Low-rank corrections to a shared reference Jacobian for GPU-parallel AC $N\!-\!1$ contingency analysis

1

Low-rank corrections to a shared reference Jacobian for GPU-parallel AC N −1 contingency analysis

arXiv:2611.03012v1 [eess.SY] 2 Nov 2026

Marc Hunkemöller , Nico Westerbeck , Joost van Dijk , Lars Schewe

and Dirk Witthaut

Abstract—Contingency analysis repeatedly solves the powerflow equations under many single-element outages, an increasingly demanding workload as operating margins tighten and renewable variability grows. Graphics processing units (GPUs) are attractive for this task, but a change in topology alters the sparsity pattern of the Jacobian, preventing the uniform batched factorization on which efficient GPU solvers rely. This paper presents a GPU-parallel Newton–Raphson solver for N−1 contingency analysis that separates the change in operating state from the change in topology. A single reference Jacobian is inverted once and reused across all contingencies. Each topology change is expressed as a low-rank Woodbury correction, while the residual state-dependent perturbation is handled by an inner iterative solve, a truncated Neumann series or the biconjugate gradient stabilized method. We analyze the spectral radius and condition number of the resulting iteration matrices and show that the reference Jacobian is an effective preconditioner. On transmission-scale benchmarks, the method attains per-topology solve times competitive with an optimized 80-core central processing unit baseline, with the low-rank correction accounting for approximately 1–2% of the workload. This makes GPU-parallel N −1 contingency analysis feasible at transmission scale without per-topology refactorization.

For accurate AC analysis, Newton–Raphson (NR) methods combined with sparse linear algebra remain the standard, owing to their robustness and fast local convergence [2], [3], [4]. In many operational settings, however, the computational bottleneck is not a single NR solve but the repeated solution of many closely related PF problems, most prominently N − 1 contingency analysis, in which every single-element outage of a system with N components must be assessed. Because the number of contingencies grows with system size, this workload, combined with the high similarity between contingency cases (only one element changes), is a natural target for parallel hardware such as graphics processing units (GPUs). However, the obstacle is structural. Varying nodal injections preserves the Jacobian’s sparsity pattern, so a single symbolic factorization can be reused. A topology change, however, alters this pattern. The resulting fill-in and irregular memory access make sparse direct solvers a poor fit for the single-instruction, multiple-data (SIMD) parallelism of GPUs, limiting their application across topologies.

Index Terms—Contingency analysis, GPU computing, low-rank update, Newton-Raphson method, parallel algorithms, power flow analysis

In this work we develop an AC PF solver that resolves this obstacle by separating the Jacobian’s dependence on the system state from its dependence on the topology. Each local topology modification is represented as a structured low-rank correction to a fixed-sparsity reference Jacobian, whose inverse is computed once and reused across all contingencies. The low-rank correction is applied at solve time through the Woodbury matrix identity, thus avoiding a refactorization, while the remaining state-dependent perturbation is resolved by an iterative inner solve, either a truncated Neumann series or biconjugate gradient stabilized (Bi-CGSTAB) method. Then, every contingency reduces to batched dense and sparse matrix products with shared sparsity pattern. Thus, the entire N −1 study maps efficiently onto GPU hardware while retaining the accuracy and robustness of a full AC NR solve.

I. I NTRODUCTION AST and reliable power-flow (PF) computation is essential for power-system operation and planning. It is used in state estimation, contingency screening and security assessment and appears as the repeated inner solve in optimization tasks such as optimal power flow and topology optimization [2], [3], [4]. These demands intensify with the energy transition: growing shares of volatile, converter-based generation, tighter operating margins from demand growth outpacing transmission expansion and more frequent corrective actions require evaluating large numbers of PF scenarios under time constraints [5], [6].

F

M.H., D.W. and L.S. acknowledge financial support from the Elia Group via the Elia Group Open Research Challenge 2025. The authors gratefully acknowledge computing time on the supercomputer JURECA [1] at Forschungszentrum Jülich under grant no. gpu-to-power. M. Hunkemöller and D. Witthaut are with the Institute of Climate and Energy Systems (ICE-1), Forschungszentrum Jülich, 52428 Jülich, Germany, and also with the Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany (e-mail: m.hunkemoeller@fz-juelich.de; d.witthaut@fz-juelich.de). L. Schewe is with the School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Edinburgh EH9 3FD, U.K. (e-mail: lars.schewe@ed.ac.uk). J. van Dijk is with TenneT TSO BV, 6812 AR Arnhem, The Netherlands (e-mail: joost.van.dijk@tennet.eu). N. Westerbeck is with Elia Group, 1000 Brussels, Belgium (e-mail: nico.westerbeck@50hertz.com).

The remainder of this paper is organized as follows. Section II reviews related work on linearized screening, lowrank modification techniques, GPU acceleration and iterative solvers for PF. Section III recalls the AC PF equations, the low-rank structure of the Jacobian under topology changes and the numerical tools we build on. Sections IV and V present the Jacobian decompositions and the resulting GPU-parallel algorithms. Section VI analyzes the spectral properties that govern their convergence and selects the best variant, which Section VII benchmarks against multi-core central processing unit (CPU) baselines on transmission-scale networks. Section VIII concludes.

2

II. L ITERATURE REVIEW A. Linearized and screening methods For computational efficiency, linearized models such as DC PF and its refinements are widely used [7], [8], with extensions incorporating losses, reactive power and voltage magnitudes [9], [10], [11]. Distribution-factor methods in particular enable extremely fast contingency and topology screening [12], [13], [14], because a topology change enters as a low-rank update to a fixed linear operator and recent formulations extend these factors to capture voltage sensitivity [15]. However, the linear approximation can break down in stressed operating regimes where the nonlinear AC effects in the contingency assessment are no longer negligible. B. Low-rank modifications and factorization updates The observation that a single-element outage is a low-rank change to the system matrix has a long history in power engineering, dating back to compensation-based techniques for handling network modifications without refactorizing [16]. In numerical linear algebra, efficient algorithms for low-rank updates of sparse LU factorizations are well established [17], [18]. Direct application of such factor-update schemes on GPUs is nontrivial, because the update algorithms pivot across rows and columns, requiring data movement in GPU memory that undermines throughput. And, for contingency analysis, the full updated factorization is not even needed, only the solution of the corresponding linear system. This motivates applying the low-rank correction at solve time rather than updating the factors, for which the Woodbury matrix identity provides the natural tool [19], [20]. In the DC approximation, low-rank updates of the power transfer distribution factors can represent line outages and substation reconfigurations without refactorization, achieving extreme throughput on GPUs [21]. This demonstrates the potential of the low-rank viewpoint for AC PF with topology changes. C. GPU acceleration for power flow GPUs have been widely used to accelerate the many PF solves required when numerous operating scenarios share a single network topology. The typical case is probabilistic power flow, where a Monte-Carlo ensemble of nodal injections, e.g., to reflect the uncertainty of renewable generation and load, is evaluated on a fixed grid [22], [23]. Time-series simulation and scenario-based planning studies fall into the same category. In all of these studies, only the injections vary while the topology and hence the Jacobian sparsity pattern, stays fixed. The workload therefore reduces to a batch of structurally identical linear solves, which maps naturally onto the GPU through batched sparse factorization or batched matrix operations [24], [25]. The situation is fundamentally different for N − 1 contingency analysis and topology optimization, where each scenario removes or modifies a network element. Such a change alters the sparsity pattern of the Jacobian, so the fixed-pattern batching of the multi-scenario case no longer applies directly

