© The Author(s), published by EDP Sciences, 2024

1 Introduction

With the large-scale development and grid-connection of new energy such as photovoltaic and wind power, the permeability of new energy in electric power systems has been continuously improved in recent years [1]. In view of this, a series of studies related to system inertia have been carried out. Among them, the quantification of the inertia level of the power system is the premise and foundation for solving this problem. It provides important guiding significance for system planning, operation scheduling, and equipment control, and so on [2, 3].

There are three types of disturbances on which the existing inertia evaluation methods are based [4]: frequency events, small disturbance events, and quasi-steady-state quasi-steady state operation. Frequency events refer to events with large frequency deviations that occur in the system. Small disturbance events refer to small amplitude disturbances or disturbances in the system, which usually occur when the system is operating in a steady state. Quasi-steady-state operation refers to the operation of the power system under relatively stable working conditions. Through exciting the transient characteristics of the system, the rotor motion equation [5–9] and electromechanical wave theory [8–11] are used for analyzing and processing the measured data to obtain the system inertia in an evaluation method based on frequency events or small disturbances. However, the above method requires external disturbances, which not only fails to track the inertia in real time but also has adverse effects on the safety and stability of the power system. In order to get rid of the above limitations, a statistical model for real-time estimating system inertia based on quasi-steady-state frequency changes is proposed in [5]. Historical data from two consecutive years is used for training to obtain a switching Markov–Gaussian model that can reflect the correlation between steady-state frequency changes and system inertia. Then the system inertia level can be obtained online through measured data. Based on the system equivalent aggregation model, the dynamic regression expansion and hybrid method are used to solve the model parameters in [6]. In [7], a statistical method of power system inertia based on network recursive observation data is proposed. The dynamic model is used to describe the inertial random process in a feasible and adaptive way and the maximum likelihood method is used to perform statistical inference on the model, so as to solve the dynamic change of the inertia of the power system.

The popularization and application of wide-area measurement devices provide favorable conditions for solving the problem of online inertia evaluation. During normal operation of the power system, there are always small disturbances caused by normal load fluctuations, power flow control, and transformer tap adjustment, resulting in random fluctuations similar to noise in state variables such as power, frequency, and voltage of the system nodes. This kind of signal contains rich information on the dynamic characteristics of the power system and can be continuously collected.

At present, identification is the most important method to determine the system parameters. Compared with the mechanism modeling method, the model can be obtained by observed data, and the actual system can be simplified without making assumptions. Therefore, the normalized and continuous tracking of the system inertia constant can be realized by using the identification method to process the collected system noise-like signals.

Aiming at the problems of quantitative evaluation and fast-tracking of inertia in new energy power systems, the inertia constant tracking method based on system identification is studied. This method is divided into two categories: non-recursive algorithm and recursive algorithm. Through analyzing and processing PMU and other measurement data and system identification, the dynamic model of the generator is established to realize the fast identification of the inertia constant of the generator, region, and system level. Theoretical analysis and experimental results show that the recursive algorithm has high identification accuracy, stable identification results, and fast identification speed. It is suitable for online tracking of the inertia constant of new energy power systems. In addition, measurement noise is taken into account, which has more engineering application value.

2 Quantitative evaluation principle of inertia constant of power system

When the generator (or equivalent synchronous machine) is disturbed, the dynamic relationship between power and frequency can be expressed by the swing equation. A single generator i, it can be described as:(1)where Δf_i is the unit value of the rotor frequency deviation, P_m,i is the unit value of mechanical power of generator i, P_e,i is the electromagnetic power of generator i, H_i and D_i are the inertia constant of the generator and damping coefficient of generator i respectively.

Similarly, for the region j in the power system, it is equivalent to the analysis of the equivalent synchronous machine, and the equation is expressed as:(2)where P_m,i, P_L,i and P_line are the per unit of the equivalent mechanical power, the active power of the internal load and the output power of the tie line in region j respectively, H_j is the inertia constant of the regional equivalent machine j, D_j is the damping coefficient and f_COI,j is the aggregated frequency per unit value of region j, which is obtained by aggregating all generator frequencies in the region.

