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
12	87.252332	True
6	90.192938	True
11	92.433671	True
1	100.972104	True
13	102.425559	True
15	103.683133	True
2	104.730022	True
4	107.617958	True
10	122.417866	True
3	92.606869	False
9	99.139406	False
14	105.959348	False
8	111.582674	False
5	112.176961	False
0	116.376136	False
7	120.110660	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
12	87.252332	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
6	90.192938	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
11	92.433671	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
1	100.972104	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
13	102.425559	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
15	103.683133	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
2	104.730022	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
4	107.617958	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
10	122.417866	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
3	92.606869	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
9	99.139406	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
14	105.959348	False	True	True	True	True	True	True	True	False	...	True	True	True	True	False	True	True	True	True	True
8	111.582674	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
5	112.176961	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
0	116.376136	False	False	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
7	120.110660	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	0.0	0.0	0.0	0.0	0.0
1	Binary	0.0	1.0	x	[0, 1]	None	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.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	-1.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]	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([12, 6, 11, 1, 13, 15, 2, 4, 10], 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()