Resiliency Module#

API reference for sdom.resiliency – the operational-resiliency evaluation module. See the user-guide page Resiliency Evaluation for the mathematical background, file-format requirements and worked examples.

Top-Level Convenience#

sdom.resiliency.evaluate_resiliency(snapshot_dir, *, inputs_dir, outage_spec, year=2030, scenario_id=1, n_hours=8760, hours=None, min_soc_per_tech=None, slack_penalty=10000.0, curtailment_penalty=0.0, formulation_overrides=None, n_workers=None, solver='highs', solver_options=None, profile_baseline=False, profile_outages=False)[source]#

End-to-end helper: load -> baseline dispatch -> outage evaluation.

Parameters:
  • snapshot_dir (str or pathlib.Path) – Snapshot directory passed to load_designed_system().

  • inputs_dir (str or pathlib.Path) – Previous-stage inputs directory.

  • outage_spec (OutageSpec) – Outage scenario specification.

  • year (int, optional) – Calendar year of the snapshot. Default 2030.

  • scenario_id (int, optional) – Scenario / Run id resolved from the snapshot CSVs. Default 1.

  • n_hours (int, optional) – Baseline-dispatch horizon length. Default 8760.

  • hours (iterable of int, optional) – Anchor hours to evaluate. None (default) evaluates every hour 1..n_hours.

  • min_soc_per_tech (dict, optional) – Operational SOC floor per storage tech (fraction of Cap_E); forwarded to both the baseline builder and the per-hour outage runner. Default None.

  • slack_penalty (float, optional) – Penalty (USD/MWh) applied to slack u[t] in the outage LP. Default 10_000.

  • curtailment_penalty (float, optional) – Penalty applied to curtailed VRE energy (USD/MWh). Default 0.

  • formulation_overrides (dict, optional) – Component formulation overrides forwarded to load_designed_system().

  • n_workers (int, optional) – Worker pool size for the per-hour evaluation. None (default) resolves to max(1, os.cpu_count() - 1) inside run_resiliency_evaluation().

  • solver (str, optional) – Pyomo solver name. "highs" first tries appsi_highs. Default "highs".

  • solver_options (dict, optional) – Solver options forwarded to both the baseline solve and every per-hour outage solve.

  • profile_baseline (bool, optional) – When True, attach a ModelInitProfiler to the baseline build/solve and print summary tables. Default False (opt-in to avoid runtime/logging overhead in the top-level helper).

  • profile_outages (bool, optional) – When True and n_workers == 1, profile every per-hour outage build. Ignored when n_workers > 1 (a per-worker summary is rarely useful). Default False (opt-in).

Returns:

Per-hour records (sorted by anchor hour) plus run metadata including n_workers_used, n_hours, solver and the outage_spec reference.

Return type:

ResiliencyResults

Examples

>>> from sdom.resiliency import OutageSpec, evaluate_resiliency
>>> spec = OutageSpec(
...     duration_hours=4,
...     recovery_hours=4,
...     outaged_assets={"imports": "all"},
... )
>>> results = evaluate_resiliency(
...     "snapshot/",
...     inputs_dir="inputs/",
...     outage_spec=spec,
...     n_hours=24,
...     hours=[1, 5, 10],
...     n_workers=1,
... )  

Data Loading#

sdom.resiliency.load_designed_system(snapshot_dir, *, inputs_dir, year=2030, scenario_id=1, formulation_overrides=None)[source]#

Load a fixed-capacity designed system from SDOM output snapshots.

Parameters:
  • snapshot_dir (str or pathlib.Path) – Directory containing OutputSummary_* and OutputSelectedVRE_* CSVs from a prior SDOM design run. OutputGeneration_* and OutputStorage_* files are ignored at this stage; files containing Phase1 in their names are also excluded.

  • inputs_dir (str or pathlib.Path) – Directory holding the previous-stage time-series and parameter CSVs (Load_hourly_*, CFSolar_*, CFWind_*, Data_Balancing_units_*, StorageData_*, Import_Cap_*, import_prices_*, Export_Cap_*, export_prices_*, fixed_dem_charges.csv, var_dem_charges.csv).

  • year (int, optional) – Calendar year used in CSV filenames. Default 2030.

  • scenario_id (int, optional) – Scenario / Run id to extract from snapshot CSVs. Default 1. If a snapshot file contains a single unique scenario, that one is used and a warning is emitted when scenario_id differs. If multiple scenarios are present, scenario_id must match one of them; otherwise ValueError is raised listing the available ids.

  • formulation_overrides (dict, optional) – Mapping {component: formulation_name} overlaid on the default formulation map.

