1. 从“烧铁淬火”到全局寻优模拟退火算法初印象如果你在数学建模、算法竞赛或者优化问题的实际工程中摸爬滚打过一阵子大概率会听过“模拟退火”这个名字。我第一次接触它是在为一个物流中心的选址问题挠头时。当时试遍了梯度下降、遗传算法效果总是不尽如人意要么陷入局部最优解出不来要么收敛速度慢得让人心焦。直到一位前辈扔过来一句“试试模拟退火吧这玩意儿有时候能给你惊喜。” 我抱着试试看的心态去研究结果发现这个算法的思想之巧妙简直是把物理世界的智慧用到了数学优化里。模拟退火算法英文是Simulated Annealing简称SA。它的核心灵感来源于冶金学中的“退火”工艺。想象一下铁匠打铁先把铁块加热到高温这时铁内部的原子处于高能、活跃的无序状态然后缓慢地、有控制地降温原子会在这个过程中逐渐找到一个能量最低、结构最稳定的排列方式从而得到一块坚韧的好钢。SA算法就是把我们要优化的目标函数比如成本、距离、误差想象成这个系统的“能量”通过模拟“加热”和“缓慢降温”的过程让解在解空间中“跳跃”有策略地接受一些暂时看起来更差的解从而有机会跳出局部最优的陷阱最终逼近全局最优解。这听起来有点反直觉对吧我们找最优解不就应该一路向下只接受更好的解吗但现实世界中的很多优化问题其解空间就像一片连绵起伏的山脉布满了深坑局部最优和高峰。如果你只肯往下走那很容易掉进一个离你起点最近的深坑里再也看不到远处那个更深的峡谷全局最优。模拟退火给你的就是一把“登山杖”允许你偶尔往上爬几步去看看山那边的风景。它特别适合解决那些解空间复杂、存在大量局部最优解的组合优化问题比如旅行商问题、调度问题、布局问题以及函数优化等。在数学建模中当你的模型目标函数复杂、非线性、不可微或者约束条件难以处理时SA往往是一个值得考虑的“暴力美学”式工具。2. 算法核心温度、扰动与Metropolis准则要理解模拟退火怎么工作得先拆解它的三个核心部件温度、新解产生扰动和接受准则。这三者共同构成了算法迭代的引擎。2.1 温度控制“探索野心”的调度器温度是SA算法中最重要的控制参数它直接决定了算法在多大程度上允许接受“坏解”。你可以把它理解为算法的“探索野心”或“随机性”。初始温度算法开始时温度设得很高。高温下系统几乎处于“熔融”状态算法会以很高的概率接受任何新产生的解无论它比当前解好还是差。这个阶段的目标是进行全局探索让解能在整个解空间里大幅跳跃避免过早地局限在某个小区域。降温策略随着迭代进行温度会按照某个退火进度表逐渐降低。最常用的是指数降温T_{k1} α * T_k其中α是一个接近1的常数比如0.95或0.99。温度降低意味着算法接受差解的概率越来越小其行为从“大胆探索”逐渐转向“精细打磨”。降温速度是关键太快了α太小如0.8容易淬火过快陷入局部最优太慢了α太大如0.999计算时间会变得无法忍受。终止温度当温度降低到一个预设的极小值比如1e-8时算法停止。此时系统几乎只接受更好的解算法行为类似于传统的局部搜索在当前位置附近进行微调。实操心得设置初始温度和降温系数是调参的艺术。一个经验法则是初始温度T0可以设置为使初始接受概率对差解的接受率在0.8左右的值。你可以通过少量试验来估算随机产生大量扰动计算目标函数差ΔE的均值然后根据Metropolis准则反推T0。降温系数α通常在[0.9, 0.999]之间对于复杂问题更慢的降温更大的α往往能得到更好的结果但需要权衡时间成本。2.2 新解产生如何在当前解的附近“蹦跶”如何从当前解x_old产生一个新解x_new这个过程称为“产生邻域解”或“扰动”。扰动的方式完全取决于你问题的具体形式。连续函数优化如果解是连续变量比如x [x1, x2, ..., xn]常见的扰动是在每个维度上加一个随机步长。例如x_new[i] x_old[i] (random() - 0.5) * step_size。这里step_size是步长通常可以与温度挂钩温度高时步长大便于大范围探索温度低时步长小便于局部精细搜索。组合优化如TSP如果解是一个排列如城市的访问顺序扰动操作可以是交换随机选择两个位置交换它们对应的城市。逆转随机选择一段子序列将其顺序反转。插入随机选择一个城市将其插入到另一个随机位置。背包问题可以随机翻转一个物品的“装入/不装入”状态。关键点扰动应该保证能够遍历整个解空间。也就是说从任何一个解出发经过有限次的扰动理论上可以到达解空间中的任何其他解。同时扰动不宜过大或过小需要与当前温度配合。2.3 Metropolis准则决定“跳不跳”的判官这是模拟退火最精髓的部分。产生新解x_new后计算目标函数值的变化ΔE E_new - E_old对于最小化问题E是目标函数值。如果ΔE 0说明新解更优那么无条件接受这个新解x_old x_new。如果ΔE 0说明新解更差。此时我们不是直接拒绝而是以一定的概率接受它。这个概率由Metropolis准则决定P exp(-ΔE / T)其中T是当前温度。这个公式非常美妙ΔE越大新解越差接受概率P越小。温度T越高接受概率P越大。高温时即使是很差的解也有不小概率被接受低温时接受差解的概率微乎其微。当T趋近于0时P也趋近于0算法退化为只接受更好解的爬山法。在实际程序中我们通过生成一个[0,1)区间的随机数rand与计算出的概率P比较。如果rand P则接受这个差解否则拒绝保留旧解。为什么这样有效接受差解给了算法跳出局部最优的能力。假设当前解在一个局部最优的“坑底”任何微小变动都会让目标函数变差ΔE0。传统的贪婪算法会困死在这里。但SA在温度还比较高的时候有机会接受一个“上山”的移动从而爬出这个坑有机会找到旁边更深的坑全局最优。3. 手把手实现一个经典的TSP问题求解程序理论说再多不如一行代码。我们用一个经典的旅行商问题作为例子来完整实现一个模拟退火算法。TSP问题描述很简单给定N个城市的坐标找一条访问每个城市恰好一次并回到起点的最短路径。3.1 问题定义与辅助函数首先我们定义城市和距离。为了简单我们在一个二维平面上随机生成城市坐标。import numpy as np import matplotlib.pyplot as plt import random import math import time # 设置随机种子确保结果可复现 np.random.seed(42) # 参数设置 num_cities 50 # 城市数量 initial_temperature 1000.0 # 初始温度 cooling_rate 0.995 # 降温系数 min_temperature 1e-8 # 终止温度 max_iterations_per_temp 100 # 每个温度下的迭代次数马尔可夫链长度 # 随机生成城市坐标 (x, y)范围在[0, 100]之间 cities np.random.rand(num_cities, 2) * 100 # 计算两个城市之间的欧氏距离 def distance(city1, city2): return np.sqrt(np.sum((city1 - city2) ** 2)) # 计算一条路径的总长度 def total_distance(path, cities): 计算给定路径顺序下的总旅行距离 total_dist 0.0 num_points len(path) for i in range(num_points): from_city cities[path[i]] # 取模运算实现闭环最后一个城市连接到第一个城市 to_city cities[path[(i 1) % num_points]] total_dist distance(from_city, to_city) return total_dist # 可视化函数 def plot_route(cities, path, title): 绘制城市和路径 plt.figure(figsize(10, 6)) # 绘制城市点 plt.scatter(cities[:, 0], cities[:, 1], cred, s50, zorder5) # 绘制路径线 for i in range(len(path)): start_idx path[i] end_idx path[(i 1) % len(path)] start_city cities[start_idx] end_city cities[end_idx] plt.plot([start_city[0], end_city[0]], [start_city[1], end_city[1]], b-, alpha0.6, linewidth1) # 标记起点 start_city cities[path[0]] plt.scatter(start_city[0], start_city[1], cgreen, s100, marker*, zorder10, labelStart) plt.title(title) plt.xlabel(X Coordinate) plt.ylabel(Y Coordinate) plt.legend() plt.grid(True, alpha0.3) plt.show()3.2 核心的模拟退火主循环接下来是算法的核心。我们从一个随机路径开始通过交换两个城市的位置来产生新解。def simulated_annealing_tsp(cities, initial_temp, cooling_rate, min_temp, iterations_per_temp): 模拟退火算法求解TSP 参数: cities: 城市坐标数组 initial_temp: 初始温度 cooling_rate: 降温系数 min_temp: 最低温度 iterations_per_temp: 每个温度下的迭代次数 返回: best_path: 找到的最佳路径 best_distance: 最佳路径长度 history: 迭代过程中的距离历史用于绘图分析 num_cities len(cities) # 1. 初始化生成一个随机路径作为当前解 current_path list(range(num_cities)) random.shuffle(current_path) # 随机打乱顺序 current_distance total_distance(current_path, cities) # 记录最优解 best_path current_path.copy() best_distance current_distance # 记录历史用于分析收敛过程 distance_history [current_distance] temperature_history [initial_temp] # 2. 开始退火过程 temperature initial_temp iteration 0 while temperature min_temp: for _ in range(iterations_per_temp): # 2.1 产生新解通过交换两个随机城市的位置 new_path current_path.copy() # 随机选择两个不同的索引 i, j random.sample(range(num_cities), 2) # 交换城市 new_path[i], new_path[j] new_path[j], new_path[i] # 计算新路径的距离 new_distance total_distance(new_path, cities) # 计算目标函数差 (我们是最小化距离) delta_distance new_distance - current_distance # 2.2 Metropolis准则判断是否接受新解 if delta_distance 0: # 新解更优直接接受 accept True else: # 新解更差以一定概率接受 acceptance_probability math.exp(-delta_distance / temperature) accept (random.random() acceptance_probability) # 2.3 更新当前解 if accept: current_path new_path current_distance new_distance # 更新全局最优解 if current_distance best_distance: best_path current_path.copy() best_distance current_distance iteration 1 distance_history.append(current_distance) temperature_history.append(temperature) # 2.4 降温 temperature * cooling_rate # 可选打印进度对于大规模问题可以每N次降温打印一次 if int(temperature * 100) % 10 0: # 粗略控制打印频率 print(fIteration: {iteration}, Temp: {temperature:.4f}, Current Dist: {current_distance:.2f}, Best Dist: {best_distance:.2f}) print(fSA completed. Total iterations: {iteration}) print(fBest distance found: {best_distance:.2f}) return best_path, best_distance, distance_history, temperature_history3.3 运行与结果分析现在让我们运行这个程序并看看效果。# 运行模拟退火算法 start_time time.time() best_path, best_distance, dist_history, temp_history simulated_annealing_tsp( cities, initial_temperature, cooling_rate, min_temperature, max_iterations_per_temp ) end_time time.time() print(f\n算法运行时间: {end_time - start_time:.2f} 秒) # 绘制最优路径 plot_route(cities, best_path, fBest TSP Route Found by SA\nTotal Distance: {best_distance:.2f}) # 绘制优化过程收敛曲线 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(dist_history, linewidth0.5) plt.title(Distance Convergence History) plt.xlabel(Iteration) plt.ylabel(Total Distance) plt.grid(True, alpha0.3) plt.subplot(1, 2, 2) plt.semilogy(temp_history, linewidth0.5) # 温度对数坐标显示更清晰 plt.title(Temperature Schedule) plt.xlabel(Iteration) plt.ylabel(Temperature (log scale)) plt.grid(True, alpha0.3) plt.tight_layout() plt.show() # 打印前10个城市的访问顺序作为示例 print(\nBest path (first 10 cities index):, best_path[:10])运行结果解读 当你运行这段代码你会看到控制台输出迭代过程最终得到一个比初始随机路径短得多的总距离。收敛曲线图会显示距离值在初期剧烈波动高温阶段大量接受差解随后波动逐渐减小最终趋于稳定低温阶段精细搜索。路径图则会直观展示出一条相对合理的、交叉较少的环游路线。避坑经验初始解的影响虽然SA理论上对初始解不敏感因为高温阶段会打乱它但一个好的初始解如最近邻法生成的路径可以显著加快收敛速度。对于时间敏感的应用值得花点心思初始化。马尔可夫链长度参数max_iterations_per_temp每个温度下的迭代次数很重要。太短系统在每个温度下来不及达到平衡状态太长计算开销大。一个实用技巧是让它与问题规模相关比如设为城市数量的若干倍如10*N或者采用自适应策略当连续拒绝的解达到一定数量时就提前进入降温步骤。邻域结构设计本例使用了最简单的“交换”操作。对于TSP逆转(2-opt)操作通常更有效它能直接消除路径交叉。你可以尝试实现2-opt作为扰动方式对比效果。4. 调参与改进让SA算法更高效实用一个基础的SA框架跑起来后你会发现效果可能时好时坏非常依赖参数设置。这部分就是算法工程师的“内功”了。下面分享几个关键的调参与改进方向。4.1 参数调优寻找属于你问题的“退火进度表”SA有四个核心参数初始温度T0、终止温度T_min、降温系数α、马尔可夫链长度L。没有放之四海而皆准的最优值必须针对具体问题调整。自适应初始温度与其拍脑袋定一个T01000不如让算法自己决定。可以这样做进行若干次比如1000次随机扰动计算所有ΔE的平均值ΔE然后根据一个期望的初始接受概率P0比如0.8来反推T0 -ΔE / ln(P0)。这样得到的T0更能匹配当前问题的能量尺度。自适应降温可以根据搜索过程动态调整降温速度。例如如果当前温度下接受率很高0.6说明温度还太高可以加快降温用更小的α如果接受率很低0.2说明降温太快了系统可能被“冻住”在次优解此时应该减慢降温速度用更大的α甚至短暂“回温”。终止条件多样化除了温度低于阈值还可以结合其他条件连续若干个温度下最优解都没有改进。达到最大迭代次数或最大运行时间。目标函数值已经达到一个可接受的阈值。4.2 邻域操作的进阶设计扰动操作的设计极大影响算法性能。好的邻域操作应该在“扰动强度”和“搜索效率”之间取得平衡。对于TSP的2-opt操作前面提到的交换操作是点操作而2-opt是边操作。它随机选择两条不相邻的边(i, i1)和(j, j1)然后逆转i1到j之间的子路径。这种操作能直接解开路径交叉是TSP求解器的标配。实现2-opt后你会发现收敛速度和最终解质量通常比简单交换好很多。混合邻域不要只使用一种扰动。可以以一定概率选择不同的邻域操作。例如80%的概率使用细粒度的“交换两个城市”20%的概率使用粗粒度的“逆转一段长路径”或“随机插入一个城市到新位置”。这有助于平衡局部搜索和全局探索。问题特定的启发式将领域知识融入扰动。例如在车辆路径问题中可以设计“将一条路线上的一个客户移到另一条路线”的扰动。这种有导向的扰动比完全随机的扰动更有效。4.3 记忆与重启避免“失忆”和“早熟”基础的SA算法只保留一个“当前解”和一个“历史最优解”。我们可以增加一些机制来提升鲁棒性。模拟回火不仅记录历史最优解还记录历史上遇到过的其他优秀解。当算法陷入停滞时可以从这些优秀解中随机选择一个作为新的当前解并适当提高温度重新开始搜索。这相当于给了算法多次“重启”的机会。带记忆的SA在每次接受新解时不仅更新当前解还以一定概率将新解存入一个“精英解池”。这个池子的大小有限采用类似LRU最近最少使用的机制进行更新。在降温后期可以从池中选取最优解作为最终输出这比单纯依赖最后一次迭代的当前解更可靠。4.4 与其他算法结合混合策略SA不排斥和其他优化方法联用往往能产生“112”的效果。SA 局部搜索在SA的每次迭代中或者当温度降到较低水平时对当前解执行一次快速的局部搜索如对于TSP进行几次贪婪的2-opt优化。这能快速提升解的质量相当于在SA的全局探索框架内嵌入了局部挖掘的能力。这种算法常被称为“模拟退火局部搜索”。SA初始化其他算法用SA快速得到一个质量不错的解然后将这个解作为更复杂算法如遗传算法、蚁群算法的初始种群之一。SA的全局探索能力可以帮助跳出局部最优为后续算法提供一个更好的起点。个人体会调参过程很像在带一个不确定性的优化系统。我的习惯是先跑一个参数范围较广的快速实验比如α从0.8到0.999T0从10到10000观察收敛曲线和最终解的质量锁定一个大概的“好区域”。然后在这个区域里进行更精细的网格搜索或随机搜索。一定要把每次实验的参数、结果包括最终解和收敛过程图记录下来久而久之你对不同规模、不同类型问题的参数敏感度就会有直觉。5. 从TSP到通用框架如何将SA应用于你的问题TSP只是一个例子。模拟退火的威力在于其框架的通用性。只要你能够定义以下三要素就可以用SA来求解你的优化问题解的表达如何用一个数据结构数组、列表、字典、树等表示你的一个候选解。目标函数如何计算一个解的好坏需要最小化或最大化的值。邻域操作如何从一个给定解通过一个小的、随机的变动产生一个新的、“相邻的”解。让我们看两个其他领域的例子体会一下这种转换。5.1 案例一函数优化寻找全局最小值假设我们想找Rastrigin函数在定义域内的全局最小值。这是一个著名的多峰测试函数有很多局部极小点。import numpy as np import matplotlib.pyplot as plt # Rastrigin 函数 def rastrigin(x): x可以是一个多维向量 A 10 n len(x) return A * n np.sum(x**2 - A * np.cos(2 * np.pi * x)) # SA求解器通用框架 def simulated_annealing_continuous(obj_func, bounds, initial_temp, cooling_rate, min_temp, iterations_per_temp): 用于连续函数优化的SA框架 参数: obj_func: 目标函数 bounds: 每个变量的上下界列表如 [(lb1, ub1), (lb2, ub2), ...] ... 其他参数同前 dim len(bounds) # 初始化当前解在边界内随机生成 current_solution np.array([np.random.uniform(low, high) for low, high in bounds]) current_value obj_func(current_solution) best_solution current_solution.copy() best_value current_value temperature initial_temp while temperature min_temp: for _ in range(iterations_per_temp): # 产生新解在当前解基础上添加随机扰动 # 扰动幅度可以与温度相关实现变步长搜索 step_size temperature / initial_temp # 一种简单的步长缩放策略 new_solution current_solution step_size * np.random.randn(dim) * (np.array([ub - lb for lb, ub in bounds])) # 确保新解在边界内越界处理反射或裁剪 for i in range(dim): lb, ub bounds[i] if new_solution[i] lb: new_solution[i] lb (lb - new_solution[i]) # 反射 if new_solution[i] ub: new_solution[i] ub # 二次越界则裁剪 elif new_solution[i] ub: new_solution[i] ub - (new_solution[i] - ub) if new_solution[i] lb: new_solution[i] lb new_value obj_func(new_solution) delta new_value - current_value if delta 0 or np.random.rand() np.exp(-delta / temperature): current_solution new_solution current_value new_value if current_value best_value: best_solution current_solution.copy() best_value current_value temperature * cooling_rate return best_solution, best_value # 定义搜索边界假设是2维问题 bounds [(-5.12, 5.12), (-5.12, 5.12)] best_sol, best_val simulated_annealing_continuous(rastrigin, bounds, initial_temperature100, cooling_rate0.99, min_temperature1e-8, iterations_per_temp100) print(f找到的最优解: {best_sol}) print(f最优值: {best_val:.6f}) print(f理论全局最优解: [0, 0], 最优值: 0)这个例子展示了如何将SA应用于连续空间。关键在于邻域操作的设计这里用了高斯随机扰动并且让步长随着温度降低而减小实现了从粗搜索到细搜索的过渡。5.2 案例二0-1背包问题背包问题给定一组物品每个物品有重量和价值在背包容量限制下选择物品使得总价值最大。def knapsack_sa(weights, values, capacity, initial_temp, cooling_rate, min_temp, iterations): 用SA解0-1背包问题 解的表达: 一个二进制列表1表示选中0表示不选。 邻域操作: 随机翻转一位即改变一个物品的选择状态如果翻转为1导致超重则拒绝该扰动。 num_items len(weights) # 辅助函数计算一个解的价值和重量 def evaluate(solution): total_value np.dot(solution, values) total_weight np.dot(solution, weights) if total_weight capacity: return -float(inf), total_weight # 不可行解价值设为负无穷 return total_value, total_weight # 生成初始可行解贪心法按价值密度排序直到装满 density [v/w for v, w in zip(values, weights)] sorted_idx np.argsort(density)[::-1] current_sol np.zeros(num_items, dtypeint) current_weight 0 for idx in sorted_idx: if current_weight weights[idx] capacity: current_sol[idx] 1 current_weight weights[idx] current_value, _ evaluate(current_sol) best_sol current_sol.copy() best_value current_value temperature initial_temp while temperature min_temp: for _ in range(iterations): # 产生新解随机翻转一位 new_sol current_sol.copy() flip_idx np.random.randint(num_items) new_sol[flip_idx] 1 - new_sol[flip_idx] # 0变11变0 new_value, new_weight evaluate(new_sol) # 如果新解不可行超重直接拒绝 if new_value -float(inf): continue delta current_value - new_value # 注意背包问题是最大化所以delta是旧值减新值 # Metropolis准则 (对于最大化问题接受更差解的概率是 exp(delta/T)) if delta 0 or np.random.rand() np.exp(delta / temperature): current_sol new_sol current_value new_value if current_value best_value: best_sol current_sol.copy() best_value current_value temperature * cooling_rate return best_sol, best_value # 示例数据 weights [2, 3, 4, 5, 9] values [3, 4, 5, 8, 10] capacity 20 best_solution, best_val knapsack_sa(weights, values, capacity, initial_temp50, cooling_rate0.95, min_temp1e-6, iterations200) print(f最优选择物品索引: {np.where(best_solution1)[0]}) print(f总价值: {best_val}) print(f总重量: {np.dot(best_solution, weights)})这个例子展示了如何将SA应用于离散组合优化。解的表达是二进制串邻域操作是单点翻转。注意我们对不可行解超重的处理直接拒绝。这是一种简单的约束处理方法。对于更复杂的约束可能需要设计专门的修复算子或者在目标函数中加入惩罚项。通用框架总结无论你的问题是什么套用SA的步骤都是一样的1) 编码你的解2) 定义如何评价解目标函数3) 定义如何产生一个相似的“邻居”解4) 设定退火进度表5) 跑起来然后耐心调参。SA的魅力就在于只要你能定义这三点它就能为你工作虽然不一定总能找到绝对最优但在很多复杂场景下它提供的是一个在合理时间内找到高质量近似解的可靠途径。 模拟退火算法：从物理退火到全局优化，附Python实现与调参指南
1. 从“烧铁淬火”到全局寻优模拟退火算法初印象如果你在数学建模、算法竞赛或者优化问题的实际工程中摸爬滚打过一阵子大概率会听过“模拟退火”这个名字。我第一次接触它是在为一个物流中心的选址问题挠头时。当时试遍了梯度下降、遗传算法效果总是不尽如人意要么陷入局部最优解出不来要么收敛速度慢得让人心焦。直到一位前辈扔过来一句“试试模拟退火吧这玩意儿有时候能给你惊喜。” 我抱着试试看的心态去研究结果发现这个算法的思想之巧妙简直是把物理世界的智慧用到了数学优化里。模拟退火算法英文是Simulated Annealing简称SA。它的核心灵感来源于冶金学中的“退火”工艺。想象一下铁匠打铁先把铁块加热到高温这时铁内部的原子处于高能、活跃的无序状态然后缓慢地、有控制地降温原子会在这个过程中逐渐找到一个能量最低、结构最稳定的排列方式从而得到一块坚韧的好钢。SA算法就是把我们要优化的目标函数比如成本、距离、误差想象成这个系统的“能量”通过模拟“加热”和“缓慢降温”的过程让解在解空间中“跳跃”有策略地接受一些暂时看起来更差的解从而有机会跳出局部最优的陷阱最终逼近全局最优解。这听起来有点反直觉对吧我们找最优解不就应该一路向下只接受更好的解吗但现实世界中的很多优化问题其解空间就像一片连绵起伏的山脉布满了深坑局部最优和高峰。如果你只肯往下走那很容易掉进一个离你起点最近的深坑里再也看不到远处那个更深的峡谷全局最优。模拟退火给你的就是一把“登山杖”允许你偶尔往上爬几步去看看山那边的风景。它特别适合解决那些解空间复杂、存在大量局部最优解的组合优化问题比如旅行商问题、调度问题、布局问题以及函数优化等。在数学建模中当你的模型目标函数复杂、非线性、不可微或者约束条件难以处理时SA往往是一个值得考虑的“暴力美学”式工具。2. 算法核心温度、扰动与Metropolis准则要理解模拟退火怎么工作得先拆解它的三个核心部件温度、新解产生扰动和接受准则。这三者共同构成了算法迭代的引擎。2.1 温度控制“探索野心”的调度器温度是SA算法中最重要的控制参数它直接决定了算法在多大程度上允许接受“坏解”。你可以把它理解为算法的“探索野心”或“随机性”。初始温度算法开始时温度设得很高。高温下系统几乎处于“熔融”状态算法会以很高的概率接受任何新产生的解无论它比当前解好还是差。这个阶段的目标是进行全局探索让解能在整个解空间里大幅跳跃避免过早地局限在某个小区域。降温策略随着迭代进行温度会按照某个退火进度表逐渐降低。最常用的是指数降温T_{k1} α * T_k其中α是一个接近1的常数比如0.95或0.99。温度降低意味着算法接受差解的概率越来越小其行为从“大胆探索”逐渐转向“精细打磨”。降温速度是关键太快了α太小如0.8容易淬火过快陷入局部最优太慢了α太大如0.999计算时间会变得无法忍受。终止温度当温度降低到一个预设的极小值比如1e-8时算法停止。此时系统几乎只接受更好的解算法行为类似于传统的局部搜索在当前位置附近进行微调。实操心得设置初始温度和降温系数是调参的艺术。一个经验法则是初始温度T0可以设置为使初始接受概率对差解的接受率在0.8左右的值。你可以通过少量试验来估算随机产生大量扰动计算目标函数差ΔE的均值然后根据Metropolis准则反推T0。降温系数α通常在[0.9, 0.999]之间对于复杂问题更慢的降温更大的α往往能得到更好的结果但需要权衡时间成本。2.2 新解产生如何在当前解的附近“蹦跶”如何从当前解x_old产生一个新解x_new这个过程称为“产生邻域解”或“扰动”。扰动的方式完全取决于你问题的具体形式。连续函数优化如果解是连续变量比如x [x1, x2, ..., xn]常见的扰动是在每个维度上加一个随机步长。例如x_new[i] x_old[i] (random() - 0.5) * step_size。这里step_size是步长通常可以与温度挂钩温度高时步长大便于大范围探索温度低时步长小便于局部精细搜索。组合优化如TSP如果解是一个排列如城市的访问顺序扰动操作可以是交换随机选择两个位置交换它们对应的城市。逆转随机选择一段子序列将其顺序反转。插入随机选择一个城市将其插入到另一个随机位置。背包问题可以随机翻转一个物品的“装入/不装入”状态。关键点扰动应该保证能够遍历整个解空间。也就是说从任何一个解出发经过有限次的扰动理论上可以到达解空间中的任何其他解。同时扰动不宜过大或过小需要与当前温度配合。2.3 Metropolis准则决定“跳不跳”的判官这是模拟退火最精髓的部分。产生新解x_new后计算目标函数值的变化ΔE E_new - E_old对于最小化问题E是目标函数值。如果ΔE 0说明新解更优那么无条件接受这个新解x_old x_new。如果ΔE 0说明新解更差。此时我们不是直接拒绝而是以一定的概率接受它。这个概率由Metropolis准则决定P exp(-ΔE / T)其中T是当前温度。这个公式非常美妙ΔE越大新解越差接受概率P越小。温度T越高接受概率P越大。高温时即使是很差的解也有不小概率被接受低温时接受差解的概率微乎其微。当T趋近于0时P也趋近于0算法退化为只接受更好解的爬山法。在实际程序中我们通过生成一个[0,1)区间的随机数rand与计算出的概率P比较。如果rand P则接受这个差解否则拒绝保留旧解。为什么这样有效接受差解给了算法跳出局部最优的能力。假设当前解在一个局部最优的“坑底”任何微小变动都会让目标函数变差ΔE0。传统的贪婪算法会困死在这里。但SA在温度还比较高的时候有机会接受一个“上山”的移动从而爬出这个坑有机会找到旁边更深的坑全局最优。3. 手把手实现一个经典的TSP问题求解程序理论说再多不如一行代码。我们用一个经典的旅行商问题作为例子来完整实现一个模拟退火算法。TSP问题描述很简单给定N个城市的坐标找一条访问每个城市恰好一次并回到起点的最短路径。3.1 问题定义与辅助函数首先我们定义城市和距离。为了简单我们在一个二维平面上随机生成城市坐标。import numpy as np import matplotlib.pyplot as plt import random import math import time # 设置随机种子确保结果可复现 np.random.seed(42) # 参数设置 num_cities 50 # 城市数量 initial_temperature 1000.0 # 初始温度 cooling_rate 0.995 # 降温系数 min_temperature 1e-8 # 终止温度 max_iterations_per_temp 100 # 每个温度下的迭代次数马尔可夫链长度 # 随机生成城市坐标 (x, y)范围在[0, 100]之间 cities np.random.rand(num_cities, 2) * 100 # 计算两个城市之间的欧氏距离 def distance(city1, city2): return np.sqrt(np.sum((city1 - city2) ** 2)) # 计算一条路径的总长度 def total_distance(path, cities): 计算给定路径顺序下的总旅行距离 total_dist 0.0 num_points len(path) for i in range(num_points): from_city cities[path[i]] # 取模运算实现闭环最后一个城市连接到第一个城市 to_city cities[path[(i 1) % num_points]] total_dist distance(from_city, to_city) return total_dist # 可视化函数 def plot_route(cities, path, title): 绘制城市和路径 plt.figure(figsize(10, 6)) # 绘制城市点 plt.scatter(cities[:, 0], cities[:, 1], cred, s50, zorder5) # 绘制路径线 for i in range(len(path)): start_idx path[i] end_idx path[(i 1) % len(path)] start_city cities[start_idx] end_city cities[end_idx] plt.plot([start_city[0], end_city[0]], [start_city[1], end_city[1]], b-, alpha0.6, linewidth1) # 标记起点 start_city cities[path[0]] plt.scatter(start_city[0], start_city[1], cgreen, s100, marker*, zorder10, labelStart) plt.title(title) plt.xlabel(X Coordinate) plt.ylabel(Y Coordinate) plt.legend() plt.grid(True, alpha0.3) plt.show()3.2 核心的模拟退火主循环接下来是算法的核心。我们从一个随机路径开始通过交换两个城市的位置来产生新解。def simulated_annealing_tsp(cities, initial_temp, cooling_rate, min_temp, iterations_per_temp): 模拟退火算法求解TSP 参数: cities: 城市坐标数组 initial_temp: 初始温度 cooling_rate: 降温系数 min_temp: 最低温度 iterations_per_temp: 每个温度下的迭代次数 返回: best_path: 找到的最佳路径 best_distance: 最佳路径长度 history: 迭代过程中的距离历史用于绘图分析 num_cities len(cities) # 1. 初始化生成一个随机路径作为当前解 current_path list(range(num_cities)) random.shuffle(current_path) # 随机打乱顺序 current_distance total_distance(current_path, cities) # 记录最优解 best_path current_path.copy() best_distance current_distance # 记录历史用于分析收敛过程 distance_history [current_distance] temperature_history [initial_temp] # 2. 开始退火过程 temperature initial_temp iteration 0 while temperature min_temp: for _ in range(iterations_per_temp): # 2.1 产生新解通过交换两个随机城市的位置 new_path current_path.copy() # 随机选择两个不同的索引 i, j random.sample(range(num_cities), 2) # 交换城市 new_path[i], new_path[j] new_path[j], new_path[i] # 计算新路径的距离 new_distance total_distance(new_path, cities) # 计算目标函数差 (我们是最小化距离) delta_distance new_distance - current_distance # 2.2 Metropolis准则判断是否接受新解 if delta_distance 0: # 新解更优直接接受 accept True else: # 新解更差以一定概率接受 acceptance_probability math.exp(-delta_distance / temperature) accept (random.random() acceptance_probability) # 2.3 更新当前解 if accept: current_path new_path current_distance new_distance # 更新全局最优解 if current_distance best_distance: best_path current_path.copy() best_distance current_distance iteration 1 distance_history.append(current_distance) temperature_history.append(temperature) # 2.4 降温 temperature * cooling_rate # 可选打印进度对于大规模问题可以每N次降温打印一次 if int(temperature * 100) % 10 0: # 粗略控制打印频率 print(fIteration: {iteration}, Temp: {temperature:.4f}, Current Dist: {current_distance:.2f}, Best Dist: {best_distance:.2f}) print(fSA completed. Total iterations: {iteration}) print(fBest distance found: {best_distance:.2f}) return best_path, best_distance, distance_history, temperature_history3.3 运行与结果分析现在让我们运行这个程序并看看效果。# 运行模拟退火算法 start_time time.time() best_path, best_distance, dist_history, temp_history simulated_annealing_tsp( cities, initial_temperature, cooling_rate, min_temperature, max_iterations_per_temp ) end_time time.time() print(f\n算法运行时间: {end_time - start_time:.2f} 秒) # 绘制最优路径 plot_route(cities, best_path, fBest TSP Route Found by SA\nTotal Distance: {best_distance:.2f}) # 绘制优化过程收敛曲线 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(dist_history, linewidth0.5) plt.title(Distance Convergence History) plt.xlabel(Iteration) plt.ylabel(Total Distance) plt.grid(True, alpha0.3) plt.subplot(1, 2, 2) plt.semilogy(temp_history, linewidth0.5) # 温度对数坐标显示更清晰 plt.title(Temperature Schedule) plt.xlabel(Iteration) plt.ylabel(Temperature (log scale)) plt.grid(True, alpha0.3) plt.tight_layout() plt.show() # 打印前10个城市的访问顺序作为示例 print(\nBest path (first 10 cities index):, best_path[:10])运行结果解读 当你运行这段代码你会看到控制台输出迭代过程最终得到一个比初始随机路径短得多的总距离。收敛曲线图会显示距离值在初期剧烈波动高温阶段大量接受差解随后波动逐渐减小最终趋于稳定低温阶段精细搜索。路径图则会直观展示出一条相对合理的、交叉较少的环游路线。避坑经验初始解的影响虽然SA理论上对初始解不敏感因为高温阶段会打乱它但一个好的初始解如最近邻法生成的路径可以显著加快收敛速度。对于时间敏感的应用值得花点心思初始化。马尔可夫链长度参数max_iterations_per_temp每个温度下的迭代次数很重要。太短系统在每个温度下来不及达到平衡状态太长计算开销大。一个实用技巧是让它与问题规模相关比如设为城市数量的若干倍如10*N或者采用自适应策略当连续拒绝的解达到一定数量时就提前进入降温步骤。邻域结构设计本例使用了最简单的“交换”操作。对于TSP逆转(2-opt)操作通常更有效它能直接消除路径交叉。你可以尝试实现2-opt作为扰动方式对比效果。4. 调参与改进让SA算法更高效实用一个基础的SA框架跑起来后你会发现效果可能时好时坏非常依赖参数设置。这部分就是算法工程师的“内功”了。下面分享几个关键的调参与改进方向。4.1 参数调优寻找属于你问题的“退火进度表”SA有四个核心参数初始温度T0、终止温度T_min、降温系数α、马尔可夫链长度L。没有放之四海而皆准的最优值必须针对具体问题调整。自适应初始温度与其拍脑袋定一个T01000不如让算法自己决定。可以这样做进行若干次比如1000次随机扰动计算所有ΔE的平均值ΔE然后根据一个期望的初始接受概率P0比如0.8来反推T0 -ΔE / ln(P0)。这样得到的T0更能匹配当前问题的能量尺度。自适应降温可以根据搜索过程动态调整降温速度。例如如果当前温度下接受率很高0.6说明温度还太高可以加快降温用更小的α如果接受率很低0.2说明降温太快了系统可能被“冻住”在次优解此时应该减慢降温速度用更大的α甚至短暂“回温”。终止条件多样化除了温度低于阈值还可以结合其他条件连续若干个温度下最优解都没有改进。达到最大迭代次数或最大运行时间。目标函数值已经达到一个可接受的阈值。4.2 邻域操作的进阶设计扰动操作的设计极大影响算法性能。好的邻域操作应该在“扰动强度”和“搜索效率”之间取得平衡。对于TSP的2-opt操作前面提到的交换操作是点操作而2-opt是边操作。它随机选择两条不相邻的边(i, i1)和(j, j1)然后逆转i1到j之间的子路径。这种操作能直接解开路径交叉是TSP求解器的标配。实现2-opt后你会发现收敛速度和最终解质量通常比简单交换好很多。混合邻域不要只使用一种扰动。可以以一定概率选择不同的邻域操作。例如80%的概率使用细粒度的“交换两个城市”20%的概率使用粗粒度的“逆转一段长路径”或“随机插入一个城市到新位置”。这有助于平衡局部搜索和全局探索。问题特定的启发式将领域知识融入扰动。例如在车辆路径问题中可以设计“将一条路线上的一个客户移到另一条路线”的扰动。这种有导向的扰动比完全随机的扰动更有效。4.3 记忆与重启避免“失忆”和“早熟”基础的SA算法只保留一个“当前解”和一个“历史最优解”。我们可以增加一些机制来提升鲁棒性。模拟回火不仅记录历史最优解还记录历史上遇到过的其他优秀解。当算法陷入停滞时可以从这些优秀解中随机选择一个作为新的当前解并适当提高温度重新开始搜索。这相当于给了算法多次“重启”的机会。带记忆的SA在每次接受新解时不仅更新当前解还以一定概率将新解存入一个“精英解池”。这个池子的大小有限采用类似LRU最近最少使用的机制进行更新。在降温后期可以从池中选取最优解作为最终输出这比单纯依赖最后一次迭代的当前解更可靠。4.4 与其他算法结合混合策略SA不排斥和其他优化方法联用往往能产生“112”的效果。SA 局部搜索在SA的每次迭代中或者当温度降到较低水平时对当前解执行一次快速的局部搜索如对于TSP进行几次贪婪的2-opt优化。这能快速提升解的质量相当于在SA的全局探索框架内嵌入了局部挖掘的能力。这种算法常被称为“模拟退火局部搜索”。SA初始化其他算法用SA快速得到一个质量不错的解然后将这个解作为更复杂算法如遗传算法、蚁群算法的初始种群之一。SA的全局探索能力可以帮助跳出局部最优为后续算法提供一个更好的起点。个人体会调参过程很像在带一个不确定性的优化系统。我的习惯是先跑一个参数范围较广的快速实验比如α从0.8到0.999T0从10到10000观察收敛曲线和最终解的质量锁定一个大概的“好区域”。然后在这个区域里进行更精细的网格搜索或随机搜索。一定要把每次实验的参数、结果包括最终解和收敛过程图记录下来久而久之你对不同规模、不同类型问题的参数敏感度就会有直觉。5. 从TSP到通用框架如何将SA应用于你的问题TSP只是一个例子。模拟退火的威力在于其框架的通用性。只要你能够定义以下三要素就可以用SA来求解你的优化问题解的表达如何用一个数据结构数组、列表、字典、树等表示你的一个候选解。目标函数如何计算一个解的好坏需要最小化或最大化的值。邻域操作如何从一个给定解通过一个小的、随机的变动产生一个新的、“相邻的”解。让我们看两个其他领域的例子体会一下这种转换。5.1 案例一函数优化寻找全局最小值假设我们想找Rastrigin函数在定义域内的全局最小值。这是一个著名的多峰测试函数有很多局部极小点。import numpy as np import matplotlib.pyplot as plt # Rastrigin 函数 def rastrigin(x): x可以是一个多维向量 A 10 n len(x) return A * n np.sum(x**2 - A * np.cos(2 * np.pi * x)) # SA求解器通用框架 def simulated_annealing_continuous(obj_func, bounds, initial_temp, cooling_rate, min_temp, iterations_per_temp): 用于连续函数优化的SA框架 参数: obj_func: 目标函数 bounds: 每个变量的上下界列表如 [(lb1, ub1), (lb2, ub2), ...] ... 其他参数同前 dim len(bounds) # 初始化当前解在边界内随机生成 current_solution np.array([np.random.uniform(low, high) for low, high in bounds]) current_value obj_func(current_solution) best_solution current_solution.copy() best_value current_value temperature initial_temp while temperature min_temp: for _ in range(iterations_per_temp): # 产生新解在当前解基础上添加随机扰动 # 扰动幅度可以与温度相关实现变步长搜索 step_size temperature / initial_temp # 一种简单的步长缩放策略 new_solution current_solution step_size * np.random.randn(dim) * (np.array([ub - lb for lb, ub in bounds])) # 确保新解在边界内越界处理反射或裁剪 for i in range(dim): lb, ub bounds[i] if new_solution[i] lb: new_solution[i] lb (lb - new_solution[i]) # 反射 if new_solution[i] ub: new_solution[i] ub # 二次越界则裁剪 elif new_solution[i] ub: new_solution[i] ub - (new_solution[i] - ub) if new_solution[i] lb: new_solution[i] lb new_value obj_func(new_solution) delta new_value - current_value if delta 0 or np.random.rand() np.exp(-delta / temperature): current_solution new_solution current_value new_value if current_value best_value: best_solution current_solution.copy() best_value current_value temperature * cooling_rate return best_solution, best_value # 定义搜索边界假设是2维问题 bounds [(-5.12, 5.12), (-5.12, 5.12)] best_sol, best_val simulated_annealing_continuous(rastrigin, bounds, initial_temperature100, cooling_rate0.99, min_temperature1e-8, iterations_per_temp100) print(f找到的最优解: {best_sol}) print(f最优值: {best_val:.6f}) print(f理论全局最优解: [0, 0], 最优值: 0)这个例子展示了如何将SA应用于连续空间。关键在于邻域操作的设计这里用了高斯随机扰动并且让步长随着温度降低而减小实现了从粗搜索到细搜索的过渡。5.2 案例二0-1背包问题背包问题给定一组物品每个物品有重量和价值在背包容量限制下选择物品使得总价值最大。def knapsack_sa(weights, values, capacity, initial_temp, cooling_rate, min_temp, iterations): 用SA解0-1背包问题 解的表达: 一个二进制列表1表示选中0表示不选。 邻域操作: 随机翻转一位即改变一个物品的选择状态如果翻转为1导致超重则拒绝该扰动。 num_items len(weights) # 辅助函数计算一个解的价值和重量 def evaluate(solution): total_value np.dot(solution, values) total_weight np.dot(solution, weights) if total_weight capacity: return -float(inf), total_weight # 不可行解价值设为负无穷 return total_value, total_weight # 生成初始可行解贪心法按价值密度排序直到装满 density [v/w for v, w in zip(values, weights)] sorted_idx np.argsort(density)[::-1] current_sol np.zeros(num_items, dtypeint) current_weight 0 for idx in sorted_idx: if current_weight weights[idx] capacity: current_sol[idx] 1 current_weight weights[idx] current_value, _ evaluate(current_sol) best_sol current_sol.copy() best_value current_value temperature initial_temp while temperature min_temp: for _ in range(iterations): # 产生新解随机翻转一位 new_sol current_sol.copy() flip_idx np.random.randint(num_items) new_sol[flip_idx] 1 - new_sol[flip_idx] # 0变11变0 new_value, new_weight evaluate(new_sol) # 如果新解不可行超重直接拒绝 if new_value -float(inf): continue delta current_value - new_value # 注意背包问题是最大化所以delta是旧值减新值 # Metropolis准则 (对于最大化问题接受更差解的概率是 exp(delta/T)) if delta 0 or np.random.rand() np.exp(delta / temperature): current_sol new_sol current_value new_value if current_value best_value: best_sol current_sol.copy() best_value current_value temperature * cooling_rate return best_sol, best_value # 示例数据 weights [2, 3, 4, 5, 9] values [3, 4, 5, 8, 10] capacity 20 best_solution, best_val knapsack_sa(weights, values, capacity, initial_temp50, cooling_rate0.95, min_temp1e-6, iterations200) print(f最优选择物品索引: {np.where(best_solution1)[0]}) print(f总价值: {best_val}) print(f总重量: {np.dot(best_solution, weights)})这个例子展示了如何将SA应用于离散组合优化。解的表达是二进制串邻域操作是单点翻转。注意我们对不可行解超重的处理直接拒绝。这是一种简单的约束处理方法。对于更复杂的约束可能需要设计专门的修复算子或者在目标函数中加入惩罚项。通用框架总结无论你的问题是什么套用SA的步骤都是一样的1) 编码你的解2) 定义如何评价解目标函数3) 定义如何产生一个相似的“邻居”解4) 设定退火进度表5) 跑起来然后耐心调参。SA的魅力就在于只要你能定义这三点它就能为你工作虽然不一定总能找到绝对最优但在很多复杂场景下它提供的是一个在合理时间内找到高质量近似解的可靠途径。