在讨论非完美算法时，模拟退火是经常被提及的算法。但多数论文对该算法的解释一般都限于与固体退火原理的类比：将固体加温至温度充分高，再让其徐徐冷却。加温时，固体内部粒子随温升变为无序状，内能增大；而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态,内能减为最小。

然而，固体退火原理对我们而言是非常陌生的。用一个陌生的原理来类比模拟退火算法，并不能帮助我们理解算法，反而带来许多更难于理解的“恐怖”术语：马尔可夫过程、Metropolis准则等等。尽管模拟退火算法的最初灵感来自固体退火过程，但能不能抛开固体退火原理，从另外的角度来理解模拟退火算法？

我们从一个简单的搜索算法出发，类比解释模拟退火算法的基本原理，并借助于实例来加深对算法的理解。

一、爬山法

爬山法是一种局部搜索算法，它运行迭代改进技术，在搜索空间中选取一个当前点。每次迭代从当前点的邻域中选取一个新点，如果这个新点的评估值优于当前点，则它就取代当前点；否则就选取邻域中的其它点来进行比较。算法一般在没有进一步改进的可能时或者运行了规定次数后终止。

显然，这种爬山法只能提供局部最优解，并且解的质量还依赖于起始点的选取。此外，也没有一个通用的过程来判断它与全局最优解之间的差距有多远，因为全局最优解是未知的。通常我们不得不从大量不同的起始点开始爬山法，并寄希望于这些起始点中总有一些会通向全局最优解。下面是算法伪代码：

void hill_climb()

{

   i = 0;             //爬山次数

   初始化best;     //全局最优解

   do {

              local = false;

              随机选取当前点Vc;

              do {

                      在Vc的邻域中选择所有新点

                      从新点集中选取评估函数eva l的最好值点Vn

                      if ( eva l(Vn) 优于 eva l(Vc))

                            Vc = Vn;

                      else

                            local = true;

              }while(local==false);

               i++;

               if ( eva l(Vc) 优于 eva l(best) )

                     best = Vc;

       }while(i < MAX);       //MAX是最大的爬山次数

}

上述算法的内层循环总是返回一个局部最优解Vc。它通过新的（随机）位置开始一次新的搜索（外层循环）来跳离局部最优。

下面对这一算法作一些修改：

（1）不是通过枚举当前点Vc邻域内的所有点来选出最好点，而是随机选取邻域中的任一点Vn;

（2）以一定概率接受这个新点，接受概率取决于这两个点的相对性能，即它们对应的评估函数值之间的差别。

通过这样的简单修改后，得到一个新算法，即随机爬山法：

void random_climb()

{

       i = 0;

       随机选择一个当前点Vc;

       do {

              从Vc的邻域中随机选择任一点Vn;

              计算p = 

              x = random[0,1];

              if (p > x)

                Vc = Vn;

              i++;

       }while(i < MAX);

}

这个算法的特点是：首先，只有一个循环，我们并不需要从不同的随机点开始重复这一迭代爬山过程；其次，所选择的新点Vn以一定概率p被接受，这意味着从当前点Vc转移到新点Vn是概率性的，新接受的点的性能有可能比当前点更糟糕。接受概率p依赖于两个点之间的质量差异，即eva l(Vc)-eva l(Vn)以及附加参数T的值。我们分析这两方面因素对p的影响。首先分析T的作用。

假设Vc和Vn的评估值分别是eva l(Vc)=108和eva l(Vn)=120，它们之间的差值等于-12，表明新点Vn比Vc要好。那么，对于不同的T值，接受这个新点的概率如何呢？下表列出了一些结果。

从表中数据可得出：T越大，Vn和Vc之间的相对质量的重要性越小。特别地，如果T很大，则接受的概率接近0.5，等同于纯随机搜索；如果T很小，则接受的概率接近1，表明Vn好于Vc时，一定被接受，随机爬山法就变成普通的爬山法了。因此，对于特定的问题，必须找到适当的T值，使之不要太大也不要太小。我们可以推知：T越大，越可能接受质量差的点，也就越有机会跳离局部搜索空间。但如果T一直保持很大的值，那么得到的近似最优解难以保证质量。

