[关闭]
@Plams 2018-06-03T13:14:58.000000Z 字数 10964 阅读 791

ACO

首先采用原始未优化的蚁群算法，求解eil51.tsp问题。
eil51.tsp 含有51个点，最优解总长度426。


一个较好的结果， 约为最优解的105%。

可以看出，仍有交叉的路径，并非最优解。

2-opt ACO

配合2-opt的local search。
一个比较好的解答, 不超过最优解的100.7%。

  1. Shortest path found:
  2. [ 4 37 10 31  0 21  1 15 49  8 29 33 20 28 19 34 35  2 27 30  7 25  6 42
  3.  23 22 47  5 26 50 45 11 46  3 17 13 24 12 40 39 18 41 43 16 36 14 44 32
  4.  38  9 48  4]
  5. Minimum distance found = 428.981647

使用邻域搜索的一个附带好处是 可以大幅度减少蚂蚁数，同时不会严重影响求解质量

接下来尝试使用2-OPT_ACO来求解一个规模稍大的问题 tsp225 (3919)

  1. Shortest path found:
  2. [ 22  21  20  19 202  18  17  16  15  14  13  12  11  10   9   8  39  40
  3.    7   6   5   4   2   0 199 197   3 196 194  42  41  43  45 193 217 192
  4.   44  47 195 191 190 223 198 132 204 189 224 206  48  49  50  56  55  54
  5.   51  52  53  57  58   1  46 188  26 187 116 186 185 184 183 119 121 120
  6.  174 118 117 114 112 113 222 115  59  60  61  62  63 111 110 109  64  65
  7.   66  67  68  69  70  71  72  73  75  74  97  98  99 100 101 102 103 219
  8.  104 105 106 107 108 123 122 170 172 181 180 179 178 177 176 175 173 171
  9.  169 168 167 124 211 213 149 150 151 148 147 146 145 152 153 154 155 156
  10.  143 144 142 200 141 138 137 139 140 158 159 160 161 162 136 135 182 134
  11.  157 212 133 214 163 164 165 166 125 126 127 128 131 221 129 210 130  85
  12.   84  81  82  83 209  86  87  88  89  90  91  92  93 220  80  79  94 208
  13.   95  96  77  78  76 201  29  28  27 203  25  33  32  34 205 216 218 215
  14.   30  31  38  37  36  35  24 207  23  22]
  15. Minimum distance found = 4028.625948

Code

原始ACO

main.py

  1. import math
  3. from aco import ACO, Graph
  4. from plot import plot
  7. def distance(city1: dict, city2: dict):
  8.     return math.sqrt((city1['x'] - city2['x']) ** 2 + (city1['y'] - city2['y']) ** 2)
  11. def main():
  12.     cities = []
  13.     points = []
  14.     with open('eil51.txt') as f:
  15.         for line in f.readlines():
  16.             city = line.split(' ')
  17.             cities.append(dict(index=int(city[0]), x=int(city[1]), y=int(city[2])))
  18.             points.append((int(city[1]), int(city[2])))
  19.     cost_matrix = []
  20.     rank = len(cities)
  21.     for i in range(rank):
  22.         row = []
  23.         for j in range(rank):
  24.             row.append(distance(cities[i], cities[j]))
  25.         cost_matrix.append(row)
  26.     aco = ACO(10, 2000, 1.0, 10.0, 0.5, 10, 2)
  27.     graph = Graph(cost_matrix, rank)
  28.     path, cost = aco.solve(graph)
  29.     print('cost: {}, path: {}'.format(cost, path))
  30.     plot(points, path)
  32. if __name__ == '__main__':
  33.     for _ in range(20):
  34.         main()

