An Analytical Method for Estimating the Maximum Penetration of DFIG Considering Frequency Stability

: With the increasing wind power in power systems and the wide application of frequency regulation technology, the accurate calculation of the limit wind power capacity in systems is critical to ensure the stability of the frequency and guide the planning of wind power sources. This paper proposes an analytical method for calculating the maximum wind generation penetration under the constraints of frequency regulation control and frequency stability taking doubly fed induction generator as an example. Firstly, the frequency-domain dynamic model of the doubly fed induction generator is established considering the supplementary frequency proportion-di ﬀ erentiation control under small disturbance. The equivalent inertia time constant of the doubly fed induction generator is calculated. On this basis, the frequency response model of the power system with the consideration of wind power integration in frequency regulation control is constructed. Then, the frequency-domain analytical solution of the system frequency is obtained. Finally, with the constraint by the steady-state deviation and dynamic change rate of the system frequency, the maximum wind generation penetration is analytically solved. The accuracy of the proposed analytical calculation method for the limit value of the percentage of wind power is veriﬁed by MATLAB / Simulink.


Background and Motivation
Wind power is generally regarded as the most promising renewable energy [1,2]. At present, the most widely used doubly fed induction generator (DFIG) has an excellent performance of capturing the maximum wind energy at real-time wind speed. However, its rotor speed is no longer coupled with the system frequency. The rotational kinetic energy contained in the rotor cannot be released, resulting in a drop in the inertia of the power system. The continuous increase of wind power capacity seriously threatens the stability of the frequency [3,4]. Moreover, frequency is one of the most critical factors restricting the acceptable capacity of wind power in the system [5,6]. For better guiding power planning and maintaining the safe and stable operation of the system, it has essential significance to estimate the maximum wind generation penetration, which is defined as the ratio of the maximum capacity of the acceptable wind power to the system capacity [7], based on the frequency constraints [8,9].

Literature Review
(1) Dynamic modeling of the power system In order to study the limit wind power capacity from the perspective of frequency stability, the prerequisite is to develop the dynamic model of the system. The typical models describing the frequency response of the power system include the full-time domain simulation model, average system frequency (ASF), and system frequency response (SFR) [10][11][12][13][14][15][16]. The full-time domain simulation model can simulate the dynamic equipment and network of the system in detail. However, it has a high dimensionality and a large amount of calculation [10,11]. Both ASF and SFR models assume that the system frequency remains the same and ignore the reactive power-voltage dynamic effects. They can intuitively reflect the relationship between active power balance and system frequency, and significantly reduce the amount of calculation [12,13]. However, the ignorance of the frequency control of wind power has made these models no longer suitable for the study of the dynamic frequency characteristics of the grid [14]. In this regard, Ghosh S et al. [15,16] considered the effects of wind power virtual inertial response and primary frequency regulation response. In [15], the inherent inertia time constant of wind power was added to the improved ASF model. However, the effect of the inertial response of DFIG is time-varying, which will bring a biased result if using an inherent constant to represent. Yan et al. [16] added the primary frequency-regulation droop control link of DFIG as a static model to the improved ASF model, which cannot fully reveal the dynamic response characteristics of the primary frequency regulation control of the DFIG.
(2) Assessment method For the measurement of the wind power capacity that can be connected to the power system, most of the current research has used simulation methods to solve the problem [17][18][19][20][21]. Zhang J. et al. [17] conducted many dynamic simulation tests based on the initial value of wind generation penetration they set at the beginning. According to various constraints, such as voltage and frequency, and the deviation of multiple indicators in each simulation, the limit value of the penetration power was continuously revised until the set constraint conditions were met. Although the penetration power can be obtained accurately, this method does not provide an analytical calculation method. The process of this simulation-based method is relatively complicated, and it must be re-operated when the system operation mode changes, which will bring a certain computational burden. In [18], the optimization methods were adopted to solve the maximum wind generation penetration with transient stability as a constraint and the maximum acceptable wind power capacity in the system as the objective. However, it ignored a large number of physical quantities and the frequency-regulation control of the DFIG, so the accuracy needs to be improved. Wang Man et al. [19] presented a novel three-stage optimization framework that targets the estimation of wind penetration limit based on the chance-constrained OPF incorporated with transient stability constraints of the power system. The analysis of dynamic model requires time series data, e.g., forecasting values of generations and demands [20]. With the popularity of wind power participating in primary frequency regulation, Shichun L. [21] established a primary frequency-regulation control response model of the DFIG and the frequency response model of the power system. On this basis, the researchers took the steady-state frequency deviation of the system and the frequency change rate at the initial moment of the disturbance as constraints and proposed a method for calculating the maximum wind generation penetration. However, the use of a higher-order wind power primary frequency-regulation model dramatically increases the computational complexity.
The above review of the literature indicates that the calculation of limit wind generation penetration is an issue that has not as yet been studied in sufficient depth. The limitations of the existing research are mainly manifested in the following aspects. Firstly, some of the previous research focused on the dynamic modeling of wind generators and primary frequency support. How to dynamically display the inertial response process and the primary frequency modulation process of the DFIG with Sustainability 2020, 12, 9850 3 of 21 additional frequency control in the system frequency dynamic modeling without greatly increasing the amount of calculations requires further research. Moreover, in the evaluation method of the maximum wind generation penetration, it is necessary to improve the accuracy and applicability of the calculation method and reduce the calculation complexity while considering the frequency regulation control of DFIG.

