模拟退火学习笔记
- 本文基于多篇博客编写，仅供学习使用，将会存在多处对其他文章内容的直接复制
模拟退火是一种随机化算法。当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，我们常使用模拟退火求解。用一句话概括：如果新状态的解更优则修改答案，否则以一定概率接受新状态。
——OI WIKI
原理
模拟退火来自冶金学的专有名词退火。退火是将材料加热后再经特定速率冷却，目的是增大晶粒的体积，并且减少晶格中的缺陷。材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动。退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置。
模拟退火的原理也和金属退火的原理近似：我们将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。
可以证明，模拟退火算法所得解依概率收敛到全局最优解。
———— 维基百科
大概操作
不妨先求最小值，这与金属退火的逻辑相合。
定义当前温度为\(T\)，新状态会由已知状态随机转移得到，而新旧状态间的能量差（目标函数之差）为\(\Delta E=f(\text{now_x})-f(\text{pre_x})\)。若新状态合法且更优\((\Delta E<0)\)显然直接转移，否则，按\(P(\Delta E)=e^{\frac{-\Delta E}{T}}\)的概率转移。显然，在\(\Delta E\)一定的情况下，\(T\)越小，发生转移的概率越小。这正好符合了人们对退火的直观印象，我们希望在退火完成（即\(T\rightarrow 0\)）时不再发生状态转移，来到终点——最优解。而\(\Delta E\)越小，说明前后两者状态相差越小，这是一个较容易接受的微扰，发生转移的概率也就越大。
随着温度的降低，跳跃越来越不随机，最优解所处的区间也越来越稳定（缩至一点）。
（此处为求最大值的模型）（图来自Wiki-StimulatedAnnealing）
若求的是最大值，已知求\(f(x)\)的最大值与求\(-f(x)\)的最小值逻辑类似，则定义\(\Delta E=-f(\text{now_x})-(-f(\text{pre_x}))=f(\text{pre_x})-f(\text{now_x})\)，若\(\Delta E<0\)（新状态的目标函数值更大）则直接转移，否则仍按\(P(\Delta E)=e^{\frac{-\Delta E}{T}}\)的概率转移。
具体实现
在这里我们采用指数降温，即\(T=k·T(d\rightarrow 1^-)\)。其他降温方式留待后续补充。
设定三个参数：初始温度\(T_0\)，降温系数\(k\)，终止温度\(T_k\)。为了不要过快结束尝试，\(T_0\)应较大，\(k\rightarrow 1^-,T_k\rightarrow 0^+\)。当\(T
我们可以从定义域中任取一个点开始降温，随机向某一个点进行跳跃（具体看代码）。
例题
Strange Fuction
输入高次函数的一次项系数，求该函数的最小值。
代码
[click]
#include
#include
#include
#include
#include
const double k=0.96;
const double T_0=10000;
const double T_k=0.001;
double y;
double min(double x,double y) {return xT_k)
{
double ux=x+t*(frand()*2-1);
if (ux<0.0||ux>100.0)
continue;
double uy=calc(ux);
double dE=uy-y;
ans=min(ans, uy);
if (dE<0.0||exp(-dE/t)>frand())//考虑是否转移
x=ux,y=uy;
t*=k;//指数级降温
}
for (int i=1;i<=1000;i++)
{
double ux=x+t*(frand()*2-1);
if (ux<0.0||ux>100.0)
continue;
double uy=calc(ux);
ans=min(ans, uy);
}
//在附近尝试寻找更精确的答案
return ans;
}
int read()
{
int res=0;
char ch=getchar();
while(!isdigit(ch))
ch=getchar();
while(isdigit(ch))
res=res*10+ch-'0',ch=getchar();
return res;
}
int main()
{
srand((unsigned int)time(0));
int Q=read();
while(Q--)
{
scanf("%lf",&y);
double ans=1e60;
for (int i=1;i<=2;i++)
ans=min(ans, stimulatedAnneal());
printf("%.4lf\n", ans);
}
// system("pause");
return 0;
}
[JSOI2016 炸弹攻击]
题意简述
给出建筑物（在二维平面上视作圆）坐标和半径，以及敌人（视为质点）的坐标，可在任意坐标放置一个任意轰炸半径的炸弹，求不伤害建筑物的情况下可攻击到的最大敌人数量。
代码
[click]
#include
#include
#include
#include
#include
const int maxn=10+5;
const int maxm=1e3+10;
const double maxTime=0.85;
const double T_0=5000;
const double T_k=1e-5;
const double k=0.996;
const double eps=1e-5;
int maxp=0;
int n,m,R;
struct building
{
int x,y,r;
}b[maxn];
struct enemy
{
int x,y;
}e[maxm];
int max(int x,int y) {return x>y?x:y;}
double min(double x,double y) {return xT_k)
{
double ux=x+t*(frand()*2-1);
double uy=y+t*(frand()*2-1);
if (!check(ux, uy))
continue;
double ur=getRadius(ux, uy);
if (ur<0)
continue;
int uf=calc(ux, uy, ur);
ans=max(uf, ans);
double dE=f-uf;
if (dE<0||exp(-dE/t-sqrt(t))>frand())
x=ux,y=uy,r=ur,f=uf;
t*=k;
}
return ans;
}
int main()
{
srand((unsigned int)time(0));
n=read(),m=read(),R=read();
for (int i=1;i<=n;i++)
b[i].x=read(),b[i].y=read(),b[i].r=read();
for (int i=1;i<=m;i++)
{
e[i].x=read(),e[i].y=read();
maxp=max(maxp, abs(e[i].x)),maxp=max(maxp, abs(e[i].y));
}
maxp+=R;
int ans=0;
// while ((double)clock()/CLOCKS_PER_SEC < maxTime)
// ans=max(ans, stimulatedAnneal());
for (int i=1;i<=8;i++)
ans=max(ans, stimulatedAnneal());
printf("%d\n", ans);
system("pause");
return 0;
}
适用范围
模拟退火算法可以高效地求解\(NP\)完全问题，如货郎担问题\((\text{Travelling Salesman Problem,TSP})\)、最大截问题\((\text{Max Cut Problem})\)、0-1背包问题\((\text{Zero One Knapsack Problem})\)、图着色问题\((\text{Graph Colouring Problem})\)等等。
一些技巧
- 模拟退火的参数难以控制，不能保证一次退火就收敛到最优值，可以在一次结束后以当前温度在得到的解附近多次随机状态，尝试得到更优的解（其过程与模拟退火相似）。
- 可以在实际操作时二分参数来确定所设置的值。
- 可以根据数据的范围设定对应的温度参数。
- 当函数的峰过多时可以分块分别进行退火再在最终几个块的答案中取最优。
- 可以用时间取代温度的参数来强制退火时长，从而确保降温不会过快也不会超时。
参考内容
- 模拟退火-OI WIKI
- 深度学习 --- 模拟退火算法详解（Simulated Annealing， SA)
- Wiki-Stimulated Annealing
- Wiki-Adaptive Stimulated Annealing