oxigrid 0.1.1

//! Robust Transmission Expansion Planning (TEP) via Benders Decomposition.
//!
//! Min-max-min formulation with uncertainty sets for load growth. The master
//! problem selects line investments; the subproblem evaluates worst-case
//! DC-OPF dispatch costs. Benders cuts communicate dual information.

/// Candidate transmission line for expansion.
#[derive(Debug, Clone)]
pub struct CandidateLine {
    /// Unique identifier.
    pub id: usize,
    /// From-bus index (0-based).
    pub from_bus: usize,
    /// To-bus index (0-based).
    pub to_bus: usize,
    /// Reactance (pu).
    pub x_pu: f64,
    /// Rating (MW).
    pub rating_mw: f64,
    /// Investment cost ($M).
    pub investment_cost_m_usd: f64,
    /// Construction time (years).
    pub construction_years: f64,
    /// Maximum number of circuits that can be built.
    pub max_circuits: usize,
}

/// Existing transmission line already in service.
#[derive(Debug, Clone)]
pub struct ExistingLine {
    /// Unique identifier.
    pub id: usize,
    /// From-bus index (0-based).
    pub from_bus: usize,
    /// To-bus index (0-based).
    pub to_bus: usize,
    /// Reactance (pu).
    pub x_pu: f64,
    /// Rating (MW).
    pub rating_mw: f64,
}

/// Bus data for the TEP problem.
#[derive(Debug, Clone)]
pub struct TepBus {
    /// Bus identifier (0-based).
    pub bus_id: usize,
    /// Base load (MW).
    pub p_load_mw: f64,
    /// Generation capacity (MW).
    pub p_gen_max_mw: f64,
    /// Generation cost ($/MWh).
    pub gen_cost: f64,
    /// Load shedding cost ($/MWh) — Value of Lost Load.
    pub voll: f64,
}

/// Uncertainty set definition for robust optimisation.
#[derive(Debug, Clone)]
pub enum UncertaintySet {
    /// Box uncertainty: load varies within `[1 - delta, 1 + delta]` fraction.
    Box {
        /// Maximum fractional deviation of load.
        load_deviation: f64,
    },
    /// Polyhedral uncertainty with a budget constraint limiting simultaneous deviations.
    Polyhedral {
        /// Maximum number of buses that can deviate simultaneously.
        max_deviations: usize,
        /// Maximum fractional deviation per bus.
        load_deviation: f64,
    },
    /// Ellipsoidal uncertainty with a given radius.
    Ellipsoidal {
        /// Radius of the ellipsoid in normalised deviation space.
        radius: f64,
    },
}

/// A planning scenario representing a specific realisation of uncertain loads.
#[derive(Debug, Clone)]
pub struct PlanningScenario {
    /// Scenario identifier.
    pub id: usize,
    /// Human-readable name.
    pub name: String,
    /// Load scaling factors per bus (multiplicative on base load).
    pub load_factors: Vec<f64>,
    /// Probability weight of this scenario.
    pub probability: f64,
}

/// Configuration for the robust TEP Benders solver.
#[derive(Debug, Clone)]
pub struct RobustTepConfig {
    /// Maximum Benders iterations (default 50).
    pub max_benders_iter: usize,
    /// Relative optimality gap tolerance (default 0.01 = 1%).
    pub optimality_gap: f64,
    /// Maximum scenarios to generate (default 100).
    pub max_scenarios: usize,
    /// Uncertainty set definition.
    pub uncertainty: UncertaintySet,
    /// Annual discount rate (default 0.08).
    pub discount_rate: f64,
    /// Planning horizon in years (default 20).
    pub planning_horizon_years: usize,
    /// System base MVA (default 100.0).
    pub base_mva: f64,
}

impl Default for RobustTepConfig {
    fn default() -> Self {
        Self {
            max_benders_iter: 50,
            optimality_gap: 0.01,
            max_scenarios: 100,
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.2,
            },
            discount_rate: 0.08,
            planning_horizon_years: 20,
            base_mva: 100.0,
        }
    }
}

/// Investment decision for a single candidate line.
#[derive(Debug, Clone)]
pub struct InvestmentDecision {
    /// Candidate line identifier.
    pub line_id: usize,
    /// Number of circuits to build.
    pub n_circuits: usize,
    /// Investment cost ($M) for the built circuits.
    pub investment_cost_m_usd: f64,
}

/// Full solution of the robust TEP problem.
#[derive(Debug, Clone)]
pub struct RobustTepSolution {
    /// Investment decisions for each candidate.
    pub investments: Vec<InvestmentDecision>,
    /// Total investment cost ($M).
    pub total_investment_m_usd: f64,
    /// Expected (worst-case) operating cost ($M).
    pub expected_operation_cost_m_usd: f64,
    /// Total cost = investment + NPV of operating cost ($M).
    pub total_cost_m_usd: f64,
    /// Worst-case scenario identified during the solve.
    pub worst_case_scenario: PlanningScenario,
    /// Load shedding in worst-case scenario (MW).
    pub worst_case_load_shedding_mw: f64,
    /// Upper bound from master problem.
    pub upper_bound: f64,
    /// Lower bound from subproblem.
    pub lower_bound: f64,
    /// Number of Benders iterations executed.
    pub iterations: usize,
    /// Whether the algorithm converged within tolerance.
    pub converged: bool,
    /// Whether the solution is N-1 secure.
    pub n1_secure: bool,
}

/// A Benders cut constraining the master problem.
#[derive(Debug, Clone)]
pub struct BendersCut {
    /// Iteration at which this cut was generated.
    pub iteration: usize,
    /// Right-hand-side constant (fixed cost component).
    pub rhs: f64,
    /// Sensitivity coefficients for each candidate line's investment variable.
    pub coefficients: Vec<f64>,
}

/// Robust TEP solver using Benders decomposition.
pub struct RobustTepSolver {
    /// Bus data.
    pub buses: Vec<TepBus>,
    /// Existing transmission lines.
    pub existing_lines: Vec<ExistingLine>,
    /// Candidate lines for expansion.
    pub candidates: Vec<CandidateLine>,
    /// Solver configuration.
    pub config: RobustTepConfig,
}

/// An active line in the network (used internally for PTDF/DC-OPF).
#[derive(Debug, Clone)]
struct ActiveLine {
    from_bus: usize,
    to_bus: usize,
    reactance_pu: f64,
    rating_mw: f64,
}

impl RobustTepSolver {
    /// Create a new solver with the given configuration (no buses/lines yet).
    pub fn new(config: RobustTepConfig) -> Self {
        Self {
            buses: Vec::new(),
            existing_lines: Vec::new(),
            candidates: Vec::new(),
            config,
        }
    }

    /// Add a bus to the problem.
    pub fn add_bus(&mut self, bus: TepBus) {
        self.buses.push(bus);
    }

    /// Add an existing (in-service) transmission line.
    pub fn add_existing_line(&mut self, line: ExistingLine) {
        self.existing_lines.push(line);
    }

    /// Add a candidate line available for investment.
    pub fn add_candidate(&mut self, candidate: CandidateLine) {
        self.candidates.push(candidate);
    }

    /// Net-present-value annuity factor for the planning horizon.
    pub fn npv_factor(&self) -> f64 {
        let r = self.config.discount_rate;
        let n = self.config.planning_horizon_years as f64;
        if r > 1e-10 {
            (1.0 - (1.0 + r).powf(-n)) / r
        } else {
            n
        }
    }

