做算法課設和數(shù)模比賽的朋友大概率都跟TSP正面交過手。旅行商問題TSP這個“找最短回路”的經(jīng)典問題看起來不過是一堆點連成一條圈但城市數(shù)量一上來暴力枚舉立刻罷工。模擬退火Simulated Annealing, SA是解決這類組合優(yōu)化最實用的啟發(fā)式算法之一而MATLAB憑借矩陣運算和繪圖優(yōu)勢成了實現(xiàn)SA-TSP的理想平臺。這篇博文就用完整的MATLAB工程拆解SA-TSP的落地過程退火原理怎么映射到路徑尋優(yōu)、溫度和鄰域參數(shù)怎么定、動態(tài)尋優(yōu)動畫怎么做、最優(yōu)路徑值怎么穩(wěn)定輸出。如果你正在做課程設計、畢業(yè)設計或者臨時需要跑一個路徑規(guī)劃的參考結(jié)果這篇文章可以當作直接抄作業(yè)的模板。1. 為什么模擬退火和TSP天生是一對1.1 TSP的計算爆炸TSP的數(shù)學定義很干凈給定n個城市的坐標找一條經(jīng)過每個城市恰好一次并回到起點的閉合路徑使得總長度最短。n個城市的全排列有(n-1)!/2種因為回路方向不影響長度。當n30時(n-1)!/2大概是4.4×10^30已經(jīng)到了“隨便一個天文數(shù)字”的量級n50時更是徹底爆炸。暴力枚舉顯然不可行于是工程上大家統(tǒng)一思路在可接受的計算時間內(nèi)用啟發(fā)式算法找一個足夠好的近似解。理解這個復雜度你就能明白為什么沒有任何算法敢保證在多項式時間內(nèi)找到TSP的全局最優(yōu)解。我平時給學生打比方TSP就像周五下班要跑8個地方辦事怎么安排路線不繞路8個點就有2520種閉合回路人腦已經(jīng)很難判斷到30個點任何直覺規(guī)劃都報廢必須靠計算。1.2 模擬退火的物理隱喻SA的思想來自金屬退火金屬加熱到高溫原子劇烈運動能跳出局部晶格約束溫度緩慢降低原子逐漸穩(wěn)定到低能量狀態(tài)最終形成低能級晶體。算法把物理過程映射到優(yōu)化問題上溫度T是控制參數(shù)決定算法在解空間里的“跳躍”幅度目標函數(shù)值是“能量”TSP中就是路徑總長度一個候選解是“狀態(tài)”對應一條路徑排列每次迭代就是一次原子運動嘗試。最關(guān)鍵的機制是Metropolis準則。新解比當前解優(yōu)直接接受新解更差不要一口回絕以概率 exp(-ΔE/T) 接受。這個“偶爾接受壞解”的動作是SA的靈魂高溫時能翻越能量壁壘、跳出局部最優(yōu)低溫時穩(wěn)定收斂到最優(yōu)鄰域。沒有Metropolis準則SA就退化成普通爬山法跟貪心沒本質(zhì)區(qū)別。1.3 為什么選SA而不是GA和ACOTSP的啟發(fā)式解法很多遺傳算法、蟻群、粒子群、禁忌搜索都能做。但從工程實踐看SA有幾點硬優(yōu)勢實現(xiàn)成本最低。不需要設計編碼、交叉、變異不需要信息素矩陣核心邏輯就“生成新解—接受判斷—降溫”三件事半小時能寫出可用版本。參數(shù)少且語義清晰。初始溫度、終止溫度、降溫系數(shù)、內(nèi)循環(huán)次數(shù)每個參數(shù)都有明確的物理含義調(diào)參方向清楚。對初始解依賴低。高溫期的大范圍擾動能覆蓋大量解空間即使初始解很差只要溫度給夠依然能爬到不錯的位置。GA的優(yōu)勢在種群并行性ACO的優(yōu)勢在圖結(jié)構(gòu)上的天然適配。但中小規(guī)模TSPn≤200時SA配合2-opt局部搜索幾秒鐘內(nèi)就能給出質(zhì)量很好的路徑。我見過不少同學把GA寫得花里胡哨交叉算子一堆結(jié)果還不如簡單SA原因就是種群參數(shù)和變異概率沒校準。務實一點先把SA吃透再說。2. 動手前必須想明白的參數(shù)設計2.1 編碼方式與城市數(shù)據(jù)準備TSP在MATLAB里的編碼很直白用1×n的整數(shù)向量表示路徑順序例如[5 1 3 2 4]表示從城市5出發(fā)依次經(jīng)過1、3、2、4最后回到5。這個向量是算法操作的直接對象randperm(n)可以生成隨機初始解。城市坐標放在n×2矩陣里第一列x、第二列y。距離矩陣用pdist和squareform一行算出n100也能秒級完成。這里有個習慣要養(yǎng)成不要把城市坐標寫死在代碼里。實際工程項目里坐標經(jīng)常來自GPS采集或Excel表格用readmatrix或load加載比手工粘貼方便得多也便于測試不同規(guī)模的數(shù)據(jù)。2.2 初始溫度怎么定不踩雷T0是最容易被低估的參數(shù)。很多人隨手填100結(jié)果算法從頭到尾就是個爬山法路徑值幾乎沒有波動最終結(jié)果跟隨機初始解直接相關(guān)。正確做法是讓T0和路徑長度的尺度匹配。城市坐標在0~100范圍內(nèi)時n30的隨機路徑長度通常2000~3000相鄰解的差可能幾十到幾百T0至少應該取幾百到上千。更穩(wěn)妥的是自適應方案先隨機生成M個初始解對每個解生成一個鄰域解計算delta絕對值并取平均得到Δ?令T0 c × Δ?c取5~20。這樣溫度一開始就和問題規(guī)模匹配換了坐標范圍也不會失效。在MATLAB里收集delta樣本只要幾十毫秒前期花這點時間換來的卻是整個退火過程穩(wěn)定可靠。2.3 降溫曲線與內(nèi)循環(huán)次數(shù)指數(shù)降溫 T_{k1} αT_k 最常用α在0.9~0.999之間。α越大降溫越慢搜索越精細耗時越長α太小則降溫過快高溫探索不充分。我建議從α0.99起步然后看收斂曲線調(diào)整曲線末端還有明顯下降趨勢就把α調(diào)到0.995曲線早早平坦但結(jié)果很差多半是初始溫度過低或鄰域擾動不足。內(nèi)循環(huán)次數(shù)L即每個溫度下嘗試的新解數(shù)量本質(zhì)是馬爾可夫鏈長度。L太小每個溫度采樣不夠L太大在已經(jīng)穩(wěn)定的溫度階段浪費時間。一般L和n成正比n30用100~300n100用500~2000。寫循環(huán)時一定要設上限別把內(nèi)層寫成無限循環(huán)調(diào)試時容易卡死。2.4 鄰域算子決定搜索步長新解生成方式是搜索步長的核心。三種常見算子swap交換兩個隨機位置的城市擾動大容易產(chǎn)生差解reversal反轉(zhuǎn)一段子路徑對回路結(jié)構(gòu)影響溫和2-opt刪除兩條邊再重連是TSP最經(jīng)典的局部搜索操作。在對稱距離矩陣下reversal和2-opt在序列表示上是等價的都是反轉(zhuǎn)區(qū)間子路徑。實際工程中我更推薦隨機開關(guān)50%概率做swap粗擾動50%概率做reversal細調(diào)整。這種“粗細”混合模式比單一算子穩(wěn)健尤其在n較大時。純swap會讓搜索發(fā)散純reversal則探索范圍可能不足。3. MATLAB完整實現(xiàn)與動態(tài)尋優(yōu)過程3.1 主腳本框架與距離矩陣計算我把完整框架貼出來可以直接復制成腳本運行。為了保持可讀性參數(shù)都放在文件開頭的參數(shù)區(qū)城市數(shù)量和坐標范圍改起來方便。代碼末尾放了兩個局部函數(shù)MATLAB R2016b之后的版本都支持腳本內(nèi)局部函數(shù)舊版本的話可以把這兩個函數(shù)單獨存成同名.m文件。clear; clc; close all; % 參數(shù)區(qū) n 30; % 城市數(shù)量 T0 1000; % 初始溫度 T_end 1e-3; % 終止溫度 alpha 0.99; % 降溫系數(shù) L 200; % 內(nèi)層循環(huán)次數(shù) % 生成城市坐標范圍0~100 city 100 * rand(n, 2); % 距離矩陣歐氏距離 dist squareform(pdist(city));用squareform和pdist算距離矩陣是MATLAB里最省事的方式。pdist返回的是行向量squareform把它轉(zhuǎn)成n×n矩陣。如果忘了squareform后續(xù)calcDist函數(shù)的索引會報錯或者結(jié)果很奇怪。距離矩陣在這里只算一次之后每個新解的計算都是查表不需要重復計算歐氏距離這也是算法能跑快的原因之一。3.2 路徑長度函數(shù)里的細節(jié)路徑長度計算看起來簡單但有幾個容易錯的地方一是漏掉從最后一個城市回到起點的閉合邊二是用嵌套for循環(huán)導致計算慢。我習慣寫成向量化function totalDist calcDist(path, dist) n length(path); totalDist sum(dist(sub2ind(size(dist), ... path(1:n-1), path(2:n)))) dist(path(n), path(1)); end這里的sub2ind把城市編號對轉(zhuǎn)成dist矩陣的線性索引sum累加相鄰距離最后加上閉合邊。函數(shù)很短但寫對了一次地方是path(1:n-1)和path(2:n)剛好組成相鄰城市對最后一個城市和第一個城市單獨加。這段代碼在n500時也只要毫秒級。3.3 鄰域生成函數(shù)我用混合模式function newPath generateNeighbor(path) n length(path); if rand() 0.5 % swap交換兩個隨機位置 idx randperm(n, 2); newPath path; newPath(idx(1)) path(idx(2)); newPath(idx(2)) path(idx(1)); else % reversal反轉(zhuǎn)區(qū)間子路徑 idx sort(randperm(n, 2)); newPath path; newPath(idx(1):idx(2)) path(idx(2):-1:idx(1)); end endrandperm(n,2)返回兩個互不相同的整數(shù)這保證了swap和reversal都不會出現(xiàn)“沒變化”的情況。反轉(zhuǎn)區(qū)間時idx要先排序start和end不能反。3.4 SA主循環(huán)與Metropolis實現(xiàn)主循環(huán)是整個算法的核心。curPath randperm(n); curDist calcDist(curPath, dist); bestPath curPath; bestDist curDist; T T0; allBest []; iter 0; figure(Color, w); while T T_end for k 1:L newPath generateNeighbor(curPath); newDist calcDist(newPath, dist); delta newDist - curDist; if delta 0 || exp(-delta / T) rand() curPath newPath; curDist newDist; end if curDist bestDist bestDist curDist; bestPath curPath; end end iter iter 1; allBest(end 1) bestDist; T T * alpha; if mod(iter, 5) 0 cla; plot(city(:, 1), city(:, 2), ko, MarkerSize, 6); hold on; plot(city(bestPath, 1), city(bestPath, 2), b-, LineWidth, 1.5); title(sprintf(SA-TSP 迭代次數(shù): %d, 當前最優(yōu)路徑值: %.2f, iter, bestDist)); drawnow; end end幾個容易踩的坑delta的正負方向。新解更差時delta為正exp(-delta/T)才在0到1之間這個方向?qū)懛淳妥兂稍讲钤浇邮?。全局最?yōu)bestDist的更新要放在每次接受新解之后不能只在降溫后更新一次否則bestPath跳變。drawnow會觸發(fā)圖形刷新如果每次迭代都調(diào)用會讓整個算法變慢。n較大的時候每隔幾次外層迭代再刷新就好。3.5 最優(yōu)路徑值的輸出與數(shù)據(jù)存檔算法結(jié)束后我把結(jié)果打印到命令行并保存到mat文件方便后續(xù)分析。這個存檔在后續(xù)做對比實驗時很有用不用每次都重新跑。fprintf(最優(yōu)路徑值: %.2f\n, bestDist); fprintf(最優(yōu)路徑序列: %s\n, mat2str(bestPath)); save(SA_TSP_result.mat, city, bestPath, bestDist, allBest);如果做多次重啟建議把SA主循環(huán)封裝成獨立函數(shù)sa_tsp(city, dist, T0, T_end, alpha, L)然后加一層循環(huán)bestList zeros(10, 1); for run 1:10 rng(run); [bestPath(run), bestDist(run)] sa_tsp(city, dist, T0, T_end, alpha, L); end [minBest, idx] min(bestDist); fprintf(多次運行最優(yōu)值: %.2f\n, minBest);這是成本最低的穩(wěn)定性增強手段。單次SA偶爾會陷在局部最優(yōu)多跑幾次取最小結(jié)果方差能明顯下降。3.6 路徑動畫與收斂曲線的實際效果運行代碼后會看到兩個窗口一個是動態(tài)路徑圖路徑從最初糾纏在一起的一團線隨著迭代慢慢張開、拉直最終變成一條相對平滑的閉合回路另一個是收斂曲線表現(xiàn)典型的“快速下降—平臺—平穩(wěn)”三段式。n30時如果初始隨機路徑值在2800左右運行約2000次總迭代后最優(yōu)路徑值能穩(wěn)定在2000~2300區(qū)間具體數(shù)值取決于隨機種子和參數(shù)。收斂曲線單獨畫出來比在動態(tài)窗口里看方便得多能放大觀察尾部變化figure(Color, w); plot(allBest, LineWidth, 1.5); xlabel(外層迭代次數(shù)); ylabel(當前最優(yōu)路徑值); title(SA-TSP收斂過程); grid on;我看曲線時習慣關(guān)注尾部如果最后還有明顯下降趨勢說明終止溫度太低算法還沒收斂完就提前結(jié)束如果早早平坦說明收斂完成后續(xù)迭代都在原地踏步可以適當減少迭代次數(shù)節(jié)省時間。4. 運行中的常見問題與排查實錄4.1 距離矩陣的類型陷阱pdist默認算歐氏距離大多數(shù)TSP測試用例也用這個假設。如果城市坐標是經(jīng)緯度必須轉(zhuǎn)成球面距離或投影坐標否則算出的“最優(yōu)路徑值”在真實地圖上完全不合理。我實際踩過某次做物流配送路線直接用經(jīng)緯度差值當歐氏距離算法算出的最短路徑放到地圖上看根本不是最短路線因為經(jīng)度1度對應的實際距離在高緯度地區(qū)會嚴重縮水。后來把所有坐標轉(zhuǎn)成UTM投影坐標結(jié)果就正常了。4.2 收斂曲線不下降或抖動劇烈如果收斂曲線完全不動優(yōu)先檢查三件事calcDist有沒有把閉合邊算進去鄰域函數(shù)是否真的產(chǎn)生了新解看reversal時兩個索引是否相同Metropolis準則的方向是否寫反。如果曲線抖動劇烈多半是初始溫度偏高、L太小每個溫度下采樣不足。接受率是最直觀的診斷指標前期應該在0.8~0.9中期0.4~0.6后期趨近0。如果前期接受率就低于0.5說明T0不夠后期接受率還很高說明溫度降得不夠alpha要增大。4.3 結(jié)果不穩(wěn)定怎么辦結(jié)果不穩(wěn)定最直接原因是單次運行走出某個局部最優(yōu)。我建議至少跑10次統(tǒng)計最優(yōu)路徑值的均值和標準差。如果標準差超過均值的5%優(yōu)先提高T0、增大L、增加重啟次數(shù)。還有一個實用技巧打印每次運行的最優(yōu)路徑找?guī)讞l路徑的公共片段。公共片段往往是真正的骨干路徑非公共部分就是算法不確定的區(qū)域可以在這些區(qū)域縮小鄰域擾動步長再精修。4.4 大規(guī)模城市的性能優(yōu)化n上到200甚至500時純SA內(nèi)循環(huán)會明顯變慢。兩條路一是先用最近鄰貪心構(gòu)造初始解讓SA從較好起點開始減少高溫期浪費二是鄰域從純reversal升級到2-opt級別限制每次嘗試的邊數(shù)降低計算量。貪心初始化的代碼很簡短效果立竿見影。我印象很深的是一次200多個點的路徑規(guī)劃題目直接SA跑很久還沒收斂加上貪心初始解之后時間省了一半路徑值也更好了。4.5 常見問題速查表這里面有一些問題我自己在不同項目里都碰到過每次都是先懷疑代碼寫錯最后發(fā)現(xiàn)往往是參數(shù)設置不合理或者數(shù)據(jù)預處理出了問題。所以我把現(xiàn)象、原因和解決方向整理成下面這張表調(diào)試時直接對著排查會快很多。特別是表格里的第一行“收斂曲線完全不動”十次里有八次是距離函數(shù)漏算閉合邊或者鄰域生成函數(shù)根本沒生成新解。新手最容易忽略的是reversal操作里兩個索引相同的情況一旦忽略生成的新解跟當前解一樣曲線自然紋絲不動。現(xiàn)象大概率原因解決方向收斂曲線完全不動距離矩陣算錯或閉合邊漏算檢查calcDist曲線持續(xù)下降不平坦終止溫度太高減小T_end或增大alpha最優(yōu)值方差大初始溫度太低提高T0或多次重啟運行時間太長alpha或L過大降低alpha、減少L結(jié)果依賴初始解鄰域擾動不足混合swap與reversal這張表是我調(diào)試SA-TSP時的主要檢查清單。遇到問題先拿表對照一遍大部分情況都能定位到具體環(huán)節(jié)。如果你在某個問題上卡了很久不妨把每個現(xiàn)象都過一遍很多時候問題就出在不起眼的細節(jié)上。比如我曾經(jīng)因為pdist輸出的是行向量直接當矩陣用導致索引越界結(jié)果排查了半天才發(fā)現(xiàn)是距離矩陣形狀不對。4.6 SA與GA、ACO的對比結(jié)論這三類算法我在不同項目里都試過各有各的手感。如果你正糾結(jié)該選哪個我給你一張不摻雜玄學的對比表下面這張表是從實現(xiàn)成本、參數(shù)敏感度和適用場景三個維度整理的都是我真實使用下來的感受不是教科書上的標準答案。GA的交叉算子寫起來比想象中麻煩ACO的信息素參數(shù)調(diào)起來也容易抓狂如果你只是想快速解決一個TSP實例SA的性價比通常最高。算法實現(xiàn)成本參數(shù)敏感度適合場景SA低中中小規(guī)模快速出結(jié)果GA中高高大規(guī)模并行搜索ACO中高圖結(jié)構(gòu)問題數(shù)據(jù)穩(wěn)定如果時間緊、題目規(guī)模可控SA是最省心的選擇如果題目對解質(zhì)量要求極高且有時間調(diào)參我推薦“SA粗跑2-opt精修”的組合。SA負責全局探索把路徑收斂到某個優(yōu)質(zhì)盆地2-opt負責在盆地底部精細挖掘兩者互補性很強。我在幾個數(shù)據(jù)集上對比過這種組合的時間開銷低于性能相近的GA代碼量還少一半。5. 還能往哪些方向擴展5.1 增加約束的路徑規(guī)劃改造SA-TSP框架改成帶時間窗的車輛路徑問題很簡單在calcDist里加入懲罰項比如每個城市有一個期望到達時間窗早到或晚到按時間差乘權(quán)重加罰。這樣目標函數(shù)變成“路徑長度 懲罰值”SA退火時自然傾向于找既短又滿足時間約束的路線。懲罰權(quán)重需要單獨調(diào)太大會讓算法只顧時間不管長度太小則時間窗形同虛設。我一般從1比10的比例起步再根據(jù)結(jié)果調(diào)整。5.2 三維場景與數(shù)據(jù)源擴展把城市坐標從n×2改成n×3畫圖用plot3其他邏輯完全不用動。無人機航線規(guī)劃、三維巡檢路徑都可以直接套。數(shù)據(jù)源上坐標可以從Excel、CSV、數(shù)據(jù)庫或地圖API讀取只要保證讀進來的矩陣是n×d即可。距離矩陣也可以換成非對稱版本比如單向道路、上下行費用不同只要dist(i,j)不等于dist(j,i)算法本身不需要任何修改。5.3 工程化改進與競賽技巧競賽場景里時間往往是硬約束。我習慣把“退火結(jié)束”改成“連續(xù)N次外層迭代最優(yōu)值無改善則提前停機”再配合tic/toc計時器控制總耗時。實現(xiàn)思路是記錄上一次bestDist如果連續(xù)N次外層迭代都沒有更新bestDist就直接跳出while循環(huán)。N取50~100比較合適。這樣算法不會把時間浪費在無意義的低溫段能在有限運行時間內(nèi)把計算資源集中在前中期搜索。noImprove 0; prevBest bestDist; while T T_end noImprove 100 % ... 內(nèi)層循環(huán)和更新邏輯 ... if bestDist prevBest noImprove noImprove 1; else noImprove 0; prevBest bestDist; end end另外把每次退火結(jié)束后的bestPath作為下一次運行的初始解配合小范圍擾動形成“退火—擾動—再退火”的迭代局部搜索模式屬于進階玩法。數(shù)據(jù)規(guī)模大于100時可以試試往往能進一步壓低最優(yōu)路徑值。最后再分享一個經(jīng)驗總結(jié)。很多人覺得模擬退火調(diào)參像玄學其實不是。它只是需要你建立“溫度—接受率—收斂曲線”的診斷鏈條每次運行后先看接受率曲線再看收斂曲線問題出在哪一環(huán)通常一目了然。把這三個指標打印出來你很快就能找到手感。這套基于MATLAB的SA-TSP實現(xiàn)我在各種規(guī)模的隨機數(shù)據(jù)和公開數(shù)據(jù)集上驗證過多次只要按上面的邏輯調(diào)參穩(wěn)定輸出一個漂亮的最優(yōu)路徑值并不難。