Lab3: Linear Programming Modeling

Case Study 1

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

• \(w\) is the number of acres allotted to wheat.
• \(b\) is the number of acres allotted to barley.

Now we can define \(\hat X = \begin{pmatrix} w \\ b \end{pmatrix}\) as the decision variable vector. Note that it must be \(\hat X \geq 0\).

We would like to maximize the profit, so we must set our objective function as follows

\[ P(w, b) = (110*1.30) w + (30*2.00) b = MAX(profit) \]

Then the expression to be maximized, that is the profit, is clearly

\[P(w, b) = 143 w + 60 b = MAX(profit) \]

which means that \[\hat P = \begin{pmatrix} 143 \\ 60 \end{pmatrix}\].

The constraints can be set in the following ways:

1. Cost for seed, fertilizer, etc. \[ 120w + 210b \leq 15000 \]
2. Storage space \[ 110w + 30b \leq 4000 \]
3. Area of farmland \[ w + b \leq 75 \]
4. Acres allotted to wheat \[ w \geq 0 \]
5. Acres allotted to barley \[ b \geq 0 \]

We can now define the coefficient matrix \[A = \begin{pmatrix} 120 & 210 \\ 110 & 30 \\ 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{pmatrix}\], and

\[B = \begin{pmatrix} 15000 \\ 4000 \\ 75 \\ 0 \\ 0 \end{pmatrix}\].

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 P so it has maximum performance.

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

Objective Function

Here are the coefficients of the decision variables:

• The net profit per bushel of wheat \(w\) is \(\$1.30\)
• The net profit per bushel of barley \(b\) is \(\$2.00\)

Therefore, the obj function is:

\[ P(w, b) = (110*1.30) w + (30*2.00) b \] \[P(w, b) = 143 w + 60 b \]

# Set the coefficients of the decision variables -> P
P <- c(143, 60)        
P

## [1] 143  60

Constraint Matrix

• First row is for the cost for seed, fertilizer, etc.
• Second row is for the storage space
• Third row is for the area of farmland
• Fourth row is for acres allotted to wheat
• Fifth row is for acres allotted to barley

# Create constraint matrix 
A <- matrix(c(120, 210,
              110, 30,
              1, 1,
              1, 0,
              0, 1), nrow=5, byrow=TRUE)
A

##      [,1] [,2]
## [1,]  120  210
## [2,]  110   30
## [3,]    1    1
## [4,]    1    0
## [5,]    0    1

# Right hand side for the constraints
B <- c(15000, 4000, 75, 0, 0)
B

## [1] 15000  4000    75     0     0

# Direction of the constraints
constranints_direction  <- c("<=", "<=", "<=", ">=", ">=")
constranints_direction

## [1] "<=" "<=" "<=" ">=" ">="

The Optimum Result

# Find the optimal solution
optimum <-  lp(direction="max",
               objective.in = P,
               const.mat = A,
               const.dir = constranints_direction,
               const.rhs = B,
               all.int = T)

str(optimum)

## List of 28
##  $ direction       : int 1
##  $ x.count         : int 2
##  $ objective       : num [1:2] 143 60
##  $ const.count     : int 5
##  $ constraints     : num [1:4, 1:5] 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 6266
##  $ solution        : num [1:2] 22 52
##  $ 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 w, b
best_sol <- optimum$solution
names(best_sol) <- c("w", "b") 
print(best_sol)

##  w  b 
## 22 52

# Check the value of objective function at optimal point
print(paste("Maksimum profit: ", optimum$objval, sep=""))

## [1] "Maksimum profit: 6266"

# Disconnect from the model and the optimum solution
rm(optimum, constranints_direction, best_sol)

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 5 constraints and 2 decision variables
lprec <- make.lp(nrow = 5, ncol = 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 P
set.objfn(lprec, P)

# 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])
add.constraint(lprec, A[5, ], ">=", B[5])

# Display the LPsolve matrix
lprec

## Model name: 
##            C1   C2             
## Maximize  143   60             
## R1          0    0  free      0
## R2          0    0  free      0
## R3          0    0  free      0
## R4          0    0  free      0
## R5          0    0  free      0
## R6        120  210    <=  15000
## R7        110   30    <=   4000
## R8          1    1    <=     75
## R9          1    0    >=      0
## R10         0    1    >=      0
## Kind      Std  Std             
## Type      Int  Int             
## Upper     Inf  Inf             
## Lower       0    0

# Solve problem
solve(lprec)

## [1] 0

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

## [1] 22 52

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

## [1] 6266

# 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)