oxigrid 0.1.2

//! Voltage Stability Assessment (VSA) — advanced indices, modal analysis, and fast screening.
//!
//! Implements:
//! - **L-index** (Kessel & Glavitsch, 1986): per-load-bus proximity-to-collapse index
//! - **Reduced Jacobian Modal Analysis**: minimum singular value and bus participation factors
//! - **Fast Voltage Stability Index (FVSI)**: per-branch line stability index
//! - **Line Stability Index (Lmn)**: per-line PSI-based index
//! - **N-1 Voltage Stability Screening**: contingency ranking by severity
//! - **Voltage Stability Margin**: continuation-based λ_max estimation

use crate::error::{OxiGridError, Result};
use crate::network::{BusType, PowerNetwork};
use crate::powerflow::{PowerFlowConfig, PowerFlowMethod, PowerFlowResult};
use num_complex::Complex64;
use serde::{Deserialize, Serialize};

// ─── L-Index (Kessel-Glavitsch) ───────────────────────────────────────────────

/// Result of the Kessel-Glavitsch L-index computation.
///
/// L_j = 0 → no load (ideal); L_j → 1 → voltage collapse imminent.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct LIndexResult {
    /// Per load-bus L-index values (length = number of PQ buses).
    pub l_index: Vec<f64>,
    /// System L-index = max(L_j) across all load buses.
    pub l_max: f64,
    /// Internal bus index with the highest L-index.
    pub critical_bus: usize,
    /// `true` when `l_max < 1.0` (power flow still feasible).
    pub stable: bool,
    /// Stability headroom: `1.0 - l_max`.
    pub stability_margin: f64,
    /// Load-bus internal indices in the same order as `l_index`.
    pub load_bus_indices: Vec<usize>,
}

/// Compute the Kessel-Glavitsch L-index from a power-flow solution and Y-bus.
///
/// # Theory
/// Partition the admittance matrix:
/// ```text
///   [I_L]   [Y_LL  Y_LG] [V_L]
///   [I_G] = [Y_GL  Y_GG] [V_G]
/// ```
/// The transfer matrix `F_LG = -Y_LL⁻¹ · Y_LG` relates generator voltages to
/// load bus voltages.  For load bus j:
/// ```text
///   L_j = |Σ_{i∈G} F_{ji} · V_i / V_j|
/// ```
/// `L_j ∈ [0, 1)` for a stable system; `L_j → 1` at voltage collapse.
pub fn compute_l_index(
    network: &PowerNetwork,
    v_mag: &[f64],
    v_ang: &[f64],
    ybus: &[Vec<Complex64>],
) -> Result<LIndexResult> {
    let n = network.buses.len();
    if n == 0 {
        return Err(OxiGridError::InvalidNetwork("empty network".into()));
    }
    if v_mag.len() < n || v_ang.len() < n {
        return Err(OxiGridError::InvalidParameter(
            "v_mag/v_ang length mismatch with bus count".into(),
        ));
    }

    // Identify generator buses (Slack + PV) and load buses (PQ)
    let gen_indices: Vec<usize> = network
        .buses
        .iter()
        .enumerate()
        .filter(|(_, b)| b.bus_type == BusType::Slack || b.bus_type == BusType::PV)
        .map(|(i, _)| i)
        .collect();

    let load_indices: Vec<usize> = network
        .buses
        .iter()
        .enumerate()
        .filter(|(_, b)| b.bus_type == BusType::PQ)
        .map(|(i, _)| i)
        .collect();

    let n_l = load_indices.len();
    let n_g = gen_indices.len();

    if n_l == 0 {
        return Ok(LIndexResult {
            l_index: vec![],
            l_max: 0.0,
            critical_bus: 0,
            stable: true,
            stability_margin: 1.0,
            load_bus_indices: vec![],
        });
    }

    // Build a mapping from internal index → row/col position in Y_LL and Y_LG
    let mut load_pos = vec![usize::MAX; n];
    for (pos, &idx) in load_indices.iter().enumerate() {
        load_pos[idx] = pos;
    }
    let mut gen_pos = vec![usize::MAX; n];
    for (pos, &idx) in gen_indices.iter().enumerate() {
        gen_pos[idx] = pos;
    }

    // Extract Y_LL (n_L × n_L) and Y_LG (n_L × n_G) from full Y-bus
    let mut y_ll: Vec<Vec<Complex64>> = vec![vec![Complex64::new(0.0, 0.0); n_l]; n_l];
    let mut y_lg: Vec<Vec<Complex64>> = vec![vec![Complex64::new(0.0, 0.0); n_g]; n_l];

    for (row_l, &i) in load_indices.iter().enumerate() {
        if i >= ybus.len() {
            continue;
        }
        for (j, &y_val) in ybus[i].iter().enumerate() {
            if load_pos[j] != usize::MAX {
                y_ll[row_l][load_pos[j]] = y_val;
            } else if gen_pos[j] != usize::MAX {
                y_lg[row_l][gen_pos[j]] = y_val;
            }
        }
    }

    // Compute F_LG = -Y_LL⁻¹ · Y_LG via Gaussian elimination (augmented matrix)
    // Augment Y_LL with Y_LG columns: [Y_LL | Y_LG]
    let rhs_cols = n_g;
    let mut aug: Vec<Vec<Complex64>> = (0..n_l)
        .map(|i| {
            let mut row = y_ll[i].clone();
            row.extend_from_slice(&y_lg[i]);
            row
        })
        .collect();

    // Forward elimination with partial pivoting
    for col in 0..n_l {
        // Find pivot row
        let pivot_row = (col..n_l)
            .max_by(|&r1, &r2| {
                aug[r1][col]
                    .norm()
                    .partial_cmp(&aug[r2][col].norm())
                    .unwrap_or(std::cmp::Ordering::Equal)
            })
            .unwrap_or(col);

        if aug[pivot_row][col].norm() < 1e-14 {
            // Singular or near-singular Y_LL — fall back to diagonal approximation
            return compute_l_index_diagonal_fallback(
                network,
                v_mag,
                v_ang,
                &load_indices,
                &gen_indices,
            );
        }

        aug.swap(col, pivot_row);

        let pivot = aug[col][col];
        #[allow(clippy::needless_range_loop)]
        for j in col..n_l + rhs_cols {
            aug[col][j] /= pivot;
        }

        for row in 0..n_l {
            if row == col {
                continue;
            }
            let factor = aug[row][col];
            #[allow(clippy::needless_range_loop)]
            for j in col..n_l + rhs_cols {
                let sub = factor * aug[col][j];
                aug[row][j] -= sub;
            }
        }
    }

    // Extract Y_LL⁻¹ · Y_LG from the augmented right-hand side
    // F_LG = -Y_LL⁻¹ · Y_LG
    let mut f_lg: Vec<Vec<Complex64>> = vec![vec![Complex64::new(0.0, 0.0); n_g]; n_l];
    for i in 0..n_l {
        for j in 0..n_g {
            f_lg[i][j] = -aug[i][n_l + j];
        }
    }

    // Compute L-index for each load bus using Kessel-Glavitsch formula:
    //   L_j = |Σ_{i∈G} F_ji · V_i / V_j  -  1|
    //
    // At no-load: Σ F_ji V_i = V_j  →  sum/V_j = 1  →  L_j = 0  (stable)
    // At collapse: L_j → 1
    let mut l_index = vec![0.0f64; n_l];
    for (j, &load_bus) in load_indices.iter().enumerate() {
        let v_j = Complex64::new(
            v_mag[load_bus] * v_ang[load_bus].cos(),
            v_mag[load_bus] * v_ang[load_bus].sin(),
        );
        if v_j.norm() < 1e-9 {
            l_index[j] = 1.0;
            continue;
        }

        let mut sum = Complex64::new(0.0, 0.0);
        for (k, &gen_bus) in gen_indices.iter().enumerate() {
            let v_i = Complex64::new(
                v_mag[gen_bus] * v_ang[gen_bus].cos(),
                v_mag[gen_bus] * v_ang[gen_bus].sin(),
            );
            sum += f_lg[j][k] * v_i / v_j;
        }
        // L_j = |F * V_G / V_L - 1|: 0 at no-load, 1 at collapse
        l_index[j] = (sum - Complex64::new(1.0, 0.0)).norm().min(2.0);
    }

    // Find maximum and critical bus
    let (max_pos, &l_max_val) = l_index
        .iter()
        .enumerate()
        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
        .unwrap_or((0, &0.0));

    let l_max = l_max_val;
    let critical_bus = load_indices.get(max_pos).copied().unwrap_or(0);

    Ok(LIndexResult {
        l_index,
        l_max,
        critical_bus,
        stable: l_max < 1.0,
        stability_margin: (1.0 - l_max).max(0.0),
        load_bus_indices: load_indices,
    })
}

