A model in CPLEX format

Let’s define a linear programming model01 in CPLEX format. The bounds are not really necessary here, but are set to illustrate how to enter them in the model.

model01 <- 
'Maximize
  cost: 3x1 + 1x2 + 3x3
Subject To
  r1: - x1 + 2x2 + x3 <= 4
  r2: 4x2 - 3x3 <= 2
  r3: x1 - 3x2 + 2x3 <= 3
Bounds
  0 <= x1 <= 99999
  0 <= x2 <= 99999
Integer
  x1 
Binaries
  x3
End'

Let’s obtain the solution for this model using solve_lp:

solve_lp(model01, type = "CPLEX_LP")

## $sol
##   x1   x2   x3 
## 4.00 1.25 1.00 
## 
## $obj
## [1] 16.25
## 
## $success
## [1] TRUE

Coding the model in MathProg

Let’s write the same model in MathProg format. The syntax is similar to CPLEX, although variable types and bounds are specified when defining variables.

model01b <- 
'var x1, >=0, <=99999, integer;
var x2, >=0, <=99999;
var x3, binary;
maximize obj: 3 * x1 + x2 + 3 * x3;
s.t. r1: - x1 + 2 * x2 + x3 <= 4;
s.t. r2: 4 * x2 - 3 * x3 <= 2;
s.t. r3: x1 - 3 * x2 + 2 * x3 <= 3;
end;'

The solution of this model is:

solve_lp(model01b, type = "MathProg")

## $sol
##   x1   x2   x3 
## 4.00 1.25 1.00 
## 
## $obj
## [1] 16.25
## 
## $success
## [1] TRUE

A formulation in MathProg

Mathprog also allows writing not only linear programs as formulations of a generic problem that can be applied to different instances. Let’s see the formulation of the transportation problem, implemented in a small instance:

transport <- 'set I;
/* canning plants */

set J;
/* markets */

param a{i in I};
/* capacity of plant i in cases */

param b{j in J};
/* demand at market j in cases */

param d{i in I, j in J};
/* distance in thousands of miles */

param f;
/* freight in dollars per case per thousand miles */

param c{i in I, j in J} := f * d[i,j] / 1000;
/* transport cost in thousands of dollars per case */

var x{i in I, j in J} >= 0;
/* shipment quantities in cases */

minimize cost: sum{i in I, j in J} c[i,j] * x[i,j];
/* total transportation costs in thousands of dollars */

s.t. supply{i in I}: sum{j in J} x[i,j] <= a[i];
/* observe supply limit at plant i */

s.t. demand{j in J}: sum{i in I} x[i,j] >= b[j];
/* satisfy demand at market j */

data;

set I := Seattle San-Diego;

set J := New-York Chicago Topeka;

param a := Seattle     350
           San-Diego   600;

param b := New-York    325
           Chicago     300
           Topeka      275;

param d :              New-York   Chicago   Topeka :=
           Seattle     2.5        1.7       1.8
           San-Diego   2.5        1.8       1.4  ;

param f := 90;

end;'

Note how do we enter the generic formulation in the first section, and the specific values of the instance in the data section.

The solution for this instance is:

solve_lp(transport, type = "MathProg")

## $sol
##   x[Seattle,New-York]    x[Seattle,Chicago]     x[Seattle,Topeka] 
##                    50                   300                     0 
## x[San-Diego,New-York]  x[San-Diego,Chicago]   x[San-Diego,Topeka] 
##                   275                     0                   275 
## 
## $obj
## [1] 153.675
## 
## $success
## [1] TRUE

The formulation of the transportation problem above uses sets to define the parameters and variables of the model. We can also use formulations with subscripts, like in this formulation of a production planning with acquisition and storage costs:

prod_plan <- 'param n;

param sini;

param prodcosts{1..n};
param demand{1..n};
param storagecosts{1..n};
param prodmax{1..n};

var q{i in 1..n}, >=0, <=prodmax[i];
var s{i in 0..n}, >=0;

minimize cost: sum{i in 1..n} prodcosts[i] * q[i] + sum{i in 1..n} storagecosts[i] * s[i];

subject to ist: s[0] = sini;
subject to dem{i in 1..n}: s[i-1] + q[i] - s[i] = demand[i];

data;

param n:= 4;

param sini := 10;

param prodcosts := 1 100
                   2 150
                   3 200
                   4 250;

param demand :=  1 300
                2 350
                3 400
                4 150;

param storagecosts := 1 50
                      2 50
                      3 50
                      4 50;

param prodmax := 1 350
                 2 350 
                 3 350
                 4 350;

