Uses of the Popov Stability Criterion for Analyzing Global Asymptotic Stability in Power System Dynamic Models

1. Introduction

Stability studies remain a crucial aspect of power systems dynamic analysis. Despite the fact that stability problems in power systems have been researched for a very long time [1,2,3], the increasing penetration of renewable energy sources, and the rapidly evolving technological advancements, give rise to more intricate and complex system dynamics, that necessitates a deeper and wider perspective regarding their stability [4,5,6]. The core of the problem lies in the fact that renewable energy sources, such as wind and solar, are inherently intermittent and often unpredictable due to their dependence on environmental conditions. This variability introduces challenges in maintaining a stable frequency within the power grid, a parameter that is crucial for the reliable operation of electrical power systems. Frequency stability is maintained by balancing the supply and demand of power within the grid. However, the intermittent nature of renewable sources can lead to significant deviations if not properly managed. These deviations can cause a range of issues, from minor inefficiencies to severe system failures, including blackouts or damage to sensitive infrastructure. Energy storage systems (ESS) have emerged as a vital solution to address these challenges, providing the necessary flexibility to absorb or release energy in response to fluctuations in supply and demand [7,8,9]. By quickly adjusting their output, ESS can help smooth out the variability in renewable energy generation, thereby stabilizing the frequency of the grid.

Many works in the recent literature study stability problems in power systems, using different approaches. This vast field of research can be divided roughly to three main approaches: Numerical analysis [10,11,12], linearization based techniques [13,14,15,16], and direct methods - which utilize at their core various Lyapunov energy functions [17,18,19]. Considering the direct methods, one may examine for instance, work [20], that presents a theoretical foundation of direct methods for both network-reduction and network-preserving power system models. The authors suggest a systematic procedure for constructing energy functions for both network-reduction and network-preserving power system models, utilizing the “boundary of stability region based controlling unstable equilibrium point” (or BCU) method for commuting the unstable equilibrium point, together with an algorithm for the numerical solution of the direct method. Finally, a practical demonstration of the proposed method for online transient stability assessments on two power systems is outlined. Moreover, paper [21] introduces a novel framework to construct the Region of Attraction (ROA) of a power system centered around a stable equilibrium by using stable state trajectories of the system dynamics. The proposed method leverages a Lyapunov function along with the Gaussian process approach. In addition, a sampling algorithm is designed to reconcile the tradeoffs between the exploitation for enlarging the ROA and the exploration for enhancing the confidence level of the sample region. The writers conduct various simulations and experimental validations to substantiate the assessment approach for the ROA of an IEEE test system with real data. It is demonstrated that the proposed approach can significantly enlarge the estimated ROA compared to previous works.

As described above, a large number of works utilize linearization techniques to study the stability of power systems. Due to the vast amount of research in this field, we naturally review only several representative papers. For example, study [22] investigates the application of a multi-variable nonlinear controller for generators. The proposed controller is based on linearization, and its main goal is to improve the transient stability and voltage regulation under post-fault conditions. Simulation results show the improvement in swing stability and system damping under large disturbances, as well as robust terminal voltage regulation. To continue this line of thinking, researchers in [23] offer a new framework to study integration methods for power systems impacted by delays. They use matrix-pencil theory for numerical stability analysis and accuracy analysis. The proposed approach covers all implicit Runge–Kutta methods for time-delay systems. The authors illustrate the idea through simulations on the IEEE 14-bus system, and consider three examples of implicit integration methods, such as “Backward Euler”, “trapezoidal”, and “2-s Radau IIA”.

As seen in the literature review above, there exist three central approaches to analyze the stability of power system dynamic models, namely - numeric simulations, linearization techniques, and direct methods that utilize a variety of Lyapunov functions. In this paper, we continue the line of thinking represented by this last approach, and propose a new criterion for finding global stability properties of equilibrium points in nonlinear power system models. Our approach is based on the Popov stability criterion, which provides exactly this kind of information. With respect to previous works, the proposed method has the following properties:

• With respect to linearization techniques, the proposed approach provides information on the global stability of an equilibrium point of a power system model, rather then on its local stability.
• With respect to numerical techniques, the proposed approach provides analytical conditions, which apply generally to the type of system under study.
• With respect to direct methods, the proposed approach provides a systematic procedure for applying the Popov stability criterion in various power system models, which has not been used before in this application, and thus enables a systematic analysis of several small-scale systems, as described in the text below.

We present first the general mathematical background, and then proceed with practical applications of this proposed approach.

2. Materials and Methods

2.1. Mathematical Background - The Popov Stability Criterion

