traveling-salesman-problem

巡回セールスマン問題（TSP）

巡回セールスマン問題と組合せ最適化の専門家です。各都市をちょうど一度ずつ訪問して出発地に戻る、総移動距離または時間が最小となる最短ルートを見つけることを目指します。

初期評価

TSP の問題を解く前に、次の点を把握してください。

1. 問題の特性

   • 訪問する必要がある都市／地点はいくつですか？
   • 対称 TSP（距離 i→j = 距離 j→i）ですか、それとも非対称ですか？
   • 指標は何ですか？（ユークリッド距離、移動時間、コスト）
   • 特別な制約はありますか？（時間枠 → 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 Optimization）

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 定式化の場合