Introduction

gridoxide is an AC power flow solver written in Rust. It solves the power flow equations for an electrical grid with the Newton-Raphson method, using a sparse Jacobian throughout — assembly, factorization, and solve.

This book covers both the method (what the equations are, how the sparse solve works, what each modeling feature changes about the equation system) and the tool (how to build it, which linear solver backends exist, how CGMES input is mapped onto the internal network model).

Where to start

• Getting Started — build the Rust project, run a solve, or pip install gridoxide and drive it from Python.
• Power Flow — the Newton-Raphson formulation, and the three modeling features that change it: reactive power limits, zero-impedance branches, and multiple islands.
• Sparse Linear Solvers — the five interchangeable linear-solver backends, and a step-by-step walkthrough of the KLU algorithm all of them are measured against.
• CGMES Data Model — reading ENTSO-E RDF/XML grid models, and how individual CIM classes map onto buses, branches, and injections.
• Reference — how gridoxide compares against five other power flow tools, where the benchmark numbers live, and the licensing of every vendored and translated piece of third-party code.

Building and Running

Building

You need the Rust toolchain — see rustup.rs for installation instructions. Then:

cargo build

For an optimized release build:

cargo build --release

A default build needs no C compiler, no system libraries, and no environment variables. Everything beyond that — the klu, pardiso, cgmes, and python features — is opt-in, and each is described on the page that covers it (backends, CGMES input, Python bindings).

Running

cargo run

Or run the built executable directly from the project root:

./target/debug/gridoxide     # debug build
./target/release/gridoxide   # release build

Testing

cargo test

Note that cargo test never compiles the python feature — the feature must not be combined with a plain cargo invocation at all (see Python bindings). Tests for the optional backends live in their own files (tests/block_jacobian_test.rs, tests/klu_jacobian_test.rs, tests/klu_native_jacobian_test.rs, tests/pardiso_jacobian_test.rs) and only run when the matching feature is enabled.

Next steps

To measure rather than just run, see Benchmarking and Profiling for the benchmark harnesses, the perf setup, and where the measured numbers are recorded.

Python Bindings

gridoxide ships as a pip-installable package exposing the solver to Python:

pip install gridoxide

Prebuilt wheels are published for Linux (x86_64), Windows, and macOS (arm64).

import gridoxide

model = gridoxide.PowerFlowModel.from_pgm_json("grid.json", backend="klu_native")
model.solve()
print(model.voltage_mag())  # per-unit magnitude, one entry per node
print(model.voltage_ang())  # angle in radians, one entry per node

Grids are loaded from power-grid-model (PGM) JSON input files. python/README.md in the repository is the package's own PyPI landing page and carries a full worked example of that input format.

API

• PowerFlowModel.from_pgm_json(path, backend="scalar", tol=1e-6, max_iter=20, s_base_va=1e6, freq_hz=50.0) — loads a PGM JSON file and builds the Y-bus admittance matrix.
• model.n_nodes — number of buses, including one virtual slack bus per active source.
• model.solve() — runs Newton-Raphson from a flat/linear-initial-guess start; raises RuntimeError if it doesn't converge within max_iter iterations.
• model.reset() — discards the cached symbolic factorization; call before the next solve() if the topology has changed.
• model.voltage_mag() / model.voltage_ang() — per-bus results in node order.

Reusing factorization across repeated solves

PowerFlowModel wraps the Rust solver::PersistentSolver — it is that API, not a reimplementation of it — so repeated .solve() calls on one model reuse the cached symbolic factorization exactly as described in Backends and Factorization Reuse. Construct one model per topology, then solve as many times as needed:

model = gridoxide.PowerFlowModel.from_pgm_json("grid.json", backend="scalar")
for scenario in scenarios:
    apply_scenario(model, scenario)  # changes p/q values only
    model.solve()
    results.append(model.voltage_mag())

Call model.reset() if the topology itself changes between solves, not just bus values.

This is what lets the benchmark suite run its whole comparison in pure Python (scripts/bench/bench_gridoxide_native.py), timing gridoxide with the same time.perf_counter()-around-a-persistent-solve-object methodology every other tool there already uses (PGM's PowerGridModel, lightsim2grid's GridModel, pandapower's net), rather than shelling out to a compiled Rust binary and parsing its stdout.

Backends available from Python

Two further backends exist in the source tree but are not in the published wheel, since they need extra system dependencies at build time — build from source with the matching Cargo feature:

• "klu" — links vendored SuiteSparse C directly (--features python,klu).
• "pardiso" — Intel oneMKL's PARDISO solver (--features python,pardiso, needs MKLROOT set).

See Backends and Factorization Reuse for what each one actually does and how they compare.

Building input grids

Two helpers produce or convert PGM JSON, so a working grid doesn't require hand-writing one. Both ship in the pip package itself (they used to live only under scripts/bench/) and are installed as console scripts:

• gridoxide.generate_grid — synthetic radial MV/LV distribution grid generator at any scale, pure stdlib, no extra dependencies:

  from gridoxide.generate_grid import generate
  generate(target_nodes=2200, seed=42, out_path="grid.json")  # ~2,600 nodes
  

  or gridoxide-generate-grid grid.json --target-nodes 2200 --seed 42.

• gridoxide.matpower (needs pip install gridoxide[matpower]) — converts a raw MATPOWER .m/.mat case into PGM JSON:

  from gridoxide.matpower import convert
  convert("case14.m", "case14.json")
  

  or gridoxide-matpower case14.m case14.json.

scripts/bench/generate_grid.py and matpower_to_pgm.py are thin CLI wrappers delegating to these same package modules, so the conversion logic only lives in one place.

There is deliberately no pandapower-based converter: it would pull in the full pandapower + power-grid-model-io dependency chain, and gridoxide.matpower already covers the same real-world test-case grids straight from their original MATPOWER sources. If you already have a pandapower.pandapowerNet and pandapower installed, scripts/bench/convert_pandapower_case.py in the main repo is a standalone (not packaged) converter.

How the extension is built

src/python.rs exposes PersistentSolver and PGM-JSON loading as gridoxide._gridoxide, a private compiled extension module built with maturin:

maturin develop --release --features python,klu

It is gated entirely behind the opt-in python Cargo feature and compiled by nothing else, so a plain cargo build/cargo test never touches it — the feature must never be combined with a plain cargo invocation.

This is a mixed Rust/Python maturin project (pyproject.toml's python-source = "python" plus module-name = "gridoxide._gridoxide"): pure-Python code lives in python/gridoxide/ and ships in the same wheel as the compiled extension, re-exported through python/gridoxide/__init__.py so callers only ever write import gridoxide.

python/tests/ holds a pytest suite (scalar/block only — no klu, matching what's published) checked against this project's own committed PGM reference fixtures, run by .github/workflows/python.yml on every push/PR. .github/workflows/pypi.yml builds wheels (Linux/Windows/macOS) plus an sdist and publishes to PyPI via trusted publishing on v* tags. The published wheel deliberately omits the klu backend — LGPL-2.1-or-later vendored SuiteSparse source, plus a C compiler and libclang needed on every target platform. See Provenance and Licensing.

The Power Flow Problem

The powerflow problem is about the calculation of voltage magnitudes and angles for all network nodes. The solution is obtained from a subset of voltages and power injections.

Power System Model

Power systems are modeled as a network of nodes (buses) and branches (lines and transformers). Power sources (generators) and sinks (loads) can be connected to the nodes. Each node in the network is fully described by the following four electrical quantities:

• \(\vert V_k \vert\): voltage magnitude
• \(\theta_k\): voltage phase angle
• \(P_k\): active power
• \(Q_k\): reactive power

There are three types of network nodes: VD (also called the slack or reference bus — types::BusType::Slack in gridoxide's own code), PV and PQ. Depending on the node type, two of the four electrical quantities are specified.

Newton Raphson

The goal is to bring a nonlinear mismatch function \(f\) to zero. The value of the mismatch function depends on a solution vector \(x\):

\[ f(x) = 0 \]

As \(f(x)\) is nonlinear, the equation system is solved iteratively using Newton-Raphson:

\[ x_{i+1} = x_i + \Delta x_i = x_i - \textbf{J}_f(x_i)^{-1} f(x_i) \]

where \(\Delta x\) is the correction of the solution vector and \(\textbf{J}_f\) is the Jacobian matrix.

Instead of computing \(\Delta x_i = - \textbf{J}_f(x_i)^{-1} f(x_i)\), the linear equation set

\[ - \textbf{J}_f(x_i) \Delta x_i = f(x_i) \]

is solved for \(\Delta x_i\).

Iterations are stopped when the mismatch is sufficiently small:

\[ f(x_i) < \epsilon \]

Powerflow Solution

The solution vector \(x\) represents the voltage \(V\) either in polar coordinates

\[ \left [ \begin{array}{c} \delta \ \vert V \vert \end{array} \right ] \]

or rectangular coordinates

\[ \left [ \begin{array}{c} V_{real} \ V_{imag} \end{array} \right ] \]

The mismatch function \(f\) represents the power mismatch

\[ \Delta S = \left [ \begin{array}{c} \Delta P \ \Delta Q \end{array} \right ] \]

or the current mismatch

\[ \Delta I = \left [ \begin{array}{c} \Delta I_{real} \ \Delta I_{imag} \end{array} \right ] \]

This results in four different formulations of the powerflow problem:

• power mismatch function and polar coordinates
• power mismatch function and rectangular coordinates
• current mismatch function and polar coordinates
• current mismatch function and rectangular coordinates

To solve the problem using Newton-Raphson, we need to formulate \(\textbf{J}_f\) and \(f\) for each powerflow problem formulation.

Powerflow with Power Mismatch Function and Polar Coordinates

The injected power at a node \(k\) is given by

\[ S_k = V_k I_k^* \]

The current injection into any node \(k\) is

\[ I_k = \sum_{j=1}^N Y_{kj} V_j \]

Substitution yields

\[ \begin{align*} S_k &= V_k \left ( \sum_{j=1}^N Y_{kj} V_j \right )^* \\ &= V_k \sum_{j=1}^N Y_{kj}^* V_j^* \end{align*} \]

\(G_{kj}\) and \(B_{kj}\) are defined as the real and imaginary part of the admittance matrix element \(Y_{kj}\), so that \(Y_{kj} = G_{kj} + jB_{kj}\). This results in

\[ \begin{align*} S_k &= V_k \sum_{j=1}^N Y_{kj}^* V_j^* \\ &= \vert V_k \vert \angle \theta_k \sum_{j=1}^N (G_{kj} + jB_{kj})^* ( \vert V_j \vert \angle \theta_j)^* \\ &= \vert V_k \vert \angle \theta_k \sum_{j=1}^N (G_{kj} - jB_{kj}) ( \vert V_j \vert \angle - \theta_j) \\ &= \sum_{j=1}^N \left \vert V_k \vert \vert V_j \vert \angle (\theta_k - \theta_j) \right (G_{kj} - jB_{kj}) \\ &= \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( cos(\theta_k - \theta_j) + jsin(\theta_k - \theta_j) \right ) (G_{kj} - jB_{kj}) \end{align*} \]

If we perform the algebraic multiplication of the two terms inside the parentheses, and collect real and imaginary parts, and recall that \(S_k = P_k + jQ_k\), we can split this into two equations: one for the real part, and one for the imaginary part.

\[ \theta_{kj} = \theta_k - \theta_j \\ P_k = \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( G_{kj}cos(\theta_{kj}) + B_{kj} sin(\theta_{kj}) \right ) \\ Q_k = \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( G_{kj}sin(\theta_{kj}) - B_{kj} cos(\theta_{kj}) \right ) \]

These are called the power flow equations.

We consider a power system network having \(N\) buses. We assume one VD bus, \(N_{PV}-1\) PV buses and \(N-N_{PV}\) PQ buses. We assume that the VD bus is numbered bus \(1\), the PV buses are numbered \(2,...,N_{PV}\), and the PQ buses are numbered \(N_{PV}+1,...,N\). We define the vector of unknowns as the composite vector of unknown angles \(\theta\) and voltage magnitudes \(\vert V \vert\):

\[ x = \left[ \begin{array}{c} \theta \\ \vert V \vert \\ \end{array} \right ] = \left[ \begin{array}{c} \theta_2 \\ \theta_{3} \\ \vdots \\ \theta_N \\ \vert V_{N_{PV+1}} \vert \\ \vert V_{N_{PV+2}} \vert \\ \vdots \\ \vert V_N \vert \end{array} \right] \]

The right-hand sides of the powerflow equations for \(P_k\) and \(Q_k\) depend on the elements of the unknown vector \(x\).

Expressing this dependency more explicitly, we rewrite these equations as

\[ \begin{align*} P_k^{spec} = P_k^{calc} (x) \Rightarrow P_k^{calc} (x) - P_k^{spec} &= 0 \quad \quad k = 2,...,N \\ Q_k^{spec} = Q_k^{calc} (x) \Rightarrow Q_k^{calc} (x) - Q_k^{spec} &= 0 \quad \quad k = N_{PV}+1,...,N \end{align*} \]

We define the mismatch \({f} (x)\) as

\[ \begin{align*} f(x) = \left [ \begin{array}{c} f_1(x) \\ \vdots \\ f_{N-1}(x) \\ ------ \\ f_N(x) \\ \vdots \\ f_{2N-N_{PV} -1}(x) \end{array} \right ] = \left [ \begin{array}{c} P_2(x) - P_2 \\ \vdots \\ P_N(x) - P_N \\ --------- \\ Q_{N_{PV}+1}(x) - Q_{N_{PV}+1} \\ \vdots \\ Q_N(x) - Q_N \end{array} \right] = \left [ \begin{array}{c} \Delta P_2 \\ \vdots \\ \Delta P_N \\ ------ \\ \Delta Q_{N_{PV}+1} \\ \vdots \\ \Delta Q_N \end{array} \right ] = 0 \end{align*} \]

That is a system of nonlinear equations. The nonlinearity stems from the fact that \(P_k\) and \(Q_k\) have terms containing products of unknowns and also terms containing trigonometric functions of unknowns.

The Jacobian matrix is obtained by taking all first-order partial derivates of the power mismatch function with respect to the voltage angles \(\theta_k\) and magnitudes \(\vert V_k \vert\):

\[ \theta_{jk} = \theta_j - \theta_k \\ \begin{align*} J_{jk}^{P \theta} &= \frac{\partial P_j (x ) } {\partial \theta_k} = \vert V_j \vert \vert V_k \vert \left ( G_{jk} sin(\theta_{jk}) - B_{jk} cos(\theta_{jk} ) \right ) \\ J_{jj}^{P \theta} &= \frac{\partial P_j(x)}{\partial \theta_j} = -Q_j (x ) - B_{jj} \vert V_j \vert ^{2} \\ J_{jk}^{Q \theta} &= \frac{\partial Q_j(x)}{\partial \theta_k} = - \vert V_j \vert \vert V_k \vert \left ( G_{jk} cos(\theta_{jk}) + B_{jk} sin(\theta_{jk}) \right ) \\ J_{jj}^{Q \theta} &= \frac{\partial Q_j(x)}{\partial \theta_j} = P_j (x ) - G_{jj} \vert V_j \vert ^{2} \\ J_{jk}^{PV} &= \frac{\partial P_j (x ) } {\partial \vert V_k \vert } = \vert V_j \vert \left ( G_{jk} cos(\theta_{jk}) + B_{jk} sin(\theta_{jk}) \right ) \\ J_{jj}^{PV} &= \frac{\partial P_j(x)}{\partial \vert V_j \vert } = \frac{P_j (x )}{\vert V_j \vert} + G_{jj} \vert V_j \vert \\ J_{jk}^{QV} &= \frac{\partial Q_j (x ) } {\partial \vert V_k \vert } = \vert V_j \vert \left ( G_{jk} sin(\theta_{jk}) - B_{jk} cos(\theta_{jk}) \right ) \\ J_{jj}^{QV} &= \frac{\partial Q_j(x)}{\partial \vert V_j \vert } = \frac{Q_j (x )}{\vert V_j \vert} - B_{jj} \vert V_j \vert \\ \end{align*} \]

The linear system of equations that is solved in every Newton iteration can be written in matrix form as follows

\[ \begin{align*} -J(x) \left [ \begin{array}{c} \Delta \theta \\ \Delta \vert V \vert \end{array} \right ] &= \left [ \begin{array}{c} \Delta P \\ \Delta Q \end{array} \right ] \\ \Rightarrow J(x) \left [ \begin{array}{c} \Delta \theta \\ \Delta \vert V \vert \end{array} \right ] &= \left [ \begin{array}{c} -\Delta P \\ -\Delta Q \end{array} \right ] \end{align*} \]

\[ \begin{align*} \left [ \begin{array}{cccccc} \frac{\partial \Delta P_2 }{\partial \theta_2} & \cdots & \frac{\partial \Delta P_2 }{\partial \theta_N} & \frac{\partial \Delta P_2 }{\partial \vert V_{N_{G+1}} \vert} & \cdots & \frac{\partial \Delta P_2 }{\partial \vert V_N \vert} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \Delta P_N }{\partial \theta_2} & \cdots & \frac{\partial \Delta P_N}{\partial \theta_N} & \frac{\partial \Delta P_N}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta P_N}{\partial \vert V_N \vert} \\ \frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_2} & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_N} & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_N \vert} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \Delta Q_N}{\partial \theta_2} & \cdots & \frac{\partial \Delta Q_N}{\partial \theta_N} & \frac{\partial \Delta Q_N}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_N}{\partial \vert V_N \vert} \end{array} \right ] \left [ \begin{array}{c} \Delta \theta_2 \\ \vdots \\ \Delta \theta_N \\ \Delta \vert V_{N_{G+1}} \vert \\ \vdots \\ \Delta \vert V_N \vert \end{array} \right ] = \left [ \begin{array}{c} -\Delta P_2 \\ \vdots \\ -\Delta P_N \\ -\Delta Q_{N_{G+1}} \\ \vdots \\ -\Delta Q_N \end{array} \right ] \end{align*} \]

