3
$\begingroup$

I am familiar with how to use ompr::MIPModel but I am trying to learn how to use MILPModel to take advantage of the model build speed. A simplified version of my model is below. I have two decision variables, x and y, binary and of equal length. I have constraints on the sum of all the x decision variables, and the sum of all the y decision variables. So far so good with MILPModel, I can build the model and solve it fast.

The problem is when I try to use the next constraint. The LHS of this constraint multiplies the x binary decision variables by a numeric column in a dataframe of the same length, then multiplies that by a matrix where the rows are equal to the length of x. Similar story in the RHS with the y variable. I then iterate this constraint 20 times to represent all the columns of the matrix.

I've used constraints similar to this many times using MIPModel, but now when I try this I get an error message, non-numeric argument to binary operator. I assume this has something to do with the colwise function, but I am completely lost on how to approach this, even after reading up on the ompr github site. Thanks in advance for any help.

add_variable(x[i], i=1:10, type='binary') %>%
add_variable(y[i], i=1:10, type='binary') %>%
add_constraint(sum_expr(x[i],i=1:10) <= 5) %>%
add_constraint(sum_expr(y[i],i=1:10) <= 3) %>%

#model builds and solves until this point...
add_constraint( 
sum_expr( x[i]* df$numeric_column[i] * matrix_a[i,j],i=1:10) <= 
sum_expr(  2* y[i]* df$numeric_column[i] * df$other_numeric_column[i] * matrix_a[i,j],i=1:10), 
j=1:20) 
$\endgroup$
  • $\begingroup$ I can confirm the error message but don't know the fix, so I contacted Dirk Schumacher, the author of OMPR. He said he'll take a look next week. $\endgroup$ – prubin Aug 9 at 13:16
  • $\begingroup$ Answered here: stackoverflow.com/questions/63294305/… $\endgroup$ – Ralph Asher Aug 9 at 13:53
  • $\begingroup$ @prubin I answered in my crosspost in SO. There were some gymnastics involved in adapting the matrix multiplication in the last constraint, to allow for vectorization. Also - I learned about OMPR from your blog a few years ago, thank you so much for what you've done for the OR community! $\endgroup$ – Ralph Asher Aug 9 at 13:54
1
$\begingroup$

Adding a follow-up on how to style matrix multiplication-type constraints and objective function values. This has been a major pain point for me, and now that I have a template, this will be a real productivity boost as I switch to being "MILPModel native".

Replace SYMPHONY with the solver of your choice...


library(tidyverse)
library(magrittr)
library(ompr)
library(ompr.roi)
library(ROI.plugin.symphony)


rm(list=ls())



set.seed(42)

#in this example, mat1 is a numeric matrix representing coefficients for 40 binary decision variables
mat1 <- matrix(ncol=10,nrow=4,runif(400))

#define this function, it will be necessary for matrix multiplication inside a MILPModel
matrix_multiplication_fcn <- function(static_matrix, row_variable, column_variable){
  vapply(seq_along(row_variable), function(k) static_matrix[row_variable[k], column_variable[k]], numeric(1L))  }

milp_model <- ompr::MILPModel() %>%
  add_variable(assign_units[rowindex,colindex], rowindex=1:4,colindex=1:10,type='binary') %>%
  #total binaries ==10
  add_constraint(sum_expr( assign_units[rowindex,colindex],rowindex=1:4,colindex=1:10 )==10 ) %>%
  
  #sum of binaries * mat1 <= 7
  add_constraint( sum_expr( ompr::colwise(
    matrix_multiplication_fcn(static_matrix=mat1,row_variable=rowindex,column_variable=colindex)) *
                              assign_units[rowindex,colindex], 
                            rowindex = 1:4, colindex = 1:10) <= 7) %>%

  #objective: maximize value
  set_objective(sum_expr(
    ompr::colwise(matrix_multiplication_fcn(static_matrix=mat1,row_variable=rowindex,column_variable=colindex)) * 
      assign_units[rowindex,colindex],
    rowindex=1:4,colindex=  1:10),sense='max')
    
milp_model_out <-  milp_model %>%   
  ompr::solve_model(with_ROI(solver = "symphony",verbosity=-2,gap_limit=0,time_limit=180, node_limit=-1,first_feasible=FALSE))


