powerio_prob/solution/
socwr.rs

1//! Solution values for the PowerModels SOCWR relaxation of AC optimal power
2//! flow.
3
4use std::collections::BTreeMap;
5use std::sync::Arc;
6
7use powerio_core::Error;
8use powerio_tx::{BalancedNetwork, BusId};
9
10use crate::ThreeWindingTransformerTerminalPower;
11use crate::diagnostics::codes;
12use crate::instance::AcOpfInstance;
13use crate::operating::row_identity;
14use crate::solution::{Producer, Residuals, Termination};
15
16fn check_length(what: &str, actual: usize, expected: usize) -> Result<(), Error> {
17    if actual == expected {
18        Ok(())
19    } else {
20        Err(Error::new(
21            &codes::BUILD_SOLUTION_SHAPE_MISMATCH,
22            format!("{what} carries {actual} values; expected {expected}"),
23        ))
24    }
25}
26
27fn check_nonnegative(what: &str, values: &[f64]) -> Result<(), Error> {
28    if let Some((position, value)) = values
29        .iter()
30        .copied()
31        .enumerate()
32        .find(|(_, value)| !value.is_finite() || *value < 0.0)
33    {
34        return Err(Error::new(
35            &codes::BUILD_SOLUTION_MULTIPLIER_INVALID,
36            format!(
37                "{what} value {position} is {value}; multipliers must be finite and nonnegative"
38            ),
39        ));
40    }
41    Ok(())
42}
43
44/// The primal quantities reported by a SOCWR AC OPF solve.
45///
46/// Bus, branch, and generator vectors follow their source network table
47/// order. Power quantities use MW, MVAr, or MVA. The three voltage product
48/// quantities use per unit squared.
49#[derive(Clone, Debug, Default, PartialEq)]
50#[non_exhaustive]
51pub struct SocwrOpfValues {
52    /// `w[i] = |V_i|²`, by bus.
53    pub bus_voltage_magnitude_squared: Vec<f64>,
54    /// `wr[e] = Re(V_from * conj(V_to))`, by oriented branch.
55    pub branch_voltage_product_real: Vec<f64>,
56    /// `wi[e] = Im(V_from * conj(V_to))`, by oriented branch.
57    pub branch_voltage_product_imaginary: Vec<f64>,
58    /// Active generator dispatch, MW.
59    pub generator_active_power: Vec<f64>,
60    /// Reactive generator dispatch, MVAr.
61    pub generator_reactive_power: Vec<f64>,
62    /// Active power into each branch at its from terminal, MW.
63    pub branch_from_active_power: Vec<f64>,
64    /// Reactive power into each branch at its from terminal, MVAr.
65    pub branch_from_reactive_power: Vec<f64>,
66    /// Active power into each branch at its to terminal, MW.
67    pub branch_to_active_power: Vec<f64>,
68    /// Reactive power into each branch at its to terminal, MVAr.
69    pub branch_to_reactive_power: Vec<f64>,
70    /// Terminal powers for each three winding transformer, in transformer
71    /// table order.
72    pub three_winding_transformer_terminal_powers: Vec<ThreeWindingTransformerTerminalPower>,
73}
74
75/// Optional dual quantities reported by a SOCWR AC OPF solve.
76///
77/// `None` means the producer did not report that column. Present bus columns
78/// follow bus table order; present branch columns follow branch table order.
79#[derive(Clone, Debug, Default, PartialEq)]
80#[non_exhaustive]
81pub struct SocwrOpfDuals {
82    /// Objective derivative per added MW of active demand.
83    pub bus_active_power_marginal: Option<Vec<f64>>,
84    /// Objective derivative per added MVAr of reactive demand.
85    pub bus_reactive_power_marginal: Option<Vec<f64>>,
86    /// Multiplier on the from terminal apparent power limit.
87    pub branch_from_thermal_limit_multiplier: Option<Vec<f64>>,
88    /// Multiplier on the to terminal apparent power limit.
89    pub branch_to_thermal_limit_multiplier: Option<Vec<f64>>,
90}
91
92#[derive(Clone, Debug)]
93struct SolutionIndex {
94    buses: BTreeMap<BusId, usize>,
95    branches: BTreeMap<String, usize>,
96    generators: BTreeMap<String, usize>,
97}
98
99impl SolutionIndex {
100    fn new(network: &BalancedNetwork) -> Result<Self, Error> {
101        let mut buses = BTreeMap::new();
102        let mut branches = BTreeMap::new();
103        let mut generators = BTreeMap::new();
104        for (position, bus) in network.buses().iter().enumerate() {
105            if buses.insert(bus.id, position).is_some() {
106                return Err(duplicate_identity("bus", &bus.id.to_string()));
107            }
108        }
109        for (position, branch) in network.branches().iter().enumerate() {
110            let id = row_identity(branch.uid.as_deref(), "branches", position);
111            if branches.insert(id.clone(), position).is_some() {
112                return Err(duplicate_identity("branch", &id));
113            }
114        }
115        for (position, generator) in network.generators().iter().enumerate() {
116            let id = row_identity(generator.uid.as_deref(), "generators", position);
117            if generators.insert(id.clone(), position).is_some() {
118                return Err(duplicate_identity("generator", &id));
119            }
120        }
121        Ok(Self {
122            buses,
123            branches,
124            generators,
125        })
126    }
127}
128
129fn duplicate_identity(component: &str, id: &str) -> Error {
130    Error::new(
131        &codes::BUILD_OPERATING_POINT_IDENTITY_UNKNOWN,
132        format!("duplicate {component} identity `{id}`"),
133    )
134}
135
136/// A solution of the PowerModels SOCWR relaxation of AC optimal power flow.
137///
138/// This is a relaxation result and its objective is a lower bound. It is not
139/// an [`crate::AcOpfSolution`] and makes no claim that its voltage products
140/// recover an AC feasible voltage phasor.
141#[derive(Clone, Debug)]
142pub struct SocwrOpfSolution {
143    instance: Arc<AcOpfInstance>,
144    termination: Termination,
145    residuals: Residuals,
146    producer: Producer,
147    values: SocwrOpfValues,
148    duals: SocwrOpfDuals,
149    objective_lower_bound: f64,
150    index: SolutionIndex,
151}
152
153impl SocwrOpfSolution {
154    pub const FORMULATION: &'static str = "socwr";
155
156    /// Construct a SOCWR solution from table ordered primal values and its
157    /// objective lower bound.
158    ///
159    /// # Errors
160    /// Any value column has a length different from its network table.
161    pub fn new(
162        instance: Arc<AcOpfInstance>,
163        termination: Termination,
164        values: SocwrOpfValues,
165        objective_lower_bound: f64,
166    ) -> Result<Self, Error> {
167        let network = instance.network();
168        let buses = network.buses().len();
169        let branches = network.branches().len();
170        let generators = network.generators().len();
171        check_length(
172            "bus voltage magnitude squared",
173            values.bus_voltage_magnitude_squared.len(),
174            buses,
175        )?;
176        check_length(
177            "branch voltage product real part",
178            values.branch_voltage_product_real.len(),
179            branches,
180        )?;
181        check_length(
182            "branch voltage product imaginary part",
183            values.branch_voltage_product_imaginary.len(),
184            branches,
185        )?;
186        check_length(
187            "generator active power",
188            values.generator_active_power.len(),
189            generators,
190        )?;
191        check_length(
192            "generator reactive power",
193            values.generator_reactive_power.len(),
194            generators,
195        )?;
196        check_length(
197            "branch from terminal active power",
198            values.branch_from_active_power.len(),
199            branches,
200        )?;
201        check_length(
202            "branch from terminal reactive power",
203            values.branch_from_reactive_power.len(),
204            branches,
205        )?;
206        check_length(
207            "branch to terminal active power",
208            values.branch_to_active_power.len(),
209            branches,
210        )?;
211        check_length(
212            "branch to terminal reactive power",
213            values.branch_to_reactive_power.len(),
214            branches,
215        )?;
216        check_length(
217            "three winding transformer terminal powers",
218            values.three_winding_transformer_terminal_powers.len(),
219            network.transformers_3w().len(),
220        )?;
221        let index = SolutionIndex::new(network)?;
222        Ok(Self {
223            instance,
224            termination,
225            residuals: Residuals::default(),
226            producer: None,
227            values,
228            duals: SocwrOpfDuals::default(),
229            objective_lower_bound,
230            index,
231        })
232    }
233
234    #[must_use]
235    pub const fn formulation(&self) -> &'static str {
236        Self::FORMULATION
237    }
238
239    #[must_use]
240    pub fn instance(&self) -> &AcOpfInstance {
241        &self.instance
242    }
243
244    #[must_use]
245    pub fn shared_instance(&self) -> Arc<AcOpfInstance> {
246        Arc::clone(&self.instance)
247    }
248
249    #[must_use]
250    pub fn network(&self) -> &BalancedNetwork {
251        self.instance.network()
252    }
253
254    #[must_use]
255    pub const fn termination(&self) -> &Termination {
256        &self.termination
257    }
258
259    #[must_use]
260    pub const fn residuals(&self) -> &Residuals {
261        &self.residuals
262    }
263
264    #[must_use]
265    pub fn producer(&self) -> Option<&str> {
266        self.producer.as_deref()
267    }
268
269    #[must_use]
270    pub const fn values(&self) -> &SocwrOpfValues {
271        &self.values
272    }
273
274    #[must_use]
275    pub const fn duals(&self) -> &SocwrOpfDuals {
276        &self.duals
277    }
278
279    #[must_use]
280    pub const fn objective_lower_bound(&self) -> f64 {
281        self.objective_lower_bound
282    }
283
284    #[must_use]
285    pub fn with_producer(mut self, producer: impl Into<String>) -> Self {
286        self.producer = Some(producer.into());
287        self
288    }
289
290    #[must_use]
291    pub const fn with_residuals(mut self, residuals: Residuals) -> Self {
292        self.residuals = residuals;
293        self
294    }
295
296    /// Attach optional dual columns.
297    ///
298    /// # Errors
299    /// A present column has the wrong length, or a thermal limit multiplier
300    /// is negative or nonfinite.
301    pub fn with_duals(mut self, duals: SocwrOpfDuals) -> Result<Self, Error> {
302        let buses = self.network().buses().len();
303        let branches = self.network().branches().len();
304        if let Some(values) = &duals.bus_active_power_marginal {
305            check_length("bus active power marginal", values.len(), buses)?;
306        }
307        if let Some(values) = &duals.bus_reactive_power_marginal {
308            check_length("bus reactive power marginal", values.len(), buses)?;
309        }
310        if let Some(values) = &duals.branch_from_thermal_limit_multiplier {
311            check_length(
312                "branch from terminal thermal limit multiplier",
313                values.len(),
314                branches,
315            )?;
316            check_nonnegative("branch from terminal thermal limit multiplier", values)?;
317        }
318        if let Some(values) = &duals.branch_to_thermal_limit_multiplier {
319            check_length(
320                "branch to terminal thermal limit multiplier",
321                values.len(),
322                branches,
323            )?;
324            check_nonnegative("branch to terminal thermal limit multiplier", values)?;
325        }
326        self.duals = duals;
327        Ok(self)
328    }
329
330    pub fn bus_order(&self) -> impl ExactSizeIterator<Item = BusId> + '_ {
331        self.network().buses().iter().map(|bus| bus.id)
332    }
333
334    pub fn branch_order(&self) -> impl ExactSizeIterator<Item = String> + '_ {
335        self.network()
336            .branches()
337            .iter()
338            .enumerate()
339            .map(|(position, branch)| row_identity(branch.uid.as_deref(), "branches", position))
340    }
341
342    pub fn generator_order(&self) -> impl ExactSizeIterator<Item = String> + '_ {
343        self.network()
344            .generators()
345            .iter()
346            .enumerate()
347            .map(|(position, generator)| {
348                row_identity(generator.uid.as_deref(), "generators", position)
349            })
350    }
351
352    #[must_use]
353    pub fn bus_voltage_magnitude_squared(&self, bus: BusId) -> Option<f64> {
354        Some(self.values.bus_voltage_magnitude_squared[*self.index.buses.get(&bus)?])
355    }
356
357    #[must_use]
358    pub fn branch_voltage_product_real(&self, branch: &str) -> Option<f64> {
359        Some(self.values.branch_voltage_product_real[*self.index.branches.get(branch)?])
360    }
361
362    #[must_use]
363    pub fn branch_voltage_product_imaginary(&self, branch: &str) -> Option<f64> {
364        Some(self.values.branch_voltage_product_imaginary[*self.index.branches.get(branch)?])
365    }
366
367    #[must_use]
368    pub fn generator_active_power(&self, generator: &str) -> Option<f64> {
369        Some(self.values.generator_active_power[*self.index.generators.get(generator)?])
370    }
371
372    #[must_use]
373    pub fn generator_reactive_power(&self, generator: &str) -> Option<f64> {
374        Some(self.values.generator_reactive_power[*self.index.generators.get(generator)?])
375    }
376
377    #[must_use]
378    pub fn branch_from_active_power(&self, branch: &str) -> Option<f64> {
379        Some(self.values.branch_from_active_power[*self.index.branches.get(branch)?])
380    }
381
382    #[must_use]
383    pub fn branch_from_reactive_power(&self, branch: &str) -> Option<f64> {
384        Some(self.values.branch_from_reactive_power[*self.index.branches.get(branch)?])
385    }
386
387    #[must_use]
388    pub fn branch_to_active_power(&self, branch: &str) -> Option<f64> {
389        Some(self.values.branch_to_active_power[*self.index.branches.get(branch)?])
390    }
391
392    #[must_use]
393    pub fn branch_to_reactive_power(&self, branch: &str) -> Option<f64> {
394        Some(self.values.branch_to_reactive_power[*self.index.branches.get(branch)?])
395    }
396}
397
398#[cfg(test)]
399mod tests {
400    use super::*;
401    use powerio_core::Source;
402
403    fn instance() -> Arc<AcOpfInstance> {
404        let path = concat!(env!("CARGO_MANIFEST_DIR"), "/../tests/data/case9.m");
405        let network = powerio_tx::parse(Source::open(path).unwrap())
406            .unwrap()
407            .into_value();
408        Arc::new(AcOpfInstance::from_network(network).unwrap())
409    }
410
411    fn values(instance: &AcOpfInstance) -> SocwrOpfValues {
412        let buses = instance.network().buses().len();
413        let branches = instance.network().branches().len();
414        let generators = instance.network().generators().len();
415        SocwrOpfValues {
416            bus_voltage_magnitude_squared: vec![1.0; buses],
417            branch_voltage_product_real: vec![0.99; branches],
418            branch_voltage_product_imaginary: vec![0.01; branches],
419            generator_active_power: vec![10.0; generators],
420            generator_reactive_power: vec![2.0; generators],
421            branch_from_active_power: vec![3.0; branches],
422            branch_from_reactive_power: vec![0.5; branches],
423            branch_to_active_power: vec![-2.9; branches],
424            branch_to_reactive_power: vec![-0.4; branches],
425            three_winding_transformer_terminal_powers: Vec::new(),
426        }
427    }
428
429    #[test]
430    fn result_is_explicitly_a_relaxation_lower_bound() {
431        let instance = instance();
432        let solution = SocwrOpfSolution::new(
433            Arc::clone(&instance),
434            Termination::Converged,
435            values(&instance),
436            5_000.0,
437        )
438        .unwrap()
439        .with_producer("test-solver");
440        assert_eq!(solution.formulation(), "socwr");
441        assert!((solution.objective_lower_bound() - 5_000.0).abs() < f64::EPSILON);
442        assert_eq!(solution.producer(), Some("test-solver"));
443        assert!(
444            (solution.bus_voltage_magnitude_squared(BusId(1)).unwrap() - 1.0).abs() < f64::EPSILON
445        );
446        let branch = solution.branch_order().next().unwrap();
447        assert!(
448            (solution.branch_voltage_product_real(&branch).unwrap() - 0.99).abs() < f64::EPSILON
449        );
450        assert!(std::ptr::eq(solution.instance(), instance.as_ref()));
451    }
452
453    #[test]
454    fn dimensions_and_multiplier_signs_are_checked() {
455        let instance = instance();
456        let mut wrong = values(&instance);
457        wrong.bus_voltage_magnitude_squared.push(1.0);
458        assert!(
459            SocwrOpfSolution::new(Arc::clone(&instance), Termination::Converged, wrong, 0.0,)
460                .is_err()
461        );
462
463        let solution = SocwrOpfSolution::new(
464            Arc::clone(&instance),
465            Termination::Converged,
466            values(&instance),
467            0.0,
468        )
469        .unwrap();
470        let branches = instance.network().branches().len();
471        assert!(
472            solution
473                .with_duals(SocwrOpfDuals {
474                    branch_from_thermal_limit_multiplier: Some(vec![-1.0; branches]),
475                    ..SocwrOpfDuals::default()
476                })
477                .is_err()
478        );
479    }
480}
powerio_prob/solution/socwr.rs

powerio_prob/solution/
socwr.rs
