Add option for MESSAGE-MAgPIE connection

message_ix/model/MESSAGE/data_load.gms

@@ -103,6 +103,33 @@ Set rating_unrated(rating) ;
 rating_unrated(rating) = yes ;
 rating_unrated('unrated') = no ;
 
+*----------------------------------------------------------------------------------------------------------------------*
+* auxiliary mappings for split between cumulative and annualland-use emission calculations                             *
+*----------------------------------------------------------------------------------------------------------------------*
+
+* remove all emission types from subset for cumulative land-use emission accounting
+emission_cumulative(emission) = no;
+* add emission types where land-use emissions should use the cumulative accounting approach
+* TO DO: to avoid having to declarde these emissions in gams code, this section should be translated into a Python process
+emission_cumulative('TCE_CO2') = yes;
+emission_cumulative('LU_CO2') = yes;
+emission_cumulative('TCE') = yes;
+* add all emission types to subset for annual land-use emission accounting,
+* then remove the emission types that have been declared in the cumulative accounting set
+emission_annual(emission) = yes;
+emission_annual(emission_cumulative) = no;
+
+*----------------------------------------------------------------------------------------------------------------------*
+* auxiliary mappings for split between cumulative and annualland-use emission calculations                             *
+*----------------------------------------------------------------------------------------------------------------------*
+
+emission_cumulative(emission) = no;
+emission_cumulative('TCE_CO2') = yes;
+emission_cumulative('LU_CO2') = yes;
+emission_cumulative('TCE') = yes;
+emission_annual(emission) = yes;
+emission_annual(emission_cumulative) = no;
+
 *----------------------------------------------------------------------------------------------------------------------*
 * assignment and computation of MESSAGE-specific auxiliary parameters                                                  *
 *----------------------------------------------------------------------------------------------------------------------*

message_ix/model/MESSAGE/model_core.gms

@@ -105,10 +105,15 @@ Positive Variables
     ACT_LO(node,tec,year_all,time)   relaxation variable for dynamic constraints on activity (downwards)
 * land-use model emulator
     LAND(node,land_scenario,year_all) relative share of land-use scenario
+    LAND_COST_NEW(node, year_all)     Land cost including debt from scenario switching
+    LAND_COST_DEBT(node, year_all,year_all2) Land cost debt from scenario switching 
+    EMISS_LU_AUX(node,emission,type_tec,year_all,year_all2)    positive emissions overshoot of historic emissions compared to chosen land scenario mix
+
 * content of storage
     STORAGE(node,tec,mode,level,commodity,year_all,time)       state of charge (SoC) of storage at each sub-annual time slice (positive)
 ;
 
+
 Variables
 * intertemporal stock variables (input or output quantity into the stock)
     STOCK_CHG(node,commodity,level,year_all,time) annual input into and output from stocks of commodities
@@ -120,12 +125,18 @@ Variables
     COST_NODAL(node, year_all)                   system costs at the node level over time
 * auxiliary variable for aggregate emissions by technology type and land-use model emulator
     EMISS(node,emission,type_tec,year_all)       aggregate emissions by technology type and land-use model emulator
+* auxiliary variable for aggregate emissions from land-use model emulator
+    EMISS_LU(node,emission,type_tec,year_all)    aggregate emissions from land-use model emulator
 * auxiliary variable for left-hand side of relations (linear constraints)
     REL(relation,node,year_all)                  auxiliary variable for left-hand side of user-defined relations
 * change in the content of storage device
     STORAGE_CHARGE(node,tec,mode,level,commodity,year_all,time)    charging of storage in each time slice (negative for discharge)
 ;
 
+
+
+
+
 ***
 * .. _section_auxiliary_variable_def:
 *
@@ -280,12 +291,18 @@ Equations
     ACTIVITY_CONSTRAINT_LO          dynamic constraint on the market penetration of a technology activity (lower bound)
     ACTIVITY_SOFT_CONSTRAINT_LO     bound on relaxation of the dynamic constraint on market penetration (lower bound)
     EMISSION_EQUIVALENCE            auxiliary equation to simplify the notation of emissions