The Popov stability criterion applies to nonlinear feedback systems that include a memory-less nonlinear component. Unlike linearization techniques, this criterion provides conditions for global asymptotic stability of an equilibrium point. Consider the following dynamic system:

$$\begin{array}{cc}\hfill \dot{x}& =Ax-b\varphi \left(\sigma \right)\hfill \\ \hfill \sigma & =c^{\mathrm{T}}x,\hfill \end{array}$$

(1)

where A is a constant matrix, b and c are constant vectors, and $\varphi \left(\sigma \right)$ is a continuous function of $\sigma$ that describes the nonlinearity. It is assumed that

• The eigenvalues of the constant matrix A are in the open left half-plane,
• $\varphi \left(0\right)=0$,
• $0<\varphi \left(\sigma \right)/\sigma \le k<\infty$, for all $\sigma \ne 0$, where k is a positive constant.

The system is described in Figure 1.

Popov’s stability criterion states that an equilibrium point of the dynamic system (1) is globally asymptotically stable if there exist real numbers $q,\eta$ such that the inequality

$$\mathrm{Re}\left\{\right(1+j\omega q\left)G\right(j\omega \left)\right\}+1/k\ge \eta >0,$$

(2)

holds for every real $\omega \ge 0$. Defining

$$\begin{array}{cc}\hfill G\left(s\right)& =c^{\mathsf{T}}{(sI-A)}^{-1}b,\hfill \\ \hfill X_p\left(\omega \right)& =\mathrm{Re}\left\{G\right(j\omega \left)\right\},\hfill \\ \hfill Y_p\left(\omega \right)& =\omega \mathrm{Im}\left\{G\right(j\omega \left)\right\},\hfill \end{array}$$

(3)

one can plot the “Popov plot”, shown in Figure 2. In this figure, if the plot is indeed in the right-half plane defined by the straight line with slope $1/q$ which intersects with x-axis at $-1/k$, then the conditions of the Popov criterion hold, and the equilibrium point of (1) is globally asymptotically stable, for every $0<\varphi \left(\sigma \right)/\sigma \le k$. We proceed by describing the implementation of this criterion in power system dynamic models.

2.2. Dynamic Stability of a Synchronous Generator, Connected to a Non-Linear Frequency Dependent Load

Consider a system consisting of a synchronous generator connected to non-linear load. We adopt an idealized model, where all the power generated is transferred without loss directly to the load. The system is described by the following differential equation:

$$\mathcal{H}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\omega +D·(\omega -{\omega }_s)=P_{\mathrm{ref}}-P_{\mathrm{load}}\left(\omega \right),$$

(4)

where $\mathcal{H}>0$ represents the moment of inertia of the rotor, $\omega$ is the electrical frequency of the rotor, $D>0$ stands for the damping coefficient, and without loss of generality, $P_{\mathrm{ref}}=P_{\mathrm{load}}\left({\omega }_s\right)$ represents the reference power, where ${\omega }_s$ is the nominal angular velocity. The power consumed by the load is a general non-linear function of the electrical angular frequency $P_{\mathrm{laod}}=P_{\mathrm{laod}}\left(\omega \right)$. For the system (4) there is at least one equilibrium point located at $\omega ={\omega }_s$. It is trivial to verify that linearizing (4) in a small neighborhood of the equilibrium point $\omega ={\omega }_s$ leads to

$$\mathcal{H}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\widehat{\omega }=-\left(D+{\left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{load}}}{\mathrm{d}\omega }}\right|}_{\omega ={\omega }_s}\right)·\widehat{\omega },$$

(5)

where $\widehat{\omega }=\omega -{\omega }_s$, so the equilibrium point $\omega ={\omega }_s$ is locally asymptotically stable if

$$-D<{\left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{load}}}{\mathrm{d}\omega }}\right|}_{\omega ={\omega }_s}.$$

(6)

We extend here this well-known result by providing a sufficient condition for global asymptotic stability of the equilibrium $\omega ={\omega }_s$, that is, we provide a condition under which any trajectory $\omega \left(t\right)$ of (4) converges to ${\omega }_s$ as $t\to \infty$, from any initial point. Our main claim consists of two parts, and it is the following one:

Theorem 1.

Consider the dynamic system (4) with a load profile $P_{\mathit{load}}\left(\omega \right)$, which is generally a nonlinear but continuous function of ω. If there exists a constant k such that

$$-D<{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}}\le k$$

(7)

for all $\omega \ne {\omega }_s$, then the equilibrium point $\omega ={\omega }_s$ is the only equilibrium point of (4), and it isglobally asymptotically stable .

