夜深人静写算法（十二）- 模拟退火

一、引例

        1、函数最值

       函数最值分为函数最大值和函数最小值，最小值即定义域内函数的最小值， 最大值即定义域内函数的最大值。函数最大（小）值的几何意义为函数图像的最高（低）点的纵坐标。

       那么，让我们来看几种简单的情况：

       1) 一次函数

图一-1-1

       画出函数图像如下：

图一-1-2

        在定义域[x1,x2]内，函数的最小值和最大值分别取在两端点上。

        2) 二次函数

图一-1-3

图一-1-4

       由于该函数存在极值点，所以定义域的不同，函数最值也不尽相同。当然也可以通过函数图像分析出来，分三种情况讨论：

       a)当定义域[x1,x2]满足x2<=x0时，最大值取在x1，最小值取在x2；

       b)当定义域[x1,x2]满足x1>=x0时，最大值取在x2，最小值取在x1；

       c)当定义域[x1,x2]满足x1<x0<x2时，最大值为两端点的小者，最小值取在x0；

      3) 周期函数

图一-1-5

图一1-6

       这个函数为标准的正弦函数，周期为2π，所以计算最值的时候可以将定义域平移到[0,2π)区间，然后分情况讨论即可（参考二次函数的情况）。

       接下来，让我们看下更加复杂的情况下的函数最值。

         4) 复合函数

图一-1-7

       这种函数乍看之下，很难想象它的函数图像，那么我们可以利用描点法绘制出它的函数图像，点越多越精确。得到了一个令人绝望的函数图像。

图一-1-8

       这种情况下，如何求它在给定的定义域下的函数最值呢？

       这就是本文今天要介绍的算法---模拟退火。 模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能，目前已在工程中得到了广泛应用，诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。

二、袋鼠跳

        1、从天而降的袋鼠

       在介绍模拟退火之前，我们先把上面提到的问题用一种有趣的方式描述出来，这样有助于理解模拟退火的算法核心。

       我们将上文的函数定义为y=f(x)，图像想象成一座山，上面有无数个连绵起伏的山峰和山谷，现在我们要想办法找到这座山的最低谷，如图二-1-1所示。

图二-1-1

       首先，我们在x坐标上随机采点，从天而降一些袋鼠。每只袋鼠有自己的思维，它们都可以各自选择是往左跳还是往右跳。

图二-1-2

       宏观来说，它们的行为是并行的（实际上，只需要用一个枚举来模拟每只袋鼠的行为，并非真正的CPU并行计算），所以我们其实只需要考虑某一只袋鼠的行为。

        2、选择性跳跃

       对于每只袋鼠，我们首先定义一个步长S，如果当前袋鼠的位置为x，那么在x-S和x+S两个位置中挑选一个相对较低的位置（即min{f(x-S),f(x+S)}），比较这个位置和当前袋鼠所在的位置。如果比当前点更低，那么毫不犹豫的跳过去，否则以一定概率选择性的跳跃。

图二-1-3

        3、退火（减少步长）

       然后，等比例减少步长S（一般是乘上一个系数，比如0.99），继续迭代这样的操作，直到步长S小于某个精度，那时候袋鼠也跳不动了。

       这样做，只要参数控制的好，袋鼠有很大几率落在最低点，但是也不一定。所以我们选了很多袋鼠，然后并行迭代后取最后落在最低点的那只即可。好了，如果能够听懂，说明你已经无形之中领悟了模拟退火的精髓。

       三、模拟退火

        1、初衷

        模拟退火算法来源于固体退火原理，是基于蒙特卡罗迭代求解策略的一种随机概率寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

        模拟退火算法从某一较高初温出发，随着温度参数的不断下降，结合概率的跳跃特性，在解空间中随机寻找目标函数的全局最优解。即在局部最优解时能够概率性的跳出，并且最终趋于全局最优。

        根据Metropolis准则，粒子在温度T时趋于平衡的概率为：

