OpenJij AdapterでTSPの経路をサンプリングする#

巡回セールスマン問題（TSP）をOMMXで表し、OpenJijの焼きなまし法で複数の経路候補を探す応用例です。都市の配置と、得られた実行可能な経路を描画します。

Instanceの構築、Preparation、SampleSetの読み方を知っている人を想定しています。基本操作を小さな例で確認したい場合は、OpenJij SA AdapterのReferenceを参照してください。以下では必要なデータとモデルをすべて作り直すため、別ページの実行状態は使いません。

実行環境とデータ#

共通のインストール手順とOpenJij Adapterの追加に従い、同じmainのPython SDKとAdapterを用意してください。経路の描画には、同じ環境にmatplotlibも追加します。

uv pip install matplotlib

TSPは、すべての都市を一度ずつ訪問して出発地へ戻る最短経路を求める問題です。ここでは、固定した乱数seedを使って10×10の領域に16都市を生成します。都市配置は再現できますが、サンプリングされる経路は実行ごとに変わる場合があります。外部データのダウンロードは不要です。

都市と距離を用意する#

from random import Random

CITY_COUNT = 16
random = Random(42)
city_points = [
    (random.uniform(0.0, 10.0), random.uniform(0.0, 10.0))
    for _ in range(CITY_COUNT)
]

%matplotlib inline
from matplotlib import pyplot as plt

fig, ax = plt.subplots(figsize=(7, 7))
city_x, city_y = zip(*city_points)
ax.scatter(city_x, city_y, s=60)
for city, point in enumerate(city_points):
    ax.annotate(str(city), point, xytext=(5, 5), textcoords="offset points")
ax.set(
    title="City locations",
    xlabel="x coordinate",
    ylabel="y coordinate",
)
ax.set_aspect("equal", adjustable="box")
ax.grid(alpha=0.2)
plt.show()

../_images/18820f642a08c89230f27ba9592b92d910a99faeb427f02c10364b61e08a9da5.png

都市\(i\)と都市\(j\)のユークリッド距離を\(d(i,j)\)とします。

def distance(left: tuple[float, float], right: tuple[float, float]) -> float:
    return ((left[0] - right[0]) ** 2 + (left[1] - right[1]) ** 2) ** 0.5


distances = [
    [distance(city_points[i], city_points[j]) for j in range(CITY_COUNT)]
    for i in range(CITY_COUNT)
]

経路の\(t\)番目に都市\(i\)を訪問するときに1となるBinary決定変数\(x_{t,i}\)を使います。目的関数は、出発地へ戻るまでの経路長です。

\[ \sum_{t=0}^{N-1} \sum_{i,j=0}^{N-1} d(i,j)x_{t,i}x_{(t+1) \bmod N,j}. \]

各位置では都市を1つだけ選び、すべての都市が1回ずつ現れる必要があります。

\[ \sum_{i=0}^{N-1}x_{t,i}=1 \quad (\forall t), \qquad \sum_{t=0}^{N-1}x_{t,i}=1 \quad (\forall i). \]

決定変数に付けた名前と添字はOMMXに保持され、後でsampleから経路を復元する際に使えます。

from ommx import DecisionVariable, Instance, Sense

route_variables = [
    [
        DecisionVariable.binary(
            city + CITY_COUNT * position,
            name="x",
            subscripts=[position, city],
        )
        for city in range(CITY_COUNT)
    ]
    for position in range(CITY_COUNT)
]

objective = sum(
    distances[i][j]
    * route_variables[position][i]
    * route_variables[(position + 1) % CITY_COUNT][j]
    for position in range(CITY_COUNT)
    for i in range(CITY_COUNT)
    for j in range(CITY_COUNT)
)

position_constraints = {
    position: (
        sum(route_variables[position][city] for city in range(CITY_COUNT)) == 1
    )
    .set_name("position")
    .add_subscripts([position])
    for position in range(CITY_COUNT)
}
city_constraints = {
    CITY_COUNT + city: (
        sum(route_variables[position][city] for position in range(CITY_COUNT))
        == 1
    )
    .set_name("city")
    .add_subscripts([city])
    for city in range(CITY_COUNT)
}

instance = Instance.from_components(
    decision_variables=[
        route_variables[position][city]
        for position in range(CITY_COUNT)
        for city in range(CITY_COUNT)
    ],
    objective=objective,
    constraints={**position_constraints, **city_constraints},
    sense=Sense.Minimize,
)

OpenJijによるサンプリング#

このモデルには制約があるため、そのままではAdapterの入力条件を満たしません。推奨Policyに、制約違反をサンプリング用の目的関数へ加えるpenaltyの重みを設定し、 prepare()で変換してから、追加変換をしないsample_without_preparation()へ渡します。 num_readsは読み出す候補の数、num_sweepsは各試行の探索量を調整する設定です。

from ommx import FixedPenaltyPreparation
from ommx_openjij_adapter import OMMXOpenJijSAAdapter

input_class = OMMXOpenJijSAAdapter.INPUT_CLASS

# OpenJij向けに推奨されるモデル変換から始めます。
policy = OMMXOpenJijSAAdapter.recommended_preparation_policy()

