NP困難な最適化問題として知られる巡回セールスマン問題を、複数の近似解法と厳密解法を使って解いてみました。
使用した解法は以下の４個です。
- 整数計画問題(MTZ制約を使用、厳密解法)
- 遺伝的アルゴリズム(近似解法)
- 焼きなまし法(近似解法)
- 2-opt(山登り法として実装、近似解法)
実装はpythonで行いました。使用したモジュールは以下の通りです。
- pulp（整数計画問題用）
- deap（遺伝的アルゴリズム）
今回使用していないがTSPを解くのに活用できるモジュールは以下の通りです。
- simanneal（焼きなまし法）
- OR-Tools（Googleの最適化ソルバー）
- scipy.spatial Convexhull（凸包を作ってくれる）
適宜ソースコードを使いながら、それぞれについて順番に解説していきます。
コードの全体はこちらにあります。
整数計画問題として厳密解を求める
今回はMTZ条件を使って実装しました。
定式化
訪問する都市の総数を$n$、都市間の距離行列を$D$とします。
都市$i$から都市$j$に向かうかどうかを表現する変数を$x_{ij}\in\{0,1\}$、都市$i$の訪問順序を表現する整数変数を$u_i\in\{0,\dots,n-1\}$とします。
この時、必要な制約式は以下のようになります。
- $\sum_{i}x_{ij}=1\quad(\forall j)$・・・移動もとの都市は一つだけ
- $\sum_{j}x_{ij}=1\quad(\forall i)$・・・移動先の都市は一つだけ
- $u_i-u_j+nx_{ij}\leq n-1$・・・巡回路除去制約（MTZ条件）
この制約のもとで、$\sum_{i,j}D_{ij}x_{ij}$を最小化します。
実装
ソースコードの一部分を抜粋します。
def __init__(self, ncity: int, D: List[float], MTZ_level: int =0) -> None:
self.ncity = ncity
self.D = D
self.cities = list(range(ncity))
# 変数を定義する
self.x = pulp.LpVariable.dicts("x", (self.cities, self.cities), cat="Binary")# 都市iから都市jに向かうかを表す0-1変数
self.u = pulp.LpVariable.dicts("u", self.cities, cat="Integer", lowBound=0, upBound=ncity-1)# 訪問順序を表す
# 問題を宣言する
self.problem = pulp.LpProblem("TSP_IP")
# 目的関数を入れる
self.problem += pulp.lpSum(self.D[i][j]*self.x[i][j] for i in self.cities for j in self.cities)
for i in self.cities:# 制約を入れる
self.problem += pulp.lpSum(self.x[i][j] for j in self.cities) == 1# 移動先は一つの都市
self.problem += pulp.lpSum(self.x[j][i] for j in self.cities) == 1# 移動元は一つの都市
self.problem += self.x[i][i] == 0# 同じ都市に止まることはNG
self.set_MTZ(level=MTZ_level)# MTZ制約を入れる
#変数が取る値を調整
self.problem += self.u[self.cities[0]] == 0
for i in self.cities[1:]:
self.problem += self.u[i] >= 1
def set_MTZ(self, level: int) -> None:
# MTZ制約を入れる
# level=0が上で紹介したもの、その他は不等式の範囲を狭くすることでより条件を強めたもの
if level == 0:
for i in self.cities:
for j in self.cities[1:]:# 部分巡回除去制約(MTZ条件)
if i == j:
continue
self.problem += self.u[i] - self.u[j] + self.ncity*self.x[i][j] <= self.ncity - 1
elif level >= 1:
for i in self.cities[1:]:
for j in self.cities[1:]:# 部分巡回除去制約ちょっと強め(MTZ条件)
if i == j:
continue
# self.problem += self.u[i] + 1 - (self.ncity - 1)*(1 - self.x[i][j]) + (self.ncity - 3)*self.x[j][i] <= self.u[j]
self.problem += self.u[i] - self.u[j] + (self.ncity - 1)*self.x[i][j] + (self.ncity - 3)*self.x[j][i] <= self.ncity - 2
if level >= 2:# 部分巡回除去制約、強め(MTZ条件)
for i in self.cities[1:]:
self.problem += 2 - self.x[0][i] + (self.ncity - 3)*self.x[i][0] <= self.u[i]
self.problem += self.u[i] <= (self.ncity - 1) - (1 - self.x[i][0]) - (self.ncity - 3)*self.x[0][i]最適化を実行する場合は`problem.solve()`とすれば良いです。
コード全体はこちら。
2-opt近傍を使って山登り法で近似解を求める
2-opt近傍とは、巡回路の中で隣接していない都市を二つ選び、その順番を入れ替えるような操作を指します。初期解をランダムに作成し、2-optによって逐次改善していきます。目的関数値が改善した場合に限り都市の入れ替えを行います。
実装
計算量を削減するため、都市の入れ替え操作による目的関数値の差分のみを見て、入れ替えを実行するかどうか決定しています。
今回は入れ替える都市をランダムに選択しています。
操作によって目的関数が改善するかどうかの判定はswap関数の中で行なっています。
コード全体はこちらです。
def swap(self, i: int, j: int) -> None:
_swap_cost = self.swap_cost(i, j)
if _swap_cost < 0:
self.current_tour[i], self.current_tour[j] = self.current_tour[j], self.current_tour[i]
self.current_obj += _swap_cost
self.best_tour = list(self.current_tour)
self.best_obj = self.current_obj
def solve_two_opt(self, MAXSTEP: int =100000, strategy: str=None) -> List[int]:
# 初期化
self.current_tour, self.current_obj = self.initial_tour(strategy=strategy)
self.best_tour, self.best_obj = list(self.current_tour), self.current_obj
for _ in range(MAXSTEP):
i, j = random.randint(0, self.ncity - 1), random.randint(0, self.ncity - 1)
self.swap(i, j)
return self.best_tour2-opt近傍を使って焼きなまし法で近似解を求める
先程の入れ替え操作を実行するかどうかの判定を確率的に行うようにしたものがこちらです。
焼きなまし法についてはこの記事で紹介しています。
実装
コード全体はこちらです。
一回の操作が終了するたびに温度パラメーターを減少させるinhomogeneousなアルゴリズムと、一定回数操作を行なった後で温度パラメーターを減少させるhomogeneousなアルゴリズムを実装しています。
def solve_simulated_annealing(self, T: float =1e5, T_final: float =1e2, C: float =0.995, MAXSTEP: int =100, strategy: str=None, _type: str ="inhom") -> List[int]:
self.current_tour, self.current_obj = self.initial_tour(strategy=strategy)
self.best_tour, self.best_obj = list(self.current_tour), self.current_obj
# homogeneous algorithm
# 収束が遅い
if _type == "hom":
while T > T_final:
for _ in range(MAXSTEP):
i, j = random.randint(0, self.ncity - 1), random.randint(0, self.ncity - 1)
self.swap(i, j, T)
if self.current_obj < self.best_obj:
self.best_tour = list(self.current_tour)
self.best_obj = self.current_obj
T *= C
# inhomogeneous algorithm
elif _type == "inhom":
while T > T_final:
i, j = random.randint(0, self.ncity - 1), random.randint(0, self.ncity - 1)
self.swap(i, j, T)
if self.current_obj < self.best_obj:
self.best_tour = list(self.current_tour)
self.best_obj = self.current_obj
T *= C
else:
raise NotImplementedError
return self.best_tour遺伝的アルゴリズムを使用する
遺伝的アルゴリズムとは、最適化問題の解を遺伝子に見立て、交配、選択、突然変異を繰り返すことでより適応度の高い（＝目的関数値のいい）解を求めようとするアルゴリズムです。
今回は各操作を以下のように行いました。
- 交配操作「二つの解の一部分のインデックスを交換する」
- 突然変異「解の要素をシャッフル（ランダムに並べ替える）」
- 選択「適応度の高い個体を選択」
モジュールはdeapを使用し、実装は公式のチュートリアルを参考に行いました。
実装
コード全体はこちらです。
注意点は、適応度を表現するための関数（下のコードだと_eval関数）の返り値がタプルでないと動作しないことです。
class TSPGA(TwoOpt):
def __init__(self, ncity: int, D: List[float], indpb: float =0.05, toursize: int =3, npop: int =300) -> None:
super().__init__(ncity, D)
creator.create("FitnessMin", base.Fitness, weights=(-1.0, ))
creator.create("Individual", list, fitness=creator.FitnessMin)
self.toolbox = base.Toolbox()
self.toolbox.register("random_init", self.initial_tour_ga, "random")
self.toolbox.register("individual", tools.initIterate, creator.Individual, self.toolbox.random_init)
self.toolbox.register("population", tools.initRepeat, list, self.toolbox.individual)
self.toolbox.register("mate", tools.cxPartialyMatched)
self.toolbox.register("mutate", tools.mutShuffleIndexes, indpb=indpb)
self.toolbox.register("select", tools.selTournament, tournsize=toursize)
self.toolbox.register("evaluate", self._eval)
self.pop = self.toolbox.population(n=npop)
def initial_tour_ga(self, strategy: str =None) -> list:
return self.initial_tour(strategy=strategy)[0]
def _eval(self, tour: List[int]) -> Tuple[float]:
return self.calc_score(tour),
def solve(self, cxpb: float =0.7, mutpb: float =0.2, ngen: int =40, hof: int=1) -> object:
if hof:
hof = tools.HallOfFame(hof)
algorithms.eaSimple(self.pop, self.toolbox, cxpb, mutpb, ngen, halloffame=hof, verbose=False)
return hof[0]
この記事は役に立ちましたか？
もし参考になりましたら、下記のボタンで教えてください。
コメント