powerio_matrix/
lindist3flow_devices.rs

1//! Device, objective, and nodal-balance preparation for LinDist3Flow.
2
3use std::collections::BTreeMap;
4
5use num_complex::Complex64;
6use powerio_dist::{Configuration, DistLoad, DistLoadVoltageModel};
7use powerio_prob::{LinDist3FlowOpfInstance, ObjectiveTerm};
8
9use crate::{
10    ConnectionPowerMap, Error, LinDist3FlowNetworkData, Result, build_lindist3flow_network_data,
11    connection_power_map, cross_voltage_coefficients, winding_voltage_coefficients,
12};
13
14/// One coefficient on a squared-voltage variable.
15#[derive(Clone, Copy, Debug, PartialEq)]
16#[non_exhaustive]
17pub struct LinDist3FlowVoltageTerm {
18    pub node: usize,
19    pub coefficient: f64,
20}
21
22/// An affine scalar over the prepared squared-voltage variables.
23#[derive(Clone, Debug, Default, PartialEq)]
24#[non_exhaustive]
25pub struct LinDist3FlowAffineExpression {
26    pub constant: f64,
27    pub voltage_terms: Vec<LinDist3FlowVoltageTerm>,
28}
29
30/// Active and reactive affine power expressions.
31#[derive(Clone, Debug, Default, PartialEq)]
32#[non_exhaustive]
33pub struct LinDist3FlowComplexAffinePower {
34    pub active: LinDist3FlowAffineExpression,
35    pub reactive: LinDist3FlowAffineExpression,
36}
37
38/// One physical load and its reference-frozen terminal allocation.
39#[derive(Clone, Debug, PartialEq)]
40#[non_exhaustive]
41pub struct LinDist3FlowLoadData {
42    pub load: String,
43    pub source_load_row: usize,
44    pub terminal_nodes: Vec<usize>,
45    pub incidence: Vec<Vec<f64>>,
46    pub terminal_power_map: ConnectionPowerMap,
47    pub channel_power: Vec<LinDist3FlowComplexAffinePower>,
48}
49
50/// One affine terminal shunt-power row.
51#[derive(Clone, Debug, PartialEq)]
52#[non_exhaustive]
53pub struct LinDist3FlowShuntData {
54    pub shunt: String,
55    pub source_shunt_row: usize,
56    pub terminal_nodes: Vec<usize>,
57    pub terminal_power: Vec<LinDist3FlowComplexAffinePower>,
58}
59
60/// Bounds, objective coefficient, and SOC data for one dispatch channel.
61#[derive(Clone, Debug, PartialEq)]
62#[non_exhaustive]
63pub struct LinDist3FlowDispatchChannel {
64    pub active_min: Option<f64>,
65    pub active_max: Option<f64>,
66    pub reactive_min: Option<f64>,
67    pub reactive_max: Option<f64>,
68    pub apparent_power_limit: Option<f64>,
69    pub current_limit: Option<f64>,
70    /// Affine live `|D v|²` used by the rotated current cone.
71    pub squared_winding_voltage: LinDist3FlowAffineExpression,
72    pub reference_winding_voltage: f64,
73    /// Linear objective coefficient on active watts.
74    pub active_objective_coefficient: f64,
75}
76
77/// One dispatchable generator and its terminal allocation.
78#[derive(Clone, Debug, PartialEq)]
79#[non_exhaustive]
80pub struct LinDist3FlowGeneratorData {
81    pub generator: String,
82    pub source_generator_row: usize,
83    pub terminal_nodes: Vec<usize>,
84    pub incidence: Vec<Vec<f64>>,
85    pub terminal_power_map: ConnectionPowerMap,
86    pub channels: Vec<LinDist3FlowDispatchChannel>,
87}
88
89/// One voltage source. Its channel powers inject directly at corresponding
90/// terminals; the squared voltages are fixed in [`LinDist3FlowNetworkData`].
91#[derive(Clone, Debug, PartialEq)]
92#[non_exhaustive]
93pub struct LinDist3FlowSourceData {
94    pub source: String,
95    pub source_source_row: usize,
96    pub terminal_nodes: Vec<usize>,
97    pub active_objective_coefficient: Vec<f64>,
98}
99
100/// Semantic scalar variables used by nodal-balance rows.
101#[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
102#[non_exhaustive]
103pub enum LinDist3FlowVariable {
104    LineActive { line: usize, conductor: usize },
105    LineReactive { line: usize, conductor: usize },
106    GeneratorActive { generator: usize, channel: usize },
107    GeneratorReactive { generator: usize, channel: usize },
108    SourceActive { source: usize, channel: usize },
109    SourceReactive { source: usize, channel: usize },
110}
111
112/// One coefficient on a semantic decision variable.
113#[derive(Clone, Debug, PartialEq)]
114#[non_exhaustive]
115pub struct LinDist3FlowBalanceTerm {
116    pub variable: LinDist3FlowVariable,
117    pub coefficient: f64,
118}
119
120/// One equality expression constrained to zero.
121#[derive(Clone, Debug, Default, PartialEq)]
122#[non_exhaustive]
123pub struct LinDist3FlowBalanceEquation {
124    pub affine: LinDist3FlowAffineExpression,
125    pub variable_terms: Vec<LinDist3FlowBalanceTerm>,
126}
127
128/// Active and reactive KCL at one retained bus terminal.
129#[derive(Clone, Debug, PartialEq)]
130#[non_exhaustive]
131pub struct LinDist3FlowBalanceRow {
132    pub node: usize,
133    pub active: LinDist3FlowBalanceEquation,
134    pub reactive: LinDist3FlowBalanceEquation,
135}
136
137/// All device rows and the fully expanded nodal balances.
138#[derive(Clone, Debug, PartialEq)]
139#[non_exhaustive]
140pub struct LinDist3FlowDeviceData {
141    pub loads: Vec<LinDist3FlowLoadData>,
142    pub shunts: Vec<LinDist3FlowShuntData>,
143    pub generators: Vec<LinDist3FlowGeneratorData>,
144    pub sources: Vec<LinDist3FlowSourceData>,
145    pub balances: Vec<LinDist3FlowBalanceRow>,
146}
147
148/// Complete solver-neutral preparation for the currently supported L3F slice.
149#[derive(Clone, Debug, PartialEq)]
150#[non_exhaustive]
151pub struct LinDist3FlowPreparation {
152    pub network: LinDist3FlowNetworkData,
153    pub devices: LinDist3FlowDeviceData,
154}
155
156fn invalid(reason: impl Into<String>) -> Error {
157    Error::InvalidLinDist3FlowCoefficients {
158        reason: reason.into(),
159    }
160}
161
162fn node_key(bus: &str, terminal: &str) -> (String, String) {
163    (bus.to_ascii_lowercase(), terminal.to_owned())
164}
165
166/// Build the physical-channel incidence matrix for a typed connection.
167///
168/// # Errors
169/// The terminal and channel arities do not describe a supported grounded-wye,
170/// single-phase, or delta connection.
171pub fn lindist3flow_connection_incidence(
172    configuration: Configuration,
173    terminal_count: usize,
174    channel_count: usize,
175) -> Result<Vec<Vec<f64>>> {
176    match configuration {
177        Configuration::SinglePhase | Configuration::Delta
178            if terminal_count == 2 && channel_count == 1 =>
179        {
180            Ok(vec![vec![1.0, -1.0]])
181        }
182        Configuration::Wye | Configuration::SinglePhase if terminal_count == channel_count => {
183            Ok((0..terminal_count)
184                .map(|row| {
185                    (0..terminal_count)
186                        .map(|column| f64::from(row == column))
187                        .collect()
188                })
189                .collect())
190        }
191        Configuration::Delta if terminal_count == 3 && channel_count == 3 => Ok(vec![
192            vec![1.0, -1.0, 0.0],
193            vec![0.0, 1.0, -1.0],
194            vec![-1.0, 0.0, 1.0],
195        ]),
196        _ => Err(invalid(format!(
197            "{configuration:?} connection with {terminal_count} terminals and {channel_count} channels is unsupported"
198        ))),
199    }
200}
201
202fn resolve_nodes(
203    positions: &BTreeMap<(String, String), usize>,
204    bus: &str,
205    terminals: &[String],
206    family: &str,
207    name: &str,
208) -> Result<Vec<usize>> {
209    terminals
210        .iter()
211        .map(|terminal| {
212            positions
213                .get(&node_key(bus, terminal))
214                .copied()
215                .ok_or_else(|| {
216                    invalid(format!(
217                        "{family} `{name}` names unknown terminal `{bus}/{terminal}`"
218                    ))
219                })
220        })
221        .collect()
222}
223
224fn references(instance: &LinDist3FlowOpfInstance, nodes: &[usize]) -> Vec<Complex64> {
225    nodes
226        .iter()
227        .map(|&node| {
228            let reference = &instance.reference().voltages[node];
229            Complex64::from_polar(reference.magnitude, reference.angle)
230        })
231        .collect()
232}
233
234fn global_affine(
235    constant: f64,
236    scale: f64,
237    local: &crate::AffineScalarCoefficients,
238    nodes: &[usize],
239) -> LinDist3FlowAffineExpression {
240    LinDist3FlowAffineExpression {
241        constant: constant + scale * local.constant,
242        voltage_terms: local
243            .coefficients
244            .iter()
245            .zip(nodes)
246            .filter_map(|(coefficient, node)| {
247                let coefficient = scale * coefficient;
248                (coefficient.abs() > f64::EPSILON).then_some(LinDist3FlowVoltageTerm {
249                    node: *node,
250                    coefficient,
251                })
252            })
253            .collect(),
254    }
255}
256
257fn channel_value(values: &[f64], channel: usize, channels: usize, label: &str) -> Result<f64> {
258    let value = match values.len() {
259        1 => values[0],
260        length if length == channels => values[channel],
261        length => {
262            return Err(invalid(format!(
263                "{label} has length {length}, expected 1 or {channels}"
264            )));
265        }
266    };
267    if !value.is_finite() {
268        return Err(invalid(format!("{label} entry {channel} is non-finite")));
269    }
270    Ok(value)
271}
272
273fn load_channel_power(
274    load: &DistLoad,
275    incidence: &[Vec<f64>],
276    reference: &[Complex64],
277    terminal_nodes: &[usize],
278    channel: usize,
279) -> Result<LinDist3FlowComplexAffinePower> {
280    let channels = load.p_nom.len();
281    let p_nom = load.p_nom[channel];
282    let q_nom = load.q_nom[channel];
283    if !p_nom.is_finite() || !q_nom.is_finite() {
284        return Err(invalid(format!(
285            "load `{}` channel {channel} nominal power is non-finite",
286            load.name
287        )));
288    }
289    if matches!(
290        load.voltage_model,
291        DistLoadVoltageModel::ConstantPower { .. }
292    ) {
293        return Ok(LinDist3FlowComplexAffinePower {
294            active: LinDist3FlowAffineExpression {
295                constant: p_nom,
296                voltage_terms: Vec::new(),
297            },
298            reactive: LinDist3FlowAffineExpression {
299                constant: q_nom,
300                voltage_terms: Vec::new(),
301            },
302        });
303    }
304    let winding = winding_voltage_coefficients(&incidence[channel], reference)?;
305    let (v_nom, alpha_z, alpha_p, beta_z, beta_p) = match &load.voltage_model {
306        DistLoadVoltageModel::ConstantImpedance { v_nom } => (
307            channel_value(v_nom, channel, channels, "load v_nom")?,
308            1.0,
309            0.0,
310            1.0,
311            0.0,
312        ),
313        DistLoadVoltageModel::Zip {
314            v_nom,
315            alpha_z,
316            alpha_i,
317            alpha_p,
318            beta_z,
319            beta_i,
320            beta_p,
321        } => {
322            let alpha_i = channel_value(alpha_i, channel, channels, "load alpha_i")?;
323            let beta_i = channel_value(beta_i, channel, channels, "load beta_i")?;
324            if alpha_i.abs() > f64::EPSILON || beta_i.abs() > f64::EPSILON {
325                return Err(invalid(format!(
326                    "load `{}` channel {channel} has a nonzero current fraction",
327                    load.name
328                )));
329            }
330            (
331                channel_value(v_nom, channel, channels, "load v_nom")?,
332                channel_value(alpha_z, channel, channels, "load alpha_z")?,
333                channel_value(alpha_p, channel, channels, "load alpha_p")?,
334                channel_value(beta_z, channel, channels, "load beta_z")?,
335                channel_value(beta_p, channel, channels, "load beta_p")?,
336            )
337        }
338        _ => {
339            return Err(invalid(format!(
340                "load `{}` voltage model is outside the supported ZP slice",
341                load.name
342            )));
343        }
344    };
345    if !v_nom.is_finite() || v_nom <= 0.0 {
346        return Err(invalid(format!(
347            "load `{}` channel {channel} v_nom must be finite and positive",
348            load.name
349        )));
350    }
351    Ok(LinDist3FlowComplexAffinePower {
352        active: global_affine(
353            p_nom * alpha_p,
354            p_nom * alpha_z / v_nom.powi(2),
355            &winding,
356            terminal_nodes,
357        ),
358        reactive: global_affine(
359            q_nom * beta_p,
360            q_nom * beta_z / v_nom.powi(2),
361            &winding,
362            terminal_nodes,
363        ),
364    })
365}
366
367fn shunt_terminal_power(
368    name: &str,
369    g: &[Vec<f64>],
370    b: &[Vec<f64>],
371    reference: &[Complex64],
372    terminal_nodes: &[usize],
373    terminal: usize,
374) -> Result<LinDist3FlowComplexAffinePower> {
375    let n = terminal_nodes.len();
376    if g.len() != n || b.len() != n {
377        return Err(invalid(format!(
378            "shunt `{name}` admittance does not have {n} rows"
379        )));
380    }
381    let mut active = LinDist3FlowAffineExpression::default();
382    let mut reactive = LinDist3FlowAffineExpression::default();
383    for other in 0..n {
384        if g[terminal].len() != n || b[terminal].len() != n {
385            return Err(invalid(format!(
386                "shunt `{name}` admittance row {terminal} does not have {n} entries"
387            )));
388        }
389        let admittance = Complex64::new(g[terminal][other], b[terminal][other]);
390        if !admittance.re.is_finite() || !admittance.im.is_finite() {
391            return Err(invalid(format!(
392                "shunt `{name}` admittance entry ({terminal}, {other}) is non-finite"
393            )));
394        }
395        let cross = cross_voltage_coefficients(reference[terminal], reference[other])?;
396        let scale = admittance.conj();
397        let constant = scale * cross.constant;
398        let own = scale * cross.coefficient_phi;
399        let other_coefficient = scale * cross.coefficient_psi;
400        active.constant += constant.re;
401        reactive.constant += constant.im;
402        active.voltage_terms.push(LinDist3FlowVoltageTerm {
403            node: terminal_nodes[terminal],
404            coefficient: own.re,
405        });
406        reactive.voltage_terms.push(LinDist3FlowVoltageTerm {
407            node: terminal_nodes[terminal],
408            coefficient: own.im,
409        });
410        active.voltage_terms.push(LinDist3FlowVoltageTerm {
411            node: terminal_nodes[other],
412            coefficient: other_coefficient.re,
413        });
414        reactive.voltage_terms.push(LinDist3FlowVoltageTerm {
415            node: terminal_nodes[other],
416            coefficient: other_coefficient.im,
417        });
418    }
419    Ok(LinDist3FlowComplexAffinePower { active, reactive })
420}
421
422fn checked_vector_value(
423    values: Option<&[f64]>,
424    channel: usize,
425    channels: usize,
426    label: &str,
427) -> Result<Option<f64>> {
428    let Some(values) = values else {
429        return Ok(None);
430    };
431    if values.len() != channels {
432        return Err(invalid(format!(
433            "{label} has length {}, expected {channels}",
434            values.len()
435        )));
436    }
437    let value = values[channel];
438    if !value.is_finite() {
439        return Err(invalid(format!("{label} entry {channel} is non-finite")));
440    }
441    Ok(Some(value))
442}
443
444fn positive_limit(
445    values: Option<&[f64]>,
446    channel: usize,
447    channels: usize,
448    label: &str,
449) -> Result<Option<f64>> {
450    let value = checked_vector_value(values, channel, channels, label)?;
451    if value.is_some_and(|value| value <= 0.0) {
452        return Err(invalid(format!("{label} entry {channel} must be positive")));
453    }
454    Ok(value)
455}
456
457fn generator_bounds(
458    nominal: f64,
459    lower: Option<&[f64]>,
460    upper: Option<&[f64]>,
461    channel: usize,
462    channels: usize,
463    label: &str,
464    selected: bool,
465) -> Result<(Option<f64>, Option<f64>)> {
466    if !selected {
467        return Ok((None, None));
468    }
469    let lower = checked_vector_value(lower, channel, channels, &format!("{label} minimum"))?;
470    let upper = checked_vector_value(upper, channel, channels, &format!("{label} maximum"))?;
471    let (lower, upper) = if lower.is_none() && upper.is_none() {
472        if !nominal.is_finite() {
473            return Err(invalid(format!("{label} nominal value is non-finite")));
474        }
475        (Some(nominal), Some(nominal))
476    } else {
477        (lower, upper)
478    };
479    if lower.zip(upper).is_some_and(|(lower, upper)| lower > upper) {
480        return Err(invalid(format!("{label} bounds are inverted")));
481    }
482    Ok((lower, upper))
483}
484
485fn cost_value(values: Option<&[f64]>, channel: usize, channels: usize, label: &str) -> Result<f64> {
486    let Some(values) = values else {
487        return Ok(0.0);
488    };
489    Ok(channel_value(values, channel, channels, label)? / 1000.0)
490}
491
492fn uses_dispatch_cost(instance: &LinDist3FlowOpfInstance) -> Result<bool> {
493    let mut enabled = false;
494    for term in instance.base_instance().objective().terms() {
495        match term {
496            ObjectiveTerm::ActivePowerDispatchCost => enabled = true,
497            ObjectiveTerm::NetworkGeneratorCost => {
498                return Err(invalid(
499                    "LinDist3Flow does not compile balanced-network generator cost curves",
500                ));
501            }
502            _ => {
503                return Err(invalid(
504                    "the LinDist3Flow objective contains an unknown term",
505                ));
506            }
507        }
508    }
509    Ok(enabled)
510}
511
512fn add_affine(
513    target: &mut LinDist3FlowAffineExpression,
514    source: &LinDist3FlowAffineExpression,
515    scale: f64,
516) {
517    target.constant += scale * source.constant;
518    target
519        .voltage_terms
520        .extend(source.voltage_terms.iter().filter_map(|term| {
521            let coefficient = scale * term.coefficient;
522            (coefficient.abs() > f64::EPSILON).then_some(LinDist3FlowVoltageTerm {
523                node: term.node,
524                coefficient,
525            })
526        }));
527}
528
529fn add_variable(
530    equation: &mut LinDist3FlowBalanceEquation,
531    variable: LinDist3FlowVariable,
532    coefficient: f64,
533) {
534    if coefficient.abs() > f64::EPSILON {
535        equation.variable_terms.push(LinDist3FlowBalanceTerm {
536            variable,
537            coefficient,
538        });
539    }
540}
541
542/// Prepare load, shunt, generator, source, objective, and KCL rows against an
543/// already prepared L3F node/line axis.
544///
545/// Every balance equation is stated as an affine expression plus semantic
546/// decision-variable terms, all constrained to zero. This keeps row numbering
547/// and cone storage decisions in the solver adapter.
548///
549/// # Errors
550/// A device has an unsupported connection or voltage model, inconsistent
551/// dimensions, an invalid bound/rating, or an unresolved terminal identity.
552#[allow(clippy::too_many_lines)]
553pub fn build_lindist3flow_device_data(
554    instance: &LinDist3FlowOpfInstance,
555    network_data: &LinDist3FlowNetworkData,
556) -> Result<LinDist3FlowDeviceData> {
557    let network = instance.network();
558    let positions = network_data
559        .nodes
560        .iter()
561        .enumerate()
562        .map(|(position, node)| (node_key(&node.node.bus, &node.node.terminal), position))
563        .collect::<BTreeMap<_, _>>();
564    let dispatch_cost = uses_dispatch_cost(instance)?;
565
566    let mut loads = Vec::with_capacity(network.loads().len());
567    for (row, load) in network.loads().iter().enumerate() {
568        if load.p_nom.len() != load.q_nom.len() || load.p_nom.is_empty() {
569            return Err(invalid(format!(
570                "load `{}` active/reactive channel counts are empty or unequal",
571                load.name
572            )));
573        }
574        let terminal_nodes = resolve_nodes(
575            &positions,
576            &load.bus,
577            &load.terminal_map,
578            "load",
579            &load.name,
580        )?;
581        let incidence = lindist3flow_connection_incidence(
582            load.configuration,
583            terminal_nodes.len(),
584            load.p_nom.len(),
585        )?;
586        let reference = references(instance, &terminal_nodes);
587        let terminal_power_map = connection_power_map(&incidence, &reference)?;
588        let channel_power = (0..load.p_nom.len())
589            .map(|channel| {
590                load_channel_power(load, &incidence, &reference, &terminal_nodes, channel)
591            })
592            .collect::<Result<Vec<_>>>()?;
593        loads.push(LinDist3FlowLoadData {
594            load: load.name.clone(),
595            source_load_row: row,
596            terminal_nodes,
597            incidence,
598            terminal_power_map,
599            channel_power,
600        });
601    }
602
603    let mut shunts = Vec::with_capacity(network.shunts().len());
604    for (row, shunt) in network.shunts().iter().enumerate() {
605        let terminal_nodes = resolve_nodes(
606            &positions,
607            &shunt.bus,
608            &shunt.terminal_map,
609            "shunt",
610            &shunt.name,
611        )?;
612        let reference = references(instance, &terminal_nodes);
613        let terminal_power = (0..terminal_nodes.len())
614            .map(|terminal| {
615                shunt_terminal_power(
616                    &shunt.name,
617                    &shunt.g,
618                    &shunt.b,
619                    &reference,
620                    &terminal_nodes,
621                    terminal,
622                )
623            })
624            .collect::<Result<Vec<_>>>()?;
625        shunts.push(LinDist3FlowShuntData {
626            shunt: shunt.name.clone(),
627            source_shunt_row: row,
628            terminal_nodes,
629            terminal_power,
630        });
631    }
632
633    let mut generators = Vec::with_capacity(network.generators().len());
634    for (row, generator) in network.generators().iter().enumerate() {
635        let channels = generator.p_nom.len();
636        if channels == 0 || generator.q_nom.len() != channels {
637            return Err(invalid(format!(
638                "generator `{}` active/reactive channel counts are empty or unequal",
639                generator.name
640            )));
641        }
642        let terminal_nodes = resolve_nodes(
643            &positions,
644            &generator.bus,
645            &generator.terminal_map,
646            "generator",
647            &generator.name,
648        )?;
649        let incidence = lindist3flow_connection_incidence(
650            generator.configuration,
651            terminal_nodes.len(),
652            channels,
653        )?;
654        let reference = references(instance, &terminal_nodes);
655        let terminal_power_map = connection_power_map(&incidence, &reference)?;
656        let capability_selected = instance
657            .base_instance()
658            .constraints()
659            .generator_capability
660            .selects(&generator.name);
661        let mut prepared_channels = Vec::with_capacity(channels);
662        for (channel, incidence_row) in incidence.iter().enumerate() {
663            let winding = winding_voltage_coefficients(incidence_row, &reference)?;
664            let squared_winding_voltage = global_affine(0.0, 1.0, &winding, &terminal_nodes);
665            let reference_winding_voltage =
666                terminal_power_map.reference_winding_voltage[channel].norm();
667            let (active_min, active_max) = generator_bounds(
668                generator.p_nom[channel],
669                generator.p_min.as_deref(),
670                generator.p_max.as_deref(),
671                channel,
672                channels,
673                "generator active power",
674                capability_selected,
675            )?;
676            let (reactive_min, reactive_max) = generator_bounds(
677                generator.q_nom[channel],
678                generator.q_min.as_deref(),
679                generator.q_max.as_deref(),
680                channel,
681                channels,
682                "generator reactive power",
683                capability_selected,
684            )?;
685            prepared_channels.push(LinDist3FlowDispatchChannel {
686                active_min,
687                active_max,
688                reactive_min,
689                reactive_max,
690                apparent_power_limit: capability_selected
691                    .then(|| {
692                        positive_limit(
693                            generator.s_max.as_deref(),
694                            channel,
695                            channels,
696                            "generator apparent-power limit",
697                        )
698                    })
699                    .transpose()?
700                    .flatten(),
701                current_limit: capability_selected
702                    .then(|| {
703                        positive_limit(
704                            generator.i_max.as_deref(),
705                            channel,
706                            channels,
707                            "generator current limit",
708                        )
709                    })
710                    .transpose()?
711                    .flatten(),
712                squared_winding_voltage,
713                reference_winding_voltage,
714                active_objective_coefficient: if dispatch_cost {
715                    cost_value(
716                        generator.cost.as_deref(),
717                        channel,
718                        channels,
719                        "generator cost",
720                    )?
721                } else {
722                    0.0
723                },
724            });
725        }
726        generators.push(LinDist3FlowGeneratorData {
727            generator: generator.name.clone(),
728            source_generator_row: row,
729            terminal_nodes,
730            incidence,
731            terminal_power_map,
732            channels: prepared_channels,
733        });
734    }
735
736    let mut sources = Vec::with_capacity(network.sources().len());
737    for (row, source) in network.sources().iter().enumerate() {
738        let terminal_nodes = resolve_nodes(
739            &positions,
740            &source.bus,
741            &source.terminal_map,
742            "voltage source",
743            &source.name,
744        )?;
745        let channels = terminal_nodes.len();
746        let active_objective_coefficient = (0..channels)
747            .map(|channel| {
748                if dispatch_cost {
749                    cost_value(
750                        source.energy_cost_rate.as_deref(),
751                        channel,
752                        channels,
753                        "voltage source energy cost",
754                    )
755                } else {
756                    Ok(0.0)
757                }
758            })
759            .collect::<Result<Vec<_>>>()?;
760        sources.push(LinDist3FlowSourceData {
761            source: source.name.clone(),
762            source_source_row: row,
763            terminal_nodes,
764            active_objective_coefficient,
765        });
766    }
767
768    let mut balances = (0..network_data.nodes.len())
769        .map(|node| LinDist3FlowBalanceRow {
770            node,
771            active: LinDist3FlowBalanceEquation::default(),
772            reactive: LinDist3FlowBalanceEquation::default(),
773        })
774        .collect::<Vec<_>>();
775    for (line_index, line) in network_data.lines.iter().enumerate() {
776        for conductor in 0..line.parent_nodes.len() {
777            add_variable(
778                &mut balances[line.parent_nodes[conductor]].active,
779                LinDist3FlowVariable::LineActive {
780                    line: line_index,
781                    conductor,
782                },
783                -1.0,
784            );
785            add_variable(
786                &mut balances[line.parent_nodes[conductor]].reactive,
787                LinDist3FlowVariable::LineReactive {
788                    line: line_index,
789                    conductor,
790                },
791                -1.0,
792            );
793            add_variable(
794                &mut balances[line.child_nodes[conductor]].active,
795                LinDist3FlowVariable::LineActive {
796                    line: line_index,
797                    conductor,
798                },
799                1.0,
800            );
801            add_variable(
802                &mut balances[line.child_nodes[conductor]].reactive,
803                LinDist3FlowVariable::LineReactive {
804                    line: line_index,
805                    conductor,
806                },
807                1.0,
808            );
809        }
810    }
811    for load in &loads {
812        for terminal in 0..load.terminal_nodes.len() {
813            let balance = &mut balances[load.terminal_nodes[terminal]];
814            for channel in 0..load.channel_power.len() {
815                let power = &load.channel_power[channel];
816                let real = load.terminal_power_map.real_part[terminal][channel];
817                let imag = load.terminal_power_map.imag_part[terminal][channel];
818                add_affine(&mut balance.active.affine, &power.active, -real);
819                add_affine(&mut balance.active.affine, &power.reactive, imag);
820                add_affine(&mut balance.reactive.affine, &power.active, -imag);
821                add_affine(&mut balance.reactive.affine, &power.reactive, -real);
822            }
823        }
824    }
825    for shunt in &shunts {
826        for (terminal, &node) in shunt.terminal_nodes.iter().enumerate() {
827            add_affine(
828                &mut balances[node].active.affine,
829                &shunt.terminal_power[terminal].active,
830                -1.0,
831            );
832            add_affine(
833                &mut balances[node].reactive.affine,
834                &shunt.terminal_power[terminal].reactive,
835                -1.0,
836            );
837        }
838    }
839    for (generator_index, generator) in generators.iter().enumerate() {
840        for terminal in 0..generator.terminal_nodes.len() {
841            let balance = &mut balances[generator.terminal_nodes[terminal]];
842            for channel in 0..generator.channels.len() {
843                let real = generator.terminal_power_map.real_part[terminal][channel];
844                let imag = generator.terminal_power_map.imag_part[terminal][channel];
845                add_variable(
846                    &mut balance.active,
847                    LinDist3FlowVariable::GeneratorActive {
848                        generator: generator_index,
849                        channel,
850                    },
851                    real,
852                );
853                add_variable(
854                    &mut balance.active,
855                    LinDist3FlowVariable::GeneratorReactive {
856                        generator: generator_index,
857                        channel,
858                    },
859                    -imag,
860                );
861                add_variable(
862                    &mut balance.reactive,
863                    LinDist3FlowVariable::GeneratorActive {
864                        generator: generator_index,
865                        channel,
866                    },
867                    imag,
868                );
869                add_variable(
870                    &mut balance.reactive,
871                    LinDist3FlowVariable::GeneratorReactive {
872                        generator: generator_index,
873                        channel,
874                    },
875                    real,
876                );
877            }
878        }
879    }
880    for (source_index, source) in sources.iter().enumerate() {
881        for (channel, &node) in source.terminal_nodes.iter().enumerate() {
882            add_variable(
883                &mut balances[node].active,
884                LinDist3FlowVariable::SourceActive {
885                    source: source_index,
886                    channel,
887                },
888                1.0,
889            );
890            add_variable(
891                &mut balances[node].reactive,
892                LinDist3FlowVariable::SourceReactive {
893                    source: source_index,
894                    channel,
895                },
896                1.0,
897            );
898        }
899    }
900
901    Ok(LinDist3FlowDeviceData {
902        loads,
903        shunts,
904        generators,
905        sources,
906        balances,
907    })
908}
909
910/// Build the complete solver-neutral preparation for the supported L3F slice.
911///
912/// # Errors
913/// As [`build_lindist3flow_network_data`] and
914/// [`build_lindist3flow_device_data`].
915pub fn build_lindist3flow_preparation(
916    instance: &LinDist3FlowOpfInstance,
917) -> Result<LinDist3FlowPreparation> {
918    let network = build_lindist3flow_network_data(instance)?;
919    let devices = build_lindist3flow_device_data(instance, &network)?;
920    Ok(LinDist3FlowPreparation { network, devices })
921}
922
923#[cfg(test)]
924mod tests {
925    use std::f64::consts::PI;
926
927    use approx::assert_relative_eq;
928    use powerio_dist::{
929        Configuration, DistBus, DistGenerator, DistLine, DistLineCode, DistLoad, DistShunt,
930        MulticonductorNetwork, VoltageSource,
931    };
932    use powerio_prob::{
933        ConstraintSelection, LinDist3FlowBuildOptions, LinDist3FlowOpfInstance, McAcOpfInstance,
934        MulticonductorActiveConstraints,
935    };
936
937    use super::*;
938
939    fn one_phase_network() -> MulticonductorNetwork {
940        let terminal = vec!["1".to_owned()];
941        let mut network = MulticonductorNetwork::named("one-phase");
942        network
943            .buses_mut()
944            .push(DistBus::new("source", terminal.clone()));
945        network
946            .buses_mut()
947            .push(DistBus::new("load", terminal.clone()));
948        network
949            .line_codes_mut()
950            .push(DistLineCode::new("one", vec![vec![0.1]], vec![vec![0.2]]));
951        network.lines_mut().push(DistLine::new(
952            "line",
953            "source",
954            "load",
955            terminal.clone(),
956            terminal.clone(),
957            "one",
958            1.0,
959        ));
960        network.sources_mut().push(VoltageSource::new(
961            "grid",
962            "source",
963            terminal,
964            vec![230.0],
965            vec![0.0],
966        ));
967        network
968    }
969
970    fn three_phase_network() -> MulticonductorNetwork {
971        let terminals = vec!["1".to_owned(), "2".to_owned(), "3".to_owned()];
972        let mut network = MulticonductorNetwork::named("three-phase");
973        network
974            .buses_mut()
975            .push(DistBus::new("source", terminals.clone()));
976        network
977            .buses_mut()
978            .push(DistBus::new("load", terminals.clone()));
979        network.line_codes_mut().push(DistLineCode::new(
980            "three",
981            vec![
982                vec![0.1, 0.0, 0.0],
983                vec![0.0, 0.1, 0.0],
984                vec![0.0, 0.0, 0.1],
985            ],
986            vec![
987                vec![0.2, 0.0, 0.0],
988                vec![0.0, 0.2, 0.0],
989                vec![0.0, 0.0, 0.2],
990            ],
991        ));
992        network.lines_mut().push(DistLine::new(
993            "line",
994            "source",
995            "load",
996            terminals.clone(),
997            terminals.clone(),
998            "three",
999            1.0,
1000        ));
1001        network.sources_mut().push(VoltageSource::new(
1002            "grid",
1003            "source",
1004            terminals,
1005            vec![230.0; 3],
1006            vec![0.0, -2.0 * PI / 3.0, 2.0 * PI / 3.0],
1007        ));
1008        network
1009    }
1010
1011    fn evaluate(expression: &LinDist3FlowAffineExpression, squared_voltage: &[f64]) -> f64 {
1012        expression.constant
1013            + expression
1014                .voltage_terms
1015                .iter()
1016                .map(|term| term.coefficient * squared_voltage[term.node])
1017                .sum::<f64>()
1018    }
1019
1020    #[test]
1021    fn connection_incidence_matches_supported_physical_channels() {
1022        assert_eq!(
1023            lindist3flow_connection_incidence(Configuration::SinglePhase, 2, 1).unwrap(),
1024            [vec![1.0, -1.0]]
1025        );
1026        assert_eq!(
1027            lindist3flow_connection_incidence(Configuration::Wye, 2, 2).unwrap(),
1028            [vec![1.0, 0.0], vec![0.0, 1.0]]
1029        );
1030        assert_eq!(
1031            lindist3flow_connection_incidence(Configuration::Delta, 3, 3).unwrap()[2],
1032            [-1.0, 0.0, 1.0]
1033        );
1034        assert!(lindist3flow_connection_incidence(Configuration::Delta, 3, 1).is_err());
1035    }
1036
1037    #[test]
1038    fn constant_power_load_and_line_flows_have_the_reference_kcl_signs() {
1039        let mut network = one_phase_network();
1040        network.loads_mut().push(DistLoad::new(
1041            "demand",
1042            "load",
1043            vec!["1".to_owned()],
1044            Configuration::Wye,
1045            vec![1_000.0],
1046            vec![200.0],
1047        ));
1048        let instance =
1049            LinDist3FlowOpfInstance::from_network(network, LinDist3FlowBuildOptions::default())
1050                .unwrap();
1051        let preparation = build_lindist3flow_preparation(&instance).unwrap();
1052
1053        let source = &preparation.devices.balances[0];
1054        assert!(
1055            source
1056                .active
1057                .variable_terms
1058                .contains(&LinDist3FlowBalanceTerm {
1059                    variable: LinDist3FlowVariable::SourceActive {
1060                        source: 0,
1061                        channel: 0,
1062                    },
1063                    coefficient: 1.0,
1064                })
1065        );
1066        assert!(
1067            source
1068                .active
1069                .variable_terms
1070                .contains(&LinDist3FlowBalanceTerm {
1071                    variable: LinDist3FlowVariable::LineActive {
1072                        line: 0,
1073                        conductor: 0,
1074                    },
1075                    coefficient: -1.0,
1076                })
1077        );
1078        let load = &preparation.devices.balances[1];
1079        assert_relative_eq!(load.active.affine.constant, -1_000.0, epsilon = 1e-12);
1080        assert_relative_eq!(load.reactive.affine.constant, -200.0, epsilon = 1e-12);
1081        assert!(
1082            load.active
1083                .variable_terms
1084                .contains(&LinDist3FlowBalanceTerm {
1085                    variable: LinDist3FlowVariable::LineActive {
1086                        line: 0,
1087                        conductor: 0,
1088                    },
1089                    coefficient: 1.0,
1090                })
1091        );
1092    }
1093
1094    #[test]
1095    fn delta_zp_load_is_affine_and_exact_at_the_reference() {
1096        let mut network = three_phase_network();
1097        let mut load = DistLoad::new(
1098            "delta",
1099            "load",
1100            vec!["1".to_owned(), "2".to_owned(), "3".to_owned()],
1101            Configuration::Delta,
1102            vec![1_000.0, 2_000.0, 3_000.0],
1103            vec![100.0, 200.0, 300.0],
1104        );
1105        load.voltage_model = DistLoadVoltageModel::Zip {
1106            v_nom: vec![230.0 * 3.0_f64.sqrt()],
1107            alpha_z: vec![0.4],
1108            alpha_i: vec![0.0],
1109            alpha_p: vec![0.6],
1110            beta_z: vec![0.25],
1111            beta_i: vec![0.0],
1112            beta_p: vec![0.75],
1113        };
1114        network.loads_mut().push(load);
1115        let instance =
1116            LinDist3FlowOpfInstance::from_network(network, LinDist3FlowBuildOptions::default())
1117                .unwrap();
1118        let preparation = build_lindist3flow_preparation(&instance).unwrap();
1119        let squared_voltage = preparation
1120            .network
1121            .nodes
1122            .iter()
1123            .map(|node| node.reference_magnitude.powi(2))
1124            .collect::<Vec<_>>();
1125        let load = &preparation.devices.loads[0];
1126
1127        for channel in 0..3 {
1128            assert_relative_eq!(
1129                evaluate(&load.channel_power[channel].active, &squared_voltage),
1130                [1_000.0, 2_000.0, 3_000.0][channel],
1131                epsilon = 1e-9
1132            );
1133            let allocation_sum = (0..3)
1134                .map(|terminal| load.terminal_power_map.matrix[terminal][channel])
1135                .sum::<Complex64>();
1136            assert_relative_eq!(allocation_sum.re, 1.0, epsilon = 1e-12);
1137            assert_relative_eq!(allocation_sum.im, 0.0, epsilon = 1e-12);
1138        }
1139    }
1140
1141    #[test]
1142    fn generator_rows_keep_bounds_cost_soc_voltage_and_terminal_map() {
1143        let mut network = one_phase_network();
1144        let mut generator = DistGenerator::new(
1145            "der",
1146            "load",
1147            vec!["1".to_owned()],
1148            Configuration::Wye,
1149            vec![500.0],
1150            vec![0.0],
1151        );
1152        generator.p_min = Some(vec![0.0]);
1153        generator.p_max = Some(vec![1_000.0]);
1154        generator.q_min = Some(vec![-300.0]);
1155        generator.q_max = Some(vec![300.0]);
1156        generator.cost = Some(vec![0.2]);
1157        generator.s_max = Some(vec![1_100.0]);
1158        generator.i_max = Some(vec![5.0]);
1159        network.generators_mut().push(generator);
1160        let instance =
1161            LinDist3FlowOpfInstance::from_network(network, LinDist3FlowBuildOptions::default())
1162                .unwrap();
1163        let preparation = build_lindist3flow_preparation(&instance).unwrap();
1164        let channel = &preparation.devices.generators[0].channels[0];
1165
1166        assert_eq!(channel.active_min, Some(0.0));
1167        assert_eq!(channel.active_max, Some(1_000.0));
1168        assert_eq!(channel.reactive_min, Some(-300.0));
1169        assert_eq!(channel.reactive_max, Some(300.0));
1170        assert_eq!(channel.apparent_power_limit, Some(1_100.0));
1171        assert_eq!(channel.current_limit, Some(5.0));
1172        assert_relative_eq!(
1173            channel.active_objective_coefficient,
1174            0.0002,
1175            epsilon = 1e-15
1176        );
1177        assert_relative_eq!(channel.reference_winding_voltage, 230.0, epsilon = 1e-12);
1178        assert_relative_eq!(
1179            evaluate(
1180                &channel.squared_winding_voltage,
1181                &[230.0_f64.powi(2), 230.0_f64.powi(2)]
1182            ),
1183            230.0_f64.powi(2),
1184            epsilon = 1e-9
1185        );
1186    }
1187
1188    #[test]
1189    fn shunt_power_enters_balance_as_an_affine_absorption() {
1190        let mut network = one_phase_network();
1191        network.shunts_mut().push(DistShunt::new(
1192            "capacitive",
1193            "load",
1194            vec!["1".to_owned()],
1195            vec![vec![0.0]],
1196            vec![vec![0.01]],
1197        ));
1198        let instance =
1199            LinDist3FlowOpfInstance::from_network(network, LinDist3FlowBuildOptions::default())
1200                .unwrap();
1201        let preparation = build_lindist3flow_preparation(&instance).unwrap();
1202        let squared_voltage = vec![230.0_f64.powi(2); 2];
1203        let shunt = &preparation.devices.shunts[0].terminal_power[0];
1204
1205        assert_relative_eq!(
1206            evaluate(&shunt.active, &squared_voltage),
1207            0.0,
1208            epsilon = 1e-9
1209        );
1210        assert_relative_eq!(
1211            evaluate(&shunt.reactive, &squared_voltage),
1212            -0.01 * 230.0_f64.powi(2),
1213            epsilon = 1e-9
1214        );
1215        assert_relative_eq!(
1216            evaluate(
1217                &preparation.devices.balances[1].reactive.affine,
1218                &squared_voltage
1219            ),
1220            0.01 * 230.0_f64.powi(2),
1221            epsilon = 1e-9
1222        );
1223    }
1224
1225    #[test]
1226    fn inactive_constraint_families_remove_only_their_numerical_limits() {
1227        let mut network = one_phase_network();
1228        network.buses_mut()[1].v_min = Some(210.0);
1229        network.lines_mut()[0].i_max = Some(vec![100.0]);
1230        let mut generator = DistGenerator::new(
1231            "der",
1232            "load",
1233            vec!["1".to_owned()],
1234            Configuration::Wye,
1235            vec![500.0],
1236            vec![0.0],
1237        );
1238        generator.p_min = Some(vec![0.0]);
1239        generator.p_max = Some(vec![1_000.0]);
1240        network.generators_mut().push(generator);
1241        let mut constraints = MulticonductorActiveConstraints::default();
1242        constraints.terminal_voltage_bounds = ConstraintSelection::None;
1243        constraints.conductor_limits = ConstraintSelection::None;
1244        constraints.generator_capability = ConstraintSelection::None;
1245        let base = McAcOpfInstance::from_network(network)
1246            .unwrap()
1247            .with_constraints(constraints);
1248        let instance =
1249            LinDist3FlowOpfInstance::from_mc_ac(base, LinDist3FlowBuildOptions::default()).unwrap();
1250        let preparation = build_lindist3flow_preparation(&instance).unwrap();
1251
1252        assert!(preparation.network.nodes[1].squared_voltage_min.is_none());
1253        assert_eq!(preparation.network.lines[0].current_limit, [None]);
1254        assert!(
1255            preparation.devices.generators[0].channels[0]
1256                .active_min
1257                .is_none()
1258        );
1259    }
1260}
powerio_matrix/lindist3flow_devices.rs

powerio_matrix/
lindist3flow_devices.rs