Proposed Method and Contributions
In this study, an analytical solution method for the maximum wind generation penetration that considers the frequency-regulation control of DFIGs and frequency stability constraints is proposed. In this regard, the main contributions of this study can be summarized as follows: (1) Contrarily to [15,16] which simply use static models to characterize the frequency regulation characteristics of DFIG, this paper calculates the equivalent inertial time constant in the form of the transfer function and adds a dynamic primary frequency-regulation link to characterize the inertial response and the primary frequency regulation process of the DFIG. This improvement can fully demonstrate the control function of DFIG in the system frequency dynamic response.
(2) By considering the participation of the DFIG in frequency regulation, the frequency response model of the power system with wind power in the time scales of the inertial response process and the primary modulation is established. The frequency-domain analytical solution of the system frequency is obtained.
(3) Based on the obtained expression of system frequency, the maximum wind generation penetration can be analytically calculated under the dynamic and static frequency stability constraints. Compared with solving through repeated simulation operations [17], this method is suitable for various situations and can fully demonstrate the control process of the DFIG while reducing the computational complexity. Compared with the optimization algorithm that ignores a large number of dynamic processes [18], this calculation can completely mathematically model the frequency dynamic response characteristics of the system, to better understand the whole process of the system frequency dynamic response.

Organization of the Paper
This paper is organized as follows: the dynamic modeling of the inertial response process of DFIG is analyzed in Section 2, and the equivalent inertia time constant of the DFIG is calculated in Section 3. In Section 4, the SFR model formulations of the power system with wind power are established. The maximum wind generation penetration considering frequency stability constraint is analytically solved in Section 5. The time-domain and assessment simulation results are provided in Section 6. Finally, the paper is concluded in Section 7.

Dynamic Modeling of Inertial Response Process of DFIG
In order to solve the problem in a targeted manner, this section establishes a simplified DFIG model suitable for studying frequency response. Considering that the inertial response has the dynamic response characteristics of the electromechanical time scale, the following assumptions were made. Firstly, it is thought that wind speed and mechanical torque remain unchanged during the inertial response dynamic process. Secondly, when the wind speed is greater than the rated wind speed for a long time, the wind turbine does not use inertial control but changes the pitch angle to realize frequency regulation due to the constraints of mechanical structure and electrical load bearing capacity. This process involves a more complex modeling process, which needs to be further discussed in future research. Therefore, only a case where the wind speed is less than the rated wind speed is discussed here. Thirdly, the secondary frequency regulation of the DFIG was disabled due to the limitation of the structure of wind turbines. In the modeling of the DFIG, the secondary frequency regulation is not considered. According to the above assumptions, the models closely related to inertial Sustainability 2020, 12, 9850 4 of 21 response include the generator electromechanical transient model, frequency regulation controller model, speed controller model, and power tracking control model.

Generator Electromechanical Transient Model
The transient electromechanical characteristics of DFIG can be represented by an equivalent mathematical model in a two-phase d-q rotating coordinate system. The flux equation of the generator can be expressed as where ψ ds , ψ qs are the d-axis and q-axis component of the stator flux linkage. ψ dr , ψ qr is the d-axis and q-axis component of the rotor flux linkage. i ds , i qs are the d-axis and q-axis component of the stator current. i dr , i qr are the d-axis and q-axis component of the rotor current. L m is the mutual inductance between the stator and rotor windings. L s , L r is the self-inductance of the stator winding and the rotor winding. ω r is the rotor speed.
Ignoring the influence of the shafting characteristics on the power grid characteristics, use [ψ ds , ψ qs , i dr , i qr , ω r ] as the state variables to establish a fifth-order transient model to show the electromagnetic torque T em directly related to the virtual inertial response: where p is the number of generator pole pairs. In a synchronous rotating coordinate system based on stator voltage vector control, ψ qs = 0. Therefore, the expression of electromagnetic torque can be simplified to where k 1 is the constant coefficient of T em . At the same time, the rotor motion equation of the generator is where H DFIG is the inherent inertia time constant of the DFIG. T mach is the mechanical torque of the generator.

Speed Controller Model
In steady state, DFIG realizes real-time speed tracking through the speed controller. When the speed of the DFIG changes, the inertial response process will be dynamically adjusted with the speed controller. The mathematical model describing the dynamic response process of the speed controller is where ∆T * vc is the output reference torque increment of the speed controller. K iv is the integral coefficient. K pv is the proportional coefficient. ∆ω r is the rate of change of the rotor speed of the DFIG.

Power Tracking Controller Model
The energy conversion efficiency of a variable speed wind turbine is determined by the wind energy utilization factor. When the DFIG is running in the maximum power tracking state in the initial state, the rotational speed usually tracks the wind speed to maintain the best tip speed ratio λ opt . At this time, the DFIG is running at the highest efficiency point, and the wind energy utilization coefficient takes the maximum value. In order to participate in power system frequency regulation, DFIGs need to reserve spare capacity in advance. For this purpose, overspeed control is usually selected to obtain spare capacity. At this time, the DFIG gives up capturing the maximum power and operates with deloaded power curves, which are pre-defined according to the requested active power output imposed by wind park supervisory control. The output power torque T* del under this condition is [22] T where k del is the derating power tracking coefficient of DFIG. When the wind speed is constant, and the system has a small disturbance, the electromagnetic torque increment ∆T * del provided by the power tracking controller is where ω r0 is the initial value of the rotor angular velocity of DFIG.

