模拟退火

U 盘坏了导致我以前的东西都丢掉了……/wq/ll，所以这大概是我新 U 盘里第一个 .md 文件（？）

模拟退火（simulate anneal，简称 SA），说白了就是一种玄学的求某个函数 (f) 最大值/最小值的算法，其中 (f) 可以是一元函数、二元函数……等等，它的定义域也可以是 (mathbb{R})，或 (mathbb{Z})，甚至可以是某个长度为 (n) 的数组。它的原理，来自于物理学给高温物体降温使其内能达到最低的原理：

一个处于很高温度的物体，现在要给它降温，使物体内能降到最低。我们常规的思维是，越快越好，让它的温度迅速地降低。

然而，实际上，过快地降温使得物体来不及有序地收缩，难以形成结晶。而结晶态，才是物体真正内能降到最低的形态。

正确的做法，是徐徐降温，也就是退火，才能使得物体的每一个粒子都有足够的时间找到自己的最佳位置并紧密有序地排列。开始温度高的时候，粒子活跃地运动并逐渐找到一个合适的状态。在这过程中温度也会越降越低，温度低下来了，那么粒子也渐渐稳定下来，相较于以前不那么活跃了。这时候就可以慢慢形成最终稳定的结晶态了。

以上文字摘自 FlashHu 大佬的博客，不要问我我不懂物理

