Travel salesman problem python

Travelling Salesman Problem (TSP) with Python

In this blog we shall discuss on the Travelling Salesman Problem (TSP) — a very famous NP-hard problem and will take a few attempts to solve it (either by considering special cases such as Bitonic TSP and solving it efficiently or by using algorithms to improve runtime, e.g., using Dynamic programming, or by using approximation algorithms, e.g., for Metric TSP and heuristics, to obtain not necessarily optimal but good enough solutions, e.g., with Simulated Annealing and Genetic Algorithms) and work on the corresponding python implementations. Few of the problems discussed here appeared as programming assignments in the Coursera course Advanced Algorithms and Complexity and some of the problem statements are taken from the course.

Improving the runtime of the Travelling Salesman Problem with Dynamic Programming

In this problem we shall deal with a classical NP-complete problem called Traveling Salesman Problem. Given a graph with weighted edges, you need to find the shortest cycle visiting each vertex exactly once. Vertices correspond to cities. Edges weights correspond to the cost (e.g., time) to get from one vertex to another one. Some vertices may not be connected by an edge in the general case.

• Here we shall use dynamic programming to solve TSP: instead of solving one problem we will solve a collection of (overlapping) subproblems.
• A subproblem refers to a partial solution
• A reasonable partial solution in case of TSP is the initial part of a cycle
• To continue building a cycle, we need to know the last vertex as well as the set of already visited vertices
• It will be convenient to assume that vertices are integers from 1 to n and that the salesman starts his trip in (and also returns back to) vertex 1.
• The following figure shows the Dynamic programming subproblems, the recurrence relation and the algorithm for TSP with DP.

Implementation tips

• In order to iterate through all subsets of , it will be helpful to notice that there is a natural one-to-one correspondence between integers in the range from 0 and 2^n − 1 and subsets of : k ↔ .
• For example, k = 1 (binary 001) corresponds to the set , where k = 5 (binary 101) corresponds to the set
• In order to find out the integer corresponding to S − (for j ∈ S), we need to flip the j-th bit of k (from 1 to 0). For this, in turn, we can compute a bitwise XOR of k and 2^j (that has 1 only in j-th position)
• In order to compute the optimal path along with the cost, we need to maintain back-pointers to store the path.