SFPD Controller Model
The existing DFIG mostly adopts the supplementary frequency proportion-differentiation (SFPD) frequency regulation control that combines droop control and virtual inertia control. The virtual inertia control adds an additional electromagnetic torque proportional to the change rate of the system frequency, which enables the DFIG to spontaneously realize the negative feedback effect of frequency change and release the kinetic energy stored in the generator rotor during the dynamic process. Therefore, the DFIG can quickly realize the virtual inertial response within a short period of time when the disturbance occurs. The droop control is the mainstream method for primary regulation, which takes the real-time frequency deviation ∆f of the system as the feedback quantity and adds an active power control signal through the proportional gain link to adjust the active output. Compared with virtual inertial control, this method mainly works on a longer time scale after the disturbance. It effectively increases the output active power of the wind turbine and restores the system frequency.
The reference value of additional electromagnetic torque ∆T* emu generated by frequency regulation control is where ∆T * emu is the output reference torque increment of the SFPD frequency regulation controller. K if is the proportional coefficient of the virtual inertia control link. K df is the proportional coefficient of the droop control link. T is the time constant of the low-pass filter. ∆ω s is the increment of the synchronous angular velocity of the system. As shown in Figure 1, by establishing the control model of each module, a simplified control model of the DFIG can be obtained. Moreover, T max and T min are the maximum and minimum reference torques that the speed controller can output. f s is the frequency of the power system [23].
Kif is the proportional coefficient of the virtual inertia control link. Kdf is the proportional coefficient of the droop control link. T is the time constant of the low-pass filter. Δωs is the increment of the synchronous angular velocity of the system. As shown in Figure 1, by establishing the control model of each module, a simplified control model of the DFIG can be obtained. Moreover, Tmax and Tmin are the maximum and minimum reference torques that the speed controller can output. fs is the frequency of the power system [23]. iv pv

Calculation of Equivalent Inertia Time Constant of DFIG
When DFIG adds SFPD frequency regulation control so that it can provide certain inertial support, the rotor speed of the DFIG is no longer decoupled from the angular velocity of the system. At this time, DFIG can support the system frequency by increasing active output according to frequency change when the system is disturbed. Furthermore, the supporting capacity can be

Calculation of Equivalent Inertia Time Constant of DFIG
When DFIG adds SFPD frequency regulation control so that it can provide certain inertial support, the rotor speed of the DFIG is no longer decoupled from the angular velocity of the system. At this time, DFIG can support the system frequency by increasing active output according to frequency change when the system is disturbed. Furthermore, the supporting capacity can be characterized by the defined equivalent inertia time constant. Its value reflects the virtual inertia response capability of the wind turbine relative to the grid frequency change.
According to [23], the variation of kinetic energy during the inertial response of DFIG is expressed by the deviation of rotor speed and inherent inertial torque and the variation of system synchronous angular velocity and the virtual moment of inertia of DFIG. Thus, an equation relationship between them can be established. The virtual moment of inertia of the DFIG can be obtained as where J equ is the virtual moment of inertia of the DFIG. ω s0 is the initial synchronous angular velocity of the system. J DFIG is the inherent moment of inertia of the DFIG. According to Equation (9) and the definition of inertia time constant, the virtual inertia time constant of the DFIG in the inertial response process can be defined as [23] where ω nom is the rated angular velocity of the DFIG. S N is the rated capacity of the DFIG. H eqDFIG quantitatively characterizes the inertial response capability of the DFIG when the system frequency fluctuates, and its size is determined by ∆ω r /∆ω s . At the same time, the DFIG mainly realizes the inertial response process by changing the electromagnetic torque, so H eqDFIG is mainly related to the electromagnetic torque increase. Because the SFPD-integrated frequency regulation controller, power tracking controller, and speed controller jointly affect the inertial response process of DFIG, the total electromagnetic torque ∆T em can be expressed as the sum of their output torque increments. By combining Equations (5)-(8), we obtain At the same time, by considering the inertial response process of a smaller time scale and ignoring the change in the rotor speed of the DFIG, the frequency-domain expression obtained from the rotor motion equation can be obtained from (4) 2H DFIG · ∆ω r s = −∆T em (12) Sustainability 2020, 12, 9850 Bringing in Equation (19), we obtain Combined with Equation (10), the frequency-domain expression of H eqDFIG can be obtained as

