Monday, April 17, 2017

MIP model in R using the OMPR package

A small model with data in R is a good opportunity to try out the OMPR modeling tool (2). Also a first time for me to see the solver Symphony (3) in action.  

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  125

How 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:

image

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    125

OMPR 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 4
treatments <- 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 <- 300

MIP 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     1

We 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         C
References
  1. 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
  2. Dirk Schumacher, OMPR: R package to model Mixed Integer Linear Programs, https://github.com/dirkschumacher/ompr
  3. Symphony MIP Solver: https://projects.coin-or.org/SYMPHONY
  4. Paul Rubin, MIP Models in R with OMPR, http://orinanobworld.blogspot.com/2016/11/mip-models-in-r-with-ompr.html