powerio_matrix/matrix/
mod.rs

1//! Sparse matrix builders for MATPOWER cases.
2//!
3//! Sign convention: the susceptance matrix has the form `B = A diag(b) Aᵀ`
4//! with the node-by-edge incidence `A` (n×m) and per-edge weight `b_e = x/(r²+x²)`
5//! (see `bprime.rs` for the entry-level form). Resulting matrices satisfy positive
6//! diagonal, negative off-diagonal, `diag = sum of |off-diagonal|` — the
7//! positive (M-matrix) Laplacian form SDDM solvers expect.
8
9mod adjacency;
10mod bdoubleprime;
11mod bprime;
12pub mod incidence;
13// The DC-OPF interior-point operators are experimental and off by default,
14// built only under `--features kkt`, which the default build and the main CI
15// jobs skip.
16#[cfg(feature = "kkt")]
17pub mod kkt;
18mod lacpf;
19pub mod laplacian;
20pub mod opf;
21pub mod sensitivity;
22pub mod triplet;
23mod ybus;
24
25#[cfg(test)]
26mod tests;
27
28pub use adjacency::build_adjacency;
29pub use bdoubleprime::build_bdoubleprime;
30pub use bprime::build_bprime;
31pub use incidence::{
32    DcConvention, IncidenceParts, build_flow_map, build_incidence, susceptance_diag,
33};
34pub use lacpf::build_lacpf;
35pub use laplacian::{
36    GroundedIndexMap, build_weighted_laplacian, ground_at, ground_at_each, reference_indicator,
37    unit_vector,
38};
39pub use opf::{BusCosts, GenCosts, OpfInstance, Units, build_opf_instance};
40pub use sensitivity::{build_lodf, build_ptdf, build_ptdf_lodf};
41pub use ybus::{YbusParts, build_ybus};
42// Crate-internal: the gridfm columnar export reuses the per-branch admittance and
43// flow kernels so its branch table and Y_bus agree with `build_ybus` by construction.
44#[cfg(feature = "gridfm")]
45pub(crate) use ybus::{YbusFlags, branch_admittance, branch_flows};
46
47use sprs::CsMat;
48
49/// Which FDPF scheme to use for B'.
50///
51/// - `Bx`: B' uses `-x / (r² + x²)` (what most modern solvers do).
52/// - `Xb`: B' uses `-1 / x` (original Stott/Alsac 1974). Requires `x ≠ 0`.
53#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, serde::Serialize, serde::Deserialize)]
54#[non_exhaustive]
55pub enum Scheme {
56    #[default]
57    Bx,
58    Xb,
59}
60
61#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
62pub struct BuildOptions {
63    pub scheme: Scheme,
64    /// Apply tap ratios when building B″ and Y-bus. (B′ always ignores taps.)
65    pub include_taps: bool,
66    /// Apply phase shifts when building Y-bus. (B′/B″ always ignore shifts.)
67    pub include_shifts: bool,
68    /// Drop branches whose `r² + x² = 0` (true) or error out (false).
69    pub skip_zero_impedance: bool,
70}
71
72impl Default for BuildOptions {
73    fn default() -> Self {
74        Self {
75            scheme: Scheme::Bx,
76            include_taps: true,
77            include_shifts: true,
78            skip_zero_impedance: true,
79        }
80    }
81}
82
83/// Which branch denominator a matrix builder uses when deciding whether a branch
84/// can contribute a finite admittance or susceptance.
85#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
86#[non_exhaustive]
87pub enum ZeroImpedanceRule {
88    /// Full series impedance, `r² + x²`. Used by Y_bus, LACPF, B' in BX mode,
89    /// and B'' in XB mode.
90    Series,
91    /// Reactance only denominator, `x`. Used by DC incidence, B' in XB mode,
92    /// and B'' in BX mode after resistance is zeroed.
93    Reactance,
94}
95
96/// Branch rows skipped because the selected builder cannot represent a zero
97/// branch denominator.
98#[derive(Debug, Clone, Default, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
99#[non_exhaustive]
100pub struct ZeroImpedanceSkips {
101    pub count: usize,
102    pub branch_indices: Vec<usize>,
103}
104
105impl ZeroImpedanceSkips {
106    pub fn new(branch_indices: Vec<usize>) -> Self {
107        Self {
108            count: branch_indices.len(),
109            branch_indices,
110        }
111    }
112
113    pub fn is_empty(&self) -> bool {
114        self.count == 0
115    }
116}
117
118/// Count in-service branch rows the given builder rule will skip. This is the
119/// shared accounting used by matrix metadata and solver property regressions.
120pub fn skipped_zero_impedance(
121    case: &crate::indexed::IndexedNetwork,
122    rule: ZeroImpedanceRule,
123) -> ZeroImpedanceSkips {
124    let branch_indices = case
125        .in_service_branches()
126        .filter_map(|(row, br)| {
127            let zero = match rule {
128                ZeroImpedanceRule::Series => br.r * br.r + br.x * br.x == 0.0,
129                ZeroImpedanceRule::Reactance => br.x == 0.0,
130            };
131            zero.then_some(row)
132        })
133        .collect();
134    ZeroImpedanceSkips::new(branch_indices)
135}
136
137/// Common stats over a sparse matrix used by the TUI and `meta.json`.
138#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
139#[non_exhaustive]
140pub struct MatrixStats {
141    pub n: usize,
142    pub nnz: usize,
143    pub min_diag: f64,
144    pub max_diag: f64,
145    /// Smallest `D_ii - sum_j |O_ij|` across all rows. Negative means
146    /// the matrix is not diagonally dominant.
147    pub min_dd_margin: f64,
148    /// Whether all off-diagonals are ≤ 0 (M-matrix sign pattern).
149    pub m_matrix_sign: bool,
150    pub frobenius_norm: f64,
151    /// Branch rows skipped under `BuildOptions::skip_zero_impedance`.
152    #[serde(default)]
153    pub skipped_zero_impedance: usize,
154    /// Source branch row indices skipped under `BuildOptions::skip_zero_impedance`.
155    #[serde(default)]
156    pub skipped_zero_impedance_branches: Vec<usize>,
157}
158
159impl MatrixStats {
160    pub fn from_csr(a: &CsMat<f64>) -> Self {
161        let n = a.rows();
162        let mut min_diag = f64::INFINITY;
163        let mut max_diag = f64::NEG_INFINITY;
164        let mut min_dd = f64::INFINITY;
165        let mut m_sign = true;
166        let mut fro_sq = 0.0_f64;
167
168        for (row_idx, row) in a.outer_iterator().enumerate() {
169            let mut diag = 0.0_f64;
170            let mut off_abs = 0.0_f64;
171            for (col, &v) in row.iter() {
172                fro_sq += v * v;
173                if col == row_idx {
174                    diag = v;
175                } else {
176                    off_abs += v.abs();
177                    if v > 0.0 {
178                        m_sign = false;
179                    }
180                }
181            }
182            min_diag = min_diag.min(diag);
183            max_diag = max_diag.max(diag);
184            min_dd = min_dd.min(diag - off_abs);
185        }
186
187        Self {
188            n,
189            nnz: a.nnz(),
190            min_diag,
191            max_diag,
192            min_dd_margin: min_dd,
193            m_matrix_sign: m_sign,
194            frobenius_norm: fro_sq.sqrt(),
195            skipped_zero_impedance: 0,
196            skipped_zero_impedance_branches: Vec::new(),
197        }
198    }
199
200    #[must_use]
201    pub fn with_zero_impedance_skips(mut self, skips: ZeroImpedanceSkips) -> Self {
202        self.skipped_zero_impedance = skips.count;
203        self.skipped_zero_impedance_branches = skips.branch_indices;
204        self
205    }
206}
207
208/// Negate every stored value of a sparse matrix in place. Used where the input
209/// is owned and consumed straight away (B″ and the `YbusB` pipeline arm), so no
210/// clone of the structure is needed.
211pub(crate) fn negate_into(mut a: CsMat<f64>) -> CsMat<f64> {
212    a.data_mut().iter_mut().for_each(|v| *v = -*v);
213    a
214}
215
216/// Whether a matrix is SDDM (symmetric diagonally dominant M-matrix).
217/// Useful as a quick sanity check before feeding it to an SDDM solver.
218pub fn sddm_check(a: &CsMat<f64>) -> bool {
219    let stats = MatrixStats::from_csr(a);
220    stats.m_matrix_sign && stats.min_dd_margin >= -1e-12 && stats.min_diag > 0.0
221}
powerio_matrix/matrix/mod.rs

powerio_matrix/matrix/
mod.rs
