Modelling techniques in OMPR

This vignettes discribes the modelling techniques available in ompr to make your life easier when developing a mixed integer programming model.

A MIP Model

You can think of a MIP Model as a big constraint maxtrix and a set of vectors. But you can also think of it as a set of decision variables, an objective function and a number of constraints as equations/inequalities. ompr implements the latter approach.

For example, Wikipedia describes the Knapsack problem like this: \[ \begin{equation*} \begin{array}{ll@{}ll} \text{max} & \displaystyle\sum\limits_{i=1}^{n} v_{i}x_{i} & &\\ \text{subject to}& \displaystyle\sum\limits_{i=1}^{n} w_{i}x_{i} \leq W, & &\\ & x_{i} \in \{0,1\}, &i=1 ,\ldots, n& \end{array} \end{equation*} \]

This is the ompr equivalent:

n <- 10; W <- 2
v <- runif(n);w <- runif(n)
model <- MIPModel() %>% 
  add_variable(x[i], i = 1:n, type = "binary") %>% 
  set_objective(sum_over(v[i] * x[i], i = 1:n)) %>% 
  add_constraint(sum_over(w[i] * x[i], i = 1:n) <= W)

The overall idea is to use modern R idioms to construct models like the one above as readable as possible directly in R. ompr will do the heavy lifting and transforms everything into matrices/vectors and pass it to your favorite solver.

Pipes

Each function in ompr creates immutable copies of the models. In addition the function interface has been designed to work with magrittr pipes. You always start with an empty model and add components to it.

MIPModel() %>% 
  add_variable(x) %>% 
  set_objective(x) %>% 
  add_constraint(x <= 1)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 1 
##   Integer: 0 
##   Binary: 0 
## Model sense: maximize 
## Constraints: 1

Variable types

Variables can be of type continuous, integer or binary.

MIPModel() %>% 
  add_variable(x, type = "integer") %>% 
  add_variable(y, type = "continuous") %>% 
  add_variable(z, type = "binary")

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 1 
##   Integer: 1 
##   Binary: 1 
## No objective function. 
## Constraints: 0

Variable bounds

Variables can have lower and upper bounds.

MIPModel() %>% 
  add_variable(x, lb = 10) %>% 
  add_variable(y, lb = 5, ub = 10)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 2 
##   Integer: 0 
##   Binary: 0 
## No objective function. 
## Constraints: 0

Indexed variables

Often when you develop a complex model you work with indexed variables. This is an important concept ompr supports.

MIPModel() %>% 
  add_variable(x[i], i = 1:10) %>%  # creates 10 decision variables
  set_objective(x[5]) %>% 
  add_constraint(x[5] <= 10)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 10 
##   Integer: 0 
##   Binary: 0 
## Model sense: maximize 
## Constraints: 1

Summation over variables

If you have indexed variables then you often want to sum over a subset of variables.

The following code creates a model with three decision variables \(x_1\), \(x_2\), \(x_3\). An objective function \(\sum_i x_i\) and one constraint \(\sum_i x_i \leq 10\).

MIPModel() %>% 
  add_variable(x[i], i = 1:3) %>% 
  set_objective(sum_over(x[i], i = 1:3)) %>% 
  add_constraint(sum_over(x[i], i = 1:3) <= 10)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 3 
##   Integer: 0 
##   Binary: 0 
## Model sense: maximize 
## Constraints: 1

Quantifiers

add_variable, add_constraint, set_bounds, sum_over all support a common quantifier interface that also supports filter expression. A more complex example will show what that means.

MIPModel() %>% 
  # Create x_{i, j} variables for all combinations of i and j where
  # i = 1:10 and j = 1:10.
  add_variable(x[i, j], type = "binary", i = 1:10, j = 1:10) %>% 
  
  # add a y_i variable for all i between 1 and 10 with i mod 2 = 0
  add_variable(y[i], type = "binary", i = 1:10, i %% 2 == 0) %>% 
  
  # we maximize all x_{i,j} where i = j + 1
  set_objective(sum_over(x[i, j], i = 1:10, j = 1:10, i == j + 1)) %>% 
  
  # for each i between 1 and 10 with i mod 2 = 0
  # we add a constraint \sum_j x_{i,j}
  add_constraint(sum_over(x[i, j], j = 1:10) <= 1, i = 1:10, i %% 2 == 0) %>% 
  
  # of course you can leave out filters or add more than 1
  add_constraint(sum_over(x[i, j], j = 1:10) <= 2, i = 1:10)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 0 