There is a normal small amplitude fluctuation in the active power frequency in the ambient excitation. Due to the limitation of the dead zone, the primary frequency modulation does not act and the mechanical power of the prime mover remains unchanged. The dynamic relationship between power and frequency of generator i and region j can be expressed by:(3) (4)

For region j, the sum of the active power of the internal load and the output power of the tie line can be regarded as the total output power of this region. Therefore, (3) and (4) can be unified as:(5)where Δf is the per unit of deviation of generator terminal bus frequency or regional equivalent aggregation frequency, ΔP is the per unit of deviation of generator output power or regional equivalent total output power.

The Laplace transform of (5) is obtained as:(6)

The inverse Laplace transform of (6) is carried out, and then the initial impulse response value is obtained:(7)

In the noise-like environment, when a large disturbance occurs in the system, the frequency of the generator (or equivalent synchronous machine) will deviate from the rated value and exceed the dead zone limit, which triggers the primary frequency modulation. At this time, the active frequency transfer function includes the inertial response and the primary frequency modulation response process G_p(s). It can be expressed as:(8)

Similarly, the inverse Laplace transform of (8) is carried out, and the initial impulse response value is solved:(9)

It can be seen from (7) and (9) that in the noise-like environment, no matter if the system has a large disturbance, the inertial constant implied in the time series of power fluctuation and frequency deviation can be determined, and it is the negative reciprocal of twice the initial impulse response, namely H =−1/2g(0₊).

3 Inertia constant tracking based on system identification

The system identification method is applied to the online tracking of power system inertia, which can be divided into the following two parts: system identification [8] and inertia constant extraction and tracking.

3.1 System identification

There are generally three steps of system identification [9]: data preprocessing, model structure selection, and model parameter identification.

3.1.1 Data preprocessing

In order to ensure the reliability of the identification results, it is necessary to preprocess the generator output power and its terminal bus frequency signal, mainly including normalization, de-mean (or de-trend), and filter processing. Considering that the inertia constant of the generator is generally longer than 2 s, the cut-off frequency of the low-pass filter f_c is set to 0.5 Hz.

3.1.2 Model structure selection

Considering the noise and external random interference in the actual environment, this paper adopts a time series model: Autoregressive Moving Average with Extra Input, ARMAX. As a typical model of the linear system, it has the advantages of being simple, fast, accurate and robust [10, 11].

The relationship between the input signal u(k) and output signal y(k) of the model ARMAX can be expressed as:(10)where , , , ε(k) is zero mean white noise of the system, n_a, n_b, n_c are the order of each delay polynomial, n_k is input–output delay, z − 1 is tge retardation factor.

After the model type is determined, it is also necessary to determine the order of each delay polynomial in the model. If the order is too small, it may not be enough to reflect the main dynamic characteristics of the system. If the order is too large, the complexity of the model will also increase, and even cause over-fitting problems. Therefore, the choice of model order needs to seek the best balance between model complexity and the ability of data set description. Akaike Information Criterion (AIC) is a standard to measure the goodness of fit of statistical models [13]. Models with small AIC values are generally selected. In the case of ARMAX model, the AIC function is defined as:(11)where N is the sample number, is model variance, k = (n_a + n_b + n_c) is the number of independent parameters of the model.

A set of model classes can be obtained by setting a lower-order range. By comparing the fitting degree of the simulation output and the actual measurement output of each model, the model with a higher fitting degree is selected as the reliable model class, and then the AIC is used to select the optimal model among the model class.

3.1.3 Identification criterion and identification algorithm

The identification criterion is used to measure how close the model is to the actual system, which is generally expressed as the error function of the model and the actual system. The formula is determined as:(12)

Identification algorithms mainly include classical system identification algorithms and modern system identification algorithms. Among them, the least squares method is a classical identification algorithm with low computational cost, including two basic structures: one is a non-recursive algorithm, which is characterized by batch processing using a batch of data to obtain the estimated value of the identification model parameters; the other is the recursive algorithm. Based on the estimated value of the model parameters at the previous moment, the estimated value is corrected according to the fresh data. As time goes by, the identification results will be continuously corrected until the identification error requirements are met.

• 1.

  Non-recursive identification algorithm based on the ARMAX model

