powerio_matrix/matrix/
mod.rs

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

powerio_matrix/matrix/
mod.rs