aco.py

  1. import random
  2. class Graph(object):
  3.     def __init__(self, cost_matrix: list, rank: int):
  4.         """
  5.         :param cost_matrix:
  6.         :param rank: rank of the cost matrix
  7.         """
  8.         self.matrix = cost_matrix
  9.         self.rank = rank
  10.         # noinspection PyUnusedLocal
  11.         self.pheromone = [[1 / (rank * rank) for j in range(rank)] for i in range(rank)]
  14. class ACO(object):
  15.     def __init__(self, ant_count: int, generations: int, alpha: float, beta: float, rho: float, q: int,
  16.                  strategy: int):
  17.         """
  18.         :param ant_count:
  19.         :param generations:
  20.         :param alpha: relative importance of pheromone
  21.         :param beta: relative importance of heuristic information
  22.         :param rho: pheromone residual coefficient
  23.         :param q: pheromone intensity
  24.         :param strategy: pheromone update strategy. 0 - ant-cycle, 1 - ant-quality, 2 - ant-density
  25.         """
  26.         self.Q = q
  27.         self.rho = rho
  28.         self.beta = beta
  29.         self.alpha = alpha
  30.         self.ant_count = ant_count
  31.         self.generations = generations
  32.         self.update_strategy = strategy
  34.     def _update_pheromone(self, graph: Graph, ants: list):
  35.         for i, row in enumerate(graph.pheromone):
  36.             for j, col in enumerate(row):
  37.                 graph.pheromone[i][j] *= self.rho
  38.                 for ant in ants:
  39.                     graph.pheromone[i][j] += ant.pheromone_delta[i][j]
  41.     # noinspection PyProtectedMember
  42.     def solve(self, graph: Graph):
  43.         """
  44.         :param graph:
  45.         """
  46.         best_cost = float('inf')
  47.         best_solution = []
  48.         for gen in range(self.generations):
  49.             # noinspection PyUnusedLocal
  50.             ants = [_Ant(self, graph) for i in range(self.ant_count)]
  51.             for ant in ants:
  52.                 for i in range(graph.rank - 1):
  53.                     ant._select_next()
  54.                 ant.total_cost += graph.matrix[ant.tabu[-1]][ant.tabu[0]]
  55.                 if ant.total_cost < best_cost:
  56.                     best_cost = ant.total_cost
  57.                     best_solution = [] + ant.tabu
  58.                 # update pheromone
  59.                 ant._update_pheromone_delta()
  60.             self._update_pheromone(graph, ants)
  61.             # print('generation #{}, best cost: {}, path: {}'.format(gen, best_cost, best_solution))
  62.         return best_solution, best_cost
  65. class _Ant(object):
  66.     def __init__(self, aco: ACO, graph: Graph):
  67.         self.colony = aco
  68.         self.graph = graph
  69.         self.total_cost = 0.0
  70.         self.tabu = []  # tabu list
  71.         self.pheromone_delta = []  # the local increase of pheromone
  72.         self.allowed = [i for i in range(graph.rank)]  # nodes which are allowed for the next selection
  73.         self.eta = [[0 if i == j else 1 / graph.matrix[i][j] for j in range(graph.rank)] for i in
  74.                     range(graph.rank)]  # heuristic information
  75.         start = random.randint(0, graph.rank - 1)  # start from any node
  76.         self.tabu.append(start)
  77.         self.current = start
  78.         self.allowed.remove(start)
  80.     def _select_next(self):
  81.         denominator = 0
  82.         for i in self.allowed:
  83.             denominator += self.graph.pheromone[self.current][i] ** self.colony.alpha * self.eta[self.current][
  84.                                                                                             i] ** self.colony.beta
  85.         # noinspection PyUnusedLocal
  86.         probabilities = [0 for i in range(self.graph.rank)]  # probabilities for moving to a node in the next step
  87.         for i in range(self.graph.rank):
  88.             try:
  89.                 self.allowed.index(i)  # test if allowed list contains i
  90.                 probabilities[i] = self.graph.pheromone[self.current][i] ** self.colony.alpha * \
  91.                     self.eta[self.current][i] ** self.colony.beta / denominator
  92.             except ValueError:
  93.                 pass  # do nothing
  94.         # select next node by probability roulette
  95.         selected = 0
  96.         rand = random.random()
  97.         for i, probability in enumerate(probabilities):
  98.             rand -= probability
  99.             if rand <= 0:
  100.                 selected = i
  101.                 break
  102.         self.allowed.remove(selected)
  103.         self.tabu.append(selected)
  104.         self.total_cost += self.graph.matrix[self.current][selected]
  105.         self.current = selected
  107.     # noinspection PyUnusedLocal
  108.     def _update_pheromone_delta(self):
  109.         self.pheromone_delta = [[0 for j in range(self.graph.rank)] for i in range(self.graph.rank)]
  110.         for _ in range(1, len(self.tabu)):
  111.             i = self.tabu[_ - 1]
  112.             j = self.tabu[_]
  113.             if self.colony.update_strategy == 1:  # ant-quality system
  114.                 self.pheromone_delta[i][j] = self.colony.Q
  115.             elif self.colony.update_strategy == 2:  # ant-density system
  116.                 # noinspection PyTypeChecker
  117.                 self.pheromone_delta[i][j] = self.colony.Q / self.graph.matrix[i][j]
  118.             else:  # ant-cycle system
  119.                 self.pheromone_delta[i][j] = self.colony.Q / self.total_cost

plot.py

  1. import operator
  3. import matplotlib.pyplot as plt
  6. def plot(points, path: list):
  7.     x = []
  8.     y = []
  9.     for point in points:
  10.         x.append(point[0])
  11.         y.append(point[1])
  12.     # noinspection PyUnusedLocal
  13.     y = list(map(operator.sub, [max(y) for i in range(len(points))], y))
  14.     plt.plot(x, y, 'co')
  16.     for _ in range(1, len(path)):
  17.         i = path[_ - 1]
  18.         j = path[_]
  19.         # noinspection PyUnresolvedReferences
  20.         plt.arrow(x[i], y[i], x[j] - x[i], y[j] - y[i])
  22.     # noinspection PyTypeChecker
  23.     plt.xlim(0, max(x) * 1.1)
  24.     # noinspection PyTypeChecker
  25.     plt.ylim(0, max(y) * 1.1)
  26.     plt.show()

