用Python模拟退火算法搞定旅行商问题（TSP）：从物理退火到代码实现的保姆级拆解

山月刀岚月刀

292人浏览 · 2026-05-30 11:41:25

山月刀岚月刀 · 2026-05-30 11:41:25 发布

从物理退火到Python实战：用模拟退火算法优雅解决旅行商问题

想象一下，你是一位物流公司的算法工程师，老板扔给你一份包含50个配送点的城市地图，要求你在半小时内规划出最短的配送路线。传统方法可能需要几天时间计算，而模拟退火算法能在喝杯咖啡的功夫给出90分以上的解决方案——这就是我们今天要探索的算法魔法。

1. 物理世界的启示：退火过程与算法思想的奇妙对应

1983年，Kirkpatrick等科学家观察到一个有趣现象：金属在加热后缓慢冷却时，原子会从高能态逐渐排列成稳定的晶体结构。这个过程被称为退火(annealing)，它启发了计算机科学家们解决复杂优化问题的思路。

温度与粒子运动的关系 ：

高温状态：金属原子剧烈运动，能量高且排列无序
缓慢降温：原子逐渐找到更低能态的位置
低温状态：原子形成稳定晶体结构，系统能量最低

把这个物理过程映射到算法中：

原子状态 → 问题的可能解
系统能量 → 目标函数值(如TSP中的路径长度)
温度参数 → 控制算法探索行为的参数

# 物理概念与算法参数的对应关系
physics_to_algorithm = {
    "原子排列": "解向量",
    "系统能量": "目标函数",
    "温度": "接受劣解的概率控制",
    "冷却速率": "参数衰减系数"
}

2. 算法核心：Metropolis准则的智慧

1953年，Metropolis提出了一个改变优化算法历史的准则：当系统从状态i转移到状态j时：

如果E(j)≤E(i)，接受新状态j
如果E(j)>E(i)，以概率P=exp[-(E(j)-E(i))/T]接受j

这个简单的规则赋予了算法"战略性冒险"的能力——允许暂时接受劣质解以避免陷入局部最优。

关键参数解析 ：

参数	物理意义	算法作用	典型取值
T₀	初始温度	控制初始探索范围	1e5-1e7
α	冷却系数	控制降温速度	0.85-0.99
Tf	终止温度	停止条件	0.1-1.0
L	马尔可夫链长	每温度迭代次数	100-1000

注意：过快的冷却(α太小)可能导致"淬火"现象，算法过早收敛到局部最优

3. TSP问题建模：从业务需求到数学表达

旅行商问题(TSP)可以表述为：给定n个城市及其相互距离，找到访问每个城市一次并返回起点的最短路径。

DFJ模型数学表达 ：

min ΣΣ dᵢⱼxᵢⱼ
约束条件：
1. 每个城市恰好进入一次：Σxᵢⱼ = 1, ∀i
2. 每个城市恰好离开一次：Σxᵢⱼ = 1, ∀j
3. 消除子回路：Σxᵢⱼ ≤ |S|-1, ∀S⊂V
4. 决策变量：xᵢⱼ ∈ {0,1}

在实际编程中，我们通常采用更简单的整数编码表示路径，如[0,3,1,2]表示0→3→1→2→0的环路。

4. Python实现：从理论到完整代码

让我们用Python实现一个完整的模拟退火算法解决TSP问题。代码分为几个关键部分：

4.1 数据准备与初始化

import numpy as np
import matplotlib.pyplot as plt

# 城市坐标数据
city_coordinates = {
    0: (60, 200), 1: (180, 200), 2: (80, 180), 
    3: (140, 180), 4: (20, 160), 5: (100, 160)
}

num_city = len(city_coordinates)

# 计算距离矩阵
def create_distance_matrix(coords):
    n = len(coords)
    dist_mat = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            if i != j:
                xi, yi = coords[i]
                xj, yj = coords[j]
                dist_mat[i][j] = np.sqrt((xi-xj)**2 + (yi-yj)**2)
    return dist_mat

distance_matrix = create_distance_matrix(city_coordinates)

4.2 核心算法实现

def simulated_annealing(cities, max_iter=1000, t0=10000, alpha=0.95):
    current_solution = np.random.permutation(len(cities))
    current_cost = calculate_cost(current_solution)
    best_solution = current_solution.copy()
    best_cost = current_cost
    
    t = t0
    costs = []
    
    for _ in range(max_iter):
        # 生成邻域解
        new_solution = generate_neighbor(current_solution)
        new_cost = calculate_cost(new_solution)
        
        # 计算成本差
        delta = new_cost - current_cost
        
        # Metropolis准则
        if delta < 0 or np.random.rand() < np.exp(-delta/t):
            current_solution = new_solution
            current_cost = new_cost
            
            if current_cost < best_cost:
                best_solution = current_solution.copy()
                best_cost = current_cost
                
        # 降温
        t *= alpha
        costs.append(best_cost)
        
    return best_solution, best_cost, costs

