# Developing Customized Branch-&-Cut algorithms¶

This chapter discusses some features of Python-MIP that allow the development of improved Branch-&-Cut algorithms by linking application specific routines to the generic algorithm included in the solver engine. We start providing an introduction to cutting planes and cut separation routines in the next section, following with a section describing how these routines can be embedded in the Branch-&-Cut solver engine using the generic cut callbacks of Python-MIP.

## Cutting Planes¶

In many applications there are strong formulations
that may include an exponential number of constraints. These formulations cannot be direct handled by the MIP Solver: entering all these constraints at once is usually not practical, except for very small instances. In the Cutting Planes [Dantz54] method the LP relaxation is solved and only constraints which are *violated* are inserted. The model is re-optimized and at each iteration a stronger formulation is obtained until no more violated inequalities are found. The problem of discovering which are the missing violated constraints is also an optimization problem (finding *the most* violated inequality) and it is called the *Separation Problem*.

As an example, consider the Traveling Salesman Problem. The compact formulation (`tsp-label`

) is a *weak* formulation: dual bounds produced
at the root node of the search tree are distant from the optimal solution cost
and improving these bounds requires a potentially intractable number of
branchings. In this case, the culprit are the sub-tour elimination constraints
linking variables \(x\) and \(y\). A much stronger TSP formulation
can be written as follows: consider a graph \(G=(N,A)\) where \(N\) is
the set of nodes and \(A\) is the set of directed edges with associated
traveling costs \(c_a \in A\). Selection of arcs is done with binary
variables \(x_a \,\,\, \forall a \in A\). Consider also that edges arriving
and leaving a node \(n\) are indicated in \(A^+_n\) and \(A^-_n\),
respectively. The complete formulation follows:

The third constraints are sub-tour elimination constraints. Since these
constraints are stated for *every subset* of nodes, the number of these
constraints is \(O(2^{|N|})\). These are the constraints that will be
separated by our cutting pane algorithm. As an example, consider the following
graph:

The optimal LP relaxation of the previous formulation without the sub-tour elimination constraints has cost 237:

As it can be seen, there are tree disconnected sub-tours. Two of these include only two nodes. Forbidding sub-tours of size 2 is quite easy: in this case we only need to include the additional constraints: \(x_{(d,e)}+x_{(e,d)}\leq 1\) and \(x_{(c,f)}+x_{(f,c)}\leq 1\).

Optimizing with these two additional constraints the objective value increases to 244 and the following new solution is generated:

Now there are sub-tours of size 3 and 4. Let’s consider the sub-tour defined by
nodes \(S=\{a,b,g\}\). The valid inequality for \(S\) is:
\(x_{(a,g)} + x_{(g,a)} + x_{(a,b)} + x_{(b,a)} + x_{(b,g)} + x_{(g,b)} \leq 2\).
Adding this cut to our model increases the objective value to 261, a significant improvement. In our example, the visual identification of the isolated subset is easy, but how to automatically identify these subsets efficiently in the general case ? Isolated subsets can be identified when a *cut* is found in the graph defined by arcs active in the unfeasible solution. To identify the most isolated subsets we just have to solve the Minimum cut problem in graphs.
In python you can use the networkx min-cut module. The following code implements a cutting plane algorithm for the asymmetric traveling
salesman problem:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | ```
from itertools import product
from networkx import minimum_cut, DiGraph
from mip import Model, xsum, BINARY, OptimizationStatus, CutType
N = ["a", "b", "c", "d", "e", "f", "g"]
A = { ("a", "d"): 56, ("d", "a"): 67, ("a", "b"): 49, ("b", "a"): 50,
("f", "c"): 35, ("g", "b"): 35, ("g", "b"): 35, ("b", "g"): 25,
("a", "c"): 80, ("c", "a"): 99, ("e", "f"): 20, ("f", "e"): 20,
("g", "e"): 38, ("e", "g"): 49, ("g", "f"): 37, ("f", "g"): 32,
("b", "e"): 21, ("e", "b"): 30, ("a", "g"): 47, ("g", "a"): 68,
("d", "c"): 37, ("c", "d"): 52, ("d", "e"): 15, ("e", "d"): 20,
("d", "b"): 39, ("b", "d"): 37, ("c", "f"): 35, }
Aout = {n: [a for a in A if a[0] == n] for n in N}
Ain = {n: [a for a in A if a[1] == n] for n in N}
m = Model()
x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A}
m.objective = xsum(c * x[a] for a, c in A.items())
for n in N:
m += xsum(x[a] for a in Aout[n]) == 1, "out({})".format(n)
m += xsum(x[a] for a in Ain[n]) == 1, "in({})".format(n)
newConstraints = True
while newConstraints:
m.optimize(relax=True)
print("status: {} objective value : {}".format(m.status, m.objective_value))
G = DiGraph()
for a in A:
G.add_edge(a[0], a[1], capacity=x[a].x)
newConstraints = False
for (n1, n2) in [(i, j) for (i, j) in product(N, N) if i != j]:
cut_value, (S, NS) = minimum_cut(G, n1, n2)
if cut_value <= 0.99:
m += (xsum(x[a] for a in A if (a[0] in S and a[1] in S)) <= len(S) - 1)
newConstraints = True
if not newConstraints and m.solver_name.lower() == "cbc":
cp = m.generate_cuts([CutType.GOMORY, CutType.MIR,
CutType.ZERO_HALF, CutType.KNAPSACK_COVER])
if cp.cuts:
m += cp
newConstraints = True
``` |

