#ToDo Bounds anpassen, Interventionsbeschreibung, Konstanten setzen, 


## die optimierung muss drei design parameter variieren k?nnen. Sleepertype, Bicycle surface type und ..vom Tunnel
## Hans Karl Bartholdtsen 373092
## Anne Knappe 377357
## Muhammad 463116

rm (list = ls())

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


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) {
  #sort the events
  events <- sort(events, decreasing = T)
  #create the distribution of the events
  distribution.events <- sapply(events, seq, from = 0, to = lifetime)
  #all events unlisted
  all.events <- melt(distribution.events)
  colnames(all.events) <- c("frequency", "event")
  #sort the events ascending
  all.events <- all.events[order(all.events$frequency),]
  #get the unique sequence of events
  unique.events <- all.events[!duplicated(all.events$frequency),]
  unique.events$event[which(unique.events$frequency == 0)] <- "DC"
  #plot the timeline
  #plot.timeline(unique.events$frequency, unique.events$event, start.Date, option.Name)
  return(unique.events)
}



#######################################################################
#######################################################################
##### COMBINE TIMELINES
##Same lifetime
##Same Starting data
#######################################################################
#######################################################################


combine.lifeTimelines <- function(product.1, product.2, product.3, p1.dur, p2.dur, p3.dur) {
  #function body
  base <- list(Names = c(), frequency = c(), duration = c())
  base$frequency <- sort(unique(c(product.1$frequency,
                                  product.2$frequency,
                                  product.3$frequency)),
                         decreasing = FALSE)
  base$Names[1] <- "DC"
  base$duration[1] <- 0
  base$Names[length(base$frequency)] <- "END"
  base$duration[length(base$frequency)] <- 0
  for(index in 2:(length(base$frequency) - 1)) {
    phase.1 <- product.1$event[which(product.1$frequency ==
                                       base$frequency[index])]
    phase.2 <- product.2$event[which(product.2$frequency ==
                                       base$frequency[index])]
    phase.3 <- product.3$event[which(product.3$frequency ==
                                       base$frequency[index])]
    if(length(phase.1) != 0 && length(phase.2) != 0 && length(phase.3) != 0) {
      base$Names[index] <- paste(phase.1, phase.2, phase.3)
      base$duration[index] <- max(p1.dur[which(names(p1.dur) == phase.1)],
                                  p2.dur[which(names(p2.dur) == phase.2)],
                                  p3.dur[which(names(p2.dur) == phase.3)])
    } else if(length(phase.1) != 0 && length(phase.2) != 0 && length(phase.3) == 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 && length(phase.3) != 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)])
    } else if(length(phase.1) != 0 && length(phase.2) == 0 && length(phase.3) != 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)])
    } else if(length(phase.1) == 0 && length(phase.2) == 0 && length(phase.3) != 0) {
      base$Names[index] <- phase.3
      base$duration[index] <- p3.dur[which(names(p3.dur) == phase.3)]
    } else if(length(phase.1) != 0 && length(phase.2) == 0 && length(phase.3) == 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 && length(phase.3) == 0) {
      base$Names[index] <- phase.2
      base$duration[index] <- p2.dur[which(names(p2.dur) == phase.2)]
    }else {
      base$duration[index] <- 0
    }
  }
  return(base)
}


par(mfrow = c(2, 1))
start.Date = "2020-01-01"
lifetime <- 70

events.Op1.Tunnel =  c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration.ev.Tunnel <- c(TR = 2, M = 3, TSR = 4, PR = 5)
design.Options.Tunnel <- list(desing.Op1.Tunnel <- "Opt 1 - Tunnel")

events.Op1.TramTrack <- c(RPS = 15, GR = 10, RPR = 20, TP = 5, CL = 15, END = lifetime)
duration.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)
design.Options.TramTrack <- list(desing.Op1.TramTrack <- "Opt 1 - TramTrack")

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 30, END = lifetime)
duration.ev.Bicycle <- c(MI = 7, SP = 4, RS = 5)
design.Options.Bicycle <- list(desing.Op1.Bicycle <- "Opt 1 - Bicycle")

#design.Options.Tunnel <- list(desing.Op1 <- "Opt1 - Tunnel")



maintenance.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date, design.Options.TramTrack)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
maintenance.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date, design.Options.Bicycle)

#######################################################################
#####Variations of maintenance strategies
#######################################################################

#################################################
##########################Strategy 1 a tram track with two underpasses -> Tunnel and Bicycle Underpass
events.Op1.Tunnel =  c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration1.ev.Tunnel <- c(TR = 2, M = 3, TSR = 4, PR = 5)
design.Options.Tunnel <- list(desing.Op1.Tunnel <- "Opt 1 - Tunnel")

