From 05d7a6daec74eeaf0670628b826bf46f288c7e75 Mon Sep 17 00:00:00 2001 From: lecfab Date: Mon, 25 Aug 2025 12:06:22 +0200 Subject: [PATCH 1/8] refactor 80_optimization/nash --- .../iea2014/declarations.gms | 8 +- modules/80_optimization/nash/bounds.gms | 21 +- modules/80_optimization/nash/datainput.gms | 115 ++-- modules/80_optimization/nash/declarations.gms | 123 ++-- modules/80_optimization/nash/equations.gms | 60 +- modules/80_optimization/nash/postsolve.gms | 586 +++++++++--------- modules/80_optimization/nash/preloop.gms | 122 ++-- modules/80_optimization/nash/presolve.gms | 36 +- modules/80_optimization/nash/realization.gms | 16 +- modules/80_optimization/nash/sets.gms | 23 +- modules/80_optimization/nash/solve.gms | 66 +- .../testOneRegi/declarations.gms | 2 +- .../80_optimization/testOneRegi/equations.gms | 2 +- .../80_optimization/testOneRegi/preloop.gms | 22 +- 14 files changed, 571 insertions(+), 631 deletions(-) diff --git a/modules/04_PE_FE_parameters/iea2014/declarations.gms b/modules/04_PE_FE_parameters/iea2014/declarations.gms index 9d075aa5a8..8f0f3d65d3 100644 --- a/modules/04_PE_FE_parameters/iea2014/declarations.gms +++ b/modules/04_PE_FE_parameters/iea2014/declarations.gms @@ -7,10 +7,10 @@ *** SOF ./modules/04_PE_FE_parameters/iea2014/declarations.gms parameter -pm_IO_input(all_regi,all_enty,all_enty,all_te) "Energy input based on IEA data" +pm_IO_input(all_regi,all_enty,all_enty,all_te) "Energy input based on IEA data" 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" @@ -20,8 +20,8 @@ pm_histfegrowth(all_regi,all_enty) "average growth p04_prodCoupleGlob(all_enty,all_enty,all_te,all_enty) "global couple products" -p04_IO_output_beforeFix(ttot,all_regi,all_enty,all_enty,all_te) "Energy output based on IEA data as read in from input data before correction from FE trajectories" -p04_IO_output_beforeFix_Total(ttot,all_regi,all_enty) "Energy output based on IEA data as read in from input data before correction from FE trajectories summed over SE" +p04_IO_output_beforeFix(ttot,all_regi,all_enty,all_enty,all_te) "Energy output based on IEA data as read in from input data before correction from FE trajectories" +p04_IO_output_beforeFix_Total(ttot,all_regi,all_enty) "Energy output based on IEA data as read in from input data before correction from FE trajectories summed over SE" ; *** EOF ./modules/04_PE_FE_parameters/iea2014/declarations.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..49d53afae0 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; +*** Negishi weights, not used in Nash +pm_w(regi) = 1; -*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 - ); +*** convergence with trade surplus thresholds +if(cm_nash_autoconverge > 0, + cm_iteration_max = 100; !! set max number of iterations + +*** default values for cm_nash_autoconverge = 1 - coarse + 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") = 350 * 12/44 / 1000; !! 350 MtCO2eq, converted into internal unit GtC + + loop(trade $ (tradePe(trade) or tradeMacro(trade)), +*** convergences thresholds - fine + p80_surplusMaxTolerance(trade) $ (cm_nash_autoconverge = 2) = 0.2 * p80_surplusMaxTolerance(trade); +*** convergences thresholds - very coarse + 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 +*** --------------- 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: +*** - if too low, markets jump far away from clearance. +*** - if too high, changes in trade patten over iterations are very slow, convergence takes many many iterations. p80_etaAdj(tradePe) = 80; p80_etaAdj("good") = 100; p80_etaAdj("perm") = 10; -*LB* parameter for nash price algorithm within the optimization. +*** LB: parameter for nash price algorithm within the optimization. +*** Multiplicator to price anticipation p80_etaXp(tradePe) = 0.1; p80_etaXp("good") = 0.1; p80_etaXp("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. +*** LB: parameter for Nash price algorithm between different iterations. These parameters are pretty sensitive: +*** - if market surpluses diverge, try higher values (up to 1). +*** - if surpluses oscillate, try lower values. + +*** short term price ajustment elasticity p80_etaST(tradePe) = 0.3; +p80_etaST("pebiolc") = 0.8; !! AJS: bio market seems to like this +p80_etaST("peur") = 0.2; !! uranium market is more sensitive, so choose lower etaST p80_etaST("good") = 0.25; p80_etaST("perm") = 0.3; +$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. -$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. +*** long term price ajustment elasticity +p80_etaLT(trade) = 0; +p80_etaLT("perm") = 0.03; -*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; +*** --------------- 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. + +*** 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,23 @@ 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 +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..9bd2a39944 100644 --- a/modules/80_optimization/nash/declarations.gms +++ b/modules/80_optimization/nash/declarations.gms @@ -7,32 +7,29 @@ *** SOF ./modules/80_optimization/nash/declarations.gms parameter -*LB* parameters for ajustments within one iteration. These cause price anticipation +*** LB: parameters for ajustments within one iteration. These cause price anticipation p80_etaXp(all_enty) "Parameter governing price anticipation on commodity markets" - -*LB* parameters for ajustments between different iterations -p80_etaLT(all_enty) "long term price ajustment elasticity " +*** 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" - -*AJS* adjustment costs between iterations p80_etaAdj(all_enty) "Adjustment costs for changes of trade pattern between iterations" -***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.", +*** 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." -p80_normalizeLT(all_enty) "Aggregated intertemporal market volume", -p80_normalize0(ttot,all_regi,all_enty) "Normalization parameter for market volume" +p80_itertemporalMarketRevenue(all_enty) "Aggregated intertemporal market volume" +p80_marketVolume(ttot,all_regi,all_enty) "Normalization parameter for market volume" -***parameter containing the respective level values from last iteration (the first set of values taken from gdx in the first iteration, respectively) +*** 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_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 [%]" -*LB* diagnostic parameters +*** 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" @@ -42,20 +39,20 @@ o80_trackSurplusSign(ttot,all_enty,iteration) "auxiliary parameter t 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, 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_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_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 +62,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 +105,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." +*** 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." ; 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 +129,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: +*** defining 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..b0187f731b 100644 --- a/modules/80_optimization/nash/equations.gms +++ b/modules/80_optimization/nash/equations.gms @@ -7,47 +7,41 @@ *** 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) 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) !! only active in 24_trade capacity mode + + 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_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) - vm_taxrev(ttot,regi)) $ (sameas(trade,"good") and ttot.val > 2005)) + / (p80_marketVolume(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) - ) -; +q80_costAdjNash(ttot,regi) $ (ttot.val >= cm_startyear) .. + vm_costAdjNash(ttot,regi) + =e= sum(trade $ (not tradeSe(trade)), + p80_etaAdj(trade) + * pm_pvp(ttot,trade) + * ((pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade))) ** 2 + / (p80_marketVolume(ttot,regi,trade) + sm_eps)); -*' link between permit budget and emission budget +*' 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)); + 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 2beabf08e9..c6479a0fdc 100644 --- a/modules/80_optimization/nash/postsolve.gms +++ b/modules/80_optimization/nash/postsolve.gms @@ -7,290 +7,263 @@ *** 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 prices that clears the markets. +*' 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_etaST_correct = p80_etaST * p80_surplus / p80_marketVolume + +*' Compute market volume for different trades (take values from last iteration for regions that were not solved optimally) +loop(ttot $ (ttot.val >= 2005), + loop(regi $ (pm_SolNonInfes(regi) = 1), + p80_marketVolume(ttot,regi,"good") = max(sm_eps, vm_cons.l(ttot,regi)); + p80_marketVolume(ttot,regi,tradePe) = max(sm_eps, (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2); + ); + loop(regi $ (pm_SolNonInfes(regi) = 0), + p80_marketVolume(ttot,regi,"good") = max(sm_eps, p80_marketVolume(ttot,regi,"good")); + p80_marketVolume(ttot,regi,tradePe) = max(sm_eps, p80_marketVolume(ttot,regi,tradePe)); ); -p80_DevPriceAnticipGlobAll(ttot) = sum(trade$(NOT tradeSe(trade)), p80_DevPriceAnticipGlob(ttot,trade)); +*** ML: normalize permit trade corrections to consumption or positive cap path instead of emissions, as those may be negative + p80_marketVolume(ttot,regi,"perm") = max(sm_eps, abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050"))); ); -*' 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) ) -; +loop(trade $ (not tradeSe(trade)), +*' Calculate residual surplus on the markets + 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 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_marketVolume) + p80_intertemporalSurplusRevenue(trade,iteration) = sum(ttot $ (ttot.val >= cm_startyear), + pm_pvp(ttot,trade) * pm_ts(ttot) * p80_surplus(ttot,trade,iteration)); + p80_itertemporalMarketRevenue(trade) = max(sm_eps, sum((ttot, regi) $ (ttot.val >= 2005), + pm_pvp(ttot,trade) * pm_ts(ttot) * p80_marketVolume(ttot,regi,trade))); + + p80_etaLT_correct(trade,iteration) = + p80_etaLT(trade) * p80_intertemporalSurplusRevenue(trade,iteration) / p80_itertemporalMarketRevenue(trade); + +*' Short term 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_etaST_correct(ttot,trade,iteration) $ (ttot.val >= 2005) = + p80_etaST(trade) * p80_surplus(ttot,trade,iteration) / max(sm_eps , sum(regi, p80_marketVolume(ttot,regi,trade))); + + p80_etaST_correct(ttot,"perm",iteration) $ (ttot.val >= 2005) = p80_etaST_correct(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_etaST_correct(ttot,tradePe(trade),iteration) $ (ttot.val >= 2005) = p80_etaST_correct(ttot,trade,iteration) + * (sm_fadeoutPriceAnticip + (1-sm_fadeoutPriceAnticip) * (pm_pvp(ttot,trade) / pm_pvp("2050",trade))); +); -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"); +*' If the surplus remains over several iterations, increase the price correction terms +p80_etaST_correct_safecopy(ttot,trade,iteration) $ (not tradeSe(trade)) = p80_etaST_correct(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); +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 -*' 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) ) ); - -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; + if(iteration.val > 15 and o80_trackSurplusSign(ttot,trade,iteration) >= 5, !! if surplus was beyond tolerance 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 > 20 and o80_trackSurplusSign(ttot,trade,iteration) >= 10, !! push stronger if previous increase did not help after a few 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 > 25 and o80_trackSurplusSign(ttot,trade,iteration) >= 15, !! push stronger if previous increase did not help after a few iterations + p80_etaST_correct(ttot,trade,iteration) = 2 * p80_etaST_correct(ttot,trade,iteration); + o80_counter_iteration_trade_ttot(ttot,trade,iteration) = 3; + ); + ); !! 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_etaLT_correct(trade,iteration) - p80_etaST_correct(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)) / sum(regi, p80_marketVolume(ttot2,regi,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_itertemporalMarketRevenue("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_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) - vm_taxrev.l(ttot,regi)) $ (ttot.val > 2005 and sameas(trade,"good")) + ) + / (p80_marketVolume(ttot,regi,trade) + sm_eps); -*** 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_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) ) ); + ); + 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"); + + +*' For display of price change p80_PriceChangePriceAnticipReg, round to 0.1% +o80_PriceChangePriceAnticipReg(ttot,trade,regi) = round(p80_PriceChangePriceAnticipReg(ttot,trade,regi), 1); -***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; +*' 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 +273,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 applied 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 "Did REMIND run sufficient iterations": currently set at 18, to allow for at least 4 iterations with EDGE-T +if(iteration.val <= 17, + s80_bool = 0; p80_messageShow("IterationNumber") = YES; ); -***additional criterion: did taxes converge? (only checked if cm_TaxConvCheck is 1) +*' additional 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? +*' additional 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? +*' additional 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 +333,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? +*' additional 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 +343,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? +*' additional 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,23 +353,28 @@ loop((t,regi,entyPe)$pm_implicitPePriceTarget(t,regi,entyPe), ); $endIf.cm_implicitPePriceTarget -*** check global budget target from core/postsolve, must be within 2 Gt of target value +*' additional criterion: global budget target from core/postsolve must be within 2 Gt of target value p80_globalBudget_absDev_iter(iteration) = sm_globalBudget_absDev; -if (abs(p80_globalBudget_absDev_iter(iteration)) gt cm_budgetCO2_absDevTol, +if (abs(p80_globalBudget_absDev_iter(iteration)) > cm_budgetCO2_absDevTol, s80_bool = 0; p80_messageShow("target") = YES; ); -*** additional criterion: if damage internalization is on, is damage iteration converged? +*' additional 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: "; @@ -529,12 +507,12 @@ display o80_trackSurplusSign, o80_SurplusOverTolerance, o80_counter_iteration_t ***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 +if( (s80_bool = 0) and (iteration.val = 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 "#### !! 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"), @@ -622,8 +600,12 @@ $endIf.cm_implicitPePriceTarget ); +***------------------------------------------------------------------------------ +*' #### Finishing or aborting +***------------------------------------------------------------------------------ + ***if all conditions are met, stop optimization. -if(s80_bool eq 1, +if(s80_bool = 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; @@ -654,25 +636,21 @@ if(s80_bool eq 1, ); -*** 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 @@ -680,13 +658,13 @@ 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, + 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."; @@ -695,18 +673,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; @@ -725,7 +701,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); ); @@ -757,8 +733,4 @@ display p80_eoMargEmiCum, p80_eoMargPermBudg, p80_eoEmiMarg, p80_eoMargAverage, ); $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..a1fe4f538c 100644 --- a/modules/80_optimization/nash/preloop.gms +++ b/modules/80_optimization/nash/preloop.gms @@ -5,97 +5,83 @@ *** | 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; + +*** 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; +*** 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 +*** using execute_loadpoint in core/preloop.gms. The price paths are the trickiest part. Try to find p80_priceXXX prices in the gdx first, +*** if that fails, fallback to price path read from input/prices_NASH.inc +loop(ttot $ (ttot.val >= 2005), + loop(trade $ (not tradeSe(trade)), + if(pm_pvp(ttot,trade) = NA or pm_pvp(ttot,trade) < 1e-12 or pm_pvp(ttot,trade) > 0.1, !! in case price paths from the gdx are not valid + 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. '; + ); + + if(pm_pvp(ttot,trade) = NA, pm_pvp(ttot,trade) = 0;); !! in case pm_pvp is not found in gdx -*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, + loop(regi, + pm_Xport0(ttot,regi,trade) = vm_Xport.l(ttot,regi,trade); + p80_Mport0(ttot,regi,trade) = vm_Mport.l(ttot,regi,trade); + + if(pm_Xport0(ttot,regi,trade) = NA, !! in case export is not found in gdx pm_Xport0(ttot,regi,trade) = 0; vm_Xport.L(ttot,regi,trade) = 0; - ); - if(p80_Mport0(ttot,regi,trade) eq NA, + ); + if(p80_Mport0(ttot,regi,trade) = NA, !! in case import is not found in gdx 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)); - -p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); - - ); ); -); - -loop(regi, - loop(tradePe, -if(p80_Mport0("2005",regi,tradePe) eq NA, p80_Mport0("2005",regi,tradePe) = 0); -);); + p80_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_marketVolume(ttot,regi,"perm") $ (ttot.val >= 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1e-6); + p80_marketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; -*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); + ); + ); ); +p80_Mport0("2005",regi,tradePe) $ (p80_Mport0("2005",regi,tradePe) = NA) = 0; -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); +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: 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(ttot $ (not sameas(ttot,"2005")), + p80_pvp_itr(ttot,tradePe,"1") $ (p80_pvp_itr(ttot,tradePe,"1") = 0) = p80_pvp_itr(ttot-1,tradePe,"1"); ); -***debug display -display pm_pvp,p80_normalize0; +*** debug display +display pm_pvp,p80_marketVolume; display pm_Xport0,p80_Mport0; display p80_surplusMaxTolerance; -*EMIOPT +*** 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 08c40b7541..e8862bdc02 100644 --- a/modules/80_optimization/nash/presolve.gms +++ b/modules/80_optimization/nash/presolve.gms @@ -5,22 +5,32 @@ *** | 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 +*** 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. +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: take care of the case of infeasible solution -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_co2eqForeign(ttot,regi) + = sum(regi2 $ (not sameas(regi,regi2)), pm_co2eq0(ttot,regi2)); !! does this interfere with the initialization in datainput? -pm_emissionsForeign(ttot,regi,enty)$((ttot.val ge 2005) and (pm_SolNonInfes(regi) eq 1)) = sum(regi2$((NOT sameas(regi,regi2))), pm_emissions0(ttot,regi2,enty)); + pm_emissionsForeign(ttot,regi,enty) + = sum(regi2 $ (not sameas(regi,regi2)), pm_emissions0(ttot,regi2,enty)); -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_fuExtrForeign(ttot,regi,enty,rlf) + = sum(regi2 $ (not sameas(regi,regi2)), vm_fuExtr.l(ttot,regi2,enty,rlf)); -display pm_capCumForeign, pm_co2eqForeign,pm_emissionsForeign; + 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, pm_emissionsForeign; *** 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 6d12f52df1..85bb20509a 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,target,regiTarget,implicitEnergyTarget,cm_implicitPriceTarget,cm_implicitPePriceTarget,damage,DevPriceAnticip, IterationNumber + infes, surplus, nonopt, taxconv, IterationNumber, anticip, target, DevPriceAnticip, regiTarget, + implicitEnergyTarget, cm_implicitPriceTarget, cm_implicitPePriceTarget, damage / 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("target") = YES; activeConvMessage80("DevPriceAnticip") = YES; @@ -49,5 +35,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..cb1d6b07eb 100644 --- a/modules/80_optimization/nash/solve.gms +++ b/modules/80_optimization/nash/solve.gms @@ -10,20 +10,18 @@ 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 +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 $ifthen.repeatNonOpt "%cm_repeatNonOpt%" == "off" - OR p80_repy(all_regi,"modelstat") eq 7 + or p80_repy(all_regi,"modelstat") = 7 $endif.repeatNonOpt - ), + ), p80_repy_thisSolitr(all_regi,solveinfo80) = 0; continue; @@ -31,12 +29,12 @@ $endif.repeatNonOpt 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 +42,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 +64,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"); ); @@ -95,18 +93,18 @@ 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) *** 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); ); @@ -118,20 +116,16 @@ $else.repeatNonOpt = 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/testOneRegi/declarations.gms b/modules/80_optimization/testOneRegi/declarations.gms index ebe385d87d..368278d692 100644 --- a/modules/80_optimization/testOneRegi/declarations.gms +++ b/modules/80_optimization/testOneRegi/declarations.gms @@ -11,7 +11,7 @@ p80_pvpFallback(ttot,all_enty) "Helper parameter. Price pat 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" +p80_marketVolume(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" ; diff --git a/modules/80_optimization/testOneRegi/equations.gms b/modules/80_optimization/testOneRegi/equations.gms index ca5777abf7..7f9aec5223 100644 --- a/modules/80_optimization/testOneRegi/equations.gms +++ b/modules/80_optimization/testOneRegi/equations.gms @@ -18,7 +18,7 @@ q80_budg_intertemp(regi).. (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) + / (p80_marketVolume(ttot,regi,trade) + 1E-6) ) ) + vm_capacityTradeBalance(ttot,regi) diff --git a/modules/80_optimization/testOneRegi/preloop.gms b/modules/80_optimization/testOneRegi/preloop.gms index f43337a8ba..ba3fac09af 100644 --- a/modules/80_optimization/testOneRegi/preloop.gms +++ b/modules/80_optimization/testOneRegi/preloop.gms @@ -20,7 +20,7 @@ ***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)), + 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: @@ -32,9 +32,9 @@ $ontext 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); + 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 @@ -48,15 +48,13 @@ $ontext ); $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_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_marketVolume(ttot,regi,"perm")$(ttot.val ge 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1E-6); + p80_marketVolume(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); - - ); - ); + p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); + ); + ); ); display "info: starting from this price path"; From 5c664e68d299fca0939f6e74d254940479563540 Mon Sep 17 00:00:00 2001 From: lecfab Date: Fri, 29 Aug 2025 13:44:47 +0200 Subject: [PATCH 2/8] nash: replacing ** by power --- modules/80_optimization/nash/declarations.gms | 5 ++-- modules/80_optimization/nash/equations.gms | 2 +- modules/80_optimization/nash/postsolve.gms | 23 ++++++++++--------- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/modules/80_optimization/nash/declarations.gms b/modules/80_optimization/nash/declarations.gms index 9bd2a39944..75c8b27e03 100644 --- a/modules/80_optimization/nash/declarations.gms +++ b/modules/80_optimization/nash/declarations.gms @@ -19,8 +19,9 @@ p80_etaAdj(all_enty) "Adjustment costs for changes of tra 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." -p80_itertemporalMarketRevenue(all_enty) "Aggregated intertemporal market volume" -p80_marketVolume(ttot,all_regi,all_enty) "Normalization parameter for market volume" +p80_marketVolume(ttot,all_regi,all_enty) "Normalization parameter for market volume [amount of trade item]" +p80_intertemporalSurplusRevenue(all_enty) "Aggregated intertemporal value of the surplus [T$2017]" +p80_itertemporalMarketRevenue(all_enty) "Aggregated intertemporal market volume [T$2017]" *** 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" diff --git a/modules/80_optimization/nash/equations.gms b/modules/80_optimization/nash/equations.gms index b0187f731b..7f2465fd80 100644 --- a/modules/80_optimization/nash/equations.gms +++ b/modules/80_optimization/nash/equations.gms @@ -34,7 +34,7 @@ q80_costAdjNash(ttot,regi) $ (ttot.val >= cm_startyear) .. =e= sum(trade $ (not tradeSe(trade)), p80_etaAdj(trade) * pm_pvp(ttot,trade) - * ((pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade))) ** 2 + * power((pm_Xport0(ttot,regi,trade) - p80_Mport0(ttot,regi,trade)) - (vm_Xport(ttot,regi,trade) - vm_Mport(ttot,regi,trade)), 2) / (p80_marketVolume(ttot,regi,trade) + sm_eps)); *' link between permit budget and emission budget diff --git a/modules/80_optimization/nash/postsolve.gms b/modules/80_optimization/nash/postsolve.gms index 53b126729a..a52b1f547a 100644 --- a/modules/80_optimization/nash/postsolve.gms +++ b/modules/80_optimization/nash/postsolve.gms @@ -12,7 +12,7 @@ *' #### Price corrections to improve convergence of next iteration ***------------------------------------------------------------------------------ -*' The objective of the nash optimisation is to find a set of prices that clears the markets. +*' The objective of the nash optimisation is to find a set of prices that clears the markets for all trade items. *' 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_etaST_correct = p80_etaST * p80_surplus / p80_marketVolume @@ -20,34 +20,35 @@ *' Compute market volume for different trades (take values from last iteration for regions that were not solved optimally) loop(ttot $ (ttot.val >= 2005), loop(regi $ (pm_SolNonInfes(regi) = 1), - p80_marketVolume(ttot,regi,"good") = max(sm_eps, vm_cons.l(ttot,regi)); - p80_marketVolume(ttot,regi,tradePe) = max(sm_eps, (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2); + p80_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); + p80_marketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; ); loop(regi $ (pm_SolNonInfes(regi) = 0), - p80_marketVolume(ttot,regi,"good") = max(sm_eps, p80_marketVolume(ttot,regi,"good")); - p80_marketVolume(ttot,regi,tradePe) = max(sm_eps, p80_marketVolume(ttot,regi,tradePe)); + p80_marketVolume(ttot,regi,"good") = p80_marketVolume(ttot,regi,"good"); + p80_marketVolume(ttot,regi,tradePe) = p80_marketVolume(ttot,regi,tradePe); ); *** ML: normalize permit trade corrections to consumption or positive cap path instead of emissions, as those may be negative - p80_marketVolume(ttot,regi,"perm") = max(sm_eps, abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050"))); + p80_marketVolume(ttot,regi,"perm") = abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050")); + p80_marketVolume(ttot,regi,trade) = max(sm_eps, p80_marketVolume(ttot,regi,trade)) !! ensure market volume is positive ); loop(trade $ (not tradeSe(trade)), -*' Calculate residual surplus on the markets +*' 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 correction takes into account the aggregated intertemporal market revenue (instead of volume) defined by +*' 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_marketVolume) - p80_intertemporalSurplusRevenue(trade,iteration) = sum(ttot $ (ttot.val >= cm_startyear), + p80_intertemporalSurplusRevenue(trade) = sum(ttot $ (ttot.val >= cm_startyear), pm_pvp(ttot,trade) * pm_ts(ttot) * p80_surplus(ttot,trade,iteration)); p80_itertemporalMarketRevenue(trade) = max(sm_eps, sum((ttot, regi) $ (ttot.val >= 2005), pm_pvp(ttot,trade) * pm_ts(ttot) * p80_marketVolume(ttot,regi,trade))); p80_etaLT_correct(trade,iteration) = - p80_etaLT(trade) * p80_intertemporalSurplusRevenue(trade,iteration) / p80_itertemporalMarketRevenue(trade); + p80_etaLT(trade) * p80_intertemporalSurplusRevenue(trade) / p80_itertemporalMarketRevenue(trade); -*' Short term correction takes into account the market surplus volume of a single time step +*' 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_etaST_correct(ttot,trade,iteration) $ (ttot.val >= 2005) = p80_etaST(trade) * p80_surplus(ttot,trade,iteration) / max(sm_eps , sum(regi, p80_marketVolume(ttot,regi,trade))); From adb8420361dce595d0f1b54f4d57b29997dfdf3c Mon Sep 17 00:00:00 2001 From: lecfab Date: Mon, 22 Dec 2025 17:31:03 +0100 Subject: [PATCH 3/8] nash refactor: proofreed all, rename "eta" variables, update testOneRegi --- core/bounds.gms | 6 +- core/declarations.gms | 2 +- main.gms | 6 +- modules/80_optimization/module.gms | 2 +- modules/80_optimization/nash/datainput.gms | 59 ++- modules/80_optimization/nash/declarations.gms | 55 +- modules/80_optimization/nash/equations.gms | 28 +- modules/80_optimization/nash/postsolve.gms | 498 +++++++++--------- modules/80_optimization/nash/preloop.gms | 66 ++- modules/80_optimization/nash/presolve.gms | 16 +- modules/80_optimization/nash/solve.gms | 15 +- .../80_optimization/testOneRegi/bounds.gms | 2 +- .../80_optimization/testOneRegi/datainput.gms | 6 +- .../testOneRegi/declarations.gms | 22 +- .../80_optimization/testOneRegi/equations.gms | 36 +- .../80_optimization/testOneRegi/postsolve.gms | 14 +- .../80_optimization/testOneRegi/preloop.gms | 74 +-- .../testOneRegi/realization.gms | 2 +- 18 files changed, 434 insertions(+), 475 deletions(-) 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/declarations.gms b/core/declarations.gms index 36ced4acbc..e1f2e98a22 100644 --- a/core/declarations.gms +++ b/core/declarations.gms @@ -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/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/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/datainput.gms b/modules/80_optimization/nash/datainput.gms index 49d53afae0..9cea72df23 100644 --- a/modules/80_optimization/nash/datainput.gms +++ b/modules/80_optimization/nash/datainput.gms @@ -6,22 +6,20 @@ *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/datainput.gms -*** Negishi weights, not used in Nash -pm_w(regi) = 1; *** convergence with trade surplus thresholds if(cm_nash_autoconverge > 0, cm_iteration_max = 100; !! set max number of iterations -*** default values for cm_nash_autoconverge = 1 - coarse - 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") = 350 * 12/44 / 1000; !! 350 MtCO2eq, converted 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)), -*** convergences thresholds - fine +*** fine convergence thresholds (cm_nash_autoconverge = 2) p80_surplusMaxTolerance(trade) $ (cm_nash_autoconverge = 2) = 0.2 * p80_surplusMaxTolerance(trade); -*** convergences thresholds - very coarse +*** coarse convergence thresholds (cm_nash_autoconverge = 3) p80_surplusMaxTolerance(trade) $ (cm_nash_autoconverge = 3) = 2 * p80_surplusMaxTolerance(trade); ); ); @@ -30,35 +28,34 @@ if(cm_nash_autoconverge > 0, *** Nash adjustment costs (default value around 150). *** Multiplicator to penalise trade patterns deviations from the last iteration. Involves a trade-off: -*** - if too low, markets jump far away from clearance. -*** - if too high, changes in trade patten over iterations are very slow, convergence takes many many iterations. -p80_etaAdj(tradePe) = 80; -p80_etaAdj("good") = 100; -p80_etaAdj("perm") = 10; +*** - 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. *** Multiplicator to price anticipation -p80_etaXp(tradePe) = 0.1; -p80_etaXp("good") = 0.1; -p80_etaXp("perm") = 0.2; +p80_priceAnticipStrength(tradePe) = 0.1; +p80_priceAnticipStrength("good") = 0.1; +p80_priceAnticipStrength("perm") = 0.2; -*** LB: parameter for Nash price algorithm between different iterations. These parameters are pretty sensitive: -*** - if market surpluses diverge, try higher values (up to 1). -*** - if surpluses oscillate, try lower values. +*** 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 *** short term price ajustment elasticity -p80_etaST(tradePe) = 0.3; -p80_etaST("pebiolc") = 0.8; !! AJS: bio market seems to like this -p80_etaST("peur") = 0.2; !! uranium market is more sensitive, so choose lower etaST -p80_etaST("good") = 0.25; -p80_etaST("perm") = 0.3; -$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. +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. *** long term price ajustment elasticity -p80_etaLT(trade) = 0; -p80_etaLT("perm") = 0.03; +p80_longTermFac(trade) = 0; +p80_longTermFac("perm") = 0.03; *** --------------- Initialise convergence process parameters --------------- @@ -67,6 +64,9 @@ s80_converged = 0; s80_fadeoutPriceAnticipStartingPeriod = 0; sm_fadeoutPriceAnticip = 1; +*** Negishi weights, not used in Nash +pm_w(regi) = 1; + *** parameter for spillover externality (aggregated productivity level) pm_cumEff(t, regi, in) = 100; @@ -101,6 +101,7 @@ $include "./modules/80_optimization/nash/input/prices_NASH.inc"; *** --------------- 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" diff --git a/modules/80_optimization/nash/declarations.gms b/modules/80_optimization/nash/declarations.gms index 75c8b27e03..b251522399 100644 --- a/modules/80_optimization/nash/declarations.gms +++ b/modules/80_optimization/nash/declarations.gms @@ -7,37 +7,38 @@ *** 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" -*** 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" -p80_etaAdj(all_enty) "Adjustment costs for changes of trade pattern between iterations" +*** 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" *** 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." +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." -p80_marketVolume(ttot,all_regi,all_enty) "Normalization parameter for market volume [amount of trade item]" -p80_intertemporalSurplusRevenue(all_enty) "Aggregated intertemporal value of the surplus [T$2017]" -p80_itertemporalMarketRevenue(all_enty) "Aggregated intertemporal market volume [T$2017]" +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]" -*** 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 [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 [%]" +*** 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 [%]" -*** 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" +*** 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" -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_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)" p80_surplusMax_iter(all_enty,iteration,tall) "Diagnostics for Nash: Worst residual market surplus until given year, absolute value [TWa, T$, GtC]" @@ -106,7 +107,7 @@ 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. +*** 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." ; @@ -118,7 +119,7 @@ q80_budgetPermRestr(all_regi) "constraints regional permit budget to given regi scalars *** 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) +*** 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 7f2465fd80..d07c6e9f17 100644 --- a/modules/80_optimization/nash/equations.gms +++ b/modules/80_optimization/nash/equations.gms @@ -8,40 +8,44 @@ *' @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") !! 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) !! only active in 24_trade capacity mode + + 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_etaXp(trade) + + 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_marketVolume(ttot,regi,trade) + sm_eps) + / p80_regiMarketVolume(ttot,regi,trade) ) ) ) ); -*' quadratic adjustment costs, penalizing deviations from the trade pattern of the last iteration. +*' 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_etaAdj(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_marketVolume(ttot,regi,trade) + sm_eps)); + * 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)); *' 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)); +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 53a5b127ce..7208f73f22 100644 --- a/modules/80_optimization/nash/postsolve.gms +++ b/modules/80_optimization/nash/postsolve.gms @@ -12,24 +12,21 @@ *' #### Price corrections to improve convergence of next iteration ***------------------------------------------------------------------------------ -*' The objective of the nash optimisation is to find a set of prices that clears the markets for all trade items. +*' 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_etaST_correct = p80_etaST * p80_surplus / p80_marketVolume +*' for example p80_shortTermCorrect = p80_shortTermFac * p80_surplus / p80_regiMarketVolume -*' Compute market volume for different trades (take values from last iteration for regions that were not solved optimally) +*' Compute market volume for different trades loop(ttot $ (ttot.val >= 2005), - loop(regi $ (pm_SolNonInfes(regi) = 1), - p80_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); - p80_marketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; +*** 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; ); - loop(regi $ (pm_SolNonInfes(regi) = 0), - p80_marketVolume(ttot,regi,"good") = p80_marketVolume(ttot,regi,"good"); - p80_marketVolume(ttot,regi,tradePe) = p80_marketVolume(ttot,regi,tradePe); - ); -*** ML: normalize permit trade corrections to consumption or positive cap path instead of emissions, as those may be negative - p80_marketVolume(ttot,regi,"perm") = abs(pm_shPerm(ttot,regi) * pm_emicapglob("2050")); - p80_marketVolume(ttot,regi,trade) = max(sm_eps, p80_marketVolume(ttot,regi,trade)) !! ensure market volume is positive + 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 ); loop(trade $ (not tradeSe(trade)), @@ -39,35 +36,35 @@ loop(trade $ (not tradeSe(trade)), + (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_marketVolume) - p80_intertemporalSurplusRevenue(trade) = sum(ttot $ (ttot.val >= cm_startyear), +*' 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_itertemporalMarketRevenue(trade) = max(sm_eps, sum((ttot, regi) $ (ttot.val >= 2005), - pm_pvp(ttot,trade) * pm_ts(ttot) * p80_marketVolume(ttot,regi,trade))); + 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_etaLT_correct(trade,iteration) = - p80_etaLT(trade) * p80_intertemporalSurplusRevenue(trade) / p80_itertemporalMarketRevenue(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_etaST_correct(ttot,trade,iteration) $ (ttot.val >= 2005) = - p80_etaST(trade) * p80_surplus(ttot,trade,iteration) / max(sm_eps , sum(regi, p80_marketVolume(ttot,regi,trade))); + p80_shortTermCorrect(ttot,trade,iteration) $ (ttot.val >= 2005) = + p80_shortTermFac(trade) * p80_surplus(ttot,trade,iteration) / p80_globMarketVolume(ttot,trade); - p80_etaST_correct(ttot,"perm",iteration) $ (ttot.val >= 2005) = p80_etaST_correct(ttot,"perm",iteration) + 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_etaST_correct(ttot,tradePe(trade),iteration) $ (ttot.val >= 2005) = p80_etaST_correct(ttot,trade,iteration) + 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 ); *' If the surplus remains over several iterations, increase the price correction terms -p80_etaST_correct_safecopy(ttot,trade,iteration) $ (not tradeSe(trade)) = p80_etaST_correct(ttot,trade,iteration); !! copy of initial values - 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 + 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 @@ -79,17 +76,11 @@ loop(ttot $ (ttot.val >= 2005), ); ); - if(iteration.val > 15 and o80_trackSurplusSign(ttot,trade,iteration) >= 5, !! if surplus was beyond tolerance 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 > 20 and o80_trackSurplusSign(ttot,trade,iteration) >= 10, !! push stronger if previous increase did not help after a few 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 > 25 and o80_trackSurplusSign(ttot,trade,iteration) >= 15, !! push stronger if previous increase did not help after a few iterations - p80_etaST_correct(ttot,trade,iteration) = 2 * p80_etaST_correct(ttot,trade,iteration); - o80_counter_iteration_trade_ttot(ttot,trade,iteration) = 3; +*** 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 @@ -98,7 +89,7 @@ loop(ttot $ (ttot.val >= 2005), *** 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_etaLT_correct(trade,iteration) - p80_etaST_correct(ttot,trade,iteration)); + 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)), @@ -114,8 +105,8 @@ p80_taxrev0(ttot,regi) $ (ttot.val >= max(2010, cm_startyear) and pm_SolNonInfes 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)) / sum(regi, p80_marketVolume(ttot2,regi,trade))); + 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"); @@ -131,7 +122,7 @@ loop(trade $ (not tradeSe(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_itertemporalMarketRevenue("good") / pm_pvp("2005","good")); +p80_defic_sum_rel(iteration) = 100 * p80_defic_sum(iteration) / (p80_itertempMarketRevenue("good") / pm_pvp("2005","good")); *** adjust parameters for next iteration @@ -170,11 +161,11 @@ pm_FEPrice_iter(iteration,t,regi,enty,sector,emiMkt) = pm_FEPrice(t,regi,enty,se loop(ttot $ (ttot.val >= 2005), loop(trade $ (not tradeSe(trade)), p80_PriceChangePriceAnticipReg(ttot,trade,regi) = 100 * - sm_fadeoutPriceAnticip * p80_etaXp(trade) + 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_marketVolume(ttot,regi,trade) + sm_eps); + / p80_regiMarketVolume(ttot,regi,trade); p80_DevPriceAnticipReg(ttot,trade,regi) = ( vm_Xport.l(ttot,regi,trade) - vm_Mport.l(ttot,regi,trade) ) @@ -274,7 +265,7 @@ loop(regi, ); ); !! loop over regi -*' criterion: are the anticipation terms sufficienctly small? (only for checking, not applied anymore) +*' criterion: are the anticipation terms sufficienctly small? (only for checking, not enforced anymore) p80_fadeoutPriceAnticip_iter(iteration) = sm_fadeoutPriceAnticip; if(sm_fadeoutPriceAnticip > cm_maxFadeOutPriceAnticip, !! s80_bool = 0; !! not an active convergence criterion anymore @@ -401,218 +392,225 @@ 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, "target"), - 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)"; - ); +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, "target"), + 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 = 0) and (iteration.val = 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, "target"), - 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; + ); +$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 "################################################################################################"; ); @@ -621,36 +619,33 @@ $endIf.cm_implicitPePriceTarget *' #### Finishing or aborting ***------------------------------------------------------------------------------ -***if all conditions are met, stop optimization. +*** if all conditions are met, stop optimization. if(s80_bool = 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 - + 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 @@ -666,7 +661,7 @@ if (cm_abortOnConsecFail > 0, ); if (smax(regi, p80_trackConsecFail(regi)) >= cm_abortOnConsecFail, - if ((s80_runInDebug = 0) and (cm_nash_mode ne 1), !! auto-start debug only if not already in debug mode + 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."; @@ -679,8 +674,7 @@ if (cm_abortOnConsecFail > 0, 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 + 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; @@ -744,9 +738,7 @@ 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 diff --git a/modules/80_optimization/nash/preloop.gms b/modules/80_optimization/nash/preloop.gms index a1fe4f538c..436d2427fe 100644 --- a/modules/80_optimization/nash/preloop.gms +++ b/modules/80_optimization/nash/preloop.gms @@ -6,7 +6,11 @@ *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/preloop.gms -*** MLB/AG: for Nash algorithm read initial price data from gdx +***------------------------------------------------------------------------------ +*** 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; @@ -15,43 +19,31 @@ 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; - -*** 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 -*** using execute_loadpoint in core/preloop.gms. The price paths are the trickiest part. Try to find p80_priceXXX prices in the gdx first, -*** if that fails, fallback to price path read from input/prices_NASH.inc loop(ttot $ (ttot.val >= 2005), loop(trade $ (not tradeSe(trade)), - if(pm_pvp(ttot,trade) = NA or pm_pvp(ttot,trade) < 1e-12 or pm_pvp(ttot,trade) > 0.1, !! in case price paths from the gdx are not valid +*** 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. '; ); - if(pm_pvp(ttot,trade) = NA, pm_pvp(ttot,trade) = 0;); !! in case pm_pvp is not found in gdx +*** 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; - loop(regi, - pm_Xport0(ttot,regi,trade) = vm_Xport.l(ttot,regi,trade); - p80_Mport0(ttot,regi,trade) = vm_Mport.l(ttot,regi,trade); - - if(pm_Xport0(ttot,regi,trade) = NA, !! in case export is not found in gdx - pm_Xport0(ttot,regi,trade) = 0; - vm_Xport.L(ttot,regi,trade) = 0; - ); - if(p80_Mport0(ttot,regi,trade) = NA, !! in case import is not found in gdx - p80_Mport0(ttot,regi,trade) = 0; - vm_Mport.L(ttot,regi,trade) = 0; - ); - p80_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); - p80_marketVolume(ttot,regi,"perm") $ (ttot.val >= 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1e-6); - p80_marketVolume(ttot,regi,tradePe) = (sum(rlf, vm_fuExtr.l(ttot,regi,tradePe,rlf)) + vm_prodPe.l(ttot,regi,tradePe)) / 2; +*** 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 - p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); - ); + p80_taxrev0(ttot,regi) = vm_taxrev.l(ttot,regi); ); ); -p80_Mport0("2005",regi,tradePe) $ (p80_Mport0("2005",regi,tradePe) = NA) = 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, @@ -60,23 +52,25 @@ if(cm_emiscen ne 1 and cm_emiscen ne 9 and smax(t, pm_pvp(t,"perm")) = 0, 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(ttot $ (not sameas(ttot,"2005")), - p80_pvp_itr(ttot,tradePe,"1") $ (p80_pvp_itr(ttot,tradePe,"1") = 0) = p80_pvp_itr(ttot-1,tradePe,"1"); +*** 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); ); +*** 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_marketVolume; +display pm_pvp,p80_regiMarketVolume; display pm_Xport0,p80_Mport0; display p80_surplusMaxTolerance; +*** 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 = 6, diff --git a/modules/80_optimization/nash/presolve.gms b/modules/80_optimization/nash/presolve.gms index 1719b3466f..cb427b69d3 100644 --- a/modules/80_optimization/nash/presolve.gms +++ b/modules/80_optimization/nash/presolve.gms @@ -6,22 +6,16 @@ *** | Contact: remind@pik-potsdam.de *** SOF ./modules/80_optimization/nash/presolve.gms -*** 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. +*** 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: take care of the case of infeasible solution + pm_cesdata("2005",regi2,in,"eff") * vm_effGr.l(ttot,regi2,in)); !! TODO: case where regi2 is infeasible - pm_co2eqForeign(ttot,regi) - = sum(regi2 $ (not sameas(regi,regi2)), pm_co2eq0(ttot,regi2)); !! does this interfere with the initialization in datainput? - - 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)); + 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. diff --git a/modules/80_optimization/nash/solve.gms b/modules/80_optimization/nash/solve.gms index cb1d6b07eb..7cd433347a 100644 --- a/modules/80_optimization/nash/solve.gms +++ b/modules/80_optimization/nash/solve.gms @@ -18,11 +18,8 @@ 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 -$ifthen.repeatNonOpt "%cm_repeatNonOpt%" == "off" - or p80_repy(all_regi,"modelstat") = 7 -$endif.repeatNonOpt +$if "%cm_repeatNonOpt%" == "off" or p80_repy(all_regi,"modelstat") = 7 ), - p80_repy_thisSolitr(all_regi,solveinfo80) = 0; continue; ); @@ -66,7 +63,7 @@ $endif.repeatNonOpt 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; @@ -90,9 +87,8 @@ 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, 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 ; @@ -109,11 +105,10 @@ loop(regi, ); *** 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 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..3d87d9f655 100644 --- a/modules/80_optimization/testOneRegi/datainput.gms +++ b/modules/80_optimization/testOneRegi/datainput.gms @@ -10,9 +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; +p80_priceAnticipStrength(tradePe) = 1; +p80_priceAnticipStrength("good") = 1; +p80_priceAnticipStrength("perm") = 1; *load fallback prices diff --git a/modules/80_optimization/testOneRegi/declarations.gms b/modules/80_optimization/testOneRegi/declarations.gms index 368278d692..24d1a88ec5 100644 --- a/modules/80_optimization/testOneRegi/declarations.gms +++ b/modules/80_optimization/testOneRegi/declarations.gms @@ -7,23 +7,23 @@ *** 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_marketVolume(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_pvpFallback(ttot,all_enty) "Helper parameter. Price path from input/prices_NASH.inc. Only used if reading prices from gdx fails." +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 7f9aec5223..196bf408e1 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_marketVolume(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) + ) ) - ); + ) + ); -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/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 ba3fac09af..416ac23f58 100644 --- a/modules/80_optimization/testOneRegi/preloop.gms +++ b/modules/80_optimization/testOneRegi/preloop.gms @@ -5,55 +5,33 @@ *** | 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; - -*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_marketVolume(ttot,regi,"good") = vm_cons.l(ttot,regi); - p80_marketVolume(ttot,regi,"perm")$(ttot.val ge 2005) = max(abs(pm_shPerm(ttot,regi) * pm_emicapglob(ttot)) , 1E-6); - p80_marketVolume(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); - ); +***------------------------------------------------------------------------------ +*** 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); ); ); 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. From 5f90bcfb250bf439f839b30e0927b1efaf937043 Mon Sep 17 00:00:00 2001 From: lecfab Date: Tue, 23 Dec 2025 15:06:46 +0100 Subject: [PATCH 4/8] capitalising s_ZJ_2_TWa and making precision consistent with other conversion factors --- core/datainput.gms | 21 +++-- core/declarations.gms | 2 +- core/loop.gms | 192 ++++++++++++++++++------------------------ 3 files changed, 93 insertions(+), 122 deletions(-) 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 e1f2e98a22..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/, 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"; ); From 9ec08db9cc1d163f83e6247b5ec8acaaecbece63 Mon Sep 17 00:00:00 2001 From: lecfab Date: Tue, 23 Dec 2025 15:07:53 +0100 Subject: [PATCH 5/8] avoid division by zero although p80_regiMarketVolume is never supposed to be zero --- modules/80_optimization/nash/equations.gms | 5 +++-- modules/80_optimization/testOneRegi/equations.gms | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/modules/80_optimization/nash/equations.gms b/modules/80_optimization/nash/equations.gms index d07c6e9f17..d92c0e1063 100644 --- a/modules/80_optimization/nash/equations.gms +++ b/modules/80_optimization/nash/equations.gms @@ -26,7 +26,7 @@ q80_budg_intertemp(regi).. + 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) + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) ) ) ) @@ -39,7 +39,8 @@ q80_costAdjNash(ttot,regi) $ (ttot.val >= cm_startyear) .. 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)); + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) + ); *' link between permit budget and emission budget q80_budgetPermRestr(regi) $ (cm_emiscen = 6) .. diff --git a/modules/80_optimization/testOneRegi/equations.gms b/modules/80_optimization/testOneRegi/equations.gms index 196bf408e1..be9e1b8b60 100644 --- a/modules/80_optimization/testOneRegi/equations.gms +++ b/modules/80_optimization/testOneRegi/equations.gms @@ -20,7 +20,7 @@ q80_budg_intertemp(regi).. (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) + / (p80_regiMarketVolume(ttot,regi,trade) + sm_eps) ) ) ) From fb94550f96e8d18f7d5ba66540e8379279c282b7 Mon Sep 17 00:00:00 2001 From: lecfab Date: Tue, 23 Dec 2025 15:08:50 +0100 Subject: [PATCH 6/8] nash refactor: fixing compile bugs --- modules/80_optimization/nash/postsolve.gms | 8 ++++---- modules/80_optimization/nash/preloop.gms | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/modules/80_optimization/nash/postsolve.gms b/modules/80_optimization/nash/postsolve.gms index 7208f73f22..5779740ccf 100644 --- a/modules/80_optimization/nash/postsolve.gms +++ b/modules/80_optimization/nash/postsolve.gms @@ -279,8 +279,8 @@ if(p80_DevPriceAnticipGlobAllMax2100Iter(iteration) > 0.1 * p80_surplusMaxTolera 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 <= 17, +*' 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; ); @@ -349,7 +349,7 @@ $endIf.cm_implicitPePriceTarget p80_globalBudget_absDev_iter(iteration) = sm_globalBudget_absDev; if (abs(p80_globalBudget_absDev_iter(iteration)) > cm_budgetCO2_absDevTol, !! check if CO2 budget met in tolerance range, s80_bool = 0; - p80_messageShow("target") = YES; + p80_messageShow("globalbudget") = YES; ); *' criterion: is the peak budget year the same as the year of maximum cumulative CO2 emissions? @@ -430,7 +430,7 @@ loop(convMessage80$(p80_messageShow(convMessage80)), display "#### Check out sm_fadeoutPriceAnticip which needs to be below cm_maxFadeOutPriceAnticip."; display sm_fadeoutPriceAnticip, cm_maxFadeOutPriceAnticip; ); - if(sameas(convMessage80, "target"), + 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"; diff --git a/modules/80_optimization/nash/preloop.gms b/modules/80_optimization/nash/preloop.gms index 436d2427fe..6a665a8211 100644 --- a/modules/80_optimization/nash/preloop.gms +++ b/modules/80_optimization/nash/preloop.gms @@ -39,7 +39,7 @@ loop(ttot $ (ttot.val >= 2005), 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_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); ); From 15c1a8e9c426b516a91d1088bad312dfe61c4d2e Mon Sep 17 00:00:00 2001 From: lecfab Date: Tue, 23 Dec 2025 15:10:32 +0100 Subject: [PATCH 7/8] nash refactor: porting modifications to testOneRegi --- modules/80_optimization/testOneRegi/datainput.gms | 4 ---- modules/80_optimization/testOneRegi/declarations.gms | 1 - modules/80_optimization/testOneRegi/not_used.txt | 3 ++- modules/80_optimization/testOneRegi/preloop.gms | 7 ++++++- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/modules/80_optimization/testOneRegi/datainput.gms b/modules/80_optimization/testOneRegi/datainput.gms index 3d87d9f655..c32d2dcee5 100644 --- a/modules/80_optimization/testOneRegi/datainput.gms +++ b/modules/80_optimization/testOneRegi/datainput.gms @@ -14,10 +14,6 @@ p80_priceAnticipStrength(tradePe) = 1; p80_priceAnticipStrength("good") = 1; p80_priceAnticipStrength("perm") = 1; -*load fallback prices - -$include "./modules/80_optimization/testOneRegi/input/prices_NASH.inc"; - *MLB 12062013* initialize learning externality (can be improved by including file) pm_capCumForeign(ttot,regi,teLearn)$(ttot.val ge 2005) = 0; diff --git a/modules/80_optimization/testOneRegi/declarations.gms b/modules/80_optimization/testOneRegi/declarations.gms index 24d1a88ec5..412364d039 100644 --- a/modules/80_optimization/testOneRegi/declarations.gms +++ b/modules/80_optimization/testOneRegi/declarations.gms @@ -8,7 +8,6 @@ parameter p80_priceAnticipStrength(all_enty) "Parameter governing price anticipation on commodity markets" -p80_pvpFallback(ttot,all_enty) "Helper parameter. Price path from input/prices_NASH.inc. Only used if reading prices from gdx fails." 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" 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/preloop.gms b/modules/80_optimization/testOneRegi/preloop.gms index 416ac23f58..16d45d54f3 100644 --- a/modules/80_optimization/testOneRegi/preloop.gms +++ b/modules/80_optimization/testOneRegi/preloop.gms @@ -29,12 +29,17 @@ loop(ttot $ (ttot.val >= 2005), 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_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: 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"; display pm_pvp; From c9dc023940e1b3c970a1ca0f20ddb14953fcf2a7 Mon Sep 17 00:00:00 2001 From: lecfab Date: Tue, 23 Dec 2025 15:11:47 +0100 Subject: [PATCH 8/8] nash refactor: removing unneeded file --- modules/80_optimization/negishi/not_used.txt | 2 + .../testOneRegi/input/prices_NASH.inc | 135 ------------------ 2 files changed, 2 insertions(+), 135 deletions(-) delete mode 100644 modules/80_optimization/testOneRegi/input/prices_NASH.inc 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/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