说人话 用 OI 语言来描述，以求函数最小值为例，为了模拟物理上的降温过程，我们要先给函数一个初始位置 (x_0)，以及一个初始温度 (T_0)（一般取 (2000sim 3000)），同时设定一个温度变化量 (Delta T)（一般取 (0.995sim 0.9999)）表示下一步温度会变为原来的 (Delta T) 倍，并在退火的过程中不断记录最优解 (x_{min})，以及当前位置 (x)（初始 (x=x_0)）。每一步对当前位置进行扰动，即随机一个新的值 (x') 并求出 (f(x')) 的值， 其中 (x') 随机的范围与 (T) 成正比。（如果 (f) 的定义域是一个数组也可以通过 swap 数组中某两个元素实现扰动）如果 (f(x')) 的值 (<) (f(x))，那么我们肯定会贪心地令 (x=x')，而对于 (f(x')ge f(x)) 的情况，如果我们不以一定概率接受它，那么我们就难以跳出局部最小值，搜寻更小的 (f(x)) 了，因此我们考虑记 (d=f(x')-f(x_{min}))，经过物理学家精密的计算我们接受这个解的概率应为 (e^{d/T})，也就是说我们随机一个 rand()，如果随机出的这个值小于 (exp(dfrac{d}{T}) imes ext{RAND_MAX})，那我们就接受 (x')，否则我们不改变 (x) 的值，如此进行下去直到 (T) 趋近于 (0)（可以设一个末温 (T_{min})，一般取 (10^{-14}sim 10^{-7})）即可那么最终得到的 (x_{min}) 就有一定概率是函数的最值点。

正确性证明？问物理学家，not me，thx

但 zszz 模拟退火仅仅只是一个玄学骗分算法（虽然它有一定概率能够取得不错的成绩），如果我们力求极致，想要通过模拟退火在一道题中拿到满分的话就要在调参上花费不少心思，这里给了一些模拟退火调参的技巧：

1. 相较于其他参数而言，(Delta T) 的设置至关重要，因为 (Delta T) 设小了就会导致降温不彻底，无法获得最小值，会 WA，设大了就会导致 TLE。对于 (T_0=3000,T_{min}=10^{-10},Delta T=0.999)，单次模拟退火的操作次数大约在 (2.5 imes 10^4)，由此可以估算模拟退火的效率以及 (Delta T) 的取值（同时也说明模拟退火只能适用于数据范围比较小的题，单次计算 (f) 的操作次数应顶多在 (10^3) 数量级）

2. 有时候一次模拟退火无法求得最小值，此时应进行多次模拟退火，最好卡个 clock()，并且最好从第二次模拟退火开始，将每次模拟退火的初始值设为上一次的最优解。

3. 调 (Delta T) 也可以用观察法，一边退火一边输出当前的温度等信息，大致感受一下解的降低速率，一般来说，如果解的降低速率比较均匀，就能较好地找到最优解，不均匀的话，就调整参数，将 (Delta T) 调大一点，速率就能减慢一点。反之同理。

4. 如果函数一直陷入局部最优解，那么可以尝试将 (d) 改为 (f(x')-f(x))，即当前解的大小减去上一次解的大小，这样如果 (f(x')<f(x_{min})) 就正常更新，否则将接受最优解概率中的 (d) 改为 (f(x')-f(x))，就可以比较有效地避免这种情况了。

5. 如果每次 WA 的点都一样，那么有可能是你模拟退火之外的其他地方写锅了，也有可能是你种子没有随机（即没有写 srand(time(0))（不过这好像跟模拟退火调参没啥关系啊/cy/cy）。

例题：

1. P1337 [JSOI2004]平衡点 / 吊打XXX

模板题，没啥好说的，正解是计算几何，但是用模拟退火实现很方便。

贴个代码：

const int MAXN=1000;
int n,x[MAXN+5],y[MAXN+5],w[MAXN+5];
double calc(double sx,double sy){
	double res=0;
	for(int i=1;i<=n;i++)
		res+=sqrt((sx-x[i])*(sx-x[i])+(sy-y[i])*(sy-y[i]))*w[i];
	return res;
}
double ansx=0,ansy=0,res;
void simulate_anneal(){
	double curx=ansx,cury=ansy,res=calc(curx,cury);
	for(double T=3000;T>1e-14;T*=0.993){
		double new_x=curx+(rand()*2-RAND_MAX)*T;
		double new_y=cury+(rand()*2-RAND_MAX)*T;
		double new_res=calc(new_x,new_y),dif=new_res-res;
		if(dif<0){
			curx=new_x;cury=new_y;
			ansx=new_x;ansy=new_y;
			res=new_res;
		} else if(rand()<RAND_MAX*exp(-dif/T)) curx=new_x,cury=new_y;
	}
}
int main(){
	scanf("%d",&n);double sx=0,sy=0;srand(time(0));
	for(int i=1;i<=n;i++) scanf("%d%d%d",&x[i],&y[i],&w[i]),sx+=x[i],sy+=y[i];
	ansx=0;ansy=0;res=calc(ansx,ansy);
	simulate_anneal();simulate_anneal();simulate_anneal();
	simulate_anneal();simulate_anneal();simulate_anneal();
	simulate_anneal();simulate_anneal();simulate_anneal();
	printf("%.3lf %.3lf
",ansx,ansy);
	return 0;
}

2. P3878 [TJOI2010]分金币

模拟退火的题感觉都没啥好说的，就是调参、调参、再调参。

这题就是上文中说的定义域是一个数组的题，每次扰动时交换两个值并计算出新的解即可实现模拟退火。

3. P5544 [JSOI2016]炸弹攻击1

还是没啥好说的……注意卡个时，还有一个注意点就是由于对于大部分 (f(x,y)) 都是 (0)，也就是说模拟退火一开始扰动幅度很大，随机到的 (f(x,y)) 都是 (0)，接受它们的概率为 (1)，就会导致到后面无法搜寻到最优解，因此这里有一个 trick，就是每次扰动从极值点开始扰动，然后一开始极值点横纵坐标设为 (m) 个敌人坐标的横纵坐标的平均值，这样就没啥问题了。

4. P3936 Coloring

代码五分钟，调参两小时，心态崩溃。。。交了整整两页半。。。

还是按照模拟退火的套路，每次随机交换两个位置上的值，按照之前的公式计算交换的概率，然后进行多次模拟退火 (+) 卡 clock()，这题 ( ext{TL}=5 ext{s}) 够你卡的了，只不过这题有不少注意点：

• 每次交换两个值对答案的贡献是可以 (mathcal O(1)) 计算的，duck 不 bee (mathcal O(nm)) 重新扫一遍矩阵，这样可以提高模拟退火的效率，增加模拟退火次数。
• 此题单纯地交换两个值比较容易陷入局部最优解，因此可以考虑按照调参套路中的第 (4) 条将 (d) 改为新随机出来的函数值减去原来的函数值计算概率

5. P2538 [SCOI2008]城堡

还是模板题，只不过之前我一直 ( ext{WA}) (52)，一开始以为是调参的问题，结果后来不论怎么调总是 (52)，因此查了下其他地方的代码，发现是最后读入的 (x) 没有加 (1)（读入下标从 (0) 开始）。。。。。。改了改就阿掉了。。。。。。

所以现在知道为啥我今天突然多了 (100) 多发提交了吧（