Online Estimation of Power System Inertia Constant Under Normal Operating Conditions

An online estimation method for the power system inertia constant under normal operating conditions is proposed. First of all, a dynamic model relating the active power to the bus frequency at the generation node is identified in the frequency domain using ambient data measured with the phasor measurement units (PMUs). Then, the inertia constant at the generation node is extracted from the unit step response of the identified model in the time domain using the swing equation. Finally, with the sliding window method and the exponential smoothing method, the estimated inertia constant is updated in real-time. Compared to the conventional methods using large disturbance data or field test data, the proposed method can estimate the inertia constant under normal operating conditions, and therefore, can provide the tracking trajectory of the power system inertia constant in real-time. The effectiveness of the proposed method is validated in the IEEE 39-bus system. The results show that the relative error of the identified inertia constant is below 5% and the identified inertia constant can be updated within 1s.


I. INTRODUCTION
In the point view of physics, the inertia of a power system is its capability to resist energy fluctuations caused by external disturbances, which, in conventional power systems, are supplied mainly by the kinetic energy stored in the rotating mass of the synchronous generators and quantified by the inertia constant [1].
With the increasing share of the power electronic inverter interfaced renewable energy sources (RESs), some changes have taken place in power system inertia. Many synchronous generators are displaced by RESs, thus leading to a persistent decrease in the conventionally available inertia resources [2]- [4]. Meanwhile, with different control strategies and parameters, various inertia suppliers such as virtual inertia control and energy storage systems are employed to The associate editor coordinating the review of this manuscript and approving it for publication was Lei Chen . improve the power system inertia [5]- [7]. Therefore, the power system inertia, as well as its quantification, not only decrease but also become time-varying [8].
In conventional power systems, the inertia constant is steady over the long term, and therefore, is chosen as a fundamental reference to the design of frequency stability controls and many of the protection relays [1]. When the penetration of RESs becomes higher, the situation is different. For instance, the system frequency may drop badly when the system suffers an active power deficiency, and consequently, the protection and control devices such as under-frequency load shedding or disconnection of generators may be triggered reluctantly [9]. Since there is no mature method to track the power system inertia in real-time, to avoid the malfunction of the protection relays and stability controllers, the transmission system operators (TSOs) must adopt more conservative operational schemes to ensure the stable operation of the power system, which inflates the power system operational costs.
Given the issues mentioned above, the need for online estimation of power system inertia (constant) is highlighted and has drawn consistent attention in recent years. According to the measurement data types, the existing estimation methods can be categorized into two groups: 1) disturbance data methods, and 2) non-disturbance data methods.
For most of the estimation methods based on disturbance measurements, the inertia (constant) is estimated using the swing equation, utilizing post-event data that record the transient active power and frequency characteristics of the generators [10]- [21]. For instance, a procedure for estimating the total inertia of the Great Britain power system was proposed in [13], which calculated the total inertia for the whole system by summing all the estimated regional inertias. An online algorithm to estimate the system inertia after a disturbance based on sliding windows of active power and frequency derivative measurements was proposed in [14]. In [17], an inertia estimation method based on the extended Kalman filter was proposed, which needs to assume the time of disturbance. Based on the transient and steady state characteristics of the frequency response after a disturbance, a method to estimate the equivalent inertia and damping constant simultaneously was proposed in [18]. In [19], an approach for online inertia estimation in the power system network with solar photovoltaic source was proposed using the synchronized measurements from PMUs. In [20], based on the measured frequency response for any arbitrary disturbance in the system, a method to estimate the available inertia in an islanded microgrid was proposed. In addition, some other disturbance data methods have recently been developed. In [21], an approach based on electromechanical wave theory was presented to identify the change of power system inertia distribution. Though both the offline and online methods based on recorded disturbances can estimate the inertia accurately, they cannot achieve continuous inertia estimation as the natural transient events are deficient, and the transient field experiments are expensive.
Compared to the disturbance data methods, estimators based on non-disturbance measurements are quite limited. In [22], a statistical model-based real-time inertia estimation method was proposed, in which the model was trained to learn the features that relate the steady-state average frequency variations and the system inertia. In [23], ambient frequency and active power data were employed to estimate the effective inertia of a power system, where a combined model of inertial response and primary control was identified, and then the inertia was extracted from the impulse response of the model. However, the inertia tracking trajectory could not be provided in real-time in [23]. In our previous work, a closed-loop identification method was proposed for the power system equivalent inertia constant online estimation, which could achieve a precise estimation result, but it required the injection of an additional probing signal [24].
In this paper, we propose a more robust online estimation method for the power system inertia constant under normal operating conditions. Compared to the conventional methods based on transient test, probing and etc., the proposed method can achieve a precise inertia constant estimation without any disturbance event and probing injection, and can provide the tracking trajectory of the inertia constant in real-time. To our knowledge, it is the first work to track the inertia constant under normal operating conditions. In relation to the existing work on inertia estimation, the contributions of this study are clarified in Fig. 1 and summarized as follows: 1) it provides a more stable and precise solution to identify power system inertia constant using ambient signals; 2) it realizes the estimation of inertia constant at different hierarchies (individual generator, area and the whole system) under normal operating conditions; 3) it realizes the real-time online tracking of inertia constant in the time scales of seconds under normal operating conditions, which can timely provide important information for stable operation of the power system.
The remainder of the paper is organized as follows: section II provides the theoretical fundamentals about power system inertia; section III introduces the proposed inertia constant online estimation method; section IV presents simulation results in the IEEE 39-bus system to verify the feasibility of the proposed method; section V concludes the paper.

