Benchmarking CPU, GPU, and ML Power-Flow Solvers from 14 to 70,000 Buses

Introduction

Power flow (PF) solving is a computational bottleneck both for grid operations and for ML pipelines that rely on PF solvers to generate training data at scale: rising renewable penetration is increasing the number of instances operators must evaluate, and recent large-scale neural surrogate efforts report data-generation budgets that are themselves a substantial compute cost.

A large and active line of machine-learning work responds to this need for speed with learned surrogate models [15, 14]. These comparisons are established against CPU-based baselines, since batched GPU-based solving has not, until recently, been practical: the sparse factorizations at the core of a Newton-Raphson solve mapped poorly onto GPU hardware at the scale of a single instance. Batched sparse solvers such as cuDSS [11] make GPU-based physics solving newly practical, which raises the question we investigate: how does the ML throughput case, established mainly against CPU baselines, hold up once that option exists?

Comparisons between ML models and a single CPU-based baseline can be misleading, as CPU solvers themselves vary substantially in performance. We observe throughput differences of up to ∼2.5–3 orders of magnitude across CPU implementations, so comparing against just the slowest or fastest baseline can yield opposite conclusions on whether ML is “needed” for speed. We built a benchmark spanning more than 3.5 orders of magnitude in grid size (14 to 70,000 buses): 7 CPU-based tools, 2 GPU-based physics solvers (including gpusim2grid, introduced here), and 3 ML surrogates, across multiple languages and academic/industrial tools – to our knowledge the first time PF computation is evaluated at this scale against this many alternatives on common ground. Code, Footer: 40th Conference on Neural Information Processing Systems (NeurIPS 2026). Workshop: AI Foundations for Power Grids. including the reference KCL evaluator (§ 2.1) used to check every solver on equal terms, is available at an anonymized repository¹.

We find that a sufficiently parallelized CPU-based solver, or gpusim2grid batch physics solver (§ 5.3), matches or exceeds the throughput of the ML surrogates that we evaluate on the power-flow task while being orders of magnitude more accurate. GENCO [14] reports genuine gains for learned surrogates on more computation-heavy tasks such as Optimal Power Flow (OPF) (§ 2); we read the two as complementary rather than contradictory. In summary, our two contributions are: (1) the benchmark suite described above, addressing the reproducibility gap left by closed-source GPU accelerators in this space; and (2) gpusim2grid itself, matching the precision of CPU-based solvers and exceeding the throughput of the ML surrogates.

Problem Setting and Related Work

Power flow formulation

We consider the classical AC power-flow problem, whose goal is to find, for each bus $i$, the complex voltage $v_i$ satisfying the equation:

where $Y \cdot v$ is the matrix product of the admittance matrix $Y$ and $v$, $^*$ is the conjugate transpose, and $\odot$ is the Hadamard product. We follow the bus-type conventions and mismatch definition of Zimmerman et al. [2011] (also used by Rivera et al. [2025], Thurner et al. [2018]).

Related work

GPU-based acceleration of power flow is an active line of work. ExaPF.jl [17] moves Newton-based AC power-flow solves onto GPU. gpusim2grid instead targets the complementary, batched regime relevant to contingency analysis and training-data generation, alongside three other efforts in this same batched, sparsity-aware GPU regime: SABLE [12], a cuDSS-based batched refactorization scheme, not yet open-sourced to our knowledge. A concurrent JAX/GMRES implementation [22] and an earlier batched sparse-LU approach [19] both lack publicly available code, to our knowledge. None release code or are evaluated against a common set of baselines – precisely the reproducibility gap our benchmark suite tries to close.

Concurrently, GENCO [14] proposes a neural solver unifying PF, OPF, and State Estimation (ST), reporting up to 30× throughput over a PF baseline in a hardware-matched CPU/GPU comparison, mainly on grids $\ge2{,}000$ buses at a looser tolerance than our CPU baseline. We do not dispute this comparison; together with our results (§ 3), it suggests learned surrogates show their clearest gains on OPF or ST rather than plain PF, where gpusim2grid is, to our knowledge, the only openly released implementation benchmarked across a comparably wide grid size and solver range.

Experiments and results

Common experimental settings

Every solver benchmarked below solves the same AC power-flow system on the same 15 grids (14 to 70,000 buses)²: 7 open-source Newton-Raphson CPU-based tools – pandapower (pp) [18], powermodels [3], exapf-CPU (exapf) [17], pypowsybl/OLF (olf) [RTE and pypowsybl contributors, 2026, LF Energy Foundation, 2026], light-sim2grid (ls) [6], power_grid_models (pgm) [21, 20], and rust-power [Cheng, 2026] – 2 open-source Newton-Raphson GPU-based physics solvers (exapf-GPU, gpusim2grid), and 3 ML surrogates (PowerFlowNet [10], GNS-S [6], CANOS-PF [13]). We used Park et al. [2026]’s data-sampling methodology by default.

