Generalized H matrix stability for power systems

Notations

$\mathcal{M}$ — matrix of generators’ inertia constants

$\mathcal{D}$ — matrix of generators’ damping coefficients

$\mathcal{U}$ — matrix mapping generator angles to relative rotor angles

$\mathcal{H}$ — matrix of partial derivatives of generator active powers w.r.t. relative rotor angles

$\mathcal{J}$ — full Jacobian matrix of active powers w.r.t. relative rotor angles

$Y$ — admittance matrix of the transmission network at the nominal frequency

$A$ — system matrix of the linearized dynamic system

$\omega$ — vector of generator electric frequencies relative to the reference frequency

$V_i$ — voltage amplitude at bus $i$

$P_i$ — active power injected at bus $i$

$X_{-}(A)$ — strictly stable subspace of matrix $A$

$X_{+}(A)$ — strictly unstable subspace of matrix $A$

$\sigma(A)$ — the set of eigenvalues of matrix $A$

Introduction

Within the field of electrical engineering, power system stability remains a primary focus of contemporary research, primarily driven by the increasing penetration of distributed energy resources (DERs) and the ongoing structural evolution of modern electrical grids. These stability challenges are exacerbated by the inherent complexity and heterogeneity of the underlying power system models [1, 2].

The stability of power system models has been researched for many years, using many methods and points of view, see for instance [3, 4, 5, 6, 7, 8, 9, 10]. The literature in this field is extensive (the aforementioned citations representing only a subset) and originates from diverse research disciplines over many years. In general, finding sufficient conditions for the stability of a power network is extremely difficult. This is because power networks are nonlinear and high-dimensional systems. Every type of analysis is based on an approximate model of the system. For instance, an individual synchronous generator can be modeled with 0, 1, or 2 damper windings, with constant field current or with varying field current, and so forth. The

REFERENCES

complexity of synchronous machine models ranges from order 9, see for instance [11 Chapter 4] or the survey [2], down to order 2, see for instance [10], [12]. Likewise, prime movers and their governors, as well as loads, can be modeled with varying complexity [9], [10], [12]–[15]. Clearly, the more complex the model of the power network, the more difficult it is to prove stability results for it. While computational methods may be used to simulate the behaviour of any specific system [16], such methods are generally non-scalable, and rarely lead to deep insights.

The approximation of time-varying phasors is commonly used to model dynamic phenomena in large-scale power systems [17]. This approach describes voltages and currents with time-varying phasor amplitudes and phases, allowing the transmission network to be modeled using power-flow equations and simplified phasor-based models for generators and loads. The core assumption is that system transients are significantly slower than the nominal system frequency [18], [19]. While this generally yields good results for slow transients in large systems, its accuracy may diminish with high renewable energy penetration or in systems with fast transients [18], [20]–[24]. Still, time-varying phasor models are very popular for several reasons. Firstly, they use intuitive quantities like amplitude, phase, active power, and reactive power, which power engineers readily understand as they stem directly from phasor notations. Secondly, these models are simpler and more compact numerically than other, more complex models, making them computationally efficient. Thirdly, time-varying phasor models are remarkably flexible, capable of describing a wide variety of sources and loads, with numerous models already documented in literature. Finally, their main assumption of slow transients holds true for most large-scale network explorations.

Within the framework of time-varying phasors, the dynamic stability of power systems has also been studied extensively. A few examples include the use of Kuramoto models [8], [12], [25]–[27], network-reduced power system models (NRPS) [5], [15], DC power flow [3], [28], linearization techniques [11], and Lyapunov functions [16], [29], [30], among others. Alternatively, predictive methods have been proposed to assess small-signal stability, particularly when system configurations are complex or unknown. For instance, [31] introduces an automated dynamic frequency scanning tool designed for the analysis of systems with unknown internal configurations. By injecting multi-sine disturbances into an electromagnetic transient environment and analyzing the resulting waveforms, the tool extracts the small-signal frequency-dependent impedance of the system at a specific operating point. This enables the characterization of the small-signal frequency response and the prediction of stability margins without requiring access to the underlying control structures or internal parameters of the components. These models serve as practical tools for real-world analysis, as demonstrated in [32]. Gaining a better understanding of the fundamental characteristics of the matrices used for stability analysis may be highly relevant for recent research focusing

on active control methods, such as [33], [34] and [35]. The properties of the $H$ matrix may establish a practical direction for such control strategies, where the synchronizing power coefficients are used as a real-time stability objective. Within this framework, virtual inertia and damping parameters might be adaptively tuned to ensure that the $H$ matrix remains positive definite. One prominent method for assessing the small-signal stability of power systems that use phasor-based dynamic models is the $H$ matrix stability criterion, as detailed for instance in [36] and [37]. The classical application of this criterion is specifically designed for a particular power system model, that incorporates time-varying phasors, uses second order models for generators, assumes a lossless transmission network, and considers constant power loads. In this context, the matrix $H$ is computed at a certain equilibrium point $x^*$, and is constructed from the partial derivatives of generator powers with respect to their relative angles. The stability of the equilibrium point is then determined as follows: If $H$ is positive definite, the equilibrium point is asymptotically stable, and if $H$ is full-rank but not positive definite, the equilibrium point is unstable (see [36]).

Recently, data-driven frameworks such as the one presented in [38], are proposed to learn power system dynamics by incorporating physical laws into the network architecture. In comparison to such data-driven approaches, the $H$ matrix stability criterion provides analytical rigor through a necessary and sufficient condition for stability (under the assumptions stated in Section II), where data-driven frameworks serve as a predictive tool, whose accuracy is dependent on training data, and can suffer from loss of previously learned dynamics when transitioning between stable and unstable regimes. The $H$ matrix method also exhibits higher online computational speed, as the assessment is performed directly without the requirement for intensive training. For example, when examining a large-scale system, such as the 2869-bus system, the stability of the equilibrium point is determined in 0.0011 seconds using the $H$ matrix stability criterion. Alternatively, the proposed data-driven method, while more efficient than traditional data-driven approaches, still requires 0.024 seconds to solve a small 4-bus system. Furthermore, the $H$ matrix stability criterion is much more interpretable than any data-driven methods.

While a detailed comparison of the proposed stability criterion to previous results may be found in Section V, the most notable advancement lies in the topology of the transmission network. Previous works, for instance [36], assume that loads are situated on the same buses as generators. The proposed updated criterion eliminates this assumption, allowing loads to be located anywhere in the transmission network. This enhancement enables the $H$ matrix criterion to be applied to any lossless transmission network model, making it more general and broadly applicable than previously available methods. The rather nontrivial fact that under the stated assumptions $H$ is symmetric and self-adjoint is also proven. Furthermore, several other relations between the eigenvalues of $H$ and the eigenvalues of the system’s Jacobian matrix $A$ are established. For instance, it is shown that, under certain conditions, the number of unstable eigenvalues of $A$ equals the number of negative eigenvalues of $H$ when algebraic multiplicities are counted. A number of case studies, both theoretical and numerical, that demonstrate the usefulness of the proposed updated criterion are also explored.

