2 Week3 - IPs, MILPs & Network Design

2.1 Lesson 1 - IPs & MILPs

2.1.1 Banner Chemicals ll

Question Info:assuming you can only produce in lots of 10 barrels.

data <- tibble(Type = c("high","Supreme","Max"),
               Additive_1  = c(3,2,300), Additive_2 = c(1,3,280),
               Profit = c(80,200,NA))

kable(data, caption = "Banner Chemicals ll")

Table 2.1: Banner Chemicals ll
Type	Additive_1	Additive_2	Profit
high	3	1	80
Supreme	2	3	200
Max	300	280

$$ \[\begin{align*} \textbf{Maximize} & \ y = \sum_{i=1}^{N} \frac{x_{i} \cdot p_{i}}{Q} \\ \textbf{subject to}\\ & \sum_{i=1}^{N} x_{i} \le \frac{P_i}{Q} \\ & \sum_{i=1}^{N} x_{i} \cdot a_{i,j} \le \frac{Max_j}{Q}, \ \ j \in M \\ & x_{i} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of product in N} \\ & j = \text{Product component in M} \\ & p_i = \text{Profit for product i} \\ & a_{i,j} = \text{Product requirement for stage j} \\ & P_i = \text{Capacity for Plant} \\ & Max_j = \text{Capacity for components} \\ & Q = \text{Lot} \\ & x_{i} = \text{Quantity of product} \\ \end{align*}\] $$

rm(x)

## Warning in rm(x): objeto 'x' não encontrado

n = 2 #Number of products
m = 2 #Number of components
P = 110 #Plant Capacity

#subsetting data to R
max = as.matrix(data[3,2:3])
a = as.matrix(data[1:2,2:3])
p = as.matrix(data[1:2,4])

Model <- function(Q){

model <- MIPModel() %>%
   #Variable of profit
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr((p[i] * x[i])/Q, i = 1:n), "max") %>%
  
  #Don't produce something you don't have demand for
  add_constraint(sum_expr(x[i], i = 1:n) <= P/Q) %>%
  
  #Dont exceed maximum required
  add_constraint(sum_expr(x[i]*a[i,j], i = 1:n) <= max[j]/Q, j = 1:m)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk", verbose = F))

return(list(Objective = result$objective_value*Q^2,
            Variables = result$solution*Q))
}

#Question 1: Result when Lot = 10

Model(10)

## $Objective
## [1] 18800
## 
## $Variables
## x[1] x[2] 
##   10   90

#Question 2: Result when Lot = 5

Model(5)

## $Objective
## [1] 19000
## 
## $Variables
## x[1] x[2] 
##   25   85

2.1.2 GOnuts

data <- tibble(Type = c("C_Ethiopia","C_Tanzania","C_Nigeria","Demand"),
               Ginko  = c(21,22.5,23,550), Kola = c(22.5,24.5,25.5,450),
               Capacity = c(425,400,750,NA))
kable(data, caption = "GOnuts")

Table 2.2: GOnuts
Type	Ginko	Kola	Capacity
C_Ethiopia	21	22	425
C_Tanzania	22	24	400
C_Nigeria	23	26	750
Demand	550	450

$$ \[\begin{align*} \textbf{Minimize} & \ y = \sum_{i=1}^{N} \sum_{j=1}^{M} c_{i,j} \cdot x_{i,j} \\ \textbf{subject to}\\ & \sum_{i=1}^{N} x_{i} \le C_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i} \cdot a_{i,j} \ge D_i, \ \ j \in M \\ & x_{i} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of product in N} \\ & j = \text{Product component in M} \\ & p_i = \text{Profit for product i} \\ & a_{i,j} = \text{Product requirement for stage j} \\ & P_i = \text{Capacity for Plant} \\ & x_{i} = \text{Quantity of product} \\ \end{align*}\] $$

rm(x)

## Warning in rm(x): objeto 'x' não encontrado

n = 3 #Number of products
m = 2 #Number of components

#subsetting data to R
D = as.matrix(data[4,2:3])
a = as.matrix(data[1:3,2:3])
C = as.matrix(data[1:3,4])

#Modeling for Lesson and question, z is for the exercise, by changing the capacity of plant's constraint

Model <- function(z = NA){
  
  if (is.na(z)) {
   z = C
  } else {
   z = as.matrix(replicate(n = n, z))
  }
  
  #Get back to model
  model <- MIPModel() %>%
   #Variable of productW
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i,j]*a[i,j] ,j = 1:m, i = 1:n), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Attend all demand
  add_constraint(sum_expr(x[i,j], j = 1:m) <= z[i], i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk", verbose = FALSE))
  
return(list(Obj_value = result$objective_value,
            Solution = result$solution))
  
}

#Question 1: Using the capacity of the lesson
Model()

## $Obj_value
## [1] 22638
## 
## $Solution
## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2] 
##      0    425    375     25    175      0

#Question 2: All plants capacity are 350
Model(350)

## $Obj_value
## [1] 22850
## 
## $Solution
## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2] 
##      0    350    250    100    300      0

2.1.3 QQ4 - Two Simple ILPs

2.1.3.1 Part l

data <- tibble(Type = c("Maximize", "C_1","C_2"),
               Var_1 = c(3,5,5), Var_2  = c(4,-3,10), RH = c(NA,8,63))
kable(data)

Type	Var_1	Var_2	RH
Maximize	3	4
C_1	5	-3	8
C_2	5	10	63

rm(x)

## Warning in rm(x): objeto 'x' não encontrado

n = 2 #Number of Variables & constraints

#subsetting data to R
a <- as.matrix(data[1,2:3])
c <- as.matrix(data[2:3,2:3])
RH <- as.matrix(data[2:3,4])

#Get back to model
model <- MIPModel() %>%
   #Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[i,j], j = 1:n) <= RH[i], i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 28

