oxigrid/powerflow/
dc_powerflow.rs

1use crate::error::{OxiGridError, Result};
2use crate::network::PowerNetwork;
3use crate::powerflow::result::BranchFlow;
4use crate::powerflow::{PowerFlowConfig, PowerFlowResult, PowerFlowSolver};
5use nalgebra::DMatrix;
6use nalgebra::DVector;
7
8pub struct DcPowerFlowSolver;
9
10impl PowerFlowSolver for DcPowerFlowSolver {
11    fn solve(&self, network: &PowerNetwork, _config: &PowerFlowConfig) -> Result<PowerFlowResult> {
12        network.validate()?;
13
14        let n = network.bus_count();
15        let slack_idx = network.slack_bus_index()?;
16
17        let mut b_prime = DMatrix::<f64>::zeros(n, n);
18
19        for branch in &network.branches {
20            if !branch.status {
21                continue;
22            }
23            let i = network.bus_index(branch.from_bus)?;
24            let j = network.bus_index(branch.to_bus)?;
25            let tap = branch.effective_tap();
26
27            let bij = 1.0 / (branch.x * tap);
28            b_prime[(i, j)] -= bij;
29            b_prime[(j, i)] -= bij;
30            b_prime[(i, i)] += bij;
31            b_prime[(j, j)] += bij;
32        }
33
34        let non_slack: Vec<usize> = (0..n).filter(|&i| i != slack_idx).collect();
35        let m = non_slack.len();
36
37        let mut b_red = DMatrix::<f64>::zeros(m, m);
38        for (ri, &i) in non_slack.iter().enumerate() {
39            for (rj, &j) in non_slack.iter().enumerate() {
40                b_red[(ri, rj)] = b_prime[(i, j)];
41            }
42        }
43
44        let (p_sched, _) = network.net_injection();
45        let p_red = DVector::from_vec(non_slack.iter().map(|&i| p_sched[i]).collect::<Vec<_>>());
46
47        let lu = b_red.lu();
48        let theta_red = lu
49            .solve(&p_red)
50            .ok_or_else(|| OxiGridError::LinearAlgebra("B' matrix is singular".to_string()))?;
51
52        let mut v_ang = vec![0.0_f64; n];
53        for (ri, &i) in non_slack.iter().enumerate() {
54            v_ang[i] = theta_red[ri];
55        }
56
57        let v_mag = vec![1.0_f64; n];
58
59        // Compute branch flows (DC: P only)
60        let mut p_inj_bus = vec![0.0_f64; n];
61        let mut branch_flows = Vec::with_capacity(network.branches.len());
62
63        for (idx, branch) in network.branches.iter().enumerate() {
64            if !branch.status {
65                branch_flows.push(BranchFlow {
66                    branch_index: idx,
67                    from_bus: branch.from_bus,
68                    to_bus: branch.to_bus,
69                    p_from_mw: 0.0,
70                    q_from_mvar: 0.0,
71                    p_to_mw: 0.0,
72                    q_to_mvar: 0.0,
73                    p_loss_mw: 0.0,
74                    q_loss_mvar: 0.0,
75                    loading_pct: 0.0,
76                });
77                continue;
78            }
79            let i = network.bus_index(branch.from_bus)?;
80            let j = network.bus_index(branch.to_bus)?;
81            let tap = branch.effective_tap();
82            let p_ij = (v_ang[i] - v_ang[j]) / (branch.x * tap);
83            let p_ij_mw = p_ij * network.base_mva;
84
85            p_inj_bus[i] += p_ij_mw;
86            p_inj_bus[j] -= p_ij_mw;
87
88            let loading = if branch.rate_a > 0.0 {
89                p_ij_mw.abs() / branch.rate_a * 100.0
90            } else {
91                0.0
92            };
93
94            branch_flows.push(BranchFlow {
95                branch_index: idx,
96                from_bus: branch.from_bus,
97                to_bus: branch.to_bus,
98                p_from_mw: p_ij_mw,
99                q_from_mvar: 0.0,
100                p_to_mw: -p_ij_mw, // lossless in DC
101                q_to_mvar: 0.0,
102                p_loss_mw: 0.0,
103                q_loss_mvar: 0.0,
104                loading_pct: loading,
105            });
106        }
107
108        Ok(PowerFlowResult {
109            voltage_magnitude: v_mag,
110            voltage_angle: v_ang,
111            p_injected: p_inj_bus,
112            q_injected: vec![0.0; n],
113            branch_flows,
114            total_p_loss_mw: 0.0,
115            total_q_loss_mvar: 0.0,
116            converged: true,
117            iterations: 1,
118            max_mismatch: 0.0,
119        })
120    }
121}
122
123#[cfg(test)]
124mod tests {
125    use super::*;
126    use crate::network::branch::Branch;
127    use crate::network::bus::{Bus, BusType};
128    use crate::network::topology::Generator;
129
130    #[test]
131    fn test_dc_2bus() {
132        let mut net = PowerNetwork::new(100.0);
133        net.buses.push(Bus::new(1, BusType::Slack));
134        let mut bus2 = Bus::new(2, BusType::PQ);
135        bus2.pd = crate::units::Power(50.0);
136        net.buses.push(bus2);
137        net.branches.push(Branch {
138            from_bus: 1,
139            to_bus: 2,
140            r: 0.01,
141            x: 0.1,
142            b: 0.0,
143            rate_a: 100.0,
144            rate_b: 100.0,
145            rate_c: 100.0,
146            tap: 0.0,
147            shift: 0.0,
148            status: true,
149        });
150        net.generators.push(Generator {
151            bus_id: 1,
152            pg: 50.0,
153            qg: 0.0,
154            qmax: 999.0,
155            qmin: -999.0,
156            vg: 1.0,
157            mbase: 100.0,
158            status: true,
159            pmax: 999.0,
160            pmin: 0.0,
161        });
162
163        let config = PowerFlowConfig::default();
164        let result = DcPowerFlowSolver.solve(&net, &config).unwrap();
165        assert!(result.converged);
166        assert!((result.voltage_angle[0]).abs() < 1e-10);
167        // Branch flow should be ~50 MW
168        assert_eq!(result.branch_flows.len(), 1);
169        assert!((result.branch_flows[0].p_from_mw - 50.0).abs() < 1.0);
170    }
171}
oxigrid/powerflow/dc_powerflow.rs

oxigrid/powerflow/
dc_powerflow.rs
