Linear Programming Modeling

Case Study

A manufacturing manager is in charge of minimizing the total costs (raw materials, labor and storage costs) of the following four months. In Table 3.1 can be found the cost of raw materials of one unit of final product, the demand of final product and the working hours available for each month. Labor costs are of 12 € per hour, and only worked hours are payed. Each unit of final product needs 30 minutes of labor. Storage costs are equal to 2 € for each unit stored at the end of the month. Any unit produced at a given month can be used to cover the demand of the same month, or be stored to cover the demand of months to come. At the beginning of month 1 there is no stock, and there are no minimum stock requirements for any month.

#Table 3.1
Month <- c("Unit cost", "Demand (units)", "Working hours available")
One   <- c( 6,100,200)
Two   <- c( 8,200,200)
Three <- c(10,150,150)
Four  <- c(12,400,150)
Table <- data.frame(Month, One, Two, Three, Four)
Table

##                     Month One Two Three Four
## 1               Unit cost   6   8    10   12
## 2          Demand (units) 100 200   150  400
## 3 Working hours available 200 200   150  150

Problem definition

First, we define the following variables :

• \(q_i\) is the quantity produced in month \(i\).
• \(s_i\) is the stock at the end of month \(i\).

We would like to minimize the total cost so we must set our objective function as follows

\[cost(s_i, q_i) = 2 s_i + 12 q_i = MIN(cost)\]

Now, we will solving the problem using lpSolve and lpSolveAPI

Solution with lpSolve

lpsolve: is callable from R via an extension or module. As such, it looks like lpsolve is fully integrated with R. Matrices can directly be transferred between R and lpsolve in both directions.

# install.packages("lpSolve")
library(lpSolve)

## Warning: package 'lpSolve' was built under R version 3.6.2

# set the coefficients of the objective function
o.func <- c(2, 12)

# set matrix corresponding to coefficients of constraint by rows
# Don't consider the non-negative constraint, it is automatically assumed
cons   <- matrix(c(100, 400,
                   200, 400,
                   150, 300,
                   400, 300), 
                 nrow = 4, byrow = T)

# set right hand side for the constraints
rhs <- c(6, 8, 10, 12)

# set direction of the constraints
c.direct <- c(">=", ">=", ">=", ">=")

# The Final value (z) : 0 = success, 2 = no feasible solution
lp("min", o.func, cons, c.direct, rhs)

## Success: the objective function is 0.1333333

# Variables final values
lp("min", o.func, cons, c.direct, rhs)$solution

## [1] 0.06666667 0.00000000

# sensitivities
lp("min", o.func, cons, c.direct, rhs, compute.sens=T)$sens.coef.from

## [1] 0 4

lp("min", o.func, cons, c.direct, rhs, compute.sens=T)$sens.coef.to

## [1] 6e+00 1e+30

# Dual Values (first dual of the constraints and then dual of the variables)
# Duals of the constraints and variables are mixed
lp("min", o.func, cons, c.direct, rhs, compute.sens=T)$duals

## [1] 0.00000000 0.00000000 0.01333333 0.00000000 0.00000000 8.00000000

So, the minimum value of z is 0.13, where \(q_i\) is 0.067, \(s_i\) is 0.00 and the sensitivity coefficient changes from 0 and 4 to 6e+00 and 1e+30.

Solution with lpSolveAPI

The lpSolveAPI R package is a second implementation of an interface of lpsolve to R. It provides an R API mirroring the lp_solve C API and hence provides a great deal more functionality but has a steeper learning curve.

# install.packages("lpSolveAPI")
library(lpSolveAPI)

## Warning: package 'lpSolveAPI' was built under R version 3.6.3

# Set 4 constraints and 2 decision variables
lprec   <- make.lp(nrow = 4, ncol = 2)

# 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 objective function coefficients vector
set.objfn(lprec, o.func)

# Add constraints
add.constraint(lprec, cons[1, ], ">=", rhs[1])
add.constraint(lprec, cons[2, ], ">=", rhs[2])
add.constraint(lprec, cons[3, ], ">=", rhs[3])
add.constraint(lprec, cons[4, ], ">=", rhs[4])

# Solve problem
solve(lprec)

## [1] 0

# Display the lpsolve matrix
lprec

## Model name: 
##             C1    C2          
## Minimize     2    12          
## R1           0     0  free   0
## R2           0     0  free   0
## R3           0     0  free   0
## R4           0     0  free   0
## R5         100   400    >=   6
## R6         200   400    >=   8
## R7         150   300    >=  10
## R8         400   300    >=  12
## Kind       Std   Std          
## Type      Real  Real          
## Upper      Inf   Inf          
## Lower        0     0

# Get the decision variables values
get.variables(lprec)

## [1] 0.06666667 0.00000000

# Get the value of the objective function
get.objective(lprec)

## [1] 0.1333333

# Note that the default boundaries on the decision variable are c(0, 0) and c(Inf, Inf)
get.bounds(lprec)

## $lower
## [1] 0 0
## 
## $upper
## [1] Inf Inf

# Boundaries can be set with following function
lpSolveAPI::set.bounds(lprec)