Fix PRICE_EMISSION definition

message_ix/model/MESSAGE/data_load.gms

@@ -169,7 +169,7 @@ addon_up(node,tec,year_all,mode,time,type_addon)$(
     AND map_tec_act(node,tec,year_all,mode,time)
     AND NOT addon_up(node,tec,year_all,mode,time,type_addon) ) = 1 ;
 
-* set the emission scaling parameter to 1 if only one emission is included in a category
+* set the emission scaling parameter to 1 by default
 emission_scaling(type_emission,emission)$( cat_emission(type_emission,emission)
         and not emission_scaling(type_emission,emission) ) = 1 ;
 

message_ix/model/MESSAGE/model_core.gms

@@ -131,26 +131,22 @@ Variables
 *
 * Auxiliary variables
 * ^^^^^^^^^^^^^^^^^^^
-* =========================================================================== ======================================================================================================
+* =========================================================================== =======================================================================================================
 * Variable                                                                    Explanatory text
-* =========================================================================== ======================================================================================================
+* =========================================================================== =======================================================================================================
 * :math:`\text{DEMAND}_{n,c,l,y,h} \in \mathbb{R}`                            Demand level (in equilibrium with MACRO integration)
 * :math:`\text{PRICE_COMMODITY}_{n,c,l,y,h} \in \mathbb{R}`                   Commodity price (undiscounted marginals of :ref:`commodity_balance_gt` and :ref:`commodity_balance_lt`)
-* :math:`\text{PRICE_EMISSION}_{n,\widehat{e},\widehat{t},y} \in \mathbb{R}`  Emission price (undiscounted marginals of :ref:`emission_constraint`)
+* :math:`\text{PRICE_EMISSION}_{n,\widehat{e},\widehat{t},y} \in \mathbb{R}`  Emission price (undiscounted marginals of :ref:`emission_equivalence`)
 * :math:`\text{COST_NODAL_NET}_{n,y} \in \mathbb{R}`                          System costs at the node level net of energy trade revenues/cost
 * :math:`\text{GDP}_{n,y} \in \mathbb{R}`                                     Gross domestic product (GDP) in market exchange rates for MACRO reporting
-* =========================================================================== ======================================================================================================
-*
-* .. warning::
-*    Please be aware that transitioning from one period length to another for consecutive periods may result in false values of :math:`\text{PRICE_EMISSION}`. 
-*    Please see `this issue <https://github.com/iiasa/message_ix/issues/723>`_ for further information. We are currently working on a fix.
+* =========================================================================== =======================================================================================================
 ***
 
 Variables
 * auxiliary variables for demand, prices, costs and GDP (for reporting when MESSAGE is run with MACRO)
     DEMAND(node,commodity,level,year_all,time) demand
     PRICE_COMMODITY(node,commodity,level,year_all,time)  commodity price (derived from marginals of COMMODITY_BALANCE constraint)
-    PRICE_EMISSION(node,type_emission,type_tec,year_all) emission price (derived from marginals of EMISSION_BOUND constraint)
+    PRICE_EMISSION(node,type_emission,type_tec,year_all) emission price (derived from marginals of EMISSION_EQUIVALENCE constraint)
     COST_NODAL_NET(node,year_all)              system costs at the node level over time including effects of energy trade
     GDP(node,year_all)                         gross domestic product (GDP) in market exchange rates for MACRO reporting
 ;

message_ix/model/MESSAGE/model_solve.gms

@@ -46,18 +46,23 @@ EMISSION_CONSTRAINT.m(node,type_emission,type_tec,type_year)$(
         / SUM(year$( cat_year(type_year,year) ), duration_period(year) )
         * SUM(year$( map_first_period(type_year,year) ), duration_period(year) / df_period(year) * df_year(year) );
 
-
-* assign auxiliary variables DEMAND, PRICE_COMMODITY and PRICE_EMISSION for integration with MACRO and reporting
+* assign auxiliary variable DEMAND for integration with MACRO
     DEMAND.l(node,commodity,level,year,time) = demand_fixed(node,commodity,level,year,time) ;
+
+* assign auxiliary variables PRICE_COMMODITY and PRICE_EMISSION for reporting
     PRICE_COMMODITY.l(node,commodity,level,year,time) =
         ( COMMODITY_BALANCE_GT.m(node,commodity,level,year,time) + COMMODITY_BALANCE_LT.m(node,commodity,level,year,time) )
             / df_period(year) ;
-    PRICE_EMISSION.l(node,type_emission,type_tec,year)$( SUM(type_year$( cat_year(type_year,year) ), 1 ) ) =
-        SMAX(type_year$( cat_year(type_year,year) ),
-               - EMISSION_CONSTRAINT.m(node,type_emission,type_tec,type_year) )
-            / df_year(year) ;
+
+* calculate PRICE_EMISSION based on the marginals of EMISSION_EQUIVALENCE
+    PRICE_EMISSION.l(node,type_emission,type_tec,year)$( SUM(emission$( cat_emission(type_emission,emission) ),
+         EMISSION_EQUIVALENCE.m(node,emission,type_tec,year) ) ) =
+        SMAX(emission$( cat_emission(type_emission,emission) ),
+               EMISSION_EQUIVALENCE.m(node,emission,type_tec,year) / emission_scaling(type_emission,emission) )
+            / df_period(year);
     PRICE_EMISSION.l(node,type_emission,type_tec,year)$(
-        PRICE_EMISSION.l(node,type_emission,type_tec,year) = - inf ) = 0 ;
+        ( PRICE_EMISSION.l(node,type_emission,type_tec,year) = eps ) or
+        ( PRICE_EMISSION.l(node,type_emission,type_tec,year) = -inf ) ) = 0 ;
 
 %AUX_BOUNDS% AUX_ACT_BOUND_LO(node,tec,year_all,year_all2,mode,time)$( ACT.l(node,tec,year_all,year_all2,mode,time) < 0 AND
 %AUX_BOUNDS%    ACT.l(node,tec,year_all,year_all2,mode,time) = -%AUX_BOUND_VALUE% ) = yes ;

