rm (list = ls())
library(timelineS)
library(lubridate)
library(ggplot2)
library(dplyr)
library(reshape2)
library(rPref)
library(mco)
library(MASS)
plot.timeline <- function(lifetime, events.name, start.Date, plot.Name) {
dataP <- data.frame(Events = events.name, Event_Dates = ymd(start.Date) +
years(lifetime))
timelineS(dataP, main = plot.Name,
labels = events.name, label.direction = “up”, label.position = 3)
}
dist.Events <- function (lifetime, events, start.Date, option.Name) {
#sort the events
events <- sort(events, decreasing = T)
#create the distribution of the events
distribution.events <- sapply(events, seq, from = 0, to = lifetime)

#all events unlisted
all.events <- melt(distribution.events)
colnames(all.events) <- c(“frequency”, “event”)
#sort the events ascending
all.events <- all.events[order(all.events$frequency),]
#get the unique sequence of events
unique.events <- all.events[!duplicated(all.events$frequency),]
unique.events$event[which(unique.events$frequency == 0)] <- “DC”
#plot the timeline
plot.timeline(unique.events$frequency, unique.events$event, start.Date, option.Name)
return(unique.events)
}
combine.lifeTimelines <- function(product.1, product.2, p1.dur, p2.dur) {
#function body
base <- list(Names = c(), frequency = c(), duration = c())
base$frequency <- sort(unique(c(product.1$frequency,
product.2$frequency)),
decreasing = FALSE)
base$Names[1] <- “DC”
base$duration[1] <- 0
base$Names[length(base$frequency)] <- “END”
base$duration[length(base$frequency)] <- 0
for(index in 2:(length(base$frequency) – 1)) {
phase.1 <- product.1$event[which(product.1$frequency ==
base$frequency[index])]
phase.2 <- product.2$event[which(product.2$frequency ==
base$frequency[index])]
if(length(phase.1) != 0 && length(phase.2) != 0) {
base$Names[index] <- paste(phase.1, phase.2)
base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)],
p2.dur[which(names(p2.dur) == phase.2)])
} else if(length(phase.1) != 0 && length(phase.2) == 0) {
base$Names[index] <- phase.1
base$duration[index] <- p1.dur[which(names(p1.dur) == phase.1)]
} else if (length(phase.1) == 0 && length(phase.2) != 0) {
base$Names[index] <- phase.2
base$duration[index] <- p2.dur[which(names(p2.dur) == phase.2)]
} else {
base$duration[index] <- 0
}
}
return(base)
}
par(mfrow = c(2, 1))
start.Date = “2020-01-01″
lifetime <- 100
events.Op1.Bridge = c(SDO = 15, M = 5, DR = 30, END = lifetime)
duration.ev.Bridge <- c(SDO = 7, M = 2, DR = 21)
events.Op1.Tunnel <- c(BR = 15, M.r = 7, AACL = 15, END = lifetime)
duration.ev.Tunnel <- c(BR = 21, M.r = 2, AACL = 14)
design.Options.Bridge <- list(desing.Op1 <- “1.Cast in place deck and steel girders”)
design.Options.Tunnel <- list(desing.Op1.Tunnel <- “Tunnel”)
maintenance.Bridge <- dist.Events(lifetime, events.Op1.Bridge, start.Date, design.Options.Bridge)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
integrated.interv <- combine.lifeTimelines(maintenance.Bridge, maintenance.Tunnel,
duration.ev.Bridge, duration.ev.Tunnel)
sum(integrated.interv$duration)
events.Op1.Bridge = c(SDO = 15, M = 5, DR = 30, END = lifetime)
duration.ev.Bridge <- c(SDO = 7, M = 2, DR = 21)
events.Op1.Tunnel <- c(BR = 17, M.r = 6, AACL = 15, END = lifetime)
duration.ev.Tunnel <- c(BR = 21, M.r = 2, AACL = 14)
design.Options.Bridge <- list(desing.Op1.Bridge <- “1.Cast in place deck and steel girders”)
design.Options.Tunnel <- list(desing.Op1.Tunnel <- “Tunnel”)
maintenance.Bridge <- dist.Events(lifetime, events.Op1.Bridge, start.Date,
design.Options.Bridge)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date,
design.Options.Tunnel)
integrated.interv <- combine.lifeTimelines(maintenance.Bridge, maintenance.Tunnel,
duration.ev.Bridge, duration.ev.Tunnel)
sum(integrated.interv$duration)
design.explore <- function(events1, events2) {
results <- c()
for(i in 1: dim(events1)[1]) {
ev1 <- unlist(events1[i, ])
dist.1 <- dist.Events(lifetime, ev1, start.Date, desing.Op1)
dur.ev1 <- ev1/2
for (j in 1: dim(events1)[1]) {
ev2 <- unlist(events2[j, ])
dist.2 <- dist.Events(lifetime, ev2, start.Date, design.Options.Tunnel)
dur.ev2 <- ev2/2
combined.lifetime <- combine.lifeTimelines(dist.1, dist.2,
dur.ev1, dur.ev2)
min.dist.int <- min(abs(combined.lifetime$frequency[1:(length(combined.lifetime$frequency) – 1)] –
combined.lifetime$frequency[2:length(combined.lifetime$frequency)]))
results <- rbind(results, c(ev1, ev2, dur = sum(combined.lifetime$duration), dist.inter = min.dist.int))
}
}
return(as.data.frame(results))
}
n.grid <- 5
events.grid.bridge <- expand.grid(SDO = sample(seq(10,20, by = 1), n.grid),
M = sample(seq(3,8,1), n.grid),
DR = sample(seq(20, 35, 1), n.grid))
events.grid.Tunnel <- expand.grid(BR = sample(seq(10, 20, by = 1), n.grid),
M.r = sample(seq(3,8,1), n.grid),
R.r = sample(seq(10, 20, 1), n.grid))
response.space <- design.explore(events.grid.bridge, events.grid.Tunnel)
p <- low(dur)* high(dist.inter)
sky <- psel(response.space, p)
pareto2 <- psel(response.space, p, top = nrow(response.space))
ggplot(response.space, aes(x = dur, y = dist.inter)) +
geom_point(shape = 21) +
geom_point(data = pareto2, size = 3, aes(color = factor(pareto2$.level)))
show_front <- function(pref) {
plot(response.space$dur, response.space$dist.inter)
sky <- psel(response.space, pref)
plot_front(response.space, pref, col = rgb(0, 0, 1))
points(sky$dur, sky$dist.inter, lwd = 3)
}
show_front(p)
p <- high(dur) * low(dist.inter)
show_front(p)
LCI.materials <- read.csv(“LCImaterials.csv”)
LCA.bridge <- function(length, width, height, thickness, girder.Option,
deck.Option, materials, dist.event) {
prefab.girder.Section <- 0.78
steel.girders.unitWeight <- 317
asphalt.Q <- length * width * thickness
#based on the use of each material, we will split the data frame
materials.split <- split(materials, materials$scope)

#new line to get the interventions
s <- summary(as.factor(dist.event$event))
s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))

