State estimation

The weighted-least-squares (WLS) state estimator fits a network state to redundant, noisy telemetry: given measurements z = h(x) + e, it finds the state x^ minimising

J(x) = sum_i (z_i - h_i(x))^2 / sigma_i^2

Unlike a power flow - which needs a complete, exactly-determined P/Q specification at every bus - the estimator takes whatever measurements exist, wherever they are, and lets redundancy fight noise. Everything documented here lives in src/GridSim.Core/Estimation/; the CLI surface is the estimate verb (user manual chapter 10).

The state vector

The state is x = [Va at every non-slack bus (radians); Vm at every bus (pu)], in model bus order (Estimation/MeasurementModel.cs:MeasurementModel). Slack-bus angles are fixed throughout as the reference, so StateSize = (n - slackCount) + n. Angle columns come first (AngleColumn(busIndex), -1 for a slack), then Vm columns (VmColumn(busIndex) = AngleCount + busIndex). These are the same conventions as the Newton power flow, so estimates are directly comparable to solver results.

The measurement model: h(x) per kind

A Measurement (Estimation/Measurement.cs) is (Kind, Value, Sigma, Bus, From, To, End, Circuit, Source). Bus measurements address the bus by id (not solver index); branch measurements address the branch by its (From, To) endpoints plus Circuit - the branch name, or the ordinal form "#k" (model branch index) - to disambiguate parallel circuits. Values travel in engineering units and are converted to internal per-unit once, at MeasurementModel construction.

Kind	h(x)	Unit conversion at ingestion
`VmPu`	`Vm[bus]`	none (pu)
`VaDeg`	`Va[bus]`	value and sigma x pi/180 (the state's angle unit)
`PInjMw`	net real injection at the bus	/ `BaseMva`
`QInjMvar`	net reactive injection	/ `BaseMva`
`PFlowMw`	real flow into the measured end	/ `BaseMva`
`QFlowMvar`	reactive flow into the measured end	/ `BaseMva`

VaDeg is a PMU/phase-angle row: h(x) = Va[bus] is linear in the state, so the row costs the Jacobian a single constant entry (1.0 in its angle column). A VaDeg row on a slack bus has no state column at all - the reference is fixed - so it is structurally empty: numerically a constant residual, legal but uninformative. A real feed (e.g. Gridradar) carries angles against its own arbitrary reference; re-reference against the slack before authoring rows - the estimator does not solve for a per-feed offset.

Injection rows are evaluated over the bus's electrical neighbourhood only (MeasurementModel.ComputeInjectionsSparse) - numerically identical to NewtonRaphsonPowerFlow.ComputeInjections (the skipped terms are exact structural zeros, and adding 0.0 is exact in IEEE 754; an agreement test pins the bit-for-bit equality), at O(nnz) instead of O(n^2). Flow rows use the same pi-model admittance stamps (yff, yft, ytf, ytt, including tap ratio and phase shift) as YBus.Build and NewtonShell.ComputeLineFlows; End selects which terminal is metered. The Jacobian partials are the standard polar power-flow forms, laid out over the full estimator state, and there is exactly one evaluation path - MeasurementModel.EvaluateJacobian<TSink> writes into either a dense array or per-row sparse lists via Estimation/JacobianSinks.cs, so the dense and sparse backends cannot drift apart.

MeasurementSet (Estimation/MeasurementSet.cs) is the versioned input document (schema gridsim-measurements/1, camelCase JSON, enums by name). Load rejects unknown schemas and horizon-checks the optional TakenUtc stamp (see chapter 8); Validate(model) throws one exception naming every offender - non-finite values, non-positive sigmas, unknown bus ids, unresolvable or ambiguous branch references, and a base-MVA mismatch - so a NaN or zero-weight row can never silently poison the gain matrix.

The phase-angle feed

Estimation/PhaseAngleFeed.cs turns per-bus voltage phase-angle observations into VaDeg rows. It exists because bus angles are the weakest-observed directions in the real GB runs: the P-rich, VAr-blind measurement layer pins them only indirectly through the P-injection/P-flow rows, so they surface as the near-singular gain pivots that drive the trust cap and the backfill machinery. A genuine per-bus angle adds a direct 1.0 pivot on that bus's angle column, turning an indirectly-observed direction into a directly-observed one, anchoring the angle profile and adding redundancy that sharpens the angle covariance diag(G^-1) (StateVaSigmaDeg).

The one physics that governs correct use: a per-bus absolute angle is phi_i(t) = theta_sys(t) + delta_i, where theta_sys is the system common-mode phase (the COI angle vs the 50 Hz frame, driven by frequency) and delta_i is the bus's spatial spread (power-flow determined). The estimator solves delta_i against a fixed slack, so re-referencing any absolute feed to slack - z_i = phi_i - phi_datum + prior = delta_i - delta_datum - cancels theta_sys exactly. Two consequences the code enforces:

Genuine per-bus differences improve the estimate (proven: injecting the true per-bus angles cuts the estimated-angle RMSE by more than half - PhaseAngleFeedTests.Genuine_per_bus_angles_sharpen_the_estimate).
A single system angle cannot - injected identically on every bus it cancels against the slack. So GDA's frequency-integrated PHASEANGLE_INFERED (a national COI series, the same for every bus; GridSim.Gda/Neso/PhaseAngleInferedReader.cs) is used only as the absolute reference frame, the true instantaneous frequency/RoCoF, and the common-mode to de-rotate an absolute per-bus feed - never authored as a per-bus row.

ToMeasurements conditions each observation in three optional steps, by mode (AutoReference (default) / AlreadyReferenced / Raw): time-project to a common instant using the observation's frequency and RoCoF (theta(T) = theta(t) + 360*(f-50)*dt + 180*rocof*dt^2), subtract an optional system common-mode (theta_sys(T) from the GDA reader), and re-reference to the slack via differences at a datum bus (the datum's model VaDeg is its prior anchor; the slack is the default datum). Slack-bus rows are dropped as uninformative. The CLI exposes it as estimate --phaseangle BUS:ANGLE:UTC:FREQ ... (--phaseangle-ref, --phaseangle-sigma, --phaseangle-at, --phaseangle-from-gda); injected rows carry the same TakenUtc horizon guard as every dated input. GDA authors no genuine per-bus angle yet - the override is the injection path, and build_estimator_measurements.py carries a documented (currently empty) _phase_angle_rows hook for a future PMU/synchrophasor feed.

