Analytic assessment of the power system frequency security

Recently, with the increased intermittent renewable energy penetration, many power grids have been incorporated into a large-scale long-distance UHV AC–DC hybrid system. In the actual grid, power disturbances, such as DC blocking faults and trip-off of wind turbines, often occur, resulting in power shortages and large frequency ﬂuctuations. However, the existing approaches to assess the system frequency security are not reliable. This paper proposes an analytic formula for the system frequency response based on a generic system frequency-response (SFR) model, which can be applied to modern large-scale power systems. First, a generic SFR model with a reasonable structure is designed according to the parameter determination strategy. Second, according to the transfer function in the model, the time-domain analytic formula of the frequency response is obtained by realizing the inverse Laplace transform. Moreover, four main indexes are established to represent the characteristics of the frequency dynamic process in different periods, and these indexes are qualitatively and quantitatively analysed. Finally, the New England 10-unit 39-bus power system and East China Power Grid are considered to demonstrate the key features of the proposed method. The results show that the proposed analytic assessment method can be effectively adopted in several applications.


INTRODUCTION
Frequency, as an essential dynamic parameter of power systems, considerably influences the power grid and users [1]. The frequency of a large-scale power grid depends on the realtime balance between power and load. Traditionally, the power system is integrated with thermal and hydro generation systems and subjected to strict frequency regulation; in this scenario, the frequency stability can be ensured. However, with the increased intermittent renewable energy penetration, the mechanical inertia of power systems has gradually decreased, which poses a threat to the system frequency stability [2]. Additionally, many power grids have gradually been incorporated to a large-scale long-distance UHV AC-DC hybrid system, leading to an increased possibility of frequency deviations caused by abnormal operation [3][4][5]. With the accelerating deployment of smart grid technologies, the traditional control methods of a power system have been enhanced by using advanced information communication technologies [6]. The precise estimation of frequency drops is the main concern for frequency control resources involved in addressing the frequency stability issue at the second or millisecond level [7]. However, the existing analysis methods for the frequency dynamic response are not sufficiently accurate and cannot ensure the effectiveness of the control measures. Therefore, establishing an approach to realize the analytic assessment for the modern large-scale power system frequency security is of significance. The accurate assessment of the system frequency security is based on the precise modelling of the system frequency response. The existing modelling methods for the system frequency response (SFR) mainly include full model time-domain simulations, linearized models, artificial intelligence techniques, and single-machine equivalent models. At present, the full model time-domain simulation method is the most widely used technique to calculate the dynamic frequency owing to its clear principles and comprehensive considerations. However, this approach is not suitable to realize online analyses owing to the large number of parameters involved and complex parameter setting. The linearized model is based on the full model time-domain simulation method with partial linearization of the model; however, this approach involves the abovementioned limitations as well [8,9]. The artificial intelligence method depends on a large amount of measured data, and thus, its popularization and application in actual power grids are limited FIGURE 1 Classic SFR model [10,11]. The single machine equivalent model can obtain the analytical solution of the frequency dynamic response and is widely applied to assess the power system frequency security [12][13][14][15][16].
The equivalent models mainly include the average systemfrequency model (ASF) [17], and SFR model [18]. The classical SFR model proposed by P. M. Anderson is based on the single-machine equivalent model and simplified reheated steam turbine-governor model. Therefore, this approach cannot be applied to modern power systems involving hydro and renewable generation. Several studies have attempted to overcome the limitation of the classical SFR model. In [19], an analytical method is proposed to aggregate the multi-machine SFR model as a single-machine model. In [20], a modified SFR model considering the available inertial and droop responses from wind farms is described. However, these two enhanced variants of the SFR model are not universal and cannot be applied to modern power systems. In [21], a generic SFR model with a reasonable structure and parameter-determination strategy is presented, which can achieve a high accuracy in lower order frameworks. The generic SFR model in [21] can well reflect the frequency dynamic response, which provides the basis for this paper.
The frequency security is usually represented in terms of certain major indexes, such as the initial rate of change, extreme value, and steady-state value. Therefore, to realize analytic assessments, it is of significance to determine the value and characteristics of the frequency security indexes. However, recent studies on frequency security are mainly limited to the use of the classical equivalent model, owing to which, the analytic assessment lacks sufficient accuracy and is unsuitable for modern power systems. In [22], an adaptive precise system frequency response (SFR) model is developed, based on which the appropriate amount of load rejections is determined to formulate the under-frequency load shedding (UFLS) plans. In [23], an analytical model for the frequency nadir prediction based on the ASF model is described.
To more accurately assess the system frequency security, this paper presents an analytic formula for the system frequency response by using the generic SFR model. Moreover, the qualitative analysis and quantitative calculations are performed for the major indexes of the frequency security. The effectiveness of the analytic assessment of the system frequency security is validated considering the actual disturbance data recorded in the East China Power Grid. The remainder of this paper is organized as follows. Section 2 presents the analytic formula of the system frequency response. Section 3 describes the analysis of the major indexes of the frequency security. Section 4 describes the validation of the analytic formula and assessment method through case studies. Finally, the concluding remarks are presented in Section 5.