（E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数）

       模拟退火算法更加直观的一种解释：想象成一个坑坑洼洼的土地，从天上丢下来一些弹性非常好的球，由于能量损失，这些球每次的弹跳高度都会下降，直到一直弹不起来。那么一开始它总能弹到一些本来不属于初始下降区域的地方，即跳出局部最优解。慢慢逼近全局最优解。

        2、算法过程

        (0) 初始化：初始温度T(充分大)，初始解状态X(是算法迭代的起点)，每个温度T下的迭代次数L；

        (1) for i = 1 to L 执行(2)(3)(4)(5)

        (2) 产生新解X′（以某种策略生成新解）；

        (3) 计算温度增量ΔT=h(X′)-h(X)，其中h(X)为估价函数；

        (4) 若ΔT<0则接受X′作为新的当前解，否则以概率e(-ΔT/T)接受X′作为新的当前解；

        (5) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法，否则当温度到达某个最小精度时结束。

        (6) 否则，逐渐减少T，然后转(1)。

        3、算法剖析

       这里的温度T就是袋鼠的跳跃步长，会随着时间的推移逐渐减少（等比下降，即退火过程）；初始解X是随机生成的，可以认为是上述问题中函数图像的x坐标，L就是袋鼠的个数。步骤(0)即初始化；步骤(1)则是并行计算每只袋鼠的跳跃情况；步骤(2)产生左右两个解，并且取更优解作为候选解X'。

        步骤(3)这一步出现了一个新的名词，估价函数。估价函数即对当前解的一个评估，从而决定是否要从当前解转向候选解。很显然，在求函数最小值的问题中，估价函数正好是函数值本身，即h(X)=f(X)。求最大值时，估价函数是函数值的相反数，即h(X)=-f(X)（原因是在该问题的前提下，估价函数一定是越小越优）。

        步骤(4)对应了袋鼠的选择性跳跃中的“选择”二字。如果能够确认更优，则直接替换当前解，否则按照上述概率进行筛选。逐渐迭代温度T，直接温度满足某个精度为止，这个很好理解，因为温度对应步长，当温度很小时，左右移动意义也不大了。

        4、算法实现

       1）类定义

[cpp] view plain copy

1. class simulatedAnnealing {  
2.     static const double minTemperature;      // 稳态（最低）温度  
3.     static const double deltaTemperature;    // 温度下降率  
4.     static const int solutionCount;          // 并行解个数  
5.     static const int candidateCount;         // 每个解的迭代次数  
6. private:  
7.     Point3D bound;  
8.     Point3D x[MAXC];  
9.     Point3DSet pointSet;  
10.     double temperature;  
11.   
12.     bool valid(const Point3D& pt);  
13.     double randIn01();  
14.     Point3D getRandomPoint();  
15.     Vector3D getRandomDirection();  
16.     Point3D getNext(const Point3D& now);  
17. public:  
18.     void start(double temper, Point3D ptBound, Point3DSet& ptSet);  
19.     double evaluateFunc(const Point3D& pt);  
20.     Point3D getSolution();  
21.     static simulatedAnnealing& Instance();  
22. };  
23. const double simulatedAnnealing::minTemperature = 1e-4;  
24. const double simulatedAnnealing::deltaTemperature = 0.95;  
25. const int simulatedAnnealing::solutionCount = 10;  
26. const int simulatedAnnealing::candidateCount = 30;  

       2）退火主流程

      三个参数temper代表退火的初始温度、ptBound代表最优解的范围、ptSet代表数据集。

      temper：初始化温度，即初温temper，直接赋值给类变量temperature即可。

      ptBound：从ptBound范围内随机挑选solutionCount个点作为初始点集，ptBound表示的是一个点，所以实际的初始点的范围就是从原点到ptBound随机取值，然而对于最优解在负数坐标的情况，需要将整个坐标系平移，即平移到所有解都位于原点之上。这里的点是个抽象的概念，可以是一维的（上文提到的x-y函数），也可以是二维的（平面坐标系上的点），也可以是三维的（空间坐标系中的点），当然也可以是更高维度的 。我写的是三维的模板，如果要转换成二维只要将第三维坐标z置0即可。

ptSet：辅助数据，用于求特定的问题，例如问题描述为求一个任意多边形的最大内接圆，那么这里任意的多边形就用这个ptSet点集来表示。

      退火主流程分成三步走：

       a）初始化参数和初始解集生成（代码中的步骤0和1）；

初始温度选择很重要，选的越大迭代时间越久越精确，会以时间作为代价。所以需要权衡两者之间的平衡。

       b）对每个初始解进行最优化筛选和替换（代码中的步骤2、3、4、5、6、7）

       minTemperature为最低温度控制，即当temperature小于等于这个温度时算法终止，一般问题如果求解精度在0.01，那么最低温度最好控制在0.001，即尽量比求解精度小1个数量级，这是为了提高准确性 。

       solutionCount作为初始解集数量，为事先定义的常量，在情况比较特殊（比如单调性很明显的单峰单谷问题）时可以取1，取得越大效率越慢，但是更加精确；相反，效率高但是精度有可能会下降。

       candidateCount为当前某个解的下一候选解的数量，即上文提到的袋鼠选择往左还是往右，对于一维的情况，即函数求极值，只需要将candidateCount设为2即可（因为只有左右两个方向），对于更高维的情况， 则需要随机方向，所以取值可能是30、50、100，具体问题具体分析，通常是实验结果。

       getNext(x[i])返回的是当前解x[i]的下一个候选解，并且用evaluateFunc计算估价函数后，取估价值最小的作为实际候选解。这个候选解需要迭代candidateCount次，选取最优，然后和当前解x[i]进行估价函数的比较。如果候选解的估价值更小，则直接替换当前解；否则，以一定概率接收候选解。注意：这里的一定概率也可能是零概率，也就是第二个分支永远跑不到，比如求单谷函数的最小值，完全没有必要接收估价更大的解，这种情况下，模拟退火算法退化成“梯度下降”算法。

       c）退火（代码中的步骤8）

       temperature 直接乘上一个常量deltaTemperature，一般取0.9左右，根据具体情况而定，越大精度越高，时间越慢；反之，则精度越低，时间越快。

