##################################################################
#GROUP C - LCA WS21/22
#ANAS AL ABAJI
#PRATHYASH KURIAN
#KYLE SONGPANYA
#KARTHEEK MYLAVARAPU
####################################################################


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)] <- "Construction"
  #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] <- "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])]
    
    #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 = "2022-01-01"
lifetime <- 30 # Since most of the systems had a lifetime of around 25

#Events and Duration of each specific intervention for every subsystem
#1. Subsystem: Offshore Piles
events.Op1.OffshorePiles = c(ReMa = 5, PiRe = 10 , END = lifetime) # ReMa = regular maintenance; PiRe = pile replacement
duration.ev.OffshorePiles <- c(ReMa = 5, PiRe = 10)# in days

#2. Subsystem: Mooring Lines
events.Op2.MooringLines = c(VI = 5,CR = 25 , END = lifetime) # VI = visual inspection;CR = complete replacement 
duration.ev.MooringLines <- c(VI = 5,CR = 10 ) # in days

#3. Subsystem: TT
events.Op3.TT = c(GI = 5, HI = 6, AI = 8, AC = 10, END = lifetime) # GI = ground inspection; HI = heat inspection; AI = air inspection; AC = Air control;
duration.ev.TT <- c(GI = 2, HI = 2, AI = 15, AC = 22)# in days

#4. Subsystem: Wind Turbine
events.Op4.WindTurbine = c(RM = 5, SSR = 4, PR = 8,  END = lifetime) # RM = regular maintenance; SSR = Support Structure Repair; PR = pile repair 
duration.ev.WindTurbine <- c(RM = 3, SSR = 2, PR = 8)# in days

#Define design options for each system
design.Options.OffshorePiles <- list(design.Op1.OffshorePiles <- "1. Offshore Piles")
design.Options.MooringLines <- list(design.Op2.MooringLines <- "2.Mooring Lines")
design.Options.TT <- list(design.Op3.TT <- "3. Transmission Tower")
design.Options.WindTurbine <- list(design.Op4.WindTurbine <- "4. Wind Turbine")

#Generate distribution of interventions over lifetime
maintenance.OffshorePiles <- dist.Events(lifetime, events.Op1.OffshorePiles, start.Date, design.Options.OffshorePiles)
maintenance.MooringLines <- dist.Events(lifetime, events.Op2.MooringLines, start.Date, design.Options.MooringLines)
maintenance.TT <- dist.Events(lifetime, events.Op3.TT, start.Date, design.Options.TT)
maintenance.WindTurbine <- dist.Events(lifetime, events.Op4.WindTurbine, start.Date, design.Options.WindTurbine)

#Combination of all four timelines
integrated.interv <- combine.lifeTimelines(maintenance.OffshorePiles, maintenance.MooringLines, maintenance.TT, maintenance.WindTurbine,
                                           duration.ev.OffshorePiles, duration.ev.MooringLines, duration.ev.TT, duration.ev.WindTurbine)

#Sum up total duration of interruptions for the set lifetime of 25 years
sum(integrated.interv$duration)

