powerio_matrix/matrix/
incidence.rs

1//! Internal DC solver factors: the signed bus by branch incidence matrix and
2//! positive branch weights.
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 the accepted in-service branches in `case.branches` order.
7
8use sprs::CsMat;
9
10use powerio_tx::BranchSusceptanceFormula;
11
12use crate::Result;
13use crate::indexed::IndexedNetwork;
14use crate::matrix::triplet::CooBuilder;
15
16use super::BuildOptions;
17
18/// The incidence factorization used by sensitivity builders.
19#[derive(Debug, Clone)]
20#[non_exhaustive]
21pub(crate) struct IncidenceParts {
22    /// Signed incidence `A`, shape `n × m`.
23    pub(crate) a: CsMat<f64>,
24    /// Positive Laplacian edge weights `b_e`, length `m`: the factor weight
25    /// a sparse solver uses (`|b|` of the PowerModels susceptance). The
26    /// signed susceptance lives on the DC operator surface.
27    pub(crate) b: Vec<f64>,
28}
29
30impl IncidenceParts {
31    #[inline]
32    pub(crate) fn n(&self) -> usize {
33        self.a.rows()
34    }
35
36    #[inline]
37    pub(crate) fn m(&self) -> usize {
38        self.a.cols()
39    }
40}
41
42/// Build the internal bus by branch incidence factor and positive weights.
43///
44/// Self-loops (from == to) are dropped. A branch whose reactance is too small
45/// to divide by has no DC susceptance the Laplacian can carry; it is skipped
46/// when `opts.skip_zero_impedance` is true and rejected with
47/// [`powerio_tx::Error::ZeroImpedance`] when it is false. A tap ratio under the same bound
48/// is [`powerio_tx::Error::DegenerateTap`] either way, as it is in Y_bus.
49pub(crate) fn build_incidence(
50    case: &IndexedNetwork,
51    formula: BranchSusceptanceFormula,
52    opts: &BuildOptions,
53) -> Result<IncidenceParts> {
54    let n = case.n();
55
56    // Pass 1: resolve and filter, fixing the column order.
57    let mut cols: Vec<Column> = Vec::new();
58    for (idx, br) in case.in_service_branches() {
59        let i = case
60            .bus_index(br.from)
61            .ok_or(powerio_tx::Error::UnknownBus {
62                bus_id: br.from,
63                element_index: idx,
64            })?;
65        let j = case.bus_index(br.to).ok_or(powerio_tx::Error::UnknownBus {
66            bus_id: br.to,
67            element_index: idx,
68        })?;
69        // Zero impedance in every sense the selected formula can act on:
70        // `x = 1e-300` gives a finite weight of `1e300` that annihilates
71        // every real branch sharing a diagonal with it. The guard reads only
72        // the denominator the selected formula divides by, so a value the
73        // formula never reads cannot reject a branch (#324): the reciprocal
74        // rules bound the reactance, and the series formula bounds the whole
75        // impedance magnitude.
76        let degenerate = match formula {
77            BranchSusceptanceFormula::SeriesSusceptance => {
78                br.r.hypot(br.x) < crate::matrix::MIN_DIVISIBLE_MAGNITUDE
79            }
80            _ => br.x.abs() < crate::matrix::MIN_DIVISIBLE_MAGNITUDE,
81        };
82        if i == j || degenerate {
83            if i != j && degenerate && !opts.skip_zero_impedance {
84                return Err(powerio_tx::Error::ZeroImpedance { row: idx }.into());
85            }
86            continue;
87        }
88        // Only `TapAdjustedReactance` divides by the tap, so only it can be
89        // bounded or rejected by one (#324); the other formulas never read
90        // the value.
91        let tap = if formula.reads_tap() {
92            br.calc_divisible_tap(idx)?
93        } else {
94            1.0
95        };
96        // The incidence parts carry the internal positive factor weight (the
97        // Laplacian edge weight a sparse solver factors); public PowerModels
98        // sign results are the DC operator surface's to emit.
99        let b_e = formula.calc_solver_edge_weight(br.r, br.x, tap);
100        // A NaN reactance slips past the guard above and poisons the whole
101        // Laplacian.
102        if !b_e.is_finite() {
103            return Err(powerio_tx::Error::NonFiniteSusceptance { row: idx }.into());
104        }
105        cols.push(Column { i, j, b_e });
106    }
107
108    // Pass 2: assemble.
109    let m = cols.len();
110    let mut a = CooBuilder::with_capacity_rect(n, m, 2 * m);
111    let mut b = Vec::with_capacity(m);
112    for (k, col) in cols.iter().enumerate() {
113        a.add(col.i, k, 1.0);
114        a.add(col.j, k, -1.0);
115        b.push(col.b_e);
116    }
117
118    Ok(IncidenceParts {
119        a: a.finish_csr(),
120        b,
121    })
122}
123
124struct Column {
125    i: usize,
126    j: usize,
127    b_e: f64,
128}
129
130/// Sparse diagonal matrix from `values` (square, `len × len`).
131pub fn calc_diagonal(values: &[f64]) -> CsMat<f64> {
132    let n = values.len();
133    let mut d = CooBuilder::with_capacity(n, n);
134    for (k, &v) in values.iter().enumerate() {
135        d.add(k, k, v);
136    }
137    d.finish_csr()
138}
139
140/// `B = diag(b)`, shape `m × m`.
141pub fn calc_susceptance_diagonal(b: &[f64]) -> CsMat<f64> {
142    calc_diagonal(b)
143}
144
145/// Calculate the branch flow matrix used by solver preparations,
146/// `Diagonal(w) * Aᵀ`, shape `m × n`, where `w` contains positive factor
147/// weights.
148///
149/// A complete affine branch flow also adds `-b * shift`; problem
150/// preparations expose that term through
151/// `DcOpfPreparation::calc_branch_flow_offset`. In the public PowerModels sign
152/// spelling the same flow is `p_branch = -Bf va + b .* shift` with negated
153/// susceptances.
154pub(crate) fn calc_solver_branch_flow_matrix(
155    bus_branch_incidence: &CsMat<f64>,
156    susceptance_magnitude: &[f64],
157) -> CsMat<f64> {
158    let diagonal = calc_susceptance_diagonal(susceptance_magnitude);
159    let transpose = bus_branch_incidence.transpose_view().to_csr();
160    &diagonal * &transpose
161}
powerio_matrix/matrix/incidence.rs

powerio_matrix/matrix/
incidence.rs