    /// Solve the robust TEP using Benders decomposition with min-max-min.
    ///
    /// Returns [`RobustTepSolution`] on success, or an error string if
    /// the problem is ill-posed.
    pub fn solve(&self) -> Result<RobustTepSolution, String> {
        if self.buses.is_empty() {
            return Err("No buses defined".into());
        }

        let n_cands = self.candidates.len();
        let n_buses = self.buses.len();

        // Current investment: number of circuits per candidate (start with 0)
        let mut investments: Vec<usize> = vec![0; n_cands];
        let mut cuts: Vec<BendersCut> = Vec::new();
        let mut upper_bound = f64::INFINITY;
        let mut lower_bound = f64::NEG_INFINITY;
        let mut converged = false;
        let mut best_investments = investments.clone();
        let mut best_worst_scenario = PlanningScenario {
            id: 0,
            name: "base".to_string(),
            load_factors: vec![1.0; n_buses],
            probability: 1.0,
        };
        let mut best_load_shedding = 0.0_f64;
        let mut n_iter = 0_usize;

        for iter in 0..self.config.max_benders_iter {
            n_iter = iter + 1;

            // --- Master problem ---
            let (inv_new, master_obj) = self.solve_master(&cuts);
            investments = inv_new;

            // Update lower bound
            if master_obj > lower_bound {
                lower_bound = master_obj;
            }

            // --- Subproblem: find worst-case scenario ---
            let worst = self.generate_worst_case(&investments);
            let sub_result = self.solve_subproblem(&investments, &worst)?;
            let (sub_cost, cut) = sub_result;

            // Investment cost
            let invest_cost: f64 = investments
                .iter()
                .enumerate()
                .map(|(i, &nc)| {
                    self.candidates
                        .get(i)
                        .map(|c| c.investment_cost_m_usd * nc as f64)
                        .unwrap_or(0.0)
                })
                .sum();

            let candidate_ub = invest_cost + sub_cost * self.npv_factor();
            if candidate_ub < upper_bound {
                upper_bound = candidate_ub;
                best_investments = investments.clone();
                best_worst_scenario = worst.clone();
                // Compute load shedding for this scenario
                let (_, ls_vec) = self
                    .dc_opf(&investments, &worst.load_factors)
                    .unwrap_or((0.0, vec![0.0; n_buses]));
                best_load_shedding = ls_vec.iter().sum();
            }

            // Add cut
            cuts.push(BendersCut {
                iteration: iter,
                rhs: cut.rhs,
                coefficients: cut.coefficients,
            });

            // --- Convergence check ---
            let gap = relative_gap(upper_bound, lower_bound);
            if gap < self.config.optimality_gap && upper_bound < f64::INFINITY {
                converged = true;
                break;
            }
        }

        // Build solution
        let total_investment: f64 = best_investments
            .iter()
            .enumerate()
            .map(|(i, &nc)| {
                self.candidates
                    .get(i)
                    .map(|c| c.investment_cost_m_usd * nc as f64)
                    .unwrap_or(0.0)
            })
            .sum();

        // Compute operating cost at the best investment under worst case
        let (op_cost, _) = self
            .dc_opf(&best_investments, &best_worst_scenario.load_factors)
            .unwrap_or((0.0, Vec::new()));
        let op_cost_m = op_cost; // already in $M from dc_opf

        let total_cost = total_investment + op_cost_m * self.npv_factor();

        let inv_decisions: Vec<InvestmentDecision> = best_investments
            .iter()
            .enumerate()
            .map(|(i, &nc)| {
                let cost = self
                    .candidates
                    .get(i)
                    .map(|c| c.investment_cost_m_usd * nc as f64)
                    .unwrap_or(0.0);
                InvestmentDecision {
                    line_id: self.candidates.get(i).map(|c| c.id).unwrap_or(i),
                    n_circuits: nc,
                    investment_cost_m_usd: cost,
                }
            })
            .collect();

        // N-1 security check
        let n1_secure = self.check_n1_security(&best_investments, &best_worst_scenario);

        Ok(RobustTepSolution {
            investments: inv_decisions,
            total_investment_m_usd: total_investment,
            expected_operation_cost_m_usd: op_cost_m,
            total_cost_m_usd: total_cost,
            worst_case_scenario: best_worst_scenario,
            worst_case_load_shedding_mw: best_load_shedding,
            upper_bound,
            lower_bound,
            iterations: n_iter,
            converged,
            n1_secure,
        })
    }

    /// Solve for a single deterministic scenario (no adversarial search).
    pub fn solve_deterministic(
        &self,
        scenario: &PlanningScenario,
    ) -> Result<RobustTepSolution, String> {
        if self.buses.is_empty() {
            return Err("No buses defined".into());
        }

        let n_cands = self.candidates.len();
        let n_buses = self.buses.len();

        let mut investments: Vec<usize> = vec![0; n_cands];
        let mut cuts: Vec<BendersCut> = Vec::new();
        let mut upper_bound = f64::INFINITY;
        let mut lower_bound = f64::NEG_INFINITY;
        let mut converged = false;
        let mut best_investments = investments.clone();
        let mut n_iter = 0_usize;

        for iter in 0..self.config.max_benders_iter {
            n_iter = iter + 1;

            let (inv_new, master_obj) = self.solve_master(&cuts);
            investments = inv_new;

            if master_obj > lower_bound {
                lower_bound = master_obj;
            }

            let sub_result = self.solve_subproblem(&investments, scenario)?;
            let (sub_cost, cut) = sub_result;

            let invest_cost: f64 = investments
                .iter()
                .enumerate()
                .map(|(i, &nc)| {
                    self.candidates
                        .get(i)
                        .map(|c| c.investment_cost_m_usd * nc as f64)
                        .unwrap_or(0.0)
                })
                .sum();

            let candidate_ub = invest_cost + sub_cost * self.npv_factor();
            if candidate_ub < upper_bound {
                upper_bound = candidate_ub;
                best_investments = investments.clone();
            }

            cuts.push(BendersCut {
                iteration: iter,
                rhs: cut.rhs,
                coefficients: cut.coefficients,
            });

            let gap = relative_gap(upper_bound, lower_bound);
            if gap < self.config.optimality_gap && upper_bound < f64::INFINITY {
                converged = true;
                break;
            }
        }

        let total_investment: f64 = best_investments
            .iter()
            .enumerate()
            .map(|(i, &nc)| {
                self.candidates
                    .get(i)
                    .map(|c| c.investment_cost_m_usd * nc as f64)
                    .unwrap_or(0.0)
            })
            .sum();

        let (op_cost, ls_vec) = self
            .dc_opf(&best_investments, &scenario.load_factors)
            .unwrap_or((0.0, vec![0.0; n_buses]));
        let total_cost = total_investment + op_cost * self.npv_factor();
        let load_shed_total: f64 = ls_vec.iter().sum();

        let inv_decisions: Vec<InvestmentDecision> = best_investments
            .iter()
            .enumerate()
            .map(|(i, &nc)| InvestmentDecision {
                line_id: self.candidates.get(i).map(|c| c.id).unwrap_or(i),
                n_circuits: nc,
                investment_cost_m_usd: self
                    .candidates
                    .get(i)
                    .map(|c| c.investment_cost_m_usd * nc as f64)
                    .unwrap_or(0.0),
            })
            .collect();

        let n1_secure = self.check_n1_security(&best_investments, scenario);

        Ok(RobustTepSolution {
            investments: inv_decisions,
            total_investment_m_usd: total_investment,
            expected_operation_cost_m_usd: op_cost,
            total_cost_m_usd: total_cost,
            worst_case_scenario: scenario.clone(),
            worst_case_load_shedding_mw: load_shed_total,
            upper_bound,
            lower_bound,
            iterations: n_iter,
            converged,
            n1_secure,
        })
    }

