模拟退火算法求解TSP：Matlab实现与调参实战

"拿到一个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？因为一条回路你从哪个城市开始走、顺时针还是逆时针，本质上是同一条路。我用一个表格把组合爆炸的直观程度拉满：

所以你要是自己写一个暴力搜索去解30个城市的TSP，估计要跑到宇宙热寂。这就是组合优化里最经典的"维度诅咒"。

1.2 精确解法的天花板在哪里

很多人会问：动态规划不是能解TSP吗？确实能，Held-Karp算法可以在O(n²·2^n)时间内解决TSP，核心思想是状态压缩DP。但你把n=30代进去，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]，随机抽到i=3、j=6，反转后变成 [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 = i+1），反转其实就是交换相邻两个城市，这也是一个合法且有效的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(i+1)); 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个城市、L=300为例，我从0.9到0.995都测过：

需要说明的是，这是某一次运行的记录，带有随机性，但趋势非常规律：α在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(i+1))，不要在循环里反复用欧氏距离公式现场算。第二，缩短L到100左右，同时把α从0.98改成0.95，降低外循环数量。第三，向量化内层循环：把2-opt操作和Metropolis判断写成批量形式，一次处理一批候选解，利用Matlab的矩阵运算优势。我自己实验下来，第三点在城市数超过100时才明显起效，小规模问题犯不上。

还有一个偏门技巧：对固定数据多次运行求最优时，可以先用大α快速跑一遍得到一个"靠谱解"，再用这个解作为下一轮SA的初始解，配合低温小范围搜索做精修，效果类似于两阶段的refinement，比单纯多跑几遍SA省时间得多。

5.4 代码层面的小坑清单

我把自己写代码时踩过的坑整理成一个速查表，照着排查基本能解决九成问题：

5.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和遗传算法混合，用种群选择替代单点搜索，效果会稳定很多。不过这些都是在敲熟基础版本之后的事了。

最后分享一个做演示的小窍门：如果你要交作业或者答辩，记得给代码加上一个总开关和必要的注释，固定随机种子、把城市数和参数独立成变量，让评审老师可以一键复现你的实验。这不只是为了"规范"，更是为了让别人在面对你的结果时，能快速建立信任——这一点在工程协作中比算法本身还重要。