{lpSolve} and {lpSolveAPI}Reference: https://www.r-bloggers.com/linear-programming-in-r/
All credits to the writer of that post. This is just me fooling around with the code they wrote.
“This is an example of linear optimization that I found in the book “Modeling and Solving Linear Programming with R” by Jose M. Sallan, Oriol Lordan and Vincenc Fernandez. The example is named “Production of two models of chairs” and can be found at page 57, section 3.5. I’m going to solve only the first point."
A company produces two models of chairs: 4P and 3P. The model 4P needs 4 legs, 1 seat and 1 back. On the other hand, the model 3P needs 3 legs and 1 seat. The company has a initial stock of 200 legs, 500 seats and 100 backs. If the company needs more legs, seats and backs, it can buy standard wood blocks, whose cost is 80 euro per block. The company can produce 10 seats, 20 legs and 2 backs from a standard wood block. The cost of producing the model 4P is 30 euro/chair, meanwhile the cost of the model 3P is 40 euro/chair. Finally, the company informs that the minimum number of chairs to produce is 1000 units per month. Define a linear programming model, which minimizes the total cost (the production costs of the two chairs, plus the buying of new wood blocks).
# Solution with lpSolve -----------Here are the coefficients of the decision variables:
Therefore, the obj function is:
\(Cost = 30*4P + 40*3P + 80*WoodenBlocks\)
## Set the coefficients of the decision variables -> C
C <- c(30, 40, 80)There is one row for each constraint, and one column for each decision variable.
First row is for the Seat constraint. It says that:
Second row is for the Legs constraint
Third row is for the Backs constraint
Fourth row is for the Min production constraint
# Create constraint martix B
A <- matrix(c(1, 1, -10,
4, 3, -20,
1, 0, -2,
1, 1, 0), nrow=4, byrow=TRUE)
# Right hand side for the constraints
B <- c(500, 200, 100, 1000)
# Direction of the constraints
constranints_direction <- c("<=", "<=", "<=", ">=")lp function# Find the optimal solution
optimum <- lp(direction="min",
objective.in = C,
const.mat = A,
const.dir = constranints_direction,
const.rhs = B,
all.int = T)
str(optimum)## List of 28
## $ direction : int 0
## $ x.count : int 3
## $ objective : num [1:3] 30 40 80
## $ const.count : int 4
## $ constraints : num [1:5, 1:4] 1 1 -10 1 500 4 3 -20 1 200 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "" "" "" "const.dir.num" ...
## .. ..$ : NULL
## $ int.count : int 3
## $ int.vec : int [1:3] 1 2 3
## $ bin.count : int 0
## $ binary.vec : int 0
## $ num.bin.solns : int 1
## $ objval : num 48680
## $ solution : num [1:3] 420 580 161
## $ 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
## - attr(*, "class")= chr "lp"# Print status: 0 = success, 2 = no feasible solution
print(optimum$status)## [1] 0# Display the optimum values for x_4p, x_3p and x_w
best_sol <- optimum$solution
names(best_sol) <- c("x_4p", "x_3p", "x_w")
print(best_sol)## x_4p x_3p x_w
## 420 580 161# Check the value of objective function at optimal point
print(paste("Total cost: ", optimum$objval, sep=""))## [1] "Total cost: 48680"#*********************************************************
# Output #
#*********************************************************
# [1] 0
# x_4p x_3p x_w
# 420 580 161
# "Total cost: 48680"
rm(optimum, constranints_direction, best_sol)
#*********************************************************# Solution with lpSolveAPI -------------------
# Let's try to solve the problem again using lpSolveAPI
# Use lpSolveAPI
require(lpSolveAPI)
# Set 4 constraints and 3 decision variables
lprec <- make.lp(nrow = 4, ncol = 3)
# Set the type of problem we are trying to solve
lp.control(lprec, sense="min")## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] -1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "minimize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"# Set type of decision variables
set.type(lprec, 1:3, type=c("integer"))
# Set objective function coefficients vector C
set.objfn(lprec, C)
# Add constraints
add.constraint(lprec, A[1, ], "<=", B[1])
add.constraint(lprec, A[2, ], "<=", B[2])
add.constraint(lprec, A[3, ], "<=", B[3])
add.constraint(lprec, A[4, ], ">=", B[4])
# Display the LPsolve matrix
lprec## Model name:
## C1 C2 C3
## Minimize 30 40 80
## R1 0 0 0 free 0
## R2 0 0 0 free 0
## R3 0 0 0 free 0
## R4 0 0 0 free 0
## R5 1 1 -10 <= 500
## R6 4 3 -20 <= 200
## R7 1 0 -2 <= 100
## R8 1 1 0 >= 1000
## Kind Std Std Std
## Type Int Int Int
## Upper Inf Inf Inf
## Lower 0 0 0# Solve problem
solve(lprec)## [1] 0# Get the decision variables values
get.variables(lprec)## [1] 420 580 161# Get the value of the objective function
get.objective(lprec)## [1] 48680# Note that the default boundaries on the decision variable are c(0, 0, 0) and c(Inf, Inf, Inf)
get.bounds(lprec)## $lower
## [1] 0 0 0
##
## $upper
## [1] Inf Inf Inf# Boundaries can be set with following function
#lpSolveAPI::set.bounds()
#*********************************************************
# Output #
#*********************************************************
# [1] 420 580 161
# [1] 48680