Return type:

DesignedSystem

Raises:
  • FileNotFoundError – When the OutputSummary_* snapshot is missing.

  • ValueError – When scenario_id cannot be resolved against the snapshot.

Examples

>>> ds = load_designed_system(
...     "Data/resiliency_eval/3MW_critical_load_24hrs_outage_24hrs_recovery",
...     inputs_dir="Data/resiliency_eval/inputs_previous_stage/Paper_PGnE/Paper",
...     year=2030,
...     scenario_id=1,
... )
>>> ds.scenario_id
1
class sdom.resiliency.DesignedSystem(storage_caps: dict[str, dict[str, float]] = <factory>, thermal_caps: dict[str, dict[str, float]] = <factory>, solar_caps: dict[str, float] = <factory>, wind_caps: dict[str, float] = <factory>, load: ~pandas.core.series.Series | None = None, cf_solar: ~pandas.core.frame.DataFrame | None = None, cf_wind: ~pandas.core.frame.DataFrame | None = None, nuclear: ~pandas.core.series.Series | None = None, hydro: ~pandas.core.series.Series | None = None, other_renewables: ~pandas.core.series.Series | None = None, import_cap: ~pandas.core.series.Series | None = None, import_price: ~pandas.core.series.Series | None = None, export_cap: ~pandas.core.series.Series | None = None, export_price: ~pandas.core.series.Series | None = None, phi_fix_t: ~pandas.core.series.Series | None = None, phi_var_t: ~pandas.core.series.Series | None = None, month_of_hour: ~pandas.core.series.Series | None = None, scenario_id: int = 1, year: int = 2030, formulation_map: dict[str, str] = <factory>)[source]#

Fixed-capacity designed system loaded from SDOM output snapshots.

Parameters:
  • storage_caps (dict) – Mapping {tech: {"Cap_Pch", "Cap_Pdis", "Cap_E", "eta_ch", "eta_dis", "soc_min_frac", "vom"}} for each storage technology with non-zero capacity. Capacities are in MW / MWh.

  • thermal_caps (dict) – Mapping {tech: {"capacity_MW", "heat_rate", "fuel_cost", "vom", "var_cost"}} for each thermal technology with non-zero capacity. var_cost = heat_rate * fuel_cost + vom.

  • solar_caps (dict) – Mapping {plant_id: capacity_MW} for selected solar plants.

  • wind_caps (dict) – Mapping {plant_id: capacity_MW} for selected wind plants.

  • load (pandas.Series) – Hourly time-series (length 8760) indexed by hour-of-year (1..8760).

  • nuclear (pandas.Series) – Hourly time-series (length 8760) indexed by hour-of-year (1..8760).

  • hydro (pandas.Series) – Hourly time-series (length 8760) indexed by hour-of-year (1..8760).

  • other_renewables (pandas.Series) – Hourly time-series (length 8760) indexed by hour-of-year (1..8760).

  • cf_solar (pandas.DataFrame) – Hourly capacity factors with columns indexed by plant id.

  • cf_wind (pandas.DataFrame) – Hourly capacity factors with columns indexed by plant id.

  • import_cap (pandas.Series) – Hourly grid-exchange capacity and price series.

  • import_price (pandas.Series) – Hourly grid-exchange capacity and price series.

  • export_cap (pandas.Series) – Hourly grid-exchange capacity and price series.

  • export_price (pandas.Series) – Hourly grid-exchange capacity and price series.

  • phi_fix_t (pandas.Series) – Hourly fixed and variable demand-charge tariffs (USD/MW or USD/MWh).

  • phi_var_t (pandas.Series) – Hourly fixed and variable demand-charge tariffs (USD/MW or USD/MWh).

  • month_of_hour (pandas.Series) – Mapping from hour-of-year (1..8760) to calendar month (1..12) used to bill demand charges per month.

  • scenario_id (int) – Scenario / Run id resolved from the snapshot CSVs.

  • year (int) – Calendar year of the snapshot.

  • formulation_map (dict) – Mapping {component: formulation_name} resolved from defaults plus user-provided overrides.

