Analytical Approach to Understanding the Effects of Implementing Fast-Frequency Response by Wind Turbines on the Short-Term Operation of Power Systems

: The signiﬁcant presence of variable-speed wind turbines in worldwide power systems has led to planners and grid operators requiring them to participate in frequency control tasks. To address this demand, a large number of wind frequency control proposals have been reported in the literature in recent years. Many of these solutions have been tested by speciﬁc experiments carried out in computer simulation environments. This paper proposes a methodology to evaluate the effects of enabling frequency support by wind turbines on the magnitudes that characterize the inertial response of a power system by using an analytical approach. The derived formulation and the illustrations are designed to provide a better understanding of both the mechanisms that determine the frequency stability indices and the improvement achieved by enabling the inertial response of wind turbines by implementing a virtual inertia-based method on the active power controllers of these machines. To facilitate the comprehension of the results obtained, the analytical approach is complemented with time-domain simulations in a predeﬁned test system implemented in MATLAB/Simulink ® . The proposed methodology achieves a generalization of the results and can be used for the assessment of any power system conﬁguration. linearized models of the speed-governor of the synchronous generation, linearization of the frequency control characteristic of wind turbines, disregard of the dynamics of wind turbine power converters and non-linearities, among others. All these practices are well justiﬁed in the literature, and they have been indicated in the corresponding part of the article’s body.


Introduction
The imminent depletion of fossil fuel reserves, the growing concern for the preservation of the environment and the constant increase in energy needs are some of the factors that have motivated the use of renewable energy sources for production of electrical energy in recent decades. From the various alternatives available for the use of non-conventional renewable energy resources, wind generation has been making its way in the electricity generation market, so that, nowadays, it presents a significant share of participation in the energy mix of many power systems in the world and a high level of technological maturity [1,2]. Among the factors that have motivated this trend in the use of this generation technology are: the cost of the energy produced, the annual availability of the resource, the efficiency in energy conversion, the geographical areas required for its exploitation and the flexibility to respond to the technical requirements of the grid operator. However, the variable nature of wind and the uncertainty associated with its prediction make wind generation management complex, compared to conventional electricity generation, making its integration into the grid somewhat difficult, even more so when the level of its participation becomes significant [3,4].
One of the main challenges the grid operator must deal with in its effort to guarantee the quality, stability and reliability standards of the system, is related to frequency control, system implemented in EMTP-RV simulator. A methodology for the analysis of frequency dynamics in large-scale power systems with a high level of wind energy penetration by means of a simplified model for wind generators and load frequency control-scheme is presented by the authors in [21]. In this work, an aggregated representation of the Spanish bulk system is used as a test bench in MATLAB/Simulink ® simulation tool. In addition, by using a similar approach to a previous work, a novel frequency sensitive-based virtual inertia control to extract the kinetic energy of the wind turbine and stored energy from the DC-link capacitor for short-term grid frequency regulation is proposed in [22]. There, the MATLAB/Simulink ® simulation results are complemented with hardware-in-the-loop simulations implemented in OPAL-RT software. An exhaustive assessment of the impact on the grid by wind power short-term frequency support by taking advantage of pitch angle de-loading [23], kinetic energy extraction [24], and wind turbine over-speeding is investigated [25]. A common aspect that characterizes all these research papers is the use of MATLAB/Simulink ® as a simulation tool, except for [25] which assessed the effectiveness of its proposal in the IEEE-39 bus system implemented in DIgSILENT/PowerFactory software. In all the works described here, very promising results are obtained that justify the participation of wind turbines in grid frequency control tasks. Conclusions from these works are based on both simulation and experimental tests under certain operating scenarios and controlled environments, designed to accurately evaluate the benefits of the proposed control strategies. To address the study of frequency stability issues in power systems, two main different approaches can be found in the literature: one based on the concept of load frequency control (LFC) [26,27], and the other based on detailed computational simulations, sometimes including hardware-in-the-loop [28], or real-time simulations [29]. Each approach has its capabilities and advantages, although one can complement the other [30]. The literature review presented here shows that there are innovative proposals that successfully address the problem of high-level penetration of wind generators when these participate in the frequency control tasks. Various computer simulation programs are used to evaluate the proposals on a test-bench designed for the purpose of applying different approaches regarding the detailed representation of the system under analysis. From this review, it has been found that: (1) none of these studies offers a detail of how the corrective actions of the frequency control support by wind generators might impact on the parameters that define the inertial characteristics of the system as a whole, and (2) the proposals for participation of wind generators in frequency control are evaluated on specific operating scenarios under controlled environments.
In the first case, the cited papers conclusions are based on the final response of the system when the frequency control actions of the wind turbines are implemented. Although, in these results the reduction of frequency nadir and rate of change of frequency (RoCoF) when enabling frequency support by wind turbines is evident, it is not possible to identify which specific characteristic of the implemented frequency control strategy in wind turbines is responsible for improving both the equivalent inertia constant as the damping constant of the system. Finding this relationship would lead to a better understanding of the effect of the implementation of fast-frequency response by wind turbines on the short-term operation of the system.
In the second case, the design, modeling, and simulation of a specific test-bench to evaluate the performance of the wind-frequency control strategies, provide a restricted vision of the problem, since it would be necessary to carry out as many simulations as so many scenarios need to be represented with different operating conditions (for example, various wind generation penetration levels and incoming wind speeds) in order to obtain more general results.
These research gaps have motivated the development of our proposal, which consists of: 1.
The proposal of a methodology for deriving the LFC-transfer function of a weak and isolated power system that incorporates the participation of wind turbines in frequency control, which is based on the theory of small-signal analysis. With the derived LFC-transfer function, the authors design an analysis tool devised for studying the frequency response of the power system under contingency by considering all the possible scenarios characterized by different wind penetration levels and wind speed conditions. The numerical evaluation of the LFC-transfer function to get the main parameters that characterize the inertial response of the whole system, instead of implementing a specific test system in simulation, make it possible to represent a vast number of contingency scenarios without these having to run one at a time in the simulation environment.

