Sufficient Conditions for Robust Frequency Stability of AC Power Systems

This paper analyses the frequency stability of ac grids in the presence of non-dispatchable generation and stochastic loads. Its main goal is to evaluate conditions in which the system is robust to large, persistent active power disturbances without recurring to time-domain simulations. Considering the ongoing energy transition to more renewable sources, defining robustness boundaries is a key topic for power system planning and operation. However, much of the research on long-term studies has not dealt with robust dynamic constraints, while short-term analyses usually depend on time-consuming simulations to evaluate nonlinearities. To bridge this gap, the authors derive an algebraic equation that provides sufficient conditions for robust frequency stability in ac power systems and a relationship among four key quantities: the maximum active power perturbation, the minimum system damping, the steady-state and the transient frequency limits. To achieve this goal, it uses a nonlinear average-model of the ac grid and Lyapunov's direct method extended by perturbation analysis requiring only limited knowledge of the system parameters. The algebraic calculations are validated using time-domain simulations of the IEEE 39-bus test system and results are compared to the traditional Swing Equation model.


I. INTRODUCTION
A NY modern society requires an energy system that is affordable, accessible, secure and sustainable. However, recent technological advancements and social changes have been challenging the limits of this balance and pushing for a complete rethinking of energy systems. Electric power generation is now: 1) becoming more decentralized; 2) taking place closer to final users; 3) undergoing a technological transformation that comprehends more diversified primary sources, new forms of energy storage and power electronics [1], [2].
In this new paradigm, planning and operation of electrical ac grids are becoming more involved. For instance, conversion from primary sources and storage is performed using not only synchronous machines but also converter-interfaced generators (CIGs). The latter can supply up to 50% of the annual demand and 100% of the hourly demand in places such as Colorado, Denmark, Hawaii, Ireland, Tasmania, Texas and South Australia [3], [4]. This introduces new dynamics into the electric power system that are "yet to be fully understood" [5] and blurry the boundaries to which well-established models can be applied [6].
Moreover, groups of interconnected loads and distributed energy resources, also known as microgrids (MGs) [7], [8], can form islands and operate independently from the bulk power system. However, this desired feature is only achieved when primary sources, storage units, power converters and control systems are properly planned, designed and operated.
From this perspective, one of the main challenges is defining sufficient conditions for stable operation of low-inertia systems in islanded mode [2]. Stability here is seen in the sense of system states remaining within acceptable ranges for given changes in inputs, i.e. input-to-state stability. In recent years, this topic has been an active area of research and references [5], [9] present extensive reviews on it. Providing proper models and clear stability criteria is essential for tasks, such as: 1) setting proper controller gains [10], [11]; 2) establishing technical constraints for dispatching [12]; 3) deciding about load shedding or generation curtailment [13].
Specifically, guaranteeing frequency stability of low-inertia systems can be a major challenge and requires proper sizing of dispatchable power sources and energy storage. Nevertheless, an inspection of recent literature reviews on ac MG planning [14], [15] shows that conditions for frequency stability are largely overlooked in the problem formulation. Specifically, most optimization algorithms assume that matching power demand with production is the only required dynamic constraint, without imposing minimum requirements on system damping or frequency containment reserves (FCR). While such results may be optimal from the affordability point of view, they could be questionable from the security aspect and may be operationally unfeasible.
On the other hand, an inspection of the literature reveals that methods for designing low-inertia power systems in the presence of stochastic loads and non-dispatchable generation This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ considering frequency stability constraints have been proposed lately. Some time ago, [16] introduced a method to optimally size and operate an energy storage system that provides FCR in an isolated power system (IPS) with wind, hydraulic, gas turbines and stochastic loads. Results were validated using a dynamic model in Matlab Simulink. In addition, [17] demonstrated how a fast-acting energy storage can increase the frequency stability of the Guadaloupe archipelago in the Caribbean and reduce the actuation of automatic load shedding using power-hardware-inthe-loop simulations. In [18], the authors also used time-domain simulations (TDSs) to study how the size and controller settings of an energy storage based on ultracapacitors impacted the frequency stability and the performance of the automatic load shedding scheme in the IPS of La Palma, Spain. Later, the same authors extended this work with field-tests results in [19]. Not least, [20] solved the unit commitment problem including FCR constraints for an IPS using mixed-integer linear programming. For that, it determined the minimum system damping empirically relying on the load-frequency sensitivity index. Moreover, frequency stability constraints have been proposed in planning and operation studies, such as optimal power flow problems [21]- [24].
Note that, in the different approaches discussed above, there is a heavy reliance on TDSs and/or model linearization to establish key boundaries for frequency stability. TDSs provide the most accurate results but are hard to integrate in optimization problems such as optimum power flow, while linearizations around an operation point do not give realistic values for frequency deviations during large active power disturbances (APDs). Taking into account this background, the main contribution of this paper is the derivation of sufficient conditions for frequency stability of an ac power system where sources and loads are subject to non-vanishing, bounded perturbations. Using a nonlinear average-model of the ac grid, Lyapunov's direct method and perturbation analysis, an algebraic equation is derived and provides a relationship among four boundaries: the maximum APD, the minimum system damping, the steady-state and the transient frequency limits. The main goal is providing a framework to implement frequency stability constraints in optimization algorithms without recurring to time-consuming simulations or unnecessary linearizations of the power balance equations.
The text is organized as follows: section II introduces a nonlinear, average model of the ac grid that is appropriate for frequency stability analysis during large APDs; section III presents the main proposition of this work followed by its mathematical proof; section IV validates the proposition using time-domain simulations of the IEEE 39-bus test system [26] and compares its results to the linear Swing Equation model; section V discusses applications of the proposition, its main limitations and directions for future research; finally, section VI presents the concluding remarks.

