add fixed point diagnostics

[luke-kiernan]

Fixed point jacobian diagnostics, in collaboration with @jmaack24 . Consider g(x) = x - J^{-1}(x)F(x): our solution of NR is a fixed point of that function. If the largest eigenvalue of the Jacobian of that function is bigger than 1, then that indicates likely trouble. Performance penalty of computing the largest eigenvalue is very minimal, only 70ms per iteration for a 8k bus system.

[luke-kiernan]

Creating a new factorization each time here is less-than-efficient. Might be not worth fixing, though--if I recall correctly, the KrylovKit eigensolve is the bottleneck here.

Also, I'm using closures and boxing here: F isn't an argument of matvec.

[Copilot]

Pull request overview

Adds a Newton fixed-point Jacobian spectral radius diagnostic to help detect likely Newton-Raphson convergence issues (based on the local fixed-point contraction property of (g(x)=x-J(x)^{-1}F(x))), along with tests and a solver option to enable per-iteration monitoring.

Changes:

• Introduces matrix-free fixed-point Jacobian spectral radius estimation using KrylovKit eigen-solves and KLU back-solves.
• Adds an ACPowerFlow option flag to enable per-iteration logging of ( \rho ), ( |F|_\infty ) location, and an estimated ( \kappa(J) ).
• Adds tests for the Hessian-vector bilinear form, the fixed-point Jacobian-vector product, and logging/flag behavior.

Reviewed changes

Copilot reviewed 6 out of 7 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
src/fixed_point_spectral_radius.jl	Implements acpf_hvvp, spectral radius estimation, and condition estimation.
src/power_flow_method.jl	Adds optional per-iteration diagnostic logging (NR/TR) and an x0 warning/log line.
src/power_flow_types.jl	Adds compute_fixed_point_spectral_radius::Bool to ACPowerFlow and getter.
src/PowerFlowData.jl	Delegating getter for the new option flag on PowerFlowData.
src/PowerFlows.jl	Imports KrylovKit and includes the new implementation file.
test/test_fixed_point_spectral_radius.jl	Adds convergence-order tests and smoke/logging tests for the new diagnostic.
Project.toml	Adds KrylovKit dependency and compat entry.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[github-actions]

Performance Results

Precompile Time

Main	This Branch	Delta
3.034 s	2.983 s	-1.7%

Solve Time

Polar AC

Test	Main	This Branch	Delta
NewtonRaphsonACPowerFlow First Solve	12.145 s	16.265 s	+33.9%
NewtonRaphsonACPowerFlow Second Solve	98.4 ms	92.2 ms	-6.3%
RobustHomotopyPowerFlow First Solve	18.227 s	18.612 s	+2.1%
RobustHomotopyPowerFlow Second Solve	10.216 s	10.5 s	+2.8%
NewtonRaphsonACPowerFlow(iwamoto) First Solve	964.3 ms	217.7 ms	-77.4%
NewtonRaphsonACPowerFlow(iwamoto) Second Solve	87.9 ms	93.0 ms	+5.8%
TrustRegionACPowerFlow(iwamoto) First Solve	1.914 s	2.114 s	+10.5%
TrustRegionACPowerFlow(iwamoto) Second Solve	97.8 ms	92.4 ms	-5.4%

Rectangular CI

Test	Main	This Branch	Delta
ACRectangularPowerFlow{NR} First Solve	5.263 s	6.398 s	+21.6%
ACRectangularPowerFlow{NR} Second Solve	43.1 ms	48.1 ms	+11.7%
ACRectangularPowerFlow{NR}(iwamoto) First Solve	164.0 ms	170.7 ms	+4.1%
ACRectangularPowerFlow{NR}(iwamoto) Second Solve	45.0 ms	45.3 ms	+0.7%
ACRectangularPowerFlow{TR} First Solve	5.71 s	7.173 s	+25.6%
ACRectangularPowerFlow{TR} Second Solve	44.8 ms	46.9 ms	+4.7%
ACRectangularPowerFlow{TR}(iwamoto_fallback) First Solve	170.4 ms	168.6 ms	-1.1%
ACRectangularPowerFlow{TR}(iwamoto_fallback) Second Solve	44.1 ms	45.6 ms	+3.5%

Mixed CPB