events.Op1.TramTrack <- c(RPS = 40, GR = 10, RPR = 20, TP = 5, CL = 15, END = lifetime)
duration1.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)
design.Options.TramTrack <- list(desing.Op1.TramTrack <- "Opt 1 - TramTrack")

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 55, END = lifetime)
duration1.ev.Bicycle <- c(MI = 2, SP = 3, RS = 8)
design.Options.Bicycle <- list(desing.Op1.Bicycle <- "Opt 1 - Bicycle")

maintenance1.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date, design.Options.TramTrack)
maintenance1.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
maintenance1.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date, design.Options.Bicycle)

#################################################
##########################Strategy 2 10 km tram track with concrete sleeper
########################### one large tunnel for cars and bicycles

duration2.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)

events.Op1.Tunnel =  c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration2.ev.Tunnel <- c(TR = 4, M = 6, TSR = 6, PR = 7)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 55, END = lifetime)
duration2.ev.Bicycle <- c(MI = 0, SP = 0, RS = 0)


maintenance2.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date, design.Options.TramTrack)
maintenance2.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
maintenance2.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date, design.Options.Bicycle)

#################################################
##########################Strategy 3 10 km tram track with timber sleeper
########################### combination of maintenance measurements of the tunnels

events.Op1.Tunnel =  c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration3.ev.Tunnel <- c(TR = 2, M = 2, TSR = 4, PR = 5)

events.Op1.TramTrack <- c(RPS = 20, GR = 10, RPR = 20, TP = 5, CL = 15, END = lifetime)
duration3.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 55, END = lifetime)
duration3.ev.Bicycle <- c(MI = 1, SP = 3, RS = 8)

maintenance3.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date, design.Options.TramTrack)
maintenance3.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
maintenance3.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date, design.Options.Bicycle)

#################################################
##########################Strategy 4 Material optimization
###########################  
events.Op1.Tunnel =  c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration4.ev.Tunnel <- c(TR = 2, M = 3, TSR = 4, PR = 5)

events.Op1.TramTrack <- c(RPS = 60, GR = 20, RPR = 20, TP = 15, CL = 15, END = lifetime)
duration4.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 55, END = lifetime)
duration4.ev.Bicycle <- c(MI = 2, SP = 3, RS = 7)

maintenance4.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date, design.Options.TramTrack)
maintenance4.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
maintenance4.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date, design.Options.Bicycle)

#######################################################################
#####Timeline integration
#######################################################################

integrated1.interv <- combine.lifeTimelines(maintenance1.Tunnel, maintenance1.TramTrack, maintenance1.Bicycle,
                                            duration1.ev.Tunnel, duration1.ev.TramTrack, duration1.ev.Bicycle)

interruption1 <- sum(integrated1.interv$duration)
interruption1
serviceLevel <- 100-interruption1/(lifetime*365)*100
serviceLevel

integrated2.interv <- combine.lifeTimelines(maintenance2.Tunnel, maintenance2.TramTrack, maintenance2.Bicycle,
                                            duration2.ev.Tunnel, duration2.ev.TramTrack, duration2.ev.Bicycle)

interruption2 <-sum(integrated2.interv$duration)
interruption2
serviceLevel <- 100-interruption2/(lifetime*365)*100
serviceLevel

integrated3.interv <- combine.lifeTimelines(maintenance3.Tunnel, maintenance3.TramTrack, maintenance3.Bicycle,
                                            duration3.ev.Tunnel, duration3.ev.TramTrack, duration3.ev.Bicycle)

interruption3 <-sum(integrated3.interv$duration)
interruption3
serviceLevel <- 100-interruption3/(lifetime*365)*100
serviceLevel

integrated4.interv <- combine.lifeTimelines(maintenance4.Tunnel, maintenance4.TramTrack, maintenance4.Bicycle,
                                            duration4.ev.Tunnel, duration4.ev.TramTrack, duration4.ev.Bicycle)

interruption4 <- sum(integrated4.interv$duration)
interruption4
serviceLevel <- 100-interruption4/(lifetime*365)*100
serviceLevel



#######################################################################
#####Automation of exploration
#######################################################################

design.explore <- function(events1, events2, events3) {
  results <- c()
  for(i in 1: dim(events1)[1]) {
    ev1 <- unlist(events1[i, ])
    dist.1 <- dist.Events(lifetime, ev1, start.Date, desing.Options.Tunnel)
    dur.ev1 <- ev1
    
    for (j in 1: dim(events2)[1]) {
      ev2 <- unlist(events2[j, ])
      dist.2 <- dist.Events(lifetime, ev2, start.Date, design.Options.TramTrack)
      dur.ev2 <- ev2
      
      for (f in 1: dim(events3)[1]) {
        ev3 <- unlist(events3[f, ])
        dist.3 <- dist.Events(lifetime, ev3, start.Date, design.Options.Bicycle)
        dur.ev3 <- ev3
        
        combined.lifetime <- combine.lifeTimelines(dist.1, dist.2, dist.3,
                                                   dur.ev1, dur.ev2, dur.ev3)
        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, dur = sum(combined.lifetime$duration), dist.inter = min.dist.int))
        
      }
    }
  }
  return(as.data.frame(results))
}

