优化算法 | 禁忌搜索算法求解TSP问题及其Python代码实现

导读

禁忌搜索算法（Tabu Search, TS）是局部搜索算法（Local Search, LS）的推广，其核心思想是通过一种“记忆”方法，即禁忌表，对已经进行的搜索过程进行记录和选择，指导下一步的搜索方向，进而避免陷入局部最优解。1986年，Glover首次提出了禁忌搜索算法概念。1997年，Glover与Laguna合著的《Tabu Search》专著出版，标志着禁忌搜索相关研究日趋完善。随后，禁忌搜索算法相关研究步入了高速发展阶段，被广泛应用到学术和工程的各个领域，并取得了较大成功。本文将介绍禁忌搜索算法的基本理论及算法流程，并以TSP问题为例进行编程求解。本文具体内容包括：

1、介绍禁忌搜索算法的基本理论；

2、介绍禁忌搜索算法的算法流程；

3、介绍基于TSP问题的禁忌搜索算法Python代码。

▲ 本文完整代码链接可在『交通与优化』公众号后台回复“TS”获取

禁忌搜索算法的原理介绍

禁忌搜索算法是一种全局逐步寻优的全局性邻域搜索算法，主要涉及以下10个关键技术：编码、初始解、评价函数、邻域结构、候选集、禁忌表、选择策略、破禁策略、停止准则以及解码。

1.1 编码

编码是将优化问题的决策变量从其解空间转换到算法所能处理的搜索空间的表示方法。常用的编码方式有两种：二进制编码和序列编码。以TSP问题为例对这两种编码方式进行说明。

（1）二进制编码

二进制编码方法是使用0，1两个基本字符将决策变量表示成二进制编码符号串。假设TSP问题中城市数量为N=4，旅行商的路线为从城市4出发，依次经过城市2、1、3，最终返回至城市4。此时的路线对应TSP问题的一个可行解（路线），即：城市4→城市2→城市1→城市3→城市4。采用二进制编码方式，可将其表示为：，如图1所示。

二进制编码方式符合计算机最小字符集原则，具有可靠性高，易于实现等优点。但也存在着可读性差等缺点。

（2）序列编码

序列编码是一种根据问题决策变量的形式，进行特殊编码的方法。以上述TSP问题为例进行说明，当旅行商路线为城市4→城市2→城市1→城市3→城市4时，采用序列编码，可将其表示为：{4, 2, 1, 3}，如图2所示。

序列编码方式具有可读性强、易于理解等优点。从图2可以看出，序列编码方式显得更为自然、合理。

1.2初始解

在寻求最优解的过程中，禁忌搜索算法是以初始解为基础进行全局搜索的。初始解可以随机产生，也可以根据经验或由其他算法产生。以TSP问题为例，对随机产生初始解进行说明：假设城市数量为4，产生初始解时，需要生成城市序列，每个序列中包含4个元素，每个元素的取值均为范围内的整数，且每个元素的取值各不相同。采用随机方法产生的初始解为{4,2,1,3}，其对应的路线为：城市4→城市2→城市1→城市3→城市4。

1.3评价函数

通常可以将目标函数作为评价函数。评价函数也可以有其他形式，可根据具体问题对评价函数进行设计。

1.4邻域结构

本节介绍TSP问题三种常见的邻域结构：元素交换，元素嵌入和元素反转。

（1）元素交换

元素交换指的是将部分元素的位置进行交换。以TSP问题为例进行说明：如图3所示，假设当前解为，选中元素{2}和{4}，交换其位置，产生新的解为：(1, 4, 3, 2)。

（2）元素嵌入

元素嵌入指的是将某部分元素序列嵌入到其他位置。以TSP问题为例进行说明：如图4所示，假设当前解为，选中序列{3, 4}，将其嵌入至元素{1}和{2}之间，产生新的解为：(1, 3, 4, 2)。

（3）元素反转

元素反转指的是当前解的序列进行反转。以TSP问题为例进行说明：如图5所示，假设当前解为，将当前解的所有元素进行反转，产生新的解为：(4, 3, 2, 1)。

