著名的模擬退火演算法，它是一種基於蒙特卡洛思想設計的近似求解最優化問題的方法。

一點歷史——如果你不感興趣，可以跳過

美國物理學家 N.Metropolis 和同仁在1953年發表研究複雜系統、計算其中能量分佈的文章，他們使用蒙特卡羅模擬法計算多分子系統中分子的能量分佈。這相當於是本文所探討之問題的開始，事實上，模擬退火中常常被提到的一個名詞就是Metropolis準則，後面我們還會介紹。

美國IBM公司物理學家 S.Kirkpatrick、C. D. Gelatt 和 M. P. Vecchi 於1983年在《Science》上發表了一篇頗具影響力的文章：《以模擬退火法進行最優化（Optimization by Simulated Annealing）》。他們借用了Metropolis等人的方法探討一種旋轉玻璃態系統（spin glass system）時，發覺其物理系統的能量和一些組合最優（combinatorial optimization）問題（著名的旅行推銷員問題TSP即是一個代表例子）的成本函式相當類似：尋求最低成本即似尋求最低能量。由此，他們發展出以 Metropolis  方法為本的一套演算法，並用其來解決組合問題等的尋求最優解。

幾乎同時，歐洲物理學家 V.Carny 也發表了幾乎相同的成果，但兩者是各自獨立發現的；只是Carny“運氣不佳”，當時沒什麼人注意到他的大作；或許可以說，《Science》雜誌行銷全球，“曝光度”很高，素負盛名，而Carny卻在另外一本發行量很小的專門學術期刊《J.Opt.Theory Appl.》發表其成果因而並未引起應有的關注。

Kirkpatrick等人受到Metropolis等人用蒙特卡羅模擬的啟發而發明了“模擬退火”這個名詞，因為它和物體退火過程相類似。尋找問題的最優解（最值）即類似尋找系統的最低能量。因此係統降溫時，能量也逐漸下降，而同樣意義地，問題的解也“下降”到最值。

一、什麼是退火——物理上的由來

在熱力學上，退火（annealing）現象指物體逐漸降溫的物理現象，溫度愈低，物體的能量狀態會低；夠低後，液體開始冷凝與結晶，在結晶狀態時，系統的能量狀態最低。大自然在緩慢降溫（亦即，退火）時，可“找到”最低能量狀態：結晶。但是，如果過程過急過快，快速降溫（亦稱「淬鍊」，quenching）時，會導致不是最低能態的非晶形。

如下圖所示，首先（左圖）物體處於非晶體狀態。我們將固體加溫至充分高（中圖），再讓其徐徐冷卻，也就退火（右圖）。加溫時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最後在常溫時達到基態，內能減為最小（此時物體以晶體形態呈現）。

似乎，大自然知道慢工出細活：緩緩降溫，使得物體分子在每一溫度時，能夠有足夠時間找到安頓位置，則逐漸地，到最後可得到最低能態，系統最安穩。

二、模擬退火（Simulate Anneal）

如果你對退火的物理意義還是暈暈的，沒關係我們還有更為簡單的理解方式。想象一下如果我們現在有下面這樣一個函式，現在想求函式的（全域性）最優解。如果採用Greedy策略，那麼從A點開始試探，如果函式值繼續減少，那麼試探過程就會繼續。而當到達點B時，顯然我們的探求過程就結束了（因為無論朝哪個方向努力，結果只會越來越大）。最終我們只能找打一個區域性最後解B。

模擬退火其實也是一種Greedy演算法，但是它的搜尋過程引入了隨機因素。模擬退火演算法以一定的概率來接受一個比當前解要差的解，因此有可能會跳出這個區域性的最優解，達到全域性的最優解。以上圖為例，模擬退火演算法在搜尋到區域性最優解B後，會以一定的概率接受向右繼續移動。也許經過幾次這樣的不是區域性最優的移動後會到達B 和C之間的峰點，於是就跳出了局部最小值B。

根據Metropolis準則，粒子在溫度T時趨於平衡的概率為exp(-ΔE/(kT))，其中E為溫度T時的內能，ΔE為其改變數,k為Boltzmann常數。Metropolis準則常表示為

