Transportation cost optimisation using OMPR for a large data set

Question

I am solving a transport optimization problem given a set of constraints. The following are the three key data sets that I have

#demand file demand - has demand(DEMAND) across 4821(DPP) sale points(D)

head(demand)
                       D        PP DEMAND                             DPP
1 ADILABAD (V) - T:11001  OPC:PACK 131.00 ADILABAD (V) - T:11001:OPC:PACK
2 ADILABAD (V) - T:13003  OPC:PACK 235.00 ADILABAD (V) - T:13003:OPC:PACK
3  ADILABAD (V) - T:2006  PPC:PACK  30.00  ADILABAD (V) - T:2006:PPC:PACK
4  ADILABAD (V) - T:4001  OPC:PACK  30.00  ADILABAD (V) - T:4001:OPC:PACK
5  ADILABAD (V) - T:7006 OPC:NPACK  34.84 ADILABAD (V) - T:7006:OPC:NPACK
6         AHMEDABAD:1001  OPC:PACK 442.10         AHMEDABAD:1001:OPC:PACK

#Capacity file cc - has capacity constraint (MaxP, MinP) across 1823 sources(SOURCE)

head(cc,4)
                     SOURCE MinP  MaxP
1 CHILAMKUR:P:OPC:NPACK:0:R  900 10806
2 CHILAMKUR:P:OPC:NPACK:0:W  900 10806
3  CHILAMKUR:P:OPC:PACK:0:R 5628 67536
4  CHILAMKUR:P:OPC:PACK:0:W 5628 67536

#LandingCost file LCMat - This is a matrix with the landing cost to deliver the product across the demand location (DPP) from a given source(SOURCE). This is an 1823 x 4821 matrix. Since the landing costs to all locations do not exist from a given source, I have replace that with a huge cost (10^6) to such DPPs.

I am using the OMPR package in R to optimize shipping material to meet the demand. This is potentially a very simple transport problem but it is taking a lot of time. I am using a 16GB ram machine

The following is the code. Could anyone guide me on what I should do better?

a = Sys.time()
grid = expand.grid(i = 1:nrow(LCMat),j = 1:ncol(LCMat))
grid_solve = grid[which(LCMat < 10^6),]
grid_notsolve = grid[which(LCMat >= 10^6),]

model <- MILPModel() %>% 
  add_variable(x[grid$i, grid$j],lb = 0, type = "continuous") %>%
  add_constraint(x[grid_notsolve$i, grid_notsolve$j] == 0) %>%
  add_constraint(sum_over(x[i,j], i = 1:nrow(LCMat)) <= demand$DEMAND[j], j = 1:ncol(LCMat)) %>%
  add_constraint(sum_over(x[i,j], j = 1:ncol(LCMat)) <= cc$MaxP[i], i = 1:nrow(LCMat)) %>%
  add_constraint(sum_over(x[i,j], j = 1:ncol(LCMat)) >= cc$MinP[i], i = 1:nrow(LCMat)) %>%
  set_objective(sum_expr(LCMat[grid_solve$i,grid_solve$j]*x[grid_solve$i,grid_solve$j]),"min")

solution = model %>% solve_model(with_ROI(solver = "glpk", verbose = TRUE))
Sys.time() - a

Answer 1

Two options to potentially speed things up:

1. Make sure you use the latest CRAN versions of ompr and listcomp.
2. Try to use filter conditions to only create/use variables that are relevant to the model, instead of adding all nrow(LCMat)*ncol(LCMat) variables and then setting (potentially) a lot of them to 0. See the code below for an example. Depending on how sparse your problem is that could help as well.

The following code takes a sparse matrix (i.e. a matrix with many 0 elements or 10^6 elements in your case) and only generates x[i,j] variables that have an entry in sparse_matrix which is greater than 0. It hopefully illustrates how to use that feature and apply it to your case.

library(ompr)
sparse_matrix <- matrix(
  c(
    1, 0, 0, 1,
    0, 1, 0, 1,
    0, 0, 0, 1,
    1, 0, 0, 0
  ), byrow = TRUE, ncol = 4
)
is_connected <- function(i, j) {
  sparse_matrix[i, j] > 0
}
n <- nrow(sparse_matrix)
m <- ncol(sparse_matrix)
model <- MIPModel() |> 
  add_variable(x[i, j], i = 1:n, j = 1:m, is_connected(i, j)) |> 
  set_objective(sum_over(x[i, j], i = 1:n, j = 1:m, is_connected(i, j))) |> 
  add_constraint(sum_over(x[i, j], i = 1:n, is_connected(i, j)) <= 1, j = 1:m)

variable_keys(model)
#> [1] "x[1,1]" "x[1,4]" "x[2,2]" "x[2,4]" "x[3,4]" "x[4,1]"

extract_constraints(model)
#> $matrix
#> 3 x 6 sparse Matrix of class "dgCMatrix"
#>                 
#> [1,] 1 . . . . 1
#> [2,] . . 1 . . .
#> [3,] . 1 . 1 1 .
#> 
#> $sense
#> [1] "<=" "<=" "<="
#> 
#> $rhs
#> [1] 1 1 1

^{Created on 2022-03-12 by the reprex package (v2.0.1)}

Answer 2

• Both OMPR and GLPK are slow for large models.
• You are duplicating sum_over(x[i,j], j = 1:ncol(LCMat)). That leads to more nonzero elements than needed. I usually try to prevent that (even at the expense of more variables).