700字范文 > 局部搜索：禁忌搜索(Tabu Search)解决TSP问题的python案例

局部搜索：禁忌搜索(Tabu Search)解决TSP问题的python案例

时间：2019-08-18 23:30:07

禁忌搜索解决TSP问题的python案例

Tabu Search

禁忌（Tabu Search）算法是一种亚启发式(meta-heuristic)随机搜索算法，它从一个初始可行解出发，选择一系列的特定搜索方向（移动）作为试探，选择实现让特定的目标函数值变化最多的移动。为了避免陷入局部最优解，TS搜索中采用了一种灵活的“记忆”技术，对已经进行的优化过程进行记录和选择，指导下一步的搜索方向，这就是Tabu表的建立。【百度百科】

局部搜索的算法研究到这里，可以总结出几个重要的点：

1.问题的表示方法（Representation of the solution）

2.评价函数（Evaluation function）

3.如何产生新的邻居解（Neighbourhood function/operator）

4.接受邻居解的规则（Acceptance criteria）

在TSP问题中，产生邻居解的方法一般是随机交换任意的两个city。有关禁忌算法的介绍，百度上搜到的讲的清楚地并不多，这里列出几篇，有助于理解，个人觉得第二篇的代码写的不是很清楚：

浅谈禁忌搜索（TabuSearch）
/yang2110862/article/details/106888667

禁忌搜索(Tabu Search)算法及python实现
/adkjb/article/details/81712969

算法中最重要的就是两个概念：禁忌列表和特赦规则，禁忌列表用来记录尝试过得解法，可以想象为有大小限制的队列，队列满了之后就对很久以前的解法就pop掉。特赦规则可以忽略禁忌表的限制而接受该解。有些算法没有特赦规则，仅从邻居解中选择不在禁忌表中的最优解。

在维基百科的介绍中，给出了不包含特赦规则的算法伪代码（对于邻居解没有排序，而是遍历找出满足条件的最优解）：

/wiki/Tabu_search

在这里安利一个Python学习案例的GitHub，里面包含了各个领域的Python基础代码，机器学习，数据结构，计算机视觉，搜索算法等，本文参考了其中的部分代码：

/TheAlgorithms/Python

结合以上内容，开始我们的正题。

TSP问题的数据来源

TSP问题是什么，就不多介绍了，关注该问题的缘由是一个YouTube视频，有关傅里叶级数的，先放个图：

/watch?v=-qgreAUpPwM

使用傅里叶级数拟合图像边缘，借助Opencv我们可以轻松的得到图像的轮廓，而且是按照顺序的轮廓，但是对于多个轮廓，怎么把它们连在一起呢？（连在一起才能仅用一组傅里叶级数把它们拟合出来）这个问题类似一笔画问题，所以想到了TSP，但是用这个算法解决这个问题貌似并不恰当（我们需要的是需要把这几个独立的轮廓连成一个周期的解），所以这里仅仅是给TSP问题提供个数据，下面的演示中我们仅仅用字母n的轮廓寻找一条近似最优解，注意坐标的转换，图片的原点在左上角。（傅里叶级数的这个视频，最终用Java和Processing库复现了，但是没有解决一笔画问题，那就多组轮廓一起画…，以后再考虑别的算法吧。）

获得图片轮廓点集，有关findContours的介绍请自行查找，在这里我们得到了所有轮廓的点集，注意下采样，否则点太多了，不适合我们的TSP问题的数据。

#得到图像轮廓,每个轮廓下采样,多个轮廓合并在一起def get_edges(img_path):img = cv2.imread(img_path)width, height = img.shape[1], img.shape[0]gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)ret, binary = cv2.threshold(gray, 127, 255, cv2.THRESH_BINARY)contours, hierarchy = cv2.findContours(binary,cv2.RETR_LIST,cv2.cv2.CHAIN_APPROX_NONE) # img = cv2.drawContours(img,contours,-1,(255,0,0),3)# cv2.imshow('contours', img)# cv2.waitKey(0)# cv2.destroyAllWindows()edges = []for contour in contours[2:3]:count = 0for point in contour:if count == 3:edges.append([point[0][0] ,height - point[0][1]])count = 0else:count += 1return edges

