- Aさん、Bさんの二人が、それぞれ重さ$w$が異なる$N$個の荷物を出来るだけ均等に分配する。但し、それぞれの荷物には持ちにくさ$D$があり、持ちにくさを考慮した上で、出来るだけ荷物を均等に分配する問題を、Dwaveのアニーリングマシンを用いて解く。
参考：東北大学大関先生講義 第2回
Keyword ：Dwave基本設定、Pyquboの使い方、罰則法
1. 基本準備/設定
1.1 準備1:量子アニーリング用モジュールのインストール
まず、初めに、以下のモジュールをインストールする。
- dwave-ocean-sdk: Dwaveアニーリングマシンのモジュール
- openjij: アニーリングシュミレーター
- pyqubo: Qubo行列作成サポートモジュール
pip install dwave-ocean-sdk
pip install openjij
pip install pyqubo
1.2 準備2: Dwaveの初期設定
続いて、Dwaveのtoken、enpoint(アクセスURL)、使用するマシンの指定を行う。実際には、関数として数値を指定するだけ。後で、Dwaveにアクセスする際に、上記の関数を指定することになる。
token = "DEV-73ac79954ef2135*************************"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage_system4.1"
更に、Dwaveマシンを指定するDwaveSamplerと、Qbitへの割り当て行うEmbeddingCompositeをインポートする
from dwave.system import DWaveSampler, EmbeddingComposite
続いて、インポートしたDWaveSamplerを用いてDwaveマシンを指定する。
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)
尚、使用するマシンの指定方法を含むDwaveの使用方法は、Qiitaの記事を参考：https://qiita.com/y_daito/items/aad74517927e393564ee
1.3 準備3: 各種モジュールのインポート
openjijからSQASampler、pyquboからArray、その他としてnumpy、pandas、matplotlib.pyplotをインポートする。
from openjij import SQASampler
from pyqubo import Array
import pandas as pd
2. 定式化
2.1 エネルギー式の整理
今回の問題は、重さ$w$がそれぞれ異なる$N$個の荷物をAさんとBさんの二人に可能な限り均等に分担する。均等であるとは、AさんとBさんの総重量がイコールになれば良いということ。アニーリングのエネルギー式は二次関数である必要がある。
そこで、Aさんの重量を$W_A$、Bさんの重量を$W_B$とすると、エネルギー式$E(x)$を以下のように出来る。
E(x) =(W_A - W_B)^2
続いて、アニーリング問題は、場合分けが0 or 1のスピンで規定されないといけない。
そこで、$i$番目の荷物をAさんがとる場合、$x_i=1$とし、Aさんがとらない場合、つまりBさんがとる場合を$x_i=0$とする。
続いて、それぞれの荷物の重さを$w_i$とする。
すると、Aさんの荷物の重量は以下となる。
W_A = \sum_{i=1}^N w_i x_i
一方で、全ての荷物の重量$W$は以下となる。
W = \sum_{i=1}^N w_i
ここで、$W_B$を$∑$で表現することも可能だが、$x_i$がAが荷物を持つかどうかの判断となっている為、新たに$x_i$を設定する必要があり面倒。そこで、$W$、$W_A$だけで、$E(x)$を表現してみる。
E(x) = (W_A - W_B)^2 \\
E(x) = (2W_A - W)^2 \\
ここで、$W_A$のみ$∑$で表記してみる。尚、$W$は定数なので$∑$を使用せず、$W$のまま記載する。
E(x) = (2\sum_{i=1}^N w_i x_i - W)^2 \\
続いて、これを展開する。$∑$の添え字はべき乗する場合、同じ添え字を使用せず、別の添え字に変更する。
E(x) = (2\sum_{i=1}^N w_i x_i - W)(2\sum_{j=1}^N w_j x_j - W) \\
E(x) = 4\sum_{i=1}^N \sum_{j=1}^N w_i w_j x_i x_j - 4 W\sum_{i=1}^N w_i x_i + W^2
尚、上記第二項は本来添え字が$i$と$j$で別れるが、いずれも独立している為、$i$でも$j$でも同じ、Aさんの総重量を示す。従い、$i$でまとめている。
2.2 コスト関数の整理
更に問題を複雑にする為に、夫々の荷物の組み合わせ事に、持ちにくさ$D_{ij}$を定義する。 ここでの添え字$i、j$は、例えば1番目と2番目などと分けることが目的。尚、$i$番目と$i$番目の荷物を同時に持つことは無いので、$D_{ii}$という場合は存在しない。
ここでコスト関数$C_1(x)$を整理すると以下となる。
C_1(x) = \sum_{i=1}^N \sum_{j=1, i≠j}^N D_{ij} x_i x_j
2.3 ハミルトニアンの作成
上記で整理したエネルギー式とコスト関数を足し合わせることでハミルトニアンを作成する。尚、コスト関数は、結果見ながら影響度合いを検討せざるを得ない。影響度を$α_1$と定義した場合、ハミルトニアンは以下となる。
H(x) = E(x) + α_1 C_1(x) \\
H(x) = 4\sum_{i=1}^N \sum_{j=1}^N w_i w_j x_i x_j - 4 W\sum_{i=1}^N w_i x_i + W^2 + α_1\sum_{i=1}^N \sum_{j=1, i≠j}^N D_{ij} x_i x_j
数式の整理は以上。
以降、数式の実装を進める。
3. 実装
3.1 問題の設定
今回分担する荷物$N$が10個とした場合。夫々の重さ$w_i$、夫々の持ちにくさ$D_{ij}$、またコスト関数の影響度$α_1$を以下の通り設定する。
N = 10
w = [1,2,3,4,5,6,7,8,9,10]
W = sum(w)
D = []
for i in range(1,N+1):
_a = [i*x for x in w]
_a[i-1] = 0
D.append(_a)
a1 = 0
3.2 Pyquboを使用せず直接的にQubo行列を表現する場合
まず始めにQubo行列を格納する箱を作成する。$N$行x$N$行の正方行列である。
Qubo = np.zeros(N**2).reshape(N,N)
Qubo
<Output>
array([[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., 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., 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., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]
Qubo行列の箱に、以下の数式の各項の結果を入力していく。
H(x) = 4\sum_{i=1}^N \sum_{j=1}^N w_i w_j x_i x_j - 4 W\sum_{i=1}^N w_i x_i + W^2 + α_1\sum_{i=1}^N \sum_{j=1, i≠j}^N D_{ij} x_i x_j
まず、第1項から、表現する。$x_i$、$x_j$以外をそのまま代入すれば良い。
for i in range(N):
for j in range(N):
Qubo[i][j] = 4 * w[i] * w[j]
Qubo
<Output>
array([[ 4., 8., 12., 16., 20., 24., 28., 32., 36., 40.],
[ 8., 16., 24., 32., 40., 48., 56., 64., 72., 80.],
[ 12., 24., 36., 48., 60., 72., 84., 96., 108., 120.],
[ 16., 32., 48., 64., 80., 96., 112., 128., 144., 160.],
[ 20., 40., 60., 80., 100., 120., 140., 160., 180., 200.],
[ 24., 48., 72., 96., 120., 144., 168., 192., 216., 240.],
[ 28., 56., 84., 112., 140., 168., 196., 224., 252., 280.],
[ 32., 64., 96., 128., 160., 192., 224., 256., 288., 320.],
[ 36., 72., 108., 144., 180., 216., 252., 288., 324., 360.],
[ 40., 80., 120., 160., 200., 240., 280., 320., 360., 400.]])
続いて、第2項。
for i in range(N):
Qubo[i][i] = Qubo[i][i] - 4*W*w[i]
Qubo
<Output>
array([[ -216., 8., 12., 16., 20., 24., 28., 32.,
36., 40.],
[ 8., -424., 24., 32., 40., 48., 56., 64.,
72., 80.],
[ 12., 24., -624., 48., 60., 72., 84., 96.,
108., 120.],
[ 16., 32., 48., -816., 80., 96., 112., 128.,
144., 160.],
[ 20., 40., 60., 80., -1000., 120., 140., 160.,
180., 200.],
[ 24., 48., 72., 96., 120., -1176., 168., 192.,
216., 240.],
[ 28., 56., 84., 112., 140., 168., -1344., 224.,
252., 280.],
[ 32., 64., 96., 128., 160., 192., 224., -1504.,
288., 320.],
[ 36., 72., 108., 144., 180., 216., 252., 288.,
-1656., 360.],
[ 40., 80., 120., 160., 200., 240., 280., 320.,
360., -1800.]])
次に、第3項であるが、ここには$x_i$も$x_j$も出てこないので、Qubo行列には影響を与えない。従い、作業無。
最後に、第4項のコスト関数。
for j in range(N):
if not i == j:
Qubo[i][j] = Qubo[i][j] + a1*D[i][j]
Qubo
array([[ -216., 8., 12., 16., 20., 24., 28., 32.,
36., 40.],
[ 8., -424., 24., 32., 40., 48., 56., 64.,
72., 80.],
[ 12., 24., -624., 48., 60., 72., 84., 96.,
108., 120.],
[ 16., 32., 48., -816., 80., 96., 112., 128.,
144., 160.],
[ 20., 40., 60., 80., -1000., 120., 140., 160.,
180., 200.],
[ 24., 48., 72., 96., 120., -1176., 168., 192.,
216., 240.],
[ 28., 56., 84., 112., 140., 168., -1344., 224.,
252., 280.],
[ 32., 64., 96., 128., 160., 192., 224., -1504.,
288., 320.],
[ 36., 72., 108., 144., 180., 216., 252., 288.,
-1656., 360.],
[ 40., 80., 120., 160., 200., 240., 280., 320.,
360., -1800.]])
以上でQubo行列の作成は完了。
最後に、DwaveマシンにQubo行列を投げてQubo行列を最小化する$x_i$の組み合わせを抽出する。
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(Qubo, num_reads=10)
print(solution.record)
<Output>
[([0, 0, 0, 1, 0, 1, 0, 1, 0, 1], -3024., 1, 0.)
([1, 0, 1, 1, 0, 0, 0, 0, 1, 1], -3024., 2, 0.)
([1, 0, 0, 0, 0, 0, 1, 0, 1, 1], -3024., 1, 0.)
([1, 0, 1, 0, 0, 1, 0, 1, 0, 1], -3024., 1, 0.)
([1, 0, 0, 1, 0, 1, 1, 0, 0, 1], -3024., 1, 0.)
([0, 1, 0, 1, 0, 1, 0, 1, 1, 0], -3016., 1, 0.)
([0, 0, 0, 1, 1, 0, 1, 0, 0, 1], -3016., 1, 0.)
([0, 1, 0, 0, 1, 0, 0, 1, 1, 0], -2976., 1, 0.)
([1, 0, 0, 1, 0, 0, 0, 0, 1, 1], -2976., 1, 0.)]
1番目、2番目、3番目の場合でのAさんの重量、Bさんの重量を表現すると以下の通り。
for i in range(5):
wa = np.dot(solution.record[i][0],w)
wb = W-wa
print(f"{i+1}番目の結果：Aさんの総重量は{wa}、Bさんの総重量は{wb}")
1番目の結果：Aさんの総重量は28、Bさんの総重量は27
2番目の結果：Aさんの総重量は27、Bさんの総重量は28
3番目の結果：Aさんの総重量は27、Bさんの総重量は28
4番目の結果：Aさんの総重量は28、Bさんの総重量は27
5番目の結果：Aさんの総重量は28、Bさんの総重量は27
3.3 Pyquboを使用してQubo行列を表現する場合
Pyquboを使用する場合、以下の数式の$x_i$の集合であるベクトル$\vec{x}$を定義することから始まる。
H(x) = (2\sum_{i=1}^N w_i x_i - W)^2 + α_1\sum_{i=1}^N \sum_{j=1, i≠j}^N D_{ij} x_i x_j \\
H(x) = 4\sum_{i=1}^N \sum_{j=1}^N w_i w_j x_i x_j - 4 W\sum_{i=1}^N w_i x_i + W^2 + α_1\sum_{i=1}^N \sum_{j=1, i≠j}^N D_{ij} x_i x_j
from pyqubo import Array
x = Array.create(name="x", shape=(N),vartype="BINARY")
<Output>
Array([Binary('x[0]'), Binary('x[1]'), Binary('x[2]'), Binary('x[3]'), Binary('x[4]'), Binary('x[5]'), Binary('x[6]'), Binary('x[7]'), Binary('x[8]'), Binary('x[9]')])
PyquboのArrayでは、0 or 1のバイナリ変数、-1 or +1のイジング変数のいずれかが選択できる。vartypeで"BINARY"or"SPIN"を選択する。
Pyquboを使用する場合、直接Qubo行列を作成するのでなく、上記の$H(x)$を先ほど定義した$\vec{x}$を用いて、数式自体を表現する。出来上がった数式からQubo行列を抜き出すという手順になる。Pyquboの最大の利点は展開せずとも計算出来る点であるが、展開しても使用可能。まずは、展開した場合で計算を行う。
以下のように一挙に表現する。
### まずハミトルニアンを準備する。
H=0
### 第1項の表現
for i in range(N):
for j in range(N):
H = H + 4*w[i]*w[j]*x[i]*x[j]
### 第2項の表現
for i in range(N):
H = H - 4*W*w[i]*x[i]
### 第3項の表現
H = H + W**2
### 第4項の表現
for i in range(N):
for j in range(N):
if not i == j:
H = H + a1*D[i][j]*x[i]*x[j]
$H(x)$の作成が完了すれば、以下の方法でQubo行列のみを抽出する。
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(Qubo, num_reads=10)
print(solution.record)
for i in range(5):
wa = np.dot(solution.record[i][0],w)
wb = W-wa
print(f"{i+1}番目の結果：Aさんの総重量は{wa}、Bさんの総重量は{wb}")
[([0, 0, 1, 0, 0, 1, 0, 1, 0, 1], -3024., 1, 0. )
([1, 0, 1, 1, 1, 1, 0, 0, 0, 1], -3016., 1, 0. )
([1, 0, 0, 0, 1, 1, 1, 0, 0, 1], -3016., 1, 0. )
([1, 1, 1, 0, 0, 1, 1, 0, 0, 1], -3016., 1, 0. )
([1, 0, 0, 0, 0, 0, 1, 1, 1, 0], -3000., 1, 0. )
([0, 0, 1, 0, 1, 0, 1, 0, 0, 1], -3000., 1, 0. )
([1, 1, 0, 0, 1, 1, 1, 1, 0, 0], -3016., 1, 0.1)
([0, 1, 0, 1, 0, 1, 0, 0, 0, 1], -2904., 1, 0. )
([1, 1, 0, 0, 0, 1, 1, 0, 1, 1], -2800., 1, 0.1)
([0, 1, 1, 1, 0, 0, 0, 0, 0, 1], -2736., 1, 0. )]
1番目の結果：Aさんの総重量は27、Bさんの総重量は28
2番目の結果：Aさんの総重量は29、Bさんの総重量は26
3番目の結果：Aさんの総重量は29、Bさんの総重量は26
4番目の結果：Aさんの総重量は29、Bさんの総重量は26
5番目の結果：Aさんの総重量は25、Bさんの総重量は30
続いて、数式を展開しない場合で$H(x)$を表現してみる。
### まずハミトルニアンを準備する。
H = 0
### 第1項の表現
H = (2*np.dot(w,x)-W)**2
### 第2項の表現
for i in range(N):
for j in range(N):
if not i == j:
H = H + a1*D[i][j]*x[i]*x[j]
続いて、Qubo行列を抽出。
model = H.compile()
Qubo, offset = model.to_qubo()
print(Qubo)
<Output>
{('x[6]', 'x[0]'): 56.0, ('x[6]', 'x[6]'): -1344.0, ('x[7]', 'x[6]'): 448.0, ('x[5]', 'x[0]'): 48.0, ('x[5]', 'x[1]'): 96.0, ('x[8]', 'x[3]'): 288.0, ('x[3]', 'x[3]'): -816.0, ('x[7]', 'x[2]'): 192.0, ('x[9]', 'x[3]'): 320.0, ('x[9]', 'x[1]'): 160.0, ('x[9]', 'x[9]'): -1800.0, ('x[9]', 'x[2]'): 240.0, ('x[5]', 'x[4]'): 240.0, ('x[4]', 'x[1]'): 80.0, ('x[1]', 'x[1]'): -424.0, ('x[2]', 'x[2]'): -624.0, ('x[5]', 'x[3]'): 192.0, ('x[7]', 'x[5]'): 384.0, ('x[4]', 'x[3]'): 160.0, ('x[6]', 'x[1]'): 112.0, ('x[4]', 'x[4]'): -1000.0, ('x[8]', 'x[5]'): 432.0, ('x[1]', 'x[0]'): 16.0, ('x[8]', 'x[2]'): 216.0, ('x[7]', 'x[3]'): 256.0, ('x[8]', 'x[8]'): -1656.0, ('x[2]', 'x[1]'): 48.0, ('x[7]', 'x[0]'): 64.0, ('x[3]', 'x[2]'): 96.0, ('x[6]', 'x[2]'): 168.0, ('x[8]', 'x[7]'): 576.0, ('x[7]', 'x[4]'): 320.0, ('x[6]', 'x[3]'): 224.0, ('x[4]', 'x[2]'): 120.0, ('x[7]', 'x[7]'): -1504.0, ('x[0]', 'x[0]'): -216.0, ('x[9]', 'x[6]'): 560.0, ('x[8]', 'x[6]'): 504.0, ('x[6]', 'x[5]'): 336.0, ('x[2]', 'x[0]'): 24.0, ('x[8]', 'x[1]'): 144.0, ('x[4]', 'x[0]'): 40.0, ('x[9]', 'x[8]'): 720.0, ('x[9]', 'x[0]'): 80.0, ('x[9]', 'x[4]'): 400.0, ('x[5]', 'x[2]'): 144.0, ('x[6]', 'x[4]'): 280.0, ('x[9]', 'x[5]'): 480.0, ('x[8]', 'x[0]'): 72.0, ('x[3]', 'x[0]'): 32.0, ('x[3]', 'x[1]'): 64.0, ('x[7]', 'x[1]'): 128.0, ('x[5]', 'x[5]'): -1176.0, ('x[9]', 'x[7]'): 640.0, ('x[8]', 'x[4]'): 360.0}
同じく、Dwaveで解いてみると。
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(Qubo, num_reads=10)
print(solution.record)
for i in range(5):
wa = np.dot(solution.record[i][0],w)
wb = W-wa
print(f"{i+1}番目の結果：Aさんの総重量は{wa}、Bさんの総重量は{wb}")
[([1, 0, 1, 1, 1, 1, 0, 1, 0, 0], -3024., 1, 0. )
([0, 1, 0, 1, 0, 1, 1, 1, 0, 0], -3024., 1, 0. )
([0, 1, 0, 1, 1, 0, 0, 1, 1, 0], -3024., 1, 0. )
([1, 0, 0, 1, 0, 1, 1, 1, 0, 0], -3016., 1, 0. )
([1, 1, 0, 1, 0, 0, 0, 0, 1, 1], -3016., 1, 0. )
([1, 0, 0, 1, 1, 0, 1, 1, 0, 0], -3000., 1, 0. )
([1, 1, 0, 0, 1, 0, 1, 0, 0, 1], -3000., 1, 0. )
([1, 1, 0, 1, 1, 1, 1, 0, 0, 0], -3000., 1, 0. )
([0, 0, 0, 0, 0, 0, 0, 1, 1, 0], -2584., 1, 0. )
([0, 0, 0, 1, 1, 1, 1, 0, 1, 0], -2976., 1, 0.1)]
1番目の結果：Aさんの総重量は27、Bさんの総重量は28
2番目の結果：Aさんの総重量は27、Bさんの総重量は28
3番目の結果：Aさんの総重量は28、Bさんの総重量は27
4番目の結果：Aさんの総重量は26、Bさんの総重量は29
5番目の結果：Aさんの総重量は26、Bさんの総重量は29
これまでと同様の回答が得られた。PyQuboを使用することで、展開せずにハミルトニアンを表現できる為、数学的にプログラミングが容易になる。
3.4 Openjijシュミレーターを使用する場合
参考まで、アニーリングシュミレーターを使用した場合のプログラミング流れを整理する。
### シュミレーターの呼び出し
sampler_jij = SQASampler()
### 計算の実施
solution = sampler_jij.sample_qubo(Qubo,num_reads=10)
print(solution.record)
### 結果の表示
for i in range(5):
wa = np.dot(solution.record[i][0],w)
wb = W-wa
print(f"{i+1}番目の結果：Aさんの総重量は{wa}、Bさんの総重量は{wb}")
[([0, 0, 0, 0, 1, 1, 1, 0, 0, 1], -3024., 1)
([1, 0, 1, 0, 0, 1, 0, 1, 1, 0], -3024., 1)
([0, 1, 1, 1, 1, 0, 1, 1, 0, 0], -3016., 1)
([1, 1, 1, 1, 0, 0, 0, 1, 1, 0], -3024., 1)
([0, 0, 1, 0, 0, 0, 1, 1, 1, 0], -3024., 1)
([1, 0, 0, 0, 0, 1, 1, 0, 1, 0], -2944., 1)
([1, 1, 1, 1, 0, 0, 0, 1, 0, 1], -3024., 1)
([0, 0, 1, 0, 1, 0, 0, 0, 1, 1], -3024., 1)
([0, 0, 0, 0, 1, 1, 1, 0, 1, 0], -3024., 1)
([0, 1, 0, 0, 0, 0, 1, 1, 0, 1], -3024., 1)]
1番目の結果：Aさんの総重量は28、Bさんの総重量は27
2番目の結果：Aさんの総重量は27、Bさんの総重量は28
3番目の結果：Aさんの総重量は29、Bさんの総重量は26
4番目の結果：Aさんの総重量は27、Bさんの総重量は28
5番目の結果：Aさんの総重量は27、Bさんの総重量は28
計算精度は、アニーリング実機と比較して低下するが、作成したQubo行列が機能しているか堂かの検証を無償で実施出来る。
おわりに
今回はアニーリングによる最適化の一番の基礎を整理した。アニーリングは0 or 1しか選択できない為、$x_i$の0はAさんが持たない、1はAさんが持つ、と表現した。しかし、これではAさん、Bさんの二人の場合しか表現できない。それでは、Aさん、Bさん、Cさんと増えてきた場合、つまりハミルトニアンが2次式以上になる場合、どのように$x_i$を設定すべきか、考えないといけない。
今後の課題: ハミルトニアンが2次式以上になる場合の対応