Mixed Integer Quadratically Constrained Programming (MIQCP) with gurobi or alabama plugin

[dsiddons]

Dear Dirk,
Sorry if this is in the wrong section... I wasn't sure if one could class this as an issue...

I am trying to use rmpk to create (what I think is...) a MIQCP. I have tried to use the "alabama" plugin, but I fear that my understanding of its controls (start values for example) are the main reason why it will not produce a solution (I have limited the number of iterations to try to get it to spit out a result to no avail!).

I then got an academic licence for "gurobi" - supposed to be THE package for QP... and using this as the solver I get:

Error: model$Q is NOT matrix, sparseMatrix or simple_triplet_matrix

There is no slot model$Q...

How can I specify a quadratic constraint? If this is not possible, how could I import the matrices and objective function into the gurobi or the ROI r framework?

I attach the dataset and the code I am using in the hopes that anyone can help!
aceites_LP.xlsx

library(readxl)
library(tidyverse)
library(rmpk)
library(ROI)
library(ROI.plugin.glpk)
library(ROI.plugin.alabama)
library(ROI.plugin.gurobi)

aceites_LP <- data.frame(read_excel("data/aceites_LP.xlsx", 
                                      sheet = "ACEITES", skip = 1))
aceites_LP$Pk_j <- factor(aceites_LP$Pk_j)
aceites_LP$Pl_j <- factor(aceites_LP$Pl_j)
  
shelves_LP <- data.frame(read_excel("data/aceites_LP.xlsx", 
                                      sheet = "SHELF", skip = 1))
  
S <- nrow(shelves_LP) # number of shelves_LP
P <- nrow(aceites_LP) # number of products we are dealing with


nextshelf <- function(S) {
  mat1 <- combn(1:S,2)[,which(combn(1:S,2)[2,]-combn(1:S,2)[1,] != 1)]
  return(mat1)
}
n = 400
control_optim <- list(trace = 1, tmax = 1, maxit = 10L)

control <- list(start = rep(0,400), lam0 = 10, sig0 = 100,  
                tol = 1e-07, max_iter = 1, verbose = TRUE, NMinit = FALSE, 
                ilack.max = 6, i.scale = 1, e.scale = 1, kkt2.check = TRUE,
                control.optim = control_optim)

#{solver <- ROI_optimizer("alabama", control = control)# All variables classes as continuous for this

solver <- ROI_optimizer("gurobi")
model <- optimization_model(solver)

x <- model$add_variable("x", i = 1:S, j = 1:P, type = "binary", lb = 0, ub = 1)
f <- model$add_variable("f", i = 1:S, j = 1:P, type = "integer", lb = 0, ub = 10)


# Objective Maximise profit number of products on shelf
model$set_objective(sum_expr(f[i,j]*colwise(aceites_LP$Pu_j[j]),i = 1:S, j = 1:P), "max")
 
# constraints

# length of products cannot exceed shelf length
model$add_constraint(sum_expr(colwise(aceites_LP$Pw_j[j])*f[i,j], j = 1:P) <= shelves_LP$Sl_i[i], i = 1:S) 

# number of products cannot exceed stock
model$add_constraint(sum_expr(f[i,j], i = 1:S) <= aceites_LP$Ps_j[j], j = 1:P) 

# total products on all shelves within maximum and minimum values
model$add_constraint(sum_expr(f[i,j], i = 1:S) <= aceites_LP$Fmax_j[j], j = 1:P) 
model$add_constraint(sum_expr(f[i,j], i = 1:S) >= aceites_LP$Fmin_j[j], j = 1:P) 

# number of shelves used for products cannot exceed minimum and maximum number of shelves
model$add_constraint(sum_expr(x[i,j], i = 1:S) <= aceites_LP$smax_j[j], j = 1:P) 
model$add_constraint(sum_expr(x[i,j], i = 1:S) >= aceites_LP$smin_j[j], j = 1:P) 

# shelves selected for products have ot have a product on it...
model$add_constraint(f[i,j] - 100*x[i,j] <= 0, j = 1:P, i = 1:S) 
model$add_constraint(f[i,j] - x[i,j] >= 0, j = 1:P, i = 1:S)

# products have to be on contiguous shelves
model$add_constraint(x[a,j] + x[b,j] <= 1, j = colwise(1:P), a = nextshelf(S)[1,], b = nextshelf(S)[2,])  

# this constraint should make that a product on more than one shelf will have the same number of products on each shelf
# this is the quadratic constraint causing issues :)
model$add_constraint(f[i,j]*sum_expr(x[i,j], i = 1:S) - sum_expr(f[i,j], i = 1:S) <= 0, i = 1:S, j = 1:P)
  
}
model$optimize()

model$get_variable_value(f[i,j]) %>% 
    arrange(i) %>% 
    pivot_wider(names_from = i,
              values_from = value) %>% 
  as.data.frame()