#########################################################################
# #Hypothetical Maintenance Scenario
# #1. Subsystem: Offshore Piles
# events.Op1.OffshorePiles = c(RM = 5, PR = 2 , END = lifetime) # RM = regular maintenance; PR = pile replacement
# duration.ev.OffshorePiles <- c(RM = 3, PR = 8)# in days
# 
# #2. Subsystem: Mooring Lines
# events.Op2.MooringLines = c(VI = 5,CR = 1 , END = lifetime) # VI = visual inspection;CR = complete replacement 
# duration.ev.MooringLines <- c(VI = 5,CR = 8 ) # in days
# 
# #3. Subsystem: TT
# events.Op3.TT = c(GI = 1, HI = 2, AI = 2, AC = 3, END = lifetime) # GI = ground inspection; HI = heat inspection; AI = air inspection; AC = Air control;
# duration.ev.TT <- c(GI = 2, HI = 2, AI = 15, AC = 22)# in days
# 
# #4. Subsystem: Wind Turbine
# events.Op4.WindTurbine = c(RM = 5, SSR = 6, PR = 3,  END = lifetime) # RM = regular maintenance; SSR = Support Structure Report; PR = pile repair 
# duration.ev.WindTurbine <- c(RM = 3, SSR = 2, PR = 8)# in days
# 
# #Define design options for each system
# design.Options.OffshorePiles <- list(design.Op1.OffshorePiles <- "1. Offshore Piles")
# design.Options.MooringLines <- list(design.Op2.MooringLines <- "2.Mooring Lines")
# design.Options.TT <- list(design.Op3.TT <- "3. Transmission Tower")
# design.Options.WindTurbine <- list(design.Op4.WindTurbine <- "4. Wind Turbine")
# 
# #Generate distribution of interventions over lifetime
# maintenance.OffshorePiles <- dist.Events(lifetime, events.Op1.OffshorePiles, start.Date, design.Options.OffshorePiles)
# maintenance.MooringLines <- dist.Events(lifetime, events.Op2.MooringLines, start.Date, design.Options.MooringLines)
# maintenance.TT <- dist.Events(lifetime, events.Op3.TT, start.Date, design.Options.TT)
# maintenance.WindTurbine <- dist.Events(lifetime, events.Op4.WindTurbine, start.Date, design.Options.WindTurbine)
# 
# #Combination of all four timelines
# integrated.interv <- combine.lifeTimelines(maintenance.OffshorePiles, maintenance.MooringLines, maintenance.TT, maintenance.WindTurbine,
#                                            duration.ev.OffshorePiles, duration.ev.MooringLines, duration.ev.TT, duration.ev.WindTurbine)
# 
# #Sum up total duration of interruptions for the set lifetime of 25 years
# sum(integrated.interv$duration)
#########################################################################
#Function to automate the combination of 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.OffshorePiles)
    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.MooringLines)
      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.TT)
        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.WindTurbine)
          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 # this can be varied however, taking larger numbers leads to increase in computational time and we don't have time lol.
#ReMa = 5, PiRe = 10
#VI = 5,CR = 25
#GI = 5, HI = 6, AI = 8, AC = 10
#RM = 5, SSR = 4, PR = 8

#Intervention matrix scenario for the OffshorePiles
# Here the stuff inside the seq part is the range of values we are considering for each type of maintenance 
events.grid.OffshorePiles <- expand.grid(ReMa = sample(seq(2,10, by = 1), n.grid),
                                 PiRe = sample(seq(5,15, by = 1), n.grid))
events.grid.OffshorePiles

#Intervention matrix scenario for the MooringLines
events.grid.MooringLines <- expand.grid(VI = sample(seq(3,10, by = 1), n.grid),
                                        CR = sample(seq(15,29, by = 1), n.grid))
events.grid.MooringLines

#Intervention matrix scenario for the TT
events.grid.TT <- expand.grid(GI = sample(seq(4,12, by = 1), n.grid),
                              HI = sample(seq(3,15, by = 1), n.grid),
                              AI = sample(seq(5,14, by = 1), n.grid),
                              AC = sample(seq(4,20, by = 1), n.grid))
events.grid.TT

#Intervention matrix scenario for the Wind Turbine
events.grid.WindTurbine <- expand.grid(RM = sample(seq(3,12, by = 1), n.grid),
                                       SSR = sample(seq(2,16, by = 1), n.grid),
                                       PR = sample(seq(6,12, by = 1), n.grid))
events.grid.WindTurbine

#Use the function design.explore() for automation of system combination
response.space <- design.explore(events.grid.OffshorePiles, events.grid.MooringLines, events.grid.TT, events.grid.WindTurbine)
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. What are our preferences? Our preference is to minimize the duration of the interruptions and to maximize the distance between interventions
p2 <- high(dur) * low(dist.inter)
show_front(p2)
sky2 <- psel(response.space, p2)


#############  2. LCI and Analysis of Integrated Civil Systems ##############
# change this to the directory where the csv file is on your PC

LCI.materials <- read.csv("groupC.csv")

#ReMa = 5, PiRe = 10
#ViIn = 5,CoRe = 25
#GI = 5, HI = 6, AiIn = 8, AC = 10
#RM = 5, SSR = 4, PR = 8
##############################  Offshore Piles (Anas)  ############################## 

