@yuichirominato 2018.06.11更新

D-WaveでVW社の交通最適化アプリケーションの実装を解く

D-Wave QUBO イジング 交通流最適化 量子アニーリング

はじめに

組合せ最適化問題を社会実装する際に量子コンピュータを活用した実装方法を確認します。使用するマシンはカナダのD-Wave社のマシンで、自社で借りているものを使用しました。

参考資料など

今回はVW社が北京で行なった交通最適化問題の論文と、その参考資料をベースとします。

Traffic flow optimization using a quantum annealer
Florian Neukart, Gabriele Compostella, Christian Seidel, David von Dollen, Sheir Yarkoni, Bob Parney
(Submitted on 4 Aug 2017 (v1), last revised 9 Aug 2017 (this version, v2))
https://arxiv.org/pdf/1708.01625.pdf
参考資料
https://www.dwavesys.com/sites/default/files/VW.pdf

どんな問題か?

北京の市内から空港までの交通混雑状況をD-Waveを使用した組合せ最適化問題で混在解消するという社会実験です。下記の図の左側の混雑状況が右側のように緩和されます。

引用:https://www.dwavesys.com/media-coverage/automotive-it-vw-cio-technology-being-readied-address-real-issues

手順

手順は下記の通りです。論文中のchaper2から引用しました。既存計算機を「古典計算機」と表現しています。

1. Classical: Pre-process map and GPS data.
(古典計算機)地図とGPSデータからデータの準備をする。

2. Classical: Identify areas where traffic congestion occurs.
(古典計算機)次に混雑が起こっている場所を特定する。

3. Classical: Determine spatially and temporally valid alternative routes for each car in the dataset, if possible.
(古典計算機)データセット内の自動車に対して代替ルートを提案する

4. Classical: Formulate the minimization problem as a QUBO (to minimize congestion in road segments on overlapping routes).
(古典計算機)混雑が緩和するような組合せ最適化問題に落とし込む。その際のQUBOと呼ばれる形式を採用する。

5. Hybrid Quantum/Classical: Find a solution that reduces congestion among route assignments in the whole traffic graph.
(古典計算機・量子コンピュータハイブリッド)問題を古典計算機による分割と量子コンピュータによる最適化を繰り返す。

6. Classical: Redistribute the cars based on the results.
(古典計算機)上記の得られた答えから自動車の位置を再配置する。

7. Iterate over steps 2 to 6 until no traffic congestion is identified.
上記のステップを混雑が緩和されるまで繰り返し計算する。

引用:https://arxiv.org/pdf/1708.01625.pdf

具体的な解き方の参考

上記にはQUBOなどの専門用語が出てきます。QUBOは問題をバイナリの{0,1}で考える手法です。下記の過去の記事などを参照してください。

イジングモデルや量子アニーリングについては過去の記事
「量子アニーリング、イジングモデルとフレームワーク」
http://blog.mdrft.com/post/6

一般的なイジングモデルのNP問題の解法はこちら
「NP問題のイジング(量子アニーリングなど向け)」
http://blog.mdrft.com/post/42

では、実際の解法の準備

論文や参考資料など色々ありますが、全体を簡単に説明しているのは下記の参考資料です。これを全体に拡張することで複雑な混雑緩和をD-Waveで行うことができます。今回は一部を抜き出して、簡単な例題で行なって見たいと思います。

「Quantum Computing at Volkswagen:
Traffic Flow Optimization using the D-Wave Quantum Annealer」
引用:https://www.dwavesys.com/sites/default/files/VW.pdf

今回は解きたい問題をバイナリで考えます。バイナリ値で考える場合、こちらはQUBOと呼ばれる規格に合わせて問題を実装します。

##(例題)自動車2台の経路最適化を行う
下記のような道に対して2台の自動車の最適化を行います。例題は上記のVW社のものをそのまま使います。

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

上記はAからBまでの経路を12本の道路で表現しています。
まずは、自動車1と自動車2に対して、AからBまでの取りうる経路を3つずつ決めます。
さらに、その3つずつの経路の中から1本経路を選びます。

上記の例では、自動車1は#1から#3の経路提案があり、#1を選んでいます。
自動車2は同様で、#2の道路を選んでいます。
ちなみに、自動車1と自動車2は全く同じ3つの経路を経路提案として持っていますが、全く違うものを使っても大丈夫です。