In addition, if $P_{\mathrm{load}}\left(\omega \right)$ is differentiable at ${\omega }_s$, and $-D>{\left.\left(\mathrm{d}{P}_{\mathit{load}}/\mathrm{d}\omega \right)\right|}_{\omega ={\omega }_s}$, or equivalently

$$-D>\underset{\omega \to {\omega }_s}{lim}{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}},$$

(8)

then the equilibrium $\omega ={\omega }_s$ isunstable .

Proof. 

We start by proving the second part of the theorem, the instability condition given in (8). Examine the linearized system described in (5). It is clear that this system has one real pole with value $-D-{\left.{\displaystyle \frac{d}{d\omega }}{P}_{\mathrm{load}}\left(\omega \right)\right|}_{\omega ={\omega }_s}$. A strictly positive pole means that the equilibrium point ${\omega }_s$ is unstable, thus, the system is unstable if

$$-D-{\left.{\displaystyle \frac{d}{d\omega }}{P}_{\mathrm{load}}\left(\omega \right)\right|}_{\omega ={\omega }_s}>0,$$

(9)

which is the same as

$${\left.{\displaystyle \frac{d}{d\omega }}{P}_{\mathrm{load}}\left(\omega \right)\right|}_{\omega ={\omega }_s}<-D,$$

(10)

and by definition of the derivative, we get

$${\left.{\displaystyle \frac{d}{d\omega }}{P}_{\mathrm{load}}\left(\omega \right)\right|}_{\omega ={\omega }_s}=\underset{\omega \to 0}{lim}{\displaystyle \frac{P_{\mathrm{load}}({\omega }_s+\omega )-P_{\mathrm{load}}\left({\omega }_s\right)}{\omega }}.$$

(11)

Now, the following equivalent expression may be used

$${\left.{\displaystyle \frac{d}{d\omega }}{P}_{\mathrm{load}}\left(\omega \right)\right|}_{\omega ={\omega }_s}=\underset{\omega \to {\omega }_s}{lim}{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{load}}\left({\omega }_s\right)}{\omega -{\omega }_s}}.$$

(12)

Substituting (12) in (10) results in

$$\underset{\omega \to {\omega }_s}{lim}{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}}<-D,$$

(13)

which completes the proof for the instability condition in (8).

Now, we proceed to prove the first claim in the main theorem, the stability condition given in (7), based on Popov’s criterion. Let us define $\varphi \left(\omega \right)=P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}$ and $\omega \leftarrow \omega -{\omega }_s$. Substitute these values in (4) to obtain

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\omega =-{\displaystyle \frac{D}{\mathcal{H}}}\omega -{\displaystyle \frac{1}{\mathcal{H}}}\varphi \left(\omega \right).$$

(14)

Assuming condition (7) is true, we now show that (2) holds. For the choice of $x=\omega ,A=-D/\mathcal{H},b=1/\mathcal{H}$ and $\sigma =1·\omega$, the transfer function of the system is

$$G\left(s\right)={\displaystyle \frac{\frac{1}{\mathcal{H}}}{s+\frac{D}{\mathcal{H}}}}.$$

(15)

Following, we shall now find $q\in \mathbb{R}$ and a small and positive constant $\eta$ that satisfy (2). Substitution of (15) in the right expression of (2) results in

$$\mathrm{Re}\left\{(1+j\omega q)G\left(j\omega \right)\right\}+{\displaystyle \frac{1}{k}}={\displaystyle \frac{\frac{1}{H}}{{\left(\frac{D}{H}\right)}^2+{\omega }^2}}\left({\displaystyle \frac{D}{H}}+{\omega }^{2}q\right)+{\displaystyle \frac{1}{k}}.$$

(16)

For the choice of $\omega =0$ we can directly choose the lower bound to be $\eta ={\displaystyle \frac{1}{D}}+{\displaystyle \frac{1}{k}}>0$. Now, by analysing the “Popov plot” we find the values of $q\in \mathbb{R}$ for which (2) holds. To plot the “Popov plot”, we rely on

$$\begin{array}{c}\hfill X_p(s=j\omega )=\mathrm{Re}\left\{G(s=j\omega )\right\}={\displaystyle \frac{D}{{\left(\frac{D}{H}\right)}^2+{\omega }^2}},\\ \hfill Y_p(s=j\omega )=\mathrm{Im}\left\{G(s=j\omega )\right\}={\displaystyle \frac{-\frac{1}{H}}{{\left(\frac{D}{H}\right)}^2+{\omega }^2}}{\omega }^2.\end{array}$$

(17)