The main result, Theorem 2, is an updated and generalized $H$ matrix stability criterion for phasor-based dynamic power grid models with locally constant power loads. Various extensions to arbitrary power loads are also derived in Section IV. In the model, each generator is represented by a swing equation, as in (2). Such a generator may also be an inverter controlled to emulate a synchronous generator, which is an often employed control strategy for grid forming inverters, see [39, 40, 15].

The contributions of the paper are summarized as follows:

1. An updated stability criterion is developed for the analysis of small-signal stability in power systems with arbitrary load distributions.

2. The symmetry of the $H$ matrix is formally established for lossless transmission networks with locally constant power loads.

3. Quantitative relationships between the eigenvalues of the $H$ matrix and the system Jacobian are derived to characterize stability and instability.

The advantage of the proposed criterion is demonstrated by several case studies, which show that correct stability assessments are maintained in large-scale power systems where alternative methodologies, such as network reduction methods, are found to be overly conservative and do not yield unequivocal results. This limitation is exemplified in the 2869-bus Pegase system, where the network reduction approaches established in [12] and [8] are unable to determine whether the system is stable. From a different perspective, although correct results are achieved by the methodology proposed in [13], through the use of the full Jacobian matrix, the process is considerably slower. Specifically, the $H$ matrix calculation is performed over 340 times faster than the Jacobian-based approach for this system. These efficiency gains are obtained because the $H$ matrix is characterized by a dimension approximately half that of the system matrix $A$, which results in substantial reductions in both computational time and memory requirements.

The paper continues as follows: In Section II, basic definitions and the problem formulation are provided. It is also proven that in the case of locally constant power loads, $H$ is symmetric. Section III contains the main result, Section IV generalizes the result for arbitrary loads, and Section V discusses the relation between the main result of this paper and related earlier results. Sections VI–VIII provide case studies. A short conference version of this material (with some proofs omitted and with different examples) is in [41].

Problem formulation

A general balanced AC power system is considered, whose dynamic behavior is governed by the quasi-stationary model,

as considered for instance in [17, 30, 36]. Note that this paper relies only on the quasi-static model of a transmission network with loads. In this model, the consumed active and reactive powers are functions of the generator voltages and phase differences only, and there are no differential equations in the description of the transmission network. The system is illustrated in Figure 1, and includes a general lossless transmission network, $G$ generators, $L$ loads and hence $N = G + L$ relevant buses. Without loss of generality, it is assumed that the generator buses are labeled by $1,\ldots,G$, and the load buses are labeled by $(G+1),\ldots,N$.

Several assumptions are employed that are commonly used in time-varying phasor models:

• The power system is a balanced three-phase system. Transients in the system are much slower than a 50/60 Hz cycle, and dynamic phenomena are represented by phasors, based on the quasi-static (or quasi-stationary) approximation [42].

• The transmission network is represented by a linear single-phase equivalent circuit, and can be described well enough by a constant and symmetric admittance matrix $Y$, which is calculated at the nominal system frequency. It is also assumed that the transmission network is lossless, such that the total generated active power equals the total active power supplied to the loads. Equivalently, the entries of $Y$ are purely imaginary¹.

• Each generator is represented by the simplest quasi-static model: a voltage source with constant amplitude behind a series impedance [17]. The generators’ synchronous impedances are modeled as part of the transmission network. The mechanical power driving the rotor is governed by a linear droop-control mechanism, where the droop power is a decreasing affine function of the rotor frequency, with no additional dynamics².

• Each load can be decomposed into a power sink, which is connected in parallel to additional reactive components (inductors and capacitors). The main assumption for each power sink is that the consumed active and reactive powers ($P_{\ell}$ and $Q_{\ell}$) are continuously differentiable.

functions of only the local load’s voltage amplitude ($V_\ell$), and

at an equilibrium point of interest $x^*$. This is the case, for instance, if each load consumes constant active and reactive powers. The additional reactive components can be regarded as part of the lossless transmission network.

The convention that all the voltage phasor angles $\delta_2,\delta_3,\ldots,\delta_N$ are measured with respect to the first generator, so that $\delta_1=0$, is used. Also, angles that differ by an integer multiple of $2\pi$ are identified. Then, the following vectors are defined:

Here, $\omega_i(t)$ is the electric frequency of generator $i$ in rad/s, $P_i(t)$ is the generator’s active power output, $\delta_i(t)$ is the relative rotor angle as defined above, $\omega_s$ is the reference frequency, and $x(t)$ is the state vector. The system is described by the following differential and algebraic equations [17]:

The matrix $M$ represents the generators’ inertia constants, and the matrix $D$ represents the generators’ damping coefficients, for instance due to the damper windings, combined with the prime movers’ droop control constants. Both matrices are symmetric and positive definite$^3$. The vector $P_{\mathrm{ref}}$ contains the reference powers for the generators. The matrix $U$ describes the relations $\dot{\delta}_i=\omega_i-\omega_1$ for $i=2,\ldots,G$, and is given by

where $\mathbf{1}_{(G-1)\times1}$ is a vector of ones, and $I_{(G-1)\times(G-1)}$ is the identity matrix. The function $f(\cdot)$ depends on the transmission network and the loads. This function is usually not available in closed form, and has to be evaluated numerically, by solving the network’s AC power-flow equations. At an equilibrium point $x^*$ (if one exists), define

Since the transmission network is lossless, it also holds that

where $P_L$ is the loads’ total active power. Due to (1), it follows that for all $k\in\{2,\ldots,G\}$,

From (4) it follows that $\partial P/\partial\delta=U^{\mathsf{T}}H$. Based on the definitions above, the Jacobian matrix of the dynamic system $\dot{x}=F(x)$ from (2) is

Before stating the main result, it is shown that, under the assumptions listed above, the matrix $H$ is symmetric:

Proof. Define $P_i$, $Q_i$ to be the active and reactive powers injected into bus $i$, and $V_i$, $\delta_i$ to be the voltage amplitude and phase at bus $i$. In addition, denote the components of $Y$ by $y_{i,n}$. These variables are linked by the standard AC power flow equations [11, 17, 36]:

for $i=1,\ldots,N$. Define $J$ to be the Jacobian of these equations at the equilibrium point $x^*$ (see for instance [17, 36]), such that

where $\mathrm{d}P_i$, $\mathrm{d}Q_i$, $\mathrm{d}\delta_i$, $\mathrm{d}V_i$ are the differentials at $x^*$, and $V_i^*$ are the voltage amplitudes at $x^*$, for $i=2,\ldots,N$. Note that $\mathrm{d}V_1$ does not appear here, since it is assumed that $\mathrm{d}V_1=0$. Now reorder the rows and columns of $J$ to obtain

with

By straightforward calculation, it can be verified that if the admittance matrix $Y$ is purely imaginary, then the matrix $J$ is symmetric (see [17 Chapt. 9], or [36 p. 117]). The matrix $\tilde{J}$ may be obtained from $J$ via a permutation matrix $S$, such that $\tilde{J}=SJS^{\mathrm T}$, so $\tilde{J}$ is also symmetric. It follows that (6) can be rewritten as

