expand violation reporting

📈 Dynamic release and power computation based on head-power-flow tables that can be easily constructed using efficiency curves specified for individual units at each plant (see offline, maximum efficiency DP model for unit commitment)

⚡ Realized hydropower and power violations determined after accounting for water availbility, pool level, and release constraints

⚡ Realized hydropower, target-power deviations, and limiting-cause diagnostics after accounting for water availability, pool level, and release constraints

⏩ Fast run times (~0.05s for an 8 reservoir network simulated over 168 hours) allow for thousands of simulations on desktop machines, paving the way for many interations within a model coupling environment, or for uncertainty analysis using probabilistic inflow forecasts

print("non-power release (Mm3/hr)", res["non_power_release"][:5])

print("storage (Mm3)",    res["storage"][:5])

print("power (MW)",       res["actual_power"][:5])

print("limiters",         res["limiting_cause"][:5])

```

#### Notes

Returned per-object results:

- Reservoir: `storage`, `pool_elevation`, `head`, `release`, `non_power_release`, `total_controlled_outflow`, `total_outflow`, `spill`, `catchment_inflow`, `total_inflow`, `target_release`, `target_power`, `actual_power`, `tailwater_elevation`, `violation_min_power_pool`, `violation_min_release`, `violation_max_release`, `violation_min_non_power_release`, `violation_max_non_power_release`, `violation_max_operating_elevation`

- Reservoir: `storage`, `pool_elevation`, `head`, `release`, `non_power_release`, `total_controlled_outflow`, `total_outflow`, `spill`, `catchment_inflow`, `total_inflow`, `target_release`, `target_power`, `actual_power`, `tailwater_elevation`, `violation_min_power_pool`, `violation_min_release`, `violation_max_release`, `violation_min_non_power_release`, `violation_max_non_power_release`, `violation_max_operating_elevation`, `violation_target_power_shortfall`, `violation_target_power_excess`, `limiting_cause`

Notes:

- Input `min_release` / `max_release` constrain total controlled outflow, not passive spill.

- Output `release` and `target_release` remain turbine-only series.

- `violation_target_power_shortfall` records unmet target power even when end-of-step constraint violations are zero.

- `violation_target_power_excess` records over-generation above the requested target, such as when minimum release constraints force extra turbine flow.

- `limiting_cause` is a per-timestep summary label such as `none`, `below_power_pool`, `water_availability`, `max_release`, `min_release`, `min_non_power_release`, `hpf_max_power`, `hpf_min_power`, or `max_operating_elevation`.

- River: `inflow`, `outflow` (with travel-time lag)

- Confluence: `inflow`, `outflow` (no lag)

const MM3_IN_M3: f64 = 1_000_000.0;

// const MEGAWATT_TO_WATT: f64 = 1_000_000.0;

/// Elevation computation based on storage.

/// Performs linear interpolation to find a y-value (elevation)

/// corresponding to an x-value (storage)

pub fn interpolate_elevation(

    storage: f64,

    storage_values: &[f64],

    elevation_values: &[f64]) -> f64 {

    elevation_values: &[f64],

) -> f64 {

    // Handle edge cases

    if storage_values.is_empty() || elevation_values.is_empty() {

        return 0.0; // Default value if tables are empty

pub fn interpolate_storage(

    elevation: f64,

    storage_values: &[f64],

    elevation_values: &[f64]) -> f64 {

    elevation_values: &[f64],

) -> f64 {

    // Handle edge cases

    if storage_values.is_empty() || elevation_values.is_empty() {

        return 0.0; // Default value if tables are empty

    storage_values[0] // Fallback

}

/// Tailwater elevation lookup from outflow-TW rating curve.

/// Performs linear interpolation to find tailwater elevation (meters)

/// corresponding to total outflow (cumecs).

    tw_elevation_values[0]

}

/// Result of computing release to meet a target power.

/// Contains both the volumetric release and the effective (possibly clamped) power.

pub struct ReleaseResult {

    hpf_power_mw: &[f64],

    hpf_flow_cumecs: &[f64],

) -> ReleaseResult {

    // target power of zero returns target release of zero

    if power_mw == 0.0 {

        return ReleaseResult { release: 0.0, effective_power: 0.0 };

        return ReleaseResult {

            release: 0.0,

            effective_power: 0.0,

        };

    }

    debug_assert_eq!(hpf_head_m.len(), hpf_power_mw.len());

    unique_p.dedup_by(|a, b| approx_eq(*a, *b));

    if unique_h.is_empty() || unique_p.is_empty() {

        return ReleaseResult { release: f64::NAN, effective_power: f64::NAN };

        return ReleaseResult {

            release: f64::NAN,

            effective_power: f64::NAN,

        };

    }

    // 2) Clamp H0 and P0 to grid bounds and find bracket indices.

    // 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 ReleaseResult { release: 0.0, effective_power: 0.0 },

        _ => {

            return ReleaseResult {

                release: 0.0,

                effective_power: 0.0,

            }

        }

    };

    let min_feasible_p = unique_p[min_idx];

    // Symmetric clamping for infeasible power targets (maintains reversibility)

    // If target power is below minimum feasible, clamp to min feasible

    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 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 ReleaseResult {

            release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,

    // 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 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 ReleaseResult {

            release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,

    let mut valid_pi_hi = pi_hi;

    // Search backward for a valid lower bound

    while valid_pi_lo > 0 &&

          !(q_at_h_lo[valid_pi_lo].is_finite() && q_at_h_hi[valid_pi_lo].is_finite()) {

    while valid_pi_lo > 0

        && !(q_at_h_lo[valid_pi_lo].is_finite() && q_at_h_hi[valid_pi_lo].is_finite())

    {

        valid_pi_lo -= 1;

    }

    // Search forward for a valid upper bound

    while valid_pi_hi < n_p - 1 &&

          !(q_at_h_lo[valid_pi_hi].is_finite() && q_at_h_hi[valid_pi_hi].is_finite()) {

    while valid_pi_hi < n_p - 1

        && !(q_at_h_lo[valid_pi_hi].is_finite() && q_at_h_hi[valid_pi_hi].is_finite())

    {

        valid_pi_hi += 1;

    }

    // Ensure we have valid corners

    if !(q_at_h_lo[valid_pi_lo].is_finite() && q_at_h_hi[valid_pi_lo].is_finite() &&

         q_at_h_lo[valid_pi_hi].is_finite() && q_at_h_hi[valid_pi_hi].is_finite()) {

    if !(q_at_h_lo[valid_pi_lo].is_finite()

        && q_at_h_hi[valid_pi_lo].is_finite()

        && q_at_h_lo[valid_pi_hi].is_finite()

        && q_at_h_hi[valid_pi_hi].is_finite())

    {

        // 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 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 ReleaseResult {

            release: q_cumecs * SECONDS_IN_HOUR / MM3_IN_M3,

        q11

    } else if approx_eq(h_lo, h_hi) {

        // linear in P only

        let u = if approx_eq(p_lo, p_hi) { 0.0 } else { (p0 - p_lo) / (p_hi - p_lo) };

        let u = if approx_eq(p_lo, p_hi) {

            0.0

        } else {

            (p0 - p_lo) / (p_hi - p_lo)

        };

        lerp(q11, q12, u)

    } else if approx_eq(p_lo, p_hi) {

        // linear in H only

        let t = if approx_eq(h_lo, h_hi) { 0.0 } else { (h0 - h_lo) / (h_hi - h_lo) };

        let t = if approx_eq(h_lo, h_hi) {

            0.0

        } else {

            (h0 - h_lo) / (h_hi - h_lo)

        };

        lerp(q11, q21, t)

    } else {

        let t = (h0 - h_lo) / (h_hi - h_lo);

    }

}

/// Thin wrapper returning only the release value (preserves existing API).

pub fn compute_release_to_meet_target_power(

    power_mw: 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

    compute_release_and_effective_power(power_mw, head_m, hpf_head_m, hpf_power_mw, hpf_flow_cumecs)

        .release

}

/// HELPER FUNCTION: compute_power_given_actual_release

///

/// Computes power (MW) given volumetric release (Mm3 over 1 hour) and head (m),

        match grid.binary_search_by(|v| v.partial_cmp(&x).unwrap()) {

            Ok(i) => (i, i),

            Err(i) => {

                if i == 0 { (0, 0) }

                else if i >= grid.len() { let j = grid.len() - 1; (j, j) }

                else { (i - 1, i) }

                if i == 0 {

                    (0, 0)

                } else if i >= grid.len() {

                    let j = grid.len() - 1;

                    (j, j)

                } else {

                    (i - 1, i)

                }

            }

        }

    };

    let h_hi = unique_h[hi_hi];

    // Compute H interpolation factor

    let t = if approx_eq(h_lo, h_hi) { 0.0 } else { (h0 - h_lo) / (h_hi - h_lo) };

    let t = if approx_eq(h_lo, h_hi) {

        0.0

    } else {

        (h0 - h_lo) / (h_hi - h_lo)

    };

    // 3) Build a single-pass lookup: collect Q at (h_lo, P) and (h_hi, P) for all P

    //    Use a hashmap-like approach with indices into unique_p

        // Find P index using binary search

        if let Ok(p_idx) = unique_p.binary_search_by(|v| {

            if approx_eq(*v, pv) { std::cmp::Ordering::Equal }

            else { v.partial_cmp(&pv).unwrap() }

            if approx_eq(*v, pv) {

                std::cmp::Ordering::Equal

            } else {

                v.partial_cmp(&pv).unwrap()

            }

        }) {

            if approx_eq(h, h_lo) {

                q_at_h_lo[p_idx] = qv;

    let u = (q0 - a) / b;

    lerp(p_lo, p_hi, u)

}

            "violation_min_power_pool", "violation_min_release",

            "violation_max_release", "violation_min_non_power_release",

            "violation_max_non_power_release", "violation_max_operating_elevation",

            "violation_target_power_shortfall", "violation_target_power_excess",

            "limiting_cause",

        ]

        for field in expected_fields:

            assert field in res, f"Missing field: {field}"

        for name, result in cumberland_results.items():

            for field, values in result.items():

                for i, v in enumerate(values):

                    if isinstance(v, (int, float)):

                        assert not math.isnan(v), f"NaN in {name}.{field}[{i}]"

        f"actual_power ({actual}) exceeds HPF max ({hpf_max}). "

        f"The simulator echoed back an impossible target_power."

    )

    assert result["violation_target_power_shortfall"][0] == absurd_power - actual

    assert result["violation_target_power_excess"][0] == 0.0

    assert result["limiting_cause"][0] == "hpf_max_power"

def test_actual_power_clamped_moderate_overshoot():