# Install required packages (run once, then comment out)
# install.packages("timelineS")
# install.packages("lubridate")
# install.packages("ggplot2")
# install.packages("dplyr")
# install.packages("reshape2")
# install.packages("rPref")
# install.packages("mco")
# install.packages("MASS")

# setwd("D:/1st semester/LCA Final project-R")
# Load required libraries
library(timelineS)
library(lubridate)
library(ggplot2)
library(dplyr)
library(reshape2)
library(rPref)
library(mco)
library(MASS)

# Clear environment
rm(list = ls())

# Define helper functions
plot.timeline <- function(lifetime, events.name, start.Date, plot.Name)  # Plot timeline
  {
  dataP <- data.frame(Events = events.name, Event_Dates = ymd(start.Date) + years(lifetime))  # Create data frame
  timelineS(dataP, main = plot.Name, labels = events.name, label.direction = "up", label.position = 3)  # Plot timeline
}

dist.Events <- function(lifetime, events, start.Date, option.Name) {
  events <- sort(events, decreasing = T)  # Sort events
  distribution.events <- sapply(events, seq, from = 0, to = lifetime) # Generate distribution
  all.events <- melt(distribution.events)  # Melt data
  colnames(all.events) <- c("frequency", "event")  # Rename columns
  all.events <- all.events[order(all.events$frequency),]# Order data
  unique.events <- all.events[!duplicated(all.events$frequency),]# Remove duplicates
  unique.events$event[which(unique.events$frequency == 0)] <- "DC" # Set start event
  #plot timeliine
  # plot.timeline(unique.events$frequency, unique.events$event, start.Date, option.Name)
  return(unique.events) 
}

