前一阵帮朋友处理一个配送路线规划需求，城市总数不到三十个，但要把一条覆盖所有站点的最短环线算出来。一开始想当然以为直接穷举就行，看到组合数之后立刻放弃了——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_{k=1}^{n-2} d(p_k, p_{k+1}) + 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(k+1)); 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(b+1: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 次取最优记录，结果如下：

可以看到，从组合 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(k+1),:)).^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 次迭代就把当前路径刷新一次，这样我能一边看算法收敛一边观察路径形态的变化。头几次运行时发现路径总是扭成奇怪的麻花状，就是通过这种动态观察找到温度衰减过快的证据的。如果你也在调参调到头大，不妨也把可视化窗口打开，看着路径一步步从混乱变规整，你往往能比纯看数字更快判断算法工作得对不对。