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

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

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

巡回セールスマン問題（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{split} \begin{align*} \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{align*} \end{split}\]

これに対応する 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 = [
    (sum(x[t][i] for i in range(N)) == 1)
        .set_id(t)  # type: ignore
        .add_name("place")
        .add_subscripts([t])
    for t in range(N)
]
time_constraint = [
    (sum(x[t][i] for t in range(N)) == 1)
        .set_id(i + N)  # type: ignore
        .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

	objective	feasible
sample_id
3	78.299937	True
1	93.674431	True
13	93.945017	True
0	94.909303	True
11	96.279679	True
15	101.301923	True
2	102.480075	True
9	102.553352	True
14	109.891757	True
8	89.173259	False
6	101.764051	False
10	106.451972	False
4	112.338591	False
7	112.676912	False
12	117.568245	False
5	132.687222	False

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

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

sample_set.summary_with_constraints

	objective	feasible	place[0]	place[1]	place[2]	place[3]	place[4]	place[5]	place[6]	place[7]	...	time[6]	time[7]	time[8]	time[9]	time[10]	time[11]	time[12]	time[13]	time[14]	time[15]
sample_id
3	78.299937	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
1	93.674431	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
13	93.945017	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
0	94.909303	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
11	96.279679	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
15	101.301923	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
2	102.480075	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
9	102.553352	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
14	109.891757	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
8	89.173259	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
6	101.764051	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
10	106.451972	False	False	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
4	112.338591	False	True	True	True	True	True	True	False	True	...	True	True	True	True	False	True	True	True	True	True
7	112.676912	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
12	117.568245	False	True	True	True	False	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
5	132.687222	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True

16 rows × 34 columns

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

sample_set.decision_variables_df.head(2)

	kind	lower	upper	name	subscripts	description	0	1	2	3	...	6	7	8	9	10	11	12	13	14	15
id
0	Binary	0.0	1.0	x	[0, 0]	None	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0
1	Binary	0.0	1.0	x	[0, 1]	None	1.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

2 rows × 22 columns

sample_set.constraints_df.head(2)

	equality	used_ids	name	subscripts	description	removed_reason	value.0	value.1	value.2	value.3	...	feasible.6	feasible.7	feasible.8	feasible.9	feasible.10	feasible.11	feasible.12	feasible.13	feasible.14	feasible.15
id
0	=0	{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...	place	[0]	None	uniform_penalty_method	0.0	0.0	0.0	0.0	...	True	True	True	True	False	True	True	True	True	True
1	=0	{16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 2...	place	[1]	None	uniform_penalty_method	0.0	0.0	0.0	0.0	...	True	True	True	True	True	True	True	True	True	True

2 rows × 38 columns

得られたサンプルを取得するには 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)]

0.0

\(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

Index([3, 1, 13, 0, 11, 15, 2, 9, 14], dtype='int64', name='sample_id')

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

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()