II. THEORETICAL FUNDAMENTALS A. INERTIAL RESPONSE AND INERTIA CONSTANT
After a disturbance, the frequency response of a traditional power system will sequentially go through three stages: the inertial response, the primary response and the secondary response [13]. During the inertial response, the rotating kinetic energy stored in the synchronous generators is released spontaneously to maintain the power balance, thus reducing the rate of change of frequency (RoCoF).
The rotational kinetic energy of a synchronous generator is usually normalized to an inertia constant, which is defined as the ratio of the stored rotational kinetic energy at the rated rotating speed to the rated capacity of the synchronous generator [1], namely: where H i is the inertia constant of generator i, E i is the rotational kinetic energy of generator i and S i is the rated capacity of generator i. Physically, H i represents the time duration to supply energy for the demand that equals to the rated capacity of the generator, without any additional mechanical input. In a power-electronics-dominated power system, the converter-interfaced RESs are initially inertia-free and the system frequency response does not have clear three-stages. However, virtual inertial response can be obtained with control of electrical converters, then converter-interfaced RESs can provide equivalent inertia for the power system. Besides, the dynamics of inertial response provided by converter-interfaced RESs can be described by a first-order differential equation [25], [26], which is similar to the swing equation of the synchronous generator. Namely, the inertial response of converter-interfaced RESs is mathematically equivalent to that of the synchronous generators [27]. Therefore, when estimating the equivalent inertia constant of the converter-interfaced RESs, it is reasonable to regard them as the equivalent synchronous generators and the proposed estimation method can be employed.
Further, in a multimachine power system, if we consider the other inertia contributors as the equivalent synchronous generators, then the equivalent inertia constant of the entire power system can be calculated as follows: where H sys is the equivalent inertia constant of the system, S sys is the rated capacity of the system and N is the number of generators.

