qsim_solvers/
dc.rs

1//! DC Power Flow Solver
2//!
3//! Linear approximation assuming:
4//! - Voltage magnitudes = 1.0 p.u.
5//! - Resistance negligible (R << X)
6//! - Small angle differences (sin θ ≈ θ)
7//!
8//! Solves: P = B × θ
9
10use nalgebra::{DMatrix, DVector};
11use qsim_core::{CoreError, Result, Solver, SolverResult, StateStore, Topology};
12
13/// DC Power Flow Solver
14///
15/// Uses B-matrix method for fast linear power flow solution.
16#[derive(Debug, Clone)]
17pub struct DcPowerFlowSolver {
18    /// Tolerance for convergence check
19    pub tolerance: f64,
20}
21
22impl DcPowerFlowSolver {
23    /// Create a new DC power flow solver
24    pub fn new() -> Self {
25        Self { tolerance: 1e-6 }
26    }
27
28    /// Create with custom tolerance
29    pub fn with_tolerance(tolerance: f64) -> Self {
30        Self { tolerance }
31    }
32
33    /// Build the B matrix from topology
34    /// 
35    /// B[i][j] = -1/X[i][j] for connected buses
36    /// B[i][i] = sum of 1/X for all branches connected to bus i
37    pub fn build_b_matrix(&self, _topology: &Topology, state: &StateStore) -> DMatrix<f64> {
38        let n = state.bus_count();
39        let mut b_matrix = DMatrix::zeros(n, n);
40        
41        // TODO: Build from actual branch data
42        // For now, return identity matrix as placeholder
43        for i in 0..n {
44            b_matrix[(i, i)] = 1.0;
45        }
46        
47        b_matrix
48    }
49}
50
51impl Default for DcPowerFlowSolver {
52    fn default() -> Self {
53        Self::new()
54    }
55}
56
57impl Solver for DcPowerFlowSolver {
58    fn name(&self) -> &'static str {
59        "DC Power Flow"
60    }
61
62    fn solve(&self, topology: &Topology, state: &mut StateStore) -> Result<SolverResult> {
63        let n = state.bus_count();
64        
65        if n == 0 {
66            return Err(CoreError::SimulationError("No buses in network".into()));
67        }
68        
69        if topology.bus_count() != n {
70            return Err(CoreError::SimulationError(
71                "Topology and state bus count mismatch".into(),
72            ));
73        }
74
75        // Build B matrix
76        let b_matrix = self.build_b_matrix(topology, state);
77        
78        // Build power injection vector (excluding slack bus)
79        let p_vector = DVector::from_vec(state.active_power.clone());
80        
81        // Solve B × θ = P
82        // For DC power flow, we exclude the slack bus (reference angle = 0)
83        // Simplified: assume bus 0 is slack
84        
85        if n > 1 {
86            // Reduced system (exclude slack bus)
87            let b_reduced = b_matrix.view((1, 1), (n - 1, n - 1)).clone_owned();
88            let p_reduced = p_vector.rows(1, n - 1).clone_owned();
89            
90            // Solve linear system
91            match b_reduced.lu().solve(&p_reduced) {
92                Some(theta) => {
93                    // Update state with computed angles
94                    state.voltage_angle[0] = 0.0; // Slack bus reference
95                    for i in 0..theta.len() {
96                        state.voltage_angle[i + 1] = theta[i];
97                    }
98                }
99                None => {
100                    return Err(CoreError::SimulationError(
101                        "B matrix is singular".into(),
102                    ));
103                }
104            }
105        }
106        
107        Ok(SolverResult::converged(1, 0.0))
108    }
109}
110
111#[cfg(test)]
112mod tests {
113    use super::*;
114
115    #[test]
116    fn test_solver_creation() {
117        let solver = DcPowerFlowSolver::new();
118        assert_eq!(solver.name(), "DC Power Flow");
119    }
120
121    #[test]
122    fn test_empty_network() {
123        let solver = DcPowerFlowSolver::new();
124        let topology = Topology::new();
125        let mut state = StateStore::new(0);
126        
127        let result = solver.solve(&topology, &mut state);
128        assert!(result.is_err());
129    }
130}