算法内容

#得到距离矩阵,对称的,对角线为1,使用欧氏距离def get_weights(citys):#[[x,y],[x,y],[x,y],....,[x,y]]def cost(city_1, city_2):return pow(pow(city_1[0] - city_2[0], 2) + pow(city_1[1] - city_2[1], 2), 0.5)weights = np.ones((len(citys), len(citys)))for i, city_from in enumerate(citys):for j, city_to in enumerate(citys):if i > j:weights[i][j] = cost(city_from, city_to)weights[j][i] = weights[i][j]return weights#随机交换两个城市, 最多产生max_num个解, 返回path和cost,eg.[[0,2,5,9,....,0],1000]def find_neighbour(path, weights, max_num):solution_neighbours = []for i in range(0, max_num):exchange = random.sample(range(1, len(path)-1), 2)temp_path = copy.deepcopy(path)temp_path[exchange[0]] = path[exchange[1]]temp_path[exchange[1]] = path[exchange[0]]if temp_path not in solution_neighbours:cost = fitness(temp_path, weights)solution_neighbours.append([temp_path, cost])return solution_neighbours#最小优化目标,计算总距离(cost), path: [2,4,5,7,1,10,7,...,2]def fitness(path, weights):cost_sum = 0for i in range(0, len(path)-1):cost_sum += weights[path[i]][path[i+1]]return cost_sum

为了避免在寻找最优解的过程中频繁计算欧式距离，所以“一劳永逸”，先生成代价矩阵，其实是个对称的矩阵，对角线用不到，随便什么数值，生成矩阵的方法，使用numpy.ones((w, h))；

邻居解就是通过随机交换path中任意两个city的位置，注意生成新的路径时，用copy.deepcopy，有关python赋值和引用的内容，请自行脑补。在这个函数中，最多生成max_num个邻居解（邻居解的优化问题也是个研究方向）；

评价函数就是计算路径的cost，没什么可以解释的。

搜索过程

def tabu_search(edge_points, weights, iter, tabu_size, neighbour_num):#绘图部分pyplot.close()fig = pyplot.figure()path_fig = fig.add_subplot(1, 2, 1)cost_fig = fig.add_subplot(1, 2, 2)path_fig.axis("equal")path_fig.set_title('Best Path')cost_fig.set_title('Current Cost')cost_fig.set_xlabel('iterations')cost_fig.set_ylabel('fitness')pyplot.subplots_adjust(wspace= 0.5)pyplot.ion()#打开交互path_fig.scatter([i[0] for i in edge_points], [i[1] for i in edge_points], s=2 ,color='red')pyplot.pause(0.001)first_path = [i for i in range(0, len(weights))]first_path.append(0)curr_solution = [first_path, fitness(first_path,weights)]best_solution = copy.deepcopy(curr_solution)tabu_list = list()cost_history = list()cost_history.append(curr_solution[1])for k in range(0, iter):#更新绘图cost_fig.plot(cost_history, color='b')path_fig.plot([edge_points[p][0] for p in best_solution[0]], [edge_points[p][1] for p in best_solution[0]], color='b', linewidth=1)pyplot.pause(0.001)path_fig.lines.pop(0)neighbour_solution = find_neighbour(curr_solution[0], weights, neighbour_num)neighbour_solution.sort(key= lambda x: x[1])best_neighbour_solution_index = 0best_neighbour_solution = neighbour_solution[best_neighbour_solution_index]found = Falsefind_count = 0while (not found) and find_count < len(neighbour_solution):tabucode = [0, 0]first_city, second_city = 0, 0for index in range(0, len(curr_solution[0])):if curr_solution[0][index] != best_neighbour_solution[0][index]:tabucode[0], tabucode[1] = curr_solution[0][index], best_neighbour_solution[0][index]break# print(tabucode)# print(tabu_list)# if tabucode not in tabu_list and tabucode.reverse() not in tabu_list:if tabucode not in tabu_list:found = True#print("接受解%d"%(find_count))#没有在禁忌表中,接受该解,可能是次优的tabu_list.append(tabucode)curr_solution = best_neighbour_solutionif best_neighbour_solution[1] < best_solution[1]:best_solution = best_neighbour_solutionelse:#特赦规则,该邻居解是目前最好的解if best_neighbour_solution[1] < best_solution[1]:curr_solution = best_solutionbest_solution = best_neighbour_solutionfound = True#print("接受特赦解")#没有最优,选择次优else:best_neighbour_solution_index += 1best_neighbour_solution = neighbour_solution[best_neighbour_solution_index]find_count += 1if len(tabu_list)>tabu_size:tabu_list.pop(0)#print("禁忌表满")cost_history.append(curr_solution[1])print("iterations: %d, current cost:%.2f, best cost:%.2f"%(k, curr_solution[1], best_solution[1]))# os.system('pause')path_fig.plot([edge_points[p][0] for p in best_solution[0]], [edge_points[p][1] for p in best_solution[0]], color='b', linewidth=1)pyplot.savefig('result.jpg')return curr_solution, best_solution