/// Diagonal approximation for L-index when Y_LL is near-singular.
fn compute_l_index_diagonal_fallback(
    _network: &PowerNetwork,
    v_mag: &[f64],
    _v_ang: &[f64],
    load_indices: &[usize],
    gen_indices: &[usize],
) -> Result<LIndexResult> {
    let n_l = load_indices.len();
    // Simple approximation: L_j ≈ 1 - V_j (normalized voltage deviation)
    // Generators are treated as ideal voltage sources at 1.0 pu
    let v_gen_avg = if gen_indices.is_empty() {
        1.0
    } else {
        gen_indices.iter().map(|&i| v_mag[i]).sum::<f64>() / gen_indices.len() as f64
    };

    // Diagonal fallback: L_j ≈ |V_gen/V_j - 1|
    // This preserves the L_j = 0 at no-load property when V_j ≈ V_gen
    let l_index: Vec<f64> = load_indices
        .iter()
        .map(|&i| {
            let vj = v_mag[i];
            if vj < 1e-9 {
                1.0
            } else {
                // At no-load V_j = V_gen → L_j = 0; at collapse V_j < V_gen → L_j > 0
                ((v_gen_avg / vj) - 1.0).abs().min(2.0)
            }
        })
        .collect();

    let (max_pos, &l_max) = l_index
        .iter()
        .enumerate()
        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
        .unwrap_or((0, &0.0));

    let critical_bus = load_indices.get(max_pos).copied().unwrap_or(0);

    // Build actual n_l-sized l_index
    let mut full_l = vec![0.0f64; n_l];
    for (i, &val) in l_index.iter().enumerate() {
        full_l[i] = val;
    }

    Ok(LIndexResult {
        l_index: full_l,
        l_max,
        critical_bus,
        stable: l_max < 1.0,
        stability_margin: (1.0 - l_max).max(0.0),
        load_bus_indices: load_indices.to_vec(),
    })
}

/// Convenience wrapper: compute L-index from a `PowerFlowResult`.
///
/// Builds the complex Y-bus internally and passes voltages from the result.
pub fn compute_l_index_from_result(
    network: &PowerNetwork,
    result: &PowerFlowResult,
) -> Result<LIndexResult> {
    let n = network.buses.len();
    let ybus_sparse = network
        .admittance_matrix()
        .map_err(|e| OxiGridError::InvalidNetwork(format!("Y-bus: {e}")))?;

    // Convert sparse Y-bus to dense Vec<Vec<Complex64>>
    let mut ybus_dense: Vec<Vec<Complex64>> = vec![vec![Complex64::new(0.0, 0.0); n]; n];
    for (&val, (i, j)) in ybus_sparse.iter() {
        if i < n && j < n {
            ybus_dense[i][j] = val;
        }
    }

    compute_l_index(
        network,
        &result.voltage_magnitude,
        &result.voltage_angle,
        &ybus_dense,
    )
}

// ─── Reduced Jacobian Modal Analysis ─────────────────────────────────────────

/// Result of modal voltage stability analysis.
///
/// The minimum singular value of the reduced Jacobian `J_R = J_QV - J_Qθ·J_Pθ⁻¹·J_PV`
/// measures proximity to singularity (voltage collapse).
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ModalVsaResult {
    /// Smallest singular value of `J_R` (→ 0 at collapse).
    pub min_singular_value: f64,
    /// Right singular vector corresponding to the minimum singular value.
    pub critical_mode: Vec<f64>,
    /// Bus participation factors (|v_i|² for each PQ bus in the critical mode).
    pub bus_participation: Vec<f64>,
    /// Branch participation factors (per branch, based on voltage difference).
    pub branch_participation: Vec<f64>,
    /// `true` when `min_singular_value > threshold` (default 0.01).
    pub stable: bool,
}

