LM performance improvements#338
Merged
Merged
Conversation
Contributor
There was a problem hiding this comment.
Pull request overview
This PR improves the performance characteristics of the Levenberg–Marquardt (LM) AC power flow solver (Issue #336) by reducing per-iteration sparse matrix allocations and moving toward reusable linear algebra workspace.
Changes:
- Introduces an
LMWorkspaceto build the augmented sparse matrix once and update Jacobian/λ entries in-place. - Simplifies the LM step computation by removing higher-order terms and updating damping using a carried-over parameter (
μ). - Minor formatting cleanup in
gradient_descent_ac_power_flow.jl.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
src/levenberg-marquardt.jl |
Adds LM workspace + in-place augmented matrix updates, changes damping/step logic, uses SPQR QR solve. |
src/gradient_descent_ac_power_flow.jl |
Removes an extra blank line. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| ) | ||
| actual_reduction = (residualSize - newResidualSize) | ||
| predicted_reduction = residualSize - dot(temp_x, temp_x) | ||
| actual_reduction = residualSize - newResidualSize |
There was a problem hiding this comment.
ρ = actual_reduction / predicted_reduction can divide by zero (or a negative/very small predicted_reduction), producing Inf/NaN and potentially accepting a bad step. Add a guard (e.g., if predicted_reduction <= 0 or not finite, treat as a rejected step and increase damping) before computing ρ.
Suggested change
| actual_reduction = residualSize - newResidualSize | |
| actual_reduction = residualSize - newResidualSize | |
| # Guard against zero/negative or non-finite predicted reduction to | |
| # avoid Inf/NaN ρ and accidentally accepting a bad step. | |
| if predicted_reduction <= 0.0 || !isfinite(predicted_reduction) | |
| # Treat as a rejected step: restore residual to match x. | |
| residual(x, time_step) | |
| return 0.0 | |
| end |
Comment on lines
125
to
+127
| @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) | ||
| λ, μ = update_damping_factor!(x, residual, J, μ, time_step, ws) |
There was a problem hiding this comment.
The LM loop only checks !isnan(λ), but with the new μ *= 4 updates, μ/λ can overflow to Inf (or become non-finite via factorization), and the loop will keep running until maxIterations. Consider switching the loop condition and convergence logic to isfinite(λ) (and/or bounding μ), and optionally restoring an explicit “give up after N rejected steps” safeguard similar to the previous maxTestλs behavior.
Contributor
Performance ResultsPrecompile Time
Solve Time
|
jd-lara
approved these changes
May 16, 2026
jd-lara
added a commit
that referenced
this pull request
May 16, 2026
… models vs. solvers (#371) * simplify hessian matrix structure initialization * add test of RH w/ headroom slack * address copilot's comments Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com> * simplify HomotopyHessian constructor; switch to IS.@assert_op Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com> * restore defensive fill-with-ones in HomotopyHessian constructor Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com> * Update org/domain/workspace references from NREL-Sienna to Sienna-Platform (#358) * Initial plan * docs: update org, docs host, and slack workspace references Agent-Logs-Url: https://github.com/Sienna-Platform/PowerFlows.jl/sessions/6568f54e-0a28-4f56-b017-353e8b212eb1 Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> * docs: fix README typo in updated contribution line Agent-Logs-Url: https://github.com/Sienna-Platform/PowerFlows.jl/sessions/6568f54e-0a28-4f56-b017-353e8b212eb1 Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> --------- Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> * reject LCC systems in RobustHomotopyPowerFlow LCCs add state variables to the Jacobian that the homotopy Hessian formulation does not account for, producing a dimension mismatch deep inside the solve. Fail fast at HomotopyHessian construction with a clear error message. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * add storage to possibly be PV (#359) * fix: update fields for tap transformers and v35 exporter default values (#360) * Bump version from 0.16.3 to 0.16.4 * Add literate.jl docs infrastructure (#356) * Add precommit format yaml * Add docs deps * Add docs preview cleanup * Add make_tutorials infrastructure * add generated files to gitignore * Translate tutorial to literate and fix type name and DC reporting bugs * Add prettyurls * Tidy dc power flow docstrings * Move LCC to explanation * Remove double execution * Compute url for tagged versions * Filter tutorial files with endswith * Add links as absolute url * Single image insert * Use non-greedy capture * Switch to abspath * Add input/output documentation * Add kwargs to base.show override * Clarify image regex * correct last org refs * Add local vs github link logic * move code for reuse * udpate LCC indexing * add rectangular residual calculator * add rectangular jacobian * export new type * add NR power flow for rectangular * add the setup code * use the new klu branch * exports and definitions * AI assisted testing * add code to catch unsopported methods * add the testing and docstring * bump psy version * apply the fixes from copilot's comment * delete plans * fix copilot suggestion * address comments * update static arrays * apply memory reduction suggestions * fix format and docs * test allocation regressions * delete plan * fix jacobians * add more testing * address PR comments and testing * fix the ZIP discrepancy * last LCC updates * Fix LCC Jacobian: correct dQ/dV, dQ/dt; add missing inverter-side entries The simplified-φ Q-derivatives in lcc_utils.jl had algebra errors (spurious V factor on the first term, sin²(φ) instead of V·sin(φ) in the second, missing sign(I_dc)). The Jacobian assembly was also missing all six bus-side entries on the inverter terminal — only the rectifier mirror was wired up — so finite-difference verification on any LCC system fails by construction. Add the missing entries (and their sparsity pattern), guard the helpers against the sin(φ)→0 clamp boundary (true derivative is 0 there), and log a warning when the arccos argument leaves [-1,1]. Re-enables verify_jacobian on the LCC test system, which was previously commented out for exactly these reasons. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Address Copilot review: clean warning string, clarify test comment Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix LCC docs inconsistency; global const for sin \phi tol * Deduplicate LCC tail residual/Jacobian; asymptotic Jacobian verification Factor the LCC tail residual rows and the LCC Jacobian tail-block scalar coefficients into shared helpers in lcc_utils.jl, consumed by both the polar and rectangular CI formulations. Mirrors the dispatch pattern of _update_ybus_lcc!. Drops the implicit allocations from polar's broadcast slicing of the LCC tail and prevents the two formulations from drifting on future LCC fixes. Stale docstring in the rect Jacobian about the inverter-side bus diagonal block (corrected upstream in be9ac82) is removed and replaced with a forward-looking note about porting the true-φ Q-derivative helpers. Replace the fixed-tolerance FD-parity Jacobian tests (rect LCC, four rect Jacobian testsets, polar verify_jacobian) with an asymptotic verifier in test/test_utils/jacobian_verification.jl that checks the linearization remainder decays as O(Δx²) under a geometric Δx-sweep (matching the methodology already used by test_homotopy_hessian.jl). Filters Δx-points that fall to the FP-noise floor and, on failure, prints per-Δx remainders, observed Taylor orders, and a worst-row symbolic-vs-FD diff for the suspect Jacobian entry. Each verification point now applies a deterministic random perturbation so we no longer evaluate at flat-start (θ=0) or polar-converged states where sin-style cross-terms silently vanish. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Polar LCC P-row Jacobian: migrate to true-φ via dP/dV, dP/dt helpers Polar's bus-P-row entries (∂P/∂V_fb, ∂P/∂t_fb, and the matching inverter-side entries) and the tail-row × bus-column entries (∂F_t/∂V) were using the α-approximation — treating φ as a constant equal to α and dropping the ∂φ/∂(V, t) chain terms. That's exact only when the LCC commutation reactance x_t is zero; for x_t > 0 it leaves out a real term. Add _calculate_dP_dV_lcc and _calculate_dP_dt_lcc parallel to the existing _calculate_dQ_d{V,t,α}_lcc helpers. Unlike the Q helpers, neither requires a sin(φ)→0 guard: differentiating cos(φ) introduces a -sin(φ) factor that cancels the 1/sin(φ) from ∂φ/∂x, leaving a regular closed form. ∂P/∂α is left on the α-approximation since the analogous cancellation makes that formula already exact. No observable behavior change in the current test suite — every LCC test fixture (case5_2_lcc.raw, the synthetic systems in test_jacobian.jl) has x_t = 0, which collapses the new term to zero. The change matters for the first system we hit with nonzero commutation reactance. Polar's tail × tail block (d_Ft_* fields of _lcc_jacobian_scalars) is still on the α-approximation: it's shared with rect via the NamedTuple, and migrating it would either touch the rect path or require splitting the helper. Left for a follow-up. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * True-φ LCC Jacobian: polar P-row clamp guard + rect migration The polar P-row entries had a latent clamp boundary bug: the α-approximation form (cos(α), sin(α)) is *algebraically equivalent* to true-φ in the interior — the ϕ-definition identity cos(α) - cos(ϕ) = x_t·I_dc/(√2·V·tap) makes the chain correction collapse into the α-trig form — but the identity *breaks at the clamp*, where ϕ is pinned and ∂ϕ/∂x = 0. The α-approximation form keeps using cos(α), sin(α) regardless, mismatching the residual (which sees cos(ϕ_clamp), sin(ϕ_clamp) = 0). Add _calculate_dP_d{V,t,α}_lcc helpers parallel to the existing _calculate_dQ_d* set; each guards sin(ϕ) → 0 so the chain term is dropped at the clamp (where ∂ϕ/∂x = 0 physically). The dP/dV and dP/dt helpers happen to be regular at the clamp (no 1/sin(ϕ)), but the guard still matters because the α-approximation form they replace uses cos(α), not cos(ϕ_clamp). The dP/dα helper has the same 1/sin(ϕ) singularity the Q helpers do. Rewire _lcc_jacobian_scalars to call these helpers for the d_Ft_* tail × tail entries — a correctness fix for both polar and rect tail × tail (the tail equations and tail columns are formulation-invariant). New test in test_jacobian.jl constructs a 3-bus system with large inverter x_t and small extinction angle to force ϕ_i = π (clamp); without the migration the asymptotic verifier reports observed Taylor order ≈ 1 on the inverter-α and tap_i columns. Then migrate rect's LCC Jacobian: rect's α-approximation is *not* algebraically equivalent to true-ϕ even in the interior (no analogous identity, because the substitution cos(α) → cos(ϕ) doesn't factor through (e, f) the way it does through V). Add _dphi_d{V,t,α}_lcc helpers exposing ∂ϕ/∂x with the same sin(ϕ) → 0 guard, and rewrite _set_entries_for_lcc_rect! to use the universal admittance form F_real = -A·u(ϕ), F_imag = -A·w(ϕ) (with the inverter sign convention absorbed into cos(ϕ_i), sin(ϕ_i)) plus the per-derivative chain terms. Two new rect tests cover nonzero xc in the interior and the inverter ϕ-clamp; both pass. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Replace stale α-approx comment in polar; expose dP_dt_* in scalars The "True-φ ∂P/∂V and ∂P/∂t for the rectifier side..." comment block in `_set_entries_for_lcc` was written in the first polar P-side migration commit and is no longer accurate: - "polar α-approximation drops the ∂φ/∂Vm chain term — which vanishes only when x_t = 0": false. The α-form is *algebraically equivalent* to true-ϕ in the interior via the ϕ-definition identity `cos(α) − cos(ϕ) = x_t·I_dc/(√2·V·tap)`; the chain term doesn't vanish, it cancels into the α-trig form. - "For ∂P/∂α the analogous cancellation makes the α-approximation exact, so those entries stay on `common_alpha_*`": also false now. The α-form is exact only in the interior; at the clamp it's wrong (the inverter-ϕ-clamp test caught it), so we migrated ∂P/∂α to the guarded `s.dP_dα_*` in the rect-migration commit. Replace with a concise comment pointing at the true-ϕ equivalence-in- the-interior and the clamp-guard rationale, and reference the test that exercises the clamp. While here: expose `dP_dt_fb`, `dP_dt_tb` from `_lcc_jacobian_scalars` alongside the existing `dP_dV_*` and `dP_dα_*`, and have polar read them all from the scalars NamedTuple instead of recomputing the same helper calls locally. No behavior change. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * rename * split abstraction between polar and rectangular * testing and performance changes * PR feedback: PV-bus Q-derivatives, noise-floor docs, cleanups Address Copilot review on PR #366: * **Bug: PV-bus LCC Q-derivatives missing** (`src/ac_power_flow_jacobian.jl`): ∂Q/∂tap and ∂Q/∂α at the LCC AC terminal were inside the `bus_type == PQ` block. For PV terminals the Q-mismatch row still depends on tap/α through Q_lcc, so those entries should be filled regardless of bus type — the V-dependent entries (∂{P,Q}/∂V and the tail-row × V chain rule) are the only ones that need PQ gating. Move bus-row × tail-column writes out of the PQ block; gate only the V-column writes. Add a "Jacobian verification with LCC at a PV terminal" test that exercises this regime. * **Noise-floor sensitivity audit** (`test/test_utils/jacobian_verification.jl`): document the analytical detection threshold for entry-level Jacobian errors: `|δ|_min ≈ c · eps · ||F0|| / Δx_max ≈ 2.2e-9 · max(||F0||,1)` with current defaults. Print the sensitivity floor in the failure diagnostic so a reader can tell what a passing test does *not* exclude. Tighten `Jv` type to `AbstractMatrix` (Copilot noted the failure path uses `Jv[row, j]` which requires `getindex`). * **Cleanups**: replace `sqrt(6)/π` with the existing `SQRT6_DIV_PI` global const in `lcc_utils.jl` (3 occurrences). Trim the long-winded comment block in polar `_set_entries_for_lcc` introduced earlier this PR. Replace stale "Both thyristor angles must be strictly positive" comment in the polar LCC verify_jacobian test with one that accurately points at the separate inverter-ϕ-clamp test for boundary coverage. Deferred to a follow-up commit: setpoint_at_rectifier=false branching (needs sparsity-structure additions on both polar and rect). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix warm start for rectangular * format change * PR comments * add more testing * wrap test in logging * Bump version from 0.16.4 to 0.17.0 * Accept narrowly-typed Dicts for ACPowerFlow solver_settings (#369) Broaden the kwarg type to AbstractDict and coerce internally to Dict{Symbol, Any}, so that callers can pass a plain literal like `Dict(:maxIterations => 50)` (inferred as Dict{Symbol, Int}) without having to spell out the value type. Dict is invariant in its parameters, so the prior `Dict{Symbol, Any}` annotation rejected these. Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix broken test * LM performance improvements (#338) * LM performance improvements * copilot code review; fix x0 kwarg * add total iterations to RH method * add method comparison script * review the PR and add rectangular testing --------- Co-authored-by: Jose Daniel Lara <jdlara@berkeley.edu> * Bring branch up to date with main; restructure power-flow docs around models vs. solvers Agent-Logs-Url: https://github.com/Sienna-Platform/PowerFlows.jl/sessions/e631d3fd-287d-4d0a-92c1-1bfbfb3a6c01 Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> --------- Co-authored-by: Luke Kiernan <lkiernan@nrel.gov> Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com> Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com> Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> Co-authored-by: Rodrigo Henríquez-Auba <rodrigomha@gmail.com> Co-authored-by: Marck Llerena V. <140716266+mcllerena@users.noreply.github.com> Co-authored-by: kdayday <kdayday@users.noreply.github.com> Co-authored-by: Jose Daniel Lara <jdlara@berkeley.edu> Co-authored-by: Luke Kiernan <86331877+luke-kiernan@users.noreply.github.com> Co-authored-by: Copilot <copilot@github.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fix up LM, in collaboration with @jmaack24. Address issue #336.
λfrom previous iterationAandb