combine.lifeTimelines <- function(product.1, product.2, p1.dur, p2.dur) # Combine two timelines
  {
  base <- list(Names = c(), frequency = c(), duration = c())
  base$frequency <- sort(unique(c(product.1$frequency, product.2$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)) # Combine timelines
    {
    phase.1 <- product.1$event[which(product.1$frequency == base$frequency[index])] # Get phase 1
    phase.2 <- product.2$event[which(product.2$frequency == base$frequency[index])] # Get phase 2
    if(length(phase.1) != 0 && length(phase.2) != 0) # Combine phases
      {
      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)
}

# Define parameters
par(mfrow = c(2,1))##############################################
start.Date = "2025-02-01"
lifetime <- 130
events.Op1.MSBuilding = c(S = 20, M = 12, R = 40, END = lifetime)
duration.ev.MSBuilding <- c(S = 45, M = 10, R = 100)
events.Op1.OSBuilding <- c(M.o = 15, PR = 30, CR = 50, END = lifetime)
duration.ev.OSBuilding <- c(M.o = 10, PR = 45, CR = 100)
design.Options.MSBuilding <- list(desing.Op1.MSBuilding <- "2. Cast in place concrete Beam and Slab")
design.Options.OSBuilding <- list(desing.Op1.OSBuilding <- "3. One-Way RC Slab")

# Generate maintenance schedules
maintenance.MSBuilding <- dist.Events(lifetime, events.Op1.MSBuilding, start.Date, design.Options.MSBuilding)
maintenance.OSBuilding <- dist.Events(lifetime, events.Op1.OSBuilding, start.Date, design.Options.OSBuilding)
#print(maintenance.MSBuilding)

# Combine timelines
integrated.interv <- combine.lifeTimelines(maintenance.MSBuilding, maintenance.OSBuilding, duration.ev.MSBuilding, duration.ev.OSBuilding)

# Calculate total interruption duration
sum(integrated.interv$duration)


# ###################### Strategy1 change frequency as (12, *, 45)  ########################################
# events.Op1.MSBuilding = c(S = 20, M = 12, R = 45, END = lifetime)
# duration.ev.MSBuilding <- c(S = 45, M = 10, R = 100)
# events.Op1.OSBuilding <- c(M.o = 12, PR = 30, CR = 45, END = lifetime)
# duration.ev.OSBuilding <- c(M.o = 10, PR = 45, CR = 100)
# design.Options.MSBuilding <- list(desing.Op1.MSBuilding <- "2. Cast in place concrete Beam and Slab")
# design.Options.OSBuilding <- list(desing.Op1.OSBuilding <- "3. One-Way RC Slab")
# 
# # Generate maintenance schedules
# maintenance.MSBuilding <- dist.Events(lifetime, events.Op1.MSBuilding, start.Date, design.Options.MSBuilding)
# maintenance.OSBuilding <- dist.Events(lifetime, events.Op1.OSBuilding, start.Date, design.Options.OSBuilding)
# #print(maintenance.MSBuilding)
# 
# # Combine timelines
# integrated.interv <- combine.lifeTimelines(maintenance.MSBuilding, maintenance.OSBuilding, duration.ev.MSBuilding, duration.ev.OSBuilding)
# 
# # Calculate total interruption duration
# sum(integrated.interv$duration)



# ###################### Strategy2 change frequency as (12, 25, *) ########################################
# events.Op1.MSBuilding = c(S = 25, M = 12, R = 40, END = lifetime)
# duration.ev.MSBuilding <- c(S = 45, M = 10, R = 100)
# events.Op1.OSBuilding <- c(M.o = 12, PR = 25, CR = 50, END = lifetime)
# duration.ev.OSBuilding <- c(M.o = 10, PR = 45, CR = 100)
# design.Options.MSBuilding <- list(desing.Op1.MSBuilding <- "2. Cast in place concrete Beam and Slab")
# design.Options.OSBuilding <- list(desing.Op1.OSBuilding <- "3. One-Way RC Slab")
# 
# # Generate maintenance schedules
# maintenance.MSBuilding <- dist.Events(lifetime, events.Op1.MSBuilding, start.Date, design.Options.MSBuilding)
# maintenance.OSBuilding <- dist.Events(lifetime, events.Op1.OSBuilding, start.Date, design.Options.OSBuilding)
# #print(maintenance.MSBuilding)
# 
# # Combine timelines
# integrated.interv <- combine.lifeTimelines(maintenance.MSBuilding, maintenance.OSBuilding, duration.ev.MSBuilding, duration.ev.OSBuilding)
# 
# # Calculate total interruption duration
# sum(integrated.interv$duration)
# 



# 
###################### Strategy3 change frequency as (12, *, 40) ########################################
events.Op1.MSBuilding = c(S = 20, M = 12, R = 40, END = lifetime)
duration.ev.MSBuilding <- c(S = 45, M = 10, R = 100)
events.Op1.OSBuilding <- c(M.o = 12, PR = 30, CR = 40, END = lifetime)
duration.ev.OSBuilding <- c(M.o = 10, PR = 45, CR = 100)
design.Options.MSBuilding <- list(desing.Op1.MSBuilding <- "2. Cast in place concrete Beam and Slab")
design.Options.OSBuilding <- list(desing.Op1.OSBuilding <- "3. One-Way RC Slab")

# Generate maintenance schedules
maintenance.MSBuilding <- dist.Events(lifetime, events.Op1.MSBuilding, start.Date, design.Options.MSBuilding)
maintenance.OSBuilding <- dist.Events(lifetime, events.Op1.OSBuilding, start.Date, design.Options.OSBuilding)
#print(maintenance.MSBuilding)

# Combine timelines
integrated.interv <- combine.lifeTimelines(maintenance.MSBuilding, maintenance.OSBuilding, duration.ev.MSBuilding, duration.ev.OSBuilding)

# Calculate total interruption duration
sum(integrated.interv$duration)
# 
# 
# 
# 
# 





















# Define design exploration function
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.MSBuilding)
    dur.ev1 <- ev1/2
    for (j in 1:dim(events1)[1]) {
      ev2 <- unlist(events2[j, ])
      dist.2 <- dist.Events(lifetime, ev2, start.Date, design.Options.OSBuilding)
      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))
}




# set.seed(200)  # You can change 200 to any fixed number
# Generate design space
n.grid <- 5 # Number of grid points

# ##strategy 0
# events.grid.MSBuilding <- expand.grid(S = sample(seq(15,25, by = 1), n.grid),
#                                   M = sample(seq(5,19,1), n.grid),
#                                   R = sample(seq(20, 60, 1), n.grid)) # Design space for MSBuilding
# events.grid.OSBuilding <- expand.grid(M.o = sample(seq(5, 25, by = 1), n.grid),
#                                 PR = sample(seq(10, 50, 1), n.grid),
#                                 CR = sample(seq(40, 60, 1), n.grid)) # Design space for OSBuilding



# ##strategy 1
# events.grid.MSBuilding <- expand.grid(S = sample(seq(15,25, by = 1), n.grid),
#                                   M = sample(seq(5,19,1), n.grid),
#                                   R = sample(seq(25, 65, 1), n.grid)) # Design space for MSBuilding
# events.grid.OSBuilding <- expand.grid(M.o = sample(seq(7, 17, by = 1), n.grid),
#                                 PR = sample(seq(10, 50, 1), n.grid),
#                                 CR = sample(seq(35, 55, 1), n.grid)) # Design space for OSBuilding

