setwd("C:\\Users\\user\\Desktop\\LCA_Projects\\LCA_Group")
#import libraries
library(timelineS) 
library(lubridate)
library(ggplot2) 
library(dplyr)
library(reshape2) 
library(rPref)
library(mco) 
library(MASS)

rm(list = ls()) #clean the environment


#Integrated Maintenance Planning
#First Step: Define Functions
#create plot function
plot.timeline <- function(lifetime, events.name, start.Date, plot.Name, event.colors){  
  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,
            label.color = event.colors) 
}

#create dist.event function
dist.Events <- function (lifetime, events, start.Date, option.Name, event.colors) {
  #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"
  
  # Generate color vector based on event names
  event.color <- event.colors[match(unique.events$event, names(event.colors))]
  
  #plot the timeline
  plot.timeline(unique.events$frequency, unique.events$event, start.Date, option.Name, event.color)
  return(unique.events)
}

# Define color vector for each event
event.colors.Buildings <- c(DC = "black", H.B = "blue", M.B = "green", F.B = "red", S.B = "purple", END = "black")
event.colors.Water <- c(DC = "black", SI.W = "lightblue", P.W = "lightgreen", W.W = "yellowgreen", T.W = "coral", END = "black")

#create function combine the life cycle timelines of individual systems
combine.lifeTimelines <- function(product.1, product.2, p1.dur, p2.dur) {
  #define vectors contain the results
  base <- list(Names = c(), frequency = c(), duration = c())
  #take the lifetime distribution of subsystems and combine them in a list with unique values
  base$frequency <- sort(unique(c(product.1$frequency,
                                  product.2$frequency)), decreasing = FALSE)
  
  #assume subsystems have same starting data & lifetime
  #initilaize and name the beginning & the end of the combined sequence ,
  #for which the duration of the intervention is 0, meaning the constrction time won't be taked into account
  #as well ass the time need to demolish the subsystems
  base$Names[1] <- "DC" #beginning of the intervention
  base$duration[1] <- 0
  base$Names[length(base$frequency)] <- "END" 
  base$duration[length(base$frequency)] <- 0
  
  #iterate over the combined list,
  #and add right names for each phase based on values we have in the unique list of events
  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])] 
    
    if(length(phase.1) != 0 && length(phase.2) != 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)])
    } 
    else if(length(phase.1) != 0 && length(phase.2) == 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) { 
      base$Names[index] <- phase.2
      base$duration[index] <- p2.dur[which(names(p2.dur) == phase.2)]
    } 
    else {
      base$duration[index] <- 0
    } }
  return(base)
}


#Second Step: Define the information for life cycle timelines
#define the start date, the duration for the hypothetical lifetime for subsystems
#and the interventions and frequency(unit: year)
par(mfrow = c(2, 1)) 
start.Date = "2024-01-01" 
lifetime <- 60

#========================================================
#           Second Maintenance Strategy
#========================================================
#define information about duration for each intervention
#keep same names for simplicity & to automate some of the tasks
#define the duration of interventions in days
# Buildings
events.Op1.Buildings = c(H.B = 2, M.B = 5, S.B = 15, END = lifetime)
duration.ev.Buildings <- c(H.B = 3, M.B = 5, S.B = 20)

# Water Supply Facility (steel-framed)
events.Op1.WaterSupplyFacility <- c(SI.W = 10, P.W = 5, W.W = 7, END = lifetime)
duration.ev.WaterSupplyFacility <- c(SI.W = 5, P.W = 3, W.W = 5)

#add the name of the design options
#use dist.events() to generate the distribution of the interventions over the lifetime
design.Options.Buildings <- list(desing.Op1.Buildings <- "Buildings")
design.Options.WaterSupplyFacility <- list(desing.Op1.WaterSupplyFacility <- "WaterSupplyFacility")

maintenance.Buildings <- dist.Events(lifetime, events.Op1.Buildings, start.Date, design.Options.Buildings, event.colors.Buildings)
maintenance.WaterSupplyFacility <- dist.Events(lifetime, events.Op1.WaterSupplyFacility, start.Date, design.Options.WaterSupplyFacility, event.colors.Water)

