traveling-salesman-problem

旅行商问题（TSP）

你是旅行商问题和组合优化领域的专家。你的目标是帮助找到最短路线，使其恰好访问每座城市一次并返回起点城市，同时最小化总行程距离或时间。

初步评估

在求解 TSP 实例之前，请先了解：

1. 问题特征

   • 需要访问多少座城市／地点？
   • 是对称 TSP（距离 i→j = 距离 j→i）还是非对称 TSP？
   • 使用什么度量？（欧几里得距离、行程时间、成本）
   • 有什么特殊约束吗？（时间窗 → 参见 vrp-time-windows）

2. 问题规模

   • 小型（< 20 座城市）：可行精确方法
   • 中型（20-100 座城市）：高级精确方法或优质启发式方法
   • 大型（100-1000 座城市）：需要元启发式方法
   • 超大型（> 1000 座城市）：需要专用算法

3. 求解要求

   • 需要经过证明的最优解吗？
   • 可接受的最优性差距是多少？（例如，与最优解相差不超过 5%）
   • 求解时间有何限制？
   • 需要单条路线还是多条路线？

4. 数据格式

   • 坐标（纬度／经度、x/y）？
   • 是否提供距离矩阵？
   • 是否需要计算距离？
   • 使用真实道路网络还是欧几里得距离？

---

数学建模

经典 TSP 建模（MTZ）

决策变量：

• x_{ij} ∈ {0,1}：若弧 (i,j) 属于路线，则为 1；否则为 0
• u_i ∈ ℝ：城市 i 在路线中的位置（用于消除子回路）

参数：

• c_{ij}：从城市 i 到城市 j 的成本／距离
• n：城市数量

目标函数：

Minimize: Σ_{i=1}^n Σ_{j=1}^n c_{ij} * x_{ij}

约束：

1. Each city has exactly one outgoing arc:
   Σ_{j=1}^n x_{ij} = 1,  ∀i

2. Each city has exactly one incoming arc:
   Σ_{i=1}^n x_{ij} = 1,  ∀j

3. Subtour elimination (Miller-Tucker-Zemlin):
   u_i - u_j + n*x_{ij} ≤ n-1,  ∀i,j ∈ {2,...,n}, i≠j

4. Variable bounds:
   x_{ij} ∈ {0,1}
   2 ≤ u_i ≤ n,  ∀i ∈ {2,...,n}

备选方案：DFJ（Dantzig-Fulkerson-Johnson）建模

子回路消除约束：

Σ_{i∈S} Σ_{j∈S} x_{ij} ≤ |S| - 1,  ∀S ⊂ V, 2 ≤ |S| ≤ n-1

此建模包含指数数量的约束，但能产生更紧的 LP 松弛。

---

精确算法

1. 分支定界与动态规划

import numpy as np
from itertools import combinations

def tsp_dynamic_programming(dist_matrix):
    """
    Held-Karp algorithm for TSP using dynamic programming

    Time complexity: O(n^2 * 2^n)
    Space complexity: O(n * 2^n)

    Optimal for small instances (n ≤ 20)

    Args:
        dist_matrix: n x n distance matrix

    Returns:
        tuple: (optimal_cost, optimal_tour)
    """
    n = len(dist_matrix)

    # C[S][i] = minimum cost path visiting all nodes in S ending at i
    C = {}

    # Initialize: paths of length 1
    for i in range(1, n):
        C[(1 << i, i)] = (dist_matrix[0][i], 0)

    # Iterate over subset sizes
    for subset_size in range(2, n):
        # Generate all subsets of size subset_size
        for subset in combinations(range(1, n), subset_size):
            # Convert to bitmask
            bits = 0
            for bit in subset:
                bits |= 1 << bit

            # Find minimum cost to reach each node in subset
            for i in subset:
                # Previous subset without node i
                prev_bits = bits & ~(1 << i)

                min_cost = float('inf')
                min_prev = None

                # Try all possible previous nodes
                for j in subset:
                    if j == i:
                        continue

                    if (prev_bits, j) in C:
                        cost = C[(prev_bits, j)][0] + dist_matrix[j][i]
                        if cost < min_cost:
                            min_cost = cost
                            min_prev = j

                if min_prev is not None:
                    C[(bits, i)] = (min_cost, min_prev)

    # Find optimal tour
    bits = (1 << n) - 2  # All nodes except 0
    min_cost = float('inf')
    last_node = None

    for i in range(1, n):
        if (bits, i) in C:
            cost = C[(bits, i)][0] + dist_matrix[i][0]
            if cost < min_cost:
                min_cost = cost
                last_node = i

    # Reconstruct tour
    tour = [0]
    bits = (1 << n) - 2

    while last_node is not None:
        tour.append(last_node)
        prev_bits = bits & ~(1 << last_node)

        if (bits, last_node) in C:
            last_node = C[(bits, last_node)][1]
            bits = prev_bits
        else:
            break

    tour.append(0)

    return min_cost, tour


