模拟退火算法实战：用Python解决旅行商问题（TSP）的完整代码解析

每次面对一个看似无解的复杂优化问题，比如规划一条覆盖几十个城市的最短路线，或者为一个庞大的数据中心寻找最优的服务器布局，那种感觉就像面对一团乱麻。传统的穷举法在问题规模稍大时就变得遥不可及，而一些贪婪算法又常常一头扎进局部最优的陷阱里出不来。这时候，一种灵感来源于物理世界的算法——模拟退火，就成了我们工具箱里一件优雅而强大的武器。它不保证找到绝对的最优解，但在有限的时间和资源内，它往往能给出一个令人惊喜的、接近最优的答案。这篇文章，就是为你，一位希望将理论算法转化为实际生产力的Python开发者，准备的一份从原理到实战的深度指南。我们将一起动手，用代码实现模拟退火算法，并把它应用在经典的旅行商问题上，看着算法如何一步步“退火”，从一团混乱中结晶出优美的路径。

1. 理解模拟退火：物理灵感与算法内核

在深入代码之前，我们有必要先抛开那些复杂的数学公式，从最直观的物理图像来理解模拟退火。想象一下，你是一位铁匠，手中有一块烧得通红的铁。此时，铁内部的原子处于高度活跃、无序的状态。你的目标是通过“退火”工艺，让这些原子重新排列，形成坚固、稳定的晶体结构。这个过程的关键在于缓慢降温。如果你把烧红的铁直接丢进冷水里淬火，原子会被“冻”在某个高能量的无序状态，材料会变得硬而脆。但如果你让它在一个高温炉中慢慢冷却，原子就有足够的时间，通过热运动的扰动，逐渐“滑入”能量更低的稳定位置。

模拟退火算法正是对这一物理过程的绝妙模拟：

• 解的状态 对应材料的微观状态（原子排列）。
• 目标函数值（如路径总长度） 对应系统的能量。我们的目标是找到能量最低（目标函数值最小）的状态。
• 温度 是算法中一个核心的控制参数。高温时，算法接受“坏解”（使能量升高的新状态）的概率大，探索能力强；随着温度降低，接受坏解的概率变小，算法逐渐稳定，倾向于在好的解附近进行局部精细搜索。

算法的灵魂在于其接受新解的 Metropolis准则。这个准则允许算法以一定的概率接受一个比当前解更差的解。正是这个机制，赋予了算法跳出局部最优陷阱的能力。接受差解的概率 p 由公式决定：

p = exp(-ΔE / T)

其中 ΔE 是新解与旧解的目标函数值之差（对于最小化问题，ΔE > 0 表示解变差），T 是当前温度。从这个公式我们可以直观地看到：

• 当 T 很高时，即使 ΔE 很大，p 也可能接近1，算法几乎“肆无忌惮”地探索解空间。
• 当 T 很低时，只有 ΔE 很小的差解才有被接受的可能，算法行为更接近传统的局部搜索。
• 当 ΔE 为负（解变好）时，p 恒大于1，意味着好解总是被接受。

注意：理解“温度”在算法中的抽象意义至关重要。它不是一个物理量，而是一个控制算法“探索”与“利用”平衡的衰减参数。初始温度设置、降温速度（退火计划）是影响算法性能的关键超参数。

为了更清晰地对比模拟退火与其他常见优化策略的核心区别，我们可以看下面这个简单的表格：

特性	模拟退火 (SA)	梯度下降法	随机搜索
核心思想	模拟物理退火过程，以概率接受差解	沿目标函数梯度反方向迭代	在解空间中完全随机采样
跳出局部最优能力	强，得益于Metropolis准则	弱，容易陷入最近的局部极小点	理论上强，但效率极低
收敛性	理论上能以概率1收敛到全局最优（无限时间）	收敛到局部最优	不保证收敛
对目标函数要求	极低，只需能计算函数值	需要可微	极低
主要超参数	初始温度、降温速率、马尔可夫链长度	学习率	采样次数
适用场景	组合优化、非凸函数、离散问题	连续可微凸/非凸问题	任何问题，但效率是瓶颈