/// Perform modal voltage stability analysis on a Newton-Raphson Jacobian.
///
/// # Arguments
/// - `jacobian` — the full `(npvpq + npq) × (npvpq + npq)` NR Jacobian as dense rows
/// - `n_pq_buses` — number of PQ buses (determines J_R extraction)
/// - `n_branches` — number of branches in the network (for branch participation sizing)
///
/// # Theory
/// The reduced Jacobian `J_R` is extracted from the 2×2 block structure:
/// ```text
///   J = [J_Pθ  J_PV]
///       [J_Qθ  J_QV]
/// ```
/// `J_R = J_QV - J_Qθ · J_Pθ⁻¹ · J_PV`
///
/// The minimum singular value is found via inverse power iteration.
pub fn modal_voltage_stability(
    jacobian: &[Vec<f64>],
    n_pq_buses: usize,
    branches_pq_indices: &[(usize, usize)],
) -> Result<ModalVsaResult> {
    if jacobian.is_empty() || n_pq_buses == 0 {
        return Err(OxiGridError::InvalidParameter(
            "empty Jacobian or no PQ buses".into(),
        ));
    }

    let n_branches = branches_pq_indices.len();
    let j_size = jacobian.len();
    let npq = n_pq_buses;

    // The NR Jacobian has dimension (npvpq + npq) × (npvpq + npq)
    // where npvpq = total non-slack buses
    // Blocks (0-indexed):
    //   J_Pt (rows 0..npvpq, cols 0..npvpq)  = ∂P/∂θ
    //   J_PV (rows 0..npvpq, cols npvpq..end) = ∂P/∂|V|
    //   J_Qt (rows npvpq..end, cols 0..npvpq) = ∂Q/∂θ
    //   J_QV (rows npvpq..end, cols npvpq..end) = ∂Q/∂|V|
    //
    // npvpq = j_size - npq
    let npvpq = j_size.saturating_sub(npq);
    if npvpq == 0 || npvpq > j_size {
        return Err(OxiGridError::InvalidParameter(
            "n_pq_buses exceeds Jacobian dimension".into(),
        ));
    }

    // Extract sub-blocks
    // J_Ptheta: npvpq × npvpq
    // J_PV:     npvpq × npq
    // J_Qtheta: npq × npvpq
    // J_QV:     npq × npq
    let mut j_pt = vec![vec![0.0f64; npvpq]; npvpq];
    let mut j_pv = vec![vec![0.0f64; npq]; npvpq];
    let mut j_qt = vec![vec![0.0f64; npvpq]; npq];
    let mut j_qv = vec![vec![0.0f64; npq]; npq];

    for i in 0..npvpq {
        if i >= jacobian.len() {
            continue;
        }
        for j in 0..npvpq {
            if j < jacobian[i].len() {
                j_pt[i][j] = jacobian[i][j];
            }
        }
        #[allow(clippy::needless_range_loop)]
        for j in 0..npq {
            let col = npvpq + j;
            if col < jacobian[i].len() {
                j_pv[i][j] = jacobian[i][col];
            }
        }
    }
    for i in 0..npq {
        let row = npvpq + i;
        if row >= jacobian.len() {
            continue;
        }
        for j in 0..npvpq {
            if j < jacobian[row].len() {
                j_qt[i][j] = jacobian[row][j];
            }
        }
        #[allow(clippy::needless_range_loop)]
        for j in 0..npq {
            let col = npvpq + j;
            if col < jacobian[row].len() {
                j_qv[i][j] = jacobian[row][col];
            }
        }
    }

    // Compute J_R = J_QV - J_Qt * J_Pt⁻¹ * J_PV
    // Step 1: Solve J_Pt * X = J_PV  (X = J_Pt⁻¹ · J_PV), shape npvpq × npq
    let x_mat = solve_linear_system_multi_rhs(&j_pt, &j_pv)?;

    // Step 2: J_Qt * X → shape npq × npq
    let j_qt_x = mat_mul(&j_qt, &x_mat, npq, npvpq, npq);

    // Step 3: J_R = J_QV - J_Qt_X
    let mut j_r = vec![vec![0.0f64; npq]; npq];
    for i in 0..npq {
        for j in 0..npq {
            j_r[i][j] = j_qv[i][j] - j_qt_x[i][j];
        }
    }

    // Compute minimum singular value via inverse power iteration
    let (min_sv, right_sv) = min_singular_value_power(&j_r, 200, 1e-8)?;

    // Bus participation factors: |v_i|² normalized
    let sv_norm_sq: f64 = right_sv.iter().map(|&x| x * x).sum();
    let bus_participation: Vec<f64> = if sv_norm_sq > 1e-14 {
        right_sv.iter().map(|&x| x * x / sv_norm_sq).collect()
    } else {
        vec![1.0 / npq as f64; npq]
    };

    // Branch participation: proportional to squared voltage-magnitude differences
    // across branch endpoints in the PQ-bus subvector.
    let branch_participation: Vec<f64> = if branches_pq_indices.is_empty() {
        Vec::new()
    } else {
        let diffs: Vec<f64> = branches_pq_indices
            .iter()
            .map(|&(fi, ti)| {
                let vf = right_sv.get(fi).copied().unwrap_or(0.0);
                let vt = right_sv.get(ti).copied().unwrap_or(0.0);
                (vf - vt).powi(2)
            })
            .collect();
        let total: f64 = diffs.iter().sum();
        if total > 1e-14 {
            diffs.iter().map(|&d| d / total).collect()
        } else {
            vec![1.0 / n_branches as f64; n_branches]
        }
    };

    Ok(ModalVsaResult {
        min_singular_value: min_sv,
        critical_mode: right_sv,
        bus_participation,
        branch_participation,
        stable: min_sv > 0.01,
    })
}

/// Find the minimum singular value of a matrix via inverse power iteration on `A^T A`.
///
/// Uses the identity: min singular value of `A` = sqrt(min eigenvalue of `A^T A`).
/// Inverse power iteration converges to the smallest eigenvalue of `A^T A`.
fn min_singular_value_power(
    matrix: &[Vec<f64>],
    max_iter: usize,
    tol: f64,
) -> Result<(f64, Vec<f64>)> {
    let m = matrix.len();
    if m == 0 {
        return Ok((0.0, vec![]));
    }
    let n = matrix[0].len();
    if n == 0 {
        return Ok((0.0, vec![]));
    }

    // Build A^T A  (n × n)
    let mut ata = vec![vec![0.0f64; n]; n];
    #[allow(clippy::needless_range_loop)]
    for i in 0..n {
        for j in 0..n {
            let mut s = 0.0f64;
            #[allow(clippy::needless_range_loop)]
            for k in 0..m {
                if i < matrix[k].len() && j < matrix[k].len() {
                    s += matrix[k][i] * matrix[k][j];
                }
            }
            ata[i][j] = s;
        }
    }

    // Regularize slightly to ensure invertibility for inverse iteration
    let regularization = 1e-12;
    #[allow(clippy::needless_range_loop)]
    for i in 0..n {
        ata[i][i] += regularization;
    }

    // Initial vector (random-like using index-based initialization)
    let mut x: Vec<f64> = (0..n).map(|i| ((i + 1) as f64).recip()).collect();
    let x_norm: f64 = x.iter().map(|v| v * v).sum::<f64>().sqrt();
    if x_norm > 1e-14 {
        for xi in &mut x {
            *xi /= x_norm;
        }
    }

    let mut eigenvalue = 0.0f64;
    let mut prev_eigenvalue = f64::INFINITY;

    for _iter in 0..max_iter {
        // Solve (A^T A) y = x  (inverse power iteration step)
        let y = solve_linear_system(&ata, &x)?;

        // Rayleigh quotient: λ ≈ x^T (A^T A) x / x^T x = 1 / (x^T y / ||y||²)
        let xy: f64 = x.iter().zip(y.iter()).map(|(xi, yi)| xi * yi).sum();
        let yy: f64 = y.iter().map(|yi| yi * yi).sum();

        eigenvalue = if xy.abs() > 1e-14 { yy / xy } else { 0.0 };

        // Normalize y to get new x
        let y_norm: f64 = yy.sqrt();
        if y_norm < 1e-14 {
            break;
        }
        x = y.iter().map(|yi| yi / y_norm).collect();

        if (eigenvalue - prev_eigenvalue).abs() < tol * eigenvalue.abs().max(1e-10) {
            break;
        }
        prev_eigenvalue = eigenvalue;
    }

    // min singular value = sqrt(min eigenvalue of A^T A)
    let min_sv = eigenvalue.max(0.0).sqrt();
    Ok((min_sv, x))
}