The non-recursive algorithm of least squares is used to solve the parameter values of ARMAX model, and (10) is rewritten as:(13)where h(k) is the data vector, θ is the parameter vector.(14)

According to (12) and (12), the criterion function is obtained as:(15)

According to the extremum principle, (15) is minimized, and the parameter vector of the optimal model needs to satisfy the following equalities:(16)

It can be seen from (16) that if the equation holds, the criterion function J(θ) takes the minimum value.

• 2.

  Recursive identification algorithm based on ARMAX model.

From , the parameter estimation at time k is determined as:(17)

If and , then:(18)

Defining that and , using the inversion formula, the least squares recursive expression of the ARMAX model is obtained as:(19)

Based on and , the first term of (17) can be expressed as:(20)

In order to make (20) approximate (17), P ⁻¹(0) approach 0, and let approach 0, the initial value is defined as:(21)where α is a sufficiently large real number, ε is a sufficiently small real vector.

As the sample data increases, the gain matrix K and the correction matrix P gradually approach zero. The correction effect of the fresh data on the parameter estimation is getting smaller and smaller, and the phenomenon of “data saturation” appears [14], resulting in large errors. Therefore, the fading memory method is used to weight the criterion function, and a forgetting factor is added to the old data to highlight the amount of information provided by the new data. The least squares recursive expression using the fading memory method is expressed as:(22)

A smaller value of λ indicates that the old data has less influence on the current parameter estimation. However, due to the existence of noise interference, when the forgetting factor is too small, the variance of parameter estimation will increase, so the general value range of λ is 0.9–1. The recursive flow chart is shown in Figure 1.

Fig. 1 Recursive algorithm flow chart.

In order to further illustrate the difference between the non-recursive algorithm and the recursive algorithm in principle, the DIGSILENT (or Powerfactory) simulation software is used to simulate and analyze the IEEE39 node system. A 100 s sample data of the No. 6 generator is collected and analyzed in a noise-like environment. A total of 500 sets of input and output data are collected. The model order was set to (n_a, n_b, n_c) = (2, 2, 3), and the delay between input and output n_k is 0.

Figure 2 is the comparison between the output results of the two identification algorithms and the measured frequency of the system. It can be seen that the fitting degree of the identification output of the recursive algorithm identification model is higher than that of the non-recursive algorithm. The reason is that the identification of the non-recursive algorithm needs to take into account the fitting effect of the entire data segment. When the model order is limited, the identification model can only fit some of the characteristics of the actual system. However, for the recursive algorithm, when the iteration reaches about 2 s, the output value of the identification model is close to the actual output value, and the parameters of the identification model are constantly approaching the true value from the initial value , as shown in Figure 3.

Fig. 2 Comparison of output results of two identification algorithms with system-measured frequency.

Fig. 3 Iterative curve of parameter vector .

3.2 Extraction and tracking of inertial constant

3.2.1 Extraction of inertial constant

The ARMAX model reflecting the main characteristics of the original generator is obtained by system identification, and the inertia constant is extracted from it. In order to avoid the problem of error accumulation caused by model reduction, the time domain impulse response curve of G(s) is directly obtained, and the minimum value in the range of 0–2 s of the curve is taken as the maximum frequency change point, and then the H_i is obtained by substituting it into (7) and (9). Therefore, the extraction steps of inertia constant can be summarized as follows: (1) The bilinear variation method is used to construct the mapping relationship between the z-plane and the s-plane, and the identification model is transformed from the discrete difference equation form G(z) to the continuous high-order transfer function form G(s). (2) The time domain impulse response curve of G(s) is obtained. (3) The minimum value of G(s) time domain impulse response curve in the range of 0–2 s is used as the initial frequency change rate to solve H_i.

3.2.2 Tracking of inertia constant

In order to monitor the inertia level of the power system in real-time, it is necessary to continuously track the extracted inertia constant H_i. The schematic diagram of the inertia constant tracking principle is shown in Figure 4. For the non-recursive algorithm, the real-time update of H_i can be realized by moving the data window method. The H_i(t) of the time period is obtained in each data window [t − T_w, t], and then the inertia constant of the next moment H_i (t + Δt) is obtained by moving the data window. For the recursive algorithm, the iterative process of the algorithm continues over time, so it is only necessary to obtain the corresponding identification model at each iteration and then extract H_i(k) to complete the real-time update of H_i. The process of the power system inertia constant tracking method based on system identification is shown in Figure 5.

