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.

• \(W\) is the area of wheat per acres
• \(B\) is the area of barley per acres

Now, we would like to maximize the profit. So, we have to set our objective function as follows:

• The profit of each bushel of wheat(W) is \(\$143\)
• The profit of each bushel of barley(B) is \(\$60\)

\(profit(W,B)=(110)(1.3)W+(30)(2)B\)

\(P=143W+60B\)

The constraints can be set in the following ways:

1. Production Cost

• \(120W+210B\le15000\)

2. Total Area in Acres

• \(W+B\le75\)

3. Storage Space

• \(110W+30B\le4000\)
• \(W\ge0\)
• \(B\ge0\)

Solution with lpSolve

library(lpSolve)

# Set the coefficients of the objective function
obj.func <- c(143, 60)

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

# Create right hand side or the result for the constraints
rhs <- c(15000, 4000, 75, 0, 0)

# Create the direction of the constraints
direc <- c("<=","<=","<=",">=",">=")

# Find the optimum solution
optimum <- lp(direction="max",
              objective.in = obj.func,
              const.mat = constr,
              const.dir = direc,
              const.rhs = 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 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 and B
solution <- optimum$solution
names(solution) <- c("W","B")
print(solution)

##  W  B 
## 22 52

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

## [1] "Total Profit: 6266"

So, the farmer should plant 22 bushels of wheat and 52 bushels of barley to maximize the profit. The profit is $6266.

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

Solution with lpSolveAPI

library(lpSolveAPI)

# Set 5 constraints and 2 decision variables
lprec <- make.lp(nrow = 0,
                 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 obj.func
set.objfn(lprec,
          obj.func)

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

# Display the matrix
lprec

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

# Solve problem
solve(lprec)

## [1] 0

# Get the wheat and barley variables value
get.variables(lprec)

## [1] 22 52

# Get the maximum profit
get.objective(lprec)

## [1] 6266

# Note that the default boundaries on the decision variables 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)

So, the farmer should plant 22 bushels of wheat and 52 bushels of barley to maximize the profit. The profit is $6266.