surge_network/network/
angle_reference.rs

1// SPDX-License-Identifier: LicenseRef-PolyForm-Noncommercial-1.0.0
2//! Shared output-angle reference conventions for power-flow solutions.
3
4use serde::{Deserialize, Serialize};
5use tracing::warn;
6
7use crate::network::model::Network;
8
9/// Reference angle convention for reported bus voltage angles.
10///
11/// This is an output/reporting convention, not a physical network property.
12/// Changing the angle reference shifts all bus angles in an island uniformly
13/// and therefore does not change branch flows or bus injections derived from
14/// angle differences.
15#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, Default)]
16pub enum AngleReference {
17    /// Preserve the initialized angle of the original reference bus.
18    #[default]
19    PreserveInitial,
20
21    /// Force the original reference bus angle to zero.
22    Zero,
23
24    /// Shift all angles so a weighted-average reference angle is zero.
25    Distributed(DistributedAngleWeight),
26}
27
28/// Weighting scheme for [`AngleReference::Distributed`].
29#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, Default)]
30#[serde(rename_all = "snake_case")]
31pub enum DistributedAngleWeight {
32    /// Weight each bus by its total active load (MW).
33    #[default]
34    LoadWeighted,
35    /// Weight each bus by its total online generation (MW).
36    GenerationWeighted,
37    /// Weight each bus by connected generator inertia (`H * MBASE`).
38    InertiaWeighted,
39}
40
41/// Apply an output angle reference convention to all buses in the network.
42pub fn apply_angle_reference(
43    angles: &mut [f64],
44    network: &Network,
45    reference_bus_idx: usize,
46    reference_angle0_rad: f64,
47    mode: AngleReference,
48) {
49    let bus_indices: Vec<usize> = (0..angles.len()).collect();
50    apply_angle_reference_subset(
51        angles,
52        network,
53        &bus_indices,
54        reference_bus_idx,
55        reference_angle0_rad,
56        mode,
57    );
58}
59
60/// Apply an output angle reference convention to one island/subset of buses.
61///
62/// The transformation is a uniform shift over `bus_indices`, so downstream
63/// branch flows within that connected component are unchanged.
64pub fn apply_angle_reference_subset(
65    angles: &mut [f64],
66    network: &Network,
67    bus_indices: &[usize],
68    reference_bus_idx: usize,
69    reference_angle0_rad: f64,
70    mode: AngleReference,
71) {
72    if bus_indices.is_empty() || reference_bus_idx >= angles.len() {
73        return;
74    }
75
76    let shift = match mode {
77        AngleReference::PreserveInitial => reference_angle0_rad - angles[reference_bus_idx],
78        AngleReference::Zero => -angles[reference_bus_idx],
79        AngleReference::Distributed(weight_mode) => {
80            match distributed_reference_angle(network, angles, bus_indices, weight_mode) {
81                Some(theta_ref) => -theta_ref,
82                None => {
83                    warn!(
84                        ?weight_mode,
85                        "distributed angle reference: all weights are zero, falling back to PreserveInitial"
86                    );
87                    reference_angle0_rad - angles[reference_bus_idx]
88                }
89            }
90        }
91    };
92
93    if shift.abs() <= 1e-12 {
94        return;
95    }
96    for &bus_idx in bus_indices {
97        if let Some(angle) = angles.get_mut(bus_idx) {
98            *angle += shift;
99        }
100    }
101}
102
103fn distributed_reference_angle(
104    network: &Network,
105    angles: &[f64],
106    bus_indices: &[usize],
107    weight_mode: DistributedAngleWeight,
108) -> Option<f64> {
109    let mut in_subset = vec![false; network.buses.len()];
110    for &bus_idx in bus_indices {
111        if let Some(flag) = in_subset.get_mut(bus_idx) {
112            *flag = true;
113        }
114    }
115
116    let bus_map = network.bus_index_map();
117    let mut weighted_angle_sum = 0.0;
118    let mut total_weight = 0.0;
119
120    match weight_mode {
121        DistributedAngleWeight::LoadWeighted => {
122            for load in &network.loads {
123                if !load.in_service {
124                    continue;
125                }
126                let Some(&bus_idx) = bus_map.get(&load.bus) else {
127                    continue;
128                };
129                if !in_subset[bus_idx] {
130                    continue;
131                }
132                let weight = load.active_power_demand_mw.abs();
133                if weight > 0.0 {
134                    weighted_angle_sum += weight * angles[bus_idx];
135                    total_weight += weight;
136                }
137            }
138        }
139        DistributedAngleWeight::GenerationWeighted => {
140            for generator in &network.generators {
141                if !generator.in_service {
142                    continue;
143                }
144                let Some(&bus_idx) = bus_map.get(&generator.bus) else {
145                    continue;
146                };
147                if !in_subset[bus_idx] {
148                    continue;
149                }
150                let weight = generator.p.abs();
151                if weight > 0.0 {
152                    weighted_angle_sum += weight * angles[bus_idx];
153                    total_weight += weight;
154                }
155            }
156        }
157        DistributedAngleWeight::InertiaWeighted => {
158            for generator in &network.generators {
159                if !generator.in_service {
160                    continue;
161                }
162                let Some(&bus_idx) = bus_map.get(&generator.bus) else {
163                    continue;
164                };
165                if !in_subset[bus_idx] {
166                    continue;
167                }
168                let weight = generator.h_inertia_s.unwrap_or(0.0) * generator.machine_base_mva;
169                if weight > 0.0 {
170                    weighted_angle_sum += weight * angles[bus_idx];
171                    total_weight += weight;
172                }
173            }
174        }
175    }
176
177    (total_weight > 1e-12).then_some(weighted_angle_sum / total_weight)
178}
179
180#[cfg(test)]
181mod tests {
182    use super::{AngleReference, DistributedAngleWeight, apply_angle_reference_subset};
183    use crate::network::{Bus, BusType, Generator, Load, Network};
184
185    fn base_network() -> Network {
186        Network {
187            buses: vec![
188                Bus {
189                    number: 1,
190                    bus_type: BusType::Slack,
191                    ..Bus::default()
192                },
193                Bus {
194                    number: 2,
195                    bus_type: BusType::PQ,
196                    ..Bus::default()
197                },
198                Bus {
199                    number: 3,
200                    bus_type: BusType::PV,
201                    ..Bus::default()
202                },
203            ],
204            loads: vec![
205                Load {
206                    bus: 1,
207                    active_power_demand_mw: 10.0,
208                    ..Load::default()
209                },
210                Load {
211                    bus: 2,
212                    active_power_demand_mw: 30.0,
213                    ..Load::default()
214                },
215            ],
216            generators: vec![
217                Generator {
218                    bus: 1,
219                    p: 50.0,
220                    ..Generator::default()
221                },
222                Generator {
223                    bus: 3,
224                    p: 150.0,
225                    h_inertia_s: Some(4.0),
226                    machine_base_mva: 100.0,
227                    ..Generator::default()
228                },
229            ],
230            ..Network::default()
231        }
232    }
233
234    #[test]
235    fn distributed_load_reference_zeroes_subset_weighted_mean() {
236        let network = base_network();
237        let mut angles = vec![0.3, -0.1, 0.8];
238        apply_angle_reference_subset(
239            &mut angles,
240            &network,
241            &[0, 1],
242            0,
243            0.0,
244            AngleReference::Distributed(DistributedAngleWeight::LoadWeighted),
245        );
246        let weighted_mean = (10.0 * angles[0] + 30.0 * angles[1]) / 40.0;
247        assert!(weighted_mean.abs() < 1e-12);
248        assert!((angles[2] - 0.8).abs() < 1e-12);
249    }
250
251    #[test]
252    fn distributed_generation_reference_uses_generation_weights() {
253        let network = base_network();
254        let mut angles = vec![0.2, 0.1, -0.4];
255        apply_angle_reference_subset(
256            &mut angles,
257            &network,
258            &[0, 2],
259            0,
260            0.0,
261            AngleReference::Distributed(DistributedAngleWeight::GenerationWeighted),
262        );
263        let weighted_mean = (50.0 * angles[0] + 150.0 * angles[2]) / 200.0;
264        assert!(weighted_mean.abs() < 1e-12);
265        assert!((angles[1] - 0.1).abs() < 1e-12);
266    }
267
268    #[test]
269    fn distributed_inertia_reference_uses_h_times_mbase() {
270        let network = base_network();
271        let mut angles = vec![0.1, 0.0, 0.7];
272        apply_angle_reference_subset(
273            &mut angles,
274            &network,
275            &[0, 2],
276            0,
277            0.0,
278            AngleReference::Distributed(DistributedAngleWeight::InertiaWeighted),
279        );
280        let weighted_mean = (4.0 * 100.0 * angles[2]) / (4.0 * 100.0);
281        assert!(weighted_mean.abs() < 1e-12);
282    }
283
284    #[test]
285    fn zero_weight_distributed_reference_falls_back_to_preserve_initial() {
286        let network = Network {
287            buses: vec![
288                Bus {
289                    number: 1,
290                    bus_type: BusType::Slack,
291                    ..Bus::default()
292                },
293                Bus {
294                    number: 2,
295                    bus_type: BusType::PQ,
296                    ..Bus::default()
297                },
298            ],
299            ..Network::default()
300        };
301        let mut angles = vec![0.6, -0.2];
302        apply_angle_reference_subset(
303            &mut angles,
304            &network,
305            &[0, 1],
306            0,
307            0.25,
308            AngleReference::Distributed(DistributedAngleWeight::LoadWeighted),
309        );
310        assert!((angles[0] - 0.25).abs() < 1e-12);
311    }
312}
surge_network/network/angle_reference.rs

surge_network/network/
angle_reference.rs