当T与Vc固定时，接受概率仅依赖于Vn的评估值。设T=12，eva l(Vc)=108。下表给出了不同的eva l(Vn)对接受概率的影响。

从表中数据可以分析得出：（1）如果Vn与Vc有相同的质量，则接受的概率是0.5，这是合理的。因为两个点质量相同，无所谓选择哪个点。（2）如果新点更差，则接受的概率小于0.5，并且其值随评估差值（正值）增大而减小。（3）如果新点更好，则接受的概率大于0.5，并且其值随评估差值（负值）增大而增大。如果新点比当前点好得多，则接受的概率接近1。

结论是：当T与Vc固定时，Vn的质量越好，接受的概率越大。

二、模拟退火算法

模拟退火算法是受上述概率公式规律的启发，做如下的操作：多次迭代搜索全局最优解，每次迭代的参数T都适当减小。开始时，T设为一个较大的值，此时接受质量差的解的概率大，使得有较多的机会跳离局部最优，这时内层循环类似于一个单纯的随机搜索。然后逐渐减小T值，使得接受质量差的解的概率越来越小，接受质量好的解的概率越来越大。到算法结束时，T值变得很小，使得模拟退火的最后阶段类似于普通的爬山法：如果新解比当前解更好，则总是接受它。下面是算法伪代码：

void sa()

{

   i = 0;

   初始化T;

   随机选择一当前点Vc;

   do {

              do {

                     在Vc邻域内选择一新点Vn;

                     if ( eva l(Vn) > eva l(Vc) )

                       Vc = Vn;

                     else {

                            p =

                            x = random[0, 1];

                            if (p > x)

                              Vc = Vn;

                     }

              }while(不满足终止条件);

              T按一定规则减小;

              i++;

       }while(不满足停机条件);

}

模拟退火算法在实现时，要确定四个特殊的条件：

（1）初始“温度”T的值是多少？初始温度T=1000，100还是其它值呢？

（2）温度应衰减多少或者多大因子？是每次减少5，还是%1？

（3）如何选取终止条件？是执行一定数目的迭代还是使用其它规则？

（4）如何选择停机条件，即“凝固”温度是多少？

三、实例

旅行商问题，即TSP问题（Travelling Salesman Problem）是著名的数学问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径长度为所有路径之中的最小值。

图1旅行商问题

一个最容易想到的方法是利用排列组合的方法把所有的路径都计算出来，并逐一比较，选出最小的路径。虽然该方法在理论上是可行的，但路径的个数等于n!。当城市个数较大时，该方法的求解时间是难以忍受的，甚至是不可能完成的。当包含20个城市时，20！= 2432902008176640000，以每秒1亿次的计算速度来估算，求解时间长达七百多年！因此，可以认为，旅行商问题的全局最优解是无法确定的，我们只可能得到近似最优解。这样，这类问题追求的目标就变成：以尽可能短的时间求得质量尽可能好的近似解。模拟退火算法能非常好地解决这一问题。

求解TSP的模拟退火算法模型可描述如下：

解空间：解空间S是遍访每个城市恰好一次的所有路经，解可以表示为{w1,w2 ,……, wn}，w1, ……, wn是1,2,……,n的一个排列，表明w1城市出发，依次经过w2, ……, wn城市，再返回w1城市。初始解可选为(1,……, n) ;

目标函数：目标函数为访问所有城市的路径总长度;

我们要求的最优路径为目标函数为最小值时对应的路径。  

新解的产生：产生新解的算法很多。不同的算法对解的质量有较显著的影响。一种简单的算法是：

随机产生1和n之间的两相异数k和m，不妨假设k<m，则将原路径

　　(w1,w2,…,wk,wk+1,…,wm,wm+1,…,wn)

　　变为新路径：

　　(w1,w2,…,wm,wk+1,…,wk,wm+1,…,wn)