Notice that by taking $\omega \to 0$ the resulting point on the “Popov plot” is $({\displaystyle \frac{H^2}{D}},0)$, and for $\omega \to \infty$ the point on the “Popov plot” is $(0,-{\displaystyle \frac{1}{\mathcal{H}}})$. Between these points the transfer function is continuous and is shaped like a quarter of a circle. By the geometric interpretation of (2), it is sufficient to find a straight line, with a slope of $1/q$ that intersects the x-axis at $-1/k$, and is located strictly to the left of the plot of the “Popov plot” of the transfer function. In this case, it is clear that for every choice of a positive slope, the “Popov plot” of the transfer function will be strictly to the right of this straight line, in particular for the value of $q=-1$. Consequently, if condition (7) holds, then there exist $q\in \mathbb{R}$ and a small positive constant $\eta$ such that (2) holds, meaning that the equilibrium is globally asymptotically stable, as desired. □

The following corollary can be used to analyze the stability of a typical small power system, in which a synchronous machine is feeding a general frequency dependant impedance $Z\left(j\omega \right)$.

Corollary 1.

If the load power is given by

$$P_{\mathrm{load}}\left(\omega \right)={\displaystyle \frac{{3\left|V\right|}^2\mathrm{Re}\left\{Z\left(j\omega \right)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}},$$

(18)

where $\left|V\right|$ is constant, $\mathrm{Re}\left\{Z\right(j\omega \left)\right\}$ is bounded from above, $\left|Z\right(j\omega \left)\right|$ is bounded from below by a positive constant, $\mathrm{Re}\left\{Z\left(j\omega \right)\right\},\mathrm{Im}\left\{Z\left(j\omega \right)\right\}\in {\mathbb{C}}^1$, and

$$D>{\displaystyle \frac{{3\left|V\right|}^2}{\omega -{\omega }_s}}\left({\displaystyle \frac{\mathrm{Re}\left\{Z\left(j{\omega }_s\right)\right\}}{|Z\left(j{\omega }_s\right){|}^2}}-{\displaystyle \frac{\mathrm{Re}\left\{Z\right(j\omega \left)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}}\right),$$

(19)

for all $\omega \ne {\omega }_s$, then the equilibrium point $\omega ={\omega }_s$ is the only equilibrium point of (4), and it is globally asymptotically stable.

Proof. 

To prove this statement, we show that it satisfies the conditions of Theorem 1. Substituting (18) and $P_{\mathrm{ref}}=P_{\mathrm{load}}\left({\omega }_s\right)$ in Corollary 1 results in

$$D>-{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}},$$

(20)

which yields directly the lower bound:

$$-D<{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}}.$$

(21)

Following, using the fact that $\left|V\right|$ is constant, $\mathrm{Re}\left\{Z\left(j\omega \right)\right\}\le k_1$, and $0<k_2\le \left|Z\left(j\omega \right)\right|$, we proceed to determine the upper bound. Examining the expression

$${\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}},$$

(22)

and substituting (18), results in

