Email             :
RPubs            : https://rpubs.com/ardifo/
Jurusan          : Statistika
Instagram      : Ardifosyaa
Address         : ARA Center, Matana University Tower
                         Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.

1 Question 1

Find the Maximum Solution

\[ Z = 4x +3y \]

Suppose that the objective function is subject to the following constraints:

\[ x \geq 0 \\ y \geq 2 \\ 2y \leq 25-x \\ 4y \leq 2x-8 \\ y \leq 2x-5 \\ \] Solve these problems using mathematical model and R packages( please explain it step-by-step)

1.1 Answer

1.1.1 Math Model

  • Problem definition

\[ Z = 4x +3y = MAX(solution) \] which means that

\[ \widehat{Z} = \binom{4}{3} \] 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} \]

1.1.2 Call lpSolve Package

library(lpSolve)                                          # call lpSolve Package
## Warning: package 'lpSolve' was built under R version 4.2.1

1.1.3 Input Coefisient From the Function

fc <- c(4,3)                                              # Input Coefisient from the function

1.1.4 Create the Matrix

fm <- matrix (c(1,0,                                      # Make the Matrix
                0,1,
                1,2,
                -2,4,
                -2,1), nrow=5, byrow=TRUE)
fm
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## [3,]    1    2
## [4,]   -2    4
## [5,]   -2    1
# right for the constraints
fv <- c(0,2,25,-8,-5)                                    
fv 
## [1]  0  2 25 -8 -5

1.1.5 Make the Directions of the Constrains

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

1.1.6 Optimal Solution

# Find the Optimal Solutioin

optimum <- lp(direction="max",
              objective.in = fc, 
              const.mat = fm,
              const.dir = constranints_direction,
              const.rhs = fv,
              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.1.7 Find the Optimum Status

# Find the optimum status 
print(optimum$status)                   # Print status: 0 = Success, 2 = no feasible solution  
## [1] 0

1.1.8 Find the Best Value For \(x\) and \(y\)

best_ASgank <- optimum$solution
names(best_ASgank) <- c("x", "y")          # Display the optimum valued for x, y
print(best_ASgank)
##  x  y 
## 21  2

1.1.9 Find the Total Profit

# Check the value of objective function at optimal point
print(paste("Maximum Solution: ", optimum$objval, sep=""))
## [1] "Maximum Solution: 90"
# Disconnect from the model and the optimum solution
rm(optimum, constranints_direction, best_ASgank)

2 Questions 2

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

2.1 Answer

2.1.1 Math Model

  • Problem Definition

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

  • \(x\) is the number of Cars A
  • \(y\) is the number of Cars B

Now we can define \(\hat{P} = \binom{x}{y}\) as the decision variable vector. Note that it must be \(\hat{P} \geq 0\) . We would li \[ P(x,y) =30000x + 45000y = MAX(profit) \] which means that \[ \hat{P} = \binom{30000}{45000} \]

the constraints can be set in the following ways: \[ 1. \ Robot : 3x + 4y \leq 30\\ 2. \ Engineer : 5x + 6y \leq 60\\ 3. \ Detailer : 1.5x +3y \leq 21\\ 4. \ Car A : x \geq 0\\ 5. \ Car B : y \geq 0 \] Now, we can 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} \]

2.1.2 Call lpSolve Package

library(lpSolve)

2.1.3 Objective Function

Here are the coefficients of the decision variables:

  • The net profit per Car A or \(x\) is \(30000\)
  • The net profit per Car B or \(y\) is \(45000\)

Therefore, the objective function is:

\[ P(x,y) = 30000x + 45000y = MAX(profit) \]

2.1.4 Input Coefficient From the Function

# Set the coefficients of the decision variables -> P
fc2 <- c(30000, 45000)

2.1.5 Create the Matrix

  • First row is for Robot
  • Second row is for Engineer
  • Third row is for Detailer
  • Fourth row is for Car A allotted
  • Fif row is for Car B allotted