// ─── Voltage Stability Margin (Continuation-based) ────────────────────────────

/// Method used to compute the voltage stability margin.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub enum VsaMethod {
    /// Continuation power flow (CPF) – most accurate.
    ContinuationPowerFlow,
    /// Direct point-of-collapse (POC) method.
    PointOfCollapse,
    /// Sensitivity-based fast screening (approximate).
    FastVoltageStability,
}

/// Voltage stability margin from operating point to voltage collapse.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct VoltageStabilityMargin {
    /// Maximum loadability (p.u. loading) where power flow converges.
    pub lambda_max: f64,
    /// Bus voltage magnitudes at the nose point (collapse).
    pub voltage_at_collapse: Vec<f64>,
    /// Bus indices where voltage drops below 0.8 p.u. at collapse.
    pub critical_buses: Vec<usize>,
    /// MW headroom = (λ_max − 1) × total base load \[MW\].
    pub loading_margin_mw: f64,
    /// Method used to estimate the margin.
    pub method: VsaMethod,
}

/// Compute the voltage stability margin by bisection on convergence.
///
/// Incrementally loads the system until power flow diverges, then bisects
/// to find λ_max with the requested tolerance.
pub fn compute_voltage_stability_margin(
    network: &PowerNetwork,
    method: VsaMethod,
    lambda_step: f64,
    tol: f64,
) -> Result<VoltageStabilityMargin> {
    let config = PowerFlowConfig {
        method: PowerFlowMethod::NewtonRaphson,
        max_iter: 100,
        tolerance: 1e-6,
        enforce_q_limits: false,
    };

    let base_load_mw = network.total_load_mw();

    match method {
        VsaMethod::FastVoltageStability => {
            // Fast sensitivity-based approximation:
            // Estimate λ_max from L-index headroom
            let base_result = network.solve_powerflow(&config).map_err(|e| {
                OxiGridError::InvalidNetwork(format!("base power flow failed: {e}"))
            })?;

            let l_idx = compute_l_index_from_result(network, &base_result)?;
            // Approximate loading margin from L-index stability margin
            let lambda_max = if l_idx.stability_margin > 1e-6 {
                l_idx.stability_margin * 2.0 // heuristic scaling
            } else {
                0.01
            };

            let v_collapse = base_result
                .voltage_magnitude
                .iter()
                .map(|&v| (v * 0.8).max(0.5))
                .collect::<Vec<_>>();

            let critical_buses = v_collapse
                .iter()
                .enumerate()
                .filter(|(_, &v)| v < 0.8)
                .map(|(i, _)| i)
                .collect();

            Ok(VoltageStabilityMargin {
                lambda_max,
                voltage_at_collapse: v_collapse,
                critical_buses,
                loading_margin_mw: (lambda_max - 1.0).max(0.0) * base_load_mw,
                method,
            })
        }
        _ => {
            // Bisection on convergence for ContinuationPowerFlow and PointOfCollapse
            let mut lo = 0.0f64;
            let mut hi = lambda_step;
            let mut v_at_hi: Vec<f64> = network.buses.iter().map(|_| 1.0).collect();

            // Find first divergence
            loop {
                let scaled = scale_network_lambda(network, hi);
                match scaled.solve_powerflow(&config) {
                    Ok(r) if r.converged => {
                        v_at_hi = r.voltage_magnitude;
                        lo = hi;
                        hi *= 2.0;
                        if hi > 20.0 {
                            break;
                        }
                    }
                    _ => break,
                }
            }

            // Bisect to find λ_max within tolerance
            for _ in 0..60 {
                if hi - lo < tol {
                    break;
                }
                let mid = 0.5 * (lo + hi);
                let scaled = scale_network_lambda(network, mid);
                match scaled.solve_powerflow(&config) {
                    Ok(r) if r.converged => {
                        v_at_hi = r.voltage_magnitude;
                        lo = mid;
                    }
                    _ => {
                        hi = mid;
                    }
                }
            }

            let lambda_max = lo;
            let critical_buses: Vec<usize> = v_at_hi
                .iter()
                .enumerate()
                .filter(|(_, &v)| v < 0.8)
                .map(|(i, _)| i)
                .collect();

            Ok(VoltageStabilityMargin {
                lambda_max,
                voltage_at_collapse: v_at_hi,
                critical_buses,
                loading_margin_mw: (lambda_max - 1.0).max(0.0) * base_load_mw,
                method,
            })
        }
    }
}

// ─── Fast Voltage Stability Index (FVSI) ─────────────────────────────────────