#result solution
result$solution

## x[1] x[2] 
##    4    4

2.1.3.2 Part 2

data <- tibble(Type = c("Maximize", "C_1","C_2","C_3"),
               Var_1 = c(2,-2,2,-2), Var_2  = c(3,1,0,2),
               State = c(NA,"<=","<=",">="), RH = c(NA,4,4,-2))
kable(data)

Type	Var_1	Var_2	State	RH
Maximize	2	3
C_1	-2	1	<=	4
C_2	2	0	<=	4
C_3	-2	2	>=	-2

rm(x)

n = 2 #Number of Variables & constraints


#subsetting data to R
a <- as.matrix(data[1,2:3])
c <- as.matrix(data[2:4,2:3])
RH <- as.matrix(data[2:4,5])


#Get back to model
model <- MIPModel() %>%
   #Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[j,i], i = 1:n) <= RH[j], j = 1:n) %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[j,i], i = 1:n) >= RH[j], j = 3)


#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 28

#result solution
result$solution

## x[1] x[2] 
##    2    8

2.1.4 GOnuts ll

data <- tibble(Type = c("C_Ethiopia","C_Tanzania","C_Nigeria","Demand"),
               Ginko  = c(21,22.5,23,550), Kola = c(22.5,24.5,25.5,450),
               Capacity = c(425,400,750,NA), Plant_Cost = c(1500,2000,3000,NA))

kable(data, caption = "GOnuts ll")

Table 2.3: GOnuts ll
Type	Ginko	Kola	Capacity	Plant_Cost
C_Ethiopia	21	22	425	1500
C_Tanzania	22	24	400	2000
C_Nigeria	23	26	750	3000
Demand	550	450

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{N} \sum_{j=1}^{M} c_{i,j} \cdot x_{i,j} + \sum_{j=1}^{M} f_j \cdot y_j \\ \textbf{subject to} \\ & \sum_{i=1}^{N} x_{i} \le C_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i} \cdot a_{i,j} \ge D_i, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} - f_j \cdot B_M \le 0\\ & x_{i} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Demand in M} \\ & p_i = \text{Profit for product i} \\ & a_{i,j} = \text{Product requirement for stage j} \\ & P_j = \text{Capacity for Plant} \\ & B_M = \text{Big number}\\ & Max_j = \text{Capacity for components} \\ & x_{i} = \text{Quantity of product} \\ \end{align*}\] $$

rm(x)

## Warning in rm(x): objeto 'x' não encontrado

n = 3 #Number of products
m = 2 #Number of components

#subsetting data to R
D = as.matrix(data[4,2:3])
a = as.matrix(data[1:3,2:3])
C = as.matrix(data[1:3,4])
Pc = as.matrix(data[1:3,5])
B = sum(D) 

Model <- function(z = NA){
  
  if (is.na(z)) {
   z = C
  } else {
   z = as.matrix(replicate(n = n, z))
  }

#Get back to model
model <- MIPModel() %>%
   #Variable of production
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%
  
  #Binary Variable
  add_variable(f[i], i = 1:n, type = "binary") %>%

   # maximize profit
  set_objective(sum_expr(x[i,j]*a[i,j] ,j = 1:m, i = 1:n) +
                sum_expr(f[i]*Pc[i], i = 1:n), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Attend all demand
  add_constraint(sum_expr(x[i,j], j = 1:m) <= z[i], i = 1:n) %>%
  
  #
  add_constraint(sum_expr(x[i,j], j = 1:m) - f[i]* B <= 0, i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk", verbose = FALSE))

return(list(Obj_value = result$objective_value,
            Solution = result$solution))

}

#Question 1: Using the capacity of the lesson
Model()

## $Obj_value
## [1] 27350
## 
## $Solution
## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2]   f[1]   f[2]   f[3] 
##      0    425      0      0    550     25      1      0      1

#Question 2: All plants capacity are 450
Model(450)

## $Obj_value
## [1] 29050
## 
## $Solution
## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2]   f[1]   f[2]   f[3] 
##      0    450    450      0    100      0      1      1      1

2.1.5 GOnuts lll

data <- tibble(Type = c("C_Ethiopia","C_Tanzania","C_Nigeria","Demand"),
               Plant_Cost = c(1500,2000,3000,NA),
               Ginko  = c(21,22.5,23,550), Kola = c(22.5,24.5,25.5,450),
               Min_Cap = c(100,250,600,NA), Max_Cap = c(425,400,750,NA))

kable(data, caption = "GOnuts lll")

Table 2.4: GOnuts lll
Type	Plant_Cost	Ginko	Kola	Min_Cap	Max_Cap
C_Ethiopia	1500	21	22	100	425
C_Tanzania	2000	22	24	250	400
C_Nigeria	3000	23	26	600	750
Demand		550	450

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{N} \sum_{j=1}^{M} c_{i,j} \cdot x_{i,j} + \sum_{j=1}^{M} f_j \cdot y_j \\ \textbf{subject to} \\ & \sum_{i=1}^{N} x_{i,j} \le C_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} \ge L_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} \cdot a_{i,j} \ge D_i, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} - f_j \cdot B_M \le 0\\ & x_{i} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Demand in M} \\ & p_i = \text{Profit for product i} \\ & a_{i,j} = \text{Product requirement for stage j} \\ & P_j = \text{Capacity for Plant} \\ & B_M = \text{Big number}\\ & Max_j = \text{Capacity for components} \\ & x_{i} = \text{Quantity of product} \\ \end{align*}\] $$

rm(x)

n = 3 #Number of products
m = 2 #Number of components

