Inertia Estimation Through Covariance Matrix

This work presents a technique to estimate on-line the inertia of a power system based on ambient measurements. The proposed technique utilizes the covariance matrix of these measurements and solves an optimization problem that fits such measurements to the synchronous machine classical model. We show that the proposed technique is adequate to accurately estimate the actual inertia of synchronous machines and also the virtual inertia provided by the controllers of converter-interfaced generators that emulate the behavior of synchronous machines. We also show that the proposed approach is able to estimate the equivalent damping of the classical synchronous machine model. This feature is exploited to estimate the droop of grid-following converters, which has a similar effect of the swing equation equivalent damping. The technique is comprehensively tested on a modified version of the IEEE 39-bus system as well as on a dynamic 1479-bus model of the all-island Irish transmission system.


A. Motivation
The inertia constant has become a volatile parameter in power systems with high penetration of renewable and nonsynchronous resources [1].Therefore, the estimation of the available inertia is a valuable information that can help system operators to ensure the security of the grid [2].This paper addresses this topic and aims at developing an on-line inertia estimation method based on ambient measurements, their stochastic behavior and a fair knowledge of the grid model.

B. Literature Review
There are mainly two methods for the on-line estimation of the inertia available in a power system: (i) methods based on the measurement of frequency and power variations after a disturbance or by injecting a probing signal and (ii) methods that utilize ambient measurements [3], [4].Both methods have advantages and shortcomings.F. Bizzarri is with Politecnico di Milano, DEIB, p.zza Leonardo da Vinci, no.32, 20133 Milano, Italy and also with the Advanced Research Center on Electronic Systems for Information and Communication Technologies E. De Castro (ARCES), University of Bologna, 41026 Bologna, Italy.(e-mail: federico.bizzarri@polimi.it).D. del Giudice, S. Grillo, D. Linaro, and A. Brambilla are with Politecnico di Milano, DEIB, p.zza Leonardo da Vinci, no.32, Milano, 20133, Italy.(e-mails: {davide.delgiudice,samuele.grillo,daniele.linaro,angelo.brambilla}@polimi.it).
Italian MIUR project PRIN 2017K4JZEE 006 funded the work of S. Grillo (partially) and D. del Giudice (totally).
The methods that are based on disturbances can be very accurate [5] provided that the on-set of the disturbance is correctly detected [6]- [9].This can be challenging in low-inertia power systems, which show a "rich" dynamical behavior [10].Moreover, the estimation cannot be performed continuously, but only following disturbances [11] [12].Methods based on signal probing must inject active power of adequate magnitude and require the installation of specific equipment or the modification of pre-existing controllers [13].
The methods that utilize ambient measurements are conceptually more challenging than those based on disturbances.The main advantage of these methods is that they employ measurements obtained in normal operating conditions and hence the estimation can be performed in a continuous fashion.The starting point is to build the dynamical model of the system, which can be either complete or reduced to a small number of coherent areas [14], [15].Then, the values of the parameters and of the observable state variables are fit to available measurements.It is worth mentioning that the number of areas and, consequently, the number of equivalent generators, whose inertia is to be estimated, is relevant.In [16], the authors raise the question whether the center of inertia (COI) of a single-machine-equivalent power system model could be identified through ambient measurements.The answer given is that the combined effects of a well-damped COI dominant mode and of a low-pass-filter action of the Ornstein-Uhlenbeck (OU) processes modelling load behavior make it practically impossible to perform this task.However, this conclusion does not hold when a multi-machine system is considered.
In the latter group of inertia estimation methods, we cite, for example, [17], where the authors obtain the fitting based on a robust Kalman filter.In [13] such a task is performed using a closed-loop identification method.In [18]- [20] inertia estimation is performed by observing the step response of the reduced-order model obtained from the full state-space model describing the grid, which is, in turn, obtained through parameter identification using ambient measurements.In [21] the power grid is reduced to a two-machine system and then the system inertia is estimated using inter-area oscillation modes.In [22], the swing equation of the area single-machineequivalent generator is transformed in Fourier domain.Then, the electromechanical oscillation, obtained by applying a fast Fourier transform to the generator estimated electrical power and rotor speed, are used to estimate the equivalent inertia.This method requires capturing the whole electromechanical oscillation trajectory.The approach proposed in [23] estimates the inertia constants of each generator by applying the regres-sion theorem to a set of OU processes.The method utilizes the covariance and correlation matrices of the state variables to estimate the state matrix of the system.Then, the inertia constants are obtained solving a least-square problem.The main limitation of [23] is the assumption that the relationship between the active power variations of synchronous machines and the voltage angle variations, as well as the active and reactive power variations of converter-interfaced generators are linear.In a similar way, in [24], the regression theorem of a multivariate OU process is applied to power systems including converter-interfaced generators (CIGs), modeled in detail with their phase locked loops (PLLs) and supposed to work at unity power factor.The estimation is carried out separately for CIGs and synchronous generators by identifying the relevant submatrices of the system state matrix.

