nick-example-ompr.R

library(dplyr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)

rb_mat <- function(r,c,prob) matrix(rbinom(r*c,1,prob),r,c)

my_A <- rb_mat(r = 1000, c = 200) # 1000 heart attacks, 200 potential locations

# number of OHCA incidents
J <- nrow(my_A)
I <- ncol(my_A)
# number of AEDs
N <- 5

MIPModel() %>%
    add_variable(x[j], j = 1:J, type = "binary") %>%
    add_variable(y[i], i = 1:I, type = "binary") %>%
    # add_variable(a[i,j], i = 1:I, j = 1:J, type = "binary") %>%
    add_variable(my_A[i,j], i = 1:I, j = 1:J, type = "binary") %>%
    set_objective(sum_exp(x[j], j = 1:J), "max") %>%
    add_constraint(sum_exp(y[i], i = 1:I), "==", N) %>%
    add_constraint(sum_exp(a[i,j] * y[i], i = 1:I, j = 1:J), ">=", x[j])

# > Error in on_element(push, inplace_update_ast, get_ast_value, element) : 
# > The expression contains a variable, that is not part of the model.

@njtierney
Variable a[i,j] is not part of the model.
However if you multiply these two together (a and y) then you will have a quadratic constraint. Ompr cannot handle that at the moment.

Here is a working example:

library(dplyr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)

rb_mat <- function(r,c,prob) matrix(rbinom(r*c,1,prob),r,c)

a <- rb_mat(r = 20, c = 20, prob = 0.3) # 1000 heart attacks, 200 potential locations

# number of OHCA incidents
J <- nrow(a)
I <- ncol(a)
# number of AEDs
N <- 5

model <- MIPModel() %>%
  add_variable(x[j], j = 1:J, type = "binary") %>%
  add_variable(y[i], i = 1:I, type = "binary") %>%
  set_objective(sum_exp(x[j], j = 1:J), "max") %>%
  add_constraint(sum_exp(y[i], i = 1:I), "==", N) %>%
  # this creates a constraint for each j
  # also you can use normal R variables as model parameters
  add_constraint(x[j], "<=", sum_exp(a[i,j] * y[i], i = 1:I), j = 1:J)
model
result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))

get_solution(result, x[i])
get_solution(result, y[i])

In addition, ompr might still be a bit slow. If you want to see where the package is going take a look at JumP. An excellent modeling package in Julia. https://jump.readthedocs.io/en/latest/

Also, since it seems you are associated with an academic institution take a look at gurobi.com. This is currently the best commercial solver and it is free for academic use as far as I know. And at least in the operations research community it is generally accepted to use commercial solvers in papers.

Another option is to use CPLEX. Also free for academic use and you can use that with omprtogether with ROI.

Thanks for this Dirk!

This reads so much better than me munging together a bunch of matrices!

I've read up on gurobi and it seems great, I'll see how I go with using open source methods and then if speed is a serious issue I'll look into things like gurobi. Getting cplex is a bit tricky for me at this stage, unfortunately. I'd be open to the idea of using Julia as well, but I'm not sure how well julia talks with R these days.

OK, so the above works great, but it only seems to work when the A matrix is square - I get an error if the A matrix is rectangular.

If we consider something where we have 20 heart attack locations, and 10 possible locations, we have a 20x10 matrix.

a <- rb_mat(r = 20, # 20 heart attacks
            c = 10, # 10 possible locations
            prob = 0.3) # only 30% values are 1

# number of OHCA incidents
J <- nrow(a)
I <- ncol(a)
# number of AEDs
N <- 2

model <- MIPModel() %>%
    add_variable(x[j], j = 1:J, type = "binary") %>%
    add_variable(y[i], i = 1:I, type = "binary") %>%
    set_objective(sum_exp(x[j], j = 1:J), "max") %>%
    add_constraint(sum_exp(y[i], i = 1:I), "==", N) %>%
    # this creates a constraint for each j
    # also you can use normal R variables as model parameters
    add_constraint(x[j], "<=", sum_exp(a[i,j] * y[i], i = 1:I), j = 1:J)

model

result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))

#> Error in acc[[key]] : no such index at level 1

I'm not sure if this is because I have not specified something properly, or if there is something internally in ompr going on?

Hi @dirkschumacher,

I'm looking to replace my current code in copertura with ompr so that I can give users the capability to select other solvers, however trying the code that previously worked above, I get an error regarding a missing constraint relation.

library(dplyr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)