B. SWING EQUATION
The dynamics between active power and frequency of a synchronous generator in a short time frame after a power mismatch can be modeled by the swing equation. For synchronous generator i, considering the damping effects, the swing equation can be written as [1] where P m,i and P e,i are the mechanical power (in p.u.) and the electrical power (in p.u.) of generator i, respectively; f r,i is the rotor electrical frequency (in p.u.) of generator i; and D i is the damping coefficient of generator i. Physically, when suffering a power mismatch, the RoCoF of the synchronous generator is constrained by the inertia constant and the damping coefficient in a short time frame, thus the frequency of the generator cannot change suddenly and the frequency stability can be improved.
As an approximation, the dynamic behavior of a certain area or the system can be represented as an equivalent synchronous generator j, which leads to the aggregated swing equation as follows: where P m,j and P e,j denotes the total mechanical power and the total electrical power of the area or the system (in p.u.), respectively; f j denotes the aggregated frequency of the area or the system (in p.u.); H j denotes the equivalent inertia constant of the area or the system; and D j denotes the total damping coefficient of the area or the system. The theoretical value of H j can be calculated using (2). Under normal operating conditions, all the variables in (3) vary around the steady-state operating point. Therefore, formula (3) can be written as the incremental formulation around the steady-state operating point as follows: Assuming P m,i to be zero and taking the Laplace transform on both sides of (5), we can reformulate the swing equation as a first-order transfer function: where f r,i represents the rotor electrical frequency deviation of generator i, P e,i represents the electrical power deviation of generator i, and G i (s) (s is the Laplace operator) is the transfer function from P e,i to f r,i , which is characterized by H i and D i . Besides, the transfer function from P e,j to f j is similar to that in (6), which is characterized by H j and D j .

III. METHODOLOGY
The methodology proposed in this paper is for power system inertia constant real-time online estimation under normal operating conditions, which is shown in Fig.2. The key procedures include signal selection and preprocessing, system identification, inertia constant extraction and inertia constant tracking. All the key procedures can be realized automatically with low computational burden, and their details are expanded in the following subsections.

A. SIGNAL SELECTION AND PREPROCESSING
According to section II-B, the dynamic model of generator i can be identified using electrical power P e,i as the input and rotor electrical frequency f r,i as the output. However, P e,i and f r,i in the real power system are difficult to measure. As a substitute, P e,i and f r,i can be approximated by the active power output P i and the frequency f i at the generator connection bus, respectively. In practice, P i and f i can be measured by PMUs installed at the generator connection bus. Therefore, we can identify the dynamic model of generator i using active power output P i as the input and bus frequency f i as the output. Namely, the following equation holds: Additionally, it is possible to identify the dynamic model of a certain area (or the system) using the total active power output and the aggregated frequency of the area (or the system) as input and output, respectively. It would be simple to obtain the total active power output by summing up all active power outputs in the area (or the system). To aggregate the frequency, however, the center of the inertia frequency is commonly used, which is an average frequency weighted by the inertia of each node and cannot be measured directly [23].
Here we propose a simplified aggregated frequency to represent the center of inertia frequency. The frequency of the area (or the system) is evaluated by a weighted average of the measured frequencies as follows: where f i denotes the frequency at the generator connection bus, and N j denotes the number of online generators in the area (or the system). Here, the weights w i are determined based on the following consideration: the larger the inertia constant of generator i, the smaller the variations of the bus frequency f i under normal operating conditions. Therefore, the weights w i are defined as the inverse of the variance of bus frequency f i , namely, Under normal operating conditions, the ambient data measured by PMUs are commonly polluted with noise, so signal preprocessing is a necessary procedure to improve identification efficiency and accuracy before running an identification algorithm [28], [29]. First, all the signals are converted into per unit values by dividing their base values. Then, all the signals are detrended by removing their mean values and prefiltered using a noncausal Butterworth low-pass filter. As the typical values for H i are in the range of 2-10 s [1], a noncausal Butterworth low-pass filter with a 0.5 Hz cut-off frequency is suitable to attenuate the higher frequency components that can impair the inertia constant estimation. Finally, the signals are downsampled to the range of 5-10 Hz to avoid the numerical problems when running the identification algorithm.