Test	Main	This Branch	Delta
ACMixedPowerFlow{NR} First Solve	5.072 s	6.036 s	+19.0%
ACMixedPowerFlow{NR} Second Solve	47.7 ms	46.5 ms	-2.6%
ACMixedPowerFlow{NR}(iwamoto) First Solve	255.5 ms	172.6 ms	-32.4%
ACMixedPowerFlow{NR}(iwamoto) Second Solve	46.8 ms	47.0 ms	+0.4%
ACMixedPowerFlow{TR} First Solve	5.159 s	6.473 s	+25.5%
ACMixedPowerFlow{TR} Second Solve	102.6 ms	49.4 ms	-51.8%
ACMixedPowerFlow{TR}(iwamoto_fallback) First Solve	173.2 ms	212.4 ms	+22.6%
ACMixedPowerFlow{TR}(iwamoto_fallback) Second Solve	47.5 ms	75.2 ms	+58.1%
ACMixedPowerFlow{LM} First Solve	7.583 s	8.036 s	+6.0%
ACMixedPowerFlow{LM} Second Solve	1.433 s	1.413 s	-1.4%

DC

Test	Main	This Branch	Delta
DCPowerFlow First Solve	4.292 s	4.169 s	-2.9%
DCPowerFlow Second Solve	18.7 ms	18.3 ms	-2.3%
PTDFDCPowerFlow First Solve	1.41 s	1.431 s	+1.5%
PTDFDCPowerFlow Second Solve	85.4 ms	85.2 ms	-0.3%
vPTDFDCPowerFlow First Solve	5.88 s	5.792 s	-1.5%
vPTDFDCPowerFlow Second Solve	5.377 s	5.415 s	+0.7%

[luke-kiernan]

Copilot review responses:

• Return type doc (line 150): fixed — now documents (ρ, info, condest).
• Math identity (line 169): fixed — rewritten as (Gv)_k = (J⁻¹ w)_k with w_k = vᵀ H_k u, matching what acpf_hvvp actually computes.
• Docstring indexing (line 20): trimmed the PQ-only sentence and added a pointer to _create_jacobian_matrix_structure, which already documents the bus-type-dependent state vector convention. Better to have one source of truth than restate it here.
• Deterministic v_init (line 239): fixed — now ones(n)/sqrt(n).
• Re-factoring KLU each call (lines 234, 237, power_flow_method.jl:613): known and noted in my own self-review. The KrylovKit eigsolve dominates wall time (verified on CATS); the per-iteration extra factorization is in the noise compared to that, and reusing KLULinSolveCache here would require plumbing the cache through the diagnostic API. Deferring to a follow-up if profiling shows it matters.

[luke-kiernan]

Reverting to draft for now--maintence state of the KrylovKit dependence is mildly iffy

[luke-kiernan]

Actually the maintenance status doesn't look too bad: KrylovKit.jl is maintained by one person, but there's commits from 3 months ago and last release was October of 2025.

ArnoldiMethod.jl is an alternative: it's used elsewhere in the Julia ecosystem and has commits from a dozen maintainers. However, the most recent commit is 2 years ago.

[luke-kiernan]

Decision: KrylovKit.jl is fine for now. We'll use this branch to prototype, the meta-algorithm of fallbacks, then open an new PR so as to avoid repeated rebases/merges. At that point we may implement our own.

[luke-kiernan]

Okay I've pushed a first draft of an adaptive method.

Strategy: start a default lightweight-ish Newton-type strategy (polar TR at the moment). Compute 3 metrics each iteration: spectral radius ρ, Jacobian condition estimate κ̂, and residual size ‖F‖_∞. If ρ > 1 for 2 consequective iterations, then switch to something more robust (polar LM at the moment). If we see less than 10% improvement for 3 iterations and the current residual is less than 10x the theoretical noise floor (κ̂·ε·‖F_0‖_∞, ε datatype precision), then we give up.

Here's what that looks like for one of the tougher time steps for CATS:

[ Info: Finding subnetworks via iterative union find
┌ Warning: Using caller-provided x0; skipping improve_x0.
└ @ PowerFlows ~/Documents/julia/Sienna/psy5/PowerFlows.jl/src/power_flow_setup.jl:295
[ Info: Initial residual size: 78.18302991625828 L2, 54.67301427002515 L∞
[ Info: x0 (time_step 22): ρ = 3.874, ‖F‖_∞ = 54.67 at bus 1951 (P), κ̂(J) = 9.271e8
[ Info: Adaptive TR iter 0: ρ = 0.6472, ‖F‖_∞ = 33.76 at bus 4209 (Q), κ̂(J) = 7.897e8
[ Info: Adaptive TR iter 1: ρ = 0.9246, ‖F‖_∞ = 8.638 at bus 8268 (Q), κ̂(J) = 1.605e9
[ Info: Adaptive TR iter 2: ρ = 17.29, ‖F‖_∞ = 0.6821 at bus 2806 (Q), κ̂(J) = 1.09e10
[ Info: Adaptive TR iter 3: ρ = 0.664, ‖F‖_∞ = 0.015 at bus 595 (P), κ̂(J) = 9.401e8
[ Info: Adaptive TR iter 4: ρ = 0.9293, ‖F‖_∞ = 0.004697 at bus 3664 (P), κ̂(J) = 1.988e9
[ Info: Adaptive TR iter 5: ρ = 808.4, ‖F‖_∞ = 0.002076 at bus 595 (P), κ̂(J) = 7.621e10
[ Info: Adaptive TR iter 6: ρ = 786.1, ‖F‖_∞ = 0.002076 at bus 595 (P), κ̂(J) = 7.621e10
[ Info: Adaptive solver: ρ ≥ 1.0 for 2 consecutive iters; escalating to LevenbergMarquardt at iter 6
[ Info: LM iter 0: ρ = 1274.0, ‖F‖_∞ = 0.001881 at bus 595 (P), κ̂(J) = 7.621e10
[ Info: LM iter 1: ρ = 0.6429, ‖F‖_∞ = 0.001881 at bus 595 (P), κ̂(J) = 2.119e9
[ Info: LM iter 2: ρ = 282.2, ‖F‖_∞ = 0.0001008 at bus 3012 (P), κ̂(J) = 4.533e10
[ Info: LM iter 3: ρ = 53.01, ‖F‖_∞ = 0.0001008 at bus 3012 (P), κ̂(J) = 1.959e10
[ Info: LM iter 4: ρ = 53.01, ‖F‖_∞ = 0.0001008 at bus 3012 (P), κ̂(J) = 1.959e10
[ Info: LM hit numerical floor: ‖F‖_∞ = 0.000101, est. floor κ̂·ε·max(1,‖F‖_init) ≈ 0.000238
[ Info: Final residual size: 0.0005938791533058349 L2, 0.00010075019468018531 L∞.
┌ Warning: The AdaptiveACPowerFlow solver stopped at the numerical floor after 11 iterations: the residual is limited by κ̂·ε·max(1,‖F‖_init) backward stability, not the requested tolerance. Reporting non-convergence.
└ @ PowerFlows ~/Documents/julia/Sienna/psy5/PowerFlows.jl/src/power_flow_method.jl:964

I'm not set on the particular methods here--seems like something we could experiment with--but I do want feedback on my approach. Do my "switching strategies" criteria make sense? Are there other cases I need to address, besides hitting the noise floor and converging?

[luke-kiernan]

Okay the noise floor logic here needs work:

1. κ̂·ε·‖F_0‖_∞ is a worst-case, where our residual is aligned with the smallest singular vector of J. In practice, that rarely happens.
2. I was using ε = eps(Float64), because that fit what I was seeing empirically. But our y bus uses 32-bit data types! So we should be using ε = eps(Float32), in which case the formula says we should loose precision much much sooner than what I'm seeing.

I'll replace it with "bail out when we stop making progress," and change the language so we're not attributing it to loss of precision.

[luke-kiernan]

I've introduced a full-fledged convergence failure/search stagnation analysis. There's now scoped enum called ACPowerFlowSolveStatus with the following values:

• CONVERGED: ‖F‖_∞ < tol.
• NON_ROOT_LOCAL_MIN: gradient ≈ 0, positive directional curvature
• NON_ROOT_SADDLE: gradient ≈ 0, negative directional curvature
• NON_ROOT_STATIONARY: gradient ≈ 0, curvature inconclusive
• SINGULAR_SUBSPACE: fixed point trapped on a near-singular subspace of J
• FOLD: ‖½H[Δx,Δx]‖/‖F‖ >> 1, quadratic curvature dominates the Newton step
• LIMIT_CYCLE: settled into a cycle of period ≥ 2.
• STAGNATED_OTHER: ‖F‖_∞ flat but no specific diagnostic threshold
• DIVERGED: LM-specific: μ hit cap or λ went non-finite.
• FAILED: hit maxIterations with no stagnation or divergence.
• RUNNING: internal sentinel; never returned to callers.

We may not want all of these long-term, but they're certainly useful for troubleshooting convergence failure and prototyping hybrid strategies.

New solver_settings:

• detect_stagnation (default true): do the convergence failure/search stagnation analysis
• bail_on_stagnation (default false, true for adaptive): bail when the search stagnates.
• monitor_ρ: per-iteration logging. Defaults to get_compute_fixed_point_spectral_radius(data)