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, product.3, product.4,
p1.dur, p2.dur, p3.dur, p4.dur) {
base <- list(Names = c(), frequency = c(), duration = c())
base$frequency <- sort(unique(c(product.1$frequency,
product.2$frequency,
product.3$frequency,
product.4$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])]
phase.3 <- product.3$event[which(product.3$frequency ==
base$frequency[index])]
phase.4 <- product.4$event[which(product.4$frequency ==
base$frequency[index])]

list <- c(phase.1, phase.2, phase.3, phase.4)
dist <- c(p1.dur, p2.dur, p3.dur, p4.dur)

l <- c()
d <- c(0)

for (i in 1:length(list)){
if (length(list[i] != 0)){
l <- append(l, list[i])
x <- dist[which(names(dist) == list[i])]
d <- append(d, x)
}
else{
d <- append(d, 0)
}
}

base$Names[index] <- paste(l, collapse = “&”)
base$duration[index] <- max(d)
}
return(base)
}

par(mfrow = c(2, 1))
start.Date = “2020-01-01″
lifetime <- 70

events.Op1.Metro = c(RR.m = 28, M.m = 5, FER.m = 15, CS.m = 50, END=lifetime)   ## Re-read!!!!
duration.ev.Metro <- c(RR.m = 7, M.m = 2, FER.m = 21)

events.Op1.Bridge = c(SDO.b = 18, M.b = 6, DR.b = 36, END = lifetime)
duration.ev.Bridge <- c(SDO.b = 4, M.b = 2, DR.b = 32)

events.Op1.Building <- c(SDR.bu = 15, M.bu = 10, DR.bu = 50, END = lifetime)
duration.ev.Building <- c(SDO.bu = 5, M.bu = 6, DR.bu = 11)

events.Op1.Water <- c(PI.w = 2, M.w = 10, PR.w = 55, END = lifetime)
duration.ev.Water <- c(PI.w = 8, M.w = 12, PR.w = 31)

design.Options.Metro <- list(desing.Op1.Metro <- “Metro”)
design.Options.Bridge <- list(desing.Op1.Bridge <- “Bridge”)
design.Options.Building <- list(desing.Op1.Building <- “Building”)
design.Options.Water <- list(desing.Op1.Water <- “Water”)

maintenance.Metro <- dist.Events(lifetime, events.Op1.Metro, start.Date,
design.Options.Metro)
maintenance.Bridge <- dist.Events(lifetime, events.Op1.Bridge, start.Date,
design.Options.Bridge)
maintenance.Building <- dist.Events(lifetime, events.Op1.Building, start.Date,
design.Options.Building)
maintenance.Water <- dist.Events(lifetime, events.Op1.Water, start.Date,
design.Options.Water)

integrated.interv <- combine.lifeTimelines(maintenance.Metro, maintenance.Bridge,
maintenance.Building, maintenance.Water,
duration.ev.Metro, duration.ev.Bridge,
duration.ev.Building, duration.ev.Water)
sum(integrated.interv$duration)

design.explore <- function(events1, events2, events3, events4) {
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.Road)
dur.ev2 <- ev2/2
for (k in 1: dim(events1)[1]) {
ev3 <- unlist(events3[k, ])
dist.3 <- dist.Events(lifetime, ev3, start.Date, design.Options.Road)
dur.ev3 <- ev3/2
for (l in 1: dim(events1)[1]) {
ev4 <- unlist(events4[k, ])
dist.4 <- dist.Events(lifetime, ev4, start.Date, design.Options.Road)
dur.ev4 <- ev4/2
combined.lifetime <- combine.lifeTimelines(dist.1, dist.2, dist.3, dist.4,
dur.ev1, dur.ev2, dur.ev3, dur.ev4)
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, ev3, ev4, dur = sum(combined.lifetime$duration), dist.inter = min.dist.int))
}
}
}
}
return(as.data.frame(results))
}

n.grid <- 2

events.Op1.Metro = c(RR.m = 28, M.m = 5, FER.m = 15, CS.m = 50, END=lifetime)   ## Re-read!!!!
duration.ev.Metro <- c(RR.m = 7, M.m = 2, FER.m = 21)

