## The Problem

Suppose a farmer has 75 acres on which to plant two crops: wheat and barley. To produce these crops, it costs the farmer (for seed, fertilizer, etc.) $120 per acre for the wheat and $210 per acre for the barley. The farmer has $15000 available for expenses. But after the harvest, the farmer must store the crops while awaiting avourable market conditions. The farmer has storage space for 4000 bushels. Each acre yields an average of 110 bushels of wheat or 30 bushels of barley. If the net profit per bushel of wheat (after all expenses have been subtracted) is $1.30 and for barley is $2.00, how should the farmer plant the 75 acres to maximize profit?

Problem definition

First, we need to translate the problem in a mathematical way. Let’s define the following variables

Let’s define the optimization function and the constraints.

Maximize : \[g = (110)(1.3)x + (30)(2.00)y = 143x + 60y\]

Subject to \[120x + 210y \leq 15000\] \[110x + 30y \leq 4000\] \[x + y \leq 75\] \[x \geq 0\] \[y \geq 0\] We will be solving this problem in both lpSolve and lpSolveAPI

R packages

Last week you already learn a little bit how to implement this is in R. There are many different packages which can solve this kind of problems but my favourites are lpSolve and lpSolveAPI which is a kind of API built on top of lpSolve. Here, I will try to use both of them for you to know more.

A nice feature about the lpSolve package is that you can specify the direction of the constraints. Indeed in our case the last constraint of minimum number of chairs produced does not fit in with the mathematical definiton of the problem that we gave in the previous paragraph. Here we can either change the signs (and therefore the inequality direction) or specify the inequality direction in lpSolve. I’ll do this second way.

We can set that all the variables should be integers by setting the argument “all.int=TRUE” in the lpSolve::lp function. Of course, lpSolve can work with both integers and real numbers. It’s also particularly important to check the status at the end of the execution: if it is 0, then a solution has been found but if it is 2 then this means that there is no feasible solution.

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. The complete interface is written in C so it has maximum performance.

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

Objective Function

Here are the coefficients of the decision variables:

  • Cost of each unit of \(x\) is $ $143$
  • Cost of each unit of \(y\) is $ $60$

Therefore, the obj function is:

\[g = (110)(1.3)x + (30)(2.00)y = 143x + 60y\]

# Set the coefficients of the decision variables -> C
obj.fun <- c(143,60)        
obj.fun
## [1] 143  60

Constraint Matrix

There is one row for each constraint, and one column for each decision variable.

  • First row is for the Seat constraint. It says that:
    • Each unit of 4P uses 1 seat from the seat inventory
    • Each unit of 3P uses 1 seat from the seat inventory
    • Each unit of wooden blocks adds 1 seat to the seat inventory
  • 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
constr <- matrix(c(120,210,
              110,30,
              1, 1), nrow=3, byrow=TRUE)
constr
##      [,1] [,2]
## [1,]  120  210
## [2,]  110   30
## [3,]    1    1
# Right hand side for the constraints
rhs <- c(15000,4000,75)
rhs
## [1] 15000  4000    75
# Direction of the constraints
const.dir  <- c("<=", "<=", ">=")
const.dir
## [1] "<=" "<=" ">="

The Optimum Result

# Find the optimal solution
optimum <-  lp(direction="max",
               obj.fun ,
               constr,
               const.dir,
               rhs,
               all.int = T)

str(optimum)
## List of 28
##  $ direction       : int 1
##  $ x.count         : int 2
##  $ objective       : num [1:2] 143 60
##  $ const.count     : int 3
##  $ constraints     : num [1:4, 1:3] 120 210 1 15000 110 30 1 4000 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "" "" "const.dir.num" "const.rhs"
##   .. ..$ : NULL
##  $ int.count       : int 2
##  $ int.vec         : int [1:2] 1 2
##  $ bin.count       : int 0
##  $ binary.vec      : int 0
##  $ num.bin.solns   : int 1
##  $ objval          : num 6460
##  $ solution        : num [1:2] 20 60
##  $ 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", "y") 
print(best_sol)
##  x  y 
## 20 60
# Check the value of objective function at optimal point
print(paste("Total cost: ", optimum$objval, sep=""))
## [1] "Total cost: 6460"

Thus, to achieve the maximum profit ($6460), the farmer should plant 20 acres pf Wheat and 60 acres of Barley.

# Disconnect from the model and the optimum solution
rm(optimum, constranints_direction, best_sol)
## Warning in rm(optimum, constranints_direction, best_sol): object
## 'constranints_direction' not found

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)

Let’s try to solve the problem again using lpSolveAPI:

# Set 4 constraints and 3 decision variables
lprec <- make.lp(0, 2)
# Set the type of problem we are trying to solve
lp.control(lprec, sense="max")
## $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] "maximize"
## 
## $simplextype
## [1] "dual"   "primal"
## 
## $timeout
## [1] 0
## 
## $verbose
## [1] "neutral"
# Set type of decision variables
set.type(lprec, 1:2, type=c("integer"))
# Set objective function coefficients vector C
set.objfn(lprec , c(143,60))
# Add constraints
add.constraint(lprec, c(120, 210), "<=", 15000)
add.constraint(lprec, c(110, 30), "<=", 4000)
add.constraint(lprec, c(1, 1), "<=", 75)
# Let's have a look at the LP problem that we have defined
lprec
## Model name: 
##            C1   C2           
## Maximize  143   60           
## R1        120  210  <=  15000
## R2        110   30  <=   4000
## R3          1    1  <=     75
## Kind      Std  Std           
## Type      Int  Int           
## Upper     Inf  Inf           
## Lower       0    0
# Solve problem
solve(lprec)
## [1] 0
# Getting the optimized profit value
get.objective(lprec)
## [1] 6266
# Finally get the value for how many acres of Wheat and Barley to be planted
get.variables(lprec)
## [1] 22 52
# 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
## 
## $upper
## [1] Inf Inf
# Boundaries can be set with following function
lpSolveAPI::set.bounds(lprec)

Thus, to achieve the maximum profit ($6266), the farmer should plant 22 acres pf Wheat and 52 acres of Barley.