n.grid <- 4

events.grid.tunnel <- expand.grid(TR = sample(seq(11, 15, by = 1), n.grid),
                                  TSR = sample(seq(22, 40, 2), n.grid),
                                  PR = sample(seq(90, 120, 3), n.grid))
events.grid.tramTrack <- expand.grid(RPS = sample(seq(20, 40, by = 1), n.grid),
                                     GR = sample(seq(30,35,1), n.grid),
                                     RPR = sample(seq(10, 20, 2), n.grid))
events.grid.bicycle <- expand.grid(SP = sample(seq(12,20,1), n.grid),
                                   RS = sample(seq(30, 60, 1), n.grid))



response.space <- design.explore(events.grid.tunnel, events.grid.tramTrack, events.grid.bicycle)
response.space

p <- low(dur)* high(dist.inter)
p
sky <- psel(response.space, p)
sky

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

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)

#######################################################################
#######################################################################
##### Inventory 
########################### All
#######################################################################
start.Date = "2020-01-01"
lifetime<-70


#######################################################################
#######################################################################
##### Inventory 
########################### Tram Track
#######################################################################



LCA.track <- function(lengthRail, widthBody, railType.Option, sleeperType.Option, lengthSleeper, heightSleeper, widthSleeper, distanceSleeper, options, materials,lifetime, dist.event){
  
  
  #new line to get the interventions
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  ###Weight of rail dependent on rail type and rail length
  selectLine <- subset(options,componentType==railType.Option)
  weightM <- selectLine$weight
  weightRail <- weightM*lengthRail
  #print(paste0("Die Schienen haben ein Gewicht von ",weightRail," kg"))
  #Replacement rail, grinding, partial replacement if high accident probability
  interventions.rail <- s2$RPR+s2$GR*0.1+ 0.30 * (if (!is.null(s2$PR)) {s2$PR}
                                                  else {0})
  
  print(interventions.rail)
  ###CO2 NOx SO2 rail
  optionComponent <- selectLine$component
  optionMaterial <- selectLine$material
  
  selectSet <- subset(materials,scope==optionComponent) 
  selectMaterial <- subset(selectSet,material==optionMaterial)
  optionDurability <- selectMaterial$serviceLife
  
  #CO2 NOx SO2 rail 
  railRow <- list(CO2 = weightRail*interventions.rail*selectMaterial$kgCO2,
                  NOx = weightRail*interventions.rail*selectMaterial$kgNOx,
                  SO2 = weightRail*interventions.rail*selectMaterial$kgSO2Equi)
  matrix <- matrix(ncol = 3, nrow=0)
  matrix <- rbind(matrix,railRow)
  
  
  CO2Rail <-weightRail*selectMaterial$kgCO2*lifetime/optionDurability
  NOxRail <- weightRail*selectMaterial$kgNOx*lifetime/optionDurability
  SO2Rail <- weightRail*selectMaterial$kgSO2Equi*lifetime/optionDurability
  
  #print(matrix)
  
  #cost rail 
  costRail <- selectMaterial$lcc*lengthRail/1000
  
  ###Weight of the sleeper
  #Volume and amount of sleeper
  dis <- distanceSleeper+widthSleeper 
  div <- lengthRail/dis
  amountSleeper = ceiling(div)
  #print(paste0("Die Schwellen haben eine Anzahl von ",amountSleeper))
  volumeSleeper <-lengthSleeper*widthSleeper*heightSleeper*amountSleeper;
  #print(paste0("Die Schwellen haben ein Volumen von ",volumeSleeper," m³"))
  #weight
  #print(sleeperType.Option)
  selectLine <- subset(options,componentType==sleeperType.Option)
  # print(selectLine)
  weightM3 <- selectLine$weight
  weightSleeper <- weightM3 * volumeSleeper
  #print(paste0("Die Schwellen haben Gewicht von ",weightSleeper," kg"))
  
  #sleeper replacement
  interventions.sleeper <- s2$RPS
  
  
  ###CO2 NOx SO2 sleeper
  optionComponent <- selectLine$component
  optionMaterial <- selectLine$material
  
  
  selectSet <- subset(materials,scope==optionComponent) 
  selectMaterial <- subset(selectSet,material==optionMaterial)
  optionDurability <- selectMaterial$serviceLife
  sleeperRow <- list(CO2 = weightSleeper*interventions.sleeper*selectMaterial$kgCO2,
                     NOx = weightSleeper*interventions.sleeper*selectMaterial$kgNOx,
                     SO2 = weightSleeper*interventions.sleeper*selectMaterial$kgSO2Equi)
  matrix <- rbind(matrix,sleeperRow)
  
  CO2Sleeper <-weightSleeper*selectMaterial$kgCO2*lifetime/optionDurability
  NOxSleeper <- weightSleeper*selectMaterial$kgNOx*lifetime/optionDurability
  SO2Sleeper <- weightSleeper*selectMaterial$kgSO2Equi*lifetime/optionDurability
  
  
  #cost sleeper 
  costSleeper <-selectMaterial$lcc*lengthRail/1000
  
  ###Weight of ballast 
  #volume widthBody*heightsleeper*length
  volumeBallast <- widthBody*heightSleeper*lengthRail-volumeSleeper
  selectLine <- subset(options,componentType=="ballast")
  weightm3Ballast <- selectLine$weight
  weightBallast <- weightm3Ballast * volumeBallast
  #print(paste0("Die Ballast hat ein Gewicht von ",weightBallast," kg"))
  
  interventions.ballast <- s2$TP
  
  ###CO2 NOx SO2 ballast
  optionComponent <- selectLine$component
  optionMaterial <- selectLine$material
  
  
  selectSet <- subset(materials,scope==optionComponent) 
  selectMaterial <- subset(selectSet,material==optionMaterial)
  optionDurability <- selectMaterial$serviceLife
  
  ballastRow <- list(CO2 = weightBallast*interventions.ballast*selectMaterial$kgCO2,
                     NOx = weightBallast*interventions.ballast*selectMaterial$kgNOx,
                     SO2 = weightBallast*interventions.ballast*selectMaterial$kgSO2Equi)
  matrix <- rbind(matrix,ballastRow)
  CO2Ballast <-weightBallast*selectMaterial$kgCO2*lifetime/optionDurability
  NOxBallast <- weightBallast*selectMaterial$kgNOx*lifetime/optionDurability
  SO2Ballast <- weightBallast*selectMaterial$kgSO2Equi*lifetime/optionDurability
  
  #cost ballast 
  costBallast <-selectMaterial$lcc * lengthRail/1000
  #print(matrix)
  
  LCA.results<-  list(CO2 = CO2Rail+CO2Sleeper+CO2Ballast,
                      NOx = NOxRail+NOxSleeper+NOxBallast,
                      SO2 = SO2Rail+SO2Sleeper+SO2Ballast
                      #Lcc = costRail+costSleeper+costBallast
                      )
  #print(LCA.results)
  return(LCA.results)
}




