An input-to-state stability approach to verify almost global stability of a synchronous-machine-infinite-bus system

1. Introduction

(a) Motivating challenges

Modern industrialized societies heavily rely on the use and steady supply of electric energy. Therefore, it is of paramount importance to guarantee the stable, reliable and efficient operation of our power systems. However, the latter are bulk, complex and highly nonlinear systems that are continuously subjected to a large variety of disturbances and contingencies [1,2]. Hence, ensuring the stability of a power system is a daunting task, which has been at the core of power system operation since its early days in the 1920s [3,4]—see [1, ch. 1] for a review of the research history on power system stability analysis.

Over the past decades, the complexity of this task has been continuously exacerbated due to the penetration of large shares of volatile renewable energy sources. While this is a very desirable development from environmental and societal factors, the reliable and efficient integration of high amounts of renewable energy sources represents a major technical challenge for the operation of power systems [5,6]. In particular, the increasing presence of such units results in power systems more frequently operating closer to their security margins [5,6]. As a consequence, the problem of power system stability can be expected to become even more relevant in the coming years. This makes it important to derive easily and quickly verifiable analytic conditions for stability, which is the topic addressed in the present work.

A further motivation for our work is the fact that—despite the rich and extensive literature on power system stability—even today many basic questions remain open. This applies to the analysis of vast networks, but also to individual network components and is explained by the complexity of their dynamics. Thus, researchers and practitioners alike are forced to invoke several assumptions simplifying the mathematical tasks of stability analysis and control design. Among the most prevalent assumptions are neglecting fast dynamics [7–9], constant voltage amplitudes and small frequency variations [1, ch. 11, 9]. Invoking such assumptions allows one to derive reduced-order synchronous generator (SG) models [1, ch. 11] and employ algebraic line models [7–9], which significantly simplifies the analysis. Unfortunately, most of the employed assumptions are not physically justifiable in generic operation scenarios. A paramount example for this is the common approximation of the motion of the machine rotor, i.e. the swing equation, in terms of mechanical and electrical power instead of their corresponding torques. However, as this approximation is only valid for small frequency variations around the nominal frequency, it is not an adequate representation of the true physical system when the latter is subjected to large disturbances [1,2,10]; see also the discussions in [11–14]. In summary, a main limitation of current power system stability analysis methods is that they are based on reduced-order models only valid within a limited range of operating conditions. A direct consequence of this is that any stability assessment relying on such reduced-order models may lead to erroneous predictions [13,14], which could lead to severe outages and—in the worst case—even to a blackout.

(b) Existing literature and contributions

In this work, we perform a global stability analysis of a classical and very well-known system called the single-machine-infinite-bus (SMIB) model [2,10]. Thereby, the main distinction to existing results [1,2] is that, motivated by the above discussion, we consider a fourth-order nonlinear SG model derived from first principles—instead of the usual simplified SG models. A thorough analysis and understanding of the SMIB scenario is a fundamental step in the development of a systematic mathematical framework for the analysis of more complex power system dynamics under less restrictive and more realistic assumptions as in [11,12] or in the classical literature [1,2,10]. Furthermore, the insights that can be gained from our analysis can provide valuable practical conclusions for the design, control and operation of SGs as well as for the implementation of an emerging strand of control concepts that aim at operating grid-connected inverters such that they mimic the dynamics of conventional SGs [15–17]. The latter objective also has motivated the related works [18,19]. The analysis in [18] proceeds along the classical lines of constructing an integro–differential equation resembling the forced pendulum equation and, subsequently, showing that the SMIB system is almost globally asymptotically stable if and only if an equilibrium of that forced pendulum equation is almost globally asymptotically stable. In [19], the same authors provide slightly simpler conditions for stability resulting from verifying if a real-valued nonlinear map defined on a finite interval is a contraction. But, as stated in [19], these conditions are hard to verify analytically. Furthermore, and perhaps more importantly, the geometric tools employed to establish the results in [18,19], do not seem to be applicable to a multi-machine power system. In [12], a scenario similar to that of the SMIB system is analysed. However, the analysis in [12] is conducted under very stringent assumptions on the specific form of the infinite bus voltage, as well as the steady-state values of the mechanical torque. A nonlinear analysis of the SMIB system is also conducted in [20], yet using a power rather than a torque balance in the swing equation and neglecting the stator dynamics. Then the SMIB model reduces to a nonlinear forced pendulum, as also studied in [21,22].

