Linear Programming Using LpSolve and LpSolveAPI
Introduction
The LpSolveAPI and lpSolve are package provides in R to solve Mixed Integer Linear Programming (MILP) with support for for pure linear, (mixed) integer/binary, semi-continuous models.
Structure of linear program model The linear programming problem consists in optimizing (that is, either minimize or maximize) the value of a linear objective function of a vector of decision variables, considering that the variables can only take the values defined by a set of linear constraints. Linear programming is a case of mathematical programming, where objective function and constraints are linear.
A formulation of a linear program in its canonical form is:
\[\begin{align*} \text{Maximize} \quad &z = c_1x_1 + c_2x_2 + \ldots + c_nx_n \\ \text{Subject to} \quad &a_{11}x_1 + a_{12}x_2 + \ldots + a_{1n}x_n \leq b_1 \\ &a_{21}x_1 + a_{22}x_2 + \ldots + a_{2n}x_n \leq b_2 \\ &\vdots \\ &a_{m1}x_1 + a_{m2}x_2 + \ldots + a_{mn}x_n \leq b_m \\ &x_1, x_2, \ldots, x_n \geq 0 \end{align*}\]The model has the following elements:
-
An objective function of n decision variables \(x_j\). Decision variables are affected by cost coeffecient \(c_j\)
-
It also encompasses a set of
m
constraints. In these constraints, a linear combination of the variables affected by coefficients \(a_{ij}\) must be less than or equal to its corresponding right-hand side value \(b_i\). Constraints may have signs greater than or equal, and equalities are also permissible. -
The bounds of the decision variables, in this case all decision variables has to be non-negative.
Problem Statement
“This example of linear optimization can be 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.
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).
lpSolve Solution
library(lpSolve)
C <- c(30, 40, 80)
# 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_direc <- c("<=", "<=", "<=", ">=")
# Find the optimal solution
optimal <- lp(direction="min",
objective.in = C,
const.mat = A,
const.dir = constranints_direc,
const.rhs = B,
all.int = T)
str(optimal)
print(optimal$status)
# Display the optimum value
best_solution <- optimal$solution
names(best_solution)
print(best_solution)
# Check the value of objective function at optimal point
print(paste("Total cost: ", optimal$objval, sep=""))
# print(best_solution)
# [1] 420 580 161
# print(paste("Total cost: ", optimum$objval, sep=""))
# "Total cost: 48680"
lpSolveAPI Solution
library(lpSolveAPI)
## Cost=30*4P + 40*3P + 80*WoodenBlocks
## Set the coefficients of the decision variables -> C
C <- c(30, 40, 80)
# 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)
# set 4 constraints
lprec <- make.lp(nrow=4, ncol=3)
lp.control(lprec, sense="min")
# 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])
lprec
# Solve the problem
solve(lprec)
# get the decision variable values
get.variables(lprec)
# get objective function
get.objective(lprec)
# get.variables(lprec)
# [1] 420 580 161
# get.objective(lprec)
# [1] 48680