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次式以上になる場合の対応