旅行商问题 (TSP)#

数学描述为，给定 $n$ 个城市和它们之间的距离 $d_{ij},\ i,j \in \{1, 2, \dots, n\}$，要求找到一条经过每个城市一次且仅一次的最短路径。

解空间 $\mathcal{X}$ 是遍历每个城市恰好一次的回路，可以表示为 $\{1, 2, \dots, n\}$ 的全部循环排列。记状态为 $x = (x_1, x_2, \dots, x_n)$，其中 $x_i$ 表示第 $i$ 个城市的编号。目标函数为路径长度，即 $f(x) = d_{x_1 x_2} + d_{x_2 x_3} + \dots + d_{x_{n-1} x_n} + d_{x_n x_1}$。

新解的产生可以通过互换操作、插入操作、逆序操作来实现。

• 互换操作：随机选择两个城市 $i$ 和 $j$，交换它们在路径中的位置。
• 插入操作：随机选择一个城市 $i$ 和一个位置 $j$，将城市 $i$ 插入到路径中的位置 $j$。
• 逆序操作：随机选择两个城市 $i$ 和 $j$，将它们之间的路径逆序。

e.g. 设当前路径为 $x = (1, 2, 3, 4, 5, 6, 7, 8, 9)$，互换操作产生的新解为 $x_{new} = (1, \underline{8}, 3, 4, 5, 6, 7, \underline{2}, 9)$，插入操作产生的新解为 $x_{new} = (1, \underline{8}, 2, 3, 4, 5, 6, 7, 9)$，逆序操作产生的新解为 $x_{new} = (1, \underline{8, 7, 6, 5, 4, 3, 2}, 9)$。

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

def generate_distance_matrix(n_cities: int, seed: int = 42) -> np.ndarray:
    """
    生成随机城市坐标和距离矩阵
    
    参数:
    - n_cities: 城市数量
    - seed: 随机种子
    
    返回:
    - distances: 距离矩阵
    - coords: 城市坐标
    """
    random.seed(seed)
    np.random.seed(seed)
    
    # 生成随机坐标
    coords = np.random.rand(n_cities, 2) * 100
    
    # 计算距离矩阵
    distances = np.zeros((n_cities, n_cities))
    for i in range(n_cities):
        for j in range(n_cities):
            if i != j:
                # 欧几里得距离
                distances[i, j] = np.sqrt(np.sum((coords[i] - coords[j])**2))
    
    return distances, coords

def calculate_path_length(path: List[int], distances: np.ndarray) -> float:
    """
    计算路径长度
    
    参数:
    - path: 路径
    - distances: 距离矩阵
    
    返回:
    - path_length: 路径长度
    """
    path_length = 0
    for i in range(len(path) - 1):
        path_length += distances[path[i], path[i+1]]
    # 加上从最后一个城市回到起点的距离
    path_length += distances[path[-1], path[0]]
    return path_length

def random_path(n_cities: int) -> List[int]:
    """
    生成随机路径
    
    参数:
    - n_cities: 城市数量
    
    返回:
    - path: 随机路径
    """
    path = list(range(n_cities))
    random.shuffle(path)
    return path

def generate_neighbor(path: List[int]) -> List[int]:
    """
    生成邻居解（使用三种操作之一：互换、插入或逆序）
    
    参数:
    - path: 当前路径
    
    返回:
    - new_path: 新路径
    """
    new_path = path.copy()
    n = len(path)
    
    # 随机选择操作类型
    operation = random.choice(["swap", "insert", "reverse"])
    
    if operation == "swap":
        # 互换操作：随机选择两个城市并交换它们的位置
        i, j = random.sample(range(n), 2)
        new_path[i], new_path[j] = new_path[j], new_path[i]
    
    elif operation == "insert":
        # 插入操作：随机选择一个城市，并将其移动到一个新位置
        i, j = random.sample(range(n), 2)
        city = new_path.pop(i)
        new_path.insert(j, city)
    
    else:  # operation == "reverse"
        # 逆序操作：随机选择两个位置，并将它们之间的路径逆序
        i, j = sorted(random.sample(range(n), 2))
        new_path[i:j+1] = reversed(new_path[i:j+1])
    
    return new_path