[dirkschumacher]

Hi David,

it should work with quadratic constraints, but I have never tried it with Gurobi. If you can reduce your problem a bit and make it reproducible without having to load data from an excel file, I am happy to look into it. See also the reprex package for more guidence.

Without having tested it I see three options:

• There is a bug in rmpk which causes a problem (I started moving towards a MOI inspired solver interface which might have broken something. Then the pandemic hit and I was occupied with other stuff :))
• Alabama has a problem with that specific model
• Gurobi's ROI plugin does not support quadratic constraints properly? I am not an academic though, so I don't have access to Gurobi. So difficult to test. But we could write a specific solver interface for Gurobi or improve the existing ROI plugin a bit.

Best,
Dirk

[dsiddons]

Dirk, Thanks for the swift reply! I shall pare down the script and use reprex. Worst case I will make the table in R. I am almost certainly doing something wrong with Alabama...

…

On Wed, 31 Mar 2021 at 09:49, Dirk Schumacher ***@***.***> wrote: Hi David, it should work with quadratic constraints, but I have never tried it with Gurobi. If you can reduce your problem a bit and make it reproducible without having to load data from an excel file, I am happy to look into it. See also the reprex <https://github.com/tidyverse/reprex> package for more guidence. Without having tested it I see three options: - There is a bug in rmpk which causes a problem (I started moving towards a MOI inspired solver interface which might have broken something. Then the pandemic hit and I was occupied with other stuff :)) - Alabama has a problem with that specific model - Gurobi's ROI plugin does not support quadratic constraints properly? I am not an academic though, so I don't have access to Gurobi. So difficult to test. But we could write a specific solver interface for Gurobi or improve the existing ROI plugin a bit. Best, Dirk — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#54 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACZYMBEVS2DXPOWIGPXMPEDTGMY7VANCNFSM42DEJDDQ> .

-- David C. Siddons Investigador Principal y Docente Escuela de Biología, Ecología y Gestión Universidad del Azuay Teléfono (593)74091000 Ext (458) http://www.uazuay.edu.ec

--

