A Spectral Model of Grid Frequency for Assessing the Impact of Inertia Response on Wind Turbine Dynamics

: The recent developments in renewable energy have led to a higher proportion of converter-connected power generation sources in the grid. Operating a high renewable energy penetration power system and ensuring the frequency stability could be challenging due to the reduced system inertia, which is usually provided by the conventional synchronous generators. Previous studies have shown the potential of wind turbines to provide an inertia response to the grid based on the measured rate of change of the grid frequency. This is achieved by controlling the kinetic energy extraction from the rotating parts by its converters. In this paper, we derive a spectral-based model of the grid frequency by analyzing historical measurements. The spectral model is then used to generate realistic, generic, and stochastic signals of the grid frequency for typical aero-elastic simulations of wind turbines. The spectral model enables the direct assessment of the additional impact of the inertia response control on wind turbines: the spectra of wind turbine output signals such as generator speed, tower base bending moment, and shaft torsional moment are calculated directly from the developed spectral model of the grid frequency and a commonly used spectral model of the turbulent wind. The calculation of output spectra is veriﬁed with non-linear time-domain simulations and spectral estimation. Based on this analysis, a notch ﬁlter is designed to signiﬁcantly alleviate the negative impact on wind turbine’s structural loads due to the inertia response with only a small reduction on the grid support.


Introduction
The power system needs to be balanced in real time to ensure a reliable supply of electricity. The imbalance between power production and consumption is eventually reflected as the grid frequency deviation from the reference value [1]. In response to climate change, many countries have stepped into the transition towards a more renewable energy-based power system [2]. Most renewable energy resources such as wind and solar are weather dependent, and renewable energy systems usually act as passive sources in the grid: they are connected to the grid through converters and do not interact with the grid frequency. Balancing the future electrical power system with a higher share of weather-dependent renewable generation is challenging: the renewable energy systems not only draw more uncertainties in the stochastic generation, but also reduce system inertia by replacing conventional synchronous generators (SGs) [3][4][5].
One advantage of modern wind turbines is that they can provide a synthetic inertia (SI) response to the grid, thus compensating part of the reduction in the system inertia due to high renewable energy penetration [6,7]. Generally, the synthetic inertia response of wind turbines can be achieved by controlling the kinetic energy extraction from the rotating parts through its converters, achieving a non-passive and grid frequency coupled power output. Typically, the amount of kinetic energy extraction can be determined by the measured rate of change of frequency (RoCoF) and can be achieved by controlling the generator torque [8,9]. Since the converter control time scale is much smaller than the inertial time scale, the fast inertia response from converter-based turbines can be considered analogous to the inherent inertia of synchronous generators. The inertia response control methods for the widely adopted modern converter-based wind turbines (Type-3 and Type-4 based on [10]) were summarized by [11].
Until today, the transmission system operator (TSO) responsible for Quebec, Canada, has been utilizing SI from wind turbines for years [12]. The recent grid code of Hydro-Québec [13] specifies the activation of SI response when a certain threshold of frequency deviation is met in under-frequency scenarios [14]. However, there is currently no TSO in Europe that has strictly required SI from wind power plants (WPPs). This is due to a large synchronous area of the Continental European (CE) grid where a large amount of SGs are still present. The European Network of Transmission System Operators for Electricity (ENTSO-E) suggests that a lowest allowable per unit system inertia at the synchronous area level needs to be considered in the future system design [15]. This will draw the necessity of providing SI from wind turbines if the allowable system inertia cannot be fulfilled by other sources.
Since the inertia response is achieved by controlling generator torque using RoCoF as the input, the stochastic frequency fluctuation needs to be studied first to understand the impact of inertia response on wind turbines. Several research papers have been carried out focusing on the frequency fluctuation characteristics and modeling for the duration above hourly power market trading [16][17][18][19]. The main conclusions from these studies are the non-Gaussian distribution of frequency deviation and the variation over the continents. Further, the long-term randomness is affected by consumer behaviors, fluctuations by weather-dependent renewable resources, and market rules [16,18]. On the other side, some authors directly studied the stochastic RoCoF behavior, which is directly linked to the control of the SI response. In [20], an algorithm is proposed to evaluate the probability of regional RoCoF, which can be used for grid operation and planning. In [21], the RoCoF probability based on measured data was analyzed, as well as the stress on wind turbine loads if SI response is required.
In [8,22], the historical measured grid frequency was applied to simulate the continuous provision of SI by wind turbines. However, measured grid frequency data are often not available. Further, traditional aero-elastic simulations of wind turbines use stochastic wind speed time series of usually 10 min based on defined spectral wind models to quantify the structural loads for certification [23]. Therefore, a method to extend this approach to quantify the impact of inertia response for wind turbines still needs to be developed. In this work, we derive a spectral-based model of the grid frequency by analyzing the historical measurement. The spectral model can then be used to generate realistic, generic, and stochastic signals of the grid frequency for typical aero-elastic simulations of wind turbines. Further, the model can be used in the frequency domain together with a wind speed spectrum to directly calculate the spectra of wind turbine output signals such as generator speed, tower base bending moment, and shaft torsional moment. This is helpful for tuning the SI response, for comparing its impact to the impact of the stochastic wind, for designing filters, and for investigating the impact of implementing SI control on turbine dynamics with or without filters.
The paper is structured as follows: Section 2 presents the data source and statistical analysis of the grid frequency. Section 3 provides the spectral estimation from the measurement data, develops the spectral model, and shows how to generate stochastic time series from the spectral model. In Section 4, the wind turbine model for assessing the inertia response impact on wind turbine dynamics is derived. A case study using a typical reference turbine model and controller, as well as the proposed spectral model are demonstrated in Section 5. Finally, Section 6 concludes this research and proposes possible future work.

