Search code examples
rmathematical-optimizationlinear-programminglpsolve

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.


Solution

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