powerio_matrix/matrix/
opf.rs

1//! DC-OPF instance data derived from the network's generators and their cost
2//! curves: cost, bounds, thermal limits, the generator→bus map, and nodal load.
3//!
4//! The paper treats generation as a nodal variable `p_g ∈ ℝⁿ`, so the
5//! canonical vectors here are bus-indexed (length `n`), formed by scattering
6//! the generator-space data through `C_g`. The generator-space vectors and
7//! `C_g` ride along so a MATPOWER-faithful per-generator formulation can be
8//! reconstructed exactly.
9
10use sprs::CsMat;
11
12use crate::indexed::IndexedNetwork;
13use crate::matrix::incidence::{IncidenceParts, diagonal};
14use crate::matrix::triplet::CooBuilder;
15use crate::{Error, Result};
16
17/// Unit system for the emitted quantities.
18#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, serde::Serialize, serde::Deserialize)]
19#[non_exhaustive]
20pub enum Units {
21    /// Power divided by `baseMVA`; cost coefficients scaled so the cost is a
22    /// function of per unit power (`q ← 2c₂·base²`, `c ← c₁·base`). Keeps the
23    /// whole instance dimensionally consistent with the per unit Laplacian. A
24    /// network from [`to_normalized`](powerio::Network::to_normalized) is already
25    /// per unit, so this leaves it unchanged (the scaling divisor is `1.0`).
26    #[default]
27    PerUnit,
28    /// Raw MATPOWER units: power in MW, cost in native `$·MWh⁻¹` coefficients.
29    Native,
30}
31
32/// Length-n bus-indexed cost and bound vectors (paper form). All share index
33/// space; each is zero at buses with no generator.
34#[derive(Debug, Clone)]
35#[non_exhaustive]
36pub struct BusCosts {
37    /// Quadratic cost diagonal `q`, `cost = ½ q p² + c p`.
38    pub q: Vec<f64>,
39    /// Linear cost `c`.
40    pub c: Vec<f64>,
41    /// Upper generation bound (summed over generators at the bus).
42    pub pmax: Vec<f64>,
43    /// Lower generation bound.
44    pub pmin: Vec<f64>,
45    /// Nodal load `p_d`.
46    pub p_d: Vec<f64>,
47}
48
49/// Generator-space provenance (length n_gen, in `C_g` column order).
50#[derive(Debug, Clone)]
51#[non_exhaustive]
52pub struct GenCosts {
53    pub q: Vec<f64>,
54    pub c: Vec<f64>,
55    pub pmax: Vec<f64>,
56    pub pmin: Vec<f64>,
57    /// Column `g` → index into the in-service generators.
58    pub gen_of_col: Vec<usize>,
59}
60
61/// Static DC-OPF instance data for a case.
62#[derive(Debug, Clone)]
63#[non_exhaustive]
64pub struct OpfInstance {
65    pub n: usize,
66    pub m: usize,
67    /// Bus-indexed cost/bounds/load (length n).
68    pub bus: BusCosts,
69    /// Generator-space provenance (length n_gen).
70    pub gen_costs: GenCosts,
71    /// Thermal limit `f̄` (`RATE_A`); `0` means unlimited per MATPOWER. Length m.
72    pub f_max: Vec<f64>,
73    /// Generator→bus incidence, `n × n_gen`, one `1` per column.
74    pub c_g: CsMat<f64>,
75}
76
77impl OpfInstance {
78    /// Number of in-service generators (the generator-space vector length).
79    #[must_use]
80    pub fn n_gen(&self) -> usize {
81        self.gen_costs.q.len()
82    }
83}
84
85/// Build the OPF instance. Errors with [`Error::NoGenerators`] if the case has
86/// no in-service generators, [`Error::MissingGenCost`] if a generator has no
87/// cost row, or [`Error::UnsupportedCostModel`] if its cost is present but not
88/// a polynomial of degree ≤ 2.
89pub fn build_opf_instance(
90    case: &IndexedNetwork,
91    incidence: &IncidenceParts,
92    units: Units,
93) -> Result<OpfInstance> {
94    let n = case.n();
95    let m = incidence.m();
96    // per_unit_base is 1.0 for an already-normalized network, so Units::PerUnit
97    // on a normalized case divides by 1 instead of scaling a second time.
98    let base = case.per_unit_base();
99
100    let p_scale = match units {
101        Units::PerUnit => 1.0 / base,
102        Units::Native => 1.0,
103    };
104    // Native cost is c₂p² + c₁p with p in MW. For per unit p (p_MW = base·p_pu),
105    // q_pu = 2c₂·base² and c_pu = c₁·base.
106    let (q_scale, c_scale) = match units {
107        Units::PerUnit => (base * base, base),
108        Units::Native => (1.0, 1.0),
109    };
110
111    let in_service: Vec<(usize, &crate::network::Generator)> = case.in_service_gens().collect();
112    let n_gen = in_service.len();
113    if n_gen == 0 {
114        return Err(Error::NoGenerators);
115    }
116
117    let mut q_gen = Vec::with_capacity(n_gen);
118    let mut c_gen = Vec::with_capacity(n_gen);
119    let mut pmax_gen = Vec::with_capacity(n_gen);
120    let mut pmin_gen = Vec::with_capacity(n_gen);
121    let mut gen_of_col = Vec::with_capacity(n_gen);
122    let mut cg = CooBuilder::with_capacity_rect(n, n_gen, n_gen);
123
124    for (col, &(gidx, generator)) in in_service.iter().enumerate() {
125        let bus = case.bus_index(generator.bus).ok_or(Error::UnknownBus {
126            bus_id: generator.bus,
127            element_index: gidx,
128        })?;
129        // Distinguish a genuinely absent cost from a present-but-unsupported
130        // one (piecewise model 1, or polynomial degree ≥ 3) so the error tells
131        // the truth about which case the file hit.
132        let cost = generator
133            .cost
134            .as_ref()
135            .ok_or(Error::MissingGenCost { gen_index: gidx })?;
136        let (q_raw, c_raw) = cost.quadratic().ok_or(Error::UnsupportedCostModel {
137            gen_index: gidx,
138            model: cost.model,
139            ncost: cost.ncost,
140        })?;
141        q_gen.push(q_raw * q_scale);
142        c_gen.push(c_raw * c_scale);
143        pmax_gen.push(generator.pmax * p_scale);
144        pmin_gen.push(generator.pmin * p_scale);
145        gen_of_col.push(gidx);
146        cg.add(bus, col, 1.0);
147    }
148    let c_g = cg.finish_csr();
149
150    let q_bus = project_gen_to_bus(&c_g, &q_gen);
151    let c_bus = project_gen_to_bus(&c_g, &c_gen);
152    let pmax_bus = project_gen_to_bus(&c_g, &pmax_gen);
153    let pmin_bus = project_gen_to_bus(&c_g, &pmin_gen);
154
155    let p_d: Vec<f64> = case.pd().iter().map(|&p| p * p_scale).collect();
156
157    let f_max: Vec<f64> = incidence
158        .branch_of_col
159        .iter()
160        .map(|&k| case.branches()[k].rate_a * p_scale)
161        .collect();
162
163    Ok(OpfInstance {
164        n,
165        m,
166        bus: BusCosts {
167            q: q_bus,
168            c: c_bus,
169            pmax: pmax_bus,
170            pmin: pmin_bus,
171            p_d,
172        },
173        gen_costs: GenCosts {
174            q: q_gen,
175            c: c_gen,
176            pmax: pmax_gen,
177            pmin: pmin_gen,
178            gen_of_col,
179        },
180        f_max,
181        c_g,
182    })
183}
184
185/// Scatter-sum a generator-space vector onto buses: `out = C_g v`. Buses with
186/// several generators get the sum; one generator per bus is exact.
187pub fn project_gen_to_bus(c_g: &CsMat<f64>, v: &[f64]) -> Vec<f64> {
188    let mut out = vec![0.0; c_g.rows()];
189    for (&val, (bus, g)) in c_g {
190        out[bus] += val * v[g];
191    }
192    out
193}
194
195/// `Q = diag(q)` as a sparse matrix — the quadratic-cost analog of
196/// [`susceptance_diag`](crate::matrix::susceptance_diag). Feeds the DC-OPF QP
197/// objective `½ pᵀ Q p + cᵀ p` consumed by the `kkt` interior-point operators.
198pub fn cost_quadratic_diag(q: &[f64]) -> CsMat<f64> {
199    diagonal(q)
200}
powerio_matrix/matrix/opf.rs

powerio_matrix/matrix/
opf.rs
