powerio_prob/
dc.rs

1use serde::{Deserialize, Serialize};
2
3use powerio::{BusId, DcConvention, IndexedNetwork};
4
5use crate::{Error, Result};
6
7use crate::{ReferenceBuses, limits, nodal};
8
9/// Unit system for power and generator cost data.
10#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
11#[non_exhaustive]
12pub enum Units {
13    /// Power is per unit. Cost coefficients are scaled for per unit power.
14    #[default]
15    PerUnit,
16    /// Power remains in the source unit, normally MW.
17    Native,
18}
19
20impl std::str::FromStr for Units {
21    type Err = String;
22
23    /// The one alias table for the bindings: `per-unit`/`perunit`/`pu` and
24    /// `native`, case insensitive, `-`/`_` ignored.
25    fn from_str(name: &str) -> std::result::Result<Self, Self::Err> {
26        match name.to_ascii_lowercase().replace(['-', '_'], "").as_str() {
27            "perunit" | "pu" => Ok(Units::PerUnit),
28            "native" => Ok(Units::Native),
29            other => Err(format!(
30                "unknown units `{other}`; expected \"per-unit\" or \"native\""
31            )),
32        }
33    }
34}
35
36impl Units {
37    /// `(power, admittance)` multipliers for source data on `base` MVA. MW
38    /// valued quantities (demand, bounds, limits, MW valued shunts) scale by
39    /// the first; per unit admittances and susceptances by the second.
40    pub(crate) fn power_scales(self, base: f64) -> (f64, f64) {
41        match self {
42            Self::PerUnit => (1.0 / base, 1.0),
43            Self::Native => (1.0, base),
44        }
45    }
46
47    /// `(quadratic, linear)` generator cost coefficient multipliers for the
48    /// same unit selection. The constant term never scales.
49    pub(crate) fn cost_scales(self, base: f64) -> (f64, f64) {
50        match self {
51            Self::PerUnit => (base * base, base),
52            Self::Native => (1.0, 1.0),
53        }
54    }
55}
56
57/// Options for DC OPF instance assembly.
58#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
59pub struct DcOpfOptions {
60    pub convention: DcConvention,
61    pub units: Units,
62    /// Skip non-self-loop branches with zero reactance. If false, assembly
63    /// returns [`powerio::Error::ZeroImpedance`].
64    pub skip_zero_impedance: bool,
65    /// Give a branch with no thermal rating the bound
66    /// [`Branch::synthesize_rate_a`](powerio::Branch::synthesize_rate_a)
67    /// states. If false, `rate_a <= 0` reaches `f_max` as zero, which reads as
68    /// unlimited.
69    pub synthesize_unrated_limits: bool,
70}
71
72impl Default for DcOpfOptions {
73    fn default() -> Self {
74        Self {
75            convention: DcConvention::default(),
76            units: Units::default(),
77            skip_zero_impedance: true,
78            synthesize_unrated_limits: false,
79        }
80    }
81}
82
83/// Generator data in generator column order.
84#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
85#[non_exhaustive]
86pub struct DcGeneratorData {
87    /// Generator column to dense bus index.
88    pub bus_of_gen: Vec<usize>,
89    /// Generator column to source generator row.
90    pub source_rows: Vec<usize>,
91    /// Quadratic objective diagonal in `0.5 * q * p^2 + c * p + c0`.
92    pub q: Vec<f64>,
93    /// Linear objective coefficient.
94    pub c: Vec<f64>,
95    /// Constant objective term. Unscaled in both unit systems: it carries no
96    /// power dimension. It does not move the argmin, but a consumer reporting
97    /// or comparing objective values needs it.
98    pub c0: Vec<f64>,
99    pub pmax: Vec<f64>,
100    pub pmin: Vec<f64>,
101}
102
103/// Branch data in active branch column order.
104#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
105#[non_exhaustive]
106pub struct DcBranchData {
107    pub from_bus: Vec<usize>,
108    pub to_bus: Vec<usize>,
109    /// Branch susceptance in the selected power unit per radian, positive for
110    /// an inductive branch.
111    pub b: Vec<f64>,
112    /// Phase shift in radians. Zero unless the convention carries phase shift
113    /// injections.
114    pub shift: Vec<f64>,
115    /// Thermal limit in the selected power unit. Zero means unlimited.
116    pub f_max: Vec<f64>,
117    /// Branch angle bounds in radians.
118    pub angle_min: Vec<f64>,
119    pub angle_max: Vec<f64>,
120    /// Branch column to source branch row.
121    pub source_rows: Vec<usize>,
122    /// Source branch rows omitted because their reactance was zero.
123    pub skipped_zero_impedance: Vec<usize>,
124}
125
126/// Generator data in dense bus order, aggregated over the generators at each
127/// bus. See [`DcOpfInstance::nodal_generator_data`].
128#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
129#[non_exhaustive]
130pub struct NodalGeneratorData {
131    pub q: Vec<f64>,
132    pub c: Vec<f64>,
133    pub c0: Vec<f64>,
134    pub pmax: Vec<f64>,
135    pub pmin: Vec<f64>,
136    /// Which buses host a generator. A bus without one has a zero range and a
137    /// zero cost, which a formulation must not read as a free generator.
138    pub has_gen: Vec<bool>,
139}
140
141/// Matrix free DC OPF input data.
142///
143/// A problem instance is complete numerical input for one problem family. It
144/// is separate from the source network, a matrix projection, a solver
145/// formulation, and a solution.
146#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
147#[non_exhaustive]
148pub struct DcOpfInstance {
149    pub name: String,
150    pub n_buses: usize,
151    pub n_source_generators: usize,
152    pub n_source_branches: usize,
153    pub base_mva: f64,
154    pub units: Units,
155    pub convention: DcConvention,
156    pub skip_zero_impedance: bool,
157    /// Dense bus index to external bus ID.
158    pub bus_ids: Vec<BusId>,
159    pub reference_buses: ReferenceBuses,
160    /// Nodal active demand in dense bus order.
161    pub p_d: Vec<f64>,
162    /// Nodal shunt conductance in dense bus order.
163    ///
164    /// The DC approximation holds the voltage magnitude at one per unit, so a
165    /// shunt draws the constant real power `g_s` and does not depend on the
166    /// angle. It belongs in the injection: the bus susceptance matrix keeps
167    /// zero row sums and carries no shunt. A nodal balance subtracts it
168    /// beside [`Self::p_d`], as MATPOWER `runpf` does.
169    pub g_s: Vec<f64>,
170    /// Nodal phase shift injection in dense bus order.
171    pub p_shift: Vec<f64>,
172    pub generators: DcGeneratorData,
173    pub branches: DcBranchData,
174}
175
176impl DcOpfInstance {
177    #[must_use]
178    pub fn n_generators(&self) -> usize {
179        self.generators.q.len()
180    }
181
182    #[must_use]
183    pub fn n_branches(&self) -> usize {
184        self.branches.b.len()
185    }
186
187    /// Project generator cost and bounds to bus space.
188    ///
189    /// The bounds at a bus are the sum of the generator bounds, which is the
190    /// range the bus total can reach. The cost curves at a bus combine by the
191    /// parallel rule `q = 1 / Σ(1/qᵢ)`, the curve that the least cost split of
192    /// the bus total follows. That combination is an approximation: it agrees
193    /// with generator space only while the split stays inside the bound of
194    /// each generator. A bus with one generator keeps that generator's own
195    /// coefficients.
196    #[must_use]
197    pub fn nodal_generator_data(&self) -> NodalGeneratorData {
198        let n = self.n_buses;
199        let generators = &self.generators;
200        let bus_of_gen = &generators.bus_of_gen;
201        let costs =
202            nodal::combine_costs(n, bus_of_gen, &generators.q, &generators.c, &generators.c0);
203        NodalGeneratorData {
204            q: costs.q,
205            c: costs.c,
206            c0: costs.c0,
207            pmax: nodal::sum_by_bus(n, bus_of_gen, &generators.pmax),
208            pmin: nodal::sum_by_bus(n, bus_of_gen, &generators.pmin),
209            has_gen: nodal::buses_with_generators(n, bus_of_gen),
210        }
211    }
212}
213
214/// Build a matrix free DC OPF instance from an indexed network.
215#[allow(clippy::too_many_lines)]
216pub fn build_dc_opf_instance(
217    case: &IndexedNetwork,
218    options: &DcOpfOptions,
219) -> Result<DcOpfInstance> {
220    case.check_reference_coverage()?;
221    case.network().check_base_mva()?;
222
223    let n_buses = case.n();
224    let base = case.per_unit_base();
225    let (p_scale, b_scale) = options.units.power_scales(base);
226    let thermal = limits::ThermalLimits {
227        synthesize_unrated: options.synthesize_unrated_limits,
228        power_scale: p_scale,
229        admittance_scale: b_scale,
230    };
231    let (q_scale, c_scale) = options.units.cost_scales(base);
232
233    let mut bus_of_gen = Vec::new();
234    let mut generator_rows = Vec::new();
235    let mut q = Vec::new();
236    let mut c = Vec::new();
237    let mut c0 = Vec::new();
238    let mut pmax = Vec::new();
239    let mut pmin = Vec::new();
240
241    for (source_row, generator) in case.in_service_gens() {
242        let bus = case
243            .bus_index(generator.bus)
244            .ok_or(powerio::Error::UnknownBus {
245                bus_id: generator.bus,
246                element_index: source_row,
247            })?;
248        let cost = generator
249            .cost
250            .as_ref()
251            .ok_or(powerio::Error::MissingGenCost {
252                gen_index: source_row,
253            })?;
254        let (q_raw, c_raw, c0_raw) = nodal::quadratic_terms(cost, source_row)?;
255        bus_of_gen.push(bus);
256        generator_rows.push(source_row);
257        q.push(q_raw * q_scale);
258        c.push(c_raw * c_scale);
259        c0.push(c0_raw);
260        pmax.push(generator.pmax * p_scale);
261        pmin.push(generator.pmin * p_scale);
262    }
263    if q.is_empty() {
264        return Err(Error::NoGenerators);
265    }
266
267    let mut from_bus = Vec::new();
268    let mut to_bus = Vec::new();
269    let mut b = Vec::new();
270    let mut shift = Vec::new();
271    let mut f_max = Vec::new();
272    let mut angle_min = Vec::new();
273    let mut angle_max = Vec::new();
274    let mut branch_rows = Vec::new();
275    let mut skipped_zero_impedance = Vec::new();
276    let mut p_shift = vec![0.0; n_buses];
277    // Dense bus order is the position order of `network().buses`.
278    let buses = &case.network().buses;
279
280    for (source_row, branch) in case.in_service_branches() {
281        let from = case
282            .bus_index(branch.from)
283            .ok_or(powerio::Error::UnknownBus {
284                bus_id: branch.from,
285                element_index: source_row,
286            })?;
287        let to = case
288            .bus_index(branch.to)
289            .ok_or(powerio::Error::UnknownBus {
290                bus_id: branch.to,
291                element_index: source_row,
292            })?;
293        if from == to {
294            // A self-loop carries no angle difference, so it contributes no
295            // DC flow, and its shift injection cancels at its own bus.
296            continue;
297        }
298        // The reactance the DC matrix builders bound, on the same rule: an
299        // `x = 1e-300` gives a finite `b = 1e300` that annihilates every real
300        // branch sharing a bus with it. Exact zero used to be the whole test.
301        if branch.x.abs() < powerio::dc::MIN_DIVISIBLE_MAGNITUDE {
302            if options.skip_zero_impedance {
303                skipped_zero_impedance.push(source_row);
304                continue;
305            }
306            return Err(powerio::Error::ZeroImpedance { row: source_row }.into());
307        }
308        let branch_b = options.convention.branch_susceptance(
309            branch.r,
310            branch.x,
311            branch.divisible_tap(source_row)?,
312        ) * b_scale;
313        if !branch_b.is_finite() {
314            return Err(powerio::Error::NonFiniteSusceptance { row: source_row }.into());
315        }
316        let shift_rad = if options.convention.includes_phase_shifts() {
317            case.angle_radians(branch.shift)
318        } else {
319            0.0
320        };
321        if shift_rad != 0.0 {
322            p_shift[from] -= branch_b * shift_rad;
323            p_shift[to] += branch_b * shift_rad;
324        }
325        let amin = case.angle_radians(branch.angmin);
326        let amax = case.angle_radians(branch.angmax);
327        from_bus.push(from);
328        to_bus.push(to);
329        b.push(branch_b);
330        shift.push(shift_rad);
331        f_max.push(thermal.of(branch, amin, amax, buses[from].vmax, buses[to].vmax));
332        angle_min.push(amin);
333        angle_max.push(amax);
334        branch_rows.push(source_row);
335    }
336
337    Ok(DcOpfInstance {
338        name: case.name().to_owned(),
339        n_buses,
340        n_source_generators: case.generators().len(),
341        n_source_branches: case.branches().len(),
342        base_mva: case.base_mva(),
343        units: options.units,
344        convention: options.convention,
345        skip_zero_impedance: options.skip_zero_impedance,
346        bus_ids: (0..n_buses).map(|index| case.bus_id(index)).collect(),
347        reference_buses: ReferenceBuses::new(case.reference_bus_indices()),
348        p_d: case.pd().iter().map(|value| value * p_scale).collect(),
349        g_s: case.gs().iter().map(|value| value * p_scale).collect(),
350        p_shift,
351        generators: DcGeneratorData {
352            bus_of_gen,
353            source_rows: generator_rows,
354            q,
355            c,
356            c0,
357            pmax,
358            pmin,
359        },
360        branches: DcBranchData {
361            from_bus,
362            to_bus,
363            b,
364            shift,
365            f_max,
366            angle_min,
367            angle_max,
368            source_rows: branch_rows,
369            skipped_zero_impedance,
370        },
371    })
372}
powerio_prob/dc.rs

powerio_prob/
dc.rs
