1use nalgebra::{DMatrix, DVector};
11use qsim_core::{CoreError, Result, Solver, SolverResult, StateStore, Topology};
12
13#[derive(Debug, Clone)]
17pub struct DcPowerFlowSolver {
18 pub tolerance: f64,
20}
21
22impl DcPowerFlowSolver {
23 pub fn new() -> Self {
25 Self { tolerance: 1e-6 }
26 }
27
28 pub fn with_tolerance(tolerance: f64) -> Self {
30 Self { tolerance }
31 }
32
33 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 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 let b_matrix = self.build_b_matrix(topology, state);
77
78 let p_vector = DVector::from_vec(state.active_power.clone());
80
81 if n > 1 {
86 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 match b_reduced.lu().solve(&p_reduced) {
92 Some(theta) => {
93 state.voltage_angle[0] = 0.0; 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}