PaperlayerView exact revision on arXiv

Source-linked AI summary

Symbolic Direct Construction of Augmented Systems for Levenberg-Marquardt Power Flow

Tianshi Cheng, Venkata Dinavahi

arXiv:2610.34509v1eess.SYcs.NAcs.PFcs.SYmath.NA

TL;DR

Ill-conditioned power-flow equations and the sparse-product cost of conventional damped normal equations motivate this work. The paper directly constructs an augmented LM system and proves it quasi-definite for positive damping. Its reported results include faster assembly, fewer factorizations for adaptive RTE solves, and a fully accounted case2383wp run.

  • Problem

    Conventional LMPF forms damped normal equations whose J^T J product introduces fill, squares the undamped spectral condition number, and adds sparse-product cost.

    Source: §I · p. 1

  • Method

    The paper directly constructs augmented LM blocks from symbolic Jacobian structure and proves the resulting system quasi-definite for positive damping.

    Source: §I · p. 1

  • Results

    Adaptive trust-region initialization solves flat-started RTE cases in 19–25 factorizations, while case2383wp takes 237 ms for 84 iterations, 84 factorizations, and 168 additional resolves.

    Source: F. LM Variants Under Factorization… · p. 6, §VIII · p. 8

  • Takeaways & Limitations

    The authors conclude that reproducible ill-conditioned power-flow benchmarks should report initialization and complete factorization and resolve accounting.

    Source: §VIII · p. 8

Abstract

from arXiv · show

Levenberg-Marquardt power flow (LMPF) addresses ill-conditioned equations and provides least-squares stationary points for infeasible operating requests. Conventional implementations form damped normal equations, where J^T J introduces fill, squares the spectral condition number of the undamped part, and adds sparse-product cost. We propose an augmented LMPF formulation with a symbolic-permutation-based Jacobian operator that fills the augmented blocks directly. We prove that the matrix is symmetric quasi-definite for positive damping, enabling fixed-order, pivot-free QDLDL factorization. This property also motivates the solver choice: KLU value-only refactorization can reuse an unsuitable pivot order, forcing fresh full factorizations on difficult LMPF trajectories. We further introduce one-probe adaptive trust-region initialization for the radius and damping. Direct assembly is 7.5-11.9x faster than the fastest coordinate-list path and 18.0-58.8x faster than cached normal-equation sparse-product formation; QDLDL reduces numerical refactorization and triangular-solve time by up to 3.40x on identical matrices. Adaptive initialization solves flat-started RTE cases in 19-25 factorizations, 7.5-11.6x faster than reproduced three-step LM. On case2383wp, our 84-iteration run performs 84 factorizations and 168 resolves in 237 ms, versus published times of 200.780-441.72 s with fewer outer iterations, corresponding to speedups of 847-1860x. These results require complete factorization and resolve accounting and reveal initialization and legacy-parameter effects. DC-initialized NR exhibits fast quadratic convergence throughout every tested solvable loading range, motivating clearer, reproducible definitions of ill-conditioned power-flow benchmarks.

I. INTRODUCTION

The paper extends symbolic direct Jacobian construction to augmented Levenberg–Marquardt power flow, while emphasizing that initialization, solution existence, and near-singularity must be distinguished when evaluating convergence. It also frames reproducible cost comparisons around the linear-algebra work hidden by outer-iteration counts.

I. INTRODUCTION

  • Earlier power-flow work addressed ill-conditioning through optimal multipliers, Newton-like alternatives, and LM damping, trust-region, or composite-step variants –.The introduction also cites LM regularization, as a basis for those power-flow applications.Source: §I · p. 1
  • Convergence studies should distinguish poor initialization outside the NR basin, loading beyond the point of solution existence, and near-singularity while a solution remains reachable.Reliable attribution of convergence or failure requires comparisons that separate these conditions.Source: §I · p. 1
  • Outer-iteration counts omit variable trial work: one accepted LM step may require several factorizations, while normal-equation products and COO operations add implementation overhead.The paper therefore calls for complete cost accounting when comparing sparse-factorization backends.Source: §I · pp. 1–2
  • The paper extends symbolic direct Jacobian construction to the augmented LM system and combines it with a verified solver class, adaptive initialization, and complete cost accounting.The introduction lists exact normal-equation equivalence and quasi-definiteness among the formulation’s results.Source: §I · p. 1