Frequency Response Model of Power System without Wind Power
When a transient disturbance occurs in the system, the frequency is dynamically adjusted and restored to a steady-state value under the regulation of synchronous unit and load. The system frequency response process has a relatively long time scale. In order to facilitate calculation, the links with a smaller time scale can be ignored, and the system model can be appropriately simplified. In [24][25][26], a simplified second-order frequency response model is proposed. The system frequency response process only considers the generator rotor and steam turbine reheater with a large inertia time constant [24], as shown in Figure 2. The model assumes that the synchronous generators are all driven by reheat turbines, and all conventional generators in the system are equivalent to one turbine generator. H is the equivalent inertia time constant of synchronous generator set after equivalent calculation. R is the droop coefficient of the governor. TR is the time constant of the turbine reheater. FH is the proportion of the output power of the high-pressure boiler. K m is the mechanical power gain factor. ∆P G is the change in the mechanical power of the generator. ∆Pa is the acceleration power assumed by the rotor. ∆P d is the disturbance power. D is the load-damping coefficient.
Sustainability 2020, 01, x FOR PEER REVIEW 7 of 21 response process only considers the generator rotor and steam turbine reheater with a large inertia time constant [24], as shown in Figure 2. The model assumes that the synchronous generators are all driven by reheat turbines, and all conventional generators in the system are equivalent to one turbine generator. H is the equivalent inertia time constant of synchronous generator set after equivalent calculation. R is the droop coefficient of the governor. TR is the time constant of the turbine reheater.
FH is the proportion of the output power of the high-pressure boiler. Km is the mechanical power gain factor. ΔPG is the change in the mechanical power of the generator. ΔPa is the acceleration power assumed by the rotor. ΔPd is the disturbance power. D is the load-damping coefficient. According to the model in Figure 2, the equation of the system frequency response can be listed as where Gs(s) is the simplified primary frequency regulation transfer function model of conventional thermal power units. Δω is the system frequency deviation per unit value. Its meaning is consistent with the above Δωs.
Usually, a step function is used to describe the instantaneous power disturbance of the system as such: where ΔPstep is the disturbance magnitude in per unit. We assume that all equations are on a According to the model in Figure 2, the equation of the system frequency response can be listed as where G s (s) is the simplified primary frequency regulation transfer function model of conventional thermal power units. ∆ω is the system frequency deviation per unit value. Its meaning is consistent with the above ∆ω s . Usually, a step function is used to describe the instantaneous power disturbance of the system as such: where ∆P step is the disturbance magnitude in per unit. We assume that all equations are on a common system base S B , which is equal to the sum of the ratings of all generating units in the system. ∆P L is the actual size of the disturbance power. By applying the model shown in Figure 2 to a power system with n synchronous generators, the inertia time constant of the equivalent machine and the droop coefficient of the governor in the model can be expressed as where S i , H i and R i are, respectively, the rated capacity, inertia time constant and droop coefficient of the i-th synchronous machine.

Frequency Response Model of Power System with Wind Power
Considering that wind power is connected to the power system, wind power and synchronous generator sets balance the load power together. The grid conditions are basically the same. The system can be set to contain n synchronous generators and m wind turbines. Define the wind generation penetration in the power system α as the ratio of grid-connected wind power capacity to the total installed capacity of the system. That is, where α is the wind generation penetration. S j is the rated capacity of the j-th wind turbine. n and m are the number of synchronous machines and wind turbines, respectively.
At this time, the inertia time constant of the equivalence machine in the model and the droop coefficient of governor in Figure 2 can be expressed as where H and R are the inertia time constant of the equivalence machine, and the droop coefficient of the governor after the wind turbines participating in the frequency control are connected to the power system. It can be seen that the wind turbines participating in frequency control are connected to the system, which will change the equivalent inertia time constant of the system and the droop coefficient of the governor.
When DFIG adopts the SFPD frequency regulation control method, the droop control link will affect the primary frequency regulation process of the system on a longer time scale. At this time, the change of speed cannot be ignored, and it is necessary to consider the dynamic process to be larger than part of the inertial response time scale. Therefore, when considering the dynamic modeling of the power system, it is necessary to establish a dynamic model of the primary frequency regulation of the DFIG. This model describes the dynamic process in which the unit adjusts the output mechanical power to suppress the frequency disturbance when the power system frequency is disturbed. It can be obtained by solving the transfer function of unit output power increment and system frequency deviation.
When the DFIG adopts the integrated controller to participate in the system frequency adjustment, a certain reserve capacity needs to be reserved. When the DFIG participates in frequency regulation by obtaining reserve capacity through overspeed control, the decrease in rotor speed during frequency regulation causes the operating point of the deloaded power curve to move downward. Part of the active power increment is used to provide steady-state power support for the system. The other part needs to be used to compensate for the power loss caused by the downward movement of the deloaded operating point in the load-shedding power tracking curve. Therefore, during the frequency disturbance accident, the electromagnetic power change of the DFIG ∆P e can be expressed as where 3k del ω r0 2 ∆ω r 2 is equivalent to the mechanical power loss caused by the lowering of the deloaded operating point during the frequency adjustment of the DFIG. At the same time, the mechanical power P m captured under deloaded operation of the DFIG can be expressed as where ρ is the air density. S is the area swept by the blade of the DFIG. v is the wind speed. C p is the wind energy utilization coefficient under deloaded operation of the DFIG.
Assuming that the wind speed remains constant during the primary frequency regulation response, the wind speed can be set as the initial steady-state wind speed v 0 . Then, P m is uniquely determined by C p . There is a strong nonlinear implicit function relationship between C p and blade tip speed ratio λ and pitch angle β. In order to simplify the order in the actual modeling calculation, the wind energy utilization coefficient C p is simplified to an explicit functional relationship between λ and β [27]. It is assumed that the blade response is not decoupled from any significant participation in the torsional dynamics. Without individual representation of the two blades, there is no need to include wind shear and tower shadow effects at a detailed level. The wind models which are available include an additive term representing the net effect of rotational non-linearities in the driving torque: where the tip speed ratio λ ωR/v. ω is the rotation speed of the wind wheel. R is the radius of the wind wheel. d% is the load reduction ratio. Assume that the wind speed v remains constant during a frequency modulation response. The value of the wind energy utilization coefficient C p is only determined by ω and β in a short time. Then, we can obtain the small-signal incremental of C p (λ, β) that can be expressed in ω r and β: ∆C p (λ, β) = ∆C p (ω r , β) = ∂C P ∂ω r ∆ω r + ∂C P ∂ω r β = (y 0 +y 1 )∆ω r + (y 2 +y 3 + y 4 )β where Sustainability 2020, 12, 9850 10 of 21 When using the speed control method, if the wind speed is less than the rated wind speed, the pitch angle will not act. Considering β, we obtain Therefore, the wind energy utilization coefficient function ∆C p can be expressed as a linear function of ∆ω r under deloaded operation of the DFIG in the simplified calculation: where ∂C p /∂ω r ∆ω r is the partial lead of the wind energy utilization coefficient function to the speed under deloaded operation of the DFIG. ∂C p /∂ω r is the coefficient of the partial derivative. When the pitch angle is 0, its value is 0.29π(1 − d%)cos(πRω r0 /15v w0 − 0.2π).
When the wind speed remains constant and the pitch angle is fixed, the rotational speed affects the wind energy utilization coefficient function, which in turn affects the mechanical power change ∆P m captured by the DFIG.
In the primary frequency regulation stage, the change in mechanical torque cannot be ignored due to the large change in speed. Similarly, combining the small-signal analysis method and the rotor motion equation of the DFIG, we obtain By incorporating Equations (18) and (24) into Equation (25) and eliminating the intermediate variable ∆ω r , the frequency regulation response model of the DFIG can be obtained as follows: Substituting the parameters of Equation (30) into Figure 2, the modified power system frequency response model, after the wind turbines participating in the frequency regulation control are connected, can be obtained, as shown in Figure 3.
By incorporating Equations (18) and (24) into Equation (25) and eliminating the intermediate variable Δωr, the frequency regulation response model of the DFIG can be obtained as follows: Substituting the parameters of Equation (30) into Figure 2, the modified power system frequency response model, after the wind turbines participating in the frequency regulation control are connected, can be obtained, as shown in Figure 3. According to Figure 3, it can be given that According to Figure 3, it can be given that Therefore, the frequency-domain expression of the deviation per unit value of frequency is