To achieve our main objective, i.e. the provision of sufficient conditions for almost global stability of the attractive equilibrium set of the SMIB system, we employ the recently developed framework of input-to-state (ISS) stability for periodic systems with multiple invariant sets [22]. Therein, the ISS approach from [23,24] is combined with the cell structure principle developed in [25–29] to derive necessary and sufficient conditions for ISS of periodic systems. This permits to relax the usual sign definiteness requirements on the Lyapunov function and its time-derivative by exploiting the periodicity of the system. This relaxation is essential to establish the main result of the present paper. More precisely, the main contributions of the paper are threefold.

• — Provide a sufficient condition for ISS of the equilibrium set of the SMIB system.

• — By using this result, establish almost global asymptotic stability of the attractive equilibrium set of the SMIB system, i.e. we show that for all initial conditions, except a set of measure zero, the solutions of the SMIB system asymptotically converge to an asymptotically stable equilibrium point.

• — Illustrate the efficiency of the proposed conditions numerically via a benchmark example taken from Anderson & Fouad [2, Table D.2].

Differently from our previous related work [30,31], in which LaSalle’s invariance principle is employed, the present ISS-based approach has the advantage of providing additional robustness guarantees with respect to exogeneous inputs. Both properties—stability and robustness—are highly desirable from a practical point of view to characterize the performance of power systems under persistent disturbances, such as load variations [32, §D]. Yet, they are also deemed to be hard to establish, in particular, in the global setting considered in this paper [32, §D].

The remainder of the paper is structured as follows. The SMIB model is introduced in §2. The stability analysis is conducted in §3. Section 4 presents a benchmark numerical example. The paper is concluded in §5 with a summary and an outlook on future work.

(i) Notions of input-to-state stability of systems with multiple invariant sets

For the notions of ISS of systems with multiple invariant sets, we follow the conventions and definitions in [22–24] for dynamical systems evolving on a Riemannian manifold of dimension n. Hence, the distance of a point from the set is denoted by , where the symbol Δ(x₁,x₂) denotes the Riemannian distance between x₁ and x₂. Furthermore, |x|=|x|_{0} denotes the usual Euclidean norm of a vector . For further details on properties of dynamical systems and, in particular, systems with multiple invariant sets the reader is referred to Angeli & Efimov [24], Guckenheimer & Holmes [33], Nitecki & Shub[34] and references therein.

2. Model of a synchronous generator connected to an infinite bus

The SMIB system is introduced in this section based on [30,31,35]. We employ a generator reference direction, i.e. current flowing out of the SG terminals is counted positively. Following [30,31], we make some standard assumptions on the SG [15,35], which are also used in [12,18,19]: the rotor is round; the machine has one pole pair per phase; there are no damper windings and no saturation effects; there are no Eddy currents; the rotor current i_f is a real constant. The latter can be achieved by choosing the excitation voltage such that i_f is kept constant [12]. Furthermore, we assume that all three-phase signals are balanced [8,35]. With regard to the SG this is equivalent to assuming a ‘perfectly built’ machine connected in star with no neutral line (see also [18,19]).

We denote the three-phase voltage at the infinite bus by

where is the root mean squared (RMS) value of the constant voltage amplitude (line-to-neutral) and with the grid frequency ω^s being a positive real constant. Furthermore, the electrical rotor angle of the SG is denoted by and the electrical frequency by . As illustrated in figure 1, δ is the angle between the axis of coil a of the SG and the d-axis. Recall that the rotor current i_f is constant by assumption and denote the peak mutual inductance by Then the three-phase electromotive force (EMF) induced in the stator is given by [15,35]

• Download figure
• Open in new tab
• Download PowerPoint

We denote the stator resistance by and the stator inductance by L=L_s−M_s, where is the self-inductance and is the mutual inductance. In the SMIB scenario, the inductance and resistance of the transmission line connecting the SG to the infinite bus are often included in the parameters L and R [2, §4.13.1]. Also note that L>0, since in practice L_s>M_s. Then the electrical equations describing the dynamics of the three-phase stator current are given by

