Feature comparison: gridoxide vs. reference tools

A survey of five independent power-flow implementations — what each actually supports (verified against source/docs, not assumed), compared against gridoxide's own current scope — used to decide what gridoxide tackles next. Three (lightsim2grid, power-grid-model, powsybl-open-loadflow) are full local checkouts under references/ (itself gitignored, hence this file living at the repo root instead); see each tool's own CLAUDE.md/README for how to consult them further. The other two, VeraGrid (the GridCal successor) and pandapower — both also used as comparison tools in scripts/bench/run_case_suite.py — aren't checked out under references/; they're verified instead by reading their installed packages' own source directly (pip install VeraGridEngine pandapower; see each package's own directory structure for the file paths cited below). This file is a snapshot, not a living document, and will drift as gridoxide and all five tools evolve.

Scope of the most recent revision. Six rows were added for analysis types the table had no row for at all — remedial action optimization, voltage stability, harmonics, reliability, protection coordination, and input-format coverage — and the single "dynamic simulation" row was split into three, because RMS, EMT and small-signal analysis are three different machineries and one checkbox understated the gap. Absence of a row is how a gap goes untracked; that was the point of the exercise.

For the added rows, the three tools checked out under references/ (lightsim2grid, power-grid-model, powsybl-open-loadflow) were searched directly, and each ❌ below names what was searched for. VeraGrid and pandapower could not be re-surveyed — the previous revision read their installed packages and neither is installed here any more — so they are marked not surveyed rather than ❌, per the convention the "Multi-island" row introduced: absence of a survey is not evidence of absence of a feature. Their cells in the three split dynamics rows carry forward the directory names the previous survey recorded, which is evidence rather than inference.

Earlier revisions. gridoxide's own column was re-verified against current source (every cell claiming support names the function or type implementing it), and the "CGMES / CIM import" row was checked across all five comparison tools by counting CGMES/CIM-named files in their installed trees. The other five tools' cells in every row predating this revision were not re-surveyed — treat them as the older snapshot.

Summary table

Per-tool notes

lightsim2grid (C++/Python, KLU-backed)

• Solvers: NR (single-slack and distributed-slack variants), Gauss-Seidel (+ "synch"), DC, fast-decoupled (XB/BX). Linear-solver backend is pluggable (Eigen SparseLU, KLU, NICSLU, CKTSO) via SolverType.
• Elements: lines, 2-winding transformers (fixed tap ratio + phase-shift angle, changeable only between solves), shunts, loads, static generators, storage, DC lines/HVDC. No 3-winding transformers, no SVC, no switches (explicit TODO in SubstationContainer).
• Its own docs/disclaimer.rst is refreshingly explicit about what it doesn't do: no Q-limit enforcement, fixed taps mid-solve, steady-state only, symmetric only.
• ContingencyAnalysis and TimeSerie batch classes reuse Ybus factorization across many solves rather than rebuilding from scratch — same idea as gridoxide's PersistentSolver, just applied to a batch-of-scenarios use case rather than only repeated single-topology solves.
• Ingests grids from pandapower and pypowsybl/IIDM directly (gridmodel/from_pandapower, gridmodel/from_pypowsybl).

power-grid-model (C++/Python)

• Calculation types: power flow (sym + asym), state estimation (WLS, with observability checks), short-circuit (IEC 60909, phase-domain). No sensitivity/OPF.
• PF solver algorithms: Newton-Raphson (default), iterative-current, linear/linear-current (auto-selected when all loads are constant-impedance).
• No PV bus type in plain power flow ("not supported yet" per its own docs) — PV-like behavior instead comes from the newer voltage_regulator component, which fixes |U| and solves for Q; q_min/q_max exist on it but the automatic PV→PQ switching is explicitly flagged as not fully implemented.
• Same hand-rolled block-sparse LU architecture gridoxide's Block backend mirrors: per-bus 2×2 real blocks for NR power flow, full pivoting within a block only (no cross-block pivoting), pivot perturbation off by default for ordinary power flow (confirmed at the newton_raphson_pf_solver.hpp call site — this is what caused the SparseMatrixErrors investigated earlier this session).
• Batch calculations reuse the prebuilt topology graph and matrix prefactorization across scenarios when only load/gen/source setpoints change (not when topology/tap/shunt status changes) — the same invariant PersistentSolver::reset() documents for gridoxide.
• TapChangingStrategy outer loop (disabled by default): any_valid_tap, min_voltage_tap, max_voltage_tap, fast_any_tap.
• validate_input_data/validate_batch_data exist but are explicitly not run automatically for performance reasons — recommended for debugging, not the hot path.

powsybl-open-loadflow (Java, RTE)

• Calculation types: AC power flow, DC power flow, sensitivity analysis (AC+DC, incl. post-contingency), security/contingency analysis (N-1/N-k, AC+DC). No short-circuit, no state estimation.
• Solvers: Newton-Raphson (primary), Newton-Krylov, fast-decoupled — all pluggable via AcSolverFactory (service-loader based, genuinely extensible). Five voltage-initialization strategies (flat, warm/previous, uniform, DC-angle-based, magnitude-based).
• Most feature-rich of the three on voltage/reactive control: automatic PV→PQ switching with reactive capability curves, remote voltage control (one generator regulating a different bus), shared voltage control among multiple controllers, a priority scheme (generators > transformers > shunts), and even secondary voltage control (research-based).
• Distributed slack: on generators, loads, or "conform" loads; manual or automatic slack-bus selection with multiple strategies (first, largest-generator, most-meshed, named); also area-interchange-based distribution.
• Genuinely modular outer-loop architecture — OuterLoop/OuterLoopContext/OuterLoopResult abstractions, extensible via ServiceLoader, with 14+ concrete outer loops (distributed slack, area-interchange, reactive limits, transformer voltage/reactive-power control, phase control, shunt voltage control, secondary voltage control, HVDC AC-emulation limits). This is the architecture responsible for nearly every "extra" feature above the bare NR solve.
• Contingency analysis performance claim is best substantiated for DC specifically (Woodbury-formula fast path, WoodburyEngine/WoodburyDcSecurityAnalysis) — AC contingency/sensitivity analysis is documented as reusing full-resolve-style computation, and its own README's "Contributing" section flags AC performance as an open area, so the tool's reputation for contingency-analysis speed is strongest for DC, not universal.
• Supports asymmetric/unbalanced modeling (LfAsym* classes) and full node-breaker topology with connectivity traversal (NodeBreakerTraverser).
• Uses powsybl-math's LUDecomposition/MatrixFactory abstraction; native KLU via JNI is the primary path (same library gridoxide's own Klu backend vendors directly).

