# Benchmark scripts
Timing and accuracy harnesses comparing gridoxide against other power-system solvers:
| 1–3 | Synthetic radial MV/LV grids | power-grid-model |
| 4 | 12 MATPOWER test cases | PGM, lightsim2grid, pypowsybl, pandapower, VeraGrid |
| 4b | Batched power flow (multi-core) | itself, thread scaling |
| 4c–4d | Block-diagonal embedding | JAX oracle, independent solves |
| 5–6 | CGMES conformance configurations | pypowsybl |
| 7 | State estimation | power-grid-model |
## Python bindings
Every script here runs purely in Python via gridoxide's native bindings (`src/python.rs`, built with
[maturin](https://www.maturin.rs/)) — no subprocess into a Rust binary, no stdout parsing.
```bash
pip install maturin
maturin develop --release --features python # scalar + block backends
maturin develop --release --features python,klu # + klu (needs a C compiler, libclang)
maturin develop --release --features python,pardiso # + pardiso (needs MKLROOT set)
```
`maturin develop` needs an active virtualenv (or `VIRTUAL_ENV=/path/to/venv`). The `python` feature is
`#[cfg]`-gated, so plain `cargo build`/`cargo test` is unaffected — but never combine `--features python`
with `cargo build --examples`/`cargo test` in one invocation: PyO3's `extension-module` makes standalone
binaries fail to link. That's an expected PyO3 constraint, not a gridoxide issue.
```python
import gridoxide
model = gridoxide.PowerFlowModel.from_pgm_json("grid.json", backend="klu")
model.solve()
print(model.voltage_mag(), model.voltage_ang())
```
`PowerFlowModel` wraps `solver::PersistentSolver`, so repeated `.solve()` calls reuse the cached symbolic
factorization — matching how every other tool here is driven (construct once, time repeated solve calls with
`time.perf_counter()`). All comparisons below are warm-vs-warm.
## 1. Generate a synthetic benchmark grid
```bash
python3 scripts/bench/generate_grid.py grid_small.json --target-nodes 200 # -> 192 nodes
python3 scripts/bench/generate_grid.py grid_medium.json --target-nodes 1500 # -> 1,003 nodes
python3 scripts/bench/generate_grid.py grid_large.json --target-nodes 2200 # -> 2,605 nodes
```
Ports PGM's own C++ benchmark generator (`tests/benchmark_cpp/fictional_grid_generator.hpp`): a radial MV
feeder with stochastically-attached LV sub-grids. Not a bit-for-bit RNG replica, but both tools read the
*same* generated JSON, so the comparison is apples-to-apples. `--target-nodes` is approximate — LV attachment
is a Bernoulli process.
## 2. Benchmark gridoxide
```bash
cargo build --release --example bench_network
./target/release/examples/bench_network grid.json [repeat-count] [backend] [mode]
python3 scripts/bench/bench_gridoxide_native.py grid.json [backend] # pure-Python, always warm
```
- `repeat-count` (default 1) — re-runs the solve, for stable averages and for `perf record` sampling.
- `backend` — `scalar` (default), `block`, `klu`, `klu_native`, `pardiso`. See `docs/src/solvers/backends.md`.
- `mode` — `cold` (default, fresh symbolic factorization per repeat) or `warm` (one `PersistentSolver` reused).
Use `warm` for any comparison against another tool.
## 3. Benchmark power-grid-model
```bash
python3 -m venv .venv-pgm
.venv-pgm/bin/pip install power-grid-model # prebuilt wheel, no C++ build
.venv-pgm/bin/python3 scripts/bench/bench_pgm.py grid.json
```
Compare gridoxide's `total (guess + NR)` against PGM's `min`/`mean`. Voltage outputs (`voltage_mag min/max`
vs `u_pu min/max`) should match closely — a large mismatch means the comparison itself is broken.
### Results — synthetic radial grids
`newton_raphson`-only ms/run (200 warm repeats) vs PGM's `mean` (5 warm runs):
| 192 | 1.28 | 0.52 | 0.44 | 0.45 | 2.06 | 0.42 |
| 1,003 | 7.86 | 2.80 | 2.52 | 2.69 | 5.06 | 0.93 |
| 2,605 | 20.97 | 6.81 | 6.04 | 6.56 | 11.73 | 2.49 |
**PGM is faster than every gridoxide backend on this topology, warm-vs-warm — a real, standing gap.**
`KluNative` tracks `Klu` closely (1.02–1.09x). `Pardiso` is 2–4.7x slower than `Klu`: its default nonsymmetric
matching/scaling preprocessing carries a roughly size-independent fixed cost that dominates at these sizes.
Contrast with §4, where gridoxide's `Klu` is often *faster* than lightsim2grid's KLU-backed C++ solver on
transmission topologies — the result depends heavily on grid topology, not implementation language.
## 4. Benchmark against MATPOWER test cases
gridoxide against five independent Newton-Raphson implementations —
[PGM](https://github.com/PowerGridModel/power-grid-model),
[lightsim2grid](https://github.com/m-mirz/lightsim2grid) (C++/KLU),
[powsybl-open-loadflow](https://github.com/powsybl/powsybl-open-loadflow) via
[pypowsybl](https://pypi.org/project/pypowsybl/), pandapower, and
[VeraGrid](https://github.com/SanPen/VeraGrid) (numba-JIT) — on the same 12 cases lightsim2grid's own
benchmark uses: `case14`, `case118`, `case_illinois200`, `case300`, `case1354pegase`, `case1888rte`,
`case2848rte`, `case2869pegase`, `case3120sp`, `case6495rte`, `case6515rte`, `case9241pegase`.
```bash
git submodule update --init tests/data/benchmark-grids
python3 -m venv .venv-case-suite
.venv-case-suite/bin/pip install maturin numpy scipy power-grid-model pandapower lightsim2grid pypowsybl \
VeraGridEngine
VIRTUAL_ENV=.venv-case-suite .venv-case-suite/bin/maturin develop --release --features python,klu
.venv-case-suite/bin/python3 scripts/bench/run_case_suite.py --python .venv-case-suite/bin/python3
```
Loops all 12 cases, converting each MATPOWER `.m` to PGM JSON on first use (cached in
`scripts/bench/.case-cache/`, gitignored). A case that fails to convert or diverges gets an explicit
`FAILED (...)` cell with the tool's real exception. `--repeat` sets timed solves per case (default 10),
`--cache-dir`/`--out` set output locations.
Each tool runs at its closest-comparable settings: lightsim2grid with `SolverType.KLU` (its own benchmark
default), pypowsybl with the same "basic" `LoadFlowParameters` powsybl's benchmark repo uses, pandapower with
its defaults, VeraGrid with `SolverType.NR` and every automatic control (tap changers, remote voltage
regulation, alternate-solver fallback) off — see `bench_veragrid.py`'s docstring.
**Data paths.** gridoxide and PGM read MATPOWER `.m` files converted by `gridoxide.matpower`
(`python/gridoxide/matpower.py`; `matpower_to_pgm.py` is a thin CLI wrapper), *not* pandapower's importer.
Three cases (`case1888rte`, `case6495rte`, `case6515rte`) will not converge through pandapower's importer at
all: it assigns each bus a physical `baseKV`, which requires keeping every impedance/tap conversion referenced
to a consistent side. powsybl's own importer sidesteps this by assigning `nominal_v = 1.0` to every bus and
encoding off-nominal ratio as `rated_u1 = ratio, rated_u2 = 1.0` — MATPOWER's formulation needs a *consistent*
voltage reference, not a physical one. `matpower_to_pgm.py` does the same. lightsim2grid and pandapower load
via `pandapower.networks.<case>()` (they need the net object); pypowsybl and VeraGrid read the same `.m` file
through their own importers. See `gridoxide.matpower`'s module docstring for the full derivation.
Generators are modelled as genuine PV buses in every tool. gridoxide parses PGM's `voltage_regulator`
component the same way PGM's own `newton_raphson_pf_solver.hpp::set_u_ref_and_bus_types` does
(`src/pgm/mod.rs::pgm_to_buses_and_branches`).
Individual tools also run standalone: `matpower_to_pgm.py`, `bench_pgm.py`, `bench_lightsim2grid.py`,
`bench_pypowsybl.py`, `bench_pandapower.py`, `bench_veragrid.py`, `convert_pandapower_case.py`.
### Timing results
Warm-run mean, 5 timed calls on one persistent model per tool, ms:
| case14 | 15 | 0.042 | 0.032 | 0.059 | 0.027 | 0.074 | 0.155 | 0.028 | 1.985 | 14.547 | 3.083 |
| case118 | 119 | 0.211 | 0.164 | 0.106 | 0.111 | 0.340 | FAILED¹ | 0.149 | 4.969 | 15.445 | 6.389 |
| case_illinois200 | 201 | 0.664 | 0.296 | 0.257 | 0.269 | 0.649 | 0.585 | 0.329 | 6.816 | 15.741 | 7.131 |
| case300 | 301 | 1.709 | 0.582 | 0.490 | 0.512 | 1.525 | FAILED¹ | 0.528 | 9.578 | 19.078 | 11.577 |
| case1354pegase | 1355 | 6.665 | 2.969 | 2.393 | 3.167 | 3.626 | FAILED² | 2.575 | 38.163 | 22.155 | 48.832 |
| case1888rte | 1889 | 10.343 | 3.951 | 3.716 | 3.429 | 4.617 | FAILED¹ | FAILED³ | FAILED⁴ | 26.073 | FAILED⁶ |
| case2848rte | 2849 | 18.166 | 6.202 | 5.255 | 5.347 | 7.213 | FAILED¹ | 7.699 | FAILED⁴ | 30.648 | 87.872 |
| case2869pegase | 2870 | 18.856 | 7.366 | 6.109 | 6.731 | 8.270 | FAILED² | 6.135 | 97.341 | 31.530 | 99.871 |
| case3120sp | 3121 | 18.793 | 6.553 | 5.323 | 6.078 | 7.108 | FAILED² | 6.297 | 81.665 | 29.024 | 97.111 |
| case6495rte | 6496 | 61.226 | 19.387 | 16.721 | 19.994 | 20.931 | FAILED² | FAILED³ | FAILED⁴ | 111.096 | FAILED⁶ |
| case6515rte | 6516 | 73.762 | 22.021 | 18.578 | 22.537 | 25.251 | FAILED² | FAILED³ | FAILED⁴ | 49.822 | FAILED⁶ |
| case9241pegase | 9242 | 110.813 | 33.835 | 29.260 | 34.623 | 35.644 | FAILED² | 27.460 | 384.882 | 75.689 | 384.093 |
¹ PGM `IterationDiverge` (20-iteration default). ² PGM `SparseMatrixError`, raised during construction.
³ lightsim2grid `ac_pf` divergence (`V.shape[0] == 0`). ⁴ pypowsybl `MAX_ITERATION_REACHED`/`Unrealistic
state` under the same "basic" parameters powsybl's own benchmark repo uses — *without* the phase-shift-zeroing
workaround its `MatpowerUtil.java` applies to these same three cases (applying it makes powsybl converge on
all three in ~4 iterations). ⁵ needs `--features python,pardiso` and `MKLROOT`. ⁶ VeraGrid
`results.converged == 0` with `SolverType.NR` and automatic controls disabled.
These timings predate the correctness fixes below and were not fully re-measured; spot-checks agreed to
run-to-run noise (`klu_native` case14/case2869pegase/case9241pegase: 0.027/6.914/35.347 vs 0.027/6.731/34.623).
The fixes add shunts to existing Y-bus diagonals, changing neither sparsity nor iteration counts. Solved
*voltages* did change — see the accuracy section.
### Reading the timing table
**gridoxide converging on the three RTE cases is not a robustness result.** Its MATPOWER→PGM conversion rounds
`transformer.clock` to the nearest 60°, and every phase shift in these cases is under 30°, so all are zeroed —
4 branches in `case1888rte`, 17 in `case6495rte`, 16 in `case6515rte`. That is precisely the workaround
powsybl's benchmark applies (footnote 4). gridoxide is solving the same easier problem, having applied the
workaround implicitly during conversion. `check_matpower_residual.py` quantifies the cost (a five-figure MVA
residual on `case1888rte`). pandapower converges via its own native path.
**Backends.** `block` beats `scalar` 2–3x on the larger cases (case9241pegase: 33.84 vs 110.81 ms) despite
`scalar` using `faer` and `block`'s LU having no partial pivoting — factoring at 2×2-block granularity halves
the elimination steps. `klu` is ~10–15% ahead of `block`. `klu_native` (`crates/gridoxide-sparse/src/klu_native/`) converges to the
same voltages as every other backend on all 12 cases and lands within 1.0–1.3x of `klu`, after `perf` traced
an earlier 1.9–2x gap to allocator churn: `kernel::refactor_block` allocated two `Vec`s per column per Newton
iteration where real KLU's `klu_refactor.c` allocates nothing. Rewriting it to mutate in place against a
reused `RefactorScratch` dropped glibc allocator functions from ~32% of `newton_raphson` self-time to
effectively absent. `pardiso` matches the other backends' voltages exactly and sits 1.2–3.2x behind `klu`,
without the clean size trend the synthetic grids show — these cases vary more in topology and conditioning.
**vs lightsim2grid.** Roughly competitive: `klu` faster on most cases (case2848rte 5.26 vs 7.70 ms),
lightsim2grid faster on a couple (case9241pegase 29.26 vs 27.46 ms). Earlier `cold` numbers made gridoxide
look 1.3–1.7x slower, which `perf` traced to redone symbolic factorization — not a solver-speed gap. On
case9241pegase, warm cut per-solve `klu` time ~45%, with COLAMD/AMD/BTF ordering about a third of a cold solve.
**Python frameworks.** pypowsybl, pandapower and VeraGrid are markedly slower than every C/Rust-backed solver
here. VeraGrid's numbers are already post-JIT-warmup; its *first* call on a large case costs seconds.
Two fixes came out of building this table:
- `block` used to be `N/A` — it panicked on any `PV` bus, since its 2×2-per-bus indexing assumed `PQ`. Fixed
by giving a `PV` block a dummy `ΔVmag = 0` row (`solver::build_jacobian_blocks`), equivalent to the scalar
backend's dimension reduction. This surfaced a latent bug in `BlockLu::refactor`: it scattered
`adj.row(perm[j])` where it needed `adj.col(perm[j])` — silently wrong on any value-asymmetric matrix (the
real Jacobian is one; every unit test happened to use symmetric off-diagonal blocks). See
`solve_asymmetric_off_diagonal_matches_dense_reference`.
- `case3120sp` diverged catastrophically on *every* backend (`voltage_mag` ≈ -1548/930 p.u.), contradicting
this table. Bisected to `network::dc_angle_guess`, a DC-power-flow angle pre-pass added for CGMES FullGrid
but wired unconditionally into `linear_initial_guess`. Its own justification never panned out (FullGrid
still doesn't converge either way) and the full `cargo test --features cgmes` suite passes identically
without it, so it was removed outright.
PGM's numbers needed a workaround to obtain at all: once the converter began writing real `q_min`/`q_max` onto
`voltage_regulator`, all 12 inputs tripped PGM's `ExperimentalFeature` error through the public
`calculate_power_flow` (version 1.13.120 never wired `experimental_features` through the public wrapper).
`bench_pgm.py` calls the private `_calculate_power_flow` with `experimental_features="enabled"`, confirmed to
reproduce the same converged voltages the public API gave before Q-limits existed.
### Accuracy results
A MATPOWER `.m` file is power-flow *input*, not a solved snapshot — but it fully specifies the *equations*,
and any correct solution must satisfy them. This is the check to trust when it disagrees with the cross-tool
table:
```bash
python3 scripts/bench/check_matpower_residual.py --all
```
Rebuilds Ybus from the `.m` file (following `makeYbus.m`) and evaluates `dS = V .* conj(Ybus @ V) -
S_specified`, checking P at every PV/PQ bus and Q at every PQ bus. Absolute, no second tool involved.
```bash
python3 scripts/bench/accuracy_case_suite.py --out accuracy.md
```
Reports every other tool's per-bus voltage-magnitude deviation from **gridoxide's own** solution, matched by
MATPOWER bus number (each tool's ID convention confirmed empirically — see the script's docstring). Needs the
same environment as `run_case_suite.py` but runs all tools in one interpreter. **Using gridoxide as the anchor
is presentational and treacherous**: its five backends agreeing bit-identically says they share a front end,
not that the front end is right, and deviations *from* an anchor cannot distinguish "every other tool is off"
from "the anchor is off." Read alongside the residual check.
| case14 | n=14 median=0.0000% p90=0.0000% max=0.0000% | n=14 median=0.0024% p90=0.0139% max=0.0674% | n=14 median=0.0024% p90=0.0139% max=0.0674% | n=14 median=0.0024% p90=0.0139% max=0.0674% | n=14 median=0.0024% p90=0.0139% max=0.0674% |
| case118 | FAILED¹ | n=118 median=0.0000% p90=0.0001% max=0.1740% | n=118 median=0.0000% p90=0.0001% max=0.1114% | n=118 median=0.0000% p90=0.0002% max=1.3665% | n=118 median=0.0000% p90=0.0001% max=0.1114% |
| case_illinois200 | n=200 median=0.0000% p90=0.0000% max=0.0000% | n=200 median=0.0148% p90=0.0536% max=0.1531% | n=200 median=0.0767% p90=0.2268% max=1.0768% | n=199 median=1.3989% p90=3.0931% max=5.8083% | n=199 median=1.3989% p90=3.0931% max=5.8083% |
| case300 | FAILED¹ | n=300 median=0.0049% p90=0.1984% max=0.7585% | n=300 median=0.0007% p90=0.1696% max=0.7585% | n=300 median=0.1741% p90=4.5251% max=10.0565% | n=300 median=0.0112% p90=4.4378% max=10.0471% |
| case1354pegase | FAILED² | n=1354 median=0.0000% p90=0.0018% max=2.7037%⁵ | n=1354 median=0.0000% p90=0.0018% max=2.7037%⁵ | n=223³ median=1.5827% p90=4.2504% max=9.9509% | n=223³ median=1.5827% p90=4.2504% max=9.9509% |
| case1888rte | FAILED¹ | FAILED | FAILED | n=1790 median=0.8890% p90=3.1720% max=26.9861% | FAILED |
| case2848rte | FAILED¹ | FAILED | n=2848 median=0.0069% p90=0.1249% max=2.8521%⁵ | n=2767 median=1.4167% p90=4.5991% max=22.2243% | n=2767 median=1.5308% p90=4.9142% max=97.9196%⁴ |
| case2869pegase | FAILED² | n=2869 median=0.0000% p90=0.0011% max=2.7529%⁵ | n=2869 median=0.0000% p90=0.0009% max=2.7529%⁵ | n=889³ median=1.8198% p90=4.7823% max=11.8474% | n=889³ median=1.8198% p90=4.7823% max=11.8474% |
| case3120sp | FAILED² | n=3120 median=0.0004% p90=0.0100% max=2.4033% | n=3120 median=0.0130% p90=0.2114% max=1.7435% | n=3119 median=1.6588% p90=5.6333% max=27.6617%⁴ | n=3119 median=1.6588% p90=5.6333% max=27.6617%⁴ |
| case6495rte | FAILED² | FAILED | FAILED | n=6490 median=1.7131% p90=4.9194% max=81.7230%⁴ | FAILED |
| case6515rte | FAILED² | FAILED | FAILED | n=6510 median=1.8431% p90=5.2888% max=81.7200%⁴ | FAILED |
| case9241pegase | FAILED² | n=9241 median=0.0000% p90=0.0008% max=2.3527%⁵ | n=9241 median=0.0000% p90=0.0008% max=2.3527%⁵ | n=9240 median=2.4184% p90=6.8006% max=29.9954%⁴ | n=9240 median=2.4184% p90=6.8006% max=29.9954%⁴ |
¹ ² Same PGM failures as the timing table.
³ pandapower/lightsim2grid's `n` on the two PEGASE reductions is far below the real bus count. Not a naming
collision and not a constant offset — pandapower's bundled `case1354pegase()`/`case2869pegase()` derive their
per-bus numbering from the 9,241-bus network differently than the vendored `.m` files do. Unresolved, and
specific to those two cases' bundled data; every other case matches at full bus count.
⁴ pandapower/lightsim2grid `max` on the largest cases (22–98%) is far outside every other tool's (≤2.9%).
Both load `pandapower.networks.<case>()` rather than the vendored `.m`, and that data provably differs
(footnote 3), so part of this is a data difference. `check_matpower_residual.py` deliberately excludes them.
⁵ These remaining 2.3–2.9% maxima are endpoints of branches whose phase shift the conversion rounds to zero,
because PGM's `transformer.clock` can't represent a continuous angle (6 such branches in `case1354pegase` up
to 66 in `case9241pegase`). `check_matpower_residual.py --zero-phase-shifts` confirms it directly: gridoxide's
residual is nonzero on exactly these 7 cases and drops to exactly 0.0000 MVA against the phase-zeroed network.
pypowsybl and VeraGrid honour the real angle, so this is gridoxide's conversion, not their solve.
**The cross-tool table missed four real bugs.** `case14` used to read `median=0.313% p90=1.788% max=2.154%`
for *every* one of the five tools to three decimals. Five independent solvers reproducing each other's
deviation that precisely is not five tools being slightly wrong — it is the anchor being wrong. An earlier
version of this section instead concluded the gap was two valid solutions from independent Newton-Raphson
implementations. `check_matpower_residual.py` settles it without reference to any tool, and found:
1. **gridoxide ignored every PGM `shunt` on the PGM-JSON path.** `PowerFlowModel.from_pgm_json` and
`examples/bench_network.rs` built Y-bus from lines and transformers only, never calling
`pgm::pgm_shunts_1ph` + `network::stamp_shunts` — which `from_cgmes` beside it does call, so the helpers
were well covered while the real entry points skipped them. Discarded 1 shunt in `case14` (19 MVAr) up to
7,327 totalling ~89.5 GVAr in `case9241pegase`.
2. **The converter dropped branch charging on the transformer path.** A branch with off-nominal ratio *and*
nonzero `b` became a PGM `transformer`, which has no charging field. `makeYbus.m` folds it into
`Ytt = ys + j·b/2`, `Yff = Ytt/|tap|²`, which decomposes exactly into two `shunt` components — now emitted,
driving `case300`'s residual from 3.2 MVAr to 0.
3. **The converter flipped the sign of every negative reactance on the transformer path.** `uk` used
`np.hypot(r, x)`, which is unsigned, so series-capacitive branches (10 in `case3120sp`, ~40 per RTE case)
came back inductive. `network::transformer_admittances_ex` already handles signed `uk` correctly — 537 MVA
of residual on `case3120sp`.
4. **The converter dropped generator reactive output at PQ-typed buses.** `sym_gen.q_specified` was hardcoded
to `0.0`. Correct at a PV bus, wrong where MATPOWER types the bus PQ and `Qg` is a fixed injection
(`case1888rte` bus 1005 at -19 MVAr, `case2848rte` 2534 and 642). Residual there was exactly `-Qg`.
With all four fixed, gridoxide satisfies the as-published equations to 0.0000 MVA on 5 of 12 cases, and on the
other 7 once the reference's phase shifts are zeroed the same way the conversion zeroes them (footnote 5).
Cross-tool medians fell ~2 orders of magnitude (case9241pegase pypowsybl/VeraGrid: `median 0.628% p90 4.008%
max 20.838%` → `median 0.0000% p90 0.0008% max 2.3527%`), and PGM now agrees to 0.0000% where it converges.
`case_illinois200` is a cautionary note: it converged in both tools to visibly different voltages, with
VeraGrid landing much closer to PGM. That was read as "this case has multiple solution basins." The real cause
was bug 1 — gridoxide discarding four bus shunts (190 MVAr). Two tools agreeing against a third is evidence
about the third.
A separate pypowsybl-side matching bug in this benchmark: pypowsybl bus IDs look like `"VL-<n>_<k>"`, but
`<n>` is *not* the MATPOWER bus number for any bus that is only ever a transformer's secondary side —
pypowsybl names a secondary voltage level after its *primary* side's bus number with an incrementing suffix
(on case14, `"VL-4_1"` is really bus 7). Fixed by reading each `Line`/`TwoWindingsTransformer` id instead
(`"LINE-<from>-<to>"`, with `"#<k>"` for parallel branches) — see `pypowsybl_bus_to_matpower_id`. Confirmed by
VeraGrid, an entirely independent importer, matching pypowsybl's post-fix deviations to 4 decimals on every
shared case.
## 4b. Batched power flow — the multi-core CPU baseline
`bench_batch.py` measures `batch::BatchSolver`: one topology, many injection scenarios across cores. This is
the time-series/QSTS/Monte-Carlo shape, and exists to be the baseline any future GPU work must beat —
beating a *single-threaded* CPU solver is not a result.
```bash
maturin develop --release --features python,klu
python3 scripts/bench/bench_batch.py .case-cache/case9241pegase.json klu 256
```
Scenarios are ±20% uniform load scalings from a fixed seed. The script asserts voltages agree across every
thread count to <1e-12, so a scaling number is never reported without its correctness check.
**AMD Ryzen 7 250 (8 physical cores / 16 SMT threads, 16 MiB shared L3), 256 scenarios, `klu`.** Steady state
(third consecutive run, stable to ~2%):
| 1 | 28.04 | 1.00x | 2.384 | 1.00x |
| 2 | 15.75 | 1.78x | 1.369 | 1.74x |
| 4 | 10.33 | 2.72x | 0.914 | 2.61x |
| 8 | 7.91 | **3.55x** | 0.680 | **3.51x** |
| 16 | 7.71 | 3.64x (SMT) | 0.721 | 3.30x (SMT, *slower*) |
**Read absolute numbers on this machine with suspicion.** It is a thermally constrained laptop APU; a cold
first run reads far better than anything reproducible (the very first 8-thread case9241pegase measurement was
3.65x and never recurred). Every table here is steady-state, and every A/B below was taken *interleaved* —
old, new, old, new — because two numbers captured minutes apart measure the heatsink, not the code.
**Re-measure on a desktop or server part before quoting any of this.**
### Precomputed Jacobian offsets
`jacobian::JacobianPattern` derives the sparsity pattern and per-nonzero recipe once per topology, then refills
one reused `Vec<f64>` per iteration instead of rebuilding a `Vec<(usize, usize, f64)>`. Every `LinearSolver`
backend already discarded the `(row, col)` half after construction, so `factor_and_solve_values` hands over
just the values. Interleaved A/B, two rounds, min-of-3 each, ms/solve:
| before | 34.33 / 33.46 | 11.47 / 11.43 | 2.560 / 2.576 | 0.717 / 0.741 |
| after | 29.95 / 30.09 | 10.63 / 10.56 | 2.395 / 2.365 | 0.690 / 0.693 |
| **gain** | **11.4%** | **7.5%** | **7.3%** | **5.1%** |
Assembly measures at ~36% of iteration time, so single-threaded on the large case
roughly a third of that stage was allocation and index rebuilding. The gain is *smaller* at 8 threads: once
all cores run, the solve is memory-stalled in the LU, so assembly is a smaller share. It is also smaller on
the small case, whose triplet array (~20k nnz, ~0.5 MB) fits in cache; case9241pegase's is ~150k nnz, ~3.6 MB,
rebuilt every iteration. Values are bit-for-bit identical — `crates/gridoxide-core/src/jacobian.rs` compares `f64::to_bits`.
### Why sub-linear scaling is not a solver defect
Both cases land at ~3.5x on 8 physical cores (44% efficiency):
| 8 physical cores (16 logical is SMT2) | 8x | 8.0x |
| All-core clock throttle (2,977 → 2,236 MHz, from `/proc/cpuinfo` under load) | x0.75 | 6.0x |
| Shared-L3 / memory-bandwidth contention | remainder | **~3.5x observed** |
The last row was isolated by running **8 single-threaded processes** instead of 8 threads: separate address
spaces share no allocator arenas, locks or false sharing, so in-process contention would show up as better
process scaling. It did not — 3.96x processes vs 3.91x threads, back to back on an equally warm machine. The
ceiling is hardware.
An earlier draft claimed scaling degrades with grid size (63% on case1354pegase vs 46% on case9241pegase) and
blamed L3 overflow. That was an artifact of comparing a cold-machine run against a warm one; under equal
thermal conditions the two scale within 0.05x. The working-set effect is real but shows up in the Jacobian
assembly gain (11.4% vs 7.3%), not in thread scaling.
**Implication for the GPU plan.** The CPU batch path is memory-bandwidth bound, not compute bound — an
argument *for* the GPU direction, since bandwidth is what datacenter GPUs have in bulk. The number to beat
here is ~127 solves/s steady-state on case9241pegase, but a publishable claim needs a *server* CPU baseline.
## 4c. JAX oracle — validating the block-diagonal embedding
`jax_oracle.py` is an independent reimplementation of batched AC power flow in JAX (f64, CPU, dense), written
to answer questions the Rust solver cannot answer about itself. It is **not** a performance prototype — it is
deliberately the slowest power flow in this repo and must never be quoted as a speed number.
```bash
python3 -m venv .venv-jax
.venv-jax/bin/pip install jax numpy maturin
VIRTUAL_ENV=$PWD/.venv-jax .venv-jax/bin/maturin develop --release --features python,klu
.venv-jax/bin/python scripts/bench/jax_oracle.py .case-cache/case118.json 8
```
It consumes gridoxide's *own* Y-bus and bus arrays (`ybus_triplets`, `bus_spec`, `initial_guess`,
`zip_term_counts`) rather than re-deriving the model. That is the point: if the oracle built its own model, a
| case | buses | 1. oracle vs klu | 2. **BDE vs independent** | 3. oracle vs BatchSolver |
|---|---|---|---|---|
| case14 | 15 | 4.4e-16 | 4.4e-16 | 6.7e-16 |
| case118 | 119 | 1.1e-15 | 1.3e-15 | 1.2e-15 |
| case_illinois200 | 201 | 6.4e-15 | 2.7e-14 | 6.5e-11 |
| case300 | 301 | 1.5e-14 | 3.5e-14 | 6.5e-14 |
| case1354pegase | 1,355 | 3.9e-14 | 1.2e-13 | 7.5e-14 |
| case1888rte | 1,889 | 1.4e-13 | 4.6e-13 | 3.1e-12 |
**Column 2 is the one that matters.** Stacking B scenarios into one block-diagonal matrix and taking a single
LU should be equivalent to B independent solves — what lets an AMD GPU path work without a batched
refactorization API. Now checked
numerically: machine-precision agreement and exactly matching per-scenario iteration counts on every case.
Scope limits: constant-power injections only (ZIP terms asserted absent), dense Jacobian so B is auto-capped
and case9241pegase is out of reach, no Q-limit enforcement or island partitioning.
index, and the oracle "converged" in 4 iterations on a one-unknown problem. `bus_spec` now returns `Vec<u32>`
and the oracle validates array shape on load. Worth recording for what it implies about check 2: in that
broken run, BDE-vs-independent *passed* — both paths shared the same wrong indexing. A self-consistency check
between two code paths cannot detect a fault in what they share, which is why check 1 must pass first.
## 4d. Block-diagonal embedding on real sparse code
`crates/gridoxide-core/src/bde.rs` stacks B scenarios' Jacobians into one block-diagonal sparse matrix and takes a single
factorization per Newton iteration. §4c established equivalence with *dense* linear algebra; this extends it
to the sparse backends the GPU path mirrors, where the claim is far less obvious — KLU applies BTF ordering,
AMD permutation and partial pivoting to the stacked matrix, none of which know about the block structure.
```bash
cargo run --release --example bde_check -- scripts/bench/.case-cache/case1354pegase.json 16
```
16 scenarios, `klu_native`, versus independent per-scenario solves:
| case118 | 119 | 2,928 | 17,456 | **0** | 0/16 | 2.6 ms | 7.1 ms |
| case300 | 301 | 8,512 | 59,968 | **0** | 0/16 | 10.2 ms | 23.5 ms |
| case1354pegase | 1,355 | 39,184 | 253,552 | **0** | 0/16 | 47.4 ms | 130.6 ms |
| case2869pegase | 2,870 | 83,664 | 586,160 | **0** | 0/16 | 113.6 ms | 310.9 ms |
**Bit-exact, not merely within tolerance**, as theory predicts: BTF finds the B disconnected components, AMD
orders within each, and no fill crosses a block, so partial pivoting cannot reach across. Iteration counts
match per scenario, so the embedding doesn't perturb convergence. `tests/bde_test.rs` additionally checks that
the stacked pattern contains no entry linking two scenarios, and that a masked scenario's block is exactly the
identity *in the same stored positions* — masking must preserve sparsity, since dropping a converged
scenario's block would invalidate the cached symbolic factorization.
**It is ~2.7x slower on a CPU, and that is expected.** One large factorization beats B small ones only on
hardware that wants wide independent work — a GPU property. `bde.rs` is an architecture validator and the
host-side half of a GPU path, **not** a CPU optimization; use `batch::BatchSolver` for real CPU work.
## 5. Cross-validate CGMES import against pypowsybl
`cross_validate_cgmes_microgrid_be.py` checks gridoxide's CGMES import + solve against pypowsybl's own, on the
same ENTSO-E MicroGrid-BE-MAS conformance files `tests/cgmes_microgrid_be_test.rs` checks against published
`SvVoltage`. That test's doc comment claimed pypowsybl "also deviates by a comparable few percent" as a manual
finding; this script computes and asserts it.
```bash
pip install pypowsybl
python3 scripts/bench/cross_validate_cgmes_microgrid_be.py
```
Runs `examples/cgmes_microgrid_be_dump.rs` as a subprocess (predates `PowerFlowModel.from_cgmes`; §6 uses the
native binding), then zips BE-MAS + boundary files for pypowsybl and solves with the same "BASIC"
`LoadFlowParameters` used everywhere here. pypowsybl's OpenLoadFlow is pinned via `slackBusesIds` to
gridoxide's own slack so angles are comparable. Worst observed deviation: 0.22% voltage / 0.07° angle;
`--tol`/`--angle-tol` default to 1% / 0.3°.
## 6. Benchmark against CGMES conformance test configurations
gridoxide (`KluNative`) against pypowsybl/powsybl-open-loadflow, both importing and solving the *same* CGMES
profile files — no conversion step on either side. pypowsybl is the only vendored reference with native CGMES
import at all (neither PGM nor lightsim2grid has a single CGMES file in its tree).
```bash
maturin develop --release --features python,cgmes
python3 scripts/bench/bench_gridoxide_cgmes.py <fixture_name> <profile.xml>...
python3 scripts/bench/bench_pypowsybl_cgmes.py <fixture_name> <profile.xml>...
```
`bench_gridoxide_cgmes.py` uses `PowerFlowModel.from_cgmes` directly; `bench_pypowsybl_cgmes.py` zips profiles
into a temp archive first (pypowsybl's importer needs one archive). Both also report deviation from the
fixture's published `SvVoltage` (parsed by the shared `cgmes_sv.py`). On gridoxide's side this is exact —
`bus_index_for_mrid` looks a bus up by `TopologicalNode` mRID. On pypowsybl's side it's a heuristic:
`match_powsybl_buses_to_tn` reconstructs the mapping from the TP profile's `ConnectivityNodeContainer`
references (pypowsybl's bus-view ID is `"<container-mRID>_<index>"`), resolving multi-node containers by
nearest-magnitude matching against the *published* voltage — kept independent of gridoxide's own solve.
### Timing results
Warm mean (5 timed calls on one persistent model) plus each side's cold (construct+solve):
| PowerFlow | 2 | 2 | 0.007 | 1.398 | 0.55 | 19.26 |
| MiniGrid | 15 | 11 | 0.019 | 1.306 (no solve)¹ | 3.86 | 65.00 (no solve)¹ |
| PST_PhaseTapChangerTable_Type3 | 2 | 2 | 0.006 | 1.386 | 0.60 | 22.96 |
| MicroGrid-BE-MAS | 13 | 7 | 0.023 | 0.872 | 1.52 | 35.86 |
| MicroGrid-Type2-HVDC-MAS | 6 | 4 | 0.010 | 0.711 (no solve)⁵ | 0.91 | 27.33 |
| SmallGrid | 163 | 120 | 0.188 | 2.555 | 67.91 | 519.63 |
| Svedala | 191 | 104 | 0.188 | 2.511 | 65.71 | 543.74 |
| RealGrid | 6,252 | 5,806 | 18.258 | 136.959 | 1,449.74 | 5,890.94 |
gridoxide's node counts run higher because `PersistentSolver::solve` keeps every disconnected switchyard stub
as its own island with a placeholder solution, while pypowsybl's `connected_component_mode=MAIN` solves only
the largest component. gridoxide is faster on every fixture where both actually solve — excluding MiniGrid and
MicroGrid-Type2-HVDC-MAS, whose `(no solve)` cells measure time to hit an iteration limit or decide there is
nothing to calculate. On the rest the ratio shrinks with size, from ~200x on the 2-bus fixtures (where fixed
per-call overhead dominates) to ~14x on SmallGrid/Svedala and ~7.5x on RealGrid. Construction dominates
pypowsybl's cold figure far more: its Java CGMES importer costs tens to hundreds of ms even on 2-bus fixtures,
while gridoxide's cold barely exceeds its warm solve at small scale (PowerFlow: 0.55 vs 0.007 ms).
These are faster than an earlier version (RealGrid warm was 24.8 ms) because removing `network::dc_angle_guess`
(§4) eliminated one extra sparse solve per `linear_initial_guess` call — ~26% on RealGrid, whose ~6,000-bus
system was the most expensive to redo. The table predates the two importer fixes below and was not
re-measured; spot-checks gave RealGrid 17.61/1,425 and Svedala 0.172/70.2 with identical node counts.
### Accuracy results
Deviation from each fixture's published `SvVoltage` (`n` = `TopologicalNode`s matched to a solved bus):
| PowerFlow | 2 | 0.0000% | 0.0001% | 0.0001% | 2 | 0.0001% | 0.0001% | 0.0001% |
| MiniGrid | 11 | 1.248% | 2.334% | 2.387% | — | NOT CONVERGED¹ | — | — |
| PST_PhaseTapChangerTable_Type3 | 2 | 0.0001% | 0.0001% | 0.0002% | 2 | 0.0002% | 0.0002% | 0.0002% |
| MicroGrid-BE-MAS | 7 | 0.527% | 2.573% | 2.612% | 7 | 0.477% | 2.612% | 2.612% |
| MicroGrid-Type2-HVDC-MAS | 2² | 3.315%² | 3.315%² | 3.315%² | — | NO CALCULATION⁵ | — | — |
| SmallGrid | 127 | 0.002% | 0.217% | 0.447% | 7³ | 0.002%³ | 2.833%³ | 2.833%³ |
| Svedala | 108 | 0.213% | 1.000% | 2.865% | 96 | 0.285% | 1.627% | 4.071% |
| RealGrid | 6,051 | 0.018% | 0.424% | 93.475%⁴ | 5,806 | 0.099% | 0.822% | 82.068%⁴ |
¹ pypowsybl hits its Newton-Raphson iteration limit on MiniGrid under flat-start "BASIC" parameters. This cell
used to read `0.000%` with the non-convergence in a footnote — actively misleading, since 0.000% is the *best
possible score in the column*. Those figures came from the two boundary buses whose voltage never moves, i.e.
comparing the input SV profile against itself. `bench_pypowsybl_cgmes.py` now prints `NOT CONVERGED` on
stdout; it previously warned only on stderr, which the `2>/dev/null` used to produce these tables discarded.
² gridoxide's `n=2` reflects that this fixture's published SV set includes `TopologicalNode`s (DC-side detail
resolved away by `cgmes_resolve_dc_converters`, `crates/gridoxide-core/src/dc.rs`) with no AC bus in the solved output — not a
matching failure. This is the one row **not** re-measured after the two fixes below: reproducing `n=2` needs
the original profile-file set, and including the `-BD-MAS` boundary files the fixture requires yields `n=4`.
Treat these four figures as provisional.
³ pypowsybl's container-based bus reconstruction breaks down on SmallGrid: only 7 of 118 candidate containers
matched, because its TP profile groups far more `TopologicalNode`s per container (up to 9) than pypowsybl's
bus view does, so the join key mismatches. A limitation of this benchmark's matching heuristic on a fixture
with 838 `Disconnector`s + 427 `Breaker`s, not a pypowsybl accuracy issue. gridoxide's `n=127` has no such
uncertainty.
⁴ **RealGrid's `max` is a defect in the fixture**, now pinned down rather than inferred from tool agreement:
```bash
python3 scripts/bench/check_cgmes_sv_consistency.py RealGrid
```
Compares each two-winding transformer's published voltage ratio against what its own `ratedU` values declare,
adjusted by the actual `RatioTapChanger.step`. On RealGrid, 83 of 1,461 transformers miss their nameplate
ratio by >5%, and the four worst by **70–87%** — LV sides of four 63/20 kV units published at 10.76–11.88 kV
where the `TopologicalNode`'s `BaseVoltage`, the end's `ratedU`, and all five `SynchronousMachine.ratedU`
values say 20 kV. No tap position produces that; those are exactly the four buses at the top of gridoxide's
error list. Excluding the 140 buses on a >5%-inconsistent transformer, RealGrid becomes `n=5911 median=0.017%
p90=0.423% max=4.318%` — the entire 93% headline is fixture data, and the median and p90 barely move.
⁵ pypowsybl doesn't solve MicroGrid-Type2-HVDC-MAS at all: `run_ac` reports `status=NO_CALCULATION,
iteration_count=0, status_text="Network has no generator with voltage control enabled"`. CGMES
`VsConverter`/`CsConverter` aren't recognized as AC voltage-controlling sources, so `get_buses()` echoes back
the state its importer initialized from the input SV profile — both the 0.000% "deviation" and the timing are
artifacts. gridoxide's DC-aware solve (`crates/gridoxide-core/src/dc.rs`) is the only one of the two that solves this fixture.
Where both tools have comparable `n` (PowerFlow, PST, MicroGrid-BE, Svedala, RealGrid), accuracy against the
published reference is close — consistent with §5, and not with either tool being systematically more correct.
**Two CGMES importer bugs, found by working outward from the worst-deviating buses:**
1. **`build_two_winding` applied only one end's structural voltage ratio.** A `PowerTransformerEnd`'s `ratedU`
need not equal its bus's `nominalVoltage`, and that difference *is* an off-nominal ideal-transformer ratio
independent of any tap changer. The importer folded in end 1's `ratedU / bus1.u_rated` but dropped end 2's.
Invisible on every hand-authored fixture (all have `ratedU == nominalVoltage`), pervasive on real data:
**666 of RealGrid's 1,509** two-winding transformers and 4 of Svedala's 53, by up to 7.2%. Fixing it cut
RealGrid's median deviation **5x, 0.088% → 0.018%**.
2. **An out-of-service `SynchronousMachine` still regulated voltage.** The machine loop gated on terminal
connectivity and `RegulatingCondEq.controlEnabled` but never `Equipment.inService`, though
`equipment_in_service` was already applied to lines, switches and shunts in the same file. All six of
Svedala's `inService=false` machines also carry `controlEnabled=true`, and one (`_f4cde1f4`, `p=q=0`)
regulates a *remote* terminal, pinning an in-service bus to 21 kV against the published 20.134 kV —
Svedala's worst bus. `tests/cgmes_svedala_test.rs`'s de-energized count moved from 78 to the correct 83.
Net: RealGrid median 0.088% → 0.018%, p90 0.525% → 0.424%; Svedala p90 halved (2.121% → 1.000%) and max fell a
third (4.351% → 2.865%). Svedala's *median* moved the other way, 0.081% → 0.213% — about a thousandth of a
per-unit on the middle bus of a 108-bus fixture. Both changes are justified by the fixture's own equipment
data and both moved the tail substantially, but the median did not improve.
## 7. State estimation against power-grid-model
Every section above times a *power flow*; this one times a *state estimation*: gridoxide's two methods against
power-grid-model's two, on the same MATPOWER cases from byte-identical measurements. PGM is the only vendored
reference that does state estimation at all, so this is a two-tool comparison.
```bash
# 1. synthesize measurements and write the augmented documents both tools read
cargo run --release --example bench_se -- case14 case118 case300 case1354pegase case2869pegase \
--emit /tmp/se_bench
# 2. time both tools on them
VIRTUAL_ENV=$PWD/.venv-ls2g .venv-ls2g/bin/maturin develop --release --features python
.venv-ls2g/bin/python scripts/bench/bench_se_pgm.py
```
The measurement set is generated once by gridoxide and written *into* the input document, so neither tool gets
a set tuned to it: a `sym_voltage_sensor` on every node and a `sym_power_sensor` on both ends of every line and
transformer, ~4x redundancy. Values come from gridoxide's power flow for a blunt reason — PGM's power flow does
not converge on any of these converted cases (on case300 the deviation *grows* to 1394 after 200 iterations).
That is worth chasing separately; it doesn't tilt this comparison, since both tools estimate from the same
numbers.
**The data is perfectly consistent, with no noise at all.** That is the right shape for timing, where cost
follows measurement count and sparsity, and the wrong shape for judging robustness. Bad-data behaviour needs a
different harness.
### Timing results
Best of 5 warm calls on one persistent model per tool, ms. Both sides amortize symbolic factorization across
those calls — PGM through `PowerGridModel`, gridoxide through the `PersistentEstimator` that
`StateEstimationModel` holds:
| case14 | 14 | 94 | 0.2 | 0.2 | **0.1** | 0.1 |
| case118 | 118 | 862 | 0.7 | 0.5 | 1.2 | **0.5** |
| case300 | 300 | 1,944 | `SparseMatrixError`¹ | **0.7** | 6.3 | 1.3 |
| case1354pegase | 1,354 | 9,318 | `SparseMatrixError`¹ | **3.4** | 23.6 | 6.1 |
| case2869pegase | 2,869 | 21,197 | `SparseMatrixError`¹ | **8.4** | 72.2 | 16.5 |
Re-measured. The previous revision of this table gave gridoxide's Newton-Raphson as 4.3 / 24.5 / 16.5
on the last three cases; 16.5 was its *iterative-linear* figure duplicated into the wrong column, and
the harness reports 72.2 there. Newton-Raphson does not get faster from 1,354 buses to 2,869 — it
refactorizes every iteration, and that cost grows. The current figures include the measurement-model
caching in `se::functional`, worth 8-16% on the larger cases.
¹ PGM's Newton-Raphson estimator raises `SparseMatrixError` ("possibly singular matrix! ... might mean the
system is not fully observable") on every case from 300 buses up, on documents its *own* iterative-linear
method estimates from the same sensors without complaint. The suggested cause does not fit: gridoxide's
Newton-Raphson converges on the same measurements to 1e-14, and gridoxide's observability analysis flags
exactly one unobservable state per case — the virtual slack bus it appends per active source (`src/pgm/mod.rs`),
which sits at index `n_nodes`, carries no sensor by construction, is pinned rather than estimated, and is not
part of PGM's problem at all. Every physical node is observable. Reproduces from a plain `PowerGridModel`
construct-and-call with no gridoxide code in the path; not investigated further.
That extra bus is also why gridoxide's bus counts run one higher than the `buses` column (`bench_se.rs` prints
301 for case300), the same way §6's node counts exceed pypowsybl's.
**Newton-Raphson.** gridoxide is the only one of the two that answers at all above 300 buses, which outweighs
any timing statement about the cases where both run. Where both run they are within 2x either way: gridoxide
faster on case14 (0.1 vs 0.2 ms), PGM faster on case118 (0.7 vs 1.2 ms).
**Iterative-linear.** PGM is faster, by a ratio remarkably stable with size — 1.6x at case300, 1.8x at
case1354pegase, 2.0x at case2869pegase. This used to read that flatness as "a constant-factor gap in
per-iteration work, not an asymptotic one". **That inference was wrong, and backwards.** A stable ratio is
equally consistent with two ratios that happen to be stable, which is what this is.
```bash
.venv-ls2g/bin/python3 scripts/bench/se_iterations.py /tmp/se_bench
```
`se_iterations.py` measures the iteration count identically for both tools — the smallest `max_iterations`
that does not fail — having first checked the two convergence criteria are the same quantity (both take
`max over buses of |Δu|`, phase-normalized, against a 1e-8 default):
| case300 | 9 | 18 | 2.0x | **0.75x** |
| case1354pegase | 10 | 17 | 1.7x | **0.99x** |
| case2869pegase | 10 | 20 | 2.0x | **0.87x** |
**gridoxide's iterations are individually cheaper than PGM's; it takes about twice as many.** Both tools
convergence-rate gap, and the linear algebra the earlier reading pointed at is already ahead.
Those counts were 29 / 28 / 33 until the under-relaxation was allowed to *recover* rather than only ever
back off — the flat-start transient alone tripped the backoff at iteration 3, after which a whole run ran
at half a step. About 40% fewer iterations and 16-22% less wall clock, which is why the table above now
reads 1.5-1.8x where it read 2.0-2.3x, and why case118 has crossed over to gridoxide being faster.
Taking the damping out entirely does not work: undamped, gridoxide's map locks into a limit cycle at
~1.7e-1 rather than converging. `docs/src/state_estimation/iterative.md` has the trace, localizes the
instability to the branch-terminal power rows (dropping them takes the undamped step from 1.7e-1 to
7.1e-4; dropping the voltage rows changes nothing), and rules out the two obvious mechanisms — the `|U|²`
weight scaling PGM deliberately omits, and the zero-injection KKT constraints it has no equivalent of.
Why PGM needs no damping at all is still open.
**Accuracy is not a differentiator on this data.** Where PGM answers, the two tools agree to 5.8e-15 V on
case14 and 3.3e-15 V on case118. Against the state the measurements were read from, gridoxide's
Newton-Raphson lands at ~1e-14 p.u. throughout and its iterative-linear at 8.7e-10 through 1.1e-6, degrading
with size — the linearization bias inherent to the method, documented in
`docs/src/state_estimation/iterative.md`, not an implementation artifact.
A caveat this table used to carry is now retired: gridoxide's column no longer includes a fresh symbolic
factorization per call. `StateEstimationModel` holds a `crate::se::nr::PersistentEstimator` whose cache
survives across `solve()` calls, so the symbolic phase is paid once per model on both sides.
### Sanity-checking against `examples/bench_se.rs`
`examples/bench_se.rs` times gridoxide's two methods against *each other* and reports 117 ms for
Newton-Raphson on case1354pegase where the table above reports 24.5 ms. Both are correct; they are not the
same problem. `bench_se.rs` also synthesizes bus-injection measurements (11,189 scalar rows vs the emitted
document's 9,318), and an injection row is a full Y-bus row rather than a two-bus branch row, so its gain
matrix is substantially denser. The emitted document drops them because PGM has no counterpart sensor.
**Compare numbers within one harness, never across the two.**
## 8. DC (Bθ) power flow and sensitivity factors
`examples/bench_dc.rs`, on the MATPOWER cases from §4, converted the same way:
```bash
python scripts/bench/matpower_to_pgm.py tests/data/benchmark-grids/matpower/case9241pegase.m /tmp/case9241pegase.json
cargo build --release --example bench_dc
./target/release/examples/bench_dc /tmp/case9241pegase.json 20 solve
./target/release/examples/bench_dc /tmp/case1354pegase.json 5 ptdf # or lodf
```
### Timing results
`solve` is one full `dc_power_flow` per repeat — branch reduction, island partition, one
factorization per island, back-substitution and branch flows. The AC column is
`examples/bench_network.rs` in `warm` mode with the `scalar` backend, i.e. the *fastest* AC
configuration, so the ratio is a conservative statement of what DC buys.
| case118 | 119 | 187 | 0.046 ms | 0.197 ms | 4.3x |
| case1354pegase | 1,355 | 1,992 | 1.268 ms | 8.550 ms | 6.7x |
| case9241pegase | 9,242 | 16,050 | 13.40 ms | 106.4 ms | 7.9x |
The gap widens with size, which is what the two methods' structure predicts: DC pays one
factorization while AC pays one *per iteration* (4–5 of them here), and the per-iteration Jacobian
refill grows with the network too.
Sensitivity sweeps produce every column once — `ptdf` over every bus, `lodf` over every branch —
against a single factorization built in `DcSensitivity::new`:
| case118 | 0.389 ms (119 cols) | 3.3 µs | 0.609 ms (177 cols) | 3.4 µs |
| case1354pegase | 30.5 ms (1,355 cols) | 23.3 µs | 45.5 ms (1,430 cols) | 31.8 µs |
| case9241pegase | — | — | 2,848 ms (14,384 cols) | 198 µs |
The case9241pegase LODF sweep produces 14,384 columns from 16,050 branches; the remaining 1,666
are radial or carry no DC flow, and return `None` without a back-substitution.
The per-column figure is the one worth quoting, and it is the entire argument for
`sparse::RealFactorization`: these are `n` triangular solves against **one** factorization. Routing
them through `solver::LinearSolver`, whose contract refactorizes on every call, would have made a
PTDF sweep `n` factorizations rather than one — at case1354pegase that is the difference between
22.5 µs and something on the order of a full DC solve, per column.
Note that neither sweep materializes a dense matrix. At case9241pegase a dense PTDF is 1.19 GB and
a dense LODF 2.06 GB; the column API exists so that cost is optional.
`linear::batch::DcBatchSolver` runs many bus-injection scenarios against **one** numeric
factorization for the whole batch — `B` depends only on topology, so a scenario changes nothing but
the right-hand side. `./target/release/examples/bench_dc <case> 2000 batch`:
| case118 | 0.046 ms | 0.0051 ms | 9.0x |
| case1354pegase | 1.268 ms | 0.0373 ms | 34x |
| case9241pegase | 13.40 ms | 0.2423 ms | 55x |
The ratio grows with size because the factorization is what gets hoisted, and it is the part that
scales worst. This is a stronger guarantee than `batch::BatchSolver` can offer on the AC side, which
reuses only the *symbolic* half — Newton's Jacobian changes numerically at every iteration of every
scenario, so the numeric factorization cannot be shared.
N-2 screening via `DcSensitivity::multi_outage_flows`, over every pair drawn from the first 200
branches (`./target/release/examples/bench_dc <case> 200 n2`). Removing `k` branches is a rank-`k`
update, so each pair costs `k` triangular solves plus one `k × k` dense solve — never a
refactorization:
| case118 | 17,391 | 15,502 | 0.0050 ms | 0.046 ms | 9.2x |
| case1354pegase | 19,900 | 9,017 | 0.0966 ms | 1.268 ms | 13x |
| case9241pegase | 19,900 | 18,432 | 0.3955 ms | 13.40 ms | 34x |
The unsolvable pairs are *breaking sets* — removing both branches disconnects the network, so no
post-outage flow exists and `multi_outage_flows` returns `None` without a solve. They are the case
single-branch screening misses entirely: two parallel lines are each individually non-radial, but
together they may be the only path.
### Accuracy results
DC has no convergence tolerance — it is a direct solve — so the meaningful quantity is the
factorization residual `max |Σ(flows out of bus) − P_bus|`, reported as
`DcSolution::max_residual`:
| case118 | 119 | 1.3e-12 |
| case1354pegase | 1,355 | 2.8e-10 |
| case9241pegase | 9,242 | 1.3e-9 |
All round-off, growing with problem size as expected. Two observations from this run worth
recording:
- **case9241pegase contains 16 branches with negative reactance** — series capacitors, which
`dc_power_flow` reports through `DcSolution::negative_reactance_branches` rather than passing over
in silence. They make `B` indefinite without making it singular. This is the case that justifies
`dc_branches` guarding the PGM-link sign on the *exact admittance constant* rather than on
`x < 0` generally: a blanket clamp would have silently corrupted 16 real branches here.
- **Residual is the early-warning signal for clamped zero-impedance branches.** A clamped branch
carries `b ≈ 2.8e5` against ordinary branches' `b ≈ 5–50`, spreading `B`'s condition number by
four orders of magnitude — and unlike Newton, DC has no iterative refinement to hide that behind.
On `tests/data/pgm/powerflow/link/dummy-test`, which is full of them, the residual still comes in
under 1e-9 (asserted in `tests/dc_powerflow_test.rs`).
### Cross-validation against MATPOWER's own DC formulation
`check_dc_matpower.py` rebuilds the DC system straight from each `.m` file following `makeBdc.m`
and `dcpf.m`, then pushes *both* angle vectors through the same `Pf = (Bf·Va + Pfinj)·baseMVA` so
the comparison lands on branch flows in MW and needs no mapping between MATPOWER's branch rows and
gridoxide's flat branch index. No second tool is used as a reference — this is the same
self-contained approach as `check_matpower_residual.py`, and it is what the in-repo oracle
(`dc_flows_match_the_lossless_ac_limit`, which checks DC against gridoxide's *own* AC branch-flow
code) structurally cannot provide.
```bash
python scripts/bench/check_dc_matpower.py --all --zero-phase-shifts --ignore-gs
```
| case14 | 14 | 20 | 2.2e-14° | 2.1e-13 MW | 147.8 MW |
| case118 | 118 | 186 | 4.0e-13° | 4.4e-12 MW | 450.0 MW |
| case_illinois200 | 200 | 245 | 3.8e-13° | 4.0e-12 MW | 371.8 MW |
| case300 | 300 | 411 | 1.6e-13° | 5.7e-12 MW | 1292.0 MW |
| case1354pegase | 1,354 | 1,991 | 2.3e-12° | 3.0e-10 MW | 1504.8 MW |
| case1888rte | 1,888 | 2,531 | 1.2e-11° | 2.9e-10 MW | 1456.7 MW |
| case2848rte | 2,848 | 3,776 | 2.1e-12° | 7.1e-11 MW | 1442.0 MW |
| case2869pegase | 2,869 | 4,582 | 6.8e-12° | 1.8e-10 MW | 1585.6 MW |
| case3120sp | 3,120 | 3,693 | 1.8e-11° | 5.1e-10 MW | 850.2 MW |
| case6495rte | 6,495 | 9,019 | 6.9e-12° | 3.7e-10 MW | 1527.0 MW |
| case6515rte | 6,515 | 9,037 | 1.4e-11° | 3.1e-10 MW | 1519.5 MW |
| case9241pegase | 9,241 | 16,049 | 1.6e-10° | 1.1e-9 MW | 1947.1 MW |
**Exact agreement on all twelve** — the worst disagreement across the suite is 1.1e-9 MW against
peak flows of ~1,900 MW, i.e. round-off. gridoxide's Bθ is the same formulation MATPOWER computes,
including the phase-shift and transformer-ratio conventions.
The two flags each isolate one known, documented model difference, and both leave a real
discrepancy visible when omitted:
- **`--ignore-gs`.** MATPOWER folds bus shunt conductance into `Pbus` as a constant real load at
|V| = 1; gridoxide's DC deliberately ignores all shunts (`docs/src/powerflow/dc.md`, Scope).
Without the flag, exactly the three cases with nonzero `ΣGs` disagree and no others — case300
(ΣGs = 1.3 MW) by 1.30 MW, case2869pegase (9.9) by 2.45 MW, case9241pegase (56.9) by 13.2 MW.
That one-to-one correspondence is what identifies the gap as *this* modeling choice rather than a
solver error, and it quantifies what adopting MATPOWER's convention would be worth.
- **`--zero-phase-shifts`.** `gridoxide.matpower` rounds `angle` to the nearest 60-degree clock
position, which for every shift in these cases means dropping it. Without the flag the
disagreement reaches **5.1e4 MW** on case1888rte — and its max Δangle comes out at 9.92°, which is
that case's known 9.95-degree shifter. This is the same conversion loss `check_matpower_residual.py`
documents for AC, and it bites far harder here, because a phase shift enters the DC right-hand
side directly rather than as a second-order term. It is a limitation of the PGM conversion, not of
the DC solver: `tests/dc_powerflow_test.rs` exercises phase shifters natively and matches the AC
limit to 1e-6.
## 9. AC N-1 contingency screening
`examples/bench_contingency.rs`, sweeping single-branch outages through
`batch::BatchSolver::solve_contingencies` and comparing against solving each outaged network
independently — the same work with no factorization reuse, which is what a naive contingency loop
does.
```bash
cargo build --release --example bench_contingency
./target/release/examples/bench_contingency /tmp/case9241pegase.json 60 klu_native
```
### Timing results
**Both sides run single-threaded**, so this isolates factorization reuse from parallelism; the
`solve_contingencies` sweep parallelizes on top of what is shown here.
| case118 | 187 (all) | 10 | 0.198 ms/outage | 0.393 ms/outage | 1.98x |
| case1354pegase | 200 | 65 | 4.823 ms/outage | 10.008 ms/outage | 2.08x |
| case9241pegase | 60 | 2 | 42.17 ms/outage | 114.50 ms/outage | 2.72x |
The ~2x floor is what the repo's own "symbolic factorization is ~45% of solve time" figure predicts,
and the ratio grows with size because that share does. The reuse here is narrower than an ordinary
batch's: each outage has its own Y-bus, so the Jacobian's cached per-nonzero recipe must be
re-analyzed every scenario (O(nonzeros), `PersistentSolver::invalidate_admittances`). What survives
is the symbolic factorization — the fill-reducing ordering and elimination structure — because
`network::build_ybus_with_outages` keeps the outaged branch's structural entries and so leaves the
sparsity pattern bit-for-bit intact.
The **severing** column counts contingencies that disconnect the network. Those cannot use that
trick, because `connected_components` reads the same structural entries and cannot tell a zeroed one
from a live one; they are detected up front with `network::structural_component_count` and fall back
to a full rebuild with a fresh factorization. That fallback is what lets an islanded contingency come
back as an honest `IslandStatus::NoReferenceBus` rather than a singular solve — and note it is not
rare: 65 of case1354pegase's first 200 branches sever something.
### Accuracy results
`tests/ac_contingency_test.rs` checks every single-branch contingency on both committed PGM fixtures
against a direct solve of the network built without that branch — same status, and voltages agreeing
to the solver's own `1e-6` tolerance. It runs at one thread as well as four: a single worker means
every scenario shares one `PersistentSolver`, which is the sharpest test of the cache-reuse claim,
and it is the configuration that caught the one real bug in this work (a `JacobianPattern` carried
across contingencies silently solves the *previous* outage's network, because the pattern caches the
admittance each entry was analyzed against).
Cross-tool comparison against another package's contingency module is not covered. The direct-solve
oracle above is stronger for the question that matters here — whether the fast path equals the slow
one — but it says nothing about agreeing with, say, lightsim2grid's `SecurityAnalysis` on the same
case.
## 10. Node-breaker topology import
`examples/bench_node_breaker.rs`, over the four genuine node-breaker configurations in
`tests/data/CGMES-Test-Configurations/v3.0/`. EQ+EQBD+SSH only — the reader does not need TP.
```bash
cargo build --release --features cgmes --example bench_node_breaker
D=tests/data/CGMES-Test-Configurations/v3.0/SmallGrid/SmallGrid-Merged
./target/release/examples/bench_node_breaker $D/SmallGrid_EQBD.xml $D/SmallGrid_EQ.xml $D/SmallGrid_SSH.xml
```
### Timing results
| MiniGrid | 103 | 90 (90) | 11 | 4.10 ms | 0.10 ms (2.5%) |
| FullGrid | 48 | 29 (28) | 8 | 7.45 ms | 0.11 ms (1.5%) |
| SmallGrid | 1,369 | 1,266 (1,203) | 115 | 59.89 ms | 3.25 ms (5.4%) |
| Svedala | 1,179 | 1,464 (988) | 128 | 57.54 ms | 2.72 ms (4.7%) |
Reading the switch graph costs 1.5–5.4% of the XML decode that has to happen anyway, and deriving a
bus view from it is well under a millisecond even at 1,464 switches. Import cost is not what makes
node-breaker support expensive.
### What it costs to keep switches
The AC unknown count under a constrained formulation is roughly `2·buses + 2·retained_edges` — two
per non-slack bus, two per retained edge. That is the number that decides whether a retention policy
can be a default:
| MiniGrid | 13 buses, ~26 unknowns | 43 buses, 30 retained, ~146 | 103 buses, 90 retained, ~386 |
| FullGrid | 22 buses, ~44 | 45 buses, 26 retained, ~142 | 48 buses, 29 retained, ~154 |
| SmallGrid | 167 buses, ~334 | 540 buses, 373 retained, ~1,826 | 1,369 buses, 1,266 retained, ~5,270 |
| Svedala | 228 buses, ~456 | 638 buses, 857 retained, ~2,990 | 1,179 buses, 1,464 retained, ~5,286 |
Two things fall out, one confirming the expected cost of `RetainAll` and one correcting a prediction.
**Confirmed: `RetainAll` really is a ~16x blowup.** SmallGrid goes from ~334 unknowns to ~5,270,
against a predicted ~5,264. `MergeAll` staying the default is not a stylistic preference.
**Corrected: `RetainAdjacentToBusbar` is not "a retained count in the tens".** That holds on the small
models (26 on FullGrid, 30 on MiniGrid) and fails on the real ones — 373 on SmallGrid and 857 on
Svedala, the latter retaining *more* switches than SmallGrid despite having fewer nodes, because
Svedala is the most switch-dense model in the tree. It still buys a 3x reduction against `RetainAll`
on SmallGrid, so it remains the right default *for a bus-breaker view*; it is just not the small
number expected, and a contingency campaign wanting a handful of switches should use
`RetentionPolicy::Explicit` rather than assume this policy is already small.
### Accuracy results
There is no tolerance here — the bus view is a partition, so it is either the exporter's or it is
not. `tests/cgmes_node_breaker_test.rs` compares the EQ+SSH-derived partition against the
`TopologicalNode` partition the same files publish, which is free ground truth requiring no fixture
authoring. See that file's header for the full table; the two results worth repeating:
- **MiniGrid agrees exactly.** 103 connectivity nodes and 90 switches reduce to precisely the 13
topological nodes the exporter published.
- **Every disagreement elsewhere is explained.** Where gridoxide merges further than the exporter
(FullGrid, SmallGrid) it is a closed switch the export left unreduced — the case
`merge_closed_switches` was written for. Where gridoxide merges *less* (34 cases across FullGrid,
SmallGrid and Svedala) the nodes are joined only by switches that are open or out of service, and
the exporter merged across them anyway. The test asserts that cause rather than the count, so it
would fail if the reader ever missed an edge.