# Example usage
if __name__ == "__main__":
    # Small 5-city example
    dist_matrix = np.array([
        [0, 10, 15, 20, 25],
        [10, 0, 35, 25, 30],
        [15, 35, 0, 30, 20],
        [20, 25, 30, 0, 15],
        [25, 30, 20, 15, 0]
    ])

    optimal_cost, optimal_tour = tsp_dynamic_programming(dist_matrix)
    print(f"Optimal Cost: {optimal_cost}")
    print(f"Optimal Tour: {optimal_tour}")

2. 使用 PuLP 的 MIP 建模

from pulp import *
import numpy as np

def tsp_mip_mtz(dist_matrix, city_names=None):
    """
    TSP using Miller-Tucker-Zemlin formulation

    Suitable for medium-sized instances (up to 50-100 cities)

    Args:
        dist_matrix: n x n distance matrix
        city_names: optional list of city names

    Returns:
        dict with optimal_cost, tour, and solve_time
    """
    n = len(dist_matrix)

    if city_names is None:
        city_names = [f"City_{i}" for i in range(n)]

    # Create problem
    prob = LpProblem("TSP_MTZ", LpMinimize)

    # Decision variables
    # x[i,j] = 1 if arc from city i to city j is in tour
    x = {}
    for i in range(n):
        for j in range(n):
            if i != j:
                x[i,j] = LpVariable(f"x_{i}_{j}", cat='Binary')

    # u[i] = position of city i in tour (for subtour elimination)
    u = {}
    for i in range(1, n):
        u[i] = LpVariable(f"u_{i}", lowBound=1, upBound=n-1, cat='Continuous')

    # Objective: Minimize total distance
    prob += lpSum([dist_matrix[i][j] * x[i,j]
                   for i in range(n) for j in range(n) if i != j]), \
            "Total_Distance"

    # Constraints

    # 1. Each city has exactly one outgoing arc
    for i in range(n):
        prob += lpSum([x[i,j] for j in range(n) if j != i]) == 1, \
                f"Out_{i}"

    # 2. Each city has exactly one incoming arc
    for j in range(n):
        prob += lpSum([x[i,j] for i in range(n) if i != j]) == 1, \
                f"In_{j}"

    # 3. Miller-Tucker-Zemlin subtour elimination
    for i in range(1, n):
        for j in range(1, n):
            if i != j:
                prob += u[i] - u[j] + n*x[i,j] <= n-1, \
                        f"MTZ_{i}_{j}"

    # Solve
    import time
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] == 'Optimal':
        tour = [0]
        current_city = 0

        for _ in range(n-1):
            for j in range(n):
                if j != current_city and x[current_city,j].varValue > 0.5:
                    tour.append(j)
                    current_city = j
                    break

        tour.append(0)  # Return to start

        return {
            'status': 'Optimal',
            'optimal_cost': value(prob.objective),
            'tour': tour,
            'tour_names': [city_names[i] for i in tour],
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'optimal_cost': None,
            'tour': None,
            'solve_time': solve_time
        }


# Example usage
dist_matrix = np.array([
    [0, 29, 20, 21, 16, 31],
    [29, 0, 15, 29, 28, 40],
    [20, 15, 0, 15, 14, 25],
    [21, 29, 15, 0, 4, 12],
    [16, 28, 14, 4, 0, 16],
    [31, 40, 25, 12, 16, 0]
])

city_names = ['Depot', 'Customer_A', 'Customer_B',
              'Customer_C', 'Customer_D', 'Customer_E']

result = tsp_mip_mtz(dist_matrix, city_names)
print(f"\nStatus: {result['status']}")
print(f"Optimal Cost: {result['optimal_cost']:.2f}")
print(f"Optimal Tour: {result['tour_names']}")
print(f"Solve Time: {result['solve_time']:.2f} seconds")

---

构造型启发式方法

1. 最近邻法

import numpy as np

def nearest_neighbor_tsp(dist_matrix, start_city=0):
    """
    Nearest Neighbor heuristic for TSP

    Time complexity: O(n^2)
    Quality: Typically within 25% of optimal

    Args:
        dist_matrix: n x n distance matrix
        start_city: starting city (default 0)

    Returns:
        tuple: (total_cost, tour)
    """
    n = len(dist_matrix)
    unvisited = set(range(n))
    tour = [start_city]
    unvisited.remove(start_city)

    current_city = start_city
    total_cost = 0

    while unvisited:
        # Find nearest unvisited city
        nearest_city = min(unvisited,
                          key=lambda city: dist_matrix[current_city][city])

        total_cost += dist_matrix[current_city][nearest_city]
        tour.append(nearest_city)
        unvisited.remove(nearest_city)
        current_city = nearest_city

    # Return to start
    total_cost += dist_matrix[current_city][start_city]
    tour.append(start_city)

    return total_cost, tour