where $\tilde{J}_{1,1}$ and $\tilde{J}_{2,2}$ are symmetric. It is well-known that $\tilde{J}_{2,2}$, being real and symmetric, can be factored like $\tilde{J}_{2,2}=W^{\mathrm T}\tilde{\Lambda}W$, where $W$ is real, $W^{\mathrm T}W=WW^{\mathrm T}=I$ (unitary matrix), and $\tilde{\Lambda}$ is diagonal, with the eigenvalues of $\tilde{J}_{2,2}$ (which are real) on its diagonal. The order of the eigenvalues in $\tilde{\Lambda}$ is arranged such that the zero eigenvalues (if they exist) are the last. Thus,

where $\Lambda$ is a real diagonal matrix of dimension $h\times h$, $h\in\{0,1,\ldots,N-G\}$, with no zeros on its main diagonal. Equation (7) is rewritten equivalently as follows:

Decompose $W=\begin{bmatrix}W_1\\W_2\end{bmatrix}$, where $W_1$ has $h$ rows. It is claimed that $W_2\tilde{J}_{1,2}^{\mathrm T}=0$. Indeed, the last $N-G-h$ rows of (8) are

Since $r_a$ is arbitrary, and $w_b=0$ due to (1), the claim follows. By symmetry, it follows that the last $N-G-h$ rows and also the last $N-G-h$ columns in the matrix on the left side of (8) are zero. Elimination of these redundant rows and columns leads to

From the last $h$ rows it is obtained that

hence

where

This matrix $\bar{H}$ is symmetric, so (9) can rewritten as

In addition, since the generators’ voltage amplitudes are constant, $(\mathrm dV_2,\ldots,\mathrm dV_G)=0$, which leads to

with a symmetric matrix $H$. ■

The main result

Consider the dynamic system (2) with two or more generators ($G\geq2$), and a symmetric matrix of partial derivatives $H$, as defined in (3), which is calculated at an equilibrium point $x^*$. The primary result is the following:

• The equilibrium point $x^*$ is asymptotically stable if the matrix $H$ is positive definite.

• The equilibrium point $x^*$ is unstable if the matrix $H$ is full-rank and not positive definite.

Theorem 2 below shows this result (implicitly), and also provides details on the relations between the eigenvalues of the matrices $A$ and $H$. Before stating this main theorem, recall some definitions and terminology from linear algebra.

Notation and terminology. Real matrices are used, but these are often applied to complex vectors (for example, to eigenvectors). For any square matrix $A$, $\sigma(A)$ denotes the set of its eigenvalues. The algebraic multiplicity of an eigenvalue $\lambda$ is the dimension of the space of all the eigenvectors and generalized eigenvectors that correspond to $\lambda$. Let $\langle\cdot,\cdot\rangle$ denote the standard inner product on $\mathbb{C}^n$: $\langle a,b\rangle=a_1\bar{b}_1+\cdots+a_n\bar{b}_n$. A real matrix $H$ of dimension $n\times n$ is called positive if $H^{\mathrm T}=H$ and for every vector $z\in\mathbb{R}^n$, $\langle Hz,z\rangle\geq0$. Equivalently, $H$ is positive if $H^{\mathrm T}=H$ (which alone implies that $\sigma(H)\subset\mathbb{R}$) and $\sigma(H)\subset[0,\infty)$. A positive $H$ is called positive definite if $0\notin\sigma(H)$ (i.e., if $H$ has full rank). Equivalently, $H$ is positive definite if $H^{\mathrm T}=H$ and for every nonzero vector $z\in\mathbb{R}^n$, $\langle Hz,z\rangle>0$.

Symbols $\mathbb{C}_-$ and $\mathbb{C}_+$ denote the open left half-plane and the open right half-plane in $\mathbb{C}$ respectively. In this paragraph, let $A$ be a square matrix with complex entries, of dimension $n\times n$. $A$ is called stable if $\sigma(A)\subset\mathbb{C}_-$ (equivalently, if $\lim_{t\to\infty}e^{At}=0$). The strictly stable subspace of $A$ is defined by

This space is the linear span of all the eigenvectors and generalized eigenvectors that correspond to an eigenvalue in $\mathbb{C}_-$. The strictly unstable subspace of $A$ is $X_+(A)=X_-(A)$. Define the notation

From $X_-(A)\cap X_+(A)=\{0\}$ it follows that

Clearly, $\operatorname{neg}(X)$ is the sum of the algebraic multiplicities of the eigenvalues of $A$ in $\mathbb{C}_{-}$, while $\operatorname{pos}(A)$ is the sum of the algebraic multiplicities of the eigenvalues of $A$ in $\mathbb{C}_{+}$. If $A$ is the Jacobian of a dynamic system at an equilibrium point $x^{*}$, then the following facts are known [43]: If $A$ is stable then $x^{*}$ is asymptotically stable. If $A$ has an eigenvalue in $\mathbb{C}_{+}$, then $x^{*}$ is not (Lyapunov) stable.

Remarks concerning the theorem:

• If $H > 0$ then $A$ is stable. This follows from 2) by considering $\operatorname{pos}(H) = G - 1$, which gives $\operatorname{neg}(A) = n$.

• If $H$ is full-rank and not positive definite ($H \not> 0$), then $A$ is unstable, and has an eigenvalue $\lambda\in\mathbb{C}_{+}$. This follows from 2), considering the case $\operatorname{neg}(H) \ge1$.

• Concerning the stability of the equilibrium point $x^{*}$ for the nonlinear system (2): If $H > 0$, then $x^{*}$ is (locally) asymptotically stable. If $\operatorname{neg}(H) > 0$, then $x^{*}$ is not even Lyapunov stable. This follows from what was said just before the theorem.

• If $H$ is full-rank, then it is the only information needed to determine the stability of $A$ and $x^{*}$.

• Physically, the matrix $H$ represents the synchronizing power coefficients of the network, hence, the condition $H > 0$ ensures that the network generates restoring torques for any rotor angle deviation.

Proof. The proof of claim 1) is trivial: If $G = 1$, then the only state variable is $\omega_{1}$, while $M$ and $D$ are positive numbers and $A = -D/M$, which is clearly stable.

Next, consider claims 3) and 4), as they are required to prove claim 2). Claim 3) is proved first: Let $\lambda$ be an eigenvalue of $A$, and let $q = \begin{bmatrix}v\\w\end{bmatrix}$ be the corresponding eigenvector. Then

In case $v = 0$, this leads to $\lambda w = 0$. Since by definition $q \ne0$, it follows that $w \ne0$, so $\lambda= 0$.

In case $v \ne0$, multiply the first equation in (10) with $\lambda$, and use the second equation to obtain

Take the inner product with $v$ to obtain a quadratic equation in $\lambda$:

Assumed that $\lambda$ is purely imaginary, so

Since $M$, $D$, and $H$ are symmetric, the imaginary part of this equation is $|\lambda|\langle Dv,v\rangle= 0$. But $D$ is positive definite and $v \ne0$, so it follows that $\lambda= 0$.

Now claim 4) is addressed. First, it is shown that if $H$ has an eigenvalue 0, then $A$ has an eigenvalue 0. To this end, let $z$ be an eigenvector of $H$ with eigenvalue 0, such that $Hz = 0$, and let $q = \begin{bmatrix}0\\z\end{bmatrix}$. Clearly $Aq = 0$, so $0 \in\sigma(A)$.