Physical solver validation

Experimental design. This first set of experiments checks that all benchmarked solvers solve the same problem, to ensure fairness. To that end, we wrote an independent script implementing the KCL evaluator of § 2.1 (detailed in § 5.2 and developed directly from the matpower reference manual [23]) to check every solver’s output under a common configuration (per-solver details in § 5.4).

We report the “mean/worst-case KCL mismatch”: with $\operatorname{mis}=v\mathbin{\odot}(Y.v)^{*}-s$, the violation vector is $[\Re(\operatorname{mis}(pv)),\Re(\operatorname{mis}(pq)),\Im(\operatorname{mis}(pq))]$; the mean/worst-case mismatch is the average/max, over the dataset, of its mean absolute value over the grid.

We sampled 9,000 independent scenarios per grid, identical across baselines, solved at tolerance $10^{-8}$ [12]; mismatch is computed over converged samples only.

Results. At 5 iterations, gpusim2grid is as precise as the reference CPU solver in all cases (9,000 scenarios $\times$ 15 grids): worst-case KCL mismatch is $2.31 \times10^{-10}$, below the $10^{-8}$ tolerance (Table 1); all reported timings in the next sections use 5 iterations for gpusim2grid. All Newton-Raphson-based solvers solve the PF problem, independently validated by our evaluator (full per-solver results in Table 3, § 5.6).

Solver	Grids (max size)	Convergence rate	Mean mismatch	Worst-case mismatch
pandapower (ref)	15 (70000)	100.00%	$5.61\mathrm{e}{-11}$	$1\mathrm{e}{-08}$
gpusim2grid (4 iter)	15 (70000)	NA	$1.05\mathrm{e}{-09}$	0.000116
gpusim2grid (5 iter)	15 (70000)	NA	$8.91\mathrm{e}{-14}$	$2.31\mathrm{e}{-10}$
gpusim2grid (6 iter)	15 (70000)	NA	$8.8\mathrm{e}{-14}$	$1.04\mathrm{e}{-10}$

Solver	#Grids (max size)	Cvg. rate	Mean mismatch	Worst-case mismatch
pp	15 (70000)	100.00%	5.61e-11	1e-08
powermodels	14 (25000)	78.57%	9.78e-13	1.26e-08
exapf	14 (25000)	100.00%	9.11e-14	8.27e-11
olf	15 (70000)	100.00%	4.4e-10	1.75e-06
ls	15 (70000)	100.00%	1.4e-13	1e-10
pgm	5 (500)	99.88%	1e-14	1.11e-12
rustpower	15 (70000)	100.00%	1.31e-10	8.65e-09
gpusim2grid (0 iter)	15 (70000)	NA	0.0142	4.39
gpusim2grid (1 iter)	15 (70000)	NA	0.00815	45.4
gpusim2grid (2 iter)	15 (70000)	NA	0.000365	4.5
gpusim2grid (3 iter)	15 (70000)	NA	3.08e-06	0.0806
gpusim2grid (4 iter)	15 (70000)	NA	1.05e-09	0.000116
gpusim2grid (5 iter)	15 (70000)	NA	8.91e-14	2.31e-10
gpusim2grid (6 iter)	15 (70000)	NA	8.8e-14	1.04e-10
gpusim2grid (7 iter)	15 (70000)	NA	8.8e-14	1.04e-10

Timings on generation and demand varying scenarios

Experimental design. This section measures each solver’s maximum throughput – Newton-Raphson on CPU, Newton-Raphson on GPU, or a learned surrogate on GPU – pushed to its limit, on a task common in the literature (e.g., [12]): varying the generation and load profile for a fixed topology (details in § 5.7). Following [14], we report the best achieved throughput for each of the 12 benchmarked methods³; Figure 1 shows the fastest method in each of the 3 categories (NR/CPU, NR/GPU, learned model).

This figure is favorable to the benchmarked learned models: a fixed architecture (the one proposed in the PF-$\Delta$ paper, trained only up to 2000 buses) is reused unmodified on every grid in our suite, without retraining or resizing for the largest cases. Our throughput comparison therefore bounds what this specific, non-adapted surrogate achieves, not necessarily what a resized or re-trained model could reach on grids far outside its training range, in either throughput or accuracy. This part of our benchmark focuses on throughput; accuracy is compared in § 3.4.