System frequency response model
The SFR model is used to average the machine dynamics in a multiple machine system to be equivalent to that of one machine [18], and the average system frequency is defined as the weighted summation of the machine speeds [17], that is: where f is the average system frequency in per-unit; f k is the frequency of machine k in per-unit; H k is the inertia constant of the machine k in seconds. As the fluctuation of the SFR is not considered in this case, the nonlinearity is not considered in the SFR model. Based on the single-machine equivalent system and simplified reheated steam turbine-governor model, P. M. Anderson and M. Mirheydar proposed a classical SFR model in [18]. The block diagram of the classical SFR model proposed in [18] is shown in Figure 1. ∆f is the frequency deviation; ∆P d is the power disturbance including the generation and load; ∆P m is the mechanical power deviation, which is positive after considering the negative feedback; ∆P a is the accelerating power; 2H is the equivalent inertia constant of the generator, in seconds; D is the equivalent damping factor; R is the droop setting of the governor; km is the mechanical power gain factor, such that km/R is the actual droop coefficient; T R is the reheat time constant, in seconds; and F H is the fraction of the total power generated by the high-pressure turbine.
Currently, various types of power generation frameworks are available. The classical SFR model that only considers reheating turbines cannot satisfy the requirements of modern power systems. The block diagram of the improved model based on the classical SFR model, named the generic SFR model, is shown in Figure 2. This model is not limited to specific types of the prime mover-governor model, and thus, this model can be applied to modern power systems involving thermal, hydro and renewable generation. In [21], a standard transfer function is proposed to describe the equivalent dynamics of the aggregate prime movergovernor.
where a i and b j are the coefficients of the transfer function.
In [21], a parameter estimation method based on the measured data is described, which can uniquely determine all the model parameters. The results of case studies indicate that the system frequency response can be well described by using a second-order transfer function G m (s). Therefore, a three-order transfer function can be considered to describe the system frequency response, which reflects the relationship between the system frequency and power disturbance, as follows: This expression can be rewritten as where A i and B j are coefficients of the system transfer function, which is determined through the parameter estimation based on measured data. These two transfer functions are identical; the first function reflects the physical characteristics, and the second function is simply a more concise form.