rb_mat <- function(r,c,prob) matrix(rbinom(r*c,1,prob),r,c)

a <- rb_mat(r = 20, c = 20, prob = 0.3) # 1000 heart attacks, 200 potential locations

# number of OHCA incidents
J <- nrow(a)
I <- ncol(a)
# number of AEDs
N <- 5

model <- MIPModel() %>%
    add_variable(x[j], j = 1:J, type = "binary") %>%
    add_variable(y[i], i = 1:I, type = "binary") %>%
    set_objective(sum_exp(x[j], j = 1:J), "max") %>%
    add_constraint(sum_exp(y[i], i = 1:I), "==", N) %>%
    # this creates a constraint for each j
    # also you can use normal R variables as model parameters
    add_constraint(x[j], "<=", sum_exp(a[i,j] * y[i], i = 1:I), j = 1:J)
#> Error in add_constraint_.optimization_model(model, lazyeval::lazy(constraint_expr),  : Does not recognize constraint expr. Missing the constraint relation

I was just wondering if you had any insights into this? I tried specifying the constraint relation, here, but then I get an error about the constraint not being well formed.

model <- MIPModel() %>%
    add_variable(x[j], j = 1:J, type = "binary") %>%
    add_variable(y[i], i = 1:I, type = "binary") %>%
    set_objective(sum_exp(x[j], j = 1:J), "max") %>%
    add_constraint(sum_exp(y[i], i = 1:I), constraint_expr = "==", N) %>%
    # this creates a constraint for each j
    # also you can use normal R variables as model parameters
    add_constraint(x[j], "<=", sum_exp(a[i,j] * y[i], i = 1:I), j = 1:J)
