powerio_matrix/
dc_operators.rs

1//! The public DC matrix operations over a [`DcPfInstance`].
2//!
3//! The instance stays matrix free; these separate operations project it. The
4//! public values carry PowerModels signs — the branch susceptance is
5//! `imag(inv(series impedance))` under the selected formula, negative
6//! for an inductive branch — and every axis carries stable element mappings:
7//! bus axes use [`BusId`] in analysis row order, and branch axes map each
8//! column to its source branch or three winding transformer winding. With
9//! voltage angles `va` in radians, the branch flow identity is
10//!
11//! ```text
12//! p_branch = -Bf * va + b .* shift
13//! ```
14//!
15//! and the nodal balance is `p_bus = -B * va + p_shift`, all in per unit on
16//! the network MVA base. The reference constrained linear system a solver
17//! factors is stated over the internal positive factor weights (the positive
18//! semidefinite Laplacian, the negation of the public susceptances), with the
19//! sign conversion confined to filling public results.
20//!
21//! Injections update in place from the instance's specifications; the network
22//! dependent matrices are built once and an operating point update never
23//! reconstructs them.
24
25use crate::{AnalysisBranchSource, SparseMatrix};
26use powerio_core::Error;
27use powerio_tx::{BranchSusceptanceFormula, BusId, IndexedNetwork};
28
29use powerio_prob::diagnostics::codes;
30use powerio_prob::{DcBusSpecification, DcPfInstance};
31
32/// The stable row identity a mapping row reports: the element uid when one
33/// exists, else `table:row`.
34fn row_identity(uid: Option<&str>, table: &str, row: usize) -> String {
35    uid.map_or_else(|| format!("{table}:{row}"), str::to_owned)
36}
37
38fn calc_incidence(bus_count: usize, endpoints: &[(usize, usize)]) -> SparseMatrix {
39    let mut incidence = crate::matrix::triplet::CooBuilder::new_rect(bus_count, endpoints.len());
40    for (column, &(from, to)) in endpoints.iter().enumerate() {
41        incidence.add(from, column, 1.0);
42        incidence.add(to, column, -1.0);
43    }
44    incidence.finish_csr()
45}
46
47/// The reference constrained linear system: the positive definite matrix a
48/// sparse solver factors, its right hand side, and the mapping from reduced
49/// rows back to bus rows.
50#[derive(Clone, Debug, PartialEq)]
51#[non_exhaustive]
52pub struct ReferenceConstrainedSystem {
53    /// The reference grounded positive semidefinite matrix, `n - r` square:
54    /// the internal positive factor weights, the negation of the public bus
55    /// susceptance matrix with the reference rows and columns removed.
56    pub matrix: SparseMatrix,
57    /// The right hand side: net injection minus phase shift injection at the
58    /// retained buses, plus coupling from eliminated reference buses, per unit.
59    pub rhs: Vec<f64>,
60    /// Reduced row to dense bus row.
61    pub retained_rows: Vec<usize>,
62}
63
64/// Build options for [`DcOperators`].
65///
66/// The default refuses a zero impedance branch, because it has no finite DC
67/// operator row. `skip_zero_impedance` drops such a branch from the operator
68/// axis instead and records its analysis row in
69/// [`DcOperators::skipped_branch_rows`], the same choice
70/// `BuildOptions::skip_zero_impedance` offers the admittance builders.
71#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
72#[non_exhaustive]
73pub struct DcOperatorOptions {
74    /// Drop a zero impedance branch instead of refusing the build.
75    pub skip_zero_impedance: bool,
76}
77
78impl DcOperatorOptions {
79    /// The default options: refuse a zero impedance branch.
80    #[must_use]
81    pub const fn new() -> Self {
82        Self {
83            skip_zero_impedance: false,
84        }
85    }
86
87    /// Set whether a zero impedance branch is dropped instead of refused.
88    #[must_use]
89    pub const fn with_skip_zero_impedance(mut self, skip: bool) -> Self {
90        self.skip_zero_impedance = skip;
91        self
92    }
93}
94
95/// DC matrix operations built once from an instance.
96#[derive(Clone, Debug)]
97pub struct DcOperators {
98    bus_ids: Vec<BusId>,
99    branch_identities: Vec<String>,
100    /// Operator column to the analysis branch row it was built from.
101    branch_rows: Vec<usize>,
102    /// Analysis branch rows dropped under `DcOperatorOptions::skip_zero_impedance`.
103    skipped_branch_rows: Vec<usize>,
104    options: DcOperatorOptions,
105    analysis_sources: Vec<AnalysisBranchSource>,
106    /// `n × m`, `+1` at the from bus and `-1` at the to bus of each branch.
107    incidence: SparseMatrix,
108    /// Public per branch susceptance, PowerModels signs.
109    branch_susceptance: Vec<f64>,
110    /// Per branch phase shift, radians (zero when the formula carries
111    /// no shift injections).
112    shift_radians: Vec<f64>,
113    /// Per column `(from row, to row)`, stored at build so no injection or
114    /// system fill rederives it from the incidence pattern.
115    endpoints: Vec<(usize, usize)>,
116    /// Net per unit injection at each bus from the instance specifications.
117    net_injection: Vec<f64>,
118    reference_rows: Vec<usize>,
119    /// The stated reference angle, radians, dense over every bus row.
120    /// Meaningful only at a row listed in `reference_rows`; every other row
121    /// holds zero and is never read.
122    reference_va_radians: Vec<f64>,
123    branch_susceptance_formula: BranchSusceptanceFormula,
124}
125
126impl DcOperators {
127    /// Build the operators. Zero impedance branches are preserved by the
128    /// instance and have no finite DC row, so they refuse the build until
129    /// resolved explicitly with
130    /// [`merge_zero_impedance_buses`](powerio_prob::merge_zero_impedance_buses); no
131    /// branch is ever skipped silently. Out of service branches and self
132    /// loops carry no operator column.
133    ///
134    /// # Errors
135    /// A zero impedance branch, a non-finite branch value, or a branch naming
136    /// an undeclared bus.
137    pub fn build(instance: &DcPfInstance) -> Result<Self, Error> {
138        Self::build_with(instance, &DcOperatorOptions::default())
139    }
140
141    /// Build the operators under explicit [`DcOperatorOptions`]. With
142    /// `skip_zero_impedance` set, a zero impedance branch is dropped from the
143    /// operator axis and listed by [`skipped_branch_rows`](Self::skipped_branch_rows)
144    /// instead of refusing the build; every other rule of [`build`](Self::build)
145    /// holds.
146    ///
147    /// # Errors
148    /// A zero impedance branch when it is not skipped, a non-finite branch
149    /// value, or a branch naming an undeclared bus.
150    // One pass over the branch table that fills every axis and operator
151    // column; splitting it would scatter the invariants the columns share.
152    #[expect(clippy::too_many_lines)]
153    pub fn build_with(instance: &DcPfInstance, options: &DcOperatorOptions) -> Result<Self, Error> {
154        let source = instance.network();
155        let view = IndexedNetwork::new(source);
156        let network = view.network();
157        let formula = instance.branch_susceptance_formula();
158        let base = network.base_mva();
159        let bus_ids: Vec<BusId> = network.buses().iter().map(|bus| bus.id).collect();
160        let row_of: std::collections::BTreeMap<BusId, usize> = bus_ids
161            .iter()
162            .enumerate()
163            .map(|(row, &id)| (id, row))
164            .collect();
165        let position_of = |bus: BusId| row_of.get(&bus).copied();
166
167        let mut branch_identities = Vec::new();
168        let mut branch_rows = Vec::new();
169        let mut skipped_branch_rows = Vec::new();
170        let mut active_analysis_sources = Vec::new();
171        let mut branch_susceptance = Vec::new();
172        let mut shift_radians = Vec::new();
173        let mut endpoints = Vec::new();
174        let analysis_sources = crate::opf::analysis_branch_sources(source);
175        for (row, branch) in network.branches().iter().enumerate() {
176            if !branch.in_service || branch.from == branch.to {
177                continue;
178            }
179            let identity = row_identity(branch.uid.as_deref(), "branches", row);
180            let (Some(from), Some(to)) = (position_of(branch.from), position_of(branch.to)) else {
181                return Err(Error::new(
182                    &codes::BUILD_INSTANCE_SHAPE_MISMATCH,
183                    format!(
184                        "branch `{identity}` names bus {} or {} the network does not declare",
185                        branch.from, branch.to
186                    ),
187                ));
188            };
189            // Only the tap-reading formula can be bounded by a tap (#324).
190            let tap = if formula.reads_tap() {
191                branch.calc_divisible_tap(row).map_err(|_| {
192                    Error::new(
193                        &codes::BUILD_OPERATOR_NOT_A_NUMBER,
194                        format!(
195                            "branch `{identity}` states a tap the selected formula cannot divide by"
196                        ),
197                    )
198                })?
199            } else {
200                1.0
201            };
202            // The same divisibility floor the other DC builders apply: a
203            // formally nonzero impedance below it yields a finite weight big
204            // enough to annihilate every real branch sharing a bus.
205            let degenerate = match formula {
206                BranchSusceptanceFormula::SeriesSusceptance => {
207                    branch.r.hypot(branch.x) < powerio_tx::dc::MIN_DIVISIBLE_MAGNITUDE
208                }
209                // Any formula that reads a reactance is bounded by it.
210                _ => branch.x.abs() < powerio_tx::dc::MIN_DIVISIBLE_MAGNITUDE,
211            };
212            if degenerate {
213                if options.skip_zero_impedance {
214                    skipped_branch_rows.push(row);
215                    continue;
216                }
217                return Err(Error::new(
218                    &codes::BUILD_OPERATOR_ZERO_IMPEDANCE,
219                    format!(
220                        "zero impedance branch `{identity}` has no finite DC operator row; resolve it explicitly with merge_zero_impedance_buses or build with skip_zero_impedance"
221                    ),
222                ));
223            }
224            let susceptance = formula.calc_branch_susceptance(branch.r, branch.x, tap);
225            if !susceptance.is_finite() {
226                return Err(Error::new(
227                    &codes::BUILD_OPERATOR_ZERO_IMPEDANCE,
228                    format!(
229                        "branch `{identity}` has no finite DC susceptance under the selected formula; resolve it explicitly with merge_zero_impedance_buses"
230                    ),
231                ));
232            }
233            let shift = branch_phase_shift_radians(formula, network.is_normalized(), branch.shift);
234            if !shift.is_finite() {
235                return Err(Error::new(
236                    &codes::BUILD_OPERATOR_NOT_A_NUMBER,
237                    format!("branch `{identity}` states a non-finite phase shift"),
238                ));
239            }
240            branch_identities.push(identity);
241            branch_rows.push(row);
242            active_analysis_sources.push(analysis_sources[row]);
243            branch_susceptance.push(susceptance);
244            shift_radians.push(shift);
245            endpoints.push((from, to));
246        }
247
248        let incidence = calc_incidence(bus_ids.len(), &endpoints);
249        let mut operators = Self {
250            bus_ids,
251            branch_identities,
252            branch_rows,
253            skipped_branch_rows,
254            options: *options,
255            analysis_sources: active_analysis_sources,
256            incidence,
257            branch_susceptance,
258            shift_radians,
259            endpoints,
260            net_injection: Vec::new(),
261            reference_rows: Vec::new(),
262            reference_va_radians: Vec::new(),
263            branch_susceptance_formula: formula,
264        };
265        operators.refresh_injections(instance, base)?;
266        Ok(operators)
267    }
268
269    /// Refresh the injection vectors and reference rows from the instance's
270    /// current specifications. The network dependent matrices are untouched:
271    /// an operating point update goes through here and reconstructs nothing.
272    ///
273    /// # Errors
274    /// A specification list whose length disagrees with the built bus axis.
275    pub fn update(&mut self, instance: &DcPfInstance) -> Result<(), Error> {
276        let base = instance.network().base_mva();
277        self.refresh_injections(instance, base)
278    }
279
280    fn refresh_injections(&mut self, instance: &DcPfInstance, base: f64) -> Result<(), Error> {
281        let source_bus_count = instance.network().buses().len();
282        if instance.specifications().len() != source_bus_count
283            || source_bus_count > self.bus_ids.len()
284        {
285            return Err(Error::new(
286                &codes::BUILD_INSTANCE_SHAPE_MISMATCH,
287                format!(
288                    "the instance states {} bus specifications for {} source buses; the operators were built over {} analysis buses",
289                    instance.specifications().len(),
290                    source_bus_count,
291                    self.bus_ids.len()
292                ),
293            ));
294        }
295        let mut net_injection = vec![0.0; self.bus_ids.len()];
296        let mut reference_va_radians = vec![0.0; self.bus_ids.len()];
297        let mut reference_rows = Vec::new();
298        for (row, specification) in instance.specifications().iter().enumerate() {
299            match *specification {
300                DcBusSpecification::NetActivePower { p_mw } => {
301                    net_injection[row] = p_mw / base;
302                }
303                DcBusSpecification::Reference { va_degrees } => {
304                    reference_rows.push(row);
305                    reference_va_radians[row] = va_degrees.to_radians();
306                }
307                _ => {}
308            }
309        }
310        self.net_injection = net_injection;
311        self.reference_rows = reference_rows;
312        self.reference_va_radians = reference_va_radians;
313        Ok(())
314    }
315
316    /// The selected branch susceptance formula.
317    #[must_use]
318    pub const fn branch_susceptance_formula(&self) -> BranchSusceptanceFormula {
319        self.branch_susceptance_formula
320    }
321
322    /// Dense bus row to bus id, the row mapping of every bus axis.
323    #[must_use]
324    pub fn bus_ids(&self) -> &[BusId] {
325        &self.bus_ids
326    }
327
328    /// Operator column to stable analysis branch identity.
329    #[must_use]
330    pub fn branch_identities(&self) -> &[String] {
331        &self.branch_identities
332    }
333
334    /// Operator column to the analysis branch row it represents: the
335    /// position in the network's branch table (three winding transformer
336    /// windings follow the branches). Out of service branches, self loops,
337    /// and skipped zero impedance branches have no column, so this is the
338    /// row selection every branch axis shares.
339    #[must_use]
340    pub fn branch_rows(&self) -> &[usize] {
341        &self.branch_rows
342    }
343
344    /// Analysis branch rows dropped under
345    /// [`DcOperatorOptions::skip_zero_impedance`], in table order. Empty
346    /// unless the option was set and a zero impedance branch existed.
347    #[must_use]
348    pub fn skipped_branch_rows(&self) -> &[usize] {
349        &self.skipped_branch_rows
350    }
351
352    /// The options the operators were built with.
353    #[must_use]
354    pub const fn options(&self) -> DcOperatorOptions {
355        self.options
356    }
357
358    /// Operator column to the source branch or three winding transformer
359    /// winding represented by that column.
360    #[must_use]
361    pub fn analysis_sources(&self) -> &[AnalysisBranchSource] {
362        &self.analysis_sources
363    }
364
365    /// The PowerModels incidence matrix `A`, `m × n`: each branch row has
366    /// `+1` at its from bus and `-1` at its to bus.
367    #[must_use]
368    pub fn calc_incidence_matrix(&self) -> SparseMatrix {
369        self.incidence.transpose_view().to_csr()
370    }
371
372    /// Calculate the per branch susceptances `b`, with PowerModels signs.
373    #[must_use]
374    pub fn calc_branch_susceptances(&self) -> &[f64] {
375        &self.branch_susceptance
376    }
377
378    /// Calculate the branch flow matrix `Bf = Diagonal(b) * A`, `m × n`.
379    #[must_use]
380    pub fn calc_branch_flow_matrix(&self) -> SparseMatrix {
381        let transpose = self.incidence.transpose_view().to_csr();
382        scale_rows(&transpose, &self.branch_susceptance)
383    }
384
385    /// Calculate the bus susceptance matrix `B = Aᵀ * Diagonal(b) * A`, `n × n`.
386    #[must_use]
387    pub fn calc_bus_susceptance_matrix(&self) -> SparseMatrix {
388        let bf = self.calc_branch_flow_matrix();
389        &self.incidence * &bf
390    }
391
392    /// The per bus net power injection the instance states, per unit on the
393    /// network MVA base, in bus row order.
394    #[must_use]
395    pub fn bus_power_injection(&self) -> &[f64] {
396        &self.net_injection
397    }
398
399    /// Calculate the branch phase shift injection `b .* shift` in branch
400    /// order.
401    #[must_use]
402    pub fn calc_branch_phase_shift_injection(&self) -> Vec<f64> {
403        self.branch_susceptance
404            .iter()
405            .zip(self.shift_radians.iter())
406            .map(|(&susceptance, &shift)| susceptance * shift)
407            .collect()
408    }
409
410    /// Calculate the bus phase shift injection
411    /// `p_shift = Aᵀ * (b .* shift)` in bus order.
412    #[must_use]
413    pub fn calc_bus_phase_shift_injection(&self) -> Vec<f64> {
414        let mut injection = vec![0.0; self.bus_ids.len()];
415        for (column, value) in self
416            .calc_branch_phase_shift_injection()
417            .into_iter()
418            .enumerate()
419        {
420            if value == 0.0 {
421                continue;
422            }
423            // Column `column` of A is `e_from - e_to`.
424            let (from, to) = self.endpoints(column);
425            injection[from] += value;
426            injection[to] -= value;
427        }
428        injection
429    }
430
431    /// Calculate DC active power flow in branch order from bus voltage angles
432    /// in radians: `p_branch = -Bf va + b .* shift`.
433    ///
434    /// # Errors
435    /// The voltage angle length differs from the bus axis length.
436    pub fn calc_branch_flow_dc(&self, voltage_angles: &[f64]) -> Result<Vec<f64>, Error> {
437        if voltage_angles.len() != self.bus_ids.len() {
438            return Err(Error::new(
439                &codes::BUILD_INSTANCE_SHAPE_MISMATCH,
440                format!(
441                    "voltage_angles has length {}; expected {} for the bus axis",
442                    voltage_angles.len(),
443                    self.bus_ids.len()
444                ),
445            ));
446        }
447        Ok(self
448            .endpoints
449            .iter()
450            .zip(self.branch_susceptance.iter())
451            .zip(self.shift_radians.iter())
452            .map(|((&(from, to), &susceptance), &shift)| {
453                -susceptance * (voltage_angles[from] - voltage_angles[to]) + susceptance * shift
454            })
455            .collect())
456    }
457
458    /// Calculate DC active power injection in bus order from bus voltage
459    /// angles in radians: `p_bus = -B va + p_shift`.
460    ///
461    /// # Errors
462    /// The voltage angle length differs from the bus axis length.
463    pub fn calc_bus_injection_dc(&self, voltage_angles: &[f64]) -> Result<Vec<f64>, Error> {
464        let branch_flows = self.calc_branch_flow_dc(voltage_angles)?;
465        let mut injections = vec![0.0; self.bus_ids.len()];
466        for (&(from, to), flow) in self.endpoints.iter().zip(branch_flows) {
467            injections[from] += flow;
468            injections[to] -= flow;
469        }
470        Ok(injections)
471    }
472
473    /// Calculate the reference constrained linear system over the internal positive
474    /// factor weights: the reference grounded positive semidefinite matrix
475    /// `L = -B` with reference rows and columns removed, and the right hand
476    /// side at the retained buses is `p - p_shift` (from `p = -B va +
477    /// p_shift`) plus the coupling carried in from every eliminated
478    /// reference bus's stated angle, radians, so `L_grounded va = rhs`
479    /// solves the stated problem with each reference bus fixed at the angle
480    /// it states rather than at zero. Sign conversion from the public
481    /// susceptances is confined to this fill.
482    ///
483    /// # Errors
484    /// An instance with no reference row.
485    pub fn calc_reference_constrained_system(&self) -> Result<ReferenceConstrainedSystem, Error> {
486        if self.reference_rows.is_empty() {
487            return Err(Error::new(
488                &codes::BUILD_INSTANCE_NO_REFERENCE_BUS,
489                "the instance states no reference bus to ground the system",
490            ));
491        }
492        let n = self.bus_ids.len();
493        let mut is_reference = vec![false; n];
494        for &row in &self.reference_rows {
495            is_reference[row] = true;
496        }
497        let mut reduced_of_full = vec![usize::MAX; n];
498        let mut retained_rows = Vec::with_capacity(n - self.reference_rows.len());
499        for (row, reduced) in reduced_of_full.iter_mut().enumerate() {
500            if !is_reference[row] {
501                *reduced = retained_rows.len();
502                retained_rows.push(row);
503            }
504        }
505
506        let mut matrix = crate::matrix::triplet::CooBuilder::new(retained_rows.len());
507        // A branch to an eliminated reference bus still carries an
508        // off-diagonal Laplacian entry; since its column is gone, that
509        // entry's contribution moves to the retained row's right hand side
510        // instead, carrying the reference bus's stated angle in.
511        let mut reference_coupling = vec![0.0; retained_rows.len()];
512        for (column, &(from, to)) in self.endpoint_table().iter().enumerate() {
513            // The positive factor weight is the negated public susceptance.
514            let weight = -self.branch_susceptance[column];
515            let (rf, rt) = (reduced_of_full[from], reduced_of_full[to]);
516            if rf != usize::MAX {
517                matrix.add(rf, rf, weight);
518            }
519            if rt != usize::MAX {
520                matrix.add(rt, rt, weight);
521            }
522            match (rf != usize::MAX, rt != usize::MAX) {
523                (true, true) => {
524                    matrix.add(rf, rt, -weight);
525                    matrix.add(rt, rf, -weight);
526                }
527                (true, false) => reference_coupling[rf] += weight * self.reference_va_radians[to],
528                (false, true) => {
529                    reference_coupling[rt] += weight * self.reference_va_radians[from];
530                }
531                (false, false) => {}
532            }
533        }
534        let shift_injection = self.calc_bus_phase_shift_injection();
535        let rhs = retained_rows
536            .iter()
537            .zip(reference_coupling.iter())
538            .map(|(&row, &coupling)| self.net_injection[row] - shift_injection[row] + coupling)
539            .collect();
540        Ok(ReferenceConstrainedSystem {
541            matrix: matrix.finish_csr(),
542            rhs,
543            retained_rows,
544        })
545    }
546
547    fn endpoints(&self, column: usize) -> (usize, usize) {
548        self.endpoints[column]
549    }
550
551    /// Endpoint rows per column, recovered from the incidence structure.
552    fn endpoint_table(&self) -> &[(usize, usize)] {
553        &self.endpoints
554    }
555}
556
557fn branch_phase_shift_radians(
558    formula: BranchSusceptanceFormula,
559    normalized: bool,
560    shift: f64,
561) -> f64 {
562    if !formula.includes_phase_shifts() {
563        0.0
564    } else if normalized {
565        shift
566    } else {
567        shift.to_radians()
568    }
569}
570
571/// `diag(values) * matrix`, scaling each row of a CSR matrix.
572fn scale_rows(matrix: &SparseMatrix, values: &[f64]) -> SparseMatrix {
573    let mut scaled = matrix.clone();
574    for (row, mut row_vec) in scaled.outer_iterator_mut().enumerate() {
575        for (_, entry) in row_vec.iter_mut() {
576            *entry *= values[row];
577        }
578    }
579    scaled
580}
powerio_matrix/dc_operators.rs

powerio_matrix/
dc_operators.rs