System frequency response formula
If a power disturbance occurs in the power grid, such as DC blocking fault, the off-grid and large-scale loads of the largescale wind turbines increase, and the real-time power disturbance ∆P d can be described as: With reference to the Laplace transform, the s-domain function of the frequency response can be expressed as By expanding the partial fraction of Equation (6), the timedomain expression of the frequency deviation ∆f can be obtained by solving the root of the characteristic equation For the characteristic equation of the transfer function of the closed-loop control system, the most common case involving three roots include a real root and a pair of conjugate complex roots. The roots for this equation can be expressed as {s 1 , s 2 , s 3 }: Based on Equation (7), the s-domain function of the frequency response can be expressed as: The coefficients C 0 , C 1 , K 1 and K 1 * (conjugate of K 1 ) in Equation (8) can be defined as follows: If K 1 = x + y j, According to Equations (11)-(13), the time-domain analytic formula of frequency response can be expressed as: As a comparison, the time-domain analytic formula of the frequency response in [18] based on the classical SFR model is described as follows: where Furthermore, Equation (15) can be rewritten as: The proposed analytic formula for the system frequency response Equation (14) contains three terms corresponding to a constant, a monotonic decay and an oscillatory decay, while the classical formula (19) contains only two terms. In contrast to the classical SFR model, the monotonic decay term is added to more accurately describe the frequency dynamic response. This expression effectively reflects the inertial response and damping effect of the system in the initial period of the system frequency response.

Frequency security indexes
Frequency security after a disturbance includes two critical aspects: transient frequency security and steady-state frequency security. In [25], a frequency security index (FSI) is proposed to evaluate the frequency security, which combines five essential characteristics of frequency profile. In [26], the time frame of frequency response is analysed, and the definitions of related indexes are given in each period. Considering both aspects (transient and steady-state) of the frequency security, this paper proposes four main indexes representing various crucial points or periods after a disturbance. According to the time sequence of the dynamic process after a disturbance, the frequency security indexes can be defined as follows: 1. Initial rate of change of frequency, S 0 : The change of frequency in a brief period (in this study, 1 s) after the disturbance is used to measure the change in the system frequency considering only the system inertial response. 2. Extreme frequency f m : Maximum or minimum frequency in the dynamic process. At the receiving-end, if the frequency nadir is less than the threshold of the under-frequency load shedding, regional load shedding is conducted, resulting in large-scale blackout.

Qualitative analysis of the frequency security indexes
Based on Equation (14), the approximate analytic formula for each frequency security index can be derived.
1. Initial rate of change of frequency, S 0 : As the initial rate of change of frequency, S 0 occurs within 1 s after the disturbance. Thus, substituting t as 0 and 1 in Equation (14) yields According to the validation of the simulation system and real system, the value of the coefficients 2 , , and T 2 is 0.13-0.15, 1.4-1.9, and more than 10, respectively; therefore, the values of sin 2 is approximately 0, and the values of cos and e −1∕T 2 are approximately 1. The analytic formula of S 0 can thus be expressed as: 2. Extreme frequency f m : We observe the time at which the maximum frequency deviation ∆f m appears and define this time as T m . Usually, T m pertains to the time 8.5-10 s after the disturbance. At T m , because T m /T 1 is large, the monotonic decay term is extremely small and can be ignored. Therefore, the frequency nadir mainly depends on the oscillatory decay term. The analytic formula for T m is: Furthermore, 3. Frequency recovery time T r : The total time from the start of the disturbance to the time at which the oscillatory decay term approaches the steady state is defined as the frequency recovery time T r . The steady state is the period in which the value of the oscillatory decay term fluctuates in a small range, and the boundary value f r is The value of the frequency recovery time T r depends on the selection of the boundary value f r , which is considered as ±0.01 Hz in this paper. The analytic formula of T r is 4. Steady-state frequency f ∞ : When the time approaches infinity, the system frequency reaches the steady state. The analytic formula for f ∞ is