C. Contribution
The approach proposed in this work falls in the group of methods that use ambient measurements.We propose to exploit the colored noise due to random fluctuations of the power consumption of the loads and their impact on bus voltages, line currents, and power flows across different areas of a power system.Then we define analytically the statistical properties of the colored noise of these electrical quantities as functions of the system parameters and of the inertia present in the system.Finally, we utilize the variance of the electrical measurements to minimize a non-linear least-squares cost function with respect to the equivalent inertia constants as well as the equivalent damping coefficients of the generators, either synchronous (see, e.g., equation (3.203) in [25]) or powerelectronic based, that are connected to the grid.The solution is the sought value of the inertia constant or of the damping coefficient that best fits the measurements and that satisfies the grid constraints.
The main benefit of the proposed approach is that it does not require a contingency or a large disturbance to occur to be able to estimate the inertia (or the damping).Only ambient measurements are needed, together with a fair knowledge of the model of the grid.The price for this is the need to take measurements in a given time period.Nevertheless, as shown in the case studies, the proposed approach is accurate and provides estimations of the inertia in the same time scale of short-term dispatch and adjustments markets.The proposed method is original, as it shows for the first time how to perform inertia estimation by resorting to the covariance matrix of the power system without needing to estimate the state of the network or to derive (or measure) the rotor speed or any other internal variable of both synchronous machines and other devices performing frequency control.

D. Organization
The remainder of the paper is structured as follows.Section II gives the theoretical mathematical background of power systems' modeling and in Section III the mathematical model is enriched with the introduction of noise injections to the power system, being this stage instrumental to model load variability and perform the estimation.Section IV presents the proposed estimation procedure.Then the proposed method is validated through numerical simulations, which are discussed in Section V.This section also discusses the challenges posed by the proposed approach and presents a solution based on a proper filtering of the measurements.Finally, conclusions are drawn in Section VI.

II. POWER SYSTEM MODEL
The proposed estimation procedure assumes that the generators connected to the grid are synchronous machines and fits the measurements to this model.In this section, we formulate the equations of the system according to this assumption.It is important to note, however, that every simulation is carried out utilizing a full-fledged model representing in detail the dynamical behavior of all generators, either synchronous or non-synchronous, and their controllers.
To describe the dynamics of the power system model (PSM), we consider the following set of differential algebraic equations (DAEs) where, assuming that M is the number of machines, the meaning of symbols in (1) is as follows: -Ω: the base synchronous frequency in rad/s; -ω(t) ∈ R M : the per-unit rotor speeds of the machines; -ω 0 ∈ R: the per-unit reference synchronous frequency; -δ(t) ∈ R M : the rotor angles of the machines; -M ∈ R M×M : a diagonal matrix whose entries model the inertia constants H of the machines.In particular, M jj = 2H jj (for j = 1, . . ., M ); -x(t) ∈ R N : the N state variables of the PSM (ω and δ excluded) that can influence the dynamics of the machines, and T ∈ R N ×N is a mass matrix; -y(t) ∈ R S : all the algebraic variables of the PSM but v and θ; -D ∈ R M×M : a diagonal matrix whose entries d jj (for j = 1, . . ., M ) model the damping factor of the machines; -v(t) ∈ R P and θ(t) ∈ R P : bus voltages and phases, respectively, where P is the number of buses; -P m (ω, x, y) ∈ R M : the mechanical power regulated by controllers depending on ω, x, and y; -P e (δ, x, v, θ, y) ∈ R M : the electrical power exchanged by machines; -F ( • ) accounts for regulators and other dynamics included in the system; and -G( • ) models algebraic constraints such as lumped models of transmission lines, transformers and static loads.Equation ( 1) can be conveniently rewritten as where and Λ is a nonsingular diagonal matrix including M .G and F are assumed to be continuously differentiable in their definition domain and matrices of their partial derivatives are referred to as If the conditions of the implicit function theorem are satisfied, (2) can be rewritten as the equilibrium points of which, say ξ o , satisfy the condition with the Jacobian matrix

