Technical manual — chapter list

Frequency dynamics

GridSim models system frequency with a centre-of-inertia (COI) swing equation: after a power imbalance deltaP, frequency moves at

df/dt = f0 / (2E) * deltaP        E = sum H*MBase  (spinning kinetic energy, MW.s)

The initial rate of change of frequency (RoCoF) is set entirely by the inertia E and the size of the loss; governor droop and load damping then arrest the fall, producing a nadir and a settling point. This chapter documents the model exactly as implemented: the batch simulator, its configuration, the stateful integrator, the pluggable inertia-model seam and the distributed multi-domain model behind it, the swing-equation inversion used for inertia monitoring, the GB frequency limits, and the NESO trace path.

All of it lives under src/GridSim.Core/Dynamics/ and src/GridSim.Core/Frequency/.

The batch simulator: FrequencyDynamics.Simulate

Dynamics/FrequencyDynamics.cs:FrequencyDynamics.Simulate takes a system inertia (MW.s), a base demand (MW), a set of FrequencyEvent(AtSeconds, DeltaPMw) power steps (negative MW = a generation loss) and an optional DynamicsConfig, and returns a FrequencyResponse. The integration is explicit forward Euler on a fixed timestep (StepSeconds, default 20 ms - deliberately the same cadence the NESO upsampler produces, so simulated and recorded traces share a timeline). Per step, at time t:

  1. Event forcing - every event with AtSeconds <= t contributes its deltaP (steps are cumulative and permanent; there is no event decay).
  2. Governor droop - a first-order lag toward the droop target -(deltaf/f0)/DroopPu * Presp, where Presp = ResponsiveFraction * baseDemandMw is the droop-responsive capacity. The lag time constant is GovernorTimeConstantS.
  3. Device fast response (FFR) - zero inside FastResponseDeadbandHz; outside it, a proportional target -(FastResponseMw/0.5)*deltaf clamped to +/-FastResponseMw (full output at 0.5 Hz deviation), lagged with time constant max(dt, FastResponseTimeConstantS).
  4. Recorded system response services - one independent first-order state per SystemResponseService(Mw, DeadbandHz, FullDeliveryDeviationHz, TimeConstantS): the target is a droop line reaching full output at the service's own delivery deviation (DC 0.5 Hz, DM/DR 0.2 Hz per the published service terms), zero inside the deadband, capped at +/-Mw, computed by the general SwingKernel.FastResponseTarget overload. These model the response the REAL system held (volumes from the EAC auction record), not a device - the replay keeps them where it strips device terms, and counterfactuals keep them too. A single droop line omits DC's low-gain (5 %) band between 0.015 and 0.2 Hz, which is acceptable for excursions well past 0.2 Hz where delivery sits on the steep segment.
  5. Load damping - -d*deltaf with d = DampingPctPerHz/100 * baseDemandMw (MW per Hz).
  6. Swing update - `rocof = f0/(2E) * (pEvent + pGov + pTsgb + pServices
    • d*deltaf), then f += rocof * dt`.

The inertia is E = max(1, systemInertiaMwS + DeviceRotationalInertiaMwS): the device's intrinsic stored energy (real (1/2)J*omega^2) is summed into E and cuts RoCoF directly, while its fast response is a separate power injection that arrests the nadir. The two device terms are never merged - the deprecated SyntheticInertiaMwS alias exists only to forward to DeviceRotationalInertiaMwS and is marked [Obsolete].

The returned FrequencyResponse carries:

Field Meaning
InitialRocofHzPerS the largest-magnitude per-step RoCoF seen (the post-event spike)
NadirHz, NadirSeconds the lowest frequency and when it occurred
SettleHz the frequency after the final update
Points the full 20 ms trajectory as FrequencyPoint(Time, Hz, RocofHzPerS)

Trajectory timestamps are synthetic, anchored at 2000-01-01T00:00Z - the simulator never reads the clock.

DynamicsConfig fields and defaults

Dynamics/FrequencyDynamics.cs:DynamicsConfig is an immutable record. Its calibration status is part of the contract and is documented on the type itself:

Field Default Status
NominalHz 50.0 fixed
DroopPu 0.05 5 % governor droop - GB Grid Code, documented, not fitted
DampingPctPerHz 1.0 load self-damping, % of demand per Hz - literature-typical, unfitted
ResponsiveFraction 0.20 share of demand carrying droop headroom - literature-typical, unfitted
GovernorTimeConstantS 8.0 primary-response lag (s) - literature-typical, unfitted
StepSeconds 0.02 integration step; matches the NESO upsampler
DurationSeconds 30.0 simulation window (the replay validator extends it, see chapter 7)
DeviceRotationalInertiaMwS 0.0 device intrinsic (1/2)J*omega^2 stored energy, summed into E
FastResponseMw 0.0 device FFR cap (MW)
FastResponseTimeConstantS 0.5 FFR ramp lag (s)
FastResponseDeadbandHz 0.015 FFR engagement deadband (Hz)
InertiaModel LumpedScalar which inertia model to run (see below)
Distributed null DistributedConfig when the distributed model is selected
SystemResponseServices [] recorded DC/DM/DR the real system held (see step 4 above); populated from the event record by GridEvent.ApplyRecordedServices
InitialFrequencyHz null recorded pre-event frequency. When set, the sim starts in equilibrium there: governor + service lag states pre-loaded to their static targets at the offset, and a constant standing-imbalance term balances them so df/dt(0) = 0 (the deficit the BM was carrying). Device terms start undeployed. Null = nominal start, all states zero

The unfitted constants (damping, responsive fraction, governor lag) set the nadir and recovery shape but not the initial RoCoF, which depends only on E and deltaP. Their uncertainty is roughly a factor of two; the scenario engine's best/worst corners bracket exactly these three. The bundled NESO sample (data/neso/frequency-sample.csv) is ambient data (+/-0.1 Hz band) - it pins the RoCoF scale but cannot calibrate the nadir constants; treat simulated nadir/settle values accordingly until a real loss-of-infeed trace is attached.

The stateful integrator: FrequencyIntegrator

Dynamics/FrequencyIntegrator.cs:FrequencyIntegrator is the running-clock twin of Simulate for the live dashboard: identical per-step physics, but advanced one Step(dtSeconds, eventMw) call at a time with governor, damping and FFR state carried between calls. An arbitrary deltat (1 ms to 20 minutes) is split into internal sub-steps of at most 20 ms (MaxInternalDt), bounded at 4,000 sub-steps per call, so explicit Euler stays stable - a huge step settles instead of blowing up.

The shared per-step terms - governor target and FFR target - are factored into Dynamics/SwingKernel.cs:SwingKernel (GovernorTarget, FastResponseTarget) so FrequencyDynamics, FrequencyIntegrator and the fleet kernel (FleetDynamics, chapter 7) cannot drift apart.

The inertia-model seam: IInertiaModel

Abstractions/IInertiaModel.cs:IInertiaModel is the swappable frequency/inertia model seam, mirroring IPowerFlowSolver:

  • Simulate(model, baseDemandMw, events, config) - returns the COI-collapsed FrequencyResponse, so every scalar consumer works with either model;
  • CoiInertiaGvaS(model) - the aggregate every scalar readout expects;
  • LastDistributedResult - the spatial detail of the most recent solve, or null for the lumped model.

Dynamics/InertiaModelSelector.cs:InertiaModelSelector.Create resolves the InertiaModel enum to an implementation. Unlike the power-flow selector there is no shape-based auto-selection: the inertia model is an explicit user choice, defaulting to InertiaModel.LumpedScalar.

Dynamics/LumpedScalarInertiaModel.cs:LumpedScalarInertiaModel is a thin adapter: FrequencyDynamics.Simulate(model.SystemInertiaGvaS * 1000, ...) - selecting it is exactly the classic behaviour.

The distributed multi-domain model

Dynamics/Distributed/ holds the experimental spatially-resolved model. It is additive and opt-in, and its outputs carry ProvenanceTag.UnvalidatedAtScale semantics: self-consistent, passing the reduction gates, but with no GB-scale field validation (that awaits PMU/DFR data). See chapter 8 for the tag.

