【啟發式算法】Python實作禁忌搜尋算法求解TSP問題

禁忌搜尋

算法參數

• 禁忌表：tabu
• 禁忌長度：tabulen
• 特赦值：spe
• 疊代次數：Times

提醒：在對清單進行操作時，注意python的指派、淺拷貝、深拷貝

算法效果分析

1. 禁忌長度越大，鄰域搜尋時間變短（因為都被禁了），程式運作時間越短，結果的可信度相對較低
2. 特赦值越大，優于全局最優解的新解容易被特赦（容易被接受），可能會導緻陷入局部最優解；而特赦值過小可能會失去接受最優解的機會

算法步驟

1. 資料初始化：計算距離矩陣；随機得到初始解（起點為0城市）；編寫距離計算函數；初始化第一代的種群

2. 禁忌搜尋算法的核心：算法疊代Times次，每次使用雙重循環周遊所有鄰域（即所有可行解），本文采用的政策是順序周遊除起點以外的所有城市，判斷合法性，兩兩交換。

   2.1. 通過鄰域移動（即兩兩交換政策）得到新解并計算新解的路徑距離

   2.2. 如果新解優于全局最優解，且其禁忌長度小于設定的特赦值，則記錄交換城市的下标、接收新解并更新全局最優解

   2.3. 如果新解隻優于鄰域中的最優解，且禁忌長度為0（不在禁忌表中），則記錄交換城市的下标、接受新解并更新鄰域中的最優解

3. 鄰域周遊完成之後，更新目前解、更新禁忌表、重置兩個交換城市的禁忌長度
4. 直到疊代Times次，算法結束

# -*- coding: utf-8 -*-
import numpy as np
import random
import math
import matplotlib.pyplot as plt
import time
import copy
"""
Desicribe:禁忌搜尋解TSP問題
"""
# 原始資料
coordinates = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
                        [880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],
                        [1605.0, 620.0], [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],
                        [725.0, 370.0], [145.0, 665.0], [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],
                        [300.0, 465.0], [520.0, 585.0], [480.0, 415.0], [835.0, 625.0], [975.0, 580.0],
                        [1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], [410.0, 250.0],
                        [420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
                        [685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],
                        [475.0, 960.0], [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],
                        [830.0, 485.0], [1170.0, 65.0], [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],
                        [1340.0, 725.0], [1740.0, 245.0]])
gl_num = coordinates.shape[0]  # 城市個數
coord_x = coordinates[:, 0]  # 城市橫坐标
coord_y = coordinates[:, 1]  # 城市縱坐标
gl_dist = np.zeros((gl_num, gl_num))  # 距離矩陣
gl_nowRoute = list(range(0, gl_num))  # 目前路徑
gl_nowDistValue = 0.0  # 目前路徑長度
bestRoute = [0] * gl_num  #最優路徑
bestDistValue = 0.0  # 最優路徑長度
# 禁忌搜尋參數
tabu = np.zeros((gl_num, gl_num))  # 禁忌表
tabulen = 8  # 禁忌長度
spe = 5  # 特赦值
Times = 100  #疊代次數


def initDist():
    """
    初始化距離矩陣，dist[i][j]：城市i與城市j之間的歐式距離
    :return:
    """
    global gl_dist
    for i in range(gl_num):
        xi = coord_x[i]
        yi = coord_y[i]
        for j in range(i, gl_num):
            xj = coord_x[j]
            yj = coord_y[j]
            gl_dist[i][j] = gl_dist[j][i] = math.sqrt((xi - xj) ** 2 + (yi - yj) ** 2)

def getDist(route):
    """
    計算路徑長度
    :param route:
    :return:
    """
    distValue = 0.0
    for i in range(0, gl_num - 1):
        distValue += gl_dist[route[i]][route[i + 1]]
    distValue += gl_dist[route[gl_num - 1]][0]  # 從最後一個城市回到起點的距離
    return distValue