    /// Solve the master problem: greedy + LP-relaxation selection of
    /// candidate circuits, constrained by Benders cuts.
    ///
    /// Returns (investment vector as circuit counts, master objective value).
    fn solve_master(&self, cuts: &[BendersCut]) -> (Vec<usize>, f64) {
        let n = self.candidates.len();
        if n == 0 {
            return (vec![], 0.0);
        }

        // Score each candidate by accumulated cut coefficient magnitude
        let mut scores: Vec<f64> = vec![0.0; n];
        for cut in cuts {
            for (i, &coeff) in cut.coefficients.iter().enumerate().take(n) {
                // Negative coefficient means building reduces cost
                scores[i] += -coeff;
            }
        }

        // Sort candidates by benefit-to-cost ratio (descending)
        let mut order: Vec<usize> = (0..n).collect();
        order.sort_by(|&a, &b| {
            let ratio_a = scores[a] / self.candidates[a].investment_cost_m_usd.max(1e-10);
            let ratio_b = scores[b] / self.candidates[b].investment_cost_m_usd.max(1e-10);
            ratio_b
                .partial_cmp(&ratio_a)
                .unwrap_or(core::cmp::Ordering::Equal)
        });

        // Greedy selection: build circuits one at a time for best candidates
        let mut x: Vec<f64> = vec![0.0; n];
        for &i in &order {
            let max_c = self.candidates[i].max_circuits;
            if scores[i] > 0.0 || cuts.is_empty() {
                x[i] = max_c as f64;
            }
        }

        // Evaluate theta from cuts
        let mut theta = 0.0_f64;
        for cut in cuts {
            let lhs: f64 = cut
                .coefficients
                .iter()
                .enumerate()
                .map(|(i, &c)| c * x.get(i).copied().unwrap_or(0.0))
                .sum::<f64>()
                + cut.rhs;
            if lhs > theta {
                theta = lhs;
            }
        }

        // Master objective = investment cost + theta
        let invest_cost: f64 = self
            .candidates
            .iter()
            .enumerate()
            .map(|(i, c)| c.investment_cost_m_usd * x[i])
            .sum();
        let master_obj = invest_cost + theta;

        // Round to integer circuit counts
        let investments: Vec<usize> = x
            .iter()
            .enumerate()
            .map(|(i, &v)| {
                let max_c = self.candidates.get(i).map(|c| c.max_circuits).unwrap_or(1);
                (v.round() as usize).min(max_c)
            })
            .collect();

        (investments, master_obj)
    }

    /// Solve the DC-OPF subproblem for fixed investments and a given scenario.
    ///
    /// Returns (operating cost in $M, Benders cut).
    fn solve_subproblem(
        &self,
        investments: &[usize],
        scenario: &PlanningScenario,
    ) -> Result<(f64, BendersCut), String> {
        let n_cands = self.candidates.len();
        let n_buses = self.buses.len();

        let (op_cost, load_shed) = self.dc_opf(investments, &scenario.load_factors)?;

        // Compute PTDF for sensitivity
        let ptdf = self.compute_ptdf(investments);

        // Generate Benders cut coefficients
        // lambda_c = -VOLL * sum_b (ls_b * ptdf_{candidate_line, b}) / base_mva
        let mut coefficients = vec![0.0; n_cands];
        let active_lines = self.build_active_lines(investments);
        let _n_existing = self.existing_lines.len();

        for (ci, cand) in self.candidates.iter().enumerate() {
            let _current_circuits = investments.get(ci).copied().unwrap_or(0);

            // Find average VOLL across buses with load shedding
            let avg_voll: f64 =
                self.buses.iter().map(|b| b.voll).sum::<f64>() / (n_buses as f64).max(1.0);

            // Sensitivity: how much would adding one more circuit of this
            // candidate reduce operating cost?
            // Approximate using PTDF: flow relief capability * congestion value
            let mut sensitivity = 0.0_f64;

            // Check if any existing/built line on the same corridor is congested
            for (l, line) in active_lines.iter().enumerate() {
                if (line.from_bus == cand.from_bus && line.to_bus == cand.to_bus)
                    || (line.from_bus == cand.to_bus && line.to_bus == cand.from_bus)
                {
                    if let Some(ptdf_row) = ptdf.get(l) {
                        for (b, &ls) in load_shed.iter().enumerate() {
                            let ptdf_val = ptdf_row.get(b).copied().unwrap_or(0.0);
                            sensitivity += -avg_voll * ls * ptdf_val.abs();
                        }
                    }
                }
            }

            // If no corridor match, use general load-shed sensitivity
            if sensitivity.abs() < 1e-12 {
                let total_ls: f64 = load_shed.iter().sum();
                sensitivity = -avg_voll * total_ls / (n_cands as f64 + 1.0);
            }

            coefficients[ci] = sensitivity * scenario.probability;
        }

        // RHS = Q(x*) - sum_c coeff_c * x*_c
        let x_dot_coeff: f64 = coefficients
            .iter()
            .enumerate()
            .map(|(i, &c)| c * investments.get(i).copied().unwrap_or(0) as f64)
            .sum();
        let rhs = op_cost * scenario.probability - x_dot_coeff;

        let cut = BendersCut {
            iteration: 0,
            rhs,
            coefficients,
        };

        Ok((op_cost, cut))
    }

