如果你接触过组合优化一定绕不开旅行商问题TSP给定若干城市坐标找一条访问所有城市恰好一次并返回起点的最短回路。问题本身通俗易懂但它背后对应的是物流配送、工厂排程、路径规划里大量真实的排序难题。我经常在项目里用模拟退火算法Simulated Annealing快速拿到一个接近最优的路线方案Matlab写起来尤其方便没有繁琐的编译链画图也顺手。这篇文章就把我从零实现模拟退火求解TSP的完整过程、参数选择经验和踩坑记录整理出来给打算在课程设计、工程实践中用这个组合的人一份可直接复制的参考。1. TSP问题、模拟退火核心思想与算法选型1.1 TSP问题的本质与工程价值TSP的数学定义不复杂给定n个城市的坐标每个城市访问且只访问一次最后回到出发城市要求总路程最短。它之所以被称为硬骨头是因为可行解的数量实在太多了。对于一个n个城市的TSP闭合回路的总数是(n-1)!/2。n等于20的时候大约是1.2×10^16种排列n到50时这个数字已经是3×10^62量级比宇宙中基本粒子的估计数量还要大好几个数量级。这种随规模呈阶乘爆炸的现象就是组合优化里典型的NP-hard问题说白了就是规模一大最优解只能靠巧劲找不能靠蛮力争。TSP从来不止是教科书里的习题。物流公司的车辆路径规划、外卖平台的配送顺序调度、PCB电路板的钻孔顺序优化、芯片布线中的引脚连接顺序、无人机航点巡检路线底层建模几乎都能落成TSP或它的变体。我接触过的真实项目里最常见的情况是不需要数学上绝对最优的路线只需要在有限时间内给出一个比人工经验好得多的方案这种场景就是启发式算法的主场。1.2 模拟退火的物理灵感与核心步骤模拟退火算法借用了金属热处理里的概念金属被加热到高温后内部原子的排列变得混乱随后缓慢降温原子会逐渐排列成能量较低的规则晶体结构如果降温太快就会形成存在大量缺陷的非晶体。算法的创始人把这种物理过程抽象成了优化方法用一个虚拟温度控制搜索过程。映射关系非常直观算法的每一个解对应一个状态目标函数值对应系统能量温度高时算法允许大范围漫游愿意接受质量差的解温度低时算法逐渐变得保守只接受改进。核心步骤其实就是四步从当前解出发通过邻域操作生成一个新解计算新解与当前解的目标函数差值按照Metropolis准则决定是否接受新解按温度调度策略降低温度并重复这个过程。用爬山来做类比可能更好懂。传统贪心算法的特点是只往高处爬结果很容易被困在一个矮山头模拟退火则允许你偶尔朝低处走几步相当于给自己一点动能越到旁边的山坳然后才有机会登上更高的山峰。温度就是这股动能的大小温度越高你越愿意走下坡路。1.3 为什么选模拟退火而不是其他算法同样能解TSP的算法不少贪心算法、遗传算法、禁忌搜索、蚁群算法都有各自的拥趸。我在实际选择时主要看三条实现成本、参数数量、对问题形态的敏感度。贪心算法最快每步选最近的城市就行但结果往往很差尤其在城市分布不规则时很容易被短视的选择带进死胡同。遗传算法效果不错但要设计编码方式、选择策略、交叉算子、变异算子深水区都在算子设计上调试成本高。禁忌搜索需要一个记忆表防止走回头路结构上比模拟退火复杂一点而且对邻域操作和禁忌长度的设定比较敏感。相比之下模拟退火的吸引力很直接框架固定、参数少、对目标函数几乎没有要求。既不需要目标函数可导也不需要凸性假设你把一个实际问题的约束写成目标函数塞进去就能跑。当你需要在短时间内搭出一个比人工强、比精确算法快、比简单贪心稳的方案时它几乎是万金油。当然它的代价也很明显随机性强、不保证最优解、需要一定的参数调校经验后面我会专门聊怎么调。2. 算法设计与关键参数2.1 解的表达、初始解与距离矩阵TSP的编码不复杂我用整数排列来表示路径例如[1 5 3 2 4]表示从城市1出发依次访问5、3、2、4最后回到1。在Matlab里一条路径就是一个1×n的行向量操作起来很方便。初始解怎么给最省事的是randperm(n)直接随机打乱简单但随机性太大。我喜欢用随机排列配合一次最近邻修正先从城市1出发每次选距离最近且未访问过的城市依次走一圈得到一个比纯随机好得多的起点。这样做的意图是让算法一开始就站在一个相对不坏的解上把宝贵的迭代预算用在优化而不是从垃圾堆里爬出来。不过也要注意别把起点做得太好否则高温阶段的随机扰动会把解推得很远反而浪费前期探索。距离矩阵建议提前一次性算好存成一个n×n矩阵。因为整个算法过程中路径代价函数会被调用成千上万次每一次都临时用坐标算欧式距离非常浪费。如果是经纬度坐标记得用球面距离公式而不是直接把经纬度当平面坐标算欧式距离那是不少地理类项目的常见错误。2.2 邻域结构与2-opt为什么有效邻域操作决定了搜索怎么从一个解移动到另一个解这是所有启发式算法的灵魂。TSP里常见的邻域操作有三类交换算子、插入算子、反转算子。交换算子随机挑两个城市换位置简单粗暴但破坏性大插入算子把一个城市抽出来插到别的位置比较温和反转算子则是把路径中一段城市倒序排列这就是TSP里大名鼎鼎的2-opt操作。2-opt之所以效果好核心在于它能消除路径交叉。两条线段如果出现了交叉交换交点两侧的访问顺序就能把路径拉直总长度必然下降。更重要的是2-opt操作一次只改动两个连接关系对路径整体结构的影响是局部的、可控的这让它特别适合在模拟退火这种需要大量试错的过程中反复使用。相比之下swap会同时切断四条连接带来的变化更剧烈容易让已经比较规整的路径变得面目全非。我的经验是主用2-opt偶尔以10%~20%的概率插入一次swap来增加多样性。只用2-opt的缺点是搜索可能陷在某个结构模式的局部区域里混一点swap可以跳出这种模式僵局。两者结合之后模拟退火在标准数据集上的表现会明显好过单用任何一种。2.3 温度调度初始温度、降温系数与终止条件温度调度是整个算法的控制器三个问题必须回答从多高的温度开始降温多快什么时候停初始温度不能拍脑袋。一个可靠的经验方法是先让算法在当前解附近随机做几百次邻域操作记录目标函数变化的绝对值取平均值记为ΔE_avg然后设定一个理想的初始接受概率P0通常取0.8到0.95用公式T0 -ΔE_avg / ln(P0)反推。这个公式直接来自Metropolis接受概率的变换确保一开始算法有七八成以上的概率接受差解足够在解空间里撒开网。我在代码里通常把这段自动估计逻辑封装成一个子函数省得每次换问题都手动调。降温曲线最常用的是等比降温T_{k1} α·T_k。α取0.90到0.99之间越大搜索越精细但耗时越长。我自己的默认值是0.99因为温度下降慢一点后期精细搜索的收益通常能把多花的时间赚回来。终止条件我一般设两个温度低于某个下限T_end比如1e-3就停或者连续多个温度周期最优解都没有变化就提前退出省点时间。2.4 Metropolis准则与随机接受逻辑模拟退火区别于爬山法的关键就在Metropolis准则新解比当前解好无条件接受新解比当前解差以概率exp(-ΔE/T)接受。这里的ΔE是新解目标函数值减去当前解目标函数值T是当前温度。温度高指数接近1差解几乎都会进温度低指数趋近于0差解基本进不来。实现起来就是一行逻辑判断if delta 0 || rand exp(-delta / T) route new_route; end很多初学者不理解为什么要接受差解这里说透如果我们永远只接受更好的解算法本质上就是爬山法一旦站在局部最优的山头上就再也下不来了。接受差解相当于允许算法走下坡路有了这个机制才有机会绕开局部最优去探索别的高峰。有一点需要留意当温度很低时exp(-delta/T)可能因为浮点数精度下溢直接变成0这是正常的不代表出错。如果你观察收敛曲线会发现低温阶段接受率趋近于0算法基本只做改进型搜索这恰恰是设计好的降温过程快结束时应有的表现。3. Matlab实现与代码逐段解析3.1 主程序框架我用一个主函数封装整个流程外部只需要传入城市坐标和参数结构体。这版代码我刻意保持了能跑与可读之间的平衡适合直接抄去改。function [best_route, best_cost] sa_tsp(city, opts) % 模拟退火求解TSP % city: n×2矩阵每行一个城市坐标 % opts: struct可选字段 % opts.T0 初始温度留空则自动估计 % opts.T_end 终止温度 % opts.alpha 降温系数 % opts.inner_iter 每个温度下的内循环次数 % 输出 % best_route 最优路径城市序号序列 % best_cost 最优路径总距离 n size(city, 1); Dist squareform(pdist(city)); % 预计算距离矩阵 route randperm(n); % 初始解简单起见用随机排列 best_route route; best_cost path_cost(route, Dist); if isempty(opts.T0) T0 auto_initial_temperature(route, Dist, 0.9); else T0 opts.T0; end T T0; T_end opts.T_end; alpha opts.alpha; inner_iter opts.inner_iter; while T T_end for k 1:inner_iter new_route neighbor(route); % 生成邻域解 new_cost path_cost(new_route, Dist); current_cost path_cost(route, Dist); delta new_cost - current_cost; if delta 0 || rand exp(-delta / T) route new_route; % 更新当前解 end % 始终保存全局最优 if path_cost(route, Dist) best_cost best_route route; best_cost path_cost(route, Dist); end end T alpha * T; % 降温 end end主循环就是两层外层是温度循环内层是在同一温度下的探索次数。内循环里每次都生成一个新解按Metropolis准则决定是否接受。注意我每次更新完当前解后都会重新算一次目标值去更新全局最优这个细节非常关键后面会专门讲为什么不能省略。3.2 代价函数与2-opt邻域实现代价函数本质上就是累加路径上相邻城市的距离并加上最后一段回起点的距离。顺着[route, route(1)]这种拼接方式写可以避免循环遍历function c path_cost(route, Dist) idx [route, route(1)]; c sum(Dist(sub2ind(size(Dist), idx(1:end-1), idx(2:end)))); end2-opt邻域操作的实现更微妙。我的写法是随机选两个位置i和j把i到j这一段反转function new_route neighbor(route) n length(route); i randi(n - 1); j i randi(n - i); new_route route; new_route(i:j) route(j:-1:i); end这段代码看起来干净但有几个边界细节必须说明。第一个是i不能取n否则j只能等于n反转一段只有一个元素等于没操作。第二个是反转段的下标闭区间问题j:-1:i生成的是从j递减到i的序列赋值给i:j是安全的。第三个是route(j:-1:i)这种写法在Matlab里天然支持不需要额外循环。我在实际项目里会把起点固定为城市1也就是初始化时先保证route(1)1并且把2-opt的随机范围限制在i≥2。这样做的原因纯粹是减少对称重复解闭合回路绕一圈回到同一个起点起点不同含义相同固定起点等于把搜索空间缩小为原来的n分之一收敛速度有明显提升。3.3 初始温度自动估计与调用示例自动估计初始温度的子函数是这么写的function T0 auto_initial_temperature(route, Dist, P0) deltas zeros(1, 300); for i 1:300 new_route neighbor(route); deltas(i) abs(path_cost(new_route, Dist) - path_cost(route, Dist)); end avg_delta mean(deltas); T0 -avg_delta / log(P0); end它先随机生成300个邻域解统计目标函数变化绝对值的平均值然后反推初始温度使初始接受概率约等于P0。每次运行的结果可能略有差异但因为取的是平均值差异不大。如果你的目标函数值本身量纲很大比如路径长度动辄上万别用固定T0这个自动估计函数会帮你省掉很多试凑的功夫。调用和可视化也很简单city 100 * rand(30, 2); opts.T0 []; opts.T_end 1e-3; opts.alpha 0.99; opts.inner_iter 300; [best_route, best_cost] sa_tsp(city, opts); figure; plot(city(best_route, 1), city(best_route, 2), o-); title(sprintf(SA最优路径总长度%.2f, best_cost));用随机生成的城市跑一下30个城市通常几秒内就能收敛到一条相当规整的路径。如果你把opts.T0留空主函数会自动去估计温度这也是一个更省心的用法。3.4 工程化细节记录最优解、耗时控制与可视化第一个工程化细节就是全局最优解必须单独保存。模拟退火在高温阶段会大量接受差解这意味着最终返回的当前解完全可能不是整个迭代过程里出现过的最优解。我见过不少初学者的程序循环结束时输出的是当前路径结果每次运行的结果波动很大。解决办法就是像上面代码那样每次接受新解后都检查一次是否更新全局最优最后返回的是best_route而不是route。第二个细节是耗时控制。Matlab里最实用的办法是用tic和toc把总运行时间打印出来再根据时间预算去调整inner_iter和alpha。比如你只有5秒预算那就把inner_iter减半或者让alpha从0.99降到0.98。不要把时间控死因为算法本身是随机的给它的预算越充裕结果越稳定。第三个细节是路径可视化。不要只在控制台打印一个总长度数字画出路径连线图能让你一眼看出算法是否生成了合理路线。如果画出来的路径有大量交叉说明搜索还没收敛需要考虑增大alpha或内循环次数如果路径已经很顺滑但长度没有下降趋势说明参数已经足够再跑下去边际收益不大了。4. 实验结果与参数敏感性分析4.1 标准测试集上的表现要验证算法写没写对最简单是拿TSPLIB的标准数据集跑一下这些数据集有官方已知最优解能直接计算误差。我用得比较多的是att48和berlin52两个数据集。att48是48个城市的实例已知最优解是33522。用上面这套模拟退火实现配合α0.995、T_end1e-3、inner_iter100*n的参数组合多次随机运行下来最好结果通常能落到33700到34500区间大概比最优解高1%到3%。带上2-opt局部搜索精修之后还有机会压到33400到33600的范围。berlin52是52个城市的实例已知最优解是7542我的多次实验里SA经常能跑到7545到7600运气好的时候能恰好碰到7542。这个结果表明SA的定位很清楚它不会给你最优解但能在几秒到几十秒内给你一个工程上完全可用的解。如果你需要绝对最优n在50以内可以试试分支定界或者动态规划但城市数量一旦过百这些精确方法基本就别想了这时候SA和它的亲戚们才是真正的主力。4.2 参数敏感性哪些参数最关键我调参时最关心的就四个参数初始温度、降温系数、内循环次数、终止温度。它们之间的关系用一个表格总结参数推荐范围影响说明初始温度T0自动估计或取P00.8~0.95反推太低导致过早收敛太高导致前期大量时间浪费在漫游上降温系数α0.90~0.99越接近1搜索越精细耗时越长越小收敛越快但容易错过最优区域内循环次数inner_iter100~1000或100*n太小每个温度下探索不足太大运行时间线性增长终止温度T_end1e-2~1e-4越低最终结果越稳定但低温阶段收益递减实际调的时候我的策略是倒着来先固定α0.99、inner_iter100*n跑一次观察收敛曲线。如果结果不满意优先加大inner_iter而不是把α改成0.999因为后者会把运行时间拉得很夸张。如果发现曲线在温度还很高时就已经平了说明前期探索不足这时应该调高初始温度或增加内循环让算法有更多机会跳出局部。4.3 多起点重启与稳定性评估模拟退火是随机算法单次运行的结果不能代表它的真实水平。我自己的标准流程是这样的先用rng(seed)固定随机数种子跑一次确认代码可复现然后更换多个种子跑10到20次记录最优、平均和最差结果最后用matlab的parfor把这些独立运行分散到多核上并行跑十几秒就能得到一组统计结果。多起点重启在很多场景下比单次拉长搜索时间更划算。假设你有20秒预算与其让算法一次性跑20秒不如让5个独立进程各跑4秒最后取最优。这样等于同时探索了5条不同的搜索轨迹有效降低了单次随机性带来的方差。5. 常见问题与排查技巧实录5.1 常见问题速查表我把平时被问得最多的现象、可能原因、解决手段整理成了速查表现象可能原因解决方向结果始终很差接近随机路线初始温度太低或降温太快用自动温度估计α调大到0.99以上前期收敛快后期纹丝不动温度降得太快搜索结构僵化增大α混入swap邻域操作每次运行结果波动极大内循环太少随机性主导增加inner_iter做多起点重启路径出现重复或乱序编码或索引逻辑错误打印route逐步检查缩小测试规模运行时间完全不可接受α太接近1或inner_iter过大调小α用时间预算控制循环最终解明显不如过程中某一步没有单独保存全局最优像上文代码一样维护best_route5.2 三个让我印象深刻的坑第一个坑是我自己刚写SA时踩的只维护当前解不维护全局最优。我在前面反复提到过高温阶段接受差解是算法能跳出局部最优的根本原因但这同时意味着当前解随时可能变得很差。如果最后直接输出当前解你看到的可能是一个被高温阶段搅乱的中等解连算法真实能力的一半都没体现出来。正确的做法是任何位置更新解之后都顺手比较一下历史最优。第二个坑是2-opt的边界下标。我早期用randi(n)去取i和j没有限制i的最大值结果经常出现单元素反转导致算法空转。更隐蔽的情况是j取到n时反转段横跨了路径的起点虽然在数学上仍然合法但如果你的代码里用了固定起点就破坏了起点约定后续所有计算都会错位。后来我把起点固定为城市1并强制i从2开始取边界问题一次性全解决了。第三个坑是提前退出的条件设置。我一开始写过一个逻辑当前温度下连续200次迭代都没有接受新解就结束。从效果看这个规则频繁让算法在高温阶段就提前退出后期低温段的精修完全没跑到。解决办法是把这个提前退出条件只作用于当前最优解连续多轮没有任何改进而不是新解都没被接受。因为你没有接受新解不代表最优解不会在下一轮被更新不能用一个太敏感的规则打断整个降温过程。5.3 提升效果的进阶技巧如果基础版本已经通了想进一步提升解质量我推荐三个方向。第一个是退火局部搜索的组合。在每轮温度结束后对当前最优解做一轮完整的2-opt局部搜索不断尝试反转路径中的每一段直到没有改进为止。这个操作相当于把模拟退火的全局搜索和2-opt的局部精修结合起来很多TSP实例在这个组合下能达到与最优解接近的结果。实现也不难写一个while循环遍历所有ij对遇到改进就执行反转直到一轮遍历无改进才停。第二个是多起点重启。我之前说过模拟退火一定不要只跑一次。用parfor并行跑10个独立流程每个流程用不同的随机种子结束后比较所有best_cost取最小。这个方法在Matlab里改动成本极低却能让结果稳定性上一个台阶。第三个是自适应降温系数。如果你观察收敛过程发现前期几乎没有下降、后期却一直在改善说明温度下降得太快反过来如果前期下降明显、后期长时间不改善说明温度下降太慢。可以根据连续几轮最优解的变化率动态调整α改善明显时让α小一点加速收敛改善停滞时让α大一点给足探索时间。这个技巧不是必需的但对追求极致解的场合很有用。最后再分享一个我自己的习惯任何启发式算法的结果都要在固定随机种子下做可复现实验并且把最优、平均、最差三个指标一起记录。算法没有绝对的好坏只有在你给定的时间预算和解质量要求下是否够用。实践里我很少迷信某一个算法能一招通吃但模拟退火确实是我工具箱里最常被拿出来用的那一个它够简单、够灵活、够稳。很多时候工程问题的关键不是算法多花哨而是你能不能在一个下午把它调通并给出可信的结论。