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