模拟退火算法

在实际日常中，人们会经常遇到如下问题：在某个给定的定义域X 内，求函数f(x) 对应的最优值。此处以最小值问题举例（最大值问题可以等价转化成最小值问题），形式化为：

min x∈X f(x). 

如果X 是离散有限取值，那么可以通过穷取法获得问题的最优解；如果X 连续，但f(x) 是凸的，那可以通过梯度下降等方法获得最优解；如果X 连续且f(x) 非凸，虽说根据已有的近似求解法能够找到问题解，可解是否是最优的还有待考量，很多时候若初始值选择的不好，非常容易陷入局部最优值。

随着日常业务场景的复杂化，第三种问题经常遇见。如何有效地避免局部最优的困扰？模拟退火算法应运而生。其实模拟退火也算是启发式算法的一种，具体学习的是冶金学中金属加热-冷却的过程。由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的，V.Čern在1985年也独立发明此演算法。

不过模拟退火算法到底是如何模拟金属退火的原理？主要是将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。若概率大于给定的阈值，则跳转到“邻居”；若概率较小，则停留在原位置不动。

一、模拟退火算法基本思想

模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。在迭代更新可行解时，以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以下图为例，假定初始解为左边蓝色点A，模拟退火算法会快速搜索到局部最优解B，但在搜索到局部最优解后，不是就此结束，而是会以一定的概率接受到左边的移动。也许经过几次这样的不是局部最优的移动后会到达全局最优点D，于是就跳出了局部最小值。

根据热力学的原理，在温度为T 时，出现能量差为dE 的降温的概率为p(dE) ，表示为:

p(dE)=exp(dEkT ). 

其中k 是波尔兹曼常数，值为k=1.3806488(13)×10 −23  ，exp 表示自然指数，且dE<0 。因此dE/kT<0 ，所以p(dE) 函数的取值范围是(0,1)。满足概率密度函数的定义。其实这条公式更直观意思就是：温度越高，出现一次能量差为p(dE) 的降温的概率就越大；温度越低，则出现降温的概率就越小。

在实际问题中，这里的“一定的概率”的计算参考了金属冶炼的退火过程。假定当前可行解为x ，迭代更新后的解为x_new ，那么对应的“能量差”定义为：

Δf=f(x_new)−f(x). 

其对应的“一定概率”为：

p(Δf)=⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ exp(−ΔfkT )exp(ΔfkT ) 最小值优化问题最大值优化问题  

注：在实际问题中，可以设定k=1 。因为kT 可以等价于一个参数 T 。如设定k=2 、T=1000 ，等于直接设定T=2000 的效果。

二、模拟退火算法描述

1. 初始化：初始温度T （充分大），温度下限T min  （充分小），初始解状态x (是算法迭代的起点)，每个T 值的迭代次数L ；
2. 对l=1,2,...,L 做第3至第6步；
3. 产生新解x_new : (x_new=x+Δx )；
4. 利计算增量Δf=f(x_new)−f(x) ，其中f(x) 为优化目标；
5. 若Δf<0 (若寻找最小值，Δf>0 )则接受x_new 作为新的当前解，否则以概率exp(−Δf/(kT)) 接受x_new 作为新的当前解；
6. 如果满足终止条件则输出当前解作为最优解，结束程序。(终止条件通常取为连续若干个新解都没有被接受时终止算法。)；
7. T 逐渐减少，且T>T min  ，然后转第2步。

三、模拟退火算法的优缺点

模拟退火算法的应用很广泛，可以高效地求解NP完全问题，如货郎担问题(Travelling Salesman Problem，简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等，但其参数难以控制，不能保证一次就收敛到最优值，一般需要多次尝试才能获得（大部分情况下还是会陷入局部最优值）。观察模拟退火算法的过程，发现其主要存在如下三个参数问题：

　　(1) 温度T的初始值设置问题
　　温度T 的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。

　　(2) 退火速度问题，即每个T 值的迭代次数
　　模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索是相当必要的，但这也需要计算时间。循环次数增加必定带来计算开销的增大。

　　(3) 温度管理问题
　　温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：

T=α×T.α∈(0,1). 

注：为了保证较大的搜索空间，α 一般取接近于1的值，如0.95、0.9。

四、模拟退火算法Python实战

经过上面理论知识的熏陶，相信大家已经对模拟退火算法有了较深入的理解，接下来通过实战再强化一下大家的认识，此处利用模拟退火算法求解如下优化问题：

minf(x)=(x−2)∗(x+3)∗(x+8)∗(x−9) 

# -*- coding: utf-8 -*-import numpy as npimport matplotlib.pyplot as pltdef inputfun(x):    return (x-2)*(x+3)*(x+8)*(x-9)initT = 1000 #初始温度minT = 1 #温度下限iterL = 1000 #每个T值的迭代次数delta = 0.95 #温度衰减系数k = 1initx = 10*(2*np.random.rand()-1)nowt = initTprint "初始解：",initxxx = np.linspace(-10,10,300)yy = inputfun(xx)plt.figure()plt.plot(xx,yy)plt.plot(initx,inputfun(initx),'o')#模拟退火算法寻找最小值过程while nowt>minT:    for i in np.arange(1,iterL,1):        funVal = inputfun(initx)        xnew = initx+(2*np.random.rand()-1)        if xnew>=-10 and xnew<=10:            funnew = inputfun(xnew)            res = funnew-funVal            if res<0:                initx = xnew            else:                p = np.exp(-(res)/(k*nowt))                if np.random.rand()<p:                    initx = xnew#            print initx-xnew#    print initx#    print nowt    nowt = nowt*deltaprint "最优解：",initxprint "最优值：",inputfun(initx)plt.plot(initx,inputfun(initx),'*r')plt.show()

可以看到，即使初始解选取的有风险，模拟退火算法经过迭代，也可以成功跳出局部最优值，最终收敛到问题的全局最优。

参考资料

1. http://wiki.mbalib.com/wiki/模拟退火算法 模拟退火算法