The SMIB model is completed with the electromechanical equations describing the rotor dynamics, given by where the electrical torque T_e is computed as [15] Furthermore, is the total moment of inertia of the rotor masses, is the damping coefficient and is the mechanical torque provided by the prime mover, which we assume is constant.

For our stability analysis, it is convenient to perform a coordinate transformation that maps sinusoidal three-phase waveforms with constant amplitude and frequency—such as (2.1)—to constant signals. This can be achieved by representing the system (2.4)–(2.6) in dq-coordinates via the dq-transformation given in appendix Aa. In this case, we choose the transformation angle

and denote the angle difference between the rotor angle δ and the dq-transformation angle φ by In dq-coordinates, the grid voltage (2.1) is given by the following constant vector [8] where the second equality follows from (2.2). Likewise, the EMF in dq-coordinates is given by where we have defined Hence, the electrical torque T_e in (2.6) is given in dq-coordinates by Furthermore, with φ given in (2.7) and T_dq(⋅) given in (.1), we have that [8, eqn. (4.8)] Consequently, by replacing the rotor angle dynamics, i.e. with the relative rotor angle dynamics, i.e. the SMIB model given by (2.4), (2.5) and (2.6) becomes in dq-coordinates where is an external perturbation which is locally essentially bounded and measurable. The signal d(t) is included to capture model uncertainties originating from possible time-varying excitation current i_f, mechanical torque T_m or voltage amplitudes v_d, v_q. The model (2.11) is used for the analysis in this paper.

Remark 2.1.

In the related works [12,18,19], the dq-transformation angle φ is chosen as the rotor angle, i.e. φ=δ. This yields e_d=0 in (2.8). Furthermore, in [18,19] a dq-transformation is chosen in which the d-axis lags the q-axis by π/2. In our notation, this coordinate frame can be obtained by setting the transformation angle to φ=δ+π. Hence, e_q=bω in [12] and e_q=−bω in [18,19]. The analysis reported in the paper can be conducted in any coordinate frame. Yet, the one employed here is independent of the rotor angle of a specific machine and, hence, seems to be more suitable to extend the results to the multi-machine case, which motivates our choice of φ.

3. Almost global stability

This section is dedicated to the analysis of almost global stability of an equilibrium point of the system (2.11). Our analysis is conducted by employing the recently proposed framework of ISS–Leonov functions [22]. To this end, we perform the following three steps. First, we construct a suitable error system in §3a. Then, we show in §3b that the set of equilibria of the error system is decomposable in the sense of [22, Definition 3]. This is a fundamental prerequisite for the main ISS analysis—the third and final step performed in §3c—which in essence consists of verifying the conditions of [22, Theorem 3]. The almost global stability claim, given in §3d, then follows in a straightforward manner from the ISS property.

For the analysis, we employ the following definition.

Definition 3.1.

An equilibrium point (modulo 2π) of the system (2.11) is almost globally asymptotically stable if it is asymptotically stable and for all initial conditions, except those contained in a set of zero Lebesgue measure, the solutions of the system converge to that equilibrium point (modulo 2π).

(a) Error coordinates

It is convenient to introduce the two constants

Note that (2.9) implies that c is non-zero if the rotor current i_f is non-zero, which is satisfied in any practical scenario.

The existence of isolated equilibria is a natural prerequisite for any stability analysis. Therefore, we make the following natural assumption that by [30, Proposition 1] implies the existence of two unique steady-state solutions (modulo 2π) of the system (2.11).

Assumption 3.2.

The parameters of the system (2.11) are such that i_f>0 and that if for all .

With assumption 3.2, we denote one of the two equilibrium points of the system (2.11) by and introduce the error states

In the error coordinates the dynamics (2.11) take the form Moreover, in the equilibrium point with we have from (2.11) that Hence, by using (3.4) the system (3.3) becomes in error coordinates the equilibrium of which is now shifted to the origin. The remainder of this section is devoted to the analysis of this system.

Remark 3.3.

As can be seen from (2.9), the signs of c and, hence, defined in (3.1) depend on the sign of the constant rotor current i_f. For the subsequent analysis, it is important to know the sign of c, as it determines which of the two equilibria of the system (2.11) is stable. Assumption 3.2 is made without loss of generality because the analysis for i_f<0 follows verbatim.

