##############################################################################

# TU Berlin - Civil and Building Systems
# Prof. Timo Hartmann, Lucian Ungureanu
# Life-Cycle Assessment - winter semester 2020/2021
#
# Project assignment 3
# 25.02.2021

##############################################################################

# Group 3

##############################################################################

##############################################################################


#install.packages("timelineS")
#install.packages("dplyr")
#install.packages("MCDA")
#install.packages("scales")
#install.packages("reshape2")
#install.packages("tibble")
#install.packages("lubridate")
#install.packages("rPref")
#install.packages("mco")
#install.packages("MASS")


rm (list = ls())

library(timelineS)
library(dplyr)
library(MCDA)
library(scales)
library(reshape2)
library(tibble)
library(lubridate)
library(ggplot2)
library(rPref)
library(mco)
library(MASS)

#--------------------------------MAINTENANCE INTERVENTIONS--------------------------------------#

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", "events")
  all.events<-all.events[order(all.events$frequency),]
  
  unique.events<-all.events[!duplicated(all.events$frequency),]
  unique.events$event[which(unique.events$frequency==0)]<-"SI"
  
  #plot the timeline
  #this needs to be deactivated for the MOO to run but activated to plot the timelines
  #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] <- "SI"
  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(1,1))
start.Date="2021-01-15"
lifetime<-60


##Elevator
events.Op1.Elevator <- c(RD = 12, RM = 15, RS = 10, RDS = 24,  END = lifetime)
duration.ev.Elevator <- c(RD = 8.8, RM = 5.7, RS = 0.1, RDS = 0.2)

##Sprinkler pipe
events.Op2.Sprinklerpipe <- c(CF=5,PR=15,END=lifetime)
duration.ev.Sprinklerpipe <- c(CF=1.75,PR=3)

##ETICS
events.Op3.ETICS<-c(C = 3, FC = 15, END = lifetime)
duration.ev.ETICS<-c(C = 1, FC = 1)

##Climate Envelope
events.Op4.Envelope = c(I = 3,  DR = 30, END = lifetime)
duration.ev.Envelope = c(I = 1, DR = 0.5)

design.Options.Elevator <- list(desing.Op1.Elevator <- "elevator_MRL")
design.Options.Sprinklerpipe<- list(desing.Op2.Sprinklerpipe <- "sprinklerpipe_PVC")
design.Options.ETICS<- list(desing.Op3.ETICS <- "ETICS_EP")
design.Options.Envelope<- list(desing.Op4.Envelope <- "Envelope_Glulam")


maintenance.Elevator <- dist.Events(lifetime, events.Op1.Elevator, start.Date, design.Options.Elevator)
maintenance.Sprinklerpipe <- dist.Events(lifetime, events.Op2.Sprinklerpipe, start.Date, design.Options.Sprinklerpipe)
maintenance.ETICS <- dist.Events(lifetime, events.Op3.ETICS, start.Date, design.Options.ETICS)
maintenance.Envelope <- dist.Events(lifetime, events.Op4.Envelope, start.Date, design.Options.Envelope)

### for plotting
#maintenance.All <- dist.Events(lifetime, c(events.Op1.Elevator, events.Op2.Sprinklerpipe, events.Op3.ETICS, events.Op4.Envelope), start.Date, "Pedestrian walkway")


integrated.interv <- combine.lifeTimelines(maintenance.Elevator, maintenance.Sprinklerpipe,maintenance.ETICS,maintenance.Envelope,
                                           duration.ev.Elevator, duration.ev.Sprinklerpipe, duration.ev.ETICS,duration.ev.Envelope)

sum(integrated.interv$duration)

print(integrated.interv$duration)
events.Op1.Elevator = c(RD = 10, RM = 15, RS = 10, RDS = 24, END = lifetime)
duration.ev.Elevator <- c(RD = 8.8, RM = 5.7, RS = 0.1, RDS = 0.2)
events.Op1.Sprinklerpipe<- c(CF=5,PR=15,END=lifetime)
duration.ev.Sprinklerpipe <- c(CF=1.75,PR=3)
events.Op3.ETICS<-c(C = 5, FC = 15, END = lifetime)
duration.ev.ETICS<-c(C = 1, FC = 1)
events.Op4.Envelope = c(I = 5,  DR = 30, END = lifetime)
duration.ev.Envelope = c(I = 1, DR = 0.5)