Now, the other direction in claim 4) is proven, showing that if $0 \in\sigma(A)$, then $0 \in\sigma(H)$. So assume that $0 \in\sigma(A)$ and let $q = \begin{bmatrix}v\\w\end{bmatrix}$ be a corresponding eigenvector. Then

Since $D$ is invertible, this leads to $w \ne0$, and

Take inner product with $Hw$:

But $D^{-1} > 0$, hence $U^{\mathsf{T}}Hw = 0$. The only solution for this equation is $Hw = 0$, so $0 \in\sigma(H)$.

To prove 2) consider $G \ge2$ and $0 \notin\sigma(H)$. It can be verified directly that the Lyapunov equation

is satisfied by

Denote $V_{-} = X_{-}(P)$ and $m = \operatorname{neg}(H) = \operatorname{neg}(P) = \dim V_{-}$. Similarly, denote $V_{+} = X_{+}(P)$, so that $\dim V_{+} = n - m$. For any nonzero $x_{0} \in V_{-}$, $\langle Px_{0},x_{0}\rangle< 0$. Denoting $x(t) = e^{At}x_{0}$, then

so that $\langle Px(t),x(t)\rangle\le\langle Px_{0},x_{0}\rangle< 0$, which shows that $x(t)$ does not converge to 0. In other words, $V_{-} \cap X_{-} = \{0\}$, which implies that $\dim V_{-} + \dim X_{-} \le n$, so that

Now, a similar argument is used for a nonzero $x_{0} \in V_{+}$, obtaining that $x(t) = e^{-At}x_{0}$ does not converge to zero, so that $V_{+} \cap X_{+} = \{0\}$. From here it is obtained that

From claim 4) it is known that $0 \notin\sigma(A)$, and hence, using claim 3), it may be concluded that $A$ has no imaginary eigenvalues. This implies that

Putting together this last equality with (12) and (13), it follows that $\dim X_{-} = n - m$ and $\dim X_{+} = m$.

Generalization to arbitrary loads

In this section the case when the loads are not assumed to take constant power is considered. An extension of Theorem 2 to the case where the assumption that the loads in the grid consume locally constant active and reactive powers is dropped as expressed by (1) is derived. The other assumptions and the notation from Sect. II are kept.

The matrix $\mathcal{H}$ can still be defined as in (3), but the formula (4) no longer holds. This time, $\mathcal{H}$ will not be symmetric, although numerical experimentation with reasonable circuits shows that it will be very nearly symmetric.

The full Jacobian matrix is introduced to describe the sensitivity of the (active) electric powers with respect to the relative rotor angles at the equilibrium point $\boldsymbol{x}^{*}$:

The last $G-1$ rows of $\mathcal{H}$ are the matrix $H$:

Following an argument very similar to the one leading to (5), it follows that in this new context (i.e., without assuming (1)), the matrix of the linearized system is

Note that if (1) holds, then (4) implies that $\mathcal{H}=U^{\mathsf{T}}H$, and the above $A$ reduces to the one from (5). Thus, the local stability of the system (2) around the equilibrium $\boldsymbol{x}^{*}$ follows if the above $A$ is stable.

Proposition 3. The matrix $A$ from (15) is stable if and only if det $A\ne0$ and the nonzero solutions of the equation

are in the open left half-plane $\mathbb{C}_{-}$.

Proof. Assume that $\det A\ne0$ and all the solutions of (16) are in $\mathbb{C}_{-}$. Let $s\in\mathbb{C}$ be an eigenvalue of $A$ and let $\boldsymbol{x}=\begin{bmatrix}\omega\\\delta\end{bmatrix}$ be a corresponding eigenvector, then $\boldsymbol{x}\ne0$, $s\ne0$ and

Multiply the first equation above with $s$ and substitute $s\delta$ from the second equation, obtaining

Since $\boldsymbol{x}\ne0$, it follows that $\omega\ne0$, so that the above equality implies (16); thus, by assumption, $s\in\mathbb{C}_{-}$. Hence, $A$ is stable. Conversely, assume that $A$ is stable, then clearly $\det A\ne0$. Let $s$ be a nonzero solution of (16). Then there exists a vector $\omega\in\mathbb{C}^{G}$ such that $\omega\ne0$ and (17) holds. Denote $\delta=\frac{1}{s}U\omega$, then it is easy to see $\boldsymbol{x}=\begin{bmatrix}\omega\\\delta\end{bmatrix}$ is an eigenvector of $A$, corresponding to the eigenvalue $s$. Since $A$ is stable, it follows that $s\in\mathbb{C}_{-}$, as claimed. ■

In the following proposition sufficient conditions for the stability of $A$ are derived.

Proposition 4. Assume that at least one of the following conditions holds:

1. Locally constant power loads, i.e., (1) holds.

2. $M^{-1}D=\gamma I$, where $\gamma$ is a positive number.

Then the matrix $A$ from (15) is stable if

Proof. First, assume that condition (1) holds. Then (from (4)) $\mathcal{H}=U^{\mathsf{T}}H$ and (from Theorem 1) $H^{\mathsf{T}}=H$. Then, the condition (18) becomes

Notice that the above implies that $\operatorname{Ker}H=\{0\}$, so that $\det H\ne0$. Denoting $\overline{M}=UM^{-1}U^{\mathsf{T}}$, it follows that $\overline{M}>0$, because $M>0$ by assumption and $\operatorname{Ker}U^{\mathsf{T}}=\{0\}$. Denote the nonzero spectrum of a square matrix $T$ by $\sigma_{0}(T)$. It is well known from linear algebra that if $S$ and $T$ are matrices such that $ST$ and $TS$ are both square matrices (possibly of different dimensions), then $\sigma_{0}(ST)=\sigma_{0}(TS)$. It follows that the condition (18) is equivalent to $\sigma_{0}(\overline{M}^{\frac{1}{2}}H\overline{M}^{\frac{1}{2}})\subset(0,\infty)$. Since the matrices $\overline{M}^{\frac{1}{2}}$ and $H$ have nonzero determinant, $\det\overline{M}^{\frac{1}{2}}H\overline{M}^{\frac{1}{2}}\ne0$, it follows that $\sigma(\overline{M}^{\frac{1}{2}}H\overline{M}^{\frac{1}{2}})\subset(0,\infty)$. Since $\overline{M}^{\frac{1}{2}}H\overline{M}^{\frac{1}{2}}$ is symmetric, this is equivalent to $H>0$. According to Theorem 2, $A$ is stable.

Now assume that condition (2) holds. Let $\begin{bmatrix}\omega\\\delta\end{bmatrix}$ be an eigenvector of $A$ corresponding to the eigenvalue $s$. Then

Multiplying the first equation by $U$ from the left and substituting the second equation results in

Assuming the special case where $M^{-1}D=\gamma I$ for some number $\gamma>0$, the equation becomes

Using the second equation in (19) again, this may be written as

Rearranging the terms yields

This equation shows that for any non-trivial eigenvector with the corresponding eigenvalue $s$,

Let $\alpha$ be any such eigenvalue of $UM^{-1}\mathcal{H}$, which means that $s^{2}+s\gamma+\alpha=0$. Assume (18), then $\alpha>0$. Since also $\gamma>0$, the roots of $s^{2}+s\gamma+\alpha=0$ are in $\mathbb{C}_{-}$, which means that $A$ is stable. ■

