再谈模拟退火

# 引言 ~~在实际生活中，~~我们常常会遇到求函数最值的问题，那怎么办呢？我们当然可以选择爬山算法，即每次在当前最优解的附近选择一个解，如果它优于最优解，就接受它，否则不接受它，并调小选择范围，寻找下一个解。在某些情况下，它是适用的，比如下图  但这个算法的劣势非常明显——它会被局限在一个局部最优解上，无法取得全局最优解，比如下图这个函数。  这时，我们就可以使用一个玄学算法——模拟退火。 # 正文 ## 简介 >模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Čern&yacute;在1985年也独立发明此演算法。模拟退火算法是解决[TSP问题](https://baike.baidu.com/item/TSP/2905216)的有效方法之一。 > >模拟退火的出发点是基于物理中固体物质的退火过程与一般[组合优化问题](https://baike.baidu.com/item/组合优化/3314860)之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。 > >### 原理 > >模拟退火的原理也和金属[退火](https://baike.baidu.com/item/退火/1039313)的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。 ——摘自百度百科~~（当然你也可以不看）~~ 简单的说，模拟退火就是在一种一定范围内求多峰函数最值的算法。它在模拟温度降低的同时得出新解，温度越高，解的变化量越大，随着温度的逐渐降低，解的变化量也渐渐变小，并越发集中在最优解附近。最后温度达到了我们设置的最低温，对应到物理学上也就是结晶了，这时，我们可以认为当前我们取得的解就是最优解，~~当然也可能不是，~~整个算法也就终止了。 ## 过程 我们先引入几个参数：当前最优解$E_0$，新解$E$，解变动量$ΔE$（**$E$与$E_0$的差**），上一个被接受的解$E_1$，初温$T_0$，末温$T_k$，当前温度$T$，温度变动量$Δ$，再引用一张非常经典的图——  这张图非常好的展现了模拟退火的运行过程，从$T_0$开始，每次乘上$Δ$得到$T$，如果$T$小于$T_k$则终止降温。$T_0$我一般设置在$1000$~$5000$左右，$Δ$则是一个略小于1的常数，而$T_k$一般设置在1e-8到1e-15之间(或者另外一个极小的数)。 在降温的同时，我们在$E_1$（不是最优解$E_0$）的基础上扰动产生新解$E$，需要注意的是扰动大小随温度的降低而变小，因为在温度高的时候，解的变化量非常大，这时的任务是在全局范围中找到最优解的大致位置，随着温度的降低，解渐渐稳定，这时的任务是确定最优解的准确位置。 但每次得出新解以后，我们以什么原则，或者说什么概率来接受它呢？ 这时就要用到[Metropolis准则](https://baike.baidu.com/item/Metropolis接受准则/14678977?fr=aladdin)。简单说来，假设我们的目标是求最小值，如果$E$与$E_0$的差，也就是$ΔE$小于$0$，我们就接受当前解，因为它优于之前的最优解嘛。而如果$ΔE$大于$0$，也就是我们遇到了一个更劣的解，我们也要以一定的概率来接受它，因为我们要找的一个多峰函数的全局最小值，因此就不能局限于一个局部的凹函数。而这个概率是$\exp (-ΔE/T)$。 我个人对于这个概率的理解是这样的：对于$ΔE$，如果它较大，即我们遇到了一个劣得多的解，那我们接受它的概率就相对较小，因为$-ΔE$较小嘛；而如果$ΔE$较小，即我们遇到了一个较劣的解，我们接受它的概率就较大，因为$-ΔE$较大。对于$T$，随着时间的增加，$T$变得越来越小，因此我们把$-ΔE$除以$T$，这样接受的概率就随着温度的降低而越来越小，因为$-ΔE$是一个负数嘛。而对于整个式子，当$T$较大的时候，我们会接受大部分的解，当$T$较小的时候，我们就只会接受$ΔE$较小的解。~~关于Metropolis准则的具体证明，过于玄学，这里就不给出了。~~当然你也可以自己试一下。如果选择接受$E$，则把$E_1$设置为$E$，然后降温并寻找下一个解。 另外还有一点，@BFLSTiger指出了一个关于$ΔE$的问题。我翻阅了其他一些文章和题解，其实对于$ΔE$来说，$E-E_0$和$E-E_1$这两种写法都是存在的。这里我也尝试了一下： [提交记录P5544（$E-E_0$的写法）](https://www.luogu.com.cn/record/29264668) [提交记录P5544（$E-E_1$的写法）](https://www.luogu.com.cn/record/29721434) 两种写法都是可以AC此题的，所以这里对此持保留态度。但对于具体的题目，还是需要具体的分析与尝试。 这里再引用一张~~很糊的~~图：  到这里我们也就知道，模拟退火算法的速度和结果受参数（$T_0$，$T_k$，$Δ$还有随机数种子）的影响非常大，是一个玄学的算法，时间复杂度也是$O (玄学)$。 ## 实现 ### 例题 接下来我们结合一道例题来讲一讲模拟退火的c++​代码实现。[**UVA10228** A Star not a Tree?](https://www.luogu.com.cn/problemnew/show/UVA10228) （这道题其实洛谷上也有） 英文题面尽管跳过，大意是给定$n$个点，求其[费马点](https://baike.baidu.com/item/托里拆利点/22667515?fromtitle=费马点&fromid=3333221&fr=aladdin)（到这$n$个点的距离最小的点）到所有点的距离和。此题各部分的代码实现都很方便，其实就是一道模板题，代码如下： ```cpp #include<iostream> #include<cstdio> #include<stdlib.h> #include<iomanip> #include<cmath> #define R register #define rep(i,a,b) for(R int i=a;i<=b;i++) #define delta 0.996 #define maxn 50005 using namespace std; inline int read() { int x=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();} while('0'<=ch&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar(); return x*f; } struct node{ double x,y; }poi[maxn]; int T,n; double ansx,ansy,ax,ay,ans,t; void clear() { ax=0,ay=0; ans=1e8; } double calculate(double x,double y) { double res=0; rep(i,1,n) { double dx=x-poi[i].x,dy=y-poi[i].y; res+=sqrt(dx*dx+dy*dy); } return res; } void simulate_anneal() { double x=ansx,y=ansy; t=3000; while(t>1e-15) { double X=x+((rand()<<1)-RAND_MAX)*t; double Y=y+((rand()<<1)-RAND_MAX)*t; double now=calculate(X,Y); double Delta=now-ans; if(Delta<0) { ansx=X,ansy=Y; x=X,y=Y; ans=now; } else if(exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y; t*=delta; } } void work() { ansx=ax/n,ansy=ay/n; simulate_anneal(); simulate_anneal(); simulate_anneal(); simulate_anneal(); simulate_anneal(); } int main() { srand(1e9+7); T=read(); rep(i,1,T) { n=read(); clear(); rep(j,1,n) { poi[j].x=read(),poi[j].y=read(); ax+=poi[j].x,ay+=poi[j].y; } work(); cout<<round(ans)<<'\n'; if(i!=T) cout<<'\n'; } return 0; } ``` 有几个注意点：坐标位置，温度和解变动量必须开成double​，一是为了确保精度，二是为了防止爆int。~~还要注意输出换行，实在很坑。~~ 但是当你愉快的写玩此题并提交以后，可能会发现你并没有AC此题。记得之前说过的吗，我们得出不一定是最优解。这时候就涉及到一个麻烦的步骤——调参。通常有以下几种调参的方式： 1. 调大初温$T_0$。 2. 调小末温$T_k$。 3. 调大温度变动量$Δ$。 4. 选取一个新的随机数种子。 5. 多跑几遍模拟退火。 6. ~~开O2~~ 其中第一，二点对于运行时间的影响不大。而第三点则非常关键，一个微调都会使运行时间和结果发生巨大变化。第五点也是一个有用的方式，一般我们跑三到五遍模拟退火，如果时间充裕，你也可以适当多跑~~一两百~~几遍。而第四点就非常看脸了，但就我个人而言，最有用的还是这句随机数种子： ```cpp srand(time(0)); ``` ### 习题 学完一个新算法以后，当然应该练习啦。其实模拟退火的主过程基本就是模板了，唯一的麻烦点是对calculate()函数和接受概率的修改，比如下题：[[JSOI2016\]炸弹攻击1](https://www.luogu.com.cn/problemnew/show/P5544) 此题的calculate()函数倒是很简单，麻烦的是修改接受概率。 题目要求的是最大值，那么$-ΔE$就成了一个正数，怎么修改呢？其实此时我们只需把这句话： ```cpp else if(exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y; ``` 中的$>$号改为$<$号就可以了，如下： ```cpp else if(exp(-Delta/t)*RAND_MAX<rand()) x=X,y=Y; ``` 而这样就很玄学了。之前我说错了，因为$-ΔE$成了一个正数，所以$\exp (-Delta/t)$必定是大于1的，也就是没有接受劣解的概率。而此题$ΔE$波动小，搜寻范围大，所以我们这样写就可以手动避免算法陷入劣解不能自拔。~~但这样写的原因是我过了此题以后才想出来的。~~ 代码就略过了，实在很简单。 ---------------------- 模拟退火的应用不仅仅是求点坐标，还可以拿来求序列。其实过程也很简单，每次随机交换序列中的两个元素就可以了，而对于网格，看作是二维序列即可。下面有一道求序列的题目：[[SCOI2008\]城堡](https://www.luogu.com.cn/problemnew/show/P2538) 读完题以后，你可能不知道此题和序列有何关系。但我们其实可以这样考虑：把所有没有城堡的城市抽象成一个序列，而序列的前$k$个城市，就是要修建城堡的城市。 而关于calculate()函数，我们可以先用floyd​算法预处理出每个城市之间的距离，在这个函数中我们只需$n^2$扫描一次，求出所有城市中离最近城堡的距离的最大值就可以了。代码如下： ```cpp #include<iostream> #include<cstdio> #include<stdlib.h> #include<cmath> #include<ctime> #define rep(i,a,b) for(register int i=a;i<=b;i++) #define maxn 500 #define inf 0x3f3f3f3f #define delta 0.996 using namespace std; inline int read() { int f=1,x=0; char ch=getchar(); while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();} while('0'<=ch&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar(); return x*f; } void write(int x) { if(x<0) x=-x,putchar('-'); if(x>9) write(x/10); putchar(x%10+'0'); } struct edge{ int a,b,next,v; }e[maxn]; int head[maxn],cnt,n,m,k,v[maxn],cas[maxn]; int dis[maxn][maxn],p[maxn],N,X[maxn],ans=1e8; double t; int calculate(int x[]) { rep(i,1,k) cas[x[i]]=1; int res=-inf; rep(i,0,n) { int minn=inf; rep(j,0,n) if(cas[j]) minn=min(minn,dis[i][j]); res=max(res,minn); } rep(i,1,k) cas[x[i]]=0; return res; } void simulate_anneal() { int a[maxn]; rep(i,1,N) a[i]=p[i]; t=5000; while(t>1e-15) { int b[maxn]; rep(i,1,N) b[i]=a[i]; int x=rand()%N+1; int y=rand()%N+1; swap(b[x],b[y]); int now=calculate(b); double Delta=now-ans; if(Delta<0) { ans=now; rep(i,1,N) p[i]=a[i]=b[i]; } else if(exp(-Delta/t)*RAND_MAX>rand()) { rep(i,1,N) a[i]=b[i]; } t*=delta; } } void work() { simulate_anneal(); simulate_anneal(); simulate_anneal(); simulate_anneal(); simulate_anneal(); } int main() { srand(time(0)); n=read(),m=read(),k=read(); n--; rep(i,0,n) X[i]=read(); rep(i,0,n) v[i]=read(); rep(i,1,m) { int a=read(); cas[a]=1; } rep(i,0,n) if(!cas[i]) p[++N]=i; rep(i,0,n) rep(j,0,n) dis[i][j]=inf; rep(i,0,n) dis[i][X[i]]=dis[X[i]][i]=min(v[i],dis[i][X[i]]),dis[i][i]=0; rep(c,0,n) rep(i,0,n) rep(j,0,n) dis[i][j]=min(dis[i][j],dis[i][c]+dis[c][j]); work(); write(ans); return 0; } ``` # 最后 要做好调参的心理准备，我在三天之内调参交了~~七页~~。 推荐几题： 1. [[JSOI2008\]球形空间产生器](https://www.luogu.com.cn/problemnew/show/P4035) 2. [[CEOI2004\]锯木厂选址](https://www.luogu.com.cn/problemnew/show/P4360) 3. [Coloring](https://www.luogu.com.cn/problemnew/show/P3936) ~~其他的习题自己找去吧。~~ update 2020.3.3 加入了$\LaTeX$数学公式渲染，并添加了一张图。 update 2020.5.1 修锅，感谢@M_sea 纠错。 update 2020.6.21 修锅，小部分改动,致谢@Caro23333。 update 2020.7.14 我又来修锅了 如仍存笔误欢迎指正。 update 2020.7.24 修了一个~~奇怪的~~锅，致谢@coolzz27和@Social_Zhao。 update 2020.8.9 作了一些补充，致谢@BFLSTiger。