### for plotting
maintenance.All <- dist.Events(lifetime, c(events.Op1.Elevator, events.Op2.Sprinklerpipe, events.Op3.ETICS, events.Op4.Envelope), start.Date, "Pedestrian walkway_combined maintenance")


design.Options.Elevator <- list(desing.Op1.Elevator <- "elevator_MRL")
design.Options.Sprinklerpipe<- list(desing.Op2.Sprinklerpipe <- "sprinklerpipe_PVC")
design.Options.ETICS<- list(desing.Op3.ETICS <- "ETICS_EP")
design.Options.Envelope<- list(desing.Op4.Envelope <- "Envelope_Glulam")
maintenance.Elevator <- dist.Events(lifetime, events.Op1.Elevator, start.Date, design.Options.Elevator)
maintenance.Sprinklerpipe <- dist.Events(lifetime, events.Op2.Sprinklerpipe, start.Date, design.Options.Sprinklerpipe)
maintenance.ETICS <- dist.Events(lifetime, events.Op3.ETICS, start.Date, design.Options.ETICS)
maintenance.Envelope <- dist.Events(lifetime, events.Op4.Envelope, start.Date, design.Options.Envelope)
integrated.interv <- combine.lifeTimelines(maintenance.Elevator, maintenance.Sprinklerpipe, maintenance.ETICS, maintenance.Envelope,
                                           duration.ev.Elevator, duration.ev.Sprinklerpipe, duration.ev.ETICS,duration.ev.Envelope)
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, design.Options.Elevator)
    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.Sprinklerpipe)
      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.ETICS)
        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.Envelope)
          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)]))
          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.grid.Elevator <- expand.grid(RD = sample(seq(8,20, by = 1), n.grid),
                                  RM = sample(seq(10,20,1), n.grid),
                                  RS = sample(seq(2, 16, 1), n.grid),
                                  RDS = sample(seq(20, 30, 1), n.grid))


events.grid.Sprinklerpipe<- expand.grid(CF = sample(seq(1, 10, by = 0.5), n.grid),
                                PR= sample(seq(7,23,1), n.grid))



events.grid.ETICS<- expand.grid(C = sample(seq(1, 5, by = 0.2), n.grid),
                                        FC= sample(seq(8,22,1), n.grid))


events.grid.Envelope<- expand.grid(I = sample(seq(1, 5, by = 0.1), n.grid),
                                DR= sample(seq(15,40,1), n.grid))


response.space <- design.explore(events.grid.Elevator, events.grid.Sprinklerpipe, events.grid.ETICS, events.grid.Envelope)
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)




#-------------------------------------- LCA ---------------------------------------------#




######################################

LCImaterials <- read.csv("Group3_LCImaterials.csv", sep = ",", header = T)

LCImaterials$material
######################################ETICS###############################