2-OPT ACO

aco_2opt.py

  1. from math import sqrt
  2. import argparse
  3. import os
  4. import numpy as np
  5. from matplotlib import pyplot as plt
  6. import time
  8. #Global variables. Still not sure what these should be
  9. SCENT_MIN = 0.1
  10. ALPHA = 1
  11. BETA = 2
  12. EVAP_COEFF = 0.01
  14. class Ant(object):
  15.     def __init__(self, node):
  16.         self.node = node
  17.         self.path = np.array([node])
  18.         self.visited = {}
  20. class Node(object):
  21.     def __init__(self, id, x, y):
  22.         self.id = int(id) - 1                               #added -1
  23.         self.x = float(x)
  24.         self.y = float(y)
  26. def build_graph(filename):
  28.     with open (filename, "r") as myfile:
  29.         data = myfile.readlines()
  31.     node_list = []
  32.     length = len(data)
  33.     print "Number of cities = %d" % length
  34.     for i in xrange(length):
  35.         line = data[i].split()
  36.         node = Node(line[0], line[2], line[1])
  37.         node_list.append(node)
  39.     distance = np.empty([length, length], dtype=float)
  40.     heuristic = np.empty([length, length], dtype=float)
  41.     scent = np.empty([length, length], dtype=float)
  42.     scent.fill(SCENT_MIN)
  44.     (distance, heuristic) = init_distances(node_list, distance)
  46.     return (node_list, distance, heuristic, scent, length)
  48. def init_distances(node_list, distances):
  50.     for node1 in node_list:
  51.         for node2 in node_list:
  52.             dist = sqrt((node1.x - node2.x)**2 + (node1.y - node2.y)**2)
  53.             distances[node1.id][node2.id] = dist                    #removed -1's
  54.     np.fill_diagonal(distances, np.nan)
  56.     heuristic = 1/distances
  57.     return (distances, heuristic)
  59. def update_probability(scent, heuristic):
  61.     probability = ((scent**ALPHA) * (heuristic**BETA)) / \
  62.             np.nansum((scent**ALPHA) * (heuristic**BETA), axis=1)[:,np.newaxis]
  65.     return probability
  67. def update_scents(ants, scent, distance):
  69.     num_nodes = len(ants[0].path)
  70.     path_travelers = {}
  72.     for ant in ants:
  73.         for i in xrange(num_nodes - 1):
  74.             src = ant.path[i]
  75.             dst = ant.path[i + 1]
  76.             if (src,dst) in path_travelers:
  77.                 path_travelers[(src,dst)] = path_travelers[(src,dst)] + 1
  78.                 path_travelers[(dst,src)] = path_travelers[(dst,src)] + 1
  79.             else:
  80.                 path_travelers[(src,dst)] = 1
  81.                 path_travelers[(dst,src)] = 1
  82.             dist = distance[src][dst]
  84.     for key, value in path_travelers.iteritems():
  85.         src = key[0]
  86.         dst = key[1]
  87.         updated_scent  = (1 - EVAP_COEFF) * scent[src][dst] + 100 * value * (1/distance[src][dst])
  88.         scent[src][dst] = max(updated_scent, SCENT_MIN)
  90.     return scent
  93. def path_distance(path, distance):
  94.     dist = 0
  95.     num_cities = len(path)
  96.     for i in xrange(num_cities - 1):
  97.         #print "i = %d" % i
  98.         src = path[i]
  99.         dst = path[i+1]
  100.         dist = dist + distance[src][dst]
  102.     return dist
  104. def update_paths(ants, probability):
  105.     num_nodes = probability.shape[0]
  106.     for ant in ants:
  107.         new_arr = np.copy(probability)
  108.         curr_node = ant.path[-1]
  109.         for i in xrange(num_nodes):
  110.             for node in ant.path:
  111.                 new_arr[curr_node][node] = 0
  113.             row_sum = np.nansum(new_arr[curr_node])
  114.             if row_sum != 0:
  115.                 new_arr[curr_node] = new_arr[curr_node] / row_sum
  116.                 next_node = np.random.choice(np.arange(0, num_nodes), p = new_arr[curr_node])
  117.                 ant.path = np.append(ant.path, next_node)
  118.             else:
  119.                 ant.path = np.append(ant.path, ant.path[0])
  121.     return ants
  123. def update_greedy_paths(ants, distance):
  125.     for ant in ants:
  126.         curr_node = ant.path[-1]
  129. def place_ants_randomly(num_ants, num_cities):
  130.     ants = np.empty(num_ants, dtype=object)
  131.     for i in xrange(num_ants):
  132.         ants[i] = Ant(np.random.random_integers(0, num_cities - 1))
  133.     return ants
  135. def shortest_path(ants, min_dist, min_path, distance):
  137.     for ant in ants:
  138.         #print ant.path
  139.         dist = 0
  140.         for i in xrange(ant.path.shape[0] - 1):
  141.             src = ant.path[i]
  142.             dst = ant.path[i + 1]
  143.             dist = dist + distance[src][dst]
  144.         if dist < min_dist:
  145.             min_path = ant.path
  146.             min_dist = dist
  147.             #print min_path
  149.     return (min_dist, min_path)
  151. def plot(nodes, path):
  152.     locations = {}
  153.     for node in nodes:
  154.         locations[node.id] = (node.x, node.y)
  155.     x = []
  156.     y = []
  157.     nums = []
  158.     for i in path:
  159.         x.append(locations[i][0])
  160.         y.append(locations[i][1])
  162.     plt.scatter(x,y)
  163.     plt.plot(x,y)
  164.     plt.xlabel("X Position")
  165.     plt.ylabel("Y Position")
  166.     plt.title("Shortest path found")
  167.     plt.hold(True)
  168.     plt.show()
  171. def two_opt_swap(route, i, k):
  172.     route_len = route.shape[0]
  173.     new_route = np.empty(route_len, dtype=int)
  174.     new_route[0] = route[0]
  176.     j = 1
  177.     for a in xrange(1, i):
  178.         new_route[j] = route[a]
  179.         j += 1
  181.     """route[i] to route[k] in reverse"""
  182.     for a in reversed(xrange(i, k+1)):
  183.         new_route[j] = route[a]
  184.         j += 1
  186.     for a in xrange(k+1, route_len):
  187.         new_route[j] = route[a]        
  188.         j += 1
  190.     return new_route
  192. def two_opt_iter(curr_path, distance):
  194.     num_nodes = curr_path.shape[0]
  195.     best_distance = path_distance(curr_path, distance)
  196.     for i in xrange(1, num_nodes - 1):
  197.         for k in xrange(i + 1, num_nodes - 1):
  198.             new_route = two_opt_swap(curr_path, i, k)
  199.             new_dist = path_distance(new_route, distance)
  200.             if (new_dist < best_distance):
  201.                 print new_dist
  202.                 return (new_route, True)
  204.     return (curr_path, False)
  207. def two_opt(curr_path, distance):
  209.     num_nodes = curr_path.shape[0]
  210.     improvement = True
  211.     while (improvement):
  212.         (curr_path, improvement) = two_opt_iter(curr_path, distance)
  214.     return curr_path
  217. def local_search(ants, distance):
  218.     """Local search on each solution"""
  219.     for ant in ants:
  220.         ant.path = two_opt(ant.path, distance)
  222.     return ants
  225. def main():
  227.     #parse arguments
  228.     parser = argparse.ArgumentParser(description='Short sample app')
  229.     parser.add_argument('-a', action='store',dest='num_ants',default=10, type=int)
  230.     parser.add_argument('-i', action='store',dest='iterations',default=10, type=int)
  231.     parser.add_argument('-f', action='store',dest='file',type=str)
  232.     args = parser.parse_args()
  235.     #Intialize graph
  236.     num_ants = args.num_ants
  237.     iterations = args.iterations
  238.     file_path = args.file
  240.     (nodes, distance, heuristic, scent, num_cities) = build_graph(file_path)
  242.     #main loop
  243.     min_path = []
  244.     min_dist = float('inf')
  246.     for a in xrange(iterations):
  248.         ants = place_ants_randomly(num_ants, num_cities)
  250.         probability = update_probability(scent, heuristic)
  252.         ants = update_paths(ants, probability)
  254.         ants = local_search(ants, distance)
  256.         (min_dist, min_path) = shortest_path(ants, min_dist, min_path, distance)
  258.         min_dist = path_distance(min_path, distance)
  260.         scent = update_scents(ants, scent, distance)
  262.     print "Shortest path found:"
  263.     print min_path
  264.     print "Minimum distance found = %f" % min_dist
  266.     plot(nodes,min_path)
  268. if __name__ == '__main__':
  269.     main()

Reference

  1. 秦东各, 王长坤. 一种基于 2-opt 算法的混合型蚁群算法[J]. 工业控制计算机, 2018, 1: 042.
  2. "高级蚁群算法求解1000以上城市的TSP问题（旅行商），附大量TSPLIB数据！" http://club.excelhome.net/thread-825313-1-1.html
  3. https://github.com/ppoffice/ant-colony-tsp
  4. https://github.com/djzurawski/aco-tsp
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注