storage_caps: dict[str, dict[str, float]]#
thermal_caps: dict[str, dict[str, float]]#
solar_caps: dict[str, float]#
wind_caps: dict[str, float]#
load: Series | None = None#
cf_solar: DataFrame | None = None#
cf_wind: DataFrame | None = None#
nuclear: Series | None = None#
hydro: Series | None = None#
other_renewables: Series | None = None#
import_cap: Series | None = None#
import_price: Series | None = None#
export_cap: Series | None = None#
export_price: Series | None = None#
phi_fix_t: Series | None = None#
phi_var_t: Series | None = None#
month_of_hour: Series | None = None#
scenario_id: int = 1#
year: int = 2030#
formulation_map: dict[str, str]#
__init__(storage_caps: dict[str, dict[str, float]] = <factory>, thermal_caps: dict[str, dict[str, float]] = <factory>, solar_caps: dict[str, float] = <factory>, wind_caps: dict[str, float] = <factory>, load: ~pandas.core.series.Series | None = None, cf_solar: ~pandas.core.frame.DataFrame | None = None, cf_wind: ~pandas.core.frame.DataFrame | None = None, nuclear: ~pandas.core.series.Series | None = None, hydro: ~pandas.core.series.Series | None = None, other_renewables: ~pandas.core.series.Series | None = None, import_cap: ~pandas.core.series.Series | None = None, import_price: ~pandas.core.series.Series | None = None, export_cap: ~pandas.core.series.Series | None = None, export_price: ~pandas.core.series.Series | None = None, phi_fix_t: ~pandas.core.series.Series | None = None, phi_var_t: ~pandas.core.series.Series | None = None, month_of_hour: ~pandas.core.series.Series | None = None, scenario_id: int = 1, year: int = 2030, formulation_map: dict[str, str] = <factory>) None#
class sdom.resiliency.BaselineState(soc_trajectory: ~pandas.core.frame.DataFrame | None = None, solver_status: str | None = None, objective_value: float | None = None, metadata: dict[str, ~typing.Any] = <factory>)[source]#

Placeholder container for baseline-dispatch outputs (Phase 2).

Parameters:
  • soc_trajectory (pandas.DataFrame, optional) – Hourly state-of-charge per storage technology (hour x tech).

  • solver_status (str, optional) – Solver termination status from the baseline run.

  • objective_value (float, optional) – Baseline objective value (USD).

  • metadata (dict, optional) – Free-form solver / run metadata.

soc_trajectory: DataFrame | None = None#
solver_status: str | None = None#
objective_value: float | None = None#
metadata: dict[str, Any]#
__init__(soc_trajectory: ~pandas.core.frame.DataFrame | None = None, solver_status: str | None = None, objective_value: float | None = None, metadata: dict[str, ~typing.Any] = <factory>) None#

Outage Specification#

class sdom.resiliency.OutageSpec(duration_hours: int, recovery_hours: int | dict[str, int], outaged_assets: dict[str, str | ~typing.Iterable], derating_factors: dict[tuple[str, str], float] = <factory>, min_soc_recovery: dict[str, float] | None = None, per_asset_durations: dict[tuple[str, str], int] = <factory>)[source]#

Specification of an outage / de-rating scenario.

Parameters:
  • duration_hours (int) – Outage duration applied to all listed assets unless overridden by per_asset_durations.

  • recovery_hours (int or dict) – Hours allowed after the outage for storage devices to recover the SOC target. A single int broadcasts to every storage technology; a dict {tech: hours} gives per-tech values.

  • outaged_assets (dict) – Mapping {component: asset_selector} where component is one of VALID_COMPONENTS and asset_selector is either the string "all" or an iterable of asset/plant IDs. The must-run components in MUST_RUN_COMPONENTS only accept "all" in this iteration; iterables raise NotImplementedError from validate().

  • derating_factors (dict, optional) – Mapping {(component, asset_id): factor} with each factor in [0, 1]. Missing assets default to 0 (full outage) when listed in outaged_assets, or 1 (no outage) otherwise.

  • min_soc_recovery (dict, optional) – Per-tech target SOC fraction at the end of the recovery window. None (default) -> baseline SOC at the end of the recovery window divided by Cap_E per tech.

  • per_asset_durations (dict, optional) – Optional per-asset overrides {(component, asset_id): hours} for duration_hours.

Raises:

