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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159