模拟退火算法python（含例程代碼與詳解)

模拟退火算法

• • 1.算法簡介
  • • 1.1 固體退火過程：
    • 1.2 Metropolis準則
    • 1.3 冷卻進度表
  • 2.算法流程
  • • 2.1 狀态産生函數（鄰域函數）
    • 2.2 狀态接受函數
    • 2.3 初溫
    • 2.4 溫度更新函數
    • 2.5 Metropolis抽樣穩定準則
    • 2.6 算法終止準則
  • 3.算法示例
  • • 3.1 TSP問題
    • 3.2 背包問題
    • 3.3 函數優化問題的求解
  • 4.算法實作
  • 5.進階

1.算法簡介

模拟退火算法（Simulated Annealing，SA）是一種模拟實體退火的過程而設計的随機優化算法，結合爬山法和随機行走算法，同時避免算法進入局部最優，早期用于組合優化，後來發展成一種通用的優化算法。它的基本思想最早在1953年就被Metropolis提出，但直到1983年Kirkpatrick等人才設計出真正意義上的模拟退火算法并進行應用。

該算法采用類似于實體退火的過程，先在一個高溫狀态下（相當于算法随機搜尋），然後逐漸退火，在每個溫度下（相當于算法的每一次狀态轉移），徐徐冷卻（相當于算法局部搜尋），最終達到實體基态（相當于算法找到最優解）。

高溫過程——增強粒子的熱運動，使其偏離平衡位置，目的是消除系統原先可能存在的非均勻态；

等溫過程——退火過程中要讓溫度慢慢降低，在每一個溫度下要達到熱平衡狀态，對于與環境換熱而溫度不變的封閉系統滿足自由能較少定律，系統狀态的自發變化總是朝自由能減少的方向進行，當自由能達到最小時，系統達到平衡态；

冷卻過程——使粒子熱運動減弱并漸趨有序，系統能量逐漸下降，進而得到低能的晶體結構。當液體凝固為固體的晶态時退火過程完成。

是以模拟退火算法從某一高溫出發，在高溫狀态下計算初始解，然後以預設的鄰域函數産生一個擾動量，進而得到新的狀态，即模拟粒子的無序運動，比較新舊狀态下的能量，即目标函數的解。如果新狀态的能量小于舊狀态，則狀态發生轉化；如果新狀态的能量大于舊狀态，則以Metropolis接受準則發生轉化。當狀态穩定後，便可以看作達到了目前狀态的最優解，便可以開始降溫，在下一個溫度繼續疊代，最終達到低溫的穩定狀态，便得到了模拟退火算法産生的結果

該算法的關鍵點如下：

1、對固體退火過程的模拟；

2、采用Metropolis接受準則；

3、用冷卻進度表控制算法程序，使算法在多項式時間裡給出一個近似解。

固體退火過程是SAA的實體背景；Metropolis接受準則使算法跳離局部最優 “險井”；而冷卻進度表的合理選擇是算法應用的前提。

下面依次介紹上述三個關鍵點

1.1 固體退火過程：

對于該部分如果不了解，可以先簡單跳過。

固體退火過程的數學表述：

在溫度T，分子停留在狀态r滿足Boltzmann機率分布

溫度低時能量低的微觀狀态機率大，溫度趨于零時，固體幾乎處于機率最大、能量最小的基态。

在同一個溫度T，標明兩個能量E1<E2，有

是以，在同一個溫度，分子停留在能量小的狀态的機率比停留在能量大的狀态的機率要大。

1.2 Metropolis準則

Metropolis準則（1953）——以機率接受新狀态，固體在恒定溫度下達到熱平衡的過程可以用Monte Carlo方法（計算機随機模拟方法）加以模拟。

Monte Carlo模拟退火過程：蒙特卡羅(Monte Carlo)方法，或稱計算機随機模拟方法，是一種基于“随機數”的計算方法。

若在溫度T，目前狀态i → 新狀态j

若Ej<Ei，則接受 j 為目前狀态；

否則，若機率 p= e x p [ − ( E j − E i ) / k B T ] exp[-(Ej-Ei)/k_BT] exp[−(Ej−Ei)/kB​T](也就是狀态j與狀态i機率的比值) 大于[0,1)區間的随機數，則仍接受狀态 j 為目前狀态；若不成立則保留狀态 i 為目前狀态。

可以看出，高溫下可接受與目前狀态能差較大的狀态為重要狀态，而在低溫下隻能接受與目前狀态能差較小的新狀态為重要狀态。在溫度趨于零時，就不再接受Ej>Ei的新狀态j了.如此反複，達到系統在此溫度下的熱平衡。這個過程稱作Metropolis抽樣過程。上面這種方法被稱為重要性采樣法，該方法以機率形式接受新狀态，避免局部最優。

