巡回セールスマン問題

巡回セールスマン問題#

巡回セールスマン問題とは、セールスマンが各都市を一度ずつ訪問し、出発点に戻る際の移動距離を最小化する問題です。本チュートリアルでは、日本の首都圏を舞台に巡回セールスマン問題をどのように解くかを説明します。

定式化#

以下では、\(n\) 都市の巡回セールスマン問題を考えることとします。また、巡回セールスマン問題はいくつかの線形定式化が知られていますが、本チュートリアルでは次のような2次定式化を採用するものとします。

決定変数#

決定変数 \(x_{i, t}\) を次のように定義するものとします。

\[\begin{split} x_{i,t} = \begin{cases} 1~\text{if the salesman visits city }i\text{ at time }t\\ 0~\text{otherwise} \end{cases} \end{split}\]
\[ i, t \in \{0, ..., n-1\} \]

すなわち、\(x_{i, t}\) はセールスマンが \(t\) 番目に都市 \(i\) を訪れた場合に1になり、そうでない時に0になるバイナリ変数です。

制約#

次に以下の2つの制約を考えます。

  1. 都市 \(i\) に訪れるのは1度だけ

\[ \sum_t x_{i, t} = 1, ~\forall i \]

  1. \(t\) 番目に訪れる都市は1つだけ

\[ \sum_i x_{i, t} = 1, ~\forall t \]

この2つの制約条件は \(x_{i, t}\) がセールスマンの都市を巡る順序を表現する上で必要なものになっています。なぜなら、セールスマンは同時に複数の都市を訪問することはできませんし、同じ都市に何度も行く必要はないからです。

以下の図は、この2つの制約条件を満たしながら \(x_{i, t}\) に値を割り当てたものになっています。このようなイメージで定式化を捉えると以降の説明が理解しやすくなるでしょう。

目的関数#

最後に目的関数を考えましょう。巡回セールスマン問題はセールスマンの移動距離を最小化する問題なので、次のように定式化することができます。

\[ \sum_{i,j,t} d_{i, j}x_{i, t}x_{j, (t+1)\mod n} \]

ここで \(d_{i, j}\) は都市 \(i\) と都市 \(j\) の距離です。 \(x_{i, t}x_{i, t+1}\) は \(t\) 番目に都市 \(i\) を訪れ、 \(t+1\) 番目に都市 \(j\) を訪れた時に値が1になり、そうでない時には値が0になります。上記の式で \(x_{j, (t+1)\mod n}\) に \(\mod n\) が含まれているのは、\(n\) 番目に訪れた都市 \(j\) から最初の都市に戻る移動を表現するためです。

\(n\)都市の巡回セールスマン問題#

これで定式化に必要なものが揃いました。\(n\) 都市の巡回セールスマン問題は以下のように数理モデルで表すことができます。

\[\begin{split} \begin{aligned} \min_x & \sum_{i, j, t} d_{i, j} x_{i,t} x_{j, (t+1) \mod n}\\ \mathrm{s.t.}~&\sum_i x_{i, t} = 1,~\forall t\\ &\sum_t x_{i, t} = 1, ~\forall i\\ &x_{i, t} \in \{0, 1\} \end{aligned} \end{split}\]

JijModelingによる定式化#

JijModelingを使って上述の数理モデルを実装してみましょう。まずは、数理モデルに現れる変数およびパラメーターを定義します。

import jijmodeling as jm

problem = jm.Problem("TSP")
N = problem.Natural("N")
d = problem.Float("d", shape=(N, N))
x = problem.BinaryVar("x", shape=(N, N))

ここで N は都市数で、d は 都市 \(i\) と都市 \(j\) の距離 \(d_{i, j}\) を要素に持つ2次元リストです。 そして、\(x_{i, t}\) に相当するバイナリ変数として x を定義しています。

目的関数の実装#

次に、目的関数を設定しましょう。

problem += jm.sum(
    jm.product(N, N, N),
    lambda i, j, t: d[i, j] * x[i, t] * x[j, (t + 1) % N],
)

制約条件の実装#

そして、2つの制約条件を定義しましょう。

problem += problem.Constraint(
    "one-city", lambda t: jm.sum(N, lambda i: x[i, t]) == 1, domain=N
)
problem += problem.Constraint(
    "one-time", lambda i: jm.sum(N, lambda t: x[i, t]) == 1, domain=N
)