and sparse direct solvers become sensitive to the resulting fillin and irregular memory access. This is addressed by regenerating and refactorizing many Jacobians in parallel on the GPU for static security analysis [26], which recovers parallelism but repeats the factorization for every topology. Compensationtype methods treat a branch outage as a low-rank update to the constant fast-decoupled matrices, reusing a single factorization across all contingencies on parallel hardware [27]. We instead separate the state and topology dependence of the full AC Jacobian and apply the low-rank correction at solve time.

D. Iterative solvers for power flow Iterative solvers are a long-standing alternative to direct factorization for the linear systems arising in each NR PF solve. Krylov-subspace methods were introduced to PF computation as a scalable alternative to LU factorization [28]. Subsequent work combined them with NR in inexact-Newton and recycling schemes, comparing inner solvers such as the generalized minimal residual (GMRES), Bi-CGSTAB and Quasi-Minimal Residual (QMR) methods and reusing preconditioner or subspace information across NR iterations [29]. Fast-decoupled and conjugate-gradient formulations that retain a largely fixed sparsity pattern have likewise been proposed, avoiding the refactorization that a topology-induced change in matrix structure would otherwise require [30], [31]. The appeal of these methods on parallel hardware is that the dominant operation is the matrix-vector product, which parallelizes far more readily than the inherently sequential triangular solves of a direct method.

E. Fair benchmarking of GPU speedups Reported GPU speedups for PF computation span nearly two orders of magnitude, largely determined by the choice of CPU baseline. High speedups are often reported relative to general-purpose open-source tools [23], while comparisons against optimized sparse solvers can yield more modest figures [32]. This mirrors the more general finding that, once both sides are properly optimized and parallelized, the GPU advantage shrinks substantially relative to the orders of magnitude frequently claimed [33]. We therefore benchmark against multi-core CPU baselines rather than a single-threaded reference seeking a CPU baseline of comparable cost and capability rather than an artificially weak one.

III. BACKGROUND A. Power-flow equations and Newton–Raphson method The AC PF equations describe the steady state of a power grid. They relate the real and reactive power injections Pi and Qi at node i to the voltage magnitudes and phase angles at all nodes. We work with logarithmic voltage variables ui = log vi , which leads to a compact Jacobian structure [34]. The phase angles are denoted as θi . With the shorthand θij =

3

θi − θj and the nodal admittance matrix entries Gij + jBij , the PF equations read X Pi = eui euj (Gij cos θij + Bij sin θij ) =: pi (x), (1) j

Qi =

X

eui euj (Gij sin θij − Bij cos θij ) =: qi (x).

(2)

j

We distinguish n PQ nodes (loads) with fixed (Pi , Qi ), m PV nodes (generators) with fixed Pi and voltage magnitude ui and one slack node with fixed voltage magnitude and θ = 0. Indexing PQ nodes by i = 1, . . . , n, PV nodes by i = n + 1, . . . , n + m and the slack node by n + m + 1, the network has N = n + m + 1 nodes and ℓ branches. Collecting the free state variables as x = (θ1 , . . . , θn+m , u1 , . . . , un )⊤ ∈ RNm with Nm = 2n + m and the corresponding power mismatches as f (x) = (p1 (x) − P1 , . . . , q1 (x) − Q1 , . . .)⊤ , the PF problem takes the compact form f (x) = 0. The NR method solves this system iteratively. Expanding f to first order about the current iterate x(σ) yields the linear update equations J (x(σ) ) ∆x(σ) = −f (x(σ) ), (σ+1)

x

(σ)

=x

+ ∆x

(3) (σ)

,

(4)

where σ is the iteration index and J ∈ RNm ×Nm is the Jacobian. To expose the structure that the low-rank correction exploits, we fix an orientation for each of the ℓ branches and define the signed node-edge incidence matrix E ∈ RN×ℓ , with Ei,a = +1 (−1) if branch a originates (terminates) at node i and 0 otherwise, together with its unsigned counterpart Di,a = |Ei,a |. Throughout, a tilde (˜·) denotes omission of the slack entry and a hat (ˆ·) omission of the slack and PV entries, applied consistently to matrices (rows) and vectors (components). In terms of these operators we define the perbranch quantities c = cos(E ⊤ θ) ⊙ exp(D ⊤ u), s = sin(E ⊤ θ) ⊙ exp(D ⊤ u), where ⊙ denotes the element-wise multiplication. With b, g ∈ Rℓ collecting the series susceptance and conductance of each branch, ordered consistently with the edge list, the four Jacobian blocks take the compact form ⊤

⊤

⊤

⊤

⊤

⊤

⊤

⊤

J pθ = Ẽ[b ⊙ c]Ẽ − D̃[g ⊙ s]Ẽ , J qθ = D̂[b ⊙ s]Ẽ + Ê[g ⊙ c]Ẽ , J pu = Ẽ[b ⊙ s]D̂ + D̃[g ⊙ c]D̂ + [2γ̃ ⊙ exp(2ũ)], J qu = Ê[g ⊙ s]D̂ − D̂[b ⊙ c]D̂ − [2β̂ ⊙ exp(2û)], (5)   J J pu with the full Jacobian being J = J pθ . Here [·] denotes qθ J qu the diagonal matrix formed from a vector and β, γ collect the diagonal susceptance and conductance entries of the bus admittance matrix. B. Low-rank Jacobian structure under topology changes We want to solve (3) in parallel for a set of grid topologies T , for instance all single-line outage scenarios in an N − 1

security assessment. The central obstacle is that topology changes alter the sparsity pattern of J , which prevents efficient batched parallel factorization with a shared symbolic structure. The key observation is that the change is of low rank. We fix a reference topology r and consider a contingency topology t ∈ T that differs from r by the outage of a single branch a = (i, j). The series susceptance and conductance of branch a change by ∆ba and ∆ga , respectively. Similarly, the diagonal susceptance and conductance entries β, γ at the terminal nodes i, j change by ∆βi , ∆βj and ∆γi , ∆γj . The resulting change in the Jacobian has a rank-four structure, J t (x) − J r (x) = Rt K t (x) R⊤ t ,

(6)