Domain assignment

Distributed/DomainAssignment.cs:DomainAssignment.Assign partitions the model's generators into speed domains - groups sharing a mechanical speed/pole ratio, treated as one coherent inertial body:

  • Distributed/PoleMap.cs:PoleMap.ForType maps the generator type string to (poles, rpm, kind): turbo-alternators (nuclear/coal/biomass/ccgt/ocgt/chp/ gas) are 2-pole 3,000 rpm synchronous; recips 4-pole; hydro/pumped 8-pole 750 rpm; storage/battery/BESS become SyntheticInertia (virtual) domains; everything else inverter-based is Pll. This mapping is a documented modelling assumption, not per-unit data - GDA holds no pole/speed field.
  • A machine with H > 0 but an unknown type defaults to a 2-pole synchronous unit (the WSCC case9 generators land here).
  • Machines are grouped by (poles, kind); each domain aggregates E_d = sum H*MBase (MW.s), rated MVA, and GVAd = sum MBase*C where the availability/coupling coefficient C = clamp(|Pg|/MBase, 0, 1) is a dispatch-derived proxy.
  • Distributed/SpeedDomain.cs:SpeedDomain.IsSwingMass - a domain participates in the electromechanical swing only if it is Synchronous with E > 0; virtual domains contribute to dominance accounting only.

The dominance read-out (Distributed/DominanceReport.cs:DominanceReport) orders domains by GVAd and computes deltaGVa = (GVAd1 - GVAd2)/GVAd2; a dominant domain exists when deltaGVa >= DominanceThreshold (default 0.10). The 0.10 default is an engineering-judgement value, not a theory-derived or fitted constant - but it is empirically anchored: over 103,433 GB samples (GDA Derived/speed_domains, 2014-2026) the dominance gap has P10 ~ 0.106, so 0.10 flags almost exactly the bottom decile of conditions (fires ~9.4% of the time; 0.08 -> 7.5%, 0.12 -> 11.4%). It is init-settable for calibration.

The multi-domain swing and its invariants

Distributed/MultiDomainSwing.cs:MultiDomainSwing.Simulate integrates one swing equation per swing-mass domain, coupled by tie synchronising power K_ij*(delta_i - delta_j). The construction guarantees two reduction properties:

  • T0 (exact reduction). With <= 1 swing mass - including DistributedConfig.SingleDomain = true - DistributedInertiaModel delegates to the lumped FrequencyDynamics, so the degenerate case is bitwise identical to the lumped model, not merely close.
  • T1 (COI preservation). The coupling matrix is symmetric, so tie power is antisymmetric and sums to zero across domains: it moves energy between domains without moving the centre of inertia. The system-level forcing (event, governor, damping, FFR) is computed from the inertia-weighted COI frequency exactly as the lumped model does, then distributed across domains in proportion to inertia - so the per-domain updates re-aggregate to the lumped COI trajectory (to rounding).

The localised loss lands on the domain hosting the largest unit by rated MVA (its credible-loss home). The device's rotational inertia is folded into the fleet proportionally so the COI total still equals the lumped E + DeviceRotationalInertiaMwS. Integration is semi-implicit (symplectic) Euler on the same 20 ms step - frequency first, then angle with the updated frequency - which stays energy-bounded under the stiff Kron-derived coupling where plain forward Euler would grow without bound.

Coupling comes from Distributed/DomainCoupling.cs:DomainCoupling.SynchronisingMatrix: the bus admittance matrix is Kron-reduced onto the generator buses and the cross-domain synchronising susceptance read off as K = BaseMva * Im(Y_reduced) at flat 1 pu voltage (anti-synchronising artefacts <= 0 are ignored). Models with no branch reactance fall back to a uniform stiffness targeting DistributedConfig.TargetInterAreaHz (default 0.5 Hz, the GB inter-area mode band).

Propagation delays