ValueError – If a derating factor is outside [0, 1], a component name is not in VALID_COMPONENTS, or a string asset selector is not "all".

duration_hours: int#
recovery_hours: int | dict[str, int]#
outaged_assets: dict[str, str | Iterable]#
derating_factors: dict[tuple[str, str], float]#
min_soc_recovery: dict[str, float] | None = None#
per_asset_durations: dict[tuple[str, str], int]#
__post_init__() None[source]#
validate(designed_system) None[source]#

Validate the specification against designed_system.

Parameters:

designed_system (DesignedSystem) – Source of truth for asset universes (storage techs, thermal plants, VRE plants).

Raises:
  • ValueError – If an asset id in outaged_assets does not match the designed system, or if recovery_hours is a dict missing techs present in the designed system.

  • NotImplementedError – If a must-run component is given a per-asset iterable.

resolve_recovery_hours(designed_system) dict[str, int][source]#

Return {tech: recovery_hours} for every storage tech.

Raises:

ValueError – If recovery_hours is a dict that does not cover every storage tech in designed_system.

resolve_duration(component: str, asset_id: str) int[source]#

Return outage duration (hours) for a given (component, asset_id).

resolve_derating(component: str, asset_id: str) float[source]#

Return the derating multiplier rho for (component, asset_id).

Returns 1.0 when the asset is not selected for outage. Returns the user-provided factor (or the default 0.0) when it is.

resolve_outaged_asset_ids(component: str, designed_system) list[str][source]#

Return the concrete list of outaged asset ids for component.

For must-run components and imports the returned list is the canonical universe (["all"] and ["grid"] respectively).

resolve_min_soc_recovery(baseline_results, designed_system, *, recovery_end_hour: dict[str, int]) dict[str, float][source]#

Return {tech: SOC_target_fraction} at the end of recovery.

If min_soc_recovery is None, fractions default to SOC_baseline[tech, recovery_end_hour[tech]] / Cap_E[tech]. Per-tech entries in min_soc_recovery override the baseline default.

Parameters:
  • baseline_results (BaselineDispatchResults) – Provides soc_trajectory (DataFrame indexed by hour with one column per tech, MWh).

  • designed_system (DesignedSystem) – Provides storage_caps[tech]["Cap_E"].

  • recovery_end_hour (dict) – Mapping {tech: hour_of_year} of the recovery end-hour (already clipped to the baseline horizon).

sdom.resiliency.VALID_COMPONENTS#

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

sdom.resiliency.MUST_RUN_COMPONENTS#

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

Baseline Dispatch (Problem B)#

sdom.resiliency.build_baseline_dispatch(designed_system, *, n_hours=8760, min_soc_per_tech=None, curtailment_penalty=0.0, formulation_overrides=None, model_name='SDOM_BaselineDispatch', profile=False)[source]#

Build the fixed-capacity annual economic-dispatch Pyomo LP.

Composes blocks for storage, thermal (per-plant), solar (per-plant), wind (per-plant), nuclear, hydro, other-renewables, load, imports (with monthly demand charges) and exports (pure-LP capacity bound) into a single LP minimizing total operational cost.

Parameters:
  • designed_system (DesignedSystem) – Output of sdom.resiliency.load_designed_system() (or a synthetic equivalent).

  • n_hours (int, optional) – Number of hours to simulate. Default 8760.

  • min_soc_per_tech (dict, optional) – {tech: fraction in [0, 1]} enforcing operational SOC floor as a fraction of Cap_E. Missing techs default to 0.0. The user dict overrides anything stored on designed_system.storage_caps.

  • curtailment_penalty (float, optional) – Penalty applied to curtailed VRE energy (USD/MWh). Default 0.0.

  • formulation_overrides (dict, optional) – Component-wise formulation overrides (currently advisory; defaults already match DesignedSystem.formulation_map).

  • model_name (str, optional) – Pyomo model name. Default "SDOM_BaselineDispatch".

  • profile (bool, optional) – When True, instrument every block-construction step with a ModelInitProfiler, attach it as model.profiler and print a summary table at the end. Default False.

Return type:

pyomo.environ.ConcreteModel

Notes

The resiliency module is self-contained and does not call the legacy objective from formulations_system.py.

sdom.resiliency.run_baseline_dispatch(model, *, solver='highs', solver_options=None, tee=False, profile=False)[source]#

Solve the baseline dispatch and collect per-hour trajectories.