where Rt ∈ R(2n+m)×4 is a sparse matrix defined with the standard basis vectors r a = ei and sa = ej as,   r a (1 : n+m) sa (1 : n+m) 0 0 Rt = , 0 0 r a (1 : n) sa (1 : n) with r a (1 : k) denoting the projection onto the first k components. And K t (x) ∈ R4×4 contains the admittance changes and depends on the state, as defined in (7). For more details on the low-rank-correction structure of the Jacobian matrix see [35]. C. Mathematical tools Three tools are central to the algorithm developed in this paper: the Woodbury matrix identity and iterative solvers based on the Neumann series and Bi-CGSTAB. 1) Woodbury matrix identity: For conformable matrices M , U , C and V , the inverse of a low-rank-corrected matrix satisfies (M + U CV )−1
−1 = M −1 − M −1 U C −1 + V M −1 U V M −1
−1 = M −1 − M −1 U I + CV M −1 U CV M −1 , (8) provided the relevant inverses exist [19]. Here, I denotes the identity matrix. This avoids inverting the full corrected matrix: if M −1 is known, the only remaining inversion is of the small matrix C −1 + V M −1 U , whose size equals the rank of the correction and is therefore cheap. In the present context, M = J r (xr ) and the correction has rank four, so the inversion reduces to a 4 × 4 system. We apply J r (xr )−1 by storing the dense inverse explicitly, computed once at the start of the algorithm. This one-time inversion carries the dominant computational cost. The subsequent application to any right-hand side is a dense matrix– 2 vector (or matrix–matrix) product in O(Nm ) operations per column, which maps efficiently onto the GPU. This choice runs counter to common practice on CPUs, where an explicit dense inverse is avoided in favor of a sparse LU factorization and triangular solves. On the GPU, however, the regular, branch-free dense products underlying the explicit dense inverse achieve far higher throughput than the irregular, sequential memory access of sparse triangular solves, particularly when many right-hand sides are processed simultaneously. The dense inverse also carries a cost: it is

4

 +∆ba ca − ∆ga sa −∆ba ca − ∆ga sa K= +∆ba sa + ∆ga ca +∆ba sa − ∆ga ca

−∆ba ca + ∆ga sa +∆ba ca + ∆ga sa −∆ba sa − ∆ga ca −∆ba sa + ∆ga ca

+∆ba sa + ∆ga ca + 2∆γi e2ui −∆ba sa + ∆ga ca +∆ga sa − ∆ba ca − 2∆βi e2ui −∆ga sa − ∆ba ca

2 no longer sparse, with a memory footprint of O(Nm ), which remains modest for the systems considered. Because the inverse will serve only as a preconditioner, it does not need to be stored in full double precision. Single precision halves its memory footprint and substantially reduces the cost of the dense products. In addition, an explicit inverse is less numerically stable than a factorization [20], more so at reduced precision and as the condition number grows with system size, thus some accuracy is lost in the preconditioner application. As it only preconditions an inner iterative solve rather than providing the final solution, this loss is tolerable. 2) Neumann series: The inverse of a matrix (I − A) can be expanded as a power series,

(I − A)−1 =

∞ X

Ak ,

(9)

k=0

provided the spectral radius satisfies ρ(A) < 1 [36]. The series converges geometrically with asymptotic factor ρ(A), so it converges rapidly when ρ ≪ 1. Thus it can be truncated after Nk terms with a truncation error that decays asymptotically as ρ(A)Nk [37], [38]. 3) Biconjugate gradient stabilized method: The biconjugate gradient stabilized (Bi-CGSTAB) method [39] solves a linear system By = b iteratively and, as a Krylov subspace method, applies to any nonsingular B without requiring a spectral radius condition. We use Bi-CGSTAB rather than the conjugate gradient (CG) method because the PF Jacobian is nonsymmetric and CG is restricted to symmetric positive-definite systems. For a symmetric positive-definite system, the convergence of CG is rigorously bounded in terms √ of the condition number √ κ alone, decreasing by a factor ( κ − 1) / ( κ + 1) per iteration [38]. For a nonsymmetric operator no such clear bound exists: the convergence of Bi-CGSTAB depends on the full spectrum [40], [41] and the condition number can only give a coarse estimate. The efficiency of Bi-CGSTAB depends strongly on preconditioning [38]. Common general-purpose preconditioners such as diagonal (Jacobi) scaling or incomplete LU factorizations [38] would require recomputation for each topology. Rather than a generic factorization-based preconditioner, we exploit the structure of the contingency problem and use a preconditioner P built from a reference system, whose sparsity pattern is shared across all topologies. Explicit forms of the preconditioner P are detailed in Sec. V and VI. The preconditioned operator is B ′ = P −1 B, and the preconditioner is effective when B ′ is close to the identity, that is, when the deviation A := I −B ′ is small in norm. The condition number κ(B ′ ) serves as a convenient scalar measure of this closeness, with κ → 1 indicating a near-identity operator.

 +∆ba sa + ∆ga ca −∆ba sa + ∆ga ca + 2∆γj e2uj    +∆ga sa − ∆ba ca −∆ga sa − ∆ba ca − 2∆βj e2uj

(7)

IV. JACOBIAN DECOMPOSITIONS The key difficulty in parallelizing the NR solve across many topologies on a GPU is that the Jacobian J t (xt ) depends on both the topology t and the state xt , which vary across scenarios. To enable efficient parallelization, we decompose J t (xt ) into a reference part that is shared across all topologies, a low-rank correction that is included with the Woodbury matrix identity and a perturbation that is handled by an iterative solver. The reference part has a fixed sparsity pattern for all t, which is the prerequisite for efficient parallelization on a GPU. We fix a reference topology r and a reference state xr , i.e., the base-case N −0 solution, and rewrite the Jacobian for topology t ∈ T as (α)

(α) J t (xt ) = J r (xr ) + Rt K t (xt ) R⊤ t + ∆J r (xt ), (10) | {z } reference

where Rt ∈ R(2n+m)×4 is the sparse topology-change matrix, (α) K t ∈ R4×4 is a small dense matrix, both similar to the definition in Sec. III-B, and ∆J (α) r (xt ) is a perturbation. As described in the following, we can make different choices for (α) the exact definition of K t and ∆J (α) r (xt ). Formulation α = 1 (baseline): The entire difference from the reference is absorbed into the perturbation: (1)

Kt

= 0,

∆J (1) r (xt ) = J t (xt ) − J r (xr ).

(11)

No low-rank correction is applied. The perturbation contains both the topology change and the state deviation and thus is generally large. Formulation α = 2 (physical method): The low-rank term captures the topology change r → t at the contingency state xt and the perturbation retains only the state-induced change in the reference topology: J t (xt ) = J r (xr ) + J t (xt ) − J r (xt ) + J r (xt ) − J r (xr ) . | {z } | {z } (2)

= Rt K t (xt ) R⊤ t

(2)

= ∆J r (xt )

(12) The perturbation ∆J (2) has the same sparsity pattern as J r r and is small whenever the contingency state xt is close to the reference state xr . Formulation α = 3 (physical method with terminal correction): This is an extension of α = 2. The entries of ∆J (2) corresponding to the terminal nodes of the failing r branch, where the impact of the outage is largest, are moved completely into the low-rank term. This yields a modified (3) K t and a perturbation ∆J (3) with smaller entries than r ∆J (2) . In principle this alters the sparsity pattern of the r perturbation. We retain the vacated positions as explicit zeros to keep the pattern fixed across all contingencies. Since the