The estimation loop

Estimation/WlsStateEstimator.cs:WlsStateEstimator.Estimate owns the loop - start point, trust-capped damping, convergence, chi-square consistency, bad-data removal - while everything representation-dependent lives behind the IWlsStep seam.

Start point. Flat start (Vm = 1, Va = the slack reference angle) by default; EstimationOptions.WarmStart seeds from a previously solved or estimated PowerFlowResult, matched by bus id (the natural chaining for per-timestamp replay estimation, and needed for real data). Slack buses always keep their own fixed angle. The start state is also cloned as the backfill prior (below).

Gauss-Newton iteration. Each pass binds a step backend (step.Begin(model, mm, w) with weights w = 1/sigma^2), gates observability at the start point, then iterates: step.TrySolve computes the undamped correction deltax; ApplyStep applies it with trust-region capping - if any angle component exceeds MaxAngleStepRad (0.2 rad) or any Vm component exceeds MaxVmStepPu (0.1 pu), the whole step is scaled down uniformly, so the direction is preserved. Caps never bind near the optimum, so converged answers (and the machine-precision recovery gates) are unchanged; they exist to stop a stressed real-data start overshooting voltage magnitudes through zero. Convergence is ||delta-x||inf < Tolerance (1e-8); the other terminal statuses are MaxIterations (30 default), NonFinite, and Unobservable. EstimationOptions.Trace (an EstimationTrace) optionally records each iteration's ||delta-x||inf and the trust scale actually applied, per pass - carried as an additive block in gridsim-estimate-export/1 and re-checked by tools/crossval/se_residual_check.py.

A singular gain mid-iteration (after the gate passed at the start point) means the state wandered somewhere the measurements no longer pin down; it is reported as Unobservable, never as NaN.

The `IWlsStep` seam and its four backends

Estimation/WlsStep.cs:IWlsStep defines what one Gauss-Newton pass needs from a representation: Begin (bind + structural work), GateObservability, TrySolve (deltax, false on a singular gain), and NormalizedResiduals (residual-covariance diagonals + criticality).

Dense LU normal equations (`DenseWlsStep`)

The at-or-below-threshold backend. Forms G = H^TWH and g = H^TW(z-h) densely (W is diagonal, so both are single passes over H; only the upper triangle is accumulated and then mirrored) and solves through the injected ILinearSolver (dense LU with partial pivoting). Its observability gate is the full analysis: island sweep plus gain-matrix rank proof.

Sparse fixed-pattern CSC with RCM-once (`SparseWlsStep`)

