CODEWITHSUNDEEP
rmathematical-optimizationlinear-programminglpsolve
coolhand
coolhand

Reputation: 2073

Applying lpSolveAPI

I have a number of variables denoting equipment I can replace. When replaced, they improve an impact metric, [I]. Each also has an associated annual cost savings [S] and replacement cost [C].

n <- 1000 # variable count

# impact
# negative for use in minimization function
I <- -rnorm(n, mean=20000, sd=8000)
# cost savings
s <- rnorm(n, mean=2500, sd=1000)
# replacement cost
c <- rnorm(n, mean=15000, sd=5000)

I want to select which components to replace to maximize the impact across the board within a budget, while still making sure the entire project (as a whole) has meets a simple payback goal.

payback_goal <- 3
budget <- 1000000

This problem is described by the equations below.

enter image description here

I'm struggling to set this up with lpSolveAPI. Specifically, I don't know how to incorporate Eq. 3.

library(lpSolveAPI)
m <- 2 # number of constraints, disregarding binary constraint set by type
my.lp <- make.lp(m, n)
set.row(my.lp, 1, c)
# i don't think this is the correct way to set up the calculation
set.row(my.lp, 2, c/s) # cost divided by savings

# create obj fn coefficients
set.objfn(my.lp, I)
set.type(my.lp, 1:n, "binary")

eq <- c("<=", "<=")
set.constr.type(my.lp, eq)
set.rhs(my.lp, c(budget, payback_goal))

solve(my.lp)
soln1 <- get.variables(my.lp)

# total cost
s1_cost <- sum(soln1 * c)
# total impact
s1_impact <- -get.objective(my.lp)
# total simple payback
s1_pb <- sum(soln1*c) / sum(soln1*s)

I've tried a somewhat hacky workaround, that assumes the total cost will be reasonably close to the budget, which makes the third equation

enter image description here

But this is just an approximation and I'd like to know how to implement it more exactly within the R code.

Upvotes: 4

Views: 183

Answers (2)

ThomasIsCoding
ThomasIsCoding

Reputation: 102770

As pointed out by @Erwin Kalvelagen, you can reformulate the constraint to make it linear.

But, you should be aware of that, the transformation is conditional on sum([X]*[S]), that means

  • If sum([X]*[S])>0, we should have sum([X]*([C]-payback*[S])) <= 0
  • If sum([X]*[S])<0, we should have sum([X]*([C]-payback*[S])) >= 0

Then, combine two cases above (compare and take the minimum), we can margin out the impact by sum([X]*[S]).

code example (in CVXR)

set.seed(0)
n <- 1000 # variable count

# impact
# negative for use in minimization function
I <- -rnorm(n, mean = 20000, sd = 8000)
# cost savings
s <- rnorm(n, mean = 2500, sd = 1000)
# replacement cost
c <- rnorm(n, mean = 15000, sd = 5000)


library(CVXR)
X <- Variable(1, n, boolean = TRUE)
payback_goal <- 3
budget <- 1000000

obj <- Minimize(X %*% I)
# case 1: conditional on X %*% s > 0
constr1 <- list(
  X %*% c <= budget,
  X %*% (c - payback_goal * s) <= 0,
  X %*% s >= 0
)
problem1 <- Problem(obj, constr1)
sol1 <- solve(problem1)

# case 2: conditional on X %*% s < 0
constr2 <- list(
  X %*% c <= budget,
  X %*% (c - payback_goal * s) >= 0,
  X %*% s <= 0
)
problem2 <- Problem(obj, constr2)
sol2 <- solve(problem2)

# compare two cases and take the minimum between them as the solution
if (round(sol1$getValue(X)) %*% I > round(sol2$getValue(X)) %*% I) {
  sol <- sol2
} else {
  sol <- sol1
}

and you will see that

> round(sol1$getValue(X)) %*% I
         [,1]
[1,] -3064996

> round(sol2$getValue(X)) %*% I
          [,1]
[1,] -287548.7

such that sol is assigned with sol1

Upvotes: 1

Erwin Kalvelagen
Erwin Kalvelagen

Reputation: 16797

You can linearize (3) as follows, using your notation:

 sum([X]*[C]) - payback*sum([X]*[S]) <= 0

That can be rewritten as;

 sum([X]*([C]-payback*[S])) <= 0

This is valid as long as sum([X]*[S])>0.

Upvotes: 5

Related Questions