Find the maximum solution to \[Z = 4x+3y\] Suppose that the objective function is subject to the following constraints: \[x\geq0 \\ y\geq2 \\ 2y\leq25-x \\ 4y\leq 2x-8 \\ y\leq 2x-5\] Solve these problems using mathematical model and R packages (please explain it step-by-step).
\[Z = 4x+3y= MAX(solution)\] which means that \[\hat Z = \begin{pmatrix} 4 \\ 3 \end{pmatrix}\].
Now, we can define the coefficient matrix \[A = \begin{pmatrix} -1 & 0 \\ 0 & -1 \\ 1 & 2 \\ -2 & 4 \\ -2 & 1 \end{pmatrix}\], and
\[B = \begin{pmatrix} 0 \\ -2 \\ 25 \\ -8 \\ -5 \end{pmatrix}\].
lpSolve## [1] 4 3
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
## [3,] 1 2
## [4,] -2 4
## [5,] -2 1
## [1] 0 2 25 -8 -5
# Direction of the constraints
constranints_direction <- c(">=", ">=", "<=", "<=", "<=")
constranints_direction## [1] ">=" ">=" "<=" "<=" "<="
# Find the optimal solution
optimum <- lp(direction="max",
objective.in = Z,
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] 4 3
## $ const.count : int 5
## $ constraints : num [1:4, 1:5] 1 0 2 0 0 1 2 2 1 2 ...
## ..- 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 90
## $ solution : num [1:2] 21 2
## $ 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 for x, y
best_sol <- optimum$solution
names(best_sol) <- c("x", "y")
print(best_sol)## x y
## 21 2
# Check the value of objective function at optimal point
print(paste("Maksimum Profit: ", optimum$objval, sep=""))## [1] "Maksimum Profit: 90"
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 <- c("Robot","Engineer","Detailer")
car_A <- c(3,5,1.5)
car_B <- c(4,6,3)
car <- data.frame(resource,car_A,car_B)
knitr::kable(car, caption = "Work count in days")| resource | car_A | car_B |
|---|---|---|
| Robot | 3.0 | 4 |
| Engineer | 5.0 | 6 |
| Detailer | 1.5 | 3 |
Car A provides 30,000 dollars profit, which Car B offers 45,000 dollars 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.
First, we need to translate the problem in a mathematical way. Let’s define the following variables
Now we can define \(\hat P = \begin{pmatrix} x \\ y \end{pmatrix}\) as the decision variable vector. Note that it must be \(\hat P \geq 0\). We would like to maximize the profit, so we must set our objective function as follows
\[ P(x, y) = 30000x + 45000y = MAX(profit) \] which means that \[\hat P = \begin{pmatrix} 30000 \\ 45000 \end{pmatrix}\].
The constraints can be set in the following ways:
We can now define the coefficient matrix \[C = \begin{pmatrix} 3 & 4 \\ 5 & 6 \\ 1.5 & 3 \\ -1 & 0 \\ 0 & -1 \end{pmatrix}\], and
\[D = \begin{pmatrix} 30 \\ 60 \\ 21 \\ 0 \\ 0 \end{pmatrix}\].
lpSolveHere are the coefficients of the decision variables:
Therefore, the obj function is:
\[ P(x, y) = 30000x + 45000y = MAX(profit) \]
## [1] 30000 45000
## [,1] [,2]
## [1,] 3.0 4
## [2,] 5.0 6
## [3,] 1.5 3
## [4,] 1.0 0
## [5,] 0.0 1
## [1] 30 60 21 0 0
# Direction of the constraints
constraints_direction <- c("<=", "<=", "<=", ">=", ">=")
constraints_direction## [1] "<=" "<=" "<=" ">=" ">="
# Find the optimal solution
optimum_car <- lp(direction="max",
objective.in = P,
const.mat = C,
const.dir = constraints_direction,
const.rhs = D,
all.int = T)
str(optimum_car)## List of 28
## $ direction : int 1
## $ x.count : int 2
## $ objective : num [1:2] 30000 45000
## $ const.count : int 5
## $ constraints : num [1:4, 1:5] 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 for x, y
best_sol <- optimum_car$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("Maksimum Profit: ", optimum_car$objval, sep=""))## [1] "Maksimum Profit: 330000"
Let say you would like to make some sausages and you have the following ingredients available:
ingredient <- c("Chicken","Wheat","Starch")
cost_dollar <- c(4.32,2.46,1.86)
availability_kg <- c(30,20,17)
sausage <- data.frame(ingredient,cost_dollar,availability_kg)
knitr::kable(sausage, caption = "Sausage recipe in $/kg")| ingredient | cost_dollar | availability_kg |
|---|---|---|
| Chicken | 4.32 | 30 |
| Wheat | 2.46 | 20 |
| Starch | 1.86 | 17 |
Assume that you will make 2 types of sausage: • Economy (>40% Chicken) • Premium (>60% Chicken) • One sausage is 50 grams (0.05 kg)
According to government regulations of Indonesia: • The most starch you can use in your sausages is 25%. • You have a contract with a butcher, and have already purchased 23 kg Chicken, that must go in your sausages. • You have a demand for 350 economy sausages and 500 premium sausages.
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
Now we can define \(\hat Cost = \begin{pmatrix} c \\ w \\ s\end{pmatrix}\) as the decision variable vector. Note that it must be \(\hat Cost \geq 0\). We would like to minimize the cost, so we must 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) = MIN(Cost) \] The constraints can be set in the following ways:
Please visualize a contour plot of the following function: \[f(x,y) = x^2+y^2+3xy\] • How coordinate descent works • How Steepest descent works • How Newton Method works • How Conjugate Gradient works
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)Coordinate descent works is simple. This method can minimize function by successively minimizing each of the individual dimension of the function in a cyclic fashion, while holding the values of the function in the other dimension fixed.
The steepest descent is a simple gradient method for optimization. This method works with a deep slow convergence towards the optimum solution, this happens because of the steps zigzagged. When certain parameters are highly correlated with each other, the steepest descent algorithm can require many steps to reach the minimum. Depending on the starting value, the steepest descent algorithm could take many steps to wind its way towards the minimum.
Newton Method works as a technique for approximating the zero generator of a function (\(F(x) = 0\)). This method uses tangents to approximate the intersection of a function graph with the x-axis.
Conjugate gradient represent a kind of steepest descent approach “with a twist”. With steepest descent, we begin our minimization of a function \(f\) starting at \(x_0\) by traveling in the direction of the negative gradient \(-f'(x_0)\). In subsequent steps, we continue to travel in the direction of the negative gradient evaluated at each successive point until convergence.
The conjugate gradient approach begins in the same manner, but diverges from steepest descent after the first step. In subsequent steps, the direction of travel must be conjugate to the direction most recently traveled.