解释再多，不如看代码，还是解释一下吧，两个主要变量：curr_solution，best_solution，分别代表当前解和历史最优解，当前解是产生邻居解参考基础，历史最优解是特赦的条件。邻居解先进行升序排列，寻找可接受解的过程可以概括为：优于历史最优解，不管有没有在禁忌表，一定接受，如果不优于历史最优解，接受不在禁忌表中的次优解，这就是内层循环的工作，找到接受解后停止。（脑残的我把best_neighbour_solution_index += 1里面的+=写错为=+，导致不能遍历…）

外层循环就是迭代次数了，对禁忌表的管理，更新绘图，打印过程，搜索结束后返回curr_solution，best_solution。

代码的剩余部分，注意把生成的轮廓随机打乱顺序。

parser = argparse.ArgumentParser()parser.add_argument("-i", "--iters", type=int, default=800, help="How many iterations the algorithm should perform")parser.add_argument("-s", "--size", type=int, default=40, help="Size of the tabu list")parser.add_argument("-n", "--num", type=int, default=200, help="Max numbers of neighbour solutions")opt = parser.parse_args()edge = get_edges('./CSDN.jpg')random.shuffle(edge)weights = get_weights(edge)print("tabu search start...\n\solution space size: %.2e, city numbers: %d\n\iterations: %d, tabulist size: %d, max generated solutions: %d"\%(factorial(len(edge)), len(edge), opt.iters, opt.size, opt.num))os.system('pause')tabu_search(edge, weights, opt.iters, opt.size, opt.num)

动态绘图

图片动态使用了pyplot.ion()

pyplot.close()fig = pyplot.figure()path_fig = fig.add_subplot(1, 2, 1)cost_fig = fig.add_subplot(1, 2, 2)path_fig.axis("equal")path_fig.set_title('Best Path')cost_fig.set_title('Current Cost')cost_fig.set_xlabel('iterations')cost_fig.set_ylabel('fitness')pyplot.subplots_adjust(wspace= 0.5)pyplot.ion()#打开交互path_fig.scatter([i[0] for i in edge_points], [i[1] for i in edge_points], s=2 ,color='red')pyplot.pause(0.001)######################更新绘图cost_fig.plot(cost_history, color='b')path_fig.plot([edge_points[p][0] for p in best_solution[0]], [edge_points[p][1] for p in best_solution[0]], color='b', linewidth=1)pyplot.pause(0.001)path_fig.lines.pop(0) #清楚上一次画的路径

两个子图，注意的地方就是pyplot.ion()，pyplot.pause(0.001)和path_fig.lines.pop(0)，最终结果：