# A Basic Branch and Bound Solver in Python using Cvxpy

Branch and bound is a useful problem solving technique. The idea is, if you have a minimization problem you want to solve, maybe there is a way to relax the constraints to an easier problem. If so, the solution of the easier problem is a lower bound on the possible solution of the hard problem. If the solution of the easier problem just so happens to also obey the more constrained hard problem, then it must also be the solution to the hard problem. You can also use the lower bound coming from a relaxed problem to prune your search tree for the hard problem. If even the relaxed problem doesn’t beat the current best found, don’t bother going down that branch.

A standard place this paradigm occurs is in mixed integer programming. The relaxation of a binary constraint (either 0 or 1) can be all the values in between (any number between 0 and 1). If this relaxed problem can be expressed in a form amenable to a solver like a linear programming solver, you can use that to power the branch and bound search, also using returned solutions for possible heuristics.

I built a basic version of this that uses cvxpy as the relaxed problem solver. Cvxpy already has much much faster mixed integer solvers baked in (which is useful to make sure mine is returning correct results), but it was an interesting exercise. The real reason I’m toying around is I kind of want the ability to add custom branching heuristics or inspect and maintain the branch and bound search tree, which you’d need to get into the more complicated guts of the solvers bound to cvxpy to get at. Julia might be a better choice.

A somewhat similar (and better) project is https://github.com/oxfordcontrol/miosqp which doesn’t use cvxpy explicitly, but does have the branch and bound control in the python layer of the solver. There are also other projects that can use fairly arbitrary solvers like Bonmin

As a toy problem I’m using a knapsack problem where we have objects of different sizes and different values. We want to maximize the value while keeping the total size under the capacity of the bag. This can be phrased linearly like so: $ \max v \cdot x$ s.t. $ \sum_i s_i x_i<= capacity $, $ x \in {0,1}$. The basic heuristic I’m using is to branch on variables that are either 0 or 1 in even the relaxed solution. The alternative branch hopefully gets pruned fast.

```
import cvxpy as cvx
import copy
from heapq import *
import numpy as np
import itertools
counter = itertools.count()
class BBTreeNode():
def __init__(self, vars = set(), constraints = [], objective=0, bool_vars=set()):
self.vars = vars
self.constraints = constraints
self.objective = objective
self.bool_vars = bool_vars
self.children = []
def buildProblem(self):
prob = cvx.Problem(cvx.Minimize(self.objective), self.constraints) #i put Minimize, just so you know that I'm assuming it
return prob
def is_integral(self):
return all([abs(v.value - 1) <= 1e-3 or abs(v.value - 0) <= 1e-3 for v in self.bool_vars])
def branch(self):
children = []
for b in [0,1]:
n1 = copy.deepcopy(self) #yeesh. Not good performance wise, but is simple implementation-wise
v = n1.heuristic() #dangerous what if they don't do the same one? I need to do it here though because I need access to copied v.
n1.constraints.append( v == b ) # add in the new binary constraint
n1.children = []
n1.bool_vars.remove(v) #remove binary constraint from bool var set
n1.vars.add(v) #and add it into var set for later inspection of answer
#self.children.append(n1) # eventually I might want to keep around the entire search tree. I messed this up though
children.append(n1)
return children
def heuristic(self):
# a basic heuristic of taking the ones it seems pretty sure about
return min([(min(1 - v.value, v.value) , i, v) for i, v in enumerate(self.bool_vars)])[2]
def bbsolve(self):
root = self
res = root.buildProblem().solve()
heap = [(res, next(counter), root)]
bestres = 1e20 # a big arbitrary initial best objective value
bestnode = root # initialize bestnode to the root
print(heap)
nodecount = 0
while len(heap) > 0:
nodecount += 1 # for statistics
print("Heap Size: ", len(heap))
_, _, node = heappop(heap)
prob = node.buildProblem()
res = prob.solve()
print("Result: ", res)
if prob.status not in ["infeasible", "unbounded"]:
if res > bestres - 1e-3: #even the relaxed problem sucks. forget about this branch then
print("Relaxed Problem Stinks. Killing this branch.")
pass
elif node.is_integral(): #if a valid solution then this is the new best
print("New Best Integral solution.")
bestres = res
bestnode = node
else: #otherwise, we're unsure if this branch holds promise. Maybe it can't actually achieve this lower bound. So branch into it
new_nodes = node.branch()
for new_node in new_nodes:
heappush(heap, (res, next(counter), new_node ) ) # using counter to avoid possible comparisons between nodes. It tie breaks
print("Nodes searched: ", nodecount)
return bestres, bestnode
# a simple knapsack problem. we'll want to minimize the total cost of having each of these items, with different sizes.
# Use a random problem instance
N = 20
prices = -np.random.rand(N)
sizes = np.random.rand(N)
print(prices)
x = cvx.Variable(N)
constraints = []
constraints += [x <= 1, 0 <= x] #The relaxation of the binary variable constraint
constraints += [sizes*x <= 5] # total size of knapsack is 5
objective = prices * x
bool_vars = {x[i] for i in range(N)}
root = BBTreeNode(constraints = constraints, objective= objective, bool_vars = bool_vars)
res, sol = root.bbsolve()
print(sorted(list([(v.name(), v.value) for v in sol.bool_vars] + [(v.name(), v.value) for v in sol.vars] ) ))
# For comparison let's do the same problem using a built in mixed integer solver.
x = cvx.Variable(N, boolean=True)
constraints = []
constraints += [x <= 1, 0 <= x]
constraints += [sizes*x <= 5]
objective = prices * x
prob = cvx.Problem(cvx.Minimize(objective),constraints)
prob.solve()
print(x.value)
```