2.
With the numerical results generated, an inertial chart is designed to show the locus of the values of the main frequency response parameters: frequency nadir and timederivative of frequency as a function of the equivalent constant of inertia of the system. The resulting chart provides an outlook of the status of the frequency resilience of the system when wind generation is incorporated in substitution of the synchronous one.

3.
Based on the symbolic representation of the derived transfer function, this work proposes a method to identify which specific inertial parameters of the system (equivalent inertia constant or damping constant) the frequency response by wind turbines really impacts. This is done with the purpose of providing a comprehensive understanding of the corrective action provided by wind generation during its contribution to smoothing grid frequency fluctuations.
This work is organized as follows: Section 2 presents the conceptual framework and the analytical approach to the problem of primary frequency control in a power system and introduces a case study. As an example of a frequency control strategy for wind turbines, Section 3 explains the theoretical foundation and mechanics of the extended optimized power point tracking (OPPT) method [21] to allow for a better understanding of Section 4, which constitutes the main part of the work. Section 4 proposes a methodology for the assessment of the effects on the magnitudes that characterize the inertial response of the whole power system when wind turbines contribute to frequency support tasks. The application of the methodology is illustrated with the implementation of the OPPT method. In addition, the most relevant results are presented here and preliminary conclusions are established. Finally, Section 5 presents the lessons learned in this study. Figure 1 illustrates the typical time scale of different physical processes involved in the operation and control of a power system [31][32][33][34]. Since the most relevant magnitudes in a frequency control study, i.e., maximum frequency deviation, maximum rate of change of frequency, and frequency deviation in quasi-steady state, take place a few seconds after a disturbance appears, the modeling and analysis efforts used in this work are focused on the time range that characterizes the primary frequency control (see the highlighted time frame in the figure).

Load Frequency Control Scheme
Within the time frame of interest shown in Figure 1, it is common practice to represent the dynamics of the power system using the LFC approach, in which the dynamics of the mechanical variables of the prime movers of conventional synchronous generation predominates over the dynamics of electromagnetic variables, whose time constants are significantly lower [35,36]. In addition, it can be verified in the same figure that the typical dynamics of the power converter interfacing the wind turbine to the grid is much faster than the mechanical dynamics of the wind turbine. From the perspective of the time frame of interest in this study, the practically instantaneous response of power converters usually leads to their dynamics being ignored in the analytical approach to modeling [33,37]. For that reason, this work focuses its efforts on modeling the power systems components within the time window shaded in Figure 1.

Load Frequency Control Scheme
Within the time frame of interest shown in Figure 1, it is common practice to represent the dynamics of the power system using the LFC approach, in which the dynamics of the mechanical variables of the prime movers of conventional synchronous generation predominates over the dynamics of electromagnetic variables, whose time constants are significantly lower [35,36]. In addition, it can be verified in the same figure that the typical dynamics of the power converter interfacing the wind turbine to the grid is much faster than the mechanical dynamics of the wind turbine. From the perspective of the time frame of interest in this study, the practically instantaneous response of power converters usually leads to their dynamics being ignored in the analytical approach to modeling [33,37]. For that reason, this work focuses its efforts on modeling the power systems components within the time window shaded in Figure 1. Figure 2 shows the typical representation of the LFC scheme for the case of n operating synchronous generators. The equivalent inertia constant of the synchronous machines is represented by Heq, and is calculated by means of Equation (1), where Hi is the individual inertia constant of each generator and wi is the participation factor, defined as the ratio of the nominal capacity of each generation unit to a certain base power. The introduction of this factor in the scheme is very useful when there is a need to represent different levels of participation of different generation agents, both in demand coverage and frequency support. In Figure 2, ΔPmi is the mechanical power deviation experienced by the speed control system of each generator in response to a deviation of its angular speed, which, for the time frame of interest, is common to all synchronous machines, Δωeq. It is worth noting that this simplification of a unique common frequency for the whole system drastically reduces the integration step in simulations and obtains analytical expressions for describing the dynamic behavior of the power system. The deviation of power demand in the system is denoted by ΔPL. Additionally, the damping constant of the supplied load is represented by Deq. Ri is the droop of each speed controller, whose dynamics is modeled by a transfer function TCi(s).
To evaluate the frequency response of the system, Δωeq, in the event of a disturbance, ΔPL, the following transfer function is derived from Figure 2:   Figure 2 shows the typical representation of the LFC scheme for the case of n operating synchronous generators. The equivalent inertia constant of the synchronous machines is represented by H eq , and is calculated by means of Equation (1), where H i is the individual inertia constant of each generator and w i is the participation factor, defined as the ratio of the nominal capacity of each generation unit to a certain base power. The introduction of this factor in the scheme is very useful when there is a need to represent different levels of participation of different generation agents, both in demand coverage and frequency support. In Figure 2, ∆P mi is the mechanical power deviation experienced by the speed control system of each generator in response to a deviation of its angular speed, which, for the time frame of interest, is common to all synchronous machines, ∆ω eq . It is worth noting that this simplification of a unique common frequency for the whole system drastically reduces the integration step in simulations and obtains analytical expressions for describing the dynamic behavior of the power system. The deviation of power demand in the system is denoted by ∆P L . Additionally, the damping constant of the supplied load is represented by D eq . R i is the droop of each speed controller, whose dynamics is modeled by a transfer function TC i (s).

Case Study
A case study for evaluating the effects of the frequency support provided by speed-wind turbines on the inertial response of a given power system for faci matches between generation and demand is presented in Figure 3. It represents To evaluate the frequency response of the system, ∆ω eq , in the event of a disturbance, ∆P L , the following transfer function is derived from Figure 2: This equation represents the LFC-transfer function of the power system.

