diff --git a/docs/make.jl b/docs/make.jl index 9122b81..d3525fe 100644 --- a/docs/make.jl +++ b/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 " => [ diff --git a/docs/src/Literal_PMDSE_PLAN.md b/docs/src/Literal_PMDSE_PLAN.md new file mode 100644 index 0000000..15b754f --- /dev/null +++ b/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. diff --git a/docs/src/literal_pmdse.md b/docs/src/literal_pmdse.md new file mode 100644 index 0000000..fbda7b1 --- /dev/null +++ b/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. diff --git a/docs/src/literal_pmdse_observability.md b/docs/src/literal_pmdse_observability.md new file mode 100644 index 0000000..93e70a7 --- /dev/null +++ b/docs/src/literal_pmdse_observability.md @@ -0,0 +1,129 @@ +# How Literal PMDSE resolves the angular & zero-sequence reference problem + +This page connects the [Literal PMDSE — explicit matrix-based estimator](@ref) +estimator to the theory of the report *"Rigorous Mathematical Analysis of the +Angular and Zero-Sequence Reference Problem in 4-Wire Unbalanced IVR Power Flow"* +(`main.tex`, TSK-871). Literal PMDSE is the explicit, matrix-based estimator that +**operationalizes** that report's conclusion and turns its central theorem into a +runnable regression test. + +## The problem (recap of `main.tex`) + +In a 4-wire explicit-neutral IVR model the state is the rectangular nodal voltage +``\mathbf{U}\in\mathbb{C}^{4n}`` and the measurement set of interest consists of +power / current-magnitude / voltage-magnitude quantities plus the explicit +grounding of the reference-bus neutral ``U_{s,n}=0``. + +**Theorem 1 (Rotational Invariance).** *Grounding only the neutral leaves the +Jacobian singular, rank-deficient by exactly one.* A global phase rotation +``\mathbf{U}\mapsto \mathbf{U}\,e^{j\alpha}`` (i) preserves the neutral constraint +(``0\cdot e^{j\alpha}=0``) and (ii) leaves every complex power injection +invariant (``S_{k,p}=U_{k,p}\,(\sum Y\,U)^\ast`` ⇒ the ``e^{j\alpha}`` and +``e^{-j\alpha}`` cancel). Hence there is a continuous locus of solutions and the +measurement Jacobian ``H`` has the null vector + +```math +\frac{\partial \mathbf{U}(\alpha)}{\partial\alpha}\Big|_{\alpha=0} + = j\,\mathbf{U}^\ast + \;\;\Longleftrightarrow\;\; + \delta\mathbf{x}_{\mathrm{rot}}=\begin{bmatrix}-\mathbf{U}^{\mathrm i}\\[2pt]\;\;\mathbf{U}^{\mathrm r}\end{bmatrix}. +``` + +Because ``H`` is not full column rank, the gain ``G=H^\top W H`` is **singular** +and the system is **unobservable**. + +**Resolution.** Grounding the neutral fixes only the *translational / zero-sequence* +datum. The *rotational / positive-sequence* datum must be fixed separately by one +phase-angle reference ``U^{\mathrm i}_{s,a}=0``. The SOTA estimator therefore +**combines** both: + +1. ground the reference-bus neutral — ``U^{\mathrm r}_{s,n}=U^{\mathrm i}_{s,n}=0``; +2. fix a single angle datum — ``U^{\mathrm i}_{s,a}=0``. + +## How Literal PMDSE implements it + +Literal PMDSE exposes the two constraints above (and the failing proposition) as +the configurable `reference` keyword of `solve_mc_se_literal` / `LiteralModel`. +The reference partition removes the fixed components as *columns* of ``H`` (they +become known parameters), so the remaining gain is full-rank **without** assuming +a balanced reference bus: + +| `main.tex` case | Literal PMDSE `reference` | what is fixed | rank of ``H`` | +|-----------------|---------------------------|---------------|---------------| +| Proposition (neutral only) | `:prop` | ``U^{\mathrm r}_{s,n}=U^{\mathrm i}_{s,n}=0`` | **deficient by 1** | +| Ultimate solution (SOTA) | `:sota` | neutral **and** ``U^{\mathrm i}_{s,a}=0`` | full | +| Conventional fixed slack | `:full_slack` | the whole reference-bus phasor | full (forces balanced RB) | + +`:sota` is the default, and — exactly as the report prescribes — it leaves the +reference-bus phase ``b``/``c`` magnitudes and angles **free state variables**, so +the unbalanced reference bus and the neutral-to-earth voltage (NEV) are estimated +rather than assumed. + +## The theorem as a regression test + +The report's negative result is encoded directly in `test/literal/test_references.jl`. +Evaluate the Jacobian ``H`` at the true state and probe it along the analytic +rotation null vector ``\delta\mathbf{x}_{\mathrm{rot}}=[-\mathbf{U}^{\mathrm i};\, +\mathbf{U}^{\mathrm r}]`` (restricted to the free coordinates): + +```julia +import PowerModelsDistributionStateEstimation as _PMDSE +import ForwardDiff, LinearAlgebra + +# rotation-invariant measurements (|U|, P, Q) on the toy 3-bus 4-wire feeder +math, pf = toy_pf() # see test/literal/helpers.jl +M = rotation_invariant_meas(math, pf) + +for reference in (:prop, :sota) + lm = _PMDSE.LiteralModel(math; reference=reference) + model = _PMDSE.build_se_model(lm, with_meas(math, M)) + xf = pf_xfree(lm, pf) + H = ForwardDiff.jacobian(x -> _PMDSE.predict(model, x), xf) + d = rotation_dir(lm, _PMDSE.expand_state(lm, xf)) # δx_rot, free coords + @show reference, LinearAlgebra.norm(H*d) / (LinearAlgebra.norm(H)*LinearAlgebra.norm(d)) +end +``` + +which yields + +``` +(:prop, 2.56e-17) # H·δx_rot ≈ 0 → rotation is UNOBSERVABLE (Theorem 1) +(:sota, 5.12e-02) # adding the angle datum makes the rotation observable +``` + +The rotation is a numerical null vector of ``H`` under `:prop` and is removed +under `:sota` — a direct, quantitative reproduction of Theorem 1 and its +resolution. + +## The toy 3-bus 4-wire example + +`main.tex` §6 analyses a 3-bus, 4-wire radial network +(``\mathbf{Z}_{abc}=\mathrm{diag}(0.01+j0.02)\,\Omega``, +``\mathbf{Z}_n=0.05+j0.01\,\Omega``, 100 kVA / 400 V) and reports +``\mathrm{rank}(\mathbf J)=23`` (singular) with neutral-only grounding, +restored to ``\mathrm{rank}(\mathbf J)=24`` once ``U^{\mathrm i}_{1,a}=0`` is +added. The same feeder is reproduced in `test/literal/helpers.jl` (`toy_math`) +and used by the observability test above, so the report's worked example and the +estimator share one network. + +## Validation against the JuMP IVREN estimator + +Beyond observability, Literal PMDSE is validated to **agree** with the +optimization-based explicit-neutral estimator (`solve_ivr_en_mc_se` / +`IVRENPowerModel`). On a noiseless, consistent measurement set with identical +references, the literal WLS reproduces the IVREN / power-flow state to machine +precision (`max|ΔU| ~ 1e-16` for both `:full_slack` and `:sota`), while +converging in ~2 Gauss–Newton iterations and running roughly two orders of +magnitude faster than the Ipopt solve (see `test/literal/test_benchmark.jl`). + +## Summary + +| concern (`main.tex`) | mechanism | Literal PMDSE | +|----------------------|-----------|---------------| +| zero-sequence / translational datum | ground reference-bus neutral | always (grounded terminals fixed to 0) | +| positive-sequence / rotational datum | fix one phase angle | `:sota` adds ``U^{\mathrm i}_{s,a}=0`` | +| unbalanced reference bus (NEV) | leave ``b,c`` free | `:sota` estimates them | +| the negative result (singular gain) | rotational invariance | reproduced by `:prop` (rank-deficient by 1) | + +See also: [Angular Reference Models](@ref) and +[Explicit Neutral Models for DSSE](@ref). diff --git a/examples/literal_pmdse_example.jl b/examples/literal_pmdse_example.jl new file mode 100644 index 0000000..45753bb --- /dev/null +++ b/examples/literal_pmdse_example.jl @@ -0,0 +1,110 @@ +################################################################################ +# Literal PMDSE — worked example # +# # +# Demonstrates the explicit, JuMP-free explicit-neutral state estimator # +# `solve_mc_se_literal` with every estimator (:wls/:wlav/:mle) and reference # +# scheme (:sota/:full_slack/:prop), result inspection, accuracy metrics and # +# the lower-level API. # +# # +# Run with a project that has Ipopt + PowerModelsDistribution(StateEstimation):# +# julia --project examples/literal_pmdse_example.jl # +################################################################################ +import Ipopt +import ForwardDiff, LinearAlgebra +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE + +ipopt = _PMDSE.optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>0, "tol"=>1e-10) + +# ---------------------------------------------------------------------------- # +# 1. Ground truth: parse a 4-wire explicit-neutral feeder and run a power flow # +# ---------------------------------------------------------------------------- # +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) + +# ---------------------------------------------------------------------------- # +# 2. Synthesize the canonical IVREN measurement set (vr/vi + crd/cid + crg/cig) # +# `actual_meas=true` keeps the exact PF values (noiseless); use a seed for # +# sampled noise (see step 6). # +# ---------------------------------------------------------------------------- # +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) # weights = 1/(rescaler·σ)² + +# ---------------------------------------------------------------------------- # +# 3. Run every estimator with the default SOTA reference # +# ---------------------------------------------------------------------------- # +println("="^78, "\n3. estimators (reference = :sota)\n", "="^78) +for est in (:wls, :wlav, :mle) + res = _PMDSE.solve_mc_se_literal(math; estimator=est, reference=:sota, + maxiter=50, tol=1e-9, verbose=false) + m = _PMDSE.accuracy_metrics(res, pf["solution"]) + println(rpad(uppercase(string(est)), 5), + " term=", res.termination, " iters=", res.iterations, + " rmse|U|=", round(m.rmse, sigdigits=3), + " max|U|=", round(m.maxerr, sigdigits=3), + " t=", round(res.solve_time, sigdigits=3), "s") +end + +# ---------------------------------------------------------------------------- # +# 4. Reference schemes # +# :sota — ground neutral + one angle datum (estimates the unbalanced RB)# +# :full_slack — fix the whole reference-bus phasor (here, to the true values) # +# :prop — neutral only -> rank-deficient (unobservable); shown for study# +# ---------------------------------------------------------------------------- # +println("\n", "="^78, "\n4. reference schemes (estimator = :wls)\n", "="^78) +rb = _PMDSE.literal_ref_bus(math) +bsol = pf["solution"]["bus"]["$rb"]; terms = math["bus"]["$rb"]["terminals"] +refvals = (Dict(t => bsol["vr"][i] for (i, t) in enumerate(terms)), + Dict(t => bsol["vi"][i] for (i, t) in enumerate(terms))) + +for (ref, rv) in ((:sota, nothing), (:full_slack, refvals)) + res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=ref, ref_values=rv) + m = _PMDSE.accuracy_metrics(res, pf["solution"]) + println(rpad(string(ref), 11), " term=", res.termination, + " max|U vs PF|=", round(m.maxerr, sigdigits=3), + " cond(G)=", round(res.gain_cond, sigdigits=3)) +end + +# ---------------------------------------------------------------------------- # +# 5. Inspect a result object # +# ---------------------------------------------------------------------------- # +println("\n", "="^78, "\n5. result inspection\n", "="^78) +res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:full_slack, ref_values=refvals) +println("termination = ", res.termination, " iterations = ", res.iterations) +println("objective = ", round(res.objective, sigdigits=3)) +println("gain rank = ", res.gain_rank) +println("loadbus |U| (per terminal incl. neutral) = ", round.(res.solution["bus"]["3"]["vm"], digits=5)) +println("per-node |ΔU| vs PF (first 6) = ", round.(_PMDSE.voltage_errors(res, pf["solution"])[1:6], sigdigits=2)) + +# ---------------------------------------------------------------------------- # +# 6. Robustness: a single gross error — WLAV vs WLS # +# ---------------------------------------------------------------------------- # +println("\n", "="^78, "\n6. one gross error: WLAV is robust, WLS is not\n", "="^78) +mathb = deepcopy(math) +mkey = first(k for (k, mm) in mathb["meas"] if mm["var"] == :vr) # corrupt a vr meas +import Distributions as _DST +mathb["meas"][mkey]["dst"][1] = _DST.Normal(_DST.mean(mathb["meas"][mkey]["dst"][1]) + 0.3, 0.005) +for est in (:wls, :wlav) + res = _PMDSE.solve_mc_se_literal(mathb; estimator=est, reference=:full_slack, ref_values=refvals) + println(rpad(uppercase(string(est)), 5), " max|U vs PF| = ", + round(_PMDSE.accuracy_metrics(res, pf["solution"]).maxerr, sigdigits=3)) +end + +# ---------------------------------------------------------------------------- # +# 7. Lower-level API: build the model once, reuse it, get the Jacobian # +# ---------------------------------------------------------------------------- # +println("\n", "="^78, "\n7. lower-level API\n", "="^78) +lm = _PMDSE.LiteralModel(math; reference=:full_slack, ref_values=refvals) +model = _PMDSE.build_se_model(lm, math; rescaler=1.0) +res = _PMDSE.solve_wls(model) +H = ForwardDiff.jacobian(x -> _PMDSE.predict(model, x), res.x_free) +println("residual rows = ", length(model.z), " free state dim = ", length(lm.free_idx)) +println("size(H) = ", size(H), " cond(HᵀWH) = ", round(res.gain_cond, sigdigits=3)) +println("\nDONE.") diff --git a/src/PowerModelsDistributionStateEstimation.jl b/src/PowerModelsDistributionStateEstimation.jl index 88efdcf..c992259 100644 --- a/src/PowerModelsDistributionStateEstimation.jl +++ b/src/PowerModelsDistributionStateEstimation.jl @@ -84,6 +84,14 @@ include("io/postprocessing.jl") include("prob/se.jl") include("prob/se_en.jl") +# Literal PMDSE: explicit matrix-based explicit-neutral state estimator (JuMP-free) +include("bare/literal_core.jl") +include("bare/measurement_model.jl") +include("bare/solve_wls.jl") +include("bare/solve_wlav.jl") +include("bare/mle.jl") +include("bare/literal_pmdse.jl") + # export include("core/export.jl") diff --git a/src/bare/literal_core.jl b/src/bare/literal_core.jl new file mode 100644 index 0000000..40ddccd --- /dev/null +++ b/src/bare/literal_core.jl @@ -0,0 +1,314 @@ +################################################################################ +# Literal PMDSE - explicit matrix-based distribution system state estimation # +# # +# literal_core.jl : node index map, nodal admittance (Y_bus) assembly, # +# reference/observability handling and flat start. # +# # +# The "literal" estimator mirrors the 4-wire IVR explicit-neutral model # +# (`IVRENPowerModel`, see src/prob/se_en.jl) but keeps only the nodal bus # +# voltages `x = [vr; vi]` as state and solves the textbook normal equations # +# (Abur & Exposito) instead of building a JuMP model. Branch / load / gen # +# currents are eliminated through the bus admittance matrix. # +################################################################################ + +"index of the neutral conductor (kept explicit in the 4-wire EN model)" +# (`_N_IDX` is already defined in the package main module) + +""" + LiteralModel + +Container with everything the measurement model and the solvers need. The full +real state is laid out as `x = [vr_1..vr_n , vi_1..vi_n]` where node `k` +(`1:n`) is the `k`-th `(bus, terminal)` pair in `nodes`. The complex bus +admittance is stored split into real/imaginary parts `Ybus = G + im*B` so that +the nodal current injection `I = Ybus*U` is evaluated with `ForwardDiff`-friendly +real matrix products + + Ir = G*vr - B*vi , Ii = B*vr + G*vi . + +# Fields +- `nodes` : ordered vector of `(bus_id, terminal)` pairs +- `node_index` : map `(bus_id, terminal) -> position in 1:n` +- `n` : number of nodes +- `G`, `B` : real / imaginary part of `Ybus` (n×n, dense Float64) +- `ref_bus` : id of the reference bus (`bus_type == 3`) +- `reference` : reference scheme `Symbol` (`:full_slack`, `:sota`, `:prop`) +- `fixed_mask` : `BitVector` of length `2n`, `true` where a state entry is fixed +- `x_fixed` : full `2n` vector holding the fixed values (0 where free) +- `free_idx` : indices (into `1:2n`) of the free state entries +- `inj_nodes` : nodes (positions) that carry a nodal-current-injection equation +""" +struct LiteralModel + math::Dict{String,Any} + nodes::Vector{Tuple{Int,Int}} + node_index::Dict{Tuple{Int,Int},Int} + n::Int + G::Matrix{Float64} + B::Matrix{Float64} + ref_bus::Int + reference::Symbol + fixed_mask::BitVector + x_fixed::Vector{Float64} + free_idx::Vector{Int} + inj_nodes::Vector{Int} +end + +"return the reference bus id (the unique `bus_type == 3` bus)" +function literal_ref_bus(math::Dict) + refs = [bus["index"] for (_, bus) in math["bus"] if bus["bus_type"] == 3] + isempty(refs) && error("Literal PMDSE: no reference bus (bus_type == 3) found") + length(refs) > 1 && @warn "Literal PMDSE: multiple reference buses found, using $(refs[1])" + return refs[1] +end + +"build the ordered `(bus, terminal)` node list and the inverse lookup" +function literal_node_map(math::Dict) + nodes = Tuple{Int,Int}[] + for b in sort(collect(keys(math["bus"])); by = x -> parse(Int, x)) + bus = math["bus"][b] + for t in bus["terminals"] + push!(nodes, (bus["index"], t)) + end + end + node_index = Dict{Tuple{Int,Int},Int}(nd => k for (k, nd) in enumerate(nodes)) + return nodes, node_index +end + +""" + build_ybus(math, node_index, n) -> Matrix{ComplexF64} + +Assemble the nodal admittance matrix of the explicit-neutral network. Every +branch contributes its full (mutually-coupled) series admittance `Y = Z^{-1}`, +`Z = br_r + im*br_x`, stamped over the `f_connections`/`t_connections` +terminals, plus its Pi-model line shunts (`g_fr+im*b_fr`, `g_to+im*b_to`). Bus +shunts (`math["shunt"]`) are added on the diagonal block. This reproduces the +admittance implied by `constraint_mc_bus_voltage_drop` / +`constraint_mc_current_balance_se(::IVRENPowerModel, …)`. +""" +function build_ybus(math::Dict, node_index::Dict{Tuple{Int,Int},Int}, n::Int) + Ybus = zeros(ComplexF64, n, n) + + for (_, br) in math["branch"] + get(br, "br_status", 1) == 0 && continue + f_bus = br["f_bus"]; t_bus = br["t_bus"] + fc = br["f_connections"]; tc = br["t_connections"] + Z = Matrix{ComplexF64}(br["br_r"] .+ im .* br["br_x"]) + Ybr = inv(Z) + m = length(fc) + Ysh_fr = haskey(br, "g_fr") ? Matrix{ComplexF64}(br["g_fr"] .+ im .* br["b_fr"]) : zeros(ComplexF64, m, m) + Ysh_to = haskey(br, "g_to") ? Matrix{ComplexF64}(br["g_to"] .+ im .* br["b_to"]) : zeros(ComplexF64, m, m) + + F = [node_index[(f_bus, c)] for c in fc] + T = [node_index[(t_bus, c)] for c in tc] + + Ybus[F, F] .+= Ybr .+ Ysh_fr + Ybus[T, T] .+= Ybr .+ Ysh_to + Ybus[F, T] .-= Ybr + Ybus[T, F] .-= Ybr + end + + for (_, sh) in get(math, "shunt", Dict{String,Any}()) + get(sh, "status", 1) == 0 && continue + sb = sh["shunt_bus"] + conns = sh["connections"] + Ysh = Matrix{ComplexF64}(sh["gs"] .+ im .* sh["bs"]) + S = [node_index[(sb, c)] for c in conns] + Ybus[S, S] .+= Ysh + end + + return Ybus +end + +""" + reference_partition(math, nodes, node_index, n, ref_bus; reference, ref_values) + +Decide which state entries are *fixed* (known parameters, removed as columns of +the Jacobian) versus *free*. Returns `(fixed_mask, x_fixed)` of length `2n`. + +Reference schemes (configurable, see §5 of the plan / `main.tex`): +- `:full_slack` : fix the whole reference-bus phasor (sanity check; forces a + balanced reference bus). +- `:sota` : ground the reference-bus neutral *and* fix one phase-angle + datum `vi[ref, first_phase] = 0` (the SOTA combination that makes the gain + full-rank without biasing the reference bus). +- `:prop` : ground only the neutral(s); **rank-deficient by one** – used to + reproduce the observability trap of `main.tex`. + +Grounded terminals (`bus["grounded"]`) are always fixed to 0 in every scheme. +`ref_values` optionally provides `(vr, vi)` per reference-bus terminal (e.g. the +true slack phasor) for `:full_slack`; otherwise a balanced phasor scaled by the +bus base voltage start is used. +""" +function reference_partition(math::Dict, nodes, node_index, n::Int, ref_bus::Int; + reference::Symbol = :sota, ref_values = nothing) + fixed_mask = falses(2n) + x_fixed = zeros(Float64, 2n) + + # 1) grounded terminals -> hard zero everywhere + for (_, bus) in math["bus"] + for (idx, t) in enumerate(bus["terminals"]) + if bus["grounded"][idx] + k = node_index[(bus["index"], t)] + fixed_mask[k] = true; x_fixed[k] = 0.0 + fixed_mask[n+k] = true; x_fixed[n+k] = 0.0 + end + end + end + + rbus = math["bus"][string(ref_bus)] + rterms = rbus["terminals"] + phases = [t for t in rterms if t != _N_IDX] + + if reference == :full_slack + # fix the full reference-bus phasor + if ref_values === nothing + va = deg2rad.([0.0, -120.0, 120.0]) + vr_ref = Dict{Int,Float64}(); vi_ref = Dict{Int,Float64}() + for (j, t) in enumerate(phases) + vr_ref[t] = cos(va[mod1(j, 3)]); vi_ref[t] = sin(va[mod1(j, 3)]) + end + for t in rterms + t == _N_IDX && (vr_ref[t] = 0.0; vi_ref[t] = 0.0) + end + else + vr_ref, vi_ref = ref_values + end + for t in rterms + k = node_index[(ref_bus, t)] + fixed_mask[k] = true; x_fixed[k] = vr_ref[t] + fixed_mask[n+k] = true; x_fixed[n+k] = vi_ref[t] + end + elseif reference == :sota + # ground reference-bus neutral + if _N_IDX in rterms + k = node_index[(ref_bus, _N_IDX)] + fixed_mask[k] = true; x_fixed[k] = 0.0 + fixed_mask[n+k] = true; x_fixed[n+k] = 0.0 + end + # fix one phase-angle datum: vi[ref, first phase] = 0 + kp = node_index[(ref_bus, first(phases))] + fixed_mask[n+kp] = true; x_fixed[n+kp] = 0.0 + elseif reference == :prop + # only ground reference-bus neutral (rank-deficient by one) + if _N_IDX in rterms + k = node_index[(ref_bus, _N_IDX)] + fixed_mask[k] = true; x_fixed[k] = 0.0 + fixed_mask[n+k] = true; x_fixed[n+k] = 0.0 + end + else + error("Literal PMDSE: unknown reference scheme :$(reference)") + end + + return fixed_mask, x_fixed +end + +"nodes (positions) that carry a nodal current-injection equation: every +ungrounded terminal that is **not** on the reference bus (the reference bus is a +free current slack)." +function injection_node_set(math::Dict, nodes, node_index, ref_bus::Int) + grounded = Set{Tuple{Int,Int}}() + for (_, bus) in math["bus"] + for (idx, t) in enumerate(bus["terminals"]) + bus["grounded"][idx] && push!(grounded, (bus["index"], t)) + end + end + inj = Int[] + for (k, (b, t)) in enumerate(nodes) + (b == ref_bus) && continue + ((b, t) in grounded) && continue + push!(inj, k) + end + return inj +end + +""" + flat_start(math, nodes; vbase) -> Vector{Float64} + +Balanced flat start: phases at `1∠{0,-120,120}`, neutral ≈ 0. When the bus data +carries `vr_start`/`vi_start` (set by `_PMD.add_start_vrvi!`) those are used. +""" +function flat_start(math::Dict, nodes, node_index, n::Int) + x = zeros(Float64, 2n) + va = deg2rad.([0.0, -120.0, 120.0]) + for (k, (b, t)) in enumerate(nodes) + bus = math["bus"][string(b)] + if haskey(bus, "vr_start") && haskey(bus, "vi_start") + idx = findfirst(isequal(t), bus["terminals"]) + x[k] = bus["vr_start"][idx] + x[n+k] = bus["vi_start"][idx] + elseif t != _N_IDX + j = findfirst(isequal(t), [1, 2, 3]) + ang = j === nothing ? 0.0 : va[j] + x[k] = cos(ang) + x[n+k] = sin(ang) + end + end + return x +end + +""" + LiteralModel(math; reference=:sota, ref_values=nothing) + +Build the full literal model from a PMD *mathematical* data dictionary +(`data_math`, as produced by `_PMD.transform_data_model(...; kron_reduce=false, +phase_project=false)`). +""" +function LiteralModel(math::Dict; reference::Symbol = :sota, ref_values = nothing) + nodes, node_index = literal_node_map(math) + n = length(nodes) + Ybus = build_ybus(math, node_index, n) + ref_bus = literal_ref_bus(math) + fixed_mask, x_fixed = reference_partition(math, nodes, node_index, n, ref_bus; + reference = reference, ref_values = ref_values) + free_idx = [i for i in 1:2n if !fixed_mask[i]] + inj_nodes = injection_node_set(math, nodes, node_index, ref_bus) + return LiteralModel(math, nodes, node_index, n, real.(Ybus), imag.(Ybus), + ref_bus, reference, fixed_mask, x_fixed, free_idx, inj_nodes) +end + +"reconstruct the full `2n` real state from the free entries" +@inline function expand_state(lm::LiteralModel, x_free::AbstractVector{T}) where {T} + x = convert(Vector{T}, lm.x_fixed) + @inbounds for (j, i) in enumerate(lm.free_idx) + x[i] = x_free[j] + end + return x +end + +"split a full `2n` real state into `(vr, vi)` complex-part vectors of length n" +@inline function vrvi(lm::LiteralModel, x::AbstractVector) + return @views x[1:lm.n], x[lm.n+1:2lm.n] +end + +"nodal current injection `I = Ybus*U` split into real/imag parts (ForwardDiff ok)" +@inline function nodal_injection(lm::LiteralModel, vr::AbstractVector, vi::AbstractVector) + Ir = lm.G * vr .- lm.B * vi + Ii = lm.B * vr .+ lm.G * vi + return Ir, Ii +end + +"initial free-state vector from the (balanced) flat start" +function flat_start_free(lm::LiteralModel) + x0 = flat_start(lm.math, lm.nodes, lm.node_index, lm.n) + return x0[lm.free_idx] +end + +""" + state_solution(lm, x_full) -> Dict + +Convert a full real state vector into a PMD-style solution dictionary +`Dict("bus" => Dict(id => Dict("vr"=>[...], "vi"=>[...])))`, with the per-bus +vectors ordered like `bus["terminals"]`. +""" +function state_solution(lm::LiteralModel, x_full::AbstractVector) + vr, vi = vrvi(lm, x_full) + busdict = Dict{String,Any}() + for (b, bus) in lm.math["bus"] + terms = bus["terminals"] + vrb = [vr[lm.node_index[(bus["index"], t)]] for t in terms] + vib = [vi[lm.node_index[(bus["index"], t)]] for t in terms] + busdict[b] = Dict{String,Any}("vr" => vrb, "vi" => vib, + "vm" => sqrt.(vrb .^ 2 .+ vib .^ 2)) + end + return Dict{String,Any}("bus" => busdict) +end diff --git a/src/bare/literal_pmdse.jl b/src/bare/literal_pmdse.jl new file mode 100644 index 0000000..c56c567 --- /dev/null +++ b/src/bare/literal_pmdse.jl @@ -0,0 +1,76 @@ +################################################################################ +# Literal PMDSE # +# literal_pmdse.jl : user-facing entry point `solve_mc_se_literal`. # +# # +# Solves the 4-wire IVR explicit-neutral DSSE with explicit matrices # +# (Gauss-Newton WLS / IRLS WLAV / general MLE) on the *same* PMD mathematical # +# data dictionary (`data_math` + `data_math["meas"]` + `data_math["se_settings"]`)# +# consumed by `solve_ivr_en_mc_se` (`IVRENPowerModel`). No JuMP is used. # +################################################################################ + +""" + solve_mc_se_literal(data; estimator=:wls, reference=:sota, ref_values=nothing, + maxiter=50, tol=1e-9, verbose=false, kwargs...) + +Run the literal (matrix-based, JuMP-free) explicit-neutral state estimator. + +# Arguments +- `data` : PMD **mathematical** data dictionary, with `data["meas"]` populated + (e.g. via `add_measurements!`) and an optional `data["se_settings"]` + (`"rescaler"`, and optionally `"reference"`). +- `estimator` : `:wls` (default), `:wlav`, or `:mle`. +- `reference` : `:sota` (ground reference-bus neutral + one angle datum, default), + `:full_slack` (fix the whole reference-bus phasor), or `:prop` (neutral only – + rank-deficient, for observability studies). +- `ref_values` : optional `(vr_dict, vi_dict)` over reference-bus terminals for + `:full_slack`. + +Returns a [`LiteralResult`](@ref) with `solution`, `objective`, `iterations`, +`solve_time`, `gain_cond`, `gain_rank` and `termination`. +""" +function solve_mc_se_literal(data::Dict{String,<:Any}; estimator::Symbol = :wls, + reference::Symbol = :sota, ref_values = nothing, + maxiter::Int = 50, tol::Float64 = 1e-9, + verbose::Bool = false, kwargs...) + haskey(data, "meas") || error("solve_mc_se_literal: data has no \"meas\" dictionary") + se = get(data, "se_settings", Dict{String,Any}()) + rescaler = Float64(get(se, "rescaler", 1.0)) + reference = Symbol(get(se, "reference", reference)) + + lm = LiteralModel(data; reference = reference, ref_values = ref_values) + model = build_se_model(lm, data; rescaler = rescaler) + + if estimator == :wls + return solve_wls(model; maxiter = maxiter, tol = tol, verbose = verbose) + elseif estimator == :wlav + return solve_wlav(model; maxiter = maxiter, tol = tol, verbose = verbose) + elseif estimator == :mle + return solve_mle(model; maxiter = maxiter, tol = tol, verbose = verbose) + else + error("solve_mc_se_literal: unknown estimator :$(estimator) (use :wls, :wlav or :mle)") + end +end + +# ---- accuracy helpers vs a reference (PF ground truth or another solution) ---- + +"per-bus complex voltage error of a literal `result` against a PMD solution dict" +function voltage_errors(result::LiteralResult, ref_solution::Dict; include_neutral::Bool = true) + lm_sol = result.solution["bus"] + errs = Float64[] + for (b, bsol) in lm_sol + refb = ref_solution["bus"][b] + 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]))) + end + end + return errs +end + +"summary accuracy metrics (RMSE, max |U| error) of a literal result vs ground truth" +function accuracy_metrics(result::LiteralResult, ref_solution::Dict; include_neutral::Bool = true) + e = voltage_errors(result, ref_solution; include_neutral = include_neutral) + return (rmse = sqrt(sum(abs2, e) / length(e)), maxerr = maximum(e), n = length(e)) +end diff --git a/src/bare/measurement_model.jl b/src/bare/measurement_model.jl new file mode 100644 index 0000000..b8b9f17 --- /dev/null +++ b/src/bare/measurement_model.jl @@ -0,0 +1,297 @@ +################################################################################ +# Literal PMDSE # +# measurement_model.jl : rectangular explicit-neutral measurement functions # +# h(x) and assembly of the residual model. # +# # +# Every supported measurement var-type (see src/core/measurement_conversion.jl # +# for the IVR/IVREN authoritative definitions) is mapped to a closure of the # +# nodal state `x = [vr; vi]`. Branch / load / gen currents are eliminated via # +# the bus admittance `Ybus` (nodal current `I = Ybus*U`) and the branch # +# admittance (branch current `I_fr = Yff*U_f + Yft*U_t`). # +# # +# EN conventions: # +# * voltage / power measurements are phase-to-neutral (subtract `_N_IDX`) # +# * residuals are formed over `setdiff(connections, _N_IDX)` (neutral skipped) # +# * a wye component's neutral current is `-sum(phase currents)` # +################################################################################ + +"per-row evaluation context shared across the residual vector" +struct HCtx{T} + vr::Vector{T} + vi::Vector{T} + Ir::Vector{T} # real part of nodal injection Ybus*U + Ii::Vector{T} # imag part of nodal injection Ybus*U +end + +""" + SEModel + +Assembled residual model. `hfun(x_free) -> Vector` returns the model prediction +of every residual row; `z` are the observed values and `w` the diagonal weights +(`w_i = 1/(rescaler*σ_i)^2`). `rowmeta` documents each row. +""" +struct SEModel + lm::LiteralModel + rowfns::Vector{Any} # each maps HCtx -> scalar + z::Vector{Float64} + w::Vector{Float64} + rowmeta::Vector{NamedTuple} + brdata::Vector{Any} # precomputed branch admittance blocks (for branch rows) +end + +nrows(m::SEModel) = length(m.z) + +# ---- helpers --------------------------------------------------------------- + +"ordered non-neutral connections of a measurement's active connection list" +_nonneutral(active) = [c for c in active if c != _N_IDX] + +"mean (value) and std (sigma) of a per-conductor distribution entry" +function _zsigma(dst_entry, σ_exact::Float64) + if isa(dst_entry, _DST.Normal) + return _DST.mean(dst_entry), _DST.std(dst_entry) + elseif isa(dst_entry, Real) + return Float64(dst_entry), σ_exact # hard / fixed value + else + return _DST.mean(dst_entry), _DST.std(dst_entry) + end +end + +"active connections of a component for a given measurement var" +function _active_connections(math::Dict, cmp::Symbol, cmp_id::Int) + if cmp == :bus + return math["bus"][string(cmp_id)]["terminals"] + elseif cmp == :load + return math["load"][string(cmp_id)]["connections"] + elseif cmp == :gen + return math["gen"][string(cmp_id)]["connections"] + elseif cmp == :branch + return math["branch"][string(cmp_id)]["f_connections"] + else + error("Literal PMDSE: unsupported measurement component $(cmp)") + end +end + +"precompute branch admittance real/imag blocks and node indices" +function _branch_blocks(lm::LiteralModel) + bd = Dict{Int,Any}() + for (i, br) in lm.math["branch"] + get(br, "br_status", 1) == 0 && continue + fc = br["f_connections"]; tc = br["t_connections"] + Z = Matrix{ComplexF64}(br["br_r"] .+ im .* br["br_x"]); Ybr = inv(Z) + m = length(fc) + Ysh = haskey(br, "g_fr") ? Matrix{ComplexF64}(br["g_fr"] .+ im .* br["b_fr"]) : zeros(ComplexF64, m, m) + Yff = Ybr .+ Ysh; Yft = .-Ybr + F = [lm.node_index[(br["f_bus"], c)] for c in fc] + T = [lm.node_index[(br["t_bus"], c)] for c in tc] + bd[br["index"]] = (Gff = real.(Yff), Bff = imag.(Yff), Gft = real.(Yft), Bft = imag.(Yft), F = F, T = T, fc = fc) + end + return bd +end + +"from-side branch current (cr,ci) for connection position `a` (1-based into fc)" +function _branch_current(bd, ctx::HCtx, a::Int) + Gff, Bff, Gft, Bft, F, T = bd.Gff, bd.Bff, bd.Gft, bd.Bft, bd.F, bd.T + cr = zero(eltype(ctx.vr)); ci = zero(eltype(ctx.vr)) + @inbounds for b in eachindex(F) + vrf = ctx.vr[F[b]]; vif = ctx.vi[F[b]]; vrt = ctx.vr[T[b]]; vit = ctx.vi[T[b]] + # I = (Gff+iBff)(vrf+i vif) + (Gft+iBft)(vrt+i vit) + cr += Gff[a, b] * vrf - Bff[a, b] * vif + Gft[a, b] * vrt - Bft[a, b] * vit + ci += Bff[a, b] * vrf + Gff[a, b] * vif + Bft[a, b] * vrt + Gft[a, b] * vit + end + return cr, ci +end + +# ---- model assembly -------------------------------------------------------- + +""" + build_se_model(lm, math; rescaler=1.0, σ_exact=1e-8, zero_inj_sigma=1e-6) + +Assemble the residual model from `math["meas"]`. Current and power *injection* +measurements (`crd,cid,crg,cig,pd,qd,pg,qg,ptot,qtot,cmd,cmg,cad,cag`) are +aggregated to the nodal injection `Ybus*U`; this is exact at the solution and +handles buses hosting several components (e.g. multiple loads). Zero-injection +pseudo-measurements (`z=0`) close the model at source-free terminals for +observability, matching the KCL equalities of `IVRENPowerModel`. +""" +function build_se_model(lm::LiteralModel, math::Dict; rescaler::Float64 = 1.0, + σ_exact::Float64 = 1e-8, zero_inj_sigma::Float64 = 1e-6, + σ_floor::Float64 = 1e-7) + rowfns = Any[]; z = Float64[]; w = Float64[]; rowmeta = NamedTuple[] + brdata = _branch_blocks(lm) + + push_row!(fn, zi, wi, meta) = (push!(rowfns, fn); push!(z, zi); push!(w, wi); push!(rowmeta, meta)) + # floor σ so that a measurement whose value is exactly 0 (e.g. vi at a + # balanced reference bus) does not produce an infinite weight. + wt(σ) = 1.0 / (rescaler * max(σ, σ_floor))^2 + + # nodal current-injection accumulators (real & imag), keyed by node position + inj_r = Dict{Int,Float64}(); inj_i = Dict{Int,Float64}() + var_r = Dict{Int,Float64}(); var_i = Dict{Int,Float64}() + touched = Set{Int}() + # nodes that physically host a gen/load (so are NOT zero-injection) + sourced = Set{Int}() + for (_, ld) in math["load"], c in ld["connections"] + push!(sourced, lm.node_index[(ld["load_bus"], c)]) + end + for (_, g) in math["gen"], c in g["connections"] + push!(sourced, lm.node_index[(g["gen_bus"], c)]) + end + + # accumulate a phase-only current measurement onto nodes (with neutral=-sum) + function accumulate_current!(bus::Int, conns, dst, part::Symbol, sgn::Float64) + nn = _nonneutral(conns) + vals = Float64[]; sgs = Float64[] + for (idx, c) in enumerate(nn) + μ, σ = _zsigma(dst[idx], σ_exact) + push!(vals, μ); push!(sgs, σ) + k = lm.node_index[(bus, c)] + d = part == :r ? inj_r : inj_i; v = part == :r ? var_r : var_i + d[k] = get(d, k, 0.0) + sgn * μ + v[k] = get(v, k, 0.0) + σ^2 + push!(touched, k) + end + if _N_IDX in conns # neutral return = -sum(phases) + k = lm.node_index[(bus, _N_IDX)] + d = part == :r ? inj_r : inj_i; v = part == :r ? var_r : var_i + d[k] = get(d, k, 0.0) + sgn * (-sum(vals)) + v[k] = get(v, k, 0.0) + sum(abs2, sgs) + push!(touched, k) + end + end + + for (m, meas) in math["meas"] + var = meas["var"]; cmp = meas["cmp"]; cmp_id = meas["cmp_id"] + dst = meas["dst"] + active = _active_connections(math, cmp, cmp_id) + nn = _nonneutral(active) + + if var in (:vr, :vi) # native voltage + for (idx, c) in enumerate(nn) + k = lm.node_index[(cmp_id, c)] + μ, σ = _zsigma(dst[idx], σ_exact) + fn = var == :vr ? (ctx -> ctx.vr[k]) : (ctx -> ctx.vi[k]) + push_row!(fn, μ, wt(σ), (m = m, var = var, node = k)) + end + elseif var in (:vm, :vmn) # |U_c - U_n| + kn = lm.node_index[(cmp_id, _N_IDX)] + for (idx, c) in enumerate(nn) + k = lm.node_index[(cmp_id, c)] + μ, σ = _zsigma(dst[idx], σ_exact) + push_row!(ctx -> sqrt((ctx.vr[k]-ctx.vr[kn])^2 + (ctx.vi[k]-ctx.vi[kn])^2), + μ, wt(σ), (m = m, var = var, node = k)) + end + elseif var == :va # ∠(U_c - U_n) + kn = lm.node_index[(cmp_id, _N_IDX)] + for (idx, c) in enumerate(nn) + k = lm.node_index[(cmp_id, c)] + μ, σ = _zsigma(dst[idx], σ_exact) + push_row!(ctx -> atan(ctx.vi[k]-ctx.vi[kn], ctx.vr[k]-ctx.vr[kn]), + μ, wt(σ), (m = m, var = var, node = k)) + end + elseif var == :vll # line-to-line |U_i - U_j| + pairs = length(nn) > 2 ? [(1, 2), (2, 3), (3, 1)] : [(1, 2)] + for (idx, (a, b)) in enumerate(pairs) + ka = lm.node_index[(cmp_id, nn[a])]; kb = lm.node_index[(cmp_id, nn[b])] + μ, σ = _zsigma(dst[idx], σ_exact) + push_row!(ctx -> sqrt((ctx.vr[ka]-ctx.vr[kb])^2 + (ctx.vi[ka]-ctx.vi[kb])^2), + μ, wt(σ), (m = m, var = var)) + end + elseif var in (:cr, :ci, :cm, :ca) # branch current (from-side) + bd = brdata[cmp_id] + for (idx, c) in enumerate(nn) + a = findfirst(isequal(c), bd.fc) + μ, σ = _zsigma(dst[idx], σ_exact) + fn = if var == :cr; ctx -> _branch_current(bd, ctx, a)[1] + elseif var == :ci; ctx -> _branch_current(bd, ctx, a)[2] + elseif var == :cm; ctx -> (cc = _branch_current(bd, ctx, a); sqrt(cc[1]^2 + cc[2]^2)) + else ctx -> (cc = _branch_current(bd, ctx, a); atan(cc[2], cc[1])) end + push_row!(fn, μ, wt(σ), (m = m, var = var)) + end + elseif var in (:p, :q) # branch power flow (EN) + bd = brdata[cmp_id] + fbus = math["branch"][string(cmp_id)]["f_bus"] + kn = lm.node_index[(fbus, _N_IDX)] + for (idx, c) in enumerate(nn) + a = findfirst(isequal(c), bd.fc); k = lm.node_index[(fbus, c)] + μ, σ = _zsigma(dst[idx], σ_exact) + fn = if var == :p + ctx -> (cc = _branch_current(bd, ctx, a); cc[1]*(ctx.vr[k]-ctx.vr[kn]) + cc[2]*(ctx.vi[k]-ctx.vi[kn])) + else + ctx -> (cc = _branch_current(bd, ctx, a); -cc[2]*(ctx.vr[k]-ctx.vr[kn]) + cc[1]*(ctx.vi[k]-ctx.vi[kn])) + end + push_row!(fn, μ, wt(σ), (m = m, var = var)) + end + elseif var in (:crd, :cid) # load current (inj sign -) + accumulate_current!(math["load"][string(cmp_id)]["load_bus"], active, dst, + var == :crd ? :r : :i, -1.0) + elseif var in (:crg, :cig) # gen current (inj sign +) + accumulate_current!(math["gen"][string(cmp_id)]["gen_bus"], active, dst, + var == :crg ? :r : :i, +1.0) + elseif var in (:pd, :qd, :pg, :qg, :ptot, :qtot) # power injection (EN, aggregated) + isgen = var in (:pg, :qg) || (var in (:ptot, :qtot) && cmp == :gen) + bus = isgen ? math["gen"][string(cmp_id)]["gen_bus"] : math["load"][string(cmp_id)]["load_bus"] + s = isgen ? 1.0 : -1.0 + kn = lm.node_index[(bus, _N_IDX)] + ispow_p = var in (:pd, :pg, :ptot) + if var in (:ptot, :qtot) + μ, σ = _zsigma(dst[1], σ_exact) + ks = [lm.node_index[(bus, c)] for c in nn] + fn = ctx -> sum( ispow_p ? + s*(ctx.Ir[k]*(ctx.vr[k]-ctx.vr[kn]) + ctx.Ii[k]*(ctx.vi[k]-ctx.vi[kn])) : + s*(-ctx.Ii[k]*(ctx.vr[k]-ctx.vr[kn]) + ctx.Ir[k]*(ctx.vi[k]-ctx.vi[kn])) for k in ks) + push_row!(fn, μ, wt(σ), (m = m, var = var)) + else + for (idx, c) in enumerate(nn) + k = lm.node_index[(bus, c)]; μ, σ = _zsigma(dst[idx], σ_exact) + fn = ispow_p ? + (ctx -> s*(ctx.Ir[k]*(ctx.vr[k]-ctx.vr[kn]) + ctx.Ii[k]*(ctx.vi[k]-ctx.vi[kn]))) : + (ctx -> s*(-ctx.Ii[k]*(ctx.vr[k]-ctx.vr[kn]) + ctx.Ir[k]*(ctx.vi[k]-ctx.vi[kn]))) + push_row!(fn, μ, wt(σ), (m = m, var = var, node = k)) + end + end + elseif var in (:cmd, :cmg, :cad, :cag) # |I| / ∠I injection (aggregated) + isgen = var in (:cmg, :cag) + bus = isgen ? math["gen"][string(cmp_id)]["gen_bus"] : math["load"][string(cmp_id)]["load_bus"] + ismag = var in (:cmd, :cmg) + for (idx, c) in enumerate(nn) + k = lm.node_index[(bus, c)]; μ, σ = _zsigma(dst[idx], σ_exact) + fn = ismag ? (ctx -> sqrt(ctx.Ir[k]^2 + ctx.Ii[k]^2)) : (ctx -> atan(ctx.Ii[k], ctx.Ir[k])) + push_row!(fn, μ, wt(σ), (m = m, var = var, node = k)) + end + else + error("Literal PMDSE: measurement var :$(var) not supported") + end + end + + # ---- nodal current-injection rows ------------------------------------- + for k in lm.inj_nodes + if k in touched + zr = inj_r[k]; wr = wt(sqrt(max(var_r[k], σ_exact^2))) + push_row!(let kk = k; ctx -> ctx.Ir[kk] end, zr, wr, (var = :inj_r, node = k)) + zi = inj_i[k]; wi = wt(sqrt(max(var_i[k], σ_exact^2))) + push_row!(let kk = k; ctx -> ctx.Ii[kk] end, zi, wi, (var = :inj_i, node = k)) + elseif !(k in sourced) + # genuine zero-injection terminal -> KCL pseudo-measurement (z = 0) + push_row!(let kk = k; ctx -> ctx.Ir[kk] end, 0.0, wt(zero_inj_sigma), (var = :zinj_r, node = k)) + push_row!(let kk = k; ctx -> ctx.Ii[kk] end, 0.0, wt(zero_inj_sigma), (var = :zinj_i, node = k)) + end + # sourced-but-untouched (e.g. power-only nodes): handled by power rows + end + + return SEModel(lm, rowfns, z, w, rowmeta, Any[]) +end + +"evaluate the full residual-model prediction vector at free state `x_free`" +function predict(m::SEModel, x_free::AbstractVector{T}) where {T} + lm = m.lm + 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) + out = Vector{T}(undef, length(m.rowfns)) + @inbounds for i in eachindex(m.rowfns) + out[i] = m.rowfns[i](ctx) + end + return out +end diff --git a/src/bare/mle.jl b/src/bare/mle.jl new file mode 100644 index 0000000..c76c10d --- /dev/null +++ b/src/bare/mle.jl @@ -0,0 +1,69 @@ +################################################################################ +# Literal PMDSE # +# mle.jl : general Maximum-Likelihood interface for the literal estimator. # +# # +# The WLS estimator is the special case of MLE under independent Gaussian # +# measurement errors. This file exposes the likelihood abstraction so that an # +# arbitrary per-measurement distribution (any `Distributions` / # +# `Polynomials`/`ExtendedBeta`, mirroring the `mle` branch of # +# `constraint_mc_residual`) can be plugged in with **no change** to the # +# measurement model `h(x)` or the references. # +# # +# Objective: maximize ℓ(x) = Σ_i logpdf(dst_i, h_i(x)) # +# solved by a (Levenberg-damped) Newton iteration with `ForwardDiff` gradient # +# and Hessian. For `dst_i = Normal(z_i, σ_i)` this is exactly WLS. # +################################################################################ + +"default per-row measurement distributions implied by the WLS weights (Gaussian)" +function gaussian_row_dsts(model::SEModel) + return [_DST.Normal(model.z[i], 1.0 / sqrt(model.w[i])) for i in eachindex(model.z)] +end + +"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 + +""" + solve_mle(model; row_dst=nothing, maxiter=50, tol=1e-9, verbose=false) + +Maximum-likelihood state estimate by **Fisher scoring**: each iteration solves +`(Hᵀ Λ H) Δ = Hᵀ g`, where `g_i = ∂logpdf(dst_i, h_i)/∂h_i` is the per-row score +and `Λ_i = −∂²logpdf/∂h_i²` the per-row information (both via `ForwardDiff`). The +step reuses the WLS QR / orthogonal fallback, so it is robust on ill-conditioned +gains. For `dst_i = Normal(z_i, σ_i)` one has `g_i = (z_i−h_i)/σ_i²`, +`Λ_i = 1/σ_i²`, so the step is exactly the WLS normal-equation step — `solve_mle` +then returns the WLS estimate (verified in the tests). When `row_dst === nothing` +the Gaussian distributions implied by the WLS weights are used; pass `row_dst` +(aligned with the residual rows) for a general likelihood. +""" +function solve_mle(model::SEModel; row_dst = nothing, maxiter::Int = 50, + tol::Float64 = 1e-9, verbose::Bool = false) + t0 = time() + lm = model.lm + dsts = row_dst === nothing ? gaussian_row_dsts(model) : row_dst + score1(i, η) = ForwardDiff.derivative(t -> loglik_row(dsts[i], t), η) + info1(i, η) = -ForwardDiff.derivative(t -> score1(i, t), η) + + x = flat_start_free(lm); term = :maxiter; iters = 0 + for it in 1:maxiter + iters = it + h = predict(model, x) + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + g = [score1(i, h[i]) for i in eachindex(h)] # score + Λ = [max(info1(i, h[i]), 1e-12) for i in eachindex(h)] # information (PSD) + F = Symmetric_full(transpose(H) * (Λ .* H)) + score = transpose(H) * g + A = sqrt.(Λ) .* H; b = g ./ sqrt.(Λ) # A'A=F, A'b=score + Δ, _ = _solve_gain(F, score, A, b) + x .+= Δ + verbose && println(" mle it $it ‖Δ‖∞=$(LinearAlgebra.norm(Δ, Inf))") + LinearAlgebra.norm(Δ, Inf) < tol && (term = :converged; break) + end + + obj = -sum(loglik_row(dsts[i], hi) for (i, hi) in enumerate(predict(model, x))) + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + gcond, grank = _gain_diag(transpose(H) * (model.w .* H)) + x_full = expand_state(lm, x) + return LiteralResult(term, iters, obj, time() - t0, gcond, grank, x, x_full, + state_solution(lm, x_full), :mle) +end diff --git a/src/bare/solve_wlav.jl b/src/bare/solve_wlav.jl new file mode 100644 index 0000000..27a626a --- /dev/null +++ b/src/bare/solve_wlav.jl @@ -0,0 +1,47 @@ +################################################################################ +# Literal PMDSE # +# solve_wlav.jl : Weighted Least Absolute Value via Iteratively Reweighted # +# Least Squares (IRLS). Minimises Σ (1/(rescaler·σ_i)) |r_i|, # +# matching the `wlav` criterion of `constraint_mc_residual`. # +################################################################################ + +""" + solve_wlav(model; maxiter=60, tol=1e-9, huber=1e-4, verbose=false) + +IRLS solution of the WLAV problem. Each outer iteration solves a Gauss-Newton +WLS step with weights `w̃_i = √w_i / max(|r_i|, huber)` so that the quadratic +surrogate `Σ w̃_i r_i²` upper-bounds the WLAV objective `Σ √w_i |r_i|` (Huber +relaxation near `r=0`). WLAV is robust to single gross errors, unlike WLS. +""" +function solve_wlav(model::SEModel; maxiter::Int = 60, tol::Float64 = 1e-9, + huber::Float64 = 1e-4, verbose::Bool = false) + t0 = time() + lm = model.lm + x = flat_start_free(lm) + wbase = model.w; z = model.z + sw = sqrt.(wbase) # √w_i (= 1/(rescaler σ_i)) + term = :maxiter; iters = 0 + + for it in 1:maxiter + iters = it + h = predict(model, x); r = z .- h + w̃ = sw ./ max.(abs.(r), huber) # IRLS weights + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + WH = w̃ .* H + G = Symmetric_full(transpose(H) * WH) + rhs = transpose(H) * (w̃ .* r) + A = sqrt.(w̃) .* H; b = sqrt.(w̃) .* r + Δ, _ = _solve_gain(G, rhs, A, b) + x .+= Δ + verbose && println(" wlav it $it ‖Δ‖∞=$(LinearAlgebra.norm(Δ, Inf))") + LinearAlgebra.norm(Δ, Inf) < tol && (term = :converged; break) + end + + h = predict(model, x); r = z .- h + obj = sum(sw .* abs.(r)) # WLAV objective Σ √w_i |r_i| + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + gcond, grank = _gain_diag(transpose(H) * (wbase .* H)) + x_full = expand_state(lm, x) + sol = state_solution(lm, x_full) + return LiteralResult(term, iters, obj, time() - t0, gcond, grank, x, x_full, sol, :wlav) +end diff --git a/src/bare/solve_wls.jl b/src/bare/solve_wls.jl new file mode 100644 index 0000000..8d51c8d --- /dev/null +++ b/src/bare/solve_wls.jl @@ -0,0 +1,98 @@ +################################################################################ +# Literal PMDSE # +# solve_wls.jl : Gauss-Newton Weighted Least Squares (normal equations with a # +# QR / orthogonal fallback, Abur & Exposito Ch. 2-3). # +################################################################################ + +"result of a literal state-estimation solve" +struct LiteralResult + termination::Symbol + iterations::Int + objective::Float64 # weighted SSR r' W r at the solution + solve_time::Float64 + gain_cond::Float64 + gain_rank::Int + x_free::Vector{Float64} + x_full::Vector{Float64} + solution::Dict{String,Any} + estimator::Symbol +end + +""" + _solve_gain(G, rhs, A, b; cond_tol) + +Solve the normal-equation step `G Δ = rhs` with `G = HᵀWH`. Tries a Cholesky +factorization first; if `G` is not positive-definite or ill-conditioned, falls +back to the numerically robust orthogonal solve `min ‖A Δ − b‖₂` with +`A = √W·H`, `b = √W·r` (a QR least-squares step). Returns `(Δ, cond, used_qr)`. +""" +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 + +"condition number / numerical rank diagnostics, guarded against LAPACK failures" +function _gain_diag(G::Matrix{Float64}) + try + sv = LinearAlgebra.svdvals(G) + smax = sv[1]; smin = sv[end] + c = smin > 0 ? smax / smin : Inf + r = sum(sv .> 1e-9 * smax) + return c, r + catch + return NaN, -1 + end +end + +""" + solve_wls(model; maxiter=50, tol=1e-9, verbose=false, cond_tol=1e12) + +Iterate the Gauss-Newton WLS step `(HᵀWH)Δ = HᵀW r`, `r = z − h(x)`, `H = +∂h/∂x` (via `ForwardDiff`), until `‖Δ‖∞ < tol`. +""" +function solve_wls(model::SEModel; maxiter::Int = 50, tol::Float64 = 1e-9, + verbose::Bool = false, cond_tol::Float64 = 1e12) + t0 = time() + lm = model.lm + x = flat_start_free(lm) + W = model.w; z = model.z + sqrtW = sqrt.(W) + term = :maxiter; iters = 0; gcond = NaN; grank = -1 + + for it in 1:maxiter + iters = it + h = predict(model, x) + r = z .- h + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + WH = W .* H + G = Symmetric_full(transpose(H) * WH) + rhs = transpose(H) * (W .* r) + A = sqrtW .* H; b = sqrtW .* r + Δ, used_qr = _solve_gain(G, rhs, A, b; cond_tol = cond_tol) + x .+= Δ + verbose && println(" it $it ‖Δ‖∞=$(LinearAlgebra.norm(Δ, Inf)) qr=$(used_qr)") + if LinearAlgebra.norm(Δ, Inf) < tol + term = :converged; break + end + end + + # diagnostics + objective at the solution + h = predict(model, x); r = z .- h + obj = sum(W .* r .^ 2) + H = ForwardDiff.jacobian(xx -> predict(model, xx), x) + G = transpose(H) * (W .* H) + gcond, grank = _gain_diag(G) + x_full = expand_state(lm, x) + sol = state_solution(lm, x_full) + return LiteralResult(term, iters, obj, time() - t0, gcond, grank, x, x_full, sol, :wls) +end + +"materialize a symmetric matrix (full, dense) from `HᵀWH`" +Symmetric_full(M) = Matrix(LinearAlgebra.Symmetric((M .+ transpose(M)) ./ 2)) diff --git a/src/core/export.jl b/src/core/export.jl index 6e50b09..4af0481 100644 --- a/src/core/export.jl +++ b/src/core/export.jl @@ -12,6 +12,10 @@ export calculate_voltage_magnitude_error export update_load_bounds!, update_voltage_bounds!, update_generator_bounds!, update_all_bounds! export ExtendedBeta +# Literal PMDSE (explicit matrix-based, JuMP-free, explicit-neutral estimator) +export solve_mc_se_literal, LiteralModel, LiteralResult, build_se_model +export solve_wls, solve_wlav, solve_mle, accuracy_metrics, voltage_errors + # so that users do not need to import JuMP to use a solver with PowerModelsDistribution import JuMP: optimizer_with_attributes export optimizer_with_attributes diff --git a/test/literal/helpers.jl b/test/literal/helpers.jl new file mode 100644 index 0000000..ba50d23 --- /dev/null +++ b/test/literal/helpers.jl @@ -0,0 +1,199 @@ +################################################################################ +# Shared helpers for the Literal PMDSE test suite. # +# (Assumes the including scope has imported `_PMD`, `_PMDSE`, `_DST`, # +# `ForwardDiff`, `LinearAlgebra`, `Ipopt` and `Test`.) # +################################################################################ + +const LIT_N = _PMDSE._N_IDX # neutral conductor index (4) + +"main.tex toy: 3-bus 4-wire feeder with realistic per-unit impedance" +function toy_math() + Zb = (400.0^2) / 100_000.0 + zp = (0.01 + 0.02im) / Zb + zm = (0.001 + 0.002im) / Zb + zn = (0.05 + 0.01im) / Zb + zpn = (0.0005 + 0.001im) / Zb + Z = Matrix{ComplexF64}(undef, 4, 4) + for i in 1:3, j in 1:3; Z[i, j] = i == j ? zp : zm; end + Z[4, 4] = zn + for i in 1:3; Z[i, 4] = zpn; Z[4, i] = zpn; end + br_r = real.(Z); br_x = imag.(Z) + va = deg2rad.([0.0, -120.0, 120.0]) + vr0 = [cos.(va)..., 0.0]; vi0 = [sin.(va)..., 0.0] + bus(idx, bt, gnd) = Dict{String,Any}("index"=>idx, "bus_i"=>idx, "bus_type"=>bt, + "terminals"=>[1,2,3,4], "grounded"=>Bool[false,false,false,gnd], + "vbase"=>0.4/sqrt(3), "vmin"=>zeros(4), "vmax"=>fill(Inf,4), + "vr_start"=>copy(vr0), "vi_start"=>copy(vi0)) + branch(idx, f, t) = Dict{String,Any}("index"=>idx, "f_bus"=>f, "t_bus"=>t, + "f_connections"=>[1,2,3,4], "t_connections"=>[1,2,3,4], + "br_r"=>br_r, "br_x"=>br_x, "br_status"=>1) + load(idx, b, conns, pd, qd) = Dict{String,Any}("index"=>idx, "load_bus"=>b, + "connections"=>conns, "configuration"=>_PMD.WYE, "pd"=>pd, "qd"=>qd, + "model"=>_PMD.POWER, "status"=>1) + Dict{String,Any}( + "bus"=>Dict("1"=>bus(1,3,true), "2"=>bus(2,1,false), "3"=>bus(3,1,false)), + "branch"=>Dict("1"=>branch(1,1,2), "2"=>branch(2,2,3)), + "load"=>Dict("1"=>load(1,2,[1,2,3,4],[0.3,0.4,0.5],[0.1,0.1,0.2]), + "2"=>load(2,3,[1,4],[0.4],[0.1])), + "gen"=>Dict("1"=>Dict{String,Any}("index"=>1,"gen_bus"=>1, + "connections"=>[1,2,3,4],"configuration"=>_PMD.WYE)), + "shunt"=>Dict{String,Any}(), "settings"=>Dict{String,Any}("sbase"=>1.0)) +end + +"balanced reference-bus phasor for a math dict, keyed by terminal" +function balanced_ref(math) + rb = _PMDSE.literal_ref_bus(math); va = deg2rad.([0.0,-120.0,120.0]) + vr = Dict{Int,Float64}(); vi = Dict{Int,Float64}() + for t in math["bus"][string(rb)]["terminals"] + vr[t] = t==LIT_N ? 0.0 : cos(va[t]); vi[t] = t==LIT_N ? 0.0 : sin(va[t]) + end + (vr, vi) +end + +"constant-power Newton PF for the toy; returns the free-state vector of `lm`" +function toy_pf_state(lm) + math = lm.math; inj = lm.inj_nodes + function comp_inj(vr, vi) + Ic = Complex.(zeros(eltype(vr), lm.n)) + for (_, ld) in math["load"] + b = ld["load_bus"]; conns = ld["connections"]; ph = _PMDSE._nonneutral(conns) + kn = LIT_N in conns ? lm.node_index[(b,LIT_N)] : 0 + for (idx, c) in enumerate(ph) + k = lm.node_index[(b,c)] + Up = vr[k]+im*vi[k]; Un = kn==0 ? 0.0+0im : vr[kn]+im*vi[kn] + Il = conj((ld["pd"][idx]+im*ld["qd"][idx])/(Up-Un)) + Ic[k] -= Il; kn != 0 && (Ic[kn] += Il) + end + end + Ic + end + function resid(xf) + x = _PMDSE.expand_state(lm, xf); vr,vi = _PMDSE.vrvi(lm,x) + Ir,Ii = _PMDSE.nodal_injection(lm, collect(vr), collect(vi)); Ic = comp_inj(collect(vr), collect(vi)) + out = similar(xf, 2*length(inj)) + for (j,k) in enumerate(inj); out[2j-1]=Ir[k]-real(Ic[k]); out[2j]=Ii[k]-imag(Ic[k]); end + out + end + xf = _PMDSE.flat_start_free(lm) + for _ in 1:60 + F = resid(xf); J = ForwardDiff.jacobian(resid, xf); Δ = J \ F; xf .-= Δ + LinearAlgebra.norm(Δ, Inf) < 1e-13 && break + end + xf +end + +"toy ground-truth PF as a PMD-style solution dict (+ per-load pd/qd)" +function toy_pf() + math = toy_math() + lmf = _PMDSE.LiteralModel(math; reference=:full_slack, ref_values=balanced_ref(math)) + xf = toy_pf_state(lmf); sol = _PMDSE.state_solution(lmf, _PMDSE.expand_state(lmf, xf)) + for (b,bsol) in sol["bus"]; bsol["terminals"] = math["bus"][b]["terminals"]; end + pf = Dict("solution"=>Dict("bus"=>sol["bus"], "load"=>Dict{String,Any}())) + for (l,ld) in math["load"]; pf["solution"]["load"][l]=Dict("pd"=>ld["pd"],"qd"=>ld["qd"]); end + math, pf +end + +"parse + IVREN power-flow an explicit-neutral DSS feeder shipped with the package" +function load_en_feeder(name) + eng = _PMD.parse_file(joinpath(_PMDSE.BASE_DIR,"test","data","three-bus-en-models",name)) + _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) + slv = _PMDSE.optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>0, "tol"=>1e-10) + pf = _PMD.solve_mc_opf(math, _PMD.IVRENPowerModel, slv) + for (b,bsol) in pf["solution"]["bus"]; bsol["terminals"]=math["bus"][b]["terminals"]; end + math, pf +end + +"write + add the canonical IVREN measurement set (vr/vi + crd/cid + crg/cig)" +function canonical_meas!(math, pf; σ=0.005) + path = joinpath(mktempdir(),"m.csv") + _PMDSE.write_measurements!(_PMD.IVRENPowerModel, math, pf, path, σ=σ) + _PMDSE.add_measurements!(math, path, actual_meas=true) + math +end + +"reference-bus true phasor from a PF dict, as (vr_dict, vi_dict) over terminals" +function ref_values_from_pf(math, pf) + rb = _PMDSE.literal_ref_bus(math) + bus = math["bus"][string(rb)]; bsol = pf["solution"]["bus"][string(rb)] + vr = Dict{Int,Float64}(); vi = Dict{Int,Float64}() + for (idx,t) in enumerate(bus["terminals"]); vr[t]=bsol["vr"][idx]; vi[t]=bsol["vi"][idx]; end + (vr, vi) +end + +"complex node voltages from a PF dict, aligned with lm.nodes" +function pf_voltages(lm, pf) + U = zeros(ComplexF64, lm.n) + for (k,(b,t)) in enumerate(lm.nodes) + bsol = pf["solution"]["bus"][string(b)]; idx=findfirst(isequal(t), lm.math["bus"][string(b)]["terminals"]) + U[k] = bsol["vr"][idx] + im*bsol["vi"][idx] + end + U +end + +"true free-state for `lm` from a PF dict" +function pf_xfree(lm, pf) + x = copy(lm.x_fixed) + for (k,(b,t)) in enumerate(lm.nodes) + bsol = pf["solution"]["bus"][string(b)]; idx=findfirst(isequal(t), lm.math["bus"][string(b)]["terminals"]) + x[k]=bsol["vr"][idx]; x[lm.n+k]=bsol["vi"][idx] + end + x[lm.free_idx] +end + +"(gen − load) current injection per node, from PF component currents" +function pf_component_injection(lm, pf) + enI(conns, cr, ci) = (np=length(cr); I=Complex{Float64}[cr[i]+im*ci[i] for i in 1:np]; + length(conns)>np && push!(I, -sum(I)); I) + Icmp = zeros(ComplexF64, lm.n) + for (g,gen) in lm.math["gen"] + Ig = enI(gen["connections"], pf["solution"]["gen"][g]["crg"], pf["solution"]["gen"][g]["cig"]) + for (idx,c) in enumerate(gen["connections"]); Icmp[lm.node_index[(gen["gen_bus"],c)]] += Ig[idx]; end + end + for (l,load) in lm.math["load"] + Il = enI(load["connections"], pf["solution"]["load"][l]["crd"], pf["solution"]["load"][l]["cid"]) + for (idx,c) in enumerate(load["connections"]); Icmp[lm.node_index[(load["load_bus"],c)]] -= Il[idx]; end + end + Icmp +end + +"a measurement dict exercising every supported var type (for the Jacobian test)" +function kitchen_sink_meas() + n3() = [_DST.Normal(0.3*j, 0.01) for j in 1:3]; n1() = [_DST.Normal(0.2, 0.01)] + M = Dict{String,Any}(); i = 0 + add(cmp, cid, var, dst) = (M["$(i+=1)"]=Dict{String,Any}("cmp"=>cmp,"cmp_id"=>cid,"var"=>var,"dst"=>dst,"crit"=>"wls")) + for v in (:vr,:vi,:vm,:va); add(:bus,3,v,n3()); end + add(:bus,3,:vll,n3()) + for v in (:cr,:ci,:cm,:ca,:p,:q); add(:branch,1,v,n3()); end + for v in (:crd,:cid,:cmd,:cad,:pd,:qd); add(:load,1,v,n3()); end + for v in (:crg,:cig,:cmg,:cag,:pg,:qg); add(:gen,1,v,n3()); end + add(:load,1,:ptot,n1()); add(:load,1,:qtot,n1()) + M +end + +"rotation-invariant measurement set (vm + pd + qd) for the observability test" +function rotation_invariant_meas(math, pf) + M = Dict{String,Any}(); i = 0 + add(cmp, cid, var, dst) = (M["$(i+=1)"]=Dict{String,Any}("cmp"=>cmp,"cmp_id"=>cid,"var"=>var,"dst"=>dst,"crit"=>"wls")) + for (b,bus) in math["bus"] + bsol = pf["solution"]["bus"][b]; nt=length(bus["terminals"]) + vmn = [sqrt((bsol["vr"][k]-bsol["vr"][end])^2+(bsol["vi"][k]-bsol["vi"][end])^2) for k in 1:nt-1] + add(:bus, bus["index"], :vm, [_DST.Normal(v,0.01) for v in vmn]) + end + for (l,load) in math["load"] + lsol = pf["solution"]["load"][l]; np=length(lsol["pd"]) + add(:load, load["index"], :pd, [_DST.Normal(lsol["pd"][k],0.01) for k in 1:np]) + add(:load, load["index"], :qd, [_DST.Normal(lsol["qd"][k],0.01) for k in 1:np]) + end + M +end + +"rotation direction δx=[-vi; vr] restricted to free coordinates" +function rotation_dir(lm, xfull) + vr, vi = _PMDSE.vrvi(lm, xfull) + d = zeros(2*lm.n); d[1:lm.n] = -vi; d[lm.n+1:2lm.n] = vr + d[lm.free_idx] +end + +with_meas(math, M) = (m = deepcopy(math); m["meas"] = M; m) diff --git a/test/literal/test_benchmark.jl b/test/literal/test_benchmark.jl new file mode 100644 index 0000000..df0776f --- /dev/null +++ b/test/literal/test_benchmark.jl @@ -0,0 +1,68 @@ +################################################################################ +# Gate 5 : benchmark the JuMP IVREN baseline against the literal estimators # +# (WLS & WLAV) over several noise seeds. Headline metrics: estimation accuracy # +# vs ground truth (RMSE / max per-phase |U| incl. neutral) and compute # +# time / iterations. Results are *compared*, not asserted equal. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +import Statistics +include("helpers.jl") + +_sig(x; d=2) = round(x; sigdigits=d) + +function ivren_maxerr(SE, pf) + e = 0.0 + for (b, bsol) in SE["solution"]["bus"], idx in 1:length(bsol["vr"]) + e = max(e, abs((bsol["vr"][idx]+im*bsol["vi"][idx]) - + (pf["solution"]["bus"][b]["vr"][idx]+im*pf["solution"]["bus"][b]["vi"][idx]))) + end + e +end + +@testset "Literal PMDSE — IVREN vs literal benchmark (gate 5)" begin + math0, pf = load_en_feeder("3bus_4wire.dss") + refvals = ref_values_from_pf(math0, pf) + σ = 0.01; seeds = [1, 2, 3, 11] + slv = _PMDSE.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-9, "print_level"=>0) + + rows = NamedTuple[] + for s in seeds + math, _ = load_en_feeder("3bus_4wire.dss") + path = joinpath(mktempdir(), "m$s.csv") + _PMDSE.write_measurements!(_PMD.IVRENPowerModel, math, pf, path, σ=σ) + _PMDSE.add_measurements!(math, path, actual_meas=false, seed=s) + math["se_settings"] = Dict{String,Any}("criterion"=>"wls", "rescaler"=>1.0) + + SE = _PMDSE.solve_ivr_en_mc_se(math, slv) + wls = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:full_slack, ref_values=refvals) + wlav = _PMDSE.solve_mc_se_literal(math; estimator=:wlav, reference=:full_slack, ref_values=refvals) + + wls_err = _PMDSE.accuracy_metrics(wls, pf["solution"]).maxerr + wlav_err = _PMDSE.accuracy_metrics(wlav, pf["solution"]).maxerr + push!(rows, (seed=s, + ivren=(ivren_maxerr(SE, pf), get(SE, "solve_time", NaN)), + wls =(wls_err, wls.solve_time, wls.iterations), + wlav=(wlav_err, wlav.solve_time, wlav.iterations))) + + @test wls.termination == :converged + @test wlav.termination == :converged + @test SE["termination_status"] ∈ [_PMDSE.LOCALLY_SOLVED, _PMDSE.ALMOST_LOCALLY_SOLVED] + # literal accuracy is comparable to the JuMP baseline (not asserted equal) + @test wls_err < 10σ + end + + println("\n IVREN vs Literal PMDSE — 3-bus 4-wire, σ=$(σ)") + println(" seed | IVREN max|U| t[s] | WLS max|U| t[s] it | WLAV max|U| t[s] it") + for r in rows + println(" $(r.seed) | ", _sig(r.ivren[1]), " ", _sig(r.ivren[2]), + " | ", _sig(r.wls[1]), " ", _sig(r.wls[2]; d=3), " ", r.wls[3], + " | ", _sig(r.wlav[1]), " ", _sig(r.wlav[2]; d=3), " ", r.wlav[3]) + end + avg(f) = Statistics.mean(f(r) for r in rows) + println(" mean | IVREN ", _sig(avg(r->r.ivren[1])), + " | WLS ", _sig(avg(r->r.wls[1])), " | WLAV ", _sig(avg(r->r.wlav[1]))) +end diff --git a/test/literal/test_hx_jacobian.jl b/test/literal/test_hx_jacobian.jl new file mode 100644 index 0000000..b62ad32 --- /dev/null +++ b/test/literal/test_hx_jacobian.jl @@ -0,0 +1,38 @@ +################################################################################ +# Gate 2 : the ForwardDiff Jacobian H = ∂h/∂x matches finite differences for # +# every supported measurement var-type. Run on the well-scaled toy feeder so # +# that the steep magnitude/angle rows are not dominated by FD truncation; a # +# relative tolerance is used because FD truncation error scales with |∂h/∂x|. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +include("helpers.jl") + +@testset "Literal PMDSE — h(x) Jacobian vs finite differences (gate 2)" begin + math = toy_math() + lm = _PMDSE.LiteralModel(math; reference=:sota) + model = _PMDSE.build_se_model(lm, with_meas(math, kitchen_sink_meas())) + + # measurement-type coverage check + vars = Set(getfield(r, :var) for r in model.rowmeta if haskey(r, :var)) + for v in (:vr,:vi,:vm,:va,:vll,:cr,:ci,:cm,:ca,:p,:q,:inj_r,:inj_i) + @test v in vars + end + + rng_states = [0.05, -0.03, 0.07, 0.02, -0.06] + maxrel = 0.0 + for s in rng_states + x = _PMDSE.flat_start_free(lm) .+ s + Jad = ForwardDiff.jacobian(xx -> _PMDSE.predict(model, xx), x) + δ = 1e-6; Jfd = similar(Jad) + for j in 1:length(x) + xp = copy(x); xp[j]+=δ; xm = copy(x); xm[j]-=δ + Jfd[:, j] = (_PMDSE.predict(model, xp) .- _PMDSE.predict(model, xm)) ./ (2δ) + end + maxrel = max(maxrel, maximum(abs.(Jad .- Jfd) ./ (1 .+ abs.(Jad)))) + end + @test maxrel < 1e-6 +end diff --git a/test/literal/test_references.jl b/test/literal/test_references.jl new file mode 100644 index 0000000..6c62686 --- /dev/null +++ b/test/literal/test_references.jl @@ -0,0 +1,33 @@ +################################################################################ +# Gate 3 : reproduce the observability theorem of main.tex. # +# With only the neutral grounded (:prop) the global phase rotation # +# δx = [−vi; vr] is a null vector of H (power / |U| measurements are invariant # +# to a global rotation) ⇒ the gain is singular. Adding one angle datum # +# (:sota) removes that null direction ⇒ the system becomes observable. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +include("helpers.jl") + +@testset "Literal PMDSE — angular reference / observability (gate 3)" begin + math, pf = toy_pf() + M = rotation_invariant_meas(math, pf) + + rot_resid = Dict{Symbol,Float64}() + for ref in (:prop, :sota) + lm = _PMDSE.LiteralModel(math; reference=ref) + model = _PMDSE.build_se_model(lm, with_meas(math, M)) + xf = pf_xfree(lm, pf) + H = ForwardDiff.jacobian(xx -> _PMDSE.predict(model, xx), xf) + d = rotation_dir(lm, _PMDSE.expand_state(lm, xf)) + rot_resid[ref] = LinearAlgebra.norm(H * d) / (LinearAlgebra.norm(H) * LinearAlgebra.norm(d)) + end + + # without the angle datum the rotation is (numerically) unobservable … + @test rot_resid[:prop] < 1e-8 + # … and with it the rotation becomes observable + @test rot_resid[:sota] > 1e-3 +end diff --git a/test/literal/test_wls_vs_ivren.jl b/test/literal/test_wls_vs_ivren.jl new file mode 100644 index 0000000..ef0b58a --- /dev/null +++ b/test/literal/test_wls_vs_ivren.jl @@ -0,0 +1,60 @@ +################################################################################ +# Gate 4 : on a noiseless, consistent measurement set with identical references, # +# the literal WLS converges to the same state as `solve_ivr_en_mc_se` # +# (IVRENPowerModel). This pins down every sign / convention. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +include("helpers.jl") + +@testset "Literal PMDSE — WLS vs IVREN (gate 4)" begin + math, pf = load_en_feeder("3bus_4wire.dss") + canonical_meas!(math, pf; σ=0.005) # noiseless: values == PF + math["se_settings"] = Dict{String,Any}("criterion"=>"wls", "rescaler"=>1.0) + refvals = ref_values_from_pf(math, pf) + + @testset "noiseless literal WLS reproduces the PF state" begin + for ref in (:full_slack, :sota) + rv = ref == :full_slack ? refvals : nothing + res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=ref, ref_values=rv) + @test res.termination == :converged + @test _PMDSE.accuracy_metrics(res, pf["solution"]).maxerr < 1e-6 + end + end + + @testset "Gaussian MLE reduces to WLS" begin + wls = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:sota) + mle = _PMDSE.solve_mc_se_literal(math; estimator=:mle, reference=:sota) + @test mle.termination == :converged + @test maximum(abs.(mle.x_free .- wls.x_free)) < 1e-8 + end + + @testset "literal WLS matches the JuMP IVREN estimator" begin + slv = _PMDSE.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-10, "print_level"=>0) + SE = _PMDSE.solve_ivr_en_mc_se(math, slv) + res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:full_slack, ref_values=refvals) + maxdiff = 0.0 + for (b, bsol) in SE["solution"]["bus"] + for idx in 1:length(bsol["vr"]) + maxdiff = max(maxdiff, abs((res.solution["bus"][b]["vr"][idx] + im*res.solution["bus"][b]["vi"][idx]) - + (bsol["vr"][idx] + im*bsol["vi"][idx]))) + end + end + @test maxdiff < 1e-4 # both recover the same physical state + end + + @testset "noisy literal WLS stays close to ground truth" begin + mathn, _ = load_en_feeder("3bus_4wire.dss") + path = joinpath(mktempdir(), "mn.csv") + _PMDSE.write_measurements!(_PMD.IVRENPowerModel, mathn, pf, path, σ=0.01) + _PMDSE.add_measurements!(mathn, path, actual_meas=false, seed=7) # sampled noise + mathn["se_settings"] = Dict{String,Any}("criterion"=>"wls", "rescaler"=>1.0) + res = _PMDSE.solve_mc_se_literal(mathn; estimator=:wls, reference=:full_slack, + ref_values=ref_values_from_pf(mathn, pf)) + @test res.termination == :converged + @test _PMDSE.accuracy_metrics(res, pf["solution"]).maxerr < 0.05 + end +end diff --git a/test/literal/test_ybus.jl b/test/literal/test_ybus.jl new file mode 100644 index 0000000..032f77e --- /dev/null +++ b/test/literal/test_ybus.jl @@ -0,0 +1,34 @@ +################################################################################ +# Gate 1 : the assembled Ybus reproduces the network physics. # +# PMD's `calc_admittance_matrix` does not support the 4-wire explicit-neutral # +# (kron_reduce=false) model, so we validate against PMD's own power-flow: # +# `Ybus·U` must equal the component current injection (Σgen − Σload) at every # +# injection node, which is exactly the KCL that `IVRENPowerModel` enforces. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +include("helpers.jl") + +@testset "Literal PMDSE — Ybus assembly (gate 1)" begin + math, pf = load_en_feeder("3bus_4wire.dss") + lm = _PMDSE.LiteralModel(math; reference=:sota) + + @test LinearAlgebra.norm(lm.G .- transpose(lm.G), Inf) < 1e-8 # Ybus symmetric + @test LinearAlgebra.norm(lm.B .- transpose(lm.B), Inf) < 1e-8 + + U = pf_voltages(lm, pf); vr = real.(U); vi = imag.(U) + Ir, Ii = _PMDSE.nodal_injection(lm, vr, vi) + Icalc = Ir .+ im .* Ii + Icmp = pf_component_injection(lm, pf) + + maxerr = maximum(abs(Icalc[k] - Icmp[k]) for k in lm.inj_nodes) + @test maxerr < 1e-6 # Ybus·U == component injection (Σgen − Σload) + + # toy network: well-scaled Ybus must also be symmetric and finite + lmt = _PMDSE.LiteralModel(toy_math(); reference=:sota) + @test all(isfinite, lmt.G) && all(isfinite, lmt.B) + @test LinearAlgebra.norm(lmt.G .- transpose(lmt.G), Inf) < 1e-10 +end diff --git a/test/runtests.jl b/test/runtests.jl index 9d00b72..855ba7f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ import PowerModels import PowerModelsDistribution as _PMD import Statistics using Test +using SafeTestsets #network and feeder from ENWL for tests ntw, fdr = 4, 2 @@ -51,6 +52,13 @@ ipopt_solver = _PMDSE.optimizer_with_attributes(Ipopt.Optimizer,"max_cpu_time" = include("with_errors.jl") end +# Literal PMDSE: explicit matrix-based explicit-neutral estimator (JuMP-free) +@safetestset "Literal PMDSE — Ybus assembly" begin include("literal/test_ybus.jl") end +@safetestset "Literal PMDSE — h(x) Jacobian" begin include("literal/test_hx_jacobian.jl") end +@safetestset "Literal PMDSE — references" begin include("literal/test_references.jl") end +@safetestset "Literal PMDSE — WLS vs IVREN" begin include("literal/test_wls_vs_ivren.jl") end +@safetestset "Literal PMDSE — benchmark" begin include("literal/test_benchmark.jl") end + ambiguities = Test.detect_ambiguities(_PMDSE); if !isempty(ambiguities) println("ambiguities detected: $ambiguities")