(b) Decomposable invariant set containing all equilibria

With assumption 3.2, the system (3.5) has a finite and disjoint set of equilibria (modulo 2π) in the sense of [22, Definition 3]. Denote this set by

The lemma below shows that—for bounded solutions— contains all α- and ω-limit sets of the system (3.5) (see [22–24] for details). This feature is instrumental to apply the ISS–Leonov framework [22] and establish our main result.

Lemma 3.4.

Consider the system (3.5) with assumption3.2. Suppose that all its solutions are bounded and that

Then the setis decomposable in the sense of [22, Definition 3] and contains all α- and ω-limit sets of the system (3.5).

Proof.

In order to show that the set contains all α-and ω-limit sets of (3.5) and that it is decomposable, note that the electrical dynamics in (3.5) takes the form

where we defined Following [30], consider the function where x is defined in (3.2), and¹

Clearly, W is continuously differentiable as a function (but not as a function , where denotes the unit circle). With (3.5), a lengthy but straightforward calculation yields

where

In order to show that indeed contains all α- and ω-limit sets of (3.5), we consider the nominal version of the system (3.5), i.e. with for all . Since R>0, the Schur complement implies that M>0 if and only if (3.7) holds. Then, we have immediately that (3.7) implies that

Consequently, by invoking LaSalle’s invariance principle [36] we conclude that all bounded solutions of the nominal system (3.5) converge to the set where From (3.8), we have that implies constant. Since also implies constant, is the only invariant set under the restriction (3.13) and it contains all α- and ω-limit sets of (3.5). Consequently, if for all then implies that is decomposable in the sense of [22, Definition 3], completing the proof. ▪

(c) A condition for input-to-state stability

Recall the function defined in (3.10) and let

where ϵ is a positive real parameter. Denote It is shown in [30, Lemma 1] that c>0 defined in (3.1) is an upper bound for i.e. the infimum in (3.15) is indeed always achievable. The lemma below is useful to establish our main ISS result.

Lemma 3.5.

LetThen for alland any

Proof.

By the definition of we have that

Hence, for any To derive the lower bound for h, we note that since completing the proof. ▪

To streamline the presentation of our result the following assumption is needed. A physical interpretation of this key assumption is given in the discussion in remark 3.10 below.

Assumption 3.6.

There exist constants χ>0 and ψ>0 satisfying

In addition, we borrow the following result from [30].

Lemma 3.7 ([30]).

Consider the system (3.5) verifying assumptions3.2and3.6. Then, all solutionsof the system (3.5) are bounded.

We are now in a position to present our ISS result.

Theorem 3.8.

Consider the system (3.5) with assumptions3.2and3.6. The system (3.5) is ISS with respect to the setdefined in (3.6).

Proof.

Note that lemma 3.7 implies that all solutions of the system (3.5) are bounded under assumption 3.6. Moreover, we see that in (3.16) reduces to (3.7) for χ=ψ=0. As for some χ>0, ψ>0 by assumption, the condition of lemma 3.4 is met. Consequently, the set is decomposable and contains all α- and ω-limit sets of the system (3.5).

The remainder of the claim is established by invoking [22, Theorem 3]. To this end, we need to construct a function V, that—along all solutions of the system (3.5)—satisfies the following two properties

and where g≥0 is a scalar, and is a continuous function satisfying for non-negative arguments. Following [22, Definition 7], such function V is called an ISS–Leonov function.²

Recall W defined in (3.9), h defined in (3.14) and consider the following ISS–Leonov function candidate

In the first step of the proof, we show that V satisfies condition (3.17). Note that lemma 3.5 implies that for all

for some real κ>0. Thus, property (3.17) is satisfied.

The second and last step of the proof consists in verifying condition (3.18). Recall M and M_d given in (3.12). Differentiating V given in (3.19) along solutions of the system (3.5) yields

