-
-
Save njtierney/f5e003bda327899f4eb621b26e0f2da5 to your computer and use it in GitHub Desktop.
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. |
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.
In the meantime the API has changed and I simplified it a bit.
sum_exp
is now calledsum_expr
and you can now write your (in)equalities directly as an R expression: e.g. (a <= b)