浅谈玄学算法——模拟退火
初级篇
本篇讲解SA的基本概念。
如果您已经了解SA的基本概念，您可以跳过这一段。
简介
模拟退火算法（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}}$
Metropolis准则(1953)—以概率接受新状态
若在温度T，当前状态i → 新状态j
若Ej<Ei，则接受 j 为当前状态；
否则，若概率 p=exp[-(Ej-Ei)/KT] 大于[0,1)区间的随机数，则仍接受状态 j 为当前状态；若不成立，则保留状态 i 为当前状态。
p=$e^{frac{-(Ej-Ei)}{kT}}$ ：在高温下，可接受与当前状态能量差较大的新状态；在低温下，只接受与当前状态能量差较小的新状态。
玻尔兹曼常量系热力学的一个基本常量，记为“k”，数值为：k=1.380649 × 10-23 J/K，玻尔兹曼常量可以推导得到：理想气体常数R等于玻尔兹曼常数乘以阿伏伽德罗常数（即R=k·NA）。
2018年11月16日，国际计量大会通过决议，1开尔文定义为“对应玻尔兹曼常数为1.380649×10-23 J/K的热力学温度” [1] 。新的定义于2019年5月20日起正式生效 [1] 。
过程
降温
模拟退火时有三个参数，分别是初始温度 $T_0$ 、降温系数 $Delta$ 、终止温度 $T_k$ 。
其中， $T_0$ 是一个比较大的数， $Delta$ 是一个略小于 $1$ 的正数， $T_k$ 是一个略大于 $0$ 的正数。
我们先让温度 $T=T_0$ ，然后每次降温时 $T=Tcdot Delta$ ，直到 $Tleq T_k$ 为止。
大致过程如下：
可以看出，随着温度的降低，解逐渐稳定下来，并逐渐集中在最优解附近。
其它
程序开始时，我们要先srand(一个常数)。这个常数可以决定分数。你可以使用233333，2147483647，甚至某个八位质数。
一遍SA往往无法跑出最优解，所以可以多跑几遍。
可以用一个全局变量记录所有跑过的SA的最优解，每次从那个最优解开始继续SA，可以减小误差。
时间复杂度
时间复杂度 $O( ext{玄学})$ 。
一般降温系数 $Delta$ 与 $1$ 的差减少一个数量级, 耗时大约多 $10$ 倍； $T_0$ 和 $T_k$ 变化一个数量级, 耗时不会变化很大。
中级篇
本篇讲解SA的实际应用。
如何调参
这是我P1337的提交记录。可以看到，不同的 $Delta$ 、 $T_0$ ， $T_k$ ，甚至 $srand()$ 和SA的次数都会影响到答案。
我们探讨一下SA的玄学调参。
Q：答案不是最优的怎么办？
A：有以下几种方法：调大 $Delta$ 、调大 $T_0$ 、调小 $T_k$ ，以及多跑几遍SA。
当您的程序跑的比较快时，可以选择多跑几遍SA，或者调大 $Delta$ ，从而增大得到最优解的概率。
调大 $T_0$ 和调小 $T_k$ 也可以，而且时间并不会增大太多。
Q：还是跑不出最优解怎么办？
A：那可能是您太非了。 $quad$ 尝试更换随机数种子，或者 $srand(rand())$ ，总之，总有可能跑出正解。
Q：我是非酋，交了两页也没有用模拟退火AC，怎么办？
A：您还是选择打正解吧。
如何生成新解
- 坐标系内：随机生成一个点，或者生成一个向量。
- 序列问题： $random\_shuffle$ 或者随机交换两个数。
- 网格问题：可以看做二维序列，每次交换两个格子即可。
例题
这里以洛谷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(1******7); srand(rand()); srand(rand()); //玄学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
",ansx,ansy);
return 0;
}这份代码跑了 $476 ms$ ，可以通过，而且时间比较充裕。
高级篇
本篇讲解SA的其它毒瘤应用，以及如何结合其它算法。
分块模拟退火
这是管理员@ComeIntoPower 在第一次审稿时要我加上去的内容。
我们发现，有时候模拟退火是不适用的，比如这种情况：
此时函数的峰特别多，所以我们要用分块模拟退火的做法。
大致算法是：将其分为几块，然后对每块跑一遍SA，最后再合并答案。
然而好像没有要用这个的题。不过万一以后考了呢【雾】
upd：这里块的数量不是 $sqrt{n}$ ，而是一个比较小的数。要根据不同的题用不同的大小。
模拟退火套其它算法
其实只是在SA的板子内，用其它算法获得答案而已。
例题
考虑用SA跑 $sum_{i=1}^n(x_i-overline x)$ 的最小值。
我们可以随机一个序列，然后用DP跑出这个序列的最小。
其它的就是板子。
事实证明我还是太非，交了一页多才过，而且还暴力枚举了一下qwq
#include <bits/stdc++.h>
#define re register
#define sqr(x) ((x)*(x))
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;
}
const double delta=0.993;
int n,m;
int a[30],s[30];
double f[30][30];
inline double DP() { //DP跑出解
memset(f,127,sizeof(f));
f[0][0]=0;
for (re int i=1;i<=n;i++) s[i]=s[i-1]+a[i];
for (re int i=1;i<=n;i++)
for (re int j=1;j<=i;j++)
for (re int k=0;k<i;k++)
f[i][j]=min(f[i][j],f[k][j-1]+sqr(s[i]-s[k]-(double)s[n]/m));
return f[n][m];
}
double ans;
inline void SA() { //模拟退火
double nowans=ans,T=2000;
while (T>1e-6) {
int x=0,y=0;
while (x==y) x=rand()%n+1,y=rand()%n+1;
swap(a[x],a[y]);
double now=DP();
double Delta=now-ans;
if (Delta<0) ans=nowans=now;
else if (exp(-Delta/T)*RAND_MAX>rand()) nowans=now;
else swap(a[x],a[y]);
T*=delta;
}
for (re int i=1;i<=n;i++)
for (re int j=1;j<=n;j++) {
swap(a[i],a[j]);
ans=min(ans,DP());
swap(a[i],a[j]);
}
}
inline void Solve() {
ans=1e10;
while ((double)clock()/CLOCKS_PER_SEC<0.8) SA(); //卡时
}
int main() {
srand(1******7); srand(rand()); srand(rand());
n=read(),m=read();
for (re int i=1;i<=n;i++) a[i]=read();
Solve();
if (ans/m<1e-6) printf("0.00
");
else printf("%.2f
",sqrt(ans/m));
return 0;
}在不会TLE的情况下尽量多地跑SA
我们知道，有一个 $clock()$ 函数，返回程序运行时间。
那么这样即可：
while ((double)clock()/CLOCKS_PER_SEC<MAX_TIME) SA();其中MAX_TIME是一个自定义的略小于1的正数，可以取0.7~0.8。
习题
自己找题去吧（逃）
反正什么求最优解的题都可以用模拟退火搞一搞