+    EMISSION_EQUIVALENCE_AUX_ANNUAL auxiliary equation calculating land-use emissions from annual scenario input
+    EMISSION_EQUIVALENCE_AUX_CUMU   auxiliary equation calculating land-use emissions from cumulative scenario input
+    EMISSION_EQUIVALENCE_AUX_CUMU_AUX auxiliary equation calculating the land-use emissions overshoot if positive compared to historic scenario mix
     EMISSION_CONSTRAINT             nodal-regional-global constraints on emissions (by category)
+    LAND_COST_CUMU                  land cost including debt from scenario switching                                 
+    LAND_COST_CUMU_DEBT             land cost debt from scenario switching
     LAND_CONSTRAINT                 constraint on total land use (linear combination of land scenarios adds up to 1)
     DYNAMIC_LAND_SCEN_CONSTRAINT_UP dynamic constraint on land scenario change (upper bound)
     DYNAMIC_LAND_SCEN_CONSTRAINT_LO dynamic constraint on land scenario change (lower bound)
     DYNAMIC_LAND_TYPE_CONSTRAINT_UP dynamic constraint on land-use change (upper bound)
     DYNAMIC_LAND_TYPE_CONSTRAINT_LO dynamic constraint on land-use change (lower bound)
+*    TAU_CONSTRAINT                  constraint on land-use intensity growth (regional tau is not allowed to shrink)
     RELATION_EQUIVALENCE            auxiliary equation to simplify the implementation of relations
     RELATION_CONSTRAINT_UP          upper bound of relations (linear constraints)
     RELATION_CONSTRAINT_LO          lower bound of relations (linear constraints)
@@ -411,8 +428,9 @@ COST_ACCOUNTING_NODAL(node, year)..
         * tax_emission(node,type_emission,type_tec,type_year)
         * EMISS(node,emission,type_tec,year) )
 * cost terms from land-use model emulator (only includes valid node-land_scenario-year combinations)
-    + SUM(land_scenario$( land_cost(node,land_scenario,year) ),
-        land_cost(node,land_scenario,year) * LAND(node,land_scenario,year) )
+*    + SUM(land_scenario$( land_cost(node,land_scenario,year) ),
+*        land_cost(node,land_scenario,year) * LAND(node,land_scenario,year) )
+    + LAND_COST_NEW(node, year)
 * cost terms associated with linear relations
     + SUM(relation$( relation_cost(relation,node,year) ),
         relation_cost(relation,node,year) * REL(relation,node,year) )
@@ -1836,8 +1854,7 @@ ACTIVITY_SOFT_CONSTRAINT_LO(node,tec,year,time)$( soft_activity_lo(node,tec,year
 *          \sum_{n^L \in N(n)} \Bigg(
 *              \sum_{t \in T(\widehat{t}),y^V \leq y,m,h }
 *                  \text{emission_factor}_{n^L,t,y^V,y,m,e} \cdot \text{ACT}_{n^L,t,y^V,y,m,h} \\
-*              + \sum_{s} \ \text{land_emission}_{n^L,s,y,e} \cdot \text{LAND}_{n^L,s,y}
-*                   \text{ if } \widehat{t} \in \widehat{T}^{LAND} \Bigg)
+*              + \text{EMISS_LU}_{n^L,e,\widehat{t},y} \text{ if } \widehat{t} \in \widehat{T}^{LAND} \Bigg)
 *
 ***
 EMISSION_EQUIVALENCE(node,emission,type_tec,year)..
@@ -1848,11 +1865,87 @@ EMISSION_EQUIVALENCE(node,emission,type_tec,year)..
         SUM((tec,vintage,mode,time)$( cat_tec(type_tec,tec)
             AND map_tec_act(location,tec,year,mode,time) AND map_tec_lifetime(location,tec,vintage,year) ),
         emission_factor(location,tec,vintage,year,mode,emission) * ACT(location,tec,vintage,year,mode,time) )
-* emissions from land use if 'type_tec' is included in the dynamic set 'type_tec_land'
-        + SUM(land_scenario$( type_tec_land(type_tec) ),
-            land_emission(location,land_scenario,year,emission) * LAND(location,land_scenario,year) )
+* emissions from land-use activity based on land scenario mix and its path dependencies 
+* for selected variables (emission_cumulative) if 'type_tec' is included in the dynamic set 'type_tec_land'
+        + EMISS_LU(location,emission,type_tec,year) $ ( type_tec_land(type_tec) )
       ) ;
 
