巡回セールスマン問題をいろいろな手法で解いてみる②

Last updated at 2024-11-10Posted at 2024-11-10

はじめに

前回記事に引き続き、巡回セールスマン問題を題材にし、最適化を行います。
巡回セールスマン問題をいろいろな手法で解いてみる①
前回記事ではpulpを用いて厳密解を確認しました。
今回は局所探索の一つの山登り法を実装して、厳密解と比較していきます。

近似解法

局所探索を実装する前に、近似解法について簡単にまとめます。
近似解法は以下の二つに大別されます。

構築型近似解法
改善型近似解法

どちらも近似解を出力しますが、アプローチが異なります。

構築型近似解法

構築型では、解を構造体としてとらえ、解要素を加えていくことで1つの解を構築していく手法です。
巡回セールスマン問題の場合ですと、初めに都市Aを訪れ、次に都市Bを訪れ、次に都市Cを訪れ...といった具合に経路を1つずつ決定していきます。
代表的なものとして、貪欲法が挙げられます。

改善型近似解法

改善型では、すでに出来上がった解の一部分を変更し、解を改善していく方法です。
巡回セールスマン問題の中でも解の構造はいくつかのパターンがあげられます。
最もメジャーなものの一つとして、訪れる都市の順番を順列で表現する方法があります。
順列の順番の一部を変更することで解を改善していきます。
局所探索法は、この中の山登り法をはじめ、
タブーサーチ、シミュレーテッドアニーリング法などが該当します。

山登り法

局所探索と山登り法

ここでは局所探索について、まとめます。
局所探索は以下のアルゴリズムの枠組みで実装されるアルゴリズムの総称です。

解を一つ生成する
現在の解の近傍操作を行い近傍解（解の一部分を変更した解）とする
一定の規則で近傍解と現在の解と入れ換える
終了条件を満たすまで 2. 以下を繰り返す

一定のルールによって、現在の解を生成された近傍解で更新していくことで、最適化を行っていきます。
中でも山登り法では生成された近傍解の中で最も評価値の高い解に更新することで最適化を行います。

山登り法の実装（交換近傍）

では、山登り法を実装していきます。
アルゴリズムの流れは以下になります。

ランダムに初期解を生成
終了条件（繰り返し数が上限になる）を満たすまで3～5から
現在の解をもとに近傍解を生成する
近傍解を評価する
現在の解を近傍解の中で最もよい解に更新する

近傍解の生成方法（近傍操作）はこれまでにいくつものパターンが考案されています。
また、近傍操作と解の表現方法は密接に関係しています。

今回は巡回セールスマン問題の解に順列表現を用います。
順列表現ではリスト構造のインデクスが若い順に巡る都市のインデクスが登録されています。
近傍操作は単純な交換近傍とします。
これはリスト構造体の2つの要素を取り出し交換するというものです。
今回は一つの解に対してランダムに30個の近傍解を生成します。