II. NONLINEAR AVERAGE MODEL OF AN AC POWER SYSTEM
When omitting network topology and voltage dynamics, the average behavior of an ac grid angular speed can be modeled using first principles as an equivalent rotating mass using Newton's second law of motion [5], [26]: where J is the equivalent moment of inertia [kgm −2 ] of the system and B(t) is a time-varying function representing the equivalent damping coefficient [Nmsrad −1 ] of the system; T S (t, ω) and T L (t, ω) are the equivalent torques from sources and loads [Nm], functions of the time t [s] and the center of inertia (COI) angular speed ω [rads −1 ]; ω s is a constant representing the synchronous angular speed [rads −1 ]. Proposition 1: The nonlinear system in eq. (1) can be expressed as:˙ where x is the state variable of the system, f (t, x) is piecewise continuous in t and represents the dynamics of the system, g(t, x) is a perturbation term that results from uncertainties and disturbances.
Proof: Eq.(1) can be rearranged into the desired form when multiplying it by ω and normalizing it with the system base power (S b ) [W]: where x = ω ω s and x = x − 1 are the normalized COI angular speed and the COI angular speed deviation from the synchronous speed

III. SUFFICIENT CONDITIONS FOR ROBUST FREQUENCY STABILITY
This section introduces the main proposition of this paper, which is an algebraic expression relating four key boundaries in an ac power system: the maximum APD u(t, x) − w(t, x) < P b [pu], the minimum system damping D min [pu], the steady-state r ss and the transient r tr frequency limits [pu]. Once three of these boundaries are set, the fourth can be unequivocally defined. While the involved mathematical proof of this proposition is presented in detail below, there will be no loss of understanding by proceeding directly to eq. (4), which summarizes the relationship between these quantities.