2. 问题定义：旅行商问题（TSP）的建模

旅行商问题是一个经典的NP-hard组合优化问题，描述非常简单：一个商人需要访问N个城市，每个城市只访问一次，最后回到起点，要求找到总距离最短的访问路线。尽管描述简单，但其解空间随着城市数量N呈阶乘级((N-1)!/2)增长，使得精确求解在N较大时变得不可能。

用模拟退火解决TSP，我们首先需要将问题“翻译”成算法能处理的形式：

1. 解的表达：一个解就是城市的一个排列（permutation）。例如，对于4个城市[A, B, C, D]，一个可能的解是 [A, C, B, D]，表示访问顺序。
2. 目标函数：计算该排列下，按顺序遍历所有城市并返回起点的总距离。我们需要一个计算距离的函数。这里我们使用欧几里得距离作为示例。
3. 邻域操作（产生新解）：这是模拟退火在TSP上的关键设计。如何从当前解产生一个“邻居”解？常用的操作有：
   • 交换（Swap）：随机选择两个位置，交换这两个位置上的城市。
   • 逆转（Reverse/2-opt）：随机选择一段子路径，将其顺序完全反转。这是TSP中非常高效的一种邻域操作。
   • 插入（Insert）：随机选择一个城市，将其插入到另一个随机位置。

在接下来的实现中，我们将主要使用交换和逆转操作来生成新解，你会发现它们对解的质量有不同影响。

我们先来搭建项目的基础结构，并实现城市坐标生成、距离计算和可视化函数。这是所有后续工作的基石。

import numpy as np
import matplotlib.pyplot as plt
import random
import math
from typing import List, Tuple
import time

class TSPProblem:
    """定义TSP问题实例"""
    def __init__(self, num_cities: int = 20, seed: int = 42):
        """
        初始化TSP问题，随机生成城市坐标。
        
        参数:
            num_cities: 城市数量
            seed: 随机种子，确保结果可复现
        """
        np.random.seed(seed)
        self.num_cities = num_cities
        # 在[0, 100]的二维平面内随机生成城市坐标
        self.coords = np.random.rand(num_cities, 2) * 100
        # 预计算距离矩阵，加速距离查询
        self._compute_distance_matrix()

    def _compute_distance_matrix(self):
        """计算并存储所有城市两两之间的欧氏距离矩阵"""
        self.dist_matrix = np.zeros((self.num_cities, self.num_cities))
        for i in range(self.num_cities):
            for j in range(i+1, self.num_cities):
                dist = np.linalg.norm(self.coords[i] - self.coords[j])
                self.dist_matrix[i][j] = dist
                self.dist_matrix[j][i] = dist # 距离矩阵是对称的

    def get_distance(self, city_a: int, city_b: int) -> float:
        """快速获取两个城市间的距离"""
        return self.dist_matrix[city_a][city_b]

    def total_distance(self, tour: List[int]) -> float:
        """
        计算给定路径序列的总距离。
        
        参数:
            tour: 城市索引的列表，例如 [0, 3, 1, 2]
        返回:
            路径的总长度
        """
        total = 0.0
        num = len(tour)
        for i in range(num):
            # 从当前城市到下一个城市的距离，最后一个城市连接到起点
            total += self.get_distance(tour[i], tour[(i + 1) % num])
        return total

    def plot_tour(self, tour: List[int], title: str = "TSP Path"):
        """可视化路径"""
        plt.figure(figsize=(10, 6))
        # 绘制城市点
        plt.scatter(self.coords[:, 0], self.coords[:, 1], c='red', s=100, zorder=5)
        for i, (x, y) in enumerate(self.coords):
            plt.text(x, y, str(i), fontsize=12, ha='center', va='center', zorder=6)

        # 绘制路径线
        tour_coords = self.coords[tour + [tour[0]]] # 使路径闭合
        plt.plot(tour_coords[:, 0], tour_coords[:, 1], 'b-', linewidth=1, alpha=0.8)
        plt.scatter(tour_coords[0, 0], tour_coords[0, 1], c='green', s=150, marker='*', zorder=7, label='Start/End')

        plt.xlabel("X Coordinate")
        plt.ylabel("Y Coordinate")
        plt.title(f"{title} | Total Distance: {self.total_distance(tour):.2f}")
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.axis('equal')
        plt.tight_layout()
        plt.show()

