From 82fd14a7000cb8a730e0b4c73f6560e87553bdc4 Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Mon, 2 Oct 2023 14:35:27 +0200 Subject: [PATCH 1/8] fix LCOE reporting script --- R/reportLCOE.R | 45 +++++++++++++++++++++------------------------ 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index a6692dd9..5a9715c3 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -64,11 +64,20 @@ reportLCOE <- function(gdx, output.type = "both"){ LCOE.out <- NULL - # read in general data + # read in general data (needed for average and marginal LCOE calculation) s_twa2mwh <- readGDX(gdx,c("sm_TWa_2_MWh","s_TWa_2_MWh","s_twa2mwh"),format="first_found") ttot <- as.numeric(readGDX(gdx,"ttot")) ttot_before2005 <- paste0("y",ttot[which(ttot <= 2000)]) ttot_from2005 <- paste0("y",ttot[which(ttot >= 2005)]) + te <- readGDX(gdx,"te") + te <- te[!te %in% c("lng_liq","gas_pipe", "lng_gas", "lng_ves", "coal_ves", "pipe_gas", "termX_lng", "termM_lng", "vess_lng")] + p_priceCO2 <- readGDX(gdx,name=c("p_priceCO2","pm_priceCO2"),format="first_found") # co2 price + + + ## equations + qm_pebal <- readGDX(gdx,name=c("q_balPe"),field="m",format="first_found") + qm_budget <- readGDX(gdx,name=c("qm_budget"),field="m",format="first_found") + ######################################################## @@ -84,8 +93,6 @@ reportLCOE <- function(gdx, output.type = "both"){ temapse <- readGDX(gdx,"en2se") temapall <- readGDX(gdx,c("en2en","temapall"),format="first_found") teall2rlf <- readGDX(gdx,c("te2rlf","teall2rlf"),format="first_found") - te <- readGDX(gdx,"te") - te <- te[!te %in% c("lng_liq","gas_pipe", "lng_gas", "lng_ves", "coal_ves", "pipe_gas", "termX_lng", "termM_lng", "vess_lng")] te2stor <- readGDX(gdx,"VRE2teStor") te2grid <- readGDX(gdx,"VRE2teGrid") teVRE <- readGDX(gdx,"teVRE") @@ -109,7 +116,6 @@ reportLCOE <- function(gdx, output.type = "both"){ pm_ts <- readGDX(gdx,"pm_ts") pm_data <- readGDX(gdx,"pm_data") pm_emifac <- readGDX(gdx,"pm_emifac", restore_zeros=F) # emission factor per technology - p_priceCO2 <- readGDX(gdx,name=c("p_priceCO2","pm_priceCO2"),format="first_found") # co2 price pm_taxemiMkt <- readGDX(gdx,"pm_taxemiMkt") # regional co2 price pm_eta_conv <- readGDX(gdx,"pm_eta_conv", restore_zeros=F) # efficiency oftechnologies with time-dependent eta pm_dataeta <- readGDX(gdx,"pm_dataeta", restore_zeros=F)# efficiency of technologies with time-independent eta @@ -126,11 +132,7 @@ reportLCOE <- function(gdx, output.type = "both"){ vm_cap <- readGDX(gdx,name=c("vm_cap"),field="l",format="first_found") vm_prodFe <- readGDX(gdx,name=c("vm_prodFe"),field="l",restore_zeros=FALSE,format="first_found") v_emiTeDetail <- readGDX(gdx,name=c("vm_emiTeDetail","v_emiTeDetail"),field="l",restore_zeros=FALSE,format="first_found") - - - ## equations - qm_pebal <- readGDX(gdx,name=c("q_balPe"),field="m",format="first_found") - qm_budget <- readGDX(gdx,name=c("qm_budget"),field="m",format="first_found") + # module-specific # amount of curtailed electricity @@ -142,13 +144,6 @@ reportLCOE <- function(gdx, output.type = "both"){ v32_curt <- 0 } - # # sensitivites (only for investment cost sensitivity runs) - # p_costFac <- readGDX(gdx,name=c("p_costFac"),react = "silent") # sensitivity factor - # if (is.null(p_costFac)) { - # p_costFac <- new.magpie(getRegions(v_directteinv),getYears(v_directteinv),getNames(p_omeg,dim=2), fill=1) - # } - - ####### calculate needed parameters ############### @@ -183,7 +178,7 @@ reportLCOE <- function(gdx, output.type = "both"){ v_investcost[,ttot_before2005,te] * dimSums(vm_deltaCap[teall2rlf][,ttot_before2005,te],dim=3.2), v_directteinv_wadj[,ttot_from2005,te] ) - + ########## calculation of LCOE of standing system ####### ######## (old annuities included) ###################### @@ -223,7 +218,7 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_inv_cost[,t0,a] <- dimSums(pm_ts[,t_id,] * te_inv_annuity[,t_id,a] * p_omeg_h[,t_id,a] ,dim=2) } # a } # t0 - + te_annual_inv_cost_wadj <- new.magpie(getRegions(te_inv_annuity_wadj[,ttot,]),getYears(te_inv_annuity_wadj[,ttot,]),magclass::getNames(te_inv_annuity_wadj[,ttot,])) # loop over ttot for(t0 in ttot){ @@ -239,7 +234,7 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_inv_cost_wadj[,t0,a] <- dimSums(pm_ts[,t_id,] * te_inv_annuity_wadj[,t_id,a] * p_omeg_h[,t_id,a] ,dim=2) } # a } # t0 - + ####### 2. sub-part: fuel cost ################################# # fuel cost = PE price * PE demand of technology @@ -457,7 +452,7 @@ reportLCOE <- function(gdx, output.type = "both"){ # convert to better dimensional format df.lcoe.avg <- as.quitte(LCOE.avg) %>% select(region, period, data, value) %>% - rename(variable = data) %>% + rename(variable = data) %>% replace(is.na(.), 0) @@ -479,11 +474,10 @@ reportLCOE <- function(gdx, output.type = "both"){ df.lcoe.avg <- df.lcoe.avg %>% mutate( unit = "US$2015/MWh") %>% - select(region, period, type, output, tech, sector, unit, cost, value) %>% - as.quitte() + select(region, period, type, output, tech, sector, unit, cost, value) # reconvert to magpie object - LCOE.avg.out <- as.magpie(df.lcoe.avg, spatial=1, temporal=2, datacol=9) + LCOE.avg.out <- suppressWarnings(as.magpie(df.lcoe.avg, spatial=1, temporal=2, datacol=9)) # bind to output file LCOE.out <- mbind(LCOE.out, LCOE.avg.out) } @@ -593,6 +587,7 @@ reportLCOE <- function(gdx, output.type = "both"){ pe2se <- readGDX(gdx,"pe2se") # pe2se technology mappings se2se <- readGDX(gdx,"se2se") # hydrogen <--> electricity technologies se2fe <- readGDX(gdx,"se2fe") # se2fe technology mappings + te <- readGDX(gdx, "te") # all technologies teStor <- readGDX(gdx, "teStor") # storage technologies for VREs teGrid <- readGDX(gdx, "teGrid") # grid technologies for VREs ccs2te <- readGDX(gdx, "ccs2te") # ccsinje technology @@ -781,6 +776,8 @@ reportLCOE <- function(gdx, output.type = "both"){ ### 7. retrieve carbon price + p_priceCO2 <- readGDX(gdx,name=c("p_priceCO2","pm_priceCO2"),format="first_found") # co2 price + df.co2price <- as.quitte(p_priceCO2) %>% select(region, period, value) %>% # where carbon price is NA, it is zero @@ -1429,7 +1426,7 @@ reportLCOE <- function(gdx, output.type = "both"){ left_join(df.CCStax, by = c("region", "period", "tech")) %>% left_join(df.CCSCost, by = c("region", "period", "tech")) %>% left_join(df.flexPriceShare, by = c("region", "period", "tech")) %>% - left_join(df.FEtax, by = c("region", "period", "output")) %>% + left_join(df.FEtax,relationship = "many-to-many", by = c("region", "period", "output")) %>% left_join(df.AddTeInvH2, by = c("region", "period", "tech", "sector")) %>% # filter to only have LCOE technologies filter( tech %in% c(te_LCOE)) From 39a128354348fa1cfb51203146caecd14d7834cf Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Mon, 2 Oct 2023 14:36:24 +0200 Subject: [PATCH 2/8] improve documentation of marginal LCOE calculation in LCOE reporting script --- R/reportLCOE.R | 46 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index 5a9715c3..062d1bce 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -513,7 +513,7 @@ reportLCOE <- function(gdx, output.type = "both"){ OMF <- NULL OMV <- NULL lifetime <- NULL - disc.fac <- NULL + annuity.fac <- NULL CapFac <- NULL fuel.price.weighted.mean <- NULL `Investment Cost` <- NULL @@ -733,13 +733,13 @@ reportLCOE <- function(gdx, output.type = "both"){ # as.magpie() # read lifetime of technology - # calculate annuity factor to annuitize CAPEX and OMF (annuity factor labeled "disc.fac") + # calculate annuity factor to annuitize CAPEX and OMF (annuity factor labeled "annuity.fac") lt <- readGDX(gdx, name="fm_dataglob", restore_zeros = F)[,,"lifetime"][,,te_LCOE_Inv][,,"lifetime"] df.lifetime <- as.quitte(lt) %>% select(all_te, value) %>% rename(tech = all_te, lifetime = value) %>% - mutate( disc.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) + mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) ### note: just take LCOE investment cost formula with fix lifetime, ignoring that in REMIND capacity is drepeciating over time ### actually you would need to divide CAPEX by CapFac * 8760 * p_omeg and then add a discounting. @@ -1059,9 +1059,9 @@ reportLCOE <- function(gdx, output.type = "both"){ # # conversion from tr USD 2005/TW to USD2015/kW # mutate(CAPEX = CAPEX *1.2 * 1e3) %>% # # calculate annuity factor for investment cost - # mutate( disc.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) %>% + # mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) %>% # # investment cost LCOE in USD/MWh - # mutate( `Investment Cost` = CAPEX * disc.fac / (CapFac*8760)*1e3) %>% + # mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% # # OMF cost LCOE in USD/MWh # mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% # # calculate pre-loss VRE LCOE as investment cost + OMF cost LCOE @@ -1185,7 +1185,7 @@ reportLCOE <- function(gdx, output.type = "both"){ # #### calculate LCOE # df.LCOE <- df.LCOE %>% # # investment cost LCOE in USD/MWh - # mutate( `Investment Cost` = CAPEX * disc.fac / (CapFac*8760)*1e3) %>% + # mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% # # OMF cost LCOE in USD/MWh # mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% # mutate( `OMV Cost` = OMV) %>% @@ -1206,7 +1206,7 @@ reportLCOE <- function(gdx, output.type = "both"){ left_join(df.gridfactor, by = "region") %>% filter(tech == "gridwind") %>% rename(gridtech = tech) %>% - mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * disc.fac) %>% + mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * annuity.fac) %>% mutate(`OMF Cost` = CAPEX*1e12*1.2/s_twa2mwh * OMF) %>% mutate(grid.cost = (`Investment Cost` + `OMF Cost`) / grid.factor) %>% full_join(teVRE.grid, by = "gridtech") %>% @@ -1247,7 +1247,7 @@ reportLCOE <- function(gdx, output.type = "both"){ left_join(df.eff, by=c("region"="region","period"="period","teStor"="tech")) %>% left_join(df.CapFac, by=c("region"="region","period"="period","teStor"="tech")) %>% # cost per MWh v32_storloss - mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * disc.fac) %>% + mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * annuity.fac) %>% mutate(`OMF Cost` = CAPEX*1e12*1.2/s_twa2mwh * OMF) %>% # cost per MWh VRE, #assuming that for one additional unit of VRE ratio of @@ -1459,28 +1459,46 @@ reportLCOE <- function(gdx, output.type = "both"){ # conversion from tr USD 2005/TW to USD2015/kW mutate(CAPEX = CAPEX *1.2 * 1e3) %>% # conversion from tr USD 2005/TWa to USD2015/MWh - mutate(OMV = OMV * 1.2 / s_twa2mwh * 1e12) %>% + mutate(OMV = OMV * 1.2 / as.numeric(s_twa2mwh) * 1e12) %>% # share of stored carbon from captured carbon is only relevant for CCS technologies, others -> 1 mutate( CO2StoreShare = ifelse(tech %in% teCCS, CO2StoreShare, 1)) %>% - # calculate discount factor for investment cost - mutate( disc.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) + # calculate annuity factor for annualizing investment cost over lifetime + mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) + #### calculate LCOE df.LCOE <- df.LCOE %>% + # investment cost LCOE in USD/MWh - mutate( `Investment Cost` = CAPEX * disc.fac / (CapFac*8760)*1e3) %>% + mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% # OMF cost LCOE in USD/MWh mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% mutate( `OMV Cost` = OMV) %>% mutate( `Fuel Cost` = fuel.price / eff) %>% + # CO2 Tax cost on SE level refer to (supply-side) emissions of pe2se technologies + # CO2 Tax cost on FE level refer to (demand-side) emissions of FE carriers mutate( `CO2 Tax Cost` = co2.price * (emiFac / eff) * CO2StoreShare) %>% mutate( `CO2 Provision Cost` = Co2.Capt.Price * co2_dem) %>% + # fuel cost of second fuel if technology has two inputs or outputs + # is positive for cost of second input + # is negative for benefit of a second output mutate( `Second Fuel Cost` = -(secfuel.prod * secfuel.price)) %>% - mutate( `VRE Storage Cost` = VREstor.cost, `Grid Cost` = grid.cost) %>% + # VRE integration cost (grid, storage, curtailment cost) + # note: marginal grid and storage cost are set to 0 for now, needs to be checked again + mutate( `VRE Storage Cost` = VREstor.cost, + `Grid Cost` = grid.cost) %>% + # curtailment cost are generation LCOE of VRE technologies of curtailed generation mutate( `Curtailment Cost` = curtShare / (1-curtShare) * (`Investment Cost` + `OMF Cost` + `OMV Cost`)) %>% + # CCS Tax cost as defined in 21_tax module mutate( `CCS Tax Cost` = CCStax.cost, `CCS Cost` = CCSCost) %>% + # Flex Tax benefit for electrolysis + # FlexPriceShare denotes share of the electricity price that electrolysis sees mutate( `Flex Tax` = -(1-FlexPriceShare) * `Fuel Cost`) %>% - mutate( `FE Tax` = FEtax, `Additional H2 t&d Cost` = AddH2TdCost) %>% + # se2fe technologies come with FE tax + mutate( `FE Tax` = FEtax, + # FE H2 has some additional t&d cost at low H2 shares (phase-in cost) in REMIND + `Additional H2 t&d Cost` = AddH2TdCost) %>% + # calculate total LCOE by adding all components mutate( `Total LCOE` = `Investment Cost` + `OMF Cost` + `OMV Cost` + `Fuel Cost` + `CO2 Tax Cost` + `CO2 Provision Cost` + `Second Fuel Cost` + `VRE Storage Cost` + `Grid Cost` + `CCS Tax Cost` + `Curtailment Cost` + `CCS Cost` + `Flex Tax` + `FE Tax` + `Additional H2 t&d Cost`) From 14545ca4da414bf86d0ca0d8f9ef8a04a083a8d2 Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Thu, 9 Nov 2023 11:32:09 +0100 Subject: [PATCH 3/8] clean up reportLCOE function --- R/reportLCOE.R | 382 ++++++++++--------------------------------------- 1 file changed, 74 insertions(+), 308 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index 062d1bce..e3c67371 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -80,12 +80,11 @@ reportLCOE <- function(gdx, output.type = "both"){ - ######################################################## - ### A) Calculation of average (standing system) LCOE ### - ######################################################## +#### A) Calculation of average (standing system) LCOE ---- + if (output.type %in% c("both", "average")) { - ####### read in needed data ####################### + # read in needed data ---- ## sets opTimeYr <- readGDX(gdx,"opTimeYr") @@ -146,7 +145,7 @@ reportLCOE <- function(gdx, output.type = "both"){ - ####### calculate needed parameters ############### + # calculates annuity cost: # annuity cost = 1/ sum_t (p_omeg(t) / 1.06^t) * direct investment cost @@ -179,15 +178,12 @@ reportLCOE <- function(gdx, output.type = "both"){ v_directteinv_wadj[,ttot_from2005,te] ) - ########## calculation of LCOE of standing system ####### - ######## (old annuities included) ###################### - + # average LCOE components ---- # calculate sub-parts of "COSTS" +# LCOE calculation only for pe2se technologies so far! - # LCOE calculation only for pe2se technologies so far! - - ####### 1. sub-part: Investment Cost ################################# +# 1. sub-part: Investment Cost ---- # AnnualInvCost(t) = sum_(td) [annuity(td) * p_pmeg(t-td+1) * deltaT], # where td in [t-lifetime,t] @@ -235,7 +231,7 @@ reportLCOE <- function(gdx, output.type = "both"){ } # a } # t0 - ####### 2. sub-part: fuel cost ################################# + # 2. sub-part: fuel cost ---- # fuel cost = PE price * PE demand of technology @@ -243,7 +239,7 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_fuel_cost[,,pe2se$all_te] <- setNames(1e+12 * qm_pebal[,ttot_from2005,pe2se$all_enty] / qm_budget[,ttot_from2005,] * setNames(vm_demPe[,,pe2se$all_te], pe2se$all_enty), pe2se$all_te) - ####### 3. sub-part: OMV cost ################################# + # 3. sub-part: OMV cost ---- # omv cost (from pm_data) * SE production @@ -252,7 +248,7 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_OMV_cost <- new.magpie(getRegions(te_inv_annuity),ttot_from2005,magclass::getNames(te_inv_annuity), fill=0) te_annual_OMV_cost[,,temapse$all_te] <- 1e+12 * collapseNames(pm_data[,,"omv"])[,,temapse$all_te] * setNames(vm_prodSe[,,temapse.names],temapse$all_te) - ####### 4. sub-part: OMF cost ################################# + # 4. sub-part: OMF cost ---- # omf cost = omf (from pm_data) * investcost (trUSD/TW) * capacity # omf in pm_data given as share of investment cost @@ -263,7 +259,7 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_OMF_cost[,,getNames(te_inv_annuity)] <- 1e+12 * collapseNames(pm_data[,,"omf"])[,,getNames(te_inv_annuity)] * v_investcost[,ttot_from2005,getNames(te_inv_annuity)] * dimSums(vm_cap[,ttot_from2005,], dim = 3.2)[,,getNames(te_inv_annuity)] - ####### 5. sub-part: storage cost (for wind, spv, csp) ################################# + # 5. sub-part: storage cost (for wind, spv, csp) ---- # storage cost = investment cost + omf cost # of corresponding storage technology ("storwind", "storspv", "storcsp") @@ -276,9 +272,9 @@ reportLCOE <- function(gdx, output.type = "both"){ te_annual_stor_cost_wadj <- new.magpie(getRegions(te_inv_annuity),ttot_from2005,magclass::getNames(te_inv_annuity), fill=0) te_annual_stor_cost_wadj[,,te2stor$all_te] <- setNames(te_annual_inv_cost_wadj[,ttot_from2005,te2stor$teStor] + te_annual_OMF_cost[,,te2stor$teStor],te2stor$all_te) - - - ####### 6. sub-part: grid cost ################################# + + + # 6. sub-part: grid cost ---- # same as for storage cost only with grid technologies: "gridwind", "gridspv", "gridcsp" # only "gridwind" technology active, wind requires 1.5 * the gridwind capacities as spv and csp @@ -306,9 +302,9 @@ reportLCOE <- function(gdx, output.type = "both"){ # see q32_limitCapTeGrid grid_factor_tech * vm_prodSe[,,te2grid$all_te] / vm_VRE_prodSe_grid - - - ####### 7. sub-part: ccs injection cost ################################# + + + # 7. sub-part: ccs injection cost ---- # same as for storage/grid but with ccs inejection technology # distributed to technolgies according to share of total captured co2 of ccs technology @@ -325,7 +321,7 @@ reportLCOE <- function(gdx, output.type = "both"){ dimSums( v_emiTeDetail[,,"cco2"], dim = 3, na.rm = T), teCCS) - ####### 8. sub-part: co2 cost ################################# + # 8. sub-part: co2 cost ---- te_annual_co2_cost <- new.magpie(getRegions(te_inv_annuity), ttot_from2005, getNames(te_inv_annuity), fill = 0) @@ -342,24 +338,24 @@ reportLCOE <- function(gdx, output.type = "both"){ ), tmp) rm(tmp) - + # regional part of the CO2 price (p47_taxCO2eq_AggFE only exits when regipol is run) if(!is.null(p47_taxCO2eq_AggFE)){ tmp <- intersect(getNames(te_inv_annuity), getNames(v_emiTeDetail[,,"co2"], dim = 3)) - + te_annual_co2_cost[,getYears(p47_taxCO2eq_AggFE),tmp] <- setNames( ( dimSums(p47_taxCO2eq_AggFE,dim=3, na.rm = TRUE) * dimSums(v_emiTeDetail[,getYears(p47_taxCO2eq_AggFE),tmp], dim = c(3.1, 3.2, 3.4), na.rm = TRUE) * 1e9 * 1e3 # T$/tC * GtC/yr * 1e9 t/Gt 1e3$/T$ = $/yr ), tmp) - + rm(tmp) } - - ####### 9. sub-part: curtailment cost ################################# - ########### (only relevant for RLDC version!) ########################## + + # 9. sub-part: curtailment cost ---- + # (only relevant for RLDC version!) # note: this step can only be done after all the other parts have already been calcuated! @@ -408,10 +404,9 @@ reportLCOE <- function(gdx, output.type = "both"){ - #################################################### - ######### calculate average LCOE ################## - #################################################### - LCOE.avg <- NULL + +# LCOE Calculation (average) ---- +LCOE.avg <- NULL # calculate standing system LCOE # divide total cost of standing system in that time step by total generation (before curtailment) in that time step @@ -482,11 +477,10 @@ reportLCOE <- function(gdx, output.type = "both"){ LCOE.out <- mbind(LCOE.out, LCOE.avg.out) } - ############################################## - ######################################################## - ### B) Calculation of marginal (new plant) LCOE ######## - ######################################################## + + #### B) Calculation of marginal (new plant) LCOE ---- + if (output.type %in% c("marginal", "both", "marginal detail")) { @@ -573,14 +567,12 @@ reportLCOE <- function(gdx, output.type = "both"){ fe2es.buildings <- pm_tau_fe_tax_ES_st <- pm_tau_fe_sub_ES_st <- NULL model <- scenario <- variable <- unit <- NULL - ############################################## + # Prepare data for LCOE calculation ---- - ##### Prepare data for LCOE calculation ####### - - ### 1. read variables, parameters sets, mappings needed from gdx + ### Read sets and mappings from GDX ---- ### technologies @@ -634,28 +626,25 @@ reportLCOE <- function(gdx, output.type = "both"){ filter(all_te %in% te_LCOE) colnames(en2en) <- c("fuel", "output", "tech") - # VRE to storage VRE2teStor <- readGDX(gdx, "VRE2teStor") %>% rename(tech = all_te) - ### time steps + ### time steps ttot <- as.numeric(readGDX(gdx,"ttot")) ttot_from2005 <- paste0("y",ttot[which(ttot >= 2005)]) ### conversion factors s_twa2mwh <- as.vector(readGDX(gdx,c("sm_TWa_2_MWh","s_TWa_2_MWh","s_twa2mwh"),format="first_found")) - # annuity factor from REMIND - p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) - - ### 2. retrieve investment and O&M cost + ### Read investment and O&M cost ---- # investment cost vm_costTeCapital <- readGDX(gdx, "vm_costTeCapital", field = "l", restore_zeros = F)[,ttot_from2005,te_LCOE_Inv] + df.CAPEX <- as.quitte(vm_costTeCapital) %>% select(region, period, all_te, value) %>% rename(tech = all_te, CAPEX = value) %>% @@ -676,7 +665,8 @@ reportLCOE <- function(gdx, output.type = "both"){ select(region, all_te, value) %>% rename(tech = all_te, OMV = value) - # 3. retrieve/calculate capacity factors + + # Read/calculate capacity factors ---- # capacity factor of non-renewables vm_capFac <- readGDX(gdx, "vm_capFac", field="l", restore_zeros = F)[,ttot_from2005,] @@ -718,7 +708,7 @@ reportLCOE <- function(gdx, output.type = "both"){ rbind(df.CapFac.ren) %>% mutate( period = as.numeric(period)) - ### 4. plant lifetime and annuity factor + ### Read plant lifetime and annuity factor ---- #discount rate r <- 0.051 @@ -741,14 +731,19 @@ reportLCOE <- function(gdx, output.type = "both"){ rename(tech = all_te, lifetime = value) %>% mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) + ### note: just take LCOE investment cost formula with fix lifetime, ignoring that in REMIND capacity is drepeciating over time ### actually you would need to divide CAPEX by CapFac * 8760 * p_omeg and then add a discounting. ### Would need to discuss again how to add the discount in such formula. - ### 6. retrieve fuel price + # # annuity factor from REMIND, + # TODO: check whether this is the same as calculated above + # p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) + + + ### Read fuel price ---- + - # retrieve capacity distribution over lifetime (fraction of capacity still standing in that year of plant lifetime) - p_omeg <- readGDX(gdx,c("pm_omeg","p_omeg"),format="first_found") # fuels to calculate price for fuels <- c("peoil","pegas","pecoal","peur", "pebiolc" , "pebios","pebioil", @@ -763,7 +758,7 @@ reportLCOE <- function(gdx, output.type = "both"){ Fuel.Price <- mbind(pm_PEPrice,pm_SEPrice )[,,fuels]*1e12/s_twa2mwh*1.2 - # Fuel price + # Fuel price for the time period for which LCOE are calculated df.Fuel.Price <- as.quitte(Fuel.Price) %>% select(region, period, all_enty, value) %>% rename(fuel = all_enty, fuel.price = value) %>% @@ -773,11 +768,10 @@ reportLCOE <- function(gdx, output.type = "both"){ tidyr::fill(fuel.price, .direction = "downup") %>% ungroup() - - - ### 7. retrieve carbon price + ### Read carbon price ---- p_priceCO2 <- readGDX(gdx,name=c("p_priceCO2","pm_priceCO2"),format="first_found") # co2 price + # carbon price for the time period for which LCOE are calculated df.co2price <- as.quitte(p_priceCO2) %>% select(region, period, value) %>% # where carbon price is NA, it is zero @@ -930,7 +924,7 @@ reportLCOE <- function(gdx, output.type = "both"){ - ### 12. calculate CO2 capture cost + ### Calculate CO2 capture cost ---- # temporary solution, take LCOD of DAC as the CO2 Capture Cost @@ -1027,7 +1021,8 @@ reportLCOE <- function(gdx, output.type = "both"){ # Note: This is assuming that the CO2 capture to storage share stays constant over time. # Still to do: calculate average over lifetime of plant - ### 14. calculate cost of second fuel (if required) +### Calculate cost of second fuel ---- + # for technologies with coupled production # (pm_prodCouple: negative values mean own consumption, positive values mean coupled product) @@ -1044,226 +1039,10 @@ reportLCOE <- function(gdx, output.type = "both"){ right_join(df.Fuel.Price, by = c("region", "fuel")) %>% rename(secfuel.price = fuel.price) - # ### VRE integration cost - # - # # first calculate LCOE of VREs before storage losses (Investment Cost + OMF cost) - # # (will be done below again, - # # this is only to calculate the integration cost which are added in the end) - # # also caculate LCOE (investment and OMF cost) of grid (only gridwind used) - # # and storage technologies - # df.LCOE.VRE.preloss <- df.CAPEX %>% - # left_join(df.OMF) %>% - # left_join(df.CapFac) %>% - # left_join(df.lifetime) %>% - # filter( tech %in% c(VRE2teStor[,1], VRE2teStor[,2], "gridwind")) %>% - # # conversion from tr USD 2005/TW to USD2015/kW - # mutate(CAPEX = CAPEX *1.2 * 1e3) %>% - # # calculate annuity factor for investment cost - # mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) %>% - # # investment cost LCOE in USD/MWh - # mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% - # # OMF cost LCOE in USD/MWh - # mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% - # # calculate pre-loss VRE LCOE as investment cost + OMF cost LCOE - # mutate( `VRE LCOE preloss` = `Investment Cost` + `OMF Cost`) %>% - # select(region, period, tech, CapFac, `VRE LCOE preloss`) - # - # - # - # # calculate marginal storloss: - # # (Marginal amount of energy that is lost in storage when - # # one more MWh of 'usable seel' (='upgraded seel that is balanced by storage') - # # is produced. Unit: MWh. Can be larger than 1!) - # - # # generation share of technology - # df.shSeEl <- read.gdx(gdx, c("vm_shSeEl","v32_shSeEl"), fields = "l", format= "first_found") %>% - # rename(shSeEl = value, tech = all_te, period = ttot, region = all_regi) - # # usable (after storage loss) SE per technology - # df.usableSeTe <- read.gdx(gdx, "vm_usableSeTe", fields = "l") %>% - # rename(usableSeTe = value, tech = all_te, period = ttot, region = all_regi) %>% - # select(-entySe) - # # usable (after storage loss) SE - # df.usableSe <- read.gdx(gdx, "vm_usableSe", fields = "l") %>% - # rename(usableSe = value, period = ttot, region = all_regi) %>% - # select(-entySe) - # # storage exponent - # df.storexp <- read.gdx(gdx, "p32_storexp") %>% - # rename(storexp = value, tech = all_te, region = all_regi) - # - # - # # questions: !! - # ## marginal storloss of usable SE partly very high (30 for wind EUR), which unit does this have - # # it is a share, right? Or does it have energy units? - # - # # calculate marginal storloss - # df.marg.storloss <- df.LCOE.VRE.preloss %>% - # filter( tech %in% VRE2teStor[,1]) %>% - # left_join(df.shSeEl) %>% - # left_join(df.usableSeTe) %>% - # left_join(df.usableSe) %>% - # left_join(df.eff %>% - # right_join(VRE2teStor) %>% - # select(-teStor) ) %>% - # left_join(df.storexp) %>% - # mutate( loss = 1-eff) %>% - # # formula derived by RP: - # # This formula is an approximation of the marginal change of storage loss - # # that was calculated outside ReMIND by hand. - # mutate( first.term = (loss * (shSeEl^storexp) + - # (usableSeTe * shSeEl^(storexp-1) * - # (1/usableSe - usableSeTe/(usableSe^2)))) / - # ((1-eff)*shSeEl^storexp)) %>% - # mutate( second.term = ( loss^2 * shSeEl^storexp * storexp * shSeEl^(storexp-1) * - # (1/usableSe - usableSeTe/(usableSe^2))) / - # (((1-eff)*shSeEl^storexp)^2)) %>% - # mutate( marg.storloss = first.term + second.term) %>% - # mutate( marg.storloss = ifelse(loss == 0, 0, marg.storloss)) %>% - # # convert storage loss from usable SE to total SE level - # mutate(marg.storloss.SE = marg.storloss / (1+marg.storloss)) %>% - # # calculate marginal LCOE of storloss - # mutate( loss.MLCOE = `VRE LCOE preloss` * marg.storloss.SE / - # ( 1 - marg.storloss.SE)) %>% - # # calculate marginal LCOE postloss - # mutate( postloss.LCOE = `VRE LCOE preloss` / - # ( 1 - marg.storloss.SE)) %>% - # # calculate postloss - preloss LCOE -> VRE Storage Loss Cost - # mutate( storloss.MLCOE = postloss.LCOE-`VRE LCOE preloss`) - # - # # VRE grid weight: how much grid capacities the respective VRE technology needs - # df.VRE2teGrid <- data.frame(tech = c("spv","csp","wind"), GridReq = c(1,1,1.5)) - # - # - # # calculate marginal integration cost (cost of building storage and grid for the marginal VRE unit) - # df.IntCost <- df.marg.storloss %>% - # select(region, period, tech, marg.storloss) %>% - # # join with MCLOE for storage and grid technologies - # # add grid capacity MLCOE - # left_join(df.LCOE.VRE.preloss %>% - # filter(tech == "gridwind") %>% - # select(-tech) %>% - # rename(gridCap.MLCOE = `VRE LCOE preloss`)) %>% - # # add storage capacity MLCOE - # left_join(df.LCOE.VRE.preloss %>% - # filter(tech %in% VRE2teStor[,2]) %>% - # rename( teStor = tech, StorCap.MLCOE = `VRE LCOE preloss`) %>% - # left_join(VRE2teStor) ) %>% - # # add grid requirement factor, current IntC impl.: wind needs 1.5*gridwind capacity than solar - # left_join(df.VRE2teGrid) %>% - # # add storage efficiencies - # left_join(df.eff %>% - # filter(tech %in% VRE2teStor[,2]) %>% - # rename(teStor = tech)) %>% - # # calculate marginal grid and storage capacity cost per unit of new VRE - # mutate( MLCOE.grid = marg.storloss / CapFac * GridReq * gridCap.MLCOE, - # MLCOE.stor = marg.storloss * eff / (1-eff) / CapFac * StorCap.MLCOE) - # - # - # - # - # - # - # # o_te_marg_storloss_usable(ttot,regi,teReNoBio) !!RP: This formula is an approximation of the marginal change of storage loss that was calculated outside ReMIND by hand. - # # = p_aux_loss(teReNoBio) - # # * (p_aux_shareseel ** p_storexp(regi,teReNoBio) - # # + ( - # # p_aux_usablese_te * p_aux_shareseel ** (p_storexp(regi,teReNoBio) -1) - # # * ( 1 / p_aux_usablese - p_aux_usablese_te/(p_aux_usablese ** 2) ) - # # ) - # # ) - # # / (1 - p_aux_loss(teReNoBio) * p_aux_shareseel ** p_storexp(regi,teReNoBio) ) - # # + ( - # # p_aux_loss(teReNoBio) ** 2 * p_aux_shareseel ** p_storexp(regi,teReNoBio) * p_storexp(regi,teReNoBio) * p_aux_shareseel ** (p_storexp(regi,teReNoBio) -1) - # # * ( 1 / p_aux_usablese - p_aux_usablese_te/(p_aux_usablese ** 2) ) - # # ) - # # / ( - # # (1 - p_aux_loss(teReNoBio) * p_aux_shareseel ** p_storexp(regi,teReNoBio) ) - # # ** 2 - # # ); - # - # - # - # #### calculate LCOE - # df.LCOE <- df.LCOE %>% - # # investment cost LCOE in USD/MWh - # mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% - # # OMF cost LCOE in USD/MWh - # mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% - # mutate( `OMV Cost` = OMV) %>% - - ### 15. calculate grid cost of VRE technologies - teVRE.grid <- data.frame(tech=c("spv","csp","wind"), gridtech = c("gridwind", - "gridwind","gridwind") ) - - p32_grid_factor <- readGDX(gdx, "p32_grid_factor") - - df.gridfactor <- as.quitte(p32_grid_factor) %>% - rename(grid.factor = value) %>% - select(region, grid.factor) - - df.gridcost <- df.CAPEX %>% - left_join(df.OMF, by = c("region", "tech")) %>% - left_join(df.lifetime, by = "tech") %>% - left_join(df.gridfactor, by = "region") %>% - filter(tech == "gridwind") %>% - rename(gridtech = tech) %>% - mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * annuity.fac) %>% - mutate(`OMF Cost` = CAPEX*1e12*1.2/s_twa2mwh * OMF) %>% - mutate(grid.cost = (`Investment Cost` + `OMF Cost`) / grid.factor) %>% - full_join(teVRE.grid, by = "gridtech") %>% - # wind requires 1.5 times the grid capacity than spv and csp - mutate(grid.cost = ifelse(tech == "wind", 1.5*grid.cost, grid.cost)) %>% - select(region, period, tech, grid.cost) - - - # Note: marginal grid cost calculation needs to be checked, see Robert's old version above, - # set to 0 meanwhile - df.gridcost <- df.gridcost %>% - mutate( grid.cost = 0) - - # 16. VRE electricity storage cost - - v32_storloss <- readGDX(gdx, "v32_storloss", field = "l") - # VRE production - vm_prodSe <- readGDX(gdx, "vm_prodSe", field = "l", restore_zeros = F)[,,c("spv","wind","csp")] - - df.prodSeVRE <- as.quitte(vm_prodSe) %>% - rename(tech = all_te, vm_prodSeVRE = value) %>% - select(region, period, tech, vm_prodSeVRE) - - df.storloss <- as.quitte(v32_storloss) %>% - rename(storloss = value, tech = all_te) %>% - select(region, period, tech, storloss) - - - # calculation following q32_limitCapTeStor - df.VREstorcost <- df.CAPEX %>% - filter(tech %in% VRE2teStor$teStor) %>% - left_join(df.OMF, by = c("region", "tech")) %>% - left_join(df.lifetime, by = "tech") %>% - rename(teStor = tech) %>% - full_join(VRE2teStor, by = "teStor") %>% - left_join(df.storloss, by = c("region", "period", "tech")) %>% - left_join(df.prodSeVRE, by = c("region", "period", "tech")) %>% - left_join(df.eff, by=c("region"="region","period"="period","teStor"="tech")) %>% - left_join(df.CapFac, by=c("region"="region","period"="period","teStor"="tech")) %>% - # cost per MWh v32_storloss - mutate(`Investment Cost` = CAPEX *1e12*1.2/s_twa2mwh * annuity.fac) %>% - mutate(`OMF Cost` = CAPEX*1e12*1.2/s_twa2mwh * OMF) %>% - # cost per MWh VRE, - #assuming that for one additional unit of VRE ratio of - #v32_storloss/vm_prodSe constant - mutate(VREstor.cost = (`Investment Cost` + `OMF Cost`) / - CapFac / eff * (1-eff) * storloss / vm_prodSeVRE) %>% - mutate( VREstor.cost = ifelse(is.na(VREstor.cost), 0, VREstor.cost)) %>% - select(region, period, tech, VREstor.cost) - - # Note: marginal storage cost calculation needs to be checked, see Robert's old version above, - # set to 0 meanwhile - df.VREstorcost <- df.VREstorcost %>% - mutate( VREstor.cost = 0) - - - ### calculate Curtailment share (required for calculating Curtailment Cost later for VREs) + + + +### Read curtailment share ----- # curtailment share is needed to calculate curtailment cost of VREs: # curtailment cost = LCOE(VRE) * curtshare/((1-curtShare)), @@ -1285,7 +1064,8 @@ reportLCOE <- function(gdx, output.type = "both"){ - ### 17. CCS tax +### Calculate CCS tax ---- + # following q21_taxrevCCS vm_co2CCS <- readGDX(gdx, "vm_co2CCS", field = "l", restore_zeros = F) sm_ccsinjecrate <- readGDX(gdx, c("sm_ccsinjecrate", "s_ccsinjecrate"), format = "first_found") @@ -1319,22 +1099,11 @@ reportLCOE <- function(gdx, output.type = "both"){ - ### 18. CO2 Storage Cost - p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) - # adding annualized capital cost and omf cost of ccsinje technology, - # multiply with vm_co2CCS_m (stored CO2 of CCS technology in GtC/TWa(output)), - # convert to USD2015/MWh - CCSCost <- vm_co2CCS_m*dimReduce(vm_costTeCapital[,,"ccsinje"]/vm_capFac[,,"ccsinje"]*(p_teAnnuity[,,"ccsinje"]+pm_data_omf[,,"ccsinje"]))*1.2*1e12/s_twa2mwh - df.CCSCost <- as.quitte(CCSCost) %>% - rename(tech = all_te, CCSCost = value) %>% - select(region, period, tech, CCSCost) +# mutate( CCSCost = 0) - # Note: CCS cost still to fix, set temporarily to 0 - df.CCSCost <- df.CCSCost %>% - mutate( CCSCost = 0) +### Read Flexibility Tax ---- - ### 19. Flexibility Tax cm_FlexTax <- readGDX(gdx, "cm_flex_tax") v32_flexPriceShare <- readGDX(gdx, "v32_flexPriceShare", field = "l", restore_zeros = F) if (is.null(v32_flexPriceShare) | is.null(cm_FlexTax)) { @@ -1351,7 +1120,7 @@ reportLCOE <- function(gdx, output.type = "both"){ rename(tech = all_te, FlexPriceShare = value) %>% select(region, period, tech, FlexPriceShare) - ### 20. Final Energy Taxes +### Read Final Energy Taxes ---- entyFe2Sector <- readGDX(gdx, "entyFe2Sector") @@ -1369,7 +1138,8 @@ reportLCOE <- function(gdx, output.type = "both"){ - ### add additional H2 transmission and distribution cost +### read H2 Phase-in cost for buildings and industry ---- + ### for buildings and industry (calculate average cost for simplicity, not marginal) @@ -1402,8 +1172,7 @@ reportLCOE <- function(gdx, output.type = "both"){ - ####################### LCOE calculation (New plant/marginal) ######################## - + # Join all data for marginal LCOE calculation ---- ### create table with all parameters needed for LCOE calculation df.LCOE <- df.CAPEX %>% @@ -1420,17 +1189,15 @@ reportLCOE <- function(gdx, output.type = "both"){ left_join(df.co2_dem, by = c("region", "period", if (module2realisation["CCU",2] == "on") "tech")) %>% left_join(df.CO2StoreShare, by = c("region", "period")) %>% left_join(df.secfuel, by = c("region", "period", "tech", "fuel")) %>% - left_join(df.gridcost, by = c("region", "period", "tech")) %>% - left_join(df.VREstorcost, by = c("region", "period", "tech")) %>% left_join(df.curtShare, by = c("region", "period", "tech")) %>% left_join(df.CCStax, by = c("region", "period", "tech")) %>% - left_join(df.CCSCost, by = c("region", "period", "tech")) %>% left_join(df.flexPriceShare, by = c("region", "period", "tech")) %>% left_join(df.FEtax,relationship = "many-to-many", by = c("region", "period", "output")) %>% left_join(df.AddTeInvH2, by = c("region", "period", "tech", "sector")) %>% # filter to only have LCOE technologies filter( tech %in% c(te_LCOE)) + # only retain unique region, period, tech combinations df.LCOE <- df.LCOE %>% unique(by=c("region", "period", "tech")) @@ -1439,7 +1206,7 @@ reportLCOE <- function(gdx, output.type = "both"){ # replace NA by 0 in certain columns # columns where NA should be replaced by 0 col.NA.zero <- c("OMF","OMV", "co2.price", "fuel.price","co2_dem","emiFac.se2fe","Co2.Capt.Price", - "secfuel.prod", "secfuel.price", "grid.cost","VREstor.cost", "curtShare","CCStax.cost","CCSCost","FEtax","AddH2TdCost") + "secfuel.prod", "secfuel.price", "curtShare","CCStax.cost","FEtax","AddH2TdCost") df.LCOE[,col.NA.zero][is.na(df.LCOE[,col.NA.zero])] <- 0 # replace NA by 1 in certain columns @@ -1466,9 +1233,10 @@ reportLCOE <- function(gdx, output.type = "both"){ mutate( annuity.fac = r * (1+r)^lifetime/(-1+(1+r)^lifetime)) - #### calculate LCOE - df.LCOE <- df.LCOE %>% +# LCOE Calculation (marginal) ---- + + df.LCOE <- df.LCOE %>% # investment cost LCOE in USD/MWh mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% # OMF cost LCOE in USD/MWh @@ -1483,14 +1251,10 @@ reportLCOE <- function(gdx, output.type = "both"){ # is positive for cost of second input # is negative for benefit of a second output mutate( `Second Fuel Cost` = -(secfuel.prod * secfuel.price)) %>% - # VRE integration cost (grid, storage, curtailment cost) - # note: marginal grid and storage cost are set to 0 for now, needs to be checked again - mutate( `VRE Storage Cost` = VREstor.cost, - `Grid Cost` = grid.cost) %>% # curtailment cost are generation LCOE of VRE technologies of curtailed generation mutate( `Curtailment Cost` = curtShare / (1-curtShare) * (`Investment Cost` + `OMF Cost` + `OMV Cost`)) %>% # CCS Tax cost as defined in 21_tax module - mutate( `CCS Tax Cost` = CCStax.cost, `CCS Cost` = CCSCost) %>% + mutate( `CCS Tax Cost` = CCStax.cost) %>% # Flex Tax benefit for electrolysis # FlexPriceShare denotes share of the electricity price that electrolysis sees mutate( `Flex Tax` = -(1-FlexPriceShare) * `Fuel Cost`) %>% @@ -1500,10 +1264,12 @@ reportLCOE <- function(gdx, output.type = "both"){ `Additional H2 t&d Cost` = AddH2TdCost) %>% # calculate total LCOE by adding all components mutate( `Total LCOE` = `Investment Cost` + `OMF Cost` + `OMV Cost` + `Fuel Cost` + `CO2 Tax Cost` + - `CO2 Provision Cost` + `Second Fuel Cost` + `VRE Storage Cost` + `Grid Cost` + `CCS Tax Cost` + `Curtailment Cost` + `CCS Cost` + `Flex Tax` + `FE Tax` + `Additional H2 t&d Cost`) + `CO2 Provision Cost` + `Second Fuel Cost` + `CCS Tax Cost` + `Curtailment Cost` + `Flex Tax` + `FE Tax` + `Additional H2 t&d Cost`) + # Levelized Cost of UE in Buildings Putty Realization ---- + # if detailed buildings module on, calculate LCOE of fe2ue technologies # following 36_modules/buildings/services_putty/presolve.gms @@ -1572,7 +1338,7 @@ reportLCOE <- function(gdx, output.type = "both"){ } - +# Transform LCOE to output format ---- # transform marginal LCOE to mif output format df.LCOE.out <- df.LCOE %>% From 547c6fdd0391ce6a3578d722aa945069f2aa7cae Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Thu, 9 Nov 2023 11:33:49 +0100 Subject: [PATCH 4/8] reintroduce calculation of fuel cost and co2 tax cost in reportLCOE based on discounted and capacity-weighted averages over the operation time of the plant --- R/reportLCOE.R | 218 +++++++++++++++++++++++++------------------------ 1 file changed, 111 insertions(+), 107 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index e3c67371..5eebfbeb 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -24,7 +24,7 @@ #' @importFrom gdx readGDX #' @importFrom magclass new.magpie dimSums getRegions getYears getNames setNames clean_magpie dimReduce as.magpie magpie_expand #' @importFrom dplyr %>% mutate select rename group_by ungroup right_join filter full_join arrange summarise -#' @importFrom quitte as.quitte overwrite getRegs +#' @importFrom quitte as.quitte overwrite getRegs getPeriods #' @importFrom tidyr spread gather expand fill @@ -779,110 +779,104 @@ LCOE.avg <- NULL mutate(value = ifelse(is.na(value), 0, value*1.2/3.66)) %>% rename(co2.price = value) - - - ### TODO: maybe deleted if not necessary anymore:this is the old section where - # capacity-weighted/discounted averages of fuel and co2 prices are calculated - # over the lifetime of the plant - - # - # # Fuel price - # df.Fuel.Price <- as.quitte(Fuel.Price) %>% - # select(region, period, all_enty, value) %>% - # rename(fuel = all_enty, fuel.price = value, opTimeYr = period) %>% - # # replace zeros by last or (if not available) next value of time series (marginals are sometimes zero if there are other constriants) - # mutate( fuel.price = ifelse(fuel.price == 0, NA, fuel.price)) %>% - # group_by(region, fuel) %>% - # tidyr::fill(fuel.price, .direction = "downup") %>% - # ungroup() - # - # - - # - # ### 8. calculate capacity distribution weighted prices over lifetimes of plant - # - # # capacity distribution over lifetime - # df.pomeg <- as.quitte(p_omeg) %>% - # rename(tech = all_te, pomeg = value) %>% - # # to save run time, only take p_omeg in ten year time steps only up to lifetimes of 50 years, - # # erase region dimension, pomeg same across all regions, only se generating technologies - # filter(opTimeYr %in% c(1,seq(5,50,5)), region %in% getRegions(vm_costTeCapital)[1] - # , tech %in% te_LCOE) %>% - # select(opTimeYr, tech, pomeg) - # - # # expanded df.omeg, with additional period dimension combined with opTimeYr dimensions, - # # period will be year of building the plant, - # # opTimeYr will be year within lifetime of plant (where only a certain fraction, pomeg, of capacity is still standing) - # df.pomeg.expand <- df.pomeg %>% - # expand(opTimeYr, period = seq(2005,2150,5)) %>% - # full_join(df.pomeg) %>% - # mutate( opTimeYr = as.numeric(opTimeYr)) %>% - # mutate( opTimeYr = ifelse(opTimeYr > 1, opTimeYr + period, period)) %>% - # left_join(en2en) %>% - # arrange(tech, period, opTimeYr) - # - - - ### 9. weight fuel price and carbon price with capacity density over liftime - - # ### TODO: Is that correct ?fuel price and co2 price with capacity-weighting and discounting - # - # # join fuel price to df.omeg and weigh fuel price by df.pomeg (pomeg = fraction of capacity still standing in that year of plant lifetime) - # df.fuel.price.weighted <- df.Fuel.Price %>% - # full_join(df.pomeg.expand) %>% - # # calculate discount factor over years at 5%/yr - # mutate( discount = 0.95^(as.numeric(opTimeYr)-as.numeric(period))) %>% - # # mean of fuel prices over first 50-years of plant lifetime weighted by pomeg and discount factor - # group_by(region, period, fuel, tech) %>% - # summarise( fuel.price.weighted.mean = mean(fuel.price * pomeg * discount)) %>% - # ungroup() - - - # ### 7. retrieve carbon price - # df.co2price <- as.quitte(p_priceCO2) %>% - # select(region, period, value) %>% - # # where carbon price is NA, it is zero - # # convert from USD2005/tC CO2 to USD2015/tCO2 - # mutate(value = ifelse(is.na(value), 0, value*1.2/3.66)) %>% - # rename(co2.price = value, opTimeYr = period) - - # - # # join carbon price to df.omeg and weigh it by df.pomeg (pomeg = fraction of capacity still standing in that year of plant lifetime) - # df.co2price.weighted <- df.co2price %>% - # full_join(df.pomeg.expand) %>% - # filter( !is.na(tech)) %>% - # # calculate discount factor over years at 5%/yr - # mutate( discount = 0.95^(as.numeric(opTimeYr)-as.numeric(period))) %>% - # # mean of carbon price over first 50-years of plant lifetime weighted by pomeg and discount factor - # group_by(region, period, tech) %>% - # summarise( co2.price.weighted.mean = mean(co2.price * pomeg * discount)) %>% - # ungroup() - # - # - # - - # # join fuel price to df.omeg and weigh fuel price by df.pomeg (pomeg = fraction of capacity still standing in that year of plant lifetime) - # df.fuel.price.weighted <- df.Fuel.Price %>% - # full_join(df.pomeg.expand) %>% - # # mean of fuel prices over first 50-years of plant lifetime weighted by pomeg - # group_by(region, period, fuel, tech) %>% - # summarise( fuel.price.weighted.mean = sum(fuel.price * pomeg / sum(pomeg))) %>% - # ungroup() - # - # # join carbon price to df.omeg and weigh it by df.pomeg (pomeg = fraction of capacity still standing in that year of plant lifetime) - # df.co2price.weighted <- df.co2price %>% - # full_join(df.pomeg.expand) %>% - # filter( !is.na(tech)) %>% - # # mean of carbon price over first 50-years of plant lifetime weighted by pomeg - # group_by(region, period, tech) %>% - # summarise( co2.price.weighted.mean = sum(co2.price * pomeg / sum(pomeg))) %>% - # ungroup() - - - - - - ### 10. get fuel conversion efficiencies + # Calculate weighted averages of fuel and carbon prices over time ---- + + # This section calculates weighted-average of the prices seen by a plant over + # its time of operation. It takes into account that fuel cost are discounted over time + # and that capacities depreciate over time such that near-term prices + # should be weighted higher in the LCOE than long-term prices. + # Weighted average fuel prices FP_avg are calculated + # from fuel prices over period t FP_t as: + # FP_avg = sum_t ( (discount_t * pomeg_t * FP_t) / (sum_t (discount_t * pomeg_t)) + # discount_t: discounting factor + # pomeg_t: depreciation factor, gives fraction of operating capacity over lifetime of plant + + + # retrieve capacity distribution over lifetime + # (fraction of capacity still standing in that year of plant lifetime) + p_omeg <- readGDX(gdx,c("pm_omeg","p_omeg"),format="first_found") + + # operation time steps of plant, only take every 5 years and limit calculation to 50 years as + # by then most of the capacity of all technologies is already retired + set_operation_period <- c(1,seq(5,50,5)) + + + # capacity distribution over lifetime + df.pomeg <- as.quitte(p_omeg) %>% + rename(tech = all_te, pomeg = value) %>% + # to save run time, only take p_omeg in ten year time steps only up to lifetimes of 50 years, + # erase region dimension, pomeg same across all regions + filter(opTimeYr %in% set_operation_period, + region %in% getRegions(vm_costTeCapital)[1], + tech %in% te_LCOE) %>% + select(opTimeYr, tech, pomeg) + + # create dataframe with all combinations of time steps "period" in which plant was built (commissioning period) + # time steps "operationPeriod" in which production takes place, year of operation "opTimeYr" in lifetime of plant + df.period_operationPeriod <- expand.grid(opTimeYr = unique(df.pomeg$opTimeYr), + period = getPeriods(df.Fuel.Price)) %>% + # operationPeriod is commissioning period + years of operation (opTimeYr) + # for comissioning period take value of first year of operation + mutate( operationPeriod = ifelse(as.numeric(opTimeYr) > 1, + as.numeric(opTimeYr) + period, + period)) %>% + # remove time steps from operationPeriod that are not remind_timesteps + filter( operationPeriod %in% unique(remind_timesteps$period)) + + # Join capacity distribution with above dataframe to get + # capacity distribution over all commissioning years (period) and + # operation years (operationPeriod) + df.pomeg.expand <- df.pomeg %>% + # only take capacity distribution of 5-year time steps and up to 50 years of operation + right_join(df.period_operationPeriod, + relationship = "many-to-many") %>% + # add energy input and output carrier dimension + left_join(en2en) + +# calculate average capacity-weighted and discounted fuel price over lifetime of plant +df.fuel.price.weighted <- df.pomeg.expand %>% + left_join(df.Fuel.Price, + by = c("operationPeriod" = "period", + "fuel" = "fuel"), + relationship = "many-to-many") %>% + # calculate discount factor + mutate( discount =1 / (1+r)^(operationPeriod - period)) %>% + # calculate weights over lifetime as operational capacity * discount factor, + # normalize by sum over all time steps for weights to add up to one + group_by(region, period, tech) %>% + mutate( weight = pomeg * discount / sum(pomeg * discount)) %>% + # calculate average weighted fuel price over lifetime as + # the sum of the capacity-discount weights + summarise( fuel.price.weighted = sum(fuel.price * weight)) %>% + ungroup() + + +df.fuel.price.weighted.check <- df.fuel.price.weighted %>% + filter(period == 2040, + region == "MEA", + tech == "refliq") + +# calculate average capacity-weighted and discounted fuel price over lifetime of plant +df.co2price.weighted <- df.pomeg.expand %>% + left_join(df.co2price, + by = c("operationPeriod" = "period"), + relationship = "many-to-many") %>% + # calculate discount factor + mutate( discount =1 / (1+r)^(operationPeriod - period)) %>% + # calculate weights over lifetime as operational capacity * discount factor, + # normalize by sum over all time steps for weights to add up to one + group_by(region, period, tech) %>% + mutate( weight = pomeg * discount / sum(pomeg * discount)) %>% + # calculate average weighted fuel price over lifetime as + # the sum of the capacity-discount weights + summarise( co2.price.weighted = sum(co2.price * weight)) %>% + ungroup() + + + + + + ### Read conversion efficiencies ---- pm_eta_conv <- readGDX(gdx,"pm_eta_conv", restore_zeros=F)[,ttot_from2005,] # efficiency oftechnologies with time-independent eta pm_dataeta <- readGDX(gdx,"pm_dataeta", restore_zeros=F)[,ttot_from2005,]# efficiency of technologies with time-dependent eta @@ -1182,6 +1176,8 @@ LCOE.avg <- NULL left_join(df.lifetime, by = "tech") %>% left_join(df.Fuel.Price, by = c("region", "period", "fuel")) %>% left_join(df.co2price, by = c("region", "period")) %>% + left_join(df.fuel.price.weighted, by = c("region", "period", "tech")) %>% + left_join(df.co2price.weighted, by = c("region", "period", "tech")) %>% left_join(df.eff, by = c("region", "period", "tech")) %>% left_join(df.emiFac, by = c("region", "tech")) %>% left_join(df.emifac.se2fe, by = c("region", "tech")) %>% @@ -1242,10 +1238,18 @@ LCOE.avg <- NULL # OMF cost LCOE in USD/MWh mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% mutate( `OMV Cost` = OMV) %>% - mutate( `Fuel Cost` = fuel.price / eff) %>% + # Fuel Cost + # # Fuel cost with fuel price in commissioning time step of plant + # mutate( `Fuel Cost` = fuel.price / eff) %>% + # Fuel cost with weighted average fuel price over plant lifetime + mutate( `Fuel Cost` = fuel.price.weighted / eff) %>% + # CO2 Tax cost # CO2 Tax cost on SE level refer to (supply-side) emissions of pe2se technologies # CO2 Tax cost on FE level refer to (demand-side) emissions of FE carriers - mutate( `CO2 Tax Cost` = co2.price * (emiFac / eff) * CO2StoreShare) %>% + # # CO2 Tax cost with co2 price in commissioning time step of plant + # mutate( `CO2 Tax Cost` = co2.price * (emiFac / eff) * CO2StoreShare) %>% + # CO2 Tax cost with weighted average co2 price over plant lifetime + mutate( `CO2 Tax Cost` = co2.price.weighted * (emiFac / eff) * CO2StoreShare) %>% mutate( `CO2 Provision Cost` = Co2.Capt.Price * co2_dem) %>% # fuel cost of second fuel if technology has two inputs or outputs # is positive for cost of second input From 282bacdadd54765f3add5b8e8f40bd23345cd777 Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Wed, 22 Nov 2023 11:03:35 +0100 Subject: [PATCH 5/8] add marginal adjustment cost LCOE calculation --- R/reportLCOE.R | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index 5eebfbeb..337f6725 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -740,6 +740,17 @@ LCOE.avg <- NULL # TODO: check whether this is the same as calculated above # p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) + ### Read marginal adjustment costs ---- + + # Read marginal adjustment cost calculated in core/postsolve.gms + # It is calculated as d(v_costInvTeAdj) / d(vm_deltaCap). + # Unit: trUSD2005/ (TW(out)/yr). + o_margAdjCostInv <- readGDX(gdx, "o_margAdjCostInv", restore_zeros = F) + + df.margAdjCostInv <- as.quitte(o_margAdjCostInv) %>% + rename(tech = all_te, + AdjCost = value) %>% + select( region, period, tech, AdjCost) ### Read fuel price ---- @@ -1174,6 +1185,7 @@ df.co2price.weighted <- df.pomeg.expand %>% left_join(df.OMV, by = c("region", "tech")) %>% left_join(df.CapFac, by = c("region", "period", "tech")) %>% left_join(df.lifetime, by = "tech") %>% + left_join(df.margAdjCostInv, by = c("region", "period", "tech")) %>% left_join(df.Fuel.Price, by = c("region", "period", "fuel")) %>% left_join(df.co2price, by = c("region", "period")) %>% left_join(df.fuel.price.weighted, by = c("region", "period", "tech")) %>% From 8dd1ab682d2fdf61b01672f8dae33fd2b1ed61e3 Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Wed, 22 Nov 2023 11:05:10 +0100 Subject: [PATCH 6/8] clean up LCOE reporting, add LCOX calculation for ccsinje technology and add two version of LCOE calculate including either fuel and carbon prices of time step or intertemporal prices --- R/reportLCOE.R | 153 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 112 insertions(+), 41 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index 337f6725..9327a271 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -5,6 +5,8 @@ #' This script calculates two different types of LCOE: average LCOE (standing system) and marginal LCOE (new plant). #' The average LCOE reflect the total cost incurred by the technology deployment in a specific time step divided by its energy output. #' The marginal LCOE estimate the per-unit lifetime cost of the output if the model added another capacity of that technology in the respective time step. +#' The marginal LCOE are calculate in two versions: 1) with fuel prices and co2 taxes of the time step for which the LCOE are calculated (time step prices), +#' or 2) with an intertemporal weighted-average of fuel price and co2 taxes over the lifetime of the plant (intertemporal prices). #' #' #' @@ -78,6 +80,9 @@ reportLCOE <- function(gdx, output.type = "both"){ qm_pebal <- readGDX(gdx,name=c("q_balPe"),field="m",format="first_found") qm_budget <- readGDX(gdx,name=c("qm_budget"),field="m",format="first_found") + ## variables + vm_prodSe <- readGDX(gdx,name=c("vm_prodSe"),field="l",restore_zeros=FALSE,format="first_found") + #### A) Calculation of average (standing system) LCOE ---- @@ -126,7 +131,6 @@ reportLCOE <- function(gdx, output.type = "both"){ vm_capEarlyReti <- readGDX(gdx,name=c("vm_capEarlyReti"),field="l",format="first_found")[,ttot,] vm_deltaCap <- readGDX(gdx,name=c("vm_deltaCap"),field="l",format="first_found")[,ttot,] vm_demPe <- readGDX(gdx,name=c("vm_demPe","v_pedem"),field="l",restore_zeros=FALSE,format="first_found") - vm_prodSe <- readGDX(gdx,name=c("vm_prodSe"),field="l",restore_zeros=FALSE,format="first_found") v_investcost <- readGDX(gdx,name=c("vm_costTeCapital","v_costTeCapital","v_investcost"),field="l",format="first_found")[,ttot,] vm_cap <- readGDX(gdx,name=c("vm_cap"),field="l",format="first_found") vm_prodFe <- readGDX(gdx,name=c("vm_prodFe"),field="l",restore_zeros=FALSE,format="first_found") @@ -606,10 +610,10 @@ LCOE.avg <- NULL # all technologies to calculate LCOE for - te_LCOE <- c(pe2se$all_te, se2se$all_te,se2fe$all_te) + te_LCOE <- c(pe2se$all_te, se2se$all_te,se2fe$all_te, "ccsinje") - # all technologies to calculate investment and O&M LCOE for (needed for CCS, storage, grid cost) - te_LCOE_Inv <- c(te_LCOE, as.vector(teStor), as.vector(teGrid), ccs2te$all_te, te[te %in% c("dac")]) + # all technologies to calculate investment, adjustment and O&M LCOE for (needed for storage, grid cost) + te_LCOE_Inv <- c(te_LCOE, as.vector(teStor), as.vector(teGrid), te[te %in% c("dac")]) # technologies to produce SE te_SE_gen <- c(pe2se$all_te, se2se$all_te) @@ -790,7 +794,7 @@ LCOE.avg <- NULL mutate(value = ifelse(is.na(value), 0, value*1.2/3.66)) %>% rename(co2.price = value) - # Calculate weighted averages of fuel and carbon prices over time ---- + # Calculate intertemporal fuel and carbon prices ---- # This section calculates weighted-average of the prices seen by a plant over # its time of operation. It takes into account that fuel cost are discounted over time @@ -891,7 +895,7 @@ df.co2price.weighted <- df.pomeg.expand %>% pm_eta_conv <- readGDX(gdx,"pm_eta_conv", restore_zeros=F)[,ttot_from2005,] # efficiency oftechnologies with time-independent eta pm_dataeta <- readGDX(gdx,"pm_dataeta", restore_zeros=F)[,ttot_from2005,]# efficiency of technologies with time-dependent eta - df.eff <- as.quitte(mbind(pm_eta_conv, pm_dataeta)) %>% + df.eff <- as.quitte(mbind(pm_eta_conv, pm_dataeta[,,setdiff(getNames(pm_dataeta),getNames(pm_eta_conv))])) %>% rename(tech = all_te, eff = value) %>% select(region, period, tech, eff) @@ -970,7 +974,7 @@ df.co2price.weighted <- df.pomeg.expand %>% add_dimension(add = "tech", nm = "dac") %>% add_dimension(add = "output", nm = "cco2") %>% add_dimension(add = "type", nm = "marginal") %>% - add_dimension(add = "sector", dim=3.4, nm = "supply-side") + add_dimension(add = "sector", dim=3.4, nm = "carbon management") # Co2 Capture price, marginal of q_balcapture, convert from tr USD 2005/GtC to USD2015/tCO2 @@ -1061,6 +1065,7 @@ df.co2price.weighted <- df.pomeg.expand %>% v32_curt <- 0 } + curt_share <- v32_curt[,,teVRE] / vm_prodSe[,,teVRE] df.curtShare <- as.quitte(curt_share) %>% @@ -1090,6 +1095,7 @@ df.co2price.weighted <- df.pomeg.expand %>% pm_eff <- mbind(pm_eta_conv, pm_dataeta[,, setdiff(getNames(pm_dataeta), getNames(pm_eta_conv)) ]) vm_co2CCS_m <- pm_emifac_cco2/pm_eff[,,getNames(pm_emifac_cco2, dim=3)]*collapseNames(p_share_carbonCapture_stor) + # calculate CCS tax markup following q21_taxrevCCS, convert to USD2015/MWh CCStax <- dimReduce(pm_data_omf[,,"ccsinje"]*vm_costTeCapital[,,"ccsinje"]*vm_co2CCS_m^2/pm_dataccs[,,"quan.1"]/pm_ccsinjecrate/s_twa2mwh*1e12*1.2) @@ -1104,9 +1110,6 @@ df.co2price.weighted <- df.pomeg.expand %>% - -# mutate( CCSCost = 0) - ### Read Flexibility Tax ---- cm_FlexTax <- readGDX(gdx, "cm_flex_tax") @@ -1127,23 +1130,36 @@ df.co2price.weighted <- df.pomeg.expand %>% ### Read Final Energy Taxes ---- + + entyFe2Sector <- readGDX(gdx, "entyFe2Sector") + sector.mapping <- c("build" = "buildings", "indst" = "industry", "trans" = "transport") - pm_tau_fe_tax <- readGDX(gdx, "pm_tau_fe_tax")[,ttot_from2005,entyFe2Sector$all_enty] - pm_tau_fe_sub <- readGDX(gdx, "pm_tau_fe_sub")[,ttot_from2005,entyFe2Sector$all_enty] + pm_tau_fe_tax <- readGDX(gdx, "pm_tau_fe_tax", restore_zeros = F) + pm_tau_fe_sub <- readGDX(gdx, "pm_tau_fe_sub", restore_zeros = F) - sector.mapping <- c("build" = "buildings", "indst" = "industry", "trans" = "transport") + df.taxrate <- as.quitte(pm_tau_fe_tax * 1.2 / s_twa2mwh * 1e12) %>% + rename(taxrate = value) + df.subrate <- as.quitte(pm_tau_fe_sub * 1.2 / s_twa2mwh * 1e12) %>% + rename(subrate = value) - df.FEtax <- as.quitte((pm_tau_fe_tax+pm_tau_fe_sub)* 1.2 / s_twa2mwh * 1e12) %>% + df.FEtax <- df.taxrate %>% + full_join(df.subrate) %>% + # set NA to 0 + mutate(subrate = ifelse(is.na(subrate),0,subrate), + taxrate = ifelse(is.na(taxrate),0,taxrate)) %>% + # taxrate + subsidy rate = net FE tax + mutate( FEtax = taxrate + subrate) %>% filter(emi_sectors %in% names(sector.mapping)) %>% revalue.levels(emi_sectors = sector.mapping) %>% - rename( sector = emi_sectors, output = all_enty, FEtax = value) + rename( sector = emi_sectors, output = all_enty) %>% + select(region,period,sector,output,FEtax) -### read H2 Phase-in cost for buildings and industry ---- +### Read H2 Phase-in cost for buildings and industry ---- ### for buildings and industry (calculate average cost for simplicity, not marginal) @@ -1200,7 +1216,7 @@ df.co2price.weighted <- df.pomeg.expand %>% left_join(df.curtShare, by = c("region", "period", "tech")) %>% left_join(df.CCStax, by = c("region", "period", "tech")) %>% left_join(df.flexPriceShare, by = c("region", "period", "tech")) %>% - left_join(df.FEtax,relationship = "many-to-many", by = c("region", "period", "output")) %>% + left_join(df.FEtax, relationship = "many-to-many", by = c("region", "period", "output")) %>% left_join(df.AddTeInvH2, by = c("region", "period", "tech", "sector")) %>% # filter to only have LCOE technologies filter( tech %in% c(te_LCOE)) @@ -1208,12 +1224,12 @@ df.co2price.weighted <- df.pomeg.expand %>% # only retain unique region, period, tech combinations df.LCOE <- df.LCOE %>% - unique(by=c("region", "period", "tech")) + unique(by=c("region", "period", "tech","sector")) # replace NA by 0 in certain columns # columns where NA should be replaced by 0 - col.NA.zero <- c("OMF","OMV", "co2.price", "fuel.price","co2_dem","emiFac.se2fe","Co2.Capt.Price", + col.NA.zero <- c("OMF","OMV", "AdjCost","co2.price","co2.price.weighted", "fuel.price","fuel.price.weighted", "co2_dem","emiFac.se2fe","Co2.Capt.Price", "secfuel.prod", "secfuel.price", "curtShare","CCStax.cost","FEtax","AddH2TdCost") df.LCOE[,col.NA.zero][is.na(df.LCOE[,col.NA.zero])] <- 0 @@ -1223,18 +1239,44 @@ df.co2price.weighted <- df.pomeg.expand %>% df.LCOE[,col.NA.one][is.na(df.LCOE[,col.NA.one])] <- 1 # replace NA for sectors by "supply-side" (for all SE generating technologies) + # Carbon management for all technologies handling CO2 df.LCOE <- df.LCOE %>% - mutate( sector = ifelse(is.na(sector), "supply-side", as.character(sector))) + mutate( sector = ifelse(tech %in% c(ccs2te$all_te), + "carbon management", + sector)) %>% + mutate( sector = ifelse(tech %in% c(pe2se$all_te, + se2se$all_te), + "supply-side", + sector)) %>% + filter( !is.na(sector)) + - ### data preparation before LCOE calculation + + + ### unit conversion and data preparation before LCOE calculation df.LCOE <- df.LCOE %>% # remove some unnecessary columns select(-model, -scenario, -variable, -unit) %>% - # unit conversions for CAPEX and OMV cost + # unit conversions for CAPEX and OMV cost, adjustment cost # conversion from tr USD 2005/TW to USD2015/kW - mutate(CAPEX = CAPEX *1.2 * 1e3) %>% + # in case of CCS technologies convert from tr USD2005/GtC to USD2015/tCO2 + mutate(CAPEX = ifelse(tech %in% ccs2te$all_te, + # CCS technology unit conversion + CAPEX *1.2*1e3/3.66, + # energy technology unit conversion + CAPEX *1.2 * 1e3), + AdjCost = ifelse(tech %in% ccs2te$all_te, + # CCS technology unit conversion + AdjCost *1.2*1e3/3.66, + # energy technology unit conversion + AdjCost *1.2 * 1e3)) %>% # conversion from tr USD 2005/TWa to USD2015/MWh - mutate(OMV = OMV * 1.2 / as.numeric(s_twa2mwh) * 1e12) %>% + # in case of CCS technologies convert from tr USD2005/GtC to USD2015/tCO2 + mutate(OMV = ifelse(tech %in% ccs2te$all_te, + # CCS technology unit conversion + OMV * 1.2 * 1e3 / 3.66, + # energy technology unit conversion + OMV * 1.2 / as.numeric(s_twa2mwh) * 1e12)) %>% # share of stored carbon from captured carbon is only relevant for CCS technologies, others -> 1 mutate( CO2StoreShare = ifelse(tech %in% teCCS, CO2StoreShare, 1)) %>% # calculate annuity factor for annualizing investment cost over lifetime @@ -1245,23 +1287,39 @@ df.co2price.weighted <- df.pomeg.expand %>% df.LCOE <- df.LCOE %>% - # investment cost LCOE in USD/MWh - mutate( `Investment Cost` = CAPEX * annuity.fac / (CapFac*8760)*1e3) %>% - # OMF cost LCOE in USD/MWh - mutate( `OMF Cost` = CAPEX * OMF / (CapFac*8760)*1e3) %>% + # investment cost LCOE in USD/MWh or USD/tCO2 + mutate( `Investment Cost` = ifelse(tech %in% ccs2te$all_te, + # CCS technology, CAPEX in USD/tCO2 + CAPEX * annuity.fac / CapFac, + # energy technology, CAPEX in USD/kW(out) + CAPEX * annuity.fac / (CapFac*8760)*1e3)) %>% + # OMF cost LCOE in USD/MWh, OMF are defined as share of CAPEX + mutate( `OMF Cost` = ifelse(tech %in% ccs2te$all_te, + # CCS technology, CAPEX in USD/tCO2 + CAPEX * OMF / CapFac, + # energy technology, CAPEX in USD/kW(out) + CAPEX * OMF / (CapFac*8760)*1e3)) %>% mutate( `OMV Cost` = OMV) %>% + # adjustment cost LCOE in USD/MWh or USD/tCO2 + # (parallel to investment cost calculation) + mutate( `Adjustment Cost` = ifelse(tech %in% ccs2te$all_te, + # CCS technology, CAPEX in USD/tCO2 + AdjCost * annuity.fac / CapFac, + # energy technology, CAPEX in USD/kW(out) + AdjCost * annuity.fac / (CapFac*8760)*1e3)) %>% # Fuel Cost - # # Fuel cost with fuel price in commissioning time step of plant - # mutate( `Fuel Cost` = fuel.price / eff) %>% + # # Fuel cost with fuel price of time step for which LCOE are calculated + mutate( `Fuel Cost (time step prices)` = fuel.price / eff) %>% # Fuel cost with weighted average fuel price over plant lifetime - mutate( `Fuel Cost` = fuel.price.weighted / eff) %>% + mutate( `Fuel Cost (intertemporal prices)` = fuel.price.weighted / eff) %>% # CO2 Tax cost # CO2 Tax cost on SE level refer to (supply-side) emissions of pe2se technologies # CO2 Tax cost on FE level refer to (demand-side) emissions of FE carriers - # # CO2 Tax cost with co2 price in commissioning time step of plant - # mutate( `CO2 Tax Cost` = co2.price * (emiFac / eff) * CO2StoreShare) %>% + # # CO2 Tax cost with co2 price of time step for which LCOE are calculated + mutate( `CO2 Tax Cost (time step prices)` = co2.price * (emiFac / eff) * CO2StoreShare) %>% # CO2 Tax cost with weighted average co2 price over plant lifetime - mutate( `CO2 Tax Cost` = co2.price.weighted * (emiFac / eff) * CO2StoreShare) %>% + # These are the CO2 tax cost that feature the total LCOE calculate below + mutate( `CO2 Tax Cost (intertemporal prices)` = co2.price.weighted * (emiFac / eff) * CO2StoreShare) %>% mutate( `CO2 Provision Cost` = Co2.Capt.Price * co2_dem) %>% # fuel cost of second fuel if technology has two inputs or outputs # is positive for cost of second input @@ -1273,17 +1331,21 @@ df.co2price.weighted <- df.pomeg.expand %>% mutate( `CCS Tax Cost` = CCStax.cost) %>% # Flex Tax benefit for electrolysis # FlexPriceShare denotes share of the electricity price that electrolysis sees - mutate( `Flex Tax` = -(1-FlexPriceShare) * `Fuel Cost`) %>% + mutate( `Flex Tax` = -(1-FlexPriceShare) * `Fuel Cost (time step prices)`) %>% # se2fe technologies come with FE tax mutate( `FE Tax` = FEtax, # FE H2 has some additional t&d cost at low H2 shares (phase-in cost) in REMIND `Additional H2 t&d Cost` = AddH2TdCost) %>% # calculate total LCOE by adding all components - mutate( `Total LCOE` = `Investment Cost` + `OMF Cost` + `OMV Cost` + `Fuel Cost` + `CO2 Tax Cost` + + mutate( + # Total LCOE with fuel cost and co2 tax cost based on fuel prices and carbon prices of time step for which LCOE are calculated + `Total LCOE (time step prices)` = `Investment Cost` + `OMF Cost` + `OMV Cost` + `Adjustment Cost` + `Fuel Cost (time step prices)` + `CO2 Tax Cost (time step prices)` + + `CO2 Provision Cost` + `Second Fuel Cost` + `CCS Tax Cost` + `Curtailment Cost` + `Flex Tax` + `FE Tax` + `Additional H2 t&d Cost`, + # Total LCOE with fuel cost and co2 tax cost based on fuel prices and carbon prices that are intertemporally weighted and averaged over the plant lifetime + `Total LCOE (intertemporal prices)` = `Investment Cost` + `OMF Cost` + `OMV Cost` + `Adjustment Cost` + `Fuel Cost (intertemporal prices)` + `CO2 Tax Cost (intertemporal prices)` + `CO2 Provision Cost` + `Second Fuel Cost` + `CCS Tax Cost` + `Curtailment Cost` + `Flex Tax` + `FE Tax` + `Additional H2 t&d Cost`) - # Levelized Cost of UE in Buildings Putty Realization ---- # if detailed buildings module on, calculate LCOE of fe2ue technologies @@ -1358,11 +1420,20 @@ df.co2price.weighted <- df.pomeg.expand %>% # transform marginal LCOE to mif output format df.LCOE.out <- df.LCOE %>% - select(region, period, tech, output, sector, `Investment Cost`, `OMF Cost`, `OMV Cost`, `Fuel Cost` , - `CO2 Tax Cost`,`CO2 Provision Cost`,`Second Fuel Cost`,`VRE Storage Cost` ,`Grid Cost`, `Curtailment Cost`, - `CCS Tax Cost`, `CCS Cost`,`Flex Tax`,`FE Tax`,`Additional H2 t&d Cost`,`Total LCOE`) %>% + select(region, period, tech, output, sector, + `Investment Cost`, `Adjustment Cost`, `OMF Cost`, `OMV Cost`, + `Fuel Cost (time step prices)` , `CO2 Tax Cost (time step prices)`, + `Fuel Cost (intertemporal prices)`, `CO2 Tax Cost (intertemporal prices)`, + `CO2 Provision Cost`,`Second Fuel Cost`, `Curtailment Cost`, + `CCS Tax Cost`, `Flex Tax`,`FE Tax`,`Additional H2 t&d Cost`, + `Total LCOE (time step prices)`, + `Total LCOE (intertemporal prices)` + ) %>% gather(cost, value, -region, -period, -tech, -output, -sector) %>% - mutate(unit = "US$2015/MWh", type="marginal") %>% + mutate(unit = ifelse(tech %in% ccs2te$all_te, + "US$2015/tCO2", + "US$2015/MWh"), + type="marginal") %>% filter( value != 0) %>% select(region, period, type, output, tech, sector, unit, cost, value) From 817b21e29506e14c3c0a206c9b1cd8cdbc2f043c Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Wed, 22 Nov 2023 13:28:42 +0100 Subject: [PATCH 7/8] fix errors of buildLibrary tests --- R/reportLCOE.R | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/R/reportLCOE.R b/R/reportLCOE.R index b5ed7dc4..bd424882 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -489,6 +489,7 @@ LCOE.avg <- NULL if (output.type %in% c("marginal", "both", "marginal detail")) { # variable definitions for dplyr operations in the following section +# avoids error of "no visible binding for global variable for X" in buildLibrary() region <- NULL period <- NULL all_te <- NULL @@ -570,6 +571,23 @@ LCOE.avg <- NULL emi_sectors <- NULL fe2es.buildings <- pm_tau_fe_tax_ES_st <- pm_tau_fe_sub_ES_st <- NULL model <- scenario <- variable <- unit <- NULL + AdjCost <- NULL + operationPeriod <- NULL + discount <- NULL + weight <- NULL + p_teAnnuity <- NULL + subrate <- NULL + taxrate <- NULL + fuel.price.weighted <- NULL + co2.price.weighted <- NULL + `Fuel Cost (time step prices)` <- NULL + `Fuel Cost (intertemporal prices)` <- NULL + `CO2 Tax Cost (time step prices)` <- NULL + `CO2 Tax Cost (intertemporal prices)` <- NULL + `Total LCOE (time step prices)` <- NULL + `Total LCOE (intertemporal prices)` <- NULL + `Adjustment Cost` <- NULL + @@ -742,7 +760,8 @@ LCOE.avg <- NULL # # annuity factor from REMIND, # TODO: check whether this is the same as calculated above - # p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) + # so far only used in levelized cost of DAC calculation below + p_teAnnuity <- readGDX(gdx, "p_teAnnuity", restore_zeros = F) ### Read marginal adjustment costs ---- @@ -836,7 +855,7 @@ LCOE.avg <- NULL as.numeric(opTimeYr) + period, period)) %>% # remove time steps from operationPeriod that are not remind_timesteps - filter( operationPeriod %in% unique(remind_timesteps$period)) + filter( operationPeriod %in% unique(quitte::remind_timesteps$period)) # Join capacity distribution with above dataframe to get # capacity distribution over all commissioning years (period) and From f07f5d0f1483f4bc97cf159f163c05c24c6b57ba Mon Sep 17 00:00:00 2001 From: Felix Schreyer Date: Wed, 22 Nov 2023 13:30:15 +0100 Subject: [PATCH 8/8] increment library version number --- .buildlibrary | 2 +- CITATION.cff | 4 ++-- DESCRIPTION | 4 ++-- NAMESPACE | 1 + README.md | 6 +++--- man/reportLCOE.Rd | 2 ++ 6 files changed, 11 insertions(+), 8 deletions(-) diff --git a/.buildlibrary b/.buildlibrary index 3fcd1617..4e8ed656 100644 --- a/.buildlibrary +++ b/.buildlibrary @@ -1,4 +1,4 @@ -ValidationKey: '220922208' +ValidationKey: '221236920' AcceptedWarnings: - 'Warning: package ''.*'' was built under R version' - 'Warning: namespace ''.*'' is not available and has been replaced' diff --git a/CITATION.cff b/CITATION.cff index 4ffb5886..d5e97235 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,8 +2,8 @@ cff-version: 1.2.0 message: If you use this software, please cite it using the metadata from this file. type: software title: 'remind2: The REMIND R package (2nd generation)' -version: 1.123.2 -date-released: '2023-11-08' +version: 1.124.0 +date-released: '2023-11-22' abstract: Contains the REMIND-specific routines for data and model output manipulation. authors: - family-names: Rodrigues diff --git a/DESCRIPTION b/DESCRIPTION index 499a113d..0870041a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Type: Package Package: remind2 Title: The REMIND R package (2nd generation) -Version: 1.123.2 -Date: 2023-11-08 +Version: 1.124.0 +Date: 2023-11-22 Authors@R: c( person("Renato", "Rodrigues", , "renato.rodrigues@pik-potsdam.de", role = c("aut", "cre")), person("Lavinia", "Baumstark", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index ade6725f..23e0a152 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -262,6 +262,7 @@ importFrom(quitte,as.quitte) importFrom(quitte,calcCumulatedDiscount) importFrom(quitte,df.2.named.vector) importFrom(quitte,getColValues) +importFrom(quitte,getPeriods) importFrom(quitte,getRegs) importFrom(quitte,getScenarios) importFrom(quitte,inline.data.frame) diff --git a/README.md b/README.md index 863a65de..c9d978c8 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # The REMIND R package (2nd generation) -R package **remind2**, version **1.123.2** +R package **remind2**, version **1.124.0** [![CRAN status](https://www.r-pkg.org/badges/version/remind2)](https://cran.r-project.org/package=remind2) [![R build status](https://github.com/pik-piam/remind2/workflows/check/badge.svg)](https://github.com/pik-piam/remind2/actions) [![codecov](https://codecov.io/gh/pik-piam/remind2/branch/master/graph/badge.svg)](https://app.codecov.io/gh/pik-piam/remind2) [![r-universe](https://pik-piam.r-universe.dev/badges/remind2)](https://pik-piam.r-universe.dev/builds) @@ -49,7 +49,7 @@ In case of questions / problems please contact Renato Rodrigues . +Rodrigues R, Baumstark L, Benke F, Dietrich J, Dirnaichner A, Führlich P, Giannousakis A, Hasse R, Hilaire J, Klein D, Koch J, Kowalczyk K, Levesque A, Malik A, Merfort A, Merfort L, Morena-Leiva S, Pehl M, Pietzcker R, Rauner S, Richters O, Rottoli M, Schötz C, Schreyer F, Siala K, Sörgel B, Spahr M, Strefler J, Verpoort P, Weigmann P (2023). _remind2: The REMIND R package (2nd generation)_. R package version 1.124.0, . A BibTeX entry for LaTeX users is @@ -58,7 +58,7 @@ A BibTeX entry for LaTeX users is title = {remind2: The REMIND R package (2nd generation)}, author = {Renato Rodrigues and Lavinia Baumstark and Falk Benke and Jan Philipp Dietrich and Alois Dirnaichner and Pascal Führlich and Anastasis Giannousakis and Robin Hasse and Jérome Hilaire and David Klein and Johannes Koch and Katarzyna Kowalczyk and Antoine Levesque and Aman Malik and Anne Merfort and Leon Merfort and Simón Morena-Leiva and Michaja Pehl and Robert Pietzcker and Sebastian Rauner and Oliver Richters and Marianna Rottoli and Christof Schötz and Felix Schreyer and Kais Siala and Björn Sörgel and Mike Spahr and Jessica Strefler and Philipp Verpoort and Pascal Weigmann}, year = {2023}, - note = {R package version 1.123.2}, + note = {R package version 1.124.0}, url = {https://github.com/pik-piam/remind2}, } ``` diff --git a/man/reportLCOE.Rd b/man/reportLCOE.Rd index 3518acb1..97f0de56 100644 --- a/man/reportLCOE.Rd +++ b/man/reportLCOE.Rd @@ -23,6 +23,8 @@ It includes most technologies that generate secondary energy and the distributio This script calculates two different types of LCOE: average LCOE (standing system) and marginal LCOE (new plant). The average LCOE reflect the total cost incurred by the technology deployment in a specific time step divided by its energy output. The marginal LCOE estimate the per-unit lifetime cost of the output if the model added another capacity of that technology in the respective time step. +The marginal LCOE are calculate in two versions: 1) with fuel prices and co2 taxes of the time step for which the LCOE are calculated (time step prices), +or 2) with an intertemporal weighted-average of fuel price and co2 taxes over the lifetime of the plant (intertemporal prices). } \examples{