JijModeling + OMMX Adapterを用いた数理モデルの構築とOpenJijによる最適化計算#
このチュートリアルでは、JijModelingを用いて数理モデルを定式化し、得られた数理モデルをQUBOに変換し、OpenJijで解くという流れを説明したいと思います。
まず初めに、必要なパッケージをインストールします。
数理モデルを簡単に構築するためのモジュールであるJijModelingと、OpenJijの入出力フォーマットと相互に変換する ommx-openjij-adapter をインストールします。
これらは、pipを使ってインストールすることができます。
# !pip install jijmodeling ommx-openjij-adapter
巡回セールスマン問題#
制約条件付き最適化問題の例として巡回セールスマン問題を解いていきたいと思います。 巡回セールスマン問題は、一人のセールスマンが決められた都市を全て一度づつ訪問し、最終的に元の都市に帰ってくる時に、都市を巡回する最短経路を求めろという問題です。
制約条件#
この問題では、セールスマンは一つの地点に一度しか訪れることができないという位置に関する制約条件と、セールスマンが一人なのである時刻では一つの都市にしか存在しないという時間に関する制約条件が存在します。
都市\(i\)を\(t\)番目に訪れるとき\(x_{i,t}=1\)、それ以外では\(x_{i,t}=0\)とするバイナリ変数を用いると、上記の二つの制約条件は、
と書くことができます。
目的関数#
巡回セールスマン問題は、都市を巡回する最短経路を求めろという問題でした。 そこで、地点\(i\)と\(j\)の間の距離を\(d_{ij}\)とすると、時刻\(t\)で都市\(i\)を訪れ、時刻\(t+1\)で都市\(j\)を訪れた時の移動距離は、
と書くことができます。 これを合計したもの、
が今回最小化したい目的関数である、合計移動距離になります。
これまでのチュートリアルで述べたようにイジング最適化を行うためには、このような制約条件を持つ数理モデルをIsingハミルトニアンやQUBOハミルトニアンに変換する必要があります。 このような作業を手で行うと面倒ですし、実際に構築した数理モデルとQUBOの間にバグが入り込む可能性があります。 そこで、このような作業を全て自動でおこなってくれるのがJijModelingです。 JijModelingを用いることで、上記のように構築した数理モデルをコーディングし、それを自動的にQUBOに変換してくれます。 ここでは、上記で説明した巡回セールスマン問題を例にとって、JijModelingの使い方について説明していきます。
JijModelingによる定式化#
まず初めに、JijModelingを用いて、問題の数式を記述していきます。 JijModelingでは、通常の数理最適化計算用のモデラーとは異なり、問題インスタンスのデータとは独立に数理モデルを構築していきます。 このように数理モデルを構築することで、数理モデルの汎用性が担保でき、かつ、紙の上で数式を書くように直感的に数式を記述することができます。 さらに、JijModelingで記述された数理モデルは、notebook上ではLaTeXで確認することができます。
ここでは、JijModelingを用いた数理モデルの構築についてひとつづつ見ていきたいと思います。
まずは、問題を記述するための変数を定義しましょう。
import jijmodeling as jm
problem = jm.Problem("TSP")
N = problem.Natural("N")
d = problem.Float("d", shape=(N, N))
x = problem.BinaryVar("x", shape=(N, N))
ここで、jm.Problemで問題を定義しており、TSPという名前をつけています。
problem.Natural、problem.Floatは定数を表現しており、ここでは都市数\(N\)と距離行列\(d\)を表現するのに用いています。
巡回セールスマン問題においては、この距離行列と都市数によってさまざまな問題が表現されることになります。
バイナリ変数を表現するのが、problem.BinaryVarです。
ここでは、\(N\times N\)のバイナリ変数を定義しています。
JijModelingでは、jm.Problemに目的関数や制約条件を追加していきます。
目的関数#
次に、目的関数を定式化しましょう。
problem += jm.sum(
jm.product(N, N, N),
lambda i, j, t: d[i, j] * x[i, t] * x[j, (t + 1) % N],
)
総和はjm.sumを用いて表現することができます。
jm.sumの最初の引数は総和をとるドメインで、TSPの目的関数では3つの添字について総和を取るので、jm.product(N, N, N)で直積集合を指定し、lambda式で各添字に対する式を記述しています。
ここまでの数理モデルをLaTeX表示を通して確認してみましょう。
problem
制約条件#
次に、制約条件を追加していきましょう。
# 制約条件1 : 位置に関するonehot制約
problem += problem.Constraint(
"one-city",
lambda t: jm.sum(N, lambda i: x[i, t]) == 1, domain=N
)
# 制約条件2 : 時間に関するonehot制約
problem += problem.Constraint(
"one-time",
lambda i: jm.sum(N, lambda t: x[i, t]) == 1, domain=N
)
problem.Constraintを用いて制約条件を表現することができます。
最初の引数は制約条件の名前、2つ目の引数が制約条件を表現した数式になります。
ここでは "one-city" と名付けられた制約により、都市 \(i\) に訪れるのは1度だけを意味する制約を定義しています。
そして"one-time" と名付けられた制約によって、\(t\) 番目に訪れる都市は1つだけを意味する制約条件を設定しています。
制約条件を追加した数理モデルを確認してみましょう。
problem
説明で用いた数理モデルと同じ数式が表示されていることがわかります。
以上で、数理モデルの構築は終わりです。 このようにJijModelingを用いると、手元の数式と見比べながら数理モデルをコーディングしていくことができます。
インスタンスの作成#
問題の数理モデルができたので、次に問題に使うインスタンスデータを作成します。 ここでは、単純な都市数7で、都市間の距離をランダムにした問題を解いていきます。
import matplotlib.pyplot as plt
import numpy as np
inst_N = 7
np.random.seed(1)
x_pos = np.random.rand(inst_N)
y_pos = np.random.rand(inst_N)
plt.plot(x_pos, y_pos, 'o')
plt.xlim(0, 1)
plt.ylim(0, 1)
(0.0, 1.0)
数理モデルを作る際に用いたproblem.Naturalやproblem.Floatにデータを代入するので、変数の名前をkeyに持つ辞書でデータを渡す必要があります。
今回の問題では、Nとdに値を渡す必要があります。
XX, XX_T = np.meshgrid(x_pos, x_pos)
YY, YY_T = np.meshgrid(y_pos, y_pos)
inst_d = np.sqrt((XX - XX_T)**2 + (YY - YY_T)**2)
instance_data = {"N": inst_N, "d": inst_d}
instance_data
{'N': 7,
'd': array([[0. , 0.30759475, 0.45952133, 0.13629233, 0.43406434,
0.35402107, 0.58040301],
[0.30759475, 0. , 0.73408488, 0.41859314, 0.64201675,
0.6567735 , 0.71897319],
[0.45952133, 0.73408488, 0. , 0.32503125, 0.20721367,
0.34684999, 0.38700806],
[0.13629233, 0.41859314, 0.32503125, 0. , 0.30817754,
0.30035264, 0.4733741 ],
[0.43406434, 0.64201675, 0.20721367, 0.30817754, 0. ,
0.48383715, 0.19690151],
[0.35402107, 0.6567735 , 0.34684999, 0.30035264, 0.48383715,
0. , 0.6801809 ],
[0.58040301, 0.71897319, 0.38700806, 0.4733741 , 0.19690151,
0.6801809 , 0. ]])}
OMMX Instanceへの変換#
数理モデルとデータの用意ができたので、次にJijModelingで記述したモデルに先ほどのinstance_dataを代入してOMMX Instanceオブジェクトを作成します。
OMMX Instanceは数理最適化におけるモデラー<>ソルバー間、モデラー<>モデラー間のデータの受け渡しを行うための共通フォーマットであるOMMXで規定された数理モデルのインスタンスを表すデータ規格に則ったオブジェクトです。JijModelingとOMMXへの変換などはJijModelingのドキュメントまたはOMMXのドキュメントを参照してください。
# instance_dataをproblemに代入してommx instanceを取得
instance = problem.eval(instance_data)
この段階では instance はまだ制約付き最適化のデータです。簡単に中身を確認してみましょう。
instance.decision_variables_df.head(5)
| kind | lower | upper | name | subscripts | description | substituted_value | parameters.subscripts | |
|---|---|---|---|---|---|---|---|---|
| id | ||||||||
| 0 | Binary | 0.0 | 1.0 | x | [0, 0] | <NA> | <NA> | [0,0] |
| 1 | Binary | 0.0 | 1.0 | x | [0, 1] | <NA> | <NA> | [0,1] |
| 2 | Binary | 0.0 | 1.0 | x | [0, 2] | <NA> | <NA> | [0,2] |
| 3 | Binary | 0.0 | 1.0 | x | [0, 3] | <NA> | <NA> | [0,3] |
| 4 | Binary | 0.0 | 1.0 | x | [0, 4] | <NA> | <NA> | [0,4] |
instance.constraints_df
| equality | type | used_ids | name | subscripts | description | |
|---|---|---|---|---|---|---|
| id | ||||||
| 0 | =0 | Linear | {0, 35, 7, 42, 14, 21, 28} | one-city | [0] | <NA> |
| 1 | =0 | Linear | {1, 36, 8, 43, 15, 22, 29} | one-city | [1] | <NA> |
| 2 | =0 | Linear | {2, 37, 9, 44, 16, 23, 30} | one-city | [2] | <NA> |
| 3 | =0 | Linear | {3, 38, 10, 45, 17, 24, 31} | one-city | [3] | <NA> |
| 4 | =0 | Linear | {32, 4, 39, 11, 46, 18, 25} | one-city | [4] | <NA> |
| 5 | =0 | Linear | {33, 5, 40, 12, 47, 19, 26} | one-city | [5] | <NA> |
| 6 | =0 | Linear | {34, 6, 41, 13, 48, 20, 27} | one-city | [6] | <NA> |
| 7 | =0 | Linear | {0, 1, 2, 3, 4, 5, 6} | one-time | [0] | <NA> |
| 8 | =0 | Linear | {7, 8, 9, 10, 11, 12, 13} | one-time | [1] | <NA> |
| 9 | =0 | Linear | {14, 15, 16, 17, 18, 19, 20} | one-time | [2] | <NA> |
| 10 | =0 | Linear | {21, 22, 23, 24, 25, 26, 27} | one-time | [3] | <NA> |
| 11 | =0 | Linear | {32, 33, 34, 28, 29, 30, 31} | one-time | [4] | <NA> |
| 12 | =0 | Linear | {35, 36, 37, 38, 39, 40, 41} | one-time | [5] | <NA> |
| 13 | =0 | Linear | {42, 43, 44, 45, 46, 47, 48} | one-time | [6] | <NA> |
バイナリ変数と制約条件が付与されていることが確認できます。
では次に、OMMX InstanceをQUBOに変換します。OMMX Instanceはペナルティ法を用いてQUBOに変換されます。
instance.to_qubo()を実行することで、QUBOに変換されます。to_quboについては to_qubo API Referenceを参照してください。ペナルティ係数の重さを調整したい場合は、uniform_penalty_weightで全てのペナルティ係数を一括で調整することができます。もしくは上記のConstraint tableのidをキーにして制約条件に対応するペナルティに対してpenalty_weights: dict[int, float] 引数を指定することで、個別にペナルティ係数を調整することができます。指定しない場合は全てのペナルティ係数が1.0になります。
qubo, const = instance.to_qubo(uniform_penalty_weight=0.8)
OpenJijによる最適化計算の実行#
ここまでで得たQUBOを用いることで、これまでと同様に最適化計算を行うことができます。
from openjij import SASampler
sampler = SASampler()
res = sampler.sample_qubo(Q=qubo, num_reads=10)
OpenJijで最適化を行なっても01の羅列が帰ってくるので、これでは読みにくいです。ommx-openjij-adapter を用いてバイナリ列を用いて制約を満たしたかの確認と、その後使いやすい形に変換しましょう。
from ommx_openjij_adapter import OMMXOpenJijSAAdapter
# AdapterはOMMX InstanceとOpenJijの間の変換を担います
adapter = OMMXOpenJijSAAdapter(instance)
# OpenJijの結果をOMMX SampleSetに変換.
# このタイミングで元の最適化問題の目的関数の計算や制約条件を満たしているかを計算しています.
sampleset = adapter.decode_to_sampleset(res)
samplesetは複数の解をもつフォーマットです。この中から最良の解を選択します。
.best_feasibleはinstanceに入っている制約条件を満たす解 (実行可能解)の中で最も目的関数が小さいものを選びます。しかし instance.to_qubo とすると制約条件が全て削除され、ペナルティとして目的関数にまとめられているため、このメソッドで実行可能解を選ぶことはできません。
制約条件を削除してペナルティとして目的関数にまとめる操作を relax (緩和) と呼びます。つまりinstanceは.to_quboを呼んだことでrelaxされています。しかし元の制約条件の満たしているかを確認したいはずです。実はrelaxして制約条件が消されたように見えても実は内部では元の制約条件を別のところに確保しているため緩和する前の状態での実行可能解を抽出する.best_feasible_unrelaxed が用意されています。QUBOにした後はこのメソッドを用いて実行可能解からベストの目的関数値を持つ解を選ぶことができます。
best_sample = sampleset.best_feasible_unrelaxed
decision_variables_dfを用いて、OpenJijからの結果を元の数理モデルの変数に戻した見やすい形で取得することができます。
df = best_sample.decision_variables_df
nonzero_keys = {tuple(row["subscripts"]) for _, row in df.iterrows() if row["name"] == "x" and row["value"] > 0.5}
nonzero_keys
{(0, 4), (1, 3), (2, 6), (3, 2), (4, 1), (5, 5), (6, 0)}
sampleset.constraints_dfには元の制約条件とその値を確認することができます。value.には制約条件の左辺引く右辺の値が格納されています。制約条件を満たしている場合は0、満たしていない場合は0以外の値が入ります。feasible.はこの制約を満たしたかというbool値が格納されています。.best_feasible_unrelaxedを用いたので全ての制約条件が満たされていることが確認できます。
またremoved_reasonにはこの制約条件がuniform_penalty_methodでペナルティとして目的関数にまとめられたということが格納されています。
sampleset.constraints_df
| equality | used_ids | name | subscripts | description | removed_reason | value.0 | value.1 | value.2 | value.3 | ... | feasible.0 | feasible.1 | feasible.2 | feasible.3 | feasible.4 | feasible.5 | feasible.6 | feasible.7 | feasible.8 | feasible.9 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| id | |||||||||||||||||||||
| 0 | =0 | {0, 35, 21, 7, 42, 28, 14} | one-city | [0] | None | uniform_penalty_method | -1.0 | 0.0 | 0.0 | 0.0 | ... | False | True | True | True | True | True | True | True | True | True |
| 1 | =0 | {1, 36, 22, 8, 43, 29, 15} | one-city | [1] | None | uniform_penalty_method | 0.0 | -1.0 | 0.0 | 0.0 | ... | True | False | True | True | False | True | True | True | True | True |
| 2 | =0 | {16, 2, 37, 23, 9, 44, 30} | one-city | [2] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | False | True | True | True | True | True |
| 3 | =0 | {17, 3, 38, 24, 10, 45, 31} | one-city | [3] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | False | True | True | True | True | True |
| 4 | =0 | {32, 18, 4, 39, 25, 11, 46} | one-city | [4] | None | uniform_penalty_method | 0.0 | 1.0 | 0.0 | 0.0 | ... | True | False | True | True | True | True | True | True | True | True |
| 5 | =0 | {33, 19, 5, 40, 26, 12, 47} | one-city | [5] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | False | True | True | True | True | False |
| 6 | =0 | {48, 34, 20, 6, 41, 27, 13} | one-city | [6] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | True | True | True | False | True | False |
| 7 | =0 | {0, 1, 2, 3, 4, 5, 6} | one-time | [0] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | True | True | True | True | True | True |
| 8 | =0 | {7, 8, 9, 10, 11, 12, 13} | one-time | [1] | None | uniform_penalty_method | -1.0 | 0.0 | 0.0 | 0.0 | ... | False | True | True | True | True | False | True | True | True | True |
| 9 | =0 | {16, 17, 18, 19, 20, 14, 15} | one-time | [2] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | True | True | True | True | True | True |
| 10 | =0 | {21, 22, 23, 24, 25, 26, 27} | one-time | [3] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | True | False | True | True | True | True |
| 11 | =0 | {32, 33, 34, 28, 29, 30, 31} | one-time | [4] | None | uniform_penalty_method | 0.0 | 0.0 | -1.0 | 0.0 | ... | True | True | False | True | True | True | True | False | True | True |
| 12 | =0 | {35, 36, 37, 38, 39, 40, 41} | one-time | [5] | None | uniform_penalty_method | 0.0 | 0.0 | 0.0 | 0.0 | ... | True | True | True | True | True | True | True | True | True | True |
| 13 | =0 | {48, 42, 43, 44, 45, 46, 47} | one-time | [6] | None | uniform_penalty_method | 0.0 | 0.0 | 1.0 | 0.0 | ... | True | True | False | True | True | True | True | True | True | True |
14 rows × 26 columns
ではこの解を可視化してみましょう。
まずはnonzero_keysを用いて巡回路を作成します。
tour = [0 for _ in range(inst_N)]
for (i, t) in nonzero_keys:
tour[t] = i
tour.append(tour[0]) # 最後の都市から最初の都市へ戻る
tour
[6, 4, 3, 1, 0, 5, 2, 6]
plt.plot(x_pos, y_pos, 'o',markersize=12)
plt.plot(x_pos[tour], y_pos[tour], 'blue')
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.show()
ommx-openjij-adapter を用いたショートカット#
上記のチュートリアルでは quboへ変換 -> OpenJijで最適化 -> ommx-openjij-adapter でデコード/評価 という3つのステップを踏みましたが、実は ommx-openjij-adapter を用いることで、これらのステップを一気に行うことができます (現在はSimulated Annealingの場合のみ)。
adapter = OMMXOpenJijSAAdapter(instance)
sampleset = adapter.sample(instance)
このようにAdapter.sampleメソッドを用いることで一行で全ての流れを実行できます。またAdapter.solveというメソッドも提供されており、こちらではsampleで複数の回を得た後に最も良い解にフィルタリングされた結果を得ることができます。