Quantitative calculation of the frequency security indexes
Assuming that the steady-state system frequency before the disturbance is f 0 , and the observation duration and interval are T and t s , respectively, according to Equation (14), the curve of the system frequency f(k) can be defined as: 1. Initial rate of change of frequency S 0 : Based on the value of the observation interval t s , points k 1 and k 2 are selected to characterize the change of frequency in 1 s after the disturbance as k 1 = 0, k 2 = 1/t s . Therefore, S 0 can be calculated as 2. Extreme frequency f m : f m can be calculated as 3. Frequency recovery time T r : The selections of the boundary value f r and observation interval t s directly affect the calculation of the frequency recovery time, and for curves with dense points, T r can be calculated as follows.
4. Steady-state frequency f ∞ : By considering the average of the last n points, f ∞ can be calculated as Considering the approximate analytic formula of each frequency security index, the qualitative analysis results can be summarized as follows.
a. The steady-state frequency depends only on the constant term, which is proportional to the coefficient C 0 . b. The extreme frequency is independent of the monotonic decay term and thus independent of coefficients C 1 and T 1 . c. The initial rate of change of frequency is related to the decay term and mainly depends on the monotonic decay term, which is proportional to the coefficient C 1 and inversely proportional to the coefficient T 1 . d. The frequency recovery time depends only on the oscillatory decay term, which is proportional to the coefficients C 2 and T 2 .

Case 1: New England 10-unit 39-bus power system
The simulation is based on the New England 10-unit 39-bus power system, including hydro, thermal, and wind power units, as shown in Figure 3. The simulation is conducted using the software PSD-BPA, which is a widely used power system simulation software developed by the China Electric Power Research Institute. The model structures of the prime mover and its governor in the hydro and thermal units are described in Sections 5.1.2 and 5.1.3 of [24], respectively. The model structures of the controller for the DFIG-based wind power generators are described in Section 6.2 of [24]. The settings of the generators are listed in Table 1. We set a load increase of 5% as the power disturbance when the system frequency is 50.00 Hz. The actual data of the system frequency response are obtained based on the detailed model. In this case, the analytic formula of the frequency response Δ f (t ) can be uniquely determined as To validate the adaptability of the formula to the distribution of load increase, we set three cases with same load increase of 5%. The total initial load in the case study is 6192.8 MW, and load increase occurs in some randomly selected buses, and the detailed amount is listed in Table 2. The analytic formula of this paper is based on the data of Case I, while Case II and Case III have completely different load distributions. Figure 4 shows that the analytic formula can fit the other two cases well. The errors of frequency security indexes listed in Table 3 are within acceptable range.
To validate the adaptability of the formula to the power disturbance, we consider three total load increase values of +2.5%, +5.0% and +7.5% and obtain the curves of the frequency, as shown in Figure 5. The errors of the frequency security indexes under different load increase values against the actual data are listed in Table 4. Table 4 and Figure 5 indicate that the out- put of the analytic formula is close to the actual system frequency response; in particular, the errors of the frequency security indexes are small. Moreover, the analytic formula exhibits a high adaptability to the different degrees of disturbances and it can thus considerably reduce the computation in practice.

Qualitative analysis
The validation demonstrates that the analytic formula is sufficiently accurate to be applied for the analytic assessment of the frequency security. Next, in order to compare with the classical SFR model, the curves of the three terms of the analytic formula are considered, and the corresponding influence on the frequency dynamic response in different periods is analysed. Figure 6 shows the characteristics of the three terms of the analytic formula and concludes that each term has a unique effect on the system frequency. Moreover, following conclusions can be derived from Figure 6: a. The monotonic decay term is effective only in the initial period of the frequency dynamic response, in which it affects the initial rate of change of frequency and later promptly decays to almost zero. b. The oscillatory decay term is a key part of the frequency dynamic process, as it determines the overall trend and notably affects the frequency nadir. This term decays to almost zero after a relatively long time, indicating the end of the frequency dynamic process. c. The constant term determines the steady-state frequency.   In addition to the direct relationship with the power disturbance ∆P d , the relationship between the frequency security indexes and formula parameters is noteworthy. To highlight the possible linear relationship between the parameters and indexes, this paper considers the extreme frequency f m and steady-state frequency f ∞ in terms of the maximum frequency deviation ∆f m and steady-state frequency deviation ∆ f ∞, respectively. As shown in Figure 7, at the original operating point, different proportional changes are introduced in the formula parameter X, as follows. According to Figure 7, the relationship between the indexes and parameters is consistent with that noted in the qualitative analysis. Moreover, the extremely small error presented in Table 5 indicates that the proposed approximate formula can be used as the analytic expression for the frequency security index.