# #strategy 2
# events.grid.MSBuilding <- expand.grid(S = sample(seq(15,35, by = 1), n.grid),
#                                       M = sample(seq(5,19,1), n.grid),
#                                       R = sample(seq(30, 50, 1), n.grid)) # Design space for MSBuilding
# events.grid.OSBuilding <- expand.grid(M.o = sample(seq(7, 17, by = 1), n.grid),
#                                       PR = sample(seq(20, 30, 1), n.grid),
#                                       CR = sample(seq(45, 55, 1), n.grid)) # Design space for OSBuilding

#strategy 3
events.grid.MSBuilding <- expand.grid(S = sample(seq(10,30, by = 1), n.grid),
                                      M = sample(seq(5,19,1), n.grid),
                                      R = sample(seq(30, 50, 1), n.grid)) # Design space for MSBuilding
events.grid.OSBuilding <- expand.grid(M.o = sample(seq(7, 17, by = 1), n.grid),
                                      PR = sample(seq(25, 35, 1), n.grid),
                                      CR = sample(seq(35, 45, 1), n.grid)) # Design space for OSBuilding











# Explore design space
response.space <- design.explore(events.grid.MSBuilding, events.grid.OSBuilding)  #design.explore is a function

# Define preferences and select best options
p <- low(dur) * high(dist.inter) # Define preferences
sky <- psel(response.space, p) # Select best options
pareto2 <- psel(response.space, p, top = nrow(response.space)) # Select Pareto-optimal solutions

# Visualize results
ggplot(response.space, aes(x = dur, y = dist.inter)) +
  geom_point(shape = 21) +
  geom_point(data = pareto2, size = 3, aes(color = factor(pareto2$.level)))


# Define function to show Pareto front
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 Pareto front
show_front(p)

p1 <- high(dur) * low(dist.inter)
show_front(p1)



# Life Cycle Inventory and Analysis (change the directory to yours)
LCI.materials <- read.csv("C:/Users/panji/Desktop/LCImaterials-2.csv")



############################define MSBuilding function####################
LCA.MSBuilding <- function(length, width, thickness, MSbeam.Option, MSslab.Option, materials, dist.event) {  
  SemiPrecast.prefab.MSbeam.Section <- 0.135
  SemiPrecast.cast_in_place.MSbeam.Section <- 0.045
  SemiPrecast.MSbeam.Section <- SemiPrecast.prefab.MSbeam.Section + SemiPrecast.cast_in_place.MSbeam.Section
  cast_in_place.MSbeam.Section <- 0.18
  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))
  # print(s2)
  # print(names(s2))
  # calculate the volume of the MSslab based on different materials strategies    
  if(MSslab.Option == "RC") {      
    if(MSbeam.Option =="none"){
      MSslab.volume <- length * width * thickness
      interventions.MSslab <- s2$DC + 0.3*s2$S + 0.7 * (if (!is.null(s2$R)) {s2$R} else {0}) + 0.15 * s2$M  #describe the reason for the weight of the maintenance
    } else if(MSbeam.Option =="RC"){
      MSslab.volume <- (length * width - 2 * (length + width) * 0.3)  * thickness
      interventions.MSslab <- s2$DC + 0.3*s2$S + 0.7 * (if (!is.null(s2$R)) {s2$R} else {0}) + 0.15 * s2$M 
    }
    
  } else if (MSslab.Option == "PRC+RC") {      
    MSslab.volume <-  0.14 * (length * width - 2 * (length + width) * 0.3)      
    interventions.MSslab <- s2$DC + 0.3*s2$S + 0.7 * (if (!is.null(s2$R)) {s2$R} else {0}) + 0.15 * s2$M  
  }
  #MSbeam options    
  if (MSbeam.Option == "PRC+RC") {      
    #get the numbers of beams      
    n1 <- round(length / 4, 0)
    n2 <- round(width / 4, 0)      
    interventions.beams <- s2$DC + 0.3*s2$S + 0.7 * (if (!is.null(s2$R)) {s2$R} else {0}) + 0.15 * s2$M
    #get the volume of the concrete for the prefab beams      
    beams.V <- (n1 * SemiPrecast.MSbeam.Section * width + n2 * SemiPrecast.MSbeam.Section * length )
  } else if (MSbeam.Option == "RC")  
  {      
    n1 <- round(length / 4, 0)
    n2 <- round(width / 4, 0)      
    beams.V <- (n1 * cast_in_place.MSbeam.Section * width + n2 * cast_in_place.MSbeam.Section * length )
    interventions.beams <- s2$DC + 0.3*s2$S + 0.7 * (if (!is.null(s2$R)) {s2$R} else {0}) + 0.15 * s2$M   
  } else if (MSbeam.Option == "none") {      
    n <- 0      
    beams.V <- 0      
    interventions.beams <- 0  
  } 
  # print(length(interventions.MSslab))
  # print(nrow(materials.split[[MSslab.Option]]))

  MSslab <- mutate(materials.split[[MSslab.Option]],  building.Q = MSslab.volume, interventions = 
                   interventions.MSslab)  
  if (!is.null(materials.split[[MSbeam.Option]])) {      
    MSbeam <- mutate(materials.split[[MSbeam.Option]], building.Q = beams.V, 
                   interventions = interventions.beams)  
    LCA.matrix <- rbind(MSslab, MSbeam)  
  } else {  
    LCA.matrix <- rbind(MSslab)    
  }   
  
  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)  
  print(LCA.matrix)
  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)) 
  print("interventions.MSslab value:")
  print(interventions.MSslab)
  return(LCA.results) 
  
}  