B. SYSTEM IDENTIFICATION
Generally, the power system is nonlinear. However, under normal operating conditions, the disturbance to the power system is small, so the nonlinear power system can be approximated by a linear state space model around the steady-state operating point. An n th other multi-input-multi-output state space model can be described as follows: where x k ∈ R n is the state vector; u k ∈ R m is the input vector; y k ∈ R l is the output vector; and A∈ R n×n , B∈ R n×m , C∈ R l×n , and D∈ R l×m are system matrixes to be identified. Additionally, ω k and ν k denote random sequences of process noises and measurement noises, respectively. Subspace identification methods are effective algorithms to identify the state space model with ambient data and can be implemented in different ways [30]. As a widely used algorithm in subspace identification [28], the N4SID (Numerical algorithm for Subspace State Space System IDentification) algorithm is used in this paper. The main feature of the N4SID algorithm is to calculate matrix k through oblique projection as follows: where / denotes the oblique projection. Then, singular value decomposition (SVD) is applied on k to determine the order VOLUME 8, 2020 of the identified model. Specifically, the order of the identified model is equal to the number of the dominant singular values of matrix k . The SVD can be partitioned into the following form: where W 1 and W 2 are the identity weighing matrixes. In (11), the insignificant singular values are neglected by removing S 2 as the dominant singular values determine the main dynamics of the system. Finally, the system matrixes A, B, C and D can be obtained by solving the linear equations. Readers can refer to [28] for the details of the N4SID algorithm. Generally, the order of the real system is rather high as it contains many complicated control systems. However, a model with a lower order is identified in this part. Though the order of the identified model is lower than the order of the real power system, it is accurate enough to capture the dynamics of the inertial response. The N4SID algorithm can search the best order for the identified model automatically after setting a range of orders from n min to n max . Empirically, n min can be set as 1 and n max can be set as 10 for inertia constant estimation.
After running the identification algorithm, the reliability of the identified model should be verified by model crossvalidation. The model cross-validation can be performed by comparing the validated outputŷ(t) and the original output y(t). To evaluate the reliability of the identified model quantitatively, the fitting ratio (FR) between the validated output y(t) and original output y(t) is defined as follows: where N is the number of samples.

C. INERTIA CONSTANT EXTRACTION
Theoretically, the state space model includes the inertia constant but as an implicit value, so further analysis should be employed for the identified model. A direct way is to extract the inertia constant from the model itself in the frequency domain after some transformation. First, the state space model can be transformed into a transfer function from P i to f i . Then, the transfer function is reduced to a first order inertia function the same as (7) and finally, the inertia constant H i is estimated together with D i . However, large error may be introduced during the order reduction process, leading to inaccurate estimation of the inertia constant. For this reason, we propose a method to extract the inertia constant from the step response of the identified model in the time domain. According to section III-B, the real power system can be approximated by a linear model and the model is reliable if the model cross-validation performs well. In other words, the identified model can be regarded as the real power system model to some extent. In this perspective, when a certain disturbance is applied to the identified model, the disturbance source can be regarded as the active power deviation P i and the corresponding response can be regarded as the bus frequency deviation f i .
If we disturb G i (s) with an unit step signal, namely, P i = −ε(t) (ε(t) is the unit step function), then f i can be expressed as (13) in Laplace domain.
The equation can be solved directly and be written as follows in the time domain: Then the slope of the unit step response at t = 0 can be calculated as: According to (15), the inertia constant H i is determined by the initial RoCoF ḟ i | t=0 , namely the initial slope of the unit step response of G i (s). Therefore, we can disturb the identified model with a unit step signal and calculate the initial slope of the corresponding response as ḟ i | t=0 . The problem turns into how to estimate ḟ i | t=0 after the disturbance as accurate as possible. As the inertial response activates immediately after the disturbance and lasts a short time, we recommend calculating ḟ i | t=0 with a 500 ms sampleby-sample sliding window 1 over 1-2 s period [13] from the unit step response of the identified model. For power systems of various sizes, the data length of the unit step response for the ḟ i | t=0 calculation can be adjusted properly according to the dynamics of the systems. During each sliding window, ḟ i | t=0 is estimated as the slope of a linear fit to the response. The maximum slope is then taken to represent the ḟ i | t=0 during the inertial response period following the disturbance.
Note that, the discrete time identified model should be converted into continuous time model before evaluating step response. Then, the continuous time model should be checked for stability in s-domain, namely, the real part of poles should be less than zero. The functions d2c, step, polyfit in MATLAB can be used in this part for inertia constant extraction.
The proposed RoCoF calculation method is also applicable for the event data of real power systems. In other words, the inertia constant can be estimated from the event data based on the proposed RoCoF calculation method and the swing equation. In this case, the event data can be treated as a complementary of ambient data.