where We need to verify that Since γ>0, the Schur complement implies that if and only if Furthermore, since χ>0, ε>0 and γ>0, applying the Schur complement also to the matrix in (3.21) implies that if and only if and Recall from lemma 3.5 that Hence, the left-hand side of (3.23) attains its maximum with respect to ε for We select this value for ε and proceed by requiring (3.23) to be satisfied with strict inequality. This ensures that there exists a for which It is then straightforward to verify that this implies that the expression on the left-hand side of (3.23) attains its maximum with respect to ϵ for For this choice of ϵ (3.23) becomes Hence, by defining the auxiliary variable ψ=((Lω^s)²+R²)γ⁻¹>0, (3.23) being satisfied with strict inequality is equivalent to positive definiteness of the matrix in (3.16). Clearly, also implies that (3.22) is satisfied. Recall that by assumption 3.6. Hence, and which is property (3.18). Consequently, with assumption 3.6, V satisfies conditions (3.17) and (3.18), which implies that V is an ISS–Leonov function for the system (3.5). Hence, by invoking [22, Theorem 3] (or, equivalently, [22, Corollary 1]) it follows that the system (3.5) is ISS with respect to the set This completes the proof. ▪

Remark 3.9.

The ISS–Leonov function V defined in (3.19) is not periodic in . Thus, it does not qualify as a usual ISS Lyapunov function [24]. Yet, this non-periodicity of V in is essential for being able to treat the constant terms in the dynamics (3.5) and to establish the ISS claim in proposition 3.8. Furthermore, in order to ensure continuity of the ISS–Leonov function in (3.19) it is necessary to consider the angle evolving in rather than, as usual, on the torus . The same applies to the function W defined in (3.9). In is not bounded a priori. While it is shown in [22, Theorem 3] that existence of an ISS–Leonov function with ensures boundedness of the distance with this does not imply boundedness of for Since boundedness of solutions is a fundamental prerequisite to invoke lemma 3.4 (that ensures decomposability of ), we need the auxiliary lemma 3.7 to establish boundedness of solutions of the system (3.5) evolving in

Remark 3.10.

By using the standard definitions for the reactance X, the conductance G and the susceptance B [2,10], i.e.

condition (3.7)—which is implied by assumption 3.6—can be rewritten as This shows that a high damping factor D and a high R/X ratio, i.e. a high electrical dissipation, are beneficial to ensure ISS of the system (3.5) with respect to the equilibrium set . On the contrary, a high value of |b|, i.e. a high excitation and consequently large EMF amplitude, deteriorate the likelihood of verifying the ISS property for the system (3.5). Both observations are sensible from a physical perspective and consistent with practical experience.

(d) Main stability result

We give now conditions for almost global asymptotic stability of one of the two equilibria (modulo 2π) of the system (2.11). This result follows as a corollary of theorem 3.8 after the relaxation of assumption 3.6—setting ψ=0—as indicated below.

Corollary 3.11.

Consider the system (3.5) with assumption3.2. Suppose thatfor alland that assumption3.6is satisfied forψ=0. Then, the equilibrium pointsatisfying (modulo 2π) with ϕdefined in (3.11) is almost globally asymptotically stable, i.e. for all initial conditions, except a set of measure zero, the solutions of the system (2.11) asymptotically converge to that equilibrium point.

Proof.

For for all we have that γ=0 in in (3.20). Hence, is equivalent to see (3.20). By following the proof of theorem 3.8, with we see that is implied by for ψ=0. As the latter holds by assumption, it follows that, see (3.25),

This shows that for all initial conditions, the solutions of the system (3.5) asymptotically converge to the set . Recall from the proof of theorem 3.8 that, under the standing assumptions, Furthermore, by the definition of in (3.10), we have that Hence, the function W defined in (3.9) has a local minimum at the origin if while it has a local maximum at Therefore, if the origin is locally asymptotically stable, while it is unstable if

To establish almost global asymptotic stability of the equilibrium point satisfying we follow [31] and invoke [37, Proposition 11]. To this end, consider the Jacobian of the system (3.5) evaluated at an equilibrium point, i.e.

By using (3.4) and (3.11), we obtain Hence, if |θ^s−ϕ|>π/2, i.e. if the considered equilibrium point is unstable. Furthermore, is a real-valued matrix of dimension 4. Consequently, if then has at least one positive real eigenvalue and from [37, Proposition 11] it follows that the region of attraction of the unstable equilibrium has zero Lebesgue measure. Hence, for all initial conditions, except a set of measure zero, the solutions of the system (3.5) asymptotically converge to the equilibrium point satisfying (modulo 2π). This completes the proof. ▪

