← 目次 ｜ ← 前: 21 ｜ 次: 23 →

22. 組合せ最適化のヒューリスティクス

テーマ: 配送ルート最適化（巡回セールスマン問題）を貪欲法＋2-opt で解く

学習点: 総当たりの計算量, 近傍探索, 局所最適とランダム再スタート, itertools.permutations, 計算時間の計測

依存: 標準ライブラリのみ ／ 難易度: 上級

実行方法

uv run 22_tsp_heuristic.py

スクリプト冒頭の PEP 723 メタデータ（# /// script）により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

営業所を出発して 10 か所の配送先をすべて 1 回ずつ回り、営業所に戻る最短ルートを探す――物流業界で日常的に現れる巡回セールスマン問題（TSP）です。このスクリプトは同じ 11 拠点の問題を 4 つの方法で解き、解の品質と計算時間を比較します。

比較するのは、(1) 全経路を調べ上げる総当たりの厳密解、(2) 「いま居る場所から最も近い未訪問地点へ行く」を繰り返す最近傍法（貪欲法）、(3) その結果を 2-opt 近傍探索で磨いた解、(4) 出発点を変えて繰り返すマルチスタート、の 4 つです。最後に総当たりの経路数が n とともに爆発する様子を示し、なぜ実務でヒューリスティクスが不可欠なのかを体感させます。

コードの読みどころ

理論的背景

TSP は NP 困難で、n 都市の巡回路は向きの対称性を考慮して (n−1)!/2 通りあります。n=12 で約 2,000 万通り、n=25 では宇宙年齢を超える計算時間になり、厳密解法は小規模問題にしか使えません。2-opt は経路中の 2 辺を繋ぎ替える局所探索（Croes 1958）で、理論保証はないものの実用上は最適解の数 % 以内に収まることが多いとされます。局所最適に陥る弱点は、初期解を変えるランダム再スタート（マルチスタート）で補います。

実行結果の見方

厳密解は距離 310.001 ですが 6.38 秒かかっています。一方、最近傍法は 0.00004 秒で解を出すものの +3.13% 長い経路です。ここに 2-opt を加えると、ほぼゼロ秒のまま厳密解と同じ 310.001 に到達しました。「最適経路」と「2-opt 解」の並びが逆順なのは、同じ巡回路を逆回りに表示しているだけで実質同一です。

最後の「計算量の壁」の表が本題です。総当たりは n=12 で 4 秒、n=15 で 2.4 時間、n=20 で 386 年、n=25 では約 20 億年と爆発します。拠点が数十を超える現実の配送計画では、厳密解をあきらめてヒューリスティクスや分枝限定法を使う理由がこの数字に集約されています。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = []
# ///
"""22: 組合せ最適化のヒューリスティクス -----------------------------------
テーマ: 配送ルート最適化（巡回セールスマン問題）を貪欲法＋2-opt で解く
学習点: 総当たりの計算量, 近傍探索, 局所最適とランダム再スタート,
        itertools.permutations, 計算時間の計測
根拠: TSP は NP困難で、n都市の総経路数は (n-1)!/2。n=12 で約2000万通り。
      2-opt は経路中の2辺を繋ぎ替える局所探索で、実用上 最適解の
      数%以内に収まることが多い（Croes 1958）。
"""
import itertools
import math
import random
import time

random.seed(3)


def distance_matrix(pts):
    n = len(pts)
    return [[math.dist(pts[i], pts[j]) for j in range(n)] for i in range(n)]


def tour_length(tour, D):
    return sum(D[tour[i]][tour[(i + 1) % len(tour)]] for i in range(len(tour)))


def brute_force(D):
    n = len(D)
    best, best_len = None, float("inf")
    for perm in itertools.permutations(range(1, n)):     # 出発点を固定
        t = (0,) + perm
        L = tour_length(t, D)
        if L < best_len:
            best, best_len = t, L
    return list(best), best_len


def nearest_neighbor(D, start=0):
    n = len(D)
    unvisited = set(range(n)) - {start}
    tour = [start]
    while unvisited:
        last = tour[-1]
        nxt = min(unvisited, key=lambda j: D[last][j])
        tour.append(nxt)
        unvisited.remove(nxt)
    return tour