VeraGrid (Python, SanPen/VeraGrid, the GridCal successor)

• By far the broadest simulation-category scope of the five — installed as the headless VeraGridEngine package (not the Qt-GUI-bundled VeraGrid package), its Simulations/ directory alone has 25+ top-level categories: beyond power flow, also OPF (linear + nonlinear AC), state estimation, short-circuit, contingency analysis, sensitivity (PTDF/LODF), continuation power flow (PV curves), stochastic/Monte Carlo analysis, reliability analysis, investment/expansion-planning evaluation, net/available transfer capacity (NTC/ATC), topology reduction, and — uniquely among all five — electro- magnetic-transient (EMT) and RMS time-domain dynamic simulation with small-signal stability analysis. pandapower rivals it in raw feature count (see below) but has no equivalent of EMT/RMS dynamic simulation at all; lightsim2grid/PGM/powsybl are themselves narrower, purpose-built power-flow-focused engines by comparison.
• Solvers: among the most pluggable of the five via SolverType — NR, Gauss-Seidel, Fast-decoupled, Levenberg-Marquardt, Iwamoto-NR, Powell's Dog Leg, HELM (Holomorphic Embedding), Decoupled-LU, plus linear/ linear-AC modes and dedicated linear/nonlinear OPF solver types, all in one PowerFlowOptions.solver_type enum — though pandapower's own algorithm parameter is comparably broad and additionally offers a backward/forward-sweep solver neither VeraGrid nor any of the other four tools here have.
• Voltage/reactive/tap control is a set of independent boolean flags on PowerFlowOptions (control_q, distributed_slack, control_remote_voltage, control_taps_modules, control_taps_phase, orthogonalize_controls) applied inside the NR iteration itself, not a modular outer-loop registry the way powsybl's OuterLoop abstraction is — closer in spirit to power-grid-model's ad hoc tap optimizer than to powsybl's extensible architecture.
• Devices include 3-winding and generic N-winding transformers, switches (with CIM/IIDM node-breaker import), HVDC lines, VSC, and UPFC, plus a ControllableShunt device (stepped Bmin/Bmax regulating a bus's voltage to a setpoint) filling the SVC role neither lightsim2grid, PGM, nor powsybl have — broad FACTS coverage, though pandapower's own dedicated create_svc/create_tcsc/create_ssc set turns out broader still (see below).
• MATPOWER import (parse_matpower_file) reads each bus's type column and each generator's Vg setpoint directly, so genuine PV-bus modeling comes for free with no PGM-voltage_regulator-style conversion step — see scripts/bench/bench_veragrid.py.
• Its own numerical kernels are numba-JIT-compiled (first call per process pays a multi-second JIT-compilation cost unrelated to the power-flow algorithm itself — scripts/bench/bench_veragrid.py's warm-up call absorbs this) and its sparse LU solve is SciPy's SuperLU (scipy.sparse.linalg._dsolve._superlu.gstrf, wrapped in a numba-jitted custom CSC type, Utils/Sparse/csc2.py) — not pluggable across multiple sparse backends the way lightsim2grid or powsybl are.
• On the 12-case real-MATPOWER benchmark (scripts/bench/README.md), converges on 9 of 12 (the same three hard RTE cases every tool but gridoxide/pandapower also fails on) and lands roughly on par with pypowsybl — markedly slower than the C/Rust-backed solvers here, consistent with being a general-purpose Python framework rather than one optimized around raw repeated-solve throughput.

pandapower (Python, e2nIEE/pandapower)

• Also very broad in scope, though — unlike VeraGrid's from-scratch simulation engines — pandapower's own numerical power-flow/OPF path is largely a thin, numba-accelerated wrapper around PYPOWER (pandapower/pypower/, itself a Python port of MATPOWER), with pandapower supplying the richer network model (switches, controllers, 3-winding transformers, FACTS devices) and everything else (contingency, timeseries, estimation, shortcircuit, diagnostic) as sibling top-level packages built on top of that core.
• Solvers (runpp(algorithm=...)): "nr" (default, PYPOWER's Newton-Raphson, numba-accelerated), Iwamoto-NR ("maybe slower... but more robust" per its own docstring), backward/forward sweep ("bfsw", specially suited to radial/weakly-meshed networks — a solver category none of the other five tools here offer), Gauss-Seidel, and two explicitly separate fast-decoupled variants ("fdbx"/"fdxb"), plus HELM.
• Switches (create_switch/create_switches, bus-bus/bus-line/bus-trafo) are core to how pandapower represents topology at all, not a bolted-on extra the way they are for some other tools here — closest in spirit to powsybl's node-breaker model among the tools surveyed.
• FACTS-device coverage (create_svc/create_tcsc/create_ssc/create_vsc*) is the broadest of the six tools surveyed, including a thyristor-controlled series capacitor (TCSC) none of the others model.
• The generic Controller/BasicCtrl framework (control/basic_controller.py, driven by run_control=True) is genuinely extensible — any third-party code can subclass Controller and register it on net.controller — closer to powsybl's outer-loop extensibility in spirit than to PGM's/VeraGrid's fixed option flags, even though the underlying convergence-loop architecture isn't identical.
• contingency module includes a run_contingency_ls2g variant that offloads the actual repeated solves to lightsim2grid for speed — a real cross-tool dependency between two of the tools surveyed here, not just a coincidental feature overlap.
• OPF is two-tiered: a native, no-external-dependency PDIPM-based AC/DC OPF (runopp/rundcopp, inherited from PYPOWER) for standard formulations, plus an optional bridge to Julia's PandaModels.jl (runpm.py) for more advanced formulations (storage, multi-stage, etc.) when that external toolchain is installed.
• diagnostic() (diagnostic/diagnostic_helpers.py) is a real, generic input-validation/consistency-check function (disconnected elements, implausible parameter values, wrong reference system, ...) — closer to PGM's validate_input_data than to VeraGrid's format-specific-only import validation.
• Also has a dedicated protection package (protection-device/relay-coordination modeling) that none of the other five tools here have any equivalent of — outside the scope of this table's rows, but worth noting as another area where pandapower's breadth exceeds a pure power-flow-engine comparison.
• This is the same pandapower already used elsewhere in this benchmark suite (bench_pandapower.py, bench_lightsim2grid.py's and lightsim2grid's own init_from_pandapower) — see Backends and Factorization Reuse and scripts/bench/README.md for its own timing numbers, where it's the only tool besides gridoxide to converge on all 12 real MATPOWER cases.

Where gridoxide already exceeds or matches

• Asymmetric power flow: already solving and tested (matches PGM/powsybl/VeraGrid/pandapower; lightsim2grid doesn't have this at all).
• 3-winding transformers: already have a passing fixture (matches PGM/powsybl/VeraGrid/pandapower; lightsim2grid doesn't have this at all).
• Factorization reuse across repeated solves (PersistentSolver): conceptually identical to what lightsim2grid's ContingencyAnalysis/TimeSerie and PGM's batch-calculation path rely on.
• Batched solving over one topology (batch::BatchSolver): the API layered on top of that reuse, and the shape time-series/QSTS and Monte Carlo runs actually need — thousands of independent scenarios over an unchanging topology, spread across cores on rayon's own pool, one cached symbolic factorization per worker. Matches PGM's batch-calculation path and lightsim2grid's TimeSerie on the injection-scenario case; still short of both on contingency, which needs per-scenario topology (see gap 2 below). (bde::solve_batch_block_diagonal stacks a batch into one block-diagonal factorization instead, validated bit-exact against independent per-scenario solves in scripts/bench/README.md §4d — but it is ~2.7x slower on a CPU and exists to validate a future GPU path's architecture, so it is not a batching capability this table should credit.)
• Block-sparse LU backend granularity: matches PGM's own per-bus 2×2 block design, and gridoxide's faer-backed solve handles pivots PGM's own hand-rolled solver refuses (no pivot perturbation) on the same real transmission-scale data — still true after the converter fixes below, which changed PGM's input but not its failure pattern (same 6 SparseMatrixError / 4 IterationDiverge cases as before).
• Sparse-solver breadth: five backends (Scalar/Block/Klu/KluNative/Pardiso — the count previously read "four" while listing five) already exceeds VeraGrid's and pandapower's single fixed-solver paths, though it's still short of lightsim2grid's runtime-pluggable KLU/Eigen/NICSLU/CKTSO selection.
• CGMES import depth: one of four tools here with any CGMES/CIM import at all, and the only one of those four that is otherwise a focused AC power-flow library rather than a general-purpose framework. On the 8 conformance configurations benchmarked in scripts/bench/README.md §6 it is faster than pypowsybl on every fixture where both actually solve, and solves MicroGrid-Type2-HVDC-MAS, which pypowsybl declines to attempt (iteration_count=0, "Network has no generator with voltage control enabled").
• Multi-island solving: solves every connected component with per-island status rather than only the main one.
• Remedial action optimization: the one row where none of the five comparison tools competes — the reference is powsybl-open-rao, a separate project. And it is the most externally-validated thing here: the reference publishes its expectations as a Cucumber suite, so 156 of its own scenarios are read unmodified from the benchmark-grids submodule and scored assertion by assertion (tests/rao_cucumber_test.rs), rather than gridoxide checking its own arithmetic against itself.
• Solution verification tooling: scripts/bench/check_matpower_residual.py checks a solved case against the MATPOWER file's own power-flow equations, and check_cgmes_sv_consistency.py checks a CGMES fixture's published SvVoltage against its own EQ/SSH data. Neither needs a second tool as a reference. This is a benchmark-harness capability, not an input-validation feature — it does not close the "Input validation" row above, which is about validating input before a solve (PGM's validate_input_data, pandapower's diagnostic()).

Identified gaps, ranked by how often reference tools flag them as important

1. Q-limit enforcement / PV→PQ switching — every one of the five either has it, half-has it, or explicitly disclaims not having it as a known limitation. gridoxide's Bus already carried q_min/q_max, unused. Done — outerloop::ReactiveLimits implements the standard MATPOWER-style one-directional PV→PQ switching outer loop, tested in tests/q_limits_test.rs across all three Jacobian backends; PgmVoltageRegulator now parses PGM's own q_min/q_max fields. Opt-in: plain newton_raphson/PersistentSolver::solve are unchanged, so no existing test/benchmark behavior shifted. It was solver::newton_raphson_enforcing_q_limits, a standalone entry point, until the outer-loop layer made it composable with the two below.

2. Contingency/N-1 batch analysis — the one gap with most of its machinery already standing: PersistentSolver's factorization reuse and batch::BatchSolver's across-scenario parallelism are exactly what lightsim2grid, powsybl, VeraGrid (the most comprehensive, with linear, nonlinear, and HELM-based variants), and pandapower (which even offloads some of its own contingency solves to lightsim2grid for speed) build this on. Done for DC, which is the mode these tools actually screen contingencies in: linear::sensitivity::DcSensitivity::outage_flows gives post-outage flows from one triangular solve with no refactorization (and reports radial branches, which have no redistribution factors), while linear::batch::DcBatchSolver runs injection scenarios against a single factorization for the whole batch — 9x to 55x per scenario against independent solves, see scripts/bench/README.md §8. N-2 and N-k are covered too, via the generalized (Woodbury) formulation in multi_outage_flows — 0.40 ms per pair on case9241pegase against a 13.4 ms re-solve — including is_breaking_set for outage sets that disconnect the network, which single-branch screening misses by construction. AC contingency is done as well, via batch::BatchSolver::solve_contingencies. The old objection recorded here — "a branch outage gives each scenario its own Y-bus and therefore its own sparsity pattern" — turned out to be half right: the Y-bus does differ per scenario, but the pattern need not, because network::build_ybus_with_outages takes a branch out while keeping its structural entries. The symbolic factorization therefore carries across a whole N-1 sweep; only the Jacobian's cached admittances are re-derived, at O(nonzeros). Measured at 2.0x (case118) to 2.7x (case9241pegase) against independent solves single-threaded, before any parallelism. Contingencies that genuinely sever the network cannot use that trick — connected_components cannot distinguish a zeroed structural entry from a live one — so they are detected up front and fall back to a full rebuild, which is what lets an islanded contingency report NoReferenceBus honestly instead of a singular solve.

3. Distributed slack — 4 of 5 tools have it (lightsim2grid, powsybl, VeraGrid, pandapower); real transmission grids split slack across several generators, and a single slack absorbing an unscheduled few hundred megawatts distorts every branch flow around it. Done — outerloop::DistributedSlack, an outer loop in the same shape as the Q-limit one: solve, compare what the slack actually produced against what its own p_spec schedules it for, move the difference onto the participating generators by normalized weight, solve again.

   Two things are worth recording. The slack's schedule is its own p_spec, a field ordinary power flow ignores for a slack bus because the slack's output is an answer rather than an input — so a document that never needed it may leave it at zero, in which case the entire output is redistributed. And the outer loop is cheap in a way the Q-limit loop is not: only p_spec changes between passes, so bus types, n_unknowns and the Jacobian's sparsity pattern all hold and PersistentSolver keeps its symbolic factorization throughout.

   Convergence is linear, not one-and-done. The first-order term cancels exactly in a single pass — distributing \(\Delta\) adds \((1-\alpha_s)\Delta\) to the other participants while the slack's own schedule rises by \(\alpha_s\Delta\) — but the leftover is the change in transmission losses from the redistributed flows, which shrinks by a few per cent per pass rather than vanishing. Measured at tolerance = 1e-8: case14_ieee, case30_ieee and case118_ieee each take seven passes, moving 2.3, 2.4 and 16.5 per-unit off their slacks. Tested in tests/distributed_slack_test.rs, with every check re-deriving injections from network::power_injections rather than reading the loop's own bookkeeping.

   newton_raphson_distributing_slack is now outerloop::DistributedSlack, one loop among the ordered list of gap 11's entry, which is what lets it run together with Q-limit enforcement — the two used to be alternatives. Both are reachable from gridoxide solve --distribute-slack --enforce-q-limits and from PowerFlowModel.solve(distribute_slack=True, enforce_q_limits=True); the "library-only" note this paragraph used to carry is out of date.

   Area interchange is done too, as outerloop::AreaInterchange. It is the same mechanism generalized from one balance target to several, which is why powsybl's AcAreaInterchangeControlOuterLoop constructs a DistributedSlackOuterLoop as its no-area fallback — and here that relationship is an assertion rather than a remark: one area with a zero target reproduces distributed slack's per-bus shifts and solved voltages exactly.

   Building it turned up something the earlier sketch of this paragraph missed. Both ends of a tie line are measured into the branch, so the positions sum to the tie's losses rather than cancelling, and a set of agreed net positions — which sums to zero — is unachievable by exactly that amount. The knobs are one aggregate schedule per area, \(N\); the conditions wanted are \(N\) positions plus the slack on its own schedule, \(N+1\). Over-determined by one, always. So the slack's own area is the dependent one: its target is not enforced, it absorbs the tie losses, and the report names it and states the residual rather than implying it.

   Both importers are wired now. CGMES supplies both halves — ControlArea.netInterchange for the schedule, TieFlow for the boundary — and membership is derived by cutting the tie branches and taking connected components, because CGMES never lists what is inside an area. MicroGrid is the tree's genuine two-area fixture and the loop puts BE on its declared −236.977 MW exactly. UCTE supplies membership from its ##Z country codes and no schedule, which is all any UCTE file states; the twelve-node case is a real four-country interconnection. Both reach gridoxide solve --area-interchange.

   Two things the fixtures corrected. A TieFlow names the terminal at the boundary — the X-node two areas' lines meet at — not inside the area, so the seed is the cut branch's far end; seeding from the named end gave five contested buses and one area owning nothing. And the sign: CGMES states an import, AreaDefinition wants an export, checked against SmallGrid's own published state sitting at 210.271 MW against its declared 210.

   IIDM is wired too, and is the only importer that states membership directly — <iidm:area> lists its voltage levels, where CGMES gives only a boundary and UCTE a country per node. The twelve-node case exists in both UCTE and IIDM and their areas agree to 1e-6 (BE +2000 MW, DE −2500, FR +1000, NL −500), by entirely separate code on each side.

   One thing real data forced. AreaDefinition::uniform participates Slack and PV buses only, inherited from distributed slack — right for frequency response, wrong for a net position, which is met by redispatch, and a generator on a fixed active set-point is exactly what gets redispatched. powsybl's own two_area_case.xiidm has both of AREA2's machines at voltageRegulatorOn="false", so under uniform that area had no participant and its −400 MW schedule was unreachable. AreaDefinition::by_generation weights by p_spec instead and the CLI uses it; AREA2 then reaches −400.000 exactly.

   Still absent: nothing in the tree states a net position for a UCTE network, and PGM JSON has no area concept at all. See Area Interchange Control.

4. DC power flow as a first-class mode — every reference tool treats this as a basic offering. Done, and in the doing it turned out the old single "DC / linear power flow" row was conflating two unrelated algorithms, which is why the summary table above now carries two. The premise that linear_initial_guess was "most of the way there" was only true for one of them: that function is power-grid-model's CalculationMethod.linear (complex, constant-admittance, keeps resistance and produces |V|), and it is now exposed as linear::impedance::linear_power_flow. What was genuinely missing was real Bθ — the thing lightsim2grid, powsybl-open-loadflow and pandapower's rundcpp mean by DC — and that is new code in linear::btheta, with both of powsybl's approximation variants, phase-shifter support, per-island factorization, and PTDF/LODF sensitivities on top. See docs/src/powerflow/dc.md and linear_impedance.md; tested in tests/dc_powerflow_test.rs and tests/dc_sensitivity_test.rs.

5. Switches, HVDC and SVC as first-class model elements — no longer absent, but reachable only through CGMES import: switching state, converter setpoints and SVC regulation are all fixed at import time, with no element in gridoxide's own network model to change between solves. That is exactly the shape lightsim2grid's disclaimer calls out for its own fixed taps, and it is what stands between the current support and the contingency/time-series work in item 2.

6. Short-circuit calculation — Done. This was previously listed in item 6 below as outside gridoxide's scope; it is not any more. crates/gridoxide-core/src/shortcircuit/ implements IEC 60909 in the phase domain: all four fault types, both voltage-scaling choices, bolted and impedance faults, several simultaneous faults, and de-energized islands, with results reported in phase quantities and symmetrical components. The premise that made it affordable is that gridoxide was already sequence-in, phase-assembled — Line3Ph carries r0/x0, transformer_seq_params returns (y0, y1, y2), and network::fortescue_to_phase was already building the phase-domain Y-bus from them — so the new code is mostly fault boundary conditions on top of tested assets.

   Building it closed three genuine gaps in the shared three-phase model, each found by a fixture rather than by inspection: transformer_seq_params supported only two winding pairs and panicked on the rest (it now implements power-grid-model's general zero-sequence algorithm, covering YNd, zigzag and the low-susceptance rule); Line3Ph had no conductance field, silently dropping the tan δ dielectric-loss term; and the three-phase path dropped link components entirely, which quietly de-energizes whatever hangs off one.

   Validated against all 15 of power-grid-model's short-circuit fixtures. Eleven match to its own 1e-8; four cannot, because power-grid-model changed its zero-sequence regularization in November 2025 and never regenerated fixtures dating from 2023 — an inconsistency inside the reference set rather than a choice gridoxide can make differently. tests/data/pgm/short_circuit/README.md and KnownDivergence in the test record which, why, and by how much.

7. Sensitivity analysis — Done for both models, and it was previously listed in the item below as outside gridoxide's scope. The DC factors (PTDF, LODF, and the Woodbury N-k generalization) were already there under item 2; what is new is ac_sensitivity::AcSensitivity, which differentiates a converged AC operating point directly.

   The premise that made it cheap is that the derivative needs no new matrix: \(dx/dp = -J^{-1}\,\partial g/\partial p\) uses the same Jacobian Newton already assembles, and for an injection the right-hand side is a unit vector, so one triangular solve answers for the whole network. The \(\partial f/\partial x\) row for a branch flow is branch_flow::terminal_flow_derivs, written for the state estimator and reused unchanged. The only genuinely new numerics are the tap derivatives, which reduce to a scaling of each admittance entry by itself.

   Both directions are offered against one factorization — forward (one solve per variable, "what does this move?") and adjoint (one solve per monitored quantity, "what would move this?") — with sparse::RealFactorization::solve_transpose reusing the same LU rather than factorizing the transpose. Variables are active and reactive injection, transformer ratio and phase-shifter angle; functions are branch P/Q at either terminal and bus voltage magnitude or angle.

   Validated against a central-difference re-solve of the full nonlinear power flow, which shares no code with the module. One finding worth recording: a tapped branch's own flow carries a direct \(\partial f/\partial p\) term on top of the chain rule, and dropping it flips the sign of that branch's own sensitivity (\(-0.276\) to \(+0.286\) on distribution-case) while leaving every other branch correct — a failure mode that looks entirely plausible in the output.

   Still absent: AC contingency sensitivities (an outage is a finite topology change, not a derivative — BatchSolver::solve_contingencies re-solves instead) and differentiation through the outer loops.

8. DC-OPF (opf::dc, opf feature) — least-cost dispatch as a convex QP, with locational marginal prices and binding-limit shadow prices, behind a solver-independent boundary. Two backends sit behind it: an in-house interior-point QP solver needing nothing installed, and HiGHS as a cross-checked reference. Validated against analytic cases, KKT certificates, pglib-opf's published DC baseline (all five cases inside 0.03%), and solver-vs-solver agreement.

   AC-OPF (opf::ac) as well: the full nonconvex problem over generator P and Q, bus voltage magnitudes and angles, with apparent-power branch limits, solved by gridoxide's own nonlinear interior-point method (opf::nlp — line search, adaptive regularization, gradual barrier) on the injection Hessians of injection_hessian. All five pglib cases match the published AC objectives to 0.001% at violations of 1e-9 or better; derivatives are checked against finite differences and prices against a numerical d(cost)/d(load). A KKT point of a nonconvex problem is locally optimal, which is what every AC-OPF tool reports. CLI (gridoxide opf --ac) and Python (ac_opf). See Optimal Power Flow.

   Still absent: taps and phase shifters as decision variables, unit commitment, and security-constrained OPF.

9. Remedial action optimization — Done, and it is not a row any of the five comparison tools fills: the reference here is powsybl-open-rao, a separate project from open-loadflow. src/rao/ carries the CRAC model and readers, the evaluation kernel, an LP over range actions, a search tree over network actions, the CASTOR perimeter decomposition, automaton simulation, MNEC soft constraints, conditional usage rules, RA usage limits and the curative stop criterion, with an AC re-validation stage on top.

   What makes it unusual among the items here is the gate: the reference ships its expectations as a Cucumber suite, so 197 of its own scenarios are read unmodified from the benchmark-grids submodule and scored — 1830 of 1862 checkable assertions, 180/186 DC, 276/282 plus 963/977 AC, 123/123 second preventive and 288/294 costly optimization, with nothing left unsupported out of what was 177 steps. Forty-odd real defects were found by it, every one internally consistent and externally wrong, and several only after a confident diagnosis had to be withdrawn. See The Remedial Action Problem for the mathematics.

   All 32 assertions that do not match are recorded disagreements rather than a backlog: each has been measured on the reference's own objective and in none is gridoxide worse — its BestTapFinder rounding a set-point blind to a virtual cost its own javadoc warns about, and two where gridoxide secures the network more cheaply than the reference asks. The last open defect closed by reading the reference's own source rather than guessing an eighth mechanism.

   Genuinely absent: loop flows, relative margins, HVDC and counter-trade range actions, the angle/voltage monitoring pass, MARMOT multi-timestamp and CNE XML export — and an integer tap, which was measured as unjustified rather than assumed: all 24 scenarios whose configuration asks for APPROXIMATED_INTEGERS already match with continuous optimization and rounding, including six the reference built as controlled pairs against their own continuous twins.

10. The analysis types with no row until this revision — voltage stability, harmonics, reliability, protection coordination, and the three dynamics rows. None of them is a refinement of something already built; each is a different question asked of the network, and they are worth ranking rather than lumping together:

    • Voltage stability / continuation power flow was the cheapest by a distance, and is now built (crates/gridoxide-core/src/continuation/) — a predictor–corrector loop around the Newton solver that already handled Q-limits, distributed slack and multi-island, answering the question nothing else in the set answers: how much further can this be loaded? Q-V curves are built too (crates/gridoxide-core/src/qv.rs), so the row is complete. The estimate held: the reusable parts were the Jacobian pattern, the reactive-limit outer loop and the sparse backends, and the genuinely new numerics are small. What the estimate missed is that two of the three outer loops must not be run inside the corrector — they move p_spec as a function of the solved state, which corrupts the tangent — and are folded into the loading direction analytically instead. dynawo-algorithms offers voltage margin and load-increase-to-collapse as an external cross-check.
    • Quasi-static time series is next, and half-built: BatchSolver already runs many scenarios over one topology in parallel. What is missing is chronology — carrying storage state of charge, tap positions and controller memory from one step to the next — plus a driver and a result writer. simbench ships the profiles.
    • RMS dynamics was the largest, and is now built (src/dynamics/) — the DAE, its integrator, discrete events, the machine/exciter/governor/stabilizer library, three readers, and the surfaces. The estimate that the corpus and methodology would come close to free held, and better than expected: Dynawo's repository ships its own solver's reference outputs alongside each example, so the external gate needed no Dynawo install at all — a sparse clone of the example files was enough. What the estimate missed is where the real risk lay. It was not the integrator, which is a hundred lines and gated by closed forms; it was the conventions. Published statements of the subtransient machine disagree about the sign of ψ_2q, and saturation is stated incompatibly by PSS/E and Dynawo. Both had to be settled by derivation and by declining to guess, not by copying a source. And the one difference no self-consistency check could ever have found — the ω ≈ 1 stator approximation, worth 0.6% of terminal power per 0.9% of speed deviation — was found by reading Dynawo's Modelica after its trajectory diverged from ours. That is precisely what an external reference is for.
    • EMT, harmonics, reliability and protection coordination stay genuinely far off. Each needs modelling gridoxide has no foundation for — sub-cycle three-phase stepping, per-frequency Y-buses, outage-rate data, relay curves — and none is currently pointed at by any plan.

    See Resources for the datasets and reference implementations each would be validated against.

11. Transformer tap control, and the outer-loop layer under it — Done. This entry used to lump both in with "materially larger undertaking or outside gridoxide's current scope", and singled out static taps as a modelling assumption four of the five comparison tools do not make. They are one job, and the second half turned out to be the load-bearing one: newton_raphson_enforcing_q_limits and newton_raphson_distributing_slack were standalone entry points, so a caller could have at most one — not an architecture nicety but two shipped features that could not run together. Both are outerloop::OuterLoop implementations now, driven over an ordered list by powsybl's own schedule, with TransformerVoltageControl and PhaseControl as the third and fourth.

    Two things the fixtures corrected, both recorded as assertions so they cannot quietly stop being true. Svedala's own published solution violates all eleven of its declared deadbands by 2–5%, so the obvious fixed-point gate — taps must not move, since SSH and SV agree — had a false premise, and a loop that left those taps alone would be the defective one. And a discrete tap often cannot reach a deadband at all: PST_PhaseTapChangerLinear_Type2 asks for zero flow within 5e-4 pu on a shifter whose closest position leaves 0.11 pu, which needs an outcome of its own rather than being reported as success.

    Found on the way: CGMES transformer ordering was nondeterministic — ends_by_pt is a HashMap and Rust randomizes its iteration order per process, so every flat branch index derived from the transformer list came out differently on each run of the same program against the same file. Nothing asserted on one, so it never failed.

    Still absent: TCSC/SSC and the wider FACTS set; powsybl's two continuous-then-round voltage strategies; reactivePower-mode tap control; shunt-section control; tap control inside a contingency or batch sweep, and in DC. And the outer-loop layer is deliberately internal — the list is built by the crate, not discovered — so the row above reads ⚠️ rather than ✅. See Outer Loops and Transformer Tap Control.

Note on state estimation

Done in both domains — snapshot and batched, with voltage, power and current sensors, symmetric and per-phase. Within that scope gridoxide matches the most capable reference tool and leads it in two places. What is left of the three gaps this note used to list is one part of one of them, at the end. An earlier version claimed parity outright, which overstated it; a later one left the summary row saying "symmetric only" for some time after the phase domain worked, which understated it in the one place a reader looks first.

se::nr::estimate is Gauss-Newton on the normal equations, validated against power-grid-model's own state-estimation fixtures (committed under tests/data/pgm/state_estimation/ with their MPL-2.0 license files): per-unit magnitudes agree to 1.5e-9 on transmission-case — measured, against a standing assertion of 1e-6, since a gate is a floor and not a claim — and every sparse backend produces the same answer, since the gain matrix is an ordinary square system. See the State Estimation chapter.

Both analyses VeraGrid is credited with above are present. Observability (se::observability::analyze) separates structural from numerical unobservability and names the buses and quantities involved, rather than only reporting that a factorization failed. Bad-data detection (se::bad_data::analyze) runs the chi-squared test and identifies culprits by largest normalized residual. Zero injections are enforced as hard equality constraints rather than as high-weight pseudo-measurements — the approach that avoids the ill-conditioning power-grid-model has two fixtures named after.

Both of power-grid-model's calculation methods are implemented and agree with each other: Newton-Raphson (se::nr) and the prefactorized iterative_linear (se::iterative), selectable per call. link is modelled now (stamped as a branch, see the zero-impedance chapter), so the fixtures using one are reachable. Pseudo-measurement augmentation — filling an unobservable region with forecast values, which VeraGrid does — is not implemented; gridoxide reports the unobservable set instead, which is the prerequisite for it.

The two leads

• Bad-data detection, which power-grid-model does not have at all. Checked against its own documentation rather than assumed: it reports a per-sensor residual and stops there — no chi-squared test, no identification of a culprit. se::bad_data::analyze does both.
• Newton-Raphson robustness. power-grid-model's Newton-Raphson estimator raises SparseMatrixError on every benchmark case from 300 buses up, on documents its own iterative-linear method estimates from the same sensors without complaint, and that gridoxide's Newton-Raphson converges on to 1e-14. See scripts/bench/README.md §7.

The gaps that remain

Measured against power-grid-model 1.13 (references/power-grid-model/), in order of how much they matter:

1. Asymmetric state estimation. Done, and reachable. What remains is narrower than the heading suggests, and is a modelling difference rather than a missing capability.

   A phase-domain document now estimates end to end. SeNetwork::from_3ph builds the measurement model for a 3N-bus network, pgm::pgm_3ph_maps supplies the object-ID maps a sensor needs to resolve against it, and pgm::measurement::measurements_from_pgm_3ph maps sensors onto phase-expanded targets. tests/se_three_phase_test.rs estimates power-grid-model's transmission-case in the phase domain and matches the answer it published for that network solved asymmetrically — all 33 phase-buses to 1e-6, with angles agreeing up to the single rotation nothing measures.

   Nothing in the estimator changed for it. Target, StateLayout, the Jacobian, the constraints, both methods and the batch solver carry over untouched, because a three-phase branch terminal is a six-coefficient CurrentFunctional where a scalar one has two. A branch is indexed 3·branch + phase to match the 3·node + phase bus convention, which is what keeps Target identical between the two domains.

   That carry-over is now gated rather than asserted. It was true, but it had been a claim about function signatures — observability, bad-data detection and the batch solver take the same arguments in either domain, so they ought to work — and the reference handling is exactly where a symmetric assumption would have hidden, since three per-phase rotations are three symmetries where StateLayout removes one. tests/se_three_phase_test.rs now checks each of them in the phase domain: observability finds every physical bus determined and only synthesized ones not; bad-data detection names the corrupted phase by a factor of five over the runner-up; and a batch of five scenarios matches a sequential loop bit for bit.

   Two things came out of writing those gates. pgm::se::se_case_from_pgm_3ph is the single entry point, because standing the phase-domain model up took eight calls in a particular order against four converters and three stamping functions — which is why the estimator worked for some time with nothing outside the test suite able to call it. gridoxide estimate --asymmetric is that door, as is StateEstimationModel.from_pgm_json(..., asymmetric=True) — the same class, with solve, solve_batch, observability and bad_data unchanged, since none of them knows which domain it is in. On a balanced document the two domains agree per phase to 1e-6, which is a cross-check neither domain could perform alone.

   And Gauss-Newton has a divergence cliff, in both domains. It takes an undamped step, so a measurement far enough from consistent throws the first one past the basin: on transmission-case the phase domain converges at a 15σ error and diverges at 20σ, the symmetric problem converges at 20σ and diverges at 50σ. It reports MaxIterations rather than a converged-looking answer, which is the part that matters — a diverged state carries residuals of 1e17, and a bad-data report built on one is nonsense that looks like a diagnosis. A line search would raise the cliff; there is none, and the behaviour is gated so that adding one will announce itself.

   The symmetric sensors that fixture carries need no conversion beyond a rotation, and it is worth recording why: a voltage reading is line-to-line over u_rated in the scalar case and line-to-neutral over u_rated/√3 here, which is the same number for a balanced set, and a power reading is a three-phase total over s_base against a per-phase value over s_base/3, likewise. So the value replicates and only the angle rotates by 0/−120/+120 — which is exactly power-grid-model's own ComplexValue<asymmetric_t> broadcast.

   Asymmetric sensors describe their three phases separately here rather than reducing to the symmetric problem: asym_voltage_sensor, asym_power_sensor and asym_current_sensor each select their own phase's reading, against a line-to-neutral voltage base and a s_base/3 power base. Checked against single-node-source-asym-voltage-sensor, where the sensor determines the answer outright and power-grid-model reports back exactly what it read.

   What is left is one modelling difference, and it is worth stating precisely. Voltage magnitudes alone do not determine a phase relationship. With flows in the set the source impedance couples the phases — a flow through it depends on all six of its phasors — but given only magnitudes the three per-phase rotations are three separate symmetries where StateLayout removes one, and the gain matrix is correctly singular. power-grid-model answers such a case because its source is a boundary condition, a fixed balanced three-phase voltage; gridoxide's is an unknown behind a synthesized impedance, the same difference that leaves SeReport::unconstrained naming a virtual bus per source on the symmetric side.

   The obvious fix is wrong, and it was tried rather than assumed. gridoxide builds that virtual bus balanced, so constraining its three angles to differ by ±120° looks like free information and removes exactly the two directions in question. It also contradicts the data: power-grid-model's own single-node-source-asym-voltage-sensor reads three phases whose sequence angles are 0.1, 0.2 and 0.3, on a node with no appliance — zero injection, so zero current through the source branch, so V_virtual = V_node exactly. The virtual bus is as unbalanced as the measurement says the node is, and the constraint moves that fixture's answer from 0.1 to 0.2. The balance is a property of the initial state gridoxide synthesizes, not of the equivalent it represents.

   So this is not a missing feature but a real limit: those two directions are undetermined, and reporting singular is the correct answer. power-grid-model answers instead because it has no source-internal bus to be undetermined about. Both directions are asserted in tests/se_three_phase_test.rs, so a change that supplies them will announce itself.

   Nothing is refused now, where four components used to be, and pgm_3ph_maps no longer returns a Result — it could not fail, and saying otherwise in a signature is its own kind of wrong. Every committed fixture builds in both domains.

   The winding-pair restriction went when transformer_seq_params grew power-grid-model's general zero-sequence algorithm. voltage_regulator is now ignored rather than refused: these maps serve state estimation and nothing else, and a regulator pins a generator's bus and lets Q float — a boundary condition for a power flow, meaningless where every bus voltage is an unknown to be recovered. Refusing it made an estimable document inestimable for a reason that could not have changed the estimate, and inconsistently, since the symmetric path read the same document without comment. And link is now modelled: an ideal connection is phase-transparent, so its sequence parameters are one branch_calc_param repeated — which is literally what power-grid-model's own link.hpp computes, calc_param_y_asym(y_link, 0.0, y_link, 0.0, 1.0). The short-circuit path had modelled it that way all along and carried its own copy of the loop; there is one now, shared, so the two cannot disagree about what a link is. Six committed state-estimation fixtures use a link and none of them could be estimated in the phase domain before.

   Two things came with it. A current sensor on a link is refused in the phase domain, as it already was symmetrically and as power-grid-model refuses it — a link's admittance is a chosen regularization constant, so the current through one is an artifact of that choice. And the Kirchhoff gate that checks the terminal functionals against the Y-bus now runs over a link-bearing document, with the link's admittance pinned separately against LINK_Y and against coupling no phase to any other — because Kirchhoff compares two descriptions built from the same blocks and would pass a link that was consistently wrong.

   And three_winding_transformer is modelled, which was the last of the four and the only one that was ever a real gap. It is three two-winding legs to a synthesized star node — the same star equivalent the symmetric path uses, three_winding_legs now assembling it for both so they cannot diverge — plus the star node itself, which is three buses like any other and injects nothing. What the phase domain adds is each leg's winding configuration, which the symmetric model never had to carry: T1 is wye_n/wye_n at clock 0 whatever winding_1 is, while T2 and T3 carry their own winding against winding_1 at a reversed clock, exactly as three_winding_transformer.hpp::convert_to_two_winding_transformers specifies. PgmThreeWindingTransformer grew winding_1/_2/_3 fields to hold them; nothing had read them before.

   Gated against power-grid-model's own asymmetric answer for a three-winding transformer over ten tap positions — the only such reference in the vendored corpus, newly committed — and, for state estimation, against its published magnitudes for the fixture whose three power sensors sit on the transformer's three sides, which also exercises resolving a side to its leg.

   What that reference cannot gate is worth stating, because the obvious reading of it is wrong. Its network is a source and a transformer with no load, so the excitation is balanced and no zero-sequence current flows in it at all — the winding configurations, which are the whole difficulty, are invisible to it. Verified rather than assumed: swapping the two sides of legs T2 and T3 leaves every number in that comparison unchanged. So the zero-sequence structure is asserted separately and directly, and that check does fail under the swap.

2. Current sensors. Done. sym_current_sensor and asym_current_sensor are read in both angle frames, on both calculation methods, checked against power-grid-model's own global-current-sensor and local-current-sensor fixtures — which are identical but for the frame and converge to visibly different states, so the distinction is genuinely exercised rather than nominally supported.

   Stored decomposed into real and imaginary components rather than as a magnitude and an angle, following power-grid-model, and for a decisive reason of gridoxide's own: arg(I) has a branch cut and gridoxide has no phase_mod_2pi anywhere, so a polar residual taken near ±π would silently chase a 2π error. |I| also has an unbounded derivative on an unloaded branch. The variance decomposition reproduces power-grid-model's second-order formula exactly.

   Two rules are enforced that power-grid-model checks only in its Python validation layer, its C++ core accepting and double-counting the mixture: a power sensor and a current sensor may not share a terminal, and two current sensors on one terminal may not disagree about the frame. A current sensor on a link is refused outright — a link's admittance is a regularization constant, so the current through one is an artifact of that choice rather than a measurement.

   One divergence worth recording: power-grid-model refuses to run at all when a global-angle sensor has no voltage angle to reference, raising NotObservableError. gridoxide reports it through ObservabilityReport::global_current_without_angle_reference instead. The state is fully determined there — determined to the wrong reference, since StateLayout pins a bus the sensor contradicts — so calling it unobservable would misname it.

3. Batch state estimation. Done. se::batch::SeBatchSolver (crates/gridoxide-core/src/se/batch.rs) estimates many scenarios over one topology and measurement structure, parallel across cores, each worker amortizing one symbolic factorization — the same shape batch::BatchSolver has for power flow, and exactly what PersistentEstimator's already-written cache-validity condition allows. MeasurementOverride varies values and sigmas and refuses to vary kind or target, mirroring BusOverride's refusal to change bus_type. Exposed as StateEstimationModel.solve_batch(scenarios, threads). Checked against power-grid-model's own sensor-update-* and unbalanced-power-measurements-* batch fixtures, and asserted bit-for-bit identical to a sequential loop at every thread count.

On speed, the iterative-linear method runs 1.6-2.0x behind power-grid-model's across an order of magnitude of problem size (scripts/bench/README.md §7). Measured rather than inferred, that is entirely an iteration-count gap: gridoxide's own iterations are 30-40% cheaper than power-grid-model's and it takes about three times as many, and undamped its map does not converge at all. See docs/src/state_estimation/iterative.md.

A fourth gap surfaced while closing the third and is now closed too: de-energized islands. power-grid-model reports a node in a component containing no source as energized: 0 with a state of exactly zero — topology decides it, and a voltage sensor on such a node is simply ignored. gridoxide had no equivalent, and the consequence was not a wrong answer but no answer: jacobian::mask_untouched pins a column nothing structurally touches, which catches a fully isolated node, but a de-energized node reached by a zero-injection constraint or by its own sensor is touched and undetermined, so the gain matrix came back singular. SeNetwork::energized now carries the same topological verdict solver::PersistentSolver has always applied on the power-flow side (network::connected_components + mark_unreferenced_islands): such a bus contributes no rows and no constraints, and is reported at zero.

Two smaller gaps have closed. A sensor on a three-winding transformer side (measured_terminal_type 6/7/8) used to return MeasurementError::UnsupportedTerminalType; it now maps to the corresponding leg's From terminal, since a three-winding transformer is already resolved into three two-winding branches around a star bus. And the estimator no longer starts flat: se::nr::linear_start carries the network's structural phase shifts, without which Gauss-Newton converges to a different stationary point on any network containing a phase-shifting transformer — reporting success, with an objective nine orders of magnitude worse than the true optimum.

One caveat on all of the above that is about evidence rather than features. Bad-data detection is exercised on deliberate single gross errors — bad_data's own a_gross_error_is_detected_and_identified, and, on a real fixture in the phase domain, bad_data_detection_finds_the_corrupted_measurement, which perturbs one phase of one node by fifteen sigmas and requires the culprit to rank first by a factor of five over the runner-up, with a_clean_phase_domain_estimate_is_not_rejected as its other half. That is a real check and a narrow one.

What is still missing is realistic noise: every benchmark and fixture here estimates from data that is perfectly consistent apart from an injected error, so nothing measures how the detector behaves against a plausible measurement set — several small errors at once, correlated errors, or a gross error on a barely-redundant measurement. Nor is lead 1 comparative: power-grid-model has no bad-data detection to compare against, so "better" here means "present". A noise-generating harness does not exist.

Note on realistic benchmark coverage

gridoxide.matpower (python/gridoxide/matpower.py — the conversion logic moved into the pip package itself; scripts/bench/matpower_to_pgm.py is now only a thin CLI wrapper around it) populates voltage_regulator.q_min/q_max from MATPOWER's gen matrix Qmax/Qmin columns (summed across every active gen at a bus, matching how p_specified is already summed). Confirmed against all 12 real benchmark cases: 11 of them have at least one PV bus whose unconstrained Q genuinely exceeds its nameplate limit (from 4 violations on the smallest case to 166 on case3120sp), and outerloop::ReactiveLimits converges on every one of them, including cases needing dozens of simultaneous PV→PQ switches across several outer iterations. MATPOWER represents "no limit" as literal +-Inf on some real cases (e.g. case9241pegase) — the converter omits the key entirely in that case rather than writing a non-standard Infinity JSON token, matching PGM's own "unset means unbounded" convention exactly.