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

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

ここでは巡回セールスマン問題を例として、問題をQUBOに変換しサンプリングを行う方法を説明します。

../_images/taraimawashi_businessman.png

たらい回しのイラスト（スーツ・男性）#

巡回セールスマン問題（TSP）は一人のセールスマンが複数の都市を順番に巡る方法を求める問題です。都市間の移動コストが与えられたときコストが最小になる経路を求めます。ここでは次の都市の配置を考えましょう

# From ulysses16.tsp in TSPLIB
ulysses16_points = [
    (38.24, 20.42),
    (39.57, 26.15),
    (40.56, 25.32),
    (36.26, 23.12),
    (33.48, 10.54),
    (37.56, 12.19),
    (38.42, 13.11),
    (37.52, 20.44),
    (41.23, 9.10),
    (41.17, 13.05),
    (36.08, -5.21),
    (38.47, 15.13),
    (38.15, 15.35),
    (37.51, 15.17),
    (35.49, 14.32),
    (39.36, 19.56),
]

都市の位置をプロットしてみましょう

%matplotlib inline
from matplotlib import pyplot as plt

x_coords, y_coords = zip(*ulysses16_points)
plt.scatter(x_coords, y_coords)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Ulysses16 Points')
plt.show()

コストとして単純に移動距離を考えましょう。\(i\)番目の都市と\(j\)番目の都市の距離 \(d(i, j)\)を計算しておきます。

def distance(x, y):
    return ((x[0] - y[0])**2 + (x[1] - y[1])**2)**0.5

# 都市の数
N = len(ulysses16_points)
# 各都市間の距離
d = [[distance(ulysses16_points[i], ulysses16_points[j]) for i in range(N)] for j in range(N)]

これを使って次のような最適化問題としてTSPを定式化します。まずある時刻 \(t\) に都市 \(i\) にいるかどうかをバイナリ変数 \(x_{t, i}\) で表します。このとき、以下の制約を満たすような \(x_{t, i}\) を求めます。するとセールスマンが移動する距離は次で与えられます：

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

ただし \(x_{t, i}\) は自由に取れるわけではなく、各時刻 \(t\) において一箇所の都市にしかいられないという制約と各都市について一度だけ訪れるという制約

\[ \sum_{i=0}^{N-1} x_{t, i} = 1, \quad \sum_{t=0}^{N-1} x_{t, i} = 1 \]

を満たす必要があります。これらを合わせてTSPは制約付き最適化問題として定式化できます

\[ \begin{aligned} \min \quad & \sum_{t=0}^{N-1} \sum_{i, j = 0}^{N-1} d(i, j) x_{t, i} x_{(t+1 \% N), j} \\ \text{s.t.} \quad & \sum_{i=0}^{N-1} x_{t, i} = 1 \quad (\forall t = 0, \ldots, N-1) \\ \quad & \sum_{t=0}^{N-1} x_{t, i} = 1 \quad (\forall i = 0, \ldots, N-1) \end{aligned} \]

これに対応する ommx.v1.Instance は次のように作成できます

from ommx.v1 import DecisionVariable, Instance

x = [[
        DecisionVariable.binary(
            i + N * t,  # 決定変数のID
            name="x",           # 決定変数の名前、解を取り出すときに使う
            subscripts=[t, i])  # 決定変数の添字、解を取り出すときに使う
        for i in range(N)
    ]
    for t in range(N)
]

objective = sum(
    d[i][j] * x[t][i] * x[(t+1) % N][j]
    for i in range(N)
    for j in range(N)
    for t in range(N)
)
place_constraint = {
    t: (sum(x[t][i] for i in range(N)) == 1)
        .add_name("place")
        .add_subscripts([t])
    for t in range(N)
}
time_constraint = {
    i + N: (sum(x[t][i] for t in range(N)) == 1)
        .add_name("time")
        .add_subscripts([i])
    for i in range(N)
}

instance = Instance.from_components(
    decision_variables=[x[t][i] for i in range(N) for t in range(N)],
    objective=objective,
    constraints={**place_constraint, **time_constraint},
    sense=Instance.MINIMIZE
)

バイナリの決定変数の作成時 DecisionVariable.binary に追加した決定変数の名前と添字は後で得られたサンプルを解釈する際に使います。

OpenJijによるサンプリング#

ommx.v1.Instance で記述されたQUBOをOpenJijを使ってサンプリングするには ommx-openjij-adapter を使います。

from ommx_openjij_adapter import OMMXOpenJijSAAdapter

sample_set = OMMXOpenJijSAAdapter.sample(instance, num_reads=16, uniform_penalty_weight=20.0)
sample_set.summary

OMMXOpenJijSAAdapter.solve は ommx.v1.SampleSet を返し、これはサンプルの値に加えて、目的関数の値や制約条件の破れを計算した値を保持しています。SampleSet.summary プロパティはこれらの要約情報を表示するためのものです。feasible はQUBOに変換する前の、元の問題に対する実行可能性を示しています。これは qubo インスタンスの removed_constraints に格納されている情報を使って計算されます。

各制約条件毎のfeasibilityを見るには summary_with_constraints プロパティを使います。

sample_set.summary_with_constraints

より詳しい情報は SampleSet.decision_variables_df 及び SampleSet.constraints_df プロパティを使って取得できます。

sample_set.decision_variables_df.head(2)

sample_set.constraints_df.head(2)

得られたサンプルを取得するには SampleSet.extract_decision_variables メソッドを使います。これは ommx.v1.DecisionVariables を作る時に登録した name と subscripts を使ってサンプルを解釈します。例えば sample_id=1 の x という名前の決定変数の値を取得するには次のようにすると dict[subscripts, value] の形で取得できます。

sample_id = 1
x = sample_set.extract_decision_variables("x", sample_id)
t = 2
i = 3
x[(t, i)]

\(x_{t, i}\)に対するサンプルが得れたのでこれをTSPのパスに変換します。これは今回の定式化自体に依存するので自分で処理を書く必要があります。

def sample_to_path(sample: dict[tuple[int, ...], float]) -> list[int]:
    path = []
    for t in range(N):
        for i in range(N):
            if sample[(t, i)] == 1:
                path.append(i)
    return path

これを表示してみましょう。まず元の問題に対してfeasibleであるサンプルのIDを取得します。

feasible_ids = sample_set.summary.query("feasible == True").index
feasible_ids

これらについて最適化された経路を表示しましょう

fig, axie = plt.subplots(3, 3, figsize=(12, 12))

for i, ax in enumerate(axie.flatten()):
    if i >= len(feasible_ids):
        break
    s = feasible_ids[i]
    x = sample_set.extract_decision_variables("x", s)
    path = sample_to_path(x)
    xs = [ulysses16_points[i][0] for i in path] + [ulysses16_points[path[0]][0]]
    ys = [ulysses16_points[i][1] for i in path] + [ulysses16_points[path[0]][1]]
    ax.plot(xs, ys, marker='o')
    ax.set_title(f"Sample {s}, objective={sample_set.objectives[s]:.2f}")

plt.tight_layout()
plt.show()
目次