1.3 冷卻進度表

在Metropolis抽樣過程中溫度T緩慢的降低，模拟退火過程就是通過T參數的變化使狀态收斂于最小能量處。因而，T參數的選擇對于算法最後的結果有很大影響。初始溫度和終止溫度設定的過低或過高都會延長搜尋時間。降溫步驟太快，往往會漏掉全局最優點，使算法收斂至局部最優點。降溫步驟太慢，則會大大延長搜尋全局最優點的計算時間，進而難以實際應用。是以選擇一個合适的降溫函數非常重要，冷卻進度表就是指的從某一高溫狀态T向低溫狀态冷卻時的降溫函數。

對實體退火過程和模拟退火算法進行一個簡單的對比：

2.算法流程

有了上面對該算法基礎的了解後，下面對算法的流程進行介紹：

可以看出該算法有内外兩次循環，内循環用來模拟同一溫度下多次的狀态轉移，也稱為抽樣穩定準則（為什麼稱為“抽樣穩定”準則，我的了解是等溫過程在鄰域選擇下一次狀态轉移，一般需要在鄰域内采樣多次，進行多次狀态轉移，進而達到這一溫度下的穩定狀态），外循環作為算法終止準則。

另外其中有3個比較重要的函數，分别為狀态産生函數（鄰域函數），狀态接受函數，退溫函數。同時另外還需注意到一個關鍵部分就是初溫、初始解的初始化。下面對上述内容展開介紹：

2.1 狀态産生函數（鄰域函數）

設計的出發點：盡可能保證産生的候選解遍布全部解空間。該函數包括兩部分的内容：産生候選解的方式、候選解産生的機率分布。前者決定由目前解産生候選解的方式，後者決定在目前解産生的候選解中選擇不同狀态的機率。

候選解的産生方式由問題的性質決定，通常在目前狀态的鄰域結構内以一定機率方式産生，而鄰域函數和機率方式可以多樣化設計，其中機率分布可以是均勻分布、正态分布、指數分布、柯西分布等。

2.2 狀态接受函數

狀态接受函數的目的是盡可能接受優化解，狀态接受函數一般以機率的方式給出，不同接受函數的差别主要在于接受機率的形式不同，設計狀态接受機率，應遵循的原則：

• 固定溫度下，接受使目标函數值下降的候選解的機率要大于使目标函數值上升的候選解的機率。
• 随溫度的下降，接受使目标函數值上升的解的機率要逐漸減小。
• 當溫度趨于零時，隻能接受目标函數值下降的解。

下面為标準形式下的狀态接受機率，通過将該機率與[0,1]之間的随機數進行比較，進而決定是否接受該候選解

2.3 初溫

實驗表明：初溫值隻要選擇充分大，獲得高品質解的機率就大！但花費計算時間增加。

是以，初溫的選擇要足夠高。初溫的确定應折衷考慮優化品質和優化效率，有以下兩種方式：

• 均勻抽樣一組狀态，選各狀态目标值的方差

• 利用經驗公式

  a.随機産生一組狀态，确定兩兩狀态間的最大目标值差，然後依據內插補點，利用一定的函數确定初溫。例如

  

  b.

  

2.4 溫度更新函數

溫度更新函數即溫度的下降方式，用于在外循環中修改溫度值，衰減量“以小為宜”。

實驗表明降溫速度越慢，獲得高品質解的幾率越大，但花費的計算時間将同時增加。是以可以使其溫度高時下降的慢些，溫度低時下降的快些。主要有以下幾種方式：

a.

Kirkpatrick首先提出對于lambd可選0.95，Johnson,Bonomi及Lutton采用[0.5,0.99].

b.Nahar及Skiscim等人把[0, t 0 t_0 t0​]劃分成K個小區間，溫度更新函數為

2.5 Metropolis抽樣穩定準則

用于決定在各溫度下産生候選解的數目( L k L_k Lk​).

具體應與問題規模成比例。實驗表明高溫時疊代次數越多越好，低溫時疊代次數可以适當減少。

在這裡引入下馬爾可夫鍊

馬氏鍊模型：

時齊算法的收斂性：

結論1 時齊模拟退火算法對應的有限狀态馬氏鍊存在平穩分布。

結論2 當溫度趨于0時，馬氏鍊以機率1收斂到最優狀态集，而收斂到非最優狀态的機率為0。

實作途徑：通過各溫度下各狀态序列無限長得以實作！

非時齊算法的收斂性：

對于非時齊算法

