add fixed point diagnostics

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"

src/PowerFlowData.jl

@@ -224,6 +224,10 @@ 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_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)

src/PowerFlows.jl

@@ -7,6 +7,8 @@ export NewtonRaphsonACPowerFlow
 export TrustRegionACPowerFlow
 export LevenbergMarquardtACPowerFlow
 export RobustHomotopyPowerFlow
+export AdaptiveACPowerFlow
+export ACPowerFlowSolveStatus
 export ACPolarPowerFlow
 export ACRectangularPowerFlow
 export ACMixedPowerFlow
@@ -37,6 +39,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 +75,9 @@ 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("jacobian_eigenvalues.jl")
+include("stagnation_diagnostics.jl")
 include("solve_ac_power_flow.jl")
 include("power_flow_setup.jl")
 include("power_flow_method.jl")

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 !isnothing(x0)
+        @warn "Using caller-provided x0; skipping improve_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

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

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,

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,

src/fixed_point_spectral_radius.jl

@@ -0,0 +1,242 @@
+"""
+    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. 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.
+"""
+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, condest::Float64)
+
+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 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 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.
+
+# 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(residual, 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)
+    # 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,
+    )
+    ρ = isempty(vals) ? NaN : abs(vals[1])
+    condest = KLU.condest(F)
+    return ρ, info, condest
+end