【启发式算法】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))