SA算法從某一個初始狀态開始後，每一步狀态轉移均是在目前狀态的鄰域中随機産生新狀态，然後以一定機率進行接受的。接受機率僅依賴于新狀态和目前狀态，并由溫度加以控制。是以，SAA對應一個馬氏鍊。

對于時齊算法，若固定每一溫度，算法均計算馬氏鍊的變化直至平穩分布，然後下降溫度。

對于非時齊算法，若無需各溫度下算法均達到平穩分布，但溫度需按照一定的速率下降。或稱非平穩馬氏鍊算法。

馬爾可夫鍊這部分我目前的了解還不是很深刻，後面有了認識後會加以補充。先記住下面的結論：

非時齊SAA：每個溫度下隻産生一個或少量候選解。

時齊算法——常用Metropolis抽樣穩定準則包括：

• 檢驗目标函數值的均值是否穩定
• 連續若幹步目标函數值的變化較小
• 按一定的步數抽樣

2.6 算法終止準則

用于決定算法何時結束。

理論上要求溫度終值趨于零，通常的做法包括：

• 設定終止溫度的閥值
• 設定外循環疊代次數(6-50)
• 算法搜尋到的最優值連續若幹步保持不變
• 檢驗系統熵是否穩定(暫未了解)

小結：

模拟退火算法基本要素和設定方法

3.算法示例

3.1 TSP問題

解空間：

目标函數：

初始解：

狀态産生函數：新解的産生

互換操作(SWAP)：随機交換兩個不同城市的位置。

逆序操作(INV)：兩個不同随機位置的城市逆序。

插入操作(INS)：随機選擇某個城市插入到不同随機位置。

初始溫度的計算：

for i=1:100
         route=randperm(CityNum);
         fval0(i)=CalDist(dislist,route);
     end
     t0=-(max(fval0)-min(fval0))/log(0.9);

           

3.2 背包問題

已知背包的裝載量為c=8，現有n=5個物品，它們的重量和價值分别是(2, 3, 5, 1, 4)和(2, 5, 8, 3, 6)。試使用模拟退火算法求解該背包問題，寫出關鍵的步驟。

求解：假設問題的一個可行解用0和1的序清單示，例如i=(1010)表示選擇第1和第3個物品，而不選擇第2和第4個物品。用模拟退火算法求解此例的關鍵過程如圖所示：

3.3 函數優化問題的求解

4.算法實作

下面以一個非常簡單的例子進行算法實作，優化f(x)= x 3 − 60 x 2 − 4 x + 6 x^3-60x^2-4x+6 x3−60x2−4x+6，求其在定義域[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)
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)  #降溫函數，也可使用T=0.9T

print(x, aimFunction(x))
           

運作結果：

最小值大概在40左右，通過求導計算得到最小值為40.033，該模拟退火算法最終的計算結果如下

為尋找在有限時間逼近全局最優的模拟退火算法，設定了許多控制算法收斂的參數。在退火過程中指定了有限的退火溫度值和在每一溫度下的轉移數目。Kirlpatrick等人在退火步驟中設定的參數如下：

（1）初始溫度值：初始溫度值T0要選的足夠高，保證模拟退火算法中所有可能的轉移都能被接受。

（2）溫度的下降：原先使用指數函數實作溫度的下降。但是這種方法使降溫幅度過小，進而延長搜尋時間。在實際中，通常使用下式：

此處λ是一小于卻接近于1的常數。λ通常的取值在0.8至0.99之間。

（3）終止溫度：如果在連續的若幹個溫度下沒有可接受的新狀态，系統當機或退火停止。

5.進階

1.算法特點

可以保證全局最優特别适合組合優化問題，可以随機選擇初始解，對問題本身沒有特别要求，不會因為問題執行個體的改變影響性能，簡單易行，通用性好。

模拟退火算法（SAA）在某一初溫下，伴随溫度參數的不斷下降，結合機率突跳特性在解空間中随機尋找目标函數的全局最優解，即在局部最優解能機率性地跳出并最終趨于全局最優。

2.模拟退火算法在求解規模較大的實際問題時，往往存在以下缺點：

（1）收斂速度比較慢。

（2）盡管理論上隻要計算時間足夠長，模拟退火法就可以保證以機率1收斂于全局最優點。但是在實際算法的實作過程中，由于計算速度和時間的限制，在優化效果和計算時間二者之間存在沖突，因而難以保證計算結果為全局最優點，優化效果不甚理想。

（3）在每一溫度下很難判定是否達到了平衡狀态。

為此，人們對模拟退火算法提出了各種各樣的改進，其中包括并行模拟退火算法、快速模拟退火算法（Cauchy機）和對模拟退火算法中各個函數和參數的重新設計等。

3.改進的可行方案

