rm (list = ls()) library(timelineS) library(lubridate) library(ggplot2) library(dplyr) library(reshape2) library(rPref) library(mco) library(MASS) # ----------------------------- # Functions # ----------------------------- plot.timeline <- function(lifetime, events.name, start.Date, plot.Name) { dataP <- data.frame( Events = events.name, Event_Dates = ymd(start.Date) + years(lifetime) ) timelineS( dataP, main = plot.Name, labels = events.name, label.direction = "up", label.position = 3 ) } dist.Events <- function(lifetime, events, start.Date, option.Name) { events <- sort(events, decreasing = TRUE) distribution.events <- sapply(events, seq, from = 0, to = lifetime) all.events <- melt(distribution.events) colnames(all.events) <- c("frequency", "event") all.events <- all.events[order(all.events$frequency), ] unique.events <- all.events[!duplicated(all.events$frequency), ] unique.events$event[unique.events$frequency == 0] <- "DC" return(unique.events) } # ----------------------------- # Integrate 3 Systems # ----------------------------- combine.lifeTimelines <- function(product.1, product.2, product.3, p1.dur, p2.dur, p3.dur) { 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])] #case 1: if both have maintenance, combine into 1 event if (sum(c(length(phase.1), length(phase.2), length(phase.3)) != 0) >= 2) { base$Names[index] <- paste(na.omit(c(phase.1, phase.2, phase.3)), collapse = ", ") base$duration[index] <- max( c( if (!is.null(phase.1)) p1.dur[which(names(p1.dur) == phase.1)] else NA, if (!is.null(phase.2)) p2.dur[which(names(p2.dur) == phase.2)] else NA, if (!is.null(phase.3)) p3.dur[which(names(p3.dur) == phase.3)] else NA ), na.rm = TRUE ) #case 2: only 1 product has an event } else if(sum(c(length(phase.1), length(phase.2), length(phase.3)) != 0) ==1) { base$Names[index] <- paste(na.omit(c(phase.1, phase.2, phase.3)), collapse = ", ") base$duration[index] <- max( c( if (!is.null(phase.1)) p1.dur[which(names(p1.dur) == phase.1)] else NA, if (!is.null(phase.2)) p2.dur[which(names(p2.dur) == phase.2)] else NA, if (!is.null(phase.3)) p3.dur[which(names(p3.dur) == phase.3)] else NA ), na.rm = TRUE ) #case 3: no products have an event } else { base$duration[index] <- 0 } } return(base) } # ----------------------------- # First Maintenance Strategy # ----------------------------- # Define parameters lifetime <- 30 start.Date <- "2025-01-01" events.Op1.Highway = c(M_h = 5, R_h = 15, O_h = 30, END = lifetime) duration.ev.Highway <- c(M_h = 2, R_h = 8, O_h = 20) events.Op1.Urban <- c(M_u = 6, R_u = 18, U_u = 30, END = lifetime) duration.ev.Urban <- c(M_u = 2, R_u = 6, U_u = 15) events.Op1.Drain <- c(C_d = 7, R_d = 30, emp = 30 , END = lifetime) duration.ev.Drain <- c(C_d = 1, R_d = 6, emp = 0) design.Options.Highway <- list(desing.Op1.Highway <- "Highway") design.Options.Urban <- list(desing.Op1.Urban <- "Urban") design.Options.Drain <- list(desing.Op1.Drain <- "Drain") # Generate timelines maintenance.H <- dist.Events(lifetime, events.Op1.Highway, start.Date, "Highway") maintenance.U <- dist.Events(lifetime, events.Op1.Urban, start.Date, "Urban") maintenance.D <- dist.Events(lifetime, events.Op1.Drain, start.Date, "Drainage") # Combine timelines integrated <- combine.lifeTimelines(maintenance.H, maintenance.U, maintenance.D, duration.ev.Highway, duration.ev.Urban, duration.ev.Drain) # Calculate total interruption duration cat("Total system downtime in days over", lifetime, "years:\n") sum(integrated$duration) # Plot timelines par(mfrow = c(3, 1)) plot.timeline(maintenance.H$frequency, maintenance.H$event, start.Date, "Highway") plot.timeline(maintenance.U$frequency, maintenance.U$event, start.Date, "Urban") plot.timeline(maintenance.D$frequency, maintenance.D$event, start.Date, "Drainage") # ----------------------------- # Second Maintenance Strategy # ----------------------------- lifetime <- 30 start.Date <- "2025-01-01" events.Op2.Highway = c(M_h = 6, R_h = 18, O_h = 30, END = lifetime) duration.ev2.Highway <- c(M_h = 2, R_h = 8, O_h = 20) events.Op2.Urban <- c(M_u = 6, R_u = 18, U_u = 30, END = lifetime) duration.ev2.Urban <- c(M_u = 2, R_u = 6, U_u = 15) events.Op2.Drain <- c(C_d = 7, R_d = 30, emp = 30 , END = lifetime) duration.ev2.Drain <- c(C_d = 1, R_d = 6, emp = 0) design.Options.Highway2 <- list(desing.Op2.Highway <- "Highway") design.Options.Urban2 <- list(desing.Op2.Urban <- "Urban") design.Options.Drain2 <- list(desing.Op2.Drain <- "Drain") maintenance.H2 <- dist.Events(lifetime, events.Op2.Highway, start.Date, "Highway") maintenance.U2 <- dist.Events(lifetime, events.Op2.Urban, start.Date, "Urban") maintenance.D2 <- dist.Events(lifetime, events.Op2.Drain, start.Date, "Drainage") integrated <- combine.lifeTimelines(maintenance.H2, maintenance.U2, maintenance.D2, duration.ev2.Highway, duration.ev2.Urban, duration.ev2.Drain) cat("Total system downtime in days over", lifetime, "years:\n") sum(integrated$duration) par(mfrow = c(3, 1)) plot.timeline(maintenance.H2$frequency, maintenance.H2$event, start.Date, "Highway") plot.timeline(maintenance.U2$frequency, maintenance.U2$event, start.Date, "Urban") plot.timeline(maintenance.D2$frequency, maintenance.D2$event, start.Date, "Drainage") # ----------------------------- # Third Maintenance Strategy # ----------------------------- lifetime <- 30 start.Date <- "2025-01-01" events.Op3.Highway = c(M_h = 6, R_h = 18, O_h = 30, END = lifetime) duration.ev3.Highway <- c(M_h = 2, R_h = 8, O_h = 20) events.Op3.Urban <- c(M_u = 7, R_u = 20, U_u = 30, END = lifetime) duration.ev3.Urban <- c(M_u = 2, R_u = 6, U_u = 15) events.Op3.Drain <- c(C_d = 5, R_d = 30, emp = 30 , END = lifetime) duration.ev3.Drain <- c(C_d = 1, R_d = 6, emp = 0) design.Options.Highway3 <- list(desing.Op3.Highway <- "Highway") design.Options.Urban3 <- list(desing.Op3.Urban <- "Urban") design.Options.Drain3 <- list(desing.Op3.Drain <- "Drain") maintenance.H3 <- dist.Events(lifetime, events.Op3.Highway, start.Date, "Highway") maintenance.U3 <- dist.Events(lifetime, events.Op3.Urban, start.Date, "Urban") maintenance.D3 <- dist.Events(lifetime, events.Op3.Drain, start.Date, "Drainage") integrated <- combine.lifeTimelines(maintenance.H3, maintenance.U3, maintenance.D3, duration.ev3.Highway, duration.ev3.Urban, duration.ev3.Drain) cat("Total system downtime in days over", lifetime, "years:\n") sum(integrated$duration) par(mfrow = c(3, 1)) plot.timeline(maintenance.H3$frequency, maintenance.H3$event, start.Date, "Highway") plot.timeline(maintenance.U3$frequency, maintenance.U3$event, start.Date, "Urban") plot.timeline(maintenance.D3$frequency, maintenance.D3$event, start.Date, "Drainage") # ----------------------------- # Fourth Maintenance Strategy # ----------------------------- lifetime <- 30 start.Date <- "2025-01-01" events.Op4.Highway = c(M_h = 6, R_h = 18, O_h = 30, END = lifetime) duration.ev4.Highway <- c(M_h = 2, R_h = 8, O_h = 20) events.Op4.Urban <- c(M_u = 6, R_u = 18, U_u = 30, END = lifetime) duration.ev4.Urban <- c(M_u = 2, R_u = 6, U_u = 15) events.Op4.Drain <- c(C_d = 6, R_d = 30, emp = 30 , END = lifetime) duration.ev4.Drain <- c(C_d = 1, R_d = 6, emp = 0) design.Options.Highway4 <- list(desing.Op4.Highway <- "Highway") design.Options.Urban4 <- list(desing.Op4.Urban <- "Urban") design.Options.Drain4 <- list(desing.Op4.Drain <- "Drain") maintenance.H4 <- dist.Events(lifetime, events.Op4.Highway, start.Date, "Highway") maintenance.U4 <- dist.Events(lifetime, events.Op4.Urban, start.Date, "Urban") maintenance.D4 <- dist.Events(lifetime, events.Op4.Drain, start.Date, "Drainage") integrated <- combine.lifeTimelines(maintenance.H4, maintenance.U4, maintenance.D4, duration.ev4.Highway, duration.ev4.Urban, duration.ev4.Drain) cat("Total system downtime in days over", lifetime, "years:\n") sum(integrated$duration) par(mfrow = c(3, 1)) plot.timeline(maintenance.H4$frequency, maintenance.H4$event, start.Date, "Highway") plot.timeline(maintenance.U4$frequency, maintenance.U4$event, start.Date, "Urban") plot.timeline(maintenance.D4$frequency, maintenance.D4$event, start.Date, "Drainage") # ----------------------------- # Scenario 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.Op1.Highway) 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.Urban) dur.ev2 <- ev2/2 for (k in 1:dim(events1)[1]) { ev3 <- unlist(events3[k, ]) dist.3 <- dist.Events(lifetime, ev3, start.Date, design.Options.Drain) dur.ev3 <- ev3/2 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)) } # Generate design space n.grid <- 2 events.grid.Highway <- expand.grid(M_h = sample(seq(3,8, by = 1), n.grid), R_h = sample(seq(10,25,1), n.grid), O_h = sample(seq(20,40,1), n.grid)) events.grid.Urban <- expand.grid(M_u = sample(seq(2,9, by = 1), n.grid), R_u = sample(seq(10,25,1), n.grid), U_u = sample(seq(20,40,1), n.grid)) events.grid.Drain <- expand.grid(C_d = sample(seq(4,9, by = 1), n.grid), R_d = sample(seq(20,40,1), n.grid), emp = sample(seq(30,30,1), n.grid)) # Explore design space response.space <- design.explore(events.grid.Highway, events.grid.Urban, events.grid.Drain) # Define preferences and select best options p <- low(dur) * high(dist.inter) sky <- psel(response.space, p) pareto2 <- psel(response.space, p, top = nrow(response.space)) # Visualize results par(mfrow = c(1, 1)) 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 with ggplot instead 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) p <- high(dur) * low(dist.inter) show_front(p) # ------------------------------- # LCA part # ------------------------------- library(dplyr) LCI.materials <- read.csv("LCImaterialss.csv", check.names = FALSE) # ---- Helpers: split UN into slab vs asphalt rows, HY into asphalt vs gravel ---- is_asphalt_row <- function(mat_name) grepl("Asphalt", mat_name, ignore.case = TRUE) is_gravel_row <- function(mat_name) grepl("^Gravel$", mat_name, ignore.case = TRUE) | grepl("Gravel", mat_name, ignore.case = TRUE) # ---- Highway: HY scope, asphalt layer + gravel base layer ---- # NEW INPUT: gravel.tk (base thickness) so gravel has its own volume LCA.highway <- function(length, width, asphalt.tk, gravel.tk, materials, dist.event) { asphalt.Q <- length * width * asphalt.tk # m3 gravel.Q <- length * width * gravel.tk # m3 materials.split <- split(materials, materials$scope) if (is.null(materials.split$HY)) stop("LCI is missing scope 'HY'") HY <- materials.split$HY # interventions for highway asphalt system (same logic as before) s <- summary(as.factor(dist.event$event)) s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors = FALSE, row.names = 1:length(s)) interventions <- s2$DC + 0.25 * (if(!is.null(s2$M.h)) s2$M.h else 0) + 0.55 * (if(!is.null(s2$R.h)) s2$R.h else 0) + 1.00 * (if(!is.null(s2$O.h)) s2$O.h else 0) # Split HY rows into gravel vs asphalt-like rows using material names HY_gravel <- HY %>% filter(is_gravel_row(material)) HY_asphalt <- HY %>% filter(!is_gravel_row(material)) if (nrow(HY_gravel) == 0) warning("HY scope has no 'Gravel' row -> gravel layer will be 0 in LCA.") if (nrow(HY_asphalt) == 0) warning("HY scope has no asphalt rows -> asphalt layer will be 0 in LCA.") # Apply different system volumes HY_gravel <- HY_gravel %>% mutate(sys.Q = gravel.Q, interventions = interventions) HY_asphalt <- HY_asphalt %>% mutate(sys.Q = asphalt.Q, interventions = interventions) LCA.matrix <- bind_rows(HY_gravel, HY_asphalt) %>% mutate( TotalMaterials.Q = quantities * sys.Q / 1000, # ton (if quantities ~ kg/m3) 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) } # ---- Urban: UN scope, slab (RC-like) + asphalt layer ---- LCA.urban <- function(length, width, slab.depth, asphalt.tk, materials, dist.event) { slab.Q <- length * width * slab.depth # m3 asphalt.Q <- length * width * asphalt.tk # m3 materials.split <- split(materials, materials$scope) if (is.null(materials.split$UN)) stop("LCI is missing scope 'UN'") UN <- materials.split$UN # interventions for slab (same as your urban logic) s <- summary(as.factor(dist.event$event)) s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors = FALSE, row.names = 1:length(s)) interventions.slab <- s2$DC + 0.25 * (if(!is.null(s2$M.u)) s2$M.u else 0) + 0.55 * (if(!is.null(s2$R.u)) s2$R.u else 0) + 1.00 * (if(!is.null(s2$U.u)) s2$U.u else 0) # asphalt in urban: keep your assumption (affected by all events) asph.interventions <- length(dist.event$event) - 1 # Split UN rows into asphalt vs slab rows using material name UN_asphalt <- UN %>% filter(is_asphalt_row(material)) UN_slab <- UN %>% filter(!is_asphalt_row(material)) if (nrow(UN_asphalt) == 0) warning("UN scope has no asphalt rows (material contains 'Asphalt') -> asphalt layer will be 0.") if (nrow(UN_slab) == 0) warning("UN scope has no slab rows -> slab layer will be 0.") UN_slab <- UN_slab %>% mutate(sys.Q = slab.Q, interventions = interventions.slab) UN_asphalt <- UN_asphalt %>% mutate(sys.Q = asphalt.Q, interventions = asph.interventions) LCA.matrix <- bind_rows(UN_slab, UN_asphalt) %>% mutate( TotalMaterials.Q = quantities * sys.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) } # ---- Drainage: DE scope (PVC etc.), kg-based ---- LCA.drain <- function(PVC.kg, materials, dist.event) { materials.split <- split(materials, materials$scope) if (is.null(materials.split$DE)) stop("LCI is missing scope 'DE'") DE <- materials.split$DE s <- summary(as.factor(dist.event$event)) s2 <- as.data.frame(rbind(Freq = s), stringsAsFactors = FALSE, row.names = 1:length(s)) interventions <- s2$DC + 0.20 * (if(!is.null(s2$C.d)) s2$C.d else 0) + 1.00 * (if(!is.null(s2$R.d)) s2$R.d else 0) LCA.matrix <- DE %>% mutate( sys.Q = PVC.kg, interventions = interventions, TotalMaterials.Q = quantities * sys.Q, # kg-based (quantities as fraction or kg/kg) 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) } # ------------------------------- # Example geometry # ------------------------------- h.length <- 1000 #1 km h.width <- 21 #6 lane h.asph.tk <- 0.12 h.gravel.tk <- 0.20 u.length <- 1000 #1 km u.width <- 14 #4 lane u.slab.depth <- 0.25 u.asph.tk <- 0.10 PVC.kg <- 25000 Option.H <- LCA.highway(h.length, h.width, h.asph.tk, h.gravel.tk, LCI.materials, maintenance.H4) Option.U <- LCA.urban(u.length, u.width, u.slab.depth, u.asph.tk, LCI.materials, maintenance.U4) Option.D <- LCA.drain(PVC.kg, LCI.materials, maintenance.D4) integrated.Design <- data.frame( Energy = Option.H$Energy + Option.U$Energy + Option.D$Energy, CO2 = Option.H$CO2 + Option.U$CO2 + Option.D$CO2, NOx = Option.H$NOx + Option.U$NOx + Option.D$NOx, SO2 = Option.H$SO2 + Option.U$SO2 + Option.D$SO2 ) Energy.costs <- 0.128 CO2.unitcost <- 26 NOx.unitCost <- 42 SO2.unitCosts <- 85 gn <- integrated.Design %>% mutate(Costs = (Energy * Energy.costs + CO2*CO2.unitcost + NOx*NOx.unitCost + SO2*SO2.unitCosts)/1e9) print(Option.H); print(Option.U); print(Option.D) print(gn) # =============================== # Multi-Objective Optimization # =============================== fitness <- function(x) { # define the output vector (same as tutorial) z <- numeric(7) # keep integer decision variables (same as tutorial) x <- round(x, 0) # --------------------------- # A) Decision variables (PROJECT-SPECIFIC, STRUCTURE SAME) # --------------------------- # Highway: SDO, M, PR # Urban: M.u, PR.u, SR # Drainage: C, Rd, Repl # + width (like tutorial has b.width) y <- expand.grid( M_h = x[1], R_h = x[2], O_h = x[3], M_u = x[4], R_u = x[5], U_u = x[6], C_d = x[7], R_d = x[8], width = x[9] ) # --------------------------- # B) Durations derived from decision variables (SAME LOGIC: ev/2) # --------------------------- dur.ev.H <- unlist(y[c("M_h","R_h","O_h")] / 2); names(dur.ev.H) <- c("M_h","R_h","O_h") dur.ev.U <- unlist(y[c("M_u","R_u","U_u")] / 2); names(dur.ev.U) <- c("M_u","R_u","U_u") dur.ev.D <- unlist(y[c("C_d","R_d")] / 2); names(dur.ev.D) <- c("C_d","R_d") # --------------------------- # C) Maintenance schedules (SAME STRUCTURE: dist.Events) # --------------------------- evH <- c(M_h = y$M_h, R_h = y$R_h, O_h = y$O_h, END = lifetime) evU <- c(M_u = y$M_u, R_u = y$R_u, U_u = y$U_u, END = lifetime) evD <- c(C_d = y$C_d, R_d = y$R_d, END = lifetime) dist.H <- dist.Events(lifetime, evH, start.Date, "Highway") dist.U <- dist.Events(lifetime, evU, start.Date, "Urban") dist.D <- dist.Events(lifetime, evD, start.Date, "Drainage") results <- combine.lifeTimelines(dist.H, dist.U, dist.D, dur.ev.H, dur.ev.U, dur.ev.D) z[1] <- sum(results[["duration"]]) z[2] <- -min(abs( results$frequency[1:(length(results$frequency) - 1)] - results$frequency[2:length(results$frequency)] )) # --------------------------- # E) LCA calls (SAME PLACE IN FLOW) # --------------------------- # LCA calls (UPDATED highway signature) width <- y$width h.length <- 1000 h.asph.tk <- 0.12 h.gravel.tk <- 0.20 # u.length <- 1000 u.slab.depth <- 0.25 u.asph.tk <- 0.10 PVC.kg <- 25000 Energy.costs <- 0.128 CO2.unitcost <- 26 NOx.unitCost <- 42 SO2.unitCosts <- 85 prod.H <- LCA.highway(h.length, width, h.asph.tk, h.gravel.tk, LCI.materials, dist.H) prod.U <- LCA.urban(u.length, width, u.slab.depth, u.asph.tk, LCI.materials, dist.U) prod.D <- LCA.drain(PVC.kg, LCI.materials, dist.D) integrated.system <- as.data.frame(list( Energy = prod.H$Energy + prod.U$Energy + prod.D$Energy, CO2 = prod.H$CO2 + prod.U$CO2 + prod.D$CO2, NOx = prod.H$NOx + prod.U$NOx + prod.D$NOx, SO2 = prod.H$SO2 + prod.U$SO2 + prod.D$SO2 )) integrated.system <- mutate(integrated.system, Costs = (Energy*Energy.costs + CO2*CO2.unitcost + NOx*NOx.unitCost + SO2*SO2.unitCosts)/10^9 ) # same objective mapping as tutorial 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) } # =============================== # NSGA-II call (SAME STRUCTURE) # =============================== r2 <- nsga2( fitness, idim = 9, odim = 7, generations = 10, popsize = 100, lower.bounds = c(3, 10, 20, 2, 10, 20, 4, 20, 14), upper.bounds = c(8, 25, 40, 9, 25, 40, 9, 40, 21) ) r2Results <- as.data.frame(r2$value) outNames <- c("duration", "interv.dist", "energy", "co2", "nox", "so2", "cost") colnames(r2Results) <- outNames pareto3 <- as.data.frame(paretoFront(r2)) colnames(pareto3) <- outNames ggplot(r2Results, aes(x = duration, y = cost)) + geom_point(shape = 21) + geom_point(data = pareto3, size = 3, color="red") + geom_line(data = pareto3, color="blue") # =============================== # PARCOORD # =============================== input.params <- round(r2$par, 0) all.results <- cbind(input.params, r2Results, r2$pareto.optimal) colnames(all.results) <- c( "M_h","R_h","O_h","M_u","R_u","U_u","C_d","R_d","width", "duration","interv.dist","energy","co2","nox","so2","cost", "pareto" ) dev.new(width = 14, height = 5) par(mfrow = c(1, 1)) par(mar = c(5, 4, 4, 1)) # bottom, left, top, right par(xpd = NA) 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))