LCA.OffshorePiles <- function(dist.event) { 

  
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  
  PR.interventions <- (if(!is.null(s2$PiRe)){s2$PiRe} else {0})

  OffshorePiles.Energy <- LCI.materials$Energy[12]*(1 + PR.interventions) # Energy is needed for the first construction and for the number of Pile replacements. Hence Energy*(1 + Number of replacements)
  OffshorePiles.CO2 <- LCI.materials$CO2[12]*(1 + PR.interventions)
  OffshorePiles.NOx <- LCI.materials$NOx[12]*(1 + PR.interventions)
  OffshorePiles.SO2 <- LCI.materials$SO2[12]*(1 + PR.interventions)
  
  OffshorePiles.LCA.results <- list(Energy = (OffshorePiles.Energy) ,   
                                    CO2 = (OffshorePiles.CO2),         
                                    NOx = (OffshorePiles.NOx),
                                    SO2 = (OffshorePiles.SO2))
  
  return(OffshorePiles.LCA.results)
  
}


##############################  MooringLines (Kyle)  ##############################

LCA.MooringLines <- function(materials,dist.event) { #
  
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  
  CR.intervals <- (if(!is.null(s2$CR)){s2$CR} else {0})
  materials.split <- split(materials, materials$Scope)
  MooringLines.Energy <- sum(materials.split$Moor[3])*(1 + CR.intervals)
  MooringLines.CO2 <- sum(materials.split$Moor[5])*(1 + CR.intervals)
  MooringLines.NOx <- sum(materials.split$Moor[6])*(1 + CR.intervals)
  MooringLines.SO2 <- sum(materials.split$Moor[7])*(1 + CR.intervals)
  
  MooringLines.LCA.results <- list(Energy = (MooringLines.Energy) ,
                                  CO2 = (MooringLines.CO2),
                                  NOx = (MooringLines.NOx),
                                  SO2 = (MooringLines.SO2))
  
  return(MooringLines.LCA.results)
  
}

##############################  TT (Kartheek)  ##############################

LCA.TT <- function(materials,dist.event) {
  
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  AI.intervals <- (if(!is.null(s2$AI)){s2$AI} else {0})
  AC.intervals <- (if(!is.null(s2$AC)){s2$AC} else {0})
  materials.split <- split(materials, materials$Scope)
  TT.Energy <- sum(materials.split$SLT[3])*(1 + (AI.intervals)*0.3 + (AC.intervals)*0.2)
  TT.CO2 <- sum(materials.split$SLT[5])*(1 + (AI.intervals)*0.3 + (AC.intervals)*0.2)
  TT.NOx <- sum(materials.split$SLT[6])*(1 + (AI.intervals)*0.3 + (AC.intervals)*0.2)
  TT.SO2 <- sum(materials.split$SLT[7])*(1 + (AI.intervals)*0.3 + (AC.intervals)*0.2)
  
  TT.LCA.results <- list(Energy = (TT.Energy) ,
                         CO2 = (TT.CO2),
                         NOx = (TT.NOx),
                         SO2 = (TT.SO2))
  
  return(TT.LCA.results)
  
}


##############################  WindTurbine (Prathyash)  ##############################

LCA.WindTurbine<-function(materials,dist.event){
  
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  PR.intervals <- (if(!is.null(s2$PR)){s2$PR} else {0})
  materials.split <- split(materials, materials$Scope)
  WindTurbine.Energy <- sum(materials.split$MonopileConcrete[3])*(1 + PR.intervals)
  WindTurbine.CO2 <- sum(materials.split$MonopileConcrete[5])*(1 + PR.intervals)
  WindTurbine.NOx <- sum(materials.split$MonopileConcrete[6])*(1 + PR.intervals)
  WindTurbine.SO2 <- sum(materials.split$MonopileConcrete[7])*(1 + PR.intervals)
  
  WindTurbine.LCA.results <- list(Energy = (WindTurbine.Energy) ,
                         CO2 = (WindTurbine.CO2),
                         NOx = (WindTurbine.NOx),
                         SO2 = (WindTurbine.SO2))
  return(WindTurbine.LCA.results)
  
}



#Application of LCA-functions to calculate the impact on environment
Option.OffshorePiles <- LCA.OffshorePiles(maintenance.OffshorePiles)
Option.MooringLines <- LCA.MooringLines(LCI.materials,maintenance.MooringLines)
Option.TT <- LCA.TT(LCI.materials,maintenance.TT)
Option.WindTurbine <- LCA.WindTurbine(LCI.materials,maintenance.WindTurbine)