# 示例：创建一个包含15个城市的问题并随机生成一条路径看看
if __name__ == "__main__":
    problem = TSPProblem(num_cities=15)
    random_tour = list(range(15))
    random.shuffle(random_tour)
    print(f"随机路径的总距离: {problem.total_distance(random_tour):.2f}")
    problem.plot_tour(random_tour, "Random Initial Tour")

运行上面的代码，你会看到一张散点图，红色的点代表城市，蓝色的线是随机生成的一条访问路径。这个距离通常很大，路径交叉严重，我们的目标就是让模拟退火算法把它“熨平”。

3. 模拟退火算法核心实现

有了问题定义和基础工具，我们现在可以构建模拟退火算法的核心引擎了。这个部分我们将把之前讨论的物理概念——温度、退火计划、Metropolis准则——一一用代码实现。

一个健壮的模拟退火实现需要考虑以下几个模块：

1. 初始解生成：一个简单的随机排列即可。
2. 邻域移动生成器：实现交换、逆转等操作。
3. 退火计划：决定温度如何随着迭代下降。常见的有线性衰减、指数衰减等。
4. 内循环（马尔可夫链长度）：在每个温度下，进行多少次邻域搜索尝试。
5. 停止准则：何时结束算法？例如温度低于阈值、连续若干代解无改进等。

下面是我们完整的模拟退火求解器实现。我添加了大量注释，并设计了灵活的接口，方便你后续调整策略和参数。