Fig. 4 Inertial constant tracking diagram.

Fig. 5 Block diagram of power system inertial constant tracking method based on system identification.

In summary, the difference between the two identification algorithms is as follows: 1) In terms of update speed, the recursive algorithm can complete an iteration and H_i update every interval of a data sample, and its update interval is shorter, generally not exceeding 1/f_c = 0.2s; the update speed of the non-recursive algorithm depends on the sliding length Δt of the window. If Δt is set too small, it means that the number of identifications required will increase inversely, resulting in an increase in the overall time-consuming ratio. Therefore, the update interval is generally not less than 1 s; 2) In terms of storage capacity, the non-recursive algorithm needs to store all the data in the data window at that time for each identification. The sample size is f_c*T_w, and the data storage capacity is large, while the recursive algorithm only needs to store the data vector h(k) at that time, and the data storage capacity is small; 3) In terms of computational complexity, if the total length of the data segment is T_total, the number of recursive algorithm operations is f_c*T_total, and the number of non-recursive algorithm operations is (T_total − T_w)/Δt. Although the recursive algorithm has more identification times and the steps of each iteration are more complicated (K(k) and P(k) need to be obtained), the non-recursive algorithm requires sufficient data to ensure its identification accuracy, so the T_w value is generally more than 50 s. With the sliding of the data window, (1 − 2Δt/T_total) × 100% of the data is reused, and the number of reuses is T_w/Δt, so the number of non-recursive algorithm operations is relatively more. In addition, the H_L dimension in the non-recursive algorithm is larger, and the matrix inversion operation is larger. Therefore, the overall computational complexity of the recursive algorithm is much smaller than that of the non-recursive algorithm, and the speed of identification is faster.

4 Simulation example analysis

In order to verify the correctness of the theoretical analysis, the time domain simulation of the IEEE39-node system is carried out by using DIGSILENT/Powerfactory simulation software. The system wiring diagram is shown in Figure 6.

Fig. 6 IEEE39 node system connection diagram.

In order to verify that the method has good adaptability under different load disturbance types and different new energy penetration scenarios, four sets of real experiments were designed to simulate the system operation scenarios under different working conditions.

Experiment I (random small disturbance experiment): Small amplitude Gaussian white noise is superimposed on the 19 loads of the system to simulate the random load fluctuation in the actual power system. The simulation time is 5 min, and the non-recursive algorithm is used. Data window length T_w = 50 s, window displacement Δt = 1 s; recursive algorithm forgetting factor λ = 1.

Experiment II (large disturbance experiment): Based on the noise environment of Experiment I, the simulation time was extended to 10 min. When t₁ = 200 s, the load of node 29 increases by 150 MW for 1 s, and when t₂ = 400 s, the load of node 29 decreases by 150 MW for 2s. The non-recursive algorithm data window length T_w = 100 s, window displacement Δt = 1 s, recursive algorithm forgetting factor λ = 0.98.

Experiment III (new energy permeability change experiment): The experimental setup is the same as Experiment I. The new energy permeability in the simulation scene is changed, and some synchronous generators in the system are replaced with wind turbines or photovoltaic power supplies without virtual inertia control.

Experiment IV (online tracking experiment): In order to verify the online tracking ability of the method, the inertia constants of No. 1, No. 2, No. 4, and No. 8 generators in the system are modified. The parameters are shown in Table 1, and the parameters of the remaining generators remain unchanged. Based on the noise environment of Experiment 1, the simulation time is extended to 60 min, and the forgetting factor is 0.99.

4.1 Random small perturbation experiment

For the non-recursive algorithm, the single identification is prone to contingency and the error is large. The moving data window method can be used to identify and extract 251 times of 5 min measured data independently, and the results are mathematically counted to reduce the error. According to the coherence of the generator after the disturbance, the system is partitioned as shown in Figure 7. The obtained generator inertia H_i is further calculated by frequency aggregation and power accumulation to obtain the inertia constant of the region/system. Figure 7 is the box plot of the multiple identification results of each generator and region/system. Figure 8 is the distribution of the identification value of the inertia constant of the No. 10 generator, and the identification results of the other generators are similar, which is not repeated here.