A. Reduced Power-Flow Equations

  • The reduced polar state contains voltage angles at nonslack buses and voltage magnitudes at PQ buses; its mismatch vector contains active-power mismatches at nonslack buses and reactive-power mismatches at PQ buses.The Jacobian is defined as J = ∂F/∂x and evaluated using the same complex-domain expressions as the production NR solver.Source: A. Reduced Power-Flow Equations · p. 2

B. Least-Squares and Damped Steps

The section formulates power flow as a least-squares problem with damped Gauss–Newton steps. It explains how damping filters singular directions and how trial steps are assessed and updated.

B. Least-Squares and Damped Steps

  • LM replaces the power-flow root equation with minimization of a least-squares objective, and its damped Gauss–Newton step uses positive damping µ and a positive diagonal metric D.The experiments use D = I in retained polar coordinates; the implementation also supports voltage- and network-dependent diagonal metrics.Source: B. Least-Squares and Damped Steps · p. 2

C. Singular-Value Filtering and Step Limits

  • Damping preserves Newton-like corrections in well-observed singular directions but attenuates small-singular-value directions; at rank deficiency the correction is minimum-norm least-squares.At large damping, the direction approaches scaled steepest descent; these limits describe the direction at a fixed initial point.Source: C. Singular-Value Filtering and Step… · p. 2
  • Trial quality is measured by actual-to-predicted reduction, whose denominator is the undamped quadratic-model decrease at the damped step; Nielsen rules update damping after acceptance or rejection.A step is accepted when ρ > 10^-4; high-quality steps reduce µ and rejected steps increase it.Source: C. Singular-Value Filtering and Step… · p. 2

D. DC Initialization

The section contrasts DC initialization with the numerical and structural costs of normal equations, then introduces an equivalent augmented system. Positive damping makes that system quasi-definite, supporting fixed symbolic ordering without numerical pivot changes.

D. DC Initialization

  • DC power flow supplies an angle initial estimate with one linear solve, while voltage magnitudes stay at specified flat-start values; transformer injection corrections can materially affect initialization quality.The DC linear solve is counted as initialization, separately from subsequent NR iterations.Source: D. DC Initialization · p. 2

A. Why Avoid Explicit Normal-Equation Formation

  • Explicitly forming J^T J adds a sparse product and possible fill, squares the undamped spectral condition number, and discards the direct-Jacobian assembly advantage.These concerns are especially relevant near a fold, where a Newton-like direction requires small µ.Source: A. Why Avoid Explicit Normal-Equation… · p. 2
  • Reusing KLU’s earlier pivot order can become unsuitable as LMPF iterates change, and the reported case6515rte control stalled without a factorization failure.Fresh full factorization of the same matrices converged, motivating QDLDL alongside its measured speed difference.Source: A. Why Avoid Explicit Normal-Equation… · p. 2

B. Exact Augmented Form

  • Eliminating the auxiliary block recovers the damped normal-equation step, while the augmented representation changes the sparse linear-system form.A fixed fill-reducing symbolic ordering can be reused for all Jacobian values and positive µ; the factorization requires a positive damping floor.Source: B. Exact Augmented Form · p. 3
  • For positive damping and positive diagonal D, the augmented matrix is symmetric quasi-definite for every real square J, with inertia (n, n, 0).The proof relies on a positive leading block µD and negative trailing block −I.Source: B. Exact Augmented Form · p. 3

C. Damping-Only Updates