Maximum Penetration of DFIG Constrained by Steady-State Frequency Deviation
The steady-state frequency deviation (SFD) can reflect the ability of the system to resist active disturbances through frequency regulation. If the SFD is large, it will not only affect the power quality but also increase the risk of system frequency instability when faced with large active power disturbances. Therefore, the power system usually stipulates that the SFD should be within ±0.2 Hz. According to the final value theorem of Laplace transform, the SFD of (28) can be obtained as where f N is the rated frequency of the system. Assuming that the boundary of the steady-state frequency deviation of the system is ξ, the SFD satisfies SFD ≤ |ξ| (32) By combining (29) and (30), the maximum penetration of DFIG under the frequency stability constraint can be obtained: where α max1 is the maximum wind generation penetration of the wind power under the constraint of steady-state frequency deviation. It can be seen that the maximum wind generation penetration is mainly affected by factors such as the droop control coefficient K df of the DFIG, the load change level ∆P L /S B , and the load reduction ratio of DFIG when constrained by the steady-state frequency deviation. The load reduction ratio of DFIG mainly affects the partial conduction of mechanical power to speed in the load reduction condition. The droop control coefficient K df has a more significant influence on the steady-state frequency deviation, and the virtual inertia control coefficient has a small impact on the steady-state frequency deviation.

Maximum Penetration of DFIG Constrained by Dynamic Rate of Change of Frequency
The rate of change of frequency (ROCOF) reflects the speed of frequency change of the power system when the system is disturbed. The frequency change rate is an important indicator of the dynamic frequency response of the system, which is directly related to whether the system frequency reaches the dynamic safety limit. The starting fixed value of the load shedding action of the existing frequency change rate relay is usually set at ±0.05 Hz/s~0.1 Hz/s. Therefore, the boundary value of the frequency change rate is η = ±0.05 Hz/s in the calculation of this article. Usually, the frequency change rate reaches its maximum value at the beginning of the frequency disturbance caused by power change. According to the Laplace initial value theorem, the maximum value of the frequency change rate can be obtained as follows: Assuming that the boundary of the system frequency change rate is η, ROCOF satisfies ROCOF ≤ η By combining (32) and (33), it can be obtained that the maximum wind generation penetration under the constraint of dynamic frequency stability satisfies It can be seen that the maximum wind generation penetration has more influence factors when the dynamic frequency change rate is constrained. It mainly includes the droop control coefficient K df , the proportional coefficient K if of the virtual inertia control link, the load change level ∆P L /S B , and the load reduction ratio of the DFIG. The proportional coefficient K if of the virtual inertia control link has a significant influence on the result.
Considering the above two situations and taking the smaller value, the maximum wind generation penetration can be obtained as follows:

Case Study
This paper builds the model shown in Figure 4 on the MATLAB/Simulink simulation platform. This model is based on the detailed nonlinear simulation model of the Simulink module and is used to verify the accuracy of Equation (31). The model is based on the IEEE 9-node model. The three synchronous generators are all thermal power units that are equivalent to several units and are merged into the wind farm on bus 9. Its initial base capacity is 100 MVA, and then it changes with simulation conditions. The wind farm is composed of several DFIGs. The power system frequency response model parameters are shown in Table 1, and the parameters of DFIG are shown in Table 2. The number of DFIG in the wind farm is determined according to the simulation scenario. Assuming that the wind speed is constant at 9m/s within the time scale of the first frequency regulation response after the disturbance, the active power balance is maintained. The remaining parameters of the system are shown in [28].