Solution Steps

1. Set the iteration counter to \(i=1\). Use the initial solution \(V_{i} = 1 \angle 0^{\circ}\)
2. Compute the mismatch vector \(f({x_i})\) using the power flow equations
3. Check the stopping criterion
   • If \(\vert \Delta P_{i} \vert < \epsilon_{P}\) for all type PQ and PV buses and
   • If \(\vert \Delta Q_{i} \vert < \epsilon_{Q}\) for all type PQ
   • Then go to step 6
   • Else, go to step 4
4. Evaluate the Jacobian matrix \(\textbf{J}_f(x_i)\) and compute \(\Delta x_i\).
5. Compute the new solution vector \(x_{i+1}\) and return to step 3.
6. Stop.

Reactive Power Limits (PV → PQ Switching)

Motivation

The plain Newton-Raphson formulation on the Powerflow page treats every PV bus as if its generator could supply or absorb any amount of reactive power needed to hold \(\vert V_k \vert\) at its setpoint. Real generators can't: each has a nameplate reactive capability \(Q_k^{min} \le Q_k \le Q_k^{max}\), usually tightest for round-rotor machines operating near their active-power limit. If the voltage setpoint on a PV bus would require \(Q_k\) outside that range, the unconstrained solution is not physically achievable — the generator saturates at its limit and the bus can no longer hold its voltage.

The standard fix is PV → PQ switching: solve normally, check every PV bus's computed \(Q_k\) against its limits, and for any bus that violates one, convert it to a PQ bus with \(Q_k\) pinned at the violated limit — trading the voltage-magnitude equation for a reactive-power equation at that bus — then re-solve. Repeat until no bus violates its limits.

What changes in the equation system

Recall from the Powerflow page that a PV bus contributes only a \(P\)-mismatch row to \(f(x)\) — its voltage magnitude \(\vert V_k \vert\) is known (the setpoint), so it isn't part of the unknown vector \(x\) and there's no \(Q\)-mismatch row for it. A PQ bus contributes both a \(P\)- and a \(Q\)-mismatch row, and its \(\vert V_k \vert\) is part of \(x\).

Switching bus \(k\) from PV to PQ therefore:

