模拟退火详解

写在前面

你有没有写不出正解而被迫写rand的经历？

你有没有提交rand数十次却每次都在二三十分上下浮动？

你有没有在提交一份rand代码之前净身更衣焚香？

如果是，那么，这款模拟退火正适合你！

本店开业活动为期一周，模拟退火买一送一，货到付款，逾期不候！（（（

模拟退火

模拟退火究竟是什么呢？

退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却。目的是降低硬度，改善切削加工性；降低残余应力，稳定尺寸，减少变形与裂纹倾向；细化晶粒，调整组织，消除组织缺陷。准确的说，退火是一种对材料的热处理工艺，包括金属材料、非金属材料。而且新材料的退火目的也与传统金属退火存在异同。——百度百科

模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。——百度百科

然而，你只需要知道——模拟退火是一种基于随机化的玄学算法！

现在，调用起你的物理知识，想一想，如果给你一块烧红的铁块，你怎样让它的内能降到尽量低？

你是不是会用水一下子将它浇灭？

然而这样做是错误的，正确的做法是，将它慢慢冷却，使得粒子渐趋有序，最终使得内能最小。

模拟退火就是借鉴了这样的原理。我们来打一个形象的比喻：

假设有一只兔子喝醉了酒，想要从山脚跳到山顶，首先，它肯定是这样跳的：

但是，由于它喝醉了，它可能跳到山的另一端：

由此，它会一直在山的左右两侧来回跳，但是可能永远也不会到达山顶。

但是，由于兔子在逐渐清醒，方向感也越来越强，所以兔子将会越来越逼近山顶。

在实际问题中，我们也不能够找到精确解，但是可以无限逼近解。因为精度的限制，可以通过题目。

那这就是模拟退火了么？

很遗憾，并不是。这种算法叫做爬山算法。当我们把爬山算法应用到下图的山脉中去：

兔子可能会在其中一个山峰兜兜转转，而不会到达最高的山峰。

因此，爬山算法适用于单峰函数。

那么，如果问题已经不是单峰函数了，我们怎样求解呢？

我们可以发现，爬山算法的劣势在于一旦进入一座山，就会一直在当前山的山峰横跳，而不会去考虑别的山峰。

那么，我们如何做到在逼近一个山峰的同时仍然去考虑别的山峰呢？

当我们对金属进行退火时，温度降低，内能减少，粒子不断稳定，粒子的移动也越来越不随机，但是仍有可能跳跃到远处。

这时，就得到了我们的模拟退火（SA）！

如果新状态的解更优则修改答案，否则以一定概率接受新状态。——OI Wiki

我们设初始温度为 $T$ ，降温系数为 $d$，我们每次让 $T=d\cdot T$，模拟降温的过程，直到温度非常接近 $0$ 为止，其中 $T$ 应该较大，$d$ 应该接近 $1$ 但小于 $1$。

那么我们怎样得到最优解呢？

我们设当前状态为 $x$，我们要通过一定的变动来得到新状态 $x^{\prime }$，这个变动值设为 $\Delta x$，则 $x^{\prime }=x+\Delta x$。按照退火原理，$\Delta x$ 应当在一个与 $T$ 成正比的范围随机选取。

我们现在考虑是否让 $x$ 变为 $x^{\prime }$。我们设当前状态和新状态的能量差为 $\Delta E$。显然，若新状态优于当前状态，即 $\Delta E<0$，肯定要接受新状态，那要是新状态劣于当前状态呢？

为了不被限制于局部最优解（即当前的山峰），我们应当按照一定的概率选择这个新状态。这个概率是 $e^{\frac{-\Delta E}{T}}$。

以上的过程，叫做Metropolis接受准则，原本是物理学中的概念。1982年，Kirkpatrick等人意识到固体退火过程与组合优化问题之间存在的类似性，应该把Metropolis准则引入到优化过程中来，从而诞生了模拟退火算法。

准确的来说，接受新状态的概率为

$$P(\Delta E)=\{\begin{array}{ll}1& \Delta E\le 0\\ e^{\frac{-\Delta E}{T}}& \Delta E>0\end{array}$$

物理的东西我也不会证明，当做结论记住就行……

这里附一张图：

可见随着温度逐渐降低，最优解越来越稳定。

模拟退火的关键问题，也是最麻烦的问题就是调参，初始温度，降温系数以及精度要求，麻烦得很。

至于怎么调吗，毕竟是玄学算法，我现在也只能凭直觉玄学调参……

例题

洛谷 P1337 [JSOI2004]平衡点 / 吊打XXX

这算是模拟退火的经典例题了，尽管它的正解是计算几何……

首先，这题需要用到一点物理知识：在平衡状态时，物体的总能量最小，即要选取一个点，使得 $\sum _{i=1}^n{d}_i{w}_i$ 最小，其中 $d_i$ 为选取的点到 $i$ 点的距离。

那么我们就可以模拟退火了！不断更新 $x$ 和 $y$ 使得能量最小，这题还算挺简单的，随便调调参就过了。

代码

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<ctime>
using namespace std;
struct node
{
	double x;
	double y;
	double w;
}
a[1010];
int n;
double ansx,ansy;
const double Down=0.998;
double energy(double x,double y)//算能量
{
	double res=0;
	for(int i=1;i<=n;i++)
	{
		double dx=x-a[i].x;
		double dy=y-a[i].y;
		res+=sqrt(dx*dx+dy*dy)*a[i].w;
	}
	return res;
}
void SA()
{
	double T=3000;//初温
	while(T>1e-15)//直到温度接近于0
	{
		double tx=ansx+(rand()*2-RAND_MAX)*T;//RAND_MAX是自带常数，为rand的最大值
		double ty=ansy+(rand()*2-RAND_MAX)*T;//rand()*2-RAND_MAX，随机产生-RAND_MAX到RAND_MAX的整数
		double d=energy(tx,ty)-energy(ansx,ansy);
		if(d<0)
		{
			ansx=tx;
			ansy=ty;
		}
		else if(exp(-d/T)*RAND_MAX>rand())//exp(x)为e的x次方，rand()/RAND_MAX随机产生(0,1]的小数，乘过去避免精度损失
		{
			ansx=tx;
			ansy=ty;
		}
		T*=Down;//降温
	}
}
int main()
{
	srand(time(0));
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		scanf("%lf%lf%lf",&a[i].x,&a[i].y,&a[i].w);
	for(int i=1;i<=n;i++)
	{
		ansx+=a[i].x;
		ansy+=a[i].y;
	}
	ansx/=(double)n;
	ansy/=(double)n;
	SA();
	printf("%.3lf %.3lf",ansx,ansy);
	return 0;
}