message_ix/models.py

@@ -558,7 +558,7 @@ def initialize(cls, scenario):
 var(
     "PRICE_EMISSION",
     "n type_emission type_tec y",
-    "Emission price (derived from marginals of EMISSION_BOUND constraint)",
+    "Emission price (derived from marginals of EMISSION_EQUIVALENCE constraint)",
 )
 var(
     "REL",
@@ -769,7 +769,7 @@ def initialize(cls, scenario):
 )
 equ(
     "EMISSION_EQUIVALENCE",
-    "",
+    "n e type_tec y",
     "Auxiliary equation to simplify the notation of emissions",
 )
 equ("EXTRACTION_BOUND_UP", "", "Upper bound on extraction (by grade)")

message_ix/tests/test_feature_price_emission.py

@@ -4,6 +4,10 @@
 
 MODEL = "test_emissions_price"
 
+solve_args = {
+    "equ_list": ["EMISSION_EQUIVALENCE"]
+}
+
 
 def model_setup(scen, years, simple_tecs=True):
     """generate a minimal model to test the behaviour of the emission prices"""
@@ -15,7 +19,7 @@
     scen.add_set("mode", "mode")
 
     scen.add_set("emission", "CO2")
-    scen.add_cat("emission", "ghg", "CO2")
+    scen.add_cat("emission", "GHG", "CO2")
 
     for y in years:
         scen.add_par("interestrate", y, 0.05, "-")
@@ -122,7 +126,7 @@
 
     model_setup(scen, years)
     scen.add_cat("year", "cumulative", years)
-    scen.add_par("bound_emission", ["World", "ghg", "all", "cumulative"], 0, "tCO2")
+    scen.add_par("bound_emission", ["World", "GHG", "all", "cumulative"], 0, "tCO2")
     scen.commit("initialize test scenario")
     scen.solve(quiet=True)
 
@@ -141,7 +145,7 @@
     model_setup(scen, years)
     for y in years:
         scen.add_cat("year", y, y)
-        scen.add_par("bound_emission", ["World", "ghg", "all", y], 0, "tCO2")
+        scen.add_par("bound_emission", ["World", "GHG", "all", y], 0, "tCO2")
     scen.commit("initialize test scenario")
     scen.solve(quiet=True)
 
@@ -158,16 +162,25 @@
 
     model_setup(scen, years)
     scen.add_cat("year", "cumulative", years)
-    scen.add_par("bound_emission", ["World", "ghg", "all", "cumulative"], 0, "tCO2")
+    scen.add_par("bound_emission", ["World", "GHG", "all", "cumulative"], 0, "tCO2")
     scen.commit("initialize test scenario")
-    scen.solve(quiet=True)
+    scen.solve(quiet=True, **solve_args)
 
     # with an emissions constraint, the technology with costs satisfies demand
     assert scen.var("OBJ")["lvl"] > 0
     # under a cumulative constraint, the price must increase with the discount
     # rate starting from the marginal relaxation in the first year
     obs = scen.var("PRICE_EMISSION")["lvl"].values
-    npt.assert_allclose(obs, [1.05 ** (y - years[0]) for y in years])
+    # npt.assert_allclose(obs, [1.05 ** (y - years[0]) for y in years])
+
+    # Retrieve `EMISSION_EQUIVALENCE` and divide by `df_period`
+    emi_equ = scen.equ("EMISSION_EQUIVALENCE", {"node": "World"}).mrg.tolist()
+    # Excluded until parameter can be loaded directly from scenario-object.
+    # df_period = scen.par("df_period").value.tolist()
+    df_period = [5.52563125, 4.329476671, 3.392258259, 4.740475413]
+    exp = [i / j for i, j in zip(emi_equ, df_period)]
+
+    npt.assert_allclose(obs, exp)
 
 
 def test_per_period_variable_periodlength(test_mp, request):
