焼きなまし法と巡回セールスマン問題

以前、[グラフ理論の基礎](https://qiita.com/maskot1977/items/e1819b7a1053eb9f7d61)や[グラフ理論の基礎をmatplotlibアニメーションで](https://qiita.com/maskot1977/items/b55c6d1c85bcf0ec0d08)などの記事を書きましたが、グラフ理論の中でNP困難とされる難解な問題の一つである「巡回セールスマン問題」の近似解法として「焼きなまし法」(Simulated Annealing)を使った方法を試してみました。

# 巡回セールスマン問題

「巡回セールスマン問題(traveling salesman problem、TSP)は、頂点の集合と各頂点間の移動コストが与えられたとき、全ての頂点をちょうど一度ずつ巡り出発地に戻る巡回路のうちで総移動コストが最小のものを求める組合せ最適化問題です。

# 焼きなまし法

焼きなまし法（Simulated Annealing、SA）は、金属工学における「焼きなまし」から名前を取った大域的最適化問題を解く方法。局所解に落ちてもすぐに脱出できるような「温度」が高い状態で探索を開始し、徐々に「温度」を下げていって大域的最適解を探します。

# 都道府県の県庁所在地データ

[グラフ理論の基礎](https://qiita.com/maskot1977/items/e1819b7a1053eb9f7d61)や[グラフ理論の基礎をmatplotlibアニメーションで](https://qiita.com/maskot1977/items/b55c6d1c85bcf0ec0d08)のデータを例に、巡回セールスマン問題を解いてみようと思います。


```python
import urllib.request
url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/location.txt'
urllib.request.urlretrieve(url, 'location.txt') # データのダウンロード
```




    ('location.txt', <http.client.HTTPMessage at 0x7f9f4e7685c0>)




```python
import pandas as pd
japan = pd.read_csv('location.txt')
japan
```


<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Town</th>
      <th>Longitude</th>
      <th>Latitude</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>Sapporo</td>
      <td>43.06417</td>
      <td>141.34694</td>
    </tr>
    <tr>
      <th>1</th>
      <td>Aomori</td>
      <td>40.82444</td>
      <td>140.74000</td>
    </tr>
    <tr>
      <th>2</th>
      <td>Morioka</td>
      <td>39.70361</td>
      <td>141.15250</td>
    </tr>
    <tr>
      <th>3</th>
      <td>Sendai</td>
      <td>38.26889</td>
      <td>140.87194</td>
    </tr>
    <tr>
      <th>4</th>
      <td>Akita</td>
      <td>39.71861</td>
      <td>140.10250</td>
    </tr>
    <tr>
      <th>5</th>
      <td>Yamagata</td>
      <td>38.24056</td>
      <td>140.36333</td>
    </tr>
    <tr>
      <th>6</th>
      <td>Fukushima</td>
      <td>37.75000</td>
      <td>140.46778</td>
    </tr>
    <tr>
      <th>7</th>
      <td>Mito</td>
      <td>36.34139</td>
      <td>140.44667</td>
    </tr>
    <tr>
      <th>8</th>
      <td>Utsunomiya</td>
      <td>36.56583</td>
      <td>139.88361</td>
    </tr>
    <tr>
      <th>9</th>
      <td>Maebashi</td>
      <td>36.39111</td>
      <td>139.06083</td>
    </tr>
    <tr>
      <th>10</th>
      <td>Saitama</td>
      <td>35.85694</td>
      <td>139.64889</td>
    </tr>
    <tr>
      <th>11</th>
      <td>Chiba</td>
      <td>35.60472</td>
      <td>140.12333</td>
    </tr>
    <tr>
      <th>12</th>
      <td>Tokyo</td>
      <td>35.68944</td>
      <td>139.69167</td>
    </tr>
    <tr>
      <th>13</th>
      <td>Yokohama</td>
      <td>35.44778</td>
      <td>139.64250</td>
    </tr>
    <tr>
      <th>14</th>
      <td>Niigata</td>
      <td>37.90222</td>
      <td>139.02361</td>
    </tr>
    <tr>
      <th>15</th>
      <td>Toyama</td>
      <td>36.69528</td>
      <td>137.21139</td>
    </tr>
    <tr>
      <th>16</th>
      <td>Kanazawa</td>
      <td>36.59444</td>
      <td>136.62556</td>
    </tr>
    <tr>
      <th>17</th>
      <td>Fukui</td>
      <td>36.06528</td>
      <td>136.22194</td>
    </tr>
    <tr>
      <th>18</th>
      <td>Kofu</td>
      <td>35.66389</td>
      <td>138.56833</td>
    </tr>
    <tr>
      <th>19</th>
      <td>Nagano</td>
      <td>36.65139</td>
      <td>138.18111</td>
    </tr>
    <tr>
      <th>20</th>
      <td>Gifu</td>
      <td>35.39111</td>
      <td>136.72222</td>
    </tr>
    <tr>
      <th>21</th>
      <td>Shizuoka</td>
      <td>34.97694</td>
      <td>138.38306</td>
    </tr>
    <tr>
      <th>22</th>
      <td>Nagoya</td>
      <td>35.18028</td>
      <td>136.90667</td>
    </tr>
    <tr>
      <th>23</th>
      <td>Tsu</td>
      <td>34.73028</td>
      <td>136.50861</td>
    </tr>
    <tr>
      <th>24</th>
      <td>Otsu</td>
      <td>35.00444</td>
      <td>135.86833</td>
    </tr>
    <tr>
      <th>25</th>
      <td>Kyoto</td>
      <td>35.02139</td>
      <td>135.75556</td>
    </tr>
    <tr>
      <th>26</th>
      <td>Osaka</td>
      <td>34.68639</td>
      <td>135.52000</td>
    </tr>
    <tr>
      <th>27</th>
      <td>Kobe</td>
      <td>34.69139</td>
      <td>135.18306</td>
    </tr>
    <tr>
      <th>28</th>
      <td>Nara</td>
      <td>34.68528</td>
      <td>135.83278</td>
    </tr>
    <tr>
      <th>29</th>
      <td>Wakayama</td>
      <td>34.22611</td>
      <td>135.16750</td>
    </tr>
    <tr>
      <th>30</th>
      <td>Tottori</td>
      <td>35.50361</td>
      <td>134.23833</td>
    </tr>
    <tr>
      <th>31</th>
      <td>Matsue</td>
      <td>35.47222</td>
      <td>133.05056</td>
    </tr>
    <tr>
      <th>32</th>
      <td>Okayama</td>
      <td>34.66167</td>
      <td>133.93500</td>
    </tr>
    <tr>
      <th>33</th>
      <td>Hiroshima</td>
      <td>34.39639</td>
      <td>132.45944</td>
    </tr>
    <tr>
      <th>34</th>
      <td>Yamaguchi</td>
      <td>34.18583</td>
      <td>131.47139</td>
    </tr>
    <tr>
      <th>35</th>
      <td>Tokushima</td>
      <td>34.06583</td>
      <td>134.55944</td>
    </tr>
    <tr>
      <th>36</th>
      <td>Takamatsu</td>
      <td>34.34028</td>
      <td>134.04333</td>
    </tr>
    <tr>
      <th>37</th>
      <td>Matsuyama</td>
      <td>33.84167</td>
      <td>132.76611</td>
    </tr>
    <tr>
      <th>38</th>
      <td>Kochi</td>
      <td>33.55972</td>
      <td>133.53111</td>
    </tr>
    <tr>
      <th>39</th>
      <td>Fukuoka</td>
      <td>33.60639</td>
      <td>130.41806</td>
    </tr>
    <tr>
      <th>40</th>
      <td>Saga</td>
      <td>33.24944</td>
      <td>130.29889</td>
    </tr>
    <tr>
      <th>41</th>
      <td>Nagasaki</td>
      <td>32.74472</td>
      <td>129.87361</td>
    </tr>
    <tr>
      <th>42</th>
      <td>Kumamoto</td>
      <td>32.78972</td>
      <td>130.74167</td>
    </tr>
    <tr>
      <th>43</th>
      <td>Oita</td>
      <td>33.23806</td>
      <td>131.61250</td>
    </tr>
    <tr>
      <th>44</th>
      <td>Miyazaki</td>
      <td>31.91111</td>
      <td>131.42389</td>
    </tr>
    <tr>
      <th>45</th>
      <td>Kagoshima</td>
      <td>31.56028</td>
      <td>130.55806</td>
    </tr>
    <tr>
      <th>46</th>
      <td>Naha</td>
      <td>26.21250</td>
      <td>127.68111</td>
    </tr>
  </tbody>
</table>
</div>


都市間の位置関係を図示するとこうなります。

```python
%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 8))
plt.scatter(japan['Latitude'], japan['Longitude'])
for city, x, y in zip(japan['Town'], japan['Latitude'], japan['Longitude']):
    plt.text(x, y, city, alpha=0.5, size=12)
plt.grid()
```

![output_3_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/23ed3f00-153c-b053-0dc6-df01d0e66bfa.png)


# simanneal パッケージ

焼きなまし法を解くためのパッケージsimannealをインストールします。


```python
!pip install simanneal
```

このパッケージには例題として巡回セールスマン問題を解くためのクラスがすでに用意してあるので、そのまま使えます。


```python
from simanneal import Annealer
class TravellingSalesmanProblem(Annealer):

    """Test annealer with a travelling salesman problem.
    """

    # pass extra data (the distance matrix) into the constructor
    def __init__(self, state, distance_matrix):
        self.distance_matrix = distance_matrix
        super(TravellingSalesmanProblem, self).__init__(state)  # important!

    def move(self):
        """Swaps two cities in the route."""
        # no efficiency gain, just proof of concept
        # demonstrates returning the delta energy (optional)
        initial_energy = self.energy()

        a = random.randint(0, len(self.state) - 1)
        b = random.randint(0, len(self.state) - 1)
        self.state[a], self.state[b] = self.state[b], self.state[a]

        return self.energy() - initial_energy

    def energy(self):
        """Calculates the length of the route."""
        e = 0
        for i in range(len(self.state)):
            e += self.distance_matrix[self.state[i-1]][self.state[i]]
        return e
```

# 距離行列

距離行列は、たとえば scipy を使ってこのように計算できますが

```python
import numpy as np
from scipy.spatial import distance

mat = japan[['Latitude', 'Longitude']].values
dist_mat = distance.cdist(mat, mat, metric='euclidean') # ユークリッド距離
```

simanneal で使うためには、このように整形しないといけないようです。

```python
distance_matrix = {}
for i, town in enumerate(japan['Town']):
    if town not in distance_matrix.keys():
        distance_matrix[town] = {}
    for j, town2 in enumerate(japan['Town']):
        distance_matrix[town][town2] = dist_mat[i][j]
```

# 初期状態のセット

都市をランダムな順番に並べて、それを「初期状態の巡回路」とします。

```python
import random
init_state = list(japan['Town'])
random.shuffle(init_state)
```

それを巡回セールスマン問題を解くクラスにセット。

```python
tsp = TravellingSalesmanProblem(init_state, distance_matrix)
```

# 計算開始

```python
tsp.set_schedule(tsp.auto(minutes=0.2))
tsp.copy_strategy = "slice"
state, e = tsp.anneal()
```

# 計算結果

巡回路が出力されました。


```python
state
```




    ['Otsu',
     'Kyoto',
     'Nara',
     'Osaka',
     'Kobe',
     'Tottori',
     'Matsue',
     'Hiroshima',
     'Yamaguchi',
     'Fukuoka',
     'Saga',
     'Nagasaki',
     'Naha',
     'Kagoshima',
     'Miyazaki',
     'Kumamoto',
     'Oita',
     'Matsuyama',
     'Kochi',
     'Okayama',
     'Takamatsu',
     'Tokushima',
     'Wakayama',
     'Tsu',
     'Gifu',
     'Nagoya',
     'Shizuoka',
     'Kofu',
     'Maebashi',
     'Niigata',
     'Akita',
     'Aomori',
     'Sapporo',
     'Morioka',
     'Sendai',
     'Yamagata',
     'Fukushima',
     'Utsunomiya',
     'Mito',
     'Chiba',
     'Yokohama',
     'Tokyo',
     'Saitama',
     'Nagano',
     'Toyama',
     'Kanazawa',
     'Fukui']


# 指定した都市からの巡回路

得られた巡回路を、たとえば、東京からスタートの巡回路にします。

```python
while state[0] != 'Tokyo':
        state = state[1:] + state[:1]  # rotate NYC to start

print()
print("%i mile route:" % e)
print(" ➞  ".join(state))
```

    
    56 mile route:
    Tokyo ➞  Saitama ➞  Nagano ➞  Toyama ➞  Kanazawa ➞  Fukui ➞  Otsu ➞  Kyoto ➞  Nara ➞  Osaka ➞  Kobe ➞  Tottori ➞  Matsue ➞  Hiroshima ➞  Yamaguchi ➞  Fukuoka ➞  Saga ➞  Nagasaki ➞  Naha ➞  Kagoshima ➞  Miyazaki ➞  Kumamoto ➞  Oita ➞  Matsuyama ➞  Kochi ➞  Okayama ➞  Takamatsu ➞  Tokushima ➞  Wakayama ➞  Tsu ➞  Gifu ➞  Nagoya ➞  Shizuoka ➞  Kofu ➞  Maebashi ➞  Niigata ➞  Akita ➞  Aomori ➞  Sapporo ➞  Morioka ➞  Sendai ➞  Yamagata ➞  Fukushima ➞  Utsunomiya ➞  Mito ➞  Chiba ➞  Yokohama

# 図示



```python
plt.figure(figsize=(10, 10))
Xs = []
Ys = []
for i in range(len(state)):
    Xs.append(list(japan[japan['Town'] == state[i]].iloc[:, 2])[0])
    Ys.append(list(japan[japan['Town'] == state[i]].iloc[:, 1])[0])

plt.plot(Xs, Ys)
for city, x, y in zip(japan['Town'], japan['Latitude'], japan['Longitude']):
    plt.text(x, y, city, alpha=0.5, size=12)
```

![output_14_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/733246d2-069f-7bbe-4040-d72fd39e48c4.png)

なるほど。最適解ではなさそうですが、それに近いものは得られているようですね。