D. INERTIA CONSTANT TRACKING
In this part, sliding window 2 method and exponential smoothing method are used to update the estimated results in realtime, thus realizing the online tracking of the inertia constant. 1,2 The sliding windows in these two parts are different. The former is for RoCoF calculation, and the latter is for inertia constant tracking.
For the sliding window method, a fixed size time window is used to estimate the inertia constant, and the window is gradually updated for the next estimation. The performance of the sliding window method is dependent on both the choice of sliding window length and the estimation refresh rate.
Generally, the estimation accuracy increases with sliding window length. However, no estimates are available during the first sliding window, and the time delay of the first estimation result can be reduced with a shorter sliding window length. The estimation refresh rate determines how fast a new estimation should be done. Though a faster refresh rate makes the estimator more responsive for real-time estimation as well as tracks change of inertia constant better, it increases computational burden, thus leading to a longer execution time.
Though the inertia constant extracted from the identified models can be updated in real-time, there may be a few low FR models that introduce inaccurate estimates, leading to large fluctuations of the inertia constant. To smooth the estimates from the low FR models, the exponential smoothing method was used, whose original form is defined as: (16) where ν t and ν t−1 are the smoothing value at time t and time t − 1, respectively; θ t is the actual value at time t; and α is the smoothing constant, ranging from 0 to 1. In addition, the inertia constant extracted from the unstable identified model is far from the real value, which should be detected and removed. In this paper, we replace the abnormal estimates from the unstable identified model at time t with the estimates at time t − 1.
Combining the methods mentioned above, the inertia constant can be updated as follows: if unstable (17) where H t and H t−1 are the smoothing value of inertia constant at time t and time t − 1, respectively; h t is the inertia constant extracted from the identified model at time t; µ t is the fitting coefficient at time t, which is equal to FR divided by 100; and k is the exponent, which can be set from 20 to 80 to get a good smoothing effect.

IV. CASE STUDY A. SYSTEM BACKGROUND
The proposed method was tested in the IEEE 39-bus system, which is a simplified model of high voltage transmission system in the northeast of the USA (New England area). The system consists of 10 generators, 39 buses, 19 loads, 34 lines and 12 transformers. The rated frequency is 60 Hz and the main voltage level is 345 kV. Generator 1 is the equivalent generator for the external power grid, and generator 2 is the balance generator. Automatic voltage regulators and governors are used for generator 2-10. The time domain simulations are carried out using Digsilent/ Powerfactory software. The single line diagram of the test system is shown in Fig. 3.

B. ONLINE INERTIA CONSTANT ESTIMATION
In this part, the inertia constant estimation accuracy of the proposed method is validated in the widely used IEEE 39-bus system mentioned above. All 19 loads in the 39-bus system are injected with Gaussian noise filtered by a low-pass filter with a cut-off frequency of 5 Hz, thus simulating the small random load fluctuations and other related small variations of real power systems during normal operating conditions. Then, the inertia constant of all the individual generators as well as the equivalent inertia constant of both the areas and the whole system is estimated using the methods introduced in section III-A to C.