1.5候选集

一般地，禁忌搜索算法不需要在每次迭代中搜索邻域中所有的解，可以只搜索邻域的一个子集，即候选集。通过这种方式，可以只对当前解邻域内部分解进行评价，有利于加快算法的搜索速度。例如：假设TSP问题中当前解为，其邻域为：，可选择其中

作为候选解集合。但是，由于没有搜索当前解的整个邻域，这种方法也可能使解的质量变差。

1.6禁忌表

禁忌表是用来存放禁忌对象的一个容器，禁忌表中的禁忌对象在解禁之前不能被再次搜索。引入禁忌表是为了防止搜索出现循环，陷入局部最优解的情况出现。禁忌表包括禁忌对象和禁忌长度两个主要指标。禁忌对象指的是禁忌表中被禁止变化的元素（即禁掉谁）。禁忌对象通常有三种：

（1）以状态本身或者状态的变化作为禁忌对象。例如：在上述TSP问题中，可将解作为禁忌对象。

（2）以状态分量以及分量的变化作为禁忌对象。例如：在上述TSP问题中，解的变化为：，可将元素{3, 4}变化的过程作为禁忌对象。

（3）以状态分量以及目标值变化作为禁忌对象。例如：在上述TSP问题中，将解(1, 2, 3, 4)对应的目标函数值44作为禁忌对象。

禁忌长度指的是搜索过程中，禁忌对象被禁忌的步长（即禁多久）。禁忌长度可以设置为固定值，也可以是动态变化的，可根据具体问题设计某种规则进行变化。例如：可将TSP问题中禁忌长度设置为固定值10。此外，亦可动态设置禁忌长度，例如：前100次迭代禁忌长度设置为5，后续迭代禁忌长度设置为10。

禁忌长度的选择会对算法的有效性产生重要影响。如果禁忌长度过短，可能会出现无法跳出局部最优解的现象。如果禁忌长度过长，可能会出现候选解全部被禁忌，搜索无法进行下去的现象。此外，禁忌长度过长，搜索时间通常也会更长。

1.7选择策略

在禁忌搜索算法中，选择的目的是从候选解集合中选择合适的解，以进行下一次迭代运算。一种常用的方法是选择目标函数值最小的解，从而使算法趋向于最优解。例如：在上述例子中，候选解集合

\{(4,2,3,1),(1,3,2,4),(1,2,4 ,3)\}

对应的目标函数值集合为，可选择目标函数值为所对应的解作为下一次迭代的当前解。

1.8破禁策略

破禁策略也叫藐视准则，是对禁忌表的适度放松。其作用是为了防止由于禁忌表的存在，而错过优异解的现象出现，同时也是为了避免候选集中所有的解均被禁忌掉，算法无法继续进行的情况出现。具体过程如下：从候选集中选择适应值最好的候选解，将其与当前最好解进行比较，若其是被禁忌的解，但是优于当前最好解，则可将其解禁，用来作为下一迭代的当前解。

1.9停止准则

与其他启发式算法类似，常用的停止准则有以下三种：

（1）设置最大迭代次数，例如：将算法的最大迭代次数设置为500。

（2）设置最好的解不再更新的次数，例如：将最好的解不再更新的次数设置为100，达到该条件后停止迭代。

（3）设置某个禁忌对象禁忌的次数，例如：假设上述TSP问题中，禁忌对象为解(1, 2, 3, 4)，当其禁忌次数达到50次时，停止迭代。

1.10解码

解码是编码的反变换过程。算法终止运行并输出最优结果后，可以进一步将最优结果解译为更为直观、易懂的表现形式。以TSP问题为例，若输出的最优解为(4, 2, 3, 1)，将其解码可得如下结果：“最优路径为：城市4→城市2→城市3→城市1→城市4”。

算法流程

禁忌搜索算法的主要步骤如下：

（1）初始化，产生初始解，清空禁忌表；

