From 248d7649db474ecbed046a6f7bb0dd015d3f5894 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 7 Apr 2026 21:08:05 -0600 Subject: [PATCH 01/14] add fixed point diagnostics --- Project.toml | 2 + src/PowerFlowData.jl | 2 + src/PowerFlows.jl | 2 + src/fixed_point_spectral_radius.jl | 245 +++++++++++++++++++++++ src/power_flow_method.jl | 38 ++++ src/power_flow_types.jl | 6 + test/test_fixed_point_spectral_radius.jl | 175 ++++++++++++++++ 7 files changed, 470 insertions(+) create mode 100644 src/fixed_point_spectral_radius.jl create mode 100644 test/test_fixed_point_spectral_radius.jl diff --git a/Project.toml b/Project.toml index c091d51fb..c40724d7a 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1" JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" KLU = "ef3ab10e-7fda-4108-b977-705223b18434" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -25,6 +26,7 @@ Dates = "1" InfrastructureSystems = "3" JSON3 = "1" KLU = "^0.6" +KrylovKit = "0.10.2" LineSearches = "7.4.0" LinearAlgebra = "1" Logging = "1" diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 8989b5da0..2100af9a0 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -224,6 +224,8 @@ get_computed_gspf(pfd::PowerFlowData) = pfd.computed_generator_slack_participati get_calculate_loss_factors(pfd::PowerFlowData) = get_calculate_loss_factors(pfd.pf) get_calculate_voltage_stability_factors(pfd::PowerFlowData) = get_calculate_voltage_stability_factors(pfd.pf) +get_compute_fixed_point_spectral_radius(pfd::PowerFlowData) = + get_compute_fixed_point_spectral_radius(pfd.pf) get_network_reductions(pfd::PowerFlowData) = get_network_reductions(pfd.pf) get_time_steps(pfd::PowerFlowData) = get_time_steps(pfd.pf) get_time_step_names(pfd::PowerFlowData) = get_time_step_names(pfd.pf) diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index d7343b77f..b4885ee4c 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -37,6 +37,7 @@ import LinearAlgebra: norm, dot, ldiv!, mul! import LinearAlgebra: norm, dot import JSON3 import KLU +import KrylovKit import SparseArrays import InfrastructureSystems as IS import PowerNetworkMatrices as PNM @@ -72,6 +73,7 @@ include("rectangular_ci_power_flow_jacobian.jl") include("mixed_cpb_setup.jl") include("mixed_cpb_power_flow_residual.jl") include("mixed_cpb_power_flow_jacobian.jl") +include("fixed_point_spectral_radius.jl") include("solve_ac_power_flow.jl") include("power_flow_setup.jl") include("power_flow_method.jl") diff --git a/src/fixed_point_spectral_radius.jl b/src/fixed_point_spectral_radius.jl new file mode 100644 index 000000000..f2cd5c4eb --- /dev/null +++ b/src/fixed_point_spectral_radius.jl @@ -0,0 +1,245 @@ +""" + acpf_hvvp(data, time_step, v, u) -> Vector{Float64} + +Compute the bilinear form of every component Hessian of the AC power flow residual +with two vectors `v` and `u`: + +```math +w[k] = v^\\top \\left(\\frac{\\partial^2 F_k}{\\partial x \\partial x^\\top}\\right) u +``` + +where ``F_k`` is the ``k``-th component of the AC power flow residual at the current +state stored in `data` (i.e., `data.bus_magnitude` and `data.bus_angles`). + +This is matrix-free: the residual's 3-tensor of second derivatives is never +materialized. State-vector indexing follows the PowerFlows.jl convention: +`x[2i-1] = Vm` and `x[2i] = Va` for PQ buses, with the same column positions used +for variable lookups via `is_PQ` / `has_θ` gates. + +LCC state variables are not yet supported; the returned vector covers only the +`2 * num_buses` bus state entries. +""" +function acpf_hvvp( + data::ACPowerFlowData, + time_step::Int, + v::AbstractVector{Float64}, + u::AbstractVector{Float64}, +) + Yb = data.power_network_matrix.data + Vm = view(data.bus_magnitude, :, time_step) + θ = view(data.bus_angles, :, time_step) + num_buses = first(size(data.bus_type)) + n = 2 * num_buses + @assert length(v) >= n && length(u) >= n + w = zeros(Float64, n) + + for i in 1:num_buses + bt_i = data.bus_type[i, time_step] + has_θi = (bt_i == PSY.ACBusTypes.PQ) || (bt_i == PSY.ACBusTypes.PV) + Pi_θiθi, Qi_θiθi = 0.0, 0.0 + Pi_Viθi, Qi_Viθi = 0.0, 0.0 + + for k in data.neighbors[i] + k == i && continue + bt_k = data.bus_type[k, time_step] + has_θk = (bt_k == PSY.ACBusTypes.PQ) || (bt_k == PSY.ACBusTypes.PV) + Gik, Bik = real(Yb[i, k]), imag(Yb[i, k]) + θik = θ[i] - θ[k] + s, c = sin(θik), cos(θik) + + if has_θi + d2P_θkθi = Vm[i] * Vm[k] * (Gik * c + Bik * s) + d2Q_θkθi = Vm[i] * Vm[k] * (Gik * s - Bik * c) + # contribution towards sum in ∂²Δ{Pᵢ,Qᵢ}/∂θᵢ∂θᵢ + Pi_θiθi -= d2P_θkθi + Qi_θiθi -= d2Q_θkθi + + if has_θk + # ∂²Δ{Pᵢ,Qᵢ}/∂θₖ∂θᵢ + # really ∂²ΔPᵢ/∂θₖ∂θᵢ * v[θₖ] * u[θᵢ] + ∂²ΔQᵢ/∂θᵢ∂θₖ * v[θᵢ] * u[θₖ] + # but ∂²ΔPᵢ/∂θₖ∂θᵢ equals ∂²ΔPᵢ/∂θᵢ∂θₖ. + w[2 * i - 1] += d2P_θkθi * (v[2 * k] * u[2 * i] + v[2 * i] * u[2 * k]) + w[2 * i] += d2Q_θkθi * (v[2 * k] * u[2 * i] + v[2 * i] * u[2 * k]) + end + end + + if bt_i == PSY.ACBusTypes.PQ + d2P_θkVi = Vm[k] * (Gik * s - Bik * c) + d2Q_θkVi = Vm[k] * (-Gik * c - Bik * s) + # contribution towards sum in ∂²Δ{Pᵢ,Qᵢ}/∂θᵢ∂Vᵢ + Pi_Viθi -= d2P_θkVi + Qi_Viθi -= d2Q_θkVi + + if has_θk + # ∂²Δ{Pᵢ,Qᵢ}/∂θₖ∂Vᵢ + w[2 * i - 1] += + d2P_θkVi * (v[2 * k] * u[2 * i - 1] + v[2 * i - 1] * u[2 * k]) + w[2 * i] += + d2Q_θkVi * (v[2 * k] * u[2 * i - 1] + v[2 * i - 1] * u[2 * k]) + end + end + + if has_θk + # ∂²Δ{Pᵢ,Qᵢ}/∂²θₖ + d2P_θkθk = Vm[i] * Vm[k] * (-Gik * c - Bik * s) + d2Q_θkθk = Vm[i] * Vm[k] * (-Gik * s + Bik * c) + w[2 * i - 1] += d2P_θkθk * v[2 * k] * u[2 * k] + w[2 * i] += d2Q_θkθk * v[2 * k] * u[2 * k] + end + + if bt_k == PSY.ACBusTypes.PQ + # ∂²Δ{Pᵢ,Qᵢ}/∂θₖ∂Vₖ + d2P_θkVk = Vm[i] * (Gik * s - Bik * c) + d2Q_θkVk = Vm[i] * (-Gik * c - Bik * s) + w[2 * i - 1] += + d2P_θkVk * (v[2 * k] * u[2 * k - 1] + v[2 * k - 1] * u[2 * k]) + w[2 * i] += + d2Q_θkVk * (v[2 * k] * u[2 * k - 1] + v[2 * k - 1] * u[2 * k]) + + if bt_i == PSY.ACBusTypes.PQ + # ∂²Δ{Pᵢ,Qᵢ}/∂Vₖ∂Vᵢ + d2P_VkVi = Gik * c + Bik * s + d2Q_VkVi = Gik * s - Bik * c + w[2 * i - 1] += + d2P_VkVi * + (v[2 * k - 1] * u[2 * i - 1] + v[2 * i - 1] * u[2 * k - 1]) + w[2 * i] += + d2Q_VkVi * + (v[2 * k - 1] * u[2 * i - 1] + v[2 * i - 1] * u[2 * k - 1]) + end + + if has_θi + # ∂²Δ{Pᵢ,Qᵢ}/∂Vₖ∂θᵢ + d2P_Vkθi = -Vm[i] * (Gik * s - Bik * c) + d2Q_Vkθi = -Vm[i] * (-Gik * c - Bik * s) + w[2 * i - 1] += + d2P_Vkθi * (v[2 * k - 1] * u[2 * i] + v[2 * i] * u[2 * k - 1]) + w[2 * i] += + d2Q_Vkθi * (v[2 * k - 1] * u[2 * i] + v[2 * i] * u[2 * k - 1]) + end + end + end + + # diagonal terms in i: accumulated sums [except ∂²Vᵢ] + if has_θi + # ∂²Δ{Pᵢ,Qᵢ}/∂²θᵢ + w[2 * i - 1] += Pi_θiθi * v[2 * i] * u[2 * i] + w[2 * i] += Qi_θiθi * v[2 * i] * u[2 * i] + end + + if bt_i == PSY.ACBusTypes.PQ + # ∂²Δ{Pᵢ,Qᵢ}/∂Vᵢ∂θᵢ + w[2 * i - 1] += + Pi_Viθi * (v[2 * i - 1] * u[2 * i] + v[2 * i] * u[2 * i - 1]) + w[2 * i] += + Qi_Viθi * (v[2 * i - 1] * u[2 * i] + v[2 * i] * u[2 * i - 1]) + + # ∂²Δ{Pᵢ,Qᵢ}/∂²Vᵢ + d2P_ViVi = 2.0 * real(Yb[i, i]) + d2Q_ViVi = -2.0 * imag(Yb[i, i]) + w[2 * i - 1] += d2P_ViVi * v[2 * i - 1] * u[2 * i - 1] + w[2 * i] += d2Q_ViVi * v[2 * i - 1] * u[2 * i - 1] + end + end + + return w +end + +""" + compute_fixed_point_spectral_radius(data, time_step; x0, tol, maxiter, krylovdim) + -> (ρ::Float64, info) + +Estimate the spectral radius of the Jacobian of the Newton fixed-point map + +```math +g(x) = x - J(x)^{-1} F(x) +``` + +at the state `x0` (default: the result of `calculate_x0(data, time_step)`). The +spectral radius ``\\rho(\\partial g / \\partial x)`` is a local convergence +diagnostic for Newton-Raphson: ``\\rho < 1`` is necessary (and locally sufficient) +for the fixed-point iteration to converge from a small enough neighborhood of +``x``. Values close to or above 1 indicate that NR is at risk of stalling or +diverging from this starting point. + +Differentiating ``g`` and using ``\\partial F / \\partial x = J`` gives the +identity + +```math +G := \\frac{\\partial g}{\\partial x} = J^{-1} \\sum_k H_k u_k, +\\qquad u = J^{-1} F(x) +``` + +where ``H_k`` is the Hessian of the ``k``-th residual component. ``G`` is never +materialized: the matvec ``G v`` is computed as ``J^{-1} \\cdot (\\sum_k H_k u_k v)`` +via two `KLU` back-solves and one [`acpf_hvvp`](@ref) call. KrylovKit's +matrix-free Lanczos / Arnoldi solver then extracts the eigenvalue of largest +magnitude. + +# Notes + +- At a true power-flow solution, ``F(x) = 0`` so ``u = 0`` and ``\\rho(G) = 0`` + trivially. The diagnostic is meaningful at *non-solution* iterates such as the + flat-start `x0`. +- LCC state variables are not yet supported. +""" +function compute_fixed_point_spectral_radius( + data::ACPowerFlowData, + time_step::Int; + x0::Union{Vector{Float64}, Nothing} = nothing, + tol::Float64 = 1e-6, + maxiter::Int = 200, + krylovdim::Int = 30, +) + if size(data.lcc.p_set, 1) > 0 + throw( + ArgumentError( + "compute_fixed_point_spectral_radius does not yet support LCC HVDC systems", + ), + ) + end + + residual = ACPowerFlowResidual(data, time_step) + jac = ACPowerFlowJacobian( + data, + residual.bus_slack_participation_factors, + residual.subnetworks, + time_step, + ) + x = isnothing(x0) ? calculate_x0(data, time_step) : copy(x0) + residual(x, time_step) + jac(time_step) + return _fixed_point_spectral_radius!( + data, residual, jac, time_step; + tol = tol, maxiter = maxiter, krylovdim = krylovdim, + ) +end + +"""In-place spectral radius computation that reuses an already-evaluated +`residual` and `jac` (both must have been called at the current state). Returns +`(ρ, info, condest)`, where `condest` is a Hager 1-norm estimate of the +condition number of `jac.Jv` computed from the same KLU factor used for the +spectral radius matvecs. Used by the per-iteration monitor inside +`_run_power_flow_method`.""" +function _fixed_point_spectral_radius!( + data::ACPowerFlowData, + residual::ACPowerFlowResidual, + jac::ACPowerFlowJacobian, + time_step::Int; + tol::Float64 = 1e-6, + maxiter::Int = 200, + krylovdim::Int = 30, +) + n = 2 * size(data.bus_type, 1) + F = KLU.klu(jac.Jv) + u = F \ copy(residual.Rv) + matvec(v::AbstractVector) = F \ acpf_hvvp(data, time_step, v, u) + v_init = randn(n) + vals, _, info = KrylovKit.eigsolve( + matvec, v_init, 1, :LM; + tol = tol, maxiter = maxiter, krylovdim = krylovdim, + ) + ρ = isempty(vals) ? NaN : abs(vals[1]) + condest = KLU.condest(F) + return ρ, info, condest +end diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index c9dbca459..7adfd1b57 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1,3 +1,21 @@ +"""Format a per-iteration diagnostic line: ρ, ‖F‖_∞ and the bus/equation where +that infinity-norm is attained, and a κ̂(J) condition estimate. All numerics +rounded to 4 significant figures.""" +function _log_diagnostics( + label::String, + ρ::Float64, + Rv::AbstractVector{Float64}, + condest::Float64, +) + abs_max, ix = findmax(abs, Rv) + pow_type = ix % 2 == 1 ? "P" : "Q" + bus_ix = div(ix + 1, 2) + sf(x) = round(x; sigdigits = 4) + @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at bus $bus_ix ($pow_type), " * + "κ̂(J) = $(sf(condest))" + return +end + """Cache for non-linear methods. # Fields @@ -585,6 +603,7 @@ function _run_power_flow_method(time_step::Int, validate_vms = validate_voltage_magnitudes i, converged = 1, false consecutive_reverts = 0 + monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) while i < maxIterations && !converged if iwamoto made_progress = _iwamoto_step( @@ -622,6 +641,10 @@ function _run_power_flow_method(time_step::Int, vm_validation_range, i, ) + if monitor_ρ + ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("NR iter $i", ρ, residual.Rv, condest) + end converged = norm(residual.Rv, Inf) < tol if !converged i += 1 @@ -680,6 +703,7 @@ function _run_power_flow_method(time_step::Int, linf = norm(residual.Rv, Inf) @debug "initially: sum of squares $(siground(residualSize)), L ∞ norm $(siground(linf)), Δ $(siground(delta))" + monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) while i < maxIterations && !converged delta = _trust_region_step( time_step, @@ -699,6 +723,10 @@ function _run_power_flow_method(time_step::Int, vm_validation_range, i, ) + if monitor_ρ + ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("TR iter $i", ρ, residual.Rv, condest) + end converged = norm(residual.Rv, Inf) < tol if !converged i += 1 @@ -802,6 +830,16 @@ function _newton_power_flow( pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol + if get_compute_fixed_point_spectral_radius(data) + ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) + _log_diagnostics("x0 (time_step $time_step)", ρ, residual.Rv, condest) + if ρ >= 1.0 + @warn "fixed-point spectral radius ρ = $(round(ρ; sigdigits = 4)) ≥ 1 " * + "at x0 (time_step $time_step); Newton-Raphson may not converge " * + "from this starting point" + end + end + i = 0 x_final = x0_init if !converged diff --git a/src/power_flow_types.jl b/src/power_flow_types.jl index a120164d9..2345bb878 100644 --- a/src/power_flow_types.jl +++ b/src/power_flow_types.jl @@ -152,6 +152,7 @@ struct ACPolarPowerFlow{ACSolver <: ACPowerFlowSolverType} <: AbstractACPowerFlo exporter::Union{Nothing, PowerFlowEvaluationModel} calculate_loss_factors::Bool calculate_voltage_stability_factors::Bool + compute_fixed_point_spectral_radius::Bool generator_slack_participation_factors::Union{ Nothing, Dict{Tuple{DataType, String}, Float64}, @@ -205,6 +206,7 @@ function ACPolarPowerFlow{ACSolver}(; exporter::Union{Nothing, PowerFlowEvaluationModel} = nothing, calculate_loss_factors::Bool = false, calculate_voltage_stability_factors::Bool = false, + compute_fixed_point_spectral_radius::Bool = false, generator_slack_participation_factors::Union{ Nothing, Dict{Tuple{DataType, String}, Float64}, @@ -234,6 +236,7 @@ function ACPolarPowerFlow{ACSolver}(; exporter, calculate_loss_factors, calculate_voltage_stability_factors, + compute_fixed_point_spectral_radius, generator_slack_participation_factors, enhanced_flat_start, robust_power_flow, @@ -275,6 +278,9 @@ get_calculate_loss_factors(pf::ACPolarPowerFlow) = pf.calculate_loss_factors get_calculate_voltage_stability_factors(::AbstractACPowerFlow) = false get_calculate_voltage_stability_factors(pf::ACPolarPowerFlow) = pf.calculate_voltage_stability_factors +get_compute_fixed_point_spectral_radius(::PowerFlowEvaluationModel) = false +get_compute_fixed_point_spectral_radius(pf::ACPolarPowerFlow) = + pf.compute_fixed_point_spectral_radius """ ACRectangularPowerFlow{ACSolver}(; kwargs...) where {ACSolver <: ACPowerFlowSolverType} diff --git a/test/test_fixed_point_spectral_radius.jl b/test/test_fixed_point_spectral_radius.jl new file mode 100644 index 000000000..8f2cc1820 --- /dev/null +++ b/test/test_fixed_point_spectral_radius.jl @@ -0,0 +1,175 @@ +@testset "test acpf_hvvp via O(ε²) finite-difference convergence" begin + # Verify acpf_hvvp via the convergence-order test: + # w_fd(ε) := [F(x+εv+εu) - F(x+εv-εu) - F(x-εv+εu) + F(x-εv-εu)] / (4ε²) + # = vᵀ H u + O(ε²) + # Halving ε should shrink the truncation error by ~4×. This is a principled + # correctness check that doesn't depend on a hand-picked atol. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PowerFlowData(pf, sys) + time_step = 1 + + residual = PF.ACPowerFlowResidual(data, time_step) + x0 = PF.calculate_x0(data, time_step) + n = length(x0) + + function residual_at(x) + residual(x, time_step) + return copy(residual.Rv) + end + + Random.seed!(42) + # normalize directions so ε is comparable across n + v = randn(n) ./ sqrt(n) + u = randn(n) ./ sqrt(n) + + function w_fd(ε) + return ( + residual_at(x0 .+ ε .* v .+ ε .* u) .- + residual_at(x0 .+ ε .* v .- ε .* u) .- + residual_at(x0 .- ε .* v .+ ε .* u) .+ + residual_at(x0 .- ε .* v .- ε .* u) + ) ./ (4 * ε^2) + end + + w_analytic = PF.acpf_hvvp(data, time_step, v, u) + residual_at(x0) # restore data state + + # Sweep ε in the regime where O(ε²) truncation dominates, well above the + # roundoff floor (~eps^(1/4) ≈ 1.2e-4 for a 4-point central scheme). + εs = [1e-2, 5e-3, 2.5e-3, 1.25e-3] + errs = [norm(w_analytic .- w_fd(ε), Inf) for ε in εs] + residual_at(x0) + + # Each halving of ε should shrink the error by ~4. Allow [3, 5] to absorb + # finite-precision noise. + ratios = errs[1:(end - 1)] ./ errs[2:end] + for r in ratios + @test 3.0 <= r <= 5.0 + end + # Sanity: errors are monotonically shrinking. + @test issorted(errs; rev = true) +end + +@testset "test fixed-point Jacobian-vector product via O(ε²) FD convergence" begin + # G(x) = ∂g/∂x where g(x) = x - J(x)⁻¹ F(x). The matrix-free product is + # G v = J⁻¹ · acpf_hvvp(v, u), u = J⁻¹ F(x) + # Verify it via central-difference convergence: + # (g(x+εv) - g(x-εv)) / (2ε) = G v + O(ε²) + # so the error halves by 4× when ε halves. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PowerFlowData(pf, sys) + time_step = 1 + + residual = PF.ACPowerFlowResidual(data, time_step) + jac = PF.ACPowerFlowJacobian( + data, + residual.bus_slack_participation_factors, + residual.subnetworks, + time_step, + ) + x0 = PF.calculate_x0(data, time_step) + n = length(x0) + + function g(x) + residual(x, time_step) + jac(time_step) + return x .- (jac.Jv \ copy(residual.Rv)) + end + + Random.seed!(7) + v = randn(n) ./ sqrt(n) + + Gv_fd(ε) = (g(x0 .+ ε .* v) .- g(x0 .- ε .* v)) ./ (2 * ε) + + # Matrix-free G·v via the same machinery the diagnostic uses. + residual(x0, time_step) + jac(time_step) + u = jac.Jv \ copy(residual.Rv) + Gv_analytic = jac.Jv \ PF.acpf_hvvp(data, time_step, v, u) + + εs = [1e-3, 5e-4, 2.5e-4, 1.25e-4] + errs = [norm(Gv_analytic .- Gv_fd(ε), Inf) for ε in εs] + g(x0) # restore + + ratios = errs[1:(end - 1)] ./ errs[2:end] + for r in ratios + @test 3.0 <= r <= 5.0 + end + @test issorted(errs; rev = true) +end + +@testset "test compute_fixed_point_spectral_radius smoke test" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PowerFlowData(pf, sys) + ρ, info, condest = PF.compute_fixed_point_spectral_radius(data, 1) + # c_sys14 is a well-conditioned 14-bus system; NR converges fast from flat + # start, so ρ should be safely below 1. + @test isfinite(ρ) + @test 0 < ρ < 1 + # condest is a 1-norm condition estimate; for a healthy 14-bus system it + # should be finite and well below the borderline-singular regime we saw + # on CATS (~1e9). + @test isfinite(condest) + @test condest > 0 + @test condest < 1e6 +end + +@testset "test compute_fixed_point_spectral_radius option flag" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; + correct_bustypes = true, + compute_fixed_point_spectral_radius = true, + ) + data = PowerFlowData(pf, sys) + # smoke test: solve runs end-to-end with the flag enabled + @test solve_power_flow!(data) +end + +@testset "test per-iteration spectral radius monitor emits log lines" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; + correct_bustypes = true, + compute_fixed_point_spectral_radius = true, + ) + data = PowerFlowData(pf, sys) + + # Capture log output through a TestLogger and verify the diagnostic lines + # appear in the expected order: a single x0 line, then NR iter lines. + test_logger = Test.TestLogger(; min_level = Logging.Info) + Logging.with_logger(test_logger) do + solve_power_flow!(data) + end + msgs = [r.message for r in test_logger.logs] + + x0_lines = filter(m -> occursin("x0 (time_step 1)", m), msgs) + iter_lines = filter(m -> occursin(r"NR iter \d+", m), msgs) + @test length(x0_lines) == 1 + @test length(iter_lines) >= 1 + # Each diagnostic line carries ρ, ‖F‖_∞, the binding bus/equation, and κ̂(J). + for line in vcat(x0_lines, iter_lines) + @test occursin("ρ = ", line) + @test occursin("‖F‖_∞ = ", line) + @test occursin(r"at bus \d+ \(P\)|at bus \d+ \(Q\)", line) + @test occursin("κ̂(J) = ", line) + end +end + +@testset "test per-iteration monitor under TrustRegionACPowerFlow" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{TrustRegionACPowerFlow}(; + correct_bustypes = true, + compute_fixed_point_spectral_radius = true, + ) + data = PowerFlowData(pf, sys) + + test_logger = Test.TestLogger(; min_level = Logging.Info) + Logging.with_logger(test_logger) do + solve_power_flow!(data) + end + msgs = [r.message for r in test_logger.logs] + iter_lines = filter(m -> occursin(r"TR iter \d+", m), msgs) + @test length(iter_lines) >= 1 +end From cbb56f5c4f364a6d3bfa7ba71248a3d7b0a75b84 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 8 Apr 2026 08:50:14 -0600 Subject: [PATCH 02/14] address copilot's comments Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fixed_point_spectral_radius.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/fixed_point_spectral_radius.jl b/src/fixed_point_spectral_radius.jl index f2cd5c4eb..c54e7e7b9 100644 --- a/src/fixed_point_spectral_radius.jl +++ b/src/fixed_point_spectral_radius.jl @@ -12,9 +12,8 @@ where ``F_k`` is the ``k``-th component of the AC power flow residual at the cur state stored in `data` (i.e., `data.bus_magnitude` and `data.bus_angles`). This is matrix-free: the residual's 3-tensor of second derivatives is never -materialized. State-vector indexing follows the PowerFlows.jl convention: -`x[2i-1] = Vm` and `x[2i] = Va` for PQ buses, with the same column positions used -for variable lookups via `is_PQ` / `has_θ` gates. +materialized. See [`_create_jacobian_matrix_structure`](@ref) for the +bus-type-dependent state vector convention. LCC state variables are not yet supported; the returned vector covers only the `2 * num_buses` bus state entries. @@ -147,7 +146,7 @@ end """ compute_fixed_point_spectral_radius(data, time_step; x0, tol, maxiter, krylovdim) - -> (ρ::Float64, info) + -> (ρ::Float64, info, condest::Float64) Estimate the spectral radius of the Jacobian of the Newton fixed-point map @@ -166,13 +165,15 @@ Differentiating ``g`` and using ``\\partial F / \\partial x = J`` gives the identity ```math -G := \\frac{\\partial g}{\\partial x} = J^{-1} \\sum_k H_k u_k, +(G v)_k = \\bigl(J^{-1} w\\bigr)_k, +\\qquad w_k = v^\\top H_k u, \\qquad u = J^{-1} F(x) ``` where ``H_k`` is the Hessian of the ``k``-th residual component. ``G`` is never -materialized: the matvec ``G v`` is computed as ``J^{-1} \\cdot (\\sum_k H_k u_k v)`` -via two `KLU` back-solves and one [`acpf_hvvp`](@ref) call. KrylovKit's +materialized: the matvec ``G v`` is computed as one [`acpf_hvvp`](@ref) call +(producing ``w``) followed by two `KLU` back-solves (one for ``u``, one for +``J^{-1} w``). KrylovKit's matrix-free Lanczos / Arnoldi solver then extracts the eigenvalue of largest magnitude. @@ -234,7 +235,8 @@ function _fixed_point_spectral_radius!( F = KLU.klu(jac.Jv) u = F \ copy(residual.Rv) matvec(v::AbstractVector) = F \ acpf_hvvp(data, time_step, v, u) - v_init = randn(n) + # Deterministic init for reproducibility across runs / CI logs. + v_init = ones(Float64, n) ./ sqrt(n) vals, _, info = KrylovKit.eigsolve( matvec, v_init, 1, :LM; tol = tol, maxiter = maxiter, krylovdim = krylovdim, From 3884b4cbb588c0d568371f15a7f4990031bde843 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Mon, 11 May 2026 11:21:38 -0600 Subject: [PATCH 03/14] robust homotopy: respect override x0 --- src/RobustHomotopy/homotopy_hessian.jl | 11 ++++++++-- src/RobustHomotopy/robust_homotopy_method.jl | 3 ++- test/test_robust_power_flow.jl | 22 ++++++++++++++++++++ 3 files changed, 33 insertions(+), 3 deletions(-) diff --git a/src/RobustHomotopy/homotopy_hessian.jl b/src/RobustHomotopy/homotopy_hessian.jl index 113ee443e..37977b212 100644 --- a/src/RobustHomotopy/homotopy_hessian.jl +++ b/src/RobustHomotopy/homotopy_hessian.jl @@ -64,8 +64,15 @@ function gradient_value!(grad::Vector{Float64}, return grad end -function homotopy_x0(data::ACPowerFlowData, time_step::Int) - x = calculate_x0(data, time_step) +function homotopy_x0(data::ACPowerFlowData, time_step::Int; + x0::Union{Vector{Float64}, Nothing} = nothing, +) + x = if OVERRIDE_x0 && !isnothing(x0) + @warn "Overriding initial guess x0." + copy(x0) + else + calculate_x0(data, time_step) + end for (bus_ix, bt) in enumerate(get_bus_type(data)[:, time_step]) # PERF: allocating if bt == PSY.ACBusTypes.PQ x[2 * bus_ix - 1] = 1.0 diff --git a/src/RobustHomotopy/robust_homotopy_method.jl b/src/RobustHomotopy/robust_homotopy_method.jl index 62b51ebed..c3193e106 100644 --- a/src/RobustHomotopy/robust_homotopy_method.jl +++ b/src/RobustHomotopy/robust_homotopy_method.jl @@ -2,10 +2,11 @@ function _newton_power_flow(pf::ACPolarPowerFlow{<:RobustHomotopyPowerFlow}, data::ACPowerFlowData, time_step::Int64; Δt_k::Float64 = DEFAULT_Δt_k, + x0::Union{Vector{Float64}, Nothing} = nothing, _ignored..., ) homHess = HomotopyHessian(data, time_step) - x = homotopy_x0(data, time_step) + x = homotopy_x0(data, time_step; x0 = x0) t_k = 0.0 # the sparse structure of the Hessian is different at t_k = 0.0 and t_k > 0.0 diff --git a/test/test_robust_power_flow.jl b/test/test_robust_power_flow.jl index 651e203d7..f563246f6 100644 --- a/test/test_robust_power_flow.jl +++ b/test/test_robust_power_flow.jl @@ -35,3 +35,25 @@ end @test isapprox(data_nr.bus_angles, data_hom.bus_angles; atol = 1e-4) @test isapprox(data_nr.bus_magnitude, data_hom.bus_magnitude; atol = 1e-6) end + +if PF.OVERRIDE_x0 + @testset "robust homotopy respects override_x0" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf_hom = ACPowerFlow{PF.RobustHomotopyPowerFlow}() + data = PowerFlowData(pf_hom, sys) + + x0 = PF.calculate_x0(data, 1) + # perturb angles so the override is detectable + for (i, bt) in enumerate(PF.get_bus_type(data)[:, 1]) + if bt == PSY.ACBusTypes.PQ + x0[2 * i] *= 1.05 # small angle perturbation + end + end + + @test_logs (:warn, r"Overriding initial guess x0.*") match_mode = :any solve_power_flow!( + data; + pf = pf_hom, + x0 = x0, + ) + end +end From 5ce1dad94be0435298d4ba5df43f05afc206a880 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 19 May 2026 10:00:02 -0600 Subject: [PATCH 04/14] post-rebase fixes --- src/fixed_point_spectral_radius.jl | 7 +------ test/test_fixed_point_spectral_radius.jl | 7 +------ 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/src/fixed_point_spectral_radius.jl b/src/fixed_point_spectral_radius.jl index c54e7e7b9..bc5a3b1c1 100644 --- a/src/fixed_point_spectral_radius.jl +++ b/src/fixed_point_spectral_radius.jl @@ -201,12 +201,7 @@ function compute_fixed_point_spectral_radius( end residual = ACPowerFlowResidual(data, time_step) - jac = ACPowerFlowJacobian( - data, - residual.bus_slack_participation_factors, - residual.subnetworks, - time_step, - ) + jac = ACPowerFlowJacobian(residual, time_step) x = isnothing(x0) ? calculate_x0(data, time_step) : copy(x0) residual(x, time_step) jac(time_step) diff --git a/test/test_fixed_point_spectral_radius.jl b/test/test_fixed_point_spectral_radius.jl index 8f2cc1820..85f15f300 100644 --- a/test/test_fixed_point_spectral_radius.jl +++ b/test/test_fixed_point_spectral_radius.jl @@ -63,12 +63,7 @@ end time_step = 1 residual = PF.ACPowerFlowResidual(data, time_step) - jac = PF.ACPowerFlowJacobian( - data, - residual.bus_slack_participation_factors, - residual.subnetworks, - time_step, - ) + jac = PF.ACPowerFlowJacobian(residual, time_step) x0 = PF.calculate_x0(data, time_step) n = length(x0) From 8f433c4cc952c6f95d07a05028fc98c37ed0fe3f Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 19 May 2026 10:01:50 -0600 Subject: [PATCH 05/14] test RH w dist slack; remove PF from test env. --- test/Project.toml | 1 - test/test_distributed_slack.jl | 7 ++----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 8db34ad89..abff44511 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -13,7 +13,6 @@ JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" PROPACK = "b169e327-5944-5131-97a6-5d3d3f0a476a" -PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" diff --git a/test/test_distributed_slack.jl b/test/test_distributed_slack.jl index 71cf886b4..2e7d74e9d 100644 --- a/test/test_distributed_slack.jl +++ b/test/test_distributed_slack.jl @@ -298,10 +298,8 @@ end end @testset "AC PF with headroom-proportional distributed slack" begin - for ACSolver in - filter(x -> !(x in (RobustHomotopyPowerFlow,)), AC_SOLVERS_TO_TEST) + for ACSolver in AC_SOLVERS_TO_TEST @testset "$(ACSolver)" begin - ACSolver == RobustHomotopyPowerFlow && continue sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") # Introduce active power imbalance so slack distribution is nontrivial @@ -344,8 +342,7 @@ end end @testset "Headroom-proportional slack: DataFrame and generator-level redistribution" begin - for ACSolver in - filter(x -> !(x in (RobustHomotopyPowerFlow,)), AC_SOLVERS_TO_TEST) + for ACSolver in AC_SOLVERS_TO_TEST @testset "$(ACSolver)" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") From a5ae37d987f9a604e15cb4e9eb18ae052e9d948e Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 19 May 2026 14:53:47 -0600 Subject: [PATCH 06/14] override x0 kwarg eliminated --- src/RobustHomotopy/homotopy_hessian.jl | 4 ++-- test/test_robust_power_flow.jl | 32 ++++++++++++-------------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/RobustHomotopy/homotopy_hessian.jl b/src/RobustHomotopy/homotopy_hessian.jl index 37977b212..293c2349c 100644 --- a/src/RobustHomotopy/homotopy_hessian.jl +++ b/src/RobustHomotopy/homotopy_hessian.jl @@ -67,8 +67,8 @@ end function homotopy_x0(data::ACPowerFlowData, time_step::Int; x0::Union{Vector{Float64}, Nothing} = nothing, ) - x = if OVERRIDE_x0 && !isnothing(x0) - @warn "Overriding initial guess x0." + x = if !isnothing(x0) + @warn "Using caller-provided x0; skipping improve_x0." copy(x0) else calculate_x0(data, time_step) diff --git a/test/test_robust_power_flow.jl b/test/test_robust_power_flow.jl index f563246f6..b2b9377f5 100644 --- a/test/test_robust_power_flow.jl +++ b/test/test_robust_power_flow.jl @@ -36,24 +36,22 @@ end @test isapprox(data_nr.bus_magnitude, data_hom.bus_magnitude; atol = 1e-6) end -if PF.OVERRIDE_x0 - @testset "robust homotopy respects override_x0" begin - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") - pf_hom = ACPowerFlow{PF.RobustHomotopyPowerFlow}() - data = PowerFlowData(pf_hom, sys) +@testset "robust homotopy respects override_x0" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf_hom = ACPowerFlow{PF.RobustHomotopyPowerFlow}() + data = PowerFlowData(pf_hom, sys) - x0 = PF.calculate_x0(data, 1) - # perturb angles so the override is detectable - for (i, bt) in enumerate(PF.get_bus_type(data)[:, 1]) - if bt == PSY.ACBusTypes.PQ - x0[2 * i] *= 1.05 # small angle perturbation - end + x0 = PF.calculate_x0(data, 1) + # perturb angles so the override is detectable + for (i, bt) in enumerate(PF.get_bus_type(data)[:, 1]) + if bt == PSY.ACBusTypes.PQ + x0[2 * i] *= 1.05 # small angle perturbation end - - @test_logs (:warn, r"Overriding initial guess x0.*") match_mode = :any solve_power_flow!( - data; - pf = pf_hom, - x0 = x0, - ) end + + @test_logs (:warn, r"Using caller-provided x0.*") match_mode = :any solve_power_flow!( + data; + pf = pf_hom, + x0 = x0, + ) end From 474e0fa99d3ec429113ccd6e575f1c572ab2295a Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 20 May 2026 10:23:08 -0600 Subject: [PATCH 07/14] add an adaptive method --- src/PowerFlows.jl | 1 + src/levenberg-marquardt.jl | 38 +++++- src/power_flow_method.jl | 229 +++++++++++++++++++++++++++++++++++-- src/power_flow_types.jl | 23 +++- 4 files changed, 276 insertions(+), 15 deletions(-) diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index b4885ee4c..43b37a070 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -7,6 +7,7 @@ export NewtonRaphsonACPowerFlow export TrustRegionACPowerFlow export LevenbergMarquardtACPowerFlow export RobustHomotopyPowerFlow +export AdaptiveACPowerFlow export ACPolarPowerFlow export ACRectangularPowerFlow export ACMixedPowerFlow diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index c7860a250..990dec1db 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -144,10 +144,11 @@ function _newton_power_flow( pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol i = 0 + floor_reached = false if !converged use_scaling = something(marquardt_scaling, _default_marquardt_scaling(pf)) ws = LMWorkspace(J.Jv; marquardt_scaling = use_scaling) - converged, i = _run_power_flow_method( + converged, i, floor_reached = _run_power_flow_method( time_step, x0, residual, @@ -160,7 +161,8 @@ function _newton_power_flow( # (or is the already-converged initial state if the loop was skipped). _finalize_formulation!(pf, data, x0, residual, time_step) return _finalize_power_flow( - converged, i, "LevenbergMarquardtACPowerFlow", residual, data, J.Jv, time_step) + converged, i, "LevenbergMarquardtACPowerFlow", residual, data, J.Jv, time_step; + floor_reached) end function _run_power_flow_method( @@ -173,6 +175,9 @@ function _run_power_flow_method( maxIterations::Int = DEFAULT_NR_MAX_ITER, tol::Float64 = DEFAULT_NR_TOL, λ_0::Float64 = DEFAULT_λ_0, + monitor_ρ::Bool = get_compute_fixed_point_spectral_radius(J.data), + F_norm_init::Union{Float64, Nothing} = nothing, + bail_at_floor::Bool = false, _ignored..., ) μ::Float64 = λ_0 @@ -182,12 +187,35 @@ function _run_power_flow_method( resSize = dot(residual.Rv, residual.Rv) linf = norm(residual.Rv, Inf) @debug "initially: sum of squares $(siground(resSize)), L ∞ norm $(siground(linf)), λ = $λ" + # Problem-scale reference for the backward-stability floor estimate. + # First Newton step has backward error ≈ u·κ·‖F_init‖, and that floor persists + # Clamped at 1.0 pu so warm-starts (tiny ‖F_init‖) don't collapse the estimate. + # The adaptive driver passes its own pre-TR ‖F‖_init for a more conservative reference. + F_scale = max(1.0, something(F_norm_init, linf)) + floor_reached = false + F_window = Float64[] while i < maxIterations && !converged && isfinite(λ) && μ < DEFAULT_μ_MAX λ, μ = update_damping_factor!(x, residual, J, μ, time_step, ws) - converged = isfinite(λ) && norm(residual.Rv, Inf) < tol + F_inf = norm(residual.Rv, Inf) + if monitor_ρ + ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("LM iter $i", ρ, residual.Rv, condest) + hit_now, floor_est = _check_numerical_floor!( + F_window, F_inf, condest, F_scale) + if hit_now && !floor_reached + @info "LM hit numerical floor: ‖F‖_∞ = $(siground(F_inf)), " * + "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" + floor_reached = true + end + if floor_reached && bail_at_floor + i += 1 + break + end + end + converged = isfinite(λ) && F_inf < tol i += 1 end - if !converged + if !converged && !floor_reached if !isfinite(λ) @error "λ is not finite ($(λ))" elseif μ >= DEFAULT_μ_MAX @@ -197,7 +225,7 @@ function _run_power_flow_method( end end - return converged, i + return converged, i, floor_reached end # LM implementation based on standard Levenberg-Marquardt method. diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 7adfd1b57..e8c6c2b12 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1,6 +1,29 @@ """Format a per-iteration diagnostic line: ρ, ‖F‖_∞ and the bus/equation where that infinity-norm is attained, and a κ̂(J) condition estimate. All numerics rounded to 4 significant figures.""" +"""Update a rolling 3-iter window of `‖F‖_∞` and decide whether the residual +has hit the backward-stability floor `κ̂·ε·F_scale`. Returns `true` only when +the current residual is within `10×` of the estimated floor AND the last 3 +iterates are within 10% of each other (genuine stagnation, not transient). +Mutates `F_window` in place. The caller is responsible for logging.""" +function _check_numerical_floor!( + F_window::Vector{Float64}, + F_inf::Float64, + condest::Float64, + F_scale::Float64, +) + push!(F_window, F_inf) + length(F_window) > 3 && popfirst!(F_window) + floor_est = condest * eps() * F_scale + if F_inf <= 10 * floor_est && length(F_window) == 3 + F_lo, F_hi = extrema(F_window) + if F_hi - F_lo < 0.1 * F_hi + return true, floor_est + end + end + return false, floor_est +end + function _log_diagnostics( label::String, ρ::Float64, @@ -598,12 +621,17 @@ function _run_power_flow_method(time_step::Int, validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, iwamoto::Bool = false, + F_norm_init::Union{Float64, Nothing} = nothing, + bail_at_floor::Bool = false, _ignored..., # absorb unknown keys from caller without error ) validate_vms = validate_voltage_magnitudes i, converged = 1, false consecutive_reverts = 0 monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) + F_scale = max(1.0, something(F_norm_init, norm(residual.Rv, Inf))) + floor_reached = false + F_window = Float64[] while i < maxIterations && !converged if iwamoto made_progress = _iwamoto_step( @@ -644,13 +672,24 @@ function _run_power_flow_method(time_step::Int, if monitor_ρ ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) _log_diagnostics("NR iter $i", ρ, residual.Rv, condest) + hit_now, floor_est = _check_numerical_floor!( + F_window, norm(residual.Rv, Inf), condest, F_scale) + if hit_now && !floor_reached + @info "NR hit numerical floor: ‖F‖_∞ = " * + "$(siground(norm(residual.Rv, Inf))), " * + "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" + floor_reached = true + end + if floor_reached && bail_at_floor + break + end end converged = norm(residual.Rv, Inf) < tol if !converged i += 1 end end - return converged, i + return converged, i, floor_reached end """Runs the full `TrustRegionNRMethod`. @@ -679,6 +718,8 @@ function _run_power_flow_method(time_step::Int, iwamoto_fallback::Bool = DEFAULT_IWAMOTO_FALLBACK, validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, + F_norm_init::Union{Float64, Nothing} = nothing, + bail_at_floor::Bool = false, _ignored..., # absorb unknown keys from caller without error ) validate_vms = validate_voltage_magnitudes @@ -704,6 +745,9 @@ function _run_power_flow_method(time_step::Int, @debug "initially: sum of squares $(siground(residualSize)), L ∞ norm $(siground(linf)), Δ $(siground(delta))" monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) + F_scale = max(1.0, something(F_norm_init, linf)) + floor_reached = false + F_window = Float64[] while i < maxIterations && !converged delta = _trust_region_step( time_step, @@ -726,17 +770,175 @@ function _run_power_flow_method(time_step::Int, if monitor_ρ ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) _log_diagnostics("TR iter $i", ρ, residual.Rv, condest) + hit_now, floor_est = _check_numerical_floor!( + F_window, norm(residual.Rv, Inf), condest, F_scale) + if hit_now && !floor_reached + @info "TR hit numerical floor: ‖F‖_∞ = " * + "$(siground(norm(residual.Rv, Inf))), " * + "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" + floor_reached = true + end + if floor_reached && bail_at_floor + break + end end converged = norm(residual.Rv, Inf) < tol if !converged i += 1 end end - return converged, i + return converged, i, floor_reached +end + +# Hardcoded escalation thresholds for AdaptiveACPowerFlow. Calibrated against +# CATS multi-period traces: TR's diverging time steps all had two consecutive +# iterates with ρ ≥ 1 by iter 3–4, while converging time steps never crossed +# ρ ≥ 1 after iter 0. +const ADAPTIVE_ρ_THRESHOLD = 1.0 +const ADAPTIVE_ρ_CONSECUTIVE = 2 + +"""Driver for the AdaptiveACPowerFlow method: runs Trust Region with a +spectral-radius watchdog, and escalates to Levenberg-Marquardt (reusing the +current iterate as the starting point) if `ρ ≥ 1` for +`$ADAPTIVE_ρ_CONSECUTIVE` consecutive iterations.""" +function _newton_power_flow( + pf::AbstractACPowerFlow{AdaptiveACPowerFlow}, + data::ACPowerFlowData, + time_step::Int64; + tol::Float64 = DEFAULT_NR_TOL, + maxIterations::Int = DEFAULT_NR_MAX_ITER, + validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, + vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, + factor::Float64 = DEFAULT_TRUST_REGION_FACTOR, + eta::Float64 = DEFAULT_TRUST_REGION_ETA, + autoscale::Bool = DEFAULT_AUTOSCALE, + iwamoto_fallback::Bool = DEFAULT_IWAMOTO_FALLBACK, + λ_0::Float64 = DEFAULT_λ_0, + marquardt_scaling::Union{Bool, Nothing} = nothing, + x0::Union{Vector{Float64}, Nothing} = nothing, + extra_kwargs..., +) + init_kwargs = if isnothing(x0) + (; validate_voltage_magnitudes, vm_validation_range) + else + (; validate_voltage_magnitudes, vm_validation_range, x0) + end + residual, J, x0_init = initialize_power_flow_variables( + pf, data, time_step; init_kwargs...) + converged = norm(residual.Rv, Inf) < tol + F_norm_init = norm(residual.Rv, Inf) + + ρ_x0, _, condest_x0 = _fixed_point_spectral_radius!(data, residual, J, time_step) + _log_diagnostics("x0 (time_step $time_step)", ρ_x0, residual.Rv, condest_x0) + + i = 0 + x_final = x0_init + escalated = false + floor_reached = false + if !converged + linSolveCache = KLULinSolveCache(J.Jv) + symbolic_factor!(linSolveCache, J.Jv) + stateVector = StateVectorCache(x0_init, residual.Rv) + + if autoscale + for k in 1:length(stateVector.x) + stateVector.d[k] = norm(view(J.Jv, :, k)) + if stateVector.d[k] == 0.0 + stateVector.d[k] = 1.0 + end + end + end + + delta = norm(stateVector.x) > 0 ? factor * norm(stateVector.x) : factor + delta_max = DEFAULT_TRUST_REGION_DELTA_MAX_FACTOR * delta + consec_bad_ρ = 0 + F_scale_tr = max(1.0, F_norm_init) + F_window_tr = Float64[] + + # ---- Phase 1: Trust Region with spectral-radius watchdog ---- + while i < maxIterations && !converged && !escalated + delta = _trust_region_step( + time_step, + stateVector, + linSolveCache, + residual, + J, + delta, + delta_max, + eta, + autoscale, + iwamoto_fallback, + ) + validate_voltage_magnitudes && _validate_state_magnitudes( + residual, + stateVector.x, + vm_validation_range, + i, + ) + ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) + _log_diagnostics("Adaptive TR iter $i", ρ, residual.Rv, condest) + + floor_reached, floor_est = _check_numerical_floor!( + F_window_tr, norm(residual.Rv, Inf), condest, F_scale_tr) + if floor_reached + @info "Adaptive TR hit numerical floor: ‖F‖_∞ = " * + "$(siground(norm(residual.Rv, Inf))), est. floor " * + "κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est)); stopping " * + "before escalation" + break + end + + consec_bad_ρ = ρ >= ADAPTIVE_ρ_THRESHOLD ? consec_bad_ρ + 1 : 0 + if consec_bad_ρ >= ADAPTIVE_ρ_CONSECUTIVE + @info "Adaptive solver: ρ ≥ $(ADAPTIVE_ρ_THRESHOLD) for " * + "$(consec_bad_ρ) consecutive iters; escalating to " * + "LevenbergMarquardt at iter $i" + escalated = true + break + end + converged = norm(residual.Rv, Inf) < tol + if !converged + i += 1 + end + end + x_final = stateVector.x + + # ---- Phase 2: Levenberg-Marquardt handoff (reuses current iterate) ---- + if escalated && !converged + use_scaling = something(marquardt_scaling, _default_marquardt_scaling(pf)) + ws = LMWorkspace(J.Jv; marquardt_scaling = use_scaling) + lm_max = max(0, maxIterations - i) + lm_converged, lm_iters, lm_floor = _run_power_flow_method( + time_step, + x_final, + residual, + J, + ws; + tol, + maxIterations = lm_max, + λ_0, + monitor_ρ = true, + F_norm_init, + bail_at_floor = true, + extra_kwargs..., + ) + converged = lm_converged + floor_reached = lm_floor + i += lm_iters + end + end + _finalize_formulation!(pf, data, x_final, residual, time_step) + return _finalize_power_flow( + converged, i, "AdaptiveACPowerFlow", residual, data, J.Jv, time_step; + floor_reached) end """Log final residual, report convergence, compute optional post-processing factors, -and return `true`/`false`. Shared by all AC power flow drivers.""" +and return `true`/`false`. Shared by all AC power flow drivers. When +`floor_reached = true`, emits a distinct message indicating the solver stopped +because the residual hit the κ̂·ε·max(1,‖F‖_init) backward-stability floor (not +convergence to `tol`); the return value is still `false` since the requested +tolerance was not met.""" function _finalize_power_flow( converged::Bool, i::Int, @@ -744,7 +946,8 @@ function _finalize_power_flow( residual::Union{ACPowerFlowResidual, ACRectangularCIResidual, ACMixedCPBResidual}, data::ACPowerFlowData, Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, - time_step::Int64, + time_step::Int64; + floor_reached::Bool = false, ) @info("Final residual size: $(norm(residual.Rv, 2)) L2, $(norm(residual.Rv, Inf)) L∞.") if converged @@ -757,7 +960,16 @@ function _finalize_power_flow( end return true end - @error("The $solver_name solver failed to converge after $i iterations.") + if floor_reached + @warn( + "The $solver_name solver stopped at the numerical floor after $i " * + "iterations: the residual is limited by κ̂·ε·max(1,‖F‖_init) " * + "backward stability, not the requested tolerance. " * + "Reporting non-convergence." + ) + else + @error("The $solver_name solver failed to converge after $i iterations.") + end return false end @@ -842,11 +1054,12 @@ function _newton_power_flow( i = 0 x_final = x0_init + floor_reached = false if !converged linSolveCache = KLULinSolveCache(J.Jv) symbolic_factor!(linSolveCache, J.Jv) stateVector = StateVectorCache(x0_init, residual.Rv) - converged, i = _run_power_flow_method( + converged, i, floor_reached = _run_power_flow_method( time_step, stateVector, linSolveCache, @@ -868,5 +1081,7 @@ function _newton_power_flow( x_final = stateVector.x end _finalize_formulation!(pf, data, x_final, residual, time_step) - return _finalize_power_flow(converged, i, string(T), residual, data, J.Jv, time_step) + return _finalize_power_flow( + converged, i, string(T), residual, data, J.Jv, time_step; + floor_reached) end diff --git a/src/power_flow_types.jl b/src/power_flow_types.jl index 2345bb878..99cb9695a 100644 --- a/src/power_flow_types.jl +++ b/src/power_flow_types.jl @@ -110,6 +110,21 @@ See also: [`ACPowerFlow`](@ref). """ struct RobustHomotopyPowerFlow <: ACPowerFlowSolverType end +""" + AdaptiveACPowerFlow <: ACPowerFlowSolverType + +An [`ACPowerFlowSolverType`](@ref) that switches strategies during the solve +based on online diagnostics: fixed-point spectral radius `ρ`, condition +estimate `κ̂(J)`, and residual `‖F‖`. Starts with a fast inner solver +(Newton-Raphson / Trust Region) and escalates to a more robust method +(Robust Homotopy) when the diagnostics indicate the inner solver is +diverging or stuck in a singular neighborhood. + +Polar formulation only. See also: [`ACPowerFlow`](@ref), +[`RobustHomotopyPowerFlow`](@ref). +""" +struct AdaptiveACPowerFlow <: ACPowerFlowSolverType end + """ ACPowerFlow{ACSolver}(; kwargs...) where {ACSolver <: ACPowerFlowSolverType} ACPowerFlow(; kwargs...) @@ -357,12 +372,13 @@ function ACRectangularPowerFlow{ACSolver}(; if ACSolver <: Union{ RobustHomotopyPowerFlow, GradientDescentACPowerFlow, + AdaptiveACPowerFlow, } throw( ArgumentError( "$(ACSolver) is not supported by ACRectangularPowerFlow. " * - "Robust Homotopy and Gradient Descent operate on the polar " * - "formulation only. Use ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}, " * + "Robust Homotopy, Gradient Descent, and Adaptive operate on the " * + "polar formulation only. Use ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}, " * "{TrustRegionACPowerFlow}, or {LevenbergMarquardtACPowerFlow}, " * "or run the solver on ACPolarPowerFlow.", ), @@ -466,11 +482,12 @@ function ACMixedPowerFlow{ACSolver}(; if ACSolver <: Union{ RobustHomotopyPowerFlow, GradientDescentACPowerFlow, + AdaptiveACPowerFlow, } throw( ArgumentError( "$(ACSolver) is not supported by ACMixedPowerFlow. " * - "Robust Homotopy and Gradient Descent do not operate on the " * + "Robust Homotopy, Gradient Descent, and Adaptive do not operate on the " * "mixed current-power formulation. Use " * "ACMixedPowerFlow{NewtonRaphsonACPowerFlow}, " * "{TrustRegionACPowerFlow}, or {LevenbergMarquardtACPowerFlow}, " * From bd2c1e6a2a4f58fa4fdc75aba4f00176dfbab4b7 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 20 May 2026 11:20:32 -0600 Subject: [PATCH 08/14] add back log diagnostics helper; make iteration not reset when switching strategies --- src/levenberg-marquardt.jl | 3 ++- src/power_flow_method.jl | 31 ++++++++++++++++--------------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index 990dec1db..3970fea44 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -178,6 +178,7 @@ function _run_power_flow_method( monitor_ρ::Bool = get_compute_fixed_point_spectral_radius(J.data), F_norm_init::Union{Float64, Nothing} = nothing, bail_at_floor::Bool = false, + iter_offset::Int = 0, _ignored..., ) μ::Float64 = λ_0 @@ -199,7 +200,7 @@ function _run_power_flow_method( F_inf = norm(residual.Rv, Inf) if monitor_ρ ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("LM iter $i", ρ, residual.Rv, condest) + _log_diagnostics("LM iter $(i + iter_offset)", ρ, residual.Rv, condest) hit_now, floor_est = _check_numerical_floor!( F_window, F_inf, condest, F_scale) if hit_now && !floor_reached diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index e8c6c2b12..1bc94c865 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1,6 +1,21 @@ """Format a per-iteration diagnostic line: ρ, ‖F‖_∞ and the bus/equation where that infinity-norm is attained, and a κ̂(J) condition estimate. All numerics rounded to 4 significant figures.""" +function _log_diagnostics( + label::String, + ρ::Float64, + Rv::AbstractVector{Float64}, + condest::Float64, +) + abs_max, ix = findmax(abs, Rv) + pow_type = ix % 2 == 1 ? "P" : "Q" + bus_ix = div(ix + 1, 2) + sf(x) = round(x; sigdigits = 4) + @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at bus $bus_ix ($pow_type), " * + "κ̂(J) = $(sf(condest))" + return +end + """Update a rolling 3-iter window of `‖F‖_∞` and decide whether the residual has hit the backward-stability floor `κ̂·ε·F_scale`. Returns `true` only when the current residual is within `10×` of the estimated floor AND the last 3 @@ -24,21 +39,6 @@ function _check_numerical_floor!( return false, floor_est end -function _log_diagnostics( - label::String, - ρ::Float64, - Rv::AbstractVector{Float64}, - condest::Float64, -) - abs_max, ix = findmax(abs, Rv) - pow_type = ix % 2 == 1 ? "P" : "Q" - bus_ix = div(ix + 1, 2) - sf(x) = round(x; sigdigits = 4) - @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at bus $bus_ix ($pow_type), " * - "κ̂(J) = $(sf(condest))" - return -end - """Cache for non-linear methods. # Fields @@ -920,6 +920,7 @@ function _newton_power_flow( monitor_ρ = true, F_norm_init, bail_at_floor = true, + iter_offset = i + 1, extra_kwargs..., ) converged = lm_converged From 1166f9c1fa3e1bd8b39fdc404021a7931f44c53a Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Fri, 22 May 2026 13:24:58 -0600 Subject: [PATCH 09/14] second pass at adaptive method --- src/PowerFlows.jl | 2 + src/ac_power_flow_jacobian.jl | 6 +- src/definitions.jl | 31 ++++ src/levenberg-marquardt.jl | 66 ++++--- src/power_flow_method.jl | 278 ++++++++++++++++++----------- src/stagnation_diagnostics.jl | 319 ++++++++++++++++++++++++++++++++++ 6 files changed, 575 insertions(+), 127 deletions(-) create mode 100644 src/stagnation_diagnostics.jl diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 43b37a070..ca3cacf76 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -8,6 +8,7 @@ export TrustRegionACPowerFlow export LevenbergMarquardtACPowerFlow export RobustHomotopyPowerFlow export AdaptiveACPowerFlow +export ACPowerFlowSolveStatus export ACPolarPowerFlow export ACRectangularPowerFlow export ACMixedPowerFlow @@ -75,6 +76,7 @@ include("mixed_cpb_setup.jl") include("mixed_cpb_power_flow_residual.jl") include("mixed_cpb_power_flow_jacobian.jl") include("fixed_point_spectral_radius.jl") +include("stagnation_diagnostics.jl") include("solve_ac_power_flow.jl") include("power_flow_setup.jl") include("power_flow_method.jl") diff --git a/src/ac_power_flow_jacobian.jl b/src/ac_power_flow_jacobian.jl index 8e46c31c7..3f3ffe31d 100644 --- a/src/ac_power_flow_jacobian.jl +++ b/src/ac_power_flow_jacobian.jl @@ -456,7 +456,7 @@ function _create_jacobian_matrix_structure( end function _set_entries_for_neighbor(::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE}, - Y_from_to::ComplexF32, + Y_from_to::YBUS_ELTYPE, Vm_from::Float64, Vm_to::Float64, θ_from_to::Float64, @@ -481,7 +481,7 @@ function _set_entries_for_neighbor(::SparseArrays.SparseMatrixCSC{Float64, J_IND end function _set_entries_for_neighbor(Jv::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE}, - Y_from_to::ComplexF32, + Y_from_to::YBUS_ELTYPE, Vm_from::Float64, Vm_to::Float64, θ_from_to::Float64, @@ -515,7 +515,7 @@ function _set_entries_for_neighbor(Jv::SparseArrays.SparseMatrixCSC{Float64, J_I end function _set_entries_for_neighbor(Jv::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE}, - Y_from_to::ComplexF32, + Y_from_to::YBUS_ELTYPE, Vm_from::Float64, Vm_to::Float64, θ_from_to::Float64, diff --git a/src/definitions.jl b/src/definitions.jl index 539dd88e5..ba33e17ac 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -67,6 +67,37 @@ const MinMax = NamedTuple{(:min, :max), Tuple{Float64, Float64}} const DEFAULT_VALIDATION_RANGE = (min = 0.5, max = 1.5) # const MAX_INDS_TO_PRINT = 10 +IS.@scoped_enum( + ACPowerFlowSolveStatus, + RUNNING = 0, + CONVERGED = 1, + NON_ROOT_LOCAL_MIN = 2, + NON_ROOT_SADDLE = 3, + NON_ROOT_STATIONARY = 4, + SINGULAR_SUBSPACE = 5, + FOLD = 6, + LIMIT_CYCLE = 7, + STAGNATED_OTHER = 8, + DIVERGED = 9, + FAILED = 10, +) +@doc """ +Concluding state of an AC power flow solve. Used to classify the reason any solver stopped. + +# 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. + """ ACPowerFlowSolveStatus + const FACTS_MODE_MAP = Dict( PSY.FACTSOperationModes.OOS => 0, PSY.FACTSOperationModes.NML => 1, diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index 3970fea44..4d78388d0 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -144,11 +144,12 @@ function _newton_power_flow( pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol i = 0 - floor_reached = false + status::ACPowerFlowSolveStatus = + converged ? ACPowerFlowSolveStatus.CONVERGED : ACPowerFlowSolveStatus.RUNNING if !converged use_scaling = something(marquardt_scaling, _default_marquardt_scaling(pf)) ws = LMWorkspace(J.Jv; marquardt_scaling = use_scaling) - converged, i, floor_reached = _run_power_flow_method( + converged, i, status = _run_power_flow_method( time_step, x0, residual, @@ -162,7 +163,7 @@ function _newton_power_flow( _finalize_formulation!(pf, data, x0, residual, time_step) return _finalize_power_flow( converged, i, "LevenbergMarquardtACPowerFlow", residual, data, J.Jv, time_step; - floor_reached) + status) end function _run_power_flow_method( @@ -176,8 +177,8 @@ function _run_power_flow_method( tol::Float64 = DEFAULT_NR_TOL, λ_0::Float64 = DEFAULT_λ_0, monitor_ρ::Bool = get_compute_fixed_point_spectral_radius(J.data), - F_norm_init::Union{Float64, Nothing} = nothing, - bail_at_floor::Bool = false, + bail_on_stagnation::Bool = false, + detect_stagnation::Bool = true, iter_offset::Int = 0, _ignored..., ) @@ -188,45 +189,60 @@ function _run_power_flow_method( resSize = dot(residual.Rv, residual.Rv) linf = norm(residual.Rv, Inf) @debug "initially: sum of squares $(siground(resSize)), L ∞ norm $(siground(linf)), λ = $λ" - # Problem-scale reference for the backward-stability floor estimate. - # First Newton step has backward error ≈ u·κ·‖F_init‖, and that floor persists - # Clamped at 1.0 pu so warm-starts (tiny ‖F_init‖) don't collapse the estimate. - # The adaptive driver passes its own pre-TR ‖F‖_init for a more conservative reference. - F_scale = max(1.0, something(F_norm_init, linf)) - floor_reached = false + status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING F_window = Float64[] + ρ_window = Float64[] + κ_window = Float64[] while i < maxIterations && !converged && isfinite(λ) && μ < DEFAULT_μ_MAX λ, μ = update_damping_factor!(x, residual, J, μ, time_step, ws) F_inf = norm(residual.Rv, Inf) + ρ_now::Union{Float64, Nothing} = nothing + condest_now::Union{Float64, Nothing} = nothing if monitor_ρ - ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("LM iter $(i + iter_offset)", ρ, residual.Rv, condest) - hit_now, floor_est = _check_numerical_floor!( - F_window, F_inf, condest, F_scale) - if hit_now && !floor_reached - @info "LM hit numerical floor: ‖F‖_∞ = $(siground(F_inf)), " * - "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" - floor_reached = true - end - if floor_reached && bail_at_floor - i += 1 - break + ρ_now, _, condest_now = + _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("LM iter $(i + iter_offset)", ρ_now, residual.Rv, condest_now) + end + if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING + kind = _check_stagnation!( + F_window, ρ_window, κ_window, F_inf, ρ_now, condest_now) + if kind === :fixed_point + msg, status = _stagnation_diagnostic( + J.data, time_step, residual.Rv, J.Jv; condest = condest_now) + @info "LM stagnated at fixed point: ‖F‖_∞ = " * + "$(siground(F_inf)), stable across " * + "$(STAGNATION_WINDOW) iterations" * msg + elseif kind === :limit_cycle + msg, status = _limit_cycle_diagnostic( + residual.Rv, ρ_window, κ_window) + @info "LM in limit cycle: ‖F‖_∞ = " * + "$(siground(F_inf)) across " * + "$(STAGNATION_WINDOW) iterations" * msg end end + if status !== ACPowerFlowSolveStatus.RUNNING && bail_on_stagnation + i += 1 + break + end converged = isfinite(λ) && F_inf < tol i += 1 end - if !converged && !floor_reached + if converged + status = ACPowerFlowSolveStatus.CONVERGED + elseif status === ACPowerFlowSolveStatus.RUNNING if !isfinite(λ) @error "λ is not finite ($(λ))" + status = ACPowerFlowSolveStatus.DIVERGED elseif μ >= DEFAULT_μ_MAX @error "The LevenbergMarquardtACPowerFlow damping factor μ hit the cap (DEFAULT_μ_MAX=$(DEFAULT_μ_MAX)) after $i iterations; aborting (likely divergence)." + status = ACPowerFlowSolveStatus.DIVERGED elseif i == maxIterations @error "The LevenbergMarquardtACPowerFlow solver didn't coverge in $maxIterations iterations." + status = ACPowerFlowSolveStatus.FAILED end end - return converged, i, floor_reached + return converged, i, status end # LM implementation based on standard Levenberg-Marquardt method. diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 1bc94c865..5e0840f8b 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -16,29 +16,6 @@ function _log_diagnostics( return end -"""Update a rolling 3-iter window of `‖F‖_∞` and decide whether the residual -has hit the backward-stability floor `κ̂·ε·F_scale`. Returns `true` only when -the current residual is within `10×` of the estimated floor AND the last 3 -iterates are within 10% of each other (genuine stagnation, not transient). -Mutates `F_window` in place. The caller is responsible for logging.""" -function _check_numerical_floor!( - F_window::Vector{Float64}, - F_inf::Float64, - condest::Float64, - F_scale::Float64, -) - push!(F_window, F_inf) - length(F_window) > 3 && popfirst!(F_window) - floor_est = condest * eps() * F_scale - if F_inf <= 10 * floor_est && length(F_window) == 3 - F_lo, F_hi = extrema(F_window) - if F_hi - F_lo < 0.1 * F_hi - return true, floor_est - end - end - return false, floor_est -end - """Cache for non-linear methods. # Fields @@ -314,6 +291,7 @@ function _trust_region_step(time_step::Int, @debug "Trust region step: ρ = $(siground(rho)), η = $(siground(eta)), ||Δx|| = $(siground(norm(stateVector.Δx_proposed)))" step_accepted = false + iwamoto_accepted = false if rho > eta # Successful iteration @debug "Step accepted: sum of squares $(siground(dot(residual.Rv, residual.Rv))), L ∞ norm $(siground(norm(residual.Rv, Inf))), Δ = $(siground(delta)), ||Δx|| = $(siground(norm(stateVector.Δx_proposed)))" @@ -332,7 +310,7 @@ function _trust_region_step(time_step::Int, # expansion logic because rho corresponds to the rejected full step. delta = min(delta / 2, delta_max) @debug "Trust region decreased (Iwamoto fallback accepted): δ $(siground(old_delta)) → $(siground(delta))" - return delta + return delta, true end else stateVector.x .-= stateVector.Δx_proposed @@ -356,7 +334,7 @@ function _trust_region_step(time_step::Int, @debug "Trust region unchanged: δ = $(siground(delta))" end delta = min(delta, delta_max) - return delta + return delta, (step_accepted || iwamoto_accepted) end """Evaluate the Iwamoto objective g(μ) = ‖(1-μ)f₀ + μ²f₁‖² expanded as @@ -621,17 +599,18 @@ function _run_power_flow_method(time_step::Int, validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, iwamoto::Bool = false, - F_norm_init::Union{Float64, Nothing} = nothing, - bail_at_floor::Bool = false, + bail_on_stagnation::Bool = false, + detect_stagnation::Bool = true, _ignored..., # absorb unknown keys from caller without error ) validate_vms = validate_voltage_magnitudes i, converged = 1, false consecutive_reverts = 0 monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) - F_scale = max(1.0, something(F_norm_init, norm(residual.Rv, Inf))) - floor_reached = false + status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING F_window = Float64[] + ρ_window = Float64[] + κ_window = Float64[] while i < maxIterations && !converged if iwamoto made_progress = _iwamoto_step( @@ -669,27 +648,45 @@ function _run_power_flow_method(time_step::Int, vm_validation_range, i, ) + ρ_now::Union{Float64, Nothing} = nothing + condest_now::Union{Float64, Nothing} = nothing if monitor_ρ - ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("NR iter $i", ρ, residual.Rv, condest) - hit_now, floor_est = _check_numerical_floor!( - F_window, norm(residual.Rv, Inf), condest, F_scale) - if hit_now && !floor_reached - @info "NR hit numerical floor: ‖F‖_∞ = " * - "$(siground(norm(residual.Rv, Inf))), " * - "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" - floor_reached = true - end - if floor_reached && bail_at_floor - break + ρ_now, _, condest_now = + _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("NR iter $i", ρ_now, residual.Rv, condest_now) + end + F_inf_now = norm(residual.Rv, Inf) + if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING + kind = _check_stagnation!( + F_window, ρ_window, κ_window, F_inf_now, ρ_now, condest_now) + if kind === :fixed_point + msg, status = _stagnation_diagnostic( + J.data, time_step, residual.Rv, J.Jv; condest = condest_now) + @info "NR stagnated at fixed point: ‖F‖_∞ = " * + "$(siground(F_inf_now)), stable across " * + "$(STAGNATION_WINDOW) iterations" * msg + elseif kind === :limit_cycle + msg, status = _limit_cycle_diagnostic( + residual.Rv, ρ_window, κ_window) + @info "NR in limit cycle: ‖F‖_∞ = " * + "$(siground(F_inf_now)) across " * + "$(STAGNATION_WINDOW) iterations" * msg end end - converged = norm(residual.Rv, Inf) < tol + if status !== ACPowerFlowSolveStatus.RUNNING && bail_on_stagnation + break + end + converged = F_inf_now < tol if !converged i += 1 end end - return converged, i, floor_reached + if converged + status = ACPowerFlowSolveStatus.CONVERGED + elseif status === ACPowerFlowSolveStatus.RUNNING + status = ACPowerFlowSolveStatus.FAILED + end + return converged, i, status end """Runs the full `TrustRegionNRMethod`. @@ -718,8 +715,8 @@ function _run_power_flow_method(time_step::Int, iwamoto_fallback::Bool = DEFAULT_IWAMOTO_FALLBACK, validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, - F_norm_init::Union{Float64, Nothing} = nothing, - bail_at_floor::Bool = false, + bail_on_stagnation::Bool = false, + detect_stagnation::Bool = true, _ignored..., # absorb unknown keys from caller without error ) validate_vms = validate_voltage_magnitudes @@ -745,11 +742,17 @@ function _run_power_flow_method(time_step::Int, @debug "initially: sum of squares $(siground(residualSize)), L ∞ norm $(siground(linf)), Δ $(siground(delta))" monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) - F_scale = max(1.0, something(F_norm_init, linf)) - floor_reached = false + status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING F_window = Float64[] + ρ_window = Float64[] + κ_window = Float64[] + # Cache ρ/κ̂ across iterations so we don't recompute when TR rejects a step + # (x and J unchanged on rejection ⇒ ρ and κ̂ unchanged). The cached values + # still feed the stagnation windows so a stretch of rejections registers. + ρ_cached::Union{Float64, Nothing} = nothing + condest_cached::Union{Float64, Nothing} = nothing while i < maxIterations && !converged - delta = _trust_region_step( + delta, accepted = _trust_region_step( time_step, stateVector, linSolveCache, @@ -768,26 +771,46 @@ function _run_power_flow_method(time_step::Int, i, ) if monitor_ρ - ρ, _, condest = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("TR iter $i", ρ, residual.Rv, condest) - hit_now, floor_est = _check_numerical_floor!( - F_window, norm(residual.Rv, Inf), condest, F_scale) - if hit_now && !floor_reached - @info "TR hit numerical floor: ‖F‖_∞ = " * - "$(siground(norm(residual.Rv, Inf))), " * - "est. floor κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est))" - floor_reached = true + if accepted + ρ_cached, _, condest_cached = + _fixed_point_spectral_radius!(J.data, residual, J, time_step) + _log_diagnostics("TR iter $i", ρ_cached, residual.Rv, condest_cached) + else + @info "TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" end - if floor_reached && bail_at_floor - break + end + F_inf_now = norm(residual.Rv, Inf) + if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING + kind = _check_stagnation!( + F_window, ρ_window, κ_window, F_inf_now, ρ_cached, condest_cached) + if kind === :fixed_point + msg, status = _stagnation_diagnostic( + J.data, time_step, residual.Rv, J.Jv; condest = condest_cached) + @info "TR stagnated at fixed point: ‖F‖_∞ = " * + "$(siground(F_inf_now)), stable across " * + "$(STAGNATION_WINDOW) iterations" * msg + elseif kind === :limit_cycle + msg, status = _limit_cycle_diagnostic( + residual.Rv, ρ_window, κ_window) + @info "TR in limit cycle: ‖F‖_∞ = " * + "$(siground(F_inf_now)) across " * + "$(STAGNATION_WINDOW) iterations" * msg end end - converged = norm(residual.Rv, Inf) < tol + if status !== ACPowerFlowSolveStatus.RUNNING && bail_on_stagnation + break + end + converged = F_inf_now < tol if !converged i += 1 end end - return converged, i, floor_reached + if converged + status = ACPowerFlowSolveStatus.CONVERGED + elseif status === ACPowerFlowSolveStatus.RUNNING + status = ACPowerFlowSolveStatus.FAILED + end + return converged, i, status end # Hardcoded escalation thresholds for AdaptiveACPowerFlow. Calibrated against @@ -816,6 +839,7 @@ function _newton_power_flow( λ_0::Float64 = DEFAULT_λ_0, marquardt_scaling::Union{Bool, Nothing} = nothing, x0::Union{Vector{Float64}, Nothing} = nothing, + detect_stagnation::Bool = true, extra_kwargs..., ) init_kwargs = if isnothing(x0) @@ -826,7 +850,6 @@ function _newton_power_flow( residual, J, x0_init = initialize_power_flow_variables( pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol - F_norm_init = norm(residual.Rv, Inf) ρ_x0, _, condest_x0 = _fixed_point_spectral_radius!(data, residual, J, time_step) _log_diagnostics("x0 (time_step $time_step)", ρ_x0, residual.Rv, condest_x0) @@ -834,7 +857,7 @@ function _newton_power_flow( i = 0 x_final = x0_init escalated = false - floor_reached = false + status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING if !converged linSolveCache = KLULinSolveCache(J.Jv) symbolic_factor!(linSolveCache, J.Jv) @@ -852,12 +875,15 @@ function _newton_power_flow( delta = norm(stateVector.x) > 0 ? factor * norm(stateVector.x) : factor delta_max = DEFAULT_TRUST_REGION_DELTA_MAX_FACTOR * delta consec_bad_ρ = 0 - F_scale_tr = max(1.0, F_norm_init) F_window_tr = Float64[] + ρ_window_tr = Float64[] + κ_window_tr = Float64[] + ρ = NaN + condest = NaN # ---- Phase 1: Trust Region with spectral-radius watchdog ---- while i < maxIterations && !converged && !escalated - delta = _trust_region_step( + delta, accepted = _trust_region_step( time_step, stateVector, linSolveCache, @@ -875,17 +901,37 @@ function _newton_power_flow( vm_validation_range, i, ) - ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) - _log_diagnostics("Adaptive TR iter $i", ρ, residual.Rv, condest) - - floor_reached, floor_est = _check_numerical_floor!( - F_window_tr, norm(residual.Rv, Inf), condest, F_scale_tr) - if floor_reached - @info "Adaptive TR hit numerical floor: ‖F‖_∞ = " * - "$(siground(norm(residual.Rv, Inf))), est. floor " * - "κ̂·ε·max(1,‖F‖_init) ≈ $(siground(floor_est)); stopping " * - "before escalation" - break + if accepted + ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) + _log_diagnostics("Adaptive TR iter $i", ρ, residual.Rv, condest) + else + @info "Adaptive TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" + end + + F_inf_now = norm(residual.Rv, Inf) + ρ_for_window = isfinite(ρ) ? ρ : nothing + κ_for_window = isfinite(condest) ? condest : nothing + if detect_stagnation + kind = _check_stagnation!( + F_window_tr, ρ_window_tr, κ_window_tr, + F_inf_now, ρ_for_window, κ_for_window) + if kind === :fixed_point + msg, status = _stagnation_diagnostic( + J.data, time_step, residual.Rv, J.Jv; condest = condest) + @info "Adaptive TR stagnated at fixed point: ‖F‖_∞ = " * + "$(siground(F_inf_now)), stable across " * + "$(STAGNATION_WINDOW) iterations; stopping before " * + "escalation" * msg + break + elseif kind === :limit_cycle + msg, status = _limit_cycle_diagnostic( + residual.Rv, ρ_window_tr, κ_window_tr) + @info "Adaptive TR in limit cycle: ‖F‖_∞ = " * + "$(siground(F_inf_now)) across " * + "$(STAGNATION_WINDOW) iterations; stopping before " * + "escalation" * msg + break + end end consec_bad_ρ = ρ >= ADAPTIVE_ρ_THRESHOLD ? consec_bad_ρ + 1 : 0 @@ -908,7 +954,7 @@ function _newton_power_flow( use_scaling = something(marquardt_scaling, _default_marquardt_scaling(pf)) ws = LMWorkspace(J.Jv; marquardt_scaling = use_scaling) lm_max = max(0, maxIterations - i) - lm_converged, lm_iters, lm_floor = _run_power_flow_method( + lm_converged, lm_iters, lm_status = _run_power_flow_method( time_step, x_final, residual, @@ -918,28 +964,34 @@ function _newton_power_flow( maxIterations = lm_max, λ_0, monitor_ρ = true, - F_norm_init, - bail_at_floor = true, + bail_on_stagnation = true, + detect_stagnation, iter_offset = i + 1, extra_kwargs..., ) converged = lm_converged - floor_reached = lm_floor + status = lm_status i += lm_iters end end + if converged + status = ACPowerFlowSolveStatus.CONVERGED + elseif status === ACPowerFlowSolveStatus.RUNNING + status = ACPowerFlowSolveStatus.FAILED + end _finalize_formulation!(pf, data, x_final, residual, time_step) return _finalize_power_flow( converged, i, "AdaptiveACPowerFlow", residual, data, J.Jv, time_step; - floor_reached) + status) end """Log final residual, report convergence, compute optional post-processing factors, -and return `true`/`false`. Shared by all AC power flow drivers. When -`floor_reached = true`, emits a distinct message indicating the solver stopped -because the residual hit the κ̂·ε·max(1,‖F‖_init) backward-stability floor (not -convergence to `tol`); the return value is still `false` since the requested -tolerance was not met.""" +and return `true`/`false`. Shared by all AC power flow drivers. `status` is a +fine-grained classification of the concluding state — one of `ACPowerFlowSolveStatus.CONVERGED`, +`ACPowerFlowSolveStatus.NON_ROOT_STATIONARY` (infeasibility), `ACPowerFlowSolveStatus.SINGULAR_SUBSPACE` (trapped on +near-singularity), `ACPowerFlowSolveStatus.FOLD` (curvature dominates), `ACPowerFlowSolveStatus.STAGNATED_OTHER`, or +`ACPowerFlowSolveStatus.FAILED` (iteration limit, no stagnation detected). The function returns +`true` iff `status == ACPowerFlowSolveStatus.CONVERGED`.""" function _finalize_power_flow( converged::Bool, i::Int, @@ -948,11 +1000,11 @@ function _finalize_power_flow( data::ACPowerFlowData, Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, time_step::Int64; - floor_reached::Bool = false, + status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.FAILED, ) @info("Final residual size: $(norm(residual.Rv, 2)) L2, $(norm(residual.Rv, Inf)) L∞.") if converged - @info("The $solver_name solver converged after $i iterations.") + @info("The $solver_name solver converged after $i iterations [status: $status].") if get_calculate_loss_factors(data) _calculate_loss_factors(data, Jv, time_step) end @@ -961,16 +1013,43 @@ function _finalize_power_flow( end return true end - if floor_reached - @warn( - "The $solver_name solver stopped at the numerical floor after $i " * - "iterations: the residual is limited by κ̂·ε·max(1,‖F‖_init) " * - "backward stability, not the requested tolerance. " * - "Reporting non-convergence." - ) + msg = if status === ACPowerFlowSolveStatus.NON_ROOT_LOCAL_MIN + "stopped at a non-root local minimum of ½‖F‖² (positive directional " * + "curvature). The power flow equations are infeasible at this " * + "operating point — consider PV→PQ reactive-limit switching, switched-" * + "shunt regulation, dispatch changes, or topology checks" + elseif status === ACPowerFlowSolveStatus.NON_ROOT_SADDLE + "stopped at a non-root saddle of ½‖F‖² (negative directional " * + "curvature). A different initial guess may find a root; consider " * + "homotopy, a DC-power-flow seed, or a perturbed warm start" + elseif status === ACPowerFlowSolveStatus.NON_ROOT_STATIONARY + "stopped at a non-root stationary point of ½‖F‖² (curvature " * + "inconclusive). The power flow equations are likely infeasible at " * + "this operating point — consider PV→PQ reactive-limit switching, " * + "switched-shunt regulation, dispatch changes, or topology checks" + elseif status === ACPowerFlowSolveStatus.SINGULAR_SUBSPACE + "stopped trapped on a near-singular subspace of J. The system is at " * + "or past a fold/loadability boundary — homotopy/continuation may be " * + "needed to track the solution branch" + elseif status === ACPowerFlowSolveStatus.FOLD + "stopped at a fold of the residual manifold: quadratic curvature " * + "dominates the Newton step, so further linear iteration cannot " * + "reduce ‖F‖. Try homotopy or a more robust formulation" + elseif status === ACPowerFlowSolveStatus.LIMIT_CYCLE + "settled into a period-≥2 limit cycle: the iteration map has no " * + "fixed point in this neighborhood (typically a fold/saddle-node " * + "bifurcation). Point-wise classifiers (local-min/saddle/fold) don't " * + "apply. Try homotopy/continuation or relax the operating point" + elseif status === ACPowerFlowSolveStatus.STAGNATED_OTHER + "stagnated for unclear reasons (‖F‖_∞ flat across 3+ iterates but no " * + "specific diagnostic threshold tripped); inspect the solver log" else - @error("The $solver_name solver failed to converge after $i iterations.") + "failed to converge" end + @warn( + "The $solver_name solver $msg after $i iterations [status: $status]. " * + "Reporting non-convergence." + ) return false end @@ -1055,12 +1134,13 @@ function _newton_power_flow( i = 0 x_final = x0_init - floor_reached = false + status::ACPowerFlowSolveStatus = + converged ? ACPowerFlowSolveStatus.CONVERGED : ACPowerFlowSolveStatus.RUNNING if !converged linSolveCache = KLULinSolveCache(J.Jv) symbolic_factor!(linSolveCache, J.Jv) stateVector = StateVectorCache(x0_init, residual.Rv) - converged, i, floor_reached = _run_power_flow_method( + converged, i, status = _run_power_flow_method( time_step, stateVector, linSolveCache, @@ -1084,5 +1164,5 @@ function _newton_power_flow( _finalize_formulation!(pf, data, x_final, residual, time_step) return _finalize_power_flow( converged, i, string(T), residual, data, J.Jv, time_step; - floor_reached) + status) end diff --git a/src/stagnation_diagnostics.jl b/src/stagnation_diagnostics.jl new file mode 100644 index 000000000..9c4486dee --- /dev/null +++ b/src/stagnation_diagnostics.jl @@ -0,0 +1,319 @@ +# Post-solve stagnation diagnostics. Self-contained: given a Jacobian, residual +# vector, and (for the Hessian terms) `data`/`time_step`, classify the cause of +# stagnation into an `ACPowerFlowSolveStatus` and produce a one-line log +# summary. Hessian-vector-vector products go through `acpf_hvvp`, which is +# polar-only today; non-polar formulations would need their own hvvp +# implementation to plug in here. + +# Classification thresholds. Tuned against CATS multi-period traces; treat as +# illustrative starting points, not load-bearing. +const STAGNATION_STATIONARY_THRESHOLD = 1e-2 # ‖JᵀF‖/(‖J‖_F·‖F‖) below → non-root stationary +const STAGNATION_FOLD_THRESHOLD = 0.5 # ‖½H[Δx,Δx]‖_∞/‖F‖_∞ above → fold +const STAGNATION_SUBSPACE_THRESHOLD = 0.9 # subspace alignment above → trapped on near-singularity +# Default number of bottom singular vectors probed by the subspace alignment. +const STAGNATION_SUBSPACE_K = 3 + +# Stability-gate parameters for `_check_stagnation!`. Two gates run in parallel: +# the **strict** gate (`ρ` and `κ̂` within ~10%) detects a true fixed point of +# the iteration map; the **loose** gate (~2×) detects period-≥2 limit cycles +# where ρ/κ̂ alternate between bounded clusters but never settle. `‖F‖_∞` uses +# the same 10% band for both — it's flat in either case. +const STAGNATION_WINDOW = 5 +const STAGNATION_F_BAND = 1.1 # max/min ‖F‖_∞ in window +const STAGNATION_ρ_BAND_STRICT = 1.1 # ρ within 10% ⇒ fixed point +const STAGNATION_κ_BAND_STRICT = 1.1 # κ̂ within 10% ⇒ fixed point +const STAGNATION_ρ_BAND_LOOSE = 2.0 # ρ within 2× ⇒ limit cycle +const STAGNATION_κ_BAND_LOOSE = 2.0 # κ̂ within 2× ⇒ limit cycle + +"""Update rolling windows of `‖F‖_∞`, `ρ`, and `κ̂(J)` and classify the recent +behavior of the iteration. Returns one of: + +- `:fixed_point` — `‖F‖`, `ρ`, `κ̂` all within tight bands across the window. + The iterate has converged (in the iteration-map sense) to a single point. + Run [`_stagnation_diagnostic`](@ref) to determine the geometric character + (local min / saddle / fold). +- `:limit_cycle` — `‖F‖` within tight band but `ρ` or `κ̂` only within the + loose band. The iterate is in a period-≥2 cycle (alternates between bounded + clusters but never settles). Typical near a fold/saddle-node bifurcation; + point-wise classifiers like gradient/curvature tests don't apply. +- `:none` — not (yet) stagnated. + +Pass `nothing` for `ρ`/`κ` when `compute_fixed_point_spectral_radius = false`; +the `‖F‖_∞` window is always required, and without `ρ`/`κ` data the helper +defaults to `:fixed_point` when `‖F‖` is flat (can't distinguish further). + +Mutates the three windows in place.""" +function _check_stagnation!( + F_window::Vector{Float64}, + ρ_window::Vector{Float64}, + κ_window::Vector{Float64}, + F_inf::Float64, + ρ::Union{Float64, Nothing}, + κ::Union{Float64, Nothing}, +) + _push_capped!(window, x) = begin + push!(window, x) + length(window) > STAGNATION_WINDOW && popfirst!(window) + end + _push_capped!(F_window, F_inf) + if !isnothing(ρ) && isfinite(ρ) + _push_capped!(ρ_window, ρ) + end + if !isnothing(κ) && isfinite(κ) + _push_capped!(κ_window, κ) + end + _stable(window, band) = + length(window) >= STAGNATION_WINDOW && + let (lo, hi) = extrema(window) + hi / max(lo, eps()) < band + end + _stable(F_window, STAGNATION_F_BAND) || return :none + # If we have no spectral data, we can only check F; report fixed_point. + have_ρ = !isempty(ρ_window) + have_κ = !isempty(κ_window) + if !have_ρ && !have_κ + return :fixed_point + end + # Loose gate first — anything failing loose isn't stagnated at all. + have_ρ && + (_stable(ρ_window, STAGNATION_ρ_BAND_LOOSE) || return :none) + have_κ && + (_stable(κ_window, STAGNATION_κ_BAND_LOOSE) || return :none) + # Strict gate: all tracked spectral quantities within tight bands. + strict_ρ = !have_ρ || _stable(ρ_window, STAGNATION_ρ_BAND_STRICT) + strict_κ = !have_κ || _stable(κ_window, STAGNATION_κ_BAND_STRICT) + return (strict_ρ && strict_κ) ? :fixed_point : :limit_cycle +end + +"""Deflated inverse iteration on `JᵀJ` (via the KLU factor) to estimate the +bottom-`k` left-singular vectors of `J`. Returns the per-vector cosine +similarities `|⟨u_i, F⟩| / ‖F‖`, the subspace alignment +`√(Σ |⟨u_i, F⟩|²) / ‖F‖`, and the corresponding right singular vectors `V` +(needed for the directional-curvature test). A subspace alignment near 1 with +individual cosines all moderate is the canonical signature of a multi- +dimensional near-singular subspace. + +`n_starts` runs the procedure with multiple random initial vectors and returns +the result with the largest subspace alignment. Mostly a stability check — +the bottom-`k` subspace is mathematically unique, so disagreement across +starts indicates inverse iteration hasn't fully converged.""" +function _residual_subspace_alignment( + Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, + Rv::AbstractVector{Float64}; + k::Int = STAGNATION_SUBSPACE_K, + n_iter::Int = 10, + n_starts::Int = 3, +) + nF = norm(Rv) + nF > 0 || return Float64[], NaN, Vector{Vector{Float64}}() + F = KLU.klu(Jv) + n = size(Jv, 1) + best_cosines = Float64[] + best_subspace = -Inf + best_V = Vector{Vector{Float64}}() + for _ in 1:n_starts + V = Vector{Vector{Float64}}() + cosines = Float64[] + for _ in 1:k + v = randn(n) + for vj in V + v .-= dot(v, vj) .* vj + end + nv = norm(v) + nv > 0 || break + v ./= nv + converged_inner = true + for _ in 1:n_iter + w = F' \ v + v = F \ w + for vj in V + v .-= dot(v, vj) .* vj + end + nv = norm(v) + if nv == 0 + converged_inner = false + break + end + v ./= nv + end + converged_inner || break + push!(V, v) + u = Jv * v + nu = norm(u) + nu > 0 || break + u ./= nu + push!(cosines, abs(dot(u, Rv)) / nF) + end + subspace = isempty(cosines) ? -Inf : sqrt(sum(c^2 for c in cosines)) + if subspace > best_subspace + best_subspace = subspace + best_cosines = cosines + best_V = V + end + end + best_subspace == -Inf && return Float64[], NaN, Vector{Vector{Float64}}() + return best_cosines, best_subspace, best_V +end + +"""Format a limit-cycle diagnostic line. Skips the point-wise classifier +(gradient / curvature / fold) since those tests are meaningless when the +iterate isn't actually at a single point. Reports the ρ and κ̂ ranges observed +across the window (which capture the cycle's amplitude) plus condest and the +top mismatches. Always returns `status = ACPowerFlowSolveStatus.LIMIT_CYCLE`.""" +function _limit_cycle_diagnostic( + Rv::AbstractVector{Float64}, + ρ_window::AbstractVector{Float64}, + κ_window::AbstractVector{Float64}, +) + sf(x) = round(x; sigdigits = 4) + ρ_part = if isempty(ρ_window) + "" + else + ρ_lo, ρ_hi = extrema(ρ_window) + ", ρ alternates in [$(sf(ρ_lo)), $(sf(ρ_hi))] (×$(sf(ρ_hi/max(ρ_lo,eps()))))" + end + κ_part = if isempty(κ_window) + "" + else + κ_lo, κ_hi = extrema(κ_window) + ", κ̂ alternates in [$(sf(κ_lo)), $(sf(κ_hi))] (×$(sf(κ_hi/max(κ_lo,eps()))))" + end + k_top = min(3, length(Rv)) + top_ix = partialsortperm(Rv, 1:k_top; by = abs, rev = true) + parts = String[] + for ix in top_ix + pow_type = ix % 2 == 1 ? "P" : "Q" + bus_ix = div(ix + 1, 2) + push!(parts, "bus $bus_ix ($pow_type) = $(sf(Rv[ix]))") + end + status = ACPowerFlowSolveStatus.LIMIT_CYCLE + msg = + ρ_part * κ_part * + "; top $k_top mismatches: " * join(parts, ", ") * + " [status: $status]" + return msg, status +end + +"""Classify why a Newton-type AC power flow solver got stuck, and produce a +log-friendly summary. Self-contained: takes only `(data, time_step, Rv, Jv)` +plus optional `condest`; no dependency on solver-internal structures. + +Returns `(message::String, status::ACPowerFlowSolveStatus)`. The status is one +of `NON_ROOT_LOCAL_MIN`, `NON_ROOT_SADDLE`, `NON_ROOT_STATIONARY`, +`SINGULAR_SUBSPACE`, `FOLD`, or `STAGNATED_OTHER`, picked in priority order +(most actionable first): + +1. **NON_ROOT_(LOCAL_MIN|SADDLE|STATIONARY)** — `‖JᵀF‖/(‖J‖_F·‖F‖) < + $(STAGNATION_STATIONARY_THRESHOLD)`: Gauss-Newton stationary point. Sign of + the directional curvature `v_minᵀ ∇²(½‖F‖²) v_min` further distinguishes: + `+` → genuine local min (infeasibility), `−` → saddle (a different x0 + might find a root), undeterminable → fall back to `NON_ROOT_STATIONARY`. +2. **SINGULAR_SUBSPACE** — subspace alignment ≥ + $(STAGNATION_SUBSPACE_THRESHOLD): `F` lies in the bottom-`k` singular + subspace of `J`; trapped at or past a fold of the state manifold. +3. **FOLD** — `‖½H[Δx,Δx]‖_∞/‖F‖_∞ ≥ $(STAGNATION_FOLD_THRESHOLD)`: a full + Newton step would not reduce `‖F‖` because nonlinearity dominates. +4. **STAGNATED_OTHER** — none of the above fire; mechanism unclear. + +The Hessian-vector-vector product is hardcoded to the polar formulation via +[`acpf_hvvp`](@ref); other formulations would need their own hvvp routine to +hook into the FOLD and curvature checks.""" +function _stagnation_diagnostic( + data::ACPowerFlowData, + time_step::Int, + Rv::AbstractVector{Float64}, + Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}; + condest::Union{Float64, Nothing} = nothing, + k::Int = STAGNATION_SUBSPACE_K, +) + sf(x) = round(x; sigdigits = 4) + cond_part = isnothing(condest) || !isfinite(condest) ? "" : + ", κ̂(J) = $(sf(condest))" + cosines, subspace, V = _residual_subspace_alignment(Jv, Rv; k = k) + align_part = if isempty(cosines) + "" + else + cos_list = join((sf(c) for c in cosines), ", ") + ", |⟨F, u_i⟩|/‖F‖ for i=1..$(length(cosines)) = ($cos_list), " * + "subspace alignment = $(sf(subspace))" + end + # Stationary-point check for ½‖F‖²: gradient is JᵀF; at a Gauss-Newton local + # minimum, JᵀF ≈ 0. Normalize by ‖J‖_F·‖F‖ (Cauchy-Schwarz upper bound), so + # the result is in [0, 1] and ≪ 1 means we're at a non-root stationary point. + JtF = Jv' * Rv + nF = norm(Rv) + nJ_F = norm(Jv.nzval) # Frobenius norm of J + grad_ratio = (nF > 0 && nJ_F > 0) ? norm(JtF) / (nJ_F * nF) : NaN + grad_part = isfinite(grad_ratio) ? + ", ‖JᵀF‖/(‖J‖_F·‖F‖) = $(sf(grad_ratio))" : "" + # Fold check: predicted residual after a full Newton step is + # F + J·Δx + ½ H[Δx,Δx,·] = ½ H[Δx,Δx,·]. Ratio ≈ 0 → Newton converges + # cleanly; ≳ 1 → curvature kills the linear step. + quad_ratio = NaN + quad_part = try + F_klu = KLU.klu(Jv) + Δx = -(F_klu \ Vector(Rv)) + h = acpf_hvvp(data, time_step, Δx, Δx) + nh = norm(h, Inf) + nFinf = norm(Rv, Inf) + quad_ratio = (nFinf > 0) ? 0.5 * nh / nFinf : NaN + isfinite(quad_ratio) ? + ", ‖½H[Δx,Δx]‖_∞/‖F‖_∞ = $(sf(quad_ratio))" : "" + catch _ + "" + end + # Second-order test along v_min (right singular vector of smallest σ). + # ∇²(½‖F‖²) = JᵀJ + Σ F_k Hᵏ; directional curvature at v_min is + # ‖J·v_min‖² + Fᵀ·H[v_min, v_min, ·]. Near singularity the first term ≈ + # σ_min² is tiny, so sign is decided by the F·H part. + → robust local min, + # − → saddle. + curv_ratio = NaN + curv_part = if !isempty(V) + try + v_min = V[1] + Jv_min = Jv * v_min + h_vv = acpf_hvvp(data, time_step, v_min, v_min) + curvature = dot(Jv_min, Jv_min) + dot(Rv, h_vv) + scale = max(dot(Jv_min, Jv_min), abs(dot(Rv, h_vv)), eps()) + curv_ratio = curvature / scale + if isfinite(curv_ratio) + ", v_minᵀ∇²(½‖F‖²)v_min (normalized) = $(sf(curv_ratio))" + else + "" + end + catch _ + "" + end + else + "" + end + # Classify + status = if isfinite(grad_ratio) && grad_ratio < STAGNATION_STATIONARY_THRESHOLD + if isfinite(curv_ratio) && curv_ratio < 0 + ACPowerFlowSolveStatus.NON_ROOT_SADDLE + elseif isfinite(curv_ratio) && curv_ratio > 0 + ACPowerFlowSolveStatus.NON_ROOT_LOCAL_MIN + else + ACPowerFlowSolveStatus.NON_ROOT_STATIONARY + end + elseif isfinite(subspace) && subspace >= STAGNATION_SUBSPACE_THRESHOLD + ACPowerFlowSolveStatus.SINGULAR_SUBSPACE + elseif isfinite(quad_ratio) && quad_ratio >= STAGNATION_FOLD_THRESHOLD + ACPowerFlowSolveStatus.FOLD + else + ACPowerFlowSolveStatus.STAGNATED_OTHER + end + k_top = min(3, length(Rv)) + top_ix = partialsortperm(Rv, 1:k_top; by = abs, rev = true) + parts = String[] + for ix in top_ix + pow_type = ix % 2 == 1 ? "P" : "Q" + bus_ix = div(ix + 1, 2) + push!(parts, "bus $bus_ix ($pow_type) = $(sf(Rv[ix]))") + end + msg = + cond_part * align_part * grad_part * quad_part * curv_part * + "; top $k_top mismatches: " * join(parts, ", ") * + " [status: $status]" + return msg, status +end From cd1be0cfa1d966dc6bc4eae50b16fe87be6155ca Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Fri, 22 May 2026 14:15:54 -0600 Subject: [PATCH 10/14] small test fixes --- src/power_flow_method.jl | 16 ++++++++++++---- test/test_dc_power_flow.jl | 2 +- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 5e0840f8b..22b47dc32 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1046,10 +1046,18 @@ function _finalize_power_flow( else "failed to converge" end - @warn( - "The $solver_name solver $msg after $i iterations [status: $status]. " * - "Reporting non-convergence." - ) + # Stagnation-classified statuses get @warn (we know what happened, just + # couldn't reach tol); plain FAILED / DIVERGED keep @error to match + # historical behavior — divergence already logged its own @error inside LM, + # but the duplicate at finalize is harmless. + full_msg = "The $solver_name solver $msg after $i iterations " * + "[status: $status]. Reporting non-convergence." + if status === ACPowerFlowSolveStatus.FAILED || + status === ACPowerFlowSolveStatus.DIVERGED + @error(full_msg) + else + @warn(full_msg) + end return false end diff --git a/test/test_dc_power_flow.jl b/test/test_dc_power_flow.jl index a313a27b7..e65261db6 100644 --- a/test/test_dc_power_flow.jl +++ b/test/test_dc_power_flow.jl @@ -124,7 +124,7 @@ end power_flow_with_units(sys, T, PSY.UnitSystem.NATURAL_UNITS) line_name2, flow_system = power_flow_with_units(sys, T, PSY.UnitSystem.SYSTEM_BASE) @test line_name == line_name2 - @test flow_natural == flow_system + @test isapprox(flow_natural, flow_system) end end From 7ce686d762ad820d5ad07a3946ae11b717703a7b Mon Sep 17 00:00:00 2001 From: Luke Kiernan <86331877+luke-kiernan@users.noreply.github.com> Date: Fri, 22 May 2026 16:23:50 -0600 Subject: [PATCH 11/14] format Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/power_flow_method.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 22b47dc32..b2ebc681b 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1050,8 +1050,9 @@ function _finalize_power_flow( # couldn't reach tol); plain FAILED / DIVERGED keep @error to match # historical behavior — divergence already logged its own @error inside LM, # but the duplicate at finalize is harmless. - full_msg = "The $solver_name solver $msg after $i iterations " * - "[status: $status]. Reporting non-convergence." + full_msg = + "The $solver_name solver $msg after $i iterations " * + "[status: $status]. Reporting non-convergence." if status === ACPowerFlowSolveStatus.FAILED || status === ACPowerFlowSolveStatus.DIVERGED @error(full_msg) From 9601bf53ec6d3348871ba5e31c3925687362c335 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 26 May 2026 16:29:28 -0600 Subject: [PATCH 12/14] Make solver diagnostics formulation-aware MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The per-iteration, stagnation, and limit-cycle diagnostics decoded residual indices assuming the polar [P,Q]-interleaved 2-per-bus layout (ix%2 for P/Q, div(ix+1,2) for the bus). That mislabels the rectangular-CI and mixed-CPB residuals, whose entries are current/voltage mismatches on variable-width per-bus blocks plus an LCC tail, and reports the wrong bus. - Add _describe_residual_entry, dispatched on residual type, that maps a global index to the correct bus (via bus_state_offset for the variable-block formulations) and labels the quantity (P/Q, ΔI_re/ΔI_im, |V|²-V_set², or an LCC equation row). - Guard the acpf_hvvp-based FOLD and directional-curvature checks in _stagnation_diagnostic to the polar formulation only; acpf_hvvp is polar-only, so on a non-polar Δx it either throws or returns a silently-wrong number. Non-polar runs now note the checks were skipped. - Thread residual/data/time_step through _log_diagnostics, _stagnation_diagnostic, and _limit_cycle_diagnostic and their call sites in the NR/TR/adaptive-TR loops and Levenberg-Marquardt. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/levenberg-marquardt.jl | 8 +- src/power_flow_method.jl | 36 ++++---- src/stagnation_diagnostics.jl | 152 ++++++++++++++++++++++++++++------ 3 files changed, 151 insertions(+), 45 deletions(-) diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index 4d78388d0..bacce1b10 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -201,20 +201,22 @@ function _run_power_flow_method( if monitor_ρ ρ_now, _, condest_now = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("LM iter $(i + iter_offset)", ρ_now, residual.Rv, condest_now) + _log_diagnostics( + "LM iter $(i + iter_offset)", ρ_now, residual, J.data, time_step, + condest_now) end if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING kind = _check_stagnation!( F_window, ρ_window, κ_window, F_inf, ρ_now, condest_now) if kind === :fixed_point msg, status = _stagnation_diagnostic( - J.data, time_step, residual.Rv, J.Jv; condest = condest_now) + residual, J.data, time_step, J.Jv; condest = condest_now) @info "LM stagnated at fixed point: ‖F‖_∞ = " * "$(siground(F_inf)), stable across " * "$(STAGNATION_WINDOW) iterations" * msg elseif kind === :limit_cycle msg, status = _limit_cycle_diagnostic( - residual.Rv, ρ_window, κ_window) + residual, J.data, time_step, ρ_window, κ_window) @info "LM in limit cycle: ‖F‖_∞ = " * "$(siground(F_inf)) across " * "$(STAGNATION_WINDOW) iterations" * msg diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index b2ebc681b..f61891a6d 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -4,14 +4,16 @@ rounded to 4 significant figures.""" function _log_diagnostics( label::String, ρ::Float64, - Rv::AbstractVector{Float64}, + residual, + data::ACPowerFlowData, + time_step::Int, condest::Float64, ) + Rv = residual.Rv abs_max, ix = findmax(abs, Rv) - pow_type = ix % 2 == 1 ? "P" : "Q" - bus_ix = div(ix + 1, 2) sf(x) = round(x; sigdigits = 4) - @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at bus $bus_ix ($pow_type), " * + @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at " * + "$(_describe_residual_entry(residual, data, time_step, ix)), " * "κ̂(J) = $(sf(condest))" return end @@ -653,7 +655,7 @@ function _run_power_flow_method(time_step::Int, if monitor_ρ ρ_now, _, condest_now = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("NR iter $i", ρ_now, residual.Rv, condest_now) + _log_diagnostics("NR iter $i", ρ_now, residual, J.data, time_step, condest_now) end F_inf_now = norm(residual.Rv, Inf) if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING @@ -661,13 +663,13 @@ function _run_power_flow_method(time_step::Int, F_window, ρ_window, κ_window, F_inf_now, ρ_now, condest_now) if kind === :fixed_point msg, status = _stagnation_diagnostic( - J.data, time_step, residual.Rv, J.Jv; condest = condest_now) + residual, J.data, time_step, J.Jv; condest = condest_now) @info "NR stagnated at fixed point: ‖F‖_∞ = " * "$(siground(F_inf_now)), stable across " * "$(STAGNATION_WINDOW) iterations" * msg elseif kind === :limit_cycle msg, status = _limit_cycle_diagnostic( - residual.Rv, ρ_window, κ_window) + residual, J.data, time_step, ρ_window, κ_window) @info "NR in limit cycle: ‖F‖_∞ = " * "$(siground(F_inf_now)) across " * "$(STAGNATION_WINDOW) iterations" * msg @@ -774,7 +776,8 @@ function _run_power_flow_method(time_step::Int, if accepted ρ_cached, _, condest_cached = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("TR iter $i", ρ_cached, residual.Rv, condest_cached) + _log_diagnostics( + "TR iter $i", ρ_cached, residual, J.data, time_step, condest_cached) else @info "TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" end @@ -785,13 +788,13 @@ function _run_power_flow_method(time_step::Int, F_window, ρ_window, κ_window, F_inf_now, ρ_cached, condest_cached) if kind === :fixed_point msg, status = _stagnation_diagnostic( - J.data, time_step, residual.Rv, J.Jv; condest = condest_cached) + residual, J.data, time_step, J.Jv; condest = condest_cached) @info "TR stagnated at fixed point: ‖F‖_∞ = " * "$(siground(F_inf_now)), stable across " * "$(STAGNATION_WINDOW) iterations" * msg elseif kind === :limit_cycle msg, status = _limit_cycle_diagnostic( - residual.Rv, ρ_window, κ_window) + residual, J.data, time_step, ρ_window, κ_window) @info "TR in limit cycle: ‖F‖_∞ = " * "$(siground(F_inf_now)) across " * "$(STAGNATION_WINDOW) iterations" * msg @@ -852,7 +855,8 @@ function _newton_power_flow( converged = norm(residual.Rv, Inf) < tol ρ_x0, _, condest_x0 = _fixed_point_spectral_radius!(data, residual, J, time_step) - _log_diagnostics("x0 (time_step $time_step)", ρ_x0, residual.Rv, condest_x0) + _log_diagnostics( + "x0 (time_step $time_step)", ρ_x0, residual, data, time_step, condest_x0) i = 0 x_final = x0_init @@ -903,7 +907,8 @@ function _newton_power_flow( ) if accepted ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) - _log_diagnostics("Adaptive TR iter $i", ρ, residual.Rv, condest) + _log_diagnostics( + "Adaptive TR iter $i", ρ, residual, J.data, time_step, condest) else @info "Adaptive TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" end @@ -917,7 +922,7 @@ function _newton_power_flow( F_inf_now, ρ_for_window, κ_for_window) if kind === :fixed_point msg, status = _stagnation_diagnostic( - J.data, time_step, residual.Rv, J.Jv; condest = condest) + residual, J.data, time_step, J.Jv; condest = condest) @info "Adaptive TR stagnated at fixed point: ‖F‖_∞ = " * "$(siground(F_inf_now)), stable across " * "$(STAGNATION_WINDOW) iterations; stopping before " * @@ -925,7 +930,7 @@ function _newton_power_flow( break elseif kind === :limit_cycle msg, status = _limit_cycle_diagnostic( - residual.Rv, ρ_window_tr, κ_window_tr) + residual, J.data, time_step, ρ_window_tr, κ_window_tr) @info "Adaptive TR in limit cycle: ‖F‖_∞ = " * "$(siground(F_inf_now)) across " * "$(STAGNATION_WINDOW) iterations; stopping before " * @@ -1133,7 +1138,8 @@ function _newton_power_flow( if get_compute_fixed_point_spectral_radius(data) ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) - _log_diagnostics("x0 (time_step $time_step)", ρ, residual.Rv, condest) + _log_diagnostics( + "x0 (time_step $time_step)", ρ, residual, data, time_step, condest) if ρ >= 1.0 @warn "fixed-point spectral radius ρ = $(round(ρ; sigdigits = 4)) ≥ 1 " * "at x0 (time_step $time_step); Newton-Raphson may not converge " * diff --git a/src/stagnation_diagnostics.jl b/src/stagnation_diagnostics.jl index 9c4486dee..76787b81a 100644 --- a/src/stagnation_diagnostics.jl +++ b/src/stagnation_diagnostics.jl @@ -13,6 +13,84 @@ const STAGNATION_SUBSPACE_THRESHOLD = 0.9 # subspace alignment above → tr # Default number of bottom singular vectors probed by the subspace alignment. const STAGNATION_SUBSPACE_K = 3 +# --------------------------------------------------------------------------- +# Formulation-aware residual-entry labeling for diagnostics. +# +# Decoding a global residual index `ix` into a readable "bus N (quantity)" / +# "LCC #k row" label must dispatch on the residual type, because the layout +# differs: +# * polar (ACPowerFlowResidual): [P₁,Q₁,…,Pₙ,Qₙ | LCC tail], fixed 2/bus. +# * rect CI (ACRectangularCIResidual): variable per-bus blocks (PQ/REF = 2 → +# ΔI real/imag; PV = 3 → ΔI real/imag + |V|²−V_set²), indexed by +# bus_state_offset. +# * mixed (ACMixedCPBResidual): 2/bus, but the two rows depend on bus type. +# Both non-polar layouts share rect's `[bus state | 4·n_LCC tail]` split. +# --------------------------------------------------------------------------- + +const _LCC_RESIDUAL_ROW_NAMES = + ("P-setpoint", "DC-line balance", "rectifier α-limit", "inverter α-limit") + +_diag_bus_number(data, bus_ix::Int) = axes(data.power_network_matrix, 1)[bus_ix] + +function _describe_lcc_residual_entry(data, tail_ix::Int) + i = div(tail_ix - 1, 4) + 1 # 1-based LCC index + row = mod1(tail_ix, 4) + (from_no, to_no) = data.lcc.arcs[i] + return "LCC #$i ($from_no→$to_no) $(_LCC_RESIDUAL_ROW_NAMES[row])" +end + +# Owning bus (matrix index) and 0-based row within its block, for the +# variable-block (rect/mixed) formulations. `bus_state_offset[b]` is the 1-based +# block start; `searchsortedlast` returns the block containing `ix`. +function _locate_variable_block(bus_state_offset::AbstractVector, ix::Int) + b = searchsortedlast(bus_state_offset, ix) + return b, ix - Int(bus_state_offset[b]) +end + +""" + _describe_residual_entry(residual, data, time_step, ix) -> String + +Formulation-aware label for global residual index `ix`: the owning bus number +(or LCC arc) and the physical quantity of that equation. Used by the +stagnation / limit-cycle / per-iteration diagnostics so they don't mislabel +non-polar residuals — whose entries are current/voltage mismatches on +variable-width blocks, not the polar fixed 2-per-bus P/Q.""" +function _describe_residual_entry(::ACPowerFlowResidual, data, time_step::Int, ix::Int) + n_bus_eqs = 2 * size(data.bus_type, 1) + ix > n_bus_eqs && return _describe_lcc_residual_entry(data, ix - n_bus_eqs) + bus_ix = div(ix + 1, 2) + qty = isodd(ix) ? "P" : "Q" + return "bus $(_diag_bus_number(data, bus_ix)) ($qty)" +end + +function _describe_residual_entry( + r::ACRectangularCIResidual, data, time_step::Int, ix::Int, +) + ix > r.total_bus_state && + return _describe_lcc_residual_entry(data, ix - r.total_bus_state) + bus_ix, row = _locate_variable_block(r.bus_state_offset, ix) + # 2-block (PQ/REF): real/imag current mismatch. 3-block (PV) adds |V|²−V_set². + qty = row == 0 ? "ΔI_re" : row == 1 ? "ΔI_im" : "|V|²−V_set²" + return "bus $(_diag_bus_number(data, bus_ix)) ($qty)" +end + +function _describe_residual_entry(r::ACMixedCPBResidual, data, time_step::Int, ix::Int) + ix > r.total_bus_state && + return _describe_lcc_residual_entry(data, ix - r.total_bus_state) + bus_ix, row = _locate_variable_block(r.bus_state_offset, ix) + bt = data.bus_type[bus_ix, time_step] + # MCPB layout: PQ = divided-current balance, imag-first; PV = power balance + # then |V|²−V_set²; REF = rect's real/imag current rows. + qty = if bt == PSY.ACBusTypes.PV + row == 0 ? "ΔP" : "|V|²−V_set²" + elseif bt == PSY.ACBusTypes.PQ + row == 0 ? "ΔI_im" : "ΔI_re" # imag-first + else # REF + row == 0 ? "ΔI_re" : "ΔI_im" + end + return "bus $(_diag_bus_number(data, bus_ix)) ($qty)" +end + # Stability-gate parameters for `_check_stagnation!`. Two gates run in parallel: # the **strict** gate (`ρ` and `κ̂` within ~10%) detects a true fixed point of # the iteration map; the **loose** gate (~2×) detects period-≥2 limit cycles @@ -161,10 +239,13 @@ iterate isn't actually at a single point. Reports the ρ and κ̂ ranges observe across the window (which capture the cycle's amplitude) plus condest and the top mismatches. Always returns `status = ACPowerFlowSolveStatus.LIMIT_CYCLE`.""" function _limit_cycle_diagnostic( - Rv::AbstractVector{Float64}, + residual, + data::ACPowerFlowData, + time_step::Int, ρ_window::AbstractVector{Float64}, κ_window::AbstractVector{Float64}, ) + Rv = residual.Rv sf(x) = round(x; sigdigits = 4) ρ_part = if isempty(ρ_window) "" @@ -180,12 +261,10 @@ function _limit_cycle_diagnostic( end k_top = min(3, length(Rv)) top_ix = partialsortperm(Rv, 1:k_top; by = abs, rev = true) - parts = String[] - for ix in top_ix - pow_type = ix % 2 == 1 ? "P" : "Q" - bus_ix = div(ix + 1, 2) - push!(parts, "bus $bus_ix ($pow_type) = $(sf(Rv[ix]))") - end + parts = [ + "$(_describe_residual_entry(residual, data, time_step, ix)) = $(sf(Rv[ix]))" + for ix in top_ix + ] status = ACPowerFlowSolveStatus.LIMIT_CYCLE msg = ρ_part * κ_part * @@ -216,16 +295,24 @@ of `NON_ROOT_LOCAL_MIN`, `NON_ROOT_SADDLE`, `NON_ROOT_STATIONARY`, 4. **STAGNATED_OTHER** — none of the above fire; mechanism unclear. The Hessian-vector-vector product is hardcoded to the polar formulation via -[`acpf_hvvp`](@ref); other formulations would need their own hvvp routine to -hook into the FOLD and curvature checks.""" +[`acpf_hvvp`](@ref), so for non-polar formulations (rectangular CI, mixed CPB) +the FOLD and directional-curvature checks are skipped entirely — they would +need a formulation-specific hvvp routine. The gradient/subspace tests, which +use only `J` and `F`, still run for every formulation.""" function _stagnation_diagnostic( + residual, data::ACPowerFlowData, time_step::Int, - Rv::AbstractVector{Float64}, Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}; condest::Union{Float64, Nothing} = nothing, k::Int = STAGNATION_SUBSPACE_K, ) + Rv = residual.Rv + # acpf_hvvp is polar-only; running it on a non-polar Δx (different basis, + # possibly different length) would either throw or — if the dimensions + # happen to coincide — return a silently wrong number. Gate the FOLD and + # directional-curvature checks on the polar formulation. + is_polar = residual isa ACPowerFlowResidual sf(x) = round(x; sigdigits = 4) cond_part = isnothing(condest) || !isfinite(condest) ? "" : ", κ̂(J) = $(sf(condest))" @@ -250,17 +337,21 @@ function _stagnation_diagnostic( # F + J·Δx + ½ H[Δx,Δx,·] = ½ H[Δx,Δx,·]. Ratio ≈ 0 → Newton converges # cleanly; ≳ 1 → curvature kills the linear step. quad_ratio = NaN - quad_part = try - F_klu = KLU.klu(Jv) - Δx = -(F_klu \ Vector(Rv)) - h = acpf_hvvp(data, time_step, Δx, Δx) - nh = norm(h, Inf) - nFinf = norm(Rv, Inf) - quad_ratio = (nFinf > 0) ? 0.5 * nh / nFinf : NaN - isfinite(quad_ratio) ? - ", ‖½H[Δx,Δx]‖_∞/‖F‖_∞ = $(sf(quad_ratio))" : "" - catch _ + quad_part = if !is_polar "" + else + try + F_klu = KLU.klu(Jv) + Δx = -(F_klu \ Vector(Rv)) + h = acpf_hvvp(data, time_step, Δx, Δx) + nh = norm(h, Inf) + nFinf = norm(Rv, Inf) + quad_ratio = (nFinf > 0) ? 0.5 * nh / nFinf : NaN + isfinite(quad_ratio) ? + ", ‖½H[Δx,Δx]‖_∞/‖F‖_∞ = $(sf(quad_ratio))" : "" + catch _ + "" + end end # Second-order test along v_min (right singular vector of smallest σ). # ∇²(½‖F‖²) = JᵀJ + Σ F_k Hᵏ; directional curvature at v_min is @@ -268,7 +359,9 @@ function _stagnation_diagnostic( # σ_min² is tiny, so sign is decided by the F·H part. + → robust local min, # − → saddle. curv_ratio = NaN - curv_part = if !isempty(V) + curv_part = if !is_polar + "" + elseif !isempty(V) try v_min = V[1] Jv_min = Jv * v_min @@ -305,14 +398,19 @@ function _stagnation_diagnostic( end k_top = min(3, length(Rv)) top_ix = partialsortperm(Rv, 1:k_top; by = abs, rev = true) - parts = String[] - for ix in top_ix - pow_type = ix % 2 == 1 ? "P" : "Q" - bus_ix = div(ix + 1, 2) - push!(parts, "bus $bus_ix ($pow_type) = $(sf(Rv[ix]))") - end + parts = [ + "$(_describe_residual_entry(residual, data, time_step, ix)) = $(sf(Rv[ix]))" + for ix in top_ix + ] + # Non-polar formulations skip the curvature/FOLD tests (acpf_hvvp is + # polar-only), so make that explicit rather than letting their absence read + # as "curvature was inconclusive". + formulation_part = + is_polar ? "" : + ", curvature/FOLD checks skipped (non-polar formulation)" msg = cond_part * align_part * grad_part * quad_part * curv_part * + formulation_part * "; top $k_top mismatches: " * join(parts, ", ") * " [status: $status]" return msg, status From b7eab224808b1a78cc80696099bb6dd0bffb8857 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 27 May 2026 08:21:07 -0600 Subject: [PATCH 13/14] Add smallest-magnitude Jacobian eigenvalue routine compute_min_jacobian_eigenvalue finds the eigenvalue of J closest to the origin via inverse iteration, reusing the KLU factor as matvec v -> J\v. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/PowerFlows.jl | 1 + src/jacobian_eigenvalues.jl | 77 +++++++++++++++++++++++++++++++ test/test_jacobian_eigenvalues.jl | 56 ++++++++++++++++++++++ 3 files changed, 134 insertions(+) create mode 100644 src/jacobian_eigenvalues.jl create mode 100644 test/test_jacobian_eigenvalues.jl diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index ca3cacf76..cdb1b57ff 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -76,6 +76,7 @@ include("mixed_cpb_setup.jl") include("mixed_cpb_power_flow_residual.jl") include("mixed_cpb_power_flow_jacobian.jl") include("fixed_point_spectral_radius.jl") +include("jacobian_eigenvalues.jl") include("stagnation_diagnostics.jl") include("solve_ac_power_flow.jl") include("power_flow_setup.jl") diff --git a/src/jacobian_eigenvalues.jl b/src/jacobian_eigenvalues.jl new file mode 100644 index 000000000..072817000 --- /dev/null +++ b/src/jacobian_eigenvalues.jl @@ -0,0 +1,77 @@ +""" + compute_min_jacobian_eigenvalue(data, time_step; x0, tol, maxiter, krylovdim) + -> (λ, info, condest::Float64) + +Estimate the eigenvalue of the AC power flow Jacobian ``J(x)`` of smallest +magnitude (the one closest to the origin) at the state `x0` (default: the result +of `calculate_x0(data, time_step)`). The eigenvalue closest to zero measures how +near ``J`` is to singular: at a saddle-node bifurcation (voltage collapse) a real +eigenvalue of ``J`` crosses zero, so a small ``|λ_min|`` flags an ill-conditioned +or near-singular operating point where Newton-type methods struggle. + +# Method + +``J`` is large, sparse, and non-symmetric, so its eigenvalues are in general +complex and the smallest-magnitude end of the spectrum is the hard end for a +direct Arnoldi iteration on ``J``. Instead we use **inverse iteration**, reusing +the sparse `KLU` factorization of ``J`` that the solver already builds: the +matvec ``v \\mapsto J^{-1} v`` is applied as a back-solve against that factor — +the inverse is *never* formed explicitly (it would be dense). KrylovKit's +matrix-free Arnoldi solver then extracts the eigenvalue ``μ`` of ``J^{-1}`` of +largest magnitude, and we report ``λ_min = 1 / μ``, the eigenvalue of ``J`` +closest to the origin. + +Returns `(λ_min, info, condest)`, where `λ_min` may be complex, `info` is the +KrylovKit convergence info, and `condest` is a Hager 1-norm estimate of the +condition number of ``J`` computed from the same `KLU` factor. + +# Notes + +- Unlike [`compute_fixed_point_spectral_radius`](@ref), this needs only the + Jacobian itself (not the residual Hessian), so it supports every AC + formulation, including LCC HVDC systems. +""" +function compute_min_jacobian_eigenvalue( + data::ACPowerFlowData, + time_step::Int; + x0::Union{Vector{Float64}, Nothing} = nothing, + tol::Float64 = 1e-6, + maxiter::Int = 200, + krylovdim::Int = 30, +) + residual = ACPowerFlowResidual(data, time_step) + jac = ACPowerFlowJacobian(residual, time_step) + x = isnothing(x0) ? calculate_x0(data, time_step) : copy(x0) + residual(x, time_step) + jac(time_step) + return _min_jacobian_eigenvalue!( + jac; + tol = tol, maxiter = maxiter, krylovdim = krylovdim, + ) +end + +"""In-place smallest-magnitude Jacobian eigenvalue computation that reuses an +already-evaluated `jac` (its `jac.Jv` must hold the Jacobian at the current +state). Factorizes `jac.Jv` once with `KLU` and runs inverse iteration; see +[`compute_min_jacobian_eigenvalue`](@ref). Returns `(λ_min, info, condest)`.""" +function _min_jacobian_eigenvalue!( + jac::ACPowerFlowJacobian; + tol::Float64 = 1e-6, + maxiter::Int = 200, + krylovdim::Int = 30, +) + n = size(jac.Jv, 1) + F = KLU.klu(jac.Jv) + # matvec v ↦ J⁻¹ v: a back-solve against the existing factor, never forming J⁻¹. + matvec(v::AbstractVector) = F \ v + # Deterministic init for reproducibility across runs / CI logs. + v_init = ones(Float64, n) ./ sqrt(n) + # Largest-magnitude eigenvalue μ of J⁻¹ ⇒ smallest-magnitude eigenvalue 1/μ of J. + vals, _, info = KrylovKit.eigsolve( + matvec, v_init, 1, :LM; + tol = tol, maxiter = maxiter, krylovdim = krylovdim, + ) + λ_min = isempty(vals) ? complex(NaN, NaN) : inv(vals[1]) + condest = KLU.condest(F) + return λ_min, info, condest +end diff --git a/test/test_jacobian_eigenvalues.jl b/test/test_jacobian_eigenvalues.jl new file mode 100644 index 000000000..88f776845 --- /dev/null +++ b/test/test_jacobian_eigenvalues.jl @@ -0,0 +1,56 @@ +@testset "test compute_min_jacobian_eigenvalue matches dense eigvals" begin + # Ground-truth check: the inverse-iteration estimate of the smallest-magnitude + # eigenvalue of J should match the dense `eigvals` of the same Jacobian to + # solver tolerance. J is non-symmetric so the eigenvalue is generally complex; + # we compare the full complex value, not just its magnitude. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PowerFlowData(pf, sys) + time_step = 1 + + # Build the Jacobian at the flat start, exactly as the routine does internally. + residual = PF.ACPowerFlowResidual(data, time_step) + jac = PF.ACPowerFlowJacobian(residual, time_step) + x0 = PF.calculate_x0(data, time_step) + residual(x0, time_step) + jac(time_step) + + ev = eigvals(Matrix(jac.Jv)) + λ_true = ev[argmin(abs.(ev))] + + λ, info, condest = PF.compute_min_jacobian_eigenvalue(data, time_step) + @test info.converged >= 1 + # Match the complex eigenvalue closest to the origin to tight relative error. + @test abs(λ - λ_true) / abs(λ_true) < 1e-8 + # 1/λ is the largest-magnitude eigenvalue of J⁻¹; sanity-check the inverse map. + @test abs(abs(λ) - minimum(abs.(ev))) < 1e-8 + + # condest is a 1-norm condition estimate from the same KLU factor; healthy + # 14-bus system should be well-conditioned. + @test isfinite(condest) + @test 0 < condest < 1e6 +end + +@testset "test compute_min_jacobian_eigenvalue accepts a custom x0" begin + # The eigenvalue is state-dependent; evaluating at a perturbed point should + # still agree with the dense ground truth built at that same point. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PowerFlowData(pf, sys) + time_step = 1 + + residual = PF.ACPowerFlowResidual(data, time_step) + jac = PF.ACPowerFlowJacobian(residual, time_step) + x0 = PF.calculate_x0(data, time_step) + Random.seed!(11) + x_pert = x0 .+ 0.05 .* randn(length(x0)) + + residual(x_pert, time_step) + jac(time_step) + ev = eigvals(Matrix(jac.Jv)) + λ_true = ev[argmin(abs.(ev))] + + λ, info, _ = PF.compute_min_jacobian_eigenvalue(data, time_step; x0 = x_pert) + @test info.converged >= 1 + @test abs(λ - λ_true) / abs(λ_true) < 1e-8 +end From 29ac70ca8d69dd30f3759f1f9a813f7199daef97 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 27 May 2026 13:48:55 -0600 Subject: [PATCH 14/14] Wire min-eigenvalue diagnostic to its own flag MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit compute_min_jacobian_eigenvalue is now a separate solver flag on all three AC formulations, decoupled from compute_fixed_point_spectral_radius. The monitor logs whichever of {ρ, κ̂, λ_min} were computed, so λ_min works for rect-CI, mixed-CPB, and LCC systems (where the polar-only spectral radius throws). Co-Authored-By: Claude Opus 4.7 (1M context) --- src/PowerFlowData.jl | 2 + src/jacobian_eigenvalues.jl | 21 +++--- src/levenberg-marquardt.jl | 12 +++- src/power_flow_method.jl | 108 ++++++++++++++++++++++++------ src/power_flow_types.jl | 14 ++++ test/test_jacobian_eigenvalues.jl | 77 +++++++++++++++++++++ 6 files changed, 199 insertions(+), 35 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 2100af9a0..3e84675ef 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -226,6 +226,8 @@ get_calculate_voltage_stability_factors(pfd::PowerFlowData) = get_calculate_voltage_stability_factors(pfd.pf) get_compute_fixed_point_spectral_radius(pfd::PowerFlowData) = get_compute_fixed_point_spectral_radius(pfd.pf) +get_compute_min_jacobian_eigenvalue(pfd::PowerFlowData) = + get_compute_min_jacobian_eigenvalue(pfd.pf) get_network_reductions(pfd::PowerFlowData) = get_network_reductions(pfd.pf) get_time_steps(pfd::PowerFlowData) = get_time_steps(pfd.pf) get_time_step_names(pfd::PowerFlowData) = get_time_step_names(pfd.pf) diff --git a/src/jacobian_eigenvalues.jl b/src/jacobian_eigenvalues.jl index 072817000..a46c26b88 100644 --- a/src/jacobian_eigenvalues.jl +++ b/src/jacobian_eigenvalues.jl @@ -11,15 +11,9 @@ or near-singular operating point where Newton-type methods struggle. # Method -``J`` is large, sparse, and non-symmetric, so its eigenvalues are in general -complex and the smallest-magnitude end of the spectrum is the hard end for a -direct Arnoldi iteration on ``J``. Instead we use **inverse iteration**, reusing -the sparse `KLU` factorization of ``J`` that the solver already builds: the +We employ **inverse iteration**, reusing a sparse `KLU` factorization of ``J``. matvec ``v \\mapsto J^{-1} v`` is applied as a back-solve against that factor — -the inverse is *never* formed explicitly (it would be dense). KrylovKit's -matrix-free Arnoldi solver then extracts the eigenvalue ``μ`` of ``J^{-1}`` of -largest magnitude, and we report ``λ_min = 1 / μ``, the eigenvalue of ``J`` -closest to the origin. +the inverse is *never* formed explicitly, as it would be dense. Returns `(λ_min, info, condest)`, where `λ_min` may be complex, `info` is the KrylovKit convergence info, and `condest` is a Hager 1-norm estimate of the @@ -39,6 +33,7 @@ function compute_min_jacobian_eigenvalue( maxiter::Int = 200, krylovdim::Int = 30, ) + # PERF: have caller pass in the already constructed J and residual. residual = ACPowerFlowResidual(data, time_step) jac = ACPowerFlowJacobian(residual, time_step) x = isnothing(x0) ? calculate_x0(data, time_step) : copy(x0) @@ -53,18 +48,20 @@ end """In-place smallest-magnitude Jacobian eigenvalue computation that reuses an already-evaluated `jac` (its `jac.Jv` must hold the Jacobian at the current state). Factorizes `jac.Jv` once with `KLU` and runs inverse iteration; see -[`compute_min_jacobian_eigenvalue`](@ref). Returns `(λ_min, info, condest)`.""" +[`compute_min_jacobian_eigenvalue`](@ref). Returns `(λ_min, info, condest)`. + +Accepts any AC formulation's Jacobian (polar, rectangular-CI, mixed-CPB) — it +only touches `jac.Jv`, which every formulation provides.""" function _min_jacobian_eigenvalue!( - jac::ACPowerFlowJacobian; + jac::Union{ACPowerFlowJacobian, ACRectangularCIJacobian, ACMixedCPBJacobian}; tol::Float64 = 1e-6, maxiter::Int = 200, krylovdim::Int = 30, ) n = size(jac.Jv, 1) + # PERF: use existing factorization. F = KLU.klu(jac.Jv) - # matvec v ↦ J⁻¹ v: a back-solve against the existing factor, never forming J⁻¹. matvec(v::AbstractVector) = F \ v - # Deterministic init for reproducibility across runs / CI logs. v_init = ones(Float64, n) ./ sqrt(n) # Largest-magnitude eigenvalue μ of J⁻¹ ⇒ smallest-magnitude eigenvalue 1/μ of J. vals, _, info = KrylovKit.eigsolve( diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index bacce1b10..be87aba8b 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -177,6 +177,7 @@ function _run_power_flow_method( tol::Float64 = DEFAULT_NR_TOL, λ_0::Float64 = DEFAULT_λ_0, monitor_ρ::Bool = get_compute_fixed_point_spectral_radius(J.data), + monitor_λ::Bool = get_compute_min_jacobian_eigenvalue(J.data), bail_on_stagnation::Bool = false, detect_stagnation::Bool = true, iter_offset::Int = 0, @@ -198,12 +199,19 @@ function _run_power_flow_method( F_inf = norm(residual.Rv, Inf) ρ_now::Union{Float64, Nothing} = nothing condest_now::Union{Float64, Nothing} = nothing + λ_min_now::Union{Number, Nothing} = nothing if monitor_ρ ρ_now, _, condest_now = _fixed_point_spectral_radius!(J.data, residual, J, time_step) + end + if monitor_λ + λ_min_now, _, condest_λ = _min_jacobian_eigenvalue!(J) + monitor_ρ || (condest_now = condest_λ) + end + if monitor_ρ || monitor_λ _log_diagnostics( - "LM iter $(i + iter_offset)", ρ_now, residual, J.data, time_step, - condest_now) + "LM iter $(i + iter_offset)", residual, J.data, time_step; + ρ = ρ_now, condest = condest_now, λ_min = λ_min_now) end if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING kind = _check_stagnation!( diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index f61891a6d..d55632cf7 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -1,20 +1,41 @@ -"""Format a per-iteration diagnostic line: ρ, ‖F‖_∞ and the bus/equation where -that infinity-norm is attained, and a κ̂(J) condition estimate. All numerics -rounded to 4 significant figures.""" +"""Format a (possibly complex) eigenvalue to 4 significant figures as `a` or +`a ± b im`.""" +function _fmt_eig(z::Number, sf) + iz = imag(z) + return if iz == 0 + "$(sf(real(z)))" + else + "$(sf(real(z))) $(iz < 0 ? "-" : "+") $(sf(abs(iz)))im" + end +end + +"""Format a per-iteration diagnostic line: always ‖F‖_∞ and the bus/equation +where that infinity-norm is attained, plus whichever of the optional diagnostics +were computed — the fixed-point spectral radius `ρ`, the condition estimate +`κ̂(J)`, and `λ_min` (the Jacobian eigenvalue closest to the origin, shown with +its magnitude). All numerics rounded to 4 significant figures. The optional +diagnostics are gated by separate solver flags, so any subset may be `nothing`.""" function _log_diagnostics( label::String, - ρ::Float64, residual, data::ACPowerFlowData, - time_step::Int, - condest::Float64, + time_step::Int; + ρ::Union{Float64, Nothing} = nothing, + condest::Union{Float64, Nothing} = nothing, + λ_min::Union{Number, Nothing} = nothing, ) Rv = residual.Rv abs_max, ix = findmax(abs, Rv) sf(x) = round(x; sigdigits = 4) - @info "$label: ρ = $(sf(ρ)), ‖F‖_∞ = $(sf(abs_max)) at " * - "$(_describe_residual_entry(residual, data, time_step, ix)), " * - "κ̂(J) = $(sf(condest))" + parts = String[] + isnothing(ρ) || push!(parts, "ρ = $(sf(ρ))") + push!(parts, + "‖F‖_∞ = $(sf(abs_max)) at " * + "$(_describe_residual_entry(residual, data, time_step, ix))") + isnothing(condest) || push!(parts, "κ̂(J) = $(sf(condest))") + isnothing(λ_min) || + push!(parts, "λ_min(J) = $(_fmt_eig(λ_min, sf)) (|λ_min| = $(sf(abs(λ_min))))") + @info "$label: " * join(parts, ", ") return end @@ -609,6 +630,7 @@ function _run_power_flow_method(time_step::Int, i, converged = 1, false consecutive_reverts = 0 monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) + monitor_λ = get_compute_min_jacobian_eigenvalue(J.data) status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING F_window = Float64[] ρ_window = Float64[] @@ -652,10 +674,20 @@ function _run_power_flow_method(time_step::Int, ) ρ_now::Union{Float64, Nothing} = nothing condest_now::Union{Float64, Nothing} = nothing + λ_min_now::Union{Number, Nothing} = nothing if monitor_ρ ρ_now, _, condest_now = _fixed_point_spectral_radius!(J.data, residual, J, time_step) - _log_diagnostics("NR iter $i", ρ_now, residual, J.data, time_step, condest_now) + end + if monitor_λ + # min-eigenvalue also yields a condest; use it only when the + # spectral-radius path (its own condest source) didn't run. + λ_min_now, _, condest_λ = _min_jacobian_eigenvalue!(J) + monitor_ρ || (condest_now = condest_λ) + end + if monitor_ρ || monitor_λ + _log_diagnostics("NR iter $i", residual, J.data, time_step; + ρ = ρ_now, condest = condest_now, λ_min = λ_min_now) end F_inf_now = norm(residual.Rv, Inf) if detect_stagnation && status === ACPowerFlowSolveStatus.RUNNING @@ -744,6 +776,7 @@ function _run_power_flow_method(time_step::Int, @debug "initially: sum of squares $(siground(residualSize)), L ∞ norm $(siground(linf)), Δ $(siground(delta))" monitor_ρ = get_compute_fixed_point_spectral_radius(J.data) + monitor_λ = get_compute_min_jacobian_eigenvalue(J.data) status::ACPowerFlowSolveStatus = ACPowerFlowSolveStatus.RUNNING F_window = Float64[] ρ_window = Float64[] @@ -753,6 +786,7 @@ function _run_power_flow_method(time_step::Int, # still feed the stagnation windows so a stretch of rejections registers. ρ_cached::Union{Float64, Nothing} = nothing condest_cached::Union{Float64, Nothing} = nothing + λ_min_cached::Union{Number, Nothing} = nothing while i < maxIterations && !converged delta, accepted = _trust_region_step( time_step, @@ -772,12 +806,19 @@ function _run_power_flow_method(time_step::Int, vm_validation_range, i, ) - if monitor_ρ + if monitor_ρ || monitor_λ if accepted - ρ_cached, _, condest_cached = - _fixed_point_spectral_radius!(J.data, residual, J, time_step) + if monitor_ρ + ρ_cached, _, condest_cached = + _fixed_point_spectral_radius!(J.data, residual, J, time_step) + end + if monitor_λ + λ_min_cached, _, condest_λ = _min_jacobian_eigenvalue!(J) + monitor_ρ || (condest_cached = condest_λ) + end _log_diagnostics( - "TR iter $i", ρ_cached, residual, J.data, time_step, condest_cached) + "TR iter $i", residual, J.data, time_step; + ρ = ρ_cached, condest = condest_cached, λ_min = λ_min_cached) else @info "TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" end @@ -855,8 +896,14 @@ function _newton_power_flow( converged = norm(residual.Rv, Inf) < tol ρ_x0, _, condest_x0 = _fixed_point_spectral_radius!(data, residual, J, time_step) + λ_min_x0 = if get_compute_min_jacobian_eigenvalue(data) + first(_min_jacobian_eigenvalue!(J)) + else + nothing + end _log_diagnostics( - "x0 (time_step $time_step)", ρ_x0, residual, data, time_step, condest_x0) + "x0 (time_step $time_step)", residual, data, time_step; + ρ = ρ_x0, condest = condest_x0, λ_min = λ_min_x0) i = 0 x_final = x0_init @@ -907,8 +954,14 @@ function _newton_power_flow( ) if accepted ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) + λ_min_now = if get_compute_min_jacobian_eigenvalue(data) + first(_min_jacobian_eigenvalue!(J)) + else + nothing + end _log_diagnostics( - "Adaptive TR iter $i", ρ, residual, J.data, time_step, condest) + "Adaptive TR iter $i", residual, J.data, time_step; + ρ = ρ, condest = condest, λ_min = λ_min_now) else @info "Adaptive TR iter $i: step rejected (state unchanged; ρ/κ̂ skipped)" end @@ -1136,12 +1189,25 @@ function _newton_power_flow( pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol - if get_compute_fixed_point_spectral_radius(data) - ρ, _, condest = _fixed_point_spectral_radius!(data, residual, J, time_step) + monitor_ρ_x0 = get_compute_fixed_point_spectral_radius(data) + monitor_λ_x0 = get_compute_min_jacobian_eigenvalue(data) + if monitor_ρ_x0 || monitor_λ_x0 + ρ_x0::Union{Float64, Nothing} = nothing + condest_x0::Union{Float64, Nothing} = nothing + λ_min_x0::Union{Number, Nothing} = nothing + if monitor_ρ_x0 + ρ_x0, _, condest_x0 = + _fixed_point_spectral_radius!(data, residual, J, time_step) + end + if monitor_λ_x0 + λ_min_x0, _, condest_λ = _min_jacobian_eigenvalue!(J) + monitor_ρ_x0 || (condest_x0 = condest_λ) + end _log_diagnostics( - "x0 (time_step $time_step)", ρ, residual, data, time_step, condest) - if ρ >= 1.0 - @warn "fixed-point spectral radius ρ = $(round(ρ; sigdigits = 4)) ≥ 1 " * + "x0 (time_step $time_step)", residual, data, time_step; + ρ = ρ_x0, condest = condest_x0, λ_min = λ_min_x0) + if !isnothing(ρ_x0) && ρ_x0 >= 1.0 + @warn "fixed-point spectral radius ρ = $(round(ρ_x0; sigdigits = 4)) ≥ 1 " * "at x0 (time_step $time_step); Newton-Raphson may not converge " * "from this starting point" end diff --git a/src/power_flow_types.jl b/src/power_flow_types.jl index 99cb9695a..fe8d2d0fa 100644 --- a/src/power_flow_types.jl +++ b/src/power_flow_types.jl @@ -168,6 +168,7 @@ struct ACPolarPowerFlow{ACSolver <: ACPowerFlowSolverType} <: AbstractACPowerFlo calculate_loss_factors::Bool calculate_voltage_stability_factors::Bool compute_fixed_point_spectral_radius::Bool + compute_min_jacobian_eigenvalue::Bool generator_slack_participation_factors::Union{ Nothing, Dict{Tuple{DataType, String}, Float64}, @@ -222,6 +223,7 @@ function ACPolarPowerFlow{ACSolver}(; calculate_loss_factors::Bool = false, calculate_voltage_stability_factors::Bool = false, compute_fixed_point_spectral_radius::Bool = false, + compute_min_jacobian_eigenvalue::Bool = false, generator_slack_participation_factors::Union{ Nothing, Dict{Tuple{DataType, String}, Float64}, @@ -252,6 +254,7 @@ function ACPolarPowerFlow{ACSolver}(; calculate_loss_factors, calculate_voltage_stability_factors, compute_fixed_point_spectral_radius, + compute_min_jacobian_eigenvalue, generator_slack_participation_factors, enhanced_flat_start, robust_power_flow, @@ -296,6 +299,11 @@ get_calculate_voltage_stability_factors(pf::ACPolarPowerFlow) = get_compute_fixed_point_spectral_radius(::PowerFlowEvaluationModel) = false get_compute_fixed_point_spectral_radius(pf::ACPolarPowerFlow) = pf.compute_fixed_point_spectral_radius +# Available on every AC formulation (and with LCC): unlike the spectral radius, +# the min-eigenvalue diagnostic only needs the Jacobian itself. +get_compute_min_jacobian_eigenvalue(::PowerFlowEvaluationModel) = false +get_compute_min_jacobian_eigenvalue(pf::AbstractACPowerFlow) = + pf.compute_min_jacobian_eigenvalue """ ACRectangularPowerFlow{ACSolver}(; kwargs...) where {ACSolver <: ACPowerFlowSolverType} @@ -343,6 +351,7 @@ struct ACRectangularPowerFlow{ACSolver <: ACPowerFlowSolverType} <: Vector{Dict{Tuple{DataType, String}, Float64}}, } enhanced_flat_start::Bool + compute_min_jacobian_eigenvalue::Bool skip_redistribution::Bool distribute_slack_proportional_to_headroom::Bool network_reductions::Vector{PNM.NetworkReduction} @@ -361,6 +370,7 @@ function ACRectangularPowerFlow{ACSolver}(; Vector{Dict{Tuple{DataType, String}, Float64}}, } = nothing, enhanced_flat_start::Bool = true, + compute_min_jacobian_eigenvalue::Bool = false, skip_redistribution::Bool = false, distribute_slack_proportional_to_headroom::Bool = false, network_reductions::Vector{PNM.NetworkReduction} = PNM.NetworkReduction[], @@ -394,6 +404,7 @@ function ACRectangularPowerFlow{ACSolver}(; exporter, generator_slack_participation_factors, enhanced_flat_start, + compute_min_jacobian_eigenvalue, skip_redistribution, distribute_slack_proportional_to_headroom, network_reductions, @@ -453,6 +464,7 @@ struct ACMixedPowerFlow{ACSolver <: ACPowerFlowSolverType} <: Vector{Dict{Tuple{DataType, String}, Float64}}, } enhanced_flat_start::Bool + compute_min_jacobian_eigenvalue::Bool skip_redistribution::Bool distribute_slack_proportional_to_headroom::Bool network_reductions::Vector{PNM.NetworkReduction} @@ -471,6 +483,7 @@ function ACMixedPowerFlow{ACSolver}(; Vector{Dict{Tuple{DataType, String}, Float64}}, } = nothing, enhanced_flat_start::Bool = true, + compute_min_jacobian_eigenvalue::Bool = false, skip_redistribution::Bool = false, distribute_slack_proportional_to_headroom::Bool = false, network_reductions::Vector{PNM.NetworkReduction} = PNM.NetworkReduction[], @@ -505,6 +518,7 @@ function ACMixedPowerFlow{ACSolver}(; exporter, generator_slack_participation_factors, enhanced_flat_start, + compute_min_jacobian_eigenvalue, skip_redistribution, distribute_slack_proportional_to_headroom, network_reductions, diff --git a/test/test_jacobian_eigenvalues.jl b/test/test_jacobian_eigenvalues.jl index 88f776845..e0b86963b 100644 --- a/test/test_jacobian_eigenvalues.jl +++ b/test/test_jacobian_eigenvalues.jl @@ -54,3 +54,80 @@ end @test info.converged >= 1 @test abs(λ - λ_true) / abs(λ_true) < 1e-8 end + +# Helper: solve under `pf` and return the per-iteration/x0 diagnostic log lines. +function _diagnostic_lines(pf, sys) + data = PowerFlowData(pf, sys) + tl = Test.TestLogger(; min_level = Logging.Info) + Logging.with_logger(tl) do + solve_power_flow!(data) + end + return [ + r.message for r in tl.logs + if occursin("x0 (time_step", r.message) || occursin(r"iter \d+", r.message) + ] +end + +@testset "compute_min_jacobian_eigenvalue flag is independent of spectral radius" begin + # With only the min-eigenvalue flag set, the per-iteration diagnostic lines + # carry λ_min(J) and κ̂(J) but NOT ρ — the (polar-only) spectral-radius path + # must not run. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; + correct_bustypes = true, compute_min_jacobian_eigenvalue = true) + lines = _diagnostic_lines(pf, sys) + @test length(lines) >= 2 + for line in lines + @test occursin("λ_min(J) = ", line) + @test occursin("|λ_min| = ", line) + @test occursin("κ̂(J) = ", line) + @test !occursin("ρ = ", line) # spectral radius must be absent + end +end + +@testset "compute_min_jacobian_eigenvalue works for the rectangular-CI formulation" begin + # The spectral-radius flag is polar-only (getter returns false for rect/mixed), + # so rect-CI never logged before. The min-eigenvalue flag is available on every + # AC formulation. + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = PF.ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}(; + correct_bustypes = true, compute_min_jacobian_eigenvalue = true) + lines = _diagnostic_lines(pf, sys) + @test length(lines) >= 2 + for line in lines + @test occursin("λ_min(J) = ", line) + @test !occursin("ρ = ", line) + end +end + +@testset "compute_min_jacobian_eigenvalue works on LCC systems (spectral radius does not)" begin + # On an LCC system the spectral-radius monitor errors (its Hessian-vector + # product covers only bus states, mismatching the LCC-augmented factor), + # whereas the min-eigenvalue diagnostic only factorizes the Jacobian and works. + sys = System(joinpath(TEST_DATA_DIR, "case5_2_lcc.raw")) + @test length(collect(get_components(TwoTerminalLCCLine, sys))) == 2 + + # spectral-radius flag: the monitor blows up on the LCC-augmented system. + data_ρ = PowerFlowData( + ACPowerFlow{NewtonRaphsonACPowerFlow}(; + compute_fixed_point_spectral_radius = true), + sys) + @test_throws Exception solve_power_flow!(data_ρ) + + # min-eigenvalue flag: solves cleanly and logs λ_min on every line. + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}(; compute_min_jacobian_eigenvalue = true) + data = PowerFlowData(pf, sys) + tl = Test.TestLogger(; min_level = Logging.Info) + converged = Logging.with_logger(tl) do + solve_power_flow!(data) + end + @test converged + lines = [ + r.message for r in tl.logs + if occursin("x0 (time_step", r.message) || occursin(r"iter \d+", r.message) + ] + @test length(lines) >= 2 + for line in lines + @test occursin("λ_min(J) = ", line) + end +end