Error in add_constraint_.optimization_model(model, lazyeval::lazy(constraint_expr),  : 
  constraint not well formed. Must be a linear (in)equality.

Sorry if I've missed something basic here, this isn't my strong suit!

In the meantime the API has changed and I simplified it a bit. sum_exp is now called sum_expr and you can now write your (in)equalities directly as an R expression: e.g. (a <= b)

rb_mat <- function(r,c,prob) matrix(rbinom(r*c,1,prob),r,c)

a <- rb_mat(r = 20, c = 20, prob = 0.3) # 1000 heart attacks, 200 potential locations

# number of OHCA incidents
J <- nrow(a)
I <- ncol(a)
# number of AEDs
N <- 5

model <- MIPModel() %>%
  add_variable(x[j], j = 1:J, type = "binary") %>%
  add_variable(y[i], i = 1:I, type = "binary") %>%
  set_objective(sum_expr(x[j], j = 1:J), "max") %>%
  add_constraint(sum_expr(y[i], i = 1:I) == N) %>%
  add_constraint(x[j] <= sum_expr(a[i,j] * y[i], i = 1:I), j = 1:J)

Hi @dirkschumacher,

Just to extend from this, I'm using some "real data" now, and trying out ompr, I get a message warning of a 6 hour wait, perhaps I am doing something wrong here?

# testing out ompr to solve the optimum covering location problem

library(maxcovr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

york_existing <- york %>% filter(grade == "I")
york_proposed <- york %>% filter(grade != "I")
crime_user <- york_crime

existing_facility_cpp <- york_existing %>%
    select(lat,long) %>%
    as.matrix()

user_cpp <- crime_user %>%
    select(lat,long) %>%
    as.matrix()

dat_nearest_dist <-
    maxcovr::nearest_facility_dist(facility = existing_facility_cpp,
                                   user = user_cpp)

# make nearest dist into dataframe
# leave only those not covered
dat_nearest_no_cov <- dat_nearest_dist %>%
    dplyr::as_data_frame() %>%
    rename(user_id = V1,
           facility_id = V2,
           distance = V3) %>%
    filter(distance > 100) # 100m would be distance_cutoff

# give user an index
user <- crime_user %>% mutate(user_id = 1:n())

# join them, to create the "not covered" set of data
user_not_covered <- dat_nearest_no_cov %>%
    left_join(user,
              by = "user_id")

# # this takes the original user list, the full set of crime, or ohcas, etc
# existing_user <- user
#
# # update user to be the new users, those who are not covered
# user <- user_not_covered

# } # end NULL

proposed_facility_cpp <- york_proposed %>%
    select(lat, long) %>%
    as.matrix()

user_cpp <- user_not_covered %>%
    select(lat, long) %>%
    as.matrix()

a <- maxcovr::binary_matrix_cpp(facility = proposed_facility_cpp,
                                user = user_cpp,
                                distance_cutoff = 100)
# number of OHCA incidents
J <- nrow(a)
I <- ncol(a)
# number of AEDs
N <- 5

library(ROI)
#> ROI: R Optimization Infrastructure
#> Registered solver plugins: nlminb, glpk.
#> Default solver: auto.
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)

model <- MIPModel() %>%
    add_variable(x[j], j = 1:J, type = "binary") %>%
    add_variable(y[i], i = 1:I, type = "binary") %>%
    set_objective(sum_expr(x[j], j = 1:J), "max") %>%
    add_constraint(sum_expr(y[i], i = 1:I) == N) %>%
    add_constraint(x[j] <= sum_expr(a[i,j] * y[i], i = 1:I), j = 1:J)

#> |                    |  0% ~6 h remaining

Hey @njtierney, what are the values of J and I?

ompr is currently not optimised for speed and thus comes with a quite high cost of abstraction. I have some ideas on how to make ompr as fast as building a matrix, but at the moment, it is quite slow.

For maxcovr your could use ROI for example and build the coefficient matrix directly.

Hi @dirkschumacher,

I'm using some "real data" and trying to use ompr to solve my simple optimization problem, I get a message warning of "Error in eval(x$expr, x$env, emptyenv()) : object 'i' not found". Can you pls help me what is wrong. Initially it did run but had many warnings. Now for some reason, I'm getting this message.

Hotspots<- Ranking[1:1453, ] # Ranking is a dataset with estimated risk values
CMF<- data.frame(id= 1:3, cmf.name= c("CMF.Passive_FLB","CMF.Passive_FLBG","CMF.FLB_FLBG"), cmf.value= c(0.39, 0.297, 0.647))
cost<- data.frame(id= 1:3, c.name= c("c.Passive_FLB","c.Passive_FLBG","c.FLB_FLBG"), c.value= c(30000, 180000, 150000))

library(ompr)
library(magrittr)

#p means crash risk after the treatment application
p<- function(i, j) {
crash.before.treatment <- Hotspots[i, ] #from /hotspots dataset
CMFj <- CMF[j, ]
output<- crash.before.treatment$WeightedEstimate*CMFj$cmf.value
return(output)
}
#e.g.
p(1, 1)

#For excluding those site-countermeasure pairs not feasible for upgrade
#It is based on the logic that:
#If crossing type is passive, treatment options are: FLB, FLBG
#If crossing type is FLB, treatment option is FLBG
#If crossing type is FLBG, no option for upgrade

y<- function(i, j){
sitei <- Hotspots[i, ] #get from hotspots dataset
CrossingType<- sitei$ProtectionType

  if (j == 1 && CrossingType == "SRCS") {
      output<- 1
      } else {
      output<- 0
      }
  
 if (j == 2 && CrossingType == "SRCS") {
    output<- 1
  } else {
    output<- 0
  }
    
  if (j == 3 && CrossingType == "FLB") {
    output<- 1
  } else {
    output<- 0
  }
  
  return(output)

}

y(1,2) # testing

cost_var<- function(j){
cj <- cost[j, ] #get from cost dataset
costj<- cj$c.value
return (costj)
}
cost_var(3)# test

n<- 1453
m<- 3
b<- 1000000

model <- MIPModel() %>%

  add_variable(x[i, j], i = 1:n, j= 1:m, type = "binary") %>%
  
  set_objective(sum_expr(p(i,j) * x[i,j], i = 1:n, j= 1:m), "max") %>%
  
  add_constraint(sum_expr(cost_var(j) * x[i,j], i = 1:n, j= 1:m)<= b) %>%
  add_constraint(sum_expr(x[i,j], j = 1:m)<= 1, i= 1:n) %>% # at most one countermeasure used
  add_constraint(x[i,j] <= y(i,j), i= 1:n, j = 1:m) # only selective countermeasures can be implemented

model

#Thanks in advance for your feedback.

Hi @dirkschumacher,
Im trying to solve a problem on optimization .We have an list of 68 raw materials in inventory stock in a vector called Inv stock .And we have the required amount of items in an array qty_req. We have a matrix of raw materials which has list of raw materials as column and rows contains the quantity of raw material req to complete a unit of item per order. That means complete raw materials required to complete each order can be determined by multiplying the qty_req vector with the raw material matrix .

Problem: We want to find qtysold as a output which should give us the appropriate qty to be manufactured instead of qty req so that we complete atleast 80% of each order and the max revenue can be obtained. We have revenue/unit for each order .

Eg of Raw material matrix is as below:

OrderNo.	A	B	C
1	0	1	1
2	2	3	0
3	1	2	3

Eg. of req_qty vector

OrderNo.	Qty req
1	2
2	1
3	3

Eg. of Inv_stock

A	B	C
4	0	4

Eg of Revenue per unit

OrderNo.	Revenue/unit
1	2400
2	4500
3	450

Total material required for manufacturing the req quantity will become req_qty* Raw_material.

Please find below the code which is not working. Please provide your suggestions .
model <- MIPModel() %>%
add_variable(Qtysold[i],i=1:n,type = "integer" ,lb=0 ) %>%
add_variable(Total_mat[i,j],i=1:n,j=1:68, type = "continuous", lb = 0) %>%
set_objective(sum_expr(Qtysold[i]*Revenueperunit[i],i=1:n), "max")%>%
add_constraint(Qtysold[i] <= qty_req[i], i = 1:n)%>%
add_constraint(Qtysold[i]*Raw_mat==Total_mat[i,j],i=1:n,j=1:68)%>%
add_constraint(sum_expr(Total_mat[i,j],i=1:n)<=Inv_Stock[1,j],j=1:68)
result <- solve_model(model, with_ROI(solver = "glpk"))

On adding the last two lines im getting errors like unexpected operator '('.

Thanks in advance.

Hi @dirkschumacher, @njtierney,
I am solving an optimization problem using the followoing code.
library(Rglpk)
library(magrittr)
library(dplyr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)
data<-read.csv('input.csv',header=TRUE,sep=",")
Inv_Stock<-read.csv('INV_DATA.csv',header = TRUE,sep = ",")
Inv_Stock<- sapply(Inv_Stock,as.matrix)
qty_req<- data[,5]
qty_req<- sapply(qty_req,as.double)

Revenueperunit<-data[,4]
Revenueperunit<-sapply(Revenueperunit,as.double)

Raw_mat<-data[,-c(1,2,3,4,5)]
Raw_mat <- sapply(Raw_mat,as.matrix)

model <- MIPModel() %>%
add_variable(Qtysold[i],i=1:n,type = "integer" ,lb=0 ) %>%
add_variable(Total_mat[i,j],i=1:n,j=1:68, type = "continuous", lb = 0) %>%
set_objective(sum_expr(Qtysold[i]*Revenueperunit[i],i=1:n), "max")%>%
add_constraint(Qtysold[i] <= qty_req[i], i = 1:n)%>%
add_constraint(rawmat[j,i]*Qtysold[i]==Total_mat[j,i])%>%
add_constraint(sum_expr(Total_mat[i,j],i=1:n)<=Inv_Stock[1,j],j=1:68)
result <- solve_model(model, with_ROI(solver = "glpk"))

Context of the problem:
A vector holding inventory of total raw materials available by colums containing the raw material name and a single row with total qty available per column.
A matrix raw material with column names of raw material and each row indicating per order number, the quantity if raw material req for unit item to be manufactured.
A vector qty req indicating per order number the number of quantity asked by a customer to manufacture.
A vector revenue/unit indicating the revenue from each unit of material per order.
Problem statement:
We want to maximize total revenue
we want to determine qtysold such that atleast 80% of the qty req per order is delivered.
when we multiply qtysold with req raw material matrix , total req raw material columnwise should not exceed the total available in inventory stock.

Error: while adding last two constraints, its not giving answer, it says'(' unexpected operator found.
Error in acc[[key]] : no such index at level 1

Hi, I am now facing issue in only last line, can anyone please tell if this is making it a quadratic constraint.

model <- MIPModel() %>%
add_variable(Qtysold[i],i=1:n,type = "integer" ,lb=0 ) %>%
add_variable(Total_mat[i,j],i=1:n,j=1:68, type = "continuous", lb = 0) %>%
set_objective(sum_expr(Qtysold[i]*Revenueperunit[i],i=1:n), "max")%>%
add_constraint(Qtysold[i] <= qty_req[i], i = 1:n)%>%
add_constraint(Qtysold[i]*Raw_mat[i,j]==Total_mat[i,j],i=1:n,j=1:68)%>%
add_constraint(colSums(Total_mat[i,j],i=1:n)<=Inv_Stock[j],j=1:68)

Error:rror in check_for_unknown_vars_impl(model, the_ast) :
The expression contains a variable that is not part of the model.