def multi_start_nearest_neighbor(dist_matrix, num_starts=None):
    """
    Run Nearest Neighbor from multiple starting cities
    and return best solution

    Args:
        dist_matrix: n x n distance matrix
        num_starts: number of different starts (default: all cities)

    Returns:
        tuple: (best_cost, best_tour)
    """
    n = len(dist_matrix)

    if num_starts is None:
        num_starts = n

    best_cost = float('inf')
    best_tour = None

    for start in range(min(num_starts, n)):
        cost, tour = nearest_neighbor_tsp(dist_matrix, start)

        if cost < best_cost:
            best_cost = cost
            best_tour = tour

    return best_cost, best_tour

2. 最小成本插入法

def cheapest_insertion_tsp(dist_matrix):
    """
    Cheapest Insertion heuristic for TSP

    Time complexity: O(n^3)
    Quality: Generally better than Nearest Neighbor

    Args:
        dist_matrix: n x n distance matrix

    Returns:
        tuple: (total_cost, tour)
    """
    n = len(dist_matrix)

    # Start with smallest edge
    min_dist = float('inf')
    best_pair = (0, 1)

    for i in range(n):
        for j in range(i+1, n):
            if dist_matrix[i][j] < min_dist:
                min_dist = dist_matrix[i][j]
                best_pair = (i, j)

    # Initial tour
    tour = list(best_pair)
    unvisited = set(range(n)) - set(tour)

    # Insert remaining cities
    while unvisited:
        best_insertion = None
        best_cost_increase = float('inf')
        best_position = None

        # Try inserting each unvisited city
        for city in unvisited:
            # Try inserting between each pair of adjacent cities
            for i in range(len(tour)):
                j = (i + 1) % len(tour)

                # Cost increase of inserting city between tour[i] and tour[j]
                cost_increase = (dist_matrix[tour[i]][city] +
                               dist_matrix[city][tour[j]] -
                               dist_matrix[tour[i]][tour[j]])

                if cost_increase < best_cost_increase:
                    best_cost_increase = cost_increase
                    best_insertion = city
                    best_position = j

        # Insert best city
        tour.insert(best_position, best_insertion)
        unvisited.remove(best_insertion)

    # Calculate total cost
    total_cost = sum(dist_matrix[tour[i]][tour[(i+1)%len(tour)]]
                    for i in range(len(tour)))

    tour.append(tour[0])  # Close tour

    return total_cost, tour


def farthest_insertion_tsp(dist_matrix):
    """
    Farthest Insertion heuristic (variant)

    Similar to cheapest insertion but selects farthest city
    from tour at each step
    """
    n = len(dist_matrix)

    # Start with two farthest cities
    max_dist = 0
    best_pair = (0, 1)

    for i in range(n):
        for j in range(i+1, n):
            if dist_matrix[i][j] > max_dist:
                max_dist = dist_matrix[i][j]
                best_pair = (i, j)

    tour = list(best_pair)
    unvisited = set(range(n)) - set(tour)

    while unvisited:
        # Find farthest city from current tour
        farthest_city = None
        max_min_dist = 0

        for city in unvisited:
            min_dist_to_tour = min(dist_matrix[city][tour_city]
                                  for tour_city in tour)
            if min_dist_to_tour > max_min_dist:
                max_min_dist = min_dist_to_tour
                farthest_city = city

        # Find best insertion position
        best_position = 0
        best_cost_increase = float('inf')

        for i in range(len(tour)):
            j = (i + 1) % len(tour)
            cost_increase = (dist_matrix[tour[i]][farthest_city] +
                           dist_matrix[farthest_city][tour[j]] -
                           dist_matrix[tour[i]][tour[j]])

            if cost_increase < best_cost_increase:
                best_cost_increase = cost_increase
                best_position = j

        tour.insert(best_position, farthest_city)
        unvisited.remove(farthest_city)

    total_cost = sum(dist_matrix[tour[i]][tour[(i+1)%len(tour)]]
                    for i in range(len(tour)))
    tour.append(tour[0])

    return total_cost, tour

---

改进型启发式方法

1. 2-Opt 局部搜索