# Define parameters for LCA

#initialize the variables we need
b.length <- 8 # units: m  
b.width <- 4 #units m  
bd.depth <- 0.12 #units m  

#vectors, which will contain different design options for beam and for slab.
beam.Options <- c("RC", "PRC+RC", "none") 
slab.Options <- c("RC", "PRC+RC")  







############################## Define OSBuilding function ############################
LCA.OSBuilding <- function(length, width, height, thickness, Joist.Option, materials, dist.event) {  
  RC.Joist.Section <- 0.0272  # Maximum dimensions: Width 16 cm and height 17 cm    
  steel.Joist.unitWeight <- 2.94  # kg/m (Density * Cross-sectional area)
  
  # Split materials by scope
  materials.split <- split(materials, materials$scope)
  
  # Extract event frequencies from dist.event
  s <- summary(as.factor(dist.event$event))  
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors = FALSE, row.names = 1:length(s))
  
  # Slab volume & interventions
  OSslab.Option <- "RCO"
  OSslab.volume <- length * width * thickness
  interventions.OSslab <- s2$DC + 0.15*s2$M.o + 0.3*if (!is.null(s2$PR)) {s2$PR} else{ 0 } + 0.7* (if(! is.null(s2$CR)) {s2$CR} else {0})
  
  # Joist options
  if (Joist.Option == "RC Joist") {    
    n <- round(width / 0.75, 0)   
    joists.volume <- n * RC.Joist.Section * length
    interventions.joists <- s2$DC + 0.15*s2$M.o + 0.3*if (!is.null(s2$PR)) {s2$PR} else{ 0 } + 0.7* (if(! is.null(s2$CR)) {s2$CR} else {0})
  } else if (Joist.Option == "Steel Joist") {      
    n <- round(width / 0.75, 0) 
    joists.volume <- n * steel.Joist.unitWeight * length
    interventions.joists <- s2$DC + 0.15*s2$M.o + 0.3*if (!is.null(s2$PR)) {s2$PR} else{ 0 } + 0.7* (if(! is.null(s2$CR)) {s2$CR} else {0})
  } else if (Joist.Option == "none") {      
    n <- 0      
    joists.volume <- 0      
    interventions.joists <- 0
  } 
  
  # Hollow block volume
  if (Joist.Option %in% c("RC Joist", "Steel Joist")) {
    Block.Volume <- length * width * 0.15
  } else {
    Block.Volume <- 0
  }
  print(length(interventions.OSslab))
  print(nrow(materials.split[[OSslab.Option]]))
  # Define materials LCA calculations
  OSslab <- mutate(materials.split$RCO, OSBuilding.Q = OSslab.volume, interventions = interventions.OSslab) 
  Block <- mutate(materials.split$Block, OSBuilding.Q = Block.Volume, interventions = 0)
  
  if (!is.null(materials.split[[Joist.Option]])) {      
    joists <- mutate(materials.split[[Joist.Option]], OSBuilding.Q = joists.volume, interventions = interventions.joists)  
    LCA.matrix <- rbind(OSslab, joists, Block)  
  } else {  
    LCA.matrix <- rbind(OSslab, Block)    
  } 
  
  # Total material calculations
  LCA.matrix <- mutate(LCA.matrix, 
                       TotalMaterials.Q = quantities * OSBuilding.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)
  
  # Summarize results
  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)  
}



