linear programming in R | limits to what it can do, or my mistake?

Question

Question for 'experts' in LP using R (using the lpSolve package, say) -- which does not apply to me for the sort of problem I describe below.

I've run any number of LP's using lpSolve in R, but all of them to date have objective and constraint functions that both contain the same variables. This lets you set up a LHS and RHS matrix/vector that are symmetrical.

But, for a problem a student posed in class, I'm stuck with how to do it in R, if its even possible (its trivial in Maxima, Maple...even using Solver in Excel, but I haven't been remotely successful in getting anything to work in R).

Suppose you have a production system that at 4 sequential time steps generate 640, 825, 580, and 925 units. At each time step, you need to decide how many of those units need to be 'quality control' (QC) checked in some fashion, subject to some constraints.

 --> at no point in time can the number of units in the system be >1000
 --> at the end of the production cycle, there can be no units left
 --> 'QC checking' costs money, varying as a function of the time step -- 35, 55, 50 and 65 for each unit, for each time step in turn. 

Objective is to minimize total cost. The total cost objective function is trivial. Let p1 = number sent out time step 1, p2 number sent out at time step 3, and so on. So, total cost function we want to minimize is simply

  cost=(35*p1)+(55*p2)+(50*p3)+(65*p4)

where p1+p2+p3+p4=(640+825+580+925)=2970 (i.e., all the products get checked). The question is, what number do you send out at each time step to minimize cost?

Where I get hung up in R is the fact that if I let t(i) be the number of products at each time step, then

 t1=640,
 t2=t1-p1+825
 t3=t2-p2+580
 t4=t3-p3+925

such that t1+t2+t3+t4=2970 (as it must), with additional constraints being

 p1<=t1, p2<=t2, p3<=t3, p4<=t4, {t1..t4}<=1000, and t4-p4=0. 

There may be algebraic ways to reduce the number of functions needed to describe the constraints, but I can't for the life of me see how I can create a coefficient matrix (typically, the LHS) since each line of said matrix, which corresponds to the constraints, needs to be a function of the unknowns in the objective function -- being, p1, p2, p3 and p4.

In Maple (for example), this is trivial:

 cost:=35*p10+55*p12+50*p14+65*p16;      
 cnsts:={t10=640,t12=t10-p10+825,t14=t12-p12+580,t16=t14-p14+925,t16-p16=0,p10<=t10,p12<=t12,p14<=t14,p16<=t16,t10<=1000,t12<=1000,t14<=1000,t16<=1000};
 Minimize(cost,cnsts,assume={nonnegative});

which yields (correctly):

 p1=640, p2=405, p3=1000, p4=925

for minimized cost of 154800.

Took only a minute to also set this up in Maxima, and using Solver in Excel. But danged if I can suss out any way to do this in R.

Pointers to the obvious welcome.

Answer 1

1) CVXR The CVXR package allows a specification similar to the one shown in the question. (SO keeps falsely claiming I used AI to generate this answer but I did not -- I wrote it myself.)

library(CVXR)

p10 <- Variable(1); p12 <- Variable(1); p14 <- Variable(1); p16 <- Variable(1)
t10 <- Variable(1); t12 <- Variable(1); t14 <- Variable(1); t16 <- Variable(1)

cost <- 35*p10+55*p12+50*p14+65*p16

cnsts <- list(t10==640,t12==t10-p10+825,t14==t12-p12+580,t16==t14-p14+925,t16-p16==0,p10<=t10,p12<=t12,p14<=t14,p16<=t16,t10<=1000,t12<=1000,t14<=1000,t16<=1000)

objective <- Minimize(cost)
prob <- Problem(objective, cnsts)
soln <- solve(prob)

soln$status
$$ [1] "optimal"

soln$value
$$ [1] 154800

soln$getValue(p10); soln$getValue(p12); soln$getValue(p14); soln$getValue(p16)
## [1] 640
## [1] 405
## [1] 1000
## [1] 925

2) lpSolve With lpSolve we have the tedious job of forming the inputs, particularly the constraint matrix, but we can do it automatically by taking the derivative wrt each variable since the problem is linear.

library(lpSolve)

# inputs from question
cnsts <- "t10==640,t12==t10-p10+825,t14==t12-p12+580,t16==t14-p14+925,t16-p16==0,p10<=t10,p12<=t12,p14<=t14,p16<=t16,t10<=1000,t12<=1000,t14<=1000,t16<=1000"
cost <- "35*p10+55*p12+50*p14+65*p16"

# cnsts.vec is cnsts split into individual constraints
cnsts.vec <- cnsts |> strsplit(",") |> unlist()

# lhs minus rhs as a language object
cnsts.lang <- cnsts.vec |> 
  sub("[<=]+(.*)", "-\\(\\1\\)", x = _) |>
  sub("\\b(\\d+)\\b", "\\1*C", x = _) |>
  sapply(str2lang)

