python實作大規模領域搜尋(LNS)求解旅行商問題(TSP)

大規模領域搜尋算法

參考《Handbook of Metaheuristics (Third Edition)》中的Large neighborhood search章節, 建議直接閱讀英文原版

LNS定義

大規模領域搜尋(LNS) 屬于超大領域搜尋(Very Large-Scale Neighborhood Search, VLNS)的一類，随着算例規模的增大，鄰域搜尋算法的鄰域規模呈指數增長或者當鄰域太大而不能在實際中明确搜尋時 (the neighborhood it searches grows exponentially with the instance size or if the neighborhood is simply too large to be searched explicitly in practice)，我們把這類鄰域搜尋算法(Neighborhood Search, NS)歸類于VLNS;

LNS領域

• 鄰域搜尋算法關鍵在于鄰域結構的選擇，即鄰域定義的方式。通常來講，鄰域越大，局部最優解就越好，獲得的全局最優解就越好。同時，鄰域越大，每次疊代搜尋鄰域所需的時間也越長。
• 在大規模鄰域搜尋算法中，鄰域由一種破壞（destroy）和一種修複（repair）算子隐式定義（the neighborhood is defined implicitly by a destroy and a repair method）。destroy算子會破壞目前解的一部分（變成不可行解），repair算子會對被破壞的解進行重建（重新變成可行解），相當于一個領域動作變換動作。破壞算子通常包含随機性的元素，以便在每次調用destroy方法時破壞解的不同部分（The destroy method typically contains an element of stochasticity such that different parts of the solution are destroyed in every invocation of the method）。
• 解 x x x的鄰域 N ( x ) N(x) N(x)就可以定義為：首先利用destroy算子破壞解 x x x，然後利用repair算子重建解 x x x，進而得到的一系列解的集合（The neighborhood N(x) of a solution x is then defined as the set of solutions that can be reached by first applying the destroy method and then the repair method）。
• LNS算法不會搜尋一個解的整個鄰域（entire neighborhood），而隻是對該鄰域進行采樣(samples)搜尋; (按照算法流程，隻定義一種detory和repair方式，雖然detory裡有随機部分，但也不包含整個領域)

LNS架構

• 僞代碼

在LNS中隻定義一種破壞和修複算子，是以破壞和修複算子的定義很重要，直接決定了能搜尋到的目前解x的領域。如果領域太小，那麼局部最優解的品質就一般，因為有些空間沒搜尋到。如果領域太大，又缺少帶點啟發式資訊方向搜尋。是以Ropke在論文《An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows》定義了多種破壞和修複算子，提出來自适應大規模領域搜尋（Adaptive Large Neighborhood Search, ALNS）的架構，有空總結寫一篇關于ALNS的文章。

旅行商問題TSP

• 旅行商問題（Travelling Salesman Problem, TSP），通俗而言它是指對于給定的一系列城市和每對城市之間的距離，找到通路每一座城市僅一次并回到起始城市的最短回路。建立不同的模組化方式會有不同的求解方式，比如Dantzig-Fulkerson-Johnson模型、Miller-Tucker-Zemlin模型、Commodity Flow、最短路徑等；(這裡不再贅述，注意TSP問題各種形式的變種)

python代碼示例及結果

• python代碼

import random
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']

class Solution:
    # 使用符号編号表示一個通路路徑route
    def __init__(self):
        self.route = []
        self.cost = 0  # 解對應的總成本

class Lns_tsp(object):
    def __init__(self, distance, num_node):
        self.distance = distance
        self.num_node = num_node

    def get_route_cost(self, route):
        # 計算成本的函數
        cost = 0
        for i in range(1, len(route)):
            cost += self.distance[route[i-1]][route[i]]
        return cost

    def destroy_operator(self, solution, num_destroy):
        # 破壞算子: 随機選擇num_destroy個不重複的破壞點（即删除num_destroy個城市）
        destroy_node_bank = []  # 儲存被删除的城市節點
        while len(destroy_node_bank) < num_destroy:
            n = random.randint(0, self.num_node-1)
            while n in destroy_node_bank:
                n = random.randint(0, self.num_node-1)
            destroy_node_bank.append(n)
            solution.route.remove(n)
        return solution, destroy_node_bank

    def repair_operator(self, solution, destroy_node_bank):
        # 修複算子: 貪婪插入，插入到成本最小的位置
        for n in destroy_node_bank:
            #計算将n插入各個位置的成本
            insert_list = np.full(len(solution.route), 0)
            for i in range(0, len(solution.route)):
                insert_list[i] = self.distance[solution.route[i-1]][n] + self.distance[n][solution.route[i]] - self.distance[solution.route[i]][solution.route[i-1]]

            greedy_index = np.where(insert_list == min(insert_list))[0][0]

            solution.route.insert(greedy_index, n)
        return solution