class SimulatedAnnealingTSP:
    """模拟退火算法求解TSP"""
    def __init__(self, problem: TSPProblem, initial_temperature: float = 1000.0,
                 cooling_rate: float = 0.995, min_temperature: float = 1e-3,
                 iterations_per_temp: int = 1000, use_2opt: bool = True):
        """
        初始化SA求解器。
        
        参数:
            problem: TSPProblem实例
            initial_temperature: 初始温度，控制早期探索性
            cooling_rate: 降温系数，每次迭代温度乘以这个系数（指数退火）
            min_temperature: 最低温度，低于此温度算法终止
            iterations_per_temp: 每个温度下的迭代次数（马尔可夫链长度）
            use_2opt: 是否在邻域操作中使用2-opt（逆转）作为主要扰动方式
        """
        self.problem = problem
        self.T_init = initial_temperature
        self.cooling_rate = cooling_rate
        self.T_min = min_temperature
        self.iter_per_temp = iterations_per_temp
        self.use_2opt = use_2opt

        # 记录历史数据用于分析
        self.best_tour_history = []
        self.best_distance_history = []
        self.current_distance_history = []
        self.temperature_history = []

    def _generate_initial_solution(self) -> List[int]:
        """生成初始解：随机排列"""
        tour = list(range(self.problem.num_cities))
        random.shuffle(tour)
        return tour

    def _generate_neighbor(self, tour: List[int]) -> List[int]:
        """
        通过邻域操作产生一个新解（邻居）。
        这里我们随机选择两种操作之一：交换两个城市，或进行2-opt逆转。
        """
        new_tour = tour.copy()
        n = len(new_tour)

        if self.use_2opt and random.random() > 0.3: # 70%的概率使用2-opt
            # 2-opt操作：随机选择两个索引i, j (i < j)，逆转i到j之间的子路径
            i = random.randint(0, n-2)
            j = random.randint(i+1, n-1)
            new_tour[i:j+1] = reversed(new_tour[i:j+1])
        else:
            # 交换操作：随机选择两个不同的位置进行交换
            i, j = random.sample(range(n), 2)
            new_tour[i], new_tour[j] = new_tour[j], new_tour[i]

        return new_tour

    def _calculate_delta_distance(self, old_tour: List[int], new_tour: List[int]) -> float:
        """
        高效计算新旧路径的距离变化量ΔE。
        注意：由于我们只做了局部改动（交换或逆转），可以只计算受影响部分距离的变化，
        而不必重新计算整条路径，这能极大提升性能。这里为清晰起见，我们使用完整计算。
        在实际高性能实现中，强烈建议实现增量计算。
        """
        old_dist = self.problem.total_distance(old_tour)
        new_dist = self.problem.total_distance(new_tour)
        return new_dist - old_dist

    def solve(self, verbose: bool = True) -> Tuple[List[int], float]:
        """
        执行模拟退火算法主循环。
        
        返回:
            best_tour: 找到的最佳路径
            best_distance: 最佳路径对应的距离
        """
        # 初始化
        current_tour = self._generate_initial_solution()
        current_distance = self.problem.total_distance(current_tour)
        best_tour = current_tour.copy()
        best_distance = current_distance

        T = self.T_init
        iteration = 0

        if verbose:
            print(f"开始模拟退火优化...")
            print(f"初始随机路径距离: {best_distance:.2f}")
            print("-" * 50)

        # 主退火循环
        while T > self.T_min:
            accepted_moves = 0
            for _ in range(self.iter_per_temp):
                # 1. 产生邻居解
                new_tour = self._generate_neighbor(current_tour)
                delta_dist = self._calculate_delta_distance(current_tour, new_tour)

                # 2. 根据Metropolis准则决定是否接受新解
                if delta_dist < 0:
                    # 新解更好，总是接受
                    accept = True
                else:
                    # 新解更差，以概率 exp(-ΔE / T) 接受
                    probability = math.exp(-delta_dist / T)
                    accept = random.random() < probability

                if accept:
                    current_tour = new_tour
                    current_distance += delta_dist # 更新当前距离
                    accepted_moves += 1

                    # 3. 更新历史最优解
                    if current_distance < best_distance:
                        best_tour = current_tour.copy()
                        best_distance = current_distance
                        if verbose and iteration % 500 == 0:
                            print(f"Iter {iteration:6d}, T={T:.4f}, Best Dist={best_distance:.2f}")

                # 记录数据
                self.current_distance_history.append(current_distance)
                self.best_distance_history.append(best_distance)
                self.temperature_history.append(T)
                iteration += 1

            # 4. 降温
            T *= self.cooling_rate

            # 可选：根据接受率动态调整马尔可夫链长度（高级技巧）
            # acceptance_rate = accepted_moves / self.iter_per_temp
            # if acceptance_rate < 0.1:
            #     # 接受率太低，可以提前升温或跳出循环
            #     pass

        if verbose:
            print("-" * 50)
            print(f"优化完成！共迭代 {iteration} 次。")
            print(f"找到的最短路径距离: {best_distance:.2f}")
            print(f"相比初始解提升: {(self.problem.total_distance(self._generate_initial_solution()) - best_distance) / best_distance * 100:.1f}% (近似)")

        self.best_tour_history = best_tour
        return best_tour, best_distance

    def plot_optimization_process(self):
        """绘制优化过程曲线，展示距离和温度随迭代次数的变化"""
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        iterations = range(len(self.best_distance_history))

        ax1.plot(iterations, self.best_distance_history, 'g-', linewidth=1.5, label='Best Distance', alpha=0.7)
        ax1.plot(iterations, self.current_distance_history, 'r-', linewidth=0.5, label='Current Distance', alpha=0.4)
        ax1.set_xlabel('Iteration')
        ax1.set_ylabel('Distance')
        ax1.set_title('Simulated Annealing Optimization Progress')
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        ax2.plot(iterations, self.temperature_history, 'b-', linewidth=1)
        ax2.set_xlabel('Iteration')
        ax2.set_ylabel('Temperature (log scale)')
        ax2.set_title('Temperature Schedule')
        ax2.set_yscale('log')
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        plt.show()

现在，让我们创建一个问题实例并用我们的求解器来跑一下，看看效果。

# 创建并解决问题
problem = TSPProblem(num_cities=25, seed=123)

