はじめに
wildcatの使い方は https://github.com/mdrft/wildcat を見てください。
こちらも参考になるかもしれません。wildcatで4x4の数独を解いてみた
また、巡回セールスマン問題の詳細はこちらの記事が参考になります。物理のいらない量子アニーリング入門
それでは標題の通り試してみます。
ハミルトニアンの作成
参考[2]から、巡回セールスマン問題のハミルトニアンは下記の様になります。
H = \sum_{t,a,b}D_{a,b}x_{t,a}x_{t+1,b} + A\sum_{t}\left(\sum_{a}x_{t,a}-1\right)^2 + A\sum_{a}\left(\sum_{t}x_{t,a}-1\right)^2
また巡回する町の数を $n$ とします。
QUBOに変形
ハミルトニアンを下記の形式に変形します。
H = \sum_{i,j}Q_{i,j}x_{i}x_{j}
ここで $H$ を見ると $x$ の一次の項が出そうですが、$x$ が0か1しか値が取らないという性質を利用して、
x_{t,a} = x_{t,a}^2
とすることで、目的の形式に持っていきます。
ステップ１
いきなりは難しいので、まず下記の形式に持っていきます。
H = \sum_{t_1,t_2,a,b}J_{t_1,t_2,a,b}x_{t_1,a}x_{t_2,b}
第一項
\begin{eqnarray}
\sum_{t,a,b}D_{a,b}x_{t,a}x_{t+1,b} = \sum_{t_1,t_2,a,b}\delta_{t_1+1,t_2}D_{a,b}x_{t_1,a}x_{t_2,b}
\end{eqnarray}
ここで $\delta_{a,b}$ はクロネッカーデルタです。
第二項
\begin{eqnarray}
A\sum_{t}\left(\sum_{a}x_{t,a}-1\right)^2 = A\sum_{t_1,t_2,a,b}\left(\delta_{t_1,t_2}-2\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + An
\end{eqnarray}
ここで $x_{t,a} = x_{t,a}^2$ を使いました。
第三項
\begin{eqnarray}
A\sum_{a}\left(\sum_{t}x_{t,a}-1\right)^2 = A\sum_{t_1,t_2,a,b}\left(\delta_{a,b}-2\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + An
\end{eqnarray}
ここで $x_{t,a} = x_{t,a}^2$ を使いました。
したがって、$H$ は下記の様になりました。
H=\sum_{t_1,t_2,a,b}\left(\delta_{t_1+1,t_2}D_{a,b}+A\delta_{t_1,t_2}+A\delta_{a,b}-4A\delta_{a,b}\delta_{t_1,t_2}\right)x_{t_1,a}x_{t_2,b} + 2An
ここで、定数項はエネルギーの基準をずらしてゼロにして考えます。あとは $J$ を計算して、それを変形して $Q$ を計算します。これで数学の計算は終わりにします。
コード
それではコードで書いていきます。
まず、ライブラリを読み込みます。そして断熱計算の設定を行います。今回はローカルで計算をします。また、計算が成功する確率を上げるため、初期の温度を1000度にしました。
from wildcat.solver.qubo_solver import QuboSolver
from wildcat.network.local_endpoint import LocalEndpoint
from wildcat.annealer.simulated.simulated_annealer import SimulatedAnnealer
from wildcat.annealer.simulated.single_spin_flip_strategy import SingleSpinFlipStrategy
from wildcat.annealer.simulated.temperature_schedule import TemperatureSchedule
from wildcat.util.matrix import hamiltonian_energy
schedule = TemperatureSchedule(initial_temperature=1000, last_temperature=0.1, scale=0.8)
strategy = SingleSpinFlipStrategy(repetition=10)
annealer = SimulatedAnnealer(schedule=schedule, strategy=strategy)
local_endpoint = LocalEndpoint(annealer=annealer)
次に参考[2]を参考に町の位置と距離の計算を定義し、また定数は適当に400としました。
#座標
positions = np.array((
(24050.0000, 123783),
(24216.6667, 123933),
(24233.3333, 123950),
(24233.3333, 124016),
(24250.0000, 123866),
(24300.0000, 123683),
(24316.6667, 123900),
(24316.6667, 124083),
(24333.3333, 123733),
))
#距離の計算
def dist(a, b):
a = np.array(a)
b = np.array(b)
return np.sqrt(((a - b)**2).sum())
#町の数
N = len(positions)
#定数A
A = 400
それでは一気に $J$ を計算していきます。
#行列Dを計算
D = np.empty((N, N), dtype = np.float64)
for i in range(N):
for j in range(N):
D[i, j] = dist(positions[i], positions[j])
#クロネッカーデルタを定義
def delta(i, j):
if(i == j):
return 1
else:
return 0
#行列Jを計算
J = np.empty((N, N, N, N), dtype = np.float64)
for a in range(N):
for t1 in range(N):
for b in range(N):
for t2 in range(N):
J[a, t1, b, t2] = delta(t1+1, t2) * D[a, b] + A * delta(t1, t2) + A * delta(a, b)\
- 4 * A * delta(a, b) * delta(t1, t2)
#定数項の計算
C = 2 * A * N
それでは式を変形して $Q$ を$J$から計算していきます。
#Qの計算
Q = np.empty((N*N,N*N), dtype = np.float64)
x = 0
for a in range(N):
for t1 in range(N):
y = 0
for b in range(N):
for t2 in range(N):
Q[x,y] = J[a,t1,b,t2]
y = y + 1
x = x + 1
$Q$ が出来たので、実際に最適化の計算を行います。
#最適化計算
solver = QuboSolver(Q)
def callback(arrangement):
e = solver.hamiltonian_energy(arrangement)
print("Energy: ", solver.hamiltonian_energy(arrangement))
print("Spins: ", arrangement)
arrangement = solver.solve(callback, endpoint=local_endpoint).result()
[Out]
Energy: -179479.61895467353
Spins: [-1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 -1 1
-1 -1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1
-1 -1 -1 1 -1 -1 -1 -1 -1]
ここでEnergyは低いほうが正解近くなります。またSpinsは量子ビットのスピンの向きを表しています。
結果をわかりやすくすると下記のようになります。
result = np.empty((N, N), dtype = np.int)
i = 0
for x in range(N):
for y in range(N):
if(arrangement[i] == 1):
result[x, y] = 1
else:
result[x, y] = 0
i += 1
print("Result: \n", result)
[Out]
Result:
[[0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 1 0 0 0]
[0 0 1 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0]
[0 0 0 1 0 0 0 0 0]]
見てみると、巡回の条件を満たしていない部分があることがわかります。
もう一度やってみると、今度は成功しました。
Result:
[[1 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 1 0 0 0]
[0 0 0 1 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 0 1 0 0]]
成功したときの距離を計算すると下記のようになります。
#距離を計算する
length = 0
for t in range(N):
for a in range(N):
for b in range(N):
if(t < N-1):
if(result[a, t] == 1 and result[b, t+1] == 1):
length += D[b, a]
else:
for c in range(N):
if(result[a, t] == 1 and result[c, 0] == 1):
length += D[a, c]
print("Length:", length)
[Out]
Length: 5203.634602384536
興味がある方はぜひ試してみてください。コードは上げておきます。https://github.com/Yokomakura/wildcat-TSP/blob/master/wildcat%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3.ipynb
急ぎ足でやったので、ミスってたらごめんなさい。