library(timelineS) library(lubridate) library(ggplot2) library(dplyr) library(reshape2) library(rPref) library(mco) library(MASS) rm (list = ls()) ############# 1. Integrated Maintenance Planning ############## #Representation of life cycle of the subsystems 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) } #Function for different specific maintenance actions of each subsystem 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) } #Definition of function to combine life cycle timelines of four subsystems + identify their overlap 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])] #all are not equal to zero if(length(phase.1) != 0 && length(phase.2) != 0 && length(phase.3) != 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.1, phase.2, phase.3, phase.4) 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)]) } #1 is zero else if(length(phase.2) != 0 && length(phase.1) == 0 && length(phase.3) != 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.2, phase.3, phase.4) base$duration[index] <- max(p2.dur[which(names(p2.dur) == phase.2)], p3.dur[which(names(p3.dur) == phase.3)], p4.dur[which(names(p4.dur) == phase.4)]) } #2 is zero else if(length(phase.1) != 0 && length(phase.2) == 0 && length(phase.3) != 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.1, phase.3, phase.4) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p3.dur[which(names(p3.dur) == phase.3)], p4.dur[which(names(p4.dur) == phase.4)]) } #3 is zero else if(length(phase.1) != 0 && length(phase.3) == 0 && length(phase.2) != 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.1, phase.2, phase.4) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p2.dur[which(names(p2.dur) == phase.2)], p4.dur[which(names(p4.dur) == phase.4)]) } #4 is zero else if(length(phase.1) != 0 && length(phase.4) == 0 && length(phase.3) != 0 && length(phase.2) != 0) { base$Names[index] <- paste(phase.1, phase.3, phase.2) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p3.dur[which(names(p3.dur) == phase.3)], p2.dur[which(names(p2.dur) == phase.2)]) } #1 and 2 are not zero else if(length(phase.1) != 0 && length(phase.2) != 0 && length(phase.3) == 0 && length(phase.4) == 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)]) } #1 and 3 are not zero else if(length(phase.1) != 0 && length(phase.3) != 0 && length(phase.2) == 0 && length(phase.4) == 0) { base$Names[index] <- paste(phase.1, phase.3) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p3.dur[which(names(p3.dur) == phase.3)]) } #1 and 4 are not zero else if(length(phase.1) != 0 && length(phase.4) != 0 && length(phase.3) == 0 && length(phase.2) == 0) { base$Names[index] <- paste(phase.1, phase.4) base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)], p4.dur[which(names(p4.dur) == phase.4)]) } #2 and 3 are not zero else if(length(phase.1) == 0 && length(phase.2) != 0 && length(phase.3) != 0 && length(phase.4) == 0) { base$Names[index] <- paste(phase.2, phase.3) base$duration[index] <- max(p2.dur[which(names(p2.dur) == phase.2)], p3.dur[which(names(p3.dur) == phase.3)]) } #2 and 4 are not zero else if(length(phase.1) == 0 && length(phase.2) != 0 && length(phase.3) == 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.2, phase.4) base$duration[index] <- max(p2.dur[which(names(p2.dur) == phase.2)], p4.dur[which(names(p4.dur) == phase.4)]) } #3 and 4 are not zero else if(length(phase.1) == 0 && length(phase.2) == 0 && length(phase.3) != 0 && length(phase.4) != 0) { base$Names[index] <- paste(phase.3, phase.4) base$duration[index] <- max(p3.dur[which(names(p3.dur) == phase.3)], p4.dur[which(names(p4.dur) == phase.4)]) } # 1 is not zero else if (length(phase.2) == 0 && length(phase.1) != 0 && length(phase.3) == 0 && length(phase.4) == 0 ) { base$Names[index] <- phase.1 base$duration[index] <- p1.dur[which(names(p1.dur) == phase.1)] } # 2 is not zero else if (length(phase.1) == 0 && length(phase.2) != 0 && length(phase.3) == 0 && length(phase.4) == 0 ) { base$Names[index] <- phase.2 base$duration[index] <- p2.dur[which(names(p2.dur) == phase.2)] } # 3 is not zero else if (length(phase.1) == 0 && length(phase.3) != 0 && length(phase.2) == 0 && length(phase.4) == 0 ) { base$Names[index] <- phase.3 base$duration[index] <- p3.dur[which(names(p3.dur) == phase.3)] } # 4 is not zero else if (length(phase.1) == 0 && length(phase.4) != 0 && length(phase.3) == 0 && length(phase.2) == 0 ) { base$Names[index] <- phase.4 base$duration[index] <- p4.dur[which(names(p4.dur) == phase.4)] } else { base$duration[index] <- 0 } } return(base) } par(mfrow = c(2, 1)) start.Date = "2020-01-01" lifetime <- 50 #chosen intervention of 50 years, since the the Heating system has a lifetime of 5 years, so the interventions can be multiplied by 10 (easier to calculate) #Events and Duration of each specific intervention for every subsystem #1. Subsystem: Frame events.Op1.Frame = c(SFO = 9, M.F = 3, F.R = 15, END = lifetime) #lifetime of subsystem of 100 years duration.ev.Frame <- c(SFO = 3.5, M.F = 1, F.R = 10.5) #2. Subsystem: Xpelair events.Op2.Xpelair = (c(FR = 1, FM = 2, END = lifetime)) #lifetime of subsystem of 5 years duration.ev.Xpelair <- (c(FR = 2, FM = 1)) #3. Subsystem: Stellio events.Op3.Stellio = c(M.S = 4, RM = 18, RMC=14, END = lifetime) #lifetime of subsystem of 30 years duration.ev.Stellio <- c(M.S = 1, RM = 15, RMC = 7) #4. Subsystem: Isolated Footing events.Op1.Isolated.Footing = c(M.I = 8, G= 12, RI = 35, END = lifetime) #lifetime of subsystem of 70 years duration.ev.Isolated.Footing <- c(M.I = 2, G = 3, RI = 28) #Define design options for each system design.Options.Frame <- list(desing.Op1.Frame <- "1.Wood frame standing concrete foundation") design.Options.Xpelair <- list(desing.Op2.Xpelair <- "2.Xpelair") design.Options.Stellio <- list(desing.Op3.Stellio <- "3. Stellio") design.Options.Isolated.Footing <- list(desing.Op1.Isolated.Footing <- "4. Isolated Footing") #Generate distribution of interventions over lifetime maintenance.Frame <- dist.Events(lifetime, events.Op1.Frame, start.Date, design.Options.Frame) maintenance.Xpelair <- dist.Events(lifetime, events.Op2.Xpelair, start.Date, design.Options.Xpelair) maintenance.Stellio <- dist.Events(lifetime, events.Op3.Stellio, start.Date, design.Options.Stellio) maintenance.Isolated.Footing <- dist.Events(lifetime, events.Op1.Isolated.Footing, start.Date, design.Options.Isolated.Footing) #Combination of all four timelines integrated.interv <- combine.lifeTimelines(maintenance.Frame, maintenance.Xpelair, maintenance.Stellio, maintenance.Isolated.Footing, duration.ev.Frame, duration.ev.Xpelair, duration.ev.Stellio, duration.ev.Isolated.Footing) #Sum up total duration of interruptions for the set lifetime of 50 years sum(integrated.interv$duration) # in case we have another maintenance strategy #events.Op1.Frame = c(SFO=18, M.F = 6, F.R = 30, END = lifetime) #duration.ev.Frame <- c(SFO = 7, M.F = 2, F.R = 21) #events.Op2.Xpelair <- c(FR= 2, FM = 2, END = lifetime) #actually every 12 and 6 months #duration.ev.Xpelair <- c(FR = 2, FM = 1) #events.Op3.Stellio = c(C = 1, M.S = 5, RM = 25, RMC=20, END = lifetime) #this system has a lifetime of 30 years #duration.ev.Stellio <- c(C = 7, M.S = 2, RM = 21, RMC = 10) #events.Op1.Isolated.Footing <- c(M.I = 8, G= 12, RI = 35, END = lifetime) #duration.ev.Isolated.Footing <- c(M.I = 7, G = 2, RI = 21) #design.Options.Frame <- list(desing.Op1.Frame <- "1.Wood frame standing concrete foundation") #design.Options.Xpelair <- list(desing.Op2.Xpelair <- "2.Xpelair") #design.Options.Stellio <- list(desing.Op3.Stellio <- "3. Stellio") #design.Options.Isolated.Footing <- list(desing.Op1.Isolated.Footing <- "4. Isolated Footing") #maintenance.Frame <- dist.Events(lifetime, events.Op1.Frame, start.Date, design.Options.Frame) #maintenance.Xpelair <- dist.Events(lifetime, events.Op2.Xpelair, start.Date, design.Options.Xpelair) #maintenance.Stellio <- dist.Events(lifetime, events.Op3.Stellio, start.Date, design.Options.Stellio) #maintenance.Isolated.Footing <- dist.Events(lifetime, events.Op1.Isolated.Footing, start.Date, design.Options.Isolated.Footing) #integrated.interv <- combine.lifeTimelines(maintenance.Frame, maintenance.Xpelair, maintenance.Stellio, maintenance.Isolated.Footing, #duration.ev.Frame, duration.ev.Xpelair, duration.ev.Stellio, duration.ev.Isolated.Footing) #sum(integrated.interv$duration) #Function to automate the combination od two systems and their maintenance plans 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, design.Options.Frame) dur.ev1 <- ev1/4 for (j in 1: dim(events2)[1]) { ev2 <- unlist(events2[j, ]) dist.2 <- dist.Events(lifetime, ev2, start.Date, design.Options.Xpelair) dur.ev2 <- ev2/4 for(k in 1: dim(events3)[1]) { ev3 <- unlist(events3[k, ]) dist.3 <- dist.Events(lifetime, ev3, start.Date, design.Options.Stellio) dur.ev3 <- ev3/4 for (l in 1: dim(events4)[1]) { ev4 <- unlist(events4[l, ]) dist.4 <- dist.Events(lifetime, ev4, start.Date, design.Options.Isolated.Footing) dur.ev4 <- ev4/4 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)] -combined.lifetime$frequency[3:length(combined.lifetime$frequency)] -combined.lifetime$frequency[4: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)) } #Create scenario space #Parameter, which controls how the parameter will vary n.grid <- 2 #Intervention matrix scenario for the Frame events.grid.Frame <- expand.grid(SFO = sample(seq(13, 23, by = 1), n.grid), M.F = sample(seq(3, 9, by = 1), n.grid), F.R = sample(seq(25, 35, by = 1), n.grid)) events.grid.Frame #Intervention matrix scenario for the Xpelair events.grid.Xpelair <- expand.grid(FR = sample(seq(1,3, by = 1), n.grid), FM = sample(seq(1,3, by = 1), n.grid)) events.grid.Xpelair #Intervention matrix scenario for the Stellio events.grid.Stellio <- expand.grid(M.S = sample(seq(3, 8,by = 1), n.grid), RM = sample(seq(20, 30, by = 1), n.grid), RMC = sample(seq(15, 25, by = 1), n.grid)) events.grid.Stellio #Intervention matrix scenario for the Isolated Footing events.grid.Isolated.Footing <- expand.grid(M.I = sample(seq(5, 11, by = 1), n.grid), G = sample(seq(7, 17, by = 1), n.grid), RI = sample(seq(30, 40, by = 1), n.grid)) events.grid.Isolated.Footing #Use the function design.explore() for automation of system combination response.space <- design.explore(events.grid.Frame, events.grid.Xpelair, events.grid.Stellio, events.grid.Isolated.Footing) response.space #Minimize the duration of the interruptions and Maximize the distance between interventions p <- low(dur)* high(dist.inter) #Evaluation of preference and SDelection of the alternative, which is the best for preference sky <- psel(response.space, p) sky #Get ranked alternatives pareto2 <- psel(response.space, p, top = nrow(response.space)) #Visualization of classification of the different subsystems ggplot(response.space, aes(x = dur, y = dist.inter)) + geom_point(shape = 21) + geom_point(data = pareto2, size = 3, aes(color = factor(pareto2$.level))) #Function to find the Pareto Frontier 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) #Visualize Pareto Front #blue line represents the Pareto frontier, dark black point represents the best alternative based on defined preferences p2 <- high(dur) * low(dist.inter) show_front(p2) sky2 <- psel(response.space, p2) ############# 2. LCI and Analysis of Integrated Civil Systems ############## LCI.materials <- read.csv("LCImaterials_Group4.csv") ############################## Frame (Alex) ############################## LCA.Frame <- function(length, width, height, Framevolume, Frameinterventions, materials) { length.Frame <- 3.2*length # length 48, total length 15m width.Frame <- 1.67*width #width 25m, total width 25m height <- 1 #units m length of a side of the Frame having a square section steel.bracing.uniWeight <- 35.5 #weight of HEA180 steel profile Frame.volume = length.Frame * width.Frame * height #wood bracing.volume = steel.bracing.uniWeight * (8*sqrt(5*5+15*15)+16*sqrt(5*5+6.25*6.25)) #units m³, cross section of bracing steel x its length in between the column and its length between the girders x its weigth interventions.Frame <- 7 intervention.bracing <- 1 #total consumption/emission Frame.energy <- as.numeric(LCI.materials$energy[2]) * Frame.volume*interventions.Frame + as.numeric(LCI.materials$energy[1])*bracing.volume*intervention.bracing Frame.CO2 <- as.numeric(LCI.materials$CO2[2])*Frame.volume*interventions.Frame + as.numeric(LCI.materials$CO2[1])*bracing.volume*intervention.bracing Frame.NOx <- as.numeric(LCI.materials$NOx[2])*Frame.volume*interventions.Frame + as.numeric(LCI.materials$NOx[1])*bracing.volume*intervention.bracing Frame.SO2 <- as.numeric(LCI.materials$SO2[2])*Frame.volume*interventions.Frame + as.numeric(LCI.materials$SO2[1])*bracing.volume*intervention.bracing ##creation of list with added up energy consumption and emission Frame.LCA.results <- list(Energy = (Frame.energy) , CO2 = (Frame.CO2), NOx = (Frame.NOx), SO2 = (Frame.SO2)) return(Frame.LCA.results) } ############################## Xpelair (Aya) ############################## LCA.Xpelair <- function(length, width, depth, filterVolume, filterInterventions,insulationVolume,insulationInterventions, materials) { # calculate the volume of the cover, steel thickness 0,75mm length.xpelair = length*0.02436 #60.9 cm, total length 25m width.xpelair = width*0.0368 # 55.2 cm, total width 15m depth.xpelair <- 28.5 #units cm filter.volume <- 1.94*2.25*0.96 #polyester filterIntervention <- 40 #multiplied by 10 because the lifetime is 5x10 = 50 years (10 times the lifetime of the system) insulationVolume <- length.xpelair * width.xpelair * 0.12 #rubber insulationIntervention=10 #multiplied by 10 because the lifetime is 5x10 = 50 years (10 times the lifetime of the system) steel.Q <- length.xpelair*width.xpelair*depth.xpelair-((length.xpelair-2*0.075)*(width.xpelair-2*0.075)*(depth.xpelair-2*0.075)) #steel #total consumption/emission Xpelair.energy <- as.numeric(LCI.materials$energy[5])*filter.volume*filterIntervention + as.numeric(LCI.materials$energy[4])*insulationVolume*insulationIntervention + as.numeric(LCI.materials$energy[3])*steel.Q Xpelair.CO2 <- as.numeric(LCI.materials$CO2[5])*filter.volume*filterIntervention + as.numeric(LCI.materials$CO2[4])*insulationVolume*insulationIntervention + as.numeric(LCI.materials$CO2[3])*steel.Q Xpelair.NOx <- as.numeric(LCI.materials$NOx[5])*filter.volume*filterIntervention + as.numeric(LCI.materials$NOx[4])*insulationVolume*insulationIntervention + as.numeric(LCI.materials$NOx[3])*steel.Q Xpelair.SO2 <- as.numeric(LCI.materials$SO2[5])*filter.volume*filterIntervention + as.numeric(LCI.materials$SO2[4])*insulationVolume*insulationIntervention + as.numeric(LCI.materials$SO2[3])*steel.Q ##creation of list with added up energy consumption and emission Xpelair.LCA.results <- list(Energy = (Xpelair.energy) , CO2 = (Xpelair.CO2), NOx = (Xpelair.NOx), SO2 = (Xpelair.SO2)) return(Xpelair.LCA.results) } ############################## Stellio (Stephi) ############################## LCA.Stellio <- function(width, oncentrator, structure, grounding, materials) { #geometrical characteristics - unit weights steel.unitweight <- 1 zinc.unitweight <- 1 PVC.unitweight <- 1 concrete.unitweight <- 1 reinforcement.uniweight <- 1 glass.unitweight <- 1 #total weights concentrator.weight = 11.0 * glass.unitweight #unit: kg structure.weight = 26 * steel.unitweight + 0.13 * zinc.unitweight #unit: kg grounding.weight = 26.0 * concrete.unitweight + 2.20 * reinforcement.uniweight #unit: kg height.Stellio <- width * 0.2 #Stellio heigth of 3m, total heigth = width of 15m Stellio.volume <- 48 * height.Stellio #interventions intervention.concentrator <- 4 intervention.structure <- 9 #total consumption/emission Stellio.energy <- as.numeric(LCI.materials$energy[10])*concentrator.weight*intervention.concentrator + as.numeric(LCI.materials$energy[8])*26 * steel.unitweight*intervention.structure + as.numeric(LCI.materials$energy[9])*0.13 * zinc.unitweight*intervention.structure + as.numeric(LCI.materials$energy[6])*26.0 * concrete.unitweight*intervention.structure + as.numeric(LCI.materials$energy[7])*2.20 * reinforcement.uniweight*intervention.structure Stellio.CO2 <- as.numeric(LCI.materials$CO2[10])*intervention.concentrator*Stellio.volume + as.numeric(LCI.materials$CO2[8])*26 * intervention.structure*Stellio.volume + as.numeric(LCI.materials$CO2[9])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$CO2[6])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$CO2[7])*intervention.structure*Stellio.volume Stellio.NOx <- as.numeric(LCI.materials$NOx[10])*concentrator.weight*intervention.concentrator*Stellio.volume + as.numeric(LCI.materials$NOx[8])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$NOx[9])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$NOx[6])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$NOx[7])*intervention.structure*Stellio.volume Stellio.SO2 <- as.numeric(LCI.materials$SO2[10])*concentrator.weight*intervention.concentrator*Stellio.volume + as.numeric(LCI.materials$SO2[8])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$SO2[9])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$SO2[6])*intervention.structure*Stellio.volume + as.numeric(LCI.materials$SO2[7])**intervention.structure*Stellio.volume Stellio.LCA.results <- list(Energy = (Stellio.energy) , CO2 = (Stellio.CO2), NOx = (Stellio.NOx), SO2 = (Stellio.SO2)) return(Stellio.LCA.results) } ############################## Isolated.Footing (Basil) ############################## LCA.Isolated.Footing<-function(length,width,design.Op,material){ Depth.footing<-1.5 Depth.PCC<-0.1 Factor.footing.rcc<-0.03 Factor.footing.PCC<-0.1 RCC.volume<-length*width*Depth.footing*Factor.footing.rcc intervention<-2 #total consumption/emission RCC RCC.energy <- as.numeric(LCI.materials$energy[11])*RCC.volume*intervention + as.numeric(LCI.materials$energy[12])*RCC.volume*intervention + as.numeric(LCI.materials$energy[13])*RCC.volume*intervention + as.numeric(LCI.materials$energy[14])*RCC.volume*intervention RCC.CO2 <- as.numeric(LCI.materials$CO2[11])*RCC.volume*intervention + as.numeric(LCI.materials$CO2[12])*RCC.volume*intervention + as.numeric(LCI.materials$CO2[13])*RCC.volume*intervention + as.numeric(LCI.materials$CO2[14])*RCC.volume*intervention RCC.NOx <- as.numeric(LCI.materials$NOx[11])*RCC.volume*intervention + as.numeric(LCI.materials$NOx[12])*RCC.volume*intervention + as.numeric(LCI.materials$NOx[13])*RCC.volume*intervention + as.numeric(LCI.materials$NOx[14])*RCC.volume*intervention RCC.SO2 <- as.numeric(LCI.materials$SO2[11])*RCC.volume*intervention + as.numeric(LCI.materials$SO2[12])*RCC.volume*intervention + as.numeric(LCI.materials$SO2[13])*RCC.volume*intervention + as.numeric(LCI.materials$SO2[14])*RCC.volume*intervention PCC.volume<-length*width*Depth.PCC*Factor.footing.PCC intervention.PCC<-2 #total consumption and emission of PCC PCC.energy<-as.numeric(LCI.materials$energy[15])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$energy[16])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$energy[17])*PCC.volume*intervention.PCC PCC.CO2<-as.numeric(LCI.materials$CO2[15])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$CO2[16])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$CO2[17])*PCC.volume*intervention.PCC PCC.NOx<-as.numeric(LCI.materials$NOx[15])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$NOx[16])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$NOx[17])*PCC.volume*intervention.PCC PCC.SO2<-as.numeric(LCI.materials$SO2[15])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$SO2[16])*PCC.volume*intervention.PCC + as.numeric(LCI.materials$SO2[17])*PCC.volume*intervention.PCC ##creation of list with added up energy consumption and emission Isolated.Footing.LCA.results <- list(Energy = (RCC.energy+PCC.energy) , CO2 = (RCC.CO2+PCC.CO2), NOx = (RCC.NOx+PCC.NOx), SO2 = (RCC.SO2+PCC.SO2)) return(Isolated.Footing.LCA.results) } #Define main variables to run the function length <- 25 # units: m width <- 15 #units: m #Application of LCA-functions to calculate the impact on environment Option.Frame <- LCA.Frame(length,width, height, Frame.volume,interventions.Frame, LCI.materials) Option.Xpelair <- LCA.Xpelair(length, width, depth, filterVolume, filterInterventions,insulationVolume,insulationInterventions, LCI.materials) Option.Stellio <- LCA.Stellio(width, concentrator.weight, structure.weight, grounding.weight, LCI.materials) Option.Isolated.Footing <- LCA.Isolated.Footing(length,width,design.Op,material) #Combination of results of all subsystems into one dataframe integrated.Design <- as.data.frame(list(Energy = Option.Frame$Energy + Option.Xpelair$Energy + Option.Stellio$Energy + Option.Isolated.Footing$Energy, CO2 = Option.Frame$CO2 + Option.Xpelair$CO2 + Option.Stellio$CO2 + Option.Isolated.Footing$CO2, NOx = Option.Frame$NOx + Option.Xpelair$NOx + Option.Stellio$NOx + Option.Isolated.Footing$NOx, SO2 = Option.Frame$SO2 + Option.Xpelair$SO2 + Option.Stellio$SO2 + Option.Isolated.Footing$SO2)) integrated.Design #Costs of the impact on the environment - Definition of unit prices Energy.costs <- 0.128 CO2.unitcost <- 26 # per metric tone NOx.unitCost <- 42 # per metric tone SO2.unitCosts <- 85 # per metric tone #Add costs to matrix of integrated designs integrated.Design <- mutate(integrated.Design, Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx * NOx.unitCost + SO2 * SO2.unitCosts)/10^9) integrated.Design ############# 3. Multi-Objective Optimization ############## ##################################### Fitness Function ######################################### fitness <- function(x) { #Definition of output dimension z <- numeric(7) x <- round(x, 0) y <- expand.grid(SFO = x[1], M.F = x[2], F.R = x[3], FR = x[4], FM = x[5], M.S = x[6], RM = x[7], RMC = x[8], M.I = x[9], G = x[10], RI = x[11]) dur.ev1 <- unlist(y[1:3] / 4) dur.ev2 <- unlist(y[4:5] / 4) dur.ev3 <- unlist(y[6:8] / 4) dur.ev4 <- unlist(y[9:11] / 4) dist.1 <- apply(y[1:3], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.2 <- apply(y[4:5], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.3 <- apply(y[6:8], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date) dist.4 <- apply(y[9:11], 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) #Add total number of interruptions and minimum distance between interventions to output vector z[1] <- sum(results[["duration"]]) z[2] <- -min(abs(results$frequency[1:(length(results$frequency) - 1)] - results$frequency[2:length(results$frequency)] - results$frequency[3:length(results$frequency)] - results$frequency[4:length(results$frequency)])) #Set static variables width <- x[12] #based on the value generated by the genetic algorithm length <- 25 # units: m materials <- LCI.materials Energy.costs <- 0.128 # per metric tone CO2.unitcost <- 26 # per metric tone NOx.unitCost <- 42 # per metric tone SO2.unitCosts <- 85 # per metric tone #Call function for LCA for each subsystem product.1.Frame <- LCA.Frame(length, width, LCI.materials, dist.1[[1]]) product.2.Xpelair <- LCA.Xpelair(length, width, LCI.materials, dist.2[[1]]) product.3.Stellio <- LCA.Stellio(length, width, LCI.materials, dist.3[[1]]) product.4.Isolated.Footing <- LCA.Isolated.Footing(length, width, LCI.materials, dist.4[[1]]) #Call function for design integration integrated.system <- as.data.frame(list(Energy = product.1.Frame$Energy + product.2.Xpelair$Energy + product.3.Stellio$Energy + product.4.Isolated.Footing$Energy, CO2 = product.1.Frame$CO2 + product.2.Xpelair$CO2 + product.3.Stellio$CO2 + product.4.Isolated.Footing$CO2, NOx = product.1.Frame$NOx + product.2.Xpelair$NOx + product.3.Stellio$NOx + product.4.Isolated.Footing$CO2, SO2 = product.1.Frame$SO2 + product.2.Xpelair$SO2 + product.3.Stellio$SO2 + product.4.Isolated.Footing$SO2)) options(scipen = 999) #Add column for cost of the impact on environment integrated.system <- mutate(integrated.system, Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx * NOx.unitCost + SO2 * SO2.unitCosts)) #Add results to output of fitness function 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) } #Define GA optimization function r2 <- nsga2(fitness, idim = 12, odim = 7, generations=10, popsize=100, lower.bounds=c(13, 3, 25, 1, 1, 3, 20, 15, 5, 7, 30, 10), upper.bounds= c(23, 9, 35, 3, 3, 8, 30, 25, 11, 17, 40, 25)) #Transform results of to data frame r2Results <- as.data.frame(r2$value) outNames <- c("duration", "interv.dist", "energy", "co2", "nox", "so2", "cost") colnames(r2Results) <- outNames #Extract values, which represent the pareto front and transform these into data frame pareto3 <- as.data.frame(paretoFront(r2)) colnames(pareto3) <- outNames #Plot the Pareto front described by total duration of intervention and costs 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") #How does accumulated impact of input parameters affect the performance criteria? input.params <- round(r2$par, 0) all.results <- cbind(input.params, r2Results, r2$pareto.optimal) colnames(all.results) <- c("SFo", "M.R", "F.R", "FR", "FM", "M.S", "RM", "RMC","M", "G", "RI", "width", "duration", "interv.dist","energy", "co2", "nox", "so2", "cost", "pareto") #interpretation #setting the number of alternatives we have y=length(pareto3$duration) y #Defining the colors of our alternatives col1=colors()[seq(2, 160, by= 160/y)] par(mfrow = c(1,1)) parcoord(all.results[, 1:20], var.label = T, col = ifelse(all.results$pareto == TRUE, col1, "skyblue2"), lty = ifelse(all.results$pareto == TRUE, 1, 3), lwd = ifelse(all.results$pareto == TRUE, 3, 1)) result_Energy = all.results$energy[c(1:y)] result_CO2 = all.results$co2[c(1:y)] result_NOx = all.results$nox[c(1:y)] result_SO2 = all.results$so2[c(1:y)] result_Cost = all.results$cost[c(1:y)] result_duration = all.results$duration[c(1:y)] result_inter_dist = -all.results$interv.dist[c(1:y)] barplot(result_Energy,col=col1, main = "Energy consumption level for differnt options", ylab = "Energy [KJ]", xlab = "alternatives" ) barplot(result_CO2,col=col1, main = "CO2 consumption level for differnt options", ylab = "CO2 [Kg]", xlab = "alternatives" ) barplot(result_NOx,col=col1, main = "NOx consumption level for differnt options", ylab = "NOx [Kg]", xlab = "alternatives" ) barplot(result_SO2,col=col1, main = "SO2 consumption level for differnt options", ylab = "SO2 [Kg]", xlab = "alternatives" ) barplot(result_Cost,col=col1, main = "cost consumption level for differnt options", ylab = "Cost [Euros]", xlab = "alternatives" ) barplot(result_duration, col = col1, main = "Maintainance duration for different options", ylab = "Duration [days]", xlab="alternatives") barplot(result_inter_dist, col = col1, main = "interval distance for different options", ylab = "interval [years]", xlab="alternatives") result_Energy_2 = 100-(all.results$energy[c(1:y)]*(100/max(result_Energy))) result_CO2_2 = 100-(all.results$co2[c(1:y)]*(100/max(result_CO2))) result_NOx_2 = 100-(all.results$nox[c(1:y)]*(100/max(result_NOx))) result_SO2_2 = 100-(all.results$so2[c(1:y)]*(100/max(result_SO2))) result_Cost_2 = 100-(all.results$cost[c(1:y)]*(100/max(result_Cost))) result_duration_2 = 100-(all.results$duration[c(1:y)]*(100/max(result_duration))) result_inter_dist_2 = (all.results$duration[c(1:y)]*(100/max(result_duration))) Total_result = data.frame(result_inter_dist_2,result_Energy_2,result_CO2_2,result_NOx_2,result_SO2_2,result_duration_2,result_Cost_2) Total_result par(mfrow = c(1,1)) parcoord(Total_result, var.label = T, col=col1,lty = 1, lwd = 3,) write.csv(pareto3, "C:/Users/alext/Documents/Masters/S1/LCA/Project3/results.csv")