Merge pull request #67 from MariannaR/tobemerged

Request to merge the new transport module EDGE-T

scripts/iterative/EDGE_transport.R

+require(data.table)
+require(gdx)
+require(gdxdt)
+require(edgeTrpLib)
+require(rmndt)
+## use cached input data for speed purpose
+require(moinput)
+setConfig(forcecache=T)
+
+
+mapspath <- function(fname){
+    file.path("../../modules/35_transport/edge_esm/input", fname)
+}
+
+datapath <- function(fname){
+    file.path("input_EDGE", fname)
+}
+
+REMINDpath <- function(fname){
+    file.path("../../", fname)
+}
+
+REMINDyears <- c(1990,
+           seq(2005, 2060, by = 5),
+           seq(2070, 2110, by = 10),
+           2130, 2150)
+
+gdx <- "input.gdx"
+if(file.exists("fulldata.gdx"))
+    gdx <- "fulldata.gdx"
+
+load("config.Rdata")
+scenario <- cfg$gms$cm_GDPscen
+EDGE_scenario <- cfg$gms$cm_EDGEtr_scen
+
+EDGEscenarios <- fread("../../modules/35_transport/edge_esm/input/EDGEscenario_description.csv")[scenario_name == EDGE_scenario]
+
+merge_traccs <<- EDGEscenarios[options == "merge_traccs", switch]
+addvintages <<- EDGEscenarios[options == "addvintages", switch]
+inconvenience <<- EDGEscenarios[options == "inconvenience", switch]
+selfmarket_taxes <<- EDGEscenarios[options == "selfmarket_taxes", switch]
+selfmarket_policypush <<- EDGEscenarios[options == "selfmarket_policypush", switch]
+selfmarket_acceptancy <<- EDGEscenarios[options == "selfmarket_acceptancy", switch]
+
+if (EDGE_scenario == "Conservative_liquids") {
+  techswitch <<- "Liquids"
+} else if (EDGE_scenario %in% c("Electricity_push", "Smart_lifestyles_Electricity_push")) {
+  techswitch <<- "BEV"
+} else if (EDGE_scenario == "Hydrogen_push") {
+  techswitch <<- "FCEV"
+} else {
+  print("You selected a not allowed scenario. Scenarios allowed are: Conservative_liquids, Hydrogen_push, Electricity_push, Smart_lifestyles_Electricity_push")
+  exit()
+}
+
+endogeff <<- EDGEscenarios[options== "endogeff", switch]
+enhancedtech <<- EDGEscenarios[options== "enhancedtech", switch]
+rebates_febates <<- EDGEscenarios[options== "rebates_febates", switch] ##NB THEY ARE ONLY IN PSI! ONLY WORKING IN EUROPE
+savetmpinput <<- FALSE
+smartlifestyle <<- EDGEscenarios[options== "smartlifestyle", switch]
+
+
+
+REMIND2ISO_MAPPING <- fread(REMINDpath(cfg$regionmapping))[, .(iso = CountryCode, region = RegionCode)]
+EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv"))
+
+
+## input data loading
+input_path = paste0("../../modules/35_transport/edge_esm/input/")
+
+inputdata = createRDS(input_path, SSP_scenario = scenario, EDGE_scenario = EDGE_scenario)
+vot_data = inputdata$vot_data
+sw_data = inputdata$sw_data
+inco_data = inputdata$inco_data
+logit_params = inputdata$logit_params
+int_dat = inputdata$int_dat
+nonfuel_costs = inputdata$nonfuel_costs
+price_nonmot = inputdata$price_nonmot
+
+## add learning optional
+setlearning = TRUE
+## add optional vintages
+addvintages = TRUE
+## optional average of prices
+average_prices = FALSE
+## inconvenience costs instead of preference factors
+inconvenience = TRUE
+
+if (setlearning | addvintages){
+  ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, REMINDyears)
+  ## select from total demand only the passenger sm
+  ES_demand = ES_demand[sector == "trn_pass",]
+}
+
+
+if (setlearning & file.exists("demand_previousiter.RDS")) {
+  ## load previous iteration number of cars
+  demand_BEVtmp = readRDS("demand_BEV.RDS")
+  ## load previous iteration demand
+  ES_demandpr = readRDS("demand_previousiter.RDS")
+  ## calculate non fuel costs and
+  nonfuel_costs = applylearning(gdx,REMINDmapping,EDGE2teESmap, demand_BEVtmp, ES_demandpr)
+  saveRDS(nonfuel_costs, "nonfuel_costs_learning.RDS")
+}
+
+## load price
+REMIND_prices <- merge_prices(
+    gdx = gdx,
+    REMINDmapping = REMIND2ISO_MAPPING,
+    REMINDyears = REMINDyears,
+    intensity_data = int_dat,
+    nonfuel_costs = nonfuel_costs)
+
+## save prices
+## read last iteration count
+keys <- c("iso", "year", "technology", "vehicle_type")
+setkeyv(REMIND_prices, keys)
+
+pfile <- "EDGE_transport_prices.rds"
+iter <- as.vector(gdxrrw::rgdx(gdx, list(name="o_iterationNumber"))$val)
+
+REMIND_prices[, iternum := iter]
+
+## save REMIND prices (before dampening)
+saveRDS(REMIND_prices, paste0("REMINDprices", iter, ".RDS"))
+
+
+if(average_prices){
+
+  if(max(unique(REMIND_prices$iternum)) >= 20 & max(unique(REMIND_prices$iternum)) <= 30){
+    old_prices <- readRDS(pfile)
+    all_prices <- rbind(old_prices, REMIND_prices)
+    setkeyv(all_prices, keys)
+    ## apply moving avg
+    REMIND_prices <- REMIND_prices[
+      all_prices[iternum >= 20, mean(tot_price), by=keys], tot_price := V1]
+    all_prices <- rbind(old_prices, REMIND_prices)
+  }else{
+    all_prices <- REMIND_prices
+  }
+  saveRDS(all_prices, pfile)
+
+  ## save REMIND prices (after dampening)
+  saveRDS(REMIND_prices,paste0("REMINDpricesDampened", iter, ".RDS"))
+
+}
+
+REMIND_prices[, "iternum" := NULL]
+
+## calculates logit
+if (inconvenience) {
+  years=copy(REMINDyears)
+
+  logit_data <- calculate_logit_inconv_endog(
+    prices= REMIND_prices[tot_price > 0],
+    vot_data = vot_data,
+    inco_data = inco_data,
+    logit_params = logit_params,
+    intensity_data = int_dat,
+    price_nonmot = price_nonmot)
+
+} else{
+
+  logit_data <- calculate_logit(
+    REMIND_prices[tot_price > 0],
+    REMIND2ISO_MAPPING,
+    vot_data = vot_data,
+    sw_data = sw_data,
+    logit_params = logit_params,
+    intensity_data = int_dat,
+    price_nonmot = price_nonmot)
+
+}
+shares <- logit_data[["share_list"]] ## shares of alternatives for each level of the logit function
+## shares$VS1_shares=shares$VS1_shares[,-c("sector","subsector_L2","subsector_L3")]
+
+mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level
+prices <- logit_data[["prices_list"]] ## prices at each level of the logit function, 1990USD/pkm
+
+if(addvintages){
+  ## calculate vintages (new shares, prices, intensity)
+  vintages = calcVint(shares = shares,
+                      totdem_regr = ES_demand,
+                      prices = prices,
+                      mj_km_data = mj_km_data,
+                      years = REMINDyears)
+
+  shares$FV_shares = vintages[["shares"]]$FV_shares
+  prices = vintages[["prices"]]
+  mj_km_data = vintages[["mj_km_data"]]
+}
+
+## use logit to calculate shares and intensities (on tech level)
+EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv"))
+
+shares_intensity_demand <- shares_intensity_and_demand(
+    logit_shares=shares,
+    MJ_km_base=mj_km_data,
+    EDGE2CESmap=EDGE2CESmap,
+    REMINDyears=REMINDyears,
+    scenario=scenario,
+    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING)
+
+demByTech <- shares_intensity_demand[["demand"]] ##in [-]
+intensity <- shares_intensity_demand[["demandI"]] ##in million pkm/EJ
+norm_demand <- shares_intensity_demand$demandF_plot_pkm ## total demand is 1, required for costs
+
+if (setlearning) {
+  demand_BEV=calc_num_vehicles( norm_dem_BEV = norm_demand[technology == "BEV" & ## battery vehicles
+                                                           subsector_L1 == "trn_pass_road_LDV_4W", ## only 4wheelers
+                                                           c("iso", "year", "sector", "vehicle_type", "demand_F") ],
+                                ES_demand = ES_demand)
+
+  ## save number of vehicles for next iteration
+  saveRDS(demand_BEV, "demand_BEV.RDS")
+  ## save the demand for next iteration renaming the column
+  setnames(ES_demand, old ="demand", new = "demandpr")
+  saveRDS(ES_demand, "demand_previousiter.RDS")
+}
+
+
+## use logit to calculate costs
+budget <- calculate_capCosts(
+    base_price=prices$base,
+    Fdemand_ES = norm_demand,
+    EDGE2CESmap = EDGE2CESmap,
+    EDGE2teESmap = EDGE2teESmap,
+    REMINDyears = REMINDyears,
+    scenario = scenario,
+    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING)
+
+## full REMIND time range for inputs
+REMINDtall <- c(seq(1900,1985,5),
+                seq(1990, 2060, by = 5),
+                seq(2070, 2110, by = 10),
+                2130, 2150)
+
+## prepare the entries to be saved in the gdx files: intensity, shares, non_fuel_price. Final entries: intensity in [trillionkm/Twa], capcost in [2005USD/trillionpkm], shares in [-]
+finalInputs <- prepare4REMIND(
+    demByTech = demByTech,
+    intensity = intensity,
+    capCost = budget,
+    EDGE2teESmap = EDGE2teESmap,
+    REMINDtall = REMINDtall,
+    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING)
+
+
+
+## add the columns of SSP scenario and EDGE scenario to the output parameters
+for (i in names(finalInputs)) {
+             finalInputs[[i]]$SSP_scenario <- scenario
+             finalInputs[[i]]$EDGE_scenario <- EDGE_scenario
+           }
+
+
+## calculate shares
+finalInputs$shFeCes = finalInputs$demByTech[, value := value/sum(value), by = c("tall", "all_regi", "all_in")]
+
+
+## CapCosts
+writegdx.parameter("p35_esCapCost.gdx", finalInputs$capCost, "p35_esCapCost",
+                   valcol="value", uelcols=c("tall", "all_regi", "SSP_scenario", "EDGE_scenario", "all_teEs"))
+
+## Intensities
+writegdx.parameter("p35_fe2es.gdx", finalInputs$intensity, "p35_fe2es",
+                   valcol="value", uelcols = c("tall", "all_regi", "SSP_scenario", "EDGE_scenario", "all_teEs"))
+
+## Shares: demand can represent the shares since it is normalized
+writegdx.parameter("p35_shFeCes.gdx", finalInputs$shFeCes, "p35_shFeCes",
+                   valcol="value",
+                   uelcols = c("tall", "all_regi", "SSP_scenario", "EDGE_scenario", "all_enty", "all_in", "all_teEs"))