/// Compute the Fast Voltage Stability Index (FVSI) for all branches.
///
/// For branch i→j:
/// ```text
///   FVSI_ij = 4 · Z² · Q_j / (V_i² · X)
/// ```
/// where `Z = |R + jX|`, `X = Im(Z)`, `Q_j` = reactive power at receiving end.
///
/// FVSI → 0: stable; FVSI → 1: voltage instability on the branch.
pub fn compute_fvsi(network: &PowerNetwork, v_mag: &[f64], v_ang: &[f64]) -> Result<Vec<f64>> {
    let n_branches = network.branches.len();
    let n_buses = network.buses.len();

    if v_mag.len() < n_buses || v_ang.len() < n_buses {
        return Err(OxiGridError::InvalidParameter(
            "v_mag/v_ang too short for bus count".into(),
        ));
    }

    let mut fvsi = Vec::with_capacity(n_branches);

    for branch in &network.branches {
        let from_idx = network
            .buses
            .iter()
            .position(|b| b.id == branch.from_bus)
            .ok_or_else(|| {
                OxiGridError::InvalidNetwork(format!("Bus {} not found", branch.from_bus))
            })?;
        let to_idx = network
            .buses
            .iter()
            .position(|b| b.id == branch.to_bus)
            .ok_or_else(|| {
                OxiGridError::InvalidNetwork(format!("Bus {} not found", branch.to_bus))
            })?;

        if !branch.status {
            fvsi.push(0.0);
            continue;
        }

        let r = branch.r;
        let x = branch.x;
        let z_sq = r * r + x * x;

        // Skip degenerate branches
        if z_sq < 1e-18 || x.abs() < 1e-12 {
            fvsi.push(0.0);
            continue;
        }

        let vi = v_mag[from_idx];
        if vi < 1e-9 {
            fvsi.push(0.0);
            continue;
        }

        // Compute reactive power at receiving (to) bus using π-model
        // Q_j = Im(V_j * I_j*) where current is estimated from voltage difference
        let vi_c = Complex64::new(vi * v_ang[from_idx].cos(), vi * v_ang[from_idx].sin());
        let vj_c = Complex64::new(
            v_mag[to_idx] * v_ang[to_idx].cos(),
            v_mag[to_idx] * v_ang[to_idx].sin(),
        );

        let tap = branch.effective_tap();
        let z_series = Complex64::new(r, x);
        let y_series = Complex64::new(1.0, 0.0) / z_series;

        // Current flowing into to-bus (receiving end)
        let i_to = y_series * (vj_c - vi_c / tap);
        let s_to = vj_c * i_to.conj();
        let q_j = s_to.im.abs() * network.base_mva;

        // FVSI = 4 * Z² * |Q_j| / (V_i² * X)
        let index = 4.0 * z_sq * q_j / (vi * vi * x.abs() * network.base_mva);
        fvsi.push(index.min(10.0)); // cap to prevent overflow
    }

    Ok(fvsi)
}

// ─── Line Stability Index (Lmn / PSI) ────────────────────────────────────────

/// Compute the Line Stability Index (Lmn) for all branches.
///
/// Per-line index based on the real power transfer:
/// ```text
///   Lmn = 4·P·Z² / (V_s²·(Z·cos(θ−δ))²)
/// ```
/// where `P` is the real power at the receiving end, `θ = ∠Z`, `δ = δ_s − δ_r`.
///
/// Voltage collapse on a line when `Lmn ≥ 1`.
pub fn compute_line_stability_index(
    network: &PowerNetwork,
    v_mag: &[f64],
    v_ang: &[f64],
) -> Result<Vec<f64>> {
    let n_branches = network.branches.len();
    let n_buses = network.buses.len();

    if v_mag.len() < n_buses || v_ang.len() < n_buses {
        return Err(OxiGridError::InvalidParameter(
            "v_mag/v_ang too short for bus count".into(),
        ));
    }

    let mut lmn = Vec::with_capacity(n_branches);

    for branch in &network.branches {
        let from_idx = network
            .buses
            .iter()
            .position(|b| b.id == branch.from_bus)
            .ok_or_else(|| {
                OxiGridError::InvalidNetwork(format!("Bus {} not found", branch.from_bus))
            })?;
        let to_idx = network
            .buses
            .iter()
            .position(|b| b.id == branch.to_bus)
            .ok_or_else(|| {
                OxiGridError::InvalidNetwork(format!("Bus {} not found", branch.to_bus))
            })?;

        if !branch.status {
            lmn.push(0.0);
            continue;
        }

        let r = branch.r;
        let x = branch.x;
        let z_sq = r * r + x * x;
        let z_mag = z_sq.sqrt();

        if z_sq < 1e-18 {
            lmn.push(0.0);
            continue;
        }

        let vs = v_mag[from_idx];
        if vs < 1e-9 {
            lmn.push(0.0);
            continue;
        }

        // Angle of impedance: θ_z = atan2(x, r)
        let theta_z = x.atan2(r);

        // Power angle difference: δ = δ_s - δ_r
        let delta = v_ang[from_idx] - v_ang[to_idx];

        // Receiving-end power from π-model
        let vi_c = Complex64::new(vs * v_ang[from_idx].cos(), vs * v_ang[from_idx].sin());
        let vj_c = Complex64::new(
            v_mag[to_idx] * v_ang[to_idx].cos(),
            v_mag[to_idx] * v_ang[to_idx].sin(),
        );

        let tap = branch.effective_tap();
        let z_series = Complex64::new(r, x);
        let y_series = Complex64::new(1.0, 0.0) / z_series;
        let i_to = y_series * (vj_c - vi_c / tap);
        let s_to = vj_c * i_to.conj();
        let p_r = s_to.re.abs() * network.base_mva;

        // Denominator: (Z * cos(θ_z - δ))²
        let cos_term = (theta_z - delta).cos();
        let denom = vs * vs * (z_mag * cos_term) * (z_mag * cos_term);

        let index = if denom > 1e-14 {
            4.0 * p_r * z_sq / denom / network.base_mva
        } else {
            10.0 // collapse indicator
        };

        lmn.push(index.min(10.0));
    }

    Ok(lmn)
}

// ─── N-1 Voltage Stability Screening ─────────────────────────────────────────

/// Configuration for N-1 voltage stability contingency screening.
pub struct VoltageStabilityScreener {
    /// Branch indices (0-based) to remove for contingency analysis.
    pub contingencies: Vec<usize>,
    /// Per-unit load increment to apply during fast check.
    pub loading_increment: f64,
    /// L-index alarm threshold (default 0.8).
    pub threshold: f64,
}

impl Default for VoltageStabilityScreener {
    fn default() -> Self {
        Self {
            contingencies: vec![],
            loading_increment: 0.05,
            threshold: 0.8,
        }
    }
}

/// Result for a single contingency in N-1 voltage stability screening.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ContingencyVsaResult {
    /// Branch index that was removed (0-based).
    pub contingency_branch: usize,
    /// Pre-contingency system L-index.
    pub l_index_pre: f64,
    /// Post-contingency system L-index.
    pub l_index_post: f64,
    /// Maximum FVSI across all branches post-contingency.
    pub fvsi_max_post: f64,
    /// `true` if any bus voltage < 0.95 p.u. post-contingency.
    pub voltage_violation: bool,
    /// Severity rank (1 = most severe, assigned after sorting).
    pub rank: usize,
}