Numerical experiments on a variety of networks with various resistive loads and $\gamma$ values, confirm that the spectrum of the matrix $UM^{-1}\mathcal{H}$ is real and positive for many stable systems as long as the condition $M^{-1}D=\gamma I$ holds. It seems that the condition $M^{-1}D = \gamma I$ implies that $\sigma(UM^{-1}H)$ is real, but there is no formal proof.

Proof. By assumption, $M^{-1}D = \gamma I$, with $\gamma> 0$. Let $s$ be the positive solution of the equation

To prove the proposition, it suffices to show that there exists $\begin{bmatrix}\omega\\\delta\end{bmatrix}$ with $\omega\in\mathbb{R}^{G}$, $\delta\in\mathbb{R}^{G-1}$, $\delta\ne0$, such that

Let $\delta\in\mathbb{R}^{G-1}$ be such that $UM^{-1}H\delta= \alpha\delta$, $\delta\ne0$. Since $U$ is onto, there exists $\omega\in\mathbb{R}^{G-1}$ such that $U\omega= s\delta$. There are many such vectors $\omega$, since $U1_G = 0$, where $1_G$ is the vector in $\mathbb{R}^{G}$ with all entries equal to 1. Thus, any value can be chosen for the first component of $\omega$, and it can be chosen such that

Equation (20) is now multiplied by $\delta$ and rewritten as follows:

Using once more that $U\omega= s\delta$, results in

The $(j-1)$-th row of the above ($j \in\{2,3,\ldots G\}$) is

Combining this with (22), yields

for all $j \in\{2,3,\ldots G\}$. This fact, together with (22), means that

This, together with $U\omega= s\delta$, means that (21) holds. ■

Comparison to related previous results

Following is a comparison to some previous works that explore the stability of similar dynamic models of multi-machine grids (each generator being described by a second order swing equation) and give algebraic criteria for local stability. The proposed method is compared with three significant references:

• Machowski, Bialek and Bumby [36] from 2012,

• Dörfler and Bullo [12] from 2012,

• Wu and Liu [13] from 1986.

Consider the stability criterion given in [36 Chapter 12]. The end result provided in [36] is similar to the one shown in the current paper: the power system dynamic model is stable if the matrix $H$ is positive definite. However, this result is obtained under more restrictive assumptions:

• The matrix $M$ of inertia constants and the matrix $D$ of damping coefficients are assumed in [36] to be diagonal and normalized, such that $D_i/M_i$ is constant (independent of $i$). This work eliminates this assumption, and it is only assumed that the matrices $D$ and $M$ are positive definite (not even diagonal).

• More importantly, careful analysis shows that in fact, the model in [36 Chapter 12] corresponds to the model used in this work without any load bus ($L=0$). This absence of load buses is not stated explicitly, and is revealed only after some reasoning. Indeed, in [36 Chapter 12] it is assumed that the network loads are represented by constant impedances, which are lumped within the transmission network model, and expressed via the network admittance matrix $Y$. It is then assumed that the transmission network, including the loads, is lossless. This assumption is more restrictive than the one used in the proposed framework, since here it is only assumed that the transmission network without the loads is lossless. This assumption is essential for ensuring the symmetry of the stability matrix, because resistive losses render the matrix $H$ non-symmetric. Thus, the assumptions in [36] enable only two types of possible loads: either the loads have no resistive components, and therefore consume no active power (so that they may be regarded as a part of the transmission network), or they are located directly on the generator buses, and therefore do not affect the system dynamics. Thus, the result in Section III may be regarded as a generalization of the one in [36 Chapter 12], since locally constant power loads are allowed to appear anywhere within the network (in other words, arbitrary $L \ge0$ are allowed instead of $L=0$).

Still as part of the comparison with [36], it is noted that the proposed stability criterion is both necessary and sufficient (under the assumptions stated in Section II), since for a full-rank matrix $H$, an equilibrium point is stable if and only if $H$ is positive definite. In [36], only the ‘if’ part is stated. Several other properties of the eigenvalues of $H$ are also established.

Next, consider the stability result in the influential paper [12]. The main result concerns the well-known NRPS model mentioned earlier, where each generator is represented by a second order swing equation, and they are interconnected with a linear, balanced and lossy electric network containing also the linear loads. The derivation of the model relies on Section 3 of [8]. The transmission network together with the loads is represented by its (symmetric) admittance matrix $Y$, which is regarded as a constant complex matrix, because the frequencies are assumed to be in a narrow range around the nominal frequency (this assumption about $Y$ is made also in [13, 36]). The idea in [12] is to consider the limit of this system when the inertia constants of the generators vanish, called the nonlinear Kuramoto (NK) model, and then show that for small enough inertia constants and for initial states in a neighborhood of a certain “preferred” equilibrium point (the equilibrium point that corresponds to the smallest relative angles), the state trajectories of the NRPS model remain close to the state trajectories of the NK model, which in turn is locally stable. A corrected, strengthened and generalized version of [12] has been published in [8], with a detailed explanation of what the improvements are. For clarity, the achieved results are compared with those in [8].

A crucial fact about the result in [8] concerning the NRPS model is that the inertia constants are assumed to be of the form $M_j = \varepsilon m_j$, where $m_j > 0$ is fixed and $\varepsilon> 0$ is a small parameter, and all the results hold for initial states in a neighborhood of the preferred equilibrium, and only for sufficiently small $\varepsilon$, the upper bound on $\varepsilon$ depending on the chosen neighborhood. Two quantities are introduced, denoted by $\Gamma_{\min}$ and $\Gamma_{\mathrm{crit}}$, which depend on the parameters of the system (but not on the inertia constants), and all the results are given under the assumption that $\Gamma_{\min} > \Gamma_{\mathrm{crit}}$. This assumption is rather restrictive, according to the authors’ experience. If this assumption is satisfied, then according to part (4) of Theorem 4.1, the preferred equilibrium point $x^*$ is exponentially attractive, i.e., all the state trajectories in certain neighborhoods of $x^*$ converge to $x^*$ at an exponential rate, if $\varepsilon$ is sufficiently small. Compared to the achieved result, the differences are: (a) linear loads, as opposed to arbitrary loads, as discussed in Section IV, (b) there is no requirement for an assumption resembling $\Gamma_{\min} > \Gamma_{\mathrm{crit}}$, (c) the inertias are not assumed to be small, (d) necessary and sufficient conditions are given (for local asymptotic stability), as opposed to sufficient conditions.

To empirically validate the theoretical differences between the proposed criterion and the approach described in [8], both methods were tested on standard IEEE benchmark test cases with constant power loads. Table 1 demonstrates the fundamental limitation of the network reduction approach: while the $H$ matrix method consistently provides definitive stability assessment across all cases, the method described in [8] becomes increasingly unreliable as system complexity grows. Indeed the sufficient condition $\Gamma_{\min} > \Gamma_{\mathrm{crit}}$ successfully guarantees synchronization for $\varepsilon$ small enough only for small, tightly-coupled networks, such as “case30” and “case30”. However, for larger systems, the method in [8] provides no result, as observed in “case118”, “case1354pegase” or “case2869pegase”. The condition $\Gamma_{\min} > \Gamma_{\mathrm{crit}}$ becomes overly conservative and impractical for large, sparse power grids. In contrast, the proposed approach maintains numerical stability and provides definitive stability assessment.

