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
6	85.937163	True
15	87.742206	True
0	92.420147	True
8	92.738861	True
4	93.314939	True
11	94.043606	True
12	94.060986	True
1	102.181668	True
14	102.867245	True
3	93.482688	False
9	98.658943	False
7	101.748679	False
5	105.570150	False
2	107.455912	False
10	109.686372	False
13	114.116205	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
6	85.937163	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
15	87.742206	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
0	92.420147	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
8	92.738861	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
4	93.314939	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
11	94.043606	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
12	94.060986	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
1	102.181668	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
14	102.867245	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
3	93.482688	False	True	True	True	True	False	True	True	True	...	True	True	True	True	False	True	True	True	True	True
9	98.658943	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
7	101.748679	False	True	True	True	False	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
5	105.570150	False	True	True	True	True	True	False	True	True	...	True	True	True	True	False	True	True	True	True	True
2	107.455912	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
10	109.686372	False	True	True	True	False	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
13	114.116205	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 及び SampleSet.constraints プロパティを使って取得できます。

sample_set.decision_variables.head(2)

	kind	lower	upper	name	subscripts	description	substituted_value	10	15	0	...	4	5	6	7	8	9	11	12	13	14
id
0	binary	0.0	1.0	x	[0, 0]	<NA>	<NA>	1.0	1.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
16	binary	0.0	1.0	x	[1, 0]	<NA>	<NA>	0.0	0.0	0.0	...	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0

2 rows × 23 columns

sample_set.constraints.head(2)

	equality	used_ids	name	subscripts	description	removed_reason	value.0	value.1	value.2	value.3	...	feasible.14	feasible.6	feasible.15	feasible.12	feasible.5	feasible.0	feasible.8	feasible.13	feasible.1	feasible.4
id
0	=0	{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...	place	[0]	<NA>	uniform_penalty_method	0.0	0.0	0.0	0.0	...	True	True	True	True	True	True	True	True	True	True
1	=0	{16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 2...	place	[1]	<NA>	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([6, 15, 0, 8, 4, 11, 12, 1, 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()