用Python玩转模拟退火算法：从物理退火到TSP求解的保姆级实战

Python实现模拟退火算法：从物理原理到旅行商问题实战

想象一下金属工匠如何打造一把完美的剑——他们将金属加热至炽热状态，然后缓慢冷却，让原子在降温过程中自然排列成最稳定的结构。这种古老的工艺启发了计算机科学家开发出模拟退火算法，一种能优雅跳出局部最优陷阱的全局优化方法。本文将带您从物理退火的直观理解出发，手把手实现解决旅行商问题(TSP)的完整Python方案。

1. 物理退火与算法思想的完美映射

固体退火过程中，温度的高低直接影响着原子的活跃程度。当金属被加热到极高温度时：

• 高温阶段：原子剧烈运动，几乎不受束缚
• 降温过程：原子逐渐找到能量更低的位置
• 最终状态：形成规则的晶体结构，达到最小内能

模拟退火算法将这一物理过程抽象为三个核心要素：

# 算法伪代码框架 def simulated_annealing(): current_solution = random_initial_solution() current_energy = calculate_energy(current_solution) T = initial_temperature while T > final_temperature: new_solution = perturb(current_solution) new_energy = calculate_energy(new_solution) if acceptance_probability(current_energy, new_energy, T) > random(): current_solution = new_solution current_energy = new_energy T = cool(T)

关键理解：算法在高温时更愿意接受劣解（探索全局），随着温度降低逐渐收敛（利用局部）

2. TSP问题建模与解表示

旅行商问题需要找到访问所有城市的最短环路。我们采用排列编码表示解：

• 每个解是城市索引的排列，如[0,3,1,2]表示0→3→1→2→0的路线
• 评价函数计算环路总距离：

def calculate_total_distance(route): total = 0 for i in range(len(route)-1): total += distance_matrix[route[i]][route[i+1]] total += distance_matrix[route[-1]][route[0]] # 返回起点 return total

距离矩阵预先计算并存储，避免重复运算：

# 生成距离矩阵示例 def create_distance_matrix(coords): n = len(coords) dist_matrix = np.zeros((n, n)) for i in range(n): for j in range(n): dx = coords[i][0] - coords[j][0] dy = coords[i][1] - coords[j][1] dist_matrix[i][j] = np.sqrt(dx**2 + dy**2) return dist_matrix

3. 算法核心实现细节

3.1 初始解生成

随机排列虽然简单，但结合贪心策略可以获得更好的起点：

def generate_initial_solution(method='random'): if method == 'random': return np.random.permutation(num_cities) elif method == 'nearest_neighbor': start = np.random.randint(num_cities) unvisited = set(range(num_cities)) - {start} route = [start] while unvisited: last = route[-1] next_city = min(unvisited, key=lambda x: distance_matrix[last][x]) route.append(next_city) unvisited.remove(next_city) return np.array(route)

3.2 邻域搜索策略

采用2-opt交换产生新解，比简单交换两个城市更有效：

def two_opt_swap(route): i, j = sorted(np.random.choice(len(route), 2, replace=False)) new_route = route.copy() new_route[i:j+1] = route[i:j+1][::-1] # 反转片段 return new_route

3.3 退火调度设计

冷却进度表对算法性能至关重要，我们实现三种降温策略：

def cooling_schedule(T, method='exponential', **params): if method == 'exponential': return params['alpha'] * T elif method == 'linear': return T - params['delta'] elif method == 'logarithmic': return T / (1 + params['beta'] * T)

典型参数设置：

• 初始温度：使初始接受概率≈0.8
• 终止温度：1e-6
• 降温系数α：0.85-0.99
• 每个温度迭代次数：100-1000

4. 完整算法实现与可视化

将各组件整合成完整算法：

def simulated_annealing_tsp(city_coords, max_iter=10000): # 初始化 dist_matrix = create_distance_matrix(city_coords) current_route = generate_initial_solution('nearest_neighbor') current_dist = calculate_total_distance(current_route) best_route, best_dist = current_route.copy(), current_dist # 温度参数 T = initial_temperature(current_dist) cooling = lambda t: 0.95 * t # 记录过程 history = [] for i in range(max_iter): # 生成新解 new_route = two_opt_swap(current_route) new_dist = calculate_total_distance(new_route) # Metropolis准则 if new_dist < current_dist or \ np.random.rand() < np.exp((current_dist - new_dist)/T): current_route, current_dist = new_route, new_dist # 更新最优解 if current_dist < best_dist: best_route, best_dist = current_route.copy(), current_dist # 降温 T = cooling(T) history.append(best_dist) # 提前终止 if T < 1e-6: break return best_route, best_dist, history

结果可视化展示优化过程：

def plot_results(route, history, city_coords): plt.figure(figsize=(12,5)) # 优化曲线 plt.subplot(121) plt.plot(history) plt.title('Optimization Progress') plt.xlabel('Iteration') plt.ylabel('Total Distance') # 路线图 plt.subplot(122) x = [city_coords[i][0] for i in route] + [city_coords[route[0]][0]] y = [city_coords[i][1] for i in route] + [city_coords[route[0]][1]] plt.plot(x, y, 'o-') plt.title('Best Route Found') plt.xlabel('X Coordinate') plt.ylabel('Y Coordinate') plt.tight_layout() plt.show()

5. 性能优化与实用技巧

5.1 加速计算的关键

• 向量化距离计算：利用NumPy广播机制

def vectorized_distance(route): shifted = np.roll(route, -1) pairs = np.column_stack((route, shifted)) return np.sum(distance_matrix[pairs[:,0], pairs[:,1]])

• 增量计算：只计算受交换影响的部分

def delta_distance(route, i, j): n = len(route) a, b = route[i], route[(i+1)%n] c, d = route[j], route[(j+1)%n] # 移除旧边，添加新边 return (distance_matrix[a][d] + distance_matrix[c][b]) - \ (distance_matrix[a][b] + distance_matrix[c][d])

5.2 参数调优指南

通过实验分析参数影响：

5.3 混合策略提升效果

结合局部搜索增强收敛性：

def local_search(route): improved = True while improved: improved = False for i in range(len(route)): for j in range(i+1, len(route)): delta = delta_distance(route, i, j) if delta < 0: route[i+1:j+1] = route[i+1:j+1][::-1] improved = True return route

在退火过程中定期调用局部搜索，可以显著提升解的质量。