元启发式算法 | 禁忌搜索(Tabu Search)解决TSP问题(Python)

[元启发式算法]禁忌搜索(Tabu Search)解决TSP问题(Python)

完整代码见github。

1.Tabu Search基本概念

禁忌搜索(Tabu Search,TS，以下简称TS） 是一种基于邻域搜索策略的元启发式算法，由Fred W. Glover 在1986提出[1]，并于1989构建[2][3]。应用于各类组合优化问题。

禁忌搜索的核心思想： 从一个初始解出发，按一系列规则对邻域进行探索，对已搜索过的途径和局部最优解进行记录（建立禁忌表，tabu list），并在进一步的迭代搜索中尽量避开这些对象（不是绝对禁止循环），减少重复搜索，从而保证对不同的有效搜索途径的探索，提高搜索效率。

禁忌搜索涉及到的几个概念。邻域（neighborhood）、禁忌表（tabu list）、禁忌长度（tabu length）、候选解（candidate）、藐视准则（aspiration criterion）(有的也叫特赦规则)

2. Tabu Search算法实现细节

以TSP（Traveling Salesman Problem）问题为例，0号为起始点，需要依次拜访1,2,…10号，最终返回0号，求最短路径。已知各个点坐标位置，各个点的位置如下图所示。解的形式表示为一个列表，如：route=[1,2,3,4,5,6,7,8,9,10].

（这里用是随机生成一系列坐标点。）

邻域（neighborhood）：邻域一般定义为由给定转化规则对给定的问题域上每结点进行转化所得到的问题域上结点的集合。那么在这里即1-10号点的排列组合。
候选解（candidate）：从邻域中选择若干个目标值或评价值最佳的邻居作为候选解。候选解集合的生成规则一定程度上决定了搜索的方向和邻域大小，十分关键。候选集过大增加计算内存和计算时间，过小容易过早陷入局部最优。候选集的选择一般由邻域中的邻居组成，可以选择所有邻居，也可以选择表现较好的邻居，还可以随机选择几个邻居。这里我们采用两两交换法则生成候选解集合。
禁忌表（tabu list）：记录最优候选解和对应元素，这些元素在下次搜索时将不会被考虑。
禁忌长度（tabu length）：禁忌表的最大长度。
藐视准则（aspiration criterion）：禁忌搜索算法中，迭代的某一步会出现候选集的某一个元素被禁止搜索，但是若解禁该元素，则会使评价函数有所改善，因此我们需要设置一个特赦规则，当满足该条件时该元素从禁忌表中跳出。
评价函数（evaluation）:评价解的好坏。在这里即路径距离。

禁忌搜索算法核心问题：
禁忌对象：禁掉谁？根据受禁对象的不同选择，可行解是一禁禁一个；还是一禁禁掉一大片。主要对禁忌范围，及搜索范围有影响
禁忌长度：禁多久？禁掉的东西什么时候放出来？禁忌长度过短，会导致循环；禁忌长度过长，会导致计算时间过长。
候选集： 邻域中可行解的选取？候选集的大小，过大增加计算内存和计算时间，过小过早陷入局部最优。

python代码实现，可直接运行。

from itertools import combinations
import os,sys,copy
import numpy as np
import time
import matplotlib.pyplot as plt
from GetData import *

class Tabu():
    def __init__(self,disMatrix,max_iters=50,maxTabuSize=10):
        """parameters definition"""
        self.disMatrix = disMatrix
        self.maxTabuSize = maxTabuSize
        self.max_iters = max_iters
        self.tabu_list=[]

    def get_route_distance(self,route):
        '''
        Description: function to calculate total distance of a route. evaluate function.
        parameters: route : list
        return : total distance : folat
        '''        
        routes = [0] + route + [0]    # add the start and end point 
        total_distance = 0
        for i,n in enumerate(routes):
            if i != 0 :
                total_distance = total_distance +  self.disMatrix[last_pos][n] 
            last_pos = n
        return total_distance

    def exchange(self,s1,s2,arr):
        """
        function to Swap positions of two elements in an arr
        Args: int,int,list
            s1 : target 1 
            s2 : target 2  
            arr : target array 
        Ouput: list
            current_list : target array 
        """
        current_list = copy.deepcopy(arr)
        index1 , index2 = current_list.index(s1) , current_list.index(s2)  # get index
        current_list[index1], current_list[index2]= arr[index2] , arr[index1]
        return current_list

    def generate_initial_solution(self,num=10,mode='greedy'):
        """
        function to get the initial solution,there two different way to generate route_init.
        Args: 
            num :  int
                the number of points 
            mode : string
                "greedy" : advance step by choosing optimal one 
                "random" : randomly generate a series number
        Ouput: list
            s_init : initial solution route_init
        """
        if mode == 'greedy':
            route_init=[0]
            for i in range(num):
                best_distance = 10000000
                for j in range(num+1):
                    if self.disMatrix[i][j] < best_distance and j not in route_init:  
                        best_distance = self.disMatrix[i][j]
                        best_candidate = j
                route_init.append(best_candidate)
            route_init.remove(0)
                            
        if mode == 'random':
            route_init = np.arange(1,num+1)  #init solution from 1 to num
            np.random.shuffle(route_init)  #shuffle the list randomly

        return list(route_init)

    def tabu_search(self,s_init):   
        """tabu search"""
        s_best = s_init 
        bestCandidate = copy.deepcopy(s_best)
        routes , temp_tabu = [] , []   # init
        routes.append(s_best)
        while(self.max_iters):
            self.max_iters -= 1 # Number of iterations
            neighbors = copy.deepcopy(s_best)
            for s in combinations(neighbors, 2):   
                sCandidate = self.exchange(s[0],s[1],neighbors)  # exchange number to generate candidates
                if s not in self.tabu_list and self.get_route_distance(sCandidate) < self.get_route_distance(bestCandidate):
                    bestCandidate = sCandidate
                    temp_tabu = s                           
            if self.get_route_distance(bestCandidate) < self.get_route_distance(s_best): # record the best solution 
                s_best = bestCandidate
            if  temp_tabu not in self.tabu_list:
                self.tabu_list.append(temp_tabu)
            if len(self.tabu_list) > self.maxTabuSize :
                self.tabu_list.pop(0)
            routes.append(bestCandidate)
        return s_best, routes

if __name__ == "__main__":
	np.random.seed(2020)
    customerNum = 10  # 定义多少个点
    data=GetData()
    tsp_data = data.generate_locations(num_points=customerNum+1,map_size=100)  #在100*100的图中，随机生成位置，customerNum+1 多一个depot点
    dismatrix = data.get_euclidean_distance_matrix(tsp_data.locations)

    """ Tabu : 
        disMatrix : the distance matrix from 0 to X , 0 represernt starting and stopping point。 
        for example:   disMatrix = [[0,3,4,...
        							 1,0,5,...
        							 3,5,0,...]]
       that means the distance from 0 to 0 is 0, from 0 to 1 is 3,... from 1 to 3 is 5....		
        max_iters : maximum iterations 
        maxTabuSize : maximum iterations 
    """
    tsp = Tabu(disMatrix=dismatrix ,max_iters=10,maxTabuSize=10) 
	# two different way to generate initial solution
	# num : the number of points   
    s_init = tsp.generate_initial_solution(num=customerNum,mode='greedy') # mode = "greedy"  or "random"
    print('init route : ' , s_init)
    print('init distance : ' , tsp.get_route_distance(s_init))

    start = time.time()
    best_route , routes = tsp.tabu_search(s_init)     # tabu search
    end = time.time()

    print('best route : ' , best_route)
    print('best best_distance : ' , tsp.get_route_distance(best_route))
    print('the time cost : ',end - start )

    # plot the result changes with iterations
    results=[]
    for i in routes:
        results.append(tsp.get_route_distance(i))    
    plt.plot(np.arange(len(results)) , results)
    plt.show()
    # plot the route
    data.plot_route(tsp_data.locations,[0]+best_route+[0])

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

(https://blog.csdn.net/DCXY71/article/details/110991670)

如需使用Solomon标准算例，只需更调用GetData.read_solomon() 函数，详细代码见github。

最后求解可以得到：solution=[2, 4, 3, 6, 1, 5, 9, 7, 8, 8] ,total_distance = 395.

3. 问题与总结

对于中小规模TSP的求解，求解速度和解的质量都不错，禁忌表长度似乎对求解速度影响不大。但随着点数增加，求解规模增大，所用时间明显增加，解的质量出现明显下降，且解的质量与给定初始解密切相关。

1.求解时间增加
这里我们采用的是combinations()函数，对列表元素两两组合，将会产生$C_n^2$种组合，当n变大时，组合数将非常可观，一般可采用随机组合的方式，随交换元素或是随机选择组合，没必要全部遍历。

            self.max_iters -= 1 # Number of iterations
            neighbors = copy.deepcopy(s_best)
            for s in combinations(neighbors, 2):  #这里将所有元素两两组合，若点数增加至1000点，每次将循环499500次
                sCandidate = self.exchange(s[0],s[1],neighbors) 
1
2
3
4

2.随着规模增加，解的质量出现下降
我们的候选解集合是通过两两交换规则产生，其实还是针对局部邻域进行贪婪的搜索，搜索范围有限，虽然算法通用易实现，且容易理解，但搜索性能完全依赖于邻域结构和初解，随着规模的增加尤其会陷入局部极小而无法保证全局优化型。为了实现全局优化，可尝试结合其他算法思想：例如以可控性概率接受劣解来跳出局部极小，如模拟退火算法；扩大邻域搜索结构，如TSP的2-opt扩展到k-opt。

3.禁忌长度
禁忌长度是禁忌搜索非常关键的参数，其大小直接影响到整个算法的搜索进程和的解的质量。在这里我们设置禁忌长度为一常数，采用先进先出(FIFO)的形式，但当问题比较复杂，规模比较大时，可以考虑让禁忌长度动态的自适应变化。

3.一些细节
设定合适的终止条件：（在这里我们设置了最大迭代次数）

1. 设置最大迭代次数。
2. 采用频率控制，当某个解、目标值或元素序列出现的频率高于给定的阈值，停止算法。
3. 如果在给定的迭代次数内，历史最优值未发生变化，可以停止算法。

还需要注意的是：

• 不要仅仅局限于将解作为禁忌对象，禁忌对象可以是其他任何合理的元素。例如，在前面代码中，我们禁忌的对象是交换的元素tabu_list=[(2, 9), (6, 10), (5, 7), (10, 4), (1, 3)]，但其实还可以将距离（价值）也作为禁忌的对象，当一个解的的价值已经位于禁忌表，那么就强迫算法选择不同价值的解，从而达到跳出局部极值的目的。
• 在实际搜索中，可能出现所有候选解都处于禁止选中的状态，为了能够让算法能够正常继续下去，就需要引入特赦准则（在我们这个算例中不明显）。当走投无路的时候（候选集全部位于禁忌表）或违反禁忌表可获得更优解，那么 特赦准则允许候选集违反禁忌规则，对象可以重新被选中。从这个角度来看，这相当于加强了对某个较优局部的搜索，以希望发现更好的解。
  
  

参考文献

[1] Fred Glover (1986). Future Paths for Integer Programming and Links to Artificial Intelligence. Computers and Operations Research. 13 (5): 533–549. doi:10.1016/0305-0548(86)90048-1.

[2] Fred Glover (1989). Tabu Search – Part 1. ORSA Journal on Computing. 1 (2): 190–206. doi:10.1287/ijoc.1.3.190.

[3] Fred Glover (1990). Tabu Search – Part 2. ORSA Journal on Computing. 2 (1): 4–32. doi:10.1287/ijoc.2.1.4.