（2）判断是否满足终止条件，若满足，则输出当前最优解并进一步执行解码操作；否则执行步骤（3）；

（3）在当前解的邻域内产生候选解集合，并计算候选解的适应值；

（4）将适应值最好的解与当前搜索得到的最好的解进行比较：如果优于当前最好解，更新其为最好的解，并作为下次迭代的当前解；否则从候选集中选择未被禁忌的最好解作为下次迭代的当前解；

（5）更新禁忌表；

（6）返回步骤（2）。

禁忌搜索算法的流程如图6所示。

Python代码实现

本节以禁忌搜索算法求解50个城市TSP问题为背景，从环境声明、主函数以及子函数定义几个方面，详细介绍Python代码实现过程。

3.1环境声明

首先，介绍本次代码实战所使用的环境，包含关于数值运算的函数库numpy、关于随机数的模块random和读取Excel表格的模块xlrd，其中random库和copy库是内建函数，xlrd库的版本为1.1.0。

1. import numpy as np  #关于数值运算的函数库
2. import random  #关于随机数的函数库
3. import xlrd #读取excel的函数库
4. import copy  #复制变量的函数
5. import matplotlib.pyplot as plt #绘图
6. from matplotlib import rcParams #绘图

接下来，为增加代码的可读性，将城市相关信息与求解结果封装为类




    
1. ############## 类声明 ###############################
2. class Node:   #将节点的所有信息打包为Node类变量
3.     def __init__(self):
4.         self.node_id = 0 #节点编号
5.         self.x = 0.0     #节点横坐标
6.         self.y = 0.0     #节点纵坐标
7. class Solution: #将求解打包为Solution类变量
8.     def __init__(self):
9.         self.cost = 0
10.         self.route = []
11.     def copy(self):
12.         solution = Solution()
13.         solution.cost = self.cost
14.         solution.route = copy.deepcopy(self.route)
15.         return solution

为增强代码的可拓展性，将城市编号及坐标信息存储在Excel文件中，读取Excel文件中的数据并存储在字典变量中。

1. ############## 数据导入 #############################
2. #将"input_node.xlsx"中的节点信息汇总为集合N
3. N={}
4. book = xlrd.open_workbook("input_node.xlsx")
5. sh = book.sheet_by_index(0)
6. for l in range(1, sh.nrows):  # read each lines
7.     node_id = str(int(sh.cell_value(l, 0)))
8.     node = Node()
9.     node.node_id = int(sh.cell_value(l, 0))
10.     node.x = float(sh.cell_value(l, 1))
11.     node.y = float(sh.cell_value(l, 2))
12.     N[l-1]=node

基于以上数据，计算城市之间的欧式距离，并存储至二维列表中。在该步骤中，输入参数为各城市的编号及坐标信息，输出结果为城市之间的欧式距离（distance）。

1. #根据节点信息计算距离矩阵distance，distance[i][j]代表从i节点到j节点的出行成本
2. distance=[[0 for j in range(sh.nrows-1)] for i in range(sh.nrows-1)]
3. for i in range(sh.nrows-1):
4.     for j in range(sh.nrows-1):
5.         distance[i][j]=((N[i].x-N[j].x)**2+(N[i].y-N[j].y)**2)**0.5

算法关键参数设置如下，包括最大迭代次数、节点数量、禁忌长度和邻域候选集的个数。

1. ############## 参数设置 ############################
2. Iter=500      #迭代次数
3. City_number = sh.nrows-1  # 城市数量
4. tabu_length = City_number*0.2  #禁忌长度
5. candidate_length = 200 #候选集的个数

分别记录当前以及历次迭代的最优值和最优路径。具体的，用history_best记录历次迭代的最优值；history_route记录当前迭代的最优路径。

1. history_best = [] #存放历史最优目标值
2. history_route=[[]]#存放历史最优路径
3. best_solution.cost = 9999#初始最优值

3.2主函数

主函数部分定义禁忌搜索算法的整个流程，算法的禁忌对象为状态分量的变化，禁忌步长为City_number*0.2，其中City_number为城市的数量。算法的终止条件为达到最大迭代次数500。具体代码如下：