##   Integer: 0 
##   Binary: 105 
## Model sense: maximize 
## Constraints: 15

Special bounds on a subset of variables

Imagine you want to model a matching problem with a single binary decision variable \(x_{i,j}\) that is \(1\) iff object \(i\) is matched to object \(j\). One constraint would be to allow matches only if \(i \neq j\). This can be modeled by a constraint or by selectively changing bounds on variables. The latter approach can be used by solvers to improve the solution process.

MIPModel() %>% 
  add_variable(x[i, j], i = 1:10, j = 1:10, 
               type = "integer", lb = 0, ub = 1) %>% 
  set_objective(sum_over(x[i, j], i = 1:10, j = 1:10)) %>% 
  add_constraint(x[i, i] == 0, i = 1:10) %>% 
  
   # this sets the ub to 0 without adding new constraints
  set_bounds(x[i, i] <= 0, i = 1:10)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 0 
##   Integer: 100 
##   Binary: 0 
## Model sense: maximize 
## Constraints: 10

External model parameters

Of course you will need external parameters for your models. You can reuse any variable defined in your R environment within the MIP Model.

n <- 5 # number of our variables
costs <- rpois(n, lambda = 3) # a cost vector
max_elements <- 3
MIPModel() %>% 
  add_variable(x[i], type = "binary", i = 1:n) %>% 
  set_objective(sum_over(costs[i] * x[i], i = 1:n)) %>% 
  add_constraint(sum_over(x[i], i = 1:n) <= max_elements)

## Mixed integer linear optimization problem
## Variables:
##   Continuous: 0 
##   Integer: 0 
##   Binary: 5 
## Model sense: maximize 
## Constraints: 1

Extract model solutions

Once you have a model, you pass it to a solver and get back a solutions. The main interface to extract variable values from a solution is the function get_solution. It returns a data.frame for indexed variable and thus makes it easy to subsequently use the value.

We use ROI and GLPK to solve it.

library(ROI)

## ROI: R Optimization Infrastructure

## Registered solver plugins: nlminb, glpk.

## Default solver: auto.

library(ROI.plugin.glpk)
library(ompr.roi)

set.seed(1)
n <- 5
weights <- matrix(rpois(n * n, 5), ncol = n, nrow = n)
result <- MIPModel() %>% 
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>% 
  set_objective(sum_over(weights[i, j] * x[i, j], i = 1:n, j = 1:n)) %>% 
  add_constraint(sum_over(x[i, j], j = 1:n) == 1, i = 1:n) %>% 
  solve_model(with_ROI("glpk", verbose = TRUE))

## <SOLVER MSG>  ----
## GLPK Simplex Optimizer 5.0
## 5 rows, 25 columns, 25 non-zeros
##       0: obj =  -0.000000000e+00 inf =   5.000e+00 (5)
##       5: obj =   2.400000000e+01 inf =   0.000e+00 (0)
## *    14: obj =   4.400000000e+01 inf =   0.000e+00 (0)
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer 5.0
## 5 rows, 25 columns, 25 non-zeros
## 25 integer variables, all of which are binary
## Integer optimization begins...
## Long-step dual simplex will be used
## +    14: mip =     not found yet <=              +inf        (1; 0)
## +    14: >>>>>   4.400000000e+01 <=   4.400000000e+01   0.0% (1; 0)
## +    14: mip =   4.400000000e+01 <=     tree is empty   0.0% (0; 1)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----

get_solution(result, x[i, j]) %>% 
  dplyr::filter(value == 1)

##   variable i j value
## 1        x 4 1     1
## 2        x 2 2     1
## 3        x 5 3     1
## 4        x 3 4     1
## 5        x 1 5     1

You can also fix certain indexes.

get_solution(result, x[2, j])

##   variable j value
## 1        x 1     0
## 2        x 2     1
## 3        x 3     0
## 4        x 4     0
## 5        x 5     0

Feedback

Do you have any questions, ideas, comments? Or did you find a mistake? Let’s discuss on Github.

Dirk Schumacher

2017-04-17