Distributed/PropagationDelay.cs:PropagationDelay.Compute supplies the inter-domain delay matrix tau. It is off by default (DistributedConfig.EnablePropagation = false) because the retardation perturbs the centre of inertia (real physics), so the machine-eps COI invariant only holds with it off. When enabled:

  • domain centroids are computed from the members' Bus.X/Y map coordinates, scaled by the GB bounding-box spans (660 km x 1,020 km);
  • the wave speed is derived from inertia density (v ~ sqrt(rho_ref/rho), reference 200 MW.s/km -> 1,000 km/s) and clamped to the published electromechanical-wave band 500-1,500 km/s (MinWaveSpeedKmPerS/MaxWaveSpeedKmPerS) - the clamp is the honest statement: only the band is asserted, not the absolute GB value;
  • tau_ij = distance/v, unless DistributedConfig.ProvidedTauSeconds supplies an externally-derived tau for the domain pair's representative buses (GDA's propagation_delays lake product, fetched over the bridge by Net/GridBridgeClient.cs:GridBridgeClient.QueryPropagation).

In the swing loop each tie reads the far domain's angle retarded by tau via a per-step angle history; tau = 0 reduces exactly to the undelayed model. A negative/NaN/infinite tau collapses to zero delay and any finite tau is clamped to the run length, so a bad value cannot index outside the history.

Distributed outputs

Distributed/DistributedResult.cs:DistributedResult rides on IInertiaModel.LastDistributedResult: per-domain DomainTraces (same shape as FrequencyResponse, one per coherent mass), the DominanceReport, a per-bus event-RoCoF map, and - when propagation is on - a PropagationReport (wave speed used, largest delay, per-domain arrival times, where arrival = the first departure of more than 1 mHz from nominal).

The authored per-asset fleet: pAssetFleet

The swing loop is additive in power, so the inverter fleet enters as one more term. Dynamics/AssetFrequencyDynamics.Simulate is a faithful copy of FrequencyDynamics.Simulate that appends, last in the imbalance sum, the net power of the authored asset fleet and the cumulative Loss-of-Mains trip deficit:

pImbalance = pEvent + pStanding + pGov + pTsgb + pServices - d*df  +  pAssetFleet - assetTripDeficit

Assets/AssetDynamicsFleet advances each bound controller (Assets/GridFollowingInverter, GridFormingInverter) against its local frequency (an IGridFrequencySampler - the COI for the lumped model). A predictor-corrector samples this step's predicted RoCoF, so a grid-forming unit's synthetic-inertia term (p = -synInertia*rocof) lowers the initial RoCoF - as a power injection, never merged into E. A grid-following unit trips on local RoCoF/vector-shift/LVRT; its lost dispatched Pg becomes deficit exactly as RocofScreening treats a lost unit.

Gated by DynamicsConfig.EnableAssetDynamics (default off). With an empty/passive fleet both added terms are 0.0, so x + 0.0 - 0.0 == x makes the trajectory bit-identical to the baseline - the reduction gate (AssetDynamicsReductionTests). The live twin is AssetAwareFrequencyIntegrator. Because every inverter control parameter is unmeasurable ([VERIFY]/CuratedPlaceholder), any real-fleet claim uses a sweep envelope (AssetEnvelope) tagged UnvalidatedAtScale, validated against a recorded trace (AssetValidation) - not a point curve. See docs/manual/18-extending-gridsim.md for authoring classes.

Inverting the swing equation: InertiaEstimator

Dynamics/InertiaEstimator.cs:InertiaEstimator recovers inertia from an observed excursion - what NESO does after every large loss:

E = f0 * |deltaP| / (2 * |RoCoF|)
  • EstimateInertiaMwS(deltaPMw, rocofHzPerS, nominalHz = 50) - signs are ignored; returns +inf when |RoCoF| < 1e-9 (no observable slope);
  • EstimateInertiaGvaS - the same divided by 1,000, the GB reporting unit;
  • EstimateHSeconds(..., baseMva) - the equivalent H = E/S on a given base;
  • PredictRocof(inertiaMwS, deltaPMw) - the forward check, with E floored at 1 MW.s.

