[OI]模拟退火小结

模拟退火

恭喜您打开玄学世界的大门（雾

为什么要写这篇文章

模拟退火其实都已经快被讲烂了……我是说前置知识。像是爬山啊物理学啊乱七八糟讲了一大堆，然后真正的模拟退火占了一点篇幅，之后就开始做题……这啥啊……

实在看不下去了，不如自己动手，模拟退火！（大雾

简单的介绍

虽然被讲烂了，但还是要提一提的。

1.爬山算法

要学模拟退火的人应该没有不会爬山算法的。本质就是贪心。毫无技巧的贪心。在某一范围内，只要找到一个更优的解就接受他，找不到更优解就退出。错误性是很显然的。

2.退火是啥

固体退火：

1.先将固体加热至溶化然后徐徐冷却。

2.退火要徐徐进行使得系统在每一温度下都达到平衡。

3.冷却时不能急剧减温

真正的模拟退火（Simulate Anneal）

引子：

假设你现在有一个箱子，水平放在空中，箱子的底部有一些凹凸不平的坑。现在你向其中扔一个小球（足够小），假设忽略弹力作用，那么这个小球会在着落点附近向更低处移动，最后停在一个坑中。除非你是欧皇，否则这个坑有很大概率不是最深的那个坑。但是现在，我们使劲晃一晃这个箱子，那么小球就有一定概率滚出当前的坑，进入另外一个坑。这个坑是不是最深的，我们也不得而知，但是可以知道的是，从理论上来说，只要你晃得次数够多，那这个小球就一定会进入最深的那个坑。

Metropoils算法

映射到题目中，假设我们可能求出的所有解组成了一个高低错落的函数图像，这时，如果有更优解，我们一定接受它，而如果找不到最优解了，说明你找到了一个局部最优解，但不知道是不是全局最优解。那么我们委屈求全，先以一定的概率（假设为p）接受一个劣解，看看是否在这个局部最优解之外有没有更有的解。这就是模拟退火（Simulate Anneal，简称SA）的主要思想。

Metropoils算法。名字听起来很niubi，但其实就是定义了上面的转移概率p。用数学式表达起来是这样的：

\[p=\begin{cases}1 & f(n+1)<f(n)\\exp(\frac{f(n+1)-f(n)}{T}) & f(n+1)>=f(n)\end{cases} \]

这样虽然也不能保证能得到最优解，但是我们可以将这一过程循环迭代多次，就可以有较高的概率得到最优值。

为什么它比纯随机算法有更高的出解概率

听起来还不错。那么，问题来了：这个东西既然是以一定概率接受，那么就是说这个东西是个随机化算法。这个算法跟我随机在区域内瞎选有什么不同吗？

显然有很大不同。首先，瞎选的成功概率太小。为什么呢？你随机选一个点的时候，明明附近就有更优的解，你却仍固执的认为你现在这个才是最优解，试问：傻不傻？

那么问题又来了：那我每次随机选一个点，找到他附近的最优解，然后退出，重新随机选一个点，是不是就比较优了？

这就是对数据进行多次爬山算法啊。但实际效果虽然不及模拟退火，但仍然很优。

但是遇到大数据，这种算法就也不实用了，因为每次选完之后就真的完了，不会对以前的解，以后的解有一点帮助（除了判断全局最优）。这显然是很浪费的。而且这种算法及其不稳定，准确的说，是根本没有稳定性可言。对于一个区间，所选的点显然是有相同的概率均分在这个区间内。但是解不是啊，肯定会有一部分局部最优解比较少，一部分比较多啊。你如果非要较真的话，就像下面的函数图像：

VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV

行吧。你直接放弃模拟退火算了。这题没法用SA做。

这种算法的不稳定性，就像是炼铁的时候急剧升高温度，又急剧降低温度，最后的结果只能是导致炼出的铁十分的脆。为了保证随机算法的稳定性，但同时又要保证它有概率找到最优解，一代代算法研究者经过仔细研究，终于在固体退火处得到了灵感，赋予了这个随机化算法另一个操作：降温。

注意到先前提到的Metropoils算法中p的定义式，是有一个\(\frac{1}{T}\)的，这里面的T就是当前的温度。随着时间的推移，温度渐渐降低，那么转移概率p也随之升高。由此可见，选择点的不稳定性随时间增长而变大，也就是说跳出当前局部最优解的概率越大。只要T的衰减速度设计足够巧妙，就有高概率的正确答案。

奉上wiki的动图：

玄学调参

现在问题全都集中在了求T的衰减函数上。这个东西可以写成\(T=\alpha * T(0<\alpha < 1)\)的形式。那么如何求\(\alpha\)？

看到这里，恭喜你，成功的真正进入了玄学调参的大门。

之所以这么讲，是因为这个参数根本就是没有固定值的。由于每道题目的函数图像不同，这个参数的取值也有相当的改变。

至此，这个算法有了两个不确定因素：迭代次数和衰减系数。

唯一的求解办法是：

多试几组极限数据，找一个跑的最快的值。

就没有什么技巧吗？

没有。

只有一点点。对于大部分题目，我们的T一般和数据范围有关。而降温系数一般在1附近，因为越往后降温越快，是成指数级增长的。0.99其实都一点也不小。实际效果上，1-7e-3，也就是0.993效果很好。需要注意的是，降温系数不要改的比0.98小，一般降温系数与1的差减少一个数量级, 时间消耗大约多 10 倍……至于迭代次数……不超时就越大越好。

如果程序一直得不到最优解，可以考虑适当增加T和衰减系数。精度不够可以减小最终温度。

