1.問題描述
旅行商問題(Travelling Salesman Problem, 簡記TSP,亦稱貨郎擔問題):設有n個城市和距離矩陣D=[dij],其中dij表示城市i到城市j的距離(i,j=1,2 … n),則問題是要找出遍訪每個城市恰好一次的一條回路並使其路徑長度為最短。
2.算法設計
對原問題進行分析,TSP的一個解可表述為一個循環排列:
Π= (Π1,Π2,Π3… Πn),即
Π1→ Π2→ … → Πn→ Π1
有(n-1)!/2 種不同方案,若使用窮舉法,當n很大時計算量是不可接受的。旅行商問題綜合了一大類組合優化問題的典型特征,屬於NP難題,不能在多項式時間內進行檢驗。若使用動態規划的方法時間復雜性和空間復雜性都保持為n的指數函數。
本次實驗利用模擬退火算法(Simulated Annealing)求解TSP問題。模擬退火算法最早由N.Metropolis等人於1953年提出,基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。該算法從某一較高初溫出發,伴隨溫度參數的不斷下降,結合概率突跳特性在解空間隨機尋找全局最優解。
退火是將固體加熱到足夠高的溫度,使分子呈隨機排列態,然后逐步降溫冷卻,最后分子以低能狀態排列,得到穩定狀態的固體。退火的過程有:
(1)加溫過程:增強粒子運動,消除系統原本可能存在的非均勻態;
(2)等溫過程:對於與環境換熱而溫度不變的封閉系統,系統狀態的自發變化總是朝向自由能減少的方向進行,當自由能達到最小時,系統平衡;
(3)冷卻過程:使粒子熱運動減弱並逐漸趨於有序,系統能量逐漸下降,從而得到低能的晶體結構。
其中,固體在恆溫下達到熱平衡的過程采用Metropolis方法進行模擬:
溫度恆定為T時,當前狀態i轉為新狀態j,如果j狀態的能量小於i,則接受狀態j為當前狀態;否則,如果概率p=exp{-(Ej-Ei)/(k*T)}大於[0,1)區間的隨機數,則仍接受狀態j為當前狀態;若不成立則保留狀態i為當前狀態。
溫度變化時,由於p=exp{-(Ej-Ei)/(k*T)},因此在高溫下,可接受當前狀態能量差較大的新狀態;在低溫下,只接受與當前狀態能量差較小的新狀態。
退火過程由冷卻進度表(Cooling Schedule)控制,包括控制參數的初值t及其衰減因子Δt、每個t值時的迭代次數L和停止條件S。
我們使用模擬退火算法實現TSP的參數選取如下:
控制參數初值:t0 = 2000
停止准則:連續2個Mapkob鏈中對路徑無任何變動(優化或惡化的)時即停止算法運行。
冷卻進度表中的控制參數t的衰減函數:α(t)=0.98 * t
Mapkob鏈長:定長20000
新解的產生:隨機交換兩個序號。任選序號u和v (u<v),將u和v的順序交換。
(Π1 …Πu-1ΠuΠu+1 …Πv-1ΠvΠv+1… Πn)
變為
(Π1 …Πu-1ΠvΠu+1 …Πv-1ΠuΠv+1… Πn)
3.程序流程
4.核心偽代碼
Begin:
INITIALIZE; { 初始化i0= Π1 …Πn, t0=2000, L= 20000 }
randomize; { 初始化隨機數種子}
Repeat
bChange:=0;
for l:=1 to L do
begin:
GENERATE; { 隨機交換兩個序號產生新的路徑}
CALCULATE(df); { 計算df= f(j)-f(i) 的值}
if ACCEPT(t, df) then { 根據Metropolis准則判斷是否接受新的路徑}
begin
f: = f + df; { 計算已接受的路徑的長度}
bChange:= 1;
end;
end;
t:=t * 0.9; { 衰減函數α(t)=0.98 * t }
if (bChange= 0) then s:=s+1 else s:=0;
until s = 2 { 停止准則為連續兩個Mapkob鏈路徑無變化}
End;
5.代碼運行及測試
運行環境:VS2013,Windows 10系統
輸入說明:
第一行輸入需要遍歷的城市數目num,之后的num行輸入鄰接矩陣用於描述城市兩兩之間的距離。
輸出說明:輸出包括最短路徑回路和最短路徑的距離。
程序運行結果:
(1)
程序輸入:
程序輸出1:
程序輸出2:
程序輸出3:
(2)
程序輸入:
程序輸出1:
程序輸出2:
程序輸出3:
由上述實驗結果可知,如果全局最優解不唯一,那么對於同一組輸入數據,多次運行得到的路徑也不相同。SA算法的效果依賴於鏈長、溫度變化范圍、控制參數等等,如果參數選取不合適將不能得到好的解。此外,實驗過程中發現,在當前的參數(t0=2000,L=20000,alpha=0.98)情況下,問題規模很小時耗費時間較長,而當問題規模較大時耗費時間明顯減少,可見模擬退火算法在大規模問題下的高效性。
6.結論
本次實驗采用模擬退火算法解決TSP問題,該算法具有超越算法本身的啟發式意義。模擬退火算法將自然物理過程中的固體退火過程與組合優化問題進行類比,成功地將退火過程遷移用於組合優化問題,不僅應用范圍廣、使用靈活,而且初始化條件少。但是收斂速度較慢,運行時間較長。模擬退火算法是在狀態空間搜索近似最優解的過程,由於搜索過程中會以一定概率接受更差的結果,所以能夠跳出局部最優狀態范圍。
源碼:
#include <iostream> #include <time.h> #include <math.h> using namespace std; #define MAX_CITY_NUM 100 #define T0 2000 #define T 1e-5 #define ALPHA 0.98 #define L 20000 struct path{ int route[MAX_CITY_NUM]; double dis; }; double w[MAX_CITY_NUM][MAX_CITY_NUM]; int num; double s = 0; path p0; void Init_City() { cout << "Please input the number of cities:" << endl; cin >> num; for (int i = 0; i < num; ++i) for (int j = 0; j < num; ++j) cin >> w[i][j]; } void Init_path() { p0.dis = 0; for (int i = 0; i < num; ++i) { p0.route[i] = i; if (i != num - 1) p0.dis += w[i][i + 1]; } p0.dis += w[num - 1][0]; } path generate(path p) { int x = 0, y = 0; while (x == y) { x = (int)(num * (rand() / (RAND_MAX + 1.0))); y = (int)(num * (rand() / (RAND_MAX + 1.0))); } path gen = p; int tmp; tmp = gen.route[x]; gen.route[x] = gen.route[y]; gen.route[y] = tmp; gen.dis = 0; for (int i = 0; i < num - 1; ++i) gen.dis += w[gen.route[i]][gen.route[i + 1]]; gen.dis += w[gen.route[num - 1]][gen.route[0]]; return gen; } void TSP_SA() { double t = T0; srand(time(NULL)); path cur = p0; path next = p0; int bChange; while (t > T) { bChange = 0; for (int i = 0; i < L; ++i) { next = generate(cur); double df = next.dis - cur.dis; if (df <= 0) { cur = next; bChange = 1; } else { double rndp = rand() / (RAND_MAX + 1.0); double eps = exp(-df / t); if (eps > rndp && eps < 1) { cur = next; bChange = 1; } } } if (cur.dis < p0.dis) p0 = cur; t *= ALPHA; if (!bChange) ++s; else s = 0; if (s == 2) break; } } int main() { Init_City(); Init_path(); TSP_SA(); cout << "The best path is:" << endl; for (int i = 0; i < num; ++i) cout << p0.route[i] << " -> "; cout << p0.route[0] << endl; cout << "The distance of the best path is: " << p0.dis << endl; system("pause"); return 0; }