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