Model Formulation: Gurobi v. OMPR

[sbmack]

Gurobi has some sample models formulated in R available on its web page. I reformulated this Gurobi model
facility.R using OMPR in order to do a simple compare and contrast between the two modeling platforms . Below is a portion of the Gurobi model followed by the OMPR equivalent using the same data.

A visual comparison of both formulations shows how much simpler the OMPR version is and also how much it resembles legacy algebraic model management systems. OMPR is not yet fully featured, but its potential as a model management platform given its efficient and intuitive design paradigm seems to be significant. Note that Gurobi could probably consider using an augmented OMPR version as its formulation engine in R.

# Copyright 2021, Gurobi Optimization, LLC
#
# Facility location: a company currently ships its product from 5 plants
# to 4 warehouses. It is considering closing some plants to reduce
# costs. What plant(s) should the company close, in order to minimize
# transportation and fixed costs?
#
# Note that this example uses lists instead of dictionaries.  Since
# it does not work with sparse data, lists are a reasonable option.
#
# Based on an example from Frontline Systems:
#   http://www.solver.com/disfacility.htm
# Used with permission.

library(Matrix)
library(gurobi)

# define primitive data
nPlants     <- 5
nWarehouses <- 4
# Warehouse demand in thousands of units
Demand      <- c(15, 18, 14, 20)
# Plant capacity in thousands of units
Capacity    <- c(20, 22, 17, 19, 18)
# Fixed costs for each plant
FixedCosts  <- c( 12000, 15000, 17000, 13000, 16000)
# Transportation costs per thousand units
TransCosts  <- c(4000, 2000, 3000, 2500, 4500,
                 2500, 2600, 3400, 3000, 4000,
                 1200, 1800, 2600, 4100, 3000,
                 2200, 2600, 3100, 3700, 3200 )

flowidx <- function(w, p) {nPlants * (w-1) + p}

# Build model
model <- list()
model$modelname <- 'facility'
model$modelsense <- 'min'

# initialize data for variables
model$lb       <- 0
model$ub       <- c(rep(1, nPlants),   rep(Inf, nPlants * nWarehouses))
model$vtype    <- c(rep('B', nPlants), rep('C', nPlants * nWarehouses))
model$obj      <- c(FixedCosts, TransCosts)
model$varnames <- c(paste0(rep('Open',nPlants),1:nPlants),
                    sprintf('Trans%d,%d',
                            c(mapply(rep,1:nWarehouses,nPlants)),
                            1:nPlants))

# build production constraint matrix
A1 <- spMatrix(nPlants, nPlants, i = c(1:nPlants), j = (1:nPlants), x = -Capacity)
A2 <- spMatrix(nPlants, nPlants * nWarehouses,
               i = c(mapply(rep, 1:nPlants, nWarehouses)),
               j = mapply(flowidx,1:nWarehouses,c(mapply(rep,1:nPlants,nWarehouses))),
               x = rep(1, nWarehouses * nPlants))
A3 <- spMatrix(nWarehouses, nPlants)
A4 <- spMatrix(nWarehouses, nPlants * nWarehouses,
               i = c(mapply(rep, 1:nWarehouses, nPlants)),
               j = mapply(flowidx,c(mapply(rep,1:nWarehouses,nPlants)),1:nPlants),
               x = rep(1, nPlants * nWarehouses))
model$A           <- rbind(cbind(A1, A2), cbind(A3, A4))
model$rhs         <- c(rep(0, nPlants),   Demand)
model$sense       <- c(rep('<', nPlants), rep('=', nWarehouses))
model$constrnames <- c(sprintf('Capacity%d',1:nPlants),
                       sprintf('Demand%d',1:nWarehouses))

# Save model
gurobi_write(model,'facilityR.lp')

This is the OMPR equivalent:

tCosts <- matrix(TransCosts, nrow = 5) # tCosts is the TransCosts vector converted to a matrix

model <- MIPModel() %>% 
  add_variable(plants[i], i = 1:nPlants, type = "binary") %>%
  add_variable(ship[i, j], i = 1:nPlants, j = 1:nWarehouses, lb = 0) %>% 
  
  
  set_objective(sum_over(FixedCosts[i] * plants[i], i = 1:nPlants) +
                  sum_over(tCosts[i, j] * ship[i, j], i = 1:nPlants, j = 1:nWarehouses), "min") %>%

  
  add_constraint(sum_over(ship[i, j], j = 1:nWarehouses) <= Capacity[i] * plants[i], i = 1:nPlants) %>% 
  add_constraint(sum_over(ship[i, j], i = 1:nPlants) >= Demand[j], j = 1:nWarehouses)

[datadrivensupplychain]

Agreed, I recently purchased gurobi for a very large MILP project I'm working. Installing ROI.plugin.gurobi was a bit of a hassle, but the ability to model algebraically in R then pass to ompr is so much easier.

[dirkschumacher]

You can also use ompr without ROI by using ompr to formulate the problem and then using the API function to create the inputs for a custom solver interface (e.g. Gurobi). See this tweet for an example.

[datadrivensupplychain]

I just tried using a MILPModel and the gurobi package MILP formulation. No luck, so I'll stay with ROI and ROI.plugin.gurobi

[sbmack]

