From 2f0b94835762b52ed4c9fb0b7829095ef6984e10 Mon Sep 17 00:00:00 2001 From: ktehranchi Date: Sat, 29 Aug 2026 15:26:57 -0700 Subject: [PATCH] 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")