Results. Figure 1 shows a large gap between the fastest CPU solver (lightsim2grid for this benchmarked task) on a single core and at maximum throughput⁴ – on average 16$\times$ faster – which matters when comparing a GPU to one CPU core rather than to the full CPU: the fastest learned model (the “Graph Neural Solver”) beats 1 CPU core of the fastest CPU baseline, but not the full CPU. The evaluated NR algorithm on the GPU (gpusim2grid), however, matches the CPU’s throughput on the largest grids.

PF-Delta and accuracy

Experimental design. Table 2 reports the Power Balance Loss [15] for pandapower and gpusim2grid on the PF-$\Delta$ datasets⁵, alongside published accuracy numbers for the best ML surrogates among PowerFlowNet, Graph Neural Solver, CANOS-PF and GENCO Base (the best one) in their best configuration, from their source publication.

Results. In this setting – most favorable to the ML surrogates (best model, evaluated on its best-performing task), matching train/test distributions – pandapower and gpusim2grid are far more accurate on case118 across $N/N$-1/$N$-2: over 1000× vs. the most accurate ML model (GENCO Base, 2.7e-3, against pandapower’s 1.1e-6). C2I results are skewed by three non-converging outliers (see table note); gpusim2grid achieves near-full convergence within the same budget.

Conclusion and discussion

Our results show a well-parallelized CPU-based solver or gpusim2grid outperforms the GNN-based ML surrogates evaluated for power flow specifically. This does not concern all learned methods (or architectures), and our findings are consistent with GENCO’s reporting higher gains on OPF and state estimation [14]. These results do not concern the use of GNNs for other tasks such as tertiary voltage control [5] or grid topology optimization [8].