# right hand side constants
const.rhs <- - sapply(cnsts.lang, function(x) eval(D(x, "C")))

# v is all variables v <- c("p10", "p12", ..., "t16")    
v <- cnsts |> chartr(",", ";", x = _) |> parse(text = _) |> 
  all.vars(unique = TRUE) |> sort()

# constraint matrix
cnsts.mat <- sapply(v, function(v) sapply(cnsts.lang, function(x) eval(D(x, v))))
rownames(cnsts.mat) <- cnsts.vec

# objective
obj <- sapply(v, function(v) D(str2lang(cost), v))

# constraint directions: vector of "<=" or "=="
cnsts.dir <- gsub("[^=<]", "", cnsts.vec)

res <- lp("min", obj, cnsts.mat, cnsts.dir, const.rhs)
str(res)

giving

List of 29
 $ direction       : int 0
 $ x.count         : int 8
 $ objective       : Named num [1:8] 35 55 50 65 0 0 0 0
  ..- attr(*, "names")= chr [1:8] "p10" "p12" "p14" "p16" ...
 $ const.count     : int 13
 $ constraints     : num [1:10, 1:13] 0 0 0 0 1 0 0 0 3 640 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "p10" "p12" "p14" "p16" ...
  .. ..$ : chr [1:13] "t10==640" "t12==t10-p10+825" "t14==t12-p12+580" "t16==t14-p14+925" ...
 $ int.count       : int 0
 $ int.vec         : int 0
 $ bin.count       : int 0
 $ binary.vec      : int 0
 $ num.bin.solns   : int 1
 $ objval          : num 154800
 $ solution        : num [1:8] 640 405 1000 925 640 825 1000 925
 $ presolve        : int 0
 $ compute.sens    : int 0
 $ sens.coef.from  : num 0
 $ sens.coef.to    : num 0
 $ duals           : num 0
 $ duals.from      : num 0
 $ duals.to        : num 0
 $ scale           : int 196
 $ use.dense       : int 0
 $ dense.col       : int 0
 $ dense.val       : num 0
 $ dense.const.nrow: int 0
 $ dense.ctr       : num 0
 $ use.rw          : int 0
 $ tmp             : chr "Nobody will ever look at this"
 $ status          : int 0
 $ timeout         : int 0
 - attr(*, "class")= chr "lp"

Answer 2

Assuming you need to solve a variable vector x of size 8, where x is defined as (p1, p2, p3, p4, t1, t2, t3, t4) since both p_1..4 and t_1..4 can be treated as variables, you could construct the following constraints:

• the equality constraint Aeq %*% x == beq
• and the inequality constraint Aineq %*% x <= bineq

to minimize the inner product C %*% x with respect the the given C coefficients c(35, 55, 50, 65, 0, 0, 0, 0)

---

Below is an implementation with CVXR

library(CVXR)
# equality constraint
Aeq <- rbind(
  c(0, 0, 0, 0, 1, 0, 0, 0),
  c(1, 0, 0, 0, -1, 1, 0, 0),
  c(0, 1, 0, 0, 0, -1, 1, 0),
  c(0, 0, 1, 0, 0, 0, -1, 1),
  c(1, 1, 1, 1, 0, 0, 0, 0),
  c(0, 0, 0, -1, 0, 0, 0, 1)
)
beq <- c(640, 825, 580, 925, 2970, 0)

# inequality constraint
Aineq <- rbind(
  c(1, 0, 0, 0, -1, 0, 0, 0),
  c(0, 1, 0, 0, 0, -1, 0, 0),
  c(0, 0, 1, 0, 0, 0, -1, 0),
  c(0, 0, 0, 1, 0, 0, 0, -1),
  c(0, 0, 0, 0, 1, 0, 0, 0),
  c(0, 0, 0, 0, 0, 1, 0, 0),
  c(0, 0, 0, 0, 0, 0, 1, 0),
  c(0, 0, 0, 0, 0, 0, 0, 1)
)
bineq <- c(0, 0, 0, 0, 1000, 1000, 1000, 1000)

# coeffs and variable
C <- matrix(c(35, 55, 50, 65, 0, 0, 0, 0), 1)
x <- Variable(8, integer = TRUE)

# objective function and constraint list
obj <- Minimize(C %*% x)
constr <- list(
  Aineq %*% x - bineq <= 0,
  Aeq %*% x - beq == 0
)

# problem formulation and solution
prob <- Problem(obj, constr)
res <- solve(prob)

and you will see

# the first 4 entries of `x`, say, p1, p2, p3, and p4
> round(res$getValue(x))[1:4]
[1]  640  405 1000  925

# optimum value
> res$value
[1] 154800