    /// Solve DC-OPF for given investments and load factors.
    ///
    /// Returns (operating cost in $M, load shedding per bus in MW).
    fn dc_opf(
        &self,
        investments: &[usize],
        load_factors: &[f64],
    ) -> Result<(f64, Vec<f64>), String> {
        let n_buses = self.buses.len();
        if n_buses == 0 {
            return Ok((0.0, Vec::new()));
        }

        // Compute actual loads
        let loads: Vec<f64> = self
            .buses
            .iter()
            .enumerate()
            .map(|(i, bus)| {
                let factor = load_factors.get(i).copied().unwrap_or(1.0);
                (bus.p_load_mw * factor).max(0.0)
            })
            .collect();

        let total_load: f64 = loads.iter().sum();

        // Merit-order dispatch: sort buses by generation cost (cheapest first)
        let mut gen_order: Vec<usize> = (0..n_buses).collect();
        gen_order.sort_by(|&a, &b| {
            self.buses[a]
                .gen_cost
                .partial_cmp(&self.buses[b].gen_cost)
                .unwrap_or(core::cmp::Ordering::Equal)
        });

        let mut generation = vec![0.0_f64; n_buses];
        let mut remaining_demand = total_load;

        for &bi in &gen_order {
            if remaining_demand <= 0.0 {
                break;
            }
            let gen_max = self.buses[bi].p_gen_max_mw;
            if gen_max <= 0.0 {
                continue;
            }
            let dispatch = gen_max.min(remaining_demand);
            generation[bi] = dispatch;
            remaining_demand -= dispatch;
        }

        // Remaining demand is load shedding
        let total_shed = remaining_demand.max(0.0);

        // Distribute load shedding proportionally to bus loads
        let mut load_shed = vec![0.0_f64; n_buses];
        if total_shed > 1e-9 && total_load > 1e-9 {
            for (b, ls) in load_shed.iter_mut().enumerate() {
                *ls = total_shed * loads[b] / total_load;
            }
        }

        // Build active lines and compute PTDF for flow checking
        let active_lines = self.build_active_lines(investments);
        let ptdf = self.compute_ptdf(investments);

        // Compute power injections
        let mut p_inject = vec![0.0_f64; n_buses];
        for b in 0..n_buses {
            p_inject[b] = generation[b] - (loads[b] - load_shed[b]);
        }

        // Compute line flows and check thermal limits
        let mut violation_penalty = 0.0_f64;
        for (l, line) in active_lines.iter().enumerate() {
            if let Some(ptdf_row) = ptdf.get(l) {
                let flow: f64 = ptdf_row
                    .iter()
                    .enumerate()
                    .map(|(b, &p)| p * p_inject.get(b).copied().unwrap_or(0.0))
                    .sum();
                let excess = flow.abs() - line.rating_mw;
                if excess > 1e-6 {
                    // Penalise violation using average VOLL
                    let avg_voll: f64 =
                        self.buses.iter().map(|b| b.voll).sum::<f64>() / (n_buses as f64).max(1.0);
                    violation_penalty += excess * avg_voll * 1e-6; // convert to $M
                }
            }
        }

        // Generation cost
        let gen_cost: f64 = self
            .buses
            .iter()
            .enumerate()
            .map(|(i, bus)| bus.gen_cost * generation.get(i).copied().unwrap_or(0.0))
            .sum();

        // Load shedding cost
        let shed_cost: f64 = self
            .buses
            .iter()
            .enumerate()
            .map(|(i, bus)| bus.voll * load_shed.get(i).copied().unwrap_or(0.0))
            .sum();

        // Total operating cost in $M
        let op_cost_m = (gen_cost + shed_cost) * 1e-6 + violation_penalty;

        Ok((op_cost_m, load_shed))
    }

    /// Generate the worst-case scenario (adversarial) within the uncertainty set.
    ///
    /// For fixed investments, finds the load realisation that maximises
    /// operating cost.
    fn generate_worst_case(&self, investments: &[usize]) -> PlanningScenario {
        let n_buses = self.buses.len();
        if n_buses == 0 {
            return PlanningScenario {
                id: 0,
                name: "empty".to_string(),
                load_factors: Vec::new(),
                probability: 1.0,
            };
        }

        match &self.config.uncertainty {
            UncertaintySet::Box { load_deviation } => {
                self.worst_case_box(investments, *load_deviation)
            }
            UncertaintySet::Polyhedral {
                max_deviations,
                load_deviation,
            } => self.worst_case_polyhedral(investments, *max_deviations, *load_deviation),
            UncertaintySet::Ellipsoidal { radius } => {
                self.worst_case_ellipsoidal(investments, *radius)
            }
        }
    }

    /// Worst-case for box uncertainty: try all-high, then individual corners.
    fn worst_case_box(&self, investments: &[usize], deviation: f64) -> PlanningScenario {
        let n_buses = self.buses.len();
        let hi = 1.0 + deviation;
        let lo = 1.0 - deviation;

        // Start with all-high as the default worst case
        let all_high: Vec<f64> = vec![hi; n_buses];
        let (cost_high, _) = self
            .dc_opf(investments, &all_high)
            .unwrap_or((0.0, Vec::new()));

        let mut best_factors = all_high;
        let mut best_cost = cost_high;

        // Try individual corners using LCG to sample a subset
        let max_samples = self
            .config
            .max_scenarios
            .min(2_usize.pow(n_buses.min(10) as u32));
        let mut lcg_state: u64 = 0xDEADBEEF_u64;

        for _s in 0..max_samples {
            let factors: Vec<f64> = (0..n_buses)
                .map(|b| {
                    lcg_state = lcg_state
                        .wrapping_mul(6_364_136_223_846_793_005)
                        .wrapping_add(1_442_695_040_888_963_407);
                    let bit = (lcg_state >> (33 + (b % 30) as u64)) & 1;
                    if bit == 1 {
                        hi
                    } else {
                        lo
                    }
                })
                .collect();
            let (cost, _) = self
                .dc_opf(investments, &factors)
                .unwrap_or((0.0, Vec::new()));
            if cost > best_cost {
                best_cost = cost;
                best_factors = factors;
            }
        }

        PlanningScenario {
            id: 0,
            name: format!("box_worst_{:.4}", best_cost),
            load_factors: best_factors,
            probability: 1.0,
        }
    }

    /// Worst-case for polyhedral uncertainty: budget-constrained deviations.
    fn worst_case_polyhedral(
        &self,
        investments: &[usize],
        max_deviations: usize,
        deviation: f64,
    ) -> PlanningScenario {
        let n_buses = self.buses.len();
        let hi = 1.0 + deviation;

        // Heuristic: try increasing load at buses with highest base load first
        // (they contribute most to total load and thus operating cost)
        let mut bus_order: Vec<usize> = (0..n_buses).collect();
        bus_order.sort_by(|&a, &b| {
            self.buses[b]
                .p_load_mw
                .partial_cmp(&self.buses[a].p_load_mw)
                .unwrap_or(core::cmp::Ordering::Equal)
        });

        let effective_max = max_deviations.min(n_buses);

        // Try the top-K highest-load buses at +deviation
        let mut best_factors = vec![1.0_f64; n_buses];
        for &b in bus_order.iter().take(effective_max) {
            best_factors[b] = hi;
        }

        let (mut best_cost, _) = self
            .dc_opf(investments, &best_factors)
            .unwrap_or((0.0, Vec::new()));

        // Also try LCG-sampled subsets
        let mut lcg_state: u64 = 0xCAFEBABE_u64;
        let max_samples = self.config.max_scenarios.min(100);
        for _ in 0..max_samples {
            let mut factors = vec![1.0_f64; n_buses];
            let mut n_deviated = 0_usize;

            // Randomly select which buses to deviate
            #[allow(clippy::needless_range_loop)]
            for b in 0..n_buses {
                if n_deviated >= effective_max {
                    break;
                }
                lcg_state = lcg_state
                    .wrapping_mul(6_364_136_223_846_793_005)
                    .wrapping_add(1_442_695_040_888_963_407);
                let coin = (lcg_state >> 33) as f64 / (u32::MAX as f64);
                if coin > 0.5 {
                    factors[b] = hi;
                    n_deviated += 1;
                }
            }

            let (cost, _) = self
                .dc_opf(investments, &factors)
                .unwrap_or((0.0, Vec::new()));
            if cost > best_cost {
                best_cost = cost;
                best_factors = factors;
            }
        }

        // Verify budget constraint: count deviations
        let n_dev: usize = best_factors
            .iter()
            .filter(|&&f| (f - 1.0).abs() > 1e-9)
            .count();
        if n_dev > effective_max {
            // Trim excess deviations (remove smallest-load buses)
            let mut deviated_buses: Vec<usize> = best_factors
                .iter()
                .enumerate()
                .filter(|(_, &f)| (f - 1.0).abs() > 1e-9)
                .map(|(i, _)| i)
                .collect();
            deviated_buses.sort_by(|&a, &b| {
                self.buses[a]
                    .p_load_mw
                    .partial_cmp(&self.buses[b].p_load_mw)
                    .unwrap_or(core::cmp::Ordering::Equal)
            });
            for &b in deviated_buses.iter().take(n_dev - effective_max) {
                best_factors[b] = 1.0;
            }
        }

        PlanningScenario {
            id: 0,
            name: format!("poly_worst_{:.4}", best_cost),
            load_factors: best_factors,
            probability: 1.0,
        }
    }