+
+* .. _equation_emission_equivalence_aux_annual:
+*
+* Equation EMISSION_EQUIVALENCE_AUX_ANNUAL
+* """""""""""""""""""""""""""""
+* This auxillary equation accounts for emissions in set emission_annual without considering a path-dependent emission debt
+*
+*
+*   .. math::
+*      \text{EMISS_LU}_{n^L,e^a,\widehat{t},y} =
+*               \sum_{s} \text{land_emission}_{n^L,s,y,e^a} \cdot \text{LAND}_{n^L,s,y}
+*
+***
+EMISSION_EQUIVALENCE_AUX_ANNUAL(location,emission,type_tec,year) $ emission_annual(emission)..
+    EMISS_LU(location,emission,type_tec,year)
+    =E=
+    SUM(land_scenario ,
+            land_emission(location,land_scenario,year,emission) * LAND(location,land_scenario,year)
+    ) ;
+
+
+* .. _equation_emission_equivalence_aux_cumu:
+*
+* Equation EMISSION_EQUIVALENCE_AUX_CUMU
+* """""""""""""""""""""""""""""
+* This auxillary equation accounts for emissions in set emission_cumulative considering a path-dependent emission debt
+* This debt is calculated by comparing the historic emissions of the current time step's land-use mix in variable LAND 
+* to the emissions of the actually chosen scenario mix of previous timesteps, including any emission debts accounted for 
+* before. To do so, first, for every previous model time step, the emissions in those steps under the current land-use 
+* scenario mix and the delta to the chosen mix in that time step are calculated. To ensure any historic emission debt is
+* only accounted for the first time step(s) it occurs, the accounted debt towards any previous time step is subtracted 
+* towards a minimum of 0.
+* 
+*
+*   .. math::
+*      \text{EMISS_LU}_{n^L,e^c,\widehat{t},y} =
+*               \sum_{s} \text{land_emission}_{n^L,s,y,e^c} \cdot \text{LAND}_{n^L,s,y} \\
+*              + \sum_{\widehat{y}} \text{EMISS_LU_AUX}_{n^L,e^c,\widehat{t},y,\widehat{y}}
+*
+***
+EMISSION_EQUIVALENCE_AUX_CUMU(location,emission,type_tec,year) $ emission_cumulative(emission)..
+    EMISS_LU(location,emission,type_tec,year)
+    =E=
+    SUM(land_scenario ,
+            land_emission(location,land_scenario,year,emission) * LAND(location,land_scenario,year) )
+    + SUM(year2, EMISS_LU_AUX(location,emission,type_tec,year, year2) )
+    ;
+
+
+* find positive emissions overshoot for history of current land scenario mix compared to mix of earlier time steps
+EMISSION_EQUIVALENCE_AUX_CUMU_AUX(location,emission,type_tec,year,year2) $ (emission_cumulative(emission) AND model_horizon(year) AND year2.pos < year.pos)..
+    EMISS_LU_AUX(location,emission,type_tec,year,year2) $ ( year2.pos < year.pos ) 
+    =G=
+    SUM(land_scenario,
+            LAND(location, land_scenario, year) * land_emission(location, land_scenario, year2, emission)
+            - LAND(location, land_scenario, year2) * land_emission(location, land_scenario, year2, emission) )
+        - SUM(year3 $ ( year3.pos < year.pos AND year2.pos < year3.pos ), EMISS_LU_AUX(location, emission,type_tec, year3, year2) ) ;
+
+* calculate scenario cost of current land scenario mix under consideration of its path dependencies and 
+* associated cost differences compared to the land scenario mix of earlier time steps
+LAND_COST_CUMU(location, year)$( model_horizon(year) )..
+       LAND_COST_NEW(location, year) =E=
+       SUM(land_scenario$( land_cost(location,land_scenario,year) ),
+       land_cost(location,land_scenario,year) * LAND(location,land_scenario,year) )
+       + SUM(year2, 
+            LAND_COST_DEBT(location, year, year2)
+            ) ; 
+
+LAND_COST_CUMU_DEBT(location, year, year2) $ (model_horizon(year) AND year2.pos < year.pos)..
+        LAND_COST_DEBT(location, year, year2) $ ( year2.pos < year.pos ) =G=
+        SUM(land_scenario$( land_cost(location,land_scenario,year) ),
+            LAND(location, land_scenario, year) * land_cost(location,land_scenario,year2)
+            - LAND(location, land_scenario, year2) * land_cost(location,land_scenario,year2) ) * df_period(year2) / df_period(year)
+        - SUM(year3 $ ( year3.pos < year.pos AND year2.pos < year3.pos ), LAND_COST_DEBT(location, year3, year2) ) ;
+
+
 ***
 * Bound on emissions
 * ^^^^^^^^^^^^^^^^^^
