引言

在实际生活中，我们常常会遇到求函数最值的问题，那怎么办呢？我们当然可以选择爬山算法，即每次在当前最优解的附近选择一个解，如果它优于最优解，就接受它，否则不接受它，并调小选择范围，寻找下一个解。在某些情况下，它是适用的，比如下图


但这个算法的劣势非常明显——它会被局限在一个局部最优解上，无法取得全局最优解，比如下图这个函数。


这时，我们就可以使用一个玄学算法——模拟退火。

正文

简介

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

原理

模拟退火的原理也和金属退火的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。

——摘自百度百科（当然你也可以不看）

简单的说，模拟退火就是在一种一定范围内求多峰函数最值的算法。它在模拟温度降低的同时得出新解，温度越高，解的变化量越大，随着温度的逐渐降低，解的变化量也渐渐变小，并越发集中在最优解附近。最后温度达到了我们设置的最低温，对应到物理学上也就是结晶了，这时，我们可以认为当前我们取得的解就是最优解，当然也可能不是，整个算法也就终止了。

过程

我们先引入几个参数：当前最优解 E₀​ ，新解 E ，解变动量 ΔE （ E 与 E0​ 的差），上一个被接受的解 E_1​ ，初温 T₀​ ，末温*T_k*​ ，当前温度 T ，温度变动量 Δ ，再引用一张非常经典的图——


这张图非常好的展现了模拟退火的运行过程，从 T\_0*T*0​ 开始，每次乘上 ΔΔ 得到 T*T* ，如果 T*T* 小于 T\_k*Tk*​ 则终止降温。 T\_0*T*0​ 我一般设置在 10001000 ~ 50005000 左右， ΔΔ 则是一个略小于1的常数，而 T\_k*Tk*​ 一般设置在 1^{-8}1−8 到 1^{-15}1−15 之间(或者另外一个极小的数)。

在降温的同时，我们在 E1​ （不是最优解 E0​ ）的基础上扰动产生新解 E ，需要注意的是扰动大小随温度的降低而变小，因为在温度高的时候，解的变化量非常大，这时的任务是在全局范围中找到最优解的大致位置，随着温度的降低，解渐渐稳定，这时的任务是确定最优解的准确位置。

但每次得出新解以后，我们以什么原则，或者说什么概率来接受它呢？

这时就要用到Metropolis准则。简单说来，假设我们的目标是求最小值，如果 E 与 E0​ 的差，也就是 ΔE 小于 0 ，我们就接受当前解，因为它优于之前的最优解嘛。而如果 ΔE 大于 0 ，也就是我们遇到了一个更劣的解，我们也要以一定的概率来接受它，因为我们要找的一个多峰函数的全局最小值，因此就不能局限于一个局部的凹函数。而这个概率是 exp(−ΔE/T) 。

我个人对于这个概率的理解是这样的：对于 ΔE ，如果它较大，即我们遇到了一个劣得多的解，那我们接受它的概率就相对较小，因为 −ΔE 较小嘛；而如果 ΔE 较小，即我们遇到了一个较劣的解，我们接受它的概率就较大，因为 −ΔE 较大。对于 T ，随着时间的增加， T 变得越来越小，因此我们把 −ΔE 除以 T ，这样接受的概率就随着温度的降低而越来越小，因为 −ΔE 是一个负数嘛。而对于整个式子，当 T 较大的时候，我们会接受大部分的解，当 T 较小的时候，我们就只会接受 ΔE 较小的解。关于Metropolis准则的具体证明，过于玄学，这里就不给出了。当然你也可以自己试一下。如果选择接受 E ，则把 E₁​ 设置为 E ，然后降温并寻找下一个解。

这里再引用一张很糊的图：


到这里我们也就知道，模拟退火算法的速度和结果受参数（ *T*₀​ ， *T_k*​ ， Δ 还有随机数种子）的影响非常大，是一个玄学的算法，时间复杂度也是 *O*(玄学) 。

实现

例题

接下来我们结合一道例题来讲一讲模拟退火的c++​代码实现。UVA10228 A Star not a Tree? （这道题其实洛谷上也有）

英文题面尽管跳过，大意是给定 n 个点，求其费马点（到这 n 个点的距离最小的点）到所有点的距离和。此题各部分的代码实现都很方便，其实就是一道模板题，代码如下：

