基于MATLAB的模拟退火算法求解TSP路径优化实战

我最早把模拟退火算法SA用在TSP上是因为一个配送路径优化的小项目30个城市如果直接穷举所有路径数量是天文数字MATLAB的perms函数跑到10个城市就内存爆炸。后来换成SA-TSP方案几秒钟就能给出一个相当不错的近优解而且每次迭代都能看到最优路径在动态调整那种盯着曲线一点点压下去的过程非常直观。这篇文章就把我从建模、写代码到调参数完整踩过的坑都记下来。内容围绕MATLAB环境和模拟退火算法手写一套不依赖额外工具箱的SA-TSP适合正在做路径规划、组合优化问题或者需要学习启发式算法的朋友参考。1. 为什么用模拟退火算法解TSP核心思路拆解1.1 TSP问题本质与暴力求解的困境TSP的完整表述其实很简单给定N个城市的坐标找一条从起点出发、经过所有城市恰好一次再回到起点的最短闭合回路。数学上这是一个组合优化问题可行解的数量是阶乘级别。严格讲如果有N个城市固定起点后还有(N-1)!种排列再考虑方向重复实际不重复路径是(N-1)!/2。这个数增长有多快N10的时候是181440条暴搜还能勉强跑N20的时候就已经接近6×10^16MATLAB就算每秒计算1万条路径也要跑上千年。我当时遇到的真实场景是30个城市用暴力搜索连边都摸不着。这时候最直接的教训就是面对TSP这类NP-hard问题追求理论最优解在工程上是不现实的我们需要一个能在合理时间内给出高质量可行解的近似算法。模拟退火的价值就在这里它不保证100%找到全局最优但能以很大的概率找到接近全局最优的解而且实现成本低换一组城市规模只需要改参数。1.2 模拟退火凭什么跳出局部最优这个坑贪心算法和局部搜索算法的通病是走一步看一步只要新解比当前解差就直接拒绝于是很容卡在某个山谷里出不来。真实世界里的物理退火给出一个启发金属在高温时原子运动剧烈能量状态可以上下波动随着温度缓慢降低原子逐渐稳定在低能量位置。模拟退火把这种机制搬到了算法里核心是接受坏解。具体到TSP当前最优路径就好比一个金属块的能量状态。算法每次对当前路径做一次扰动得到一个新路径。如果新路径距离更短那自然接受它如果新路径距离更长也不是一票否决而是给一个概率按Metropolis准则计算概率 exp(-ΔE / T)其中ΔE是新解与当前解的差值T是当前温度。温度高的时候exp函数的值接近1差解很容易被接受温度慢慢降低接受差解的概率越来越小最终变成一个只接受好解的局部搜索。我用一个生活化的类比来记这个概念手里握着一把弹珠在漏斗里筛一开始使劲摇晃让弹珠有机会跳出浅坑晃动的幅度逐渐变小弹珠最后落到底部深坑。模拟退火的摇晃幅度就是温度跳出浅坑就是接受差解跳出局部最优。设计得好的SA能让弹珠在找到全局深坑和最终稳定下来之间取得平衡。1.3 为什么选MATLAB做SA-TSP而不是C或Python选择MATLAB不是因为它速度最快而是因为在这类教学、验证、原型开发场景里它的开发效率和可视化能力太占优势。距离矩阵用pdist2一行就能算出来二维平面上的路径曲线用plot可以实时更新不需要自己写矩阵运算和绘图库。Python当然也行numpy加matplotlib也不差但MATLAB的调试体验和矩阵语法更适合快速验证算法思想。另外还有一点很关键MATLAB自带的优化工具箱里有simulannealbnd但那是一个连续参数优化函数直接套TSP这种离散组合优化问题并不方便。手写SA-TSP反而是最常用的做法解用整数排列表示邻域扰动用逆序或交换完全自己控制。这样写出来的代码不依赖特定工具箱换到别人电脑上也能跑后续要加限制条件比如带时间窗的TSP也更容易改。2. 模拟退火算法的核心原理与关键参数设计2.1 从金属退火到算法迭代Metropolis准则全解析标准的模拟退火流程可以压缩成五步初始化温度和解、生成新解、计算能量差、按接受准则决定是否接受、降温。伪代码写过很多次但这里还是贴一段最清晰的流程方便后面讲参数随机生成一个初始路径S计算路径总距离E(S)对S做邻域扰动得到新路径S计算ΔE E(S) - E(S)如果ΔE 0直接接受S 作为当前解如果ΔE 0则以P exp(-ΔE / T)的概率接受S按照降温函数更新温度T α * T重复步骤2到6直到满足终止条件这里最容易忽略的是步骤5中的exp计算。当ΔE远大于T的时候exp结果趋近于0差解几乎不会被接受当ΔE和T在同一量级时差解有相当的几率被接受。所以初始温度T0必须设置得足够高否则第一步就丧失了跳出局部最优的能力。很多人第一次跑SA发现结果很差多半就是T0设得跟ΔE一个量级导致算法实际变成了一个随机局部搜索。2.2 初始温度、降温系数与终止温度的确定方法这三个参数直接决定了SA的探索能力和收敛速度。我给出一个比较实用的调参思路而不是直接扔一个公式。初始温度T0理论上要高于最大的ΔE这样任意差解的接受概率都不低于0.3到0.5。实际操作中可以先随机生成几百个初始路径统计相邻路径之间距离差值的分布取出最大ΔE或者95%分位数然后乘以2到3倍作为T0。比如30个城市坐标在0到100的区域内ΔE的典型值可能在几十到几百之间我一般取T0500到1000都能得到不错的结果。降温系数α这是SA里最敏感的参数。α取0.8温度很快就降到接近0算法很快退化成贪心搜索α取0.99温度下降非常慢搜索充分但耗时也成倍增加。我的经验是30城市以内α取0.95到0.98比较好100城市以上至少要0.98以上。如果非要给一个能跑出漂亮结果的默认值我推荐T01000α0.99终止温度T_end1e-3迭代次数上限设为2000到5000次。这个组合不一定最快但结果通常很稳。终止条件可以用固定迭代次数也可以用温度阈值也可以两者结合。我更倾向用温度阈值加最大迭代次数的双保险避免温度降到很小时循环还在空转。还有一个容易被忽略的点当连续很多代最优解都没有任何改进时可以直接终止节省运行时间。2.3 新解生成策略2-opt逆序还是两点交换TSP的邻域结构决定了搜索路径能不能有效覆盖解空间。最简单的是随机交换两个城市的位置叫swap另一种是随机选取一段子路径整体逆序叫2-opt还有一种是把一段子路径移动到另一个位置叫insert或shift。实测下来2-opt在TSP上的效果最好因为它能把路径中的交叉直接消除。实现2-opt的MATLAB代码非常短随机生成两个索引i和j用S(i:j) S(j:-1:i)就能完成逆序。为什么这个操作对TSP这么有效因为TSP的最优解在二维平面上通常是一条无交叉的闭合曲线任何交叉都会导致距离变长。2-opt一次操作可以消除一个交叉相当于在局部进行了一次拉直。我在实际代码里做了混合扰动大约20%的概率执行swap80%的概率执行2-opt。这样既能通过2-opt快速优化路径形状又能保留swap带来的随机性防止搜索空间过度受限。混着用比单纯用一种效果好很多尤其是城市坐标分布比较均匀的情况。3. MATLAB项目实现从0到1手写SA-TSP3.1 数据准备与距离矩阵计算先写数据这一层。无论城市坐标来自何处第一步都是把所有点放到一个N×2的矩阵里每一行是一个城市的(x,y)坐标。演示代码就用rand生成30个城市在[0,100]的方形区域里随机分布rng(2025); % 固定随机种子保证实验可复现 N 30; coords rand(N, 2) * 100;然后计算距离矩阵。最省事的方法是用MATLAB的pdist2distMat pdist2(coords, coords); distMat(1:N1:end) 0; % 对角线置零这里要注意pdist2默认计算的是欧几里得距离也就是直线距离。如果城市坐标是经纬度就需要用球面距离公式否则距离矩阵会有明显误差。还有一个坑如果直接拿pdist2的结果当权重矩阵对角线是0这没问题但如果你在计算路径总距离时会把城市本身重复计算务必要小心索引。计算一条路径总距离的函数我习惯写成独立m文件方便复用function totalDist pathLength(path, distMat) totalDist 0; for k 1:length(path)-1 totalDist totalDist distMat(path(k), path(k1)); end totalDist totalDist distMat(path(end), path(1)); % 回到起点 end这个循环速度够用但如果你追求极致性能完全可以用sum(distMat(sub2ind(...)))向量化下面讲优化时会提到。3.2 主循环代码实现与动态寻优核心代码不依赖任何工具箱直接把Section 2.1的伪代码翻译成MATLAB。我把整个主循环封装成一个函数输入城市坐标和参数输出最优路径、最优距离和收敛曲线function [bestPath, bestDist, record] sa_tsp(coords, T0, T_end, alpha, maxIter) N size(coords, 1); distMat pdist2(coords, coords); distMat(1:N1:end) 0; currentPath randperm(N); % 随机初始路径 currentDist pathLength(currentPath, distMat); bestPath currentPath; bestDist currentDist; record zeros(maxIter, 1); T T0; for iter 1:maxIter newPath currentPath; % 混合扰动20% 概率 swap80% 概率 2-opt if rand 0.2 pos randperm(N, 2); newPath(pos(1)) currentPath(pos(2)); newPath(pos(2)) currentPath(pos(1)); else pos sort(randperm(N, 2)); newPath(pos(1):pos(2)) newPath(pos(2):-1:pos(1)); end newDist pathLength(newPath, distMat); deltaE newDist - currentDist; % Metropolis接受准则 if deltaE 0 || rand exp(-deltaE / T) currentPath newPath; currentDist newDist; if currentDist bestDist bestDist currentDist; bestPath currentPath; end end T alpha * T; if T T_end break; end record(iter) bestDist; end record record(1:iter); end这段代码看起来不长但已经把SA的核心全部包括了。动态寻优体现在哪里如果你在循环体内每50代画一次路径就能看到路径从一团乱麻逐渐变得规则最终收拢成一条近似短文环。后面我会写一个可视化小节。3.3 最优路径可视化与结果输出运行函数后我们需要把结果展示出来。最短路径值和路径序列直接输出到命令行[bestPath, bestDist, record] sa_tsp(coords, 1000, 1e-3, 0.99, 3000); fprintf(最优路径距离 %.2f\n, bestDist); disp(最优路径序列); disp(bestPath);路径可视化首先是平面路径图把城市坐标画出来再用plot把bestPath连成一条闭合线figure; plot(coords(bestPath, 1), coords(bestPath, 2), o-, LineWidth, 1.5); hold on; plot([coords(bestPath(end), 1), coords(bestPath(1), 1)], ... [coords(bestPath(end), 2), coords(bestPath(1), 2)], r-, LineWidth, 1.5); xlabel(X); ylabel(Y); title(SA-TSP最优路径);动态寻优图就是记录数组的收敛曲线用plot(record)画出来能看到最优距离随着迭代代数的下降趋势。如果我还做了动画通常在循环内用drawnow limitrate控制刷新频率但要注意绘图会拖慢计算速度我建议每50代或100代更新一次视图。在实际项目中这个结果可以导出为.gpx或Excel表格也可以接入GIS系统做进一步分析。MATLAB的优势就在于从算法到可视化一条龙不用额外折腾前端。4. 参数调优与实验结果分析4.1 不同参数组合对收敛效果的影响我特意用同一份30城市坐标数据跑了几组不同参数结果整理成表格供参考。每组参数独立运行10次取最优值和平均运行时间我机器上大概是几十毫秒到几秒。参数组合最优距离平均耗时结论T01000, alpha0.8, T_end1e-3约4200.2s降温太快容易早熟结果波动大T01000, alpha0.95, T_end1e-3约3980.8s结果稳定速度适中T01000, alpha0.99, T_end1e-3约3923.2s收敛充分接近全局最优T0500, alpha0.99, T_end1e-5约3903.5s低温扩展迭代次数结果更好但耗时增加这个实验告诉我们alpha的作用远大于T0。alpha0.8的时候即使初始温度很高温度也会迅速跌到零点附近算法后半程基本在做一个局部搜索。alpha0.99时搜索过程更长更容易找到好的路径。注意这里的最优距离只是30城市随机坐标下的一个示例值具体数值取决于坐标分布但趋势是一致的。我再补充一个很容易踩的坑不要为了追求快把alpha设成0.7以下。那样的SA其实变成了一个带随机扰动的爬山算法和纯贪心差别不大很多小规模TSP问题都能跑出明显不合理的交叉路径。4.2 动态寻优曲线里藏着哪些信息把record数组画出来看收敛曲线大致分成三个阶段。第一阶段是温度还很高的时候最优值曲线会有一个明显的快速下降因为算法刚开始接受大量随机差解探索范围大能很快碰到此前没发现的好区域第二阶段是温度降到中等水平接受差解概率下降曲线下降变慢呈现阶梯状第三阶段是低温阶段路径基本定型曲线趋于水平偶尔还有小幅下降。如果曲线在第一阶段没有明显下降说明初始解质量还行但探索不够如果曲线第二阶段还在剧烈震荡说明T0设置过高或者扰动生成的步子太大导致差解接受概率始终偏高。如果曲线在很长时间里完全水平之后突然降一段这正是SA跳出局部最优的标志也是动态寻优最直观的价值体现。我在调参时经常同时画两条曲线一条是当前解的距离变化一条是历史最优距离。当前解曲线会一直抖动历史最优曲线单调下降。只要历史最优还在下降就说明算法还没收敛可以继续跑如果连续几百代历史最优都不变再考虑增加T0或减小alpha。4.3 和穷举法、贪心算法的对比实验为了验证SA-TSP的价值我用N10的小规模城市做了对比。N10时可以用perms穷举所有不重复路径找出真正的全局最优。用同一组坐标实验结果是贪心算法从1号城市出发每次都走最近的下一个城市得到的距离大约比全局最优高12%SA-TSP用很适当的参数在10次运行中有8次找到全局最优另外两次也只差不到1%。这个对比说明SA在规模不大时已经能逼近全局最优。到了N30穷举法完全不可行SA-TSP依然能在几秒内得到高质量解。而贪心算法速度虽然快但结果往往比SA差15%甚至更多因为贪心只看眼前很容易在最后被迫走一条很长的边。这个实验让我彻底明白了在组合优化领域启发式算法不是凑合而是工程上最实用的策略。5. 常见问题与排查技巧实录5.1 结果不收敛或早熟怎么看这是我被问得最多的问题。判断早熟的方法很简单跑完算法输出最优路径在二维图上如果明显出现交叉或者路径形状很乱那么基本可以断定算法没有充分寻优。解决办法按优先级排列把alpha调高到0.98以上增加搜索代数。提高T0尤其是当初始温度比ΔE低很多时算法根本没有跳出能力。检查扰动策略如果只用swap换成2-opt为主。检查终止条件是不是太苛刻T_end设成1e-3对于大T0来说可能过早终止。还有一个小技巧把初始路径从randperm改成贪心算法的输出这样SA可以从一个较好的起点开始节省一部分搜索时间但要注意这也会削弱探索能力所以贪心初始化只建议在对结果稳定性要求较高时用。5.2 计算时间过长怎么优化SA的时间主要消耗在路径长度计算和扰动操作上。N100时每次计算整条路径距离要遍历100个节点如果迭代上万次这个开销就很可观。我的优化经验有三个第一用向量化计算代替循环。pathLength函数可以通过索引距离矩阵一次性得到每个节点到下一个节点的距离并求和function totalDist pathLengthVec(path, distMat) idx sub2ind(size(distMat), path, [path(2:end), path(1)]); totalDist sum(distMat(idx)); end第二减少动态绘图的频率。实时绘图非常耗时把所有绘图操作移到循环结束之后或者每50代记录一次绘图数据速度能提升好几倍。第三将终止条件改为连续无改进次数上限。如果300代内最优解都没有更新大概率已经收敛直接跳出循环。这个策略在α很小的早期版本里特别有效。5.3 动态寻优曲线抖动剧烈的原因与处理收敛曲线抖动剧烈说明算法在高温度阶段接受了大量差解这是设计使然不算bug。但如果你发现抖动在整个运行过程中一直存在甚至在低温阶段还是高幅度震荡那就是参数出了问题。一种常见情况是T0设得过大比如10000以上导致前期几百代几乎接受所有新解完全像随机游走。这时可以缩短高温阶段让T0落在ΔE分布的1到2倍范围内即可。另一种情况是扰动太激进比如2-opt逆序的长度经常等于N/2以上会让路径结构被频繁打乱难以收敛。我一般限制逆序长度不超过N/2用if diff N/2, continue; end之类的方式过滤掉过大的扰动。还有一个和动态寻优相关的技巧始终保留历史最优解。不要只维护当前解。在高温阶段当前解可能跑得很差但历史最优已经被记录下来输出时用历史最优而不是当前解。这个最简单也最容易被忽略。6. 进阶扩展与实际工程经验6.1 从TSP到带约束路径优化SA-TSP的逻辑完全可以扩展到更实际的场景比如带时间窗的车辆路径问题VRPTW、多旅行商问题MTSP。只需要修改两点一是解的表达方式比如MTSP需要在排列中插入分隔符二是目标函数把时间窗惩罚、车辆数等约束作为额外项加到总代价里。模拟退火最友好的地方就是目标函数可以写得非常灵活不需要像线性规划那样严格构造。我做过一个带时间窗的变种目标函数改为总行驶距离 惩罚系数 × 超时总量惩罚系数可以看成另一个温度控制量一开始放大惩罚后面逐渐调整权重效果相当不错。6.2 多次运行取优与随机种子管理由于SA是随机算法单次运行存在不确定性工程上很容易出现两次结果差异大的情况。我习惯的做法是写一个外层脚本连续跑10次或20次每次都换随机种子记录每次的最优距离最后取最小值、均值和方差。标准差能告诉你算法的稳定性如果标准差太大说明参数还需要调整。设置全局随机种子rng(固定值)可以让结果可复现这对调试和写文档非常有用。我在项目服务器上跑批处理时会为每次实验分配不同的种子确保结果的独立性和可对比性。6.3 给初学者的三个实操建议第一先不要急着优化性能把算法跑通的代码控制在100行以内所有参数都用明文常量方便观察。第二务必画收敛曲线和路径图肉眼观察比任何指标都直接。第三改参数时一次只改一个否则出了问题根本不知道是哪个参数引起的。我自己在调试阶段还有一个坚持很久的习惯把每次实验的参数和结果记录到一个表格里用日期命名。不要靠记忆因为alpha从0.95改成0.96效果可能差别不大但记录多了就能看出模式。比如我记录下18组参数后才确定30城市下alpha0.99是最合适的设定。最后再分享一个小技巧在算法主循环里加一个lastImproved变量记录历史最优解更新的迭代数。当循环结束还没有找到更好的解时如果lastImproved距离当前迭代已经超过500代可以直接跳出。这个小改动让我的SA-TSP在几十个测试用例上平均快了30%同时几乎不损失解质量。模拟退火这东西看着简单但把细节抠到位效果能差出好几个档次。