JijModelingを用いた数理モデルの構築とOpenJijでの最適化計算#
このチュートリアルでは、JijModelingを用いて数理モデルを定式化し、得られた数理モデルをQUBOに変換し、OpenJijで解くという流れを説明したいと思います。
まず初めに、必要なパッケージをインストールします。
数理モデルを簡単に構築するためのモジュールであるJijModelingとJijModelingで記述された数理モデルをQUBOに変換するためのモジュールであるjijmodeling-transpilerをインストールします。
これらは、pip
を使ってインストールすることができます。
!pip install jijmodeling
!pip install jijmodeling-transpiler
Collecting jijmodeling
Downloading jijmodeling-1.4.1-cp39-cp39-manylinux_2_28_x86_64.whl.metadata (12 kB)
Requirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling) (1.26.4)
Collecting pandas (from jijmodeling)
Downloading pandas-2.2.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (19 kB)
Collecting orjson<4.0.0,>=3.8.0 (from jijmodeling)
Downloading orjson-3.10.5-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (49 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 49.7/49.7 kB 16.4 MB/s eta 0:00:00
?25hRequirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pandas->jijmodeling) (2.9.0.post0)
Collecting pytz>=2020.1 (from pandas->jijmodeling)
Downloading pytz-2024.1-py2.py3-none-any.whl.metadata (22 kB)
Collecting tzdata>=2022.7 (from pandas->jijmodeling)
Downloading tzdata-2024.1-py2.py3-none-any.whl.metadata (1.4 kB)
Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from python-dateutil>=2.8.2->pandas->jijmodeling) (1.16.0)
Downloading jijmodeling-1.4.1-cp39-cp39-manylinux_2_28_x86_64.whl (2.0 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.0/2.0 MB 97.5 MB/s eta 0:00:00
?25hDownloading orjson-3.10.5-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (144 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 144.8/144.8 kB 42.9 MB/s eta 0:00:00
?25hDownloading pandas-2.2.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (13.1 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 13.1/13.1 MB 89.8 MB/s eta 0:00:00
?25hDownloading pytz-2024.1-py2.py3-none-any.whl (505 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 505.5/505.5 kB 78.8 MB/s eta 0:00:00
?25hDownloading tzdata-2024.1-py2.py3-none-any.whl (345 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 345.4/345.4 kB 68.2 MB/s eta 0:00:00
?25hInstalling collected packages: pytz, tzdata, orjson, pandas, jijmodeling
Successfully installed jijmodeling-1.4.1 orjson-3.10.5 pandas-2.2.2 pytz-2024.1 tzdata-2024.1
Collecting jijmodeling-transpiler
Downloading jijmodeling_transpiler-0.6.15-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl.metadata (3.9 kB)
Requirement already satisfied: jijmodeling<2.0.0,>=1.0.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (1.4.1)
Requirement already satisfied: numpy<1.27.0,>=1.17.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (1.26.4)
Collecting typeguard (from jijmodeling-transpiler)
Downloading typeguard-4.3.0-py3-none-any.whl.metadata (3.7 kB)
Collecting pydantic (from jijmodeling-transpiler)
Downloading pydantic-2.7.4-py3-none-any.whl.metadata (109 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 109.4/109.4 kB 22.5 MB/s eta 0:00:00
?25hCollecting pyqubo (from jijmodeling-transpiler)
Downloading pyqubo-1.4.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (7.3 kB)
Collecting mip (from jijmodeling-transpiler)
Downloading mip-1.15.0-py3-none-any.whl.metadata (21 kB)
Requirement already satisfied: dimod in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (0.12.15)
Requirement already satisfied: pandas in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling<2.0.0,>=1.0.0->jijmodeling-transpiler) (2.2.2)
Requirement already satisfied: orjson<4.0.0,>=3.8.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from jijmodeling<2.0.0,>=1.0.0->jijmodeling-transpiler) (3.10.5)
Collecting cffi==1.15.* (from mip->jijmodeling-transpiler)
Downloading cffi-1.15.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (1.1 kB)
Collecting pycparser (from cffi==1.15.*->mip->jijmodeling-transpiler)
Downloading pycparser-2.22-py3-none-any.whl.metadata (943 bytes)
Collecting annotated-types>=0.4.0 (from pydantic->jijmodeling-transpiler)
Downloading annotated_types-0.7.0-py3-none-any.whl.metadata (15 kB)
Collecting pydantic-core==2.18.4 (from pydantic->jijmodeling-transpiler)
Downloading pydantic_core-2.18.4-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (6.5 kB)
Requirement already satisfied: typing-extensions>=4.6.1 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pydantic->jijmodeling-transpiler) (4.12.2)
Collecting dwave-neal>=0.5.7 (from pyqubo->jijmodeling-transpiler)
Downloading dwave_neal-0.6.0-py3-none-any.whl.metadata (3.0 kB)
Collecting Deprecated>=1.2.12 (from pyqubo->jijmodeling-transpiler)
Downloading Deprecated-1.2.14-py2.py3-none-any.whl.metadata (5.4 kB)
Requirement already satisfied: six>=1.15.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pyqubo->jijmodeling-transpiler) (1.16.0)
Requirement already satisfied: importlib-metadata>=3.6 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from typeguard->jijmodeling-transpiler) (7.2.0)
Collecting wrapt<2,>=1.10 (from Deprecated>=1.2.12->pyqubo->jijmodeling-transpiler)
Downloading wrapt-1.16.0-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (6.6 kB)
Collecting dwave-samplers<2.0.0,>=1.0.0 (from dwave-neal>=0.5.7->pyqubo->jijmodeling-transpiler)
Downloading dwave_samplers-1.2.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (9.8 kB)
Requirement already satisfied: zipp>=0.5 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from importlib-metadata>=3.6->typeguard->jijmodeling-transpiler) (3.19.2)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pandas->jijmodeling<2.0.0,>=1.0.0->jijmodeling-transpiler) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pandas->jijmodeling<2.0.0,>=1.0.0->jijmodeling-transpiler) (2024.1)
Requirement already satisfied: tzdata>=2022.7 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from pandas->jijmodeling<2.0.0,>=1.0.0->jijmodeling-transpiler) (2024.1)
Collecting networkx>=2.4.0 (from dwave-samplers<2.0.0,>=1.0.0->dwave-neal>=0.5.7->pyqubo->jijmodeling-transpiler)
Downloading networkx-3.2.1-py3-none-any.whl.metadata (5.2 kB)
Downloading jijmodeling_transpiler-0.6.15-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl (7.6 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.6/7.6 MB 104.4 MB/s eta 0:00:00
?25hDownloading mip-1.15.0-py3-none-any.whl (15.3 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 15.3/15.3 MB 101.3 MB/s eta 0:00:00
?25hDownloading cffi-1.15.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (441 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 441.2/441.2 kB 76.3 MB/s eta 0:00:00
?25hDownloading pydantic-2.7.4-py3-none-any.whl (409 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 409.0/409.0 kB 75.3 MB/s eta 0:00:00
?25hDownloading pydantic_core-2.18.4-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (2.0 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.0/2.0 MB 106.5 MB/s eta 0:00:00
?25hDownloading pyqubo-1.4.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (245 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 245.5/245.5 kB 58.4 MB/s eta 0:00:00
?25hDownloading typeguard-4.3.0-py3-none-any.whl (35 kB)
Downloading annotated_types-0.7.0-py3-none-any.whl (13 kB)
Downloading Deprecated-1.2.14-py2.py3-none-any.whl (9.6 kB)
Downloading dwave_neal-0.6.0-py3-none-any.whl (8.7 kB)
Downloading dwave_samplers-1.2.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (6.7 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 6.7/6.7 MB 107.6 MB/s eta 0:00:00
?25hDownloading wrapt-1.16.0-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (80 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 80.1/80.1 kB 24.5 MB/s eta 0:00:00
?25hDownloading pycparser-2.22-py3-none-any.whl (117 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 117.6/117.6 kB 36.0 MB/s eta 0:00:00
?25hDownloading networkx-3.2.1-py3-none-any.whl (1.6 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.6/1.6 MB 104.2 MB/s eta 0:00:00
?25hInstalling collected packages: wrapt, pydantic-core, pycparser, networkx, annotated-types, typeguard, pydantic, dwave-samplers, Deprecated, cffi, mip, dwave-neal, pyqubo, jijmodeling-transpiler
Successfully installed Deprecated-1.2.14 annotated-types-0.7.0 cffi-1.15.1 dwave-neal-0.6.0 dwave-samplers-1.2.0 jijmodeling-transpiler-0.6.15 mip-1.15.0 networkx-3.2.1 pycparser-2.22 pydantic-2.7.4 pydantic-core-2.18.4 pyqubo-1.4.0 typeguard-4.3.0 wrapt-1.16.0
巡回セールスマン問題#
制約条件付き最適化問題の例として巡回セールスマン問題を解いていきたいと思います。 巡回セールスマン問題は、一人のセールスマンが決められた都市を全て一度づつ訪問し、最終的に元の都市に帰ってくる時に、都市を巡回する最短経路を求めろという問題です。
制約条件#
この問題では、セールスマンは一つの地点に一度しか訪れることができないという位置に関する制約条件と、セールスマンが一人なのである時刻では一つの都市にしか存在しないという時間に関する制約条件が存在します。
番目に都市を訪れるとき、それ以外ではとするバイナリ変数を用いると、上記の二つの制約条件は、
と書くことができます。
目的関数#
巡回セールスマン問題は、都市を巡回する最短経路を求めろという問題でした。 そこで、地点との間の距離をとすると、時刻で都市を訪れ、時刻で都市を訪れた時の移動距離は、
と書くことができます。 これを合計したもの、
が今回最小化したい目的関数である、合計移動距離になります。
これまでのチュートリアルで述べたようにイジング最適化を行うためには、このような制約条件を持つ数理モデルをIsingハミルトニアンやQUBOハミルトニアンに変換する必要があります。 このような作業を手で行うと面倒ですし、実際に構築した数理モデルとQUBOの間にバグが入り込む可能性があります。 そこで、このような作業を全て自動でおこなってくれるのがJijModelingです。 JijModelingを用いることで、上記のように構築した数理モデルをコーディングし、それを自動的にQUBOに変換してくれます。 ここでは、上記で説明した巡回セールスマン問題を例にとって、JijModelingの使い方について説明していきます。
JijModelingを用いた巡回セールスマン問題の数理モデルの構築#
まず初めに、JijModelingを用いて、問題の数式を記述していきます。 JijModelingでは、通常の数理最適化計算用のモデラーとは異なり、問題インスタンスのデータとは独立に数理モデルを構築していきます。 このように数理モデルを構築することで、数理モデルの汎用性が担保でき、かつ、紙の上で数式を書くように直感的に数式を記述することができます。 さらに、JijModelingで記述された数理モデルは、notebook上ではLaTeXで確認することができます。
ここでは、JijModelingを用いた数理モデルの構築についてひとつづつ見ていきたいと思います。
まずは、問題を記述するための変数と定数を表現しましょう。
import jijmodeling as jm
dist = jm.Placeholder("dist", ndim=2)
N = jm.Placeholder("N")
x = jm.BinaryVar("x", shape=(N, N))
i = jm.Element("i", belong_to=(0,N))
j = jm.Element("j", belong_to=(0,N))
t = jm.Element("t", belong_to=(0,N))
ここで、jm.Placeholder
は定数を表現しており、ここでは距離行列と都市数を表現するのに用いています。
巡回セールスマン問題においては、この距離行列と都市数によってさまざまな問題が表現されることになります。
バイナリ変数を表現するのが、jm.BinaryVar
です。
ここでは、のバイナリ変数を定義しています。
次に、総和などで添字の範囲を表現するためにjm.Element
を用いて、i
,j
,t
という添字を定義しています。
JijModelingでは、jm.Problem
インスタンスを作成し、それに目的関数や制約条件を追加していきます。
では次に、定義した変数を用いて目的関数を定義していきます。
problem = jm.Problem("TSP")
problem += jm.sum([t, i, j], dist[i, j] * x[t, i] * x[(t + 1) % N, j])
problem
総和はjm.sum
を用いて表現することができます。
jm.sum
の最初の引数は総和をとる添字で、TSPの目的関数では、3つの添字について総和を取るので、それらの添字(jm.element
)をリストで渡しています。
次に、制約条件を追加していきます。
# 制約条件1 : 位置に関するonehot制約
problem += jm.Constraint(
"onehot_location",
x[:, i].sum() == 1,
forall=i,
)
# 制約条件2 : 時間に関するonehot制約
problem += jm.Constraint(
"onehot_time",
x[t, :].sum() == 1,
forall=t,
)
jm.Constraint
を用いて制約条件表現することができます。
最初の引数は制約条件の名前、2つ目の引数が制約条件を表現した数式になります。
この制約条件には、オプション引数として、forall
というものがあります。
これは、数理モデルにおいて、「任意のについて」や、「」というように表現されているものをJijModelingで表現するための引数です。
また、制約条件の中に現れているx[:,i].sum()
という記法は、jm.sum(t,x[t,i])
という書き方の糖衣構文になっています。
最後に、今記述した数理モデルを確認してみましょう。
problem
説明で用いた数理モデルと同じ数式が表示されていることがわかります。
以上で、数理モデルの構築は終わりです。 このようにJijModelingを用いると、手元の数式と見比べながら数理モデルをコーディングしていくことができます。
問題データの作成#
問題の数理モデルができたので、次に問題に使うデータを作成します。 ここでは、単純な都市数10で、都市間の距離をランダムにした問題を解いていきます。
import matplotlib.pyplot as plt
import numpy as np
inst_N = 5
np.random.seed(3)
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)
数理モデルを作る際に用いたjm.Placeholder
にデータを代入するので、Placeholderの名前をkeyに持つ辞書でデータを渡す必要があります。
今回の問題では、N
とdist
に値を渡す必要があります。
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, "dist": inst_d}
instance_data
{'N': 5,
'dist': array([[0. , 0.7866063 , 0.73643374, 0.84577089, 0.56967619],
[0.7866063 , 0. , 0.4251585 , 0.21078131, 0.36540009],
[0.73643374, 0.4251585 , 0. , 0.26950348, 0.64576184],
[0.84577089, 0.21078131, 0.26950348, 0. , 0.54552992],
[0.56967619, 0.36540009, 0.64576184, 0.54552992, 0. ]])}
JijModeling-Transpilerを用いた数理モデルのQUBOへの変換#
数理モデルとデータの用意ができたので、次にJijModeling-Transpilerを用いて数理モデルとデータをQUBOに変換します。
まずはjijmodeling_transpiler
をインポートしましょう。
import jijmodeling_transpiler as jmt
QUBOへの変換は、以下の手順で行うことができます。
# compile
compiled_model = jmt.core.compile_model(problem, instance_data, {})
# get qubo model
pubo_builder = jmt.core.pubo.transpile_to_pubo(compiled_model=compiled_model, relax_method=jmt.core.pubo.RelaxationMethod.AugmentedLagrangian)
qubo, const = pubo_builder.get_qubo_dict(multipliers={})
jmt.core.compile_model
を用いて、JijModelingで作成したproblem
と用意したinstance_data
をコンパイルします。
コンパイルされたモデルをjmt.core.pubo.transpile_to_pubo
に渡すことで、polynomial unconstrained binary optimization (PUBO: 多項式の制約なしバイナリ最適化問題)へトランスパイルします。
トランスパイルされた結果から.get_qubo_dict
メソッドを用いることで、QUBOの情報を辞書型で得ることができます。
OpenJijを用いた最適化の実行#
ここまでで得たQUBOを用いることで、これまでと同様に最適化計算を行うことができます。
import openjij as oj
sampler = oj.SASampler()
res = sampler.sample_qubo(Q=qubo, num_reads=1)
JijModeling-Transpilerの機能を用いれば、最適化によって得られた結果をより解析しやすい形にすることができます。
decode_from_openjij
機能を使ってみましょう。
sampleset = jmt.core.pubo.decode_from_openjij(res, pubo_builder, compiled_model)
sampleset
SampleSet(record=Record(solution={'x': [(([1, 0, 2, 3, 4], [1, 4, 0, 2, 3]), [1.0, 1.0, 1.0, 1.0, 1.0], (5, 5))]}, num_occurrences=[1]), evaluation=Evaluation(energy=[-6.8035391666144855], objective=[2.703473530206421], constraint_violations={"onehot_location": [0.0], "onehot_time": [0.0]}, constraint_forall={"onehot_location": [[0], [1], [2], [3], [4]], "onehot_time": [[0], [1], [2], [3], [4]]}, constraint_values=[{"onehot_location": [0.0, 0.0, 0.0, 0.0, 0.0], "onehot_time": [0.0, 0.0, 0.0, 0.0, 0.0]}], penalty={}), measuring_time=MeasuringTime(solve=SolvingTime(preprocess=None, solve=None, postprocess=None), system=SystemTime(post_problem_and_instance_data=None, request_queue=None, fetch_problem_and_instance_data=None, fetch_result=None, deserialize_solution=None), total=None), metadata={})
得られた解はsampleset.record
のsolution
の中に入っています。
この中に結果は疎行列の形で入っています。
大事なのが最初の二つの要素で、ひとつ目が行列の中のindexそして、二つ目がそのindexにおける値が入っています。
バイナリ変数の場合には、1となった値のみが入っているので、通常、値には1しか入っていません。
sparse_index,value,_ = sampleset.record.solution['x'][0]
sparse_index
([1, 0, 2, 3, 4], [1, 4, 0, 2, 3])
sampleset.evaluation
の中には、最適化によって得られた評価値として、エネルギーや目的関数、そして、制約条件の破れが入っています。
ここで、制約条件を満たした解が得られたかをチェックしてみます。
sampleset.evaluation.constraint_violations
{'onehot_location': array([0.]), 'onehot_time': array([0.])}
制約条件を満たした解が得られているようです。
そこで、この解を可視化してみましょう。 ルートをプロットする場合には、時間に関するindexを用いて、都市のindexをソートしてあげることで、どの順番に都市を回るかという順序を得ることができます。
time_indices, city_indices = zip(*sorted(zip(*sparse_index)))
time_indices, city_indices
((0, 1, 2, 3, 4), (4, 1, 0, 2, 3))
plt.plot(x_pos, y_pos, 'o',markersize=12)
plt.xlim(0, 1)
plt.ylim(0, 1)
for i, city_index in enumerate(city_indices[:-1]):
next_city_index = city_indices[i+1]
plt.plot([x_pos[city_index],x_pos[next_city_index ]],[y_pos[city_index],y_pos[next_city_index ]],c = "blue")
plt.plot([x_pos[city_indices[-1]],x_pos[city_indices[0]]],[y_pos[city_indices[-1]],y_pos[city_indices[0]]],c = "blue")
[<matplotlib.lines.Line2D at 0x7fb25c3ce610>]