巡回セールスマン問題 ― 距離が最小の回路を探す

巡回セールスマン問題 ― 距離が最小の回路を探す#

この章で学ぶこと#

すべての都市を1回ずつ回って出発点に戻る、一番短い経路を探します。

  • add_circuit() で「1本の巡回ルート」を作る

  • model.minimize() で距離の合計を最小にする

  • 目的関数の改善をコールバックで記録する

  • 制約を足して「線分を必ず通る」など、変わった TSP も解く

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

複数の都市を1回ずつ訪れて出発点に戻るとき、移動距離が最小になる順番を求める問題です。

TSP は NP困難で、都市数 \(N\) が増えると経路の数は \((N-1)!/2\) のように爆発します。

都市数

経路の数（およそ）

10

18万

15

4億

20

1.2兆

総当たりでは太刀打ちできないので、ソルバーの出番です。CP-SAT では「各都市から1本ずつ出て、1本ずつ入る」という回路の制約と、距離の最小化を組み合わせます。

from matplotlib import pyplot as plt
import numpy as np
from scipy import spatial
from matplotlib.collections import LineCollection
from ortools.sat.python import cp_model

テストデータ#

ランダムに100点を作り、すべての点の間の距離を計算します。CP-SAT は整数しか扱えないので、距離は1000倍して整数にします。

n = 100
np.random.seed(42)
points = np.random.uniform(0, 100, size=(n, 2))

# 距離行列（1000倍して整数にする）
dist_matrix = (spatial.distance.squareform(spatial.distance.pdist(points)) * 1000).astype(np.int32)

fig, ax = plt.subplots()
ax.scatter(points[:, 0], points[:, 1]);
_images/509d12f031dd3c68589d18f64a0540b01b0b1fb8471dec4c5600553dfd55b585.png

改善の履歴を記録する#

CpSolverSolutionCallback を継承すると、解が見つかるたびに呼ばれる関数を書けます。ここでは目的関数（総距離）の推移を記録します。

class ObjectiveTracker(cp_model.CpSolverSolutionCallback):
    """探索中に見つかった目的関数の値を記録する"""

    def __init__(self):
        super().__init__()
        self.history = []

    def on_solution_callback(self):
        self.history.append(self.objective_value)

解く#

solve_points_tsp() の流れは次のとおりです。

  1. 距離行列から、(i, j) → 距離 の辞書 edges を作る（i == j は除く）

  2. 各辺にブール変数を作る（真ならその道を通る）

  3. add_circuit() で「1本の閉路」を強制する

  4. LinearExpr.weighted_sum() で総距離を計算し、model.minimize() で最小化する

  5. 見つかった辺から、経路の座標を取り出す

max_time を過ぎると探索を打ち切ります。そのときは「最適とは限らないが、実行可能な解」が返ります（状態が OPTIMAL か FEASIBLE かで分かります）。

def solve_points_tsp(points, dist_matrix, log_history=False, max_time=10):
    """点の集合を回る最短の巡回ルートを求める"""
    edges = {}
    for (i, j), d in np.ndenumerate(dist_matrix):
        if i != j:
            edges[i, j] = int(d)

    model = cp_model.CpModel()
    variables = {key: model.new_bool_var(f"e{i}") for i, key in enumerate(edges)}
    model.add_circuit([(s, e, v) for (s, e), v in variables.items()])

    # 通る辺の距離の合計を最小にする
    length = cp_model.LinearExpr.weighted_sum(list(variables.values()), list(edges.values()))
    model.minimize(length)

    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = max_time
    if log_history:
        callback = ObjectiveTracker()
        status = solver.solve(model, callback)
        print("status:", solver.status_name(status))
        print("距離の推移:", callback.history)
    else:
        status = solver.solve(model)
        print("status:", solver.status_name(status))

    selected = points[np.array([key for key, v in variables.items() if solver.value(v)])]
    return selected

100点で解いて、選ばれた辺を線で結びます。

selected_edges = solve_points_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(points[:, 0], points[:, 1])
ax.add_collection(LineCollection(selected_edges));

status: OPTIMAL
距離の推移: [5094854.0, 779053.0, 766522.0, 766177.0, 765060.0, 764047.0, 763658.0, 759054.0, 749670.0, 747664.0, 745684.0]
_images/2491c9fde3465e6db7b76f64883ea16c6e11a4c39f6278a1a8fe0f2940263c1a.png

線分の TSP#

「点」ではなく「線分」を巡回する問題を考えます。50本の線分があり、各線分のどちらかの端から入って、もう一方の端から出る、という制約を足します。