def cop(a, b):  # 把b數組的值指派a數組，深拷貝
    for i in range(gl_num):
        a[i] = b[i]


def init():
    """
    初始化初始解和初始路徑長度
    :return:
    """
    global gl_nowRoute, gl_nowDistValue, bestDistValue, bestRoute
    initDist()
    temp_list = gl_nowRoute[1:]  # 資料切片，左閉右開，擷取除起點終點以外的所有元素
    np.random.shuffle(temp_list)
    gl_nowRoute = gl_nowRoute[:1] + temp_list  # 切片的拼接

    # gl_nowRoute = [0, 21, 30, 17, 2, 16, 20, 41, 6, 1, 29, 22, 19, 49, 28, 15, 45, 43, 33, 38, 39, 37, 14, 4, 5, 3, 24, 11, 27, 26, 25, 46, 12, 13, 51, 10, 50, 32, 42, 9, 8, 7, 40, 18, 44, 31, 48, 35, 47, 23, 36, 34]
    gl_nowDistValue = getDist(gl_nowRoute)
    print("初始路徑：", gl_nowRoute)
    print("初始路徑長度：", gl_nowDistValue)
    cop(bestRoute, gl_nowRoute)
    bestDistValue = gl_nowDistValue  # 将初始解視作最優路徑和最優解


def solve():
    global gl_nowRoute, bestRoute, bestDistValue
    tempResult = [0] * gl_num  # 中間變量記錄交換結果
    index1 = 0
    index2 = 0  # 記錄交換城市的下标
    tempRoute = copy.deepcopy(gl_nowRoute)
    tempDistValue = gl_nowDistValue  # 暫存鄰域中的最優路徑與最優路徑長度
    # 搜尋所有鄰域
    for i in range(1, gl_num):  # 順序周遊除起點以外的所有城市，判斷合法性，兩兩交換
        for j in range(1, gl_num):
            if (i + j) >= gl_num:
                break
            if i == j:
                continue
            cop(tempResult, gl_nowRoute)
            tempResult[i], tempResult[i + j] = tempResult[i + j], tempResult[i]  # 交換
            tempValue = getDist(tempResult)
            # 如果優于全局最優且禁忌長度小于特赦值
            if (tempValue <= bestDistValue) & (tabu[i][i + j] < spe):
                # 接收新解并更新全局最優解
                cop(bestRoute, tempResult)
                bestDistValue = tempValue
                index1 = i
                index2 = i + j
                cop(tempRoute, tempResult)
                tempDistValue = tempValue
            # 如果優于鄰域中的最優解且禁忌長度為0則
            elif (tabu[i][i + j] == 0) & (tempValue < tempDistValue):
                # 接受新解
                cop(tempRoute, tempResult)
                tempDistValue = tempValue
                index1 = i
                index2 = i + j

    cop(gl_nowRoute, tempRoute)  # 更新目前解
    for i in range(gl_num):  # 更新禁忌表
        for j in range(gl_num):
            if tabu[i][j] > 0:
                tabu[i][j] -= 1
    tabu[index1][index2] = tabulen  # 重置兩個交換城市的禁忌長度


def draw(route):
    """
    畫出最優路徑
    :param route:
    :return:
    """
    x = [0] * (gl_num + 1)
    y = [0] * (gl_num + 1)
    for i in range(gl_num):
        index = route[i]
        x[i] = coordinates[index][0]
        y[i] = coordinates[index][1]
    x[gl_num] = coordinates[0][0]
    y[gl_num] = coordinates[0][1]
    plt.plot(x, y, c='r', marker='*')
    plt.show()


if __name__ == "__main__":
    start = time.time()
    init()
    for i in range(Times):
        solve()
    end = time.time()
    print("最優路徑：", bestRoute)
    print("最優路徑長度：", bestDistValue)
    draw(bestRoute)
    print('Running time: %s Seconds' % (end - start))