Statistical Analysis of Grid Frequency
In this section, we describe the source and analyze the stationarity of the data used.

Data Source Description
The frequency data were measured at the campus of Flensburg University of Applied Sciences (synchronized with the CE grid), using a Dewetron 2010 measurement computer and the measurement software Dewesoft 7.1 [24]. The grid frequency was computed by measuring the grid voltage using a high-voltage module with a sampling frequency of 50 kHz. The computation of the grid frequency was done by fitting sine waves into measured time traces of the grid voltage. The frequency of the best fit sinusoidal wave was then recorded as the grid frequency. An overview of the measurement coverage is summarized in Table 1, where the valid days mean the number of days with complete measurement and the coverage means the number of valid days divided by the total number of days in a month. As for the excluded days, the measurement was not complete due to measuring system failure, and these days were not used for later analysis.

Stationarity Analysis
A random process such as the grid frequency fluctuation in this work can be considered as stationary if the statistical properties calculated over a time interval do not vary significantly from one interval to the next [25]. The spectral model should be derived on the time interval by which stationarity approximation can be met. If the sample time interval is too short, non-stationarity caused by fluctuations of low frequency can distort the spectral estimation results.
In terms of the atmospheric turbulence spectrum, the interval is chosen between 10 min and 30 min to meet the stationarity [26], and 10 min is specified by the IEC 61400-1 standard [23] to access the turbine load caused by turbulent wind flow. This is however not proper for grid frequency analysis, since the 15 min imbalance settlement period, the hourly clearing, and human daily routines can influence the power balancing and therefore the frequency fluctuation significantly. An analysis interval shorter than these marketand human behavior-related periods can make the data highly non-stationary. Similar to the study by [16], we used 24 h as the sample time interval, above which the data can be considered as stationary [19].
After ignoring the incomplete data period, for each day i, the meanf grid,i , variance σ 2 grid,i , and auto-correlation R grid,i (j) of the data are calculated by: where j is a discrete displacement and N = 527,346 is the total number of measurement data in each day given by a sampling period of ∆t = 0.1638 s. In the auto-correlation function R grid,i (j), the actual time lag is τ = j∆t. For simplicity and conformity in units, we often used the standard deviation σ grid,i , which is the square root of the variance. The daily mean and standard deviation of the grid frequency measured in 2018 are shown in Figure 1. For both the mean and standard deviation, there was no conclusive trend that changed following the days of the year, and they only varied around the yearly mean values, which were approximately 50 Hz and 0.02 Hz, respectively. This implies the stationarity of the data analyzed with a daily time interval [19].
For each day, the auto-correlation was also calculated, and the mean and variation are shown in Figure 2. The peak drawn by the auto-correlation function denotes the similarity between the shifted data and the original data. It can be seen that the cycles due to the hourly market clearing were pronounced, followed by the 30 min continuous intra-day trading [27]. Small peaks due to the 15 min imbalance settlement period were also clearly seen in the auto-correlation. From the confidence intervals, it can be concluded that the auto-correlation functions within the year were well centered due to the narrow 68% interval. Furthermore, the distribution of the auto-correlation was balanced because of the good agreement between the mean and median curves. Overall, the similarity of the daily auto-correlation over a full year justifies the daily analysis.

