量子コンピュータデータサイエンティスト講座：経路最適化

はじめに

量子コンピュータとモビリティを語る際には外せない交通最適化、経路最適化の問題を解いてみたいとおもいます。

参考

今回は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

手順

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

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

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

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

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

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

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

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

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

例題

１、１台の車に2つのルート選択を許容する（古典）
２、現在の選択ルートから混雑状況を計算（古典）
３、混雑度を最小にするようにルート選択を最小化（量子）

スタートAからゴールBまで12の道路が通し番号0から11まで付けられています。自動車１と自動車２があり、それぞれ指定されたルートを取ります。ルートはそれぞれ2つずつ提案されており、そのうちの一つを取ります。今回は一番混雑度が低くなるようなルートを選択します。今回は最小値問題を解くアルゴリズムを利用し、混雑状況とルート選択の制約を同時に満たすような計算を学んでいきます。

car1
route1-1(q0)：s0,s3,s6,s9
route1–2 (q1)：s0,s3,s8,s11

car2
route2–1(q2):s0,s3,s8,s11
route2–2(q3):s2,s7,s10,s11

混雑度を計算

まず、道路の混雑度を上記の提案ルートを利用して求めます。それぞれの道の混雑を量子ビットを使って表現します。道s0から順番にコストを足し合わせていきます。コストはその道の登場する量子ビットを足し合わせて２乗します。途中 $0^2=0$ または、 $1^2=1$ より次数を下げています。

h=(q_0+q_1+q_2)^2+(q_0+q_1+q_2)^2+q_0^2+q_0^2+(q_1+q_2)^2+(q_1+q_2+q_3)^2+q_3^2+q_3^2+q_3^2\\ =4q_0+4q_1+4q_2+4q_3+4q_0q_1+4q_0q_2+8q_1q_2+2q_1q_3+2q_2q_3

これがそのまま定式化（ハミルトニアン）になります。

制約条件を課す（ソフト制約）

しかしここでは、ルート選択は１つだけを選ぶわけですから、 $q_0$ か $q_1$ の片方は1で、もう片方は0となるのが条件としてついてきます。また、 $q_2$ と $q_3$ も同様です。今回は、まず、この制約条件をソフト制約といって、式の形で実装してみます。
ここで、ソフト制約を課すには、

(q_0+q_1-1)^2 + (q_2+q_3-1)^2 = q_0^2 + q_1^2 + 2q_0q_1 - 2q_0 - 2q_1 + 1 + q_2^2 + q_3^2 + 2q_2q_3 - 2q_2 - 2q_3 + 1 = 2q_0q_1 - q_0 - q_1 + 2q_2q_3 - q_2 - q_3 + 2

という式を準備します。途中 $0^2=0$ または、 $1^2=1$ より次数を下げています。この式は $q_0$ か $q_1$ のどちらかが1の時かつ、 $q_2$ もしくは $q_3$ のどちらかが1の時にこの式は最小値をとります。今回利用するのは最小値を求める計算をベースにして計算を行いますので、最小値を利用して制限をかけます。

混雑度と制約の式を組み合わせる

そして、上記の混雑度の計算の式と制約条件の式を組み合わせます。二つの式はそのままでは組み合わせることができず、今回はソフト制約の式のほうにラグランジュパラメータと呼ばれる調整変数をつけます。
$$
h=4q_0+4q_1+4q_2+4q_3+4q_0q_1+4q_0q_2+8q_1q_2+2q_1q_3+2q_2q_3 + M*(2q_0q_1 - q_0 - q_1 + 2q_2q_3 - q_2 - q_3)
$$
上記、ソフト制約の定数は最小値問題では除いても影響がないので、省いています。
$$
h = (4-M)q_0 + (4-M)q_1 + (4-M)q_2 + (4-M)q_3 + (4+2M)q_0q_1 + 4q_0q_2 + 8q_1q_2 + 2q_1q_3 + (2+2M)q_2q_3
$$
式を整理するとこのようになりました。これで定式化は完了です。このMをパラメータとして調整しながら、制約を満たしながら、式を最小にする組合せを探します。