# calculate the volume of the deck based on different materials strategies
if(deck.Option == “RC”) {
deck.volume <- length * width * height
interventions.deck <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
} else if (deck.Option == “PC”) {
deck.volume <- 10.2*length
interventions.deck <- s2$DC + s2$BR + 0.25 * s2$AACL + 0.35 * (if (!is.null(s2$PF))
{s2$PF} else {0})

}
#girder options
if (girder.Option == “steel”) {
#get the numbers of girders
n <-round(width / 3, 0)
interventions.girders <- s2$DC + if (!is.null(s2$DR)) {s2$DR} else{ 0 } + 0.55 * (if(!is.null(s2$PR)) {s2$PR} else {0})

#get the volume of the concrete for the prefab girders
girders.V <- n * steel.girders.unitWeight * length
} else if (girder.Option == “C”) {
n <- 0
girders.V <- 7.85*0.05*0.65*(11.3-10)*6*666
interventions.girders <- s2$DC + if (!is.null(s2$BR)) {s2$BR} else{ 0 } + 0.55 * (if(!is.null(s2$AACL)) {s2$AACL} else {0})

}

asph.interventions = length(dist.event$event)-1
asphalt <- mutate(materials.split$asphalt, bridge.Q = asphalt.Q, interventions =
asph.interventions)
deck <- mutate(materials.split[[deck.Option]], bridge.Q = deck.volume, interventions =
interventions.deck)

girders <- mutate(materials.split[[girder.Option]], bridge.Q = girders.V, interventions =
interventions.girders)
LCA.matrix <- rbind(deck, girders, asphalt)
LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * bridge.Q / 1000,
materials.LC = TotalMaterials.Q * interventions,
Energy.LC = materials.LC * energy,
CO2.LC = materials.LC * CO2 * 1000,
NOx.LC = materials.LC * NOx * 1000,
SO2.LC = materials.LC * SO2 * 1000)
LCA.results <- list(Energy = sum(LCA.matrix$Energy.LC),
CO2 = sum(LCA.matrix$CO2.LC),
NOx = sum(LCA.matrix$NOx.LC),
SO2 = sum(LCA.matrix$SO2.LC))