<pre class="EnlighterJSRAW" data-enlighter-group="" data-enlighter-highlight="" data-enlighter-language="cpp" data-enlighter-linenumbers="" data-enlighter-lineoffset="" data-enlighter-theme="" data-enlighter-title="">#include<iostream>
#include<cstdio>
#include<stdlib.h>
#include<iomanip>
#include<cmath>
#define R register
#define rep(i,a,b) for(R int i=a;i<=b;i++)
#define delta 0.996
#define maxn 50005
using namespace std;

inline int read() {
    int x=0,f=1;
    char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();}
    while('0'<=ch&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
    return x*f;
}

struct node{
    double x,y;
}poi[maxn];
int T,n;
double ansx,ansy,ax,ay,ans,t;

void clear() {
    ax=0,ay=0;
    ans=1e8;
}

double calculate(double x,double y) {
    double res=0;
    rep(i,1,n) {
        double dx=x-poi[i].x,dy=y-poi[i].y;
        res+=sqrt(dx*dx+dy*dy);
    }
    return res;
}

void simulate_anneal() {
    double x=ansx,y=ansy;
    t=3000;
    while(t>1e-15) {
        double X=x+((rand()<<1)-RAND_MAX)*t;
        double Y=y+((rand()<<1)-RAND_MAX)*t;
        double now=calculate(X,Y);
        double Delta=now-ans;
        if(Delta<0) {
            ansx=X,ansy=Y;
            x=X,y=Y;
            ans=now;
        } else if(exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y;
        t*=delta;
    }
}

void work() {
    ansx=ax/n,ansy=ay/n;
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
}

int main() {
    srand(1e9+7);
    T=read();
    rep(i,1,T) {
        n=read();
        clear();
        rep(j,1,n) {
            poi[j].x=read(),poi[j].y=read();
            ax+=poi[j].x,ay+=poi[j].y;
        }
        work();
        cout<<round(ans)<<'\n';
        if(i!=T) cout<<'\n';
    }
    return 0;
}

有几个注意点：坐标位置，温度和解变动量必须开成double​，一是为了确保精度，二是为了防止爆int。还要注意输出换行，实在很坑。

但是当你愉快的写玩此题并提交以后，可能会发现你并没有AC此题。记得之前说过的吗，我们得出不一定是最优解。这时候就涉及到一个麻烦的步骤——调参。通常有以下几种调参的方式：

1. 调大初温 T₀​ 。
2. 调小末温 T_k​ 。
3. 调大温度变动量 Δ 。
4. 选取一个新的随机数种子。
5. 多跑几遍模拟退火。
6. 开O2

其中第一，二点对于运行时间的影响不大。而第三点则非常关键，一个微调都会使运行时间和结果发生巨大变化。第五点也是一个有用的方式，一般我们跑三到五遍模拟退火，如果时间充裕，你也可以适当多跑一两百几遍。而第四点就非常看脸了，你当然可以选择某个恶臭的八位质数，但就我个人而言，最有用的还是这句随机数种子：

<pre class="EnlighterJSRAW" data-enlighter-group="" data-enlighter-highlight="" data-enlighter-language="cpp" data-enlighter-linenumbers="" data-enlighter-lineoffset="" data-enlighter-theme="" data-enlighter-title="">srand(time(0));

习题

学完一个新算法以后，当然应该练习啦。其实模拟退火的主过程基本就是模板了，唯一的麻烦点是对calculate()函数和接受概率的修改，比如下题：[JSOI2016]炸弹攻击1

此题的calculate()函数倒是很简单，麻烦的是修改接受概率。

题目要求的是最大值，那么 -ΔE−ΔE 就成了一个正数，怎么修改呢？其实此时我们只需把这句话：

<pre class="EnlighterJSRAW" data-enlighter-group="" data-enlighter-highlight="" data-enlighter-language="cpp" data-enlighter-linenumbers="" data-enlighter-lineoffset="" data-enlighter-theme="" data-enlighter-title="">else if(exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y;

中的 >> 号改为 << 号就可以了，如下：

<pre class="EnlighterJSRAW" data-enlighter-group="" data-enlighter-highlight="" data-enlighter-language="cpp" data-enlighter-linenumbers="" data-enlighter-lineoffset="" data-enlighter-theme="" data-enlighter-title="">else if(exp(-Delta/t)*RAND_MAX<rand()) x=X,y=Y;

而这样就很玄学了。之前我说错了，因为 −ΔE 成了一个正数，所以 exp(−Delta/t) 必定是大于1的，也就是没有接受劣解的概率。而此题 ΔE 波动小，搜寻范围大，所以我们这样写就可以手动避免算法陷入劣解不能自拔。但这样写的原因是我过了此题以后才想出来的。

代码就略过了，实在很简单。

---

模拟退火的应用不仅仅是求点坐标，还可以拿来求序列。其实过程也很简单，每次随机交换序列中的两个元素就可以了，而对于网格，看作是二维序列即可。下面有一道求序列的题目：[SCOI2008]城堡

读完题以后，你可能不知道此题和序列有何关系。但我们其实可以这样考虑：把所有没有城堡的城市抽象成一个序列，而序列的前 k 个城市，就是要修建城堡的城市。

而关于calculate()函数，我们可以先用floyd​算法预处理出每个城市之间的距离，在这个函数中我们只需 n² 扫描一次，求出所有城市中离最近城堡的距离的最大值就可以了。代码如下：

<pre class="EnlighterJSRAW" data-enlighter-group="" data-enlighter-highlight="" data-enlighter-language="cpp" data-enlighter-linenumbers="" data-enlighter-lineoffset="" data-enlighter-theme="" data-enlighter-title="">#include<iostream>
#include<cstdio>
#include<stdlib.h>
#include<cmath>
#include<ctime>
#define rep(i,a,b) for(register int i=a;i<=b;i++)
#define maxn 500
#define inf 0x3f3f3f3f
#define delta 0.996
using namespace std;

inline int read() {
    int f=1,x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();}
    while('0'<=ch&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
    return x*f;
} 

void write(int x) {
    if(x<0) x=-x,putchar('-');
    if(x>9) write(x/10);
    putchar(x%10+'0');
}

struct edge{
    int a,b,next,v;
}e[maxn];
int head[maxn],cnt,n,m,k,v[maxn],cas[maxn];
int dis[maxn][maxn],p[maxn],N,X[maxn],ans=1e8;
double t;

int calculate(int x[]) {
    rep(i,1,k) cas[x[i]]=1;
    int res=-inf;
    rep(i,0,n) {
        int minn=inf;
        rep(j,0,n) if(cas[j]) minn=min(minn,dis[i][j]);
        res=max(res,minn);
    }
    rep(i,1,k) cas[x[i]]=0;
    return res;
}

void simulate_anneal() {
    int a[maxn];
    rep(i,1,N) a[i]=p[i];
    t=5000;
    while(t>1e-15) {
        int b[maxn];
        rep(i,1,N) b[i]=a[i];
        int x=rand()%N+1;
        int y=rand()%N+1;
        swap(b[x],b[y]);
        int now=calculate(b);
        double Delta=now-ans;
        if(Delta<0) {
            ans=now;
            rep(i,1,N) p[i]=a[i]=b[i];
        } else if(exp(-Delta/t)*RAND_MAX>rand()) {
            rep(i,1,N) a[i]=b[i];
        }
        t*=delta;
    }
}

void work() {
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
}

int main() {
    srand(time(0));
    n=read(),m=read(),k=read();
    n--;
    rep(i,0,n) X[i]=read();
    rep(i,0,n) v[i]=read();
    rep(i,1,m) {
        int a=read();
        cas[a]=1;
    }
    rep(i,0,n) if(!cas[i]) p[++N]=i;
    rep(i,0,n) rep(j,0,n) dis[i][j]=inf;
    rep(i,0,n) dis[i][X[i]]=dis[X[i]][i]=min(v[i],dis[i][X[i]]),dis[i][i]=0;
    rep(c,0,n) rep(i,0,n) rep(j,0,n)
    dis[i][j]=min(dis[i][j],dis[i][c]+dis[c][j]);
    work();
    write(ans);
    return 0;
}

最后

要做好调参的心理准备，我在三天之内调参交了七页。

推荐几题：

1. [JSOI2008]球形空间产生器
2. [CEOI2004]锯木厂选址
3. Coloring

其他的习题自己找去吧。

update 2020.3.3 加入了 LATE​X 数学公式渲染，并添加了一张图。

update 2020.5.1 修锅，感谢@M_sea 纠错。

update 2020.6.21 修锅，小部分改动,致谢@Caro23333。