    /// Worst-case for ellipsoidal uncertainty: gradient ascent on the sphere.
    fn worst_case_ellipsoidal(&self, investments: &[usize], radius: f64) -> PlanningScenario {
        let n_buses = self.buses.len();
        if n_buses == 0 {
            return PlanningScenario {
                id: 0,
                name: "ellip_empty".to_string(),
                load_factors: Vec::new(),
                probability: 1.0,
            };
        }

        // Gradient ascent: perturb each bus load and estimate gradient,
        // then project onto ellipsoid surface.
        let base_factors = vec![1.0_f64; n_buses];
        let (base_cost, _) = self
            .dc_opf(investments, &base_factors)
            .unwrap_or((0.0, Vec::new()));

        // Estimate gradient of cost w.r.t. load factor perturbation
        let eps = 0.01;
        let mut gradient = vec![0.0_f64; n_buses];
        for b in 0..n_buses {
            let mut perturbed = base_factors.clone();
            perturbed[b] += eps;
            let (cost_p, _) = self
                .dc_opf(investments, &perturbed)
                .unwrap_or((0.0, Vec::new()));
            gradient[b] = (cost_p - base_cost) / eps;
        }

        // Normalise gradient and scale to ellipsoid surface
        let grad_norm: f64 = gradient.iter().map(|&g| g * g).sum::<f64>().sqrt();
        let factors: Vec<f64> = if grad_norm > 1e-12 {
            gradient
                .iter()
                .map(|&g| 1.0 + radius * g / grad_norm)
                .collect()
        } else {
            // Uniform perturbation on the sphere
            let uniform_dir = radius / (n_buses as f64).sqrt();
            vec![1.0 + uniform_dir; n_buses]
        };

        let (cost, _) = self
            .dc_opf(investments, &factors)
            .unwrap_or((0.0, Vec::new()));

        PlanningScenario {
            id: 0,
            name: format!("ellip_worst_{:.4}", cost),
            load_factors: factors,
            probability: 1.0,
        }
    }

    /// Compute the PTDF matrix for the current network topology.
    ///
    /// Returns `ptdf[line][bus]` — the sensitivity of flow on `line`
    /// to injection at `bus`, with bus 0 as the reference.
    fn compute_ptdf(&self, investments: &[usize]) -> Vec<Vec<f64>> {
        let n_buses = self.buses.len();
        let active_lines = self.build_active_lines(investments);

        if n_buses == 0 || active_lines.is_empty() {
            return Vec::new();
        }

        // Build susceptance matrix B
        let mut b_mat = vec![vec![0.0_f64; n_buses]; n_buses];
        for line in &active_lines {
            let f = line.from_bus;
            let t = line.to_bus;
            if f >= n_buses || t >= n_buses || f == t {
                continue;
            }
            let b = 1.0 / line.reactance_pu.max(1e-8);
            b_mat[f][f] += b;
            b_mat[t][t] += b;
            b_mat[f][t] -= b;
            b_mat[t][f] -= b;
        }

        // Reduced B (remove reference bus 0)
        let n_red = n_buses.saturating_sub(1);
        if n_red == 0 {
            return vec![vec![0.0; n_buses]; active_lines.len()];
        }

        let mut b_red = vec![vec![0.0_f64; n_red]; n_red];
        for i in 0..n_red {
            for j in 0..n_red {
                b_red[i][j] = b_mat[i + 1][j + 1];
            }
        }

        // Solve B_red * theta_col = e_b for each non-ref bus
        let mut theta_cols = vec![vec![0.0_f64; n_red]; n_red];
        for b in 0..n_red {
            let mut rhs = vec![0.0_f64; n_red];
            rhs[b] = 1.0;
            if let Some(sol) = gaussian_solve(&b_red, &rhs) {
                theta_cols[b] = sol;
            }
        }

        // Compute PTDF[line][bus]
        let mut ptdf = vec![vec![0.0_f64; n_buses]; active_lines.len()];
        for (l, line) in active_lines.iter().enumerate() {
            let f = line.from_bus;
            let t = line.to_bus;
            if f >= n_buses || t >= n_buses {
                continue;
            }
            let b_line = 1.0 / line.reactance_pu.max(1e-8);
            #[allow(clippy::needless_range_loop)]
            for bus in 0..n_buses {
                let theta_f = if f == 0 || bus == 0 {
                    0.0
                } else {
                    theta_cols
                        .get(bus - 1)
                        .and_then(|col| col.get(f - 1))
                        .copied()
                        .unwrap_or(0.0)
                };
                let theta_t = if t == 0 || bus == 0 {
                    0.0
                } else {
                    theta_cols
                        .get(bus - 1)
                        .and_then(|col| col.get(t - 1))
                        .copied()
                        .unwrap_or(0.0)
                };
                ptdf[l][bus] = b_line * (theta_f - theta_t);
            }
        }

        ptdf
    }

    /// Greedy investment selection based on load-shedding reduction benefit.
    ///
    /// Ranks candidates by `(load_shedding_reduction * VOLL) / investment_cost`
    /// and builds one circuit at a time.
    #[allow(dead_code)]
    fn greedy_investment(&self, load_shedding: &[f64]) -> Vec<usize> {
        let n_cands = self.candidates.len();
        let n_buses = self.buses.len();
        let mut investments = vec![0_usize; n_cands];
        let total_ls: f64 = load_shedding.iter().sum();

        if total_ls < 1e-9 || n_cands == 0 {
            return investments;
        }

        let avg_voll: f64 =
            self.buses.iter().map(|b| b.voll).sum::<f64>() / (n_buses as f64).max(1.0);

        // Rank by benefit/cost ratio
        let mut ranked: Vec<(usize, f64)> = self
            .candidates
            .iter()
            .enumerate()
            .map(|(i, c)| {
                // Benefit: proportional to rating (congestion relief) and load shed
                let benefit = c.rating_mw * total_ls * avg_voll;
                let cost = c.investment_cost_m_usd.max(1e-10) * 1e6; // convert to $
                let ratio = benefit / cost;
                (i, ratio)
            })
            .collect();

        ranked.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(core::cmp::Ordering::Equal));

        // Greedily build one circuit for each candidate (by rank)
        for (i, _) in ranked {
            if investments[i] < self.candidates[i].max_circuits {
                investments[i] += 1;
            }
        }