Derivation of the Spectral Model
In this section, the measurement data are used to estimate the power spectral density (PSD) of the grid frequency. A spectral model is then developed and fit to the estimated spectrum. The analytical PSD of RoCoF is also derived based on the spectral model of the grid frequency. Finally, the generation of stochastic grid frequency and RoCoF time series from the spectral model are discussed.

Estimation of the Grid Frequency Spectrum
The power spectral density can either be calculated via the Fourier transform of the auto-correlation function or estimated via the finite Fourier transform (FFT) [28] of the time domain records. The latter approach, more specifically Bartlett's method [29], was used in this work, that is: in which F {} and F * {} denote the FFT and the conjugate of the FFT, respectively. However, due to the non-continuous one year record, only the one-day-data blocks with effective continuous records were used. From each one-day-data block, the first 2 19 data points were processed to ensure an FFT length of 2 19 , which resulted in a faster calculation and avoided zero-padding. Therefore, only 2 19 × ∆t = 23.86 h for each day was used in the following analysis. This gave a discrete frequency vector f i range from ∆ f e = 1/(2 19 ∆t) Hz to half of the sampling frequency f M = 1/(2∆t) = 3.05 Hz, with a constant frequency step ∆ f e . Before applying the FFT, each data block was normalized to have unit variance. The mean estimated PSD was finally calculated by taking the sample mean of all the PSDs estimated from each data block.
The mean estimated PSD of the grid frequency is shown by Figure 3. There were obvious spectral peaks in frequencies that corresponded to 1 h, 30 min, and 15 min, which were in accordance with the auto-correlation analysis in the previous section. In addition, the peak at 1/60 Hz and also other spikes might be due to electrical or magnetic-related effects. The electrical-magnetic illusions from the grid frequency measurements were discussed in [30]. These effects could result in unnecessary inertia responses since they are most likely not related to power imbalances. Except for the spikes, there were some wave shapes in the frequency range higher than 0.1 Hz, which can be the result of local oscillations of SGs.