test case	buses	num. of gen/loads	Proposed method	method in [8, 12]	$\Gamma_{\min}/\Gamma_{\mathrm{crit}}$
9	9	3/6	stable	$\Gamma_{\min} > \Gamma_{\mathrm{crit}}$, stable for $\varepsilon$ small enough	4.50
14	14	5/9	stable	$\Gamma_{\min} > \Gamma_{\mathrm{crit}}$, stable for $\varepsilon$ small enough	2.65
30	30	6/24	stable	$\Gamma_{\min} > \Gamma_{\mathrm{crit}}$, stable for $\varepsilon$ small enough	1.18
118	118	54/64	stable	$\Gamma_{\min} < \Gamma_{\mathrm{crit}}$	0.18
1354pegase	1354	260/1094	stable	$\Gamma_{\min} < \Gamma_{\mathrm{crit}}$	0.12
2869pegase	2869	510/2359	stable	$\Gamma_{\min} < \Gamma_{\mathrm{crit}}$	0.04

If $\Gamma_{\min} \le\Gamma_{\mathrm{crit}}$, then [8, 12] give no stability result.

Finally, the results are compared with [13], which at its time has been a significant advance in the field, by proposing to look at power system models using systems of differential and algebraic equations and finding an algebraic condition for local stability. Their model is similar to the one used in this paper, with a lossless transmission network and with constant generator terminal voltage amplitudes, but also more complicated: [13] distinguishes between the generator internal buses (where the synchronous internal voltage is) and the generator terminal buses (In this work, only the internal buses are considered, while the terminal buses and the stator impedance are treated as part of the transmission network). Paper [13] assumes that the loads consume active power that is an affine function of the local frequency (equation (34)) and reactive power that is a constant plus the reactive power consumption of a shunt capacitor (equation (35)). They also assume that the reactive power output at each generator internal bus is constant. This is not intuitive, since the generator has to maintain constant terminal voltage amplitude and a given active power, and this seems to be one constraint too many on each generator.

In addition, [13] treats every voltage (at every node that is not a generator terminal) as a state variable. Moreover, they treat the phase angle at every node, except one, that is the reference, as a state variable. This leads to more than two times more state variables than what is needed, see equation (28). With these state variables, the system is described by a collection of differential and algebraic equations, (36) and (37). Under an additional mild assumption (A7), the reference [13] computes the linearization of the aforementioned differential-algebraic system around an equilibrium point $x$ (equations (41) and (42)) and defines a square matrix $\tilde{F}(x)$, which is proved to be symmetric (Lemma 1). The matrix $\tilde{F}(x)$ is of dimension $n-1$, which is simply $2G+L-1$. The main result of [13] (Theorem 1) has some formal resemblance with the main result of this work, namely, the equilibrium point is (locally) stable if and only if $\tilde{F}(x) > 0$. However, the dimension of matrix $\tilde{F}(x)$ is at least twice the dimension of the matrix $H$, and its definition is more complex. In addition, there is no information in [13] about the number of unstable eigenvalues of the linearized system, in case it is unstable.

Case study 1: Loads located on the same buses as generators

A number of case studies follow, demonstrating the usefulness of the stability criterion in Theorem 2. The first case study is a famous abstract power system model, in which the system is divided into “areas”, and each area includes one generator and one load. This implies that the loads are located on the same buses as the generators, see Figure 2.

This simplified model has been studied extensively in both the power systems and control communities for many years. In the control community, this model is sometimes described as an instance of the NRPS model, or analyzed in the context of a Kuramoto model [12, 25]. In the power systems community, this model appears in classic textbooks, for instance [36]. Since the mathematical properties of this model are already well known [12, 27, 36], the main objective in this section is not to learn something new about it, but rather to illustrate and confirm the proposed stability criterion for this particular class of models.

For this system, the function $f(\cdot)$ from (2) can be expressed by simple formulas:

for $i=1,\ldots,G$, where $y_{i,n}$ are elements of the admittance matrix, $V_i$ are the voltage amplitudes of the generators, and $P_{L,i}$ are the active powers of the loads. For the more specific case of a lossless transmission network that is primarily inductive⁴, these equations may be also written as

where $a_{i,n}=a_{n,i}=V_iV_n|y_{i,n}|$. Computing the derivatives of these equations yields

where $i,k\in\{2,\ldots,G\}$, and the star notation marks the values at the equilibrium point. The equality $a_{i,k}=a_{k,i}$ implies $h_{i,k}=h_{k,i}$, so the matrix $H$ is symmetric, as expected.

According to the proposed stability criterion, for $x^*$ to be asymptotically stable, the matrix $H$ must be positive definite. To this end, consider the quadratic form $z^{\mathsf{T}}Hz$, where $z=(z_2,\ldots,z_G)$ is any real vector:

Define $b_i=a_{i,1}\cos(\delta_i^*)$ to obtain $h_{i,i}=b_i-\sum_{\substack{k=2\\k\ne i}}^{G}h_{i,k}$, and substitute to get

If this expression is $>0$ for every $z\ne0$, then $H>0$.

For $G\ge2$, a network is called connected⁵ if for every vector $z=(z_2,\ldots,z_G)^{\mathsf{T}}\ne0$, it is true that

Assume that indeed $G\ge2$ and the network is connected, as above. In this case, it can be shown that

for every $z\ne0$. When comparing (26) to (24), it is evident that since the coefficients $a_{i,k}$ are all non-negative, the expression in (24) is strictly positive for every $z\ne0$, if $\cos(\delta_i^*)>0$, $\cos(\delta_i^*-\delta_k^*)>0$, for $i,k=2,\ldots,G$. As a result $H$ is positive definite, and the equilibrium point $x^*$ is asymptotically stable, if

More formally, consider an undirected weighted graph with adjacency matrix $A$ (see, for instance, [44]), whose elements are $a_{i,k}\ge0$ for $i,k=1,\ldots,G$, with $i\ne k$, and its main diagonal is zero. This graph is connected if and only if (25) holds for every vector $z=(z_2,\ldots,z_G)^{\mathsf{T}}\ne0$. Here is a quick informal proof for the “if” part of the claim: Assume that (25) does not hold. By negation that there is a group of nodes that form an island within the network, in the sense that they are connected to each other, but none of them is connected directly or indirectly to the node 1. Let $B$ be the set of indices of the nodes within this group. Define a vector $z$ as follows: $z_i=1$ if $i\in B$ and $z_i=0$ otherwise, for $i=2,\ldots,G$. Clearly this vector zeros the sum in (25), which contradicts the definition of a connected graph. Therefore, a group of nodes that form an island cannot exist. The “only if” part of the claim can be proved based on a similar reasoning.

for $i,k=2,\ldots,G$. Other stable equilibrium points can be obtained by shifting any angle by $2\pi$. Note that [8] and [12] refer to the same set of angles as in (27), using the notation $\delta\in\Delta_{\mathrm{grnd}}(\pi/2)$. The resulting region is plotted for the case of $G=3$ in Figure 3.

