前言:

模拟退火网上已经有很多大佬讲解的很清楚也很详细了,所以本篇暂时只是我个人的记录(所以也写上了”随笔”两字),本篇末尾我会放大佬的链接,有需要的小伙伴还是看大佬的吧🤣

当然很多题目上,退火终究不是正解,而且也有点玄学,还是慎用为好。

【转载说明】本文优先发布于我的个人博客www.226yzy.com ,转载请注明出处并注明作者:星空下的YZY。

本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0许可协议。

算法简介

该部分内容参考了百度百科的描述

模拟退火是一种随机化算法【所以些许玄学,毕竟有用到随机数🤣】。模拟退火算法赋予了搜索过程一种时变且最终趋于零的概率突跳性,从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

该部分内容参考了OI-wiki的描述

接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S,否则以概率exp(-ΔT/T)接受S′作为新的当前解S。(百度百科)

先用一句话概括:如果新状态的解更优则修改答案,否则以一定概率接受新状态。

我们定义当前温度为 $T$,新状态与已知状态(由已知状态通过随机的方式得到)之间的能量(值)差为 $\Delta E$($\Delta E\geqslant 0$),则发生状态转移(修改最优解)的概率为

注意:我们有时为了使得到的解更有质量,会在模拟退火结束后,以当前温度在得到的解附近多次随机状态,尝试得到更优的解(其过程与模拟退火相似)。【确实,这样被玄学坑的概率也许可以减低😂】

伪代码

以《算法竞赛入门到进阶》/罗勇军、郭卫斌著 中的描述为参考

1
2
3
4
5
6
7
8
9
10
11
12
13
14
eps=1e-8;	//终止温度,接近于0,用于控制精度
T=100; //初始温度,应为较高的温度
delta=0.98 //降温系数,控制火的快慢,应小于1
g(x); //状态x时的评价函数,自定义
now,next; //当前状态和新状态
while(T>eps){ //如果温度未降到eps
g(now);g(next); //计算当前状态和新状态时的函数值
dE=g(now)-g(next); 当前状态和新状态的函数值的差
if(dE>=0) //新状态更优,接受新状态
now=next;
else if(exp(dE/T)>rand()) //如果状态更差,在一定概率上接受它【至于这个exp(dE/T),看上面吧】
now=next;
T*=delta; //降温,模拟退火过程
}

另外是百度百科上的伪代码

1
2
3
4
5
6
7
8
9
s:=s0;e:=E(s)//设定目前状态为s0,其能量E(s0)
k:=0//评估次数k
while k<kmax and e>emax//若还有时间(评估次数k还不到kmax)且结果还不够好(能量e不够低)则:
sn:=neighbour(s)//随机选取一临近状态sn
en:=Esn)//sn的能量为E(sn)
if random()<P(e,en,temp(k/kmax)) then//决定是否移至临近状态sn
s:=sn; e:=en//移至临近状态sn
k:=k+1//评估完成,次数k加一
returns//回转状态s

例题

模拟退火的典型应用有函数最值问题、TSP旅行商问题、最小圆覆盖、最小球覆盖等

Problem - 2899 (hdu.edu.cn)

以《算法竞赛入门到进阶》/罗勇军、郭卫斌著 中的代码为参考

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include<bits/stdc++.h>
using namespace std;
const double eps = 1e-8; //终止温度
double y;
double func(double x) { //计算函数值
return 6*pow(x,7.0)+8*pow(x,6.0)+7*pow(x,3.0)+5*pow(x,2.0)-y*x;
}
double solve() {
double T = 100; //初始温度
double delta = 0.98; //降温系数
double x = 50.0; //x的初始值
double now = func(x); //计算初始函数值
double ans = now; //返回值
while(T > eps) { //eps是终止温度
int f[2] = {1, -1};
double newx = x + f[rand()%2]*T; //按概率改变x,随T的降温而减少
if(newx >= 0 && newx <= 100) {
double next = func(newx);
ans = min(ans, next);
if(now - next > eps) {x = newx; now = next;} //更新x
}
T *= delta;
}
return ans;
}
int main(){
int cas;
scanf("%d", &cas);
while(cas--){
scanf("%lf", &y);
printf("%.4f\n", solve());
}
return 0;
}

上面的一些参数你根据题目适当修改也是可以的

当然、这道题我也尝试过可以用二分或三分做🤣

二分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
double y;
int f[2]={1,-1};
double ff(double x) {
return 6*pow(x,7.0)+8*pow(x,6.0)+7*pow(x,3.0)+5*pow(x,2.0)-y*x;
}
double ffc(double x) {
return 42*pow(x,6.0)+48*pow(x,5.0)+21*pow(x,2.0)+10*x-y;
}
double solve(){
double l=0.0,r=100.0;
while(r-l>eps){
double mid=(l+r)/2.0;
if(ffc(mid)<eps)l=mid;
else r=mid;
}
return ff(l);
}
int main(){
int t;
scanf("%d",&t);
while(t--){
scanf("%lf",&y);
printf("%.4f\n", solve());
}
return 0;
}