def two_opt(tour, D, max_pass=100):
    """辺 (i-1,i) と (j,j+1) を反転して繋ぎ替える局所探索。"""
    best = tour[:]
    improved = True
    n = len(best)
    while improved and max_pass > 0:
        improved = False
        max_pass -= 1
        for i in range(1, n - 1):
            for j in range(i + 1, n):
                a, b = best[i - 1], best[i]
                c, d = best[j], best[(j + 1) % n]
                if a == c or b == d:
                    continue
                delta = (D[a][c] + D[b][d]) - (D[a][b] + D[c][d])
                if delta < -1e-12:
                    best[i:j + 1] = reversed(best[i:j + 1])
                    improved = True
    return best


def main() -> None:
    n = 11
    pts = [(random.uniform(0, 100), random.uniform(0, 100)) for _ in range(n)]
    D = distance_matrix(pts)
    print(f"配送先 {n} 拠点（1拠点目が営業所）")
    print(f"総経路数 = (n-1)!/2 = {math.factorial(n-1)//2:,} 通り\n")

    t0 = time.perf_counter()
    opt, opt_len = brute_force(D)
    t_bf = time.perf_counter() - t0
    print(f"[厳密解 総当たり] 距離 {opt_len:8.3f}  ({t_bf:.2f} 秒)")

    t0 = time.perf_counter()
    nn = nearest_neighbor(D)
    nn_len = tour_length(nn, D)
    print(f"[貪欲法 最近傍]   距離 {nn_len:8.3f}  "
          f"(厳密解比 +{nn_len/opt_len - 1:.2%}, {time.perf_counter()-t0:.5f} 秒)")

    t0 = time.perf_counter()
    imp = two_opt(nn, D)
    imp_len = tour_length(imp, D)
    print(f"[貪欲法 + 2-opt]  距離 {imp_len:8.3f}  "
          f"(厳密解比 +{imp_len/opt_len - 1:.2%}, {time.perf_counter()-t0:.5f} 秒)")

    t0 = time.perf_counter()
    best, best_len = None, float("inf")
    for s in range(n):                       # 出発点を変えたマルチスタート
        t = two_opt(nearest_neighbor(D, s), D)
        L = tour_length(t, D)
        if L < best_len:
            best, best_len = t, L
    print(f"[マルチスタート]  距離 {best_len:8.3f}  "
          f"(厳密解比 +{best_len/opt_len - 1:.2%}, {time.perf_counter()-t0:.5f} 秒)")

    print(f"\n最適経路: {' -> '.join(map(str, opt))} -> {opt[0]}")
    print(f"2-opt解 : {' -> '.join(map(str, best))} -> {best[0]}")

    print("\n[計算量の壁] 総当たりの経路数")
    for k in [8, 10, 12, 15, 20, 25]:
        c = math.factorial(k - 1) // 2
        secs = c / 5e6                       # 1秒500万経路と仮定
        unit = (f"{secs:.1f}秒" if secs < 60 else
                f"{secs/60:.1f}分" if secs < 3600 else
                f"{secs/3600:.1f}時間" if secs < 86400 * 365 else
                f"{secs/(86400*365):.3g}年")
        print(f"  n={k:>3}: {c:>25,} 通り  概算計算時間 {unit}")
    print("  ※ だからこそ実務ではヒューリスティクスや分枝限定法が使われる。")


if __name__ == "__main__":
    main()

実行結果

配送先 11 拠点（1拠点目が営業所）
総経路数 = (n-1)!/2 = 1,814,400 通り

[厳密解 総当たり] 距離  310.001  (6.38 秒)
[貪欲法 最近傍]   距離  319.693  (厳密解比 +3.13%, 0.00004 秒)
[貪欲法 + 2-opt]  距離  310.001  (厳密解比 +0.00%, 0.00004 秒)
[マルチスタート]  距離  310.001  (厳密解比 +0.00%, 0.00046 秒)

最適経路: 0 -> 3 -> 1 -> 9 -> 8 -> 5 -> 6 -> 7 -> 10 -> 2 -> 4 -> 0
2-opt解 : 0 -> 4 -> 2 -> 10 -> 7 -> 6 -> 5 -> 8 -> 9 -> 1 -> 3 -> 0

[計算量の壁] 総当たりの経路数
  n=  8:                     2,520 通り  概算計算時間 0.0秒
  n= 10:                   181,440 通り  概算計算時間 0.0秒
  n= 12:                19,958,400 通り  概算計算時間 4.0秒
  n= 15:            43,589,145,600 通り  概算計算時間 2.4時間
  n= 20:    60,822,550,204,416,000 通り  概算計算時間 386年
  n= 25: 310,224,200,866,619,719,680,000 通り  概算計算時間 1.97e+09年
  ※ だからこそ実務ではヒューリスティクスや分枝限定法が使われる。

← 目次 ｜ ← 前: 21 ｜ 次: 23 →