Modeling of the Grid Frequency Spectrum
Based on the estimation from the previous section, the PSD of the grid frequency generally showed a decreasing trend by two slopes in the logarithmic scale, as shown by Figure 3. The PSD model of grid frequency is developed and fit in this section, which is able to represent the general trend of this signal in the frequency domain.
A function that has a linear shape in the logarithmic scale can be written as: where d = 1 s is used to get a dimensionless frequency, or: where lg() = log 10 () denotes the common logarithm operation, lg(a) is the intercept at lg( f d) = 0 or f = 1 Hz, and b is a slope constant. However, the grid frequency PSD generally drops with two slopes, which is not possible to represent using only one constant b in the model. Thus, a function that combines two linear lines in different frequency ranges is introduced, that is: where b 1 and b 2 are constants determining the slopes and a 1 and a 2 are constants determining the intercepts. Here, a 1 ( f d) b 1 and a 2 ( f d) b 2 represent the two lines in the logarithmic scale, as shown in Figure 4. T( f ) is a transition function that allows S grid ( f ) to follow the two lines in different frequency ranges. It is expressed as: where c 1 and c 2 are constant frequencies determining the transition interval. With r = 10 lg(r) , the transition function has fixed values at c 1 and c 2 :  Table 2.
A general explanation to the expression is that T( f ) equals T(c 1 ) when f = c 1 , and it further approximates one for smaller f . T( f ) is equal to T(c 2 ) when f = c 2 and further approximates zero for larger f .
These characteristics allow the model S grid ( f ) to approximate two lines in different frequency ranges. By adjusting the constant frequencies c 1 and c 2 , the transition interval can be manipulated.
In order to fit the model with the PSD estimation results, the optimization problem is formulated as: The target is to minimize the error between the mean estimated PSD and PSD model on a logarithmic scale. The logarithmic scale is chosen to have same clarity in both low-and high-frequency ranges. The logarithmic scale error is divided by the discrete frequency f i to have equivalent weightings of the cost function. If 1/ f i is not considered, the optimization result tends to focus on better fitting on a high-frequency range due to the very dense discrete frequencies. The range of the discrete frequencies in optimization Equation (11) should be selected based on the specific application. For this work, a minimal frequency of f m was chosen to comply with the turbulence spectrum (10 min) used in the aero-elastic simulation, so f m then becomes the first discrete frequency larger than 1/(10 min) = 1/600 Hz in the frequency vector f i resulting from Section 3.1. The maximal frequency is just half of the sampling frequency f M discussed in Section 3.1, which is limited by the measuring devices.
Before applying the optimization, the spikes in the estimated PSD were removed because modeling these electromagnetic dynamics in the spectral model was not desired. The spike filtering can be accomplished by a simple median filter (done using "medfilt1" by MATLAB 2018b, MathWorks, Natick, MA, USA) [31]. Based on the discussions above, the fit parameters for the spectral model are summarized in Table 2. The estimated spectrum without the spikes and the fitted spectrum are shown in Figure 3, where the fitted model is able to represent the overall trend. Table 2. Summary of the fitted parameters for the spectral model.
Once the spectral model for the grid frequency is derived, it is not difficult to write out a spectral model for the RoCoF, that is: The reason for deriving the RoCoF spectrum is that one can also directly generate a time series of the RoCoF using the method to be discussed in the next section.

Generation of Stochastic Grid Frequency Time Series
The stochastic process theory has been applied to evaluate mechanics and structural loads for decades [32]. In wind turbine design, the spectral model is applied to simulate turbulent wind to assess the safety and reliability of the design [23]. With the requirement of a continuous inertia response from turbines, the stochastic frequency fluctuation can also be simulated to analyze the additional impact on wind turbines.
Based on [32], a common and simple way to simulate a time-domain signal from its PSD is to represent the signal by a Fourier series: where f j is the jth discrete frequency of the time series, K is the number of discrete frequencies,f grid is the mean grid frequency of the simulated time interval, φ j is the uniformly distributed random phase between zero to 2π, and A j is the amplitude determined by: with ∆ f g the frequency step for the generated time series, σ grid the target standard deviation of the generated time series, and σ model the standard deviation determined by the model: With the model of the PSD, the time series can be efficiently calculated using the fast inverse Fourier transform [28]. The discrete frequency for generating the time series ranges from ∆ f g to f K , with a step length of ∆ f g . It is worth mentioning that f K = 1/(2∆t g ), where ∆t g is the sampling time step of the generated time series. For the simulation in Section 5, f K was chosen to be 2 14 /600/2 = 13.65 Hz to ensure a good PSD estimation from the simulation results. However, the spectral model was not fit in the frequency range higher than 3.05 Hz; we still assumed the model could describe the PSD in this range. The assumption however did not influence the assessment for this study, which will further be discussed later in Section 5. Figure 5 compares the generated time series of f grid and RoCoF with one set of 10 min measurements. The time series were generated using the developed spectral model with random phase angles and scaled to have identical variance with the measured time series.