まず線分を作ります。

n = 50
np.random.seed(42)
start_points = np.random.uniform(0, 100, size=(n, 2))
length = np.random.uniform(5, 10, size=n)
angle = np.random.uniform(0, 2 * np.pi, size=n)
end_points = start_points + np.c_[length * np.cos(angle), length * np.sin(angle)]
segments = np.column_stack([start_points, end_points])
points = segments.reshape(-1, 2)

fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2)))

<matplotlib.collections.LineCollection at 0x1f0b7154cd0>
_images/683b1b657fe2a87859c236faf2ef7a8c3225dbf292a8cda3bde040f1b31accd6.png

制約が弱いと、線分を通らない#

まず「同じ線分の両端の距離を 0 にする」だけを試します。線分の上を移動するコストは 0 になりますが、それでも全部の線分を通るとは限りません。コスト 0 の道を使わずに、近道をしたほうが総距離が小さくなるからです。

次の図では、赤い線が元の線分、黒い線が選ばれたルートです。通っていない線分があることを確認してください。

dist_matrix = (spatial.distance.squareform(spatial.distance.pdist(segments.reshape(-1, 2))) * 1000).astype(np.int32)
start_index = np.arange(0, n * 2, 2)
end_index = start_index + 1
# 線分の両端の距離を 0 にする
dist_matrix[start_index, end_index] = 0
dist_matrix[end_index, start_index] = 0

selected_edges = solve_points_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2), linewidth=5, color='red', alpha=0.6))
ax.add_collection(LineCollection(selected_edges, linewidth=2))

status: FEASIBLE
距離の推移: [5215210.0, 5208664.0, 947648.0, 930508.0, 746063.0, 691436.0, 527524.0, 461945.0]

<matplotlib.collections.LineCollection at 0x1f0b5fc2fd0>
_images/707604b3a1628b9d09c68417439abbf74fc104237ecbc05c402a7968b93d2233.png

「必ず通る」を制約で足す#

線分 (s, e) について、s → e か e → s のどちらかを必ず通るように add_exactly_one() で制約します。これで、すべての線分を通るルートになります。

def solve_segments_tsp(points, dist_matrix, log_history=False, max_time=10):
    """すべての線分を通る最短ルートを求める"""
    edges = {}
    for (i, j), d in np.ndenumerate(dist_matrix):
        if i != j:
            edges[i, j] = int(d)

    n = dist_matrix.shape[0]
    start_index = np.arange(0, n, 2)
    end_index = start_index + 1

    model = cp_model.CpModel()
    variables = {key: model.new_bool_var(f"e{i}") for i, key in enumerate(edges)}
    model.add_circuit([(s, e, v) for (s, e), v in variables.items()])

    # 線分ごとに、どちらかの向きを必ず通る
    for s, e in zip(start_index, end_index):
        model.add_exactly_one([variables[s, e], variables[e, s]])

    length = cp_model.LinearExpr.weighted_sum(list(variables.values()), list(edges.values()))
    model.minimize(length)

    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = max_time
    if log_history:
        callback = ObjectiveTracker()
        status = solver.solve(model, callback)
        print("status:", solver.status_name(status))
        print("距離の推移:", callback.history)
    else:
        status = solver.solve(model)
        print("status:", solver.status_name(status))

    selected = points[np.array([key for key, v in variables.items() if solver.value(v)])]
    return selected

これで、すべての線分を通るルートが得られます。

selected_edges = solve_segments_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2), linewidth=5, color='red', alpha=0.6))
ax.add_collection(LineCollection(selected_edges, linewidth=2))

status: FEASIBLE
距離の推移: [949519.0, 827425.0, 698362.0, 621424.0, 618375.0, 617604.0, 550116.0, 490336.0, 488825.0, 486431.0]

<matplotlib.collections.LineCollection at 0x1f0b7285bd0>
_images/1b5db3123756d8678afccc4149932b61c36e5bf6dc6b7faad576bc01de549ced.png

まとめ#

  • TSP は「add_circuit() ＋ minimize(距離の合計)」で書ける

  • 距離は整数にする（1000倍など）。CP-SAT は小数を扱えない

  • 解の改善を追いたいときは CpSolverSolutionCallback を継承する

  • max_time_in_seconds で時間を区切れる。時間切れのときの解は FEASIBLE（最適とは限らない）

  • 「必ず通る」などの条件は、辺の変数に add_exactly_one などを足すだけで表現できる

  • 制約が弱いと、意図した経路にならない（線分の例）。条件は明示的に書くこと