Fix actual_power not clamped to HPF table bounds (513c58bb) · Commits · HydroCHiPPs / Powersheds

src/helpers.rs

+42 −10

Original line number	Diff line number	Diff line
		@@ -81,10 +81,18 @@ pub fn interpolate_storage(
		}


		/// HELPER FUNCTION: compute_release_to_meet_target_power
		/// Result of computing release to meet a target power.
		/// Contains both the volumetric release and the effective (possibly clamped) power.
		pub struct ReleaseResult {
		pub release: f64,
		pub effective_power: f64,
		}

		/// HELPER FUNCTION: compute_release_and_effective_power
		///
		/// Computes volumetric release (Mm3 over 1 hour) required to meet target power,
		/// using bilinear interpolation on a regular (H, P) grid with Q values.
		/// Also returns the effective power, which may differ from target if clamped to HPF bounds.
		///
		/// Inputs:
		/// - `power_mw`: target power (MW)
		@@ -111,17 +119,17 @@ pub fn interpolate_storage(
		/// Notes:
		/// - Converts cumecs to Mm3 over one hour using existing constants:
		/// `SECONDS_IN_HOUR` and `MM3_IN_M3`
		pub fn compute_release_to_meet_target_power(
		pub fn compute_release_and_effective_power(
		power_mw: f64,
		head_m: f64,
		hpf_head_m: &[f64],
		hpf_power_mw: &[f64],
		hpf_flow_cumecs: &[f64],
		) -> f64 {
		) -> ReleaseResult {

		// target power of zero returns target release of zero
		if power_mw == 0.0 {
		return 0.0;
		return ReleaseResult { release: 0.0, effective_power: 0.0 };
		}

		debug_assert_eq!(hpf_head_m.len(), hpf_power_mw.len());
		@@ -150,7 +158,7 @@ pub fn compute_release_to_meet_target_power(
		unique_p.dedup_by(\|a, b\| approx_eq(a, b));

		if unique_h.is_empty() \|\| unique_p.is_empty() {
		return f64::NAN;
		return ReleaseResult { release: f64::NAN, effective_power: f64::NAN };
		}

		// 2) Clamp H0 and P0 to grid bounds and find bracket indices.
		@@ -224,7 +232,7 @@ pub fn compute_release_to_meet_target_power(
		// If no feasible power levels exist at this head, return 0 (can't operate)
		let (min_idx, max_idx) = match (min_feasible_p_idx, max_feasible_p_idx) {
		(Some(min), Some(max)) => (min, max),
		_ => return 0.0,
		_ => return ReleaseResult { release: 0.0, effective_power: 0.0 },
		};

		let min_feasible_p = unique_p[min_idx];
		@@ -235,14 +243,20 @@ pub fn compute_release_to_meet_target_power(
		if power_mw < min_feasible_p {
		let t = if approx_eq(h_lo, h_hi) { 0.0 } else { (h0 - h_lo) / (h_hi - h_lo) };
		let q_cumecs = lerp(q_at_h_lo[min_idx], q_at_h_hi[min_idx], t);
		return q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3;
		return ReleaseResult {
		release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,
		effective_power: min_feasible_p,
		};
		}

		// If target power is above maximum feasible, clamp to max feasible
		if power_mw > max_feasible_p {
		let t = if approx_eq(h_lo, h_hi) { 0.0 } else { (h0 - h_lo) / (h_hi - h_lo) };
		let q_cumecs = lerp(q_at_h_lo[max_idx], q_at_h_hi[max_idx], t);
		return q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3;
		return ReleaseResult {
		release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,
		effective_power: max_feasible_p,
		};
		}

		// Target power is within feasible range - find bracketing P indices
		@@ -270,7 +284,10 @@ pub fn compute_release_to_meet_target_power(
		// Fallback: use min feasible power
		let t = if approx_eq(h_lo, h_hi) { 0.0 } else { (h0 - h_lo) / (h_hi - h_lo) };
		let q_cumecs = lerp(q_at_h_lo[min_idx], q_at_h_hi[min_idx], t);
		return q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3;
		return ReleaseResult {
		release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,
		effective_power: min_feasible_p,
		};
		}

		let p_lo = unique_p[valid_pi_lo];
		@@ -302,7 +319,22 @@ pub fn compute_release_to_meet_target_power(

		// 5) Convert cumecs → Mm3 over 1 hour.
		let volumetric_release_mm3 = q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3;
		volumetric_release_mm3
		ReleaseResult {
		release: volumetric_release_mm3,
		effective_power: power_mw,
		}
		}


		/// Thin wrapper returning only the release value (preserves existing API).
		pub fn compute_release_to_meet_target_power(
		power_mw: f64,
		head_m: f64,
		hpf_head_m: &[f64],
		hpf_power_mw: &[f64],
		hpf_flow_cumecs: &[f64],
		) -> f64 {
		compute_release_and_effective_power(power_mw, head_m, hpf_head_m, hpf_power_mw, hpf_flow_cumecs).release
		}

src/lib.rs

+5 −2

Original line number	Diff line number	Diff line
		@@ -8,6 +8,7 @@ mod helpers;
		use helpers::interpolate_elevation;
		use helpers::interpolate_storage;
		use helpers::compute_release_to_meet_target_power;
		use helpers::compute_release_and_effective_power;
		use helpers::compute_power_given_actual_release;


		@@ -146,13 +147,15 @@ fn simulate_timestep(
		let head = previous_pool_elevation - reservoir.tailwater_elevation;

		// Compute the target release given the power generation desired
		let target_release = compute_release_to_meet_target_power(
		let release_result = compute_release_and_effective_power(
		target_power,
		head,
		&reservoir.hpf_h,
		&reservoir.hpf_p,
		&reservoir.hpf_q
		);
		let target_release = release_result.release;
		let effective_power = release_result.effective_power;

		// Compute available water given inflow and previous storage...
		// ... then determine actual release and resulting storage
		@@ -180,7 +183,7 @@ fn simulate_timestep(
		)
		} else {

		target_power
		effective_power

		};

tests/fixtures/cumberland_baseline.json

+5 −5

Original line number	Diff line number	Diff line
		{
		"WolfCreek": {
		"final_storage": 5086.765250757235,
		"sum_actual_power": 14509.02,
		"sum_actual_power": 14541.0,
		"sum_spill": 0.0,
		"mean_pool_elevation": 220.33568136836752
		},
		@@ -19,7 +19,7 @@
		},
		"CenterHill": {
		"final_storage": 2719.2148359068806,
		"sum_actual_power": 3635.94,
		"sum_actual_power": 3962.94,
		"sum_spill": 0.0,
		"mean_pool_elevation": 210.18284075212327
		},
		@@ -37,13 +37,13 @@
		},
		"Cheatham": {
		"final_storage": 594.2905274140278,
		"sum_actual_power": 2270.82,
		"sum_actual_power": 1563.74,
		"sum_spill": 0.0,
		"mean_pool_elevation": 123.3456746352141
		},
		"Barkley": {
		"final_storage": 1423.766365644019,
		"sum_actual_power": 5460.54,
		"sum_actual_power": 4733.753333333333,
		"sum_spill": 0.0,
		"mean_pool_elevation": 110.13269757075129
		}

tests/test_power_clamping.py

0 → 100644

+63 −0

Original line number	Diff line number	Diff line
		"""
		Test that actual power is properly clamped to the HPF table maximum.

		The HPF table defines a finite range of achievable power at each head.
		If a user requests a target power above the HPF maximum, actual_power
		must be clamped to the HPF maximum — not echo back the impossible target.
		"""

		import powersheds
		from tests.helpers import make_reservoir, make_cascade, run_single_reservoir


		# Default HPF table: 2 heads (40, 60) x 3 powers (0, 50, 100)
		# Max power in the table is 100 MW.
		# Default reservoir: tailwater=150, initial_pool_elevation=200 → head≈50m


		def test_actual_power_clamped_to_hpf_max():
		"""Requesting power far above HPF max should NOT produce actual_power > max."""
		absurd_power = 99999.0
		result = run_single_reservoir(
		n_hours=1,
		target_power=[absurd_power],
		max_release=100.0, # high enough to not constrain the release
		catchment_inflow=[10.0], # plenty of water
		)
		actual = result["actual_power"][0]
		hpf_max = 100.0 # max power in default HPF table

		assert actual <= hpf_max, (
		f"actual_power ({actual}) exceeds HPF max ({hpf_max}). "
		f"The simulator echoed back an impossible target_power."
		)


		def test_actual_power_clamped_moderate_overshoot():
		"""Requesting power slightly above HPF max should clamp to max."""
		result = run_single_reservoir(
		n_hours=1,
		target_power=[150.0], # 50% above HPF max of 100
		max_release=100.0,
		catchment_inflow=[10.0],
		)
		actual = result["actual_power"][0]
		hpf_max = 100.0

		assert actual <= hpf_max, (
		f"actual_power ({actual}) exceeds HPF max ({hpf_max})."
		)


		def test_within_range_power_unchanged():
		"""Target power within the HPF table should be achievable as-is."""
		result = run_single_reservoir(
		n_hours=1,
		target_power=[50.0],
		max_release=100.0,
		catchment_inflow=[10.0],
		)
		actual = result["actual_power"][0]
		assert abs(actual - 50.0) < 1.0, (
		f"Expected actual_power ≈ 50 MW, got {actual}"
		)