This is at least solving the problem fairly quickly. It needs better heuristics and to be sped up, which is possible in lots of ways. I was not trying to avoid all performance optimizations. It takes maybe 5 seconds, whereas the cvxpy solver is almost instantaneous.

```
Nodes searched: 67
[('var0[0]', 0.9999999958228145), ('var0[10]', -1.2718338055950193e-08), ('var0[11]', -1.3726395012104872e-08), ('var0[12]', 0.9999999982326986), ('var0[13]', 0.9999999973744331), ('var0[14]', 0.9999999988156902), ('var0[15]', -1.1908085711772973e-08), ('var0[16]', 0.9999999903780872), ('var0[17]', 0.9999999863334883), ('var0[18]', -1.1481655920777931e-08), ('var0[19]', 0.9999999996667646), ('var0[1]', 0.9999999969549299), ('var0[2]', 0.9999999979596141), ('var0[3]', -9.282428548104736e-09), ('var0[4]', -1.1378022795740783e-08), ('var0[5]', 0.9999999868240312), ('var0[6]', 0.9999999995068807), ('var0[7]', 0.9999999995399617), ('var0[8]', 0.9999999859520627), ('var0[9]', 0.9999999948062767)]
[ 1.00000000e+00 1.00000000e+00 1.00000000e+00 -1.44435650e-12
-1.88491321e-12 1.00000000e+00 1.00000000e+00 1.00000000e+00
1.00000000e+00 1.00000000e+00 -7.11338729e-13 1.99240081e-13
1.00000000e+00 1.00000000e+00 1.00000000e+00 -1.48697107e-12
1.00000000e+00 1.00000000e+00 -1.75111698e-12 1.00000000e+00]
```

Edit : I should investigate the Parameter functionality of cvxpy. That would make a make faster version than the one above based on deepcopy. If you made the upper and lower vectors on the binary variables parameters, you could restrict the interval to 0/1 without rebuilding the problem or copying everything.

```
#rough sketch
b = cvx.Variable(N)
u = cvx.Parameter(N)
u.value = np.ones(N)
l = cvx.Parameter(N)
l.value = np.zeros(N)
constraints += [b <= u, l <= b]
# change l.value and u.value in search loop.
```