This tutorial was written by Michael Schmitt, a student of the course "Discrete Optimization".
====== Gomory Cuts ======
In this Tutorial, I describe how to solve an ILP using the gomory cut approach implemented in polymake.
===== Preamble =====
This preamble must be integrated in every script file, so that we can use the polytope functionalities of polymake.
> use application "polytope";
===== 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($denom * $x);
> my $remainder = $nom % $denom;
> if($remainder < 0){
> $remainder += $denom;
> }
> return new Rational($remainder, $denom);
> }
===== Main Routine =====
The main Routine of this script gets invoked with the following parameters:
* The Matrix of equalities as an Matrix
* The Matrix of inequalites as an Matrix
* The Objective function as a list of Rationals.
> sub gomory {
>
> my $A = shift;
> my $ineqs = shift;
> my @obj_list = @_;
> my $obj = new Vector(@obj_list);
>
> my $p = new Polytope;
> $p->EQUATIONS = $A;
> $p->INEQUALITIES = $ineqs;
>
>
> ### Feasibility check
> # Since all the inequalities that we add are valid for the integer hull
> # of the polytope, we can conclude that the integer hull is empty,
> # when the examined polytope in any iteration cycle is empty.
>
> if(!$p->FEASIBLE){
> print "The integer hull of P is empty.\n";
> return;
> }
>
>
> ### Solving the relaxation
>
> $p->LP = new LinearProgram(LINEAR_OBJECTIVE=>$obj);
>
> my $x = $p->LP->MAXIMAL_VERTEX;
> print "Optimal solution of the relaxation: ", $x,"\n";
>
>
> my @B = ();
> my @N = ();
> for(my $k = 1; $k < $x->dim(); $k++){
> if($x->[$k] > 0){
> push @B, $k;
> }else{
> push @N, $k;
> }
> }
>
>
> ### Check for integer solution
> # If we found an integer solution in the linear relaxation,
> # it's the solution of the ILP aswell.
> # If not, we remember the coordinate that is not integer for later steps.
>
> my $i = -1;
> for(my $k=0; $k<@B; $k++){
> if(f($x->[$B[$k]]) != 0){
> $i = $k;
> }
> }
>
> if($i == -1){
> print "Found optimal integer solution: ", $x, " -> ", $p->LP->MAXIMAL_VALUE,"\n";
> return;
> }
>
>
> ### Compute optimal basis
> # In this example, I only implemented the computation of a Basis
> # in the non-degenerated case (that is, where every Basis Variables are strictly positive.
>
> my @A_b = ();
> foreach my $b (@B){
> push @A_b, $A->col($b);
> }
>
> my $A_b = transpose(new Matrix(@A_b));
> my $A_inv = inv($A_b);
>
>
>
> ### Compute new equality
> # Here is where the remembered non-integer variable comes into play,
> # to compute the Gomory cut itself.
>
> my @new_row;
> my @new_col;
>
> for(my $j=0; $j<$A->rows(); $j++){
> $new_col[$j] = 0;
> }
>
> $new_row[0] = f($A_inv->row($i)*($A->col(0)));
> foreach my $j (@N){
> $new_row[$j] = (-1)*f((-1)*$A_inv->row($i)*$A->col($j));
> }
> foreach my $j (@B){
> $new_row[$j] = 0;
> }
> push @new_row, 1;
> my $new_row_vec = new Vector(@new_row);
> my $new_col_vec = new Vector(@new_col);
>
> print "New equality: $new_row_vec\n";
>
>
> ### Compute new inequality
> # The new inequality demands, that the newly introduced slack variable will be positive.
>
> my @new_ineq_col;
> for(my $j=0; $j<$ineqs->rows(); $j++){
> $new_ineq_col[$j] = 0;
> }
>
> my @new_ineq_row = @new_ineq_col;
> #push @new_ineq_row, 0; #this leads to dimension mismatches
> push @new_ineq_row, 1;
>
> my $new_ineq_row_vec = new Vector(@new_ineq_row);
> my $new_ineq_col_vec = new Vector(@new_ineq_col);
>
>
> ### Construct the new matrices
> # The two matrices constructed are the equality matrix A and the inequality matrix
> # that restricts the variables to be positive.
> # In both cases we add a new column on the right side and a new row at the bottom.
>
> my $new_A = new Matrix(($A|$new_col_vec) / $new_row_vec);
> my $new_ineqs = new Matrix(($ineqs|$new_ineq_col_vec) / $new_ineq_row_vec);
>
> push @obj_list, 0;
>
> gomory($new_A, $new_ineqs, @obj_list);
> }
==== First call of the Function ====
Finally we create two example matrices and an example objective function, to call the main routine.
> $A = new Matrix([
> [-14, 7, -2, 1, 0, 0],
> [-3, 0, 1, 0, 1, 0],
> [-3, 2, -2, 0, 0, 1]]);
>
> $ineqs = new Matrix([
> [0, 1, 0, 0, 0, 0],
> [0, 0, 1, 0, 0, 0],
> [0, 0, 0, 1, 0, 0],
> [0, 0, 0, 0, 1, 0],
> [0, 0, 0, 0, 0, 1]]);
>
> @objective = (0,4,-1,0,0,0);
>
> gomory($A, $ineqs, @objective);
Optimal solution of the relaxation: 1 20/7 3 0 0 23/7
New equality: 5/7 0 0 -2/7 -4/7 0 1
Optimal solution of the relaxation: 1 5/2 7/4 0 5/4 3/2 0
New equality: 1/2 0 0 0 0 0 -1/2 1
Optimal solution of the relaxation: 1 2 1/2 1 5/2 0 1 0
New equality: 1/2 0 0 0 0 -1/2 0 0 1
Optimal solution of the relaxation: 1 2 1 2 2 1 1 0 0
Found optimal integer solution: 1 2 1 2 2 1 1 0 0 -> 7