Branch & Cut |
Travelling Salesperson Problem (TSP)
- We have a network of cities that we want to visit.
- We start in one city and want to visit all other cities.
- We need to visit each city only once and must return to the starting city at the end.
- We want to find the shortest route.
Example: The Odyssey
Odysseus has to visit 22 places in the Mediterranean region
latlon = [
(38.24, 20.42),
(39.57, 26.15),
(40.56, 25.32),
(36.26, 23.12),
(33.48, 10.54),
(37.56, 12.19),
(38.42, 13.11),
(37.52, 20.44),
(41.23, 09.10),
(41.17, 13.05),
(36.08, -5.21),
(38.47, 15.13),
(38.15, 15.35),
(37.51, 15.17),
(35.49, 14.32),
(39.36, 19.56),
(38.09, 24.36),
(36.09, 23.00),
(40.44, 13.57),
(40.33, 14.15),
(40.37, 14.23),
(37.57, 22.56),
]
We can model this as a graph $G=(V, E)$.
- Each node $i \in V$ represents a location.
- The edge costs $c_{(i, j)}$ represent the distances.
- In the Odysseus example, we can take the straight-line distance between the locations as the distance.
- Here, the graph would then be fully connected.
from geopy.distance import distance
distance(latlon[0], latlon[1]).km
Which variables make sense to model this as an IP?
As with other graph problems, a binary variable for each edge: \[ x_{ij} \in \{0, 1\} \quad \forall (i, j) \in E \]
For each edge, we know the distance, so we want to minimize the total distance.
\[ \min \sum_{(i, j) \in E} c_{ij} x_{ij} \]
What constraints do we need?
Each node is left exactly once:
\[ \sum_{(i,j) \in E} x_{ij} = 1 \quad \forall i \in V \]
Each node is visited exactly once:
\[ \sum_{(j, i) \in E} x_{ji} = 1 \quad \forall i \in V \]
Are we done with this already?
import mip, itertools
model = mip.Model()
model.objective.sense = mip.MINIMIZE
variables = {}
for p1, p2 in itertools.product(latlon, repeat=2):
if p1 != p2:
variables[p1, p2] = model.add_var(var_type=mip.BINARY, obj= distance(p1, p2).km)
for p in latlon:
outflow = sum(variables[p, p2] for p2 in latlon if p != p2)
model.add_constr(outflow == 1)
inflow = sum(variables[p2, p] for p2 in latlon if p != p2)
model.add_constr(inflow == 1)
model.optimize()
The solution consists of small subtours.
So, in order to find a correct solution, we need to do something else.
Idea: For each node, count at which position it appears in the route.
- We choose an arbitrary node $v$ as the first one.
- Every other node is assigned a variable indicating its position:
\[ y_i \geq 0 \quad \forall i \in V \setminus \{v\} \]
In order for the successor of $v$ to be assigned the correct value:
\[ y_i \geq x_{vi} \quad \forall i \in V \setminus \{v\} \]
If $j$ is the successor of $i$, the following should hold:
\[ y_j \geq y_i + 1 \]However, we need to deactivate this inequality if $j$ is not actually the successor of $i$:
\[ y_j \geq y_i + 1 {\color{red} - M \cdot (1 - x_{ij})} \]No variable can take a value larger than $|V|$, therefore we can set $M = |V|$:
\[ y_j \geq y_i + 1 - |V| \cdot (1 - x_{ij}) \\[2mm] \forall i, j \in V \setminus \{v\}, i \neq j \]
\[ \begin{align*} \min \quad & \sum_{(i, j) \in E} c_{ij} x_{ij} \\ \text{s.t} \quad & \sum_{(i,j) \in E} x_{ij} = 1 && \forall i \in V \\ & \sum_{(j, i) \in E} x_{ji} = 1 && \forall i \in V \\ & y_j \geq y_i + 1 - |V| \cdot (1 - x_{ij}) && \forall i, j \in V \setminus \{v\}, i \neq j \\ & x_{ij} \in \{0, 1\} && \forall (i, j) \in E \\ & y_i \geq 0 && \forall i \in V \setminus \{v\} \end{align*} \]
order_vars = {p: model.add_var(lb=0) for p in latlon[1:]}
for p1, p2 in itertools.product(latlon, latlon[1:]):
if p1 == latlon[0]:
model.add_constr( order_vars[p2] - variables[p1,p2] >= 0 )
elif p1 != p2:
model.add_constr( order_vars[p2] >= order_vars[p1] + 1 - len(latlon) * (1 - variables[p1,p2]) )
Now we have the correct solution.
Disadvantage of this model: The Big-M constraints slow down the model.
Can we find better constraints to exclude subtours?
How many selected edges are there within the nodes of a subtour?
How many should there have been?
For each subset of nodes, the following must hold:
\[ \sum_{\substack{(i, j) \in E\\ i \in S, j \in S}} x_{ij} \leq |S| - 1 \quad \forall S \subset V, |S| \geq 2 \]
These are good constraints.
But unfortunately too many:
Number of subsets of $V$ is $2^{|V|}$
We can try to exclude all 2-subtours:
\[ x_{i,j} + x_{j,i} \leq 1 \quad \forall i, j \in V, i \neq j \]
for p1, p2 in itertools.product(latlon, repeat=2):
if p1 != p2:
model.add_constr(variables[p1, p2] + variables[p2, p1] <= 1)
The solution is at least better than before with 2 sub-routes.
How do we prevent larger sub-tours?
Idea: Only generate constraints when necessary.
We can pass a callback to the solver, which generates new constraints for a current solution.
We need to figure out how to find a subtour.
Which nodes are those of a subtour?
Those on one side of a cut with a value less than 1.
So we are looking for a minimum cut (regardless of which nodes) with a value less than 1.
One strategy: Look for the minimum s-t-cut between $v$ and all the other nodes.
Advantage: we can find the average even when the solution is still fractional.
The ConstrsGenerator defines the callback:
class SubTourCutGenerator(mip.ConstrsGenerator):
def generate_constrs(self, node_m: mip.Model, depth: int = 0, npass: int = 0):
t_vars = node_m.translate(variables)
Through preprocessing, the variables in the node are modified and need to be translated back.
Finding a new constraint:
g = nx.DiGraph()
for (p1, p2), v in t_vars.items():
g.add_edge(p1, p2, capacity=v.x)
s = latlon[0]
for t in latlon[1:]:
val, (S, NS) = nx.minimum_cut(g, s, t)
if val < 0.99:
node_m.add_constr(
sum(t_vars[p1, p2] for p1, p2 in itertools.product(S, repeat=2) if p1 != p2) <= len(S) - 1)
node_m.add_constr(
sum(t_vars[p1, p2] for p1, p2 in itertools.product(NS, repeat=2) if p1 != p2) <= len(NS) - 1)
Now we need to register the callback in the model.
model.cuts_generator = SubTourCutGenerator()
model.lazy_constrs_generator = SubTourCutGenerator()
It also makes sense to exclude the 2-edge subtours directly with the ConstrsGenerator
Each callback invocation takes time, so we save some of them.
Now we have an optimal solution and are significantly faster than with the auxiliary variables.
There are two types of generated constraints.
- Cuts:
cuts_generator - Lazy Constraints:
lazy_constrs_generator
Cuts
Optional constraints that enhance the LP relaxation.
Called with fractional solutions.
Lazy Constraints
Relevant constraints that could render a solution infeasible.
We are called with integer solutions.
Disadvantage of Lazy Constraints:
The solver can do less preprocessing
as it does not know the entire model.
May potentially cost performance.
In our case, we can have the best of both worlds.
- Helper variables to ensure the model is complete from the start.
- Cut Generator to generate better constraints if needed.
- No lazy constraints.
With this code we can already solve quite large TSPs optimally in acceptable time:
Generating constraints on demand can often bring performance benefits.
Even if the number of constraints is not exponential but only cubic, for example.
There is an analogous concept for generating variables only when needed:
Column Generation
However, this is a bit more cumbersome.