diff --git a/core/equations.gms b/core/equations.gms index 678b791c..c7eb26ac 100644 --- a/core/equations.gms +++ b/core/equations.gms @@ -31,7 +31,7 @@ qDummyObjPGALL(allCy,YTIME)$(TIME(YTIME) and runCy(allCy)).. 1e-5 * SQR( i04MatFacPlaAvailCap(allCy,PGALL,YTIME) - i04MatFacPlaAvailCap(allCy,PGALL,YTIME-1) ) + - 0 * SQR( + 1e-4 * SQR( i04ScaleEndogScrap(allCy,PGALL,YTIME) - i04ScaleEndogScrap(allCy,PGALL,YTIME-1) ) ) / CARD(PGALL); @@ -71,7 +71,7 @@ qDummyObjINDDOMShares(allCy,YTIME,DSBS)$(TIME(YTIME) and runCy(allCy) and (INDDO ) / SUM(EFS$SECtoEF(DSBS,EFS), 1) )$((INDSE(DSBS) or sameas("NEN",DSBS) or sameas("PCH",DSBS)) and t02FinalEnergyINDSE(allCy,DSBS,YTIME)) + SUM(ITECH$SECTTECH(DSBS,ITECH), - 1e-5 * SQR( + 1e-3 * SQR( i02ScaleEndogScrap(allCy,DSBS,ITECH,YTIME) - i02ScaleEndogScrap(allCy,DSBS,ITECH,YTIME-1) ) ) / card(ITECH); diff --git a/core/input.gms b/core/input.gms index a5731f64..2fba7c80 100644 --- a/core/input.gms +++ b/core/input.gms @@ -100,18 +100,18 @@ imDisc(runCy,"PC",YTIME) = 0.11; * FIXME: Drive the emission factors with mrprom * author=giannou parameter iCo2EmiFacAllSbs(EF) "CO2 emission factors (kgCO2/kgoe fuel burned)" / -CRO 2.76 -LGN 4.15330622, -HCL 3.941453651, +CRO 3.2 +LGN 4.2, +HCL 4.2, *SLD 4.438008647, -GSL 2.872144882, -GDO 3.068924588, -LPG 2.612562612, -KRS 2.964253636, -RFO 3.207089028, -OLQ 3.207089028, -NGS 2.336234395, -OGS 2.336234395, +GSL 3.2, +GDO 3.2, +LPG 3.2, +KRS 3.2, +RFO 3.2, +OLQ 3.2, +NGS 2.5, +OGS 2.5, BMSWAS 0/; *--- imCo2EmiFac(runCy,SBS,EF,YTIME)$(not (sameas("NEN",SBS) or sameas("PCH",SBS))) = iCo2EmiFacAllSbs(EF); @@ -177,7 +177,6 @@ imFuelPrice(runCy,TRANSE,"RFO",YTIME) = imFuelPrice(runCy,"BU","RFO",YTIME); imFuelPrice(runCy,TRANSE,"OGS",YTIME) = imFuelPrice(runCy,TRANSE,"NGS",YTIME); imFuelPrice(runCy,TRANSE,"OLQ",YTIME) = imFuelPrice(runCy,TRANSE,"GDO",YTIME); imFuelPrice(runCy,TRANSE,"H2F",YTIME) = 2 * imFuelPrice(runCy,TRANSE,"H2F",YTIME); -imFuelPrice(runCy,"PA","H2F",YTIME) = 2 * imFuelPrice(runCy,"PA","KRS",YTIME); imFuelPrice(runCy,"ICT",EFS,YTIME)$SECtoEF("ICT",EFS) = imFuelPrice(runCy,"SE",EFS,YTIME); *--- table imPriceFuelsIntBase(WEF,YTIME) "International Fuel Prices USED IN BASELINE SCENARIO ($2015/toe)" @@ -387,8 +386,8 @@ SE.TNGS 0.2244 6.8 20 0.88 SE.TOGS 0.2244 10.88 20 0.8 *SE.PGTSOL 0.86224 1.36 20 0.85 SE.TBMSWAS 0.323544 10.88 20 0.5 -SE.TELC 0.3 8.976 12 0.85 -SE.THEATPUMP 0.432 12.9254 20 1.848 +SE.TELC 0.3 8.976 12 0.97 +SE.THEATPUMP 0.432 12.9254 20 3.2 SE.TSOL 0.432 12.9254 20 1 SE.TGEO 0.432 12.9254 20 0.5 AG.THCL 0.323544 10.88 20 0.7 @@ -403,7 +402,7 @@ AG.TNGS 0.2244 6.8 20 0.88 AG.TOGS 0.2244 10.88 20 0.8 *AG.PGTSOL 0.86224 1.36 20 0.85 AG.TBMSWAS 0.323544 10.88 20 0.5 -AG.TELC 0.3 8.976 12 0.85 +AG.TELC 0.3 8.976 12 0.9 AG.THEATPUMP 0.432 12.9254 20 1.848 AG.TSOL 0.432 12.9254 20 1 AG.TGEO 0.432 12.9254 20 0.5 @@ -419,8 +418,8 @@ HOU.TNGS 0.2244 6.8 20 0.88 HOU.TOGS 0.2244 10.88 20 0.8 *HOU.PGTSOL 0.86224 1.36 20 0.85 HOU.TBMSWAS 0.323544 10.88 20 0.5 -HOU.TELC 0.3 8.976 12 0.85 -HOU.THEATPUMP 0.432 12.9254 20 1.848 +HOU.TELC 0.3 8.976 12 0.97 +HOU.THEATPUMP 0.432 12.9254 20 3.2 HOU.TSOL 0.432 12.9254 20 1 HOU.TGEO 0.432 12.9254 20 0.5 ; @@ -696,8 +695,8 @@ $offtext $ELSE.calib variable imMatrFactor(allCy,DSBS,TECH,YTIME) "Maturity factor per technology and subsector for all countries (1)"; -imMatrFactor.LO(runCy,DSBS,TECH,YTIME) = 1e-6; -imMatrFactor.UP(runCy,DSBS,TECH,YTIME) = 10; +imMatrFactor.LO(runCy,DSBS,TECH,YTIME) = 1e-2; +imMatrFactor.UP(runCy,DSBS,TECH,YTIME) = 1; imMatrFactor.L(runCy,DSBS,TECH,YTIME) = iMatrFactorData(runCy,DSBS,TECH,YTIME); imMatrFactor.FX(runCy,DSBS,TECH,YTIME)$(not (sameas(DSBS,"PC") or sameas(DSBS,"PB") or sameas(DSBS,"GU") or INDDOM(DSBS) or sameas("NEN",DSBS) or sameas("PCH",DSBS))) = iMatrFactorData(runCy,DSBS,TECH,YTIME); imMatrFactor.FX(runCy,DSBS,TECH,YTIME)$(sameas(DSBS,"AG") and not EU28(runCy)) = iMatrFactorData(runCy,DSBS,TECH,YTIME); @@ -712,6 +711,7 @@ parameters !!imFacSubsiCapCostSupply(SSBS,STECH) !!State subsidy (%) factor in technology capex (supply side) !!imGrantCapCostSupply(SSBS,STECH) !!State granting in technology capex (supply side) imCapCostTechMin(allCy,DSBS,TECH,YTIME) !!Factor for the minimum capex of a demand technology after the state subsidy +!!#UPT imCostCapTechDisc(YTIME) !!Discount rate for capital costs of power generation technologies ; $ontext @@ -820,9 +820,10 @@ imUsfEneConvSubTech(runCy,DOMSE,TECH,YTIME) = imDataDomTech(DOMSE,TECH,"USC"); imFixOMCostTech(runCy,NENSE,TECH,YTIME)= imDataNonEneSec(NENSE,TECH,"FC"); imVarCostTech(runCy,NENSE,TECH,YTIME) = imDataNonEneSec(NENSE,TECH,"VC"); imUsfEneConvSubTech(runCy,NENSE,TECH,YTIME) = imDataNonEneSec(NENSE,TECH,"USC"); -imUsfEneConvSubTech(runCy,"BU","TH2F",YTIME) = 0.7; -imUsfEneConvSubTech(runCy,"BU","TNGS",YTIME) = 0.6; +imUsfEneConvSubTech(runCy,"BU","TH2F",YTIME) = 0.8; +imUsfEneConvSubTech(runCy,"BU","TNGS",YTIME) = 0.5; imUsfEneConvSubTech(runCy,"BU","TGSL",YTIME) = 0.5; +imCapCostTech(runCy,"BU",TECH,YTIME)$SECTTECH("BU",TECH) = imCapCostTech(runCy,"GN","TGDO",YTIME); imCapCostTech(runCy,"BU","TH2F",YTIME) = 1.5 * imCapCostTech(runCy,"BU","TGDO",YTIME); *--- ** CDR @@ -868,4 +869,8 @@ imPlantEffByType(runCy,STECH,"effHeat",YTIME)$(not PGALL(STECH))= imPlantEffByTy *--- ** Conversion of GW mean power into TWh/y, depending on whether it's a leap year smGwToTwhPerYear(YTIME) = 8.76 + 0.024 $ (mod(YTIME.val,4) = 0 and mod (YTIME.val,100) <> 0); -*-- \ No newline at end of file +*-- +!!#UPT imCostCapTechDisc(YTIME) = 0; +!!#UPT imCostCapTechDisc(YTIME)$(ord(YTIME) = 20) = 0.75; +!!#UPT imCostCapTechDisc(YTIME)$(ord(YTIME) > 20 and ord(YTIME) <= 40) = 0.75 + (ord(YTIME) - 20) * (0.5 - 0.75) / (40 - 20); +!!#UPT imCostCapTechDisc(YTIME)$(ord(YTIME) > 40) = 0.5; \ No newline at end of file diff --git a/core/preloop.gms b/core/preloop.gms index ace4ebd8..126899e3 100644 --- a/core/preloop.gms +++ b/core/preloop.gms @@ -47,7 +47,7 @@ VmPriceFuelSubsecCarVal.L(runCy,SBS,EF,YTIME)$SECtoEF(SBS,EF) = 1; $IFTHEN %softLinkMAgPIE% == on VmPriceFuelSubsecCarVal.FX(runCy,SBS,"BMSWAS",YTIME)$(An(YTIME)) = iPricesMagpie(runCy,SBS,YTIME); $ENDIF -VmPriceFuelSubsecCarVal.FX(runCy,SBS,EF,YTIME)$(SECtoEF(SBS,EF) and not sameas("NUC",EF) and DATAY(YTIME)) = imFuelPrice(runCy,SBS,EF,YTIME); +VmPriceFuelSubsecCarVal.FX(runCy,SBS,EF,YTIME)$(SECtoEF(SBS,EF) and not sameas("NUC",EF) and not sameas("H2F",EF) and DATAY(YTIME)) = imFuelPrice(runCy,SBS,EF,YTIME); * Alternative fuel prices are set explicitly below instead of using ALTMAP * FIXME: VmPriceFuelSubsecCarVal (NUC/MET/ETH/BGDO) should be computed endogenously after startYear, and with mrprom before startYear * author=giannou @@ -60,6 +60,7 @@ VmPriceFuelSubsecCarVal.FX(runCy,"H2P",EF,YTIME)$(SECtoEF("H2P",EF)$DATAY(YTIME) VmPriceFuelSubsecCarVal.FX(runCy,"STEAMP",EF,YTIME)$(SECtoEF("STEAMP",EF)$DATAY(YTIME)) = imFuelPrice(runCy,"PG",EF,YTIME); VmPriceFuelSubsecCarVal.FX(runCy,SBS,"CRO",YTIME) = imFuelPrice(runCy,SBS,"CRO",YTIME); VmPriceFuelSubsecCarVal.FX(runCy,SBS,"STE",YTIME)$(SECtoEF(SBS,"STE") and DATAY(YTIME)) = imFuelPrice(runCy,"OI","ELC",YTIME); +VmPriceFuelSubsecCarVal.FX(runCy,SBS,"H2F",YTIME)$(SECtoEF(SBS,"H2F") and DATAY(YTIME)) = 1.5 * imFuelPrice(runCy,"OI","ELC",YTIME); *--- VmPriceFuelAvgSub.LO(runCy,DSBS,YTIME) = 0; VmPriceFuelAvgSub.L(runCy,DSBS,YTIME) = 1; diff --git a/modules/01_Transport/simple/equations.gms b/modules/01_Transport/simple/equations.gms index 6e0e45a1..e1b2531f 100644 --- a/modules/01_Transport/simple/equations.gms +++ b/modules/01_Transport/simple/equations.gms @@ -286,12 +286,12 @@ Q01PremScrp(allCy,TRANSE,TTECH,YTIME)$(TIME(YTIME)$SECTTECH(TRANSE,TTECH)$runCy( V01PremScrp(allCy,TRANSE,TTECH,YTIME) =E= 1 - - (V01CostFuel(allCy,TRANSE,TTECH,YTIME-1)) ** (-2) / + (V01CostFuel(allCy,TRANSE,TTECH,YTIME-1)) ** (-1) / ( - V01CostFuel(allCy,TRANSE,TTECH,YTIME-1) ** (-2) + - 0.05*i01PremScrpFac(allCy,TRANSE,TTECH,YTIME) * + V01CostFuel(allCy,TRANSE,TTECH,YTIME-1) ** (-1) + + 0.02*i01PremScrpFac(allCy,TRANSE,TTECH,YTIME) * SUM(TTECH2$(not sameas(TTECH2,TTECH) and SECTTECH(TRANSE,TTECH2)), - V01CostTranspPerMeanConsSize(allCy,TRANSE,TTECH2,YTIME-1) ** (-2) + V01CostTranspPerMeanConsSize(allCy,TRANSE,TTECH2,YTIME-1) ** (-1) ) ); diff --git a/modules/01_Transport/simple/input.gms b/modules/01_Transport/simple/input.gms index 946644c6..37cfa123 100644 --- a/modules/01_Transport/simple/input.gms +++ b/modules/01_Transport/simple/input.gms @@ -230,6 +230,7 @@ i01AvgVehCapLoadFac(runCy,TRANSE,TRANSUSE,YTIME) = i01CapDataLoadFacEachTransp(T ** Transport Sector i01TechLft(runCy,TRANSE,TTECH,YTIME) = imDataTransTech(TRANSE,TTECH,"LFT",YTIME); i01TechLft(runCy,TRANSE,TTECH,YTIME) = 20; +i01TechLft(runCy,DOMSE,"TELC",YTIME) = 20; *--- ** Industrial Sector i01TechLft(runCy,INDSE,ITECH,YTIME) = imDataIndTechnology(INDSE,ITECH,"LFT"); diff --git a/modules/02_Industry/technology/equations.gms b/modules/02_Industry/technology/equations.gms index 7c507d66..1e66b54b 100644 --- a/modules/02_Industry/technology/equations.gms +++ b/modules/02_Industry/technology/equations.gms @@ -46,13 +46,13 @@ Q02RatioRem(allCy,DSBS,ITECH,YTIME)$(TIME(YTIME)$(SECTTECH(DSBS,ITECH) and (INDD Q02PremScrpIndu(allCy,DSBS,ITECH,YTIME)$(TIME(YTIME)$(SECTTECH(DSBS,ITECH) and (INDDOM(DSBS) or NENSE(DSBS)))$runCy(allCy)).. V02PremScrpIndu(allCy,DSBS,ITECH,YTIME) =E= - 1 - (V02VarCostTech(allCy,DSBS,ITECH,YTIME-1) * i02util(allCy,DSBS,ITECH,YTIME-1) + 1e-3) ** (-2) / + 1 - (V02VarCostTech(allCy,DSBS,ITECH,YTIME-1) * i02util(allCy,DSBS,ITECH,YTIME-1) + 1e-3) ** (-1) / ( - (V02VarCostTech(allCy,DSBS,ITECH,YTIME-1) * i02util(allCy,DSBS,ITECH,YTIME-1) + 1e-3) ** (-2) + - i02ScaleEndogScrap(allCy,DSBS,ITECH,YTIME) * + (V02VarCostTech(allCy,DSBS,ITECH,YTIME-1) * i02util(allCy,DSBS,ITECH,YTIME-1) + 1e-3) ** (-1) + + 0.01 * i02ScaleEndogScrap(allCy,DSBS,ITECH,YTIME) * sum(ITECH2$(not sameas(ITECH2,ITECH) and SECTTECH(DSBS,ITECH2)), - V02CostTech(allCy,DSBS,ITECH2,YTIME-1) + V02VarCostTech(allCy,DSBS,ITECH2,YTIME-1) - ) ** (-2) + (V02CostTech(allCy,DSBS,ITECH2,YTIME-1) + 1e-3) ** (-1) + ) ); *'NEW EQUATION' - kind of --> substitutes Q02ConsRemSubEquipSubSec diff --git a/modules/02_Industry/technology/preloop.gms b/modules/02_Industry/technology/preloop.gms index a1a7523e..5a4a0b9a 100644 --- a/modules/02_Industry/technology/preloop.gms +++ b/modules/02_Industry/technology/preloop.gms @@ -66,10 +66,14 @@ V02VarCostTech.FX(runCy,DSBS,ITECH,YTIME)$(DATAY(YTIME) and (INDDOM(DSBS) or NEN ( sum(EF$ITECHtoEF(ITECH,EF), i02ShareBlend(runCy,DSBS,ITECH,EF,YTIME) * - VmPriceFuelSubsecCarVal.L(runCy,DSBS,EF,YTIME) + - imCO2CaptRateIndustry(runCy,ITECH,YTIME) * VmCstCO2SeqCsts.L(runCy,YTIME) * 1e-3 * (imCo2EmiFac(runCy,DSBS,EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) + - (1-imCO2CaptRateIndustry(runCy,ITECH,YTIME)) * 1e-3 * (imCo2EmiFac(runCy,DSBS,EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * - (sum(NAP$NAPtoALLSBS(NAP,"PG"), VmCarVal.L(runCy,NAP,YTIME))) + ( + ( + VmPriceFuelSubsecCarVal.L(runCy,DSBS,EF,YTIME) + + imCO2CaptRateIndustry(runCy,ITECH,YTIME) * VmCstCO2SeqCsts.L(runCy,YTIME) * 1e-3 * imCo2EmiFac(runCy,DSBS,EF,YTIME) + + (1-imCO2CaptRateIndustry(runCy,ITECH,YTIME)) * 1e-3 * imCo2EmiFac(runCy,DSBS,EF,YTIME) * + sum(NAP$NAPtoALLSBS(NAP,"PG"), VmCarVal.L(runCy,NAP,YTIME)) + ) + ) ) + imVarCostTech(runCy,DSBS,ITECH,YTIME) / sUnitToKUnit ) / imUsfEneConvSubTech(runCy,DSBS,ITECH,YTIME); diff --git a/modules/03_RestOfEnergy/legacy/equations.gms b/modules/03_RestOfEnergy/legacy/equations.gms index 775c1988..8086ad49 100644 --- a/modules/03_RestOfEnergy/legacy/equations.gms +++ b/modules/03_RestOfEnergy/legacy/equations.gms @@ -17,18 +17,13 @@ Q03LossesDistr(allCy,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))).. VmLossesDistr(allCy,EFS,YTIME) =E= + imRateLossesFinCons(allCy,EFS,YTIME) * ( - imRateLossesFinCons(allCy,EFS,YTIME) * - ( - SUM(DSBS,VmFinalEnergy(allCy,DSBS,EFS,YTIME)) + - V03ProdPrimary(allCy,EFS,YTIME)$sameas(EFS,"CRO") - ) - )$(not H2EF(EFS)) + + SUM(DSBS,VmFinalEnergy(allCy,DSBS,EFS,YTIME)) + + V03ProdPrimary(allCy,EFS,YTIME)$sameas(EFS,"CRO") + ) !! FIXME: Do we need to add LQD,GAS,SLD here too? - ( - 0!!VmDemTotH2(allCy,YTIME) - - !!sum(SBS$SECtoEF(SBS,"H2F"), VmDemSecH2(allCy,SBS,YTIME)) - )$H2EF(EFS); +; $ontext *' The equation calculates the refineries' capacity for a given scenario and year. @@ -65,7 +60,8 @@ Q03InpTotTransf(allCy,SSBS,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))$SECtoEF(SSBS,E ( i03InputEffSupply(allCy,SSBS,EFS,"%fBaseY%") * SUM(EFS2, V03OutTotTransf(allCy,SSBS,EFS2,YTIME)) - )$(sameas(SSBS,"GAS") or sameas(SSBS,"SLD") or sameas(SSBS,"LQD")); + )$(sameas(SSBS,"GAS") or sameas(SSBS,"SLD") or sameas(SSBS,"LQD")) + + (SUM(H2TECH$H2TECHtoFEEDSTOCK(H2TECH,EFS),VmProdH2(allCy,H2TECH,YTIME) * i05InputOverOutH2ProdFeed(allCy,H2TECH,EFS,YTIME)))$sameas("H2P",SSBS); *' The equation calculates the total transformation output for a specific energy branch in a given scenario and year. *' The result is obtained by summing the transformation outputs from different sources, including thermal power stations, District Heating Plants, @@ -175,7 +171,7 @@ Q03ConsFiEneSec(allCy,SSBS,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))).. V03ProdPrimary(allCy,EFS2,YTIME)$(not PGRENEF(EFS2)) ) )$(not sameas("H2P",SSBS)) + - VmConsFuelH2Prod(allCy,EFS,YTIME)$sameas("H2P",SSBS); + (SUM(H2TECH$H2TECHtoENERGY(H2TECH,EFS),VmProdH2(allCy,H2TECH,YTIME) * i05InputOverOutH2ProdEnergy(allCy,H2TECH,EFS,YTIME)))$sameas("H2P",SSBS); Q03FinalEnergy(allCy,DSBS,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))$(SECtoEF(DSBS,EFS))$(not sameas("ICT",DSBS))).. VmFinalEnergy(allCy,DSBS,EFS,YTIME) diff --git a/modules/03_RestOfEnergy/legacy/input.gms b/modules/03_RestOfEnergy/legacy/input.gms index 367806bf..a26fc63c 100644 --- a/modules/03_RestOfEnergy/legacy/input.gms +++ b/modules/03_RestOfEnergy/legacy/input.gms @@ -98,6 +98,7 @@ imRateLossesFinCons(runCy,EFS,YTIME) = (sum(DSBS,imFuelCons(runCy,DSBS,EFS,YTIME)) + i03PrimProd(runCy,"CRO",YTIME)$sameas("CRO",EFS)) ]$(sum(DSBS,imFuelCons(runCy,DSBS,EFS,YTIME)) + i03PrimProd(runCy,"CRO",YTIME)$sameas("CRO",EFS)); imRateLossesFinCons(runCy,EFS,YTIME)$AN(YTIME) = imRateLossesFinCons(runCy,EFS,"%fBaseY%"); +imRateLossesFinCons(runCy,"H2F",YTIME) = imRateLossesFinCons(runCy,"ELC",YTIME); *--- i03RatioPrimaryFuels(runCy,EFS,YTIME)$DATAY(YTIME) = ( diff --git a/modules/04_PowerGeneration/simple/equations.gms b/modules/04_PowerGeneration/simple/equations.gms index 5c6451bc..2a38b98a 100644 --- a/modules/04_PowerGeneration/simple/equations.gms +++ b/modules/04_PowerGeneration/simple/equations.gms @@ -53,7 +53,9 @@ Q04CostCapTech(allCy,PGALL,YTIME)$(time(YTIME) $runCy(allCy)).. V04CostCapTech(allCy,PGALL,YTIME) =E= V04CapexRESRate(allCy,PGALL,YTIME) * V04CapexFixCostPG(allCy,PGALL,YTIME) / - (i04AvailRate(allCy,PGALL,YTIME) * smGwToTwhPerYear(YTIME) * 1000); + (i04AvailRate(allCy,PGALL,YTIME) * smGwToTwhPerYear(YTIME) * 1000) + !!#UPT * (imCostCapTechDisc(YTIME)$(sameas(PGALL,"ATHBMSCCS")) + 1$(not sameas(PGALL,"ATHBMSCCS"))) + ; *' Compute the variable cost of each power plant technology for every region, *' By utilizing the gross cost, fuel prices, CO2 emission factors & capture, and plant efficiency. diff --git a/modules/04_PowerGeneration/simple/input.gms b/modules/04_PowerGeneration/simple/input.gms index dda34954..182bb74d 100644 --- a/modules/04_PowerGeneration/simple/input.gms +++ b/modules/04_PowerGeneration/simple/input.gms @@ -100,6 +100,7 @@ $ELSE.calib parameter i04MatFacPlaAvailCap(allCy,PGALL,YTIME) "Maturity factor related to plant available capacity (1)"; parameter i04ScaleEndogScrap(allCy,PGALL,YTIME) "Scale parameter for endogenous scrapping applied to the sum of full costs (1)"; i04MatFacPlaAvailCap(runCy,PGALL,YTIME) = iMatFacPlaAvailCapData(runCy,PGALL,YTIME); +i04MatFacPlaAvailCap(runCy,"ATHBMSCCS",YTIME)$(ord(YTIME) > 22) = 0.1; i04ScaleEndogScrap(runCy,PGALL,YTIME) = iScaleEndogScrapData(runCy,PGALL,YTIME); $ENDIF.calib *--- diff --git a/modules/05_Hydrogen/legacy/declarations.gms b/modules/05_Hydrogen/legacy/declarations.gms index 5e061276..f8ea8320 100644 --- a/modules/05_Hydrogen/legacy/declarations.gms +++ b/modules/05_Hydrogen/legacy/declarations.gms @@ -2,21 +2,16 @@ *' @code Variables -V05GapShareH2Tech1(allCy, H2TECH, YTIME) "Shares of H2 production technologies in new market competition 1" -V05GapShareH2Tech2(allCy, H2TECH, YTIME) "Shares of H2 production technologies in new market competition 2" +V05GapShareH2Tech(allCy, H2TECH, YTIME) "Shares of H2 production technologies in new market competition 1" V05CapScrapH2ProdTech(allCy, H2TECH, YTIME) "Decommissioning of capacity by H2 production technology" V05PremRepH2Prod(allCy, H2TECH, YTIME) "Premature replacement of H2 production technologies" V05ScrapLftH2Prod(allCy, H2TECH, YTIME) "Scrapping of equipment due to lifetime (normal scrapping)" V05DemGapH2(allCy, YTIME) "Demand for H2 to be covered by new equipment in mtoe" V05CostProdH2Tech(allCy, H2TECH, YTIME) "Hydrogen production cost per technology in US$2015 per toe of hydrogen" V05CostVarProdH2Tech(allCy, H2TECH, YTIME) "Variable cost (including fuel cost) for hydrogen production by technology in US$2015 per toe" -V05ShareCCSH2Prod(allCy, H2TECH, YTIME) "Share of CCS technology in the decision tree between CCS and no CCS" -V05ShareNoCCSH2Prod(allCy, H2TECH, YTIME) "Share of technology without CCS in the decision tree between CCS and no CCS" -V05AcceptCCSH2Tech(allCy, YTIME) "Acceptance of investment in CCS technologies" -V05CostProdCCSNoCCSH2Prod(allCy, H2TECH, YTIME) "Production cost of the composite technology with and without CCS in Euro per toe" VmCostAvgProdH2(allCy, YTIME) "Average production cost of hydrogen in Euro per toe" V05CaptRateH2(allCy,H2TECH,YTIME) - +V05UtilRate(allCy,H2TECH,YTIME) $ontext *' **Infrastructure Variables** V05H2InfrArea(allCy, YTIME) "Number of stylised areas covered by H2 infrastructure" @@ -39,28 +34,21 @@ $offtext *' **Interdependent Variables** VmDemTotH2(allCy, YTIME) "Hydrogen production requirement in Mtoe for meeting final demand" VmProdH2(allCy, H2TECH, YTIME) "Hydrogen Production by technology in Mtoe" -VmConsFuelTechH2Prod(allCy, H2TECH, EF, YTIME) "Fuel consumption by hydrogen production technology in Mtoe" -VmDemSecH2(allCy, SBS, YTIME) "Demand for H2 by sector in mtoe" +VmCapH2(allCy, H2TECH, YTIME) "Hydrogen capacity by technology in kw-output" VmCostAvgProdH2(allCy, YTIME) "Average production cost of hydrogen in Euro per toe" -VmConsFuelH2Prod(allCy, EF, YTIME) "Total fuel consumption for hydrogen production in Mtoe" ; Equations -Q05GapShareH2Tech1(allCy, H2TECH, YTIME) "Equation for calculating the shares of technologies in hydrogen gap using Weibull equations 1" -Q05GapShareH2Tech2(allCy, H2TECH, YTIME) "Equation for calculating the shares of technologies in hydrogen gap using Weibull equations 2" +Q05GapShareH2Tech(allCy, H2TECH, YTIME) "Equation for calculating the shares of technologies in hydrogen gap using Weibull equations 1" Q05CapScrapH2ProdTech(allCy, H2TECH, YTIME) "Equation for decommissioning of capacity by H2 production technology" Q05PremRepH2Prod(allCy, H2TECH, YTIME) "Equation for premature replacement of H2 production technologies" Q05ScrapLftH2Prod(allCy, H2TECH, YTIME) "Equation for scrapping of equipment due to lifetime (normal scrapping)" Q05DemGapH2(allCy, YTIME) "Equation for gap in hydrogen demand" Q05CostProdH2Tech(allCy, H2TECH, YTIME) "Equation for hydrogen production cost per technology" Q05CostVarProdH2Tech(allCy, H2TECH, YTIME) "Equation for variable cost (including fuel cost) for hydrogen production by technology in Euro per toe" -Q05ShareCCSH2Prod(allCy, H2TECH, YTIME) "Equation for share of CCS technology in the decision tree between CCS and no CCS" -Q05ShareNoCCSH2Prod(allCy, H2TECH, YTIME) "Equation for share of technology without CCS in the decision tree between CCS and no CCS" -Q05AcceptCCSH2Tech(allCy, YTIME) "Equation for acceptance in CCS technologies" -Q05ConsFuelH2Prod(allCy, EF, YTIME) "Equation for total fuel consumption for hydrogen production" -Q05CostProdCCSNoCCSH2Prod(allCy, H2TECH, YTIME) "Equation for calculating the production cost of the composite technology with and without CCS" Q05CostAvgProdH2(allCy, YTIME) "Equation for average production cost of hydrogen in Euro per toe" Q05CaptRateH2(allCy,H2TECH,YTIME) +Q05UtilRate(allCy,H2TECH,YTIME) $ontext *' **Infrastructure Equations** Q05H2InfrArea(allCy, YTIME) "Equation for infrastructure area" @@ -79,9 +67,8 @@ Q05CostTotH2(allCy, SBS, YTIME) "Equation of total hydrogen co $offtext *' **Interdependent Equations** Q05DemTotH2(allCy, YTIME) "Equation for total hydrogen demand in a country in Mtoe" -Q05ProdH2(allCy, H2TECH, YTIME) "Equation for H2 production by technology" -Q05ConsFuelTechH2Prod(allCy, H2TECH, EF, YTIME) "Equation for fuel consumption by technology for hydrogen production" -Q05DemSecH2(allCy, SBS, YTIME) "Equation for demand of H2 by sector in mtoe" +Q05ProdH2(allCy,H2TECH,YTIME) "Equation for H2 production by technology" +Q05CapH2(allCy,H2TECH,YTIME) "Equation for H2 capacity by technology" ; Scalars diff --git a/modules/05_Hydrogen/legacy/equations.gms b/modules/05_Hydrogen/legacy/equations.gms index f3a48c07..abf7e7cc 100644 --- a/modules/05_Hydrogen/legacy/equations.gms +++ b/modules/05_Hydrogen/legacy/equations.gms @@ -7,43 +7,13 @@ Q05DemTotH2(allCy,YTIME)$(TIME(YTIME)$(runCy(allCy))).. VmDemTotH2(allCy,YTIME) =E= V03ConsGrssInl(allCy,"H2F",YTIME) - VmImpNetEneBrnch(allCy,"H2F",YTIME); -$ontext - sum(SBS$SECtoEF(SBS,"H2F"), - VmDemSecH2(allCy,SBS, YTIME) / - prod(INFRTECH$H2INFRSBS(INFRTECH,SBS), - i05EffH2Transp(allCy,INFRTECH,YTIME)* - (1-i05ConsSelfH2Transp(allCy,INFRTECH,YTIME)) - ) - ) !! increase the demand due to transportation losses -$offtext - -*' This equation calculates the sectoral hydrogen demand (VmDemSecH2) for each demand subsector (DSBS), year, and region. -*' It sums up hydrogen consumption from both industrial/tertiary sectors (using VmConsFuel) and transport sectors (using VmDemFinEneTranspPerFuel), -*' ensuring each subsector receives only the relevant demand. -Q05DemSecH2(allCy,SBS,YTIME)$(TIME(YTIME)$(runCy(allCy))).. - VmDemSecH2(allCy,SBS,YTIME) - =E= - 0;!! Unnecessary variable? author = mmadianos - !!sum(INDDOM$SAMEAS(INDDOM,SBS), VmConsFuel(allCy,INDDOM,"H2F",YTIME)) + - !!sum(TRANSE$SAMEAS(TRANSE,SBS), VmDemFinEneTranspPerFuel(allCy,TRANSE,"H2F",YTIME)) + - !!VmConsFuelCDRProd(allCy,"H2F",YTIME)$sameas("DAC",SBS) + - !!VmConsFuelCDRProd(allCy,"H2F",YTIME)$sameas("EW",SBS) + - !!VmConsFuelElecProd(allCy,"H2F",YTIME)$sameas("PG",SBS); *' This equation defines the amount of hydrogen production capacity that is scrapped due to the expiration of the useful life of plants. *' It considers the remaining lifetime of hydrogen production facilities and the impact of past production gaps. Q05ScrapLftH2Prod(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. - V05ScrapLftH2Prod(allCy,H2TECH,YTIME) - =E= - (1/i05ProdLftH2(H2TECH,YTIME))$(ord(YTIME)>14+i05ProdLftH2(H2TECH,YTIME)) -$ontext - ( - V05GapShareH2Tech1(allCy,H2TECH,YTIME-i05ProdLftH2(H2TECH,YTIME)) * - V05DemGapH2(allCy,YTIME-i05ProdLftH2(H2TECH,YTIME)) / - (VmProdH2(allCy,H2TECH,YTIME-1) + 1e-6) - )$(ord(YTIME)>14+i05ProdLftH2(H2TECH,YTIME)) -$offtext -; + V05ScrapLftH2Prod(allCy,H2TECH,YTIME) + =E= + (1/i05ProdLftH2(H2TECH,YTIME))$(ord(YTIME)>14+i05ProdLftH2(H2TECH,YTIME)); *' This equation models the premature replacement of hydrogen production capacity. It adjusts for the need to replace aging *' or inefficient hydrogen production technologies before their expected end of life based on economic factors such as cost, @@ -53,16 +23,16 @@ Q05PremRepH2Prod(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))$H2TECHPM(H2TECH =E= V05CostVarProdH2Tech(allCy,H2TECH,YTIME-1)**(-i05WBLGammaH2Prod(allCy,YTIME)) / ( - iWBLPremRepH2Prod(allCy,H2TECH,YTIME) * + 0.03 * iWBLPremRepH2Prod(allCy,H2TECH,YTIME) * ( sum(H2TECH2$(not sameas(H2TECH,H2TECH2)), - V05CostProdH2Tech(allCy,H2TECH2,YTIME-1) + V05CostProdH2Tech(allCy,H2TECH2,YTIME-1) ** (-i05WBLGammaH2Prod(allCy,YTIME)) !!V05GapShareH2Tech1(allCy,H2TECH2,YTIME)* !!(1/i05AvailH2Prod(allCy,H2TECH,YTIME)* !!V05CostProdH2Tech(allCy,H2TECH2,YTIME) + !!(1-1/i05AvailH2Prod(allCy,H2TECH,YTIME)) * V05CostVarProdH2Tech(allCy,H2TECH2,YTIME)) ) - )**(-i05WBLGammaH2Prod(allCy,YTIME)) + + ) + V05CostVarProdH2Tech(allCy,H2TECH,YTIME-1)**(-i05WBLGammaH2Prod(allCy,YTIME)) ); @@ -77,23 +47,25 @@ Q05CapScrapH2ProdTech(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. *' The hydrogen demand gap equation defines the difference between the total hydrogen demand (calculated in Q05DemTotH2) and *' the actual hydrogen production capacity. It ensures that the gap value is non-negative, preventing overproduction or underproduction of hydrogen. -Q05DemGapH2(allCy, YTIME)$(TIME(YTIME)$(runCy(allCy))).. - V05DemGapH2(allCy, YTIME) - =E= - ( - VmDemTotH2(allCy,YTIME) - - sum(H2TECH, - (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * - VmProdH2(allCy,H2TECH,YTIME-1) - ) + - SQRT( SQR( - VmDemTotH2(allCy,YTIME) - - sum(H2TECH, - (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * - VmProdH2(allCy,H2TECH,YTIME-1) - ) - )) )/2 -; +Q05DemGapH2(allCy,YTIME)$(TIME(YTIME)$(runCy(allCy))).. + V05DemGapH2(allCy,YTIME) + =E= + ( + VmDemTotH2(allCy,YTIME) - + sum(H2TECH, + (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * + VmCapH2(allCy,H2TECH,YTIME-1) * + i05AvailH2Prod(allCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe + ) + + SQRT(SQR( + VmDemTotH2(allCy,YTIME) - + sum(H2TECH, + (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * + VmCapH2(allCy,H2TECH,YTIME-1) * + i05AvailH2Prod(allCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe + ) + )) + ) / 2 + 1e-6; *' This equation calculates the production costs of hydrogen, including both fixed costs (e.g., capital investment) *' and variable costs (e.g., operational expenses). The costs are typically differentiated by hydrogen production @@ -103,122 +75,89 @@ Q05CostProdH2Tech(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. =E= ( imDisc(allCy,"H2P",YTIME) * - exp(imDisc(allCy,"H2P",YTIME)* i05ProdLftH2(H2TECH,YTIME)) / + exp(imDisc(allCy,"H2P",YTIME) * i05ProdLftH2(H2TECH,YTIME)) / (exp(imDisc(allCy,"H2P",YTIME) * i05ProdLftH2(H2TECH,YTIME))-1) * ( i05CostCapH2Prod(allCy,H2TECH,YTIME) + - i05CostFOMH2Prod(allCy,H2TECH,YTIME) + - i05CostVOMH2Prod(allCy,H2TECH,YTIME) + i05CostFOMH2Prod(allCy,H2TECH,YTIME) ) + - (V04CapexFixCostPG(allCy,"PGSOL",YTIME))$sameas(H2TECH,"wes") + - (V04CapexFixCostPG(allCy,"PGAWNO",YTIME))$sameas(H2TECH,"wew") + V04CapexFixCostPG(allCy,"PGSOL",YTIME)$sameas(H2TECH,"wes") + + V04CapexFixCostPG(allCy,"PGAWNO",YTIME)$sameas(H2TECH,"wew") ) / (i05AvailH2Prod(allCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe) + V05CostVarProdH2Tech(allCy,H2TECH,YTIME); - + *' This equation models the variable costs associated with hydrogen production, factoring in fuel prices (e.g., electricity or natural gas), *' CO₂ emission costs, and the efficiency of the production technology. This helps to understand the fluctuating costs based on market conditions. Q05CostVarProdH2Tech(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. V05CostVarProdH2Tech(allCy,H2TECH,YTIME) =E= - sum(EF$H2TECHEFtoEF(H2TECH,EF), - ( - VmPriceFuelSubsecCarVal(allCy,"H2P",EF,YTIME) * 1e3 + - V05CaptRateH2(allCy,H2TECH,YTIME) * (imCo2EmiFac(allCy,"H2P",EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * VmCstCO2SeqCsts(allCy,YTIME) + - (1-V05CaptRateH2(allCy,H2TECH,YTIME)) * imCo2EmiFac(allCy,"H2P",EF,YTIME) * sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal(allCy,NAP,YTIME)) - ) - )$(not H2TECHREN(H2TECH)) / i05EffH2Prod(allCy,H2TECH,YTIME) + - (i04VarCost("PGSOL",YTIME) / (smTWhToMtoe))$(sameas(H2TECH,"wes")) + - (i04VarCost("PGAWNO",YTIME) / (smTWhToMtoe))$(sameas(H2TECH,"wew")) -; - -*' This equation models the acceptance of carbon capture and storage (CCS) technologies in hydrogen production. -*' It evaluates the economic feasibility of adding CCS to the hydrogen production process, considering cost, -*' environmental policies, and technology readiness. -Q05AcceptCCSH2Tech(allCy,YTIME)$(TIME(YTIME)$(runCy(allCy))).. - V05AcceptCCSH2Tech(allCy,YTIME) - =E= - i05WBLGammaH2Prod(allCy,YTIME)*2 + - EXP(-0.06*((sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal(allCy,NAP,YTIME -1))))) -; - -*' This equation determines the share of hydrogen produced using CCS technologies compared to those produced without CCS. -*' The share is calculated based on relative costs, technological feasibility, and policy incentives supporting CCS. -Q05ShareCCSH2Prod(allCy,H2TECH,YTIME)$(TIME(YTIME) $H2CCS(H2TECH) $(runCy(allCy))).. - V05ShareCCSH2Prod(allCy,H2TECH,YTIME) - =E= - 1.5 * - iWBLShareH2Prod(allCy,H2TECH,YTIME) * - (V05CostProdH2Tech(allCy,H2TECH,YTIME-1) + 1e-3)**(-V05AcceptCCSH2Tech(allCy,YTIME)) / - ( - 1.5 * - iWBLShareH2Prod(allCy,H2TECH,YTIME) * - (V05CostProdH2Tech(allCy,H2TECH,YTIME-1) + 1e-3)**(-V05AcceptCCSH2Tech(allCy,YTIME)) + - - sum(H2TECH2$H2CCS_NOCCS(H2TECH,H2TECH2), - - 1 * - iWBLShareH2Prod(allCy,H2TECH2,YTIME) * - (V05CostProdH2Tech(allCy,H2TECH,YTIME-1) + 1e-3)**(-V05AcceptCCSH2Tech(allCy,YTIME))) - ) -; - -*' Similar to Q05ShareCCSH2Prod, this equation models the share of hydrogen produced without CCS technologies. -*' It calculates the proportion of production from non-CCS methods like electrolysis or SMR without CO₂ capture. -Q05ShareNoCCSH2Prod(allCy,H2TECH,YTIME)$(TIME(YTIME) $H2NOCCS(H2TECH) $(runCy(allCy))).. - V05ShareNoCCSH2Prod(allCy,H2TECH,YTIME) - =E= - 1 - sum(H2TECH2$H2CCS_NOCCS(H2TECH2,H2TECH) , V05ShareCCSH2Prod(allCy,H2TECH2,YTIME) ) -; - -*' This equation computes the weighted average production cost of hydrogen, incorporating both CCS and non-CCS production methods. -*' It provides an overall cost perspective, helping to assess which production methods dominate the market based on cost-efficiency. -Q05CostProdCCSNoCCSH2Prod(allCy,H2TECH,YTIME)$(TIME(YTIME) $H2NOCCS(H2TECH) $(runCy(allCy))) .. - V05CostProdCCSNoCCSH2Prod(allCy,H2TECH,YTIME) - =E= - V05ShareNoCCSH2Prod(allCy,H2TECH,YTIME)*V05CostProdH2Tech(allCy,H2TECH,YTIME)+ - sum(H2CCS$H2CCS_NOCCS(H2CCS,H2TECH), V05ShareCCSH2Prod(allCy,H2CCS,YTIME)*V05CostProdH2Tech(allCy,H2CCS,YTIME)) -; - -*' This equation calculates the market share of different hydrogen production technologies, considering factors -*' like cost competitiveness, policy support, and fuel availability. It adjusts market shares based on technological -*' performance and shifting cost dynamics. -Q05GapShareH2Tech2(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. - V05GapShareH2Tech2(allCy,H2TECH,YTIME) - =E= ( - iWBLShareH2Prod(allCy,H2TECH,YTIME) * - ( - V05CostProdH2Tech(allCy,H2TECH,YTIME-1)$(not H2NOCCS(H2TECH)) + - V05CostProdCCSNoCCSH2Prod(allCy,H2TECH,YTIME-1)$H2NOCCS(H2TECH) - )**(-i05WBLGammaH2Prod(allCy,YTIME)) / - sum(H2TECH2$(not H2CCS(H2TECH2)), - iWBLShareH2Prod(allCy,H2TECH2,YTIME) * + sum(EF$(H2TECHtoFEEDSTOCK(H2TECH,EF) or H2TECHtoENERGY(H2TECH,EF)), ( - V05CostProdH2Tech(allCy,H2TECH2,YTIME-1)$(not H2NOCCS(H2TECH2)) + - V05CostProdCCSNoCCSH2Prod(allCy,H2TECH2,YTIME-1)$H2NOCCS(H2TECH2) - )**(-i05WBLGammaH2Prod(allCy,YTIME)) - ) - )$(not H2CCS(H2TECH)) + - sum(H2NOCCS$H2CCS_NOCCS(H2TECH,H2NOCCS), V05GapShareH2Tech2(allCy,H2NOCCS,YTIME))$H2CCS(H2TECH); + VmPriceFuelSubsecCarVal(allCy,"H2P",EF,YTIME) * 1e3 + + V05CaptRateH2(allCy,H2TECH,YTIME) * (imCo2EmiFac(allCy,"H2P",EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * VmCstCO2SeqCsts(allCy,YTIME) + + (1-V05CaptRateH2(allCy,H2TECH,YTIME)) * imCo2EmiFac(allCy,"H2P",EF,YTIME) * sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal(allCy,NAP,YTIME)) - + V05CaptRateH2(allCy,H2TECH,YTIME) * (4.17$sameas("BMSWAS", EF)) * sum(NAP$NAPtoALLSBS(NAP,"H2P"), VmCarVal(allCy,NAP,YTIME)) + ) * (i05InputOverOutH2ProdFeed(allCy,H2TECH,EF,YTIME) + i05InputOverOutH2ProdEnergy(allCy,H2TECH,EF,YTIME)) + ) + + i05CostVOMH2Prod(allCy,H2TECH,YTIME) + + (i04VarCost("PGSOL",YTIME) / smTWhToMtoe)$sameas(H2TECH,"wes") + + (i04VarCost("PGAWNO",YTIME) / smTWhToMtoe)$sameas(H2TECH,"wew") + + + SQRT(SQR( + sum(EF$(H2TECHtoFEEDSTOCK(H2TECH,EF) or H2TECHtoENERGY(H2TECH,EF)), + ( + VmPriceFuelSubsecCarVal(allCy,"H2P",EF,YTIME) * 1e3 + + V05CaptRateH2(allCy,H2TECH,YTIME) * (imCo2EmiFac(allCy,"H2P",EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * VmCstCO2SeqCsts(allCy,YTIME) + + (1-V05CaptRateH2(allCy,H2TECH,YTIME)) * imCo2EmiFac(allCy,"H2P",EF,YTIME) * sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal(allCy,NAP,YTIME)) - + V05CaptRateH2(allCy,H2TECH,YTIME) * (4.17$sameas("BMSWAS", EF)) * sum(NAP$NAPtoALLSBS(NAP,"H2P"), VmCarVal(allCy,NAP,YTIME)) + ) * (i05InputOverOutH2ProdFeed(allCy,H2TECH,EF,YTIME) + i05InputOverOutH2ProdEnergy(allCy,H2TECH,EF,YTIME)) + ) + + i05CostVOMH2Prod(allCy,H2TECH,YTIME) + + (i04VarCost("PGSOL",YTIME) / smTWhToMtoe)$sameas(H2TECH,"wes") + + (i04VarCost("PGAWNO",YTIME) / smTWhToMtoe)$sameas(H2TECH,"wew") + )) + ) / 2 + 1e-3; *' This equation further adjusts the market share of hydrogen technologies, particularly considering *' the relative competitiveness between CCS and non-CCS technologies. It helps to model the transition *' between different production technologies over time. -Q05GapShareH2Tech1(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. - V05GapShareH2Tech1(allCy,H2TECH,YTIME) +Q05GapShareH2Tech(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. + V05GapShareH2Tech(allCy,H2TECH,YTIME) =E= - V05GapShareH2Tech2(allCy,H2TECH,YTIME)$((not H2CCS(H2TECH)) $(not H2NOCCS(H2TECH))) + - V05GapShareH2Tech2(allCy,H2TECH,YTIME)*(V05ShareCCSH2Prod(allCy,H2TECH,YTIME)$H2CCS(H2TECH) + V05ShareNoCCSH2Prod(allCy,H2TECH,YTIME)$H2NOCCS(H2TECH)) -; + i05MatFacH2(allCy,H2TECH,YTIME) * + V05CostProdH2Tech(allCy,H2TECH,YTIME-1) ** (-2) / + SUM(H2TECH2, + i05MatFacH2(allCy,H2TECH2,YTIME) * + V05CostProdH2Tech(allCy,H2TECH2,YTIME-1) ** (-2) + ); + +*' This equation defines the actual hydrogen production levels, considering both scrapped capacity and the demand gap. +*' It allocates production from different technologies to meet the overall demand, adjusting for changes in capacity and technology availability. +Q05CapH2(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. + VmCapH2(allCy,H2TECH,YTIME) + =E= + (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * VmCapH2(allCy,H2TECH,YTIME-1) + + V05GapShareH2Tech(allCy,H2TECH,YTIME) * V05DemGapH2(allCy,YTIME) / (i05AvailH2Prod(allCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe); + +Q05UtilRate(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. + V05UtilRate(allCy,H2TECH,YTIME) + =E= + VmDemTotH2(allCy,YTIME) / + SUM(H2TECH2, + VmCapH2(allCy,H2TECH2,YTIME) * + i05AvailH2Prod(allCy,H2TECH2,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe + ); *' This equation defines the actual hydrogen production levels, considering both scrapped capacity and the demand gap. *' It allocates production from different technologies to meet the overall demand, adjusting for changes in capacity and technology availability. Q05ProdH2(allCy,H2TECH,YTIME)$(TIME(YTIME)$(runCy(allCy))).. VmProdH2(allCy,H2TECH,YTIME) =E= - (1-V05CapScrapH2ProdTech(allCy,H2TECH,YTIME)) * VmProdH2(allCy,H2TECH,YTIME-1) + - V05GapShareH2Tech1(allCy,H2TECH,YTIME) * V05DemGapH2(allCy,YTIME); + V05UtilRate(allCy,H2TECH,YTIME) * + VmCapH2(allCy,H2TECH,YTIME) * + i05AvailH2Prod(allCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe; *' This equation calculates the average cost of hydrogen production across all technologies in the system. *' It accounts for varying costs of different technologies (e.g., electrolysis vs. SMR) to provide an overall assessment of hydrogen production cost. @@ -231,24 +170,6 @@ Q05CostAvgProdH2(allCy,YTIME)$(TIME(YTIME)$(runCy(allCy))).. ) / sum(H2TECH,VmProdH2(allCy,H2TECH,YTIME) + 1e-6); -*' This equation calculates the fuel consumption for each hydrogen production technology, considering -*' the efficiency of the technology and the amount of fuel required for producing a unit of hydrogen. -*' It provides insight into fuel demand for hydrogen production. -Q05ConsFuelTechH2Prod(allCy,H2TECH,EF,YTIME)$(TIME(YTIME) $H2TECHEFtoEF(H2TECH,EF) $(runCy(allCy))).. - VmConsFuelTechH2Prod(allCy,H2TECH,EF,YTIME) - =E= -* VmConsFuelTechH2Prod(allCy,H2TECH,EF,YTIME-1)+ - VmProdH2(allCy,H2TECH,YTIME)/i05EffH2Prod(allCy,H2TECH,YTIME)!!- -* (VmProdH2(allCy,H2TECH,YTIME-1)/i05EffH2Prod(allCy,H2TECH,YTIME-1)) -; - -*' This equation aggregates the total fuel consumption across all hydrogen production technologies in the system, -*' summing up the fuel requirements from all sources. It helps track the total fuel demand for hydrogen production. -Q05ConsFuelH2Prod(allCy,EF,YTIME)$(TIME(YTIME)$H2PRODEF(EF)$(runCy(allCy))).. - VmConsFuelH2Prod(allCy,EF,YTIME) - =E= - sum(H2TECH$H2TECHEFtoEF(H2TECH,EF),VmConsFuelTechH2Prod(allCy,H2TECH,EF,YTIME)) -; $ontext !! diff --git a/modules/05_Hydrogen/legacy/input.gms b/modules/05_Hydrogen/legacy/input.gms index e86dbc09..8038abfc 100644 --- a/modules/05_Hydrogen/legacy/input.gms +++ b/modules/05_Hydrogen/legacy/input.gms @@ -2,7 +2,7 @@ *' @code *--- -table i05H2Production(ECONCHARHY,H2TECH,YTIME) "Data for Hydrogen production" +table i05H2Production(ECONCHARHY2,H2TECH,YTIME) "Data for Hydrogen production" $ondelim $include"./iH2Production.csv" $offdelim @@ -36,7 +36,8 @@ i05CostCapH2Prod(allCy,H2TECH,YTIME) "Capital cost of hydrogen production i05CostFOMH2Prod(allCy,H2TECH,YTIME) "Fixed operating and maintenance costs of hydrogen production technologies in US$2015 per kW output H2" i05CostVOMH2Prod(allCy,H2TECH,YTIME) "Variable operating and maintenance costs of hydrogen production technologies in US$2015 per kW output H2" i05AvailH2Prod(allCy,H2TECH,YTIME) "Availability of hydrogen production technologies" -i05EffH2Prod(allCy,H2TECH,YTIME) "Efficiency of hydrogen production technologies" +i05InputOverOutH2ProdFeed(allCy,H2TECH,EF,YTIME) +i05InputOverOutH2ProdEnergy(allCy,H2TECH,EF,YTIME) i05CostInvH2Transp(allCy,INFRTECH,YTIME) "Investment cost of infrastructure technology" !! - Turnpike pipeline in Euro per km !! - Low pressure urban pipeline in Euro per km @@ -56,6 +57,7 @@ i05EffNetH2Transp(allCy,INFRTECH,YTIME) "Total efficiency of the distributio i05CostAvgWeight(allCy,YTIME) "Weight for pricing in average cost or in marginal cost" iWBLShareH2Prod(allCy,H2TECH,YTIME) "Maturity factors for H2 technologies" iWBLPremRepH2Prod(allCy,H2TECH,YTIME) "Maturity factors for premature replacement of H2 technologies" +i05MatFacH2(allCy,H2TECH,YTIME) ; *--- iWBLShareH2Prod(runCy,H2TECH,YTIME) = iTechShareH2Prod(H2TECH,YTIME); @@ -67,8 +69,6 @@ i05ProdLftH2("wes",YTIME) = i05H2Production("LFT","weg",YTIME); i05ProdLftH2("wew",YTIME) = i05H2Production("LFT","weg",YTIME); *--- i05CaptRateH2Prod(H2TECH) = i05H2Production("CR",H2TECH,"%fBaseY%"); -i05CaptRateH2Prod("wes") = i05CaptRateH2Prod("weg"); -i05CaptRateH2Prod("wew") = i05CaptRateH2Prod("weg"); i05CaptRateH2Prod(H2TECH)$(not H2CCS(H2TECH)) = 0; *--- i05H2Adopt(runCy,"b",YTIME) = i05H2Parameters(runCy,"B"); @@ -77,6 +77,8 @@ i05H2Adopt(runCy,"mid",YTIME) = i05H2Parameters(runCy,"mid"); i05TranspLftH2(INFRTECH,YTIME) = i05H2InfrCapCosts("LFT",INFRTECH,YTIME); *--- i05CostCapH2Prod(runCy,H2TECH,YTIME) = i05H2Production("IC",H2TECH,YTIME); +i05CostCapH2Prod("CHA","cgf",YTIME) = 0.5 * i05CostCapH2Prod("CHA","cgf",YTIME); +i05CostCapH2Prod("CHA","cgs",YTIME) = 0.54 * i05CostCapH2Prod("CHA","cgs",YTIME); i05CostCapH2Prod(runCy,"wes",YTIME) = i05H2Production("IC","weg",YTIME); i05CostCapH2Prod(runCy,"wew",YTIME) = i05H2Production("IC","weg",YTIME); *--- @@ -89,12 +91,14 @@ i05CostVOMH2Prod(runCy,"wes",YTIME) = i05H2Production("VC","weg",YTIME); i05CostVOMH2Prod(runCy,"wew",YTIME) = i05H2Production("VC","weg",YTIME); *--- i05AvailH2Prod(runCy,H2TECH,YTIME) = i05H2Production("AVAIL",H2TECH,YTIME); -i05AvailH2Prod(runCy,"wes",YTIME) = min(i05AvailH2Prod(runCy,"weg",YTIME),i04AvailRate(runCy,"PGSOL",YTIME)); -i05AvailH2Prod(runCy,"wew",YTIME) = min(i05AvailH2Prod(runCy,"weg",YTIME),i04AvailRate(runCy,"PGAWNO",YTIME)); +i05AvailH2Prod(runCy,"wes",YTIME) = i04AvailRate(runCy,"PGSOL",YTIME); +i05AvailH2Prod(runCy,"wew",YTIME) = i04AvailRate(runCy,"PGAWNO",YTIME); *--- -i05EffH2Prod(runCy,H2TECH,YTIME) = i05H2Production("EFF",H2TECH,YTIME); -i05EffH2Prod(runCy,"wes",YTIME) = i05H2Production("EFF","weg",YTIME); -i05EffH2Prod(runCy,"wew",YTIME) = i05H2Production("EFF","weg",YTIME); +i05InputOverOutH2ProdFeed(runCy,H2TECH,EFS,YTIME)$H2TECHtoFEEDSTOCK(H2TECH,EFS) = 0.7 * i05H2Production("INOUT_HEAT",H2TECH,YTIME); +i05InputOverOutH2ProdEnergy(runCy,H2TECH,EFS,YTIME)$(H2TECHtoENERGY(H2TECH,EFS) and not sameas("ELC",EFS)) = 0.3 * i05H2Production("INOUT_HEAT",H2TECH,YTIME); +i05InputOverOutH2ProdEnergy(runCy,H2TECH,EFS,YTIME)$(H2TECHtoENERGY(H2TECH,EFS) and sameas("ELC",EFS)) = i05H2Production("INOUT_ELC",H2TECH,YTIME); +i05InputOverOutH2ProdFeed(runCy,"wes",EFS,YTIME)$H2TECHtoFEEDSTOCK("wes",EFS) = i05InputOverOutH2ProdEnergy(runCy,"weg","ELC",YTIME); +i05InputOverOutH2ProdFeed(runCy,"wew",EFS,YTIME)$H2TECHtoFEEDSTOCK("wew",EFS) = i05InputOverOutH2ProdEnergy(runCy,"weg","ELC",YTIME); *--- i05CostInvH2Transp(runCy,INFRTECH,YTIME) = i05H2InfrCapCosts("IC",INFRTECH,YTIME); *--- @@ -118,7 +122,7 @@ i05HabAreaCountry(runCy) = i05H2Parameters(runCy,"AREA"); *--- i05EffNetH2Transp(runCy,INFRTECH,YTIME) = i05EffH2Transp(runCy,INFRTECH,YTIME)*(1-i05ConsSelfH2Transp(runCy,INFRTECH,YTIME)); *--- -iWBLPremRepH2Prod(runCy,H2TECH,YTIME) = 0.1 ; +iWBLPremRepH2Prod(runCy,H2TECH,YTIME) = 0.01 ; *--- loop H2EFFLOOP do loop INFRTECH2$H2NETWORK(INFRTECH2,H2EFFLOOP) do @@ -130,4 +134,7 @@ i05CostAvgWeight(runCy,YTIME) = 1; loop YTIME$(An(YTIME)) do i05CostAvgWeight(runCy,YTIME) = -1/19+i05CostAvgWeight(runCy,YTIME-1); endloop; -*--- \ No newline at end of file +*--- +i05MatFacH2(allCy,H2TECH,YTIME) = 1; +i05MatFacH2(allCy,H2TECH,YTIME)$(H2CCS(H2TECH) and ord(YTIME) <= 25) = 0.5; +i05MatFacH2(allCy,H2TECH,YTIME)$(H2CCS(H2TECH) and ord(YTIME) <= 20) = 0; diff --git a/modules/05_Hydrogen/legacy/postsolve.gms b/modules/05_Hydrogen/legacy/postsolve.gms index 96537ac4..e03b5e4c 100644 --- a/modules/05_Hydrogen/legacy/postsolve.gms +++ b/modules/05_Hydrogen/legacy/postsolve.gms @@ -5,14 +5,12 @@ * Hydrogen Module *--- -VmConsFuelTechH2Prod.FX(runCyL,H2TECH,EF,YTIME)$TIME(YTIME) = VmConsFuelTechH2Prod.L(runCyL,H2TECH,EF,YTIME)$TIME(YTIME); -V05GapShareH2Tech1.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = V05GapShareH2Tech1.L(runCyL,H2TECH,YTIME)$TIME(YTIME); VmProdH2.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = VmProdH2.L(runCyL,H2TECH,YTIME)$TIME(YTIME); +VmCapH2.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = VmCapH2.L(runCyL,H2TECH,YTIME)$TIME(YTIME); V05DemGapH2.FX(runCyL,YTIME)$TIME(YTIME) = V05DemGapH2.L(runCyL,YTIME)$TIME(YTIME); VmCostAvgProdH2.FX(runCyL,YTIME)$TIME(YTIME) = VmCostAvgProdH2.L(runCyL,YTIME)$TIME(YTIME); V05CaptRateH2.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = V05CaptRateH2.L(runCyL,H2TECH,YTIME)$TIME(YTIME); V05CostProdH2Tech.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = V05CostProdH2Tech.L(runCyL,H2TECH,YTIME)$TIME(YTIME); -V05CostProdCCSNoCCSH2Prod.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = V05CostProdCCSNoCCSH2Prod.L(runCyL,H2TECH,YTIME)$TIME(YTIME); V05CostVarProdH2Tech.FX(runCyL,H2TECH,YTIME)$TIME(YTIME) = V05CostVarProdH2Tech.L(runCyL,H2TECH,YTIME)$TIME(YTIME); *V05DelivH2InfrTech.FX(runCyL,INFRTECH,YTIME)$TIME(YTIME) = V05DelivH2InfrTech.L(runCyL,INFRTECH,YTIME)$TIME(YTIME); *--- \ No newline at end of file diff --git a/modules/05_Hydrogen/legacy/preloop.gms b/modules/05_Hydrogen/legacy/preloop.gms index e3d36975..5d0154a9 100644 --- a/modules/05_Hydrogen/legacy/preloop.gms +++ b/modules/05_Hydrogen/legacy/preloop.gms @@ -9,8 +9,8 @@ *V05CostTotH2.FX(runCy,INDDOM,YTIME)$(not An(YTIME)) = imFuelPrice(runCy,INDDOM,"STE1AH2F",YTIME)$(not An(YTIME)); *display V05CostTotH2.L; *--- -V05GapShareH2Tech1.UP(runCy,H2TECH,YTIME) = 1; -V05GapShareH2Tech1.LO(runCy,H2TECH,YTIME) = 0; +V05GapShareH2Tech.UP(runCy,H2TECH,YTIME) = 1; +V05GapShareH2Tech.LO(runCy,H2TECH,YTIME) = 0; *--- V05DemGapH2.LO(runCy,YTIME) = 0; V05DemGapH2.L(runCy,YTIME) = 10; @@ -20,23 +20,23 @@ VmDemTotH2.LO(runCy,YTIME) = 0; VmDemTotH2.L(runCy,YTIME) = (i03DataGrossInlCons(runCy,"H2F","%fBaseY%") - imFuelTrade(runCy,"IMPORTS","H2F","%fBaseY%") + imFuelTrade(runCy,"EXPORTS","H2F","%fBaseY%")) + 1; VmDemTotH2.FX(runCy,YTIME)$DATAY(YTIME) = (i03DataGrossInlCons(runCy,"H2F",YTIME) - imFuelTrade(runCy,"IMPORTS","H2F",YTIME) + imFuelTrade(runCy,"EXPORTS","H2F",YTIME)); *--- -VmProdH2.LO(runCy,H2TECH, YTIME) = 0; -VmProdH2.L(runCy,H2TECH, YTIME) = 0.5; -VmProdH2.FX(runCy,H2TECH, YTIME)$DATAY(YTIME) = 0; +VmProdH2.LO(runCy,H2TECH,YTIME) = 0; +VmProdH2.L(runCy,H2TECH,YTIME) = 1; +VmProdH2.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = 0; *--- -*VmConsFuelTechH2Prod.L(runCy,H2TECH,EF,YTIME)$(not An(YTIME)$H2TECHEFtoEF(H2TECH,EF)) = 0; -VmConsFuelTechH2Prod.FX(runCy,H2TECH,EF,"%fBaseY%")$(H2TECHEFtoEF(H2TECH,EF)) = (VmProdH2.L(runCy,H2TECH,"%fBaseY%")/i05EffH2Prod(runCy,H2TECH,"%fBaseY%")); -display i05EffH2Prod; -display VmConsFuelTechH2Prod.L; +V05UtilRate.LO(runCy,H2TECH,YTIME) = 0; +V05UtilRate.L(runCy,H2TECH,YTIME) = 1; +V05UtilRate.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = 1; +*--- +VmCapH2.LO(runCy,H2TECH,YTIME) = 0; +VmCapH2.L(runCy,H2TECH,YTIME) = 1; +VmCapH2.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = 0; *--- *V05DelivH2InfrTech.L(runCy,INFRTECH,YTIME) = 2; *V05DelivH2InfrTech.FX(runCy,INFRTECH,YTIME)$(not An(YTIME)) = 1e-5; *V05DelivH2InfrTech.FX(runCy,INFRTECH,"%fBaseY%") = 0; *display V05DelivH2InfrTech.L; *--- -V05GapShareH2Tech2.LO(runCy,H2TECH,YTIME) = 0; -V05GapShareH2Tech2.UP(runCy,H2TECH,YTIME) = 1; -*--- V05CapScrapH2ProdTech.LO(runCy,H2TECH,YTIME) = 0; V05CapScrapH2ProdTech.UP(runCy,H2TECH,YTIME) = 1; *--- @@ -50,12 +50,15 @@ V05ScrapLftH2Prod.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = 1/i05ProdLftH2(H2TECH,YT V05CostVarProdH2Tech.LO(runCy,H2TECH,YTIME) = 0; V05CostVarProdH2Tech.L(runCy,H2TECH,YTIME) = 2; V05CostVarProdH2Tech.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = -sum(EF$H2TECHEFtoEF(H2TECH,EF), - VmPriceFuelSubsecCarVal.L(runCy,"H2P",EF,YTIME) * 1e3 + - V05CaptRateH2.L(runCy,H2TECH,YTIME) * (imCo2EmiFac(runCy,"H2P",EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * VmCstCO2SeqCsts.L(runCy,YTIME) + - (1-V05CaptRateH2.L(runCy,H2TECH,YTIME)) * (imCo2EmiFac(runCy,"H2P",EF,YTIME)) * - sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal.L(runCy,NAP,YTIME)) -)$(not H2TECHREN(H2TECH)) / i05EffH2Prod(runCy,H2TECH,YTIME) + +sum(EF$(H2TECHtoFEEDSTOCK(H2TECH,EF) or H2TECHtoENERGY(H2TECH,EF)), + ( + VmPriceFuelSubsecCarVal.L(runCy,"H2P",EF,YTIME) * 1e3 + + V05CaptRateH2.L(runCy,H2TECH,YTIME) * (imCo2EmiFac(runCy,"H2P",EF,YTIME) + 4.17$(sameas("BMSWAS", EF))) * VmCstCO2SeqCsts.L(runCy,YTIME) + + (1-V05CaptRateH2.L(runCy,H2TECH,YTIME)) * (imCo2EmiFac(runCy,"H2P",EF,YTIME)) * + sum(NAP$NAPtoALLSBS(NAP,"H2P"),VmCarVal.L(runCy,NAP,YTIME)) + ) * (i05InputOverOutH2ProdFeed(runCy,H2TECH,EF,YTIME) + i05InputOverOutH2ProdEnergy(runCy,H2TECH,EF,YTIME)) +) + +i05CostVOMH2Prod(runCy,H2TECH,YTIME) + (i04VarCost("PGSOL",YTIME) / (smTWhToMtoe))$(sameas(H2TECH,"wes")) + (i04VarCost("PGAWNO",YTIME) / (smTWhToMtoe))$(sameas(H2TECH,"wew")); *--- @@ -68,8 +71,7 @@ V05CostProdH2Tech.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = (exp(imDisc(runCy,"H2P",YTIME) * i05ProdLftH2(H2TECH,YTIME))-1) * ( i05CostCapH2Prod(runCy,H2TECH,YTIME) + - i05CostFOMH2Prod(runCy,H2TECH,YTIME) + - i05CostVOMH2Prod(runCy,H2TECH,YTIME) + i05CostFOMH2Prod(runCy,H2TECH,YTIME) ) + V04CapexFixCostPG.L(runCy,"PGSOL",YTIME)$sameas(H2TECH,"wes") + V04CapexFixCostPG.L(runCy,"PGAWNO",YTIME)$sameas(H2TECH,"wew") @@ -78,27 +80,9 @@ V05CostProdH2Tech.FX(runCy,H2TECH,YTIME)$DATAY(YTIME) = (i05AvailH2Prod(runCy,H2TECH,YTIME) * smGwToTwhPerYear(YTIME) * smTWhToMtoe) + V05CostVarProdH2Tech.L(runCy,H2TECH,YTIME); *--- -V05ShareCCSH2Prod.LO(runCy,H2TECH,YTIME) = 0; -V05ShareCCSH2Prod.UP(runCy,H2TECH,YTIME) = 1; -*--- -V05ShareNoCCSH2Prod.LO(runCy,H2TECH,YTIME) = 0; -V05ShareNoCCSH2Prod.UP(runCy,H2TECH,YTIME) = 1; -*--- -VmConsFuelH2Prod.FX(runCy,EF,YTIME)$(DATAY(YTIME) and H2PRODEF(EF)) = sum(H2TECH$H2TECHEFtoEF(H2TECH,EF),VmConsFuelTechH2Prod.L(runCy,H2TECH,EF,YTIME)); -VmConsFuelH2Prod.FX(runCy,EF,YTIME)$(not H2PRODEF(EF)) = 0; -*--- -V05CostProdCCSNoCCSH2Prod.LO(runCy,H2TECH,YTIME) = epsilon6; -V05CostProdCCSNoCCSH2Prod.L(runCy,H2TECH,YTIME) = 2; -*--- VmCostAvgProdH2.LO(runCy,YTIME) = 0; VmCostAvgProdH2.L(runCy,YTIME) = 1; -VmCostAvgProdH2.FX(runCy,YTIME)$DATAY(YTIME) = -sum(H2TECH, - (VmProdH2.L(runCy,H2TECH,YTIME) + 1e-6) * - V05CostProdH2Tech.L(runCy,H2TECH,YTIME) -) / -sum(H2TECH,VmProdH2.L(runCy,H2TECH,YTIME) + 1e-6); - +VmCostAvgProdH2.FX(runCy,YTIME)$DATAY(YTIME) = V05CostProdH2Tech.L(runCy,"gsr",YTIME); *--- *V05InvNewReqH2Infra.L(runCy,INFRTECH,YTIME) = 2; diff --git a/modules/05_Hydrogen/legacy/sets.gms b/modules/05_Hydrogen/legacy/sets.gms index 2ef68bdb..a6f0f525 100644 --- a/modules/05_Hydrogen/legacy/sets.gms +++ b/modules/05_Hydrogen/legacy/sets.gms @@ -90,18 +90,25 @@ LPIPU SSGG / *--- -H2TECHEFtoEF(H2TECH,EF) "Mapping between production technologies and fuels" +H2TECHtoFEEDSTOCK(H2TECH,EF) "Mapping between production technologies and feedstock fuels for processes" / (gsr,gss).ngs !! ,smr (cgf,cgs).hcl (bgfls,bgfl).BMSWAS !! bpy,bgfs, *sht.SOL *(nht,wen).NUC -weg.ELC -wes.ELC -wew.ELC +wes.SOL +wew.WND *(opo,ops).RFO / + +H2TECHtoENERGY(H2TECH,EF) "Mapping between production technologies and fuels for energy use (combustion)" +/ +(gsr,gss).(ngs,elc) +(cgf,cgs).(hcl,elc) +(bgfls,bgfl).(BMSWAS,elc) +weg.ELC +/ *--- $ontext H2TECHtoPGALL(H2TECH,PGALL) "Mapping between hydrogen production technologies and power generation technologies used for water electrolysis" @@ -200,6 +207,27 @@ B mid CR / + +ECONCHARHY2 "Technical - Economic characteristics for demand technologies Hydrogen" +/ +IC +FC +VC +INOUT_ELC +INOUT_HEAT +SELF +AVAIL +LFT +H2KMTOE +mpips +lpipu +mpipu +AREA +MAXAREA +B +mid +CR +/ *--- INFRTECHLAB(INFRTECH,ECONCHARHY) / diff --git a/modules/06_CO2/legacy/declarations.gms b/modules/06_CO2/legacy/declarations.gms index ecd33962..b80e99fb 100644 --- a/modules/06_CO2/legacy/declarations.gms +++ b/modules/06_CO2/legacy/declarations.gms @@ -44,6 +44,6 @@ Scalars *' Proposed values for S06EmissPercCDR between 0.005 - 0.02. *' S06CapFacMinNewCDR is responsible for the minimum deployment of CDR technologies. In ambitious scenarios, this reflects the post-net-zero phase. *' Proposed values for S06CapFacMinNewCDR between 0.005 - 0.025. -S06EmissPercCDR "The percentage of emissions that needs to be captured by new CDR equipment" /0.003/ -S06CapFacMinNewCDR "The minimum level of CDR capacity expansion as a percentage of last year's capacity" /0.015/ +S06EmissPercCDR "The percentage of emissions that needs to be captured by new CDR equipment" /0.004/ +S06CapFacMinNewCDR "The minimum level of CDR capacity expansion as a percentage of last year's capacity" /0.03/ ; \ No newline at end of file diff --git a/modules/06_CO2/legacy/equations.gms b/modules/06_CO2/legacy/equations.gms index eb2b2b84..a37bc6b0 100644 --- a/modules/06_CO2/legacy/equations.gms +++ b/modules/06_CO2/legacy/equations.gms @@ -20,8 +20,9 @@ Q06CO2CaptureCCS(allCy,SBS,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))$SECtoEF(SBS,EF imPlantEffByType(allCy,PGALL,"effELC",YTIME) * V04CO2CaptRate(allCy,PGALL,YTIME) )$sameas("PG", SBS) + - sum(H2TECH$H2TECHEFtoEF(H2TECH,EFS), - VmConsFuelTechH2Prod(allCy,H2TECH,EFS,YTIME) * + sum(H2TECH$(H2TECHtoFEEDSTOCK(H2TECH,EFS) or H2TECHtoENERGY(H2TECH,EFS)), + VmProdH2(allCy,H2TECH,YTIME) * + (i05InputOverOutH2ProdFeed(allCy,H2TECH,EFS,YTIME) + i05InputOverOutH2ProdEnergy(allCy,H2TECH,EFS,YTIME)) * V05CaptRateH2(allCy,H2TECH,YTIME) )$sameas("H2P", SBS) + sum(DSBS$sameas(DSBS,SBS), diff --git a/modules/07_Emissions/legacy/declarations.gms b/modules/07_Emissions/legacy/declarations.gms index de8944ef..a5e292bb 100644 --- a/modules/07_Emissions/legacy/declarations.gms +++ b/modules/07_Emissions/legacy/declarations.gms @@ -20,6 +20,7 @@ Q07CostAbateBySrcRegTim(E07SrcMacAbate, allCy, YTIME) "Calculate total abate Q07EmiActBySrcRegTim(E07SrcMacAbate, allCy, YTIME) "Calculate remaining actual emissions" Q07GrossEmissCO2Demand(allCy,DSBS,YTIME) "Calculate gross emissions of demand subsectors" Q07EmissionsNet(allCy,YTIME) "Calculate net emissions after abatement" +Q07GrossEmissCO2Processes(allCy,SSBS,YTIME) ; Variables @@ -33,4 +34,5 @@ V07EmiActBySrcRegTim(E07SrcMacAbate,allCy,YTIME) "Actual emissions" V07CostAbateBySrcRegTim(E07SrcMacAbate,allCy,YTIME) "Total abatement cost" V07GrossEmissCO2Demand(allCy,DSBS,YTIME) "Gross emissions of demand subsectors" V07EmissionsNet(allCy,YTIME) "Net emissions after abatement" +V07GrossEmissCO2Processes(allCy,SSBS,YTIME) *; \ No newline at end of file diff --git a/modules/07_Emissions/legacy/equations.gms b/modules/07_Emissions/legacy/equations.gms index a82e49a9..6ac08480 100644 --- a/modules/07_Emissions/legacy/equations.gms +++ b/modules/07_Emissions/legacy/equations.gms @@ -74,16 +74,24 @@ Q07GrossEmissCO2Supply(allCy,SSBS,YTIME)$(TIME(YTIME)$runCy(allCy)).. =E= SUM(EFS, ( - V03InpTotTransf(allCy,SSBS,EFS,YTIME)$SSBSEMIT(SSBS) + + V03InpTotTransf(allCy,SSBS,EFS,YTIME)$SSBSENERGY(SSBS) + VmConsFiEneSec(allCy,SSBS,EFS,YTIME) ) * imCo2EmiFac(allCy,SSBS,EFS,YTIME) ); - + +Q07GrossEmissCO2Processes(allCy,SSBS,YTIME)$(TIME(YTIME)$runCy(allCy)).. + V07GrossEmissCO2Processes(allCy,SSBS,YTIME) + =E= + SUM(EFS, + V03InpTotTransf(allCy,SSBS,EFS,YTIME) * + imCo2EmiFac(allCy,SSBS,EFS,YTIME) + )$sameas(SSBS,"H2P"); + *' This equation calculates the total absolute abatement of non-CO2 emissions for a specific source, country, and time period. *' The determination is based on the Marginal Abatement Cost (MAC) curves, the exogenous carbon price, and specific unit conversion factors. The equation *' identifies the maximum abatement potential by scanning the MAC curve steps and selecting the highest reduction level where the implementation cost is less than or *' equal to the adjusted carbon price. This ensures that the model adopts all abatement measures that are economically viable given the current carbon price. -Q07RedAbsBySrcRegTim(E07SrcMacAbate, allCy, YTIME)$(TIME(YTIME)$(runCy(allCy))).. +Q07RedAbsBySrcRegTim(E07SrcMacAbate, allCy, YTIME)$(TIME(YTIME)$(runCy(allCy))$(ord(YTIME) > 17)).. V07RedAbsBySrcRegTim(E07SrcMacAbate, allCy, YTIME) =E= smax(E07MAC$(p07MacCost(E07MAC) <= iCarbValYrExog(allCy, YTIME) * p07UnitConvFactor(E07SrcMacAbate)), diff --git a/modules/07_Emissions/legacy/preloop.gms b/modules/07_Emissions/legacy/preloop.gms index 411b7171..a7a9567d 100644 --- a/modules/07_Emissions/legacy/preloop.gms +++ b/modules/07_Emissions/legacy/preloop.gms @@ -7,23 +7,10 @@ V07GrossEmissCO2Supply.LO(runCy,SSBS,YTIME) = 0; V07GrossEmissCO2Supply.FX(runCy,"H2INFR",YTIME) = 0; V07GrossEmissCO2Supply.FX(runCy,SSBS,YTIME)$DATAY(YTIME) = SUM(EFS, - ( - (-i03InpTotTransfProcess(runCy,SSBS,EFS,YTIME))$SSBSEMIT(SSBS) + - i03DataOwnConsEne(runCy,SSBS,EFS,YTIME) - - SUM(CCS$PGALLtoEF(CCS,EFS), - SUM(PGEF$sameas(PGEF,EFS), - i04ShareFuels(runCy,CCS,PGEF)) * - VmProdElec.L(runCy,CCS,YTIME) * smTWhToMtoe / - imPlantEffByType(runCy,CCS,"effELC",YTIME) * - V04CO2CaptRate.L(runCy,CCS,YTIME) - )$sameas("PG",SSBS) - - SUM(H2CCS$H2TECHEFtoEF(H2CCS,EFS), - VmProdH2.L(runCy,H2CCS,YTIME) / - i05EffH2Prod(runCy,H2CCS,YTIME) * - V05CaptRateH2.L(runCy,H2CCS,YTIME) - )$sameas("H2P",SSBS) - ) * - imCo2EmiFac(runCy,"PG",EFS,YTIME) + ( + V03InpTotTransf.L(runCy,SSBS,EFS,YTIME)$SSBSENERGY(SSBS) + + VmConsFiEneSec.L(runCy,SSBS,EFS,YTIME) + ) * imCo2EmiFac(runCy,SSBS,EFS,YTIME) ); *--- V07GrossEmissCO2Demand.LO(runCy,DSBS,YTIME) = 0; diff --git a/modules/07_Emissions/legacy/sets.gms b/modules/07_Emissions/legacy/sets.gms index 52768778..127cf2c4 100644 --- a/modules/07_Emissions/legacy/sets.gms +++ b/modules/07_Emissions/legacy/sets.gms @@ -3,7 +3,7 @@ sets *--- -SSBSEMIT(SSBS) Supply Subsectors emitting inputs /PG,H2P,CHP,STEAMP/ +SSBSENERGY(SSBS) Supply Subsectors with energy processes /PG,CHP,STEAMP/ E07MAC "Cost categories for Marginal abatement costs curves (MACC) -2010$/tC for CH4,N20 and 2005$/tC for F-gases" / 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, diff --git a/modules/08_Prices/legacy/declarations.gms b/modules/08_Prices/legacy/declarations.gms index 7329ea78..86681daf 100644 --- a/modules/08_Prices/legacy/declarations.gms +++ b/modules/08_Prices/legacy/declarations.gms @@ -18,6 +18,7 @@ Q08PriceFuelSubsecCarVal(allCy,SBS,EF,YTIME) "Compute fuel prices Q08PriceFuelAvgSub(allCy,DSBS,YTIME) "Compute average fuel price per subsector" *Q08PriceFuelSubsecCHP(allCy,DSBS,EF,YTIME) "Compute fuel prices per subsector and fuel especially for chp plants" Q08PriceElecInd(allCy,TCHP,YTIME) "Compute electricity industry prices" +Q08PriceCarbon(allCy,SBS,EFS,YTIME) ; Parameters @@ -41,6 +42,7 @@ VmPriceFuelSubsecCarVal(allCy,SBS,EF,YTIME) "Fuel prices per subs VmPriceFuelAvgSub(allCy,DSBS,YTIME) "Average fuel prices per subsector (k$2015/toe)" * VmPriceFuelSubsecCHP(allCy,DSBS,EF,YTIME) "Fuel prices per subsector and fuel for CHP plants (kUS$2015/toe)" VmPriceElecInd(allCy,TCHP,YTIME) "Electricity index - a function of industry price (1)" +VmPriceCarbon(allCy,SBS,EFS,YTIME) *' *** Miscellaneous *V08FuelPriSubNoCarb(allCy,SBS,EF,YTIME) "Fuel prices per subsector and fuel without carbon value (kUS$2015/toe)" diff --git a/modules/08_Prices/legacy/equations.gms b/modules/08_Prices/legacy/equations.gms index f2b160b3..28ee829a 100644 --- a/modules/08_Prices/legacy/equations.gms +++ b/modules/08_Prices/legacy/equations.gms @@ -52,6 +52,7 @@ $ENDIF.magpieQuantityEquation Q08BmswasPriceFactor(allCy,YTIME)$(TIME(YTIME) $runCy(allCy)).. V08BmswasPriceFactor(allCy,YTIME) =E= +( $IFTHEN.mode %bmswasPriceMode% == curve $IFTHEN.emulatorCurve %landUseEmulator% == globiom ( 1e-3 + sum(activeGlobiomScen, i08BmswasSupplyCoefGlobiom(activeGlobiomScen,allCy,"a",YTIME)) @@ -78,9 +79,12 @@ $ENDIF.emulatorCurve $ELSEIF.mode %bmswasPriceMode% == softfx VmPriceFuelSubsecCarVal(allCy,"PG","BMSWAS",YTIME) / VmPriceFuelSubsecCarVal(allCy,"PG","BMSWAS",YTIME-1) $ELSE.mode - 1 +1 $ENDIF.mode - ; +) * +EXP(1.5 * (SUM(runCy2,V03ProdPrimary(runCy2,"BMSWAS",YTIME-1)) / (140 * 23.88458966275)) ** 3) / +EXP(1.5 * (SUM(runCy2,V03ProdPrimary(runCy2,"BMSWAS",YTIME-2)) / (140 * 23.88458966275)) ** 3) +; Q08PriceFuelSubsecCarVal(allCy,SBS,EFS,YTIME)$(SECtoEF(SBS,EFS) $(not sameas("CRO",EFS)) $TIME(YTIME) $IFTHEN %softLinkMAgPIE% == on @@ -98,17 +102,13 @@ $ENDIF.magpiePriceDomain *' crude-oil pass-through, and the change in carbon cost. VmPriceFuelSubsecCarVal(allCy,SBS,EFS,YTIME-1) * (1 + (VmCostPowGenAvgLng(allCy,YTIME-1) / VmCostPowGenAvgLng(allCy,YTIME-2) - 1)$sameas("ELC",EFS)) * - (1 + (VmCostAvgProdH2(allCy,YTIME-1) / VmCostAvgProdH2(allCy,YTIME-2) - 1)$sameas("H2F",EFS)) * + (1 + ((VmCostAvgProdH2(allCy,YTIME-1) / VmCostAvgProdH2(allCy,YTIME-2)) ** 0.7 - 1)$sameas("H2F",EFS)) * (1 + (VmCostAvgProdSte(allCy,YTIME-1) / VmCostAvgProdSte(allCy,YTIME-2) - 1)$sameas("STE",EFS)) * (1 + (V08BmswasPriceFactor(allCy,YTIME) ** i08PriceTransElast(EFS,"BMSWAS") - 1)$(BIOFUELS(EFS) or sameas("BMSWAS",EFS))) * (1 + ((VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME) / VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME-1)) ** i08PriceTransElast(EFS,"CRO") - 1)$sameas("NGS",EFS)) * (1 + ((VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME) / VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME-1)) ** i08PriceTransElast(EFS,"CRO") - 1)$SECtoEFPROD("LQD",EFS)) * (1 + ((VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME) / VmPriceFuelSubsecCarVal(allCy,SBS,"CRO",YTIME-1)) ** i08PriceTransElast(EFS,"CRO") - 1)$(sameas("HCL",EFS) or sameas("LGN",EFS))) + - 1e-3 * ( - VmCarVal(allCy,"TRADE",YTIME) * imCo2EmiFac(allCy,SBS,EFS,YTIME) - - VmCarVal(allCy,"TRADE",YTIME-1) * imCo2EmiFac(allCy,SBS,EFS,YTIME-1) - )$DSBS(SBS) - ; + VmPriceCarbon(allCy,SBS,EFS,YTIME) - VmPriceCarbon(allCy,SBS,EFS,YTIME-1); $IFTHEN.magpiePriceEquation "%bmswasPriceMode%" == "curve" $IFTHEN.magpiePriceEquationSource "%landUseEmulator%" == "magpie" @@ -132,6 +132,13 @@ V08PriceFuelSepCarbonWght(allCy,DSBS,EF,YTIME) SUM(EFS2$SECtoEF(DSBS,EFS2), (VmFinalEnergy(allCy,DSBS,EFS2,YTIME) - V02FinalElecNonSubIndTert(allCy,DSBS,YTIME)$ELCEF(EFS2)) + 1e-6) ); +Q08PriceCarbon(allCy,SBS,EFS,YTIME)$(TIME(YTIME)$(runCy(allCy))).. + VmPriceCarbon(allCy,SBS,EFS,YTIME) + =E= + 1e-3 * ( + VmCarVal(allCy,"TRADE",YTIME)$(INDSE1(SBS) or ((DOMSE1(SBS) or TRANS1(SBS) or sameas("BU", SBS)) and ord(YTIME) > 17)) + ) * imCo2EmiFac(allCy,SBS,EFS,YTIME); + *' The equation calculates the average fuel price per subsector. These average prices are used to further compute electricity prices in industry *' (using the OI "other industry" avg price), as well as the aggregate fuel demand (of substitutable fuels) per subsector. *' In the transport sector they feed into the calculation of the activity levels. diff --git a/modules/08_Prices/legacy/input.gms b/modules/08_Prices/legacy/input.gms index 402a5e78..f9451840 100644 --- a/modules/08_Prices/legacy/input.gms +++ b/modules/08_Prices/legacy/input.gms @@ -12,6 +12,7 @@ $ondelim $include "iPriceTransElast.csv" $offdelim ; +i08PriceTransElast(EFS,"CRO")$SECtoEFPROD("LQD",EFS) = 0.65; *--- $IFTHEN %softLinkMAgPIE% == on table iPricesMagpie(allCy,SBS,YTIME) "Prices of biomass per subsector (k$2015/toe)" diff --git a/modules/08_Prices/legacy/postsolve.gms b/modules/08_Prices/legacy/postsolve.gms index 49c545be..aa85fa27 100644 --- a/modules/08_Prices/legacy/postsolve.gms +++ b/modules/08_Prices/legacy/postsolve.gms @@ -8,6 +8,7 @@ VmPriceFuelAvgSub.FX(runCyL,DSBS,YTIME)$TIME(YTIME) = VmPriceFuelAvgSub.L(runCyL VmPriceFuelSubsecCarVal.FX(runCyL,SBS,EF,YTIME)$TIME(YTIME) = VmPriceFuelSubsecCarVal.L(runCyL,SBS,EF,YTIME)$TIME(YTIME); VmPriceElecInd.FX(runCyL,TCHP,YTIME)$TIME(YTIME) = VmPriceElecInd.L(runCyL,TCHP,YTIME)$TIME(YTIME); V08PriceFuelSepCarbonWght.FX(runCyL,DSBS,EF,YTIME)$TIME(YTIME) = V08PriceFuelSepCarbonWght.L(runCyL,DSBS,EF,YTIME)$TIME(YTIME); +VmPriceCarbon.FX(runCyL,SBS,EFS,YTIME)$TIME(YTIME) = VmPriceCarbon.L(runCyL,SBS,EFS,YTIME)$TIME(YTIME); *--- *' Land-use emulator emission accounting (landEmiMode == curve only) *' diff --git a/modules/08_Prices/legacy/preloop.gms b/modules/08_Prices/legacy/preloop.gms index 95509505..80a323c9 100644 --- a/modules/08_Prices/legacy/preloop.gms +++ b/modules/08_Prices/legacy/preloop.gms @@ -20,6 +20,9 @@ $offtext V08BmswasPriceFactor.LO(runCy,YTIME) = 0; V08BmswasPriceFactor.L(runCy,YTIME) = 1; *--- +VmPriceCarbon.LO(runCy,SBS,EFS,YTIME) = 0; +VmPriceCarbon.FX(runCy,SBS,EFS,YTIME)$DATAY(YTIME) = 1e-3 * iCarbValYrExog(runCy,YTIME)$INDSE1(SBS) * imCo2EmiFac(runCy,SBS,EFS,YTIME); +*--- $IFTHEN %landEmiMode% == curve * Both emulator backends use the same native-MAgPIE AFOLU history on DATAY. * TIME values are calculated by the selected backend in postsolve. diff --git a/scripts/tasks/findCarbonPrice.R b/scripts/tasks/findCarbonPrice.R index 9e4e1279..0d6ef401 100644 --- a/scripts/tasks/findCarbonPrice.R +++ b/scripts/tasks/findCarbonPrice.R @@ -76,7 +76,9 @@ suppressPackageStartupMessages({ inputCsvPath <- "data/iEnvPolicies.csv" outputCsvPath <- "data/iEnvPolicies_updated.csv" backupCsvPath <- "data/iEnvPolicies_backup.csv" +lastTestedCsvPath <- "data/iEnvPolicies_last_tested.csv" # last policy tested before a failure workDir <- getwd() +lastTestedPolicy <- NULL # most recent policy table sent to OPEN-PROM # One-time backup of the original canonical file if (file.exists(inputCsvPath) && !file.exists(backupCsvPath)) { @@ -100,10 +102,10 @@ readEnvPolicies <- function(csvPath = inputCsvPath) { list(envWide = envWide, envLong = envLong, yearCols = yearCols) } -applyAlpha <- function(envWide, yearCols, alpha, targetRegion) { +applyAlpha <- function(envWide, yearCols, alpha, targetRegion, fromYear = changeCarbonPriceFromYear) { x <- data.table::copy(envWide) - yearColsFuture <- yearCols[as.integer(yearCols) >= changeCarbonPriceFromYear] + yearColsFuture <- yearCols[as.integer(yearCols) >= fromYear] if (is.null(targetRegion)) { # GLOBAL: apply to all regions @@ -119,10 +121,11 @@ applyAlpha <- function(envWide, yearCols, alpha, targetRegion) { } writeFinalPolicyFiles <- function(envWide, yearCols, alphaFinal, region, + fromYear = changeCarbonPriceFromYear, canonicalPath = file.path("data","iEnvPolicies.csv"), updatedPath = file.path("data","iEnvPolicies_updated.csv"), alsoTimestamped = TRUE) { - envFinal <- applyAlpha(envWide, yearCols, alphaFinal, region) + envFinal <- applyAlpha(envWide, yearCols, alphaFinal, region, fromYear) dir.create(dirname(canonicalPath), showWarnings = FALSE, recursive = TRUE) dir.create(dirname(updatedPath), showWarnings = FALSE, recursive = TRUE) fwrite(envFinal, canonicalPath, na = "NA") @@ -162,6 +165,7 @@ run_gams <- function(gms = "main.gms", # Runs OPEN-PROM for a given alpha and returns the tracked emissions value. emissionsOPENPROM <- function(envWide, yearCols, alpha, targetRegion, targetYear, + fromYear = changeCarbonPriceFromYear, dataDir = "data", gms = "main.gms", gamsArgs = GAMSCmdArgs, @@ -172,7 +176,9 @@ emissionsOPENPROM <- function(envWide, yearCols, alpha, targetRegion, targetYear if (file.exists(backupCsvPath)) file.copy(backupCsvPath, canonicalCsv, overwrite = TRUE) }, add = TRUE) - fwrite(applyAlpha(envWide, yearCols, alpha, targetRegion), canonicalCsv, na = "NA") + # Remember the exact policy table sent to GAMS so it can be recovered if this run fails. + lastTestedPolicy <<- applyAlpha(envWide, yearCols, alpha, targetRegion, fromYear) + fwrite(lastTestedPolicy, canonicalCsv, na = "NA") ok <- run_gams(gms = gms, args = gamsArgs, log = log, echo_on_success = echo_on_success) if (!ok) stop("OPEN-PROM run failed for alpha=", alpha) @@ -226,35 +232,58 @@ alphaSeedLinear <- function(alpha0, E0, alphar, Er, Etarget, warn = TRUE, stopIf alpha0 + (Etarget - E0) * (alphar - alpha0) / (Er - E0) } -autoBracketFromSeed <- function(seedAlpha, budgetTarget, envWide, yearCols, targetRegion, targetYear, +autoBracketFromSeed <- function(seedAlpha, budgetTarget, envWide, yearCols, targetRegion, targetYear, + fromYear = changeCarbonPriceFromYear, minAlpha = 0.0, maxAlpha = 5.0, expandFactor = 1.35, maxProbes = 20, verbose = TRUE) { - if (verbose) message(sprintf("Seeding bracket near alpha ≈ %.4f", seedAlpha)) - probe <- function(a) emissionsOPENPROM(envWide, yearCols, a, targetRegion, targetYear) + probe <- function(a) emissionsOPENPROM(envWide, yearCols, a, targetRegion, targetYear, fromYear) + + # --- Feasibility probe: test the maximum carbon price first --- + # If even the highest allowed price cannot pull emissions down to the target, the + # target is unreachable — stop now (one run) instead of climbing toward it probe by probe. + if (verbose) message(sprintf("Feasibility probe at maxAlpha=%.4f", maxAlpha)) + Emax <- probe(maxAlpha) + if (verbose) message(sprintf("maxAlpha: alpha=%.4f -> E=%.6f (target=%.6f)", maxAlpha, Emax, budgetTarget)) + if (Emax > budgetTarget) { + stop(sprintf( + "maxAlpha too low: at alpha=%.4f emissions are %.6f, still above target %.6f. Increase maxAlpha.", + maxAlpha, Emax, budgetTarget)) + } + if (verbose) message(sprintf("maxAlpha feasible; seeding bracket near alpha %.4f", seedAlpha)) Eseed <- probe(seedAlpha) if (verbose) message(sprintf("Seed: alpha=%.4f → E=%.6f (target=%.6f)", seedAlpha, Eseed, budgetTarget)) if (Eseed <= budgetTarget) { - aU <- seedAlpha; EU <- Eseed - aL <- max(minAlpha, seedAlpha / expandFactor); tries <- 0 - repeat { - EL <- probe(aL); tries <- tries + 1 - if (verbose) message(sprintf("Down probe: alpha=%.4f → E=%.6f", aL, EL)) - if (EL > budgetTarget || tries >= maxProbes || aL <= minAlpha + 1e-9) break - aL <- max(minAlpha, aL / expandFactor) + # Seed already meets the target. Check whether unchanged prices (alpha = 0) meet it + # too — if so there is nothing to optimize, so keep the original prices and skip. + E0 <- probe(0) + if (verbose) message(sprintf("No-change probe: alpha=0.0000 -> E=%.6f (target=%.6f)", E0, budgetTarget)) + if (E0 <= budgetTarget) { + if (verbose) message("Region already meets target with unchanged carbon prices; skipping optimization.") + return(list(alreadyMet = TRUE, alpha = 0, + lowerAlpha = 0, upperAlpha = 0, EL = E0, EU = E0)) } - if (EL <= budgetTarget) stop("No failing lower bound; decrease minAlpha or revisit monotonicity.") + # alpha = 0 fails but the seed passes → [0, seed] brackets the target. + aL <- 0; EL <- E0 + aU <- seedAlpha; EU <- Eseed } else { + # Seed exceeds the target → expand upward toward the (already feasible) maxAlpha. aL <- seedAlpha; EL <- Eseed aU <- min(maxAlpha, seedAlpha * expandFactor); tries <- 0 repeat { + if (aU >= maxAlpha - 1e-9) { aU <- maxAlpha; EU <- Emax; break } # reuse feasibility probe EU <- probe(aU); tries <- tries + 1 if (verbose) message(sprintf("Up probe: alpha=%.4f → E=%.6f", aU, EU)) - if (EU <= budgetTarget || tries >= maxProbes || aU >= maxAlpha - 1e-9) break + if (EU <= budgetTarget || tries >= maxProbes) break + # Still above target: this probe is a tighter lower bound than the seed, so keep + # it instead of leaving aL stuck at seedAlpha (e.g. bracket [0.4, 1.6], not [0.1, 1.6]). + aL <- aU; EL <- EU aU <- min(maxAlpha, aU * expandFactor) } - if (EU > budgetTarget) stop("No passing upper bound; increase maxAlpha or revisit monotonicity.") + # Probes exhausted with the last one still above target: it is a valid lower bound too, + # so promote it before falling back to the known-feasible max. + if (EU > budgetTarget) { aL <- aU; EL <- EU; aU <- maxAlpha; EU <- Emax } } list(lowerAlpha = aL, upperAlpha = aU, EL = EL, EU = EU) } @@ -262,19 +291,20 @@ autoBracketFromSeed <- function(seedAlpha, budgetTarget, envWide, yearCols, targ findAlphaForBudget <- function(envWide, yearCols, budgetTarget, lowerAlpha, upperAlpha, eLow = NULL, eHigh = NULL, targetRegion, targetYear, + fromYear = changeCarbonPriceFromYear, tolAlphaRel = 1e-3, tolEmisAbs = 1e-3, maxIter = 60, verbose = TRUE, writeFinalCsv = TRUE) { # Evaluate bounds if not already provided by autoBracketFromSeed. - if (is.null(eLow)) eLow <- emissionsOPENPROM(envWide, yearCols, lowerAlpha, targetRegion, targetYear) - if (is.null(eHigh)) eHigh <- emissionsOPENPROM(envWide, yearCols, upperAlpha, targetRegion, targetYear) + if (is.null(eLow)) eLow <- emissionsOPENPROM(envWide, yearCols, lowerAlpha, targetRegion, targetYear, fromYear) + if (is.null(eHigh)) eHigh <- emissionsOPENPROM(envWide, yearCols, upperAlpha, targetRegion, targetYear, fromYear) if (verbose) message(sprintf( "Initial: aL=%.6f -> E=%.6f; aU=%.6f -> E=%.6f; target=%.6f", lowerAlpha, eLow, upperAlpha, eHigh, budgetTarget)) if (eLow <= budgetTarget) { - if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, lowerAlpha, targetRegion) + if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, lowerAlpha, targetRegion, fromYear) return(list(alpha = lowerAlpha, emissions = eLow, converged = TRUE, iters = 0)) } if (eHigh > budgetTarget) stop("upperAlpha still exceeds budget. Increase it or check monotonicity.") @@ -302,12 +332,12 @@ findAlphaForBudget <- function(envWide, yearCols, budgetTarget, } prevAM <- aM - emisM <- emissionsOPENPROM(envWide, yearCols, aM, targetRegion, targetYear) + emisM <- emissionsOPENPROM(envWide, yearCols, aM, targetRegion, targetYear, fromYear) if (verbose) message(sprintf("Iter %02d: aM=%.6f -> E=%.6f (target=%.6f)", it, aM, emisM, budgetTarget)) if (abs(emisM - budgetTarget) < tolEmisAbs || abs(aU - aL) / max(1.0, abs(aM)) < tolAlphaRel) { if (verbose) message("Converged.") - if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, aM, targetRegion) + if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, aM, targetRegion, fromYear) return(list(alpha = aM, emissions = emisM, converged = TRUE, iters = it)) } @@ -319,7 +349,7 @@ findAlphaForBudget <- function(envWide, yearCols, budgetTarget, } warning("Max iterations reached without strict tolerance convergence.") - if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, aU, targetRegion) + if (writeFinalCsv) writeFinalPolicyFiles(envWide, yearCols, aU, targetRegion, fromYear) list(alpha = aU, emissions = emisU, converged = FALSE, iters = it) } configureGamsFile <- function(gmsPath, targetRegion) { @@ -348,9 +378,23 @@ extractEmissions <- function(dataMagpie) { # Run # ---------------------------- start_time <- Sys.time() -GAMSCmdArgs <- c("--DevMode=0", "--GenerateInput=off", "lo=4", "idir=./data", "--CountrySolveMode=parallel") -selectedYear <- 2100 -changeCarbonPriceFromYear <- 2026 +selectedYear <- 2050 # default target year, overridable per region via targetList `year` +changeCarbonPriceFromYear <- 2026 # default first year the price is scaled, overridable per region via targetList `fromYear` + +# Model scenario passed to GAMS via --fScenario (overrides $evalGlobal fScenario in main.gms): +# 0 = No carbon price, 1 = NPi_Default, 2 = 1.5C, 3 = 2C +selectedScenario <- 2 + +# --fEndY caps the solve horizon at selectedYear instead of always running to 2100, so +# shortening a run is just a matter of lowering selectedYear (use --fEndY, never +# fEndHorizon, which triggers domain-violation errors). --fScenario selects the scenario. +GAMSCmdArgs <- c("--DevMode=0", "--GenerateInput=off", "lo=4", "idir=./data", + "--CountrySolveMode=parallel", + paste0("--fEndY=", selectedYear), + paste0("--fScenario=", selectedScenario)) + +# Keep an unmodified template of the GAMS args so we can substitute per-region end years +GAMSCmdArgsTemplate <- GAMSCmdArgs # --- Emissions variable to track --- # Set emissionsVariable to any variable name returned by reportEmissions(). @@ -360,8 +404,8 @@ changeCarbonPriceFromYear <- 2026 # "Emissions|CO2|Cumulated.Gt CO2" * 1000 -> Mt CO2 (cumulated) # "Emissions|CO2.Mt CO2/yr" * 1 -> Mt CO2/yr # "Emissions|Kyoto Gases.Mt CO2-equiv/yr" * 1 -> Mt CO2-equiv/yr -emissionsVariable <- "Emissions|CO2|Cumulated.Gt CO2" -emissionsScale <- 1000 # Gt -> Mt +emissionsVariable <- "Emissions|Kyoto Gases.Mt CO2-equiv/yr" +emissionsScale <- 1 # Gt -> Mt # EU27 member regions — share a single carbon price in iEnvPolicies.csv. # Never optimised individually; always solved as one aggregated group. @@ -370,7 +414,14 @@ EU27_REGIONS <- c("AUT","BEL","BGR","CYP","CZE","DEU","DNK","ESP","EST", "LVA","MLT","NLD","POL","PRT","ROU","SVK","SVN","SWE") # --- Target list --- -# Each entry is a named budget in the unit of emissionsVariable * emissionsScale. +# Each entry is one of: +# - a scalar numeric budget (backwards compatible): targetList$REGION = BUDGET +# - a named list: targetList$REGION = list(budget = BUDGET, year = YYYY, fromYear = YYYY) +# * budget = emissions budget the region must reach. +# * year = year by which the budget must be met (optional; defaults to selectedYear). +# * fromYear = first year whose carbon price is scaled (optional; defaults to +# changeCarbonPriceFromYear). Set it per region to change the price +# from a different start year than the global default. # Special keys: # "EU27" -> shared alpha applied to all 27 EU member rows; emissions summed over members. # "WORLD" -> alpha applied to all region rows; emissions summed globally. @@ -379,20 +430,21 @@ EU27_REGIONS <- c("AUT","BEL","BGR","CYP","CZE","DEU","DNK","ESP","EST", # # Current unit: cumulated Mt CO2 (Emissions|CO2|Cumulated.Gt CO2 * 1000) targetList <- list( - #"WORLD" = 1257571, # optional: comment out to skip world run - "EU27" = 92917, # sum of all 27 EU member budgets - "CAZ" = 22517, - "CHA" = 307728, - "GBR" = 14692, - "IND" = 223711, - "JPN" = 29285, - "LAM" = 117008, - "MEA" = 102029, - "NEU" = 22276, - "OAS" = 220646, - "REF" = 64216, - "SSA" = 175989, - "USA" = 100855 + # Examples (mix-and-match supported): + # WORLD = 1257571, # optional: comment out to skip world run + EU27 = list(budget = 0, year = 2050), + CAZ = list(budget = 0, year = 2050), + # CHA = list(budget = 13447, year = 2050), + GBR = list(budget = 0, year = 2050), + IND = list(budget = 0, year = 2070), + JPN = list(budget = 0, year = 2050), + LAM = list(budget = 703, year = 2050), + MEA = list(budget = 3245, year = 2060), + NEU = list(budget = 101, year = 2050), + OAS = list(budget = 1186, year = 2060), + REF = list(budget = 355, year = 2060), + SSA = list(budget = 2238, year = 2050), + USA = list(budget = 0, year = 2050) # numeric form still supported: interpreted as budget with fallback year = selectedYear ) logFilePath <- "Carbon_price_optimization.log" @@ -421,7 +473,20 @@ resultsLog <- list() for (regName in names(targetList)) { - bg <- targetList[[regName]] + lastTestedPolicy <- NULL # don't carry a previous region's tested policy into this run + entry <- targetList[[regName]] + # Support two formats: numeric (budget only) or list(budget=..., year=..., fromYear=...) + if (is.list(entry) && !is.null(entry$budget)) { + bg <- as.numeric(entry$budget) + regionTargetYear <- if (!is.null(entry$year)) as.integer(entry$year) else selectedYear + regionFromYear <- if (!is.null(entry$fromYear)) as.integer(entry$fromYear) else changeCarbonPriceFromYear + } else if (is.numeric(entry) && length(entry) == 1) { + bg <- as.numeric(entry) + regionTargetYear <- selectedYear + regionFromYear <- changeCarbonPriceFromYear + } else { + stop(sprintf("Invalid targetList entry for '%s' — must be numeric or list(budget=..., year=..., fromYear=...)", regName)) + } if (regName == "WORLD") { actualRegion <- NULL # NULL -> alpha applied to all region rows @@ -434,7 +499,7 @@ for (regName in names(targetList)) { displayName <- regName } - message(sprintf("\n--- Optimizing %s (Target: %.4f) ---", displayName, bg)) + message(sprintf("\n--- Optimizing %s (Target: %.4f, Year: %d, From: %d) ---", displayName, bg, regionTargetYear, regionFromYear)) skipRegion <- FALSE tryCatch({ @@ -444,49 +509,68 @@ for (regName in names(targetList)) { configureGamsFile("main.gms", actualRegion) } + # Ensure GAMS runs use the region-specific solve horizon (end year) + GAMSCmdArgs <- GAMSCmdArgsTemplate + i_endy <- grep("^--fEndY=", GAMSCmdArgs) + if (length(i_endy)) GAMSCmdArgs[i_endy] <- paste0("--fEndY=", regionTargetYear) else GAMSCmdArgs <- c(GAMSCmdArgs, paste0("--fEndY=", regionTargetYear)) + brkt <- autoBracketFromSeed( seedAlpha = 0.1, budgetTarget = bg, envWide = currentEnvWide, yearCols = yearCols, targetRegion = actualRegion, - targetYear = selectedYear, + targetYear = regionTargetYear, + fromYear = regionFromYear, minAlpha = -0.5, # Allow price reduction up to -50% if needed - maxAlpha = 10.0, # Allow up to +1000% increase - expandFactor = 3.0, + maxAlpha = 40, # Allow up to +1000% increase + expandFactor = 4.0, maxProbes = 7, verbose = TRUE ) - # C. Solve - solveResult <- findAlphaForBudget( - envWide = currentEnvWide, - yearCols = yearCols, - budgetTarget = bg, - targetRegion = actualRegion, # Passes NULL if global - targetYear = selectedYear, - lowerAlpha = brkt$lowerAlpha, - upperAlpha = brkt$upperAlpha, - eLow = brkt$EL, - eHigh = brkt$EU, - tolAlphaRel = 1e-2, - tolEmisAbs = 1e+1, - maxIter = 60, - verbose = TRUE, - writeFinalCsv = FALSE - ) - - finalAlpha <- solveResult$alpha - message(sprintf(" -> Converged %s: Alpha=%.3f", displayName, finalAlpha)) + # C. Solve (skip the root-find when the region already meets its target unchanged) + if (isTRUE(brkt$alreadyMet)) { + finalAlpha <- 0 + message(sprintf(" -> %s already meets target; carbon prices left unchanged (Alpha=0).", displayName)) + } else { + solveResult <- findAlphaForBudget( + envWide = currentEnvWide, + yearCols = yearCols, + budgetTarget = bg, + targetRegion = actualRegion, # Passes NULL if global + targetYear = regionTargetYear, + fromYear = regionFromYear, + lowerAlpha = brkt$lowerAlpha, + upperAlpha = brkt$upperAlpha, + eLow = brkt$EL, + eHigh = brkt$EU, + tolAlphaRel = 1e-2, + tolEmisAbs = 1e+1, + maxIter = 60, + verbose = TRUE, + writeFinalCsv = FALSE + ) + + finalAlpha <- solveResult$alpha + message(sprintf(" -> Converged %s: Alpha=%.3f", displayName, finalAlpha)) + } # Apply converged alpha and persist as the new baseline for subsequent regions - currentEnvWide <- applyAlpha(currentEnvWide, yearCols, finalAlpha, actualRegion) + currentEnvWide <- applyAlpha(currentEnvWide, yearCols, finalAlpha, actualRegion, regionFromYear) fwrite(currentEnvWide, inputCsvPath, na = "NA") file.copy(inputCsvPath, backupCsvPath, overwrite = TRUE) resultsLog[[regName]] <- list(status = "OK", alpha = finalAlpha) }, error = function(e) { message(sprintf(" !! FAILURE for %s: %s", displayName, e$message)) + # Preserve the carbon prices that were being tested when the run failed, in a + # region-specific file, so each failing region's last-tested policy survives. + if (!is.null(lastTestedPolicy)) { + regionLastTestedCsvPath <- sub("\\.csv$", paste0("_", regName, ".csv"), lastTestedCsvPath) + fwrite(lastTestedPolicy, regionLastTestedCsvPath, na = "NA") + message(sprintf(" -> Saved last-tested policy to %s", regionLastTestedCsvPath)) + } message(" -> Reverting to last good state and skipping.") if (file.exists(backupCsvPath)) file.copy(backupCsvPath, inputCsvPath, overwrite = TRUE) resultsLog[[regName]] <<- list(status = "FAILED", error = e$message)