はじめに

こんにちは，株式会社Ridge-iエンジニアの中村です．本記事では量子アニーリング手法を用いた最適化についてご紹介します．

量子アニーリングでは，イジングモデルという統計力学の理論モデルを用いて最適化計算を行うため，通常の数理モデルの定式化からQUBOと呼ばれる式へ変換する必要があります．私自身QUBOについてはあまり知らなかったため，これを機に調べてみました．またスケジュールタスクに量子アニーリングを適用した論文 "Applying Quantum Annealing for Shift Scheduling Problem for Call Centers" の問題設定を模擬して実装してみました．本記事が最適化・量子アニーリングに興味がある方にとって情報収集の一助となれば幸いです．

最適化とは

最適化とは，与えられた制約条件の下で，目的関数の値を最小(もしくは最大)にする解を求める問題のことをいいます[1]．最適化の問題はいくつか種類に分類することができますが，今回は組み合わせ最適化(離散最適化)について考えます．組み合わせ最適化の有名な例ですと，巡回セールスマン問題やナップサック問題があります．

巡回セールスマン問題は「1人のセールスマンがある都市を出発した後，目的地となる複数の都市を経由し帰ってくる場合に，どのルートが最短となるかを求める問題」です．ただし以下の制約を満たす必要があります．

ある都市を出発し，すべての都市を経由して出発の都市に戻る
1度訪れた都市を経由して，別の都市に向かうのは禁止

上記の条件の中で，最短で回ることができるルートを求めることになります．

組み合わせ最適化問題では，要素の数(上記の巡回セールスマン問題でいえば訪問都市数)に対して解の選択肢が指数関数的に増大する「組み合わせ爆発」が起こることが知られています．仮に都市の数が20の場合，組み合わせの数は20!=2432902008176640000通りとなり，100MIPSで演算をするコンピューターで約39年かかることになります[2]．

これから紹介する量子アニーリング(Quantum Annealing)[3]は，このような組み合わせ最適化問題を現実時間で解くことが期待されています．

イジングモデルとQUBO式

イジングモデル

量子アニーリングは，イジングモデルと呼ばれる統計力学の理論モデルを用いて組み合わせ最適化問題を考えるため，はじめにイジングモデルについて説明します．イジングモデルは強磁性体の振る舞いを調べるために考案された模型の1つです．内容としては， $d$ 次元立法格子を考えて，それぞれの原子が±1の2値を取るスピンを持っている場合，スピン配置でのエネルギー $H$ は以下の式で表すことができるというものです．

$\displaystyle H = \sum\_{i>j} J\_{ij} σ_i σ_j + \sum h_i σ_i \tag{1}$

ここで $σ_i \in \{1, -1\}$ は $i$ 番目のスピンを表すバイナリ変数です． $J_{ij}$ は $i$ と $j$ のスピン間相互作用を表しており， $h_i$ は各スピンが{1, -1}のどちらかを取りたがる性質の強さを表しています[4]．

具体的に $d=1$ 次元の場合を考えると下記の図のようになります．ここでは簡単のため $J_{ij} = -J$ ， $h_{i} = 0$ とします．

上向きを $+1$ ，下向きを $-1$ とした場合，エネルギーは以下のように計算されます．

$\displaystyle H = J + J + (-J) + (-J) + (-J) = -J$

計算式からもわかるように $J_{ij}=-J$ の場合，隣り合うスピンが同じ方向を向いている方がエネルギーが低くなります．

上記で示した $H$ は最適化問題におけるコスト関数に対応しており，イジングモデルのエネルギー $H$ を最小化するようなスピンの配列を求める問題に置き換えて解くのが，イジングモデルを用いた最適化計算になります．

QUBO(Quadratic Unconstrained Binary Optimization)式について

イジングモデルを用いた組み合わせ最適化問題を解くために，以下のようにスピン $σ_i \in \{1, -1\}$ を $x_i \in \{0, 1 \}$ に変換します．

$\displaystyle x_i = \frac{\sigma_i + 1}{2} \tag{2}$

$(2)$ 式を用いて， $(1)$ 式のイジングモデルで表されるエネルギー $H$ を書き換えると以下のようになります．

$\displaystyle H = \sum\_{i < j}(Q\_{ij}x_ix_j) + \sum\_{i}(Q\_{ii}x_i) \tag{3}$

$(3)$ 式はQUBO(Quadratic Unconstrained Binary Optimization)式と呼ばれています．イジングモデルを用いた組み合わせ最適化問題ではQUBO式が用いられます．

