Technical manual — chapter list

Physics conformance and the constrained estimator

This chapter documents the physics-law audit engine, the ten laws it scores and which of them gate a state, the physics-constrained state estimator built on that gate, and the two instruments (a C# verb and an independent NumPy re-derivation) that keep the audit honest.

The motivating problem: the WLS state estimator could produce physically impossible states. On real GB data the undamped Gauss-Newton solve "walked Vm through zero", and a VAr-blind measurement set drags Vm to collapse. Those were patched reactively (trust-region damping, a voltage-schedule prior, the zero-injection KKT). This subsystem generalises those patches into one idea: the physical laws define the admissible region, and the estimator is constrained to live inside it.

The audit engine (Core/Physics/)

PhysicsLaws.Audit(model, result, tol, synchronism?) re-derives the network laws from the raw branch data by an independent per-branch pi-model path - it accumulates each branch's S_from/S_to and the bus shunt draw, and never calls IPowerFlowSolver. So a pass is genuine corroboration that the solved voltages satisfy the network equations, not a check of the solver against itself. The per-branch loop lives once in PhysicsAuditContext; the laws read its shared arrays.

Audit returns a PhysicsLawReport (Physics/PhysicsLawReport.cs): a LawResult per law (number, name, statement, LawStatus, residual, unit, tolerance, detail) plus Gate(hardLaws) - the acceptance verdict.

The ten laws, and what is checkable

GridSim is a positive-sequence phasor model plus an electromechanical frequency layer. That makes the circuit laws and AC synchronism directly checkable; it is not an electromagnetic-field model, so the field laws are embodied constitutively rather than verified. The audit states this honestly rather than fabricating a residual.

# Law Status class What is checked
1 Conservation of energy hard independent total dissipation == reported loss
2 Ohm's law hard branch flow y·dV matches the reported line flow
3 Joule heating hard branch-loss sum consistent with the reported total
4 Faraday induction constitutive transformer tap / gen EMF; no field to check
5 Ampere's law (passivity) reported passive branches dissipate: I²R >= 0
6 Kirchhoff current hard injection balance at constrained buses
7 Kirchhoff voltage structural nodal potentials are single-valued (loops vanish)
8 Complex power S=VI* hard the reactive component of the same balance
9 AC synchronism dynamics rotor-angle spread / frequency excursion margin
10 Gauss / Coulomb not modelled field/geometry, outside a phasor model

PhysicsCheckRegistry.HardGateLaws = {1, 2, 3, 6, 8} - the laws whose violation is physically impossible. A breach of any of these means the state is not a real network state, so the constrained estimator rejects it. Laws 1, 6 and 8 are three views of the single power-balance residual and are reported as such - no fabricated independence.

The Joule-vs-Ampere passivity split (a real finding)

The corpus sweep flagged case145, case588 and case9241pegase on Joule heating with large residuals. They converge and conserve energy exactly, but carry a branch with negative real loss - their data has a negative-resistance branch, a network-reduction/equivalencing artefact. A negative I²R is a model passivity property, not a bad state, and the estimator cannot fix the model. So Joule (hard) checks only loss consistency, and passivity moved to Ampere (law 5), reported when violated but not hard-gated. The gate then correctly passes these states while surfacing the non-passive branch. This is the audit doing its job across the corpus.

The constrained estimator (EstimationMethod.Feasible)

Three rings around the existing Gauss-Newton loop keep the estimate inside the admissible region. Every ring is inert for the other methods, so Wls, Huber, Lav and Constrained are byte-for-byte unchanged.

  • Ring A - interior guard (WlsStateEstimator.ApplyStep). Extends the existing trust-region damping with a fraction-to-boundary scale and a hard clamp so a voltage magnitude can never reach zero. This is the direct fix for "walked Vm through zero"; ApplyStep receives the feasibility spec only for the feasible method (null otherwise).
  • Ring B - active-set KKT (FeasibleWlsStep). Composes ConstrainedDenseWlsStep (reusing its zero-injection equalities, observability and covariance) and adds, each iteration, a linearised constraint row for every currently-violated state-coupled inequality: generator Q-capability Qg(x) in [Qmin,Qmax] (Jacobian = the injection-Q partials MeasurementModel.InjectionQRow), across-branch angle difference |theta_f - theta_t| <= limit (the two +-1 angle entries), and, opt-in, branch thermal rating |S_end| <= RateA. The thermal limit is nonlinear, so it is enforced through its squared form S^2 = P^2 + Q^2 with the row 2P*dP + 2Q*dQ from MeasurementModel.BranchFlowRow (the branch's own pi-model stamps, no copy-drift); only the more-loaded end is held (the two ends are near-collinear, so both would singularise the KKT), a small hysteresis band keeps the constraint active at the boundary, and its RHS is clamped so it can only hold or pull an over-rating flow down, never push a slack branch up onto its limit. Thermal enforcement is off by default - unlike a generator's reactive capability (a device limit), a thermal rating is operational and a real grid legitimately runs over it, so forcing the estimate onto the rating would distort a true over-rating state rather than report it (the overload is always surfaced on LimitViolations). The active set re-detects each iteration and converges with the state. Load-bearing invariant: with no equality and no active inequality the augmented system is the free normal equations, so the step matches DenseWlsStep to solver tolerance - the feasible method costs nothing on an already-feasible state. Dense only, under the same 3,000-bus fail-loud guard as Constrained.
  • Ring C - acceptance gate (WlsStateEstimator, post-convergence). Audits the converged state (evaluated at x-hat via ToPowerFlowResult - pure evaluation, no forward inference) and, on a hard-law breach, demotes it to EstimationStatus.Infeasible rather than shipping a fantasy. Operational limits (PhysicsFeasibility: generator Q-capability slack-exempt, thermal ratings, voltage band) are always reported on EstimationResult.LimitViolations but demote only under FeasibilityConstraints.RejectOnLimitViolation (a real grid legitimately runs branches over rating and buses outside the nominal band).

Constraints, spec and gate live in Estimation/FeasibilityConstraints.cs; EstimationStatus.Infeasible is the distinct loud terminal state.

Instruments

  • gridsim physics-check <case> - audit one solved state; prints the per-law table and the gate verdict (exit 1 on gate fail).
  • gridsim conformance - sweep the audit across the whole case corpus, writing a case × law matrix (results/conformance/conformance.{md,csv} + run-metadata.json, RunInfo-stamped, git-commit per CSV row). Every converged benchmark case is physically admissible.
  • Additive physics block in gridsim-solve-export/1 and gridsim-estimate-export/1 (schema unchanged; readers ignore unknown fields): the per-law report and gate verdict, so an estimated or presolved state carries its own conformance evidence.
  • tools/crossval/physics_audit.py - the primary instrument. Rebuilds Ybus and recomputes laws 1/3/6/8 in NumPy (reusing residual_check.build_ybus), independent of the C# math, and checks the exported gate verdict agrees with the independent recomputation. Verified agreeing to machine epsilon on real exports, including the negative-resistance case145 (both gates PASS; the passivity artefact is reported, non-gating). PhysicsExportContractTests locks the export block shape in CI.

Reproducing

gridsim physics-check case118                      # one state
gridsim conformance                                # the whole corpus -> a matrix
gridsim crossval case145 --tol 1e-4                # produce a solve-export
python tools/crossval/physics_audit.py results/crossval-case145/solve-export.json --tol 1e-4
gridsim estimate <case> --method feasible          # the physics-constrained estimate

Experimental models (gridsim experimental, Physics/Experimental/)

Opt-in, provenance-tagged, unvalidated-at-scale models live behind the experimental verb and off every validated path. Each follows the PropagationDelay discipline: a PROVENANCE block, off by default, and a degenerate limit that collapses to established physics (pinned by a test).

  • experimental reciprocity <case> (ReciprocityCheck, law 11) - passive-network reciprocity Y = Y^T. Phase-shifting transformers are the deliberate non-reciprocal exception (excluded); an asymmetry over the rest is a model/build defect, not an impossible state, so it never gates. Collapses to a trivial pass on a phase-shifter-free network.
  • experimental wave <case> (ContinuumWaveModel) - the grid as a damped-wave medium. The linearised swing on the network is a wave equation M d2(delta)/dt2 + D d(delta)/dt + K_L delta = P (inertia M = the 1/2 J w^2 density, damping D = the viscosity analogue, synchronising graph-Laplacian K_L = the elastic coupling). Its checkable content is the MODAL structure: the eigenvalues of K_L phi = w^2 M phi (via the self-contained SymmetricEigen Jacobi solver) are the network's inter-area electromechanical normal modes (~0.1-2 Hz), and the null mode is the centre of inertia. It is a linearisation of the validated nonlinear ClassicalMultiMachine swing - oscillation frequencies and mode shapes, not large-disturbance transients or turbulence - and collapses to the SMIB analytic frequency w_n = sqrt(w_s P_s / 2H) in its degenerate limit (test-pinned against the same oracle ClassicalMultiMachineTests uses). The dominant-mode damping ratio is the Reynolds analogue: < 1 underdamped (wave-like), >= 1 overdamped (diffusive). Needs per-machine inertia (fails loud through ClassicalInit on a case with no H).

On real GB data

The experimental and analysis models run end-to-end on the real GB spine (gb-spine, 316 real transmission buses assembled from NESO/DNO data, with the real ETYS boundary set). Headline results from one operating point:

  • Boundary attribution (experimental boundaries): B6 Scotland-England (SCOTEX) carries ~3,780 MW (limit 4,900), B8 North-to-Midlands (FLOWSTH) ~8,620 MW, SEIMP importing ~2,790 MW - each decomposed into the specific generators supplying the transfer. This is the north-to-south story, quantified per generator.
  • Carbon tracing (experimental carbon): a flow-based system intensity of ~181 gCO2/kWh with per-load intensities to ~335 gCO2/kWh, from the real per-generator fuel types.
  • DC-OPF (dcopf): ~53 GW of demand met at ~£257k/h as the provably cost-optimal security-constrained dispatch - a bounded-variable LP (balance + a two-sided thermal limit per rated line) solved to its global optimum by the self-contained BoundedSimplex. Cheapest (wind) plant is dispatched first, then backed off exactly where the network requires costlier local generation to hold every line inside its rating (the earlier marginal one-line-at-a-time redispatch under-secured, so this correct full-LP figure is higher than the old heuristic's).
  • Lake conformance (tools/crossval/lake_conformance.py): the physical-law audit over the stored parquet lake found 93/98 sampled gb-spine periods physically admissible to machine epsilon (~1e-9 MW), and flagged 5 anomalous periods - the sweep doing its job across years of real GB history. NOTE: the lake states are demand-scaled per period, so the audit compares S = V conj(Y V) to the STORED injection, not the static case's Pg-Pd.

The full gb-full case (3,539 real transmission buses) also loads and audits end-to-end: it converges and passes the hard-gate laws to machine epsilon (energy/KCL/S=VI* residuals ~4e-9, ~524 MW total losses). An earlier stale export carried two generator names with unescaped inner quotes that broke the JSON loader; both the C# GdaModelWriter and the Python writer.py emit through a JSON serializer that escapes correctly, so a regeneration is clean and the stale file was repaired in place.

Invariants a change here must preserve

  • Hard-gate laws {1,2,3,6,8} are the definition of "physically possible"; do not add operational preferences (band, ratings) to that set - they are reported, not gated.
  • The feasible method must remain inert for already-feasible states (FeasibleWlsStep == DenseWlsStep to 1e-12 with nothing binding) and for every non-feasible method (ApplyStep with feas == null is the old path).
  • The audit stays an independent recompute (per-branch pi-model, never the solver); the Python instrument is the cross-language backstop.
  • No forward inference: the gate consumes a ToPowerFlowResult projection, never a solver call.
Protections Protection assessment (a passive relay overlay)