import math
import random
import time
from matplotlib import pyplot as plt
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 _ 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:
    #昇順
    def initial_route(self,instance):
        route = list(range(instance.n))
        return route
    def cal_cost(self,route,instance):
        sum_cost = 0
        n = instance.n
        for i in range(n - 1):
            sum_cost += instance.cost[route[i]][route[i + 1]]
        sum_cost += instance.cost[route[n - 1]][route[0]]
        return sum_cost
    def get_cost(self, index_i, index_j, route, instance):
        if index_i == instance.n:
            index_i = 0
        if index_j == instance.n:
            index_j = 0
        return instance.cost[route[index_i]][route[index_j]]
    #交換近傍
    def exchange_neighbor(self, index_i, index_j, route, instance):
        n = instance.n
        ### i,j:都市のID
        ### index_i,index_j:route上の都市i,jのインデクス
        ### delta:近傍操作による変化量
        delta = 0
        ### 隣接する場合
        if(index_i == index_j - 1 or index_i == index_j + n - 1):
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
        elif(index_i == index_j + 1 or index_j == index_i + n - 1):
            delta -= (self.get_cost(index_j - 1, index_j, route, instance) + self.get_cost(index_i, index_i + 1, route, instance))
            delta += (self.get_cost(index_j - 1, index_i, route, instance) + self.get_cost(index_j, index_i + 1, route, instance))
        ### 隣接しない場合
        else:
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_i, index_i + 1, route, instance))
            delta += (self.get_cost(index_j - 1, index_i, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
            delta -= (self.get_cost(index_j - 1, index_j, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_j, index_i + 1, route, instance))
        return delta
    def exchange_route(self, route, index_i, index_j):
        route[index_i], route[index_j] = route[index_j], route[index_i]
        return route
    #山登り探索
    def local_search(self, ini_route, iteration, instance, neighbor_ope, change_route, word):
        best = 999999
        local = self.cal_cost(ini_route,instance)
        best = min(best, local)
        #### search
        route = ini_route.copy()
        best_route = ini_route.copy()
        neighbors = dict()
        iteration_values = [0]*iteration
        fig_index = 0
        gen_fig(instance, best_route, 'Local'+word+'_fig/'+f'{fig_index:04}'+".png")
        for k in range(iteration):
            #近傍解を列挙
            neighbors = []
            for _ in range(30):
                index_i, index_j =random.sample(range(instance.n), 2)
                neighbors.append([neighbor_ope(index_i, index_j, route, instance), index_i, index_j])
            delta, index_i, index_j = neighbors[min(range(len(neighbors)), key=neighbors.__getitem__)]
            #最もよい近傍に遷移
            if(delta < 0):
                #解の更新
                route = change_route(route, index_i, index_j)
                local += delta
                best = local
                best_route = route.copy()
                fig_index += 1
                gen_fig(instance, best_route, 'Local'+word+'_fig/'+f'{fig_index:04}'+".png")
            iteration_values[k] = best
        result = dict()
        result['tour'] = str(best_route)
        result['cost'] = best
        result['iteration_values'] = iteration_values
        return result

    #結果のプロット
    def plot(self,results,iteration):
        fig, ax = plt.subplots()
        x = list(range(iteration))
        for i in range(len(results)):
            y = results[i]['iteration_values']
            ax.plot(x, y, label = "local search"+str(i))
        ax.set_xlabel('iteration')  # x軸ラベル
        ax.set_ylabel('cost')  # y軸ラベル
        ax.set_title('Shifting costs of solutions') # グラフタイトル
        ax.legend(loc = 0)
        plt.savefig("result.png")
        plt.clf()

def gen_fig(instance, path, fig_name):
    G = nx.Graph()
    G.add_nodes_from(range(instance.n))
    for i in range(instance.n):
        G.add_edge(path[i-1], path[i])
    nx.draw_networkx(G, pos=instance.points, node_color="c")
    plt.savefig(fig_name)
    plt.clf()

def main():
    problem_size = 30
    instance = Instance(problem_size)
    iteration = 300
    solve = Solve()
    ini_route1 = solve.initial_route(instance)
    ini_routes = [ini_route1]
    for ini_route in ini_routes:
        t = time.time()
        local_result1 = solve.local_search(ini_route, iteration, instance, solve.exchange_neighbor, solve.exchange_route, 'Exchange')
        print("交換近傍を用いた局所探索")
        print('目的関数値 : '+str(local_result1['cost']))
        print(time.time()-t,'(s)')
        print()
        
        solve.plot([local_result1],iteration)

if __name__ == '__main__':
    main()

計算結果（交換近傍）

今回のアルゴリズムで生成された近似解は526でした。
前回のpulpで生成した解（432）と比較すると、94のギャップがあります。

前回と同じく、解を可視化してみましょう。

解の遷移の様子も可視化してみましょう。

解は少しずつ改善されており、最適化はされているようです。
しかしながら、生成される解の精度はまだまだといった印象。

せっかくなので、ほかの近傍操作も実装して比較してみましょう。

山登り法の実装（2-opt近傍）

ここでは2-opt近傍を使った山登り法の実装を行います。
2-opt近傍の操作では、順列表現の解に対して、二つの要素を抽出し、その区間の要素すべての順序を反転させる操作になります。
こちらでもは一つの解に対してランダムに30個の近傍解を生成します。
実装では近傍操作とルートの更新で下記のコードを使用します。

class Solve:
    #2-opt近傍
    def two_opt_neighbor(self, index_i, index_j, route, instance):
        ### i,j:都市のID
        ### index_i,index_j:route上の都市i,jのインデクス
        ### index_i < index_j
        if index_i > index_j:
            index_i, index_j = index_j, index_i
        delta = 0
        if not (index_i == 0 and index_j == instance.n-1):
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
        return delta
    def two_opt_route(self, route, index_i, index_j):
        if index_i > index_j:
            index_i, index_j = index_j, index_i
        route[index_i:index_j+1] = reversed(route[index_i:index_j+1])
        return route

計算結果（2-opt近傍）

計算結果は450でした。

解を可視化するとこのような結果です。

遷移の様子はこちら。

まずまずですね。

おわりに

今回は局所探索法の一つである山登り法を実装しました。
その中で２つの近傍操作を実装し、それぞれの解をpulpで生成した厳密解と比較しました。
今回のケースでは2-opt近傍によるアプローチのほうが良い結果を示すようです。
ただ、厳密解と比較するとまだギャップがあります。
次回は今回実装した山登り法をベースにタブーサーチを実装します。

コード全体

長くなりますが、下記しておきます。

import math
import random
import time
from matplotlib import pyplot as plt
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 _ 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:
    #昇順
    def initial_route(self,instance):
        route = list(range(instance.n))
        return route
    def cal_cost(self,route,instance):
        sum_cost = 0
        n = instance.n
        for i in range(n - 1):
            sum_cost += instance.cost[route[i]][route[i + 1]]
        sum_cost += instance.cost[route[n - 1]][route[0]]
        return sum_cost
    def get_cost(self, index_i, index_j, route, instance):
        if index_i == instance.n:
            index_i = 0
        if index_j == instance.n:
            index_j = 0
        return instance.cost[route[index_i]][route[index_j]]
    #交換近傍
    def exchange_neighbor(self, index_i, index_j, route, instance):
        n = instance.n
        ### i,j:都市のID
        ### index_i,index_j:route上の都市i,jのインデクス
        ### delta:近傍操作による変化量
        delta = 0
        ### 隣接する場合
        if(index_i == index_j - 1 or index_i == index_j + n - 1):
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
        elif(index_i == index_j + 1 or index_j == index_i + n - 1):
            delta -= (self.get_cost(index_j - 1, index_j, route, instance) + self.get_cost(index_i, index_i + 1, route, instance))
            delta += (self.get_cost(index_j - 1, index_i, route, instance) + self.get_cost(index_j, index_i + 1, route, instance))
        ### 隣接しない場合
        else:
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_i, index_i + 1, route, instance))
            delta += (self.get_cost(index_j - 1, index_i, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
            delta -= (self.get_cost(index_j - 1, index_j, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_j, index_i + 1, route, instance))
        return delta
    def exchange_route(self, route, index_i, index_j):
        route[index_i], route[index_j] = route[index_j], route[index_i]
        return route
    #山登り探索
    def local_search(self, ini_route, iteration, instance, neighbor_ope, change_route, word):
        best = 999999
        local = self.cal_cost(ini_route,instance)
        best = min(best, local)
        #### search
        route = ini_route.copy()
        best_route = ini_route.copy()
        neighbors = dict()
        iteration_values = [0]*iteration
        fig_index = 0
        gen_fig(instance, best_route, 'Local'+word+'_fig/'+f'{fig_index:04}'+".png")
        for k in range(iteration):
            #近傍解を列挙
            neighbors = []
            for _ in range(30):
                index_i, index_j =random.sample(range(instance.n), 2)
                neighbors.append([neighbor_ope(index_i, index_j, route, instance), index_i, index_j])
            delta, index_i, index_j = neighbors[min(range(len(neighbors)), key=neighbors.__getitem__)]
            #最もよい近傍に遷移
            if(delta < 0):
                #解の更新
                route = change_route(route, index_i, index_j)
                local += delta
                best = local
                best_route = route.copy()
                fig_index += 1
                gen_fig(instance, best_route, 'Local'+word+'_fig/'+f'{fig_index:04}'+".png")
            iteration_values[k] = best
        result = dict()
        result['tour'] = str(best_route)
        result['cost'] = best
        result['iteration_values'] = iteration_values
        return result

    #2-opt近傍
    def two_opt_neighbor(self, index_i, index_j, route, instance):
        ### i,j:都市のID
        ### index_i,index_j:route上の都市i,jのインデクス
        ### index_i < index_j
        if index_i > index_j:
            index_i, index_j = index_j, index_i
        delta = 0
        if not (index_i == 0 and index_j == instance.n-1):
            delta -= (self.get_cost(index_i - 1, index_i, route, instance) + self.get_cost(index_j, index_j + 1, route, instance))
            delta += (self.get_cost(index_i - 1, index_j, route, instance) + self.get_cost(index_i, index_j + 1, route, instance))
        return delta
    def two_opt_route(self, route, index_i, index_j):
        if index_i > index_j:
            index_i, index_j = index_j, index_i
        route[index_i:index_j+1] = reversed(route[index_i:index_j+1])
        return route

    #結果のプロット
    def plot(self,results,iteration):
        fig, ax = plt.subplots()
        x = list(range(iteration))
        for i in range(len(results)):
            y = results[i]['iteration_values']
            ax.plot(x, y, label = "local search"+str(i))
        ax.set_xlabel('iteration')  # x軸ラベル
        ax.set_ylabel('cost')  # y軸ラベル
        ax.set_title('Shifting costs of solutions') # グラフタイトル
        ax.legend(loc = 0)
        plt.savefig("result.png")
        plt.clf()

def gen_fig(instance, path, fig_name):
    G = nx.Graph()
    G.add_nodes_from(range(instance.n))
    for i in range(instance.n):
        G.add_edge(path[i-1], path[i])
    nx.draw_networkx(G, pos=instance.points, node_color="c")
    plt.savefig(fig_name)
    plt.clf()

def main():
    problem_size = 30
    instance = Instance(problem_size)
    iteration = 300
    solve = Solve()
    ini_route1 = solve.initial_route(instance)
    ini_routes = [ini_route1]
    for ini_route in ini_routes:
        t = time.time()
        local_result1 = solve.local_search(ini_route, iteration, instance, solve.exchange_neighbor, solve.exchange_route, 'Exchange')
        print("交換近傍を用いた局所探索")
        print('目的関数値 : '+str(local_result1['cost']))
        print(time.time()-t,'(s)')
        print()
        
        t = time.time()
        local_result2 = solve.local_search(ini_route, iteration, instance, solve.two_opt_neighbor, solve.two_opt_route, 'TwoOpt')
        print("2-opt近傍を用いた局所探索")
        print('目的関数値 : '+str(local_result2['cost']))
        print(time.time()-t,'(s)')
        print()

        solve.plot([local_result1, local_result2],iteration)

if __name__ == '__main__':
    main()

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up