ここでは "one-city" と名付けられている制約により、\(t\) 番目に訪れる都市は1つだけを意味する制約を定義しています。 そして"one-time" と名付けられている制約によって、都市 \(i\) に訪れるのは1度だけを意味する制約条件を設定しています。

これで数理モデルの実装は完了です。正しく数理モデルが実装されているかをLaTeX表示を通して確認してみましょう。

problem
\[\begin{split}\begin{array}{rl} \text{Problem}\colon &\text{TSP}\\\displaystyle \min &\displaystyle \sum _{i=0}^{N-1}{\sum _{j=0}^{N-1}{\sum _{t=0}^{N-1}{{d}_{i,j}\cdot {x}_{i,t}\cdot {x}_{j,\left(t+1\right)\bmod N}}}}\\&\\\text{s.t.}&\\&\begin{aligned} \text{one-city}&\quad \displaystyle \sum _{i=0}^{N-1}{{x}_{i,t}}=1\quad \forall t\;\text{s.t.}\;t\in \left\{0,\ldots ,N-1\right\}\\\text{one-time}&\quad \displaystyle \sum _{t=0}^{N-1}{{x}_{i,t}}=1\quad \forall i\;\text{s.t.}\;i\in \left\{0,\ldots ,N-1\right\}\end{aligned} \\&\\\text{where}&\\&\text{Decision Variables:}\\&\qquad \begin{alignedat}{2}x&\in \mathop{\mathrm{Array}}\left[N\times N;\left\{0, 1\right\}\right]&\quad &2\text{-dim binary variable}\\\end{alignedat}\\&\\&\text{Placeholders:}\\&\qquad \begin{alignedat}{2}d&\in \mathop{\mathrm{Array}}\left[N\times N;\mathbb{R}\right]&\quad &2\text{-dimensional array of placeholders with elements in }\mathbb{R}\\N&\in \mathbb{N}&\quad &\text{A scalar placeholder in }\mathbb{N}\\\end{alignedat}\end{array} \end{split}\]

インスタンスデータの準備#

次はパラメーターに入力するインスタンスデータを準備します。 ここでは日本の首都圏（関東平野）の都県を回る問題を考えてみましょう。

import numpy as np

# set the name list of traveling points
points = ['茨城県', '栃木県', '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県']

# set the latitudes and longitudes of traveling points
latlng_list = [
    [36.2869536, 140.4703384],
    [36.6782167, 139.8096549],
    [36.52198, 139.033483],
    [35.9754168, 139.4160114],
    [35.549399, 140.2647303],
    [35.6768601, 139.7638947],
    [35.4342935, 139.374753]
]

# make distance matrix
inst_N = len(points)
inst_d = np.zeros((inst_N, inst_N))
for i in range(inst_N):
    for j in range(inst_N):
        a = np.array(latlng_list[i])
        b = np.array(latlng_list[j])
        inst_d[i][j] = np.linalg.norm(a-b)

# normalize distance matrix
inst_d = (inst_d-inst_d.min()) / (inst_d.max()-inst_d.min())

geo_data = {'points': points, 'latlng_list': latlng_list}
instance_data = {'N': inst_N, 'd': inst_d}

JijZeptSolverで解く#

jijzept_solverを用いて、巡回セールスマン問題を解きましょう。

import jijzept_solver

instance = problem.eval(instance_data)
solution = jijzept_solver.solve(instance, solve_limit_sec=1.0)

解の可視化#

得られた解を用いて、セールスマンが都市を訪問する順序を可視化してみましょう。

import matplotlib.pyplot as plt

df = solution.decision_variables_df
x_indices = df[df["value"] == 1.0]["subscripts"]
tour = [index for index, _ in sorted(x_indices, key=lambda x: x[1])]
tour.append(tour[0])

position = np.array(geo_data['latlng_list']).T[[1, 0]]

plt.axes().set_aspect('equal')
plt.plot(*position, "o")
plt.plot(position[0][tour], position[1][tour], "-")
plt.show()

print(np.array(geo_data['points'])[tour])
../_images/bae14b1cc1854db0fcccf87e1bbf81a9c6def7a8ee8edb104aec4867fa9e0b85.png 
['埼玉県' '群馬県' '栃木県' '茨城県' '千葉県' '東京都' '神奈川県' '埼玉県']
目次