events.Op1.Bridge = c(SDO.b = 18, M.b = 6, DR.b = 36, END = lifetime)
duration.ev.Bridge <- c(SDO.b = 4, M.b = 2, DR.b = 32)

events.Op1.Building <- c(SDR.bu = 15, M.bu = 10, DR.bu = 50, END = lifetime)
duration.ev.Building <- c(SDO.bu = 5, M.bu = 6, DR.bu = 11)

events.Op1.Water <- c(PI.w = 2, M.w = 10, PR.w = 55, END = lifetime)
duration.ev.Water <- c(PI.w = 8, M.w = 12, PR.w = 31)

events.grid.metro <- expand.grid(SDO.m = sample(seq(10,20, by = 1), n.grid),
M.m = sample(seq(3,8,1), n.grid),
DR.m = sample(seq(20, 35, 1), n.grid))

events.grid.bridge <- expand.grid(SDO.br = sample(seq(12, 25, by = 1), n.grid),
M.br = sample(seq(5,10,1), n.grid),
R.br = sample(seq(27, 37, 1), n.grid))

events.grid.building <- expand.grid(SDO.bu = sample(seq(10,20, by = 1), n.grid),
M.bu = sample(seq(3,8,1), n.grid),
DR.bu = sample(seq(20, 35, 1), n.grid))

events.grid.water <- expand.grid(SDO.w = sample(seq(12, 25, by = 1), n.grid),
M.w = sample(seq(5,10,1), n.grid),
R.w = sample(seq(27, 37, 1), n.grid))

response.space <- design.explore(events.grid.metro, events.grid.bridge, events.grid.building, events.grid.water)

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)

## Metro

LCI.materials <- read.csv(“LCI_Materials-5.csv”, sep=”;”, dec=”,”)

LCA.station <- function(length, superstructure.Option, materials, dist.event) {
substructure.section <- 0.506
conv.superstructure.section <- 0.46575
FRP_conv.superstructure.section <- 0.46575
prefab.superstructure.section <- 0.3795
FRP_prefab.superstructure.section <- 0.3795
rail.unitWeight <- 54.77

materials.split <- split(materials, materials$scope)

s <- summary(as.factor(dist.event$event))
s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors = F, row.names = 1:length(s))

if(superstructure.Option == “CRC”) {
superstructure.volume <- length * conv.superstructure.section * 2 # multiplication by “2” is for each track line direction
}
else if (superstructure.Option == “PRC”) {
superstructure.volume <- length * prefab.superstructure.section * 2
}
else if (superstructure.Option == “CFRC”) {
superstructure.volume <- length * FRP_conv.superstructure.section * 2
}
else if (superstructure.Option == “PFRC”) {
superstructure.volume <- length * FRP_prefab.superstructure.section * 2

}

interventions.sup <- s2$DC + s2$M + (if (!is.null(s2$PR)) {s2$PR} else {0}) +
(if (!is.null(s2$FPR)) {s2$FPR} else {0})

substructure.volume <- length * substructure.section * 2
interventions.sub <- s2$DC + s2$M + (if (!is.null(s2$CS)) {s2$CS} else {0})

rail.weight <- length * rail.unitWeight * 2 * 2 #multiplication by 4 is because there are 2 rail on each direction
interventions.rail <- s2$DC + s2$RR + s2$M + (if (!is.null(s2$FER)) {s2$FER} else {0})

rail <- mutate(materials.split$Rail, metro.Q = rail.weight, interventions = interventions.rail)
superstructure <- mutate(materials.split[[superstructure.Option]], metro.Q = superstructure.volume, interventions = interventions.sup)
substructure <- mutate(materials.split[[“CRC”]], metro.Q = substructure.volume, interventions = interventions.sub)

LCA.matrix <- rbind(substructure, superstructure, rail)

LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * metro.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,
TotalCost.LC = materials.LC * UnitPrice * 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),
Cost = sum(LCA.matrix$TotalCost.LC))
return(LCA.results)
}

metro.length <- 80

# chosen.Superstructure.Option <- readline(prompt = “Enter metro station ooption index (1-4): “)
superstructure.Options <- c(“CRC”, “PRC”, “CFRC”, “PFRC”)