#######################################################################
#######################################################################
##### Inventory 
########################### Tunnel
#######################################################################


LCA.tunnel <- function(length, width, height, thickness, girder.Option, deck.Option, materials, dist.event) {
  
  
  #new line to get the interventions
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  
  prefab.girder.Section <- 0.78
  steel.girders.unitWeight <- 317 #the weight for HEM800 steel profile
  asphalt.Q <- length * width * thickness
  materials.split <- split(materials, materials$scope)
  
  # calculate the volume of the deck based on different materials strategies
  if(deck.Option == "PCS") {
    deck.volume <- length * width * height
    interventions.deck <- 2.5
  } else if (deck.Option == "CIS") {
    deck.volume <- 0.5 * length * width * height
    interventions.deck <- 2
  } else if (deck.Option == "FRP") {
    deck.volume <- 0.2 * length * width * height
    interventions.deck <- 1
  }
  
  #girder options
  if (girder.Option == "RB") {
    #get the numbers of girders
    n <-round(width / 3.75, 0)
    interventions.girders <- 2
    #get the volume of the concrete for the prefab girders
    girders.V <- n * prefab.girder.Section * length
  } else if (girder.Option == "RBS") {
    n <- round(width / 3, 0)
    girders.V <- n * steel.girders.unitWeight * length
    interventions.girders <- 2
  } else if (girder.Option == "none") {
    n <- 0
    girders.V <- 0
    interventions.girders <- 0
  }
  #print(materials.split$asphalt)
  
  asphalt <- mutate(materials.split$asphalt, bridge.Q = asphalt.Q, interventions = 12)
  #print("test")
  #print(deck.Option)
  #print(materials.split$deck.Option)
  #print(interventions.deck)
  #print(deck.volume)
  #print("s: ")
  #print(s)
  #print("s2: ")
  #print(s2)
  #print("s2TR: ")
  #print(s2$TR)
  #print(s2$M*0.2)
  s2$M <- s2$DC
  s2$TR <- s2$M
  s2$TSR <- s2$TR
  s2$PR <- s2$TSR
  interventions.Tunnel <- s2$TR + s2$M*0.2 + s2$TSR*0.35+s2$PR*0.1
  #print(interventions.Tunnel)
  deck <- mutate(materials.split[[deck.Option]], bridge.Q = deck.volume, interventions = interventions.Tunnel)
  #print("test2")
  if (!is.null(materials.split[[girder.Option]])) {
    girders <- mutate(materials.split[[girder.Option]], bridge.Q = girders.V, interventions = interventions.girders)
    LCA.matrix <- rbind(deck, girders, asphalt)
  } else {
    LCA.matrix <- rbind(deck, asphalt)
  }
  
  LCA.matrix <- mutate(LCA.matrix, TotalMaterials.Q = quantities * bridge.Q / 1000,
                       materials.LC = TotalMaterials.Q * interventions,
                       Energy.LC = materials.LC * energy,
                       CO2.LC = materials.LC * CO2,
                       NOx.LC = materials.LC * NOx,
                       SO2.LC = materials.LC * SO2)
  
  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)
}