Parameters:
  • model (pyomo.environ.ConcreteModel) – Model returned by build_baseline_dispatch().

  • solver (str, optional) – Pyomo solver name. "highs" first tries appsi_highs then falls back to highs. Default "highs".

  • solver_options (dict, optional) – Extra options passed to solver.solve(..., options=...).

  • tee (bool, optional) – Stream solver output to the console. Default False.

  • profile (bool, optional) – When True, instrument the solve / extraction stages with a ModelInitProfiler and print a summary table. Default False.

Returns:

Hourly dispatch trajectories plus solver metadata.

Return type:

BaselineDispatchResults

Raises:

AttributeError – If model was not produced by build_baseline_dispatch().

class sdom.resiliency.BaselineDispatchResults(soc_trajectory: ~pandas.core.frame.DataFrame | None = None, pcha_trajectory: ~pandas.core.frame.DataFrame | None = None, pdis_trajectory: ~pandas.core.frame.DataFrame | None = None, pthermal_trajectory: ~pandas.core.frame.DataFrame | None = None, psolar_trajectory: ~pandas.core.frame.DataFrame | None = None, pwind_trajectory: ~pandas.core.frame.DataFrame | None = None, pimp: ~pandas.core.series.Series | None = None, pexp: ~pandas.core.series.Series | None = None, nuclear: ~pandas.core.series.Series | None = None, hydro: ~pandas.core.series.Series | None = None, other_renewables: ~pandas.core.series.Series | None = None, load: ~pandas.core.series.Series | None = None, month_of_hour: ~pandas.core.series.Series | None = None, objective_value: float | None = None, solver_status: str | None = None, metadata: dict[str, ~typing.Any] = <factory>)[source]#

Trajectories and metadata produced by run_baseline_dispatch().

Parameters:
  • soc_trajectory (pandas.DataFrame) – Hourly state-of-charge per storage technology, indexed by hour and with one column per tech (MWh).

  • pcha_trajectory (pandas.DataFrame) – Hourly charge / discharge per storage tech (MW).

  • pdis_trajectory (pandas.DataFrame) – Hourly charge / discharge per storage tech (MW).

  • pthermal_trajectory (pandas.DataFrame) – Hourly thermal dispatch per balancing-unit Plant_id (MW). Empty DataFrame when no thermal units survive the snapshot filter.

  • psolar_trajectory (pandas.DataFrame) – Hourly dispatched solar / wind power per plant id (MW).

  • pwind_trajectory (pandas.DataFrame) – Hourly dispatched solar / wind power per plant id (MW).

  • pimp (pandas.Series) – Hourly imports / exports (MW).

  • pexp (pandas.Series) – Hourly imports / exports (MW).

  • nuclear (pandas.Series) – Hourly time-series parameters echoed from the input system (MW).

  • hydro (pandas.Series) – Hourly time-series parameters echoed from the input system (MW).

  • other_renewables (pandas.Series) – Hourly time-series parameters echoed from the input system (MW).

  • load (pandas.Series) – Hourly time-series parameters echoed from the input system (MW).

  • month_of_hour (pandas.Series) – Hour -> month mapping used by the demand-charge billing.

  • objective_value (float) – Operational objective value (USD).

  • solver_status (str) – Solver termination condition (e.g. "optimal").

  • metadata (dict, optional) – Free-form solver / run metadata.

soc_trajectory: DataFrame | None = None#
pcha_trajectory: DataFrame | None = None#
pdis_trajectory: DataFrame | None = None#
pthermal_trajectory: DataFrame | None = None#
psolar_trajectory: DataFrame | None = None#
pwind_trajectory: DataFrame | None = None#
pimp: Series | None = None#
pexp: Series | None = None#
nuclear: Series | None = None#
hydro: Series | None = None#
other_renewables: Series | None = None#
load: Series | None = None#
month_of_hour: Series | None = None#
objective_value: float | None = None#
solver_status: str | None = None#
metadata: dict[str, Any]#
__init__(soc_trajectory: ~pandas.core.frame.DataFrame | None = None, pcha_trajectory: ~pandas.core.frame.DataFrame | None = None, pdis_trajectory: ~pandas.core.frame.DataFrame | None = None, pthermal_trajectory: ~pandas.core.frame.DataFrame | None = None, psolar_trajectory: ~pandas.core.frame.DataFrame | None = None, pwind_trajectory: ~pandas.core.frame.DataFrame | None = None, pimp: ~pandas.core.series.Series | None = None, pexp: ~pandas.core.series.Series | None = None, nuclear: ~pandas.core.series.Series | None = None, hydro: ~pandas.core.series.Series | None = None, other_renewables: ~pandas.core.series.Series | None = None, load: ~pandas.core.series.Series | None = None, month_of_hour: ~pandas.core.series.Series | None = None, objective_value: float | None = None, solver_status: str | None = None, metadata: dict[str, ~typing.Any] = <factory>) None#