来看下算法实现：

[cpp] view plain copy

1. void simulatedAnnealing::start(double temper, Point3D ptBound, Point3DSet& ptSet) {  
2.     // 0.初始化温度  
3.     temperature = temper;  
4.     bound = ptBound;  
5.     pointSet = ptSet;  
6.     int i, j;  
7.   
8.     // 1.随机生成solutionCount个初始解  
9.     for(i = 0; i < solutionCount; ++i) {  
10.         x[i] = getRandomPoint();  
11.     }  
12.   
13.     while (temperature > minTemperature) {  
14.         // 2.对每个当前解进行最优化选择  
15.         for(i = 0; i < solutionCount; ++i) {  
16.             double nextEval = INF;  
17.             Point3D nextOpt;  
18.             // 3.对于每个当前解，随机选取附近的candidateCount个点，并且将最优的那个解保留  
19.             for(j = 0; j < candidateCount; ++j) {  
20.                 Point3D next = getNext(x[i]);  
21.                 if(!valid(next)) {  
22.                     continue;  
23.                 }  
24.                 double Eval = evaluateFunc(next);  
25.                 if(Eval < nextEval) {  
26.                     nextEval = Eval;  
27.                     nextOpt = next;  
28.                 }  
29.             }  
30.   
31.             // 4.没有生成可行解  
32.             if(nextEval >= INF)  
33.                 continue;  
34.   
35.             // 5.计算生成的最优解和原来的解进行比较  
36.             double deltaEval = evaluateFunc(nextOpt) - evaluateFunc(x[i]);  
37.             if(deltaEval < 0) {  
38.                 // 6.比原来的解更优，直接替换  
39.                 x[i] = nextOpt;  
40.             }else {  
41.                 // 7.没有原来的解优，则以一定概率进行接收  
42.                 // 这个概率上限会越来越小，直到最后趋近于0  
43.                 // 理论上，这个分支也可能不考虑  
44.                 if( randIn01() < exp(-deltaEval/temperature) ) {  
45.                     x[i] = nextOpt;  
46.                 }      
47.             }  
48.         }  
49.         // 8.退火  
50.         temperature *= deltaTemperature;  
51.     }  
52. }  

       3）估计函数

模拟退火中的估价函数，类似A*，即通过某个解估算出一个价值（整数或者浮点数皆可），我们定义为估价函数越小，解越优。不同问题，估价函数的实现也不尽相同，这个在下一章节中会详细介绍。

       4）其它的一些小函数就不介绍了，完整C++代码，请参见 ： 模拟退火算法C++实现

        5、算法复杂度

       最后，分析一下实际退火过程的时间复杂度。首先最外层迭代，从最高温度到最低温度经过一个等比式的迭代，令初温为T，最低温度TMin，温度下降率DT，下降次数为N，那么有：

       从而得出N的最大值为：

       然后来看每次迭代，解集数量定义为L，每个解的候选解数量为C，每次计算估价函数的时间复杂度为H，则得到算法的时间复杂度为：

      得到算法复杂度后，我们就可以根据实际情况，设定L和C的值，权衡时间和精度进行算法求解了。

      四、模拟退火算法举例

      最后，我们对这个算法展开一些问题的探索，加深对算法本身的理解。

      1、函数最值

      以函数最小值为例，回到一开始提到的那个复合函数：

      【例题1】给定x1,x2，求该函数在区间[x1, x2]上的的最小值，精确到小数点后2位。

       直接套用模拟退火，初始温度可以取函数定义域的区间长度x2-x1，最小温度0.001，温度下降率0.95，初始解集10个，候选解集2个，估价函数直接取函数值，实现如下：

[cpp] view plain copy

1. double simulatedAnnealing::evaluateFunc(double x) {  
2.     double sx = sin(x);        // sin(x)计算比较费时，预先计算出来  
3.     return sx*sx - 2*sx + x*sx;  
4. }  

        2、曲面最值

      【例题2】给定一个二次曲面方程，保证这个曲面有一个最低点，求这个最低点的坐标。

