捕食搜索算法（PS）-求解TSP问题

捕食搜索算法（PS）

PS算法是一种用于解决组合优化问题的模拟动物模拟动物捕食行为的空间搜索策略。

算法背景简介

【1】在没有发现猎物和猎物痕迹时，在整个空间内沿一定的方向快速的寻找猎物；
【2】一旦发现猎物或者发现猎物痕迹，捕猎者立即改变自己的运动方式，减慢速度，在该附近区域进行集中搜索，持续不断的接近猎物。
【3】在搜寻一段时间没有找到猎物后，捕猎者将放弃这种集中的区域，而继续在整个空间寻找猎物，即所谓的区域限制的搜索策略。
！！！隐含应用：猎物在搜索空间内是聚集的，若为随机分布则算法可能失效

算法优势

区域限制搜索策略可以很好的权衡全面和集中搜索，实现搜索广度和深度的结合。
文献指出：该算法适合求解TSP,GB(二分图)和 Flow-shop Sequencing等组合优化问题。

核心要义

使用限制（Restriction）来表征较优解的临近域大小，注意！！！限制级别越高，越宽松，搜索范围更大。以此来实现搜索空间的增大和减小，平衡算法的探索和开发能力。
具体实现上：当局部搜索陷入局部最优点时，通过拉高限制级别（lhighthreshold），使得算法可以接受当前解向较差的解的移动，进而在较差的解附近再开始新一轮的搜索。（实践中发现这里得避免算法在这两个“相似”的解附近循环跳跃，陷入死循环）

具体解释一下如何用限制来调控搜索范围

对于极小化目标而言，限制【0】到限制【3】，目标函数值依次增大（即该级别的阈值），在当限制级别L较小时，比如为0，在空间中搜索可以接受的解的空间就在图中最下面包含1的那一块区域，当在这个限制级别（L=0）下搜索若干次无果后，令L=L+1(此时L=1)，可搜索的范围变成了包含{1，2，10，12}所在的更大的区域，然后依次增大限制级别（即限制级别越大，搜索区域越宽松）

如何跳出局部最优点？

当L增大到一个阈值Lthreshold时，意味着算法已经在所限制的区域内搜索了跟多很多次而找不到改进解了，所以将L设置为一个相对较高的值LhighThreshold，这个阈值表示，就算搜索邻域得到一个较差解的适应值，，算法还是允许当前解向该差解移动，然后从该差解附近开始搜索，以此跳出原来所限制的较小的区域。PS:在第三步中若限制级别L很大，就可以接受较差解，跳出局部最优

关于算法中参数的设定（具体看代码）

代码演示