# Initialize variables
slab.length <- 10 
slab.width <- 6 
slab.thickness <- 0.1 
Block.height <- 0.15 

# Design options
Joist.Options <- c("RC Joist", "Steel Joist", "none")  



# Perform LCA for MSBuilding and OSBuilding
Option.MSBuilding <- LCA.MSBuilding(b.length, b.width, bd.depth, beam.Options[1], slab.Options[1], LCI.materials, maintenance.MSBuilding)
Option.OSBuilding <- LCA.OSBuilding(slab.length, slab.width, Block.height, slab.thickness, Joist.Options[3], LCI.materials, maintenance.OSBuilding)

# Combine LCA results
integrated.Design <- as.data.frame(list(Energy = Option.MSBuilding$Energy + Option.OSBuilding$Energy,
                                        CO2 = Option.MSBuilding$CO2 + Option.OSBuilding$CO2,
                                        NOx = Option.MSBuilding$NOx + Option.OSBuilding$NOx,
                                        SO2 = Option.MSBuilding$SO2 + Option.OSBuilding$SO2))

# Define unit costs
Energy.costs <- 0.04328
CO2.unitcost <- 160
NOx.unitCost <- 1857.5
SO2.unitCosts <- 7843.9

# 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########################################
fitness <- function(x) {
  z <- numeric(7)
  
  x <- round(x, 0)
  y <- expand.grid(S = x[1], M = x[2], R = x[3], M.o = x[4], PR = x[5], CR = 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)
  dist.2 <- apply(y[4:6], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  results <- combine.lifeTimelines(dist.1[[1]], dist.2[[1]], dur.ev1, dur.ev2)
  
  z[1] <- sum(results[["duration"]])
  z[2] <- -min(abs(results$frequency[1:(length(results$frequency) - 1)] - results$frequency[2:length(results$frequency)]))
  
  Distance <- x[7]
  #initialize the variables we need in MSbuilding
  b.length <- 8 # units: m  
  b.width <- 4 #units m  
  bd.depth <- 0.12 #units m  
  
  #vectors, which will contain different design options for beam and for slab.
  beam.Options <- c("RC", "PRC+RC", "none") 
  slab.Options <- c("RC", "PRC+RC")  
  
  # Initialize variables in OSBuilding
  slab.length <- 10 
  slab.width <- 6 
  slab.thickness <- 0.1 
  Block.height <- 0.15 
  
  # Design options
  Joist.Options <- c("RC Joist", "Steel Joist", "none")  
  
  materials <- LCI.materials
  
  product.1.MSBuilding <- LCA.MSBuilding(b.length, b.width, bd.depth, beam.Options[1], slab.Options[1], LCI.materials, dist.1[[1]])
  product.2.OSBuilding <- LCA.OSBuilding(slab.length, slab.width, Block.height, slab.thickness, Joist.Options[3], LCI.materials, dist.2[[1]])
  
  integrated.system <- as.data.frame(list(Energy = product.1.MSBuilding$Energy + product.2.OSBuilding$Energy,
                                          CO2 = product.1.MSBuilding$CO2 + product.2.OSBuilding$CO2,
                                          NOx = product.1.MSBuilding$NOx + product.2.OSBuilding$NOx,
                                          SO2 = product.1.MSBuilding$SO2 + product.2.OSBuilding$SO2))
  
  integrated.system <- mutate(integrated.system, 
                              Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx * NOx.unitCost + SO2 * SO2.unitCosts)/10^9)
  
  z[3] <- integrated.system$Energy
  z[4] <- integrated.system$CO2
  z[5] <- integrated.system$NOx
  z[6] <- integrated.system$SO2
  z[7] <- integrated.system$Costs
  
  return(z)
}


# set.seed(200)  # You can change 200 to any fixed number
# Run NSGA-II optimization
r2 <- nsga2(fitness, idim = 7, odim = 7,
            generations = 10, popsize = 100,
            lower.bounds = c(15, 2, 25, 6, 30, 90, 16),##7th parameter is the distance between 2 buildings
            upper.bounds = c(25, 8, 35, 14, 60, 180, 30))##7th parameter is the distance between 2 buildings

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

# Extract Pareto front
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("S", "M", "R", "M.o", "PR", "CR",
                           "distance", "duration", "interv.dist",
                           "energy", "co2", "nox", "so2", "cost",
                           "pareto")

# Create parallel coordinate plot
par(mfrow = c(1,1))
parcoord(all.results[, 1:14], 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")
}

# Save results
write.csv(pareto3, "pareto_optimal_solutions.csv", row.names = FALSE)