scripts/output/single/EDGETransportValidation.R

+require(rmarkdown)
+
+## run EDGE transport validation output if required
+if(cfg$gms$transport == "edge_esm"){
+  md_template <- "EDGETransportReport.Rmd"
+  file.copy(file.path("./scripts/output/single/notebook_templates", md_template), outputdir)
+  rmarkdown::render(path(outputdir, md_template), output_format="pdf_document")
+  if(cfg$gms$c_keep_iteration_gdxes == 1){
+    md_template <- "EDGETransportMultiIterationAnalysis.Rmd"
+    file.copy(file.path("./scripts/output/single/notebook_templates", md_template), outputdir)
+    rmarkdown::render(path(outputdir, md_template), output_format="pdf_document")
+  }
+}

scripts/output/single/notebook_templates/EDGETransportMultiIterationAnalysis.Rmd

+---
+title: "Analysis of the transport module for a run with resolved iteration domain"
+output: html_document
+---
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+require(ggplot2)
+require(lusweave)
+require(rmndt)
+
+
+knitr::opts_chunk$set(fig.width=12, fig.height=12) 
+
+folder <- "./"
+files <- list.files(path = folder, pattern = "fulldata_[0-9]+\\.gdx")
+
+year_toplot = 2050
+iter_toplot = 25
+maxiter = 100
+
+print(paste0("Year: ", year_toplot))
+print(paste0("Iteration: ", iter_toplot))
+
+```
+
+```{r, echo=FALSE}
+addyrs <- function(dt, yrcol="ttot"){
+    dt[, year := as.numeric(get(yrcol))][, (yrcol) := NULL]
+    return(dt)
+}
+
+get_trp_shares <- function(gdx){
+    tes <- readGDX(gdx, "teEs_dyn35")
+    
+    shares_data <- readgdx(gdx, "pm_shFeCes")[all_teEs %in% tes]
+    shares_data <- addyrs(shares_data)
+    
+    setnames(shares_data, "value", "tech_share")
+    return(shares_data)
+}
+
+get_fuel_prices <- function(gdx){
+  REMINDyears <- c(1990,
+                   seq(2005, 2060, by = 5),
+                   seq(2070, 2110, by = 10),
+                   2130, 2150)
+  ## report prices from REMIND gdx in 2005$/MJ
+  
+  tdptwyr2dpgj <- 31.71  #TerraDollar per TWyear to Dollar per GJ
+  CONV_2005USD_1990USD <- 0.67
+  startyear <- 2020
+  
+  ## load entries from the gdx
+  fety <- readGDX(gdx, c("entyFe", "fety"), format = "first_found")
+  
+  budget.m <- readGDX(gdx, name = "qm_budget", types = "equations", field = "m",
+                      format = "first_found")[, REMINDyears[REMINDyears >= startyear],]  # Alternative: calcPrice
+  
+  interpolate_first_timesteps <- function(obj){
+    ## interpolate values for 1990, 2005 and 2010
+    obj = time_interpolate(obj, c(1990, seq(2005, startyear, 5)),
+                           integrate_interpolated_years=T,
+                           extrapolation_type = "constant")
+    return(obj)
+  }
+  
+  budget.m <- interpolate_first_timesteps(budget.m)
+  budget.m <- lowpass(budget.m)
+  
+  
+  bal_eq <- "qm_balFeForCesAndEs"
+  febal.m <- readGDX(gdx, name = bal_eq, types = "equations",
+                      field = "m", format = "first_found")[, REMINDyears[REMINDyears>=startyear], fety]
+    
+  if(any(febal.m > 0)){
+    sprintf("Found positive marginals on %s. We correct this, but the issue should be adressed.", bal_eq)
+    febal.m[febal.m > 0] <- -1e-10
+  }
+    
+  febal.m <- interpolate_first_timesteps(febal.m)
+    
+  ## in some regions and time steps, 0 final energy demand for an entry could give problems
+  tmp <- magclass::setNames(lowpass(lowpass(febal.m[, , "fegat"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fegat")
+    
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "feelt"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "feelt"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "feh2t"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "feh2t"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "fedie"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fedie"))
+    
+  tmp <- mbind(tmp, magclass::setNames(lowpass(lowpass(febal.m[, , "fepet"]))/(budget.m + 1e-10) * tdptwyr2dpgj, "fepet"))
+                              
+  tmp <- magpie2dt(tmp, regioncol = "all_regi", yearcol = "year", datacols = "all_enty")
+
+  setnames(tmp, old = "value", new = "fe_price")
+  return(tmp)
+}
+
+get_demand <- function(gdx){
+    ces_nodes <- readGDX(gdx, "ppfen_dyn35")
+    demand <- readgdx(gdx, "vm_cesIO")[all_in %in% ces_nodes]
+    demand <- addyrs(demand, "tall")
+    setnames(demand, "value", "CES_demand")
+    return(demand)
+}
+```
+
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+sorted_files <- paste0(folder, "fulldata_", 1:length(files), ".gdx")
+
+shares <- lapply(sorted_files, get_trp_shares)
+prices <- lapply(sorted_files, get_fuel_prices)
+demand <- lapply(sorted_files, get_demand)
+
+library(stringr)
+# lapply(files, function(fname){str_extract(fname, "[0-9]+")})
+
+for(fname in files){
+    idx <- as.numeric(str_extract(fname, "[0-9]+"))
+    shares[[idx]][, iter := idx]
+    prices[[idx]][, iter := idx]
+    demand[[idx]][, iter := idx]
+}
+
+shares <- rbindlist(shares)
+prices <- rbindlist(prices)
+demand <- rbindlist(demand)
+
+## we are interested only in years up to 2100
+shares = shares[year <= 2100]
+prices = prices[year <= 2100]
+demand = demand[year <= 2100]
+```
+
+##Prices and shares in the iteration domain
+
+```{r, echo=FALSE}
+toplot <- prices[year == year_toplot  & iter < maxiter]
+
+ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2, scales = "free" )
+# +
+#   ylim(c(0,50))
+```
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-5,year_toplot,year_toplot+5)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-5, ", ", year_toplot, ", ", year_toplot+5)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-5,year_toplot)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-5, ", ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year %in% c(year_toplot-10,year_toplot-5,year_toplot)]
+# toplot <- toplot[,.(fe_price=mean(fe_price)), by=c("all_regi","all_enty","iter")]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Iteration Domain, average of ", year_toplot-10, ", ", year_toplot-5, ", ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2, scales = "free")
+
+
+```
+
+
+```{r, echo=FALSE}
+# toplot <- prices[year == year_toplot & all_regi == "EUR"]
+# 
+# ggplot(aes(iter, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line(size = 2) +
+#     labs(x = "iteration #", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("European Fuel Prices in the Iteration Domain, in ", year_toplot)) +
+#    # facet_wrap(~ all_regi, ncol=2)+
+#   theme_light()+
+#   scale_color_brewer(palette = "Set1") 
+```
+
+
+```{r, echo =FALSE}
+# toplot <- shares[year == year_toplot & iter < maxiter]
+# 
+# ggplot(aes(iter, tech_share, colour=all_teEs, shares=all_teEs), data=toplot) +
+#     geom_line() +
+#     labs(x = "iteration #", y = "Tech Share", color = "LDV Tech Type", title=paste0("Tech Shares in the Iteration Domain, in ", year_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2)
+```
+
+```{r, echo=FALSE}
+
+toplot <- demand[year == year_toplot & all_in %in% c("entrp_pass_sm", "entrp_pass_lo")  & iter < maxiter]
+
+ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "CES node aggr. pass. demand (trillion pkm)", color = "CES category", title=paste0("Aggregated Passenger Tranport Demand in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2)
+```
+
+
+```{r, echo=FALSE}
+
+# toplot <- demand[year == year_toplot & all_in %in% c("entrp_pass_sm", "entrp_pass_lo") & all_regi == "EUR"]
+# 
+# ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+#     geom_line(size = 2) +
+#     labs(x = "iteration #", y = "CES node aggr. pass. demand (trillion pkm)", color = "CES category", title=paste0("Europe: aggregated Passenger Tranport Demand in the Iteration Domain, in ", year_toplot)) +
+#    # facet_wrap(~ all_regi, ncol=2)+
+#   theme_light()
+```
+
+
+```{r, echo=FALSE}
+toplot <- demand[year == year_toplot & all_in %in% c("entrp_frgt_sm", "entrp_frgt_lo")  & iter < maxiter]
+ggplot(aes(iter, CES_demand, colour=all_in, shares=all_in), data=toplot) +
+    geom_line() +
+    labs(x = "iteration #", y = "CES node aggr. freight demand (trillion tkm)", color = "CES category", title= paste0("Aggregated Freight Transport Demand in the Iteration Domain, in ", year_toplot)) +
+    facet_wrap(~ all_regi, ncol=2)
+```
+
+##Prices in the time domain
+
+```{r, echo=FALSE}
+# toplot <- prices[iter == iter_toplot]
+# 
+# ggplot(aes(year, fe_price, group=all_enty, color=all_enty), data=toplot) +
+#     geom_line() +
+#     labs(x = "Year", y = "US$ 2015/GJ", color = "Fuel Type", title=paste0("Fuel Prices in the Time Domain, in iteration ", iter_toplot)) +
+#     facet_wrap(~ all_regi, ncol=2)
+# +
+#   ylim(c(0,50))
+```
+
+```{r, echo=FALSE}
+
+enties = c("fegat", "fepet", "fedie", "feelt", "feh2t")
+
+for (enty in enties) {
+  toplot <- prices[all_enty == enty  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+
+  plot=ggplot(aes(year, fe_price, group=interaction(iter, all_enty), color=iter), data=toplot) +
+    geom_line() +
+    labs(x = "Year", y = "US$ 2015/GJ", color = "Iteration number", title=paste0("Fuel Prices in the Time Domain, across iterations, for ", enty)) +
+    facet_wrap(~ all_regi, scales = "free", ncol = 2)
+
+  print(plot)
+}
+
+```
+
+```{r, echo =FALSE}
+
+
+enty_vals = unique(prices$all_enty)
+toplot <- prices[year>=2005 & iter < maxiter]
+toplot[, iter := as.numeric(iter)]
+toplot[, min := min(fe_price), by = c("all_regi", "all_enty", "year")]
+toplot[, max := max(fe_price), by = c("all_regi", "all_enty", "year")]
+
+plot= ggplot()+
+  geom_ribbon(data = toplot, aes(x=year, group = all_enty, fill = all_enty, ymin = min, ymax = max), alpha = 0.5)+
+  facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+  theme_minimal() +
+  facet_wrap(~ all_regi,ncol = 2) +
+  geom_line(data = toplot[iter == max(iter)], aes(year, fe_price, group=all_enty, color = all_enty), linetype = "dashed")+
+    labs(x = "Year", y = "US$ 2015/GJ", fill = "Fuel Type", color = "Final iteration", title=paste0("Fuel Prices in the Time Domain, across iterations"))
+
+
+print(plot)
+```
+
+
+## Shares in the time domain
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- shares[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+
+  plot=ggplot() +
+    geom_line(aes(year, tech_share, group=interaction(iter, all_teEs), color=all_teEs, alpha = iter, linetype ="Intermediate Iter." ), data=toplot) +
+    geom_line(data = toplot[iter == max(iter)], aes(year, tech_share, group=all_teEs, linetype ="Last Iter.", color ="Last Iter."), color="black")+
+    scale_linetype_manual("Iteration Status",values=c("Intermediate Iter." = 1, "Last Iter."=2))+
+    labs(x = "Year", y = "[-]", color = "Fuel Type", title=paste0("Tech shares in the Time Domain, across iterations, for ", in_val)) +
+    facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+    theme_minimal() +
+    facet_wrap(~ all_regi,ncol = 2)
+
+  print(plot)
+  
+}
+```
+
+
+
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- shares[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+  toplot[, min := min(tech_share), by = c("all_regi", "all_teEs", "year")]
+  toplot[, max := max(tech_share), by = c("all_regi", "all_teEs", "year")]
+
+ plot= ggplot()+
+  geom_ribbon(data = toplot, aes(x=year, group = all_teEs, fill = all_teEs, ymin = min, ymax = max), alpha = 0.5)+
+  facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+  theme_minimal() +
+  facet_wrap(~ all_regi,ncol = 2) +
+  geom_line(data = toplot[iter == max(iter)], aes(year, tech_share, group=all_teEs, color = all_teEs), linetype = "dashed")+
+    labs(x = "Year", y = "[-]", color = "Fuel Type", title=paste0("Tech shares in the Time Domain, across iterations, for ", in_val))
+
+
+  print(plot)
+  
+}
+```
+
+
+
+```{r, echo =FALSE}
+
+
+in_vals = unique(shares$all_in)
+
+for (in_val in in_vals) {
+  toplot <- demand[year>=2005 & all_in == in_val  & iter < maxiter]
+  toplot[, iter := as.numeric(iter)]
+  if (in_val %in% c("entrp_pass_sm","entrp_pass_lo")) {
+    ylabel = "[billion pkm]"
+  } else {
+    ylabel = "[billion tkm]"
+  }
+
+ plot= ggplot()+
+   geom_line(data = toplot, aes(x=year, y = CES_demand, group = iter, color = all_in, alpha = iter, linetype ="Intermediate Iter." ))+
+   geom_line(data = toplot[iter == max(iter)], aes(year, CES_demand, group=all_in, linetype ="Last Iter.", color ="Last Iter."), color="black")+
+   facet_wrap(~ all_regi, scales = "free", ncol = 2)+
+   theme_minimal() +
+   facet_wrap(~ all_regi,ncol = 2) +
+   scale_linetype_manual("Iteration Status",values=c("Intermediate Iter." = 1, "Last Iter."=2))+
+    labs(x = "Year", y = ylabel, color = "Transport Type", title=paste0("CES node value in Time Domain, across iterations, for ", in_val))
+
+
+  print(plot)
+  
+}
+```
+## A quality measure of convergence
+
+For starters, let us implement something rather simple:
+
+$$Q{r,s}=\sum_{t}\frac{\sum_{i=m}^N [s_{t,i}-Avg_i(s_{t,i})]^2}{Avg^2_i(s_{t,i})}$$
+
+
+for a variable $s$ over times $t$ with iteration index $i$ (first m iterations being ignored) and where $Avg_i$ denotes the average over the iteration index dimension. The region index $r$ is omitted on the right hand side.
+
+
+```{r, echo =FALSE}
+quality <- function(arr, m=3){
+    # we ignore the first m iterations
+    sum((arr[-(1:m)] - ave(arr[-(1:m)]))^2/ave(arr[-(1:m)])^2)
+}
+
+shares[, quality := quality(.SD$tech_share), by=c("all_regi", "all_teEs", "year")]
+top <- head(shares[order(-quality), max(quality), by=c("all_regi", "all_teEs", "year")], 15)
+top
+
+```
+
+```{r, echo=FALSE}
+ggplot(aes(iter, tech_share), data=shares[all_regi == top[1]$all_regi & all_teEs == top[1]$all_teEs & year == top[1]$year]) +
+    geom_line() +
+    labs(x = "iteration #", y = "Tech Share", color = "year", title=paste(top[1]$all_te, "Shares in", top[1]$all_regi, "in the Iteration Domain, in ", top[1]$year))
+```
+