Metropolis準則表明，在溫度為T時，出現能量差為dE的降溫的概率為P(dE)，表示為：P(dE) = exp( dE/(kT) )。其中k是一個常數，exp表示自然指數，且dE<0。所以P和T正相關。這條公式就表示：溫度越高，出現一次能量差為dE的降溫的概率就越大；溫度越低，則出現降溫的概率就越小。又由於dE總是小於0（因為退火的過程是溫度逐漸下降的過程），因此dE/kT < 0 ，所以P(dE)的函式取值範圍是(0,1) 。隨著溫度T的降低，P(dE)會逐漸降低。 我們將一次向較差解的移動看做一次溫度跳變過程，我們以概率P(dE)來接受這樣的移動。也就是說，在用固體退火模擬組合優化問題，將內能E模擬為目標函式值 f，溫度T演化成控制引數 t，即得到解組合優化問題的模擬退火演演算法：由初始解 i 和控制引數初值 t 開始，對當前解重複“產生新解→計算目標函式差→接受或丟棄”的迭代，並逐步衰減 t 值，演算法終止時的當前解即為所得近似最優解，這是基於蒙特卡羅迭代求解法的一種啟發式隨機搜尋過程。退火過程由冷卻進度表(Cooling Schedule)控制，包括控制引數的初值 t 及其衰減因子Δt 、每個 t 值時的迭代次數L和停止條件S。

總結起來就是：

• 若f( Y(i+1) ) <= f( Y(i) )  (即移動後得到更優解)，則總是接受該移動；
• 若f( Y(i+1) ) > f( Y(i) )  (即移動後的解比當前解要差)，則以一定的概率接受移動，而且這個概率隨著時間推移逐漸降低（逐漸降低才能趨向穩定）相當於上圖中，從B移向BC之間的小波峰時，每次右移（即接受一個更糟糕值）的概率在逐漸降低。如果這個坡特別長，那麼很有可能最終我們並不會翻過這個坡。如果它不太長，這很有可能會翻過它，這取決於衰減 t 值的設定。

關於普通Greedy演算法與模擬退火，有一個有趣的比喻：

• • 普通Greedy演算法：兔子朝著比現在低的地方跳去。它找到了不遠處的最低的山谷。但是這座山谷不一定最低的。這就是普通Greedy演算法，它不能保證區域性最優值就是全域性最優值。
  • 模擬退火：兔子喝醉了。它隨機地跳了很長時間。這期間，它可能走向低處，也可能踏入平地。但是，它漸漸清醒了並朝最低的方向跳去。這就是模擬退火。

      通過一個例項來程式設計演示模擬退火的執行。特別地，我們這裡所採用的例項是著名的“旅行商問題”（TSP，Traveling Salesman Problem），它是哈密爾頓迴路的一個例項化問題，也是最早被提出的NP問題之一。

     TSP是一個最常被用來解釋模擬退火用法的問題，因為這個問題比較有名，我們這裡不贅言重述，下面直接給出C++實現的程式碼：

#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>

#define N     30      //城市數量
#define T     3000    //初始溫度
#define EPS   1e-8    //終止溫度
#define DELTA 0.98    //溫度衰減率

#define LIMIT 1000   //概率選擇上限
#define OLOOP 20    //外迴圈次數
#define ILOOP 100   //內迴圈次數

using namespace std;

//定義路線結構體
struct Path
{
    int citys[N];
    double len;
};

//定義城市點座標
struct Point
{
    double x, y;
};

Path bestPath;        //記錄最優路徑
Point p[N];       //每個城市的座標
double w[N][N];   //兩兩城市之間路徑長度
int nCase;        //測試次數