#same with MIPModel

mip_model <- ompr::MIPModel() %>%
  add_variable(assign_units[rowindex,colindex], rowindex=1:4,colindex=1:10,type='binary') %>%
  #total binaries ==10
  add_constraint(sum_expr( assign_units[rowindex,colindex],rowindex=1:4,colindex=1:10 )==10 ) %>%
  
  #sum of binaries * mat1 <= 7
  add_constraint( sum_expr( mat1[rowindex,colindex] * assign_units[rowindex,colindex], 
                            rowindex = 1:4, colindex = 1:10) <= 7) %>%
  
  #objective: maximize value
  set_objective(sum_expr(mat1[rowindex,colindex] * assign_units[rowindex,colindex],
    rowindex=1:4,colindex=  1:10),sense='max')

mip_model_out <-  mip_model %>%   
  ompr::solve_model(with_ROI(solver = "symphony",verbosity=-2,gap_limit=0,time_limit=180, node_limit=-1,first_feasible=FALSE))

#compare results, both should be slightly above 7 (thus slightly violating the constraint), I assume due to floating point issues
mip_model_out
milp_model_out
```
| improve this answer | |
$\endgroup$
0
$\begingroup$

Self-answered by OP from Stack Overflow, which I'll be posting here on their behalf as a community wiki to get this question off the unanswered list.


Figured it out. To use matrix algebra in a constraint requires a little bit of acrobatics. Good luck figuring out how to use matrix algebra in the objective function, should you need to.

Example comparing MIPModel to MILPModel is below.

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

rm(list=ls())

numvec1 <- runif(10)
numvec2 <- runif(10)
matrix_a <- matrix(nrow=10,ncol=20,data=runif(10*20))

my_mip_model <- MIPModel() %>%
 
 add_variable(x[i], i=1:10, type='binary') %>%
 add_variable(y[i], i=1:10, type='binary') %>%
 add_constraint(sum_expr(numvec1[i]*x[i],i=1:10) <= 5) %>%
 add_constraint(sum_expr(2*y[i],i=1:10) <= 3) %>%
 
 add_constraint( 
   sum_expr( x[i]* numvec1[i] * matrix_a[i,j],i=1:10) <= 
     sum_expr(  2* y[i]* numvec1[i] * numvec2[i] * matrix_a[i,j],i=1:10), 
   j=1:20) %>%
 set_objective( sum_expr(3*x[i]*numvec1[i],i=1:10),sense='max')

my_mip_model_solve <- my_mip_model %>% solve_model(with_ROI(solver='glpk'))



#functionally equivalent using MILPmodel----

my_milp_model <- MILPModel() %>%

 add_variable(x[i], i=1:10, type='binary') %>%
 add_variable(y[i], i=1:10, type='binary') %>%
 add_constraint(sum_expr( colwise(numvec1[i]) * x[i],i=1:10) <= 5) %>%
 add_constraint(sum_expr( colwise(2) * y[i],i=1:10) <= 3)  %>%
 set_objective(sum_expr( colwise(3*numvec1[i]) * x[i],i=1:10),sense='max')

#now to add the matrix constraints, add a loop on the matrix column index j.
#with MIPModel we could just iterate on j in a single constraint, but here it appears
#we need to add the same constraint multiple times, and use the value of j to
#calculate the indices in as.numeric(matrix_a) that we want to use.
for(j in 1:ncol(matrix_a)){
 
 my_milp_model %<>% add_constraint(
   
   sum_expr( x[i]* colwise(numvec1[i] *
     as.numeric(matrix_a)[(i + (nrow(matrix_a)*j -nrow(matrix_a)))]),i=1:10) <= 
       sum_expr(  y[i]* colwise(2* numvec1[i] * numvec2[i] * 
   as.numeric(matrix_a)[(i + (nrow(matrix_a)*j -nrow(matrix_a)))]) ,i=1:10) ) 
   
}

my_milp_model_solve <- my_milp_model %>% solve_model(with_ROI(solver='glpk'))

#objective value and results should be equal...
my_mip_model_solve
my_milp_model_solve

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.