Harmonics

The harmonic-analysis subsystem lives in src/GridSim.Core/Harmonics/. It answers a different question from the fundamental power flow - not "what is the 50 Hz operating point" but "what does the network do to the non-sinusoidal currents that power electronics inject into it" - and it answers it with the same machinery and the same discipline. The physics is grounded: the frequency-dependent admittance Y(h), the per-order penetration solve, the KCL that closes it, and the sequence rotation are all derived from the network model already in hand. The one assumed input is the harmonic-source spectrum - the current a converter, drive or rectifier emits at each order - and that input is provenance-tagged so an assumed-typical spectrum (CuratedPlaceholder + [VERIFY]) and a metered one (Recorded) are the same shape and upgrade in place. This chapter describes the design, the math as implemented, and the narrative the subsystem exists to tell.

The grounded engine and the swappable-source seam

Everything downstream of the source spectrum is as grounded as the fundamental power flow. Y(h) is built from the same branch impedances, shunts and transformer taps the fundamental YBus uses, scaled by the physics of frequency (below); the penetration solve is Kirchhoff's current law at each order (Y(h) V(h) = I(h), the same linear system the Newton step solves once linearised); the indices, sequence and compliance checks are arithmetic over the result. None of it is assumed. The seam is drawn at exactly one place: the harmonic-source spectrum, the injected current phasor I_h per bus per order. That is the only quantity the engine cannot derive from the network model, because it is a property of the connected equipment, not the grid.

The source set is therefore a first-class, versioned, provenance-tagged input (HarmonicSourceSet, loaded from a gridsim-harmonics/1 document by HarmonicSourceSet.Load). Each per-bus spectrum carries a ProvenanceTag:

CuratedPlaceholder (paired with a [VERIFY] note) - an assumed-typical spectrum synthesised from the built-in device library (a six-pulse drive, a PV inverter, an EV charger, ...). It is a defensible engineering default, not a measurement, and it is labelled as such everywhere it surfaces.
Recorded - a spectrum measured at the point of connection by a power-quality meter (GDA's PQ feed, or a third party's). This is a fact.
DerivedFromRecorded - a measured spectrum transformed (rescaled to a different operating point, aggregated over a feeder) but traceable to metered data.

Because the tag travels with the spectrum and every output echoes it, a study built on placeholders and a study built on metered data are structurally identical; swapping one bus from assumed to measured is an in-place upgrade that sharpens exactly that bus's contribution without touching the engine. This is the same discipline the evidence engine applies to event data (08-provenance-and-invariants.md): the physics is never in doubt, the assumptions are always named.

The complex-solver seam

Harmonic penetration is a complex linear solve - Y(h) and the injection I(h) are complex, the unknown V(h) is complex - so the subsystem introduces its own solver interface, IComplexLinearSolver, parallel to the real ISparseLinearSolver of the power flow (04-sparse-and-gpu-numerics.md). Three implementations sit behind it, cross-validated against each other so the seam is safe to switch:

Implementation	Method	Role
`DenseComplexLuSolver`	Dense complex LU with partial pivoting	The oracle: exact, small-case reference every other solver is diffed against.
`RealExpandedComplexSolver`	The `[G -B; B G]` real 2n-expansion, solved on the real `SparseLu`	The default: reuses the grounded, RCM-ordered sparse-direct factorisation the power flow already trusts.
`SparseComplexBiCgStabSolver`	Native complex preconditioned BiCGSTAB	The scalable native-complex path where the expansion's 2x width is unwelcome.

The RealExpandedComplexSolver is the default because it is the least new code on the most tested foundation. A complex system (G + jB)(Vr + jVi) = (Ir + jIi) is equivalent to the real 2n-by-2n block system

[ G  -B ] [ Vr ]   [ Ir ]
[ B   G ] [ Vi ] = [ Ii ]

which the existing real Gilbert-Peierls SparseLu factorises directly - same fill-reducing ordering, same stiffness-proof direct solve, no new numeric kernel. The dense oracle exists to prove the expansion (and the native complex BiCGSTAB) return the same phasors to machine epsilon on every case where all three run; the native solver exists so a large case need not pay the 2x dimension when an iterative complex method converges. Selection is by case size and solver availability, exactly as the real path chooses dense vs sparse Newton.

The frequency-dependent admittance `Y(h)`

FrequencyDependentYBus builds the network admittance matrix at harmonic order h (frequency h*f0). It is not the fundamental YBus scaled by a scalar; each element scales by its own physics:

Series reactance scales with frequency: X -> hX (an inductor's impedance is jhwL), while capacitive branch susceptance scales the same way, B -> hB (jhwC). This is the dominant effect and the reason resonances exist at all - inductive and capacitive reactances that cancel at one order do not cancel at h.
Skin effect raises series resistance: R(h) grows with frequency as current crowds to the conductor surface. The model is applied per equipment class - a distinct R(h) law for cable vs overhead line - because their skin behaviour differs materially, and each branch carries its own conductor type.
Shunts are frequency-dependent: fixed shunt reactors and capacitor banks scale their susceptance with h; this is where capacitor banks turn from benign reactive support at 50 Hz into the capacitive half of a parallel resonance at higher orders.
Optional load-shunt damping: loads can be represented as a frequency-dependent shunt admittance, which damps resonances - the loads are part of the grid's harmonic sink. This is optional because it is a modelling choice with its own assumptions, off by default so the undamped (worst-case) network is the baseline.
The distributed-parameter (long-line) correction: for electrically long branches at high orders the lumped pi-model understates the impedance. The hyperbolic correction Z' = Z * sinh(gl)/(gl) (with propagation constant g and length l) applies the exact distributed-parameter solution of the transmission-line equations, and the matching tanh(gl/2)/(gl/2) factor on the shunt legs. At low orders gl is small and sinh(gl)/(gl) -> 1, so the correction vanishes; at high orders on long lines it is the difference between a plausible and an implausible resonant frequency.

The load-bearing invariant: at h = 1 every one of these reductions is the identity, and FrequencyDependentYBus returns exactly the fundamental YBus. hX = X, hB = B, R(1) is the DC/fundamental resistance, sinh(gl)/(gl) evaluated at the fundamental gl matches the lumped model to the same order the power flow already uses. This is checked directly and is what lets the harmonic engine claim the same grounding as the fundamental one: it is the same matrix, continued to h > 1 by the physics of frequency.

Linear penetration, the impedance scan, and the sequence rotation

HarmonicPenetration is the core solve: for each order h present in the source set, build Y(h), assemble the injection vector I(h) from the provenance-tagged spectra, and solve Y(h) V(h) = I(h) through the complex solver seam. The result is the harmonic voltage phasor at every bus at every order - the penetration of each source's emission into the network. It is linear and superposition-based: each order is independent, because a linear network does not couple orders (the coupled, nonlinear case is below).

HarmonicImpedanceScan computes the driving-point (self) impedance Z(h) = 1 / Y(h)_ii seen looking into a chosen bus, swept over a fine grid of orders (or fractional orders). Peaks in |Z(h)| are parallel resonances - the frequencies at which a small injected current produces a large voltage, because network inductance and capacitance resonate. The scan identifies each resonance's order, magnitude and sharpness; a sharp, tall peak near a characteristic order (the 5th, 7th, 11th, 13th of six-pulse converters) is the dangerous case, because that is where real injection lands.

HarmonicSequence applies the phase-sequence rotation. In a balanced three-phase system a harmonic of order h rotates as a positive-, negative- or zero-sequence set according to h mod 3: h mod 3 == 1 -> positive (the 1st, 7th, 13th), h mod 3 == 2 -> negative (the 5th, 11th), h mod 3 == 0 -> zero (the triplens: 3rd, 9th, 15th). This is why the 5th appears as a negative-sequence quantity and the triplens sum in the neutral. The rotation is a property of the order, not the network, and it drives how each order couples through transformers (a delta winding traps zero-sequence triplens) and how it must be handled in the three-phase penetration below.

Indices, compliance, and the harmonic-damping-loss metric

HarmonicIndices reduces the per-order voltage spectrum to the standard engineering scalars:

THD (total harmonic distortion) per bus - the RMS of the harmonic content as a fraction of the fundamental, the headline number.
K-factor - the transformer de-rating weight sum (h^2 * I_h^2) / sum I_h^2, which quantifies the extra eddy-current heating harmonics impose on a transformer.
Damping sharpness - the resonance quality factor derived from the impedance scan: how narrow and how tall each |Z(h)| peak is. A high sharpness means a lightly damped, high-Q resonance.

HarmonicCompliance checks the computed distortion against the two limit regimes GridSim carries: IEEE 519 (the individual-harmonic and total-distortion voltage limits by system voltage level) and the GB Engineering Recommendation G5/5 planning levels. It reports pass/margin per bus per limit, so a study ends in a verdict, not just a spectrum.

The on-mission metric is harmonic-damping-loss: the resonance sharpness weighed against synchronous penetration - how far into the network the resonance reaches and how lightly it is damped as the synchronous fleet (the grid's harmonic sink) shrinks. It is the harmonic analogue of the inertia headroom the frequency subsystem tracks, and it is the number that carries the narrative in the closing section.

Three-phase and unbalanced harmonics

ThreePhaseHarmonicPenetration runs the penetration on the multi-phase network model (Distribution/MultiPhaseModel) rather than the positive-sequence equivalent, so it captures unbalance and the sequence structure explicitly. On a balanced network the results collapse to what HarmonicSequence predicts - the triplens are pure zero-sequence, the 5th is pure negative-sequence - which is the cross-check that the three-phase build is correct. On an unbalanced network (the realistic distribution case) the orders mix sequences, and the three-phase solve is the only one that gets the neutral currents and the per-phase distortion right. This is the harmonic side of the distribution and multi-phase work in 18-real-time-distribution-and-dynamic-estimation.md.

Coupled iterative harmonic power flow

Linear penetration treats each source as a fixed current injection. A real converter's emission depends on the voltage at its terminals - a nonlinear, order-coupling relationship - and where that matters, CoupledHarmonicPowerFlow iterates. Its structure:

The converter is modelled as a Norton equivalent at each order: a current source in parallel with an admittance. The admittance is folded into Y(h) (it stiffens the bus, and correctly damps or sharpens the local resonance), and the current source is the nonlinear injection.
The nonlinear injection is iterated: solve the penetration with the current estimate of the injections, recompute each converter's emission from the updated terminal voltage through its coupling model, and repeat to a fixed point.

The coupling behaviour is supplied through IHarmonicCouplingModel, with NortonCouplingModel the built-in implementation. Note the seam: the Norton admittance and the voltage-to-emission law are the one assumed part of the coupled solve - they describe the converter's internal control, which the grid model does not contain - so NortonCouplingModel's parameters are CuratedPlaceholder until a specific converter's data is supplied, exactly like the source spectra. The network physics folding it in is grounded; the converter model folded in is the assumption, and it is labelled.

Harmonic state estimation

HarmonicStateEstimator closes the loop the same way the fundamental estimator does (06-state-estimation.md): given sparse measurements, recover the drivers. It runs a linear complex weighted least-squares per order - because at a fixed order the network is linear, the harmonic estimation problem is a single complex WLS solve, not the iterative Gauss-Newton of the nonlinear fundamental case. From a set of sparse HarmonicVoltageMeasurements (metered harmonic voltage phasors at a handful of buses) it estimates the source injections I(h) that best explain them under Y(h). This is the estimator identity carried into the harmonic domain: the grounded Y(h) plus a few measurements pins down the many unknown injections, turning a sparse PQ-meter deployment into a network-wide picture of who is emitting what. It is the mechanism by which Recorded provenance propagates - a measured voltage at one bus upgrades the inferred spectra at the sources that feed it.

The narrative: the retiring harmonic sink

Synchronous machines are the grid's harmonic sink. Their low, largely resistive sub-transient impedance damps resonances and provides a stiff path that absorbs harmonic current. As they retire and are replaced by converter- interfaced generation, two things happen at once: the sources of harmonics multiply, and the damping that used to absorb them disappears. The impedance scan's peaks grow taller and narrower - resonances sharpen - and the harmonic-damping-loss metric rises. This is the under-told half of the low-inertia story: the same fleet retirement that erodes frequency inertia (05-frequency-dynamics.md) also erodes harmonic damping, and a low-inertia grid is a lightly-damped, resonance-prone grid as well as a fast-RoCoF one. The harmonic subsystem exists to make that second consequence as measurable as the first: run the same case with the synchronous fleet at 260, 155, 120, 102, 50 GVA.s of inertia and watch the resonances sharpen alongside the RoCoF.

Future extension: supraharmonics

The subsystem models the classical harmonic range (integer orders up to the low tens). The supraharmonic band (2-150 kHz), where switching-frequency emissions of modern converters live, is a documented future extension: it needs the same Y(h) continued to far higher frequencies (where the distributed-parameter correction and skin model dominate everything) and a source library at those frequencies. The seams - IComplexLinearSolver, FrequencyDependentYBus, the provenance-tagged source set - are built to carry it without redesign.