return(LCA.results)
}
b.length <- 500 # units: m
b.width <- 15 #units: m
bd.depth <- 0.35 #units: m
asphalt.tk <- 0.12 #units: m
t.depth <- 0.35 #units: m
t.lenght <- 1000 #m 100 meters on each side of the bridge
girder.Options <- c(“steel”, “C”)
deck.options <- c(“RC”, “PC”)
Option.bridge <- LCA.bridge(b.length, b.width, bd.depth, asphalt.tk,
girder.Options[1], deck.options[1], LCI.materials, maintenance.Bridge)
Option.Tunnel <- LCA.bridge(t.lenght, b.width, t.depth, asphalt.tk,
girder.Options[2], deck.options[2], LCI.materials, maintenance.Tunnel)
integrated.Design <- as.data.frame(list(Energy = Option.bridge$Energy + Option.Tunnel$Energy,
CO2 = Option.bridge$CO2 + Option.Tunnel$CO2,
NOx = Option.bridge$NOx + Option.Tunnel$NOx,
SO2 = Option.bridge$SO2 + Option.Tunnel$SO2))
Energy.costs <- 0.128
CO2.unitcost <- 26 # per metric tone
NOx.unitCost <- 42 # per metric tone
SO2.unitCosts <- 85 # per metric tone
integrated.Design <- mutate(integrated.Design,
Costs = (Energy * Energy.costs +
CO2*CO2.unitcost +
NOx * NOx.unitCost +
SO2 * SO2.unitCosts)/10^9)
fitness <- function(x) {
#define the output dimension
z <- numeric(7)
#function body
x <- round(x, 0)
y <- expand.grid(SDO = x[1], M = x[2], DR = x[3], BR = x[4], M.r = x[5], AACL = x[6])

dur.ev1 <- unlist(y[1:3] / 2)
dur.ev2 <- unlist(y[4:6] / 2)

dist.1 <- apply(y[1:3], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.Bridge)
dist.2 <- apply(y[4:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.Tunnel)
results <- combine.lifeTimelines(dist.1[[1]], dist.2[[1]], dur.ev1, dur.ev2)
z[1] <- sum(results[[“duration”]])
z[2] <- -min(abs(results$frequency[1:(length(results$frequency) – 1)] –
results$frequency[2:length(results$frequency)]))
b.width <- x[7]
b.length <- 25 # units: m
bd.depth <- 0.35 #units m
asphalt.tk <- 0.12 #units m
t.depth <- 0.35 #units: m
t.lenght <- 1000 #m 100 meters on each side of the bridge
materials <- LCI.materials
girder.Options <- c(“steel”, “C”)
deck.options <- c(“RC”, “PC”)
Energy.costs <- 0.128
CO2.unitcost <- 26
NOx.unitCost <- 42
SO2.unitCosts <- 85
product.1.bridge <- LCA.bridge(b.length, b.width, bd.depth, asphalt.tk,
girder.Options[1], deck.options[1], LCI.materials, dist.1[[1]])

product.2.Tunnel <- LCA.bridge(t.lenght, b.width, t.depth, asphalt.tk,
girder.Options[2], deck.options[2], LCI.materials, dist.2[[1]])

integrated.system <- as.data.frame(list(Energy = product.1.bridge$Energy +
product.2.Tunnel$Energy,
CO2 = product.1.bridge$CO2 + product.2.Tunnel$CO2,
NOx = product.1.bridge$NOx + product.2.Tunnel$NOx,
SO2 = product.1.bridge$SO2 + product.2.Tunnel$SO2))
integrated.system <- mutate(integrated.system,
Costs = (Energy * Energy.costs +
CO2*CO2.unitcost +
NOx * NOx.unitCost +
SO2 * SO2.unitCosts)/10^9)
z[3] <- integrated.system$Energy
z[4] <- integrated.system$CO2
z[5] <- integrated.system$NOx
z[6] <- integrated.system$SO2
z[7] <- integrated.system$Costs
return(z)
}
r2 <- nsga2(fitness, idim = 7, odim = 7,
generations=5, popsize=4,
lower.bounds=c(10, 3, 20, 10, 10, 3, 25),
upper.bounds= c(20, 8, 35, 20, 20, 8, 45))
r2 <- nsga2(fitness, idim = 7, odim = 7,
generations=10, popsize=100,
lower.bounds=c(10, 3, 20, 12, 4, 27, 25),
upper.bounds= c(20, 8, 35, 25, 10, 37, 45))
r2Results <- as.data.frame(r2$value)
outNames <- c(“duration”, “interv.dist”, “energy”, “co2″, “nox”, “so2″, “cost” )
colnames(r2Results) <- outNames
pareto3 <- as.data.frame(paretoFront(r2))
colnames(pareto3) <- outNames
ggplot(r2Results, aes(x = duration, y = cost)) + geom_point(shape = 21) +
geom_point(data = pareto3, size = 3, color=”red”) +
geom_line(data = pareto3, color=”blue”)
input.params <- round(r2$par, 0)
all.results <- cbind(input.params, r2Results, r2$pareto.optimal)
colnames(all.results) <- c(“SDO”, “M”, “DR”,
“BR”, “M.r”, “AACL”,
“b.width”, “duration”, “interv.dist”,
“energy”, “co2″, “nox”, “so2″, “cost”,
“pareto”)
par(mfrow = c(1,1))
parcoord(all.results[, 1:14], var.label = T,
col = ifelse(all.results$pareto == TRUE, “indianred”, “skyblue2″),
lty = ifelse(all.results$pareto == TRUE, 1, 3),
lwd = ifelse(all.results$pareto == TRUE, 3, 1))