ここで準備するのは、自動車1に対して提案された経路をそれぞれ、$Q_{11},Q_{12},Q_{13}$とし、自動車2に対して提案は$Q_{21},Q_{22},Q_{23}$とします。

また、それぞれの経路提案の中に具体的にどのセグメント(道路)が含まれているかを書き出しておきます。

混雑状況を計算する

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

次に混雑を計算します。ここでは、先ほどの経路提案でセグメントが何回含まれているかをカウントし、それを二乗することでコストを算出しています。

たとえば、セグメントのs0は提案$Q_{11},Q_{12},Q_{21},Q_{22}$に含まれます。コスト関数は上記のように$(Q_{11}+Q_{12}+Q_{21}+Q_{22})^2$のようになります。実際には、$Q_{11}とQ_{22}$が選択されているので、s0の経路のコストは、$(1+0+0+1)^2=4$となります。経路が選ばれていれば1、選ばれていなければ0です。

すべてのセグメントの混雑のコストを足し合わせると、この場合には12になります。

ルートをずらして見て、コストを再確認する

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

ルートを実際にずらして選択して見た場合に、コストが下がるかどうかを確認します。
上記の図から見ると自動車2がルートの$Q_{23}$を選択した場合、全体のコストは8に低下しています。このように、コストの下がるようなルートの組み合わせをD-Waveを使って計算します。

制約条件

さて、これでコストを選んでいけば解けそうですが、イジングマシンの場合にはそう簡単には行きません。
ここで制約条件と呼ばれる、本来満たすべき量子ビット間の関係性を表す式が必要になります。イジングで上記の問題を解こうとするとき、制約をかけない場合には、1台の自動車が複数の経路を選択してしまうことが起こり得てしまいます。この制約条件をイジングマシンでは明示的に数式で表現をし、コスト関数に導入する必要があります。それが下記です。

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

この場合、自動車1と自動車2がそれぞれルートを1つだけ選択するという条件を付与する必要があります。その条件の作り方は、$Q_{11}もしくはQ_{12}もしくはQ_{13}$のうち、どれか1つが一となればいいという条件を最小値問題で表すとき、$(Q_{11}+Q_{12}+Q_{13}-1)^2$とします。上記条件を満たすとき、この式は0となり、満たさないときは最小値になりません。イジングマシンは最小値を取ろうとするので、これをうまく導入することで制約条件を満たすことができます。

この際、元のコスト関数との影響率を調整するハイパーパラメータKを導入します。

条件を満たすQUBOmatrix

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

参考資料によると上記の条件を満たすQUBOのmatrixは上記のようになるようです。あとは自分でKを調整して、色々解いて見ることになります。

本当に解けるのか細かく見ていく

まずは混雑状況を計算する組合せ最適化問題の式を作ります。

引用:https://www.dwavesys.com/sites/default/files/VW.pdf

こちらを引用して、下記の方にExampleとして式が載っています、これらを式展開してQUBOmatrixというものを作ります。便宜上$Q_{11}=q_0,Q_{12}=q_1,Q_{13}=q_2,Q_{21}=q_3,Q_{22}=q_4,Q_{23}=q_5$とすると、全体の道の混雑を表すコスト関数(ハミルトニアン)は、上記を書き換えて、