Damping retries update only the leading diagonal entries after the initial assembly, allowing reuse of the Jacobian blocks and sparsity pattern. This makes factorization and resolve counts observable measures of nested solver work.

C. Damping-Only Updates

  • After one assembly, changing µ updates n scalar diagonal slots while reusing J, J^T, the −I block, and the sparsity pattern.Each damping retry refactorizes the same augmented pattern and computes a candidate step.Source: C. Damping-Only Updates · p. 3
  • Composite methods can issue triangular solves for additional right-hand sides while retaining the same factorization.The passage identifies these solves alongside damping retries as directly observable implementation work.Source: C. Damping-Only Updates · p. 3
  • Factorization and triangular-resolve counts expose nested work that outer-iteration counts alone do not show.This accounting follows from retries and additional right-hand sides within an iteration.Source: C. Damping-Only Updates · p. 3

A. Admittance-Prefix Structure

Sorted Ybus columns provide prefix-based cuts for the reduced Jacobian, while symbolic segment lengths and diagonal ranks support direct address calculation during filling.

A. Admittance-Prefix Structure

  • Angle and PQ-magnitude Jacobian columns use retained active- and reactive-power row prefixes from Ybus, making the reduced pattern a type-restricted view of its graph.Symbolic analysis also supplies the local diagonal rank within each column.Source: A. Admittance-Prefix Structure · p. 3
  • Symbolic prefix sums locate reduced-Jacobian columns and their structural transpose; the augmented pattern places these graphs off-diagonal and appends two diagonal blocks.No numerical values are needed to construct the pattern.Source: A. Admittance-Prefix Structure · p. 3

B. Addresses Derived Inside the Fill Loop

  • Each fill-loop iteration derives destination addresses from Ybus starts, retained segment lengths, and diagonal ranks, replacing address searches with integer additions.In the upper-triangle layout, the diagonal correction is addressed by adding the local rank to the derived residual-column base.Source: B. Addresses Derived Inside the Fill… · pp. 3–4
  • The full-layout column path derives starts for all four Jacobian quadrants from shared block bases and local offsets, replacing persistent quadrant start arrays.Shared graph bases and segment lengths keep global layout and local writes consistent by construction.Source: B. Addresses Derived Inside the Fill… · p. 4

C. Column Fill and Direct Transpose Fill

The implementation fills augmented Jacobian blocks directly, including a transpose-oriented route for QDLDL, and uses a trust-region LM loop that retries damping values on a fixed pattern.

C. Column Fill and Direct Transpose Fill

  • Column-oriented kernels write Jacobian values directly to final positions; the QDLDL upper-triangle route computes J^T into residual columns without materializing a numerical J block.The symbolic transpose supplies row-ordered admittance entries, and phase-shifting transformers preserve the symmetric pattern despite numerically unsymmetric Ybus values.Source: C. Column Fill and Direct Transpose Fill · p. 4

D. Storage and Retry Complexity

  • The upper-triangle augmented pattern stores nnz(J) + 2n entries versus 2 nnz(J) + 2n for the full layout; numerical fill is Θ(M) without sorting or dynamic insertion.A damping retry changes n diagonal values; symbolic analysis and ordering are reused for a topology and bus-type partition.Source: D. Storage and Retry Complexity · p. 4

A. Trust-Region LM

  • Trust-region LM bounds the damped step by ∥δ∥D ≤ ∆, while damping shapes the trial step and the radius independently limits its length.Each outer iteration can try multiple damping values, refactorizing the same augmented pattern; oversized candidates trigger a doubled damping value, while evaluated rejections halve the radius and double damping.Source: A. Trust-Region LM · p. 4

B. Three-Step Composite LM

The paper compares a three-step composite LM implementation with adaptive trust-region initialization, which probes the starting point to choose whether to retire or initialize the radius.

B. Three-Step Composite LM

  • Three-step LM evaluates three corrections from the same state, Jacobian, and residual using one factorization and two additional triangular resolves.The fresh-J control instead reassembles and refactorizes at every correction.Source: B. Three-Step Composite LM · p. 4

