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

powerio_matrix/matrix/
incidence.rs