LCA.ETICS <- function(length, width, materials, ETICSarea, ETICSinterventions) 
{
  ETICS.area <- length * 3.14 * width
  EPSfactor<- 3.5
  baseadweight<- 5 
  rbase.factor <- 4
  meshfactor <- 1
  renderfactor <- 4
  primerfactor <- 0.25
  fcfactor <- 2.90
  #Interventions
  re.intervention <- 5
  pri.intervention <-5
  fc.intervention <-5
  
  
  #ENERGY AND EMISSIONS CALCULATION
  Etics.Energy <- LCImaterials$energy[1]*ETICS.area*baseadweight+ 
    LCImaterials$energy[2]*ETICS.area*EPSfactor+
    LCImaterials$energy[3]*ETICS.area*rbase.factor+
    LCImaterials$energy[4]*ETICS.area*meshfactor+
    LCImaterials$energy[5]*ETICS.area*renderfactor*re.intervention+
    LCImaterials$energy[6]*ETICS.area*primerfactor*pri.intervention+
    LCImaterials$energy[7]*ETICS.area*fcfactor*fc.intervention
  
  Etics.CO2 <- LCImaterials$CO2[1]*ETICS.area*baseadweight+ 
    LCImaterials$CO2[2]*ETICS.area*EPSfactor+
    LCImaterials$CO2[3]*ETICS.area*rbase.factor+
    LCImaterials$CO2[4]*ETICS.area*meshfactor+
    LCImaterials$CO2[5]*ETICS.area*renderfactor*re.intervention+
    LCImaterials$CO2[6]*ETICS.area*primerfactor*pri.intervention+
    LCImaterials$CO2[7]*ETICS.area*fcfactor*fc.intervention
  
  Etics.SO2 <- LCImaterials$SO2[1]*ETICS.area*baseadweight+ 
    LCImaterials$SO2[2]*ETICS.area*EPSfactor+
    LCImaterials$SO2[3]*ETICS.area*rbase.factor+
    LCImaterials$SO2[4]*ETICS.area*meshfactor+
    LCImaterials$SO2[5]*ETICS.area*renderfactor*re.intervention+
    LCImaterials$SO2[6]*ETICS.area*primerfactor*pri.intervention+
    LCImaterials$SO2[7]*ETICS.area*fcfactor*fc.intervention
  
  #LIST WITH ADDED ENERGY AND EMISSIONS
  ETICS.LCA.results <- list(Energy = (Etics.Energy) ,   
                            CO2 = (Etics.CO2),         
                            SO2 = (Etics.SO2))
  
  
  return(ETICS.LCA.results)
} 

######################################SPRINKLER PIPE###############################
LCA.Sprinklerpipe <- function(length, width, materials, pipelength, sprinklerpipeinterventions)
{  
  
  pipelength <-((length*width)/5.76)*2.4
  sprinklerpipeinterventions<-2
    
    
    
  #ENERGY AND EMISSIONS CALCULATION
  Sprinklerpipe.Energy <- LCImaterials$energy[9]*pipelength*sprinklerpipeinterventions+ 
    LCImaterials$energy[10]*pipelength*sprinklerpipeinterventions+
    LCImaterials$energy[11]*pipelength*sprinklerpipeinterventions
    
  
  Sprinklerpipe.CO2 <- LCImaterials$CO2[9]*pipelength*sprinklerpipeinterventions+ 
    LCImaterials$CO2[10]*pipelength*sprinklerpipeinterventions+
    LCImaterials$CO2[11]*pipelength*sprinklerpipeinterventions
  
  Sprinklerpipe.SO2 <- LCImaterials$SO2[9]*pipelength*sprinklerpipeinterventions+ 
    LCImaterials$SO2[10]*pipelength*sprinklerpipeinterventions+
    LCImaterials$SO2[11]*pipelength*sprinklerpipeinterventions
  
  
  #LIST WITH ADDED ENERGY AND EMISSIONS
  Sprinklerpipe.LCA.results <- list(Energy = (Sprinklerpipe.Energy) ,   
                            CO2 = (Sprinklerpipe.CO2),         
                            SO2 = (Sprinklerpipe.SO2))
  
  
  return(Sprinklerpipe.LCA.results)
}


######################################ELEVATOR###############################

LCA.Elevator <- function(numberofelevators,materials)
  
{
  numberofelevators <- 2
  
  Elevator.Energy <- numberofelevators*(LCImaterials$energy[12] + LCImaterials$energy[13])
  Elevator.CO2 <- numberofelevators*(LCImaterials$CO2[12] + LCImaterials$energy[13])
  Elevator.SO2 <- numberofelevators*(LCImaterials$SO2[12] + LCImaterials$energy[13])
  
  Elevator.LCA.results <- list(Energy = (Elevator.Energy) ,   
                                    CO2 = (Elevator.CO2),         
                                    SO2 = (Elevator.SO2))
  return(Elevator.LCA.results)
}