III. INCLUSION OF NOISE
We assume stochastic noise is injected into the PSM as where η is a vector of Z independent Ornstein-Uhlenbeck (OU) processes [26], and Ξ ∈ R (2P +S)×Z is a constant matrix.
The OU processes are defined through a set of stochastic differential equations (SDEs), as follows where the drift Υ ∈ R Z×Z and diffusion Σ ∈ R Z×Z are diagonal matrices with positive entries, W t ∈ R Z is a vector of Z Wiener processes, and the differentials rather than time derivatives are utilized to account for the idiosyncrasies of the Wiener processes.The OU processes are characterized by a mean-reversion property and show bound standard deviation.Moreover, these processes show a spectrum that is an accurate model of the stochastic variability of power loads [27]- [31].
In (6), we assume that noise injection models the effect of the consumption randomness of Z loads.Note that, in (6), the vector of stochastic processes η appears only in the algebraic equations.This is assumed for simplicity but without lack of generality.The interested reader can find for instance in [29] a stochastic PSM model where the noise perturbs both differential and algebraic equations.
We also assume that the action of η can be safely modeled as a small-signal perturbation around an equilibrium point.Hence, the solution of ( 6) can be written as that is, the effects of perturbations are assumed as additive.It is thus possible to linearise (6) around ξ o and obtain the random ordinary differential equation [32] where the F ξ , F ζ , and The small-signal variations of the algebraic variables are given by Augmenting ( 8) with (7) allows obtaining the linear SDEs (in narrow sense [26]) governing the overall dynamics of the linearised noisy PSM.Adopting the typical formalism of the SDEs, the complete set of linearised SDEs reads The solution of (10) with a normally distributed (or constant) initial condition is a (2M + N + Z)−dimensional Gaussian stochastic process.

IV. PROPOSED TECHNIQUE TO ESTIMATE THE INERTIA
The mean and the covariance matrix of process ( 10) can be derived by solving two sets of linear ordinary differential equations (ODEs) [26].Since (10) is stable under the hypothesis that ( 3) is stable at ξ o , these ODEs reveal that, at steady state, the mean of X t is zero, and its K Xt covariance matrix derives from the solution of the following Lyapunov equation The diagonal elements of K Xt are the (steady-state) variances of the components of the X t process.In particular, the last Z diagonal elements, corresponding to the η sub-vector of X t , are the (steady-state) variances of the Z independent OU processes introduced in (7).Hence, these terms can be written as σ 2 z /(2υ z ), for z = 1, ..., Z, where σ z and υ z are the diagonal elements of Σ and Υ, respectively [26].The remaining 2M + N diagonal elements of K Xt refer to the ξ η sub-vector of X t and are influenced by η through the submatrix (10).According to (9), ζ η can be written as a linear combination of the entries of X t through E, thus being a multidimensional Gaussian process, too.Hence, it is possible to write the K ζη covariance matrix of the small-signal algebraic variables as the quadratic form [33], [34] The proposed approach exploits (12) and the fact that K ζη (and, thus, the variance of the algebraic variables) depends through K Xt and A on the elements of Λ, a subset of which are the inertia constants of the synchronous machines.
Let us assume one wants to estimate a subset M of these inertia constants and let Z be a set of PMU measurements of bus voltages and line currents flowing through a given number of transmission lines and/or transformers.
The above derivations show that, if the model and the parameters of the grid are known, one can compute the σ 2 Z variances of the measured quantities, which are the diagonal elements of the K ζη matrix.However, it is also possible to solve an inverse problem: knowing through measurements the variance of the elements of Z, the values of the elements of M can be determined by minimizing the following non-linear least-squares cost function (see also Fig. 1) We note that the procedure above is agnostic w.r.t. to the devices that are connected to the grid.That is, the estimation of the inertia is obtained by fitting the synchronous machine model to the measurements.If there exist non-synchronous devices, such as converter-interfaced generation, their "equivalent" inertia constants can be obtained using (13) regardless of the fact that the converter is set up to be grid forming or grid following.The procedure is also agnostic in terms of the parameters to be estimated, provided that such parameters do have an effect on the frequency variations of the system.In particular, we use this property to estimate also the equivalent damping coefficients of the devices, which can be attained by minimizing a cost function for a subset of damping coefficients D, in a similar way as for M in (13),