System parameter influence analysis
The qualitative analysis can help determine the formula parameters, validate the accuracy of the proposed approximate formula and clarify the relationship between the frequency security indexes and formula parameters. Because the system parameters with a clear physical meaning correspond to a higher controllability compared with the determined formula parameters, it is necessary to analyse the influence of the system parameters on the frequency security. According to the s-domain function of the frequency response (6), the formula parameters can be uniquely determined in terms of the coefficients of the fourthorder system transfer function, which are related to parameters in the generic SFR model. Among the parameters, three system parameters can be expressed as follows.
Therefore, a clear relationship exists between the formula parameters and system parameters. With reference to the qualitative analysis, Figure 8 shows the influence of the three system parameters on the frequency security indexes, which can be summarized as follows: a. The system inertial response is reflected in the monotonic decay term because the system inertia H considerably influences the coefficients C 1 and T 1 . As shown in Figure 7a, this parameter reduces the initial rate of change of frequency by 10%.  b. The system damping and system primary frequency modulation are reflected in the oscillatory decay term and constant term. The increase in the system damping K D adversely influences all the frequency security indexes by more than 10%.
In particular, the variation in the maximum frequency deviation and frequency recovery time can reach 40% owing to the notable influence of the parameter on the coefficients C 2 and T 2.
c. The primary frequency regulation capability of the unit K G considerably influences the steady-state frequency deviation by up to 10%, similar to the system damping K D . In addition, the influence of this parameter on the maximum frequency deviation is only 5%, mainly affected by the dynamic characteristics of the governor.

Case 2: East china power grid
This section describes the analysis of the practical application of the proposed approach in the East China Power Grid. The East China Power Grid is a regional power grid with the largest electric load in China, which provides electricity for Shanghai city, Jiangsu province, Zhejiang province, Fujian province, and Anhui province. We   The Prony analysis can be conducted to fit a reduced-order model to a high-order system. Prony algorithm is used to extract the amplitude, phase, frequency and damping factor of the signal [27,28]. A discrete time signal can be described as where A i is the amplitude, θ i is the phase, f i is the oscillation frequency, α i is the damping factor, N is the number of data, p is the order, and △t is the sample interval. Considering the measured data, the Prony analysis is adopted to validate the composition of the frequency response analytic formula. Table 6 lists the main components based on the energy of each mode in the descending order. Modes 1 and 2 occupy the major proportion, corresponding to the monotonic and oscillatory decay terms, respectively. Therefore, the analytic formula of the system frequency response contains three terms, the constant term and the monotonic and oscillatory decay terms. Compared with the classical SFR model, the analytic formula based on the generic SFR model can more accurately reflect the system frequency response.
Using the actual data measured, the analytic formula of the frequency response Δ f (t ) can be uniquely determined: As shown in Figure 9, the output of the analytic formula is close to the actual system frequency responses. Furthermore, the qualitative calculation of the frequency security indexes and errors can be realized. As shown in Table 7, the listed relative errors are small and acceptable.

CONCLUSION
A method to realize the analytic assessment of the power system frequency security is developed. An analytic formula of the system frequency response is obtained based on a generic SFR model, which consists of three components: the constant term and the monotonic and oscillatory decay terms. Compared to the existing approach, the proposed formula with the monotonic decay term can more accurately reflect the system frequency response, especially in the initial period. Moreover, considering the characteristics of the frequency dynamic process in different periods, four main indexes are identified to realize the assessment of the frequency security. Furthermore, the qualitative analysis of the frequency security index is conducted considering the formula parameters, system parameters and other influencing factors. In addition, considering the actual frequency data, the quantitative calculation method for the frequency security index is established. Case studies on an IEEE test system and a practical power system demonstrate that the proposed analytic assessment of the frequency security corresponds to a reduced computational time and higher accuracy.