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
11use crate::indexed::IndexedNetwork;
12use crate::matrix::triplet::CooBuilder;
13use crate::{Error, Result};
14
15use super::{BuildOptions, ZeroImpedanceSkips};
16
17/// DC susceptance convention for `b_e` and the Laplacian.
18#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, serde::Serialize, serde::Deserialize)]
19#[non_exhaustive]
20pub enum DcConvention {
21    /// `b_e = 1/x_e`; taps and phase shifts ignored: the textbook DC power
22    /// flow weight.
23    #[default]
24    PaperPure,
25    /// `b_e = 1/(x_e·τ_e)` with a phase shift injection vector, matching
26    /// MATPOWER `makeBdc`.
27    Matpower,
28}
29
30/// The incidence factorization of a case under one DC convention.
31#[derive(Debug, Clone)]
32#[non_exhaustive]
33pub struct IncidenceParts {
34    /// Signed incidence `A`, shape `n × m`.
35    pub a: CsMat<f64>,
36    /// Branch susceptances `b_e`, length `m`.
37    pub b: Vec<f64>,
38    /// Phase shift bus injection, length `n`. All zeros unless the MATPOWER
39    /// convention is used and shifters are present.
40    pub p_shift: Vec<f64>,
41    /// Column `k` → index into `case.branches`.
42    pub branch_of_col: Vec<usize>,
43    /// In-service branch rows skipped because their DC denominator is zero.
44    pub skipped_zero_impedance: ZeroImpedanceSkips,
45}
46
47impl IncidenceParts {
48    #[inline]
49    pub fn n(&self) -> usize {
50        self.a.rows()
51    }
52
53    #[inline]
54    pub fn m(&self) -> usize {
55        self.a.cols()
56    }
57}
58
59/// Build `A`, `b`, the phase shift injection, and the column→branch map.
60///
61/// Self-loops (from == to) are dropped. Branches with `x == 0` have no finite DC
62/// susceptance; they are skipped when `opts.skip_zero_impedance` is true and
63/// rejected with [`Error::ZeroImpedance`] when it is false.
64pub fn build_incidence(
65    case: &IndexedNetwork,
66    conv: DcConvention,
67    opts: &BuildOptions,
68) -> Result<IncidenceParts> {
69    let n = case.n();
70
71    // Pass 1: resolve and filter, fixing the column order.
72    let mut cols: Vec<Column> = Vec::new();
73    let mut skipped_zero_impedance = Vec::new();
74    for (idx, br) in case.in_service_branches() {
75        let i = case.bus_index(br.from).ok_or(Error::UnknownBus {
76            bus_id: br.from,
77            element_index: idx,
78        })?;
79        let j = case.bus_index(br.to).ok_or(Error::UnknownBus {
80            bus_id: br.to,
81            element_index: idx,
82        })?;
83        if i == j || br.x == 0.0 {
84            if i != j && br.x == 0.0 {
85                if !opts.skip_zero_impedance {
86                    return Err(Error::ZeroImpedance { row: idx });
87                }
88                skipped_zero_impedance.push(idx);
89            }
90            continue;
91        }
92        let b_e = match conv {
93            DcConvention::PaperPure => 1.0 / br.x,
94            DcConvention::Matpower => 1.0 / (br.x * br.effective_tap()),
95        };
96        // A NaN reactance slips past the `x == 0.0` guard above, and a
97        // denormal `x` yields Inf; either poisons the whole Laplacian.
98        if !b_e.is_finite() {
99            return Err(Error::NonFiniteSusceptance { row: idx });
100        }
101        let shift_rad = match conv {
102            DcConvention::PaperPure => 0.0,
103            // angle_radians, not to_radians: a normalized network's shift is
104            // already in radians, so converting again would double-scale it.
105            DcConvention::Matpower => case.angle_radians(br.shift),
106        };
107        cols.push(Column {
108            i,
109            j,
110            b_e,
111            shift_rad,
112            branch: idx,
113        });
114    }
115
116    // Pass 2: assemble.
117    let m = cols.len();
118    let mut a = CooBuilder::with_capacity_rect(n, m, 2 * m);
119    let mut b = Vec::with_capacity(m);
120    let mut p_shift = vec![0.0; n];
121    let mut branch_of_col = Vec::with_capacity(m);
122    for (k, col) in cols.iter().enumerate() {
123        a.add(col.i, k, 1.0);
124        a.add(col.j, k, -1.0);
125        b.push(col.b_e);
126        branch_of_col.push(col.branch);
127        if col.shift_rad != 0.0 {
128            // MATPOWER makeBdc: Pbusinj = (Cf − Ct)ᵀ (b ⊙ (−shift)). Column k
129            // of (Cf − Ct) is e_from − e_to.
130            p_shift[col.i] -= col.b_e * col.shift_rad;
131            p_shift[col.j] += col.b_e * col.shift_rad;
132        }
133    }
134
135    Ok(IncidenceParts {
136        a: a.finish_csr(),
137        b,
138        p_shift,
139        branch_of_col,
140        skipped_zero_impedance: ZeroImpedanceSkips::new(skipped_zero_impedance),
141    })
142}
143
144struct Column {
145    i: usize,
146    j: usize,
147    b_e: f64,
148    shift_rad: f64,
149    branch: usize,
150}
151
152/// Sparse diagonal matrix from `values` (square, `len × len`).
153pub fn diagonal(values: &[f64]) -> CsMat<f64> {
154    let n = values.len();
155    let mut d = CooBuilder::with_capacity(n, n);
156    for (k, &v) in values.iter().enumerate() {
157        d.add(k, k, v);
158    }
159    d.finish_csr()
160}
161
162/// `B = diag(b)`, shape `m × m`.
163pub fn susceptance_diag(b: &[f64]) -> CsMat<f64> {
164    diagonal(b)
165}
166
167/// The flow map `B Aᵀ`, shape `m × n`: `f = (B Aᵀ) θ`.
168pub fn build_flow_map(a: &CsMat<f64>, b: &[f64]) -> CsMat<f64> {
169    let d = susceptance_diag(b);
170    let at = a.transpose_view().to_csr();
171    &d * &at
172}
powerio_matrix/matrix/incidence.rs

powerio_matrix/matrix/
incidence.rs