Lines 6-13 are the input data. Nodes are labeled with letters in a list
`N`

and a dictionary `A`

is used to store the weighted directed
graph. Lines 14 and 15 store output and input arcs per node. The mapping of
binary variables \(x_a\) to arcs is made also using a dictionary in line
18. Line 20 sets the objective function and the following tree lines include
constraints enforcing one entering and one leaving arc to be selected for each
node. Line 29 will only solve the LP relaxation and the separation routine can
be executed. Our separation routine is executed for each pair or nodes at line
38 and whenever a disconnected subset is found the violated inequality is
generated and included at line 40. The process repeats while new violated
inequalities are generated.

Python-MIP also supports the automatic generation of cutting planes, i.e., cutting planes that can be generated for any model just considering integrality constraints. Line 43 triggers the generation of these cutting planes with the method `generate_cuts()`

when our sub-tour elimination procedure does not finds violated sub-tour elimination inequalities anymore.

## Cut Callback¶

The cutting plane method has some limitations: even though the first rounds of
cuts improve significantly the lower bound, the overall number of iterations
needed to obtain the optimal integer solution may be too large. Better results
can be obtained with the Branch-&-Cut algorithm, where cut generation is
*combined* with branching. If you have an algorithm like the one included in
the previous Section to separate inequalities for your application you can
combine it with the complete BC algorithm implemented in the solver engine
using *callbacks*. Cut generation callbacks (CGC) are called at each node of
the search tree where a fractional solution is found. Cuts are generated in the callback and returned to the MIP solver engine which adds these cuts to the *Cut Pool*. These cuts are merged with the cuts generated with the solver
builtin cut generators and a *subset* of these cuts is included to the LP
relaxation model. Please note that in the Branch-&-Cut algorithm context cuts
are *optional* components and only those that are classified as *good* cuts by
the solver engine will be accepted, i.e., cuts that are too dense and/or have
a small violation could be discarded, since the cost of solving a much larger
linear program may not be worth the resulting bound improvement.

When using cut callbacks be sure that cuts are used only to *improve* the LP
relaxation but not to *define* feasible solutions, which need to be defined by
the initial formulation. In other words, the initial model without cuts may be
*weak* but needs to be *complete* 1. In the case of TSP, we can include the weak sub-tour elimination constraints presented in `tsp-label`

in
the initial model and then add the stronger sub-tour elimination constraints
presented in the previous section as cuts.

In Python-MIP, CGC are implemented extending the `ConstrsGenerator`

class. The following example implements the previous cut separation algorithm as a `ConstrsGenerator`

class and includes it as a cut generator for the branch-and-cut solver engine. The method that needs to be implemented in this class is the `generate_constrs()`

procedure. This method receives as parameter the object `model`

of type `Model`

. This object must be used to query the fractional values of the model `vars`

, using the `x`

property. Other model properties can be queried, such as the problem constraints (`constrs`

). Please note that, depending on which solver engine you use, some variables/constraints from the original model may have been removed in the pre-processing phase. Thus, direct references to the original problem variables may be invalid. The method `translate()`

(line 15) translates references of variables from the original model to references of variables in the model received in the callback procedure. Whenever a violated inequality is discovered, it can be added to the model using the `+=`

operator (line 31). In our example, we temporarily store the generated cuts in a `CutPool`

