This tutorial was written by Michael Schmitt, a student of the course “Discrete Optimization”.
In this Tutorial, I like to describe, how one can implement the Branch & Bound Algorithm for polytopes using polymake.
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.
> declare $z = -10000; > declare $x_best;
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); > }
It follows the main routine, that will call itself recursively. The parameters that are passed to this function are:
> sub branchnbound { > > > # The following line gets the parameters. > > my $objective = shift; > my @ineqs = @_; > > my $mat = new Matrix<Rational>(@ineqs); > #print "Current Inequalities: \n", $mat, "\n"; > > my $p = new Polytope<Rational>; > $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. > > if($ctx <= $z){ > print "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 ... > > 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($objective, @new_ineqs_1); > print "Forget x$i <= $x_i_down\n"; > > print "Set x$i >= $x_i_up\n"; > branchnbound($objective, @new_ineqs_2); > print "Forget x$i >= $x_i_up\n"; > > }
Finally, the first call of the Routine:
> @ineqs = ([14,-7,2],[3,0,-1],[3,-2,2],[0,1,0],[0,0,1]); > $objective = [0,4,-1]; > branchnbound($objective, @ineqs); > > print "The solution is: $x_best with the optimal value $z\n"; Relaxed solution: 1 20/7 3 -> 59/7 Set x1 <= 2 Relaxed solution: 1 2 1/2 -> 15/2 Set x2 <= 0 Relaxed solution: 1 3/2 0 -> 6 Set x1 <= 1 Relaxed solution: 1 1 0 -> 4 Solution is Integer. Solution is better then the last optimum. Forget x1 <= 1 Set x1 >= 2 EMPTY! Forget x1 >= 2 Forget x2 <= 0 Set x2 >= 1 Relaxed solution: 1 2 1 -> 7 Solution is Integer. Solution is better then the last optimum. Forget x2 >= 1 Forget x1 <= 2 Set x1 >= 3 EMPTY! Forget x1 >= 3 The solution is: 1 2 1 with the optimal value 7