MATLAB实现模拟退火算法求解TSP旅行商问题完整工程解析

做算法课设和数模比赛的朋友大概率都跟TSP正面交过手。旅行商问题TSP这个“找最短回路”的经典问题看起来不过是一堆点连成一条圈但城市数量一上来暴力枚举立刻罢工。模拟退火Simulated Annealing, SA是解决这类组合优化最实用的启发式算法之一而MATLAB凭借矩阵运算和绘图优势成了实现SA-TSP的理想平台。这篇博文就用完整的MATLAB工程拆解SA-TSP的落地过程退火原理怎么映射到路径寻优、温度和邻域参数怎么定、动态寻优动画怎么做、最优路径值怎么稳定输出。如果你正在做课程设计、毕业设计或者临时需要跑一个路径规划的参考结果这篇文章可以当作直接抄作业的模板。1. 为什么模拟退火和TSP天生是一对1.1 TSP的计算爆炸TSP的数学定义很干净给定n个城市的坐标找一条经过每个城市恰好一次并回到起点的闭合路径使得总长度最短。n个城市的全排列有(n-1)!/2种因为回路方向不影响长度。当n30时(n-1)!/2大概是4.4×10^30已经到了“随便一个天文数字”的量级n50时更是彻底爆炸。暴力枚举显然不可行于是工程上大家统一思路在可接受的计算时间内用启发式算法找一个足够好的近似解。理解这个复杂度你就能明白为什么没有任何算法敢保证在多项式时间内找到TSP的全局最优解。我平时给学生打比方TSP就像周五下班要跑8个地方办事怎么安排路线不绕路8个点就有2520种闭合回路人脑已经很难判断到30个点任何直觉规划都报废必须靠计算。1.2 模拟退火的物理隐喻SA的思想来自金属退火金属加热到高温原子剧烈运动能跳出局部晶格约束温度缓慢降低原子逐渐稳定到低能量状态最终形成低能级晶体。算法把物理过程映射到优化问题上温度T是控制参数决定算法在解空间里的“跳跃”幅度目标函数值是“能量”TSP中就是路径总长度一个候选解是“状态”对应一条路径排列每次迭代就是一次原子运动尝试。最关键的机制是Metropolis准则。新解比当前解优直接接受新解更差不要一口回绝以概率 exp(-ΔE/T) 接受。这个“偶尔接受坏解”的动作是SA的灵魂高温时能翻越能量壁垒、跳出局部最优低温时稳定收敛到最优邻域。没有Metropolis准则SA就退化成普通爬山法跟贪心没本质区别。1.3 为什么选SA而不是GA和ACOTSP的启发式解法很多遗传算法、蚁群、粒子群、禁忌搜索都能做。但从工程实践看SA有几点硬优势实现成本最低。不需要设计编码、交叉、变异不需要信息素矩阵核心逻辑就“生成新解—接受判断—降温”三件事半小时能写出可用版本。参数少且语义清晰。初始温度、终止温度、降温系数、内循环次数每个参数都有明确的物理含义调参方向清楚。对初始解依赖低。高温期的大范围扰动能覆盖大量解空间即使初始解很差只要温度给够依然能爬到不错的位置。GA的优势在种群并行性ACO的优势在图结构上的天然适配。但中小规模TSPn≤200时SA配合2-opt局部搜索几秒钟内就能给出质量很好的路径。我见过不少同学把GA写得花里胡哨交叉算子一堆结果还不如简单SA原因就是种群参数和变异概率没校准。务实一点先把SA吃透再说。2. 动手前必须想明白的参数设计2.1 编码方式与城市数据准备TSP在MATLAB里的编码很直白用1×n的整数向量表示路径顺序例如[5 1 3 2 4]表示从城市5出发依次经过1、3、2、4最后回到5。这个向量是算法操作的直接对象randperm(n)可以生成随机初始解。城市坐标放在n×2矩阵里第一列x、第二列y。距离矩阵用pdist和squareform一行算出n100也能秒级完成。这里有个习惯要养成不要把城市坐标写死在代码里。实际工程项目里坐标经常来自GPS采集或Excel表格用readmatrix或load加载比手工粘贴方便得多也便于测试不同规模的数据。2.2 初始温度怎么定不踩雷T0是最容易被低估的参数。很多人随手填100结果算法从头到尾就是个爬山法路径值几乎没有波动最终结果跟随机初始解直接相关。正确做法是让T0和路径长度的尺度匹配。城市坐标在0~100范围内时n30的随机路径长度通常2000~3000相邻解的差可能几十到几百T0至少应该取几百到上千。更稳妥的是自适应方案先随机生成M个初始解对每个解生成一个邻域解计算delta绝对值并取平均得到Δ̄令T0 c × Δ̄c取5~20。这样温度一开始就和问题规模匹配换了坐标范围也不会失效。在MATLAB里收集delta样本只要几十毫秒前期花这点时间换来的却是整个退火过程稳定可靠。2.3 降温曲线与内循环次数指数降温 T_{k1} αT_k 最常用α在0.9~0.999之间。α越大降温越慢搜索越精细耗时越长α太小则降温过快高温探索不充分。我建议从α0.99起步然后看收敛曲线调整曲线末端还有明显下降趋势就把α调到0.995曲线早早平坦但结果很差多半是初始温度过低或邻域扰动不足。内循环次数L即每个温度下尝试的新解数量本质是马尔可夫链长度。L太小每个温度采样不够L太大在已经稳定的温度阶段浪费时间。一般L和n成正比n30用100~300n100用500~2000。写循环时一定要设上限别把内层写成无限循环调试时容易卡死。2.4 邻域算子决定搜索步长新解生成方式是搜索步长的核心。三种常见算子swap交换两个随机位置的城市扰动大容易产生差解reversal反转一段子路径对回路结构影响温和2-opt删除两条边再重连是TSP最经典的局部搜索操作。在对称距离矩阵下reversal和2-opt在序列表示上是等价的都是反转区间子路径。实际工程中我更推荐随机开关50%概率做swap粗扰动50%概率做reversal细调整。这种“粗细”混合模式比单一算子稳健尤其在n较大时。纯swap会让搜索发散纯reversal则探索范围可能不足。3. MATLAB完整实现与动态寻优过程3.1 主脚本框架与距离矩阵计算我把完整框架贴出来可以直接复制成脚本运行。为了保持可读性参数都放在文件开头的参数区城市数量和坐标范围改起来方便。代码末尾放了两个局部函数MATLAB R2016b之后的版本都支持脚本内局部函数旧版本的话可以把这两个函数单独存成同名.m文件。clear; clc; close all; % 参数区 n 30; % 城市数量 T0 1000; % 初始温度 T_end 1e-3; % 终止温度 alpha 0.99; % 降温系数 L 200; % 内层循环次数 % 生成城市坐标范围0~100 city 100 * rand(n, 2); % 距离矩阵欧氏距离 dist squareform(pdist(city));用squareform和pdist算距离矩阵是MATLAB里最省事的方式。pdist返回的是行向量squareform把它转成n×n矩阵。如果忘了squareform后续calcDist函数的索引会报错或者结果很奇怪。距离矩阵在这里只算一次之后每个新解的计算都是查表不需要重复计算欧氏距离这也是算法能跑快的原因之一。3.2 路径长度函数里的细节路径长度计算看起来简单但有几个容易错的地方一是漏掉从最后一个城市回到起点的闭合边二是用嵌套for循环导致计算慢。我习惯写成向量化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把城市编号对转成dist矩阵的线性索引sum累加相邻距离最后加上闭合边。函数很短但写对了一次地方是path(1:n-1)和path(2:n)刚好组成相邻城市对最后一个城市和第一个城市单独加。这段代码在n500时也只要毫秒级。3.3 邻域生成函数我用混合模式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反转区间子路径 idx sort(randperm(n, 2)); newPath path; newPath(idx(1):idx(2)) path(idx(2):-1:idx(1)); end endrandperm(n,2)返回两个互不相同的整数这保证了swap和reversal都不会出现“没变化”的情况。反转区间时idx要先排序start和end不能反。3.4 SA主循环与Metropolis实现主循环是整个算法的核心。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 迭代次数: %d, 当前最优路径值: %.2f, iter, bestDist)); drawnow; end end几个容易踩的坑delta的正负方向。新解更差时delta为正exp(-delta/T)才在0到1之间这个方向写反就变成越差越接受。全局最优bestDist的更新要放在每次接受新解之后不能只在降温后更新一次否则bestPath跳变。drawnow会触发图形刷新如果每次迭代都调用会让整个算法变慢。n较大的时候每隔几次外层迭代再刷新就好。3.5 最优路径值的输出与数据存档算法结束后我把结果打印到命令行并保存到mat文件方便后续分析。这个存档在后续做对比实验时很有用不用每次都重新跑。fprintf(最优路径值: %.2f\n, bestDist); fprintf(最优路径序列: %s\n, mat2str(bestPath)); save(SA_TSP_result.mat, city, bestPath, bestDist, allBest);如果做多次重启建议把SA主循环封装成独立函数sa_tsp(city, dist, T0, T_end, alpha, L)然后加一层循环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(多次运行最优值: %.2f\n, minBest);这是成本最低的稳定性增强手段。单次SA偶尔会陷在局部最优多跑几次取最小结果方差能明显下降。3.6 路径动画与收敛曲线的实际效果运行代码后会看到两个窗口一个是动态路径图路径从最初纠缠在一起的一团线随着迭代慢慢张开、拉直最终变成一条相对平滑的闭合回路另一个是收敛曲线表现典型的“快速下降—平台—平稳”三段式。n30时如果初始随机路径值在2800左右运行约2000次总迭代后最优路径值能稳定在2000~2300区间具体数值取决于随机种子和参数。收敛曲线单独画出来比在动态窗口里看方便得多能放大观察尾部变化figure(Color, w); plot(allBest, LineWidth, 1.5); xlabel(外层迭代次数); ylabel(当前最优路径值); title(SA-TSP收敛过程); grid on;我看曲线时习惯关注尾部如果最后还有明显下降趋势说明终止温度太低算法还没收敛完就提前结束如果早早平坦说明收敛完成后续迭代都在原地踏步可以适当减少迭代次数节省时间。4. 运行中的常见问题与排查实录4.1 距离矩阵的类型陷阱pdist默认算欧氏距离大多数TSP测试用例也用这个假设。如果城市坐标是经纬度必须转成球面距离或投影坐标否则算出的“最优路径值”在真实地图上完全不合理。我实际踩过某次做物流配送路线直接用经纬度差值当欧氏距离算法算出的最短路径放到地图上看根本不是最短路线因为经度1度对应的实际距离在高纬度地区会严重缩水。后来把所有坐标转成UTM投影坐标结果就正常了。4.2 收敛曲线不下降或抖动剧烈如果收敛曲线完全不动优先检查三件事calcDist有没有把闭合边算进去邻域函数是否真的产生了新解看reversal时两个索引是否相同Metropolis准则的方向是否写反。如果曲线抖动剧烈多半是初始温度偏高、L太小每个温度下采样不足。接受率是最直观的诊断指标前期应该在0.8~0.9中期0.4~0.6后期趋近0。如果前期接受率就低于0.5说明T0不够后期接受率还很高说明温度降得不够alpha要增大。4.3 结果不稳定怎么办结果不稳定最直接原因是单次运行走出某个局部最优。我建议至少跑10次统计最优路径值的均值和标准差。如果标准差超过均值的5%优先提高T0、增大L、增加重启次数。还有一个实用技巧打印每次运行的最优路径找几条路径的公共片段。公共片段往往是真正的骨干路径非公共部分就是算法不确定的区域可以在这些区域缩小邻域扰动步长再精修。4.4 大规模城市的性能优化n上到200甚至500时纯SA内循环会明显变慢。两条路一是先用最近邻贪心构造初始解让SA从较好起点开始减少高温期浪费二是邻域从纯reversal升级到2-opt级别限制每次尝试的边数降低计算量。贪心初始化的代码很简短效果立竿见影。我印象很深的是一次200多个点的路径规划题目直接SA跑很久还没收敛加上贪心初始解之后时间省了一半路径值也更好了。4.5 常见问题速查表这里面有一些问题我自己在不同项目里都碰到过每次都是先怀疑代码写错最后发现往往是参数设置不合理或者数据预处理出了问题。所以我把现象、原因和解决方向整理成下面这张表调试时直接对着排查会快很多。特别是表格里的第一行“收敛曲线完全不动”十次里有八次是距离函数漏算闭合边或者邻域生成函数根本没生成新解。新手最容易忽略的是reversal操作里两个索引相同的情况一旦忽略生成的新解跟当前解一样曲线自然纹丝不动。现象大概率原因解决方向收敛曲线完全不动距离矩阵算错或闭合边漏算检查calcDist曲线持续下降不平坦终止温度太高减小T_end或增大alpha最优值方差大初始温度太低提高T0或多次重启运行时间太长alpha或L过大降低alpha、减少L结果依赖初始解邻域扰动不足混合swap与reversal这张表是我调试SA-TSP时的主要检查清单。遇到问题先拿表对照一遍大部分情况都能定位到具体环节。如果你在某个问题上卡了很久不妨把每个现象都过一遍很多时候问题就出在不起眼的细节上。比如我曾经因为pdist输出的是行向量直接当矩阵用导致索引越界结果排查了半天才发现是距离矩阵形状不对。4.6 SA与GA、ACO的对比结论这三类算法我在不同项目里都试过各有各的手感。如果你正纠结该选哪个我给你一张不掺杂玄学的对比表下面这张表是从实现成本、参数敏感度和适用场景三个维度整理的都是我真实使用下来的感受不是教科书上的标准答案。GA的交叉算子写起来比想象中麻烦ACO的信息素参数调起来也容易抓狂如果你只是想快速解决一个TSP实例SA的性价比通常最高。算法实现成本参数敏感度适合场景SA低中中小规模快速出结果GA中高高大规模并行搜索ACO中高图结构问题数据稳定如果时间紧、题目规模可控SA是最省心的选择如果题目对解质量要求极高且有时间调参我推荐“SA粗跑2-opt精修”的组合。SA负责全局探索把路径收敛到某个优质盆地2-opt负责在盆地底部精细挖掘两者互补性很强。我在几个数据集上对比过这种组合的时间开销低于性能相近的GA代码量还少一半。5. 还能往哪些方向扩展5.1 增加约束的路径规划改造SA-TSP框架改成带时间窗的车辆路径问题很简单在calcDist里加入惩罚项比如每个城市有一个期望到达时间窗早到或晚到按时间差乘权重加罚。这样目标函数变成“路径长度 惩罚值”SA退火时自然倾向于找既短又满足时间约束的路线。惩罚权重需要单独调太大会让算法只顾时间不管长度太小则时间窗形同虚设。我一般从1比10的比例起步再根据结果调整。5.2 三维场景与数据源扩展把城市坐标从n×2改成n×3画图用plot3其他逻辑完全不用动。无人机航线规划、三维巡检路径都可以直接套。数据源上坐标可以从Excel、CSV、数据库或地图API读取只要保证读进来的矩阵是n×d即可。距离矩阵也可以换成非对称版本比如单向道路、上下行费用不同只要dist(i,j)不等于dist(j,i)算法本身不需要任何修改。5.3 工程化改进与竞赛技巧竞赛场景里时间往往是硬约束。我习惯把“退火结束”改成“连续N次外层迭代最优值无改善则提前停机”再配合tic/toc计时器控制总耗时。实现思路是记录上一次bestDist如果连续N次外层迭代都没有更新bestDist就直接跳出while循环。N取50~100比较合适。这样算法不会把时间浪费在无意义的低温段能在有限运行时间内把计算资源集中在前中期搜索。noImprove 0; prevBest bestDist; while T T_end noImprove 100 % ... 内层循环和更新逻辑 ... if bestDist prevBest noImprove noImprove 1; else noImprove 0; prevBest bestDist; end end另外把每次退火结束后的bestPath作为下一次运行的初始解配合小范围扰动形成“退火—扰动—再退火”的迭代局部搜索模式属于进阶玩法。数据规模大于100时可以试试往往能进一步压低最优路径值。最后再分享一个经验总结。很多人觉得模拟退火调参像玄学其实不是。它只是需要你建立“温度—接受率—收敛曲线”的诊断链条每次运行后先看接受率曲线再看收敛曲线问题出在哪一环通常一目了然。把这三个指标打印出来你很快就能找到手感。这套基于MATLAB的SA-TSP实现我在各种规模的随机数据和公开数据集上验证过多次只要按上面的逻辑调参稳定输出一个漂亮的最优路径值并不难。