1. ############## 主函数 ############################
2. for k in range(Iter):
3.     candidate,candidate_distance,swap_position = get_candidate(current_solution.route,distance,City_number)
4.     min_index = np.argmin(candidate_distance)
5.     #判断邻域动作是否在禁忌表中，若不在禁忌表中则判断候选集中最小值与全局最优值的大小
6.     if swap_position[min_index] not in tabu_list:#若邻域动作不在禁忌表中
7.         if candidate_distance[min_index] #若候选集中最小值小于最优值，则更新当前解与最优解
8.             best_solution.cost = candidate_distance[min_index]
9.             best_solution.route = candidate[min_index].copy()
10.             current_solution.route = candidate[min_index].copy()
11.             current_solution.cost = candidate_distance[min_index]
12.             tabu_list.append(swap_position[min_index])#更新禁忌表
13.             history_best.append(best_solution.cost)
14.             history_route.append(best_solution.route)
15.         else


    
:#若候选集中最小值大于最优值，则只更新当前解
16.             current_solution.route = candidate[min_index].copy()
17.             current_solution.cost = candidate_distance[min_index]
18.             tabu_list.append(swap_position[min_index])
19.             history_best.append(best_solution.cost)
20.             history_route.append(best_solution.route)
21.     else:#若邻域动作在禁忌表中，判断候选集中最小值与全局最优值的关系
22.         if candidate_distance[min_index] #若候选集中最小值小于最优值，则此邻域动作解禁，并更新全局最最优解与当前解
23.             del_list.append(tabu_list[tabu_list.index(swap_position[min_index])])
24.             del tabu_list[tabu_list.index(swap_position[min_index])]
25.             best_solution.cost = candidate_distance[min_index]
26.             best_solution.route = candidate[min_index].copy()
27.             current_solution.route = candidate[min_index].copy()
28.             current_solution.cost = candidate_distance[min_index]
29.             history_best.append(best_solution.cost)
30.         else:#若候选集中最小值大于最优值，则将此解设置为最大，不再搜索，更新当前解
31.             candidate_distance[min_index] = 99999
32.             current_solution.route = candidate[min_index].copy()
33.             current_solution.cost = candidate_distance[min_index]
34.             tabu_list.append(swap_position[min_index])#更新禁忌表
35.             history_best.append(best_solution.cost)
36.             history_route.append(best_solution.route)
37.     if len(tabu_list) > tabu_length:#禁忌表移动
38.         del tabu_list[0]

3.3子函数

子函数部分包括3个关键技术。

（1）关键技术1：计算目标函数值。

1. ############## 计算路径成本消耗 ############################
2. def route_cost(route,distance,City_number):
3.     cost = 0
4.     for i in range(City_number - 1):
5.         cost += distance[int(route[i])][int(route[i + 1])]
6.     cost += distance[int(route[City_number - 1])][int(route[0])]  # 加上返回起点的距离
7.     return cost

（2）关键技术2：实现两点交换的邻域结构，即：2-opt。

1. ############## 邻域操作 ############################
2. def swap(route,route_i,route_j):
3.     new_route=copy.deepcopy(route)
4.     while route_i5.         temp=new_route[route_i]
6.         new_route[route_i]=new_route[route_j]
7.         new_route[route_j]=temp
8.         route_i+=1
9.         route_j-=1    
10.     return new_route

（3）关键技术3：根据邻域结构随机生成候选解集合。

1. ############## 产生邻域候选解集合############################
2. def get_candidate(route,distance,City_number):
3.     swap_position = []
4.     i = 0
5.     while i 6.         current = random.sample(range(0,City_number),2)
7.         if current not in swap_position:
8.             swap_position.append(current)
9.             candidate[i] = swap(route,current[0],current[1])
10.             candidate_distance[i] = route_cost(candidate[i],distance,City_number)
11.             i += 1
12.     return candidate,candidate_distance,swap_position