## 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.

## 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

- 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