powerio_matrix/matrix/
incidence.rs

1//! DC network primitives: the signed incidence matrix `A`, branch
2//! susceptances `b`, the flow map `B Aᵀ`, and the phase shift injection.
3//!
4//! Edge orientation is fixed to MATPOWER's from→to: column `e` of `A` has
5//! `+1` at the from bus (tail) and `−1` at the to bus (head). Columns run
6//! over in-service branches in `case.branches` order; `branch_of_col` maps a
7//! column back to its source branch index.
8
9use sprs::CsMat;
10
11pub use powerio::DcConvention;
12
13use crate::Result;
14use crate::indexed::IndexedNetwork;
15use crate::matrix::triplet::CooBuilder;
16
17use super::{BuildOptions, ZeroImpedanceSkips};
18
19/// The incidence factorization of a case under one DC convention.
20#[derive(Debug, Clone)]
21#[non_exhaustive]
22pub struct IncidenceParts {
23    /// Signed incidence `A`, shape `n × m`.
24    pub a: CsMat<f64>,
25    /// Branch susceptances `b_e`, length `m`.
26    pub b: Vec<f64>,
27    /// Phase shift bus injection, length `n`. All zeros unless the MATPOWER
28    /// convention is used and shifters are present.
29    pub p_shift: Vec<f64>,
30    /// Column `k` → index into `case.branches`.
31    pub branch_of_col: Vec<usize>,
32    /// In-service branch rows skipped because their DC denominator is zero.
33    pub skipped_zero_impedance: ZeroImpedanceSkips,
34}
35
36impl IncidenceParts {
37    #[inline]
38    pub fn n(&self) -> usize {
39        self.a.rows()
40    }
41
42    #[inline]
43    pub fn m(&self) -> usize {
44        self.a.cols()
45    }
46}
47
48/// Build `A`, `b`, the phase shift injection, and the column→branch map.
49///
50/// Self-loops (from == to) are dropped. A branch whose reactance is too small
51/// to divide by has no DC susceptance the Laplacian can carry; it is skipped
52/// when `opts.skip_zero_impedance` is true and rejected with
53/// [`powerio::Error::ZeroImpedance`] when it is false. A tap ratio under the same bound
54/// is [`powerio::Error::DegenerateTap`] either way, as it is in Y_bus.
55pub fn build_incidence(
56    case: &IndexedNetwork,
57    conv: DcConvention,
58    opts: &BuildOptions,
59) -> Result<IncidenceParts> {
60    let n = case.n();
61
62    // Pass 1: resolve and filter, fixing the column order.
63    let mut cols: Vec<Column> = Vec::new();
64    let mut skipped_zero_impedance = Vec::new();
65    for (idx, br) in case.in_service_branches() {
66        let i = case.bus_index(br.from).ok_or(powerio::Error::UnknownBus {
67            bus_id: br.from,
68            element_index: idx,
69        })?;
70        let j = case.bus_index(br.to).ok_or(powerio::Error::UnknownBus {
71            bus_id: br.to,
72            element_index: idx,
73        })?;
74        // Zero impedance in every sense the builder can act on: `x = 1e-300`
75        // gives a finite `b = 1e300` that annihilates every real branch sharing
76        // a diagonal with it. Exact zero used to be the whole test.
77        let degenerate_x = br.x.abs() < crate::matrix::MIN_DIVISIBLE_MAGNITUDE;
78        if i == j || degenerate_x {
79            if i != j && degenerate_x {
80                if !opts.skip_zero_impedance {
81                    return Err(powerio::Error::ZeroImpedance { row: idx }.into());
82                }
83                skipped_zero_impedance.push(idx);
84            }
85            continue;
86        }
87        // `Matpower` divides the susceptance by the tap, so it is bounded here
88        // by the same rule Y_bus and the instance builders apply.
89        let b_e = conv.branch_susceptance(br.r, br.x, br.divisible_tap(idx)?);
90        // A NaN reactance slips past the guard above and poisons the whole
91        // Laplacian.
92        if !b_e.is_finite() {
93            return Err(powerio::Error::NonFiniteSusceptance { row: idx }.into());
94        }
95        // angle_radians, not to_radians: a normalized network's shift is
96        // already in radians, so converting again would double-scale it.
97        let shift_rad = if conv.includes_phase_shifts() {
98            case.angle_radians(br.shift)
99        } else {
100            0.0
101        };
102        cols.push(Column {
103            i,
104            j,
105            b_e,
106            shift_rad,
107            branch: idx,
108        });
109    }
110
111    // Pass 2: assemble.
112    let m = cols.len();
113    let mut a = CooBuilder::with_capacity_rect(n, m, 2 * m);
114    let mut b = Vec::with_capacity(m);
115    let mut p_shift = vec![0.0; n];
116    let mut branch_of_col = Vec::with_capacity(m);
117    for (k, col) in cols.iter().enumerate() {
118        a.add(col.i, k, 1.0);
119        a.add(col.j, k, -1.0);
120        b.push(col.b_e);
121        branch_of_col.push(col.branch);
122        if col.shift_rad != 0.0 {
123            // MATPOWER makeBdc: Pbusinj = (Cf − Ct)ᵀ (b ⊙ (−shift)). Column k
124            // of (Cf − Ct) is e_from − e_to.
125            p_shift[col.i] -= col.b_e * col.shift_rad;
126            p_shift[col.j] += col.b_e * col.shift_rad;
127        }
128    }
129
130    Ok(IncidenceParts {
131        a: a.finish_csr(),
132        b,
133        p_shift,
134        branch_of_col,
135        skipped_zero_impedance: ZeroImpedanceSkips::new(skipped_zero_impedance),
136    })
137}
138
139struct Column {
140    i: usize,
141    j: usize,
142    b_e: f64,
143    shift_rad: f64,
144    branch: usize,
145}
146
147/// Sparse diagonal matrix from `values` (square, `len × len`).
148pub fn diagonal(values: &[f64]) -> CsMat<f64> {
149    let n = values.len();
150    let mut d = CooBuilder::with_capacity(n, n);
151    for (k, &v) in values.iter().enumerate() {
152        d.add(k, k, v);
153    }
154    d.finish_csr()
155}
156
157/// `B = diag(b)`, shape `m × m`.
158pub fn susceptance_diag(b: &[f64]) -> CsMat<f64> {
159    diagonal(b)
160}
161
162/// The flow map `B Aᵀ`, shape `m × n`: `f = (B Aᵀ) θ`.
163pub fn build_flow_map(a: &CsMat<f64>, b: &[f64]) -> CsMat<f64> {
164    let d = susceptance_diag(b);
165    let at = a.transpose_view().to_csr();
166    &d * &at
167}
powerio_matrix/matrix/incidence.rs

powerio_matrix/matrix/
incidence.rs
