diff --git a/.claude/settings.local.json b/.claude/settings.local.json new file mode 100644 index 0000000000000000000000000000000000000000..3efe268a8760deccd017111ea03eec9aa0a3c7f1 --- /dev/null +++ b/.claude/settings.local.json @@ -0,0 +1,82 @@ +{ + "permissions": { + "allow": [ + "Bash(dir:*)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" log --oneline -20)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" config --list)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" status)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" diff)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" log --oneline -5)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" add CLAUDE.md)", + "Bash(git -C \"C:\\\\Users\\\\snk\\\\powersheds\" commit -m \"$\\(cat <<''EOF''\nadd CLAUDE.md for Claude Code guidance\n\nProvides build commands, architecture overview, and key concepts\nfor future Claude Code instances working in this repository.\n\nCo-Authored-By: Claude Opus 4.5 \nEOF\n\\)\")", + "Bash(.venv/Scripts/activate)", + "Bash(maturin develop:*)", + "Bash(.venv/Scripts/python:*)", + "Bash(.venv/Scripts/pip install:*)", + "Bash(.venv/Scripts/maturin develop:*)", + "Bash(uv run maturin:*)", + "Bash(uv run python:*)", + "Bash(uv run pytest:*)", + "Bash(uv add:*)", + "Bash(uv sync)", + "Bash(uv pip install:*)", + "Bash(git stash:*)", + "Bash(git checkout:*)", + "Bash(git add:*)", + "Bash(git commit:*)", + "Bash(cargo clean:*)", + "Bash(ls:*)", + "Bash(cargo build:*)", + "Bash(git -C C:/Users/snk/powersheds status)", + "Bash(git -C C:/Users/snk/powersheds log --oneline -5)", + "Bash(git -C C:/Users/snk/powersheds diff src/helpers.rs)", + "Bash(git -C C:/Users/snk/powersheds diff src/lib.rs)", + "Bash(git -C C:/Users/snk/powersheds add src/helpers.rs src/lib.rs tests/test_hpf.py)", + "Bash(git -C C:/Users/snk/powersheds commit -m \"$\\(cat <<''EOF''\nHandle infeasible power targets with symmetric clamping\n\nWhen target power falls outside the feasible range at a given head,\nthe HPF interpolation now clamps to the nearest feasible boundary\n\\(min or max\\) rather than returning NaN. This maintains roundtrip\nconsistency between P→Q and Q→P functions.\n\nChanges:\n- src/helpers.rs: Implement symmetric clamping for infeasible power,\n use linear scan for P index lookup \\(binary search unreliable with\n approx_eq\\), remove duplicate actual_power calculation in lib.rs\n- tests/test_hpf.py: Update tolerances for flat regions \\(~2 MW corner,\n ~1.5 MW edge\\), filter clamping cases from roundtrip tests, increase\n performance thresholds for linear scan algorithm\n\nFixes NaN values in target_release for Cumberland example simulation.\n\nCo-Authored-By: Claude Opus 4.5 \nEOF\n\\)\")", + "Bash(git -C C:/Users/snk/powersheds push)", + "Bash(python:*)", + "Bash(python3:*)", + "Bash(find:*)", + "Bash(.venvScriptsactivate.ps1)", + "Bash(.\\\\.venv\\\\Scripts\\\\python.exe:*)", + "Bash(.venv/Scripts/python.exe:*)", + "Bash(git rm:*)", + "Bash(git push:*)", + "Bash(.venv\\\\Scripts\\\\python.exe:*)", + "Bash(C:Userssnkpowersheds.venvScriptspython.exe examples/plot_network.py examples/Cumberland/cascade_config.yaml)", + "Bash(cmd /c \"C:\\\\Users\\\\snk\\\\powersheds\\\\.venv\\\\Scripts\\\\activate.bat && maturin develop\")", + "Bash(cmd /c \"C:\\\\Users\\\\snk\\\\powersheds\\\\.venv\\\\Scripts\\\\activate.bat && maturin develop 2>&1\")", + "Bash(C:Userssnkpowersheds.venvScriptsmaturin.exe develop)", + "Bash(powershell -Command:*)", + "WebSearch", + "Bash(..venvScriptsactivate.ps1)", + "Bash(source:*)", + "Bash(jupyter nbconvert:*)", + "Bash(head:*)", + "Bash(. .venv/Scripts/activate)", + "Bash(done)", + "Bash(grep:*)", + "Bash(\"../../../.venv/Scripts/python.exe\" -c:*)", + "Bash(.venv/Scripts/maturin.exe:*)", + "Bash(tasklist:*)", + "Bash(taskkill:*)", + "Bash(cmd //c \"taskkill /F /PID 48996\")", + "Bash(cmd //c:*)", + "Bash(wc:*)", + "Bash(for f in /c/Users/snk/powersheds/examples/private/Kennebec/storage_HWelevation_tables/*.csv)", + "Bash(do)", + "Bash(basename=\"$f##*/\")", + "Bash(basename=\"$basename%.csv\")", + "Bash(if [ ! -f \"/c/Users/snk/powersheds/examples/private/Kennebec/outflow_TWelevation_tables/$basename.csv\" ])", + "Bash(then)", + "Bash(echo:*)", + "Bash(fi:*)", + "Bash(md5sum:*)", + "Bash(xargs:*)", + "Bash(\"C:/Users/snk/powersheds/.venv/Scripts/python\" -m jupyter nbconvert --to notebook --execute experiments.ipynb --output executed_experiments.ipynb --ExecutePreprocessor.timeout=600)", + "Bash(\"C:/Users/snk/powersheds/.venv/Scripts/python\" \"C:/Users/snk/powersheds/examples/private/Kennebec/_test_exp.py\")", + "Bash(\"C:/Users/snk/powersheds/.venv/Scripts/python\" -c \"import powersheds; print\\(powersheds.__file__\\)\")", + "Bash(Rscript:*)" + ] + } +} diff --git a/README.md b/README.md index b19b1d6cd78b1fda12bb7c81da3a14efda81c3cc..612d19205b4b30f5cfc774f81324f8719d86769e 100644 --- a/README.md +++ b/README.md @@ -81,46 +81,74 @@ Dependencies (pandas, numpy, scipy, matplotlib, PyYAML, pyarrow, etc.) are decla ```python import powersheds - -# Minimal single-reservoir cascade (toy numbers for illustration) -cascade_data = { - "reservoirs": {}, - "rivers": {}, - "confluences": {} -} +from dataclasses import dataclass, field + +@dataclass +class ReservoirData: + object_type: str + simulation_order: int + downstream_object: str + capacity: float + initial_pool_elevation: float + min_power_pool: list[float] + tailwater_elevation: float + max_release: list[float] + min_release: list[float] + set_storage: list[float] + set_elevation: list[float] + catchment_inflow: list[float] + target_power: list[float] + hpf_h: list[float] + hpf_p: list[float] + hpf_q: list[float] + set_tw_outflow: list[float] = field(default_factory=list) + set_tw_elevation: list[float] = field(default_factory=list) + max_operating_elevation: list[float] = field(default_factory=list) + +@dataclass +class CascadeData: + reservoirs: dict[str, ReservoirData] + rivers: dict[str, object] + confluences: dict[str, object] # Define a reservoir with 24 hours of inflow and power targets -cascade_data["reservoirs"]["DemoReservoir"] = { - "object_type": "reservoir", - "simulation_order": 1, - "downstream_object": "NA", # terminal object - "capacity": 500.0, # Mm3 - "initial_pool_elevation": 200.0, - "min_power_pool": 190.0, - "tailwater_elevation": 150.0, - "max_release": 2.0, # Mm3/hr - "min_release": 0.0, # Mm3/hr - # storage–elevation (SET) curve - "set_storage": [0.0, 250.0, 500.0], - "set_elevation": [180.0, 195.0, 210.0], - # hourly inputs - "catchment_inflow": [0.2] * 24, - "target_power": [50.0] * 24, # MW - # head–power–flow (HPF) grid: cartesian product of H and P - # H (m): [40, 60] - # P (MW): [ 0, 50, 100] - # Q (m3/s): row-major values for each (H,P) pair - "hpf_h": [40, 40, 40, 60, 60, 60], - "hpf_p": [ 0, 50,100, 0, 50,100], - "hpf_q": [ 0, 60,120, 0, 45, 90], -} +cascade_data = CascadeData( + reservoirs={ + "DemoReservoir": ReservoirData( + object_type="reservoir", + simulation_order=1, + downstream_object="NA", # terminal object + capacity=500.0, # Mm3 + initial_pool_elevation=200.0, + min_power_pool=[190.0] * 24, + tailwater_elevation=150.0, + max_release=[2.0] * 24, # Mm3/hr + min_release=[0.0] * 24, # Mm3/hr + # storage–elevation (SET) curve + set_storage=[0.0, 250.0, 500.0], + set_elevation=[180.0, 195.0, 210.0], + # hourly inputs + catchment_inflow=[0.2] * 24, + target_power=[50.0] * 24, # MW + # head–power–flow (HPF) grid: cartesian product of H and P + # H (m): [40, 60] + # P (MW): [ 0, 50, 100] + # Q (m3/s): row-major values for each (H,P) pair + hpf_h=[40, 40, 40, 60, 60, 60], + hpf_p=[0, 50, 100, 0, 50, 100], + hpf_q=[0, 60, 120, 0, 45, 90], + ) + }, + rivers={}, + confluences={}, +) results = powersheds.simulate_cascade(cascade_data) res = results["DemoReservoir"] -print("release (Mm3/hr)", res.release[:5]) -print("storage (Mm3)", res.storage[:5]) -print("power (MW)", res.actual_power[:5]) +print("release (Mm3/hr)", res["release"][:5]) +print("storage (Mm3)", res["storage"][:5]) +print("power (MW)", res["actual_power"][:5]) ``` #### Notes @@ -133,7 +161,7 @@ print("power (MW)", res.actual_power[:5]) ## 📚 Data Model -Pass a single `cascade_data` mapping with three sections: `reservoirs`, `rivers`, and `confluences`. Each contains a dictionary of named objects. +Pass a single `cascade_data` object with `reservoirs`, `rivers`, and `confluences` attributes. A lightweight Python `dataclass` works well; each section should be a dictionary of named objects with matching attributes. Reservoir @@ -141,9 +169,12 @@ Reservoir - `simulation_order`: integer; lower runs earlier - `downstream_object`: name of downstream object or "NA" - `capacity`: reservoir storage capacity (Mm3) -- `initial_pool_elevation` and `min_power_pool` (m) +- `initial_pool_elevation` (m) +- `min_power_pool`: list[float] length T (m) - `tailwater_elevation` (m) -- `max_release` and `min_release` (Mm3/hr) +- `max_release`: list[float] length T (Mm3/hr) +- `min_release`: list[float] length T (Mm3/hr) +- `max_operating_elevation`: optional list[float] length T, or empty - `set_storage` and `set_elevation`: arrays for the storage–elevation curve - `catchment_inflow`: list[float] length T - `target_power`: list[float] length T (MW) @@ -170,12 +201,12 @@ Confluence Python entry point ```python -powersheds.simulate_cascade(cascade_data: dict) -> dict[str, Result] +powersheds.simulate_cascade(cascade_data) -> dict[str, Result] ``` Returned per-object results: -- Reservoir: `storage`, `pool_elevation`, `head`, `release`, `spill`, `catchment_inflow`, `total_inflow`, `target_release`, `target_power`, `actual_power` +- Reservoir: `storage`, `pool_elevation`, `head`, `release`, `spill`, `catchment_inflow`, `total_inflow`, `target_release`, `target_power`, `actual_power`, `tailwater_elevation`, `violation_min_power_pool`, `violation_min_release`, `violation_max_release`, `violation_max_operating_elevation` - River: `inflow`, `outflow` (with travel-time lag) - Confluence: `inflow`, `outflow` (no lag) diff --git a/examples/Cumberland/cumberland_diagnostics.html b/examples/Cumberland/cumberland_diagnostics.html new file mode 100644 index 0000000000000000000000000000000000000000..31b34a631a720cbead424d55526ecec6e7b722bb --- /dev/null +++ b/examples/Cumberland/cumberland_diagnostics.html @@ -0,0 +1,3888 @@ + + + +
+
+ + \ No newline at end of file diff --git a/examples/Cumberland/notebook.ipynb b/examples/Cumberland/notebook.ipynb index 330a29cd9717490cbf731ac23a2061516e6cdb68..5982aed037014350cd105ba498920e31a410d093 100644 --- a/examples/Cumberland/notebook.ipynb +++ b/examples/Cumberland/notebook.ipynb @@ -24,14 +24,14 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "from dataclasses import dataclass, field\n\n@dataclass\nclass ReservoirData:\n object_type: str\n capacity: float\n initial_pool_elevation: float\n min_power_pool: float\n set_storage: list[float]\n set_elevation: list[float]\n hpf_h: list[float]\n hpf_p: list[float]\n hpf_q: list[float]\n tailwater_elevation: float\n max_release: float\n min_release: float\n catchment_inflow: list[float]\n target_power: list[float]\n simulation_order: int\n downstream_object: str\n set_tw_outflow: list[float] = field(default_factory=list)\n set_tw_elevation: list[float] = field(default_factory=list)\n\n@dataclass\nclass RiverData:\n object_type: str\n simulation_order: int\n downstream_object: str\n lag: int\n legacy_flows: list[float]\n\n@dataclass\nclass ConfluenceData:\n object_type: str\n simulation_order: int\n downstream_object: str\n\n@dataclass\nclass CascadeData:\n reservoirs: dict[str, ReservoirData]\n rivers: dict[str, RiverData]\n confluences: dict[str, ConfluenceData]" + "source": "from dataclasses import dataclass, field\n\n@dataclass\nclass ReservoirData:\n object_type: str\n capacity: float\n initial_pool_elevation: float\n min_power_pool: list[float]\n set_storage: list[float]\n set_elevation: list[float]\n hpf_h: list[float]\n hpf_p: list[float]\n hpf_q: list[float]\n tailwater_elevation: float\n max_release: list[float]\n min_release: list[float]\n catchment_inflow: list[float]\n target_power: list[float]\n simulation_order: int\n downstream_object: str\n set_tw_outflow: list[float] = field(default_factory=list)\n set_tw_elevation: list[float] = field(default_factory=list)\n max_operating_elevation: list[float] = field(default_factory=list)\n\n@dataclass\nclass RiverData:\n object_type: str\n simulation_order: int\n downstream_object: str\n lag: int\n legacy_flows: list[float]\n\n@dataclass\nclass ConfluenceData:\n object_type: str\n simulation_order: int\n downstream_object: str\n\n@dataclass\nclass CascadeData:\n reservoirs: dict[str, ReservoirData]\n rivers: dict[str, RiverData]\n confluences: dict[str, ConfluenceData]" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "CUMBERLAND_DIR = Path(\".\").resolve()\n\nwith open(CUMBERLAND_DIR / \"cascade_config.yaml\") as f:\n config_dict = yaml.safe_load(f)\n\nreservoir_dict = {}\nriver_dict = {}\nconfluence_dict = {}\n\nfor name, specs in config_dict.items():\n if specs[\"object_type\"] == \"reservoir\":\n ts = pd.read_csv(CUMBERLAND_DIR / \"time_series\" / f\"{name}.csv\")\n se = pd.read_csv(CUMBERLAND_DIR / \"storage_HWelevation_tables\" / f\"{name}.csv\")\n hpf = pd.read_parquet(CUMBERLAND_DIR / \"head_power_flow_tables\" / name)\n tw_path = CUMBERLAND_DIR / \"outflow_TWelevation_tables\" / f\"{name}.csv\"\n if tw_path.exists():\n tw = pd.read_csv(tw_path)\n set_tw_outflow = tw[\"outflow_cumecs\"].tolist()\n set_tw_elevation = tw[\"elevation_m\"].tolist()\n else:\n set_tw_outflow, set_tw_elevation = [], []\n reservoir_dict[name] = ReservoirData(\n **specs,\n catchment_inflow=ts[\"catchment_inflow\"].tolist(),\n target_power=ts[\"target_power\"].tolist(),\n set_storage=se[\"storage_Mm3\"].tolist(),\n set_elevation=se[\"elevation_m\"].tolist(),\n hpf_h=hpf[\"H_m\"].astype(float).tolist(),\n hpf_p=hpf[\"P_MW\"].astype(float).tolist(),\n hpf_q=hpf[\"Q_cumecs\"].astype(float).tolist(),\n set_tw_outflow=set_tw_outflow,\n set_tw_elevation=set_tw_elevation,\n )\n elif specs[\"object_type\"] == \"river\":\n legacy = pd.read_csv(CUMBERLAND_DIR / \"time_series\" / f\"{name}.csv\")\n river_dict[name] = RiverData(\n **specs,\n legacy_flows=legacy[\"legacy_flow\"].tolist(),\n )\n elif specs[\"object_type\"] == \"confluence\":\n confluence_dict[name] = ConfluenceData(**specs)\n\ncascade_data = CascadeData(\n reservoirs=reservoir_dict,\n rivers=river_dict,\n confluences=confluence_dict,\n)\n\nprint(f\"Loaded {len(reservoir_dict)} reservoirs, {len(river_dict)} rivers, {len(confluence_dict)} confluences\")" + "source": "CUMBERLAND_DIR = Path(\".\").resolve()\n\nwith open(CUMBERLAND_DIR / \"cascade_config.yaml\") as f:\n config_dict = yaml.safe_load(f)\n\nreservoir_dict = {}\nriver_dict = {}\nconfluence_dict = {}\n\nfor name, specs in config_dict.items():\n if specs[\"object_type\"] == \"reservoir\":\n ts = pd.read_csv(CUMBERLAND_DIR / \"time_series\" / f\"{name}.csv\")\n se = pd.read_csv(CUMBERLAND_DIR / \"storage_HWelevation_tables\" / f\"{name}.csv\")\n hpf = pd.read_parquet(CUMBERLAND_DIR / \"head_power_flow_tables\" / name)\n tw_path = CUMBERLAND_DIR / \"outflow_TWelevation_tables\" / f\"{name}.csv\"\n if tw_path.exists():\n tw = pd.read_csv(tw_path)\n set_tw_outflow = tw[\"outflow_cumecs\"].tolist()\n set_tw_elevation = tw[\"elevation_m\"].tolist()\n else:\n set_tw_outflow, set_tw_elevation = [], []\n\n n_hours = len(ts)\n # Convert scalar YAML constraint values to time series\n specs = dict(specs) # copy to avoid mutating config_dict\n specs[\"min_power_pool\"] = [specs[\"min_power_pool\"]] * n_hours\n specs[\"max_release\"] = [specs[\"max_release\"]] * n_hours\n specs[\"min_release\"] = [specs[\"min_release\"]] * n_hours\n\n reservoir_dict[name] = ReservoirData(\n **specs,\n catchment_inflow=ts[\"catchment_inflow\"].tolist(),\n target_power=ts[\"target_power\"].tolist(),\n set_storage=se[\"storage_Mm3\"].tolist(),\n set_elevation=se[\"elevation_m\"].tolist(),\n hpf_h=hpf[\"H_m\"].astype(float).tolist(),\n hpf_p=hpf[\"P_MW\"].astype(float).tolist(),\n hpf_q=hpf[\"Q_cumecs\"].astype(float).tolist(),\n set_tw_outflow=set_tw_outflow,\n set_tw_elevation=set_tw_elevation,\n )\n elif specs[\"object_type\"] == \"river\":\n legacy = pd.read_csv(CUMBERLAND_DIR / \"time_series\" / f\"{name}.csv\")\n river_dict[name] = RiverData(\n **specs,\n legacy_flows=legacy[\"legacy_flow\"].tolist(),\n )\n elif specs[\"object_type\"] == \"confluence\":\n confluence_dict[name] = ConfluenceData(**specs)\n\ncascade_data = CascadeData(\n reservoirs=reservoir_dict,\n rivers=river_dict,\n confluences=confluence_dict,\n)\n\nprint(f\"Loaded {len(reservoir_dict)} reservoirs, {len(river_dict)} rivers, {len(confluence_dict)} confluences\")" }, { "cell_type": "markdown", diff --git a/python/powersheds/diagnostics.py b/python/powersheds/diagnostics.py index bd151711ce61b870e783d355c81e034533ab6005..90bd93625bd706da6be7f842daf61cdd036f0b50 100644 --- a/python/powersheds/diagnostics.py +++ b/python/powersheds/diagnostics.py @@ -55,6 +55,8 @@ _COLORS = { "infeasible_fill": "rgba(238,102,119,0.25)", "release_band": "rgba(160,160,160,0.15)", "release_bound": "rgba(140,140,140,0.4)", + "elevation_band": "rgba(102,204,238,0.18)", + "elevation_bound": "rgba(102,204,238,0.45)", } # Style constants for the matplotlib network diagram @@ -552,16 +554,19 @@ def _build_flow_traces(data, x, visible, release_constraints=None): )) if release_constraints is not None: min_rel, max_rel = release_constraints + # Normalize to lists (support both scalar and list constraints) + min_rel_y = min_rel if isinstance(min_rel, list) else [min_rel] * len(x) + max_rel_y = max_rel if isinstance(max_rel, list) else [max_rel] * len(x) # Upper bound (invisible line for fill anchor) traces.append(go.Scatter( - x=x, y=[max_rel] * len(x), + x=x, y=max_rel_y, line=dict(color=_COLORS["release_bound"], width=1, dash="dot"), mode="lines", showlegend=False, hoverinfo="skip", visible=visible, legendgroup="flows", )) # Lower bound with fill up to max_release traces.append(go.Scatter( - x=x, y=[min_rel] * len(x), + x=x, y=min_rel_y, fill="tonexty", fillcolor=_COLORS["release_band"], line=dict(color=_COLORS["release_bound"], width=1, dash="dot"), @@ -572,26 +577,29 @@ def _build_flow_traces(data, x, visible, release_constraints=None): return traces -def _build_elevation_traces(data, x, visible, show_storage): +def _build_elevation_traces( + data, + x, + visible, + show_storage, + elevation_constraints=None, +): """Build elevation panel traces for a single reservoir.""" traces = [] common = dict(visible=visible, legendgroup="elev") + if elevation_constraints is not None: + traces.extend( + _build_elevation_constraint_traces( + x, visible, elevation_constraints + ) + ) + traces.append(go.Scatter( x=x, y=data["pool_elevation"], name="Pool elevation", line=dict(color=_COLORS["pool_elevation"], width=1.5), showlegend=visible, **common, )) - traces.append(go.Scatter( - x=x, y=data["tailwater_elevation"], name="Tailwater elevation", - line=dict(color=_COLORS["tailwater_elevation"], width=1.5), - showlegend=visible, **common, - )) - traces.append(go.Scatter( - x=x, y=data["head"], name="Head", - line=dict(color=_COLORS["head"], width=1.5), - showlegend=visible, **common, - )) if show_storage: traces.append(go.Scatter( x=x, y=data["storage"], name="Storage", @@ -601,6 +609,31 @@ def _build_elevation_traces(data, x, visible, show_storage): return traces +def _build_elevation_constraint_traces(x, visible, elevation_constraints): + """Build elevation operating-band traces for a single reservoir.""" + min_elev, max_elev = elevation_constraints + min_y = min_elev if isinstance(min_elev, list) else [min_elev] * len(x) + max_y = max_elev if isinstance(max_elev, list) else [max_elev] * len(x) + + return [ + go.Scatter( + x=x, y=max_y, + line=dict(color=_COLORS["elevation_bound"], width=1, dash="dot"), + mode="lines", showlegend=False, hoverinfo="skip", + visible=visible, legendgroup="elev", + ), + go.Scatter( + x=x, y=min_y, + fill="tonexty", + fillcolor=_COLORS["elevation_band"], + line=dict(color=_COLORS["elevation_bound"], width=1, dash="dot"), + mode="lines", name="Operating band", + showlegend=visible, + visible=visible, legendgroup="elev", + ), + ] + + def _build_power_traces(data, x, visible): """Build power panel traces with infeasible generation shading.""" traces = [] @@ -621,6 +654,7 @@ def _build_power_traces(data, x, visible): target = data["target_power"] actual = data["actual_power"] fill_min = [min(a, t) for a, t in zip(actual, target)] + shortfall = [max(t - a, 0.0) for a, t in zip(actual, target)] # Upper boundary (target) — invisible line traces.append(go.Scatter( @@ -636,6 +670,13 @@ def _build_power_traces(data, x, visible): fillcolor=_COLORS["infeasible_fill"], line=dict(width=0), mode="lines", name="Infeasible", + customdata=np.column_stack([target, actual, shortfall]), + hovertemplate=( + "Target: %{customdata[0]:.3f} MW
" + "Actual: %{customdata[1]:.3f} MW
" + "Infeasible: %{customdata[2]:.3f} MW" + "" + ), showlegend=visible, visible=visible, legendgroup="power", )) @@ -655,6 +696,8 @@ def cascade_diagnostics( height: int = 900, width: int = 1100, show_storage: bool = False, + start_date: str | None = None, + window_days: int | None = None, ) -> go.Figure: """Multi-panel interactive time series diagnostic plot. @@ -672,7 +715,8 @@ def cascade_diagnostics( reservoirs found in *results*. cascade_data : CascadeData, optional The cascade configuration. When provided, min/max release - constraint bands are drawn on the Flows panel. + constraint bands are drawn on the Flows panel and an elevation + operating band is drawn on the Elevation panel. time_index : list or DatetimeIndex, optional X-axis values. Defaults to integer hour indices. title : str @@ -681,6 +725,12 @@ def cascade_diagnostics( Figure dimensions in pixels. show_storage : bool Include a storage trace on the elevation panel. + start_date : str, optional + ISO-format date string (e.g. ``"2021-04-01"``) for the start of a + focus window. Requires *time_index* to be provided. + window_days : int, optional + Length of the focus window in days. Defaults to 14 if *start_date* + is given but *window_days* is omitted. """ res_names = reservoir_names or _detect_reservoirs(results) if not res_names: @@ -689,17 +739,66 @@ def cascade_diagnostics( n_hours = len(results[res_names[0]]["actual_power"]) x = list(time_index) if time_index is not None else list(range(n_hours)) - # Extract release constraints from cascade_data when available + # --- Optional date-window focus --- + if start_date is not None: + if time_index is None: + raise ValueError( + "time_index is required when start_date is specified." + ) + import pandas as pd + from datetime import timedelta + + ti = pd.DatetimeIndex(time_index) + start_ts = pd.Timestamp(start_date) + if ti.tz is not None and start_ts.tz is None: + start_ts = start_ts.tz_localize(ti.tz) + days = window_days if window_days is not None else 14 + end_ts = start_ts + timedelta(days=days) + mask = (ti >= start_ts) & (ti < end_ts) + idx = np.where(mask)[0] + if len(idx) == 0: + raise ValueError( + f"No data in window [{start_ts}, {end_ts}). " + f"time_index spans {ti[0]} to {ti[-1]}." + ) + i_start, i_end = int(idx[0]), int(idx[-1]) + 1 + x = x[i_start:i_end] + + # Deep-slice every list value in results for each reservoir + sliced_results = {} + for name in res_names: + sliced_results[name] = { + k: v[i_start:i_end] if isinstance(v, list) else v + for k, v in results[name].items() + } + results = sliced_results + + # Extract flow/elevation constraints from cascade_data when available has_constraints = cascade_data is not None release_constraints = {} + elevation_constraints = {} if has_constraints: for name in res_names: res = cascade_data.reservoirs.get(name) if res is not None: - release_constraints[name] = ( - getattr(res, "min_release", 0.0), - getattr(res, "max_release", float("inf")), - ) + min_rel = getattr(res, "min_release", [0.0]) + max_rel = getattr(res, "max_release", [float("inf")]) + min_elev = getattr(res, "min_power_pool", None) + max_elev = getattr(res, "max_operating_elevation", []) + if not max_elev: + max_elev = [max(res.set_elevation)] * len(x) + if start_date is not None: + if isinstance(min_rel, list): + min_rel = min_rel[i_start:i_end] + if isinstance(max_rel, list): + max_rel = max_rel[i_start:i_end] + if isinstance(min_elev, list): + min_elev = min_elev[i_start:i_end] + if isinstance(max_elev, list): + max_elev = max_elev[i_start:i_end] + release_constraints[name] = (min_rel, max_rel) + if min_elev is not None: + elevation_constraints[name] = (min_elev, max_elev) # ----- Create subplots ----- fig = make_subplots( @@ -712,7 +811,7 @@ def cascade_diagnostics( # ----- Trace counts ----- n_flow = 6 if has_constraints else 4 - n_elev = 3 + int(show_storage) + n_elev = (3 if has_constraints else 1) + int(show_storage) n_power = 4 n_per_res = n_flow + n_elev + n_power @@ -730,13 +829,16 @@ def cascade_diagnostics( ) if name in release_constraints: min_r, max_r = release_constraints[name] - fr = _merged_range(fr, [min_r, max_r]) + lo = min(min_r) if isinstance(min_r, list) else min_r + hi = max(max_r) if isinstance(max_r, list) else max_r + fr = _merged_range(fr, [lo, hi]) flow_ranges[name] = fr - er = _merged_range( - _compute_axis_ranges(d, "pool_elevation"), - _compute_axis_ranges(d, "tailwater_elevation"), - _compute_axis_ranges(d, "head"), - ) + er = _compute_axis_ranges(d, "pool_elevation") + if name in elevation_constraints: + min_elev, max_elev = elevation_constraints[name] + lo = min(min_elev) if isinstance(min_elev, list) else min_elev + hi = max(max_elev) if isinstance(max_elev, list) else max_elev + er = _merged_range(er, [lo, hi]) if show_storage: er = _merged_range(er, _compute_axis_ranges(d, "storage")) elev_ranges[name] = er @@ -750,9 +852,12 @@ def cascade_diagnostics( vis = i == 0 d = results[name] rc = release_constraints.get(name) if has_constraints else None + ec = elevation_constraints.get(name) if has_constraints else None for t in _build_flow_traces(d, x, vis, release_constraints=rc): fig.add_trace(t, row=1, col=1) - for t in _build_elevation_traces(d, x, vis, show_storage): + for t in _build_elevation_traces( + d, x, vis, show_storage, elevation_constraints=ec + ): fig.add_trace(t, row=2, col=1) for t in _build_power_traces(d, x, vis): fig.add_trace(t, row=3, col=1) diff --git a/src/lib.rs b/src/lib.rs index b4ea99a1a1f17a83882d13deb75a58d2606f2762..191566b7590904e63dd3b792e8e07fdbcc8cd330 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -39,12 +39,12 @@ struct ReservoirData { object_type: String, capacity: f64, initial_pool_elevation: f64, - min_power_pool: f64, + min_power_pool: Vec, set_storage: Vec, set_elevation: Vec, tailwater_elevation: f64, - max_release: f64, - min_release: f64, + max_release: Vec, + min_release: Vec, catchment_inflow: Vec, target_power: Vec, simulation_order: i32, @@ -54,6 +54,7 @@ struct ReservoirData { hpf_q: Vec, set_tw_outflow: Vec, set_tw_elevation: Vec, + max_operating_elevation: Vec, } // Set River data struct @@ -91,6 +92,10 @@ struct ReservoirState { spill: f64, actual_power: f64, tailwater_elev: f64, + violation_min_power_pool: f64, + violation_min_release: f64, + violation_max_release: f64, + violation_max_operating_elevation: f64, } // Define results structures for each object type @@ -107,6 +112,10 @@ struct ReservoirResults { target_power: Vec, actual_power: Vec, tailwater_elevation: Vec, + violation_min_power_pool: Vec, + violation_min_release: Vec, + violation_max_release: Vec, + violation_max_operating_elevation: Vec, } #[derive(FromPyObject, IntoPyObject)] @@ -137,8 +146,19 @@ fn simulate_timestep( target_power: f64, previous_storage: f64, previous_total_outflow: f64, + t: usize, ) -> ReservoirState { + // Sanitize constraints: NaN → "no constraint" sentinel + let min_power_pool = if reservoir.min_power_pool[t].is_nan() { f64::NEG_INFINITY } else { reservoir.min_power_pool[t] }; + let min_release = if reservoir.min_release[t].is_nan() { 0.0 } else { reservoir.min_release[t] }; + let max_release = if reservoir.max_release[t].is_nan() { f64::INFINITY } else { reservoir.max_release[t] }; + let max_op_elev = if !reservoir.max_operating_elevation.is_empty() && !reservoir.max_operating_elevation[t].is_nan() { + Some(reservoir.max_operating_elevation[t]) + } else { + None + }; + // Use previous (initial) storage to estimate pool elevation let previous_pool_elevation = interpolate_elevation( previous_storage, @@ -147,11 +167,23 @@ fn simulate_timestep( ); // Determine whether the pool has dropped below the power pool minimum - let below_power_pool = previous_pool_elevation < reservoir.min_power_pool; + let below_power_pool = previous_pool_elevation < min_power_pool; // Compute available water given inflow and previous storage let available_water = previous_storage + inflow; + // Determine effective max storage: min(capacity, max_operating_storage) + let effective_capacity = if let Some(moe) = max_op_elev { + let max_op_storage = interpolate_storage( + moe, + &reservoir.set_storage, + &reservoir.set_elevation, + ); + reservoir.capacity.min(max_op_storage) + } else { + reservoir.capacity + }; + // Branch: dynamic tailwater (non-empty TW table) vs static tailwater let has_dynamic_tw = !reservoir.set_tw_outflow.is_empty(); @@ -191,14 +223,14 @@ fn simulate_timestep( 0.0 } else { target_release - .max(reservoir.min_release) - .min(reservoir.max_release) + .max(min_release) + .min(max_release) .min(available_water) }; let potential_storage = available_water - release; - spill = if potential_storage > reservoir.capacity { - potential_storage - reservoir.capacity + spill = if potential_storage > effective_capacity { + potential_storage - effective_capacity } else { 0.0 }; @@ -236,8 +268,8 @@ fn simulate_timestep( 0.0 } else { target_release - .max(reservoir.min_release) - .min(reservoir.max_release) + .max(min_release) + .min(max_release) .min(available_water) }; @@ -255,8 +287,8 @@ fn simulate_timestep( }; let potential_storage = available_water - release; - let (storage, spill) = if potential_storage > reservoir.capacity { - (reservoir.capacity, potential_storage - reservoir.capacity) + let (storage, spill) = if potential_storage > effective_capacity { + (effective_capacity, potential_storage - effective_capacity) } else { (potential_storage, 0.0) }; @@ -267,6 +299,16 @@ fn simulate_timestep( &reservoir.set_elevation, ); + // Compute violations (NaN constraints → 0.0 violation) + let violation_min_power_pool = if reservoir.min_power_pool[t].is_nan() { 0.0 } else { (min_power_pool - pool_elevation).max(0.0) }; + let violation_min_release = if reservoir.min_release[t].is_nan() { 0.0 } else { (min_release - release).max(0.0) }; + let violation_max_release = if reservoir.max_release[t].is_nan() { 0.0 } else { (release - max_release).max(0.0) }; + let violation_max_operating_elevation = if let Some(moe) = max_op_elev { + (pool_elevation - moe).max(0.0) + } else { + 0.0 + }; + ReservoirState { storage, pool_elevation, @@ -276,6 +318,10 @@ fn simulate_timestep( spill, actual_power, tailwater_elev: tw, + violation_min_power_pool, + violation_min_release, + violation_max_release, + violation_max_operating_elevation, } } else { // --- Legacy path: static tailwater --- @@ -295,8 +341,8 @@ fn simulate_timestep( 0.0 } else { target_release - .max(reservoir.min_release) - .min(reservoir.max_release) + .max(min_release) + .min(max_release) .min(available_water) }; @@ -314,8 +360,8 @@ fn simulate_timestep( }; let potential_storage = available_water - release; - let (storage, spill) = if potential_storage > reservoir.capacity { - (reservoir.capacity, potential_storage - reservoir.capacity) + let (storage, spill) = if potential_storage > effective_capacity { + (effective_capacity, potential_storage - effective_capacity) } else { (potential_storage, 0.0) }; @@ -326,6 +372,16 @@ fn simulate_timestep( &reservoir.set_elevation, ); + // Compute violations (NaN constraints → 0.0 violation) + let violation_min_power_pool = if reservoir.min_power_pool[t].is_nan() { 0.0 } else { (min_power_pool - pool_elevation).max(0.0) }; + let violation_min_release = if reservoir.min_release[t].is_nan() { 0.0 } else { (min_release - release).max(0.0) }; + let violation_max_release = if reservoir.max_release[t].is_nan() { 0.0 } else { (release - max_release).max(0.0) }; + let violation_max_operating_elevation = if let Some(moe) = max_op_elev { + (pool_elevation - moe).max(0.0) + } else { + 0.0 + }; + ReservoirState { storage, pool_elevation, @@ -335,6 +391,10 @@ fn simulate_timestep( spill, actual_power, tailwater_elev: reservoir.tailwater_elevation, + violation_min_power_pool, + violation_min_release, + violation_max_release, + violation_max_operating_elevation, } } } @@ -349,6 +409,25 @@ fn simulate_cascade( .ok_or_else(|| PyValueError::new_err("No reservoirs provided"))?; let n = first_reservoir.catchment_inflow.len(); + // Validate constraint time series lengths + for (name, reservoir) in &cascade_data.reservoirs { + if reservoir.min_power_pool.len() != n + || reservoir.max_release.len() != n + || reservoir.min_release.len() != n + { + return Err(PyValueError::new_err(format!( + "Reservoir '{}': constraint time series must have length {}", name, n + ))); + } + if !reservoir.max_operating_elevation.is_empty() + && reservoir.max_operating_elevation.len() != n + { + return Err(PyValueError::new_err(format!( + "Reservoir '{}': max_operating_elevation must be empty or length {}", name, n + ))); + } + } + let mut results = HashMap::new(); let mut ordered_objects = Vec::new(); @@ -387,6 +466,10 @@ fn simulate_cascade( target_power: reservoir.target_power.clone(), actual_power: vec![0.0; n], tailwater_elevation: vec![0.0; n], + violation_min_power_pool: vec![0.0; n], + violation_min_release: vec![0.0; n], + violation_max_release: vec![0.0; n], + violation_max_operating_elevation: vec![0.0; n], })); }, ObjectType::River => { @@ -480,6 +563,7 @@ fn simulate_cascade( reservoir.target_power[t], previous_storage, previous_total_outflow_cumecs, + t, ); current_results.storage[t] = state.storage; @@ -491,6 +575,10 @@ fn simulate_cascade( current_results.total_inflow[t] = current_inflow; current_results.actual_power[t] = state.actual_power; current_results.tailwater_elevation[t] = state.tailwater_elev; + current_results.violation_min_power_pool[t] = state.violation_min_power_pool; + current_results.violation_min_release[t] = state.violation_min_release; + current_results.violation_max_release[t] = state.violation_max_release; + current_results.violation_max_operating_elevation[t] = state.violation_max_operating_elevation; // Update previous outflow for next timestep (in Mm3/hr) prev_outflow.insert(name.to_string(), state.release + state.spill); diff --git a/tests/conftest.py b/tests/conftest.py index e1e9b40391b2a14c7f3501b479bec6b04c033c8f..26a1e4c2c28639b75eeb61340e8d22990426f5ea 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -58,6 +58,13 @@ def cumberland_cascade(): else: set_tw_outflow, set_tw_elevation = [], [] + n_hours = len(ts) + # Convert scalar YAML constraint values to time series + specs = dict(specs) # copy to avoid mutating config_dict + specs["min_power_pool"] = [specs["min_power_pool"]] * n_hours + specs["max_release"] = [specs["max_release"]] * n_hours + specs["min_release"] = [specs["min_release"]] * n_hours + reservoir_dict[name] = ReservoirData( **specs, catchment_inflow=ts["catchment_inflow"].tolist(), diff --git a/tests/helpers.py b/tests/helpers.py index cb3d11f8c36df3bb8175c988fdc860d1ec4258ac..238f13b571bee9d0369eb4825698f32267f5f513 100644 --- a/tests/helpers.py +++ b/tests/helpers.py @@ -19,21 +19,22 @@ class ReservoirData: object_type: str capacity: float initial_pool_elevation: float - min_power_pool: float + min_power_pool: list set_storage: list set_elevation: list hpf_h: list hpf_p: list hpf_q: list tailwater_elevation: float - max_release: float - min_release: float + max_release: list + min_release: list catchment_inflow: list target_power: list simulation_order: int downstream_object: str set_tw_outflow: list = field(default_factory=list) set_tw_elevation: list = field(default_factory=list) + max_operating_elevation: list = field(default_factory=list) @dataclass @@ -134,10 +135,10 @@ def make_reservoir(n_hours=24, **overrides): downstream_object="NA", capacity=500.0, initial_pool_elevation=200.0, - min_power_pool=190.0, + min_power_pool=[190.0] * n_hours, tailwater_elevation=150.0, - max_release=2.0, - min_release=0.0, + max_release=[2.0] * n_hours, + min_release=[0.0] * n_hours, set_storage=DEFAULT_SET_STORAGE[:], set_elevation=DEFAULT_SET_ELEVATION[:], catchment_inflow=[0.2] * n_hours, @@ -147,6 +148,7 @@ def make_reservoir(n_hours=24, **overrides): hpf_q=DEFAULT_HPF_Q[:], set_tw_outflow=[], set_tw_elevation=[], + max_operating_elevation=[], ) defaults.update(overrides) return ReservoirData(**defaults) diff --git a/tests/test_cascade.py b/tests/test_cascade.py index 42c736adc7a583702ae2b607dfdde7b3094dcbad..88ece9f5ef309c0ba4a7a72265f52c1036b79263 100644 --- a/tests/test_cascade.py +++ b/tests/test_cascade.py @@ -37,7 +37,9 @@ class TestOutputStructure: expected_fields = [ "storage", "pool_elevation", "head", "release", "spill", "catchment_inflow", "total_inflow", "target_release", - "target_power", "actual_power", + "target_power", "actual_power", "tailwater_elevation", + "violation_min_power_pool", "violation_min_release", + "violation_max_release", "violation_max_operating_elevation", ] for field in expected_fields: assert field in res, f"Missing field: {field}" @@ -111,7 +113,6 @@ class TestTwoReservoirChain: downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -155,7 +156,6 @@ class TestSimulationOrder: downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade(reservoirs={"ResA": res_a, "ResB": res_b}) @@ -186,7 +186,6 @@ class TestDownstreamRouting: n_hours=n, simulation_order=3, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -235,7 +234,6 @@ class TestDownstreamRouting: n_hours=n, simulation_order=6, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -292,12 +290,13 @@ class TestPhysicalInvariants: When below_power_pool, release is forced to 0 regardless of min_release. """ + n = 24 result = run_single_reservoir( - n_hours=24, - min_release=0.05, - max_release=1.5, + n_hours=n, + min_release=[0.05] * n, + max_release=[1.5] * n, ) - for t in range(24): + for t in range(n): r = result["release"][t] if r > 0: assert r >= 0.05 - 1e-9, f"Release below min at t={t}: {r}" @@ -311,7 +310,6 @@ class TestPhysicalInvariants: capacity=350.0, catchment_inflow=[5.0] * 48, target_power=[0.0] * 48, - min_release=0.0, ) for t in range(48): if result["spill"][t] > 1e-9: diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index d5d7780522a6e24396f0682a0fe3f1978bff236d..d8d6c1efeaf7452f43b3eb7a030f13a188609663 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -33,7 +33,6 @@ class TestElevationAtBreakpoints: initial_pool_elevation=195.0, catchment_inflow=[0.0], target_power=[0.0], - min_release=0.0, ) assert result["pool_elevation"][0] == pytest.approx(195.0, abs=0.01) @@ -47,7 +46,6 @@ class TestElevationAtBreakpoints: initial_pool_elevation=180.0, catchment_inflow=[0.0], target_power=[0.0], - min_release=0.0, ) assert result["pool_elevation"][0] == pytest.approx(180.0, abs=0.01) @@ -62,7 +60,6 @@ class TestElevationAtBreakpoints: initial_pool_elevation=210.0, catchment_inflow=[100.0], target_power=[0.0], - min_release=0.0, ) assert result["pool_elevation"][0] == pytest.approx(210.0, abs=0.01) assert result["spill"][0] > 0 # confirms overflow happened @@ -82,7 +79,6 @@ class TestInterpolationMidpoints: initial_pool_elevation=195.0, catchment_inflow=[0.0, 125.0], target_power=[0.0, 0.0], - min_release=0.0, ) # t=0: storage=250, elev=195 assert result["pool_elevation"][0] == pytest.approx(195.0, abs=0.01) @@ -110,8 +106,7 @@ class TestRoundtrip: initial_pool_elevation=elev, catchment_inflow=[0.0], target_power=[0.0], - min_release=0.0, - ) + ) assert result["pool_elevation"][0] == pytest.approx(elev, abs=0.01), \ f"Roundtrip failed for elevation {elev}: got {result['pool_elevation'][0]}" @@ -131,7 +126,6 @@ class TestStorageTracking: initial_pool_elevation=195.0, # storage=250 catchment_inflow=[inflow] * n, target_power=[0.0] * n, - min_release=0.0, ) init_storage = expected_storage(195.0) # 250.0 for t in range(n): @@ -151,7 +145,6 @@ class TestStorageTracking: initial_pool_elevation=195.0, catchment_inflow=[inflow] * n, target_power=[0.0] * n, - min_release=0.0, ) for t in range(n): storage = result["storage"][t] @@ -185,7 +178,6 @@ class TestMonotonicity: initial_pool_elevation=185.0, catchment_inflow=[2.0] * n, target_power=[0.0] * n, - min_release=0.0, ) for t in range(1, n): assert result["pool_elevation"][t] >= result["pool_elevation"][t - 1] - 1e-9, \ @@ -199,7 +191,6 @@ class TestMonotonicity: initial_pool_elevation=205.0, catchment_inflow=[0.0] * n, target_power=[50.0] * n, - min_release=0.0, ) for t in range(1, n): assert result["pool_elevation"][t] <= result["pool_elevation"][t - 1] + 1e-9, \ diff --git a/tests/test_power_clamping.py b/tests/test_power_clamping.py index 73c47caa735f05aba1804a0823e0e8928160b6e9..ea9039da0774d3a19e0fe62e6fa526032e414de3 100644 --- a/tests/test_power_clamping.py +++ b/tests/test_power_clamping.py @@ -21,7 +21,7 @@ def test_actual_power_clamped_to_hpf_max(): result = run_single_reservoir( n_hours=1, target_power=[absurd_power], - max_release=100.0, # high enough to not constrain the release + max_release=[100.0], # high enough to not constrain the release catchment_inflow=[10.0], # plenty of water ) actual = result["actual_power"][0] @@ -38,7 +38,7 @@ def test_actual_power_clamped_moderate_overshoot(): result = run_single_reservoir( n_hours=1, target_power=[150.0], # 50% above HPF max of 100 - max_release=100.0, + max_release=[100.0], catchment_inflow=[10.0], ) actual = result["actual_power"][0] @@ -54,7 +54,7 @@ def test_within_range_power_unchanged(): result = run_single_reservoir( n_hours=1, target_power=[50.0], - max_release=100.0, + max_release=[100.0], catchment_inflow=[10.0], ) actual = result["actual_power"][0] diff --git a/tests/test_regression.py b/tests/test_regression.py index 1a3643a042a79fbb1d5a18a972d1254053cae5f7..4587bb4ecc5a8c502af0add13b0a6e75844ab397 100644 --- a/tests/test_regression.py +++ b/tests/test_regression.py @@ -74,10 +74,10 @@ class TestSmokeTest: downstream_object="NA", capacity=500.0, initial_pool_elevation=200.0, - min_power_pool=190.0, + min_power_pool=[190.0] * 24, tailwater_elevation=150.0, - max_release=2.0, - min_release=0.0, + max_release=[2.0] * 24, + min_release=[0.0] * 24, set_storage=[0.0, 250.0, 500.0], set_elevation=[180.0, 195.0, 210.0], catchment_inflow=[0.2] * 24, diff --git a/tests/test_river_confluence.py b/tests/test_river_confluence.py index 5f6aa5583ac97780dd843740b6781ec6681a6d54..2d047902d346b19558be1d77c5b4485320b0f7b2 100644 --- a/tests/test_river_confluence.py +++ b/tests/test_river_confluence.py @@ -41,7 +41,6 @@ class TestRiverLagRouting: n_hours=n, simulation_order=3, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -80,7 +79,6 @@ class TestRiverLagRouting: n_hours=n, simulation_order=3, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -110,7 +108,6 @@ class TestRiverLagRouting: n_hours=n, simulation_order=3, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -190,7 +187,6 @@ class TestConfluenceRouting: n_hours=n, simulation_order=6, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( @@ -226,7 +222,6 @@ class TestConfluenceRouting: n_hours=n, simulation_order=4, downstream_object="NA", catchment_inflow=[0.0] * n, target_power=[0.0] * n, - min_release=0.0, ) results = powersheds.simulate_cascade( make_cascade( diff --git a/tests/test_simulate_timestep.py b/tests/test_simulate_timestep.py index c3720def4738b6330003ad135b1e4c43455ccbd9..94218f3667914950ff72b0cf4b229f15ad8bc3e0 100644 --- a/tests/test_simulate_timestep.py +++ b/tests/test_simulate_timestep.py @@ -83,7 +83,7 @@ class TestReleaseConstraints: """ result = run_single_reservoir( n_hours=1, - max_release=0.1, + max_release=[0.1], target_power=[100.0], ) assert result["release"][0] == pytest.approx(0.1, abs=1e-9) @@ -100,7 +100,7 @@ class TestReleaseConstraints: """ result = run_single_reservoir( n_hours=1, - min_release=0.1, + min_release=[0.1], target_power=[5.0], ) assert result["release"][0] == pytest.approx(0.1, abs=1e-9) @@ -133,7 +133,7 @@ class TestInsufficientWater: result = run_single_reservoir( n_hours=1, initial_pool_elevation=180.01, - min_power_pool=179.0, + min_power_pool=[179.0], catchment_inflow=[0.0], target_power=[50.0], ) @@ -160,7 +160,7 @@ class TestSpill: capacity=350.0, catchment_inflow=[50.0], target_power=[0.0], - min_release=0.0, + min_release=[0.0], ) init_storage = expected_storage(200.0) expected_spill = init_storage + 50.0 - 350.0 diff --git a/tests/test_time_varying_constraints.py b/tests/test_time_varying_constraints.py new file mode 100644 index 0000000000000000000000000000000000000000..a73a219e1d0bf85fac8ba220dfecd38c9f567ae5 --- /dev/null +++ b/tests/test_time_varying_constraints.py @@ -0,0 +1,462 @@ +""" +Time-Varying Constraints & Max Operating Elevation Tests + +Tests for: +- Time-varying min_power_pool, max_release, min_release +- Max operating elevation constraint (spill trigger) +- Violation output time series +""" + +import pytest + +import powersheds +from helpers import ( + expected_elevation, + expected_storage, + make_cascade, + make_reservoir, + run_single_reservoir, +) + + +class TestTimeVaryingMaxRelease: + """Tests for time-varying max_release constraint.""" + + def test_max_release_ramps_down(self): + """Max release ramps from 2.0 to 0.05 over 12 hours. + + Release should track the max_release constraint when target_release + exceeds the max. + """ + n = 12 + # Large target power -> large target release (above any max_release) + max_rel = [2.0 - i * (1.95 / 11) for i in range(n)] + result = run_single_reservoir( + n_hours=n, + max_release=max_rel, + target_power=[100.0] * n, # high target -> wants lots of release + ) + for t in range(n): + assert result["release"][t] <= max_rel[t] + 1e-9, ( + f"Release {result['release'][t]} exceeds max_release {max_rel[t]} at t={t}" + ) + + +class TestTimeVaryingMinRelease: + """Tests for time-varying min_release constraint.""" + + def test_min_release_increases_midway(self): + """Min release jumps from 0 to 0.15 at t=6. + + With low target power, release should be raised to min_release + in the second half. + """ + n = 12 + min_rel = [0.0] * 6 + [0.15] * 6 + result = run_single_reservoir( + n_hours=n, + min_release=min_rel, + target_power=[5.0] * n, # low target -> wants small release + ) + for t in range(6, n): + assert result["release"][t] >= 0.15 - 1e-9, ( + f"Release {result['release'][t]} below min_release 0.15 at t={t}" + ) + + +class TestTimeVaryingMinPowerPool: + """Tests for time-varying min_power_pool constraint.""" + + def test_min_power_pool_change_affects_release(self): + """Pool starts at 191 (above 190 threshold). At t=6, threshold rises + to 195 (above pool), forcing release=0 for remaining hours. + + initial_pool_elevation=191 -> storage=183.33. + With some target power, pool elevation will drift slightly as + release occurs. When min_power_pool jumps to 195, previous_pool_elevation + will be below 195, so release goes to 0. + """ + n = 12 + min_pp = [190.0] * 6 + [195.0] * 6 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=191.0, + min_power_pool=min_pp, + catchment_inflow=[0.0] * n, + target_power=[50.0] * n, + ) + # In the second half, release should be 0 because pool is below + # the raised min_power_pool=195 + for t in range(7, n): # t=7 onward (t=6 uses previous pool from t=5) + assert result["release"][t] == 0.0, ( + f"Expected release=0 when below min_power_pool at t={t}, " + f"got {result['release'][t]}" + ) + + +class TestMaxOperatingElevationSpill: + """Tests for max_operating_elevation constraint triggering spill.""" + + def test_spill_at_max_operating_elevation(self): + """Set max_operating_elevation=200 (below capacity elevation=210). + High inflow with no power target -> storage builds up. + Spill should occur once pool reaches max_operating_elevation. + """ + n = 24 + # Storage at elevation 200 = (200-180)/30*500 = 333.33 Mm3 + max_op_storage = expected_storage(200.0) + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=198.0, + catchment_inflow=[10.0] * n, + target_power=[0.0] * n, + max_operating_elevation=[200.0] * n, + ) + # At some point, spill should occur + has_spill = any(s > 1e-9 for s in result["spill"]) + assert has_spill, "Expected spill when pool reaches max_operating_elevation" + + # When spill occurs, storage should be at max_op_storage (not capacity) + for t in range(n): + if result["spill"][t] > 1e-9: + assert result["storage"][t] == pytest.approx( + max_op_storage, abs=0.1 + ), ( + f"Storage should be at max_operating level when spilling, " + f"t={t}: storage={result['storage'][t]}, expected={max_op_storage}" + ) + + def test_no_constraint_when_empty(self): + """Empty max_operating_elevation -> same as legacy (only capacity spill).""" + n = 12 + result_with = run_single_reservoir( + n_hours=n, + max_operating_elevation=[], + ) + result_without = run_single_reservoir( + n_hours=n, + ) + for t in range(n): + assert result_with["storage"][t] == pytest.approx( + result_without["storage"][t], abs=1e-9 + ), f"Results differ at t={t} — empty max_operating_elevation should be identical" + + def test_max_operating_tighter_than_capacity(self): + """max_operating_elevation < capacity_elevation -> spill triggers earlier. + + capacity=500 -> elevation=210. Set max_operating_elevation=195 -> + max_op_storage=250. With large inflow, spill should start at 250 not 500. + """ + n = 12 + max_op_storage = expected_storage(195.0) # 250.0 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=194.0, + catchment_inflow=[50.0] * n, + target_power=[0.0] * n, + max_operating_elevation=[195.0] * n, + ) + # Spill should occur and storage should cap at max_op_storage + for t in range(n): + assert result["storage"][t] <= max_op_storage + 0.1, ( + f"Storage {result['storage'][t]} exceeds max_operating storage " + f"{max_op_storage} at t={t}" + ) + + def test_spill_does_not_constrain_release(self): + """When pool is above max_operating_elevation, excess goes to spill. + Release stays within max_release; total outflow (release + spill) may + exceed max_release. + """ + n = 6 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=205.0, + catchment_inflow=[20.0] * n, + target_power=[50.0] * n, + max_release=[0.5] * n, + max_operating_elevation=[195.0] * n, + ) + for t in range(n): + # Release should respect max_release + assert result["release"][t] <= 0.5 + 1e-9, ( + f"Release {result['release'][t]} exceeds max_release at t={t}" + ) + # But spill can add to total outflow + if result["spill"][t] > 1e-9: + total_outflow = result["release"][t] + result["spill"][t] + assert total_outflow > 0.5, ( + f"Expected total outflow > max_release due to spill at t={t}" + ) + + +class TestMassBalanceWithConstraints: + """Mass balance must hold with all new constraints active.""" + + def test_mass_balance_with_max_operating_elevation(self): + """Mass balance: prev_storage + inflow == storage + release + spill.""" + n = 24 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=198.0, + catchment_inflow=[5.0] * n, + target_power=[30.0] * n, + max_operating_elevation=[200.0] * n, + ) + init_storage = expected_storage(198.0) + for t in range(n): + prev = init_storage if t == 0 else result["storage"][t - 1] + balance = ( + prev + result["total_inflow"][t] + - result["storage"][t] + - result["release"][t] + - result["spill"][t] + ) + assert balance == pytest.approx(0.0, abs=1e-9), ( + f"Mass balance violated at t={t}: imbalance={balance}" + ) + + +class TestViolationOutputs: + """Tests for violation time series in output.""" + + def test_violations_zero_when_unconstrained(self): + """Under normal operation, all violations should be 0.""" + result = run_single_reservoir(n_hours=12) + for t in range(12): + assert result["violation_min_power_pool"][t] == 0.0, f"t={t}" + assert result["violation_min_release"][t] == 0.0, f"t={t}" + assert result["violation_max_release"][t] == 0.0, f"t={t}" + assert result["violation_max_operating_elevation"][t] == 0.0, f"t={t}" + + def test_violation_min_power_pool_magnitude(self): + """When pool drops below min_power_pool, violation = difference. + + Set min_power_pool=200, but pool will be lower due to draining. + """ + n = 12 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=198.0, + min_power_pool=[200.0] * n, + catchment_inflow=[0.0] * n, + target_power=[0.0] * n, + ) + # Pool elevation starts at ~198 (with zero release/inflow, stays ~198) + # min_power_pool=200, so violation = 200 - 198 = 2.0 + for t in range(n): + elev = result["pool_elevation"][t] + expected_violation = max(0.0, 200.0 - elev) + assert result["violation_min_power_pool"][t] == pytest.approx( + expected_violation, abs=0.01 + ), f"t={t}: elev={elev}, violation={result['violation_min_power_pool'][t]}" + + def test_violation_min_release_magnitude(self): + """When release < min_release (below_power_pool forces release=0), + violation = min_release - 0 = min_release. + """ + n = 4 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=185.0, # below min_power_pool=190 + min_release=[0.1] * n, + catchment_inflow=[0.5] * n, + target_power=[50.0] * n, + ) + # Release is forced to 0 because below_power_pool + for t in range(n): + if result["release"][t] == 0.0: + assert result["violation_min_release"][t] == pytest.approx( + 0.1, abs=1e-9 + ), f"t={t}" + + def test_violation_max_operating_elevation_magnitude(self): + """Pool above max_operating_elevation should report violation magnitude. + + With max_operating_elevation constraint active, pool is capped, so + violation should be 0 when the constraint fully works. But if we set + a very tight constraint, the pool might still end slightly above due + to interpolation precision. + + Actually, the spill mechanism should prevent pool from exceeding + max_operating_elevation, so the violation should be 0 or near-zero + when the constraint is active. + """ + n = 12 + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=198.0, + catchment_inflow=[10.0] * n, + target_power=[0.0] * n, + max_operating_elevation=[200.0] * n, + ) + # Since we spill excess, pool should not exceed max_operating_elevation + for t in range(n): + assert result["violation_max_operating_elevation"][t] == pytest.approx( + 0.0, abs=0.01 + ), ( + f"t={t}: pool={result['pool_elevation'][t]}, " + f"violation={result['violation_max_operating_elevation'][t]}" + ) + + def test_mixed_violations(self): + """Scenario with multiple constraint activations. + + Phase 1 (t=0-5): normal operation, no violations + Phase 2 (t=6-11): min_power_pool raised high -> below_power_pool + -> release=0 -> min_release violation + """ + n = 12 + min_pp = [190.0] * 6 + [210.0] * 6 # raise threshold above max pool + min_rel = [0.0] * 6 + [0.1] * 6 + result = run_single_reservoir( + n_hours=n, + min_power_pool=min_pp, + min_release=min_rel, + catchment_inflow=[0.2] * n, + target_power=[50.0] * n, + ) + # In phase 2, release should be 0 (below_power_pool) with min_release violation + for t in range(7, n): # t=7+ uses previous pool from t=6+ + assert result["release"][t] == 0.0, f"t={t}" + assert result["violation_min_release"][t] == pytest.approx( + 0.1, abs=1e-9 + ), f"t={t}" + + +class TestConstraintLengthValidation: + """Tests that constraint time series length mismatches raise errors.""" + + def test_wrong_length_min_power_pool_raises(self): + """min_power_pool length != n_hours should raise ValueError.""" + res = make_reservoir(n_hours=24) + res.min_power_pool = [190.0] * 10 # wrong length + cascade = make_cascade(reservoirs={"Res": res}) + with pytest.raises(Exception, match="constraint time series must have length"): + powersheds.simulate_cascade(cascade) + + def test_wrong_length_max_operating_elevation_raises(self): + """Non-empty max_operating_elevation with wrong length should raise.""" + res = make_reservoir(n_hours=24) + res.max_operating_elevation = [200.0] * 10 # wrong length + cascade = make_cascade(reservoirs={"Res": res}) + with pytest.raises(Exception, match="max_operating_elevation must be empty or length"): + powersheds.simulate_cascade(cascade) + + +class TestNaNConstraints: + """Tests that NaN constraint values are treated as 'no constraint'.""" + + def test_nan_max_release_means_no_limit(self): + """NaN in max_release should behave like no upper bound on release.""" + n = 6 + nan = float("nan") + result_nan = run_single_reservoir( + n_hours=n, + max_release=[nan] * n, + target_power=[100.0] * n, + ) + result_high = run_single_reservoir( + n_hours=n, + max_release=[1e6] * n, + target_power=[100.0] * n, + ) + for t in range(n): + assert result_nan["release"][t] == pytest.approx( + result_high["release"][t], abs=1e-9 + ), f"t={t}: NaN max_release should not constrain release" + + def test_nan_min_release_means_no_floor(self): + """NaN in min_release should behave like min_release=0.""" + n = 6 + nan = float("nan") + result_nan = run_single_reservoir( + n_hours=n, + min_release=[nan] * n, + target_power=[5.0] * n, + ) + result_zero = run_single_reservoir( + n_hours=n, + min_release=[0.0] * n, + target_power=[5.0] * n, + ) + for t in range(n): + assert result_nan["release"][t] == pytest.approx( + result_zero["release"][t], abs=1e-9 + ), f"t={t}: NaN min_release should behave like 0" + + def test_nan_min_power_pool_means_no_check(self): + """NaN in min_power_pool should never trigger below-power-pool. + + With a very low pool elevation (185) that is normally below + min_power_pool=190, NaN should allow generation to proceed. + """ + n = 4 + nan = float("nan") + result = run_single_reservoir( + n_hours=n, + initial_pool_elevation=185.0, + min_power_pool=[nan] * n, + catchment_inflow=[0.5] * n, + target_power=[50.0] * n, + ) + # Should NOT be forced to zero — NaN means no power pool check + for t in range(n): + assert result["release"][t] > 0, ( + f"t={t}: NaN min_power_pool should not suppress release" + ) + + def test_nan_max_operating_elevation_means_no_spill_cap(self): + """NaN in max_operating_elevation should behave like no constraint + (only capacity triggers spill).""" + n = 12 + nan = float("nan") + result_nan = run_single_reservoir( + n_hours=n, + max_operating_elevation=[nan] * n, + ) + result_empty = run_single_reservoir( + n_hours=n, + max_operating_elevation=[], + ) + for t in range(n): + assert result_nan["storage"][t] == pytest.approx( + result_empty["storage"][t], abs=1e-9 + ), f"t={t}: NaN max_operating_elevation should match empty" + + def test_nan_violations_are_zero(self): + """All violation outputs should be 0.0 when the constraint is NaN.""" + n = 6 + nan = float("nan") + result = run_single_reservoir( + n_hours=n, + min_power_pool=[nan] * n, + min_release=[nan] * n, + max_release=[nan] * n, + max_operating_elevation=[nan] * n, + ) + for t in range(n): + assert result["violation_min_power_pool"][t] == 0.0, f"t={t}" + assert result["violation_min_release"][t] == 0.0, f"t={t}" + assert result["violation_max_release"][t] == 0.0, f"t={t}" + assert result["violation_max_operating_elevation"][t] == 0.0, f"t={t}" + + def test_mixed_nan_and_real_constraints(self): + """NaN at some timesteps, real values at others. + + First 3 hours: max_release=NaN (no limit), last 3: max_release=0.05. + """ + n = 6 + nan = float("nan") + max_rel = [nan, nan, nan, 0.05, 0.05, 0.05] + result = run_single_reservoir( + n_hours=n, + max_release=max_rel, + target_power=[100.0] * n, + ) + # First 3 hours: release unconstrained by max_release + # Last 3 hours: release clamped to 0.05 + for t in range(3, n): + assert result["release"][t] <= 0.05 + 1e-9, ( + f"t={t}: release should be clamped to 0.05" + ) diff --git a/uv.lock b/uv.lock index 9f4582464ef075c9e88b3e7484e2d64c92dff606..564504c1ee0a381e061dfbb18782d79bca1c57aa 100644 --- a/uv.lock +++ b/uv.lock @@ -495,6 +495,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/af/33/ee4519fa02ed11a94aef9559552f3b17bb863f2ecfe1a35dc7f548cde231/matplotlib_inline-0.2.1-py3-none-any.whl", hash = "sha256:d56ce5156ba6085e00a9d54fead6ed29a9c47e215cd1bba2e976ef39f5710a76", size = 9516, upload-time = "2025-10-23T09:00:20.675Z" }, ] +[[package]] +name = "narwhals" +version = "2.16.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/fc/6f/713be67779028d482c6e0f2dde5bc430021b2578a4808c1c9f6d7ad48257/narwhals-2.16.0.tar.gz", hash = "sha256:155bb45132b370941ba0396d123cf9ed192bf25f39c4cea726f2da422ca4e145", size = 618268, upload-time = "2026-02-02T10:31:00.545Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/03/cc/7cb74758e6df95e0c4e1253f203b6dd7f348bf2f29cf89e9210a2416d535/narwhals-2.16.0-py3-none-any.whl", hash = "sha256:846f1fd7093ac69d63526e50732033e86c30ea0026a44d9b23991010c7d1485d", size = 443951, upload-time = "2026-02-02T10:30:58.635Z" }, +] + [[package]] name = "nest-asyncio" version = "1.6.0" @@ -681,6 +690,19 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/cb/28/3bfe2fa5a7b9c46fe7e13c97bda14c895fb10fa2ebf1d0abb90e0cea7ee1/platformdirs-4.5.1-py3-none-any.whl", hash = "sha256:d03afa3963c806a9bed9d5125c8f4cb2fdaf74a55ab60e5d59b3fde758104d31", size = 18731, upload-time = "2025-12-05T13:52:56.823Z" }, ] +[[package]] +name = "plotly" +version = "6.5.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "narwhals" }, + { name = "packaging" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/e3/4f/8a10a9b9f5192cb6fdef62f1d77fa7d834190b2c50c0cd256bd62879212b/plotly-6.5.2.tar.gz", hash = "sha256:7478555be0198562d1435dee4c308268187553cc15516a2f4dd034453699e393", size = 7015695, upload-time = "2026-01-14T21:26:51.222Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/8a/67/f95b5460f127840310d2187f916cf0023b5875c0717fdf893f71e1325e87/plotly-6.5.2-py3-none-any.whl", hash = "sha256:91757653bd9c550eeea2fa2404dba6b85d1e366d54804c340b2c874e5a7eb4a4", size = 9895973, upload-time = "2026-01-14T21:26:47.135Z" }, +] + [[package]] name = "pluggy" version = "1.6.0" @@ -703,6 +725,11 @@ dependencies = [ { name = "scipy" }, ] +[package.optional-dependencies] +viz = [ + { name = "plotly" }, +] + [package.dev-dependencies] dev = [ { name = "pytest" }, @@ -714,10 +741,12 @@ requires-dist = [ { name = "matplotlib", specifier = ">=3.10" }, { name = "numpy", specifier = ">=2.2,<2.3" }, { name = "pandas", specifier = ">=2.2,<2.3" }, + { name = "plotly", marker = "extra == 'viz'", specifier = ">=5.18" }, { name = "pyarrow", specifier = "==21.0" }, { name = "pyyaml", specifier = "==6.0.2" }, { name = "scipy", specifier = ">=1.15" }, ] +provides-extras = ["viz"] [package.metadata.requires-dev] dev = [{ name = "pytest", specifier = ">=9.0.2" }]