Online inertia estimation for power systems with high penetration of RES using recursive parameters estimation

NewZealandAidProgrammeNewZealandScholarshipsforinternationaltertiarystudents Abstract Low and time-changing inertia values due to the high percentage of renewable energy sources (RESs) can cause stability problems in power systems due to rapid frequency instabilities. Inertia monitoring will assist operators to apply suitable actions and more proper control methods to alleviate stability issues. Therefore, this paper proposes an online method to estimate the total inertia of a network using a recursive least-squares approach. The proposed method uses network measurements with a non-recursive sys-tem identiﬁcation approach to initially estimate the network hypothesis model. Then, the recursive method is used together with time changing measurements to recursively estimate model parameters and extract online inertia estimates of the network. During the estimation, the method does not need to store previous data after each sample step; therefore, the computation burden is signiﬁcantly reduced. More importantly, the technique incorporates the use of available electromechanical oscillation modes in the system, which are linked with system parameters, to determine the network inertia estimates. The applicability of the proposed method has been validated by numerical simulations of the IEEE 39-bus network and the aggregated New Zealand network with its actual inertia data.


Background and motivation
Inertia has a great impact on system dynamics.Inertia takes part in deciding the ability of the network to retain stability when exposed to power imbalances.Conventionally, inertia in power systems has been provided mainly by synchronous generators [1].However, due to the increasing renewable energy sources (RESs) penetration level to power networks, some synchronous generators are replaced by converter-based RESs.Consequently, the inertia constant in power system has reduced.For instance, [2] presents the inertia decrease for different countries between the year 1996 and 2016.The study shows a decrease of inertia value in Denmark by 60%.On the other hand, [3] estimated the inertia value of the UK network to be 9 s by 2008 and anticipated to be 3 s, around 67% decrease, by 2020.These are remarkable decreases of inertia in power systems.Therefore, the study of the power system's inertia is becoming vital.
While inertia is only essential in the first instants after the events of power imbalance, its influence is critical in the network as it suppresses frequency deviations under sudden power imbalances.When the inertia of the network decreases, the standard indicator is the rate of change of frequency (RoCoF), which increases.If the RoCoF increases beyond a critical value of a specific network, RoCoF protective relays will operate to isolate a generating unit from the network.Isolating a generating unit from the network may result in cascading failures and possibly blackout.To address the problem of the increased RoCoF because of the high penetration of RESs, some countries such as Denmark, Germany and UK have updated their critical RoCoF relay settings from 0.5 to 2.5 Hz/s to accommodate high RoCoF in their networks [4].Therefore, the challenges related to reduced inertia in power systems, especially during contingencies, should be given research attention [1,5].
It is traditionally known that for a network with multisynchronous generators, its system equivalent inertia constant can be obtained by using Equation (1) [6].
where H Total and S Total are the total inertia constant and the capacity of the entire system, respectively; H i and S i are the inertia constant and the capacity of the ith synchronous generator, respectively; N is the number of generators in the network [5,7].It is apparent that, for the conventional power system, the system inertia could be calculated by using Equation (1).However, Equation (1) may be unsuccessful at giving an accurate estimate of the actual network inertia in the presence of a high percentage of RESs.The RESs generation units have non-synchronous inertias called synthetic inertias.The existence of synthetic inertia in the power system can further complicate the use of Equation (1) in inertia calculation.Failure to estimate the network reduced inertia might result in difficulties in both network operation and control schemes planning [8,9].An Australian blackout in 2018 as well as UK's and the Island of Tenerife's blackouts in 2019 are good examples of the latest challenges related to a failure of estimating the values of network inertia [10].Other problems related to the high percentage of RESs in power systems are presented in [11][12][13].Some research works such as in [14,15] have addressed the problem of low inertia in power systems by introducing synthetic inertia as a control strategy to replace reduced conventional synchronous generator's inertia.As most synthetic inertias are dependent on stochastic RESs, the future system inertia expected to be a time-varying quantity in power systems.The high probability of inertia becoming a time-varying quantity due to a high percentage of RESs in the network prompts more consideration of power system stability and reliability.As a result, for the safety and reliability of modern and future power systems with low and time-varying values of inertia, online estimation and monitoring of inertia values as well as evaluation of frequency response in power systems are necessary.Prior awareness of inertia values in power systems will help in planning for frequency response in power systems.Besides, the deployment of proper frequency containment services depending on the level of inertia at a given time in power systems can be planned in advance.Consequently, online inertia estimation techniques need to operate in real-time to quickly quantify the risks of power systems blackouts.