Option.metro <- LCA.station(metro.length, superstructure.Options[2], LCI.materials, maintenance.Metro)

## Bridge

LCA.bridge <- function(length, width, height, thickness, girder.Option, deck.Option, materials, dist.event) {
steel.girders.unitWeight <- 317 #the weight for HEM800 steel profile
asphalt.Q <- length * width * thickness

materials.split <- split(materials, materials$scope)

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 == “steel”) {
deck.volume <-    0.1 * 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 == “rsteel”) {
deck.volume <- 0.1 * length * width * height
interventions.deck <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
}
#girder options
if (girder.Option == “steel”) {
n <- round(width / 3, 0)
girders.V <- n * steel.girders.unitWeight * length
interventions.girders <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})

} else if (girder.Option == “rsteel”) {
n <- round(width / 3, 0)
girders.V <- n * steel.girders.unitWeight * length
interventions.girders <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
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)

# bridge.Q in m^3 or t (steel)
LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * bridge.Q,
materials.LC = TotalMaterials.Q * interventions,
Energy.LC = materials.LC * energy,
CO2.LC = materials.LC * CO2,
Cost.LC = materials.LC * cost)
LCA.results <- list(Energy = sum(LCA.matrix$Energy.LC),
CO2 = sum(LCA.matrix$CO2.LC),
Cost = sum(LCA.matrix$Cost.LC))
return(LCA.results)
}

bridge.length <- 30 # units: m
bridge.width <- 10 #units m
bridge.depth <- 0.25 #units m
bridge_asphalt.tk <- 0.12 #units m

girder.Options <- c(“steel”, “rsteel”)
deck.options <- c(“RC”, “steel”, “rsteel”)

#chosen.Bridge_deck.Option <- readline(prompt = “Enter bridge deck option index (1-3): “)
#chosen.Bridge_girder.Option <- readline(prompt = “Enter bridge girder option index (1-2): “)

Option.bridge <- LCA.bridge(bridge.length, bridge.width, bridge.depth, bridge_asphalt.tk,
girder.Options[2], deck.options[2], LCI.materials, maintenance.Bridge)

## Building

LCA.pile <- function(length, width, depth, th, pile.Option, materials, dist.event)
{
concpile.crs <- 0.018
concpile.cap.section <- 16
timber.crs <- 0.19
steel.crs <- 49.50
cap.Q <- concpile.cap.section * th

materials.split <- split(materials, materials$scope)

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 pile based on different materials strategies
if(pile.Option == “RC”) {
x <- round(width / 5)
y <- round(length / 5)
n <- x * y
pile.volume <- n * concpile.crs * depth
interventions.deck <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
} else if (pile.Option == “steel”) {
x <- round(width / 5)
y <- round(length / 5)
n <- x * y
pile.volume <- n * timber.crs * depth
interventions.deck <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
} else if (pile.Option == “wood”) {
x <- round(width / 5)
y <- round(length / 5)
n <- x * y
pile.volume <- n * steel.crs * depth
interventions.deck <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
}
cap.interventions = length(dist.event$event)-1
cap <- mutate(materials.split$RC, foun.Q = cap.Q, interventions = cap.interventions)
pile <- mutate(materials.split[[pile.Option]], foun.Q = pile.volume, interventions =
interventions.deck)
LCA.matrix <- rbind(pile, cap)
LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * foun.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)
}

p.length <- 16 # units: m
p.width <- 15 #units m
p.depth <- 50 #units m
p.th <- 0.12 #units m

#chosen.pile.Option <- readline(prompt = “Enter pile option index (1-3): “)
pile.Option <- c(“RC”, “steel”, “wood”)

Option.building <- LCA.pile(p.length, p.width, p.depth, p.th, pile.Option[2], LCI.materials, maintenance.Building)