######################################CLIMATE###############################

#ENERGY AND EMISSIONS CALCULATION

LCA.Envelope <- function(numberofgirders,Envelopeinterventions,materials)
  
{
  

  
  numberofgirders <- 8
  Envelopeinterventions<-2
  
  
  
  Envelope.Energy <- numberofgirders*(LCImaterials$energy[14]*Envelopeinterventions)
  Envelope.CO2 <- numberofgirders*(LCImaterials$CO2[14]*Envelopeinterventions)
  Envelope.SO2 <- numberofgirders*(LCImaterials$SO2[14]*Envelopeinterventions)
  
  Envelope.LCA.results <- list(Energy = (Envelope.Energy) ,   
                               CO2 = (Envelope.CO2),         
                               SO2 = (Envelope.SO2))
  return(Envelope.LCA.results)
}
    
  
#Walkway Dimensions
w.length <- 28 # units: m
w.width <- 15 #units: m
w.height <- 21 #units: m


Option.Elevator<- LCA.Elevator (numberofelevators,LCImaterials)
Option.Sprinklerpipe <- LCA.Sprinklerpipe(w.length, w.width, LCImaterials, pipelength, sprinklerpipeinterventions)
Option.ETICS <- LCA.ETICS(w.length, w.width, LCImaterials, ETICSarea, ETICSinterventions) 
Option.Envelope <- LCA.Envelope(numberofgirders,Envelopeinterventions,LCImaterials)

integrated.Design <- as.data.frame(list(Energy = Option.Elevator$Energy + Option.Sprinklerpipe$Energy+Option.ETICS$Energy+Option.Envelope$Energy,
                                        CO2 = Option.Elevator$CO2 + Option.Sprinklerpipe$CO2+Option.ETICS$CO2+Option.Envelope$CO2,
                                        SO2 = Option.Elevator$SO2 + Option.Sprinklerpipe$SO2+Option.ETICS$SO2+Option.Envelope$SO2))
Energy.costs <- 1.08#0.128
CO2.unitcost <- 0.09#26 # per metric tonne
SO2.unitCosts <- 15.1#85 # per metric tonne

integrated.Design <- mutate(integrated.Design, Costs = (Energy * Energy.costs + CO2*CO2.unitcost +
                                                         SO2 * SO2.unitCosts))
integrated.Design

pieplotdf <- data.frame(system = c("Elevator", "Sprinklerpipe", "ETICS", "Envelope"), 
           energy_MJ = c(Option.Elevator$Energy, Option.Sprinklerpipe$Energy, Option.ETICS$Energy, Option.Envelope$Energy),
           CO2_kg = c(Option.Elevator$CO2, Option.Sprinklerpipe$CO2, Option.ETICS$CO2, Option.Envelope$CO2),
           SO2_kg = c(Option.Elevator$SO2, Option.Sprinklerpipe$SO2, Option.ETICS$SO2, Option.Envelope$SO2), 
           cost_EUR = c(Option.Elevator$Energy*Energy.costs+Option.Elevator$CO2*CO2.unitcost+Option.Elevator$SO2*SO2.unitCosts,
                    Option.Sprinklerpipe$Energy*Energy.costs+Option.Sprinklerpipe$CO2*CO2.unitcost+Option.Sprinklerpipe$SO2*SO2.unitCosts,
                    Option.ETICS$Energy*Energy.costs+Option.ETICS$CO2*CO2.unitcost+Option.ETICS$SO2*SO2.unitCosts,
                    Option.Envelope$Energy*Energy.costs+Option.Envelope$CO2*CO2.unitcost+Option.Envelope$SO2*SO2.unitCosts)
           )

ggplot(pieplotdf, aes(x="", y=energy_MJ, fill=system)) +
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0) +
  labs(x = NULL, y = NULL, fill = NULL, title = "energy") +
  theme_void()

ggplot(pieplotdf, aes(x="", y=CO2_kg, fill=system)) +
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0) +
  labs(x = NULL, y = NULL, fill = NULL, title = "CO2") +
  theme_void()