Proposition 2:
Let u(t, x) − w(t, x) < P b be the maximum APD of the system defined in eq. (2) in which M is a positive constant, D(t) is piecewise continuous in the time t and lower-bounded by a positive constant D min . If where r tr , r ss are the maximum allowed values of x during the transient and steady-state periods, respectively.
Then, for all x(t 0 ) < r ss and r ss < r tr , the solution x of eq. (2) satisfies and some finite T that represents the duration of the transient period.
Proof: Let V( x) = M 2 x 2 be an input-to-state Lyapunov function candidate for the system in eq. (2). Note that V( x) is: 1) a piecewise, continuous differentiable function on R; 2) class K ∞ according to definition 4.2 in [27]; 3) upper-and lower-bounded by M 2 x 2 . Moreover, its time derivative is given by: Assuming that x ∈ D and using the Cauchy-Schwarz's inequality, it holds that: As W 3 ( x) is a continuous positive definite function on R and ρ( z ) is class K, all conditions of the input-to-state theorem in Appendix C are fulfilled. Hence, the system defined in eq. (2) is input-to-state stable on D.
Now, let us analyze the unforced system, i.e. the case where u(t, x) − w(t, x) = 0. The inequalities: With the previous results, it is possible to apply lemma 9.2 in [27]. In this context, when assuming that: for all t ≥ 0, all x ∈ D and 0 < r ss < r tr . Hence, if the maximum APD satisfies the upper-bound: Then, for all x(t 0 ) < r ss and for some finite T , the solution x(t) of eq. (2) satisfies: where the minimum convergence rate of x is: In summary, the steps of the proof are the following: 1) Show that the system in eq. (2) is input-to-state stable on D using the input-to-state theorem in Appendix C.
r Use the Cauchy-Schwarz's inequality to show that 1 x+1 ≤ 1 1−r tr when x ∈ D. 2) Using Theorem 4.10 in [27], show that x = 0 is an exponentially stable equilibrium point of the unforced system, i.e. g(t, x) = 0 in eq. (2). 3) Apply Lemma 9.2 in [27] assuming that: is the perturbation term of the system in eq. (2) r The state must be bounded by r ss , i.e. r t < t 0 is the steady-state period before the APD u(t, x) − w(t, x) < P b appears at t = t 0 . The system is in an equilibrium point x(t 0 ) whose norm is upper-bounded by r ss . r t 0 ≤ t < t 0 + T is the transient period after the APD is applied to the system at t = t 0 in which the state norm r t ≥ t 0 + T is the steady-state period after the system settles down on the new equilibrium point x(t 0 + T ), whose norm is upper-bounded by r ss . Remark 2.2: The proof assumes that x always remains in D.
In other words, Proposion 2 is only valid when x lies between −r tr and r tr . Though, note that this is a reasonable assumption. In ac power systems, when x / ∈ D, load shedding or generation curtailment will take place to reduce the power imbalance and bring it back to D.

Remark 2.3: The inequality
1 x+1 ≤ 1 1−r tr and, hence, Proposion 2 are valid when 0 < r tr < 1, i.e. 0 < ω < 2ω s . For most ac power systems, this is a reasonable assumption. Remark 2.4: As seen in eq. (5), γ and hence the maximum rate of change of frequency (RoCoF) is dependent on D min and M . However, the equilibrium point of the system is independent of M , as highlighted in eq. (4). This statement is formally justified with an equilibrium analysis in Appendix A.
Remark 2.5: From a practical perspective, θ = r ss r tr is a safety margin. The closer to 1 its value is, the closer to marginal stability the system will be after an APD of amplitude P b . When ignoring the term (1 − r tr ), eq. (4) reduces to P b = D min r ss which is the expression retrieved from the linear Swing Equation [28], [29].
Remark 2.6: In the problem formulation, time delays in actuators of primary frequency controllers are neglected. While they affect the dynamics of the system and important benchmarks such as the RoCoF, the frequency nadir and zenith, time delays have no influence in the new equilibrium point of the system after a large APD. Appendix B formally justifies this statement.

IV. VALIDATION USING THE IEEE 39-BUS TEST SYSTEM
At first, the advantages of the proposed model over the linear Swing Equation model may not be obvious. Therefore, this section presents study cases of large APDs in the IEEE 39-bus test system [25]. The aim is to compare results of rms TDSs in DIgSILENT PowerFactory 2020 SP2 with algebraic calculations using eq. (4) and the reference equation P b = D min r ss derived from the linear Swing Equation model. Fig. 1 shows the single-line diagram of the IEEE 39-bus test system and Table I and II introduce relevant system parameters used in the calculations. For all simulations and calculations presented herein, S b = 100 MW. Further details required to reproduce the results in this section are available in [30], which contains the model file and its thorough documentation including a description of modeling differences when compared to [25].

A. Load Damping
The damping coefficient D introduced in eq. (2) can be subdivided into: where D L is the load damping and i D i is the sum of the frequency droop of all primary frequency controllers. In real  systems, D L is typically defined either introducing APDs or applying statistical methods such as system identification to historical data [20].
For the study case, the term i D i = 471.75 is directly retrieved from Table II. However, for the initial conditions shown in Table I, the up-regulation capacity of G5 is extremely limited and this unit will saturate for very small net-load increases. Therefore, its damping is ignored, i.e. i D i = 446.25, and G5 is considered a constant active power source at rated capacity during a net-load increase. To maintain the energy balance, the up-regulation capacity of G5 is subtracted from the APD for all calculations. With the assumptions above, D L is obtained by suddenly increasing bus 4 load by 100 MW (20% step) and analyzing the load-frequency sensitivity index. Fig. 2 shows the COI frequency deviation of the TDSs in PowerFactory and a new equilibrium at −0.002017 pu. Substituting its absolute value as r ss in the linear and eq. (4) models, assuming that r ss = r tr (i.e. no safety margin) and P b = 100−2 100 = 0.98 pu: where D l L is the load damping calculated using the linear model and D p L using eq. (4). Note that D p L is almost 2.5% higher than D l L .