Verification of the Calculation Method for Equivalent Inertia Time Constant of DFIG
In the above IEEE9-node simulation system, a sudden load increase of 20 MW was set at the position shown in Figure 4 at t = 20 s to simulate the system frequency transient when the active power was short. Through simulation, the ratio of Δωr to Δωs (Δωr/Δωs) during the inertial response of the DFIG is obtained, and then the simulation value of HeqDFIG is obtained.
In order to verify the correctness of the calculation results, the simulation results of HeqDFIG(t) in

Verification of the Calculation Method for Equivalent Inertia Time Constant of DFIG
In the above IEEE9-node simulation system, a sudden load increase of 20 MW was set at the position shown in Figure 4 at t = 20 s to simulate the system frequency transient when the active power was short. Through simulation, the ratio of ∆ω r to ∆ω s (∆ω r /∆ω s ) during the inertial response of the DFIG is obtained, and then the simulation value of H eqDFIG is obtained.
In order to verify the correctness of the calculation results, the simulation results of H eqDFIG (t) in the process of system frequency transients caused by active power shortage are compared with the calculation results. The simulation conditions are K df = 6, K if = 6, T = 0.1 s. Figure 5 is the curve of the inertial response H eqDFIG produced by the DFIG when the power is insufficient when the wind speed is 8 m/s and the time of failure t = 20 s. It can be seen that the calculated value is consistent with the simulated value, which proves the correctness of the calculated result. Unlike ordinary synchronous generators, which have a constant inertia time constant, the equivalent inertia time constant of the DFIG has time-varying characteristics.  Because the change rate of the system frequency is tremendous when the system disturbance occurs, the HeqDFIG at the initial moment of the inertial response quickly rises to a larger initial value. At that time, the DFIG releases the rotor kinetic energy, and the rotor speed drops. Moreover, when the change speed of the system frequency is higher than the change speed of the rotor speed of the DFIG, Δωr/Δωs decreases. Therefore, HeqDFIG gradually decreases from a larger initial value. Then, the inertial response of the SFPD controller slows down the system frequency change rate, and DFIG is still in the rotor kinetic energy release stage, so Δωr continues to increase. Then, Δωr/Δωs increases instead, so HeqDFIG has a rising process. Then, the control effect of the speed controller becomes more evident with the increase of Δωr, and the DFIG starts to restore the speed. As the system frequency gradually recovers, the inertial response process also ends. HeqDFIG transitions to a minimal steadystate value after small fluctuations.

Verification of Power System Frequency Response Model
This section uses the detailed nonlinear simulation model based on the Simulink module as shown in Figure 4 to verify the accuracy of the established reduced-order power system frequency response model. The reduced-order power system frequency response model parameters of Figure 4 are calculated as shown in Table 3. The wind farm contains 90 DFIGs, which is equivalent to a DFIG with a rated capacity of 150 MW. The wind generation penetration in the system is 33.33%. The loadshedding ratio of DFIG is 10%. The wind speed is 9m/s. Table 3. Power system frequency response model parameters.

Classification
Parameter Value

Synchronous generator
Droop coefficient of steam turbine R 0.05 The time constant of steam turbine reheater TR 10 s The proportion of output power of high-pressure boiler FH 0.28 The mechanical power gain factor Km 1
The inherent inertia time constant of the rotor H 5.04s The partial conductance of electromagnetic power to speed3kdelω 2 ro_del 0.78 The partial conduction of mechanical power to speed −0.3572 We supposed that the system load Load1 power suddenly increases at t = 50 s. The parameter Because the change rate of the system frequency is tremendous when the system disturbance occurs, the H eqDFIG at the initial moment of the inertial response quickly rises to a larger initial value. At that time, the DFIG releases the rotor kinetic energy, and the rotor speed drops. Moreover, when the change speed of the system frequency is higher than the change speed of the rotor speed of the DFIG, ∆ω r /∆ω s decreases. Therefore, H eqDFIG gradually decreases from a larger initial value. Then, the inertial response of the SFPD controller slows down the system frequency change rate, and DFIG is still in the rotor kinetic energy release stage, so ∆ω r continues to increase. Then, ∆ω r /∆ω s increases instead, so H eqDFIG has a rising process. Then, the control effect of the speed controller becomes more evident with the increase of ∆ω r , and the DFIG starts to restore the speed. As the system frequency gradually recovers, the inertial response process also ends. H eqDFIG transitions to a minimal steady-state value after small fluctuations.

Verification of Power System Frequency Response Model
This section uses the detailed nonlinear simulation model based on the Simulink module as shown in Figure 4 to verify the accuracy of the established reduced-order power system frequency response model. The reduced-order power system frequency response model parameters of Figure 4 are calculated as shown in Table 3. The wind farm contains 90 DFIGs, which is equivalent to a DFIG with a rated capacity of 150 MW. The wind generation penetration in the system is 33.33%. The load-shedding ratio of DFIG is 10%. The wind speed is 9 m/s. We supposed that the system load Load1 power suddenly increases at t = 50 s. The parameter values of the load sudden-increased power and SFPD frequency regulation control parameters were respectively changed for comparison and simulation. The simulation case includes three scenarios shown in Table 4. Scenario 1 and scenario 2 maintain the SFPD frequency regulation control parameters K if and K df unchanged and change the load sudden-increased power ∆P L . Scenario 3 and scenario 2 maintain the load burst power ∆P L unchanged and change the SFPD frequency regulation control parameters K if and K df unchanged. Table 3. Power system frequency response model parameters.