def two_opt(tour, dist_matrix):
    """
    2-opt improvement heuristic

    Iteratively removes two edges and reconnects the tour
    in a different way if it improves the solution

    Args:
        tour: initial tour (list of city indices)
        dist_matrix: n x n distance matrix

    Returns:
        tuple: (improved_cost, improved_tour)
    """
    n = len(tour) - 1  # Exclude duplicate last city
    improved = True
    best_tour = tour[:-1]  # Remove last city (duplicate of first)

    while improved:
        improved = False

        for i in range(1, n - 1):
            for j in range(i + 1, n):
                # Current edges: (tour[i-1], tour[i]) and (tour[j], tour[j+1])
                # New edges: (tour[i-1], tour[j]) and (tour[i], tour[j+1])

                current_cost = (dist_matrix[best_tour[i-1]][best_tour[i]] +
                               dist_matrix[best_tour[j]][best_tour[(j+1)%n]])

                new_cost = (dist_matrix[best_tour[i-1]][best_tour[j]] +
                           dist_matrix[best_tour[i]][best_tour[(j+1)%n]])

                if new_cost < current_cost:
                    # Reverse the segment between i and j
                    best_tour[i:j+1] = reversed(best_tour[i:j+1])
                    improved = True
                    break

            if improved:
                break

    # Calculate total cost
    total_cost = sum(dist_matrix[best_tour[i]][best_tour[(i+1)%n]]
                    for i in range(n))
    total_cost += dist_matrix[best_tour[n-1]][best_tour[0]]

    best_tour.append(best_tour[0])  # Close tour

    return total_cost, best_tour


def two_opt_optimized(tour, dist_matrix, max_iterations=1000):
    """
    Optimized 2-opt with early stopping
    """
    n = len(tour) - 1
    best_tour = tour[:-1]

    def calculate_tour_cost(t):
        return sum(dist_matrix[t[i]][t[(i+1)%len(t)]] for i in range(len(t)))

    best_cost = calculate_tour_cost(best_tour)
    iterations = 0
    no_improvement_count = 0

    while iterations < max_iterations and no_improvement_count < 50:
        iterations += 1
        improved = False

        for i in range(1, n - 1):
            for j in range(i + 1, n):
                # Calculate delta (change in cost)
                delta = (dist_matrix[best_tour[i-1]][best_tour[j]] +
                        dist_matrix[best_tour[i]][best_tour[(j+1)%n]] -
                        dist_matrix[best_tour[i-1]][best_tour[i]] -
                        dist_matrix[best_tour[j]][best_tour[(j+1)%n]])

                if delta < -1e-10:  # Improvement found
                    best_tour[i:j+1] = reversed(best_tour[i:j+1])
                    best_cost += delta
                    improved = True
                    no_improvement_count = 0
                    break

            if improved:
                break

        if not improved:
            no_improvement_count += 1

    best_tour.append(best_tour[0])
    return best_cost, best_tour

2. 3-Opt 局部搜索

def three_opt(tour, dist_matrix, max_iterations=100):
    """
    3-opt improvement heuristic

    More powerful than 2-opt but slower
    Considers removing 3 edges and reconnecting

    Args:
        tour: initial tour
        dist_matrix: n x n distance matrix
        max_iterations: maximum iterations

    Returns:
        tuple: (improved_cost, improved_tour)
    """
    n = len(tour) - 1
    best_tour = tour[:-1]

    def calculate_cost(t):
        return sum(dist_matrix[t[i]][t[(i+1)%len(t)]] for i in range(len(t)))

    best_cost = calculate_cost(best_tour)

    for iteration in range(max_iterations):
        improved = False

        for i in range(n - 2):
            for j in range(i + 2, n - 1):
                for k in range(j + 2, n):
                    # Try all possible 3-opt reconnections
                    # There are 8 ways to reconnect 3 segments

                    # Current tour segments: A-B-C where
                    # A = tour[0:i+1]
                    # B = tour[i+1:j+1]
                    # C = tour[j+1:k+1]
                    # D = tour[k+1:]

                    # Try different reconnections
                    segments = [
                        best_tour[:i+1],
                        best_tour[i+1:j+1],
                        best_tour[j+1:k+1],
                        best_tour[k+1:]
                    ]

                    # Case 1: A + reverse(B) + C + D
                    new_tour = segments[0] + segments[1][::-1] + segments[2] + segments[3]
                    new_cost = calculate_cost(new_tour)

                    if new_cost < best_cost:
                        best_tour = new_tour
                        best_cost = new_cost
                        improved = True
                        break

                if improved:
                    break

            if improved:
                break

        if not improved:
            break

    best_tour.append(best_tour[0])
    return best_cost, best_tour

3. Or-Opt