1.2
Literature review A power system with a significant percentage of converterbased RESs has low inertia that may be inadequate to immediately and effectively respond to system dynamics when the network is exposed to power imbalances.The low inertia system is prone to the high RoCoF during contingencies and, hence, large frequency deviations.The high RoCoF may lead to frequency instabilities and consequently blackouts.Similarly, since most of RESs integrated into modern power systems are weather dependent, the number of committed synchronous generators over time will be affected by weather conditions.This tendency, consequently, leads to challenges in the planning, operation and control of networks [8,[16][17][18].Subsequently, maintaining the stability of the grid is problematic as the penetration levels of the converter-based RESs keep increasing in the network.
The decrease in the percentage of inertial machines, which is a result of the increase in the percentage of stochastic RESs in the power system needs more research.The research is dedicated to ensuring that the stability and flexibility of the modern and future power grid are maintained [16,17].For the safe and reliable operation of modern and future power systems with low and time-changing values of inertia, online estimation and tracking of inertia in power systems are essential.Besides, real-time assessment of the frequency response in power systems is also important.By successfully assessing the frequency response in real-time, appropriate measures such as control schemes can be planned [19].Generally, quick and continuous information of the inertia value in the network will help power system operators (PSOs) to plan and act either before the contingency or very rapidly after the contingency with appropriate measures.
Even though inertia estimation methods have been implemented in the traditional network, the methods may not be suitable for modern and complex networks.The methods used traditionally include the well-established equation expressed by Equation ( 1), probing test [20,21] and transient test [22]; all these techniques estimate the inertia constants of the networks predominantly driven by synchronous generators.As synthetic inertia is becoming a reality in power systems, hence making the total inertia a time-varying quantity, traditional methods would not work in estimating variable inertia in modern power systems.To address inertia estimation in modern networks, new methods have been proposed recently [5,19,[23][24][25].However, the majority of them are still dependent on Equation (1) and based on historical data analysis as well as large data sets after frequency events.Furthermore, the system inertia estimation in [26] is conducted by using an equivalent system model.Unfortunately, the model requires a large amount of recorded data, which slows down the entire process.Therefore, this approach is not suitable for online inertia estimation.Similarly, [23] uses historical data and the Bayesian approach in the methods to estimate the aggregated inertia of the network.But the approach is prone to computational problems owing to using a large and complex computational model that stores large historic data during computation.Therefore, it is unsuitable for online inertia estimation.Besides, [27] presents a method that has a drawback of downtime revising and updating the entire system model when considering the addition or exclusion of a generator or load in the system.Due to this drawback, the method is irrelevant for online inertia estimation.A downside of decreased precision due to phase step issues as clarified in [28] makes the methodology for inertia estimation in [29] lack credibility for online inertia estimation.
Ref. [30] proposes an online inertia estimation technique, which overlooks the variable inertia due to the high percentage of RESs in modern networks.Furthermore, the method proposed in [31] is also built on the Bayesian approach, which is associated with high computation burden as analysed in [32].The computation burden associated with the Bayesian approach makes the method impractical for online inertia tracking in power systems.Alternatively, [32] proposes two inertia estimation techniques for estimating synchronous and nonsynchronous generators' inertias.The techniques are mainly developed based on the approach in [33] and [34].However, the techniques are limited to individual devices only.
The challenges revealed in the analysed literature above lead to the development of a compact online inertia estimation and monitoring technique.The developed technique should consider time-varying inertia for the entire network, should not suffer from phase step issues and should have minimum computation burden.

Novelty and organisation
To overcome the above-mentioned vulnerabilities and to provide an accurate and reliable online inertia estimation and tracking, this paper proposes a recursive parameter estimation and online inertia tracking in power systems.The proposed technique has the following merits: • Reduces computation burden significantly as it does not need to store previous data after each sample step during estimation.This paper covers the gaps that exist in most proposed estimation methods such as having high computation burdens, considering only electrical model information and leaving behind electromechanical oscillations, which limit most of the inertia estimation methods to only post-event analysis.The remaining parts of this paper are structured as follows.Theoretical background is presented in Section 2. In this section, the role of inertia in frequency stability is discussed.Furthermore, dynamic modes and eigen structures in relation to inertia extraction in power systems are introduced.Section 3 describes the proposed online inertia estimation method based on the recursive leastsquares approach.Moreover, the application of the proposed technique on the case study network is presented and examined in Section 4. The application and validation of the proposed technique on the real network data is examined in Section 5. Lastly, Section 6 presents the conclusion and briefly underlines the future work.

Inertia in power system dynamics
Inertia has a critical role in maintaining the stability of the network after power imbalances [1,35].Inertia determines how fast or slow the frequency can change after power imbalances in the network [5,23,36].To describe the behaviour of frequency response in networks, the traditional swing equation is very important.When the traditional swing equation is simplified and normalised with respect to the unit values as depicted in [9], the transfer function (TF) of the system expressed in Equation (2) can be obtained.
where G (s) is the TF from ΔP e to Δ f , ΔP e is the active electrical power deviation, Δ f is the frequency deviation of the rotor of the generator, H is the inertia constant, D is the damping coefficient and s is the Laplace transform operator.As the PMUs are employed to measure the network's parameters, a system estimation technique can be employed to estimate the TF of the network.Then, network parameters such as inertia can be estimated.
For a network with several synchronous generators, the total inertia constant of the entire network can be obtained by employing Equation (1).Also, a network with several synchronous generators can be combined and assumed as a single generator network [19].To achieve the combined network, all generators' output active powers are aggregated to give a single value of active power.The frequencies of different buses can be averaged at a selected bus to have the so-called centre of inertia (COI) frequency.Traditionally, the COI frequency of the aggregated network is given by Equation (3) [37].
where, f COI is centre of inertia frequency of the network, H i and f i is inertia constant and frequency of the i th synchronous generator, respectively, and N is the total number of synchronous generators in the network.
For traditional networks whose parameters of the installed generators are known in advance, the approach in Equation ( 3) can be implemented when considering a multi-generators network.However, due to the nature of the modern grid with significant RESs penetration and synthetic inertia inclusion, the inertia of some components in the network may not be known in advance.The advanced wide-area measurement systems (WAMS) that are equipped in modern networks can provide measurements of different units in the network, which can facilitate real-time tracking of each component added in the network [38].The information of each component can be obtained and immediately incorporated into the network.Therefore, the sum of active power deviations and the frequency responses, which are tracked in real-time, are used as inputs for the proposed adaptive algorithm to calculate COI frequency in actual networks.The algorithm in this proposed method is designed in such a way it adaptively adjusts to incorporate reduced inertia in the network.

Dynamic modes and eigen structure analysis
An interesting interpretation of system dynamic behaviour can be achieved by examining the eigenvalues of a dynamic matrix A of the state-space model of a power system [39].As matrix A describes the dynamics of the system as characterised by its eigenvalue, the system's inertia constant can be extracted from this matrix.To do this, the dynamic modes of the system are found and then the eigenvalues of matrix A, which are poles of the system, are calculated [40,41].The general form of a dynamic system can be expressed as in Equation ( 4).
where  represents the vector of the state variable and A is the state matrix.If it is assumed that matrix A is diagonalisable with eigenvalue decomposition, then the estimated state variable can be given as in Equation ( 5).
x = Λ Γm where  is a set of functions obtained from the system data, which physically represents the oscillation of the system, Γm represents the estimated row vectors containing the temporary coefficient evaluated at each observation, and ing of empirical Ritz eigenvalues  j of the dynamic matrix A [42].Therefore, the estimated state variable x can then be expanded in a linear combination of modal components as in Equation ( 6).
x ≈ where ã j is a set of temporal amplitudes,  j is a set of dynamic modes, and  j is the associated eigenvalues of the dynamic matrix A of the system model.In this way, the eigenvalues and eigenvectors of the state matrices can then be found.Since there is a cross-coupling between the state variables and the dynamic matrix, it is reasonable to assume that, for n th order system, the homogeneous response of each of n state variables x i (t ) is a weighted sum of n exponential components given by a derivative of a state variable in Equation ( 7) [43].
If a set of m i j and  j can be found that satisfy Equation ( 7), the assumed exponential form is a solution to the homogeneous state equation.In this case,  j is the j th eigenvalue and m j is the corresponding eigenvector of a given square dynamic matrix.Since system dynamics correlate with eigenvalue and eigenvector, any change in dynamic parameters can indirectly change the eigenvalue and eigenvector through the state variable.By determining the modes of the system and analysing the eigenvalues of the dynamic matrix A, network parameters such as inertia can be determined.

THE PROPOSED ONLINE INERTIA ESTIMATION TECHNIQUE
The proposed online inertia estimation method is summarised in the algorithm flow diagram presented in Figure 1.It should be noted that PMUs are used to measure power system data.The recorded data is passed through a non-causal Butterworth lowpass filter with 0.5 Hz cut-off frequency to attenuate the higher frequency components that may impair the inertia estimation.The filtered recorded data is then used to estimate the timevariant TF model of the discrete-time process as the hypothesis model of the system.From the hypothesis model identified, the recursive least-squares method is applied to recursively estimate the parameters of the model.A model reduction is then performed, followed by a transformation to the state-space model.From the state space representation, the dynamic matrix A can be obtained from which the eigenvalues and eigenvectors can be attained to extract the inertia estimate of the system.

Identification of the hypothesis model
The PMU network measurements are categorised as input vectors (u 0 , u 1 , u 2 , … .um ) and output vectors (y 0 , y 1 , y 2 , … .ym ).The PMU measurements are evenly sampled with the sampling period T s = 1/f s , where f s is the sampling frequency.Considering a time t at a discrete-time index k, which corresponds to the maximum number of data points, t k is given by t k = kT s .The obtained data vectors are used to attain the estimated hypothesis model of the network.For the identification procedure, a non-recursive least squares approach is employed to obtain a system hypothesis model.The hypothesis model is represented as a time-variant parametric system in the form of G (s) in Equation (8).
where b k (for k = 0 to n-1) and a m (for m = 0 to n) define model parameters in the numerator and the denominator, respectively.n − 1 and n are the highest orders of the operator "s" in the numerator and the denominator, respectively.
To use discrete data sets for model identification, the s domain TF Equation ( 8) must be transformed to z domain as given in Equation (9).
where a 1 … .a n , b 1 … .b m are system parameters, n and m are real numbers, while z is a forward shift operator.For a decreasing time index (k-1), different data samples can be incorporated into Equation ( 9) to estimate the parameters a i and b i of the hypothesis TF.However, Equation ( 9) must be in the difference equation form as given by Equation ( 10) [44].
where y(k) and u(k-1) are the output and input discrete data samples at time indexes t = k and t = (k-1), respectively.Equation (10) can, therefore, be reordered to get Equation (11).
From Equation ( 9), m and n determine the number of model parameters for the numerator and the denominator, respectively.For N data points available in the experiment but only l data points to be used, i.e., (l < N and l > n), then, general matrices Equations ( 12), ( 13) and ( 14) are defined for the vector equations representation of measurements data and system parameters, respectively.
Combining Equations ( 12), ( 13) and ( 14) to represent Equation (11) with vectors of several data points (l), Equation ( 15) is obtained.Equation ( 15) connects several discrete data sets of the network to estimate  system parameters.
It is also possible for N data sets to be used in Equation (15).For this case, when l = N, Equation ( 16) is obtained.Equation ( 16) is comparable to Equation (15) with an exception that all data sets are used in Equation (16).When more data sets are used, it makes Equation ( 16) to be overdetermined, and therefore, the system parameters identification is simplified [44].

Recursive model parameter estimation
Before the recursive least square approach is applied, the nonrecursive least square method is used to obtain the hypothesis model of the network.It should be noted that the non-recursive method is used with offline measurements of the network to estimate the network's model and its order dimension.Then, the recursive least square approach is employed for online inertia estimation.For online inertia estimation in this proposed approach, only parameters of the model should be recursively estimated and updated concurrently with the measurement process.After each sample step or a certain number of sample steps, the parameter estimates should be available.The recursive approach is appropriate in online estimation as it reduces the computation burden and provides an update of the parameter estimates after each sample step.The previous data do not need to be stored.However, the recursive least squares approach is derived based on the non-recursive least-squares method.The parameter estimate Θ(k) for the non-recursive method of least squares for the sample step k is given as Equation ( 17) [45]. where where  represents the data matrix of time sampled input and output measurement of the system as in Equation ( 21).
Likewise, the parameter estimates Θ(k + 1) for the sample step k + 1 can be given as in Equation (22).
where θ(k + 1) is the new parameter estimate, θ(k) is the old parameter estimate and P(k + 1)(k + 1) is the correction vector.From Equation (26), y(k = 1)is the new measurement and  T (k + 1) θ(k) is the predicted measurements based on the last parameter estimate.
From Equations ( 25) and ( 26), a recursive formulation of the estimation problem can be realised.To calculate P −1 (k + 1) recursively as per Equation ( 24), it needs one matrix inversion per update step.However, to prevent computation burden by inversion of a matrix each time step, Equation ( 24) can be rewritten as Equation (27).
Since the term ( T (k + 1)P(k)(k + 1) + 1) is a scalar quantity only, there is no need to invert a full matrix any longer.
When Equation ( 28) is combined with Equation ( 25), they return a recursive method of least squares as Equation (29).
where  (k) is the correction vector given by Equation (30).
From Equation ( 27), it follows that Equation (31) can be obtained.
The recursive method of least squares is, therefore, given by the three equations above, which need to be performed in the sequence of Equations ( 30), ( 29) and (31).
It should be noted that only the network's model parameters, not the whole model of the network, is updated recursively.If the model of the entire network is updated using a recursive approach, then the computation burden is significantly experienced in the proposed method.To avoid the computation burden, only the parameters of the hypothesis model that show a significant effect in the estimation process are updated recursively.The order dimension of the identified hypothesis model is kept constant throughout the simulation time.

Model reduction
Big dataset-driven system identification methods are associated with huge and complex models.Besides, the obtained models have very high degrees of freedom that represent the network.The high order models obtained, which are also complex to analyse, result in high computational burden algorithms.Due to these reasons, the resulted high order models need to be reduced.However, the application of reduction methods should be done in such a way it maintains the essential dynamics of the system.To achieve low order models from high order models, the singular value decomposition (SVD) technique can be applied.To obtain the decomposed model, the SVD technique is applied to the oblique projection vector of the model  i , which is defined in Equation ( 32) [46]. where As a result, SVD is given as Equation (33).
where W 1 and W 2 are the identity weighting matrices, U and V are orthogonal matrices, while matrix  is the diagonal matrix with real entries.Likewise, SVD can be presented as a matrix multiplication by Equation (34).The decomposed order of the system is chosen as  1 provided  1 major singular values are identified in matrix  [46].
where U 1 and U 2 are the components of U , while V 1 and V 2 are the components of V .Then, the TF model of the network at each sample time step is transformed into a state-space representation.From the state-space model, the dynamic matrix A of the system that contains the dynamic behaviour of the system can be extracted.

Inertia extraction
The dynamic matrix A can be represented by an estimated general form of a dynamic system as in Equation (4).Then, the aggregated system dynamic behaviour is estimated by the swing Equation (35).
2H ω + D = P m − P e (35) Since the mechanical power is constant during the transient state, the linearized form of the swing equation can be represented as Equation (36).
For any output signal y(t ), the system model comprises a sequence of the dynamic behaviour of  modes.This sequence can be approximated in a linear model around a stable operating point where TF is written as Equation (38).
where j = 1, 2, … , m;i = 1, 2, … , ; with m being the number of outputs and  is the number of modes.R i is the residual, which is directly coupled with the amplitude of each mode to generate the output signal.When the output signal of the system is sampled at a constant sampling rate Δt , the sampled signal is represented as Equation (39).
where k denotes the samples, and z i =   i Δt is the discretization of model variable z(t ), with  =  + j  (eigenvalue).The set of snapshot matrix of A representing the dynamic modes of the system is given by Equation (40).
Since the states of the system in the time domain can be estimated as a linear combination of the terms  i e  i t , it can be seen that eigenvalues  i can be linked to the physical parameters of the system by a linear combination factor (LCF)  i [43].Considering that dynamic modes are observed during the dynamic response of the system and each observed variable is recorded by N samples with interval Δt , the dynamic modes and eigenvalues can be used to approximate the dynamic characteristics of the system.
To estimate the inertia, the frequency and power deviations are recorded as input and output of the reduced power system model in samples to construct the snapshot matrix.Considering that the LCFs for the frequency deviation and active power variation for i th dynamic modes are  i, Δ f and  i, ΔP e , respectively, and  i is the i th eigenvalue, the time progressions of the dynamic mode reconstructions for the frequency and power deviations can be given as Equation (41) [47].
where  i represents initial value coefficient corresponding to the i th eigenvalue.Inserting Equation (41) in the linearized swing Equation ( 36) and neglecting the damping, the linearized dynamic equation is reconstructed by the dynamic mode as presented in Equation (42).
i e  i t b i  i,ΔP e (42) Rewriting Equation ( 42) for arbitrary t , e t can be eliminated from Equation (42) to obtain Equation (43).
The effective inertia of the system H e can therefore be obtained by solving Equation (43).Since the equation is linear and does not contain the derivative of the rotor speed, the proposed methodology can accommodate large-scale power systems with high dimension.
Also, to avoid some eigenvalue computation problems such as singularity in state matrices and difficulty in finding a complete set of unstable poles of the state matrices in analysing big power systems, the proposed algorithm is incorporated with the following additional attributes: • The algorithm extracts only significant eigenvalues.
• The algorithm incorporates subspace accelerated Rayleigh quotient iteration to speed up eigenvalue computation, improve convergence, and compute poorly damped eigenvalues.• The algorithm is designed to avoid any singularity in state matrices.
Figure 2 presents the generalized flow diagram of the proposed method.Likewise, Tables 1 and 2 provide additional clarity on the algorithm of the proposed technique.

SIMULATION CASE STUDY: IEEE 39 BUS NETWORK
In this section, a modified IEEE 39-bus system is considered to demonstrate the applicability and accuracy of the proposed method.The system is a typical interconnected power system, which can be used to validate a proposed technique in power system studies.The network has 10 generators, 12 transformers and 19 loads.Besides, the network is interconnected by 34 lines as depicted in Figure 3, which includes RESs.DIgSI-LENT PowerFactory 2019 is used to model, test and simulate the dynamic behaviour of the system.

Data pre-processing and model preparation
The case study network is used to obtain data sets to be used for the proposed online inertia constant estimation method.The real theoretical inertia constant value (H) is entered for Estimate the hypothesis of the system TF model using (10) to (17) applying the non-recursive least-squares approach.3. Start a recursive least squares method to estimate the parameters of the model by running ( 31), ( 30) and ( 32) sequentially.4. Perform the SVD of the weighted oblique projection to determine the order by inspecting the singular values in S using: W 1  i W 2 = USV T 5. Transform TF to state-space model and extract system dynamic matrix A.

Part II: Eigenvalue and inertia extraction
1. Determine the eigenvectors and eigenvalues from the obtained dynamic matrix A of the system.
Determine the quantitative relationship between the eigenvalue of the dynamic matrix A of the model and the linear estimated swing equation.
i e  i t b i  i, ΔP e 3. Extract the estimated inertia of the system by defining the dynamic mode corresponding to the eigenvalue that corresponds to the inertia of the system.
2H  i, Δ f b = − i, ΔP e b 0 To perform these steps, the algorithm to estimate the total inertia constant for the aggregated network is presented on Tables 1 and 2 and summarize in Figure 2.
For the aggregated network, permanent magnet synchronous generator (PMSG) wind power plants with synthetic inertia controls are also integrated into the model to replace some of the synchronous generators.The PMSG wind power plants with synthetic inertia controls are integrated step by step to emulate different levels of penetration of RESs into the grid as presented in the modified IEEE 39-bus network in Figure 3.
For each level of RESs penetration, the algorithm is used to estimate the available inertia in the system.
To capture the network transients, root mean square/electromagnetic transient (RMS/EMT) is used as the simulation type in PowerFactory.For the modes and eigenstructure analysis, Model/Eigenvalue analysis is used as the simulation type in PowerFactory.The case study network is aggregated to perform as a single generator network.Furthermore, all the PMUs measurements from all generator buses are aggregated to perform as a single machine infinity bus network.
Before performing the algorithm for online effective inertia estimation of the network, it is essential to construct the operational structure of the system and to determine the required measurements and the COI where the effective inertia can be FIGURE 4 Step response of the IEEE 39-bus network model showing frequencies at different generator buses and COI frequency estimated.Equation ( 3) is used to locate the relatively acceptable COI, and frequency is measured from this COI.By definition, the COI frequency is the frequency of the bus at or around the centre of the network relative to other buses of the network [37].For the aggregated IEEE 39-bus network, Figure 4 presents the frequency of bus 14 as the selected best bus to represent the COI frequency using Equation ( 3), where all frequency measurements for effective inertia estimation are taken.
Based on power variations at different buses, which further cause frequency responses of the network, the data are collected as presented in Figure 5. Before the model identification step, the data is pre-processed to improve the identification process and the estimation performance.In this case, the algorithm part I procedure of the proposed method, as presented in Table 1, is performed to estimate the hypothesis model and recursive parameters of the system model.The model is crossvalidated using data collected in a different time span.The crossvalidation is intended to authenticate the relevance of the model in representing the dynamics of the system.The cross-validation for the model is presented in Figure 6, which gives a fitting ratio (FR) of 96%.

Online inertia constant estimation for each generator
Based on the recorded rotor speed response at each generator as shown in Figure 7, the input snapshot matrix for the dynamic modes of generator G01 and the dynamic matrix A extraction as presented in Equation ( 4) is constructed.The examples of extracted dynamic modes of G01 are presented in Figure 8.The eigenvalues are extracted from the dynamic matrix A as presented in Figure 9.According to the suggested approach in [39], estimating inertia for each generator requires at least a fourthorder estimated state-space model.This means that at least four distinctive eigenvalues must be attained.Figure 9 shows a reasonable number of obtained stable pair of conjugate and real eigenvalues.By using the FR comparison, the results show that   the extracted eigenvalues and eigenvectors are acceptable by achieving the FR of more than 95% [5].In this simulation, the inertia values are provided by each synchronous generator.The proposed algorithm is intended to estimate the inertias from different generators.The comparison between the actual inertia and the average estimated inertia results in an estimation error.Table 3 presents the average estimated inertia for each generator in the network with the associated percentage error for each estimation.Figure 10 presents pictorial comparisons of the actual and the average estimated inertia for each generator in the network.The corresponding error graph for each generator is depicted in Figure 11.The worst-case is observed at bus 38, where G09 is located.At this generator, the algorithm returns the inertia estimation with a percentage error of 3.57, which is the worst-case compared to the rest of the generators' estimations.

Online inertia constant estimation for aggregated network
The proposed inertia estimation technique is tested on the aggregated IEEE 39-bus network to further verify its applica-  bility and accuracy.Initially, the equivalent inertia (H eq ) is calculated using Equation ( 3) for all generators in the network as presented in Table 4.
To incorporate RESs in the network, PMSG wind power plants are used to replace the conventional synchronous generators in the network.Based on the topology of the IEEE 39-bus network, a modified system is generated in steps by increasing the percentage of RESs in the network from 0% to 65%.The first scenario considers a purely conventional network driven by traditional synchronous generators only.Then, wind generators replace synchronous generators G08 and G09 at buses 37 and 38, respectively, to incorporate the second scenario of 10% RESs penetration in the network.The third scenario is the replacement of generators G07, G08, G09 and G10 for penetration of 20% RESs in the network.The fourth scenario is the replacement of generators G02, G03, G04, G05, G06, G07, G08, G09 and G10 to obtain a penetration of 40% in the network.Furthermore, the fifth scenario of 60% penetration is obtained by replacing G01 in the network.The last scenario is obtained by replacing G01, G02 and G05 to achieve a 65% penetration of RESs in the network.
For each penetration of RESs in the network, the average RoCoF is obtained.Table 5 presents different average RoCoF values as recorded for different levels of RESs penetration in the case study network used for simulations.
The electromechanical response of each of the modified network scenario is stimulated by power variations to record the frequency response.The measured data of power variations and frequency response of each scenario as presented in Figure 5 are employed in the proposed algorithm to estimate the total inertia of the study case network with the corresponding percentage of RESs penetration.The effective inertia value corresponding to every level of RESs penetration in the network is presented in Table 6. Figure 12 compares the actual inertia in the case study network to the average inertia estimated by the proposed method.The corresponding errors are presented in Figure 13.

Online inertia constant tracking for aggregated network
In support of the applicability of the proposed method, Figure 14 presents the simulation results of time-varying inertia together with a tracking estimation trajectory.The tracking tra-FIGURE 12 Actual inertia compared to average estimated inertia for various percentages of RES in the aggregated case study network FIGURE 13 Actual inertia and corresponding estimation error for the aggregated IEEE 39-bus network with different RESs penetrations jectory is obtained using the proposed method explained in Section 3. Using the recursive model parameters estimation to track the inertia of the aggregated network, it takes a total of 8 s to run a 200 s simulation with a sampling time of 0.01 s.On the other hand, a 1 s simulation of the entire hypothesis model with a sampling time of 0.01 s updated using the recursive method takes a total time of 5.6 h to complete.For these simulations, an office computer (Intel(R) Core (TM) i7-9700 CPU @ 3 GHz, 32GB RAM) is used to simulate both models.A massive difference in computation time can be noted between the two approaches.Hence, the significant reduction in computation burden using the proposed method is verified.
It can also be noted from the simulation result on Figure 14 that the tracked inertia, which is obtained by using Equation (1) for each step of inertia reduction, exhibits a step change for each variation.It should be noted that even after the penetration  14 that the tracking trajectory has an average speed of 1.5 s in tracking the actual inertia changes.This means, after each step change in the actual inertia, the proposed method takes an average of 1.5 s to estimate the new value of inertia change in the network.This is a good speed of estimation given the fact that weather changes might not result in a step change of the actual inertia in the network.The step changes of inertia provided in the simulation are faster than weather dependent inertia changes.For this reason, the proposed method can easily track the weather-dependent time-changing inertia in power systems.

VALIDATION USING REAL MEASUREMENTS DATA
The effectiveness of the proposed technique is tested using data from the New Zealand network.The New Zealand network is one of the smallest networks worldwide [48].It includes North and South islands, which are interlinked with an HVDC line.The network consists of hydro generators, Combined Cycle Gas Turbines (CCGTs) and geothermal stations.The respective ranges of inertias for the generators are 2 to 4 s, 5 to 6 s and 3 to 6 s.In addition, the network comprises wind generation units that do not supply any inertia to the network.Therefore, more integration of wind power units into power systems, which may replace some synchronous generators, reduces the total inertia of the network.
The wind power contribution in the current status of the New Zealand network is small as it contributes less than 5% of the entire installed capacity.The total inertia value of the network is calculated by using Equation (3) as presented in Table 7. Figure 15 shows the normal record of the network frequency  Using Equation (3), the total inertia value of the network was calculated to be 8.02 s.This value was obtained with a wind power contribution of 0.9% in the network.This contribution is too low to make a sound effect of the total inertia variation in the network.The system bases power and frequency used for the calculation of the total inertia constant value in the network are 3000 MVA and 50 Hz, respectively.Table 6 also presents different information for different types of generators at the time of data recording.The field measured data is pre-processed to enhance identification by improving efficiency and accuracy.To remove noise from the data, the signals are passed through an order of ten non-causal, zero-phase low-pass filter set at 3 Hz and then down-sampled to four samples per second [50].The data is down-sampled to increase the estimation accuracy of modes of the generators.
The recorded active power variation and the frequency responses are sampled to represent input data u(t ) and output data y(t ), respectively.Then, the non-recursive least-squares approach steps 1 and 2 from the algorithm in Table 1 are used to estimate the model of the network.Then, the comparison between the estimated and the actual models is done by calculating their FR.For this test network, the resulted FR is 85%.This FR is good, satisfactory and acceptable for the network data used [44].Figure 16 represents the graphs of the actual and estimated models of the network with the FR of 85%.
The proposed recursive least squares method is used for parameters estimation of the estimated model and recursively estimating and extracting the inertia value of the network.Using the proposed method, the estimation of the average inertia constant ( H ) of the network is found to be 7.5 s.Comparing the actual total inertia constant H with estimated inertia as presented in Table 6, the percentage error is computed to be 8.54%.As the percentage error is less than 10%, this confirms and validates the applicability of the proposed method.

CONCLUSION AND FUTURE WORK
An online inertia estimation technique is presented in this paper.First, the inertia estimation is carried out for each generator in the network.Second, inertia estimation is also carried out for different levels of inertia constant in the aggregated case study network.The different levels of inertia constant presented are a result of different percentage of RESs in the network.Further, consistent estimates using the proposed online inertia estima-tion technique are obtained, which are in the range of 0.74% to 3.57% for individual generators in the network.Besides, the error range of 2.39% to 3.80% is observed for the equivalent inertia of the aggregated network.Finally, an error of 8.54% is obtained for the actual data from the real network.Besides, the proposed method proves to be effective in online estimating the power system inertia.On top of that, since the method does not need to store previous data after each sample step, the computation burden is significantly reduced in the proposed technique.More importantly, the technique incorporates the use of available electromechanical oscillation modes in the system, which can be linked to parameters of the power system for online estimation of the network inertia.Future work of this research is to extend the online inertia estimation method while incorporating coordinated synthetic inertia in the network as a support for reduced inertia in future power systems.Likewise, due to the high penetration of weather-dependent generation units, the study of weather forecast effects in the equivalent inertia of the network is important to be incorporated.In the same way, how fast the estimated inertia can be communicated to PSOs for fast stability control of the network will be looked at.

FIGURE 1
FIGURE 1 Summarised diagram of the proposed method

FIGURE 2 TABLE 1
FIGURE 2 Flow diagram of the proposed recursive algorithm for online inertia estimation

FIGURE 3 A
FIGURE 3 A modified single line diagram representing the IEEE 39-bus network

FIGURE 5 FIGURE 6
FIGURE 5 Simulation results of case study network indicating trajectories of active power and COI frequency

FIGURE 7 FIGURE 8 FIGURE 9
FIGURE 7The time-domain trajectory of rotor speed for generators

FIGURE 11
FIGURE 11 Actual inertia and corresponding estimation error for different generators in the case study network

FIGURE 14
FIGURE 14 Online inertia constant tracking for the aggregated IEEE 39-bus network with different RESs penetrations

FIGURE 15
FIGURE 15 Logged frequency fluctuations at the COI of the New Zealand network

FIGURE 16
FIGURE16 An FR of 85% to justify the estimated model when compared to the actual New Zealand network model

TABLE 3
Inertia estimation for different generators

TABLE 4
Generators' actual inertias and the total inertia constant for the case study network

TABLE 5
Average RoCoF values for various percentages of RESs penetration

TABLE 6
Actual inertia compared to average estimated inertia for various percentages of RES in the case study network

TABLE 7
New Zealand network power generating units with related power capacities and inertia constants