模拟退火学习笔记
- 本文基于多篇博客编写，仅供学习使用，将会存在多处对其他文章内容的直接复制
模拟退火是一种随机化算法。当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，我们常使用模拟退火求解。用一句话概括：如果新状态的解更优则修改答案，否则以一定概率接受新状态。
——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<T_k\)时退火结束，当前的最优解（注意不是当前解，因为有概率导致解不够好却转移）即为所求的答案。
我们可以从定义域中任取一个点开始降温，随机向某一个点进行跳跃（具体看代码）。
例题
Strange Fuction
输入高次函数的一次项系数，求该函数的最小值。
代码
[click]
#include <cstdio>
#include <cctype>
#include <cstdlib>
#include <ctime>
#include <cmath>
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 x<y?x:y;}
double frand() {return (double)rand()/RAND_MAX;}//随机获得一个[0, 1]内的浮点数
double calc(double x)//计算目标函数
{
return 6.0*pow(x, 7)+8.0*pow(x, 6)+7.0*pow(x, 3)+5.0*x*x-y*x;
}
double stimulatedAnneal()
{
double x=frand()*100,y=calc(x);//随机选一个当前状态
double t=T_0,ans=y;
//设置初始状态
while(t>T_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 <cstdio>
#include <cctype>
#include <cmath>
#include <cstdlib>
#include <ctime>
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;
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 x<y?x:y;}
int read()
{
int res=0,p=1;
char ch=getchar();
while(!isdigit(ch))
{
if (ch=='-')
p=-1;
ch=getchar();
}
while(isdigit(ch))
res=res*10+ch-'0',ch=getchar();
return res*p;
}
double frand() {return (double)rand()/RAND_MAX;}//随机获得一个[0, 1]上的浮点数
double getDis(double x,double y,double a,double b) {return sqrt((a-x)*(a-x)+(b-y)*(b-y));}//返回欧拉距离
double getRadius(double x,double y)//计算坐标对应最大的轰炸半径
{
double r=R;
for (int i=1;i<=n;i++)
r=min(r, getDis(x, y, b[i].x, b[i].y)-b[i].r);
return r;
}
bool check(double x,double y) {return fabs(x)<=maxp&&fabs(y)<=maxp;}//一个小的剪枝优化
int calc(double x,double y,double r)
{
int res=0;
for (int i=1;i<=m;i++)
res+=(getDis(x, y, e[i].x, e[i].y)<=r);
return res;
}
int stimulatedAnneal()
{
double t=2*maxp;//根据数据范围而改变的初温
double x=maxp*(frand()*2-1),y=maxp*(frand()*2-1),r=getRadius(x, y);//随机获得初态
int f=calc(x, y, r),ans=f;
while(t>T_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})\)等等。
一些技巧
- 模拟退火的参数难以控制，不能保证一次退火就收敛到最优值，可以在一次结束后以当前温度在得到的解附近多次随机状态，尝试得到更优的解（其过程与模拟退火相似）。（对于函数值是离散的情况可能作用不大）
- 可以在实际操作时二分参数来确定所设置的值。
- 可以根据数据的范围设定对应的温度参数。
- 当函数的峰过多时可以分块分别进行退火再在最终几个块的答案中取最优。
- 可以用时间取代降火的次数来确保在不超时的情况下尽可能多地降火，增大获得正确答案的概率。