Classification
Parameter Value

Synchronous generator
Droop coefficient of steam turbine R 0.05 The time constant of steam turbine reheater T R 10 s The proportion of output power of high-pressure boiler F H 0.28 The mechanical power gain factor K m 1

DFIG
The initial speed at 10% load reduction ω r0_del 1.0656 p.u. The inherent inertia time constant of the rotor H 5.04 s The partial conductance of electromagnetic power to speed3k del ω 2 ro_del 0.78 The partial conduction of mechanical power to speed −0.3572 In scenario 1, the time, frequency and steady-state frequency deviation Based on the above analysis, it can be seen that the nonlinear simulation model based on Simulink and the frequency response model based on Figure 3 have a high degree of agreement in different scenarios. By comparing Scenario 2 and Scenario 3, it can be seen that when the frequency regulation control parameter is fixed, the system frequency deviation increases when the load sudden-increased power increases. By comparing Scenario 1 and Scenario 2, it can be seen that when the load sudden-increased power remains unchanged, the steady-state deviation of the system frequency is reduced after DFIG participates in the frequency regulation compared to when DFIG does not add SFPD frequency regulation control. Moreover, the frequency drop speed is reduced.
In scenario 1, the time, frequency and steady-state frequency deviation of the system frequency reaching the lowest point in the Simulink nonlinear simulation model are 54.23 s, 49.657 Hz and −0.178 Hz, respectively. The active power increment of DFIG is 0.0054 p.u. Based on Equation (37), the analytical calculation model system frequency reaches the lowest point time, frequency and steady-state frequency deviation are 54 s, 49.661 Hz and −0.18 Hz, respectively, and the wind turbine active power increment is 0.0057 p.u. Combined with Figure 6a,b, it can be seen that the frequency response characteristics obtained from the Simulink nonlinear simulation model and the analytical calculation model based on (37)

Verification of the Maximum Wind Generation Penetration under Different K df
We supposed that the system load Load1 power suddenly increases at t = 50 s. The parameter values of the load sudden-increased power and SFPD frequency regulation control parameters were respectively changed for comparison and simulation. The simulation case includes three scenarios shown in Table 5. The maximum wind generation penetration in the system can be calculated by Equation (31) are 16.67%, 18.09%, and 19.77%, respectively. The detailed nonlinear simulation model based on Simulink verifies the accuracy of the calculated maximum wind generation penetration through time-domain simulation. When the maximum wind generation penetration is 16.67%, 18.09%, and 19.77%, the number of wind farms in the time domain simulation model shown in Figure 4 contains 38, 41, and 44 DFIGs, respectively. The power of Load4 is 20.5 MW, 21.9 MW, and 23.25 MW, respectively. When t = 50 s, the load sudden change power is 24.255 MW, 24.545 MW, and 24.878 MW, respectively. Compare the time-domain simulation of three extreme scenarios of the wind power ratio. Figure 7 shows the system frequency response curve. Under different primary frequency regulation control parameters, the system frequency stabilizes back to around 49.8 Hz after transient disturbance. The steady-state frequency deviation is restricted to the maximum allowable deviation value of 0.2 Hz, and the frequency change rate is restricted to 0.05 Hz/s. The simulation results verify the accuracy of the calculated maximum wind generation penetration. By comparing the simulation results under different frequency-regulation control coefficient scenarios K df , it can be seen that the maximum wind generation penetration is smaller without droop frequency-regulation control of DFIG(K df = 0). With the addition of droop control to wind power, the maximum wind generation penetration is increased. Moreover, as the droop control coefficient K df of DFIG increases, the maximum wind generation penetration grows. This is due to the increase in the active power contribution of DFIGs with droop control during the primary frequency regulation stage during the frequency disturbance of the system. When K df increases from 0 to 10, the maximum wind generation penetration increases by about 3%.

Verification of the Maximum Wind Generation Penetration under Different Kif
We supposed that the system load Load1 power suddenly increases at t = 50s. The parameter values of the load sudden-increased power and SFPD frequency regulation control parameters were respectively changed for comparison and simulation. The simulation case includes three scenarios shown in Table 6. The maximum wind generation penetration that can be calculated by Equation (31) are 17.21%, 18.01%, and 18.99%, respectively. The detailed nonlinear simulation model based on the Simulink module verifies the accuracy of the calculated maximum wind generation penetration through time-domain simulation. When the maximum wind generation penetration is 17.21%, 18.01%, and 18.99%, the numbers of DFIGs in the time domain simulation model, shown in Figure 5, is 39, 41, and 42, respectively. The power of Load4 is 20.8 MW, 21.05 MW, and 22.95 MW. When t = 50 s, the sudden load power is 21.88 MW, 24.5 MW, and 24.52 MW. Time-domain simulation is performed on several extreme scenarios of the wind power ratio. Figure 8 shows the system frequency response curve. In different Kif scenarios, the system frequency stabilized back to around 49.8 Hz after transient disturbance. The steady-state frequency deviation is restricted to the maximum allowable deviation value of 0.2 Hz, and the