$
\begin{eqnarray*}
H&=&(q_0+q_1+q_3+q_4)^2+(q_0+q_1+q_3+q_4)^2+(q_0+q_3)^2+(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2\\
&=&2(q_0+q_1+q_3+q_4)^2+2(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+3(q_2+q_5)^2\\
&=&4q_0^2+4q_0q_1+8q_0q_3+4q_0q_4+4q_1^2+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2^2 +2q_2q_4+8q_2q_5+4q_3^2+4q_3q_4+4q_4^2+2q_4q_5+4q_5^2
\end{eqnarray*}
$

ここで、量子ビットはバイナリ値なので、$x_i^2 = x_i$が成り立つことから上記式は下記のように展開できる。

$
\begin{eqnarray*}
H&=&4q_0+4q_0q_1+8q_0q_3+4q_0q_4+4q_1+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2 +2q_2q_4+8q_2q_5+4q_3+4q_3q_4+4q_4+2q_4q_5+4q_5
\end{eqnarray*}
$

これを係数を見やすくQUBOmatrixの形に直すと、

$
\begin{array}
– & q_0 & q_1 & q_2 & q_3 & q_4 & q_5\\
q_0 & 4 & 4 & 0 & 8 & 4 & 0\\
q_1 & & 4 & 2 & 4 & 8 & 2\\
q_2 & & & 4 & 0 & 2 & 8\\
q_3 & & & & 4 & 4 & 0\\
q_4 & & & & & 4 & 2\\
q_5 & & & & & & 4
\end{array}
$

これはまだ制約条件を満たしていないので、制約条件をここに足し合わせていく。

制約条件の展開

制約条件は前述から、

$
K(q_0+q_1+q_2-1)^2+K(q_3+q_4+q_5-1)^2
$

これを展開して、

$
Kq_0^2+2Kq_0q_1+2Kq_0q_2-2Kq_0+Kq_1^2+2Kq_1q_2-2Kq_1+Kq_2^2-2Kq_2+Kq_3^2+2Kq_3q_4+2Kq_3q_5-2Kq_3+Kq_4^2+2Kq_4q_5-2Kq_4+Kq_5^2-2Kq_5+2K
$

さらに$q_i^2 = q_i$より、

$
-Kq_0+2Kq_0q_1+2Kq_0q_2-Kq_1+2Kq_1q_2-Kq_2-Kq_3+2Kq_3q_4+2Kq_3q_5-Kq_4+2Kq_4q_5-Kq_5+2K
$

こちらの係数を先ほどのQUBOmatrixに導入して、

$
\begin{array}
– & q_0 & q_1 & q_2 & q_3 & q_4 & q_5\\
q_0 & 4-K & 4+2K & 2K & 8 & 4 & 0\\
q_1 & & 4-K & 2+2K & 4 & 8 & 2\\
q_2 & & & 4-K & 0 & 2 & 8\\
q_3 & & & & 4-K & 4+2K & 2K\\
q_4 & & & & & 4-K & 2+2K\\
q_5 & & & & & & 4-K
\end{array}
$

となりました。

QUBOからイジングへ変換

QUBOではバイナリで考えましたが、実際のD-Waveは{-1,1}のイジングスピンで考えますので変換をします。

下記の変換ルールです。

$Jij=\frac{1}{4}Q_{i,j}\\h_i=\frac{1}{2}Q_{i,i}+\frac{1}{4}\sum_{i<j}Q_{i,j}$

こちらを適用すると、さきほどのQUBOmatrixはイジングへと変換され下記のようになります。

$
\begin{array}
– & q_0 & q_1 & q_2 & q_3 & q_4 & q_5\\
q_0 & \frac{12+K}{2} & \frac{2+K}{2} & \frac{K}{2} & 2 & 1 & 0\\
q_1 & & \frac{14+K}{2} & \frac{1+K}{2} & 1 & 2 & \frac{1}{2}\\
q_2 & & & \frac{10+K}{2} & 0 & \frac{1}{2} & 2\\
q_3 & & & & \frac{12+K}{2} & \frac{2+K}{2} & \frac{K}{2}\\
q_4 & & & & & \frac{14+K}{2} & \frac{1+K}{2}\\
q_5 & & & & & & \frac{10+K}{2}
\end{array}
$

こちらをD-Waveにいれていきます。

Kを調整して答えを出して見る

上記の6量子ビットの完全結合は17量子ビットで作れました。

上記の太線は量子ビットのコピーです。
この回路に上記を代入します。

たとえば、K=1として、全体を1/10してみて、D-Waveに入力して見ると、

実際にエネルギーが最小な状態で、$q_2$と$q_3$が選ばれました。

実際に、$q_2$に対応するセグメントはs2,s7,s10,s11、$q_3$に対応するセグメントはs0,s3,s6,s9と重なる道路がなく、混雑状況として最小になっているのが確認できました。

Kの値の調整はやはりコツが少し必要でした。
以上で、実際にD-Waveで交通混雑の緩和式を最小値問題で解くことができました。

あわせて読みたい