时间戳

首页 > 技术文章 > 模拟退火（速通教程）

模拟退火（速通教程）

autoloop 2021-08-21 15:13 原文

速通模拟退火（Simulated Annealing）

作者：郭程(judgemeido@163.com)
博客园、公众号、github同步更新。

如果您喜欢我的作品，可以到我的github上给我一个star。也欢迎关注我的博客园。

github地址：https://github.com/AutoLoop255/acm

博客园地址：https://www.cnblogs.com/autoloop

唯一微信公众号：

公众号二维码

本教程为原创教程，转载请标明出处，禁止商业用途。

本教程属于快速入门（Quick-Start)系列，不对具体退火细节深究，只论算法竞赛题目研究。

1.是什么？为什么？怎么做？

Q1:模拟退火是什么算法？

模拟退火是模拟物理上退火方法，通过N次迭代（退火），逼近函数的上的一个最值（最大或者最小值）。

比如逼近这个函数的最大值C点：

image-20210820234047376

Q2:模拟退火为什么可行？

讨论这个问题需要理解一下物理原型是怎么样的，也就是原来是怎么“退火”的：

模拟退火算法的思想借鉴于固体的退火原理，当固体的温度很高的时候，内能比较大，固体的内部粒子处于快速无序运动，当温度慢慢降低的过程中，固体的内能减小，粒子的慢慢趋于有序，最终，当固体处于常温时，内能达到最小，此时，粒子最为稳定。

注意标粗字体：

  1. 温度高->运动速度快（温度低->运动速度慢）
  2. 温度是缓慢（想象成特别慢的那种）降低的
  3. 温度基本不再变化后，趋于有序(最后内能达到最小，也就是接近最优)

我们通过模拟这个操作，使得我们需要的答案“趋于有序”，也就是最靠近需要的值（最值）。

Q3:怎么做？

大方向

首先，理解一下大方向：

模拟退火就是一种循环算法。

  1. 我们先设定一个初始的温度\(T\)（这个温度会比较高，比如2000）
  2. 每次循环都退火一次。（具体怎么操作后面详解）
  3. 然后降低\(T\)的温度，我们通过让\(T\)和一个“降温系数”\(\Delta T\)（一个接近1的小数，比如\(0.99\))相乘，达到慢慢降低温度的效果，直到接近于0（我们用\(eps\)来代表一个接近0的数(比如0.00001)，只要\(T<eps\)就可以退出循环了）

所以总的来说，用伪代码表示退火的流程是这样的：

double T = 2000; //代表开始的温度
double dT = 0.99; //代表系数delta T
double eps = 1e-14; //相当于0.0000000000000001
while(T > eps) {
    //--------------
    //这里是每一次退火的操作
	//--------------
    T = T * dT; //温度每次下降一点点， T * 0.99
}

退火详解

我们要求解的答案无非是两个：自变量\(x\)和对应函数的最大值\(f(x)\)

