14
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

wildcatで巡回セールスマン問題を試してみた

Posted at

#はじめに
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

とすることで、目的の形式に持っていきます。
##ステップ1
いきなりは難しいので、まず下記の形式に持っていきます。

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
急ぎ足でやったので、ミスってたらごめんなさい。
#参考

  1. wildcatで4x4の数独を解いてみた
  2. 物理のいらない量子アニーリング入門
14
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
14
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?