From f5eb10f0f1bde3bde17af13aaf38f9262a1cf787 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:09:46 -0700 Subject: [PATCH 01/23] Add SERVM load-allocation weights artifact (CPUC CA demand, PR 1) Builds a fractional per-(SERVM region, cluster bus) allocation table so CPUC SERVM demand for the six California load regions can be disaggregated onto the simplified network. The shares are derived by composing the two busmaps (base bus -> substation -> {simpl} cluster bus) against elec_base_network.nc rather than reading a bus column of the clustered network: aggregate_to_substations drops `balancing_area` from the buses, and a cluster can straddle two SERVM regions anyway (LA County holds both LDWP and CISO-SCE buses). The base network is the last stage that still carries both `balancing_area` and the population-based `load_weight`, so the composition recovers where each base bus's weight lands. A straddling cluster therefore appears once per region it overlaps, each row carrying only that region's share, and the factors sum to 1.0 within a region. Balancing areas that are deliberately out of scope (CISO-VEA, Nevada footprint) keep a row with an empty region in servm_region_map.csv, so an unknown balancing area introduced by upstream relabeling still hard-fails while the excluded one is dropped with a logged load share. Co-Authored-By: Claude Fable 5 --- workflow/repo_data/CPUC/servm_region_map.csv | 9 + workflow/rules/build_electricity.smk | 20 ++ workflow/scripts/build_servm_load_weights.py | 206 ++++++++++++++++++ .../test/test_build_servm_load_weights.py | 188 ++++++++++++++++ 4 files changed, 423 insertions(+) create mode 100644 workflow/repo_data/CPUC/servm_region_map.csv create mode 100644 workflow/scripts/build_servm_load_weights.py create mode 100644 workflow/scripts/test/test_build_servm_load_weights.py diff --git a/workflow/repo_data/CPUC/servm_region_map.csv b/workflow/repo_data/CPUC/servm_region_map.csv new file mode 100644 index 00000000..32ac294f --- /dev/null +++ b/workflow/repo_data/CPUC/servm_region_map.csv @@ -0,0 +1,9 @@ +balancing_area,servm_region,note +CISO-PGAE,PGE, +CISO-SCE,SCE, +CISO-SDGE,SDGE, +IID,IID, +LDWP,LADWP, +BANC,NCNC, +TIDC,NCNC, +CISO-VEA,,Nevada footprint; excluded from CA-only networks diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk index daf03dc6..ecd5fa0d 100644 --- a/workflow/rules/build_electricity.smk +++ b/workflow/rules/build_electricity.smk @@ -847,6 +847,26 @@ rule cluster_resources: "../scripts/cluster_simpl.py" +rule build_servm_load_weights: + input: + network=NETWORKS + "{interconnect}/elec_base_network.nc", + busmap_b=BUSMAPS + "{interconnect}/busmap_b.csv", + busmap_s=BUSMAPS + "{interconnect}/busmap_s{simpl}.csv", + region_map="repo_data/CPUC/servm_region_map.csv", + output: + weights=DEMAND + "{interconnect}/servm_load_weights_s{simpl}.csv", + log: + LOGS + "{interconnect}/build_servm_load_weights_s{simpl}.log", + threads: 1 + resources: + mem_mb=lambda wildcards, input, attempt: (input.size // 200000) * attempt * 2, + walltime=config_provider( + "walltime", "build_servm_load_weights", default="00:20:00" + ), + script: + "../scripts/build_servm_load_weights.py" + + rule cluster_network: params: cluster_network=config_provider("clustering", "cluster_network"), diff --git a/workflow/scripts/build_servm_load_weights.py b/workflow/scripts/build_servm_load_weights.py new file mode 100644 index 00000000..3073f67c --- /dev/null +++ b/workflow/scripts/build_servm_load_weights.py @@ -0,0 +1,206 @@ +# BY PyPSA-USA Authors +"""Fractional load-allocation factors mapping SERVM load regions onto clustered buses. + +The CPUC SERVM demand dataset is published for six California load regions (IID, +LADWP, NCNC, PGE, SCE, SDGE). Attaching it to the model needs, per region, the +share of that region's demand that lands on each bus of the simplified network. + +The share cannot be read off the network the demand is attached to: +``aggregate_to_substations`` drops ``balancing_area`` from the buses, and a single +substation/cluster bus may straddle two SERVM regions anyway (LA County holds both +LDWP and CISO-SCE buses). Both facts are handled by going back to +``elec_base_network.nc`` — the last network that still carries ``balancing_area`` +and the population-based ``load_weight`` per bus — and composing the two busmaps +(base bus -> substation -> ``{simpl}`` cluster bus) to find where each base bus's +weight ends up. A straddling cluster therefore appears once per SERVM region it +overlaps, each row carrying only that region's share. + +The resulting long-format table (``bus``, ``servm_region``, ``laf``) sums to 1.0 +within every region, so it is scale-free: it disaggregates whatever regional +demand the downstream rule reads. +""" + +import logging + +import pandas as pd +import pypsa +from _helpers import configure_logging + +logger = logging.getLogger(__name__) + +SERVM_REGIONS = ("IID", "LADWP", "NCNC", "PGE", "SCE", "SDGE") + + +def load_servm_region_map(path: str) -> pd.Series: + """ + Balancing-area -> SERVM region lookup read from ``servm_region_map.csv``. + + Balancing areas that are deliberately out of scope (CISO-VEA, whose footprint + is in Nevada) keep their row with a null region rather than being dropped, so + that :func:`build_load_weights` can tell a known-but-excluded balancing area + from one that upstream relabeling has made unknown. + """ + df = pd.read_csv(path, dtype=str) + ba = df["balancing_area"].str.strip() + region = df["servm_region"].str.strip().replace("", None) + + duplicated = ba[ba.duplicated()].unique() + if len(duplicated): + raise ValueError(f"Duplicate balancing areas in {path}: {list(duplicated)}") + + return pd.Series(region.to_numpy(), index=pd.Index(ba, name="balancing_area"), name="servm_region") + + +def compose_busmaps(busmap_b: pd.Series, busmap_s: pd.Series) -> pd.Series: + """ + Chain the base-bus -> substation and substation -> cluster busmaps. + + ``busmap_b`` is written by ``aggregate_to_substations`` (index: base network + bus, value: ``sub_id``) and ``busmap_s`` by ``cluster_simpl`` (index: + ``sub_id``, value: ``cluster_bus``). ``cluster_simpl`` runs for every + ``{simpl}`` value, emitting an identity map when the wildcard is empty, so the + second leg is always available. Substation ids round-trip through CSV as + integers, hence the explicit string normalization on both keys and values. + + Base buses whose substation is absent from ``busmap_s`` map to NaN; whether + that is fatal depends on whether they carry demand weight, which is decided in + :func:`build_load_weights`. + """ + busmap_b = busmap_b.astype(str) + busmap_b.index = busmap_b.index.astype(str) + busmap_s = busmap_s.astype(str) + busmap_s.index = busmap_s.index.astype(str) + + composed = busmap_b.map(busmap_s) + composed.name = "cluster_bus" + composed.index.name = "bus" + + unmapped = int(composed.isna().sum()) + if unmapped: + logger.info( + "%d of %d base buses have no cluster bus (substation missing from the simpl busmap).", + unmapped, + len(composed), + ) + return composed + + +def build_load_weights( + buses: pd.DataFrame, + region_map: pd.Series, + busmap: pd.Series, +) -> pd.DataFrame: + """ + Per-(SERVM region, cluster bus) load-allocation factors in long format. + + ``buses`` are the base-network buses, which must still carry + ``balancing_area`` and ``load_weight``. Returns columns ``bus``, + ``servm_region``, ``laf`` with ``laf`` summing to 1.0 within each region. + """ + missing_cols = {"balancing_area", "load_weight"}.difference(buses.columns) + if missing_cols: + raise ValueError( + f"Base network buses are missing {sorted(missing_cols)}. SERVM weights must be built from " + "elec_base_network.nc, the last network that carries the balancing area and load weight.", + ) + + # the index name is dropped so the added "bus" column cannot collide with it + df = pd.DataFrame( + { + "balancing_area": buses.balancing_area.fillna("").astype(str).str.strip(), + "load_weight": pd.to_numeric(buses.load_weight, errors="coerce").fillna(0.0), + }, + ).rename_axis(None) + weighted = df.load_weight > 0 + + blank_ba = weighted & (df.balancing_area == "") + if blank_ba.any(): + offenders = df.index[blank_ba] + raise ValueError( + f"{blank_ba.sum()} buses carry demand weight but have no balancing area, so their load " + f"cannot be assigned to a SERVM region: {list(offenders[:10])}", + ) + + unknown_ba = weighted & ~df.balancing_area.isin(region_map.index) + if unknown_ba.any(): + raise ValueError( + f"Balancing areas carrying demand weight are absent from the SERVM region map: " + f"{sorted(df.balancing_area[unknown_ba].unique())}. Add them to " + "repo_data/CPUC/servm_region_map.csv (with an empty region to exclude them).", + ) + + excluded_bas = region_map.index[region_map.isna()] + excluded = df.balancing_area.isin(excluded_bas) & weighted + total_weight = df.load_weight.sum() + if excluded.any() and total_weight > 0: + logger.info( + "Dropping %d buses in balancing areas excluded from SERVM (%s), carrying %.3f%% of total bus load weight.", + int(excluded.sum()), + ", ".join(sorted(df.balancing_area[excluded].unique())), + 100 * df.load_weight[excluded].sum() / total_weight, + ) + + df["servm_region"] = df.balancing_area.map(region_map) + df = df[df.servm_region.notna() & weighted].copy() + + df["bus"] = df.index.to_series().astype(str).map(busmap) + unmapped = df.bus.isna() + if unmapped.any(): + offenders = df.index[unmapped] + raise ValueError( + f"{unmapped.sum()} buses carrying SERVM demand weight are absent from the composed busmap: " + f"{list(offenders[:10])}", + ) + + weights = df.groupby(["servm_region", "bus"], as_index=False).load_weight.sum() + weights["laf"] = weights.load_weight / weights.groupby("servm_region").load_weight.transform("sum") + weights = weights[["bus", "servm_region", "laf"]].sort_values(["servm_region", "bus"]).reset_index(drop=True) + + if weights.empty: + raise ValueError( + "No bus carries load weight for any SERVM region. Check that the network covers California " + "and that its balancing areas match repo_data/CPUC/servm_region_map.csv.", + ) + + region_sums = weights.groupby("servm_region").laf.sum() + off = region_sums[(region_sums - 1.0).abs() > 1e-9] + if len(off): + raise ValueError(f"Load allocation factors do not sum to 1 per SERVM region: {off.to_dict()}") + + absent = [r for r in SERVM_REGIONS if r not in region_sums.index] + if absent: + logger.warning( + "No buses found for SERVM regions %s; their demand cannot be allocated in this network.", + ", ".join(absent), + ) + + return weights + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "build_servm_load_weights", + interconnect="western", + simpl="", + ) + configure_logging(snakemake) + + n = pypsa.Network(snakemake.input.network) + + busmap_b = pd.read_csv(snakemake.input.busmap_b, index_col=0, dtype=str).iloc[:, 0] + busmap_s = pd.read_csv(snakemake.input.busmap_s, index_col=0, dtype=str).iloc[:, 0] + busmap = compose_busmaps(busmap_b, busmap_s) + + region_map = load_servm_region_map(snakemake.input.region_map) + weights = build_load_weights(n.buses, region_map, busmap) + + logger.info( + "Built %d (SERVM region, bus) allocation weights across %d regions and %d buses.", + len(weights), + weights.servm_region.nunique(), + weights.bus.nunique(), + ) + weights.to_csv(snakemake.output.weights, index=False) diff --git a/workflow/scripts/test/test_build_servm_load_weights.py b/workflow/scripts/test/test_build_servm_load_weights.py new file mode 100644 index 00000000..d2e7f2f4 --- /dev/null +++ b/workflow/scripts/test/test_build_servm_load_weights.py @@ -0,0 +1,188 @@ +"""Unit tests for build_servm_load_weights helpers.""" + +import os +import sys + +import pandas as pd +import pytest + +sys.path.append(os.path.join(os.path.dirname(__file__), "..")) + +from build_servm_load_weights import ( + build_load_weights, + compose_busmaps, + load_servm_region_map, +) + + +@pytest.fixture +def region_map(): + """Mirror of repo_data/CPUC/servm_region_map.csv, CISO-VEA excluded.""" + return pd.Series( + { + "CISO-PGAE": "PGE", + "CISO-SCE": "SCE", + "CISO-SDGE": "SDGE", + "IID": "IID", + "LDWP": "LADWP", + "BANC": "NCNC", + "TIDC": "NCNC", + "CISO-VEA": None, + }, + name="servm_region", + ).rename_axis("balancing_area") + + +def make_buses(records): + """Base-network buses from (bus, balancing_area, load_weight) triples.""" + df = pd.DataFrame(records, columns=["bus", "balancing_area", "load_weight"]) + return df.set_index("bus") + + +def test_load_region_map_keeps_excluded_ba_as_null(tmp_path): + path = tmp_path / "servm_region_map.csv" + path.write_text( + "balancing_area,servm_region,note\nCISO-PGAE,PGE,\nCISO-VEA,,Nevada footprint\n", + ) + mapping = load_servm_region_map(str(path)) + assert mapping["CISO-PGAE"] == "PGE" + # excluded, not unknown: the row survives so the BA stays distinguishable + assert "CISO-VEA" in mapping.index + assert pd.isna(mapping["CISO-VEA"]) + + +def test_compose_busmaps_chains_base_to_cluster(): + busmap_b = pd.Series({"b1": "10", "b2": "10", "b3": "20"}, name="sub_id") + busmap_s = pd.Series({10: "c1", 20: "c2"}, name="cluster_bus") + + composed = compose_busmaps(busmap_b, busmap_s) + + assert composed.to_dict() == {"b1": "c1", "b2": "c1", "b3": "c2"} + + +def test_laf_sums_to_one_per_region(region_map): + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", "CISO-PGAE", 10.0), + ("b3", "CISO-SCE", 60.0), + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c2", "b3": "c3"}) + + weights = build_load_weights(buses, region_map, busmap) + + assert set(weights.columns) == {"bus", "servm_region", "laf"} + per_region = weights.groupby("servm_region").laf.sum() + assert per_region.round(12).eq(1.0).all() + pge = weights[weights.servm_region == "PGE"].set_index("bus").laf + assert pge["c1"] == pytest.approx(0.75) + assert pge["c2"] == pytest.approx(0.25) + + +def test_straddling_cluster_splits_across_regions(region_map): + """A cluster spanning LADWP and SCE gets one row per region, each at laf 1.""" + buses = make_buses( + [ + ("b1", "LDWP", 30.0), + ("b2", "CISO-SCE", 70.0), + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c1"}) + + weights = build_load_weights(buses, region_map, busmap) + + assert len(weights) == 2 + assert set(weights.bus) == {"c1"} + assert weights.set_index("servm_region").laf.to_dict() == {"LADWP": 1.0, "SCE": 1.0} + + +def test_unmapped_ba_with_load_raises(region_map): + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", "CISO-PGE", 20.0), # hypothetical upstream relabeling + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c2"}) + + with pytest.raises(ValueError, match="CISO-PGE"): + build_load_weights(buses, region_map, busmap) + + +def test_blank_ba_with_load_raises(region_map): + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", None, 20.0), + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c2"}) + + with pytest.raises(ValueError, match="b2"): + build_load_weights(buses, region_map, busmap) + + +def test_blank_ba_without_load_is_dropped(region_map): + """Offshore buses carry no balancing area but also no weight.""" + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", "Offshore", 0.0), + ("b3", None, 0.0), + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c1", "b3": "c1"}) + + weights = build_load_weights(buses, region_map, busmap) + + assert weights.laf.tolist() == [1.0] + + +def test_vea_dropped_without_error(region_map): + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", "CISO-PGAE", 10.0), + ("b3", "CISO-VEA", 60.0), + ], + ) + busmap = pd.Series({"b1": "c1", "b2": "c2", "b3": "c3"}) + + weights = build_load_weights(buses, region_map, busmap) + + assert "c3" not in set(weights.bus) + # the excluded weight leaves PGE's internal shares untouched + assert weights.set_index("bus").laf.to_dict() == {"c1": pytest.approx(0.75), "c2": pytest.approx(0.25)} + + +def test_missing_region_warns_not_raises(region_map, caplog): + buses = make_buses([("b1", "CISO-PGAE", 30.0)]) + busmap = pd.Series({"b1": "c1"}) + + with caplog.at_level("WARNING"): + weights = build_load_weights(buses, region_map, busmap) + + assert weights.servm_region.tolist() == ["PGE"] + assert "IID" in caplog.text + + +def test_bus_missing_from_busmap_raises(region_map): + buses = make_buses( + [ + ("b1", "CISO-PGAE", 30.0), + ("b2", "CISO-SCE", 20.0), + ], + ) + busmap = pd.Series({"b1": "c1"}) + + with pytest.raises(ValueError, match="b2"): + build_load_weights(buses, region_map, busmap) + + +def test_empty_result_raises(region_map): + buses = make_buses([("b1", "CISO-VEA", 30.0)]) + busmap = pd.Series({"b1": "c1"}) + + with pytest.raises(ValueError, match="No bus carries load weight"): + build_load_weights(buses, region_map, busmap) From 907906d458b9e97bfe84d57a2a52d4b24fbe4143 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:10:33 -0700 Subject: [PATCH 02/23] Wire aggregate interface transmission limits (RESOLVE CAISO import cap, PR 3) Rewires the two dead config keys `model_topology.interface_transmission_limits` and `electricity.transmission_interface_limits`: the RESOLVE interface table is now applied as a per-snapshot MW cap on the aggregate flow across each interface. Scope is the electricity import/export Links added by add_extra_components, so the constraint is inert when trade is disabled. Region_2 entries that are inside the network contribute no trade links, so internal AC lines escape the cap; that understatement is documented in the module, not corrected here. Co-Authored-By: Claude Fable 5 --- workflow/rules/solve_electricity.smk | 1 + workflow/rules/validate.smk | 1 + workflow/scripts/opts/interfaces.py | 103 ++++++++ workflow/scripts/solve_network.py | 5 + workflow/scripts/test/test_interfaces.py | 287 +++++++++++++++++++++++ 5 files changed, 397 insertions(+) create mode 100644 workflow/scripts/opts/interfaces.py create mode 100644 workflow/scripts/test/test_interfaces.py 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-")] From 2f0b94835762b52ed4c9fb0b7829095ef6984e10 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:26:57 -0700 Subject: [PATCH 03/23] Add CPUC SERVM demand profile for California (PR 2) Adds `electricity.demand.profile: servm`, wiring CPUC SERVM 2026 hourly load for the six California load regions (IID, LADWP, NCNC, PGE, SCE, SDGE) through the existing read/write demand strategies. Retrieve: `retrieve_cpuc_servm_load` pulls one ~118 MB CSV per forecast year (2026-2045) from files.cpuc.ca.gov, and `retrieve_cpuc_baseline_generators` pulls the CAISO baseline generator list, both through a thin `retrieve_cpuc_data.py` downloader patterned on `retrieve_eer_data.py`. ReadServm parses the three-row header positionally: the seventh index column ("Hour of Day") carries the stray ('Region', 'Unit Type') labels, so the calendar block cannot be found by testing the upper header levels for blankness. Every published component is kept on the `subsector` index level so the zonal artifact stays component-resolved; only `Net Load` reaches the model, selected by a new `ReadStrategy.default_subsector` hook that `Context` applies so `main()` stays profile-agnostic. Files are indexed by the forecast year parsed from each basename rather than by list order. Strips are fixed PST (UTC-8, no DST -- verified against the BTMPV solar-noon centroid, Dec 12.52 vs Jul 12.68) and are rolled to UTC with the same `np.roll` convention ReadEer uses for CST. Two calendar misalignments are accepted and documented on `ReadServm._assign_snapshots`: 1. SERVM lays its hours on a synthetic Monday-start calendar, so weekday and weekend hours do not line up with the real weekdays of the planning horizon. 2. For a leap *weather* year the strip contains February 29 and omits December 31, while the model snapshots do the opposite, so every hour after February lands one calendar day earlier than it sat in the source file. The 8760-hour strip is therefore mapped positionally onto the network's own per-period snapshots rather than onto a synthesised `date_range`, which would run a day short of December 31 for the leap planning horizons (2028/2032/2040) because `get_snapshots` drops February 29 from them. WriteServm disaggregates by matrix product against the PR-1 weights table pivoted to (region x bus), rather than the base class's one-zone-per-bus mapping: a cluster bus can straddle two SERVM regions (LA County holds both LDWP and CISO-SCE buses), and the product gives such a bus the sum of its share of every region it overlaps. `build_electrical_demand` gains a second output, the reader's zonal component-resolved demand written before disaggregation. It is produced for every profile (for efs/eer it is simply the single subsector="all" slice), and is captured off the `Context` rather than re-read. Config: `electricity.demand.scenario.servm_weather_years` selects the weather year from the stacked 2000-2024 record. A single-entry list is deterministic; multiple entries raise NotImplementedError pending stochastic scenarios. A value differing from the top-level `renewable_weather_years` logs a warning. Co-Authored-By: Claude Fable 5 --- docs/source/configtables/electricity.csv | 3 +- workflow/repo_data/config/config.default.yaml | 11 +- workflow/rules/build_electricity.smk | 39 +- workflow/rules/retrieve.smk | 34 ++ workflow/scripts/build_demand.py | 491 +++++++++++++++++- workflow/scripts/retrieve_cpuc_data.py | 24 + .../scripts/test/test_build_demand_servm.py | 418 +++++++++++++++ .../test/test_build_demand_servm_write.py | 139 +++++ 8 files changed, 1148 insertions(+), 11 deletions(-) create mode 100644 workflow/scripts/retrieve_cpuc_data.py create mode 100644 workflow/scripts/test/test_build_demand_servm.py create mode 100644 workflow/scripts/test/test_build_demand_servm_write.py diff --git a/docs/source/configtables/electricity.csv b/docs/source/configtables/electricity.csv index af570d63..50eb3437 100644 --- a/docs/source/configtables/electricity.csv +++ b/docs/source/configtables/electricity.csv @@ -33,12 +33,13 @@ gaslimit,MWh thermal,float,"Cap on annual gas-fired primary energy from gas carr ,,, demand:,,, -- bus_allocation,--,"One of {``population``, ``breakthrough``}","How zone-level demand is distributed to individual buses. ``population`` (default) weights buses by 2020 Decennial Census county populations (split evenly across each county's substations, then each substation's buses). ``breakthrough`` uses the legacy nominal-demand column (``Pd``) from the 2016-vintage Breakthrough Energy grid model." --- profile,--,"One of {``efs``, ``eia``, ``eer``}","Datasource for electrical load data. ``EFS`` pulls future state level electrical demand data. ``EIA`` pulls historical balancing level electrical demand data. ``EER`` pulls future state-level profiles from the EER dataset; when selected, ``planning_horizons`` must be one of 2021, 2025, 2030, 2035, 2040, 2045, or 2050 and ``renewable_weather_years`` must contain exactly one year from 2007-2013 or 2016-2023." +-- profile,--,"One of {``efs``, ``eia``, ``eer``, ``servm``}","Datasource for electrical load data. ``EFS`` pulls future state level electrical demand data. ``EIA`` pulls historical balancing level electrical demand data. ``EER`` pulls future state-level profiles from the EER dataset; when selected, ``planning_horizons`` must be one of 2021, 2025, 2030, 2035, 2040, 2045, or 2050 and ``renewable_weather_years`` must contain exactly one year from 2007-2013 or 2016-2023. ``SERVM`` pulls CPUC SERVM hourly load for the six California load regions (California models only); when selected, ``planning_horizons`` must be one of 2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, or 2045." -- scenario:,,, -- -- efs_case,--,"One of {``reference``, ``medium``, ``high``}",(UNDER DEVELOPMENT) Extracts EFS data according to level of adoption -- -- efs_speed,--,"One of {``slow``, ``moderate``, ``fast``}",(UNDER DEVELOPMENT) Extracts EFS data according to speed of electrification -- -- eer_file,--,"One of {``demand_EER2025_100by2050.h5``, ``demand_EER2025_Baseline_AEO2023.h5``, ``demand_EER2025_IRAlow.h5``}",Selects the EER demand dataset file to download and use when ``profile`` is ``eer``. -- -- aeo,--,One of the AEO scenarios `here `_,(UNDER DEVELOPMENT) Scales future demand according to the AEO scenario +-- -- servm_weather_years,--,List with exactly one year from 2000-2024,"Weather year drawn from the stacked SERVM record when ``profile`` is ``servm``. Multiple entries are reserved for stochastic scenarios and currently raise ``NotImplementedError``. Keep equal to the top-level ``renewable_weather_years`` so load and renewable profiles share a weather year; a mismatch logs a warning." ,,, demand_response:,,,Settings to activate and configure demand response -- shift,per_unit,"float {0 <=, >= 1} or 'inf'",Allowable load to be shifted per snapshot. Set to 0 to turn off demand response. Set to 'inf' to not enforce capacity limits. diff --git a/workflow/repo_data/config/config.default.yaml b/workflow/repo_data/config/config.default.yaml index 74a5d984..99fcc23d 100644 --- a/workflow/repo_data/config/config.default.yaml +++ b/workflow/repo_data/config/config.default.yaml @@ -163,13 +163,20 @@ electricity: # ---------------- Demand ---------------- demand: - profile: efs # demand time series source; efs (EIA-EFS) | eia (historical actuals) | eer (EER scenarios) + profile: efs # demand time series source; efs (EIA-EFS) | eia (historical actuals) | eer (EER scenarios) | servm (CPUC SERVM, California only) bus_allocation: population # per-bus demand weight; population (2020 census counties) | breakthrough (legacy BE Pd) - scenario: # EFS/EER scenario knobs (ignored if profile=eia) + scenario: # EFS/EER/SERVM scenario knobs (ignored if profile=eia) efs_case: reference # reference | medium | high efs_speed: moderate # slow | moderate | rapid eer_file: demand_EER2025_100by2050.h5 # EER h5 (profile=eer): demand_EER2025_100by2050 | demand_EER2025_Baseline_AEO2023 | demand_EER2025_IRAlow aeo: reference # AEO (EIA Annual Energy Outlook) scaling case; reference | high | low + # ---- profile: servm (CPUC SERVM 2026, California regions) ---- + # Weather year drawn from the stacked 2000-2024 SERVM record. Single-entry + # list = deterministic; multiple entries are reserved for stochastic + # scenarios (phase 3) and currently raise NotImplementedError. + # RECOMMENDED: keep equal to top-level `renewable_weather_years` so load + # and renewable profiles share a weather year (mismatch logs a warning). + servm_weather_years: [2019] demand_response: # price-responsive shiftable load; 0 = disabled shift: 0 # fraction of hourly load that can shift in time diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk index ecd5fa0d..37dec6fc 100644 --- a/workflow/rules/build_electricity.smk +++ b/workflow/rules/build_electricity.smk @@ -315,6 +315,15 @@ def eer_demand_file(): return filename +SERVM_LOAD_FILE = "cpuc/servm/HourlyLoad_CA_Regions_V2025E_2224_Mon_{year}.csv" + + +def servm_demand_files(): + """One CPUC SERVM hourly-load file per planning horizon.""" + horizons = sorted(set(config["scenario"]["planning_horizons"])) + return [DATA + SERVM_LOAD_FILE.format(year=year) for year in horizons] + + def demand_raw_data(wildcards): # get profile to use end_use = wildcards.end_use @@ -346,6 +355,8 @@ def demand_raw_data(wildcards): return DATA + f"nrel_efs/EFSLoadProfile_{efs_case}_{efs_speed}.csv" elif profile == "eer": return DATA + f"eer/{eer_demand_file()}" + elif profile == "servm": + return servm_demand_files() elif profile == "ferc": return [ DATA + "pudl/out_ferc714__hourly_estimated_state_demand.parquet", @@ -381,13 +392,23 @@ def demand_raw_data(wildcards): def demand_disaggregate_data(wildcards): - """CLIU county-level industrial loads are the only disaggregation input. + """Extra per-profile input needed to spread zonal demand over buses. - All other end uses disaggregate by population and need no extra file. + CLIU county-level industrial loads serve the industry end use; the SERVM + profile needs its precomputed (region, bus) allocation weights. Everything + else disaggregates by population and needs no extra file. """ - if wildcards.end_use != "industry": - return [] - return DATA + "industry_load/2014_update_20170910-0116.csv" + if wildcards.end_use == "industry": + return DATA + "industry_load/2014_update_20170910-0116.csv" + if ( + wildcards.end_use == "power" + and config["electricity"]["demand"]["profile"] == "servm" + ): + return ( + DEMAND + + f"{wildcards.interconnect}/servm_load_weights_s{wildcards.simpl}.csv" + ) + return [] def demand_scaling_data(wildcards): @@ -410,6 +431,8 @@ def demand_scaling_data(wildcards): return [] elif profile == "eer": return [] + elif profile == "servm": + return [] else: return "" @@ -423,14 +446,20 @@ rule build_electrical_demand: profile_year=pd.to_datetime(config["snapshots"]["start"]).year, planning_horizons=config["scenario"]["planning_horizons"], renewable_weather_years=config["renewable_weather_years"], + servm_weather_years=config["electricity"]["demand"]["scenario"].get( + "servm_weather_years", [] + ), snapshots=config["snapshots"], pudl_path=config_provider("pudl_path"), input: network=NETWORKS + "{interconnect}/elec_s{simpl}.nc", demand_files=demand_raw_data, + dissagregate_files=demand_disaggregate_data, demand_scaling_file=demand_scaling_data, output: elec_demand=DEMAND + "{interconnect}/{end_use}_electricity_s{simpl}.csv", + zonal_components=DEMAND + + "{interconnect}/{end_use}_zonal_components_s{simpl}.parquet", log: LOGS + "{interconnect}/{end_use}_build_demand_s{simpl}.log", benchmark: diff --git a/workflow/rules/retrieve.smk b/workflow/rules/retrieve.smk index 38cda565..290fa245 100644 --- a/workflow/rules/retrieve.smk +++ b/workflow/rules/retrieve.smk @@ -103,6 +103,40 @@ rule retrieve_eer_demand_data: "../scripts/retrieve_eer_data.py" +CPUC_SERVM_URL = "https://files.cpuc.ca.gov/energy/modeling/2026_servm_updates/" + + +rule retrieve_cpuc_servm_load: + wildcard_constraints: + servm_year="2026|2028|2030|2032|2035|2037|2040|2042|2045", + params: + url=lambda wildcards: CPUC_SERVM_URL + + f"HourlyLoad_CA_Regions_V2025E_2224_Mon_{wildcards.servm_year}.csv", + output: + DATA + "cpuc/servm/HourlyLoad_CA_Regions_V2025E_2224_Mon_{servm_year}.csv", + resources: + mem_mb=5000, + log: + "logs/retrieve/retrieve_cpuc_servm_load_{servm_year}.log", + retries: 2 + script: + "../scripts/retrieve_cpuc_data.py" + + +rule retrieve_cpuc_baseline_generators: + params: + url=CPUC_SERVM_URL + "BaselineGeneratorList_CAISO.xlsx", + output: + DATA + "cpuc/BaselineGeneratorList_CAISO.xlsx", + resources: + mem_mb=5000, + log: + "logs/retrieve/retrieve_cpuc_baseline_generators.log", + retries: 2 + script: + "../scripts/retrieve_cpuc_data.py" + + sector_datafiles = [ # heating sector "population/DECENNIALDHC2020.P1-Data.csv", diff --git a/workflow/scripts/build_demand.py b/workflow/scripts/build_demand.py index f8c3f946..cc5bd9ab 100644 --- a/workflow/scripts/build_demand.py +++ b/workflow/scripts/build_demand.py @@ -5,6 +5,7 @@ import calendar import logging +import re import sys from abc import ABC, abstractmethod from pathlib import Path @@ -36,6 +37,7 @@ def __init__(self, read_strategy, write_strategy) -> None: """(read_strategy: ReadStrategy, write_strategy: WriteStrategy).""" self._read_strategy = read_strategy self._write_strategy = write_strategy + self._zonal_demand = None @property def read_strategy(self): # returns ReadStrategy: @@ -47,17 +49,46 @@ def write_strategy(self): # returns WriteStrategy: """The Context maintains a reference to the Strategy objects.""" return self._write_strategy + @property + def zonal_demand(self): + """ + The zonal (pre-disaggregation) demand read during the last prepare_* call. + + This is the reader's own output: hourly demand indexed by + (snapshot, sector, subsector, fuel) with one column per source zone. + It is kept so the rule can persist a component-resolved zonal artifact + without re-reading the (large) source files. ``None`` until a + ``prepare_*`` method has run. + """ + return self._zonal_demand + def _read(self) -> pd.DataFrame: """Delegate reading to the strategy.""" - return self._read_strategy.read_demand() + demand = self._read_strategy.read_demand() + self._zonal_demand = demand + return demand def _write(self, demand: pd.DataFrame, zone: str, **kwargs) -> pd.DataFrame: """Delegate writing to the strategy.""" return self._write_strategy.dissagregate_demand(demand, zone, **kwargs) + def _apply_default_subsector(self, kwargs: dict) -> dict: + """Let a read strategy nominate which subsector is the modeled load. + + Readers that resolve demand into several components (SERVM) keep every + component on the ``subsector`` index level, but only one of them is the + load the model should see. Defaulting it here keeps ``main()`` generic; + an explicit ``subsector=`` argument still wins. + """ + default_subsector = getattr(self._read_strategy, "default_subsector", None) + if default_subsector is not None: + kwargs.setdefault("subsector", default_subsector) + return kwargs + def prepare_demand(self, **kwargs) -> pd.DataFrame: """Read in and dissagregate demand.""" demand = self._read() + kwargs = self._apply_default_subsector(kwargs) return self._write(demand, self._read_strategy.zone, **kwargs) def prepare_multiple_demands( @@ -71,6 +102,7 @@ def prepare_multiple_demands( fuels = [fuels] demand = self._read() + kwargs = self._apply_default_subsector(kwargs) data = {} for fuel in fuels: @@ -96,6 +128,11 @@ class ReadStrategy(ABC): of some algorithm. """ + # Subsector holding the modeled load, for readers that resolve demand into + # several components. ``None`` means the reader emits a single component + # and the caller should not filter on the subsector level. + default_subsector: ClassVar[str | None] = None + def __init__(self, filepath: str | list[str] | None = None) -> None: self.filepath = filepath @@ -531,6 +568,339 @@ def _format_data(self, data: dict[int, pd.DataFrame]) -> pd.DataFrame: return df +class ReadServm(ReadStrategy): + """Reads CPUC SERVM hourly load for the six California load regions. + + The CPUC publishes one CSV per forecast year + (``HourlyLoad_CA_Regions_V2025E_2224_Mon_{year}.csv``). Each file stacks 25 + weather years (2000-2024) of a full year of hourly values for every + (region, component) pair, in fixed Pacific Standard Time with no DST. + + Only ``Net Load`` is the load the model dispatches against, but every + published component (``Load``, ``BTMPV``, ``EV``, ``DATA_CEN``, ...) is + carried through on the ``subsector`` index level so the zonal artifact + stays component-resolved. ``default_subsector`` selects ``Net Load`` at + disaggregation time. + """ + + MODEL_YEARS: ClassVar[tuple[int, ...]] = ( + 2026, + 2028, + 2030, + 2032, + 2035, + 2037, + 2040, + 2042, + 2045, + ) + WEATHER_YEARS: ClassVar[tuple[int, ...]] = tuple(range(2000, 2025)) + REGIONS: ClassVar[tuple[str, ...]] = ("IID", "LADWP", "NCNC", "PGE", "SCE", "SDGE") + + # First seven columns are the (Weather Year, Season, Month, Day, + # Day of Month, Hour, Hour of Day) calendar block. They must be selected + # positionally: "Hour of Day" carries the stray upper-level labels + # ('Region', 'Unit Type'), so testing the upper levels for blankness + # misses it. + INDEX_COLUMNS: ClassVar[int] = 7 + INDEX_NAMES: ClassVar[tuple[str, ...]] = ( + "Weather Year", + "Season", + "Month", + "Day", + "Day of Month", + "Hour", + "Hour of Day", + ) + WEATHER_YEAR_COLUMN: ClassVar[str] = "Weather Year" + LOAD_COMPONENT: ClassVar[str] = "Net Load" + default_subsector: ClassVar[str | None] = "Net Load" + + HOURS_PER_YEAR: ClassVar[int] = const.HOURS_PER_YEAR + # SERVM strips are fixed PST (UTC-8) with no daylight-saving transition. + # Verified empirically against the BTMPV solar-noon centroid, which sits at + # hour 12.5 in December and 12.7 in July - a DST-observing series would move + # by a full hour between the two. + PST_TO_UTC_SHIFT: ClassVar[int] = 8 + + # The forecast year is the trailing "_YYYY" of the basename. Anchoring on + # the extension keeps the vintage tag ("V2025E") out of the match. + FILENAME_YEAR: ClassVar[str] = r"_(\d{4})\.csv$" + + def __init__( + self, + filepath: str | list[str] | None = None, + planning_horizons: list[int] | None = None, + servm_weather_years: list[int] | None = None, + renewable_weather_years: list[int] | None = None, + snapshots: pd.MultiIndex | pd.DatetimeIndex | None = None, + ) -> None: + super().__init__(filepath) + self._zone = "servm" + self.planning_horizons = self._validate_planning_horizons(planning_horizons) + self.weather_year = self._validate_weather_years( + servm_weather_years, + renewable_weather_years, + ) + self.period_snapshots = self._validate_snapshots(snapshots) + self.files = self._index_files_by_year(filepath) + + @property + def zone(self): # noqa: D102 + return self._zone + + @classmethod + def _validate_planning_horizons( + cls, + planning_horizons: list[int] | None, + ) -> list[int]: + if not planning_horizons: + raise ValueError("SERVM demand requires scenario.planning_horizons.") + + years = [int(year) for year in planning_horizons] + invalid_years = sorted(set(years) - set(cls.MODEL_YEARS)) + if invalid_years: + raise ValueError( + f"SERVM demand supports planning_horizons {cls.MODEL_YEARS}; " + f"received unsupported year(s): {invalid_years}.", + ) + return years + + @classmethod + def _validate_weather_years( + cls, + servm_weather_years: list[int] | None, + renewable_weather_years: list[int] | None = None, + ) -> int: + if not servm_weather_years: + raise ValueError( + "SERVM demand requires electricity.demand.scenario.servm_weather_years with exactly one weather year.", + ) + + years = [int(year) for year in servm_weather_years] + if len(years) > 1: + raise NotImplementedError( + "multiple electricity.demand.scenario.servm_weather_years requires " + "stochastic scenarios (phase 3); the electrical demand output path is " + f"not weather-year specific, so only one entry can be built. Received {years}.", + ) + + weather_year = years[0] + if weather_year not in cls.WEATHER_YEARS: + raise ValueError( + f"SERVM demand supports weather years {cls.WEATHER_YEARS}; received {weather_year}.", + ) + + if renewable_weather_years: + renewable = {int(year) for year in renewable_weather_years} + if renewable != {weather_year}: + logger.warning( + "SERVM weather year %s does not match renewable_weather_years %s. " + "Load and renewable profiles will be drawn from different weather " + "years; set them equal unless the mismatch is intentional.", + weather_year, + sorted(renewable), + ) + + return weather_year + + def _validate_snapshots( + self, + snapshots: pd.MultiIndex | pd.DatetimeIndex | None, + ) -> dict[int, pd.DatetimeIndex]: + """Group the network snapshots by investment period. + + SERVM strips carry no usable absolute calendar of their own (see + :meth:`_assign_snapshots`), so the model's own snapshots are the only + source of timestamps and are therefore required. + """ + if snapshots is None or len(snapshots) == 0: + raise ValueError( + "SERVM demand requires the network snapshots to map its hourly strips onto; none were provided.", + ) + + if isinstance(snapshots, pd.MultiIndex): + periods = np.asarray(snapshots.get_level_values(0)) + stamps = pd.DatetimeIndex(snapshots.get_level_values(-1)) + else: + stamps = pd.DatetimeIndex(snapshots) + periods = np.asarray(stamps.year) + + by_period = {} + for period in pd.unique(periods): + by_period[int(period)] = stamps[periods == period] + + missing = sorted(set(self.planning_horizons) - set(by_period)) + if missing: + raise ValueError( + f"Network snapshots contain no timesteps for planning horizon(s) {missing}; " + f"found periods {sorted(by_period)}.", + ) + return by_period + + def _index_files_by_year(self, filepath: str | list[str] | None) -> dict[int, str]: + """Map each input file to the forecast year parsed out of its basename. + + Indexing on the parsed year rather than on list order keeps the reader + correct however snakemake happens to order ``demand_files``. + """ + if not filepath: + raise ValueError("Must provide filepath(s) for SERVM data.") + + files = [filepath] if isinstance(filepath, str) else list(filepath) + + indexed: dict[int, str] = {} + for f in files: + match = re.search(self.FILENAME_YEAR, Path(f).name) + if not match: + raise ValueError( + f"Cannot parse a forecast year out of SERVM filename '{Path(f).name}'; " + "expected it to end in '_YYYY.csv'.", + ) + year = int(match.group(1)) + if year in indexed and indexed[year] != f: + raise ValueError( + f"Two SERVM files claim forecast year {year}: {indexed[year]} and {f}.", + ) + indexed[year] = f + + missing = sorted(set(self.planning_horizons) - set(indexed)) + if missing: + raise ValueError( + f"No SERVM load file provided for planning horizon(s) {missing}; " + f"the supplied files cover {sorted(indexed)}.", + ) + return indexed + + def _read_data(self) -> dict[int, pd.DataFrame]: + """Reads SERVM profiles for each requested model year.""" + logger.info( + f"Building Load Data using CPUC SERVM demand for weather year {self.weather_year}", + ) + return {year: self._read_model_year(year) for year in self.planning_horizons} + + @staticmethod + def _normalize_column(column: tuple) -> tuple[str, str, str]: + """Blank out pandas' placeholder labels for empty header cells.""" + return tuple("" if str(level).startswith("Unnamed:") else str(level).strip() for level in column) + + def _split_header(self, raw: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + """Split the three-row header into the calendar block and the data block.""" + index_block = raw.iloc[:, : self.INDEX_COLUMNS].copy() + index_names = [str(column[-1]).strip() for column in index_block.columns] + if tuple(index_names) != self.INDEX_NAMES: + raise ValueError( + f"SERVM index columns changed: expected {self.INDEX_NAMES}, found {tuple(index_names)}.", + ) + index_block.columns = index_names + + data_block = raw.iloc[:, self.INDEX_COLUMNS :].copy() + data_block.columns = pd.MultiIndex.from_tuples( + [self._normalize_column(column) for column in data_block.columns], + names=["region", "unit_type", "component"], + ) + return index_block, data_block + + def _validate_components(self, columns: pd.MultiIndex, filepath: str) -> None: + """Every modeled region must publish the modeled load component.""" + available = set(zip(columns.get_level_values("region"), columns.get_level_values("component"), strict=True)) + missing = [ + (region, self.LOAD_COMPONENT) for region in self.REGIONS if (region, self.LOAD_COMPONENT) not in available + ] + if missing: + raise ValueError( + f"SERVM file '{filepath}' is missing the '{self.LOAD_COMPONENT}' column for " + f"region(s) {[region for region, _ in missing]}. The published column layout " + "has changed; the reader must be updated before these results can be trusted.", + ) + + def _read_model_year(self, model_year: int) -> pd.DataFrame: + """Read one forecast year and cut out the configured weather year.""" + filepath = self.files[model_year] + raw = pd.read_csv(filepath, header=[0, 1, 2], low_memory=False) + index_block, data_block = self._split_header(raw) + self._validate_components(data_block.columns, filepath) + + weather_years = pd.to_numeric( + index_block[self.WEATHER_YEAR_COLUMN], + errors="coerce", + ) + block = data_block.loc[(weather_years == self.weather_year).to_numpy()] + if len(block) != self.HOURS_PER_YEAR: + raise ValueError( + f"SERVM file '{filepath}' holds {len(block)} rows for weather year " + f"{self.weather_year}; expected {self.HOURS_PER_YEAR}.", + ) + + block = block.astype(float) + rolled = pd.DataFrame( + np.roll(block.to_numpy(), self.PST_TO_UTC_SHIFT, axis=0), + columns=block.columns, + ) + rolled.index = self._assign_snapshots(model_year) + return self._to_long_format(rolled) + + def _assign_snapshots(self, model_year: int) -> pd.DatetimeIndex: + """Map the hourly strip positionally onto this period's snapshots. + + The strip is assigned by position rather than by rebuilding a calendar, + because neither calendar is the model's. ``get_snapshots`` drops + February 29 from leap planning horizons, so a synthesised + ``date_range(f"{year}-01-01", periods=HOURS_PER_YEAR)`` would run one day + short of the network's own December 31 for 2028/2032/2040. Two calendar + misalignments are accepted as a consequence, and are immaterial for an + hourly capacity-expansion model: + + 1. SERVM lays its hours on a synthetic Monday-start calendar, so + weekday-versus-weekend hours do not line up with the real weekdays of + the planning horizon. + 2. For a leap *weather* year the strip contains February 29 and omits + December 31, while the model snapshots do the opposite. Every hour + after February therefore lands one calendar day earlier than it sat + in the source file. + """ + snapshots = self.period_snapshots[model_year] + if len(snapshots) != self.HOURS_PER_YEAR: + raise ValueError( + f"Planning horizon {model_year} has {len(snapshots)} snapshots; SERVM " + f"demand needs a full {self.HOURS_PER_YEAR}-hour year to map onto.", + ) + return pd.DatetimeIndex(snapshots, name="snapshot") + + def _to_long_format(self, df: pd.DataFrame) -> pd.DataFrame: + """Move the component level onto the index, keeping regions as columns.""" + components = list(dict.fromkeys(df.columns.get_level_values("component"))) + + frames = {} + for component in components: + subset = df.loc[:, df.columns.get_level_values("component") == component] + subset.columns = pd.Index(subset.columns.get_level_values("region"), name=None) + if subset.columns.has_duplicates: + duplicates = sorted(subset.columns[subset.columns.duplicated()].unique()) + raise ValueError( + f"SERVM component '{component}' appears more than once for region(s) {duplicates}.", + ) + frames[component] = subset + + # Components missing for a region (EV and friends exist only for + # PGE/SCE/SDGE) align to NaN, which is the honest reading of "not + # published" and never reaches the model: only LOAD_COMPONENT, present + # everywhere, is disaggregated. + long = pd.concat(frames, names=["subsector"]) + return long.reorder_levels(["snapshot", "subsector"]).sort_index() + + def _format_data(self, data: dict[int, pd.DataFrame]) -> pd.DataFrame: + """Formats raw SERVM data to the demand strategy contract.""" + df = pd.concat(data.values()).reset_index() + # snapshots are taken from the network, so they are datetimes already; + # _format_snapshot_index() cannot be used here because MultiIndex + # set_levels() requires unique level values. + df["snapshot"] = pd.to_datetime(df["snapshot"]) + df["sector"] = "all" + df["fuel"] = "electricity" + return df.set_index(["snapshot", "sector", "subsector", "fuel"]).sort_index() + + class ReadEulp(ReadStrategy): """Reads in End Use Load Profile data.""" @@ -1553,7 +1923,7 @@ def dissagregate_demand( df: pd.DataFrame Demand dataframe zone: str - Zones of demand ('ba', 'state', 'reeds') + Zones of demand ('ba', 'state', 'reeds', 'servm') sector: Optional[str | List[str]] = None, Sectors to group subsector: Optional[str | List[str]] = None, @@ -1575,7 +1945,7 @@ def dissagregate_demand( """ # 'state' is states based on power regions # 'full_state' is actual geographic boundaries - assert zone in ("ba", "state", "reeds") + assert zone in ("ba", "state", "reeds", "servm") self._check_datastructure(df) # get zone area demand for specific sector and fuel @@ -1765,6 +2135,94 @@ def _get_load_allocation_factor( return bus_load.load_weight / zone_loads +class WriteServm(WritePopulation): + """ + Disaggregates SERVM regional demand with the precomputed allocation weights. + + ``build_servm_load_weights`` writes a long table of (bus, servm_region, laf) + shares that already sum to 1.0 within every region. Because a cluster bus can + straddle two SERVM regions, a bus may appear under more than one region, so + the allocation cannot be expressed as the one-zone-per-bus mapping the base + class builds. Pivoting the table to a (region x bus) matrix and taking the + matrix product against the (snapshot x region) demand handles straddling + buses exactly: each bus receives the sum of its share of every region it + overlaps. + """ + + def __init__(self, n: pypsa.Network, filepath: str) -> None: + super().__init__(n) + self.filepath = filepath + self.weights = self._read_weights(filepath) + + def _read_weights(self, filepath: str) -> pd.DataFrame: + """Read the weights table and pivot it to a (region x bus) matrix.""" + if isinstance(filepath, list | tuple): + if len(filepath) != 1: + raise ValueError( + f"SERVM disaggregation needs exactly one weights file; received {list(filepath)}.", + ) + filepath = filepath[0] + + df = pd.read_csv(filepath) + missing_columns = {"bus", "servm_region", "laf"}.difference(df.columns) + if missing_columns: + raise ValueError( + f"SERVM weights file '{filepath}' is missing column(s) {sorted(missing_columns)}.", + ) + + df["bus"] = df.bus.astype(str) + unknown = sorted(set(df.bus) - set(self.n.buses.index.astype(str))) + if unknown: + raise ValueError( + f"SERVM weights file '{filepath}' allocates demand to {len(unknown)} bus(es) that " + f"are not in the network: {unknown[:10]}. The weights were built against a " + "different network than the one demand is being attached to.", + ) + + weights = df.pivot_table( + index="servm_region", + columns="bus", + values="laf", + aggfunc="sum", + fill_value=0.0, + ) + logger.info( + "Allocating SERVM demand over %d buses from %d regions.", + weights.shape[1], + weights.shape[0], + ) + return weights.astype(float) + + def dissagregate_demand( + self, + df: pd.DataFrame, + zone: str, + sector: str | list[str] | None = None, + subsector: str | list[str] | None = None, + fuel: str | list[str] | None = None, + sns: pd.DatetimeIndex | None = None, + ) -> pd.DataFrame: + """Allocate regional demand to buses via the weights matrix product.""" + assert zone == "servm" + self._check_datastructure(df) + + demand = self._filter_demand(df, sector, subsector, fuel, sns) + demand = self._group_demand(demand) + if demand.empty: + demand = self._make_empty_demand(columns=df.columns) + demand = demand.astype(float) + + weights = self.weights.reindex(index=demand.columns, fill_value=0.0) + unweighted = weights.index[~weights.index.isin(self.weights.index)] + if len(unweighted): + logger.warning( + "No bus weights found for SERVM region(s) %s; their demand is dropped.", + sorted(unweighted), + ) + + return demand.dot(weights) + + class WriteIndustrial(WriteStrategy): """ Based on county level energy use from 2014. @@ -2310,6 +2768,11 @@ def get_demand_params( scaling_method = "aeo_electricity" elif demand_profile == "eer": scaling_method = None + elif demand_profile == "servm": + # SERVM publishes one file per forecast year, so no scaling is + # needed; its regions are their own disaggregation zone. + demand_disaggregation = "servm" + scaling_method = None else: logger.warning( f"No scaling method available for {demand_profile} profile. Setting to 'aeo_electricity'", @@ -2495,6 +2958,16 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]: ) sns = n.snapshots.get_level_values(1) + elif demand_profile == "servm": + reader = ReadServm( + demand_files, + planning_horizons=planning_horizons, + servm_weather_years=snakemake.params.get("servm_weather_years", None), + renewable_weather_years=snakemake.params.get("renewable_weather_years", None), + snapshots=n.snapshots, + ) + sns = n.snapshots.get_level_values(1) + elif demand_profile == "ferc": assert profile_year in range(2018, 2024) @@ -2555,6 +3028,8 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]: elif demand_disaggregation == "cliu": cliu_file = snakemake.input.dissagregate_files writer = WriteIndustrial(n, cliu_file) + elif demand_disaggregation == "servm": + writer = WriteServm(n, snakemake.input.dissagregate_files) else: raise NotImplementedError @@ -2574,6 +3049,16 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]: sns=sns, ) # dict[str, pd.DataFrame] + # persist the reader's own zonal, component-resolved demand before it is + # disaggregated onto buses. For single-component profiles (efs, eer, ...) + # this is simply the subsector='all' slice, so no profile needs special + # casing here. + zonal_output = snakemake.output.get("zonal_components", None) + if zonal_output: + zonal_demand = demand_converter.zonal_demand + assert zonal_demand is not None, "no zonal demand captured during read" + zonal_demand.astype(float).round(4).to_parquet(zonal_output) + # scale demand and align snapshots. this is outside the main read/write # strategy as extra arguments are required to fill in data if scaling_method == "aeo_electricity": diff --git a/workflow/scripts/retrieve_cpuc_data.py b/workflow/scripts/retrieve_cpuc_data.py new file mode 100644 index 00000000..ea9d1627 --- /dev/null +++ b/workflow/scripts/retrieve_cpuc_data.py @@ -0,0 +1,24 @@ +"""Retrieve CPUC modeling datasets (SERVM hourly load, baseline generator list).""" + +import logging +from pathlib import Path + +from _helpers import configure_logging, progress_retrieve + +logger = logging.getLogger(__name__) + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake("retrieve_cpuc_servm_load", servm_year="2026") + + configure_logging(snakemake) + + output = Path(snakemake.output[0]) + output.parent.mkdir(parents=True, exist_ok=True) + + url = snakemake.params.url + logger.info(f"Downloading CPUC data from '{url}'.") + progress_retrieve(url, output) diff --git a/workflow/scripts/test/test_build_demand_servm.py b/workflow/scripts/test/test_build_demand_servm.py new file mode 100644 index 00000000..d6640eab --- /dev/null +++ b/workflow/scripts/test/test_build_demand_servm.py @@ -0,0 +1,418 @@ +"""Tests for the CPUC SERVM demand profile reader.""" + +import logging + +import numpy as np +import pandas as pd +import pytest +from _helpers import get_multiindex_snapshots +from build_demand import Context, ReadServm + +HOURS = ReadServm.HOURS_PER_YEAR + +# Component layout copied from the published file. IID/LADWP/NCNC carry the six +# universal components; PGE/SCE/SDGE additionally carry the EV, BTM-storage, +# climate-change and data-centre components. +SIMPLE_COMPONENTS = ("Load", "Modified Load", "Net Load", "BTMPV", "AAEE", "AAFS") +SIMPLE_UNITS = ("", "", "", "R", "R", "R") +FULL_COMPONENTS = ( + "Load", + "Modified Load", + "Net Load", + "BTMPV", + "AAEE", + "EV", + "BTMStorageShapeDischarge", + "AAFS", + "BTMStorageShapeCharge", + "CLIM_CHG_nou", + "CLIM_CHG_gen", + "DATA_CEN", +) +FULL_UNITS = ("", "", "", "R", "R", "R", "R", "R", "R", "R", "R", "R") + +REGION_LAYOUT = ( + ("IID", SIMPLE_COMPONENTS, SIMPLE_UNITS), + ("LADWP", SIMPLE_COMPONENTS, SIMPLE_UNITS), + ("NCNC", SIMPLE_COMPONENTS, SIMPLE_UNITS), + ("PGE", FULL_COMPONENTS, FULL_UNITS), + ("SCE", FULL_COMPONENTS, FULL_UNITS), + ("SDGE", FULL_COMPONENTS, FULL_UNITS), +) + +# Literal first two header rows of the real file, for the leading calendar block. +# The seventh column carries the stray ('Region', 'Unit Type') labels. +INDEX_HEADER_0 = ["", "", "", "", "", "", "Region"] +INDEX_HEADER_1 = ["", "", "", "", "", "", "Unit Type"] +INDEX_HEADER_2 = [ + "Weather Year", + "Season", + "Month", + "Day", + "Day of Month", + "Hour", + "Hour of Day", +] + + +def _column_layout(drop=()): + """(region, unit, component) triples in published order, minus `drop`.""" + columns = [] + for region, components, units in REGION_LAYOUT: + for component, unit in zip(components, units, strict=True): + if (region, component) in drop: + continue + columns.append((region, unit, component)) + return columns + + +def _value(column_index: int, weather_year: int, hour: int, offset: int = 0) -> int: + """Deterministic, collision-free cell value.""" + return offset + column_index * 100_000 + (weather_year - 2000) * 10_000 + hour + + +def write_servm_fixture( + path, + weather_years=(2000, 2001), + offset: int = 0, + drop=(), +): + """Write a synthetic SERVM CSV carrying the real three-row header quirks.""" + columns = _column_layout(drop=drop) + n_columns = len(columns) + + header_lines = [ + ",".join(INDEX_HEADER_0 + [region for region, _, _ in columns]), + ",".join(INDEX_HEADER_1 + [unit for _, unit, _ in columns]), + ",".join(INDEX_HEADER_2 + [component for _, _, component in columns]), + ] + + hours = np.arange(HOURS) + blocks = [] + for weather_year in weather_years: + index_block = pd.DataFrame( + { + "Weather Year": weather_year, + "Season": "Winter", + "Month": 1, + "Day": 1, + "Day of Month": 1, + "Hour": hours + 1, + "Hour of Day": (hours % 24) + 1, + }, + ) + data = np.empty((HOURS, n_columns), dtype=np.int64) + for column_index in range(n_columns): + data[:, column_index] = _value(column_index, weather_year, hours, offset) + blocks.append(pd.concat([index_block, pd.DataFrame(data)], axis=1)) + + body = pd.concat(blocks, ignore_index=True) + with open(path, "w") as f: + f.write("\n".join(header_lines) + "\n") + body.to_csv(f, header=False, index=False, lineterminator="\n") + return path + + +def make_snapshots(planning_horizons, base_year=2019): + """Snapshots exactly as the pipeline builds them.""" + return get_multiindex_snapshots( + { + "start": f"{base_year}-01-01 00:00", + "end": f"{base_year}-12-31 23:00", + "inclusive": "both", + }, + planning_horizons, + ) + + +def column_index_of(region: str, component: str, drop=()) -> int: + columns = _column_layout(drop=drop) + return columns.index( + next(c for c in columns if c[0] == region and c[2] == component), + ) + + +@pytest.fixture(scope="module") +def servm_file(tmp_path_factory): + path = tmp_path_factory.mktemp("servm") / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + return str(write_servm_fixture(path)) + + +@pytest.fixture(scope="module") +def servm_demand(servm_file): + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=make_snapshots([2028]), + ) + return reader.read_demand() + + +def test_multiheader_parse_recovers_region_and_component(servm_demand): + """The three-row header resolves to regions on columns, components on subsector.""" + assert list(servm_demand.columns) == list(ReadServm.REGIONS) + assert servm_demand.index.names == ["snapshot", "sector", "subsector", "fuel"] + + subsectors = set(servm_demand.index.get_level_values("subsector")) + assert "Net Load" in subsectors + assert set(servm_demand.index.get_level_values("sector")) == {"all"} + assert set(servm_demand.index.get_level_values("fuel")) == {"electricity"} + assert len(servm_demand) == HOURS * len(set(FULL_COMPONENTS) | set(SIMPLE_COMPONENTS)) + + +def test_weather_year_filter_selects_correct_block(servm_file): + """The chosen weather year selects its own 8760-row block, not the first one.""" + snapshots = make_snapshots([2028]) + column = column_index_of("SCE", "Net Load") + + values = {} + for weather_year in (2000, 2001): + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[weather_year], + snapshots=snapshots, + ) + demand = reader.read_demand() + # hour 0 of the strip lands PST_TO_UTC_SHIFT hours into the year + stamp = snapshots.get_level_values(1)[ReadServm.PST_TO_UTC_SHIFT] + values[weather_year] = demand.loc[(stamp, "all", "Net Load", "electricity"), "SCE"] + + assert values[2000] == _value(column, 2000, 0) + assert values[2001] == _value(column, 2001, 0) + assert values[2000] != values[2001] + + +def test_all_components_preserved_on_subsector_level(servm_demand): + """Every published component survives; EV stays absent where CPUC omits it.""" + subsectors = set(servm_demand.index.get_level_values("subsector")) + assert subsectors == set(SIMPLE_COMPONENTS) | set(FULL_COMPONENTS) + + ev = servm_demand.xs("EV", level="subsector") + assert ev[["PGE", "SCE", "SDGE"]].notna().all().all() + assert ev[["IID", "LADWP", "NCNC"]].isna().all().all() + + net_load = servm_demand.xs("Net Load", level="subsector") + assert net_load.notna().all().all() + + +def test_net_load_is_the_default_subsector(servm_file): + """Context filters to Net Load without main() knowing about SERVM.""" + assert ReadServm.default_subsector == "Net Load" + + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=make_snapshots([2028]), + ) + + class _RecordingWriter: + def __init__(self): + self.kwargs = None + + def dissagregate_demand(self, df, zone, **kwargs): + self.kwargs = dict(kwargs, zone=zone) + return df + + writer = _RecordingWriter() + Context(reader, writer).prepare_demand() + + assert writer.kwargs["subsector"] == "Net Load" + assert writer.kwargs["zone"] == "servm" + + +def test_explicit_subsector_overrides_the_default(servm_file): + """An explicit subsector= argument still wins over the reader's default.""" + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=make_snapshots([2028]), + ) + + class _RecordingWriter: + def __init__(self): + self.kwargs = None + + def dissagregate_demand(self, df, zone, **kwargs): + self.kwargs = kwargs + return df + + writer = _RecordingWriter() + Context(reader, writer).prepare_demand(subsector="BTMPV") + + assert writer.kwargs["subsector"] == "BTMPV" + + +def test_snapshots_align_with_leap_model_year(servm_file): + """A leap planning horizon keeps the network's own (Feb-29-free) snapshots.""" + snapshots = make_snapshots([2028]) + timesteps = snapshots.get_level_values(1) + assert len(timesteps) == HOURS # the pipeline drops Feb 29 + + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=snapshots, + ) + demand = reader.read_demand() + + stamps = pd.DatetimeIndex(demand.index.get_level_values("snapshot")) + assert stamps.year.unique().tolist() == [2028] + assert not ((stamps.month == 2) & (stamps.day == 29)).any() + assert stamps.max() == pd.Timestamp("2028-12-31 23:00") + + net_load = demand.xs("Net Load", level="subsector") + got = pd.DatetimeIndex(net_load.index.get_level_values("snapshot")) + pd.testing.assert_index_equal(got, timesteps, check_names=False) + + +def test_pst_shift_rolls_by_eight(servm_file): + """PST (UTC-8) with no DST: hour 0 of the strip is 08:00 UTC.""" + snapshots = make_snapshots([2028]) + timesteps = snapshots.get_level_values(1) + column = column_index_of("IID", "Net Load") + + reader = ReadServm( + servm_file, + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=snapshots, + ) + demand = reader.read_demand().xs("Net Load", level="subsector") + + shift = ReadServm.PST_TO_UTC_SHIFT + assert demand.loc[(timesteps[shift], "all", "electricity"), "IID"] == _value(column, 2000, 0) + # the tail of the strip wraps onto the first hours of the year + assert demand.loc[(timesteps[0], "all", "electricity"), "IID"] == _value(column, 2000, HOURS - shift) + + +def test_rejects_unsupported_planning_horizon(tmp_path): + with pytest.raises(ValueError, match="unsupported year"): + ReadServm( + str(tmp_path / "HourlyLoad_2027.csv"), + planning_horizons=[2027], + servm_weather_years=[2019], + snapshots=make_snapshots([2027]), + ) + + +def test_rejects_unsupported_weather_year(tmp_path): + with pytest.raises(ValueError, match="supports weather years"): + ReadServm( + str(tmp_path / "HourlyLoad_2028.csv"), + planning_horizons=[2028], + servm_weather_years=[1999], + snapshots=make_snapshots([2028]), + ) + + +def test_multiple_weather_years_raises_not_implemented(tmp_path): + with pytest.raises(NotImplementedError, match="stochastic scenarios"): + ReadServm( + str(tmp_path / "HourlyLoad_2028.csv"), + planning_horizons=[2028], + servm_weather_years=[2018, 2019], + snapshots=make_snapshots([2028]), + ) + + +def test_mismatched_renewable_weather_years_warns(tmp_path, caplog): + path = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + path.touch() + + with caplog.at_level(logging.WARNING, logger="build_demand"): + ReadServm( + str(path), + planning_horizons=[2028], + servm_weather_years=[2019], + renewable_weather_years=[2012], + snapshots=make_snapshots([2028]), + ) + + assert "does not match renewable_weather_years" in caplog.text + + +def test_matched_renewable_weather_years_do_not_warn(tmp_path, caplog): + path = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + path.touch() + + with caplog.at_level(logging.WARNING, logger="build_demand"): + ReadServm( + str(path), + planning_horizons=[2028], + servm_weather_years=[2019], + renewable_weather_years=[2019], + snapshots=make_snapshots([2028]), + ) + + assert "does not match renewable_weather_years" not in caplog.text + + +def test_missing_net_load_column_raises(tmp_path): + """A CPUC layout change that drops a region's Net Load must fail loudly.""" + path = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + write_servm_fixture(path, weather_years=(2000,), drop=(("SDGE", "Net Load"),)) + + reader = ReadServm( + str(path), + planning_horizons=[2028], + servm_weather_years=[2000], + snapshots=make_snapshots([2028]), + ) + with pytest.raises(ValueError, match="missing the 'Net Load' column"): + reader.read_demand() + + +def test_files_indexed_by_basename_year_not_order(tmp_path): + """Files are matched to horizons by their filename year, whatever the order.""" + file_2026 = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2026.csv" + file_2028 = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + write_servm_fixture(file_2026, weather_years=(2000,), offset=0) + write_servm_fixture(file_2028, weather_years=(2000,), offset=1_000_000_000) + + # deliberately reversed relative to the planning horizons + reader = ReadServm( + [str(file_2028), str(file_2026)], + planning_horizons=[2026, 2028], + servm_weather_years=[2000], + snapshots=make_snapshots([2026, 2028]), + ) + assert reader.files == {2026: str(file_2026), 2028: str(file_2028)} + + demand = reader.read_demand().xs("Net Load", level="subsector") + column = column_index_of("PGE", "Net Load") + shift = ReadServm.PST_TO_UTC_SHIFT + + hour_2026 = pd.Timestamp("2026-01-01 00:00") + pd.Timedelta(hours=shift) + hour_2028 = pd.Timestamp("2028-01-01 00:00") + pd.Timedelta(hours=shift) + assert demand.loc[(hour_2026, "all", "electricity"), "PGE"] == _value(column, 2000, 0) + assert demand.loc[(hour_2028, "all", "electricity"), "PGE"] == _value(column, 2000, 0, offset=1_000_000_000) + + +def test_missing_file_for_planning_horizon_raises(tmp_path): + path = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2026.csv" + path.touch() + + with pytest.raises(ValueError, match=r"No SERVM load file provided for planning horizon"): + ReadServm( + str(path), + planning_horizons=[2026, 2028], + servm_weather_years=[2019], + snapshots=make_snapshots([2026, 2028]), + ) + + +def test_snapshots_are_required(tmp_path): + path = tmp_path / "HourlyLoad_CA_Regions_V2025E_2224_Mon_2028.csv" + path.touch() + + with pytest.raises(ValueError, match="requires the network snapshots"): + ReadServm( + str(path), + planning_horizons=[2028], + servm_weather_years=[2019], + ) diff --git a/workflow/scripts/test/test_build_demand_servm_write.py b/workflow/scripts/test/test_build_demand_servm_write.py new file mode 100644 index 00000000..2c26b596 --- /dev/null +++ b/workflow/scripts/test/test_build_demand_servm_write.py @@ -0,0 +1,139 @@ +"""Tests for the SERVM demand write (disaggregation) strategy.""" + +import numpy as np +import pandas as pd +import pypsa +import pytest +from _helpers import get_multiindex_snapshots +from build_demand import WriteServm + +REGIONS = ("PGE", "SCE") + + +def make_network(buses=("bus_a", "bus_b", "bus_c")): + n = pypsa.Network() + n.snapshots = get_multiindex_snapshots( + {"start": "2028-01-01 00:00", "end": "2028-01-01 03:00", "inclusive": "both"}, + [2028], + ) + n.set_investment_periods(periods=[2028]) + for bus in buses: + n.add("Bus", bus) + return n + + +def write_weights(path, rows): + """rows: iterable of (bus, servm_region, laf).""" + pd.DataFrame(rows, columns=["bus", "servm_region", "laf"]).to_csv(path, index=False) + return str(path) + + +def make_demand(n, values): + """Zonal demand frame on the reader's 4-level contract. + + ``values`` maps region -> list of hourly values for the Net Load component. + A second (BTMPV) component is always added so the subsector filter is + actually exercised. + """ + snapshots = n.snapshots.get_level_values(1) + frames = [] + for subsector, scale in (("Net Load", 1.0), ("BTMPV", 100.0)): + frame = pd.DataFrame( + {region: np.asarray(series, dtype=float) * scale for region, series in values.items()}, + index=snapshots, + ) + frame.index.name = "snapshot" + frame["sector"] = "all" + frame["subsector"] = subsector + frame["fuel"] = "electricity" + frames.append(frame.set_index(["sector", "subsector", "fuel"], append=True)) + return pd.concat(frames).sort_index() + + +def test_writeservm_matrix_product_matches_manual(tmp_path): + """Bus load is the region-share weighted sum of regional demand.""" + n = make_network() + weights_file = write_weights( + tmp_path / "weights.csv", + [ + ("bus_a", "PGE", 0.75), + ("bus_b", "PGE", 0.25), + ("bus_c", "SCE", 1.0), + ], + ) + + demand = make_demand(n, {"PGE": [100, 200, 300, 400], "SCE": [10, 20, 30, 40]}) + + writer = WriteServm(n, weights_file) + result = writer.dissagregate_demand(demand, "servm", subsector="Net Load") + + assert list(result.columns) == ["bus_a", "bus_b", "bus_c"] + np.testing.assert_allclose(result["bus_a"], [75, 150, 225, 300]) + np.testing.assert_allclose(result["bus_b"], [25, 50, 75, 100]) + np.testing.assert_allclose(result["bus_c"], [10, 20, 30, 40]) + + # the BTMPV component (100x) must not have leaked into the modeled load + assert result.to_numpy().sum() == pytest.approx( + demand.xs("Net Load", level="subsector").to_numpy().sum(), + ) + + +def test_straddling_bus_receives_sum_of_both_regions(tmp_path): + """A cluster spanning two SERVM regions collects a share of each.""" + n = make_network(buses=("bus_a", "bus_b")) + weights_file = write_weights( + tmp_path / "weights.csv", + [ + ("bus_a", "PGE", 0.6), + ("bus_b", "PGE", 0.4), + ("bus_a", "SCE", 0.1), # bus_a straddles PGE and SCE + ("bus_b", "SCE", 0.9), + ], + ) + + demand = make_demand(n, {"PGE": [100, 100, 100, 100], "SCE": [50, 50, 50, 50]}) + + writer = WriteServm(n, weights_file) + result = writer.dissagregate_demand(demand, "servm", subsector="Net Load") + + np.testing.assert_allclose(result["bus_a"], [65.0] * 4) # 0.6*100 + 0.1*50 + np.testing.assert_allclose(result["bus_b"], [85.0] * 4) # 0.4*100 + 0.9*50 + # nothing is created or lost + np.testing.assert_allclose(result.sum(axis=1), [150.0] * 4) + + +def test_weights_bus_not_in_network_raises(tmp_path): + """Weights built against a different network must fail loudly.""" + n = make_network(buses=("bus_a",)) + weights_file = write_weights( + tmp_path / "weights.csv", + [("bus_a", "PGE", 0.5), ("bus_missing", "PGE", 0.5)], + ) + + with pytest.raises(ValueError, match="not in the network"): + WriteServm(n, weights_file) + + +def test_region_without_weights_is_dropped_with_warning(tmp_path, caplog): + """Demand for a region absent from the weights table cannot be allocated.""" + n = make_network(buses=("bus_a",)) + weights_file = write_weights(tmp_path / "weights.csv", [("bus_a", "PGE", 1.0)]) + + demand = make_demand(n, {"PGE": [100] * 4, "SCE": [50] * 4}) + + writer = WriteServm(n, weights_file) + with caplog.at_level("WARNING", logger="build_demand"): + result = writer.dissagregate_demand(demand, "servm", subsector="Net Load") + + np.testing.assert_allclose(result["bus_a"], [100.0] * 4) + assert "No bus weights found for SERVM region(s)" in caplog.text + + +def test_wrong_zone_is_rejected(tmp_path): + n = make_network(buses=("bus_a",)) + weights_file = write_weights(tmp_path / "weights.csv", [("bus_a", "PGE", 1.0)]) + demand = make_demand(n, {"PGE": [100] * 4}) + + writer = WriteServm(n, weights_file) + with pytest.raises(AssertionError): + writer.dissagregate_demand(demand, "state", subsector="Net Load") From 099977c25a3fe3be158313a8da0ec69dcab67a90 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:44:22 -0700 Subject: [PATCH 04/23] Add maintained California SERVM configuration, docs, and phase-2/3 stubs (PR 4) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Ties the SERVM demand profile (PR 2) and the aggregate interface transmission limits (PR 3) together into a runnable, maintained California entry point. New `workflow/repo_data/config/config.california.yaml`: CA-only Western model on CPUC SERVM 2026 demand, at REeDS-zone resolution (clusters must be 4 — the ReEDS backbone cannot cluster below the four CA zones p8/p9/p10/p11), with planning horizons restricted to SERVM forecast years (2030/2035/2040/2045), the RESOLVE CAISO interface caps switched on, and imports/exports enabled at a 25%-of-demand annual volume limit. A clearly-marked commented alternative carries the county-resolution variant (clusters: 58, simpl: 'county'); the county NARIS flowgate file is selected automatically by `add_extra_components` from `topological_boundaries`. Phase-2 hook `conventional.ambient_derate` is stubbed in both the California and default configs and guarded in `add_electricity`: enabling it raises NotImplementedError rather than silently doing nothing. The comments record that when implemented it REPLACES the EIA-860 seasonal derate and must never stack with it or with a UCAP-derated capacity credit. `retrieve.smk` reserves `retrieve_cpuc_thermal_derate` for the profiles it will need. Tests: `config.test.california.yaml` is a SERVM variant of the CA test harness config. It keeps a full 8760-hour snapshot year — `ReadServm._assign_snapshots` maps its hourly strips positionally onto the network's own per-period snapshots and raises unless the horizon carries exactly 8760, so a truncated window cannot build SERVM demand. Three dry-run cases are added (california -> cluster_network, california -> solve_network, and a county-mode case driven by `--config` overrides), and `tests/integration/test_servm_demand_artifacts.py` asserts weight conservation, demand/network alignment, and the PST->UTC roll against a `--until add_demand` build. Docs: a SERVM section in data-demand.md (source URL pattern, the nine forecast years, the six regions and their balancing-area mapping including the CISO-VEA exclusion, Net-Load-only dispatch with the component-resolved zonal artifact, weather-year semantics, fixed PST, and both calendar caveats); the interface constraint documented in model-constraints.md and data-transmission.md including the p8 gap; the `transmission_interface_limits` config row corrected from "not currently consumed"; a California example in config-spatial.md; and a release-notes entry. Co-Authored-By: Claude Fable 5 --- docs/source/config-configuration.md | 9 +- docs/source/config-spatial.md | 43 ++ docs/source/configtables/electricity.csv | 4 +- docs/source/data-demand.md | 100 +++- docs/source/data-transmission.md | 31 ++ docs/source/datatables/demand.csv | 1 + docs/source/model-constraints.md | 36 ++ docs/source/release-notes.md | 30 ++ tests/integration/conftest.py | 106 ++++ .../test_servm_demand_artifacts.py | 146 ++++++ tests/static/test_dag_dryrun.py | 52 +- .../repo_data/config/config.california.yaml | 487 ++++++++++++++++++ workflow/repo_data/config/config.default.yaml | 14 + .../config/config.test.california.yaml | 135 +++++ workflow/rules/retrieve.smk | 6 + workflow/scripts/add_electricity.py | 17 + 16 files changed, 1195 insertions(+), 22 deletions(-) create mode 100644 tests/integration/test_servm_demand_artifacts.py create mode 100644 workflow/repo_data/config/config.california.yaml create mode 100644 workflow/repo_data/config/config.test.california.yaml diff --git a/docs/source/config-configuration.md b/docs/source/config-configuration.md index f6f6b38d..411912f0 100644 --- a/docs/source/config-configuration.md +++ b/docs/source/config-configuration.md @@ -50,8 +50,13 @@ network is aggregated to. `transmission_network` chooses between the ReEDS zonal TAMU synthetic nodal network; `topological_boundaries` sets the zone type used after clustering (county, REeDS zone, state, or balancing area). Use `include` to subset the modeled footprint to specific zones, states, or balancing authorities (mixed zone types are not supported), and -`aggregate` to pre-aggregate buses into larger regions. `interface_transmission_limits` applies -NARIS2024 inter-regional transfer capacity limits and requires the ReEDS backbone. +`aggregate` to pre-aggregate buses into larger regions. `interface_transmission_limits` switches +on the aggregate inter-regional transfer limits read from +`electricity: transmission_interface_limits` and requires the ReEDS backbone; the limits are +applied in `solve_network` as a per-snapshot cap on the total flow across each interface, and +constrain only the import/export links, so they are inert unless `electricity: imports` or +`electricity: exports` is enabled. See {ref}`spatial` for a worked California example +(`workflow/repo_data/config/config.california.yaml`) at both REeDS-zone and county resolution. ```{eval-rst} .. literalinclude:: ../../workflow/repo_data/config/config.default.yaml diff --git a/docs/source/config-spatial.md b/docs/source/config-spatial.md index 04fe7f71..f94deb9b 100644 --- a/docs/source/config-spatial.md +++ b/docs/source/config-spatial.md @@ -35,6 +35,49 @@ model_topology: Alternatively, you can use the code reeds_state: 'CA' option to achieve the same result by specifying the entire state. +A complete, maintained California configuration ships with the repository as +`workflow/repo_data/config/config.california.yaml`. It pairs the footprint below with CPUC +SERVM demand, the RESOLVE CAISO interface limits, and enabled imports/exports: + +```yaml +scenario: + interconnect: [western] + planning_horizons: [2030, 2035, 2040, 2045] # CPUC SERVM forecast years + clusters: [4] # p8, p9, p10, p11 + simpl: [75] + +model_topology: + transmission_network: 'reeds' + topological_boundaries: 'reeds_zone' + interface_transmission_limits: true # RESOLVE CAISO interface caps + include: + reeds_state: ['CA'] +``` + +`clusters` is pinned to 4 because the ReEDS zonal backbone cannot be clustered below the +number of zones in the footprint, and California holds exactly four. + +To run the same footprint at **county resolution**, switch `topological_boundaries` to +`county` and raise `clusters` to 58 — the number of California counties, and the number of +`p06xxx` nodes in the county NARIS interface table. `simpl: ['county']` selects the +county-FIPS fast path in `cluster_simpl`, so the resource layer is built directly on county +boundaries (a numeric `simpl` of at least 58 also works if you want more resource zones than +transmission nodes): + +```yaml +scenario: + clusters: [58] # 58 California counties + simpl: ['county'] # county-FIPS fast path + +model_topology: + topological_boundaries: 'county' +``` + +`add_extra_components` switches from the balancing-area NARIS flowgate file to the county one +(`transmission_capacity_init_AC_county_NARIS2024.csv`) automatically when +`topological_boundaries` is `county` — no other key needs changing. `config.california.yaml` +carries this same block as a commented alternative. + In addition to filtering by `reeds_zone` and `reeds_state`, you can filter by `reeds_ba`, `trans_reg`, and `nerc_reg` shown graphically below. diff --git a/docs/source/configtables/electricity.csv b/docs/source/configtables/electricity.csv index 50eb3437..9cbee3fd 100644 --- a/docs/source/configtables/electricity.csv +++ b/docs/source/configtables/electricity.csv @@ -24,7 +24,7 @@ erm:,,,Energy Reserve Margin settings (used when ERM opt is enabled). Ensures su regional_Co2_limits,--,path,"CSV of per-region CO2 caps in tCO2/yr (``config/policy_constraints/regional_Co2_limits.csv``). Enforced when the ``REM`` keyword is present in the ``{opts}`` wildcard." technology_capacity_targets,--,path,"CSV of forced minimum/maximum capacity builds by technology and region (``config/policy_constraints/technology_capacity_targets.csv``). Enforced when the ``TCT`` keyword is present in the ``{opts}`` wildcard." portfolio_standards,--,path,"CSV of RPS/CES clean-energy fractions by region (``config/policy_constraints/portfolio_standards.csv``). Enforced when the ``RPS`` keyword is present in the ``{opts}`` wildcard (covers both RPS and CES targets)." -transmission_interface_limits,--,path,"CSV of MW limits on flows across inter-regional transmission interfaces (``config/policy_constraints/transmission_interface_limits.csv``), paired with ``model_topology: interface_transmission_limits``. Reserved setting — not currently consumed by the workflow." +transmission_interface_limits,--,path,"CSV of MW limits on flows across inter-regional transmission interfaces (``config/policy_constraints/transmission_interface_limits.csv``, columns ``interface, region_1, region_2, flow_12, flow_21``). Applied by ``solve_network`` when ``model_topology: interface_transmission_limits`` is ``true``, as a per-snapshot cap on the **aggregate** flow across each interface rather than a path-by-path limit: ``flow_12`` caps exports out of ``region_1`` into ``region_2``, ``flow_21`` caps imports in the opposite direction. Only the import/export ``Link`` components added by ``add_extra_components`` are constrained, so the caps are inert when ``imports``/``exports`` are disabled." ,,, co2limit_enable,bool,true or false,"Switch to activate the system-wide CO2 cap below. Optional; defaults to false when unset. Can also be set via the ``Co2L`` keyword in the ``{opts}`` wildcard." co2limit,:math:`t_{CO_2}/a`,float,"System-wide cap on annual CO2 emissions, added as a global constraint in ``prepare_network``. Only applied when ``co2limit_enable`` is true." @@ -33,7 +33,7 @@ gaslimit,MWh thermal,float,"Cap on annual gas-fired primary energy from gas carr ,,, demand:,,, -- bus_allocation,--,"One of {``population``, ``breakthrough``}","How zone-level demand is distributed to individual buses. ``population`` (default) weights buses by 2020 Decennial Census county populations (split evenly across each county's substations, then each substation's buses). ``breakthrough`` uses the legacy nominal-demand column (``Pd``) from the 2016-vintage Breakthrough Energy grid model." --- profile,--,"One of {``efs``, ``eia``, ``eer``, ``servm``}","Datasource for electrical load data. ``EFS`` pulls future state level electrical demand data. ``EIA`` pulls historical balancing level electrical demand data. ``EER`` pulls future state-level profiles from the EER dataset; when selected, ``planning_horizons`` must be one of 2021, 2025, 2030, 2035, 2040, 2045, or 2050 and ``renewable_weather_years`` must contain exactly one year from 2007-2013 or 2016-2023. ``SERVM`` pulls CPUC SERVM hourly load for the six California load regions (California models only); when selected, ``planning_horizons`` must be one of 2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, or 2045." +-- profile,--,"One of {``efs``, ``eia``, ``eer``, ``servm``}","Datasource for electrical load data. ``EFS`` pulls future state level electrical demand data. ``EIA`` pulls historical balancing level electrical demand data. ``EER`` pulls future state-level profiles from the EER dataset; when selected, ``planning_horizons`` must be one of 2021, 2025, 2030, 2035, 2040, 2045, or 2050 and ``renewable_weather_years`` must contain exactly one year from 2007-2013 or 2016-2023. ``SERVM`` pulls CPUC SERVM hourly load for the six California load regions (California models only); when selected, ``planning_horizons`` must be one of 2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, or 2045, and each horizon must carry a full 8760-hour snapshot year. Only the ``Net Load`` component is dispatched against; the full component split is written to ``power_zonal_components_s{simpl}.parquet``. See the SERVM section of the demand data page." -- scenario:,,, -- -- efs_case,--,"One of {``reference``, ``medium``, ``high``}",(UNDER DEVELOPMENT) Extracts EFS data according to level of adoption -- -- efs_speed,--,"One of {``slow``, ``moderate``, ``fast``}",(UNDER DEVELOPMENT) Extracts EFS data according to speed of electrification diff --git a/docs/source/data-demand.md b/docs/source/data-demand.md index 89be90f8..36dc28fe 100644 --- a/docs/source/data-demand.md +++ b/docs/source/data-demand.md @@ -21,6 +21,98 @@ the model years 2021, 2025, 2030, 2035, 2040, 2045, and 2050, and each profile i historical weather year, so `renewable_weather_years` must contain exactly one year from 2007-2013 or 2016-2023. +(servm-demand)= +### CPUC SERVM (California) + +`profile: servm` uses the hourly load forecast the California Public Utilities Commission +publishes for its 2026 Integrated Resource Planning cycle, produced with the SERVM +production-cost model. It is a **California-only** dataset — use it with a footprint scoped +to California (`model_topology: include: reeds_state: ['CA']`); the maintained entry point is +`workflow/repo_data/config/config.california.yaml`. + +The workflow retrieves one CSV per forecast year from + +``` +https://files.cpuc.ca.gov/energy/modeling/2026_servm_updates/HourlyLoad_CA_Regions_V2025E_2224_Mon_{year}.csv +``` + +(~118 MB each, via `retrieve_cpuc_servm_load`). **Nine forecast years are published: 2026, +2028, 2030, 2032, 2035, 2037, 2040, 2042, and 2045.** Unlike EFS, SERVM demand is *not* +interpolated or AEO-scaled between published years, so `scenario: planning_horizons` must be +drawn from that set — any other year raises in `ReadServm`. + +#### Regions + +SERVM reports six California load regions. Each maps onto the balancing areas PyPSA-USA +carries on its buses (`workflow/repo_data/CPUC/servm_region_map.csv`): + +| SERVM region | PyPSA-USA balancing area(s) | Notes | +| --- | --- | --- | +| `PGE` | `CISO-PGAE` | PG&E CAISO footprint | +| `SCE` | `CISO-SCE` | Includes Valley Electric Association's California load, per the CPUC data dictionary | +| `SDGE` | `CISO-SDGE` | San Diego Gas & Electric | +| `IID` | `IID` | Imperial Irrigation District | +| `LADWP` | `LDWP` | Los Angeles Department of Water and Power | +| `NCNC` | `BANC` + `TIDC` | Northern California non-CAISO: Balancing Authority of Northern California and Turlock Irrigation District | + +`CISO-VEA` — the Valley Electric Association balancing area — is a Nevada footprint and is +deliberately **excluded** from California-only networks. It carries an empty region in the +mapping file, so its (small) load share is dropped with a log message, while any *unknown* +balancing area introduced by upstream relabeling still hard-fails. + +Because the SERVM regions are balancing areas rather than states, the demand is disaggregated +with a purpose-built weights table (`build_servm_load_weights`) instead of the generic +state/BA path: a cluster bus can straddle two SERVM regions (Los Angeles County holds both +`LDWP` and `CISO-SCE` buses), so a bus receives the sum of its share of every region it +overlaps. The underlying per-bus weights are still the `bus_allocation` weights described +below. + +#### Components + +Each file publishes several load components per region (`Load`, `BTMPV`, `EV`, `DATA_CEN`, ...) +alongside `Net Load`. **Only `Net Load` is dispatched against by the model** — it is what +remains after behind-the-meter PV and other embedded resources. Every published component is +nonetheless carried through on the `subsector` index level and written to the zonal artifact +`resources//demand/{interconnect}/power_zonal_components_s{simpl}.parquet`, so the +component split stays available for reporting. Components that exist for only some regions +(`EV` and friends are published for PGE/SCE/SDGE only) align to `NaN` there. + +#### Weather years + +Each forecast-year file stacks 25 weather years (2000-2024) of a full hourly year. +`electricity: demand: scenario: servm_weather_years:` selects which one to use. It takes a +list with **exactly one** entry; multiple entries are reserved for stochastic scenarios +(phase 3) and currently raise `NotImplementedError`, because the demand output path is not +weather-year specific. + +**Set `servm_weather_years` equal to the top-level `renewable_weather_years`.** Drawing load +and wind/solar profiles from different weather years decorrelates them and will understate +both the peak-net-load and the flexibility need. A mismatch is permitted but logs a warning. + +#### Timezone and calendar caveats + +SERVM strips are in **fixed Pacific Standard Time (UTC−8) with no daylight-saving +transition**. This is not stated in the source files; it was verified empirically from the +behind-the-meter PV solar-noon centroid, which sits at hour 12.52 in December and 12.68 in +July — a DST-observing series would move by a full hour between the two. The strips are +rolled forward 8 hours to UTC before being attached. + +Two calendar misalignments are accepted, and are immaterial for an hourly +capacity-expansion model, but matter if you compare hour-for-hour against another source: + +1. **Monday-start synthetic calendar.** SERVM lays each year's 8760 hours on a synthetic + calendar that starts on a Monday, so weekday-versus-weekend hours do not line up with the + real weekdays of the planning horizon. +2. **Leap weather years.** For a leap *weather* year the SERVM strip contains February 29 and + omits December 31, while PyPSA-USA's snapshots do the opposite (`get_snapshots` drops + February 29 from leap planning horizons). Every hour after February therefore lands one + calendar day earlier than it sat in the source file. + +The strip is mapped **positionally** onto the network's own per-period snapshots rather than +onto a synthesised `date_range` — the latter would run a day short of December 31 for the +leap planning horizons (2028, 2032, 2040). As a consequence each planning horizon must carry +exactly 8760 snapshots; a truncated snapshot window cannot be used with `profile: servm`. + ## Demand Disaggregation All of the demand sources above arrive at a coarser resolution than the network: EIA930 @@ -68,8 +160,9 @@ horizons setting, and the electricity demand setting. If conducting historical s user must select a planning horizon in the past (2018-2023) and set `profile: eia`. If conducting forward-looking planning cases the user must set a future planning horizon — -2030, 2040, or 2050 with `profile: efs`, or any of 2021, 2025, 2030, 2035, 2040, 2045, and -2050 with `profile: eer`. +2030, 2040, or 2050 with `profile: efs`; any of 2021, 2025, 2030, 2035, 2040, 2045, and +2050 with `profile: eer`; or any of 2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, and 2045 +with `profile: servm` (California only). For planning horizons between the EFS data years, PyPSA-USA implements a scaling factor that interpolates between future years or scales historical demand using forecasts from the Annual @@ -81,11 +174,12 @@ scenario: electricity: demand: - profile: efs # efs, eia, eer + profile: efs # efs, eia, eer, servm scenario: efs_case: reference # reference, medium, high efs_speed: moderate # slow, moderate, rapid eer_file: demand_EER2025_100by2050.h5 # used when profile: eer + servm_weather_years: [2019] # used when profile: servm; exactly one year, 2000-2024 aeo: reference ``` diff --git a/docs/source/data-transmission.md b/docs/source/data-transmission.md index 4930c6df..63761b9f 100644 --- a/docs/source/data-transmission.md +++ b/docs/source/data-transmission.md @@ -38,6 +38,37 @@ While representative of the US electricity system, the TAMU network is synthetic See the [Spatial Configuration](./config-spatial.md) page for information on how to choose between networks. ``` +## Interface Transmission Limits + +The path-by-path ratings above are complemented by **interface** limits: aggregate MW caps on +the total simultaneous flow across a bundle of paths. PyPSA-USA ships the CPUC RESOLVE +interface table at `config/policy_constraints/transmission_interface_limits.csv`, which rates +the CAISO import/export capability against the rest of WECC: + +| Interface | `region_1` (inside) | `region_2` (outside) | `flow_12` (MW) | `flow_21` (MW) | +| --- | --- | --- | --- | --- | +| `CA_NW` | p9, p10, p11 | p2, p5, p6, p7, p8 | 3,592 | 9,269 | +| `CA_SW` | p9, p10, p11 | p12, p13, p25, p27, p28, p30 | 10,901 | 10,463 | +| `CAISO_Imports` | p9, p10, p11 | all of the above | 9,728 | 10,208 | + +**Flow orientation:** `flow_12` is the cap on flow *out of* `region_1` (exports), `flow_21` the +cap on flow *into* `region_1` (imports). Enable the table with +`model_topology: interface_transmission_limits: true`; the constraint formulation is described +in [Model Constraints](./model-constraints.md#interface-transmission-limits). + +```{warning} +The interface caps are applied to the virtual `imports` / `exports` links created by +`add_extra_components`, so they bind only when `electricity: imports` / `electricity: exports` +are enabled, and they only see flow that crosses the boundary of the modeled footprint. + +`p8` (northeastern California) appears in the `region_2` list of every RESOLVE row but is +itself a California zone. In a California-only model it is therefore *inside* the network, and +the internal `p8`-`p9` AC corridor — about 300 MW in the ReEDS/NARIS balancing-area table — +carries no trade links and escapes the `CAISO_Imports` cap. Simultaneous CAISO imports are +understated by roughly that amount. This gap is documented rather than corrected: closing it +would require constraining internal AC lines alongside the trade links. +``` + (transmission-data)= ### Data ```{eval-rst} diff --git a/docs/source/datatables/demand.csv b/docs/source/datatables/demand.csv index 18b33f81..ef55c509 100644 --- a/docs/source/datatables/demand.csv +++ b/docs/source/datatables/demand.csv @@ -1,3 +1,4 @@ Characteristic,Data Source,Spatial Scale,Temporal Scale Historical Demand,GridEmissions (EIA930),Balancing Area,Hourly (2018- 2023) Future Demand,NREL Electrification Futures Study (EFS),States,"Hourly (2030, 2040, 2050)" +Future Demand (California),CPUC SERVM 2026 IRP Hourly Load,Six California load regions (PGE/SCE/SDGE/IID/LADWP/NCNC),"Hourly (2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, 2045; weather years 2000-2024)" diff --git a/docs/source/model-constraints.md b/docs/source/model-constraints.md index ad510003..42443b67 100644 --- a/docs/source/model-constraints.md +++ b/docs/source/model-constraints.md @@ -52,6 +52,7 @@ Regions used by the policy constraints may be specified as state codes, ReEDS zo | [Bidirectional link coupling](#bidirectional-link-coupling) | Equal capacity expansion of paired forward/reverse links | always active | [bidirectional_link.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/opts/bidirectional_link.py) | | [Demand-response capacity](#demand-response-capacity) | Shifted load bounded by a fixed share of nominal load per bus and snapshot | `electricity: demand_response: shift` | [sector.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/opts/sector.py) | | [Import/export volume limits](#import-and-export-volume-limits) | Traded energy bounded by a share of demand per balancing period | `electricity: imports/exports: volume_limit` | [interchange.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/opts/interchange.py) | +| [Interface transmission limits](#interface-transmission-limits) | Aggregate MW cap on the total flow across a bundle of transmission paths | `model_topology: interface_transmission_limits`; `electricity: transmission_interface_limits` | [interfaces.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/opts/interfaces.py) | | [National emission cap](#national-emission-cap-co2l) | System-wide CO2 cap via PyPSA `GlobalConstraint` | `Co2L` opts token; `electricity: co2limit` | [prepare_network.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/prepare_network.py) | | [Natural gas limit](#natural-gas-limit-ch4l) | Cap on annual gas-fired primary energy | `CH4L` opts token; `electricity: gaslimit` | [prepare_network.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/prepare_network.py) | | [Emission pricing](#emission-pricing-ep) | CO2 price added to marginal costs (objective, not a constraint) | `Ep` opts token; `costs: emission_prices` | [prepare_network.py](https://github.com/PyPSA/pypsa-usa/blob/master/workflow/scripts/prepare_network.py) | @@ -307,6 +308,41 @@ percent, and {math}`d_t` total AC load in period {math}`\tau`. In sector studies measured as the flow into the end-use sectors and the bound becomes a linear constraint in both trade and demand variables. +(interface-transmission-limits)= +## Interface transmission limits + +A transmission *interface* is a bundle of paths between two groups of regions that is rated in +aggregate rather than path-by-path — CAISO's simultaneous import capability being the canonical +example. Setting `model_topology: interface_transmission_limits: true` reads the interface table +at `electricity: transmission_interface_limits` (columns +`interface, region_1, region_2, flow_12, flow_21`) and adds one per-snapshot constraint per +interface and direction: + +**Trigger:** `model_topology: interface_transmission_limits: true`. + +\begin{align*} + &\ \hspace{1cm} \sum_{\ell \in I^{\rightarrow}} p_{\ell,t} \;\leq\; F_{12} + \hspace{0.5cm} \forall_t + \hspace{1cm} + \sum_{\ell \in I^{\leftarrow}} p_{\ell,t} \;\leq\; F_{21} + \hspace{0.5cm} \forall_t +\end{align*} + +where {math}`F_{12}` (`flow_12`) caps flow **out of** `region_1` into `region_2` and +{math}`F_{21}` (`flow_21`) caps flow in the opposite direction. {math}`I^{\rightarrow}` and +{math}`I^{\leftarrow}` are selected by bus membership and carrier, never by link name. + +```{important} +Only the `imports` / `exports` links added by `add_extra_components` are constrained, so the +limits are a **no-op when `electricity: imports` and `electricity: exports` are both disabled**. +A `region_2` entry that is itself inside the modeled footprint contributes no trade links, so +internal AC lines between it and `region_1` escape the cap. In the shipped +`CAISO_Imports` row this applies to `p8`, which is a California zone: in a California-only run +the internal `p8`-`p9` corridor (~300 MW in the ReEDS/NARIS balancing-area table) is not +counted against the CAISO import cap, understating simultaneous imports by roughly that +corridor's rating. This is documented, not corrected. +``` + (national-emission-cap-co2l)= ## National emission cap (Co2L) diff --git a/docs/source/release-notes.md b/docs/source/release-notes.md index 4e2715a0..14dfe52b 100644 --- a/docs/source/release-notes.md +++ b/docs/source/release-notes.md @@ -49,6 +49,36 @@ The full engineering change-log, including per-change expected effects on model results, is maintained in the repository at [`docs/CHANGELOG-v1-epic.md`](https://github.com/PyPSA/pypsa-usa/blob/master/docs/CHANGELOG-v1-epic.md). +### California / CPUC SERVM + +- **New demand source `electricity: demand: profile: servm`** — CPUC SERVM 2026 IRP hourly + load for the six California load regions (PGE, SCE, SDGE, IID, LADWP, NCNC), retrieved + per forecast year from files.cpuc.ca.gov. Nine forecast years are published (2026, 2028, + 2030, 2032, 2035, 2037, 2040, 2042, 2045) and `planning_horizons` is restricted to them. + `electricity: demand: scenario: servm_weather_years` picks one weather year out of the + stacked 2000-2024 record. Only `Net Load` is dispatched; the full component split is + written to a new component-resolved zonal artifact + (`power_zonal_components_s{simpl}.parquet`), which is now produced for every demand + profile. See {ref}`servm-demand`. +- **SERVM load-allocation weights** — a new `build_servm_load_weights` rule composes the + base→substation→cluster busmaps into a fractional `(SERVM region, bus)` table, so a cluster + that straddles two regions (Los Angeles County holds both LDWP and CISO-SCE buses) receives + the sum of its share of each. +- **Interface transmission limits are live.** `model_topology: interface_transmission_limits` + and `electricity: transmission_interface_limits` were previously dead keys. They now apply + the RESOLVE interface table as a per-snapshot cap on the *aggregate* flow across each + interface. The constraint scopes to the import/export links, so it is inert when trade is + disabled; the resulting understatement for `region_2` entries inside the footprint (notably + `p8` in California-only runs) is documented in {doc}`data-transmission`. +- **New maintained config `config.california.yaml`** — a runnable California-only model on + SERVM demand with the CAISO interface caps and imports/exports enabled, at REeDS-zone + resolution (`clusters: 4`) with a commented county-resolution alternative + (`clusters: 58`, `simpl: county`). +- **Phase-2 hook `conventional: ambient_derate`** — reserved for CPUC SERVM unit-specific + ambient-temperature derates. It is not implemented; enabling it raises `NotImplementedError` + in `add_electricity`. When it lands it replaces the EIA-860 seasonal derate rather than + stacking on it. + ### Documentation - New Model Description section ({doc}`model-workflow`, {doc}`model-components`, diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index a2500105..ec9c4429 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -119,6 +119,112 @@ def busmap_s(self) -> Path: return self.base / "busmaps" / self.interconnect / f"busmap_s{self.simpl}.csv" +@dataclass(frozen=True) +class ServmArtifacts: + """Paths to the CPUC SERVM demand artifacts produced by the Tier B build. + + The ``interconnect`` and ``simpl`` defaults MUST match + ``workflow/repo_data/config/config.test.california.yaml``'s ``scenario`` + section. If you change one, change both. + """ + + run_name: str + base: Path # resources/{run_name}/ + interconnect: str = "western" + simpl: str = "20" + + @property + def elec_s(self) -> Path: + """Simplified network the weights and demand were built against.""" + return self.base / "networks" / self.interconnect / f"elec_s{self.simpl}.nc" + + @property + def elec_s_dem(self) -> Path: + """Simplified network with SERVM demand attached.""" + return self.base / "networks" / self.interconnect / f"elec_s{self.simpl}_dem.nc" + + @property + def weights(self) -> Path: + """(bus, servm_region, laf) allocation table from build_servm_load_weights.""" + return self.base / "demand" / self.interconnect / f"servm_load_weights_s{self.simpl}.csv" + + @property + def demand(self) -> Path: + """Per-bus hourly demand CSV written by build_electrical_demand.""" + return self.base / "demand" / self.interconnect / f"power_electricity_s{self.simpl}.csv" + + @property + def zonal_components(self) -> Path: + """Component-resolved zonal demand parquet (pre-disaggregation).""" + return self.base / "demand" / self.interconnect / f"power_zonal_components_s{self.simpl}.parquet" + + +def _run_snakemake(configfile: str, until: str, run_name: str) -> None: + """Run one ``snakemake --until `` build against ``workflow/``.""" + cmd = [ + "snakemake", + "--until", + until, + "--configfile", + configfile, + "--config", + f"run={{name: '{run_name}', shared_cutouts: true}}", + "-j", + str(os.cpu_count() or 2), + # Force greedy scheduler to avoid the ILP scheduler's cbc dependency + # (cbc is shipped non-executable in some envs and causes PermissionError). + "--scheduler", + "greedy", + "--quiet", + ] + try: + result = subprocess.run( + cmd, + cwd=WORKFLOW_DIR, + capture_output=True, + text=True, + timeout=600, + ) + except subprocess.TimeoutExpired as e: + raw_stderr = e.stderr + if isinstance(raw_stderr, bytes): + raw_stderr = raw_stderr.decode("utf-8", errors="replace") + stderr_tail = "\n".join(raw_stderr.splitlines()[-100:]) if raw_stderr else "" + pytest.fail( + f"snakemake build timed out after {e.timeout}s\nstderr (last 100 lines):\n{stderr_tail}", + ) + if result.returncode != 0: + pytest.fail( + f"snakemake build failed (exit {result.returncode}):\n" + f"stderr (last 100 lines):\n" + "\n".join(result.stderr.splitlines()[-100:]), + ) + + +@pytest.fixture(scope="session") +def servm_built(tmp_path_factory) -> ServmArtifacts: + """Run ``snakemake --until add_demand`` on the SERVM test config once per session. + + Stops at ``add_demand`` rather than ``cluster_network``: the SERVM + artifacts under test (weights, per-bus demand CSV, zonal components) are all + produced at or before that stage, and the extra clustering work is already + covered by the ``built`` fixture. + + Skips like ``built`` when the data dirs are missing. This build also + downloads one ~118 MB CPUC SERVM load file the first time it runs. + """ + missing = [d for d in DATA_DIRS if not d.exists()] + if missing: + pytest.skip( + "Integration tests require populated data dirs; missing: " + ", ".join(str(m) for m in missing), + ) + run_name = f"pytest_servm_{tmp_path_factory.mktemp('servm_run').name}" + _run_snakemake("config/config.test.california.yaml", "add_demand", run_name) + return ServmArtifacts( + run_name=run_name, + base=WORKFLOW_DIR / "resources" / run_name, + ) + + @pytest.fixture(scope="session") def built(tmp_path_factory) -> BuiltArtifacts: """Run ``snakemake --until cluster_network`` once per session and expose artifact paths. diff --git a/tests/integration/test_servm_demand_artifacts.py b/tests/integration/test_servm_demand_artifacts.py new file mode 100644 index 00000000..73d45138 --- /dev/null +++ b/tests/integration/test_servm_demand_artifacts.py @@ -0,0 +1,146 @@ +"""Tier B — assert shape and orientation of the CPUC SERVM demand artifacts. + +Built from ``workflow/repo_data/config/config.test.california.yaml`` by the +``servm_built`` fixture (``snakemake --until add_demand``). These are structural +checks — conservation of the allocation weights, alignment of the per-bus demand +CSV with the network it was built against, and the PST->UTC roll. No numerical +regression against the CPUC reference case here; that is the (separate) +benchmark rule ``run.benchmark_cpuc`` reserves. +""" + +from __future__ import annotations + +import pandas as pd +import pypsa +import pytest + +pytestmark = pytest.mark.integration + +# The six CPUC SERVM California load regions (ReadServm.REGIONS). +SERVM_REGIONS = {"IID", "LADWP", "NCNC", "PGE", "SCE", "SDGE"} + + +class TestServmLoadWeights: + """Assertions on servm_load_weights_s{simpl}.csv (build_servm_load_weights).""" + + def test_file_exists(self, servm_built): + """The weights table was produced by snakemake.""" + assert servm_built.weights.exists(), f"missing {servm_built.weights}" + + def test_columns(self, servm_built): + """The table carries the three columns WriteServm reads.""" + weights = pd.read_csv(servm_built.weights) + assert {"bus", "servm_region", "laf"}.issubset(weights.columns) + + def test_laf_sums_to_one_per_region(self, servm_built): + """Allocation factors conserve each region's demand exactly.""" + weights = pd.read_csv(servm_built.weights) + totals = weights.groupby("servm_region")["laf"].sum() + assert not totals.empty, "no SERVM regions in the weights table" + pd.testing.assert_series_equal( + totals, + pd.Series(1.0, index=totals.index, name="laf"), + check_exact=False, + rtol=1e-6, + ) + + def test_regions_are_known(self, servm_built): + """Only the six published SERVM regions appear (CISO-VEA is excluded).""" + weights = pd.read_csv(servm_built.weights) + unknown = set(weights.servm_region.unique()) - SERVM_REGIONS + assert not unknown, f"unexpected SERVM region(s) {sorted(unknown)}" + + def test_buses_exist_in_network(self, servm_built): + """Every weighted bus is a bus of the network the weights were built for.""" + weights = pd.read_csv(servm_built.weights) + n = pypsa.Network(str(servm_built.elec_s)) + unknown = set(weights.bus.astype(str)) - set(n.buses.index.astype(str)) + assert not unknown, f"weights reference {len(unknown)} bus(es) not in elec_s{servm_built.simpl}.nc" + + +class TestServmDemand: + """Assertions on power_electricity_s{simpl}.csv (build_electrical_demand).""" + + def test_file_exists(self, servm_built): + """The per-bus demand CSV was produced by snakemake.""" + assert servm_built.demand.exists(), f"missing {servm_built.demand}" + + def test_columns_are_network_buses(self, servm_built): + """Demand columns are buses of the network, and there is at least one.""" + demand = pd.read_csv(servm_built.demand, index_col=0) + n = pypsa.Network(str(servm_built.elec_s)) + assert len(demand.columns) > 0 + unknown = set(demand.columns.astype(str)) - set(n.buses.index.astype(str)) + assert not unknown, f"demand columns are not network buses: {sorted(unknown)[:10]}" + + def test_row_count_matches_snapshots(self, servm_built): + """One row per network snapshot.""" + demand = pd.read_csv(servm_built.demand, index_col=0) + n = pypsa.Network(str(servm_built.elec_s)) + assert len(demand) == len(n.snapshots), ( + f"demand has {len(demand)} rows, network has {len(n.snapshots)} snapshots" + ) + + def test_demand_is_positive(self, servm_built): + """No NaNs, no negative load, and a non-trivial total.""" + demand = pd.read_csv(servm_built.demand, index_col=0) + assert not demand.isna().any().any() + assert (demand.to_numpy() >= 0).all(), "SERVM Net Load produced negative demand" + assert demand.to_numpy().sum() > 0 + + def test_attached_to_network(self, servm_built): + """add_demand attached the SERVM load to elec_s{simpl}_dem.nc.""" + n = pypsa.Network(str(servm_built.elec_s_dem)) + assert len(n.loads) > 0 + assert not n.loads_t.p_set.isna().any().any() + assert n.loads_t.p_set.to_numpy().sum() > 0 + + def test_peak_hour_consistent_with_pst(self, servm_built): + """The annual peak lands in the CA afternoon/evening once rolled to UTC. + + ``ReadServm`` rolls the fixed-PST strip forward by 8 hours, so the + snapshot index is effectively UTC. A California system peak sits in the + late afternoon local time; 15:00-20:00 PST maps to UTC hours 23-04. The + band is deliberately wide — it is testing the direction and magnitude of + the roll, not the exact peak hour. (Peak *date* is not asserted: SERVM + lays its hours on a synthetic Monday-start calendar, and a leap weather + year shifts every post-February hour by one calendar day.) + """ + demand = pd.read_csv(servm_built.demand, index_col=0, parse_dates=True) + peak_hour = demand.sum(axis=1).idxmax().hour + assert peak_hour in { + 23, + 0, + 1, + 2, + 3, + 4, + }, f"annual peak at UTC hour {peak_hour}; expected 23-04 for a PST-rolled California profile" + + +class TestServmZonalComponents: + """Assertions on power_zonal_components_s{simpl}.parquet.""" + + def test_file_exists(self, servm_built): + """The component-resolved zonal artifact was produced.""" + assert servm_built.zonal_components.exists(), f"missing {servm_built.zonal_components}" + + def test_regions_are_columns(self, servm_built): + """Columns are the SERVM regions; every region carries nonzero energy.""" + zonal = pd.read_parquet(servm_built.zonal_components) + assert set(zonal.columns) <= SERVM_REGIONS, f"unexpected column(s) {sorted(set(zonal.columns) - SERVM_REGIONS)}" + net_load = zonal.xs("Net Load", level="subsector") + assert (net_load.sum() > 0).all(), ( + f"zero annual energy in region(s) {sorted(net_load.columns[net_load.sum() <= 0])}" + ) + + def test_net_load_component_present(self, servm_built): + """``Net Load`` — the only component the model dispatches against — is kept.""" + zonal = pd.read_parquet(servm_built.zonal_components) + assert "Net Load" in zonal.index.get_level_values("subsector") + + def test_components_beyond_net_load_are_kept(self, servm_built): + """The zonal artifact stays component-resolved (BTMPV, EV, ... survive).""" + zonal = pd.read_parquet(servm_built.zonal_components) + components = set(zonal.index.get_level_values("subsector")) + assert components - {"Net Load"}, "zonal artifact collapsed to Net Load only" diff --git a/tests/static/test_dag_dryrun.py b/tests/static/test_dag_dryrun.py index 26310624..96459580 100644 --- a/tests/static/test_dag_dryrun.py +++ b/tests/static/test_dag_dryrun.py @@ -44,33 +44,55 @@ def _seed_runtime_configs(): shutil.copy2(src, dst) +# County-resolution overrides for the California config. `config.california.yaml` +# ships this same block commented out; snakemake's ``--config`` performs a +# recursive dict update, so only the listed sub-keys are replaced. California has +# 58 counties, and ``simpl: county`` selects the county-FIPS fast path in +# ``cluster_simpl``. +CALIFORNIA_COUNTY_OVERRIDE = [ + "scenario={clusters: [58], simpl: ['county']}", + "model_topology={topological_boundaries: 'county'}", +] + + @pytest.mark.fast @pytest.mark.parametrize( - "configfile,target", + "configfile,target,overrides", [ - ("config/config.tutorial.yaml", "cluster_network"), - ("config/config.tutorial.yaml", "solve_network"), - ("config/config.default.yaml", "cluster_network"), + ("config/config.tutorial.yaml", "cluster_network", []), + ("config/config.tutorial.yaml", "solve_network", []), + ("config/config.default.yaml", "cluster_network", []), + ("config/config.california.yaml", "cluster_network", []), + ("config/config.california.yaml", "solve_network", []), + ( + "config/config.california.yaml", + "cluster_network", + CALIFORNIA_COUNTY_OVERRIDE, + ), ], + ids=lambda v: "+".join(v) if isinstance(v, list) else v, ) -def test_snakemake_dryrun_resolves(configfile, target): +def test_snakemake_dryrun_resolves(configfile, target, overrides): + cmd = [ + "snakemake", + "-n", + "--configfile", + configfile, + "--until", + target, + "--quiet", + ] + if overrides: + cmd += ["--config", *overrides] result = subprocess.run( - [ - "snakemake", - "-n", - "--configfile", - configfile, - "--until", - target, - "--quiet", - ], + cmd, cwd=WORKFLOW_DIR, capture_output=True, text=True, timeout=60, ) assert result.returncode == 0, ( - f"snakemake -n failed for {configfile} --until {target}\n" + f"snakemake -n failed for {configfile} --until {target} (overrides={overrides})\n" f"stderr:\n{result.stderr}\n" f"stdout (last 50 lines):\n" + "\n".join(result.stdout.splitlines()[-50:]) ) diff --git a/workflow/repo_data/config/config.california.yaml b/workflow/repo_data/config/config.california.yaml new file mode 100644 index 00000000..0a80ae24 --- /dev/null +++ b/workflow/repo_data/config/config.california.yaml @@ -0,0 +1,487 @@ +# ==================================================================== +# PyPSA-USA — California Configuration (CPUC SERVM demand) +# ==================================================================== +# A maintained, runnable California-only capacity-expansion model. +# +# What it is: +# * Footprint : the four California REeDS zones (p8, p9, p10, p11), +# carved out of the Western Interconnection. +# * Demand : CPUC SERVM 2026 IRP hourly load for the six California +# load regions (IID, LADWP, NCNC, PGE, SCE, SDGE), +# disaggregated onto buses with the population-derived +# SERVM allocation weights (`build_servm_load_weights`). +# * Trade : imports/exports to the rest of WECC are enabled and +# capped in aggregate by the RESOLVE CAISO interface +# limits (`model_topology.interface_transmission_limits`). +# * Horizons : 2030 / 2035 / 2040 / 2045, solved with perfect foresight. +# +# How to run it: +# cd workflow +# uv run snakemake -j1 --configfile config/config.california.yaml +# # data model only (no solve): +# uv run snakemake data_model -j1 --configfile config/config.california.yaml +# +# The `retrieve_cpuc_servm_load` rule pulls one ~118 MB CSV per planning +# horizon from files.cpuc.ca.gov the first time this config is built. +# +# Option reference: docs/source/config-configuration.md +# Demand reference: docs/source/data-demand.md (SERVM section) +# Spatial reference: docs/source/config-spatial.md +# ==================================================================== + + +# ==================================================================== +# RUN — run identity and shared-resource toggles +# ==================================================================== +run: + name: "california" # resources/california/ and results/california/ + disable_progressbar: false + shared_resources: false # isolate resources per run.name + shared_cutouts: true # cutouts are large; share them across runs + validation: false + # Marks this run as a CPUC benchmark case. The rule that consumes the flag + # (comparison of model output against the CPUC SERVM/RESOLVE reference case) + # lands in a separate PR; the key is inert until then and is safe to leave on. + benchmark_cpuc: true + + +# ==================================================================== +# RENEWABLE DATASET — capacity-factor source for wind & solar +# ==================================================================== +renewable: + dataset: godeeep # atlite | godeeep + + +# ==================================================================== +# SCENARIO — workflow wildcards and planning horizons +# ==================================================================== +# planning_horizons MUST be drawn from the nine CPUC SERVM forecast years +# {2026, 2028, 2030, 2032, 2035, 2037, 2040, 2042, 2045}; any other year +# raises in `ReadServm`. One SERVM load file is downloaded per horizon. +scenario: + interconnect: [western] # CA is carved out of the Western Interconnection + planning_horizons: [2030, 2035, 2040, 2045] # SERVM forecast years only + clusters: [4] # p8, p9, p10, p11 — the four CA REeDS zones + simpl: [75] # resource resolution before transmission clustering + ll: [v1.0] # line-limit scenario; today's transfer capacity + opts: [REM-3h] # REM = regional CO2 limits, 3h = 3-hour resolution + scope: "total" # urban | rural | total + sector: "" # electricity-only + +foresight: 'perfect' # all horizons solved in one monolithic problem + + +# ==================================================================== +# MODEL TOPOLOGY — transmission backbone and zonal aggregation +# ==================================================================== +# `clusters` is pinned to 4 because the ReEDS zonal backbone cannot be +# clustered below the number of zones in the footprint, and California +# holds exactly four (p8, p9, p10, p11). +model_topology: + transmission_network: 'reeds' # ReEDS/NARIS zonal backbone + topological_boundaries: 'reeds_zone' # county | reeds_zone | state + interface_transmission_limits: true # apply the RESOLVE CAISO interface caps + include: + reeds_state: ['CA'] # equivalent to reeds_zone: ['p8','p9','p10','p11'] + aggregate: {} + +# -------------------------------------------------------------------- +# ALTERNATIVE: county-resolution California +# -------------------------------------------------------------------- +# To run California at county resolution instead of REeDS-zone resolution, +# replace the `scenario` and `model_topology` values above with the block +# below. California has 58 counties, so `clusters` must be 58 — the county +# NARIS interface table (transmission_capacity_init_AC_county_NARIS2024.csv) +# carries exactly 58 `p06xxx` nodes. `simpl: ['county']` uses the county-FIPS +# fast path in `cluster_simpl` so the resource layer is built directly on +# county boundaries; a numeric `simpl` >= 58 also works if you want more +# resource zones than transmission nodes. +# +# `add_extra_components` swaps to the county NARIS flowgate file automatically +# when `topological_boundaries: 'county'` — no other key needs changing. +# +# scenario: +# interconnect: [western] +# planning_horizons: [2030, 2035, 2040, 2045] +# clusters: [58] # 58 California counties +# simpl: ['county'] # county-FIPS fast path in cluster_simpl +# ll: [v1.0] +# opts: [REM-3h] +# scope: "total" +# sector: "" +# +# model_topology: +# transmission_network: 'reeds' +# topological_boundaries: 'county' +# interface_transmission_limits: true +# include: +# reeds_state: ['CA'] +# aggregate: {} +# +# See docs/source/config-spatial.md for the county workflow and the +# minimum-cluster table. +# -------------------------------------------------------------------- + + +# ==================================================================== +# ENABLE — top-level feature flags +# ==================================================================== +enable: + build_cutout: false # consume the prebuilt atlite cutout + + +# ==================================================================== +# SNAPSHOTS & WEATHER YEARS — temporal scope +# ==================================================================== +# `renewable_weather_years` and `electricity.demand.scenario.servm_weather_years` +# are deliberately the same year: SERVM publishes 25 stacked weather years +# (2000-2024) per forecast year, and drawing load and renewable profiles from +# different weather years decorrelates load from wind/solar. A mismatch is +# allowed but logs a warning in `build_electrical_demand`. +renewable_weather_years: [2019] + +snapshots: + start: "2019-01-01" + end: "2020-01-01" + inclusive: "left" # a full 8760-hour year + +renewable_scenarios: ["rcp85cooler"] # GODEEEP climate scenario + +renewable_snapshots: + start_month: 1 + start_day: 1 + end_month: 12 + end_day: 31 + end_inclusive: true + + +# ==================================================================== +# ELECTRICITY — generators, storage, demand, reserves, trade +# ==================================================================== +electricity: + conventional_carriers: [nuclear, oil, OCGT, CCGT, coal, geothermal, biomass, waste] + renewable_carriers: [onwind, offwind_floating, solar, hydro] + retirement: economic + + extendable_carriers: + Generator: [solar, onwind, offwind_floating, OCGT, CCGT, CCGT-95CCS, nuclear, hydrogen_ct] + StorageUnit: [4hr_battery_storage, 8hr_battery_storage] + Store: [] + Link: [] + + # ---------------- Reserves & resource adequacy ---------------- + SAFE_reservemargin: 0.14 + SAFE_regional_reservemargins: 'config/policy_constraints/SAFE_regional_prm.csv' + + operational_reserve: + activate: false + epsilon_load: 0.02 + epsilon_vres: 0.02 + contingency: 4000 + + erm: + all: 0.15 + + # ---------------- Policy & emissions ---------------- + regional_Co2_limits: 'config/policy_constraints/regional_Co2_limits.csv' + technology_capacity_targets: 'config/policy_constraints/technology_capacity_targets.csv' + portfolio_standards: 'config/policy_constraints/portfolio_standards.csv' + # Aggregate MW caps on the CAISO import/export interfaces (RESOLVE), applied + # when model_topology.interface_transmission_limits is true. + transmission_interface_limits: 'config/policy_constraints/transmission_interface_limits.csv' + + # ---------------- Demand ---------------- + demand: + profile: servm # CPUC SERVM 2026 IRP hourly load (California only) + bus_allocation: population # 2020 Census county populations + scenario: + # Weather year drawn from the stacked 2000-2024 SERVM record. Exactly one + # entry: multiple entries are reserved for stochastic scenarios (phase 3) + # and currently raise NotImplementedError. + servm_weather_years: [2019] + aeo: reference # unused by SERVM (one file per forecast year) + + demand_response: + shift: 0 + marginal_cost: 999999 + + # ---------------- Inter-regional trade ---------------- + # California is not an island: the rest of WECC is represented as virtual + # import/export links whose aggregate flow is capped by the RESOLVE CAISO + # interfaces above. `volume_limit: 25` caps annual imported/exported energy + # at 25% of total demand, roughly CAISO's historical net-import share. + imports: + enable: true + costs: wholesale # priced at monthly wholesale hub prices + co2_emissions: 0.428 # tCO2/MWh assigned to imported energy + capacity_limit: true # cap import power at historical maxima + volume_limit: 25 # % of total demand per balancing period + balancing_period: year + + exports: + enable: true + costs: wholesale + capacity_limit: true + volume_limit: 25 + balancing_period: year + + +# ==================================================================== +# CONVENTIONAL — unit-commitment and fuel-price overrides +# ==================================================================== +conventional: + unit_commitment: false + must_run: false + dynamic_fuel_price: + enable: false + pudl: true + wholesale: true + + # ---- PHASE 2 (NOT IMPLEMENTED) ---------------------------------- + # Ambient-temperature capacity derates for thermal units, taken from the + # CPUC SERVM unit-specific hourly derate profiles. `enable: true` currently + # raises NotImplementedError in add_electricity — the CPUC derate profiles + # are not retrieved or applied yet (see the reserved + # `retrieve_cpuc_thermal_derate` rule in workflow/rules/retrieve.smk). + # + # When it lands it REPLACES the EIA-860 summer/winter seasonal derate + # (`replaces_seasonal_derate: true`) rather than stacking on top of it. + # Stacking an ambient derate on a seasonal derate, or on a UCAP-derated + # capacity credit, double-counts the same thermal deficiency and must never + # be done. + ambient_derate: + enable: false # PHASE 2 — leave false + source: cpuc_servm # CPUC SERVM unit-specific hourly derates + replaces_seasonal_derate: true # never stack with the EIA-860 seasonal derate + + +# ==================================================================== +# LINES — AC transmission lines +# ==================================================================== +lines: + s_max_pu: 0.7 + s_nom_max: .inf + max_extension: 20000 + length_factor: 1.25 + + +# ==================================================================== +# LINKS — HVDC and controllable links +# ==================================================================== +links: + p_max_pu: 1.0 + p_nom_max: .inf + max_extension: 20000 + + +# ==================================================================== +# CO2 — sequestration storage and pipeline transport +# ==================================================================== +co2: + storage: false + network: + enable: false + capital_cost: 2736000 + marginal_cost: 4 + lifetime: 40 + discount_rate: 0.07 + + +# ==================================================================== +# DAC — Direct Air Capture +# ==================================================================== +dac: + enable: false + granularity: "node" + capital_cost: 6000000 + electricity_input: 2.5 + lifetime: 20 + discount_rate: 0.07 + + +# ==================================================================== +# COSTS — capex/opex scenarios and policy incentives +# ==================================================================== +costs: + atb: + model_case: "Market" + scenario: "Moderate" + overrides: + + aeo: + scenario: "reference" + + social_discount_rate: 0.02 + ng_fuel_year: 2019 # CAISO NG price vintage year + + emission_prices: + enable: false + co2: 0. + co2_monthly_prices: false + + ptc_modifier: + onwind: 27.50 + biomass: 27.50 + + itc_modifier: + solar: 0.3 + offwind: 0.3 + offwind_floating: 0.3 + EGS: 0.3 + geothermal: 0.3 + SMR: 0.3 + nuclear: 0.3 + hydro: 0.3 + 2hr_battery_storage: 0.3 + 4hr_battery_storage: 0.3 + 6hr_battery_storage: 0.3 + 8hr_battery_storage: 0.3 + 10hr_battery_storage: 0.3 + 8hr_PHS: 0.3 + 10hr_PHS: 0.3 + 12hr_PHS: 0.3 + + min_year: + hydrogen_ct: 2040 + max_growth: + + +# ==================================================================== +# CLUSTERING — pre-cluster (simpl) and final (cluster) settings +# ==================================================================== +clustering: + simplify_network: + weighting_strategy: demand-capacity + algorithm: kmeans + cluster_network: + weighting_strategy: demand-capacity + algorithm: kmeans + exclude_carriers: [] + consider_efficiency_classes: false + + aggregation_strategies: + generators: + build_year: 'capacity_weighted_average' + lifetime: 'capacity_weighted_average' + start_up_cost: 'capacity_weighted_average' + min_up_time: 'capacity_weighted_average' + min_down_time: 'capacity_weighted_average' + ramp_limit_up: max + ramp_limit_down: max + committable: any + vom_cost: mean + fuel_cost: mean + heat_rate: mean + + temporal: + resolution_elec: false # the 3h in `opts` does the temporal reduction + resolution_sector: false + +focus_weights: + + +# ==================================================================== +# SOLVING — solver selection, options, and per-solver tunings +# ==================================================================== +solving: + options: + load_shedding: false + clip_p_max_pu: 1.e-2 + noisy_costs: true + skip_iterations: true + rolling_horizon: false + seed: 123 + track_iterations: false + min_iterations: 4 + max_iterations: 6 + transmission_losses: 2 + linearized_unit_commitment: true + horizon: 8760 + assign_all_duals: true + + solver: + name: gurobi + options: gurobi-default + + solver_options: + highs-default: + # https://ergo-code.github.io/HiGHS/options/definitions.html + threads: 4 + solver: "ipm" + run_crossover: "off" + small_matrix_value: 1e-6 + large_matrix_value: 1e9 + primal_feasibility_tolerance: 1e-5 + dual_feasibility_tolerance: 1e-5 + ipm_optimality_tolerance: 1e-4 + parallel: "on" + random_seed: 123 + gurobi-default: + threads: 8 + method: 2 + crossover: 0 + BarHomogeneous: 1 + BarConvTol: 1.e-5 + OptimalityTol: 1.e-4 + FeasibilityTol: 1.e-3 + ScaleFlag: 1 + Seed: 123 + AggFill: 0 + PreDual: 0 + GURO_PAR_BARDENSETHRESH: 200 + gurobi-numeric-focus: + name: gurobi + NumericFocus: 3 + method: 2 + crossover: 0 + BarHomogeneous: 1 + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 + ObjScale: -0.5 + threads: 8 + Seed: 123 + gurobi-fallback: + name: gurobi + crossover: 0 + method: 2 + BarHomogeneous: 1 + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-5 + OptimalityTol: 1.e-5 + Seed: 123 + threads: 8 + cplex-default: + threads: 4 + lpmethod: 4 + solutiontype: 2 + barrier.convergetol: 1.e-5 + feasopt.tolerance: 1.e-6 + cbc-default: {} + glpk-default: {} + + mem: 30000 + walltime: "12:00:00" + + +# ==================================================================== +# WALLTIME — per-rule HPC walltime overrides (only used by SLURM driver) +# ==================================================================== +walltime: + build_renewable_profiles: '04:00:00' + build_fuel_prices: '00:20:00' + add_demand: '02:00:00' + add_electricity: '04:00:00' + aggregate_to_substations: '02:00:00' + cluster_resources: '05:00:00' + cluster_network: '04:00:00' + solve_network: '20:00:00' + + +# ==================================================================== +# CUSTOM FILES — bring-your-own network or costs +# ==================================================================== +custom_files: + activate: false + files_path: '' + network_name: '' diff --git a/workflow/repo_data/config/config.default.yaml b/workflow/repo_data/config/config.default.yaml index 99fcc23d..4af707e9 100644 --- a/workflow/repo_data/config/config.default.yaml +++ b/workflow/repo_data/config/config.default.yaml @@ -215,6 +215,20 @@ conventional: pudl: true # PUDL receipts-based fuel cost time series (monthly resolution) wholesale: true # CAISO/wholesale natural-gas hub prices (overrides PUDL for NG) + # PHASE 2 (NOT IMPLEMENTED): ambient-temperature capacity derates for thermal + # units, from the CPUC SERVM unit-specific hourly derate profiles. Setting + # `enable: true` raises NotImplementedError in add_electricity — the profiles + # are not retrieved or applied yet (see the reserved + # `retrieve_cpuc_thermal_derate` rule in workflow/rules/retrieve.smk). + # When implemented it REPLACES the EIA-860 summer/winter seasonal derate + # rather than stacking on it; stacking an ambient derate on a seasonal derate + # (or on a UCAP-derated capacity credit) double-counts the same thermal + # deficiency and must never be done. + ambient_derate: + enable: false # PHASE 2 — leave false + source: cpuc_servm # CPUC SERVM unit-specific hourly derates + replaces_seasonal_derate: true # never stack with the EIA-860 seasonal derate + # ==================================================================== # LINES — AC transmission lines diff --git a/workflow/repo_data/config/config.test.california.yaml b/workflow/repo_data/config/config.test.california.yaml new file mode 100644 index 00000000..66648931 --- /dev/null +++ b/workflow/repo_data/config/config.test.california.yaml @@ -0,0 +1,135 @@ +# PyPSA-USA SERVM test config — used by tests/integration to build the +# CPUC SERVM demand artifacts on a minimal CA-only Western slice. +# DO NOT use as a user-facing entry point; see config.california.yaml +# for the maintained California model. +# +# NOTE ON SNAPSHOTS: unlike config.test.yaml this config keeps a full +# 8760-hour year. `ReadServm._assign_snapshots` maps each SERVM hourly +# strip positionally onto the network's own per-period snapshots and +# raises unless the planning horizon carries exactly 8760 of them, so a +# truncated snapshot window cannot build SERVM demand. Runtime is kept +# down by the small `simpl`/`clusters` values instead. +# +# NOTE ON DATA: building this config downloads one ~118 MB CPUC SERVM +# load file (the 2030 forecast year) from files.cpuc.ca.gov. +run: + name: "test-california" + disable_progressbar: true + shared_resources: false + shared_cutouts: true + validation: false + +foresight: 'perfect' + +scenario: + interconnect: [western] + clusters: [4m] + simpl: [20] + opts: [REM-3h] + ll: [v1.0] + scope: "total" + sector: "" + planning_horizons: [2030] # must be a CPUC SERVM forecast year + +model_topology: + transmission_network: 'reeds' + topological_boundaries: 'reeds_zone' + interface_transmission_limits: true + include: + reeds_state: ['CA'] + aggregate: {} + +enable: + build_cutout: false + +renewable_weather_years: [2019] + +snapshots: + start: "2019-01-01" + end: "2020-01-01" + inclusive: "left" + +electricity: + conventional_carriers: [nuclear, oil, OCGT, CCGT, coal, geothermal, biomass, waste] + renewable_carriers: [onwind, solar, hydro] + retirement: economic + extendable_carriers: + Generator: [solar, onwind, OCGT, CCGT] + StorageUnit: [4hr_battery_storage] + Store: [] + Link: [] + transmission_interface_limits: 'config/policy_constraints/transmission_interface_limits.csv' + demand: + profile: servm + bus_allocation: population + scenario: + servm_weather_years: [2019] + aeo: reference + imports: + enable: true + costs: wholesale + co2_emissions: 0.428 + capacity_limit: true + volume_limit: 25 + balancing_period: year + exports: + enable: true + costs: wholesale + capacity_limit: true + volume_limit: 25 + balancing_period: year + +conventional: + unit_commitment: false + must_run: false + dynamic_fuel_price: + enable: false + pudl: true + wholesale: true + ambient_derate: + enable: false # PHASE 2 — raises NotImplementedError when true + +lines: + s_max_pu: 0.7 + s_nom_max: .inf + max_extension: 20000 + length_factor: 1.25 + +links: + p_max_pu: 1.0 + p_nom_max: .inf + max_extension: 20000 + +costs: + ng_fuel_year: 2019 + +clustering: + simplify_network: + algorithm: kmeans + cluster_network: + algorithm: kmeans + exclude_carriers: [] + consider_efficiency_classes: false + aggregation_strategies: + generators: + build_year: 'capacity_weighted_average' + lifetime: 'capacity_weighted_average' + start_up_cost: 'capacity_weighted_average' + min_up_time: 'capacity_weighted_average' + min_down_time: 'capacity_weighted_average' + ramp_limit_up: max + ramp_limit_down: max + committable: any + vom_cost: mean + fuel_cost: mean + heat_rate: mean + temporal: + resolution_elec: false + resolution_sector: false + +focus_weights: + +custom_files: + activate: false + files_path: '' + network_name: '' diff --git a/workflow/rules/retrieve.smk b/workflow/rules/retrieve.smk index 290fa245..2e67902b 100644 --- a/workflow/rules/retrieve.smk +++ b/workflow/rules/retrieve.smk @@ -137,6 +137,12 @@ rule retrieve_cpuc_baseline_generators: "../scripts/retrieve_cpuc_data.py" +# RESERVED (phase 2): `retrieve_cpuc_thermal_derate` will pull the CPUC SERVM +# unit-specific ambient-temperature derate profiles that back the +# `conventional.ambient_derate` config hook. Until it lands, enabling that key +# raises NotImplementedError in add_electricity. + + sector_datafiles = [ # heating sector "population/DECENNIALDHC2020.P1-Data.csv", diff --git a/workflow/scripts/add_electricity.py b/workflow/scripts/add_electricity.py index 56324251..c6aade7b 100755 --- a/workflow/scripts/add_electricity.py +++ b/workflow/scripts/add_electricity.py @@ -988,6 +988,22 @@ def broadcast_investment_horizons_index(n: pypsa.Network, df: pd.DataFrame): return df +def check_ambient_derate_not_enabled(conventional: dict) -> None: + """Guard the phase-2 ``conventional.ambient_derate`` hook. + + Ambient-temperature derates are meant to REPLACE the EIA-860 seasonal + derate applied by :func:`apply_seasonal_capacity_derates`, never to stack + on top of it (or on top of a UCAP-derated capacity credit) — stacking + double-counts the same thermal deficiency. Until the CPUC SERVM + unit-specific hourly derate profiles are retrieved and applied, enabling + the key is a hard error rather than a silent no-op. + """ + if (conventional or {}).get("ambient_derate", {}).get("enable", False): + raise NotImplementedError( + "conventional.ambient_derate is a phase-2 hook; CPUC unit-specific hourly derate profiles are not yet retrieved or applied.", + ) + + def apply_seasonal_capacity_derates( n: pypsa.Network, plants: pd.DataFrame, @@ -1255,6 +1271,7 @@ def main(snakemake): renewable_carriers, unit_commitment=params.conventional["unit_commitment"], ) + check_ambient_derate_not_enabled(params.conventional) apply_seasonal_capacity_derates( n, plants, From 1e9fcdc6d9fc21ff5ea3aea9188adf87f2e92741 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:49:50 -0700 Subject: [PATCH 05/23] Add CPUC Baseline Generator List capacity benchmark (PR 5) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Compares installed capacity between the CPUC Baseline Generator List and the model's generator fleet, aggregated by SERVM benchmark region and technology category, per planning horizon. Emits a long-format comparison CSV and a per-horizon deviation heatmap under the run's figures/benchmark/ directory. Gated by the new run.benchmark_cpuc flag (default false), which is independent of the demand profile: the fleet benchmark is useful for any California run. Region attribution: EIA reports every CAISO plant under the single BA code CISO, so powerplants.csv cannot separate PGE from SCE from SDGE. The benchmark therefore runs at the coarsest resolution both sides support — CAISO (= PGE + SCE + SDGE), LADWP (LDWP), IID, NCNC (BANC + TIDC) — with the collapse stated in repo_data/CPUC/servm_benchmark_regions.csv. The model side is restricted to state == CA so EIA's CISO code does not drag in the Nevada CISO-VEA footprint. Technology attribution: repo_data/CPUC/servm_tech_map.csv maps both SERVM tech categories and PyPSA carriers into a shared compare_category. Anything unmapped on either side becomes an explicit UNMAPPED: row rather than being dropped. Vintage/retirement filtering is shared by both sides and mirrors add_electricity.load_powerplants, so the benchmark measures the fleet the model actually builds. Co-Authored-By: Claude Fable 5 --- docs/source/configtables/run.csv | 1 + workflow/Snakefile | 23 + .../CPUC/servm_benchmark_regions.csv | 8 + workflow/repo_data/CPUC/servm_tech_map.csv | 52 ++ workflow/repo_data/config/config.default.yaml | 1 + workflow/rules/validate.smk | 29 + workflow/scripts/benchmark_cpuc_baseline.py | 626 ++++++++++++++++++ .../test/test_benchmark_cpuc_baseline.py | 474 +++++++++++++ 8 files changed, 1214 insertions(+) create mode 100644 workflow/repo_data/CPUC/servm_benchmark_regions.csv create mode 100644 workflow/repo_data/CPUC/servm_tech_map.csv create mode 100644 workflow/scripts/benchmark_cpuc_baseline.py create mode 100644 workflow/scripts/test/test_benchmark_cpuc_baseline.py diff --git a/docs/source/configtables/run.csv b/docs/source/configtables/run.csv index 0cce5f09..4986acc8 100644 --- a/docs/source/configtables/run.csv +++ b/docs/source/configtables/run.csv @@ -4,3 +4,4 @@ disable_progrssbar,bool,"{true, false}",Switch to select whether progressbar sho shared_resources,bool,"{true, false}",Switch to select whether resources should be shared across runs. shared_cutouts,bool,"{true, false}",Switch to select whether cutouts should be shared across runs. validation,bool,"{true, false}",Switch to enable back-casting validation plotting +benchmark_cpuc,bool,"{true, false}",Switch to enable the CPUC Baseline Generator List capacity benchmark (California runs) diff --git a/workflow/Snakefile b/workflow/Snakefile index 809f06e6..35c8c0a5 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -258,6 +258,28 @@ def validation_figures(wildcards): return [] +BENCHMARK_CPUC_OUTPUTS = [ + "cpuc_capacity_benchmark.csv", + "cpuc_capacity_deviation.pdf", +] + + +def benchmark_figures(wildcards): + """CPUC Baseline Generator List capacity benchmark, opt-in via run.benchmark_cpuc. + + Gated on its own flag rather than on the demand profile: the fleet + benchmark is useful for any California run, whatever drives its load. + """ + if not config["run"].get("benchmark_cpuc", False): + return [] + return expand( + RESULTS + + "{interconnect}/figures/s{simpl}_cluster_{clusters}/l{ll}_{opts}_{sector}/benchmark/{figure}", + **config["scenario"], + figure=BENCHMARK_CPUC_OUTPUTS, + ) + + rule all: input: #"repo_data/dag.jpg", @@ -265,6 +287,7 @@ rule all: electricity_figures, sector_figures, validation_figures, + benchmark_figures, rule data_model: diff --git a/workflow/repo_data/CPUC/servm_benchmark_regions.csv b/workflow/repo_data/CPUC/servm_benchmark_regions.csv new file mode 100644 index 00000000..c273e14a --- /dev/null +++ b/workflow/repo_data/CPUC/servm_benchmark_regions.csv @@ -0,0 +1,8 @@ +servm_region,eia_ba_code,benchmark_region,note +PGE,CISO,CAISO,EIA reports every CAISO plant under the single BA code CISO; PGE/SCE/SDGE cannot be split without geography +SCE,CISO,CAISO,rolled into CAISO for the same reason +SDGE,CISO,CAISO,rolled into CAISO for the same reason +LADWP,LDWP,LADWP, +IID,IID,IID, +NCNC,BANC,NCNC,Balancing Authority of Northern California +NCNC,TIDC,NCNC,Turlock Irrigation District diff --git a/workflow/repo_data/CPUC/servm_tech_map.csv b/workflow/repo_data/CPUC/servm_tech_map.csv new file mode 100644 index 00000000..5f8bd6f8 --- /dev/null +++ b/workflow/repo_data/CPUC/servm_tech_map.csv @@ -0,0 +1,52 @@ +side,source_category,compare_category,note +cpuc,CC,Gas CC, +cpuc,CT,Gas CT/ICE/Steam, +cpuc,ICE,Gas CT/ICE/Steam, +cpuc,Steam,Gas CT/ICE/Steam,CA gas steam boilers; EIA maps 'Natural Gas Steam Turbine' to OCGT so it cannot be split out on the model side +cpuc,Cogen,Gas Cogen/CHP,CPUC-only bucket; EIA/PyPSA carriers have no CHP concept so these units land in CCGT/OCGT on the model side +cpuc,Coal,Coal, +cpuc,Nuclear,Nuclear, +cpuc,Geothermal,Geothermal, +cpuc,Biomass/Wood,Biomass, +cpuc,Biogas,Biomass, +cpuc,Hydro,Hydro, +cpuc,PSH,PSH, +cpuc,Wind,Wind, +cpuc,OOS_Wind,Wind,out-of-state wind contracted to a CA region; absent from the CA subset of the 2026 baseline +cpuc,Solar_1Axis,Solar, +cpuc,Solar_Fixed,Solar, +cpuc,Solar_2Axis,Solar, +cpuc,Solar_Thermal,Solar, +cpuc,Paired_Solar_1Axis,Solar, +cpuc,Paired_Solar_Fixed,Solar, +cpuc,Hybrid_Solar_1Axis,Solar, +cpuc,Hybrid_Solar_Fixed,Solar, +cpuc,Battery_4h,Battery, +cpuc,Paired_BattStorage,Battery, +cpuc,Hybrid_BattStorage,Battery, +cpuc,DR,Demand Response,CPUC-only; PyPSA-USA has no demand-response resource +cpuc,Pumping_Load,Pumping Load,CPUC-only; modelled as load rather than as a resource in PyPSA-USA +pypsa,CCGT,Gas CC, +pypsa,CCGT-95CCS,Gas CC, +pypsa,OCGT,Gas CT/ICE/Steam,EIA 'Natural Gas Steam Turbine' / 'Combustion Turbine' / 'Internal Combustion Engine' all map here +pypsa,hydrogen_ct,Gas CT/ICE/Steam, +pypsa,coal,Coal, +pypsa,nuclear,Nuclear, +pypsa,geothermal,Geothermal, +pypsa,EGS,Geothermal, +pypsa,biomass,Biomass, +pypsa,waste,Biomass,municipal solid waste; CPUC folds it into Biogas +pypsa,hydro,Hydro,pumped storage is split off before mapping via prime_mover_code == PS +pypsa,PHS,PSH, +pypsa,onwind,Wind, +pypsa,offwind,Wind, +pypsa,offwind_floating,Wind, +pypsa,solar,Solar, +pypsa,battery,Battery, +pypsa,2hr_battery_storage,Battery, +pypsa,4hr_battery_storage,Battery, +pypsa,6hr_battery_storage,Battery, +pypsa,8hr_battery_storage,Battery, +pypsa,10hr_battery_storage,Battery, +pypsa,oil,Oil,model-only; CPUC has no oil bucket for CA +pypsa,other,Other,"model-only catch-all (flywheels, other gases, compressed air)" diff --git a/workflow/repo_data/config/config.default.yaml b/workflow/repo_data/config/config.default.yaml index 99fcc23d..30de5ef2 100644 --- a/workflow/repo_data/config/config.default.yaml +++ b/workflow/repo_data/config/config.default.yaml @@ -34,6 +34,7 @@ run: shared_resources: false # true = reuse resources/ across all runs (DAG-keyed); false = isolate per run.name shared_cutouts: true # true = share atlite cutouts across runs (recommended; cutouts are large and slow) validation: false # true = wire in back-casting validation plots (historical-only) + benchmark_cpuc: false # true = compare installed capacity against the CPUC Baseline Generator List (California runs) # ==================================================================== # RENEWABLE DATASET — capacity-factor source for wind & solar diff --git a/workflow/rules/validate.smk b/workflow/rules/validate.smk index a6b968bb..98518516 100644 --- a/workflow/rules/validate.smk +++ b/workflow/rules/validate.smk @@ -66,3 +66,32 @@ rule plot_validation_figures: mem_mb=5000, script: "../scripts/plot_validation_production.py" + + +# Compares installed capacity against the CPUC Baseline Generator List, +# aggregated by SERVM benchmark region and technology, per planning horizon. +# Gated by run.benchmark_cpuc (see benchmark_figures in the Snakefile). +rule benchmark_cpuc_baseline: + params: + planning_horizons=config_provider("scenario", "planning_horizons"), + input: + network=NETWORKS + + "{interconnect}/elec_s{simpl}_c{clusters}_ec_l{ll}_{opts}_{sector}.nc", + powerplants="resources/powerplants/powerplants.csv", + cpuc_baseline=DATA + "cpuc/BaselineGeneratorList_CAISO.xlsx", + region_map="repo_data/CPUC/servm_benchmark_regions.csv", + tech_map="repo_data/CPUC/servm_tech_map.csv", + output: + comparison=RESULTS + + "{interconnect}/figures/s{simpl}_cluster_{clusters}/l{ll}_{opts}_{sector}/benchmark/cpuc_capacity_benchmark.csv", + heatmap=RESULTS + + "{interconnect}/figures/s{simpl}_cluster_{clusters}/l{ll}_{opts}_{sector}/benchmark/cpuc_capacity_deviation.pdf", + log: + LOGS + + "benchmark_cpuc_baseline/{interconnect}/elec_s{simpl}_c{clusters}_l{ll}_{opts}_{sector}.log", + threads: 1 + resources: + walltime="00:20:00", + mem_mb=5000, + script: + "../scripts/benchmark_cpuc_baseline.py" diff --git a/workflow/scripts/benchmark_cpuc_baseline.py b/workflow/scripts/benchmark_cpuc_baseline.py new file mode 100644 index 00000000..779711db --- /dev/null +++ b/workflow/scripts/benchmark_cpuc_baseline.py @@ -0,0 +1,626 @@ +# BY PyPSA-USA Authors +"""Benchmark the modelled California fleet against the CPUC Baseline Generator List. + +The CPUC publishes, alongside the SERVM hourly load, a unit-level baseline +generator list covering the whole WECC. Its California rows are the reference +fleet the CPUC's own SERVM/RESOLVE runs start from, so comparing installed +capacity against PyPSA-USA's fleet is a direct check on whether the model starts +a California study from a plausible resource mix. + +Two axes have to be reconciled before the sides are comparable. + +Region +------ +The CPUC tags every unit with a ``SERVM Region`` (PGE, SCE, SDGE, IID, LADWP, +NCNC). PyPSA-USA's ``powerplants.csv`` carries ``balancing_authority_code_eia``, +which is the *EIA* balancing-authority code — and EIA reports every CAISO plant +under the single code ``CISO`` with no sub-BA breakdown. There is therefore no +column on the model side that separates PGE from SCE from SDGE, and inventing +one (point-in-polygon against the clustered BA shapes) would need a +``resources/`` artifact that this rule does not consume. + +So the benchmark is run at the coarsest region resolution both sides support: + +* ``CAISO`` — CPUC PGE + SCE + SDGE vs. model BA ``CISO`` (the CPUC data + dictionary defines CAISO as exactly PGE + SCE + SDGE) +* ``LADWP`` — CPUC LADWP vs. model BA ``LDWP`` +* ``IID`` — CPUC IID vs. model BA ``IID`` +* ``NCNC`` — CPUC NCNC vs. model BAs ``BANC`` + ``TIDC`` + +The mapping lives in ``repo_data/CPUC/servm_benchmark_regions.csv`` so the +collapse is explicit and editable rather than hard-coded here. The model side is +additionally restricted to ``state == "CA"``: EIA's ``CISO`` also covers the +Valley Electric (CISO-VEA) footprint, which sits in Nevada and is excluded from +the SERVM California regions. + +Technology +---------- +SERVM splits technology far finer than PyPSA carriers do (nine solar buckets, a +paired/hybrid storage split, a CHP bucket). Both sides are therefore mapped into +a common ``compare_category`` by ``repo_data/CPUC/servm_tech_map.csv``, which +carries a ``side`` column (``cpuc`` / ``pypsa``) so one file documents both +directions. Anything unmapped on either side becomes ``UNMAPPED:`` and is +reported as its own row — a category is never silently dropped, because a silent +drop is exactly how a benchmark starts agreeing with itself. + +Two mappings deserve to be flagged when reading the output: + +* ``Gas Cogen/CHP`` is CPUC-only. EIA technology descriptions have no CHP + concept, so California gas cogeneration lands in ``CCGT``/``OCGT`` by prime + mover. Expect the model to be short in ``Gas Cogen/CHP`` and long in the two + gas buckets by roughly the same amount. +* ``Demand Response`` and ``Pumping Load`` are CPUC-only resources with no + PyPSA-USA counterpart; they are kept as rows with ``model_mw == 0``. + +Vintage and retirement +---------------------- +Both sides are filtered by the same rule (:func:`filter_active`): in service by +December 31 of the horizon, and not retired by then. The model side reproduces +``add_electricity.load_powerplants`` exactly (see +:func:`prepare_model_plants`) so this benchmark measures the fleet the model +actually builds, not a differently-filtered idealisation of it. +""" + +from __future__ import annotations + +import logging + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pypsa +import seaborn as sns +from _helpers import configure_logging + +logger = logging.getLogger(__name__) + +sns.set_theme("paper", style="whitegrid") + +DPI = 300 + +CPUC_SHEET = "BaselineGeneratorList" +#: The workbook's first row is blank; the real header is the second one. +CPUC_HEADER_ROW = 1 + +#: SERVM regions that make up California. Every other ``SERVM Region`` value in +#: the workbook is a non-CA WECC balancing area. +CA_SERVM_REGIONS = ("IID", "LADWP", "NCNC", "PGE", "SCE", "SDGE") + +CPUC_CAPACITY_COL = "Capmax MW" +CPUC_REGION_COL = "SERVM Region" +CPUC_TECH_COL = "SERVM Tech Category" +CPUC_INSV_COL = "Insvdt" +CPUC_RETIRE_COL = "RetireDate" + +#: Sentinel used by :func:`filter_active` for "no retirement date on record". +NEVER_RETIRES = pd.Timestamp("2262-01-01") + +UNMAPPED_PREFIX = "UNMAPPED:" + +#: Pseudo-region prefix for :func:`reconcile_totals`' informational rows. +RECONCILE_PREFIX = "RECONCILE:" + + +# --------------------------------------------------------------------------- # +# inputs +# --------------------------------------------------------------------------- # + + +def read_cpuc_baseline(xlsx_path: str) -> pd.DataFrame: + """ + Read the California rows of the CPUC Baseline Generator List. + + The workbook opens with a blank row, so the header is row 2 (``header=1``). + The published sheet covers all of WECC; only units whose ``SERVM Region`` is + one of :data:`CA_SERVM_REGIONS` are kept. + + ``Capmax MW`` is the nameplate capacity of the unit as SERVM sees it. Paired + and hybrid resources are published as separate rows per component (a + ``Paired_Solar_1Axis`` row and a ``Paired_BattStorage`` row), so summing + ``Capmax MW`` by technology counts each component once and never + double-counts the point of interconnection. + """ + df = pd.read_excel(xlsx_path, sheet_name=CPUC_SHEET, header=CPUC_HEADER_ROW) + + missing = [ + c + for c in (CPUC_CAPACITY_COL, CPUC_REGION_COL, CPUC_TECH_COL, CPUC_INSV_COL, CPUC_RETIRE_COL) + if c not in df.columns + ] + if missing: + raise ValueError( + f"{xlsx_path} sheet '{CPUC_SHEET}' is missing {missing}. The CPUC has changed the " + "workbook layout; update benchmark_cpuc_baseline.py's column constants.", + ) + + df = df.copy() + df[CPUC_REGION_COL] = df[CPUC_REGION_COL].astype("object").where(df[CPUC_REGION_COL].notna()).str.strip() + df[CPUC_TECH_COL] = df[CPUC_TECH_COL].astype("object").where(df[CPUC_TECH_COL].notna()).str.strip() + df[CPUC_CAPACITY_COL] = pd.to_numeric(df[CPUC_CAPACITY_COL], errors="coerce") + for col in (CPUC_INSV_COL, CPUC_RETIRE_COL): + df[col] = pd.to_datetime(df[col], errors="coerce") + + ca = df[df[CPUC_REGION_COL].isin(CA_SERVM_REGIONS)].copy() + if ca.empty: + raise ValueError( + f"No rows in {xlsx_path} carry a California SERVM Region ({', '.join(CA_SERVM_REGIONS)}). " + "Either the wrong sheet was read or the region labels have changed.", + ) + + logger.info( + "Read %d CPUC baseline units, %d of them in California across %d SERVM regions.", + len(df), + len(ca), + ca[CPUC_REGION_COL].nunique(), + ) + return ca + + +def load_tech_map(path: str) -> dict[str, dict[str, str]]: + """ + ``{side: {source_category: compare_category}}`` read from ``servm_tech_map.csv``. + + One file holds both directions, keyed by a ``side`` column of ``cpuc`` or + ``pypsa``, so the two halves of a comparison category can never drift apart + across separate files. + """ + df = pd.read_csv(path, dtype=str) + for col in ("side", "source_category", "compare_category"): + if col not in df.columns: + raise ValueError( + f"{path} must have columns side, source_category, compare_category; got {list(df.columns)}", + ) + df[col] = df[col].str.strip() + + unknown_sides = sorted(set(df.side.unique()) - {"cpuc", "pypsa"}) + if unknown_sides: + raise ValueError(f"{path} has unknown side values {unknown_sides}; expected 'cpuc' or 'pypsa'.") + + out: dict[str, dict[str, str]] = {} + for side, grp in df.groupby("side"): + duplicated = grp.source_category[grp.source_category.duplicated()].unique() + if len(duplicated): + raise ValueError(f"{path} maps {list(duplicated)} more than once on side '{side}'.") + out[side] = dict(zip(grp.source_category, grp.compare_category)) + + for side in ("cpuc", "pypsa"): + out.setdefault(side, {}) + return out + + +def load_benchmark_regions(path: str) -> tuple[pd.Series, pd.Series]: + """ + Return ``(servm_region -> benchmark_region, eia_ba_code -> benchmark_region)``. + + Both lookups come from the same file so the CAISO collapse (PGE/SCE/SDGE and + EIA ``CISO`` landing on one ``CAISO`` row) is stated once. + """ + df = pd.read_csv(path, dtype=str) + for col in ("servm_region", "eia_ba_code", "benchmark_region"): + if col not in df.columns: + raise ValueError( + f"{path} must have columns servm_region, eia_ba_code, benchmark_region; got {list(df.columns)}", + ) + df[col] = df[col].str.strip() + + servm = df.drop_duplicates("servm_region").set_index("servm_region").benchmark_region + dupes = servm.index[servm.index.duplicated()].unique() + if len(dupes): + raise ValueError(f"{path} maps SERVM regions {list(dupes)} to more than one benchmark region.") + + ba = df.drop_duplicates("eia_ba_code") + conflicting = ba.groupby("eia_ba_code").benchmark_region.nunique() + conflicting = conflicting[conflicting > 1] + if len(conflicting): + raise ValueError(f"{path} maps EIA BA codes {list(conflicting.index)} to more than one benchmark region.") + + return servm, ba.set_index("eia_ba_code").benchmark_region + + +# --------------------------------------------------------------------------- # +# shared filtering / mapping +# --------------------------------------------------------------------------- # + + +def filter_active( + df: pd.DataFrame, + horizon: int, + insv_col: str, + retire_col: str, +) -> pd.DataFrame: + """ + Units in service by the end of ``horizon`` and not yet retired. + + A unit counts if it enters service on or before December 31 of the horizon + year and its retirement year is strictly greater than the horizon — the same + inequality ``add_electricity.load_powerplants`` applies to the model fleet, + so the two sides of the benchmark are filtered identically. + + A missing retirement date means "never retires" (the workbook leaves the + field empty for units with no announced retirement; it also uses a + ``2050-12-01`` placeholder for many, which needs no special handling since a + 2050 retirement genuinely keeps the unit alive through any earlier horizon). + A missing in-service date is treated as already in service: the workbook is a + list of the *existing* baseline fleet, so an absent date is a data gap rather + than a future project. + """ + insv = pd.to_datetime(df[insv_col], errors="coerce") + retire = pd.to_datetime(df[retire_col], errors="coerce").fillna(NEVER_RETIRES) + + n_missing_insv = int(insv.isna().sum()) + if n_missing_insv: + logger.warning( + "%d of %d units have no in-service date; treating them as already in service.", + n_missing_insv, + len(df), + ) + + in_service = insv.isna() | (insv <= pd.Timestamp(year=horizon, month=12, day=31)) + not_retired = retire.dt.year > horizon + return df[in_service & not_retired] + + +def map_compare_category(values: pd.Series, mapping: dict[str, str]) -> pd.Series: + """ + Map source technology labels onto comparison categories. + + Anything absent from ``mapping`` becomes ``UNMAPPED: