1"""Example that solves the Traveling Salesman Problem using the simple compact 2formulation presented in Miller, C.E., Tucker, A.W and Zemlin, R.A. "Integer 3Programming Formulation of Traveling Salesman Problems". Journal of the ACM 47(4). 1960.""" 5 6from itertools import product 7from sys import stdout as out 8from mip import Model, xsum, minimize, BINARY 9import random as rnd 10from typing import List 11 12# names of places to visit 13places = ['Antwerp', 'Bruges', 'C-Mine', 'Dinant', 'Ghent', 14 'Grand-Place de Bruxelles', 'Hasselt', 'Leuven', 15 'Mechelen', 'Mons', 'Montagne de Bueren', 'Namur', 16 'Remouchamps', 'Waterloo'] 17 18# distances in an upper triangular matrix 19dists = [[83, 81, 113, 52, 42, 73, 44, 23, 91, 105, 90, 124, 57], 20 [161, 160, 39, 89, 151, 110, 90, 99, 177, 143, 193, 100], 21 [90, 125, 82, 13, 57, 71, 123, 38, 72, 59, 82], 22 [123, 77, 81, 71, 91, 72, 64, 24, 62, 63], 23 [51, 114, 72, 54, 69, 139, 105, 155, 62], 24 [70, 25, 22, 52, 90, 56, 105, 16], 25 [45, 61, 111, 36, 61, 57, 70], 26 [23, 71, 67, 48, 85, 29], 27 [74, 89, 69, 107, 36], 28 [117, 65, 125, 43], 29 [54, 22, 84], 30 [60, 44], 31 [97], 32 []] 33 34# number of nodes and list of vertices 35n, V = len(dists), range(len(dists)) 36 37# distances matrix 38c = [[0 if i == j 39 else dists[i][j-i-1] if j > i 40 else dists[j][i-j-1] 41 for j in V] for i in V] 42 43model = Model() 44 45# binary variables indicating if arc (i,j) is used on the route or not 46x = [[model.add_var(var_type=BINARY) for j in V] for i in V] 47 48# continuous variable to prevent subtours: each city will have a 49# different sequential id in the planned route except the first one 50y = [model.add_var() for i in V] 51 52# objective function: minimize the distance 53model.objective = minimize(xsum(c[i][j]*x[i][j] for i in V for j in V)) 54 55# constraint : leave each city only once 56for i in V: 57 model += xsum(x[i][j] for j in set(V) - {i}) == 1 58 59# constraint : enter each city only once 60for i in V: 61 model += xsum(x[j][i] for j in set(V) - {i}) == 1 62 63# subtour elimination 64for (i, j) in set(product(set(V) - {0}, set(V) - {0})): 65 model += y[i] - (n+1)*x[i][j] >= y[j]-n 66 67# running a best insertion heuristic to obtain an initial feasible solution: 68# test every node j not yet inserted in the route at every intermediate 69# position p and select the pair (j, p) that results in the smallest cost 70# increase 71seq = [0, max((c[0][j], j) for j in V)[1]] + [0] 72Vout = set(V)-set(seq) 73while Vout: 74 (j, p) = min([(c[seq[p]][j] + c[j][seq[p+1]], (j, p)) for j, p in 75 product(Vout, range(len(seq)-1))])[1] 76 77 seq = seq[:p+1]+[j]+seq[p+1:] 78 assert(seq[-1] == 0) 79 Vout = Vout - {j} 80cost = sum(c[seq[i]][seq[i+1]] for i in range(len(seq)-1)) 81print('route with cost %g built' % cost) 82 83 84# function to evaluate the cost of swapping two positions in a route in 85# constant time 86def delta(d: List[List[float]], S: List[int], p1: int, p2: int) -> float: 87 p1, p2 = min(p1, p2), max(p1, p2) 88 e1, e2 = S[p1], S[p2] 89 if p1 == p2: 90 return 0 91 elif abs(p1-p2) == 1: 92 return ((d[S[p1-1]][e2] + d[e2][e1] + d[e1][S[p2+1]]) 93 - (d[S[p1-1]][e1] + d[e1][e2] + d[e2][S[p2+1]])) 94 else: 95 return ( 96 (d[S[p1-1]][e2] + d[e2][S[p1+1]] + d[S[p2-1]][e1] + d[e1][S[p2+1]]) 97 - (d[S[p1-1]][e1] + d[e1][S[p1+1]] + d[S[p2-1]][e2] + d[e2][S[p2+1]])) 98 99 100# applying the Late Acceptance Hill Climbing 101rnd.seed(0) 102L = [cost for i in range(50)] 103sl, cur_cost, best = seq.copy(), cost, cost 104for it in range(int(1e7)): 105 (i, j) = rnd.randint(1, len(sl)-2), rnd.randint(1, len(sl)-2) 106 dlt = delta(c, sl, i, j) 107 if cur_cost + dlt <= L[it % len(L)]: 108 sl[i], sl[j], cur_cost = sl[j], sl[i], cur_cost + dlt 109 if cur_cost < best: 110 seq, best = sl.copy(), cur_cost 111 L[it % len(L)] = cur_cost 112 113print('improved cost %g' % best) 114 115model.start = [(x[seq[i]][seq[i+1]], 1) for i in range(len(seq)-1)] 116# optimizing 117model.optimize(max_seconds=30) 118 119# checking if a solution was found 120if model.num_solutions: 121 out.write('route with total distance %g found: %s' 122 % (model.objective_value, places[0])) 123 nc = 0 124 while True: 125 nc = [i for i in V if x[nc][i].x >= 0.99][0] 126 out.write(' -> %s' % places[nc]) 127 if nc == 0: 128 break 129 out.write('\n') 130model.check_optimization_results() 131