5

Algorithm 1 Woodbury-first, single NR step using Neumann approximation, all t ∈ T in parallel. For the Bi-CGSTAB variant replace the Neumann inner loop (steps 5–8) by a preconditioned Bi-CGSTAB solve. Require: reference inverse J r (xr )−1 (precomputed once) 1: for t ∈ T in parallel do ▷ parallel over topologies 2: compute residual f t (xt ), perturbation ∆J (α) r (xt ) (α) and low-rank factors Rt , K t 3: S t ← J −1 ▷ only once r Rt 4: z 0 ← −f t (xt ), ∆xt ← 0 5: for k = 0, 1, . . . , Nk − 1 do ▷ Neumann series 6: y k ← P −1 ▷ via Woodbury t zk 7: ∆xt ← ∆xt + y k 8: z k+1 ← − ∆J (α) r (xt ) y k

matrix enters only through matrix products, the stored zeros carry no cost beyond storage. V. PARALLEL SOLVING ALGORITHM Given the decomposition (10), the NR update ∆xt is determined by h i (α) (α) J r (xr ) + Rt K t R⊤ ∆xt = −f t (xt ). (13) t + ∆J r The three summands suggest two natural strategies for organizing the solve. One may first apply the Woodbury matrix iden(α) tity to invert J r (xr ) + Rt K t (xt )R⊤ t exactly and then treat (α) ∆J r (xt ) as a perturbation handled by an iterative solver. We refer to this as the Woodbury-first ordering. Alternatively, one may first handle ∆J (α) r (xt ) iteratively using J r (xr ) as a preconditioner and then apply the Woodbury correction. We refer to this as the iterative-first ordering. The two orderings differ in which part of the operator is inverted exactly and which is left to the iterative solver and hence in the iteration matrix that governs convergence. We derive both orderings below and analyze their spectral properties in Sec. VI. In each ordering, the reference inverse J r (xr )−1 is computed once and reused across all topologies and NR iterations, so that the per-topology factorization of conventional PF is replaced by this shared inversion, a low-rank Woodbury correction, and an iterative solve. A. Woodbury-first For the Woodbury-first ordering, the reference Jacobian and the low-rank correction are grouped into the matrix that (α) is inverted exactly, P t = J r (xr ) + Rt K t (xt )R⊤ t . Left−1 multiplying the NR update equation by P t gives h i (α) I + P −1 ∆J (x ) ∆xt = −P −1 (14) t t r t f t (xt ), which is solved with the Neumann series approximation ! NX k −1  k −1 (α) ∆xt = − −P t ∆J r (xt ) P −1 t f t (xt ). (15)

Algorithm 2 Iterative-first, single NR step using Neumann approximation, all t ∈ T in parallel. For the Bi-CGSTAB variant replace the Neumann inner loop (steps 5–8) by a preconditioned Bi-CGSTAB solve. Require: reference inverse J −1 r (xr ) (precomputed once) 1: for t ∈ T in parallel do ▷ parallel over topologies 2: compute residual f t (xt ), perturbation ∆J (α) r (xt ) (α) and low-rank factors R , K t t   3: stack RHS Z 0 ← −f t (xt ) Rt ▷ 5 columns 4: G←0 ▷ accumulated Ξ−1 t [ −f t | Rt ] 5: for k = 0, 1, . . . , Nk − 1 do ▷ Neumann series 6: Y k ← J −1 ▷ dense GEMM r Zk 7: G←G+Yk 8: Z k+1 ← − ∆J (α) ▷ sparse SpMM r (xt ) Y k 9: g t ← G[:,1] , W t ← G[:,2:5] ▷ split columns   (α) (α) 10: solve 4 × 4: I + K t R⊤ W z = K t R⊤ t t t t gt 11: ∆xt ← g t − W t z t ▷ Woodbury combine of P −1 to a vector v is computed via the Woodbury matrix t identity, −1 v P −1 t v = J r (xr ) h i−1 (α) ⊤ (α) −1 − S t I + K t Rt S t K t R⊤ v, t J r (xr )

where S t = J r (xr )−1 Rt . The only remaining solve is the 4× (α) 4 system in I + K t R⊤ t S t , negligible compared to the other steps. The complete Woodbury-first algorithm for a single NR step is given in Algorithm 1. Steps 5–8 realize the Neumann series, with the inner application of P −1 in step 6 evaluated through the Woodbury t identity (16). The Bi-CGSTAB variant instead solves the via Eq. (16). preconditioned system (14), again applying P −1 t B. Iterative-first Now, instead we only use the reference Jacobian J r (xr ) as a preconditioner and group the reference Jacobian and the perturbation together, Ξt = J r (xr ) + ∆J (α) r (xt ), so that the update equation reads h i (α) Ξt + Rt K t (xt )R⊤ t ∆xt = −f t (xt ).

(17)

(18)

The Woodbury identity is now applied using the full inverse of Ξt , −1 ∆xt = −Ξ−1 t f t (xt ) + Ξt Rt h i−1 (19) (α) (α) −1 −1 × I + K t R⊤ K t R⊤ t Ξt Rt t Ξt f t (xt ). −1 In practice, the terms Ξ−1 t f t (xt ) and Ξt Rt are not computed separately but stacked column-wise into a single righthand side Z t = [ −f t (xt ) | Rt ] and solved together by a single batched Neumann iteration,  −1 (α) Ξ−1 Zt t Z t = J r (xr ) + ∆J r (xt )

k=0

Similarly, Eq. (14) can be viewed as preconditioned system and solved with Bi-CGSTAB. In both cases, each application

(16)

=

NX k −1 

k −1 −(J r (xr ))−1 ∆J (α) (J r (xr )) Z t , r (xt )

k=0

6

α=3

ρ(A )

100

I

101

101 I

W)

10−2

100

α=3

100

10−1

10−1

102

102 I

ρ(A )

W)

α=2

10−1

ρ(A )

W)

α=2

10−2