$${\displaystyle \frac{{3\left|V\right|}^2\left(\frac{\mathrm{Re}\left\{Z\left(j{\omega }_s\right)\right\}}{|Z\left(j{\omega }_s\right){|}^2}-\frac{\mathrm{Re}\left\{Z\right(j\omega \left)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}\right)}{\omega -{\omega }_s}}.$$

(23)

For $\omega \to 0$ and $\omega \to \infty$ the expression (23) is well defined, and the upper and lower bounds can be calculated directly. Let us denote them by $c_l,c_u$ respectively. In the regions $[0,{\omega }_s)$ and $({\omega }_s,\infty )$ the expression (23) is continuous, thus there are no additional extremum points that can be found there. Finally, consider the limit

$$\underset{\omega \to {\omega }_s}{lim}\left({\displaystyle \frac{\mathrm{Re}\left\{Z\left(j{\omega }_s\right)\right\}}{|Z\left(j{\omega }_s\right){|}^2}}-{\displaystyle \frac{\mathrm{Re}\left\{Z\right(j\omega \left)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}}\right){\displaystyle \frac{1}{\omega -{\omega }_s}}.$$

(24)

We now show that this limit exists, and is finite. That is, our goal is to prove that

$$\underset{\omega \to {\omega }_s}{lim}\left({\displaystyle \frac{Re\left\{Z\left(j{\omega }_s\right)\right\}}{|Z\left(j{\omega }_s\right){|}^2}}-{\displaystyle \frac{\mathrm{Re}\left\{Z\right(j\omega \left)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}}\right){\displaystyle \frac{1}{\omega -{\omega }_s}}=C.$$

(25)

We rely on the following assumptions presented in Theorem 1:

• There exists $k_1$ such that $\mathrm{Re}\left\{Z\left(j\omega \right)\right\}\le k_1$,
• There exists $k_2$ such that $\left|Z\left(j\omega \right)\right|\ge k_2$,
• $\mathrm{Re}\left\{Z\left(j\omega \right)\right\},\mathrm{Im}\left\{Z\left(j\omega \right)\right\}\in {\mathbb{C}}^1$.

Examining the limit (24) again, application of L’Hopital’s rule results in

$$\underset{\omega \to {\omega }_s}{lim}{\displaystyle \frac{\frac{\mathrm{d}}{\mathrm{d}\omega }\left(\frac{\mathrm{Re}\left\{Z\right(j\omega \left)\right\}}{{\left|Z\left(j\omega \right)\right|}^2}\right)}{\frac{\mathrm{d}}{\mathrm{d}\omega }(\omega -{\omega }_s)}}=\underset{\omega \to {\omega }_s}{lim}{\displaystyle \frac{\left(\frac{\mathrm{d}}{\mathrm{d}\omega }\mathrm{Re}\left\{Z\left(j\omega \right)\right\}\right){\left|Z\left(j\omega \right)\right|}^2-\mathrm{Re}\left\{Z\left(j\omega \right)\right\}\left(\frac{\mathrm{d}}{\mathrm{d}\omega }{\left|Z\left(j\omega \right)\right|}^2\right)}{{\left|Z\left(j\omega \right)\right|}^4}}.$$

(26)

Since $\mathrm{Re}\left\{Z\left(j\omega \right)\right\},\mathrm{Im}\left\{Z\left(j\omega \right)\right\}\in {\mathbb{C}}^1$, we have that $\left|Z\left(j\omega \right)\right|={\left(\mathrm{Re}\left\{Z\right(j\omega \left)\right\}\right)}^2+{\left(\mathrm{Im}\left\{Z\right(j\omega \left)\right\}\right)}^2\in {\mathbb{C}}^1$. Hence, the following limits are well defined:

• ${lim}_{\omega \to {\omega }_s}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\omega }}\mathrm{Re}\left\{Z\left(j\omega \right)\right\}=c_1$,
• ${lim}_{\omega \to {\omega }_s}{\left|Z\left(j\omega \right)\right|}^2=c_2$,
• ${lim}_{\omega \to {\omega }_s}\mathrm{Re}\left\{Z\left(j\omega \right)\right\}=c_3$,
• ${lim}_{\omega \to {\omega }_s}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\omega }}{\left|Z\left(j\omega \right)\right|}^2=c_4$.

Therefore, the original limit we are interested in is equal to

$${\displaystyle \frac{c_1{c}_2-c_3{c}_4}{c_2^4}}\equiv C.$$

(27)

Note that, depending on the system’s parameters, it is possible to get $C>0$ or $C<0$. In each case, the value of C is a candidate for upper or lower limit respectively. w.l.o.g. assuming $C>0$ we deduce it is a candidate for the upper bound. Consequently, (22) is bounded by

$$-D<{\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}}\le k,$$

(28)

where $k=max\{c_u,C\}$ and $D=c_l$. Thereby, the conditions of Theorem 1 holds in this unique case, meaning that $\omega ={\omega }_s$ is indeed the only equilibrium point of (4), and it is globally asymptotically stable. □

2.3. Synchronous Machine Driving a Resistive-Inductive Load

In this case study, we analyze a system in which the nonlinear load has a resistive-inductive behavior. The main claim used for the analysis is Corollary 1. Our aim is to find the minimal value of D which bounds from above the expression (19), which is derived from the active power of the system. For this type of model, the active power is given by

$$P_{\mathrm{load}}=3{\left|V\right|}^2·{\displaystyle \frac{R}{R^2+{\left(\omega L\right)}^2}}.$$

(29)

Examine the expression

$${\displaystyle \frac{P_{\mathrm{load}}-P_{\mathrm{ref}}}{\omega -{\omega }_s}}.$$

(30)

Substitution of (29) in (30) leads to

$${\displaystyle \frac{{3\left|V\right|}^2}{\omega -{\omega }_s}}\left({\displaystyle \frac{R}{R^2+{\left({\omega }_{s}L\right)}^2}}-{\displaystyle \frac{R}{R^2+{\left(\omega L\right)}^2}}\right).$$

(31)

Re-writing this equation results in

$${\displaystyle \frac{1}{\omega -{\omega }_s}}·{\displaystyle \frac{{\omega }^2-{\omega }_s^2}{(R^2+{\omega }_s^2{L}^2)·(R^2+{\omega }^2{L}^2)}}·3{\left|V\right|}^{2}R{L}^2,$$

