40.014 Engineering Systems Architecture - Homework 2
Name: Ong Jun Sheng (Jarvin)
Student ID: 1002040
Date: 14/04/18
Problem 1 – Aquavista [50 points + 10 bonus points]
Task 1
Duo reservoir will be used more as it has a lower cost of transport per unit volume of water. A good water use strategy will be to use as much water from the Duo reservoir as possible and allocating the minimum amount of water possible to the Uno reservoir, constrained by the water demand, electricity demand and the reservoir capacities.
Task 2
min \(8x_1+3x_2\) st
\(x_1+x_2≥300\) (Water demand) \(9x_1+5x_2≥2000\) (Electricity demand)
\(x_1≤250\) (Uno capacity)
\(x_2≤200\) (Duo capacity) \(x_1, x_2 ≥ 0\)
Task 3
library(lpSolve)
constraint_coefficients <- matrix(c(1, 1,
9, 5,
1, 0,
0, 1,
1, 0,
0, 1),
nrow = 6, byrow = TRUE)
result_p1t3 <- lp(objective.in = c(8, 3),
const.mat = constraint_coefficients,
const.dir = c('>=', '>=', '<=', '<=', '>=', '>='),
const.rhs = c(300, 2000, 250, 200, 0, 0))
print(paste("Optimal value of objective function:", result_p1t3$objval))## [1] "Optimal value of objective function: 1488.88888888889"print(paste("x1 =", result_p1t3$solution[1]))## [1] "x1 = 111.111111111111"print(paste("x2 =", result_p1t3$solution[2]))## [1] "x2 = 200"The solution matches my anticipation as the water used from the Duo reservoir is at its upper bound of 200, while the water used from the Uno reservoir is the minimum required to fulfil the electricity demand constraint.
Task 4
Uno reservoir will be used more as it has a negative cost. A good water use strategy will be to use as much water as possible from the Uno reservoir and allocating the minimum amount of water possible to be used from the Duo reservoir, subject to water demand, electricity demand and the reservoir capacities constraints.
Task 5
constraint_coefficients <- matrix(c(1, 1,
9, 5,
1, 0,
0, 1,
1, 0,
0, 1),
nrow = 6, byrow = TRUE)
result_p1t5 <- lp(objective.in = c(-4, 3),
const.mat = constraint_coefficients,
const.dir = c('>=', '>=', '<=', '<=', '>=', '>='),
const.rhs = c(300, 2000, 250, 200, 0, 0))
print(paste("Optimal value of objective function:", result_p1t5$objval))## [1] "Optimal value of objective function: -850"print(paste("x1 =", result_p1t5$solution[1]))## [1] "x1 = 250"print(paste("x2 =", result_p1t5$solution[2]))## [1] "x2 = 49.9999999999999"Task 6
wpi_p1t3 <- -4*result_p1t3$solution[1] + 3*result_p1t3$solution[2]
cost_p1t5 <- 8*result_p1t5$solution[1] + 3*result_p1t5$solution[2]
print(paste("Water pollution index for Task 3 =", wpi_p1t3))## [1] "Water pollution index for Task 3 = 155.555555555555"print(paste("Daily cost of water supply for Task 5 =", cost_p1t5))## [1] "Daily cost of water supply for Task 5 = 2150"Task 7
print(paste("Range for daily cost of water supply:",
result_p1t3$objval,
"to",
cost_p1t5))## [1] "Range for daily cost of water supply: 1488.88888888889 to 2150"print(paste("Range for water pollution index:",
result_p1t5$objval,
"to",
wpi_p1t3))## [1] "Range for water pollution index: -850 to 155.555555555555"I would suggest varying epsilon from 1489 to 2150.
Task 8
epsilon_values <- seq(1489, 2150, length.out = 5)
constraint_coefficients <- matrix(c(1, 1,
9, 5,
1, 0,
0, 1,
1, 0,
0, 1,
8, 3),
nrow = 7,
byrow = TRUE)
result_p1t8 <- data.frame()
for(epsilon_value in epsilon_values) {
result <- lp(objective.in = c(-4, 3),
const.mat = constraint_coefficients,
const.dir = c(">=", ">=", "<=", "<=", ">=", ">=", "<="),
const.rhs = c(300, 2000, 250, 200, 0, 0, epsilon_value))
# uno[epsilon_value] <- result$solution[1]
# duo[epsilon_value] <- result$solution[2]
# cost[epsilon_value] <- 8*result$solution[1] + 3*result$solution[2]
# wpi[epsilon_value] <- result$objval
result_row <- data.frame(epsilon_value,
result$solution[1],
result$solution[2],
8*result$solution[1] + 3*result$solution[2],
result$objval)
result_p1t8 <- rbind(result_p1t8, result_row)
}
colnames(result_p1t8) <- c("epsilon",
"uno",
"duo",
"cost",
"wpi")
print(result_p1t8)## epsilon uno duo cost wpi
## 1 1489.00 111.1538 199.9231 1489.00 155.1538
## 2 1654.25 150.8500 149.1500 1654.25 -155.9500
## 3 1819.50 183.9000 116.1000 1819.50 -387.3000
## 4 1984.75 216.9500 83.0500 1984.75 -618.6500
## 5 2150.00 250.0000 50.0000 2150.00 -850.0000library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.3.3Decision Space
ggplot(result_p1t8, aes(x = uno, y = duo)) +
labs(x = "Water used from Uno Reservoir",
y = "Water used from Duo Reservoir",
title = "Decision Space (Water used from Uno and Duo Reservoirs)") +
geom_point()
Objective Space
ggplot(result_p1t8, aes(x = cost, y = wpi)) +
labs(x = "Cost",
y = "Water Pollution Index",
title = "Objective Space (Cost vs Water Pollution Index)") +
geom_point()
Problem 2
Task 1
library(goalprog)
coefficients_p2t1 <- matrix(c(4, 2,
4, 2,
1, 0,
0, 1),
nrow=4, byrow=TRUE)
targets_p2t1 <- c(48, 32, 7, 10)
objective_p2t1 <- c(1, 2, 3, 4)
priority_p2t1 <- c(1, 2, 3, 4)
positive_deviations_p2t1 <- c(0, 2, 0, 0)
negative_deviations_p2t1 <- c(1, 1, 5, 5)
achievements_p2t1 <- data.frame(objective = objective_p2t1,
priority = priority_p2t1,
p = positive_deviations_p2t1,
n = negative_deviations_p2t1)
solution_p2t1 <- llgp(coefficients = coefficients_p2t1,
targets = targets_p2t1,
achievements = achievements_p2t1)
print(paste("Decision Variables:",
"x1 =",
solution_p2t1$out$x[1],
", x2 =",
solution_p2t1$out$x[2]))## [1] "Decision Variables: x1 = 7 , x2 = 10"print(paste("Objective Function Value:",
"G1 =",
solution_p2t1$out$f[1,],
", G2 =",
solution_p2t1$out$f[2,],
", G3 =",
solution_p2t1$out$f[3,],
", G4 =",
solution_p2t1$out$f[4,]))## [1] "Objective Function Value: G1 = 48 , G2 = 48 , G3 = 7 , G4 = 10"No, not all the goals are met. There is a negative deviation for G4, which means the demand for Wonder is not met. There is also a positive deviation for G3, which means the production of Marvel exceeds the goal.
Task 2
coefficients_p2t2 <- matrix(c(4, 2,
4, 2,
1, 0,
0, 1),
nrow=4, byrow=TRUE)
targets_p2t2 <- c(48, 32, 7, 10)
objective_p2t2 <- c(1, 2, 3, 4)
priority_p2t2 <- c(1, 1, 2, 3)
positive_deviations_p2t2 <- c(0, 2, 0, 0)
negative_deviations_p2t2 <- c(1, 1, 5, 5)
achievements_p2t2 <- data.frame(objective = objective_p2t2,
priority = priority_p2t2,
p = positive_deviations_p2t2,
n = negative_deviations_p2t2)
solution_p2t2 <- llgp(coefficients = coefficients_p2t2,
targets = targets_p2t2,
achievements = achievements_p2t2)
print(paste("Decision Variables:",
"x1 =",
solution_p2t2$out$x[1],
", x2 =",
solution_p2t2$out$x[2]))## [1] "Decision Variables: x1 = 7 , x2 = 2"print(paste("Objective Function Value:",
"G1 =",
solution_p2t2$out$f[1,],
", G2 =",
solution_p2t2$out$f[2,],
", G3 =",
solution_p2t2$out$f[3,],
", G4 =",
solution_p2t2$out$f[4,]))## [1] "Objective Function Value: G1 = 32 , G2 = 32 , G3 = 7 , G4 = 2"Task 3
Objective Function Values
obj_val_p2t3 <- data.frame(G1 = c(50, 0),
G2 = c(50, 0),
G3 = c(10, 0),
G4 = c(10, 0))
obj_val_p2t3 <- rbind(obj_val_p2t3,
as.vector(solution_p2t1$out$f),
as.vector(solution_p2t2$out$f))
library(fmsb)
radarchart(obj_val_p2t3, title = "Objective Function Values (Task 1 Solution vs Task 2 Solution)")
legend(1, 1, lty=c(1,1), lwd=c(2.5,2.5,2.5),
col=c("black","red"),
c("Task 1 Solution","Task 2 Solution"),
cex=0.8)
Deviation of Objective Function Values
obj_val_p2t3x <- data.frame(G1 = c(50, 0),
G2 = c(50, 0),
G3 = c(10, 0),
G4 = c(10, 0))
obj_val_p2t3x <- rbind(obj_val_p2t3x,
as.vector(abs(solution_p2t1$out$f - solution_p2t1$out$b)),
as.vector(abs(solution_p2t2$out$f - solution_p2t2$out$b)))
radarchart(obj_val_p2t3x, title = "Deviation of Objective Function Values\n(Task 1 Solution vs Task 2 Solution)")
legend(1, 1, lty=c(1,1), lwd=c(2.5,2.5,2.5),
col=c("black","red"),
c("Task 1","Task 2"),
cex=0.8)
Problem 3
Task 1
Visualizing the function
min_p3t1 <- -10; max_p3t1 <- 10
objective_function_p3t1 <- function(x) cos(x) - 3*exp(-(x-0.2)^4)
curve(objective_function_p3t1,
from = min_p3t1,
to = max_p3t1,
n = 1000,
ylab = "f(x)")
Using an initial population of 5:
library(GA)## Loading required package: foreach## Loading required package: iterators## Package 'GA' version 3.0.2## Type 'citation("GA")' for citing this R package in publications.ga_p3t1 <- ga(type = "real-valued",
fitness = function(x) -objective_function_p3t1(x),
min = min_p3t1, max = max_p3t1, popSize = 5, maxiter = 25)## Warning in ga(type = "real-valued", fitness = function(x) -
## objective_function_p3t1(x), : The population size is less than 10.curve(objective_function_p3t1,
from = min_p3t1,
to = max_p3t1,
n = 1000,
ylab = "f(x)")
points(ga_p3t1@solution, -ga_p3t1@fitnessValue, col = 2, pch = 19)
GA does not find the global minimum with an initial population of 5.
Using an initial population of 20:
ga_p3t1x <- ga(type = "real-valued",
fitness = function(x) -objective_function_p3t1(x),
min = min_p3t1, max = max_p3t1, popSize = 20, maxiter = 25)
curve(objective_function_p3t1,
from = min_p3t1,
to = max_p3t1,
n = 1000,
ylab = "f(x)")
points(ga_p3t1x@solution, -ga_p3t1x@fitnessValue, col = 2, pch = 19)
GA produces the global minimum with an initial population of 20
Task 2
Bukin_function_n.6 <- function(x, y) {
100*sqrt(abs(y-0.01*x^2)) + 0.01*abs(x+10)
}
ga_p3t2 <- ga(type = "real-valued",
fitness = function(x) -Bukin_function_n.6(x[1],x[2]),
min = c(-15, -3),
max = c(-5, 3))
x_range_p3t2 <- seq(-15, -5, by = 0.1)
y_range_p3t2 <- seq(-5, 3, by = 0.1)
f_p3t2 <- outer(x_range_p3t2, y_range_p3t2, Bukin_function_n.6)
filled.contour(x_range_p3t2,
y_range_p3t2,
f_p3t2,
color.palette = bl2gr.colors,
plot.axes = {
axis(1); axis(2);
points(ga_p3t2@solution[,1],
ga_p3t2@solution[,2],
pch = 3, cex = 2, col = "white", lwd = 2);
points(-10,
1,
pch = 19, cex = 2, col = "red", lwd = 2)
})
Red point is the global minimum. White “+” is the approximated minimizer.