Python — 模擬退火演算法(Simulated Annealing, SA)求解旅行推銷員問題(Traveling Salesman Problem, TSP)
模擬退火演算法(Simulated Annealing, SA)是S.Kirkpatrick, C. D. Gelatt和M.P.Vecchi等人於1983年所提出的通用概率演算法,用來在一個範圍空間內搜尋問題的最佳解。是基於Monte-Carlo迭代求解策略的一種隨機尋優算法。
何謂模擬退火演算法(Simulated Annealing, SA)?
相較於其他演算法來說,模擬退火演算法屬於入門級的演算法,容易理解及應用,普遍用以求解最佳化問題與旅行推銷員問題。
其靈感來自於冶煉金屬的降溫過程,於此過程中熱運動趨近於穩定,等同於優化過程中解的收斂。在以下本文中,針對SA的算法原理除了會說明其流程步驟之外,也會說明如何使用python實現演算法。
何謂旅行推銷員問題(Traveling Salesman Problem, TSP)?
此部分主要來自搜尋維基百科,在此只會針對問題定義上大致做一個講解,若想詳細了解可以點擊以下連結,但相信大部分會點擊此篇文章的讀者,應都已具備相對應的觀念。
TSP定義為,給定一系列城市和每對城市之間的距離,求解存取每一座城市一次並回到起始城市的最短路徑。它是組合最佳化中的一個NP困難問題,在作業研究和理論電腦科學中非常重要。
在本文你將會學會以下技巧:
1.培養分析問題的能力。
2.以Python實現簡單的演算法。
3.了解旅行推銷員問題。
4.使用模擬退火演算法解決旅行推銷員問題。
SA演算原理及流程
首先初始化基本參數,決定起始溫度 t0、終止溫度 tmin、溫度衰退係數、初始解,並透過城市點 N 生成相對應的距離矩陣(Distance Matrix),於每次降溫過程中,重複 k 次運算以決定該溫度內的最適解。
在不改變起始和結尾城市的情況下,每次以交換 2 個城市順序來產生一新路徑,並計算新路徑的總距離,再和原距離相減以得到兩者之間的差異量diff。
因TSP本身為最小化問題,如果diff<0則直接替換新解,反之若是>0則以exp(-diff/t)的概率接受其成為最適解;在此 t 是當前溫度,原理為通過在溫度高時較有可能接受差解和溫度低時幾乎不接受差解從而避免陷入局部最優的情況。每個溫度內需完成 k 次運算,溫度從t0降到tmin,方可達到終止條件。
SA偽代碼如下:
求解TSP問題簡述
SA算是一種通用的優化演算法,理論上具有概率的全局優化性能,所以在以下求解的流程採用方法為,優先生成一距離矩陣(distance matrix),此矩陣是求解最佳路徑中,所需點中彼此間的距離,並以此矩陣為基礎,在每個點不重複的情況下,求出一總距離最短的路徑。
其整體流程撇除演算法的核心步驟,整體還是和筆者一篇文章,人工蜂群演算法求解旅行推銷員問題,原理上大致相同,有興趣的讀者可以參考看看,裡面針對如何使用實數解求解TSP有詳細的分析過程。
以下正文開始~
基本參數設定
- D:10。城市點數量,以下完整路徑會以長度(D+1)表示。
- N: 10。和 D 相同為城市點數量,根據 N 決定生成 N*N 的距離矩陣。
- t0:500。起始溫度。
- tmin:0.1。最低溫度,意味著迭代結束的溫度點 。
- k:40。每次降溫階段中內循環的運算次數。
- coolnum : 0.98。於每次迭代過程中的衰退速率。
決定城市數量並生成距離矩陣(Distance Matrix)
首先,決定本次要求解的城市點數,參數 N 決定城市點為10,意指在10個城市間選擇出距離最短的路徑解法。可以看到以下是 10x10 的距離矩陣,分別代表第 i 個城市到第 j 個城市的距離,以及相對應的第 j 個城市到第 i 個城市的距離,且每個城市點對點間距離限制設定在 10~50 之間。(此距離限制可自行設定)
路徑問題初始化
本文求解TSP的限制因子為,透過上述參數 D,可知本次試驗的總城市數,但本文求解最後會回到一開始出發的城市。所以可以知道,透過上述距離矩陣(Distance Matrix)可知此次試驗城市點有10個,但是不管從第幾個城市出發最後都會回到那個出發的城市,因此實際上在算路徑時會有11個點。
以下為初始化路徑範例,隨機的結果為從第 4 個城市出發,最後一樣會回到第 4 個城市,配合距離矩陣(Distance Matrix)可知初始路徑總長為234。
產生新路徑
主要是在不改變起始和結尾城市的情況下,每次交換 2 個城市的先後順序來產生一新路徑。
def Inversion2(Route0,time):
Copy = Route0.copy()
place = random.sample([1,2,3,4,5,6,7,8,9],2)
Copy[place[0]],Copy[place[1]] = Copy[place[1]],Copy[place[0]]
return Copy
若是考慮的城市點較多,可採用隨機因子較大的方法生成一新路徑。
以下範例是以內循環運算次數 k 中的偶數次和奇數次分成不同的交換方法為主,其交換方法內容可以自訂。偶數次採隨機交換兩點;奇數次為二分法,設定一半機率交換前半部路徑的城市點,另一半交換後半部路徑的城市點。
def Inversion(Route0,time):
Copy = Route0.copy()
if time % 2 == 0:
place = random.sample([1,2,3,4,5,6,7,8,9],2)
else:
Randnum = random.randint(0,1)
if Randnum == 0:
place = random.sample([1,2,3,4],2)
else:
place = random.sample([5,6,7,8,9],2)
Copy[place[0]],Copy[place[1]] = Copy[place[1]],Copy[place[0]]
return Copy
Method of exhaustion
為判斷是否能求出最佳解,採用窮舉法(Method of exhaustion),將頭尾都是 4的路徑中全部都嘗試一遍。
在此不用特別設定頭尾的城市點,我寫的代碼會根據頭尾的城市而限制要窮舉的組合。
得出最佳路徑和最短距離為:
將上述提及參數實際丟入SA中運算得到的最佳解如下:
將迭代流程視覺化如下:
可看出這次試驗成功,最佳結果有迭代到窮舉法的最佳解。
完整代碼分享
結語
相較其他演算法,SA算是較易理解的演算法,入門實踐也較為簡單。透過練習這類簡單的演算法,除了可增加自己思考能力外,也可以增加 coding 實戰的能力,最後謝謝大家願意看完這篇文章,謝謝。