#=====================-------------- GROUP 5 --------------===================== # Selda Demirel - 458462 # Vincent Paul - 365951 # Mudder Marcel Nasser - 370032 # Jamila Loutfi - 355035 # Padmanabha Saikia - 456008 # #-----------------------------------Installing packages # install.packages("timelineS") # install.packages("dplyr") # install.packages("MCDA") # install.packages("ggplot2") # install.packages("scales") # install.packages("reshape2") # install.packages("tibble") # install.packages("lubridate") rm (list = ls()) cat("\014") library(timelineS) library(lubridate) library(ggplot2) library(dplyr) library(reshape2) library(rPref) library(mco) library(MASS) library(MCDA) library(scales) library(tibble) library("rstudioapi") #===================================================== 1. Integrated Maintenance Planning ===================================================== 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) { events <- sort(events, decreasing = T) distribution.events <- sapply(events, seq, from = 0, to = lifetime) all.events <- melt(distribution.events) colnames(all.events) <- c("frequency", "event") all.events <- all.events[order(all.events$frequency),] unique.events <- all.events[!duplicated(all.events$frequency),] unique.events$event[which(unique.events$frequency == 0)] <- "Construction" return(unique.events) } combine.lifeTimelines <- function(product.1, product.2,product.3,product.4,product.5, p1.dur, p2.dur, p3.dur, p4.dur, p5.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, product.5$frequency)), decreasing = FALSE) base$Names[1] <- "Construction" 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])] phase.5 <- product.5$event[which(product.5$frequency == base$frequency[index])] if(length(phase.1) != 0 && length(phase.2) != 0 && length(phase.3) != 0 && length(phase.4) != 0&& length(phase.5) != 0) { base$Names[index] <- paste(phase.1, phase.2, phase.3, phase.4, phase.5) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p2.dur[which(names(p2.dur) == phase.2)], p3.dur[which(names(p3.dur) == phase.3)], p4.dur[which(names(p4.dur) == phase.4)], p5.dur[which(names(p5.dur) == phase.5)]) } else if(length(phase.1) != 0 && length(phase.2) == 0 && length(phase.3) == 0 && length(phase.4) == 0&& length(phase.5) == 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 && length(phase.3) == 0 && length(phase.4) == 0&& length(phase.5) == 0) { base$Names[index] <- phase.2 base$duration[index] <- p2.dur[which(names(p2.dur) == phase.2)] } else if (length(phase.1) == 0 && length(phase.2) == 0 && length(phase.3) != 0 && length(phase.4) == 0&& length(phase.5) == 0) { base$Names[index] <- phase.3 base$duration[index] <- p3.dur[which(names(p3.dur) == phase.3)] } else if (length(phase.1) == 0 && length(phase.2) == 0 && length(phase.3) == 0 && length(phase.4) != 0&& length(phase.5) == 0) { base$Names[index] <- phase.4 base$duration[index] <- p4.dur[which(names(p4.dur) == phase.4)] } else if (length(phase.1) == 0 && length(phase.2) == 0 && length(phase.3) == 0 && length(phase.4) == 0&& length(phase.5) != 0) { base$Names[index] <- phase.5 base$duration[index] <- p5.dur[which(names(p5.dur) == phase.5)] } else { base$duration[index] <- 0 } } return(base) } par(mfrow = c(2, 1)) start.Date = "2021-01-01" lifetime <- 100 events.Op1.Bridge = c(PD_M = 5, PD_R = 35, END = lifetime) duration.ev.Bridge <- c(PD_M = 10, PD_R = 60) events.Op1.Pylon <- c(P_M = 25, P_R = 50, C_A = 10, C_R = 20, END = lifetime) duration.ev.Pylon <- c(P_M = 8, P_R = 40, C_A = 2, C_R = 10) events.Op1.Tunnel <- c(SC_M = 25, SC_I = 5,END = lifetime) duration.ev.Tunnel <- c(SC_M = 20, SC_I = 5) events.Op1.Pavement <- c(AS_R = 20, AIB_R = 50, AS_M = 8, AIB_M = 25, END = lifetime) duration.ev.Pavement <- c(AS_R = 7, AIB_R = 14, AS_M = 2, AIB_M = 2) events.Op1.Barrier <- c(B_CBM = 10, B_RCB = 20, B_GBM = 12, B_RGB = 22, END = lifetime) duration.ev.Barrier <- c(B_CBM = 1, B_RCB = 2, B_GBM = 2, B_RGB = 7) design.Options.Bridge <- list(design.Op1.Bridge <- "Steel girder & precast concrete deck") design.Options.Pylon <- list(design.Op1.Pylon <- "Reinforced concrete cross section") design.Options.Tunnel <- list(design.Op1.Tunnel <- "Steel-Concrete Hybrid deck") design.Options.Pavement <- list(design.Op1.Pavement <- "Recycling Asphalt") design.Options.Barrier <- list(design.Op1.Barrier <- "GFRC Guard Barriers and Steel Crash Barrier") maintenance.Bridge <- dist.Events(lifetime, events.Op1.Bridge, start.Date, design.Options.Bridge) maintenance.Pylon <- dist.Events(lifetime, events.Op1.Pylon, start.Date, design.Options.Pylon) maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel) maintenance.Pavement <- dist.Events(lifetime, events.Op1.Pavement, start.Date, design.Options.Pavement) maintenance.Barrier <- dist.Events(lifetime, events.Op1.Barrier, start.Date, design.Options.Barrier) integrated.interv <- combine.lifeTimelines(maintenance.Bridge, maintenance.Pylon, maintenance.Tunnel, maintenance.Pavement, maintenance.Barrier, duration.ev.Bridge, duration.ev.Pylon, duration.ev.Tunnel, duration.ev.Pavement, duration.ev.Barrier) sum(integrated.interv$duration) getmode <- function(v) { uniqv <- unique(v) uniqv[which.max(tabulate(match(v, uniqv)))] } design.explore <- function(events1, events2, events3,events4,events5) { results <- c() for(i in 1: dim(events1)[1]) { ev1 <- unlist(events1[i, ]) dist.1 <- dist.Events(lifetime, ev1, start.Date, design.Options.Bridge) 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.Pylon) 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.Tunnel) dur.ev3 <- ev3/2 for (l in 1: dim(events1)[1]) { ev4 <- unlist(events4[l, ]) dist.4 <- dist.Events(lifetime, ev4, start.Date, design.Options.Pavement) dur.ev4 <- ev4/2 for (m in 1: dim(events1)[1]) { ev5 <- unlist(events5[m, ]) dist.5 <- dist.Events(lifetime, ev5, start.Date, design.Options.Barrier) dur.ev5 <- ev5/2 combined.lifetime <- combine.lifeTimelines(dist.1, dist.2,dist.3, dist.4,dist.5, dur.ev1, dur.ev2, dur.ev3, dur.ev4, dur.ev5) min.dist.int <- getmode(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, ev5, dur = sum(combined.lifetime$duration), dist.inter = min.dist.int)) } } } } } return(as.data.frame(results)) } n.grid <- 3 events.grid.Bridge <- expand.grid(PD_M = sample(seq(3,28, by = 1), n.grid), PD_R = sample(seq(30,40,1), n.grid)) events.grid.Pylon <- expand.grid(P_M = sample(seq(20, 29, by = 1), n.grid), P_R = sample(seq(43,58,1), n.grid), C_A = sample(seq(7, 23, 1), n.grid), C_R = sample(seq(16, 25, 1), n.grid)) events.grid.Tunnel <- expand.grid(SC_M = sample(seq(5, 25, by = 1), n.grid), SC_I = sample(seq(2,15,1), n.grid)) events.grid.Pavement <- expand.grid(AS_R = sample(seq(15, 25, by = 1), n.grid), AIB_R = sample(seq(40,60,1), n.grid), AS_M = sample(seq(5, 21, 1), n.grid), AIB_M = sample(seq(20,30,1), n.grid)) events.grid.Barrier <- expand.grid(B_CBM = sample(seq(5, 43, by = 1), n.grid), B_RCB = sample(seq(4,25,1), n.grid), B_GMB = sample(seq(5,38,1), n.grid), B_RGB = sample(seq(4,30,1), n.grid)) response.space <- design.explore(events.grid.Bridge, events.grid.Pylon, events.grid.Tunnel, events.grid.Pavement, events.grid.Barrier) 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) #worst opt show_front(p) #=================================== 2. Life Cycle Inventory and Analysis of Integrated Civil Systems ================================== setwd(dirname(getActiveDocumentContext()$path)) LCI.materials <- read.csv("LCImaterialsGroup5.csv") #______________________________________parameter CB.Option <- 'Steel' GB.Option <- 'GFRC' length_bridge = 200 width_bridge = 15 height_deck = 0.3 # thicknes_asphalt = 0.12 girder_Option = 'steel' deck_Option = 'PRC' thicknes_section = 0.04 diameter_Rebar = 0.025 amount_Rebar_nos = 24 diameterSRF_tie = 0.01 distanceSRF_tie = 0.1 name_des_opt = 'RConcrete' H_Slab = 0.2 A_Beam = 0.0156 BeamMaterial_Tun = 4 Reinforcment_Tun = 1 Roh_Tun = 1 length_tunnel = 100 length_all = length_tunnel+length_bridge road_surface = "asphaltSurface" road_base = "asphaltBase" thickness_asphaltSurface <- 0.04 #units m thickness_asphaltBase <- 0.18 #units m LCA.Bridge <- function(length, width, height, thickness, girder.Option, deck.Option, materials, dist.event) { #15 deck.section <- 3 steel.girders.unitWeight <- 8713.5 steel.deck.unitWeight <- 9420 # asphalt.Q <- length * width * thickness materials.split <- split(materials, materials$scope) ###NEW s <- summary(as.factor(dist.event$event)) s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s)) if (deck.Option == "steel") { deck.volume <- length * steel.deck.unitWeight interventions.deck <- 1 } else if(deck.Option == "RC") { deck.volume <- length * width * height interventions.deck <- 2.5 } else if (deck.Option == "PRC") { deck.volume <- 0.67 * length * width * height interventions.deck <- s2$Construction + 0.88*s2$PD_M + 0.13*s2$PD_R } if (girder.Option == "steel") { n <- 1 girders.V <- n * steel.girders.unitWeight * length interventions.girders <- 2 } # asphalt <- mutate(materials.split$asphalt, bridge.Q = asphalt.Q, interventions = 12) deck <- mutate(materials.split[[deck.Option]], bridge.Q = deck.volume, interventions = interventions.deck) if (!is.null(materials.split[[girder.Option]])) { girders <- mutate(materials.split[[girder.Option]], bridge.Q = girders.V, interventions = interventions.girders) LCA.matrix <- rbind(deck, girders) } else { LCA.matrix <- rbind(deck) } #22 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) #23 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) } LCA.Pylon <- function(length, width, thickness, diameter, amountRF, diameterSRF, distanceSRF, name, LCI.materials) { dens_steel = 7.850 # density of steel in t/m³ dens_concrete = 2.4 # density of concrete in t/m³ row_steel = 1 row_RF = 2 row_concrete = 3 if (name == "Steel") { # Design Option 1 selected # Calculation of cross section area for a hollow steel cross section volume <- length*width*width-length*(width-thickness)*(width-thickness) # 5% more for longitudinale ribs and other stiffening elements percent <- 5 volumeSteel <- (1 + percent/100) * volume # Per m3 energyPerm3 <- LCI.materials[row_steel,4] CO2Perm3 <- LCI.materials[row_steel,5] NOXPerm3 <- LCI.materials[row_steel,6] SO2Perm3 <- LCI.materials[row_steel,7] # Calculation of the energy amount (Life Stages: A1-A3) Energy <- volumeSteel * dens_steel * energyPerm3 CO2 <- volumeSteel * CO2Perm3 NOx <- volumeSteel * NOXPerm3 SO2 <- volumeSteel * SO2Perm3 } else if (name == "RConcrete") { # Design Option 2 selected # REINFORCEMENT crossareaSteel <- amountRF * (diameter/2)^2 * pi volumeSteelLong <- crossareaSteel * length # 10mm/10mm # Concrete Covering covering <- 50*1e-3 # in [m] # Adding Stirrup Reinforcement stirrupRFVolume <- (width-covering)*4*pi*(diameterSRF/2)^2 SRFVolume <- stirrupRFVolume*(1/distanceSRF)*length # CONCRETE volumeConcrete <- length * width * width - volumeSteelLong - SRFVolume # Per m3 for Reinforcement energyPerm3_S <- LCI.materials[row_RF,4] CO2Perm3_S <- LCI.materials[row_RF,5] NOxPerm3_S <- LCI.materials[row_RF,6] SO2Perm3_S <- LCI.materials[row_RF,7] # Per m3 for Concrete energyPerm3_C <- LCI.materials[row_concrete,4] CO2Perm3_C <- LCI.materials[row_concrete,5] NOxPerm3_C <- LCI.materials[row_concrete,6] SO2Perm3_C <- LCI.materials[row_concrete,7] # Calculation of the energy amount (Life Stages: A1-A3) Energy <- (volumeSteelLong+SRFVolume) * dens_steel * energyPerm3_S + volumeConcrete * dens_concrete * energyPerm3_C CO2 <- (volumeSteelLong+SRFVolume) * CO2Perm3_S + volumeConcrete * CO2Perm3_C NOx <- (volumeSteelLong+SRFVolume) * NOxPerm3_S + volumeConcrete * NOxPerm3_C SO2 <- (volumeSteelLong+SRFVolume) * SO2Perm3_S + volumeConcrete * SO2Perm3_C } else if (name == "SComposite") { # Design Option 3 selected # HEB 800 A_HEB800 <- 334*1e-4 # in [m²] volumeSteel <- A_HEB800 * length # in [m³] # REINFORCEMENT # Assumption: Since HEB800 will also carry a lot of stress -> diameter of reinforcement can be less diameter <- 16*1e-3 crossareaSteel <- amountRF * (diameter/2)^2 * pi volumeSteelLong <- crossareaSteel * length # Concrete Covering covering <- 5*1e-2 # in [m] # Adding Stirrup Reinforcement stirrupRFVolume <- (width-covering)*4*pi*(diameterSRF/2)^2 SRFVolume <- stirrupRFVolume*(1/distanceSRF)*length # Volume of Concrete volumeConcrete <- length * (width * width - A_HEB800 - volumeSteelLong - SRFVolume) # Per m3 energyPerm3_S <- LCI.materials[row_steel,4] CO2Perm3_S <- LCI.materials[row_steel,5] NOxPerm3_S <- LCI.materials[row_steel,6] SO2Perm3_S <- LCI.materials[row_steel,7] # Per m3 energyPerm3_RF <- LCI.materials[row_RF,4] CO2Perm3_RF <- LCI.materials[row_RF,5] NOxPerm3_RF <- LCI.materials[row_RF,6] SO2Perm3_RF <- LCI.materials[row_RF,7] # Per m3 energyPerm3_C <- LCI.materials[row_concrete,4] CO2Perm3_C <- LCI.materials[row_concrete,5] NOxPerm3_C <- LCI.materials[row_concrete,6] SO2Perm3_C <- LCI.materials[row_concrete,7] # Calculation of the energy amount (Life Stages: A1-A3) Energy <- volumeSteel * dens_steel * energyPerm3_S + (volumeSteelLong + SRFVolume) * dens_steel * energyPerm3_RF + volumeConcrete * dens_concrete * energyPerm3_C CO2 <- volumeSteel * CO2Perm3_S + (volumeSteelLong + SRFVolume) * CO2Perm3_RF + volumeConcrete * CO2Perm3_C NOx <- volumeSteel * NOxPerm3_S + (volumeSteelLong + SRFVolume) * NOxPerm3_RF + volumeConcrete * NOxPerm3_C SO2 <- volumeSteel * SO2Perm3_S + (volumeSteelLong + SRFVolume) * SO2Perm3_RF + volumeConcrete * SO2Perm3_C } else { # in case that not the correct name has been given to the function -> print Error and Stop stop("Selected material does not exist!") } LCA.results <- list (Energy = Energy, CO2 = CO2, NOx = NOx, SO2 = SO2) return(LCA.results) } LCA.Tunnel <- function(HSlab,ABeam,BeamMaterial,Reinforcment,Roh,length_tun,width) { hight_tun=3.80 # height of tunnel amount = length_tun+1 # amount of steel beams VolumeBeam <- ABeam*(2*hight_tun+width)*amount VolumeSlab <- HSlab*length_tun*(width+2*hight_tun) LCA.results <- list (Energy = VolumeBeam*LCI.materials[BeamMaterial,4]+VolumeSlab*LCI.materials[2,4]+VolumeBeam*Roh*LCI.materials[Reinforcment,4], CO2 = VolumeBeam*LCI.materials[BeamMaterial,5]+VolumeSlab*LCI.materials[2,5]+VolumeBeam*Roh*LCI.materials[Reinforcment,5], SO2 = VolumeBeam*LCI.materials[BeamMaterial,6]+VolumeSlab*LCI.materials[2,6]+VolumeBeam*Roh*LCI.materials[Reinforcment,6], NOx = VolumeBeam*LCI.materials[BeamMaterial,7]+VolumeSlab*LCI.materials[2,7]+VolumeBeam*Roh*LCI.materials[Reinforcment,7]) return (LCA.results) } LCA.Pavement <- function(length, width, roadsurface.Option, roadbase.option,thickness.asphaltSurface,thickness.asphaltBase, materials, dist.event) { materials.split <- split(materials, materials$scope) ##NEW 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 asphalt Surface if(roadsurface.Option == "asphaltSurface") { roadsurface.volume <- length * width * thickness.asphaltSurface interventions.roadsurface <- s2$Construction + 0.29*s2$AS_R + 0.79*s2$AS_M } # calc. vol. asphalt Base if (roadbase.option == "asphaltBase"){ roadbase.volume <- length * width * thickness.asphaltBase interventions.roadbase <- s2$Construction + 0.33*s2$AIB_R + 0.71*s2$AIB_M } surface <- mutate(materials.split[[roadsurface.Option]], pave.Q = roadsurface.volume, interventions = interventions.roadsurface) base <- mutate(materials.split[[roadbase.option]], pave.Q = roadbase.volume, interventions = interventions.roadbase) LCA.matrix <- rbind(surface,base) LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * pave.Q, materials.LC = TotalMaterials.Q * interventions, Energy.LC = materials.LC * energy, CO2.LC = materials.LC * CO2 , NOx.LC = materials.LC * NOx , SO2.LC = materials.LC * SO2 ) LCA.results <- list(Energy = sum(LCA.matrix$Energy.LC*2,5), CO2 = sum(LCA.matrix$CO2.LC*1000), NOx = sum(LCA.matrix$NOx.LC*1000), SO2 = sum(LCA.matrix$SO2.LC*1000)) return(LCA.results) } LCA.Barrier <- function(CB.Option, GB.Option, materials,dist.event) { unit.quantity <- function(component,material){ if (component=='CB'& material=='PVC'){ areaExtrudable=0.110911;noOfExElementsPerMeter=1;lenghtNonExtrudable=0;areaNonExtrudable=0;noOfNonExElementsPerMeter=0;SteelFactor=1;PVCFactor=1500; } else if (component=='CB'& material=='Steel'){ areaExtrudable=0.000660;noOfExElementsPerMeter=1;lenghtNonExtrudable=0.500;areaNonExtrudable=0.000660;noOfNonExElementsPerMeter=2;SteelFactor=7850;PVCFactor=1; } else if (component=='GB'& material=='Steel'){ areaExtrudable=0.000212;noOfExElementsPerMeter=3;lenghtNonExtrudable=2.200;areaNonExtrudable=0.001323;noOfNonExElementsPerMeter=5;SteelFactor=7850;PVCFactor=1; } else if (component=='GB'& material=='GFRC'){ areaExtrudable=0.498840;noOfExElementsPerMeter=1;lenghtNonExtrudable=0;areaNonExtrudable=0;noOfNonExElementsPerMeter=0;SteelFactor=1;PVCFactor=1; } else if (component=='GB'& material=='RCC'){ areaExtrudable=0.757490;noOfExElementsPerMeter=1;lenghtNonExtrudable=0;areaNonExtrudable=0;noOfNonExElementsPerMeter=0;SteelFactor=1;PVCFactor=1; } else{ areaExtrudable=0;noOfExElementsPerMeter=0;lenghtNonExtrudable=0;areaNonExtrudable=0;noOfNonExElementsPerMeter=0;SteelFactor=1;PVCFactor=1; } UnitQuantity <- (1*areaExtrudable*noOfExElementsPerMeter + lenghtNonExtrudable*areaNonExtrudable*noOfNonExElementsPerMeter)*SteelFactor*PVCFactor return(UnitQuantity) } 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(GB.Option == "Steel") { GB.quantity <- unit.quantity("GB",GB.Option) interventions.GB <- 5 } else if (GB.Option == "GFRC") { GB.quantity <- unit.quantity("GB",GB.Option) interventions.GB <- s2$Construction + 0.67*s2$B_GBM + 0.33*s2$B_RGB } else if (GB.Option == "RCC") { GB.quantity <- unit.quantity("GB",GB.Option) interventions.GB <- 2.5 } if (CB.Option == "PVC") { CB.quantity <- unit.quantity("CB",CB.Option) interventions.CB <- 3.33 } else if (CB.Option == "Steel") { CB.quantity <- unit.quantity("CB",CB.Option) interventions.CB <- s2$Construction + 0.67*s2$B_CBM + 0.33*s2$B_RCB } else if (CB.Option == "none") { CB.quantity <- 0 interventions.CB <- 0 } GB <- mutate(materials.split[[GB.Option]], viaduct.Q = GB.quantity, interventions = interventions.GB) ###when CB needed if (!is.null(materials.split[[CB.Option]])) { CB <- mutate(materials.split[[CB.Option]], viaduct.Q = CB.quantity, interventions = interventions.CB) LCA.matrix <- rbind(GB, CB) } ###when CB NOT needed else { LCA.matrix <- rbind(GB) } LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * viaduct.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) } Option.Bridge <- LCA.Bridge(length_bridge, width_bridge, height_deck, thicknes_asphalt, girder_Option, deck_Option, LCI.materials, maintenance.Bridge) Option.Pylon <- LCA.Pylon(length_bridge, width_bridge, thicknes_section, diameter_Rebar, amount_Rebar_nos, diameterSRF_tie, distanceSRF_tie, name_des_opt, LCI.materials) Option.Tunnel <- LCA.Tunnel(H_Slab,A_Beam,BeamMaterial_Tun, Reinforcment_Tun, Roh_Tun, length_tunnel, width_bridge) Option.Pavement <- LCA.Pavement(length_all, width_bridge, road_surface, road_base,thickness_asphaltSurface, thickness_asphaltBase, LCI.materials, maintenance.Pavement) Option.Barrier <- LCA.Barrier(CB.Option,GB.Option,LCI.materials, maintenance.Barrier) integrated.Design <- as.data.frame(list(Energy = Option.Bridge$Energy + Option.Pylon$Energy+ Option.Tunnel$Energy+ Option.Pavement$Energy+ Option.Barrier$Energy, CO2 = Option.Bridge$CO2 + Option.Pylon$CO2+ Option.Tunnel$CO2+ Option.Pavement$CO2+ Option.Barrier$CO2, NOx = Option.Bridge$NOx + Option.Pylon$NOx+ Option.Tunnel$NOx+ Option.Pavement$NOx+ Option.Barrier$NOx, SO2 = Option.Bridge$SO2 + Option.Pylon$SO2+ Option.Tunnel$SO2+ Option.Pavement$SO2+ Option.Barrier$SO2)) #=================================== Contribution of subsystems(Additional) ================================== # results_contri <- Option.Bridge # results_contri <- rbind.data.frame(results_contri,Option.Pylon,Option.Tunnel,Option.Pavement,Option.Barrier) # # row.names(results_contri) <-c('Bridge','Pylon','Tunnel',"Pavement","Barrier") # # # barplot(main="Energy",results_contri$Energy, names.arg = row.names(results_contri)) # barplot(main="CO2",results_contri$CO2, names.arg = row.names(results_contri)) # barplot(main="NOx",results_contri$NOx, names.arg = row.names(results_contri)) # barplot(main="SO2",results_contri$SO2, names.arg = row.names(results_contri)) 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) #============================================== 3. Multi-Objective Optimization ============================================ fitness <- function(x) { z <- numeric(7) x <- round(x, 0) y <- expand.grid(PD_M= x[1], PD_R = x[2], P_M = x[3], P_R = x[4], C_A = x[5], C_R = x[6], SC_M = x[7], SC_I = x[8], AS_R = x[9], AIB_R = x[10], AS_M = x[11], AIB_M = x[12], B_CBM = x[13], B_RCB = x[14], B_GMB = x[15], B_RGB = x[16]) dur.ev1 <- unlist(y[1:4] / 2) dur.ev2 <- unlist(y[5:6] / 2) dur.ev3 <- unlist(y[7:8] / 2) dur.ev4 <- unlist(y[9:12] / 2) dur.ev5 <- unlist(y[13:16] / 2) dist.1 <- apply(y[1:4], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.2 <- apply(y[5:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.3 <- apply(y[7:8], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.4 <- apply(y[9:12], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.5 <- apply(y[13:16], 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]],dist.5[[1]], dur.ev1, dur.ev2, dur.ev3, dur.ev4, dur.ev5) z[1] <- sum(results[["duration"]]) z[2] <- -getmode(abs(results$frequency[1:(length(results$frequency) - 1)] - results$frequency[2:length(results$frequency)])) CB.Option <- 'Steel' GB.Option <- 'GFRC' length_bridge = x[17] width_bridge = x[18] height_deck = 0.3 thicknes_asphalt = 0.12 girder_Option = 'steel' deck_Option = 'PRC' thicknes_section = 0.04 diameter_Rebar = 0.025 amount_Rebar_nos = 24 diameterSRF_tie = 0.01 distanceSRF_tie = 0.1 name_des_opt = 'RConcrete' H_Slab = 0.2 A_Beam = 0.0156 BeamMaterial_Tun = 4 Reinforcment_Tun = 1 Roh_Tun = 1 length_tunnel=x[19] length_all = length_bridge + length_tunnel road_surface = "asphaltSurface" road_base = "asphaltBase" thickness_asphaltSurface <- 0.04 #units m thickness_asphaltBase <- 0.18 #units m Energy.costs <- 0.128 CO2.unitcost <- 26 NOx.unitCost <- 42 SO2.unitCosts <- 85 product.1.Bridge <- LCA.Bridge(length_bridge, width_bridge, height_deck, thicknes_asphalt, girder_Option, deck_Option, LCI.materials, maintenance.Bridge) product.2.Pylon <- LCA.Pylon(length_bridge, width_bridge, thicknes_section, diameter_Rebar, amount_Rebar_nos, diameterSRF_tie, distanceSRF_tie, name_des_opt, LCI.materials) product.3.Tunnel <- LCA.Tunnel(H_Slab,A_Beam,BeamMaterial_Tun, Reinforcment_Tun, Roh_Tun, length_tunnel, width_bridge) product.4.Pavement <- LCA.Pavement(length_all, width_bridge, road_surface, road_base,thickness_asphaltSurface, thickness_asphaltBase, LCI.materials,maintenance.Pavement) product.5.Barrier <- LCA.Barrier(CB.Option,GB.Option,LCI.materials,maintenance.Barrier) integrated.system <- as.data.frame(list(Energy = product.1.Bridge$Energy + product.2.Pylon$Energy + product.3.Tunnel$Energy + product.4.Pavement$Energy + product.5.Barrier$Energy, CO2 = product.1.Bridge$CO2 + product.2.Pylon$CO2 + product.3.Tunnel$CO2 + product.4.Pavement$CO2 + product.5.Barrier$CO2, NOx = product.1.Bridge$NOx + product.2.Pylon$NOx + product.3.Tunnel$NOx + product.4.Pavement$NOx + product.5.Barrier$NOx, SO2 = product.1.Bridge$SO2 + product.2.Pylon$SO2 + product.3.Tunnel$SO2 + product.4.Pavement$SO2 + product.5.Barrier$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 = 19, odim = 7, generations=10, popsize=100, lower.bounds=c(5,10,12,15,10,20,10,7,15,30,5,20,5,10,5,17, 180, 12, 80), upper.bounds=c(8,40,29,58,13,25,30,10,20,50,8,25,10,20,10,20, 220, 17, 120)) 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('PD-M','PD-R','P-M','P-R', 'C-A','C-R','SC-M','SC-I', 'AS-R','AIB-R','AS-M','AIB-M', 'B-CBM','B-RCB','B-GBM','B-RGB','B_Len','B_Wid','T_len', "duration", "interv.dist", "energy", "co2", "nox", "so2", "cost", "pareto") par(mfrow = c(1,1)) parcoord(all.results[, 1:24], 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)) #============================================== 4. AHP(Additional) ============================================ # # results_DurInt = pareto3[1:5,1:2] # # # row.names(results_DurInt) <-c('Strategy 1','Strategy 2','Strategy 3',"Strategy 4","Strategy 5") # # results_AHP = pareto3[1:5,3:6] # # # row.names(results_AHP) <-c('Strategy 1','Strategy 2','Strategy 3',"Strategy 4","Strategy 5") # # write.csv(results_AHP,"results_AHP.csv", row.names = FALSE) # # # barplot(main="Duration",results_DurInt$duration, names.arg = row.names(results_DurInt)) # barplot(main="Average Distance between Intervention",results_DurInt$interv.dist, names.arg = row.names(results_DurInt)) # barplot(main="Energy",results_AHP$energy, names.arg = row.names(results_AHP)) # barplot(main="CO2",results_AHP$co2, names.arg = row.names(results_AHP)) # barplot(main="NOx",results_AHP$nox, names.arg = row.names(results_AHP)) # barplot(main="SO2",results_AHP$so2, names.arg = row.names(results_AHP)) # # results_AHP.df = results_AHP %>% # rownames_to_column(var = "option") %>% # mutate_each(funs(rescale), -option) %>% # melt(id.vars = c('option'), measure.vars = colnames(results_AHP)) %>% # arrange(option) # results_AHP.df %>% # ggplot(aes(x=option, y=value, group=variable, color=variable)) + # geom_polygon(fill=NA) + # coord_polar() + theme_bw() + facet_wrap(~ variable) + # theme(axis.text.x = element_text(size = 5)) # # # alternatives <- rownames(results_AHP) # indicators <- colnames(results_AHP) # # energy <- t(matrix(c( 1, 2, 1/5, 1/4, 1/3, # 1/2, 1, 1/5, 1/5, 1/3, # 5, 5, 1, 2, 4, # 4, 5, 1/2, 1, 3, # 3, 3, 1/4, 1/3, 1), # nrow=length(alternatives),ncol=length(alternatives), # dimnames = list(alternatives, alternatives))) # # energy == t(1/energy) # # CO2 <- t(matrix(c(1, 1/3, 1/2, 1/5, 1/5, # 3, 1, 2, 1/3, 1/2, # 2, 1/2, 1, 1/5, 1/4, # 5, 3, 5, 1, 2, # 5, 2, 4, 1/2, 1), # nrow=length(alternatives),ncol=length(alternatives), # dimnames = list(alternatives, alternatives))) # CO2 == t(1/CO2) # # NOx <- t(matrix(c(1, 1/4, 2, 1/3, 1/5, # 4, 1, 4, 2, 1/2, # 1/2, 1/4, 1, 1/3, 1/5, # 3, 1/2, 3, 1, 1/2, # 5, 2, 5, 2, 1), # nrow=length(alternatives),ncol=length(alternatives), # dimnames = list(alternatives, alternatives))) # NOx == t(1/NOx) # # SO2 <- t(matrix(c(1, 1/2, 1/2, 1/5, 1/4, # 2, 1, 2, 1/2, 1/2, # 2, 1/2, 1, 1/4, 1/3, # 5, 2, 4, 1, 2, # 4, 2, 3, 1/2, 1), # nrow=length(alternatives),ncol=length(alternatives), # dimnames = list(alternatives, alternatives))) # SO2 == t(1/SO2) # # APCL <- list(energy=energy, # CO2=CO2, # NOx=NOx, # SO2 = SO2) # # CWPC <- t(matrix(c( 1, 2, 4, 3, # 1/2, 1, 3, 2, # 1/4, 1/3, 1, 1/2, # 1/3, 1/2, 2, 1 ), # nrow=length(APCL),ncol=length(APCL), # dimnames = list(names(APCL), names(APCL)))) # CWPC == t(1/CWPC) # # # results_AHP.AHP <- AHP(CWPC, APCL) # # results_AHP.AHP # # plot.new() # piepercent<- round(100*results_AHP.AHP, 1) # pie(results_AHP.AHP, labels = piepercent, # col = c("lightblue", "mistyrose", "lightcyan","lavender", # "cornsilk", "thistle4", "steelblue4"), # radius = 1, # main = "Ranking of the design options using AHP") # legend("bottomleft", # names(results_AHP.AHP), cex = 0.8, # fill = c("lightblue", "mistyrose", "lightcyan","lavender", # "cornsilk", "thistle4", "steelblue4"))