那么模拟退火的是怎么做到的呢？【下面车速很快，请系好安全带！】

  1. 我们先随机找一点\(x_0\) ，不论是哪个点都可以，随机！（不超过定义域就行）。

    这个点作为我们的初始值（相当于物体里的一个粒子）

    image-20210821003424926

  2. 再找到一点\(f(x_0)\),来代表\(x_0\)所对应的函数值

    image-20210821003548607

  3. 现在正式开始退火！

    刚才我们说了\(x_0\)相当于是一个粒子，所以我们会进行一个无序运动，也就是向左或者向右随机移动

    是的，是随机移动，可能向左，也可能向右，但是请记住一个关键点：移动的幅度和当前的温度\(T\)有关。

    温度\(T\)越大，移动的幅度越大。温度\(T\)越小，移动的幅度就越小。这是在模拟粒子无序运动的状态。

  4. 接受(Accept)更"好"的状态

    image-20210821004445418

    假设我们移动到了\(x_1\)处，那么这个点对应的\(f(x_1)\)很明显答案是优于（大于）当前的\(f(x_0)\)的

    image-20210821005349066

    因此我们将答案进行更新。也就是将初始值进行替换：\(x_0=x_1,f(x_0)=f(x_1)\)。这是一种贪心的思想。

  5. 以一定概率接受(Accept)更差的状态

    这是退火最精彩的部分。

    为什么我们要接受一个更加差的状态呢？因为可能在一个较差的状态旁边会出现一个更加高的山峰

    image-20210821010525132

    如果我们鼠目寸光，只盯着右半区，很容易随着温度的下降、左右跳转幅度的减小而迷失自己，最后困死在小山丘中。

    而我们如果找到了左边山峰的低点，以一定的概率接受了它（概率大小和温度以及当前的值的关键程度有关），会在跳转幅度减少之前，尽可能找到最优点。

    那么我们以多少的概率去接受它呢？我们用一个公式表示（这个公式我们只需记住，这是科学家推导出来的结论）：

    \[\Huge e^{\frac{\Delta f}{kT}} \]

    别慌！很简单！我们来理解一下这里面的变量：

    1. \(e\)是自然对数，约等于2.71。我们可以把右上角这一坨值\(\frac{\Delta f}{kT}\)看成一个整体\(x\):

      \(e^x\)的图形画出来是这样的：

      image-20210821012147778

      因为我们想要函数\(e^x\)来代表一个概率值，所以我们只需要关注x为负数的部分即可：

      负数部分的值域是在\((0,1)\)开区间内，x越小，越接近0，越大越靠近1。

      因为在0到1之间，所以这个值相当于是概率了。比如\(e^x=0.97\)，那么我们接受的概率就是\(97\%\)

      而正数部分的值域会大于1，也就是说概率会超过\(100\%\)，所以会一定选（其实是上一种找到更优的情况）

    2. \(kT\)

      \(k\)其实是个物理学常数，我们在代码中不会用到。

      \(T\)很简单，就是当前的温度。所以实际上这个分母就是\(T\)，\(k\)当做\(1\)使用。

    3. \(\Delta f\)

      我们着重讲一下什么是\(\Delta f\)。

      其实从前面的函数\(e^x\)中可以发现，\(\Delta f\)必须是个负数！

      我们想要函数\(e^x\)来代表一个概率值，一定要让它的值域属于\((0,1)\)，所以\(\frac{\Delta f}{kT}\)必须是个负数。但是\(kT\)在我们的模拟中一定是正数，那么\(\Delta f\)必须是个负数！

      其实\(\Delta f\)就是当前解的函数值与目标解函数值之差，$\Delta f= - \left | f(x_0)-f(x_1) \right | $，并且一定是个负数。这个需要具体问题具体分析。

      比如现在我们求一个函数的最大值，那么如果\(f(x_0) < f(x_1)\)了，那就说明结果变好了，我们肯定选择它（见第4点）

      如果\(f(x_0) > f(x_1)\)，那就说明结果变差了，我们需要概率选择它，因此\(\Delta f=-(f(x_0) - f(x_1))\)

      image-20210821005349066

      所以总结一下就是：

      • 随机后的函数值如果结果更好，我们一定选择它(即\(x_0=x_1,f(x_0)=f(x_1)\))
      • 随机后的函数值如果结果更差，我们以\(\LARGE e^{\frac{\Delta f}{kT}}\)的概率接受它

2021年8月21日注：

关于接受概率\(e^x\)和\(\Delta f\)以及\(kT\)的关系：

当\(kT\)越大时(温度越高)，由于\(\Delta f\)是一个负数，所以算出来的值会越大。比如-1 / 1000 等于 -0.001，很明显-0.001 > -1，由于\(e^x\)是个单调递增函数，所以温度越高时，接受的概率就越大。

\(\Delta f\)越小，说明答案越糟糕，那么接受的概率就越小。

伪代码流程

注：对代码中的函数作出解释：

①对rand()函数

  1. rand()函数可以默认拿到[0,32767]内的随机整数
  2. RAND_MAX = 32767，可以看作常量。本质是宏定义: #define RAND_MAX 32767
  3. rand() * 2 的范围是[0,32767 * 2]
  4. rand() * 2 - RAND_MAX 的范围是[-32767, 32767]

②对exp()函数

  1. exp(x)代表\(e^x\)