#subsetting data to R
D = as.matrix(data[4,3:4])
a = as.matrix(data[1:3,3:4])
L = as.matrix(data[1:3,5])
C = as.matrix(data[1:3,6])
Pc = as.matrix(data[1:3,2])
B = sum(D)

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%
  
  #Binary Variable
  add_variable(f[i], i = 1:n, type = "binary") %>%

  # maximize profit
  set_objective(sum_expr(x[i,j]*a[i,j] ,j = 1:m, i = 1:n) +
                sum_expr(f[i]*Pc[i], i = 1:n), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Attend all demand
  add_constraint(sum_expr(x[i,j], j = 1:m) <= C[i], i = 1:n) %>%
  
  #
  add_constraint(sum_expr(x[i,j], j = 1:m) >= L[i]*f[i], i = 1:n) %>%
  
  #
  add_constraint(sum_expr(x[i,j], j = 1:m) - f[i]* B <= 0, i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 27425

#result solution
result$solution

## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2]   f[1]   f[2]   f[3] 
##      0    400      0      0    550     50      1      0      1

2.1.6 GOnuts lV

Operate with 2 Plants at maximum constraint.

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{N} \sum_{j=1}^{M} c_{i,j} \cdot x_{i,j} + \sum_{j=1}^{M} f_j \cdot y_j \\ \textbf{subject to} \\ & \sum_{i=1}^{N} x_{i,j} \le C_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} \ge L_j, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} \cdot a_{i,j} \ge D_i, \ \ j \in M \\ & \sum_{i=1}^{N} x_{i,j} - f_j \cdot B_M \le 0\\ & \sum_{j=1}^{M} f_j \le K \\ & x_{i} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Demand in M} \\ & p_i = \text{Profit for product i} \\ & a_{i,j} = \text{Product requirement for stage j} \\ & P_j = \text{Capacity for Plant} \\ & B_M = \text{Big number}\\ & Max_j = \text{Capacity for components} \\ & x_{i} = \text{Quantity of product} \\ & K = \text{Maximum of Plants to be open} \\ \end{align*}\] $$

rm(x)

## Warning in rm(x): objeto 'x' não encontrado

n = 3 #Number of products
m = 2 #Number of components
K = 2

#subsetting data to R
D = as.matrix(data[4,3:4])
a = as.matrix(data[1:3,3:4])
L = as.matrix(data[1:3,5])
C = as.matrix(data[1:3,6])
Pc = as.matrix(data[1:3,2])
B = sum(D)

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%
  
  #Binary Variable
  add_variable(f[i], i = 1:n, type = "binary") %>%

  # maximize profit
  set_objective(sum_expr(x[i,j]*a[i,j] ,j = 1:m, i = 1:n) +
                sum_expr(f[i]*Pc[i], i = 1:n), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Attend all demand
  add_constraint(sum_expr(x[i,j], j = 1:m) <= C[i], i = 1:n) %>%
  
  #TODO
  add_constraint(sum_expr(f[j], j = 1:m) <= K) %>%
  
  #TODO
  add_constraint(sum_expr(x[i,j], j = 1:m) >= L[i]*f[i], i = 1:n) %>%
  
  #TODO
  add_constraint(sum_expr(x[i,j], j = 1:m) - f[i]* B <= 0, i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 27425

#result solution
result$solution

## x[1,1] x[1,2] x[2,1] x[2,2] x[3,1] x[3,2]   f[1]   f[2]   f[3] 
##      0    400      0      0    550     50      1      0      1

2.2 Recitations

2.2.1 Integer Linear Programming

data <- tibble(Type = c("Maximize", "C_1","C_2","C_3"),
               Var_1 = c(8,1,3,1), Var_2  = c(20,1,2,3), RH = c(NA,11,30,28))

kable(data)

Type	Var_1	Var_2	RH
Maximize	8	20
C_1	1	1	11
C_2	3	2	30
C_3	1	3	28

rm(x)

n = 2 #Number of Variables
m = 3 #Number of constraints

#subsetting data to R
a <- as.matrix(data[1,2:3])
c <- as.matrix(data[2:4,2:3])
RH <- as.matrix(data[2:4,4])

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%

   # maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[j,i], i = 1:n) <= RH[j], j = 1:m)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 188

#result solution
result$solution

## x[1] x[2] 
##    1    9

2.2.2 Capital Budgeting Problem

data <- tibble(Project = c("P1","P2","P3","P4","P5","Funds"),
               Exp_1 = c(5,4,3,7,8,25), Exp_2 = c(1,7,9,4,6,25),
               Exp_3 = c(8,10,2,1,10,25), Returns = c(20,40,20,15,30,NA))


kable(data, caption = "Capital Budgeting Problem")

Table 2.5: Capital Budgeting Problem
Project	Exp_1	Exp_2	Exp_3	Returns
P1	5	1	8	20
P2	4	7	10	40
P3	3	9	2	20
P4	7	4	1	15
P5	8	6	10	30
Funds	25	25	25

rm(x)

n = 5 #Number of Projects
m = 3 #Number of Years

#subsetting data to R
a <- as.matrix(data[1:5,5])
c <- as.matrix(data[1:5,2:4])
RH <- as.matrix(data[6,2:4])

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i], i = 1:n, type = "binary", lb = 0) %>%

  # maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[i,j], i = 1:n) <= RH[j], j = 1:m)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 95

#result solution
result$solution

## x[1] x[2] x[3] x[4] x[5] 
##    1    1    1    1    0

2.2.3 Fixed Charge Problem

data <- tibble(Type = c("Shirts","Short","Pants","RH"),
               Labor = c(3,2,6,150), Cloth = c(4,3,4,160),
               Sales = c(12,8,15,NA), Fix_Cost = c(200,150,100,NA),
               Var_Cost = c(6,4,8,NA))

kable(data, caption = "Fixed Charge Problem")

Table 2.6: Fixed Charge Problem
Type	Labor	Cloth	Sales	Fix_Cost	Var_Cost
Shirts	3	4	12	200	6
Short	2	3	8	150	4
Pants	6	4	15	100	8
RH	150	160