## Water Distribution System
LCA.pipe <- function(length, time.length, externalDiameter, internalDiameter, pipe.Option, soil.option, sources, dist.event) {

pipe.frequency <- 398.97 # value of corrosion mass loss rate

# pipe option
# calculate corrosion rate

if(pipe.Option == “Cast Iron”) {
pipe.density <- 7850 #unit kg/m3
pipe.initial.mass <- 47.34*length #unit kg (weight/length = 31.24 kg/m)
pipe.modulus <- 206000 # unit MPa
pipe.poisson <- 0.3
interventions.pipe <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})

} else if (pipe.Option == “Ductile Iron”) {
pipe.density <- 7850 #unit kg/m3
pipe.initial.mass <- 47.34*length #unit kg (weight/length = 31.24 kg/m)
pipe.modulus <- 206000 # unit MPa
pipe.poisson <- 0.3
interventions.pipe <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})

} else if (pipe.Option == “PVC”) {
pipe.density <- 1450 #unit kg/m3
pipe.initial.mass <- 25.44*length #unit kg (weight/length = 16.21 kg/m)    pipe.modulus <- 3275 # unit MPa
pipe.poisson <- 0.3
interventions.pipe <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})

} else if (pipe.Option == “Concrete”) {
pipe.external.diameter <-  750 #unit mm
pipe.density <- 2500 #unit kg/m3
pipe.initial.mass <- 227.5*length #unit kg (weight/length = 16.21 kg/m)
pipe.modulus <- 43800 # unit MPa
pipe.poisson <- 0.2
interventions.pipe <- s2$DC + s2$DR + 0.25 * s2$SDO + 0.35 * (if (!is.null(s2$PR)) {s2$PR}
else {0})
}

pipe.mass.loss.rate <- (1/(2*time.length))*(externalDiameter-sqrt(((4*pipe.frequency*pipe.initial.mass)/(pi*pipe.density))-(internalDiameter**2)))
#  pipe.mass.loss.rate <- abs(pipe.mass.loss.rate)

log_OTC <- filter(sources, type == soil.option)
log_OTC <- select(log_OTC, OTC)

log_SMZ <- filter(sources, type == soil.option)
log_SMZ <- select(log_SMZ, SMZ)

log_NOR <- filter(sources, type == soil.option)
log_NOR <- select(log_NOR, NOR)

LCA.matrix <- cbind(pipe.mass.loss.rate, log_OTC, log_SMZ, log_NOR)

LCA_result <- list(pipe.mass.loss.rate,
log_OTC,
log_SMZ,
log_NOR)

return(LCA_result)
}

pipe.length <- 10 # units: m
pipe.time <- 20 #units years
pipe.external.diameter <-  0.400 #units: m
pipe.internal.diameter <- 0.360 #unit: m

# chosen.pipe.Option <- readline(prompt = “Enter pipe option index (1-4): “)
pipe.Options <- c(“Cast Iron”, “Ductile Iron”, “PVC”, “Concrete”)
#pipe.density <- c(7850, 7850, 1450, 2500)
# chosen.soil.Option <- readline(prompt = “Enter soil option index (1-6): “)
soil.option <- c(“Fluvo-aquic soil”, “Red soil”, “Yellow soil”, “Bog soil”, “Brown soil”, “Latosolic”)

Option.water <- LCA.pipe(p.length, p.time, p.external.diameter, p.internal.diameter, pipe.Options[2],
soil.option[2], LCI.materials, maintenance.Water)

### combined

