improved the robustness of the peakBudget-finding algorithm for iterative...

improved the robustness of the peakBudget-finding algorithm for iterative target_adj = 9 and of diffCurvPhaseIn2Lin. Also updated the description in default.cfg

config/default.cfg

@@ -206,7 +206,10 @@ cfg$gms$banking  <- "off"           # def = off
 # * (ExogSameAsPrevious): uses the identical carbon prices of a previous run, so far requires manual updating of the input file - see outcommented text in scripts/start_functions.R
 # * (temperatureNotToExceed): [test and verify before using it!] Find the optimal carbon carbon tax (set cm_emiscen=9) for the not-to-exceed temperature target defined by cm_carbonprice_temperatureLimit.
 # * (NDC2constant): linearly phase in global constant price from NDC prices (default 2020-2040 phase-in)
-# * (diffPhaseIn2Constant): linearly phase in global constant price, with starting values differentiated by GDP/cap (Robert's proposal)
+# * (diffCurvPhaseIn2Lin): [REMIND 2.1 default for validation peakBudget runs, in combination with "iterative_target_adj" = 9] curved convergence of CO2 prices between regions until cm_CO2priceRegConvEndYr; developed countries have linear path from 0 in 2010 through cm_co2_tax_2020 in 2020; 
+# * (diffPhaseIn2Constant): !experimental! linearly phase in global constant price, with starting values differentiated by GDP/cap 
+# * 
+
 
 cfg$gms$carbonprice <- "none"           # def = none
 
@@ -660,6 +663,7 @@ cfg$RunsUsingTHISgdxAsBAU <- NA
 # (5): iterative adjustment V for both budget and tax runs based on CO2 emissions 2011-2100, see core/postsolve.gms for direct algorithms of adjustment
 # (6): iterative adjustment VI for both budget and tax runs based on peak CO2 emissions budget, without changing temporal profile (i.e. with overshoot), see core/postsolve.gms for direct algorithms of adjustment
 # (7): iterative adjustment VII for tax runs based on peak CO2 emissions, with change of temporal tax profile after time of peak budget, aiming for net-zero thereafter, see core/postsolve.gms for direct algorithms of adjustment
+# (9): [REMIND 2.1 default for validation peakBudget runs, in combination with carbonprice = diffCurvPhaseIn2Lin) ] aims at achieving c_budgetCO2 as peakBudget. Scales up and down the original CO2 price path before the peaking year; after the peaking year annual increase by cm_taxCO2inc_after_peakBudgYr. Automatically shifts cm_peakBudgYr to find the correct peaking year for a given .
 #c_budgetCO2 "cumulative CO2|Total emissions budget 2011-2100, Gt CO2"
 # budgets from 1.5°C report, chapter 2, Table 2.2: 67% 1.5°C GSAT budget 2018-2100 = 420 Gt CO2, 67% 2.0°C GSAT budget 2018-2100 = 1170
 # add emissions 2011-2018 = 290 Gt CO2, subtract emissions from feedbacks (e.g. thawing permafrost) 100 Gt CO2

core/declarations.gms

@@ -53,14 +53,16 @@ o_diff_to_Budg(iteration)                             "Difference between actual
 o_totCO2emi_peakBudgYr(iteration)                     "Total CO2 emissions in the peakBudgYr"
 o_peakBudgYr_Itr(iteration)                           "Year in which the CO2 budget is supposed to peak. Is changed in iterative_target_adjust = 9"
 o_factorRescale_taxCO2_afterPeakBudgYr(iteration)     "Multiplicative factor for rescaling the CO2 price in the year after peakBudgYr - only needed if flip-flopping of peakBudgYr occurs"