(32)

or alternatively

$${\displaystyle \frac{\omega +{\omega }_s}{R^2+{\omega }^2{L}^2}}·{\displaystyle \frac{{3\left|V\right|}^{2}R{L}^2}{R^2+{\omega }_s^2{L}^2}}.$$

(33)

Next, we aim to find the minimal value that bounds (33) from above. For this purpose, notice that the expression given in (33) consists of continuous functions of $\omega$. Hence, the derivative may be calculated directly and is given by

$${\displaystyle \frac{{3\left|V\right|}^{2}R}{R^2+{\omega }_s^2{L}^2}}·\left[{\displaystyle \frac{-{\omega }^2-2{\omega }_s\omega +{\left(\frac{R}{L}\right)}^2}{{\left({\omega }^2+{\left(\frac{R}{L}\right)}^2\right)}^2}}\right],$$

(34)

which critical points are

$$r_{1,2}=-{\omega }_s\pm \sqrt{{\omega }_s^2+{\left({\displaystyle \frac{R}{L}}\right)}^2}.$$

(35)

Since $\omega >0$, the only valid root is $r_2$, so to analyze (34) we examine its behavior close to the critical point $r_2$. Given that ${\displaystyle \frac{{\left|V\right|}^{2}R}{R^2+{\omega }_s^2{L}^2}}>0$ and ${\left({\omega }^2+{\left({\displaystyle \frac{R}{L}}\right)}^2\right)}^2>0$, it is sufficient to look at the sign of the expression $-{\omega }^2-2{\omega }_s\omega +{{\displaystyle \frac{R}{L}}}^2$. When $\omega \to r_2^{-}$, the derivative (34) is negative, and when $\omega \to r_2^{+}$, the expression (34) is positive, meaning that $r_2$ is a maximum point. Consequently, the upper bound may be defined as

$$D={\displaystyle \frac{\sqrt{{\omega }_s^2+{\left(\frac{R}{L}\right)}^2}}{R^2+{\left(-{\omega }_s+\sqrt{{\omega }_s^2+{\left(\frac{R}{L}\right)}^2}\right)}^2·L^2}}·{\displaystyle \frac{{3\left|V\right|}^{2}R{L}^2}{R^2+{\omega }_s^2{L}^2}}.$$

(36)

This results in a strict upper bound for the expression given in (33). Indeed, the condition given in (19) holds, therefore, ${\omega }_s$ is the only global asymptotic equilibrium point of the given system, as desired, and a simplified expression is

$$D={\displaystyle \frac{\sqrt{{\omega }_s^2+{\left(\frac{R}{L}\right)}^2}}{{\frac{R}{L}}^2+{\left(-{\omega }_s+\sqrt{{\omega }_s^2+{\left(\frac{R}{L}\right)}^2}\right)}^2}}·{\displaystyle \frac{{3\left|V\right|}^{2}R}{R^2+{\omega }_s^2{L}^2}}.$$

(37)

Notice that the active power at the operating point is given by $P^{*}={\displaystyle \frac{{3\left|V\right|}^{2}R}{R^2+{\omega }_s^2{L}^2}}$ and the reactive power at the operating point is defined as $Q^{*}={\displaystyle \frac{{3\left|V\right|}^2{\omega }_{s}L}{R^2+{\omega }_s^2{L}^2}}$. We therefore have ${\displaystyle \frac{R}{L}}={\omega }_S·{\displaystyle \frac{P^{*}}{Q^{*}}}$. Substituting these terms in (37) produces