# Penaltyの大きさはapplication固有なので、呼び出し側で選びます。
policy.fixed_penalty = (
    FixedPenaltyPreparation.uniform_penalty_method_with_fixed_weight(
        weight=10.0,
    )
)
instance.prepare(input_class, policy)

# 各readが異なるランダムな軌道をたどれるよう、seedは指定しません。
sample_set = OMMXOpenJijSAAdapter.sample_without_preparation(
    instance,
    num_reads=32,
    num_sweeps=2000,
)

Instanceの変換Policy#

このAdapterのINPUT_CLASSは、制約がなく、使用する変数がBinaryで、目的関数が多項式の最小化問題です。多項式の次数を2次以下に限定してはいません。このTSPの例では、等式制約を penaltyへ変換すると2次以下の式になりますが、OpenJijを使うために必ずQUBO（制約なしBinary2次最小化問題）へ変換するという意味ではありません。

推奨Policyはこのモデルに必要な固定penaltyの重みを指定しないため、呼び出し側で設定します。重みはOpenJijの探索設定ではなく、OMMXによるモデル変換の設定です。すべてのモデルに適する共通の値ではありません。

この例では最大の辺コストが約9.5なので、同程度の10.0を出発点にしています。これは十分な重みを保証する公式ではありません。また、入力条件を満たすモデルへ変換できても、サンプラーが元の制約を満たす候補を返すとは限りません。結果のfeasibleと候補の多様性を確認しながら、モデルごとに重みや探索設定を調整してください。重みを変える場合は、モデルを作るセルから実行し直します。変換済みのInstanceはすでに入力条件を満たすため、Policyを変えてprepare()を呼ぶだけでは、適用済みの重みは置き換わりません。

prepare()は同じinstanceを更新します。編集したPolicyを適用済みなので、そのInstanceを sample_without_preparation()へ渡します。sample()は毎回新しい推奨Policyを取得するため、ここで編集したPolicyを引き継ぎません。一般的な仕組みはprepare()、変換の正確な条件はprepare()と uniform_penalty_method_with_fixed_weight()を参照してください。

結果の表示#

summary = sample_set.summary
summary.head(10)

	objective	feasible
sample_id
9	70.293540	True
11	70.349960	True
27	73.998265	True
28	77.165848	True
4	77.960204	True
31	78.769024	True
14	80.251552	True
2	81.954426	True
20	81.954426	True
22	81.954426	True

SampleSet.summaryは、各候補の目的関数値と実行可能性を持つpandas DataFrameです。 OMMXはこのpenalty変換で取り除いた制約も評価するため、feasibleから各候補がTSPの制約を満たすか確認できます。objectiveは元の目的関数の式で評価した値で、feasible=Trueの行では経路長を表します。固定penaltyはOpenJijの探索に使う目的値に影響しますが、この報告値には含まれません。この表はOpenJij内部のエネルギー（raw energy）を表示したものではありません。

表は実行可能な候補を先頭に、その中では経路長の短い順に並びます。得られた候補の中で最良でも、全経路の中で最短であることを保証するものではありません。

登録した変数名と添字を使うと、候補から\(x_{t,i}\)を取り出して経路へ戻せます。以下では SampleSet全体から最良の実行可能な候補を選びます。結果の読み方はSampleSetも参照してください。

def sample_to_route(sample: dict[tuple[int, ...], float]) -> list[int]:
    return [
        next(
            city
            for city in range(CITY_COUNT)
            if sample[(position, city)] > 0.5
        )
        for position in range(CITY_COUNT)
    ]


if not sample_set.feasible_ids():
    raise RuntimeError(
        "Feasibleな経路が得られませんでした。read数を増やすかpenaltyを見直してください。"
    )

best_sample_id = sample_set.best_feasible_id
best_objective = sample_set.objectives[best_sample_id]
best_sample = sample_set.extract_decision_variables("x", best_sample_id)
sampled_route = sample_to_route(best_sample)
sampled_route

[9, 0, 11, 1, 12, 2, 13, 3, 15, 4, 6, 5, 8, 7, 10, 14]

route_cycle = sampled_route + [sampled_route[0]]

fig, ax = plt.subplots(figsize=(7, 7))
route_x = [city_points[city][0] for city in route_cycle]
route_y = [city_points[city][1] for city in route_cycle]
ax.plot(route_x, route_y, color="tab:blue", alpha=0.6)
ax.scatter(city_x, city_y, s=60, color="tab:blue")

for start, end in zip(route_cycle, route_cycle[1:]):
    ax.annotate(
        "",
        xy=city_points[end],
        xytext=city_points[start],
        arrowprops={"arrowstyle": "->", "color": "tab:blue", "alpha": 0.7},
    )
for city, point in enumerate(city_points):
    ax.annotate(str(city), point, xytext=(5, 5), textcoords="offset points")

ax.set(
    title=(
        f"Best feasible sampled route\n"
        f"sample_id={best_sample_id}, distance={best_objective:.2f}"
    ),
    xlabel="x coordinate",
    ylabel="y coordinate",
)
ax.set_aspect("equal", adjustable="box")
ax.grid(alpha=0.2)
plt.show()

../_images/11b4697ee3637038fa3f0d3b93da8588bf0f4082d5998613047203cafa41b739.png

OpenJij AdapterでTSPの経路をサンプリングする

目次