#combine the life cycle timelines of individual systems
integrated.interv <- combine.lifeTimelines(maintenance.Buildings,  maintenance.WaterSupplyFacility,
                                           duration.ev.Buildings, duration.ev.WaterSupplyFacility)

#analyze the total duration of the interventions
cat("Total Duration (Second Maintenance Strategy):", sum(integrated.interv$duration), "\n")


#========================================================
#           Third Maintenance Strategy
#========================================================
#define information about duration for each intervention
#keep same names for simplicity & to automate some of the tasks
#define the duration of interventions in days
#Buildings
events.Op1.Buildings = c(H.B = 2, M.B = 4, S.B = 20, END = lifetime)
duration.ev.Buildings <- c(H.B = 3, M.B = 5, S.B = 20)

# Water Supply Facility (steel-framed)
events.Op1.WaterSupplyFacility <- c(SI.W = 20, P.W = 4, W.W = 6, END = lifetime)
duration.ev.WaterSupplyFacility <- c(SI.W = 5, P.W = 3, W.W = 5)

#add the name of the design options
#use dist.events() to generate the distribution of the interventions over the lifetime
design.Options.Buildings <- list(desing.Op1.Buildings <- "Buildings")
design.Options.WaterSupplyFacility <- list(desing.Op1.WaterSupplyFacility <- "WaterSupplyFacility")

maintenance.Buildings <- dist.Events(lifetime, events.Op1.Buildings, start.Date, design.Options.Buildings, event.colors.Buildings)
maintenance.WaterSupplyFacility <- dist.Events(lifetime, events.Op1.WaterSupplyFacility, start.Date, design.Options.WaterSupplyFacility, event.colors.Water)

#combine the life cycle timelines of individual systems
integrated.interv <- combine.lifeTimelines(maintenance.Buildings,  maintenance.WaterSupplyFacility,
                                           duration.ev.Buildings, duration.ev.WaterSupplyFacility)

#analyze the total duration of the interventions
cat("Total Duration (Third Maintenance Strategy):", sum(integrated.interv$duration), "\n")

#-----------------------------------------
#create function to explore more scenarios
#-----------------------------------------

design.explore <- function(events1, events2) {
  results <- c()
  for(i in 1:dim(events1)[1]) {
    ev1 <- unlist(events1[i, ])
    dist.1 <- dist.Events(lifetime, ev1, start.Date, design.Options.Buildings, event.colors.Buildings)
    dur.ev1 <- ev1/2
    for (j in 1:dim(events2)[1]) {
      ev2 <- unlist(events2[j, ])
      dist.2 <- dist.Events(lifetime, ev2, start.Date, design.Options.WaterSupplyFacility, event.colors.Water)
      dur.ev2 <- ev2/2

      combined.lifetime <- combine.lifeTimelines(dist.1, dist.2,
                                                 dur.ev1, dur.ev2)
      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, dur = sum(combined.lifetime$duration), dist.inter = min.dist.int))
    }
  }
  return(as.data.frame(results))
}
#-----------------------------------------

#create parameter allows us to control how the parameters will vary
n.grid <- 5
#generate the intervention matrix scenarios for the buildings
events.grid.Buildings <- expand.grid(H.B = sample(seq(2,10, by = 1), n.grid),
                                     M.B = sample(seq(5,15,1), n.grid),
                                     S.B = sample(seq(15,30,1), n.grid))

#generate the design space for each variable which is a sequence of numbers within the interval [min, max]
#use the function sample(x, n) to get a number of n random samples
#expand.grid() will take the values of each variable,
#and will combine these into a matrix, making sure that the all possible combinations are covered

events.grid.WaterSupplyFacility <- expand.grid(SI.W = sample(seq(15,30,1), n.grid),
                                               P.W = sample(seq(1,10,1), n.grid),
                                               W.W = sample(seq(6,15,1), n.grid))
response.space <- design.explore(events.grid.Buildings, events.grid.WaterSupplyFacility)


#visualize and select the maxima based on specific preferences
##preference is to minimize the duration of the interruptions
##and to maximize the distance between interventions