Whatever works. But I'm curious, is the fail in model formulation using MILPModel or in solver execution with a gurobi call from MILPModel?

As I noted above, generating constraints using matrix/array data structures can be confusing with MILPModel. While the updated, (development) version of MIPModel retains the simpler syntax of the original. So supersedes MILPModel.

That said, if your failure mode is in formulation suggest refactoring with the development version of MIPModel and trying that. Good Luck...

[datadrivensupplychain]

I’ve formulated a MILPModel using the CRAN version of ompr. (Very large optimization models, so using MILPModel is so much faster than MIPModel, at least with the CRAN version). I cannot extract the exact arguments from my MILPModel object, to then pass directly to the gurobi package. Right now I’m on the critical path with my project, but I’ll dig into using the dev version of MIPModel once the project closes in a few weeks. From: SteveM ***@***.***> Sent: Monday, May 02, 2022 12:25 PM To: dirkschumacher/ompr ***@***.***> Cc: datadrivensupplychain ***@***.***>; Comment ***@***.***> Subject: Re: [dirkschumacher/ompr] Model Formulation: Gurobi v. OMPR (Discussion #421) Whatever works. But I'm curious, is the fail in model formulation using MILPModel or in solver execution with a gurobi call from MILPModel? As I noted above, generating constraints using matrix/array data structures can be confusing with MILPModel. While the updated, (development) version of MIPModel retains the simpler syntax of the original. So supersedes MILPModel. That said, if your failure mode is in formulation suggest refactoring with the development version of MIPModel and trying that. Good Luck... — Reply to this email directly, view it on GitHub <#421 (reply in thread)> , or unsubscribe <https://github.com/notifications/unsubscribe-auth/ASKSE64WN6YMUI4C3TSXQH3VIAFXZANCNFSM5TBAN4HQ> . You are receiving this because you commented. <https://github.com/notifications/beacon/ASKSE66Y2EDH5DZ5DVWDPGTVIAFXZA5CNFSM5TBAN4H2YY3PNVWWK3TUL52HS4DFWFCGS43DOVZXG2LPNZBW63LNMVXHJKTDN5WW2ZLOORPWSZGOAAUNE5A.gif> Message ID: ***@***.*** ***@***.***> >

[sbmack]

Another obvious question is if your MILPModel is generating your model input matrix correctly. A crude way to examine a formulated model is to write out it out to as text file:

wd <- getwd()
file_out <- paste(wd, "test.txt", sep = "/")
ROI_write(as_ROI_model(model), file_out, "lp_cplex")

Unlike commercial algebraic modeling systems ROI_write does not use the actually variable names. But if you use reduced index sets to generate only a few constraints of each type you can visually examine them in the text output to see you are getting what you want.

[datadrivensupplychain]

@sbmack Do you know if it's possible to use a Gurobi solution and restart the optimization from a current solution? Sometimes I'll run a model for X hours just to see what will result; but then later want to use the solution as a starting point to further run the optimization.

[sbmack]

I haven't used Gurobi. But Gurobi documentation demonstrates the functionality for doing a warm start:
https://www.gurobi.com/documentation/9.5/examples/mip_starts.html

So if you can set that control parameter in your ROI.plugin.gurobi function call, it should work.

[harryprince]

I am working with ompr.roi and gurobi, when the number of variables in my problem goes to 1851360+, the adding the constraints will be very slow, does anyone can help me?

mdl <- MILPModel() %>%
   add_variable(x[i, k], i = 1:n, k = 1:m, type = "binary") %>% # binary selection for each variable
   add_variable(z[i, j, k], i = 1:n,j = 1:n, k = 1:m, type = "binary") %>% # binary selection for each variable pair
   set_objective(sum_over(x[i, k] * ord_w[i, k] , i = 1:n, k = 1:m), "max") %>% # more order 
   add_constraint(sum_over(con_wsp[i, j] * x[j, k] ,  j = 1:n) - x[i,k] >= 0,  i = 1:n, k = 1:m) %>%  # connectivity constraint
   add_constraint(sum_over(z[i, j, k], i = 1:n,j = 1:n) >= sum_over(x[i,k] + x[j,k] - 1,i=1:n,j=1:n),k = 1:m) %>%  # make zijk = xik * xjk
   add_constraint(sum_over(z[i, j, k], i = 1:n,j = 1:n) <= sum_over(x[i,k] ,i=1:n),k = 1:m) %>% 
   add_constraint(sum_over(z[i, j, k], i = 1:n,j = 1:n) <= sum_over(x[j,k] ,j=1:n),k = 1:m)

[sbmack]

I am assuming that your constraint matrix is sparse or very sparse with only a small percentage of non-zero coefficients in the constraint matrix.

An obvious tactic is to pre-process your data into sparse format and then pass that to ompr. That should significantly reduce the execution time of model creation. I played around with a sparse formulation but unfortunately have not fully developed a completed example. If you check the documentation of the ROI package: https://cran.r-project.org/web/packages/ROI/ROI.pdf you can see how he does it. Although I haven't tried to map his sparse construct to ompr.

Also note that MILPModel has been superseded by updates to MIPModel. So MIPModel is again the default version of the platform to use. The latest updates to the development version can be installed from github:

remotes::install_github("dirkschumacher/ompr")
remotes::install_github("dirkschumacher/ompr.roi")