1) ONLINE INERTIA CONSTANT ESTIMATION OF THE INDIVIDUAL GENERATOR
The experiment of a single estimation can be carried out with the following procedures: Step 1: The bus frequency f i at the generator connection bus and the active power output P i are measured by PMUs. The length of the measurements is 200 s and the sampling rate is 100 Hz. The first 100 s signals are used for model identification, and the last 100 s signals are used for model cross-validation.
Step 2: Data preprocessing is employed to the signals. The bus frequency f i and the active power output P i are converted into per unit values by dividing their base values 60 Hz and 100 MW, respectively; the trend is removed; the high-frequency components are filtered using a noncausal low-pass filter with a cut-off frequency of 0.5 Hz; and the sampling rate is decreased from 100 Hz to 5 Hz to avoid numerical unstable problems. Fig. 4 shows the preprocessed active power and bus frequency of generator 2. The fluctuations of both the active power and the bus frequency under normal operating conditions is rather small. The preprocessed measurements of other VOLUME 8, 2020 generators are not provided here since they are similar to that in Fig. 4.
Step 3: With the processed active power P i and bus frequency f i , the state space model with P i as input and f i as output is identified using N4SID. Then the model cross-validation is carried out to evaluate the reliability of the identified model.  It indicates a reliable model is identified as the validated output fits the original output well.
Step 4: The inertia constant of the generator is extracted from the identified model using the method proposed in section III-C. The identified model is disturbed with a unit step signal; the unit step response of the identified model is collected; the initial RoCoF is calculated from the unit step response of the identified model; and the inertia constant is extracted using (15). Fig. 6 shows the unit step response of the identified model of generator 2 as well as the calculation of initial RoCoF  during inertial response, where the first 0.5 s of the result is enlarged. Table 1 shows the detailed results of the inertia constant estimation with a rated apparent power (S b ) of 100 MW. In Table 1, the estimated inertia constant (H est ) is very close to the real inertia constant (H ref ), and the largest relative error for the estimated inertia constant (H est ) is below 5%, thus the feasibility of the proposed method to estimate inertia constant of individual generator is verified.
Generally, the results from a single estimation are uncertain. To evaluate the accuracy of the proposed method and to further reduce the relative error by averaging the results from multiple estimations, the sliding window method is used to obtain various estimated results from one measurement period. In our simulations, a 10-minute measurement with a window length of 200 s and a refresh rate of 1 s is used; in total, 401 inertia constant estimations are carried out. For each estimation, steps 1-4 are employed sequentially; then, all the estimated results are obtained. After all of the estimations are completed, the results from the identified model with an FR above 95% are selected for statistical analysis, which means the point estimation and the 95% confidence interval (CI) estimation are carried out. Fig. 7 shows the normal distribution behavior of the estimated inertia constant (H est ) of generator 2. It can see that 106 inertia constant estimates are selected for statistical analysis. The distribution of the estimated inertia constant from other generators is not provided here due to the limited space. Table 2 shows the statistical analysis results of the estimated inertia constant. In Table 2, we can see that the average of the

2) ONLINE EQUIVALENT INERTIA CONSTANT ESTIMATION OF THE AREA AND THE WHOLE SYSTEM
Considering the coherence of the generators after a large disturbance, the 39-bus system can be divided into 4 areas, as shown in Fig. 2 [32]. Then, the equivalent inertia constant of the area and the whole system is estimated following steps 1 to 4. Note that the equivalent inertia constant of the area and the whole system is estimated by summing the active power outputs and aggregating the bus frequencies.
To evaluate the accuracy of the developed method and to further reduce the relative error, the sliding window method is used to obtain various estimates from one measurement period. In our simulations, a 10-minute measurement with a window length of 200 s and a refresh rate of 1 s is used; in total, 401 inertia constant estimations are carried out. In each estimation, steps 1 to 4 mentioned above are employed sequentially. After all of the estimations are completed, the results from the identified model with an FR above 95% are selected for statistical analysis.   Table 3 shows the statistical analysis results for the estimations. In Table 3, we can find out that the average of the estimated inertia constant (H avg ) is nearly the same as the real inertia con-

3) COMPARISON WITH THE DISTURBANCE DATA METHOD
The disturbance data method in [13] is carried out to compare with the method proposed in this paper.
Step 1: The active power of Load 15 is increased from 320 MW to 1820MW suddenly, resulting in an active power deficiency in the system. In this case, the active power imbalance P is equal to 15 (in p.u.).
Step 2: The transient frequency following this event is measured. A five-order low pass Butterworth filter with a cutoff frequency of 0.5 Hz is used to isolate the dominant system inertial response. The processed transient frequency is shown in Fig. 9.  Step 3: The RoCoF is calculated using a 500-ms sampleby-sample sliding window, over a 2 s period following the event. The maximum value is then taken to represent the RoCoF following the event before the primary frequency response starts to take effect.
Step 4: With the method in [13], the equivalent inertia constant of the whole system is estimated as 809.98 s. The estimate value is close to the real value of 782.7 s and the relative error is 3.49%.
Besides, the relative error of the estimated inertia constant from the disturbance data method and that from the method proposed in this paper are both below 5%, thus verifying the effectiveness of the proposed method further.