object (line 25) to discard repeated cuts that eventually are found.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 | ```
from typing import List, Tuple
from random import seed, randint
from itertools import product
from math import sqrt
import networkx as nx
from mip import Model, xsum, BINARY, minimize, ConstrsGenerator, CutPool
class SubTourCutGenerator(ConstrsGenerator):
"""Class to generate cutting planes for the TSP"""
def __init__(self, Fl: List[Tuple[int, int]], x_, V_):
self.F, self.x, self.V = Fl, x_, V_
def generate_constrs(self, model: Model):
xf, V_, cp, G = model.translate(self.x), self.V, CutPool(), nx.DiGraph()
for (u, v) in [(k, l) for (k, l) in product(V_, V_) if k != l and xf[k][l]]:
G.add_edge(u, v, capacity=xf[u][v].x)
for (u, v) in F:
val, (S, NS) = nx.minimum_cut(G, u, v)
if val <= 0.99:
aInS = [(xf[i][j], xf[i][j].x)
for (i, j) in product(V_, V_) if i != j and xf[i][j] and i in S and j in S]
if sum(f for v, f in aInS) >= (len(S)-1)+1e-4:
cut = xsum(1.0*v for v, fm in aInS) <= len(S)-1
cp.add(cut)
if len(cp.cuts) > 256:
for cut in cp.cuts:
model += cut
return
for cut in cp.cuts:
model += cut
n = 30 # number of points
V = set(range(n))
seed(0)
p = [(randint(1, 100), randint(1, 100)) for i in V] # coordinates
Arcs = [(i, j) for (i, j) in product(V, V) if i != j]
# distance matrix
c = [[round(sqrt((p[i][0]-p[j][0])**2 + (p[i][1]-p[j][1])**2)) for j in V] for i in V]
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 first one
y = [model.add_var() for i in V]
# objective function: minimize the distance
model.objective = minimize(xsum(c[i][j]*x[i][j] for (i, j) in Arcs))
# constraint : leave each city only once
for i in V:
model += xsum(x[i][j] for j in V - {i}) == 1
# constraint : enter each city only once
for i in V:
model += xsum(x[j][i] for j in V - {i}) == 1
# (weak) subtour elimination
# subtour elimination
for (i, j) in product(V - {0}, V - {0}):
if i != j:
model += y[i] - (n+1)*x[i][j] >= y[j]-n
# no subtours of size 2
for (i, j) in Arcs:
model += x[i][j] + x[j][i] <= 1
# computing farthest point for each point, these will be checked first for
# isolated subtours
F, G = [], nx.DiGraph()
for (i, j) in Arcs:
G.add_edge(i, j, weight=c[i][j])
for i in V:
P, D = nx.dijkstra_predecessor_and_distance(G, source=i)
DS = list(D.items())
DS.sort(key=lambda x: x[1])
F.append((i, DS[-1][0]))
model.cuts_generator = SubTourCutGenerator(F, x, V)
model.optimize()
print('status: %s route length %g' % (model.status, model.objective_value))
``` |

## Lazy Constraints¶

Python-MIP also supports the use of constraint generators to produce *lazy
constraints*. Lazy constraints are dynamically generated, just as cutting
planes, with the difference that lazy constraints are also applied to *integer
solutions*. They should be used when the initial formulation is *incomplete*.
In the case of our previous TSP example, this approach allow us to use in the
initial formulation only the degree constraints and add all required sub-tour
elimination constraints on demand. Auxiliary variables \(y\) would also not
be necessary. The lazy constraints TSP example is exactly as the cut generator
callback example with the difference that, besides starting with a smaller
formulation, we have to inform that the constraint generator will be used to generate lazy constraints using the model property `lazy_constrs_generator`

.

```
...
model.cuts_generator = SubTourCutGenerator(F, x, V)
model.lazy_constrs_generator = SubTourCutGenerator(F, x, V)
model.optimize()
...
```

## Providing initial feasible solutions¶

The Branch-&-Cut algorithm usually executes faster with the availability of an integer feasible solution: an upper bound for the solution cost improves its ability of pruning branches in the search tree and this solution is also used in local search MIP heuristics. MIP solvers employ several heuristics for the automatically production of these solutions but they do not always succeed.

If you have some problem specific heuristic which can produce an initial
feasible solution for your application then you can inform this solution to the
MIP solver using the `start`

model property. Let’s
consider our TSP application (`tsp-label`

). If the graph is complete,
i.e. distances are available for each pair of cities, then *any* permutation
\(\Pi=(\pi_1,\ldots,\pi_n)\) of the cities \(N\) can be used as an
initial feasible solution. This solution has exactly \(|N|\) \(x\)
variables equal to one indicating the selected arcs: \(((\pi_1,\pi_2),
(\pi_2,\pi_3), \ldots, (\pi_{n-1},\pi_{n}), (\pi_{n},\pi_{1}))\). Even though
this solution is obvious for the modeler, which knows that binary variables of
this model refer to arcs in a TSP graph, this solution is not obvious for the
MIP solver, which only sees variables and a constraint matrix. The following
example enters an initial random permutation of cities as initial feasible
solution for our TSP example, considering an instance with `n`

cities,
and a model `model`

with references to variables stored in a matrix
`x[0,...,n-1][0,..,n-1]`

:

1 2 3 4 | ```
from random import shuffle
S=[i for i in range(n)]
shuffle(S)
model.start = [(x[S[k-1]][S[k]], 1.0) for k in range(n)]
``` |

The previous example can be integrated in our TSP example (`tsp-label`

)
by inserting these lines before the `model.optimize()`

call. Initial
feasible solutions are informed in a list (line 4) of `(var, value)`

pairs. Please note that only the original non-zero problem variables need to be
informed, i.e., the solver will automatically compute the values of the
auxiliary \(y\) variables which are used only to eliminate sub-tours.

Footnotes

- 1
If you want to initally enter an incomplete formulation than see the next sub-section on Lazy-Constraints.