diff --git a/message_ix/model/MESSAGE-MACRO_run.gms b/message_ix/model/MESSAGE-MACRO_run.gms index 36cb49217..7c7060bd9 100644 --- a/message_ix/model/MESSAGE-MACRO_run.gms +++ b/message_ix/model/MESSAGE-MACRO_run.gms @@ -169,6 +169,8 @@ DISPLAY enestart, eneprice, total_cost ; * solve MACRO model * *----------------------------------------------------------------------------------------------------------------------* +total_cost(node_macro, year) = COST_NODAL_NET.L(node_macro, year) / 1000 ; + $INCLUDE MACRO/macro_solve.gms *----------------------------------------------------------------------------------------------------------------------* @@ -204,7 +206,7 @@ COST_NODAL_NET.L(node_macro,year) = * EMISS.L(node_macro,emission,type_tec,year) ) ; -GDP.L(node_macro,year) = (I.L(node_macro,year) + C.L(node_macro,year) + EC.L(node_macro,year)) * 1000 ; +GDP.L(node_macro,year) = (I.L(node_macro,year) + C.L(node_macro,year) + EC.L(node_macro,year) - ( trade_cost(node_macro, year) * 1E-6 ))*1000; * calculate convergence level (maximum absolute scaling factor minus 1 across all regions, sectors, and years) max_adjustment_pos = smax((node_macro,sector,year)$( NOT macro_base_period(year) AND demand_scale(node_macro,sector,year) > 1), diff --git a/message_ix/model/MESSAGE/model_solve.gms b/message_ix/model/MESSAGE/model_solve.gms index 93cdacd47..000c85eaa 100644 --- a/message_ix/model/MESSAGE/model_solve.gms +++ b/message_ix/model/MESSAGE/model_solve.gms @@ -14,7 +14,7 @@ if (%foresight% = 0, * This is the standard option; the GAMS global variable ``%foresight%=0`` by default. * * .. math:: -* \min_x \text{OBJ} = \sum_{y \in Y} \text{OBJ}_y(x_y) +* \min_x OBJ = \sum_{y \in Y} OBJ_y(x_y) *** * reset year in case it was set by MACRO to include the base year before @@ -64,6 +64,7 @@ EMISSION_CONSTRAINT.m(node,type_emission,type_tec,type_year)$( ( 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 ; %AUX_BOUNDS% AUX_ACT_BOUND_UP(node,tec,year_all,year_all2,mode,time)$( ACT.l(node,tec,year_all,year_all2,mode,time) > 0 AND @@ -82,10 +83,10 @@ else * Loop over :math:`\hat{y} \in Y`, solving * * .. math:: -* \min_x \ \text{OBJ} = \sum_{y \in \hat{Y}(\hat{y})} \text{OBJ}_y(x_y) \\ +* \min_x \ OBJ = \sum_{y \in \hat{Y}(\hat{y})} OBJ_y(x_y) \\ * \text{s.t. } x_{y'} = x_{y'}^* \quad \forall \ y' < y * -* where :math:`\hat{Y}(\hat{y}) = \{y \in Y | \ |\hat{y}| - |y| < \text{optimization_horizon} \}` and +* where :math:`\hat{Y}(\hat{y}) = \{y \in Y | \ |\hat{y}| - |y| < optimization\_horizon \}` and * :math:`x_{y'}^*` is the optimal value of :math:`x_{y'}` in iteration :math:`|y'|` of the iterative loop. * * The advantage of this implementation is that there is no need to 'store' the optimal values of all decision diff --git a/message_ix/models.py b/message_ix/models.py index 057f50e94..d411577f2 100644 --- a/message_ix/models.py +++ b/message_ix/models.py @@ -21,6 +21,8 @@ "lpmethod": 4, "threads": 4, "epopt": 1e-6, + "scaind": -1, + "memoryemphasis": "true", } #: Common dimension name abbreviations mapped to tuples with: @@ -176,7 +178,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 @@ -212,7 +214,7 @@ def run(self, scenario): optfile.write_text("\n".join(lines)) log.info(f"Use CPLEX options {self.cplex_opts}") - self.cplex_opts.update({"barcrossalg": 2}) + self.cplex_opts.update({"predual": 1}) optfile2 = Path(self.model_dir).joinpath("cplex.op2") lines2 = ("{} = {}".format(*kv) for kv in self.cplex_opts.items()) optfile2.write_text("\n".join(lines2)) @@ -560,6 +562,11 @@ def initialize(cls, scenario): "n type_emission type_tec y", "Emission price (derived from marginals of EMISSION_EQUIVALENCE constraint)", ) +var( + "PRICE_EMISSION_NEW", + "n type_emission type_tec y", + "TEMPORARY test for Emission price fix", +) var( "REL", "r nr yr", diff --git a/message_ix/tests/feature_price_emission.ipynb b/message_ix/tests/feature_price_emission.ipynb new file mode 100644 index 000000000..2826007f0 --- /dev/null +++ b/message_ix/tests/feature_price_emission.ipynb @@ -0,0 +1,327 @@ +{ + "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", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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 = 1\n", + "# cumulative_bound = 500\n", + "# cumulative_bound = 550\n", + "years = [2020, 2030, 2040, 2050]\n", + "# years = [2020, 2025, 2030, 2040, 2050]\n", + "# years = [2020, 2030, 2040, 2045, 2050]\n", + "\n", + "filters = {\"node\": \"World\"}\n", + "\n", + "model_setup(scen=scen, years=years, simple_tecs=False)" + ] + }, + { + "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.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": null, + "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.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": null, + "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.par(\"emission_factor\")" + ] + }, + { + "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_cumulative_bound.par(\"bound_emission\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scen_tax = Scenario(\n", + " mp,\n", + " MODEL,\n", + " scenario=\"tax_many_tecs\",\n", + " version=\"new\",\n", + ")\n", + "model_setup(scen_tax, years, simple_tecs=False)\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\"))\n", + "print(scenario_cumulative_bound.var(\"EMISS\"))\n", + "print(scen_tax.var(\"EMISS\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(scenario_cumulative_bound.var(\"ACT\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(scen_tax.var(\"ACT\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound.par(\"demand\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scenario_cumulative_bound.equ(\"EMISSION_EQUIVALENCE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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": null, + "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 8e885b0c4..6a655fd93 100644 --- a/message_ix/tests/test_feature_price_emission.py +++ b/message_ix/tests/test_feature_price_emission.py @@ -1,109 +1,143 @@ +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"]} +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") 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, - ), + make_df("emission_factor", technology="dirty_tec", emission="CO2", unit="tCO2"), ) # 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", - time="year", - unit="USD/GWa", - technology="clean_tec", - **common, - ), + make_df("var_cost", technology="clean_tec", time="year", unit="USD/GWa"), ) -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": 10, + "inv_cost": 500, + "fix_cost": 30, + "var_cost": 30, + "bound_activity_up": 100, + "lifetime": 20, + }, + "tec2": { + "emission_factor": -10, + "inv_cost": 1200, + "fix_cost": 10, + "var_cost": 0, + "bound_activity_up": 40, + "lifetime": 20, + }, + "tec3": { + "emission_factor": -12, + "inv_cost": 1300, + "fix_cost": 12, + "var_cost": 0, + "bound_activity_up": 30, + "lifetime": 20, + }, + "tec4": { + "emission_factor": -14, + "inv_cost": 1400, + "fix_cost": 14, + "var_cost": 0, + "bound_activity_up": 20, + "lifetime": 20, + }, + "tec5": { + "emission_factor": -16, + "inv_cost": 1500, + "fix_cost": 16, + "var_cost": 0, + "bound_activity_up": 10, + "lifetime": 20, + }, + } + year_df = scen.vintage_and_active_years() + vintage_years, act_years = year_df["year_vtg"], year_df["year_act"] + + 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("var_cost", tec_specs + ["year"], tecs[t][1], "USD/GWa") + scen.add_par("emission_factor", tec_specs + ["co2"], tecs[t][0], "tCO2") + scen.add_par( + "bound_activity_up", ["node", t, y, "mode", "year"], tecs[t][2], "GWa" + ) def test_no_constraint(test_mp, request): @@ -133,7 +167,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): @@ -230,40 +264,233 @@ def test_custom_type_variable_periodlength(test_mp, request): npt.assert_allclose(obs, exp) -def test_price_duality(test_mp, request): +@pytest.mark.parametrize("cumulative_bound", [0.25, 0.5, 0.75]) +def test_price_duality(test_mp, request, cumulative_bound): 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" - ) + scen = Scenario(test_mp, MODEL, "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" +def test_price_duality(test_mp, request, cumulative_bound, years, tag): + # 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" + ) + if cumulative: + scen.add_cat("year", "cumulative", years) + scen.add_par( + "bound_emission", + make_df( + "bound_emission", + type_year="cumulative", + **bound_emission_base, + ), ) - model_setup(tax_scen, years, simple_tecs=False) + else: for y in years: - tax_scen.add_cat("year", y, y) + scen.add_cat("year", y, y) + scen.add_par( + "bound_emission", + make_df("bound_emission", type_year=y, **bound_emission_base), + ) + scen.commit("initialize test scenario") + scen.solve(quiet=True, **solve_args) - # 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) + # ---------------------------------------------------------- + # 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) + 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" + 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) + + # -------------------------------------------------------- + # 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) + + # 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) + + +# idea: try running the same with the westeros scenario +# Let's not do that, though, because make_westeros has the westeros years hardcoded +# Instead, try to ensure this model setup follows westeros as a sanity check +# @pytest.mark.parametrize( +# "cumulative_bound, years", +# [ +# (400, [2020, 2030, 2040, 2050]), +# (400, [2020, 2025, 2030, 2040, 2050]), +# (500, [2020, 2030, 2040, 2050]), +# (500, [2020, 2025, 2030, 2040, 2050]), +# (600, [2020, 2030, 2040, 2050]), +# (600, [2020, 2025, 2030, 2040, 2050]), +# ], +# ) +# def test_price_duality_westeros(test_mp, request, cumulative_bound, years): +# baseline = make_westeros( +# mp=test_mp, model_horizon=years, emissions=True, request=request +# ) +# year_df = baseline.vintage_and_active_years() +# vintage_years, act_years = year_df["year_vtg"], year_df["year_act"] +# country = "Westeros" +# test_mp.add_unit("MtCO2") +# emission_factor = make_df( +# "emission_factor", +# node_loc=country, +# year_vtg=vintage_years, +# year_act=act_years, +# mode="standard", +# unit="tCO2/kWa", +# technology="coal_ppl", +# emission="CO2", +# value=7.4, +# ) + +# # Define cumulate emission bound and solve that Scenario +# scen_cumulative_bound = baseline.clone( +# MODEL, +# "emission_bound", +# "introducing an upper bound on emissions", +# keep_solution=False, +# ) +# scen_cumulative_bound.check_out() +# scen_cumulative_bound.add_set("emission", "CO2") +# scen_cumulative_bound.add_cat("emission", "GHG", "CO2") +# scen_cumulative_bound.add_par("emission_factor", emission_factor) +# scen_cumulative_bound.add_par( +# "bound_emission", +# [country, "GHG", "all", "cumulative"], +# value=cumulative_bound, +# unit="MtCO2", +# ) +# scen_cumulative_bound.commit( +# comment="Introducing emissions and setting an upper bound" +# ) +# scen_cumulative_bound.set_as_default() +# scen_cumulative_bound.solve(quiet=True, **solve_args) + +# # ---------------------------------------------------------- +# # Run scenario with `tax_emission` based on `PRICE_EMISSION` +# # from cumulative constraint scenario. +# # ---------------------------------------------------------- +# scen_tax = baseline.clone( +# MODEL, +# "tax_emission", +# "introducing a fixed tax on emissions", +# keep_solution=False, +# ) +# scen_tax.check_out() +# scen_cumulative_bound.add_set("emission", "CO2") +# scen_cumulative_bound.add_cat("emission", "GHG", "CO2") +# scen_cumulative_bound.add_par("emission_factor", emission_factor) + +# # use emission prices from cumulative-constraint scenario as taxes +# taxes = scen_cumulative_bound.var("PRICE_EMISSION").rename( +# columns={"year": "type_year", "lvl": "value"} +# ) +# taxes["unit"] = "USD/tCO2" +# 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_cumulative_bound.var("PRICE_EMISSION")) +# print(scen_tax.var("PRICE_EMISSION")) +# print(scen_tax.par("tax_emission")) +# print(scen_cumulative_bound.var("EMISS")) +# print(scen_tax.var("EMISS")) + +# # check emissions are close between cumulative and tax scenarios +# filters = {"node": "World"} +# emiss = scen_cumulative_bound.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_cumulative_bound.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) diff --git a/tutorial/westeros/westeros_emissions_taxes.ipynb b/tutorial/westeros/westeros_emissions_taxes.ipynb index fa4257cb5..429223a3b 100644 --- a/tutorial/westeros/westeros_emissions_taxes.ipynb +++ b/tutorial/westeros/westeros_emissions_taxes.ipynb @@ -16,6 +16,15 @@ "- You have run the tutorial on introducing emissions (``westeros_emissions_bounds.ipynb``)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, { "cell_type": "code", "execution_count": null, @@ -138,6 +147,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",