C. Adaptive TR: One-Probe Initialization

  • Adaptive TR changes only TR initialization and uses one probe to assess merit retention before selecting damped Newton or initializing the trust-region radius and damping.A prescribed radius that is too small can slow a healthy start, while one that is too large can admit a damaging first step from a poor flat start.Source: C. Adaptive TR: One-Probe Initialization · pp. 4–5
  • Healthy starts had probe retention r = 0.003–0.14, while diverging stressed flat starts had r = 200–104; reported adaptive-TR runs used κ = 5.The initialization probe costs one factorization and one residual evaluation per solve.Source: C. Adaptive TR: One-Probe Initialization · p. 5

D. Assembly-Path Terminology

The assembly comparison measures complete LM solves and finds direct augmented assembly faster than both the fastest COO path and cached normal-equation formation. Controls compare identical iterates and backend inputs, while convergence tests distinguish initialization failure, solution existence, and conditioning.

E. Experimental Controls

  • Assembly variants execute identical iterates and share Jacobian kernels; backend comparisons use the same full augmented matrix and verify equal iteration and factorization counts.Timings are medians of seven rounds after warmup, with complete optimized-release runs reported for algorithm-level wall-clock tables.Source: E. Experimental Controls · p. 5
  • Flat and DC starts use identical network data and tolerances, while loading scans follow one fixed injection direction and continuation brackets the fold along that direction.Residual and merit distinguish solved points from nonzero-residual stationary points; Jacobian conditioning is assessed at a specified solution.Source: E. Experimental Controls · p. 5

A. Assembly and End-to-End Speedups

  • 7.5–11.9× assembly speedup over the fastest COO path and 18.0–58.8× over NE-cached accompanied end-to-end speedups of 1.15×–1.43× and 1.28×–3.57×, respectively.All paths included every damping trial and factorization; AUG-upper had the lowest end-to-end time on IEEE 39 and IEEE 118, while AUG-operator led on PEGASE 9241.Source: A. Assembly and End-to-End Speedups · p. 5
  • Construction paths produced identical accepted-step and factorization counts, and backend solutions agreed to floating-point precision; on PEGASE 9241, QDLDL and KLU differed by max |∆V| = 6.2 × 10−15.The row-oriented baseline reproduced the production column Jacobian bitwise on IEEE 39.Source: A. Assembly and End-to-End Speedups · p. 5

B. Factorization Backend on Identical Matrices

The backend comparison favors QDLDL on identical augmented matrices, while the RTE and case2383wp experiments distinguish initialization sensitivity from solution limits and report costs with factorization and resolve counts.

B. Factorization Backend on Identical Matrices

  • QDLDL is faster on every tested system, reaching a 3.40× KLU/QDLDL ratio on PEGASE 9241; value-only KLU refactorization stalled on case6515rte, requiring fresh full factorizations to converge.Fresh KLU used 44 full factorizations and 2.242 s of linear-solve time; augmented QDLDL used 0.312 s on the same accepted-step sequence.Source: B. Factorization Backend on Identical… · pp. 5–6

C. RTE Systems: Initialization Versus Existence

  • Flat-start NR failed at base loading in all three stressed RTE systems, while dc-started NR converged in 6–8 iterations and at every tested loading point before the continuation fold.Its first failure occurred only after the loading boundary was crossed; the dc angle correction was 0.05–0.08 radians.Source: C. RTE Systems: Initialization Versus… · p. 6