-o_delay_increase_peakBudgYear(iteration)              "Switch that tracks if flip-flopping of peakBudgYr happened. Starts an inner loop to try and overcome this"
-o_reached_until2150pricepath(iteration)               "Switch that tracks if the inner loop of increasing the CO2 price AFTER peakBudgYr goes beyond the initial trajectory"
-p_taxCO2eq_until2150(ttot,all_regi)                  "CO2 price trajectory continued until 2150 - as if there was no change in trajectory after peakBudgYr. Needed to recalculate CO2 price trajectory after peakBudgYr was shifted right"
+o_delay_increase_peakBudgYear(iteration)              "Counter that tracks if flip-flopping of peakBudgYr happened. Starts an inner loop to try and overcome this"
+o_reached_until2150pricepath(iteration)               "Counter that tracks if the inner loop of increasing the CO2 price AFTER peakBudgYr goes beyond the initial trajectory"
+p_taxCO2eq_until2150(ttot,all_regi)                   "CO2 price trajectory continued until 2150 - as if there was no change in trajectory after peakBudgYr. Needed to recalculate CO2 price trajectory after peakBudgYr was shifted right"
 o_totCO2emi_allYrs(ttot,iteration)                    "Global CO2 emissions over time and iterations. Needed to check the procedure to find the peakBudgYr"
 o_change_totCO2emi_peakBudgYr                         "Measure for how much the CO2 emissions change around the peakBudgYr"
 p_factorRescale_taxCO2(iteration)                     "Multiplicative factor for rescaling the CO2 price to reach the target"
 p_factorRescale_taxCO2_Funneled(iteration)            "Multiplicative factor for rescaling the CO2 price to reach the target - limited by an iteration-dependent funnel"
 o_taxCO2eq_Itr_1regi(ttot,iteration)                  "CO2 taxed in the last region, tracked over iterations for debugging" 
+o_pkBudgYr_flipflop(iteration)                        "Counter that tracks if flipfloping of cm_peakBudgYr occured in the last iterations"
+o_taxCO2eq_afterPeakShiftLoop_Itr_1regi(ttot, iteration) "CO2 taxed in the last region, after the loop that shifts peakBudgYr, tracked over iterations for debugging"
 
 ***----------------------------------------------------------------------------------------
 ***-----------------------------------------------ESM module-------------------------------

core/postsolve.gms

@@ -295,8 +295,7 @@ display p_actualbudgetco2;
 *** if tax after cm_peakBudgYr is higher than normal increase rate (exceeding a 20% tolerance): shift right
 		elseif( ( sum(regi, sum(t2$(t2.val = pm_ttot_val(ttot+1)),pm_taxCO2eq(t2,regi))) > sum(regi,sum(t2$(t2.val = pm_ttot_val(ttot+1)),p_taxCO2eq_until2150(t2,regi)))*1.2 ) AND (cm_peakBudgYr < 2100) ), !! if peaking time would be after 2100, keep 2100 budget year
           if(  (iteration.val > 2) AND ( o_peakBudgYr_Itr(iteration - 1) > o_peakBudgYr_Itr(iteration) ) AND ( o_peakBudgYr_Itr(iteration - 2) = o_peakBudgYr_Itr(iteration) ) , !! if the target year was just shifted left after being shifted right
-            o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration); !! don't shift right again immediately, but go into a different loop:
-            o_delay_increase_peakBudgYear(iteration) = 1;
+            o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration); !! don't shift right again immediately
           else
 		    display "shift peakBudgYr right";
             o_peakBudgYr_Itr(iteration+1) =  pm_ttot_val(ttot + 1);  !! ttot+1 is the new peakBudgYr
