Optimization from Excel to R

Question

I am struggling with solving below optimization in R. It was easy for me to do with excel solver but I am not able to do same problem in R. I don’t have much experience in optimization in R Problem is to assign when a particular activity should be a done over the time period. A1, A2,…, A12 are 12 activities which are done in a field. 0 stands for not assigned, 1 stands for assigned.

Constraints are -

• A2 should happen after A1 is completed, A3 should occur after A2 is complete and so on
• A3, A5, A6 - only two of these activities can happen on a particular day
• Each activity should happen in a field ( in excel solver, this constraint is defined by taking row sum equals to one

It would be a great help, if anyone can help in solving this problem. I have attached image for better understanding of this problem. Also, please let me know if there is any site or book, rich in examples on "how to solve optimization problems in R". Thanking you in anticipation.

Answer 1

You can solve your problem using the lpSolve package.

You will require a cost vector and constraints information. Constraint information in fed into the function in the following structure:

• lhs: a "left hand side" matrix of coefficients, one for each decision variable
• dir: a "direction", namely <, <=, ==, >=, >
• rhs: a "right hand side" as numeric value

To build your list of constraints, I find it useful to think of each decision you could possibly take as an X (becomes a column in the lhs table) and each constraint you would like to define as a separate equation (becomes a row in the lhs table, with a corresponding value in dir and rhs, respectively)

Let us start with all possible decisions:

library(tidyverse)
library(stringr)

# What are the decision variables? ----

# Which action to take
actions <- str_c('A',seq(1:12) %>% formatC(width = 2, flag = '0'))
actions
#[1] "A01" "A02" "A03" "A04" "A05" "A06" "A07" "A08" "A09" "A10" "A11" "A12"

# When to take it
timings <- str_c('T',seq(1:12) %>% formatC(width = 2, flag = '0'))
timings
#[1] "T01" "T02" "T03" "T04" "T05" "T06" "T07" "T08" "T09" "T10" "T11" "T12"

# List of all possible decisions is this:
decisions <- expand.grid(actions, timings)

# Convert it to a vector
decision_variables <- str_c(decisions[,1], '_', decisions[,2])

# You also need a cost vector. 
# We'll use a value increasing as a function of timings, 
# as this will penalize "late" actions?
cost <- rep(seq(1:length(timings)), length(actions)) %>% sort

Each element of decision_variables is one possible action (i.e. take an action at a given time. Now we can start narrowing down the options available to the solver by introducing constraints.

---

First type of constraint: Each option should only be taken once! (This is actually your third, but I start with this as it's the simplest one) We can formulate that like this:

# Create a matrix with one column per possible decision
# and one row per action (for now)
lhs <- matrix(0, 
  nrow = length(actions), 
  ncol = length(decision_variables), 
  dimnames = list(
    actions, 
    decision_variables))

# Each action should only be taken once!
for (i in 1:length(actions)) {
  # Which fields does an action occur in?
  this_action <- str_detect(colnames(lhs), actions[i])
  # Set their coefficients to 1
  lhs[i,this_action] <- 1
}
# create corresponding dir and rhs values
dir <- rep('==', length(actions))
rhs <- rep(1, length(actions))

You can see that we set the coefficient of all X's (decisions) that contain the action in question to 1. In our final solution, each X will take the value of either 0 or 1. If X is zero, the coefficient will be irrelevant. If X is 1, the coefficient will be added to the sum of the lhs and compared using dir to the rhs value.

Here, our constraint is that the sum of coefficient * X == 1 for each of the constraints we just introduced. Coefficients are 1 for all possible decisions containing a given action. Ergo, a solution is only valid if any given action is only taken once.

---

Second constraint: Only two of c('A03', 'A05', 'A06') should co-occur on a given day.

Again, we generate one line per constraint. In this case, I think we need one constraint per day. We append the values we generate to the already existing lhs, dir, and rhs variables:

# only one of A3, A5, A6 at any given time.
# One constraint for each timestep
for (j in timings) {
  lhs <- rbind(lhs, ifelse(str_detect(decision_variables, paste0('A0[356]{1}_',j)), 1, 0))
  dir <- c(dir, '<=')
  rhs <- c(rhs, 2)
}

---

Placeholder for third constraint

---

Presto, we've formulated our problem. Now let lpSolve crunch the numbers! You can feed our problem into the algorithm like this:

library(lpSolve)

# Run lpSolve to find best solution
solution <- lp(
  # maximise or minimise the objective function?
  direction = 'min',
  # coefficients of each variable
  objective.in = cost,
  const.mat = lhs, 
  const.dir = dir, 
  const.rhs = rhs)

# Extract the values of X for the best solution:
print(solution$solution)

# Convert it into ta matrix of the format you are familiar with
matrix(solution$solution, 
       nrow = length(timings), 
       ncol = length(actions),
       dimnames = list(actions, timings))

Does this do what you need?

Any questions?