#define the preference
p <- low(dur)* high(dist.inter)
#evaluate the preference we define and to select that alternative which is the best for our preference
sky <- psel(response.space, p)

#get alternatives ranked
pareto2 <- psel(response.space, p, top = nrow(response.space))
#visualization
ggplot(response.space, aes(x = dur, y = dist.inter)) + geom_point(shape = 21) +
  geom_point(data = pareto2, size = 3, aes(color = factor(.level)))

#create function to find the Pareto frontier for different preferences
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) }

#visualize the Pareto frontier of out preference p
show_front(p)
#easily simulate different alternatives for preferences
p <- high(dur) * low(dist.inter)
#Show Pareto front
show_front(p)

print(paste("end:", Sys.time()))

#Section 5
#explore the idea of multi objective optimization using algorithms
#based on predefined design space will automatically look for optimal solutions

#Life Cycle Inventory and Analysis of Integrated Civil Systems
#the same life cycle inventory we create during previous assignment
LCI.materials <- read.csv("project_LCImaterials.csv")

#update the function to remove the options which we are no longer going to be used
#(setting different intervention senarios)
#and to consider the impact of the interventions based on their distribution over lifetime
LCA.building <- function(length, width, height, story, column_section, beam_section, structure.Option, materials, dist.event) {
  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))

  #Column
  if (structure.Option == "RC") {
    n <- round(width / 6, 0)*round(length / 6, 0) #every 6m establish a RC column
    interventions.column <- s2$DC + if (!is.null(s2$DC)) {s2$DC} else {0} + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    column.V <- n * column_section * height * story
    # print("Debug: Checking interventions.column before mutate()")
    # print(interventions.column)
  } else if (structure.Option == "Timber") {
    n <- round(width / 6, 0)*round(length / 6, 0) #every 6m establish a timber column
    interventions.column <- s2$DC + if (!is.null(s2$DC)) {s2$DC} else {0} + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    column.V <- n * column_section * height * story
    # print("Debug: Checking interventions.column before mutate()")
    # print(interventions.column)
  } else if (structure.Option == "Steel") {
    n <- round(width / 10, 0)*round(length / 10, 0) #every 10m establish a steel column
    interventions.column <- s2$DC + if (!is.null(s2$DC)) {s2$DC} else {0} + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    column.V <- n * column_section * height * story
    # print("Debug: Checking interventions.column before mutate()")
    # print(interventions.column)
  } else if (structure.Option == "Steel_Frame") {
    n <- round(width / 10, 0)*round(length / 10, 0) #every 10m establish a steel column
    interventions.column <- s2$DC + if (!is.null(s2$DC)) {s2$DC} else {0} + if (!is.null(s2$SI.W)) {s2$SI.W} else {0}
    column.V <- n * column_section * height * story
    # print("Debug: Checking interventions.column before mutate()")
    # print(interventions.column)
  }

  #Beam
  if (structure.Option == "RC") {
    n <- round(width / 6, 0)*round(length / 6, 0) #every 6m establish a RC beam
    interventions.beam <- s2$DC + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    beam.V <- n * height * story
  } else if (structure.Option == "Timber") {
    n <- round(width / 6, 0)*round(length / 6, 0) #every 6m establish a timber beam
    interventions.beam <- s2$DC + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    beam.V <- n * beam_section * height * story
  } else if (structure.Option == "Steel") {
    n <- round(width / 10, 0)*round(length / 10, 0) #every 10m establish a steel beam
    interventions.beam <- s2$DC + if (!is.null(s2$S.B)) {s2$S.B} else {0}
    beam.V <- n * beam_section * height * story
  } else if (structure.Option == "Steel_Frame") {
    n <- round(width / 10, 0)*round(length / 10, 0) #every 10m establish a steel beam
    interventions.beam <- s2$DC + if(! is.null(s2$SI.W)) {s2$SI.W} else {0}
    beam.V <- n * beam_section * height * story
  }

  #Slab
  if (structure.Option == "RC") {
    slab_thickness <- 0.15
    interventions.slab <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    slab.V <- length*width*slab_thickness * story
  } else if (structure.Option == "Timber") {
    slab_thickness <- 0.03
    interventions.slab <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    slab.V <- length*width*slab_thickness * story
  } else if (structure.Option == "Steel") {
    slab_thickness <- 0.01
    interventions.slab <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    slab.V <- length*width*slab_thickness * story
  } else if (structure.Option == "Steel_Frame") {
    slab_thickness <- 0.01
    interventions.slab <- s2$DC + if(! is.null(s2$SI.W)) {s2$SI.W} else {0} + 0.15 * if(! is.null(s2$P.W)) {s2$P.W} else {0}
    slab.V <- length*width*slab_thickness * story
  }

  #Facade
  if (structure.Option == "RC") {
    wall_thickness <- 0.15
    interventions.wall <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    wall.V <- length*width*wall_thickness * story
  } else if (structure.Option == "Timber") {
    wall_thickness <- 0.03
    interventions.wall <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    wall.V <- length*width*wall_thickness * story
  } else if (structure.Option == "Steel") {
    wall_thickness <- 0.01
    interventions.wall <- s2$DC + 0.15 * if (!is.null(s2$H.B)) {s2$H.B} else {0} + 0.15 * if(! is.null(s2$M.B)) {s2$M.B} else {0}
    wall.V <- length*width*wall_thickness * story
  } else if (structure.Option == "Steel_Frame") {
    wall_thickness <- 0.01
    interventions.wall <- s2$DC + if(! is.null(s2$SI.W)) {s2$SI.W} else {0} + 0.15 * if(! is.null(s2$P.W)) {s2$P.W} else {0}
    wall.V <- length*width*wall_thickness * story
  }

  #why there is a null?? is it because there could be no girder??
  if (!is.null(materials.split[[structure.Option]])) {
    column <- mutate(materials.split[[structure.Option]], building.Q = column.V, interventions = interventions.column)
    beam <- mutate(materials.split[[structure.Option]], building.Q = beam.V, interventions = interventions.beam)
    slab <- mutate(materials.split[[structure.Option]], building.Q = slab.V, interventions = interventions.slab)
    Wall <- mutate(materials.split[[structure.Option]], building.Q = wall.V, interventions = interventions.wall)
    LCA.matrix <- rbind(column, beam, slab, Wall)
  } else {
    LCA.matrix <- rbind(column, beam, slab, Wall)
  }

  LCA.matrix <- mutate(LCA.matrix,
                       TotalMaterials.Q = quantities * building.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)
}

