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

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

巡回セールスマン問題

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

焼きなまし法

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

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

グラフ理論の基礎やグラフ理論の基礎をmatplotlibアニメーションでのデータを例に、巡回セールスマン問題を解いてみようと思います。

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>)

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

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

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

%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()

simanneal パッケージ

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

!pip install simanneal

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

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 を使ってこのように計算できますが

import numpy as np
from scipy.spatial import distance

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

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

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]

初期状態のセット

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

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

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

tsp = TravellingSalesmanProblem(init_state, distance_matrix)

計算開始

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

計算結果

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

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']

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

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

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

図示

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)

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