★ FEATURED ARTICLE 模拟退火是我最早接触的一类元启发式算法当年在 MATLAB 里用它求解 TSP旅行商问题时从照着别人的代码改、到自己把整个流程吃透中间踩的坑确实不少。这篇文章我想把模拟退火的核心思路、TSP 建模方式、一套可以直接跑的 MATLAB 代码以及调参和排查故障的经验一次讲清楚。适合谁看呢刚学智能优化算法但被公式劝退的同学需要交课程设计又不知道从哪下手的初学者还有想在项目里快速验证组合优化思路的工程师——我觉得都值得花十分钟读完。1. 先搞懂模拟退火在干什么1.1 从金属退火到优化算法“退火”这个词来自材料科学。金属在高温下加热内部原子获得足够能量后可以自由移动整个结构处于一种比较无序的高能状态如果这时候缓慢降温原子就有充足时间重新排列最终形成规则且能量较低的晶格结构。反过来如果冷却速度太快原子还没排好就凝固了金属内部会残留缺陷性能自然变差。1983 年有人把这个物理过程抽象成了一套优化算法也就是模拟退火。核心思想特别直接把优化问题的一个候选解看成金属的“原子排列状态”把解的目标函数值看成系统的“能量”然后引入一个叫“温度”的参数来控制搜索行为。温度高时搜索非常激进连明显变差的解都接受温度慢慢下降接受差解的概率越来越低最终温度趋近于零算法退化成只接受更优解相当于在做一个精细的局部搜索。这个算法最迷人的地方在于它不要求目标函数可导、不要求问题有特殊结构只要你能写清楚“什么是解”和“怎么比较两个解的优劣”就能跑起来。TSP 恰好就是这样一类问题。1.2 跳出局部最优的底气Metropolis 准则模拟退火和普通爬山法最大的区别藏在一个叫 Metropolis 准则的接受策略里。假设当前解对应的目标函数值是 E候选新解的值是 E_new定义能量差 Δ E_new - E。如果 Δ 0说明新解更优直接接受如果 Δ ≥ 0也就是新解变差了算法并不会直接扔掉而是以概率p exp(-Δ / T)来决定是否接受。你可以把它理解为大雾天在山里找最低点高温时你步子迈得很大即使走到一个低洼处再往上爬也不怕因为你知道自己有足够能量翻过附近的山脊温度降下来之后你变得谨慎几乎只往明显更低的地方移动。所谓“跳出局部最优”靠的就是在高温阶段容忍那些暂时变差、但可能通向更好区域的中间解。我个人在实际调试中有一个体会很多同学把模拟退火跑出糟糕结果根本原因不是代码写错而是温度 T 和目标函数差值的数量级不匹配导致 exp(-Δ/T) 要么永远接近 1、要么永远接近 0接受准则彻底失效。这个问题后面讲参数的时候会重点展开。2. TSP 问题的建模思路2.1 旅行商问题的本质与难度TSP 的描述非常简单给定 n 个城市的位置找到一条经过所有城市、每个城市恰好访问一次、最后返回出发城市的最短回路。这个问题的关键是城市访问顺序决定了总路程而可能的顺序数量大得惊人。n 个城市的无方向回路数量是 (n-1)!/2。当 n30 时这个数字大约在 10 的 30 次方量级用普通电脑做全枚举跑到宇宙寿命结束也算不完。正因为这种“看起来简单、算起来爆炸”的特性TSP 成了各种元启发式算法的标准练武场——遗传算法、粒子群、蚁群、模拟退火几乎都要拿它做基准测试也方便和其他算法公开对比。在课程设计或者小规模验证中通常用 20 到 50 个城市比较常见。城市坐标可以直接随机生成也可以从公开的 TSP 基准测试集里选一组有现成的最优解可以对照评判算法好坏更有说服力。2.2 解的编码与总路程计算用模拟退火求 TSP第一步是定义“解”长什么样。最自然的方案是排列编码用一条向量表示城市访问顺序。比如[3, 7, 1, 5, 2, 6, 4]含义就是从城市 3 出发依次访问 7、1、5、2、6、4最后回到 3。这里每个数字必须恰好出现一次这是 TSP 解的硬性约束。总路程则是相邻城市间距离的累加和总距离 Σ d(path(i), path(i1)) d(path(n), path(1))如果用欧氏平面坐标城市 i 和城市 j 的距离就是它们坐标点之间的直线距离。MATLAB 里计算距离矩阵最省事的是用工具包函数但为了不依赖额外工具箱我更喜欢用双循环自己算一遍反正城市数量不大预处理一次、后面反复查表即可。这一步把“计算任意两条边长度”的复杂度降到了 O(1)对后续算法性能很关键。3. 直接能跑的 MATLAB 实现3.1 主程序框架温度循环加内层迭代模拟退火的标准框架是两层循环外层控制温度衰减内层在同一个温度下做多轮邻域搜索内层次数也叫 Markov 链长度。下面是完整的主程序城市数量 30参数都是经过实测比较稳妥的组合。%% 模拟退火求解TSP主程序 clear; clc; close all; % 随机生成30个城市坐标范围[0,100]x[0,100] city_num 30; cities 100 * rand(city_num, 2); % 预计算距离矩阵 dist_matrix calc_dist_matrix(cities); % 参数设置 T0 100; % 初始温度 T_end 1e-3; % 截止温度 alpha 0.98; % 降温系数 markov_len 200; % 每个温度下迭代次数 % 初始化随机排列作为初始解 current_path randperm(city_num); current_dist calc_total_dist(current_path, dist_matrix); best_path current_path; best_dist current_dist; T T0; iter 0; record []; while T T_end for i 1:markov_len % 2-opt产生新解翻转一段子路径 pos1 randi(city_num); pos2 randi(city_num); if pos1 pos2 tmp pos1; pos1 pos2; pos2 tmp; end if pos1 pos2 continue; end new_path current_path; new_path(pos1:pos2) fliplr(new_path(pos1:pos2)); new_dist calc_total_dist(new_path, dist_matrix); delta new_dist - current_dist; if delta 0 || rand exp(-delta / T) current_path new_path; current_dist new_dist; end if current_dist best_dist best_dist current_dist; best_path current_path; end iter iter 1; record(iter) best_dist; end T T * alpha; end代码里有一个细节内层每次迭代都更新最优解而不是只在温度循环结束时更新。这样即使温度降到很低、算法卡在某个区域我们依然保住了整个搜索历史上最好的结果。这种“当前解”和“历史最优解”分离的策略是模拟退火实践里一个非常重要的经验能显著提升最终输出质量。3.2 距离矩阵与总路程计算函数主程序用到了两个辅助函数。第一个是距离矩阵计算不依赖任何工具箱function dist_mat calc_dist_matrix(cities) n size(cities, 1); dist_mat zeros(n, n); for i 1:n for j i1:n d norm(cities(i,:) - cities(j,:)); dist_mat(i,j) d; dist_mat(j,i) d; end end end第二个是总路程计算输入一条路径和距离矩阵返回回路总长度function total calc_total_dist(path, dist_matrix) n length(path); total 0; for i 1:n-1 total total dist_matrix(path(i), path(i1)); end total total dist_matrix(path(n), path(1)); % 回到起点 end这两个函数都很朴素但把它们的逻辑写清楚之后后续换数据集、换目标函数都很容易。我在实际教学中见过不少同学把距离计算写进主程序里导致代码又长又难调试拆成独立函数之后主程序的逻辑一目了然排查问题也方便。3.3 邻域操作2-opt 为什么是首选模拟退火在“当前解”附近生成新解的操作叫邻域操作常见的有三种交换选两个位置互换城市、插入把一个城市挪到另一个位置、2-opt选一段子路径反转其顺序。我个人强烈推荐 2-opt尤其在城市数量超过 20 时。原因很简单对于 TSP 这类由连接顺序决定成本的问题2-opt 一次性改变了整段路径的连接关系搜索效率远高于只交换两个城市的局部微调。你可以观察一个事实——最优路径的图几乎不会出现交叉线段而 2-opt 操作天然就有消除交叉的能力。相比之下交换操作的搜索步长太小收敛起来明显更慢。当然交换操作也有其价值实现最简单适合先把整体流程跑通在一些特别看重速度、不需要高精度的场景下够用。我的建议是初学阶段先用交换跑通理解 SA 的整体流程之后立刻换成 2-opt你会明显看到收敛速度和最终精度的提升。4. 参数调节的核心经验4.1 初温 T0 怎么设初温是整个算法里最关键的参数因为它直接决定搜索前期的“探索能力”。一个实用原则是初始阶段对坏解的接受率应该保持在 80% 到 90% 左右。如果 T0 太低算法很快就丧失了跳出能力会早早掉进局部最优如果 T0 太高前期大量接受差解纯粹浪费计算量。实操中怎么估计 T0可以随机生成若干对相邻解算出它们的能量差 Δ然后取平均 Δ再把 T0 设为这个平均值的 5 到 10 倍。对 30 个坐标范围在 0 到 100 的随机城市一次 2-opt 翻转通常产生几十到上百的 Δ所以 T0 取 100 左右在大多数情况下是合理的。初温情况表现特征处理建议T0 过低收敛曲线很快变平结果明显差提高 T0或按平均 Δ 估算T0 合适前期快速下降中后期稳步收敛保持当前量级T0 过高前段长时间几乎不下降浪费时间降低 T0 或减少 Markov 链长度4.2 衰减系数 alpha 与 Markov 链长度温度衰减最常用的是乘法降温T(k1) α × T(k)。α 一般取 0.9 到 0.99α 越接近 1降温越慢、搜索越充分、耗时也越长。Markov 链长度则是每个温度下的内层迭代次数这两个参数是配合着调的α 取大时Markov 链可以适当短α 取小时需要更长的 Markov 链来补偿。总计算量可以粗略估算为外循环次数 × 内循环次数。例如 30 个城市T0100T_end1e-3α0.98降温到截止温度大约需要 568 次外循环每次 200 次内迭代总迭代约 11 万次。这个量级在 MATLAB 里跑起来非常快几秒钟就能完成。注意不要同时把 α 设成 0.999 又把 Markov 链长度设成 1000虽然结果会很好但耗时可能是原来的几十倍。要结合问题规模控制总计算量先小规模网格搜索再放大。4.3 初始解的选取随机还是贪心常见的初始化方式是randperm随机排列简单、无偏。但如果你希望在同样时间内得到更稳定的结果可以先用贪心法构造一个初始解从某个城市出发每次都选择距离最近且尚未访问的城市。这样初始解就比较接近局部最优退火过程相当于在好解附近做精修。我个人的经验是贪心初始化加稍大的 α在中等规模问题上表现非常稳定但它也有隐患——贪心解本身落在某个局部区域可能限制了前期全局探索的广度。所以在实际项目中我更常用的策略是“多次随机初始化取最优结果”。这也侧面说明模拟退火是个随机算法单次运行结果不能代表算法真实水平后面会专门讲怎么评估。5. 实验结果分析与可视化5.1 收敛曲线怎么看主程序里record数组记录了每个迭代周期之后的历史最优长度画出来就是收敛曲线。一个健康的收敛曲线通常是这样的前段急剧下降中段缓慢下降末段基本水平。如果曲线前段很快就水平了而且最终值明显不理想多半是初温太低或者降温太快。画图时如果用普通plot发现前段数值跨度太大、细节看不清可以改用semilogy把纵轴改成对数坐标收敛过程会展示得更清楚。下面这段代码可以直接接在主程序后面figure; plot(cities(best_path, 1), cities(best_path, 2), o-, LineWidth, 1.5); hold on; plot(cities(best_path(1), 1), cities(best_path(1), 2), rp, MarkerSize, 12); title([最佳路径长度: , num2str(best_dist)]); grid on; figure; semilogy(record, LineWidth, 1.5); title(收敛曲线); xlabel(迭代次数); ylabel(最优路径长度); grid on;5.2 多次运行取统计结果模拟退火是随机算法每次运行结果都有波动。正确的评估姿势是同一组参数跑 10 次记录最好值、最差值、平均值和标准差。平均值反映算法的整体水平标准差反映稳定性最差值则能帮你判断有没有偶尔坍缩到很差的解。我在调参时通常固定一个公共数据集然后先跑一组参数 10 次用平均长度作为主要指标。只比较单次最好值容易被运气干扰尤其是在课程设计答辩场景下你能拿出一张统计表比单纯展示一个最优结果有说服力得多。5.3 与贪心、随机搜索对比为了验证模拟退火确实有效可以在相同数据集上跑随机搜索和贪心算法做一个朴素的对比。随机搜索就是不停地随机生成排列、记录最小值贪心就是之前提到的最近邻构造。三者对比下来你会看到随机搜索在 30 个城市场景下表现最差贪心初期表现还行但很快止步于局部最优模拟退火在足够迭代次数下通常能稳定优于这两者。这个对比的意义不只是为了证明 SA 更强更重要的是让你直观理解“探索”的价值。贪心解很多时候已经不错了但差的那几个百分点恰恰是智能优化算法发挥作用的空间。6. 常见问题与排查技巧实录6.1 结果一直停在某个值不动这是最常遇到的问题。可能的原因有两类一类是 T0 太低算法在前期就没有足够概率接受差解失去了跳出局部最优的能力另一类是邻域操作太弱比如只用交换两个城市单个操作对路径结构改动太小导致搜索步长不足。我的排查顺序是先看收敛曲线如果曲线非常早地进入平台期说明初温或者降温速率有问题如果曲线持续下降但速度极慢再怀疑邻域操作。把 2-opt 换上去、T0 调到合适量级大多数“卡死”情况都能解决。6.2 运行时间太长怎么优化很多初学者把时间浪费在重复计算上。比如 2-opt 翻转一段路径之后没有必要重新计算整条路径的总距离这个操作的代价是 O(n)实际上翻转 pos1 到 pos2 这一段受影响的只有两条边原来连接 (pos1-1, pos1) 和 (pos2, pos21) 的边翻转后变成了 (pos1-1, pos2) 和 (pos1, pos21)其余内部边的长度和保持不变。因此可以直接做增量评估n length(path); i0 pos1 - 1; j0 pos2 1; if i0 1, i0 n; end if j0 n, j0 1; end delta - dist_matrix(path(i0), path(pos1)) ... - dist_matrix(path(pos2), path(j0)) ... dist_matrix(path(i0), path(pos2)) ... dist_matrix(path(pos1), path(j0));把原来的calc_total_dist换成这种增量计算代入delta即可速度能提升一个数量级。需要注意边界条件当 pos1 等于 1 时pos1-1 要回绕到 n当 pos2 等于 n 时pos21 要回绕到 1。另外如果翻转段覆盖了整个环pos11 且 pos2n此时路径其实没有变化应当跳过。6.3 接受率异常的表现与处理如果你在算法开头统计接受率会发现一些典型症状。初温太低时从第一步开始就几乎不接受坏解算法直接退化成爬山法初温太高时前期几乎所有坏解都被接受收敛曲线的第一段平得像一条直线等温度降到合理范围时宝贵的前期迭代已经浪费掉了。想做得更系统可以在主循环里加一个accept_count计数器统计每个温度下接受坏解的次数与总迭代次数的比值打印出来观察。理想情况下这个接受率应从高到低平滑下降一开始 0.8 左右后期趋近于 0。6.4 路径图出现明显交叉线如果你画出最优路径发现里面还有明显的交叉线段这说明收敛还不彻底。对于平面欧氏距离的 TSP最优解不会包含交叉边因为任何交叉都可以通过 2-opt 翻折变成更短的无交叉路径。看到交叉通常意味着迭代次数不足、α 还不够大或者初始解太差。碰到这种情况我一般会先把迭代次数增加 2 到 3 倍看改善幅度如果改善很小再回头检查 T0 和 α 是否匹配。记住一个原则调参的顺序是先保证“搜索充分”再追求“计算高效”顺序反了很容易浪费大量时间。最后再分享一个实际管用的小技巧跑完一轮模拟退火之后把得到的最优路径作为初始解用小的正温度再跑一轮短退火或者直接在这个解附近做 2-opt 局部搜索。这种“退火出来再打磨”的两阶段思路经常能用很小的额外开销把结果再提升一截——我自己在多次实验里都验证过效果比盲目加大迭代次数要划算得多。
阅读完成 · 觉得有帮助？