#extract the frequency of different events repetition in s2 variable
#using these values to calculate the impact of the interventions over the lifecycle

# Define parameters for LCA
ob.length <- 24 # units: m
ob.width <- 30 #units: m
ob.height <- 4 #units: m
ob.story <- 30

sd.length <- 24 # units: m
sd.width <- 30 #units: m
sd.height <- 4 #units: m
sd.story <- 15

rb.length <- 30 # units: m
rb.width <- 60 #units: m
rb.height <- 4 #units: m
rb.story <- 3

wf.length <- 40 # units: m
wf.width <- 50 #units: m
wf.height <- 6 #units: m
wf.story <- 1

structures.Options <- c("RC", "Timber", "Steel", "Steel_Frame")

#call function to calculate the environmental impact
# Perform LCA for buildings
Option.OfficeBuilding <- LCA.building(ob.length, ob.width, ob.height, ob.story, 0.25, 0.12, structures.Options[1], LCI.materials, maintenance.Buildings)
Option.StudentDormitory <- LCA.building(sd.length, sd.width, sd.height, sd.story, 0.25, 0.12, structures.Options[2], LCI.materials, maintenance.Buildings)
Option.ResidentialBuilding <- LCA.building(rb.length, rb.width, rb.height, rb.story, 0.06, 0.03, structures.Options[3], LCI.materials, maintenance.Buildings)
Option.WaterSupplyFacility <- LCA.building(wf.length, wf.width, wf.height, wf.story, 0.04, 0.03, structures.Options[4], LCI.materials, maintenance.WaterSupplyFacility)


