本篇博文部分搬运于M_sea，旨在参考并给出清晰系统的讲解。原文参见luogu日报。
简介
模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。
模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。
——百度百科
简单说，模拟退火是一种随机化算法。当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，我们常使用模拟退火求解。
它与爬山算法最大的不同是，在寻找到一个局部最优解时，赋予了它一个跳出去的概率，也就有更大的机会能找到全局最优解。
原理
模拟退火的原理也和金属退火的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。
——百度百科
要将模拟退火，首先要知道金属退火。
简单来说，就是将金属加热到一定温度，保持足够时间，然后以适宜速度冷却。
那么对应到OI上，就是每次随机出一个新解，如果这个解更优，则接受它，否则以一个与温度和与最优解的差相关的概率接受它。
设这个新的解与最优解的差为 \Delta E，温度为 T， k为一个随机数，那么这个概率为 e^{\frac{\Delta E}{kT}}
过程
降温
模拟退火时有三个参数，分别是初始温度 T_0、降温系数 \Delta、终止温度T_k。
其中， T_0是一个比较大的数，\Delta是一个略小于 1的正数， T_k是一个略大于 0的正数。
我们先让温度 T=T_0，然后每次降温时 T=T⋅Δ，直到 T≤T_k为止。
大致过程如下：
可以看出，随着温度的降低，解逐渐稳定下来，并逐渐集中在最优解附近。
其它
程序开始时，我们要先srand(一个常数)。这个常数可以决定分数。你可以使用233333，2147483647，甚至某个八位质数。
一遍SA往往无法跑出最优解，所以可以多跑几遍。
可以用一个全局变量记录所有跑过的SA的最优解，每次从那个最优解开始继续SA，可以减小误差。
时间复杂度
时间复杂度O(\text{玄学})。
一般降温系数 \Delta与 1的差减少一个数量级, 耗时大约多 10倍； T_0和 T_k变化一个数量级, 耗时不会变化很大。
实战
这里以洛谷1337 [JSOI2004]平衡点 / 吊打XXX为例，讲解SA的实际应用。
题目要使整个系统的能量最小。那么我们只要用SA跑出这个最小值即可。
#include <bits/stdc++.h>
#define re register
using namespace std;
inline int read() { //读入优化
int X=0,w=1; char c=getchar();
while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
while (c>='0'&&c<='9') X=(X<<3)+(X<<1)+c-'0',c=getchar();
return X*w;
}
struct node { int x,y,w; };
node a[1010];
int n,sx,sy;
double ansx,ansy; //全局最优解的坐标
double ans=1e18,t; //全局最优解、温度
const double delta=0.993; //降温系数
inline double calc_energy(double x,double y) { //计算整个系统的能量
double rt=0;
for (re int i=1;i<=n;i++) {
double deltax=x-a[i].x,deltay=y-a[i].y;
rt+=sqrt(deltax*deltax+deltay*deltay)*a[i].w;
}
return rt;
}
inline void simulate_anneal() { //SA主过程
double x=ansx,y=ansy;
t=2000; //初始温度
while (t>1e-14) {
double X=x+((rand()<<1)-RAND_MAX)*t;
double Y=y+((rand()<<1)-RAND_MAX)*t; //得出一个新的坐标
double now=calc_energy(X,Y);
double Delta=now-ans;
if (Delta<0) { //接受
x=X,y=Y;
ansx=x,ansy=y,ans=now;
}
else if (exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y; //以一个概率接受
t*=delta;
}
}
inline void Solve() { //多跑几遍SA，减小误差
ansx=(double)sx/n,ansy=(double)sy/n; //从平均值开始更容易接近最优解
simulate_anneal();
simulate_anneal();
simulate_anneal();
}
int main() {
srand(19260817);//玄学srand
n=read();
for (re int i=1;i<=n;i++) {
a[i].x=read(),a[i].y=read(),a[i].w=read();
sx+=a[i].x,sy+=a[i].y;
}
Solve();
printf("%.3f %.3f\n",ansx,ansy);
return 0;
}
如何调参
不同的 \Delta、 T_0， T_k，甚至 srand()和SA的次数都会影响到答案。
我们探讨一下SA的玄学调参。
Q：答案不是最优的怎么办？
A：有以下几种方法：调大 \Delta、调整 T_0和 T_k，以及多跑几遍SA。
当您的程序跑的比较快时，可以选择多跑几遍SA，或者调大 \Delta，从而增大得到最优解的概率。
T_0和 T_k的值需要根据值域对应调整，从而保证每次随机出的变化幅度不会跑太远。
Q：还是跑不出最优解怎么办？
A：那可能是您太非了。尝试更换随机数种子，或者 srand(rand())，总之，总有可能跑出正解。
Q：我是非酋，交了两页也没有用模拟退火AC，怎么办？
A：您还是选择打正解吧。
如何生成新解
- 坐标系内：随机生成一个点，或者生成一个向量。
- 序列问题： random_shuffle或者随机交换两个数。
- 网格问题：可以看做二维序列，每次交换两个格子即可。
0 条评论