V. CASE STUDIES
In this Section, we consider two grids, the well known IEEE 39-bus system and a 1479-bus dynamic model of the all-island Irish transmission system AIITS.The IEEE 39bus system allows us illustrating the various features and challenges of the proposed method.With this aim we first consider the conventional system with synchronous machines and their regulators (Section V-A1).Then, we evaluate the behavior of the proposed method when the system includes grid-forming (Section V-A2) and grid-following converters (Section V-A3).The all-island Irish system, on the other hand, serves to demonstrate the robustness of the proposed method when applied to a real-world complex network.

A. IEEE 39-bus System
We use as a benchmark the IEEE 39-bus power system shown in Fig. 2 [35].This is a simplified model of the transmission system in the New England region in the Northeastern United States and is composed of 10 generators, 34 lines, 19 loads, and 12 transformers.The network models and parameters are those adopted and directly extracted from its DIgSILENT PowerFactory implementation [36] and derived from [37].The G 1 generator, which represents the connection of the New England System to the rest of the US and Canadian grid, is modeled with a constant excitation [37], viz.there is no automatic voltage regulator (AVR) connected to it.On the contrary, the other generators are equipped with AVRs (specifically, rotating excitation systems of IEEE Type 1 according to [37]).Governors are considered as IEEE Type G1 (steam turbine) for G 2 -G 9 , and IEEE Type G3 (hydro turbine) for G 10 .Finally, in Table I we report the value of inertia constants and power rating of the 10 generators.
The active power consumption of each of the Z = 19 loads is altered by one of the η z (t) OU processes (z = 1, . . ., Z) in (10), which has zero mean, and variance σ 2 z /(2υ z ).All loads are assumed to incorporate stochastic power fluctuations [29]  where (for z = 1, . . ., Z) P L0 z is the nominal active power of the load, V 0z is the load voltage rating, V z (t) is the bus voltage at which the load is connected and γ governs the dependence of the load on bus voltage.In the simulation below, we assume υ z = 0.5 and σ z is defined in such a way that the standard deviation of η z (t) is 5% of P L0 .The zero mean implies that stochastic load power fluctuations do not perturb, on average, the operating point of the system.Furthermore, we assumed γ = 0.
To carry out the simulations discussed below, the numerical integration of the multi-dimensional OU process in (6) was based on the numerical scheme proposed by Gillespie in [38].Furthermore, the second-order trapezoidal implicit weak scheme for stochastic differential equations with colored noise, available in the simulator PAN [39]- [41], was adopted [42].

