diff --git a/message_ix/model/MESSAGE/data_load.gms b/message_ix/model/MESSAGE/data_load.gms index 1f094303b..d31306c54 100644 --- a/message_ix/model/MESSAGE/data_load.gms +++ b/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 ; diff --git a/message_ix/model/MESSAGE/model_core.gms b/message_ix/model/MESSAGE/model_core.gms index e549c7fb1..5c7b3afb3 100644 --- a/message_ix/model/MESSAGE/model_core.gms +++ b/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{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 `_ for further information. We are currently working on a fix. +* ======================================================================================= ======================================================================================================= +* 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_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 +* ======================================================================================= ======================================================================================================= *** 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 ; diff --git a/message_ix/model/MESSAGE/model_solve.gms b/message_ix/model/MESSAGE/model_solve.gms index 6fda3e23f..93cdacd47 100644 --- a/message_ix/model/MESSAGE/model_solve.gms +++ b/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 ; diff --git a/message_ix/models.py b/message_ix/models.py index 5c9253f71..b3fec43f1 100644 --- a/message_ix/models.py +++ b/message_ix/models.py @@ -4,7 +4,7 @@ from dataclasses import InitVar, dataclass, field from functools import partial from pathlib import Path -from typing import Mapping, MutableMapping, Optional, Tuple +from typing import MutableMapping, Optional, Tuple from warnings import warn import ixmp.model.gams @@ -175,7 +175,7 @@ class GAMSModel(ixmp.model.gams.GAMSModel): #: Mapping from model item (equation, parameter, set, or variable) names to #: :class:`.Item` describing the item. - items: Mapping[str, Item] + items: MutableMapping[str, Item] def __init__(self, name=None, **model_options): # Update the default options with any user-provided options @@ -557,7 +557,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", @@ -768,7 +768,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)") diff --git a/message_ix/tests/feature_price_emission.ipynb b/message_ix/tests/feature_price_emission.ipynb new file mode 100644 index 000000000..b0fcfea8b --- /dev/null +++ b/message_ix/tests/feature_price_emission.ipynb @@ -0,0 +1,394 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from ixmp import Platform\n", + "\n", + "from message_ix import Scenario\n", + "\n", + "MODEL = \"test_emissions_price\"\n", + "\n", + "mp = Platform(\"test_feature_price_emission\")\n", + "mp.add_unit(\"MtCO2\")\n", + "mp.add_unit(\"tCO2/kWa\")\n", + "mp.add_unit(\"USD/kW\")\n", + "\n", + "scen = Scenario(\n", + " mp,\n", + " MODEL,\n", + " scenario=\"many_tecs\",\n", + " version=\"new\",\n", + ")\n", + "# scen = Scenario(\n", + "# mp,\n", + "# MODEL,\n", + "# scenario=\"two_tecs\",\n", + "# version=\"new\",\n", + "# )" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from message_ix.tests.test_feature_price_emission import model_setup\n", + "\n", + "# 450 is too much; this bound will not affect the first model year\n", + "# (EMISS is not reduced from without cumulative bound to with it)\n", + "# In a model year without binding bound, no price_emission is produced ->\n", + "# there is one value missing for price_emission in period-specific bound\n", + "cumulative_bound = 40\n", + "# cumulative_bound = 500\n", + "# cumulative_bound = 550\n", + "# years = [2020, 2030, 2040, 2050]\n", + "years = [2020, 2025, 2030, 2040, 2045, 2050]\n", + "# years = [2020, 2030, 2040, 2045, 2050]\n", + "\n", + "filters = {\"node\": \"World\"}\n", + "\n", + "model_setup(scen=scen, years=years, simple_tecs=False)\n", + "# model_setup(scen=scen, years=years)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from message_ix.tests.test_feature_price_emission import solve_args\n", + "\n", + "scen.commit(\"initialize test scenario\")\n", + "scen.solve(quiet=True, **solve_args)\n", + "scen.var(\"EMISS\", filters)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound = scen.clone(\n", + " MODEL,\n", + " \"cumulative_emission_bound\",\n", + " \"introducing a cumulative emissions bound\",\n", + " keep_solution=False,\n", + ")\n", + "# scenario_cumulative_bound = scen.clone(\n", + "# MODEL,\n", + "# \"cumulative_emission_bound_two_tecs\",\n", + "# \"introducing a cumulative emissions bound\",\n", + "# keep_solution=False,\n", + "# )\n", + "scenario_cumulative_bound.check_out()\n", + "\n", + "scenario_cumulative_bound.add_cat(\"year\", \"cumulative\", years)\n", + "scenario_cumulative_bound.add_par(\n", + " \"bound_emission\",\n", + " [\"World\", \"GHG\", \"all\", \"cumulative\"],\n", + " cumulative_bound,\n", + " \"MtCO2\",\n", + ")\n", + "scenario_cumulative_bound.commit(\"initialize test scenario\")\n", + "scenario_cumulative_bound.solve(quiet=True, **solve_args)\n", + "scenario_cumulative_bound.var(\"EMISS\", filters)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "emiss = scenario_cumulative_bound.var(\"EMISS\", filters).set_index(\"year\").lvl\n", + "price_emission = (\n", + " scenario_cumulative_bound.var(\"PRICE_EMISSION\", filters).set_index(\"year\").lvl\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# --------------------------------------------------------\n", + "# Run scenario with annual-emission bound based on `EMISS`\n", + "# from cumulative constraint scenario.\n", + "# --------------------------------------------------------\n", + "\n", + "scenario_period_bound = scen.clone(\n", + " MODEL,\n", + " \"period_bound_many_tecs\",\n", + " \"introducing a period-specific emission_bound\",\n", + " keep_solution=False,\n", + ")\n", + "# scenario_period_bound = scen.clone(\n", + "# MODEL,\n", + "# \"period_bound_two_tecs\",\n", + "# \"introducing a period-specific emission_bound\",\n", + "# keep_solution=False,\n", + "# )\n", + "scenario_period_bound.check_out()\n", + "for year in years:\n", + " scenario_period_bound.add_cat(\"year\", year, year)\n", + "\n", + "# use emissions from cumulative-constraint scenario as period-emission bounds\n", + "emiss_period_bound = (\n", + " scenario_cumulative_bound.var(\"EMISS\", {\"node\": \"World\"})\n", + " .rename(columns={\"year\": \"type_year\", \"lvl\": \"value\"})\n", + " .drop(\"emission\", axis=1)\n", + ")\n", + "emiss_period_bound[\"type_emission\"] = \"GHG\"\n", + "emiss_period_bound[\"unit\"] = \"MtCO2\"\n", + "scenario_period_bound.add_par(\"bound_emission\", emiss_period_bound)\n", + "scenario_period_bound.commit(\"initialize test scenario for periodic emission bound\")\n", + "scenario_period_bound.solve(quiet=True, **solve_args)\n", + "scenario_period_bound.var(\"EMISS\", filters)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy.testing as npt\n", + "\n", + "# check -emissions are close between cumulative and yearly-bound scenarios\n", + "emiss_period_bound = scenario_period_bound.var(\"EMISS\", filters).set_index(\"year\").lvl\n", + "npt.assert_allclose(emiss, emiss_period_bound)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound.equ(\"EMISSION_EQUIVALENCE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(price_emission)\n", + "scenario_period_bound.var(\"PRICE_EMISSION\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# check \"PRICE_EMISSION\" is close between cumulative- and yearly-bound scenarios\n", + "price_emission_period_bound = (\n", + " scenario_period_bound.var(\"PRICE_EMISSION\", filters).set_index(\"year\").lvl\n", + ")\n", + "npt.assert_allclose(price_emission, price_emission_period_bound)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_period_bound.equ(\"EMISSION_EQUIVALENCE\")\n", + "# NOTE: mrg == df_period for two_tec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound.par(\"bound_emission\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_period_bound.par(\"bound_emission\")" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "scen_tax = Scenario(\n", + " mp,\n", + " MODEL,\n", + " scenario=\"tax_many_tecs\",\n", + " version=\"new\",\n", + ")\n", + "# scen_tax = Scenario(\n", + "# mp,\n", + "# MODEL,\n", + "# scenario=\"tax_two_tecs\",\n", + "# version=\"new\",\n", + "# )\n", + "model_setup(scen_tax, years, simple_tecs=False)\n", + "# model_setup(scen_tax, years)\n", + "\n", + "for year in years:\n", + " scen_tax.add_cat(\"year\", year, year)\n", + "# use emission prices from cumulative-constraint scenario as taxes\n", + "taxes = scenario_cumulative_bound.var(\"PRICE_EMISSION\").rename(\n", + " columns={\"year\": \"type_year\", \"lvl\": \"value\"}\n", + ")\n", + "taxes[\"unit\"] = \"USD/tCO2\"\n", + "taxes[\"node\"] = \"node\"\n", + "scen_tax.add_par(\"tax_emission\", taxes)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scen_tax.commit(\"initialize test scenario for taxes\")\n", + "scen_tax.solve(quiet=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(scenario_cumulative_bound.var(\"PRICE_EMISSION\"))\n", + "print(scen_tax.var(\"PRICE_EMISSION\"))\n", + "print(scen_tax.par(\"tax_emission\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(scenario_cumulative_bound.var(\"EMISS\"))\n", + "# print(scenario_period_bound.var(\"EMISS\"))\n", + "print(scen_tax.var(\"EMISS\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cumulative_activity = scenario_cumulative_bound.var(\"ACT\")\n", + "print(cumulative_activity.loc[cumulative_activity.lvl > 0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "period_activity = scenario_period_bound.var(\"ACT\")\n", + "print(period_activity.loc[period_activity.lvl > 0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tax_activity = scen_tax.var(\"ACT\")\n", + "print(tax_activity.loc[tax_activity.lvl > 0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scen_tax.par(\"bound_emission\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound.par(\"demand\")" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "price_emission_tax = scen_tax.var(\"PRICE_EMISSION\").set_index(\"year\").lvl\n", + "npt.assert_allclose(price_emission, price_emission_tax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# check emissions are close between cumulative and tax scenarios\n", + "emiss_tax = scen_tax.var(\"EMISS\", filters).set_index(\"year\").lvl\n", + "npt.assert_allclose(emiss, emiss_tax, rtol=0.05)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "mp.close_db()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mix312", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/message_ix/tests/test_feature_price_emission.py b/message_ix/tests/test_feature_price_emission.py index b790965dc..be9d23b11 100644 --- a/message_ix/tests/test_feature_price_emission.py +++ b/message_ix/tests/test_feature_price_emission.py @@ -1,107 +1,248 @@ +from typing import List + import numpy.testing as npt +import pytest +# from message_ix.testing import make_westeros from message_ix import Scenario, make_df MODEL = "test_emissions_price" +solve_args = { + "equ_list": ["EMISSION_EQUIVALENCE"], + # "par_list": ["df_period", "df_year", "levelized_cost"], + # At the moment, it is not possible to retrieve auxilliary parameters that + # are created in GAMS back to ixmp. The default "par_list" is coded in Java + # backened, and specifying a list here does not add/modify the default list. +} + +interest_rate = 0.05 + -def model_setup(scen, years, simple_tecs=True): +def model_setup(scen: Scenario, years: List[int], simple_tecs=True) -> None: """generate a minimal model to test the behaviour of the emission prices""" scen.add_spatial_sets({"country": "node"}) scen.add_set("commodity", "comm") scen.add_set("level", "level") - scen.add_set("year", years) - + scen.add_horizon(years) 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, "-") - scen.add_par("demand", ["node", "comm", "level", y, "year"], 1, "GWa") + scen.add_par("interestrate", y, interest_rate, "-") + scen.add_par( + "demand", + make_df( + "demand", + node="node", + commodity="comm", + level="level", + time="year", + year=y, + value=60 + 2 * (y - years[0]), + unit="GWa", + ), + ) - if simple_tecs: - add_two_tecs(scen, years) - else: - add_many_tecs(scen, years) + add_two_tecs(scen, years) if simple_tecs else add_many_tecs(scen, years) -def add_two_tecs(scen, years): +def add_two_tecs(scen: Scenario, years: List[int]) -> None: """add two technologies to the scenario""" scen.add_set("technology", ["dirty_tec", "clean_tec"]) - common = dict(node_loc="node", year_vtg=years, year_act=years, value=1, mode="mode") + common_base = dict( + node_loc="node", year_vtg=years, year_act=years, mode="mode", value=1 + ) + common_output = dict( + node_dest="node", + commodity="comm", + level="level", + time="year", + time_dest="year", + unit="GWa", + ) # the dirty technology is free (no costs) but has emissions scen.add_par( "output", - make_df( - "output", - node_dest="node", - technology="dirty_tec", - commodity="comm", - level="level", - time="year", - time_dest="year", - unit="GWa", - **common, - ), + make_df("output", technology="dirty_tec", **common_base, **common_output), ) scen.add_par( "emission_factor", make_df( "emission_factor", - unit="tCO2", technology="dirty_tec", emission="CO2", - **common, + unit="tCO2", + **common_base, ), ) # the clean technology has variable costs but no emissions scen.add_par( "output", - make_df( - "output", - node_dest="node", - technology="clean_tec", - commodity="comm", - level="level", - time="year", - time_dest="year", - unit="GWa", - **common, - ), + make_df("output", technology="clean_tec", **common_base, **common_output), ) scen.add_par( "var_cost", make_df( "var_cost", + technology="clean_tec", time="year", unit="USD/GWa", - technology="clean_tec", - **common, + **common_base, ), ) -def add_many_tecs(scen, years, n=50): +def add_many_tecs(scen: Scenario, years: List[int], n=50) -> None: """add a range of dirty-to-clean technologies to the scenario""" - output_specs = ["node", "comm", "level", "year", "year"] - - for i in range(n + 1): - t = "tec{}".format(i) + # Add some hardcoded tecs for temporary testing + tecs: dict[str, dict] = { + "tec1": { + "emission_factor": 7.4, + "inv_cost": 500, + "fix_cost": 30, + "var_cost": 30, + "bound_activity_up": 100, + "lifetime": 20, + }, + "tec2": { + "emission_factor": 0, + "inv_cost": 1500, + "fix_cost": 10, + "var_cost": 0, + "bound_activity_up": 100, + "lifetime": 20, + }, + "tec3": { + "emission_factor": 2.5, + "inv_cost": 700, + "fix_cost": 30, + "var_cost": 30, + "bound_activity_up": 1000, + "lifetime": 20, + }, + "tec4": { + "emission_factor": 0, + "inv_cost": 2000, + "fix_cost": 5, + "var_cost": 0, + "bound_activity_up": 10, + "lifetime": 40, + }, + "tec5": { + "emission_factor": 10, + "inv_cost": 400, + "fix_cost": 30, + "var_cost": 30, + "bound_activity_up": 100, + "lifetime": 20, + }, + } + year_df = scen.vintage_and_active_years() + vintage_years, act_years = year_df["year_vtg"], year_df["year_act"] + mp = scen.platform + mp.add_unit("USD/kW") + mp.add_unit("tCO2/kWa") + + for t in tecs: scen.add_set("technology", t) - for y in years: - tec_specs = ["node", t, y, y, "mode"] - # variable costs grow quadratically over technologies - # to get rid of the curse of linearity - c = (10 * i / n) ** 2 * (1.045) ** (y - years[0]) - e = 1 - i / n - scen.add_par("output", tec_specs + output_specs, 1, "GWa") - scen.add_par("var_cost", tec_specs + ["year"], c, "USD/GWa") - scen.add_par("emission_factor", tec_specs + ["CO2"], e, "tCO2") + scen.add_par( + "output", + make_df( + "output", + node_loc="node", + node_dest="node", + technology=t, + year_vtg=vintage_years, + year_act=act_years, + mode="mode", + commodity="comm", + level="level", + time="year", + time_dest="year", + value=1, + unit="GWa", + ), + ) + scen.add_par( + "technical_lifetime", + make_df( + "technical_lifetime", + node_loc="node", + year_vtg=vintage_years, + unit="y", + technology=t, + value=tecs[t]["lifetime"], + ), + ) + scen.add_par( + "var_cost", + make_df( + "var_cost", + node_loc="node", + year_vtg=vintage_years, + year_act=act_years, + mode="mode", + time="year", + unit="USD/kWa", + technology=t, + value=tecs[t]["var_cost"], + ), + ) + scen.add_par( + "inv_cost", + make_df( + "inv_cost", + node_loc="node", + year_vtg=vintage_years, + unit="USD/kW", + technology=t, + value=tecs[t]["inv_cost"], + ), + ) + scen.add_par( + "fix_cost", + make_df( + "fix_cost", + node_loc="node", + year_vtg=vintage_years, + year_act=act_years, + unit="USD/kWa", + technology=t, + value=tecs[t]["fix_cost"], + ), + ) + scen.add_par( + "emission_factor", + make_df( + "emission_factor", + node_loc="node", + year_vtg=vintage_years, + year_act=act_years, + unit="tCO2/kWa", + technology=t, + mode="mode", + value=tecs[t]["emission_factor"], + emission="CO2", + ), + ) + scen.add_par( + "bound_activity_up", + make_df( + "bound_activity_up", + node_loc="node", + technology=t, + year_act=act_years, + mode="mode", + time="year", + value=tecs[t]["bound_activity_up"], + unit="GWa", + ), + ) def test_no_constraint(test_mp, request): @@ -122,7 +263,7 @@ def test_cumulative_equidistant(test_mp, request): 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) @@ -131,7 +272,7 @@ def test_cumulative_equidistant(test_mp, request): # 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 + interest_rate) ** (y - years[0]) for y in years]) def test_per_period_equidistant(test_mp, request): @@ -141,7 +282,7 @@ def test_per_period_equidistant(test_mp, request): 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 +299,24 @@ def test_cumulative_variable_periodlength(test_mp, request): 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]) + + # 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 +326,7 @@ def test_per_period_variable_periodlength(test_mp, request): 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,53 +344,156 @@ def test_custom_type_variable_periodlength(test_mp, request): 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]) - - -def test_price_duality(test_mp, request): - years = [2020, 2025, 2030, 2040, 2050] - for c in [0.25, 0.5, 0.75]: - # set up a scenario for cumulative constraints - scen = Scenario( - test_mp, MODEL, scenario=request.node.name + "_cum_many_tecs", version="new" - ) - 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" - ) - scen.commit("initialize test scenario") - scen.solve(quiet=True) - # set up a new scenario with emissions taxes - tax_scen = Scenario( - test_mp, MODEL, scenario=request.node.name + "_tax_many_tecs", version="new" - ) - model_setup(tax_scen, years, simple_tecs=False) - for y in years: - tax_scen.add_cat("year", y, y) + # 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) + + +# TODO as per Oliver's description: +# - Only need cumulative (no need to parametrize) +@pytest.mark.parametrize( + "bound, years", + [ + (10, [2010, 2020, 2025, 2030, 2035, 2040, 2050]), + (10, [2010, 2015, 2020, 2030, 2040, 2045, 2050]), + (40, [2010, 2020, 2025, 2030, 2035, 2040, 2050]), + (40, [2010, 2015, 2020, 2030, 2040, 2045, 2050]), + (75, [2010, 2020, 2025, 2030, 2035, 2040, 2050]), + (75, [2010, 2015, 2020, 2030, 2040, 2045, 2050]), + ], +) +def test_price_duality(test_mp, request, bound, years): + # set up a scenario for cumulative constraints + scen = Scenario( + test_mp, + MODEL, + scenario=f"{request.node.name}_cum_many_tecs", + version="new", + ) + model_setup(scen, years, simple_tecs=False) + bound_emission_base = dict( + node="World", type_emission="GHG", type_tec="all", value=bound, unit="tCO2" + ) + # We use a `cumulative` base + scen.add_cat("year", "cumulative", years) + scen.add_par( + "bound_emission", + make_df( + "bound_emission", + type_year="cumulative", + **bound_emission_base, + ), + ) + scen.commit("initialize test scenario") + scen.solve(quiet=True, **solve_args) + + # ---------------------------------------------------------- + # Run scenario with `tax_emission` based on `PRICE_EMISSION` + # from cumulative constraint scenario. + # ---------------------------------------------------------- + scen_tax = Scenario( + test_mp, + MODEL, + scenario=f"{request.node.name}_tax_many_tecs", + version="new", + ) + model_setup(scen_tax, years, simple_tecs=False) + # TODO Why do we do this? We applied the bound as cumulative before, but the tax + # should be per-period? + for year in years: + scen_tax.add_cat("year", year, year) + + # use emission prices from cumulative-constraint scenario as taxes + taxes = scen.var("PRICE_EMISSION").rename( + columns={"year": "type_year", "lvl": "value"} + ) + taxes["unit"] = "USD/tCO2" + # TODO Check: bound is on "World", taxes are on "node" + taxes["node"] = "node" + scen_tax.add_par("tax_emission", taxes) + scen_tax.commit("initialize test scenario for taxes") + scen_tax.solve(quiet=True, **solve_args) + # scen_tax.solve(**solve_args) + + print(scen.var("PRICE_EMISSION")) + print(scen_tax.var("PRICE_EMISSION")) + print(scen_tax.par("tax_emission")) + print(scen.var("EMISS")) + print(scen_tax.var("EMISS")) + + # check emissions are close between cumulative and tax scenarios + filters = {"node": "World"} + emiss = scen.var("EMISS", filters).set_index("year").lvl + emiss_tax = scen_tax.var("EMISS", filters).set_index("year").lvl + npt.assert_allclose(emiss, emiss_tax, rtol=0.05) + + # check "PRICE_EMISSION" is close between cumulative and tax scenarios + filters = {"node": "World"} + price_emission = scen.var("PRICE_EMISSION", filters).set_index("year").lvl + price_emission_tax = scen_tax.var("PRICE_EMISSION", filters).set_index("year").lvl + npt.assert_allclose(price_emission, price_emission_tax) + + # check "tax_emission" and "PRICE_EMISSION" are equal in tax-scenario + filters = {"node": "World"} + tax_emission = scen_tax.par("tax_emission", filters).set_index("type_year").value + price_emission_tax = scen_tax.var("PRICE_EMISSION", filters).set_index("year").lvl + npt.assert_allclose(tax_emission, price_emission_tax) + + # -------------------------------------------------------- + # Run scenario with annual-emission bound based on `EMISS` + # from cumulative constraint scenario. + # -------------------------------------------------------- + + scenario_period_bound = Scenario( + test_mp, + MODEL, + f"{request.node.name}_period_bound_many_tecs", + version="new", + ) + model_setup(scenario_period_bound, years, simple_tecs=False) + for year in years: + scenario_period_bound.add_cat("year", year, year) + + # use emissions from cumulative-constraint scenario as period-emission bounds + emiss_period_bound = ( + scen.var("EMISS", {"node": "World"}) + .rename(columns={"year": "type_year", "lvl": "value"}) + .drop("emission", axis=1) + ) + emiss_period_bound["type_emission"] = "GHG" + emiss_period_bound["unit"] = "tCO2" + # TODO: see above, _per_period_: bound_emission added for every single year. Does + # it work like this (for all at once), too? + scenario_period_bound.add_par("bound_emission", emiss_period_bound) + scenario_period_bound.commit("initialize test scenario for periodic emission bound") + scenario_period_bound.solve(quiet=True, **solve_args) + + # check -emissions are close between cumulative and yearly-bound scenarios + emiss_period_bound = ( + scenario_period_bound.var("EMISS", filters).set_index("year").lvl + ) + npt.assert_allclose(emiss, emiss_period_bound) - # use emission prices from cumulative-constraint scenario as taxes - taxes = scen.var("PRICE_EMISSION").rename( - columns={"year": "type_year", "lvl": "value"} - ) - taxes["unit"] = "USD/tCO2" - tax_scen.add_par("tax_emission", taxes) - tax_scen.commit("initialize test scenario for taxes") - tax_scen.solve(quiet=True) - - # check that emissions are close between cumulative and tax scenario - filters = {"node": "World"} - emiss = scen.var("EMISS", filters).set_index("year").lvl - emiss_tax = tax_scen.var("EMISS", filters).set_index("year").lvl - npt.assert_allclose(emiss, emiss_tax, rtol=0.20) + # check "PRICE_EMISSION" is close between cumulative- and yearly-bound scenarios + price_emission_period_bound = ( + scenario_period_bound.var("PRICE_EMISSION", filters).set_index("year").lvl + ) + npt.assert_allclose(price_emission, price_emission_period_bound) diff --git a/message_ix/tests/test_report.py b/message_ix/tests/test_report.py index 745cf5c6d..e4c19a8b1 100644 --- a/message_ix/tests/test_report.py +++ b/message_ix/tests/test_report.py @@ -70,11 +70,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") diff --git a/tutorial/westeros/westeros_emissions_taxes.ipynb b/tutorial/westeros/westeros_emissions_taxes.ipynb index fd532291e..f3ccfaf07 100644 --- a/tutorial/westeros/westeros_emissions_taxes.ipynb +++ b/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. " ] @@ -138,6 +138,8 @@ "metadata": {}, "outputs": [], "source": [ + "from message_ix.util import make_df\n", + "\n", "horizon = [700, 710, 720]\n", "\n", "bd_emission = message_ix.make_df(\n",