Outage Dispatch (Problem O)#

sdom.resiliency.build_outage_dispatch(baseline_results, *, start_hour, outage_spec, designed_system=None, slack_penalty=10000.0, curtailment_penalty=0.0, min_soc_per_tech=None, n_hours=8760, model_name='SDOM_OutageDispatch', profile=False)[source]#

Build the per-hour outage economic-dispatch Pyomo LP anchored at start_hour.

Parameters:
  • baseline_results (BaselineDispatchResults) – Output of sdom.resiliency.run_baseline_dispatch(). Used to seed the initial SOC and (optionally) the recovery SOC target. May carry the originating DesignedSystem under baseline_results.metadata["designed_system"].

  • start_hour (int) – Anchor hour h. The outage horizon is [h, h + duration + max_recovery - 1], clipped to [1, n_hours].

  • outage_spec (OutageSpec) – Outage / de-rating specification.

  • designed_system (DesignedSystem, optional) – Source of truth for capacities and time series. If None, taken from baseline_results.metadata["designed_system"].

  • slack_penalty (float, optional) – Penalty \(\pi^{slack}\) (USD/MWh) on unserved-energy slack. Default 10_000.0.

  • curtailment_penalty (float, optional) – Penalty on curtailed VRE energy (USD/MWh). Default 0.0.

  • min_soc_per_tech (dict, optional) – Operational SOC floor per storage tech (fraction of Cap_E). Same semantics as build_baseline_dispatch().

  • n_hours (int, optional) – Length of the baseline horizon used for end-of-year clipping. Default 8760.

  • model_name (str, optional) – Pyomo model name. Default "SDOM_OutageDispatch".

  • profile (bool, optional) – When True, instrument the build stages with a ModelInitProfiler, attach it as model.profiler and print a summary table. Default False. Note: enabling this from inside a parallel ProcessPoolExecutor will produce one summary per worker on each spawned process, which is rarely useful. Prefer profile only in serial runs.

Returns:

A Pyomo LP exposing model.h (hour set), model.u (slack, NonNegativeReals), the standard dispatch sub-blocks (storage, thermal, solar, wind, imports, exports) without demand-charge variables, and an objective that minimises operational cost plus slack and curtailment penalties.

Return type:

pyomo.environ.ConcreteModel

Raises:
  • ValueError – If designed_system is not provided and cannot be recovered from baseline_results.metadata, or if validation of outage_spec against designed_system fails.

  • TypeError – If baseline_results or designed_system is the wrong type.

Notes

Initial SOC is seeded with pyomo.environ.Var.fix() rather than via an equality constraint, which keeps the LP tighter and avoids one extra row per storage technology.

Imports Formulation with Demand Charges#

The opt-in ImportsWithDemandChargesFormulation block builder used by the baseline LP. Pure linear program; not registered in io_manager.py.

sdom.resiliency.add_imports_with_demand_charges(model, *, import_cap: Series, import_price: Series, phi_fix_t: Series, phi_var_t: Series, month_of_hour: Series, block_name: str = 'imports')[source]#

Attach the ImportsWithDemandChargesFormulation block to model.

The block adds hourly import variables Pimp[t], monthly fixed and variable demand-charge variables D_fix[m] / D_var[m], the capacity bound, the demand-charge linking inequalities, and an additive cost expression total_cost_expr suitable for inclusion in any objective.

Parameters:
  • model (pyomo.environ.ConcreteModel) – Host model. Must already define model.h (a Pyomo Set of hour indices).

  • import_cap (pandas.Series) – Hourly import capacity \(\overline{P}^{imp}_{t}\) (MW), indexed by hour matching model.h.

  • import_price (pandas.Series) – Hourly import energy price \(c^{imp}_{t}\) (USD/MWh).

  • phi_fix_t (pandas.Series) – Hourly fixed demand-charge tariff \(\phi^{fix}_{t}\) (USD/MW). Must be constant within each calendar month.

  • phi_var_t (pandas.Series) – Hourly variable (time-of-use) demand-charge tariff \(\phi^{var}_{t}\) (USD/MW).

  • month_of_hour (pandas.Series) – Mapping hour -> month integer.

  • block_name (str, optional) – Name of the sub-block attached to model. Default "imports".

