Linear programming consists of optimizing a linear function subject to a set of linear constraints. If all or some of the decision variables are integer, it is called integer and mixed integer linear programming, respectively. Linear programming is a powerful optimization technique, and the first to consider when tackling combinatorial problems.

In this post, I will show how to solve linear programming models in R using the `Rglpk`

package. This is a R API for **GLPK**, the **G**NU **L**inear **P**rogramming **K**it. It is an open-source set of linear programming routines written in ANSI C and organized in the form of a callable library.

Let’s see how can we solve a linear programming formulation written in CPLEX or MathProg formats. We need two steps for doing this:

- Reading the model from a text file with
`Rglpk_read_file`

. - Solving the model with
`Rglpk_solve_LP`

. The inputs of`Rglpk_solve_LP`

are the outputs of`Rglpk_read_file`

.

To make the task easier, I have written a `solve_lp`

function that does the following:

- Creates a text file with the content of
`model`

using`writeLines`

. - Reads the text file of
`type`

MathProg or CPLEX_LP with`Rglpk_read_file`

, and saves the result in the object`lp_model`

. - Obtains the solution
`lp_solution`

from the elements of`lp_model`

with`Rglpk_solve_LP`

. - Removes the text file with
`unlink`

. - If the solution is feasible, returns the solution and
`success = TRUE`

. If not, returns`success = FALSE`

. The solution is a list including a named vector`sol`

, the value of the objective function`obj`

and the`success`

variable.

```
solve_lp <- function(model, type = "MathProg"){
writeLines(model, "model.mod")
lp_model <- Rglpk_read_file("model.mod", type = type)
lp_solution <- Rglpk_solve_LP(obj = lp_model$objective,
mat = lp_model$constraints[[1]],
dir = lp_model$constraints[[2]],
rhs = lp_model$constraints[[3]],
types = lp_model$types,
bounds = lp_model$bounds,
max = lp_model$maximum)
unlink("model.mod")
if(lp_solution$status == 0){
sol <- lp_solution$solution
names(sol) <- attr(lp_model, "objective_vars_names")
obj = lp_solution$optimum
return(list(sol = sol, obj = obj, success = TRUE))
}else
return(success = FALSE)
}
```

## 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
```