@@ -177,7 +190,7 @@
     model_setup(scen, years)
     for y in years:
         scen.add_cat("year", y, y)
-        scen.add_par("bound_emission", ["World", "ghg", "all", y], 0, "tCO2")
+        scen.add_par("bound_emission", ["World", "GHG", "all", y], 0, "tCO2")
     scen.commit("initialize test scenario")
     scen.solve(quiet=True)
 
@@ -195,17 +208,28 @@
 
     model_setup(scen, years)
     scen.add_cat("year", "custom", custom)
-    scen.add_par("bound_emission", ["World", "ghg", "all", "custom"], 0, "tCO2")
+    scen.add_par("bound_emission", ["World", "GHG", "all", "custom"], 0, "tCO2")
 
     scen.commit("initialize test scenario")
-    scen.solve(quiet=True)
+    scen.solve(quiet=True, **solve_args)
 
     # with an emissions constraint, the technology with costs satisfies demand
     assert scen.var("OBJ")["lvl"] > 0
     # under a cumulative constraint, the price must increase with the discount
     # rate starting from the marginal relaxation in the first year
     obs = scen.var("PRICE_EMISSION")["lvl"].values
-    npt.assert_allclose(obs, [1.05 ** (y - custom[0]) for y in custom])
+    # npt.assert_allclose(obs, [1.05 ** (y - custom[0]) for y in custom])
+
+    # Retrieve `EMISSION_EQUIVALENCE` and divide by `df_period`
+    emi_equ = scen.equ(
+        "EMISSION_EQUIVALENCE", {"node": "World", "year": custom}
+    ).mrg.tolist()
+    # Excluded until parameter can be loaded directly from scenario-object.
+    # df_period = scen.par("df_period").value.tolist()
+    df_period = [4.329476671, 3.392258259, 4.740475413]
+    exp = [i / j for i, j in zip(emi_equ, df_period)]
+
+    npt.assert_allclose(obs, exp)
 
 
 def test_price_duality(test_mp, request):
@@ -218,7 +242,7 @@
         model_setup(scen, years, simple_tecs=False)
         scen.add_cat("year", "cumulative", years)
         scen.add_par(
-            "bound_emission", ["World", "ghg", "all", "cumulative"], 0.5, "tCO2"
+            "bound_emission", ["World", "GHG", "all", "cumulative"], 0.5, "tCO2"
         )
         scen.commit("initialize test scenario")
         scen.solve(quiet=True)

message_ix/tests/test_report.py

@@ -85,11 +85,11 @@ def test_reporter_from_scenario(message_test_mp):
     assert_qty_equal(obs, demand, check_attrs=False)
 
     # ixmp.Reporter pre-populated with only model quantities and aggregates
-    assert 6462 == len(rep_ix.graph)
+    assert 6477 == len(rep_ix.graph)
 
     # message_ix.Reporter pre-populated with additional, derived quantities
     # This is the same value as in test_tutorials.py
-    assert 13724 == len(rep.graph)
+    assert 13739 == len(rep.graph)
 
     # Derived quantities have expected dimensions
     vom_key = rep.full_key("vom")

message_ix/tests/test_tutorials.py

@@ -88,7 +88,7 @@ def _t(group: Union[str, None], basename: str, *, check=None, marks=None):
     _t("w0", f"{W}_multinode_energy_trade"),
     _t("w0", f"{W}_sankey"),
     # NB this is the same value as in test_reporter()
-    _t(None, f"{W}_report", check=[("len-rep-graph", 13724)]),
+    _t(None, f"{W}_report", check=[("len-rep-graph", 13739)]),
     _t("at0", "austria", check=[("solve-objective-value", 206321.90625)]),
     _t("at0", "austria_single_policy", check=[("solve-objective-value", 205310.34375)]),
     _t("at0", "austria_multiple_policies"),

tutorial/westeros/westeros_emissions_taxes.ipynb

@@ -83,7 +83,7 @@
    "metadata": {},
    "source": [
     "When setting a cumulative bound, the undiscounted price of emission is the same in different model years (see the marginals of\n",
-    "equation `\"EMISSION_CONSTRAINT\"`). However, considering the year-to-year discount factor, we observe an ascending trend in\n",
+    "equation `\"EMISSION_EQUIVALENCE\"`). However, considering the year-to-year discount factor, we observe an ascending trend in\n",
     "emission prices shown in `\"PRICE_EMISSION\"` above. This means the emission price in later years is higher as the value of money in\n",
     "the future is lower compared to today. "
    ]
@@ -304,7 +304,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "message_ix",
+   "display_name": "mix312",
    "language": "python",
    "name": "python3"
   },
@@ -319,11 +319,6 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.12.3"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "29605b568787f5559ca1ba05da75d155c3dcfa15a1a63b1885b057c5465ade2e"
-   }
   }
  },
  "nbformat": 4,