1"""Example of a Branch-and-Cut implementation for the Traveling Salesman
2Problem. Initially a compact (weak) formulation is created. This formulation is
3dynamically improved with sub-tour elimination inequalities, using a
4CutsGenerator implementation that repeatedly solves the min-cut problem between
5pairs of distant nodes to find isolated sets of nodes. Cut generation is called
6from the solver engine whenever a fractional solution is generated."""
7
8from typing import List, Tuple
9from random import seed, randint
10from itertools import product
11from math import sqrt
12import networkx as nx
13from mip import Model, xsum, BINARY, minimize, ConstrsGenerator, CutPool
14
15
16class SubTourCutGenerator(ConstrsGenerator):
17    """Class to generate cutting planes for the TSP"""
18    def __init__(self, Fl: List[Tuple[int, int]], x_, V_):
19        self.F, self.x, self.V = Fl, x_, V_
20
21    def generate_constrs(self, model: Model, depth: int = 0, npass: int = 0):
22        xf, V_, cp, G = model.translate(self.x), self.V, CutPool(), nx.DiGraph()
23        for (u, v) in [(k, l) for (k, l) in product(V_, V_) if k != l and xf[k][l]]:
24            G.add_edge(u, v, capacity=xf[u][v].x)
25        for (u, v) in F:
26            val, (S, NS) = nx.minimum_cut(G, u, v)
27            if val <= 0.99:
28                aInS = [(xf[i][j], xf[i][j].x)
29                        for (i, j) in product(V_, V_) if i != j and xf[i][j] and i in S and j in S]
30                if sum(f for v, f in aInS) >= (len(S)-1)+1e-4:
31                    cut = xsum(1.0*v for v, fm in aInS) <= len(S)-1
32                    cp.add(cut)
33                    if len(cp.cuts) > 256:
34                        for cut in cp.cuts:
35                            model += cut
36                        return
37        for cut in cp.cuts:
38            model += cut
39
40
41n = 30  # number of points
42V = set(range(n))
43seed(0)
44p = [(randint(1, 100), randint(1, 100)) for i in V]  # coordinates
45Arcs = [(i, j) for (i, j) in product(V, V) if i != j]
46
47# distance matrix
48c = [[round(sqrt((p[i][0]-p[j][0])**2 + (p[i][1]-p[j][1])**2)) for j in V] for i in V]
49
50model = Model()
51
52# binary variables indicating if arc (i,j) is used on the route or not
53x = [[model.add_var(var_type=BINARY) for j in V] for i in V]
54
55# continuous variable to prevent subtours: each city will have a
56# different sequential id in the planned route except the first one
57y = [model.add_var() for i in V]
58
59# objective function: minimize the distance
60model.objective = minimize(xsum(c[i][j]*x[i][j] for (i, j) in Arcs))
61
62# constraint : leave each city only once
63for i in V:
64    model += xsum(x[i][j] for j in V - {i}) == 1
65
66# constraint : enter each city only once
67for i in V:
68    model += xsum(x[j][i] for j in V - {i}) == 1
69
70# (weak) subtour elimination
71# subtour elimination
72for (i, j) in product(V - {0}, V - {0}):
73    if i != j:
74        model += y[i] - (n+1)*x[i][j] >= y[j]-n
75
76# no subtours of size 2
77for (i, j) in Arcs:
78    model += x[i][j] + x[j][i] <= 1
79
80# computing farthest point for each point, these will be checked first for
81# isolated subtours
82F, G = [], nx.DiGraph()
83for (i, j) in Arcs:
84    G.add_edge(i, j, weight=c[i][j])
85for i in V:
86    P, D = nx.dijkstra_predecessor_and_distance(G, source=i)
87    DS = list(D.items())
88    DS.sort(key=lambda x: x[1])
89    F.append((i, DS[-1][0]))
90
91model.cuts_generator = SubTourCutGenerator(F, x, V)
92model.optimize()
93
94print('status: %s route length %g' % (model.status, model.objective_value))
95
96# sanity check
97model.check_optimization_results()
98