③关于exp((df - f) / T) * RAND_MAX > rand()

  1. 目的是要概率接受，但是\(e^x\)是个准确值，所以从理论上我们可以生成一个（0,1）的随机数，如果\(e^x\)比(0,1)这个随机数要大，那么我们就接受。
  2. 但是由于rand()只能产生[0,32767]内的随机整数，化成小数太过麻烦。所以我们可以把左边乘以RAND_MAX（也就是把概率同时扩大32767倍），效果就等同于\(e^x\)比(0,1)了。
double T = 2000; //代表开始的温度
double dT = 0.99; //代表系数delta T
double eps = 1e-14; //相当于0.0000000000000001

//用自变量计算函数值,这里可能存在多个自变量对应一个函数值的情况，比如f(x,y)
double func(int x, ... ) {
    //这里是对函数值进行计算
    double ans = .......
    return ans;
}
//原始值
double x = rand(); //x0取随机值
double f = func(x,...); //通过自变量算出f(x0)的值
while(T > eps) {
    //--------------
    //这里是每一次退火的操作
    
    //x1可以左右随机移动，幅度和温度T正相关，所以*T
    //注意这里移动可以左右移动，但是也可以单向移动
    //关于rand()详细见开头注的①
    double dx = x + (2*rand() - RAND_MAX) * T; 
    
    //让x落在定义域内，如果没在里面，就重新随机。题目有要求需要写，否则不用写
    // ================
    while(x > ? || x < ? ...) {
        double dx = x + (2*rand() - RAND_MAX) * T; 
    }
    // ================
    
    //求出f(x1)的值
    double df = func(dx);
    //这里需要具体问题具体分析，我们要接受更加优秀的情况。可能是df < f(比如求最小值）
    if(f < df) {
        f = df; x = dx;  [...,y = dy;] // 接受，替换值，如果多个自变量，那么都替换
    }
    //否则概率接受，注意这里df-f也要具体问题具体分析。
    //详细见开头注的②③
    else if(exp((df - f) / T) * RAND_MAX > rand()) {
        f = df; x = dx;  [...y = dy;] // 接受，替换值，如果多个自变量，那么都替换
    }
	//--------------
    T = T * dT; //温度每次下降一点点， T * 0.99
}
//最后输出靠近最优的自变量x值，和函数值f(x)
cout << x << " " << f << endl;

2.通过模拟退火算出\(\sqrt{n}\)的值

思路是这样的：我们试图通过退火找出一个值\(x_0\)，使得\({x_0}^2\)的值更加接近于\({\sqrt{n}}^2\)。（记住退火是让一个随机值去逼近最后的答案）

因为\({x_0}^2\)的值更加接近于\({\sqrt{n}}^2\)。因此\(x_0\)值就更加接近于\(\sqrt{n}\)

  1. 所以我们需要一个函数\(f(x)\)，算出\({x_0}^2\)和\({\sqrt{n}}^2\)的接近程度，那么毋庸置疑，我们需要算出他们绝对值的差。

    \[f(x)=\left |x^2 - n \right | \]

    也就是说我们的函数表达式就有了

    //n代表我们最后函数要逼近的值
double n;
//x表示我们随机产生的那个数的平方和n的靠近程度
double func(double x) {
    return fabs(x * x - n);
}

  2. 写出退火函数SA()

    double T = 20000; //初始温度，初始温度主要是看题目开始需要跳转的幅度。
double dT = 0.993; //变化率，这里需要速度稍微慢一点，写0.995 或者 0.997都可以，但是越靠近1速度就越慢 
const double eps = 1e-14; //10的-14次方已经非常小了，写这个大概率没问题
void SA() {
    //首先随机生成一个点x0，这里我用0代替。
    double x = 0;
    //算出x平方和n的差距f(x0)
    double f = func(x);
    while(T > eps) {
        //这里x0既可以变小，也可以变大，所以我们正负都要进行一个跳转,算出变换后的点dx
        double dx = x + (2 * rand() - RAND_MAX) * T;
        //但是请注意，dx很明显要保证 >= 0才行，因为算术平方根的定义域是>=0，因此小于0就重新随机
        while(dx < 0) dx = x + (2 * rand() - RAND_MAX) * T;
        //算出变换后的点dx的平方和n的差距，f(dx)
        double df = func(dx);
        //这里就是关键的地方了，很明显我们需要算出来的function值越小，自变量x更加接近那个根号值。
        //所以如果新来的值df 比 f更小，我们百分百接受
        if(df < f) {
            //注意更新所有变量
            f = df; x = dx;
        }
        //否则我们概率接受，这里的需要写 f - df了，因为这样才是负值。负值说明我们并不是贪心接受的，他是不太好的值。
        else if(exp((f - df)/T) * RAND_MAX > rand()) {
            //注意更新所有变量
            f = df; x = dx;
        }
        //温度下降一下
        T *= dT;
    }
    printf("%.8lf",x);
}

