1. 需求拆解TSP 问题的规模与模拟退火的位置前段时间接到一个路径规划的实验需求几十个配送点位要求规划一条经过每个点位恰好一次、最后回到起点的闭合回路并且总距离最短。这类问题平时听着简单真正动手才发现搜索空间大得吓人。我用 Python 徒手实现了一遍模拟退火算法把 TSP旅行商问题从头到尾跑通。这篇文章就把需求拆解、算法选型、代码实现、调参踩坑一次性讲清楚。适合算法入门阶段想动手写优化算法的同学、准备面试手写题的人以及想把模拟退火框架迁移到自己项目里的开发者。1.1 为什么说 TSP 是穷举没法解决的TSP 的形式化定义其实很简单给定 n 个点的坐标找到一条闭合路径让每个点都被访问一次并且路径总长度最小。坑人的地方在于搜索空间。n 个点的排列数是 n!因为闭合回路可以把起点挪到任意位置也可以正着走反着走实际不同解大约是 (n-1)!/2。n20 的时候这个数字已经超过 2.43×10^18。就算你每秒能评估一百万条路径也要七万多年才能把所有路径看一遍。这还只是 20 个点。真实场景里动辄 50 个、100 个点你可以直接死心任何穷举思路都不可行必须用启发式方法在“够好”的解附近搜索。所以 TSP 才成了元启发式算法的经典试验田它简单到任何人五分钟就能理解问题又难到可以拿来做各种智能优化算法的 benchmark。你写的模拟退火、遗传算法、蚁群哪怕只是调整了参数都可以用同样一份城市数据直接比出优劣。1.2 为什么最终选了模拟退火而不是爬山、遗传或精确算法我最早想的是先试试简单方案最近邻贪心每次找离当前点最近且没走过的点。这个方法非常快但结果短视经常在后半段被迫走很长的边路径交叉严重。然后我也想过用遗传算法编码、选择、交叉、变异、锦标赛、精英保留参数比模拟退火多一倍调起来很痛苦。精确算法更不用提动态规划在 n 小于 20 时确实能求出严格最优城市一多内存和时间一起爆炸。最终选了模拟退火核心原因是它在“实现成本”和“解质量”之间平衡得最好。爬山法只会接受改进解一旦爬到某个局部最优就彻底被困住模拟退火的灵魂在于温度高的时候允许接受一些变差的解这样就有机会从局部最优的坑里爬出去随着温度逐步降低接受变差解的概率越来越小最后收敛到某个稳定区域。还有一个很实际的理由模拟退火的框架和 TSP 的邻域结构配合得非常自然。TSP 的经典邻域操作 2-opt 只涉及两条边的变化可以写出 O(1) 的增量式距离差计算这让迭代几十万次也变得可行。第三方优化库里虽然有现成的退火函数但它面向的是连续参数空间用来做离散组合优化还要自己实现邻域和回调函数反而绕远路。自己徒手写一遍你能精确控制每一步在干什么后面扩展到带时间窗、带容量约束的问题时也更容易加约束判断。1.3 宏观流程从城市坐标到最终路径整条流程可以拆成六步先建立整体印象后面分节展开准备城市坐标数据计算两两之间的欧氏距离得到距离矩阵生成初始解通常有随机排列和最近邻两种方式设定模拟退火参数初始温度、降温系数、终止温度、内循环长度在内循环里反复做 2-opt 扰动计算新解距离差按模拟退火接受准则决定是否替换当前解每次外循环结束温度乘上降温系数重复直到温度低于阈值输出历史最优路径并可视化查看路径形状和收敛曲线这个流程写熟之后你会发现它不只适用于 TSP。任何组合优化问题只要能定义解的表达、能定义目标函数、能定义一个邻域操作模拟退火这套“高温乱跳、低温精修”的思路都可以搬过去用。2. 核心环节设计解表达、目标函数、邻域与温度策略很多人写模拟退火容易卡住不是卡在算法本身而是卡在最基础的这两件事上解该怎么表示邻域该怎么生成。这两件事定好了主循环反而只是几十行代码的事。我按自己习惯的顺序讲从数据结构一路讲到温度曲线。2.1 解应该如何表达从“城市顺序”到“闭合回路”TSP 的解我直接用 Python list 表示里面按访问顺序存每个城市的编号。比如 tour [0, 3, 1, 2] 就表示从城市 0 出发先去 3再去 1再去 2最后从 2 绕回 0。这里容易忽略的是“闭合”这两个字。很多初学者只算相邻两点的距离之和忘了把最后一段 tour[-1] 回 tour[0] 的距离加上。在代码里距离函数要把这条闭合边算进去才算完整否则你会得到一个“开放路径”的最小值和真正的 TSP 根本不是一个问题。还有一个值得注意的细节闭合环路上起点其实是无关紧要的。tour [0, 3, 1, 2] 和 tour [3, 1, 2, 0] 是同一条回路反向走也是同一条回路。模拟退火在搜索过程中没必要刻意消除这种等价性搜索算法自己会处理你只需要心里清楚表面上看状态数是 n!因为闭合性和方向性实际有效解要少一个数量级。我不推荐用“城市对”的边集来表达解因为保证不了边集能构成一条哈密顿回路检查合法性会很麻烦。老老实实用一个顺序数组任何排列天然合法生成新解的操作也直观。2.2 目标函数与距离矩阵目标函数就是总路径长度f(tour) Σ_{k0}^{n-1} dist[tour[k]][tour[(k1) mod n]]距离我用欧氏距离也就是地图上两点的直线距离。因为坐标是随机的欧氏距离没有方向和道路的限制是 TSP 最常见的假设之一。如果你手里是经纬度坐标可能需要调用地理距离公式如果是道路网络就得用真实路网距离了。实现上有一个决定性能的细节必须提前算好 n×n 的距离矩阵而不是在每次评估时现算 sqrt。原因是模拟退火一次完整运行可能要评估几十万次候选解每次都现算平方根的开销完全没法接受。距离矩阵本质上是一个静态的查找表它把“计算”变成了“查值”。网络上有现成坐标数据集比如某几个城市间的标准 TSP 实例你可以直接拿来和别人的结果对比自己做实验时我习惯用随机坐标方便复现也方便调参。2.3 邻域生成策略交换、插入、2-opt 反转模拟退火的“搜索方向”靠邻域操作来推动。TSP 常用的有三种交换swap随机选两个城市交换它们的位置插入insertion随机选一个城市插入到另一个位置2-opt 反转随机选两个下标 i、j把两者之间的路径段落反转我一开始觉得交换最简单代码一行就能写完。但实验下来效果最差因为它会让路径变得“破碎”产生很多交叉需要更多轮迭代才能恢复。插入策略生成解的多样性还可以但同样不够精细。2-opt 是 TSP 里最经典的邻域操作它的含义是删掉两条旧边 (a,b) 和 (c,d)补上两条新边 (a,c) 和 (b,d)同时把 b 到 c 之间的城市顺序反转。这会让交叉的路径自然被解开而且一次只影响两条边配合增量式距离计算性能极高。选 2-opt 而不是其它策略本质是因为它对应了 TSP 路径的一个基本局部操作把一条自我交叉的路线改成不交叉的路线。TSP 最优解的曲线通常是不自交的所以 2-opt 可以非常高效地修正路径形状。2.4 温度调度与接受准则模拟退火的接受准则长这样如果新解变好也就是 Δ 0无条件接受如果新解变差也就是 Δ 0以 p exp(-Δ/T) 的概率接受。这个概率公式就是玻尔兹曼分布的形式。温度 T 高的时候exp(-Δ/T) 接近 1变差的解大概率被接受温度低的时候exp(-Δ/T) 接近 0只有很小的概率接受变差。这就像系统还“热”的时候允许乱跳冷下来后只敢小范围修整。初始温度怎么设是个经验问题。我常用的办法是先随机生成 100 个邻域操作统计距离差的平均值 avgΔ再按期望接受率 0.8 到 0.9 反推初始温度。公式是 T0 -avgΔ / ln(accept_rate)ln(0.85) 约等于 -0.1625所以 T0 大约是 avgΔ 的 6 倍。这样做比随手写一个 100 要可靠得多能保证系统开局有足够活力。降温系数我用几何降温T ← α × T。α 的取值强烈影响效果和耗时。α 太小比如 0.9温度从 100 降到 1 只要几十轮搜索远没充分就冻住了α 太大比如 0.9999效果可能很好但运行时间很长。一般城市规模 30 到 100 的时候α 取 0.99 到 0.998 都在合理区间。终止温度设成 1e-3 就够因为到这个量级 exp(-Δ/T) 几乎都是 0系统已经退化成纯爬山法再降下去意义不大。3. 从零开始写 Python 实现到核心环节了。我尽量把代码写得直白不做花哨优化但增量式距离计算一定要保留因为它直接影响你能不能跑 50 个城市以上的例子。3.1 准备数据与验证思路用 numpy 生成随机城市坐标然后计算距离矩阵。np.random.default_rng 比全局 np.random 更可控支持单独指定种子跑不同实验互不干扰。import numpy as np def gen_cities(n, seed42): rng np.random.default_rng(seed) return rng.random((n, 2)) def distance_matrix(cities): diff cities[:, None, :] - cities[None, :, :] dist np.sqrt((diff ** 2).sum(axis-1)) return dist cities gen_cities(32, seed42) dist distance_matrix(cities)n32 是一个很舒服的实验规模既有一定的搜索难度跑起来又不会太慢肉眼也能看出路径形状好坏。距离矩阵用广播一次性算好这里的 diff 是 (n, n, 2) 的三维数组稍微费内存但 n1000 以内都不是问题。验证思路我建议这样先造几个排成直线的点比如 (0,0), (1,0), (2,0), (3,0)那么最优闭合回路必然是 0→1→2→3→0总距离是 2×(n-1)。如果程序在这个数据上算出了别的结果说明距离计算或者路径表示一定有 bug。这个验证只要十秒钟值得在写主循环前做一次。3.2 初始化随机排列与最近邻模拟退火对初始解不是特别敏感但一个合理的初始解能省下大量早期探索。最近邻算法简单粗暴从城市 0 出发每次去最近的未访问城市直到全部走完。def nearest_neighbor_tour(dist, start0): n len(dist) unvisited set(range(n)) unvisited.remove(start) tour [start] cur start while unvisited: nxt min(unvisited, keylambda x: dist[cur][x]) tour.append(nxt) unvisited.remove(nxt) cur nxt return tour我实际对比过对同样一批 32 个随机点随机初始解的总距离大约是 12 到 15最近邻初始解通常能压到 6 到 8差距接近一倍。虽然模拟退火最终能修正初始解的质量但初始解越差需要的高温探索期就越长。如果你想观察模拟退火“从很差到很好”的完整过程建议用随机排列如果你想快速得到不错的结果用最近邻初始化更省时间。3.3 核心循环2-opt 增量式模拟退火主循环之前先写两个辅助函数一个是全量计算路径距离用于初始化和最终打印另一个是增量式计算 2-opt 反转后的距离差。2-opt 增量式的核心思想是反转 tour 中从 i1 到 j 的这一段改变的边只有两条。原本连接 a 和 b 的边、连接 c 和 d 的边被删除替换成 a 到 c、b 到 d 的两条新边。于是距离差就是delta dist[a][c] dist[b][d] - dist[a][b] - dist[c][d]这个 O(1) 的 delta 是整个程序性能的关键。如果不做增量式每次重新算整条路径距离是 O(n)内外循环一叠加60 万次评估就会慢到不可接受。def tour_distance(tour, dist): total 0.0 for k in range(len(tour)): total dist[tour[k]][tour[(k 1) % len(tour)]] return total def two_opt_delta(tour, i, j, dist): n len(tour) a tour[i] b tour[(i 1) % n] c tour[j] d tour[(j 1) % n] return dist[a][c] dist[b][d] - dist[a][b] - dist[c][d]主循环这么写def simulated_annealing(dist, n, init_T100.0, alpha0.995, T_min1e-3, L200, seed42, use_nnFalse): rng np.random.default_rng(seed) if use_nn: current_tour nearest_neighbor_tour(dist) else: current_tour list(rng.permutation(n)) current_dist tour_distance(current_tour, dist) best_tour current_tour[:] best_dist current_dist T init_T history [] while T T_min: for _ in range(L): i int(rng.integers(0, n)) j int(rng.integers(0, n)) if i j: continue if i j: i, j j, i if (i 1) % n j or (j 1) % n i: continue delta two_opt_delta(current_tour, i, j, dist) if delta 0 or rng.random() np.exp(-delta / T): current_tour[i 1: j 1] current_tour[i 1: j 1][::-1] current_dist delta if current_dist best_dist: best_dist current_dist best_tour current_tour[:] history.append(best_dist) T * alpha return best_tour, best_dist, history几个边界条件一定要处理。第一i 和 j 相等时没有意义。第二如果 i1 等于 j说明要反转的部分只有一个城市反转等于没动。第三如果 j1 取模后等于 i说明选的区间覆盖了几乎整条回路反转后等价于整条反向走白白浪费一次迭代。这三个条件过滤掉无效操作不仅节约时间还能避免一些隐蔽的路径错乱。增量式更新 current_dist 是整个程序最容易被忽略的细节很多人接受到新解后习惯重新调用 tour_distance 再算一次。我强烈建议不要这样写。delta 已经算出来了反转操作真的被执行时当前距离就是 current_dist delta直接加就好少掉一次 O(n) 遍历。3.4 可视化与结果解读光看一个最终距离数字还不够直观画出路径和收敛曲线才能确认识别算法到底做了什么。import matplotlib.pyplot as plt def plot_result(cities, tour, history): loop tour [tour[0]] fig, axes plt.subplots(1, 2, figsize(12, 4.5)) axes[0].plot(cities[loop, 0], cities[loop, 1], o-, markersize4) for idx, (x, y) in enumerate(cities): axes[0].annotate(str(idx), (x, y), fontsize8, textcoordsoffset points, xytext(4, 4)) axes[0].set_title(best tour) axes[1].plot(history) axes[1].set_title(best distance) axes[1].set_xlabel(outer iteration) axes[1].set_ylabel(distance) plt.tight_layout() plt.show()我跑的一组固定种子实验里32 个随机点的模拟退火结果随机初始解大约从 13.2 降到 5.1最近邻初始化大约从 7.8 降到 4.9两条收敛曲线都能看到明显的快速下降段和后续平缓抖动段。最终路径画出来没有交叉连线分布得很均匀。从收敛曲线上你能直接判断温度策略是否合理如果曲线从很早就变成一根直线说明系统早就冻住了大概率是初始温度太低或者降温太快如果曲线一直大幅震荡说明温度太高收敛不足需要调低初始温度或者多跑几轮。4. 调参与排查把模拟退火调到能用的细节写完第一版代码能跑出结果只是个开始。真正决定模拟退火好不好用的是你愿不愿意花时间调参数。这一节集中讲我实际踩过的坑和常用的调试手段。4.1 收敛太快或者陷入局部最优怎么排查模拟退火最常见的问题就是结果很差路径里还有大量交叉。原因基本逃不出这三个初始温度太低、降温系数太小、内循环长度太短。我把现象和排查方向整理成了一张表方便你照着一项项对照。现象可能原因处理办法结果路径有大量交叉迭代不充分搜索没铺开提高初始温度或把 alpha 从 0.995 调到 0.998高温期接受率就偏低初始温度太低按平均 delta 和期望接受率 0.85 反推 T0低温期还频繁接受劣化解终止温度太高或 L 太长调小 T_min或者减小 L运行时间过长L 太大或 alpha 太接近 1开启增量式 delta缩小 L 到 100每次运行结果差异巨大随机种子不同导致随机性过大调试期固定 seed最终对比跑多轮取最优第三个现象很多人没意识到。如果每个温度循环的长度 L 特别大低温阶段即使温度已经很低系统也在大量接受劣化解整体收敛就会变慢。温度和循环长度是配套调整的不要单独猛调一个。4.2 运行性能从重新求距离到 O(1) 增量更新我最初实现为了省事每次生成新解后直接调用 tour_distance 算一遍。32 个城市还能等一换到 100 个城市就明显卡顿。原因很简单一次 2-opt 反转虽然只影响两条边但全量重算会遍历整条路径的所有边复杂度是 O(n)。假设外循环 2000 轮每轮 200 次扰动一共 40 万次评估单次 O(n) 的代价乘上 n100就是 4000 万次距离读取。而用增量式计算每次只查距离矩阵的四个元素做一次浮点加减复杂度 O(1)差异能有几十倍。我自己做过的统计对比同一个 100 城市实例全量重算版本要跑接近 30 秒增量式版本三秒不到。所以 two_opt_delta 这个函数表面看起来只是个小优化实际上它决定了这种模拟退火流程能不能扩展到实际规模。如果你以后要把 n 加到 500 甚至 1000这个增量式的思路也依然适用。4.3 三个实战小技巧第一个技巧用直线数据验证正确性。把所有城市放在同一条直线上例如坐标是 (0,0), (1,0), (2,0), (3,0), (4,0)那么唯一的最优解只能是沿着直线走一个来回总距离是 2×(n-1)。用这个测试用例跑一遍如果程序输出的不是这个值那一定不是模拟退火调参的问题而是代码本身有 bug趁早修不要带病调参。第二个技巧调试期固定随机种子正式实验期再放开。固定种子能保证每次运行结果可复现改一个参数后你看到的效果差异就是真实差异而不是随机波动。我在调参时会先用 seed42 跑通确认无误后换用不同种子跑多轮取最低的那条路径作为最终输出。第三个技巧在循环里记录接受率。你可以在每个温度区间统计“接受劣化解的次数 / 总尝试次数”接受率太高说明系统太热接受率趋近于零说明已经冻结。这个指标比只看最终路径更容易定位问题。我在自己的代码里习惯把每次外循环的 T、best_dist、accept_count 都输出调试效率能提升一个档次。最后再分享一个让我印象深刻的教训。有一次我实现带时间窗的变种问题时图省事没有单独写增量式更新而是沿用了旧的 delta 公式结果边界情况完全对不上运行结果一直比全量计算差。后来我把两种结果逐项比对才发现问题出在反转区间和闭合边界的相互作用上。从那以后凡是修改邻域操作我一定会先用“直线数据 全量计算”做基准对比再放心调参。希望这篇实战过程也能帮你少走这些弯路。