______________________________              Universidad del Azuay           (http://www.uazuay.edu.ec <http://www.uazuay.edu.ec>)

[dsiddons]

Dear Dirk,

Here is what is hopefully a completely self-sufficient code example:

library(tidyverse)
library(ompr) ## needed for "colwise"
library(rmpk)
#> 
#> Attaching package: 'rmpk'
#> The following objects are masked from 'package:ompr':
#> 
#>     MIPModel, sum_expr
library(ROI)
#> ROI: R Optimization Infrastructure
#> Registered solver plugins: nlminb, alabama, glpk, gurobi, ipop, nloptr.bobyqa, nloptr.crs2lm, nloptr.direct, nloptr.directL, nloptr.lbfgs, nloptr.neldermead, nloptr.newuoa, nloptr.sbplx, nloptr.stogo, nloptr.tnewton, nloptr.varmetric, nloptr.cobyla, nloptr.mma, nloptr.auglag, nloptr.isres, nloptr.slsqp, osqp, qpoases, quadprog, scs.
#> Default solver: auto.
library(ROI.plugin.alabama)
library(ROI.plugin.gurobi)

Pw_j <- c(14,10,12,9,7,8,8,8,8,7,7,7,9.5,9.5,10,9,7,7,7,9,6,6,7,6,8,7.5,7,6,6,5,5,5,15,15,15,14,15,4,8.3,7)
Pu_j <- c(2.1412,1.5395,0.8054,0.4917,0.4261,0.2767,0.5045,0.6213,0.4076,0.168,1.0817,2.1441,0.67,0.6078,0.909,1.09,0.3411,0.5237,0.2934,0.4846,0.3483,0.329,0.47,0.3242,0.2438,0.527,0.3119,1.5972,0.1977,1.2392,0.2472,0.4171,0.3265,0.5035,0.5694,0.2284,0.1932,0.21,1.8789,0.5757)
Fmin_j <- rep(1,40)
smin_j <- rep(1,40)
Fmax_j <- c(2,4,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,5)
smax_j <- rep(2,40)
Sl_i <- c(151,200,180,218)

P = 40
S = 4

nextshelf <- function(S) {
  mat1 <- combn(1:S,2)[,which(combn(1:S,2)[2,]-combn(1:S,2)[1,] != 1)]
  return(mat1)
}
n = 320

control_optim <- list(trace = 1, tmax = 1, maxit = 10L)
control <- list(start = rep(0,n), lam0 = 10, sig0 = 100,  
                tol = 1e-07, max_iter = 1, verbose = TRUE, NMinit = FALSE, 
                ilack.max = 6, i.scale = 1, e.scale = 1, kkt2.check = TRUE,
                control.optim = control_optim)

#{solver <- ROI_optimizer("alabama", control = control)# All variables classes as continuous for this
{solver <- ROI_optimizer("gurobi")
model <- optimization_model(solver)

x <- model$add_variable("x", i = 1:S, j = 1:P, type = "binary", lb = 0, ub = 1)
f <- model$add_variable("f", i = 1:S, j = 1:P, type = "integer", lb = 0, ub = 10)

# Objective Maximise profit number of products on shelf
model$set_objective(sum_expr(f[i,j]*colwise(Pu_j[j]),i = 1:S, j = 1:P), "max")
 
# constraints

# length of products cannot exceed shelf length
model$add_constraint(sum_expr(colwise(Pw_j[j])*f[i,j], j = 1:P) <= Sl_i[i], i = 1:S) 
# total products on all shelves within maximum and minimum values
model$add_constraint(sum_expr(f[i,j], i = 1:S) <= Fmax_j[j], j = 1:P) 
model$add_constraint(sum_expr(f[i,j], i = 1:S) >= Fmin_j[j], j = 1:P) 
# number of shelves used for products cannot exceed minimum and maximum number of shelves
model$add_constraint(sum_expr(x[i,j], i = 1:S) <= smax_j[j], j = 1:P) 
model$add_constraint(sum_expr(x[i,j], i = 1:S) >= smin_j[j], j = 1:P) 
# shelves selected for products have ot have a product on it...
model$add_constraint(f[i,j] - 100*x[i,j] <= 0, j = 1:P, i = 1:S) 
model$add_constraint(f[i,j] - x[i,j] >= 0, j = 1:P, i = 1:S)
# products have to be on contiguous shelves
model$add_constraint(x[a,j] + x[b,j] <= 1, j = colwise(1:P), a = nextshelf(S)[1,], b = nextshelf(S)[2,])  

# this constraint should make that a product on more than one shelf will have the same number of products
# on each shelf

# model$add_constraint(f[i,j]*sum_expr(x[i,j], i = 1:S) - sum_expr(f[i,j], i = 1:S) <= 0, i = 1:S, j = 1:P)
  
}
model$optimize()

model$get_variable_value(f[i,j]) %>% 
arrange(i)%>% 
  
  pivot_wider(names_from = i,
              values_from = value) %>% 
  as.data.frame()
#>    name  j 1 2 3 4
#> 1     f 33 1 0 0 0
#> 2     f 38 0 0 0 1
#> 3     f 27 0 1 0 0
#> 4     f 19 0 1 0 0
#> 5     f  8 1 1 0 0
#> 6     f 28 0 0 7 1
#> 7     f 11 3 5 0 0
#> 8     f 20 0 1 0 0
#> 9     f 21 0 0 1 0
#> 10    f  4 0 1 0 0
#> 11    f 16 0 8 0 0
#> 12    f 25 0 0 1 0
#> 13    f 32 0 0 1 0
#> 14    f 37 1 0 0 0
#> 15    f  6 1 0 0 0
#> 16    f 36 1 0 0 0
#> 17    f 15 0 0 0 7
#> 18    f  1 2 0 0 0
#> 19    f  7 0 0 1 0
#> 20    f 12 0 0 8 0
#> 21    f 22 0 0 1 0
#> 22    f 35 1 0 0 0
#> 23    f 39 0 0 3 5
#> 24    f 18 0 1 0 0
#> 25    f 34 0 1 0 0
#> 26    f  3 0 0 1 0
#> 27    f 10 1 0 0 0
#> 28    f 26 0 1 0 0
#> 29    f 17 0 1 0 0
#> 30    f 29 0 0 1 0
#> 31    f  9 1 0 0 0
#> 32    f 13 0 0 0 1
#> 33    f 14 0 1 0 0
#> 34    f 24 0 0 1 0
#> 35    f  5 1 0 0 0
#> 36    f 30 0 0 0 8
#> 37    f  2 0 0 0 4
#> 38    f 23 0 0 0 1
#> 39    f 31 1 0 0 0
#> 40    f 40 0 1 0 0

^{Created on 2021-04-05 by the reprex package (v1.0.0)}

If I add the quadratic constraint I get:

# this constraint should make that a product on more than one shelf will have the same number of products
# on each shelf

model$add_constraint(f[i,j]*sum_expr(x[i,j], i = 1:S) - sum_expr(f[i,j], i = 1:S) <= 0, i = 1:S, j = 1:P)
  
}
model$optimize()
#> Error: model$Q is NOT matrix, sparseMatrix or simple_triplet_matrix

and with alabama... It gets stuck somewhere and will not produce a solution

[dirkschumacher]

Weird. I cannot reproduce the error, but I also do not have Gurobi. Alabma is started without a modelling error.
PS: you do not need the colwise here.