def or_opt(tour, dist_matrix):
    """
    Or-opt improvement heuristic

    Relocates sequences of 1, 2, or 3 consecutive cities

    Args:
        tour: initial tour
        dist_matrix: n x n distance matrix

    Returns:
        tuple: (improved_cost, improved_tour)
    """
    n = len(tour) - 1
    best_tour = tour[:-1]
    improved = True

    def calculate_cost(t):
        return sum(dist_matrix[t[i]][t[(i+1)%len(t)]] for i in range(len(t)))

    best_cost = calculate_cost(best_tour)

    while improved:
        improved = False

        # Try relocating sequences of length 1, 2, and 3
        for seq_length in [1, 2, 3]:
            if seq_length >= n:
                continue

            for i in range(n):
                if i + seq_length > n:
                    continue

                # Extract sequence
                sequence = best_tour[i:i+seq_length]

                # Try inserting at all other positions
                for j in range(n):
                    if j >= i and j < i + seq_length:
                        continue

                    # Create new tour with sequence relocated
                    new_tour = best_tour[:i] + best_tour[i+seq_length:]
                    new_tour = new_tour[:j] + sequence + new_tour[j:]

                    new_cost = calculate_cost(new_tour)

                    if new_cost < best_cost - 1e-10:
                        best_tour = new_tour
                        best_cost = new_cost
                        improved = True
                        break

                if improved:
                    break

            if improved:
                break

    best_tour.append(best_tour[0])
    return best_cost, best_tour

---

元启发式方法

1. 模拟退火

import random
import math

def simulated_annealing_tsp(dist_matrix, initial_temp=1000,
                           cooling_rate=0.995, max_iterations=10000):
    """
    Simulated Annealing for TSP

    Probabilistically accepts worse solutions to escape local optima

    Args:
        dist_matrix: n x n distance matrix
        initial_temp: starting temperature
        cooling_rate: temperature reduction factor
        max_iterations: maximum iterations

    Returns:
        tuple: (best_cost, best_tour)
    """
    n = len(dist_matrix)

    # Generate initial solution
    current_tour = list(range(n))
    random.shuffle(current_tour)
    current_tour.append(current_tour[0])

    def calculate_cost(tour):
        return sum(dist_matrix[tour[i]][tour[i+1]]
                  for i in range(len(tour)-1))

    current_cost = calculate_cost(current_tour)
    best_tour = current_tour.copy()
    best_cost = current_cost

    temperature = initial_temp

    for iteration in range(max_iterations):
        # Generate neighbor solution (2-opt move)
        i = random.randint(1, n - 2)
        j = random.randint(i + 1, n - 1)

        new_tour = current_tour.copy()
        new_tour[i:j+1] = reversed(new_tour[i:j+1])

        new_cost = calculate_cost(new_tour)
        delta = new_cost - current_cost

        # Accept or reject new solution
        if delta < 0 or random.random() < math.exp(-delta / temperature):
            current_tour = new_tour
            current_cost = new_cost

            # Update best solution
            if current_cost < best_cost:
                best_tour = current_tour.copy()
                best_cost = current_cost

        # Cool down
        temperature *= cooling_rate

        if temperature < 0.01:
            break

    return best_cost, best_tour

2. 遗传算法

def genetic_algorithm_tsp(dist_matrix, population_size=100,
                         generations=500, mutation_rate=0.01):
    """
    Genetic Algorithm for TSP

    Uses order crossover (OX) and swap mutation

    Args:
        dist_matrix: n x n distance matrix
        population_size: number of individuals
        generations: number of generations
        mutation_rate: probability of mutation

    Returns:
        tuple: (best_cost, best_tour)
    """
    n = len(dist_matrix)

    def calculate_fitness(tour):
        cost = sum(dist_matrix[tour[i]][tour[(i+1)%n]] for i in range(n))
        return 1.0 / (1.0 + cost)  # Higher fitness = better tour

    def create_individual():
        tour = list(range(n))
        random.shuffle(tour)
        return tour

    def order_crossover(parent1, parent2):
        """Order crossover (OX)"""
        size = len(parent1)
        start, end = sorted(random.sample(range(size), 2))

        child = [-1] * size
        child[start:end] = parent1[start:end]

        # Fill remaining positions from parent2
        pos = end
        for city in parent2[end:] + parent2[:end]:
            if city not in child:
                if pos >= size:
                    pos = 0
                child[pos] = city
                pos += 1

        return child

    def mutate(tour):
        """Swap mutation"""
        if random.random() < mutation_rate:
            i, j = random.sample(range(len(tour)), 2)
            tour[i], tour[j] = tour[j], tour[i]
        return tour

    def tournament_selection(population, fitnesses, k=3):
        """Tournament selection"""
        selected = random.sample(list(zip(population, fitnesses)), k)
        return max(selected, key=lambda x: x[1])[0]

    # Initialize population
    population = [create_individual() for _ in range(population_size)]

    best_tour = None
    best_cost = float('inf')

    for generation in range(generations):
        # Calculate fitness
        fitnesses = [calculate_fitness(ind) for ind in population]

        # Track best
        for tour, fitness in zip(population, fitnesses):
            cost = sum(dist_matrix[tour[i]][tour[(i+1)%n]] for i in range(n))
            if cost < best_cost:
                best_cost = cost
                best_tour = tour.copy()

        # Create next generation
        new_population = []

        # Elitism: keep best individuals
        elite_count = int(0.1 * population_size)
        elite_indices = sorted(range(len(fitnesses)),
                              key=lambda i: fitnesses[i],
                              reverse=True)[:elite_count]
        new_population = [population[i].copy() for i in elite_indices]

        # Crossover and mutation
        while len(new_population) < population_size:
            parent1 = tournament_selection(population, fitnesses)
            parent2 = tournament_selection(population, fitnesses)

            child = order_crossover(parent1, parent2)
            child = mutate(child)

            new_population.append(child)

        population = new_population

    best_tour.append(best_tour[0])
    return best_cost, best_tour