#combine the effect of both into one dataframe
# Combine LCA results
integrated.Design <- as.data.frame(list(Energy = Option.OfficeBuilding$Energy + Option.StudentDormitory$Energy
                                        + Option.ResidentialBuilding$Energy + Option.WaterSupplyFacility$Energy,
                                        CO2 = Option.OfficeBuilding$CO2 + Option.StudentDormitory$CO2
                                        + Option.ResidentialBuilding$CO2 + Option.WaterSupplyFacility$CO2,
                                        NOx = Option.OfficeBuilding$NOx + Option.StudentDormitory$NOx
                                        + Option.ResidentialBuilding$NOx + Option.WaterSupplyFacility$NOx,
                                        SO2 = Option.OfficeBuilding$SO2 + Option.StudentDormitory$SO2
                                        + Option.ResidentialBuilding$SO2 + Option.WaterSupplyFacility$SO2))

# Define unit costs
#CO2: carbon tax, others:costs for air purification
Energy.costs <- 1.08
CO2.unitcost <- 37
NOx.unitCost <- 9696
SO2.unitCosts <- 19773

# Calculate total costs
integrated.Design <- mutate(integrated.Design,
                            Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx * NOx.unitCost + SO2 * SO2.unitCosts)/10^9)


#Multi-Objective Optimization
#to implement a genetic algorithm for the multi criteria optimization
#define the function to be minimized, the number of input variables, the numbers of output measures
#as well as the lower and upper bounds for the input variables

#input variables: intervention duration in a specific range
#               ,building length, width, height, story
#               , as well thickness of the asphalt

#output variables: total duration of the interventions
#                 ,the minimum distance between two consecutive interventions
#                 ,energy, CO2, NOx, SO2, and Costs

#Based on aforementioned variables
#Define the function to be minimized
fitness <- function(x) {
  #define the output dimension
  z <- numeric(7)

  #distribute the interventions over the lifetime
  x <- round(x, 0)
  y <- expand.grid(H.B = x[1], M.B = x[2], S.B = x[3],
                   SI.W = x[4], P.W = x[5], W.W = x[6])

  dur.ev1 <- unlist(y[1:3] / 2)
  dur.ev2 <- unlist(y[4:6] / 2)
  dist.1 <- apply(y[1:3], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date, event.colors = event.colors.Buildings)
  dist.2 <- apply(y[4:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date, event.colors = event.colors.Water)
  results <- combine.lifeTimelines(dist.1[[1]], dist.2[[1]], dur.ev1, dur.ev2)

  #expected output is represented by the total number of interruption and the minimum distance between interventions
  z[1] <- sum(results[["duration"]])
  z[2] <- -min(abs(results$frequency[1:(length(results$frequency) - 1)] - results$frequency[2:length(results$frequency)]))

  #Cuz the function we are going to use for optimization aim to minimize a specific function
  #, we need to handle those objectives which we want to maximize.
  #To do this, we will negate the results, so the maximum value will become the minimum one
  #z[2] contains the results for the distance between two successive interventions and the defined objective for this is maximize()
  #Because of this, we negated the results.

  #define those variables which are static
  #The increase in job opportunities not only leads to a rise in office building floor area and local housing demand but also creates an incentive for incoming students to relocate to the area for education.
  #Based on this interaction pattern, we assume that for every 100% increase in office building capacity, the floor area of student dormitories and residential buildings will increase by 20% and 50%, respectively.
  ob.length <-x[7] # units: m
  ob.width <- 30 #units: m
  ob.height <- 4 #units: m
  ob.story <- 30

  sd.length <-24*(ob.length/24*0.2) # units: m
  sd.width <- 30 #units: m
  sd.height <- 4 #units: m
  sd.story <- 15

  rb.length <-30*(ob.length/24*0.5) # units: m
  rb.width <- 60 #units: m
  rb.height <- 4 #units: m
  rb.story <- 3

  wf.length <-40 # units: m
  wf.width <- 50 #units: m
  wf.height <- 6 #units: m
  wf.story <- 1

  materials <- LCI.materials
  structures.Options <- c("RC", "Timber", "Steel", "Steel_Frame")

  Energy.costs <- 1.08
  CO2.unitcost <- 37
  NOx.unitCost <- 9696
  SO2.unitCosts <- 19773

  #call the function we defined for life cycle analysis
  product.1.OfficeBuilding <- LCA.building(ob.length, ob.width, ob.height, ob.story, 0.25, 0.12, structures.Options[1], LCI.materials, dist.1[[1]])
  product.2.StudentDormitory <- LCA.building(sd.length, sd.width, sd.height, sd.story, 0.25, 0.12, structures.Options[2], LCI.materials, dist.1[[1]])
  product.3.ResidentialBuilding <- LCA.building(rb.length, rb.width, rb.height, rb.story, 0.06, 0.03, structures.Options[3], LCI.materials, dist.1[[1]])
  product.4.WaterSupplyFacility <- LCA.building(wf.length, wf.width, wf.height, wf.story, 0.04, 0.03, structures.Options[4], LCI.materials, dist.2[[1]])

  integrated.system <- as.data.frame(list(Energy = product.1.OfficeBuilding$Energy + product.2.StudentDormitory$Energy
                                                 + product.3.ResidentialBuilding$Energy + product.4.WaterSupplyFacility$Energy,
                                          CO2 = product.1.OfficeBuilding$CO2 + product.2.StudentDormitory$CO2
                                              + product.3.ResidentialBuilding$CO2 + product.4.WaterSupplyFacility$CO2,
                                          NOx = product.1.OfficeBuilding$NOx + product.2.StudentDormitory$NOx
                                              + product.3.ResidentialBuilding$NOx + product.4.WaterSupplyFacility$NOx,
                                          SO2 = product.1.OfficeBuilding$SO2 + product.2.StudentDormitory$SO2
                                              + product.3.ResidentialBuilding$SO2 + product.4.WaterSupplyFacility$SO2))

  integrated.system <- mutate(integrated.system,
                              Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx * NOx.unitCost + SO2 * SO2.unitCosts)/10^9)

  # take the results obtained and we will feed these to the output variable
  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)
}
# # Run NSGA-II optimization
r2 <- nsga2(fitness, idim = 7, odim = 7,
            generations = 10, popsize = 100,
            lower.bounds = c(2, 5, 15, 15, 1, 5, 24),
            upper.bounds = c(10, 15, 30, 30, 10, 15, 60))