图4-2-1

       这个问题其实就是求一个坐标，这个坐标在曲面上，且坐标的纵坐标分量最小（即y坐标）。因为有曲面方程，所以如果x和z的值确定，那么求解一个一元二次方程就能算出y，所以起初我们可以任意随机一个曲面上的点，然后通过模拟退火不断将点逼近到最低点。估价函数的实现比较简单：

[cpp] view plain copy

1. double simulatedAnnealing::evaluateFunc(const Point3D& pt) {  
2.     return pt.y;  
3. }  

       3、最 远 的最 近 点

      【例题3】给定一个矩形区域W,H(W,H<=10000) 和 N(N<=1000)个危险点，求矩形内一个点，这个点离最近的那个危险点最远（要求精确到小数点后1位）。如图红色点代表危险点，蓝色点代表离最近的危险点最远的那个最优解所在的位置。

图四-3-1

      初始化温度T=max(W,H),随机初始化一个矩形内的点作为初始解，将初始解作为当前解进行迭代计算。

       每次以当前解(图中的黑点)为圆心,温度T为半径生成一些方向随机的候选点集(图中天蓝和深蓝的点)，然后选择一个离所有危险点(图中红色点)最小距离最大的点(图中深蓝色点)作为当前解的候选替代者（暂时不进行替换）。

图四-3-2

       模拟退火算法中，估价函数的值越小越优，而这题求的是最大距离，也就是离最近的那个危险点距离更大的点更优，所以估价函数需要在求出最小距离后加上一个相反数。

       将候选最优解和当前解的估价函数作差，如果估价差值△E<0，说明候选最优解比当前解更优，则用候选最优解替换当前解；否则，随机一个0-1的随机数p，如果p<e(-△E/T)，则用候选最优解替换当前解，否则不进行替换。

       每次计算完毕，则将温度下降一定比率，直到温度退化到一定精度，则算法结束。

      这个算法中，温度T为每次选择候选点集的偏移半径；每个候选点到最近危险点的距离的相反数作为估价函数；估价函数差值作为退火时的能量差值作为新解取舍的判断依据。

      当然，初始化生成的位置可能不够好，导致最后陷入局部最优解，所以，我们可以随机生成一个初始解集合进行“并行”（并非真正的并行）计算。

       参数设定：最小温度0.001，温度下降率0.95，初始解集5个，候选解集30个，估价函数实现如下：

[cpp] view plain copy

1. double simulatedAnnealing::evaluateFunc(const Point3D& pt) {  
2.     double minDist = INF;  
3.     for(int i = 0; i < pointSet.n; ++i) {  
4.         double dist = (pointSet.p[i] - pt).len();  
5.         if(dist < minDist)   
6.             minDist = dist;  
7.     }  
8.     return - minDist;  
9. }  

      其中pointSet就是上文我们提到的辅助数据，用来记录题目中那些“危险点”。

       4、最近的最远点

      【例题4】 给定一个矩形区域W,H(W,H<=10000) 和 N(N<=1000)个点，求矩形内一个点，这个点离最远的那个点最近（要求精确到小数点后1位）。

这个问题和例题3，正好相反，同样可以利用模拟退火的性质求解，唯一不同的就是估价函数不需要取相反数，并且选择点的时候是选取距离最大的那个点。

       5、任意多边形的最大内接圆

      【例题5】给定一个任意简单多边形（可能是凹多边形）以及一个半径为r的圆，问这个圆能否放入这个多边形内。

图四-5-1

      作为一个判定性问题，同样可以用模拟退火来进行判定。首先，初始化一些初始点集，对于每个点，以它到所有多边形线段的距离的最小值的相反数作为估价函数。因为我们想要的是让这个点到多边形所有线段的最小距离最大，所以估价函数需要取相反数。仔细一想，这个问题和“最远的最近点”是一样的，唯一不同的是它有跳出条件，即当这个最小距离已经大于等于r时，可以跳出迭代。

      需要注意的是，圆的圆心需要在多边形内，否则可能出现以图四-5-2的非法情况。

图四-5-2

      五、模拟退火相关题集整理

PKU 1379 Run Away

寻找最远的最近点

PKU 2069 Super Star

寻找最近的最远点（3维）

PKU 2420 A Star not a Tree?

寻找到所有点距离最小的点

PKU 2600 Gemetrical dreams

模拟向量旋转

HDU 3202 Circle and Triangle

三角形和圆的最大相交面积（难题）

HDU 3644 A Chocolate Manufacturer's Problem

任意简单多边形的最大内接圆

HDU 3932 Groundhog Build Home

寻找最近的最远点

HDU 4717 The Moving Points

寻找距离最小的最远点对

HDU 5017 Ellipsoid

曲面最小距离