The following python code shows an implementation of the above algorithm.

 
import numpy as np from itertools import combinationsdef TSP(G): n = len(G) C = [[np.inf for _ in range(n)] for __ in range(1  1 for size in range(1, n): for S in combinations(range(1, n), size): S = (0,) + S k = sum([1 ][j] all_index = (1 

  
The following animation shows how the least cost solution cycle is computed with the DP for a graph with 4 vertices. Notice that in order to represent C(S,i) from the algorithm, the vertices that belong to the set S are colored with red circles, the vertex i where the path that traverses through all the nodes in S ends at is marked with a red double-circle.
  
The next animation also shows how the DP table gets updated. The DP table for a graph with 4 nodes will be of size 2⁴ X 4, since there are 2⁴=16 subsets of the vertex set V= and a path going through a subset of the vertices in V may end in any of the 4 vertex.
  
The transposed DP table is shown in the next animation, here the columns correspond to the subset of the vertices and rows correspond to the vertex the TSP ends at.

 

  
The following animation shows how the least cost solution cycle is computed with the DP for a graph with 5 nodes.

 

The following animation / figure shows the TSP optimal path is computed for increasing number of nodes (where the weights for the input graphs are randomly generated) and the exponential increase in the time taken.
   
Solving TSP with Integer Linear Program


Solving with the mip package using the following python code, produces the output shown by the following animation, for a graph with randomly generated edge-weights.
 
from mip import Model, xsum, minimize, BINARY def TSP_ILP(G): start = time() V1 = range(len(G)) n, V = len(G), set(V1) model = Model() # binary variables indicating if arc (i,j) is used # on the route or not x = [[model.add_var(var_type=BINARY) for j in V] for i in V] # continuous variable to prevent subtours: each city will have a # different sequential id in the planned route except the 1st one y = [model.add_var() for i in V] # objective function: minimize the distance model.objective = minimize(xsum(G[i][j]*x[i][j] \ for i in V for j in V)) # constraint : leave each city only once for i in V: model += xsum(x[i][j] for j in V - ) == 1 # constraint : enter each city only once for i in V: model += xsum(x[j][i] for j in V - ) == 1 # subtour elimination for (i, j) in product(V - , V - ): if i != j: model += y[i] - (n+1)*x[i][j] >= y[j]-n # optimizing model.optimize() # checking if a solution was found if model.num_solutions: print('Total distance <>'.format(model.objective_value)) nc = 0 # cycle starts from vertex 0 cycle = [nc] while True: nc = [i for i in V if x[nc][i].x >= 0.99][0] cycle.append(nc) if nc == 0: break return (model.objective_value, cycle)
 
The constraint to prevent the subtours to appear in the solution is necessary, if we run without the constraint, we get a solution with subtours instead of a single cycle going through all the nodes, as shown below:
 
  
 
Comparing with Dynamic programming based solution, we can see that ILP is much more efficient for higher n values.
 
 
 
Bitonic TSP
 
 
 
The following python code snippet implements the above DP algorithm.
 
 
def dist(P, i, j): return np.sqrt((P[i][0]-P[j][0])**2+(P[i][1]-P[j][1])**2) def BTSP(P): n = len(P) D = np.ones((n,n))*np.inf path = np.ones((n,n), dtype=int)*(-1) D[n-2,n-1] = dist(P, n-2, n-1) path[n-2,n-1] = n-1 for i in range(n-3,-1,-1): m = np.inf for k in range(i+2,n): if m > D[i+1,k] + dist(P,i,k): m, mk = D[i+1,k] + dist(P,i,k), k D[i,i+1] = m path[i,i+1] = mk for j in range(i+2,n): D[i,j] = D[i+1,j] + dist(P,i,i+1) path[i,j] = i+1 D[0,0] = D[0,1] + dist(P,0,1) path[0,0] = 1 return D, path def get_tsp_path(path, i, j, n): if n < 0: return [] if i 
The following animation shows how the DP table is computed and the optimal path for Bitonic TSP is constructed. It also shows the final optimal path.
   
2-OPT Approximation Algorithm for Metric TSP

The next code snippet implements the above 2-OPT approximation algorithm.
 
import numpy as np import queue def dfs(adj, x): visited = [False]*len(adj) stack = [x] visited[x] = True path = [] while len(stack) > 0: u = stack.pop(-1) path.append(u) for v in adj[u]: if not visited[v]: stack.append(v) visited[v] = True return pathdef mst(adj): inf = np.inf c = [inf]*n s = 0 c[s] = 0 visited = [False]*n parent = [None]*n h = queue.PriorityQueue() for v in range(n): h.put((c[v], v)) edges = [] while not h.empty(): w, u = h.get() if visited[u]: continue visited[u] = True if parent[u] != None: edges.append((parent[u], u)) for v in range(n): if v == u: continue if (not visited[v]) and (c[v] > adj[u][v]): c[v] = adj[u][v] parent[v] = u h.put((c[v], v)) adj = [[] for _ in range(n)] for i in range(n): if parent[i] != None: adj[parent[i]].append(i) path = dfs(adj, 0) path += [path[0]] return path
 
The following animation shows the TSP path computed with the above approximation algorithm and compares with the OPT path computed using ILP for 20 points on 2D plane. The MST is computed with Prim’s algorithm.
 
TSP with Simulated Annealing
 
 
 
The following python code snippet shows how to implement the Simulated Annealing to solve TSP, here G represents the adjacency matrix of the input graph.
 
 
def TSP_SA(G): s = list(range(len(G))) c = cost(G, s) ntrial = 1 T = 30 alpha = 0.99 while ntrial 
The following animations show how the algorithm works:
The following animation shows the TSP path computed with SA for 100 points in 2D.
TSP with Genetic Algorithm

Here in the following implementation of the above algorithm we shall have the following assumptions:

  
We shall assume the crossover rate is 1.0, i.e., all individuals in a population participate in crossover. The mutation probability to be used is 0.1.
  
With each crossover operation between two parent chromosomes, couple of children are generated, cant just swap portions of parents chromosomes, need to be careful to make sure that the offspring represents valid TSP path.
  
Mutation is similar to swap operation implemented earlier.
  
For each generation we shall keep a constant k=20 (or 30) chromosomes (representing candidate solutions for TSP).
  
The fitness function will be the cost of the TSP path represented by each chromosome. Hence, we want to minimize the value of the fitness function — i.e., less the value of a chromosome, more fit is it to survive.
  
We shall use rank selection, i.e., after crossover and mutation, only the top k fittest offspring (i.e., with least fitness function value) will survive for the next generation.
  
The following python code shows the implementation of the above algorithm with the above assumptions.

 
import numpy as np def do_crossover(s1, s2, m): s1, s2 = s1.copy(), s2.copy() c1 = s2.copy() for i in range(m, len(s1)): c1.remove(s1[i]) for i in range(m, len(s1)): c1.append(s1[i]) c2 = s1.copy() for i in range(m, len(s2)): c2.remove(s2[i]) for i in range(m, len(s2)): c2.append(s2[i]) return (c1, c2) def do_mutation(s, m, n): i, j = min(m, n), max(m, n) s1 = s.copy() while i < j: s1[i], s1[j] = s1[j], s1[i] i += 1 j -= 1 return s1 def compute_fitness(G, s): l = 0 for i in range(len(s)-1): l += G[s[i]][s[i+1]] l += G[s[len(s)-1]][s[0]] return l def get_elite(G, gen, k): gen = sorted(gen, key=lambda s: compute_fitness(G, s)) return gen[:k] def TSP_GA(G, k=20, ntrial = 200): n_p = k mutation_prob = 0.1 gen = [] path = list(range(len(G))) while len(gen) < n_p: path1 = path.copy() np.random.shuffle(path1) if not path1 in gen: gen.append(path1) for trial in range(ntrial): gen = get_elite(G, gen, k) gen_costs = [(round(compute_fitness(G, s),3), s) \ for s in gen] next_gen = [] for i in range(len(gen)): for j in range(i+1, len(gen)): c1, c2 = do_crossover(gen[i], gen[j], \ np.random.randint(0, len(gen[i]))) next_gen.append(c1) next_gen.append(c2) if np.random.rand() < mutation_prob: m = np.random.randint(0, len(gen[i])) while True: n = np.random.randint(0, len(gen[i])) if m != n: break c = do_mutation(gen[i], m, n) next_gen.append(c) gen = next_gen
 
The following animation shows the TSP path computed with GA for 100 points in 2D.
 
Источник