#Combination of results of all subsystems into one dataframe
integrated.Design <- as.data.frame(list(Energy = Option.OffshorePiles$Energy + Option.MooringLines$Energy + Option.TT$Energy + Option.WindTurbine$Energy,
                                        CO2 = Option.OffshorePiles$CO2 + Option.MooringLines$CO2 + Option.TT$CO2 + Option.WindTurbine$CO2,
                                        NOx = Option.OffshorePiles$NOx + Option.MooringLines$NOx + Option.TT$NOx + Option.WindTurbine$NOx,
                                        SO2 = Option.OffshorePiles$SO2 + Option.MooringLines$SO2 + Option.TT$SO2 + Option.WindTurbine$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 #########################################

#ReMa = 5, PiRe = 10
#VI = 5,CR = 25
#GI = 5, HI = 6, AI = 8, AC = 10
#RM = 5, SSR = 4, PR = 8

fitness <- function(x) {
  #Definition of output dimension
  z <- numeric(7)
  x <- round(x, 0)
  y <- expand.grid(ReMa = x[1], PiRe = x[2], VI = x[3],CR = x[4],GI = x[5],HI = x[6], AI = x[7], AC = x[8],RM = x[9], SSR = x[10], PR = x[11])
  
  dur.ev1 <- unlist(y[1:2] / 4)
  dur.ev2 <- unlist(y[3:4] / 4)
  dur.ev3 <- unlist(y[5:8] / 4)
  dur.ev4 <- unlist(y[9:11] / 4)
  
  dist.1 <- apply(y[1:2], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  dist.2 <- apply(y[3:4], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  dist.3 <- apply(y[5: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)]))
  
  #Set static variables
  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.OffshorePiles <- LCA.OffshorePiles(dist.1[[1]])
  
  product.2.MooringLines <- LCA.MooringLines(LCI.materials,dist.2[[1]])
  
  product.3.TT <- LCA.TT(LCI.materials,dist.3[[1]])
  
  product.4.WindTurbine <- LCA.WindTurbine(LCI.materials,dist.4[[1]])
  
  #Call function for design integration
  integrated.system <- as.data.frame(list(Energy = product.1.OffshorePiles$Energy + product.2.MooringLines$Energy + product.3.TT$Energy + product.4.WindTurbine$Energy,
                                          CO2 = product.1.OffshorePiles$CO2 + product.2.MooringLines$CO2 + product.3.TT$CO2 + product.4.WindTurbine$CO2,
                                          NOx = product.1.OffshorePiles$NOx + product.2.MooringLines$NOx + product.3.TT$NOx + product.4.WindTurbine$NOx,
                                          SO2 = product.1.OffshorePiles$SO2 + product.2.MooringLines$SO2 + product.3.TT$SO2 + product.4.WindTurbine$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)
}
#ReMa = x[1], PiRe = x[2], VI = x[3],CR = x[4],GI = x[5],HI = x[6], AI = x[7], AC = x[8],RM = x[9], SSR = x[10], PR = x[11]
#Define GA optimization function
r2 <- nsga2(fitness, idim = 11, odim = 7,
            generations=10, popsize=100,
            lower.bounds=c(2, 5, 3,15,4,3,5,4,3,2,6),
            upper.bounds= c(10,15,10,29,12,15,14,20,12,16,12))

r2$value
r2$par
r2$pareto.optimal
which(r2$pareto.optimal==TRUE)
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")


#How does accumulated impact of input parameters affect the performance criteria?
#ReMa = x[1], PiRe = x[2], VI = x[3],CR = x[4],GI = x[5],HI = x[6], AI = x[7], AC = x[8],RM = x[9], SSR = x[10], PR = x[11]
input.params <- round(r2$par, 0)
all.results <- cbind(input.params, r2Results, r2$pareto.optimal)
colnames(all.results) <- c("ReMa", "PiRe", "VI", "CR", "GI", "HI", "AI", "AC","RM", "SSR", "PR",
                           "duration", "interv.dist","energy", "co2", "nox", "so2", "cost", "pareto")
#interpretation
#setting the number of alternatives we have
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))