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")