3. 蚁群优化

def ant_colony_optimization_tsp(dist_matrix, n_ants=20, n_iterations=100,
                                alpha=1.0, beta=2.0, evaporation=0.5, Q=100):
    """
    Ant Colony Optimization for TSP

    Args:
        dist_matrix: n x n distance matrix
        n_ants: number of ants
        n_iterations: number of iterations
        alpha: pheromone importance
        beta: distance importance
        evaporation: pheromone evaporation rate
        Q: pheromone deposit factor

    Returns:
        tuple: (best_cost, best_tour)
    """
    n = len(dist_matrix)

    # Initialize pheromone matrix
    pheromone = np.ones((n, n)) / n

    best_tour = None
    best_cost = float('inf')

    for iteration in range(n_iterations):
        all_tours = []
        all_costs = []

        # Each ant constructs a tour
        for ant in range(n_ants):
            tour = [0]  # Start from city 0
            unvisited = set(range(1, n))

            current_city = 0

            while unvisited:
                # Calculate probabilities for next city
                probabilities = []

                for city in unvisited:
                    tau = pheromone[current_city][city] ** alpha
                    eta = (1.0 / dist_matrix[current_city][city]) ** beta
                    probabilities.append(tau * eta)

                # Normalize probabilities
                total = sum(probabilities)
                probabilities = [p / total for p in probabilities]

                # Select next city
                next_city = np.random.choice(list(unvisited), p=probabilities)

                tour.append(next_city)
                unvisited.remove(next_city)
                current_city = next_city

            tour.append(0)  # Return to start

            # Calculate tour cost
            cost = sum(dist_matrix[tour[i]][tour[i+1]]
                      for i in range(len(tour)-1))

            all_tours.append(tour)
            all_costs.append(cost)

            # Update best solution
            if cost < best_cost:
                best_cost = cost
                best_tour = tour.copy()

        # Evaporate pheromone
        pheromone *= (1 - evaporation)

        # Deposit pheromone
        for tour, cost in zip(all_tours, all_costs):
            deposit = Q / cost

            for i in range(len(tour) - 1):
                city1, city2 = tour[i], tour[i+1]
                pheromone[city1][city2] += deposit
                pheromone[city2][city1] += deposit

    return best_cost, best_tour

---

使用 OR-Tools（Google 优化工具）

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

def solve_tsp_ortools(dist_matrix, time_limit_seconds=30):
    """
    Solve TSP using Google OR-Tools

    Very efficient for practical instances

    Args:
        dist_matrix: n x n distance matrix
        time_limit_seconds: time limit for solver

    Returns:
        dict with solution details
    """
    n = len(dist_matrix)

    # Create routing index manager
    manager = pywrapcp.RoutingIndexManager(n, 1, 0)

    # Create routing model
    routing = pywrapcp.RoutingModel(manager)

    # Create distance callback
    def distance_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return int(dist_matrix[from_node][to_node])

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Set search parameters
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.seconds = time_limit_seconds

    # Solve
    solution = routing.SolveWithParameters(search_parameters)

    if solution:
        tour = []
        index = routing.Start(0)

        while not routing.IsEnd(index):
            tour.append(manager.IndexToNode(index))
            index = solution.Value(routing.NextVar(index))

        tour.append(manager.IndexToNode(index))

        return {
            'status': 'Optimal' if solution.status == 1 else 'Feasible',
            'optimal_cost': solution.ObjectiveValue(),
            'tour': tour,
            'computation_time': routing.solver().WallTime() / 1000.0
        }
    else:
        return {
            'status': 'No solution found',
            'optimal_cost': None,
            'tour': None
        }