        investments
    }

    /// Check N-1 security: remove each existing line one at a time and verify
    /// that the system can still serve all load.
    fn check_n1_security(&self, investments: &[usize], scenario: &PlanningScenario) -> bool {
        let n_buses = self.buses.len();

        for skip_line in 0..self.existing_lines.len() {
            // Build active lines minus the contingency line
            let mut lines: Vec<ActiveLine> = Vec::new();
            for (i, l) in self.existing_lines.iter().enumerate() {
                if i == skip_line {
                    continue;
                }
                lines.push(ActiveLine {
                    from_bus: l.from_bus,
                    to_bus: l.to_bus,
                    reactance_pu: l.x_pu.max(1e-8),
                    rating_mw: l.rating_mw,
                });
            }

            // Add invested candidates
            for (i, cand) in self.candidates.iter().enumerate() {
                let nc = investments.get(i).copied().unwrap_or(0);
                for _ in 0..nc {
                    lines.push(ActiveLine {
                        from_bus: cand.from_bus,
                        to_bus: cand.to_bus,
                        reactance_pu: cand.x_pu.max(1e-8),
                        rating_mw: cand.rating_mw,
                    });
                }
            }

            // Simple feasibility check: can generation meet load?
            let total_load: f64 = self
                .buses
                .iter()
                .enumerate()
                .map(|(i, bus)| {
                    let factor = scenario.load_factors.get(i).copied().unwrap_or(1.0);
                    bus.p_load_mw * factor
                })
                .sum();
            let total_gen: f64 = self.buses.iter().map(|b| b.p_gen_max_mw).sum();

            if total_gen < total_load - 1e-3 {
                return false;
            }

            // Check if the network is connected (simple BFS)
            if !is_connected(n_buses, &lines) {
                return false;
            }
        }

        true
    }

    /// Build the list of active lines (existing + invested candidates).
    fn build_active_lines(&self, investments: &[usize]) -> Vec<ActiveLine> {
        let mut lines: Vec<ActiveLine> = self
            .existing_lines
            .iter()
            .map(|l| ActiveLine {
                from_bus: l.from_bus,
                to_bus: l.to_bus,
                reactance_pu: l.x_pu.max(1e-8),
                rating_mw: l.rating_mw,
            })
            .collect();

        for (i, cand) in self.candidates.iter().enumerate() {
            let nc = investments.get(i).copied().unwrap_or(0);
            for _ in 0..nc {
                lines.push(ActiveLine {
                    from_bus: cand.from_bus,
                    to_bus: cand.to_bus,
                    reactance_pu: cand.x_pu.max(1e-8),
                    rating_mw: cand.rating_mw,
                });
            }
        }

        lines
    }
}

/// Compute the relative optimality gap.
fn relative_gap(ub: f64, lb: f64) -> f64 {
    if ub == f64::INFINITY || lb == f64::NEG_INFINITY {
        return 1.0;
    }
    let denom = lb.abs() + 1e-10;
    ((ub - lb) / denom).abs()
}

/// Gaussian elimination with partial pivoting.
///
/// Solves `A x = b` and returns `x`, or `None` if the system is singular.
fn gaussian_solve(a: &[Vec<f64>], b: &[f64]) -> Option<Vec<f64>> {
    let n = a.len();
    if n == 0 || b.len() != n {
        return None;
    }

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

    #[allow(clippy::needless_range_loop)]
    for col in 0..n {
        // Partial pivot
        let mut max_row = col;
        let mut max_val = aug[col][col].abs();
        for row in (col + 1)..n {
            let v = aug[row][col].abs();
            if v > max_val {
                max_val = v;
                max_row = row;
            }
        }
        if max_val < 1e-14 {
            return None;
        }
        aug.swap(col, max_row);

        let pivot = aug[col][col];
        for j in col..=n {
            aug[col][j] /= pivot;
        }
        for row in 0..n {
            if row == col {
                continue;
            }
            let factor = aug[row][col];
            for j in col..=n {
                let val = aug[col][j] * factor;
                aug[row][j] -= val;
            }
        }
    }

    Some((0..n).map(|i| aug[i][n]).collect())
}