三分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
double y;
int f[2]={1,-1};
double ff(double x) {
return 6*pow(x,7.0)+8*pow(x,6.0)+7*pow(x,3.0)+5*pow(x,2.0)-y*x;
}
double solve(){
double T=50;
double delta=0.98;
double x=50;
double now=ff(x);
double ans=now;
while(T>eps){
double newx=x+f[rand()%2]*T;
if(newx>=0&&newx<=100){
double next=ff(newx);
ans=min(ans,next);
if(now-next>eps){
x=newx;
now=next;
}
}
T*=delta;
}
return ans;
}
int main(){
int t;
scanf("%d",&t);
while(t--){
scanf("%lf",&y);
printf("%.4f\n", solve());
}
return 0;
}

Problem - 3007 (hdu.edu.cn)

(最小覆盖圆)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
double ansx,ansy,ansr;
struct Point{
double x,y;
}p[5010];
Point maxp;
int n,f[2] = {1, -1};
double maxR(double x,double y){
double r=0;
for(int i=0;i<n;i++){
r=max(pow(x-p[i].x,2)+pow(y-p[i].y,2),r);
}
return sqrt(r);
}
void solve(){
double T=100;
double delta=0.98;
ansr=maxR(ansx,ansy);
while(T>eps){
double newx=ansx+f[rand()%2]*T/10;
double newy=ansy+f[rand()%2]*T/10;
double newr=maxR(newx,newy);
if(ansr>newr){
ansr=newr;
ansx=newx;
ansy=newy;
}
T*=delta;
}
}
int main(){
while(~scanf("%d",&n),n){
ansx=0;ansy=0;
for(int i=0;i<n;i++){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
solve();
printf("%.2lf %.2lf %.2lf\n",ansx,ansy,ansr);
}
}

这道题数据量不大,我用模拟退火是156MS

当然,正解肯定是几何法啦,46MS

所以,模拟退火还是要谨慎使用,参数要自己控制好😂

几何法

以《算法竞赛入门到进阶》/罗勇军、郭卫斌著 中的代码为参考

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const int maxn = 505;
int sgn(double x) {
if (fabs(x) < eps)return 0;
else return x < 0 ? -1 : 1;
}
struct Point {
double x, y;
};
double Distance(Point A, Point B) { return hypot(A.x - B.x, A.y - B.y);}
//求三角形abc的外接圆的圆心
Point circle_center(const Point a,const Point b,const Point c){
Point center;
double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2;
double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2;
double d = a1 * b2 - a2 * b1;
center.x = a.x + (c1 * b2 - c2 * b1) / d;
center.y = a.y + (a1 * c2 - a2 * c1) / d;
return center;
}
//求最小覆盖圆,返回圆心c,半径r
void min_cover_circle(Point *p,int n,Point &c,double &r){
random_shuffle(p, p + n); //随机函数,打乱所有点
c = p[0];r = 0; //从第一个点p0开始,圆心为p0,半径为0
for (int i = 1; i < n;i++){ //扩展所有点
if(sgn(Distance(p[i],c) - r) > 0){ //点pi在圆外
c = p[i];r = 0; //重新设置圆心为pi,半径为0
for (int j = 0; j < i;j ++){ //重新检测前面所有点
if(sgn(Distance(p[j],c) - r) > 0){ //两点定圆
c.x = (p[i].x + p[j].x) / 2;
c.y = (p[i].y + p[j].y) / 2;
r = Distance(p[j], c);
for (int k = 0; k < j;k ++){
if(sgn(Distance(p[k],c) - r) > 0){ //两点定圆不行就是三点定圆
c = circle_center(p[i], p[j], p[k]);
r = Distance(p[i], c);
}
}
}
}
}
}
}
int main(){
int n;
Point p[maxn];
Point c;double r;
while(~scanf("%d", & n) && n){
for (int i = 0; i < n;i ++) scanf("%lf %lf",&p[i].x,&p[i].y);
min_cover_circle(p, n, c, r);
printf("%.2f %.2f %.2f\n", c.x, c.y, r);
}
return 0;
}

random_shuffle(p, p + n);这个不加这个题也能过,加这个是为了随机打乱,使得从概率角度将该程序的复杂度接近O(n)【代码中有3层循环,如果不随机(也就是万一出现比较恶心的样例)是O(n³)】

其他题目

后续再补上代码

P1337JSOI2004]平衡点 / 吊打XXX - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

1379 — Run Away (poj.org)

2069 — Super Star (poj.org)

待续。。。

参考文献及推荐

【朝夕的ACM笔记】随机化-爬山与模拟退火(完全通俗解释) - 知乎 (zhihu.com)

模拟退火算法_百度百科 (baidu.com)

模拟退火 - OI Wiki (oi-wiki.org)

最后

暂时就写到这了,部分内容还不完善

欢迎访问我的小破站https://www.226yzy.com/ 或者GitHub版镜像 https://226yzy.github.io/ 或Gitee版镜像https://yzy226.gitee.io/

我的Github:226YZY (星空下的YZY) (github.com)

【转载说明】本文优先发布于我的个人博客www.226yzy.com ,转载请注明出处并注明作者:星空下的YZY。

本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0许可协议。