Case study 2: Stability under the DC power-flow approximation

The “DC power flow” approximate model of the power grid is obtained when using a small angle approximation of (23):

for $i=1\ldots G$, where $a_{i,n}=a_{n,i}=V_iV_n\lvert y_{i,n}\rvert$, and $P_{L,i}$ is the equivalent load connected at bus $i$. Straight-forward differentiation yields

where $i,k=2,\ldots,G$. Again, since $a_{i,k}=a_{k,i}$ it follows that $h_{i,k}=h_{k,i}$, so $H$ is symmetric. Consider now the quadratic form $z^{\mathsf{T}}Hz$, where $z=(z_2,\ldots,z_G)^{\mathsf{T}}$ is any real vector. This can be expressed as

If $G\ge2$, and the transmission network is connected (as defined in Section VI), then the expression in (28) is strictly positive for every $z\ne0$, which means that $H$ is positive definite. This follows from the fact that all the coefficients $a_{i,k}$ are non-negative. Because $H$ is positive definite, the equilibrium point $x^*$ is asymptotically stable, and since the dynamic system is linear in this example, the equilibrium point $x^*$ is also globally asymptotically stable.

Case study 3: Numerical results for standard ``MATPOWER'' case studies

In this section the proposed stability criterion is applied to several standard case-study networks, which are taken from the software package MatPower [45] with double-precision arithmetic. The considered case-studies are “case9”, “case14”, “case30”, “case118”, and “case1354pegase”, and “case2869pegase”. For these case-studies the matrix $H$ is computed numerically, and it is shown that in all cases $H$ is symmetric and positive-definite at the default equilibrium point $x^*$, implying that $x^*$ is asymptotically stable. The matrix $H$ is also used to partially explain what specific mathematical property of these networks makes them stable. This is done by showing that most of these networks lead to a matrix $H$ that is diagonally dominant, see explanations below.

The computation of the $H$ matrix follows this procedure: First, each transmission line’s series resistance and bus shunt conductance are set to zero (BR_R = 0, GS = 0) to ensure the lossless network assumption. The equilibrium point $x^*$ is obtained using MatPower’s Newton-Raphson power flow solver (runpf function) with default convergence tolerances ($10^{-8}$ for power mismatches). All generator buses are set as $V$-$\delta$ (reference) buses to maintain constant voltage amplitudes.

The $H$ matrix elements are computed using finite differences: each generator angle $\delta_k$ ($k=2,\ldots,G$) is perturbed by $\varepsilon=10^{-3}$ radians, the power flow is re-solved, and the partial derivatives are estimated as $h_{i,k}=\partial P_i/\partial\delta_k\approx(P_i^{\mu}-P_i^*)/\varepsilon$, where $P_i^{\mu}$ and $P_i^*$ are the updated and equilibrium generator powers, respectively. This perturbation size provides sufficient numerical accuracy while avoiding finite-precision errors. The symmetry of $H$ is verified by checking that $\lVert H-H^{\mathsf{T}}\rVert_\infty/\lVert H\rVert_\infty<10^{-4}$. For systems with fewer than 200 generators, positive definiteness is verified using the $\operatorname{eig}()$ function, for larger systems, Cholesky factorization is used as a computationally efficient alternative.

For each case-study the matrix $H$ is computed numerically, based on the following procedure: first, the series resistance in every branch and the shunt conductance on every bus are zeroed, to make the transmission network lossless. At the second stage, the default equilibrium point $x^*$ is computed by running MatPower built-in functions. A loop is then executed in which each $\delta_k$ is slightly perturbed by the substitution $\delta_k\to\delta_k+\varepsilon$, for $k=2,\ldots,G$, with a default value of $\varepsilon=10^{-3}$. For each such change the power flow is solved again, to compute the updated generator powers $P_i^{\mu}$. The partial derivatives are then estimated using the approximation $h_{i,k} = \partial P_i/\partial\delta_k = (P_i^{u} - P_i^{*})/e$. After the matrix $H$ is computed, its symmetry is verified. This numeric process is very fast, and takes a few seconds at most for all the examined case-studies.

The equilibrium point $x^{*}$ is calculated using the Newton Raphson power flow solver provided in the MATPOWER software package. A convergence tolerance of $10^{-8}$ for power mismatches is utilized to ensure numerical precision. For the assessment of stability, the eigenvalues of the matrix $H$ are computed using the standard eig function for smaller power systems. In contrast, the Cholesky factorization is employed for larger networks such as the 2869-bus Pegase system to determine positive definiteness with high computational efficiency.

For the first case-study, “case9”, the resulting matrix is

This matrix is symmetric, with eigenvalues 14.52 and 4.24, and therefore it is positive-definite, which means that the default equilibrium point $x^{*}$ is asymptotically stable.

Another reason that makes this matrix positive definite is that it is strictly diagonally dominant (SDD). Recall that a real matrix $H$ is SDD if $|h_{i,i}| > \sum_{j \ne i} |h_{i,j}|$. If a real symmetric matrix has this property SDD, and the entries on its main diagonal are positive, then this matrix is also positive-definite ($H > 0$), see [46 Theorem 6.2.27].

For “case9”, the matrix $H$ shown above is obviously SDD, so $H > 0$. This may be qualitatively explained as follows: if for every generator “the change in local power is dominant”, in the sense that $|h_{i,i}| > \sum_{j \ne i} |h_{i,j}|$, then $H$ is positive definite, and the equilibrium point $x^{*}$ is stable. Intuitively, this property will be true for many power systems, and may partially explain, based on the proposed theorem, what makes typical power systems stable.

For the next case-study, “case14”, a similar behaviour is observed. The matrix $H$ is given by

The eigenvalues are 59.49, 15.02, 10.66 and 3.88, so this matrix is positive-definite, and $x^{*}$ is stable. As observed previously, the matrix $H$ is SDD, which explains why the equilibrium point is stable.

For “case30”, the matrix $H$ is given by

The eigenvalues are 41.25, 16.14, 9.94, 6.46 and 2.06, so this matrix is positive-definite, and $x^{*}$ is stable. The matrix $H$ is again SDD. For “case118” the matrix $H$ is schematically plotted in Figure VIII. As shown in the figure, $H$ has a dominant main diagonal, but in this case it is not SDD, so the intuition provided above is not valid. Still, the eigenvalues of $H$ are all positive, so this matrix is positive-definite, and $x^{*}$ is stable.

Vector (PDF)

For “case1354pegase” the matrix $H$ is schematically plotted in Figure 5, in logarithmic scale. As in the previous case, $H$ has a dominant main diagonal, but is not SDD. Nevertheless, it can be verified using Cholesky decomposition that this matrix is positive-definite, and $x^{*}$ is stable.

Vector (PDF)

Again, the “case2869pegase” (see Figure 6) $H$ has a dominant main diagonal, but is not SDD. Nevertheless, it can be verified using Cholesky decomposition that this matrix is positive-definite, and $x^{*}$ is stable. A summary of properties of the matrix $H$ for the above case-studies is in Table 2.

Vector (PDF)

