Simulating the GB power system frequency during underfrequency events 2018–19

Lightning hit a transmission power line outside London, England on 9 August 2019. There followed a loss of power from a cascade of generator outages that exceeded contingency reserves, leading to an exceptional fall in grid frequency causing wide-spread transport disruptions and the disconnection of over 1 m households. Simulating such events typically involves a system of differential equations representing the overall generation and load present at the time. A standard model based on the swing equation assumes unlimited capacity in aggregated resources, and the availability of these services throughout the duration of the frequency deviation. In simulating the effect of outages on the GB Grid frequency on 2019/8/9, the effect of limiting these services to the capacity of resources engaged during the event is examined. It is shown that by taking these refinements into account the timing and extent of the frequency nadir can be successfully estimated. Insight is gained into the responses of various grid characteristics and how they interact with unplanned generation imbalances. Using this adapted model, further events on the GB grid are examined to validate the influence of these features. With the model’s effectiveness validated, novel mitigating measures to preserve the stability of a low-inertia grid can be evaluated.


Introduction
The largest outage on the Great Britain (GB) grid in recent years occurred on Friday August 9, 2019 (referred hereafter as the 2019/8/9 event).Lightning struck an overhead power line between Eaton Socon and Wymondley at 16:52.33 BST [1].Two major generation units went offline almost immediately, followed by a cascade of outages that led to a cumulative power loss of over 1900MW .This exceeded the capacity of the reserves held to maintain the integrity of supply.As a result, the power frequency, normally 50Hz , dropped to 48.8Hz for the first time in over a dec- ade.This triggered exceptional measures intended to preserve the overall network stability leading to 1.15 m households being disconnected.
Conventional thermal generation has traditionally provided resistance to instability through inertia: synchronous turbines converting between mechanical and electrical energy to dampen the effects of sudden drops and surges in generation.However, stability of the power supply on the GB grid has been evolving in recent years due to the emergence of renewable generation as a significant contributor to the energy mix [2].These generation sources add to the power capacity of the network but not significantly to inertia, leading to a decrease in the aggregate system inertia over time.
This paper sets out a frequency model to examine current and future grid stability using the GB network as a basis.The model 1. uses a dynamical system to represent the principal influences on the change in frequency of the grid power supply during a transient event.2. incorporates (a) Frequency Response (FR) that reacts to changes in power and frequency; (b) Frequency Response that responds over different timescales; (c) Frequency Response of limited duration; (d) the response of grid load to changes in frequency.3. is verified by the 2019/8/9 event simulation and is then applied to a further 26 under frequency events recorded National Grid Electricity System Operator (ESO) in the period 2018-19 [5].
The model can be used to evaluate variations in Frequency Response services and to simulate future mitigation measures that may smooth the transition to a stable, lowinertia grid.It will be used as the basis of a networked system, allowing the examination of inertia distribution and Frequency Response in current or future scenarios, in addition to the influence of increased Inverter-based resource (IBR) generation.The organization of this paper is as follows.Section 2 sets out the dynamic model used to simulate the 2019/8/9 event and the configuration of parameters necessary to recreate the grid's state following the lightning strike.Section 3 compares the simulated frequency with the recorded trace for the event.Section 4 applies the model to a suite of other events over the period 2018-19 to ensure that it is generalisable, and Sect. 5 discusses the results of simulating these events using the model.In Sect.6 conclusions are drawn from the results of the study and further avenues for research are outlined.
2 Modelling frequency response in 2019/8/9 event Examining the frequency trace of the 2019/8/9 event (Fig. 1) the sequence can be divided into five stages, each beginning with an occurrence on the grid.
1. 16:52:34 Outages at Little Barford power station and Hornsea windfarm follow the lightning strike, with further disconnections from Loss-of-Mains protections at embedded generation.2. 16:52:58 Primary Frequency Response (PFR) halts the drop at 49.1 Hz , and the frequency starts to recover.3. 16:53:30 A second outage at Little Barford precipitates a further fall in frequency below 49 Hz , triggering frequency-limit disconnections in embedded generation.4. 16:53:48 Low Frequency Demand Disconnection (LFDD) arrests the fall in frequency, which overcomes a further outage at Little Barford to recover. 5. 16:55:00 Frequency recovers to levels above the 49.5 Hz statutory threshold [7].
The modelling objective is to devise a high-level, lightweight, configurable model to simulate this and other events.To achieve this, an aggregate single bus DC model of the system centre of inertia is adopted, which has been shown to adequately Fig. 1 2019/8/9 event trace from 1 s historic dataset provided by ESO [9].Events marking the beginning of stages 1-5 are indicated as numbered lines represent frequency fluctuations [8] [6].In addition, as the model is intended to examine only significant perturbations, the effects of voltage, line losses, turbine response and deadbands are neglected as having a lesser effect than that of frequency response services.
The standard model consists of a 3-dimensional system of ordinary differential equations, based on the model detailed in [3].It has as its foundation the standard swing equation for u, the per unit frequency deviation from the nominal (see Table 3) given by where H SYS is the system inertia, p PFR is the Primary Frequency Response (PFR), p SFR the Secondary Frequency Response (SFR) and p DR is the change in demand in response to the frequency deviation.Here p L (t) is the per unit load imbalance at time t.Each of these factors are detailed in the following sections, which first describe the standard model from [3], and then set out the changes in the model that have been introduced to improve the agreement between simulation and actual frequency for the 2019/8/9 event.