The estimate is clean because the very first RoCoF sample precedes any governor or fast-response action: at t = 0 the only imbalance is deltaP itself. The replay validator uses this inversion as its inertia cross-check (chapter 7).

GB limits and assessment: FrequencyLimits

Dynamics/FrequencyLimits.cs:FrequencyLimits holds the operating thresholds and two classifiers.

Constant Value Meaning
OperationalLowHz / OperationalHighHz 49.8 / 50.2 NESO's operational band
StatutoryLowHz / StatutoryHighHz 49.5 / 50.5 statutory +/-0.5 Hz
LfddFirstStageHz 48.8 low-frequency demand disconnection, first stage
RocofLimitLegacyHzPerS 0.125 pre-ALoMCP loss-of-mains relay setting
RocofLimitModernHzPerS 1.0 GB relays as retuned by ~2022
RocofLimitHzPerS 0.125 back-compatible alias - pinned to legacy
ActiveRocofLimitHzPerS 1.0 (settable static) the limit the "active" flag scores against

RoCoF assessment is dual-threshold: every assessment reports the legacy and modern breaches independently, plus an "active" flag scored against ActiveRocofLimitHzPerS (settable, defaults to modern). The retained RocofExceedsLimit/RocofLimitExceeded fields keep legacy (0.125 Hz/s) semantics for back-compatibility. BreachesRocof compares with a 1e-9 Hz/s tolerance.

  • Assess(FrequencyResponse) -> FrequencyAssessment: the headline RoCoF/ nadir/settle, the limit flags (BelowOperational, BelowStatutory, LfddAtRisk at nadir <= 48.8 Hz) and a one-line summary chosen from a severity ladder (LFDD > statutory > active RoCoF > legacy RoCoF > operational > within limits).
  • AssessTrace(IReadOnlyList<FrequencySample>) -> FrequencyTraceReport: scans a measured trace - extremes, sample counts below the operational and statutory floors, whether 48.8 Hz was reached, and the worst sample-to-sample RoCoF |deltaf/deltat| against the same three thresholds.

The NESO trace path and RocofUpsampler

The measured-frequency types live in Frequency/FrequencyTypes.cs: FrequencySample(Time, Hz) is one recorded reading; FrequencyPoint(Time, Hz, RocofHzPerS) is a point on the fine timeline.

Frequency/NesoFrequencyCsv.cs:NesoFrequencyCsv.Read streams a NESO system-frequency CSV (dtm,f) row by row, so a month of 1-second data (~2.7 M rows) reads without loading into memory. Unparseable lines are skipped. Two details are load-bearing:

  • Timestamps land in UTC. Offset-carrying stamps are converted to UTC and naive stamps are assumed UTC (AssumeUniversal | AdjustToUniversal). Without this, the replay validator's RMSE alignment - which compares recorded times to the UTC event instant by raw ticks - drifts by the local offset during BST, enough to demote a Reproduced baseline to Marginal.
  • Every row passes ForwardInferenceGuard.RejectFuture - the reader is structurally incapable of ingesting a recent or future sample (chapter 8).

Frequency/RocofUpsampler.cs:RocofUpsampler refines the coarse trace:

  • Upsample(samples, subSamplesPerInterval = 50) - between consecutive samples it computes the interval RoCoF (f2-f1)/deltat and fills the gap by linear interpolation, f(t) = f1 + RoCoF*t, emitting one FrequencyPoint per sub-sample (50 per 1-second interval -> 20 ms cadence, the resolution the dynamics model and an FFR controller want). A closing point repeats the final sample with RoCoF 0.
  • Rocof(samples) - the per-interval RoCoF series, stamped at the later sample.
  • Summarize(samples) -> FrequencyStats (count, min/mean/max, worst |RoCoF| and when).

Consumers of this path: the CLI freq verb (summary + limits scan + an upsample demonstration), the CLI replay verb (upsampled 20 ms playback via Frequency/ReplayController.cs), and ReplayValidator's data-gated per-sample RMSE, which upsamples the recorded trace to the simulation timestep before comparing (see chapter 7).

See also

Sparse and GPU numerics State estimation