Find the maximum solution to \[Z=4x+3y\]
Suppose that the objective function is subject to the following constraints: \[\begin{eqnarray} x &\geq& {0}\\ y &\geq& {2}\\ 2y &\leq& {25-x}\\ 4y &\leq& {2x-8}\\ y &\leq& {2x-5} \end{eqnarray}\] Solve these problems using mathematical model and R packages (please explain it step-by-step).
We would like to maximize the solution, so we must set our objective function as follows \[Z = 4x+y = MAX(solution)\] which means that \[\hat {Z} = \begin{pmatrix} 4 \\ 3 \end{pmatrix}\]
Now, We can define the coefficient matrix \[A=\begin{pmatrix}
0 & -1 \\
1 & 2 \\
-2 & 4 \\
-2 & 1
\end{pmatrix}\] and, \[B=\begin{pmatrix}
-2 \\
25 \\
-8 \\
-5
\end{pmatrix}\] Now, we will solving the problem using lpSolve.
lpSolve# install.packages("lpSolve")
library(lpSolve)
# set the coefficients of the objective function
o.func <- c(4, 3)
# set matrix corresponding to coefficients of constraint by rows
# Don't consider the non-negative constraint, it is automatically assumed
cons <- matrix(c( 0,-1,
1, 2,
-2, 4,
-2, 1),
nrow = 4, byrow = T)
# set right hand side for the constraints
rhs <- c(-2, 25, -8, -5)
# set direction of the constraints
c.direct <- c(">=", "<=", "<=", "<=")
# Find the optimal solution
optimum <- lp(direction="max",
objective.in = o.func,
const.mat = cons,
const.dir = c.direct,
const.rhs = rhs,
all.int = T)
str(optimum)## List of 28
## $ direction : int 1
## $ x.count : int 2
## $ objective : num [1:2] 4 3
## $ const.count : int 4
## $ constraints : num [1:4, 1:4] 0 -1 2 -2 1 2 1 25 -2 4 ...
## ..- 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 100
## $ solution : num [1:2] 25 0
## $ 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"
## [1] 0
# Display the optimum values
best_sol <- optimum$solution
names(best_sol) <- c("x", "y")
print(best_sol)## x y
## 25 0
# Check the value of objective function at optimal point
print(paste("Maximum: ", optimum$objval, sep=""))## [1] "Maximum: 100"
So, the maximum solution is 100.
Let say you are working as consultant for a boutique car manufacturer, producing luxury cars. They run on one-month (30 days) cycles, we have one cycle to show we can provide value. There is one robot, 2 engineers and one detailer in the factory. The detailer has some holiday off, so only has 21 days available.
The 2 cars need different time with each resource:
| Resource Time | Car A | Car B |
|---|---|---|
| Robot | 3 days | 4 days |
| Engineer | 5 days | 6 days |
| Detailer | 1.5 days | 3 days |
Car A provides $30,000 profit, whilst Car B offers $45,000 profit. At the moment, they produce 4 of each car per month, for $300,000 profit. Not bad at all, but we think we can do better for them.
We would like to maximize the solution, so we must set our objective function as follows \[P = 30000x+45000y = MAX(profit)\] which means that \[\hat {P} = \begin{pmatrix} 30000 \\ 45000 \end{pmatrix}\]
The constraints can be set in the following ways:
Now, We can define the coefficient matrix \[A=\begin{pmatrix}
3 & 4 \\
5 & 6 \\
1.5& 3
\end{pmatrix}\] and, \[B=\begin{pmatrix}
30 \\
60 \\
21
\end{pmatrix}\] Now, we will solving the problem using lpSolve.
lpSolve# install.packages("lpSolve")
library(lpSolve)
# set the coefficients of the objective function
o.func <- c(30000, 45000)
# set matrix corresponding to coefficients of constraint by rows
# Don't consider the non-negative constraint, it is automatically assumed
cons <- matrix(c( 3, 4,
5, 6,
1.5, 3),
nrow = 3, byrow = T)
# set right hand side for the constraints
rhs <- c(30, 60, 21)
# set direction of the constraints
c.direct <- c("<=", "<=", "<=")
# Find the optimal solution
optimum <- lp(direction="max",
objective.in = o.func,
const.mat = cons,
const.dir = c.direct,
const.rhs = rhs,
all.int = T)
str(optimum)## List of 28
## $ direction : int 1
## $ x.count : int 2
## $ objective : num [1:2] 30000 45000
## $ const.count : int 3
## $ constraints : num [1:4, 1:3] 3 4 1 30 5 6 1 60 1.5 3 ...
## ..- 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 330000
## $ solution : num [1:2] 2 6
## $ 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"
## [1] 0
# Display the optimum values
best_sol <- optimum$solution
names(best_sol) <- c("x", "y")
print(best_sol)## x y
## 2 6
# Check the value of objective function at optimal point
print(paste("Maximum: ", optimum$objval, sep=""))## [1] "Maximum: 330000"
So, the maximum profit is $330000.
Let say you would like to make some sausages and you have the following ingredients available:
| Ingredient | Cost($/kg) | Availability(kg) |
|---|---|---|
| Chicken | 4.32 | 30 |
| Wheat | 2.46 | 20 |
| Starch | 1.86 | 17 |
Assume that you will make 2 types of sausage:
According to government regulations of Indonesia:
So, please figure out how to optimize the cost effectively to blend your sausages.
First, we need to translate the problem in a mathematical way. Let’s define the following variables
Then, we set our objective function as follows \[Cost = 4.32(c_e + c_p) + 2.46(w_e + w_p) + 1.86(s_e + s_p)\]
The constraints can be set in the following ways: \[\begin{eqnarray} c_e + c_p &\leq& {30}\\ w_e + w_p &\leq& {20}\\ s_e + s_p &\leq& {17}\\ c_e &\geq& 0.4 (c_e + w_e + s_e)\\ c_p &\geq& 0.6 (c_p + w_p + s_p)\\ s_c &\geq& 0.25 (c_e + w_e + s_e)\\ s_p &\geq& 0.25 (c_e + w_e + s_e)\\ c_e + w_e + s_e &=& 350 × 0.05\\ c_p + w_p + s_p &=& 500 × 0.05 \end{eqnarray}\]
Please visualize a contour plot of the following function: \[f(x,y)=x^2+y^2+3xy\]
f <- function(x, y) {x^2 + y^2 + 3*x*y}
n <- 30
xpts <- seq(-1.5, 1.5, len = n)
ypts <- seq(-1.5, 1.5, len = n)
gr <- expand.grid(x = xpts, y = ypts)
feval <- with(gr, matrix(f(x, y), nrow = n, ncol = n))
par(mar = c(5, 4, 1, 1))
contour(xpts, ypts, feval, nlevels = 20, xlab = "x", ylab = "y")
points(-1, -1, pch = 19, cex = 2)
abline(h = -1)