# 初始化求解器，参数可以调整
solver = SimulatedAnnealingTSP(problem,
                               initial_temperature=1000,
                               cooling_rate=0.995,
                               min_temperature=1e-3,
                               iterations_per_temp=1000,
                               use_2opt=True)

start_time = time.time()
best_tour, best_dist = solver.solve(verbose=True)
end_time = time.time()

print(f"\n优化耗时: {end_time - start_time:.2f} 秒")

# 可视化最终结果
problem.plot_tour(best_tour, title="Optimized Tour by Simulated Annealing")

# 查看优化过程
solver.plot_optimization_process()

运行这段代码，你会在控制台看到算法迭代的过程，最终输出优化后的路径长度和耗时。两张图会分别展示优化后的路径（应该比最初的随机路径整洁、交叉少得多）以及算法搜索过程中“当前解距离”和“历史最优解距离”的变化曲线。你会注意到，当前解距离（红线）上下跳动得很厉害，尤其是在高温阶段，这正是算法在探索解空间的表现。而历史最优解距离（绿线）则呈现一个阶梯式下降的趋势。

4. 参数调优与高级技巧：让算法更高效

模拟退火算法好用，但其性能极大地依赖于参数设置。一套糟糕的参数可能让算法在解空间里盲目游荡很久也找不到好解。这一节，我们就像调试精密仪器一样，来聊聊如何调整这些“旋钮”。

4.1 核心参数解析与调优指南

• 初始温度 (initial_temperature)：

  • 作用：决定算法初期的“探索野心”。温度太高，初期会接受大量差解，搜索随机，收敛慢；温度太低，则过早陷入局部搜索，可能跳不出局部最优。
  • 调优方法：一个经验法则是，让初始温度下，接受差解的概率在一个较高的水平（例如0.8以上）。可以通过一个小实验来设定：随机产生大量邻域移动，计算目标函数变化的平均值 ΔE_avg，然后根据 T0 = -ΔE_avg / ln(p) 来估算，其中 p 是你期望的初始接受概率。

• 降温速率 (cooling_rate)：

  • 作用：控制温度下降的速度，是“退火计划”的核心。指数退火 T_{k+1} = α * T_k 是最常用的方式。
  • 调优方法：α 通常取值在 [0.9, 0.999] 之间。值越接近1，降温越慢，在每个温度下搜索得越充分，但耗时也越长。对于复杂问题，需要更慢的退火（更大的α）。可以尝试 0.99, 0.995, 0.999 等值。

• 每个温度的迭代次数 (iterations_per_temp)：

  • 作用：也称为马尔可夫链长度。它决定了在每个温度下，算法尝试搜索的邻域解的数量。
  • 调优方法：通常与问题规模相关。对于TSP，可以设置为城市数量的若干倍（如 100*n 到 1000*n）。太短则搜索不充分，太长则增加不必要的计算。一个高级技巧是使其与温度挂钩，高温时短一些（快速探索），低温时长一些（精细搜索）。

• 终止温度 (min_temperature)：

  • 作用：算法停止的条件之一。当温度低于此阈值时，算法几乎不再接受差解，可以终止。
  • 调优方法：通常设为一个很小的正数，如 1e-3, 1e-5。也可以结合其他终止条件，如连续若干代最优解无改进。

为了直观感受参数的影响，我们可以设计一个小实验，对比不同降温速率的效果：

def compare_cooling_rates(problem, rates=[0.99, 0.995, 0.999]):
    """比较不同降温速率对最终结果的影响"""
    results = {}
    for rate in rates:
        print(f"\n测试降温速率 α = {rate}")
        solver = SimulatedAnnealingTSP(problem,
                                       initial_temperature=1000,
                                       cooling_rate=rate,
                                       min_temperature=1e-3,
                                       iterations_per_temp=500) # 固定迭代次数
        start = time.time()
        best_tour, best_dist = solver.solve(verbose=False)
        elapsed = time.time() - start
        results[rate] = {'distance': best_dist, 'time': elapsed}
        print(f"  最短距离: {best_dist:.2f}, 耗时: {elapsed:.2f}s")

    # 简单结果对比
    print("\n=== 结果对比 ===")
    for rate, data in results.items():
        print(f"α={rate}: 距离 {data['distance']:.2f}, 时间 {data['time']:.2f}s")

