★ FEATURED ARTICLE 拿到一个30个城市的TSP实例时我第一反应是大骂自己手贱——明明知道旅行商问题是个NP难问题还是忍不住想跑一遍精确解。30个城市的路径总数大约是2.65×10^32穷举一下就秒懂什么叫组合爆炸。这种情况下模拟退火算法几乎是性价比最高的入场选手。这篇文章我从TSP问题本身的难处讲起把模拟退火的核心机制拆成物理直觉、参数逻辑和Matlab代码三块最后再把我调参时踩过的坑、实测的实验记录一起整理出来。无论你是课程设计遇到算法题还是数模竞赛需要快速解决问题这篇都可以让你直接跑通一套能出图、能分析的方案。1. 先把TSP为什么难这个问题说透1.1 TSP的本质一个推销员的回家之路旅行商问题的描述非常简单一个推销员从任意一座城市出发要遍历所有n个给定的城市每个城市只去一次最后返回出发城市求一条总路程最短的闭合回路。数学上就是给定n个城市的坐标或者两两之间的距离矩阵寻找一个城市排列π (π₁, π₂, …, π_n)使得总距离D Σ d(πᵢ, πᵢ₊₁) d(π_n, π₁)最小。注意最后一项路径必须闭合这也是很多人初写代码时最容易漏掉的地方。这个问题的难点不在于描述而在于解的数量。当n不太大时所有可能的回路数是(n-1)!/2——为什么除以2因为一条回路你从哪个城市开始走、顺时针还是逆时针本质上是同一条路。我用一个表格把组合爆炸的直观程度拉满城市数 n路径总数 (n-1)!/2直观感受512手算都行10181440穷举勉强可行15约4.36×10^10计算机也很吃力20约6.08×10^16一秒算一百万个也要两千年30约2.65×10^32彻底放弃穷举所以你要是自己写一个暴力搜索去解30个城市的TSP估计要跑到宇宙热寂。这就是组合优化里最经典的维度诅咒。1.2 精确解法的天花板在哪里很多人会问动态规划不是能解TSP吗确实能Held-Karp算法可以在O(n²·2^n)时间内解决TSP核心思想是状态压缩DP。但你把n30代进去2^30 ≈ 10亿再乘以n²900运算次数接近万亿量级而且需要开一个2^n大小的状态表内存吃紧。更尴尬的是n到50、100之后这个复杂度直接崩掉。于是工程界转而依赖两类方法一类是近似算法和启发式算法比如最近邻法、2-opt局部搜索、Christofides算法速度快但结果不保证最优另一类是元启发式算法包括遗传算法、蚁群算法、模拟退火算法、禁忌搜索等。它们在有限时间内找到足够好解这件事上非常实用。我个人的经验是如果你只是想快速处理几十个城市规模的路径规划模拟退火算法是性价比之王——代码量不到一百行不用额外装工具箱调参数的空间也直观还能把收敛过程画得特别漂亮。2. 模拟退火算法从金属车间到组合优化2.1 物理灵感为什么缓慢降温能获得好晶体模拟退火算法Simulated Annealing, SA的思想源头是金属热处理工艺。金属在高温下原子运动剧烈、结构松散如果骤然冷却淬火内部会留下大量缺陷材料变脆如果让温度缓慢下降原子有足够时间重排就能形成能量更低的规则晶格材料更强韧。1983年Kirkpatrick等人在Science上发表了那篇著名的论文把这种物理过程迁移到组合优化目标函数值就相当于能量一个可行解相当于原子的一个排列状态而一个控制参数T就相当于温度。算法在高温时允许解大范围跳动随着温度一步步降低解逐渐稳定到某个低能量状态。这个迁移真正精妙的地方在于物理退火能到达全局最低能量状态是因为原子在高温下获得了克服局部势垒的能量而算法里的温度恰好赋予了较差解一丝存活机会让整个搜索过程有机会从局部最优的陷阱里爬出来。2.2 Metropolis准则以概率投靠更差的解模拟退火的核心决策逻辑叫Metropolis准则用一句话概括就是新解比当前解好无条件接受新解比当前解差扔硬币决定——但是硬币的偏置程度由温度控制。数学表达如下若ΔE E_new - E_cur 0接受新解否则接受概率为 P exp(-ΔE / T)。这里T就是当前温度。生活化理解就是你在找一条更短的出差路线老板给你提了一个新方案。如果新方案确实更短直接换如果新方案更绕你也不立刻否定而是想着先记下这个方向万一下一步能拐到更好的路上呢。温度高的时候你更愿意绕远路探索温度低了就老老实实走眼前的最优路线。Metropolis准则里的exp(-ΔE/T)很有意思。当温度很高时即使ΔE很大概率也接近1算法几乎是随机游走当温度很低时只有ΔE很小的差解才有概率被接受算法逐步收敛到精修模式。这个平滑过渡正是SA比单纯爬山法稳健的原因。2.3 落地一个SA求解器必须定义三件事纸上谈兵没用要写代码前先想清楚三件事解空间与目标函数TSP的解是一个排列目标函数是闭合回路总距离。邻域结构如何从当前解生成一个新解。TSP里最常用的是2-opt反转每次随机选一段路径把它倒序相当于把两条交叉边拆掉重连。邻域结构直接决定算法探索能力的上限2-opt简单高效对平面TSP特别友好。冷却计划包括初始温度T0、降温函数T α·T、终止温度T_end、以及每个温度下的内循环次数L也叫Markov链长度。冷却计划决定搜索节奏太急容易早熟太慢浪费时间。这三件事一头一尾决定了算法的行为风格。邻域结构是招式冷却计划是内功Metropolis准则是心法。后面所有的调参技巧本质上都是在调这三者的配合。3. Matlab实现从零写一个能出图的SA-TSP求解器3.1 数据准备城市坐标与距离矩阵先准备好测试数据。为了让结果可复现我习惯在代码开头就用rng固定随机种子这样每次运行跑出来的优化过程完全一致调参时能隔离随机因素。clear; clc; close all; rng(42); % 固定随机种子方便复现 n 30; % 城市数量 cities 100 * rand(n, 2); % 在[0,100]×[0,100]区域内随机生成城市坐标 % 计算距离矩阵 D zeros(n, n); for i 1:n for j 1:n D(i, j) sqrt(sum((cities(i,:) - cities(j,:)).^2)); end end如果你想用更Matlab化的写法距离矩阵可以直接一行搞定D squareform(pdist(cities));但注意squareform返回的是压缩后的行向量要用就得配合squareform完整函数处理新手容易晕所以我后面代码里都用双层循环版的D虽然看着啰嗦但语义清晰。3.2 路径表示与2-opt邻域操作TSP的路径我用一个1×n的行向量保存比如path [3 7 1 4 2 5 6]含义是从城市3出发依次经过7、1、4、2、5、6再回到3。用randperm(n)可以一步生成随机排列这就是当前解的起点。2-opt操作是整个算法里最核心的小动作。它的物理意义很直观当前路径上有两条边是交叉的或者至少是不顺路的你选出两个位置i和j把中间的城市顺序完全反转这样原来的两条边被拆掉换成两条新边路径总长度通常会有明显下降。我来看一个例子假设有一条路径 [1 2 3 4 5 6 7 8]随机抽到i3、j6反转后变成 [1 2 6 5 4 3 7 8]。注意我们不是交换两个城市而是把一整段倒序这样能在一瞬间消除多个路径交叉。% 2-opt邻域操作随机选两段反转 i randi([1, n]); j randi([1, n]); if i j [i, j] deal(j, i); end if i j continue; % 选的同一个位置没有意义跳过 end new_path cur_path; new_path(i:j) cur_path(j:-1:i); % 反转区间这里有两个坑我必须提醒第一i j之后要交换否则切片new_path(i:j)会变成空集操作第二如果i和j之间只隔一个城市j i1反转其实就是交换相邻两个城市这也是一个合法且有效的2-opt操作不用额外剔除。但如果i j反转后路径完全不变既是浪费时间又会让后面的Metropolis判断变得无意义所以我选择continue跳过。3.3 主循环外层降温内层搜索接下来是完整的模拟退火主流程。我把目标函数也单独写成一个函数这样代码结构更干净。% 模拟退火参数 T0 224; % 初始温度后面讲怎么估计的 T_end 1e-3; % 终止温度 alpha 0.98; % 降温系数 L 300; % 每个温度下的内循环次数 T T0; % 初始化解 cur_path randperm(n); cur_dist total_dist(cur_path, D); best_path cur_path; best_dist cur_dist; % 记录收敛过程 history []; while T T_end for k 1:L new_path cur_path; i randi([1, n]); j randi([1, n]); if i j, [i, j] deal(j, i); end if i j, continue; end new_path(i:j) cur_path(j:-1:i); new_dist total_dist(new_path, D); delta new_dist - cur_dist; if delta 0 || exp(-delta / T) rand cur_path new_path; cur_dist new_dist; if cur_dist best_dist best_dist cur_dist; best_path cur_path; end end end T T * alpha; history [history; best_dist]; end % 绘图 figure; plot(cities(best_path([1:end 1]), 1), cities(best_path([1:end 1]), 2), b-o, LineWidth, 1.2); hold on; plot(cities(:, 1), cities(:, 2), ro, MarkerSize, 8, LineWidth, 1.5); title(sprintf(SA-TSP 最优路径, 总距离 %.2f, best_dist)); xlabel(X坐标); ylabel(Y坐标); axis equal; grid on; figure; semilogy(history, LineWidth, 1.2); xlabel(降温次数); ylabel(当前最优距离); title(收敛曲线);这里可能有人问为什么接受条件写成exp(-delta / T) rand而不是rand exp(-delta / T)两者完全等价纯属个人习惯。真正需要注意的是当delta为正且delta/T很大的时候exp函数的结果在Matlab中会直接下溢为0此时exp(...) rand恒为假正好不会接受差解逻辑上没毛病。不过如果你对数值稳定性有执念可以加一层保护if delta 0 || (delta 50*T exp(-delta / T) rand)这个写法避免了大数指数计算的浪费不过Matlab里exp参数小于-700左右才会下溢一般场景用不上这么保守。3.4 计算总距离的辅助函数上面代码里调用的total_dist是一个局部函数Matlab脚本里定义局部函数必须放在脚本文件末尾这是新手最容易踩的坑——把函数写在了脚本开头然后报了一堆“Function definitions in a script must appear at the end of the file”之类的错误。function d total_dist(path, D) % path: 1×n 城市排列 % D: n×n 距离矩阵 n length(path); d 0; for i 1:n-1 d d D(path(i), path(i1)); end d d D(path(n), path(1)); % 闭合回路 end我一开始写这段代码时漏掉了最后一行结果算法跑出来说是300多实际路径苦不堪言直到把迭代过程画出来才发现回路根本没闭合推销员飞回了起点。如果你用向量化写法total_dist也可以写成d sum(D(sub2ind([n n], path, [path(2:end) path(1)])))但为了可读性循环版本更适合作业和演示场景。4. 参数调优别让算法玄学化4.1 初始温度到底取多少初始温度决定了算法刚开始时的胆量。温度太高前期几百步全是纯随机乱跳温度太低开局就陷入局部最优。一个非常实用的估计方法是随机采样一定数量的2-opt邻域操作统计ΔE的平均值再用期望初始接受率反推T0。假设我随机做了100次2-opt得到平均ΔE ΔE_avg 110具体数值取决于实例我希望初始接受率P0 0.8也就是80%的差解会被接受那么利用Metropolis公式P0 exp(-ΔE_avg / T0)反解得 T0 -ΔE_avg / ln(P0)。把数字代进去T0 -110 / ln(0.8) 110 / 0.223 ≈ 493。当然这个期望接受率本身是个近似实践中你可以在0.7到0.95之间选。我文章里的示例代码用了224对应的P0大约0.61对30个城市已经够用了你可以根据实际效果调整。写一个小脚本就能估出这个数不必每次靠猜deltas zeros(1, 100); for s 1:100 tmp_path cur_path; ii randi(n); jj randi(n); if ii jj, [ii,jj] deal(jj,ii); end if ii jj, continue; end tmp_path(ii:jj) tmp_path(jj:-1:ii); deltas(s) total_dist(tmp_path, D) - cur_dist; end T0 -mean(deltas(deltas 0)) / log(0.8);注意只统计正ΔE负的不需要接受概率。4.2 降温速率α的秘密降温系数α是0到1之间的数每次外循环结束时把温度乘上α。α越接近1降温越慢搜索越充分但运行时间线性增加α太小如0.9温度断崖式下跌算法很快就变成纯爬山典型症状是结果很差、路径上有明显的交叉。从网格化的视角看温度是视野半径初期大视野找方向后期小视野精修。α等于每次缩小视野的比例。以30个城市、L300为例我从0.9到0.995都测过降温系数 α最优距离单次运行外循环次数相对效果0.90483.294差早熟明显0.95415.6181一般仍有交叉0.98398.7377好路径干净0.995395.31077略好但耗时加倍需要说明的是这是某一次运行的记录带有随机性但趋势非常规律α在0.95以上时结果明显改善0.98附近已经足够好再往上性价比就不划算了。用0.98相当于在100×100的坐标范围里跑了几百次外循环足够让一个30城市的实例收敛到很接近最优的值。4.3 内循环长度L取多少合适内循环L决定每个温度下要尝试多少个候选解。L太小当前温度下探索不足温度就降下去了L太大后面几百度几乎在做重复功。经验法则L可以取城市数的5到20倍也就是30个城市用150到600。我在示例里用300跑完大概半分钟左右节奏刚好。有一个很容易被忽略的细节随着温度降低接受率越来越低内循环里大量时间在拒绝新解。如果你真的想省时间可以加一个提前退出机制——如果连续若干次内循环都没有接受任何新解就直接跳去降温。不过这样做会让算法的Markov链性质变弱作业演示可以严谨研究就算了。4.4 实验记录与收敛曲线判读我固定随机种子rng(42)跑完后初始随机路径长度大概是1480左右最终最优距离落在398到410这个区间。收敛曲线呈现非常典型的快速下降长尾精修形态前50次降温就把距离从1480压到500以内后面三百多次只是把400往395磨。这其实正是模拟退火的特性——大尺度路径调整必须靠高温期完成后期低温只是局部打磨。如果你看到收敛曲线后期还有明显的阶梯状下跌说明当前邻域结构还能发现明显改进空间可以适当提高α或者加长L如果曲线在某个值上彻底平掉再也降不动多半是卡在局部最优了可以考虑重退火温度跳回较高值重新搜索。5. 常见问题与排查技巧实录5.1 每次运行结果都不同而且差距很大模拟退火算法本质是随机算法每次运行结果不同是正常的。但如果你发现两次最优距离能差出30%以上说明参数有问题。排查思路如下先固定随机种子确认代码本身没有隐藏问题打印每个降温阶段的最优距离看算法是不是在前几十次降温就锁死了一个坏解如果锁死把初始温度加倍或者让α更接近1。我自己的标准流程是同一个参数跑10次记录最优、最差、中位数。如果10次最优距离的极差超过5%就调大α或者L如果10次结果都很接近说明参数稳健可以放心用。5.2 路径图上有明显交叉结果却不再下降路径交叉是2-opt这个邻域操作理论上可以消除的但如果算法还留有大量交叉就停止改进几乎可以断定温度场出了问题。最常见的原因是初始温度偏低让算法过早就失去了接受差解的能力其次是降温系数太小温度掉太快。这里有一个我经常用的排查技巧把每次接受差解的概率打点画出来。如果接受概率在循环进行到10%时就掉到0.01以下那后面90%的迭代都是在做无效爬山。平衡状态应该是高温期接受概率0.6~0.8中期逐步过渡到0.1末期趋近0。5.3 运行时间太长怎么优化优化SA的运行时间有几个立竿见影的思路。第一把total_dist里的距离计算用查表代替——距离矩阵D提前算好每次只查D(path(i), path(i1))不要在循环里反复用欧氏距离公式现场算。第二缩短L到100左右同时把α从0.98改成0.95降低外循环数量。第三向量化内层循环把2-opt操作和Metropolis判断写成批量形式一次处理一批候选解利用Matlab的矩阵运算优势。我自己实验下来第三点在城市数超过100时才明显起效小规模问题犯不上。还有一个偏门技巧对固定数据多次运行求最优时可以先用大α快速跑一遍得到一个靠谱解再用这个解作为下一轮SA的初始解配合低温小范围搜索做精修效果类似于两阶段的refinement比单纯多跑几遍SA省时间得多。5.4 代码层面的小坑清单我把自己写代码时踩过的坑整理成一个速查表照着排查基本能解决九成问题症状原因修复脚本报函数定义位置错误局部函数写在脚本开头把函数移动到脚本末尾2-opt操作后路径不变i与j相等切片反转无效果添加if i j continuei j时切片为空随机生成的两个位置未排序先交换让i j距离计算偏小忘了闭合回路项D(path(n), path(1))补上最后一段exp结果为NaNdelta和T计算中出现非数值检查路径索引是否越界结果不收敛初始温度太低或α太小用4.1的公式重新估T05.5 一个关于2-opt方向的进阶提醒2-opt反转只是最简单的邻域操作。实际操作中我发现对于城市数超过50的实例纯2-opt的邻域集合太大了随机采样的效率会下降此时可以考虑2-opt的快速变体只随机选一个位置i然后在距离矩阵中找一个最近的未访问邻接点构建新的边。当然这就是另一个话题了。对TSP入门和课程设计来说老老实实用2-opt反转配合均匀采样完全够用。还有一个容易被忽略的细节2-opt最开始的灵感是消除路径交叉但随机反转一段路径可能不一定会减少交叉甚至会增加交叉——这没关系因为Metropolis准则和温度机制会处理掉这种后退。不要因为看到某一步路径变差就手动干预要相信整个退火过程的统计力量。6. 写在最后的个人体会模拟退火这个算法我前前后后用过它解决TSP、排课调度、参数拟合和基站选址问题。说句实在话它很少是最优算法但它几乎总是一个够用算法。尤其是TSP这种解空间离散、目标函数计算便宜的问题SA的性价比极高一下午就能从零写完代码不需要额外工具箱调参空间直观画出来的收敛曲线还很漂亮。调参经验上我最大的收获是初始温度不是玄学而是有物理意义的——它决定你对差解的容忍度降温系数不是越大越好而是要和内循环长度配合着看。理解这两点之后你就不再是瞎试参数而是在控制一个模拟的冶金过程。如果你有兴趣继续往下扩展建议试试三个方向一是把2-opt升级成3-opt甚至Lin-Kernighan邻域结构更强了解的质量会上一个台阶二是加入重退火策略每过一段时间把温度拉回一个中间值避开长期陷入同一片局部最优三是把SA和遗传算法混合用种群选择替代单点搜索效果会稳定很多。不过这些都是在敲熟基础版本之后的事了。最后分享一个做演示的小窍门如果你要交作业或者答辩记得给代码加上一个总开关和必要的注释固定随机种子、把城市数和参数独立成变量让评审老师可以一键复现你的实验。这不只是为了规范更是为了让别人在面对你的结果时能快速建立信任——这一点在工程协作中比算法本身还重要。
阅读完成 · 觉得有帮助？