C. ONLINE INERTIA CONSTANT TRACKING
In this part, the inertia constant online tracking capability of the proposed method is validated in the IEEE 39-bus system. Here, each synchronous generator of the 39-bus system is regarded as an equivalent generator, which is aggregated by all the generators of a certain area. The inertia constant of the synchronous generator of the 39-bus system is changed to simulate the time-varying equivalent inertia constant caused by the switching of the generators in the area. Considering a measurement length of 120 min, the inertia constant of all generators remains at H 1 during the first 40 min, then changes to H 2 and remains there from 40 to 80 min. The inertia constant stays at H 3 during the last 40 min. H 1 , H 2 and H 3 can be regarded as the inertia constant vectors here that represent the inertia constant of all generators in different periods.  Table 4 shows the inertia constant of all generators in different periods.
All the loads in IEEE 39-bus system are modeled to be random, the same as section IV-B. Then, the inertia constant of all generators in the 39-bus system is tracked using the methods introduced in section III-A to D. The active power output P i and bus frequency f i with a 120-min length are measured by PMUs. The sliding window method with a window length of 100 s and a refresh rate of 1 s is employed. For each sliding window, steps 2-4 in section IV-B-1) are employed sequentially and then, the estimates are smoothed using the exponential smoothing method in (17). Once the measurements from the sliding window are obtained, the inertia constant will be estimated immediately. The execution time of a single estimation is less than 1 s, which means that the inertia constant can be updated before new samples are completely collected, leading to real-time online tracking of the inertia constant. Inertia constant tracking trajectories of generator 7 and generator 8 are shown in Fig. 10 and Fig. 11, respectively. The inertia constant tracking trajectories of other generators are not provided here for the limited space. Both Fig. 10 and Fig. 11 show that the inertia constant tracking trajectories fluctuate around the reference trajectories in a small range, meaning that the inertia constant can be tracked accurately in real-time using the proposed method.
According to section III-D, the execution time of the procedures depends on the choice of sliding window length and estimation refresh rate. Table 5 shows the execution time for different settings of sliding window length and refresh rate with a 120-min measurement. The simulations are carried out in a regular office laptop with Inter(R) Core(TM) i7-7700HQ and CPU at 2.80 GHz. The results of the first four settings show that the execution time is inversely proportional to the refresh rate as fewer estimations are done with a slower refresh rate. The execution time of the fifth setting is comparable to the execution time of the first setting, indicating that a larger sliding window length will increase the computational burden. Additionally, processing one measurement period takes approximately 6.22-31.36% of its length using the laptop introduced above. On average, it takes 3.73-18.81 s to analyze 1 min of the measurements. In addition, a sliding window length of 100-200 s with a refresh rate of 1-2 s is responsive enough to obtain reliably accurate results in our simulations.
For various power systems, the sliding window length can be adjusted properly until an acceptable dynamic model can be identified and the faster refresh rated can be determined to track the inertia constant better using a more powerful computer. All the analyses above confirm that the proposed method is suitable for power system inertia constant real-time online tracking.

V. CONCLUSION
This paper proposed an online estimation method for power system inertia constant under normal operating conditions. First, the dynamic model between active power output and bus frequency measured by PMUs was identified using the subspace identification method. Then, the inertia constant was extracted from the unit step response of the identified model in the time domain. Finally, the sliding window method and the exponential smoothing method were used to update the inertia constant in real-time.
The proposed method was tested in the IEEE 39-bus system. The results confirmed that different hierarchical power system inertia constant (individual generator, area and the whole system) could be estimated with high accuracy using ambient data, which could overcome the disadvantages of the disturbance data methods. Moreover, by using the proposed method, the inertia constant could be updated on a time scale of seconds, and the inertia constant tracking trajectories could also be provided in real-time.