1) Conventional Power System:
The objective of the first scenario considered in this case study is to estimate the inertia constants of the 10 generators of the IEEE 39-bus system across a time period spanning one day.The target is to obtain an estimation of the inertia constants every 15 min (estimation window) of the working day.To this end, we compute the variance of the currents flowing through the generators, using a 15 min-long moving window.The values of variance are the elements of the M set introduced in Sec.IV.The sampling period is 25 ms, i.e., a sampling rate of 40 Hz. 2e also assume that the matrices F ξ , F ζ , G ξ , G ζ are constant, i.e., that the parameters of the IEEE 39-bus system do not vary during the simulation.This assumption simplifies the analysis but is not a binding requirement of the proposed procedure.If the operating point and/or the topology of the grid change, one just needs to update those matrices.
The H inertia constant estimates were obtained with the MATLAB Global Optimization Toolbox and the fgoalattain function.The search interval of the optimization procedure was lower-bounded to zero since inertia constants cannot be negative.We performed 20 independent trials each one being a 24 h simulation of the grid and an optimization problem is solved every 15 min of simulated time, thus collecting 1920 estimates per synchronous generator.The violin plots3 in Fig. 3 (one for each generator) summarize the obtained results.The median of the estimated value of the inertia constants (black solid dots in Fig. 3) are in good agreement with the corresponding nominal values (empty squared markers in Fig. 3) listed in Table I.Despite the amplitude of the interquartile range (IQR) intervals (magenta vertical bars in Fig. 3) are reasonably narrow, the upper and lower adjacent values (green solid dots in Fig. 3) are quite far from the median Fig. 3. Violin plots of the inertia estimations of each generator of the IEEE 39-bus system.The black solid circle markers correspond to the median value of the optimization results, whereas the empty square markers represent the nominal values of the inertia constants.The magenta bars represent the IQR, viz. the spread difference between the 75 th and 25 th percentiles of the data.The green solid circle markers represent the upper adjacent value (i.e., the largest observation that is less than or equal to the third quartile plus 1.5×IQR) and the lower adjacent value (i.e., the smallest observation that is greater than or equal to the first quartile minus 1.5 × IQR).Time window: 15 min.Optimization performed on the unfiltered currents.
values, and the same holds for the tails of the distribution of the estimated values.The maximum value of the relative percentage error, computed as , is obtained for G 5 and amounts to 3.87%.The 0.44% minimum value is achieved for G 3 .
The quality of the results can be improved by filtering the considered currents before computing their moving variance.This is done by resorting to a steep bandpass filter acting from 0.1 Hz to 1.5 Hz.This bandwidth was chosen to remove the low frequency contribution of the OU process.This is done since these components account for the very slow fluctuations of the currents that heavily affect the trend of their moving variance w.r.t. the steady-state values of the variance itself.The effect of this filter is shown in Fig. 4 for the i 8 current (i.e., the current flowing through G 8 ) in terms of the absolute value of the relative percentage error between the moving variance and its steady-state value over 24 h in one of the 20 trials.In particular, the black and red curves refer to the unfiltered and the filtered case, respectively, and the improvement over the unfiltered case is evident.
The state equations governing the dynamics of the filters were added to the set of ODEs reported in (3) and this allowed to derive the steady-state variance of the filters output solving (11).The inertia constants of the IEEE 39-bus generators was estimated from the filtered currents and the results are reported in Fig. 5.The effect of this signal processing is a significant reduction of the amplitude of the interval spanned by the upper and lower adjacent values, thus guaranteeing a more reliable estimate of the inertia constants.The ǫ H % relative percentage error spans from 0.14% (for G 2 ) to 1.78% (for G 9 ).These results illustrate well that the proposed approach accurately estimates the inertia constants of the ten synchronous generators of the system.Fig. 4. The E % percentage relative error of the moving variance of the i 8 current, calculated over a sliding window of 15 min, with respect to the steady-state variance of i 8 .The black and the red curves refer to the unfiltered and filtered case, respectively.Fig. 5. Violin plots of the inertia estimations of each generator of the IEEE 39-bus system.The optimization was performed on the filtered currents, on a time window of 15 min.Note the different scale of the vertical axes w.r.t.Fig. 3.