# Example: Solve TSP with coordinates
def solve_tsp_from_coordinates(coordinates, time_limit=30):
    """
    Solve TSP given (x, y) coordinates

    Args:
        coordinates: list of (x, y) tuples
        time_limit: time limit in seconds

    Returns:
        solution dictionary
    """
    import math

    n = len(coordinates)

    # Calculate Euclidean distance matrix
    dist_matrix = np.zeros((n, n))

    for i in range(n):
        for j in range(n):
            if i != j:
                dx = coordinates[i][0] - coordinates[j][0]
                dy = coordinates[i][1] - coordinates[j][1]
                dist_matrix[i][j] = math.sqrt(dx*dx + dy*dy)

    return solve_tsp_ortools(dist_matrix, time_limit)


# Example usage
coordinates = [
    (0, 0),    # Depot
    (2, 4),    # City 1
    (5, 1),    # City 2
    (8, 3),    # City 3
    (6, 7),    # City 4
    (3, 6),    # City 5
    (7, 9),    # City 6
    (1, 8)     # City 7
]

result = solve_tsp_from_coordinates(coordinates)
print(f"Status: {result['status']}")
print(f"Optimal Cost: {result['optimal_cost']:.2f}")
print(f"Optimal Tour: {result['tour']}")
print(f"Computation Time: {result['computation_time']:.2f} seconds")

---

完整求解框架

class TSPSolver:
    """
    Complete TSP solver with multiple algorithms
    """

    def __init__(self, dist_matrix, city_names=None):
        self.dist_matrix = np.array(dist_matrix)
        self.n = len(dist_matrix)

        if city_names is None:
            self.city_names = [f"City_{i}" for i in range(self.n)]
        else:
            self.city_names = city_names

    def solve(self, method='ortools', **kwargs):
        """
        Solve TSP using specified method

        Methods:
        - 'ortools': Google OR-Tools (recommended)
        - 'dynamic_programming': Exact (small instances only)
        - 'nearest_neighbor': Fast heuristic
        - 'cheapest_insertion': Good heuristic
        - '2opt': Local search improvement
        - 'simulated_annealing': Metaheuristic
        - 'genetic': Genetic algorithm
        - 'aco': Ant colony optimization
        """

        if method == 'ortools':
            result = solve_tsp_ortools(self.dist_matrix,
                                      kwargs.get('time_limit', 30))

        elif method == 'dynamic_programming':
            if self.n > 20:
                raise ValueError("Dynamic programming only for n ≤ 20")
            cost, tour = tsp_dynamic_programming(self.dist_matrix)
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Optimal'}

        elif method == 'nearest_neighbor':
            cost, tour = multi_start_nearest_neighbor(
                self.dist_matrix,
                kwargs.get('num_starts', self.n))
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Heuristic'}

        elif method == 'cheapest_insertion':
            cost, tour = cheapest_insertion_tsp(self.dist_matrix)
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Heuristic'}

        elif method == '2opt':
            # Need initial solution
            _, initial_tour = nearest_neighbor_tsp(self.dist_matrix)
            cost, tour = two_opt_optimized(initial_tour, self.dist_matrix)
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Local Search'}

        elif method == 'simulated_annealing':
            cost, tour = simulated_annealing_tsp(
                self.dist_matrix,
                kwargs.get('initial_temp', 1000),
                kwargs.get('cooling_rate', 0.995),
                kwargs.get('max_iterations', 10000))
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Metaheuristic'}

        elif method == 'genetic':
            cost, tour = genetic_algorithm_tsp(
                self.dist_matrix,
                kwargs.get('population_size', 100),
                kwargs.get('generations', 500),
                kwargs.get('mutation_rate', 0.01))
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Metaheuristic'}

        elif method == 'aco':
            cost, tour = ant_colony_optimization_tsp(
                self.dist_matrix,
                kwargs.get('n_ants', 20),
                kwargs.get('n_iterations', 100))
            result = {'optimal_cost': cost, 'tour': tour, 'status': 'Metaheuristic'}

        else:
            raise ValueError(f"Unknown method: {method}")

        # Add city names to result
        if result['tour']:
            result['tour_names'] = [self.city_names[i] for i in result['tour']]

        return result

    def compare_methods(self, methods=['nearest_neighbor', 'cheapest_insertion',
                                      '2opt', 'ortools']):
        """
        Compare multiple solution methods
        """
        import time

        results = []

        for method in methods:
            print(f"Solving with {method}...")
            start_time = time.time()

            try:
                result = self.solve(method)
                solve_time = time.time() - start_time

                results.append({
                    'method': method,
                    'cost': result['optimal_cost'],
                    'tour': result['tour'],
                    'status': result['status'],
                    'time': solve_time
                })
            except Exception as e:
                print(f"  Error: {e}")

        # Create comparison dataframe
        import pandas as pd
        df = pd.DataFrame([{
            'Method': r['method'],
            'Cost': r['cost'],
            'Status': r['status'],
            'Time (s)': f"{r['time']:.3f}"
        } for r in results])

        # Add % from best
        best_cost = df['Cost'].min()
        df['Gap %'] = ((df['Cost'] - best_cost) / best_cost * 100).round(2)

        return df