B. Generator 10 Outage
Next, let us compare the results of TDSs with the linear Swing Equation and proposed models for an outage of G10, which is a hydro unit providing most of the FCR. For that, let assume that r tr = 0.02 pu, which is the value for maximum instantaneous frequency deviation allowed in the Nordic system, as defined in article 127 and annex III of the European grid code [31]. The latter also imposes a maximum steady-state frequency deviation of 0.01 pu, which will be considered the target for r ss .
Under these assumptions, let us consider that G10's initial injection of 250 MW and damping of 212.5 pu are lost. where r l ss is the steady-state frequency calculated using the linear model and r p ss using eq. (4). Fig. 3 compares the two calculated limits with the TDSs and shows that the linear Swing Equation model gives optimistic results which are 0.82% below the TDS's new equilibrium at −0.009147 pu. On the other hand, eq. (4) model provides pessimistic results which are 0.84% above the new equilibrium. Regardless, the three obtained references are below the adopted target for r ss and the test system can robustly operate for net-load variations of up to 250 MW without G10.

C. Islanding
In this TDS, the test system is islanded from the bulk grid by disconnecting G1. As before, results are compared to the proposed and linear models.
Under the same assumptions from section IV-B, consider that G1's initial injection is lost. Note that no damping is provided by the interconnection. Thus, i D i = 446.25 and P b = 1000−2 100 = Fig. 4. COI frequency after the system is islanded.
9.98 pu. Substituting in the linear and eq. (4) models: 9.98 39.62 + 446. 25 = 0.02054 This indicates that additional measures for frequency stability are required, such as reducing the initial injection of G1 or increasing FCR. The reader is invited to recalculate r p ss using r tr larger than the new equilibrium to confirm that the proposed model will produce pessimistic values as long as r ss < r tr .
Note that a recalculation of D or P b using the linear model does not guarantee frequency stability, as it produces optimistic values. On the contrary, eq. (4) produces pessimistic values whenever r ss < r tr and provides a stability certificate against the nonlinearities introduced in eq. (2).

V. DISCUSSION
At this point, the utility of Proposion 2 may be more apparent. By deriving sufficient conditions for robust frequency stability of ac power systems operating under large APDs, it defines an algebraic relationship between four system boundaries: the maximum APD P b , the minimum equivalent damping coefficient D min , the transient r tr and r ss steady-state frequency limits. This relationship is summarized in eq. (4) and allows calculation of one of these boundaries once the other three are predefined without resorting to time-consuming TDSs.
As shown in section IV, the proposed algebraic equation produces more realistic forecasts of the COI frequency deviation under large APDs than the linear Swing Equation model. The proposed boundary is robust to nonlinearities and includes a safety factor that depends on the relationship between the maximum instantaneous and steady-state frequency deviations allowed by grid codes or industry standards. One may argue that deviations corrected by eq. (4) are minimal in large ac grids, where frequency deviations and nonlinearities are minimal. Though, that might not be the case in low-inertia systems, where frequency excursions could be considerable under large APDs. Therefore, the use of eq. (4) may be preferred when an algebraic equation must evaluate the frequency stability of generalized ac grids. Note the model described in section II and the theory applied in section III can be used to study frequency stability of MGs dominated by CIGs if primary controllers use droop control or equivalent strategies [32].
Naturally, TDSs still provide the most accurate results and include assessment of voltage and rotor angle stability, which are ignored in the proposed model. However, the algebraic expression given in eq. (4) allows easy integration of frequency stability criteria in long-term optimization algorithms, such as sizing of primary sources and energy storage, uncertain unit commitment and economic dispatch. For instance, some of the authors are contributing to the project HES-OFF [33] in which a software tool is being developed to reduce greenhouse gases emissions of offshore oil and gas installations by integrating a wind farm and hydrogen-based energy storage in their power generation systems. In HES-OFF, eq. (4) is integrated into a multi-objective, long-term optimization procedure that considers frequency stability criteria in the design stage of the ac power system. The procedure is described in detail in [34], [35]. The proposed model reduced the total optimization time by one order of magnitude when compared to the algorithm recurring to TDSs.
Once again, it is important to emphasize that while criteria offered by Proposion 2 are sufficient from a frequency stability perspective, they are only required conditions for the overall stability of an ac power system. The model developed in section II ignores network topology and voltage dynamics, therefore rotor angle and voltage stability criteria are not considered. Evidently, analysis of these dynamics requires detailed models of synchronous machines and CIGs, as described in [36]. To make the problem more tractable, it might be interesting to apply the standard singular perturbation model described in [27] and divide the analysis in two time scales: slow (frequency) and fast (rotor angle and voltage) dynamics. This approach is one possible direction for future research.
Another limitation of this work is assuming that r tr is not violated. As commented in remark 2.6, time delays of primary frequency controllers influence the RoCoF, the frequency nadir and zenith. In low-inertia systems, this might require the reduction of controller delays, as discussed in Appendix B. On one hand, the time-delay of FCR is a variable that typically cannot be influenced. On the other hand, this scenario is changing by the introduction of fast frequency reserves (FFR). For instance, in the Nordic synchronous areas, market solutions for the procurement of FFR have been recently introduced [37], [38]. Consequently, how to optimally size FFR and FCR simultaneously to achieve frequency stability at a limit r ss without infringing a limit r tr during transients is another interesting topic for future research, in which the framing presented in this paper might be useful.

