diff --git a/.gitignore b/.gitignore index 7c316c72b..0c6454293 100644 --- a/.gitignore +++ b/.gitignore @@ -38,6 +38,7 @@ docs/site/ *.ipynb Manifest.toml +scripts/benchmarks/*.csv .vscode *.h5 data diff --git a/scripts/benchmarks/Project.toml b/scripts/benchmarks/Project.toml new file mode 100644 index 000000000..74a706228 --- /dev/null +++ b/scripts/benchmarks/Project.toml @@ -0,0 +1,6 @@ +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" +PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" +PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" diff --git a/scripts/benchmarks/method_comparison.jl b/scripts/benchmarks/method_comparison.jl new file mode 100644 index 000000000..01d3453ed --- /dev/null +++ b/scripts/benchmarks/method_comparison.jl @@ -0,0 +1,423 @@ +""" +Head-to-head comparison of all AC power flow methods across grid sizes and +initial-condition quality. Captures convergence, iteration count, initial / +final residuals, and wall-clock time. Set `OVERWRITE_NON_CONVERGED = false` +in `src/definitions.jl`, else final residual size will show NaN whenever +convergence fails. + +Perturbation model: for each trial, the true NR solution is perturbed by adding +`K * scale * randn()` independently to each state variable, where `scale` depends +on the physical quantity: `VM_SCALE` for voltage magnitudes (PQ buses), `VA_SCALE` +for voltage angles (PQ and PV buses), and `PQ_SCALE` for active/reactive power +(REF and PV buses). The same K value does NOT produce the same total perturbation +norm across systems: the expected L2 norm of the perturbation is `~K sqrt(n)` for a +system of size `n`, with a constant that depends on the bus type distribution. +**IMPORTANT**: K values are comparable across solvers on the same system, but not +across systems. + +Usage: + julia --project=scripts/benchmarks scripts/benchmarks/method_comparison.jl +""" + +using PowerFlows, PowerSystemCaseBuilder, PowerSystems +using Logging, LinearAlgebra, Random, DataFrames, Printf, Statistics, CSV + +const PSB = PowerSystemCaseBuilder +const PSY = PowerSystems +const PF = PowerFlows + +# ───────────────────────────────────────────────────────────────────────────── +# Custom logger that intercepts the @info messages emitted by the solver +# and extracts iteration count and residual norms. +# ───────────────────────────────────────────────────────────────────────────── +mutable struct MetricsCapture + iterations::Int + initial_residual_L2::Float64 + initial_residual_Linf::Float64 + final_residual_L2::Float64 + final_residual_Linf::Float64 + converged_from_log::Union{Bool, Nothing} +end +MetricsCapture() = MetricsCapture(0, NaN, NaN, NaN, NaN, nothing) + +struct CapturingLogger <: AbstractLogger + metrics::MetricsCapture + min_level::LogLevel +end +CapturingLogger(m::MetricsCapture) = CapturingLogger(m, Logging.Debug) + +Logging.min_enabled_level(l::CapturingLogger) = l.min_level +Logging.shouldlog(l::CapturingLogger, level, _module, group, id) = true +Logging.catch_exceptions(::CapturingLogger) = true + +function Logging.handle_message(l::CapturingLogger, level, message, _module, group, id, + filepath, line; kwargs...) + msg = string(message) + + # "Initial residual size: 1.23 L2, 4.56 L∞" + m = match(r"Initial residual size:\s*([\d.eE+-]+)\s*L2,\s*([\d.eE+-]+)\s*L", msg) + if m !== nothing + l.metrics.initial_residual_L2 = parse(Float64, m.captures[1]) + l.metrics.initial_residual_Linf = parse(Float64, m.captures[2]) + end + + # "Final residual size: 1.23 L2, 4.56 L∞." + m = match(r"Final residual size:\s*([\d.eE+-]+)\s*L2,\s*([\d.eE+-]+)\s*L", msg) + if m !== nothing + l.metrics.final_residual_L2 = parse(Float64, m.captures[1]) + l.metrics.final_residual_Linf = parse(Float64, m.captures[2]) + end + + # "The NewtonRaphsonACPowerFlow solver converged after 5 iterations." + m = match(r"solver converged after (\d+) iteration", msg) + if m !== nothing + l.metrics.iterations = parse(Int, m.captures[1]) + l.metrics.converged_from_log = true + end + + # "The … solver failed to converge after N iterations." + m = match(r"solver failed to converge after (\d+) iteration", msg) + if m !== nothing + l.metrics.iterations = parse(Int, m.captures[1]) + l.metrics.converged_from_log = false + elseif occursin("failed to converge", msg) + l.metrics.converged_from_log = false + end + + return nothing +end + +# ───────────────────────────────────────────────────────────────────────────── +# Configuration +# ───────────────────────────────────────────────────────────────────────────── +const SYSTEMS = [ + (PSB.PSITestSystems, "c_sys14", "14-bus", Dict(:add_forecasts => false)), + (PSB.MatpowerTestSystems, "matpower_ACTIVSg2000_sys", "2000-bus", Dict{Symbol, Any}()), + (PSB.MatpowerTestSystems, "matpower_ACTIVSg10k_sys", "10k-bus", Dict{Symbol, Any}()), +] + +const SOLVERS = [ + (PF.NewtonRaphsonACPowerFlow, "NR", Dict{Symbol, Any}()), + (PF.NewtonRaphsonACPowerFlow, "NR+Iwamoto", Dict{Symbol, Any}(:iwamoto => true)), + (PF.TrustRegionACPowerFlow, "TR+Iwamoto", Dict{Symbol, Any}(:iwamoto => true)), + (PF.LevenbergMarquardtACPowerFlow, "LM", Dict{Symbol, Any}()), + (PF.RobustHomotopyPowerFlow, "Homotopy", Dict{Symbol, Any}()), + (PF.GradientDescentACPowerFlow, "GradDescent", Dict{Symbol, Any}()), +] + +# Perturbation magnitudes (multiplied by randn to produce x0 offsets). +# Two special cases run once (deterministic): +# "system data" — use system setpoints as x0 (no override, whatever's in the System object) +# "flat start" — true flat start: Vm=1.0 pu, Va=0.0 rad for all buses +const PERTURBATIONS = [ + (NaN, "system data"), # default initialization (no x0 override) + (-1.0, "flat start"), # true flat start: Vm=1, Va=0 (sentinel value) + (0.001, "K=0.001"), + (0.005, "K=0.005"), + (0.01, "K=0.01"), + (0.025, "K=0.025"), + (0.05, "K=0.05"), + (0.075, "K=0.075"), + (0.1, "K=0.1"), + (0.15, "K=0.15"), + (0.2, "K=0.2"), + (0.25, "K=0.25"), + (0.3, "K=0.3"), +] + +const N_TRIALS = 20 # repeats per (system, solver, perturbation) triple + +# ───────────────────────────────────────────────────────────────────────────── +# Helpers +# ───────────────────────────────────────────────────────────────────────────── + +"""Suppress all logs below `Error` for the duration of `f()`.""" +function quietly(f) + with_logger(SimpleLogger(stderr, Logging.Error)) do + f() + end +end + +"""Get the true solution for a system by solving with NR from flat start.""" +function get_true_solution(sys) + pf = ACPowerFlow(; correct_bustypes = true) + data = quietly(() -> PF.PowerFlowData(pf, sys)) + n = size(PF.get_bus_type(data), 1) + bus_types = PF.get_bus_type(data)[:, 1] + quietly(() -> PF.solve_power_flow!(data; pf = pf)) + x_solved = zeros(2n) + PF.update_state!(x_solved, data, 1) + return x_solved, n, bus_types +end + +# Base perturbation scales per physical quantity. +const VM_SCALE = 0.1 # pu — voltage magnitude +const VA_SCALE = 0.2 # rad — voltage angle (~11 degrees) +const PQ_SCALE = 0.1 # pu — active/reactive power + +"""Build a perturbed x0 from x_solved. All state variables are perturbed: +Vm and Va at PQ buses, Q and Va at PV buses, P and Q at REF buses. +K is a multiplier on the base scales.""" +function perturb_x0(x_solved, bus_types, K, rng) + n = length(bus_types) + x0 = copy(x_solved) + for i in 1:n + bt = bus_types[i] + if bt == PSY.ACBusTypes.PQ + x0[2i - 1] += K * VM_SCALE * randn(rng) # Vm + x0[2i] += K * VA_SCALE * randn(rng) # Va + elseif bt == PSY.ACBusTypes.PV + x0[2i - 1] += K * PQ_SCALE * randn(rng) # Q + x0[2i] += K * VA_SCALE * randn(rng) # Va + elseif bt == PSY.ACBusTypes.REF + x0[2i - 1] += K * PQ_SCALE * randn(rng) # P + x0[2i] += K * PQ_SCALE * randn(rng) # Q + end + end + return x0 +end + +"""Run a single benchmark trial. Returns a NamedTuple of metrics.""" +function run_trial(sys, solver_type, solver_settings, x_solved, n, bus_types, K; + seed::Int = 42) + # Build PF object + pf = ACPowerFlow{solver_type}(; + correct_bustypes = true, + solver_settings = solver_settings, + ) + data = quietly(() -> PF.PowerFlowData(pf, sys)) + + # Build x0 + kwargs = Dict{Symbol, Any}() + if isnan(K) + # "system data" — no x0 override, use whatever's in the System object + elseif K < 0.0 + # "flat start" — true flat start: Vm=1.0, Va=0.0 for all buses + x0 = zeros(2n) + for i in 1:length(bus_types) + x0[2i - 1] = 1.0 # Vm = 1.0 pu (or P/Q = 1.0 for REF/PV, but solver overwrites) + x0[2i] = 0.0 # Va = 0.0 rad + end + kwargs[:x0] = x0 + else + rng = Random.MersenneTwister(seed) + x0 = perturb_x0(x_solved, bus_types, K, rng) + kwargs[:x0] = x0 + end + + # Compute initial power flow residual at the actual starting point, + # consistently across all methods (including homotopy). + residual_obj = PF.ACPowerFlowResidual(data, 1) + if haskey(kwargs, :x0) + residual_obj(kwargs[:x0], 1) + else + x0_default = PF.calculate_x0(data, 1) + residual_obj(x0_default, 1) + end + init_res_L2 = norm(residual_obj.Rv, 2) + init_res_Linf = norm(residual_obj.Rv, Inf) + + # Solve under capturing logger + metrics = MetricsCapture() + logger = CapturingLogger(metrics) + local converged::Bool + local elapsed::Float64 + + elapsed = @elapsed begin + converged = with_logger(logger) do + PF.solve_power_flow!(data; pf = pf, kwargs...) + end + end + + # Extract final state and compare to reference solution + x_final = zeros(2n) + PF.update_state!(x_final, data, 1) + solution_dist_Linf = norm(x_final .- x_solved, Inf) + solution_dist_L2 = norm(x_final .- x_solved, 2) + + # Compute final power flow residual directly — don't rely on log capture, + # which may be missing (homotopy) or absent on non-convergence. + residual_obj(x_final, 1) + final_res_L2 = norm(residual_obj.Rv, 2) + final_res_Linf = norm(residual_obj.Rv, Inf) + + iters = metrics.iterations + + # Residual reduction ratio: how much did the residual shrink? + res_reduction = + (init_res_Linf > 0 && final_res_Linf > 0) ? + init_res_Linf / final_res_Linf : NaN + + return ( + converged = converged, + iterations = iters, + initial_residual_L2 = init_res_L2, + initial_residual_Linf = init_res_Linf, + final_residual_L2 = final_res_L2, + final_residual_Linf = final_res_Linf, + elapsed = elapsed, + solution_dist_Linf = solution_dist_Linf, + solution_dist_L2 = solution_dist_L2, + res_reduction = res_reduction, + ) +end + +# ───────────────────────────────────────────────────────────────────────────── +# Main +# ───────────────────────────────────────────────────────────────────────────── +function main() + results = DataFrame(; + system = String[], + solver = String[], + perturbation = String[], + trial = Int[], + converged = Bool[], + iterations = Int[], + initial_res_L2 = Float64[], + initial_res_Linf = Float64[], + final_res_L2 = Float64[], + final_res_Linf = Float64[], + elapsed_s = Float64[], + sol_dist_Linf = Float64[], + sol_dist_L2 = Float64[], + res_reduction = Float64[], + ) + + for (group, sysname, label, build_kwargs) in SYSTEMS + println("\n", "="^70) + println("System: $label ($sysname)") + println("="^70) + + sys = quietly(() -> PSB.build_system(group, sysname; build_kwargs...)) + x_solved, n, bus_types = get_true_solution(sys) + println(" Buses: $n, true solution obtained.\n") + + for (solver_type, solver_label, solver_settings) in SOLVERS + gave_up = false + for (K, perturb_label) in PERTURBATIONS + if gave_up + println( + " %-14s %-12s ⏭ skipped (solver gave up at lower K)" |> + s -> @sprintf( + " %-14s %-12s ⏭ skipped (solver gave up at lower K)", + solver_label, perturb_label), + ) + continue + end + + n_converged = 0 + n_trials = (isnan(K) || K < 0.0) ? 1 : N_TRIALS + for trial in 1:n_trials + seed = trial + + local result + try + result = run_trial( + sys, solver_type, solver_settings, + x_solved, n, bus_types, K; seed = seed, + ) + catch e + if e isa InterruptException + println( + "\n ⏭ Interrupted, skipping remaining trials for $solver_label / $perturb_label", + ) + break + end + @warn "FAILED: $label / $solver_label / $perturb_label / trial $trial" exception = + (e, catch_backtrace()) + result = ( + converged = false, iterations = 0, + initial_residual_L2 = NaN, initial_residual_Linf = NaN, + final_residual_L2 = NaN, final_residual_Linf = NaN, + elapsed = NaN, + solution_dist_Linf = NaN, solution_dist_L2 = NaN, + res_reduction = NaN, + ) + end + + n_converged += result.converged + + push!( + results, + ( + label, solver_label, perturb_label, trial, + result.converged, result.iterations, + result.initial_residual_L2, result.initial_residual_Linf, + result.final_residual_L2, result.final_residual_Linf, + result.elapsed, + result.solution_dist_Linf, result.solution_dist_L2, + result.res_reduction, + ), + ) + + status = result.converged ? "OK" : "FAIL" + dist_flag = result.solution_dist_Linf > 1e-4 ? " ⚠DIFF" : "" + reduc_str = if isnan(result.res_reduction) + "N/A" + else + @sprintf("%.1e", result.res_reduction) + end + @printf( + " %-14s %-12s trial %d %4s %3d iter %.4f s init‖r‖∞=%.2e final‖r‖∞=%.2e ‖Δx‖∞=%.2e reduce=%s%s\n", + solver_label, perturb_label, trial, status, + result.iterations, result.elapsed, + result.initial_residual_Linf, result.final_residual_Linf, + result.solution_dist_Linf, reduc_str, dist_flag) + end + + # If no trials converged at this K, skip all higher K values. + # Exception: GradDescent is not expected to converge; skip the early exit. + if n_converged == 0 && !isnan(K) && K > 0.0 && + solver_type != PF.GradientDescentACPowerFlow + gave_up = true + println( + " %-14s %-12s → 0/$N_TRIALS converged, skipping higher K" |> + s -> + @sprintf(" %-14s %-12s → 0/%d converged, skipping higher K", + solver_label, perturb_label, N_TRIALS), + ) + end + end + end + end + + # ── Summary table ────────────────────────────────────────────────────── + println("\n\n", "="^70) + println("SUMMARY (median over $N_TRIALS trials)") + println("="^70) + + summary = combine( + groupby(results, [:system, :solver, :perturbation]), + :converged => (x -> all(x)) => :all_converged, + :iterations => median => :med_iter, + :elapsed_s => median => :med_time, + :initial_res_Linf => median => :med_init_Linf, + :final_res_Linf => median => :med_final_Linf, + :sol_dist_Linf => maximum => :max_sol_dist, + :res_reduction => median => :med_reduction, + ) + + @printf( + "\n%-10s %-14s %-12s Conv Iter Time(s) Init‖r‖∞ Final‖r‖∞ max‖Δx‖∞ Reduction\n", + "System", "Solver", "Perturbation") + println("-"^120) + for row in eachrow(summary) + conv = row.all_converged ? "YES" : "NO" + reduc_str = + isnan(row.med_reduction) ? " N/A" : @sprintf("%10.1e", row.med_reduction) + @printf( + "%-10s %-14s %-12s %-4s %4.0f %8.4f %10.2e %10.2e %10.2e %s\n", + row.system, row.solver, row.perturbation, conv, + row.med_iter, row.med_time, row.med_init_Linf, row.med_final_Linf, + row.max_sol_dist, reduc_str) + end + + # Save full results + outfile = joinpath(@__DIR__, "method_comparison_results.csv") + CSV.write(outfile, results) + println("\nFull results saved to $outfile") + + return results +end + +main() diff --git a/src/RobustHomotopy/robust_homotopy_method.jl b/src/RobustHomotopy/robust_homotopy_method.jl index ae598b99f..62b51ebed 100644 --- a/src/RobustHomotopy/robust_homotopy_method.jl +++ b/src/RobustHomotopy/robust_homotopy_method.jl @@ -18,17 +18,25 @@ function _newton_power_flow(pf::ACPolarPowerFlow{<:RobustHomotopyPowerFlow}, symbolic_factor!(hSolver, homHess.Hv) success = true + total_iters = 0 while true # go onto next t_k even if search doesn't terminate within max iterations. - converged_t_k, _ = _second_order_newton(homHess, t_k, time_step, x, hSolver) + converged_t_k, iters = _second_order_newton(homHess, t_k, time_step, x, hSolver) + total_iters += iters if t_k == 1.0 success = converged_t_k break end t_k = min(t_k + Δt_k, 1.0) end + r_L2 = norm(homHess.pfResidual.Rv, 2) + r_Linf = norm(homHess.pfResidual.Rv, Inf) + @info("Final residual size: $(r_L2) L2, $(r_Linf) L∞.") if !success - @warn "RobustHomotopyPowerFlow failed to find a solution" + @error( + "The RobustHomotopyPowerFlow solver failed to converge after $total_iters iterations." + ) else + @info("The RobustHomotopyPowerFlow solver converged after $total_iters iterations.") if get_calculate_loss_factors(data) _calculate_loss_factors(data, homHess.J.Jv, time_step) end diff --git a/src/definitions.jl b/src/definitions.jl index 3923b6ec3..bafa4783a 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -39,8 +39,9 @@ const DEFAULT_IWAMOTO_FALLBACK = true # when a trust region step is rejected, tr const PF_MAX_LOG = 10 # only used for Levenberg-Maquardt const DEFAULT_λ_0 = 1e-5 -# input is mix of powers (100 MW), voltages (0.8-1.2), and angles (-π/4 to π/4). -const DEFAULT_MAX_TEST_λs = 50 # give up after increasing damping factor 50 times. +# Upper bound on the LM damping factor μ. μ only grows on rejected steps, so +# hitting this cap is a divergence signal — the solver aborts with an error. +const DEFAULT_μ_MAX = 1e8 const DEFAULT_Δt_k = 0.2 diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index 00dd9af33..b04e24d76 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -1,8 +1,87 @@ -"""Driver for the LevenbergMarquardtACPowerFlow method: sets up the data -structures (e.g. residual), runs the power flow method via calling `_run_power_flow_method` +"""Pre-allocated workspace for the Levenberg-Marquardt solver. + +Holds the augmented matrix `[J; √λ·I]` with a fixed sparsity pattern, +a mapping to update its entries in-place, and a cached QR factorization.""" +mutable struct LMWorkspace + A::SparseMatrixCSC{Float64, Int64} + # PERF: if we're doing dense access, I'd expect a bitmask to be faster. + # Or only store λ_diag_indices (sparse) and infer J entries as the complement. + # Indices into A.nzval for the J block entries (same order as J.Jv.nzval) + j_nzval_indices::Vector{Int} + # Indices into A.nzval for the √λ diagonal entries (length n) + λ_diag_indices::Vector{Int} + # PERF: have not investigated other factorization methods. SPQR does not allow + # for in-place updates or re-use of symbolic factorization, but with non + # square matrices--J^T J form is more unstable--we don't have a lot of options. + # Cached QR factorization + F::SparseArrays.SPQR.QRSparse{Float64, Int64} + # Preallocated augmented RHS [-Rv; 0] (length m + n); bottom n stay zero. + b::Vector{Float64} +end + +"""Build the augmented matrix `[J; I]` once, recording which `A.nzval` entries +correspond to J values vs the λ diagonal.""" +function LMWorkspace(Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}) + m, n = size(Jv) + + # Convert J to Int64 indices for SPQR compatibility, then vcat with identity. + Jv64 = SparseMatrixCSC{Float64, Int64}( + Jv.m, Jv.n, + Vector{Int64}(Jv.colptr), + Vector{Int64}(Jv.rowval), + copy(Jv.nzval), + ) + Iλ = sparse(Int64.(1:n), Int64.(1:n), ones(n), n, n) + A = vcat(Jv64, Iλ) + + # Identify which A.nzval entries come from J vs the diagonal. + j_nzval_indices = Vector{Int}(undef, length(Jv.nzval)) + λ_diag_indices = Vector{Int}(undef, n) + + j_idx = 0 + for col in 1:n + for a_idx in A.colptr[col]:(A.colptr[col + 1] - 1) + row = A.rowval[a_idx] + if row <= m + j_idx += 1 + j_nzval_indices[j_idx] = a_idx + elseif row == m + col + λ_diag_indices[col] = a_idx + end + end + end + @assert j_idx == length(Jv.nzval) "Expected $(length(Jv.nzval)) J entries, found $j_idx" + + b = zeros(m + n) + F = LinearAlgebra.qr(A) + + return LMWorkspace(A, j_nzval_indices, λ_diag_indices, F, b) +end + +"""Copy current Jacobian values into the augmented matrix.""" +function copy_jacobian!(ws::LMWorkspace, Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}) + nzv = Jv.nzval + for (i, a_idx) in enumerate(ws.j_nzval_indices) + ws.A.nzval[a_idx] = nzv[i] + end + return +end + +"""Update the √λ diagonal and re-factorize.""" +function update_lambda!(ws::LMWorkspace, λ::Float64) + sqrtλ = sqrt(λ) + for a_idx in ws.λ_diag_indices + ws.A.nzval[a_idx] = sqrtλ + end + ws.F = LinearAlgebra.qr(ws.A) + return +end + +"""Driver for the LevenbergMarquardtACPowerFlow method: sets up the data +structures (e.g. residual), runs the power flow method via calling `_run_power_flow_method` on them, then handles post-processing (e.g. loss factors).""" function _newton_power_flow( - pf::ACPolarPowerFlow{LevenbergMarquardtACPowerFlow}, + pf::AbstractACPowerFlow{LevenbergMarquardtACPowerFlow}, data::ACPowerFlowData, time_step::Int64; tol::Float64 = DEFAULT_NR_TOL, @@ -10,23 +89,32 @@ function _newton_power_flow( validate_voltage_magnitudes::Bool = DEFAULT_VALIDATE_VOLTAGES, vm_validation_range::MinMax = DEFAULT_VALIDATION_RANGE, λ_0::Float64 = DEFAULT_λ_0, - maxTestλs::Int = DEFAULT_MAX_TEST_λs, + x0::Union{Vector{Float64}, Nothing} = nothing, _ignored..., ) + init_kwargs = if isnothing(x0) + (; validate_voltage_magnitudes, vm_validation_range) + else + (; validate_voltage_magnitudes, vm_validation_range, x0) + end residual, J, x0 = initialize_power_flow_variables( - pf, data, time_step; - validate_voltage_magnitudes, vm_validation_range) + pf, data, time_step; init_kwargs...) converged = norm(residual.Rv, Inf) < tol i = 0 if !converged + ws = LMWorkspace(J.Jv) converged, i = _run_power_flow_method( time_step, x0, residual, - J; - tol, maxIterations, λ_0, maxTestλs, + J, + ws; + tol, maxIterations, λ_0, ) end + # x0 was mutated in place to the converged state by _run_power_flow_method + # (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) end @@ -34,84 +122,84 @@ end function _run_power_flow_method( time_step::Int, x::Vector{Float64}, - residual::ACPowerFlowResidual, - J::ACPowerFlowJacobian; + residual::Union{ACPowerFlowResidual, ACRectangularCIResidual}, + J::Union{ACPowerFlowJacobian, ACRectangularCIJacobian}, + ws::LMWorkspace; maxIterations::Int = DEFAULT_NR_MAX_ITER, tol::Float64 = DEFAULT_NR_TOL, λ_0::Float64 = DEFAULT_λ_0, - maxTestλs::Int = DEFAULT_MAX_TEST_λs, _ignored..., ) - λ::Float64 = λ_0 + μ::Float64 = λ_0 + λ::Float64 = 0.0 i, converged = 0, false residual(x, time_step) resSize = dot(residual.Rv, residual.Rv) linf = norm(residual.Rv, Inf) @debug "initially: sum of squares $(siground(resSize)), L ∞ norm $(siground(linf)), λ = $λ" - while i < maxIterations && !converged && !isnan(λ) - λ = update_damping_factor!(x, residual, J, time_step, maxTestλs) - converged = !isnan(λ) && norm(residual.Rv, Inf) < tol + while i < maxIterations && !converged && isfinite(λ) && μ < DEFAULT_μ_MAX + λ, μ = update_damping_factor!(x, residual, J, μ, time_step, ws) + converged = isfinite(λ) && norm(residual.Rv, Inf) < tol i += 1 end - if isnan(λ) - @error "λ is NaN" - elseif i == maxIterations - @error "The LevenbergMarquardtACPowerFlow solver didn't coverge in $maxIterations iterations." + if !converged + if !isfinite(λ) + @error "λ is not finite ($(λ))" + elseif μ >= DEFAULT_μ_MAX + @error "The LevenbergMarquardtACPowerFlow damping factor μ hit the cap (DEFAULT_μ_MAX=$(DEFAULT_μ_MAX)) after $i iterations; aborting (likely divergence)." + elseif i == maxIterations + @error "The LevenbergMarquardtACPowerFlow solver didn't coverge in $maxIterations iterations." + end end return converged, i end -# LM implementation and parameters values largely based upon this paper: -# # https://www.sciencedirect.com/science/article/pii/S0142061518336342 +# LM implementation based on standard Levenberg-Marquardt method. +# See Nocedal & Wright (2006), sections 10.3 and 11.2. + +"""Compute one LM trial step. Assumes `residual` and `J` are already evaluated +at `x` by the caller. Returns the gain ratio ρ.""" function compute_error( x::Vector{Float64}, - residual::ACPowerFlowResidual, - J::ACPowerFlowJacobian, + residual::Union{ACPowerFlowResidual, ACRectangularCIResidual}, + J::Union{ACPowerFlowJacobian, ACRectangularCIJacobian}, λ::Float64, time_step::Int, residualSize::Float64, + ws::LMWorkspace, ) - residual(x, time_step) # M(x_c) - J(time_step) - - n = size(J.Jv, 2) - Iλ = sparse(1:n, 1:n, sqrt(λ) .* ones(n), n, n) # less error-prone compared to A = vcat(J.Jv, sqrt(λ) * sparse(LinearAlgebra.I, size(J.Jv))) - A = [J.Jv; Iλ] + # Update augmented matrix with current J values and λ, then factorize once. + copy_jacobian!(ws, J.Jv) + update_lambda!(ws, λ) - b_x = vcat(-residual.Rv, zeros(size(J.Jv, 2))) - Δx = A \ b_x + m = length(residual.Rv) + @assert m == length(ws.b) - size(J.Jv, 2) "residual/J size mismatch vs preallocated LM buffer (m=$m, buf=$(length(ws.b)), n=$(size(J.Jv, 2)))" + @views ws.b[1:m] .= .-residual.Rv # bottom n entries stay zero from construction + Δx = ws.F \ ws.b temp_x = residual.Rv .+ J.Jv * Δx x_trial = x .+ Δx - residual(x_trial, time_step) # M(y_c) - - b_y = vcat(-residual.Rv, zeros(size(J.Jv, 2))) - Δy = A \ b_y - temp_y = residual.Rv .+ J.Jv * Δy - - newResidualSize_y = dot(residual.Rv, residual.Rv) - residual(x_trial .+ Δy, time_step) # M(x_c+Δx+Δy) + residual(x_trial, time_step) # M(x_c + Δx) newResidualSize = dot(residual.Rv, residual.Rv) - b_z = vcat(-residual.Rv, zeros(size(J.Jv, 2))) - Δz = A \ b_z - temp_z = residual.Rv .+ J.Jv * Δz - newResidualSize_z = dot(residual.Rv, residual.Rv) + predicted_reduction = residualSize - dot(temp_x, temp_x) + actual_reduction = residualSize - newResidualSize - residual(x_trial .+ Δy .+ Δz, time_step) # M(x_c+Δx+Δy+Δz) - newResidualSize = dot(residual.Rv, residual.Rv) + # Guard against zero/negative predicted reduction. + if predicted_reduction <= 0.0 || !isfinite(predicted_reduction) + residual(x, time_step) + return 0.0 + end - predicted_reduction = ( - residualSize - dot(temp_x, temp_x) + newResidualSize_y - dot(temp_y, temp_y) + - newResidualSize_z - dot(temp_z, temp_z) - ) - actual_reduction = (residualSize - newResidualSize) ρ = actual_reduction / predicted_reduction if ρ > 1e-4 - x .+= (Δx .+ Δy .+ Δz) + x .+= Δx + else + # Bad step: restore data state to match x (not x_trial). + residual(x, time_step) end return ρ @@ -119,39 +207,26 @@ end function update_damping_factor!( x::Vector{Float64}, - residual::ACPowerFlowResidual, - J::ACPowerFlowJacobian, + residual::Union{ACPowerFlowResidual, ACRectangularCIResidual}, + J::Union{ACPowerFlowJacobian, ACRectangularCIJacobian}, + μ::Float64, time_step::Int, - maxTestλs::Int, + ws::LMWorkspace, ) residual(x, time_step) residualSize = dot(residual.Rv, residual.Rv) J(time_step) - # Now update \lambda - test_lambda::Int = 1 - # https://www.sciencedirect.com/science/article/pii/S0142061518336342 - λ = DEFAULT_λ_0 * sqrt(residualSize) - while test_lambda < maxTestλs - ρ = compute_error(x, residual, J, λ, time_step, residualSize) - if ρ > 0.75 # good step - factor = 4.0 - λ_temp = λ / factor - λ = max(λ_temp, 1.0e-8) - break - elseif ρ >= 0.25 # okay step - break - else # bad step - factor = 4.0 - λ *= factor - end - - test_lambda += 1 - end - if test_lambda == maxTestλs - @error "Unable to improve: gave up after increasing damping factor $maxTestλs times." - λ = NaN + λ = μ * sqrt(residualSize) + ρ = compute_error(x, residual, J, λ, time_step, residualSize, ws) + coef = 4.0 + if ρ > 0.75 + μ = max(μ / coef, 1e-8) + elseif ρ >= 0.25 + # intentional no-op + else + μ = min(μ * coef, DEFAULT_μ_MAX) end - return λ + return (λ, μ) end diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 06f85b99f..88c81b213 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -725,7 +725,7 @@ function _finalize_power_flow( end return true end - @error("The $solver_name solver failed to converge.") + @error("The $solver_name solver failed to converge after $i iterations.") return false end diff --git a/src/power_flow_types.jl b/src/power_flow_types.jl index 2dace7dec..a4115e28f 100644 --- a/src/power_flow_types.jl +++ b/src/power_flow_types.jl @@ -81,9 +81,12 @@ struct TrustRegionACPowerFlow <: ACPowerFlowSolverType end An [`ACPowerFlowSolverType`](@ref) corresponding to the [Levenberg-Marquardt](https://en.wikipedia.org/wiki/Levenberg–Marquardt_algorithm) iterative method. This is more robust than the basic Newton-Raphson method, but also more computationally -intensive. Due to the difficulty of tuning meta parameters, this method may occasionally +intensive. Due to the difficulty of tuning meta parameters, this method may occasionally fail to converge where other methods would succeed. +Works with both the polar ([`ACPolarPowerFlow`](@ref)) and rectangular +current-injection ([`ACRectangularPowerFlow`](@ref)) formulations. + See also: [`ACPowerFlow`](@ref). """ struct LevenbergMarquardtACPowerFlow <: ACPowerFlowSolverType end @@ -278,10 +281,10 @@ State per bus: PQ `(eᵢ, fᵢ)`, PV `(eᵢ, fᵢ, Qᵢ)`, REF `(P_genᵢ, Q_gen `ΔIᵢ = I_specᵢ − Y_bus·V`. Off-diagonal Jacobian blocks ≡ Y_bus 2×2 real blocks and are constant across iterations. -`ACSolver` defaults to [`NewtonRaphsonACPowerFlow`](@ref); only -[`NewtonRaphsonACPowerFlow`](@ref) and [`TrustRegionACPowerFlow`](@ref) are -supported. Levenberg-Marquardt, Robust Homotopy, and Gradient Descent operate -on the polar formulation only and are rejected at construction. +`ACSolver` defaults to [`NewtonRaphsonACPowerFlow`](@ref). Supported solvers: +[`NewtonRaphsonACPowerFlow`](@ref), [`TrustRegionACPowerFlow`](@ref), and +[`LevenbergMarquardtACPowerFlow`](@ref). Robust Homotopy and Gradient Descent +operate on the polar formulation only and are rejected at construction. Unlike [`ACPolarPowerFlow`](@ref), this model has no `calculate_voltage_stability_factors`, `calculate_loss_factors`, or @@ -339,19 +342,16 @@ function ACRectangularPowerFlow{ACSolver}(; solver_settings::Dict{Symbol, Any} = Dict{Symbol, Any}(), ) where {ACSolver <: ACPowerFlowSolverType} if ACSolver <: Union{ - LevenbergMarquardtACPowerFlow, RobustHomotopyPowerFlow, GradientDescentACPowerFlow, } throw( ArgumentError( "$(ACSolver) is not supported by ACRectangularPowerFlow. " * - "Levenberg-Marquardt, Robust Homotopy, and Gradient Descent " * - "operate on the polar formulation only. Use " * - "ACRectangularPowerFlow{NewtonRaphsonACPowerFlow} or " * - "{TrustRegionACPowerFlow}, or run the solver on " * - "ACPolarPowerFlow. (Rectangular LM is tracked as a separate " * - "follow-up project.)", + "Robust Homotopy and Gradient Descent operate on the polar " * + "formulation only. Use ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}, " * + "{TrustRegionACPowerFlow}, or {LevenbergMarquardtACPowerFlow}, " * + "or run the solver on ACPolarPowerFlow.", ), ) end diff --git a/test/test_rectangular_ci_power_flow.jl b/test/test_rectangular_ci_power_flow.jl index 44d4a6fad..80528157c 100644 --- a/test/test_rectangular_ci_power_flow.jl +++ b/test/test_rectangular_ci_power_flow.jl @@ -32,6 +32,20 @@ end ) end +@testset "Rectangular CI Power Flow (LM): non-convergence returns missing" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) + # maxIterations = 1 from flat start cannot converge c_sys14; the solver must + # report non-convergence (results = missing) rather than error or hang. + pf = ACRectangularPowerFlow{LevenbergMarquardtACPowerFlow}(; + solver_settings = merge(_rect_pf_settings(), + Dict{Symbol, Any}(:maxIterations => 1))) + @test_logs( + (:error, r".*solver failed to converge"), + match_mode = :any, + @test ismissing(solve_power_flow(pf, sys)) + ) +end + @testset "Rectangular CI Power Flow: parity with polar NR" begin fixtures = [ ("c_sys5", false), @@ -75,9 +89,11 @@ end @test ACPolarPowerFlow{NewtonRaphsonACPowerFlow}(; calculate_voltage_stability_factors = true) isa ACPolarPowerFlow # Polar-only solvers rejected at construction. - @test_throws ArgumentError ACRectangularPowerFlow{LevenbergMarquardtACPowerFlow}() @test_throws ArgumentError ACRectangularPowerFlow{RobustHomotopyPowerFlow}() @test_throws ArgumentError ACRectangularPowerFlow{GradientDescentACPowerFlow}() + # Levenberg-Marquardt is now supported on the rectangular formulation. + @test ACRectangularPowerFlow{LevenbergMarquardtACPowerFlow}() isa + ACRectangularPowerFlow end @testset "Rectangular CI Power Flow: step strategy variants" begin @@ -93,6 +109,8 @@ end ("Trust Region", TrustRegionACPowerFlow, Dict{Symbol, Any}()), ("TR + Iwamoto FB", TrustRegionACPowerFlow, Dict{Symbol, Any}(:iwamoto_fallback => true)), + ("Levenberg-Marquardt", LevenbergMarquardtACPowerFlow, + Dict{Symbol, Any}()), ] for (name, with_forecasts) in fixtures @testset "$name" begin @@ -123,6 +141,65 @@ end end end +@testset "Rectangular CI: LM matches polar LM" begin + for name in ("c_sys5", "c_sys14") + @testset "$name" begin + sys = PSB.build_system(PSB.PSITestSystems, name; add_forecasts = false) + pf_polar = ACPolarPowerFlow{LevenbergMarquardtACPowerFlow}() + res_polar = solve_power_flow(pf_polar, deepcopy(sys)) + pf_rect = ACRectangularPowerFlow{LevenbergMarquardtACPowerFlow}(; + solver_settings = _rect_pf_settings()) + res_rect = solve_power_flow(pf_rect, deepcopy(sys)) + @test res_rect !== missing + @test maximum( + abs.(res_polar["bus_results"].Vm - res_rect["bus_results"].Vm), + ) < 1e-7 + @test maximum( + abs.(res_polar["bus_results"].θ - res_rect["bus_results"].θ), + ) < 1e-7 + end + end +end + +@testset "ACTIVSg2000 (LM): polar and rectangular match polar NR" begin + sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_ACTIVSg2000_sys") + PSY.set_units_base_system!(sys, "SYSTEM_BASE") + + # Tight tolerance + generous iteration budget: LM refactorizes the sparse QR + # every iteration and needs more iterations than NR on a 2000-bus system. + lm_settings = Dict{Symbol, Any}(:tol => 1e-10, :maxIterations => 100) + ref_settings = Dict{Symbol, Any}(:tol => 1e-10) + + # Reference: Newton-Raphson on the polar formulation (trusted for ACTIVSg2000 + # elsewhere in the suite). + pf_ref = ACPolarPowerFlow{NewtonRaphsonACPowerFlow}(; + correct_bustypes = true, solver_settings = ref_settings) + res_ref = solve_power_flow(pf_ref, sys) + + pf_lm_polar = ACPolarPowerFlow{LevenbergMarquardtACPowerFlow}(; + correct_bustypes = true, solver_settings = lm_settings) + res_lm_polar = solve_power_flow(pf_lm_polar, sys) + + pf_lm_rect = ACRectangularPowerFlow{LevenbergMarquardtACPowerFlow}(; + correct_bustypes = true, + solver_settings = merge(_rect_pf_settings(), lm_settings)) + res_lm_rect = solve_power_flow(pf_lm_rect, sys) + + @test res_lm_polar !== missing + @test res_lm_rect !== missing + + # Polar LM and rectangular LM both reproduce the NR reference solution. + @test norm(res_lm_polar["bus_results"].Vm .- res_ref["bus_results"].Vm, Inf) < 1e-5 + @test norm(res_lm_polar["bus_results"].θ .- res_ref["bus_results"].θ, Inf) < 1e-5 + @test norm(res_lm_rect["bus_results"].Vm .- res_ref["bus_results"].Vm, Inf) < 1e-5 + @test norm(res_lm_rect["bus_results"].θ .- res_ref["bus_results"].θ, Inf) < 1e-5 + # Polar LM and rectangular LM agree with each other. + @test norm( + res_lm_polar["bus_results"].Vm .- res_lm_rect["bus_results"].Vm, Inf) < 1e-5 + @test norm( + res_lm_polar["bus_results"].θ .- res_lm_rect["bus_results"].θ, Inf) < 1e-5 +end + @testset "Rectangular CI: multi-period previous-solution warm start" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") pf = ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}(;