diff --git a/core/bounds.gms b/core/bounds.gms index 5a618be807..91f3c51168 100755 --- a/core/bounds.gms +++ b/core/bounds.gms @@ -326,8 +326,8 @@ if(c_ccsinjecratescen > 0, *** Potential of EU27 regions is pooled and redistributed according to GDP (Only upper limit for 2030) *** Norway and UK announced to store CO2 for EU27 countries. So 50% of Norway and UK potential in 2030 is attributed to EU27-Pool if(not cm_emiscen = 1, !! cm_emiscen 1 = BAU - vm_co2CCS.lo(t,regi,"cco2","ico2","ccsinje","1") $ (t.val <= 2030) = s_MtCO2_2_GtC * p_boundCapCCS(t,regi,"operational") $ (t.val <= 2030); - vm_co2CCS.up(t,regi,"cco2","ico2","ccsinje","1") $ (t.val <= 2030) = s_MtCO2_2_GtC * ( + vm_co2CCS.lo(t,regi,"cco2","ico2","ccsinje","1") $ (t.val <= 2030) = sm_MtCO2_2_GtC * p_boundCapCCS(t,regi,"operational") $ (t.val <= 2030); + vm_co2CCS.up(t,regi,"cco2","ico2","ccsinje","1") $ (t.val <= 2030) = sm_MtCO2_2_GtC * ( p_boundCapCCS(t,regi,"operational") $ (t.val <= 2030) + p_boundCapCCS(t,regi,"construction") $ (t.val <= 2030) + p_boundCapCCS(t,regi,"planned") $ (t.val <= 2030) * c_fracRealfromAnnouncedCCScap2030); @@ -372,7 +372,7 @@ loop(regi $ (p_boundCapCCSindicator(regi) = 0), *** Limit REMINDs ability to vent captured CO2 to 1 MtCO2 per yr per region. This happens otherwise to a great extend in stringent climate *** policy scenarios if CCS and CCU capacities are limited in early years, to lower overall adjustment costs of capture technologies. *** In reality, people don't have perfect foresight and without storage or usage capacities, no capture facilities will be built. -v_co2capturevalve.up(t,regi) = 1 * s_MtCO2_2_GtC; +v_co2capturevalve.up(t,regi) = 1 * sm_MtCO2_2_GtC; *** ================================================================== diff --git a/core/datainput.gms b/core/datainput.gms index cb059124b6..756524df48 100644 --- a/core/datainput.gms +++ b/core/datainput.gms @@ -561,15 +561,14 @@ $offdelim *** carbon intensities of coal, oil, and gas *** emissions factor of primary energy fossil fuels -pm_cintraw("pecoal") = 26.1 / s_zj_2_twa; -pm_cintraw("peoil") = 20.0 / s_zj_2_twa; -pm_cintraw("pegas") = 15.0 / s_zj_2_twa; +pm_cintraw("pecoal") = 26.1 / s_ZJ_2_TWa; +pm_cintraw("peoil") = 20.0 / s_ZJ_2_TWa; +pm_cintraw("pegas") = 15.0 / s_ZJ_2_TWa; $ifthen.tech_CO2capturerate not "%c_tech_CO2capturerate%" == "off" p_PECarriers_CarbonContent(peFos)=pm_cintraw(peFos); *** From conversation: 25 GtC/ZJ is the assumed carbon content of PE biomass (makes default bioh2c capture rate 90%) -*** Convert to GtC/TWa -p_PECarriers_CarbonContent("pebiolc")=25 / s_zj_2_twa; +p_PECarriers_CarbonContent("pebiolc") = 25 / s_ZJ_2_TWa; !! convert 25 GtC/ZJ to GtC/TWa loop(pe2se(entyPe,entySe,te)$(p_tech_co2capturerate(te)), if(p_tech_co2capturerate(te) gt 0, if(p_tech_co2capturerate(te) ge 1, @@ -577,7 +576,7 @@ loop(pe2se(entyPe,entySe,te)$(p_tech_co2capturerate(te)), ); *** Alter CO2 capture rate in f_dataemiglob *** f_dataemiglob is given in GtC/ZJ - f_dataemiglob(entyPe,entySe,te,"cco2") = p_tech_co2capturerate(te) * p_PECarriers_CarbonContent(entyPe) * s_zj_2_twa; + f_dataemiglob(entyPe,entySe,te,"cco2") = p_tech_co2capturerate(te) * p_PECarriers_CarbonContent(entyPe) * s_ZJ_2_TWa; if(sameAs(entyPe,"pebiolc"), f_dataemiglob(entyPe,entySe,te,"co2") = -f_dataemiglob(entyPe,entySe,te,"cco2") ; else @@ -602,8 +601,8 @@ $ifthen "%c_SSP_forcing_adjust%" == "forcing_SSP5" $endif ); *nb* specific emissions of transformation technologies (co2 in gtc/zj -> conv. gtc/twyr): -f_dataemiglob(enty,enty2,te,"co2")$pe2se(enty,enty2,te) = 1/s_zj_2_twa * f_dataemiglob(enty,enty2,te,"co2"); -f_dataemiglob(enty,enty2,te,"cco2") = 1/s_zj_2_twa * f_dataemiglob(enty,enty2,te,"cco2"); +f_dataemiglob(enty,enty2,te,"co2")$pe2se(enty,enty2,te) = 1/s_ZJ_2_TWa * f_dataemiglob(enty,enty2,te,"co2"); +f_dataemiglob(enty,enty2,te,"cco2") = 1/s_ZJ_2_TWa * f_dataemiglob(enty,enty2,te,"cco2"); table f_dataetaglob(tall,all_te) "global eta data" $include "./core/input/generisdata_varying_eta.prn" @@ -1545,9 +1544,9 @@ $include "./core/input/f_nechem_emissionFactors.cs4r" $offdelim /; -pm_emifacNonEnergy(ttot,regi,"sesofos", "fesos","indst","co2") = f_nechem_emissionFactors(ttot,regi,"solids") / s_zj_2_twa; -pm_emifacNonEnergy(ttot,regi,"seliqfos","fehos","indst","co2") = f_nechem_emissionFactors(ttot,regi,"liquids") / s_zj_2_twa; -pm_emifacNonEnergy(ttot,regi,"segafos", "fegas","indst","co2") = f_nechem_emissionFactors(ttot,regi,"gases") / s_zj_2_twa; +pm_emifacNonEnergy(ttot,regi,"sesofos", "fesos","indst","co2") = f_nechem_emissionFactors(ttot,regi,"solids") / s_ZJ_2_TWa; +pm_emifacNonEnergy(ttot,regi,"seliqfos","fehos","indst","co2") = f_nechem_emissionFactors(ttot,regi,"liquids") / s_ZJ_2_TWa; +pm_emifacNonEnergy(ttot,regi,"segafos", "fegas","indst","co2") = f_nechem_emissionFactors(ttot,regi,"gases") / s_ZJ_2_TWa; ***------ Read in projections for incineration rates of plastic waste--- *** "incineration rates [fraction]" diff --git a/core/declarations.gms b/core/declarations.gms index 36ced4acbc..22cb65f9e8 100644 --- a/core/declarations.gms +++ b/core/declarations.gms @@ -613,7 +613,7 @@ sm_trillion_2_non "trillion to non" /1e+12/, *** energy units pm_conv_TWa_EJ "conversion from TWa to EJ" /31.536/, -s_zj_2_twa "zeta joule to tw year" /31.7098/, +s_ZJ_2_TWa "zeta joule to tw year" /31.71/, sm_EJ_2_TWa "multiplicative factor to convert from EJ to TWa" /31.71e-03/, sm_GJ_2_TWa "multiplicative factor to convert from GJ to TWa" /31.71e-12/, sm_TWa_2_TWh "tera Watt year to Tera Watt hour" /8.76e+3/, @@ -628,7 +628,7 @@ sm_c_2_co2 "conversion from c to co2" /3.666666 s_NO2_2_N "convert NO2 to N [14 / (14 + 2 * 16)]" / .304 / sm_tgn_2_pgc "conversion factor 100-yr GWP from TgN to PgCeq" sm_tgch4_2_pgc "conversion factor 100-yr GWP from TgCH4 to PgCeq" -s_MtCO2_2_GtC "conversion factor from MtCO2 to native REMIND emission unit GtC" /2.727e-04/ +sm_MtCO2_2_GtC "conversion factor from MtCO2 to native REMIND emission unit GtC" /2.727e-04/ s_MtCH4_2_TWa "Energy content of methane. MtCH4 --> TWa: 1 MtCH4 = 1.23 * 10^6 toe * 42 GJ/toe * 10^-9 EJ/GJ * 1 TWa/31.536 EJ = 0.001638 TWa (BP statistical review)" /0.001638/ s_gwpCH4 "Global Warming Potentials of CH4, AR5 WG1 CH08 Table 8.7" /28/ s_gwpN2O "Global Warming Potentials of N2O, AR5 WG1 CH08 Table 8.7" /265/ diff --git a/core/loop.gms b/core/loop.gms index c9edb955e3..a477513b8a 100644 --- a/core/loop.gms +++ b/core/loop.gms @@ -6,9 +6,9 @@ *** | Contact: remind@pik-potsdam.de *** SOF ./core/loop.gms -*-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** solveoptions -*-------------------------------------------------------------------------- +***------------------------------------------------------------------- option limcol = 0; option limrow = 0; hybrid.optfile = 1; @@ -21,158 +21,130 @@ o_modelstat = 100; $ifthen.calibrate "%CES_parameters%" == "calibrate" !! CES_parameters $ifthen.subsectors "%industry%" == "subsectors" !! industry -!! Calibrating industry/subsectors lead to random infeasibilities on the order -!! of 1e-15. Relaxing this attribute a little solves this problem. +*** Calibrating industry/subsectors lead to random infeasibilities on the order of 1e-15. +*** Relaxing this attribute a little solves this problem. hybrid.tolinfeas = 1e-14; $endif.subsectors $endif.calibrate ***------------------------------------------------------------------- -*** read GDX +*** read GDX ***------------------------------------------------------------------- -*** load start gdx - execute_loadpoint "input"; -***-------------------------------------------------------------------------- -*** start iteration loop -***-------------------------------------------------------------------------- - -LOOP(iteration $(ord(iteration)<(cm_iteration_max+1)), +***------------------------------------------------------------------- +*** start iteration loop +***------------------------------------------------------------------- +loop(iteration $ (iteration.val <= cm_iteration_max), + if(iteration.val = cm_iteration_max, + option solprint = on + ); - IF(ord(iteration)>(cm_iteration_max-1), - OPTION solprint=on - ); -*-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** BOUNDS -*-------------------------------------------------------------------------- +***------------------------------------------------------------------- $include "./core/bounds.gms"; $batinclude "./modules/include.gms" bounds -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** PRESOLVE -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- $include "./core/presolve.gms"; $batinclude "./modules/include.gms" presolve -*** Fixing information (.L, .FX and .M) from run to be fixed to is read in from input_ref.gdx (t < cm_startyear) -*** happens via submit.R script (files levs.gms, fixings.gms, margs.gms) -*** submit.R looks for the unique string in the following line and replaces it with the offlisting include into the full.gms at this position -***cb20140305readinpositionforfixingfiles +*** When there is a reference run, input_ref.gdx contains the necessary fixing information (.L, .FX and .M) for t < cm_startyear +*** Script submit.R creates reference files (levs.gms, fixings.gms, margs.gms) and inludes them in full.gms by replacing the following line +*** cb20140305readinpositionforfixingfiles -*** In case of fixing, fix to prices from input_ref.gdx (t < cm_startyear). -*** Parameters are not automatically treated by the fixing mechanism above. -if( (cm_startyear gt 2005), +*** Also fix prices, which are not automatically treated by the fixing mechanism above. + if(cm_startyear > 2005, Execute_Loadpoint "input_ref" p_pvpRef = pm_pvp; - pm_pvp(ttot,trade)$( (ttot.val ge 2005) and (ttot.val lt cm_startyear) and (NOT tradeSe(trade))) = p_pvpRef(ttot,trade); -); + pm_pvp(ttot,trade) $ (ttot.val >= 2005 and ttot.val < cm_startyear and not tradeSe(trade)) = p_pvpRef(ttot,trade); + ); -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** SOLVE -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** Set options for debugging -if (cm_nash_mode eq 1, - option - solprint = on - limcol = 2147483647 - limrow = 2147483647 - ; -); - + if(cm_nash_mode = 1, + option + solprint = on + limcol = 2147483647 + limrow = 2147483647 + ; + ); o_modelstat = 100; -loop(sol_itr$(sol_itr.val <= cm_solver_try_max), - if(o_modelstat ne 2, +loop(sol_itr $ (sol_itr.val <= cm_solver_try_max), + if(o_modelstat ne 2, $batinclude "./modules/include.gms" solve - ) + ) ); !! end of sol_itr loop, when o_modelstat is not equal to 2 -***--------------------------------------------------------- -*** Track of changes between iterations -***--------------------------------------------------------- -loop(entyPe$(NOT sameas(entyPe,"peur")), - o_negitr_cumulative_peprod(iteration,entyPe) = 0.031536 - * sum(regi, - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), vm_prodPe.l(ttot,regi,entyPe) * pm_ts(ttot) ) - + sum(ttot$(ttot.val eq 2005), vm_prodPe.l(ttot,regi,entyPe) * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), vm_prodPe.l(ttot,regi,entyPe) * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) - ); -); -o_negitr_cumulative_peprod(iteration,"peur") = -sum(regi, - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), sum(pe2rlf("peur",rlf), 0.4102 * vm_prodPe.l(ttot,regi,"peur") * pm_ts(ttot) ) ) - + sum(ttot$(ttot.val eq 2005), 0.4102 * vm_prodPe.l(ttot,regi,"peur") * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), 0.4102 * vm_prodPe.l(ttot,regi,"peur") * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) -); -o_negitr_cumulative_CO2_emineg_co2luc(iteration) = -sum(regi, - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2luc") * pm_ts(ttot) ) - + sum(ttot$(ttot.val eq 2005), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2luc") * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2luc") * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) -); +***------------------------------------------------------------------- +*** Track of changes between iterations +***------------------------------------------------------------------- +o_negitr_cumulative_peprod(iteration,entyPe) = + ((1 / s_ZJ_2_TWa) $ (not sameas(entyPe,"peur")) + 0.4102 $ (sameas(entyPe,"peur"))) !! conversion from TWa (or Mt uranium) to ZJ + * sum(regi, + sum(ttot $ (ttot.val < 2100 and ttot.val > 2005), vm_prodPe.l(ttot,regi,entyPe) * pm_ts(ttot) ) + + sum(ttot $ (ttot.val = 2005), vm_prodPe.l(ttot,regi,entyPe) * pm_ts(ttot) * 0.5 ) + + sum(ttot $ (ttot.val = 2100), vm_prodPe.l(ttot,regi,entyPe) * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5 ) + ); -o_negitr_cumulative_CO2_emineg_cement(iteration) = -sum(regi, - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2cement_process") * pm_ts(ttot) ) - + sum(ttot$(ttot.val eq 2005), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2cement_process") * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), 3.6667 * vm_emiMacSector.l(ttot,regi,"co2cement_process") * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) -); -o_negitr_cumulative_CO2_emieng_seq(iteration) - = - 3.6667 +o_negitr_cumulative_CO2_emineg_co2luc(iteration) = sm_c_2_co2 !! conversion from carbon to CO2 + * sum(regi, + sum(ttot $ (ttot.val < 2100 and ttot.val > 2005), vm_emiMacSector.l(ttot,regi,"co2luc") * pm_ts(ttot) ) + + sum(ttot $ (ttot.val = 2005), vm_emiMacSector.l(ttot,regi,"co2luc") * pm_ts(ttot) * 0.5 ) + + sum(ttot $ (ttot.val = 2100), vm_emiMacSector.l(ttot,regi,"co2luc") * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5 ) + ); + +o_negitr_cumulative_CO2_emineg_cement(iteration) = sm_c_2_co2 !! conversion from carbon to CO2 * sum(regi, - sum((ttot,emi2te(enty,enty2,te,"cco2"))$( ttot.val gt 2005 AND ttot.val lt 2100 ), - vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") - * pm_ts(ttot) - ) - + sum((ttot,emi2te(enty,enty2,te,"cco2"))$( ttot.val eq 2005 ), - vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") - * pm_ts(ttot) - / 2 - ) - + sum((ttot,emi2te(enty,enty2,te,"cco2"))$( ttot.val eq 2100 ), - vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") - * (pm_ttot_val(ttot) - pm_ttot_val(ttot-1)) - / 2 - ) - ) -; + sum(ttot $ (ttot.val < 2100 and ttot.val > 2005), vm_emiMacSector.l(ttot,regi,"co2cement_process") * pm_ts(ttot) ) + + sum(ttot $ (ttot.val = 2005), vm_emiMacSector.l(ttot,regi,"co2cement_process") * pm_ts(ttot) * 0.5 ) + + sum(ttot $ (ttot.val = 2100), vm_emiMacSector.l(ttot,regi,"co2cement_process") * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5 ) + ); + +o_negitr_cumulative_CO2_emieng_seq(iteration) = sm_c_2_co2 !! conversion from carbon to CO2 + * sum((regi,emi2te(enty,enty2,te,"cco2")), + sum(ttot $ ( ttot.val > 2005 and ttot.val < 2100 ), vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") * pm_ts(ttot)) + + sum(ttot $ ( ttot.val = 2005 ), vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") * pm_ts(ttot) * 0.5) + + sum(ttot $ ( ttot.val = 2100 ), vm_emiTeDetail.l(ttot,regi,enty,enty2,te,"cco2") * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5) + ); + o_negitr_disc_cons_dr5_reg(iteration,regi) = - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * pm_ts(ttot) ) - + sum(ttot$(ttot.val eq 2005), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) -; + sum(ttot $ ( (ttot.val < 2100) and (ttot.val > 2005)), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * pm_ts(ttot) ) + + sum(ttot $ (ttot.val = 2005), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * pm_ts(ttot) * 0.5 ) + + sum(ttot $ (ttot.val = 2100), vm_cons.l(ttot,regi) * (0.95 ** (pm_ttot_val(ttot) - s_t_start)) * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5 ); + o_negitr_disc_cons_drInt_reg(iteration,regi) = - sum(ttot$( (ttot.val lt 2100) AND (ttot.val gt 2005)), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi)/ (qm_budget.m("2005",regi) + 1.e-8) * pm_ts(ttot) ) - + sum(ttot$(ttot.val eq 2005), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi)/ (qm_budget.m("2005",regi) + 1.e-8) * pm_ts(ttot) * 0.5 ) - + sum(ttot$(ttot.val eq 2100), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi)/ (qm_budget.m("2005",regi) + 1.e-8) * ( pm_ttot_val(ttot)- pm_ttot_val(ttot-1) ) * 0.5 ) -; + sum(ttot $ ( (ttot.val < 2100) and (ttot.val > 2005)), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi)/ (qm_budget.m("2005",regi) + sm_eps) * pm_ts(ttot) ) + + sum(ttot $ (ttot.val = 2005), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi) / (qm_budget.m("2005",regi) + sm_eps) * pm_ts(ttot) * 0.5 ) + + sum(ttot $ (ttot.val = 2100), vm_cons.l(ttot,regi) * qm_budget.m(ttot,regi) / (qm_budget.m("2005",regi) + sm_eps) * ( pm_ttot_val(ttot) - pm_ttot_val(ttot-1) ) * 0.5 ); -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- *** POSTSOLVE -***-------------------------------------------------------------------------- +***------------------------------------------------------------------- $include "./core/postsolve.gms"; $batinclude "./modules/include.gms" postsolve -*-------------------------------------------------------------------------- -*** save gdx -*-------------------------------------------------------------------------- -*** write the fulldata.gdx file after each optimal iteration -*** In Nash status 7 is considered optimal in that respect (see definition of -*** o_modelstat in solve.gms) +***------------------------------------------------------------------- +*** save gdx +***------------------------------------------------------------------- +*** Write the fulldata.gdx file after each optimal iteration +*** In Nash, status 7 is considered optimal (see definition of o_modelstat in solve.gms) logfile.nr = 1; -if (o_modelstat le 2, +if(o_modelstat <= 2, execute_unload "fulldata"; - !! retain gdxes of intermediate iterations by copying them using shell - !! commands - if (c_keep_iteration_gdxes eq 1, + if(c_keep_iteration_gdxes = 1, !! save gdx of intermediate iterations using shell command "copy" put_utility logfile, "shell" / "cp fulldata.gdx fulldata_" iteration.val:0:0 ".gdx"; ); else execute_unload "non_optimal"; - if (c_keep_iteration_gdxes eq 1, + if(c_keep_iteration_gdxes = 1, put_utility logfile, "shell" / "cp non_optimal.gdx non_optimal_" iteration.val:0:0 ".gdx"; ); diff --git a/main.gms b/main.gms index 384ff266e5..807659d214 100755 --- a/main.gms +++ b/main.gms @@ -504,9 +504,9 @@ parameter ; cm_nash_autoconverge = 1; !! def = 1 !! regexp = [0-3] *' * (0): manually set number of iterations by adjusting cm_iteration_max -*' * (1): run until solution is sufficiently converged - coarse tolerances, quick solution. ! do not use in production runs ! -*' * (2): run until solution is sufficiently converged - fine tolerances, for production runs. -*' * (3): run until solution is sufficiently converged using very relaxed targets - very coarse tolerances, two times higher than option 1. ! do not use in production runs ! +*' * (1): run until solution is sufficiently converged: default tolerance, for quick solution, not recommended in production runs. +*' * (2): run until solution is sufficiently converged: fine tolerance (0.2x option 1), for production runs. +*' * (3): run until solution is sufficiently converged: coarse tolerance (2x option 1), do not use in production runs! *' parameter cm_emiscen "policy scenario choice" diff --git a/modules/04_PE_FE_parameters/iea2014/declarations.gms b/modules/04_PE_FE_parameters/iea2014/declarations.gms index fc903823eb..0397bc9571 100644 --- a/modules/04_PE_FE_parameters/iea2014/declarations.gms +++ b/modules/04_PE_FE_parameters/iea2014/declarations.gms @@ -11,7 +11,7 @@ pm_IO_input(all_regi,all_enty,all_enty,all_te) "Energy input ba pm_IO_output(tall,all_regi,all_enty,all_enty,all_te) "Historical energy output per technology based on IEA data [TWa]" p04_IO_output(all_regi,all_enty,all_enty,all_te) "Energy output based on IEA data" p04_x_enty2te_dyn04(all_regi,all_enty,all_enty,all_te,all_te) "parameter for the allocation of energy flow to technologies" -pm_prodCouple(all_regi,all_enty,all_enty,all_te,all_enty) "own consumption" +pm_prodCouple(all_regi,all_enty,all_enty,all_te,all_enty) "share of couple production. Value 0.2 means 20EJ of couple product for 100EJ of main product; value -0.03 means 3EJ of own consumption for 100EJ of main product." p04_aux_data(all_regi,char, all_te) "auxiliary parameter to store the initial mix0 and eta values for gas electricity before splitting it to ngcc and ngt (needed as long as calibration routine sets ngt to 0)" p04_shareNGTinGas(all_regi) "Share of ngt in electricity produced from gas" pm_fuExtrOwnCons(all_regi, all_enty, all_enty) "energy own consumption in the extraction sector with first enty being the output produced and the second enty being the input required" diff --git a/modules/80_optimization/module.gms b/modules/80_optimization/module.gms index 41b38c0f34..1017982a41 100644 --- a/modules/80_optimization/module.gms +++ b/modules/80_optimization/module.gms @@ -11,7 +11,7 @@ *' @description The optimization module gives the opportunity to choose different solution algorithms. *' *' -*' @authors Anastasis Giannousakis, Marian Leimbach, Lavinia Baumstark, Gunnar Luderer, Anselm Schultes +*' @authors Anastasis Giannousakis, Marian Leimbach, Lavinia Baumstark, Gunnar Luderer, Anselm Schultes, Fabrice Lécuyer *###################### R SECTION START (MODULETYPES) ########################## $Ifi "%optimization%" == "nash" $include "./modules/80_optimization/nash/realization.gms" diff --git a/modules/80_optimization/nash/bounds.gms b/modules/80_optimization/nash/bounds.gms index d1d069a654..1211a26140 100644 --- a/modules/80_optimization/nash/bounds.gms +++ b/modules/80_optimization/nash/bounds.gms @@ -5,19 +5,18 @@ *** | REMIND License Exception, version 1.0 (see LICENSE file). *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/bounds.gms -*override fixing of vm_perm set in 41_emicapregi/none/bounds.gms -*mlb 20150609* allow negative permits to compute efficient solution + +if(cm_emiscen = 6, !! budget run $ifthen.emiopt %emicapregi% == 'none' -if(cm_emiscen eq 6, - vm_perm.lo(t,regi) = -10; - vm_perm.up(t,regi) = 10000; -); + vm_perm.up(t,regi) = 10000; !! override fixing of vm_perm set in 41_emicapregi/none/bounds.gms + vm_perm.lo(t,regi) = -10; !! mlb 20150609: allow negative permits to compute efficient solution $endif.emiopt -*ML* in nash with permit allocation only total budgets are meaningful; allowing permit trade -*only for the initial policy period avoids indeterminacy, hence numerical problems -loop(ttot$(ttot.val ne cm_startyear), - vm_Xport.fx(ttot,regi,"perm")$(cm_emiscen eq 6) = 0; - vm_Mport.fx(ttot,regi,"perm")$(cm_emiscen eq 6) = 0 ; +*** ML: in nash with permit allocation only total budgets are meaningful; +*** allowing permit trade only for the initial policy period avoids indeterminacy, hence numerical problems + loop(ttot $ (ttot.val ne cm_startyear), + vm_Xport.fx(ttot,regi,"perm") = 0; + vm_Mport.fx(ttot,regi,"perm") = 0; + ); ); *** EOF ./modules/80_optimization/nash/bounds.gms diff --git a/modules/80_optimization/nash/datainput.gms b/modules/80_optimization/nash/datainput.gms index 2588cb4055..9cea72df23 100644 --- a/modules/80_optimization/nash/datainput.gms +++ b/modules/80_optimization/nash/datainput.gms @@ -5,76 +5,75 @@ *** | REMIND License Exception, version 1.0 (see LICENSE file). *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/datainput.gms -pm_w(regi) = 1; -o80_trackSurplusSign(ttot,trade,iteration)$(NOT tradeSe(trade)) = 0; -*MLB 20130920* initialization only -pm_cumEff(t, regi, in) = 100; +*** convergence with trade surplus thresholds +if(cm_nash_autoconverge > 0, + cm_iteration_max = 100; !! set max number of iterations -*MLB 20140109* initialization of climate externality is sensitive -pm_co2eqForeign(t, regi) = (1 - pm_shPerm(t,regi)) * pm_emicapglob(t); - -***convergence mode -if(cm_nash_autoconverge gt 0, - -*** set max number of iterations -cm_iteration_max = 100; - - if(cm_nash_autoconverge eq 1, -***convergences thresholds - coarse - p80_surplusMaxTolerance(tradePe) = 1.5 * sm_EJ_2_TWa; !! convert EJ/yr into internal unit TWa - p80_surplusMaxTolerance("good") = 100/1000; !! in internal unit, trillion Dollar - p80_surplusMaxTolerance("perm") = 300 * 12/44 / 1000; !! convert MtCO2eq into internal unit GtC - ); - if(cm_nash_autoconverge eq 2, -***convergences thresholds - fine - p80_surplusMaxTolerance(tradePe) = 0.3 * sm_EJ_2_TWa; !! convert EJ/yr into internal unit TWa - p80_surplusMaxTolerance("good") = 20/1000; !! in internal unit, trillion Dollar - p80_surplusMaxTolerance("perm") = 70 * 12/44 / 1000 ; !! convert MtCO2eq into internal unit GtC - ); - if(cm_nash_autoconverge eq 3, -***convergences thresholds - very coarse - p80_surplusMaxTolerance(tradePe) = 2* 1.5 * sm_EJ_2_TWa; !! convert EJ/yr into internal unit TWa - p80_surplusMaxTolerance("good") = 2* 100/1000; !! in internal unit, trillion Dollar - p80_surplusMaxTolerance("perm") = 2* 300 * 12/44 / 1000; !! convert MtCO2eq into internal unit GtC - ); +*** default values (cm_nash_autoconverge = 1) + p80_surplusMaxTolerance(tradePe) = 1.5 * sm_EJ_2_TWa; !! 1.5 EJ/yr, converted into internal unit TWa + p80_surplusMaxTolerance("good") = 0.1; !! 0.1 trillion Dollar + p80_surplusMaxTolerance("perm") = 300 * sm_MtCO2_2_GtC; !! 300 MtCO2eq, converted into internal unit GtC + + loop(trade $ (tradePe(trade) or tradeMacro(trade)), +*** fine convergence thresholds (cm_nash_autoconverge = 2) + p80_surplusMaxTolerance(trade) $ (cm_nash_autoconverge = 2) = 0.2 * p80_surplusMaxTolerance(trade); +*** coarse convergence thresholds (cm_nash_autoconverge = 3) + p80_surplusMaxTolerance(trade) $ (cm_nash_autoconverge = 3) = 2 * p80_surplusMaxTolerance(trade); + ); ); - -*Nash adjustment costs. Involves a trade-off: If set too low, markets jump far away from clearance. Set too high, changes in trade patten over iterations are very slow, convergence takes many many iterations. Default value around 150 -p80_etaAdj(tradePe) = 80; -p80_etaAdj("good") = 100; -p80_etaAdj("perm") = 10; +*** --------------- Technical parameters for Nash price algorithm --------------- + +*** Nash adjustment costs (default value around 150). +*** Multiplicator to penalise trade patterns deviations from the last iteration. Involves a trade-off: +*** - when too low, markets are unstable across iterations and cannot clear +*** - when too high, changes in trade patten over iterations are very slow and convergence takes longer +p80_tradeAdj(tradePe) = 80; +p80_tradeAdj("good") = 100; +p80_tradeAdj("perm") = 10; -*LB* parameter for nash price algorithm within the optimization. -p80_etaXp(tradePe) = 0.1; -p80_etaXp("good") = 0.1; -p80_etaXp("perm") = 0.2; +*** Multiplicator to price anticipation +p80_priceAnticipStrength(tradePe) = 0.1; +p80_priceAnticipStrength("good") = 0.1; +p80_priceAnticipStrength("perm") = 0.2; -*LB* parameter for Nash price algorithm between different iterations -p80_etaLT(trade) = 0; -p80_etaLT("perm") = 0.03; -***These parameters are pretty sensitive. If market surpluses diverge, try higher values (up to 1). If surpluses oscillate, try lower values. -p80_etaST(tradePe) = 0.3; -p80_etaST("good") = 0.25; -p80_etaST("perm") = 0.3; +*** Price adjustment factor across iterations (value between zero and one). These parameters are sensitive: +*** - when too low, market surplus may diverge +*** - when too high, market surplus may oscillate -$ifi %banking% == "banking" p80_etaST("perm") = 0.2; !! in banking mode, the permit market reacts more sensitively. -$ifi %emicapregi% == "budget" p80_etaST("perm") = 0.25; !! in budget mode, the permit market reacts more sensitively. +*** short term price ajustment elasticity +p80_shortTermFac(tradePe) = 0.3; +p80_shortTermFac("pebiolc") = 0.8; !! AJS: bio market seems to like this +p80_shortTermFac("peur") = 0.2; !! uranium market is more sensitive, so choose lower etaST +p80_shortTermFac("good") = 0.25; +p80_shortTermFac("perm") = 0.3; +$ifi %banking% == "banking" p80_shortTermFac("perm") = 0.2; !! in banking mode, the permit market reacts more sensitively. +$ifi %emicapregi% == "budget" p80_shortTermFac("perm") = 0.25; !! in budget mode, the permit market reacts more sensitively. -*AJS* bio market seems to like this: -p80_etaST("pebiolc") = 0.8; -***peur market is more sensitive, so choose lower etaST -p80_etaST("peur") = 0.2; +*** long term price ajustment elasticity +p80_longTermFac(trade) = 0; +p80_longTermFac("perm") = 0.03; + +*** --------------- Initialise convergence process parameters --------------- s80_converged = 0; -***initialize some convergence process parameters s80_fadeoutPriceAnticipStartingPeriod = 0; sm_fadeoutPriceAnticip = 1; -*AJS*technical stuff. We want GAMS to import values for the following variables/parameters from the gdx, it would not do that unless you set them a (any) value beforehand. + +*** Negishi weights, not used in Nash +pm_w(regi) = 1; + +*** parameter for spillover externality (aggregated productivity level) +pm_cumEff(t, regi, in) = 100; + +*** auxiliary parameter to track how long the surplus for an item (ttot, trade) had the same sign over iterations +o80_trackSurplusSign(ttot,trade,iteration) $ (not tradeSe(trade)) = 0; + +*** AJS: need to set (any) values for the following variables/parameters so that GAMS imports them from the gdx pm_pvp(ttot,trade)$(ttot.val ge 2005) = NA; p80_pvpFallback(ttot,trade)$(ttot.val ge 2005) = NA; pm_Xport0(ttot,regi,trade)$(ttot.val ge 2005) = NA; @@ -95,21 +94,24 @@ pm_capCumForeign(ttot,regi,teLearn)$(ttot.val ge 2005)=0; qm_co2eqCum.m(regi) = 0; q80_budgetPermRestr.m(regi) = 0; -***read in price paths as fallback option -***p80_pvpFallback(ttot,trade) = 0; +*** read in price paths as fallback option +*** p80_pvpFallback(ttot,trade) = 0; $include "./modules/80_optimization/nash/input/prices_NASH.inc"; -***EMIOPT------------------------------------------------------------------------------ -if ( cm_emiscen eq 6, +*** --------------- Emissions initialisation --------------- +*** emissions, which are part of the climate policy, of other regions (nash relevant) +*** MLB 20140109: initialisation of climate externality is sensitive +*** may be overwritten in nash/presolve +pm_co2eqForeign(t, regi) = (1 - pm_shPerm(t,regi)) * pm_emicapglob(t); !! (1 - emission permit shares) * global emission cap + $ifthen.emiopt %emicapregi% == "none" -*AJS* initialize +if (cm_emiscen = 6, !! budget p80_eoMargEmiCum(regi) = 0; p80_eoMargPermBudg(regi) = 0; -*** ML 20150609 * initialization of permit budget shares (emiopt version, no permit trade) -*** convergence sensitive to initial allocation -*** ML 20161808 * If you change pm_shPerm, ensure that sum(regi, pm_shPerm) = 1; otherwise, a mismatch between global +*** ML 20150609: initialisation of permit budget shares (emiopt version, no permit trade), convergence sensitive to initial allocation +*** ML 20161808: If you change pm_shPerm, ensure that sum(regi, pm_shPerm) = 1; otherwise, a mismatch between global *** and regional budgets will likely disturb results in runs with iterative adjustment - pm_shPerm("2050",regi) = pm_pop("2050",regi)/ sum(regi2, pm_pop("2050",regi2) ); -$endif.emiopt + pm_shPerm("2050",regi) = pm_pop("2050",regi) / sum(regi2, pm_pop("2050",regi2)); ); +$endif.emiopt *** EOF ./modules/80_optimization/nash/datainput.gms diff --git a/modules/80_optimization/nash/declarations.gms b/modules/80_optimization/nash/declarations.gms index acc7f6d463..b251522399 100644 --- a/modules/80_optimization/nash/declarations.gms +++ b/modules/80_optimization/nash/declarations.gms @@ -7,55 +7,54 @@ *** SOF ./modules/80_optimization/nash/declarations.gms parameter -*LB* parameters for ajustments within one iteration. These cause price anticipation -p80_etaXp(all_enty) "Parameter governing price anticipation on commodity markets" +*** parameters for ajustments within one iteration. These cause price anticipation +p80_priceAnticipStrength(all_enty) "Parameter governing price anticipation on commodity markets" +*** parameters for ajustments between differ ent iterations +p80_longTermFac(all_enty) "Long term price ajustment elasticity" +p80_shortTermFac(all_enty) "Short term price ajustment elasticity" +p80_tradeAdj(all_enty) "Adjustment costs for changes of trade pattern between iterations" -*LB* parameters for ajustments between different iterations -p80_etaLT(all_enty) "long term price ajustment elasticity " -p80_etaST(all_enty) "short term price ajustment elasticity" +*** prices and market volumes +p80_pvp_itr(ttot,all_enty,iteration) "Price on commodity markets per iteration" +p80_pvpFallback(ttot,all_enty) "Helper parameter. Price path from input/prices_NASH.inc. Only used if reading prices from gdx fails." -*AJS* adjustment costs between iterations -p80_etaAdj(all_enty) "Adjustment costs for changes of trade pattern between iterations" +p80_regiMarketVolume(ttot,all_regi,all_enty) "Regional market volume of a trade item, used for normalisation [amount of trade item]" +p80_globMarketVolume(ttot,all_enty) "Global market volume of a trade item, used for normalisation [amount of trade item]" +p80_intertempSurplusRevenue(all_enty) "Aggregated intertemporal value of the surplus [T$2017]" +p80_itertempMarketRevenue(all_enty) "Aggregated intertemporal market volume [T$2017]" -***prices -p80_pvp_itr(ttot,all_enty,iteration) "Price on commodity markets per iteration", -p80_pvpFallback(ttot,all_enty) "Helper parameter. Price path from input/prices_NASH.inc. Only used if reading prices from gdx fails.", +*** parameter containing the respective leve l values from last iteration (the first set of values taken from gdx in the first iteration, respectively) +p80_Mport0(tall,all_regi,all_enty) "Imports in last iteration" +p80_surplus(tall,all_enty,iteration) "Surplus on commodity market" +p80_defic_trade(all_enty) "Surplus in monetary terms over all times on commodity markets [T$2017]" +p80_defic_sum(iteration) "Surplus in monetary terms over all times on all commodity markets combined [T$2017] (NOTE: to compare this number with the Negishi defic_sum, divide by around 100)" +p80_defic_sum_rel(iteration) "Surplus monetary value over all times on all commodity markets combined, normalized to consumption [%]" -p80_normalizeLT(all_enty) "Aggregated intertemporal market volume", -p80_normalize0(ttot,all_regi,all_enty) "Normalization parameter for market volume" +*** diagnostic parameters +p80_longTermCorrect(all_enty,iteration) "Long term price correction factor in percent" +p80_shortTermCorrect(tall,all_enty,iteration) "Short term price correction factor in percent" -***parameter containing the respective level values from last iteration (the first set of values taken from gdx in the first iteration, respectively) -p80_Mport0(tall,all_regi,all_enty) "Imports in last iteration" -p80_surplus(tall,all_enty,iteration) "Surplus on commodity market", -p80_defic_trade(all_enty) "Surplus in monetary terms over all times on commodity markets [trillion US$2017]", -p80_defic_sum(iteration) "Surplus in monetary terms over all times on all commodity markets combined [trillion US$2017] (NOTE: to compare this number with the Negishi defic_sum, divide by around 100)", -p80_defic_sum_rel(iteration) "Surplus monetary value over all times on all commodity markets combined, normalized to consumption [%]", +p80_shortTermCorrect_safecopy(tall,all_enty,iteration) "auxiliary parameter to remember short term price correction factor in percent, before new convergence adjustments" +o80_counter_iteration_trade_ttot(ttot,all_enty,iteration) "auxiliary parameter to display in which iteration and for which item (ttot, trade) additional convergence measures were taken" +o80_trackSurplusSign(ttot,all_enty,iteration) "auxiliary parameter to track how long the surplus for an item (ttot, trade) had the same sign over iterations" +o80_SurplusOverTolerance(ttot,all_enty,iteration) "auxiliary parameter to track in which iterations which item surpassed the tolerance (positive/negative)" -*LB* diagnostic parameters -p80_etaLT_correct(all_enty,iteration) "long term price correction factor in percent" -p80_etaST_correct(tall,all_enty,iteration) "short term price correction factor in percent" -p80_etaST_correct_safecopy(tall,all_enty,iteration) "auxiliary parameter to remember short term price correction factor in percent, before new convergence adjustments" -o80_counter_iteration_trade_ttot(ttot,all_enty,iteration) "auxiliary parameter to display in which iteration and for which item (ttot, trade) additional convergence measures were taken" -o80_trackSurplusSign(ttot,all_enty,iteration) "auxiliary parameter to track how long the surplus for an item (ttot, trade) had the same sign over iterations" -o80_SurplusOverTolerance(ttot,all_enty,iteration) "auxiliary parameter to track in which iterations which item surpassed the tolerance (positive/negative)" +p80_surplusMax_iter(all_enty,iteration,tall) "Diagnostics for Nash: Worst residual market surplus until given year, absolute value [TWa, T$, GtC]" +p80_surplusMax2100(all_enty) "Worst residual market surplus until 2100, absolute value [TWa, T$, GtC]" +p80_surplusMaxRel(all_enty,iteration,tall) "Diagnostics for Nash: Worst residual market surplus until given year, in per cent." +p80_surplusMaxTolerance(all_enty) "maximum tolerable residual value of absolute market surplus in 2100." - -p80_surplusMax_iter(all_enty,iteration,tall) "Diagnostics for Nash: Worst residual market surplus until given year, absolute value. [TWa, trillion Dollar, GtC]" -p80_surplusMax2100(all_enty) "Worst residual market surplus until 2100, absolute value. [TWa, trillion Dollar, GtC]" -p80_surplusMaxRel(all_enty,iteration,tall) "Diagnostics for Nash: Worst residual market surplus until given year, in per cent." -p80_surplusMaxTolerance(all_enty) "maximum tolerable residual value of absolute market surplus in 2100." - -p80_taxrev0(tall,all_regi) "vm_taxrev from last iteration" -p80_taxrev_agg(tall,iteration) "vm_taxrev globally from last iteration" +p80_taxrev0(tall,all_regi) "vm_taxrev from last iteration" +p80_taxrev_agg(tall,iteration) "vm_taxrev globally from last iteration" p80_handle(all_regi) "parallel mode handle parameter" p80_repy(all_regi,solveinfo80) "summary report from solver " p80_repy_iteration(all_regi,solveinfo80,iteration) "summary report from solver in iteration" p80_repyLastOptim(all_regi,solveinfo80) "p80_repy from last iteration" -p80_repy_thisSolitr(all_regi,solveinfo80) "p80_repy from the current solitr - only shows results for regions that were run in this solItr" +p80_repy_thisSolitr(all_regi,solveinfo80) "p80_repy from the current solitr - only shows results for regions that were run in this solItr" p80_repy_nashitr_solitr(all_regi,solveinfo80,iteration,sol_itr) "summary report from solver in nash iteration and solver iteration" p80_messageFailedMarket(tall,all_enty) "nash display helper" p80_messageShow(convMessage80) "nash display helper" @@ -65,42 +64,42 @@ p80_curracc(ttot,all_regi) "current account" pm_cumEff(tall,all_regi,all_in) "parameter for spillover externality (aggregated productivity level)" -p80_PriceChangePriceAnticipReg(ttot,all_enty,all_regi) "Price change of a trade good due to the price anticipation effect. [Percent]" -o80_PriceChangePriceAnticipReg(ttot,all_enty,all_regi) "only for display: Price change of a trade good due to price anticipation. If nothing is displayed, all values are <0.1%. [Percent, rounded to 0.1%]" -o80_PriceChangePriceAnticipRegMaxIter(ttot,iteration) "only for display: Largest absolute value of o80_PriceChangePriceAnticipReg until 2100/2150, tracked over iteration. [Percent, rounded to 0.1%]" -p80_DevPriceAnticipReg(ttot,all_enty,all_regi) "Deviation of the yearly monetary export/import expenditure due to price change anticipation effect. [trillion Dollar]" -p80_DevPriceAnticipGlob(ttot,all_enty) "Global sum of p80_DevPriceAnticipReg. [trillion Dollar]" -p80_DevPriceAnticipGlobIter(ttot,all_enty,iteration) "Track p80_DevPriceAnticipGlob over iterations. [trillion Dollar]" -p80_DevPriceAnticipGlobAll(ttot) "p80_DevPriceAnticipGlob summed over all trade goods. [trillion Dollar]" -p80_DevPriceAnticipGlobMax(ttot,all_enty) "Max of p80_DevPriceAnticipGlob until the given year. [trillion Dollar]" -p80_DevPriceAnticipGlobAllMax(ttot) "Max of p80_DevPriceAnticipGlobAll until the given year. [trillion Dollar]" -p80_DevPriceAnticipGlobMax2100Iter(all_enty,iteration) "Track the 2100 value of p80_DevPriceAnticipGlobMax over iterations. [trillion Dollar]" -p80_DevPriceAnticipGlobAllMax2100Iter(iteration) "Track the 2100 value of p80_DevPriceAnticipGlobAllMax over iterations. [trillion Dollar]" - -*EMIOPT relevant -p80_eoMargPermBudg(all_regi) "marginal of permit budget restriction" -p80_eoMargEmiCum(all_regi) "marginal of cumulative emissions" - -p80_eoMargAverage "global average of marginals from nash budget equation" -p80_eoMargDiff(all_regi) "scaled deviation of regional marginals from global average" -p80_eoDeltaEmibudget "total change in permit budget" -p80_eoEmiMarg(all_regi) "weighted marginal utility of emissions" -p80_eoWeights(all_regi) "welfare weights" - -p80_eoMargDiffItr(all_regi,iteration) "scaled deviation of regional marginals from global average" -p80_eoEmibudget1RegItr(all_regi,iteration) "corrected regional permit budgets" -p80_eoEmibudgetDiffAbs(iteration) "convergence indicator" -p80_count "count regions with feasible solution" - -p80_SolNonOpt(all_regi) "solve status" +p80_PriceChangePriceAnticipReg(ttot,all_enty,all_regi) "Price change of a trade good due to the price anticipation effect [Percent]" +o80_PriceChangePriceAnticipReg(ttot,all_enty,all_regi) "only for display: Price change of a trade good due to price anticipation. If nothing is displayed, all values are <0.1% [Percent, rounded to 0.1%]" +o80_PriceChangePriceAnticipRegMaxIter(ttot,iteration) "only for display: Largest absolute value of o80_PriceChangePriceAnticipReg until 2100/2150, tracked over iteration [Percent, rounded to 0.1%]" +p80_DevPriceAnticipReg(ttot,all_enty,all_regi) "Deviation of the yearly monetary export/import expenditure due to price change anticipation effect [T$]" +p80_DevPriceAnticipGlob(ttot,all_enty) "Global sum of p80_DevPriceAnticipReg [T$]" +p80_DevPriceAnticipGlobIter(ttot,all_enty,iteration) "Track p80_DevPriceAnticipGlob over iterations [T$]" +p80_DevPriceAnticipGlobAll(ttot) "p80_DevPriceAnticipGlob summed over all trade goods [T$]" +p80_DevPriceAnticipGlobMax(ttot,all_enty) "Max of p80_DevPriceAnticipGlob until the given year [T$]" +p80_DevPriceAnticipGlobAllMax(ttot) "Max of p80_DevPriceAnticipGlobAll until the given year [T$]" +p80_DevPriceAnticipGlobMax2100Iter(all_enty,iteration) "Track the 2100 value of p80_DevPriceAnticipGlobMax over iterations [T$]" +p80_DevPriceAnticipGlobAllMax2100Iter(iteration) "Track the 2100 value of p80_DevPriceAnticipGlobAllMax over iterations [T$]" + +*** EMIOPT relevant +p80_eoMargPermBudg(all_regi) "marginal of permit budget restriction" +p80_eoMargEmiCum(all_regi) "marginal of cumulative emissions" + +p80_eoMargAverage "global average of marginals from nash budget equation" +p80_eoMargDiff(all_regi) "scaled deviation of regional marginals from global average" +p80_eoDeltaEmibudget "total change in permit budget" +p80_eoEmiMarg(all_regi) "weighted marginal utility of emissions" +p80_eoWeights(all_regi) "welfare weights" + +p80_eoMargDiffItr(all_regi,iteration) "scaled deviation of regional marginals from global average" +p80_eoEmibudget1RegItr(all_regi,iteration) "corrected regional permit budgets" +p80_eoEmibudgetDiffAbs(iteration) "convergence indicator" +p80_count "count regions with feasible solution" + +p80_SolNonOpt(all_regi) "solve status" pm_fuExtrForeign(ttot,all_regi,all_enty,rlf) "foreign fuel extraction" -p80_convNashTaxrev_iter(iteration,ttot,all_regi) "deviation of tax revenue relative to GDP per iteration, thus 0.01 means 1 percent [1]" -p80_convNashObjVal_iter(iteration,all_regi) "deviation of objective value to objective value from last iteration per iteration" +p80_convNashTaxrev_iter(iteration,ttot,all_regi) "deviation of tax revenue relative to GDP per iteration, thus 0.01 means 1 percent [1]" +p80_convNashObjVal_iter(iteration,all_regi) "deviation of objective value to objective value from last iteration per iteration" p80_fadeoutPriceAnticip_iter(iteration) "Helper parameter, describes fadeout of price anticipation during iterations per iteration" $ifthen.cm_implicitQttyTarget not "%cm_implicitQttyTarget%" == "off" -p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) "deviation of current iteration quantity target from target per iteration - relative for total targets, absolute (= share points) for share targets" +p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) "deviation of current iteration quantity target from target per iteration - relative for total targets, absolute (= share points) for share targets" $endif.cm_implicitQttyTarget p80_globalBudget_absDev_iter(iteration) "absolute deviation of global cumulated CO2 emissions budget from target budget" p80_sccConvergenceMaxDeviation_iter(iteration) "max deviation of SCC from last iteration per iteration [percent]" @@ -108,17 +107,19 @@ p80_gmt_conv_iter(iteration) "global mean temperature conve ; positive variable -*AJS* Adjustment costs for Nash trade algorithm. Only non-zero in the Nash_test realization of 80_optimization module. -vm_costAdjNash(ttot,all_regi) "Adjustment costs for deviation from the trade structure of the last iteration." +*** adjustment costs for Nash trade algorithm (only non-zero in the nash realization of 80_optimization module) +vm_costAdjNash(ttot,all_regi) "Adjustment costs for deviation from the trade structure of the last iteration." ; equations -q80_budg_intertemp(all_regi) "interemporal trade balance (Nash mode only)" -q80_costAdjNash(ttot,all_regi) "calculate Nash adjustment costs (of anticipation of the difference in trade pattern, compared to the last iteration), combined for all markets" -q80_budgetPermRestr(all_regi) "constraints regional permit budget to given regional emission budget"; +q80_budg_intertemp(all_regi) "interemporal trade balance (Nash mode only)" +q80_costAdjNash(ttot,all_regi) "calculate Nash adjustment costs (of anticipation of the difference in trade pattern, compared to the last iteration), combined for all markets" +q80_budgetPermRestr(all_regi) "constraints regional permit budget to given regional emission budget" +; scalars -***convergence criteria. if met, the optimization is stopped. Feel free to adjust these to your needs. Denote maximum tolerable deviation from market clearance.(the one for goods is given in million US$2017/yr, the resources in EJ/yr) +*** Convergence criteria: if met, the optimization is stopped. Feel free to adjust these to your needs. +*** Denote maximum tolerable deviation from market clearance (goods in M$2017/yr, resources in EJ/yr) sm_fadeoutPriceAnticip "Helper parameter, describes fadeout of price anticipation during iterations" s80_fadeoutPriceAnticipStartingPeriod "Helper parameter, denotes iteration in which price anticipation fadeout starts" s80_dummy "dummy scalar" @@ -130,7 +131,7 @@ s80_converged "if nash converged, this is 1" s80_runInDebug "Is 1 if regions stayed infeasible in nash and start in debug mode automatically following the parallel mode" /0/ ; -*' defining specific output formats: +*** specific output formats option p80_DevPriceAnticipGlobAll:3:0:1; option p80_DevPriceAnticipGlobAllMax:3:0:1; option o80_PriceChangePriceAnticipReg:1:2:1; diff --git a/modules/80_optimization/nash/equations.gms b/modules/80_optimization/nash/equations.gms index 8bc6b50d25..d92c0e1063 100644 --- a/modules/80_optimization/nash/equations.gms +++ b/modules/80_optimization/nash/equations.gms @@ -7,47 +7,46 @@ *** SOF ./modules/80_optimization/nash/equations.gms *' @equations -*' For Nash solution: intertemporal trade balance must be zero (couple in agricultural trade costs: pvp deflator * net export) + +*' Intertemporal trade balance must be zero (couple in agricultural trade costs). Terms are: +*' 1. foreign assets in 2005 +*' 2. sum over all time steps of (time step duration * net exports in this time step): +*' a) net agricultural exports +*' b) difference of net exports compared to previous iteration, adjusted for price anticipation q80_budg_intertemp(regi).. -0 =e= pm_nfa_start(regi) * pm_pvp("2005","good") - + SUM(ttot$(ttot.val ge 2005), - pm_ts(ttot) - * ( - SUM(trade$(NOT tradeSe(trade) and NOT tradeCap(trade)), - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) * pm_pvp(ttot,trade) - * ( 1 + sm_fadeoutPriceAnticip*p80_etaXp(trade) - * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) - - p80_taxrev0(ttot,regi)$(ttot.val gt 2005)$(sameas(trade,"good")) + vm_taxrev(ttot,regi)$(ttot.val gt 2005)$(sameas(trade,"good")) - ) - / (p80_normalize0(ttot,regi,trade) + sm_eps) - ) + 0 =e= + pm_nfa_start(regi) * pm_pvp("2005","good") !! net value foreign assets in 2005 ++ sum(ttot $ (ttot.val >= 2005), + pm_ts(ttot) !! duration of the time step (average between previous and next time steps) + * ( pm_pvp(ttot,"good") * pm_NXagr(ttot,regi) !! net value agricultural exports + + vm_capacityTradeBalance(ttot,regi) !! trade balance for 24_trade capacity realisation + + sum(trade $ (not tradeSe(trade) and not tradeCap(trade)), + (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) * pm_pvp(ttot,trade) !! net value of exports + * ( 1 + + sm_fadeoutPriceAnticip * p80_priceAnticipStrength(trade) + * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) + - (p80_taxrev0(ttot,regi) - vm_taxrev(ttot,regi)) $ (sameas(trade,"good") and ttot.val > 2005)) + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) ) - + vm_capacityTradeBalance(ttot,regi) - + pm_pvp(ttot,"good") * pm_NXagr(ttot,regi) ) - ); + ) + ); -*' quadratic adjustment costs, penalizing deviations from the trade pattern of the last iteration. -q80_costAdjNash(ttot,regi)$( ttot.val ge cm_startyear ) .. - vm_costAdjNash(ttot,regi) - =e= sum(trade$(NOT tradeSe(trade)), - pm_pvp(ttot,trade) - * p80_etaAdj(trade) - * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) - ) - * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) - ) - / (p80_normalize0(ttot,regi,trade) + sm_eps) - ) -; +*' Quadratic adjustment costs penaliseing deviations of trade volumes from previous iteration. +q80_costAdjNash(ttot,regi) $ (ttot.val >= cm_startyear) .. + vm_costAdjNash(ttot,regi) + =e= sum(trade $ (not tradeSe(trade)), + p80_tradeAdj(trade) + * pm_pvp(ttot,trade) + * power((pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)), 2) + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) + ); -*' link between permit budget and emission budget -q80_budgetPermRestr(regi)$(cm_emiscen=6) .. - sum(ttot$(ttot.val lt sm_endBudgetCO2eq and ttot.val ge cm_startyear), pm_ts(ttot)* vm_perm(ttot,regi)) - + sum(ttot$(ttot.val eq sm_endBudgetCO2eq),pm_ts(ttot)/2 * (vm_perm(ttot,regi))) =l= - pm_budgetCO2eq(regi) - sum(ttot $((ttot.val ge 2005) and (ttot.val lt cm_startyear)), pm_ts(ttot)* vm_co2eq(ttot,regi)); +*' link between permit budget and emission budget +q80_budgetPermRestr(regi) $ (cm_emiscen = 6) .. + sum(ttot$(ttot.val < sm_endBudgetCO2eq and ttot.val >= cm_startyear), pm_ts(ttot) * vm_perm(ttot,regi)) + + sum(ttot$(ttot.val = sm_endBudgetCO2eq), pm_ts(ttot) / 2 * vm_perm(ttot,regi)) =l= + pm_budgetCO2eq(regi) - sum(ttot $((ttot.val >= 2005) and (ttot.val < cm_startyear)), pm_ts(ttot) * vm_co2eq(ttot,regi)); *' @stop *** EOF ./modules/80_optimization/nash/equations.gms diff --git a/modules/80_optimization/nash/postsolve.gms b/modules/80_optimization/nash/postsolve.gms index 982ef409f6..5779740ccf 100644 --- a/modules/80_optimization/nash/postsolve.gms +++ b/modules/80_optimization/nash/postsolve.gms @@ -7,290 +7,255 @@ *** SOF ./modules/80_optimization/nash/postsolve.gms - -*ML*2015-02-04* calculate current account -*LB* needed for decomposition script -p80_curracc(ttot, regi) = SUM(trade$(NOT tradeSe(trade)), pm_pvp(ttot,trade)/ max(pm_pvp(ttot,"good"),sm_eps) * (vm_Xport.l(ttot,regi,trade)- vm_Mport.l(ttot,regi,trade)) ); - -p80_taxrev0(ttot,regi)$( (ttot.val ge max(2010,cm_startyear)) and (pm_SolNonInfes(regi) eq 1) ) = vm_taxrev.l(ttot,regi); - -*AJS*update normalization paramaters, take values from last iteration for regions that were not solved optimally -p80_normalize0(ttot,regi,"good")$(ttot.val ge 2005) = max(vm_cons.l(ttot,regi)$(pm_SolNonInfes(regi) eq 1) + p80_normalize0(ttot,regi,"good")$(pm_SolNonInfes(regi) eq 0),sm_eps); -*ML*normalize permit trade corrections to consumption or positive cap path instead of emissions, as those may be negative -p80_normalize0(ttot,regi,"perm")$(ttot.val ge 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050")) , sm_eps); -p80_normalize0(ttot,regi,tradePe)$(ttot.val ge 2005) = max(0.5 * (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe))$(pm_SolNonInfes(regi) eq 1) - + p80_normalize0(ttot,regi,tradePe)$(pm_SolNonInfes(regi) eq 0) ,sm_eps); - - -***calculate residual surplus on the markets -loop(ttot$(ttot.val ge 2005), - loop(trade$(NOT tradeSe(trade)), - p80_surplus(ttot,trade,iteration) = sum(regi, (vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade))$(pm_SolNonInfes(regi) eq 1) - + (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade) )$(pm_SolNonInfes(regi) eq 0) ); - ); -); - -*' calculate both the size of the price change due to the price change anticipation effect in percent, as well as -*' the deviation of the yearly monetary export/import expenditure due to the price change anticipation effect: -loop(ttot$(ttot.val ge 2005), - loop(trade$(NOT tradeSe(trade)), - loop(regi, - p80_PriceChangePriceAnticipReg(ttot,trade,regi) = 100 * - ( sm_fadeoutPriceAnticip*p80_etaXp(trade) - * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade)) - - p80_taxrev0(ttot,regi)$(ttot.val gt 2005)$(sameas(trade,"good")) + vm_taxrev.l(ttot,regi)$(ttot.val gt 2005)$(sameas(trade,"good")) - ) - / (p80_normalize0(ttot,regi,trade) + sm_eps) - ) - ; - p80_DevPriceAnticipReg(ttot,trade,regi) = - ( vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade) ) - * pm_pvp(ttot,trade) / pm_pvp(ttot,"good") - * p80_PriceChangePriceAnticipReg(ttot,trade,regi) - ; - ); - p80_DevPriceAnticipGlob(ttot,trade) = sum(regi, abs( p80_DevPriceAnticipReg(ttot,trade,regi) ) ); +***------------------------------------------------------------------------------ +*' #### Price corrections to improve convergence of next iteration +***------------------------------------------------------------------------------ + +*' The objective of the nash optimisation is to find a set of global prices so that the market for each trade item clears. +*' If a certain trade has a positive surplus, it means that the global price is too high and should be reduced to obtain convergence. +*' The price corrections are given by price reduction = price adjustment elasticity * market surplus / market volume +*' for example p80_shortTermCorrect = p80_shortTermFac * p80_surplus / p80_regiMarketVolume + +*' Compute market volume for different trades +loop(ttot $ (ttot.val >= 2005), +*** calculation has to match initialisation in nash/preloop + loop(regi $ (pm_SolNonInfes(regi) = 1), !! only for optimally solved regions, others keep values from last iteration + p80_regiMarketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_regiMarketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; ); -p80_DevPriceAnticipGlobAll(ttot) = sum(trade$(NOT tradeSe(trade)), p80_DevPriceAnticipGlob(ttot,trade)); + p80_regiMarketVolume(ttot,regi,"perm") = abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050")); !! has to be positive for normalisation + p80_regiMarketVolume(ttot,regi,trade) = max(sm_eps, p80_regiMarketVolume(ttot,regi,trade)); !! ensure market volume is positive + p80_globMarketVolume(ttot,trade) = sum(regi, p80_regiMarketVolume(ttot,regi,trade)); !! compute global market volume ); -*' calculate maximum of p80_DevPriceAnticipGlob -p80_DevPriceAnticipGlobMax(ttot,trade)$((ttot.val ge cm_startyear) AND (NOT tradeSe(trade))) = - smax(ttot2$(ttot2.val ge cm_startyear AND ttot2.val le ttot.val), p80_DevPriceAnticipGlob(ttot2,trade) ) -; -*' calculate maximum of p80_DevPriceAnticipGlobAll -p80_DevPriceAnticipGlobAllMax(ttot)$((ttot.val ge cm_startyear)) = - smax(ttot2$(ttot2.val ge cm_startyear AND ttot2.val le ttot.val), p80_DevPriceAnticipGlobAll(ttot2) ) -; - -p80_DevPriceAnticipGlobIter(ttot,trade,iteration)$((ttot.val ge cm_startyear) AND (NOT tradeSe(trade))) = p80_DevPriceAnticipGlob(ttot,trade); -p80_DevPriceAnticipGlobMax2100Iter(trade,iteration)$(NOT tradeSe(trade)) = p80_DevPriceAnticipGlobMax("2100",trade); -p80_DevPriceAnticipGlobAllMax2100Iter(iteration) = p80_DevPriceAnticipGlobAllMax("2100"); - +loop(trade $ (not tradeSe(trade)), +*' Calculate residual surplus on the markets: trade items that are exported but not imported + p80_surplus(ttot,trade,iteration) $ (ttot.val >= 2005) = sum(regi, + (vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade)) $ (pm_SolNonInfes(regi) = 1) + + (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade) ) $ (pm_SolNonInfes(regi) = 0) ); + +*' Long term price correction takes into account the aggregated intertemporal market revenue (instead of volume) defined by +*' market revenue = price * duration * market yearly volume (or pm_pvp * pm_ts * p80_regiMarketVolume) + p80_intertempSurplusRevenue(trade) = sum(ttot $ (ttot.val >= cm_startyear), + pm_pvp(ttot,trade) * pm_ts(ttot) * p80_surplus(ttot,trade,iteration)); + p80_itertempMarketRevenue(trade) = max(sm_eps, sum((ttot, regi) $ (ttot.val >= 2005), + pm_pvp(ttot,trade) * pm_ts(ttot) * p80_regiMarketVolume(ttot,regi,trade))); + + p80_longTermCorrect(trade,iteration) = + p80_longTermFac(trade) * p80_intertempSurplusRevenue(trade) / p80_itertempMarketRevenue(trade); + +*' Short term price correction takes into account the market surplus volume of a single time step +*' For permit and primary energy trade, price anticipation impacts the price correction + p80_shortTermCorrect(ttot,trade,iteration) $ (ttot.val >= 2005) = + p80_shortTermFac(trade) * p80_surplus(ttot,trade,iteration) / p80_globMarketVolume(ttot,trade); + + p80_shortTermCorrect(ttot,"perm",iteration) $ (ttot.val >= 2005) = p80_shortTermCorrect(ttot,"perm",iteration) + * ((1-sm_fadeoutPriceAnticip) + sm_fadeoutPriceAnticip * sqrt(pm_pvp(ttot,"good") / pm_pvp("2100","good"))) + * (sm_fadeoutPriceAnticip + (1-sm_fadeoutPriceAnticip) * (pm_pvp(ttot,"good") / pm_pvp("2040","good"))); + p80_shortTermCorrect(ttot,tradePe(trade),iteration) $ (ttot.val >= 2005) = p80_shortTermCorrect(ttot,trade,iteration) + * (sm_fadeoutPriceAnticip + (1-sm_fadeoutPriceAnticip) * (pm_pvp(ttot,trade) / pm_pvp("2050",trade))); + + p80_shortTermCorrect_safecopy(ttot,trade,iteration) = p80_shortTermCorrect(ttot,trade,iteration); !! copy of initial values +); -*' For display of price change p80_PriceChangePriceAnticipReg, round to 0.1% -*o80_PriceChangePriceAnticipReg(ttot,trade,regi) = p80_PriceChangePriceAnticipReg(ttot,trade,regi); -o80_PriceChangePriceAnticipReg(ttot,trade,regi) = round(p80_PriceChangePriceAnticipReg(ttot,trade,regi),1); -*' determine largest price change in p80_PriceChangePriceAnticipReg -o80_PriceChangePriceAnticipRegMaxIter("2100",iteration) = smax( (ttot,trade,regi)$(ttot.val le 2100) , abs(o80_PriceChangePriceAnticipReg(ttot,trade,regi) ) ); -o80_PriceChangePriceAnticipRegMaxIter("2150",iteration) = smax( (ttot,trade,regi)$(ttot.val ge 2110) , abs(o80_PriceChangePriceAnticipReg(ttot,trade,regi) ) ); +*' If the surplus remains over several iterations, increase the price correction terms +loop(ttot $ (ttot.val >= 2005), + loop(trade $ (tradePe(trade) or sameas(trade,"good")), + if(abs(p80_surplus(ttot,trade,iteration)) <= p80_surplusMaxTolerance(trade), + o80_trackSurplusSign(ttot,trade,iteration) = 0; !! reset counter if surplus is within tolerance range + else !! if surplus is outside tolerance range + o80_SurplusOverTolerance(ttot,trade,iteration) = Sign(p80_surplus(ttot,trade,iteration)); !! track the sign of the surplus -display - p80_DevPriceAnticipGlob, - p80_DevPriceAnticipGlobMax, - p80_DevPriceAnticipGlobAllMax, - p80_DevPriceAnticipGlobMax2100Iter, - p80_DevPriceAnticipGlobAllMax2100Iter, - p80_DevPriceAnticipGlobAll, - o80_PriceChangePriceAnticipReg - o80_PriceChangePriceAnticipRegMaxIter -; - -***calculate aggregated intertemporal market volumes - used in calculation of price corrections later on -loop(trade$(NOT tradeSe(trade)), - p80_normalizeLT(trade) = sum(ttot$(ttot.val ge 2005), sum(regi, pm_pvp(ttot,trade) * pm_ts(ttot) * p80_normalize0(ttot,regi,trade) )); - if (p80_normalizeLT(trade) = 0, p80_normalizeLT(trade) = sm_eps); - ); - -*LB* calculate price correction terms -p80_etaLT_correct(trade,iteration)$(NOT tradeSe(trade)) = - p80_etaLT(trade) * - sum(ttot2$(ttot2.val ge cm_startyear), pm_pvp(ttot2,trade) * pm_ts(ttot2) * p80_surplus(ttot2,trade,iteration) ) - / p80_normalizeLT(trade); - -p80_etaST_correct(ttot,trade,iteration)$((ttot.val ge 2005) AND (NOT tradeSe(trade))) = - p80_etaST(trade) - * ( ( (1-sm_fadeoutPriceAnticip) + sm_fadeoutPriceAnticip * sqrt(pm_pvp(ttot,"good")/pm_pvp("2100","good")) )$(sameas(trade,"perm")) + 1$(NOT sameas(trade,"perm")) ) - * ((sm_fadeoutPriceAnticip + (1-sm_fadeoutPriceAnticip) * (pm_pvp(ttot,"good")/pm_pvp('2040',"good")) )$(sameas(trade,"perm")) + 1$(NOT sameas(trade,"perm")) ) - * ((sm_fadeoutPriceAnticip + (1-sm_fadeoutPriceAnticip) * (pm_pvp(ttot,trade)/pm_pvp('2050',trade)) )$(tradePe(trade)) + 1$(NOT tradePe(trade)) ) - * p80_surplus(ttot,trade,iteration) - / max(sm_eps , sum(regi, p80_normalize0(ttot,regi,trade))); - -*RP* add a stronger push to the price adjustment if convergence doesn't happen for an extended amount of iterations: -p80_etaST_correct_safecopy(ttot,trade,iteration)$(NOT tradeSe(trade)) = p80_etaST_correct(ttot,trade,iteration); !! first make a copy of the initial adjustment values - -*RP* track sign of the surplus -if(iteration.val > 2, - loop(ttot$(ttot.val ge 2005), - loop(trade$(tradePe(trade) OR sameas(trade,"good") ), - if( abs(p80_surplus(ttot,trade,iteration)) gt p80_surplusMaxTolerance(trade), - o80_SurplusOverTolerance(ttot,trade,iteration) = Sign(p80_surplus(ttot,trade,iteration) ); + if(iteration.val > 2, + if(Sign(p80_surplus(ttot,trade,iteration)) = Sign(p80_surplus(ttot,trade,iteration-1)), !! if deviation is in the same direction as in previous iteration + o80_trackSurplusSign(ttot,trade,iteration) = 1 + o80_trackSurplusSign(ttot,trade,iteration-1); + else + o80_trackSurplusSign(ttot,trade,iteration) = 0; !! reset counter if sign changed ); - ); - ); -); + ); -*RP* track continued surplusses with the same sign (to show where convergence is too slow) -if(iteration.val > 2, - loop(ttot$(ttot.val ge 2005), - loop(trade$(tradePe(trade) OR sameas(trade,"good") ), - if( ( Sign(p80_surplus(ttot,trade,iteration) ) eq Sign(p80_surplus(ttot,trade,iteration-1) ) ) - AND ( abs(p80_surplus(ttot,trade,iteration)) gt p80_surplusMaxTolerance(trade) ) , - o80_trackSurplusSign(ttot,trade,iteration) = o80_trackSurplusSign(ttot,trade,iteration-1) +1; - else - o80_trackSurplusSign(ttot,trade,iteration) = 0; +*** After 10 iterations, if the surplus keeps the same sign, apply a multiplicative factor on price correction to speed up convergence +*** The multiplicative factor follows this linear curve: https://www.desmos.com/calculator/xb1uv6ly9h + if(iteration.val > 10 + o80_trackSurplusSign(ttot,trade,iteration), + o80_counter_iteration_trade_ttot(ttot,trade,iteration) = max(1, o80_trackSurplusSign(ttot,trade,iteration) - 3); + p80_shortTermCorrect(ttot,trade,iteration) = p80_shortTermCorrect(ttot,trade,iteration) * o80_counter_iteration_trade_ttot(ttot,trade,iteration); ); - ); + ); !! if surplus is outside tolerance range + ); !! trade +); !! ttot + + +*** calculate prices for next iteration, prevent prices from turning negative by limiting extreme prices corrections +p80_pvp_itr(ttot,trade,iteration+1) $ (ttot.val >= cm_startyear and not tradeSe(trade)) = + pm_pvp(ttot,trade) * max(0.05, 1 - p80_longTermCorrect(trade,iteration) - p80_shortTermCorrect(ttot,trade,iteration)); + +*** AJS: feed updated prices and quantities into the next iteration, ML: adjustments in case of infeasibilities (increase import) +loop(trade $ (not tradeSe(trade)), + loop(ttot $ (ttot.val >= cm_startyear), + pm_pvp(ttot,trade) = p80_pvp_itr(ttot,trade,iteration+1); + pm_Xport0(ttot,regi,trade) $ (pm_SolNonInfes(regi) = 1) = vm_Xport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) $ (pm_SolNonInfes(regi) = 1) = vm_Mport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) $ (pm_SolNonInfes(regi) = 0) = vm_Mport.l(ttot,regi,trade) * 1.2; ); ); +p80_taxrev0(ttot,regi) $ (ttot.val >= max(2010, cm_startyear) and pm_SolNonInfes(regi) = 1) = vm_taxrev.l(ttot,regi); -if(iteration.val > 15, - loop(ttot$(ttot.val ge 2005), - loop(trade$(tradePe(trade) OR sameas(trade,"good")), - if( abs(p80_surplus(ttot,trade,iteration)) gt p80_surplusMaxTolerance(trade) , - if( ( abs( sum(iteration2$( (iteration2.val le iteration.val) AND (iteration2.val ge (iteration.val - 4))), - p80_surplus(ttot,trade,iteration2) !! this sum should ensure the additional price adjustment only happens if the surplus was always off the same sign - ) - ) ge ( 5 * p80_surplusMaxTolerance(trade) ) ) AND ( o80_trackSurplusSign(ttot,trade,iteration) ge 5 ) , !! check if surplus was out of the target range for 5 consecutive iterations - p80_etaST_correct(ttot,trade,iteration) = 4 * p80_etaST_correct(ttot,trade,iteration); - o80_counter_iteration_trade_ttot(ttot,trade,iteration) = 1; - - if(iteration.val gt 20, !! only start checking if a stronger push is necessary a few iterations later, so that step 1 could potentially show an effect - if( ( abs( sum(iteration2$( (iteration2.val le iteration.val) AND (iteration2.val ge (iteration.val - 9))), - p80_surplus(ttot,trade,iteration2) - ) - ) ge ( 10 * p80_surplusMaxTolerance(trade)) ) AND ( o80_trackSurplusSign(ttot,trade,iteration) ge 10 ), !! check if surplus was out of the target range for 10 consecutive iterations - p80_etaST_correct(ttot,trade,iteration) = 2 * p80_etaST_correct(ttot,trade,iteration); - o80_counter_iteration_trade_ttot(ttot,trade,iteration) = 2; - - if(iteration.val gt 25, !! only start checking if a stronger push is necessary a few iterations later, so that step 1&2 could potentially show an effect - if( ( abs( sum(iteration2$( (iteration2.val le iteration.val) AND (iteration2.val ge (iteration.val - 14))), - p80_surplus(ttot,trade,iteration2) - ) - ) ge ( 15 * p80_surplusMaxTolerance(trade)) ) AND ( o80_trackSurplusSign(ttot,trade,iteration) ge 15 ), !! check if surplus was out of the target range for 15 consecutive iterations - p80_etaST_correct(ttot,trade,iteration) = 2 * p80_etaST_correct(ttot,trade,iteration); - o80_counter_iteration_trade_ttot(ttot,trade,iteration) = 3; - ); - ); - ); - ); - ); - ); - ); !! trade - ); !! ttot -); !! iteration>15 - - -***calculate prices for next iteration -p80_pvp_itr(ttot,trade,iteration+1)$((ttot.val ge cm_startyear) AND (NOT tradeSe(trade))) = - pm_pvp(ttot,trade) - * max(0.05, !! prevent prices from turning negative by limiting extreme prices corrections - (1 - p80_etaLT_correct(trade,iteration) - - p80_etaST_correct(ttot,trade,iteration) - ) +loop(trade $ (not tradeSe(trade)), +*** AJS: calculate maximum residual surplus on markets, absolute and relative + loop(ttot $ (ttot.val >= cm_startyear), + p80_surplusMax_iter(trade,iteration,ttot) = smax(ttot2 $ (ttot2.val >= cm_startyear and ttot2.val <= ttot.val), abs(p80_surplus(ttot2,trade,iteration))); + p80_surplusMaxRel(trade,iteration,ttot) = 100 * smax(ttot2 $ (ttot2.val >= cm_startyear and ttot2.val <= ttot.val), abs(p80_surplus(ttot2,trade,iteration)) / p80_globMarketVolume(ttot2,trade)); + ); + p80_surplusMax2100(trade) = p80_surplusMax_iter(trade,iteration,"2100"); + +*** convergence indicators + p80_defic_trade(trade) = 1 / pm_pvp("2005","good") * + sum(ttot $ (ttot.val >= 2005), + pm_ts(ttot) * ( + abs(p80_surplus(ttot,trade,iteration)) * pm_pvp(ttot,trade) + + sum(regi, abs(p80_taxrev0(ttot,regi)) * pm_pvp(ttot,"good")) $ (sameas(trade,"good") and ttot.val >= max(2010,cm_startyear) ) + + sum(regi, abs(vm_costAdjNash.l(ttot,regi)) * pm_pvp(ttot,"good")) $ (sameas(trade,"good") and (ttot.val >= 2005) ) ) - ; - -*AJS* feed updated prices and quantities into the next iteration: -*ML* adjustments in case of infeasibilities (increase import) -loop(trade$(NOT tradeSe(trade)), - loop(regi, - loop(ttot$(ttot.val ge cm_startyear), - pm_pvp(ttot,trade) = p80_pvp_itr(ttot,trade,iteration+1); - pm_Xport0(ttot,regi,trade)$(pm_SolNonInfes(regi) eq 1) = vm_Xport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade)$(pm_SolNonInfes(regi) eq 1) = vm_Mport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade)$(pm_SolNonInfes(regi) eq 0) = 1.2 * vm_Mport.l(ttot,regi,trade); - ); ); ); +p80_defic_sum("1") = 1; +p80_defic_sum(iteration) = sum(trade $ (not tradeSe(trade)), p80_defic_trade(trade)); +p80_defic_sum_rel(iteration) = 100 * p80_defic_sum(iteration) / (p80_itertempMarketRevenue("good") / pm_pvp("2005","good")); -***some diagnostic output: -p80_taxrev_agg(ttot,iteration)$(ttot.val ge 2005) = sum(regi,vm_taxrev.l(ttot,regi)); -*AJS* calculate maximum residual surplusses on markets -p80_surplusMax_iter(trade,iteration,ttot)$((ttot.val ge cm_startyear) AND (NOT tradeSe(trade))) = smax(ttot2$(ttot2.val ge cm_startyear AND ttot2.val le ttot.val), abs(p80_surplus(ttot2,trade,iteration))); +*** adjust parameters for next iteration +*** Decide on when to fade out price anticipation terms if markets are reasonably cleared (doing this too early leads to diverging markets) +if( smax(tradePe, p80_surplusMax_iter(tradePe,iteration,"2150")) < 0.5 + and p80_surplusMax_iter("good",iteration,"2150") < 1 + and p80_surplusMax_iter("perm",iteration,"2150") < 1 + and s80_fadeoutPriceAnticipStartingPeriod = 0, !! as long as we are not fading out already + s80_fadeoutPriceAnticipStartingPeriod = iteration.val; +); -***from this, relative residual surplusses. -p80_surplusMaxRel(trade,iteration,ttot)$((ttot.val ge cm_startyear) AND (NOT tradeSe(trade))) = 100 * smax(ttot2$(ttot2.val ge cm_startyear AND ttot2.val le ttot.val), abs(p80_surplus(ttot2,trade,iteration)) / sum(regi, p80_normalize0(ttot2,regi,trade))); +if(s80_fadeoutPriceAnticipStartingPeriod > 0, + sm_fadeoutPriceAnticip = 0.7 ** (iteration.val - s80_fadeoutPriceAnticipStartingPeriod + 1); +); +display s80_fadeoutPriceAnticipStartingPeriod, sm_fadeoutPriceAnticip; -p80_surplusMax2100(trade)$(NOT tradeSe(trade)) = p80_surplusMax_iter(trade,iteration,"2100"); +***------------------------------------------------------------------------------ +*' #### Output and monitoring +***------------------------------------------------------------------------------ -***convergence indicators -loop(trade$(NOT tradeSe(trade)), - p80_defic_trade(trade) = 1/pm_pvp("2005","good") * - sum(ttot$(ttot.val ge 2005), pm_ts(ttot) * ( - abs(p80_surplus(ttot,trade,iteration)) * pm_pvp(ttot,trade) - + sum(regi, abs(p80_taxrev0(ttot,regi)) * pm_pvp(ttot,"good"))$(sameas(trade,"good") and (ttot.val ge max(2010,cm_startyear)) ) - + sum(regi, abs(vm_costAdjNash.l(ttot,regi)) * pm_pvp(ttot,"good"))$(sameas(trade,"good") and (ttot.val ge 2005) ) +*** ML 2015-02-04: calculate current account, LB: needed for cost decomposition script +p80_curracc(ttot, regi) = sum(trade $ (not tradeSe(trade)), + pm_pvp(ttot,trade) / max(pm_pvp(ttot,"good"), sm_eps) * (vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade))); - ) - ); -); -p80_defic_sum("1") = 1; -p80_defic_sum(iteration) = sum(trade$(NOT tradeSe(trade)), p80_defic_trade(trade)); -p80_defic_sum_rel(iteration) = 100 * p80_defic_sum(iteration) / (p80_normalizeLT("good")/pm_pvp("2005","good")); +*** diagnostic output: vm_taxrev globally from last iteration +p80_taxrev_agg(ttot,iteration) $ (ttot.val >= 2005) = sum(regi, vm_taxrev.l(ttot,regi)); +*** save all FE prices across sectors and markets [tr$2005/TWa] across iterations +pm_FEPrice_iter(iteration,t,regi,enty,sector,emiMkt) = pm_FEPrice(t,regi,enty,sector,emiMkt); -***adjust parameters for next iteration -***Decide on when to fade out price anticipation terms (doing this too early leads to diverging markets) -***if markets are reasonably cleared -if( (smax(tradePe,p80_surplusMax_iter(tradePe,iteration,'2150')) lt (10 * 0.05)) - AND ( p80_surplusMax_iter("good",iteration,'2150') lt (10 * 0.1) ) !! - AND ( p80_surplusMax_iter("perm",iteration,'2150') lt (5 * 0.2) ) - AND (s80_fadeoutPriceAnticipStartingPeriod eq 0), !! as long as we are not fading out already - s80_fadeoutPriceAnticipStartingPeriod = iteration.val; -); -***if thats the case, then start to fade out anticipation terms - begin second phase -if(s80_fadeoutPriceAnticipStartingPeriod ne 0, - sm_fadeoutPriceAnticip = 0.7**(iteration.val - s80_fadeoutPriceAnticipStartingPeriod + 1); -); -display s80_fadeoutPriceAnticipStartingPeriod, sm_fadeoutPriceAnticip; +*' calculate both the size of the price change due to the price change anticipation effect in percent +*' and the deviation of the yearly monetary export/import expenditure due to the price change anticipation effect +loop(ttot $ (ttot.val >= 2005), + loop(trade $ (not tradeSe(trade)), + p80_PriceChangePriceAnticipReg(ttot,trade,regi) = 100 * + sm_fadeoutPriceAnticip * p80_priceAnticipStrength(trade) + * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade)) + - (p80_taxrev0(ttot,regi) - vm_taxrev.l(ttot,regi)) $ (ttot.val > 2005 and sameas(trade,"good")) + ) + / p80_regiMarketVolume(ttot,regi,trade); + + p80_DevPriceAnticipReg(ttot,trade,regi) = + ( vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade) ) + * pm_pvp(ttot,trade) / pm_pvp(ttot,"good") + * p80_PriceChangePriceAnticipReg(ttot,trade,regi); -*** Save FE prices in each iteration for easier monitoring -pm_FEPrice_iter(iteration,t,regi,enty,sector,emiMkt) = - pm_FEPrice(t,regi,enty,sector,emiMkt); + p80_DevPriceAnticipGlob(ttot,trade) = sum(regi, abs( p80_DevPriceAnticipReg(ttot,trade,regi) ) ); + ); + p80_DevPriceAnticipGlobAll(ttot) = sum(trade $ (not tradeSe(trade)), p80_DevPriceAnticipGlob(ttot,trade)); +); !! end of ttot loop + + +*' calculate maximum of p80_DevPriceAnticipGlob +p80_DevPriceAnticipGlobMax(ttot,trade) $ (ttot.val >= cm_startyear and not tradeSe(trade)) = + smax(ttot2 $ (ttot2.val >= cm_startyear and ttot2.val <= ttot.val), p80_DevPriceAnticipGlob(ttot2,trade) ); +*' calculate maximum of p80_DevPriceAnticipGlobAll +p80_DevPriceAnticipGlobAllMax(ttot) $ (ttot.val >= cm_startyear) = + smax(ttot2 $ (ttot2.val >= cm_startyear and ttot2.val <= ttot.val), p80_DevPriceAnticipGlobAll(ttot2) ); + +p80_DevPriceAnticipGlobIter(ttot,trade,iteration) $ (ttot.val >= cm_startyear and not tradeSe(trade)) = p80_DevPriceAnticipGlob(ttot,trade); +p80_DevPriceAnticipGlobMax2100Iter(trade,iteration) $ (not tradeSe(trade)) = p80_DevPriceAnticipGlobMax("2100",trade); +p80_DevPriceAnticipGlobAllMax2100Iter(iteration) = p80_DevPriceAnticipGlobAllMax("2100"); -***Decide, on whether to end iterating now. if any of the following criteria (contained in the set convMessage80(surplus,infes,nonopt)) is not met, s80_bool is set to 0, and the convergence process is NOT stopped -***reset some indicators -s80_bool=1; +*' For display of price change p80_PriceChangePriceAnticipReg, round to 0.1% +o80_PriceChangePriceAnticipReg(ttot,trade,regi) = round(p80_PriceChangePriceAnticipReg(ttot,trade,regi), 1); + +*' determine largest price change in p80_PriceChangePriceAnticipReg +o80_PriceChangePriceAnticipRegMaxIter("2100",iteration) = smax( (ttot,trade,regi) $ (ttot.val <= 2100) , abs(o80_PriceChangePriceAnticipReg(ttot,trade,regi) ) ); +o80_PriceChangePriceAnticipRegMaxIter("2150",iteration) = smax( (ttot,trade,regi) $ (ttot.val >= 2110) , abs(o80_PriceChangePriceAnticipReg(ttot,trade,regi) ) ); + +display + p80_DevPriceAnticipGlob, + p80_DevPriceAnticipGlobMax, + p80_DevPriceAnticipGlobAllMax, + p80_DevPriceAnticipGlobMax2100Iter, + p80_DevPriceAnticipGlobAllMax2100Iter, + p80_DevPriceAnticipGlobAll, + o80_PriceChangePriceAnticipReg + o80_PriceChangePriceAnticipRegMaxIter +; + + +***------------------------------------------------------------------------------ +*' #### Convergence criteria +***------------------------------------------------------------------------------ +*** Iterations continue if any of the convergence criteria is not met +*** Criteria are contained in convMessage80(surplus,infes,nonopt), and continuing means that s80_bool = 0 + +*** reset indicators +s80_bool = 1; !! assume all the criteria are met p80_messageShow(convMessage80) = NO; p80_messageFailedMarket(ttot,all_enty) = NO; -***criterion ""surplus": are we converged yet? -loop(trade$(NOT tradeSe(trade)), - if(p80_surplusMax_iter(trade,iteration,"2100") gt p80_surplusMaxTolerance(trade), - s80_bool=0; - p80_messageShow("surplus") = YES; - loop(ttot$((ttot.val ge cm_startyear) and (ttot.val le 2100)), - if( (abs(p80_surplus(ttot,trade,iteration)) gt p80_surplusMaxTolerance(trade) ), - p80_messageFailedMarket(ttot,trade) = YES; - ); +*' criterion "surplus": did trade converge yet? +loop(trade $ (not tradeSe(trade)), + if(p80_surplusMax_iter(trade,iteration,"2100") > p80_surplusMaxTolerance(trade), + s80_bool = 0; + p80_messageShow("surplus") = YES; + loop(ttot $ (ttot.val >= cm_startyear and ttot.val <= 2100), + if(abs(p80_surplus(ttot,trade,iteration)) > p80_surplusMaxTolerance(trade), + p80_messageFailedMarket(ttot,trade) = YES; ); + ); ); - if(p80_surplusMax_iter(trade,iteration,"2150") gt 10 * p80_surplusMaxTolerance(trade), - s80_bool=0; - p80_messageShow("surplus") = YES; - loop(ttot$((ttot.val ge cm_startyear) and (ttot.val gt 2100)), - if( (abs(p80_surplus(ttot,trade,iteration)) gt p80_surplusMaxTolerance(trade) ), - p80_messageFailedMarket(ttot,trade) = YES; - ); + if(p80_surplusMax_iter(trade,iteration,"2150") > 10 * p80_surplusMaxTolerance(trade), + s80_bool = 0; + p80_messageShow("surplus") = YES; + loop(ttot $ (ttot.val >= cm_startyear and ttot.val > 2100), + if(abs(p80_surplus(ttot,trade,iteration)) > p80_surplusMaxTolerance(trade), + p80_messageFailedMarket(ttot,trade) = YES; ); - ); + ); + ); ); -*** critertion "infes": is any region neither optimal nor intermediate non-optimal -> then it is infeasible +*' critertion "infes": is any region neither optimal nor intermediate non-optimal? then it is infeasible loop(regi, - if( (p80_repy(regi,'modelstat') ne 2) and (p80_repy(regi,'modelstat') ne 7), !! 2 is optimal, 7 nonopt, + if(p80_repy(regi,'modelstat') ne 2 and p80_repy(regi,'modelstat') ne 7, !! 2 is optimal, 7 nonopt, s80_bool = 0; p80_messageShow("infes") = YES; ); -*** critertion "nonopt": The next lines are a workaround for the status 7 -*** problem. If the objective value does not differ too much from the last known -*** optimal solution, accept this solution as if it were optimal. + +*' critertion "nonopt": +*** The next lines are a workaround for the status 7 problem +*** If the objective value does not differ too much from the last known optimal solution, accept this solution as if it were optimal p80_convNashObjVal_iter(iteration,regi) = p80_repy(regi,'objval') - p80_repyLastOptim(regi,'objval'); - if (1 le iteration.val, - !! no last iteration if this is the first; NA value in p80_repyLastOptim is - !! sticky, so test this separately - if ( p80_repy(regi,'modelstat') eq 7 - !! The 1E-4 are quite arbitrary. One should do more research on how - !! the solution differs over iteration when status 7 occurs. - AND p80_convNashObjVal_iter(iteration,regi) lt - 1e-4, + if(1 <= iteration.val, !! no last iteration if this is the first; NA value in p80_repyLastOptim is sticky, so test this separately + if(p80_repy(regi,'modelstat') = 7 and p80_convNashObjVal_iter(iteration,regi) < - 1e-4, + !! 1E-4 is quite arbitrary. One should do more research on how the solution differs over iteration when status 7 occurs. s80_bool = 0; p80_messageShow("nonopt") = YES; display "Not all regions were status 2 in the last iteration. The deviation of the objective function from the last optimal solution is too large to be accepted:"; @@ -300,58 +265,58 @@ loop(regi, ); ); !! loop over regi -*** criterion only for checking, not applied anymore: are the anticipation terms sufficienctly small? +*' criterion: are the anticipation terms sufficienctly small? (only for checking, not enforced anymore) p80_fadeoutPriceAnticip_iter(iteration) = sm_fadeoutPriceAnticip; -if(sm_fadeoutPriceAnticip gt cm_maxFadeOutPriceAnticip, -*** s80_bool = 0; !! not an active convergence criterion anymore +if(sm_fadeoutPriceAnticip > cm_maxFadeOutPriceAnticip, + !! s80_bool = 0; !! not an active convergence criterion anymore p80_messageShow("anticip") = YES; ); *' criterion "Deviation due to price anticipation": are the resulting deviations sufficiently small? -*' compare to 1/10th of the cutoff for goods imbalance -if(p80_DevPriceAnticipGlobAllMax2100Iter(iteration) gt 0.1 * p80_surplusMaxTolerance("good"), - s80_bool=0; +*** compare to 1/10th of the cutoff for goods imbalance +if(p80_DevPriceAnticipGlobAllMax2100Iter(iteration) > 0.1 * p80_surplusMaxTolerance("good"), + s80_bool = 0; p80_messageShow("DevPriceAnticip") = YES; ); -*' criterion "Did REMIND run sufficient iterations (currently set at 18, to allow for at least 4 iterations with EDGE-T) -if( (iteration.val le 17), - s80_bool=0; +*' criterion "Were there enough iterations": currently set to allow for at least 4 iterations with EDGE-T +if(iteration.val < cm_startIter_EDGET + 4, + s80_bool = 0; p80_messageShow("IterationNumber") = YES; ); -***additional criterion: did taxes converge? (only checked if cm_TaxConvCheck is 1) +*' criterion: did taxes converge? (only checked if cm_TaxConvCheck is 1) p80_convNashTaxrev_iter(iteration,t,regi) = 0; loop(regi, - loop(t, - p80_convNashTaxrev_iter(iteration,t,regi) = vm_taxrev.l(t,regi) / vm_cesIO.l(t,regi,"inco"); - if (cm_TaxConvCheck eq 1, - if( abs(p80_convNashTaxrev_iter(iteration,t,regi)) gt 0.001, - s80_bool = 0; - p80_messageShow("taxconv") = YES; - ); - ); + loop(t, + p80_convNashTaxrev_iter(iteration,t,regi) = vm_taxrev.l(t,regi) / vm_cesIO.l(t,regi,"inco"); + if (cm_TaxConvCheck = 1, + if(abs(p80_convNashTaxrev_iter(iteration,t,regi)) > 0.001, + s80_bool = 0; + p80_messageShow("taxconv") = YES; + ); ); + ); ); -*** additional criterion: Were regional climate targets reached? +*' criterion: Were regional climate targets reached? $ifthen.emiMkt not "%cm_emiMktTarget%" == "off" -loop((ttot,ttot2,ext_regi,emiMktExt)$pm_emiMktTarget_dev(ttot,ttot2,ext_regi,emiMktExt), - if(NOT(pm_allTargetsConverged(ext_regi) eq 1), +loop((ttot,ttot2,ext_regi,emiMktExt) $ pm_emiMktTarget_dev(ttot,ttot2,ext_regi,emiMktExt), + if(not(pm_allTargetsConverged(ext_regi) = 1), s80_bool = 0; p80_messageShow("regiTarget") = YES; ); ); $endif.emiMkt -*** additional criterion: Were the quantity targets reached by implicit taxes and/or subsidies? +*' criterion: Were the quantity targets reached by implicit taxes and/or subsidies? $ifthen.cm_implicitQttyTarget not "%cm_implicitQttyTarget%" == "off" p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) = pm_implicitQttyTarget_dev(ttot,ext_regi,qttyTarget,qttyTargetGroup); -loop((ttot,ext_regi,taxType,targetType,qttyTarget,qttyTargetGroup)$pm_implicitQttyTarget(ttot,ext_regi,taxType,targetType,qttyTarget,qttyTargetGroup), - if(abs(p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup)) gt cm_implicitQttyTarget_tolerance, - if(NOT ((sameas(taxType,"tax") and p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) lt 0) - OR (sameas(taxType,"sub") and p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) gt 0)), - if(NOT(pm_implicitQttyTarget_isLimited(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) eq 1), !!no tax update either by reaching target or due to tax changes not affecting quantitties +loop((ttot,ext_regi,taxType,targetType,qttyTarget,qttyTargetGroup) $ pm_implicitQttyTarget(ttot,ext_regi,taxType,targetType,qttyTarget,qttyTargetGroup), + if(abs(p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup)) > cm_implicitQttyTarget_tolerance, + if(not ((sameas(taxType,"tax") and p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) < 0) + or (sameas(taxType,"sub") and p80_implicitQttyTarget_dev_iter(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) > 0)), + if(not(pm_implicitQttyTarget_isLimited(iteration,ttot,ext_regi,qttyTarget,qttyTargetGroup) = 1), !! no tax update either by reaching target or due to tax changes not affecting quantities s80_bool = 0; p80_messageShow("implicitEnergyTarget") = YES; ); @@ -360,9 +325,9 @@ loop((ttot,ext_regi,taxType,targetType,qttyTarget,qttyTargetGroup)$pm_implicitQt ); $endif.cm_implicitQttyTarget -*** additional criterion: Were FE price targets reached by implicit taxes and/or subsidies? +*' criterion: Were FE price targets reached by implicit taxes and/or subsidies? $ifthen.cm_implicitPriceTarget not "%cm_implicitPriceTarget%" == "off" -loop((t,regi,entyFe,entySe,sector)$pm_implicitPriceTarget(t,regi,entyFe,entySe,sector), +loop((t,regi,entyFe,entySe,sector) $ pm_implicitPriceTarget(t,regi,entyFe,entySe,sector), if((pm_implicitPrice_NotConv(regi,sector,entyFe,entySe,t)), s80_bool = 0; p80_messageShow("cm_implicitPriceTarget") = YES; @@ -370,9 +335,9 @@ loop((t,regi,entyFe,entySe,sector)$pm_implicitPriceTarget(t,regi,entyFe,entySe,s ); $endIf.cm_implicitPriceTarget -*** additional criterion: Were PE price targets reached by implicit taxes and/or subsidies? +*' criterion: Were PE price targets reached by implicit taxes and/or subsidies? $ifthen.cm_implicitPePriceTarget not "%cm_implicitPePriceTarget%" == "off" -loop((t,regi,entyPe)$pm_implicitPePriceTarget(t,regi,entyPe), +loop((t,regi,entyPe) $ pm_implicitPePriceTarget(t,regi,entyPe), if((pm_implicitPePrice_NotConv(regi,entyPe,t)), s80_bool = 0; p80_messageShow("cm_implicitPePriceTarget") = YES; @@ -380,41 +345,44 @@ loop((t,regi,entyPe)$pm_implicitPePriceTarget(t,regi,entyPe), ); $endIf.cm_implicitPePriceTarget -*** check global budget target from core/postsolve, must be within cm_budgetCO2_absDevTol (default 2 Gt) of target value +*' criterion: global budget target from core/postsolve must be within cm_budgetCO2_absDevTol (default 2Gt) of target value p80_globalBudget_absDev_iter(iteration) = sm_globalBudget_absDev; -if ( abs(p80_globalBudget_absDev_iter(iteration)) gt cm_budgetCO2_absDevTol , !! check if CO2 budget met in tolerance range, +if (abs(p80_globalBudget_absDev_iter(iteration)) > cm_budgetCO2_absDevTol, !! check if CO2 budget met in tolerance range, s80_bool = 0; p80_messageShow("globalbudget") = YES; ); - +*' criterion: is the peak budget year the same as the year of maximum cumulative CO2 emissions? $ifthen.carbonprice %carbonprice% == "functionalForm" *** check whether cm_peakBudgYr corresponds to year of maximum cumulative CO2 emissions -if ( ( cm_iterative_target_adj eq 9 - AND cm_peakBudgYr ne sm_peakBudgYr_check ), +if (cm_iterative_target_adj eq 9 and cm_peakBudgYr ne sm_peakBudgYr_check, s80_bool = 0; p80_messageShow("peakbudgyr") = YES; ); -*** Check whether difference in cumulative emissions between both time steps is greater than sm_peakbudget_diff_tolerance -if ( ( cm_iterative_target_adj eq 9 - AND abs(sm_peakbudget_diff) gt sm_peakbudget_diff_tolerance), +*** check whether difference in cumulative emissions between both time steps is within sm_peakbudget_diff_tolerance +if (cm_iterative_target_adj eq 9 and abs(sm_peakbudget_diff) > sm_peakbudget_diff_tolerance, s80_bool = 0; p80_messageShow("peakbudget") = YES; ); $endIf.carbonprice -*** additional criterion: if damage internalization is on, is damage iteration converged? +*' criterion: if damage internalization is on, is damage iteration converged? p80_sccConvergenceMaxDeviation_iter(iteration) = pm_sccConvergenceMaxDeviation; p80_gmt_conv_iter(iteration) = pm_gmt_conv; $ifthen.internalizeDamages not "%internalizeDamages%" == "off" - if(p80_sccConvergenceMaxDeviation_iter(iteration) gt cm_sccConvergence OR p80_gmt_conv_iter(iteration) gt cm_tempConvergence, - s80_bool = 0; - p80_messageShow("damage") = YES; - ); + if(p80_sccConvergenceMaxDeviation_iter(iteration) > cm_sccConvergence or p80_gmt_conv_iter(iteration) > cm_tempConvergence, + s80_bool = 0; + p80_messageShow("damage") = YES; + ); $endIf.internalizeDamages + +***------------------------------------------------------------------------------ +*' #### Display convergence diagnostics in full.lst +***------------------------------------------------------------------------------ + display "####"; display "Convergence diagnostics"; display "Iteration number: "; @@ -424,293 +392,277 @@ option decimals = 3; display "In the following you find some diagnostics on whether the model converged in this iteration: "; -display "solvestat and modelstat parameters: "; +display "Solvestat and modelstat parameters: "; display p80_repy; -display "trade convergence indicators"; +display "Trade convergence indicators:"; display p80_surplusMaxTolerance, p80_surplusMax2100; display p80_defic_trade, p80_defic_sum,p80_defic_sum_rel; -display "Reasons for non-convergence in this iteration (if not yet converged)"; - - loop(convMessage80$(p80_messageShow(convMessage80)), - if(sameas(convMessage80, "infes"), - display "#### 1.) Infeasibilities found in at least some regions in the last iteration. Please check parameter p80_repy for details. "; - display "#### Try a different gdx, or re-run the optimization with cm_nash_mode set to debug in order to debug the infes."; - ); - if(sameas(convMessage80, "surplus"), - display "#### 2.) Some markets failed to reach a residual surplus below the prescribed threshold. "; - display "#### In the following, the offending markets are indicated by a 1:"; - OPTION decimals = 0; - display p80_messageFailedMarket; - OPTION decimals = 3; - display "#### You will find detailed trade convergence indicators below, search for p80_defic_trade"; - ); - if(sameas(convMessage80, "nonopt"), - display "#### 3.) Found a feasible, but non-optimal solution. This is the infamous status-7 problem: "; - display "#### We can't accept this solution, because it is non-optimal, and, in addition, too far away from the last known optimal solution. "; - display "#### Just trying a different gdx may help."; - ); - if(sameas(convMessage80, "taxconv"), - display "#### 4.) Taxes did not converge in all regions and time steps. Absolute level of tax revenue must be smaller than 0.1 percent of GDP. Check p80_convNashTaxrev_iter below."; - ); - if(sameas(convMessage80, "DevPriceAnticip"), - display "#### 5.) The total monetary value of the price anticipation term times the traded amount are larger than the goods imbalance threshold * 0.1"; - display "#### Check out p80_DevPriceAnticipGlobAllMax2100Iter, which needs to be below 0.1 * the threshold for goods imbalance, p80_surplusMaxTolerance"; - ); - if(sameas(convMessage80, "anticip"), - display "#### 5b.) only for checking, not anymore a criterion that stops convergence: The fadeout price anticipation terms are not sufficiently small."; - display "#### Check out sm_fadeoutPriceAnticip which needs to be below cm_maxFadeOutPriceAnticip."; - display sm_fadeoutPriceAnticip, cm_maxFadeOutPriceAnticip; - ); - if(sameas(convMessage80, "globalbudget"), - display "#### 6.) A global climate target has not been reached yet."; - display "#### check sm_globalBudget_absDev for the deviation from the global target CO2 budget (convergence criterion defined via cm_budgetCO2_absDevTol [default = 2 Gt CO2]), as well as"; - display "#### pm_taxCO2eq_iter (regional CO2 tax path tracked over iterations [T$/GtC]) and"; - display "#### pm_taxCO2eq_anchor_iterationdiff (difference in global anchor carbon price to the last iteration [T$/GtC]) in diagnostics section below."; - display sm_globalBudget_absDev; - ); -$ifthen.carbonprice %carbonprice% == "functionalForm" - if(sameas(convMessage80, "peakbudgyr"), - display "#### 6.) Years are different: cm_peakBudgYr is not equal to sm_peakBudgYr_check."; - display cm_peakBudgYr; - ); - if(sameas(convMessage80, "peakbudget"), - display "#### 6.) PeakBudget not reached: sm_peakbudget_diff is greater than sm_peakbudget_diff_tolerance."; - display sm_peakbudget_diff; - ); -$endIf.carbonprice - if(sameas(convMessage80, "IterationNumber"), - display "#### 0.) REMIND did not run sufficient iterations (currently set at 18, to allow for at least 4 iterations with EDGE-T)"; - ); +display "Reasons for non-convergence in this iteration (if not yet converged):"; +loop(convMessage80$(p80_messageShow(convMessage80)), + if(sameas(convMessage80, "infes"), + display "#### 1.) Infeasibilities found in at least some regions in the last iteration. Please check parameter p80_repy for details. "; + display "#### Try a different gdx, or re-run the optimization with cm_nash_mode set to debug in order to debug the infes."; + ); + if(sameas(convMessage80, "surplus"), + display "#### 2.) Some markets failed to reach a residual surplus below the prescribed threshold. "; + display "#### In the following, the offending markets are indicated by a 1:"; + option decimals = 0; + display p80_messageFailedMarket; + option decimals = 3; + display "#### You will find detailed trade convergence indicators below, search for p80_defic_trade"; + ); + if(sameas(convMessage80, "nonopt"), + display "#### 3.) Found a feasible, but non-optimal solution. This is the infamous status-7 problem: "; + display "#### We can't accept this solution, because it is non-optimal, and, in addition, too far away from the last known optimal solution. "; + display "#### Just trying a different gdx may help."; + ); + if(sameas(convMessage80, "taxconv"), + display "#### 4.) Taxes did not converge in all regions and time steps. Absolute level of tax revenue must be smaller than 0.1 percent of GDP. Check p80_convNashTaxrev_iter below."; + ); + if(sameas(convMessage80, "DevPriceAnticip"), + display "#### 5.) The total monetary value of the price anticipation term times the traded amount are larger than the goods imbalance threshold * 0.1"; + display "#### Check out p80_DevPriceAnticipGlobAllMax2100Iter, which needs to be below 0.1 * the threshold for goods imbalance, p80_surplusMaxTolerance"; + ); + if(sameas(convMessage80, "anticip"), + display "#### 5b.) only for checking, not anymore a criterion that stops convergence: The fadeout price anticipation terms are not sufficiently small."; + display "#### Check out sm_fadeoutPriceAnticip which needs to be below cm_maxFadeOutPriceAnticip."; + display sm_fadeoutPriceAnticip, cm_maxFadeOutPriceAnticip; + ); + if(sameas(convMessage80, "globalbudget"), + display "#### 6.) A global climate target has not been reached yet."; + display "#### check sm_globalBudget_absDev for the deviation from the global target CO2 budget (convergence criterion defined via cm_budgetCO2_absDevTol [default = 2 Gt CO2]), as well as"; + display "#### pm_taxCO2eq_iter (regional CO2 tax path tracked over iterations [T$/GtC]) and"; + display "#### pm_taxCO2eq_anchor_iterationdiff (difference in global anchor carbon price to the last iteration [T$/GtC]) in diagnostics section below."; + display sm_globalBudget_absDev; + ); + if(sameas(convMessage80, "IterationNumber"), + display "#### 0.) REMIND did not run sufficient iterations (currently set at 18, to allow for at least 4 iterations with EDGE-T)"; + ); $ifthen.emiMkt not "%cm_emiMktTarget%" == "off" - if(sameas(convMessage80, "regiTarget"), - display "#### 7) A regional climate target has not been reached yet."; - display "#### Check out the pm_emiMktTarget_dev parameter of 47_regipol module."; - display "#### For budget targets, the parameter gives the percentage deviation of current emissions in relation to the target value."; - display "#### For yearly targets, the parameter gives the current emissions minus the target value in relative terms to the 2005 emissions."; - display "#### The deviation must to be less than pm_emiMktTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of 2005 emissions to reach convergence."; - display pm_emiMktTarget_tolerance, pm_emiMktTarget_dev, pm_factorRescaleemiMktCO2Tax, pm_emiMktCurrent, pm_emiMktTarget, pm_emiMktRefYear; - display pm_emiMktTarget_dev_iter; - display pm_taxemiMkt_iteration; - ); + if(sameas(convMessage80, "regiTarget"), + display "#### 7) A regional climate target has not been reached yet."; + display "#### Check out the pm_emiMktTarget_dev parameter of 47_regipol module."; + display "#### For budget targets, the parameter gives the percentage deviation of current emissions in relation to the target value."; + display "#### For yearly targets, the parameter gives the current emissions minus the target value in relative terms to the 2005 emissions."; + display "#### The deviation must to be less than pm_emiMktTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of 2005 emissions to reach convergence."; + display pm_emiMktTarget_tolerance, pm_emiMktTarget_dev, pm_factorRescaleemiMktCO2Tax, pm_emiMktCurrent, pm_emiMktTarget, pm_emiMktRefYear; + display pm_emiMktTarget_dev_iter; + display pm_taxemiMkt_iteration; + ); $endif.emiMkt $ifthen.cm_implicitQttyTarget not "%cm_implicitQttyTarget%" == "off" - if(sameas(convMessage80, "implicitEnergyTarget"), - display "#### 10) A quantity target has not been reached yet."; - display "#### Check out the pm_implicitQttyTarget_dev parameter of 47_regipol module."; - display "#### The relative deviation must to be less than cm_implicitQttyTarget_tolerance, which is 1 percent by default."; - display "#### For taxes, this means every value > +0.01, while for subsidies everything < -0.01 is problematic in the following lines."; - display cm_implicitQttyTarget_tolerance, pm_implicitQttyTarget_dev; - ); + if(sameas(convMessage80, "implicitEnergyTarget"), + display "#### 10) A quantity target has not been reached yet."; + display "#### Check out the pm_implicitQttyTarget_dev parameter of 47_regipol module."; + display "#### The relative deviation must to be less than cm_implicitQttyTarget_tolerance, which is 1 percent by default."; + display "#### For taxes, this means every value > +0.01, while for subsidies everything < -0.01 is problematic in the following lines."; + display cm_implicitQttyTarget_tolerance, pm_implicitQttyTarget_dev; + ); $endif.cm_implicitQttyTarget $ifthen.cm_implicitPriceTarget not "%cm_implicitPriceTarget%" == "off" - if(sameas(convMessage80, "cm_implicitPriceTarget"), - display "#### 11) A price target has not been reached yet."; - display "#### Check out below the pm_implicitPrice_NotConv parameter values for non convergence cases."; - display "#### Deviations must be lower than 5%."; - display "#### The pm_implicitPrice_ignConv stores the cases disconsidered in the convergence check."; - display pm_implicitPrice_NotConv, pm_implicitPrice_ignConv; - ); + if(sameas(convMessage80, "cm_implicitPriceTarget"), + display "#### 11) A price target has not been reached yet."; + display "#### Check out below the pm_implicitPrice_NotConv parameter values for non convergence cases."; + display "#### Deviations must be lower than 5%."; + display "#### The pm_implicitPrice_ignConv stores the cases disconsidered in the convergence check."; + display pm_implicitPrice_NotConv, pm_implicitPrice_ignConv; + ); $endIf.cm_implicitPriceTarget $ifthen.cm_implicitPePriceTarget not "%cm_implicitPePriceTarget%" == "off" - if(sameas(convMessage80, "cm_implicitPePriceTarget"), - display "#### 11) A primary energy price target has not been reached yet."; - display "#### Check out below the pm_implicitPePrice_NotConv parameter values for non convergence cases."; - display "#### Deviations must be lower than 5%."; - display "#### The pm_implicitPePrice_ignConv stores the cases disconsidered in the convergence check."; - display pm_implicitPePrice_NotConv, pm_implicitPePrice_ignConv; - ); + if(sameas(convMessage80, "cm_implicitPePriceTarget"), + display "#### 11) A primary energy price target has not been reached yet."; + display "#### Check out below the pm_implicitPePrice_NotConv parameter values for non convergence cases."; + display "#### Deviations must be lower than 5%."; + display "#### The pm_implicitPePrice_ignConv stores the cases disconsidered in the convergence check."; + display pm_implicitPePrice_NotConv, pm_implicitPePrice_ignConv; + ); $endIf.cm_implicitPePriceTarget $ifthen.internalizeDamages not "%internalizeDamages%" == "off" if(sameas(convMessage80,"damage"), - display "#### 11) The damage iteration did not converge."; - display "#### Check out below the values for pm_gmt_conv and pm_sccConvergenceMaxDeviation." - display "#### They should be below 0.05." - display pm_gmt_conv, pm_sccConvergenceMaxDeviation; + display "#### 11) The damage iteration did not converge."; + display "#### Check out below the values for pm_gmt_conv and pm_sccConvergenceMaxDeviation." + display "#### They should be below 0.05." + display pm_gmt_conv, pm_sccConvergenceMaxDeviation; ); $endIf.internalizeDamages - ); +); display "See the indicators below to dig deeper on the respective reasons of non-convergence: " -display "tax convergence indicators"; +display "Tax convergence indicators:"; display p80_convNashTaxrev_iter; -display "detailed trade convergence indicators"; +display "Detailed trade convergence indicators:"; display p80_defic_trade, p80_defic_sum,p80_defic_sum_rel; -OPTION decimals = 7; +option decimals = 7; display p80_surplus; -OPTION decimals = 3; +option decimals = 3; -display "Carbon tax tracked over iterations of 45_carbonprice/functionalForm/postsolve"; +display "Carbon tax tracked over iterations of 45_carbonprice/functionalForm/postsolve:"; display pm_taxCO2eq_iter; -display "Carbon tax difference to last iteration for global targets of 45_carbonprice/functionalForm/postsolve"; +display "Carbon tax difference to last iteration for global targets of 45_carbonprice/functionalForm/postsolve:"; display pm_taxCO2eq_anchor_iterationdiff; - -*RP* display effect of additional convergence push -display "display effect of additional convergence push"; -display o80_trackSurplusSign, o80_SurplusOverTolerance, o80_counter_iteration_trade_ttot, p80_etaST_correct_safecopy,p80_etaST_correct,p80_pvp_itr; - - - -***end with failure message if max number of iterations is reached w/o convergence: -if( (s80_bool eq 0) and (iteration.val eq cm_iteration_max), !! reached max number of iteration, still no convergence - OPTION decimals = 3; - display "################################################################################################"; - display "#################################### Nash Solution Report ####################################"; - display "################################################################################################"; - display "#### !! Nash did NOT converge within the maximum number of iterations allowed !!" - display "#### The reasons for failing to successfully converge are:" - loop(convMessage80$(p80_messageShow(convMessage80)), - if(sameas(convMessage80, "infes"), - display "####"; - display "#### 1.) Infeasibilities found in at least some regions in the last iteration. Plase check parameter p80_repy for details. "; - display "#### Try a different gdx, or re-run the optimization with cm_nash_mode set to debug in order to debug the infes."; - display p80_repy; - ); - if(sameas(convMessage80 , "surplus"), - display "####"; - display "#### 2.) Some markets failed to reach a residual surplus below the prescribed threshold. "; - display "#### You may try less stringent convergence target (a lower cm_nash_autoconverge), or a different gdx. "; - display "#### In the following, the offending markets are indicated by a 1:"; - OPTION decimals = 0; - display p80_messageFailedMarket; - OPTION decimals = 3; - ); - if(sameas(convMessage80, "nonopt"), - display "####"; - display "#### 3.) Found a feasible, but non-optimal solution. This is the infamous status-7 problem: "; - display "#### We can't accept this solution, because it is non-optimal, and too far away from the last known optimal solution. "; - display "#### Just trying a different gdx may help."; - ); - if(sameas(convMessage80, "taxconv"), - display "####"; - display "#### 4.) Taxes did not converge in all regions and time steps. Absolut level of tax revenue must be smaller than 0.1 percent of GDP. Check p80_convNashTaxrev_iter."; - ); - if(sameas(convMessage80, "anticip"), - display "#### 5.) The fadeout price anticipation terms are not sufficiently small."; - ); - if(sameas(convMessage80, "globalbudget"), - display "#### 6.) A global climate target has not been reached yet."; - display "#### check sm_globalBudget_absDev for the deviation from the global target CO2 budget (convergence criterion defined via cm_budgetCO2_absDevTol [default = 2 Gt CO2]), as well as"; - display "#### pm_taxCO2eq_iter (regional CO2 tax path tracked over iterations [T$/GtC]) and"; - display "#### pm_taxCO2eq_anchor_iterationdiff (difference in global anchor carbon price to the last iteration [T$/GtC]) in diagnostics section below."; - display sm_globalBudget_absDev; - ); +display "Display effect of additional convergence push:"; +display o80_trackSurplusSign, o80_SurplusOverTolerance, o80_counter_iteration_trade_ttot, p80_shortTermCorrect_safecopy,p80_shortTermCorrect,p80_pvp_itr; + + + +*** end with failure message if max number of iterations is reached without convergence +if(s80_bool = 0 and iteration.val = cm_iteration_max, + option decimals = 3; + display "################################################################################################"; + display "#################################### Nash Solution Report ####################################"; + display "################################################################################################"; + display "#### !! Nash did NOT converge within the maximum number of iterations allowed !!" + display "#### The reasons for failing to successfully converge are:" + loop(convMessage80$(p80_messageShow(convMessage80)), + if(sameas(convMessage80, "infes"), + display "####"; + display "#### 1.) Infeasibilities found in at least some regions in the last iteration. Plase check parameter p80_repy for details. "; + display "#### Try a different gdx, or re-run the optimization with cm_nash_mode set to debug in order to debug the infes."; + display p80_repy; + ); + if(sameas(convMessage80 , "surplus"), + display "####"; + display "#### 2.) Some markets failed to reach a residual surplus below the prescribed threshold. "; + display "#### You may try less stringent convergence target (a lower cm_nash_autoconverge), or a different gdx. "; + display "#### In the following, the offending markets are indicated by a 1:"; + option decimals = 0; + display p80_messageFailedMarket; + option decimals = 3; + ); + if(sameas(convMessage80, "nonopt"), + display "####"; + display "#### 3.) Found a feasible, but non-optimal solution. This is the infamous status-7 problem: "; + display "#### We can't accept this solution, because it is non-optimal, and too far away from the last known optimal solution. "; + display "#### Just trying a different gdx may help."; + ); + if(sameas(convMessage80, "taxconv"), + display "####"; + display "#### 4.) Taxes did not converge in all regions and time steps. Absolut level of tax revenue must be smaller than 0.1 percent of GDP. Check p80_convNashTaxrev_iter."; + ); + if(sameas(convMessage80, "anticip"), + display "#### 5.) The fadeout price anticipation terms are not sufficiently small."; + ); + if(sameas(convMessage80, "globalbudget"), + display "#### 6.) A global climate target has not been reached yet."; + display "#### check sm_globalBudget_absDev for the deviation from the global target CO2 budget (convergence criterion defined via cm_budgetCO2_absDevTol [default = 2 Gt CO2]), as well as"; + display "#### pm_taxCO2eq_iter (regional CO2 tax path tracked over iterations [T$/GtC]) and"; + display "#### pm_taxCO2eq_anchor_iterationdiff (difference in global anchor carbon price to the last iteration [T$/GtC]) in diagnostics section below."; + display sm_globalBudget_absDev; + ); $ifthen.carbonprice %carbonprice% == "functionalForm" - if(sameas(convMessage80, "peakbudgyr"), - display "#### 6.) Years are different: cm_peakBudgYr is not equal to sm_peakBudgYr_check."; - display cm_peakBudgYr; - ); - if(sameas(convMessage80, "peakbudget"), - display "#### 6.) PeakBudget not reached: sm_peakbudget_diff is greater than sm_peakbudget_diff_tolerance."; - display sm_peakbudget_diff; - ); + if(sameas(convMessage80, "peakbudgyr"), + display "#### 6.) Years are different: cm_peakBudgYr is not equal to sm_peakBudgYr_check."; + display cm_peakBudgYr; + ); + if(sameas(convMessage80, "peakbudget"), + display "#### 6.) PeakBudget not reached: sm_peakbudget_diff is greater than sm_peakbudget_diff_tolerance."; + display sm_peakbudget_diff; + ); $endIf.carbonprice $ifthen.emiMkt not "%cm_emiMktTarget%" == "off" - if(sameas(convMessage80, "regiTarget"), - display "#### 7) A regional climate target has not been reached yet."; - display "#### Check out the pm_emiMktTarget_dev parameter of 47_regipol module."; - display "#### For budget targets, the parameter gives the percentage deviation of current emissions in relation to the target value."; - display "#### For yearly targets, the parameter gives the current emissions minus the target value in relative terms to the 2005 emissions."; - display "#### The deviation must to be less than pm_emiMktTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of 2005 emissions to reach convergence."; - display pm_emiMktTarget_tolerance, pm_emiMktTarget_dev, pm_factorRescaleemiMktCO2Tax, pm_emiMktCurrent, pm_emiMktTarget, pm_emiMktRefYear; - display pm_emiMktTarget_dev_iter; - display pm_taxemiMkt_iteration; - ); + if(sameas(convMessage80, "regiTarget"), + display "#### 7) A regional climate target has not been reached yet."; + display "#### Check out the pm_emiMktTarget_dev parameter of 47_regipol module."; + display "#### For budget targets, the parameter gives the percentage deviation of current emissions in relation to the target value."; + display "#### For yearly targets, the parameter gives the current emissions minus the target value in relative terms to the 2005 emissions."; + display "#### The deviation must to be less than pm_emiMktTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of 2005 emissions to reach convergence."; + display pm_emiMktTarget_tolerance, pm_emiMktTarget_dev, pm_factorRescaleemiMktCO2Tax, pm_emiMktCurrent, pm_emiMktTarget, pm_emiMktRefYear; + display pm_emiMktTarget_dev_iter; + display pm_taxemiMkt_iteration; + ); $endif.emiMkt $ifthen.cm_implicitQttyTarget not "%cm_implicitQttyTarget%" == "off" - if(sameas(convMessage80, "implicitEnergyTarget"), - display "#### 10) A quantity target has not been reached yet."; - display "#### Check out the pm_implicitQttyTarget_dev parameter of 47_regipol module."; - display "#### The deviation must to be less than cm_implicitQttyTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of the defined target."; - display cm_implicitQttyTarget_tolerance, pm_implicitQttyTarget_dev; - ); + if(sameas(convMessage80, "implicitEnergyTarget"), + display "#### 10) A quantity target has not been reached yet."; + display "#### Check out the pm_implicitQttyTarget_dev parameter of 47_regipol module."; + display "#### The deviation must to be less than cm_implicitQttyTarget_tolerance. By default within 1%, i.e. in between -0.01 and 0.01 of the defined target."; + display cm_implicitQttyTarget_tolerance, pm_implicitQttyTarget_dev; + ); $endif.cm_implicitQttyTarget $ifthen.cm_implicitPriceTarget not "%cm_implicitPriceTarget%" == "off" - if(sameas(convMessage80, "cm_implicitPriceTarget"), - display "#### 11) A final energy price target has not been reached yet."; - display "#### Check out below the pm_implicitPrice_NotConv parameter values for non convergence cases."; - display "#### Deviations must be lower than 5%."; - display "#### The pm_implicitPrice_ignConv stores the cases disconsidered in the convergence check."; - display pm_implicitPrice_NotConv, pm_implicitPrice_ignConv; - ); + if(sameas(convMessage80, "cm_implicitPriceTarget"), + display "#### 11) A final energy price target has not been reached yet."; + display "#### Check out below the pm_implicitPrice_NotConv parameter values for non convergence cases."; + display "#### Deviations must be lower than 5%."; + display "#### The pm_implicitPrice_ignConv stores the cases disconsidered in the convergence check."; + display pm_implicitPrice_NotConv, pm_implicitPrice_ignConv; + ); $endIf.cm_implicitPriceTarget $ifthen.cm_implicitPePriceTarget not "%cm_implicitPePriceTarget%" == "off" - if(sameas(convMessage80, "cm_implicitPePriceTarget"), - display "#### 11) A primary energy price target has not been reached yet."; - display "#### Check out below the pm_implicitPePrice_NotConv parameter values for non convergence cases."; - display "#### Deviations must be lower than 5%."; - display "#### The pm_implicitPePrice_ignConv stores the cases disconsidered in the convergence check."; - display pm_implicitPePrice_NotConv, pm_implicitPePrice_ignConv; - ); + if(sameas(convMessage80, "cm_implicitPePriceTarget"), + display "#### 11) A primary energy price target has not been reached yet."; + display "#### Check out below the pm_implicitPePrice_NotConv parameter values for non convergence cases."; + display "#### Deviations must be lower than 5%."; + display "#### The pm_implicitPePrice_ignConv stores the cases disconsidered in the convergence check."; + display pm_implicitPePrice_NotConv, pm_implicitPePrice_ignConv; + ); $endIf.cm_implicitPePriceTarget - ); - display "#### Info: These residual market surplusses in current monetary values are:"; - display p80_defic_trade; - display "#### The sum of those, normalized to the total consumption, given in percent is: "; - display p80_defic_sum_rel; + ); - display "################################################################################################"; - display "################################################################################################"; + display "#### Info: These residual market surplusses in current monetary values are:"; + display p80_defic_trade; + display "#### The sum of those, normalized to the total consumption, given in percent is: "; + display p80_defic_sum_rel; + display "################################################################################################"; + display "################################################################################################"; ); -***if all conditions are met, stop optimization. -if(s80_bool eq 1, -***in automatic mode, set iteration_max such that no next iteration takes place - if(cm_nash_autoconverge ne 0, - cm_iteration_max = iteration.val - 1; - ); - OPTION decimals = 3; - s80_numberIterations = cm_iteration_max + 1; - display "######################################################################################################"; - display "Run converged!!"; - display "#### Nash Solution Report"; - display "#### Convergence threshold reached within ",s80_numberIterations, "iterations."; - display "############"; - display "Model solution parameters of last iteration"; - display p80_repy; - display "#### Residual market surpluses in 2100 are:"; - display p80_surplusMax2100; - display "#### This meets the prescribed tolerance requirements of: "; - display p80_surplusMaxTolerance; - display "#### Info: These residual market surplusses in monetary are :"; - display p80_defic_trade; - display "#### Info: And the sum of those (equivalent to Negishi's defic_sum):"; - display p80_defic_sum; - display "#### This value in percent of the NPV of consumption is: "; - display p80_defic_sum_rel; - display "############"; - display "######################################################################################################"; - OPTION decimals = 3; - s80_converged = 1; !! set machine-readable status parameter - +***------------------------------------------------------------------------------ +*' #### Finishing or aborting +***------------------------------------------------------------------------------ + +*** if all conditions are met, stop optimization. +if(s80_bool = 1, + s80_converged = 1; !! set machine-readable status parameter + s80_numberIterations = iteration.val; + cm_iteration_max $ (cm_nash_autoconverge > 0) = iteration.val - 1; !! set iteration_max such that no next iteration takes place + + option decimals = 3; + display "######################################################################################################"; + display "Run converged!!"; + display "#### Nash Solution Report"; + display "#### Convergence threshold reached within ", s80_numberIterations, "iterations."; + display "############"; + display "Model solution parameters of last iteration"; + display p80_repy; + display "#### Residual market surpluses in 2100 are:"; + display p80_surplusMax2100; + display "#### This meets the prescribed tolerance requirements of: "; + display p80_surplusMaxTolerance; + display "#### Info: These residual market surplusses in monetary are :"; + display p80_defic_trade; + display "#### Info: And the sum of those (equivalent to Negishi's defic_sum):"; + display p80_defic_sum; + display "#### This value in percent of the NPV of consumption is: "; + display p80_defic_sum_rel; + display "############"; + display "######################################################################################################"; + option decimals = 3; ); -*** check if any region has failed to solve consecutively for -*** cm_abortOnConsecFail times -if (cm_abortOnConsecFail gt 0, +*** check if any region has failed to solve consecutively for cm_abortOnConsecFail times +if (cm_abortOnConsecFail > 0, loop (regi, - if ( ( p80_repy_iteration(regi,"solvestat",iteration) eq 1 - AND p80_repy_iteration(regi,"modelstat",iteration) eq 2) - OR ( p80_repy_iteration(regi,"solvestat",iteration) eq 4 - AND p80_repy_iteration(regi,"modelstat",iteration) eq 7), + if ( (p80_repy_iteration(regi,"solvestat",iteration) = 1 and p80_repy_iteration(regi,"modelstat",iteration) = 2) + or (p80_repy_iteration(regi,"solvestat",iteration) = 4 and p80_repy_iteration(regi,"modelstat",iteration) = 7), !! region was solved successfully p80_trackConsecFail(regi) = 0; - else - !! region failed to solve + else !! region failed to solve p80_trackConsecFail(regi) = p80_trackConsecFail(regi) + 1; ); ); if (smax(regi, p80_trackConsecFail(regi)) >= cm_abortOnConsecFail, - if ((s80_runInDebug eq 0) AND (cm_nash_mode ne 1), !! auto-start debug only if not already in debug mode - if (sum(regi, pm_SolNonInfes(regi) ne 0) eq 0, !! if all regions are infeasible debug makes no sense + if (s80_runInDebug = 0 and cm_nash_mode ne 1, !! auto-start debug only if not already in debug mode + if (sum(regi, pm_SolNonInfes(regi) ne 0) = 0, !! if all regions are infeasible debug makes no sense execute_unload "abort.gdx"; abort "Run was aborted because the maximum number of consecutive failures was reached in at least one region! No debug started since all regions are infeasible."; else !! start debug mode only if at leat one region was feasible @@ -718,13 +670,12 @@ if (cm_abortOnConsecFail gt 0, cm_nash_mode = 1; display "Starting nash in debug mode after maximum number of consecutive failures was reached in at least one region."; ); - else !! s80_runInDebug eq 1 AND/OR cm_nash_mode eq 1 + else !! s80_runInDebug = 1 and/or cm_nash_mode = 1 execute_unload "abort.gdx"; abort "After debug mode run was aborted because the maximum number of consecutive failures was still reached in at least one region!"; ); - else - !! Set nash mode back to parallel because all regions got feasible after they have been automatically restarted as debug - if (s80_runInDebug eq 1, + else !! set nash mode back to parallel because all regions got feasible after they have been automatically restarted as debug + if (s80_runInDebug = 1, s80_runInDebug = 0; cm_nash_mode = 2; display "Set nash mode back to parallel after regions got feasible in auto-debug mode."; @@ -733,18 +684,16 @@ if (cm_abortOnConsecFail gt 0, ); -***-------------------------- -*** EMIOPT implementation -***-------------------------- +***------------------------------------------------------------------------------ +*' #### EMIOPT nash algorithm implementation +***------------------------------------------------------------------------------ $ifthen.emiopt %emicapregi% == 'none' -if(cm_emiscen eq 6, - -*** nash emiopt algorithm +if(cm_emiscen = 6, !! budget *** we iteratively reach the point where these two marginals are equal for each region by adjusting regional permit budgets: *** marginal of cumulative emissions: -p80_eoMargEmiCum(regi) = 5*(abs(qm_co2eqCum.m(regi)))$(pm_SolNonInfes(regi) eq 1); +p80_eoMargEmiCum(regi) = 5*(abs(qm_co2eqCum.m(regi))) $ (pm_SolNonInfes(regi) = 1); *** marginal of permit budget : -p80_eoMargPermBudg(regi) = 5*(abs(q80_budgetPermRestr.m(regi)))$(pm_SolNonInfes(regi) eq 1); +p80_eoMargPermBudg(regi) = 5*(abs(q80_budgetPermRestr.m(regi))) $ (pm_SolNonInfes(regi) = 1); display pm_budgetCO2eq; @@ -763,7 +712,7 @@ p80_count=0; p80_count = smax(regi, p80_eoEmiMarg(regi)); loop(regi, *** dealing with infeasibles - if ((pm_SolNonInfes(regi) eq 0), + if ((pm_SolNonInfes(regi) = 0), p80_eoEmiMarg(regi) = p80_count; else p80_eoEmiMarg(regi) = p80_eoEmiMarg(regi); ); @@ -789,14 +738,8 @@ p80_eoMargDiffItr(regi,iteration) = p80_eoMargDiff(regi); p80_eoEmibudgetDiffAbs(iteration) = sum(regi, abs(p80_eoMargDiff(regi) * p80_eoDeltaEmibudget) ); option decimals = 5; -display p80_eoMargEmiCum, p80_eoMargPermBudg, p80_eoEmiMarg, p80_eoMargAverage, p80_eoMargDiff, p80_eoDeltaEmibudget, p80_eoWeights,p80_eoEmibudget1RegItr -; - +display p80_eoMargEmiCum, p80_eoMargPermBudg, p80_eoEmiMarg, p80_eoMargAverage, p80_eoMargDiff, p80_eoDeltaEmibudget, p80_eoWeights,p80_eoEmibudget1RegItr; ); $endif.emiopt - - - - *** EOF ./modules/80_optimization/nash/postsolve.gms diff --git a/modules/80_optimization/nash/preloop.gms b/modules/80_optimization/nash/preloop.gms index 9cf0c48be1..6a665a8211 100644 --- a/modules/80_optimization/nash/preloop.gms +++ b/modules/80_optimization/nash/preloop.gms @@ -5,97 +5,77 @@ *** | REMIND License Exception, version 1.0 (see LICENSE file). *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/preloop.gms -*MLB/AG* for Nash algorithm read initial price data from gdx - Execute_Loadpoint 'input' pm_pvp = pm_pvp; - Execute_Loadpoint 'input' vm_Xport.l = vm_Xport.l; - Execute_Loadpoint 'input' vm_Mport.l = vm_Mport.l; - Execute_Loadpoint 'input' vm_cons.l = vm_cons.l; - Execute_Loadpoint 'input' vm_taxrev.l = vm_taxrev.l; - Execute_Loadpoint 'input' vm_fuExtr.l = vm_fuExtr.l; - Execute_Loadpoint 'input' vm_prodPe.l = vm_prodPe.l; -*** assign fake values for p80_repyLastOptim which gets initialised in the loop -p80_repyLastOptim(regi,solveinfo80) = NA; +***------------------------------------------------------------------------------ +*** Use gdx to initialise prices, trade volumes etc +***------------------------------------------------------------------------------ +*** read data from gdx with explicit instruction execute_loadpoint +Execute_Loadpoint 'input' pm_pvp = pm_pvp; +Execute_Loadpoint 'input' vm_Xport.l = vm_Xport.l; +Execute_Loadpoint 'input' vm_Mport.l = vm_Mport.l; +Execute_Loadpoint 'input' vm_cons.l = vm_cons.l; +Execute_Loadpoint 'input' vm_taxrev.l = vm_taxrev.l; +Execute_Loadpoint 'input' vm_fuExtr.l = vm_fuExtr.l; +Execute_Loadpoint 'input' vm_prodPe.l = vm_prodPe.l; -*AJS* initialize starting points for prices, trade volumes etc. from gdx. -***in order to read parameters like p80_priceXXX from a gdx, instead of only variables , we have to explicitly instruct gams to do so in the execute_loadpoint command in core/preloop.gms. -***the price paths are the trickiest part. Try to find p80_priceXXX prices in the gdx fist, if that fails, fallback to price path read from input/prices_NASH.inc -loop(ttot$(ttot.val ge 2005), - loop(trade$(NOT tradeSe(trade)), - if((pm_pvp(ttot,trade) eq NA) OR (pm_pvp(ttot,trade) lt 1E-12) OR (pm_pvp(ttot,trade) gt 0.1) , -***in case we have not been able to read price paths from the gdx, or these price are zero (or eps), fall back to price paths in file input/prices_NASH.inc: - pm_pvp(ttot,trade) = p80_pvpFallback(ttot,trade); - display 'Nash: Info: Could not load useful initial price from gdx, falling back to the one found in input/prices_NASH.inc. This should not be a problem, the runs can stil converge. '; - ); -***in case pm_pvp is not found in gdx: - if(pm_pvp(ttot,trade) eq NA, - pm_pvp(ttot,trade) = 0; - ); - loop(regi, - pm_Xport0(ttot,regi,trade) = vm_Xport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade) = vm_Mport.l(ttot,regi,trade); - -***in case xport/mport is not found in gdx: - if(pm_Xport0(ttot,regi,trade) eq NA, - pm_Xport0(ttot,regi,trade) = 0; - vm_Xport.L(ttot,regi,trade) = 0; - ); - if(p80_Mport0(ttot,regi,trade) eq NA, - p80_Mport0(ttot,regi,trade) = 0; - vm_Mport.L(ttot,regi,trade) = 0; - ); - p80_normalize0(ttot,regi,"good") = vm_cons.l(ttot,regi); -*** p80_normalize0(ttot,regi,"perm") = vm_cons.l(ttot,regi); - p80_normalize0(ttot,regi,"perm")$(ttot.val ge 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1E-6); - p80_normalize0(ttot,regi,tradePe) = 0.5 * (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)); +loop(ttot $ (ttot.val >= 2005), + loop(trade $ (not tradeSe(trade)), +*** in case price paths from the gdx are not valid, fall back to prices read from input/prices_NASH.inc + if(pm_pvp(ttot,trade) = NA or pm_pvp(ttot,trade) < 1e-12 or pm_pvp(ttot,trade) > 0.1, + pm_pvp(ttot,trade) = p80_pvpFallback(ttot,trade); + display 'Nash: Info: Could not load useful initial price from gdx, falling back to the one found in input/prices_NASH.inc. This should not be a problem, the runs can stil converge. '; + ); -p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); +*** initialise trade volumes from gdx, and set missing values to zero + pm_Xport0(ttot,regi,trade) = vm_Xport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) = vm_Mport.l(ttot,regi,trade); + pm_Xport0(ttot,regi,trade) $ (pm_Xport0(ttot,regi,trade) = NA) = 0; + vm_Xport.l(ttot,regi,trade) $ (pm_Xport0(ttot,regi,trade) = NA) = 0; + p80_Mport0(ttot,regi,trade) $ (p80_Mport0(ttot,regi,trade) = NA) = 0; + vm_Mport.l(ttot,regi,trade) $ (p80_Mport0(ttot,regi,trade) = NA) = 0; - ); - ); -); +*** initialise market volume for different trades; has to match calculation in nash/postsolve + p80_regiMarketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_regiMarketVolume(ttot,regi,"perm") = abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)); + p80_regiMarketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; + p80_regiMarketVolume(ttot,regi,trade) = max(sm_eps, p80_regiMarketVolume(ttot,regi,trade)); !! ensure market volume is positive -loop(regi, - loop(tradePe, -if(p80_Mport0("2005",regi,tradePe) eq NA, p80_Mport0("2005",regi,tradePe) = 0); -);); - -*AJS* starting policy runs from permit prices that are all zero doesnot work. start from 30$ price path instead -if((cm_emiscen ne 1) and (cm_emiscen ne 9) and (smax(t,pm_pvp(t,"perm"))) eq 0, - loop(ttot$(ttot.val ge 2005), -***this is a 30$/tCo2eq in 2020 trajectory: - pm_pvp(ttot,"perm") = 0.11*1.05**(ttot.val-2020) * pm_pvp(ttot,"good"); - ); - pm_pvp("2005","perm")=0; + p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); + ); ); -if((cm_emiscen eq 1) or (cm_emiscen eq 9), !! if there is no period trade, set the price to zero. - pm_pvp(ttot,"perm")=0; +*** AJS: starting policy runs from permit prices that are all zero does not work; start from 30$ price path instead +if(cm_emiscen ne 1 and cm_emiscen ne 9 and smax(t, pm_pvp(t,"perm")) = 0, + pm_pvp("2005","perm") = 0; + loop(ttot $ (ttot.val > 2005), + pm_pvp(ttot,"perm") = 0.11 * 1.05**(ttot.val - 2020) * pm_pvp(ttot,"good"); !! this is a 30$/tCo2eq in 2020 trajectory + ); ); +*** if there is no permit trade, set the price to zero +pm_pvp(ttot,"perm") $ (cm_emiscen = 1 or cm_emiscen = 9) = 0; - -p80_pvp_itr(ttot,trade,"1")$(NOT tradeSe(trade)) = pm_pvp(ttot,trade); - -*AJS* Take care of resource prices that were imported as zero (as seen for 2150, peur), as they cause problems in the covergence process. Default to last periods price: -loop(tradePe, - loop(ttot$(NOT sameas(ttot,'2005')), - if(p80_pvp_itr(ttot,tradePe,"1") eq 0, - p80_pvp_itr(ttot,tradePe,"1") = p80_pvp_itr(ttot-1,tradePe,"1")$(NOT sameas(ttot,'2005')); - ); - ); +*** AJS: for resource with price zero, fall back to previous period in order to avoid convergence problems (seen for peur in 2150) +loop(ttot $ (ttot.val > 2005), + pm_pvp(ttot,tradePe) $ (pm_pvp(ttot,tradePe) = 0) = pm_pvp(ttot-1,tradePe); ); -***debug display -display pm_pvp,p80_normalize0; +*** save prices of the first iteration +p80_pvp_itr(ttot,trade,"1") $ (not tradeSe(trade)) = pm_pvp(ttot,trade); + +*** debug display +display pm_pvp,p80_regiMarketVolume; display pm_Xport0,p80_Mport0; display p80_surplusMaxTolerance; -*EMIOPT +*** assign fake values for p80_repyLastOptim which gets initialised in nash/solve +p80_repyLastOptim(regi,solveinfo80) = NA; + +*** EMIOPT $ifthen.emiopt %emicapregi% == 'none' -if(cm_emiscen eq 6, -pm_budgetCO2eq(regi) = pm_shPerm("2050",regi) * sm_budgetCO2eqGlob; -display pm_shPerm, sm_budgetCO2eqGlob, pm_budgetCO2eq; +if(cm_emiscen = 6, + pm_budgetCO2eq(regi) = pm_shPerm("2050",regi) * sm_budgetCO2eqGlob; + display pm_shPerm, sm_budgetCO2eqGlob, pm_budgetCO2eq; ); $endif.emiopt diff --git a/modules/80_optimization/nash/presolve.gms b/modules/80_optimization/nash/presolve.gms index 42598ca19a..cb427b69d3 100644 --- a/modules/80_optimization/nash/presolve.gms +++ b/modules/80_optimization/nash/presolve.gms @@ -5,20 +5,24 @@ *** | REMIND License Exception, version 1.0 (see LICENSE file). *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/presolve.gms -*MLB* 12062013* update of learning externality -*MLB 20130920* update of spillover eternality -*ML+LB 20140107* update of climate externality (carbon market) -*LB* 20140506 update of climate externality (link to the climate module) -*AJS* TODO: this is non-optimal: in summing over regi2, we'd like to keep the last known optimal value. To achieve that, param0 should not be overwritten in core/presolve.gms - we'd have to make p80_repynoninfes an interface. -pm_capCumForeign(ttot,regi,teLearn)$((ttot.val ge 2005) and (pm_SolNonInfes(regi) eq 1) ) = sum(regi2$((NOT sameas(regi,regi2))), pm_capCum0(ttot,regi2,teLearn)); -pm_capCumForeign(ttot,regi,teLearn)$((ttot.val ge 2025) and (pm_SolNonInfes(regi) eq 1) and (cm_LearningSpillover eq 0)) = sum(regi2$((NOT sameas(regi,regi2))), pm_capCum0("2020",regi2,teLearn)); -* if cm_LearningSpillover is 0, the foreign capacity in technology learning is fixed to the level of 2020. This simulates a world of protectionism with noch no further foreign technology learning spillover. -pm_cumEff(ttot,regi, in)$(ttot.val ge 2005 and pm_SolNonInfes(regi) eq 1) = sum( regi2$(pm_SolNonInfes(regi2) eq 1), (pm_cesdata("2005",regi2,in,"eff") * vm_effGr.l(ttot,regi2,in))) - (pm_cesdata("2005",regi,in,"eff") * vm_effGr.l(ttot,regi,in)); !! TODO: take care of the case of infeasible solution +*** Compute sum of foreign variables needed for regional optimisation (use the last known optimal value) +loop(regi $ (pm_SolNonInfes(regi) = 1), + loop(ttot $ (ttot.val >= 2005), + pm_cumEff(ttot,regi,in) + = sum(regi2 $ (not sameas(regi,regi2) and pm_SolNonInfes(regi2) = 1), + pm_cesdata("2005",regi2,in,"eff") * vm_effGr.l(ttot,regi2,in)); !! TODO: case where regi2 is infeasible -pm_co2eqForeign(ttot,regi)$((ttot.val ge 2005) and (pm_SolNonInfes(regi) eq 1)) = sum(regi2$((NOT sameas(regi,regi2))), pm_co2eq0(ttot,regi2)); !! does this interfere with the initialization in datainput? - -pm_fuExtrForeign(ttot,regi,enty,rlf)$((ttot.val ge 2005) and (pm_SolNonInfes(regi) eq 1)) = sum(regi2$((NOT sameas(regi,regi2))), vm_fuExtr.l(ttot,regi2,enty,rlf)); + pm_co2eqForeign(ttot,regi) = sum(regi2 $ (not sameas(regi,regi2)), pm_co2eq0(ttot,regi2)); + pm_fuExtrForeign(ttot,regi,enty,rlf) = sum(regi2 $ (not sameas(regi,regi2)), vm_fuExtr.l(ttot,regi2,enty,rlf)); + pm_capCumForeign(ttot,regi,teLearn) = sum(regi2 $ (not sameas(regi,regi2)), pm_capCum0(ttot,regi2,teLearn)); + ); + +*** If cm_LearningSpillover is 0, the foreign capacity in technology learning is fixed to the level of 2020. +*** This simulates a world of protectionism with no further foreign technology learning spillover. + pm_capCumForeign(ttot,regi,teLearn) $ (ttot.val >= 2025 and cm_LearningSpillover = 0) + = sum(regi2 $ (not sameas(regi,regi2)), pm_capCum0("2020",regi2,teLearn)); +); display pm_capCumForeign, pm_co2eqForeign; *** EOF ./modules/80_optimization/nash/presolve.gms diff --git a/modules/80_optimization/nash/realization.gms b/modules/80_optimization/nash/realization.gms index ede59aab2b..12533d68c8 100644 --- a/modules/80_optimization/nash/realization.gms +++ b/modules/80_optimization/nash/realization.gms @@ -9,15 +9,17 @@ *' @description *' Unlike in Negishi-mode, each region forms its own optimization problem in Nash mode. *' Regions trade on goods and resource markets, but market-clearing conditions are not part of the optimization itself. -*' Instead, the Nash-algorithm iteratively computes solutions for all regions including their trade patterns, and adjusts prices such that the surplus on global markets vanishes. -*' Initial values for trade patterns, prices etc. are taken from the gdx (input.gdx). +*' Instead, the Nash-algorithm iteratively computes solutions for all regions including their trade patterns, and adjusts prices +*' such that the surplus on global markets vanishes. Initial values for trade patterns, prices etc. are taken from the gdx (input.gdx). *' -*' Potential benefits of a Nash-solution are a massive reduction in run-time (convergence within a few hours), and more flexibility in treating inter-regional externalities. -*' Learning-by-doing technologies (learnte) are included by default and cause an inter-regional spill-over. This causes a welfare difference between the solution in Nash- and Negishi-mode. -*' In Nash-mode, a subsidy on the investment cost of learning technologies can be used to internalize this spill-over externality. This subsidy is implemented in the module 22_subsidizeLearning. +*' Potential benefits of a Nash-solution are a massive reduction in run-time (convergence within a few hours), and more flexibility in treating +*' inter-regional externalities. Learning-by-doing technologies (teLearn) are included by default and cause an inter-regional spill-over. +*' This causes a welfare difference between the solution in Nash- and Negishi-mode. +*' In Nash-mode, a subsidy on the investment cost of learning technologies can be used to internalize this spill-over externality. +*' This subsidy is implemented in the module 22_subsidizeLearning. *' -*' Without internalizing the learning-by-doing spill-over due to the global learning curve, Nash and Negishi solution differ. This is the case in the default setting of the corresponding module: -*' cfg$gms$subsidizeLearning <- "off" +*' Without internalizing the learning-by-doing spill-over due to the global learning curve, Nash and Negishi solution differ. +*' This is the case in the default setting of the corresponding module: cfg$gms$subsidizeLearning <- "off" *' In Nash-mode, the subsidy internalizing this externality can be calculated. *' When activated by cfg$gms$subsidizeLearning <- "globallyOptimal" the Nash solution should be equivalent to the Negishi solution. diff --git a/modules/80_optimization/nash/sets.gms b/modules/80_optimization/nash/sets.gms index 9fdaaf9e15..79cceacefb 100644 --- a/modules/80_optimization/nash/sets.gms +++ b/modules/80_optimization/nash/sets.gms @@ -7,39 +7,25 @@ *** SOF ./modules/80_optimization/nash/sets.gms sets -learnte_dyn80(all_te) "learnte for nash" -/ - spv "solar photovoltaic" - csp "concentrating solar power" - windon "wind onshore power converters" - windoff "wind offshore power converters" - storspv "storage technology for spv" - storcsp "storage technology for csp" - storwindon "storage technology for wind onshore" - storwindoff "storage technology for wind offshore" -/, - solveinfo80 "Nash solution stats" / -solvestat, modelstat, resusd, objval + solvestat, modelstat, resusd, objval / convMessage80 "contains all convergence criteria" / -infes,surplus,nonopt,taxconv,anticip,globalbudget,peakbudgyr,peakbudget,regiTarget,implicitEnergyTarget,cm_implicitPriceTarget,cm_implicitPePriceTarget,damage,DevPriceAnticip, IterationNumber + infes, surplus, nonopt, taxconv, anticip, globalbudget, peakbudgyr, peakbudget, regiTarget, + implicitEnergyTarget, cm_implicitPriceTarget, cm_implicitPePriceTarget, damage, DevPriceAnticip, IterationNumber / activeConvMessage80(convMessage80) "all active convergence criterias" / / ; -teLearn(learnte_dyn80) = YES; - activeConvMessage80("infes") = YES; activeConvMessage80("surplus") = YES; activeConvMessage80("nonopt") = YES; activeConvMessage80("IterationNumber") = YES; - -if (cm_TaxConvCheck eq 1, activeConvMessage80("taxconv") = YES;); +if(cm_TaxConvCheck eq 1, activeConvMessage80("taxconv") = YES;); ***activeConvMessage80("anticip") = YES; activeConvMessage80("globalbudget") = YES; $if %carbonprice% == "functionalForm" activeConvMessage80("peakbudgyr") = YES; @@ -51,5 +37,4 @@ $if not "%cm_implicitPriceTarget%" == "off" activeConvMessage80("cm_implicitPric $if not "%cm_implicitPePriceTarget%" == "off" activeConvMessage80("cm_implicitPePriceTarget") = YES; $if not "%internalizeDamages%" == "off" activeConvMessage80("damage") = YES; -display teLearn; *** EOF ./modules/80_optimization/nash/sets.gms diff --git a/modules/80_optimization/nash/solve.gms b/modules/80_optimization/nash/solve.gms index 45e2e9f9eb..7cd433347a 100644 --- a/modules/80_optimization/nash/solve.gms +++ b/modules/80_optimization/nash/solve.gms @@ -10,33 +10,28 @@ regi(all_regi) = NO; hybrid.solvelink = 3; !! activate multiple-CPU mode for GAMS hybrid.optfile = 9; -if(cm_nash_mode eq 1, +if(cm_nash_mode = 1, hybrid.solvelink = 0; !! activate single-CPU mode for GAMS ); -loop (all_regi, - !! only solve for regions that do not have a valid solution from the - !! last solver iteration - if ( ( sol_itr.val gt 1 - OR s80_runInDebug eq 1) - AND ( p80_repy(all_regi,"modelstat") eq 2 -$ifthen.repeatNonOpt "%cm_repeatNonOpt%" == "off" - OR p80_repy(all_regi,"modelstat") eq 7 -$endif.repeatNonOpt - ), - +loop(all_regi, +*** only solve for regions that do not have a valid solution from the last solver iteration + if( ( sol_itr.val > 1 or s80_runInDebug = 1) + and ( p80_repy(all_regi,"modelstat") = 2 +$if "%cm_repeatNonOpt%" == "off" or p80_repy(all_regi,"modelstat") = 7 + ), p80_repy_thisSolitr(all_regi,solveinfo80) = 0; continue; ); regi(all_regi) = YES; - if (execError > 0, + if(execError > 0, execute_unload "abort.gdx"; abort "at least one execution error occured, possibly in the loop"; ); - if (cm_keep_presolve_gdxes eq 1, + if(cm_keep_presolve_gdxes = 1, execute_unload "presolve_nash.gdx"; sm_tmp = logfile.nr; sm_tmp2 = logfile.nd; @@ -44,20 +39,20 @@ $endif.repeatNonOpt logfile.nd = 0; put_utility logfile, "shell" / "mv presolve_nash.gdx presolve_nash_" all_regi.tl "_CES-" - sm_CES_calibration_iteration "_Nash-" iteration.val "_Sol-" sol_itr.val - ".gdx"; + sm_CES_calibration_iteration "_Nash-" iteration.val "_Sol-" sol_itr.val + ".gdx"; logfile.nr = sm_tmp; logfile.nd = sm_tmp2; ); solve hybrid using nlp maximizing vm_welfareGlob; - if(cm_nash_mode eq 1, + if(cm_nash_mode = 1, !! regions run in a sequence p80_repy_thisSolitr(all_regi,"solvestat") = hybrid.solvestat; p80_repy_thisSolitr(all_regi,"modelstat") = hybrid.modelstat; p80_repy_thisSolitr(all_regi,"resusd") = hybrid.resusd; p80_repy_thisSolitr(all_regi,"objval") = hybrid.objval; - if (p80_repy_thisSolitr(all_regi,"modelstat") eq 2, + if(p80_repy_thisSolitr(all_regi,"modelstat") = 2, p80_repyLastOptim(all_regi,"objval") = p80_repy_thisSolitr(all_regi,"objval"); ); ); @@ -66,15 +61,15 @@ $endif.repeatNonOpt p80_handle(all_regi) = hybrid.handle; ); !! close regi loop -if(cm_nash_mode eq 2, +if(cm_nash_mode = 2, !! regions run in parallel repeat - loop (all_regi$handlecollect(p80_handle(all_regi)), + loop(all_regi $ handlecollect(p80_handle(all_regi)), p80_repy_thisSolitr(all_regi,"solvestat") = hybrid.solvestat; p80_repy_thisSolitr(all_regi,"modelstat") = hybrid.modelstat; p80_repy_thisSolitr(all_regi,"resusd") = hybrid.resusd; p80_repy_thisSolitr(all_regi,"objval") = hybrid.objval; - if (p80_repy_thisSolitr(all_regi,"modelstat") eq 2, + if(p80_repy_thisSolitr(all_regi,"modelstat") = 2, p80_repyLastOptim(all_regi,"objval") = p80_repy_thisSolitr(all_regi,"objval"); ); @@ -92,46 +87,40 @@ pm_SolNonInfes(regi) = 0; p80_SolNonOpt(regi) = 0; putclose foo_msg; -*** This putclose serves to make foo_msg the last "active" put file, and thus makes GAMS use the foo_msg formating (namely F-format, not scientific E-format) +*** This putclose makes foo_msg the last "active" put file so that GAMS uses the foo_msg formating (namely F-format, not scientific E-format) *** Otherwise, the following put messages will try to write modelstat in scientif format, throwing errors because of insufficient space - -loop (regi, - if( (p80_repy_thisSolitr(regi,"solvestat") > 0) , +loop(regi, + if(p80_repy_thisSolitr(regi,"solvestat") > 0, put_utility foo_msg "msg" / "Solitr:" sol_itr.tl:2:0 " " regi.tl:4:0 " updated. Modstat new " p80_repy_thisSolitr(regi,"modelstat"):2:0 ", old " p80_repy(regi,"modelstat"):2:0 "; Resusd new" p80_repy_thisSolitr(regi,"resusd"):5:0 ", old" p80_repy(regi,"resusd"):5:0 "; Obj new" p80_repy_thisSolitr(regi,"objval"):7:3 ", old" p80_repy(regi,"objval"):7:3 ; p80_repy(regi,solveinfo80) = p80_repy_thisSolitr(regi,solveinfo80); !! copy info from this Solitr into p80_repy else put_utility foo_msg "msg" / "Solitr:" sol_itr.tl:2:0 " " regi.tl:4:0 " not updated. Modstat new " p80_repy_thisSolitr(regi,"modelstat"):2:0 ", old " p80_repy(regi,"modelstat"):2:0 "; Resusd new" p80_repy_thisSolitr(regi,"resusd"):5:0 ", old" p80_repy(regi,"resusd"):5:0 "; Obj new" p80_repy_thisSolitr(regi,"objval"):7:3 ", old" p80_repy(regi,"objval"):7:3 ; ); - if (p80_repy(regi,"modelstat") eq 2 OR p80_repy(regi,"modelstat") eq 7, + if(p80_repy(regi,"modelstat") = 2 or p80_repy(regi,"modelstat") = 7, pm_SolNonInfes(regi) = 1; ); - if (p80_repy(regi,"modelstat") eq 7, + if(p80_repy(regi,"modelstat") = 7, p80_SolNonOpt(regi) = 1); ); *** set o_modelstat to the highest value across all regions -o_modelstat $ifthen.repeatNonOpt "%cm_repeatNonOpt%" == "off" - = smax(regi, p80_repy(regi,"modelstat")$(p80_repy(regi,"modelstat") ne 7)); !! ignoring status 7 + o_modelstat = smax(regi, p80_repy(regi,"modelstat") $ (p80_repy(regi,"modelstat") ne 7)); !! ignoring status 7 $else.repeatNonOpt - = smax(regi, p80_repy(regi,"modelstat")); !! also taking into account status 7 + o_modelstat = smax(regi, p80_repy(regi,"modelstat")); !! also taking into account status 7 $endif.repeatNonOpt -!! add information if this region was solved in this iteration -p80_repy_iteration(regi,solveinfo80,iteration)$( - p80_repy_thisSolitr(regi,solveinfo80) ) - !! store sum of resusd for all sol_itrs - = ( p80_repy_iteration(regi,solveinfo80,iteration) - + p80_repy_thisSolitr(regi,solveinfo80)$( - p80_repy_thisSolitr(regi,solveinfo80) ne NA ) - )$( sameas(solveinfo80,"resusd") ) - + p80_repy_thisSolitr(regi,solveinfo80)$( NOT sameas(solveinfo80,"resusd") ); - -!! add information if this region was solved in this iteration -p80_repy_nashitr_solitr(regi,solveinfo80,iteration,sol_itr)$( - p80_repy_thisSolitr(regi,solveinfo80) ) - = p80_repy_thisSolitr(regi,solveinfo80); +*** add information if this region was solved in this iteration +loop((regi,solveinfo80) $ p80_repy_thisSolitr(regi,solveinfo80), + p80_repy_iteration(regi,solveinfo80,iteration) !! store sum of resusd for all sol_itrs + = ( p80_repy_iteration(regi,solveinfo80,iteration) + + p80_repy_thisSolitr(regi,solveinfo80) $ (p80_repy_thisSolitr(regi,solveinfo80) ne NA) + ) $ sameas(solveinfo80,"resusd") + + p80_repy_thisSolitr(regi,solveinfo80) $ (not sameas(solveinfo80,"resusd")); + + p80_repy_nashitr_solitr(regi,solveinfo80,iteration,sol_itr) = p80_repy_thisSolitr(regi,solveinfo80); +); put_utility "msg" / "Solve overview: The following are the results for iteration " iteration.tl:3:0 " , sol_itr " sol_itr.tl:3:0 ; display o_modelstat; diff --git a/modules/80_optimization/negishi/not_used.txt b/modules/80_optimization/negishi/not_used.txt index 3e2a255d84..ace563bb1f 100644 --- a/modules/80_optimization/negishi/not_used.txt +++ b/modules/80_optimization/negishi/not_used.txt @@ -74,3 +74,5 @@ sm_peakbudget_diff_tolerance,input,not needed cm_iterative_target_adj,input,not needed sm_peakBudgYr_check,input,not needed sm_peakbudget_diff,input,not needed +sm_MtCO2_2_GtC,input,not needed +cm_startIter_EDGET,input,not needed diff --git a/modules/80_optimization/testOneRegi/bounds.gms b/modules/80_optimization/testOneRegi/bounds.gms index f759d4feb7..0b63376568 100644 --- a/modules/80_optimization/testOneRegi/bounds.gms +++ b/modules/80_optimization/testOneRegi/bounds.gms @@ -6,6 +6,6 @@ *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/testOneRegi/bounds.gms -*AJS* Adjustment costs are only relevant for the Nash_test realization +*AJS* Adjustment costs are only relevant for the nash realization vm_costAdjNash.fx(t,regi) = 0; *** EOF ./modules/80_optimization/testOneRegi/bounds.gms diff --git a/modules/80_optimization/testOneRegi/datainput.gms b/modules/80_optimization/testOneRegi/datainput.gms index 6ea1dea919..c32d2dcee5 100644 --- a/modules/80_optimization/testOneRegi/datainput.gms +++ b/modules/80_optimization/testOneRegi/datainput.gms @@ -10,13 +10,9 @@ pm_w(regi) = 1; *LB* initialize price parameter, import from gdx in preloop pm_pvp(ttot,trade)$((ttot.val ge 2005) AND (NOT tradeSe(trade))) = 0; -p80_etaXp(tradePe) = 1; -p80_etaXp("good") = 1; -p80_etaXp("perm") = 1; - -*load fallback prices - -$include "./modules/80_optimization/testOneRegi/input/prices_NASH.inc"; +p80_priceAnticipStrength(tradePe) = 1; +p80_priceAnticipStrength("good") = 1; +p80_priceAnticipStrength("perm") = 1; *MLB 12062013* initialize learning externality (can be improved by including file) diff --git a/modules/80_optimization/testOneRegi/declarations.gms b/modules/80_optimization/testOneRegi/declarations.gms index ebe385d87d..412364d039 100644 --- a/modules/80_optimization/testOneRegi/declarations.gms +++ b/modules/80_optimization/testOneRegi/declarations.gms @@ -7,23 +7,22 @@ *** SOF ./modules/80_optimization/testOneRegi/declarations.gms parameter -p80_pvpFallback(ttot,all_enty) "Helper parameter. Price path from input/prices_NASH.inc. Only used if reading prices from gdx fails.", -p80_etaXp(all_enty) "Parameter governing price anticipation on commodity markets", -p80_taxrev0(tall,all_regi) "???" -p80_Mport0(tall,all_regi,all_enty) "Imports in last iteration" -p80_normalize0(ttot,all_regi,all_enty) "Normalization parameter for market volume" -pm_cumEff(tall,all_regi,all_in) "parameter for spillover externality (aggregated productivity level)" -pm_fuExtrForeign(ttot,all_regi,all_enty,rlf) "foreign fuel extraction" +p80_priceAnticipStrength(all_enty) "Parameter governing price anticipation on commodity markets" +p80_regiMarketVolume(ttot,all_regi,all_enty) "Regional market volume of a trade item, used for normalisation [amount of trade item]" +p80_Mport0(tall,all_regi,all_enty) "Imports in last iteration" +p80_taxrev0(tall,all_regi) "vm_taxrev from last iteration" +pm_cumEff(tall,all_regi,all_in) "parameter for spillover externality (aggregated productivity level)" +pm_fuExtrForeign(ttot,all_regi,all_enty,rlf) "foreign fuel extraction" ; positive variable -*AJS* Adjustment costs for Nash trade algorithm. Only non-zero in the Nash_test realization of 80_optimization module. -vm_costAdjNash(ttot,all_regi) "Adjustment costs for deviation from the trade structure of the last iteration." +*** adjustment costs for Nash trade algorithm (only non-zero in the nash realization of 80_optimization module) +vm_costAdjNash(ttot,all_regi) "Adjustment costs for deviation from the trade structure of the last iteration." ; equations -q80_budg_intertemp(all_regi) "interemporal trade balance (Nash mode only)" -q80_costAdjNash(ttot,all_regi) "plays a dummy role for now, allowing fixing to Nash GDX files" +q80_budg_intertemp(all_regi) "interemporal trade balance (Nash mode only)" +q80_costAdjNash(ttot,all_regi) "calculate Nash adjustment costs; no role in testOneRegi mode" ; *** EOF ./modules/80_optimization/testOneRegi/declarations.gms diff --git a/modules/80_optimization/testOneRegi/equations.gms b/modules/80_optimization/testOneRegi/equations.gms index ca5777abf7..be9e1b8b60 100644 --- a/modules/80_optimization/testOneRegi/equations.gms +++ b/modules/80_optimization/testOneRegi/equations.gms @@ -8,24 +8,26 @@ *' @equations -*' for Nash solution: intertemporal trade balance must be zero -q80_budg_intertemp(regi).. -0 =e= - SUM(ttot$(ttot.val ge 2005), - pm_ts(ttot) - * ( - SUM(trade$(NOT tradeSe(trade)), - (vm_Xport(ttot,regi,trade)-vm_Mport(ttot,regi,trade)) * pm_pvp(ttot,trade) - * ( 1 + p80_etaXp(trade) - * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) ) - / (p80_normalize0(ttot,regi,trade) + 1E-6) - ) - ) - + vm_capacityTradeBalance(ttot,regi) +*' Intertemporal trade balance must be zero. Terms are: +*' sum over all time steps of (time step duration * net exports in this time step): +*' difference of net exports compared to previous iteration, adjusted for price anticipation +q80_budg_intertemp(regi).. + 0 =e= + sum(ttot $ (ttot.val >= 2005), + pm_ts(ttot) !! duration of the time step (average between previous and next time steps) + * ( vm_capacityTradeBalance(ttot,regi) !! trade balance for 24_trade capacity realisation + + sum(trade $ (not tradeSe(trade) and not tradeCap(trade)), + (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) * pm_pvp(ttot,trade) !! net value of exports + * ( 1 + p80_priceAnticipStrength(trade) + * ( (pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)) ) + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) + ) ) - ); + ) + ); -q80_costAdjNash(ttot,regi)$(ttot.val ge cm_startyear).. -vm_costAdjNash(ttot,regi) =e= 0; + +q80_costAdjNash(ttot,regi) $ (ttot.val ge cm_startyear).. + vm_costAdjNash(ttot,regi) =e= 0; *** EOF ./modules/80_optimization/testOneRegi/equations.gms diff --git a/modules/80_optimization/testOneRegi/input/prices_NASH.inc b/modules/80_optimization/testOneRegi/input/prices_NASH.inc deleted file mode 100644 index c95c3bc655..0000000000 --- a/modules/80_optimization/testOneRegi/input/prices_NASH.inc +++ /dev/null @@ -1,135 +0,0 @@ -*** SOF ./modules/80_optimization/testOneRegi/input/prices_NASH.inc -p80_pvpFallback("2010","peoil")= 0.00228271; -p80_pvpFallback("2015","peoil")= 0.00185881; -p80_pvpFallback("2020","peoil")= 0.00163106; -p80_pvpFallback("2025","peoil")= 0.00132324; -p80_pvpFallback("2030","peoil")= 0.00118051; -p80_pvpFallback("2035","peoil")= 0.00085660; -p80_pvpFallback("2040","peoil")= 0.00067428; -p80_pvpFallback("2045","peoil")= 0.00052235; -p80_pvpFallback("2050","peoil")= 0.00044664; -p80_pvpFallback("2055","peoil")= 0.00034830; -p80_pvpFallback("2060","peoil")= 0.00027406; -p80_pvpFallback("2070","peoil")= 0.00021467; -p80_pvpFallback("2080","peoil")= 0.00014628; -p80_pvpFallback("2090","peoil")= 0.00010435; -p80_pvpFallback("2100","peoil")= 0.00007633; -p80_pvpFallback("2110","peoil")= 0.00005228; -p80_pvpFallback("2130","peoil")= 0.00003439; -p80_pvpFallback("2150","peoil")= 0.00001313; -p80_pvpFallback("2010","pegas")= 0.00094487; -p80_pvpFallback("2015","pegas")= 0.00071265; -p80_pvpFallback("2020","pegas")= 0.00052809; -p80_pvpFallback("2025","pegas")= 0.00041664; -p80_pvpFallback("2030","pegas")= 0.00034030; -p80_pvpFallback("2035","pegas")= 0.00026381; -p80_pvpFallback("2040","pegas")= 0.00020924; -p80_pvpFallback("2045","pegas")= 0.00017118; -p80_pvpFallback("2050","pegas")= 0.00014598; -p80_pvpFallback("2055","pegas")= 0.00012458; -p80_pvpFallback("2060","pegas")= 0.00009798; -p80_pvpFallback("2070","pegas")= 0.00008733; -p80_pvpFallback("2080","pegas")= 0.00006769; -p80_pvpFallback("2090","pegas")= 0.00004785; -p80_pvpFallback("2100","pegas")= 0.00003874; -p80_pvpFallback("2110","pegas")= 0.00002666; -p80_pvpFallback("2130","pegas")= 0.00001820; -p80_pvpFallback("2150","pegas")= 0.00000776; -p80_pvpFallback("2010","pecoal")= 0.00039189; -p80_pvpFallback("2015","pecoal")= 0.00032563; -p80_pvpFallback("2020","pecoal")= 0.00026147; -p80_pvpFallback("2025","pecoal")= 0.00021817; -p80_pvpFallback("2030","pecoal")= 0.00018079; -p80_pvpFallback("2035","pecoal")= 0.00014859; -p80_pvpFallback("2040","pecoal")= 0.00011972; -p80_pvpFallback("2045","pecoal")= 0.00010619; -p80_pvpFallback("2050","pecoal")= 0.00009678; -p80_pvpFallback("2055","pecoal")= 0.00007992; -p80_pvpFallback("2060","pecoal")= 0.00006888; -p80_pvpFallback("2070","pecoal")= 0.00005399; -p80_pvpFallback("2080","pecoal")= 0.00004017; -p80_pvpFallback("2090","pecoal")= 0.00003208; -p80_pvpFallback("2100","pecoal")= 0.00002555; -p80_pvpFallback("2110","pecoal")= 0.00001804; -p80_pvpFallback("2130","pecoal")= 0.00001224; -p80_pvpFallback("2150","pecoal")= 0.00000487; -p80_pvpFallback("2010","peur")= 0.00057982; -p80_pvpFallback("2015","peur")= 0.00047069; -p80_pvpFallback("2020","peur")= 0.00040449; -p80_pvpFallback("2025","peur")= 0.00035830; -p80_pvpFallback("2030","peur")= 0.00032322; -p80_pvpFallback("2035","peur")= 0.00029497; -p80_pvpFallback("2040","peur")= 0.00027181; -p80_pvpFallback("2045","peur")= 0.00025285; -p80_pvpFallback("2050","peur")= 0.00023723; -p80_pvpFallback("2055","peur")= 0.00022431; -p80_pvpFallback("2060","peur")= 0.00021287; -p80_pvpFallback("2070","peur")= 0.00020701; -p80_pvpFallback("2080","peur")= 0.00018858; -p80_pvpFallback("2090","peur")= 0.00017120; -p80_pvpFallback("2100","peur")= 0.00015159; -p80_pvpFallback("2110","peur")= 0.00013052; -p80_pvpFallback("2130","peur")= 0.00010186; -p80_pvpFallback("2150","peur")= 0.00009891; -p80_pvpFallback("2010","pebiolc")= 0.00094461; -p80_pvpFallback("2015","pebiolc")= 0.00063663; -p80_pvpFallback("2020","pebiolc")= 0.00048667; -p80_pvpFallback("2025","pebiolc")= 0.00036216; -p80_pvpFallback("2030","pebiolc")= 0.00027056; -p80_pvpFallback("2035","pebiolc")= 0.00021136; -p80_pvpFallback("2040","pebiolc")= 0.00016270; -p80_pvpFallback("2045","pebiolc")= 0.00017002; -p80_pvpFallback("2050","pebiolc")= 0.00013628; -p80_pvpFallback("2055","pebiolc")= 0.00010850; -p80_pvpFallback("2060","pebiolc")= 0.00008683; -p80_pvpFallback("2070","pebiolc")= 0.00007830; -p80_pvpFallback("2080","pebiolc")= 0.00004762; -p80_pvpFallback("2090","pebiolc")= 0.00004069; -p80_pvpFallback("2100","pebiolc")= 0.00003999; -p80_pvpFallback("2110","pebiolc")= 0.00001986; -p80_pvpFallback("2130","pebiolc")= 0.00001570; -p80_pvpFallback("2150","pebiolc")= 0.00000581; -p80_pvpFallback("2010","good")= 0.01110502; -p80_pvpFallback("2015","good")= 0.00747063; -p80_pvpFallback("2020","good")= 0.00539279; -p80_pvpFallback("2025","good")= 0.00396793; -p80_pvpFallback("2030","good")= 0.00295582; -p80_pvpFallback("2035","good")= 0.00221891; -p80_pvpFallback("2040","good")= 0.00167476; -p80_pvpFallback("2045","good")= 0.00127277; -p80_pvpFallback("2050","good")= 0.00097038; -p80_pvpFallback("2055","good")= 0.00074049; -p80_pvpFallback("2060","good")= 0.00051893; -p80_pvpFallback("2070","good")= 0.00036242; -p80_pvpFallback("2080","good")= 0.00021508; -p80_pvpFallback("2090","good")= 0.00013270; -p80_pvpFallback("2100","good")= 0.00008364; -p80_pvpFallback("2110","good")= 0.00004530; -p80_pvpFallback("2130","good")= 0.00002002; -p80_pvpFallback("2150","good")= 0.00000548; -p80_pvpFallback("2010","perm")= 0.00000000; -p80_pvpFallback("2015","perm")= 0.00000000; -p80_pvpFallback("2020","perm")= 0.00000000; -p80_pvpFallback("2025","perm")= 0.00000000; -p80_pvpFallback("2030","perm")= 0.00000000; -p80_pvpFallback("2035","perm")= 0.00000000; -p80_pvpFallback("2040","perm")= 0.00000000; -p80_pvpFallback("2045","perm")= 0.00000000; -p80_pvpFallback("2050","perm")= 0.00000000; -p80_pvpFallback("2055","perm")= 0.00000000; -p80_pvpFallback("2060","perm")= 0.00000000; -p80_pvpFallback("2070","perm")= 0.00000000; -p80_pvpFallback("2080","perm")= 0.00000000; -p80_pvpFallback("2090","perm")= 0.00000000; -p80_pvpFallback("2100","perm")= 0.00000000; -p80_pvpFallback("2110","perm")= 0.00000000; -p80_pvpFallback("2130","perm")= 0.00000000; -p80_pvpFallback("2150","perm")= 0.00000000; -p80_pvpFallback("2005","peoil")= 0.00000000; -p80_pvpFallback("2005","pegas")= 0.00000000; -p80_pvpFallback("2005","pecoal")= 0.00000000; -p80_pvpFallback("2005","peur")= 0.00056481; -p80_pvpFallback("2005","pebiolc")= 0.00221526; -p80_pvpFallback("2005","good")= 0.01274436; -p80_pvpFallback("2005","perm")= 0.00000000; -*** EOF ./modules/80_optimization/testOneRegi/input/prices_NASH.inc diff --git a/modules/80_optimization/testOneRegi/not_used.txt b/modules/80_optimization/testOneRegi/not_used.txt index 7dbe76ea04..ef87c65622 100644 --- a/modules/80_optimization/testOneRegi/not_used.txt +++ b/modules/80_optimization/testOneRegi/not_used.txt @@ -18,7 +18,6 @@ cm_nash_mode, switch, not needed cm_nash_autoconverge, switch, ??? pm_emiMktTarget_tolerance, switch, not needed cm_implicitQttyTarget_tolerance,switch, not needed -sm_eps, scalar, ??? pm_co2eq0, parameter, ??? pm_capCum0, parameter, ??? pm_NXagr, parameter, ??? @@ -69,3 +68,5 @@ sm_peakbudget_diff_tolerance,input,not needed cm_iterative_target_adj,input,not needed sm_peakBudgYr_check,input,not needed sm_peakbudget_diff,input,not needed +sm_MtCO2_2_GtC,input,not needed +cm_startIter_EDGET,input,not needed diff --git a/modules/80_optimization/testOneRegi/postsolve.gms b/modules/80_optimization/testOneRegi/postsolve.gms index e7275a371a..431d1c762f 100644 --- a/modules/80_optimization/testOneRegi/postsolve.gms +++ b/modules/80_optimization/testOneRegi/postsolve.gms @@ -8,14 +8,12 @@ o_iterationNumber = iteration.val; -*AJS* feed updated prices and quantities into the next iteration: -loop(trade$(NOT tradeSe(trade)), - loop(regi, - loop(ttot$(ttot.val ge cm_startyear), - pm_Xport0(ttot,regi,trade)$(pm_SolNonInfes(regi) eq 1) = vm_Xport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade)$(pm_SolNonInfes(regi) eq 1) = vm_Mport.l(ttot,regi,trade); - ); - ); +*** AJS: feed updated prices and quantities into the next iteration, ML: adjustments in case of infeasibilities (increase import) +loop(trade $ (not tradeSe(trade)), + loop(ttot $ (ttot.val >= cm_startyear), + pm_Xport0(ttot,regi,trade) $ (pm_SolNonInfes(regi) = 1) = vm_Xport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) $ (pm_SolNonInfes(regi) = 1) = vm_Mport.l(ttot,regi,trade); + ); ); *** Save FE prices in each iteration for easier monitoring diff --git a/modules/80_optimization/testOneRegi/preloop.gms b/modules/80_optimization/testOneRegi/preloop.gms index f43337a8ba..16d45d54f3 100644 --- a/modules/80_optimization/testOneRegi/preloop.gms +++ b/modules/80_optimization/testOneRegi/preloop.gms @@ -5,58 +5,39 @@ *** | REMIND License Exception, version 1.0 (see LICENSE file). *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/testOneRegi/preloop.gms -*MLB/AG* for testOneRegi algorithm read initial price data from gdx - Execute_Loadpoint 'input' pm_pvp = pm_pvp; - Execute_Loadpoint 'input' vm_Xport.l = vm_Xport.l; - Execute_Loadpoint 'input' vm_Mport.l = vm_Mport.l; - Execute_Loadpoint 'input' vm_cons.l = vm_cons.l; - Execute_Loadpoint 'input' vm_taxrev.l = vm_taxrev.l; - Execute_Loadpoint 'input' vm_fuExtr.l = vm_fuExtr.l; - Execute_Loadpoint 'input' vm_prodPe.l = vm_prodPe.l; - Execute_Loadpoint 'input' pm_capCumForeign = pm_capCumForeign; +***------------------------------------------------------------------------------ +*** Use gdx to initialise prices, trade volumes etc +***------------------------------------------------------------------------------ + +*** read data from gdx with explicit instruction execute_loadpoint +Execute_Loadpoint 'input' pm_pvp = pm_pvp; +Execute_Loadpoint 'input' vm_Xport.l = vm_Xport.l; +Execute_Loadpoint 'input' vm_Mport.l = vm_Mport.l; +Execute_Loadpoint 'input' vm_cons.l = vm_cons.l; +Execute_Loadpoint 'input' vm_taxrev.l = vm_taxrev.l; +Execute_Loadpoint 'input' vm_fuExtr.l = vm_fuExtr.l; +Execute_Loadpoint 'input' vm_prodPe.l = vm_prodPe.l; +Execute_Loadpoint 'input' pm_capCumForeign = pm_capCumForeign; + +loop(ttot $ (ttot.val >= 2005), + loop(trade $ (not tradeSe(trade)), +*** initialise trade volumes from gdx, and set missing values to zero + pm_Xport0(ttot,regi,trade) = vm_Xport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) = vm_Mport.l(ttot,regi,trade); + + p80_regiMarketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_regiMarketVolume(ttot,regi,"perm") = abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)); + p80_regiMarketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; + p80_regiMarketVolume(ttot,regi,trade) = max(sm_eps, p80_regiMarketVolume(ttot,regi,trade)); !! ensure market volume is positive + + p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); + ); +); -*AJS* initialize starting points for prices, trade volumes etc. from gdx. -***in order to read parameters like p80_priceXXX from a gdx, instead of only variables , we have to explicitly instruct gams to do so in the execute_loadpoint command in core/preloop.gms. -***the price paths are the trickiest part. Try to find p80_priceXXX prices in the gdx fist, if that fails, fallback to price path read from input/prices_NASH.inc -loop(ttot$(ttot.val ge 2005), - loop(trade$(NOT tradeSe(trade)), -$ontext - if((pm_pvp(ttot,trade) eq NA) OR (pm_pvp(ttot,trade) lt 1E-12) OR (pm_pvp(ttot,trade) gt 0.1) , -***in case we have not been able to read price paths from the gdx, or these price are zero (or eps), fall back to price paths in file input/prices_NASH.inc: - pm_pvp(ttot,trade) = p80_pvpFallback(ttot,trade); - display 'Nash: Info: Could not load useful initial price from gdx, falling back to the one found in input/prices_NASH.inc. This should not be a problem, the runs can stil converge. '; - ); -***in case pm_pvp is not found in gdx: - if(pm_pvp(ttot,trade) eq NA, - pm_pvp(ttot,trade) = 0; - ); -$offtext - loop(regi, - pm_Xport0(ttot,regi,trade)$(NOT tradeSe(trade)) = vm_Xport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade)$(NOT tradeSe(trade)) = vm_Mport.l(ttot,regi,trade); - -***in case xport/mport is not found in gdx: -$ontext - if(pm_Xport0(ttot,regi,trade) eq NA, - pm_Xport0(ttot,regi,trade) = 0; - vm_Xport.L(ttot,regi,trade) = 0; - ); - if(p80_Mport0(ttot,regi,trade) eq NA, - p80_Mport0(ttot,regi,trade) = 0; - vm_Mport.L(ttot,regi,trade) = 0; - ); -$offtext - - p80_normalize0(ttot,regi,"good") = vm_cons.l(ttot,regi); -*** p80_normalize0(ttot,regi,"perm") = vm_cons.l(ttot,regi); - p80_normalize0(ttot,regi,"perm")$(ttot.val ge 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1E-6); - p80_normalize0(ttot,regi,tradePe) = 0.5 * (sum(rlf,vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)); - -p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); - - ); - ); +*** AJS: for resource with price zero, fall back to previous period in order to avoid convergence problems (seen for peur in 2150) +loop(ttot $ (ttot.val > 2005), + pm_pvp(ttot,tradePe) $ (pm_pvp(ttot,tradePe) = 0) = pm_pvp(ttot-1,tradePe); ); display "info: starting from this price path"; diff --git a/modules/80_optimization/testOneRegi/realization.gms b/modules/80_optimization/testOneRegi/realization.gms index d45e43c66d..ef4f653779 100644 --- a/modules/80_optimization/testOneRegi/realization.gms +++ b/modules/80_optimization/testOneRegi/realization.gms @@ -9,7 +9,7 @@ *' @description This is a reduced model version, only containing one region. *' It is equivalent to the Nash realization with just one region. Prices of resources and goods are exogenously fixed to the values taken from the gdx. *' -*' Run time is some minutes. +*' Run time is a few minutes. *' *' @limitations This realization is only useful for testing purposes or for regions with no influence on the global prices of traded goods.