# Example usage
if __name__ == "__main__":
    # Example distance matrix
    dist_matrix = [
        [0, 29, 20, 21, 16, 31, 100],
        [29, 0, 15, 29, 28, 40, 72],
        [20, 15, 0, 15, 14, 25, 81],
        [21, 29, 15, 0, 4, 12, 92],
        [16, 28, 14, 4, 0, 16, 94],
        [31, 40, 25, 12, 16, 0, 95],
        [100, 72, 81, 92, 94, 95, 0]
    ]

    city_names = ['Depot', 'Customer A', 'Customer B', 'Customer C',
                  'Customer D', 'Customer E', 'Customer F']

    solver = TSPSolver(dist_matrix, city_names)

    # Solve with OR-Tools
    result = solver.solve('ortools', time_limit=10)
    print(f"\nOR-Tools Solution:")
    print(f"Cost: {result['optimal_cost']:.2f}")
    print(f"Tour: {' → '.join(result['tour_names'])}")

    # Compare methods
    print("\n" + "="*60)
    print("Comparing Methods:")
    print("="*60)
    comparison = solver.compare_methods()
    print(comparison.to_string(index=False))

---

工具与库

Python 库

优化：

• OR-Tools (Google)：工业级路径规划求解器
• PuLP：MIP 建模
• Pyomo：高级优化建模
• NetworkX：图算法
• python-tsp：专用 TSP 库

启发式方法：

• scikit-opt：遗传算法、模拟退火
• pymoo：多目标优化

可视化：

• matplotlib：基础绘图
• plotly：交互式可视化
• folium：基于地图的可视化

商业求解器

• Gurobi：最先进的 MIP 求解器
• CPLEX：IBM 优化求解器
• Concorde：最佳精确 TSP 求解器（专用）

在线工具

• TSPLIB：标准基准实例
• Concorde TSP Solver：免费在线求解器
• NEOS Server：远程优化服务

---

常见挑战与解决方案

挑战：问题规模过大，无法使用精确方法

问题：

• 城市数量超过 100 座时，精确方法会变得不切实际
• 需要快速获得优质解

解决方案：

• 使用设有时间限制的 OR-Tools
• 多次运行最近邻法 + 2-opt
• 元启发式方法（SA、GA、ACO）
• 问题分解（先聚类，再规划路线）

挑战：非对称 TSP

问题：

• 距离 i→j ≠ 距离 j→i（例如，单行道）
• 需要不同的建模方式

解决方案：

• 针对非对称情况修改 MTZ 建模
• OR-Tools 可自然处理非对称情况
• 启发式方法易于调整

挑战：现实约束

问题：

• 时间窗口、容量、多辆车
• 纯 TSP 不够用

解决方案：

• 有时间约束时，请参阅 vrp-time-windows
• 有容量约束时，请参阅 vehicle-routing-problem
• 有多个配送中心时，请参阅 multi-depot-vrp

挑战：动态 TSP

问题：

• 执行过程中会有新客户加入
• 需要重新优化

解决方案：

• 使用快速启发式算法（最近邻、插入法）
• 实施重新优化触发机制
• 确保部分路线仍然可行

---

输出格式

TSP 解决方案报告

问题摘要：

• 城市数量：25
• 距离度量：欧几里得距离
• 求解方法：使用 Guided Local Search 的 OR-Tools
• 计算时间：2.3 秒

解的质量：

指标	值
总距离	427.85 km
城市数量	25
解的状态	最优
最优性差距	0%

最优路线：

Depot → City_A → City_F → City_C → ... → City_Z → Depot

路线可视化： [显示路线的地图或图表]

统计数据：

• 连续城市之间的平均距离：17.1 km
• 最长路段：City_A → City_B（45.2 km）
• 最短路段：City_F → City_G（3.8 km）

---

需要询问的问题

如果需要更多背景信息：

1. 需要访问多少个城市或地点？
2. 你有坐标或距离矩阵吗？
3. 距离是否对称（i→j = j→i）？
4. 你需要已证明的最优解，还是可接受一个较好的解？
5. 可接受的计算时间是多少？
6. 是否有其他约束（时间窗口、容量）？
7. 这是一次性求解，还是需要重复求解的问题？
8. 你能使用商业求解器（Gurobi、CPLEX）吗？

---

相关技能

• vehicle-routing-problem：适用于有容量约束的多车辆问题
• vrp-time-windows：适用于有时间约束的 TSP
• pickup-delivery-problem：适用于取货和送货配对问题
• route-optimization：适用于实际路线规划
• network-flow-optimization：适用于底层图论
• metaheuristic-optimization：适用于高级求解方法
• optimization-modeling：适用于通用 MIP 建模