powerio_matrix/matrix/
mod.rs

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

powerio_matrix/matrix/
mod.rs