rm(x)

n = 3 #Number of Products

#subsetting data to R
a <- as.matrix(data[1:3,2:3])
RH <- as.matrix(data[4,2:3])
p <- as.matrix(data[1:3,4])
c <- as.matrix(data[1:3,5:6])
M <- 100

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%
  
  # Variable of fixed cost
  add_variable(y[i], i = 1:n, type = "binary", lb = 0) %>%

  # maximize profit
  set_objective(sum_expr(x[i]*(p[i]-c[i,2]) - y[i]*c[i,1], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*a[i,j], i = 1:n) <= RH[j], j = 1:2) %>%
  
  #Flow Constraint
  add_constraint(x[i] <= y[i]*M, i = 1:n)


#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 75

#result solution
result$solution

## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0    0   25    0    0    1

2.2.4 Minimum Amount Problem

data <- tibble(Type = c("Steel","Labor","Profit"),
               Compact = c(1.5,30,2000), Midsize = c(3,25,3000),
               Large = c(5,40,4000), Capacity = c(6000,60000,NA))

kable(data, caption = "Minimum Amount Problem")

Table 2.7: Minimum Amount Problem
Type	Compact	Midsize	Large	Capacity
Steel	1.5	3	5	6000
Labor	30.0	25	40	60000
Profit	2000.0	3000	4000

rm(x)

n = 3 #Number of cars
m = 2 #Number of components (labor & Steel)

#subsetting data to R
a <- as.matrix(data[1:2,2:4])
RH <- as.matrix(data[1:2,5])
p <- as.matrix(data[3,2:4])
M <- 10000
Min <- 1000

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%
  
  # Variable of minimum amount condition
  add_variable(y[i], i = 1:n, type = "binary", lb = 0) %>%

  # maximize profit
  set_objective(sum_expr(x[i]*p[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*a[j,i], i = 1:n) <= RH[j], j = 1:m) %>%
  
  #Flow Constraint
  add_constraint(x[i] <= M*y[i], i = 1:n) %>%
  
  #Min quantity
  add_constraint(x[i] >= Min*y[i], i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 6000000

#result solution
result$solution

## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0 2000    0    0    1    0

2.2.5 Transportation Problem

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{m} \sum_{i=1}^{N} c_{i,j} \cdot x_{i,j} \\ \textbf{subject to} \\ & \sum_{i=1}^{N} x_{i,j} \ge D_j, \ \ j \in M \\ & \sum_{i=1}^{M} x_{i,j} \le S_i \ \ i \in N \\ & x_{i,j} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Demand in M} \\ & c_{i,j} = \text{Transportation cost from Plant to Demand} \\ & x_{i,j} = \text{Quantity of product} \\ & D_j = \text{Demand of cities} \\ & S_j = \text{Plant Capacity} \end{align*}\] $$

data <- tibble(From = c("Plant 1","Plant 2","Plant 3", "Demand"),
               City_1 = c(8,9,14,45), City_2 = c(6,12,9,20),
               City_3 = c(10,13,16,30), City_4 = c(9,7,5,30),
               Supply = c(35,50,40,NA))

kable(data, caption = "Transportation Problem")

Table 2.8: Transportation Problem
From	City_1	City_2	City_3	City_4	Supply
Plant 1	8	6	10	9	35
Plant 2	9	12	13	7	50
Plant 3	14	9	16	5	40
Demand	45	20	30	30

rm(x)

n = 3 #Number of Plants
m = 4 #Number of cities

#subsetting data to R
c <- as.matrix(data[1:3,2:5])
S <- as.matrix(data[1:3,6])
D <- as.matrix(data[4,2:5])


#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%

  # maximize profit
  set_objective(sum_expr(x[i,j]*c[i,j], i = 1:n, j = 1:m), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Flow Constraint
  add_constraint(sum_expr(x[i,j], j = 1:m) <= S[i], i = 1:n)


#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))


#result solution
Result <- cbind(data[1:n,1],
                matrix(result$solution, nrow = n, ncol = m,
                dimnames = list(NULL, colnames(data[,2:5])))) %>%
          mutate(Supplied = colSums(data[1:n,2:m]), Supply = data$Supply[1:n]) %>%
          bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

kable(Result, caption = paste0("Optimal Value: ", result$objective_value))

Table 2.9: Optimal Value: 1020
From	City_1	City_2	City_3	City_4	Supplied	Supply
Plant 1	0	0	5	10	31	35
Plant 2	10	45	0	0	27	50
Plant 3	25	0	0	30	39	40
Total	35	45	5	40	97	125

2.2.6 Transhipment Problem

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{M} \sum_{i=1}^{N} c_{i,j} \cdot x_{i,j} + \sum_{i=1}^{N} c_{j,k} \cdot x_{j,k} \\ \textbf{subject to} \\ & \sum_{i=1}^{M} x_{i,j} \le S_i, \ \ j \forall M \\ & \sum_{i=1}^{N} x_{j,k} \ge D_j, \ \ k \forall K \\ & \sum_{i=1}^{N} x_{i,j} \le T_j, \ \ j \forall M \\ & \sum_{i=1}^{N} x_{i,j} - \sum_{j=1}^{M} x_{j,k} = 0, \ \ k \forall K \\ & x_{i,j,k} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Transhipment in M} \\ & k = \text{Demand in K} \\ & c_{i,j} = \text{Inflow Transportation cost} \\ & c_{j,k} = \text{Outflow Transportation cost} \\ & x_{i,j} = \text{inflow of product} \\ & x_{j,k} = \text{Outflow of product} \\ & D_k = \text{Demand of cities} \\ & S_i = \text{Plant Capacity} \\ & T_j = \text{Transhipment Capacity} \end{align*}\] $$