VI. CONCLUSION
Using Lyanpunov's direct method extended by perturbation theory, this article proposed sufficient conditions for frequency stability in ac power systems in the presence of active power disturbances originated by non-dispatchable sources and stochastic loads. By doing that, it defined a clear relationship between four system boundaries: the maximum active power disturbance, the minimum equivalent damping coefficient, the transient and steady-state frequency limits. Hence, when three of these boundaries are defined, the fourth is obtained without recurring to time-consuming simulations. The latter feature allows its integration into long-term optimization algorithms such as: sizing of primary sources and energy storage, uncertain unit commitment and economic dispatch. The proposition is validated using time-domain simulations of the IEEE 39-bus test system and compared to the results of the linear Swing Equation model. Moreover, limitations and possible directions for future research are also discussed.

APPENDIX A EQUILIBRIUM ANALYSIS
Let us define the admissible equilibrium set of the system in eq. (2) as Solving in x , u , w : Note that if u = w (i.e. no active power imbalance), eq. (7) has two solutions: x = 0 and x = −1. The latter is not in D, so only the former is valid. In conclusion, the equilibrium around x = 0 is independent of M and affected only by D and u − w , as highlighted in eq. (4) and (7).

APPENDIX B THE EFFECTS OF TIME DELAYS
Let us assume that primary frequency controllers are affected by time delay in their actuators and/or measurement and these can be approximated by a first-order system as suggested in [39]. In this case, the average model of an ac power system can be modeled by the following cascade system: The admissible equilibrium set of the cascade system is In the equilibrium (i.e.ẏ i = 0), the output of primary frequency controllers is y i = −D i x. As consequence, the equilibrium analysis of the complete system is reduced to solving in x , u , w the following equation: which produces the same results from Appendix A when considering D = D L + i D i . In conclusion, the equilibrium around x = 0 is independent of the time delays i and affected only by the total system damping D and the active power imbalance u − w . By induction, it is also possible to realize that time delays modeled by higher-order functions such as those presented in [40] will render the same results. On the other hand, time delays may affect the stability of the cascade system. Let us consider: Note thatV is a quadratic form, henceV < 0 ⇐⇒ Q > 0, which is indeed the condition for exponential stability of the cascade system. Using the argument that Q sym = 1 2 (Q + Q ) and Kron reduction, it is possible to show that the cascade system is stable if and only if max < M D . However, a proper stability analysis must also consider local and inter-area active power oscillations and rotor angle stability, which is assumed as pre-requisite here. This demands a proper distribution of damping and inertia among primary frequency controllers. An extensive discussion about this topic is found in [36].

APPENDIX C INPUT-TO-STATE STABILITY THEOREM
Theorem 1: 4.19 in [27] Let V : [0, ∞) × R n → R be a continuously differentiable function such that where α 1 , α 2 are class K ∞ functions, ρ is a class K function, and W 3 (x) is a continuous positive definite function on R n . Then, the systeṁ is input-to-state stable with γ = α −1 1 • α 2 • ρ, where f : [0, ∞) × R n × R m → R n is piecewise continuous in t and locally Lipschitz in x and u, and the input u(t) is a piecewise continuous, bounded function of t for all t ≥ 0.