剩下你要做的就是选好参数，然后祈祷你的人品足够好到A掉这道题……

讲了这些，模拟退火的基本思路应该比较清楚了，他比其他算法的优越性也比较清楚了。

不过还有一个最重要的问题：

代码实现

Talk is cheap, show me the code!

再找一道经典例题好了：

吊打XXX

Description

gty又虐了一场比赛，被虐的蒟蒻们决定吊打gty。gty见大势不好机智的分出了n个分身，但还是被人多势众的蒟蒻抓住了。蒟蒻们将n个gty吊在n根绳子上，每根绳子穿过天台的一个洞。这n根绳子有一个公共的绳结x。吊好gty后蒟蒻们发现由于每个gty重力不同，绳结x在移动。蒟蒻wangxz脑洞大开的决定计算出x最后停留处的坐标，由于他太弱了决定向你求助。不计摩擦，不计能量损失，由于gty足够矮所以不会掉到地上。

Input

输入第一行为一个正整数n(1<=n<=10000)，表示gty的数目。
接下来n行,每行三个整数xi，yi，wi，表示第i个gty的横坐标，纵坐标和重力。
对于20%的数据，gty排列成一条直线。
对于50%的数据，1<=n<=1000。
对于100%的数据，1<=n<=10000,-100000<=xi,yi<=100000

Output

输出1行两个浮点数（保留到小数点后3位），表示最终x的横、纵坐标。

Sample Input

3
0 0 1
0 2 1
1 1 1

Sample Output

0.577 1.000

题目要求的就是让总重力势能最小。

具体的解法去看题解。这里只是展示一下SA的写法。

暂且不管解法，模拟退火部分应该还是挺好理解的。

此处有一个小技巧：（来自@rvaluedalao的博客）

找下一个解的时候有一个提高精度的小技巧: 根据当前温度决定差值的范围. 这样在降温即将结束接近最优解的时候可以有更大的概率更精确地命中最优解.

具体做法就是使用一个产生 [0,1] 随机实数的函数，将随机区间转为[-1,1]后乘上T作为差值，（也就是生成一个[-T,T]的随机值作为差值

不过实际操作的时候我们较少直接输出最终解, 而是选择在模拟退火的过程中单独维护一个解, 只在遇到更优解的时候将其更新, 增加正确率.

#include <cstring>
#include <cstdlib>
#include <cfloat>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cctype>
#include <map>
#include <queue>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;

typedef double db;
const int MAXN = 10010;
int n, px[MAXN], py[MAXN], w[MAXN];
double minE = DBL_MAX;
db ansx, ansy;

template <typename _Tp>
inline void read(_Tp &x) {
    char ch = getchar( ); bool f = 0; x = 0;
    while (!isdigit(ch)) { if (ch == '-') f = 1; ch = getchar( ); }
    while (isdigit(ch)) { x = x * 10 + ch - '0'; ch = getchar( ); }
    x = f ? -x : x;
}

db Rnd( ) { return (db)(rand( )) / (db)(RAND_MAX); }

db Dist(db x1, db y1, db x2, db y2) {
	return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

bool Accept(db delta, db tmp) {
	return delta < 0 || Rnd( ) < exp(-delta / tmp);
}
	
db Calc(db _x, db _y) {
	db tmp = 0;
	for (int i = 1; i <= n; ++ i) 
		tmp += Dist(_x, _y, (db)px[i], (db)py[i]) * (db)w[i];
	if (tmp < minE) {
		ansx = _x;
		ansy = _y;
		minE = tmp;
	}
	return tmp;
}
	
void SA(db _x, db _y, db T, db alpha, db end) {
	db tmp = T, X = _x, Y = _y, nextX, nextY;
	db now = Calc(X, Y), nxt;
	while (tmp > end) {
		nextX = X + tmp * (Rnd( ) * 2 - 1);
		nextY = Y + tmp * (Rnd( ) * 2 - 1);
		nxt = Calc(nextX, nextY);
		if (Accept(nxt - now, tmp)) {
			X = nextX;
			Y = nextY;
			now = nxt;
		}
		tmp *= alpha;
	}
	db rndx, rndy;
	for (int i = 1; i <= 1000; ++ i) {
		rndx = (ansx + tmp * (Rnd( ) * 2 - 1));
		rndy = (ansy + tmp * (Rnd( ) * 2 - 1));
		Calc(rndx, rndy);
	}
}
	
int main( ) {
	srand(time(NULL));
	read(n);
	db _x = 0.0, _y = 0.0;
	for (int i = 1; i <= n; ++ i) {
		read(px[i]), read(py[i]), read(w[i]);
		_x += (db)px[i];
		_y += (db)py[i];
	}
	_x /= n, _y /= n;
	SA(_x, _y, 1e5, 1-7e-3, 1e-3);
	printf("%.3lf %.3lf", ansx, ansy);
	return 0;
}

觉得不保险可以把模拟退火中最后的循环变大一点，还可以进行多次SA，（SA(ansx,ansy,1e5,1-7e-3,1e-3)）前提是别T掉。

值得一提的是，关于随机种子的话，其实并不推荐使用srand(time(NULL))，更推荐的是找一个实际效果比较好的定值，srand(x)。这个x的要求是用它产生的随机数要足够均匀，模拟退火的理论正确率才会上升。如果你打了好几道随机化的题目，发现有一个数srand它后正确率很高，那么恭喜你，它就是你以后的幸运数字了。

不过最好不要用广为人知的数比如某八位质数，容易被hack。

类似的随机算法还有蚁群算法，遗传算法等，想要了解玄学的同学可以自行百度。