#TODO Might be an error on the lecture slide W3-Advanced Optimization on slide 16, the Transshipment Problem. It states the flow constraint as \sum_i x_{i,j} = \sum_i x_{j,i}. This would mean the distance from the plant to CD is the same from the CD to POS. It seems it should be \sum_i x_{i,j} = \sum_j x_{j,k} where k is the node of demand.

#TODO This could be simpler

#Variable Names
Ps <- c("P_1","P_2","P_3","P_4","P_5")
Tsn <- c("TS_1","TS_2")
Tso <- c("TO_1","TO_2")
Dn <- c("D_1","D_2","D_3","D_4","D_5")

#Values
S <- matrix(c(200,300,100,150,220), nrow = 5, dimnames = list(Ps,"Plant_Cap"))
Ts <- matrix(c(450,300), nrow = 1, dimnames = list("Transh_Cap",Tsn))
D <- matrix(c(150,100,110,200,180), nrow = 5, dimnames = list(Dn,"Demand"))
inc <- matrix(c(30,50,23,66,35,14,70,12,65,70), nrow = 5, ncol = 2, byrow = T,
              dimnames = list(Ps,Tsn))
outc <- matrix(c(12,25,22,40,41,65,22,23,12,15), nrow = 2, ncol = 5, byrow = T,
               dimnames = list(Tso,Dn))


#Visualizing dataset
data <- as.tibble(S) %>%
        cbind(inc,t(outc),D)

kable(data, caption = "Transhipment Problem")

Table 2.10: Transhipment Problem
	Plant_Cap	TS_1	TS_2	TO_1	TO_2	Demand
P_1	200	30	50	12	65	150
P_2	300	23	66	25	22	100
P_3	100	35	14	22	23	110
P_4	150	70	12	40	12	200
P_5	220	65	70	41	15	180

rm(x)

n = 5 #Number of Plants
t = 2 #Number of transhipments (CDs)
m = 5 #Number of cities (POS)

#Get back to model
model <- MIPModel() %>%
   #Variable of inflow
  add_variable(x[i,j], i = 1:n, j = 1:t, type = "integer", lb = 0) %>%
  
   #Variable of outflow
  add_variable(u[j,k], j = 1:t, k = 1:m, type = "integer", lb = 0) %>%

   #minimize distance cost
  set_objective(sum_expr(x[i,j]*inc[i,j],  i = 1:n, j = 1:t) + 
                sum_expr(u[j,k]*outc[j,k], j = 1:t, k = 1:m), "min") %>%
  
  #Supply Constraint
  add_constraint(sum_expr(x[i,j], j = 1:t) <= S[i], i = 1:n) %>%
  
  #Attend all demand
  add_constraint(sum_expr(u[j,k], j = 1:t) >= D[k], k = 1:m) %>%
  
  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) <= Ts[j], j = 1:t) %>%
  
  #Flow constraint
  add_constraint(sum_expr(x[i,j], i = 1:n) - sum_expr(u[j,k], k = 1:m) == 0,
                 j = 1:t)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))


#Result

#Inflow
a <- as.tibble(matrix(get_solution(result, x[i,j])$value,
              nrow = n, ncol = t, dimnames = list(Ps,Tsn))) %>%
     mutate(Plant = Ps, Production = rowSums(.)) %>%
     select(Plant, -Production,  everything(), Production) %>%
     cbind(Plant_Cap = data$Plant_Cap) %>%
     bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

kable(a, rownames = F, caption = paste0("Optimal Value: ", result$objective_value))

Table 2.11: Optimal Value: 30700
Plant	TS_1	TS_2	Production	Plant_Cap
P_1	140	50	190	200
P_2	300	0	300	300
P_3	0	100	100	100
P_4	0	150	150	150
P_5	0	0	0	220
Total	440	300	740	970

#Outflow
b <- as.tibble(matrix(get_solution(result, u[j,k])$value, nrow = t, ncol = m,
                        dimnames = list(Tso,Dn))) %>%
     mutate(DC = Tso, demand = rowSums(.), DC_Cap = Ts) %>%
     select(DC,-demand, everything(), demand) %>%
     bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

kable(b, row.names = F, caption = paste0("Optimal Value: ", result$objective_value))

Table 2.11: Optimal Value: 30700
DC	D_1	D_2	D_3	D_4	D_5	demand	DC_Cap
TO_1	150	100	110	0	80	440	450
TO_2	0	0	0	200	100	300	450
Total	150	100	110	200	180	740	1500

2.3 Practice Problems

2.3.1 Production Scheduling

data <- tibble(Type = c("Maximize", "Electrician","Mechanic"),
               Prod_X = c(1.5,1.2,0.5), Prod_Y  = c(2,1,1.6), RH = c(NA,7.6,6.5))

kable(data, caption = "Production Scheduling")

Table 2.12: Production Scheduling
Type	Prod_X	Prod_Y	RH
Maximize	1.5	2.0
Electrician	1.2	1.0	7.6
Mechanic	0.5	1.6	6.5

rm(x)

n = 2 #Number of Products
m = 2 #Number of constraints

#subsetting data to R
a <- as.matrix(data[1,2:3])
c <- as.matrix(data[2:3,2:3])
RH <- as.matrix(data[2:3,4])

#Get back to model
model <- MIPModel() %>%
  
   #Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[j,i], i = 1:n) <= RH[j], j = 1:m)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 10

#result solution
result$solution

## x[1] x[2] 
##    3    3

2.3.2 Capital Budgeting Problem

data <- tibble(Project = c("P1","P2","P3","P4","P5","Funds"),
               Exp_1 = c(12,7,9,2,4,25), Exp_2 = c(8,5,6,4,6,20),
               Exp_3 = c(4,1,3,8,10,20), Returns = c(30,20,25,23,27,NA))

kable(data, caption = "Capital Budgeting Problem")