运行这个比较函数，你会发现 α=0.999 的慢速退火通常能找到更好的解，但花费的时间也显著更长。这体现了优化中永恒的 “效果-效率”权衡。

4.2 邻域操作的进阶选择

我们之前实现了交换和2-opt两种邻域操作。实际上，针对TSP，学术界和工业界有更多高效的邻域结构：

• 3-opt：比2-opt更复杂的操作，通过断开路径的三条边并重新连接，能产生更大的扰动，跳出更深局部最优的能力更强，但每次评估的计算量也更大。
• Lin-Kernighan (LK) 启发式：这是一种非常强大的局部搜索启发式，可以看作是2-opt和3-opt的智能组合与迭代。它被广泛用于TSP的高性能求解器中。将模拟退火与LK结合（例如，以一定概率在SA的每次迭代中调用LK进行深度局部搜索），可以极大提升解的质量。

下面是一个概念性的代码片段，展示如何将更复杂的邻域操作集成到我们的框架中：

class AdvancedSATSPSolver(SimulatedAnnealingTSP):
    """集成更多高级邻域操作的SA求解器"""
    def _generate_neighbor_advanced(self, tour: List[int]) -> List[int]:
        """以不同概率尝试多种邻域操作"""
        rand_val = random.random()
        new_tour = tour.copy()
        n = len(new_tour)

        if rand_val < 0.6: # 60%概率使用2-opt
            i = random.randint(0, n-2)
            j = random.randint(i+1, n-1)
            new_tour[i:j+1] = reversed(new_tour[i:j+1])
        elif rand_val < 0.9: # 30%概率使用交换
            i, j = random.sample(range(n), 2)
            new_tour[i], new_tour[j] = new_tour[j], new_tour[i]
        else: # 10%概率尝试插入操作
            i = random.randint(0, n-1)
            city = new_tour.pop(i)
            j = random.randint(0, n-1)
            new_tour.insert(j, city)
        return new_tour

    # 可以在这里重写 _generate_neighbor 方法，调用上面的高级版本

4.3 自适应退火与重启策略

为了让算法更智能，我们可以引入一些自适应机制：

• 自适应退火计划：根据当前搜索状态动态调整降温速率。例如，如果当前温度下接受率一直很高，说明温度可能还太高，可以加快降温；如果接受率骤降，可以适当减缓降温甚至短暂“回温”。
• 重启策略：当算法在低温下陷入停滞（最优解长时间不更新）时，可以保存当前最优解，然后从一个新的随机解（或对当前最优解施加一个较大扰动）开始，并重置到一个中等温度，重新进行退火。这能有效增加找到全局最优的概率。

# 一个简单的重启策略伪代码示例
def solve_with_restart(self, max_restarts=5):
    global_best_tour = None
    global_best_dist = float('inf')

    for restart in range(max_restarts):
        # 每次重启，可以稍微扰动上次的全局最优解作为初始解，而非完全随机
        if global_best_tour:
            current_tour = self._perturb_solution(global_best_tour, strength=0.1)
        else:
            current_tour = self._generate_initial_solution()

        # 用稍低的初始温度开始新一轮退火
        T = self.T_init / (restart + 1)
        # ... 执行标准SA循环 ...

        # 更新全局最优
        if current_best_dist < global_best_dist:
            global_best_dist = current_best_dist
            global_best_tour = current_best_tour.copy()

        print(f"重启 {restart+1}/{max_restarts} 完成，当前全局最优: {global_best_dist:.2f}")

    return global_best_tour, global_best_dist

参数调优没有银弹，最好的方法是对你的特定问题（城市分布、规模）进行多次实验，记录不同参数组合下的结果（最终距离、运行时间、收敛曲线），然后根据你的需求（是追求极致解质量，还是要求快速得到一个可接受的解）来选择最合适的配置。这个过程本身，也充满了探索的乐趣。