2) Inclusion of Grid-Forming Converters:
This section illustrates the behavior of the proposed estimation technique for a scenario where the IEEE 39-bus system is modified to include a variable generation/load profile and grid-forming converters.These devices are assumed to mimic the behavior of synchronous machines and are thus expected to provide a virtual inertia to the system.The controllers of grid-forming converters are based on [43].The objective of this section is twofold: (i) showing that the proposed inertia estimation approach works correctly under variable generation/load conditions and (ii) illustrating the ability of the proposed estimation technique to capture the effect of grid-forming converters.
For what regards point (i), we first modified the IEEE 39bus system by mimicking a typical daily variation of the loads.To do so, we multiplied the power absorbed by the loads and the (active) power generated by the synchronous machines by a coefficient λ.To alter the nominal value of the power loads, λ was varied continuously in time according to the profile shown in the upper panel of Fig. 6.This profile was sampled every 15 min: this same time basis was used to update the set points of the generators.The effect of energy production by Fig. 6.Upper panel: time evolution of the λ coefficient used to overload the IEEE 39-bus system.Lower panel: the overall active power supplied by the solar plants that were added to the IEEE 39-bus system.three (aggregated) solar plants connected to bus 8 , bus 24 , and bus 27 was emulated by reducing the absorbed active power of the loads connected at those buses.Load power was decreased by the same active power supplied by those (aggregated) solar plants (each one supplying one third of the electrical power reported in the lower panel of Fig. 6).To balance generation, the power supplied by these plants was subtracted from the power of synchronous generators.In particular, the active power supplied by G 3 , G 7 , and G 8 was varied in proportion to their power ratings.Concerning point (ii), the power rating of each one of the virtual synchronous generators (i.e., the solar plants connected to bus 8 , bus 24 , and bus 27 ) is 100 MW.We regulated these generators so that each of them provided, through the controllers of their inverters, a total virtual inertia equivalent to H = 10 s only from 8:15 am to 5:30 pm.
Results are shown in Fig. 7, using again violin plots.The left panel refers to the time interval in which the solar plants are not connected and thus the inertia constants of their gridforming converters is zero.The right panel corresponds to the 8:15 am -5:30 pm time interval.In the high-inertia case, the ǫ H % percentage error for GFM 1 , GFM 2 , and GFM 3 is equal to 0.29%, 1.23%, and 0.60%, respectively.In the lowinertia case, ǫ H % is almost the same for the three grid-forming converters and amounts to 0.02%.
The estimation is accurate even though we used only the measurements of the current flowing through the lines connecting bus 1 − bus 39 , bus 9 − bus 39 , bus 3 − bus 4 , bus 16 − bus 17 , and bus 14 − bus 15 , (i.e., the branches indicated with dashed lines in Fig. 2), which are far away from the buses connecting the solar power plants.As reported in the literature, these lines split the IEEE 39-bus system in areas.
The measurements used in this case highlight that we do not necessarily need to put PMUs in front of synchronous generators or grid-forming/following converters.However, we do not have a systematic method that clearly indicates the optimal locations where PMUs should be connected to provide measurements used during the estimation process.This is an open issue that we will cope with in future work.
Moreover, we observe the birth of additional "swing equation" modes due to the presence of this kind of the virtual Fig. 7. Violin plots of the inertia constants estimates of the grid-forming converters GFM 1 , GFM 2 , and GFM 3 (all of which were set to 10 s from 8:15 am to 5:30 pm).Left and right panels refer to the low-and high-inertia cases, respectively.The time window is 15 min.
synchronous generators [44].These, in fact, induce inter-area oscillations that are clearly identified by peaks in the frequency responses by frequency scan.Modifying the inertia constant of these elements changes significantly the variance of the measured electrical quantities and, thus, helps increasing the accuracy of the proposed technique.
Figure 8 shows the estimation of the inertia of the system across the full day of measurements.This illustrates the effectiveness of the proposed technique to continuously estimate the inertia by stochastic fluctuations.In Fig. 8, the black traces are single trials while the red one is the median value over 50 trials.All realisations give accurate results with only a few time instants being affected by a relatively higher error.
3) Inclusion of Grid-Following (GFL) Converters: In the third and last scenario, we discuss the effect of the frequency droop control of grid-following converters on the estimation of the inertia with the proposed technique.As is well-known, the frequency droop control of the converter is equivalent to a second order model, and is thus similar to that of synchronous machines, except for the fact that it shows a large damping and approximately zero inertia [45].In this scenario, thus, we repeat the simulations of the modified IEEE 39-bus system with inclusion of the same solar power plants considered in the previous section.However, assume that these plants are now equipped with conventional frequency droop controllers.
Table II shows that the estimation of the damping/droop and of the inertia are strongly correlated.An increase in the droop coefficients of the frequency controllers leads to a lower inertia constant estimation.This result had be expected as a lower droop makes the response of the controller slower and thus its time scale tends to overlap with that of the inertial response.The estimation of the inertia constant and of the damping coefficient are performed independently, since two separate least-squares problems are solved: one fits the inertia constant and a second one fits the damping coefficient.
The main conclusion that can be drawn from the results shown in Table II is that the value of the inertia constant alone is not sufficient to define the ability of a device to regulate the Fig. 8. From top to bottom: inertia constant estimates over time of the gridforming converters GFM 1 , GFM 2 , and GFM 3 (all of which were set to 10 s from 8:15 am to 5:30 pm).The results of the 50 estimation procedures are depicted in black.In red we report the median value computed every 15 min over the 50 trials.frequency.If only the inertia constant is estimated, in fact, one would conclude that the case with higher droop coefficient is the one with lower frequency containment, which is not.
On the other hand, if both parameters, namely inertia and damping, are estimated, then one can conclude, correctly, that the system with higher droop provides "more" fast frequency control than inertial response.On the other hand, lower droop values increase the ability of the system to have an inertial response but lead, as expected, to a lower capability of providing a service as fast frequency controllers.