イジングモデルを用いた最適化ロジック

上記は組み合わせ最適化問題をイジングモデルとして定式化する内容でしたが，実際に最小化されたコスト関数を得る方法について紹介します．

焼きなまし法(Simulated Annealing)

焼きなまし(Annealing)とは，温度を少しずつ下げて欠損のない均一な物質を作る手法のことです．これを最適化問題に当てはめて考案されたメタヒューリスティックアルゴリズムが焼きなまし法(Simulated Annealing)[5]になります．温度 $T$ において，あるスピンの配列から1つのスピンを変えた場合を考えます．エネルギーが下がる場合は新しいスピンの配列を取り入れ，上がる場合は以下の確率 $P$ で新しいスピン配列を取り入れます．

$\displaystyle P = exp(-ΔH/T)$

ここで $ΔH(= H_{new} - H_{old})$ はスピン配列を変えた後から変える前のエネルギーを引いたエネルギー差になります．エネルギーが上がった場合でも確率的に新しいスピン配列を受け入れることで，グローバルな解を見つけることができます(この式から，温度 $T$ が高く，エネルギー差が小さいほど新しいスピン配列を受け入れやすくなる)．このような確率的な振る舞いは物理的な熱揺らぎを模擬しています．

量子アニーリングとSimulated Quantum Annealing

焼きなまし法は，熱揺らぎを模擬してグローバルな解を見つけるアルゴリズムでした．ミクロな世界での物理現象を説明する量子力学では，熱揺らぎとは別に量子揺らぎという確率的な効果が存在します．熱揺らぎの代わりに，量子揺らぎを用いて最適化計算を行うのが量子アニーリングです．量子アニーリングのエネルギーは以下のように記述されます[6]．

$\displaystyle H(t) = H\_{ising} + \Gamma (t)H’$

$t$ は時刻で， $\Gamma (t)$ は初期時刻では非常に大きく， $t$ が大きいときは影響が小さくなるように設定されます． $H'$ は取りうる全てのスピン配置を表しており，時間と共に $H_{ising}$ のスピン状態を求めることができます．

量子アニーリングでは純粋な量子効果を用いるため，従来のコンピューターで計算を行うことは難しいです．そこで，Simulated Quantum Annealingと呼ばれる手法を用います．Simulated Quantum Annealingでは，量子効果を模擬した $H(t)$ に対して，焼きなまし法を行います．

実践 : シフトスケジューリング問題

ここではシフトスケジューリング問題をQUBO式として定式化し，量子アニーリングとSimulated Quantum Annealingで解を求めます．シフトスケジューリング問題はコールセンターのアルバイトなどのシフトの組み合わせを最適化する組み合わせ最適化問題のことです．

今回は，量子アニーリングがスケジュール問題に対して有用であることを示した論文 Applying Quantum Annealing for Shift Scheduling Problem for Call Centers[7] を参考に問題設定を作成しました．また論文内で紹介されている定式化手法やペナルティ係数(後述)の決定方法を実践しました．ソースコードはGithubのRepositoryをご参考ください．

量子アニーリングはD-waveマシンを使用し，Simulated Quantum Annealingはアニーリングシミュレータとしてソフトウェア実装をした OpenJij [8] を用いました．ただしD-waveマシンでの計算は利用制限(アカウント作成から1ヶ月は1分のみ使用可能という制限)があるため，主な計算はOpenJijを用いて行い，お試しでD-waveマシンを使用しました．最後に，実際にD-waveマシンで計算を行うためのセットアップ方法もご紹介します(具体的な計算手順は省略)．

定式化

論文の問題設定は以下です．

シフトスケジュールを作成する
1日はいくつかのタームに分かれている
タームごとに何人必要か決められているので，なるべく必要人数に合わせて作業員を出勤させたい
作業員ごとに希望出勤回数を出しており，なるべく希望にあった回数だけ働かせる
希望出勤回数とは別に，休み希望リストがあり，休み希望の日にち・タームは出勤させてはならない
作業員の中にはグループが存在しており，同じグループの作業員は一緒に出勤しなければならない

以上の設定から定数，変数を定義し，制約条件と目的関数を定式化します．

定数

