pub fn run_power_flow(
buses: Vec<Bus>,
lines: &[Line],
transformers: &[Transformer],
shunts: &[ShuntAdm],
taps: TapData<'_>,
opts: PowerFlowOptions,
) -> PowerFlowReportExpand description
Runs a power flow by the method named in opts.
Takes branch lists rather than a finished Y-bus because
PowerFlowMethod::Dc is formulated per branch — it needs each branch’s
own x and tap, which a Y-bus has already summed away — and never
builds one at all. The other two methods assemble their Y-bus here, from
these same lists plus shunts.
run_power_flow_analysis and run_power_flow_analysis_from_ybus are
unchanged and remain the shortest path to an ordinary Newton solve;
PowerFlowOptions::default() here reproduces exactly what they do.
Examples found in repository?
examples/iidm_area_probe.rs (lines 17-28)
2fn main() {
3 for path in std::env::args().skip(1) {
4 let net = match gridoxide::iidm::read(&path) {
5 Ok(n) => n,
6 Err(e) => {
7 println!("{path}: {e}");
8 continue;
9 }
10 };
11 let a = &net.areas;
12 println!("=== {path}");
13 println!(" {:?}", a.report);
14 if a.ids.is_empty() {
15 continue;
16 }
17 let report = gridoxide::run_power_flow(
18 net.buses.clone(),
19 &net.lines,
20 &net.transformers,
21 &net.shunts,
22 gridoxide::TapData::none(),
23 gridoxide::PowerFlowOptions {
24 tol: 1e-10,
25 max_iter: 40,
26 ..Default::default()
27 },
28 );
29 let x = gridoxide::outerloop::AreaInterchange::measure(
30 &report.buses,
31 &net.lines,
32 &net.transformers,
33 &a.of_bus,
34 a.ids.len(),
35 );
36 println!(" solve {:?}", report.stats.status);
37 for (i, (id, name)) in a.ids.iter().enumerate() {
38 let n = a.of_bus.iter().filter(|z| **z == Some(i)).count();
39 println!(
40 " {id:<6} {name:<28} {n:>3} bus(es) scheduled export {:>9.3} MW actual {:>9.3} MW",
41 a.targets[i] * net.base_mva,
42 x[i] * net.base_mva
43 );
44 }
45 }
46}More examples
examples/ucte_area_probe.rs (lines 7-18)
2fn main() {
3 for path in std::env::args().skip(1) {
4 let net = gridoxide::ucte::read(&path).expect("import");
5 let (of_bus, names) = net.country_areas();
6 let (lines, transformers) = net.as_switched();
7 let report = gridoxide::run_power_flow(
8 net.buses.clone(),
9 &lines,
10 &transformers,
11 &net.shunts,
12 gridoxide::TapData::none(),
13 gridoxide::PowerFlowOptions {
14 tol: 1e-10,
15 max_iter: 40,
16 ..Default::default()
17 },
18 );
19 let x = gridoxide::outerloop::AreaInterchange::measure(
20 &report.buses,
21 &lines,
22 &transformers,
23 &of_bus,
24 names.len(),
25 );
26 println!(
27 "{path}: {} areas {names:?}, solve {:?}",
28 names.len(),
29 report.stats.status
30 );
31 for (i, n) in names.iter().enumerate() {
32 let buses = of_bus.iter().filter(|a| **a == Some(i)).count();
33 println!(
34 " {n}: {buses:>3} bus(es), exporting {:>9.3} MW",
35 x[i] * net.base_mva
36 );
37 }
38 println!(
39 " unassigned buses: {}",
40 of_bus.iter().filter(|a| a.is_none()).count()
41 );
42 println!(
43 " sum of positions (= tie losses): {:.3} MW",
44 x.iter().sum::<f64>() * net.base_mva
45 );
46 }
47}examples/tap_control_probe.rs (lines 27-44)
6fn main() {
7 let root =
8 Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/CGMES-Test-Configurations/v3.0");
9 for name in std::env::args().skip(1) {
10 let dir = root.join(&name);
11 let mut paths: Vec<PathBuf> = std::fs::read_dir(&dir)
12 .unwrap_or_else(|e| panic!("{}: {e}", dir.display()))
13 .filter_map(|e| e.ok().map(|e| e.path()))
14 .filter(|p| p.extension().is_some_and(|e| e == "xml"))
15 .collect();
16 paths.sort();
17 let refs: Vec<&Path> = paths.iter().map(|p| p.as_path()).collect();
18 let ds = gridoxide::cgmes::load_profiles(&refs).expect("load");
19 let net = gridoxide::cgmes::cgmes_to_network(&ds, 100e6).expect("convert");
20
21 let before: Vec<Option<i32>> = net
22 .tap_changers
23 .iter()
24 .map(|c| c.as_ref().map(|c| c.position))
25 .collect();
26
27 let report = gridoxide::run_power_flow(
28 net.buses.clone(),
29 &net.lines,
30 &net.transformers,
31 &net.shunts,
32 gridoxide::TapData {
33 changers: &net.tap_changers,
34 regulation: &net.regulation,
35 },
36 gridoxide::PowerFlowOptions {
37 control_taps: true,
38 enforce_q_limits: true,
39 max_outer_iter: 40,
40 tol: 1e-8,
41 max_iter: 30,
42 ..Default::default()
43 },
44 );
45 let outer = report.outer.as_ref().expect("loops ran");
46 println!("=== {name}: {:?}", report.stats.status);
47 println!(" {:?}", outer.report.as_ref().unwrap());
48 for c in &outer.taps {
49 println!(
50 " tf {:>3} {} -> {} steps {} reversals {} {:?}",
51 c.transformer,
52 c.initial_position,
53 c.final_position,
54 c.steps_moved,
55 c.direction_changes,
56 c.outcome
57 );
58 }
59 let moved: Vec<usize> = outer
60 .changers
61 .iter()
62 .enumerate()
63 .filter(|(i, c)| c.as_ref().map(|c| c.position) != before[*i])
64 .map(|(i, _)| i)
65 .collect();
66 println!(" transformers whose tap moved: {moved:?}");
67 }
68}examples/cgmes_area_probe.rs (lines 25-36)
4fn main() {
5 let root =
6 Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/CGMES-Test-Configurations/v3.0");
7 for name in std::env::args().skip(1) {
8 let dir = root.join(&name);
9 if !dir.is_dir() {
10 println!("{name}: absent");
11 continue;
12 }
13 let mut paths: Vec<PathBuf> = std::fs::read_dir(&dir)
14 .unwrap()
15 .filter_map(|e| e.ok().map(|e| e.path()))
16 .filter(|p| p.extension().is_some_and(|x| x == "xml"))
17 .collect();
18 paths.sort();
19 let refs: Vec<&Path> = paths.iter().map(|p| p.as_path()).collect();
20 let ds = gridoxide::cgmes::load_profiles(&refs).expect("load");
21 let net = gridoxide::cgmes::cgmes_to_network(&ds, 100e6).expect("convert");
22 let areas = gridoxide::cgmes::cgmes_control_areas(&ds, &net, 100e6).expect("areas");
23 println!("=== {name}");
24 println!(" {:?}", areas.report);
25 let report = gridoxide::run_power_flow(
26 net.buses.clone(),
27 &net.lines,
28 &net.transformers,
29 &net.shunts,
30 gridoxide::TapData::none(),
31 gridoxide::PowerFlowOptions {
32 tol: 1e-8,
33 max_iter: 40,
34 ..Default::default()
35 },
36 );
37 let x = gridoxide::outerloop::AreaInterchange::measure(
38 &report.buses,
39 &net.lines,
40 &net.transformers,
41 &areas.of_bus,
42 areas.ids.len(),
43 );
44 println!(" solve {:?}", report.stats.status);
45 // What sits at the buses no area claims — the boundary X-nodes.
46 let unassigned: Vec<usize> = areas
47 .of_bus
48 .iter()
49 .enumerate()
50 .filter(|(_, a)| a.is_none())
51 .map(|(i, _)| i)
52 .collect();
53 let inj: f64 = unassigned.iter().map(|&i| report.buses[i].p_spec).sum();
54 println!(
55 " {} bus(es) in no area, injecting {:.3} MW between them",
56 unassigned.len(),
57 inj * 100.0
58 );
59 println!(
60 " sum of area positions: {:.3} MW",
61 x.iter().sum::<f64>() * 100.0
62 );
63 for (i, (mrid, nm)) in areas.ids.iter().enumerate() {
64 let buses = areas.of_bus.iter().filter(|a| **a == Some(i)).count();
65 println!(
66 " {:<4} ({}) {:>4} bus(es) scheduled export {:>10.3} MW actual {:>10.3} MW",
67 nm,
68 &mrid[..8.min(mrid.len())],
69 buses,
70 areas.targets[i] * 100.0,
71 x[i] * 100.0
72 );
73 }
74 }
75}examples/tap_deadband_probe.rs (lines 47-58)
9fn main() {
10 let root =
11 Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/CGMES-Test-Configurations/v3.0");
12 for name in std::env::args().skip(1) {
13 let dir = root.join(&name);
14 let mut paths: Vec<PathBuf> = std::fs::read_dir(&dir)
15 .unwrap()
16 .filter_map(|e| e.ok().map(|e| e.path()))
17 .filter(|p| p.extension().is_some_and(|e| e == "xml"))
18 .collect();
19 paths.sort();
20 let refs: Vec<&Path> = paths.iter().map(|p| p.as_path()).collect();
21 let ds = gridoxide::cgmes::load_profiles(&refs).expect("load");
22 let mut net = gridoxide::cgmes::cgmes_to_network(&ds, 100e6).expect("convert");
23 // FullGrid is an HVDC model; without resolving its converters into
24 // AC-side injections the AC solve has nothing driving those buses.
25 match gridoxide::cgmes::cgmes_resolve_dc_converters(&ds, &mut net.buses, 100e6) {
26 Ok(_) => println!(" dc: resolved"),
27 Err(e) => println!(" dc: {e}"),
28 }
29 let index = gridoxide::cgmes::cgmes_topological_node_bus_index(&ds).expect("index");
30
31 // The published SvVoltage per bus, in per unit.
32 let mut published = vec![f64::NAN; net.buses.len()];
33 for mrid in &ds.by_type["SvVoltage"] {
34 let sv = &ds.entries[mrid];
35 let (Some(tn), Some(v)) = (
36 sv.get_ref("SvVoltage.TopologicalNode"),
37 sv.get_f64("SvVoltage.v"),
38 ) else {
39 continue;
40 };
41 if let Some(&b) = index.get(tn) {
42 published[b] = v * 1e3 / net.buses[b].u_rated;
43 }
44 }
45
46 // gridoxide's own solve at the positions the document names.
47 let report = gridoxide::run_power_flow(
48 net.buses.clone(),
49 &net.lines,
50 &net.transformers,
51 &net.shunts,
52 gridoxide::TapData::none(),
53 gridoxide::PowerFlowOptions {
54 tol: 1e-8,
55 max_iter: 30,
56 ..Default::default()
57 },
58 );
59 // And the same, after letting the tap loops run.
60 let after = gridoxide::run_power_flow(
61 net.buses.clone(),
62 &net.lines,
63 &net.transformers,
64 &net.shunts,
65 gridoxide::TapData {
66 changers: &net.tap_changers,
67 regulation: &net.regulation,
68 },
69 gridoxide::PowerFlowOptions {
70 control_taps: true,
71 max_outer_iter: 60,
72 tol: 1e-8,
73 max_iter: 30,
74 ..Default::default()
75 },
76 );
77 println!(
78 "=== {name}: base solve {:?}, controlled solve {:?}",
79 report.stats.status, after.stats.status
80 );
81 println!(
82 "{:>4} {:>9} {:>9} {:>9} {:>9} {:>9} {:>9}",
83 "bus", "target", "half-db", "gridoxide", "|err|", "published", "|err|"
84 );
85 for r in &net.regulation {
86 if r.mode != gridoxide::outerloop::RegulationMode::Voltage {
87 continue;
88 }
89 let mine = report.buses[r.controlled_bus].voltage_mag;
90 let theirs = published[r.controlled_bus];
91 println!(
92 "{:>4} {:>9.5} {:>9.5} {:>9.5} {:>9.5}{} {:>9.5} {:>9.5}{}",
93 r.controlled_bus,
94 r.target,
95 r.deadband / 2.0,
96 mine,
97 (mine - r.target).abs(),
98 if (mine - r.target).abs() > r.deadband / 2.0 {
99 " OUT"
100 } else {
101 " in"
102 },
103 theirs,
104 (theirs - r.target).abs(),
105 if (theirs - r.target).abs() > r.deadband / 2.0 {
106 " OUT"
107 } else {
108 " in"
109 },
110 );
111 let controlled = after.buses[r.controlled_bus].voltage_mag;
112 println!(
113 " -> after tap control: {controlled:.5} |err| {:.5}{}",
114 (controlled - r.target).abs(),
115 if (controlled - r.target).abs() > r.deadband / 2.0 {
116 " OUT"
117 } else {
118 " in"
119 },
120 );
121 }
122 }
123}examples/dispatch_probe.rs (lines 26-38)
6fn main() {
7 let path = std::env::args().nth(1).expect("path to a CGMES directory");
8 let dir = std::path::Path::new(&path);
9 let mut paths: Vec<std::path::PathBuf> = std::fs::read_dir(dir)
10 .expect("dir")
11 .filter_map(|e| e.ok().map(|e| e.path()))
12 .filter(|p| p.extension().is_some_and(|x| x == "xml"))
13 .collect();
14 paths.sort();
15 let refs: Vec<&std::path::Path> = paths.iter().map(|p| p.as_path()).collect();
16 let ds = gridoxide::cgmes::load_profiles(&refs).expect("load");
17 let net = gridoxide::cgmes::cgmes_to_network(&ds, 100e6).expect("convert");
18 let vc = &net.voltage_control;
19 println!(
20 "{} regulated bus(es), {} shared, {} regulating machine(s)",
21 vc.regulated_buses,
22 vc.shared_buses,
23 vc.machines.len()
24 );
25
26 let report = gridoxide::run_power_flow(
27 net.buses.clone(),
28 &net.lines,
29 &net.transformers,
30 &net.shunts,
31 gridoxide::TapData::none(),
32 PowerFlowOptions {
33 enforce_q_limits: true,
34 max_iter: 60,
35 max_outer_iter: 60,
36 ..Default::default()
37 },
38 );
39 println!("solve: {:?}", report.stats.status);
40
41 let mut ybus = build_ybus(net.buses.len(), &net.lines, &net.transformers);
42 stamp_shunts(&mut ybus, &net.shunts);
43 let ybus = ybus.finish();
44 let d = allocate(&vc.machines, &vc.nonregulating_q, &report.buses, &ybus);
45
46 // Group by controlled bus, show the shared ones.
47 let mut by_bus: std::collections::BTreeMap<usize, Vec<&gridoxide::dispatch::MachineDispatch>> =
48 Default::default();
49 for bd in &d {
50 for m in &bd.machines {
51 by_bus.entry(m.controls_bus).or_default().push(m);
52 }
53 }
54 let shared: Vec<_> = by_bus.iter().filter(|(_, v)| v.len() > 1).collect();
55 println!("{} bus(es) with >1 machine\n", shared.len());
56
57 for (bus, ms) in shared.iter().take(4) {
58 let refs: Vec<&gridoxide::types::RegulatingMachine> = vc
59 .machines
60 .iter()
61 .filter(|m| m.controls_bus == **bus)
62 .collect();
63 let (_, basis) = reactive_keys(&refs);
64 let total: f64 = ms.iter().map(|m| m.q).sum();
65 println!(
66 "bus {bus}: {} machines, keys from {:?}, total Q = {:+.6} pu",
67 ms.len(),
68 basis,
69 total
70 );
71 for m in ms.iter() {
72 println!(
73 " {:<22} q={:+.6} share={:.3} range=[{:+.4},{:+.4}]{}",
74 &m.id[..m.id.len().min(22)],
75 m.q,
76 m.share,
77 m.q_min,
78 m.q_max,
79 if m.at_limit { " AT LIMIT" } else { "" }
80 );
81 }
82 println!();
83 }
84
85 // Invariants worth seeing before writing a test around them.
86 // Compare against the fixture's own published per-terminal solution.
87 let mut published: std::collections::HashMap<String, f64> = Default::default();
88 if let Some(mrids) = ds.by_type.get("Terminal") {
89 let mut eq_of: std::collections::HashMap<&str, &str> = Default::default();
90 for m in mrids {
91 if let Some(e) = ds.entries[m].get_ref("Terminal.ConductingEquipment") {
92 eq_of.insert(m.as_str(), e);
93 }
94 }
95 if let Some(fl) = ds.by_type.get("SvPowerFlow") {
96 for m in fl {
97 let sv = &ds.entries[m];
98 let (Some(t), Some(q)) = (
99 sv.get_ref("SvPowerFlow.Terminal"),
100 sv.get_f64("SvPowerFlow.q"),
101 ) else {
102 continue;
103 };
104 if let Some(e) = eq_of.get(t) {
105 published.insert((*e).to_string(), -q / 100.0);
106 }
107 }
108 }
109 }
110 let mut matched = 0usize;
111 let mut worst_m = 0.0f64;
112 let mut sum_abs = 0.0f64;
113 let mut sum_err = 0.0f64;
114 for bd in &d {
115 for m in &bd.machines {
116 if let Some(&exp) = published.get(&m.id) {
117 matched += 1;
118 worst_m = worst_m.max((m.q - exp).abs());
119 sum_abs += exp.abs();
120 sum_err += (m.q - exp).abs();
121 }
122 }
123 }
124 println!(
125 "published per-machine Q: matched {matched}/{} worst |diff| = {worst_m:.4} pu \
126 mean |diff| / mean |published| = {:.1}%",
127 d.iter().map(|b| b.machines.len()).sum::<usize>(),
128 if sum_abs > 0.0 {
129 100.0 * sum_err / sum_abs
130 } else {
131 0.0
132 }
133 );
134
135 // Decompose: does the *bus total* agree (the solve) or the *split* (the rule)?
136 let (mut bus_err, mut bus_abs) = (0.0f64, 0.0f64);
137 let (mut split_err, mut split_abs) = (0.0f64, 0.0f64);
138 let (mut shared_split_err, mut shared_split_abs) = (0.0f64, 0.0f64);
139 let mut worst_bus = (0.0f64, 0usize);
140 for bd in &d {
141 let pub_sum: f64 = bd
142 .machines
143 .iter()
144 .filter_map(|m| published.get(&m.id))
145 .sum();
146 let n_pub = bd
147 .machines
148 .iter()
149 .filter(|m| published.contains_key(&m.id))
150 .count();
151 if n_pub != bd.machines.len() {
152 continue;
153 }
154 bus_err += (bd.attributed - pub_sum).abs();
155 bus_abs += pub_sum.abs();
156 if (bd.attributed - pub_sum).abs() > worst_bus.0 {
157 worst_bus = ((bd.attributed - pub_sum).abs(), bd.bus);
158 }
159 // Split quality, measured against the *published* total so the solve
160 // difference is factored out.
161 let (keys, _) = reactive_keys(
162 &vc.machines
163 .iter()
164 .filter(|m| m.controls_bus == bd.bus)
165 .collect::<Vec<_>>(),
166 );
167 for (k, m) in bd.machines.iter().enumerate() {
168 let exp = published[&m.id];
169 let ours_on_pub_total = pub_sum * keys[k];
170 split_err += (ours_on_pub_total - exp).abs();
171 split_abs += exp.abs();
172 if bd.machines.len() > 1 {
173 shared_split_err += (ours_on_pub_total - exp).abs();
174 shared_split_abs += exp.abs();
175 }
176 }
177 }
178 println!(
179 " bus totals (the solve): mean rel err {:.1}% worst {:.4} pu at bus {}",
180 100.0 * bus_err / bus_abs.max(1e-12),
181 worst_bus.0,
182 worst_bus.1
183 );
184 println!(
185 " the split (the rule) : mean rel err {:.1}% (shared buses only: {:.1}%)",
186 100.0 * split_err / split_abs.max(1e-12),
187 100.0 * shared_split_err / shared_split_abs.max(1e-12)
188 );
189
190 // Basis distribution over shared buses, and unattributed at a physical threshold.
191 let mut basis_counts: std::collections::BTreeMap<String, usize> = Default::default();
192 for bd in &d {
193 if bd.machines.len() > 1 {
194 *basis_counts.entry(format!("{:?}", bd.basis)).or_default() += 1;
195 }
196 }
197 println!("shared-bus key basis: {basis_counts:?}");
198 for thr in [1e-9f64, 1e-7, 1e-6, 1e-5] {
199 let n = d.iter().filter(|b| b.unattributed.abs() > thr).count();
200 let t: f64 = d
201 .iter()
202 .filter(|b| b.unattributed.abs() > thr)
203 .map(|b| b.unattributed.abs())
204 .sum();
205 println!(" unattributed > {thr:.0e}: {n} bus(es), {t:.4} pu");
206 }
207 // Which shared buses fall back, and why?
208 for bd in d
209 .iter()
210 .filter(|b| b.machines.len() > 1 && format!("{:?}", b.basis) != "Capability")
211 {
212 println!(" FALLBACK bus {} -> {:?}", bd.bus, bd.basis);
213 for m in &bd.machines {
214 println!(
215 " {} range=[{:+.6},{:+.6}] width={:.6}",
216 &m.id[..12],
217 m.q_min,
218 m.q_max,
219 m.q_max - m.q_min
220 );
221 }
222 }
223
224 let (_, qall) = gridoxide::network::power_injections(&report.buses, &ybus);
225 let mut worst: Vec<(f64, usize, usize, usize, f64, f64)> = Vec::new();
226 let mut violations = 0;
227 for (bus, ms) in &by_bus {
228 let total: f64 = ms.iter().map(|m| m.q).sum();
229 let expected = qall[*bus] - vc.nonregulating_q[*bus];
230 let pinned = ms.iter().filter(|m| m.at_limit).count();
231 worst.push((
232 (total - expected).abs(),
233 *bus,
234 ms.len(),
235 pinned,
236 total,
237 expected,
238 ));
239 for m in ms.iter() {
240 if m.q > m.q_max + 1e-9 || m.q < m.q_min - 1e-9 {
241 violations += 1;
242 }
243 }
244 }
245 let n_bad = worst.iter().filter(|w| w.0 > 1e-9).count();
246 let total_bad: f64 = worst.iter().map(|w| w.0).sum();
247 // How many regulated buses have reactive load co-located with the machine?
248 let with_load = by_bus
249 .keys()
250 .filter(|b| vc.nonregulating_q[**b].abs() > 1e-9)
251 .count();
252 println!(
253 "regulated buses: {} with co-located reactive injection: {with_load} \
254 unattributable: {n_bad} total unattributed = {total_bad:.4} pu",
255 by_bus.len()
256 );
257 worst.sort_by(|a, b| b.0.total_cmp(&a.0));
258 println!("limit violations = {violations}");
259 println!("worst mismatches (err, bus, machines, pinned, sum, expected):");
260 for w in worst.iter().take(6) {
261 let bt = report.buses[w.1].bus_type;
262 println!(
263 " {:.3e} bus {:4} n={} pinned={} sum={:+.6} expected={:+.6} type={:?} qlim=[{:+.4},{:+.4}]",
264 w.0, w.1, w.2, w.3, w.4, w.5, bt, report.buses[w.1].q_min, report.buses[w.1].q_max
265 );
266 }
267}