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

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

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

../_images/taraimawashi_businessman.png — 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()

../_images/7f4dc275af06d3e6dae6479969a286991bfc44cfd951418ee0e8b732b9db3189.png

コストとして単純に移動距離を考えましょう。 $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, \sum_{t = 0}^{N - 1} x_{t, i} = 1

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

\begin{array}{r} \begin{aligned} min & \sum_{t = 0}^{N - 1} \sum_{i, j = 0}^{N - 1} d (i, j) x_{t, i} x_{(t + 1 % N), j} \\ s.t. & \sum_{i = 0}^{N - 1} x_{t, i} = 1 (\forall t = 0, \dots, N - 1) \\ \sum_{t = 0}^{N - 1} x_{t, i} = 1 (\forall i = 0, \dots, N - 1) \end{aligned} \end{array}

これに対応する 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
5	93.532998	True
1	95.126642	True
7	95.590577	True
13	95.861402	True
10	102.935178	True
2	105.138658	True
12	108.956331	True
11	115.814312	True
9	89.359943	False
6	100.341992	False
4	101.166567	False
15	101.252804	False
0	107.179374	False
14	109.782167	False
3	111.387205	False
8	121.947213	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
5	93.532998	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
1	95.126642	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
7	95.590577	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
13	95.861402	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
10	102.935178	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
2	105.138658	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
12	108.956331	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
11	115.814312	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
9	89.359943	False	True	True	True	True	True	True	True	False	...	True	True	True	True	False	True	True	True	True	True
6	100.341992	False	True	True	False	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
4	101.166567	False	False	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
15	101.252804	False	True	True	True	True	True	False	True	True	...	True	True	True	True	False	True	True	True	True	True
0	107.179374	False	True	True	True	False	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
14	109.782167	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
3	111.387205	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
8	121.947213	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	0	1	2	...	6	7	8	9	10	11	12	13	14	15
id
0	binary	0.0	1.0	x	[0, 0]	<NA>	<NA>	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
16	binary	0.0	1.0	x	[1, 0]	<NA>	<NA>	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	1.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.0	feasible.5	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	False
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([5, 1, 7, 13, 10, 2, 12, 11], 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()

../_images/9d93c2eb602cceaf18ba3de243b2ec93c79d8690eefa92dbc5b41831a3533b54.png

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

Contents

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

OpenJijによるサンプリング#