{作業員1, …, 作業員 }
- 今回は作業員が6人いると設定
{1日, …, 日}
- 今回は7日間あると設定
{ターム1, ..., ターム }
- 今回は3タームあると設定
: 日のタームに必要なブース数
- 今回は1日2ブース(2人)必要と設定
: 作業員が日のタームへ割り当て可能かを表すバイナリ変数
- 実際の設定はGithubのConfigファイルに記載していますが，出力されるスケジュールに割り当て不可日を赤枠にしています
: 作業員aの希望割り当てシフト数
- 実際の設定はGithubのConfigファイルに記載しています
: グループに作業員が所属しているかいないかのバイナリ
- 今回は作業員 1,2 はグループ 1，作業員 3,4 はグループ 2，作業員 5,6 はグループ 3 と設定

変数

$x_{adt} \in \{0, 1\}$ : 作業員 $a$ が $d$ 日目のターム $t$ に割り振られたか

制約条件

作業員が希望していないターム( )にはシフトを割り振らない
- $r_{adt} = 0$ の場合

$\displaystyle x\_{adt} = 0 \tag{4}$

同じグループに割り振られた作業員全員が同じシフトに割り振られる

$\displaystyle \sum_{a} x_{adt} \times G_{ga} = (\sum_{a}G_{ga} \rm{or} \;0) \tag{5}$

2つ目の制約式について補足します．下図( [7] の論文より抜粋) のようにグループ $g = 1$ に作業員 $a=1,2$ が，グループ $g = 2$ に作業員 $a=3,4$ が，グループ $g = 3$ に作業員 $a=5,6$ が属している場合を考えます．ある日 $d$ のあるターム $t$ で，作業員 $1,2,6$ が出勤する( $x_{1dt},x_{2dt},x_{6dt}=1$ )とした場合，グループ $g=3$ は $\sum G_{3a} \ne \sum x_{adt} \times \sum G_{3a}$ となり制約を満たしません．つまり2つ目の制約式は，同じグループの作業員は同じシフトに割り当てられるという条件を反映していることになります．

目的関数

希望ブース数と割り当てられる作業員数が近い

$\displaystyle \sum_{dt}(\sum_{a}x_{adt} - S_{dt}) \tag{6}$

実際に割り当てる作業員数と希望シフト数の差が近い

$\displaystyle \sum_{a}(\sum_{dt}x_{adt} - R_{a}) \tag{7}$

QUBOの定式化

量子アニーリングを適用するためにQUBO式を構築する必要があります．式 $(4)$ ~ $(7)$ を使って，以下のようにQUBO式を定義します．

$\displaystyle H = \lambda_{i} H_{i} + \lambda_{j} H_{j} + W_{k} H_{k} + W_{l} H_{l} \tag{8}$

$H_{i}, H_{j}, H_{k}, H_{l}$ は式 $(6),$ $(7)$ の目的関数と，式 $(4),$ $(5)$ の制約条件を表しています．

$\displaystyle H_{i} = \sum_{dt}(\sum_{a}x_{adt} - S_{dt})^2 \tag{9}$

$\displaystyle H\_{j} = \sum_{a}(\sum_{dt}x_{adt} - R_{a})^2 \tag{10}$

$\displaystyle H_{k} = (1 - r_{adt})x_{adt} \tag{11}$

$\displaystyle H_{l} = \sum_{dtg} \{ (\sum_{a} G_{ga} - \sum_{a} x_{adt} G_{ga}) \sum_{a} x_{adt} G_{ga} \} \tag{12}$

$\lambda_{i}$ , $\lambda_{j}$ は目的関数 $H_{i}$ , $H_{j}$ の係数で， $H_{i}$ の値が小さい場合は希望ブース数と割り当てられる作業員数が近いことを表し， $H_{j}$ が小さい場合は実際に割り当てる作業員数と希望シフト数の差が近いことを表します．つまり係数に差をつけることで，求解における目的関数の優先度を調整することができます．今回は簡単のため $\lambda_{i} = \lambda_{j} = 1$ としました．

同様に $W_k$ , $W_l$ は制約項 $H_{k}$ , $H_{l}$ のペナルティ係数です．制約を満たしていれば制約項は $0$ となり，制約違反する場合は正の値を取るように定義されています．ペナルティ係数の大きさは，目的関数の値を十分に超えるような値に設定設定することで制約違反が起こらないような解を求めることができます．

ペナルティ係数のチューニング

QUBO式に落とし込んだ後は制約項のペナルティ係数を決めます．この節では $H$ を energy， $W_{k} = w_{desire},\; W_{l} = w_{group}$ と呼ぶことにします．論文では2つの係数をグリットサーチを用いて計算を行います．同じ係数の組み合わせで複数回計算を行って，制約を満たす解(実行可能解)が何回出現したかをカウントし，ヒートマップを用いて可視化します．上記の手法で実験を行うと以下のようなヒートマップが得られました．今回は同じ係数の組み合わせで100回計算行って確率を算出しています．