Returns:

The block that was attached to model.

Return type:

pyomo.environ.Block

Raises:

AttributeError – If model.h is not present.

Notes

The contributed cost is

\[Z_{imp,dc} = \sum_t c^{imp}_t\, p^{imp}_t + \sum_m \left( D^{fix}_m + D^{var}_m \right),\]

with linking constraints \(D^{k}_{m} \ge \phi^{k}_{t}\, p^{imp}_{t}\) for all \(t \in \mathcal{T}_m\) and \(k \in \{fix, var\}\).

Examples

>>> import pyomo.environ as pyo
>>> import pandas as pd
>>> m = pyo.ConcreteModel()
>>> m.h = pyo.RangeSet(1, 24)
>>> idx = range(1, 25)
>>> add_imports_with_demand_charges(  
...     m,
...     import_cap=pd.Series(100.0, index=idx),
...     import_price=pd.Series(1.0, index=idx),
...     phi_fix_t=pd.Series(50.0, index=idx),
...     phi_var_t=pd.Series(2.0, index=idx),
...     month_of_hour=pd.Series(1, index=idx),
... )

Parallel Runner#

sdom.resiliency.run_resiliency_evaluation(baseline_results, *, outage_spec, designed_system=None, hours=None, slack_penalty=10000.0, curtailment_penalty=0.0, min_soc_per_tech=None, n_hours=8760, n_workers=None, solver='highs', solver_options=None, profile_outages=False)[source]#

Run the per-hour outage evaluation in parallel and aggregate metrics.

For every anchor hour h in hours the runner builds the short-horizon outage LP via build_outage_dispatch(), solves it with HiGHS, and records EUE, USE_hours, max_unserved_MW and the objective value. Per-hour problems are independent and are solved in parallel via concurrent.futures.ProcessPoolExecutor when n_workers > 1; n_workers == 1 (or a single-hour run) falls back to an in-process serial loop.

Parameters:
  • baseline_results (BaselineDispatchResults) – Output of run_baseline_dispatch(). Must carry an SOC trajectory; if its metadata does not contain a "designed_system" entry, designed_system must be supplied.

  • outage_spec (OutageSpec) – Outage / de-rating specification, broadcast to every anchor hour.

  • designed_system (DesignedSystem, optional) – Source of truth for capacities and time series. Takes precedence over baseline_results.metadata['designed_system'] when both are provided.

  • hours (iterable of int, optional) – Anchor hours to evaluate. None (default) evaluates every hour 1..n_hours. Duplicates are removed; the result is sorted.

  • slack_penalty (float, optional) – Penalty (USD/MWh) applied to slack u[t]. Default 10_000.

  • curtailment_penalty (float, optional) – Penalty applied to curtailed VRE energy (USD/MWh). Default 0.

  • min_soc_per_tech (dict, optional) – Operational SOC floor per storage tech (fraction of Cap_E); forwarded to the outage builder.

  • n_hours (int, optional) – Length of the baseline horizon used for end-of-year clipping. Default 8760.

  • n_workers (int, optional) – Worker pool size. None -> max(1, os.cpu_count() - 1); 1 forces serial mode. Always clamped to len(hours).

  • solver (str, optional) – Pyomo solver name. "highs" first tries appsi_highs. Default "highs".

  • solver_options (dict, optional) – Solver options forwarded to solver.solve(..., options=...).

  • profile_outages (bool, optional) – When True and n_workers == 1, profile every per-hour outage build via ModelInitProfiler. Ignored (with a warning) when n_workers > 1 because each worker would emit its own summary on a separate process. Default False.

Returns:

Per-hour records (sorted by start_hour) plus run metadata including n_workers_used and a reference to outage_spec.

Return type:

ResiliencyResults

Raises:

Notes

Failure isolation: if any worker raises an exception, that hour’s record is marked solver_status="error" and its error_message field carries the formatted traceback. Other hours continue normally.

Results Container#

