模擬退火演算法(SimulatedAnnealing)是基於Monte-Carlo迭代求解策略的一種隨機尋優演算法,主要用於組合優化問題的求解。

假設現在有這麼一個函式：

現要求其在[0,100]範圍內的最小值，如果不求導計算，可能第一反應都是窮舉法，把範圍內每個值都算一遍再比較大小。如果求的是整數範圍，那麼要算100遍，但是如果要精確到小數後8位，則要算10000000000次，即便使用計算機依然是一個龐大的運算過程。而優化問題中很多都類似於問題，無法用窮舉法解出答案，我們叫這類問題為NP難問題（可檢視維基百科：NP-hard），於是，有人提出了爬山法：

但是這個方法的缺點在於最優解的產生依賴於最初值的選取，無法解決非凸函式，即容易收斂於區域性最優解：

同時，也無法解決有平臺的函式：

於是，Kirkpatrick等提出了模擬退火演算法，它是一種啟發式搜尋演算法，即按照預定的控制策略進行搜尋，在搜尋過程中獲取的中間資訊將用來改進控制策略

1.模擬退火演算法的原理

1.1概念

模擬退火演算法的思想借鑑於固體的退火原理，當固體的溫度很高的時候，內能比較大，固體的內部粒子處於快速無序運動，當溫度慢慢降低的過程中，固體的內能減小，粒子的慢慢趨於有序，最終，當固體處於常溫時，內能達到最小，此時，粒子最為穩定。模擬退火演算法便是基於這樣的原理設計而成。

模擬退火演算法從某一高溫出發，在高溫狀態下計算初始解，然後以預設的鄰域函式產生一個擾動量，從而得到新的狀態，即模擬粒子的無序運動，比較新舊狀態下的能量，即目標函式的解。如果新狀態的能量小於舊狀態，則狀態發生轉化；如果新狀態的能量大於舊狀態，則以一定的概率準則發生轉化。當狀態穩定後，便可以看作達到了當前狀態的最優解，便可以開始降溫，在下一個溫度繼續迭代，最終達到低溫的穩定狀態，便得到了模擬退火演算法產生的結果。

1.2狀態空間與鄰域函式

狀態空間也稱為搜尋空間，它由經過編碼的可行解的集合所組成。而鄰域函式應儘可能滿足產生的候選解遍佈全部狀態空間。其通常由產生候選解的方式和候選解產生的概率分佈組成。候選解一般按照某一概率密度函式對解空間進行隨機取樣獲得，而概率分佈可以為均勻分佈、正態分佈、指數分佈等。

1.3狀態轉移概率（Metropolis準則）

狀態轉移概率是指從一個狀態轉換成另一個狀態的概率，模擬退火演算法中一般採用Metropolis準則，具體如下：

其與當前溫度引數T有關，隨溫度的下降而減小。

1.4冷卻進度表

冷卻進度表是指從某一高溫狀態T向低溫狀態冷卻時的降溫函式,設時刻的溫度為T(t)，則經典模擬退火演算法的降溫方式為：

而快速模擬退火演算法的降溫方式為：

另外還有幾種常用的降溫函式個人感覺大同小異，就不過多介紹了。

1.5初始溫度

一般來說，初始溫度越大，獲得高質量解的機率越大，但是花費的時間也會隨之增加，因此，初溫的確定應該同時考慮計算效率與優化質量，常用的方法包括：

(1)均勻抽樣一組狀態，以各狀態目標值的方差為初溫。

(2)隨機產生一組狀態，確定亮亮狀態間的最大目標值差，然後根據差值，利用一定的函式確定初溫，如：

其中Pr為初始接受概率。

(3)根據經驗公式給出

1.6迴圈終止準則

內迴圈終止準則：

（1）檢驗目標函式的均值是否穩定

（2）連續若干步的目標值變化較小

（3）按一定的步數進行抽樣

外迴圈終止準則

（1）設定終止溫度

（2）設定外迴圈迭代次數

（3）演算法搜尋到的最優值連續若干步保持不變

（4）檢驗系統熵是否穩定

Python實現過程：

下面便通過python求解開頭提到的問題，首先定義函式，然後通過pyplot看看函式在[0,100]上的大致影象：

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import math

#define aim function
def aimFunction(x):
    y=x**3-60*x**2-4*x+6
    return y
x=[i/10 for i in range(1000)]
y=[0 for i in range(1000)]
for i in range(1000):
    y[i]=aimFunction(x[i])

plt.plot(x,y)

可以看到最小值大概在40左右，通過求導計算得到最小值為40.033。

接下來便構造SA模型：

定義初溫、低溫閾值並通過隨機得到初始x，同時定義時刻t。通過均勻分佈構造鄰域函式，同時設定內迴圈次數為50次，降溫函式使用

程式碼實現如下：

T=1000 #initiate temperature
Tmin=10 #minimum value of terperature
x=np.random.uniform(low=0,high=100)#initiate x
k=50 #times of internal circulation 
y=0#initiate result
t=0#time
while T>=Tmin:
    for i in range(k):
        #calculate y
        y=aimFunction(x)
        #generate a new x in the neighboorhood of x by transform function
        xNew=x+np.random.uniform(low=-0.055,high=0.055)*T
        if (0<=xNew and xNew<=100):
            yNew=aimFunction(xNew)
            if yNew-y<0:
                x=xNew
            else:
                #metropolis principle
                p=math.exp(-(yNew-y)/T)
                r=np.random.uniform(low=0,high=1)
                if r<p:
                    x=xNew
    t+=1
    print t
    T=1000/(1+t)
    
print x,aimFunction(x)

經過迴圈輸出x與y，結果如下：

多次執行：

可以看到SA演算法很好的逼近了最優解。

Reference：