4. Numerical example

The analysis is illustrated on a numerical benchmark example with data taken from [2, Table D.2, H16]. All system parameters are given in per unit (pu) including normalized time (for ω_Base=2π60 rad s⁻¹, table 1). Furthermore, based on [2, Example 5.1], we assume a line resistance of R_ℓ=0.02 pu and a line inductance of L_ℓ=0.4 pu. The nominal damping coefficient is given by D_nom=3⋅2 pu [2, Table D.2].³

Condition (3.16) is a linear matrix inequality with decision variables χ>0 and ψ>0. Hence, its feasibility can be verified efficiently using standard software tools, such as Yalmip [38]. For the subsequent experiments we take D as additional decision variable and seek to determine the minimal damping coefficient required for our condition (3.16) to be satisfied. Thereby, following our previous work [30,31], we consider a large variety of operating scenarios. More precisely, we modify the machine loading by varying T_m, as well as the excitation by varying i_f. All other parameters are kept constant. Furthermore, we define the error coordinates x in (3.2) with respect to the stable equilibrium point, i.e. satisfying .

The resulting values for are shown in figure 2a. Since the machine torque T_e is normalized with respect to the single-phase power base, the rated torque corresponds to T^s_e=3.0 pu. Also, usually it is only possible to load large SGs at their rated power if the amplitude of the EMF e_abc is higher than the rated voltage [39, ch. 5]. This explains the chosen ranges of values for E^s and T^s_e in figure 2a. By theorem 3.8 the SMIB system (3.5) is ISS with respect to the corresponding equilibrium set (if it exists) for any damping coefficient . Furthermore, by corollary 3.11 we conclude that—for for all t≥0—the equilibrium point with (modulo 2π) with ϕ defined in (3.11) is almost globally asymptotically stable. Both properties are further illustrated in figure 2b for D=62.4D_nom, T^s_e=3 pu and |E^s|=1.566 pu. The bold lines correspond to the case while the dashed lines correspond to the case . It can be seen from the simulation results that the trajectories of the perturbed system (i.e. ) remain close to those of the unperturbed system (i.e. ), as predicted by the ISS property.

• Download figure
• Open in new tab
• Download PowerPoint

Moreover, in line with the observations made in [30,31], figure 2a shows that the required damping increases with the machine loading and the excitation magnitude. This is consistent with the well-known power-angle characteristic that can be employed for reduced-order SG models and states that the likelihood of an SG to undergo instability after a change in load increases with the SG being operated closer to its generation limit [1,35].

5. Conclusion

A global stability analysis of the equilibria of a realistic SMIB model has been presented. The main result is established by using the powerful property of ISS, which is verified using an ISS–Leonov function for periodic systems with multiple invariant sets—an approach first advocated in [22]. Compared with the usual Lyapunov functions, ISS–Leonov functions have the advantage that the the conditions of sign definiteness imposed are obviated. This key feature turns out to be essential for the analysis of the SMIB system carried out in this paper. The conservativeness of the derived conditions has been assessed via extensive numerical evaluations on a benchmark problem.

Our current research is directed towards the development of transient stability analysis methods for power system models with time-dependent transmission lines and physically more realistic SG models (e.g. considering torque and voltage control loops as well as saturation effects), which is a natural extension of the presented results. Given the ‘scalable’ nature of the analysis tools employed here this seems a feasible—albeit difficult—task. Another objective is to further reduce the conservativeness of the conditions by improving the construction of the employed Lyapunov function.

Data accessibility

This article has no supporting data.

Authors' contributions

The research design, theoretical results and numerics were performed by J.S. and supported by D.E. with R.O. and N.B. supervising the work. All authors contributed to editing the manuscript.

Competing interests

We declare we have no competing interests.

Funding

The project leading to this manuscript has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 734832, the Government of Russian Federation (grant 074-U01) and the Ministry of Education and Science of Russian Federation (project 14.Z50.31.0031).

Footnotes

One contribution of 13 to a theme issue ‘Energy management: flexibility, risk and optimization’.