Problem Description
From (1) (see references below):
I have 4 stores (1,2,3,4) and I can apply 3 treatments (A,B,C) to each of the 4 stores. Each treatment has its own cost and profit.
The matrix is as follows:
Store Treatment Cost Profit
1 A 50 100
1 B 100 200
1 C 75 50
2 A 25 25
2 B 150 0
2 C 50 25
3 A 100 300
3 B 125 250
3 C 75 275
4 A 25 25
4 B 50 75
4 C 75 125How can I maximize the profit having a constraint on maximum cost in R? In reality, I have a large number of stores and hence cannot be done manually. Each store can get only 1 treatment.
Thanks in advance.
Optimization Model
The optimization model for this problem can look like:
Loading the data
Loading the data into a data frame is easy with read.table:
df<-read.table(text="
Store Treatment Cost Profit
1 A 50 100
1 B 100 200
1 C 75 50
2 A 25 25
2 B 150 0
2 C 50 25
3 A 100 300
3 B 125 250
3 C 75 275
4 A 25 25
4 B 50 75
4 C 75 125
",header=T)
df## Store Treatment Cost Profit
## 1 1 A 50 100
## 2 1 B 100 200
## 3 1 C 75 50
## 4 2 A 25 25
## 5 2 B 150 0
## 6 2 C 50 25
## 7 3 A 100 300
## 8 3 B 125 250
## 9 3 C 75 275
## 10 4 A 25 25
## 11 4 B 50 75
## 12 4 C 75 125OMPR Model
I am going to try to model this problem with OMPR (2) and solve it with the Symphony MIP solver (3).
First we need a bunch of packages:
library(dplyr)
library(tidyr)
library(ROI)
library(ROI.plugin.symphony)
library(ompr)
library(ompr.roi)Extract data
OMPR is using numbers as indices so it makes sense to extract the cost and profit data as matrices. We do this as follows:
# extract some basic data
stores<-unique(df$Store)
stores## [1] 1 2 3 4treatments <- levels(df$Treatment)
treatments## [1] "A" "B" "C"num_treatments <- length(treatments)
# cost matrix
cost <- as.matrix(spread(subset(df,select=c(Store,Treatment,Cost)),Treatment,Cost)[,-1])
cost## A B C
## 1 50 100 75
## 2 25 150 50
## 3 100 125 75
## 4 25 50 75# profit matrix
profit <- as.matrix(spread(subset(df,select=c(Store,Treatment,Profit)),Treatment,Profit)[,-1])
profit## A B C
## 1 100 200 50
## 2 25 0 25
## 3 300 250 275
## 4 25 75 125# maximum cost allowed
max_cost <- 300MIP Model
The OMPR package allows us to specify linear equations algebraically. Most LP/MIP solvers in R use a matrix notation, which often is more difficult to use. Here we go:
m <- MIPModel() %>%
add_variable(x[i,j], i=stores, j=1:num_treatments, type="binary") %>%
add_constraint(sum_expr(x[i,j], j=1:num_treatments)<=1, i=stores) %>%
add_constraint(sum_expr(cost[i,j]*x[i,j], i=stores, j=1:num_treatments) <= max_cost) %>%
set_objective(sum_expr(profit[i,j]*x[i,j], i=stores, j=1:num_treatments),"max") %>%
solve_model(with_ROI(solver = "symphony",verbosity=1))##
## Starting Preprocessing...
## Preprocessing finished...
## with no modifications...
## Problem has
## 5 constraints
## 12 variables
## 24 nonzero coefficients
##
## Total Presolve Time: 0.000000...
##
## Solving...
##
## granularity set at 25.000000
## solving root lp relaxation
## The LP value is: -650.000 [0,5]
##
##
## ****** Found Better Feasible Solution !
## ****** Cost: -650.000000
##
##
## ****************************************************
## * Optimal Solution Found *
## * Now displaying stats and best solution found... *
## ****************************************************
##
## ====================== Misc Timing =========================
## Problem IO 0.000
## ======================= CP Timing ===========================
## Cut Pool 0.000
## ====================== LP/CG Timing =========================
## LP Solution Time 0.000
## LP Setup Time 0.000
## Variable Fixing 0.000
## Pricing 0.000
## Strong Branching 0.000
## Separation 0.000
## Primal Heuristics 0.000
## Communication 0.000
## Total User Time 0.000
## Total Wallclock Time 0.000
##
## ====================== Statistics =========================
## Number of created nodes : 1
## Number of analyzed nodes: 1
## Depth of tree: 0
## Size of the tree: 1
## Number of solutions found: 1
## Number of solutions in pool: 1
## Number of Chains: 1
## Number of Diving Halts: 0
## Number of cuts in cut pool: 0
## Lower Bound in Root: -650.000
##
## ======================= LP Solver =========================
## Number of times LP solver called: 1
## Number of calls from feasibility pump: 0
## Number of calls from strong branching: 0
## Number of solutions found by LP solve: 1
## Number of bounds changed by strong branching: 0
## Number of nodes pruned by strong branching: 0
## Number of bounds changed by branching presolver: 0
## Number of nodes pruned by branching presolver: 0
##
## ==================== Primal Heuristics ====================
## Time #Called #Solutions
## Rounding I 0.00
## Rounding II 0.00
## Diving 0.00
## Feasibility Pump 0.00
## Local Search 0.00 1 0
## Restricted Search 0.00
## Rins Search 0.00
## Local Branching 0.00
##
## =========================== Cuts ==========================
## Accepted: 0
## Added to LPs: 0
## Deleted from LPs: 0
## Removed because of bad coeffs: 0
## Removed because of duplicacy: 0
## Insufficiently violated: 0
## In root: 0
##
## Time in cut generation: 0.00
## Time in checking quality and adding: 0.00
##
## Time #Called In Root Total
## Gomory 0.00
## Knapsack 0.00
## Clique 0.00
## Probing 0.00
## Flowcover 0.00
## Twomir 0.00
## Oddhole 0.00
## Mir 0.00
## Rounding 0.00
## LandP-I 0.00
## LandP-II 0.00
## Redsplit 0.00
##
## ===========================================================
## Solution Found: Node 0, Level 0
## Solution Cost: -650.0000000000
## +++++++++++++++++++++++++++++++++++++++++++++++++++
## User indices and values of nonzeros in the solution
## +++++++++++++++++++++++++++++++++++++++++++++++++++
## 1 1.0000000000
## 2 1.0000000000
## 4 1.0000000000
## 11 1.0000000000
## Solution
A “raw" solution can be generated using:
get_solution(m,x[i, j]) %>%
filter(value > 0)## variable i j value
## 1 x 2 1 1
## 2 x 3 1 1
## 3 x 1 2 1
## 4 x 4 3 1We want to clean this up a bit:
cat("Status:",solver_status(m),"\n")
cat("Objective:",objective_value(m),"\n")
get_solution(m,x[i, j]) %>%
filter(value > 0) %>%
mutate(Treatment = treatments[j],Store = i) %>%
select(Store,Treatment)## Status: optimal
## Objective: 650
## Store Treatment
## 1 2 A
## 2 3 A
## 3 1 B
## 4 4 CReferences
- The original problem is from: Mixed linear integer optimization - Optimize Profit based on different treatments to all the stores in R, http://stackoverflow.com/questions/43446194/optimize-profit-based-on-different-treatments-to-all-the-stores-in-r
- Dirk Schumacher, OMPR: R package to model Mixed Integer Linear Programs, https://github.com/dirkschumacher/ompr
- Symphony MIP Solver: https://projects.coin-or.org/SYMPHONY
- Paul Rubin, MIP Models in R with OMPR, http://orinanobworld.blogspot.com/2016/11/mip-models-in-r-with-ompr.html
No comments:
Post a Comment