#######################################################################
#######################################################################
##### Inventory 
########################### Bicycle Underpass
#######################################################################


LCA.bicycle <- function(material.Option, materials,lifetime,  dist.event) {
  # function’s body
  s <- summary(as.factor(dist.event$event))
  s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors=F, row.names = 1:length(s))
  durability <-90
  selectLine <- subset(materials,material==material.Option)
  amount <- selectLine$amount;
  interventions.rail <- s2$MI+s2$SP*0.5+s2$RS*0.2
  LCA.results <- list(CO2 = amount*interventions.rail*selectLine$CO2*lifetime/durability,
                      NOx = amount*interventions.rail**selectLine$NOX*lifetime/durability,
                      SO2 = amount*interventions.rail**selectLine$SO2*lifetime/durability)
  return(LCA.results);
}




#######################################################################
#######################################################################
##### Inventory 
########################### Integration
####################################################################### 


#########Tunnel
LCI.materials <- read.csv("MaterialsTunnel.csv")
tunnelsupp.Options <- c("RB", "RBS", "none")
tunnelseg.options <- c("PCS", "CIS", "SBS")

#########Track
LCI.materialsTrack <- read.csv("materials.csv")
LCO.options <- read.csv("options.csv")
sleeperType.Options <- c("timber", "concrete", "steel","plastic")
railType.Options <- c("vignole rail", "grooved rail")

########Bicycle

mat.Options <- c("Asphalt", "Concrete","Tiles")
LCI.materialsBicy <- read.csv("bicycleMaterials.csv")
design.Options.Bicycle <- list(desing.Op1.Bicycle <- "Opt 1 - Bicycle")

########################### Strategy 1
########################### 10 km tram track with concrete sleeper and two underpasses --> 
###########################tunnel with reinforced concrete segments and a bicycle underpass with concrete road
#######################################################################

events.Op1.Tunnel = c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration.ev.Tunnel <- c(TR = 2, M = 3, TSR = 4, PR = 5)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)

events.Op1.TramTrack <- c(RPS = 15, GR = 10, RPR = 20, TP = 5, CL = 15, END = lifetime)
duration.ev.TramTrack <- c(RPS = 2, GR = 1, RPR = 1, TP = 1, CL = 1)
maintenance.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date,design.Options.TramTrack)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 55, END = lifetime)
duration.ev.Bicycle <- c(MI = 2, SP = 3, RS = 8)
maintenance.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date,design.Options.Bicycle)

b.length <- 16 # units: m
b.width <- 15 #units m
bd.depth <- 0.25 #units m
asphalt.tk <- 0.12 #units m

Option.track <- LCA.track(10, 3.5, railType.Options[1], sleeperType.Options[2], 2.6, 0.16, 0.26, 0.6, LCO.options, LCI.materialsTrack, 50, maintenance.TramTrack)
Option.bicycle <- LCA.bicycle(mat.Options[1], LCI.materialsBicy, 60, maintenance.Bicycle)
Option.tunnel <- LCA.tunnel(b.length, b.width, bd.depth, asphalt.tk, tunnelsupp.Options[1], tunnelseg.options[1], LCI.materials, maintenance.Tunnel)


integrated.Design1 <- as.data.frame(list(x =1,
                                         CO2 = Option.tunnel$CO2 + Option.track$CO2+Option.bicycle$CO2,
                                         NOx = Option.tunnel$NOx + Option.track$NOx+Option.bicycle$NOx,
                                         SO2 = Option.tunnel$SO2 + Option.track$SO2+Option.bicycle$SO2))
integrated.Design1


########################### Strategy 2
########################### 10 km tram track with concrete sleeper
########################### one large tunnel for cars and bicycles
#######################################################################

b.length <- 16 # units: m
b.width <- 20 #units m higher width because tunnel and bicycle underpass together
bd.depth <- 0.25 #units m
asphalt.tk <- 0.12 #units m
events.Op1.Tunnel = c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration.ev.Tunnel <- c(TR = 8, M = 3, TSR = 4, PR = 2)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)
railType.Options[0]
Option.track <- LCA.track(10, 3.5, railType.Options[1], sleeperType.Options[1], 2.6, 0.16, 0.26, 0.6, LCO.options, LCI.materialsTrack, 50, maintenance.TramTrack)
#Option.bicycle <- LCA.bicycle(mat.Options[1], LCI.materialsBicy, 60, maintenance.Bicycle)
Option.tunnel <- LCA.tunnel(b.length, b.width, bd.depth, asphalt.tk, tunnelsupp.Options[1], tunnelseg.options[1], LCI.materials, maintenance.Tunnel)