@@ -1914,6 +2007,7 @@ EMISSION_CONSTRAINT(node,type_emission,type_tec,type_year)$is_bound_emission(nod
 LAND_CONSTRAINT(node,year)$( SUM(land_scenario$( map_land(node,land_scenario,year) ), 1 ) ) ..
     SUM(land_scenario$( map_land(node,land_scenario,year) ), LAND(node,land_scenario,year) ) =E= 1 ;
 
+
 ***
 * Dynamic constraints on land use
 * ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -2055,7 +2149,7 @@ DYNAMIC_LAND_TYPE_CONSTRAINT_LO(node,year,land_type)$( is_dynamic_land_lo(node,y
         land_use(node,land_scenario,year,land_type) * LAND(node,land_scenario,year) ) =G=
 * initial 'new' land used for that type (compounded over the duration of the period)
         - initial_land_lo(node,year,land_type) * (
-            ( ( POWER( 1 + growth_land_up(node,year,land_type) , duration_period(year) ) - 1 )
+            ( ( POWER( 1 + growth_land_lo(node,year,land_type) , duration_period(year) ) - 1 )
                 / growth_land_lo(node,year,land_type) )$( growth_land_lo(node,year,land_type) )
               + ( duration_period(year) )$( NOT growth_land_lo(node,year,land_type) )
             )
@@ -2073,6 +2167,39 @@ DYNAMIC_LAND_TYPE_CONSTRAINT_LO(node,year,land_type)$( is_dynamic_land_lo(node,y
 %SLACK_LAND_TYPE_LO% - SLACK_LAND_TYPE_LO(node,year,land_type)
 ;
 
+
+***
+*
+* To ensure correct replication of MAgPIE behavior, the following constraint ensures that 
+* land-use intensity cannot decrease over time
+*
+* .. _equation_tau_constraint:
+*
+* Equation TAU_CONSTRAINT
+* """"""""""""""""""""""""""""""""""""""""
+*
+*  .. math::
+*     \sum_{s \in S} \text{land_output Tau}_{n,s,y,l,h} &\cdot \text{LAND}_{n,s,y}
+*         \geq \sum_{s \in S} \big( \text{land_output Tau}_{n,s,y-1,l,h}
+*                            & \quad \quad \cdot \big( \text{LAND}_{n,s,y-1} + \text{historical_land}_{n,s,y-1} \big) \big) 
+*
+***
+
+TAU_CONSTRAINT(node, year, level, time) .. 
+   SUM(land_scenario$( map_land(node,land_scenario,year) ),
+        land_output(node, land_scenario, year, "Landuse intensity indicator Tau", level, time) 
+          * LAND(node, land_scenario, year)
+        ) =G=
+    SUM((year_all2)$( seq_period(year_all2,year) ),
+        SUM(land_scenario$( map_land(node,land_scenario,year) ),
+            land_output(node, land_scenario, year_all2, "Landuse intensity indicator Tau", level, time) 
+              * ( LAND(node, land_scenario, year_all2) $ ( model_horizon(year_all2) ) 
+                  + historical_land(node,land_scenario,year_all2) )
+            )
+      )
+;
+
+
 *----------------------------------------------------------------------------------------------------------------------*
 ***
 * .. _section_of_generic_relations:

message_ix/model/MESSAGE/sets_maps_def.gms

@@ -148,6 +148,8 @@ Sets
     tec             technologies
     mode            modes of operation
     emission        greenhouse gases - pollutants - etc.
+    emission_annual (emission) greenhouse gases - pollutants - etc. from annual inputs
+    emission_cumulative (emission) greenhouse gases - pollutants - etc. from cumulative inputs 
     land_scenario   scenarios of land use (for land-use model emulator)
     land_type       types of land use
     year_all        years (over entire model horizon)