模擬退火（速通教程）

速通模擬退火（Simulated Annealing）

作者：郭程([e-mail])

博客園、公眾號、github同步更新。

如果您喜歡我的作品，可以到我的github上給我一個star。也歡迎關注我的博客園。

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

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

唯一微信公眾號：

本教程為原創教程，轉載請標明出處，禁止商業用途。

本教程屬於快速入門（Quick-Start)系列，不對具體退火細節深究，只論算法競賽題目研究。

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

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

模擬退火是模擬物理上退火方法，通過N次迭代（退火），逼近函數的上的一個最值（最大或者最小值）。

比如逼近這個函數的最大值C點：

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\) ，不論是哪個點都可以，隨機！（不超過定義域就行）。

   這個點作為我們的初始值（相當於物體里的一個粒子）

   

2. 再找到一點\(f(x_0)\),來代表\(x_0\)所對應的函數值

   

3. 現在正式開始退火！

   剛才我們說了\(x_0\)相當於是一個粒子，所以我們會進行一個無序運動，也就是向左或者向右隨機移動

   是的，是隨機移動，可能向左，也可能向右，但是請記住一個關鍵點：移動的幅度和當前的溫度\(T\)有關。

   溫度\(T\)越大，移動的幅度越大。溫度\(T\)越小，移動的幅度就越小。這是在模擬粒子無序運動的狀態。

4. 接受(Accept)更"好"的狀態

   

   假設我們移動到了\(x_1\)處，那么這個點對應的\(f(x_1)\)很明顯答案是優於（大於）當前的\(f(x_0)\)的

   

   因此我們將答案進行更新。也就是將初始值進行替換：\(x_0=x_1,f(x_0)=f(x_1)\)。這是一種貪心的思想。

5. 以一定概率接受(Accept)更差的狀態

   這是退火最精彩的部分。

   為什么我們要接受一個更加差的狀態呢？因為可能在一個較差的狀態旁邊會出現一個更加高的山峰

   

   如果我們鼠目寸光，只盯着右半區，很容易隨着溫度的下降、左右跳轉幅度的減小而迷失自己，最后困死在小山丘中。

   而我們如果找到了左邊山峰的低點，以一定的概率接受了它（概率大小和溫度以及當前的值的關鍵程度有關），會在跳轉幅度減少之前，盡可能找到最優點。

   那么我們以多少的概率去接受它呢？我們用一個公式表示（這個公式我們只需記住，這是科學家推導出來的結論）：

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

   別慌！很簡單！我們來理解一下這里面的變量：

   1. \(e\)是自然對數，約等於2.71。我們可以把右上角這一坨值\(\frac{\Delta f}{kT}\)看成一個整體\(x\):

      \(e^x\)的圖形畫出來是這樣的：

      

      因為我們想要函數\(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))\)

      

      所以總結一下就是：

      • 隨機后的函數值如果結果更好，我們一定選擇它(即\(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 = (2*rand() - RAND_MAX) * T; 
    
    //讓x落在定義域內，如果沒在里面，就重新隨機。題目有要求需要寫，否則不用寫
    // ================
    while(x > ? || x < ? ...) {
        double dx = (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();
}