diff --git a/workflow/rules/solve_electricity.smk b/workflow/rules/solve_electricity.smk index e9dacff9..86a2ea1f 100644 --- a/workflow/rules/solve_electricity.smk +++ b/workflow/rules/solve_electricity.smk @@ -26,6 +26,7 @@ rule solve_network: safer_reeds="config/policy_constraints/reeds/prm_annual.csv", rps_reeds="config/policy_constraints/reeds/rps_fraction.csv", ces_reeds="config/policy_constraints/reeds/ces_fraction.csv", + interface_limits="config/policy_constraints/transmission_interface_limits.csv", pop_layout=pop_layout_input, ev_policy=ev_policy_input, output: diff --git a/workflow/rules/validate.smk b/workflow/rules/validate.smk index a6b968bb..1c12a2d1 100644 --- a/workflow/rules/validate.smk +++ b/workflow/rules/validate.smk @@ -9,6 +9,7 @@ rule solve_network_validation: safer_reeds="config/policy_constraints/reeds/prm_annual.csv", rps_reeds="config/policy_constraints/reeds/rps_fraction.csv", ces_reeds="config/policy_constraints/reeds/ces_fraction.csv", + interface_limits="config/policy_constraints/transmission_interface_limits.csv", output: network=RESULTS + "{interconnect}/networks/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}_{sector}_operations.nc", diff --git a/workflow/scripts/opts/interfaces.py b/workflow/scripts/opts/interfaces.py new file mode 100644 index 00000000..4dd59c2f --- /dev/null +++ b/workflow/scripts/opts/interfaces.py @@ -0,0 +1,103 @@ +"""Adds aggregate inter-regional transmission interface limits (RESOLVE/NARIS). + +An interface is a bundle of transmission paths between two groups of regions, +capped in aggregate rather than path-by-path. The limits are read from a CSV +with the columns ``interface, region_1, region_2, flow_12, flow_21``, where +``flow_12`` is the MW cap on flow out of ``region_1`` into ``region_2`` and +``flow_21`` the cap on flow in the opposite direction, for example:: + + interface,region_1,region_2,flow_12,flow_21,Notes + CAISO_Imports,"p9, p10, p11","p2, p5, p6, ...",9728,10208,RESOLVE + +The caps are applied to the import/export ``Link`` components created by +``add_extra_components.add_elec_imports_exports`` and are therefore a no-op +when ``electricity.imports``/``electricity.exports`` are disabled. +""" + +import logging + +import pandas as pd +import pypsa +from opts._helpers import get_region_buses + +logger = logging.getLogger(__name__) + +TRADE_CARRIERS = ("imports", "exports") + + +def _parse_regions(cell: str) -> list[str]: + """Split a comma separated region cell into a list of region names.""" + return [region.strip() for region in str(cell).split(",") if region.strip()] + + +def _boundary_links( + n: pypsa.Network, + inside_regions: list[str], + outside_regions: list[str], + direction: str, +) -> pd.Index: + """Get the trade links crossing an interface. + + Links are selected by bus membership and carrier, never by parsing link + names. Imports run from an external ``{zone}_imports`` bus into a bus inside + ``inside_regions``; exports run the other way into a ``{zone}_exports`` bus. + """ + if direction not in TRADE_CARRIERS: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + links = n.links[n.links.carrier == direction] + if links.empty: + return links.index + + # The external trade buses carry the *outside* zone name in their `country` + # field, which `get_region_buses` also matches on, so drop them here. + inside_buses = get_region_buses(n, inside_regions) + inside_buses = inside_buses[~inside_buses.carrier.isin(TRADE_CARRIERS)] + + external_names = {f"{zone}_{direction}" for zone in outside_regions} + external_buses = n.buses[ + (n.buses.carrier == direction) & (n.buses.index.isin(external_names) | n.buses.country.isin(outside_regions)) + ] + + if direction == "imports": + crossing = links.bus0.isin(external_buses.index) & links.bus1.isin(inside_buses.index) + else: + crossing = links.bus0.isin(inside_buses.index) & links.bus1.isin(external_buses.index) + + return links[crossing].index + + +def add_interface_transmission_limits(n: pypsa.Network, limits_csv_path: str) -> None: + """Cap the aggregate per-snapshot flow across each transmission interface. + + ``flow_21`` limits total imports into ``region_1``, ``flow_12`` total + exports out of it. Rows without any matching link are skipped. + + Note that only the import/export links are constrained. Region_2 entries + that are inside the network (e.g. `p8` in a California-only run is itself a + California zone) contribute no trade links, so internal AC lines such as + p8-p9 escape the cap. The resulting understatement is documented, not + corrected here. + """ + limits = pd.read_csv(limits_csv_path) + + for _, row in limits.iterrows(): + region_1 = _parse_regions(row.region_1) + region_2 = _parse_regions(row.region_2) + + for direction, cap in (("imports", row.flow_21), ("exports", row.flow_12)): + if pd.isna(cap): + continue + + links = _boundary_links(n, region_1, region_2, direction) + if links.empty: + logger.info(f"No {direction} links cross interface {row.interface}; skipping limit") + continue + + lhs = n.model["Link-p"].sel(name=links.tolist()).sum("name") + + n.model.add_constraints( + lhs <= float(cap), + name=f"interface_limit-{row.interface}-{direction}", + ) + logger.info(f"Added {direction} limit of {cap} MW on interface {row.interface} over {len(links)} links") diff --git a/workflow/scripts/solve_network.py b/workflow/scripts/solve_network.py index b28538c5..9858d665 100644 --- a/workflow/scripts/solve_network.py +++ b/workflow/scripts/solve_network.py @@ -37,6 +37,7 @@ from constants import HOURS_PER_YEAR from opts.bidirectional_link import add_bidirectional_link_constraints from opts.interchange import add_interchange_constraints +from opts.interfaces import add_interface_transmission_limits from opts.land import add_land_use_constraints from opts.policy import ( add_regional_co2limit, @@ -202,6 +203,10 @@ def extra_functionality(n, snapshots): if config["electricity"].get("exports", {}).get("volume_limit", False): add_interchange_constraints(n, config, "exports", sector_enabled) + # Apply aggregate interface transmission limits if configured + if config["model_topology"].get("interface_transmission_limits", False): + add_interface_transmission_limits(n, global_snakemake.input.interface_limits) + # Apply sector-specific constraints if sector is enabled if sector_enabled: # Heat pump constraints diff --git a/workflow/scripts/test/test_interfaces.py b/workflow/scripts/test/test_interfaces.py new file mode 100644 index 00000000..91b0f9c5 --- /dev/null +++ b/workflow/scripts/test/test_interfaces.py @@ -0,0 +1,287 @@ +""" +Test the aggregate transmission interface limits. + +This module contains tests for the RESOLVE/NARIS style interface constraints +applied to the electricity import/export links in PyPSA-USA. +""" + +import os +import sys + +import pandas as pd +import pypsa +import pytest + +sys.path.append(os.path.join(os.path.dirname(__file__), "..")) + +from _helpers import get_multiindex_snapshots +from opts.interfaces import ( + _boundary_links, + _parse_regions, + add_interface_transmission_limits, +) + +TOL = 1e-4 + +SHIPPED_LIMITS = os.path.join( + os.path.dirname(__file__), + "../../repo_data/config/policy_constraints/transmission_interface_limits.csv", +) + +# Fixtures + + +@pytest.fixture +def interface_network(): + """ + Build a small network with electricity import and export links. + + Mirrors the conventions of ``add_extra_components.add_elec_imports_exports``: + external buses are named ``{zone}_imports`` / ``{zone}_exports`` with the + matching carrier, import links run from the external bus into the model and + export links run the other way. + """ + n = pypsa.Network() + + n.snapshots = get_multiindex_snapshots( + sns_config={"start": "2030-01-01 00:00", "end": "2030-01-01 03:00", "inclusive": "both"}, + invest_periods=[2030], + ) + n.set_investment_periods(periods=[2030]) + + for carrier in ("AC", "gas", "imports", "exports"): + n.add("Carrier", carrier, co2_emissions=0) + + # Buses inside the model + n.add( + "Bus", + ["CA_Z1", "TX_Z1"], + carrier="AC", + country="US", + interconnect="western", + nerc_reg=["WECC", "WECC"], + reeds_state=["CA", "TX"], + reeds_zone=["CA_Z1", "TX_Z1"], + ) + + # External trade buses. `_add_import_export_buses` stamps the outside zone + # name onto `country`, which is what makes the carrier filter necessary. + n.add( + "Bus", + ["p2_imports", "p5_imports"], + carrier="imports", + country=["p2", "p5"], + interconnect="western", + ) + n.add( + "Bus", + "p2_exports", + carrier="exports", + country="p2", + interconnect="western", + ) + + # Trade links + n.add( + "Link", + "CA_Z1_p2_imports", + bus0="p2_imports", + bus1="CA_Z1", + carrier="imports", + p_nom=500, + marginal_cost=0, + ) + n.add( + "Link", + "CA_Z1_p5_imports", + bus0="p5_imports", + bus1="CA_Z1", + carrier="imports", + p_nom=500, + marginal_cost=0, + ) + # Decoy: an import link landing outside region_1 + n.add( + "Link", + "TX_Z1_p2_imports", + bus0="p2_imports", + bus1="TX_Z1", + carrier="imports", + p_nom=500, + marginal_cost=0, + ) + n.add( + "Link", + "CA_Z1_p2_exports", + bus0="CA_Z1", + bus1="p2_exports", + carrier="exports", + p_nom=500, + marginal_cost=0, + ) + # Decoy: an internal AC link between two in-model buses + n.add( + "Link", + "CA_Z1_TX_Z1", + bus0="CA_Z1", + bus1="TX_Z1", + carrier="AC", + p_nom=500, + ) + + # Generation and demand: imports are cheap, local gas is not + n.add("Generator", "import_p2", bus="p2_imports", carrier="imports", p_nom=1000, marginal_cost=1) + n.add("Generator", "import_p5", bus="p5_imports", carrier="imports", p_nom=1000, marginal_cost=1) + n.add("Generator", "gas_ca", bus="CA_Z1", carrier="gas", p_nom=1000, marginal_cost=50) + n.add("Generator", "gas_tx", bus="TX_Z1", carrier="gas", p_nom=1000, marginal_cost=50) + # Expensive backup behind the export bus, so an export cap stays feasible + n.add("Generator", "gas_p2", bus="p2_exports", carrier="gas", p_nom=1000, marginal_cost=100) + + n.add("Load", "load_ca", bus="CA_Z1", carrier="AC", p_set=pd.Series(300.0, index=n.snapshots)) + n.add("Load", "load_tx", bus="TX_Z1", carrier="AC", p_set=pd.Series(200.0, index=n.snapshots)) + n.add("Load", "load_p2", bus="p2_exports", carrier="exports", p_set=pd.Series(400.0, index=n.snapshots)) + + return n + + +def write_limits(tmp_path, rows): + """Write an interface limits CSV and return its path.""" + path = tmp_path / "transmission_interface_limits.csv" + pd.DataFrame(rows).to_csv(path, index=False) + return str(path) + + +# Tests + + +def test_parse_regions_handles_whitespace(): + assert _parse_regions("p9, p10,p11 ") == ["p9", "p10", "p11"] + assert _parse_regions("p9") == ["p9"] + assert _parse_regions("") == [] + + +def test_boundary_links_selects_only_crossing_links(interface_network): + n = interface_network + + imports = _boundary_links(n, ["CA_Z1"], ["p2", "p5"], "imports") + assert sorted(imports) == ["CA_Z1_p2_imports", "CA_Z1_p5_imports"] + + exports = _boundary_links(n, ["CA_Z1"], ["p2", "p5"], "exports") + assert sorted(exports) == ["CA_Z1_p2_exports"] + + # Only the named outside zones count + assert sorted(_boundary_links(n, ["CA_Z1"], ["p5"], "imports")) == ["CA_Z1_p5_imports"] + + # Region_1 can be given as any of the labels get_region_buses matches on + assert sorted(_boundary_links(n, ["CA"], ["p2", "p5"], "imports")) == [ + "CA_Z1_p2_imports", + "CA_Z1_p5_imports", + ] + + with pytest.raises(ValueError): + _boundary_links(n, ["CA_Z1"], ["p2"], "both") + + +def test_no_matching_links_is_a_noop_not_an_error(interface_network, tmp_path): + n = interface_network + limits = write_limits( + tmp_path, + [ + # Neither region exists in this network + {"interface": "NW_SW", "region_1": "p30", "region_2": "p33", "flow_12": 100, "flow_21": 100}, + ], + ) + + def extra_functionality(n, sns): + add_interface_transmission_limits(n, limits) + # the shipped RESOLVE table names ReEDS zones absent from this network + add_interface_transmission_limits(n, SHIPPED_LIMITS) + + n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) + + assert not [c for c in n.model.constraints if c.startswith("interface_limit-")] + + +def test_import_cap_binds(interface_network, tmp_path): + n = interface_network + cap = 100.0 + limits = write_limits( + tmp_path, + [ + { + "interface": "CAISO_Imports", + "region_1": "CA_Z1", + "region_2": "p2, p5", + "flow_12": 1e6, + "flow_21": cap, + }, + ], + ) + + def extra_functionality(n, sns): + add_interface_transmission_limits(n, limits) + + n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) + + assert "interface_limit-CAISO_Imports-imports" in n.model.constraints + + flow = n.links_t.p0[["CA_Z1_p2_imports", "CA_Z1_p5_imports"]].sum(axis=1) + assert (flow <= cap + TOL).all() + assert flow.max() >= cap - TOL, "import cap should bind in at least one snapshot" + + # The decoy import link into TX is outside the interface and stays free + assert n.links_t.p0["TX_Z1_p2_imports"].max() > cap + TOL + + +def test_export_cap_uses_flow_12(interface_network, tmp_path): + n = interface_network + cap = 150.0 + limits = write_limits( + tmp_path, + [ + { + "interface": "CAISO_Exports", + "region_1": "CA_Z1", + "region_2": "p2, p5", + "flow_12": cap, + "flow_21": 1e6, + }, + ], + ) + + def extra_functionality(n, sns): + add_interface_transmission_limits(n, limits) + + n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) + + assert "interface_limit-CAISO_Exports-exports" in n.model.constraints + + flow = n.links_t.p0["CA_Z1_p2_exports"] + assert (flow <= cap + TOL).all() + assert flow.max() >= cap - TOL, "export cap should bind in at least one snapshot" + + +def test_disabled_flag_adds_no_constraints(interface_network, tmp_path): + n = interface_network + limits = write_limits( + tmp_path, + [ + { + "interface": "CAISO_Imports", + "region_1": "CA_Z1", + "region_2": "p2, p5", + "flow_12": 100, + "flow_21": 100, + }, + ], + ) + config = {"model_topology": {"interface_transmission_limits": False}} + + def extra_functionality(n, sns): + # mirrors the gate in solve_network.extra_functionality + if config["model_topology"].get("interface_transmission_limits", False): + add_interface_transmission_limits(n, limits) + + n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) + + assert not [c for c in n.model.constraints if c.startswith("interface_limit-")]