integrated.Design <- as.data.frame(list(Energy = Option.metro$Energy + Option.bridge$Energy + Option.building$Energy + Option.water$Energy,
CO2 = Option.metro$CO2 + Option.bridge$CO2 + Option.building$CO2 + Option.water$CO2,
NOx = Option.metro$NOx + Option.bridge$NOx + Option.building$NOx + Option.water$NOx,
SO2 = Option.metro$SO2 + Option.bridge$SO2 + Option.building$SO2 + Option.water$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)

###

##  Multi-Objective Optimization

fitness <- function(x) {

z <- numeric(7)

y <- expand.grid(SDO.m = x[1], M.m = x[2], DR.m = x[3], SDO.br = x[4], M.br = x[5], R.br = x[6],
SDO.bu = x[7], M.bu = x[8], DR.bu = x[9], SDO.w = x[10], M.w = x[11], R.w = x[12])
dur.ev1 <- unlist(y[1:3] / 2)
dur.ev2 <- unlist(y[4:6] / 2)
dur.ev3 <- unlist(y[7:9] / 2)
dur.ev4 <- unlist(y[10:12] / 2)

dist.1 <- apply(y[1:3], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
dist.2 <- apply(y[4:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
dist.3 <- apply(y[7:9], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
dist.4 <- apply(y[10:12], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)

results <- combine.lifeTimelines(dist.1[[1]], dist.2[[1]], dist.3[[1]], dist.4[[1]],
dur.ev1, dur.ev2, dur.ev3, dur.ev4)

z[1] <- sum(results[[“duration”]])

z[2] <- -min(abs(results$frequency[1:(length(results$frequency) – 1)] –
results$frequency[2:length(results$frequency)]))

# Metro features
metro.length <- 80

chosen.Superstructure.Option <- readline(prompt = “Enter metro station ooption index (1-4): “)
superstructure.Options <- c(“CRC”, “PRC”, “CFRC”, “PFRC”)

product.1.metro <- LCA.station(metro.length, superstructure.Options[chosen.Superstructure.Option],
LCI.materials, dist.1[[1]])

# Bridge features
bridge.length <- 30 # units: m
bridge.width <- 10 #units m
bridge.depth <- 0.25 #units m
bridge_asphalt.tk <- 0.12 #units m

girder.Options <- c(“steel”, “rsteel”)
deck.options <- c(“RC”, “steel”, “rsteel”)

chosen.Bridge_deck.Option <- readline(prompt = “Enter bridge deck option index (1-3): “)
chosen.Bridge_girder.Option <- readline(prompt = “Enter bridge girder option index (1-2): “)

product.2.bridge <- LCA.bridge(bridge.length, bridge.width, bridge.depth, bridge_asphalt.tk,
girder.Options[chosen.Bridge_girder.Option],
deck.options[chosen.Bridge_deck.Option], LCI.materials, dist.2[[1]])

# Building features

p.length <- 16 # units: m
p.width <- 15 #units m
p.depth <- 50 #units m
p.th <- 0.12 #units m

chosen.pile.Option <- readline(prompt = “Enter pile option index (1-3): “)
pile.Option <- c(“RC”, “steel”, “wood”)

product.3.building <- LCA.pile(p.length, p.width, p.depth, p.th, pile.Option[chosen.pile.Option], LCI.materials, dist.3[[1]])

# Water_Distribution features
pipe.length <- 10 # units: m
pipe.time <- 20 #units years
pipe.external.diameter <-  0.400 #units: m
pipe.internal.diameter <- 0.360 #unit: m

chosen.pipe.Option <- readline(prompt = “Enter pipe option index (1-4): “)
pipe.Options <- c(“Cast Iron”, “Ductile Iron”, “PVC”, “Concrete”)
#pipe.density <- c(7850, 7850, 1450, 2500)
chosen.soil.Option <- readline(prompt = “Enter soil option index (1-6): “)
soil.option <- c(“Fluvo-aquic soil”, “Red soil”, “Yellow soil”, “Bog soil”, “Brown soil”, “Latosolic”)

product.4.water <- LCA.pipe(p.length, p.time, p.external.diameter, p.internal.diameter, pipe.Options[chosen.pipe.Option],
soil.option[chosen.soil.Option], LCI.materials, dist.4[[1]])

# Unit Prices
Energy.costs <- 0.128
CO2.unitcost <- 26
NOx.unitCost <- 42
SO2.unitCosts <- 85

## Combination

integrated.system <- as.data.frame(list(Energy = Option.metro$Energy + Option.bridge$Energy + Option.building$Energy + Option.water$Energy,
CO2 = Option.metro$CO2 + Option.bridge$CO2 + Option.building$CO2 + Option.water$CO2,
NOx = Option.metro$NOx + Option.bridge$NOx + Option.building$NOx + Option.water$NOx,
SO2 = Option.metro$SO2 + Option.bridge$SO2 + Option.building$SO2 + Option.water$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, 12, 4, 27, 25),
upper.bounds= c(20, 8, 35, 25, 10, 37, 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”,
“SDO.r”, “M.r”, “R.r”,
“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))