# -*- coding: utf-8 -*-#import package
import copy
import math
import randomimport matplotlib as mpl
import matplotlib.pyplot as pltmpl.rcParams['font.sans-serif'] = ['KaiTi']
mpl.rcParams['font.serif'] = ['KaiTi']#import datasets
#城市坐标  31个城市
def tsp_dataset(filename):location = []with open(filename) as f:for _ in range(7):f.readline()while True:content = f.readline() if content.strip('\n') == 'EOF':breakline_list = content.strip('\n').split()location.append((float(line_list[1]),float(line_list[2])))return location#内部加载城市数据
# location = [(1304,2312),(3639,1315),(4177,2244),(3712,1399),
#             (3488,1535),(3326,1556),(3238,1229),(4196,1004),
#             (4312,790),(4386,570),(3007,1970),(2562,1756),
#             (2788,1491),(2381,1676),(1332,695),
#             (3715,1678),(3918,2179),(4061,2370),
#             (3780,2212),(3676,2578),(4029,2838),
#             (4263,2931),(3429,1908),(3507,2367),
#             (3394,2643),(3439,3201),(2935,3240),
#             (3140,3550),(2545,2357),(2778,2826),
#             (2370,2975)
# ]
#外部加载城市数据
location = tsp_dataset('berlin52.txt')#distance matrix
n = len(location)
distance = [[0 for col in range(n)] for raw in range(n)]
#example:
# distance = [[0,0,0],
#            [0,0,0],
#            [0,0,0]]
for i in range(n):for j in range(n):distance[i][j] = math.sqrt((location[i][0]-location[j][0])**2+(location[i][1]-location[j][1])**2)#function
def sequence_generate():#random initial solutiongen_seq = random.sample(range(1,n+1),n)#1,2,3,4,……,n !!!no 0return gen_seqdef cacl_cost(seq):#calculate the cost of the seqseq_cost = 0for j in range(0,len(seq)-1):cost = distance[seq[j]-1][seq[j+1]-1]#recode from 0 rowseq_cost += costcost2 = distance[seq[0]-1][seq[-1]-1]#come back to startCityseq_cost += cost2return seq_costdef adjacent_seq(seq):'''#create a new solution basing on current seq 2-opt#不能产生和seq一摸一样的解！！！问题出在列表的切片上'''while True:exchange_number = random.sample(range(0,n),2)exchange_number.sort()if exchange_number[0] == (exchange_number[1]-1):#两个位置相邻continueelse:exchange_seq = seq[int(exchange_number[0]):int(exchange_number[1])]exchange_seq.reverse()new_seq = copy.deepcopy(seq)new_seq[int(exchange_number[0]):int(exchange_number[1])] = exchange_seqbreakreturn new_seqdef evalution_domain(seq , s_number):'''create s_number adjacent_seq\ns_number:在当前解附近搜索次数，每次搜索会得到一个解'''evalution_seq = [[0 for rol in range(n)] for raw in range(s_number)]for i in range(0,s_number):evalution_seq[i] = adjacent_seq(seq)return evalution_seqdef restriction_level(seq):'''calculate restrictionLevel in  the algorithm\n限制等级：numlevle+1'''restrict_domin = evalution_domain(seq,numlevle)restrict_cost = [0 for i in range(numlevle)]for i in range(numlevle):restrict_cost[i] = cacl_cost(restrict_domin[i])restrict_level = copy.copy(restrict_cost)restrict_level.sort()restrict_level.insert(0,cacl_cost(seq))return restrict_level#paramas setting
evalution_number = int(0.5*(n**2-n))#5% of the neighborhoods of seq 局部搜索的大小
numlevle = n #restriction level usually set to 'n'
cthreshold = 3*n#3*n
lthreshold = round(n/20)#n/20
lhighthreshold = n-round(n/20)#n-[n/20]best_seq = sequence_generate()
best_cost = cacl_cost(best_seq)
print('random initial solution:')
print(best_seq)
print('initial solution cost:')
print(best_cost)l = 0
counter = 0
res_level = restriction_level(best_seq)
solution = best_seq
cost_trend = [best_cost]#PSA algorithm
while l < n:# if l == 30:#         print('test')if l == lthreshold:#陷入局部最优解了，放开搜索区域l = lhighthresholdwhile True:#在当前解附近搜索evalution_number次，找到最好的一个解ns_domain = evalution_domain(solution,evalution_number)ns_cost = [0 for i in range(evalution_number)]for ns in range(0,evalution_number):ns_cost[ns] = cacl_cost(ns_domain[ns])proposal_cost = min(ns_cost)proposal_index = ns_cost.index(proposal_cost)proposal = ns_domain[proposal_index]if proposal_cost < res_level[l]:solution = proposalif proposal_cost < best_cost:best_seq = solutionbest_cost = proposal_costres_level = restriction_level(best_seq)l = 0counter = 0cost_trend.append(best_cost)else:counter +=1else:counter +=1if counter > cthreshold:counter = 0l += 1break#result
print('PSA SEARCH solution:')
print(best_seq)
print('PSA SEARCH cost:')
print(best_cost)#绘制行驶路线图
best_seq.append(best_seq[0])
x = [location[index-1][0] for index in best_seq]
y = [location[index-1][1] for index in best_seq]
plt.figure(figsize=(15,10))
plt.plot(x,y,'o')
plt.plot(x,y,linewidth=1,color='red')
plt.plot(x[0],y[0],'v',markersize=20)
plt.title('PSA_TSP')
plt.show()plt.figure(figsize=(15,10))
plt.plot(range(len(cost_trend)),cost_trend,lineWidth=2)
plt.show()
'''
random initial solution:
[27, 43, 49, 28, 24, 10, 31, 38, 42, 47, 14, 36, 17, 2, 34, 3, 18, 1, 39, 7, 6, 35, 25, 50, 41, 48, 46, 5, 13, 22, 33, 30, 16, 23, 26, 8, 20, 15, 32, 40, 21, 51, 12, 37, 19, 4, 29, 11, 44, 9, 45]
initial solution cost:
29242.35004468245
PSA SEARCH solution:
[24, 11, 27, 25, 46, 13, 51, 12, 26, 10, 50, 32, 42, 3, 5, 14, 4, 23, 47, 37, 36, 39, 38, 33, 34, 35, 48, 31, 44, 18, 9, 8, 7, 40, 2, 17, 16, 20, 41, 6, 1, 29, 28, 49, 19, 22, 30, 21, 43, 15, 45]
PSA SEARCH cost:
8345.526874164043
'''

运行结果

random initial solution:
[27, 43, 49, 28, 24, 10, 31, 38, 42, 47, 14, 36, 17, 2, 34, 3, 18, 1, 39, 7, 6, 35, 25, 50, 41, 48, 46, 5, 13, 22, 33, 30, 16, 23, 26, 8, 20, 15, 32, 40, 21, 51, 12, 37, 19, 4, 29, 11, 44, 9, 45]
initial solution cost:
29242.35004468245
PSA SEARCH solution:
[24, 11, 27, 25, 46, 13, 51, 12, 26, 10, 50, 32, 42, 3, 5, 14, 4, 23, 47, 37, 36, 39, 38, 33, 34, 35, 48, 31, 44, 18, 9, 8, 7, 40, 2, 17, 16, 20, 41, 6, 1, 29, 28, 49, 19, 22, 30, 21, 43, 15, 45]
PSA SEARCH cost:
8345.526874164043

参考资料：
【1】汪定伟, 王俊伟, 王洪峰, 等. 智能优化方法[M]. Gao deng jiao yu chu ban she, 2007.
【2】论文算法复现—智能调度(1) | 捕食搜索算法
【3】实验数据集-berlin52.tsp.gz