The gb-full-scale backend, selected when n > SparseBusThreshold (SolverSelector.DefaultSparseThreshold = 3,000 buses). H is held as per-measurement sparse rows. The key property: the gain's sparsity pattern is structural - it depends on measurement topology, not the operating point - so Begin does all structural work exactly once per bad-data pass:

MeasurementModel.StructuralColumns lists every state column a row can ever write (injection rows: the bus's electrical neighbourhood; flow rows: four columns; voltage rows: one);
the gain pattern G = H^TWH is built once as the union of each row's column-pair outer product; because G is symmetric, the CSR arrays double as the CSC arrays;
the fill-reducing ordering (SparseLu.Ordering, reverse Cuthill-McKee) is computed once;
a per-row slot map records the value-array index of each column pair (binary search once, not per iteration).

Each iteration then only refills values and refactors numerically (SparseLu.FactorPrepared, Gilbert-Peierls). This removes the per-iteration CSR rebuild + re-sort + re-ordering the first implementation paid.

Its observability gate is the island sweep only - the O(s^3) rank proof does not scale - so ObservabilityReport.Rank is reported as -1 (not computed) and numerical singularity surfaces through the factorisation returning null.

Constrained KKT (`ConstrainedDenseWlsStep`, `--method constrained`)

Rows tagged Source == "virtual-zero-injection" are promoted from soft pseudo-measurements to hard equalities via the Hachtel-style augmented (KKT) system

[ G_f  C^T ] [deltax]   [ g_f ]
[ C    0  ] [lambda ] = [ r_c ]

where G_f/g_f are the free rows' normal equations and C the constraint Jacobian. A pure switching node's zero net injection is a network fact, not telemetry - the constrained method honours it exactly instead of "within 10 MVAr", and conditions better than the huge-weight trick. The indefinite KKT matrix goes through the same partial-pivot LU. The observability gate uses the combined information (free rows at their weights, constraints at a strong gate-only weight 1e6): rank deficiency there is rank deficiency of [H_f; C]. For identification, the free rows use the plain covariance with constraints as very-high-weight rows (a documented approximation); constraint rows are marked critical outright with NaN normalized residuals - a network fact is never removable telemetry.

The backend is dense only: above the sparse threshold the estimator throws NotSupportedException - it fails loud rather than silently running unconstrained.

The IRLS outer loop (`--method huber` / `--method lav`)

The robust methods reuse the dense/sparse steps but wrap the Gauss-Newton inner loop in an iteratively-reweighted outer loop (up to 10 outer passes). After each converged inner solve, effective weights are re-derived from the standardized residuals r~_i = |z_i - h_i|/sigma_i at the solution:

Method	Weight factor on w_i
Huber (k = 1.345)	1 for r~ <= k; k/r~ beyond
LAV	1/max(r~, 1e-3)

The loop stops when the largest relative weight-factor change drops below 1e-3. Outliers are down-weighted smoothly, never removed - robust to several simultaneous gross errors where the sequential removal loop is slow or fragile. Two conventions follow: the chi-square statistic is computed against the original sigmas (not the IRLS weights) so it stays comparable across methods, and the removal loop is skipped entirely - the down-weighting already did the job.

Observability

Estimation/ObservabilityAnalysis.cs implements two checks:

Island sweep (AnalyzeTopology, both backends): connected components of the branch graph that no measurement touches at all are unobservable by construction; their bus ids are listed in UnmeasuredIslandBuses.
Gain rank (Analyze, dense backends): pivoted row-echelon elimination of G = H^TWH with a relative pivot floor (norm(G) x 1e-12, the same doctrine as the LU solvers). Columns that never yield a usable pivot are exactly the free state directions; they are named ("Va@bus 12"/"Vm@bus 12") in UnobservableStates with the affected bus ids in UnobservableBuses.

The report's Observable flag is the gate: when false, the run terminates with EstimationStatus.Unobservable and the reported state is the start point, not an estimate. On the sparse backend Rank = -1 by convention (not computed); singularity is caught by the factorisation instead.

The chi-square consistency gate

With redundancy (DegreesOfFreedom = m - StateSize > 0), the converged objective J(x^) is tested against the chi-square quantile at ChiSquareConfidence (0.95): ChiSquarePassed = J <= chi-square(dof, p). Estimation/ChiSquare.cs computes the quantile dependency-free with the Wilson-Hilferty cube approximation

chi-square(k, p) ~ k * (1 - 2/(9k) + z_p*sqrt(2/(9k)))^3

using Acklam's rational approximation for the standard-normal quantile z_p (relative error < 1.15e-9). Accuracy is well under 1 % for k >= 3 - ample for a detection gate. With dof <= 0 the threshold is NaN and the test is vacuously true: an undetectable gross error is still possible then, and the result says so through DegreesOfFreedom.

When the test passes and EstimationOptions.LeanResiduals is set, the m residual-covariance solves are skipped (the removal loop would never run; normalized residuals are then NaN in the residual table). The default keeps the fully-populated diagnostics.

Bad data: largest normalized residual with critical protection

When chi-square fails, Estimation/BadDataAnalysis.cs computes, at the converged state, the residual-covariance diagonal

Omega_i_i = sigma_i^2 - h_i^T * G^-1 * h_i         r^N_i = |r_i| / sqrtOmega_i_i

(one gain factorisation, then one back-substitution per measurement; the sparse path reuses one Gilbert-Peierls factorisation for all m solves). A measurement with Omega_i_i <= sigma_i^2*1e-8 is critical: it is the only observation of some state direction, its residual is structurally zero, and removing it would make the system unobservable - it is flagged (IsCritical, r^N = NaN) and never removed.

The removal loop then identifies the largest r^N among removable rows - not critical, not a prior-backfill pseudo - above NormalizedResidualThreshold (3.0, the classical 3sigma rule), records it in Removed with its pass number, excludes it, and re-estimates warm from the current state. The loop is capped at MaxBadDataRemovals (20). If chi-square fails but nothing is identifiable/removable, the failure is reported honestly (ChiSquarePassed = false) rather than force-fitted.

Prior backfill (`BackfillRemovedFromPrior`, off by default)

Observability restoration on removal: when the loop removes a real-power measurement (PInjMw/PFlowMw), a pseudo-measurement is substituted, evaluated from the start state (the warm start / measured-injection solved prior - for --prior=solve pipelines the best model-consistent guess available): value = h_prior x BaseMva, with a deliberately looser sigma max(3*sigma_removed, 5%*|value|), tagged Source = "removal-backfill-prior" (WlsStateEstimator.BackfillSource) and itself never removable. Without this, removing the only informative constraint in a weakly-redundant region (a radial GSP's sole injection) leaves the state free to run away from every remaining measurement; with it, the fit stays anchored to the physically-consistent prior.

Reactive kinds are deliberately not backfilled. In the GB measurement layer no VAr telemetry exists: the Q rows are themselves pseudos, and the solved prior embeds the measured Q - so backfilling a systematically defective Q row re-injects its own defect and churns the removal loop to its cap on otherwise-healthy periods. The option is opt-in because it changes what "removed" means (excluded -> replaced-by-prior).

The measured-injection prior

Estimation/MeasuredInjectionPrior.cs turns a measurement set into a measured-injection power flow solved by the validated Newton engine, whose solution serves as both the Gauss-Newton warm start and the value source for voltage priors. Per-bus net injections are read off the set's PInjMw/ QInjMvar rows (duplicates summed, unmeasured buses zero). Buses injecting at least PvThresholdMw (50 MW) become PV at PvSetpointPu (1.01) with a +/-0.6*capacity var range - a generating site holds its voltage; a pure-PQ conversion sags unrealistically. Everything else becomes PQ with the measured injection as negative load; the slack keeps its reference role and absorbs the residual imbalance. ApplyVoltagePrior then replaces "pseudo-voltage-schedule" rows with the solved voltage profile at sigma 0.02 pu, retagged "prior-solved" - the prior is now physics, not a flat number.

The synthetic sampler

Estimation/SyntheticMeasurements.cs:SyntheticMeasurements.Sample draws z = h(x_true) + N(0, sigma)*NoiseScale from a solved power-flow truth - the estimator's test fixture and demo feed (estimate ... --synth). It is deterministic per options: one seeded RNG (default seed 20260713), buses and branches walked strictly in model order, Box-Muller Gaussians, and NoiseScale = 0 draws nothing from the RNG (the machine-precision recovery fixture - exact values under the declared sigmas). Defaults: sigma(Vm) 0.004 pu, sigma(P/Q inj) 1.0 MW, sigma(P/Q flow) 1.6 MW; coverage fractions select meters by a deterministic stride over model order; flow rows are emitted for both branch ends. Parallel branches are keyed by the ordinal circuit form "#k" - real cases (gb-full SGT banks) carry duplicate names on parallel branches, so a name is only a safe key when unique. The written set carries a fixed, safely-historical TakenUtc (2024-03-01T12:00Z): the sampler never reads the clock.

Results and trace

Estimation/EstimationResult.cs:EstimationResult carries: Status/ Converged (true only for Converged - every other status means the estimate is untrustworthy), Iterations (across all passes), MaxStateUpdate (final ||delta-x||inf), BusIds/Vm/VaDeg (degrees, slack-referenced), ObjectiveJ, DegreesOfFreedom, ChiSquareThreshold, ChiSquarePassed, Residuals (per measurement: predicted, residual in engineering units, normalized residual or NaN, criticality), Removed (in removal order, with r^N and pass), and Observability.

ToPowerFlowResult(model) projects the estimate into the solved-state shape the rest of GridSim consumes - injections, losses and line flows recomputed at x^ through the same assembly path as the power-flow solvers - so replay, UI and export plumbing work on estimated states unchanged (this is how estimated GB states stream over the bridge as source=estimated).

The GDA-side sigma model

Real GB measurement sets are authored in the GDA data lake by two builders and consumed here via MeasurementSet.Load. The sigmas are a documented, deliberately simple model; the source tag on every row names its composition (sigma parts combine root-sum-square).

Stage 1 - D:\Work\GDA\v1\NetworkModel\build_gsp_measurements.py aggregates metered per-unit generation onto NESO GSPs through the BMU crosswalk:

Quantity	Sigma	Notes
B1610 per-GSP generation	`max(2 MW, 2%*\|MW\|) / min_confidence`	settlement-grade metering; inflated by the weakest crosswalk confidence in the aggregate (`provenance="observed"`)
PN fallback generation	`max(10 MW, 10%*\|MW\|) / min_confidence`	Physical Notifications are submitted intent, not metering (`provenance="pn"`); used when the day has no (or a partial, < 40-period) B1610 file; extends estimable history to 2016-03
Unallocated (gsp=NULL) bucket	30 %	a known-missing bucket, not a measurement - never smeared across the network
Per-PES-group net demand	`max(10 MW, 3%*\|MW\|)`	from the supplier settlement legs of the same B1610 read; not emitted on PN days

Stage 2 - D:\Work\GDA\v1\NetworkModel\build_estimator_measurements.py composes the gridsim-measurements/1 set for one settlement period on gb-spine/gb-full. Per-bus PInjMw = metered generation + interconnector flow

scaled folded LTDS demand; the sigma is the RSS of the parts:

Component	Sigma	Source tag fragment
Metered generation	carried from stage 1 (summed in quadrature)	`b1610` / `pn`
Interconnector flow	`max(5 MW, 2%*\|MW\|)`	`ic`
Demand pseudo, national-scaled	`max(5 MW, 25%*\|P_load\|)`	`pseudo-demand`
Demand pseudo, group-anchored	`max(5 MW, 15%*\|P_load\|)`	`pseudo-demand-group`
`QInjMvar` demand pseudo	`max(10 MVAr, 30%*\|Q_load\|)`	`pseudo-demand`
Switched-compensation widening	`+ (range/sqrt12)^2` in quadrature	`...+comp-range`
Zero-centred reactive pseudo (metered bus, no demand row)	`max(50, 0.3*\|P\|)`	`pseudo-reactive`
Virtual zero injection (pure switching nodes)	10 MW (P), 20 MVAr (Q)	`virtual-zero-injection`
Slack voltage reference	1.02 pu +/- 0.005	`reference`
Voltage-schedule prior, every other bus	1.01 pu +/- 0.03	`pseudo-voltage-schedule`

The two demand-pseudo fractions encode where the information actually is: a group-anchored pseudo (15 %) scales the bus's folded LTDS demand to its PES group's measured net offtake - the group total is telemetry (3 %), only the within-group shape is prior - while a national-scaled pseudo (25 %) rests on one national demand scalar and a uniform shape. The comp-range widening declares real uncertainty: a switched MSC/reactor/SVC's state is unobservable from settlement data, so a compensated site's Q pseudo carries the ETYS device range as extra sigma (uniform-distribution standard deviation = range/sqrt12). The reactive layer is entirely pseudo - no VAr telemetry exists in the lake - which is precisely why the estimator's prior-backfill option excludes Q kinds, and why the virtual-zero-injection rows are what the constrained method promotes to hard equalities.

Both builders are horizon-guarded on the GDA side, and the set is horizon-checked again by MeasurementSet.Load - the double-guarding described in chapter 8.