Fig. 7 Identification results of inertia constant of each generator and each area/system.

Fig. 8 Normal distribution of inertial constant identification value for generator No. 10.

It can be seen that the inertia constant identification results obtained by multiple identification extraction are basically normal distribution, and some generator identification values have a small amount of abnormal data, which are eliminated by the absolute median deviation algorithm. Similarly, the identification data segment that gradually tends to be stable 3 min after the recursive algorithm is intercepted is processed in the same way, and the identification results are shown in Table 2.

The error range of the non-recursive algorithm is ±0.90%–9.91%, and the standard error is 0.029–0.316. The error range of the recursive algorithm is ±0.03%–8.93%, and the standard error is 0.005–0.086. The error range and standard error are smaller than the former. Therefore, the inertial constant identification results using the recursive algorithm are more accurate and stable.

4.2 Large perturbation experiment

In Experiment 2, the identification results of the inertial constants of 10 generators under the two algorithms are obtained, as shown in Figure 9.

Fig. 9 No. 10. Identification curve of inertial constant under large disturbance. (a) Hi identification curve based on non-recursive algorithm. (b) Hi identification curve based on recursive algorithm.

For the non-recursive algorithm, if the number of samples in the first 100 s does not constitute a complete data window, there is no identification result. When two large disturbances occur in the system, the identification results fluctuate obviously, and the fluctuation is more obvious when T_w = 100 s after the large disturbance. The reason is that the amplitude of input and output data jump caused by large disturbance is larger. As the data window moves, it appears at both ends of the data window, and the change of identification parameters caused by it is more obvious than when the data jump appears in the middle of the data window. Two large disturbances were applied in the experiment, and the identification results of the non-recursive algorithm showed four large fluctuations. Table 3 shows the error analysis of the maximum identification results of the two fluctuation amplitudes. It can be seen that the maximum fluctuation error is 10.29%, and the results are more secure. For the recursive algorithm, when the system has a large disturbance, there are two obvious fluctuations in the identification value, and then with the iteration of the algorithm, it quickly returns to the true value.

4.3 New energy penetration rate change experiment

The overall inertia level of the system decreases with the increasing penetration rate of new energy sources. The results of system inertia constant identification under different new energy penetration scenarios are shown in Table 4. It is clear that both algorithms have high discrimination accuracy and can realize the quantitative assessment of the inertia level of the high percentage new energy power system.

4.4 Online tracking experiment

Figure 10 shows the online tracking results of the inertia constants of generator 1, generator 5, and generator 8, in which the parameters of generator 5 are kept unchanged as a control group, and the tracking result is a curve with small fluctuation amplitude. The tracking results of the rest of the generators are similar, and the tracking results can effectively track the changes in the actual reference curve, which will not be repeated here.

Fig. 10 Online tracking result of inertia constant.

Fig. 11 Tracking results of Hi for generators No. 8 with different.

As can be seen from Figure 10, the inertia constant tracking methods based on two different algorithms have a good tracking effect. When the system parameters do not change, the identification value fluctuates slightly around the actual value; when the system parameters change suddenly, both algorithms can realize the fast-tracking of inertia constants. Although the two algorithms do not differ much in tracking numerical results, there are large differences in the following aspects:

1. Storage Volume. The non-recursive algorithms are independent of each other for each identification, and the amount of data storage required for each identification is f_c * T_w = 250 sets; whereas the recursive identification only needs to store the data vector h(k) at that moment in time, and the amount of storage is equal to the total order of the ARMAX model, namely, (n_a + n_b + n_c). The model orders of different generators in this algorithm are slightly different, but the total order is generally not more than 18, so the storage cost of executing the recursive algorithm is very low, about 7.2% of the non-recursive algorithm, which is helpful to reduce the communication cost of data transmission and improve the overall operation speed;

