模拟退火算法详解：从物理直觉到Python求解TSP

前阵子做项目连续两个晚上栽在同一个问题上目标函数有一堆局部最优贪心算法一进去就出不来换个起点换个效果跑十次能给你十种答案。最后把我从这种泥潭里拉出来的就是模拟退火算法Simulated Annealing简称SA。这个名字听起来挺玄乎其实就是从一个非常朴素的物理现象里偷师学来的优化策略。它能容忍往坏处走这正是它能跳出局部最优的关键。这篇文章我打算抛开教科书式的推导从一个工程师的角度把模拟退火讲透物理直觉、数学逻辑、参数调优再配一个完整的Python实现——用SA求解旅行商问题TSP代码可以直接复制运行看完就能套用到你自己的优化场景里。1. 先从一块烧红的铁说起模拟退火背后的物理直觉1.1 金属退火与优化问题的对应关系搞金属材料的人都清楚一个工序把金属加热到高温然后让它在环境里缓慢冷却。这个过程叫退火。高温状态下金属内部的原子动能很大排列状态非常混乱随着温度逐步降低原子热运动变缓会慢慢趋于能量最低的晶格结构。温度降得越慢最终得到的晶体结构就越完美内能也越低。1983年Kirkpatrick等人把这个现象抽象成了优化算法。你看物理世界和数学优化问题之间几乎可以一一对应物理退火优化问题原子的一个排列状态一个可行解物体的内能目标函数值当前温度控制参数T缓慢降温迭代中逐步降低T内能最小的晶格结构全局最优解这个对应关系不是硬凑的。物理系统在退火过程中遵循热力学规律算法则用概率分布来模拟原子运动。在高温阶段允许解的能量剧烈波动对应算法大幅探索解空间低温阶段波动变小对应算法逐渐收敛到低能区。1.2 为什么爬山法会输给一个愿意走回头路的算法先看爬山法它是很多人的第一反应每次只往目标函数更优的方向走没有更好的相邻解就停止搜索。这个方法思路很直但缺陷致命。设想一个双峰地形你站在一个低坡上前方不远处是个小山峰再远处是更高的主峰。从低坡出发爬山法会顺着山坡一口气爬到小山峰顶然后发现四周都是下坡于是宣布我已经找到最优。问题是小山峰只是局部最优。模拟退火为什么能破解这个局面因为它允许算法在某个温度下以一定概率接受一个目标函数值更差的解——说白了它允许自己下坡。当程序站在小山峰上时它有机会向下走一段钻进山坳再重新往上爬最终翻过另一侧的高峰。这种用暂时的退步换取长远的进步的思路在传统局部搜索里是根本没有的。这就引出整个算法最核心的一个概率判定准则。1.3 Metropolis准则差解也有机会活下来设当前解的目标函数值是E_old通过邻域操作产生的新解目标函数值是E_new令差值ΔE E_new - E_old。如果新解更优即ΔE 0没说的直接接受。如果新解更差即ΔE 0也不直接拒绝而是以概率P接受它P exp(-ΔE / T)T是当前的温度控制参数。这个公式有两个关键特征T很大的时候指数部分接近0P接近1差解几乎都能被接受算法行为接近随机搜索。T逐渐降低P越来越小接受差解的概率越来越低。最终T趋于0时算法退化成爬山法。有一个特别直观的生活类比你现在月薪8000摆在面前一份月薪7500的工作如果你刚毕业、处在高温度状态你可能会想试两个月看看发展空间但如果你已经在这个行业积累十年、当前状态稳定也就是低温度状态大概率不会为了少500块钱挪窝。模拟退火就是这么决定要不要接受一个不如当前的解。算法靠这一个小小的概率公式就同时拿到了前期全局探索和后期局部精化的能力。它既不是纯随机撞运气也不是死脑筋只往高处走。2. 算法主流程与四个关键参数手写之前必须想明白的事很多教程把模拟退火讲成一个很短的伪代码但实际应用时坑几乎全藏在参数的选取上。主流程我写下来不超过十行真正的难点在动手之前就必须想清楚的那几个决策。2.1 核心循环到底在做什么模拟退火的主流程可以概括成下面几步初始化一个解x确定初始温度T0降温系数alpha每个温度下的迭代次数L终止温度T_min。如果T T_min进入当前温度的等温循环在当前解的邻域里生成一个新解x。计算ΔE f(x) - f(x)。若ΔE 0接受x为新解否则生成随机数r random()若r exp(-ΔE / T)也接受x。循环L次。记录到目前为止的最优解。降温T T * alpha。回到第2步直到满足终止条件。用代码抽象一下就是current init_solution() T T0 while T T_min: for _ in range(L): new neighbor(current) delta objective(new) - objective(current) if delta 0 or random.random() math.exp(-delta / T): current new T * alpha外层循环控制降温内层循环保证同一个温度下有足够多的采样次数让系统接近当前温度下的热平衡。这两个循环对应了物理退火的两个阶段温度降到多少以及在这个温度下停留多久。而这两件事分别由参数控制。2.2 初始温度 T0给多少才算烧透初始温度决定了算法最前期的探索强度。如果T0太小Metropolis准则从一开始就近乎严格拒绝差解算法直接退化成爬山法允许下坡的能力形同虚设。反过来T0太大前期接受几乎所有解看起来是在全局搜索实际上是在做大量无方向性的随机游走白白消耗算力。实际工程里我通常不直接拍脑袋定T0而是先做一个小实验随机生成一批解统计相邻解之间目标函数差值的量级用这个量级反推T0。假设采样得到的平均|ΔE|大约是mean_delta想让初始接受率大约为p0例如0.8可以取T0 ≈ -mean_delta / ln(p0)举个例子如果mean_delta 50希望初始接受率0.8ln(0.8)约等于-0.223那么T0大约给224。这个温度下算法前期既不会像醉汉一样完全乱走又保留了充分的跳跃能力。2.3 降温系数alpha与每温迭代次数L火候要慢停留要够温度控制里有两个概念大家容易混淆降温和适应。降温的速率由alpha决定每个温度下停留多久由L决定。物理退火里温度下降太快叫淬火金属会变脆算法里温度降太快也会得到一个脆弱的结果——特别容易卡在局部最优。我常用的经验值alpha取0.9到0.999之间。0.9降温非常快适合快速验证代码逻辑0.99以上是长期优化任务的标配0.999适合大规模组合优化问题。L取50到500之间。太短系统在当前温度下还没稳定就仓促降温等温过程形同虚设太长算力开销大尤其当目标函数计算本身就很贵的时候。温度下降序列的选择也是可以做的。除了几何降温T T * alpha还有线性降温T T - step等。实践下来几何降温最实用因为它前期温度高时下降快能快速穿过纯随机阶段到后期温度低时下降放缓把宝贵算力留给精细搜索符合先粗后精的调优思路。2.4 邻域算子与终止条件这两个环节最容易被低估参数设好了容易在邻域算子这里栽跟头。邻域算子决定了新解从哪里来。连续优化问题里最常用的是高斯扰动x x N(0, σ)σ取多大直接影响搜索效果。σ太大新解跑到天边概率接受机制很容易变成纯随机σ太小每一步都缩在小范围内低温阶段根本出不了包围圈。我会让σ与温度联动温度高时用较大的扰动温度低时自动收窄类似先广撒网、后精准捕捞。组合优化问题则要针对问题设计邻域。比如TSP里最经典的是2-opt算子、交换算子、插入算子后面实战部分展开细说。终止条件常见的三种温度降到阈值、连续若干轮最优解没有提升、达到总迭代次数上限。个人建议别只依赖单一条件。温度阈值设定得太小后期纯属空转设得太大又可能太早结束。我通常三选二组合T_min设成初始温度的万分之一级别同时加一个连续200轮最优解无提升就提前终止的兜底。3. Python实现用模拟退火求解旅行商问题TSP3.1 为什么选TSP当示例旅行商问题Traveling Salesman Problem是组合优化领域最经典的问题之一给定若干个城市和两两之间的距离找一条经过所有城市且每个城市只经过一次的最短回路。这个问题的直观性极强路径长度一眼就能看出来好坏同时它又是典型的NP难问题当城市数量到20个以上时穷举所有排列已经完全不现实。用这个例子讲模拟退火既能展示算法如何解决组合爆炸问题又方便把每次求解的过程可视化。3.2 城市数据与目标函数的准备先用numpy随机生成20个城市坐标范围取0到100。目标函数就很简单给定一条访问顺序计算整条回路的欧氏距离总和。注意一定不要漏了最后从终点返回起点的距离这个细节很多人一开始会忽略。import math import random import numpy as np def generate_cities(num_cities, seed42): rng np.random.default_rng(seed) return rng.random((num_cities, 2)) * 100 def route_distance(cities, route): ordered cities[route] diff np.diff(ordered, axis0) total np.sum(np.sqrt(np.sum(diff ** 2, axis1))) total np.sqrt(np.sum((ordered[-1] - ordered[0]) ** 2)) return total这段代码里route是城市索引的排列例如[0, 3, 1, 5, ...]。np.diff算的是相邻城市之间的坐标差用勾股定理算出距离最后单独加上最后一个城市回到第一个城市的那一段。3.3 两种邻域操作交换与2-opt逆序邻域操作是模拟退火的手脚它的设计直接决定搜索效率。我用两种常用的算子交换算子随机挑两个位置把这两个位置上的城市互换。2-opt算子随机挑两个位置把两个位置之间的那一段访问顺序整体反转。2-opt是TSP里非常经典的局部改善技巧它在路径不交叉方面有天然优势比交换算子更适合作退火过程中的邻域操作。下面的实现用的是2-opt代码更简洁实际效果也明显好。def two_opt_neighbor(route): new_route route.copy() i, j sorted(random.sample(range(len(route)), 2)) new_route[i:j 1] new_route[i:j 1][::-1] return new_route随机取两个索引后排序保证i小于j。切片反转之后直接赋值回去速度快而且不会破坏城市覆盖的完整性。3.4 完整可运行的Python代码把上面的函数整到一起再写主循环。为了让大家能直接看到收敛趋势我在每次降温后都记录一次当前最短路长最后画出来。import math import random import numpy as np import matplotlib.pyplot as plt def generate_cities(num_cities, seed42): rng np.random.default_rng(seed) return rng.random((num_cities, 2)) * 100 def route_distance(cities, route): ordered cities[route] diff np.diff(ordered, axis0) total np.sum(np.sqrt(np.sum(diff ** 2, axis1))) total np.sqrt(np.sum((ordered[-1] - ordered[0]) ** 2)) return total def two_opt_neighbor(route): new_route route.copy() i, j sorted(random.sample(range(len(route)), 2)) new_route[i:j 1] new_route[i:j 1][::-1] return new_route def simulated_annealing(cities, T0100, alpha0.995, T_min1e-3, max_iter_per_temp200, seed0): random.seed(seed) n len(cities) current list(range(n)) current_dist route_distance(cities, current) best current.copy() best_dist current_dist T T0 history [] while T T_min: for _ in range(max_iter_per_temp): new two_opt_neighbor(current) new_dist route_distance(cities, new) delta new_dist - current_dist if delta 0 or random.random() math.exp(-delta / T): current new current_dist new_dist if current_dist best_dist: best current.copy() best_dist current_dist history.append(best_dist) T * alpha return best, best_dist, history if __name__ __main__: cities generate_cities(20) best_route, best_length, history simulated_annealing(cities) print(最短路径长度:, round(best_length, 4)) print(最优路径:, best_route) plt.plot(history) plt.xlabel(降温轮次) plt.ylabel(当前最优路径长度) plt.title(模拟退火求解TSP收敛曲线) plt.show()运行这个脚本需要先安装numpy和matplotlibpip install numpy matplotlib3.5 代码里最容易写错的三处细节第一处delta的符号。模拟退火习惯约定最小化问题所以delta new_dist - current_dist。如果你在最大化问题里套用同一个公式必须记得取反。搞反了的话算法会专门挑更差的解效果完全相反。第二处更新best的时候要存路径的副本。Python里赋值是用引用直接把best current后面current一变best也跟着变最后输出的最优路径根本不是你记录的那个值。我用best current.copy()就是防这个坑。第三处判断接受条件时注意exp(-delta / T)的数值稳定性。当T非常小且delta不为0时exp的参数可能是一个很大的负数结果直接下溢为0这没问题但如果delta是负值程序走不进这个分支因为delta 0的分支已经先处理了。理论上不会出现exp接受更优解的情况但如果代码写得不小心把两个分支的顺序弄反就可能出事故。4. 实测结果、调参经验与三个值得尝试的进阶方向4.1 收敛曲线不会骗人一次运行的完整过程代码默认是固定随机种子在我本机运行的结果比较稳定20个城市随机生成的初始路径总长度大约在800出头模拟退火收敛完以后路径长度落在370到400这个区间。如果你用自己的机器跑数值不会完全一样因为城市坐标虽然固定了但算法内部随机过程仍然会有细微差别整体量级和收敛趋势是基本一致的。把history画出来能看到一个典型的收敛过程开始几轮温度高路径长度下降极快从800多一路掉到500左右中间阶段曲线平滑地往下走到后面温度已经很低曲线几乎变成一条平线。这种形态说明参数设置基本合理——前期探索充分中期稳步收敛后期没有明显波动浪费算力。如果你画的曲线在中后段还是剧烈锯齿状说明温度降得太慢或者接受率一直偏高如果前期下降之后就直接躺平温度还很高就不再变化那多半是初始温度不够或邻域算子力度太小。4.2 关键参数对结果的影响实测对照表下面这个表总结了我自己调参时总结出来的规律。注意这里的数值是20个城市左右的TSP问题尺度不同问题规模需要放缩但趋势方向是一致的。参数常用范围设得偏大设得偏小初始温度T0按ΔE量级估算前期随机游走占比过大有效搜索时间被压缩高温阶段缺失直接进入局部精化容易早熟降温系数alpha0.95 ~ 0.999降温慢耗时成倍增加降温快搜索粗糙结果稳定性差每温迭代次数L50 ~ 500计算开销大收益边际递减等温不充分温度下降时系统还没稳定最小温度T_min1e-2 ~ 1e-6过早结束低温精修不足后期长时间空转改进微乎其微还有一个容易被误解的点参数是联动的不是孤立的。alpha取0.999时L可以适当取小一些alpha取0.95时L最好大一些。物理退火讲究慢降温、快平衡和快降温、慢平衡之间的取舍算法也是同一个道理。4.3 进阶方向重退火、并行SA与局部搜索精炼基础版模拟退火跑通之后想在实际项目里达到更好的效果通常有三条路可以走。第一是重退火reheating。当连续很多轮温度下降时最优解都没有任何更新说明算法可能被困在某片区域。这时候人为把温度重新抬高到中间水平给系统一次重新加热的机会相当于换一个能量状态重新开始搜索。这个策略在复杂问题里非常有效尤其适合目标函数存在大量平台区域的场景。第二是并行模拟退火。单次退火能探索的空间有限做成并行版本则是同时从多个初始解出发各自独立跑退火再定期把全局最优广播给所有线程参考。这样既增加了多样性又提高了找到全局最优的概率。Python里可以用multiprocessing或concurrent.futures实现思路不复杂效果立竿见影。第三是退火结束后接局部搜索精炼。模拟退火本身擅长找到一片好区域不擅长在这片区域里抠到极限精度。退火结束之后对结果再跑一遍2-opt直到无法改进往往能再把目标函数往下压一截。工程上这种元启发式粗搜索 局部精搜的组合才是标准打法。我个人在实际优化项目里最常用的组合是多个初始解并行退火得到候选解后统一做一次局部精炼最后选目标值最好的结果。这样既能保证全局探索的广度又能保证最终解的质量稳定。最后再分享一个小技巧如果在调参时不确定当前温度设置是否合适很值得在循环里实时统计接受率——接受率一直偏高说明温度太高前期探索过度接受率骤降说明温度降得太快后半程几乎只在做爬山。接受率从高到低平滑过渡才是一张健康的退火曲线。这个指标比单纯看最终结果更能帮你诊断算法也更适合快速定位问题出在哪个环节。