{ompr} 1.0 release + mini roadmap

[dirkschumacher]

Hey everyone,

I worked hard in the recent weeks to publish a new version of {ompr} together with a new version of {ompr.roi} 🙌. {ompr.roi} will land on CRAN in the coming days, once the new {ompr} package has propagated through CRAN mirrors.

What is new?

The most important change is a redesign of the MIPModel backend, based on what I learned building rmpk. The new backend is now easier to maintain, faster for larger models, has fewer bugs (hopefully) and some new features. It should be backwards compatible. Should you notice any problem with your existing code, please let me know.

sum_over

sum_expr is now called sum_over. The reason is that sum_expr does not create an expression anymore as in the original MIPModel version, it rather sums over variables. In fact, it is now a regular function. sum_expr will stay in the package for now to not break anything.

sum_over(i, i = 1:10, i < 5)
#> [1] 10

set_bounds supports (in)equalities

Inspired by @sbmack's issue #365, MIPModel now supports (in)equalities for setting bounds which gives you more options to specify variable bounds.

library(ompr)
model <- MIPModel() |>
  add_variable(x[i], i = 1:2) |>
  set_bounds(x[i] <= i, i = 1:2) |>
  set_bounds(x[i] >= 0, i = 1:2)
variable_bounds(model)
#> $lower
#> [1] 0 0
#> 
#> $upper
#> [1] 1 2

All quantifiers are now managed by the package {listcomp}

I extracted the quantifier logic out of {rmpk} and built {listcomp}. This gives you now a consistent way of specifying quantifiers that support dependent indexes and {rlangs}' quotation syntax. It works the same for add_variable, add_constraint, set_bounds and sum_over.

library(ompr)
custom_indexes <- 10:15
bnd <- 15
MIPModel() |>
  add_variable(x[i, j], i = 1:5, j = custom_indexes[i], i + j < !!bnd) |> 
  variable_keys()
#> [1] "x[1,10]" "x[2,11]"

You can now access additional solver specific information from {ompr.roi}

{ROI} has two generic status codes 0 and 1. Where 0 means "success" and 1 otherwise. These will now be mapped to "success" or "error" in {ompr.roi}. However the solver itself reveals more information that can be accessed using additional_solver_output(solution) (available in v1.0.1). For example, when you set a time limit in GLPK, ROI returns code=1, but the actual solver code from GLPK is GLP_FEAS. You can now access that for a specific solver.

library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)
set.seed(42)
max_bins <- 10
bin_size <- 3
n <- 10
weights <- runif(n, max = bin_size)
model <- MIPModel() |>
  add_variable(y[i], i = 1:max_bins, type = "binary") |>
  add_variable(x[i, j], i = 1:max_bins, j = 1:n, type = "binary") |>
  set_objective(sum_expr(y[i], i = 1:max_bins), "min") |>
  add_constraint(sum_expr(weights[j] * x[i, j], j = 1:n) <= y[i] * bin_size, i = 1:max_bins) |>
  add_constraint(sum_expr(x[i, j], i = 1:max_bins) == 1, j = 1:n)
result <- solve_model(model, with_ROI("glpk",
  control = list(
    tm_limit = 1000,
    presolve = FALSE, verbose = TRUE
  )
))
#> <SOLVER MSG>  ----
#> GLPK Simplex Optimizer, v4.65
#> 20 rows, 110 columns, 210 non-zeros
#>       0: obj =   0.000000000e+00 inf =   1.000e+01 (10)
#>      32: obj =   6.362622072e+00 inf =   0.000e+00 (0)
#> *    36: obj =   6.362622072e+00 inf =   0.000e+00 (0)
#> OPTIMAL LP SOLUTION FOUND
#> GLPK Integer Optimizer, v4.65
#> 20 rows, 110 columns, 210 non-zeros
#> 110 integer variables, all of which are binary
#> Integer optimization begins...
#> Long-step dual simplex will be used
#> +    36: mip =     not found yet >=              -inf        (1; 0)
#> Solution found by heuristic: 8
#> + 65050: mip =   8.000000000e+00 >=   7.000000000e+00  12.5% (1399; 24849)
#> TIME LIMIT EXCEEDED; SEARCH TERMINATED
#> <!SOLVER MSG> ----
solver_status(result)
#> [1] "error"
additional_solver_output(result)$ROI$status
#> $code
#> [1] 1
#> 
#> $msg
#>   solver glpk
#>     code 2
#>   symbol GLP_FEAS
#>  message Solution is feasible.
#> roi_code 1

Take a look at the full release-notes for {ompr} and {ompr.roi} for all the details.

If you fnd any bugs, please report issues and post your question in the Discussion section of the repository.

What is next?

I do not have a big master plan, but here are some ideas for the future:

• With the new backend, it would be possible to extend {ompr} to quadratic programs (both in the objective function and the constraints).
• Although {ompr} allows different external solver packages, in practice only {ompr.roi} is really used. We could merge {ompr.roi} with {ompr} and offer {ROI} out of the box.
• MILPModel turned out to be too complicated for complex use-cases and it produces a lot of problems, as far as I can tell. As the new MIPModel is reasonable fast and could be made even faster, I strongly tend towards removing MILPModel from the package at some point in the future.
• {rmpk} will be separate for now, with some innovations coming back into {ompr}.

As I am developing all of that in my free time, I don't want to promise anything. But starting this year I went full-time freelance, which might give me more opportunities to work on {ompr} and related packages. You can follow me on Twitter, if you want to get regular updates on {ompr} et al. Or even work with me on optimization projects.