プログラミング

blueqatを利用して、上記の式をプログラムに落としていきます。

import numpy as np
import networkx as nx
from blueqat.wq import *
import matplotlib.pyplot as plt

%matplotlib inline

まずはコスト関数をのほうを実装してみます。

#initialize
A = np.zeros((4,4))

#weight on distance
A[0][0] = A[1][1] = A[2][2] =A[3][3] = 4

A[0][1] = 4
A[0][2] = 4
A[1][2] = 8
A[1][3] = 2
A[2][3] = 2

A

array([[4., 4., 4., 0.],
       [0., 4., 8., 2.],
       [0., 0., 4., 2.],
       [0., 0., 0., 4.]])

#initialize
B = np.zeros((4,4))

#weight on distance
B[0][0] = B[1][1] = B[2][2] = B[3][3] = -1

B[0][1] = 2
B[0][2] = 2

B

array([[-1.,  2.,  2.,  0.],
       [ 0., -1.,  0.,  0.],
       [ 0.,  0., -1.,  0.],
       [ 0.,  0.,  0., -1.]])

この二つを接続します。

#set a QUBO
M = 1
qubo = A+B*M

#show as network
G = nx.from_numpy_matrix(qubo)
nx.draw_networkx(G)
plt.show()

<Figure size 432x288 with 1 Axes>

#solved on blueqat
M = 5
qubo = A+B*M

a = Opt()
a.qubo = qubo
res = a.run()

res

[1, 0, 0, 1]

このように答えが出ました。答えは、 $q_0=1$ かつ $q_3=1$ なので、この二つのルートが選択されました。それによってコストが最小になります。

ハード制約

上記はソフト制約を実行しました。量子もつれを利用してハード制約を実行してみます。ハード制約ではラグランジュパラメータの調整が不要になります。

import blueqat
from blueqat import vqe,Circuit
from blueqat.pauli import X, Y, Z, I

#initialize
A = np.zeros((4,4))

#weight on distance
A[0][0] = A[1][1] = A[2][2] =A[3][3] = 4

A[0][1] = 4
A[0][2] = 4
A[1][2] = 8
A[1][3] = 2
A[2][3] = 2

A

array([[4., 4., 4., 0.],
       [0., 4., 8., 2.],
       [0., 0., 4., 2.],
       [0., 0., 0., 4.]])

まずはコスト関数を準備しました、ソフト制約の場合と同じです。そして、それをゲートマシンで計算できるようにします。今回また、シミュレーションでは量子断熱ステップというステップ数を決める必要がありますが、step=2としてみます。これは精度に影響します。

q = pauli(A)
step = 2

次に、今回のハード制約は量子もつれと呼ばれるものでつくることができます。mixerと初期状態の両方を自分で準備する必要があります。mixerは今回は詳しい説明は省略しますが、二つの量子ビットを入れ替えるXYmixerを利用します。これで、0と1の量子ビットの入れ替え、2と3の量子ビットの入れ替えを実現します。

また、初期状態はXYmixerに対応した形で、量子ビットが01と10でもつれるように作られます。

#mixer and init state
mixer = 0.5*X[0]*X[1] + 0.5*Y[0]*Y[1] + 0.5*X[2]*X[3] + 0.5*Y[2]*Y[3]
init = Circuit().h[0].cx[0,1].x[0].h[2].cx[2,3].x[2]

そしてこれを実行します。

result = vqe.Vqe(vqe.QaoaAnsatz(q, step, init, mixer)).run()
resv = result.most_common(1)[0][0]
print(resv)

(1, 0, 0, 1)

答えが出ました。1001という答えがきちんとでました。このようにハード制約を利用して経路最適化問題を解くことができました。

まとめ

ソフト制約とハード制約をみました。使い方を確認して身につけましょう。以上です。