禁忌搜索(Tabu Search)算法及python实现

版权声明：本文为博主原创文章，博客地址：/adkjb，未经博主允许。

禁忌搜索(Tabu Search,TS,又称禁忌搜寻法）是一种现代启发式算法，由美国科罗拉多大学教授Fred Glover在1986年左右提出的，是一个用来跳脱局部最优解的搜索方法。其先创立一个初始化的方案；基于此，算法“移动”到一相邻的方案。经过许多连续的移动过程，提高解的质量。 [ 百度百科 ]

由于网上提供的代码大多是C和Java的，这里笔者使用Python实现该算法

【目录】

-TS算法原理详解

-举例详述TS算法过程

-Python实现TS算法

-总结

TS算法原理详解

邻域

对于组合优化问题，给定任意可行解x，x∈D，D是决策变量的定义域，对于D上的一个映射：N：x∈D→N(x)∈2(D) 其中2(D)表示D的所有子集组成的集合，N(x)成为x的一个邻域，y∈N(x)称为x的一个邻居。

候选集合

候选集合一般由邻域中的邻居组成，可以将某解的所有邻居作为候选集合，也可以通过最优提取，也可以随机提取，例如某一问题的初始解是[1,2,3]，若通过两两交换法则生成候选集合，则可以是[1,3,2],[2,1,3],[3,2,1]中的一个或几个。

禁忌表

禁忌表包括禁忌对象和禁忌长度。由于在每次对当前解的搜索中，需要避免一些重复的步骤，因此将某些元素放入禁忌表中，这些元素在下次搜索时将不会被考虑，这些被禁止搜索的元素就是禁忌对象；

禁忌长度则是禁忌表所能接受的最多禁忌对象的数量，若设置的太多则可能会造成耗时较长或者算法停止，若太少则会造成重复搜索。

评价函数

用来评价当前解的好坏，TSP问题中是总旅程距离。

特赦规则

禁忌搜索算法中，迭代的某一步会出现候选集的某一个元素被禁止搜索，但是若解禁该元素，则会使评价函数有所改善，因此我们需要设置一个特赦规则，当满足该条件时该元素从禁忌表中跳出。

终止规则

一般当两次迭代得到的局部最优解不再变化，或者两次最优解的评价函数差别不大，或者迭代n次之后停止迭代，通常选择第三种方法。

举例详述TS算法过程

现有一架飞机，从A点出发，需要经过B,C,D,E,F之后返回A点，且每个点只能经过一次，最后返回A点，求最短路径。

该问题是一个Hamilton回路问题，其中起点和终点已经固定，因此我们可以将解形式记为，例如【A,D,C,F,E,A】,每次只需变换中间两个元素即可，现在我们将禁忌长度设置为2，候选集合长度定义为4，迭代次数为100，通过以下步骤能使读者更清洗的了解TS算法的步骤。

给定任意初始解 x1=【A,D,C,F,E,A】f(x1)=10，历史最优为10

我们发现对x1交换D和E时，f最优，此时x2=【A,E,C,F,D,A】 f(x2)=6，历史最优为6，将D-E放入禁忌表中

我们发现对x2交换C和D时，f最优，此时x3=【A,E,D,F,C,A】 f(x3)=5，历史最优为5，将D-C放入禁忌表中

此时我们发现对x3交换D和C时最优，但是由于D-C已经在禁忌表中，因此我们退而求其次，对x3交换F和D，此时x4=【A,E,F,D,C,A】 f(x4)=10，历史最优为5，将F-D放入禁忌表中，由于禁忌长度为2，因此将最先放入禁忌表中的D-E移除禁忌表

此时我们发现对x4交换D和C时最优，虽然D-C已经在禁忌表中，但是f(D-C)<历史最优5，因此满足特赦规则，现在将D-C移除禁忌表，此时x5=【A,E,F,C,D,A】 f(x5)=4，历史最优为4，然后再将D-C放入禁忌表

依次迭代下去，当迭代次数超过100时停止迭代，历史最优值即为输出解

Python实现TS算法

此处以一个简单版的TSP问题为例

现有29个城市，提供了经纬度，一架飞机需要从1号城市出发，途径剩下的28个城市然后返回1号城市，求最短距离

数据如下

解的形式为例如[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]，这去除了开头结尾的0，即1号城市，只用计算中间的航程即可

import randomimport datetimecity = []eg=for line in open("tsp.txt"):place,lon,lat = line.strip().split(" ")city.extend([(place,(lon,lat))]) #导入城市的坐标def printtravel(vec):print(city[0],city[vec[0]])for i in range(len(vec)-1):print(city[vec[i]],city[vec[i+1]])print(city[vec[i+1]],city[0]) #打印结果函数 eg=[i for i in range(1,29)] #一个例子printtravel(eg)

【下面是打印出来的路径，即1→2→3→4→…→29→1】

(‘1’, (‘1150.0’, ‘1760.0’)) (‘2’, (‘630.0’, ‘1660.0’))

(‘2’, (‘630.0’, ‘1660.0’)) (‘3’, (‘40.0’, ‘2090.0’))

(‘3’, (‘40.0’, ‘2090.0’)) (‘4’, (‘750.0’, ‘1100.0’))

(‘4’, (‘750.0’, ‘1100.0’)) (‘5’, (‘750.0’, ‘2030.0’))

(‘5’, (‘750.0’, ‘2030.0’)) (‘6’, (‘1030.0’, ‘2070.0’))

(‘6’, (‘1030.0’, ‘2070.0’)) (‘7’, (‘1650.0’, ‘650.0’))

(‘7’, (‘1650.0’, ‘650.0’)) (‘8’, (‘1490.0’, ‘1630.0’))

(‘8’, (‘1490.0’, ‘1630.0’)) (‘9’, (‘790.0’, ‘2260.0’))

(‘9’, (‘790.0’, ‘2260.0’)) (‘10’, (‘710.0’, ‘1310.0’))

(‘10’, (‘710.0’, ‘1310.0’)) (‘11’, (‘840.0’, ‘550.0’))

(‘11’, (‘840.0’, ‘550.0’)) (‘12’, (‘1170.0’, ‘2300.0’))

(‘12’, (‘1170.0’, ‘2300.0’)) (‘13’, (‘970.0’, ‘1340.0’))

(‘13’, (‘970.0’, ‘1340.0’)) (‘14’, (‘510.0’, ‘700.0’))

(‘14’, (‘510.0’, ‘700.0’)) (‘15’, (‘750.0’, ‘900.0’))

(‘15’, (‘750.0’, ‘900.0’)) (‘16’, (‘1280.0’, ‘1200.0’))

(‘16’, (‘1280.0’, ‘1200.0’)) (‘17’, (‘230.0’, ‘590.0’))

(‘17’, (‘230.0’, ‘590.0’)) (‘18’, (‘460.0’, ‘860.0’))

(‘18’, (‘460.0’, ‘860.0’)) (‘19’, (‘1040.0’, ‘950.0’))

(‘19’, (‘1040.0’, ‘950.0’)) (‘20’, (‘590.0’, ‘1390.0’))

(‘20’, (‘590.0’, ‘1390.0’)) (‘21’, (‘830.0’, ‘1770.0’))

(‘21’, (‘830.0’, ‘1770.0’)) (‘22’, (‘490.0’, ‘500.0’))

(‘22’, (‘490.0’, ‘500.0’)) (‘23’, (‘1840.0’, ‘1240.0’))

(‘23’, (‘1840.0’, ‘1240.0’)) (‘24’, (‘1260.0’, ‘1500.0’))

(‘24’, (‘1260.0’, ‘1500.0’)) (‘25’, (‘1280.0’, ‘790.0’))

(‘25’, (‘1280.0’, ‘790.0’)) (‘26’, (‘490.0’, ‘2130.0’))

(‘26’, (‘490.0’, ‘2130.0’)) (‘27’, (‘1460.0’, ‘1420.0’))

(‘27’, (‘1460.0’, ‘1420.0’)) (‘28’, (‘1260.0’, ‘1910.0’))

(‘28’, (‘1260.0’, ‘1910.0’)) (‘29’, (‘360.0’, ‘1980.0’))

(‘29’, (‘360.0’, ‘1980.0’)) (‘1’, (‘1150.0’, ‘1760.0’))

def costroad(road):cost = ((float(city[0][1][0])-float(city[road[0]][1][0]))**2+(float(city[0][1][1])-float(city[road[0]][1][1]))**2)**0.5for i in range(len(road)-1):cost = cost+((float(city[road[i+1]][1][0])-float(city[road[i]][1][0]))**2+(float(city[road[i+1]][1][1])-float(city[road[i]][1][1]))**2)**0.5cost=cost+((float(city[road[-1]][1][0])-float(city[0][1][0]))**2+(float(city[-1][1][1])-float(city[0][1][1]))**2)**0.5return(cost) #计算所求解的距离，这里为了简单，视作二位平面上的点，使用了欧式距离costroad(eg) #计算上例中的总路程

【结果】

25814.877362907795

def tabusearch(diedaitimes,cacu_time,tabu_length,origin_times,costf,printf):s1=datetime.datetime.now() #获取运行前的时间print("The program now is executing...")def pan_move(move_step,tabu_move): #判断移动是否在禁忌区域中，如果是返回True和该点索引，否则返回False和0if move_step in tabu_move:index = tabu_move.index(move_step)return(True,index)else:return(False,0)def pan_cost(cost,tabu_cost,t): #判断该移动是否比禁忌区域中该移动小，如果小则返回True，否则返回Falseif cost<tabu_cost[t]:return(True)else:return(False) def add_tabu(cost,move,tabu_cost,tabu_move,t): #为禁忌区域添加移动和成本，若超过T则剔除最先进入的禁忌tabu_cost.append(cost)tabu_move.append(move)if len(tabu_cost)>t:del tabu_cost[0]if len(tabu_move)>t:del tabu_move[0]return(tabu_cost,tabu_move)def cacu(vec,t): #为每一个初始解计算t次vec_set = []m_set = []cost_set = []h = []for i in range(t):v,m,c,h = move(vec,h)vec_set.append(v)m_set.append(m)cost_set.append(c)return(vec_set,m_set,cost_set)def cacu_tiqu(v1,m1,c1): #从上述t次筛选最小的解向量，移动和成本t = c1.index(min(c1))v_max = v1[t]m_max = m1[t]c_max = c1[t]return(v_max,m_max,c_max)def move(vec,h): #输出移动后的向量，和成本i = 1while i==1:m = random.sample(vec,2)m.sort()if m not in h:h.append(m)vec_copy = vec[:]vec_copy[vec_copy.index(m[0])] = m[1]vec_copy[vec_copy.index(m[1])] = m[0]cost = costf(vec_copy)i = 0return(vec_copy,m,cost,h) finall_road = []finall_cost = []for t1 in range(origin_times):road = [i for i in range(1,29)]random.shuffle(road)tabu_cost = []tabu_move = []for t in range(diedaitimes):i = 0while i==0:v1,m1,c1 = cacu(road,cacu_time)v_m,m_m,c_m = cacu_tiqu(v1,m1,c1)key1 = pan_move(m_m,tabu_move)if key1[0]:if pan_cost(c_m,tabu_cost,key1[1]):road = v_mfinall_road.append(road)finall_cost.append(c_m)tabu_cost,tabu_move = add_tabu(c_m,m_m,tabu_cost,tabu_move,tabu_length)i=1else:v1.remove(v_m)m1.remove(m_m)c1.remove(c_m)if len(v1)==0:i = 1else:tabu_cost,tabu_move = add_tabu(c_m,m_m,tabu_cost,tabu_move,tabu_length)road = v_mfinall_road.append(road)finall_cost.append(c_m)i = 1index = finall_cost.index(min(finall_cost))s2 = datetime.datetime.now()print("Successfully execute!,the program has executed for "+str((s2-s1).seconds)+" seconds!")return(finall_road[index],min(finall_cost),printf(finall_road[index]))

tabusearch(diedaitimes=100,cacu_time=100,tabu_length=10,origin_times=100,costf=costroad,printf=printtravel)

其中diedaitimes为每一个初始解的迭代次数，cacu_time为候选集合长度，tabu_length为禁忌长度，origin_times为整个程序循环次数，可以理解为使用不同个初始解，costf为成本函数，printtravel为打印结果函数，运行上一条程序可得结果：

The program now is executing…

Successfully execute!,the program has executed for 71 seconds!

(‘1’, (‘1150.0’, ‘1760.0’)) (‘6’, (‘1030.0’, ‘2070.0’))

(‘6’, (‘1030.0’, ‘2070.0’)) (‘13’, (‘970.0’, ‘1340.0’))

(‘13’, (‘970.0’, ‘1340.0’)) (‘15’, (‘750.0’, ‘900.0’))

(‘15’, (‘750.0’, ‘900.0’)) (‘4’, (‘750.0’, ‘1100.0’))

(‘4’, (‘750.0’, ‘1100.0’)) (‘10’, (‘710.0’, ‘1310.0’))

(‘10’, (‘710.0’, ‘1310.0’)) (‘2’, (‘630.0’, ‘1660.0’))

(‘2’, (‘630.0’, ‘1660.0’)) (‘26’, (‘490.0’, ‘2130.0’))

(‘26’, (‘490.0’, ‘2130.0’)) (‘3’, (‘40.0’, ‘2090.0’))

(‘3’, (‘40.0’, ‘2090.0’)) (‘29’, (‘360.0’, ‘1980.0’))

(‘29’, (‘360.0’, ‘1980.0’)) (‘5’, (‘750.0’, ‘2030.0’))

(‘5’, (‘750.0’, ‘2030.0’)) (‘9’, (‘790.0’, ‘2260.0’))

(‘9’, (‘790.0’, ‘2260.0’)) (‘12’, (‘1170.0’, ‘2300.0’))

(‘12’, (‘1170.0’, ‘2300.0’)) (‘21’, (‘830.0’, ‘1770.0’))

(‘21’, (‘830.0’, ‘1770.0’)) (‘20’, (‘590.0’, ‘1390.0’))

(‘20’, (‘590.0’, ‘1390.0’)) (‘18’, (‘460.0’, ‘860.0’))

(‘18’, (‘460.0’, ‘860.0’)) (‘14’, (‘510.0’, ‘700.0’))

(‘14’, (‘510.0’, ‘700.0’)) (‘17’, (‘230.0’, ‘590.0’))

(‘17’, (‘230.0’, ‘590.0’)) (‘22’, (‘490.0’, ‘500.0’))

(‘22’, (‘490.0’, ‘500.0’)) (‘11’, (‘840.0’, ‘550.0’))

(‘11’, (‘840.0’, ‘550.0’)) (‘7’, (‘1650.0’, ‘650.0’))

(‘7’, (‘1650.0’, ‘650.0’)) (‘25’, (‘1280.0’, ‘790.0’))

(‘25’, (‘1280.0’, ‘790.0’)) (‘19’, (‘1040.0’, ‘950.0’))

(‘19’, (‘1040.0’, ‘950.0’)) (‘16’, (‘1280.0’, ‘1200.0’))

(‘16’, (‘1280.0’, ‘1200.0’)) (‘23’, (‘1840.0’, ‘1240.0’))

(‘23’, (‘1840.0’, ‘1240.0’)) (‘8’, (‘1490.0’, ‘1630.0’))

(‘8’, (‘1490.0’, ‘1630.0’)) (‘27’, (‘1460.0’, ‘1420.0’))

(‘27’, (‘1460.0’, ‘1420.0’)) (‘24’, (‘1260.0’, ‘1500.0’))

(‘24’, (‘1260.0’, ‘1500.0’)) (‘28’, (‘1260.0’, ‘1910.0’))

(‘28’, (‘1260.0’, ‘1910.0’)) (‘1’, (‘1150.0’, ‘1760.0’))

([5,12,14,3,9,1,25,2,28,4,8,11,20,19,17,13,16,21,10,6,24,18,15,22,7,26,23,27],

11356.612831196624,None)

所得到的路径如下图，总航程约为11356

下面是路径，黑色为起点（这个是用R绘制的，就不上代码了）

总结

从本例出发，28个城市进行排列，共有28！≈3*10^29种组合，笔者使用暴力枚举测试了一百万个组合，耗费76秒得到里程为16270的一条航程，相比本文讨论的禁忌搜索算法，仅使用71秒就得到了11356的一条行程，因此具有非常好的效果，虽然从图像来看好像不是最优路径，但是通过调整参数会逐渐接近最优解。

数据集链接：/s/1DOUCHkgX_8KoxUtrmV-aTw