All the learned methods benchmarked are GNN-based, and other architectures, e.g. transformers, might trade off speed and accuracy differently. We also report a single best-observed throughput per solver/grid rather than repeated-run statistics, following common practice at this scale (e.g. Puech et al. [ [14], Park et al. [ [12]); repeating our 630M solves multiple times each was not tractable within our compute budget. Our benchmark was also limited by the machines we had available, (e.g., a more recent or larger CPU or GPU). The synthetic data used here also does not fully capture real-grid complexity; evaluating on GENCO’s realistic dataset is an interesting future direction.

LLM usage disclosure. This manuscript’s text was drafted and revised with the assistance of a large language model (Claude, Anthropic), which also helped write the data-analysis and figure-generation scripts used to visualize the benchmark results. Beyond manuscript preparation, the model was also used as part of the research method itself: it assisted the authors in implementing the gpusim2grid solver (§ 5.3) – the batching strategy and its CUDA/cuDSS-based implementation – under the authors’ direction, review, and testing at every stage. It was not involved in the benchmark’s experimental design (choice of solvers, grids, sampling protocol, or metrics) nor in producing any reported numerical result: all timing runs were launched by the authors, accuracy figures were generated by the authors’ own pipeline and independently verified against the reference KCL evaluator of § 2.1. All LLM-assisted content, including code, was reviewed by the authors, who take full responsibility for the manuscript’s accuracy and integrity.

AI4PowerGrids Domain Checklist

1. Physics feasibility. All results reported for every solver use the full AC power-flow equations (Eq. (1)), not a DC or linearized approximation. Feasibility is assessed with an independent reference KCL evaluator (§ 2.1, detailed in § 5.2), reporting mean and worst-case constraint-violation statistics rather than aggregate error alone (Table 1, full per-solver results in Table 3).

2. Out-of-distribution evaluation. No solver benchmarked here is trained by us, and the physics-based solvers require no training at all, so “out-of-distribution” in the usual train/test sense does not directly apply to them. That said, the benchmark itself is a structural distribution shift: all methods, including the published ML surrogates, are evaluated across 15 grids spanning 3.5 orders of magnitude in size (14 to 70,000 buses, § 3.1), and, for the PF-$\Delta$ comparison, across four topology/operating regimes – normal ($N$), single- and double-contingency ($N$-1, $N$-2), and close-to-infeasible (C2I) – none of which match the ML surrogates’ own training conditions (§ 3.4, Table 2). We report this as the closest available proxy for OOD evaluation given that no method here was trained by us.

Model	$N$	$N$-1	$N$-2	Close-to-infeasible (C2I)
pandapower (ref.)	$1.1\mathrm{e}{-6}$	$8.5\mathrm{e}{-7}$	$7.5\mathrm{e}{-7}$	$7.2\mathrm{e}{-3}$*
gpusim2grid (ours)	$4.9\mathrm{e}{-14}$	$2.3\mathrm{e}{-8}$	$3.1\mathrm{e}{-10}$	$4.5\mathrm{e}{-4}$
GENCO Base†	$2.7\mathrm{e}{-3} \pm2.1\mathrm{e}{-3}$	$4.5\mathrm{e}{-3} \pm3.4\mathrm{e}{-3}$	$5.8\mathrm{e}{-3} \pm4.3\mathrm{e}{-3}$	$2.2\mathrm{e}{-1} \pm8\mathrm{e}{-2}$

3. Failure modes. Convergence rate per solver is reported in Table 1, with the complete per-solver, per-grid breakdown and discussion of failure modes in Table 3 (§ 5.6), including power grid models’ failure to scale beyond small grids and the GPU NaN-mismatch cases tied to a cuDSS beta issue. Diverging cases are explicitly retained (not filtered out) in Table 2 and discussed in the accompanying results paragraph (§ 3.4).

4. System-level effect. Not applicable. This paper benchmarks standalone AC power-flow solving – a computational building block used by every grid operator – and makes no claim about downstream effects when PF is embedded in a larger control loop (e.g., contingency screening, OPF, or real-time dispatch). We state this explicitly as a scope limitation in § 4 rather than leaving it unaddressed.

5. Data and code. Code for gpu-0DCD is already available at an anonymized repository (https://anonymous.4open.science/r/gpusimulator-0DCD), which will also host the independent KCL evaluator and the full benchmark harness (physics-based and ML solvers) by the time of review. We do not plan to release the sampled scenario data itself: the sampling methodology is simple enough to reproduce locally, and regeneration is cheaper than storage or transfer. As a reference point, extrapolating storage costs from the GridFM case14 dataset (https://huggingface.co/datasets/gridfm/pf_small_case14_ieee/tree/main, 400MB for 200,000 instances) to the 8,928,571 case14 scenarios solved for this paper would require roughly 178GB; regenerating that same data on demand takes about 1 minute on a single core and under 3 seconds using multiple cores with the fastest solver benchmarked here.

Supplementary material

Power grids used for this experiment

We tried to use available grids in the literature that were used in published benchmarks, for example in [12] or [14], while also using the smallest and largest possible grids. We finally used:

• case14, case_ieee30, case57, case118 from [23] have also been used by the PF-$\Delta$ paper

• case1354pegase, case9241pegase from [Fliscounakis et al., 2013], case3120sp from [23] and case_ACTIVSg10k from [1] were used in SABLE

• case_ACTIVSg500, case_ACTIVSg2000, case_ACTIVSg25k, case_ACTIVSg70k from [1], case2869pegase, case6515rte from [Fliscounakis et al., 2013] and case3375wp from [23]

They cover a huge range of sizes and complexity, spanning from the “toy model” case14 up to the ACTIVSg70k test case, built to represent the eastern United States (52 areas in total) with 70 thousand buses.

Independent KCL evaluation metrics

To make sure that we were comparing solvers solving the exact same problem, we built a dedicated script that tests whether or not a given solution (provided by any solver) meets Kirchhoff’s Current Law (KCL).

This validation takes as input a matpower case file (.m file), which is the file most of the benchmarked solvers can read.

It then builds the $Y$ matrix from the bus (columns GS and BS) (adding the shunt diagonal entries) and the branch matrix (both diagonal and non-diagonal entries) using the BR_B, BR_X, BR_R, TAP and SHIFT columns (only for lines that are actually connected BR_STATUS $\ne0$), as described in the MATPOWER manual, Section “3 Modeling”, pages 22 and following, especially subsection “3.2 Branches”.

The $s$ vector is built similarly from the same manual, especially subsection “3.3 Generators” and “3.4 Loads” (pages 25 and following).

Finally, the type of each bus (“slack”, “pv” or “pq”) is not read directly from the matpower bus data, because that is not what most of the tested CPU baselines do today. To determine the type of a bus (“slack”, “pv” or “pq”), the following rule is used:

• if it is tagged as “slack” in the matpower bus data, then it is the slack; we do not change that

• a bus is said to be “pv” if it has a connected generator (GEN_STATUS $\ne0$), discarding the input matpower data for this field

• if a bus is neither the “slack” nor “pv”, it is declared “pq”.

Dive into gpusim2grid

gpusim2grid is built on top of CUDA and cuDSS [11], NVIDIA’s batched sparse direct-solver library, chosen for being, to our knowledge, the most mature batched sparse direct-solver library currently available for GPU. The trade-off is that cuDSS is CUDA-only, so gpusim2grid as implemented does not run on AMD ROCm hardware, unlike exapf-GPU in our benchmark; targeting a ROCm-compatible batched sparse solver, should one reach comparable maturity, is straightforward in principle but out of scope here.

In particular the gpusim2grid model implements a full AC powerflow on the GPU. A converged grid (on the CPU) is uploaded once to the GPU, which then:

• extracts from this converged grid all the elements required to solve the system 1: $Y$, $s$, “slack”, “pv” and “pq” indices

• then extracts the data required to build the “Jacobian” $J$ of the system⁶ and represent it in cuDSS

• stacks multiple versions (exactly batch-size many) of $Y$, $J$ and $v$, giving $Y$, $J$ and $v$. This can be done in an “embarrassingly parallel” manner because only the values of $Y$ and $J$ need to be stacked; the sparsity pattern is the same across the whole batch.

Then the “solve” part will take place, all on the GPU, without any CPU-GPU transfer during computation. For a given number of iterations (user specified, uniform across the elements in the batch):

1. the values of the stacked $J$ are updated from $v$ and $Y$ for all systems in the batch, as if it were a single (block-diagonal) matrix.

2. cuDSS then solves the stacked linear system $J.dx = \text{mis}$ for the whole batch at once

3. the data for each individual system in the batch is read back, which is used to build the updated (stacked) $v$ complex vector at each bus

4. the KCL mismatch is computed for all systems in the batch, giving the vector mis used to restart at step 1.

This way of solving the powerflow with the Newton-Raphson algorithm is beneficial (in terms of computation time) if enough data is used at step 2 above (the most time-consuming step); otherwise, the setup cost might decrease throughput. This can explain why, in Figure 1, the gpusim2grid method does not perform as well as other cpu baselines on the 4 or 5 smallest grids: the GPU might not be saturated in these cases.

With this formulation, this model can already handle, in a single batch:

• varying active generation and demand values

• varying generation voltage setpoint

• varying reactive demand

• disconnection of an arbitrary number of branches (known as $n-k$), including power line or transformer disconnections that would island part of the grid (in this case the islanded part is not solved)

This is why, for the experiments of § 3.4 we could not run the whole dataset in one pass. For now there is no possibility to batch generator disconnection directly (as it would change the “pv” to “pq” attribution of some buses). We had two solutions to handle this:

• Either we kept “everything in one batch” approach, and in this case we had to slightly modify the input data to model the transformer connecting the generator to the bus. This could be done in a way to make the modified model and the new model identical, ensuring the KCL evaluator would still report full convergence for gpusim2grid and then use the gpusim2grid possibility to disconnect these transformers (that would then also disconnect the generator behind them) in the batch.

• or we could more simply run multiple instances of the gpusim2grid simulator, each initialized from the state where the generator status are fixed. We then run all the data in "one batch" for every generator disconnection pattern in the test sets. This simpler approach has been preferred for this paper, though it makes the gpusim2grid much slower in this setting: each instance of gpusim2grid needs to copy $Y$, $s$, build the jacobian etc. which (for small grid like case118) takes a lot of time when compared to the actual powerflow solves.

What could be possible (but is not yet implemented) is to also include, for a given batch, if known before computation, more profound topological modifications (bus splitting/merging, moving an element to a different bus) and changes to bus status (e.g., “pv” to “pq” or conversely). This could be done without changing the core of the architecture and will be part of future developments.

An ability that might also interest current state-of-the-art data-processing pipelines e.g. in [14] is to slightly change the physical parameters of the power lines. This has not been done yet but could be addressed with, in some cases, parameter sampling directly on the GPU for an even higher throughput.

What is not yet addressed is the capacity to batch data from multiple grids of different sizes in the same call, or to modify parameters mid-solve (e.g., disconnecting a line or a generator between iterations 2 and 3).

Physical solvers configuration used

Having all 7 solvers solve exactly the same AC power flow problem has not been an easy task. We had to tweak the default configuration of some of them to make sure every element was represented exactly as in Equation (1).

We think that finding these exact parameters is one contribution of this paper. All the code for the evaluation of all models is provided on GitHub⁷. We encourage anyone interested in the details to check there for the actual implementation. We only expose here a summary of what has been done:

• pandapower: has been one of the most difficult to use because it uses, by default, a different representation for transformers and changes the order of the elements (the $k^{\mathrm{th}}$ bus in the matpower case file is not necessarily the $k^{\mathrm{th}}$ in pandapower). To solve the first issue, we model transformers explicitly as an “impedance” pandapower element (rather than as transformers, as done by default). We also wrote a script to “align” the IDs of the buses, loads, and generators with the matpower data.

• powermodels: actually one of the easiest to use. It is directly “aligned” with matpower data by default. We compute powerflows here with the compute_ac_pf function (rather than the costlier solve_ac_pf) for grids of all sizes (the GENCO paper instead used Ipopt, a costlier solver, for the largest grids, 2k and 10k buses). For these largest grids, like Puech et al. [ [14], we noticed some “unrealistic” solutions, though still marked as valid by PowerModels. We still report computation time for these cases, since they still meet the KCL.

• exapf: also easy to use; the default importer matches the “matpower equations” bit for bit.

• pypowsybl: this CPU solver only imports the matpower “.mat” format (not .m), so we used a converter (via a local matpower installation under the free Octave language) to convert the data. It also required changing the order of the elements (especially buses and generators) by relying on the imported names. We also had to use dedicated powerflow parameters to make pypowsybl run the exact powerflow at the given tolerance: by default, it solves a more complex problem that accounts for generator reactive limits or slack distribution, for example. These more realistic features are out of scope for this benchmark and were switched off.

• lightsim2grid: also among the easiest. Matpower data is read natively and the element order is preserved.

• power grid models (called “pgm” in the scripts, to distinguish it from the similarly named “powermodels”) was also quite challenging. It was primarily developed for distribution grids (DSO), whereas matpower data is primarily transmission-focused (TSO). We wrote a full converter from matpower data to power grid models, partly relying on the pandapower converter (after the pandapower fixes above) and forcing some parameters to specific values (e.g., for the slack, or disabling generator reactive limits). Being designed mainly for the smaller grids typical of DSOs, this model does not work for the largest grids ($\ge2000$ buses); we plan to raise this on the Power Grid Model GitHub.

• rustpower: has no direct interface for reading the matpower format, but accepts lower-level inputs directly ($Y$, $s$, “pv”, “pq”), which we extracted from lightsim2grid as the best-suited source. We also found (and verified in the rustpower code) that we needed to pass the transpose of the $Y$ matrix for results to match matpower’s⁸. We intend to open an issue on their GitHub for this finding.

Solver names and conventions

When trying to benchmark all these solvers, we faced an issue: how to name all the different solvers. Finding a good name for “powermodels” and “power grid models”, or “pandapower” (often abbreviated “pp”) and “pypowsybl” (also abbreviated “pp"), can be confusing. We therefore adopted the following convention:

• pp: pandapower

• powermodels: powermodels

• exapf: exapf, in CPU mode only

• olf: pypowsybl, with its default “open load flow” (abbreviated OLF) powerflow solver

• ls: lightsim2grid

• pgm: power grid models, its standard abbreviation

• rustpower: rustpower

• exapf-gpu: exapf, in GPU mode only

• gpusim2grid: the proposed Newton-Raphson-based model running on GPU

• gns: the “graph neural solver” model [6], with implementation and dimensions taken from the PF-$\Delta$ codebase

• pf-net: the “powerflownet” model [10], with implementation and dimensions taken from the PF-$\Delta$ codebase

• canos: the “canos-pf” model [13], with implementation and dimensions taken from the PF-$\Delta$ codebase

The coloring of the graphs also has a meaning:

• brown to yellow color corresponds to Newton-Raphson-based models running on CPU

• burgundy to red corresponds to Newton-Raphson-based models running on the GPU

• dark blue to light blue corresponds to machine learning-based models

Full results of the KCL mismatch for each tested tool

Table 3 reveals several noteworthy patterns.

First, the parameters found in § 5.4 allow all solvers to solve the exact same problem: all solvers agree within tolerance. The only number $> 1\mathrm{e}{-8}$ here is in the “olf” row (Open Load Flow), for the worst-case mismatch.

Open Load Flow’s (the default solver of pypowsybl) worst-case mismatch is $1.75\mathrm{e}{-06}$ because it does not use the $\lVert\cdot\rVert_{\infty}$ norm internally but the $\lVert\cdot\rVert_{2}$ norm; this is why some values can occasionally exceed $1\mathrm{e}{-8}$, even though the mean mismatch is always well below the tolerance ($4.4\mathrm{e}{-10} \ll1\mathrm{e}{-8}$).

PowerModels’ convergence rate is relatively low compared to the others. This is why, in a real pipeline, the GENCO authors [14] preferred to use the more robust (but slower) mode of PowerModels to generate their data. We instead used the “faster mode” here, since our aim was not to generate data; had we needed to generate real power-flow data, we would have made the same choice as the GENCO authors and switched to the more robust Ipopt-based mode of PowerModels.

PGM (power grid models) does not converge for grids above 500 buses. It was primarily designed to model distribution grids, which have different properties than the transmission grids tested here; this likely explains why it diverges on the largest test cases. We intend to inform the PGM authors about this fact in the coming weeks.

The “convergence rate” of gpusim2grid is “not applicable": as specified in § 5.3, this model’s stopping criterion is not based on accuracy (unlike the other baselines) but on the number of Newton-Raphson iterations (here varying from 0 to 7). As already detailed in the main paper, 5 iterations are enough in this setting to achieve an error lower than $1e - 8$. Though gpusim2grid, like any physical solver, can detect internally the KCL mismatch (from its own representation of the system) and could raise "non-converging flags" if needed.

What is not reported in this table, however, is that sometimes, for some grids, some batch of gpusim2grid reports a “NaN” error. After a thorough investigation, we found that the sparse linear solver sometimes fails silently to find a good solution for the linear system and returns an unrealistically high value while reporting “converged,” causing subsequent iterations to overflow to floating-point “NaN". This has already been reported to NVIDIA. More specific details, with examples of the misbehaving data, will be released alongside the paper’s code for NVIDIA to have more data to work on if interested.

We think that exapf, in GPU mode, suffers from the same type of issue. A more thorough analysis must be carried out, but Figure 2 is consistent with this explanation.

A final note on the “arbitrary” sample size of 9,000 chosen here: it was chosen so that the reported metric (mean mismatch) would be representative enough of the accuracy of the powerflows. We chose 9,000 and then checked with a Monte-Carlo-like estimate that the reported metric had indeed “converged” to its true value. Such graphs were produced for all 15 grids and both metrics; one example can be found in Figure 3.

Deep dive into the loads / generation scenarios

The goal of this experiment was to benchmark the computation time of each solver when pushed to its limit. To that end, we decided to benchmark them on a rather simple task: one grid, many scenarios of demand and generation.

The process we wanted to benchmark is, for a given solver, a given grid, and a given scenario (loaded in main computer memory, RAM):

• update the input of the solver with the new data

• solve the powerflow

As in Puech et al. [ [14] (Figure 7, page 15), we did not want to benchmark the time to load the data from disk or to initialize the solver into a usable state (e.g., loading its weights from the hard drive, compiling it, etc.).

This principle has some consequences:

• In multi-threaded or multi-process experiments, the solver state often has to be initialized separately in each worker. Whenever possible, this initialization time is excluded from the measurements, as it would be negligible in operational settings where it is paid only once and then amortized over many simulations, for example when evaluating power flows under uncertainty in generation and demand.

• As an additional safeguard, we chose workloads so that each core spent about 30 seconds in the solve phase and processed a minimum number of power flows. This amortizes solver initialization or copy overhead over at least 30 seconds of actual computation, making its impact on the reported results negligible. As a result, the fastest solvers on the smallest grids required running several million power flows, likely far more than would arise in real grid operations.

• When dealing with a GPU model, we explicitly wanted our measurements to account for the time to transfer data between the GPU and main memory, in both directions. Indeed, in an operational process, data sitting still on the GPU has close to no use. It should be transferred back to the RAM for use by human operators or automatic processes alike.

• For JIT-compiled models (pandapower for Python, the Julia-based ones, or gpusim2grid) we first ran an untimed evaluation, then started the timing.

• We “warmed up” the GPU for all GPU-based solvers, for example making sure that all the shared libraries were loaded in RAM and that the model was “transferred” to the GPU, etc.

Machine characteristics

We use three machines, each pairing a CPU and a GPU (4 GPUs of the stated model per machine):

• l20 (2×AMD EPYC 9755, 256 cores, ~3 TB RAM, NVIDIA L20 48GB VRAM)

• t4 (2×Intel Xeon Silver 4214R, 24 cores, ~400 GB RAM, NVIDIA Tesla T4 16GB VRAM)

• rtx6k (2×AMD EPYC 9135, 32 cores, ~800 GB RAM, NVIDIA RTX PRO 6000 Blackwell 96GB VRAM).

All machines were on premise (not on the cloud), running RHEL 8.10 (Ootpa). Each machine uses “hyperthreading”, so the available number of logical cores is double the number of physical cores (512 for the “l20” machine, 48 for the “t4” machine, and 64 for the “rtx6k” machine).

On “t4” and “rtx6k”, the Python version used was 3.12, and it was 3.13 on “l20”. The Julia version was 1.12.6 across all machines.

The (main) packages used were the same across all machines:

• pandapower: 3.5.4

• powermodels: 0.21.6 (julia)

• exapf: 0.13.0 (julia)

• pypowsybl / olf: 1.16.1

• lightsim2grid: 1.0.0

• pgm: 1.13.141

• rustpower: 0.5.0

• pytorch: 2.12.1+cu130

• torch geometric: 2.8.0.post1

In principle, all experiments were run on all 3 machines. We only report results for the “rtx6k” machine because we think it’s the “most standard” machine of the 3. “t4” is rather old and approaching end of life, the “l20” has 512 available CPU threads (2 times 256 real physical cores) which is not the most common hardware.

We still wanted to make as many runs on the 3 different machines as possible to ensure that the results we have are robust to different hardware specification.

Highest throughput for CPU solvers

To obtain the highest throughput for the CPU solvers, we ran all 7 solvers using different core counts (from 1 up to 1024 for some experiments) for each grid. If the solver provides a native parallelism feature, we used it; otherwise, we fell back to the native “parallel” tool available in the language (i.e., multiprocessing or multithreading for Python, or native threads in Julia). We ran most of our experiments on 3 machines at our disposal to generate figures similar to Figure 4. In parallel experiments, we disable internal threading in the benchmarked solvers to avoid thread oversubscription. This matters when a method relies on libraries such as OpenMP, BLAS, or MKL for linear algebra. As a safeguard, we set OMP_NUM_THREADS, OPENBLAS_NUM_THREADS, MKL_NUM_THREADS, NUMEXPR_NUM_THREADS, VECLIB_MAXIMUM_THREADS, and RAYON_NUM_THREADS to 1.

Each of these figures is then summarized per grid size to give the maximum throughput per method across all grids, yielding a figure similar to Figure 5.

We now have 3 different plots, one for each machine, presenting the maximum throughput for each of the 7 CPU solvers benchmarked; results are available in Figure 6.

We then report, in the main Figure 1, the lightsim2grid curve, since it is the one with the highest measured throughput for the selected “rtx6k” machine.

Making this figure requires approximately 4700 CPU-hours (≈195 core-days) for the final runs. Most of the time was spent on powermodels (the slowest of the baselines in this setting). 630M power flows have been run for this experiment, most of them by pgm and lightsim2grid on smallest grids: the fastest solver in this benchmark.

Caveats We want to emphasize some limits of this benchmark methodology. We did not try to optimize each solver for speed so as to achieve its absolute best performance, which would have been mandatory if we wanted to rank solvers based on their throughput, which is not our goal here. Our approach was centered around a “best effort” strategy: solvers were optimized for speed using easily available/documented features, as we believe is common practice in the “machine learning for power systems” community.

Some solvers (pgm, lightsim2grid) natively offer features matching almost 100% of our benchmark here (scenarios where loads and generation change but not the topology). We decided to use this capability for these two solvers. Our benchmark does not try to rank these CPU solvers, but rather to show what is possible with a CPU; this is why we used the native feature whenever available.

Conversely, some other solvers (e.g., pandapower, powermodels, pypowsybl/OLF) “start from scratch” each time a new scenario is computed, for example by rebuilding the entire $Y$ matrix and $s$ vector. This might be rather inefficient, but this absolutely does not mean these solvers are slow in general; it only means they have not been optimized for the specific task used in this benchmark. They could shine on other types of tasks not benchmarked here, for example contingency analysis.

We were more familiar with Python than with Julia. Following our “best effort” principle, we did not try to make our Python code overly complex, though it might naturally be more efficient simply because of our familiarity with the language. This would have been an issue had we tried to rank CPU solvers, which has never been the goal of this benchmark.

Highest throughput for GPU / ML solvers

We adopted a similar methodology to find the best throughput of the methods using GPUs (Newton-Raphson-based solvers and machine-learning surrogates), but instead of varying the number of CPU cores, we varied the batch size.

For learned models, we have not tried to retrain them. We only used published architectures and their hyperparameters (number of “layers”, size of each “layer”, etc.). We are only concerned with evaluation throughput and this is why we also deactivated the computation of gradients for these models.

Usually, we generated a dataset consisting of 16,384 rows (for each grid size). Then, on this dataset, we ran the same model for batch sizes 1, 2, 4, ..., 16384. To save computation time, we added two break conditions: a) if there are more than 100 batches, only the first 100 batches are simulated and the average is reported over these 100 batches; b) if a model reports an “Out of Memory” error on the device for a given batch size, the larger batches are not simulated.

By diving into the computation time, we also found some optimizations that can be made to PF-Δ’s useful InMemoryDataset, which we will share with the community, notably regarding faster preprocessing in the coming weeks.

First, we produced figures for a given method on a given grid, as shown in Figure 7.

Then we aggregate the results for a given method across all computed grids:

We now have 3 different plots, one for each machine, presenting the maximum throughput for each of the 5 GPU solvers benchmarked (2 Newton-Raphson and 3 machine-learning); results are available in Figure 9.

We then pick the best performing method in each category from this plot to include it in the Figure 1 of the main paper: Graph Neural Solver is the fastest learned model and gpusim2grid is the fastest GPU solver available.

In Figure 9 we can see the phenomenon described in § 5.6 regarding exapf (and most likely cuDSS instability). As shown in Figure 2, for some grids (like case14, case_ieee30, and case57) cuDSS fails to solve the systems. This limits the batch size to only 256 (larger sizes diverge for these grids, for reasons that remain unclear), which is not enough to saturate the GPU on these small grids, leading to poor throughput for exapf on the 3 smallest grids.

We can also notice a small artifact of our benchmark: for the smallest cases, too few scenarios were simulated for gpusim2grid to saturate the GPU. As a result, the curve is flat at the beginning of the plot, whereas it should decrease linearly on this log-log scale.