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)