class sdom.resiliency.ResiliencyResults(per_hour: ~pandas.core.frame.DataFrame, metadata: dict[str, ~typing.Any] = <factory>)[source]#

Per-hour outage outcomes (lightweight, Phase 5).

Aggregate metrics (LOLP, LOLE, percentiles) and plotting are added in Phase 6.

Parameters:
  • per_hour (pandas.DataFrame) – Indexed by hour (anchor start_hour). Columns include ["EUE", "USE_hours", "max_unserved_MW", "objective_value", "solver_status", "solve_time_s", "truncated", "error_message"].

  • metadata (dict) – Free-form run metadata. Conventionally includes {"n_workers_used", "outage_spec", "n_hours", "solver"}.

per_hour: DataFrame#
metadata: dict[str, Any]#
to_dataframe() DataFrame[source]#

Return per_hour with the index promoted to a hour column.

Returns:

A copy of per_hour with hour as a regular column, sorted by hour.

Return type:

pandas.DataFrame

eue_total() float[source]#

Return the sum of per-hour expected unserved energy (MWh).

Return type:

float

metrics(*, level: str = 'aggregate')[source]#

Aggregate or per-hour resiliency metrics.

Parameters:

level ({"aggregate", "per_hour"}, optional) – "aggregate" (default) returns a dict of scalar metrics computed over the evaluated hours (errored hours excluded). "per_hour" returns a copy of per_hour with hour promoted to a column.

Return type:

dict or pandas.DataFrame

Raises:

ValueError – If level is not one of the supported values.

Notes

Aggregate metrics exclude rows with solver_status == "error"; the count of excluded rows is reported as n_errors.

lolp() float[source]#

Return the loss-of-load probability across evaluated hours.

Return type:

float

lole() float[source]#

Return the loss-of-load expectation (mean USE hours per scenario).

Return type:

float

eue(*, p: float | None = None) float[source]#

Return the mean EUE or an empirical percentile of EUE.

Parameters:

p (float, optional) – Quantile in (0, 1). Default None returns the mean EUE.

Return type:

float

Raises:

ValueError – If p is provided and not in (0, 1).

save(path: str | Path | None = None) Path[source]#

Persist per-hour records and aggregate metrics to disk.

Parameters:

path (str or pathlib.Path, optional) – Output directory. Default: ./results_resiliency/ relative to the current working directory. The directory is created if it does not exist.

Returns:

The directory the artifacts were written to.

Return type:

pathlib.Path

Raises:

ImportError – If no Parquet engine (pyarrow or fastparquet) is available.

Notes

Writes two files to path:

  • per_hour.parquet - the per-hour DataFrame.

  • summary.json - aggregate metrics + JSON-safe metadata.

classmethod load(path: str | Path) ResiliencyResults[source]#

Load a previously-saved ResiliencyResults from path.

Parameters:

path (str or pathlib.Path) – Directory that previously received save().

Return type:

ResiliencyResults

Raises:

FileNotFoundError – If per_hour.parquet or summary.json is missing.

__init__(per_hour: ~pandas.core.frame.DataFrame, metadata: dict[str, ~typing.Any] = <factory>) None#

Plotting#

sdom.resiliency.plot_metric_distribution(results: ResiliencyResults, *, metric: str = 'EUE', kind: str = 'hist', ax: Axes | None = None, **plot_kwargs: Any) Axes[source]#

Plot the empirical distribution of a per-hour metric.

Parameters:
  • results (ResiliencyResults) – Container produced by sdom.resiliency.run_resiliency_evaluation().

  • metric (str, optional) – Numeric column of results.per_hour to plot. Default "EUE".

  • kind ({"hist", "ecdf", "exceedance"}, optional) –

    Plot style. Default "hist".

    • "hist" - histogram of metric values.

    • "ecdf" - empirical CDF, monotonically non-decreasing in [0, 1].

    • "exceedance" - exceedance curve 1 - ECDF, monotonically non-increasing.

  • ax (matplotlib.axes.Axes, optional) – Existing axes to draw on. A new figure/axes is created when None.

  • **plot_kwargs – Forwarded to the underlying matplotlib call (ax.hist for kind="hist"; ax.plot otherwise).

Return type:

matplotlib.axes.Axes

Raises:
  • ImportError – If matplotlib is not importable.

  • ValueError – If kind is not a supported value or metric is not a numeric column of results.per_hour.