Table 2.13: Capital Budgeting Problem
Project	Exp_1	Exp_2	Exp_3	Returns
P1	12	8	4	30
P2	7	5	1	20
P3	9	6	3	25
P4	2	4	8	23
P5	4	6	10	27
Funds	25	20	20

rm(x)

n = 5 #Number of Projects
m = 3 #Number of Years

#subsetting data to R
a <- as.matrix(data[1:5,5])
c <- as.matrix(data[1:5,2:4])
RH <- as.matrix(data[6,2:4])

#Get back to model
model <- MIPModel() %>%
   # Variable of production
  add_variable(x[i], i = 1:n, type = "binary", lb = 0) %>%

   # maximize profit
  set_objective(sum_expr(x[i]*a[i], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*c[i,j], i = 1:n) <= RH[j], j = 1:m)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 82

#result solution
result$solution

## x[1] x[2] x[3] x[4] x[5] 
##    1    0    1    0    1

2.3.3 Fixed Charge Problem

data <- tibble(Type = c("Homer","Turner","Sargent","RH"),
               Labor = c(3,2,6,150), Glass = c(4,3,4,160),
               Sales = c(12,8,15,NA), Fix_Cost = c(200,150,100,NA),
               Var_Cost = c(6,4,8,NA))

kable(data, caption = "Fixed Charge Problem")

Table 2.14: Fixed Charge Problem
Type	Labor	Glass	Sales	Fix_Cost	Var_Cost
Homer	3	4	12	200	6
Turner	2	3	8	150	4
Sargent	6	4	15	100	8
RH	150	160

rm(x,y)

n = 3 #Number of Products

#subsetting data to R
a <- as.matrix(data[1:3,2:3])
RH <- as.matrix(data[4,2:3])
p <- as.matrix(data[1:3,4])
c <- as.matrix(data[1:3,5:6])
M <- 100

#Get back to model
model <- MIPModel() %>%
  # Variable of production
  add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%
  
   #Variable of fixed cost
  add_variable(y[i], i = 1:n, type = "binary", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i]*(p[i]-c[i,2]) - y[i]*c[i,1], i = 1:n), "max") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i]*a[i,j], i = 1:n) <= RH[j], j = 1:2) %>%
  
  #Flow Constraint
  add_constraint(x[i] <= y[i]*M, i = 1:n)


#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

#Optimal Value
result$objective_value

## [1] 75

#result solution
result$solution

## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0    0   25    0    0    1

2.3.4 Minimums and Maximums

This exercise has some of changing conditions that we can set as a funcion: -Component Capacity -Minimum Demand -Maximum Demand -Fixed Cost

data <- tibble(Type = c("Labor","Balloons","Cost","Price"),
               Standard = c(3,2,2,3), Joyful = c(5,5,3,5),
               Fabulous = c(8,8,4,7))


kable(data, caption = "Minimums and Maximums")

Table 2.15: Minimums and Maximums
Type	Standard	Joyful	Fabulous
Labor	3	5	8
Balloons	2	5	8
Cost	2	3	4
Price	3	5	7

#TODO Resultados não estão batendo

rm(x)

n = 3 #Number of products
m = 2 #Number of components (labor & Baloons)

#Fixed data
a <- as.matrix(data[1:2,2:4])
c <- as.matrix(data[3,2:4])
p <- as.matrix(data[4,2:4])


 #Variables of questions index
 # Component Capacity = RH 1:m
 # Minimum Demand = Min 1:n
 # Maximum Demand = Max  1:n
 # Fixed Cost = Fc 1:n

#Get back to model
Model <- function(RH = NA, Min = NA, Max = NA, Fc = NA ){
  


  if (is.na(RH))  {RH = as.matrix(replicate(n = m, 9999))} else {RH = as.matrix(RH)}
  if (is.na(Min)) {Min = as.matrix(replicate(n = n, 0))} else {Min = as.matrix(Min)}
  if (is.na(Max)) {Max = as.matrix(replicate(n = n, 9999))} else {Max = as.matrix(Max)}
  if (is.na(Fc))  {Fc = as.matrix(replicate(n = n, 0))} else {Fc = as.matrix(Fc)}
  
  
  

model <- MIPModel() %>%
  
         #Variable of production
        add_variable(x[i], i = 1:n, type = "integer", lb = 0) %>%
        
         #Variable of minimum amount condition
        add_variable(y[i], i = 1:n, type = "binary", lb = 0) %>%
      
         #maximize profit
        set_objective(sum_expr(x[i]*(p[i]-c[i]) - y[i]*Fc[i], i = 1:n),"max") %>%
      
        #Attend all demand
        add_constraint(sum_expr(x[i]*a[j,i], i = 1:n) <= RH[j], j = 1:m) %>%
        
        #Min quantity
        add_constraint(x[i] >= Min[i]*y[i], i = 1:n) %>%
          
        #Max quantity
        add_constraint(x[i] <= Max[i]*y[i], i = 1:n)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))

return(list(Solution = result$solution,
            Obj_Val = result$objective_value))
  
}


Model(c(325,400))

## $Solution
## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0   65    0    0    1    0 
## 
## $Obj_Val
## [1] 130

Model(RH = c(350,400),Min = c(10,10,10))

## $Solution
## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0   70    0    0    1    0 
## 
## $Obj_Val
## [1] 140

Model(RH = c(350,400),Min = c(10,10,10), Max = c(40,40,40))

## $Solution
## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0   38   20    0    1    1 
## 
## $Obj_Val
## [1] 136

Model(RH = c(350,400),Min = c(10,10,10), Max = c(40,40,40), Fc = c(10,5,1))

## $Solution
## x[1] x[2] x[3] y[1] y[2] y[3] 
##    0   38   20    0    1    1 
## 
## $Obj_Val
## [1] 130