• adds \(\vert V_k \vert\) to the unknown vector \(x\) (it's no longer fixed at the setpoint — the voltage is now free to float),
• adds a \(Q\)-mismatch row \(Q_k^{calc}(x) - Q_k^{spec}\) to \(f(x)\), with \(Q_k^{spec}\) pinned at whichever limit was violated (\(Q_k^{min}\) or \(Q_k^{max}\)) instead of the bus's original \(q\_spec\),
• leaves that bus's \(P\)-mismatch row untouched — the same equation, unaffected by bus type.

This changes \(n_{unknowns}\) itself, so the Jacobian's dimensions grow by one row and one column per switched bus. Any cached symbolic factorization (fill-reducing ordering) computed for the old sparsity pattern is invalid once a switch happens and must be redone from scratch.

Switching strategies

That equation-system change is common to every implementation. What differs is when the check runs and whether a switch can ever be undone — three genuinely different designs, in increasing order of bookkeeping.

1. One-directional outer-loop switching

Solve to full convergence, check limits, switch every violating bus to PQ, re-solve; repeat until a pass switches nothing. A bus that has been switched to PQ stays PQ for the rest of the call, even if a later pass's solution would put its \(Q\) back within limits.

The appeal is that the anti-oscillation problem simply doesn't arise: with switching one-way, the set of PQ buses grows monotonically, so the outer loop terminates on its own. The cost is that a bus which was only transiently out of range — say during an early pass, still far from the real solution — is committed to PQ permanently, giving a slightly different, slightly more conservative answer than a scheme that could release it.

2. Bidirectional switching with a switch-count cap

The same outer loop, except a bus moved to PQ may move back to PV in a later pass if its computed \(Q\) has returned within limits — recovering exactly the transiently-violating case strategy 1 gives up on. This reintroduces the risk strategy 1 avoids: a bus can flip between PV and PQ indefinitely across passes, since each type change moves the solution enough to justify the opposite change next time. The standard mitigation is a hard cap on how many times any one bus may switch, after which it is forced to stay wherever it last was.

That cap is real bookkeeping — a per-bus counter carried across outer passes — and it makes the result mildly path-dependent (which bus hits its cap first depends on iteration order). In exchange, no bus is written off on the strength of one early, poorly-conditioned iterate.

3. Mid-iteration switching

Rather than a separate outer loop that fully re-solves after every switch, the check-and-switch step is folded directly into the Newton iteration: once the mismatch drops below some threshold, bus types are switched and the same iteration sequence continues instead of a fresh solve starting. This saves the redundant re-convergence strategies 1 and 2 pay for — the iterate is already near the solution when the switch happens, so it doesn't have to be rediscovered — at the cost of mutating the equation system underneath a running Newton iteration, so the Jacobian's dimensions (and any cached factorization) change mid-flight.

Published versions of this technique also widen the switch-back criterion: "On PV-PQ Bus Type Switching Logic in Power Flow Computation" (Jinquan Zhao) bases the decision to release a bus on comparing its voltage against its setpoint in addition to comparing its \(Q\) against its limits, not on \(Q\) alone.

Where this fits in gridoxide today

solver::newton_raphson_enforcing_q_limits() implements strategy 1, as an outer loop around a solver::PersistentSolver:

1. Solve to convergence with every PV bus free (PersistentSolver::solve). If this doesn't converge, stop and return that status — there's no point checking limits on a non-converged solution.
2. Compute every bus's actual \(Q_k\) (network::power_injections).
3. For each PV bus whose \(Q_k\) violates q_min/q_max: switch bus_type to PQ and pin q_spec to the violated limit.
4. If no bus was switched this pass, the solution is self-consistent — stop and return Converged.
5. Otherwise, reset the cached factorization (PersistentSolver::reset) and go back to step 1.
6. If max_outer_iter outer passes pass without stabilizing, stop and return MaxIterationsReached.

Strategies 2 and 3 both add real complexity (anti-oscillation bookkeeping, or reworking the iteration loop itself) to recover a case gridoxide's simpler design just accepts as a one-way commitment — a deliberate scope trade-off, not an oversight.

One further simplification on top of the one-directional rule: bus voltages carry over between outer passes rather than resetting to a flat start. Since only one bus's type changes per pass, the previous pass's converged state is normally very close to the next equilibrium, so re-solving after a switch typically takes only a few extra Newton iterations, not a full fresh solve.

Scope: net injection, not per-device limits

Bus::q_min/q_max bound the bus's net reactive injection — the same aggregate quantity q_spec already represents when multiple loads/generators share a node (see pgm::PgmVoltageRegulator's own q_min/q_max parsing, which sums each active generator's limits at a bus exactly the way p_spec is already summed). Pinning q_spec to a violated limit only achieves that limit exactly if the bus carries no voltage-dependent ZIP-model load terms (Bus::zip_terms) — true for the common case of a bus that's purely a PV/generator connection, not true in general if a voltage-dependent load is co-located there.

Validated against real data

tests/q_limits_test.rs checks the mechanics directly: a 3-bus fixture with one PV bus whose q_min is set tight enough to force a violation, confirmed to switch to PQ, pin q_spec at exactly q_min, and converge to the same voltages across every Jacobian backend (Scalar, Block, Klu) — as well as a control case with a loose q_min that shouldn't trigger any switch at all, converging to the same answer as plain unconstrained Newton-Raphson.

Beyond the unit fixture, this has been exercised against all 12 real MATPOWER benchmark cases (gridoxide.matpower — python/gridoxide/matpower.py, with scripts/bench/matpower_to_pgm.py as a thin CLI wrapper around it — populates voltage_regulator.q_min/q_max from each case's own gen matrix): 11 of the 12 have at least one PV bus whose unconstrained \(Q\) genuinely exceeds its nameplate limit — from 4 violations on the smallest case up to 166 simultaneous violations on case3120sp.

Tool reference

Ideal Switches and Zero-Impedance Branches

Motivation

Real networks contain connections that are, by design, not really "branches" in the sense the Powerflow page assumes: breakers, disconnectors, bus-bar couplers, jumpers — elements meant to either tie two points electrically together with (idealized) zero impedance, or fully separate them, with no impedance value in between and no partial state.

Modeling one as an ordinary branch runs straight into the admittance formulation the rest of that page relies on. A branch's contribution to the bus admittance matrix is built from \(Y = 1/Z\); with \(Z = 0\) that's an infinite entry, not a large-but-finite one. And this isn't just a numerical inconvenience to work around — a closed zero-impedance connection isn't naturally a "branch" at all. Its physical meaning is an equality constraint, \(V_i = V_j\), not a current flowing in proportion to a voltage difference. That's a fundamentally different kind of equation from every other term in \(f(x)\), which is exactly why real tools solve this with dedicated mechanisms rather than a variant of the ordinary branch stamp.

Three approaches

1. Topological reduction (merge before you formulate)

The most direct fix is to never let the zero-impedance connection reach the admittance formulation in the first place. Build a graph of nodes connected by closed zero-impedance edges, find its connected components, and treat each component as a single node in the actual power-flow model — an open switch is simply not an edge, so it doesn't merge anything.

This adds no new equations, no new unknowns, and no numerical stiffness. The costs are that it happens before the solver ever runs, as a separate graph pass over the input model, and that the two original terminals lose their distinct identity — nothing downstream can distinguish "flow through the left side of the coupler" from "flow through the right side" once they've become one node. Switching state also can't change between solves without redoing the reduction and rebuilding the model.

2. Large-admittance regularization

A second approach keeps the zero-impedance connection as an ordinary branch, but assigns it a series admittance many orders of magnitude larger than anything else in the network (say 10⁶ S in a 10 kV network) and zero shunt admittance. The same nodal equations that already exist for every other branch then force near-equality of the two bus voltages as a natural numerical consequence — no new equation type, no graph pass, just an extreme parameter value. Open/closed status needs no special handling either: it's the same connection-status flag every other branch already has.

The trade-off is conditioning: introducing one admittance value five or six orders of magnitude larger than the rest of the matrix widens its dynamic range substantially, which is exactly the kind of thing a direct sparse solver's pivoting has to work harder to stay accurate through as a network grows. That's a real cost, paid deliberately in exchange for needing no separate topological pre-processing stage at all — the ideal connection is just a branch, as far as the rest of the solver is concerned.

3. Equality-constrained augmented system

The third approach takes the constraint interpretation from the Motivation section literally: instead of merging the two nodes away or approximating the constraint with an extreme admittance, add the equation \(V_i = V_j\) (and, in polar form, \(\theta_i = \theta_j\)) directly into the Newton-Raphson system as its own row. Since adding an equation without adding an unknown would leave the system over-determined, a pair of new "dummy" variables is introduced alongside it — one at each of the two buses, equal and opposite, representing the (otherwise unmodeled) power flowing through the ideal connection to make the constraint hold. The system stays exactly square: one new equation, one new pair of unknowns.

A complication this approach has to handle and the other two don't: a loop of zero-impedance branches produces a linearly dependent constraint set (the last edge's constraint is implied by the others), which would make the augmented system singular. The standard fix is to compute a minimum spanning tree over each connected group of zero-impedance branches and constrain only the tree edges, leaving the redundant non-tree edges as inactive constraints that can be reactivated if the tree has to be re-routed after a later switching operation.

Compared to the other two, this is real bookkeeping — the spanning-tree maintenance in particular has no analogue in either simpler approach — but it's also the only one of the three that keeps both original terminals numerically distinct after the fact, which matters if anything downstream needs to report or reason about flow through each side of the connection separately.

Where this fits in gridoxide today

gridoxide's solver core has no node-breaker layer: every Bus reaching network::build_ybus is already a fully-resolved electrical node. Topology resolution is a precondition of the Bus/Line/Transformer model, so the choice above is made once, at import.

The rule: identity, not impedance

Merge only when the element has no identity in the output model. Otherwise it stays a branch.

The deciding evidence is what the fixtures assert. All four upstream power-flow cases containing a link publish that link's own current, and vision-validation-network publishes its full p/q/s at both ends. Two state-estimation fixtures publish per-node injections either side of a link. Merging deletes the branch those numbers describe, so it is unavailable wherever they are asserted — branch satisfies all eight, merge satisfies two.

Note this is decided by what the element is, not by what is attached to it. An earlier version of this rule asked whether the endpoints carried appliances, and vision-validation-network's link 74 refutes it: the link joins two nodes with no appliances at all — apparently the safest possible merge — and power-grid-model still reports 1.964 MW flowing through it.

The threshold row exists because a connection can be electrically zero-impedance without being declared as one. ill-conditioned-by-line-meshed carries a line at 7.07e-9 p.u.; only a value-based test catches that, which is powsybl's framing rather than power-grid-model's.

The numbers, and why they were measured

Both constants were chosen by sweeping, not derived, because the two calculation types pull in opposite directions:

• Power flow wants a link stiff. The drop across it is \(\Delta V = I/y\), and the fixtures check node voltages and the link current at 1e-5 relative.
• State estimation wants it soft. \(G = H^{T}WH\) squares the admittance, so power-grid-model's 1e8 becomes 1e16 in the gain matrix.

The window is about one order of magnitude wide. That narrowness is the argument for treating these as regularization parameters with measured values rather than physical constants: a network far outside these fixtures' power scale may need them re-measured, and if no value serves both, approach 3 is the exit — it imposes \(V_i = V_j\) exactly, with no large number anywhere.

topology::ZERO_IMPEDANCE_THRESHOLD is 1e-7 p.u. on the same basis: measured across all 86 branches in the committed PGM fixtures and every CGMES one, it sits above the single pathological line at 7.07e-9 and below the smallest legitimate branch (1.0e-6 in PGM, 2.92e-6 in CGMES), so it disturbs nothing currently modelled.

Detection and treatment are separate numbers. The threshold only decides whether a branch is an ideal connection; a branch it catches is clamped to IDEAL_CONNECTION_Y, the same admittance a declared link gets. One number cannot do both jobs: it has to sit below every legitimate branch, yet clamping merely to that level leaves \(|Y|\) as high as 1e7 — inside the range measured as singular for state estimation, and 35x stiffer than a link. Separating them satisfies both, and says the right thing besides: a line that short is an undeclared link, so it should be treated as one.

powsybl's own threshold is 1e-8, but its default treatment is the equality-constrained formulation rather than a clamp, so it never has to reconcile the two roles in a single value.

Consequences of the merge, where it is used

CGMES inherits approach 1's limitations directly: switching state cannot change between solves without re-importing, and no per-side flow is reportable across a merged switch. There is also an empirical reason that side merges rather than stamping branches — it was tried, and the AC Newton-Raphson solve diverged on FullGrid with 20-odd such branches active at once, which is exactly the conditioning cost approach 2 carries.

Tool reference

powsybl uses both approaches 1 and 3, at different layers and for different reasons: closed switches that its bus/branch view merges away never reach the solver at all, while retained switches that survive into a node/breaker solve get the augmented-system treatment.

Multi-Island Power Flow

Motivation

Real networks aren't always one connected system. Switched-out feeder sections, decommissioned equipment, or genuinely separate synchronous areas can leave a Y-bus with more than one disconnected component. Before this feature, newton_raphson solved the whole bus list as a single system regardless: if any disconnected component had no slack/reference bus in it, that component's Jacobian rows were structurally singular (no reference angle to anchor them against), and the entire solve failed — even when every other component was perfectly solvable on its own.

The concepts

Supporting disconnected networks means answering four separate questions. They're largely independent — a tool's answer to one doesn't determine its answer to the next.

1. What counts as "one island"?

The obvious definition is a connected component of the network graph. But there are two defensible graph definitions once DC links exist:

• Connected component — traverse AC branches and DC/HVDC links.
• Synchronous component — traverse AC branches only.

Two AC areas linked only by an HVDC tie are one connected component but two separate synchronous components. That's physically the right distinction, because a DC link carries power without providing angle/frequency synchronization between its two sides: each side needs its own angle reference, so each side is its own power-flow problem.

A tool with no DC modeling at all doesn't need the distinction — the two notions coincide for it — but the seam is worth keeping conceptually, since retrofitting it after the fact is harder than carrying a second, currently-redundant notion of "component".

2. Which components actually get solved?

Two answers in practice:

• Main component only — find the largest component, solve it, drop the rest. Cheap and matches the common case where everything but the main grid is modeling debris, but the dropped buses silently produce no result at all.
• Every component — solve each one independently. Every bus in the input gets an answer, at the cost of doing work on components a user might not care about.

Note that this choice is separable from the classification itself: the classification (what the components are) can run unconditionally while only the selection of which to solve is configurable.

3. What happens to a component with no reference bus?

A component containing no slack/source bus has no anchor for voltage angle — it's not a solvable power-flow problem, no matter how well-formed the rest of the network is. Three distinct responses:

• Fail the whole solve — structurally correct in that the system really is singular, but it destroys the results for every other component too. This is the behavior the motivation above describes as the actual bug.
• Emit a null/zero placeholder — mark the component's buses with fixed \(V = 0,\ P = Q = 0\) so they contribute nothing to the shared system, and report them as unsolved. There's no principled way to fabricate a reference phasor for a genuinely sourceless island, so none is invented.
• Drop the component entirely — the natural consequence of "main component only"; the buses just aren't in any output.

The mirror-image case is a component with more than one reference bus. That's physically over-determined — two independently fixed slack phasors in one electrically connected component — but it is not necessarily numerically singular, so Newton-Raphson will happily converge to a result satisfying neither slack's true power balance. A mismatch-based convergence check can't detect this after the fact; it has to be decided up front from the topology.

4. How is the result reported?

A single overall status (Converged/Diverged) can't describe a network where one island converged and another didn't. Per-island reporting means one status per component, so a caller can tell "the part I care about converged" from "everything failed".

There's a real limitation to be honest about here: if every component is solved in one shared sparse factorization (as opposed to genuinely separate solves), a singularity detected on the result vector can't be attributed to a specific component — the factorization is one object. Per-island status after a shared solve is therefore best-effort for that particular failure mode, and exact only for the ones decided from topology up front.

Where this fits in gridoxide today

gridoxide solves every solvable component in one call — no mode to choose, no opt-in; this is simply what run_power_flow_analysis/run_power_flow_analysis_from_ybus do. Sourceless components get the zero-voltage placeholder, and each component gets its own status.

The key architectural fact that makes this cheap rather than a solver rewrite: network::YBusSparse already exposes everything needed (n(), row(i)) to build a connectivity graph with zero new plumbing, and every Jacobian backend already excludes BusType::Slack buses from its unknowns with no assumption that there's exactly one Slack bus anywhere in the whole bus list. So marking a sourceless component's buses as fixed Slack placeholders makes the existing newton_raphson/PersistentSolver correctly solve every other component in the same shared call — mathematically equivalent, iteration for iteration, to solving each disconnected component independently, since there's no Jacobian coupling between them. Zero changes to the solver algorithm itself were needed; this is a data-preparation pass before the existing solve and a status-reporting pass after it.

gridoxide has no DC-side traversal in this pass today, so its notion of "component" is both notions at once (§1) — src/dc.rs's HVDC converters are resolved into AC-side injections before build_ybus runs, so a DC link is not an edge in the graph this pass walks, making its components synchronous components by construction.

The pipeline

The partitioning pipeline lives inside solver::PersistentSolver::solve / solver::newton_raphson_with_backend / solver::newton_raphson themselves — not in run_power_flow_analysis_from_ybus. That means every production entry point shares one mechanism: run_power_flow_analysis_from_ybus is a thin convenience wrapper (build a Y-bus, run linear_initial_guess, construct a one-shot PersistentSolver, call .solve()), and a caller who reaches for PersistentSolver, newton_raphson_with_backend, or newton_raphson directly — to reuse a cached factorization, pick a specific backend, or set a custom tol/max_iter — gets the exact same island handling for free, with no separate opt-in and no risk of silently bypassing it. All three return Vec<solver::IslandReport> (previously SolveStatus/()).

batch::BatchSolver inherits it too, since each of its workers calls PersistentSolver::solve per scenario — batching changes how many solves run and on which thread, not what a solve does.

One deliberate exception: bde::solve_batch_block_diagonal does not partition islands. It stacks every scenario's Jacobian into one block-diagonal sparse matrix to validate the GPU path's architecture (see scripts/bench/README.md §4d), and is an architecture validator rather than a production entry point — its own module documentation says so and points CPU callers at BatchSolver. Worth stating here rather than only there, because "every entry point partitions islands" is exactly the kind of invariant that gets relied on later. If BDE ever becomes a production path, island partitioning has to be added to it: a sourceless component would otherwise make its whole block singular, and because the blocks are structurally disjoint (§4 below), that block alone — not the batch — would be the part that fails.

Each of those functions does, internally:

1. network::connected_components — an iterative DFS over YBusSparse's own adjacency, generic over the finished Y-bus (so it works uniformly for native JSON, PGM-JSON, and CGMES input with no format-specific code).
2. network::classify — counts each component's existing Slack buses:
   • exactly one → normal, solvable.
   • zero → network::mark_unreferenced_islands pins every member bus to a fixed V = 0, P = Q = 0 placeholder. There is no principled way to fabricate a reference voltage/angle for a genuinely sourceless island, so none is attempted — every unit/sign-convention bug this project has actually fixed got fixed by matching verified physical or reference-implementation behavior, never by guessing, and inventing a slack here would repeat that mistake class.
   • more than one → left untouched in the shared solve, and reported as AmbiguousReferenceBus. This verdict is decided once, up front, and is never later overwritten by a numerically-convergent-looking mismatch check — for the reason §3 gives.
3. One shared solve (via whichever JacobianBackend the caller picked) across the whole (now correctly-classified) bus set.
4. solver::finish_island_reports recovers each component's own status after the shared solve, restricting the same power_injections/effective_injection mismatch computation newton_raphson_cached's own convergence check already uses to just that component's bus indices.

newton_raphson_enforcing_q_limits composes on top the same way it always did: each of its outer PV→PQ-switching passes calls PersistentSolver::solve (now returning Vec<IslandReport>) and returns that same Vec<IslandReport> — from its last pass once Q-limits have stabilized, or immediately if some island's own status is Singular/ MaxIterationsReached.

PowerFlowReport and IslandStatus

#![allow(unused)]
fn main() {
pub struct PowerFlowReport {
    pub buses: Vec<Bus>,
    pub islands: Vec<solver::IslandReport>,
    /// Iteration count and per-iteration convergence trace. The solve loop
    /// itself prints nothing — see `solver::SolveStats`.
    pub stats: solver::SolveStats,
}

pub enum IslandStatus {
    Converged,
    MaxIterationsReached,
    Singular,
    NoReferenceBus,        // no Slack bus at all — placeholder values only, never solved
    AmbiguousReferenceBus, // more than one Slack bus — solved, but not necessarily meaningful
}
}

A real limitation, stated plainly rather than papered over

Singular cannot, in general, be attributed to a specific island — the shared-solve caveat from §4, concretely. Every backend solves one combined sparse factorization, and singularity is detected via a finiteness check on the result vector (sparse.rs), not a specific failing pivot or row. When the overall solve reports Singular, every component whose own post-hoc mismatch is still above tolerance gets marked Singular too — best-effort, not precise. A component whose own mismatch is already below tolerance is marked Converged regardless of the overall status, since every backend returns Singular before applying that iteration's update (confirmed directly in both the Scalar and Block backends), so buses always holds the last fully-applied, fully-finite state — never a partially-updated or NaN-poisoned one — making this post-hoc check meaningful even after an overall Singular result.

One existing "sourceless placeholder" precedent

cgmes.rs's de-energized-bus block — uses CGMES's own TopologicalIsland membership (a semantic guarantee from the standard: "only energised TopologicalNode-s shall be part of the topological island"), which encodes real domain knowledge the raw admittance graph doesn't have — not just a bypass-caller accommodation, unlike the PGM case above. Kept unconditionally: per the FullGrid diagnosis in Motivation, if the generic pass replaced this block, FullGrid's real converter gap (EquivalentBranch unimplemented) would get silently repackaged as 4 spurious NoReferenceBus islands instead of surfacing as the connectivity bug it actually is.

The generic pass layers on top of CGMES's own block, unconditionally, for every entry point — a harmless, idempotent no-op for input that's already de-energization-resolved by CGMES's own TopologicalIsland logic, and the mechanism that gives CGMES richer per-island reporting once a file genuinely declares more than one TopologicalIsland.

Validated

tests/multi_island_test.rs covers, entirely through the public API (no internal-only test module needed, since PersistentSolver::solve is itself public and does the full partitioning/reporting internally): two well-formed islands matching independent single-island solves exactly; a no-slack island correctly placeholder-reported without poisoning a well-formed neighbor (the literal repro of the original bug); an ambiguous (two-slack) island likewise not poisoning a neighbor; all three non-trivial statuses coexisting in one call, in connected_components's bus-index-ascending order; and an easy-converging island correctly reported Converged even while a deliberately hard one drags the overall solve to MaxIterationsReached — exercised by calling PersistentSolver::solve directly with a custom low max_iter, since run_power_flow_analysis_from_ybus's own tol/max_iter are fixed.

Tool reference

Neither tool auto-detects a mode from the network's actual shape: powsybl's DEFAULT_COMPONENT_MODE is a fixed constant, and PGM has no mode at all. The classification itself always runs unconditionally; only selection is ever configurable.

The State Estimation Problem

Motivation

Power flow answers a planning question: given the injections at every bus, what are the voltages? Every input is known exactly, and there is one right answer.

Operations does not work that way. What a control room has is a few hundred telemetered readings, each with its own error, arriving from equipment of varying quality — some quantities measured twice, many not measured at all, and a handful simply wrong because a transducer has drifted or a sign convention was mis-wired years ago. No state reproduces all of them, because they contradict each other.

State estimation asks the question that fits that data: which grid state best explains these measurements? It is the entry point to essentially every operational application, because those applications need a complete, consistent snapshot of the network and telemetry never provides one.

Weighted least squares

Let \(x\) be the state — every bus's voltage magnitude and angle — and \(z\) the vector of measurements. Each measurement has a measurement function \(h_i(x)\): what that sensor would read if the grid were in state \(x\). A voltage sensor's is trivial (\(h = \vert V_k \vert\)); a branch power sensor's is the terminal flow; a bus injection's is the familiar power-flow expression.

The estimate is the state minimizing the weighted sum of squared disagreements:

\[ J(x) = \frac{1}{2} \left( z - h(x) \right)^{T} W \left( z - h(x) \right) \]

with \(W = \Sigma^{-1}\) diagonal, \(W_{ii} = 1/\sigma_i^2\). Weighting by inverse variance is what makes the estimate maximum likelihood under independent Gaussian errors, and it is why a sensor's declared \(\sigma\) matters as much as its value: a reading trusted to 0.1% pulls a hundred times harder than one trusted to 1%.

Setting \(\partial J / \partial x = 0\) and linearizing gives the normal equations, solved repeatedly until the step vanishes:

\[ G , \Delta x = H^{T} W r, \qquad G = H^{T} W H, \qquad r = z - h(x) \]

where \(H = \partial h / \partial x\) is the measurement Jacobian. This is Gauss-Newton.

How this differs from the power-flow Newton loop

The two loops look alike and are not the same.

The residual does not go to zero. solver::newton_raphson drives a mismatch to zero: an exact solution exists and Newton converges onto it. Gauss-Newton has no such target — the measurements are inconsistent by construction, so \(r\) stays nonzero at the optimum. Convergence is therefore tested on the size of the step, never on the residual. A converged estimate with a large \(J(x)\) has not failed to converge; it means the data disagrees with itself, which is what bad-data analysis is for.

There are no PV buses. A generator's voltage is a quantity to be estimated, not asserted. So every bus contributes a magnitude unknown, and every bus but the angle reference contributes an angle:

\[ x = \left[ \theta_0 \ldots \theta_{N-1} \text{ except the reference}, ; V_0 \ldots V_{N-1} \right] \]

giving \(n = 2N - 1\), against power flow's \(n_{angle} + n_{pq}\). Confusing the two layouts is the easiest way to produce a Jacobian that is subtly and consistently wrong, so se::jacobian::StateLayout owns the mapping and nothing else indexes by hand.

The angle reference is conditional. A network measured only in magnitudes and powers is invariant under a global phase shift, so a reference must be pinned or \(G\) is singular however many measurements there are. But that invariance disappears the moment a phasor measurement supplies an absolute angle — and pinning a reference on top of one is a false constraint that rotates the whole estimate away from the data. gridoxide pins a reference exactly when no angle is measured.

Why the normal equations, and what they cost

\(H\) is rectangular (\(m \times n\)), and every sparse backend gridoxide has is square — they are power flow's Jacobian solvers. Forming \(G = H^{T} W H\) makes the system square and symmetric, so Scalar, Block, Klu, KluNative and Pardiso all carry state estimation with no new solver abstraction. That is the single most important reuse in the design.

The price is conditioning: \(G\) squares \(H\)'s condition number. Two consequences follow, and both shape the implementation:

• Zero injections are enforced as hard constraints rather than as very-high-weight pseudo-measurements, so the weights never span more orders of magnitude than the physics requires.
• An orthogonal (QR) formulation, which never forms \(G\), remains the escape hatch if a real network ever proves too ill-conditioned for this path.

Validation

gridoxide is checked against power-grid-model's own state-estimation fixtures, committed under tests/data/pgm/state_estimation/ with their MPL-2.0 license files. On transmission-case — 11 buses, 4 transformers, 59 measurements — the per-unit magnitudes agree to 1.5 × 10⁻⁹.

Angles are checked in two regimes, because only one of them has an absolute answer. Where an angle is measured, absolute angles must match. Where none is, the estimate is invariant under a global rotation and power-grid-model's own fixtures do not agree with each other on the convention (transmission-case reports its source node at exactly 0, 1os2msr-no-angle reports its source node at −0.0130). The test then requires gridoxide's angles to match up to one constant shared by every bus — a stronger check than it sounds, since a wrong estimate produces per-bus errors rather than a uniform offset.

Two methods

gridoxide implements both of power-grid-model's calculation methods, and they agree with each other bus by bus on its fixtures.

The linearized method is faster per iteration and needs more of them. It is power-grid-model's default; gridoxide's default is Newton-Raphson, on the grounds that a library should be exact unless asked otherwise.

Running it

From the shell:

gridoxide estimate tests/data/pgm/state_estimation/1os2msr/input.json

From Python:

import gridoxide

model = gridoxide.StateEstimationModel.from_pgm_json("grid_with_sensors.json")
model.solve()
print(model.voltage_mag())          # per-unit, one entry per bus
print(model.observability())        # what the measurements fail to determine
print(model.bad_data())             # chi-squared, p-value, ranked suspects

Note the input needs sensors, but does not need p_specified on its loads or u_ref on its sources — those are quantities state estimation solves for rather than inputs it consumes.

Measurements and What They Mean

Everything on this page is about turning sensors into the rows of \(z\) and \(H\). It is the part of state estimation with the least mathematics and the most opportunity for a silent, confident, wrong answer — a flipped sign or a misattributed quantity does not fail, it converges to something plausible.

One measurement is one scalar

A Measurement is a single scalar observation, not a sensor. Newton-Raphson WLS treats \(\sigma_P\) and \(\sigma_Q\) as independent, so a power sensor becomes two rows; a voltage sensor becomes one or two depending on whether it carries an angle. A voltage sensor that does carry an angle is a phasor (PMU) measurement, and it is what makes the global phase determinate.

Aggregation runs in two opposite directions

Redundant sensors are merged before estimation rather than passed through as extra rows, and the two cases move opposite ways:

Several sensors observing one quantity merge by inverse variance:

\[ z = \frac{\sum_k z_k / \sigma_k^2}{\sum_k 1 / \sigma_k^2}, \qquad \sigma^2 = \frac{1}{\sum_k 1 / \sigma_k^2} \]

The result is more certain than any input. A quantity you measure twice is one you know better.

Several appliances making up one bus injection sum, and their variances sum with them:

\[ S = \sum_k S_k, \qquad \sigma^2 = \sum_k \sigma_k^2 \]

The result is less certain than its parts, because independent errors accumulate.

Conflating these is easy and costly. An early version of gridoxide's aggregation merged per bus rather than per appliance, which turned two sensors watching one load into two loads and doubled the injection. The ordering that works is: merge per appliance first, then sum across appliances at the bus.

Sign conventions are not uniform

Taken from power-grid-model's reference-direction rules, which differ by sensor type:

Branch measurements pass through unchanged, since that is already the convention branch_flow::terminal_flow uses. Loads and shunts are negated to become injections.

Where gridoxide's model differs from power-grid-model's

This is the subtlety that most affects correctness, and it is not a sign issue — it is a difference in what the quantity is.

power-grid-model treats a source and a shunt as appliances at a node, so their power counts toward that node's injection. gridoxide models both structurally: a source becomes a virtual slack bus feeding through an impedance branch, and a shunt becomes a Y-bus diagonal entry. Neither therefore appears in network::power_injections at that bus at all — by Kirchhoff's current law, the net injection at a source node with no load is zero, because the source's power arrives through a branch that is part of the network.

Each of the three needs its own measurement function:

Using the plain bus injection for any of them is wrong by construction. This was found rather than predicted: tests/measurement_residual_test.rs evaluates every measurement function at the state power-grid-model published and reported a 63-sigma disagreement on exactly that quantity, with the model saying 0 and the sensor saying 2.4 p.u.

Zero-injection buses

A bus with no load and no generator injects exactly nothing. That is a property of the network, not an observation of it: no sensor, no noise, no uncertainty.

The common shortcut is to feed it in as a pseudo-measurement of zero with a very small \(\sigma\). It works, and it is why so many estimators are described as ill-conditioned — the weight matrix then spans many orders of magnitude and \(G\) squares that spread. power-grid-model ships fixtures named ill-conditioned-by-line-meshed and ill-conditioned-by-link-meshed for precisely this failure.

gridoxide enforces it as a hard equality constraint instead, via the Lagrangian stationarity conditions:

\[ \begin{bmatrix} G & C^{T} \\ C & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \lambda \end{bmatrix} = \begin{bmatrix} H^{T} W r \\ -c(x) \end{bmatrix} \]

The augmented matrix is symmetric indefinite rather than positive definite, which rules out a Cholesky-style solver — but every gridoxide backend is a general sparse LU, so this costs assembly work only. A constraint row turns out to be the same bus-injection partials an injection measurement would produce; the difference between a constraint and a measurement is entirely in how the system consumes the row.

Two details worth knowing:

• Which buses qualify is read off the input document, not off Bus::p_spec. A state-estimation document leaves p_specified unset, so an unmeasured load looks exactly like zero injection in the converted network while being nothing of the sort. Getting this backwards would constrain real loads to zero.
• Sources and shunts do not disqualify a bus, because gridoxide models both structurally and neither appears in that bus's injection. The virtual slack buses are excluded, since that is where a source's unknown power enters.

The fixture that pins this down is node-injection-sensor-and-zero-injection: an injection sensor reading 0.1 p.u. on a node with no appliance attached. power-grid-model requires at least one appliance for such a sensor to mean anything, so it overrides the reading and reports both buses at exactly 1.0∠0. Without the constraint, weighted least squares fits the sensor perfectly by driving that node to \(\sqrt{2} \angle -45°\) — a 41% overvoltage on a bus with nothing connected to it, and an objective of \(4 \times 10^{-28}\). By its own criterion, a flawless answer.

The Iterative-Linear Method

power-grid-model's default calculation method, and gridoxide's optional one. It trades exactness for a much cheaper iteration.

The trick

In terms of the complex voltage vector \(\underline{U}\), the awkward part of the measurement model is that power is bilinear in voltage: \(S = U \overline{I}\). Current is not — \(I = YU\) is perfectly linear. So if the measurements were currents rather than powers, the whole estimation problem would be a single linear least-squares solve.

So they are converted. A power measurement \(S\) at a terminal whose voltage is \(U\) becomes a current measurement

\[ \underline{I} = \overline{\left( S / U \right)} \]

using the voltage from the previous iterate. A magnitude-only voltage measurement becomes a phasor by borrowing the previous iterate's angle. Each pass redoes the conversion with better voltages, and the linearization error shrinks as the angles settle.

Zero-injection buses are markedly simpler here than in the nonlinear formulation: a bus with no appliance carries no injected current, so \((YU)_i = 0\) exactly, with no linearization involved. They enter the same augmented KKT system in complex arithmetic.

Where the speed comes from

The system matrix is built from admittances and measurement weights. Admittances are fixed. The weights are nearly fixed — converting \(S\) to \(I\) divides by \(U\), so a current's variance scales as \(1/\vert U \vert^2\) — and the method assumes \(\vert U \vert\) constant for that purpose.

Granting that assumption, the matrix never changes. It is factorized once, and each iteration moves only the right-hand side. Newton-Raphson, by contrast, rebuilds and refactorizes its Jacobian every pass.

What it costs

The constant-\(\vert U \vert\) assumption, plus the interpretation of power measurements as current measurements, means this method's optimum is not exactly the weighted-least-squares optimum. power-grid-model documents the same caveat and recommends its Newton-Raphson method when precision matters. gridoxide keeps both, and defaults to Newton-Raphson on the grounds that a library should be exact unless asked otherwise.

There is also a capability difference worth knowing: a lone \(P\) or \(Q\) measurement, without its partner on the same target, is dropped. Half a complex power cannot be converted into a current. The Newton-Raphson path handles the two independently and keeps it.

Under-relaxation, and why it is safe here

The linearization inverts the voltage, and that can oscillate: an estimate that is too low produces a current that is too high, which pushes the next estimate too high, and back. The iteration then sits in a two-cycle indefinitely rather than converging.

This is not hypothetical. power-grid-model's transmission-case fixture does exactly that — the step parks at \(5.5 \times 10^{-2}\) and stays there for hundreds of iterations, neither converging nor diverging:

iteration:  1     2     3     5     10    20    50    100
step:       6.9e-2 5.2e-2 5.6e-2 5.6e-2 5.6e-2 5.5e-2 5.4e-2 5.3e-2

gridoxide backs the step off — halving a relaxation factor whenever an iteration fails to improve on the last, and letting it grow again by 1.2 whenever progress resumes — which breaks the cycle:

iteration:  1     2     3     5     10    20
step:       6.9e-2 5.2e-2 2.8e-2 6.7e-5 2.8e-6 5.9e-9

Damping is a suspicious remedy in general, because it can make a solver look convergent while walking toward the wrong answer. It is safe here for a specific reason: it changes the path and never the fixed point. iterative::tests::the_true_state_is_a_fixed_point pins that down directly — seeded at a state satisfying the measurements, one iteration leaves the residuals at zero and moves the state only by the global rotation the estimator normalizes away. So a damped run converges to the same state an undamped one would reach if it got there at all.

Relaxation engages whenever an iteration fails to improve on the last by 10%, which is not the same thing as "only on a pathological problem" — a flat start alone can trip it, and on case1354pegase it does, at iteration 3. That is exactly why the factor has to be able to recover: see Letting the relaxation recover below, where a version that could only descend turned out to be spending whole runs at half a step because of one transient.

Measured cost

examples/bench_se.rs runs both methods over the MATPOWER benchmark grids, with measurements synthesized from a converged power flow — a voltage magnitude at every bus and flows at both ends of every branch, giving a redundancy around 4x.

Both the benefit and the cost show up plainly. Newton-Raphson takes exactly six iterations at every scale — quadratic convergence — while the linearized method takes 16 to 45. It still wins on wall clock because each of its iterations is so much cheaper, and the margin widens with size, from 5.5x at 119 buses to 11.4x at 2,870: prefactorization amortizes better the larger the matrix.

The accuracy trade is equally visible. Against the state the measurements were read from, Newton-Raphson lands at ~1e-14 throughout — the data is perfectly consistent here, so it is limited only by arithmetic. The linearized method runs 8.7e-10, 1.6e-8, 8.9e-8, 5.6e-7, 1.1e-6 down that same column: still far below any real measurement noise, but degrading with size rather than holding constant. That is the linearization and the constant-|U| weighting showing through, and it is why Newton-Raphson remains the default.

Where the remaining time goes

Profiling was worth doing before optimizing, because it redirected the obvious plan. On case1354pegase the two methods spend their time quite differently:

The tempting optimization is a symmetric factorization: G is symmetric and a general sparse LU is being used on it, so in principle half of Newton-Raphson's 78% is recoverable. Two measurements argue against it.

First, zero-injection buses are 9-31% of buses on these grids, so the constrained KKT system is the normal case, not the exception, and it is symmetric indefinite — needing Bunch-Kaufman or a quasi-definite regularization rather than a Cholesky, which is two code paths and a perturbation of the constraints.

Second and more decisively, it would speed up the wrong method. Newton-Raphson is where gridoxide is already ahead — power-grid-model's own Newton-Raphson state estimator raises SparseMatrixError on every case from 300 buses up, where gridoxide's converges. The remaining gap is on this method — 1.5 ms against 0.7 ms on case300, 7.2 ms against 3.4 ms on case1354pegase, 18.3 ms against 8.4 ms on case2869pegase — and there the factorization is not the bottleneck at all. §7 of scripts/bench/README.md has the full table.

So the open lead for this method is its convergence rate, not its linear algebra: it needs 45 iterations here where Newton-Raphson needs 6.

Measured against power-grid-model's own iteration count

That last paragraph was a conclusion about gridoxide on its own. Comparing the same method in the two tools sharpens it into something more useful, and reverses the reading §7 of scripts/bench/README.md used to give.

Neither tool reports an iteration count through its public API, so scripts/bench/se_iterations.py obtains it the same way for both: the smallest max_iterations budget that does not fail, found by bisection. Before comparing anything it is worth knowing the two criteria are the same quantity, and they are — power-grid-model's iterate_unknown returns max over buses of |u_new − u_old| phase-normalized and loops while (max_dev > err_tol); gridoxide's raw_step is .map(|(a, b)| (a - b).norm()).fold(0.0, f64::max) over the same normalized voltages, and both default to 1e-8. The one asymmetry is that gridoxide tests raw_step × relaxation, where power-grid-model has no relaxation at all and always takes the full step.

On the documents examples/bench_se.rs --emit writes:

(gridoxide's counts here are the ones that stood before Letting the relaxation recover, which cuts them by about 40%.)

gridoxide's iterations are individually cheaper than power-grid-model's — by 30-40% on every case above 100 buses — and it takes about three times as many of them. The flat ~2x total is the product of those two effects pulling against each other, not evidence of a constant-factor gap in the per-iteration work. Reading a stable ratio as a per-iteration constant was the error; a stable ratio is equally consistent with two ratios that happen to be stable.

The control that makes this comparison mean anything is that both tools reach the same answer: max |Δu| between their solutions is 7.6e-9 to 5.2e-7 across these cases, i.e. agreement at their shared tolerance. Same problem, same optimum, different paths to it.

So the lever is convergence rate, confirmed rather than inferred — and specifically not the linear algebra, which is already ahead.

Letting the relaxation recover

The section above establishes two things that pull in opposite directions: the damping cannot be removed, because the undamped map does not converge at all; and the damping is what makes gridoxide take three times power-grid-model's iterations. Both are true, and the way between them is that the factor could only ever go down.

The rule was: halve whenever an iteration fails to improve on the last by 10%. Nothing ever restored it. And the flat-start transient alone trips that test — at iteration 3 on case1354pegase, before the iteration has settled into anything — so a run would spend its whole length at half a step because of one early stumble that had nothing to do with the instability the damping exists for.

Adding one clause fixes it: when the step is shrinking, by more than 25%, grow the factor by 1.2 again, capped at 1. The result on the benchmark documents:

About 40% fewer iterations, and 16-22% less wall clock — measured interleaved, old build and new, two rounds. The two do not match because a fixed setup-and-factorization cost does not shrink with the iteration count; that it is now a larger share of the total is the point of having cut the rest. The iteration-count gap to power-grid-model falls from 2.0-3.3x to 1.2-2.0x, and on case118 gridoxide is now the faster of the two outright.

The growth rate is measured rather than derived, and the measurement is the interesting part. Forcing a constant relaxation shows this map's optimum sits near 0.7 — 23 iterations on case300 against 35 at 0.5 — with the map going unstable just above it: 0.8 costs 36 iterations and 0.9 costs 80. So there is a narrow good range whose location depends on the network, which is an argument for hunting for it adaptively rather than for hardcoding 0.7. Among growth factors that all converge, 1.2 has the best worst case: 1.3 costs case2869pegase 30 iterations and 1.4 costs case14 seventy, where 1.2 needs 20 and 19.

The answer is untouched, as it must be — the_true_state_is_a_fixed_point already pins that damping changes the path and never the destination, and the estimates agree to 2.5e-9 across the change, which is below the 1e-8 the iteration is asked for.

This does not explain why power-grid-model needs no damping at all. That remains open, and the next section narrows it.

The relaxation is load-bearing, not overhead

The obvious next move from that table is to take the damping out, since power-grid-model manages without it. That does not work, and finding out why moves the problem somewhere more interesting.

Forced undamped on these same documents, the iteration does not converge slowly — it does not converge at all. The step falls for three iterations and then locks onto a constant:

iter        1        2        3        4       ...      58       59       60
step   6.09e-1  2.27e-1  1.77e-1  1.76e-1     ...  1.748610e-1  1.748610e-1  1.748610e-1

Constant to seven digits for fifty-odd iterations is a limit cycle, not a floor — the step is not shrinking at all, so the iterate is circulating rather than settling. That under-relaxation at 0.5 restores convergence places the dominant eigenvalue near −1: damping maps λ to 0.5 + 0.5λ, which sends λ ≈ −1 to ≈ 0. power-grid-model's map, on the identical document with the identical measurements and no damping at all, is stable there.

The rows responsible are the branch-terminal power ones. Dropping them from the undamped run drops the step from 1.7e-1 to 7.1e-4 — three orders of magnitude — while dropping the voltage rows instead changes nothing at all. That fits their shape: a branch row converts its reading with I = conj(S/U_at), so its right-hand side depends on the inverse of a voltage the same row solves for, and on these documents the branch rows outweigh the voltage rows by four orders of magnitude.

Two further candidate mechanisms were tested and both ruled out:

• The |U|² weight scaling. gridoxide scales a power row's weight by the reference bus's starting |U|², where power-grid-model deliberately does not — its iterative_linear_se_solver.hpp says so directly ("the variance is not scaled as an approximation"). Removing the scaling moves the limit cycle's amplitude (1.75e-1 → 1.81e-1 on case300) and does not remove it.
• The zero-injection KKT constraints. These are 66 to 869 buses on the cases above, and power-grid-model has no equivalent augmentation. Dropping them entirely shrinks the amplitude to 5.7e-2 and still leaves a limit cycle.

Whatever destabilizes the map is therefore in its core — the power-to-current conversion, the phase normalization, or the interaction between them — rather than in either of the two places gridoxide visibly departs from power-grid-model. That is the open question, and it is worth answering before anyone reaches for an accelerator: a scheme layered on an unstable map inherits the instability.

Note this is a different measurement set from the one the next section traces. These documents carry voltage magnitudes and branch flows only; examples/se_converge.rs synthesizes bus injections too, and its undamped run decays rather than cycling. Both are real, and the difference between them is itself a clue about which rows drive the instability.

What actually ends the iteration

Tracing that convergence rate is worth doing before trying to improve it, because it is not one phenomenon but two, and the second one is not a rate at all.

examples/se_converge.rs re-runs the estimate at increasing max_iter and prints the step, its ratio to the previous step, and the error against the state the measurements were read from:

cargo run --release --example se_converge case1354pegase

Forced undamped, case1354pegase decays geometrically at a strikingly stable ratio — 0.757, holding to three digits from iteration 30 through iteration 53. A single dominant mode, textbook material for Aitken or Anderson acceleration. Then it stops:

iter        50       51       52       53      ...      197      198      199      200
step    6.3e-6   4.7e-6   3.6e-6   2.7e-6      ...   2.5e-7   3.4e-7   2.9e-7   7.3e-8

Out to 200 iterations the raw step never reaches the 1e-8 tolerance; it settles on a floor around 1e-7 and bounces. Forcing any constant relaxation between 0.4 and 1.0 gives the same picture.

What ends the default run is the relaxation ratchet:

iter      1        2        3     ...      39       40       41       42       43       44       45
raw    1.19e0   1.66e0   1.72e0   ...   3.2e-7   3.7e-7   2.1e-7   2.5e-7   2.6e-7   4.7e-7   1.9e-7
relax    1.0      1.0      0.5    ...      0.5     0.25     0.25    0.125   0.0625  0.03125  0.03125

Relaxation engages at iteration 3, in the flat-start transient, and holds at 0.5 through iteration 39. Then the raw step stops falling, the 10%-improvement test fails four times over the next five iterations, and the factor halves to 1/32. Convergence is declared at iteration 45 on a reported step of \(1.86 \times 10^{-7} \times 0.03125 = 5.8 \times 10^{-9}\).

That reported step is honest about what it measures — the state genuinely moved that little, because that is how far damping let it move. But it is not the iteration reaching 1e-8. The tolerance is met by the damping factor.

The answer is unaffected, and this is the part worth being clear about. The error against the true state bottoms out at 5.5e-7 by iteration 36 — the same 5.6e-7 this method scores on case1354pegase in the table above, i.e. its own linearization bias, the floor no amount of iterating can go below. Iterations 37 to 45 buy no accuracy in any case.

Two consequences for whoever picks this up next:

• Accelerating the 0.757 mode wins iterations 3 through 39 and nothing after. The floor is a separate phenomenon and would still be there.
• A step floor near 1e-7 on a gain matrix whose weights span 1e4 to 1e6, with exact KKT constraints augmented in, is the shape of linear-solve accuracy rather than of the fixed-point map — worth testing with iterative refinement or a scaling pass before assuming the iteration is at fault.

An earlier version of this section reported roughly eighty iterations, with relaxation holding at 1.0 until iteration 78. Neither figure reproduces at HEAD; the traces above are what the current code does.

Agreement with Newton-Raphson

On every fixture gridoxide checks, the two methods agree to \(10^{-6}\) per-unit, bus by bus — a stricter statement than each agreeing with power-grid-model, since it holds at every bus rather than only where an expected value was published.

Selecting it:

gridoxide estimate grid.json --iterative-linear

model = gridoxide.StateEstimationModel.from_pgm_json(
    "grid.json", method="iterative_linear", max_iter=100
)

The larger iteration budget is deliberate: this method converges linearly where Newton-Raphson converges quadratically. Spending more iterations to make each one far cheaper is the entire point.

Observability and Bad Data

A converged estimate is not a correct one. Weighted least squares will absorb a broken sensor by bending the state toward it, and will happily return an answer for a network where half the quantities were never determined at all. Both failures look exactly like success. These two analyses are what distinguish them.

Observability: what the measurements determine

An unobservable system is not a numerical accident, and treating it as one produces the least useful diagnostic available — a factorization that failed, with no indication why. The useful answer names the unknowns nobody is watching.

There are two distinct failures.

Structural. A column of \(H\) that is identically zero: no measurement function mentions that unknown at all. The estimator has to find these anyway, since they make \(G\) singular on their own. In gridoxide they are usually the virtual slack buses synthesized per source — a network power-grid-model considers fully observable can still leave gridoxide with surplus unknowns, because PGM has no such bus in its state space.

Numerical. A column that is present but linearly dependent on the others. Two branch-flow measurements on a radial feeder with no injection measurement between them constrain the same combination of angles twice and leave another combination free. Nothing about the sparsity pattern gives this away; it takes a rank computation.

\(G = H^{T} W H\) is symmetric positive semidefinite by construction, and the standard rank-revealing factorization for that class is a Cholesky with symmetric diagonal pivoting: take the largest remaining diagonal at each step, stop when it falls below tolerance, and the number of steps is the rank. gridoxide uses faer's blocked implementation for the factorization but makes the rank decision itself — faer stops at \(\varepsilon n\), around \(10^{-16}\) relative, which is the threshold for "not numerically positive definite". Observability wants a much looser one: a direction determined only at the \(10^{-12}\) level is not usefully determined, and calling it observable would be the more misleading answer.

The cost is that \(G\) is densified, making this \(O(n^3)\) time and \(O(n^2)\) memory. The analysis refuses above DENSE_LIMIT and says so, rather than quietly allocating gigabytes on a transmission grid; a sparse rank-revealing method is the natural follow-up.

Unknowns pinned this way are held at their starting value rather than estimated, and reported. That is deliberate: an unknown nothing observes cannot be recovered, and moving it would be inventing information. It also improves the failure mode — an under-measured system now solves its observable part and names the rest, instead of failing outright.

Bad data: whether the measurements agree

Two questions, answered separately.

Is there bad data at all?

Under the assumption that each error is independent and normal with the declared \(\sigma\), the objective is chi-squared distributed:

\[ J = r^{T} W r \sim \chi^2(m - n + k) \]

with \(m\) measurements, \(n\) estimated state variables and \(k\) equality constraints — each constraint gives a degree of freedom back, and a pinned unobservable variable never consumed one. A \(J\) far out in the tail says the residuals are too large to be noise.

Rejection says something is wrong, not which measurement. A test that does not reject is also not a clean bill of health: a single moderate error, or several that partly cancel, can sit inside the threshold.

Which measurement?

The largest normalized residual. Raw residuals are not comparable — 0.1 is enormous against a \(\sigma\) of 0.001 and negligible against a \(\sigma\) of 1 — and dividing by \(\sigma\) alone is still not enough, because a redundantly measured quantity spreads its error across its neighbours and so under-shows in its own residual. The right scale is the residual's own standard deviation:

\[ r_i^{N} = \frac{\vert r_i \vert}{\sqrt{\Omega_{ii}}}, \qquad \Omega = R - H G^{-1} H^{T} \]

conventionally compared against 3. \(\Omega\) is computed through the augmented system, so the zero-injection constraints count: a constrained estimate has less freedom to move, which changes how much of an error surfaces in the residual rather than in the state.

\(\Omega_{ii}\) costs one linear solve per measurement, so the full diagonal costs more than the estimate itself. gridoxide shortlists candidates by the cheap proxy \(\vert r_i \vert / \sigma_i\) and computes \(\Omega_{ii}\) exactly for the worst 20, configurable. This is an approximation with a real failure mode, not a free shortcut: because \(\Omega_{ii}\) varies per measurement, the true worst can in principle sit outside the shortlist. Raising the limit to \(m\) removes the doubt at the corresponding cost.

A measurement whose residual has no variance at all is critical in the usual terminology — the estimate is forced to reproduce it exactly, so an error in it cannot be detected by any amount of analysis. Those are skipped rather than assigned a meaningless normalized residual.

What the fixtures show

Run over power-grid-model's own state-estimation fixtures:

The last one's worst suspect is its injection sensor at a normalized residual of exactly 100.00 — the 100-sigma conflict that fixture is built around, recovered as a number.

One caveat when reading those figures: the consistent fixtures produce \(\chi^2 \approx 0\) rather than \(\chi^2 \approx \text{dof}\), because power-grid-model generated their readings from the true state without adding noise. There is nothing for a correct estimate to disagree with. Real telemetry would sit near its degrees of freedom, and a near-zero statistic on real data would itself be suspicious — it would suggest the declared sigmas are far too large.

Backends and Factorization Reuse

The Y-bus admittance matrix and Newton-Raphson Jacobian are stored and factored as sparse matrices, not dense ones. This page covers why that matters, how factorization work is reused across solves, and what the five interchangeable linear-solver backends are.

Why sparse

At 2,605 nodes the original dense solver was over 100,000x slower than power-grid-model; the sparse rewrite closed that to roughly an order of magnitude on a cold, single-shot solve. The underlying grid is sparse — each bus only connects to a handful of neighbors — so a dense representation was doing asymptotically unnecessary work regardless of how fast the constant-factor arithmetic was.

Two things make this work, rather than just "swap in a sparse matrix type":

• Sparse-aware assembly. Both the Jacobian build (solver::build_jacobian_triplets) and the linear initial-guess warm start (network::linear_initial_guess) walk each bus's actual admittance neighbors (via network::YBusSparse::row) instead of looping over every possible bus pair. The O(n²)/O(m²) assembly cost has to go too, or a sparse solve alone doesn't fix the bottleneck.
• Symbolic factorization reuse. A Newton-Raphson Jacobian has the same sparsity pattern every iteration — same bus topology, only numeric values change — so solver::newton_raphson computes the symbolic factorization (fill-reducing ordering) once and reuses it for a cheap numeric-only refactorization on each iteration (sparse::RealSparseSystem), mirroring what PGM's own solver does internally.

Inside KLU walks through exactly what "symbolic factorization" and "numeric-only refactorization" mean, step by step, on a real Jacobian.

Reusing factorization across repeated solves

A single newton_raphson/newton_raphson_with_backend call reuses its symbolic factorization across its own NR iterations, but starts cold on every call — re-deriving the fill-reducing ordering from scratch. That is fine for a genuinely one-off solve and wasteful for anything that solves the same topology repeatedly: a time series, a batch of scenarios, contingency analysis. In those, only bus values (p_spec, q_spec, voltage guess) change between calls, not the topology.

solver::PersistentSolver extends the reuse across calls. Construct one per topology, then solve as many times as needed; only the first call pays for symbolic factorization.

#![allow(unused)]
fn main() {
use gridoxide::solver::{JacobianBackend, PersistentSolver};

let mut solver = PersistentSolver::new(JacobianBackend::Klu);
for scenario in scenarios {
    apply_scenario(&mut buses, scenario); // changes p_spec/q_spec only
    solver.solve(&mut buses, &ybus, 1e-6, 20);
}
}

Call .reset() (or construct a new solver) if the topology itself changes between solves. This is a meaningful win on real-world grids, since a cold solve otherwise redoes COLAMD/AMD/BTF ordering from scratch every call. examples/bench_network.rs exposes it as an optional warm mode; its default cold mode still measures "N independent flat-start solves with no shared state," a different and also legitimate number. See Benchmarking for the measured warm-vs-cold figures.

The backend interface

src/sparse.rs is the thin backend wrapper around faer, and intentionally the only file that imports faer types directly, so a different sparse-solver backend can be swapped in behind the same interface without touching the rest of the codebase.

solver::newton_raphson always uses the default Scalar (faer-backed) path. newton_raphson_with_backend additionally accepts four alternatives via solver::JacobianBackend:

Block

src/block_sparse.rs, no extra build requirements. Groups each bus's own (angle, voltage-magnitude) unknowns into one dense 2×2 block, mirroring power-grid-model's block-per-bus matrix structure, with a hand-written Gilbert-Peierls sparse block LU (block_sparse::BlockLu). Symmetric power flow only. Consistently faster than Scalar.

Klu

src/sparse_klu.rs, needs cargo build --features klu. The same scalar Jacobian as Scalar, solved by SuiteSparse's KLU instead of faer, vendored and compiled from source (vendor/suitesparse/) rather than depending on a third-party Rust wrapper crate. Needs a C compiler and libclang (for bindgen) at build time.

KLU and BTF (one of KLU's own dependencies) are LGPL-2.1-or-later. See Provenance and Licensing.

KluNative

src/klu_native/, no extra build requirements. A from-scratch Rust translation of the same KLU algorithm Klu links over FFI: BTF block-triangular preprocessing, per-block AMD ordering, a partial-pivoting Gilbert-Peierls LU kernel with Eisenstat-Liu pruning, and cheap numeric-only refactorization — all faithfully ported, not a simplified reimplementation. No C compiler or libclang needed, so unlike Klu it is always built.

Validated end-to-end against real KLU on all 13 real MATPOWER benchmark cases (14 to 9,241 buses): identical iteration counts and identical converged voltages on every case. One known, documented gap — row scaling (klu_native::scale) is ported and independently tested but not yet wired into the factor/refactor path (see src/klu_native/mod.rs's module doc comment). That is a numerical-stability preconditioning step, not a correctness one, so it doesn't affect the results above.

Inside KLU is a walkthrough of this port specifically.

Pardiso

src/sparse_pardiso.rs, needs cargo build --features pardiso and MKLROOT set at build time. The same scalar Jacobian as Scalar, solved by Intel oneMKL's PARDISO direct solver. Unlike Klu, nothing is vendored — MKL is proprietary (Intel Simplified Software License, not OSS), so this only dynamically links a locally-installed oneMKL (libmkl_rt.so, discovered via MKLROOT, e.g. source /opt/intel/oneapi/setvars.sh) and generates FFI bindings via bindgen against that install's own mkl_pardiso.h. No MKL header or source is copied into this repo.

PARDISO's C API is one function called repeatedly with different phase values against a persistent opaque handle, rather than KLU's separate analyze/factor/refactor/solve functions. mtype = 11 (real, nonsymmetric) and iparm[34] = 1 (0-based indexing) are the two settings that matter for matching gridoxide's CSR/CSC conventions.

Not built or tested in CI — no CI runner has MKL installed — so this is a local/manual-verification-only backend.

How they compare

All four alternatives are strictly parallel to Scalar, not replacements. A bug in any of them can't affect newton_raphson's default behavior, and every existing test keeps using Scalar unless it explicitly opts into a different backend.

All five backends produce identical converged voltages at every scale — these are purely performance comparisons, not correctness trade-offs. In rough terms:

• Block, Klu, and KluNative are all meaningfully faster than Scalar.
• Klu and KluNative land close to each other, slightly ahead of Block.
• Pardiso carries a largely size-independent fixed setup cost from its default matching/scaling preprocessing, making it the slowest backend at small problem sizes — even behind Scalar — though it scales better than Scalar as node count grows.

PGM is clearly faster than any gridoxide backend on synthetic radial-distribution/LV topology, a real, standing gap this project hasn't closed. That gap doesn't hold universally, though: on real-world transmission-topology grids, gridoxide's Klu backend is frequently faster than lightsim2grid's own KLU-backed C++ solver. The comparison depends on topology, not just implementation language.

Benchmarking points at the full measured numbers, exact ratios, and how to reproduce them.

Inside KLU: the Sparse Solve, Step by Step

Every Newton-Raphson iteration ends in the same place: solve \(\textbf{J} \Delta x = f(x)\) for the correction \(\Delta x\). For anything bigger than a toy network, that linear solve dominates the runtime, and the Jacobian is sparse — a bus only couples to its immediate neighbours, so a 22×22 Jacobian from an 11-node grid has 124 nonzeros out of 484 possible entries (74 % empty).

KLU is the sparse LU solver gridoxide uses for that step. It was designed by Tim Davis and Ekanathan Palamadai specifically for circuit-simulation matrices — very sparse, unsymmetric, with a strong zero-free diagonal — which describes a power-flow Jacobian almost exactly. gridoxide has two KLU backends (JacobianBackend::Klu, the vendored C via FFI, and JacobianBackend::KluNative, a from-scratch Rust port of the same algorithm in src/klu_native/), and both run the pipeline described below.

This page walks that pipeline one phase at a time, on matrices small enough to check by hand. Every number shown was produced by running the actual src/klu_native/ code on the matrix shown.

The pipeline at a glance

KLU splits the work into a symbolic phase that depends only on where the nonzeros are, and a numeric phase that depends on their values:

The split is the whole point. Phases 1–2 are the expensive combinatorial work, and in a power flow the Jacobian's pattern never changes between Newton iterations — only its numbers do. So gridoxide runs phases 0–2 once per topology (cached in PersistentSolver) and only phases 5 and 4 per iteration. Phase 5 exists precisely because it can skip everything phase 3 does symbolically.

Step 0: triplets to CSC

jacobian::JacobianPattern derives the Jacobian's sparsity pattern once per topology and hands the backend (row, col, value) triplets in whatever order its assembly walk happens to produce them — once, at LinearSolver::new. (solver::jacobian_triplets_reference/build_jacobian_triplets still assemble the same thing the naive way, but only under #[cfg(test)], as the oracle JacobianPattern is checked against bit-for-bit.) KLU wants compressed sparse column (CSC): one array of column start offsets, one of row indices sorted within each column, one of values.

Take this 4×4 matrix, which we'll reuse for the next several steps:

       c0  c1  c2  c3
r0      .   1   .   3
r1      .   .   2   .
r2      .   4   2   1
r3      5   .   1   .

As triplets and then as CSC (build_csc_structure in klu_native/mod.rs):

col_ptr = [0, 1, 3, 6, 8]     // column j occupies row_idx[col_ptr[j] .. col_ptr[j+1]]
row_idx = [3,  0, 2,  1, 2, 3,  0, 2]
values  = [5,  1, 4,  2, 2, 1,  3, 1]
          |c0|  c1  |    c2    |  c3 |

build_csc_structure also merges duplicate (row, col) pairs and returns a groups mapping from each CSC slot back to the triplets that fed it. That mapping is what makes cheap refactorization possible later: on the next Newton iteration, the same assembly order maps to the same CSC slots, so new values can be packed in without re-deriving the structure. It is also why later iterations hand over only the values — JacobianPattern refills one reused Vec<f64> positionally matching that first triplet order, and LinearSolver::factor_and_solve_values packs it straight into the cached slots, with no (row, col) pair rebuilt after the first iteration.

Step 1: BTF — finding block triangular structure

If a matrix can be permuted to block upper-triangular form, you never have to factor it as one piece. You factor each diagonal block independently and stitch the results together with back substitution. Blocks are smaller, so there's less fill and less work; and a 1×1 block is just a division.

BTF gets there in two moves.

1a. Maximum transversal: get a zero-free diagonal

First find a matching — a set of nonzeros, no two sharing a row or column, as large as possible. This is bipartite matching between rows and columns, solved by augmenting paths (btf/maxtrans.rs, a port of Duff's MC21 algorithm).

For the matrix above, the matching found is:

MAXTRANS: match = [1, 2, 3, 0]   nmatch = 4

Read it as "row i is matched to column match[i]": row 0 ↔ column 1 (value 1), row 1 ↔ column 2 (value 2), row 2 ↔ column 3 (value 1), row 3 ↔ column 0 (value 5). All four rows are matched (nmatch == n), so the matrix has full structural rank and can be permuted to have no zeros on its diagonal.

If nmatch < n, the matrix is structurally singular — no permutation can fill the diagonal — and btf_order completes the permutation arbitrarily so later phases still have a bijection to work with, marking the fake entries with the flip sentinel (klu_native/types.rs). The factorization then fails on a zero pivot, which is exactly how gridoxide reports SolveStatus::Singular for a disconnected island with no reference bus.

1b. Strongly connected components: find the blocks

With the matching in hand, build a directed graph: one node per column, and an edge \(j \to k\) whenever column \(j\) has a nonzero in the row that is matched to column \(k\). Its strongly connected components are exactly the irreducible diagonal blocks, and a topological order of the components gives the block ordering. btf/strongcomp.rs is Tarjan's algorithm, iterative rather than recursive.

For our matrix, the edges are:

c0 → c0                  (row 3 ↔ c0: the matched diagonal entry)
c1 → c3                  (row 2 ↔ c3)
c2 → c3, c2 → c0         (row 2 ↔ c3, row 3 ↔ c0)
c3 → c1                  (row 0 ↔ c1)

c1 → c3 → c1 is a cycle, so {c1, c3} is one component. {c0} and {c2} are singletons. And because c2 has edges out to both other components and none coming back, it must be ordered last. That is what the code returns:

BTF: p = [0, 2, 3, 1]   q = [1, 3, 0, 2]   r = [0, 2, 3, 4]

p and q are the row and column permutations; r gives the block boundaries — block 0 spans positions 0..2, block 1 is position 2, block 2 is position 3. Applying them:

                q0=c1 q1=c3 | q2=c0 | q3=c2
   p0=r0   [     1     3    |   .   |   .    ]
   p1=r2   [     4     1    |   .   |   2    ]
           [ --------------  -------  ------ ]
   p2=r3   [     .     .    |   5   |   1    ]
           [ --------------  -------  ------ ]
   p3=r1   [     .     .    |   .   |   2    ]

Block upper-triangular, with a 2×2 block and two 1×1 blocks. Instead of one 4×4 LU, KLU will do one 2×2 LU and two divisions. The two entries above the diagonal blocks (the 2 and the 1 in the last column) are stored separately, as the off-diagonal part — factor.rs keeps them in their own CSC arrays (off_p/off_i/off_x) because they're not part of any block's LU; they only appear during the solve.

What BTF does on a real power-flow Jacobian: nothing

Worth saying plainly, because it's the common case:

3-bus network.json:        n=3    nnz=9     btf_blocks=1   sizes=[3]
PGM transmission-case:     n=22   nnz=124   btf_blocks=1   sizes=[22]

A connected AC network has an irreducible Jacobian — every bus reaches every other bus through the network graph, so the whole matrix is one strongly connected component and BTF returns a single block. BTF earns its keep on circuit matrices with genuine one-way structure (a driving stage feeding a following stage that doesn't feed back).

The one power-flow case where it does something: a network with several islands that each have a reference bus. PersistentSolver::solve classifies islands but still hands the whole bus list to a single Newton solve (see Multi-Island Power Flow), so the Jacobian is block diagonal — no unknown in one island appears in any equation of another. The column digraph of step 1b then has no edges crossing island boundaries, so no strongly connected component can span two islands, and BTF necessarily recovers one block (or more) per island and factors them independently. Islands with no reference bus never reach this point: mark_unreferenced_islands converts their buses to slack first, removing them from the unknown set entirely.

BTF costs one near-linear pass, so KLU runs it unconditionally rather than trying to predict which case it's in.

Step 2: AMD — reordering to limit fill-in

Inside each diagonal block, the ordering still matters enormously, because of fill-in: entries that are zero in \(A\) but nonzero in \(L\) or \(U\). Fill is what makes a sparse factorization degenerate into a dense one.

The canonical demonstration is a star (an arrow matrix) — one hub node connected to everything, which in power system terms is a substation busbar with many feeders:

       c0  c1  c2  c3  c4
r0      5   1   1   1   1
r1      1   2   .   .   .
r2      1   .   2   .   .
r3      1   .   .   2   .
r4      1   .   .   .   2

Eliminate the hub (node 0) first and every pair of leaves becomes coupled — the remaining 4×4 block fills in completely. Eliminate the leaves first and nothing fills at all, because no two leaves are adjacent. AMD finds the second order:

AMD star perm = [4, 3, 2, 1, 0]         // hub eliminated last

star, natural order  [0,1,2,3,4]:  nnz(L)=10  nnz(U)=10  total LU nnz = 25   (fully dense)
star, AMD order      [4,3,2,1,0]:  nnz(L)=4   nnz(U)=4   total LU nnz = 13   (zero fill)

25 versus 13 on a 5×5 — and the gap widens quadratically with size. AMD (klu_native/amd/) is a greedy heuristic: repeatedly eliminate the node of smallest approximate degree, where "approximate" is the trick that makes it fast — it bounds each node's degree using quotient-graph element absorption rather than recomputing it exactly.

Two details of how KLU uses it:

• AMD orders the symmetrized pattern \(A + A^T\), so it produces one permutation applied to both rows and columns of the block. That's why analyze.rs applies pblk to p and q identically.
• Blocks of size ≤ 3 skip AMD and keep their natural order (analyze_worker's own threshold, matched exactly in analyze.rs) — for a 3×3 there is nothing to gain.

On the real Jacobians, this is where the savings actually come from:

PGM transmission-case (n=22, nnz=124):
    LU nnz with BTF+AMD = 124      (zero fill)
    LU nnz natural order = 308     (184 fill entries)

AMD ordered that Jacobian so well that the factorization produces no fill whatsoever — the LU factors have exactly as many nonzeros as the matrix. Natural order would have produced 2.5× more.

The 3-bus case is n=3, so it takes the natural-order path and both counts are 9 (a 3×3 Jacobian is dense anyway).

Step 3: numeric factorization — Gilbert-Peierls with partial pivoting

Now the values matter. KLU factors each diagonal block with a left-looking LU: column \(k\) of \(L\) and \(U\) is computed completely before column \(k+1\) is touched, using only columns to its left. kernel.rs::factor_block does this for one block.

Take this block (already ordered — pretend BTF and AMD have run):

       c0  c1  c2  c3
r0      2   .   1   .
r1      1   3   .   .
r2      .   1   4   1
r3      .   .   1   5

Each column goes through four sub-steps.

3a. Symbolic: which rows will this column touch?

Before computing anything, KLU determines the pattern of column \(k\) by a depth-first search. The rule (Gilbert-Peierls): the nonzero pattern of column \(k\) of \(L\) and \(U\) is the set of nodes reachable from the nonzero rows of \(A(:,k)\) in the directed graph of the already-computed \(L\).

Column 2 of the example shows why this matters. Its input entries are rows 0, 2, 3 — nothing in row 1. But row 0 is already pivotal (it was column 0's pivot), and column 0 of \(L\) has an entry in row 1. So the DFS reaches row 1 anyway, and column 2 gets an entry there that \(A\) never had. That's fill-in, predicted symbolically before a single flop:

U col 2: [(0, 1.0), (1, -0.5)]
                     ^^^^ fill: A(1,2) = 0

3b. Numeric: sparse triangular solve

With the pattern known, scatter \(A(:,k)\) into a dense workspace and run the updates in topological order — for each pivotal row \(j\) in the pattern, subtract \(x_j \cdot L(:,j)\). Column 2 again: \(x_0 = 1\), and \(L(1,0) = 0.5\), so \(x_1 \mathrel{-}= 0.5 \cdot 1 = -0.5\); then \(L(2,1) = 1/3\) gives \(x_2 = 4 + 1/6 = 4.1\overline{6}\).

3c. Pivot: pick the diagonal if you can live with it

Everything still below the diagonal is a pivot candidate. KLU's rule (lpivot) is diagonal preference with a threshold: take the entry the ordering intended for the diagonal, provided it's at least tol times the largest candidate in the column; otherwise take the largest.

tol defaults to 0.001 (klu_defaults.c, mirrored in types.rs), which is deliberately loose — KLU assumes its input has a strong diagonal and would rather preserve the fill-reducing ordering than chase the last digit of stability. Two 2×2 matrices show both sides:

A = [[0.5, 1],       diag candidate 0.5, column max 1.0
     [1.0, 1]]       0.5 ≥ 0.001 × 1.0  →  keep the diagonal
                     p = [0, 1]   udiag = [0.5, -1.0]

B = [[1e-6, 1],      diag candidate 1e-6, column max 1.0
     [1.0,  1]]      1e-6 < 0.001 × 1.0  →  pivot to row 1
                     p = [1, 0]   udiag = [1.0, 0.999999]

In the first case a dense LU with ordinary partial pivoting would have swapped rows (1.0 > 0.5); KLU does not, because 0.5 is good enough and keeping the row order keeps the sparsity. In the second case the diagonal is hopeless and it swaps — note the resulting \(L(1,0) = 10^{-6}\), tiny, which is the whole point of pivoting.

3d. Prune

After each pivot, prune applies Eisenstat-Liu symmetric pruning: once a column of \(L\) is known to be "covered" by a symmetric counterpart in \(U\), the DFS in step 3a no longer needs to scan past a certain point in it. This changes nothing about the result — it only shortens future searches, and it is why KLU's symbolic step stays near-linear instead of degrading on later columns.

The finished factors

For the 4×4 block above, no pivoting was needed (p = [0,1,2,3]) and the result is:

udiag  = [2, 3, 4.166666666666667, 4.76]

L col 0: [(1, 0.5)]                     U col 0: []
L col 1: [(2, 0.3333333333333333)]      U col 1: []
L col 2: [(3, 0.24)]                    U col 2: [(0, 1.0), (1, -0.5)]
L col 3: []                             U col 3: [(2, 1.0)]

i.e.

\[ L = \begin{bmatrix} 1 & & & \\ 0.5 & 1 & & \\ & \tfrac13 & 1 & \\ & & 0.24 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & & 1 & \\ & 3 & -0.5 & \\ & & 4.1\overline{6} & 1 \\ & & & 4.76 \end{bmatrix} \]

Checking against a dense solve with \(b = [1,2,3,4]^T\): both give \(x = [2/7,\ 4/7,\ 3/7,\ 5/7]\).

Note that \(L\)'s unit diagonal is never stored, and \(U\)'s diagonal lives in its own udiag array — that separation is what lets the solve divide without hunting for the diagonal entry inside a sparse column.

Step 4: solve — permutations, then blocks in reverse

With the factors in hand, solve.rs computes

\[ x = Q \left( (LU + \text{Off})^{-1} , P , b \right) \]

in four moves. Back to the BTF example from step 1, with \(b = [7, 4, 13, 10]^T\):

1. Permute the right-hand side. \(P b\) reorders b by p = [0,2,3,1]:

Pb = [b0, b2, b3, b1] = [7, 13, 10, 4]

2. Solve the blocks in reverse order. This is the part BTF bought us. An earlier block's rows may reference a later block's columns (that's what "upper block-triangular" means), so the last block must be solved first:

block 2  (1×1):   2·y3 = 4                        →  y3 = 2
block 1  (1×1):   5·y2 = 10 − 1·y3 = 8            →  y2 = 1.6
                        ^^^^^^^ off-diagonal entry, now that y3 is known
block 0  (2×2):   [1 3; 4 1]·[y0;y1] = [7; 13 − 2·y3] = [7; 9]

The subtractions are exactly the off-diagonal arrays from step 1 (off_x = [2.0, 1.0], both in the last permuted column) being applied as each block's solution becomes available.

3. Forward/back substitution inside each block. Block 0's factors are \(L = [[1,0],[4,1]]\), \(U = [[1,3],[0,-11]]\) (udiag = [1.0, -11.0]):

forward (Lz = rhs):   z0 = 7,   z1 = 9 − 4·7 = −19
back    (Uy = z):     y1 = −19 / −11 = 1.727273,   y0 = (7 − 3·1.727273) / 1 = 1.818182

4. Permute back. \(x = Q y\), i.e. x[q[k]] = y[k] with q = [1,3,0,2]:

x[1] = y0 = 1.818182
x[3] = y1 = 1.727273
x[0] = y2 = 1.6
x[2] = y3 = 2.0

x = [1.6, 1.8181818181818183, 2.0, 1.7272727272727273]

which matches a dense Gaussian-elimination solve of the original unpermuted matrix to the last digit but one ([1.6, 1.8181818181818181, 2.0, 1.7272727272727273] — the difference is one ulp of round-off, from a genuinely different order of operations).

Step 5: refactor — the one that runs every Newton iteration

Between two Newton iterations the Jacobian's values change but its pattern does not. Refactor exploits that as hard as it can: same BTF blocks, same AMD ordering, same pivot choices, same L/U sparsity pattern — only the numbers are recomputed.

That means refactor skips all of step 3a (no DFS, no reachability search) and all of step 3c (no pivot search). Real KLU's own klu_refactor.c never calls dfs/lsolve_symbolic at all, and neither does the port. The stored pattern of each \(U\) column is already in a valid topological order — it came from the original DFS — so the elimination can just walk it.

Keeping the earlier matrix's pattern and changing its values:

       c0  c1  c2   c3          before:  c0  c1  c2  c3
r0      .   1   .   3.5                   .   1   .   3
r1      .   .   3    .                    .   .   2   .
r2      .   4  2.5   2                    .   4   2   1
r3      5   .  1.5   .                    5   .   1   .

solve 1 (original values):  x = [1.6, 1.8181818181818183, 2.0, 1.7272727272727273]
solve 2 (new values):       x = [1.6, 1.6527777777777786, 1.3333333333333333, 1.5277777777777777]
     dense cross-check:         [1.6, 1.6527777777777781, 1.3333333333333333, 1.5277777777777777]

KluNativeSystem::factor_and_solve_values is the entry point the Newton loop uses: it packs the new values into the cached CSC slots (via groups from step 0), refactors in place, and solves. factor_and_solve is the same path for a caller that still holds full triplets — it reads nothing but their value component. The port's refactor_block_in_place overwrites the existing value fields rather than rebuilding Vecs — profiling showed those per-column allocations were the dominant reason the Rust backend ran ~2× slower than the C one.

The risk of reusing pivots is real but bounded: if the new values make an old pivot choice unstable, the result degrades rather than being caught by a fresh pivot search. factor_and_solve guards against the extreme case by rejecting a non-finite solution, and returns None so the caller can report SolveStatus::Singular.

The whole pipeline on a real Jacobian

The 3-bus network in tests/data/network.json (slack + PV + PQ) gives 2 angle unknowns and 1 magnitude unknown, so a 3×3 Jacobian. At the first Newton iteration from a linear initial guess (vm = [1.06, 1.04, 1.00307], va = [0, 0, −0.05160]):

J = [ 21.834690  −5.298690  −1.462838 ]
    [ −5.119349   9.032698   2.323778 ]
    [  2.005352  −3.538289   8.503522 ]

mismatch = [0.26866152, 0.00368908, 0.00153712]

Through the pipeline:

• BTF: one block of size 3 (r = [0, 3]) — connected network, irreducible.
• AMD: skipped, nk ≤ 3 → natural order. p = q = [0, 1, 2].
• Factor: no pivoting needed (pnum = [0, 1, 2]) — the Jacobian's diagonal dominates, which is the property KLU's loose tol is built around.

k=0  udiag = 21.834690   L = [(1, −0.234459), (2, 0.091842)]   U = []
k=1  udiag =  7.790370   L = [(2, −0.391720)]                  U = [(0, −5.298690)]
k=2  udiag =  9.413792   L = []                                U = [(0, −1.462838), (1, 1.980802)]

• Solve: \(\Delta x = [0.014383106, 0.008478649, 0.000316791]\) — matching a dense solve exactly. The first two entries are angle corrections in radians; the third is a voltage-magnitude correction in per unit.
• Refactor: iteration 2 rebuilds J with the updated voltages and reuses everything above.

What this port deliberately leaves out

src/klu_native/ implements KLU for exactly the configuration gridoxide uses, and the omissions are documented rather than silent:

• Row scaling is ported but not wired in. scale.rs implements both of KLU's variants and is differentially tested against the C, but factor/refactor currently run as if scaling were disabled. Scaling changes which candidate pivots partial pivoting compares — a stability preconditioner, not a correctness requirement — and a per-unit power-flow Jacobian is not pathologically scaled. The differential tests (unscaled Rust vs. scaled C, same matrices) agree to 1e-8.
• Real f64 only, single right-hand side, int32-range indices. No complex arithmetic, no batched multi-RHS, no DLONG. Newton-Raphson needs one real solve per iteration.
• AMD only. COLAMD, user-supplied orderings, and the user-callback ordering are all reachable in real KLU but never selected by gridoxide, so Options has no ordering field to select them.
• No maxwork limit on the maximum transversal — it defaults to "no limit" upstream and is never overridden, so the port always runs the matching to completion.

Where to look in the code

src/klu_native/PROVENANCE.md maps each file to its upstream source and license; the vendored C it was ported from lives in vendor/suitesparse/.

Reading CGMES Input

cargo build --features cgmes builds src/cgmes.rs, a third network-input path alongside the native JSON format and PGM-JSON, reading CGMES (Common Grid Model Exchange Standard) RDF/XML — the IEC 61970/61968 interchange format ENTSO-E and TSOs use.

It is built on cimoxide — a separate Rust project by the same author — for RDF/XML decoding, via its cimoxide-decoder/cimoxide-structs crates, pulled in under their shorter former names (see Provenance and Licensing). The feature is opt-in since some users only need JSON input and shouldn't pay for cimdecoder's dependency tree or build time.

#![allow(unused)]
fn main() {
use gridoxide::cgmes::{load_profiles, cgmes_to_buses_and_branches};
use gridoxide::network::{build_ybus, stamp_shunts};
use gridoxide::run_power_flow_analysis_from_ybus;

let ds = load_profiles(&[&eq_path, &ssh_path, &tp_path, &sv_path])?;
let (buses, lines, transformers, shunts) = cgmes_to_buses_and_branches(&ds, 100e6)?;
let mut ybus = build_ybus(buses.len(), &lines, &transformers);
stamp_shunts(&mut ybus, &shunts);
let result = run_power_flow_analysis_from_ybus(buses, ybus);
}

What the importer expects

The standard EQ+SSH+TP+SV "solved case" profile bundle:

• TP is required. TopologicalNode is used directly as gridoxide's Bus, so switch-state topology processing is assumed already resolved upstream. See Ideal Switches and Zero-Impedance Branches for what that resolution involves and how cgmes::merge_closed_switches handles the node-breaker case.
• SV must carry a populated TopologicalIsland.AngleRefTopologicalNode, used as the slack bus. See Multi-Island Power Flow for how reference buses are picked per island.

What is mapped

Loads — EnergyConsumer, ConformLoad, NonConformLoad, EquivalentInjection, ExternalNetworkInjection, and AsynchronousMachine. The last is converted like a plain load, with both P and Q negated.

Branches — ACLineSegment and SeriesCompensator, including ACLineSegment.gch, real shunt conductance, not just bch's reactive charging.

Transformers — 2- and 3-winding PowerTransformers, with RatioTapChanger (including its optional RatioTapChangerTable per-step override, falling back to the linear stepVoltageIncrement formula when absent) and all four PhaseTapChanger variants: Linear, Symmetrical, Asymmetrical, and Tabular.

Shunts — LinearShuntCompensator and NonlinearShuntCompensator.

Voltage-controlled buses — SynchronousMachine plus RegulatingControl, and the same mechanism for StaticVarCompensator and ExternalNetworkInjection, minus the active-power term for the former.

Validation

Validated end-to-end against four ENTSO-E conformance cases, with fixtures referenced via a git submodule (see tests/data/cgmes/README.md):

MicroGrid-BE-MAS and MiniGrid converge cleanly but match their own published SV voltages only within a few percent. That gap was cross-checked (for MicroGrid-BE-MAS) against pypowsybl's own independent CGMES import and AC load flow on the same case, which shows a comparable deviation from the same published values (scripts/bench/cross_validate_cgmes_microgrid_be.py) — confirming it is inherent to solving a boundary-truncated area file with fixed-injection equivalents, not a correctness bug. One known, documented limitation contributes: types::Line has no tap ratio, so it can't absorb the small nominal-voltage mismatch CGMES explicitly allows at boundary tie points.

Not built or tested in CI — the same local/manual-verification posture as klu and pardiso.

The per-class pages

The remaining pages in this section each take one CIM class or attribute that needed real modeling work, and follow the same structure: why it matters, the concepts and formulas involved, where it sits in gridoxide today, and how other tools handle it.

StaticVarCompensator

Motivation

A Static Var Compensator (SVC) is a shunt-connected, power-electronically-controlled reactive power source: instead of a fixed shunt capacitor/reactor bank, it continuously adjusts its own reactive injection to hold the voltage at its connection point near a setpoint, within a capacitive/inductive rating. Electrically it behaves like a voltage-controlled bus for load-flow purposes — the same \(\vert V_k \vert\)-known, \(Q_k\)-unknown PV formulation the Powerflow page already describes for a generator, just without any active-power term.

CGMES represents one as its own StaticVarCompensator class (a RegulatingCondEq, the same base every SynchronousMachine and ExternalNetworkInjection also derive from), carrying:

• capacitiveRating / inductiveRating — the SVC's reactive range,
• slope — a droop coefficient for voltage-vs-reactive-power regulation,
• q (SSH) — a starting/fallback reactive power value,
• an optional RegulatingControl reference — the same voltage-mode-target mechanism SynchronousMachine uses.

The concepts

1. Ratings are reactances, not powers

capacitiveRating/inductiveRating read like power quantities ("...at maximum capacitive reactive power") but are documented, and universally treated by real tools, as reactance ratings in ohms. They have to be converted to a susceptance, and from there to a reactive-power rating, rather than used directly as a Mvar value:

\[ B = \frac{1}{X_{rating}}, \qquad Q \approx V^2 B \]

The deciding evidence is that powsybl-core's CGMES importer computes exactly 1 / rating to get the susceptance it stores. If capacitiveRating were already a Mvar quantity, taking its reciprocal to get a susceptance would make no dimensional sense. (See "A bug this caught" below for what happens if you believe the doc text instead.)

An absent or zero rating conventionally means unlimited rather than zero — mapped to \(\pm\infty\) (or ±Double.MAX_VALUE), the same "no rating means no limit" convention tap-changer xMin/xMax uses.

2. Three regulation behaviors, in increasing fidelity

Hard voltage pin. An SVC that is actively regulating in voltage mode is treated exactly like a PV bus: its controlled bus's \(\vert V \vert\) is fixed at the target and its \(Q\) is the free variable, clamped to \([B_{min}V^2,\ B_{max}V^2]\). This is the same mechanism as a generator's voltage control — the SVC contributes no \(P\) term, which is the only structural difference.

Droop (slope). A real SVC doesn't hold voltage exactly; it regulates along a droop characteristic, so that absorbing more reactive power comes with a slightly lower terminal voltage. The linearized form folds directly into the voltage equation as an extra term rather than being a post-hoc correction:

\[ V + \text{slope} \cdot Q_{SVC} = V_{target} \]

Standby / dead-band. An SVC may sit idle as a fixed susceptance while voltage stays inside a dead-band, only entering active regulation when voltage leaves it. Because that decision depends on the solved voltage, it can't be made before the solve — it needs an outer loop that toggles the bus between a fixed-susceptance PQ shunt and an active PV pin between passes, the same architectural pattern the Reactive Power Limits page describes for Q-limit switching, applied here to decide whether to regulate at all rather than how far a limit was exceeded.

A tool that implements only the hard pin still solves the common case correctly; droop and standby refine it.

3. The regulated bus need not be the SVC's own bus

RegulatingControl.Terminal is independent of the equipment's own terminal, so an SVC can regulate a remote bus. Two consequences for anything converting one: the bus whose voltage gets pinned is the control terminal's bus, while the per-unit base for converting the ohm ratings into a Q limit is the SVC's own physical bus's rated voltage. The two differ whenever regulation is remote.

4. Not regulating? Fall back to a fixed injection

An SVC with controlEnabled = false, a disabled RegulatingControl, or no RegulatingControl at all is not a voltage-controlled bus. The fallback is the SSH q value as a plain fixed reactive injection — an ordinary PQ contribution.

Where this fits in gridoxide today

src/cgmes.rs's StaticVarCompensator conversion (added alongside the ACLineSegment.gch fix — see the Shunt Conductance page — after this fixture's own SVC was found to be silently dropped entirely) mirrors SynchronousMachine's existing RegulatingControl-driven PV upgrade:

• If the SVC has an enabled, voltage-mode RegulatingControl, the controlled bus (concept 3 — possibly remote) is promoted from PQ to PV and pinned to the target voltage: concept 2's hard pin.
• Otherwise it falls back to concept 4's fixed Q injection from the SSH q value, using the same sign convention SynchronousMachine.q already established empirically (no negation — see that code's own comment for why the doc text alone isn't trustworthy here).

q_min/q_max follow concept 1 — a reactance rating converted to a per-unit reactive-power limit, \(Q \approx V^2 B \approx B_{pu}\) at \(V \approx 1\) pu (the same flat-voltage approximation SynchronousMachine's own min_q/max_q already make), anchored to the SVC's own physical bus's u_rated:

#![allow(unused)]
fn main() {
let z_base = own_bus.map(|b| buses[b].u_rated * buses[b].u_rated / s_base_va);
buses[controlled_bus].q_min = match (sc.inductive_rating, z_base) {
    (Some(x), Some(zb)) if x != 0.0 => zb / x,
    _ => -f64::INFINITY,
};
buses[controlled_bus].q_max = match (sc.capacitive_rating, z_base) {
    (Some(x), Some(zb)) if x != 0.0 => zb / x,
    _ => f64::INFINITY,
};
}

Two deliberate simplifications versus concept 2's fuller model, both consistent with gridoxide's existing scope elsewhere:

• No droop/slope. StaticVarCompensator.slope is not read at all — every regulating SVC is a hard voltage pin, the same simplification gridoxide already makes for SynchronousMachine.
• No standby/monitoring mode. A non-regulating SVC falls back to its fixed q injection permanently; there's no outer loop that would later switch it back into regulation if voltage left some dead-band, since gridoxide's plain newton_raphson doesn't run any outer loop for SVCs at all (only for PV→PQ switching, and only when newton_raphson_enforcing_q_limits is used instead of the default solver).

Tool reference

power-grid-model and lightsim2grid have no SVC concept at all: PGM's only shunt-connected component is Shunt, a fixed admittance, and lightsim2grid's ShuntContainer is a fixed injection stamped straight into the Y-bus diagonal. Both are consistent with their domains — SVCs are a transmission-level device.

Line Shunt Conductance (ACLineSegment.gch)

Motivation

The Powerflow page's Y-bus construction models each line as a π-equivalent: a series impedance plus a shunt admittance split evenly across both ends. Almost every power-flow tool's documentation talks about that shunt term purely in terms of susceptance — line charging capacitance, the reactive effect of a long line's own capacitance to ground — because for the overwhelming majority of real lines the shunt's real part (conductance, representing corona loss or leakage) is negligible or simply zero.

CGMES's ACLineSegment schema doesn't assume that: alongside bch (susceptance) it defines gch (conductance) as a first-class, independent field, "of the entire line section" like every other ACLineSegment electrical attribute. Most real conformance fixtures leave it at zero, but not all — ENTSO-E's own MicroGrid-BE-MAS fixture has two lines (BE-Line_6, BE-Line_2) with non-negligible gch, together worth several MW of real power at nominal voltage.

The concepts

Every tool that models AC lines at all has shunt conductance somewhere in its branch admittance — this isn't a case of some supporting it and others not. What differs is how it's parameterized, and how visible a distinct "conductance" concept is in the data model. Three variants, all carrying the same physical information:

1. A direct conductance field, per line or per end

The most explicit form: a siemens value stored next to the susceptance. Two sub-variants matter for conversion:

• One value for the whole line — CGMES's own form, gch "of the entire line section". A converter must split it across the π-model's two ends, conventionally evenly: \(g_1 = g_2 = g_{ch}/2\), exactly the way \(b_{ch}\) is already split.
• One value per end (g1/g2) — strictly more general, since it can represent an asymmetric line whose two ends carry different shunts. A whole-line value maps into it trivially by halving; the reverse direction loses information unless the two ends happen to be equal.

2. One complex per-end shunt admittance

Instead of naming conductance separately, store a single complex number per end, \(h = g + jb\), and stamp it directly into the Y-bus diagonal: \(y_{11} = y_s + h_{or}\). The conductance is the real part — genuinely present and genuinely solved, just never given its own name. A grep for "conductance" in such a codebase finds nothing, which is a naming fact, not a capability one.

3. Derived from capacitance and a loss tangent

A physically-motivated alternative parameterization: store the shunt capacitance \(c_1\) and a dielectric loss tangent \(\tan\delta_1\), and derive the complex shunt from both:

\[ y_{shunt} = \omega c_1 \tan\delta_1 + j,\omega c_1 \]

so \(g = \omega c_1 \tan\delta_1\). Same information, expressed the way a cable datasheet expresses it. Converting into this form from a raw siemens pair means backing out the tangent as \(\tan\delta_1 = g_{ch}/b_{ch}\) rather than mapping the two fields across directly — and a converter whose source format has no loss-tangent equivalent has no way to produce a nonzero conductance at all.

Where this fits in gridoxide today

Before this fix, types::Line had a b_shunt field and nothing else — src/cgmes.rs's ACLineSegment conversion read bch but never gch, even though its own comment already documented gch as one of the fields "of the entire line section" (a stale comment that got ahead of the code, not the other way around). network::build_ybus had no way to stamp a real shunt term even if the conversion had wanted to.

The fix uses concept 1's whole-line form with the even split, since that's the same convention gridoxide's own bch handling already used and types::Line already has half-open-line self-loop folding logic that only needed a second field threaded through it:

#![allow(unused)]
fn main() {
pub struct Line {
    pub from: usize,
    pub to: usize,
    pub r: f64,
    pub x: f64,
    pub b_shunt: f64, // total line charging
    #[serde(default)]
    pub g_shunt: f64, // total shunt conductance (CGMES ACLineSegment.gch; usually 0)
}
}

#![allow(unused)]
fn main() {
// build_ybus: split shunt admittance (conductance + susceptance) equally to both ends of line
let y_shunt_half = Complex::new(ln.g_shunt / 2.0, ln.b_shunt / 2.0);
y.add(ln.from, ln.from, y_line + y_shunt_half);
y.add(ln.to, ln.to, y_line + y_shunt_half);
}

#[serde(default)] keeps every existing native-JSON and PGM-JSON fixture working unchanged. The PGM importer always sets g_shunt: 0.0: PGM's own line schema is concept 3, and PgmLine has no tan1 field to derive a conductance from — the same "no data, so no effect" stance the rest of that converter takes.

Why this mattered more than "a couple of MW out of a large network"

The missing MW didn't just make voltages a little off everywhere — it concentrated almost entirely into one bus's angle, and nowhere else. BE-Line_6/BE-Line_2 feed directly into the substation hosting MicroGrid-BE-MAS's StaticVarCompensator (see the StaticVarCompensator page), a voltage-magnitude-pinned bus. A pinned bus can absorb a reactive-power mismatch by adjusting its own Q injection, but it has no equivalent slack for an active-power one — so the several MW this fix restores had, before the fix, nowhere to go but that one bus's angle. Cross-validated against pypowsybl's own independent CGMES import (scripts/bench/cross_validate_cgmes_microgrid_be.py, with both tools pinned to the same reference bus so their angles are directly comparable): worst angle deviation across the whole fixture dropped from 0.34° to 0.07° once gch was included — a five-fold improvement concentrated almost entirely at that one substation, exactly where the missing real power was actually flowing in.

Tool reference

PhaseTapChangerLinear

Motivation

A phase-shifting transformer (PST) uses its tap changer to shift the voltage angle across the transformer, not (only) its magnitude — a way to directly control how much active power flows through a particular path in a meshed network, independent of the voltage-magnitude control an ordinary tap changer provides. CGMES models four distinct ways a phase tap changer's per-step behavior can be specified, each its own concrete class: PhaseTapChangerSymmetrical and PhaseTapChangerAsymmetrical (trigonometric formulas), PhaseTapChangerTabular (an explicit per-step lookup table), and PhaseTapChangerLinear — the simplest of the four.

PhaseTapChangerLinear is exactly what its name says: ratio is always exactly 1.0 (a pure phase shifter, no magnitude change at all), and angle is linear in tap step — \(\alpha = (\text{step} - \text{neutralStep}) \cdot \text{stepPhaseShiftIncrement}\) — a mathematical approximation of a real PST's behavior, per the CIM class's own doc comment, rather than a physically derived model like Symmetrical/Asymmetrical's trigonometric curves. ENTSO-E's conformance suite has two dedicated test configurations for exactly this class, PST_PhaseTapChangerLinear_Type1/_Type2 — without support for it, neither could be solved correctly at all.

The concepts

1. Per-step behavior: formula vs. table, normalized at import

The four CGMES classes are four specifications of the same underlying thing: for a given tap position, what complex ratio \(\rho e^{j\alpha}\) (and what r/x/g/b deviation) does the transformer have? The distinction only needs to exist while reading CGMES. Once each class's own rule is evaluated, everything downstream can work from a uniform flat representation — a per-step table of \({\rho, \alpha, r, x, g, b}\) plus a scalar tap position selecting a row, with the originating class kept only as a property for round-trip export.

For PhaseTapChangerLinear the rule is the two lines from the Motivation: \(\rho = 1\) and \(\alpha = (\text{step} - \text{neutralStep}) \cdot \text{stepPhaseShiftIncrement}\), evaluated for every step from lowStep to highStep.

A converter that only ever needs one step's effect — the current SSH position, because it never exports a step table — can evaluate the same rule for that single step and skip building the table at all.

2. Reactance varies with tap position

A phase shifter's leakage reactance is not constant across its tap range. When xMin/xMax are present, \(x\) is interpolated between them by a trigonometric rule shared between the Linear and Symmetrical cases:

\[ x(\alpha) = x_{min} + (x_{max} - x_{min}) \left( \frac{\sin(\alpha/2)}{\sin(\alpha_{max}/2)} \right)^{2} \]

where \(\alpha_{max}\) is the largest angle reachable across the changer's own step range. This value supersedes the PowerTransformerEnd.x from the EQ profile — the reason phase tap changers need an x-override path at all, where an ordinary ratio tap changer usually doesn't.

3. Who moves the tap: fixed input vs. control outer loop

Regardless of which CGMES class produced it, the phase-shift angle is a fixed parameter within any single Newton-Raphson solve — it is not an unknown, and appears in the Jacobian only if a tool wants sensitivities with respect to it. Two designs then exist for choosing that parameter's value:

• Fixed at conversion. The tap position from the input snapshot (CGMES SSH step) is baked in once before the solve, and only an explicit external call changes it between solves.
• Moved by an outer loop. Solve, check whether the controlled branch's active-power or current flow matches its target, adjust the tap position, re-solve — the same architectural pattern the Reactive Power Limits page describes for Q-limit switching. Since tap positions are discrete, an "incremental" refinement computes a continuous \(dP/d\alpha_1\) sensitivity from the Jacobian to estimate how many positions to move at once, bounded to prevent oscillation.

This choice is orthogonal to concepts 1 and 2: a tool can model all four CGMES flavors faithfully and still never move a tap, or move taps aggressively while supporting only one flavor.

Where this fits in gridoxide today

gridoxide takes concept 3's fixed-input option: tap position is chosen once at conversion time from the current SSH step, and there's no outer loop anywhere in the solver that moves a tap to hit an active-power or current target. Every tap changer conversion (RatioTapChanger, all four PhaseTapChanger variants) computes the complex ratio/angle for whatever step CGMES's SSH profile already says it's at and bakes that into the transformer's tap field before Newton-Raphson ever runs.

phase_tap_linear (src/cgmes.rs) is concepts 1 and 2 for the single current step — no step table is built, since nothing in gridoxide exports one:

#![allow(unused)]
fn main() {
fn phase_tap_linear(ptc: &cimstructs::PhaseTapChangerLinear, mrid: &str, xtx: f64) -> Result<TapEffect, CgmesError> {
    ...
    let angle_rad_at = |s: f64| -> f64 { ((s - neutral) * inc_deg).to_radians() };

    let alpha = angle_rad_at(step);
    let tap = Complex::from_polar(1.0, alpha);

    let alpha_max = (low..=high).map(|s| angle_rad_at(s as f64)).fold(f64::MIN, f64::max);
    let x_override = match (x_min_max(ptc.x_min, ptc.x_max, xtx), alpha_max != 0.0) {
        (Some((x_min, x_max)), true) => {
            let ratio = (alpha / 2.0).sin() / (alpha_max / 2.0).sin();
            Some(x_min + (x_max - x_min) * ratio * ratio)
        }
        (Some(_), false) => Some(0.0),
        (None, _) => None,
    };

    Ok(TapEffect { tap, x_override })
}
}

The x-interpolation helper is shared with phase_tap_symmetrical, matching concept 2's "one rule for both classes" — generalized to take raw xMin/xMax rather than a PhaseTapChangerNonLinear reference, since PhaseTapChangerLinear is a structurally distinct, shallower CGMES class carrying its own same-named fields, not a sibling subtype of Symmetrical/Asymmetrical.

Tool reference

RatioTapChanger.RatioTapChangerTable

Motivation

An ordinary (voltage-magnitude-only, non-phase-shifting) tap changer's simplest CGMES representation is RatioTapChanger's own stepVoltageIncrement — a single percentage, giving a ratio linear in tap step: \(\text{ratio} = 1 + (\text{step} - \text{neutralStep}) \cdot \text{stepVoltageIncrement}/100\). Real transformers don't always tap this uniformly, though — non-uniform winding turns per step, or step- dependent leakage reactance, are both real physical effects a single linear coefficient can't capture. CGMES accommodates this with an optional RatioTapChanger.RatioTapChangerTable reference: a plain RatioTapChanger may carry this alongside its own stepVoltageIncrement, pointing to a table of RatioTapChangerTablePoint rows giving each step's ratio (and optionally r/x/g/b deviation) explicitly, rather than via the linear formula.

ENTSO-E's Svedala conformance fixture (national-scale substation-area model, 53 PowerTransformers) uses this extensively: all 11 of its RatioTapChangers reference a table.

The concepts

1. An optional table alongside a formula — not a separate class

This is architecturally different from PhaseTapChangerTabular, and the difference drives everything else on this page. PhaseTapChangerTabular is its own distinct CGMES class: if its table lookup fails there is nothing else to fall back on, so a failed lookup is a genuine data error. RatioTapChangerTable is just an optional reference that a RatioTapChanger may or may not carry — the same object always has stepVoltageIncrement sitting right there.

The consequence is a precedence rule rather than an error path: try the table, fall back to the formula when the reference is absent, the table is empty, or the table is invalid. A missing or unusable table is a normal, expected condition, not a failure.

2. What a table point can override

Each RatioTapChangerTablePoint carries ratio plus optional r, x, g, b deviations (in percent) for its step — the same TapChangerTablePoint shape PhaseTapChangerTabular's own points use. A tool can therefore honor anywhere from just the ratio up to the full impedance deviation, which is a scope choice independent of concept 1's precedence rule.

3. Full step table vs. current step only

Whether a converter needs to materialize the whole table depends on what it does downstream. A tool that exports models, or that moves taps during the solve, needs every step. A tool that only ever solves the snapshot it was handed needs exactly one row — the one matching the SSH step — and can look it up directly.

Where this fits in gridoxide today

ratio_tap_table (src/cgmes.rs) implements concept 1's precedence rule, with concept 3's single-step lookup (gridoxide only ever needs this step's effect, since it neither exports step tables nor moves taps — see the PhaseTapChangerLinear page for the latter):

#![allow(unused)]
fn main() {
fn ratio_tap_table(ds: &CimDataset, table_mrid: &str, step: i64, xtx: f64) -> Option<TapEffect> {
    for pt_mrid in by_type(ds, "RatioTapChangerTablePoint") {
        let pt: &cimstructs::RatioTapChangerTablePoint = get(ds, pt_mrid)?;
        let Some(owner) = &pt.ratio_tap_changer_table else { continue };
        if owner.mrid != *table_mrid || pt.base.step != Some(step) {
            continue;
        }
        let ratio = pt.base.ratio.unwrap_or(1.0);
        let x_pct = pt.base.x.unwrap_or(0.0);
        return Some(TapEffect { tap: Complex::new(ratio, 0.0), x_override: Some(xtx * (1.0 + x_pct / 100.0)) });
    }
    None
}
}

Returning None (rather than an error) when the table or a matching point for the current step isn't found is the deliberate difference from phase_tap_tabular's own hard-error behavior — exactly concept 1: the caller falls back to the linear stepVoltageIncrement formula instead of treating a missing/invalid table as a data error.

On concept 2, gridoxide reads only ratio and x. r/g/b deviations have no representation in TapEffect at all ({ tap: Complex<f64>, x_override: Option<f64> } — no r_override, no g/b override), a pre-existing simplification across every tap-changer conversion in this file, not something specific to this feature.

Validated against Svedala — but not a strong before/after signal on this particular fixture

Solving ENTSO-E's Svedala conformance case end-to-end (191 buses after this fixture's own 3-winding transformer star-bus synthesis on a different fixture with the same feature) converges cleanly in 6 Newton-Raphson iterations, matching 108 published SvVoltage values with a mean absolute error of ~0.53% and a worst case of ~4.35%.

Disabling the table lookup (forcing every RatioTapChanger back onto its linear fallback, as a direct A/B comparison) barely moves either number — this fixture's own tables happen to have ratios that are exactly linear already (confirmed directly against the raw EQ XML: table step 1 gives ratio 0.88, matching \(1 + (1-13)\times 1/100 = 0.88\) exactly, all the way through step 25). So Svedala doesn't demonstrate a large accuracy win from this feature specifically — its real value is CIM-spec correctness for the (real, if less common) case where a table genuinely diverges from the linear approximation, which this particular fixture just doesn't happen to exercise.

Tool reference

Of the tools surveyed, only these two have a CGMES RatioTapChangerTable path at all: power-grid-model and lightsim2grid have no CGMES import, and VeraGrid's and pandapower's CIM importers weren't surveyed for this specific reference.

ExternalNetworkInjection

Motivation

Any CGMES sub-model that covers less than the full interconnected synchronous area — which is to say, almost every real CGMES file, since even a large national TSO's own model is a part of Continental Europe's actual network — needs some way to represent what lies beyond its own boundary. CGMES has (at least) two distinct classes for this: EquivalentInjection and ExternalNetworkInjection. Both carry P/Q; the difference is that ExternalNetworkInjection additionally derives from RegulatingCondEq (the same base StaticVarCompensator and SynchronousMachine share), so it can also carry a RegulatingControl and behave like a voltage-regulated source, not just a fixed injection. ENTSO-E's MiniGrid conformance fixture uses it this way for its two external-grid connection points ("Q1"/"Q2").

The concepts

1. An external injection is a generator, not a new element type

Structurally, "the rest of the interconnected system" behaves exactly like a generator at the boundary bus: it injects P and Q, it may hold a voltage setpoint, and it has (possibly unbounded) reactive limits. Nothing about it needs a distinct element in the internal network model — the same RegulatingControl-driven PV upgrade a SynchronousMachine gets applies unchanged, with the originating CGMES class kept only as a tag for round-trip export.

This applies to EquivalentInjection and ExternalNetworkInjection equally: the two CGMES classes need not map to two internal types.

2. Sign convention: both P and Q are negated

CGMES's SSH power-flow values for these classes are given in the load sign convention — power drawn from the network at the terminal — so converting to an injection means negating both:

\[ P_{inj} = -p_{SSH}, \qquad Q_{inj} = -q_{SSH} \]

This is worth stating explicitly because it is not universal across CGMES's regulating-equipment classes: SynchronousMachine's own q follows a different rule in practice (no negation), so "it derives from RegulatingCondEq, therefore it signs like a machine" is exactly the wrong inference. The right grouping is by conceptual role: an external injection stands in for a network, not for a physical rotating machine.

3. Boundary-point folding applies to only one of the two classes

An EquivalentInjection sitting exactly at a network boundary point can be folded into the virtual generation/load of the boundary line that crosses it, rather than becoming a standalone element. ExternalNetworkInjection has no equivalent path, because it isn't itself a boundary-defining class — it is always a standalone element at an ordinary bus. This is the one real structural difference between the two, and it only matters to a tool that models boundary lines as first-class objects.

4. Under-specified limits interact with slack selection

CGMES often leaves an external injection's active-power limits unspecified, and the natural default is "unbounded" (\(\pm\infty\), or ±Double.MAX_VALUE). Any heuristic that picks a slack bus by generator size then has to cope with a generator whose stated maxP is absurd — the usual treatment is a plausibility threshold above which a candidate is excluded. So an under-specified external injection tends to be filtered out of slack consideration rather than favored, purely as a side effect of the defaulting.

Where this fits in gridoxide today

The conversion (src/cgmes.rs) applies concept 1 by mirroring StaticVarCompensator's block almost exactly — same RegulatingCondEq/RegulatingControl-driven PV-bus upgrade, same fallback to a fixed injection when not actively regulating — with concept 2's sign convention:

#![allow(unused)]
fn main() {
for mrid in by_type(ds, "ExternalNetworkInjection") {
    let eni: &cimstructs::ExternalNetworkInjection = require(ds, mrid, "ExternalNetworkInjection", mrid, "(self)")?;
    ...
    buses[bus].p_spec += -eni.p.unwrap_or(0.0) * 1e6 / s_base_va;
    buses[bus].q_spec += -eni.q.unwrap_or(0.0) * 1e6 / s_base_va;
    ...
}

Both P and Q are negated here — EquivalentInjection's convention, not SynchronousMachine's Q exception (see the StaticVarCompensator page for why that exception exists and why it doesn't extend here).

Concepts 3 and 4 don't arise: gridoxide has no boundary-line object for an injection to fold into, and its slack bus comes from CGMES's own TopologicalIsland.AngleRefTopologicalNode (see Multi-Island Power Flow) rather than from any generator-size heuristic, so an unbounded maxP has nothing to distort.

Tool reference

powsybl-open-loadflow's only indirect sensitivity is concept 4: ExternalNetworkInjectionConversion defaults maxP to ±Double.MAX_VALUE, and LargestGeneratorSlackBusSelector filters out generators whose maxP exceeds a plausibility threshold when choosing a slack candidate.

power-grid-model and lightsim2grid have no CGMES import, so neither has an equivalent concept.

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. gridoxide's own column was re-verified against current source (every cell claiming support names the function or type implementing it), and the new "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 pre-existing row were not re-surveyed and are carried over from the previous revision — treat them as the older snapshot. The new "Multi-island" row marks the four tools not checked as not surveyed rather than ❌, since absence of a survey is not evidence of absence of the feature.

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.
• 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 — solver::newton_raphson_enforcing_q_limits 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.
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. What remains is the part the batch fast path deliberately excludes: a branch outage gives each scenario its own Y-bus and therefore its own sparsity pattern, so Scenario::branch_outages is currently a documented error rather than a solve.
3. Distributed slack — 4 of 5 tools have it (only lightsim2grid, powsybl, VeraGrid, and pandapower); real transmission grids often split slack across several generators.
4. DC power flow as a first-class mode — cheap, since linear_initial_guess is most of the way there already; every reference tool treats this as a basic offering.
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. Everything else in the table (TCSC/SSC and the wider FACTS set, short-circuit, sensitivity/OPF, outer-loop/controller architecture as a general extensibility mechanism, VeraGrid's unique EMT/RMS dynamic simulation, pandapower's protection-device modeling) — real capabilities elsewhere, but either a materially larger undertaking or outside gridoxide's current scope as a focused AC power-flow library.

Note on state estimation

Done for symmetric estimation — snapshot and batched, with voltage, power and current sensors — and for asymmetric networks driven by symmetric sensors. 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.

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, 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. Substantially done, and what remains is narrower than the heading suggests.

   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 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.

   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.

   pgm_3ph_maps refuses the components the three-phase conversion does not model — link, three_winding_transformer, voltage_regulator, and any transformer winding pair outside Dyn and YNyn — with a typed error rather than dropping them silently or, in the last case, panicking from inside transformer_seq_params.

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 (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: every benchmark and fixture here estimates from data that is either perfectly consistent or hand-authored. Nothing in this repo generates realistically noisy or corrupted measurements, so gridoxide's bad-data advantage — lead 1 — has never actually been measured against anything. Bad-data behaviour needs a harness that does not exist yet.

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 newton_raphson_enforcing_q_limits 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.

Benchmarking and Profiling

scripts/bench/README.md is the single source of truth for every benchmark number in this project. This page is a map into it, not a copy of it — the numbers live there, next to the scripts that produce them, so they can be updated in one place when re-measured.

Profiling

For profiling with perf, set:

sysctl kernel.perf_event_paranoid=1

What is measured, and where

scripts/bench/README.md is organized as a numbered sequence of benchmarks:

The two benchmark shapes

Synthetic radial distribution grid (§1–3). examples/bench_network.rs and scripts/bench/bench_gridoxide_native.py time gridoxide against power-grid-model on generated MV/LV topology at controllable scale. Its cold mode measures N independent flat-start solves with no shared state; the optional warm mode measures repeated solves through a PersistentSolver — see Backends and Factorization Reuse.

Real power-system test cases (§4). Twelve real IEEE/MATPOWER grids, 14 to 9,241 buses, comparing gridoxide against five independent solvers: power-grid-model, lightsim2grid, RTE's powsybl-open-loadflow (via pypowsybl), pandapower's default solver, and VeraGrid.

Two results from that second benchmark are worth stating here because they shaped the code:

• gridoxide and pandapower's own native path are the only two of the six that converge on all 12 cases. The other four each fail on a subset of the same handful of genuinely hard cases (RTE's own real production grids), confirmed by cross-checking against powsybl-open-loadflow directly — not a gridoxide gap.
• Compared warm-vs-warm, Klu is frequently faster than lightsim2grid's own KLU-backed C++ solver on this real transmission-topology data, even though PGM still clearly beats every gridoxide backend on the synthetic radial-distribution topology. The comparison genuinely depends on grid topology, not just implementation language.

Provenance and Licensing

gridoxide's own code is licensed under Apache-2.0 (LICENSE). Several pieces of third-party code are vendored, translated, or linked, and each carries its own terms. This page is the summary; the authoritative per-file detail lives in the PROVENANCE.md files kept alongside the code they describe.

The crate-wide license field

Cargo.toml's license field is:

Apache-2.0 AND BSD-3-Clause AND LGPL-2.1-or-later

That is accurate for every default cargo build, not just an opt-in one — because src/klu_native/ is always built, with no feature gate.

Always built: src/klu_native/

src/klu_native/ is a from-scratch Rust translation of vendored SuiteSparse AMD, BTF, and KLU C source. A close translation of licensed source is reasonably a derivative work regardless of implementation language, so it carries forward its upstream license — and that is not one license here:

src/klu_native/PROVENANCE.md has the exact file-by-file mapping back to the upstream C, including which specific functions each Rust file was ported from and how shared header material was classified. See Inside KLU for what the ported algorithm actually does.

Opt-in: --features klu

Building with cargo build --features klu additionally compiles the vendored SuiteSparse C itself into the binary via FFI. vendor/suitesparse/ is a partial vendoring of SuiteSparse at tag v7.12.2 (commit 42151688813c45846a597edcb601435a0e38f3dd, 2026-02-10) — only the Source/ and Include/ subdirectories of five packages (SuiteSparse_config, AMD, COLAMD, BTF, KLU), each keeping its own Doc/License.txt:

This adds no license beyond what is already listed above, but it does add LGPL's relinking obligations for anyone distributing a binary built with that feature. The klu-dynamic sub-feature exists for that case: it links a system-installed libklu.so instead of statically linking the vendored copy. See vendor/suitesparse/PROVENANCE.md for exactly what was and was not vendored, and how to update to a newer SuiteSparse release.

Opt-in: --features pardiso

A separate case from all of the above. It dynamically links a locally-installed Intel oneMKL (libmkl_rt.so) at build and run time, under Intel's own Simplified Software License — not LGPL, not OSS, and not vendored or redistributed by this repo in any form. No MKL header or source is copied in; bindgen only reads the local install's own mkl_pardiso.h at build time to generate FFI bindings.

Because nothing MKL-derived is ever copied into or shipped by this crate, Cargo.toml's license field does not change for this feature. Anyone who builds with --features pardiso and distributes the resulting binary is responsible for their own compliance with Intel's oneMKL redistribution terms — this project doesn't audit that on their behalf.

Opt-in: --features cgmes — the cimoxide dependency

The optional cgmes feature (src/cgmes.rs) depends on cimoxide's decoder and generated-structs crates, a separate Rust project by the same author providing CGMES RDF/XML decoding into typed CIM structs.

• Upstream: https://github.com/m-mirz/cimoxide
• Crates: cimoxide-decoder and cimoxide-structs, both 0.3.0

cimoxide is Apache-2.0, matching gridoxide's own license — no licensing mismatch from this dependency.

The package = rename

Cargo.toml pulls both in under their old, shorter names:

cimdecoder = { package = "cimoxide-decoder", version = "0.3.0", optional = true }
cimstructs = { package = "cimoxide-structs", version = "0.3.0", optional = true }

The crates.io names are cimoxide--prefixed, but src/cgmes.rs and the CGMES tests were written against cimdecoder/cimstructs, so package = keeps every existing use path valid rather than renaming 33 references for no behavioural gain.

Formerly a git dependency

Until cimoxide's crates.io release, these were a git dependency pinned to a vendor/gridoxide branch — main plus one commit force-adding the code generator's normally-gitignored output, since cimstructs's source is produced by cargo run -p cimgen from ENTSO-E's RDF/SHACL schemas and so does not exist on a fresh main clone.

Updating

To pick up a newer cimoxide schema or generator change:

1. Release the new version from the cimoxide repo (make generate, then make build && make test before publishing — cimstructs's generated source is what a schema change moves).
2. Bump the version = value in this repo's Cargo.toml for the cimdecoder and cimstructs dependencies, and update the version noted above.
3. Run cargo build --features cgmes and cargo test --features cgmes to confirm the converter still matches tests/cgmes_microgrid_be_test.rs's expectations — CGMES field names or shapes could in principle change between schema versions. The fixture-backed CGMES tests silently skip when tests/data/CGMES-Test-Configurations isn't checked out, so git submodule update --init tests/data/CGMES-Test-Configurations first, or the run proves nothing beyond "it compiles".

Test fixtures

The CGMES conformance fixtures under tests/data/cgmes/ are referenced via a git submodule rather than committed, because of their own licensing — see tests/data/cgmes/README.md.