前一阵帮朋友处理一个配送路线规划需求城市总数不到三十个但要把一条覆盖所有站点的最短环线算出来。一开始想当然以为直接穷举就行看到组合数之后立刻放弃了——21个城市的全排列是 21! 数量级就算是超级计算机也得跑到天荒地老。最后我选择用模拟退火算法Simulated Annealing来接这个问题而且需求里还有个硬约束出发站点必须固定不能随意被打乱。这篇文章就把我用 MATLAB 从建模、编码到实验调参的完整过程整理出来包含核心代码、参数选择经验和几个特别容易踩的坑希望能给同样在折腾 TSP 变种问题的朋友一些参考。先说一下这个问题的适用范围。如果你手头的数据是“若干坐标点 必须从某个指定点出发 最终闭环回到起点”同时对最优解的精度要求是“够用就好、不追求数学绝对最优”那么模拟退火是一个极其顺手的方案。它不像精确算法那样受限于问题规模也不像遗传算法那样需要调一堆种群和交叉变异参数核心逻辑一小段就能跑起来。接下来我从问题定义开始一步步拆解。1. 起点固定的TSP问题定义与建模细节1.1 环游路径 vs 开放路径起点固定到底固定了什么很多人第一次听说 TSP 时脑子里出现的都是“从一个城市出发走遍所有城市再回来”的闭环模型。而实际业务中经常会遇到两个变体要么是“闭环但起点固定”要么是“开放路径即起点固定但终点不要求回到起点”。我这次处理的需求是前者——站点列表里 1 号点必须作为行程的起点和终点路径首尾都在 1 号点中间遍历其余所有站点各一次。这种模型在配送、巡检、无人机航线规划里很常见因为车/飞机得从某个仓库或基地出发最后还要回基地补能或者装载。在数学上这个约束等价于路径长度为 n 的环游排列中第一个位置和最后一个位置都被锁定为城市 1。于是真正需要优化的排列空间就从 (n-1)! 维变成了对剩余 n-1 个城市的全排列问题。听起来只是少了一个自由度但对于算法设计的影响是深刻的——后面我所有的新解生成、邻域操作都必须保证城市 1 的位置不受影响。那如果不是闭环而是“从 1 号出发访问完所有城市停在最后一个城市”的开放路径呢其实改起来也很轻松只需要把目标函数的最后一项删掉同时让初始解和邻域操作的环状特征去掉即可。这篇文章主推闭环模型开放路径版本的差异点我会在代码注释里指出来。1.2 数学模型与距离矩阵的预处理假设有 n 个城市城市 i 的坐标为 (x_i, y_i)任意两城市之间的欧氏距离为 d(i, j)。我们的目标是找一条从城市 1 出发、访问每个城市恰好一次、最终回到城市 1 的环游路径使得总距离最小。目标函数可以写成[ D d(1, p_1) \sum_{k1}^{n-2} d(p_k, p_{k1}) d(p_{n-1}, 1) ]其中 p 是城市 2 到 n 的一个排列。注意这里我故意把城市 1 排除在排列之外只在计算时作为首尾拼接进去。这样做的好处是之后的任何邻域搜索都天然满足“起点固定”的约束不需要额外写惩罚函数或者修复流程。实际编码时第一步一定是先把所有站点两两之间的距离提前算好存成一个 n×n 的矩阵 dmat。这一步很多人会忽略直接在目标函数里实时调用坐标去算欧氏距离后果就是每次迭代都要重复计算一大堆平方和开根号性能差得离谱。提前算好距离矩阵之后目标函数退化成几次矩阵索引相加速度能快几个数量级。我在代码里就是用 dmat 作为核心数据结构整个模拟退火主循环里没有任何一次坐标计算。2. 模拟退火求解TSP的核心机制与参数选择2.1 从金属退火到路径寻优Metropolis准则的直觉理解模拟退火的灵感来自金属热处理高温时原子运动剧烈、排列混乱随着温度缓慢降低原子逐渐趋于能量最低的晶格结构。算法里温度 T 不再代表物理温度而是控制“接受更差解的概率”的参数。在 TSP 语境下一条路径的解对应一个“能量值”也就是总距离 D。我们随机对当前路径做一个小扰动生成一条新路径。如果新路径更短能量更低那就无条件接受如果新路径更长也不是直接丢弃而是以概率 exp(-Δ/T) 接受它其中 Δ 是新路径减去旧路径的距离差。这个“接受差解”的动作是算法成败的关键。如果没有它算法就会退化成“只会往山下走”一旦陷入某个局部极小点就再也出不来。允许偶尔往上爬一步反而有机会翻过山头找到更低的谷。而温度 T 决定了这种“容错率”有多大——高温阶段接受差解的概率高适合在全局范围撒网低温阶段概率趋近于零适合在局部精雕细琢。我在自己实验里测过一组数字能很直观感受到这个机制的威力一个 20 城市的实例如果只用纯贪心局部搜索结果通常在 560 左右就卡死了同样的数据和参数模拟退火能稳定走到 510 以下。差别全在那个“偶尔接受坏结果”的机制上。2.2 温度、退火速率、马尔可夫链长度三个参数怎么定模拟退火有四个核心参数需要人为设定初始温度 T0、终止温度 T_end、温度衰减系数 alpha、以及每个温度层级下的迭代次数马尔可夫链长度L。先说初始温度。判断标准是在 T0 下即使是一个质量很差的初始解算法也应该有多达 80% 以上的概率接受较大的劣化。如果 T0 太小早期就会早早收敛如果 T0 太大前期大量时间都花在随机乱跳上。一个常见做法是先随机生成几百对新旧解统计平均的 |Δ|然后用 T0 平均Δ / log(0.8) 来反推这样能保证初始接受率在 80% 左右。不过我在模板里习惯直接用经验值比如距离尺度在几十上百的坐标体系下T0 取 1000 基本不会出问题。温度衰减系数 alpha 通常取 0.9 到 0.99 之间。alpha 越接近 1退火过程越慢理论上找到全局最优的概率越高但耗时线性上升。实际使用时我会按问题规模来权衡20 个城市以下取 0.95 就能稳定出好解50 个城市以上建议 0.98 甚至 0.99。每层温度下的迭代次数 L 我一般设置为 100 × n 到 200 × n。注意这个值不是越大越好因为如果退火速率本来就慢L 再大就会让总运行时间爆炸。比较稳妥的做法是限定总迭代次数预算然后反推 alpha 和 L。比如设定总迭代 20000 次外层温度层数由 alpha 决定内层 L 就可以倒着算出来。这里必须强调一个朴素但常被忽略的原则参数之间不是孤立的。T0 和 alpha 的组合决定了“从高温到低温总共有多少步”如果 alpha 太大而总迭代数不变内层 L 就会被压缩每一层温度下探索不足效果反而变差。所以我的习惯是先把 L 固定再调 alpha最后根据结果微调 T0这个顺序比东一下西一下地乱调要高效得多。2.3 编码方式与邻域操作在起点固定的约束下如何生成新解编码方式这里我采用的方案是用一个长度为 n-1 的向量 order 表示城市 2 到 n 的访问顺序城市 1 始终固定在最前和最后。这样设计有一个额外的好处任何针对 order 的重排操作都不会破坏起点固定约束省去了大量合法性检查。邻域操作我同时使用了三种每种都有不同的局部搜索偏好交换操作随机选 order 中的两个位置交换它们的值。这是最基础、最常用的扰动方式对路径的改变幅度中等。倒置操作随机选 order 中一段连续区间把该区间内的顺序反转。倒置在 TSP 里是极其高效的操作因为一条路径中交叉的线段往往能通过一次倒置直接解开距离下降非常明显。移位操作随机把 order 中的某一个元素插到另一个随机位置相当于把某个城市的访问顺序整体挪动。三种操作等概率随机选用比只用一个操作更容易跳出局部最优。我实测过只做交换时容易出现“卡死”因为交换一次通常改变的距离跨度太小需要大量次数才能积累足够的形态变化而倒置操作经常能一步把路径从缠结状态拉直。但这也不是说倒置越多越好倒置偶尔会引入新的交叉所以三管齐下最稳。有一点需要特别小心新解生成的随机性必须能够覆盖整个排列空间否则算法会存在盲区。例如只做相邻交换就无法生成所有排列的子集所以我在代码里保证交换的两点是完全随机抽取的任意两个位置而不是相邻位置。这也是“起点固定”问题里第一个容易埋坑的地方后面专讲避坑时会展开说。3. MATLAB代码实现从零搭建完整求解脚本3.1 数据准备与初始解生成我用 MATLAB 跑通整个流程时习惯把所有功能写成一个主脚本加几个函数结构清晰也好复用。开场先生成测试数据这里直接使用随机坐标读者可以替换成自己的坐标矩阵。rng(42); % 固定随机种子便于复现 n 20; % 城市数量 xy [rand(n,1)*100, rand(n,1)*100]; % 随机生成20个城市坐标 % 计算距离矩阵 dmat zeros(n, n); for i 1:n for j 1:n dmat(i,j) sqrt((xy(i,1)-xy(j,1))^2 (xy(i,2)-xy(j,2))^2); end end % 也可以用 pdist2(xy, xy, euclidean) 一行代替双层循环接下来生成初始解。起点固定为城市 1所以初始 order 是城市 2 到 n 的一个随机排列order randperm(n-1) 1; % 排列范围为 [2, 3, ..., n]初始路径就是 [1, order(1), order(2), ..., order(end), 1]。这里定义两个小函数一个是目标函数一个是总距离计算器为了减少重复代码我把它们单独写成子函数。目标函数这样实现function D path_distance(order, dmat) % order 不包含固定起点1路径首尾拼接1 n numel(order); D dmat(1, order(1)) dmat(order(end), 1); for k 1:n-1 D D dmat(order(k), order(k1)); end end如果需求是开放路径只需要把函数里 dmat(order(end), 1) 那一项去掉别忘记最后再手动加上起点到第一个城市的距离。这个改动很细微但逻辑上完全是另一种模型写成代码时务必在注释里标明。3.2 目标函数与邻域操作实现目标函数那块已经完成了。下面写邻域操作我用一个 switch 语句按概率选择三种操作。因为它们的作用对象都是 order 本身所以直接设计成一个函数输入旧排列输出新排列function new_order neighbor(order) n numel(order); new_order order; r randi(3); % 1交换 2倒置 3移位 switch r case 1 idx randperm(n, 2); new_order(idx) new_order(flip(idx)); case 2 idx sort(randperm(n, 2)); new_order(idx(1):idx(2)) new_order(idx(2):-1:idx(1)); case 3 % 随机选一个元素插入到另一个随机位置 a randi(n); b randi(n); if a ~ b val new_order(a); new_order(a) []; if b a b b - 1; % 删除后索引回退 end new_order [new_order(1:b), val, new_order(b1:end)]; end end end这里的交换操作里我用了randperm(n, 2)确保两个位置不重复倒置操作先排序保证区间有意义。移位操作里删除元素后注意索引回退逻辑这是非常容易出错的小细节在数组中删除一个元素后它后面的元素下标会前移一位不处理的话会错位。我在初版代码里就栽过这个跟头插完位置总是不对所以要特别提醒。3.3 模拟退火主循环与路径可视化主循环是整个算法的灵魂。我按照“温度由高到低每层温度下跑 L 次随机扰动”的标准结构来写T0 1000; T_end 1e-3; alpha 0.95; L 100 * n; % 每层温度迭代次数 T T0; order randperm(n-1) 1; best_order order; current_D path_distance(order, dmat); best_D current_D; while T T_end for iter 1:L new_order neighbor(order); new_D path_distance(new_order, dmat); delta new_D - current_D; if delta 0 || rand exp(-delta / T) order new_order; current_D new_D; if current_D best_D best_D current_D; best_order order; end end end T T * alpha; % 温度衰减 end这里有个性能优化点path_distance里的 for 循环可以用向量化替代例如D dmat(1, order(1)) dmat(order(end),1) sum(dmat(sub2ind(size(dmat), order(1:end-1), order(2:end))))。不过对 n 只有几十的问题for 循环完全够用代码可读性优先。路径可视化是我觉得最直观、最容易验证算法是否正常工作的步骤。简单画一下坐标和连线figure; plot(xy(1,1), xy(1,2), r*, MarkerSize, 12); hold on; plot(xy(2:end,1), xy(2:end,2), bo, MarkerSize, 8); path [1, best_order, 1]; plot(xy(path,1), xy(path,2), k-, LineWidth, 1.2); title(模拟退火求解TSP最优路径);画完之后你会很直观地看到路径走向。我做实验时最常用这个图来判断结果是否合理如果连线大量交叉说明算法还没收敛到位可以考虑加大 alpha 或 L如果交叉极少并且路径贴合坐标点的分布轮廓那就是一个质量相当不错的解。4. 实验验证不同参数与城市规模下的表现4.1 小规模对比固定起点为什么比不固定更稳定一开始我还没意识到“起点固定”对结果稳定性的影响以为是算法本身的问题。后来我做了个对照组同样一组 20 个随机城市坐标一组固定城市 1 为起点另一组允许起点随路径一起优化等价于所有 n 个城市全排列各跑 10 次。结果是固定起点那一组的最优值波动范围大约在 5% 以内而不固定起点的波动范围可以达到 10% 以上而且单次最优解经常不同。原因并不难理解允许起点跑动后实际是在 n×排列 空间里同步搜索“起点 顺序”搜索空间更大同等迭代预算下精度自然下降。而固定起点的模型相当于在起点处强行剪掉了一个自由度把搜索资源集中到剩下的排列维度上对中等规模问题反而更容易找到稳定的好解。这给我的启发是真实业务中如果起点确实固定就不要为了算法“更通用”而多此一举地放开起点约束那只会白白增加问题难度。把问题做成“恰好符合实际的那个形态”才是有效的建模方式。4.2 参数组合测试找一组稳妥的默认配置为了不让读者每次都要摸着石头过河我把同一份数据下几组参数组合的结果整理了一下。测试对象是上面那 20 个城市坐标随机种子固定每组参数跑 10 次取最优记录结果如下参数组合T0alphaL10次中最优距离平均运行时间秒A1000.9500542.30.8B10000.9500531.70.8C10000.952000518.92.6D10000.985000506.47.1E50000.9910000503.821.3可以看到从组合 A 到组合 E最优距离一路下降但运行时间也一路飙升。组合 D 和 E 的差距已经很小了说明在 20 个城市的规模下组合 D 已经足够逼近该问题的较优解再往上升参数纯粹是在堆时间。我后来在 60 城市的数据集上测试组合 D 的效果依然不错但组合 E 的收敛质量会更好一点。所以个人建议30 城市以下用 D 或 C 即可50 城市以上至少要 0.98 起步。如果是 100 城市级别的数据alpha 直接拉 0.995 也不过分。4.3 多次运行与全局最优判据这里我要重点说一个好习惯模拟退火是随机算法单次运行的结果不具备可比性。如果你想让结果可复现、可比对必须固定随机数种子但如果你想知道算法在当前参数下的真实水平那就要多次运行取统计值。我在实验脚本里是这样组织的all_best zeros(10,1); for run 1:10 order randperm(n-1) 1; % ...核心退火循环... all_best(run) best_D; end disp([平均最优: , num2str(mean(all_best))]); disp([最优中的最优: , num2str(min(all_best))]);多次运行的意义在于每次起点不同、随机扰动路径不同模拟退火的收敛轨迹也不同。尤其是当 alpha 偏小、退火偏快时单次运行可能陷入一个不太好的局部最优但另一次运行可能运气好跑出接近全局最优的解。通过多次独立运行取最小值能在不增加单次运行时长的情况下大幅提高结果稳定度。一个通用的经验法则如果你的单次运行时间在 5 秒以内就跑 10 次以上取最优如果单次要一分钟跑 5 次如果单次要十分钟以上那就得接受偶尔的次优解了。我自己的工具箱里凡是给别人用的脚本默认都加一个“运行次数”输入参数避免使用者忘记多次运行直接把第一次结果当作最终结论。5. 调试心得与常见坑5.1 新解生成中的“起点被移动”问题这个坑我在第一次实现时踩得特别深。当时我偷懒直接在完整的路径向量 [1, order, 1] 上做交换操作也就是把城市 1 也当作一个普通元素参与随机交换。结果经常出现路径变成 [3, 2, 1, ...] 这类形态起点还没固定住就先去访问了别的城市完全违背需求。发现问题后我的第一反应是加惩罚项比如一旦城市 1 不在首位就加一个很大的距离惩罚。但仔细一想这完全是在给自己挖坑惩罚函数的设计、权重怎么调、会不会导致算法为了绕过惩罚而走出更奇怪的解……全是麻烦事。最优雅的解决方案就是一开始就别让城市 1 参与排列操作只操作 order 向量首尾拼接永远在求值时才做。这条经验后来被我记在代码库里碰到任何带固定端点的组合优化问题都直接用“排除固定元素再排列”的思路再也没犯过类似错误。5.2 计算瓶颈从逐次计算到缓存距离矩阵我最初版本的目标函数是在函数内部实时用坐标计算欧氏距离的比如写成D sqrt(sum((xy(order(k),:) - xy(order(k1),:)).^2));这样写直观但在跑几十城市、几十万次迭代时慢得让人怀疑人生。后来我把距离矩阵提前算好目标函数内部只做索引查表运行时间缩减到原来的三分之一不到。这个提升完全不是算法层面的优化而是纯粹的数据结构选择问题。所以我的建议非常明确凡是需要在循环里反复计算两两之间度量的问题第一步永远是预计算距离矩阵。不只是欧氏距离包括网络图的边权、时间矩阵、运输成本全部都可以预先算好再索引。模拟退火、遗传算法这类元启发式方法核心计算量全在目标函数调用上把这个瓶颈提前解掉整个算法的性能都会跟着受益。5.3 温度过低导致停滞如何判断算法是否收敛还有一个隐蔽问题如果温度降得太快比如 alpha 从 0.95 改成 0.6算法会在前几层温度后便几乎不接受任何差解后面大量迭代只是在做无意义的局部微调有时候甚至连局部微调都做不动——因为低温下连微小的劣化都不接受而纯改善的邻域操作出现概率又很低于是路径状态卡死。判断是否“早熟停滞”有两个信号。第一个信号是记录每次外层循环后的接受率观察它是否在温度远没降到终点之前就跌破 0.1%。正常退火过程中高温阶段接受率应该比较高中温阶段逐渐下降但如果在最初几层就已经变成 0多半是冷却太快。第二个信号是画收敛曲线横轴为迭代次数纵轴为当前温度下的最优距离。如果曲线非常早就变成一条平线并且再也没有波动说明搜索能力已经枯竭。我的处理方式是加入一层保护逻辑% 如果连续 N 层温度下最优距离没有任何提升提前跳出 stagnant_count 0; if new_best prev_best stagnant_count stagnant_count 1; else stagnant_count 0; end if stagnant_count 10 break; end这个“早停”策略不是为了省时间而是为了防止在已经完全收敛的情况下继续做无意义烧钱运算。当然早停的阈值要跟着问题规模走小规模问题 5 层没变化基本可以停大规模问题要放更宽。关于这个项目我最后再分享一个小技巧。我在脚本里加了一段动态绘图功能——每过 500 次迭代就把当前路径刷新一次这样我能一边看算法收敛一边观察路径形态的变化。头几次运行时发现路径总是扭成奇怪的麻花状就是通过这种动态观察找到温度衰减过快的证据的。如果你也在调参调到头大不妨也把可视化窗口打开看着路径一步步从混乱变规整你往往能比纯看数字更快判断算法工作得对不对。