def simulated_annealing_tsp(distances: np.ndarray, 
                           initial_temp: float = 1000.0, 
                           cooling_rate: float = 0.995, 
                           min_temp: float = 1e-8,
                           iterations_per_temp: int = 100) -> Tuple[List[int], float, List[float], List[float]]:
    """
    使用模拟退火算法解决旅行商问题
    
    参数:
    - distances: 城市间距离矩阵
    - initial_temp: 初始温度
    - cooling_rate: 冷却率
    - min_temp: 最小温度
    - iterations_per_temp: 每个温度下的迭代次数
    
    返回:
    - best_path: 找到的最佳路径
    - best_length: 最佳路径长度
    - length_history: 路径长度历史
    - temp_history: 温度历史
    """
    n_cities = distances.shape[0]
    
    # 初始化随机路径
    current_path = random_path(n_cities)
    current_length = calculate_path_length(current_path, distances)
    
    best_path = current_path.copy()
    best_length = current_length
    
    temp = initial_temp
    
    # 存储历史数据用于可视化
    length_history = [current_length]
    temp_history = [temp]
    
    # 主循环
    while temp > min_temp:
        for _ in range(iterations_per_temp):
            # 生成新解
            new_path = generate_neighbor(current_path)
            new_length = calculate_path_length(new_path, distances)
            
            # 计算能量差
            delta = new_length - current_length
            
            # Metropolis准则：接受更好的解或以一定概率接受较差的解
            if delta < 0 or random.random() < math.exp(-delta / temp):
                current_path = new_path
                current_length = new_length
                
                # 更新全局最优解
                if current_length < best_length:
                    best_path = current_path.copy()
                    best_length = current_length
            
            length_history.append(current_length)
            temp_history.append(temp)
        
        # 降温
        temp *= cooling_rate
    
    return best_path, best_length, length_history, temp_history

def plot_tsp_solution(coords, path, title="TSP Solution"):
    """绘制TSP解决方案"""
    plt.figure(figsize=(10, 8))
    
    # 绘制城市
    plt.scatter(coords[:, 0], coords[:, 1], c='red', s=50)
    
    # 绘制路径
    for i in range(len(path)):
        j = (i + 1) % len(path)
        plt.plot([coords[path[i], 0], coords[path[j], 0]], 
                 [coords[path[i], 1], coords[path[j], 1]], 'b-')
    
    # 添加城市编号标签
    for i, coord in enumerate(coords):
        plt.text(coord[0] + 1, coord[1] + 1, str(i), fontsize=10)
    
    plt.title(title)
    plt.grid(True)
    plt.tight_layout()
    plt.show()

def plot_convergence(length_history, temp_history):
    """绘制收敛过程"""
    plt.figure(figsize=(14, 6))
    
    # 绘制路径长度变化
    plt.subplot(1, 2, 1)
    plt.plot(length_history)
    plt.xlabel('迭代次数')
    plt.ylabel('路径长度')
    plt.title('路径长度随迭代变化')
    plt.grid(True)
    
    # 绘制温度变化
    plt.subplot(1, 2, 2)
    plt.plot(temp_history)
    plt.xlabel('迭代次数')
    plt.ylabel('温度')
    plt.title('温度随迭代变化')
    plt.yscale('log')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

# 主程序
if __name__ == "__main__":
    # 生成问题实例（20个城市）
    n_cities = 20
    distances, coords = generate_distance_matrix(n_cities)
    
    # 运行模拟退火算法
    initial_temp = 100.0
    cooling_rate = 0.995
    iterations_per_temp = 50
    
    print("开始求解TSP问题...")
    best_path, best_length, length_history, temp_history = simulated_annealing_tsp(
        distances, initial_temp, cooling_rate, min_temp=0.01, iterations_per_temp=iterations_per_temp
    )
    
    print(f"找到的最佳路径: {best_path}")
    print(f"最佳路径长度: {best_length:.2f}")
    
    # 可视化结果
    plot_tsp_solution(coords, best_path, f"TSP Solution - Length: {best_length:.2f}")
    plot_convergence(length_history, temp_history)