From f5eb10f0f1bde3bde17af13aaf38f9262a1cf787 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:09:46 -0700 Subject: [PATCH] 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)