最后贴上完整代码和注释供大家调试。

#include <bits/stdc++.h>
using namespace std;
//n代表我们最后函数要逼近的值
double n;
//x表示我们随机产生的那个数的平方和n的靠近程度
double func(double x) {
    return fabs(x * x - n);
}
double T = 20000; //初始温度，初始温度主要是看题目开始需要跳转的幅度。
double dT = 0.993; //变化率，这里需要速度稍微慢一点，写0.995 或者 0.997都可以，但是越靠近1速度就越慢 
const double eps = 1e-14; //10的-14次方已经非常小了，写这个大概率没问题
void SA() {
    //首先随机生成一个点x0，这里我用0代替。
    double x = 0;
    //算出x平方和n的差距f(x0)
    double f = func(x);
    while(T > eps) {
        //这里x0既可以变小，也可以变大，所以我们正负都要进行一个跳转,算出变换后的点dx
        double dx = x + (2 * rand() - RAND_MAX) * T;
        //但是请注意，dx很明显要保证 >= 0才行，因为算术平方根的定义域是>=0，因此小于0就重新随机
        while(dx < 0) dx = x + (2 * rand() - RAND_MAX) * T;
        //算出变换后的点dx的平方和n的差距，f(dx)
        double df = func(dx);
        //这里就是关键的地方了，很明显我们需要算出来的function值越小，自变量x更加接近那个根号值。
        //所以如果新来的值df 比 f更小，我们百分百接受
        if(df < f) {
            //注意更新所有变量
            f = df; x = dx;
        }
        //否则我们概率接受，这里的需要写 f - df了，因为这样才是负值。负值说明我们并不是贪心接受的，他是不太好的值。
        else if(exp((f - df)/T) * RAND_MAX > rand()) {
            //注意更新所有变量
            f = df; x = dx;
        }
        //温度下降一下
        T *= dT;
    }
    printf("%.8lf",x);
}
int main() 
{
    cin >> n;
    SA();
}

3.例题POJ - 2420

题目描述

给出平面上N（N<=100）个点，你需要找到一个这样的点，使得这个点到N个点的距离之和尽可能小。输出这个最小的距离和（四舍五入到最近的整数）。

Input输入

第一行N，接下来N行每行两个整数，表示N个点

Output输出

一行一个正整数，最小的距离和。

Sample Input样例输入

4
0 0
0 10000
10000 10000
10000 0

Sample Output样例输出

28284

题目解析

学过平面几何的一定知道，假定两点的坐标为\(a(x_1,y_1),b(x_2,y_2)\)，那么他们之间的距离\(d\)可以用勾股定理计算：

\[d = \sqrt{(x_1-x_2)^2-(y_1-y_2)^2} \]

所以我们的函数\(func()\)就是：求出一个随机点\(A=(x_0,y_0)\)到\(N\)个点的距离之和\(D\)，就是把这个点A和所有N个点的距离相加，即：

\[\large D=\sum_{i=1}^{N} {\sqrt{(x_0-x_i)^2-(y_0-y_i)^2}} \]

我们将N个点用结构体存起来

const int N = 1e5 + 5;
struct Point {
    double x, y;
}e[N];

//main函数中输入N个点
int main() 
{
    cin >> n;
    for(int i = 1; i <= n; i++) 
        scanf("%lf%lf", &e[i].x, &e[i].y);
}

那么我们最关键的函数（求出一个随机点\(A=(x_0,y_0)\)到\(N\)个点的距离之和\(D\)）就是：

//输入我们随机生成的点A 坐标 x0, y0
double func(double x, double y) {
    double ans = 0;
    for(int i = 1; i <= n; i++) {
        double xx = e[i].x, yy = e[i].y; // xx 就是 xi ，yy 就是 yi
        double p = fabs(e[i].x - x); // x0 - xi
        double q = fabs(e[i].y - y); // y0 - yi
        ans += (sqrt((p * p + q * q))); // ans加上到随机点A到点i的距离
    }
    return ans;
}