$$D={\displaystyle \frac{\sqrt{{\omega }_s^2+{\omega }_s^2·{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}{{\omega }_s^2·{\frac{P^{*}}{Q^{*}}}^2+{\left(-{\omega }_s+\sqrt{{\omega }_s^2+{\omega }_s^2·{\left(\frac{P^{*}}{Q^{*}}\right)}^2}\right)}^2}}·P^{*},$$

(38)

which is the same as

$$D={\displaystyle \frac{\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}{{\frac{P^{*}}{Q^{*}}}^2+{\left(-1+\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}\right)}^2}}·{\displaystyle \frac{P^{*}}{{\omega }_s}}.$$

(39)

Re-arranging the terms in the expression above, yields

$$D={\displaystyle \frac{\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}{{\frac{P^{*}}{Q^{*}}}^2+1-2\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}+1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}·{\displaystyle \frac{P^{*}}{{\omega }_s}},$$

(40)

or

$$D={\displaystyle \frac{\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}{{\frac{P^{*}}{Q^{*}}}^2+1-\sqrt{1+{\left(\frac{P^{*}}{Q^{*}}\right)}^2}}}·{\displaystyle \frac{P^{*}}{2{\omega }_s}},$$

(41)

which may be written as

$$D={\displaystyle \frac{1}{\sqrt{{\frac{P^{*}}{Q^{*}}}^2+1}-1}}·{\displaystyle \frac{P^{*}}{2{\omega }_s}},$$

(42)

thereby leading to

$$D={\displaystyle \frac{Q}{\sqrt{P^{*}+{Q^{*}}^2}-Q}}·{\displaystyle \frac{P^{*}}{2{\omega }_s}}.$$

(43)

In conclusion, using the fact that ${\left|S\right|}^2=P^2+Q^2$, the bound on D may be presented in the following way

$$D={\displaystyle \frac{Q}{\left|S\right|-Q}}·{\displaystyle \frac{P^{*}}{2{\omega }_s}}.$$

(44)

2.4. Synchronous Machine Driving a Lossless Synchronous Motor, with a Quasi-Linear Mechanical Load

Examine the following non-linear load function:

$$P_{\mathrm{load}}\left(\omega \right)=\tau \left(\omega \right)·\omega ,$$

(45)

where $\tau \left(\omega \right)$ is given by

$$\left\{\begin{array}{cc}\alpha ·\omega ,\hfill & \omega \le {\omega }_0\hfill \\ \alpha ·{\omega }_0,\hfill & \omega >{\omega }_0\hfill \end{array}\right.,$$

(46)

where $\alpha ={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s^2}}$.

To prove that $\omega ={\omega }_s$ is indeed the global asymptotic stable point, we aim to prove that the conditions of Theorem 1 hold. We start by analyzing separately the domains $[0,{\omega }_0]$ and $({\omega }_0,\infty )$. Define two constants for each domain $c_l\ge 0,c_u>0$ which denote the upper and lower bound of $\tau \left(\omega \right)$ accordingly. Following, the global upper and lower bounds are given by $k=max\{c_u^1,c_u^2\}$ and $D=min\{c_l^1,c_l^2\}$, where the superscript $i\in \{1,2\}$ denotes the domain. Consider the domain $[0,{\omega }_0]$, the expression that we are interested in is given by

$${\displaystyle \frac{P_{\mathrm{load}}\left(\omega \right)-P_{\mathrm{ref}}}{\omega -{\omega }_s}}.$$

(47)

Alternatively, we have

$${\displaystyle \frac{\frac{P_{\mathrm{ref}}}{{\omega }_s^2}·{\omega }^2-P_{\mathrm{ref}}}{\omega -{\omega }_s}},$$

(48)

which is equal to

$${\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{{\left(\frac{\omega }{{\omega }_s}\right)}^2-1}{\frac{\omega }{{\omega }_s}-1}}.$$

(49)

This can be further simplified as follows:

$${\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·\left({\displaystyle \frac{\omega }{{\omega }_s}}+1\right).$$

(50)

Using the extreme value theorem, since the expression (47) is continuous on a closed and bounded interval $[0,{\omega }_0]$, it must attain a minimum and maximum values there. Moreover, since the function is monotonic, it is assured it will attain them only once on this interval. Hence, (50) is bounded from below by $c_l^1=0$, and since it is monotonically increasing, the maximum value is obtained at ${\omega }_0$ and equals to $c_u^1={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1)$. Next, looking at the domain $({\omega }_0,\infty )$, substituting (45) in the expression (47) equals to

$${\displaystyle \frac{\frac{P_{\mathrm{ref}}}{{\omega }_s^2}·{\omega }_0·\omega -P_{\mathrm{ref}}}{\omega -{\omega }_s}},$$

(51)

which yields

$${\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{\frac{{\omega }_0}{{\omega }_s}·\frac{\omega }{{\omega }_s}-1}{\frac{\omega }{{\omega }_s}-1}}.$$

(52)

Note that ${\displaystyle \frac{{\omega }_0}{{\omega }_s}}>1$ since ${\omega }_0>{\omega }_s$. Hence, the function (52) is continuous and strictly monotonically decreasing on the interval $({\omega }_0,\infty )$, which leads to a lower bound

$$c_l^2={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{{\left(\frac{{\omega }_0}{{\omega }_s}\right)}^2-1}{\frac{{\omega }_0}{{\omega }_s}-1}},$$

(53)

which reduces to

$$c_l^2={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·\left({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1\right).$$

(54)

The upper bound on this interval can be calculated by examining the limit $\omega \to \infty$, which is

$$\underset{\omega \to \infty }{lim}\left({\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{\frac{{\omega }_0}{{\omega }_s}·\frac{\omega }{{\omega }_s}-1}{\frac{\omega }{{\omega }_s}-1}}\right)={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{{\omega }_0}{{\omega }_s}}.$$

(55)

Consequently, we can choose $c_u^2={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{{\omega }_0}{{\omega }_s}}$ as the upper bound on the current interval. Summarizing, we can globally bound (47) from below by

$$\begin{array}{cc}\hfill D& =max\left\{0,{\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·\left({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1\right)\right\}\hfill \\ \hfill & ={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·\left({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1\right),\hfill \end{array}$$

(56)

where the global upper bound is

$$\begin{array}{cc}\hfill k_0& =max\left\{{\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1),{\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·{\displaystyle \frac{{\omega }_0}{{\omega }_s}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{P_{\mathrm{ref}}}{{\omega }_s}}·\left({\displaystyle \frac{{\omega }_0}{{\omega }_s}}+1\right).\hfill \end{array}$$

(57)

This proves that the conditions of Theorem 1 hold, and therefore $\omega ={\omega }_s$ is indeed the global asymptotic equilibrium point, as desired.

3. Numeric Results

To demonstrate this result, we consider here the classical problem of managing an ideal grid connected storage device, and focus on the behavior of the poles given different values of the damping coefficient of the generator. Consider a microgrid system comprising a grid-connected ideal storage device and a photovoltaic (PV) panel, as illustrated in Figure 3.

The storage device is charged by both the grid and the PV, and supplies an aggregated load characterized by its active power consumption. The load’s active power demand, represented by a continuous, positive function $P_a:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$, is defined over a finite time interval $[0,T]$, where T is a known constant. The power generated by the PV panel is modeled as a piecewise continuous, non-negative function $P_{pv}:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$. The net power consumption of the load is given by $P_L\left(t\right)=P_a\left(t\right)-P_{pv}\left(t\right)$. The rate of charge or discharge of the storage device is described by $P_s\left(t\right)=P_g\left(t\right)-P_L\left(t\right)$, where $P_g:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ is the power supplied by the grid, and $P_s:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ is the power flowing into the storage device. Additionally, the quantities $E_g\left(t\right),E_L\left(t\right)$, and $E_s\left(t\right)$ represent the generated energy, the load energy, and the stored energy, respectively. These energy values are derived from their corresponding power functions via the relation $E\left(t\right)={\int }_0^{T}P\left(\tau \right)\mathrm{d}\tau$.

To verify our finding we performed a MATLAB simulation, using Simulink to model the storage device as an R-L circuit, using the idea presented in [24] and the analysis presented in Section 2.3. The solar generation and the consumption profiles are presented in Figure 4.

The critical value of the bound D which defines the stability condition is $D_{th}=95.402$. Performing the simulation for different values of D, shows clearly that for all values of D which are above the threshold value $D_{th}$, the system is stable and the frequency converges to the globally asymptotically stable point ${\omega }_s\approx 314$. On the contrary, for all values of D that are beneath the threshold $D_{th}$ the system diverges. For values of D in which the system is stable, the scheduling policy of the storage device is presented in Figure 5.

Examining Figure 6, it can be observed how for different values of D, the frequency either stabilizes on the unique global stable point ${\omega }_s$, or diverges. Figure 7 presents how the poles move as the values of the damping coefficient D and the resistance R changes, thus showing the transition of the system from a stable to an unstable state. We examine the following values of the resistance $R=[0.1,0.5,1,5,10,100]$ in units of $\Omega$, and the damping coefficient values tested are $D=[0.01,0.03,0.04,0.06,0.5,1,10,D_{th},D_{th}+10]$.

4. Discussion & Conclusions

Stability problems in power systems are typically explored using one of three main approaches: numerical methods, linearization techniques, or direct methods, which utilize Lyapunov energy functions. This paper belongs to the third category, and highlights the usefulness of the Popov stability criterion in the analysis of nonlinear power system models. The main advantage of this criterion is that it provides conditions for global asymptotic stability of an equilibrium point. This stands in contrast to linearization techniques, which only provides conditions for local stability. We show a general method to apply this stability criterion, and examine its uses in several specific applications. More specifically, we examine small power systems such as synchronous machine driving a resistive-inductive load and a synchronous machine driving a lossless synchronous motor, with a quasi-linear mechanical load. Following we rely on a simplified model of a storage device that is powered by a synchronous machine, to analyze the frequency stability for varying system parameters. For instance, we evaluate for which practical parameters of the storage and the synchronous machine, such as the resistance and damping coefficient, the equilibrium is globally asymptotically stable.