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/// Common stats over a sparse matrix used by the TUI and `meta.json`.
84#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
85#[non_exhaustive]
86pub struct MatrixStats {
87    pub n: usize,
88    pub nnz: usize,
89    pub min_diag: f64,
90    pub max_diag: f64,
91    /// Smallest `D_ii - sum_j |O_ij|` across all rows. Negative means
92    /// the matrix is not diagonally dominant.
93    pub min_dd_margin: f64,
94    /// Whether all off-diagonals are ≤ 0 (M-matrix sign pattern).
95    pub m_matrix_sign: bool,
96    pub frobenius_norm: f64,
97}
98
99impl MatrixStats {
100    pub fn from_csr(a: &CsMat<f64>) -> Self {
101        let n = a.rows();
102        let mut min_diag = f64::INFINITY;
103        let mut max_diag = f64::NEG_INFINITY;
104        let mut min_dd = f64::INFINITY;
105        let mut m_sign = true;
106        let mut fro_sq = 0.0_f64;
107
108        for (row_idx, row) in a.outer_iterator().enumerate() {
109            let mut diag = 0.0_f64;
110            let mut off_abs = 0.0_f64;
111            for (col, &v) in row.iter() {
112                fro_sq += v * v;
113                if col == row_idx {
114                    diag = v;
115                } else {
116                    off_abs += v.abs();
117                    if v > 0.0 {
118                        m_sign = false;
119                    }
120                }
121            }
122            min_diag = min_diag.min(diag);
123            max_diag = max_diag.max(diag);
124            min_dd = min_dd.min(diag - off_abs);
125        }
126
127        Self {
128            n,
129            nnz: a.nnz(),
130            min_diag,
131            max_diag,
132            min_dd_margin: min_dd,
133            m_matrix_sign: m_sign,
134            frobenius_norm: fro_sq.sqrt(),
135        }
136    }
137}
138
139/// Negate every stored value of a sparse matrix in place. Used where the input
140/// is owned and consumed straight away (B″ and the `YbusB` pipeline arm), so no
141/// clone of the structure is needed.
142pub(crate) fn negate_into(mut a: CsMat<f64>) -> CsMat<f64> {
143    a.data_mut().iter_mut().for_each(|v| *v = -*v);
144    a
145}
146
147/// Whether a matrix is SDDM (symmetric diagonally dominant M-matrix).
148/// Useful as a quick sanity check before feeding it to an SDDM solver.
149pub fn sddm_check(a: &CsMat<f64>) -> bool {
150    let stats = MatrixStats::from_csr(a);
151    stats.m_matrix_sign && stats.min_dd_margin >= -1e-12 && stats.min_diag > 0.0
152}
powerio_matrix/matrix/mod.rs

powerio_matrix/matrix/
mod.rs
