Search code examples
roptimizationbinaryintegerlpsolve

Integer optimization with binary choice variable in R

I am working on an optimization problem in which my choice variable (to build or not) is a binary variable. The photo attached shows the optimal build choices (solved using Solver in Excel).

I am trying to replicate this problem in R. Minimizing the sum of expected value (EV) with the budget constraint ≤ 30 and build choices as the choice variable.

The problem is that EV changes with the given combination of build decisions. The general formula for it is:

EV = probability of fire (p)*Damage(when there is a fire) + (1-p)*Damage (no fire) + Cost (K).

Damage is 0 when there is no fire. High damage when there is a fire and the plot does not have a tower and low damage when the plot has a tower built on it.

So EV = p*Damage(fire state) + Cost

Excel sheet of problem in hand.

Plot i: Ten plots of land (size or other attributes do not matter). H: High damage. High damage only if there is no build decision. L: Low damage. Low damage if we chose to build a tower in a given plot. p: Probability of fire. K: Cost of building watch tower. Given B profile: B is a binary variable. B=1 signifies we chose to built in the given plot. EV(B=1): Expected value if we chose to build. EV(B=0): Expected value if we do not build. EV: Final expected value. If we build, EV = EV(B=1) else EV = EV(B=0). Budget constraint: If we build, it is equal to the cost of building, else it is 0.

Solution

  • You can express your problem as binary linear one by having two variables for each row. Those pair of varibles are connected in following way:

    • their costs are EV1 and EV0
    • their budgets are K and 0
    • their sum is exactly 1

    data

    data <- data.frame(
  H = c(50, 55, 58, 98, 60, 78, 88, 69, 78, 63),
  L = c(20, 11, 5, 12, 15, 22, 32, 15, 8, 17),
  p = c(0.1, 0.2, 0.05, 0.07, 0.5, 0.15, 0.09, 0.02, 0.01, 0.15),
  K = c(6, 8, 6, 5, 4, 8, 9, 7, 4, 2)
)

data$EV1 <- data$L * data$p + data$K
data$EV0 <- data$H * data$p

    solution

    library(lpSolve)

n <- nrow(data)

total_budget <- 30

res <- lp(
  direction = "min",
  objective = c(data$EV1, data$EV0),
  const.mat = rbind(c(data$K, rep(0, n)), cbind(diag(n), diag(n))),
  const.dir = c("<=", rep("=", n)),
  const.rhs = c(total_budget, rep(1, n)),
  all.bin = TRUE
)

res$solution[seq_len(n)]
[1] 0 1 0 1 1 1 0 0 0 1
res$objval
[1] 61.37
sum(data$K[res$solution[seq_len(n)] == 1])
[1] 27

    where total_budget is upper limit on your budget.