はじめに

巡回セールスマン問題を例題として、いろいろな解法を試してみたので、備忘録がてら、いくつか記事を作成してみます。

第1回は

巡回セールスマン問題の紹介
pulpによる厳密解の導出
までをまとめます。

巡回セールスマン問題とは

巡回セールスマン問題（TSP：traveling salesman problem）は、複数の決められた地点をセールスマンが各１度ずつ巡回して元に戻る場合、どのように回れば最も移動距離が少ないかという問題です。
この問題は計算複雑性理論においてNP困難に属する問題で、一般的に問題規模が大きくなると計算量が爆発的に大きくなる問題です。

定式化

さて、この問題を数式モデルとして定式化してみましょう。
まずは以下を定義します。

集合
- $N$：都市の集合
定数
- $c_{ij}$：都市 $i(i \in N)$から都市 $j(j \in N)$ $(i \neq j)$ の移動に要する移動コスト
決定変数
- $x_{ij}$：都市 $i$から都市 $j$に移動する場合 1 ，それ以外の場合 0 の値をとるバイナリ決定変数

定式化は以下のようになります。

$$
min \sum_{i}^N \sum_{j}^N c_{ij} x_{ij} \label{a} \tag{1}
$$
$$
s.t. \sum_{i}^N x_{ij} = 1, \ \ \ \ \ i = 1,2,...,n　\label{b} \tag{2}
$$
$$
s.t. \sum_{i}^N x_{ij} = 1, \ \ \ \ \ j = 1,2,...,n \label{c} \tag{3}
$$
$$
u_i+1-M(1-x_{ij}) \leq u_j, \ \ \ \ \ \forall i,j \label{d} \tag{4}
$$
$$
1 \leq u_j \leq n-1 \label{e} \tag{5}
$$
$$
x_{ij} \in {0,1}, \ \ \ \ \ \forall i,j \label{f} \tag{6}
$$

$(\ref{a})$式は目的関数が巡回路の距離の総和の最小化であることを示します。
$(\ref{b})$式は都市$i$を到着するのは 1 回だけであることを示します。
同様に、$(\ref{c})$式は都市$j$を出発するのは 1 回だけであることを示します。
$(\ref{d})$$(\ref{e})$式は部分巡回路を取り除くための制約です。
$(\ref{f})$式は決定変数 𝑥𝑖𝑗が離散変数条件であることを示します。

pulpによる最適化モデリング

さて、この定式化を使って、実際に問題を解いてみましょう。
最適化において、解法は大別して、厳密解法と近似解法に分けられます。
厳密解法はその問題で最も良い答えを保証して出力しますが、問題によっては解くのに時間がかかってしまいます。
一般的に問題規模が大きいほど答えの導出に時間がかかるとされます。
一方で、近似解法を用いると、少ない時間でそこそこの答えを出力することができます。
今回は無料で使用できるpyrhonライブラリのpulpを使い厳密解を出力します。
なお、pulpについてやモデリング方法の詳細は下記を参考にしました。
最適化におけるPython（PuLP版）

プログラム

都市数は30個に設定し、各都市はランダムで配置しました。
各都市間の距離はユークリッド距離で与えます。
せっかくですので、出力結果はnetworkXライブラリを使って画像出力しましょう。

main.py

import math
import random
import time
import pulp
from pulp import LpProblem, LpMinimize, LpVariable, lpSum, LpBinary, value
from matplotlib import pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np

class Instance:
    def __init__(self,n):
        random.seed(111)
        self.n = n
        self.points = []
        for i in range(n):
            self.points.append([random.randint(0,100),random.randint(0,100)])
        self.cost = [[0] * n for i in range(n)]
        for i in range(n):
            for j in range(i+1,n):
                # #ユークリッド距離
                self.cost[i][j] = int(math.sqrt(math.fabs(self.points[i][0]-self.points[j][0])**2 + math.fabs(self.points[i][1]-self.points[j][1])**2))
                self.cost[j][i] = int(math.sqrt(math.fabs(self.points[i][0]-self.points[j][0])**2 + math.fabs(self.points[i][1]-self.points[j][1])**2))

class Solve:
    #pulpを用いた厳密解法
    def Pulp(self,instance):
        n = instance.n
        cost = np.array(instance.cost)
        M = 10000
        # Create the problem
        tsp = LpProblem("TSP", LpMinimize)
        # Decision variables
        x = [[] for _ in range(n)]
        for i in range(n):
            for j in range(n):
                if i != j:
                    x[i].append(LpVariable(f"x_{i}_{j}", cat=LpBinary))
                else:
                    x[i].append(0)
        u = [LpVariable(f"u_{i}", lowBound=0, upBound=n-1) for i in range(n)]
        # Objective function
        tsp += lpSum(cost[i][j] * x[i][j] for i in range(n) for j in range(n))
        # Assignment constraints
        for i in range(n):
            tsp += lpSum(x[i][j] for j in range(n)) == 1
        for j in range(n):
            tsp += lpSum(x[i][j] for i in range(n)) == 1
        # Subtour-breaking constraints
        for i in range(n):
            for j in range(1, n):
                if i != j:
                    tsp += u[i] + 1 - u[j] <= M * (1 - x[i][j])
        # Solve the problem
        tsp.solve(pulp.PULP_CBC_CMD())
        result = dict()
        print("厳密解：目的関数値",value(tsp.objective))
        
        G = nx.Graph()
        G.add_nodes_from(range(instance.n))
        for i in range(n):
            for j in range(n):
                if i != j:
                    if value(x[i][j]) == 1:
                        G.add_edge(i, j)
        nx.draw_networkx(G, pos=instance.points, node_color="c")
        plt.savefig('exact_sol.png')
        plt.clf()
        return result
        

def main():
    problem_size = 30
    instance = Instance(problem_size)
    solve = Solve()
    
    t = time.time()
    print("pulpによる厳密解計算")
    pulp_result = solve.Pulp(instance)
    print(time.time()-t,'(s)')
    print()

if __name__ == '__main__':
    main()

計算結果

こちらのプログラムを実行した場合、目的関数が432.0の解が出力されます。
なお、この解は次のような順路です。

いい感じですね！

終わりに

今回は巡回セールスマン問題の定式化を紹介し、pulpによって問題を解いてみました。
次回は巡回セールスマン問題に対して、代表的なメタ・ヒューリスティクスの一つである局所探索を実装します。