1. 项目缘起从“旅行推销员”到代码实践最近在整理一些经典的组合优化问题案例TSP旅行商问题自然是绕不开的一座大山。这问题说起来简单一个推销员要去N个城市推销商品每个城市去一次且仅一次最后回到起点怎么走总路程最短但就是这个看似简单的描述让无数数学家和程序员“头秃”了几十年。它属于NP-hard问题城市数量一多暴力枚举所有路径(N-1)! 条在有限时间内基本不可能。所以大家转向了各种启发式算法试图在可接受的时间内找到一个“还不错”的解。在众多启发式算法里模拟退火Simulated Annealing, SA一直是我的心头好。它不像遗传算法那样需要设计复杂的交叉、变异操作也不像蚁群算法那样参数调起来让人眼花缭乱。模拟退火的核心思想非常物理、非常直观模仿金属退火过程通过控制“温度”这个参数让搜索过程既有“跳出局部最优”的探索能力又有“收敛到优质解”的开发能力。这次我就把手头一个用Python实现模拟退火求解TSP的代码拿出来结合我自己的调试和优化经验从头到尾捋一遍。你会发现实现一个能跑的SA框架可能就几十行代码但真想让它跑得好、解得快里面的门道可不少。2. 模拟退火算法核心不只是“概率接受差解”很多人对模拟退火的初印象就是“以一定概率接受更差的解从而跳出局部最优”。这个说法没错但太笼统了。要真正用好它必须理解其背后的三个核心机制以及它们是如何在代码中具体体现的。2.1 能量函数与邻域结构定义你的“问题世界”在模拟退火的语境里“能量”就是我们要优化的目标函数值。对于TSP能量就是路径的总长度。计算总长度是基础操作这里有个小技巧通常我们会预先计算好所有城市两两之间的距离矩阵dist_matrix。这样在评估一条路径的能量时只需要做N次查表和加法而不是每次临时计算距离能极大提升效率尤其是城市数量较多时。import numpy as np import math def calc_distance(city1, city2): 计算两个城市间的欧氏距离 return math.sqrt((city1[0] - city2[0])**2 (city1[1] - city2[1])**2) def create_dist_matrix(cities): 创建距离矩阵 n len(cities) dist_matrix np.zeros((n, n)) for i in range(n): for j in range(i1, n): dist calc_distance(cities[i], cities[j]) dist_matrix[i][j] dist_matrix[j][i] dist return dist_matrix def total_distance(path, dist_matrix): 计算一条路径的总距离能量 total 0.0 n len(path) for i in range(n): total dist_matrix[path[i]][path[(i1) % n]] # 最后一个城市回到起点 return total比能量函数更重要的是“邻域结构”它定义了如何从当前解产生一个新解。在TSP中最常用的邻域操作有交换Swap随机选择路径中的两个位置交换它们对应的城市。逆转Reverse/2-opt随机选择路径中的一段子路径将其顺序完全颠倒。插入Insert随机选择一个城市将其插入到路径的另一个随机位置。我的经验是对于中小规模的TSP城市数1002-opt操作通常效果最好。因为它一次性能改变路径中多个边的连接扰动更大更容易跳出局部最优的“小水坑”。而交换操作扰动较小在低温阶段进行微调时可能更有用。在实际代码中我常常将两者结合高温时多用2-opt进行大胆探索低温时引入交换进行精细调整。def generate_new_path_2opt(old_path): 使用2-opt片段逆转产生新路径 n len(old_path) new_path old_path.copy() # 随机选择两个不同的索引并确保 i j i, j sorted(np.random.choice(n, 2, replaceFalse)) # 逆转 i 到 j 之间的片段 new_path[i:j1] new_path[i:j1][::-1] return new_path def generate_new_path_swap(old_path): 使用交换操作产生新路径 n len(old_path) new_path old_path.copy() i, j np.random.choice(n, 2, replaceFalse) new_path[i], new_path[j] new_path[j], new_path[i] return new_path2.2 退火计划表控制搜索的“节奏感”这是模拟退火的精髓所在直接决定了算法的性能和最终解的质量。它主要包含四个参数初始温度T_init温度太高算法几乎完全随机搜索效率低下温度太低又容易过早陷入局部最优。一个经验公式是T_init -ΔE_avg / ln(P_init)其中ΔE_avg是随机产生一批新解时能量差新解能量-旧解能量的平均值P_init是你期望在初始时接受差解的概率比如0.8。实践中如果嫌麻烦也可以根据目标函数值的数量级进行估算比如设为目标函数值范围的若干倍。终止温度T_end通常设为一个非常接近0的正数比如1e-7。当温度低于此值时算法停止此时它基本只接受更好的解相当于一个局部搜索。温度衰减系数alpha最常见的衰减方式是T_new alpha * T_old其中alpha是一个略小于1的数如0.95到0.99。alpha越大降温越慢搜索越充分但耗时也越长。我个人的习惯是从0.99开始尝试。马尔可夫链长度L在每个温度下迭代的次数。太短搜索不充分太长浪费时间。一个常见的策略是L 100 * NN为城市数或者设置一个固定值如1000-5000并通过实验调整。class SimulatedAnnealingTSP: def __init__(self, cities, dist_matrix): self.cities cities self.dist_matrix dist_matrix self.n len(cities) self.best_path None self.best_energy float(inf) self.history {temp: [], energy: [], best_energy: []} def solve(self, T_init1000, T_end1e-7, alpha0.99, L2000): # 初始化路径随机排列 current_path np.random.permutation(self.n).tolist() current_energy total_distance(current_path, self.dist_matrix) self.best_path current_path.copy() self.best_energy current_energy T T_init while T T_end: for _ in range(L): # 以一定概率选择不同的邻域操作 if np.random.random() 0.7: # 70%概率使用2-opt new_path generate_new_path_2opt(current_path) else: # 30%概率使用交换 new_path generate_new_path_swap(current_path) new_energy total_distance(new_path, self.dist_matrix) delta_e new_energy - current_energy # Metropolis准则判断是否接受新解 if delta_e 0 or np.random.random() math.exp(-delta_e / T): current_path, current_energy new_path, new_energy # 更新历史最优解 if current_energy self.best_energy: self.best_path current_path.copy() self.best_energy current_energy # 记录当前温度下的状态用于分析 self.history[temp].append(T) self.history[energy].append(current_energy) self.history[best_energy].append(self.best_energy) # 降温 T * alpha return self.best_path, self.best_energy2.3 Metropolis准则算法跳出能力的“灵魂”接受差解的概率由P exp(-ΔE / T)决定。这是整个算法能跳出局部最优的关键。当 ΔE 0新解更好exp(-ΔE / T)大于1所以一定接受。这是“下山”过程。当 ΔE 0新解更差以概率P接受。温度T很高时即使ΔE很大P也可能不小算法有勇气跳到更远的地方随着T降低接受差解的概率越来越小算法越来越“保守”最终稳定在一个希望是全局或优质的局部最优解附近。这里有一个极易忽略的坑ΔE和T的量级必须匹配如果ΔE的典型值是几百万而T从1000开始衰减那么-ΔE/T会是一个非常巨大的负数导致exp(-ΔE/T)在绝大多数情况下计算结果为0浮点数下溢算法实际上失去了接受差解的能力退化成普通的局部搜索。因此初始温度的设置必须参考目标函数值的尺度。如果路径总距离在10^4量级T_init设在10^3量级可能就太小了。3. Python实现中的性能陷阱与优化技巧把算法思路翻译成Python代码不难但写出高效、健壮的代码需要一些技巧。下面是我在实现过程中踩过的一些坑和总结的优化经验。3.1 距离矩阵与向量化操作前面提到了预计算距离矩阵这是最重要的优化没有之一。避免了在能量评估的循环中重复调用math.sqrt和乘法运算。对于N个城市能量计算复杂度从 O(N²) 降到了 O(N)。更进一步我们可以利用NumPy的向量化操作来加速能量计算。虽然对于单次路径评估循环和向量化差别不大但在每个温度下要进行L次评估累积起来就很可观了。def total_distance_vectorized(path, dist_matrix): 使用numpy向量化操作计算总距离稍快 # 将路径索引转换为numpy数组以便高级索引 idx np.array(path) # 利用roll操作获取下一个城市的索引 next_idx np.roll(idx, -1) # 使用高级索引一次性获取所有距离并求和 return np.sum(dist_matrix[idx, next_idx])注意np.roll会产生一个新数组对于超大规模路径其开销也需要考虑。但在大多数情况下向量化版本更具可读性和一定的速度优势。3.2 路径表示与邻域操作的效率路径通常用一个列表或numpy数组表示存储城市的访问顺序。在进行邻域操作如2-opt时要特别注意避免不必要的完整列表拷贝。我最初写的generate_new_path_2opt函数是new_path old_path.copy()然后进行切片逆转。这对于Python列表是可行的因为切片操作会创建新列表。但如果old_path是numpy数组直接切片赋值new_path[i:j1] old_path[i:j1][::-1]是原地操作的一部分但最开始的copy()仍然是必须的否则会修改原路径破坏算法状态。一个更极致的优化是在高温、大量接受差解的阶段可以尝试不总是创建完整的新路径副本而是记录对当前路径的“差分”修改并在能量计算时只更新受影响的部分距离。但这会大大增加代码复杂度除非面对城市数量极大1000的情况否则收益可能不如优化其他部分明显。3.3 随机数生成与随机性控制模拟退火依赖随机数进行邻域扰动和Metropolis判断。使用np.random模块比Python内置的random模块更快尤其是在需要生成大量随机数时。另外固定随机种子对于调试和结果复现至关重要。在开发阶段设置np.random.seed(42)可以让每次运行都产生相同的随机序列这样当你修改了某个参数比如降温系数观察到的效果变化才是真实的而不是随机性带来的噪声。def solve(self, T_init1000, T_end1e-7, alpha0.99, L2000, seedNone): if seed is not None: np.random.seed(seed) # ... 其余求解代码 ...在最终多次运行取最优解时再去掉固定的种子或者使用不同的种子运行多次。4. 参数调优实战如何让算法“跑得又好又快”模拟退火的参数没有银弹需要针对具体问题进行调整。以下是我常用的调优流程和策略。4.1 初始温度的自动化估计手动拍一个初始温度很麻烦。我们可以实现一个简单的自适应方法来估计T_init。def estimate_initial_temperature(cities, dist_matrix, num_samples100, initial_accept_prob0.8): 通过采样估计初始温度 n len(cities) current_path np.random.permutation(n).tolist() current_energy total_distance(current_path, dist_matrix) delta_es [] for _ in range(num_samples): new_path generate_new_path_2opt(current_path) new_energy total_distance(new_path, dist_matrix) delta_e new_energy - current_energy if delta_e 0: # 只收集变差的能量差 delta_es.append(delta_e) # 更新当前路径继续采样 current_path, current_energy new_path, new_energy if delta_es: avg_delta_e np.mean(delta_es) # 根据公式 T -ΔE_avg / ln(P_init) T_init -avg_delta_e / math.log(initial_accept_prob) else: # 如果采样中全是更优解说明初始路径很差可以设一个较大的默认值 T_init 1000 * current_energy / n # 一个经验公式 return max(T_init, 1.0) # 确保温度为正这个方法通过随机采样一批状态转移计算变差的ΔE的平均值然后反推出能让你期望的初始接受概率如0.8成立的温度。虽然不精确但比盲目猜测要好得多。4.2 降温策略与停止准则的变体除了等比降温还有其它策略线性降温T_new T_old - dT。降温速度恒定但需要精心选择dT。自适应降温根据当前解的接受率来调整降温速度。例如如果当前温度下的接受率很高说明还没充分搜索可以慢点降温如果接受率很低说明已经接近稳定可以加快降温。实现起来稍复杂但有时效果更好。停止准则也可以更智能连续若干温度最优解未改进如果最优解连续K个温度都没有更新可以提前终止。能量变化率过低监控当前解能量的变化如果变化微乎其微也可以停止。def solve_adaptive(self, T_initNone, L2000, no_improve_limit50): if T_init is None: T_init estimate_initial_temperature(self.cities, self.dist_matrix) T T_init current_path np.random.permutation(self.n).tolist() current_energy total_distance(current_path, self.dist_matrix) self.best_path current_path.copy() self.best_energy current_energy no_improve_count 0 while no_improve_count no_improve_limit: accepted 0 for _ in range(L): # ... 生成新解并判断是否接受 ... if accepted_this_move: accepted 1 accept_rate accepted / L # 自适应降温接受率高则慢降接受率低则快降 if accept_rate 0.6: alpha 0.98 # 慢降 elif accept_rate 0.2: alpha 0.90 # 快降 else: alpha 0.95 # 中速降 T * alpha # 检查最优解是否更新 if current_energy self.best_energy: self.best_energy current_energy self.best_path current_path.copy() no_improve_count 0 else: no_improve_count 1 if T 1e-10: # 绝对温度下限 break return self.best_path, self.best_energy4.3 马尔可夫链长度与迭代平衡L的设置需要权衡。我的经验法则是与问题规模相关L k * N其中k在50到200之间。城市越多每个温度下需要更多的尝试来探索状态空间。与降温系数配合如果alpha很大如0.99降温慢每个温度下可以设置较小的L如100*N因为总迭代次数多。如果alpha较小如0.90降温快每个温度下应设置较大的L如500*N以确保在每个温度下都能充分搜索。动态调整也可以让L随着温度降低而增加。高温时进行粗搜索L可以小一些低温时进行精细搜索L增大。一个简单的动态调整可以是L_current int(L_init * (1 math.log(1 T_init / T)))这样温度越低链长越长。5. 结果可视化与算法诊断算法跑完了怎么知道它运行得好不好光看一个最终路径长度是不够的。可视化是强大的诊断工具。5.1 绘制优化过程曲线绘制能量当前解距离和最优能量随迭代或温度下降的曲线可以直观看到算法的收敛过程。import matplotlib.pyplot as plt def plot_optimization_history(sa_solver): 绘制优化历史曲线 history sa_solver.history iterations range(len(history[energy])) plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(iterations, history[energy], b-, alpha0.6, labelCurrent Energy) plt.plot(iterations, history[best_energy], r-, linewidth2, labelBest Energy) plt.xlabel(Iteration (Temperature Step)) plt.ylabel(Total Distance) plt.title(Energy Convergence) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.subplot(1, 2, 2) plt.semilogy(iterations, history[temp], g-) plt.xlabel(Iteration (Temperature Step)) plt.ylabel(Temperature (log scale)) plt.title(Temperature Schedule) plt.grid(True, linestyle--, alpha0.5) plt.tight_layout() plt.show()从曲线中我们可以看出能量曲线是否平稳下降如果当前能量曲线剧烈震荡说明温度可能还太高或者邻域操作扰动太大。最优能量曲线是否在持续改进如果在很长一段迭代后最优解都没变化可能陷入了局部最优需要考虑增加初始温度或调整邻域操作。降温曲线是否合理是否符合预期的衰减速度。5.2 绘制最终路径图将城市坐标和最终找到的最优路径画出来是最直接的成果展示。def plot_path(cities, path, titleBest TSP Path Found): 绘制TSP路径图 cities_arr np.array(cities) path_arr np.array(path [path[0]]) # 闭合路径 plt.figure(figsize(10, 8)) plt.scatter(cities_arr[:, 0], cities_arr[:, 1], cred, s100, zorder5) plt.plot(cities_arr[path_arr, 0], cities_arr[path_arr, 1], b-, linewidth1.5, zorder4) # 标注城市编号 for i, (x, y) in enumerate(cities): plt.text(x, y, str(i), fontsize12, hacenter, vacenter, colorwhite, zorder6) plt.xlabel(X Coordinate) plt.ylabel(Y Coordinate) plt.title(title) plt.axis(equal) plt.grid(True, linestyle--, alpha0.5) plt.show()通过看图可以快速判断解的质量路径是否有明显的交叉是否绕了远路这能给你直观的反馈帮助你调整邻域操作例如2-opt操作的一个重要特性就是消除路径交叉。5.3 与基准问题对比如果你知道所求解的TSP实例的最优解或已知最优解例如TSPLIB中的标准问题可以将算法结果与最优解进行对比计算近似比(你的解长度 / 最优解长度)这是衡量算法性能的客观指标。即使不知道最优解也可以运行多次模拟退火使用不同随机种子统计解的平均值、标准差和最好值评估算法的稳定性和鲁棒性。def run_multiple_trials(cities, num_trials10): 多次运行SA统计结果 dist_matrix create_dist_matrix(cities) results [] best_of_all None best_energy_all float(inf) for trial in range(num_trials): solver SimulatedAnnealingTSP(cities, dist_matrix) best_path, best_energy solver.solve(T_init1000, seedtrial) # 使用trial作为种子 results.append(best_energy) if best_energy best_energy_all: best_energy_all best_energy best_of_all best_path.copy() print(fTrial {trial1}: Best Energy {best_energy:.2f}) results_arr np.array(results) print(f\n--- Summary after {num_trials} trials ---) print(fBest: {results_arr.min():.2f}) print(fWorst: {results_arr.max():.2f}) print(fAverage: {results_arr.mean():.2f}) print(fStd Dev: {results_arr.std():.2f}) return best_of_all, best_energy_all, results_arr多次运行能帮你确认找到的好解是运气还是算法参数设置得当。如果结果方差很大说明算法对初始状态敏感可能需要增加初始温度或马尔可夫链长度来增强探索能力。6. 超越基础进阶优化思路当你掌握了基本的模拟退火实现后可以尝试以下进阶策略来进一步提升解的质量。6.1 混合策略SA与局部搜索结合模拟退火擅长全局探索但在低温末期的局部开发能力可能不如一些专门的局部搜索算法如2-opt局部搜索、Lin-Kernighan等。一个常见的混合策略是先运行完整的模拟退火过程得到一个较好的解。将这个解作为初始解运行一个贪婪的局部搜索例如反复尝试所有可能的2-opt交换只要能使路径变短就接受直到无法再改进。这相当于用SA进行“粗调”用局部搜索进行“精修”。在代码实现上可以在SA的solve方法返回前增加一个局部搜索的步骤。def local_search_2opt(path, dist_matrix): 对给定路径进行贪婪的2-opt局部搜索 n len(path) improved True best_path path.copy() best_energy total_distance(best_path, dist_matrix) while improved: improved False for i in range(n): for j in range(i2, n): # 确保片段长度至少为2 # 尝试交换边 (i, i1) 和 (j, j1) # 计算交换后距离的变化量避免重复计算整条路径 # 这里简化处理直接构造新路径并计算能量 new_path best_path.copy() # 逆转 i1 到 j 的片段 new_path[i1:j1] new_path[i1:j1][::-1] new_energy total_distance(new_path, dist_matrix) if new_energy best_energy: best_path, best_energy new_path, new_energy improved True break # 找到改进就跳出内层循环重新开始扫描 if improved: break return best_path, best_energy # 在SA求解后调用 best_path_sa, best_energy_sa sa_solver.solve() best_path_final, best_energy_final local_search_2opt(best_path_sa, dist_matrix)6.2 并行化与多起点策略模拟退火的内循环每个温度下的L次迭代是顺序的但我们可以从两个层面并行多起点并行同时从多个不同的随机初始路径开始运行独立的SA过程最后取所有结果中的最优解。这可以充分利用多核CPU并且由于初始状态的随机性更有可能找到全局最优解。可以用Python的multiprocessing或concurrent.futures模块实现。并行尝试邻域在每个温度下可以同时生成多个候选新解例如用多个线程或进程然后并行计算它们的能量最后统一进行Metropolis判断。但这涉及到状态同步实现起来比多起点策略复杂。多起点策略实现简单收益明显是我最推荐的进阶实践。from concurrent.futures import ProcessPoolExecutor def run_sa_for_seed(seed): 单个种子下的SA运行任务 np.random.seed(seed) solver SimulatedAnnealingTSP(cities, dist_matrix) path, energy solver.solve(T_init1000, T_end1e-7, alpha0.995, L1500) return path, energy def parallel_sa(cities, num_processes4, num_trials_per_process3): 并行多起点SA dist_matrix create_dist_matrix(cities) seeds list(range(num_processes * num_trials_per_process)) # 生成不同的种子 best_global_energy float(inf) best_global_path None with ProcessPoolExecutor(max_workersnum_processes) as executor: futures [executor.submit(run_sa_for_seed, seed) for seed in seeds] for future in concurrent.futures.as_completed(futures): path, energy future.result() if energy best_global_energy: best_global_energy energy best_global_path path return best_global_path, best_global_energy6.3 针对TSP的特殊邻域操作除了通用的交换和逆转TSP还有一些更高效的专用邻域操作可以在SA的框架内使用3-opt断开路径的三条边然后以另一种方式重新连接。它比2-opt的搜索空间更大扰动更强适合在高温阶段使用。双桥移动Double-Bridge Move一种特殊的4-opt移动能产生非常大的扰动常用于跳出非常深的局部最优是许多元启发式算法如迭代局部搜索ILS中的“抖动”操作。实现这些操作稍微复杂一些需要仔细处理路径片段的切割和重组但它们能显著提升算法对复杂解空间的探索能力。