powerio_matrix/matrix/
bprime.rs

1//! Fast Decoupled Power Flow (FDPF) B' matrix — shuntless.
2//!
3//! Per Stott & Alsac (1974), B' is the susceptance Laplacian with all
4//! shunts removed and tap ratios / phase shifts ignored:
5//!
6//! - Off-diagonal `B'_ij = -x / (r² + x²)`  (BX scheme; default)
7//!   or `B'_ij = -1 / x`              (XB scheme)
8//! - Diagonal     `B'_ii = sum_j |B'_ij|`
9//!
10//! Result: positive diag, negative off-diag, diag = sum of |off-diag| — the
11//! positive (M-matrix) Laplacian convention SDDM solvers expect.
12
13use sprs::CsMat;
14
15use crate::indexed::IndexedNetwork;
16use crate::{Error, Result};
17
18use super::{BuildOptions, Scheme, triplet::CooBuilder};
19
20pub fn build_bprime(case: &IndexedNetwork, opts: &BuildOptions) -> Result<CsMat<f64>> {
21    let n = case.n();
22    let mut coo = CooBuilder::with_capacity(n, 4 * case.branches().len() + n);
23
24    for (row_idx, br) in case.in_service_branches() {
25        let i = case.bus_index(br.from).ok_or(Error::UnknownBus {
26            bus_id: br.from,
27            element_index: row_idx,
28        })?;
29        let j = case.bus_index(br.to).ok_or(Error::UnknownBus {
30            bus_id: br.to,
31            element_index: row_idx,
32        })?;
33
34        let b_off = match opts.scheme {
35            Scheme::Bx => {
36                let denom = br.r * br.r + br.x * br.x;
37                if denom == 0.0 {
38                    if opts.skip_zero_impedance {
39                        continue;
40                    }
41                    return Err(Error::ZeroImpedance { row: row_idx });
42                }
43                -br.x / denom
44            }
45            Scheme::Xb => {
46                if br.x == 0.0 {
47                    if opts.skip_zero_impedance {
48                        continue;
49                    }
50                    return Err(Error::ZeroImpedance { row: row_idx });
51                }
52                -1.0 / br.x
53            }
54        };
55
56        // A NaN/Inf reactance (the MATPOWER tokenizer accepts `NaN`/`Inf`) slips
57        // past the `== 0.0` checks above and would write a non-finite entry that
58        // silently poisons MatrixStats / sddm_check. Reject it loudly instead.
59        if !b_off.is_finite() {
60            return Err(Error::NonFiniteSusceptance { row: row_idx });
61        }
62
63        if i == j {
64            // self-loop: contributes only as a shunt, has no place in B'
65            continue;
66        }
67
68        coo.add_sym(i, j, b_off);
69        coo.add(i, i, -b_off);
70        coo.add(j, j, -b_off);
71    }
72
73    Ok(coo.finish_csr())
74}
powerio_matrix/matrix/bprime.rs

powerio_matrix/matrix/
bprime.rs
