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”.

# 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<Integer>($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<Rational>
- The Matrix of inequalites as an Matrix<Rational>
- The Objective function as a list of Rationals.

sub gomory { my $p = new Polytope<Rational>; my $A = shift; my $ineqs = shift; my @obj_list = @_; my $obj = new Vector<Rational>(@obj_list); $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<Rational>(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<Rational>(@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<Rational>(@new_row); my $new_col_vec = new Vector<Rational>(@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; push @new_ineq_row, 1; my $new_ineq_row_vec = new Vector<Rational>(@new_ineq_row); my $new_ineq_col_vec = new Vector<Rational>(@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<Rational>(($A|$new_col_vec) / $new_row_vec); my $new_ineqs = new Matrix<Rational>(($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.

my $A = new Matrix<Rational>(<<"."); -14 7 -2 1 0 0 -3 0 1 0 1 0 -3 2 -2 0 0 1 . my $ineqs = new Matrix<Rational>(<<"."); 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 . my @objective = (0,4,-1,0,0,0); gomory($A, $ineqs, @objective);