integrated.Design2 <- as.data.frame(list(x =2,
                                         CO2 = Option.tunnel$CO2 + Option.track$CO2,
                                         NOx = Option.tunnel$NOx + Option.track$NOx,
                                         SO2 = Option.tunnel$SO2 + Option.track$SO2))

########################### Strategy 3
###########################  Combine measures at tunnel and bicycle
########################### optimize the track input values --> wooden sleeper, small body
#######################################################################

b.length <- 16 # units: m
b.width <- 15 #units m higher width because tunnel and bicycle underpass together
bd.depth <- 0.2 #units m
asphalt.tk <- 0.12 #units m


events.Op1.Tunnel = c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration.ev.Tunnel <- c(TR = 2, M = 5, TSR = 4, PR = 5)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 30, END = lifetime)
duration.ev.Bicycle <- c(MI = 0, SP = 3, RS = 5)
maintenance.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date,design.Options.Bicycle)

Option.track <- LCA.track(10, 3, railType.Options[1], sleeperType.Options[3], 2.6, 0.16, 0.26, 0.6, LCO.options, LCI.materialsTrack, 50, maintenance.TramTrack)
Option.bicycle <- LCA.bicycle(mat.Options[1], LCI.materialsBicy, 60, maintenance.Bicycle)
Option.tunnel <- LCA.tunnel(b.length, b.width, bd.depth, asphalt.tk, tunnelsupp.Options[1], tunnelseg.options[1], LCI.materials, maintenance.Tunnel)
Option.track
Option.bicycle
Option.track

integrated.Design3 <- as.data.frame(list(x =3,
                                         CO2 = Option.tunnel$CO2 + Option.track$CO2+Option.bicycle$CO2,
                                         NOx = Option.tunnel$NOx + Option.track$NOx+Option.bicycle$NOx,
                                         SO2 = Option.tunnel$SO2 + Option.track$SO2+Option.bicycle$SO2))



########################### Strategy 4
###########################  Combine measures at tunnel and bicycle
########################### the tunnel is now a Net Element 
########################### the tunnel i now a Net Element -> that means the track is in the tunnel and a road is above 
#######################################################################

b.length <- 30 # units: m
b.width <- 15 #units m higher width because tunnel and bicycle underpass together
bd.depth <- 0.2 #units m
asphalt.tk <- 0.12 #units m

events.Op1.Tunnel = c(TR = 15, M = 10, TSR = 40, PR = 50, END = lifetime)
duration.ev.Tunnel <- c(TR = 7, M = 2, TSR = 3, PR = 1)
maintenance.Tunnel <- dist.Events(lifetime, events.Op1.Tunnel, start.Date, design.Options.Tunnel)

events.Op1.TramTrack <- c(RPS = 15, GR = 10, RPR = 20, TP = 5, CL = 15, END = lifetime)
duration.ev.TramTrack <- c(RPS = 5, GR = 3, RPR = 2, TP = 5, CL = 2)
maintenance.TramTrack <- dist.Events(lifetime, events.Op1.TramTrack, start.Date,design.Options.TramTrack)

events.Op1.Bicycle <- c(MI = 10, SP = 15, RS = 30, END = lifetime)
duration.ev.Bicycle <- c(MI = 2, SP = 3, RS = 8)
maintenance.Bicycle <- dist.Events(lifetime, events.Op1.Bicycle, start.Date,design.Options.Bicycle)


Option.track <- LCA.track(10, 3, railType.Options[1], sleeperType.Options[2], 2.6, 0.16, 0.26, 0.6, LCO.options, LCI.materialsTrack, 50, maintenance.TramTrack)
Option.bicycle <- LCA.bicycle(mat.Options[1], LCI.materialsBicy, 60, maintenance.Bicycle)
Option.tunnel <- LCA.tunnel(b.length, b.width, bd.depth, asphalt.tk, tunnelsupp.Options[1], tunnelseg.options[1], LCI.materials, maintenance.Tunnel)
Option.tunnel
Option.track
Option.bicycle

integrated.Design4 <- as.data.frame(list(x =4,
                                         CO2 = Option.tunnel$CO2 + Option.track$CO2 + Option.bicycle$CO2,
                                         NOx = Option.tunnel$NOx + Option.track$NOx + Option.bicycle$NOx,
                                         SO2 = Option.tunnel$SO2 + Option.track$SO2 + Option.bicycle$SO2))

integrated.Design1 <- mutate(integrated.Design1,
                             Costs = (  CO2 * 26 +
                                          NOx * 42 +
                                          SO2 * 85)/10^3)

integrated.Design2 <- mutate(integrated.Design2,
                             Costs = (  CO2 * 26 +
                                          NOx * 42 +
                                          SO2 * 85)/10^3)

integrated.Design3 <- mutate(integrated.Design3,
                             Costs = (  CO2 * 26 +
                                          NOx * 42 +
                                          SO2 * 85)/10^3)