2.3.5 Transportation Problem

Formulation already seen on Recitations

data <- tibble(From = c("Chicago","Atlanta","Denver", "Demand"),
               Boston = c(1.04,1.23,1.92,11000), Seattle = c(1.27,1.93,0.94,6300),
               Tampa = c(1.22,0.60,1.03,7400), Supply = c(10000,10000,10000, NA))

kable(data, caption = "Transportation Problem")

Table 2.16: Transportation Problem
From	Boston	Seattle	Tampa	Supply
Chicago	1.0	1.27	1.2	10000
Atlanta	1.2	1.93	0.6	10000
Denver	1.9	0.94	1.0	10000
Demand	11000.0	6300.00	7400.0

rm(x)

n = 3 #Number of Plants
m = 3 #Number of cities

#subsetting data to R
c <- as.matrix(data[1:3,2:4])
S <- as.matrix(data[1:3,5])
D <- as.matrix(data[4,2:4])


#Get back to model
model <- MIPModel() %>%
   #Variable of production
  add_variable(x[i,j], i = 1:n, j = 1:m, type = "integer", lb = 0) %>%

   #maximize profit
  set_objective(sum_expr(x[i,j]*c[i,j], i = 1:n, j = 1:m), "min") %>%

  #Attend all demand
  add_constraint(sum_expr(x[i,j], i = 1:n) >= D[j], j = 1:m) %>%
  
  #Flow Constraint
  add_constraint(sum_expr(x[i,j], j = 1:m) <= S[i],  i = 1:n)


#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))


#result solution
Result <- cbind(data[1:n,1],
                matrix(result$solution, nrow = n, ncol = m,
                dimnames = list(NULL, colnames(data[,2:4])))) %>%
          mutate(Supplied = colSums(data[1:n,2:4]), Supply = data$Supply[1:n]) %>%
          bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))


kable(Result, caption = paste0("Optimal Value: ", result$objective_value))

Table 2.17: Optimal Value: 21992
From	Boston	Seattle	Tampa	Supplied	Supply
Chicago	10000	1000	0	4.2	10000
Atlanta	0	0	6300	4.1	10000
Denver	0	7400	0	2.8	10000
Total	10000	8400	6300	11.2	30000

2.3.6 Transshipment Problem

Exercise exactly the same as the one on Recitations

2.3.7 Facility Problem

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{N} \sum_{j=1}^{M} c_{i} \cdot x_{i} \\ \textbf{subject to} \\ & \sum_{i=1}^{N} x_i = 2 \\ & x_{i} \ge 0, \ x \in \mathbb{B} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & x_{i} = \text{Binary decision} \\ & c_i = \text{Cost to transport materials} \\ \end{align*}\] $$

data <- tibble(Type = c("WH1","WH2","WH3","WH4","WH5"),
               Total_cost = c(116200,98400,107800,86000,104800))

kable(data, caption = "Facility Problem")

Table 2.18: Facility Problem
Type	Total_cost
WH1	116200
WH2	98400
WH3	107800
WH4	86000
WH5	104800

rm(x)

n = 5 #Number of WH

#subsetting data to R
c = as.matrix(data[,2])


#Get back to model
model <- MIPModel() %>%

         #Variable of production
        add_variable(x[i], i = 1:n, type = "binary", lb = 0) %>%
      
         #maximize profit
        set_objective(sum_expr(x[i]*c[i], i = 1:n), "min") %>%
      
        #Open only 2 WH
        add_constraint(sum_expr(x[i], i = 1:n) == 2)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk", verbose = FALSE))


#Results
result$objective_value

## [1] 184400

result$solution## Transhipment Problem

## x[1] x[2] x[3] x[4] x[5] 
##    0    1    0    1    0

This problem was called as facilty location problem, but the problem is mainly based on transhipment without CD capacity, so the lowest total cost CD will be choosen

$$ \[\begin{align*} \textbf{Minimize} & \ z = \sum_{i=1}^{M} \sum_{i=1}^{N} c_{i,j} \cdot x_{i,j} + \sum_{j=1}^{N} \sum_{k=1}^{K} c_{j,k} \cdot x_{j,k} + \sum_{i=1}^{N} y_j \cdot f_y \\ \textbf{subject to} \\ & \sum_{i=1}^{M} x_{i,j} \le S_i, \ \ j \forall M \\ & \sum_{i=1}^{N} x_{j,k} \ge D_j, \ \ k \forall K \\ & \sum_{i=1}^{N} x_{i,j} \le T_j, \ \ j \forall M \\ & \sum_{i=1}^{N} x_{i,j} - \sum_{j=1}^{M} x_{j,k} = 0, \ \ k \forall K \\ & x_{i,j,k} \ge 0, \ x \in \mathbb{Z} \\ \textbf{Where} \\ & i = \text{Type of Plant in N} \\ & j = \text{Transhipment in M} \\ & k = \text{Demand in K} \\ & c_{i,j} = \text{Inflow Transportation cost} \\ & c_{j,k} = \text{Outflow Transportation cost} \\ & x_{i,j} = \text{inflow of product} \\ & x_{j,k} = \text{Outflow of product} \\ & D_k = \text{Demand of cities} \\ & S_i = \text{Plant Capacity} \\ & T_j = \text{Transhipment Capacity} \end{align*}\] $$

data <- WH_Cst %>%
        inner_join(INB_Cost, by = c("OUT tsp" = "INB tsp")) %>%
        inner_join(Out_cost, by = c("OUT tsp" = "OUT tsp")) %>%
        bind_rows(Demand) %>%
        mutate(`OUT tsp`=replace(`OUT tsp`, is.na(`OUT tsp`), "Demand")) %>%
        rename("Type" = `OUT tsp`)

rm(INB_Cost,OUT,WH_Cst,Out_cost,Demand)

