Literal PMDSE: explicit matrix-based explicit-neutral DSSE (WLS / WLAV / MLE)#2
Open
MohamedNumair wants to merge 6 commits into
Open
Literal PMDSE: explicit matrix-based explicit-neutral DSSE (WLS / WLAV / MLE)#2MohamedNumair wants to merge 6 commits into
MohamedNumair wants to merge 6 commits into
Conversation
Textbook, JuMP-free state estimator that consumes the same data_math + meas dictionaries as the IVR explicit-neutral estimator (IVRENPowerModel) but solves the DSSE with explicit H/W/Gain matrices and Gauss-Newton instead of JuMP. - literal_core.jl: (bus,terminal) index map, nodal Ybus assembly (series + Pi-shunts + bus shunts), configurable references (:sota/:full_slack/:prop), flat start, solution dict. - measurement_model.jl: rectangular EN h(x) closures for every IVR/IVREN measurement var-type; component->nodal aggregation with neutral = -sum(phases) and zero-injection pseudo-measurements for observability. - solve_wls.jl: Gauss-Newton WLS, normal equations with QR/orthogonal fallback. - solve_wlav.jl: IRLS WLAV (robust to gross errors). - mle.jl: general log-likelihood interface (Newton); Gaussian reduces to WLS. - literal_pmdse.jl: solve_mc_se_literal wrapper, LiteralResult, accuracy helpers. Exports solve_mc_se_literal and registers the bare/ includes. https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
- test_ybus.jl (gate 1): Ybus*U reproduces PMD nodal current balance
- test_hx_jacobian.jl(gate 2): AD H vs finite differences, every measurement type
- test_references.jl (gate 3): rotation unobservable w/o angle datum, observable
with it (reproduces the main.tex observability proof)
- test_wls_vs_ivren.jl(gate 4): noiseless literal WLS == IVREN/PF state (~1e-16)
- test_benchmark.jl (gate 5): IVREN vs literal WLS/WLAV accuracy + time/iters
Wired into runtests.jl via SafeTestsets. docs/src/Literal_PMDSE_PLAN.md documents
the design, modelling decisions, sign conventions and the benchmark table.
https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
Printf is not a declared package/test dependency; replace @printf table formatting with round()+string interpolation so the suite stays green under Pkg.test(). Statistics (a package dependency) is retained. https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
There was a problem hiding this comment.
Pull request overview
This PR adds a new “Literal PMDSE” explicit-neutral (4-wire) distribution-system state estimation path that is matrix-based and JuMP-free, intended to mirror the existing IVREN formulation while solving WLS/WLAV/MLE via explicit measurement model construction and ForwardDiff Jacobians.
Changes:
- Adds a new literal estimator stack (
LiteralModel, measurement model assembly, WLS/WLAV solvers, and a general MLE interface) and exports the new public entrypoints. - Introduces a dedicated
test/literal/suite (5 “validation gates”) and wires it intotest/runtests.jlviaSafeTestsets. - Adds a planning/technical design document under
docs/src/.
Reviewed changes
Copilot reviewed 16 out of 16 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
| test/runtests.jl | Adds SafeTestsets and runs the new literal estimator gate tests. |
| test/literal/test_ybus.jl | Gate 1: validates Ybus assembly against PMD power-flow current balance. |
| test/literal/test_hx_jacobian.jl | Gate 2: validates AD Jacobian vs finite differences across measurement types. |
| test/literal/test_references.jl | Gate 3: validates observability behavior under reference schemes. |
| test/literal/test_wls_vs_ivren.jl | Gate 4: compares literal WLS vs PF/IVREN solutions (noiseless + noisy). |
| test/literal/test_benchmark.jl | Gate 5: benchmark comparison of IVREN vs literal WLS/WLAV. |
| test/literal/helpers.jl | Shared helper utilities for toy feeder, PF, measurements, and comparisons. |
| src/PowerModelsDistributionStateEstimation.jl | Includes the new literal estimator source files into the module. |
| src/core/export.jl | Exports the literal estimator APIs and helper utilities. |
| src/bare/literal_core.jl | Implements node indexing, Ybus assembly, reference handling, flat start, and solution materialization. |
| src/bare/measurement_model.jl | Implements h(x) row closures and model assembly (SEModel) for supported measurement types. |
| src/bare/solve_wls.jl | Implements Gauss–Newton WLS and defines LiteralResult. |
| src/bare/solve_wlav.jl | Implements IRLS-based WLAV. |
| src/bare/mle.jl | Implements a general MLE solver (damped Newton) and Gaussian default mapping. |
| src/bare/literal_pmdse.jl | Adds solve_mc_se_literal wrapper plus accuracy/voltage error helpers. |
| docs/src/Literal_PMDSE_PLAN.md | Design/plan document describing formulation, assumptions, and validation gates. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+29
to
+39
| function _solve_gain(G::Matrix{Float64}, rhs::Vector{Float64}, | ||
| A::Matrix{Float64}, b::Vector{Float64}; cond_tol::Float64 = 1e12) | ||
| if all(isfinite, G) && all(isfinite, rhs) | ||
| ch = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(G); check = false) | ||
| if LinearAlgebra.issuccess(ch) | ||
| Δ = ch \ rhs | ||
| all(isfinite, Δ) && return Δ, false | ||
| end | ||
| end | ||
| return A \ b, true # QR least squares on √W·H (orthogonal estimator) | ||
| end |
Comment on lines
+281
to
+283
|
|
||
| return SEModel(lm, rowfns, z, w, rowmeta, Any[]) | ||
| end |
Comment on lines
+288
to
+291
| x = expand_state(lm, x_free) | ||
| vr, vi = vrvi(lm, x) | ||
| Ir, Ii = nodal_injection(lm, collect(vr), collect(vi)) | ||
| ctx = HCtx{T}(collect(vr), collect(vi), Ir, Ii) |
Comment on lines
+22
to
+25
| "log-likelihood of one row given the model prediction `ĥ` (Gaussian / general)" | ||
| loglik_row(dst::_DST.Distribution, ĥ) = _DST.logpdf(dst, ĥ) | ||
| loglik_row(dst::Real, ĥ) = -((ĥ - dst)^2) # degenerate hard value | ||
|
|
Comment on lines
+10
to
+13
| iterations::Int | ||
| objective::Float64 # weighted SSR r' W r at the solution | ||
| solve_time::Float64 | ||
| gain_cond::Float64 |
Comment on lines
+62
to
+66
| terms = haskey(refb, "terminals") ? refb["terminals"] : 1:length(bsol["vr"]) | ||
| for idx in 1:length(bsol["vr"]) | ||
| !include_neutral && idx == _N_IDX && continue | ||
| push!(errs, abs((bsol["vr"][idx] + im * bsol["vi"][idx]) - | ||
| (refb["vr"][idx] + im * refb["vi"][idx]))) |
Replace the full-Newton/Levenberg negloglik minimisation with Fisher scoring: each step solves (Hᵀ Λ H) Δ = Hᵀ g with per-row score g_i = ∂logpdf/∂h_i and information Λ_i = -∂²logpdf/∂h_i² (via ForwardDiff), reusing the WLS QR/orthogonal fallback. For Gaussian dst this is exactly the WLS normal-equation step, so MLE now converges in ~2 iterations to the WLS estimate even on the ill-conditioned (:sota) gain where the previous Newton path hit maxiter. Adds an MLE==WLS test. https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
- docs/src/literal_pmdse.md: full usage guide for solve_mc_se_literal (every argument, the :wls/:wlav/:mle estimators, the :sota/:full_slack/:prop reference schemes, LiteralResult fields, accuracy helpers, lower-level API). - docs/src/literal_pmdse_observability.md: how Literal PMDSE resolves the angular & zero-sequence reference problem of main.tex (Theorem 1 -> :prop rank deficiency; the SOTA fix -> :sota), with the rotation null-vector regression test and the toy 3-bus 4-wire example. - examples/literal_pmdse_example.jl: end-to-end worked example (estimators, references, result inspection, robustness, lower-level API). - docs/make.jl: add both pages to the manual navigation. https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
- Point @ref links at the actual page headers ("Measurement Conversion") and drop the unresolved IVRENPowerModel symbol ref / fragile em-dash slug. - Describe :mle as Fisher scoring (matches the implementation). https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
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.
Literal PMDSE
A textbook, JuMP-free distribution-system state estimator that consumes the same
data_math+data_math["meas"]+data_math["se_settings"]dictionaries as the JuMP IVR explicit-neutral estimator (solve_ivr_en_mc_se/IVRENPowerModel), but solves the DSSE by forming the measurement modelz = h(x) + e, the JacobianH = ∂h/∂x(viaForwardDiff) and iterating the normal equations(HᵀWH)Δx = HᵀW r(Abur & Expósito).Design
src/bare/literal_core.jl—(bus,terminal)index map, nodalYbusassembly (mutually-coupled series + Π line-shunts + bus shunts), configurable references, flat start.src/bare/measurement_model.jl— rectangular explicit-neutralh(x)closures for every IVR/IVREN measurement var-type; component→nodal aggregation with a wye neutral injection= −Σ(phase injections)and zero-injection pseudo-measurements for observability.src/bare/solve_wls.jl— Gauss–Newton WLS, normal equations with a QR / orthogonal fallback when the gain is ill-conditioned.src/bare/solve_wlav.jl— IRLS WLAV (robust to gross errors).src/bare/mle.jl— general log-likelihood interface (damped Newton); Gaussian reduces to WLS (verified ~1e-16).src/bare/literal_pmdse.jl—solve_mc_se_literalwrapper +LiteralResult+ accuracy helpers.References / observability (reproduces
main.tex)bus_type == 3is the reference bus. Schemes::sota(ground neutral + one angle datumvi[ref,phase₁]=0— full-rank without forcing a balanced reference bus),:full_slack(fix the whole reference-bus phasor),:prop(neutral only — rank-deficient by one). The rotationδx=[−vi;vr]is a null vector ofHunder:propand becomes observable under:sota.Validation gates (all pass)
Ybus·Ureproduces PMD's nodal current balance~1e-10Hvs finite differences, every measurement type<1e-6‖H·δx_rot‖ ~1e-17vs5e-2~1e-16(:full_slack&:sota);<1e-4vs JuMP IVRENWired into
test/runtests.jlviaSafeTestsets(test/literal/).Benchmark — 3-bus 4-wire, σ = 0.01, 4 noise seeds
Literal Gauss–Newton WLS converges in 2 iterations / ~1 ms, ≈500× faster than the JuMP/Ipopt IVREN solve, with comparable accuracy. Results are compared, not asserted equal.
Sign / convention decisions (verified against IVREN)
Ybus·U = Σ crg − Σ crd(gen+, load−).p_c = s·(Ir_c·Δvr_c + Ii_c·Δvi_c),q_c = s·(−Ii_c·Δvr_c + Ir_c·Δvi_c),Δphase-to-neutral,s = +1gen /−1load.I_fr = (Y_series+Y_sh_fr)·U_f − Y_series·U_t.Assumptions / limitations (flagged)
Ybus(the EN test feeders use a stiff voltage source); branches + shunts are supported.calc_admittance_matrixdoes not support the 4-wire (kron_reduce=false) case, so gate 1 validatesYbusagainst PMD's power-flow physics.See
docs/src/Literal_PMDSE_PLAN.mdfor full details.https://claude.ai/code/session_01WqnbPpV1xjaYnLarHUFWrF
Generated by Claude Code