Assessment of the Inertia Response Impact
Based on the previous sections, the stochastic grid frequency time series can be generated and used to assess the inertia response impact on wind turbines through time-domain simulations. Here, we propose a frequency-domain method, which allows assessing the impact of the synthetic inertia response more directly. For this method, we first derived a simplified low-order wind turbine model and a baseline wind turbine controller extended with synthetic inertia control. The models were then combined and linearized. Finally, the wind turbine output spectra were derived.

Wind Turbine Model
The reduced wind turbine model [33], which represents the overall dynamic and the main motions, can be divided into the aerodynamics, drive-train dynamics, and tower fore-aft dynamics. Figure 6 shows a sketch of the reduced wind turbine model. The aerodynamics are modeled by: where M a and F a are the aerodynamic torque and thrust force, respectively, ρ is the air density, R is the rotor radius, Ω r is the rotor side rotational speed, and c P and c T are the power and thrust coefficients determined by the tip speed ratio λ and pitch angle θ. The relative rotor equivalent wind speed v rel is simplify calculated by superposition of the rotor effective wind speed v 0 (REWS) and the tower top speedẋ T .

Drive-Train Dynamics
The three degree-of-freedom (DOFs) drive-train model used here was based on a two coupled mass-spring-damper systems [34,35] modeling the dynamics of the rotor speed Ω r , the generator speed Ω g , and the shaft torsional angle φ T : Here, i gb is the gear-box ratio, J r and J g are the rotor-side and generator-side equivalent inertia, M g refers to the generator torque, and M r is the rotor-side shaft torsional moment.
Further, k eS and c eS are the equivalent drive-train torsional stiffness and damping constant, respectively.

Tower Fore-Aft Dynamics
The tower top fore-aft dynamics were modeled by a mass-spring-damper system with two dynamic states: Here, m eT , k eT , and c eT are the tower equivalent modal mass, damping, and stiffness, respectively. Further, x 0T is the static tower top position due to the gravity of the overhanging rotor mass. With the turbine hub height z H , the tower base fore-aft bending moment is:

Controller Model
A typical variable-speed wind turbine is controlled by a blade pitch and generator torque controller. The generator torque controller is here extended by the synthetic inertia. A baseline collective blade pitch control is achieved by a proportional integral (PI) controller [36]: where Ω g,ref is the speed control reference, δ PI is the integral of the error from the PI controller, k p the proportional gain, and T I the integrator time constant. The pitch controller is only active above rated wind speed, and k p and T I are scheduled to have a constant closed-loop behavior.
The main idea of SI from wind turbines is to emulate the inherent inertia response naturally existing in a synchronous generator [6]. The inertia power is calculated by: where H is the target inertia response constant to emulate, P rated is the turbine rated power, and f grid is the grid frequency. The baseline generator torque controller uses the generator speed and blade pitch as the input and returns a torque reference value, and it can simply be expressed as M gc (Ω g , θ). The details of the control algorithm can be found in [36]. The total generator torque control signal with inertia response included then is: with η g the generator efficiency.

Linearized Closed-Loop Model
For each operating point, the models listed above can be linearized and organized in the following state-space representation [37]: where x is the vector with n state variables, u is the vector with p input variables, and y is the vector with q output variables. A is the state matrix containing n × n time-invariant factors, B is the n × p input matrix, C is the q × n output matrix, and D is a q × p feedthrough matrix. For the wind turbine model we used in this work, the vectors are: Since the grid frequency deviation from the nominal value is relatively small, f grid was assumed to be the constant nominal frequency when deriving the linearized model. The extended expressions of the matrix A, B, C, and D are given in Appendix A.
In the frequency domain, the matrix of transfer functions can be calculated by: where I is an identity matrix and s denotes the complex frequency of the Laplace transformation.

