# 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 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. In the Cutting Planes [DFJ54]
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
involving 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, s 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 ?
A subset is a *cut* in a Graph. To identify the most isolated subset 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 | ```
from mip.model import Model, xsum
from mip.constants import BINARY
from itertools import product
from networkx import minimum_cut, DiGraph
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
m.relax()
while newConstraints:
m.optimize()
print('objective value : {}'.format(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
``` |

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 21 sets the objective function and the
following tree lines include constraints enforcing one entering and one
leaving arc to be selected for each node. On line 28 we relax the
integrality constraints of variables so that the optimization performed in
line 31 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 40 and whenever a disconnected subset is found the violated inequality
is generated and included at line 42. The process repeats while new
violated inequalities are generated.

## 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 in 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*. 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_cuts()`

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 by pre-processing. Thus, direct references
to the original problem variables may be invalid. Since for variables that
remain in this model Python-MIP ensures that the names from the original
variables are preserved, it is a good practice to query again the variables
list. In our example, the relationship of the variables with the arcs of the
input graph can be inferred by examining variable names which are in the format
“`x(i,j)`

” (lines 15-17). Whenever a violated inequality is discovered,
it can be added to the solver’s engine cut pool using the
`add_cut()`

`Model`

method (lines 33 and 35).
In our example, we temporarily store the generated cuts in our
`CutPool`

object (line 30).

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 | ```
from sys import argv
from typing import List, Tuple
import networkx as nx
from tspdata import TSPData
from mip.model import Model, xsum, BINARY
from mip.callbacks import ConstrsGenerator, CutPool
class SubTourCutGenerator(ConstrsGenerator):
def __init__(self, Fl: List[Tuple[int, int]]):
self.F = Fl
def generate_constrs(self, model: Model):
G = nx.DiGraph()
r = [(v, v.x) for v in model.vars if v.name.startswith('x(')]
U = [int(v.name.split('(')[1].split(',')[0]) for v, f in r]
V = [int(v.name.split(')')[0].split(',')[1]) for v, f in r]
cp = CutPool()
for i in range(len(U)):
G.add_edge(U[i], V[i], capacity=r[i][1])
for (u, v) in F:
if u not in U or v not in V:
continue
val, (S, NS) = nx.minimum_cut(G, u, v)
if val <= 0.99:
arcsInS = [(v, f) for i, (v, f) in enumerate(r)
if U[i] in S and V[i] in S]
if sum(f for v, f in arcsInS) >= (len(S)-1)+1e-4:
cut = xsum(1.0*v for v, fm in arcsInS) <= 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
return
inst = TSPData(argv[1])
n, d = inst.n, inst.d
model = Model()
x = [[model.add_var(name='x({},{})'.format(i, j),
var_type=BINARY) for j in range(n)] for i in range(n)]
y = [model.add_var(name='y({})'.format(i),
lb=0.0, ub=n) for i in range(n)]
model.objective = xsum(d[i][j] * x[i][j] for j in range(n) for i in range(n))
for i in range(n):
model += xsum(x[j][i] for j in range(n) if j != i) == 1
model += xsum(x[i][j] for j in range(n) if j != i) == 1
for (i, j) in [(i, j) for (i, j) in
product(range(1, n), range(1, n)) if i != j]:
model += y[i] - (n + 1) * x[i][j] >= y[j] - n
F = []
for i in range(n):
(md, dp) = (0, -1)
for j in [k for k in range(n) if k != i]:
if d[i][j] > md:
(md, dp) = (d[i][j], j)
F.append((i, dp))
m.cuts_generator = SubTourCutGenerator(F)
model.optimize()
arcs = [(i, j) for i in range(n) for j in range(n) if x[i][j].x >= 0.99]
print('optimal route : {}'.format(arcs))
``` |

## Lazy Constraints¶

Python-MIP also supports the use of cut 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 exaclty as the cut generator
callback example with the difference that, besides starting with a smaller
formulation, we have to inform that the cut generator will be used to generate
lazy constraints:

1 2 3 4 5 | ```
...
m.constrs_generator = SubTourCutGenerator(F)
m.constrs_generator.lazy_constraints = True
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.