end;
'

solve_lp(prod_plan, type = "MathProg")

## $sol
## q[1] q[2] q[3] q[4] s[1] s[2] s[3] s[4] s[0] 
##  340  350  350  150   50   50    0    0   10 
## 
## $obj
## [1] 199000
## 
## $success
## [1] TRUE

Entering model parameters from the R environment

Rglpk can be seen as application programming interface (PI) to use GLPK through R. This means that the linear programming is solved using the C libraries of GLPK, returning the solution to the R environment. Instead of a CPLEX or Mathprog input, we can define the inputs for the API as R objects directly.

Let’s do that with the transportation problem of the previous section. I will define a tp_lp function for solving the transportation problem entering data from a named list with elements:

• supply or maximum capacity at each origin.
• demand from each destination.
• a costs matrix with cost coefficients.

Let’s store the parameters of the model built in MathProg format in a named list with this structure:

tp_sample <- list(supply = c(350, 600),
                  demand = c(325, 300, 275),
                  costs = matrix(c(2.5, 1.7, 1.8, 2.5, 1.8, 1.4)*90/1000, nrow = 2, ncol = 3, byrow = TRUE))

The function for solving the transportation problem with this input is:

lp_tp <- function(tp){
  
  m <- length(tp$supply)
  n <- length(tp$demand)
  
  var_names <- character(n*m)
  
  for(i in 1:m)
    for(j in 1:n)
      var_names[(i-1)*n + j] <- paste0("x_", i, "_", j)
  
  obj <- c(t(tp$costs))
  rhs <- c(tp$supply, tp$demand)
  dir <- c(rep("<=", m), rep(">=", n))
  const <- matrix(0, n+m, n*m)
  
  for(i in 1:m){
    for(j in 1:n){
      const[i, n*(i-1) + j] <- 1
      const[m+j, n*(i-1) + j] <- 1
    }
  }
  
  lp_solution <- Rglpk_solve_LP(obj = obj, 
                                mat = const,
                                dir = dir,
                                rhs = rhs,
                                max = FALSE)
  
  if(lp_solution$status == 0){
    sol <- lp_solution$solution
    names(sol) <- var_names
    obj = lp_solution$optimum
    return(list(sol = sol, obj = obj, success = TRUE))
  }else
    return(success = FALSE)
  
}

Let’s apply the function to our data, to check if the results is the same as in the previous section:

lp_tp(tp_sample)

## $sol
## x_1_1 x_1_2 x_1_3 x_2_1 x_2_2 x_2_3 
##    50   300     0   275     0   275 
## 
## $obj
## [1] 153.675
## 
## $success
## [1] TRUE

The obtained solution is the same as the previous section, but now with generic variable names generated inside lp_tp.

Specifying linear programming models in R

The Rglpk package allows specifying linear programming models in a variety of formats using the Rglpk_read_file function. The output of this function is a list of R objects that can be used as input to the Rglpk_solve_LP function of the same package, or to other solvers.

The MathProg language for defining linear programming models seems the best choice for Rglpk. It allows the straight modeling of the CPLEX format, and separating the generic formulation from the specific values of an instance of the problem. Another alternative is defining the elements of the model as vectors and matrices, and use them as input of Rglpk_solve_LP.

Acknowledgements

Thanks to Oriol Lordan for the model.mod trick, and to Vicenc Fernandez for the first example.

References

• GLPK Wikibook https://en.wikibooks.org/wiki/GLPK
• GLPK (GNU Linear Programming Kit) https://www.gnu.org/software/glpk/
• Sallan, J. M.; Lordan, O. & Fernandez, V. (2016). Modeling and solving linear programming with R. Omnia Books. https://www.omniascience.com/books/index.php/scholar/catalog/book/34

Session info

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_ES.UTF-8        LC_COLLATE=es_ES.UTF-8    
##  [5] LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=es_ES.UTF-8   
##  [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] Rglpk_0.6-4 slam_0.1-49
## 
## loaded via a namespace (and not attached):
##  [1] bookdown_0.24     digest_0.6.27     R6_2.5.0          jsonlite_1.7.2   
##  [5] magrittr_2.0.1    evaluate_0.14     blogdown_1.5      stringi_1.7.3    
##  [9] rlang_0.4.12      jquerylib_0.1.4   bslib_0.2.5.1     rmarkdown_2.9    
## [13] tools_4.1.2       stringr_1.4.0     xfun_0.23         yaml_2.2.1       
## [17] compiler_4.1.2    htmltools_0.5.1.1 knitr_1.33        sass_0.4.0