Filter
An adjustable notch filter for the grid frequency was used in this work, with the following transfer function: where f notch is the notch frequency in Hz, b notch is the bandwidth in Hz, and d notch is the depth of the notch filter at the notch frequency. The benefits of using the notch filter for SI control are discussed in Section 5.

Wind Turbine Output Spectra
For a single turbine, we can assume that the stochastic rotor effective wind at the local site is uncorrelated with the stochastic grid frequency fluctuation, which is caused by the overall power imbalance. Furthermore, the effect of the inertia response from a single turbine should not change the overall grid frequency. Thus, the calculation of the PSD of the output variables can be simplified as a special case of mutually uncorrelated inputs [25].
The PSD of output variables, in our case the generator speed, tower base fore-aft bending moment and shaft torsional moment can thus be calculated with: In Equations (37)-(39), G ij ( f ) is the element in the transfer function matrix G(s) from input i to output j. Further, S v 0 is the PSD of the rotor effective wind speed v 0 , which can be derived from a spectral model of the turbulence. In this case, we used the IEC Kaimal model [23]. The details regarding the REWS spectra were discussed in [33].
If a filter is applied to the RoCoF, the output PSD of rotor speed becomes: The only difference is additionally multiplying the input PSD by the square of the transfer function |G notch ( f )| 2 . Similarly, the output PSD for other variables can be obtained by multiplying the input PSD by |G notch ( f )| 2 in Equations (38) and (39).

Case Study
In this section, we chose an above-rated operating point and used the methods presented in previous sections to assess the impact of SI response on wind turbine dynamics. The turbine model was established based on the 5MW reference turbine [36] and linearized at 20 m/s. An invariant inertia response constant H = 12 s was considered. The spectral model of the RoCoF derived in Section 3.2 together with the spectral model of the REWS introduced by [33] were used. The impact of the SI was assessed in the frequency domain and the time domain. The results were then compared and discussed.

Frequency-Domain Analysis
In a first analysis, the linearized model and controller were used to calculate the transfer function from Equation (35). In Figure 7, the magnitude of the transfer functions from every single input to every single output is displayed. The peak at 0.32 Hz corresponds to the natural frequency of the tower bending, and the peak at 2.22 Hz is the natural frequency of the shaft torsion. The simulations in Section 5.2 were carried out with a time step of ∆t sim = 600/2 14 s ≈ 0.0366 s to avoid zero-padding in the PSD estimation in Section 5.3 and to be between five and 10 times faster then the fastest dynamics, leading to a Nyquist frequency of 1/(2∆t sim ) ≈ 13.65 Hz, outside of the frequency domain used for the fitting of the spectral model in Section 3. The frequencies higher than 3.05 Hz (maximum frequency used to fit the spectral model) were well damped. Therefore, frequencies higher than 3.05 Hz and thus the limitations of the fitting process would not have large impact on the PSD of the output variables. In the next step, the PSD of the RoCoF was calculated based on the spectral model derived in Section 3.2. Based on the study by [21], the standard deviation of RoCoF observed from CE grid ranges from 0.0035 Hz/s to 0.0055 Hz/s, depending on the renewable energy penetration. To emphasize the dynamic impact of SI, the PSD of the RoCoF was scaled in this case study to have a standard deviation of 0.04 Hz/s on a 10 min scale, thus around 10 times higher than the usual value, indicating a grid with very high renewable energy penetration and low system inertia. The resulting spectrum is displayed together with the spectrum of the rotor-effective wind speed v 0 in the top row of Figure 8.
Then, the PSDs of the output variables were calculated based on Equations (37)-(40) and shown in the left and middle column of Figure 8. The PSDs of the output variables are separated to clearly see the independent contribution of each single input variable. Compared with the contribution from the turbulent wind, the contribution of the RoCoF to the PSD of rotational speed and tower fore-aft bending moment was significantly smaller. The RoCoF mainly impacted the PSD of the shaft torsional moment near the natural frequency of the shaft. Therefore, a notch filter was designed in the following steps to filter the RoCoF around the shaft torsion natural frequency based on Equation (36). The notch filter was designed to have a damping effect down to 10% at the notch frequency f notch = 2.22 Hz, which resulted roughly in a reduction of the spectrum of the generator speed and shaft torsional moment to 1% at the notch frequency. The width b notch was chosen to be 0.6 Hz to also reduce the spectrum to a similar level at the surrounding frequencies, while keeping the time delay introduced by the filtering small. Figure 9 shows the magnitude, phase, and time delay characteristics of the designed notch filter.
Finally, the notch filter was applied following Equation (40), and the effects are shown in the right column of Figure 8: the PSD of the RoCoF was filtered out at the notch frequency range, and the resulting spectral peaks for generator speed and shaft torsional moment were alleviated significantly down to roughly the expected 1%.
The peaks at the tower's natural frequency caused by the ReCoF were much smaller than the ones caused by the turbulent wind. Therefore, no additional notch filter was used for the tower natural frequency to avoid further reducing the inertia response, since the inertia power was proportional to the RoCoF based on Equation (28).