fm2 <- matrix(c(3, 4,
              5, 6,
              1.5, 3), nrow=3, byrow=TRUE)
fm2
##      [,1] [,2]
## [1,]  3.0    4
## [2,]  5.0    6
## [3,]  1.5    3
fv2 <- c(30, 60, 21)         # Right hand side for the constraints
fv2
## [1] 30 60 21

2.1.6 Create the Direction of the Constraints

constranints_direction <- c("<=", "<=", "<=")
constranints_direction
## [1] "<=" "<=" "<="

2.1.7 Optimum Solution

# Find the Optimal Solution

optimum_car <- lp(direction="max",
              objective.in = fc2, 
              const.mat = fm2,
              const.dir = constranints_direction,
              const.rhs = fv2,
              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 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"

2.1.8 Find the Optimum Status

# Find the optimum status 
print(optimum_car$status)                   # Print status: 0 = Success, 2 = no feasible solution  
## [1] 0

2.1.9 Find the Best Value for \(x\) and \(y\)

best_sol <- optimum_car$solution
names(best_sol) <- c("x", "y")          # Display the optimum valued for x, y
print(best_sol)
## x y 
## 2 6

2.1.10 Find the Total Profit

# Check the value of objective function at optimal point
print(paste("Maximum Profit: ", optimum_car$objval, sep=""))
## [1] "Maximum Profit: 330000"
# Disconnect from the model and the optimum solution
rm(optimum_car, constranints_direction, best_sol)

3 Questions 3

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:

  • 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.

3.1 Answer

The ingredient that I have

ingredient     <- c("Chicken","Wheat","Starch") 
Cost           <- c(4.32, 2.46, 1.86)
Availability   <- c(30,20,17)
Sausage        <- data.frame(ingredient,Cost,Availability)
Sausage
##   ingredient Cost Availability
## 1    Chicken 4.32           30
## 2      Wheat 2.46           20
## 3     Starch 1.86           17

So, Now we need to make the following variables above in Math way,

  • \(es\) = Ecoomy Sausage
  • \(ps\) = Premium Sausage
  • \(c\) = Chicken
  • \(w\) = Wheat
  • \(s\) = Starch

Now, we can define \(\hat{C}ost=\begin{pmatrix} c \\ w \\ s \end{pmatrix}\) as the decision variable vector. Note that it mus \(\hat{C}ost \leq 0\).

We would like to minimize the cost, so we must set our objective function as follows

\[ MIN(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:

  1.  Chicken : \(c_e + c_p \leq 53\)
  2.  Wheat : \(w_e + w_p \leq 20\)
  3.  Starch : \(S_e + S_p \leq 17\)
  4.  Chicken in economy sausage : \(C_e \geq 0.4(c_e + w_e + s_p)\)
  5.  Chicken in economy sausage : \(c_p \geq 0.6(c_p + w_p + s_p)\)
  6.  Starch in economy sausage : \(s_e \geq 0.25(c_e + w_e + s_e)\)
  7.  Starch in premium sausage : \(s_p \geq 0.25(c_p + w_e + s_p)\)
  8.  Profit economy sausage : \(c_e + w_e + s_e = 350x0.05\)
  9.  Profit premium sausage : \(c_p + w_p + s_p = 500x0.05\)

4 Questions 4

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

4.1 Answer

4.1.1 Contour Plot

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, nlevel = 20, xlab = "x", ylab = "y")
points(-1, -1, pch = 19, cex = 2)
abline(h = -1)

4.1.2 How Coordinate Descent Works

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.

4.1.3 How Steepest Descent Works

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 steeps to reach the minimum. Depending on the starting value, the steepest descent algorithm could take many steps to wind it’s way towards the minimum.

4.1.4 How Newton Method Works

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.

4.1.5 How Conjugate Gradient Works

Conjugate gradient represent a kind of steepest descent approach “with a twist”. with steepest descent, we begin our minimization of function \(f\) starting at \(x_0\) y 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 he same manner, but diverges from steepest descent after the first step. In subsequent steps, the direction most recently traveled.