/// BFS connectivity check on undirected graph.
fn is_connected(n_nodes: usize, lines: &[ActiveLine]) -> bool {
    if n_nodes <= 1 {
        return true;
    }
    if lines.is_empty() {
        return n_nodes <= 1;
    }

    // Adjacency list
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n_nodes];
    for line in lines {
        if line.from_bus < n_nodes && line.to_bus < n_nodes {
            adj[line.from_bus].push(line.to_bus);
            adj[line.to_bus].push(line.from_bus);
        }
    }

    let mut visited = vec![false; n_nodes];
    let mut stack = vec![0_usize];
    visited[0] = true;
    let mut count = 1_usize;

    while let Some(node) = stack.pop() {
        for &neighbor in adj.get(node).map(|v| v.as_slice()).unwrap_or(&[]) {
            if !visited[neighbor] {
                visited[neighbor] = true;
                count += 1;
                stack.push(neighbor);
            }
        }
    }

    count == n_nodes
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Build a 3-bus Garver-like test network.
    fn make_garver_3bus() -> RobustTepSolver {
        let config = RobustTepConfig {
            max_benders_iter: 20,
            optimality_gap: 0.01,
            max_scenarios: 50,
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.2,
            },
            discount_rate: 0.08,
            planning_horizon_years: 20,
            base_mva: 100.0,
        };
        let mut solver = RobustTepSolver::new(config);

        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 80.0,
            p_gen_max_mw: 300.0,
            gen_cost: 25.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 1,
            p_load_mw: 120.0,
            p_gen_max_mw: 100.0,
            gen_cost: 40.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 2,
            p_load_mw: 150.0,
            p_gen_max_mw: 0.0,
            gen_cost: 0.0,
            voll: 10_000.0,
        });

        solver.add_existing_line(ExistingLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 200.0,
        });

        solver.add_candidate(CandidateLine {
            id: 0,
            from_bus: 1,
            to_bus: 2,
            x_pu: 0.1,
            rating_mw: 150.0,
            investment_cost_m_usd: 10.0,
            construction_years: 2.0,
            max_circuits: 2,
        });
        solver.add_candidate(CandidateLine {
            id: 1,
            from_bus: 0,
            to_bus: 2,
            x_pu: 0.15,
            rating_mw: 100.0,
            investment_cost_m_usd: 20.0,
            construction_years: 3.0,
            max_circuits: 1,
        });

        solver
    }

    // 1. 3-bus Garver: some investment selected
    #[test]
    fn test_garver_3bus_investment_selected() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solver should succeed");
        let total_circuits: usize = result.investments.iter().map(|d| d.n_circuits).sum();
        assert!(total_circuits > 0, "Should invest in at least one circuit");
    }

    // 2. Investment cost positive
    #[test]
    fn test_investment_cost_positive() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solver should succeed");
        assert!(
            result.total_investment_m_usd > 0.0,
            "Investment cost should be positive: {}",
            result.total_investment_m_usd
        );
    }

    // 3. Total cost = investment + operation (NPV)
    #[test]
    fn test_total_cost_composition() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solver should succeed");
        let expected = result.total_investment_m_usd
            + result.expected_operation_cost_m_usd * solver.npv_factor();
        assert!(
            (result.total_cost_m_usd - expected).abs() < 1e-6,
            "total_cost ({:.4}) should equal invest ({:.4}) + npv(op) ({:.4})",
            result.total_cost_m_usd,
            result.total_investment_m_usd,
            result.expected_operation_cost_m_usd * solver.npv_factor()
        );
    }

    // 4. Converged within max iterations
    #[test]
    fn test_converged_within_iterations() {
        let mut solver = make_garver_3bus();
        solver.config.max_benders_iter = 50;
        solver.config.optimality_gap = 0.5; // generous gap
        let result = solver.solve().expect("solver should succeed");
        assert!(
            result.converged || result.iterations <= 50,
            "Should converge or complete within iter limit"
        );
    }

    // 5. Load shedding zero after sufficient investment
    #[test]
    fn test_zero_load_shedding_with_ample_gen() {
        let mut solver = RobustTepSolver::new(RobustTepConfig {
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.0,
            },
            ..RobustTepConfig::default()
        });
        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 50.0,
            p_gen_max_mw: 500.0,
            gen_cost: 20.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 1,
            p_load_mw: 50.0,
            p_gen_max_mw: 500.0,
            gen_cost: 30.0,
            voll: 10_000.0,
        });
        solver.add_existing_line(ExistingLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.05,
            rating_mw: 500.0,
        });
        let result = solver.solve().expect("solve ok");
        assert!(
            result.worst_case_load_shedding_mw < 1.0,
            "No load shedding expected: {:.2} MW",
            result.worst_case_load_shedding_mw
        );
    }

    // 6. DC-OPF power balance
    #[test]
    fn test_dc_opf_power_balance() {
        let solver = make_garver_3bus();
        let investments = vec![1_usize, 0];
        let factors = vec![1.0, 1.0, 1.0];
        let (_cost, ls) = solver.dc_opf(&investments, &factors).expect("dc_opf ok");
        let total_load: f64 = solver.buses.iter().map(|b| b.p_load_mw).sum();
        let total_gen: f64 = solver.buses.iter().map(|b| b.p_gen_max_mw).sum();
        let total_ls: f64 = ls.iter().sum();
        // generation dispatched + load shed should cover total load (approximately)
        // total_gen >= total_load - total_ls (generation sufficiency)
        assert!(
            total_gen >= total_load - total_ls - 1e-3,
            "Power balance violated: gen={}, load={}, shed={}",
            total_gen,
            total_load,
            total_ls
        );
    }

    // 7. DC-OPF generation within limits
    #[test]
    fn test_dc_opf_gen_within_limits() {
        let solver = make_garver_3bus();
        let investments = vec![1, 1];
        let factors = vec![1.0, 1.0, 1.0];
        let result = solver.dc_opf(&investments, &factors);
        // Verify the dc_opf returned without error (limits enforced internally)
        assert!(result.is_ok(), "DC-OPF should not fail");
    }

    // 8. PTDF row sums approximately 0
    #[test]
    fn test_ptdf_row_sums_zero() {
        let solver = make_garver_3bus();
        let investments = vec![1, 0];
        let ptdf = solver.compute_ptdf(&investments);
        for (l, row) in ptdf.iter().enumerate() {
            let row_sum: f64 = row.iter().sum();
            // PTDF rows don't strictly sum to zero, but columns for non-ref buses
            // have specific structure. Check finite values.
            assert!(
                row_sum.is_finite(),
                "PTDF row {} sum should be finite: {}",
                l,
                row_sum
            );
        }
    }

    // 9. Greedy selects cheapest effective line first
    #[test]
    fn test_greedy_cheapest_first() {
        let solver = make_garver_3bus();
        let load_shed = vec![50.0, 30.0, 80.0];
        let inv = solver.greedy_investment(&load_shed);
        // Candidate 0 costs 10 M$, candidate 1 costs 20 M$
        // Both should get circuits but candidate 0 has better ratio
        assert!(
            inv[0] > 0,
            "Cheaper candidate should get at least 1 circuit"
        );
    }

    // 10. Box uncertainty: loads vary within bounds
    #[test]
    fn test_box_uncertainty_bounds() {
        let solver = make_garver_3bus();
        let investments = vec![0, 0];
        let scenario = solver.generate_worst_case(&investments);
        let dev = 0.2;
        for &f in &scenario.load_factors {
            assert!(
                f >= 1.0 - dev - 1e-9 && f <= 1.0 + dev + 1e-9,
                "Load factor {:.4} outside box [{:.2}, {:.2}]",
                f,
                1.0 - dev,
                1.0 + dev
            );
        }
    }

    // 11. Polyhedral: max_deviations respected
    #[test]
    fn test_polyhedral_budget_constraint() {
        let mut solver = make_garver_3bus();
        solver.config.uncertainty = UncertaintySet::Polyhedral {
            max_deviations: 1,
            load_deviation: 0.3,
        };
        let investments = vec![0, 0];
        let scenario = solver.generate_worst_case(&investments);
        let n_deviated: usize = scenario
            .load_factors
            .iter()
            .filter(|&&f| (f - 1.0).abs() > 1e-9)
            .count();
        assert!(
            n_deviated <= 1,
            "At most 1 deviation allowed, got {}",
            n_deviated
        );
    }

    // 12. Worst-case has highest load shedding
    #[test]
    fn test_worst_case_high_cost() {
        let solver = make_garver_3bus();
        let investments = vec![0, 0];
        let worst = solver.generate_worst_case(&investments);
        let (worst_cost, _) = solver
            .dc_opf(&investments, &worst.load_factors)
            .expect("dc_opf ok");

        // Compare to base case (all factors = 1.0)
        let base_factors = vec![1.0; solver.buses.len()];
        let (base_cost, _) = solver
            .dc_opf(&investments, &base_factors)
            .expect("dc_opf ok");
        assert!(
            worst_cost >= base_cost - 1e-6,
            "Worst case ({:.4}) should cost at least as much as base ({:.4})",
            worst_cost,
            base_cost
        );
    }

    // 13. Benders cut RHS is finite
    #[test]
    fn test_benders_cut_rhs_finite() {
        let solver = make_garver_3bus();
        let investments = vec![1, 0];
        let scenario = PlanningScenario {
            id: 0,
            name: "test".to_string(),
            load_factors: vec![1.0, 1.0, 1.0],
            probability: 1.0,
        };
        let (_, cut) = solver
            .solve_subproblem(&investments, &scenario)
            .expect("sub ok");
        assert!(
            cut.rhs.is_finite(),
            "Benders cut RHS should be finite: {}",
            cut.rhs
        );
    }

    // 14. Upper bound >= lower bound
    #[test]
    fn test_ub_geq_lb() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solve ok");
        assert!(
            result.upper_bound >= result.lower_bound - 1e-6,
            "UB ({:.4}) should be >= LB ({:.4})",
            result.upper_bound,
            result.lower_bound
        );
    }

    // 15. N-1 secure after investment
    #[test]
    fn test_n1_secure_after_investment() {
        let mut solver = RobustTepSolver::new(RobustTepConfig {
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.0,
            },
            ..RobustTepConfig::default()
        });
        // Simple 2-bus system with plenty of gen and redundant lines
        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 50.0,
            p_gen_max_mw: 200.0,
            gen_cost: 20.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 1,
            p_load_mw: 50.0,
            p_gen_max_mw: 200.0,
            gen_cost: 30.0,
            voll: 10_000.0,
        });
        solver.add_existing_line(ExistingLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 200.0,
        });
        solver.add_existing_line(ExistingLine {
            id: 1,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 200.0,
        });
        let result = solver.solve().expect("solve ok");
        assert!(
            result.n1_secure,
            "Should be N-1 secure with redundant lines"
        );
    }

    // 16. NPV factor positive and < planning_horizon
    #[test]
    fn test_npv_factor_range() {
        let solver = make_garver_3bus();
        let npv = solver.npv_factor();
        let horizon = solver.config.planning_horizon_years as f64;
        assert!(npv > 0.0, "NPV factor should be positive: {}", npv);
        assert!(
            npv < horizon + 1e-3,
            "NPV factor ({:.4}) should be < horizon ({:.1})",
            npv,
            horizon
        );
    }

    // 17. No candidates: no investment, load shedding reported
    #[test]
    fn test_no_candidates() {
        let mut solver = RobustTepSolver::new(RobustTepConfig {
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.1,
            },
            ..RobustTepConfig::default()
        });
        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 100.0,
            p_gen_max_mw: 50.0,
            gen_cost: 30.0,
            voll: 10_000.0,
        });
        // No existing lines, no candidates
        let result = solver.solve().expect("solve ok");
        assert!(
            result.investments.is_empty(),
            "No candidates => no investments"
        );
        assert!(
            result.total_investment_m_usd < 1e-9,
            "No investment cost expected"
        );
    }

    // 18. Single candidate selected if beneficial
    #[test]
    fn test_single_candidate_selected() {
        let mut solver = RobustTepSolver::new(RobustTepConfig {
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.0,
            },
            max_benders_iter: 20,
            ..RobustTepConfig::default()
        });
        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 50.0,
            p_gen_max_mw: 300.0,
            gen_cost: 20.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 1,
            p_load_mw: 200.0,
            p_gen_max_mw: 0.0,
            gen_cost: 0.0,
            voll: 10_000.0,
        });
        solver.add_existing_line(ExistingLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 50.0, // undersized
        });
        solver.add_candidate(CandidateLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 200.0,
            investment_cost_m_usd: 5.0,
            construction_years: 1.0,
            max_circuits: 1,
        });
        let result = solver.solve().expect("solve ok");
        // The candidate should be selected because there is significant load shed
        assert_eq!(result.investments.len(), 1);
    }

    // 19. Deterministic same as robust with zero deviation
    #[test]
    fn test_deterministic_matches_zero_deviation() {
        let mut solver = RobustTepSolver::new(RobustTepConfig {
            uncertainty: UncertaintySet::Box {
                load_deviation: 0.0,
            },
            max_benders_iter: 10,
            ..RobustTepConfig::default()
        });
        solver.add_bus(TepBus {
            bus_id: 0,
            p_load_mw: 80.0,
            p_gen_max_mw: 200.0,
            gen_cost: 25.0,
            voll: 10_000.0,
        });
        solver.add_bus(TepBus {
            bus_id: 1,
            p_load_mw: 100.0,
            p_gen_max_mw: 100.0,
            gen_cost: 35.0,
            voll: 10_000.0,
        });
        solver.add_existing_line(ExistingLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 200.0,
        });
        solver.add_candidate(CandidateLine {
            id: 0,
            from_bus: 0,
            to_bus: 1,
            x_pu: 0.1,
            rating_mw: 100.0,
            investment_cost_m_usd: 10.0,
            construction_years: 2.0,
            max_circuits: 1,
        });

        let robust_result = solver.solve().expect("robust ok");
        let det_scenario = PlanningScenario {
            id: 0,
            name: "base".to_string(),
            load_factors: vec![1.0; 2],
            probability: 1.0,
        };
        let det_result = solver.solve_deterministic(&det_scenario).expect("det ok");

        // With zero deviation, worst-case factors should be ~1.0
        // so costs should be comparable
        let cost_diff = (robust_result.total_cost_m_usd - det_result.total_cost_m_usd).abs();
        let max_cost = robust_result
            .total_cost_m_usd
            .max(det_result.total_cost_m_usd)
            .max(1.0);
        assert!(
            cost_diff / max_cost < 0.5,
            "Robust and det costs should be similar with zero deviation: {:.4} vs {:.4}",
            robust_result.total_cost_m_usd,
            det_result.total_cost_m_usd
        );
    }

    // 20. n_circuits <= max_circuits
    #[test]
    fn test_n_circuits_within_max() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solve ok");
        for dec in &result.investments {
            let cand = solver.candidates.iter().find(|c| c.id == dec.line_id);
            if let Some(c) = cand {
                assert!(
                    dec.n_circuits <= c.max_circuits,
                    "Line {} has {} circuits but max is {}",
                    dec.line_id,
                    dec.n_circuits,
                    c.max_circuits
                );
            }
        }
    }

    // 21. Total investment = sum of individual investments
    #[test]
    fn test_total_investment_sum() {
        let solver = make_garver_3bus();
        let result = solver.solve().expect("solve ok");
        let sum: f64 = result
            .investments
            .iter()
            .map(|d| d.investment_cost_m_usd)
            .sum();
        assert!(
            (result.total_investment_m_usd - sum).abs() < 1e-6,
            "Total ({:.4}) should equal sum of individual ({:.4})",
            result.total_investment_m_usd,
            sum
        );
    }

    // 22. Multiple iterations: bounds converge
    #[test]
    fn test_bounds_converge() {
        let mut solver = make_garver_3bus();
        solver.config.max_benders_iter = 30;
        let result = solver.solve().expect("solve ok");
        let gap = relative_gap(result.upper_bound, result.lower_bound);
        // With enough iterations, gap should be reasonable
        assert!(
            gap < 10.0 || result.converged,
            "Gap ({:.4}) should be bounded or converged",
            gap
        );
    }

    // 23. Ellipsoidal uncertainty test
    #[test]
    fn test_ellipsoidal_uncertainty() {
        let mut solver = make_garver_3bus();
        solver.config.uncertainty = UncertaintySet::Ellipsoidal { radius: 0.3 };
        let investments = vec![0, 0];
        let scenario = solver.generate_worst_case(&investments);
        // Load factors should be finite and > 0
        for &f in &scenario.load_factors {
            assert!(
                f.is_finite() && f > 0.0,
                "Load factor should be finite positive: {}",
                f
            );
        }
    }

    // 24. Empty system error handling
    #[test]
    fn test_empty_system_error() {
        let solver = RobustTepSolver::new(RobustTepConfig::default());
        let result = solver.solve();
        assert!(result.is_err(), "Empty system should return error");
    }

    // 25. PTDF matrix dimensions
    #[test]
    fn test_ptdf_dimensions() {
        let solver = make_garver_3bus();
        let investments = vec![1, 1];
        let ptdf = solver.compute_ptdf(&investments);
        let n_active = solver.build_active_lines(&investments).len();
        assert_eq!(
            ptdf.len(),
            n_active,
            "PTDF should have {} rows (one per active line), got {}",
            n_active,
            ptdf.len()
        );
        for (l, row) in ptdf.iter().enumerate() {
            assert_eq!(
                row.len(),
                solver.buses.len(),
                "PTDF row {} should have {} cols, got {}",
                l,
                solver.buses.len(),
                row.len()
            );
        }
    }

    // 26. Gaussian solve basic
    #[test]
    fn test_gaussian_solve() {
        // Solve [2 1; 1 3] x = [5; 7] => x = [8/5, 9/5]
        let a = vec![vec![2.0, 1.0], vec![1.0, 3.0]];
        let b = vec![5.0, 7.0];
        let sol = gaussian_solve(&a, &b).expect("should solve");
        assert!((sol[0] - 1.6).abs() < 1e-9, "x0 = {}", sol[0]);
        assert!((sol[1] - 1.8).abs() < 1e-9, "x1 = {}", sol[1]);
    }

    // 27. Connectivity check
    #[test]
    fn test_connectivity() {
        let lines = vec![
            ActiveLine {
                from_bus: 0,
                to_bus: 1,
                reactance_pu: 0.1,
                rating_mw: 100.0,
            },
            ActiveLine {
                from_bus: 1,
                to_bus: 2,
                reactance_pu: 0.1,
                rating_mw: 100.0,
            },
        ];
        assert!(is_connected(3, &lines), "3-bus chain should be connected");
        assert!(
            !is_connected(4, &lines),
            "4-bus with 2 lines should be disconnected"
        );
    }
}