kable(data, caption = "Facility Location Problem") %>%
  kable_styling() %>%
  scroll_box(width = "100%")

Table 2.19: Facility Location Problem
Type	Variable costs	Fixed costs	Plant 1	Plant 2	RW 1	RW 2	RW 3	RW 4	RW 5	RW 6	RW 7	RW 8	RW 9
CW 1	8	100000	77	54	89	78	69	95	85	73	88	73	86
CW 2	10	90000	48	33	90	91	63	61	94	58	89	90	55
CW 3	25	80000	41	21	81	82	71	98	62	71	85	75	72
CW 4	30	70000	57	94	73	97	55	75	86	54	62	99	59
CW 5	50	60000	94	44	85	98	93	82	87	58	98	72	50
Demand					140	180	240	210	175	130	320	280	160

kable(Plant_Cst, caption = "Plant Cost")

Table 2.19: Plant Cost
Plant	Variable costs	Capacity
Plant 1	8	1500
Plant 2	13	500

Let’s set the data organized for this model, this could’ve been done before easily when the data was stratified, but it will be cool to learn a bit more of data wrangling

Plants <- grep("Plant*", names(data), value = T)
Demand <- grep("RW*", names(data), value = T)
CD <- grep("CW*", data$Type, value = T)

var_Dc <- as.vector(as.matrix(data[complete.cases(data),"Variable costs"]))
var_Pc <- as.matrix(Plant_Cst$`Variable costs`)

 inc <- t(as.matrix(data[complete.cases(data),Plants])) * var_Dc
 outc <- as.matrix(data[data$Type %in% CD,Demand]) * var_Dc
inc <- t(as.matrix(data[complete.cases(data),Plants]))
outc <- as.matrix(data[data$Type %in% CD,Demand])
S <- as.matrix(Plant_Cst$Capacity)
D <- as.matrix(data[data$Type %in% "Demand",Demand])
Fc <- as.matrix(data[complete.cases(data),"Fixed costs"])

M <- sum(D)*5 #Big Number for flow constraint

rm(x,u,z)

n = length(Plants) #Number of Plants
t = length(CD) #Number of transhipments (CDs)
m = length(Demand) #Number of cities (POS)

#Get back to model
model <- MIPModel() %>%
  
  # Variable of inflow
  add_variable(x[i,j], i = 1:n, j = 1:t, type = "integer", lb = 0) %>%
  
   #Variable of outflow
  add_variable(u[j,k], j = 1:t, k = 1:m, type = "integer", lb = 0) %>%
  
  #Variable of DC to incur a fixed cost
  add_variable(z[j], j = 1:t, type = "binary") %>%

   #minimize total cost
  set_objective(sum_expr(x[i,j] * inc[i,j],  i = 1:n, j = 1:t) +  #Inbound Distance cost
                sum_expr(u[j,k] * outc[j,k], j = 1:t, k = 1:m) +  #Outbound Distance cost
                sum_expr(x[i,j] * var_Dc[j], i = 1:n, j = 1:t) +  #DC variable cost
                sum_expr(x[i,j] * var_Pc[i], i = 1:n, j = 1:t) +  #Plant variable cost
                sum_expr(z[j] * Fc[j], j = 1:t), "min") %>%       #Plant Fixed cost
  
  #Plant Capacity
  add_constraint(sum_expr(x[i,j], j = 1:t) <= S[i], i = 1:n) %>%
  
  #Attend all demand
  add_constraint(sum_expr(u[j,k], j = 1:t) >= D[k], k = 1:m) %>%
  
  #Binary flow constraint
  add_constraint(sum_expr(x[i,j], i = 1:n) <= z[j]*M, j = 1:t) %>%
  
  #Only one DC
  add_constraint(sum_expr(z[j], j = 1:t) <= 1) %>%
  
  #Flow constraint
  add_constraint(sum_expr(x[i,j], i = 1:n) - sum_expr(u[j,k], k = 1:m) == 0,
                 j = 1:t)

#Solve
result <- solve_model(model, with_ROI(solver = "glpk"))


 #Which central warehouse minimizes all cost?
kable(matrix(get_solution(result, x[i,j])$value, nrow = n, ncol = n, dimnames = list(NULL,Plants)))

Plant 1	Plant 2
0	1335
0	500

kable(matrix(get_solution(result, u[j,k])$value, nrow = t, ncol = m, dimnames = list(CD,Demand)),
      caption = paste0("Optimal Value: ", result$objective_value))

Table 2.20: Optimal Value: 349490
	RW 1	RW 2	RW 3	RW 4	RW 5	RW 6	RW 7	RW 8	RW 9
CW 1	0	0	0	0	0	0	0	0	0
CW 2	140	180	240	210	175	130	320	280	160
CW 3	0	0	0	0	0	0	0	0	0
CW 4	0	0	0	0	0	0	0	0	0
CW 5	0	0	0	0	0	0	0	0	0

1 Week2 - Unconstrained & Constrained Optimization

3 Week4 - Algorithms & Approximations

	Plant_Cap	TS_1	TS_2	TO_1	TO_2	Demand
P_1	200	30	50	12	65	150
P_2	300	23	66	25	22	100
P_3	100	35	14	22	23	110
P_4	150	70	12	40	12	200
P_5	220	65	70	41	15	180

	Plant_Cap	TS_1	TS_2	TO_1	TO_2	Demand
P_1	200	30	50	12	65	150
P_2	300	23	66	25	22	100
P_3	100	35	14	22	23	110
P_4	150	70	12	40	12	200
P_5	220	65	70	41	15	180

	Plant_Cap	TS_1	TS_2	TO_1	TO_2	Demand
P_1	200	30	50	12	65	150
P_2	300	23	66	25	22	100
P_3	100	35	14	22	23	110
P_4	150	70	12	40	12	200
P_5	220	65	70	41	15	180