From 685968a59919d97ea5188ddcef82b7d5b3924865 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 08:40:15 +0000 Subject: [PATCH 1/2] Literal PMDSE: add PowerGridModel WLS solve options, validated against PGM MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implement the two PowerGridModel (PGM) state-estimation solve options on the literal (JuMP-free) estimator and make them general — the same four-wire, explicit-neutral code path that handles the EN feeders, with the single-phase PGM networks as the one-conductor / no-neutral special case: - src/bare/pgm_se.jl - build_se_atoms: turn data["meas"] (the package measurement dict) into a flat list of measurement atoms over the nodal rectangular state x=[vr;vi], with phase-to-neutral voltages/powers, neutral injection = -sum(phase currents), appliance-injection aggregation, Kalman-combined duplicate sensors, and zero-injection KCL pseudo-measurements. - solve_se_il (:iterative_linear, PGM default): re-linearise measurements to complex currents/voltage-phasors each iteration, re-solve a constant linear WLS; slack-angle gauge added automatically when no angle is measured. - solve_se_nr (:newton_raphson): Gauss-Newton on the nonlinear WLS, warm -started from one iterative_linear solve (flat start is degenerate for power-/magnitude-only systems). - solve_mc_se_literal gains a `method` kwarg selecting the PGM solve option; default behaviour (in-house Gauss-Newton WLS) is unchanged. Validation (TDD): - test/literal/test_pgm_benchmark.jl (gate 6): reproduce PowerGridModel's own state-estimation examples (single-node, 1os2msr +/- angle, single-line) — node voltages and line flows match PGM to ~1e-7 for both methods; the two methods agree. Golden data generated from PGM in test/literal/pgm/ (generate_golden.py). - test/literal/test_pgm_en_generality.jl (gate 7): the same solve options run unchanged on the 3-bus 4-wire explicit-neutral feeder, recover the IVREN power flow (neutral included) and agree with solve_wls. Per-unit bases match PGM (V_base = u_rated, S_base = base_power_3p = 1e6). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DnxB7KtuKVUYCK4fAQyahK --- CHANGELOG.md | 9 + docs/src/Literal_PMDSE_PLAN.md | 32 + docs/src/literal_pmdse.md | 22 +- src/PowerModelsDistributionStateEstimation.jl | 1 + src/bare/literal_pmdse.jl | 30 +- src/bare/pgm_se.jl | 603 ++++++++++++++++++ src/core/export.jl | 2 + test/literal/pgm/README.md | 34 + test/literal/pgm/generate_golden.py | 140 ++++ test/literal/pgm/golden.json | 538 ++++++++++++++++ test/literal/pgm/pgm_cases_data.jl | 65 ++ test/literal/pgm_cases.jl | 195 ++++++ test/literal/test_pgm_benchmark.jl | 93 +++ test/literal/test_pgm_en_generality.jl | 47 ++ test/runtests.jl | 2 + 15 files changed, 1807 insertions(+), 6 deletions(-) create mode 100644 src/bare/pgm_se.jl create mode 100644 test/literal/pgm/README.md create mode 100644 test/literal/pgm/generate_golden.py create mode 100644 test/literal/pgm/golden.json create mode 100644 test/literal/pgm/pgm_cases_data.jl create mode 100644 test/literal/pgm_cases.jl create mode 100644 test/literal/test_pgm_benchmark.jl create mode 100644 test/literal/test_pgm_en_generality.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index fec4b76..21d1e35 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,12 @@ +## unreleased +- Literal PMDSE: add the two PowerGridModel WLS solve options + (`solve_mc_se_literal(...; estimator=:wls, method=:iterative_linear|:newton_raphson)`), + generalised to the four-wire explicit-neutral model (`src/bare/pgm_se.jl`). + - validated against PowerGridModel's own state-estimation examples + (`test/literal/test_pgm_benchmark.jl`, golden data in `test/literal/pgm/`) + - shown to run unchanged on the explicit-neutral feeder + (`test/literal/test_pgm_en_generality.jl`) + ## v0.8.0 - @MohamedNumair adds: - explicit nuetral modelling diff --git a/docs/src/Literal_PMDSE_PLAN.md b/docs/src/Literal_PMDSE_PLAN.md index 15b754f..42fe7f2 100644 --- a/docs/src/Literal_PMDSE_PLAN.md +++ b/docs/src/Literal_PMDSE_PLAN.md @@ -16,6 +16,36 @@ res = solve_mc_se_literal(data_math; estimator = :wlav, reference = :full_slack, res = solve_mc_se_literal(data_math; estimator = :mle) # Gaussian MLE == WLS ``` +## PowerGridModel solve options + +In addition to the in-house Gauss–Newton WLS, the two +[PowerGridModel](https://github.com/PowerGridModel/power-grid-model) (PGM) +state-estimation *solve options* are reproduced and made **general** (the same +four-wire, explicit-neutral code path, not PGM's Kron-reduced single phase): + +```julia +res = solve_mc_se_literal(data_math; estimator = :wls, method = :iterative_linear) # PGM default +res = solve_mc_se_literal(data_math; estimator = :wls, method = :newton_raphson) +``` + +* `:iterative_linear` ([`solve_se_il`](@ref)) — every measurement is linearised + to a complex current / voltage-phasor at the previous voltages; the (constant) + linear WLS matrix is re-solved each iteration. Re-linearised RHS only, so it is + cheap; the slack-angle gauge `Im(U_ref)=0` is added automatically when no + voltage angle/phasor is measured. +* `:newton_raphson` ([`solve_se_nr`](@ref)) — Gauss–Newton on the nonlinear WLS, + warm-started from one `iterative_linear` solve (the flat start is degenerate for + power-/magnitude-only systems with no phasor measurement). + +Both share one measurement assembly ([`build_se_atoms`](@ref)) that consumes the +*same* `data["meas"]` dictionary as the JuMP estimators, and both converge to the +identical WLS optimum. PGM normalises measurements per-unit, so the benchmark +data is built per-unit (`V_base = u_rated`, `S_base = base_power_3p = 1e6`). + +`src/bare/pgm_se.jl` holds the implementation; `test/literal/pgm/` holds the +PGM-derived golden data (regenerate with `test/literal/pgm/generate_golden.py`). + + ## Files ``` @@ -93,6 +123,8 @@ test/literal/ | 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 | +| 6 | PGM solve options (`:iterative_linear`/`:newton_raphson`) reproduce PowerGridModel node voltages + line flows on its own SE examples | max\|U\| `<1e-3 V` (≈1e-7 typical), both methods | +| 7 | the PGM solve options run unchanged on the four-wire EN feeder (recover the IVREN PF, agree with `solve_wls`) | max\|U\| `<1e-4` | ### Benchmark (3-bus 4-wire, σ = 0.01, 4 noise seeds) diff --git a/docs/src/literal_pmdse.md b/docs/src/literal_pmdse.md index fbda7b1..899dc03 100644 --- a/docs/src/literal_pmdse.md +++ b/docs/src/literal_pmdse.md @@ -48,7 +48,7 @@ res = _PMDSE.solve_mc_se_literal(math; estimator=:wls, reference=:sota) ## `solve_mc_se_literal` ```julia -solve_mc_se_literal(data; estimator = :wls, reference = :sota, +solve_mc_se_literal(data; estimator = :wls, method = nothing, reference = :sota, ref_values = nothing, maxiter = 50, tol = 1e-9, verbose = false) ``` @@ -57,6 +57,7 @@ solve_mc_se_literal(data; estimator = :wls, reference = :sota, |----------|--------|---------| | `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). | +| `method` | `nothing` (default), `:iterative_linear`, `:newton_raphson` | for `estimator = :wls`, selects a PowerGridModel solve option (see below); `nothing` keeps the in-house Gauss–Newton WLS. | | `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. | @@ -75,6 +76,25 @@ optional `"reference"` entry overrides the `reference` keyword. | `: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`. | +### PowerGridModel solve options (`method`) + +With `estimator = :wls` the two +[PowerGridModel](https://github.com/PowerGridModel/power-grid-model) WLS *solve +options* can be selected through `method`. Both are general (four-wire, +explicit-neutral capable — the single-phase PGM network is the one-conductor, +no-neutral special case) and are validated against PowerGridModel's own +state-estimation examples (test gates 6–7). + +| `method` | algorithm | notes | +|----------|-----------|-------| +| `:iterative_linear` | measurements linearised to complex currents / voltage-phasors at the previous voltages; a constant linear WLS system re-solved each iteration | PowerGridModel's default; the slack-angle gauge is pinned automatically when no voltage angle is measured. | +| `:newton_raphson` | Gauss–Newton on the nonlinear WLS, warm-started from one `iterative_linear` solve | identical optimum to `:iterative_linear`; warm-start avoids the flat-start degeneracy of power-/magnitude-only systems. | + +```julia +res = solve_mc_se_literal(data; estimator = :wls, method = :iterative_linear) +res = solve_mc_se_literal(data; estimator = :wls, method = :newton_raphson) +``` + ```julia res_wls = _PMDSE.solve_mc_se_literal(math; estimator=:wls) res_wlav = _PMDSE.solve_mc_se_literal(math; estimator=:wlav) # rejects gross errors diff --git a/src/PowerModelsDistributionStateEstimation.jl b/src/PowerModelsDistributionStateEstimation.jl index c992259..85ac927 100644 --- a/src/PowerModelsDistributionStateEstimation.jl +++ b/src/PowerModelsDistributionStateEstimation.jl @@ -88,6 +88,7 @@ include("prob/se_en.jl") include("bare/literal_core.jl") include("bare/measurement_model.jl") include("bare/solve_wls.jl") +include("bare/pgm_se.jl") include("bare/solve_wlav.jl") include("bare/mle.jl") include("bare/literal_pmdse.jl") diff --git a/src/bare/literal_pmdse.jl b/src/bare/literal_pmdse.jl index c56c567..1482004 100644 --- a/src/bare/literal_pmdse.jl +++ b/src/bare/literal_pmdse.jl @@ -9,8 +9,8 @@ ################################################################################ """ - solve_mc_se_literal(data; estimator=:wls, reference=:sota, ref_values=nothing, - maxiter=50, tol=1e-9, verbose=false, kwargs...) + solve_mc_se_literal(data; estimator=:wls, method=nothing, reference=:sota, + ref_values=nothing, maxiter=50, tol=1e-9, verbose=false, kwargs...) Run the literal (matrix-based, JuMP-free) explicit-neutral state estimator. @@ -19,9 +19,19 @@ Run the literal (matrix-based, JuMP-free) explicit-neutral state estimator. (e.g. via `add_measurements!`) and an optional `data["se_settings"]` (`"rescaler"`, and optionally `"reference"`). - `estimator` : `:wls` (default), `:wlav`, or `:mle`. +- `method` : selects the WLS *solve option* when `estimator == :wls`: + * `nothing` (default) — the in-house Gauss–Newton WLS (`solve_wls`). + * `:iterative_linear` — the PowerGridModel default: measurements are linearised + to complex currents/voltage-phasors each iteration and a constant linear WLS + system is re-solved ([`solve_se_il`](@ref)). + * `:newton_raphson` — PowerGridModel's nonlinear Gauss–Newton, warm-started + from one `iterative_linear` solve ([`solve_se_nr`](@ref)). + Both PGM methods converge to the same WLS optimum; they are validated against + PowerGridModel's own state-estimation examples (see `test/literal`). - `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). + `:full_slack` (fix the whole reference-bus phasor), or `:prop` (fix only the + grounded/neutral terminals — the slack-angle gauge is then taken from a voltage + phasor measurement, or pinned automatically when only magnitudes are measured). - `ref_values` : optional `(vr_dict, vi_dict)` over reference-bus terminals for `:full_slack`. @@ -29,6 +39,7 @@ 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, + method::Union{Nothing,Symbol} = nothing, reference::Symbol = :sota, ref_values = nothing, maxiter::Int = 50, tol::Float64 = 1e-9, verbose::Bool = false, kwargs...) @@ -36,10 +47,19 @@ function solve_mc_se_literal(data::Dict{String,<:Any}; estimator::Symbol = :wls, se = get(data, "se_settings", Dict{String,Any}()) rescaler = Float64(get(se, "rescaler", 1.0)) reference = Symbol(get(se, "reference", reference)) + method = method === nothing ? (haskey(se, "method") ? Symbol(se["method"]) : nothing) : method lm = LiteralModel(data; reference = reference, ref_values = ref_values) - model = build_se_model(lm, data; rescaler = rescaler) + # PowerGridModel solve options (atom-based, four-wire general) + if estimator == :wls && method in (:iterative_linear, :newton_raphson) + atoms = build_se_atoms(lm, data; rescaler = rescaler) + return method == :iterative_linear ? + solve_se_il(lm, atoms; maxiter = maxiter, tol = tol, verbose = verbose) : + solve_se_nr(lm, atoms; maxiter = maxiter, tol = tol, verbose = verbose) + end + + model = build_se_model(lm, data; rescaler = rescaler) if estimator == :wls return solve_wls(model; maxiter = maxiter, tol = tol, verbose = verbose) elseif estimator == :wlav diff --git a/src/bare/pgm_se.jl b/src/bare/pgm_se.jl new file mode 100644 index 0000000..b1a62b7 --- /dev/null +++ b/src/bare/pgm_se.jl @@ -0,0 +1,603 @@ +################################################################################ +# Literal PMDSE # +# pgm_se.jl : the two PowerGridModel (PGM) state-estimation *solve options* # +# re-implemented on the literal nodal model and made GENERAL # +# (four-wire, explicit-neutral capable). # +# # +# PGM (https://github.com/PowerGridModel/power-grid-model) solves the # +# weighted-least-squares DSSE with two interchangeable algorithms: # +# # +# • iterative_linear (PGM default) — every measurement is *linearised* to an # +# equivalent complex-current / voltage-phasor measurement using the voltage # +# angle from the previous iteration; a (constant) linear WLS system is then # +# re-solved until the voltages stop moving. # +# • newton_raphson — the nonlinear WLS is solved directly by Gauss–Newton on # +# the real measurement model `z = h(x)+e`. # +# # +# Both minimise ½ (h(x)−z)ᵀ W (h(x)−z), W = diag(1/σ²), # +# i.e. they converge to the *same* state; PGM ships them as two methods and # +# this file mirrors that. Unlike PGM (which is single-phase / Kron-reduced) # +# the implementation here keeps the full per-(bus,terminal) rectangular state # +# `x=[vr;vi]`, so it is general: an explicit neutral conductor is supported # +# (phase-to-neutral voltages/powers, neutral injection `= −Σ phase currents`), # +# and the single-phase PGM networks are simply the special case of one # +# conductor per bus with no neutral. # +# # +# The measurement model is shared between the two methods through a list of # +# `SEAtom`s, each carrying the (real/imag parts of the) complex coefficient # +# vectors of the two linear functionals it needs: # +# ΔU(x) = cu·U (phase-to-neutral voltage) # +# I(x) = ci·U (a branch or nodal-injection current, via Y) # +# From these every measurement type (voltage phasor/magnitude, branch power, # +# branch current, injection power/current, zero-injection) is expressed. # +################################################################################ + +"weight given to a zero-injection KCL pseudo-measurement (1/σ²)" +const _ZINJ_W = 1.0e10 + +""" + SEAtom + +One scalar measurement (or one real/imaginary component of one complex +measurement) of the WLS problem, stored as the real and imaginary parts of the +complex coefficient vectors of its two linear functionals over the nodal complex +voltage `U` (length `n`): + +* `cu = cur + im*cui` : the phase-to-neutral voltage `ΔU = cu·U` + (used by voltage measurements and as the `U` in `S = ΔU·conj(I)`). +* `ci = cir + im*cii` : a current `I = ci·U` (a branch current row, or a + `Y_bus` row for a nodal injection). + +`kind`: +* `:vre` / `:vim` — measured real / imaginary part of `ΔU` (a voltage phasor + component); one residual row, `z = z1`. +* `:vmag` — measured magnitude `|ΔU|`; one (nonlinear) residual row for + Newton–Raphson, two (linearised phasor) rows for iterative-linear. +* `:power` — measured `(P,Q) = z1,z2` of `S = sign·ΔU·conj(I)`; two rows. +* `:cinj` — measured complex current `I = z1+im*z2`; two rows. +* `:zinj` — zero-injection constraint `I = 0`; two rows. +""" +Base.@kwdef struct SEAtom + kind::Symbol + cur::Vector{Float64} = Float64[] + cui::Vector{Float64} = Float64[] + cir::Vector{Float64} = Float64[] + cii::Vector{Float64} = Float64[] + z1::Float64 = 0.0 + z2::Float64 = 0.0 + sign::Float64 = 1.0 + sigma::Float64 = 1.0 # σ of a voltage / current measurement + sigp::Float64 = 1.0 # σ of the active-power component + sigq::Float64 = 1.0 # σ of the reactive-power component + meta::NamedTuple = NamedTuple() +end + +# ---- complex-coefficient helpers ------------------------------------------ + +"`(re,im)` of `Σ_k (cr[k]+im*ci[k])·(vr[k]+im*vi[k])` (ForwardDiff-friendly)" +@inline function _cval(cr::Vector{Float64}, ci::Vector{Float64}, vr, vi) + re = zero(eltype(vr)); im_ = zero(eltype(vr)) + @inbounds for k in eachindex(cr) + re += cr[k] * vr[k] - ci[k] * vi[k] + im_ += ci[k] * vr[k] + cr[k] * vi[k] + end + return re, im_ +end + +"same as `_cval` but for a fixed `ComplexF64` voltage vector `U`" +@inline function _cvalU(cr::Vector{Float64}, ci::Vector{Float64}, U::Vector{ComplexF64}) + s = zero(ComplexF64) + @inbounds for k in eachindex(cr) + s += (cr[k] + im * ci[k]) * U[k] + end + return s +end + +"unit complex coefficient `e_k` (length n)" +_unit_coeff(k::Int, n::Int) = (v = zeros(Float64, n); v[k] = 1.0; v) + +# ---- branch admittance blocks (from & to side) ---------------------------- + +""" + _pgm_branch_blocks(lm) -> Dict + +Per-branch full 2×2 (block) admittance: `Yff,Yft,Ytf,Ytt` (complex matrices over +the branch conductors), the from/to node positions `F,T`, the conductor list, +and the from/to neutral node positions (`Fn`/`Tn`, `nothing` if no neutral). +""" +function _pgm_branch_blocks(lm::LiteralModel) + bd = Dict{Int,Any}() + for (_, 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_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 = [lm.node_index[(br["f_bus"], c)] for c in fc] + T = [lm.node_index[(br["t_bus"], c)] for c in tc] + Fn = haskey(lm.node_index, (br["f_bus"], _N_IDX)) ? lm.node_index[(br["f_bus"], _N_IDX)] : nothing + Tn = haskey(lm.node_index, (br["t_bus"], _N_IDX)) ? lm.node_index[(br["t_bus"], _N_IDX)] : nothing + bd[br["index"]] = (Yff = Ybr .+ Ysh_fr, Yft = .-Ybr, Ytf = .-Ybr, Ytt = Ybr .+ Ysh_to, + F = F, T = T, fc = fc, tc = tc, + f_bus = br["f_bus"], t_bus = br["t_bus"], Fn = Fn, Tn = Tn) + end + return bd +end + +"neutral node position of `bus` (or `nothing` when the bus has no neutral)" +function _neutral_pos(lm::LiteralModel, bus::Int) + haskey(lm.node_index, (bus, _N_IDX)) ? lm.node_index[(bus, _N_IDX)] : nothing +end + +"phase-to-neutral coefficient `e_kc − e_kn` (or `e_kc` if no neutral) of length n" +function _delta_coeff(lm::LiteralModel, bus::Int, c::Int) + n = lm.n + cu = _unit_coeff(lm.node_index[(bus, c)], n) + kn = _neutral_pos(lm, bus) + kn !== nothing && c != _N_IDX && (cu[kn] -= 1.0) + return cu +end + +""" + _kalman_combine(a, b) -> per-phase combined `Normal` distributions + +Inverse-variance ("Kalman") combine of two per-conductor measurement vectors of +the *same* physical quantity (PGM §3.3): `1/σ² = Σ 1/σ_k²`, +`μ = σ² Σ μ_k/σ_k²`. A disabled sensor (`σ → ∞`) contributes nothing. +""" +function _kalman_combine(a, b) + return [begin + μa, σa = _zsigma(a[i], 1.0e-8); μb, σb = _zsigma(b[i], 1.0e-8) + ia = 1.0 / σa^2; ib = 1.0 / σb^2; iv = ia + ib + iv > 0 ? _DST.Normal((μa * ia + μb * ib) / iv, sqrt(1.0 / iv)) : _DST.Normal(μa, Inf) + end for i in eachindex(a)] +end + +"store `dst` under `dict[key][part]`, Kalman-combining a repeated measurement" +function _accum_meas!(dict, key, part, dst) + g = get!(dict, key, Dict{Symbol,Any}()) + g[part] = haskey(g, part) ? _kalman_combine(g[part], dst) : dst + return nothing +end + +"node position of the reference bus' first non-neutral terminal" +function _ref_phase_pos(lm::LiteralModel) + terms = lm.math["bus"][string(lm.ref_bus)]["terminals"] + ph = first(t for t in terms if t != _N_IDX) + return lm.node_index[(lm.ref_bus, ph)] +end + +# ---- measurement → atoms --------------------------------------------------- + +""" + build_se_atoms(lm, math; rescaler=1.0, σ_exact=1e-8, zinj_sigma=1e-5) + +Translate `math["meas"]` (the PMD-style measurement dictionary, same one the +JuMP estimators consume) into the flat list of [`SEAtom`](@ref)s used by both +PGM solve options. Appliance (`load`/`gen`) injection measurements are +aggregated onto the bus, exactly as PGM aggregates appliance powers into the bus +injection; current injections add the neutral return `= −Σ phase currents`. +Genuine zero-injection terminals get a `z=0` KCL pseudo-measurement. + +Recognised `var`s (per non-neutral conductor unless noted): +`:vr`,`:vi` (phasor components) · `:vm`,`:vmn` (magnitude) · +`:p`,`:q` (branch power, `meta["side"]` ∈ `(:from,:to)`, default `:from`) · +`:cr`,`:ci` (branch current) · `:pd`,`:qd`,`:pg`,`:qg` (injection power) · +`:crd`,`:cid`,`:crg`,`:cig` (injection current). +""" +function build_se_atoms(lm::LiteralModel, math::Dict; rescaler::Float64 = 1.0, + σ_exact::Float64 = 1.0e-8, zinj_sigma::Float64 = 1.0e-5) + n = lm.n + Ybus = complex.(lm.G, lm.B) + bd = _pgm_branch_blocks(lm) + atoms = SEAtom[] + rs(σ) = rescaler * σ + + # accumulators for appliance injection aggregation, keyed by node position + pw_S = Dict{Int,ComplexF64}(); pw_vp = Dict{Int,Float64}(); pw_vq = Dict{Int,Float64}() + cu_I = Dict{Int,ComplexF64}(); cu_v = Dict{Int,Float64}() + touched = Set{Int}() + sourced = Set{Int}() + for (_, ld) in get(math, "load", Dict{String,Any}()), c in ld["connections"] + push!(sourced, lm.node_index[(ld["load_bus"], c)]) + end + for (_, g) in get(math, "gen", Dict{String,Any}()), c in g["connections"] + push!(sourced, lm.node_index[(g["gen_bus"], c)]) + end + + # paired power / current branch measurements, keyed by (cmp_id, side, part) + bpow = Dict{Tuple{Int,Symbol},Dict{Symbol,Any}}() + bcur = Dict{Tuple{Int,Symbol},Dict{Symbol,Any}}() + # paired injection measurements, keyed by (cmp, cmp_id) + ipow = Dict{Tuple{Symbol,Int},Dict{Symbol,Any}}() + icur = Dict{Tuple{Symbol,Int},Dict{Symbol,Any}}() + + for (mid, meas) in math["meas"] + var = meas["var"]; cmp = meas["cmp"]; cmp_id = meas["cmp_id"]; dst = meas["dst"] + side = Symbol(get(meas, "side", :from)) + if var in (:vr, :vi) # voltage phasor component + bus = cmp_id + for (idx, c) in enumerate(_nonneutral(_active_connections(math, cmp, cmp_id))) + μ, σ = _zsigma(dst[idx], σ_exact) + cu = _delta_coeff(lm, bus, c) + push!(atoms, SEAtom(kind = (var == :vr ? :vre : :vim), + cur = real(cu), cui = imag(cu), + z1 = μ, sigma = rs(σ), meta = (mid = mid, var = var, c = c))) + end + elseif var in (:vm, :vmn) # voltage magnitude |ΔU| + bus = cmp_id + for (idx, c) in enumerate(_nonneutral(_active_connections(math, cmp, cmp_id))) + μ, σ = _zsigma(dst[idx], σ_exact) + cu = _delta_coeff(lm, bus, c) + push!(atoms, SEAtom(kind = :vmag, cur = real(cu), cui = imag(cu), + z1 = μ, sigma = rs(σ), meta = (mid = mid, var = var, c = c))) + end + elseif var == :p + _accum_meas!(bpow, (cmp_id, side), :p, dst) + elseif var == :q + _accum_meas!(bpow, (cmp_id, side), :q, dst) + elseif var == :cr + _accum_meas!(bcur, (cmp_id, side), :r, dst) + elseif var == :ci + _accum_meas!(bcur, (cmp_id, side), :i, dst) + elseif var in (:pd, :pg) + _accum_meas!(ipow, (cmp, cmp_id), :p, dst) + elseif var in (:qd, :qg) + _accum_meas!(ipow, (cmp, cmp_id), :q, dst) + elseif var in (:crd, :crg) + _accum_meas!(icur, (cmp, cmp_id), :r, dst) + elseif var in (:cid, :cig) + _accum_meas!(icur, (cmp, cmp_id), :i, dst) + else + error("Literal PMDSE (PGM): measurement var :$(var) not supported by the PGM solve options") + end + end + + # ---- branch power measurements ---------------------------------------- + for ((cmp_id, side), pq) in bpow + haskey(pq, :p) && haskey(pq, :q) || error("Literal PMDSE (PGM): branch $(cmp_id) $(side) power needs both :p and :q") + b = bd[cmp_id] + Yself = side == :from ? b.Yff : b.Ytt + Yother = side == :from ? b.Yft : b.Ytf + Fpos = side == :from ? b.F : b.T + Opos = side == :from ? b.T : b.F + sbus = side == :from ? b.f_bus : b.t_bus + for (idx, c) in enumerate(_nonneutral(b.fc)) + a = findfirst(isequal(c), b.fc) + ci = zeros(ComplexF64, n) + @inbounds for q in eachindex(Fpos); ci[Fpos[q]] += Yself[a, q]; end + @inbounds for q in eachindex(Opos); ci[Opos[q]] += Yother[a, q]; end + cu = _delta_coeff(lm, sbus, c) + μp, σp = _zsigma(pq[:p][idx], σ_exact); μq, σq = _zsigma(pq[:q][idx], σ_exact) + push!(atoms, SEAtom(kind = :power, cur = real(cu), cui = imag(cu), + cir = real(ci), cii = imag(ci), z1 = μp, z2 = μq, sign = 1.0, + sigp = rs(σp), sigq = rs(σq), meta = (var = :p, branch = cmp_id, side = side, c = c))) + end + end + + # ---- branch current measurements -------------------------------------- + for ((cmp_id, side), ri) in bcur + haskey(ri, :r) && haskey(ri, :i) || error("Literal PMDSE (PGM): branch $(cmp_id) $(side) current needs both :cr and :ci") + b = bd[cmp_id] + Yself = side == :from ? b.Yff : b.Ytt + Yother = side == :from ? b.Yft : b.Ytf + Fpos = side == :from ? b.F : b.T + Opos = side == :from ? b.T : b.F + for (idx, c) in enumerate(_nonneutral(b.fc)) + a = findfirst(isequal(c), b.fc) + ci = zeros(ComplexF64, n) + @inbounds for q in eachindex(Fpos); ci[Fpos[q]] += Yself[a, q]; end + @inbounds for q in eachindex(Opos); ci[Opos[q]] += Yother[a, q]; end + μr, σr = _zsigma(ri[:r][idx], σ_exact); μi, _ = _zsigma(ri[:i][idx], σ_exact) + push!(atoms, SEAtom(kind = :cinj, cir = real(ci), cii = imag(ci), + z1 = μr, z2 = μi, sigma = rs(σr), meta = (var = :cr, branch = cmp_id, side = side, c = c))) + end + end + + # ---- injection power measurements (aggregate per node-conductor) ------- + for ((cmp, cmp_id), pq) in ipow + haskey(pq, :p) && haskey(pq, :q) || error("Literal PMDSE (PGM): injection $(cmp) $(cmp_id) power needs both p and q") + isgen = cmp == :gen + bus = isgen ? math["gen"][string(cmp_id)]["gen_bus"] : math["load"][string(cmp_id)]["load_bus"] + conns = isgen ? math["gen"][string(cmp_id)]["connections"] : math["load"][string(cmp_id)]["connections"] + s = isgen ? 1.0 : -1.0 + for (idx, c) in enumerate(_nonneutral(conns)) + k = lm.node_index[(bus, c)] + μp, σp = _zsigma(pq[:p][idx], σ_exact); μq, σq = _zsigma(pq[:q][idx], σ_exact) + pw_S[k] = get(pw_S, k, 0.0 + 0im) + s * (μp + im * μq) + pw_vp[k] = get(pw_vp, k, 0.0) + (rs(σp))^2 + pw_vq[k] = get(pw_vq, k, 0.0) + (rs(σq))^2 + push!(touched, k) + end + end + for (k, S) in pw_S + cu = _unit_coeff(k, n) + # phase-to-neutral: subtract the neutral of the bus owning node k + (bus, c) = lm.nodes[k]; kn = _neutral_pos(lm, bus) + kn !== nothing && c != _N_IDX && (cu[kn] -= 1.0) + ci = collect(@view Ybus[k, :]) + push!(atoms, SEAtom(kind = :power, cur = real(cu), cui = imag(cu), + cir = real(ci), cii = imag(ci), z1 = real(S), z2 = imag(S), sign = 1.0, + sigp = sqrt(pw_vp[k]), sigq = sqrt(pw_vq[k]), meta = (var = :pinj, node = k))) + end + + # ---- injection current measurements (aggregate + neutral return) ------ + for ((cmp, cmp_id), ri) in icur + haskey(ri, :r) && haskey(ri, :i) || error("Literal PMDSE (PGM): injection $(cmp) $(cmp_id) current needs both cr and ci") + isgen = cmp == :gen + bus = isgen ? math["gen"][string(cmp_id)]["gen_bus"] : math["load"][string(cmp_id)]["load_bus"] + conns = isgen ? math["gen"][string(cmp_id)]["connections"] : math["load"][string(cmp_id)]["connections"] + s = isgen ? 1.0 : -1.0 + Isum = 0.0 + 0im; vsum = 0.0 + for (idx, c) in enumerate(_nonneutral(conns)) + k = lm.node_index[(bus, c)] + μr, σr = _zsigma(ri[:r][idx], σ_exact); μi, _ = _zsigma(ri[:i][idx], σ_exact) + Iph = s * (μr + im * μi) + cu_I[k] = get(cu_I, k, 0.0 + 0im) + Iph + cu_v[k] = get(cu_v, k, 0.0) + (rs(σr))^2 + push!(touched, k) + Isum += Iph; vsum += (rs(σr))^2 + end + if _N_IDX in conns # neutral return = −Σ phase currents + kn = lm.node_index[(bus, _N_IDX)] + cu_I[kn] = get(cu_I, kn, 0.0 + 0im) - Isum + cu_v[kn] = get(cu_v, kn, 0.0) + vsum + push!(touched, kn) + end + end + for (k, I) in cu_I + ci = collect(@view Ybus[k, :]) + push!(atoms, SEAtom(kind = :cinj, cir = real(ci), cii = imag(ci), + z1 = real(I), z2 = imag(I), sigma = sqrt(cu_v[k]), meta = (var = :cinj, node = k))) + end + + # ---- zero-injection KCL pseudo-measurements --------------------------- + for k in lm.inj_nodes + (k in touched || k in sourced) && continue + ci = collect(@view Ybus[k, :]) + push!(atoms, SEAtom(kind = :zinj, cir = real(ci), cii = imag(ci), + sigma = zinj_sigma, meta = (var = :zinj, node = k))) + end + + return atoms +end + +# ---- shared reference / gauge handling ------------------------------------ + +"the free indices and fixed state for a solve, adding the slack-angle gauge +`Im(U_ref)=0` when the measurements carry no absolute angle information" +function _gauge(lm::LiteralModel, atoms::Vector{SEAtom}) + has_angle = any(a -> a.kind == :vim, atoms) + fixed = copy(lm.fixed_mask) + if !has_angle + fixed[lm.n + _ref_phase_pos(lm)] = true # pin the global rotation + end + free = [i for i in 1:2lm.n if !fixed[i]] + return free, copy(lm.x_fixed), has_angle +end + +"reconstruct the full `2n` state (ForwardDiff-friendly) from the free entries" +function _expand_free(xfix::Vector{Float64}, free::Vector{Int}, xf::AbstractVector{T}) where {T} + x = convert(Vector{T}, xfix) + @inbounds for (j, i) in enumerate(free); x[i] = xf[j]; end + return x +end + +"complex flat-start voltages for `lm` (balanced phasors, neutral ≈ 0)" +function _flat_U(lm::LiteralModel) + x = flat_start(lm.math, lm.nodes, lm.node_index, lm.n) + return complex.(x[1:lm.n], x[lm.n+1:2lm.n]) +end + +# ---- iterative-linear method ---------------------------------------------- + +""" + solve_se_il(lm, atoms; maxiter=100, tol=1e-9, verbose=false) + +PowerGridModel `iterative_linear` state estimation on the literal model. Each +iteration re-linearises the power/current and voltage-magnitude measurements at +the latest voltages (only the right-hand side changes; the measurement matrix is +constant) and solves the linear WLS `min ‖√W (A x − b)‖₂` for the new +rectangular state `x=[vr;vi]`, until `max_i |Uᵢ−Uᵢ_prev| < tol`. +""" +function solve_se_il(lm::LiteralModel, atoms::Vector{SEAtom}; maxiter::Int = 100, + tol::Float64 = 1.0e-9, verbose::Bool = false) + t0 = time() + n = lm.n + free, xfix, _ = _gauge(lm, atoms) + U = _flat_U(lm) + # seed fixed nodes (e.g. a `:full_slack` reference bus) to their known phasor + @inbounds for i in 1:n + lm.fixed_mask[i] && (U[i] = xfix[i] + im * xfix[n+i]) + end + term = :maxiter; iters = 0 + + # the measurement matrix A (constant) — rhs/weights are rebuilt per iteration + rows = Vector{Float64}[] + for a in atoms + if a.kind == :vre + push!(rows, vcat(a.cur, -a.cui)) + elseif a.kind == :vim + push!(rows, vcat(a.cui, a.cur)) + elseif a.kind == :vmag + push!(rows, vcat(a.cur, -a.cui)); push!(rows, vcat(a.cui, a.cur)) + elseif a.kind in (:power, :cinj, :zinj) + push!(rows, vcat(a.cir, -a.cii)); push!(rows, vcat(a.cii, a.cir)) + end + end + A = permutedims(reduce(hcat, rows)) # m × 2n + Afree = A[:, free] + Afix_b = A * xfix # fixed-column contribution + w = _il_weights(atoms) + + for it in 1:maxiter + iters = it + b = _il_rhs(atoms, U) + beff = b .- Afix_b + sw = sqrt.(w) + xf = (sw .* Afree) \ (sw .* beff) + xfull = copy(xfix); @inbounds for (j, i) in enumerate(free); xfull[i] = xf[j]; end + Unew = complex.(xfull[1:n], xfull[n+1:2n]) + dev = maximum(abs.(Unew .- U)) + U = Unew + verbose && println(" il it $it max|ΔU|=$(dev)") + if dev < tol; term = :converged; break; end + end + + x_full = vcat(real.(U), imag.(U)) + obj, gcond, grank = _atoms_diagnostics(atoms, lm, x_full, free, xfix) + sol = state_solution(lm, x_full) + return LiteralResult(term, iters, obj, time() - t0, gcond, grank, + x_full[free], x_full, sol, :iterative_linear) +end + +"per-row weights of the iterative-linear measurement matrix (constant)" +function _il_weights(atoms::Vector{SEAtom}) + w = Float64[] + for a in atoms + if a.kind in (:vre, :vim) + push!(w, 1.0 / a.sigma^2) + elseif a.kind == :vmag + push!(w, 1.0 / a.sigma^2); push!(w, 1.0 / a.sigma^2) + elseif a.kind == :power + wS = 1.0 / (a.sigp^2 + a.sigq^2) # PGM: σ_S² = σ_P² + σ_Q² + push!(w, wS); push!(w, wS) + elseif a.kind == :cinj + push!(w, 1.0 / a.sigma^2); push!(w, 1.0 / a.sigma^2) + elseif a.kind == :zinj + push!(w, _ZINJ_W); push!(w, _ZINJ_W) + end + end + return w +end + +"right-hand side of the iterative-linear system, re-linearised at voltages `U`" +function _il_rhs(atoms::Vector{SEAtom}, U::Vector{ComplexF64}) + b = Float64[] + for a in atoms + if a.kind in (:vre, :vim) + push!(b, a.z1) + elseif a.kind == :vmag + du = _cvalU(a.cur, a.cui, U) + tgt = a.z1 * (abs(du) > 0 ? du / abs(du) : 1.0 + 0im) # use previous angle + push!(b, real(tgt)); push!(b, imag(tgt)) + elseif a.kind == :power + du = _cvalU(a.cur, a.cui, U) + Imeas = conj(a.sign * (a.z1 + im * a.z2) / du) # S → equivalent current + push!(b, real(Imeas)); push!(b, imag(Imeas)) + elseif a.kind == :cinj + push!(b, a.z1); push!(b, a.z2) + elseif a.kind == :zinj + push!(b, 0.0); push!(b, 0.0) + end + end + return b +end + +# ---- newton-raphson method ------------------------------------------------ + +"the nonlinear measurement prediction `h(x)` over the free state (ForwardDiff ok)" +function _nr_predict(atoms::Vector{SEAtom}, n::Int, xfix::Vector{Float64}, + free::Vector{Int}, xf::AbstractVector{T}) where {T} + x = _expand_free(xfix, free, xf) + vr = @view x[1:n]; vi = @view x[n+1:2n] + h = T[] + for a in atoms + if a.kind == :vre + re, _ = _cval(a.cur, a.cui, vr, vi); push!(h, re) + elseif a.kind == :vim + _, im_ = _cval(a.cur, a.cui, vr, vi); push!(h, im_) + elseif a.kind == :vmag + re, im_ = _cval(a.cur, a.cui, vr, vi); push!(h, sqrt(re^2 + im_^2)) + elseif a.kind == :power + dur, dui = _cval(a.cur, a.cui, vr, vi); ir, ii = _cval(a.cir, a.cii, vr, vi) + push!(h, a.sign * (dur * ir + dui * ii)); push!(h, a.sign * (dui * ir - dur * ii)) + elseif a.kind in (:cinj, :zinj) + ir, ii = _cval(a.cir, a.cii, vr, vi); push!(h, ir); push!(h, ii) + end + end + return h +end + +"measured values `z` and weights `w` aligned with `_nr_predict` rows" +function _nr_zw(atoms::Vector{SEAtom}) + z = Float64[]; w = Float64[] + for a in atoms + if a.kind in (:vre, :vim) + push!(z, a.z1); push!(w, 1.0 / a.sigma^2) + elseif a.kind == :vmag + push!(z, a.z1); push!(w, 1.0 / a.sigma^2) + elseif a.kind == :power + push!(z, a.z1); push!(w, 1.0 / a.sigp^2) + push!(z, a.z2); push!(w, 1.0 / a.sigq^2) + elseif a.kind == :cinj + push!(z, a.z1); push!(w, 1.0 / a.sigma^2) + push!(z, a.z2); push!(w, 1.0 / a.sigma^2) + elseif a.kind == :zinj + push!(z, 0.0); push!(w, _ZINJ_W); push!(z, 0.0); push!(w, _ZINJ_W) + end + end + return z, w +end + +""" + solve_se_nr(lm, atoms; maxiter=50, tol=1e-9, warm=true, verbose=false) + +PowerGridModel `newton_raphson` state estimation: Gauss–Newton on the nonlinear +WLS `min ½‖√W (z − h(x))‖₂²`, `H = ∂h/∂x` via `ForwardDiff`, with the same +QR/orthogonal step as `solve_wls`. By default the iteration is *warm-started* +from one `iterative_linear` solve (`warm=true`); the flat start is degenerate for +power-/magnitude-only systems with no phasor measurement (the angle Jacobian +vanishes), so warm-starting makes Newton–Raphson as robust as PGM's augmented +formulation while converging to the identical WLS optimum. +""" +function solve_se_nr(lm::LiteralModel, atoms::Vector{SEAtom}; maxiter::Int = 50, + tol::Float64 = 1.0e-9, warm::Bool = true, verbose::Bool = false) + t0 = time() + n = lm.n + free, xfix, _ = _gauge(lm, atoms) + z, w = _nr_zw(atoms) + sqrtW = sqrt.(w) + + if warm + x_full0 = solve_se_il(lm, atoms; maxiter = maxiter, tol = tol).x_full + xf = x_full0[free] + else + xf = _flat_U(lm) |> U -> vcat(real.(U), imag.(U))[free] + end + + term = :maxiter; iters = 0 + for it in 1:maxiter + iters = it + h = _nr_predict(atoms, n, xfix, free, xf) + r = z .- h + H = ForwardDiff.jacobian(xx -> _nr_predict(atoms, n, xfix, free, xx), xf) + G = Symmetric_full(transpose(H) * (w .* H)) + rhs = transpose(H) * (w .* r) + Δ, _ = _solve_gain(G, rhs, sqrtW .* H, sqrtW .* r) + xf = xf .+ Δ + verbose && println(" nr it $it ‖Δ‖∞=$(LinearAlgebra.norm(Δ, Inf))") + if LinearAlgebra.norm(Δ, Inf) < tol; term = :converged; break; end + end + + x_full = _expand_free(xfix, free, xf) + h = _nr_predict(atoms, n, xfix, free, xf); r = z .- h + obj = sum(w .* r .^ 2) + H = ForwardDiff.jacobian(xx -> _nr_predict(atoms, n, xfix, free, xx), xf) + gcond, grank = _gain_diag(transpose(H) * (w .* H)) + sol = state_solution(lm, x_full) + return LiteralResult(term, iters, obj, time() - t0, gcond, grank, xf, x_full, sol, :newton_raphson) +end + +"weighted SSR + gain diagnostics of an atom set at a full state" +function _atoms_diagnostics(atoms::Vector{SEAtom}, lm::LiteralModel, x_full::Vector{Float64}, + free::Vector{Int}, xfix::Vector{Float64}) + z, w = _nr_zw(atoms) + xf = x_full[free] + h = _nr_predict(atoms, lm.n, xfix, free, xf) + obj = sum(w .* (z .- h) .^ 2) + H = ForwardDiff.jacobian(xx -> _nr_predict(atoms, lm.n, xfix, free, xx), xf) + gcond, grank = _gain_diag(transpose(H) * (w .* H)) + return obj, gcond, grank +end diff --git a/src/core/export.jl b/src/core/export.jl index 4af0481..62820cf 100644 --- a/src/core/export.jl +++ b/src/core/export.jl @@ -15,6 +15,8 @@ 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 +# PowerGridModel solve options (four-wire general) + measurement assembly +export build_se_atoms, solve_se_il, solve_se_nr, SEAtom # so that users do not need to import JuMP to use a solver with PowerModelsDistribution import JuMP: optimizer_with_attributes diff --git a/test/literal/pgm/README.md b/test/literal/pgm/README.md new file mode 100644 index 0000000..d148b68 --- /dev/null +++ b/test/literal/pgm/README.md @@ -0,0 +1,34 @@ +# PowerGridModel benchmark for the Literal PMDSE solve options + +This folder pins the Literal PMDSE PowerGridModel (PGM) solve options +(`:iterative_linear`, `:newton_raphson`) to PGM's own state-estimation results. + +* `generate_golden.py` — runs [power-grid-model](https://github.com/PowerGridModel/power-grid-model) + on a handful of its symmetric state-estimation validation cases + (`tests/data/state_estimation/*`) and writes the golden data. +* `golden.json` — machine-readable network + sensors + the node voltages PGM + computes for each calculation method (informational / diffable). +* `pgm_cases_data.jl` — the same data as a `const PGM_CASES = [...]` literal; its + contents are embedded verbatim in `../pgm_cases.jl` (the test does **not** read + this file at runtime, so the suite needs neither Python nor PGM installed). + +## Cases + +| name | PGM folder | nodes | methods | what it exercises | +|---|---|---|---|---| +| `single_node` | `single-node-source-sym-voltage-sensor` | 1 | il, nr | trivial voltage-phasor pin | +| `os1_2msr` | `1os2msr` | 3 | il, nr | redundant set, branch from/to power, source/load injection, a disabled (σ=1e15) sensor, Y-bus + to-side back-calc | +| `os1_2msr_no_angle` | `1os2msr-no-angle` | 3 | il, nr | magnitude-only voltages → slack-angle gauge | +| `single_line_load` | `single-line-load-il` | 2 | il | inconsistent data (weighting matters), Kalman-combined duplicate load sensor | + +The single-phase PGM network is the special case of the four-wire literal model +with one conductor per bus and no neutral; per-unit bases match PGM +(`V_base = u_rated`, `S_base = base_power_3p = 1e6`). + +## Regenerate + +```bash +pip install power-grid-model numpy +python3 generate_golden.py /path/to/power-grid-model # or set PGM_REPO +# then paste pgm_cases_data.jl into the PGM_CASES block of ../pgm_cases.jl +``` diff --git a/test/literal/pgm/generate_golden.py b/test/literal/pgm/generate_golden.py new file mode 100644 index 0000000..5c20974 --- /dev/null +++ b/test/literal/pgm/generate_golden.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 +"""Regenerate the PowerGridModel (PGM) benchmark golden data for Literal PMDSE. + +Reads PGM's own state-estimation validation cases, runs PGM, and emits + + golden.json - machine-readable network + sensors + PGM node voltages + pgm_cases_data.jl - the `const PGM_CASES = [...]` block embedded by + test/literal/pgm_cases.jl + +Usage +----- + pip install power-grid-model numpy + python3 generate_golden.py /path/to/power-grid-model + +If the path is omitted, the env var PGM_REPO or ../../../../power-grid-model is +tried. Only the symmetric (single-phase) cases below are exported; they are the +ones the literal estimator is benchmarked against (the four-wire generality is +covered separately by the explicit-neutral feeder tests). +""" +import json, os, sys +import numpy as np +from power_grid_model import PowerGridModel, CalculationMethod +from power_grid_model.utils import json_deserialize + +# (julia name, PGM folder under tests/data/state_estimation, methods) +CASES = [ + ("single_node", "single-node-source-sym-voltage-sensor", ["il", "nr"]), + ("os1_2msr", "1os2msr", ["il", "nr"]), + ("os1_2msr_no_angle", "1os2msr-no-angle", ["il", "nr"]), + ("single_line_load", "single-line-load-il", ["il"]), +] +METH = {"il": CalculationMethod.iterative_linear, "nr": CalculationMethod.newton_raphson} +MTERM = {0: "branch_from", 1: "branch_to", 2: "source", 4: "load", 5: "gen"} +FREQ = 50.0 + + +def pgm_repo(): + if len(sys.argv) > 1: + return sys.argv[1] + here = os.path.dirname(os.path.abspath(__file__)) + return os.environ.get("PGM_REPO", os.path.join(here, "..", "..", "..", "..", "power-grid-model")) + + +def case_dict(base, name, folder, methods): + with open(f"{base}/tests/data/state_estimation/{folder}/input.json") as f: + inp = json_deserialize(f.read()) + d = dict(name=name, folder=folder, methods=methods) + d["nodes"] = [dict(id=int(n["id"]), u_rated=float(n["u_rated"])) for n in inp["node"]] + d["lines"] = [dict(id=int(l["id"]), f=int(l["from_node"]), t=int(l["to_node"]), + r1=float(l["r1"]), x1=float(l["x1"]), c1=float(l["c1"]), tan1=float(l["tan1"]), + fs=int(l["from_status"]), ts=int(l["to_status"])) for l in inp.get("line", [])] + d["source"] = dict(id=int(inp["source"][0]["id"]), node=int(inp["source"][0]["node"])) + d["loads"] = [dict(id=int(l["id"]), node=int(l["node"])) for l in inp.get("sym_load", [])] + appnode = {int(inp["source"][0]["id"]): ("source", int(inp["source"][0]["node"]))} + for l in inp.get("sym_load", []): + appnode[int(l["id"])] = ("load", int(l["node"])) + d["vsensors"] = [] + for s in inp.get("sym_voltage_sensor", []): + ang = float(s["u_angle_measured"]); has = not np.isnan(ang) + d["vsensors"].append(dict(node=int(s["measured_object"]), u=float(s["u_measured"]), + ang=(ang if has else None), sigma=float(s["u_sigma"]))) + d["psensors"] = [] + for s in inp.get("sym_power_sensor", []): + tt = int(s["measured_terminal_type"]); kind = MTERM[tt] + e = dict(kind=kind, p=float(s["p_measured"]), q=float(s["q_measured"]), sigma=float(s["power_sigma"])) + if kind in ("branch_from", "branch_to"): + e["branch"] = int(s["measured_object"]) + else: + e["node"] = appnode[int(s["measured_object"])][1]; e["app"] = appnode[int(s["measured_object"])][0] + d["psensors"].append(e) + model = PowerGridModel(inp, system_frequency=FREQ) + exp = {} + for m in methods: + out = model.calculate_state_estimation(symmetric=True, error_tolerance=1e-12, + max_iterations=100, calculation_method=METH[m]) + nd = out["node"] + nodes = {int(r["id"]): dict(u=float(r["u"]), + u_angle=float(r["u_angle"]) if "u_angle" in nd.dtype.names else 0.0) for r in nd} + lines = {} + if "line" in out: + for r in out["line"]: + lines[int(r["id"])] = dict(p_from=float(r["p_from"]), q_from=float(r["q_from"]), + p_to=float(r["p_to"]), q_to=float(r["q_to"])) + exp[m] = dict(nodes=nodes, lines=lines) + d["expected"] = exp + return d + + +def emit_julia(cases): + jn = lambda x: repr(float(x)) + ja = lambda a: "nothing" if a is None else repr(float(a)) + o = ["# AUTO-GENERATED by test/literal/pgm/generate_golden.py — DO NOT EDIT BY HAND.", + "# Golden state-estimation cases exported from power-grid-model (PGM),", + "# tests/data/state_estimation/*. PGM-computed node voltages per method.", + "", "const PGM_CASES = ["] + for c in cases: + o.append(" (") + o.append(f" name = {json.dumps(c['name'])},") + o.append(f" folder = {json.dumps(c['folder'])},") + o.append(" methods = [" + ", ".join(":" + m for m in c["methods"]) + "],") + o.append(" nodes = [" + ", ".join(f"(id={n['id']}, u_rated={jn(n['u_rated'])})" for n in c["nodes"]) + "],") + o.append(" lines = [" + ", ".join( + f"(id={l['id']}, f={l['f']}, t={l['t']}, r1={jn(l['r1'])}, x1={jn(l['x1'])}, c1={jn(l['c1'])}, tan1={jn(l['tan1'])}, fs={l['fs']}, ts={l['ts']})" + for l in c["lines"]) + "],") + o.append(f" source = (id={c['source']['id']}, node={c['source']['node']}),") + o.append(" loads = [" + ", ".join(f"(id={l['id']}, node={l['node']})" for l in c["loads"]) + "],") + o.append(" vsensors = [" + ", ".join( + f"(node={v['node']}, u={jn(v['u'])}, ang={ja(v['ang'])}, sigma={jn(v['sigma'])})" for v in c["vsensors"]) + "],") + ps = [] + for p in c["psensors"]: + loc = f"branch={p['branch']}" if "branch" in p else f"node={p['node']}, app={json.dumps(p['app'])}" + ps.append(f"(kind={json.dumps(p['kind'])}, {loc}, p={jn(p['p'])}, q={jn(p['q'])}, sigma={jn(p['sigma'])})") + o.append(" psensors = [" + ", ".join(ps) + "],") + o.append(" expected = Dict(") + for m, e in c["expected"].items(): + ne = ", ".join(f"{nid} => (u={jn(v['u'])}, u_angle={jn(v['u_angle'])})" for nid, v in e["nodes"].items()) + le = ", ".join(f"{lid} => (p_from={jn(v['p_from'])}, q_from={jn(v['q_from'])}, p_to={jn(v['p_to'])}, q_to={jn(v['q_to'])})" + for lid, v in e["lines"].items()) + o.append(f" {json.dumps(m)} => (nodes = Dict({ne}), lines = Dict({le})),") + o.append(" ),") + o.append(" ),") + o.append("]") + return "\n".join(o) + "\n" + + +def main(): + base = pgm_repo() + if not os.path.isdir(f"{base}/tests/data/state_estimation"): + sys.exit(f"power-grid-model repo not found at {base!r}; pass its path as argv[1] or set PGM_REPO") + cases = [case_dict(base, *c) for c in CASES] + here = os.path.dirname(os.path.abspath(__file__)) + with open(os.path.join(here, "golden.json"), "w") as f: + json.dump(cases, f, indent=1) + with open(os.path.join(here, "pgm_cases_data.jl"), "w") as f: + f.write(emit_julia(cases)) + print(f"wrote golden.json and pgm_cases_data.jl for {len(cases)} cases from {base}") + + +if __name__ == "__main__": + main() diff --git a/test/literal/pgm/golden.json b/test/literal/pgm/golden.json new file mode 100644 index 0000000..fb42746 --- /dev/null +++ b/test/literal/pgm/golden.json @@ -0,0 +1,538 @@ +[ + { + "name": "single_node", + "folder": "single-node-source-sym-voltage-sensor", + "methods": [ + "il", + "nr" + ], + "nodes": [ + { + "id": 1, + "u_rated": 10000.0 + } + ], + "lines": [], + "source": { + "id": 2, + "node": 1 + }, + "loads": [], + "vsensors": [ + { + "node": 1, + "u": 12345.0, + "ang": 0.1, + "sigma": 100.0 + } + ], + "psensors": [], + "expected": { + "il": { + "nodes": { + "1": { + "u": 12345.0, + "u_angle": 0.1 + } + }, + "lines": {} + }, + "nr": { + "nodes": { + "1": { + "u": 12345.0, + "u_angle": 0.09999999999999999 + } + }, + "lines": {} + } + } + }, + { + "name": "os1_2msr", + "folder": "1os2msr", + "methods": [ + "il", + "nr" + ], + "nodes": [ + { + "id": 1, + "u_rated": 10500.0 + }, + { + "id": 2, + "u_rated": 10500.0 + }, + { + "id": 3, + "u_rated": 10500.0 + } + ], + "lines": [ + { + "id": 4, + "f": 1, + "t": 2, + "r1": 0.11, + "x1": 0.12, + "c1": 4.1380285203892784e-05, + "tan1": 0.1076923076923077, + "fs": 1, + "ts": 1 + }, + { + "id": 5, + "f": 2, + "t": 3, + "r1": 0.15, + "x1": 0.16, + "c1": 5.411268065124442e-05, + "tan1": 0.10588235294117646, + "fs": 1, + "ts": 1 + } + ], + "source": { + "id": 6, + "node": 1 + }, + "loads": [ + { + "id": 7, + "node": 2 + }, + { + "id": 8, + "node": 3 + } + ], + "vsensors": [ + { + "node": 1, + "u": 10751.072595758282, + "ang": -0.013054638926306409, + "sigma": 105.0 + }, + { + "node": 2, + "u": 10752.698591183394, + "ang": -0.01763734945972652, + "sigma": 105.0 + }, + { + "node": 3, + "u": 10748.320749959701, + "ang": -0.02018233075947404, + "sigma": 100.0 + } + ], + "psensors": [ + { + "kind": "branch_from", + "p": 2412359.2976399013, + "q": -3024028.886598367, + "sigma": 37916.0, + "branch": 4 + }, + { + "kind": "branch_to", + "p": -2000000.0, + "q": 1000000.0, + "sigma": 1000000000000000.0, + "branch": 4 + }, + { + "kind": "branch_from", + "p": 1230426.390004009, + "q": -1742195.1033582848, + "sigma": 20878.0, + "branch": 5 + }, + { + "kind": "branch_to", + "p": -1019999.9999999485, + "q": -219999.99999999927, + "sigma": 10435.0, + "branch": 5 + }, + { + "kind": "source", + "p": 2412359.297639887, + "q": -3024028.8865982923, + "sigma": 38009.0, + "node": 1, + "app": "source" + }, + { + "kind": "load", + "p": 1010000.0, + "q": 210000.0, + "sigma": 10316.0, + "node": 2, + "app": "load" + }, + { + "kind": "load", + "p": 1020000.0, + "q": 220000.0, + "sigma": 10435.0, + "node": 3, + "app": "load" + } + ], + "expected": { + "il": { + "nodes": { + "1": { + "u": 10751.072595669397, + "u_angle": -0.01305463892705287 + }, + "2": { + "u": 10752.698591094362, + "u_angle": -0.017637349460495298 + }, + "3": { + "u": 10748.32074987062, + "u_angle": -0.020182330760257754 + } + }, + "lines": { + "4": { + "p_from": 2412359.297617528, + "q_from": -3024028.8865525783, + "p_to": -2240426.3899840107, + "q_to": 1532195.1033373913 + }, + "5": { + "p_from": 1230426.389992248, + "q_from": -1742195.1033316434, + "p_to": -1019999.9999917165, + "q_to": -219999.99999405813 + } + } + }, + "nr": { + "nodes": { + "1": { + "u": 10751.072595758176, + "u_angle": -0.013054638926306369 + }, + "2": { + "u": 10752.698591183289, + "u_angle": -0.01763734945972652 + }, + "3": { + "u": 10748.3207499596, + "u_angle": -0.020182330759474076 + } + }, + "lines": { + "4": { + "p_from": 2412359.2976397616, + "q_from": -3024028.886598335, + "p_to": -2240426.390003599, + "q_to": 1532195.1033583526 + }, + "5": { + "p_from": 1230426.3900039378, + "q_from": -1742195.1033583248, + "p_to": -1019999.9999998793, + "q_to": -219999.99999994013 + } + } + } + } + }, + { + "name": "os1_2msr_no_angle", + "folder": "1os2msr-no-angle", + "methods": [ + "il", + "nr" + ], + "nodes": [ + { + "id": 1, + "u_rated": 10500.0 + }, + { + "id": 2, + "u_rated": 10500.0 + }, + { + "id": 3, + "u_rated": 10500.0 + } + ], + "lines": [ + { + "id": 4, + "f": 1, + "t": 2, + "r1": 0.11, + "x1": 0.12, + "c1": 4.1380285203892784e-05, + "tan1": 0.1076923076923077, + "fs": 1, + "ts": 1 + }, + { + "id": 5, + "f": 2, + "t": 3, + "r1": 0.15, + "x1": 0.16, + "c1": 5.411268065124442e-05, + "tan1": 0.10588235294117646, + "fs": 1, + "ts": 1 + } + ], + "source": { + "id": 6, + "node": 1 + }, + "loads": [ + { + "id": 7, + "node": 2 + }, + { + "id": 8, + "node": 3 + } + ], + "vsensors": [ + { + "node": 1, + "u": 10751.072595758282, + "ang": null, + "sigma": 105.0 + }, + { + "node": 2, + "u": 10752.698591183394, + "ang": null, + "sigma": 105.0 + }, + { + "node": 3, + "u": 10748.320749959701, + "ang": null, + "sigma": 100.0 + } + ], + "psensors": [ + { + "kind": "branch_from", + "p": 2412359.2976399013, + "q": -3024028.886598367, + "sigma": 37916.0, + "branch": 4 + }, + { + "kind": "branch_to", + "p": -2240426.3900036793, + "q": 1532195.1033583684, + "sigma": 26809.0, + "branch": 4 + }, + { + "kind": "branch_from", + "p": 1230426.390004009, + "q": -1742195.1033582848, + "sigma": 20878.0, + "branch": 5 + }, + { + "kind": "branch_to", + "p": -1019999.9999999485, + "q": -219999.99999999927, + "sigma": 10435.0, + "branch": 5 + }, + { + "kind": "source", + "p": 2412359.297639887, + "q": -3024028.8865982923, + "sigma": 38009.0, + "node": 1, + "app": "source" + }, + { + "kind": "load", + "p": 1010000.0, + "q": 210000.0, + "sigma": 10316.0, + "node": 2, + "app": "load" + }, + { + "kind": "load", + "p": 1020000.0, + "q": 220000.0, + "sigma": 10435.0, + "node": 3, + "app": "load" + } + ], + "expected": { + "il": { + "nodes": { + "1": { + "u": 10751.072595583171, + "u_angle": 0.0 + }, + "2": { + "u": 10752.698591008038, + "u_angle": -0.004582710533466773 + }, + "3": { + "u": 10748.32074978426, + "u_angle": -0.007127691833244886 + } + }, + "lines": { + "4": { + "p_from": 2412359.2975955657, + "q_from": -3024028.886511328, + "p_to": -2240426.3899647417, + "q_to": 1532195.103320241 + }, + "5": { + "p_from": 1230426.3899807737, + "q_from": -1742195.1033068383, + "p_to": -1019999.9999835168, + "q_to": -219999.99998724568 + } + } + }, + "nr": { + "nodes": { + "1": { + "u": 10751.072595758145, + "u_angle": 0.0 + }, + "2": { + "u": 10752.69859118326, + "u_angle": -0.0045827105334202535 + }, + "3": { + "u": 10748.32074995957, + "u_angle": -0.007127691833167822 + } + }, + "lines": { + "4": { + "p_from": 2412359.297639833, + "q_from": -3024028.886598318, + "p_to": -2240426.390003647, + "q_to": 1532195.1033584457 + }, + "5": { + "p_from": 1230426.3900038258, + "q_from": -1742195.1033584112, + "p_to": -1019999.9999998506, + "q_to": -219999.99999981755 + } + } + } + } + }, + { + "name": "single_line_load", + "folder": "single-line-load-il", + "methods": [ + "il" + ], + "nodes": [ + { + "id": 11, + "u_rated": 10500.0 + }, + { + "id": 12, + "u_rated": 10500.0 + } + ], + "lines": [ + { + "id": 21, + "f": 11, + "t": 12, + "r1": 0.1, + "x1": 0.0, + "c1": 0.0, + "tan1": 0.0, + "fs": 1, + "ts": 1 + } + ], + "source": { + "id": 41, + "node": 11 + }, + "loads": [ + { + "id": 31, + "node": 12 + } + ], + "vsensors": [ + { + "node": 11, + "u": 10500.0, + "ang": null, + "sigma": 105.0 + } + ], + "psensors": [ + { + "kind": "branch_from", + "p": 1000000.0, + "q": 0.0, + "sigma": 1000.0, + "branch": 21 + }, + { + "kind": "load", + "p": 2000000.0, + "q": 0.0, + "sigma": 1000.0, + "node": 12, + "app": "load" + }, + { + "kind": "load", + "p": 2000000.0, + "q": 0.0, + "sigma": 1e+200, + "node": 12, + "app": "load" + } + ], + "expected": { + "il": { + "nodes": { + "11": { + "u": 10499.999990211563, + "u_angle": 0.0 + }, + "12": { + "u": 10485.70128888785, + "u_angle": 0.0 + } + }, + "lines": { + "21": { + "p_from": 1501363.6375900682, + "q_from": 0.0, + "p_to": -1499319.1089946209, + "q_to": -0.0 + } + } + } + } + } +] \ No newline at end of file diff --git a/test/literal/pgm/pgm_cases_data.jl b/test/literal/pgm/pgm_cases_data.jl new file mode 100644 index 0000000..0ef3a49 --- /dev/null +++ b/test/literal/pgm/pgm_cases_data.jl @@ -0,0 +1,65 @@ +# AUTO-GENERATED by test/literal/pgm/generate_golden.py — DO NOT EDIT BY HAND. +# Golden state-estimation cases exported from power-grid-model (PGM), +# tests/data/state_estimation/*. PGM-computed node voltages per method. + +const PGM_CASES = [ + ( + name = "single_node", + folder = "single-node-source-sym-voltage-sensor", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10000.0)], + lines = [], + source = (id=2, node=1), + loads = [], + vsensors = [(node=1, u=12345.0, ang=0.1, sigma=100.0)], + psensors = [], + expected = Dict( + "il" => (nodes = Dict(1 => (u=12345.0, u_angle=0.1)), lines = Dict()), + "nr" => (nodes = Dict(1 => (u=12345.0, u_angle=0.09999999999999999)), lines = Dict()), + ), + ), + ( + name = "os1_2msr", + folder = "1os2msr", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10500.0), (id=2, u_rated=10500.0), (id=3, u_rated=10500.0)], + lines = [(id=4, f=1, t=2, r1=0.11, x1=0.12, c1=4.1380285203892784e-05, tan1=0.1076923076923077, fs=1, ts=1), (id=5, f=2, t=3, r1=0.15, x1=0.16, c1=5.411268065124442e-05, tan1=0.10588235294117646, fs=1, ts=1)], + source = (id=6, node=1), + loads = [(id=7, node=2), (id=8, node=3)], + vsensors = [(node=1, u=10751.072595758282, ang=-0.013054638926306409, sigma=105.0), (node=2, u=10752.698591183394, ang=-0.01763734945972652, sigma=105.0), (node=3, u=10748.320749959701, ang=-0.02018233075947404, sigma=100.0)], + psensors = [(kind="branch_from", branch=4, p=2412359.2976399013, q=-3024028.886598367, sigma=37916.0), (kind="branch_to", branch=4, p=-2000000.0, q=1000000.0, sigma=1000000000000000.0), (kind="branch_from", branch=5, p=1230426.390004009, q=-1742195.1033582848, sigma=20878.0), (kind="branch_to", branch=5, p=-1019999.9999999485, q=-219999.99999999927, sigma=10435.0), (kind="source", node=1, app="source", p=2412359.297639887, q=-3024028.8865982923, sigma=38009.0), (kind="load", node=2, app="load", p=1010000.0, q=210000.0, sigma=10316.0), (kind="load", node=3, app="load", p=1020000.0, q=220000.0, sigma=10435.0)], + expected = Dict( + "il" => (nodes = Dict(1 => (u=10751.072595669397, u_angle=-0.01305463892705287), 2 => (u=10752.698591094362, u_angle=-0.017637349460495298), 3 => (u=10748.32074987062, u_angle=-0.020182330760257754)), lines = Dict(4 => (p_from=2412359.297617528, q_from=-3024028.8865525783, p_to=-2240426.3899840107, q_to=1532195.1033373913), 5 => (p_from=1230426.389992248, q_from=-1742195.1033316434, p_to=-1019999.9999917165, q_to=-219999.99999405813))), + "nr" => (nodes = Dict(1 => (u=10751.072595758176, u_angle=-0.013054638926306369), 2 => (u=10752.698591183289, u_angle=-0.01763734945972652), 3 => (u=10748.3207499596, u_angle=-0.020182330759474076)), lines = Dict(4 => (p_from=2412359.2976397616, q_from=-3024028.886598335, p_to=-2240426.390003599, q_to=1532195.1033583526), 5 => (p_from=1230426.3900039378, q_from=-1742195.1033583248, p_to=-1019999.9999998793, q_to=-219999.99999994013))), + ), + ), + ( + name = "os1_2msr_no_angle", + folder = "1os2msr-no-angle", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10500.0), (id=2, u_rated=10500.0), (id=3, u_rated=10500.0)], + lines = [(id=4, f=1, t=2, r1=0.11, x1=0.12, c1=4.1380285203892784e-05, tan1=0.1076923076923077, fs=1, ts=1), (id=5, f=2, t=3, r1=0.15, x1=0.16, c1=5.411268065124442e-05, tan1=0.10588235294117646, fs=1, ts=1)], + source = (id=6, node=1), + loads = [(id=7, node=2), (id=8, node=3)], + vsensors = [(node=1, u=10751.072595758282, ang=nothing, sigma=105.0), (node=2, u=10752.698591183394, ang=nothing, sigma=105.0), (node=3, u=10748.320749959701, ang=nothing, sigma=100.0)], + psensors = [(kind="branch_from", branch=4, p=2412359.2976399013, q=-3024028.886598367, sigma=37916.0), (kind="branch_to", branch=4, p=-2240426.3900036793, q=1532195.1033583684, sigma=26809.0), (kind="branch_from", branch=5, p=1230426.390004009, q=-1742195.1033582848, sigma=20878.0), (kind="branch_to", branch=5, p=-1019999.9999999485, q=-219999.99999999927, sigma=10435.0), (kind="source", node=1, app="source", p=2412359.297639887, q=-3024028.8865982923, sigma=38009.0), (kind="load", node=2, app="load", p=1010000.0, q=210000.0, sigma=10316.0), (kind="load", node=3, app="load", p=1020000.0, q=220000.0, sigma=10435.0)], + expected = Dict( + "il" => (nodes = Dict(1 => (u=10751.072595583171, u_angle=0.0), 2 => (u=10752.698591008038, u_angle=-0.004582710533466773), 3 => (u=10748.32074978426, u_angle=-0.007127691833244886)), lines = Dict(4 => (p_from=2412359.2975955657, q_from=-3024028.886511328, p_to=-2240426.3899647417, q_to=1532195.103320241), 5 => (p_from=1230426.3899807737, q_from=-1742195.1033068383, p_to=-1019999.9999835168, q_to=-219999.99998724568))), + "nr" => (nodes = Dict(1 => (u=10751.072595758145, u_angle=0.0), 2 => (u=10752.69859118326, u_angle=-0.0045827105334202535), 3 => (u=10748.32074995957, u_angle=-0.007127691833167822)), lines = Dict(4 => (p_from=2412359.297639833, q_from=-3024028.886598318, p_to=-2240426.390003647, q_to=1532195.1033584457), 5 => (p_from=1230426.3900038258, q_from=-1742195.1033584112, p_to=-1019999.9999998506, q_to=-219999.99999981755))), + ), + ), + ( + name = "single_line_load", + folder = "single-line-load-il", + methods = [:il], + nodes = [(id=11, u_rated=10500.0), (id=12, u_rated=10500.0)], + lines = [(id=21, f=11, t=12, r1=0.1, x1=0.0, c1=0.0, tan1=0.0, fs=1, ts=1)], + source = (id=41, node=11), + loads = [(id=31, node=12)], + vsensors = [(node=11, u=10500.0, ang=nothing, sigma=105.0)], + psensors = [(kind="branch_from", branch=21, p=1000000.0, q=0.0, sigma=1000.0), (kind="load", node=12, app="load", p=2000000.0, q=0.0, sigma=1000.0), (kind="load", node=12, app="load", p=2000000.0, q=0.0, sigma=1e+200)], + expected = Dict( + "il" => (nodes = Dict(11 => (u=10499.999990211563, u_angle=0.0), 12 => (u=10485.70128888785, u_angle=0.0)), lines = Dict(21 => (p_from=1501363.6375900682, q_from=0.0, p_to=-1499319.1089946209, q_to=-0.0))), + ), + ), +] diff --git a/test/literal/pgm_cases.jl b/test/literal/pgm_cases.jl new file mode 100644 index 0000000..8d80439 --- /dev/null +++ b/test/literal/pgm_cases.jl @@ -0,0 +1,195 @@ +################################################################################ +# PowerGridModel (PGM) benchmark cases for the Literal PMDSE PGM solve options. # +# # +# Each entry is a PGM symmetric (single-phase positive-sequence) state- # +# estimation case taken verbatim from # +# power-grid-model/tests/data/state_estimation// # +# together with the node voltages PGM computes for each calculation method. # +# `build_pgm_math` turns the PGM network into the *same* PMD mathematical data # +# dictionary the package consumes (per-unit, V_base = u_rated, S_base = 1e6 = # +# PGM `base_power_3p`); `build_pgm_meas!` writes `data["meas"]` using the # +# package measurement vocabulary. The single-phase PGM network is the special # +# case of the four-wire literal model with one conductor per bus and no neutral.# +# # +# Regenerate the golden numbers (PGM_CASES below) with # +# python3 test/literal/pgm/generate_golden.py (needs PGM + numpy)# +# then paste test/literal/pgm/pgm_cases_data.jl in place of the block below. # +################################################################################ + +import Distributions as _DST +import PowerModelsDistribution as _PMD + +const PGM_FREQ = 50.0 +const PGM_SBASE = 1.0e6 # PGM common/base_power_3p (symmetric calculation) + +# AUTO-GENERATED by test/literal/pgm/generate_golden.py — DO NOT EDIT BY HAND. +# Golden state-estimation cases exported from power-grid-model (PGM) v1.12, +# tests/data/state_estimation/*. Each case carries the PGM network + sensors and +# the PGM-computed node voltages (per calculation method) used as the benchmark. + +const PGM_CASES = [ + ( + name = "single_node", + folder = "single-node-source-sym-voltage-sensor", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10000.0)], + lines = [], + source = (id=2, node=1), + loads = [], + vsensors = [(node=1, u=12345.0, ang=0.1, sigma=100.0)], + psensors = [], + expected = Dict( + "il" => (nodes = Dict(1 => (u=12345.0, u_angle=0.1)), lines = Dict()), + "nr" => (nodes = Dict(1 => (u=12345.0, u_angle=0.09999999999999999)), lines = Dict()), + ), + ), + ( + name = "os1_2msr", + folder = "1os2msr", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10500.0), (id=2, u_rated=10500.0), (id=3, u_rated=10500.0)], + lines = [(id=4, f=1, t=2, r1=0.11, x1=0.12, c1=4.1380285203892784e-05, tan1=0.1076923076923077, fs=1, ts=1), (id=5, f=2, t=3, r1=0.15, x1=0.16, c1=5.411268065124442e-05, tan1=0.10588235294117646, fs=1, ts=1)], + source = (id=6, node=1), + loads = [(id=7, node=2), (id=8, node=3)], + vsensors = [(node=1, u=10751.072595758282, ang=-0.013054638926306409, sigma=105.0), (node=2, u=10752.698591183394, ang=-0.01763734945972652, sigma=105.0), (node=3, u=10748.320749959701, ang=-0.02018233075947404, sigma=100.0)], + psensors = [(kind="branch_from", branch=4, p=2412359.2976399013, q=-3024028.886598367, sigma=37916.0), (kind="branch_to", branch=4, p=-2000000.0, q=1000000.0, sigma=1000000000000000.0), (kind="branch_from", branch=5, p=1230426.390004009, q=-1742195.1033582848, sigma=20878.0), (kind="branch_to", branch=5, p=-1019999.9999999485, q=-219999.99999999927, sigma=10435.0), (kind="source", node=1, app="source", p=2412359.297639887, q=-3024028.8865982923, sigma=38009.0), (kind="load", node=2, app="load", p=1010000.0, q=210000.0, sigma=10316.0), (kind="load", node=3, app="load", p=1020000.0, q=220000.0, sigma=10435.0)], + expected = Dict( + "il" => (nodes = Dict(1 => (u=10751.072595669397, u_angle=-0.01305463892705287), 2 => (u=10752.698591094362, u_angle=-0.017637349460495298), 3 => (u=10748.32074987062, u_angle=-0.020182330760257754)), lines = Dict(4 => (p_from=2412359.297617528, q_from=-3024028.8865525783, p_to=-2240426.3899840107, q_to=1532195.1033373913), 5 => (p_from=1230426.389992248, q_from=-1742195.1033316434, p_to=-1019999.9999917165, q_to=-219999.99999405813))), + "nr" => (nodes = Dict(1 => (u=10751.072595758176, u_angle=-0.013054638926306369), 2 => (u=10752.698591183289, u_angle=-0.01763734945972652), 3 => (u=10748.3207499596, u_angle=-0.020182330759474076)), lines = Dict(4 => (p_from=2412359.2976397616, q_from=-3024028.886598335, p_to=-2240426.390003599, q_to=1532195.1033583526), 5 => (p_from=1230426.3900039378, q_from=-1742195.1033583248, p_to=-1019999.9999998793, q_to=-219999.99999994013))), + ), + ), + ( + name = "os1_2msr_no_angle", + folder = "1os2msr-no-angle", + methods = [:il, :nr], + nodes = [(id=1, u_rated=10500.0), (id=2, u_rated=10500.0), (id=3, u_rated=10500.0)], + lines = [(id=4, f=1, t=2, r1=0.11, x1=0.12, c1=4.1380285203892784e-05, tan1=0.1076923076923077, fs=1, ts=1), (id=5, f=2, t=3, r1=0.15, x1=0.16, c1=5.411268065124442e-05, tan1=0.10588235294117646, fs=1, ts=1)], + source = (id=6, node=1), + loads = [(id=7, node=2), (id=8, node=3)], + vsensors = [(node=1, u=10751.072595758282, ang=nothing, sigma=105.0), (node=2, u=10752.698591183394, ang=nothing, sigma=105.0), (node=3, u=10748.320749959701, ang=nothing, sigma=100.0)], + psensors = [(kind="branch_from", branch=4, p=2412359.2976399013, q=-3024028.886598367, sigma=37916.0), (kind="branch_to", branch=4, p=-2240426.3900036793, q=1532195.1033583684, sigma=26809.0), (kind="branch_from", branch=5, p=1230426.390004009, q=-1742195.1033582848, sigma=20878.0), (kind="branch_to", branch=5, p=-1019999.9999999485, q=-219999.99999999927, sigma=10435.0), (kind="source", node=1, app="source", p=2412359.297639887, q=-3024028.8865982923, sigma=38009.0), (kind="load", node=2, app="load", p=1010000.0, q=210000.0, sigma=10316.0), (kind="load", node=3, app="load", p=1020000.0, q=220000.0, sigma=10435.0)], + expected = Dict( + "il" => (nodes = Dict(1 => (u=10751.072595583171, u_angle=0.0), 2 => (u=10752.698591008038, u_angle=-0.004582710533466773), 3 => (u=10748.32074978426, u_angle=-0.007127691833244886)), lines = Dict(4 => (p_from=2412359.2975955657, q_from=-3024028.886511328, p_to=-2240426.3899647417, q_to=1532195.103320241), 5 => (p_from=1230426.3899807737, q_from=-1742195.1033068383, p_to=-1019999.9999835168, q_to=-219999.99998724568))), + "nr" => (nodes = Dict(1 => (u=10751.072595758145, u_angle=0.0), 2 => (u=10752.69859118326, u_angle=-0.0045827105334202535), 3 => (u=10748.32074995957, u_angle=-0.007127691833167822)), lines = Dict(4 => (p_from=2412359.297639833, q_from=-3024028.886598318, p_to=-2240426.390003647, q_to=1532195.1033584457), 5 => (p_from=1230426.3900038258, q_from=-1742195.1033584112, p_to=-1019999.9999998506, q_to=-219999.99999981755))), + ), + ), + ( + name = "single_line_load", + folder = "single-line-load-il", + methods = [:il], + nodes = [(id=11, u_rated=10500.0), (id=12, u_rated=10500.0)], + lines = [(id=21, f=11, t=12, r1=0.1, x1=0.0, c1=0.0, tan1=0.0, fs=1, ts=1)], + source = (id=41, node=11), + loads = [(id=31, node=12)], + vsensors = [(node=11, u=10500.0, ang=nothing, sigma=105.0)], + psensors = [(kind="branch_from", branch=21, p=1000000.0, q=0.0, sigma=1000.0), (kind="load", node=12, app="load", p=2000000.0, q=0.0, sigma=1000.0), (kind="load", node=12, app="load", p=2000000.0, q=0.0, sigma=1e+200)], + expected = Dict( + "il" => (nodes = Dict(11 => (u=10499.999990211563, u_angle=0.0), 12 => (u=10485.70128888785, u_angle=0.0)), lines = Dict(21 => (p_from=1501363.6375900682, q_from=0.0, p_to=-1499319.1089946209, q_to=-0.0))), + ), + ), +] + +# --------------------------------------------------------------------------- +# PGM network -> PMD mathematical data dictionary (per-unit) +# --------------------------------------------------------------------------- + +"per-unit PMD `data_math` dict equivalent to a PGM symmetric case (1 conductor)" +function build_pgm_math(case) + sbase = PGM_SBASE + vbase = Dict(n.id => n.u_rated for n in case.nodes) + ω = 2π * PGM_FREQ + busd = Dict{String,Any}(); brd = Dict{String,Any}() + loadd = Dict{String,Any}(); gend = Dict{String,Any}() + + for n in case.nodes + bt = n.id == case.source.node ? 3 : 1 + busd[string(n.id)] = Dict{String,Any}( + "index" => n.id, "bus_i" => n.id, "bus_type" => bt, + "terminals" => [1], "grounded" => Bool[false], + "vbase" => n.u_rated, "vmin" => [0.0], "vmax" => [Inf], + "vr_start" => [1.0], "vi_start" => [0.0]) + end + for l in case.lines + zb = vbase[l.f]^2 / sbase # both ends share u_rated here + b1 = ω * l.c1; g1 = l.tan1 * b1 + gsh = (g1 * zb) / 2; bsh = (b1 * zb) / 2 # π-model, half each end + brd[string(l.id)] = Dict{String,Any}( + "index" => l.id, "f_bus" => l.f, "t_bus" => l.t, + "f_connections" => [1], "t_connections" => [1], + "br_r" => reshape([l.r1 / zb], 1, 1), "br_x" => reshape([l.x1 / zb], 1, 1), + "g_fr" => reshape([gsh], 1, 1), "b_fr" => reshape([bsh], 1, 1), + "g_to" => reshape([gsh], 1, 1), "b_to" => reshape([bsh], 1, 1), + "br_status" => (l.fs == 1 && l.ts == 1) ? 1 : 0) + end + # PGM source acts like a generator (defines where injection is allowed) + gend[string(case.source.id)] = Dict{String,Any}( + "index" => case.source.id, "gen_bus" => case.source.node, + "connections" => [1], "configuration" => _PMD.WYE) + for ld in case.loads + loadd[string(ld.id)] = Dict{String,Any}( + "index" => ld.id, "load_bus" => ld.node, "connections" => [1], + "configuration" => _PMD.WYE, "pd" => [0.0], "qd" => [0.0], + "model" => _PMD.POWER, "status" => 1) + end + return Dict{String,Any}("bus" => busd, "branch" => brd, "load" => loadd, "gen" => gend, + "shunt" => Dict{String,Any}(), "settings" => Dict{String,Any}("sbase" => sbase), + "per_unit" => true) +end + +"map an appliance node id to its PGM source/load math id for this case" +function _pgm_app_id(case, node, app) + app == "source" && return case.source.id + for ld in case.loads + ld.node == node && return ld.id + end + error("PGM case $(case.name): no $(app) appliance at node $(node)") +end + +"write `data[\"meas\"]` (per-unit) from the PGM sensors, in package vocabulary" +function build_pgm_meas!(math, case) + sbase = PGM_SBASE + vbase = Dict(n.id => n.u_rated for n in case.nodes) + M = Dict{String,Any}(); k = 0 + add(d) = (k += 1; M[string(k)] = d) + N(μ, σ) = _DST.Normal(μ, σ) + + for v in case.vsensors + vb = vbase[v.node] + if v.ang === nothing # magnitude only + add(Dict{String,Any}("cmp" => :bus, "cmp_id" => v.node, "var" => :vm, + "dst" => [N(v.u / vb, v.sigma / vb)])) + else # voltage phasor -> vr & vi + add(Dict{String,Any}("cmp" => :bus, "cmp_id" => v.node, "var" => :vr, + "dst" => [N(v.u * cos(v.ang) / vb, v.sigma / vb)])) + add(Dict{String,Any}("cmp" => :bus, "cmp_id" => v.node, "var" => :vi, + "dst" => [N(v.u * sin(v.ang) / vb, v.sigma / vb)])) + end + end + for p in case.psensors + if p.kind in ("branch_from", "branch_to") + side = p.kind == "branch_from" ? :from : :to + add(Dict{String,Any}("cmp" => :branch, "cmp_id" => p.branch, "var" => :p, "side" => side, + "dst" => [N(p.p / sbase, p.sigma / sbase)])) + add(Dict{String,Any}("cmp" => :branch, "cmp_id" => p.branch, "var" => :q, "side" => side, + "dst" => [N(p.q / sbase, p.sigma / sbase)])) + else # source / load / gen injection + isgen = p.app in ("source", "gen") + cmp = isgen ? :gen : :load + aid = _pgm_app_id(case, p.node, p.app) + pv, qv = isgen ? (:pg, :qg) : (:pd, :qd) + add(Dict{String,Any}("cmp" => cmp, "cmp_id" => aid, "var" => pv, + "dst" => [N(p.p / sbase, p.sigma / sbase)])) + add(Dict{String,Any}("cmp" => cmp, "cmp_id" => aid, "var" => qv, + "dst" => [N(p.q / sbase, p.sigma / sbase)])) + end + end + math["meas"] = M + return math +end + +"complex node voltage (SI volts) of a literal result for PGM node id `nid`" +function pgm_result_voltage(case, res, nid) + b = res.solution["bus"][string(nid)] + vbase = first(n.u_rated for n in case.nodes if n.id == nid) + return (b["vr"][1] + im * b["vi"][1]) * vbase +end diff --git a/test/literal/test_pgm_benchmark.jl b/test/literal/test_pgm_benchmark.jl new file mode 100644 index 0000000..34a58c9 --- /dev/null +++ b/test/literal/test_pgm_benchmark.jl @@ -0,0 +1,93 @@ +################################################################################ +# Gate 6 : validate the Literal PMDSE PowerGridModel solve options # +# (`:iterative_linear`, `:newton_raphson`) against PowerGridModel itself. # +# # +# Every case in `pgm_cases.jl` is a PGM state-estimation example whose node # +# voltages were computed by PGM (the golden `expected` field). Here the same # +# network + sensors are fed to the literal estimator and the converged node # +# voltages (and line flows, hence the Y-bus and to-side branch back-calc) are # +# required to match PGM to tolerance — for *both* PGM methods, on the same # +# four-wire-general code path that also handles the explicit-neutral feeders. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Distributions as _DST +import LinearAlgebra +include("pgm_cases.jl") + +const _PGM_METHOD = Dict(:il => :iterative_linear, :nr => :newton_raphson) + +"line (p_from,q_from,p_to,q_to) in SI from a literal solution, via the π-model" +function pgm_line_flows(case, res) + sbase = PGM_SBASE; ω = 2π * PGM_FREQ + vbase = Dict(n.id => n.u_rated for n in case.nodes) + U(nid) = (b = res.solution["bus"][string(nid)]; (b["vr"][1] + im * b["vi"][1])) # p.u. + flows = Dict{Int,NamedTuple}() + for l in case.lines + zb = vbase[l.f]^2 / sbase + ys = 1.0 / ((l.r1 + im * l.x1) / zb) + b1 = ω * l.c1; g1 = l.tan1 * b1; ysh = ((g1 + im * b1) * zb) / 2 + Uf = U(l.f); Ut = U(l.t) + Ifr = (ys + ysh) * Uf - ys * Ut + Ito = (ys + ysh) * Ut - ys * Uf + Sf = Uf * conj(Ifr) * sbase; St = Ut * conj(Ito) * sbase + flows[l.id] = (p_from = real(Sf), q_from = imag(Sf), p_to = real(St), q_to = imag(St)) + end + return flows +end + +@testset "Literal PMDSE — PowerGridModel benchmark (gate 6)" begin + for case in PGM_CASES + @testset "$(case.name)" begin + for m in case.methods + math = build_pgm_math(case) + build_pgm_meas!(math, case) + res = _PMDSE.solve_mc_se_literal(math; estimator = :wls, + method = _PGM_METHOD[m], reference = :prop, maxiter = 100, tol = 1e-10) + @test res.termination == :converged + + ex = case.expected[string(m)] + maxe = 0.0 + for (nid, ev) in ex.nodes + Uest = pgm_result_voltage(case, res, nid) + Ugold = ev.u * exp(im * ev.u_angle) + maxe = max(maxe, abs(Uest - Ugold)) + @test abs(Uest - Ugold) < 1e-3 # SI volts (~1e-7 typical) + @test isapprox(abs(Uest), ev.u; rtol = 1e-6) + end + # Y-bus + to-side back-calculation: line flows must match PGM + if !isempty(case.lines) + fl = pgm_line_flows(case, res) + for (lid, ev) in ex.lines + f = fl[lid] + @test isapprox(f.p_from, ev.p_from; rtol = 1e-5, atol = 1.0) + @test isapprox(f.q_from, ev.q_from; rtol = 1e-5, atol = 1.0) + @test isapprox(f.p_to, ev.p_to; rtol = 1e-5, atol = 1.0) + @test isapprox(f.q_to, ev.q_to; rtol = 1e-5, atol = 1.0) + end + end + @info "PGM benchmark" case = case.name method = _PGM_METHOD[m] max_U_err_V = maxe iters = res.iterations + end + end + end + + # The two PGM methods must converge to the same WLS optimum (PGM ships both) + @testset "iterative_linear == newton_raphson" begin + for case in PGM_CASES + (:il in case.methods && :nr in case.methods) || continue + math = build_pgm_math(case); build_pgm_meas!(math, case) + il = _PMDSE.solve_mc_se_literal(math; estimator = :wls, method = :iterative_linear, + reference = :prop, maxiter = 100, tol = 1e-11) + math2 = build_pgm_math(case); build_pgm_meas!(math2, case) + nr = _PMDSE.solve_mc_se_literal(math2; estimator = :wls, method = :newton_raphson, + reference = :prop, maxiter = 100, tol = 1e-11) + d = 0.0 + for n in case.nodes + d = max(d, abs(pgm_result_voltage(case, il, n.id) - + pgm_result_voltage(case, nr, n.id))) + end + @test d < 1e-2 # same state to well under a volt on a ~10 kV system + end + end +end diff --git a/test/literal/test_pgm_en_generality.jl b/test/literal/test_pgm_en_generality.jl new file mode 100644 index 0000000..ab97c49 --- /dev/null +++ b/test/literal/test_pgm_en_generality.jl @@ -0,0 +1,47 @@ +################################################################################ +# Gate 7 : the PowerGridModel solve options are GENERAL — they run unchanged on # +# the four-wire, explicit-neutral feeder, not only on PGM's Kron-reduced # +# single-phase networks. # +# # +# The same `:iterative_linear` / `:newton_raphson` code path is fed the # +# canonical explicit-neutral IVREN measurement set (rectangular voltages + # +# per-conductor current injections, neutral kept explicit) on the 3-bus 4-wire # +# feeder. It must (a) converge, (b) recover the IVREN power-flow state, and # +# (c) agree with the gate-4-validated `solve_wls`. This exercises the neutral # +# handling (phase-to-neutral voltages, neutral injection `= −Σ phase currents`) # +# that PGM — being single-phase — never needs. # +################################################################################ +using Test +import PowerModelsDistribution as _PMD +import PowerModelsDistributionStateEstimation as _PMDSE +import Ipopt, ForwardDiff, LinearAlgebra +import Distributions as _DST +include("helpers.jl") + +@testset "Literal PMDSE — four-wire EN generality of PGM solve options (gate 7)" begin + math, pf = load_en_feeder("3bus_4wire.dss") + canonical_meas!(math, pf; σ = 0.005) # noiseless EN current/voltage set + math["se_settings"] = Dict{String,Any}("criterion" => "wls", "rescaler" => 1.0) + refvals = ref_values_from_pf(math, pf) + + base = _PMDSE.solve_mc_se_literal(math; estimator = :wls, reference = :full_slack, ref_values = refvals) + @test base.termination == :converged + + for method in (:iterative_linear, :newton_raphson) + res = _PMDSE.solve_mc_se_literal(math; estimator = :wls, method = method, + reference = :full_slack, ref_values = refvals, + maxiter = 200, tol = 1e-9) + @test res.termination == :converged + # (b) recovers the explicit-neutral PF state, neutral included + @test _PMDSE.accuracy_metrics(res, pf["solution"]; include_neutral = true).maxerr < 1e-4 + # (c) agrees with the gate-4-validated Gauss–Newton solve_wls + agree = 0.0 + for (b, bsol) in base.solution["bus"], idx in eachindex(bsol["vr"]) + agree = max(agree, abs((res.solution["bus"][b]["vr"][idx] + im * res.solution["bus"][b]["vi"][idx]) - + (bsol["vr"][idx] + im * bsol["vi"][idx]))) + end + @test agree < 1e-4 + @info "EN generality" method = method maxerr_vs_PF = + _PMDSE.accuracy_metrics(res, pf["solution"]).maxerr agree_vs_wls = agree iters = res.iterations + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 855ba7f..90cf872 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,6 +58,8 @@ 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 +@safetestset "Literal PMDSE — PGM benchmark" begin include("literal/test_pgm_benchmark.jl") end +@safetestset "Literal PMDSE — PGM EN generality" begin include("literal/test_pgm_en_generality.jl") end ambiguities = Test.detect_ambiguities(_PMDSE); if !isempty(ambiguities) From 90624a712a7f50096231caad661e825fe5cd24c0 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 12:33:37 +0000 Subject: [PATCH 2/2] Literal PMDSE: add verbose PowerGridModel tutorial (entry point) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A step-by-step, matrix-by-matrix walkthrough of the two PGM solve options as the entry point to the branch: - examples/literal_pmdse_pgm_tutorial.jl: runnable, heavily-commented script that builds everything from scratch and prints every intermediate object — per-unit data dict, Y_bus, the measurement atoms, the iterative-linear measurement matrix A and weights, the per-iteration |ΔU|, the Newton-Raphson gain HᵀWH, the comparison to PowerGridModel, then the same code on a four-wire explicit-neutral network (4×4 coupled admittance, neutral return = -Σ phase currents), and the reference schemes / robust estimators. - docs/src/literal_pmdse_tutorial.md: the narrated version with the real captured matrices; registered in docs/make.jl. Verified: the script and all new code run in Julia (1os2msr matches PGM to 6.5e-8 V for both methods; the four-wire EN case recovers the power-flow truth, neutral included, to 5.6e-16). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DnxB7KtuKVUYCK4fAQyahK --- docs/make.jl | 1 + docs/src/literal_pmdse_tutorial.md | 270 ++++++++++++++++++++++++ examples/literal_pmdse_pgm_tutorial.jl | 281 +++++++++++++++++++++++++ 3 files changed, 552 insertions(+) create mode 100644 docs/src/literal_pmdse_tutorial.md create mode 100644 examples/literal_pmdse_pgm_tutorial.jl diff --git a/docs/make.jl b/docs/make.jl index d3525fe..28467ae 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,6 +15,7 @@ makedocs( "Angular Reference Models" => "angular_ref.md", "Explicit Neutral Models" => "explicit_neutral_models.md", "Literal PMDSE (matrix-based)" => "literal_pmdse.md", + "Literal PMDSE: PowerGridModel Tutorial" => "literal_pmdse_tutorial.md", "Literal PMDSE: Reference & Observability" => "literal_pmdse_observability.md", "Bad Data Detection and Identification" => "bad_data.md", ], diff --git a/docs/src/literal_pmdse_tutorial.md b/docs/src/literal_pmdse_tutorial.md new file mode 100644 index 0000000..19ed4fb --- /dev/null +++ b/docs/src/literal_pmdse_tutorial.md @@ -0,0 +1,270 @@ +# Literal PMDSE tutorial — the PowerGridModel solve options, step by step + +This is the recommended **entry point** to the Literal PMDSE estimator. It walks, +matrix by matrix, through the two [PowerGridModel](https://github.com/PowerGridModel/power-grid-model) +(PGM) state-estimation *solve options* reproduced here — + +* `:iterative_linear` (PGM's default) and +* `:newton_raphson` + +— first on a small **single-phase PGM benchmark** (so we can check every number +against PowerGridModel itself), and then on a **four-wire, explicit-neutral** +network (which PGM cannot represent) using the *same* code. Every intermediate +object is printed: the per-unit data dictionary, the nodal admittance matrix +`Y_bus`, the measurement model, the iterative-linear measurement matrix, the +per-iteration voltage updates, and the Newton–Raphson gain. + +The runnable companion script is `examples/literal_pmdse_pgm_tutorial.jl`: + +```julia +julia --project examples/literal_pmdse_pgm_tutorial.jl +``` + +All matrices below are the *actual* output of that script. + +--- + +## Part A — a 3-bus single-phase PGM example + +We reproduce the `1os2msr` case from +`power-grid-model/tests/data/state_estimation/1os2msr`: three 10.5 kV nodes, two +lines, a source at node 1, and loads at nodes 2 and 3, with three voltage-phasor +sensors and seven power sensors. + +### A.1 Per-unit bases + +PGM solves in SI but normalises measurements internally; Literal PMDSE (like +PowerModelsDistribution) works in per-unit. We pick the **same bases PGM uses** +so the weighting — and therefore the estimate — matches exactly: + +``` +V_base = u_rated = 10500 V S_base = base_power_3p = 1e6 VA +Z_base = V_base² / S_base ≈ 110.25 Ω +``` + +A line's PGM parameters `(r1, x1, c1, tan1)` become a PMD branch with +`br_r = r1/Z_base`, `br_x = x1/Z_base`, and a Π-shunt split half to each end, +`b = ω·c1·Z_base/2`, `g = tan1·b`. + +### A.2 The nodal admittance matrix + +`LiteralModel` assembles `Y_bus` (here, one conductor per bus; node order +`[(1,1), (2,1), (3,1)]`). The nodal current injection is `I = Y_bus · U`: + +``` +Y_bus [p.u.] (3×3) + 457.72-498.53im -457.64+499.25im 0.0+0.0im + -457.64+499.25im 801.63-864.33im -343.81+366.74im + 0.0+0.0im -343.81+366.74im 343.91-365.8im +``` + +(The large magnitudes come from the small per-unit line impedances, +`r ≈ 1e-3` p.u.) + +### A.3 The measurement model — one "atom" per scalar measurement + +`build_se_atoms(lm, math)` turns `data["meas"]` (the same dictionary the JuMP +estimators consume) into a flat list of atoms. Each atom carries the real/imag +parts of two complex coefficient vectors over the nodal voltage `U`: + +* `ΔU(x) = cu·U` — a phase-to-neutral voltage (for voltage measurements and as the + `U` in `S = ΔU·conj(I)`), +* `I(x) = ci·U` — a branch current row, or a `Y_bus` row for a nodal injection. + +For `1os2msr` we get 13 atoms: + +``` +vim z=-0.0181 σ=0.01 (voltage angle/imag, node 1) +vre z=1.0238 σ=0.01 (voltage real, node 1) +vim z=-0.0134 σ=0.01 (node 2) +vim z=-0.0207 σ=0.0095 (node 3) +vre z=1.0234 σ=0.0095 (node 3) +vre z=1.0239 σ=0.01 (node 2) +power P=-2.0 Q=1.0 σ=1.0e9 (branch 4, to-side — disabled: σ huge) +power P=1.2304 Q=-1.7422 σ=0.0209 (branch 5, from) +power P=2.4124 Q=-3.024 σ=0.0379 (branch 4, from) +power P=-1.02 Q=-0.22 σ=0.0104 (branch 5, to) +power P=-1.01 Q=-0.21 σ=0.0103 (load injection, node 2) +power P=-1.02 Q=-0.22 σ=0.0104 (load injection, node 3) +power P=2.4124 Q=-3.024 σ=0.038 (source injection, node 1) +``` + +A voltage phasor sensor becomes two atoms (`:vre`, `:vim`); a power sensor becomes +one `:power` atom (`P`, `Q` with independent σ). Multiple sensors on the same +appliance are Kalman-combined; the disabled to-side sensor (σ = 1e9 p.u.) carries +essentially zero weight. + +### A.4 iterative_linear — the linear system it solves + +Each atom contributes one or two **real, linear** rows in the state +`x = [vr₁ vr₂ vr₃ vi₁ vi₂ vi₃]`. Stacking them gives the (constant) measurement +matrix `A` (20 rows here). The first eight rows: + +``` +A (first 8 of 20 rows) + 0.0 0.0 0.0 0.0 1.0 0.0 # vim node1 -> vi₁ + 1.0 0.0 0.0 -0.0 -0.0 -0.0 # vre node1 -> vr₁ + 0.0 0.0 0.0 1.0 0.0 0.0 # vim node2 + 0.0 0.0 0.0 0.0 0.0 1.0 # vim node3 + 0.0 0.0 1.0 -0.0 -0.0 -0.0 # vre node3 + 0.0 1.0 0.0 -0.0 -0.0 -0.0 # vre node2 + -457.64 457.72 0.0 -499.25 498.53 -0.0 # Re branch-4 current + 499.25 -498.53 0.0 -457.64 457.72 0.0 # Im branch-4 current +``` + +The voltage rows are trivial selectors; the power/current rows are rows of +`Y_bus`. `iterative_linear` keeps `A` and the weights `W = diag(1/σ²)` constant +and only re-builds the right-hand side `b` each iteration (the power +measurements are linearised to equivalent currents `I = conj(S/ΔU)` using the +previous voltages), solving `min ‖√W (A x − b)‖₂` each time. + +### A.5 Watch it converge + +`max_i |Uᵢ − Uᵢ_prev|` shrinks geometrically (ratio ≈ 0.23): + +``` +it 1 max|ΔU| = 3.3e-2 +it 2 = 6.7e-3 +it 3 = 1.5e-3 +it 5 = 7.8e-5 +it 10 = 4.8e-8 +it 16 = 6.7e-12 -> converged, weighted SSR = 0.0 +``` + +(The data is consistent — generated from a power flow — so the residual sum is +machine-zero; the disabled sensor is the only nonzero residual.) + +### A.6 newton_raphson — the Gauss–Newton gain + +`:newton_raphson` solves the same WLS by Gauss–Newton on the nonlinear model +`z = h(x)`, with `H = ∂h/∂x` (via `ForwardDiff`). It is warm-started from one +`iterative_linear` solve, so it converges in a single refining step here. The +normal-equation gain `G = HᵀWH` (6×6): + +``` +G = HᵀWH (entries ~1e10 — the per-unit impedances make G ill-conditioned; + the QR / orthogonal fallback in the solver handles it) + 5.19e9 -8.53e9 3.35e9 7.21e6 4.54e7 -3.76e7 + -8.53e9 1.98e10 -1.13e10 -4.53e7 -7.73e6 4.28e7 + 3.35e9 -1.13e10 7.93e9 3.76e7 -3.20e7 -1.03e7 + 7.21e6 -4.53e7 3.76e7 5.19e9 -8.54e9 3.35e9 + 4.54e7 -7.73e6 -3.20e7 -8.54e9 1.99e10 -1.13e10 + -3.76e7 4.28e7 -1.03e7 3.35e9 -1.13e10 7.97e9 +``` + +### A.7 Compare to PowerGridModel + +Both methods reproduce PowerGridModel's published `sym_output.json` voltages: + +| node | PGM \|U\| [V] | iterative_linear \|U\| | newton_raphson \|U\| | max \|U\| error | +|---|---|---|---|---| +| 1 | 10751.073 | 10751.073 | 10751.073 | 6.5e-8 V | +| 2 | 10752.699 | 10752.699 | 10752.699 | 6.5e-8 V | +| 3 | 10748.321 | 10748.321 | 10748.321 | 6.4e-8 V | + +The line flows (including the to-side back-calculation) and the disabled +sensor's residual also match PGM — this is gate 6 of the test suite +(`test/literal/test_pgm_benchmark.jl`). + +--- + +## Part B — the same solver, four-wire with an explicit neutral + +PGM is single-phase (Kron-reduced). Literal PMDSE keeps the neutral explicit: +every bus has terminals `[1,2,3,4]`, voltages are phase-to-neutral, and a wye +load's neutral current is `−Σ(phase currents)`. The *same* `build_se_atoms` / +`solve_se_*` code handles it — the single-phase case above is just the +one-conductor, no-neutral special case. + +### B.1 The 4×4 mutually-coupled branch + +``` +Y_branch = Z⁻¹ [p.u.] (row/col 4 = neutral) + 8.07-16.1im -1.93+3.9im -1.93+3.9im -0.71+0.98im + -1.93+3.9im 8.07-16.1im -1.93+3.9im -0.71+0.98im + -1.93+3.9im -1.93+3.9im 8.07-16.1im -0.71+0.98im + -0.71+0.98im -0.71+0.98im -0.71+0.98im 7.09-9.79im +``` + +### B.2 The neutral return current + +With per-phase current-injection measurements at the load, the estimator +synthesises the **neutral** nodal-injection atom as minus the sum of the phase +currents: + +``` +cinj at node (2,1) Ir=-0.298 Ii= 0.102 +cinj at node (2,2) Ir= 0.311 Ii= 0.292 +cinj at node (2,3) Ir= 0.106 Ii=-0.523 +cinj at node (2,4) Ir=-0.120 Ii= 0.129 # = −(−0.298+0.311+0.106) , −(0.102+0.292−0.523) +``` + +### B.3 Recovery, neutral included + +Both methods recover the power-flow truth to machine precision: + +``` +ITERATIVE_LINEAR: converged in 2 it, max|U − truth| = 5.55e-16 +NEWTON_RAPHSON: converged in 1 it, max|U − truth| = 5.66e-16 +bus-2 voltages |U| per terminal [1,2,3,N] = [0.99191, 0.98553, 0.98281, 0.01322] +``` + +The non-zero neutral voltage (`0.01322` p.u.) is exactly the unbalance effect a +four-wire study needs and that a Kron-reduced single-phase model cannot show. +This is gate 7 of the test suite (`test/literal/test_pgm_en_generality.jl`). + +--- + +## Part C — reference schemes and the robust estimators + +The reference scheme fixes the rotational gauge of the WLS problem: + +| `reference` | effect | +|---|---| +| `:prop` | fix only grounded terminals; the absolute angle is taken from a voltage-phasor measurement, or pinned automatically (`Im U_ref = 0`) when only magnitudes are measured. | +| `:sota` | additionally ground the reference-bus neutral and fix one angle datum. | +| `:full_slack` | fix the whole reference-bus phasor (e.g. to a known slack). | + +``` +reference schemes (method=:iterative_linear): :prop -> converged (13 it) + :sota -> converged (13 it) +``` + +Besides the two PGM solve options, the original JuMP-free **Gauss–Newton** +estimators remain available (leave `method` unset). They target the +explicit-neutral IVREN model: + +| `estimator` | meaning | +|---|---| +| `:wls` | weighted least squares (Gauss–Newton with a QR fallback) | +| `:wlav` | weighted least *absolute value* (IRLS) — robust to a single gross error | +| `:mle` | general maximum likelihood (Fisher scoring) — reduces to `:wls` for Gaussian noise | + +``` +Gauss-Newton estimators on the four-wire EN network: :wls / :wlav / :mle -> all converged +``` + +--- + +## API summary + +```julia +using PowerModelsDistributionStateEstimation + +# the two PowerGridModel solve options (general — single-phase or four-wire EN): +res = solve_mc_se_literal(data; estimator=:wls, method=:iterative_linear) +res = solve_mc_se_literal(data; estimator=:wls, method=:newton_raphson) + +# the original Gauss-Newton estimators (explicit-neutral IVREN model): +res = solve_mc_se_literal(data; estimator=:wls) # or :wlav, :mle + +# lower-level introspection: +lm = LiteralModel(data; reference=:prop) +atoms = build_se_atoms(lm, data) # the measurement model +res = solve_se_il(lm, atoms) # iterative_linear +res = solve_se_nr(lm, atoms) # newton_raphson +``` + +`res` is a [`LiteralResult`](@ref): `res.solution["bus"][id]["vr"|"vi"|"vm"]`, +`res.termination`, `res.iterations`, `res.objective`, `res.gain_cond`, +`res.gain_rank`. See [Literal PMDSE (matrix-based)](@ref) for the full reference. diff --git a/examples/literal_pmdse_pgm_tutorial.jl b/examples/literal_pmdse_pgm_tutorial.jl new file mode 100644 index 0000000..d157df0 --- /dev/null +++ b/examples/literal_pmdse_pgm_tutorial.jl @@ -0,0 +1,281 @@ +################################################################################ +# Literal PMDSE — PowerGridModel tutorial (verbose entry point) # +# # +# A heavily-commented, matrix-by-matrix walk through the two PowerGridModel # +# (PGM) state-estimation solve options reproduced in Literal PMDSE: # +# :iterative_linear and :newton_raphson # +# # +# We build everything from scratch so you can see every intermediate object: # +# the per-unit data dictionary, the nodal admittance matrix Y_bus, the # +# measurement model (one "atom" per scalar measurement), the linear # +# measurement matrix A and weights W used by iterative-linear, the per- # +# iteration voltage updates, the Newton-Raphson Jacobian H and gain HᵀWH, and # +# finally a comparison against PowerGridModel's own published results. # +# # +# Part A : a 3-bus single-phase PGM example (power-grid-model 1os2msr case) # +# Part B : the SAME code on a four-wire, explicit-neutral network # +# Part C : the other reference schemes and estimators # +# # +# Run (needs the package instantiated): # +# julia --project examples/literal_pmdse_pgm_tutorial.jl # +################################################################################ +import PowerModelsDistributionStateEstimation as _PMDSE +import PowerModelsDistribution as _PMD +import Distributions as _DST +import LinearAlgebra +const PMDSE = _PMDSE + +hr(t) = println("\n", "="^78, "\n ", t, "\n", "="^78) +cfmt(z) = string(round(real(z); digits=2), imag(z) ≥ 0 ? "+" : "", round(imag(z); digits=2), "im") +function showmat(name, M; cplx = eltype(M) <: Complex) + println("\n", name, " (", size(M, 1), "×", size(M, 2), ")") + for i in 1:size(M, 1) + print(" ") + for j in 1:size(M, 2) + print(lpad(cplx ? cfmt(M[i, j]) : string(round(M[i, j]; digits=4)), cplx ? 16 : 10), " ") + end + println() + end +end + +############################################################################ A +hr("PART A — 3-bus single-phase PGM example (the `1os2msr` validation case)") + +# --- A.1 the network, exactly as in power-grid-model/tests/data/state_estimation/1os2msr +# 3 nodes @ 10.5 kV, 2 lines, a source at node 1, loads at nodes 2 & 3. +# PGM works in SI; Literal PMDSE (like PMD) works in per-unit. We pick the +# SAME bases PGM uses internally: V_base = u_rated, S_base = base_power_3p = 1e6. +u_rated = 10.5e3 +Sbase = 1.0e6 +Zbase = u_rated^2 / Sbase # ≈ 110.25 Ω +ω = 2π * 50.0 +println("\nper-unit bases: V_base = $(u_rated) V S_base = $(Sbase) VA Z_base = $(round(Zbase;digits=3)) Ω") + +# PGM line parameters (Ω, F): r1, x1, c1, tan1 — π-model, shunt split half each end. +lines = [(id=4, f=1, t=2, r1=0.11, x1=0.12, c1=4.1380285203892784e-5, tan1=0.1076923076923077), + (id=5, f=2, t=3, r1=0.15, x1=0.16, c1=5.411268065124442e-5, tan1=0.10588235294117646)] + +# --- A.2 build the PMD *mathematical* data dictionary (per-unit, 1 conductor/bus) +function build_branch(l) + b1 = ω * l.c1; g1 = l.tan1 * b1 # shunt: b = ωC, g = tan·b + gsh = (g1 * Zbase) / 2; bsh = (b1 * Zbase) / 2 # half to each π end, in p.u. + Dict{String,Any}("index"=>l.id, "f_bus"=>l.f, "t_bus"=>l.t, + "f_connections"=>[1], "t_connections"=>[1], + "br_r"=>reshape([l.r1/Zbase],1,1), "br_x"=>reshape([l.x1/Zbase],1,1), + "g_fr"=>reshape([gsh],1,1), "b_fr"=>reshape([bsh],1,1), + "g_to"=>reshape([gsh],1,1), "b_to"=>reshape([bsh],1,1), "br_status"=>1) +end +bus(i, bt) = Dict{String,Any}("index"=>i, "bus_i"=>i, "bus_type"=>bt, "terminals"=>[1], + "grounded"=>Bool[false], "vbase"=>u_rated, "vmin"=>[0.0], "vmax"=>[Inf], + "vr_start"=>[1.0], "vi_start"=>[0.0]) +math = Dict{String,Any}( + "bus" => Dict("1"=>bus(1,3), "2"=>bus(2,1), "3"=>bus(3,1)), # node 1 = reference + "branch" => Dict(string(l.id)=>build_branch(l) for l in lines), + "gen" => Dict("6"=>Dict{String,Any}("index"=>6,"gen_bus"=>1,"connections"=>[1],"configuration"=>_PMD.WYE)), + "load" => Dict("7"=>Dict{String,Any}("index"=>7,"load_bus"=>2,"connections"=>[1],"configuration"=>_PMD.WYE,"pd"=>[0.0],"qd"=>[0.0],"model"=>_PMD.POWER,"status"=>1), + "8"=>Dict{String,Any}("index"=>8,"load_bus"=>3,"connections"=>[1],"configuration"=>_PMD.WYE,"pd"=>[0.0],"qd"=>[0.0],"model"=>_PMD.POWER,"status"=>1)), + "shunt" => Dict{String,Any}(), "settings"=>Dict{String,Any}("sbase"=>Sbase)) + +# --- A.3 the nodal admittance matrix Y_bus (the heart of the network model) +lm = PMDSE.LiteralModel(math; reference = :prop) # :prop = fix only grounded terminals +println("\nnode order (bus, terminal): ", lm.nodes) +Ybus = complex.(lm.G, lm.B) +showmat("Y_bus [p.u.] (I = Y_bus · U)", Ybus) + +# --- A.4 the measurements, in real units (volts / watts / vars), as PGM gives them +# voltage sensors carry a magnitude AND an angle (a phasor) -> rectangular vr/vi; +# power sensors sit on a branch terminal (from/to) or an appliance (source/load). +V(node,u,θ,σ) = (cmp=:bus, id=node, vars=[(:vr, u*cos(θ)), (:vi, u*sin(θ))], σ=σ) +vsens = [V(1, 10751.072595758282, -0.013054638926306409, 105.0), + V(2, 10752.698591183394, -0.017637349459726520, 105.0), + V(3, 10748.320749959701, -0.020182330759474040, 100.0)] +# (kind, object, side/appliance, P, Q, σ) +psens = [(:branch, 4, :from, 2412359.2976399013, -3024028.886598367, 37916.0), + (:branch, 4, :to, -2000000.0, 1000000.0, 1e15), # disabled (σ huge) + (:branch, 5, :from, 1230426.390004009, -1742195.1033582848, 20878.0), + (:branch, 5, :to, -1019999.9999999485, -219999.99999999927, 10435.0), + (:source, 6, :_, 2412359.297639887, -3024028.8865982923, 38009.0), + (:load, 7, :_, 1010000.0, 210000.0, 10316.0), + (:load, 8, :_, 1020000.0, 220000.0, 10435.0)] + +# write data["meas"] in the package vocabulary (per-unit), exactly as +# `add_measurements!` would produce it. +N(μ,σ) = _DST.Normal(μ, σ) +meas = Dict{String,Any}(); mid = Ref(0) +addmeas!(d) = (mid[] += 1; meas[string(mid[])] = d) +for s in vsens, (var, val) in s.vars + addmeas!(Dict{String,Any}("cmp"=>:bus, "cmp_id"=>s.id, "var"=>var, "dst"=>[N(val/u_rated, s.σ/u_rated)])) +end +for (kind, obj, side, P, Q, σ) in psens + if kind == :branch + addmeas!(Dict{String,Any}("cmp"=>:branch,"cmp_id"=>obj,"var"=>:p,"side"=>side,"dst"=>[N(P/Sbase, σ/Sbase)])) + addmeas!(Dict{String,Any}("cmp"=>:branch,"cmp_id"=>obj,"var"=>:q,"side"=>side,"dst"=>[N(Q/Sbase, σ/Sbase)])) + else + cmp, pv, qv = kind == :source ? (:gen, :pg, :qg) : (:load, :pd, :qd) + addmeas!(Dict{String,Any}("cmp"=>cmp,"cmp_id"=>obj,"var"=>pv,"dst"=>[N(P/Sbase, σ/Sbase)])) + addmeas!(Dict{String,Any}("cmp"=>cmp,"cmp_id"=>obj,"var"=>qv,"dst"=>[N(Q/Sbase, σ/Sbase)])) + end +end +math["meas"] = meas +println("\nbuilt ", length(meas), " scalar measurement entries (var=:vr/:vi/:p/:q/:pd/:qd/:pg/:qg)") + +# --- A.5 the measurement model: each measurement becomes an "atom" +# ΔU(x) = cu·U (a phase-to-neutral voltage) I(x) = ci·U (a branch/injection current) +atoms = PMDSE.build_se_atoms(lm, math; rescaler = 1.0) +println("\nmeasurement atoms (", length(atoms), "):") +for a in atoms + tag = a.kind in (:vre,:vim) ? "z=$(round(a.z1;digits=4)) σ=$(round(a.sigma;digits=4))" : + a.kind == :power ? "P=$(round(a.z1;digits=3)) Q=$(round(a.z2;digits=3)) σ=$(round(a.sigp;digits=4))" : + a.kind == :cinj ? "Ir=$(round(a.z1;digits=3)) Ii=$(round(a.z2;digits=3))" : "" + println(" ", rpad(string(a.kind),6), " ", tag) +end + +# --- A.6 iterative-linear: reconstruct the linear system A·x = b it solves +# (this is exactly what `solve_se_il` builds; we rebuild it here to see it). +function il_matrix(atoms, n) + rows = Vector{Float64}[] + for a in atoms + if a.kind == :vre ; push!(rows, vcat(a.cur, -a.cui)) + elseif a.kind == :vim ; push!(rows, vcat(a.cui, a.cur)) + elseif a.kind == :vmag ; push!(rows, vcat(a.cur, -a.cui)); push!(rows, vcat(a.cui, a.cur)) + else push!(rows, vcat(a.cir, -a.cii)); push!(rows, vcat(a.cii, a.cir)) + end + end + return permutedims(reduce(hcat, rows)) +end +A = il_matrix(atoms, lm.n) +showmat("iterative-linear measurement matrix A (rows = real/imag of every measurement; cols = [vr₁ vr₂ vr₃ vi₁ vi₂ vi₃])", A[1:min(8,end), :]) +println(" … (", size(A,1), " rows total)") + +hr("A.7 — run iterative_linear (watch max|ΔU| shrink each iteration)") +res_il = PMDSE.solve_mc_se_literal(math; estimator = :wls, method = :iterative_linear, + reference = :prop, maxiter = 50, tol = 1e-10, verbose = true) +println("\nconverged: ", res_il.termination, " in ", res_il.iterations, " iterations,", + " weighted SSR = ", round(res_il.objective; sigdigits = 4)) + +hr("A.8 — run newton_raphson and inspect the Gauss-Newton gain matrix") +res_nr = PMDSE.solve_mc_se_literal(math; estimator = :wls, method = :newton_raphson, + reference = :prop, maxiter = 50, tol = 1e-10) +# reconstruct H and the gain G = HᵀWH at the solution (the normal-equation matrix) +free, xfix, _ = PMDSE._gauge(lm, atoms); xf = res_nr.x_full[free] +H = _PMDSE.ForwardDiff.jacobian(xx -> PMDSE._nr_predict(atoms, lm.n, xfix, free, xx), xf) +z, w = PMDSE._nr_zw(atoms); G = transpose(H) * (w .* H) +println("Newton-Raphson: ", res_nr.termination, " in ", res_nr.iterations, " iteration(s)") +showmat("gain G = HᵀWH (note the ~1e10 entries: the per-unit line impedances make G ill-conditioned;\n the QR/orthogonal fallback in solve_wls/solve_se_nr handles it)", G) + +hr("A.9 — compare the estimate to PowerGridModel's published result") +# PGM node voltages for 1os2msr (tests/data/state_estimation/1os2msr/sym_output.json) +pgm = Dict(1 => (u=10751.072595758282, θ=-0.013054638926306409), + 2 => (u=10752.698591183394, θ=-0.017637349459726520), + 3 => (u=10748.320749959701, θ=-0.020182330759474040)) +println(rpad("node",6), rpad("PGM |U| [V]",16), rpad("IL |U| [V]",16), rpad("NR |U| [V]",16), "max |U| err [V]") +voltage(res, nid) = (b = res.solution["bus"][string(nid)]; (b["vr"][1] + im*b["vi"][1]) * u_rated) +for nid in 1:3 + Uil = voltage(res_il, nid); Unr = voltage(res_nr, nid); Ug = pgm[nid].u * exp(im*pgm[nid].θ) + println(rpad(nid,6), rpad(round(pgm[nid].u;digits=3),16), rpad(round(abs(Uil);digits=3),16), + rpad(round(abs(Unr);digits=3),16), round(max(abs(Uil-Ug), abs(Unr-Ug)); sigdigits=2)) +end + +############################################################################ B +hr("PART B — the SAME solver on a FOUR-WIRE, explicit-neutral network") +# PGM is single-phase (Kron-reduced). Literal PMDSE keeps the neutral explicit: +# every bus has terminals [1,2,3,4]; voltages are phase-to-neutral; a wye load's +# neutral current is −Σ(phase currents). Same `build_se_atoms`/`solve_se_*` code. + +va = deg2rad.([0.0, -120.0, 120.0]) +zp, zm, zn, zpn = 0.03+0.06im, 0.01+0.02im, 0.05+0.07im, 0.005+0.01im +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 +showmat("4×4 mutually-coupled series admittance Y_branch = Z⁻¹ [p.u.] (row/col 4 = neutral)", inv(Z)) + +busN(i,bt,gnd) = Dict{String,Any}("index"=>i,"bus_i"=>i,"bus_type"=>bt,"terminals"=>[1,2,3,4], + "grounded"=>Bool[false,false,false,gnd],"vbase"=>1.0,"vmin"=>zeros(4),"vmax"=>fill(Inf,4), + "vr_start"=>[cos.(va)...,0.0],"vi_start"=>[sin.(va)...,0.0]) +mathEN = Dict{String,Any}( + "bus"=>Dict("1"=>busN(1,3,true), "2"=>busN(2,1,false)), + "branch"=>Dict("1"=>Dict{String,Any}("index"=>1,"f_bus"=>1,"t_bus"=>2,"f_connections"=>[1,2,3,4], + "t_connections"=>[1,2,3,4],"br_r"=>real.(Z),"br_x"=>imag.(Z),"br_status"=>1)), + "gen"=>Dict("1"=>Dict{String,Any}("index"=>1,"gen_bus"=>1,"connections"=>[1,2,3,4],"configuration"=>_PMD.WYE)), + "load"=>Dict("1"=>Dict{String,Any}("index"=>1,"load_bus"=>2,"connections"=>[1,2,3,4], + "configuration"=>_PMD.WYE,"pd"=>[0.3,0.4,0.5],"qd"=>[0.1,0.12,0.15],"model"=>_PMD.POWER,"status"=>1)), + "shunt"=>Dict{String,Any}(), "settings"=>Dict{String,Any}("sbase"=>1.0)) + +# Ground truth by a quick Newton power flow (so we can check the estimator recovers it) +lmEN = PMDSE.LiteralModel(mathEN; reference = :prop); ni = lmEN.node_index; n = lmEN.n +Yen = complex.(lmEN.G, lmEN.B); Sload = [0.3+0.1im, 0.4+0.12im, 0.5+0.15im] +Uslk = ComplexF64[exp(im*va[1]), exp(im*va[2]), exp(im*va[3]), 0.0] +U = ComplexF64[(b == 1 ? Uslk[t] : (t < 4 ? exp(im*va[t]) : 0.0)) for (b,t) in lmEN.nodes] +fixed = [ni[(1,t)] for t in 1:4]; freeN = [i for i in 1:n if !(i in fixed)] +resid(x) = (Uu = copy(U); for (j,i) in enumerate(freeN); Uu[i] = x[j] + im*x[length(freeN)+j]; end; + I = Yen*Uu; Ic = zeros(ComplexF64, n); + for ph in 1:3; kc = ni[(2,ph)]; kn = ni[(2,4)]; Il = conj(Sload[ph]/(Uu[kc]-Uu[kn])); Ic[kc]-=Il; Ic[kn]+=Il; end; + vcat(real.((Yen*Uu .- Ic)[freeN]), imag.((Yen*Uu .- Ic)[freeN]))) +let x = vcat(real.(U[freeN]), imag.(U[freeN])) + for _ in 1:60 + r = resid(x); J = zeros(length(x), length(x)); h = 1e-7 + for j in eachindex(x); xp = copy(x); xp[j] += h; J[:,j] = (resid(xp) .- r) ./ h; end + d = J \ r; x .-= d; LinearAlgebra.norm(d, Inf) < 1e-13 && break + end + for (j,i) in enumerate(freeN); U[i] = x[j] + im*x[length(freeN)+j]; end +end +Utrue = copy(U) + +# measurements: voltage phasor at the slack + per-phase current injections at the load +mEN = Dict{String,Any}(); kc = Ref(0); add!(d) = (kc[] += 1; mEN[string(kc[])] = d) +add!(Dict{String,Any}("cmp"=>:bus,"cmp_id"=>1,"var"=>:vr,"dst"=>[N(real(Utrue[ni[(1,p)]]),1e-4) for p in 1:3])) +add!(Dict{String,Any}("cmp"=>:bus,"cmp_id"=>1,"var"=>:vi,"dst"=>[N(imag(Utrue[ni[(1,p)]]),1e-4) for p in 1:3])) +crd = [conj(Sload[ph] / (Utrue[ni[(2,ph)]] - Utrue[ni[(2,4)]])) for ph in 1:3] # load currents +add!(Dict{String,Any}("cmp"=>:load,"cmp_id"=>1,"var"=>:crd,"dst"=>[N(real(crd[p]),1e-4) for p in 1:3])) +add!(Dict{String,Any}("cmp"=>:load,"cmp_id"=>1,"var"=>:cid,"dst"=>[N(imag(crd[p]),1e-4) for p in 1:3])) +mathEN["meas"] = mEN + +ref = (Dict(t=>real(Utrue[ni[(1,t)]]) for t in 1:4), Dict(t=>imag(Utrue[ni[(1,t)]]) for t in 1:4)) +atomsEN = PMDSE.build_se_atoms(PMDSE.LiteralModel(mathEN; reference=:full_slack, ref_values=ref), mathEN) +println("\ncurrent-injection atoms at the load bus (terminal 4 is the neutral return = −Σ phases):") +for a in atomsEN + a.kind == :cinj || continue + nd = haskey(a.meta, :node) ? PMDSE.LiteralModel(mathEN; reference=:full_slack, ref_values=ref).nodes[a.meta.node] : "" + println(" cinj at node ", nd, " Ir=", round(a.z1; digits=3), " Ii=", round(a.z2; digits=3)) +end + +for method in (:iterative_linear, :newton_raphson) + r = PMDSE.solve_mc_se_literal(mathEN; estimator=:wls, method=method, + reference=:full_slack, ref_values=ref, maxiter=100, tol=1e-12) + err = maximum(abs((r.solution["bus"][string(b)]["vr"][findfirst(==(t), mathEN["bus"][string(b)]["terminals"])] + + im*r.solution["bus"][string(b)]["vi"][findfirst(==(t), mathEN["bus"][string(b)]["terminals"])]) - + Utrue[ni[(b,t)]]) for (b,t) in lmEN.nodes) + println("\n", uppercase(string(method)), ": ", r.termination, " in ", r.iterations, + " it, max|U − truth| = ", round(err; sigdigits=3)) + b2 = r.solution["bus"]["2"] + println(" bus-2 voltages |U| per terminal [1,2,3,N] = ", round.(b2["vm"]; digits=5), + " (neutral is non-zero: the 4-wire effect PGM cannot represent)") +end + +############################################################################ C +hr("PART C — reference/gauge schemes and the robust estimators") +# The PGM solve options (`method=:iterative_linear|:newton_raphson`) are the +# general path and run on either network. The reference scheme fixes the gauge: +# :prop — fix only grounded terminals; the rotation is taken from a +# voltage-phasor measurement, or pinned automatically (Im U_ref=0) +# when only magnitudes are measured. +# :sota — additionally ground the reference-bus neutral + one angle datum. +# :full_slack — fix the whole reference-bus phasor. +println("reference / gauge schemes on the single-phase network (method=:iterative_linear):") +for r in (:prop, :sota) + res = PMDSE.solve_mc_se_literal(math; estimator=:wls, method=:iterative_linear, reference=r) + println(" ", rpad(string(r),11), " term=", res.termination, " iters=", res.iterations) +end + +# The Gauss-Newton estimators :wls/:wlav/:mle (the original JuMP-free path, +# selected by leaving `method` unset) target the explicit-neutral IVREN model, +# so we exercise them on the four-wire network from Part B. :wlav (least +# absolute value) is robust to a single gross error; :mle reduces to :wls for +# Gaussian noise. +println("\nGauss-Newton estimators on the four-wire EN network (method unset): :wls / :wlav / :mle") +for est in (:wls, :wlav, :mle) + res = PMDSE.solve_mc_se_literal(mathEN; estimator=est, reference=:full_slack, ref_values=ref) + println(" ", rpad(string(est),5), " term=", res.termination, " iters=", res.iterations, + " SSR=", round(res.objective; sigdigits=3)) +end +println("\nDONE. See docs/src/literal_pmdse_tutorial.md for the narrated version.")