Verification of the Maximum Wind Generation Penetration under Different K if
We supposed that the system load Load1 power suddenly increases at t = 50s. The parameter values of the load sudden-increased power and SFPD frequency regulation control parameters were respectively changed for comparison and simulation. The simulation case includes three scenarios shown in Table 6. The maximum wind generation penetration that can be calculated by Equation (31) are 17.21%, 18.01%, and 18.99%, respectively. The detailed nonlinear simulation model based on the Simulink module verifies the accuracy of the calculated maximum wind generation penetration through time-domain simulation. When the maximum wind generation penetration is 17.21%, 18.01%, and 18.99%, the numbers of DFIGs in the time domain simulation model, shown in Figure 5, is 39, 41, and 42, respectively. The power of Load4 is 20.8 MW, 21.05 MW, and 22.95 MW. When t = 50 s, the sudden load power is 21.88 MW, 24.5 MW, and 24.52 MW. Time-domain simulation is performed on several extreme scenarios of the wind power ratio. Figure 8 shows the system frequency response curve. In different K if scenarios, the system frequency stabilized back to around 49.8 Hz after transient disturbance. The steady-state frequency deviation is restricted to the maximum allowable deviation value of 0.2 Hz, and the frequency change rate is restricted to 0.05 Hz/s. The above simulation results verify the accuracy of the calculated maximum wind generation penetration of wind power. By comparing the simulation results under different virtual inertia control coefficient scenarios K df , it can be seen that the maximum wind generation penetration is smaller without virtual inertia control of DFIG (K df = 0). With the addition of virtual inertia control to wind power, the maximum wind generation penetration has been increased, and as the wind power virtual inertia control coefficient K if increases, the maximum wind generation penetration grows. This is because the DFIG can release more rotor kinetic energy during the frequency disturbance of the system, which increases the contribution to the inertial response process of the system. When K if increases from 0 to 10, the maximum wind generation penetration increases by about 1.9%. Based on this, it can be concluded that the additional virtual inertia control of DFIG will promote wind power consumption to a certain extent. during the frequency disturbance of the system, which increases the contribution to the inertial response process of the system. When Kif increases from 0 to 10, the maximum wind generation penetration increases by about 1.9%. Based on this, it can be concluded that the additional virtual inertia control of DFIG will promote wind power consumption to a certain extent.

Conclusions
This paper studies the maximum wind generation penetration under wind power frequencyregulation control and system frequency constraints. Firstly, the calculation method of the equivalent inertia time constant of the DFIG with frequency regulation control is given. Unlike the invariable inertia time constant that synchronous generators have, the equivalent inertia time constant of the DFIG has a multistage time-varying characteristic expressed as a transfer function. Then, a frequency response model of the power system with SFPD control of DFIG is established. The influence of wind power connected to the system on frequency stability is mainly due to the addition of a primary frequency regulation control loop while affecting the equivalent inertia time constant of the system. Next, the analytical method is used to solve the maximum wind generation penetration under the frequency constraint. The main factors affecting the maximum wind generation penetration are the frequency regulation control coefficient of the DFIG, the level of load change, and the proportion of deloaded. Finally, the simulation analysis verifies the correctness of the results and proves that the participation of wind power in frequency regulation control will effectively promote wind power consumption. Compared with wind power that does not participate in frequency regulation, the maximum wind generation penetration when wind power is controlled by additional frequency regulation is effectively increased. Moreover, when the frequency regulation control parameter droop control gain and virtual inertia control gain of the DFIG are larger, the maximum wind generation penetration is higher. How to further apply the method to more different types of renewable energy such as photovoltaics and other types of wind turbines will be the focus of future work.

Conclusions
This paper studies the maximum wind generation penetration under wind power frequency-regulation control and system frequency constraints. Firstly, the calculation method of the equivalent inertia time constant of the DFIG with frequency regulation control is given. Unlike the invariable inertia time constant that synchronous generators have, the equivalent inertia time constant of the DFIG has a multistage time-varying characteristic expressed as a transfer function. Then, a frequency response model of the power system with SFPD control of DFIG is established. The influence of wind power connected to the system on frequency stability is mainly due to the addition of a primary frequency regulation control loop while affecting the equivalent inertia time constant of the system. Next, the analytical method is used to solve the maximum wind generation penetration under the frequency constraint. The main factors affecting the maximum wind generation penetration are the frequency regulation control coefficient of the DFIG, the level of load change, and the proportion of deloaded. Finally, the simulation analysis verifies the correctness of the results and proves that the participation of wind power in frequency regulation control will effectively promote wind power consumption. Compared with wind power that does not participate in frequency regulation, the maximum wind generation penetration when wind power is controlled by additional frequency regulation is effectively increased. Moreover, when the frequency regulation control parameter droop control gain and virtual inertia control gain of the DFIG are larger, the maximum wind generation penetration is higher. How to further apply the method to more different types of renewable energy such as photovoltaics and other types of wind turbines will be the focus of future work.
Author Contributions: The authors confirm their contributions to the paper as follows: M.Q. and F.T. proposed the idea and wrote the paper; F.L. and D.L. revised the manuscript; N.D. and B.H. reviewed the results and approved the final version of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This paper study was funded by the National Natural Science Fund Program of China (grant number 51977157).