This tutorial is probably also available as a Jupyter notebook in the `demo`

folder in the polymake source and on github.

Different versions of this tutorial: latest release, release 3.3, release 3.2, nightly master

This tutorial was written by Michael Schmitt, a student of the course “Discrete Optimization”.

# Branch & Bound

In this Tutorial, I like to describe, how one can implement the Branch & Bound Algorithm for polytopes using polymake.

## Preamble

This preamble must be integrated in every script file, so that we can use the polytope functionalities of polymake.

use application "polytope";

At first we create some global variables that keep track over the best integer solution, that has been discovered so far.

$::z = -10000; $::x_best;

## Subfunction

I introduce a subroutine, that returns the non-integer part of some rational number.

sub f { my $x=shift; my $denom = denominator($x); my $nom = convert_to<Integer>($denom * $x); my $remainder = $nom % $denom; if($remainder < 0){ $remainder += $denom; } return new Rational($remainder, $denom); }

## Main routine

It follows the main routine, that will call itself recursively. The parameters that are passed to this function are:

- The objective function as an Vector<Rational>
- A list of inequalities that define the currently examined polytope

sub branchnbound { my $p = new Array<Rational>;

The following line gets the parameters.

my $objective = shift; my @ineqs = @_; my $mat = new Matrix<Rational>(@ineqs); #print "Current Inequalities: \n", $mat, "\n"; $p->INEQUALITIES = $mat;

### Feasibility check

Let's check, whether the polytope is empty or not. If it is empty, we can stop examining this branch of the B&B Tree.

if(!$p->FEASIBLE){ print "EMPTY!\n"; return; }

### Solving the relaxation

We look for the optimal solution of the linear relaxation.

$p->LP = new LinearProgram<Rational>(LINEAR_OBJECTIVE=>$objective); my $x = $p->LP->MAXIMAL_VERTEX; my $ctx = $p->LP->MAXIMAL_VALUE; print "Relaxed solution: ", $x, " -> ", $ctx,"\n";

### Pruning

And if that is already worse then the best integer solution so far, we can prune this branch.

# Pruning if($ctx <= $::z){ "is pruned, because <= $::z\n"; return 0; }

### Check for integer solution

In the following loop, we find the first coordinate of the solution vector, that is not integer.

my $i = 0; for(my $k=1; $k<$x->dim(); $k++){ if(f($x->[$k]) != 0){ $i = $k; } }

If we didn't find any, we have found an integer solution and possibly a new best solution. In the case of an integer solution, we can stop here.

if($i == 0){ print "Solution is Integer.\n"; # Update the best solution if($ctx > $::z){ $::x_best = $x; $::z = $ctx; print "Solution is better then the last optimum.\n"; } return; }

### Branching

Now lets create some new inequalities. At first round the not integer coordinate of x up and down.

#Round up and down: my $x_i_down = $x->[$i] - f($x->[$i]); my $x_i_up = $x_i_down + 1;

To generate P_1, add the inequality x_i ⇐ x_i_down.

my @ineq_to_add_1 = ($x_i_down); for (my $k = 1; $k<$x->dim(); $k++){ if($k == $i){ push @ineq_to_add_1, -1; }else{ push @ineq_to_add_1, 0; } }

To generate P_2, add the inequality x_i >= x_i_up.

my @ineq_to_add_2 = (-$x_i_up); for (my $k = 1; $k<$x->dim(); $k++){ if($k == $i){ push @ineq_to_add_2, 1; }else{ push @ineq_to_add_2, 0; } }

Now add the new inequalities to the old ones …

@ineq_to_add_1 = new Array<Rational>(@ineq_to_add_1); @ineq_to_add_2 = new Array<Rational>(@ineq_to_add_2); my @new_ineqs_1 = @ineqs; push @new_ineqs_1, @ineq_to_add_1; my @new_ineqs_2 = @ineqs; push @new_ineqs_2, @ineq_to_add_2;

And do the recursive calls.

Here I added some printouts, that give an inorder traversal of the B&B tree and the edges that are used for that.

print "Set x$i <= $x_i_down\n"; branchnbound(@new_ineqs_1); print "Forget x$i <= $x_i_down\n"; print "Set x$i >= $x_i_up\n"; branchnbound(@new_ineqs_2); print "Forget x$i >= $x_i_up\n"; }

Finally, the first call of the Routine:

my @ineqs = ([14,-7,2],[3,0,-1],[3,-2,2],[0,1,0],[0,0,1]); my $objective = [0,4,-1]; branchnbound($objective, @ineqs); print "The solution is: $::x_best with the optimal value $::z\n";