ρ(A

case118

α=1

10−1

ρ(A

case9

ρ(A

case6ww

100

100 10−2

1

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

100 10−2

1

(α)

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

(α)

Fig. 1. Spectral radius ρ of the iteration matrix, AW = −P −1 and AI = −J r (xr )−1 ∆J r , for decompositions α = 1, 2, 3 shown for all t ∆J r N −1 contingencies and Newton–Raphson iterations for test cases case6ww, case9 and case118 (rows). Left three columns: Woodbury-first ordering. Right two columns: iterative-first ordering. Since α = 1 yields identical results for both orderings, it is shown only in the left block. The dashed line marks ρ = 1, above which the Neumann series diverges. α=1

α=2

α=3

α=2

I

3 2

3 2

I

κ(I − A )

1

101

101

100

100

103

103 κ(I − A )

102

I

W)

κ(I − A W)

κ(I − A

case9

1

case118

α=3

4 κ(I − A )

W)

κ(I − A

case6ww

4

101 100

1

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

102 101 100

1

2 3 4 5 NR Iteration

6

1

2 3 4 5 NR Iteration

6

Fig. 2. Condition number κ of the preconditioned system I − A, with AW and AI as defined in (21) and (22), for the same combinations and panel ordering as Fig. 1.

so that each Neumann term reduces to a dense product with the shared inverse J r (xr )−1 and a sparse product with ∆J (α) r (xt ), both applied to all columns of Z t at once. This batching, together with the parallelization over topologies t ∈ T , is what maps the solve efficiently onto the GPU.

scribed above and steps 9–11 perform the rank-four Woodbury step. The Bi-CGSTAB variant replaces the inner loop (steps 5–8) by a batched preconditioned Bi-CGSTAB solve with the preconditioner J r (xr ).

(1)

Note that with a hot-start xt = xr , the perturbation ∆J (α) r (xt ) vanishes at the first NR iteration for α = 2, 3. Thus, for both orderings the first NR step reduces to a single exact Woodbury correction with no inner iterations. The complete iterative-first algorithm is given in Algorithm 2. Steps 5–8 realize the batched Neumann iteration de-

VI. S PECTRAL RADIUS AND CONDITION NUMBER Both orderings reduce the Newton step to a preconditioned system of the same structure but with a different preconditioner P , which we keep generic in the following derivation and specialize to each ordering afterward. Writing the iteratively

7

solved part of the system as P ∆xt + ∆J (α) r (xt ) ∆xt = −f t (xt ) and left-multiplying by P −1 gives

where A is the iteration matrix whose spectral properties govern convergence. The two orderings correspond to two choices of P . In the Woodbury-first ordering the precondi(α) tioner P t = J r (xr ) + Rt K t (xt )R⊤ t is applied exactly via the Woodbury identity (8), giving the iteration matrix (α) AW = −P −1 t ∆J r (xt ),

(21)

(α) so that I − AW = I + P −1 t ∆J r (xt ) recovers the preconditioned system of Eq. (14). In the iterative-first ordering the reference Jacobian J r (xr ) alone serves as preconditioner, giving the iteration matrix

AI = −J r (xr )−1 ∆J (α) r (xt ),

ρ(A )

Neumann

I

10−2 t

A := −P −1 ∆J (α) r (xt ), (20)

Residual ‖f ‖∞

(I − A) ∆xt = −P −1 f t ,

101

0.08

2.0

0.06

10−8

0.04

10−11

0.02

0

1

2 3 4 5 NR iteration

6

I

2.2

10−5

10−14

κ(1 − A )

Bi-CGSTAB

0.10

1.8 1.6 1.4 1.2

0

1

2 3 4 5 NR iteration

6

Fig. 3. Residual ∥f t ∥∞ versus NR iteration for the Neumann (left) and Bi-CGSTAB (right) inner solves on case118, at inner-iteration count Nk = 5. Each line is one N−1 contingency, colored by the spectral radius ρ(AI ) (left) and condition number κ(I − AI ) (right) of its iteration matrix. The dashed line marks the convergence threshold ∥f t ∥∞ < 10−9 . Residuals plateau near machine precision (∼ 10−14 ).

(22)

and leaving the low-rank term to a subsequent Woodbury correction. The efficiency of both iterative solvers is governed by the spectral properties of A: the Neumann series requires ρ(A) < 1, while for the preconditioned Bi-CGSTAB the preconditioner quality is studied through the condition number κ(I − A) of the preconditioned system, keeping in mind that this gives only a coarse estimate of the convergence. We therefore study ρ(A) and κ(I − A) for all three decompositions α = 1, 2, 3 combined with both orderings, on the standard benchmark cases case6ww, case9 and case118, imported via pandapower [42], [43]. This spectral analysis is restricted to small systems because it requires computing eigenvalues and singular values of the iteration matrix for every contingency and NR iteration, whose costs scale as 3 ), which becomes intractable for large grids. To characO(Nm terize the iteration matrix independently of any particular inner iterative solver, the update step itself is computed with an exact linear solve of the full contingency Jacobian J t (xt ), then the spectral quantities are evaluated at each true NR iterate. The spectral radii across all N −1 contingencies and NR iterations are shown in Fig. 1 and the corresponding condition numbers in Fig. 2. For α = 1 we only show one ordering since both orderings coincide. For α = 2, 3 the perturbation ∆J (α) r (xt ) vanishes at the first NR iteration (Sec. V) due to the hot start with xr , so no box appears at the first iteration. The spectral radii for α = 1 cluster just below unity across all NR iterations, confirming that the baseline decomposition is unsuitable for a fast converging Neumann solve. Both α = 2 and α = 3 reduce ρ substantially, with α = 3 achieving the lowest values across both orderings. The iterative-first ordering yields systematically lower spectral radii than the Woodbury-first ordering for the same α. A small number of outlier contingencies with ρ > 1 persist for the Woodbury-first ordering, particularly in case118, corresponding to cases where the state deviation introduced due to the topology change is large and better handled after the iterative inner solver. The condition numbers exhibit the same qualitative behavior as the spectral radii: α = 1 yields poorly conditioned systems, while α = 2 and α = 3 lower the bulk of the condition-number

distribution, with α = 3 giving the strongest reduction, while, again, the iterative-first ordering outperforms Woodbury-first for a given α. The lowest condition numbers are achieved by α = 3 with the iterative-first ordering. This agreement across the two metrics is central to the results of this analysis. A small ρ ensures rapid convergence of the Neumann series, while a small κ reflects a wellconditioned system for the Bi-CGSTAB iteration. The fact that both are minimized by the same configuration means a single choice optimizes both iterative solvers and thus constitutes the best preconditioning strategy. We therefore adopt α = 3 with the iterative-first ordering for all subsequent numerical experiments. To confirm that these spectral quantities translate into real achieved convergence, Fig. 3 shows the residual trajectories of the actual Neumann and Bi-CGSTAB inner solves on case118, at the same inner-iteration count Nk = 5 adopted for the performance evaluation in Sec. VII. The spectral quantities are evaluated along the exact NR trajectory (cf. Fig. 1 and Fig. 2), whereas the residuals are those of the actual inner solver. The trajectories are colored by ρ(AI ) (Neumann) and κ(I − AI ) (Bi-CGSTAB). Both patterns confirm the spectral analysis: the contingencies with the largest spectral radii (Neumann) and the largest condition numbers (Bi-CGSTAB) converge visibly slower, while the remainder reach near machine precision (residual ∼ 10−14 ) within fewer iterations. The rate separation is most pronounced for the Neumann series, whose geometric convergence factor ρ directly sets the per-iteration reduction. For Bi-CGSTAB, where no single scalar fully characterizes the convergence (see Sec. III-C), the separation by condition number is less sharp. Nonetheless, a trend remains visible, confirming that κ provides a useful, though approximate, indicator of convergence behavior. VII. C OMPUTATIONAL PERFORMANCE Having identified α = 3 with the iterative-first ordering as the most effective configuration in Sec. VI, we now evaluate its computational performance. The algorithm is implemented in Julia [44] using the CUDA.jl package [45] and benchmarked on an NVIDIA H100 GPU. As CPU baselines we use

8

27%

1%

13% 17% 15% 9% 23% 13% 1000

2000 3000 Per-operation wall time (ms)

4000

Wall time (s)

2%

4% 11%

Fig. 4. Whole solver wall-time composition on case2869pegase (NVIDIA H100, (Nn , Nk ) = (5, 5)), for the Neumann and Bi-CGSTAB inner solvers. Costs are grouped and split at the divider into one-time work (left: data loading, allocation and the reference Jacobian inversion) and per-iteration work (right). Percentages give each group’s share of the total wall time, the bold number is the total wall time. Because the breakdown is measured with per-operation synchronization, the segments sum to the bar total, which exceeds the overlapped production wall time reported in Sec. VII-C.

LightSim2Grid on a single i7-13700H core [46], and 80-core implementations of LightSim2Grid and pandapower [42] on an AMD EPYC 9V84, with numbers taken from [47]. All results are reported for the case2869pegase test case [48], [49] with 2869 buses representing a realistic transmissionscale N−1 study. Consistent with using the reference Jacobian as preconditioner, we hot-start the NR iteration from the reference state xr , the N − 0 solution. The implementation used to produce the results in this paper is available at [50]. A. Per-iteration runtime breakdown To identify where the runtime is spent and to guide further optimization, Fig. 4 decomposes the total wall time of a complete N − 1 solve into grouped components, separated into one-time and per-iteration costs. Two points should be noted. First, the breakdown is a profiling run in which each operation is measured under synchronization. The segments therefore sum to a total, but that total exceeds the overlapped wall time of the production solve (Sec. VII-C), where the GPU scheduler overlaps independent operations. Second, the reported time is the total across all contingencies, not the pertopology time used in the CPU comparison. The breakdown is thus meant to show the relative distribution of work, not absolute performance. The one-time costs are dominated by setup and allocation. 3 The O(Nm ) inversion of the reference Jacobian accounts for only a small share, as J r (xr )−1 is computed once and reused across every contingency and NR iteration. We apply it through a Tensor-Core dense matrix–matrix product (GEMM) in the reduced-precision TF32 format (TensorFloat32), roughly doubling the throughput of the dominant dense product, while the outer Newton iteration and the low-rank correction remain in double precision (FP64). The preconditioner is stored once in single precision (FP32) and at each application its operands are cast down from FP64 and the result cast back up, a type conversion visible in the breakdown. The TF32 GEMM nonetheless outperforms the equivalent FP32 or FP64 product, so this overhead is recovered.

Bi-CGSTAB LS2G CPU 80C PP CPU 80C

Nk = 1 Nk = 2 Nk = 3 Nk = 4 Nk = 5 Nk = 6

1.5 1.0 0.5 0.0

5000

Residual

3% 1% 6% 3%

22%

5021 ms

16%

2.0

Converged (%)

0

14%

Neumann

SpMM GEMM Jacobian Other NR ops Vec ops Woodbury Type conversion

2434 ms

Neumann Bi-CGSTAB

Setup & allocation Reference inversion

min-max 1-99% 5-95% 11.5-88.5% 25-75% median

100 10−3 10−6 10−9 10−12 100 75 50 25 0

0

1 2 3 4 5 NR iterations Nn

6

0

1 2 3 4 5 NR iterations Nn

6

Fig. 5. Wall time (top row), residual distribution for Nk = 5 (middle row) and percentage of converged topologies (bottom row) for Neumann (left column) and Bi-CGSTAB (right column) as inner solvers on case2869pegase, for varying numbers of NR iterations Nn (horizontal axis). In the top and bottom rows, colors denote the inner-solver iteration count Nk and in the middle row, the shaded bands show percentile ranges of the final residual (see legend), at the fixed Nk = 5 with the dotted line marking the convergence threshold on the residual. In the top row, dashed (dotted) horizontal lines mark the LightSim2Grid (pandapower) 80-core CPU baselines.

Apart from the type conversions, the per-iteration cost is dominated by the linear algebra: the GEMM applying the reference inverse and the sparse matrix–matrix product (SpMM) applying the perturbation together form the largest share for both solvers. Bi-CGSTAB additionally spends a comparable amount in its vector recurrences (dot products and updates), which the Neumann series does not require. The Jacobian assembly is a small to moderate per-iteration contributor (4.34% and 9.37%, respectively), and the residual and edge-term evaluations are smaller still, summarized under the label Other NR steps, together with the calculation of (3) matrices K t (xt ) and ∆J (3) r (xt ), the state vector update itself, among other small contributors. In contrast to the whole solver workload, the rank-four Woodbury correction is negligible, contributing only 2.06% of the runtime for Neumann and 1.14% for Bi-CGSTAB. B. Iteration count trade-off Before comparing against CPU baselines, we determine the operating point (Nn , Nk ), i.e., the number of outer NR iterations Nn and inner solver (Neumann or Bi-CGSTAB) iterations Nk , that achieves acceptable convergence. Convergence for a contingency case is defined as a maximum residual

9

∥f t (xt )∥ < 10−9 . Figure 5 shows the wall time (top row), the distribution of the residual using median, quartiles and outer percentiles (middle row) and the percentage of converged cases (bottom row). Three observations follow. First, for a sufficiently accurate inner solve (large Nk , best seen for Bi-CGSTAB), four outer NR iterations suffice to converge the large majority of contingency cases. This confirms the N−1 states are close enough to the reference state for the outer loop to converge rapidly. To ensure convergence with some margin, we set (Nn , Nk ) = (5, 5) in the following. Second, at (Nn , Nk ) = (5, 5) the Neumann solver reaches wall times competitive with the LightSim2Grid 80-core CPU baseline, but does not converge for all cases: it reaches roughly 90 % coverage, as a fraction of contingencies have iterationmatrix spectral radii exceeding unity, for which the Neumann series diverges. This makes the Neumann solver suitable for fast screening, where incomplete coverage is acceptable, but not for exact N −1 assessment. Third, the Bi-CGSTAB solver reaches about 99 % convergence across the N − 1 cases, matching the coverage of the CPU baseline, but at wall times exceeding the 80core LightSim2Grid baseline. This is due to Bi-CGSTAB’s higher numerical cost per iteration, from the additional matrix products and vector operations within each step. Across all configurations, both solvers stay below the pandapower 80core benchmark. C. Comparison to CPU baselines Table I reports the per-topology time on case2869pegase, introduced above, together with two further test cases, case118 and case1354pegase [48], [49], at the fixed operating point (Nn , Nk ) = (5, 5). We report per-topology rather than total wall time so that grids of different size can be compared on a common basis. Additionally, the GPU solvers and the 80-core CPU baselines are evaluated on the exact same set of contingencies. Both measurements cover the same stages of the computation: the timed region is the solver itself, the repeated linear solves, and excludes the one-time setup of reading the grid data, allocation and, for the GPU, host–device transfer. This overhead differs between implementations and can be reduced substantially in a production setting [47]. The dense reference inverse’s memory footprint is modest: at single precision it is about 104 MB for case2869pegase, and even for the much larger case9241pegase it remains around 1 GB. The constraint on how many contingencies can be solved in parallel is therefore the batched working set, not the inverse. The GPU Neumann solver undercuts the LightSim2Grid 80-core per-topology time on case1354pegase (no 80-core LightSim2Grid value is available for case118), and stays close to it on case2869pegase. The Bi-CGSTAB solver remains above the LightSim2Grid 80-core baseline across all cases but consistently undercuts the pandapower 80-core benchmark. The two baselines differ in intent: LightSim2Grid is a performance-optimized C++ PF backend and thus the stronger reference, whereas pandapower is a broader, general-purpose Python tool not specialized for high-throughput contingency analysis.

TABLE I P ER - TOPOLOGY WALL TIME (ms) FOR N −1 SECURITY ASSESSMENT AT (Nn , Nk ) = (5, 5) ON THREE TEST GRIDS . T HE GPU N EUMANN AND B I -CGSTAB VARIANTS (NVIDIA H100) ARE COMPARED AGAINST L IGHT S IM 2G RID (LS2G) AND PANDAPOWER ON 80 CPU CORES AND LS2G ON A SINGLE CORE . C ORE COUNTS ARE INDICATED AS (80 C )/(1 C ). T HE FASTEST RESULT PER CASE IS HIGHLIGHTED IN BOLD . CPU BASELINES ARE TAKEN FROM [47] (80- CORE ) AND [46] ( SINGLE - CORE ).

# N −1 cases GPU Neumann GPU Bi-CGSTAB LS2G (80c) pandapower (80c) LS2G (1c)

case118 177

case1354pegase 1430

case2869pegase 3804

0.097 0.162 – 1.53 0.047

0.092 0.185 0.203 0.434 1.06

0.155 0.342 0.126 0.557 2.35

We also include a single-core baseline. It performs best only on the small case118 grid, where parallelization overhead does not yet pay off. On the two larger grids it is clearly outpaced by both the GPU and the 80-core CPU solvers. VIII. C ONCLUSION We have presented a GPU-parallel Newton–Raphson solver for N −1 contingency analysis that separates the power-flow state from the network topology, allowing a single reference factorization to be reused across all contingencies. The reference Jacobian serves as an effective preconditioner for the iterative inner solve, while the topology change is absorbed by a low-rank Woodbury correction, removing the need to refactorize the Jacobian for every contingency. A spectral analysis confirms the preconditioner’s effectiveness; both the Neumann series and Bi-CGSTAB converge in few iterations, with Neumann faster and Bi-CGSTAB handling the outlier cases where the series diverges. GPU parallelization of batched contingency analysis is effective but does not render optimized CPU methods obsolete: a single GPU proves comparable to 80 highly optimized CPU cores. The appropriate framing is thus one of complementary hardware rather than order-of-magnitude GPU superiority. Several directions follow. The runtime breakdown points to the precision conversion and the Bi-CGSTAB vector recurrences as the main optimization targets. Second, fast contingency screening, a cheaper, low-accuracy pass to pre-filter the a priori unknown set of critical contingencies, could further reduce the overall cost. Finally, the method extends naturally to further contingency types, such as generator and busbar outages. The latter remove an entire node and thus yield higher-rank corrections. R EFERENCES [1] Jülich Supercomputing Centre, “JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre,” Journal of large-scale research facilities, vol. 7, no. A182, 2021. [Online]. Available: http://dx.doi.org/ 10.17815/jlsrf-7-182 [2] J. J. Grainger and W. D. Stevenson Jr., Power System Analysis. New York: McGraw-Hill, 1994. [3] A. J. Wood, B. F. Wollenberg, and G. B. Sheblé, Power Generation, Operation and Control, 3rd ed. New York: John Wiley & Sons, 2014. [4] P. S. Kundur and O. P. Malik, Power system stability and control. McGraw-Hill Education, 2022.

10

[5] J. W. Simpson-Porco, F. Dörfler, and F. Bullo, “Voltage collapse in complex power grids,” Nat. Commun., vol. 7, no. 1, p. 10790, 2016. [6] D. Witthaut et al, “Collective nonlinear dynamics and self-organization in decentralized power grids,” Reviews of modern physics, vol. 94, no. 1, p. 015005, 2022. [7] B. Stott, J. Jardim, and O. Alsac, “DC power flow revisited,” IEEE Trans. Power Syst., vol. 24, no. 3, p. 1290, 2009. [8] K. Purchala et al, “Usefulness of DC power flow for active power flow analysis,” in IEEE Power Engineering Society General Meeting, 2005. IEEE, 2005, pp. 454–459. [9] C. Coffrin, P. Van Hentenryck, and R. Bent, “Approximating line losses and apparent power in ac power flow linearizations,” in Power and Energy Society General Meeting, 2012 IEEE. IEEE, 2012, pp. 1–8. [10] Z. Li, J. Yu, and Q. Wu, “Approximate linear power flow using logarithmic transform of voltage magnitudes with reactive power and transmission loss consideration,” IEEE Trans. Power Syst., vol. 33, no. 4, pp. 4593–4603, 2018. [11] Z. Yang et al, “A general formulation of linear power flow models: basic theory and error analysis,” IEEE Trans. Power Syst., vol. 34, no. 2, pp. 1315–1324, 2019. [12] T. Guler, G. Gross, and M. Liu, “Generalized line outage distribution factors,” IEEE Transactions on Power Systems, vol. 22, no. 2, pp. 879– 881, 2007. [13] H. Ronellenfitsch et al, “Dual theory of transmission line outages,” IEEE transactions on power systems, vol. 32, no. 5, pp. 4060–4068, 2017. [14] J. van Dijk et al, “Unified algebraic derivation of distribution factors in linear power flow,” arXiv preprint arXiv:2412.16164, 2024. [15] M. Titz, D. Witthaut, J. van Dijk, B. Petrick, and N. Westerbeck, “Voltage-sensitive distribution factors for contingency analysis and topology optimization,” 2025. [Online]. Available: https://arxiv.org/abs/ 2509.19976 [16] O. Alsac, B. Stott, and W. F. Tinney, “Sparsity-oriented compensation methods for modified network solutions,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-102, no. 5, pp. 1050–1060, 1983. [17] P. Stange, A. Griewank, and M. Bollhöfer, “On the efficient update of rectangular LU factorizations subject to low rank modifications,” https: //d-nb.info/1248217799/34, 2006. [18] M. A. Saunders et al, “LUSOL: Sparse LU for Ax= b,” https:// web.stanford.edu/group/SOL/software/lusol/, 2009. [19] W. W. Hager, “Updating the inverse of a matrix,” SIAM Review, vol. 31, no. 2, pp. 221–239, 1989. [Online]. Available: https: //doi.org/10.1137/1031049 [20] N. J. Higham, Accuracy and stability of numerical algorithms. SIAM, 2002. [21] N. Westerbeck, J. van Dijk, J. Viebahn, C. Merz, and D. Witthaut, “Accelerated dc loadflow solver for topology optimization,” in 2025 IEEE Kiel PowerTech, 2025, pp. 1–8. [22] G. Zhou, R. Bo, L. Chien, X. Zhang, S. Yang, and D. Su, “Gpuaccelerated algorithm for online probabilistic power flow,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 1132–1135, 2018. [23] Z. Wang, S. Wende-von Berg, and M. Braun, “Fast parallel newton–raphson power flow solver for large number of system calculations with cpu and gpu,” Sustainable Energy, Grids and Networks, vol. 27, p. 100483, 2021. [Online]. Available: https: //www.sciencedirect.com/science/article/pii/S2352467721000540 [24] X. Su, C. He, T. Liu, and L. Wu, “Full parallel power flow solution: A gpu-cpu-based vectorization parallelization and sparse techniques for newton–raphson implementation,” IEEE Transactions on Smart Grid, vol. 11, no. 3, pp. 1833–1844, 2020. [25] A. Ahmadi, M. C. Smith, E. R. Collins, V. Dargahi, and S. Jin, “Fast Newton-Raphson power flow analysis based on sparse techniques and parallel processing,” IEEE Transactions on Power Systems, vol. 37, no. 3, pp. 1695–1705, 2022. [26] G. Zhou, Y. Feng, R. Bo, L. Chien, X. Zhang, Y. Lang, Y. Jia, and Z. Chen, “GPU-accelerated batch-ACPF solution for N-1 static security analysis,” IEEE Transactions on Smart Grid, vol. 8, no. 3, pp. 1406– 1416, 2017. [27] S. Huang and V. Dinavahi, “Real-time contingency analysis on massively parallel architectures with compensation method,” IEEE Access, vol. 6, pp. 44 519–44 530, 2018. [28] A. Semlyen, “Fundamental concepts of a Krylov subspace power flow methodology,” IEEE Transactions on Power Systems, vol. 11, no. 3, pp. 1528–1537, 1996. [29] E. F. Yetkin and O. Ceylan, “Recycling Newton–Krylov algorithm for efficient solution of large scale power systems,” International Journal of Electrical Power & Energy Systems, vol. 144, p. 108559, 2023.

[Online]. Available: https://www.sciencedirect.com/science/article/pii/ S0142061522005555 [30] X. Meng and Z. Yun, “GPU-based N-1 AC power flow with an improved-convergence fast decoupled method considering memory limitations in large power grids,” in 2023 International Conference on Power System Technology (PowerCon), 2023, pp. 1–6. [31] M. Fu et al, “GPU-based N-1 static security analysis algorithm with preconditioned conjugate gradient method,” IEEE Access, vol. 8, pp. 124 066–124 075, 2020. [32] M. D’orto, S. Sjöblom, L. S. Chien, L. Axner, and J. Gong, “Comparing different approaches for solving large scale power-flow problems with the Newton-Raphson method,” IEEE Access, vol. 9, pp. 56 604–56 615, 2021. [33] V. W. Lee, C. Kim, J. Chhugani, M. Deisher, D. Kim, A. D. Nguyen, N. Satish, M. Smelyanskiy, S. Chennupaty, P. Hammarlund, R. Singhal, and P. Dubey, “Debunking the 100x GPU vs. CPU myth: an evaluation of throughput computing on CPU and GPU,” in Proceedings of the 37th Annual International Symposium on Computer Architecture, ser. ISCA ’10. New York, NY, USA: Association for Computing Machinery, 2010, p. 451–460. [Online]. Available: https://doi.org/10.1145/1815961.1816021 [34] F. Milano, “Complex frequency,” IEEE Transactions on Power Systems, vol. 37, no. 2, pp. 1230–1240, 2022. [35] M. Hunkemöller, D. Witthaut, L. Schewe, J. van Dijk, M. Farber, and N. Westerbeck, “A GPU-based Newton–Raphson algorithm for power flow under topology changes,” in 2026 Open Source Modelling and Simulation of Energy Systems (OSMSES), 2026, pp. 1–6. [36] R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. Cambridge University Press, 2012. [37] R. S. Varga, Matrix Iterative Analysis, 2nd ed., ser. Springer Series in Computational Mathematics. Berlin: Springer, 2000, vol. 27. [38] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. Society for Industrial and Applied Mathematics, 2003. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/1.9780898718003 [39] H. A. van der Vorst, “Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992. [Online]. Available: https://doi.org/10.1137/0913035 [40] A. Greenbaum, V. Pták, and Z. Strakoš, “Any nonincreasing convergence curve is possible for GMRES,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 3, pp. 465–469, 1996. [Online]. Available: https://doi.org/10.1137/S0895479894275030 [41] L. N. Trefethen and M. Embree, Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton, NJ: Princeton University Press, 2005. [42] L. Thurner, A. Scheidler, F. Schäfer, J.-H. Menke, J. Dollichon, F. Meier, S. Meinecke, and M. Braun, “Pandapower—an open-source python tool for convenient modeling, analysis, and optimization of electric power systems,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6510– 6521, 2018. [43] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power systems research and education,” IEEE Transactions on Power Systems, vol. 26, no. 1, pp. 12–19, 2011. [44] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017. [Online]. Available: https://doi.org/10.1137/141000671 [45] T. Besard, C. Foket, and B. De Sutter, “Effective extensible programming: Unleashing julia on gpus,” IEEE Transactions on Parallel and Distributed Systems, vol. 30, no. 4, pp. 827–841, 2019. [46] B. Donnot, “Lightsim2grid - A C++ backend targeting the Grid2Op platform. ,” https://GitHub.com/Grid2Op/lightsim2grid, 2020. [47] M. Farber, “Massively parallel GPU-based sparse linear solver for large scale power system problems,” Master’s thesis, Technische Universität Berlin, 2025. [48] C. Josz, S. Fliscounakis, J. Maeght, and P. Panciatici, “AC power flow data in matpower and qcqp format: itesla, rte snapshots, and pegase,” 2016. [Online]. Available: https://arxiv.org/abs/1603.01533 [49] S. Fliscounakis et al, “Contingency ranking with respect to overloads in very large power systems taking into account uncertainty, preventive, and corrective actions,” IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 4909–4917, 2013. [50] M. Hunkemöller, “Code: Low-rank corrections to a shared reference Jacobian for GPU-parallel AC N-1 contingency analysis,” Zenodo, Oct. 2026. [Online]. Available: https://doi.org/10.5281/zenodo.23662128

Related documents

• An empirical comparison of the performance of classical power indices#242057
• Suggested Modifications to IEEE 2800-2022 Requirements for Inverter-Based Resources Under Unbalanced Faults: Preprint#1481086
• Interactive Power Flow in the Browser#1157193
• Hybrid machine learning framework for electrical grid stability prediction using comparative model evaluation.#1251920
• Modular Model Order Reduction for Inverter-Based Power Systems with Stability and Accuracy Guarantees#1186791