上述变换方法就是将k和m对应的两个城市在路径序列中交换位置,称为2-opt映射。根据上述描述，模拟退火算法求解TSP问题的流程框图如下：

图2 TSP模拟退火流程图

参考源程序：

/*

模拟退火算法解决TSP问题

输入格式(tsp.in)：

第1行:1个整数N，表示城市的数量

第2..N+1行：每行有2个空格分开的整数x,y，第i+1行的x,y表示城市i的坐标

*/

#include <iostream>

#include <cmath>

using namespace std;

const int MAXN = 501;           //最大城市数500

const double INIT_T = 300;      //初始温度

const double RATE = 0.95;        //温度衰减率

const double FINAL_T = 1E-10;    //凝固温度

const int IN_LOOP = 13000;      //内层循环次数

const int OUT_LOOP = 100;       //外层循环次数

const int P_LIMIT = 10000;       //概率选择次数

struct path {                   //路径类型

    int City[MAXN];             //依次遍历的城市的序号

    double Length;              //所有城市的路径总长度

};

int N;                          //城市数量

double D[MAXN][MAXN];        //任意两个城市之间的距离

path bestpath;                   //最优的遍历路径

inline double dist(int x1, int y1, int x2, int y2) //计算两点距离

{

    return sqrt(double((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)));

}

inline double totaldist(path p)        //计算路径P的总长度

{

    int i;

    double ret = 0;

    for (i=1; i<N; i++)

      ret += D[p.City[i]][p.City[i+1]];

    return ret;

}

void init()                     //读入数据，并初始化

{

    int C[MAXN][2];                 //城市的坐标

    int i, j, k;

    freopen("tsp.in", "r", stdin);

    freopen("tsp.out", "w", stdout);

    scanf("%d", &N);

    for (i=1; i<=N; i++)

      scanf("%d%d", &C[i][0], &C[i][1]);

    for (i=1; i<N; i++)     //计算任两个城市之间的长度

      for (j=i+1; j<=N; j++)

      {

          D[i][j] = D[j][i] = dist(C[i][0], C[i][1], C[j][0], C[j][1]);

      }

    for (i=1; i<=N; i++)    //最优解的初始状态

      bestpath.City[i] = i;

    bestpath.Length = totaldist(bestpath);

    srand((unsigned)time(NULL));

}

path getnext(path p)   //新解产生函数

{

    int x, y;

    path ret;

    ret = p;

    do {

        x = rand() % N + 1;

        y = rand() % N + 1;

    }while(x == y);

    swap(ret.City[x], ret.City[y]);

    ret.Length = totaldist(ret);

    return ret;

}

void sa()

{

    double T;               //温度

    path newpath, curpath;    //当前路径和新路径

    int i, P_t=0, A_t=0;

    double delta;

    T = INIT_T;

    curpath = bestpath;

    while(true)

    {

        for (i=1; i<=IN_LOOP; i++)

        {

            newpath = getnext(curpath);

            delta = newpath.Length - curpath.Length;

            if (delta < 0.0)

            {

                curpath = newpath;

                P_t = 0;

                A_t = 0;

            }

            else

            {

                double rnd = rand()%10000 /10000.0;

                double p = exp(-delta/T);

                if (p > rnd)

                    curpath = newpath;

                P_t++;

            }

            if (P_t >=P_LIMIT)

            {

                A_t++;

                break;

            }

        }

        if (curpath.Length<bestpath.Length)

          bestpath = curpath;

        if ( A_t >= OUT_LOOP || T < FINAL_T) break;

        T = T * RATE;

    }

}

int main()

{

    init();

    printf("Initial Length: %.4f\n", bestpath.Length);

    sa();

    printf("Approx. Min Length: %.4f\n", bestpath.Length);

    return 0;

}

四、算例实测

上图左边是50个城市TSP问题的初始状态，右边是用模拟退火算法计算完成后的状态。计算时间大约是3.6秒。每次计算出的最短路径长度都有细微差别，但相较于穷举算法与模拟退火算法之间的巨大时间差异，这点差别是完全可以接受的。因此，模拟退火算法确实是一种优秀的近似计算算法。