#create plots from data
#transform results to dataframe
r2Results <- as.data.frame(r2$value)
outNames <- c("duration", "interv.dist", "energy", "co2", "nox", "so2", "cost")
colnames(r2Results) <- outNames

#Extract values of Pareto front and transform these into a dataframe
pareto3 <- as.data.frame(paretoFront(r2))
colnames(pareto3) <- outNames

# Plot Pareto front (duration vs cost)
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") +
  labs(title = "Pareto Front: Duration vs Cost",
       x = "Total Duration of Interventions",
       y = "Total Cost")

# Prepare data for parallel coordinate plot
input.params <- round(r2$par, 0)
all.results <- cbind(input.params, r2Results, r2$pareto.optimal)
colnames(all.results) <- c("H.B", "M.B", "S.B",
                           "SI.W", "P.W", "W.W",
                           "ob.length",
                           "duration", "interv.dist",
                           "energy", "co2", "nox", "so2", "cost",
                           "pareto")

# Create parallel coordinate plot
par(mfrow = c(1,1))
parcoord(all.results[, 1:15], var.label = TRUE,
         col = ifelse(all.results$pareto == TRUE, "indianred", "skyblue2"),
         lty = ifelse(all.results$pareto == TRUE, 1, 3),
         lwd = ifelse(all.results$pareto == TRUE, 3, 1))
title("Parallel Coordinate Plot of Design Parameters and Objectives")

# Print summary of Pareto-optimal solutions
cat("\nSummary of Pareto-optimal solutions:\n")
print(summary(pareto3))

# Print best solution for each objective
cat("\nBest solutions for each objective:\n")
for(obj in outNames) {
  if(obj == "interv.dist") {
    best <- pareto3[which.max(pareto3[[obj]]),]
  } else {
    best <- pareto3[which.min(pareto3[[obj]]),]
  }
  cat(paste0("Best ", obj, ":\n"))
  print(best)
  cat("\n")
}