integrated.Design4 <- mutate(integrated.Design4,
                             Costs = (  CO2 * 26 +
                                          NOx * 42 +
                                          SO2 * 85)/10^3)
integrated.Design <- rbind(integrated.Design1,integrated.Design2, integrated.Design3, integrated.Design4)



###########################  Visualization 
###########################  Boxplots
#######################################################################
integrated.Design4
test <-unlist(integrated.Design1)
test
results <- as.data.frame(rbind(Option1 = unlist(integrated.Design1),
                               Option2 = unlist(integrated.Design2),
                               Option3 = unlist(integrated.Design3),
                               Option4 = unlist(integrated.Design4)))

barplot(results$CO2, names.arg = row.names(results), main ="CO2")
barplot(results$NOx, names.arg = row.names(results), main ="NOx")
barplot(results$SO2, names.arg = row.names(results), main ="SO2")
barplot(results$Costs, names.arg = row.names(results), main ="Costs")



###########################  Visualization 
###########################  radar diagram
#######################################################################
results$x <- NULL
results.df = results %>%
  rownames_to_column(var = "option") %>%
  mutate_each(funs(rescale), -option) %>%
  melt(id.vars = c('option'), measure.vars = colnames(results)) %>%
  arrange(option)
results.df %>%
  ggplot(aes(x=option, y=value, group=variable, color=variable)) +
  geom_polygon(fill=NA) +
  coord_polar() + theme_bw() + facet_wrap(~ variable) +
  theme(axis.text.x = element_text(size = 5))
###########################  Visualization 
###########################  Graphs
#######################################################################




data <- data.frame(x = c(1,2,3,4),
                   y1 = c(integrated.Design1$CO2/10, integrated.Design1$NOx  , integrated.Design1$SO2,integrated.Design1$Costs),
                   y2 = c(integrated.Design2$CO2/10, integrated.Design2$NOx , integrated.Design2$SO2,integrated.Design2$Costs),
                   y3 = c(integrated.Design3$CO2/10, integrated.Design3$NOx , integrated.Design3$SO2,integrated.Design3$Costs),
                   y4 = c(integrated.Design4$CO2/10, integrated.Design4$NOx , integrated.Design4$SO2,integrated.Design4$Costs))
data
data_ggp <- data.frame(x = data$x,
                       y = c(data$y1, data$y2, data$y3, data$y4),
                       Scenario = c(rep("Scenario 1", nrow(data)),
                                    rep("Scenario 2", nrow(data)),
                                    rep("Scenario 3", nrow(data)),
                                    rep("Scenario 4", nrow(data))))

ggp <- ggplot(data_ggp, aes(x, y, col = Scenario)) + xlab("1 = CO2 , 2 = NOx , 3 = SO2, 4 = Costs")  +     
  ylab("kg CO2/Nox/SO2 (CO2*10)")   +      
  geom_line()
ggp  


#############################################################################################################################
############################# MULTI-OBJECTIVE OPTIMIZATION
#############################################################################################################################
## function 