Time-Domain Analysis
As for the time-domain simulation, the stochastic input time series were generated based on the approach shown in Section 3.3. The two input time series REWS and RoCoF were generated as independent and thus uncorrelated variables. Then, the wind turbine was simulated (nonlinear ODEs were solved by "ode4" (Runge-Kutta fourth order) by MATLAB 2018b, MathWorks, Natick, MA, USA.) over 600 s using the input time series and the non-linear turbine model and nonlinear controller from Sections 4.1 and 4.2, respectively.
The following three simulations were carried out: a baseline case where no SI response was considered and two further simulations where the SI response was included without and with filtering of the RoCoF signal. The first 30 s are shown in Figure 10. As predicted by the frequency analysis, the impacts of SI on rotor speed and tower fore-aft bending moment were neglectable, and the output time series were almost identical and thus overlapping in the plots. The main impact of the SI without the filter was the additional variation in the rotor side shaft torsional moment close to the natural frequency of 2.2 Hz, as indicated by the orange line. Applying the notch filter damped this variation, and the resulting shaft torsional moment (blue line) was close to the one of the baseline case (red line). However, applying the notch filter also reduced the inertia response power P SI from the turbine.
To further investigate the dynamic impact due to SI response, we defined 16 different initial random seeds and performed 16 time-domain simulations with different randomly generated input time series, following the Monte Carlo method. For each seed, the time series of the outputs variables from the simulations were collected, and the corresponding standard deviations and damage-equivalent loads (DEL) were calculated, using a Wöhler exponent of four, 2 × 10 6 as a reference number of cycles, and a lifetime of 20 years. The inertia response energy E SI was obtained by integrating the absolute value of inertia power P SI with respect to time. The sample mean values from the different initial seeds were calculated and summarized in Table 3. The statistic further illustrated that SI mainly influenced the shaft loads. Applying the notch filter generally led to a compromise between decreasing loads and reducing the inertia power. The relative change in the shaft DELs compared to the no-SI scenario was reduced from 236.56% to 48.53% by applying the notch filter, while the inertia energy provided only reduced by 2.85%. However, a more detailed study including dynamic grid simulation and an aero-elastic model of the wind turbine was necessary to finally evaluate the impact of the notch filter.

Comparison of the Results and Discussion
The PSDs from the frequency analysis were then compared to the PSDs estimated from time-domain simulations; see Figure 11. Here, the PSDs from time series were estimated by taking the sample mean of 16 realizations by different random seeds. The good agreement validated the approach we applied. The plots again proved that the additional impact on shaft load due to SI inertia response could be alleviated by introducing a notch filter. The frequency-domain model for the linear system could be applied together with the spectral model of input variables to effectively calculate the spectra of outputs variables. It was evident to directly use the model-based approach to analyze the spectra of output variables without running additional time-domain simulations. The plots also show that the linearization of the used turbine model at the chosen operating point provided good approximations to non-linear simulations.