case:	9-bus	14-bus	30-bus	118-bus	2869pegase
$\mathcal{H} > 0$	✓	✓	✓	✓	✓
$\mathcal{H}$ is SDD	✓	✓	✓	X	X

The system stability was examined under scenarios typical for inverter-based generation. The computations involved modifying the standard test cases to include random active power injections and reactive power deviations, examining the behavior of distributed renewable resources. In all examined scenarios the $H > 0$ and $H$ is SDD.

As detailed in Section IV, the stability analysis presented in this paper extends beyond the constant power load assumption to more general load models. When the damping and inertia matrices satisfy $M^{-1}D = \gamma I$ with $\gamma> 0$, Proposition 4 provides an alternative stability criterion through the spectrum of $UM^{-1}\mathcal{H}$. As noted previously, a predominant finding based on extensive numerical experiments is that the spectrum of $UM^{-1}\mathcal{H}$ is real and positive for all the stable systems that were examined.

To validate Proposition 4, the same MatPower case studies have been re-examined with modified load characteristics. Rather than maintaining constant power consumption, the loads now exhibit voltage-dependent behavior while $M^{-1}D = \gamma I$ is preserved. The computational procedure remains similar to the constant power case, with finite difference calculations of the $\mathcal{H}$ matrix elements following generator angle perturbations, and resistive loads (in addition to constant power loads) are present on each load bus.

For the “case9” network under these modified conditions, the eigenvalues of $UM^{-1}\mathcal{H}$ are 12.72 and 14.52. Despite the $\mathcal{H}$ matrix losing positive definiteness (it is no longer symmetric), the system remains asymptotically stable. The Jacobian matrix A exhibits eigenvalues at $-1.0$ and complex pairs at $-0.5 \pm3.53i$ and $-0.5 \pm3.78i$, confirming stability.

The “case14” and “case30” networks display consistent stability under both criteria. For “case14”, the eigenvalues of $UM^{-1}\mathcal{H}$ are 85.33, 22.86, 12.93, and 7.31, while $\mathcal{H}$ is nearly self-adjoint, with eigenvalues 59.47, 15.02, 10.65, and 3.88. Similarly, for the “case30” the resulting eigenvalues of $UM^{-1}\mathcal{H}$ are 70.93, 16.62, 9.95, 7.11, and 6.32. $\mathcal{H}$ is once again nearly self-adjoint, and it has positive eigenvalues.

These numerical experiments confirm that Proposition 4 encompasses a broader class of power system models than Theorem 2. The stability condition based on $UM^{-1}\mathcal{H}$ applies to systems that satisfy $M^{-1}D = \gamma I$ (which is a common assumption, as discussed in [47]) regardless of whether loads satisfy the constant power assumption.

Discussion and Conclusions

This paper extends existing knowledge on the small-signal stability of time-varying phasor models of multi-machine AC power systems, offering a generalized $\mathcal{H}$ matrix stability criterion. For two or more generators with a full-rank $\mathcal{H}$, an equilibrium point $x^*$ is asymptotically stable if $\mathcal{H} > 0$, and unstable otherwise. This result either differs from or generalizes existing literature, as detailed in § 7. Most notably, the proposed criterion enables constant power loads to be located anywhere within the network, even without explicit analytic power flow equations. The main theorem is broad enough to predict stability in other classic models, such as certain “NRPS” model cases (Section VI) and under the DC power-flow approximation (Section VII). Another small advantage of the proposed theorem is that it facilitates faster stability assessments in large-scale power systems. This is because determining if $\mathcal{H} > 0$ is numerically simpler than computing all eigenvalues of the system matrix A. This efficiency stems from two factors: first, positive definite testing is generally more efficient in terms of run-time and memory, with methods like the Cholesky factorization being particularly effective for large, sparse matrices, as demonstrated in Table 3. Second, the $\mathcal{H}$ matrix’s dimension, $(G-1) \times(G-1)$, is about half of that of A, allowing numerical algorithms to operate faster on $\mathcal{H}$. Based on several case studies conducted, the advantages of the proposed stability criterion are demonstrated. As established in the introduction, the generalized $\mathcal{H}$ matrix stability criterion accurately identifies asymptotic stability while some alternative methodologies do not yield unequivocal results. Moreover, in contrast to approaches that rely on the full Jacobian matrix, using the $\mathcal{H}$ matrix, which is approximately half the size of the system Jacobian, leads to decreased computational time and lower memory usage.

Case	Dimensions	Comp. Time (s)	Memory (bytes)	Dimensions	Comp. Time (s)	Memory (bytes)
Case9	$5 \times5$	0.0001	200	$2 \times2$	0.0001	32
Case14	$9 \times9$	0.0001	648	$4 \times4$	0.0001	128
Case30	$11 \times11$	0.0001	968	$5 \times5$	0.0001	200
Case118	$107 \times107$	0.0027	91,592	$53 \times53$	0.0001	22,472
Case2869pegase	$1019 \times1019$	0.3756	8,306,888	$509 \times509$	0.0011	2,072,648
Case1354pegase	$519 \times519$	0.0846	2,154,888	$259 \times259$	0.0004	536,648

Elinor Ginzburg-Ganz

ELINOR GINZBURG-GANZ received her BSc from the Faculty of Electrical and Computer Engineering at the Technion in 2023. She is currently pursuing her PhD on a direct track in the Faculty of Electrical and Computer Engineering, also at the Technion. Her research interests include energy and power systems, optimal control, and machine learning. Ms. Ginzburg Ganz was the recipient of The Vatat Student Research Award in the Field of Data-Sciences in 2024, and The Freud Award for Excellence in Research in 2025.

Yoash Levron

YOASH LEVRON (SM’18) is an associate professor at the Technion—Israel Institute of Technology and a faculty member since 2014. His research interests are power systems, focusing especially on dynamics, control and artificial intelligence, and power electronics. He has published over 100 papers in professional journals and won several awards, among them the Masao Horiba award, best-of-the-best paper award in the IEEE PES General Meeting Conference, the best paper award (second place) in the IEEE Transactions on Power Electronic, the Uzi and Michal Halevy Award for Innovative Applied Engineering, the Technion Commendation for Excellence in Teaching, and the Israel Defence Prize.

Juri Belikov

JURI BELIKOV (SM’19) received his B.S. degree in mathematics from Tallinn University, and his M.S. and Ph.D. degrees in computer and systems engineering from Tallinn University of Technology in 2006, 2008, and 2012, respectively. In the years 2015-2017 he was a Postdoctoral Fellow in the Faculty of Mechanical Engineering and Andrew & Erna Viterbi Faculty of Electrical Engineering, Technion—Israel Institute of Technology. He is currently an associate professor in the Department of Software Science, Tallinn University of Technology. His main research interests lie on the edge of nonlinear control theory, power systems, and computer science.

George Weiss

GEORGE WEISS received the MEng degree in control engineering from the Polytechnic Institute of Bucharest, Romania, in 1981, and the PhD degree in applied mathematics from Weizmann Institute, Rehovot, Israel, in 1989. He has worked in many places, including the University of Exeter, UK, and Imperial College London, UK. His research interests include distributed parameter systems, operator semigroups, power electronics, and the transient stability of the power grid.