Case Study
A case study for evaluating the effects of the frequency support provided by variable speed-wind turbines on the inertial response of a given power system for facing mismatches between generation and demand is presented in Figure 3. It represents a power system (with a rated frequency of 50 Hz) dominated by thermal and wind generation within a frequency control area. This is a typical case of many island electric systems, in which the excursions of the frequency (either under normal conditions, due to intermittent injection of power from non-conventional renewable energy sources, or in contingency situations, due to the sudden loss of generation units or consumption centers) are very pronounced compared to a continental system [38,39]. In the example illustrated in Figure 3, the thermal generation is represented by the linearized model of the prime mover of a steam turbine and its speed governor [35]. The input, ∆ω eq , is applied to a droop characteristic, R 1 , in order to generate signals to be applied to the first-order turbine model with time constant τ G (stage of admission of steam to the turbine). Then, the resulting signal is input to a block that represents the dynamics of each turbine stage: from high to lower pressure stages with contribution gain F HP and the time constants T HP , T RH , and T CH . More details of this tandem-compound configuration of the steam turbine can be consulted in [35]. Finally, the complete thermal generation model generates a control signal, ∆P m1 , intended to compensate the power mismatch in the LFC-scheme.  In this example, the thermal unit has an inertia constant H1 = 3.6 s and a participation factor w1 = 0.8 (where the base power is the total generation capacity of the system). Accordingly, the participation of wind generation is 20% (wWT = 0.2). It is also assumed that the wind speed is constant, v = 9.6 m/s, and that the active power of wind turbines is dispatched using the Maximum Power Point Tracking (MPPT) criterion. For the frequency control study, a contingency event characterized by a sudden load increase (for example, due to the loss of part of the operating generation) ΔPL = 1% is analyzed. This situation  In the same figure, a block can be seen implementing the mathematical model used to represent the dynamics of the wind generation (variable speed-wind turbine, VSWT) within the LFC scheme. This model, which can be consulted in detail in [40], outputs an active power deviation signal, ∆P g , according to the incoming wind speed v. By default, the wind turbine model dispatches its active power according to an optimal power conversion criterion, and it is not expected that it provides frequency support to the system in any case, unless a wind-frequency control feature is enabled. The latter is discussed extensively in Section 3. The values of all the parameters of these two models are indicated in the Appendices A and B.
In this example, the thermal unit has an inertia constant H 1 = 3.6 s and a participation factor w 1 = 0.8 (where the base power is the total generation capacity of the system). Accordingly, the participation of wind generation is 20% (w WT = 0.2). It is also assumed that the wind speed is constant, v = 9.6 m/s, and that the active power of wind turbines is dispatched using the Maximum Power Point Tracking (MPPT) criterion. For the frequency control study, a contingency event characterized by a sudden load increase (for example, due to the loss of part of the operating generation) ∆P L = 1% is analyzed. This situation provokes a grid sub-frequency, through which it is possible to evaluate the inertial response of the system and the effectiveness of the frequency recovery tasks carried out by the synchronous generation. Figure 4 shows some results from simulations (with MATLAB/Simulink ® ) of the LFC scheme of Figure 3 under the operational situation previously described. It illustrates the time domain behavior of the following variables of interest: deviation of the mechanical power experienced by the primary mover of the synchronous generator, ∆P m1 , active power dispatched by the wind turbine, P g , frequency of the system, f s , and the absolute value of its time derivative, |ROCOF|.

Analytical Approach to the Problem
The simple experiment carried out in the previous subsection has been designed t The simulation results show a typical behavior of these variables in response to a sudden power mismatch: the grid frequency (defined by the synchronous generator speed, Figure 4c) suffers a decay that triggers the frequency recovery task, carried out by the speed controllers of thermal groups (Figure 4a), leading to a maximum deviation, ∆f max = 0.083 Hz, and eventually a new equilibrium point with a quasi-steady deviation ∆f qs = 0.03 Hz. Figure 4d shows the evolution of the absolute value of the rate of change of frequency (RoCoF), where its maximum value is |ROCOF| max = 0.0837 Hz/s. Since the power injected into the grid by the wind turbine ( Figure 4b) depends solely on the incident wind speed (v = 9.6 m/s), which is assumed constant in this study, no response is expected, in terms of power, to contribute to the frequency recovery process. This situation highlights one of the challenges posed by the incorporation of wind power displacing conventional synchronous generation, as the latter is responsible for guaranteeing the frequency stability of the grid. To maintain the inertial characteristics of the system, wind turbines need to contribute to mitigating these negative effects by providing an effective inertial response.