D. Case2383wp: Conditioning and Reported Cost

  • NR converged on case2383wp in four iterations from both starts; its estimated condition number rose from 6.0 × 10^5 at the base solution to 2.8 × 10^6 near the fold, where genuine near-singularity appears.The reported historical divergence was explained by parameter settings in the MATPOWER case definition.Source: D. Case2383wp: Conditioning and… · p. 6
  • On case2383wp, TS used 84 factorizations and 168 resolves in 237 ms, while adaptive TR used 10 factorizations in 30.3 ms.TS performs one factorization and two additional triangular resolves per outer iteration; adaptive TR required 9 iterations.Source: D. Case2383wp: Conditioning and… · p. 6
  • Our TS run was 1.86 × 10^3 and 847 times faster than published TSLM and FTSLM times, respectively, despite more outer iterations; timing ratios also include implementation and hardware effects.The published times were 441.72 s and 200.780 s, with 11 and 5 outer iterations; complete factorization, resolve, assembly, and solver measurements support the comparison.Source: D. Case2383wp: Conditioning and… · p. 6

E. Tripathy 11-Bus System: A Nonzero Merit Floor

The Tripathy 11-bus nominal request lies beyond the feasible branch, where LM returns a least-squares stationary point. Across LM variants, factorization accounting distinguishes adaptive trust-region performance from composite-method work that outer-iteration counts can obscure.

E. Tripathy 11-Bus System: A Nonzero Merit Floor

  • For the Tripathy 11-bus system, continuation places the fold at α ≈ 0.643, below nominal loading; LM reaches the same least-squares stationary point from both starts.The nonzero residual and continuation result place the nominal specification beyond the feasible branch.Source: E. Tripathy 11-Bus System: A Nonzero… · p. 6

F. LM Variants Under Factorization Accounting

  • Adaptive TR converged from RTE flat starts in 19–25 factorizations, with wall-time speedups of 7.5×, 10.9×, and 11.6× over TS.Residual-scaled composite LM required 152–245 factorizations; fresh-Jacobian corrections changed accepted-iteration counts by less than 2%.Source: F. LM Variants Under Factorization… · p. 6
  • A single-step residual-scaled variant reported 15 outer iterations but executed 1,214 factorizations through damping retries, showing why outer-loop counts alone can misrank LM work.A composite outer step can evaluate three corrections using one shared factorization and additional triangular resolves.Source: F. LM Variants Under Factorization… · pp. 6–7

VII. DISCUSSION

The discussion separates initialization failure from solution nonexistence and near-singularity, and identifies LM’s role in difficult starts and infeasible requests. The conclusion reports faster assembly and factorization, while emphasizing complete work accounting and reproducible benchmark definitions.

A. A Reproducible Meaning of Ill-Conditioned Power Flow

  • Flat-start failure records the outcome from one initial point, and controlled dc initialization or corrected case parameters can change it.A reproducible benchmark should specify the network version and parameters, loading direction, initialization, convergence tolerance, dc+NR control, continuation boundary, and factorization and resolve counts.Source: A. A Reproducible Meaning of… · pp. 7–8

B. The Appropriate Role of LM Power Flow

  • Trust-region LM is useful for difficult flat starts, while NR remains simpler whenever dc+NR converges; beyond the feasible region, LM returns a least-squares stationary point.Continuation locates the loading boundary, and the LM residual describes the outcome beyond it.Source: B. The Appropriate Role of LM Power Flow · p. 8

VIII. CONCLUSION

  • AUG-operator achieves 7.5–11.9× assembly speedups over the fastest COO control and 18.0–58.8× over cached normal-equation sparse-product formation; QDLDL reduces numerical refactorization plus triangular-solve time by up to 3.40× on identical augmented systems.Source: §VIII · p. 8
  • 84 TS iterations with 84 factorizations and 168 additional resolves take 237 ms on case2383wp; adaptive flat-started RTE solves use 19–25 factorizations and run 7.5–11.6× faster than reproduced TS.Outer-iteration counts alone give an incomplete measure of implementation cost.Source: §VIII · p. 8
  • DC+NR converges at every tested RTE loading point before the continuation fold and solves case2383wp in four iterations.Benchmarks should report case parameters, initialization, loading, tolerances, factorizations, and resolves.Source: §VIII · p. 8

Related papers

Loading 2610.34509v1…