然后套用上面的公式，我们写出SA()模拟退火

double T = 2000; //初始温度
double dT = 0.993; //变化率，这里需要速度稍微慢一点，写0.995 或者 0.997都可以，但是越靠近1速度就越慢 
const double eps = 1e-14; //10的-14次方已经非常小了，写这个大概率没问题
void SA() {
    //随机生成点A(x0,y0)，这里我直接赋值为了0。当然也可以写rand()，差别都不大。
    double x = 0, y = 0;
    //生成的点A(x0,y0)对应的函数值f(x0, y0)
    double f = func(x, y);
    while(T > eps) {
        //这里对x和y同时进行一个偏移，很明显在一个平面中，上下左右都可以移动，所以我们选择了rand() * 2 - RAND_MAX
        //对于每个点，我们的偏移幅度都要 * T,温度越高，偏移量越大
        double dx = x + (rand() * 2 - RAND_MAX) * T;
        double dy = y + (rand() * 2 - RAND_MAX) * T;
        //算出偏移后的点B(dx, dy)对应的函数值f(dx, dy)
        double df = func(dx, dy);
        //这里就是关键的地方了，很明显我们需要算出来的function值越小，就更加优秀
        //所以如果新来的值df 比 f更小，我们百分百接受
        if(df < f) {
            //注意更新所有变量
            f = df; y = dy; x = dx;
        }
        //否则我们概率接受，这里的需要写 f - df了，因为这样才是负值
        else if(exp((f - df) / T) * RAND_MAX > rand()) {
            //注意更新所有变量
            f = df; y = dy; x = dx;
        }
        //降低温度
        T = T * dT;
    }
    //输出答案
    printf("%.f\n", f);
}

最后贴上完全版代码

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <ctime>
using namespace std;
const int N = 1e5 + 5;
double dT = 0.993;
double T = 2000;
const double eps = 1e-14;  
struct vega {
    double x, y;
}e[N];
int n;
//输入我们随机生成的点A 坐标 x0, y0
double func(double x, double y) {
    double ans = 0;
    for(int i = 1; i <= n; i++) {
        double xx = e[i].x, yy = e[i].y; // xx 就是 xi ，yy 就是 yi
        double p = fabs(e[i].x - x); // x0 - xi
        double q = fabs(e[i].y - y); // y0 - yi
        ans += (sqrt((p * p + q * q))); // ans加上到随机点A到点i的距离
    }
    return ans;
}
void SA() {
    //随机生成点A(x0,y0)，这里我直接赋值为了0。当然也可以写rand()，差别都不大。
    double x = 0, y = 0;
    //生成的点A(x0,y0)对应的函数值f(x0, y0)
    double f = func(x, y);
    while(T > eps) {
        //这里对x和y同时进行一个偏移，很明显在一个平面中，上下左右都可以移动，所以我们选择了rand() * 2 - RAND_MAX
        //对于每个点，我们的偏移幅度都要 * T,温度越高，偏移量越大
        double dx = x + (rand() * 2 - RAND_MAX) * T;
        double dy = y + (rand() * 2 - RAND_MAX) * T;
        //算出偏移后的点B(dx, dy)对应的函数值f(dx, dy)
        double df = func(dx, dy);
        //这里就是关键的地方了，很明显我们需要算出来的function值越小，就更加优秀
        //所以如果新来的值df 比 f更小，我们百分百接受
        if(df < f) {
            //注意更新所有变量
            f = df; y = dy; x = dx;
        }
        //否则我们概率接受，这里的需要写 f - df了，因为这样才是负值
        else if(exp((f - df) / T) * RAND_MAX > rand()) {
            //注意更新所有变量
            f = df; y = dy; x = dx;
        }
        //降低温度
        T = T * dT;
    }
    //输出答案
    printf("%.f\n", f);
}
int main()
{
    cin >> n;
    for(int i = 1; i <= n; i++) {
        scanf("%lf%lf", &e[i].x, &e[i].y);
    }
    SA();
}

推荐阅读