@@ -422,6 +421,13 @@ if(cm_iterative_target_adj eq 9,
 *** ----B1: check if cm_peakBudgYr should be shifted left or right: 
     if( abs(o_diff_to_Budg(iteration)) < 20,                      !! only think about shifting peakBudgYr if the budget is close enough to target budget
       display "close enough to target budget to check timing of peak year";
+	  
+	  !!  check if the target year was just shifted back left after being shifted right before
+	  if ( (iteration.val > 2) AND ( o_peakBudgYr_Itr(iteration - 1) > o_peakBudgYr_Itr(iteration) ) AND ( o_peakBudgYr_Itr(iteration - 2) = o_peakBudgYr_Itr(iteration) ),
+	    o_pkBudgYr_flipflop(iteration) = 1; 
+        display "flipflop observed (before loop)";
+	  );
+	  
       loop(ttot$(ttot.val = cm_peakBudgYr),                               !! look at the peak timing
         if(  ( (o_totCO2emi_peakBudgYr(iteration) < -(0.1 + o_change_totCO2emi_peakBudgYr(iteration)) ) AND (cm_peakBudgYr > 2040) ), !! no peaking time before 2040
           display "shift peakBudgYr left";
@@ -429,9 +435,13 @@ if(cm_iterative_target_adj eq 9,
           pm_taxCO2eq(t,regi)$(t.val gt pm_ttot_val(ttot - 1)) = p_taxCO2eq_until2150(ttot-1,regi) + (t.val - pm_ttot_val(ttot - 1)) * cm_taxCO2inc_after_peakBudgYr * sm_DptCO2_2_TDpGtC;  !! increase by cm_taxCO2inc_after_peakBudgYr per year after peakBudgYr
         
 		elseif ( ( o_totCO2emi_peakBudgYr(iteration) > (0.1 + o_change_totCO2emi_peakBudgYr(iteration)) ) AND (cm_peakBudgYr < 2100) ), !! if peaking time would be after 2100, keep 2100 budget year
-          if(  (iteration.val > 2) AND ( o_peakBudgYr_Itr(iteration - 1) > o_peakBudgYr_Itr(iteration) ) AND ( o_peakBudgYr_Itr(iteration - 2) = o_peakBudgYr_Itr(iteration) ) , !! if the target year was just shifted left after being shifted right
+          if(  (o_pkBudgYr_flipflop(iteration) eq 1), !! if the target year was just shifted left after being shifted right, and would now be shifted right again
+            display "peakBudgYr was left, right, left and is now supposed to be shifted right again -> flipflop, thus go into separate loop";
             o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration); !! don't shift right again immediately, but go into a different loop:
             o_delay_increase_peakBudgYear(iteration) = 1;
+		  elseif ( o_delay_increase_peakBudgYear(iteration) eq 1 ),
+		    display "still in separate loop trying to resolve flip-flop behavior";
+			o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration); !! keep current peakBudgYr,
           else
 		    display "shift peakBudgYr right";
             o_peakBudgYr_Itr(iteration+1) =  pm_ttot_val(ttot + 1);  !! ttot+1 is the new peakBudgYr
@@ -456,27 +466,36 @@ if(cm_iterative_target_adj eq 9,
       display "not shifting peakBudgYr right, instead adjusting CO2 price for following year";
       loop(ttot$(ttot.val eq cm_peakBudgYr),  !! set ttot to the current peakBudgYr 
         loop(t2$(t2.val eq pm_ttot_val(ttot+1)),  !! set t2 to the following time step
-          o_factorRescale_taxCO2_afterPeakBudgYr(iteration) = 1 + max(sum(regi2,vm_emiAll.l(t2,regi2,"co2"))/sum(regi2,vm_emiAll.l("2015",regi2,"co2")),-0.75) ; !! inspired by Christoph. This value is 1 if emissions are 0.
+          o_factorRescale_taxCO2_afterPeakBudgYr(iteration) = 1 + max(sum(regi2,vm_emiAll.l(ttot,regi2,"co2"))/sum(regi2,vm_emiAll.l("2015",regi2,"co2")),-0.75) ; 
+		  !! this was inspired by Christoph's approach. This value is 1 if emissions in the peakBudgYr are 0; goes down to 0.25 if emissions are <0 and approaching the size of 2015 emissions, and > 1 if emissions > 0. 
           
 		  !! in case the normal linear extension still is not enough to get emissions to 0 after the peakBudgYr, shift peakBudgYr right again:
-          if( ( o_reached_until2150pricepath(iteration-1) eq 1 ) AND ( (o_factorRescale_taxCO2_afterPeakBudgYr(iteration) / p_factorRescale_taxCO2_Funneled(iteration)) > 1), 
-            display "price in following year reached original path and is still not enough -> shift peakBudgYr to right";
-            o_delay_increase_peakBudgYear(iteration) = 0;
+          if( ( o_reached_until2150pricepath(iteration-1) eq 1 ) AND ( o_totCO2emi_peakBudgYr(iteration) > (0.1 + o_change_totCO2emi_peakBudgYr(iteration)) ), 
+            display "price in following year reached original path in previous iteration and is still not enough -> shift peakBudgYr to right";
+            o_delay_increase_peakBudgYear(iteration+1) = 0;  !! probably is not necessary
             o_reached_until2150pricepath(iteration) = 0;
             o_peakBudgYr_Itr(iteration+1) = t2.val;        !! shift PeakBudgYear to the following time step
             cm_peakBudgYr = o_peakBudgYr_Itr(iteration+1);
             pm_taxCO2eq(t2,regi) = p_taxCO2eq_until2150(t2,regi) ;  !! set CO2 price in t2 to value in the "continuous path"
       
             display cm_peakBudgYr;
+		  elseif ( ( o_reached_until2150pricepath(iteration-1) eq 1 ) AND ( o_totCO2emi_peakBudgYr(iteration) < (0.1 + o_change_totCO2emi_peakBudgYr(iteration)) ) ), 
+            display "New intermediate price in timestep after cm_peakBudgYr is sufficient to stabilize peaking year - go back to normal loop";	
+			o_delay_increase_peakBudgYear(iteration+1) = 0;  !! probably is not necessary
+            o_reached_until2150pricepath(iteration) = 0;
+			o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration);  
+            cm_peakBudgYr = o_peakBudgYr_Itr(iteration+1);
           else      !! either didn't reach the continued "until2150"-price path in last iteration, or the increase was high enough to get emissions to 0. 
 		            !! in this case, keep PeakBudgYr, and adjust the price in the year after the peakBudgYr to get emissions close to 0,
+			o_delay_increase_peakBudgYear(iteration+1) = 1; !! make sure next iteration peakBudgYr is not shifted right again
 		    o_peakBudgYr_Itr(iteration+1) = o_peakBudgYr_Itr(iteration);
             pm_taxCO2eq(t2,regi) = max(pm_taxCO2eq(ttot,regi), !! at least as high as the price in the peakBudgYr
                                        pm_taxCO2eq(t2,regi) * (o_factorRescale_taxCO2_afterPeakBudgYr(iteration) / p_factorRescale_taxCO2_Funneled(iteration) ) !! the full path was already rescaled by p_factorRescale_taxCO2_Funneled, so adjust the second rescaling
                                    );
-            loop(regi,                   !! this loop necessary to allow the <-comparison in the next if statement
+            loop(regi,                   !! this loop is necessary to allow the <-comparison in the next if statement
               if( p_taxCO2eq_until2150(t2,regi) < pm_taxCO2eq(t2,regi) ,   !! check if new price would be higher than the price if the peakBudgYr would be one timestep later 
-                pm_taxCO2eq(t2,regi) = p_taxCO2eq_until2150(t2,regi);
+                display "price increase reached price from path with cm_peakBudgYr one timestep later - downscale to 99%"; 
+				pm_taxCO2eq(t2,regi) = 0.99 * p_taxCO2eq_until2150(t2,regi); !! reduce the new CO2 price to 99% of the price that it would be if the peaking year was one timestep later. The next iteration will show if this is enough, otherwise cm_peakBudgYr will be shifted right 
                 o_reached_until2150pricepath(iteration) = 1;             !! upward CO2 price correction reached the continued price path - check in next iteration if this is high enough.  
               );
             );
@@ -490,7 +509,12 @@ if(cm_iterative_target_adj eq 9,
       cm_peakBudgYr = o_peakBudgYr_Itr(iteration+1);  !! this has to happen outside the loop, otherwise the loop condition might be true twice
     ); !! if o_delay_increase_peakBudgYear(iteration) = 1,   !! if there was a flip-floping in the previous iterations, try to solve this
 	
-    display o_delay_increase_peakBudgYear, o_reached_until2150pricepath, pm_taxCO2eq, o_peakBudgYr_Itr;
+	
+	loop(regi, !! not a nice solution to having only the price of one regi display (for better visibility), but this way it overwrites again and again until the value from the last regi remain
+	    o_taxCO2eq_afterPeakShiftLoop_Itr_1regi(t,iteration+1) = pm_taxCO2eq(t,regi); 
+	);
+	
+    display o_delay_increase_peakBudgYear, o_reached_until2150pricepath, pm_taxCO2eq, o_peakBudgYr_Itr, o_taxCO2eq_afterPeakShiftLoop_Itr_1regi, o_pkBudgYr_flipflop;
   ); !! if cm_emiscen eq 9,
 );   !! if cm_iterative_target_adj eq 8,
 

modules/45_carbonprice/diffCurvPhaseIn2Lin/datainput.gms

@@ -47,9 +47,11 @@ loop(ttot$((ttot.val ge cm_startyear) AND (ttot.val le cm_CO2priceRegConvEndYr)
    min(1,
        max(0, 
 	        p45_phasein_2025ratio(regi) + (1 - p45_phasein_2025ratio(regi)) 
-			                               * ( ( (ttot.val - 2025) + (cm_CO2priceRegConvEndYr - 2025) * 0.1 ) 
-                                               / ( (cm_CO2priceRegConvEndYr - 2025) * 1.1 ) 
-											 ) ** 2
+			                               * Power( 
+										       ( (ttot.val - 2025) + (cm_CO2priceRegConvEndYr - 2025) * 0.1 ) 
+                                               / ( (cm_CO2priceRegConvEndYr - 2025) * 1.1 )
+											   , 2
+											 ) !! use Power instead of ** to allow ttot be smaller than 2025, and thus the base to be negative
        )											 
    );
 );