（1）設計合适的狀态産生函數；

（2）設計高效的退火曆程；

（3）避免狀态的迂回搜尋；

（4）采用并行搜尋結構；

（5）避免陷入局部極小，改進對溫度的控制方式；

（6）選擇合适的初始狀态；

（7）設計合适的算法終止準則。

改進的方式：增加某些新的環節

（1）增加升溫或重升溫過程，避免陷入局部極小；

（2）增加記憶功能（記憶“Best so far”狀态）；

（3）增加補充搜尋過程（以最優結果為初始解）；

（4）對每一目前狀态，采用多次搜尋政策，以機率接受區域内的最優狀态；

（5）結合其它搜尋機制的算法；

（6）上述各方法的綜合。

改進的思路

（1）記錄“Best so far”狀态，并即時更新；

（2）設定雙門檻值，使得在盡量保持最優性的前提下減少計算量，即在各溫度下目前狀态連續 m1 步保持不變則認為Metropolis抽樣穩定，若連續 m2 次退溫過程中所得最優解不變則認為算法收斂。

改進的退火過程

（1）給定初溫t0，随機産生初始狀态s，令初始最優解s*=s，目前狀态為s(0)=s，i=p=0；

（2）令t=ti，以t，s和s(i)調用改進的抽樣過程，傳回其所得最優解s’和目前狀态s’(k)，令目前狀态s(i)=s’(k)；

（3）判斷C(s*)<C(s*’)? 若是，則令p=p+1；否則，令s*=s*’，p=0；

（4）退溫ti+1=update(ti)，令i=i+1；

（5）判斷p>m2? 若是，則轉第(6)步；否則，傳回第(2)步；

（6）以最優解s*作為最終解輸出，停止算法。

改進的抽樣過程

（1）令k=0時的初始目前狀态為s’(0)=s(i)，q=0；

（2）由狀态s通過狀态産生函數産生新狀态s’，計算增量∆C’=C(s’)-C(s)；

（3）若∆C’<0，則接受s’作為目前解，并判斷C(s*’)>C(s’)? 若是，則令s*’=s’，q=0；否則，令q=q+1。若∆C’>0，則以機率exp(-∆C’/t)接受s’作為下一目前狀态；

（4）令k=k+1，判斷q>m1? 若是，則轉第(5)步；否則，傳回第(2)步；

（5）将目前最優解s*’和目前狀态s’(k)傳回改進退火過程。

時間不變的噪聲算法（TINA Time-invariant noise algorithm）

狀态産生函數中擾動強度不随時間改變，而是和能量大小相關，能量大的擾動大，能量小的擾動小，能量為零，擾動也為零，算法停止。

單調升溫(Monotonic temperature rising) SA

在算法退火後期，溫度很低且陷入局部極小解的時，算法很難跳出。是以，可以适當重新提高溫度，促使算法跳出。

記憶指導SA(Simulated Annealing with Memmory Guidance ,簡記為SAMG)

增加一個記憶裝置，存儲算法計算過程産生的最好的解，以這個解為最終解。

自适應SA算法

根據鄰域搜尋進展的回報資訊, 自适應确定溫度變化和鄰域搜尋強度

特點：

1. 退火過程中溫度參數變化符合幅值遞減的下降總趨勢, 但不排除局部升溫的可能, 以保證尋求到合适的溫度序列, 避免陷入局部最優;
2. 算法的終止條件依據退火溫度和鄰域搜尋進展狀态設計;
3. 每一溫度下算法的疊代次數随溫度下降而遞增, 鄰域搜尋強度依其對目标函數的貢獻動态配置設定;
4. 溫度變化、鄰域搜尋和終止條件的控制機制由算法過程自動觸發。

并行性

操作并行性：各個環節同時處理；

程序并行性：同時多個算法運作；

空間并行性：解空間分解分别處理，最終組合。

4.傳統優化算法與全局優化算法的對比：

全局優化算法：

1、不依賴于初始條件；

2、不與求解空間有緊密關系，對解域，無可微或連續的要求。求解穩健，但收斂速度慢。能獲得全局最優。适合于求解空間不知的情況

傳統優化算法：

1、依賴于初始條件。

2、與求解空間有緊密關系，促使較快地收斂到局部解，但同時對解域有限制，如可微或連續。利用這些限制，收斂快。

3、有些方法，如Davison-Fletcher-Powell直接依賴于至少一階導數；共轭梯度法隐含地依賴于梯度。

參考連結：

https://blog.csdn.net/WFRainn/article/details/80303138

https://www.cnblogs.com/ranjiewen/p/6084052.html

https://blog.csdn.net/To_be_to_thought/article/details/81914417