impl VoltageStabilityScreener {
    /// Run N-1 screening.
    ///
    /// For each contingency: remove branch, run power flow, compute L-index and FVSI.
    /// Results are ranked by post-contingency L-index (worst first).
    pub fn screen(
        &self,
        network: &PowerNetwork,
        base_flow_result: &PowerFlowResult,
    ) -> Result<Vec<ContingencyVsaResult>> {
        let config = PowerFlowConfig {
            method: PowerFlowMethod::NewtonRaphson,
            max_iter: 100,
            tolerance: 1e-6,
            enforce_q_limits: false,
        };

        // Compute pre-contingency L-index
        let l_pre = compute_l_index_from_result(network, base_flow_result)
            .map(|r| r.l_max)
            .unwrap_or(0.0);

        let contingency_list: Vec<usize> = if self.contingencies.is_empty() {
            (0..network.branches.len()).collect()
        } else {
            self.contingencies.clone()
        };

        let mut results: Vec<ContingencyVsaResult> = Vec::new();

        for &branch_idx in &contingency_list {
            if branch_idx >= network.branches.len() {
                continue;
            }

            // Create post-contingency network (branch removed)
            let mut post_net = network.clone();
            post_net.branches[branch_idx].status = false;

            // Run post-contingency power flow
            let post_result = match post_net.solve_powerflow(&config) {
                Ok(r) if r.converged => r,
                _ => {
                    // Diverged: extreme case
                    results.push(ContingencyVsaResult {
                        contingency_branch: branch_idx,
                        l_index_pre: l_pre,
                        l_index_post: 1.0,
                        fvsi_max_post: 10.0,
                        voltage_violation: true,
                        rank: 0,
                    });
                    continue;
                }
            };

            // Compute post L-index
            let l_post = compute_l_index_from_result(&post_net, &post_result)
                .map(|r| r.l_max)
                .unwrap_or(0.0);

            // Compute post FVSI
            let fvsi_max = compute_fvsi(
                &post_net,
                &post_result.voltage_magnitude,
                &post_result.voltage_angle,
            )
            .map(|v| v.iter().cloned().fold(0.0f64, f64::max))
            .unwrap_or(0.0);

            // Check voltage violations (< 0.95 p.u.)
            let voltage_violation = post_result.voltage_magnitude.iter().any(|&v| v < 0.95);

            results.push(ContingencyVsaResult {
                contingency_branch: branch_idx,
                l_index_pre: l_pre,
                l_index_post: l_post,
                fvsi_max_post: fvsi_max,
                voltage_violation,
                rank: 0,
            });
        }

        // Sort by post-contingency L-index descending (worst first)
        results.sort_by(|a, b| {
            b.l_index_post
                .partial_cmp(&a.l_index_post)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        // Assign ranks
        for (rank, result) in results.iter_mut().enumerate() {
            result.rank = rank + 1;
        }

        Ok(results)
    }
}

// ─── Linear algebra helpers ───────────────────────────────────────────────────

/// Solve Ax = b with Gaussian elimination and partial pivoting.
fn solve_linear_system(a: &[Vec<f64>], b: &[f64]) -> Result<Vec<f64>> {
    let n = a.len();
    if n == 0 {
        return Ok(vec![]);
    }
    if b.len() != n {
        return Err(OxiGridError::LinearAlgebra(
            "dimension mismatch in solve_linear_system".into(),
        ));
    }

    // Build augmented matrix [A | b]
    let mut aug: Vec<Vec<f64>> = a
        .iter()
        .zip(b.iter())
        .map(|(row, &rhs)| {
            let mut r = row.clone();
            r.push(rhs);
            r
        })
        .collect();

    for col in 0..n {
        // Partial pivoting
        let pivot_row = (col..n)
            .max_by(|&r1, &r2| {
                aug[r1][col]
                    .abs()
                    .partial_cmp(&aug[r2][col].abs())
                    .unwrap_or(std::cmp::Ordering::Equal)
            })
            .unwrap_or(col);

        if aug[pivot_row][col].abs() < 1e-14 {
            return Err(OxiGridError::LinearAlgebra(
                "singular matrix in solve_linear_system".into(),
            ));
        }

        aug.swap(col, pivot_row);

        let pivot = aug[col][col];
        #[allow(clippy::needless_range_loop)]
        for j in col..=n {
            aug[col][j] /= pivot;
        }

        for row in 0..n {
            if row == col {
                continue;
            }
            let factor = aug[row][col];
            #[allow(clippy::needless_range_loop)]
            for j in col..=n {
                aug[row][j] -= factor * aug[col][j];
            }
        }
    }

    Ok(aug.iter().map(|row| row[n]).collect())
}

/// Solve A·X = B for multiple right-hand sides simultaneously.
///
/// Returns X as a matrix (rows = unknowns, cols = RHS columns).
fn solve_linear_system_multi_rhs(a: &[Vec<f64>], b: &[Vec<f64>]) -> Result<Vec<Vec<f64>>> {
    let n = a.len();
    let m = b.first().map(|r| r.len()).unwrap_or(0);

    if n == 0 || m == 0 {
        return Ok(vec![vec![0.0f64; m]; n]);
    }

    if b.len() != n {
        return Err(OxiGridError::LinearAlgebra(
            "RHS row count mismatch in solve_linear_system_multi_rhs".into(),
        ));
    }

    // Build augmented matrix [A | B]
    let mut aug: Vec<Vec<f64>> = a
        .iter()
        .zip(b.iter())
        .map(|(row_a, row_b)| {
            let mut r = row_a.clone();
            r.extend_from_slice(row_b);
            r
        })
        .collect();

    for col in 0..n {
        // Partial pivoting
        let pivot_row = (col..n)
            .max_by(|&r1, &r2| {
                aug[r1][col]
                    .abs()
                    .partial_cmp(&aug[r2][col].abs())
                    .unwrap_or(std::cmp::Ordering::Equal)
            })
            .unwrap_or(col);

        if aug[pivot_row][col].abs() < 1e-14 {
            // Near-singular: return zeros for this block
            return Ok(vec![vec![0.0f64; m]; n]);
        }

        aug.swap(col, pivot_row);

        let pivot = aug[col][col];
        #[allow(clippy::needless_range_loop)]
        for j in col..n + m {
            aug[col][j] /= pivot;
        }

        for row in 0..n {
            if row == col {
                continue;
            }
            let factor = aug[row][col];
            #[allow(clippy::needless_range_loop)]
            for j in col..n + m {
                aug[row][j] -= factor * aug[col][j];
            }
        }
    }

    // Extract solution columns
    let mut x = vec![vec![0.0f64; m]; n];
    for i in 0..n {
        for j in 0..m {
            x[i][j] = aug[i][n + j];
        }
    }
    Ok(x)
}

/// Dense matrix multiplication: C = A · B, where A is (m × k) and B is (k × n).
fn mat_mul(
    a: &[Vec<f64>],
    b: &[Vec<f64>],
    rows_a: usize,
    cols_a: usize,
    cols_b: usize,
) -> Vec<Vec<f64>> {
    let mut c = vec![vec![0.0f64; cols_b]; rows_a];
    for i in 0..rows_a {
        if i >= a.len() {
            continue;
        }
        for k in 0..cols_a {
            if k >= a[i].len() || k >= b.len() {
                continue;
            }
            let a_ik = a[i][k];
            for j in 0..cols_b {
                if j < b[k].len() {
                    c[i][j] += a_ik * b[k][j];
                }
            }
        }
    }
    c
}

/// Scale all loads and non-slack generators by factor (1 + lambda).
fn scale_network_lambda(net: &PowerNetwork, lambda: f64) -> PowerNetwork {
    let mut scaled = net.clone();
    let factor = 1.0 + lambda;
    for bus in &mut scaled.buses {
        bus.pd.0 *= factor;
        bus.qd.0 *= factor;
    }
    if let Ok(slack_idx) = scaled.slack_bus_index() {
        for gen in &mut scaled.generators {
            if gen.bus_id != slack_idx {
                gen.pg *= factor;
            }
        }
    }
    scaled
}

// ─── Tests ────────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;
    use crate::network::topology::Generator;
    use crate::network::{Branch, Bus, BusType, PowerNetwork};
    use crate::powerflow::PowerFlowConfig;

    fn load_ieee14() -> PowerNetwork {
        let path = concat!(env!("CARGO_MANIFEST_DIR"), "/tests/data/ieee14.m");
        PowerNetwork::from_matpower(path).expect("ieee14 parse")
    }

    /// Build a minimal 2-bus network (infinite bus + load bus) for formula verification.
    fn two_bus_network() -> PowerNetwork {
        let mut net = PowerNetwork::new(100.0);

        // Bus 1: slack (infinite bus)
        let mut b1 = Bus::new(1, BusType::Slack);
        b1.vm = 1.05;
        b1.va = 0.0;
        net.buses.push(b1);

        // Bus 2: PQ load bus
        let mut b2 = Bus::new(2, BusType::PQ);
        b2.vm = 1.0;
        b2.va = 0.0;
        b2.pd = crate::units::Power(50.0);
        b2.qd = crate::units::ReactivePower(20.0);
        net.buses.push(b2);

        // Branch 1-2: r=0.01, x=0.1, b=0
        net.branches.push(Branch {
            from_bus: 1,
            to_bus: 2,
            r: 0.01,
            x: 0.1,
            b: 0.0,
            rate_a: 0.0,
            rate_b: 0.0,
            rate_c: 0.0,
            tap: 0.0,
            shift: 0.0,
            status: true,
        });

        // Slack generator
        net.generators.push(Generator {
            bus_id: 1,
            pg: 50.0,
            qg: 0.0,
            qmax: 200.0,
            qmin: -200.0,
            vg: 1.05,
            mbase: 100.0,
            status: true,
            pmax: 500.0,
            pmin: 0.0,
        });

        net
    }

    #[test]
    fn test_l_index_base_case() {
        let net = load_ieee14();
        let config = PowerFlowConfig::default();
        let result = net.solve_powerflow(&config).expect("power flow");
        let li = compute_l_index_from_result(&net, &result).expect("L-index");
        // IEEE 14-bus at base loading: well-loaded but stable
        assert!(
            li.l_max >= 0.0,
            "L-max should be non-negative: {}",
            li.l_max
        );
        assert!(
            li.l_max < 0.5,
            "L-max should be < 0.5 at base loading: {:.4}",
            li.l_max
        );
        assert!(li.stable, "System should be stable at base loading");
        assert!(
            li.stability_margin > 0.5,
            "Margin should be > 0.5: {:.4}",
            li.stability_margin
        );
    }

    #[test]
    fn test_l_index_at_collapse() {
        let net = load_ieee14();
        // Scale loads to heavy loading (3x base)
        let mut heavy = net.clone();
        for bus in &mut heavy.buses {
            bus.pd.0 *= 3.0;
            bus.qd.0 *= 3.0;
        }
        let config = PowerFlowConfig {
            method: PowerFlowMethod::NewtonRaphson,
            max_iter: 100,
            tolerance: 1e-6,
            enforce_q_limits: false,
        };
        // Power flow may or may not converge under heavy loading
        if let Ok(result) = heavy.solve_powerflow(&config) {
            if result.converged {
                let li = compute_l_index_from_result(&heavy, &result).expect("L-index heavy");
                // Expect higher L-index under heavy load
                assert!(li.l_max > 0.0, "L-max should be positive under heavy load");
            }
        }
        // Test passes regardless (heavy loading may cause divergence)
    }

    #[test]
    fn test_fvsi_all_positive() {
        let net = load_ieee14();
        let config = PowerFlowConfig::default();
        let result = net.solve_powerflow(&config).expect("power flow");
        let fvsi =
            compute_fvsi(&net, &result.voltage_magnitude, &result.voltage_angle).expect("FVSI");

        for (i, &f) in fvsi.iter().enumerate() {
            assert!(f >= 0.0, "FVSI[{}] = {:.6} should be non-negative", i, f);
        }
    }

    #[test]
    fn test_fvsi_increases_with_load() {
        let net = load_ieee14();
        let config = PowerFlowConfig {
            method: PowerFlowMethod::NewtonRaphson,
            max_iter: 100,
            tolerance: 1e-6,
            enforce_q_limits: false,
        };

        // Base case
        let base_result = net.solve_powerflow(&config).expect("base PF");
        let fvsi_base = compute_fvsi(
            &net,
            &base_result.voltage_magnitude,
            &base_result.voltage_angle,
        )
        .expect("FVSI base");

        // Increased load (1.5x)
        let mut heavy = net.clone();
        for bus in &mut heavy.buses {
            bus.pd.0 *= 1.5;
            bus.qd.0 *= 1.5;
        }
        if let Ok(heavy_result) = heavy.solve_powerflow(&config) {
            if heavy_result.converged {
                let fvsi_heavy = compute_fvsi(
                    &heavy,
                    &heavy_result.voltage_magnitude,
                    &heavy_result.voltage_angle,
                )
                .expect("FVSI heavy");

                // Sum of FVSI should increase under heavier load
                let sum_base: f64 = fvsi_base.iter().sum();
                let sum_heavy: f64 = fvsi_heavy.iter().sum();
                assert!(
                    sum_heavy >= sum_base,
                    "FVSI sum should increase with load: base={:.4}, heavy={:.4}",
                    sum_base,
                    sum_heavy
                );
            }
        }
    }

    #[test]
    fn test_modal_vsa_positive_sv() {
        let net = load_ieee14();
        let config = PowerFlowConfig::default();
        let result = net.solve_powerflow(&config).expect("power flow");

        // Build a simple diagonal Jacobian proxy for testing
        // (3x3 PQ-submatrix identity-like)
        let npq = net
            .buses
            .iter()
            .filter(|b| b.bus_type == BusType::PQ)
            .count();
        let npvpq = net
            .buses
            .iter()
            .filter(|b| b.bus_type != BusType::Slack)
            .count();
        let j_size = npvpq + npq;

        // Build approximate Jacobian: diagonal-dominant matrix
        let mut jac = vec![vec![0.0f64; j_size]; j_size];
        let v = &result.voltage_magnitude;
        for i in 0..j_size.min(v.len()) {
            jac[i][i] = v[i] * 2.0 + 0.1;
            if i + 1 < j_size {
                jac[i][i + 1] = -0.05;
                jac[i + 1][i] = -0.05;
            }
        }

        // Map branch endpoints to PQ-bus subvector indices.
        let pq_ids: Vec<usize> = net
            .buses
            .iter()
            .filter(|b| b.bus_type == BusType::PQ)
            .map(|b| b.id)
            .collect();
        let pq_idx = |id: usize| pq_ids.iter().position(|&p| p == id);
        let branches_pq: Vec<(usize, usize)> = net
            .branches
            .iter()
            .filter_map(|br| Some((pq_idx(br.from_bus)?, pq_idx(br.to_bus)?)))
            .collect();

        let modal = modal_voltage_stability(&jac, npq, &branches_pq).expect("modal VSA");

        assert!(
            modal.min_singular_value > 0.0,
            "Min singular value should be positive: {:.6}",
            modal.min_singular_value
        );
        assert!(modal.stable, "Stable system should have stable=true");
    }

    #[test]
    fn test_branch_participation_normalized() {
        // 3-bus fully meshed (bus 0 slack, buses 1 and 2 PQ, branches 0-1, 0-2, 1-2)
        let jacobian = vec![
            vec![4.0, -2.0, -2.0, 0.0],
            vec![-2.0, 4.0, 0.0, -2.0],
            vec![-2.0, 0.0, 4.0, 0.0],
            vec![0.0, -2.0, 0.0, 4.0],
        ];
        // Branches between PQ bus indices: (0,1), (0,1) double, (0,0) self-loop as edge case
        let branches = vec![(0_usize, 1_usize), (1_usize, 0_usize)];
        let result = modal_voltage_stability(&jacobian, 2, &branches)
            .expect("modal VSA branch participation");
        let sum: f64 = result.branch_participation.iter().sum();
        assert!(
            (sum - 1.0).abs() < 1e-10 || result.branch_participation.is_empty(),
            "branch participation must sum to 1 or be empty; got {}",
            sum
        );
        assert!(
            !result.branch_participation.iter().all(|&p| p == 0.0),
            "at least one branch participation factor must be non-zero"
        );
    }

    #[test]
    fn test_line_stability_index_range() {
        let net = load_ieee14();
        let config = PowerFlowConfig::default();
        let result = net.solve_powerflow(&config).expect("power flow");
        let lmn =
            compute_line_stability_index(&net, &result.voltage_magnitude, &result.voltage_angle)
                .expect("Lmn");

        // Lmn should be non-negative and finite (capped at 10.0) for all branches.
        // Note: Lmn can exceed 1.0 for realistic loading conditions due to the
        // angle-sensitive formula. The index identifies the RELATIVE ranking of
        // line stress, not an absolute [0,1] bound.
        for (i, &val) in lmn.iter().enumerate() {
            assert!(val >= 0.0, "Lmn[{}] = {:.6} should be non-negative", i, val);
            assert!(val.is_finite(), "Lmn[{}] = {:.6} should be finite", i, val);
            assert!(
                val <= 10.0 + 1e-6,
                "Lmn[{}] = {:.6} should be bounded by the cap of 10.0",
                i,
                val
            );
        }
        // Lmn should vary across branches (not all zero)
        let max_lmn = lmn.iter().cloned().fold(0.0f64, f64::max);
        assert!(max_lmn > 0.0, "Max Lmn should be positive: {:.4}", max_lmn);
    }

    #[test]
    fn test_n1_screening_ranks_by_severity() {
        let net = load_ieee14();
        let config = PowerFlowConfig::default();
        let base_result = net.solve_powerflow(&config).expect("power flow");

        let screener = VoltageStabilityScreener {
            contingencies: (0..5).collect(), // test first 5 branches
            loading_increment: 0.05,
            threshold: 0.8,
        };

        let results = screener.screen(&net, &base_result).expect("N-1 screening");

        // Results should be sorted by l_index_post descending
        for window in results.windows(2) {
            assert!(
                window[0].l_index_post >= window[1].l_index_post,
                "Results not sorted: {:.4} < {:.4}",
                window[0].l_index_post,
                window[1].l_index_post
            );
        }

        // Rank should be 1-based and ascending
        for (i, r) in results.iter().enumerate() {
            assert_eq!(r.rank, i + 1, "Rank mismatch at position {}", i);
        }
    }

    #[test]
    fn test_l_index_balanced_network() {
        // 2-bus infinite bus system: verify formula implementation
        let net = two_bus_network();
        let config = PowerFlowConfig {
            method: PowerFlowMethod::NewtonRaphson,
            max_iter: 100,
            tolerance: 1e-6,
            enforce_q_limits: false,
        };

        let result = net.solve_powerflow(&config).expect("2-bus power flow");
        assert!(result.converged, "2-bus PF should converge");

        let li = compute_l_index_from_result(&net, &result).expect("L-index 2-bus");

        // Only 1 load bus (bus 2), so l_index has length 1
        assert_eq!(
            li.l_index.len(),
            1,
            "Should have exactly 1 load bus L-index"
        );
        // L_j should be in [0, 1)
        let lj = li.l_index[0];
        assert!(
            (0.0..1.0).contains(&lj),
            "L_j = {:.4} should be in [0, 1)",
            lj
        );
        // For a lightly loaded 2-bus system, L should be small
        assert!(
            lj < 0.5,
            "L_j = {:.4} should be < 0.5 for lightly loaded 2-bus",
            lj
        );
    }

    #[test]
    fn test_voltage_stability_margin_fast() {
        let net = load_ieee14();
        let margin =
            compute_voltage_stability_margin(&net, VsaMethod::FastVoltageStability, 0.1, 0.01)
                .expect("VSM fast");

        assert!(
            margin.lambda_max > 0.0,
            "Lambda max should be positive: {:.4}",
            margin.lambda_max
        );
    }

    #[test]
    fn test_voltage_stability_margin_bisection() {
        let net = load_ieee14();
        let margin =
            compute_voltage_stability_margin(&net, VsaMethod::ContinuationPowerFlow, 0.2, 0.05)
                .expect("VSM bisection");

        assert!(
            margin.lambda_max > 0.0,
            "Lambda max should be positive: {:.4}",
            margin.lambda_max
        );
        assert!(
            margin.lambda_max < 20.0,
            "Lambda max should be reasonable: {:.4}",
            margin.lambda_max
        );
    }
}