def plot_best_vales_iteration(best_values_record):
    # 繪制最優解随着疊代變化的趨勢
    plt.figure()
    plt.plot([i+1 for i in range(len(best_values_record))], best_values_record)
    plt.xlabel('疊代次數')
    plt.ylabel('最優值')
    plt.show()

def plot_route(route, city_location):
    plt.figure()
    # 繪制散點
    x = np.array(city_location)[:, 0]  # 橫坐标
    y = np.array(city_location)[:, 1]  # 縱坐标
    plt.scatter(x, y, color='r')
    # 繪制城市編号
    for i, txt in enumerate(range(1, len(city_location) + 1)):
        plt.annotate(txt, (x[i], y[i]))
    # 繪制方向
    x0 = x[route]
    y0 = y[route]
    for i in range(len(city_location) - 1):
        plt.quiver(x0[i], y0[i], x0[i + 1] - x0[i], y0[i + 1] - y0[i], color='b', width=0.005, angles='xy', scale=1,
                   scale_units='xy')
    plt.quiver(x0[-1], y0[-1], x0[0] - x0[-1], y0[0] - y0[-1], color='b', width=0.005, angles='xy', scale=1,
               scale_units='xy')

    plt.xlabel('橫坐标')
    plt.ylabel('縱坐标')
    plt.show()

if __name__ == '__main__':
    ############## 算例和參數設定 ############################
    # 城市節點的位置資訊，一行代表一個城市的橫坐标及縱坐标
    city_location = [[ 94,  99],
           [ 66,  67],
           [ 14,  78],
           [ 95,  56],
           [ 68,   9],
           [ 26,  20],
           [ 51,  67],
           [ 39,  39],
           [  5,  55],
           [ 12,  33],
           [ 55,  85],
           [ 98,  46],
           [ 36,  39],
           [ 65, 100],
           [ 57,  89],
           [ 88,  24],
           [ 53,  96],
           [ 91,  41],
           [ 32,  69],
           [ 38,  38],
           [ 38,  39],
           [ 85, 100],
           [  7,  37],
           [ 85,  96],
           [ 89,  48],
           [ 85,  35],
           [ 32,  29],
           [ 31,  25],
           [ 20,  17],
           [ 75,  21],
           [ 74,  29],
           [  6,  32],
           [ 20,  81],
           [ 62,   1],
           [ 11,  48],
           [  1,  69],
           [ 99,  70],
           [ 20,  27],
           [ 25,  42],
           [  6,  31],
           [ 78,  24],
           [ 42,  39],
           [ 83,  30],
           [ 94,  10],
           [ 90,  37],
           [ 76,  73],
           [  9,  56],
           [ 39,  33],
           [ 74,  15],
           [ 77,  14]]

    num_node = len(city_location)  # 城市節點的數量
    iter_num = 300  # 疊代次數
    random.seed(3)  # 随機種子
    num_destroy = int(num_node*0.2) # 破壞程度

    # 計算距離成本矩陣 distance, 直接使用歐式距離
    distance = np.full((num_node, num_node), 0)
    for i in range(num_node):
        for j in range(num_node):
            distance[i][j] = ((city_location[i][0]-city_location[j][0])**2+(city_location[i][1]-city_location[j][1])**2)**0.5

    ############## 産生初始解 ############################
    solution = Solution()
    solution.route = [i for i in range(num_node)]  # 按照節點編号依次相連構成初始解也可随機産生
    lns = Lns_tsp(distance, num_node)
    solution.cost = lns.get_route_cost(solution.route) # 計算初始解對應的目标成本
    best_solution = deepcopy(solution)  # 初始化最優解=初始解
    best_values_record = [0 for i in range(iter_num)]  # 初始化儲存最優解的集合

    ############## 執行LNS ############################
    for n_gen in range(iter_num):
        tem_solution = deepcopy(solution)
        # 執行破壞修複算子，得到臨時解
        tem_solution, destroy_node_bank = lns.destroy_operator(tem_solution, num_destroy)
        tem_solution = lns.repair_operator(tem_solution, destroy_node_bank)
        # 計算臨時解的目标值
        tem_solution.cost = lns.get_route_cost(tem_solution.route)

        # 接受标準：如果臨時解比目前解好，直接接受；且更新最優解
        if tem_solution.cost < best_solution.cost:
            solution = deepcopy(tem_solution)
            best_solution = deepcopy(tem_solution)
        best_values_record[n_gen] = best_solution.cost

    ############## 繪制結果 ############################
    plot_best_vales_iteration(best_values_record)
    plot_route(best_solution.route, city_location)
           

• 結果