2. Computation Volume. A total of (3600 – T_w)/Δt = 3550 identifications are performed in this algorithm, of which 97.2% of the frequency-active data are processed times, for a total of 8.875 × 10⁵ sets of operations; the recursive algorithm iteratively per-forms 18,000 recursive identifications, with the amount of data processed in each iteration being equal to the number of dimensions of the parameter vector h(k), namely, (n_a + n_b + n_c), and the total number of operations does not exceed 3.24 × 10⁵ groups, which is 36.5% of the non-recursive algorithm. In addition, the H_L dimension is larger in the non-recursive algorithm, and the computation volume increases sharply in the matrix inverse operation. Therefore, the recursive algorithm is smaller than the non-recursive algorithm, which is conducive to reducing the computational energy consumption shortening the computation cycle, and further accelerating the overall identification speed;

3. Influencing factors. With the sampling frequency f_c and the model order (n_a, n_b, n_c) unchanged, the identification accuracy of the non-recursive algorithm is related to the sample size of the batch, namely, to the length of the data window T_w: the larger T_w is, the more stable the identification results are. Accordingly, when the system parameters change, the longer the tracking transition time will be, and the recognition curve changes slowly and uniformly, which cannot realize the fast-tracking of inertia constants. For the recursive algorithm, its tracking effect depends on the forgetting factor λ: the larger λ is, the slower the old data is forgotten, the smoother the tracking curve is, and correspondingly the longer the tracking transition time is, as shown in Figure 10.

In order to compare the computing speed of the two algorithms under different parameter configurations, two computers are used to perform MATLAB calculations: device 1 processor is Inter Core i5-5200U, 2.20 GHz, with 4 G of RAM, and device 2 processor is AMD Ryzen 70.1’5800H, 3.20 GHz, with 16 G of RAM, and the results of the computing speed comparison are shown in Table 5.

1. The running speed of the non-recursive algorithm is mainly related to the update step Δt, with less correlation with the data window length T_w; when T_w is of the same order of magnitude, the execution time of the two algorithms is not much different;

2. The running speed of the recursive algorithm is mainly determined by Δt, which is numerically equal to the inverse of the sampling frequency f_c, and it is about 0.1–0.2 s;

3. After replacing the device with a better hardware condition 2, the non-recursive algorithm runs about 4 times faster and the recursive algorithm runs about 8–10 times faster. The running speed of the recursive algorithm is much higher than that of the non-recursive algorithm with a lower update speed, and the improvement of hardware conditions is more obvious for the recursive algorithm. Therefore, in practical engineering applications, the algorithm language and hardware conditions (processor, CPU, memory, etc.) can be improved to further accelerate the tracking accuracy and running speed of the inertia identification algorithm.

5 Conclusion

Aiming at the inertia online tracking of the new energy power system, this paper realizes the online tracking of the inertia constant of the power system through the identification of the model system and the extraction and tracking of the inertia constant. Moreover, the recursive algorithm and non-recursive algorithm are compared and analyzed from the principles and implementation methods. The results show that the recursive algorithm has high identification accuracy, stable identification results, and lower storage and operation costs. It can effectively improve the identification and tracking speed of the algorithm and is suitable for real-time monitoring and rapid sensing of inertia of new energy power systems with large numbers of nodes and complex structures.

Conflict of interest

The authors declare no conflict of interest.

Acknowledgments

This work was supported by the Science and Technology Project of Hebei Electric Power Company (kj2022-059).

References

All Tables

All Figures

	Fig. 1 Recursive algorithm flow chart.
In the text

	Fig. 2 Comparison of output results of two identification algorithms with system-measured frequency.
In the text

	Fig. 3 Iterative curve of parameter vector .
In the text

	Fig. 4 Inertial constant tracking diagram.
In the text

	Fig. 5 Block diagram of power system inertial constant tracking method based on system identification.
In the text

	Fig. 6 IEEE39 node system connection diagram.
In the text

	Fig. 7 Identification results of inertia constant of each generator and each area/system.
In the text

	Fig. 8 Normal distribution of inertial constant identification value for generator No. 10.
In the text

	Fig. 9 No. 10. Identification curve of inertial constant under large disturbance. (a) Hi identification curve based on non-recursive algorithm. (b) Hi identification curve based on recursive algorithm.
In the text

	Fig. 10 Online tracking result of inertia constant.
In the text

	Fig. 11 Tracking results of Hi for generators No. 8 with different.
In the text