# 评价函数
def calculate_cost(solution):
    total = 0
    for i in range(len(solution)-1):
        total += distance_matrix[solution[i]][solution[i+1]]
    total += distance_matrix[solution[-1]][solution[0]]
    return total

# 邻域生成(2-opt交换)
def generate_neighbor(solution):
    new_solution = solution.copy()
    i, j = np.random.choice(len(solution), 2, replace=False)
    new_solution[i], new_solution[j] = new_solution[j], new_solution[i]
    return new_solution

4.3 结果可视化与分析

# 运行算法
best_route, best_cost, cost_history = simulated_annealing(city_coordinates)

# 绘制优化过程
plt.plot(cost_history)
plt.title('Optimization Process')
plt.xlabel('Iteration')
plt.ylabel('Total Distance')
plt.show()

# 绘制最优路径
def plot_route(route, coords):
    x = [coords[i][0] for i in route] + [coords[route[0]][0]]
    y = [coords[i][1] for i in route] + [coords[route[0]][1]]
    
    plt.figure(figsize=(10,6))
    plt.plot(x, y, 'o-')
    for i, (xi, yi) in enumerate(coords.values()):
        plt.text(xi, yi, str(i))
    plt.title(f'Best Route Found (Total Distance: {best_cost:.2f})')
    plt.show()

plot_route(best_route, city_coordinates)

5. 高级技巧与实战建议

在实际应用中，我们还需要考虑以下优化点：

参数调优策略 ：

初始温度：应该足够高，使得初始接受概率≈1
- 可通过计算初始随机解的成本方差来估计
冷却计划：尝试不同的降温策略
- 指数冷却：T = αT
- 对数冷却：T = T₀/log(1+k)
终止条件：除了温度，还可结合最大迭代次数或无改进次数

算法改进方向 ：

自适应温度调节：根据接受率动态调整降温速率
混合邻域搜索：结合2-opt、3-opt等不同邻域结构
并行化实现：同时探索多个解空间区域

实用建议：对于50个城市以上的问题，建议将最大迭代次数设置为至少10000次，并考虑使用更高效的邻域生成方法

6. 与其他优化算法的对比

模拟退火并非唯一解决TSP的方法，下面是比较几种常见算法：

算法	优点	缺点	适用场景
模拟退火	避免局部最优，实现简单	参数敏感，收敛慢	中小规模问题(≤100城市)
遗传算法	并行搜索，鲁棒性强	编码复杂，早熟收敛	中等规模问题
蚁群算法	正反馈机制，分布式计算	参数多，计算量大	路径特征明显的问题
精确算法	保证最优解	计算复杂度高	小规模问题(≤20城市)

在最近的一个物流配送项目中，我们对100个配送点测试发现：

模拟退火在5分钟内找到了比遗传算法更好的解
但最终采用了混合策略：用模拟退火生成初始解，再用2-opt局部优化

7. 常见问题与解决方案

Q1：算法陷入局部最优怎么办？

增加初始温度
尝试更慢的冷却速率
引入重启机制(温度回升)

Q2：如何评估解的质量？

与已知最优解比较(对标准测试集)
多次运行看结果稳定性
计算gap = (found - best_known)/best_known

Q3：处理大规模TSP的建议？

分治策略：先聚类再分别求解
使用更高效的邻域结构(如3-opt)
考虑基于深度学习的现代方法

# 示例：带重启机制的改进版本
def advanced_SA(cities, max_restarts=3):
    best_global = None
    best_cost = float('inf')
    
    for _ in range(max_restarts):
        solution, cost, _ = simulated_annealing(cities)
        if cost < best_cost:
            best_global = solution
            best_cost = cost
            
    return best_global, best_cost

8. 实际应用案例扩展

让我们看一个更真实的案例——某电商公司在双十一期间的配送路线优化：

业务约束 ：

50个配送点，包含优先级(生鲜优先)
车辆容量限制
时间窗口约束

算法调整 ：

修改评价函数，加入惩罚项：

def calculate_cost(solution):
    distance = 0
    penalty = 0
    
    # 计算基本距离
    for i in range(len(solution)-1):
        distance += dist_matrix[solution[i]][solution[i+1]]
    
    # 添加时间窗和容量惩罚
    if violates_time_window(solution):
        penalty += 10000
    if violates_capacity(solution):
        penalty += 5000
        
    return distance + penalty