Literal PMDSE: explicit matrix-based explicit-neutral DSSE (WLS / WLAV / MLE)

docs/make.jl

@@ -14,6 +14,8 @@ makedocs(
                             "State Estimation Criteria"               => "se_criteria.md",
                             "Angular Reference Models"                => "angular_ref.md",
                             "Explicit Neutral Models"                 => "explicit_neutral_models.md",
+                            "Literal PMDSE (matrix-based)"            => "literal_pmdse.md",
+                            "Literal PMDSE: Reference & Observability" => "literal_pmdse_observability.md",
                             "Bad Data Detection and Identification"   => "bad_data.md",
                             ],
              "Library " => [

docs/src/Literal_PMDSE_PLAN.md

@@ -0,0 +1,130 @@
+# Literal PMDSE — explicit matrix-based explicit-neutral DSSE
+
+**Status: complete — all 5 validation gates pass.**
+
+`Literal PMDSE` is a textbook, JuMP-free distribution-system state estimator that
+consumes the *same* PMD mathematical dictionaries (`data_math`,
+`data_math["meas"]`, `data_math["se_settings"]`) as the JuMP IVR explicit-neutral
+estimator (`solve_ivr_en_mc_se` / `IVRENPowerModel`, `src/prob/se_en.jl`) but
+solves the DSSE by forming the measurement model `z = h(x) + e`, the Jacobian
+`H = ∂h/∂x` (via `ForwardDiff`) and iterating the normal equations
+`(HᵀWH)Δx = HᵀW r` (Abur & Expósito, *Power System State Estimation*).
+
+```julia
+res = solve_mc_se_literal(data_math; estimator = :wls,  reference = :sota)
+res = solve_mc_se_literal(data_math; estimator = :wlav, reference = :full_slack, ref_values = …)
+res = solve_mc_se_literal(data_math; estimator = :mle)        # Gaussian MLE == WLS
+```
+
+## Files
+
+```
+src/bare/
+  literal_core.jl       # node index map, Ybus assembly, references, flat start, solution dict
+  measurement_model.jl  # rectangular EN h(x) closures for every measurement var-type
+  solve_wls.jl          # Gauss–Newton WLS (normal equations + 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
+test/literal/
+  helpers.jl            # toy network, toy PF, feeder loaders, measurement builders
+  test_ybus.jl          # gate 1
+  test_hx_jacobian.jl   # gate 2
+  test_references.jl    # gate 3
+  test_wls_vs_ivren.jl  # gate 4
+  test_benchmark.jl     # gate 5
+```
+
+## Formulation (nodal, rectangular, explicit-neutral)
+
+* **State** `x = [vr; vi]` over every `(bus, terminal)` pair, terminals including
+  the neutral `_N_IDX = 4`.  Branch / load / gen currents are eliminated.
+* **Admittance** `Ybus = G + jB`: each branch contributes its full mutually
+  coupled series admittance `Y = (br_r + j·br_x)^{-1}` stamped over
+  `f_connections`/`t_connections`, plus its Π line-shunts (`g_fr,b_fr,g_to,b_to`)
+  and any bus shunts.  Validated to reproduce PMD's power-flow current balance
+  (`Ybus·U == Σgen − Σload` at every node).
+* **Nodal injection** `I = Ybus·U`, evaluated `Ir = G·vr − B·vi`,
+  `Ii = B·vr + G·vi` (ForwardDiff-friendly).
+* **Measurement functions** (EN, phase-to-neutral where relevant):
+  `vr,vi` identity; `vm/vmn = |U_c−U_n|`; `va = ∠(U_c−U_n)`; `vll = |U_i−U_j|`;
+  branch `cr,ci,cm,ca` from `Y_ff·U_f + Y_ft·U_t`; branch power `p,q` (EN);
+  current injections `crd,cid,crg,cig` and powers `pd,qd,pg,qg,ptot,qtot` and
+  `cmd,cmg,cad,cag` via the nodal injection `Ybus·U`.
+
+## Key modelling decisions
+
+1. **Component → nodal aggregation.**  Current/power *injection* measurements are
+   mapped to the nodal injection `Ybus·U`.  When several components share a
+   bus terminal (e.g. the 4 loads on one bus of `3bus_4wire`), their measured
+   currents are summed into one nodal residual — exact at the solution and the
+   only consistent choice once currents are eliminated.  A wye component's
+   **neutral injection is `−Σ(phase injections)`**.
+2. **Zero-injection pseudo-measurements** (`z = 0`, high weight) close the model
+   at source-free terminals, reproducing the KCL equalities of `IVRENPowerModel`
+   and making intermediate buses observable.
+3. **References (configurable).**  `bus_type == 3` is the reference bus.
+   * `:sota` (default) — ground the reference-bus neutral **and** fix one angle
+     datum `vi[ref, phase₁] = 0`.  This is the SOTA combination from `main.tex`
+     that makes the gain full-rank without forcing a balanced reference bus.
+   * `:full_slack` — fix the whole reference-bus phasor (used for the apples-to-
+     apples gate-4 comparison against IVREN).
+   * `:prop` — neutral only; **rank-deficient by one** (the observability trap).
+   Grounded terminals (`bus["grounded"]`) are always fixed to 0.
+4. **Numerics.**  Normal equations with a Cholesky solve, falling back to the
+   orthogonal QR least-squares on `√W·H` when the gain is ill-conditioned (the
+   per-unit impedances of `3bus_4wire` give `cond(G) ~ 1e17`, yet the QR fallback
+   recovers the state to machine precision).
+5. **Weights.**  `W = diag(1/(rescaler·σ)²)` with `σ = std(dst)`; `Float64`
+   `dst` entries are treated as near-exact (floored σ); a σ-floor prevents an
+   infinite weight when a measured value is exactly 0 (e.g. `vi` at a balanced
+   reference bus).
+6. **MLE-ready.**  `solve_mle` maximises `Σ logpdf(dst_i, h_i(x))` by a damped
+   Newton iteration; for Gaussian `dst` it is identical to WLS (verified to
+   ~1e-16).  An arbitrary `Distributions`/`Polynomials`/`ExtendedBeta` can be
+   supplied with no change to `h(x)` or the references.
+
+## Validation gates (all pass)
+
+| # | Gate | Result |
+|---|------|--------|
+| 1 | `Ybus·U` reproduces PMD's nodal current balance | max err `~1e-10` |
+| 2 | AD `H` vs finite differences, every measurement type | rel err `<1e-6` |
+| 3 | rotation unobservable w/o angle datum (`:prop`), observable with it (`:sota`) | `‖H·δx_rot‖ ~1e-17` vs `5e-2` |
+| 4 | noiseless literal WLS == IVREN / PF state | max\|U\| `~1e-16` (`:full_slack` & `:sota`); `<1e-4` vs JuMP IVREN |
+| 5 | benchmark IVREN vs literal WLS/WLAV (accuracy + time/iters) | see table |
+
+### Benchmark (3-bus 4-wire, σ = 0.01, 4 noise seeds)
+
+| seed | IVREN max\|U\| / t[s] | WLS max\|U\| / t[s] / it | WLAV max\|U\| / t[s] / it |
+|------|----------------------|--------------------------|---------------------------|
+| 1 | 6.9e-4 / 0.52 | 2.8e-4 / 0.001 / 2 | 2.5e-4 / 1.45 / 7 |
+| 2 | 8.6e-4 / 0.45 | 7.6e-4 / 0.001 / 2 | 6.8e-4 / 0.001 / 8 |
+| 3 | 3.6e-4 / 0.45 | 1.4e-3 / 0.001 / 2 | 1.3e-3 / 0.001 / 6 |
+| 11 | 4.4e-4 / 0.49 | 2.1e-4 / 0.001 / 2 | 2.0e-4 / 0.001 / 6 |
+| **mean** | **5.9e-4** | **6.6e-4** | **6.1e-4** |
+
+The literal Gauss–Newton WLS converges in **2 iterations / ~1 ms**, roughly
+**500× faster** than the JuMP/Ipopt IVREN solve, with comparable accuracy;
+WLAV is marginally more accurate (robustness).  Results are *compared*, not
+asserted equal.
+
+## Sign / convention notes (verified against IVREN)
+
+* Nodal injection `Ybus·U = Σ crg − Σ crd` (gen `+`, load `−`).
+* Power injection `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 = +1` gen / `−1` load.
+* Branch from-current `I_fr = (Y_series+Y_sh_fr)·U_f − Y_series·U_t`; branch flow
+  `p = cr·Δvr_f + ci·Δvi_f` (EN).
+
+## Limitations / future work
+
+* Transformers are not yet folded into `Ybus` (the explicit-neutral test feeders
+  use a stiff voltage source, no in-service transformer); branches + shunts are
+  supported.
+* Power-*only* injection buses rely on the power residuals for observability of
+  that node; the canonical IVREN set uses current injections.
+* `calc_admittance_matrix` in PMD does not support the 4-wire (kron_reduce=false)
+  case, so gate 1 validates `Ybus` against PMD's power-flow physics instead.
+* Optional analytic Jacobian blocks (unit-tested against AD) could replace
+  `ForwardDiff` for speed on large feeders.

docs/src/literal_pmdse.md

@@ -0,0 +1,174 @@
+# Literal PMDSE — explicit matrix-based estimator
+
+As of version 0.8.0, PMDSE ships **Literal PMDSE**: a textbook, **JuMP-free**
+distribution-system state estimator that mirrors the 4-wire IVR explicit-neutral
+model (`IVRENPowerModel`, see [Explicit Neutral Models for DSSE](@ref))
+but solves the DSSE with **explicit matrices** — the measurement model
+`z = h(x) + e`, the Jacobian `H = ∂h/∂x`, the gain `G = HᵀWH` and the normal
+equations `G Δx = HᵀW r` — instead of handing a nonlinear program to a solver
+such as Ipopt.
+
+It consumes the **same** PMD mathematical data dictionary as
+`solve_ivr_en_mc_se` (`data_math`, `data_math["meas"]`, `data_math["se_settings"]`),
+so it is a drop-in alternative whenever you want a fast, transparent Gauss–Newton
+estimate, a robust WLAV estimate, or a general maximum-likelihood estimate, and
+full control over the angular/zero-sequence reference (see
+[How Literal PMDSE resolves the angular & zero-sequence reference problem](@ref)).
+
+## Quick start
+
+```julia
+import Ipopt
+import PowerModelsDistribution as _PMD
+import PowerModelsDistributionStateEstimation as _PMDSE
+
+# --- parse a 4-wire explicit-neutral feeder and build ground truth ----------
+eng = _PMD.parse_file(joinpath(_PMDSE.BASE_DIR, "test", "data",
+                               "three-bus-en-models", "3bus_4wire.dss"))
+_PMD.transform_loops!(eng)
+_PMD.remove_all_bounds!(eng)
+math = _PMD.transform_data_model(eng, kron_reduce=false, phase_project=false)
+_PMD.add_start_vrvi!(math)
+pf = _PMD.solve_mc_opf(math, _PMD.IVRENPowerModel, Ipopt.Optimizer)
+
+# --- synthesize the canonical IVREN measurement set (vr/vi + crd/cid + …) ---
+msr = joinpath(mktempdir(), "msr.csv")
+_PMDSE.write_measurements!(_PMD.IVRENPowerModel, math, pf, msr, σ=0.005)
+_PMDSE.add_measurements!(math, msr, actual_meas=true)
+math["se_settings"] = Dict{String,Any}("rescaler" => 1.0)
+
+# --- run the literal estimator (no JuMP / no solver object needed) -----------
+res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:sota)
+
+@show res.termination          # :converged
+@show res.iterations           # ~2-3 Gauss-Newton iterations
+@show res.solution["bus"]["3"]["vm"]   # estimated phase-to-ground magnitudes
+```
+
+## `solve_mc_se_literal`
+
+```julia
+solve_mc_se_literal(data; estimator = :wls, reference = :sota,
+                    ref_values = nothing, maxiter = 50, tol = 1e-9,
+                    verbose = false)
+```
+
+| argument | values | meaning |
+|----------|--------|---------|
+| `data` | `Dict` | PMD **mathematical** dictionary with `data["meas"]` populated (see [Measurement Conversion](@ref)) and an optional `data["se_settings"]`. |
+| `estimator` | `:wls` (default), `:wlav`, `:mle` | which estimator to run (see below). |
+| `reference` | `:sota` (default), `:full_slack`, `:prop` | reference / observability scheme (see below). |
+| `ref_values` | `(vr::Dict, vi::Dict)` or `nothing` | reference-bus phasor (keyed by terminal) for `:full_slack`. |
+| `maxiter` | `Int` | maximum Gauss-Newton / IRLS / Newton iterations. |
+| `tol` | `Float64` | convergence tolerance on `‖Δx‖∞`. |
+| `verbose` | `Bool` | print per-iteration progress. |
+
+`data["se_settings"]` is honoured: `"rescaler"` scales the weights
+(`W = diag(1/(rescaler·σ)²)`, consistent with `constraint_mc_residual`), and an
+optional `"reference"` entry overrides the `reference` keyword.
+
+### Estimators
+
+| `estimator` | method | use when |
+|-------------|--------|----------|
+| `:wls`  | Gauss–Newton **Weighted Least Squares** (normal equations with a QR / orthogonal fallback for ill-conditioned gains) | the default; fastest, optimal for Gaussian noise. |
+| `:wlav` | **Weighted Least Absolute Value** via IRLS | robustness to a single gross/bad measurement. |
+| `:mle`  | general **Maximum Likelihood** by Fisher scoring on `Σ logpdf(dstᵢ, hᵢ(x))` | non-Gaussian measurement errors; reduces **exactly** to WLS for Gaussian `dst`. |
+
+```julia
+res_wls  = _PMDSE.solve_mc_se_literal(math; estimator=:wls)
+res_wlav = _PMDSE.solve_mc_se_literal(math; estimator=:wlav)   # rejects gross errors
+res_mle  = _PMDSE.solve_mc_se_literal(math; estimator=:mle)    # == WLS for Gaussian dst
+```
+
+### Reference schemes
+
+The reference bus is the bus with `bus_type == 3`. Grounded terminals
+(`bus["grounded"]`) are always fixed to `0`. On top of that:
+
+| `reference` | fixes | rank | use |
+|-------------|-------|------|-----|
+| `:sota` (default) | reference-bus neutral `vr=vi=0` **and** one angle datum `vi[ref, phase₁]=0` | full | the SOTA combination — observable **without** forcing a balanced reference bus. |
+| `:full_slack` | the whole reference-bus phasor (`ref_values`, or balanced if omitted) | full | apples-to-apples comparison against a fixed-slack solver. |
+| `:prop` | reference-bus neutral only | **deficient by 1** | the "neutral-only" proposition — *unobservable*; for observability studies. |
+
+```julia
+# SOTA: ground neutral + one angle datum (estimates the unbalanced reference bus)
+res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:sota)
+
+# Fix the full reference-bus phasor to known (e.g. measured / true) values
+rb   = _PMDSE.literal_ref_bus(math)
+bsol = pf["solution"]["bus"]["$rb"]; terms = math["bus"]["$rb"]["terminals"]
+rv   = (Dict(t=>bsol["vr"][i] for (i,t) in enumerate(terms)),
+        Dict(t=>bsol["vi"][i] for (i,t) in enumerate(terms)))
+res  = _PMDSE.solve_mc_se_literal(math; reference=:full_slack, ref_values=rv)
+```
+
+See [How Literal PMDSE resolves the angular & zero-sequence reference problem](@ref)
+for the theory behind `:prop` vs `:sota`.
+
+## The result object
+
+`solve_mc_se_literal` returns a `LiteralResult`:
+
+| field | description |
+|-------|-------------|
+| `termination` | `:converged` or `:maxiter` |
+| `iterations`  | number of iterations |
+| `objective`   | weighted SSR `rᵀWr` (WLS/MLE) or `Σ√wᵢ|rᵢ|` (WLAV) at the solution |
+| `solve_time`  | wall-clock seconds |
+| `gain_cond`, `gain_rank` | condition number / numerical rank of `G = HᵀWH` (health check) |
+| `solution`    | `Dict("bus" => id => Dict("vr","vi","vm"))`, vectors ordered like `bus["terminals"]` |
+| `x_free`, `x_full` | the estimated state (free entries / full `[vr; vi]`) |
+| `estimator`   | `:wls`, `:wlav` or `:mle` |
+
+### Accuracy vs ground truth
+
+```julia
+m = _PMDSE.accuracy_metrics(res, pf["solution"])      # (rmse, maxerr, n) over all |U|
+m = _PMDSE.accuracy_metrics(res, pf["solution"]; include_neutral=false)  # phases only
+e = _PMDSE.voltage_errors(res, pf["solution"])        # per-node |ΔU| vector (incl. NEV)
+```
+
+## Supported measurements
+
+Literal PMDSE implements `h(x)` for every measurement var-type that
+`IVRENPowerModel` supports (see [Measurement Conversion](@ref)):
+
+* **native**: `:vr, :vi, :cr, :ci, :crd, :cid, :crg, :cig`
+* **voltage**: `:vm`/`:vmn` (phase-to-neutral), `:va`, `:vll`
+* **branch**: `:cm, :ca, :p, :q`
+* **injection**: `:pd, :qd, :pg, :qg, :ptot, :qtot, :cmd, :cmg, :cad, :cag`
+
+Voltage and power measurements follow the EN phase-to-neutral convention and the
+neutral residual row is skipped (`setdiff(conns, _N_IDX)`), exactly as in
+`constraint_mc_residual`. Current/power *injection* measurements are mapped to the
+nodal injection `Ybus·U`; when several components share a bus terminal (e.g.
+multiple loads on one bus) their measurements are aggregated into one nodal
+residual — exact at the solution.
+
+## Lower-level API (advanced)
+
+For custom workflows you can build and reuse the pieces directly:
+
+```julia
+lm    = _PMDSE.LiteralModel(math; reference=:sota)   # node map, Ybus, references
+model = _PMDSE.build_se_model(lm, math; rescaler=1.0)# residual rows z, w, h-closures
+res   = _PMDSE.solve_wls(model)                       # or solve_wlav / solve_mle
+
+H = ForwardDiff.jacobian(x -> _PMDSE.predict(model, x), res.x_free)  # the Jacobian
+```
+
+`build_se_model` accepts `rescaler`, `σ_floor` (prevents an infinite weight when a
+measured value is exactly 0, e.g. `vi` at a balanced reference bus) and
+`zero_inj_sigma` (weight of the zero-injection pseudo-measurements).
+
+## Notes & limitations
+
+* Transformers are not yet folded into `Ybus` (the explicit-neutral test feeders
+  use a stiff voltage source); mutually-coupled branches and shunts are supported.
+* PMD's `calc_admittance_matrix` does not support the 4-wire
+  (`kron_reduce=false`) case, so `Ybus` is validated against PMD's power-flow
+  current balance in the test suite.
+* The Jacobian is assembled with `ForwardDiff`, keeping the estimation algebra
+  fully explicit while avoiding hand-derivative bugs.