ggplot(pieplotdf, aes(x="", y=SO2_kg, fill=system)) +
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0) +
  labs(x = NULL, y = NULL, fill = NULL, title = "SO2") +
  theme_void()

ggplot(pieplotdf, aes(x="", y=cost_EUR, fill=system)) +
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0) +
  labs(x = NULL, y = NULL, fill = NULL, title = "cost") +
  theme_void()


#-------------------------------------- MDO ---------------------------------------------#
##################################### Fitness Function #########################################

fitness <- function(x) {
  #define the output dimension
  z <- numeric(6)
  #function body
  
  x <- round(x, 0)
  y <- expand.grid(RD = x[1], RM = x[2], RS = x[3], RDS = x[4],CF = x[5], PR = x[6],C=x[7], FC=x[8], I=x[9],DR=x[10])
  
  dur.ev1 <- unlist(y[1:4] / 4)
  dur.ev2 <- unlist(y[5:6] / 4)
  dur.ev3 <- unlist(y[7:8] / 4)
  dur.ev4 <- unlist(y[9:10] / 4)
  
  dist.1 <- apply(y[1:4], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.Elevator)
  dist.2 <- apply(y[5:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.Sprinklerpipe)
  dist.3 <- apply(y[7:8], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.ETICS)
  dist.4 <- apply(y[9:10], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date,design.Options.Envelope)
  
  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)]))
  
  

  w.length <- x[11] # units: m
  w.width <- x[12] #units: m
  w.height <- 21 #units: m
  materials <- LCImaterials
  Energy.costs <- 0.128
  CO2.unitcost <- 26
  SO2.unitCosts <- 85
  
  product.1.Elevator <- LCA.Elevator(numberofelevators, LCImaterials)
  
  product.2.Sprinklerpipe <- LCA.Sprinklerpipe(w.length, w.width, LCImaterials, dist.2[[1]])
  
  product.3.ETICS <- LCA.ETICS(w.length, w.width, LCImaterials, dist.3[[1]])
  
  product.4.Envelope <- LCA.Envelope(numberofgirders,dist.4[[1]],LCImaterials)
  
  
  integrated.system <- as.data.frame(list(Energy = product.1.Elevator$Energy + product.2.Sprinklerpipe$Energy+product.3.ETICS$Energy + product.4.Envelope$Energy,
                                          CO2 = product.1.Elevator$CO2 + product.2.Sprinklerpipe$CO2 + product.3.ETICS$CO2+ product.4.Envelope$CO2,
                                          SO2 = product.1.Elevator$SO2 + product.2.Sprinklerpipe$SO2 + product.3.ETICS$CO2+ product.4.Envelope$SO2))
  integrated.system <- mutate(integrated.system,
                              Costs = (Energy * Energy.costs +
                                         CO2*CO2.unitcost +
                                         SO2 * SO2.unitCosts))
  
  z[3] <- integrated.system$Energy
  z[4] <- integrated.system$CO2
  z[5] <- integrated.system$SO2
  z[6] <- integrated.system$Costs
  
  
  return(z)
}


r2 <- nsga2(fitness, idim = 12, odim = 6,
            generations=10, popsize=100,
            lower.bounds=c(8,10,2,20,1,7,1,8,1,15,20,10),
            upper.bounds= c(20,20,16,30,10,23,5,22,5,40,35,20))

r2$value
r2$par
r2$pareto.optimal
which(r2$pareto.optimal==TRUE)

#r2Results[1:4,]


r2Results <- as.data.frame(r2$value)
outNames <- c("duration", "interv.dist", "energy", "co2", "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("RD", "RM", "RS", "RDS","CF", "PR", "C", "FC","I","DR","w.length","w.width",
                           "duration", "interv.dist","energy", "co2", "so2", "cost","pareto")
r2$value
r2$par
r2$pareto.optimal
which(r2$pareto.optimal==TRUE)
#all.results[1:4,]


par(mfrow = c(1,1))
parcoord(all.results[, 1:19], 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))