double dist(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

void GetDist(Point p[], int n)
{
    for(int i = 0; i < n; i++)
        for(int j = i + 1; j < n; j++)
            w[i][j] = w[j][i] = dist(p[i], p[j]);
}

void Input(Point p[], int &n)
{
    scanf("%d", &n);
    for(int i = 0; i < n; i++)
        scanf("%lf %lf", &p[i].x, &p[i].y);
}

void Init(int n)
{
    nCase = 0;
    bestPath.len = 0;
    for(int i = 0; i < n; i++)
    {
        bestPath.citys[i] = i;
        if(i != n - 1)
        {
            printf("%d--->", i);
            bestPath.len += w[i][i + 1];
        }
        else
            printf("%d\n", i);
    }
    printf("\nInit path length is : %.3lf\n", bestPath.len);
    printf("-----------------------------------\n\n");
}

void Print(Path t, int n)
{
    printf("Path is : ");
    for(int i = 0; i < n; i++)
    {
        if(i != n - 1)
            printf("%d-->", t.citys[i]);
        else
            printf("%d\n", t.citys[i]);
    }
    printf("\nThe path length is : %.3lf\n", t.len);
    printf("-----------------------------------\n\n");
}

Path GetNext(Path p, int n)
{
    Path ans = p;
    int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
    int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    while(x == y)
    {
        x = (int)(n * (rand() / (RAND_MAX + 1.0)));
        y = (int)(n * (rand() / (RAND_MAX + 1.0)));
    }
    swap(ans.citys[x], ans.citys[y]);
    ans.len = 0;
    for(int i = 0; i < n - 1; i++)
        ans.len += w[ans.citys[i]][ans.citys[i + 1]];
    cout << "nCase = " << nCase << endl;
    Print(ans, n);
    nCase++;
    return ans;
}

void SA(int n)
{
    double t = T;
    srand((unsigned)(time(NULL)));
    Path curPath = bestPath;
    Path newPath = bestPath;
    int P_L = 0;
    int P_F = 0;
    while(1)       //外迴圈，主要更新引數t，模擬退火過程
    {
        for(int i = 0; i < ILOOP; i++)    //內迴圈，尋找在一定溫度下的最優值
        {
            newPath = GetNext(curPath, n);
            double dE = newPath.len - curPath.len;
            if(dE < 0)   //如果找到更優值，直接更新
            {
                curPath = newPath;
                P_L = 0;
                P_F = 0;
            }
            else
            {
                double rd = rand() / (RAND_MAX + 1.0);
                //如果找到比當前更差的解，以一定概率接受該解，並且這個概率會越來越小
                if(exp(dE / t) > rd && exp(dE / t) < 1)
                    curPath = newPath;
                P_L++;
            }
            if(P_L > LIMIT)
            {
                P_F++;
                break;
            }
        }
        if(curPath.len < bestPath.len)
            bestPath = curPath;
        if(P_F > OLOOP || t < EPS)
            break;
        t *= DELTA;
    }
}

int main(int argc, const char * argv[]) {

    freopen("TSP.data", "r", stdin);
    int n;
    Input(p, n);
    GetDist(p, n);
    Init(n);
    SA(n);
    Print(bestPath, n);
    printf("Total test times is : %d\n", nCase);
    return 0;
}

注意由於是基於蒙特卡洛的方法，所以上面程式碼每次得出的結果並不完全一致。你可以通過增加迭代的次數來獲得一個更優的結果。

我們這裡需要說明的是，在之前的文章裡，我們用求最小值的例子來解釋模擬退火的執行：如果新一輪的計算結果更前一輪之結果更小，那麼我們就接受它，否則就以一個概率來拒絕或接受它，而這個拒絕的概率會隨著溫度的降低（也即是迭代次數的增加）而變大（也就是接受的概率會越來越小）。

但現在我們面對一個TSP問題，我們如何定義或者說如何獲取下一輪將要被考察的哈密爾頓路徑呢？在一元函式最小值的例子中，下一輪就是指向左或者向右移動一小段距離。而在TSP問題中，我們可以採用的方式其實是很多的。上面程式碼中GetNext()函式所採用的方式是隨機交換兩個城市在路徑中的順序。例如當前路徑為 A->B->C->D->A，那麼下一次路徑就可能是A->D->C->B->A，即交換B和D。而在文獻【3】中，作者取樣的程式碼如下（我們擷取一個片段，完整程式碼請參考原文）：

public class Tour{
    ... ...
    // Creates a random individual
    public void generateIndividual() {
        // Loop through all our destination cities and add them to our tour
        for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) {
          setCity(cityIndex, TourManager.getCity(cityIndex));
        }
        // Randomly reorder the tour
        Collections.shuffle(tour);
    }
    ... ...
}

    可見作者的方法是把上一輪路徑做一個隨機的重排（這顯然也是一種策略）。

TSP.data的資料格式如下，第一行的數字表示一個有多少座城市，第2至最後一行，每行有兩個數字表示，城市的座標（平面直角座標系）。例如：  6  20 80  16 84  23 66  62 90  11 9  35 28  最後請讀者自行編譯執行程式並觀察分析輸出的結果。