Analytical Approach to the Problem
The simple experiment carried out in the previous subsection has been designed to warn the reader about the need to involve non-synchronous generators (e.g., variable speed-wind turbines), whose integration in an electrical system implies a displacement of conventional synchronous generation that is detrimental to the equivalent inertial response of the system. In order to broaden the results obtained up to this point and to draw more general conclusions, this problem is now approached analytically, taking the transfer function of the previous power system as a starting point. This function is presented in Equation (3), obtained from the definition expressed in Equation (2) and particularized for the system outlined in Figure 3. As it can be seen in the resulting expression, the response of the system frequency, ∆ω s , to a disturbance, ∆P L , is by no means conditioned by any parameter or variable related to the incorporated wind generation. This is because, as demonstrated in the previous example, the criterion with which variable speed-wind turbines deliver their active power makes a natural inertial response from this machine impossible. Therefore, it is to be expected that the higher the level of participation of these generators, the poorer the inertial response that the resulting conventional synchronous generation can provide to face power imbalances.
From the expression in Equation (3), it is possible to systematically analyze the frequency response of the system to a certain disturbance with different levels of participation of synchronous generation in frequency control tasks. For example, the response to a ∆P L = 0.01 pu step (as defined in Section 2.2) is evaluated below for all possible values of w 1 . The numerical calculation of the step response is focused on the evaluation of the following quantities: • Peak value of the step response, max{∆ω s (t)}: for calculation of the maximum frequency deviation, ∆f max . • Steady state value of the step response, lim t→∞ {∆ω s (t)}: for calculation of the quasisteady state frequency deviation, ∆f qs . • Maximum absolute value of d∆ω s /dt: for calculation of |RoCoF| max .
From this numerical procedure, the set of curves illustrated in Figure 5 is obtained, where the different participation factors, w 1 , are expressed in terms of the equivalent inertia constant, H eq , normalized for the base case (w 1 = 0.8). In this example, 100 possible values of w 1 have been used to evaluate Equation (3), which means that the generated chart represents an equal number of simulated scenarios. In this chart, each simulation result is represented by one round marker placed over the corresponding curve. For this example, the chart contains a total of 300 numerical results (∆f max , ∆f qs , and |ROCOF| max ) for 100 simulated scenarios. If the user requires a chart with a higher resolution, a finer variation of the values of participation factors, w 1 , must be done. The inertial chart of Figure 5 allows for visualizing the effect of conventional synchronous generation displacement on the parameters that characterize the inertial response of the system within the time frame of interest. As a fact, the inertial characteristics of the system under study tend to degrade as operating synchronous generation is displaced. Therefore, non-synchronous generators should participate in the frequency recovery tasks in a similar way as the displaced conventional synchronous units did. In Figure 5 it can be seen that, for the base case (in which thermal generation is 80% of the total installed capacity), the numerical results obtained from evaluation of Equation (3) are very close to those achieved by simulation in Figure 4.
values of w1 have been used to evaluate Equation (3), which means that the ge chart represents an equal number of simulated scenarios. In this chart, each sim result is represented by one round marker placed over the corresponding curve. example, the chart contains a total of 300 numerical results (∆fmax, ∆fqs, and |ROC for 100 simulated scenarios. If the user requires a chart with a higher resolution variation of the values of participation factors, w1, must be done. The inertial char ure 5 allows for visualizing the effect of conventional synchronous generation d ment on the parameters that characterize the inertial response of the system wi time frame of interest. As a fact, the inertial characteristics of the system under stu to degrade as operating synchronous generation is displaced. Therefore, non-s nous generators should participate in the frequency recovery tasks in a similar wa displaced conventional synchronous units did. In Figure 5 it can be seen that, for case (in which thermal generation is 80% of the total installed capacity), the nu results obtained from evaluation of Equation (3) are very close to those achieved ulation in Figure 4.

Method Based on the Concept of Virtual Inertia
Making use of the significant amount of inertial resources available in a wind several works (e.g., [7,[21][22][23]) address the issue of providing an effective inertial r from this machine, so that it participates in the grid frequency control. This is a through additional control loops that complement the main MPPT controllers of t turbine. These auxiliary control loops introduce additional active power set po pending on the system requirements, as a function of the frequency deviation t

Method Based on the Concept of Virtual Inertia
Making use of the significant amount of inertial resources available in a wind turbine, several works (e.g., [7,[21][22][23]) address the issue of providing an effective inertial response from this machine, so that it participates in the grid frequency control. This is achieved through additional control loops that complement the main MPPT controllers of the wind turbine. These auxiliary control loops introduce additional active power set points, depending on the system requirements, as a function of the frequency deviation that may exist, thereby enabling an inertial response from the wind turbine. At this point, it must be taken into account that this inertial response would be forced, controlled, virtual, or synthetic, but energetically supported by a real physical property of the wind turbine: its natural inertia.
To fulfill the central objective of this work, the virtual inertia control strategy proposed by the authors in [21], called the extended optimized power point tracking (OPPT) method, whose mechanics are explained below, has been considered. Generally, the active power dispatched by variable speed wind turbines follows a set point signal, P MPPT , which is applied to the corresponding control system implemented in the machine. Equation (4) shows the definition of this signal for the optimization zone, where: f Kopt is a multiplication factor (equal to unity, by default), K opt is the optimization constant of the turbine (which depends on aerodynamic and constructive characteristics of the rotor), and ω g is the shaft speed. If this expression is plotted on the P-ω plane of the wind turbine, taking into account the technical limits of the machine, the locus shown in Figure 6 (blue line) is obtained, known as the MPPT curve. The control function of the extended OPPT method is shown in Equation (5), where: ω g0 is the rotor speed prior to the requirement of frequency support, k vir is a design constant of the method, W vir is a weighting factor, and H WT is the inertia constant of the wind turbine. This expression, derived in detail in [21], has two components that determine the magnitude of the wind inertial response sought: one dependent on the grid frequency deviation, ∆f s , and another dependent on the time derivative of frequency, df s /dt. where exist, thereby enabling an inertial response from the wind turbine. At this point, it must be taken into account that this inertial response would be forced, controlled, virtual, or synthetic, but energetically supported by a real physical property of the wind turbine: its natural inertia.
To fulfill the central objective of this work, the virtual inertia control strategy proposed by the authors in [21], called the extended optimized power point tracking (OPPT) method, whose mechanics are explained below, has been considered. Generally, the active power dispatched by variable speed wind turbines follows a set point signal, PMPPT, which is applied to the corresponding control system implemented in the machine. Equation (4) shows the definition of this signal for the optimization zone, where: fKopt is a multiplication factor (equal to unity, by default), Kopt is the optimization constant of the turbine (which depends on aerodynamic and constructive characteristics of the rotor), and ωg is the shaft speed. If this expression is plotted on the P-ω plane of the wind turbine, taking into account the technical limits of the machine, the locus shown in Figure 6 (blue line) is obtained, known as the MPPT curve.
[pu] Kopt opt g P f K (4) Figure 6. Mechanics of the method based on virtual inertia control.
The control function of the extended OPPT method is shown in Equation (5), where: ωg0 is the rotor speed prior to the requirement of frequency support, kvir is a design constant of the method, Wvir is a weighting factor, and HWT is the inertia constant of the wind turbine. This expression, derived in detail in [21], has two components that determine the magnitude of the wind inertial response sought: one dependent on the grid frequency deviation, Δfs, and another dependent on the time derivative of frequency, dfs/dt.
To explain the inertia emulation process of the wind turbine achieved by implementing the extended OPPT method (whose implementation scheme in the active power control system of the wind turbine is illustrated in Figure 7), let us return to Figure 6. Starting from the hypothesis that the wind turbine is driven by a constant wind speed (v = 9.6 m/s), under normal operating conditions, the generator will deliver an active power PA, associated with a speed ωg0, according to the MPPT control. If, for example, the power system experiences an underfrequency event caused by a sudden decompensation between generated and consumed power, the presence of non-zero values of Δfs and dfs/dt in Equation (5) will cause the corresponding implementation scheme to start the inertial response emulation process by moving the MPPT curve upwards, in order to drive the operation of the wind turbine towards point B of its P-ω characteristics. To explain the inertia emulation process of the wind turbine achieved by implementing the extended OPPT method (whose implementation scheme in the active power control system of the wind turbine is illustrated in Figure 7), let us return to Figure 6. Starting from the hypothesis that the wind turbine is driven by a constant wind speed (v = 9.6 m/s), under normal operating conditions, the generator will deliver an active power P A , associated with a speed ω g0 , according to the MPPT control. If, for example, the power system experiences an underfrequency event caused by a sudden decompensation between generated and consumed power, the presence of non-zero values of ∆f s and df s /dt in Equation (5) will cause the corresponding implementation scheme to start the inertial response emulation process by moving the MPPT curve upwards, in order to drive the operation of the wind turbine towards point B of its P-ω characteristics.
At this point, the value adopted by P B will depend on the severity of the frequency event, measured in terms of the magnitude of ∆f s and df s /dt. Under this operational situation, the value of the active power extracted from the electric generator will exceed the value of the mechanical power delivered by the turbine, causing an internal power imbalance, which will lead to the deceleration of the rotating masses of the wind turbine. This use of kinetic energy, seen as a decrease in the speed of the turbine-generator set, moves the operating point from B to C, following the path imposed by the MPPT curve shifted to the left. At point C, the active power produced by the generator, P C , equals again the mechanical power given by the turbine, at a new equilibrium condition. The excursion of the active power of the wind turbine along the path A→B→C makes its contribution to the frequency support of the system feasible. Once this support is finished, the method forces the turbine to recover an optimal operating condition (point A), thus preparing it to meet future requirements. Similarly, in case of an over-frequency event in the system, the inertial support provided by the turbine will follow the path A→D→E→A by shifting the MPPT curve downwards.  Figure 8 shows the simulation results of implementing the scheme of Figure 7 in the wind turbines that integrate the power system under study. The benefits of using the kinetic energy available in a wind turbine to contribute to the grid frequency support can be appreciated. First, the inertial response of wind generation to the frequency drop (Figure 8b) makes the mechanical power control action from thermal generation somewhat smoother and of a lower magnitude (Figure 8a). The combined control action of wind and thermal generation in frequency support gives, in this case, a 24% reduction of the maximum frequency deviation, compared to the base case (Figure 8c). Additionally, a 19% reduction in the maximum magnitude of RoCoF is obtained (Figure 8d). Here, it is important to note that the results presented are due to a particular wind generation situation. The effectiveness of the virtual inertia algorithm (extended OPPT-method) in reducing ∆f max and |RoCoF| max will depend, on the one hand, on the magnitude of the disturbance (measured in terms of the state of ∆f s and df s /dt), and on the other hand, on the level of kinetic energy stored in the machine before the request for the frequency support service (see the influence of the term ω g0 in Equation (5)). Since the shaft speed, ω g0 , is determined by pre-existing wind conditions, a more general analysis is required. This is discussed in detail in the next section. At this point, the value adopted by PB will depend on the severity of the frequency event, measured in terms of the magnitude of Δfs and dfs/dt. Under this operational situation, the value of the active power extracted from the electric generator will exceed the value of the mechanical power delivered by the turbine, causing an internal power imbalance, which will lead to the deceleration of the rotating masses of the wind turbine. This use of kinetic energy, seen as a decrease in the speed of the turbine-generator set, moves the operating point from B to C, following the path imposed by the MPPT curve shifted to the left. At point C, the active power produced by the generator, PC, equals again the mechanical power given by the turbine, at a new equilibrium condition. The excursion of the active power of the wind turbine along the path A→B→C makes its contribution to the frequency support of the system feasible. Once this support is finished, the method forces the turbine to recover an optimal operating condition (point A), thus preparing it to meet future requirements. Similarly, in case of an over-frequency event in the system, the inertial support provided by the turbine will follow the path A→D→E→A by shifting the MPPT curve downwards.  Figure 8 shows the simulation results of implementing the scheme of Figure 7 in the wind turbines that integrate the power system under study. The benefits of using the kinetic energy available in a wind turbine to contribute to the grid frequency support can be appreciated. First, the inertial response of wind generation to the frequency drop (Figure 8b) makes the mechanical power control action from thermal generation somewhat smoother and of a lower magnitude (Figure 8a). The combined control action of wind and thermal generation in frequency support gives, in this case, a 24% reduction of the maximum frequency deviation, compared to the base case (Figure 8c). Additionally, a 19% reduction in the maximum magnitude of RoCoF is obtained (Figure 8d). Here, it is important to note that the results presented are due to a particular wind generation situation. The effectiveness of the virtual inertia algorithm (extended OPPT-method) in reducing kinetic energy stored in the machine before the request for the frequency support service (see the influence of the term ωg0 in Equation (5)). Since the shaft speed, ωg0, is determined by pre-existing wind conditions, a more general analysis is required. This is discussed in detail in the next section.

Linearization of the Extended OPPT Method
In order to study the effects that frequency support from wind generation has on the inertial characteristics of the system, similar to the analysis presented in Section 2.3, it is necessary to include the particular control function of the frequency control method used, in this case the extended OPPT one in Equation (5), in the transfer function of Equation (3) that defines the dynamics of the system within the time frame of interest to this work. To get a rational transfer function, in the first place, it is necessary to carry out a linearization process.
Considering the term dependent on Δfs in Equation (5): Since, in frequency stability studies, small disturbances are considered (that is, Δfs ≈ 0), it is possible to linearize the previous expression as follows:

Linearization of the Extended OPPT Method
In order to study the effects that frequency support from wind generation has on the inertial characteristics of the system, similar to the analysis presented in Section 2.3, it is necessary to include the particular control function of the frequency control method used, in this case the extended OPPT one in Equation (5), in the transfer function of Equation (3) that defines the dynamics of the system within the time frame of interest to this work. To get a rational transfer function, in the first place, it is necessary to carry out a linearization process.
Considering the term dependent on ∆f s in Equation (5): Since, in frequency stability studies, small disturbances are considered (that is, ∆f s ≈ 0), it is possible to linearize the previous expression as follows: Energies 2021, 14, 3660

of 22
Then, as explained in Section 3.1, an increase in the active power dispatched by the wind turbine, by shifting the MPPT curve upwards (from point A to point B in Figure 6), can be formulated as follows: By substituting Equation (7) in Equation (8), and making ∆f s[pu] = ∆ω s[pu] : Now, the term dependent on df s /dt in Equation (5) is obtained from the definition in [21]: With some algebra in Equation (10) and expressing ω s in terms of its deviation from the nominal value ω s0 : Considering small disturbances (∆ω s ≈ 0): ∆ω 2 s << 2∆ω s ω s0 , and the previous equation can be reduced to: Then, adding Equations (9) and (12) to obtain the control function of the linearized extended OPPT method and taking into account the weighting factor of the frequency derivative dependent term, W vir , in the Laplace domain: Or, in compact form: Next, in order to evaluate the numerical accuracy of the linearization process with respect to the control method outlined in Figure 7 and implemented on the diagram in Figure 3, the control function in Equation (14) is included in the transfer function of the power system in Equation (3), finally being: Synchronous generation speed governor Wind turbine frequency control (extended OPPT − method) (15) where: α = α s 1/T wo + s At this point, it is essential to discuss the scalability and the complexity of the developed model in larger power systems with their generation units disaggregated. In order to obtain a general equation for implementing the proposal on a multimachine power system, let us take the complete transfer function in Equation (15) and pay particular attention to the denominator terms. The contribution to the grid-frequency support from synchronous generator and wind turbine, SG 1 (s) and WT 1 (s) terms, respectively, can be extracted from Equation (15) and defined as follows: According to the representation of the LFC-scheme for a multimachine power system illustrated in Figure 2, and expressed in Equations (1) and (2), the denominator term corresponding to the synchronous generation speed governor model and its inertia in Equation (16) can be extended to the case of n-synchronous generation units, as below: where, H i , R i , w i and TC i are the parameters and functions of the ith synchronous generator unit (similar to those introduced in Section 2.1). Now, if we consider the denominator term corresponding to the wind turbine equipped with a frequency control characteristic (extended OPPT-method, in this study), as in Equation (17), it can be extended to the case of N-wind turbines as: where, α' i , β' i , and w WTi are functions and parameters of the ith wind turbine that incorporates a frequency response characteristic. Note that, according to the definition of α and β in Equation (13), α i and β i will depend on the pre-disturbance rotor speed ω g0i , which is a function of the incoming wind speed to the ith wind turbine, v i (see Appendix D). If the user is interested in disabling the participation of a specific wind turbine in frequency support, it is enough to make the corresponding functions α' i and β' i equal to zero. Finally, the general form of the transfer function for the implementation of the proposed method in case of multimachine power systems is: Now, returning to our case of study, note that, in Equation (15), the transfer functions of the high pass and low pass filters (from the ∆f s and df s /dt dependent loops of Figure 7) have been included. The final form of this transfer function is the one that was sought to fulfill the main objective of this research.
Finally, the control action of the extended OPPT method can be expressed as follows: It is important to note that, since the linearized model, which represents the inertial response by the wind generation, that has been incorporated into the transfer function of the power system Equation (15) does not take into account the dynamics of the wind turbine's mechanical system, its accuracy will be limited to a time window of the order of the wind turbine's natural inertia constant, H WT (which for the model considered, corresponds to a value of 5.29 s). This can be verified in Figure 9, where the inertial response offered by the wind turbine for the studied contingency scenario (∆P L = 1%) has been plotted using two approaches: the exact one, by implementing Equation (5) in the electromechanical model of the wind turbine (Figure 3), and the linearized one, by means of the direct application of Equation (21) in the LFC-transfer function of the power system (Equation (15)). The values for the numerical evaluation of Equation (21) can be found in the Appendices B-D. ' ' It is important to note that, since the linearized model, which represents the inertia response by the wind generation, that has been incorporated into the transfer function o the power system Equation (15) does not take into account the dynamics of the wind tur bine's mechanical system, its accuracy will be limited to a time window of the order of th wind turbine's natural inertia constant, HWT (which for the model considered, correspond to a value of 5.29 s). This can be verified in Figure 9, where the inertial response offere by the wind turbine for the studied contingency scenario (∆PL = 1%) has been plotted usin two approaches: the exact one, by implementing Equation (5) in the electromechanica model of the wind turbine (Figure 3), and the linearized one, by means of the direct appl cation of Equation (21) in the LFC-transfer function of the power system (Equation (15) The values for the numerical evaluation of Equation (21) can be found in the Appendixe B-D. Note that the error when using the linearized approach does not exceed 0.4% in an case. This is reflected in the high numerical accuracy of the results in terms of frequenc and its derivative ( Figure 10  Note that the error when using the linearized approach does not exceed 0.4% in any case. This is reflected in the high numerical accuracy of the results in terms of frequency and its derivative ( Figure 10  ' ' It is important to note that, since the linearized model, which represents the inertia response by the wind generation, that has been incorporated into the transfer function o the power system Equation (15) does not take into account the dynamics of the wind tur bine's mechanical system, its accuracy will be limited to a time window of the order of the wind turbine's natural inertia constant, HWT (which for the model considered, correspond to a value of 5.29 s). This can be verified in Figure 9, where the inertial response offered by the wind turbine for the studied contingency scenario (∆PL = 1%) has been plotted using two approaches: the exact one, by implementing Equation (5) in the electromechanica model of the wind turbine (Figure 3), and the linearized one, by means of the direct appli cation of Equation (21) in the LFC-transfer function of the power system (Equation (15)) The values for the numerical evaluation of Equation (21) can be found in the Appendixe B-D.

Evaluation of the Effect of Wind Generation Contribution to Frequency Support on the Inertial Characteristics of the System
After evaluating the numerical accuracy of the approximate approach, this section follows a procedure similar to that described in Section 2.3 (which led to the results of Figure 5) to evaluate the transfer function in Equation (15), considering the same contingency event for different incident wind conditions and different values of the participation factor, w 1 . As a result, the inertial charts in Figure 11 are obtained, which show the effect of the displacement of the spinning synchronous generation on the ∆f max and |RoCoF| max indicators when the contribution of wind generators to frequency support is enabled. Since 100 different values of w 1 were assigned, each of the curves represents 100 simulated scenarios. This figure shows the corrective action of wind turbines if the extended OPPT method is implemented. The participation of these generation units in the frequency recovery contributes to reducing the level of degradation of the inertial characteristics of the system as a whole, the benefit being much more pronounced, the higher the wind speed is. The fact that the increase in wind generation penetration has the effect of shifting the curves upwards means that the contribution of wind generators to frequency support makes their integration (displacing synchronous generators) less harmful, in terms of the parameters of interest: ∆f max and |RoCoF| max .

Evaluation of the Effect of Wind Generation Contribution to Frequency Support on the Inertial Characteristics of the System
After evaluating the numerical accuracy of the approximate approach, this section follows a procedure similar to that described in Section 2.3 (which led to the results of Figure 5) to evaluate the transfer function in Equation (15), considering the same contingency event for different incident wind conditions and different values of the participation factor, w1. As a result, the inertial charts in Figure 11 are obtained, which show the effect of the displacement of the spinning synchronous generation on the Δfmax and |RoCoF|max indicators when the contribution of wind generators to frequency support is enabled. Since 100 different values of w1 were assigned, each of the curves represents 100 simulated scenarios. This figure shows the corrective action of wind turbines if the extended OPPT method is implemented. The participation of these generation units in the frequency recovery contributes to reducing the level of degradation of the inertial characteristics of the system as a whole, the benefit being much more pronounced, the higher the wind speed is. The fact that the increase in wind generation penetration has the effect of shifting the curves upwards means that the contribution of wind generators to frequency support makes their integration (displacing synchronous generators) less harmful in terms of the parameters of interest: Δfmax and |RoCoF|max. To better understand the corrective effect that the implementation of the wind turbine frequency control strategy has on the inertial response of the system, one can take a look at the denominator of the transfer function in Equation (15). By regrouping terms, it is possible to write the following expressions: From Equations (22) and (23) it is inferred that the loop of the extended OPPT method dependent on dfs/dt directly influences the equivalent inertia constant of the system, HT (influence of the term β'), whereas the loop dependent on Δfs affects the equivalent damping constant of the system, DT (effect of the term α'). Therefore, when the extended OPPT To better understand the corrective effect that the implementation of the wind turbine frequency control strategy has on the inertial response of the system, one can take a look at the denominator of the transfer function in Equation (15). By regrouping terms, it is possible to write the following expressions: From Equations (22) and (23) it is inferred that the loop of the extended OPPT method dependent on df s /dt directly influences the equivalent inertia constant of the system, H T (influence of the term β'), whereas the loop dependent on ∆f s affects the equivalent damping constant of the system, D T (effect of the term α'). Therefore, when the extended OPPT control scheme starts operating during the frequency event, it forces the wind turbine to produce variations in the active power injected into the grid, in such a way that the system perceives a variation in its H T and D T parameters, although temporarily and while the control action implemented in the wind turbine lasts.
To provide a notion of this phenomenon in the time domain, the case study presented in Section 3.2 will be considered again. According to Equation (21), the dynamics of coefficients α' and β' can be calculated by using Equations (24) and (25), respectively, using the numerical results obtained when evaluating Equation (15).
Finally, by replacing the results of this calculation in expressions Equations (22) and (23), it is possible to obtain the dynamics of H T and D T (Figure 12), whose precision, once again, is limited to the time interval covered by H WT (the unshaded time window). Note in this figure that the participation of wind generation in the frequency support of the system has an amplifying effect of the synchronous equivalent constant of inertia of the system, H T , even as to double the pre-disturbance value at a specific moment. Theoretically, it is expected that, once the wind power participation has ended, H T will stabilize at a new value, defined only by the spinning synchronous generation.
Energies 2021, 14,3660 control scheme starts operating during the frequency event, it forces the wind tu produce variations in the active power injected into the grid, in such a way that the perceives a variation in its HT and DT parameters, although temporarily and w control action implemented in the wind turbine lasts.
To provide a notion of this phenomenon in the time domain, the case study pr in Section 3.2 will be considered again. According to Equation (21), the dynamics ficients α' and β' can be calculated by using Equations (24) and (25), respectively the numerical results obtained when evaluating Equation (15).
Finally, by replacing the results of this calculation in expressions Equations ( (23), it is possible to obtain the dynamics of HT and DT (Figure 12), whose precisio again, is limited to the time interval covered by HWT (the unshaded time window in this figure that the participation of wind generation in the frequency suppor system has an amplifying effect of the synchronous equivalent constant of inerti system, HT, even as to double the pre-disturbance value at a specific moment. T cally, it is expected that, once the wind power participation has ended, HT will stab a new value, defined only by the spinning synchronous generation. The main effect of the transient increase of HT on the inertial characteristic system is the reduction of |RoCoF|max with respect to the base case, as already sh Figures 8d and 11b. Something similar happens with the dynamics of DT. The mented virtual inertia method introduces a sudden and significant increase in the alent damping constant of the load, a fact that contributes to the reduction of th (Figures 8c and 11a). In this case, the tendency of DT to recover its original Deq v the wind turbine contribution to frequency control ends, is evident.

Conclusions
This paper proposes a methodology for the analysis of the effects of the contr The main effect of the transient increase of H T on the inertial characteristics of the system is the reduction of |RoCoF| max with respect to the base case, as already shown in Figures 8d and 11b. Something similar happens with the dynamics of D T . The implemented virtual inertia method introduces a sudden and significant increase in the equivalent damping constant of the load, a fact that contributes to the reduction of the ∆f max (Figures 8c and 11a). In this case, the tendency of D T to recover its original D eq value, as the wind turbine contribution to frequency control ends, is evident.

Conclusions
This paper proposes a methodology for the analysis of the effects of the contribution of wind generation to frequency control tasks on the inertial characteristics of power systems. This problem is approached from an analytical perspective through the study of the transfer function that represents the dynamics of the system frequency in response to power disturbances in the primary frequency control stage. This analytical approach allows the researcher to go from implementing a specific test system in a time-domain simulation to directly studying the results generated by the numerical evaluation of the wind-penetrated power system transfer function, in order to get the main parameters that characterize the inertial response of the whole system. In this way, it is possible to represent a vast number of contingency scenarios without these having to be run one at a time in the simulation environment, as has been a common practice in the literature. With these results, an inertial chart that represents the locus of the frequency nadir and its maximum time derivative, as a function of the wind-penetration level and the incoming wind speed, is constructed to show the improvement of the frequency recovery capability of the whole system offered by the participation of the wind generation in the frequency control task. As seen by the example plotted in Figure 11, 100 possible scenarios of wind generation penetration in the system were represented, each one evaluated for five different wind conditions and one base case (600 scenarios to quantify the frequency nadir and 600 to evaluate RoCoF; 1200 numerical results in total). These results were generated by the proposed analysis tool in a single run.
Furthermore, with the symbolic form of the derived transfer function it is possible to identify the influence of implementing virtual inertia in wind turbines on the equivalent system inertia constant and load damping constant (Equations (22) and (23)). Both parameters are decisive in characterizing the equivalent inertial response of the power system as a whole, measured in terms of the maximum frequency deviation and its maximum rate of change. In the particular case used as an example illustrating the proposed methodology, the results of this analysis ( Figure 12) show the benefits achieved by promoting the participation of wind generators in this task for the operation and control of a power system, reducing the negative impact of their massive integration in terms of an adequate frequency response to power imbalances: on the one hand, the equivalent constant of inertia of the system is transiently boosted (Figure 12a), passing from the pre-disturbance value 2.88 s to 5.22 s in the fast-frequency response stage. This improvement of the inertial characteristics caused by wind generation leads to a decrease in the frequency nadir by about 24% to the base case (Figure 8c). On the other hand, with the participation of wind turbine in the frequency control tasks, the instantaneous load damping constant is rapidly increased after the disturbance appears, passing from a steady-state value of 0.5 pu to 2.368 pu (Figure 12b). This effect explains the mitigation of the frequency change rate (Figure 8d), where the maximum absolute value of RoCoF reduces by about 19% if it is compared with the base case. Without Equations (22) and (23), derived through the proposed methodology, it would have been complicated to identify the influence of the fast-frequency response from wind turbines on the system parameters H eq and D eq . This application justifies the benefits of the analytical approach used for the design of the proposal.
It also constitutes an analytical tool that allows for obtaining general conclusions and having a first vision about the behavior of the system under contingency and how different control parameters affect this behavior. For example, it can be used as a tuning tool for the virtual inertia controllers of wind turbines. It can also be used to understand the interactions with other elements connected to the network through power electronics. Of course, it is a high-level analysis that needs to be complemented with detailed studies for each specific case.
Finally, it is important to note that the spirit of this paper is to present to the scientific community the formulation and benefits of the proposed methodology. Thus, in order to achieve a final model that can be run with a low computational burden without punishing the numerical accuracy of the tool, some simplifications have been introduced: use of linearized models of the speed-governor of the synchronous generation, linearization of the frequency control characteristic of wind turbines, disregard of the dynamics of wind turbine power converters and non-linearities, among others. All these practices are well justified in the literature, and they have been indicated in the corresponding part of the article's body.
Author Contributions: D.O.: conceptualization, methodology, software, validation, investigation, resources, data curation, writing-original draft, writing-review and editing, visualization, funding acquisition. S.M.: conceptualization, methodology, validation, resources, writing-review and editing, supervision, project administration, funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.