この図から $w_{desire}$ , $w_{group}$ で実行可能解の確率が均衡の取れている値を選択します．今回は $w'_{desire} = 1.5$ , $w'_{group}=5.0$ としました．

上記で決めた係数は比率とし，次に基準値(base)となる値を求めます．つまり $w_{desire} = w'_{desire} \times \rm{base}$ , $w_{group} = w'_{group} \times \rm{base}$ として計算を行います．

baseを $0\sim 2$ の範囲で変化させ，再度実行可能解の確率を算出すると以下のようになります．右図は100個の解における energy の最小・平均・最大です．

各base値に対して100回実行した際の実行可能解が得られる確率と，エネルギーの最小・平均・最大

実行可能解の確率が 1 付近で安定し，かつ energy が安定している base 値を選択します．今回はbase $=1.4$ としました．よって係数は以下のようになります．

$\displaystyle w\_{desire} = w'\_{desire} \times \rm{base} = 1.5 \times 1.4 \tag{13}$

$\displaystyle w\_{group} = w'\_{group} \times \rm{base} = 5.0 \times 1.4 \tag{14}$

結果

式 $(13), (14)$ のペナルティ係数を用いて最適化計算を行うと以下のようなスケジュールが得られました．縦軸は各作業員，横軸は時間軸である日にちとタームを表しています(d1_t2 の場合，1日目の2タームを表す)．

今回設定した作業員の割り当てられないターム $r_{adt}$ は赤枠で，制約1つ目は満たせていることが確認できました．また作業員 $\{1,2\}, \{3,4\}, \{5,6\}$ が同じグループと設定しているため，制約2つ目についても満たせていることが確認できました．

これまではOpenJijを使って Simulated Quantum Anneaing で計算を行いましたが，D-waveマシンを使用して計算も行いました．しかし上記で決めたペナルティ係数では制約を満たさない解が出てしまったので，手動でいくつかのパターンを試してみた結果，以下のような制約を満たす実行可能解が得られました．

上記のスケジュールは，目的関数の1つ目である "希望ブース数と割り当てられる作業員数が近い(各タームで2人割り当てられるのが理想)" を大きく逸脱する結果となっています．D-waveで計算を行う際は，別途ペナルティ係数の調整を行うとより質の高い解(制約を満たし，目的関数が小さい=Energyが小さい)が得られるかもしれません．

量子アニーリングは問題設定をQUBO式に落として係数を決定する必要がありますが，他の計算手法に比べて比較的早く解を求めることが可能です．論文では，他のメタヒューリスティックアルゴリズムや，数理最適化ソルバーを使った場合との求解速度の比較を行っています．D-waveを用いた量子アニーリング手法は他と比べて速度が早く，変数が増えるほどその差は顕著になっています．

今回使用した量子アニーリングと Simulated Quantum Annealing で計算速度を比べると，本タスクでは量子アニーリングの方が100倍計算速度が早かったです．

タスク	量子アニーリング	Simulated Quantum Annealing
スケジュール算出	約2 秒	約0.019 秒 (19 m秒)

補足 : D-Waveマシンでの実行準備

Simulated Quantum Annealingの場合はローカルPCで実行可能ですが，D-waveマシンで計算を行う場合はアカウントを作成し，計算を投げる必要があります．

ここではD-waveマシンで計算を行うためのセットアップ方法について紹介します． - アカウント作成方法とAPI TOKENの取得 - D-Wave Leap(https://cloud.dwavesys.com/leap/)にアクセスし，Sign upする - Sign up後，メールからログインする - 自身のDashboard画面から，API TOKENを取得する - ライブラリ用意 - D-Waveマシンを使用するソフトウェアがまとまっている開発キットをインストールします

pip install dwave-ocean-sdk

この先の計算方法については長くなるためここでは省略します．詳細はGithubをご覧ください．

まとめ

本記事では，イジングモデルから始まり量子アニーリングを用いて組み合わせ最適化問題を解く方法をご紹介しました．量子アニーリングは，組み合わせ最適化問題をより高速に解く手法としてさらなる普及が期待されているので，ご興味のある方は調べてみてください．

さいごに

Ridge-iでは様々なポジションで積極採用中です．カジュアル面談も可能ですのでご興味がある方は是非ご連絡ください．

Ridge-i.com

Ridge-institute R&D Blog

量子アニーリングを用いた最適化問題の定式化方法について