Optimization week 3

1 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?

1.1 Problem definition

LPSOLVE

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

• \(X\) is the number of acress allotted to wheat.
• \(Y\) is the number of acress allotted to barley.

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

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

\[profit(X,Y)= (110*1.30)X+ (30*2.00)Y = MAX(profit) \] Then, the expression is like \[profit(X,Y)= 143X+ 60Y = MAX(profit) \] The constraints can be set in the following ways like: 1.Cost \[120X +210Y\leq 15000\] 2.Storage space \[110X+30Y\leq4000\] 3.Area \[X+Y\leq 75\] 4.Acres allotted to X or wheat \[X \geq 0 \] 5.Acres allotted to Y or barley \[Y \geq 0 \] which means that \[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}\].

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

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

1.1.1 Objective Function

Here are the coefficients of the decision variables:

• The net profit of \(X\) is \(\1.30\)
• The net profit of \(Y\) is \(\2.00\)

Therefore, the obj function is:

\[profit(X,Y)= (110*1.30)X + (30*2.00)Y\]

\[profit(X,Y)= 143X+ 60Y\]

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

## [1] 143  60

##CONSTRAINT MATRIX The constraints can be set in the following ways like: Row: 1.Cost 2.Storage space 3.Area 4.Acres allotted to X or wheat 5.Acres allotted to Y or barley

# Create constraint martix B
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  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] "<=" "<=" "<=" ">=" ">="

1.1.2 The Optimum Result

# Find the optimal solution
optimum <-  lp(direction="max",
               objective.in = C,
               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 X, Y
best_sol <- optimum$solution
names(best_sol) <- c("X", "Y") 
print(best_sol)

##  X  Y 
## 22 52

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

## [1] "Total profit: 6266"

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

Solve the problem again using lpSolveAPI

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

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

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