fitness <- function(x){
  # dimension of the output vector
  z<- numeric(6)
  
  #start if function body
  # ABSCHNITT MUSS NOCH UPDATET WERDEN FUNKTIONEN DIE DIE INTERVENTIONS BESCHREIBEN
  
  # M?glichkeitsraum aufspannen 
  x <- round(x, 0)
  #print("Hier folgt X: ")
  #print(x)
  y <- expand.grid(RPS = x[1], GR = x[2], RPR = x[3], TP = x[4], CL = x[5],
                   TR = x[6], M = x[7], TSR = x[8], PR = x[9],
                   MI = x[10], SP = x[11], RS = x[12])
  
  dur.ev1 <- unlist(y[1:5])
  dur.ev2 <- unlist(y[6:9])
  dur.ev3 <- unlist(y[10:12])
  
  dist.1 <- apply(y[1:5], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  dist.2 <- apply(y[6:9], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  dist.3 <- apply(y[10:12], 1, FUN = dist.Events, lifetime = lifetime, start.Date = start.Date)
  results <- combine.lifeTimelines(dist.1[[1]], dist.2[[1]], dist.3[[1]], dur.ev1, dur.ev2, dur.ev3)
  #print(dist.3)
  
  resultsFrame <- as.data.frame(results)
  print(resultsFrame)
  #print(results)
  #bis hierhin soll z1 und z2 bestimmen. 
  
  #set all constant vaLUES for inventory  ABSCHNITT MUSS NOCH GEMACHT WERDEN
  
  #der eine Wert der variiert
  distanceSleeper <- x[13] #was 0.6 bounded now between 0.2 and 1.0
  widthBody <-x[14]
  asphalt.tk <- x[15] #units m
  
  b.length <- 16 # units: m
  b.width <- 15 #units m
  bd.depth <- 0.25 #units m
  materials <- LCI.materials

  Energy.costs <- 0.128
  CO2.unitcost <- 26
  NOx.unitCost <- 42
  SO2.unitCosts <- 85
  
  maintenance.TramTrack <- dist.1[[1]]
  maintenance.Tunnel <- dist.2[[1]]
  maintenance.Bicycle <- dist.3[[1]]
 
  
  #FUNKTION F?R INVENTORY WIEVIEL WOVON BENUTZT WIRD
  
  #wir brauchen hier stattdessen LCA.track, LCA.bicycle und LCA.tunnel mit den richtigen Inputwerten von dar?ber.
  
  product.1.track <- LCA.track(10, widthBody, railType.Options[1], sleeperType.Options[2], 2.6, 0.16, 0.26, distanceSleeper, LCO.options, LCI.materialsTrack, 50, maintenance.TramTrack )
  product.2.tunnel <- LCA.tunnel(b.length, b.width, bd.depth, asphalt.tk, tunnelsupp.Options[1], tunnelseg.options[1], LCI.materials, maintenance.Tunnel)
  product.3.bicycle <- LCA.bicycle(mat.Options[1],LCI.materialsBicy, 70, maintenance.Bicycle)
  
  #die Emissionen aus den drei Subsystemen addieren
  integrated.system <- as.data.frame(list(
  #                                        Energy = product.1.track$Energy + product.2.tunnel$Energy + product.3.bicycle$Energy,
                                          CO2 = product.1.track$CO2 + product.2.tunnel$CO2 + product.3.bicycle$CO2,
                                          NOx = product.1.track$NOx + product.2.tunnel$NOx + product.3.bicycle$NOx,
                                          SO2 = product.1.track$SO2 + product.2.tunnel$SO2 + product.3.bicycle$SO2))
  # Kosten als Funktion der Emissionen
  integrated.system <- mutate(integrated.system,
                              Costs = (
  #                              Energy * Energy.costs +
                                         CO2*CO2.unitcost +
                                         NOx * NOx.unitCost +
                                         SO2 * SO2.unitCosts)/10^9)
  
  #end of function body
  
  #outputvector with total duration of interventions, min distance between interventions, energy, CO2, NOx, SO2, Costs
  z[1] <- sum(results[["duration"]])
  z[2] <- -min(abs(results$frequency[1:(length(results$frequency) - 1)] -results$frequency[2:length(results$frequency)]))
  #z[3] <- integrated.system$Energy
  z[3] <- integrated.system$CO2
  z[4] <- integrated.system$NOx
  z[5] <- integrated.system$SO2
  z[6] <- integrated.system$Costs
  
  return(z)
}


#function that optimizes our system
#GENERATION und POPSIZE f?r final run erh?hen

roptimized <- nsga2(fitness, idim = 15, odim = 6,
                          generations=30, popsize=100,
                          lower.bounds=c(  5, 2,  4,  1, 25,  10,  19,  5, 14,  5, 12, 27,  0.2, 3.5, 0.12),
                          upper.bounds= c(10, 6, 10,  6, 35,  20,  20, 15, 15, 10, 25, 37,  1, 6, 0.24))



#graphical representation
roptimizedInputs <- as.data.frame(roptimized$par)

roptimizedResults <- as.data.frame(roptimized$value)
outNames <- c("totalduration", "interv.dist.","co2","nox","so2", "cost")
colnames(roptimizedResults) <- outNames

roptimizedResults

pareto3 <- as.data.frame(paretoFront(roptimized))
colnames(pareto3) <- outNames
pareto3


#quick fix more than one data point
#roptimizedResults[4,3] <-702144760
#roptimizedResults[4,6] <- 19.00229
#pareto3[4,3] <- 702144760
#pareto3[4,6] <- 19.00229


#name der Bezeichnung muss ohne Leerzeichen sein.
# zwei Ziele auf den beiden Achsen
ggplot(roptimizedResults, aes(x = cost , y = totalduration)) + geom_point(shape = 21) +geom_point(data = pareto3, size = 3, color="red") + geom_line(data = pareto3, color="blue")


#Header fehlt.
# zwei Ziele auf den beiden Achsen
ggplot(roptimizedResults, aes(x = cost , y = interv.dist.)) + geom_point(shape = 21) +geom_point(data = pareto3, size = 3, color="red") + geom_line(data = pareto3, color="blue")


# accumulated impacts
input.params <- round(roptimized$par, 0)
all.results <- cbind(input.params, roptimizedResults, roptimized$pareto.optimal)
colnames(all.results) <- c("RPS","GR","RPR","TP","CL", "TR", "M", "TSR", "PR","MI", "SP", "RS","sl.Dist","w.Body","asphaltTk","tot.dur","int.dist.","co2","nox","so2", "cost","pareto")

all.results[,c(7,15)] <-list(NULL)

#presentation of the 13 input variables
par(mfrow = c(1,1))
parcoord(all.results[, 1:13], 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))

#presentation of the 6 output variables
par(mfrow = c(1,1))
parcoord(all.results[, 14: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))

#presentation of all variables
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))