Load imbalance p L (t)
From the time of the event t = 0 at 16:52:34, the imbalances ΔP L (t) (in MW) occur- ring up to time t are aggregated (Table 1) to arrive at the total load imbalance P L (t) .In (1) this is expressed as p L (t) , a per unit value of the total demand at the event time -approximately 29 GW [1].
Following the example in the ESO model [10] the generation reduction at Stage 3 from 30s to 60s after the initial outage combined with the second Little Barford disconnection.The ESO records that 350 MW of LFDD occurred at 16:53:48 [10].This is treated as a net reduction in demand and deducted from the load imbalance from this time on. (1)

Inertia H SYS
The effect of inertia dominates in the initial seconds after the event, before Primary Frequency Response (PFR) and Secondary Frequency Response (SFR) have ramped up.At this point the standard model (1) at time t = 1 reduces to For the 2019/8/9 event p L (1) = 1561 MW 29 GW and the largest per second Rate of Change of Frequency (RoCoF) from the 1 s data at the time of the event [9] is ̇f (1) = 0.151 Hz s −1 .Given u = (f − f 0 )∕f 0 , where f is the frequency in Hz , f 0 = 50 Hz , so that ̇u = ̇f ∕f 0 , and the result H SYS ≃ 9.1 s is obtained.

Primary frequency response
PFR is intended to halt and stabilise the frequency before longer-term FR resources are ramped up (Fig. 2).The PFR value p PFR is governed by the follow- ing differential equation where p is the PFR time constant.The default response rate is p = 10 s (Table 3).This value is the principal determinant of the extent of the first nadir.For the 2019/8/9 event, it is adjusted to p = 13.3 s to align the simulation with the event curve at this point.This adjustment could be attributed to a delay in activation that is consistent with dynamic PFR provision.
All generators connected to the GB grid must be in a position to provide PFR, including IBRs [7].In the model PFR is aggregated into the single variable p PFR .Static PFR and Enhanced PFR are activated when the frequency crosses defined thresholds and have ramp times of 1s.These thresholds are crossed in the immediate aftermath of the events examined in this paper.These two services can therefore be adequately simulated by being aggregated into the overall PFR response.
The standard PFR model assumes unlimited capacity and indefinite availability, resulting in the allocation of resources significantly in excess of those that would be available.Equation ( 3) is therefore modified to equation ( 4) to take into account both the PFR capacity present at the time of the event (Table 2) and a limit on the time that PFR resources are at full capacity.
In (4) the total capacity of PFR is PPFR calculated in Table 2 as 831 MW , converted to the per-unit value.
(2) ̇u( 1 t PFR is the time from the start of the event until dynamic, short-term PFR is disabled.The nominal value is 60s (Fig. 2).As with the ramp rate, this is modified heuristically so as to match the observed behaviour of the event frequency trace.In the case of the 2019/8/9 event it is set to 45s.The parameter PFR is the proportion of total PFR which is static, remaining in service after dynamic PFR is terminated.The reported figure (Table 2) is PFR = 51% ( 424 MW of 831 MW ).To match the event trace, we must assume that a greater proportion of PFR remains in service after t PFR , so that PFR is set to 65%.

Secondary frequency response
For the standard model ( 1), the value for SFR, p SFR , is determined by the equation where K i is the secondary control gain with default value K i = 0.05 s −1 (see Table 3) and is the integral of the per-unit frequency deviation, so that θ = u .The parameter K i is calibrated so that the slope of the recovery in frequency between Stages 2 and 3 is matched.With PFR at capacity, it is the growth in SFR at this time that alters the net imbalance on the RHS of (1).The default value of 0.05 s −1 is thereby modified to 0.013 s −1 .As the capacity of SFR is not reached over the timescale of the simulation it does not need to be accounted for.

Demand response
Demand varies due to frequency, where some devices such as synchronous motors use slightly less power when frequency is low [11].This Demand Response (DR) reduces the load from grid-connected rotating machinery.This effect is assumed to correspond to be 2.5% Hz −1 [12], so that As this is an effect that is not directly measurable, this assumption is the subject of debate [13]; nevertheless it has been included in the model.

Turbine governor and Control Room response
The generator turbine mechanical response to frequency imbalances was implemented according to [3].With limited capacity its effects during frequency excursions were found to be reproducible through an increase in the capacity of other FR measures and it was therefore omitted from the model.
Tertiary response measures, which are activated manually, are not detailed in the event reports from the ESO.It is reported in [1] that the frequency returned to 50 Hz by 16:57:15 as a result of FR measures and 1240 MW of Control Room actions.These Control Room actions are reported as instructions occurring after the ( 5) last outage at Little Barford at 16:53:58, the last event modelled in the simulation, between Stages 4 and 5.The inclusion of these actions would see the modelled trace return to the nominal frequency more rapidly than with the support of SFR alone.Indeed the frequency after the event exceeded 50 Hz for a time, which may be a con- sequence to the Control Room interventions.Fig. 2 Characteristic frequency trace during a system imbalance (based on [5]) showing a steady state variation within the statutory limits ( 49.8 Hz to 50.2 Hz ) followed by a sudden fall as a result of a disturbance.Recovery in frequency after the nadir plateaus after 60 s, returning to steady state after 30 min Parameter and variable values are summarised in Table 3. PFR increases to capacity pPFR between Stages 2 and 3, decreasing after t PFR to the level of static PFR (  PFR pPFR ) before Stage 3. SFR increases gradually over the 3 min of the simulation.DR inversely tracks changes in frequency, decreasing as it recovers.
The DR magnitude is consistent with the expected value (c.630 MW at 49.1 Hz [14]).The Net Imbalance is the summed RHS powers in the swing equation, (1).The model simulation trace (Fig. 4) approximates the extent of initial nadir level (Stage 2) while being offset slightly in terms of its timing, validating the chosen value for p and the calculated H SYS .Frequency recovers from Stage 2, the slope of recovery determined by the value for K i setting the ramp rate for SFR.From Fig. 3 it can be seen that this is due to the net imbalance being above zero at Stage 2, until the timeout of dynamic PFR.Its decline is precipitated by the second Barford outage at Stage 3, falling rapidly as PFR capacity is reached and SFR insufficient to make up the deficit.Grid frequency is therefore falling.This second fall is closely matched (Stage 3), including the timing and depth of the second nadir (Stage 4).
LFDD at Stage 4, incorporated in the p L (t) time series (Table 1), returns the net imbalance to positive, and thereafter the frequency recovers to the statutory Overall the trace compares favourably with that of the ESO Frequency Simulation Engine [10].

Other GB grid underfrequency events 2018-19
The ESO published details of 34 over-and underfrequency events on the GB grid during the period 2018-19 [5].The information provided included the magnitude and location of the disturbances causing a frequency deviation and the grid characteristics at the event time.Using the corresponding 1s frequency data at the event times [15] for comparison, it is possible with this information to apply the model that simulates the 2019/8/9 event to the simulation of these other events.
As 25 underfrequency events other than that on 2019/8/9 involved a single outage, p L (t) is static, with the event occurring at t = 0 .The loss is expressed as a proportion of the generation present.The event time is taken to be where the RoCoF is greatest, indicated by an abrupt change from normal noisy variation (Fig. 2).One of the events reported by ESO has been omitted as having the highest proportion of inertial generation, 94.3% , with the result that no variation in FR affects the frequency simulation trajectory, suggesting that the deficit caused by the outage was remedied by standard turbine response of the generation active at the time.Section 2.2 sets out the process to calculate H SYS , using the total inertial and non-inertial generation at event time given in the ESO report.
Whereas the coefficients are arbitrary, as knowledge of the grid condition at the event time is limited, it is possible to fine-tune the model based on the common characteristics of perturbations.In calibrating the model, the PFR ramp rate, p , is calibrated by matching the simulation with the frequency nadir and SFR response factor K i is set so that the recovery path of the event trace after the nadir is matched by the simulated frequency.It was demonstrated in subsection 2.4 how, at this stage in the event aftermath, SFR is the principal determinant of the frequency's path up to this point, once all other factors have been taken into account.Each parameter is calibrated based on real-world measurements.Aggregate inertia is determined by the reported RoCoF, and validated against the initial drop in frequency.PFR is calibrated to achieve the correct nadir.During the initial phase of recovery, SFR is the only force changing at the time of calibration, and is therefore set to align with the event trace path for this stage.
The timing of dynamic PFR cut-off, which is necessary to accurately reproduce the general recovery trajectory, is clearly visible in all cases.In each of the event traces there is a clearly distinguishable point where the frequency changes, usually between 30s to 60s after perturbation.The frequency's trajectory alters, demarcating the boundary between different FR configurations, from which the value for t PFR is taken.The path of the frequency after this point is determined by the reserves available to support the recovery in frequency.With SFR set, the remaining variable is the amount of static PFR active after dynamic PFR has stopped.In modelling these under-frequency events it is assumed that the Fig. 4 Model simulation of 2019/8/9 event frequency trace compared to the actual frequency trace from [1] capacity of PFR is comparable to that available for the 2019/8/9 event, which is reasonable as, under the Grid Code [16], provision is to be made for PFR to mitigate the effects of an N − 1 event, a single generator outage of 1 GW .It is also assumed that there was the same proportion of delivery at the time of these events as was the case on 2019/8/9 (Table 2), taking into account the underperformance of FR providers [17], so that PFR capacity is 831 MW PPFR = 831MW .The value PFR for the proportion of total PFR that is static is determined by the behaviour of the frequency after the disconnection of dynamic PFR.With all other FR measures fixed or previously configured, the criterion that decides whether the frequency increases or decreases at this point is whether there is sufficient FR capacity to make the Net Imbalance positive or negative.This determines the sign of the RHS of (1) and thereby the increase or decrease in frequency.Varying PFR varies the amount of static PFR so that the path of the frequency can be matched.

24 underfrequency event simulations
In all instances of the 24 event simulations (Figs.5,7,9), the initial path of the frequency after the outage matches that of the event trace, indicating an acceptable approximation of H SYS .The p , as measured by matching the simulation nadir with that of the event, is in the range of 4.5sto21.5s,apart from two instances where variations in p have no effect on the nadir.These cases have the lowest loss per unit of generation at the time of the event, suggesting that normal governor response performed the role of frequency response in recovering from the shortfall.With the exception of these two cases, the total duration of dynamic PFR provision, t PFR varied from 23.1s to 59.5s (Figs.6,8,10), in agreement with the standard model (Fig. 2).In none of the illustrated cases was the capacity of PFR, pPFR , reached.In most cases it was possible to match the recovery in frequency from the nadir by calibrating K i .From the default of 0.05 s −1 , this had a range 0.0015 s −1 to 0.06 s −1 .In the remaining cases where frequency recovers rapidly, it is possible that the origin of the outage was remedied before FR ramped up sufficiently to affect frequency.
The threshold between dynamic and static PFR outlined in Sect.2.3 was readily identifiable in the majority of cases.From this point, three characteristic trajectories were identifiable, depending on the Net Imbalance after dynamic PFR had been disabled: 1. Smooth recovery (Fig. 5), where the Net Imbalance was greater than zero, so that available FR is sufficient to attain a stable frequency (Fig. 6).2. A kink in the path (Fig. 7), where the Net Imbalance dips below zero momentarily, preventing the frequency from falling significantly, but delaying recovery until SFR increases further (Fig. 8).3. A drop in frequency after the cut-off of dynamic PFR, to levels near to or exceeding the initial nadir (Fig. 9), where the Net Imbalance is markedly negative at t PFR and FR is insufficient to stabilise the drop in frequency until SFR increases to the point where it makes the Net Imbalance positive again (Fig. 10).
Matching each of these trajectories is done by the calibration of PFR to quantify the amount of static PFR available after the cut-off at t PFR , and thereby the FR resources to balance the deficit.In the 2018-19 event simulations, this value ranged from 10% to 91% , a significant variation in grid behaviour.Correspond- ence between the simulation and event traces reduce after the dynamic PFR cutoff point.Whereas the 2019/8/9 event was well documented in reports by ESO and Ofgem (the energy regulator) [1][10] [17], there are no details about secondary outages in the ESO report on frequency events.Such incidents were a significant contributor to the 2019/8/9 incident magnitude, and during the period 2018-19 many of the same Loss-of-Mains protection systems would have been in place, responding to sudden variations in voltage and frequency.It is therefore only possible to simulate the overall path of the frequency without this information.

Conclusions
In this paper the behaviour of grid frequency during the 2019/8/9 event and other underfrequency events on the GB grid in the period 2018-19 is simulated.For this a 3-dimensional dynamic model which accounts for characteristics of grid FR is devised.From this work the following conclusions can be drawn:   1. Frequency transients and FR services can be adequately modelled using a dynamic system so that the behaviour of the frequency in response to a Net Imbalance can be examined.2. In modelling major transient events, the capacity of FR reserves must be taken into account.Once this level is reached, FR must remain static for the remainder of the simulation, or for the contracted duration if this applies.3. The limit on the duration of the availability of dynamic PFR should be allowed for, and the proportion of total PFR available after this ends calibrated to the trajectory of frequency recovery.4. The inclusion of DR is important for accurate transient modelling, but the basis for simulating its frequency sensitivity during excursions should be scrutinised.
The new contributions of this paper are that an accurate and efficient modelling of frequency trace for the 2019-08-09 event can be achieved, once detailed information on the sequence of power outages and frequency response are available.The model is also useful in analysing other frequency events even when such information is not available, and can recreate event traces outside the boundaries of normal operation.It is also sufficiently abstract and configurable to represent other grids, and it can therefore be adapted and extended to form the basis of a wider investigation of the effects and mitigation of frequency events, in both current scenarios and planning for projected future configurations.Although the study is preliminary, and more complicated simulation models can reproduce a more accurate trace, the advantage of this model is its highlevel simplicity, focussing on the significant forces impacting frequency during a Fig. 10 Imbalance and Frequency Response graph for 2019/10/03 event.Net Imbalance is markedly negative at t PFR causing the frequency to fall a second time until SFR increases sufficiently contingency, and demonstrating that a low-dimensional model can compare with complicated, highly computationally intensive models.It is therefore amenable to be extended to examine future energy scenarios cost-effectively without being tied down to specific grid characteristics.The work is of interest to others in the field who are also interested in scenario-based research, where computationally expensive models would be prohibitive.
While it has been shown that a reasonable approximation of an event frequency trace can be achieved using realistic parameters, and the model's generality has been demonstrated by examining other under-frequency excursions on the GB grid, the model would benefit from extended analysis.Greater detail on the underfrequency events would allow a more in depth study of FR services.PFR and SFR disaggregation into their static and dynamic components (with varying time constants and service durations) would result in a closer correspondence between the simulation and the actual grid configuration.This had been investigated in preliminary work but was found to be unnecessary for a first approximation of frequency behaviour.The contribution of deadbands and generator turbine response to the accuracy of an event trace, also not included in the model after experimentation, should be re-examined.
The effects of geographic generation distribution, inertia, and FR should be examined, as well as testing the model on events in other grids.This will involve adapting the grid configuration on these grids at the time of these events and a comparison of the results.Effects such as inertia, DR and FR are configuration-dependent, and their effects on obtaining good simulation results in a variety of scenarios should be investigated.
To address the single-bus model's shortcomings and facilitate further investigations of the issues described above, the future path of this work is to expand the model to represent a network.This will extend the model functionality to facilitate the investigation of local variations in inertia and optimal remedial strategies for frequency response distribution.It will also allow the study to anticipate IBR generation, with the consequence of local and aggregate changes in inertia and increasing sensitivity to perturbations.

3 2019/ 8 / 9
event model simulation The change in system imbalance and frequency responses during the simulation of the first 3 min of the 2019/8/9 event are shown in Fig 3.

Fig. 5
Fig.5 Event trace for frequency event on 2019/7/4, showing a smooth recovery path.The shaded area represents the deadband for Dynamic Containment[18]

Fig. 6 Fig. 7
Fig. 6 Imbalance and Frequency Response graph for 2019/7/4 event.PFR = 50% is assumed and Net Imbalance remains positive after the initial recovery above zero

Fig. 8 Fig. 9
Fig. 8 Imbalance and Frequency Response graph for 2019/5/31 event.Net Imbalance is briefly negative at t PFR , synchronous with the temporary drop in frequency recovery

Table 3
PFRStatic proportion of pPFR after t PFR % 50 threshold by Stage 5 as SFR increases.The initial path to recovery (Stage 5) is approximated, after which tertiary measures are reported to have intervened, correcting the imbalance.