Conclusions
This paper was motivated by the current trend and demand for inertia response from wind turbines. The paper first summarized the statistical properties of the European grid frequency measured over one year. Based on this analysis, a spectral model for the grid frequency and its derivative were introduced, which were capable of describing the general characteristic of the grid frequency spectrum. In a case study, this spectral model was used together with a linearized wind turbine and controller model to efficiently assess the impact of synthetic inertia response on wind turbine dynamics. Based on this analysis, the synthetic inertia response was extended to significantly reduce the impact on the structural loads. Thus, the main contributions of this work were: (1) statistical analysis of the European grid frequency, (2) derivation of a spectral model for the grid frequency, (3) assessment of the impact of synthetic inertia response on wind turbines dynamics, (4) and extension of the synthetic inertia response to reduce the main impact on wind turbine loads. The main conclusions for all four topics are summarized in the following paragraphs.
The statistical analysis suggested that the daily grid frequency measurement can be considered stationary. Extra peaks were observed in the auto-correlation function for time lags corresponding to the hourly market clearing, 30 min continuous intraday trading, and the 15 min imbalance settlement period.
The spectral model for the grid frequency was based on the observation that the spectrum of the grid frequency generally had a descending trend with two main slopes. This shape was represented by the developed spectral model, which had six parameters. The parameters were obtained by a fitting process in this work, but can also be extracted visually from the spectrum. A spectral model of the rate of change of frequency (RoCoF) was derived by simple calculations. With the spectral model, stochastic time-domain signals of the grid frequency and for the RoCoF could be generated, scaled to specific conditions, and integrated into simulations of wind turbines including aero-elastic tools.
For the assessment of the impact of synthetic inertia response on wind turbine dynamics, a reduced-order wind turbine model was presented. The turbine model can be combined with the spectral models of wind speed and RoCoF to derive the spectra of load-related signals such as the generator speed, shaft torsion moment, and tower fore-aft bending moment. It was verified that this model can produce accurate spectra directly without the need for running simulations. The case study revealed that even for very high grid fluctuations and a strong synthetic inertia response, the impact of the RoCoF on the generator speed and tower base bending moment was very small in comparison with the impact of the wind fluctuations. However, a significant increase was observed for the shaft torsional moment at the natural frequency. The results were confirmed by simulations where the additional shaft loads were increased by over 200%, while the extra impact on generator speed variations and tower load were both below 0.5%.
As an extension of the synthetic inertia response, a notch filter was proposed to reduce the impact on the shaft loads without a large impact on the synthetic inertia response. For the case study, the additional shaft loads could be considerably reduced down to around 50% in comparison to the no-filtering scenario, while roughly only 3% of the synthetic inertia response was sacrificed. This would make the approach much more likely to be applied to real wind turbines.
This work represents a first study in this field and will benefit from further work. Due to the limitation of the grid frequency measuring device, the proposed spectral model was only able to reflect the spectrum in a specific frequency range in Central Europe. A comparison with grid frequency data from various synchronized areas will help to further parameterize our model. Additionally, a multi-resolution analysis using the fast wavelet transform could replace the FFT to capture more detailed patterns of the grid frequency. Furthermore, the method to assess the impact on wind turbine dynamics focused on above-rated wind speed conditions where reduced linearized models provided good approximations. For assessing the inertia response in more detail, the method should be validated against aero-elastic simulations. Lastly, further studies with dynamic grid models should be carried out in the future to find a good compromise for the notch filter parameters to balance load reduction and grid support.

Abbreviations
The following abbreviations are used in this paper:

Appendix A
The state space matrices A, B, and C of the linearized model can be calculated by the partial derivatives as shown by: D is a 3 × 2 matrix with all elements being zeros.