OMMX AdapterでQUBOからサンプリングする#

このチュートリアルでは、巡回セールスマン問題（TSP）を例に、OMMXのモデルをOpenJijでサンプリングするまでの一連の流れを説明します。都市の配置とサンプリングされた経路も描画し、得られた結果をその場で確認します。

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/d0fe16ebf6a41beb14144702e8b4be56fb0b3296c8da573e59476605b27ad010.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によるサンプリング#

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#

OpenJijが直接受け付けるのは、制約なしのBinary最小化modelです。このTSP modelには制約があり、penalty magnitudeには安全な共通defaultがありません。そのためAdapterの推奨変換 Policyをcustomizeし、Instanceに対してprepare()を呼び出します。推奨Policyでは、特殊制約のlowering、 optimization senseの正規化、Integer slackの追加、Integer変数のencodingなど、OpenJijで一般的に必要となる変換を有効にします。

制約を取り除くには固定penaltyが必要です。その大きさはOpenJij samplerのparameterではなく、 OMMXによるmodel変換の設定であり、すべてのmodelに安全な値はありません。大きな値ほど feasibleなsampleを得やすくなりますが、有限の値でfeasibilityが保証されるわけではありません。

この例では最大の辺コストが約9.5なので、同程度の10.0を出発点にしています。これは十分なpenaltyを保証する公式ではありません。feasibilityとsampleの多様性を確認しながら、 modelごとに値を調整してください。samplingは確率的なので、notebookを再実行すると異なる経路が得られる場合があります。

prepare()はinstanceをin-placeで更新します。後の変換が失敗した場合、それより前に完了した変換はinstanceに残ります。その後の preparation-freeな sample_without_preparation()呼び出しでは、 InstanceがexactなOpenJij inputであることを検査します。

結果の表示#

summary = sample_set.summary
summary.head(10)

	objective	feasible
sample_id
4	48.414094	True
2	65.546568	True
13	65.712669	True
6	65.823287	True
3	66.033051	True
8	66.058482	True
0	66.106397	True
5	67.296505	True
1	68.817908	True
16	70.114236	True

SampleSet.summary は、各sampleの目的関数値とfeasibilityを持つpandas DataFrameです。 OMMXはpreparation中に取り除かれた制約を保持するため、feasible から各sampleがTSPの制約を満たすか確認できます。objectiveはすべての行で元の経路長です。固定penaltyはOpenJijの sampling energyに影響しますが、この列には含まれません。表はfeasibleなsampleを先頭に、その中では目的関数値の順に並びます。

登録した変数名と添字を使うと、sampleから\(x_{t,i}\)を取り出して経路へ戻せます。ここでは summaryに表示された中から最良のfeasible sampleを選びます。

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

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

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/65bc47906a0b63a1099f137d40af51dac5ec4ea6fa4e388d5f978c418a006618.png

OMMX AdapterでQUBOからサンプリングする

目次

OMMX AdapterでQUBOからサンプリングする#

OpenJijによるサンプリング#

Instanceの変換Policy#

結果の表示#