B. All-Island Irish Transmission System (AIITS)
In this section, the proposed method is applied to a realworld complex system.The model of the AIITS considered here consists of 1479 buses, 1851 transmission lines or transformers, 245 loads, 22 conventional synchronous power plants with AVRs and turbine governors, 6 power system stabilizers, and 169 wind power plants.Note that the secondary frequency control of the AIITS is implemented manually and, thus, no AGC is considered in the model.Both load fluctuations and wind speeds are modeled as OU processes.The resulting set of DAEs for the AIITS includes 2118 state variables (663 of which are stochastic processes) and 6175 algebraic variables.Both the active and the reactive power of the Z = 245 loads is altered by an OU process characterised by v = 0.1, while σ is defined in such a way that the standard deviation of those processes is 0.5% of the nominal active or reactive power of each load.Finally, the standard deviation and drift term of the OU processes modeling the wind speeds are set to 0.002 and 0.02, respectively.The time domain simulations of the AIITS are obtained with the software tool Dome [46].
To test the applicability of the proposed method, we estimated the H inertia constants of 11 conventional synchronous power plants chosen to be representative of non-homogeneous H constants and power ratings (see Table III).An optimization problem was solved every 15 min of simulated time by collecting 100 estimates per synchronous generator.Making reference to the diagram in Figure 1, the subset Z contains only the magnitude of the voltage at the buses where the generators are connected.Hence we did not resort to the phase of such voltages, to any current of the network, or to any rotor speed of the synchronous generators.
The ǫ H % relative percentage error spans from 0.22% (for G 7 ) to 5.42% (for G 10 ).The quality of the results, measured through ǫ H % , is similar to that obtained for the IEEE 39bus system.For the sake of readability, the violin plots in The distance of the upper and lower adjacent values (green solid dots in Figure 9) from the expected normalized value of H is lower than 0.22 for all the 11 synchronous generators, which confirms the robusteness of the proposed method when applied to a real-world complex system with several hundreds of buses and dynamic devices.

VI. CONCLUSIONS
This work proposes a technique to estimate the inertia based on the variance of measurements.Assuming a given model of the power system and the dynamics of the generators, the technique fits this model to the measurements by solving a leastsquares problem.Simulation results show that the proposed technique, given an adequate filtering of the measurements, is accurate even in a real-world complex system and can take into account both conventional and virtual synchronous machines, as well as the effect of frequency droop controllers.However, the correct interpretation of the effect of droop controllers is consistent only if also the damping/droop coefficients are estimated, which the proposed technique duly allows to obtain.
In particular, we believe that the proposed technique can be useful to reward the frequency support provided by nonsynchronous devices.In turn, we recommend that the reward should be evaluated based on the effect that the controllers of these devices have on the system, not on the actual implementation of the control itself.
Future work will focus on further investigating the theoretical and practical aspects of the proposed technique.For example, relevant questions are how to improve the accuracy of the estimated inertia constants through (i) an optimal number and location of the set of measurements and (ii) of the filtering of the measurements.

Fig. 1 .
Fig. 1.The flow-diagram of the proposed inertia estimation procedure.

Fig. 2 .
Fig. 2. The single-line diagram of the IEEE 39-bus power system.

Fig. 9
report the statistics of the H estimated inertia values normalized w.r.t. the H nominal inertia of each generator.So doing, for each violin plot in Fig. the reference value is 1.

Fig. 9 .
Fig. 9. Violin plots of the inertia constants estimates for some generators of the AIITS, normalized w.r.t.their effective value.

TABLE I SYNCHRONOUS
GENERATORS H AND P RAT

TABLE II ESTIMATION
OF THE EQUIVALENT INERTIA AND OF THE DAMPING/DROOP FOR THE GRID-FOLLOWING CONVERTERS IN THE MODIFIED IEEE 39-BUS

TABLE III SYNCHRONOUS
GENERATORS H AND P RAT OF THE AIITS