Rotor-angle versus voltage instability in the third-order model for synchronous generators

TL;DR

This paper analyzes the stability interplay between rotor-angle and voltage in power systems using a third-order generator model, providing stability conditions and a global stability map for a two-bus system.

Contribution

It offers necessary and sufficient stability conditions and analytical insights into the stability routes of the third-order generator model.

Findings

Derived stability conditions for the third-order model

Identified routes to rotor-angle and voltage instability

Created a global stability map for a two-bus system

Abstract

We investigate the interplay of rotor-angle and voltage stability in electric power systems. To this end, we carry out a local stability analysis of the third-order model which entails the classical power-swing equations and the voltage dynamics. We provide necessary and sufficient stability conditions and investigate different routes to instability. For the special case of a two-bus system we analytically derive a global stability map.

Equations

$$\dot{\delta }$$

$$M\dot{\omega }$$

$$T_d^{\prime }{\dot{E}}_q^{\prime }$$

$$P_j^{el}$$

$$I_{d,j}$$

$${\dot{\delta }}_j$$

$$M_j{\dot{\omega }}_j$$

$$T_j{\dot{E}}_j$$

$${\dot{\omega }}_j={\dot{E}}_j=0,{\dot{\delta }}_j=\Omega \text{mbox}forallj=1,\dots ,N.$$

$$\Omega$$

$$0$$

$${\delta }_j(t)={\delta }_j^{\circ }+{\xi }_j(t),{\omega }_j(t)={\omega }_j^{\circ }+{\nu }_j(t),E_j(t)=E_j^{\circ }+{ϵ}_j(t).$$

$${\dot{\xi }}_j$$

$$M_j{\dot{\nu }}_j$$

$$T_j{\dot{ϵ}}_j$$

$${\Lambda }_{j,\ell }$$

$$A_{j,\ell }$$

$$H_{j,\ell }$$

$$\begin{pmatrix}\dot{\boldsymbol{\xi}}\\ \dot{\boldsymbol{\nu}}\\ \dot{\boldsymbol{\epsilon}}\end{pmatrix}\!=\!\underbrace{\begin{pmatrix}\boldsymbol{0}&\mbox{$\mbox{1}\!\mbox{l}\;$}&\boldsymbol{0}\\ -\boldsymbol{M}^{-1}\boldsymbol{\Lambda}&-\boldsymbol{M}^{-1}\boldsymbol{D}&\boldsymbol{M}^{-1}\boldsymbol{A}^{\top}\\ \boldsymbol{T}^{-1}\boldsymbol{X}\boldsymbol{A}&\boldsymbol{0}&\boldsymbol{T}^{-1}(\boldsymbol{X}\boldsymbol{H}-\mbox{$\mbox{1}\!\mbox{l}\;$})\\ \end{pmatrix}}_{\displaystyle=:\boldsymbol{J}}\!\!\begin{pmatrix}\boldsymbol{\xi}\\ \boldsymbol{\nu}\\ \boldsymbol{\epsilon}\end{pmatrix}.$$

$${\mathcal{D}}_{\perp }=\left\{(\mathbf{\xi },\mathbf{\nu },\mathbf{ϵ}{)}^{\top }\in {\mathbb{R}}^{3N}|1^{\top }\mathbf{\xi }=0\right\}.$$

$${\tilde{\mathcal{D}}}_{\perp }=\left\{(\mathbf{\xi },\mathbf{ϵ}{)}^{\top }\in {\mathbb{R}}^{2N}|1^{\top }\mathbf{\xi }=0\right\},$$

$${\bar{\mathcal{D}}}_{\perp }=\left\{{\mathbf{\xi }}^{\top }\in {\mathbb{R}}^N|1^{\top }\mathbf{\xi }=0\right\}.$$

$${\mu }_1=0,\Re ({\mu }_2)\le \Re ({\mu }_3)\le \cdots \le \Re ({\mu }_{3N}).$$

$$\mathbf{J}\left(\begin{array}{c}\mathbf{\xi }\\ \mathbf{\nu }\\ \mathbf{ϵ}\end{array}\right)=\mu \left(\begin{array}{c}\mathbf{\xi }\\ \mathbf{\nu }\\ \mathbf{ϵ}\end{array}\right).$$

$$\mathbf{\nu }$$

$$-\mathbf{\Lambda }\mathbf{\xi }-\mathbf{D}\mathbf{\nu }+{\mathbf{A}}^{\top }\mathbf{ϵ}$$

$$\displaystyle\boldsymbol{X}\boldsymbol{A}\boldsymbol{\xi}+(\boldsymbol{X}\boldsymbol{H}-\mbox{$\mbox{1}\!\mbox{l}\;$})\boldsymbol{\epsilon}$$

$$-\mathbf{\Lambda }\mathbf{\xi }-\mu \mathbf{D}\mathbf{\xi }+{\mathbf{A}}^{\top }\mathbf{ϵ}$$

$$\mathbf{A}\mathbf{\xi }+(\mathbf{H}-{\mathbf{X}}^{-1})\mathbf{ϵ}$$

$$\begin{gathered} {\mu }^2⟨\mathbf{\xi },\mathbf{M}\mathbf{\xi }⟩+\mu ⟨\mathbf{\xi },\mathbf{D}\mathbf{\xi }⟩-{\mu }^{\ast }⟨\mathbf{ϵ},\mathbf{T}{\mathbf{X}}^{-1}\mathbf{ϵ}⟩ \\ +⟨\mathbf{\xi },\mathbf{\Lambda }\mathbf{\xi }⟩+⟨\mathbf{ϵ},(\mathbf{H}-{\mathbf{X}}^{-1})\mathbf{ϵ}⟩=0, \end{gathered}$$

$$2\Re(\mu)\Im(\mu)\langle\boldsymbol{\xi},\boldsymbol{M}\boldsymbol{\xi}\rangle+\Im(\mu)\big{[}\langle\boldsymbol{\xi},\boldsymbol{D}\boldsymbol{\xi}\rangle+\langle\boldsymbol{\epsilon},\boldsymbol{T}\boldsymbol{X}^{-1}\boldsymbol{\epsilon}\rangle\big{]}=0.$$

$$\Re (\mu )=-\frac{⟨\mathbf{ϵ},\mathbf{T}{\mathbf{X}}^{-1}\mathbf{ϵ}⟩+⟨\mathbf{\xi },\mathbf{D}\mathbf{\xi }⟩}{2⟨\mathbf{\xi },\mathbf{M}\mathbf{\xi }⟩}.$$

$$\mathbf{\Xi }=\left(\begin{array}{cc}-\mathbf{\Lambda }& {\mathbf{A}}^{\top }\\ \mathbf{A}& \mathbf{H}-{\mathbf{X}}^{-1}\end{array}\right)$$

$$(\mathbf{\Xi }+\mathbf{C}(\mu ))\left(\begin{array}{c}\mathbf{\xi }\\ \mathbf{ϵ}\end{array}\right)=0,$$

$$\mathbf{C}(\mu )=\left(\begin{array}{cc}-{\mu }^2\mathbf{M}-\mu \mathbf{D}& 0\\ 0& -\mu {\mathbf{X}}^{-1}\mathbf{T}\end{array}\right).$$

$$\mathbf{\Xi }\left(\begin{array}{c}\mathbf{\xi }\\ \mathbf{ϵ}\end{array}\right)=\beta \left(\begin{array}{c}\mathbf{\xi }\\ \mathbf{ϵ}\end{array}\right)$$

Full text

A reliable supply of electric power requires a stable operation of the electric power grid. Thousands of generators must run in a synchronous state with fixed voltage magnitudes and fixed relative phases. The ongoing transition to a renewable power system challenges the stability as line loads and temporal fluctuations increase. Maintaining a secure supply thus requires a detailed understanding of power system dynamics and stability. Among various models describing the dynamics of synchronous generators, analytic results are available mainly for the simplest second-order model which describes only the dynamics of nodal frequencies and voltage phase angles. In this article we analyze the stability of the third order model including the transient dynamics of voltage magnitudes. Within this model we provide analytical insights into the interplay of voltage and rotor-angle dynamics and characterize possible sources of instability. We provide novel stability criteria and support our studies with the analysis of a network of two coupled nodes, where a full analytic solution for the equilibria is obtained and a bifurcation analysis is performed.

I Introduction

A stable supply of electric power is essential for the economy, industry and our daily life. The ongoing transition to a renewable generation challenges the stability of power grids in several ways R. Sims et al. (2011). Electric power has to be transmitted over large distances, leading to high transmission line loads at peak times Pesch, Allelein, and Hake (2014); Witthaut et al. (2016). Wind and solar power generation fluctuates on various time scales, requiring more flexibility and challenging dynamic stability Anvari et al. (2016); Schmietendorf, Peinke, and Kamps (2017); Schäfer et al. (2017, 2018). Furthermore, the effective inertia of the grid decreases such that power fluctuations have a larger impact on system stability Ulbig, Borsche, and Andersson (2014).

A reliable power supply requires a detailed understanding of power system dynamics and stability. Numerical studies are carried out routinely at different levels of modelling detail (see Refs. Machowski, Bialek, and Bumby, 2008; Sauer and Pai, 1998 for a comparison of different models). These studies provide a concrete stability assessment for one given power grid or components. For instance, the performance of different models has been evaluated for the Western System Coordinating Council System (WSCC) and the New England & New York system in Ref. Weckesser, Jóhannsson, and Østergaard, 2013. Analytic studies into the mathematical structure of the problem have been obtained mainly for second-order models based on the power-swing equations Dörfler and Bullo (2010); Dörfler, Chertkov, and Bullo (2013); Motter et al. (2013); Schäfer et al. (2015). These models describe only the dynamics of nodal frequencies and rotor angles, assuming the voltage magnitudes to be constant in time Machowski, Bialek, and Bumby (2008); Bergen and Hill (1981); Nishikawa and Motter (2015). Voltage stability is usually investigated numerically, see Ref. Simpson-Porco, Dörfler, and Bullo, 2016 and references therein.

The present paper aims at providing analytical insights into the interplay of voltage and rotor-angle dynamics in electric power systems Kundur et al. (2004). We analyze the stability of the third-order model of synchronous generators, including the transient voltage along the $q$-axis, which has been studied so far mainly computationally Machowski, Bialek, and Bumby (2008); Sauer and Pai (1998); Schmietendorf et al. (2014); Ma et al. (2016); Ghahremani, Karrari, and Malik (2008); Dehghani and Nikravesh (2008); Karrari and Malik (2004). We provide an analytical decomposition of the Jacobian into the frequency and voltage subsystems, which gives rise to a novel stability criterion (cf. Proposition 1), and a characterization of possible sources of instability. For the most elementary network of two coupled nodes a full analytic solution of the equilibria is obtained and a bifurcation analysis is performed.

The remainder of the paper is structured as follows: Section II recalls the third-order model. In Section III we perform a local stability analysis via the Jacobian linearization. Section IV discusses routes to instability, while in Section V draws upon the example of a two bus system. The paper ends with conclusions in Section VI.

II The third-order model and its equilibria

The third-order or one-axis model describes the transient dynamics of synchronous machines Machowski, Bialek, and Bumby (2008); Sauer and Pai (1998), in particular the

• •

power angle $\delta(t)$ relative to the grid reference frame,

• •

the angular frequency $\omega(t)=\dot{\delta}$ relative to the grid reference frame and

• •

the transient voltage $E^{\prime}_{q}(t)$ in the $q$-direction of a co-rotating frame of reference.

The third-order model does not cover subtransient effects and it assumes that the transient voltage in the $d$-direction of the co-rotating frame vanishes. The equations of motion for one machine are given by Machowski, Bialek, and Bumby (2008)

[TABLE]

where the dot denotes differentiation with respect to time. Here, the symbol $P^{m}$ denotes the effective mechanical input power of the machine and $E^{f}$ the internal voltage or field flux. $P^{el}$ stands for electrical power out-flow. The parameters $D$ and $M$ denote the damping and the inertia of the mechanical motion and $T^{\prime}_{d}$ the relaxation time of the transient voltage dynamics. The voltage dynamics further depends on the difference of the static ($X_{d}$) and transient ($X^{\prime}_{d}$) reactances along the $d$-axis, where $X_{d}-X_{d}^{\prime}>0$ in general, and the current along $d$-axis $I_{d}$.

In this article we consider an extended grid consisting of several synchronous machines labeled by $j=1,2,\ldots,N$. Neglecting transmission line losses, the electric power exchanged with the grid and the current at the $j$th machine read Schmietendorf et al. (2014)

[TABLE]

where the $E^{\prime}_{q,j}$ and $\delta_{j}$ are the transient voltage and the power angle of the $j$th machine, the parameter $B_{jk}\geq 0$ denotes the susceptance of the transmission line $(j,k)$ and $B_{jj}\leq 0$ denotes the shunt susceptance of the $j$th node.

Using equations (2), the equations of motion assume a particularly simple form Schmietendorf et al. (2014); Ma et al. (2016); Schmietendorf, Peinke, and Kamps (2017). For the sake of notational convenience we drop the prime as well the subscripts $d$ and $q$ in the following and obtain

[TABLE]

Many studies consider a grid consisting of synchronous generators and ohmic loads, which can then be eliminated using a Kron reduction Dörfler and Bullo (2013). The reduced system consists of the generator nodes only and the parameters $B_{jk}$ and $P^{m}_{j}$ represent effective values characterizing the reduced network.

The negligence of line losses is a common simplification in power grid stability assessment, as the ohmic resistance is typically much smaller than the susceptance in high-voltage power transmission grids. This assumption is not valid in distribution grids where resistance and susceptance are comparable. Furthermore, losses are expected to become more important when the transmitted power is large, i.e. at the border of the stability region.

Stationary operation of a power grid corresponds to a state with constant voltages and perfect phase-synchronization, i.e. a point in configuration space where all $E_{j}$, $\omega_{j}$ and $\delta_{j}-\delta_{\ell}$ are constant in time. The latter condition requires that all nodes rotate at the same frequency $\delta_{j}(t)=\Omega t+\delta_{j}^{\circ}$ for all $j=1,\ldots,N$, such that we obtain the conditions

[TABLE]

Strictly speaking, we are searching for a stable limit cycle, but all points along the cycle are physically equivalent. We can focus on any point on the cycle as a representative for the equivalence class and call this an equilibrium in the following. Subsequently, the superscript $\cdot^{\circ}$ is used to denote the values of the rotor phase angle, frequency and voltage in this equilibrium state. Perturbations along the cycle, where we add or subtract a global phase shift $\delta$ from all phases $\delta_{j}$ simultaneously, do not affect phase synchronization and thus can be excluded from the stability analysis. We will make this precise in Definition 1.

For the third-order model (3) an equilibrium state of the power grid is thus given by the nonlinear algebraic equations

[TABLE]

We note that many equilibria, stable and unstable, can exist in networks with sufficiently complex topology Manik, Timme, and Witthaut (2017).

III Linear Stability analysis

It is well-understood in the literature Strogatz (2001); Kundur et al. (2004) that local stability properties of an equilibrium, i.e. stability with respect to small perturbations, can be evaluated by linearizing the equations of motion (3). For linear stability analysis we introduce perturbations $\xi$, $\nu$ and $\epsilon$:

[TABLE]

The main question is then whether the perturbations $\xi$, $\nu$ and $\epsilon$ grow or decay over time. If all perturbations decay (exponentially), the equilibrium is said to be ‘linearly’ stable. In the literature, this property is sometimes also called ‘local (asymptotic/exponential) stability’ or ‘small system stability’, cf. Refs. Strogatz, 2001; Khalil, 2002. We also refer to Ref. Kundur et al., 2004 for an overview of stability notions for power systems. Substituting the ansatz from above into (3), transferring to a frame of reference rotating with the frequency $\Omega$ and keeping only terms linear in $\xi_{j}$, $\nu_{j}$ and $\epsilon_{j}$ yields

[TABLE]

whereby $\boldsymbol{\Lambda},\boldsymbol{A},\boldsymbol{H}\in\mathbb{R}^{N\times N}$ are given by

[TABLE]

Furthermore, we define the diagonal matrices $\boldsymbol{M}$, $\boldsymbol{D}$, $\boldsymbol{X}$ and $\boldsymbol{T}$ (all in $\mathbb{R}^{N\times N}$) with elements $M_{j}$, $D_{j}$, $(X_{j}-X^{\prime}_{j})$ and $T_{j}$ for $j=1,\ldots,N$, respectively. We note that all these elements are strictly positive.

Now we can recast (6) into matrix form, defining the vectors $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{N})^{\top}$, $\boldsymbol{\nu}=(\nu_{1},\ldots,\nu_{N})^{\top}$ and $\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{N})^{\top}$, where the superscript $\top$ denotes the transpose of a matrix or vector. We then obtain

[TABLE]

An equilibrium is dynamically stable if small perturbations are exponentially damped. This is the case if the real part of all relevant eigenvalues of the Jacobian matrix $\boldsymbol{J}$ are strictly smaller than zero Strogatz (2001). If an eigenvalue $\mu$ has a positive real part then the corresponding eigenmode grows exponentially as $e^{\Re(\mu)t}$ and the system is linearly unstable. If the real part equals zero then the question for local stability cannot be decided using the linearization approach and a full nonlinear treatment based on a center-manifold approximation is necessary, see Ref. Strogatz, 2001. Typically, power systems are nonlinearly unstable in this case Manik et al. (2014).

We note that $\boldsymbol{J}$ always has one eigenvalue $\mu_{1}=0$ with the eigenvector $\left(\boldsymbol{\xi}~{}\boldsymbol{\nu}~{}\boldsymbol{\epsilon}\right)^{\top}=\left(\boldsymbol{1}~{}\boldsymbol{0}~{}\boldsymbol{0}\right)^{\top}.$ This corresponds to a global shift of the machines’ phase angles which has no physical significance. We thus exclude this trivial eigenmode from the stability analysis in the following, i.e. we consider only perturbations from the orthogonal complement

[TABLE]

Similarly, we define the projection onto the angle and voltage subspaces (omitting frequency)

[TABLE]

and the angle subspace (omitting frequency and voltage)

[TABLE]

Furthermore, we fix the ordering of all eigenvalues such that

[TABLE]

We then have the following consistent definition of linear stability in all directions transversal to the limit cycle (cf. Ref. Strogatz, 2001).

Definition 1**.**

The equilibrium $(\delta^{\circ}_{j},\omega^{\circ}_{j},E^{\circ}_{j})$ is linearly stable if $\Re(\mu_{n})<0$ for all eigenvalues $n=2,\ldots,3N$ of the Jacobian matrix $\boldsymbol{J}$ defined in (14).

We note that an eigenvalue $\mu_{n}$ can be complex, but then also the complex conjugate $\mu_{n}^{*}$ is an eigenvalue. Complex eigenvalues correspond to oscillatory modes, i.e. the variables oscillate around the equilibrium with a growing or shrinking amplitude. Indeed, the following lemma shows that only damped oscillations are possible.

Lemma 1**.**

If an eigenvalue is complex $\Im(\mu_{n})\neq 0$, then its real part is strictly negative $\Re(\mu_{n})<0$.

Proof.

Recall the eigenvalue problem for the Jacobian

[TABLE]

Decomposition yields

[TABLE]

Substituting (19a) into (19b) and multiplying (19c) with $\boldsymbol{X}^{-1}$ yields

[TABLE]

Multiplying the two last equations with the Hermitian conjugates $\boldsymbol{\xi}^{\dagger}$ and $\boldsymbol{\epsilon}^{\dagger}$ from the left, respectively, and equating $(\boldsymbol{\xi}^{\dagger}\boldsymbol{A}^{\top}\boldsymbol{\epsilon})^{\dagger}=\boldsymbol{\epsilon}^{\dagger}\boldsymbol{A}\boldsymbol{\xi}$, one obtains

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the standard scalar product on $\mathbb{C}^{N}$. All matrices in these expressions are Hermitian such that all scalar products are real. Thus we can easily divide (21) into real and imaginary part obtaining

[TABLE]

This condition can be satisfied in two ways:

• (1)

$\Im(\mu)=0$. In this case we have an non-oscillatory mode and $\boldsymbol{\xi}$ and $\boldsymbol{\epsilon}$ can be chosen real. This case is analyzed in Lemma 2.

• (2)

If $\Im(\mu)\neq 0$ we have a pair of complex conjugate eigenvalues and the corresponding eigenmode describes a damped or amplified oscillation. In this case we can divide equation (22) by $\Im(\mu)$ and solve for $\Re(\mu)$ with the result

[TABLE]

As all matrices $\boldsymbol{T}$, $\boldsymbol{X}$, $\boldsymbol{D}$ and $\boldsymbol{M}$ are diagonal with only positive entries we obtain that $\Im(\mu)\neq 0~{}\Rightarrow~{}\Re(\mu)<0.$ ∎

As all oscillatory modes are damped, we can concentrate on real eigenvalues, $\mu\in\mathbb{R}$, in the following.

Lemma 2**.**

The equilibrium $(\delta^{\circ}_{j},\omega^{\circ}_{j},E^{\circ}_{j})$ is linearly stable if and only if the matrix

[TABLE]

is negative definite on $\tilde{\mathcal{D}}_{\perp}$.

Proof.

By contradiction. We first consider the case that $\boldsymbol{\Xi}$ is negative definite and show that this implies $\mu<0$. We start from (20) for the eigenvalues $\mu$ and eigenvectors of $\boldsymbol{J}$, which is rearranged into matrix form:

[TABLE]

where

[TABLE]

Now let $\boldsymbol{\Xi}$ be negative definite on $\tilde{\mathcal{D}}_{\perp}$. If we assume that $\mu\geq 0$ is an eigenvalue, the matrix $\boldsymbol{C}(\mu)$ is negative semi-definite. Thus the sum of the two matrices $\boldsymbol{\Xi}+\boldsymbol{C}(\mu)$ is also negative definite such that equation (25) cannot be satisfied. Thus an eigenvalue with $\mu\geq 0$ does not exist and the system is linearly stable.

Second, consider the case that the matrix $\boldsymbol{\Xi}$ is not negative definite on $\tilde{\mathcal{D}}_{\perp}$. We then find an eigenvalue $\beta\geq 0$ with corresponding eigenvector $(\boldsymbol{\xi},\boldsymbol{\epsilon})^{\top}$. We can assume that $\boldsymbol{\epsilon}\neq\boldsymbol{0}$ because otherwise

[TABLE]

would imply that $\boldsymbol{\Lambda}\boldsymbol{\xi}=\boldsymbol{0}$ and it would mean that either $\boldsymbol{\xi}=\boldsymbol{0}$ or $\boldsymbol{\xi}\propto\boldsymbol{1}$. The first option is impossible because this would imply $(\boldsymbol{\xi},\boldsymbol{\epsilon})^{\top}=\boldsymbol{0}$, the second is ruled out by the assumption that $(\boldsymbol{\xi},\boldsymbol{\epsilon})^{\top}\in\tilde{\mathcal{D}}_{\perp}$.

Then we evaluate the expression

[TABLE]

whereby the last inequality follows from $\boldsymbol{X},\boldsymbol{T}$ being positive definite. Thus the Jacobian is not negative definite on $\mathcal{D}_{\perp}$ and the equilibrium is not linearly stable. ∎

IV Routes to instability

We can obtain further insights into the stability properties by decomposing the system dynamics into its rotor-angle and voltage parts (cf. Refs. Kundur et al., 2004; Vournas, Sauer, and Pai, 1996). We start with the following proposition, which comes in two versions:

Proposition 1** (Sufficient and necessary stability conditions).**

• I.

The equilibrium $(\delta^{\circ}_{j},\omega^{\circ}_{j},E^{\circ}_{j})$ is linearly stable if and only if (a) the matrix $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ and (b) the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}+\boldsymbol{A}\boldsymbol{\Lambda}^{+}\boldsymbol{A}^{\top}$ is negative definite, where $\boldsymbol{\Lambda}^{+}$ is the Moore-Penrose pseudoinverse of $\boldsymbol{\Lambda}$.

• II.

The equilibrium $(\delta^{\circ}_{j},\omega^{\circ}_{j},E^{\circ}_{j})$ is linearly stable if and only if (a) the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$ is negative definite and (b) the matrix $\boldsymbol{\Lambda}+\boldsymbol{A}^{\top}(\boldsymbol{H}-\boldsymbol{X}^{-1})^{-1}\boldsymbol{A}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$.

Proof.

The results follow from Lemma 2 using the Schur complement Zhang (2006). We demonstrate the proof for Part I. The proof for criterion II is obtained analogously.

First, on $\tilde{\mathcal{D}}_{\perp}$, the matrix $\boldsymbol{\Xi}$ can be decomposed into

[TABLE]

with

[TABLE]

We define

[TABLE]

Lemma 3 given in the Appendix shows that

[TABLE]

According to Lemma 2, an equilibrium is stable if and only if

[TABLE]

By the use of equations (29) and (32) the above condition is equivalent to

[TABLE]

That is, the equilibrium $(\delta^{\circ}_{j},\omega^{\circ}_{j},E^{\circ}_{j})$ is stable if and only if $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ and $\boldsymbol{H}-\boldsymbol{X}^{-1}+\boldsymbol{A}\boldsymbol{\Lambda}^{+}\boldsymbol{A}^{\top}$ is negative definite. ∎

We are especially interested in classifying possible routes to instability in the third-order model, i.e. we want to understand how a stable equilibrium can be lost if the system parameters are varied. Proposition 1 provides the basis for this task as it shows where instabilities emerge. Depending on which of the criteria in the corollary is violated, the instability can be attributed to the angular or the voltage subsystem or both. This analysis will furthermore reveal stabilizing factors of the third-order model.

IV.1 Rotor-angle stability

An equilibrium becomes unstable when the Condition I(a) in Proposition 1 is violated, i.e. when the matrix $\boldsymbol{\Lambda}$ is no longer positive definite on $\bar{\mathcal{D}}_{\perp}$. This corresponds to a pure angular instability which can be seen as follows. If we artificially fix the voltages of all nodes we arrive at a second-order model, commonly referred to as classical model Machowski, Bialek, and Bumby (2008), structure-preserving model Bergen and Hill (1981) or oscillator model Filatrella, Nielsen, and Pedersen (2008); Rohden et al. (2012) in the literature. In this case we have $\boldsymbol{\epsilon}\equiv 0$ and the linearized equations of motion (14) reduce to

[TABLE]

This system is stable if and only if $\boldsymbol{\Lambda}$ is positive definite (cf. Ref. Manik et al., 2014). Thus Condition I(a) includes a condition for the stability of the rotor-angle subsystem alone.

For a system of two machines $j=1,2$ connected by one transmission line this stability criterion is easily understood. The matrix $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ if and only if the phase difference satisfies $\cos(\delta_{2}^{\circ}-\delta_{1}^{\circ})>0$. This is a well-known criterion for stability extensively discussed in the literature Machowski, Bialek, and Bumby (2008); Sauer and Pai (1998); Chen et al. (2016a). Stability is lost if the real power transmitted between the two nodes,

[TABLE]

is strongly increased. Typically, the phase difference $\delta_{1}^{\circ}-\delta_{2}^{\circ}$ grows with $P_{1}$ until it reaches $\pi/2$ where the grid becomes unstable.

An example of such a pure angle instability is shown in Fig. 1. For $X_{j}-X_{j}^{\prime}=0$ the voltages are fixed as $E_{j}^{\circ}=E_{j}^{f}$ such that we can restrict our analysis to the stability of the rotor-angle subsystem. Stable operation is possible as long as $P_{1}<P_{\rm max}=B_{1,2}E_{1}^{f}E_{2}^{f}$ as elaborated above. As $P_{1}\rightarrow P_{\rm max}$, the second eigenvalue of the Laplacian $\boldsymbol{\Lambda}$ approaches zero and the stable equilibrium is lost in a saddle-node bifurcation as studied in detail in Refs. Manik et al., 2014; Ma et al., 2016.

This result can be generalized to larger networks in the following way. We assume that the network is connected, otherwise we can consider every connected component separately. A transmission line $(i,j)$ is called a ‘critical line’ if the angle difference across $(i,j)$ satisfies

[TABLE]

otherwise it is called a non-critical line. We then have the following corollary.

Corollary 1**.**

If there is no critical line in the network, then the matrix $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$.

This result follows from the fact that $\boldsymbol{\Lambda}$ is a Laplacian matrix of a graph with weights $w_{ij}=E_{i}^{\circ}E_{j}^{\circ}B_{i,j}\cos(\delta_{j}^{\circ}-\delta_{i}^{\circ})$. If no critical line exists, $\boldsymbol{\Lambda}$ is a conventional Laplacian, which is extensively studied in the literature Merris (1994); Fiedler (1973); Godsil and Royle (2001). In particular, all eigenvalues are positive except for one zero eigenvalue corresponding to the eigenvector $\boldsymbol{1}$, i.e. a global shift of the voltage phase angles which has no physical significance as discussed above.

The general case of a signed Laplacian with possibly negative weights $w_{ij}$ has been studied in detail only recently Manik et al. (2014); Chen et al. (2016a, b); Song, Hill, and Liu (2015). Typically only few critical lines are present such that it is insightful to analyze how stability depends on the properties of these edges Song, Hill, and Liu (2015). One can draw an analogy to resistor networks, where $1/w_{ij}$ is taken as the direct resistance of each edge in the network. The effective resistance between a pair of nodes is then defined as the voltage drop inducing a unit current between the respective nodes.

If only one critical line is present, its direct resistance is negative. However, the Laplacian matrix remains positive semidefinite as long as the effective resistance of the critical line remains positive Zelazo and Bürger (2014).

This result has been generalized to the case of several critical lines as follows. Let $G_{-}$ be the subgraph containing all vertices but only the critical edges and let $F_{-}$ be a spanning forest of this subgraph. We denote the node-edge incidence matrix of the spanning forest as $\boldsymbol{D}_{-}$. Then the matrix

[TABLE]

is an effective resistance matrix: the diagonal elements are the effective resistances of the respective edges and the non-diagonal elements are mutual effective resistances. One then finds:

Corollary 2**.**

The signed Laplacian $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ if and only if the matrix $\boldsymbol{\Gamma}_{-}$ is positive definiteChen et al. (2016b).

IV.2 Voltage stability

Condition II (a) in Proposition 1 describes the stability of the voltage subsystem alone. To see this, we artificially fix the phase angles and frequencies, i.e. we set $\boldsymbol{\xi}\equiv\boldsymbol{\nu}\equiv 0$ and consider only perturbations of the voltages $\boldsymbol{\epsilon}$. The linearized equations of motion (14) then reduce to

[TABLE]

This system is stable if and only if the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$ is negative definite. If this condition is violated we thus face a pure voltage instability. Rotor angles and frequencies and phase angles will be affected by such an instability, too, but the origin lies in the voltage subsystem.

To obtain a deeper understanding of this condition, we first consider a system of only two machines $j=1,2$ connected by one transmission line. Then we have

[TABLE]

The cosine is positive since otherwise the rotor-angle instability would arise; $X_{j}-X_{j}^{\prime}$ and $B_{j,\ell}$ are positive for $j\neq\ell$ while $B_{j,j}\leq 0$. Applying Silvester’s criterion we find that $\boldsymbol{H}-\boldsymbol{X}^{-1}$ is negative definite if and only if

[TABLE]

For applications it is desirable to find sufficient stability conditions which depend purely on the machine parameters and not on the state variables. Using the bound $\cos^{2}(\delta_{1}^{\circ}-\delta_{2}^{\circ})\leq 1$, a sufficient condition for voltage stability is obtained:

[TABLE]

This condition is satisfied if

[TABLE]

Thus, voltage stability is threatened if the difference between the static and transient reactances $(X_{j}-X_{j}^{\prime})$ becomes too large. An example for this route to instability is shown in Fig. 2.

In more detail, we consider two identical machines running idle, i.e. no power is transmitted $P_{1}=P_{2}=0$. In this case the sufficient conditions (34) for stability are also necessary conditions and the stable equilibrium is lost at the critical value

[TABLE]

When $X-X^{\prime}$ is increased above this critical value, an eigenvalue of the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$ crosses zero. Then also an eigenvalue of the Jacobian $\boldsymbol{J}$ of the full dynamical system crosses zero and the equilibrium becomes unstable, as described by Proposition 1.

This result can be generalized to larger systems comprised of more than two machines. Next, we present a sufficient and a necessary conditions for the stability of the voltage subsystem alone. Again, we will support our proposal for small values of $(X_{j}-X_{j}^{\prime})$ for ensuring voltage stability.

Corollary 3**.**

If, for all nodes $j=1,\ldots,N$,

[TABLE]

then the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$ is negative definite.

Proof.

By applying Gershgorin’s circle theorem Gerschgorin (1931) to the general form of the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$, the following condition for eigenvalues $\lambda_{j}$ is obtained:

[TABLE]

with the radius of the disk $\sum_{\ell\neq j}^{N}B_{j,\ell}\cos{(\delta_{j}^{\circ}-\delta_{l}^{\circ})}$ centered at $B_{j,j}-(X_{j}-X_{j}^{\prime})^{-1}$. The matrix is negative definite if and only if all the eigenvalues lie in the left half of the complex plane, which is guaranteed if

[TABLE]

Using the bound $\cos{(\delta_{j}^{\circ}-\delta_{l}^{\circ})}\leq 1$, a sufficient condition for negative definiteness is obtained as

[TABLE]

This concludes the proof. ∎

Corollary 4**.**

If for any subset of nodes $\mathcal{S}\subset\{1,2,\ldots,N\}$,

[TABLE]

then the matrix $\boldsymbol{H}-\boldsymbol{X}^{-1}$ is not negative definite and the equilibrium is linearly unstable.

Proof.

This result follows from evaluating the expression $\boldsymbol{x}^{\top}(\boldsymbol{H}-\boldsymbol{X}^{-1})\boldsymbol{x}$ for a trial vector $\boldsymbol{x}\in\mathbb{R}^{N}$ with entries $x_{j}=1\;\forall j\in\mathcal{S}$ and $x_{j}=0\;\forall j\notin\mathcal{S}$. ∎

IV.3 Mixed instability

The interplay of voltage and angle dynamics can lead to a third type of instability. A genuine mixed instability is observed if the Conditions I (a) and II (a) in Proposition 1 are satisfied such that no ‘pure’ angle or voltage instability occur, but one of the Conditions I (b) or II (b) is violated. An example of a mixed instability is shown in Fig. 3. The matrices quantifying voltage stability ($\boldsymbol{H}-\boldsymbol{X}^{-1}$) and angle stability ($\boldsymbol{\Lambda}$) remain negative and positive definite, respectively, but still stability is lost in a saddle node bifurcation when the transmitted power is increased.

It should be noted that a mixed instability is not exceptional, but the typical case. In the case of two coupled machines a pure angle instability is observed only if $X_{j}-X_{j}^{\prime}=0$ for $j=1,2$ (cf. Fig. 1). If $X_{j}-X_{j}^{\prime}>0$ an increase of the transmitted power can only lead to a mixed instability.

The emergence of mixed instabilities demand more rigid requirements for stability as discussed above. Based on proposition 1, we derive necessary and sufficient conditions in terms of network connectivity measures.

Definition 2**.**

Given the ordering (18), the second smallest eigenvalue $\lambda_{2}$ of the Laplacian matrix $\boldsymbol{\Lambda}$ is known as the Fiedler value or algebraic connectivity of a networkFiedler (1973). The associated eigenvector $\boldsymbol{v}_{F}$ is the Fiedler vectorFiedler (1975).

The algebraic connectivity $\lambda_{2}$ is a measure of the connectivity of a graph, embodying its topological structure and connectedness. For a conventional graph with positive edge weights the algebraic connectivity is greater than [math] if and only if the graph is connected (cf. 1).

In the limit $(X_{j}-X^{\prime}_{j})\equiv 0$ (no voltage dynamics) a necessary and sufficient condition for stability is that the Laplacian $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ which is equivalent to

[TABLE]

as discussed in detail in section IV.1. We can extend this necessary condition for the algebraic condition for small $(X_{j}-X^{\prime}_{j})$ by means of a Taylor expansion. To leading order, the necessary connectivity always increases with $(X_{j}-X^{\prime}_{j})$.

Corollary 5**.**

A necessary condition for the stability of an equilibrium point is given by

[TABLE]

where $\boldsymbol{v}_{F}$ denotes the Fiedler vector for $(X_{j}-X^{\prime}_{j})\equiv 0$.

Proof.

We denote by $\boldsymbol{v}_{F}$ the normalized Fiedler vector at $(X_{j}-X^{\prime}_{j})\equiv 0$ and by $\boldsymbol{v}^{\prime}_{F}$ the actual normalized Fiedler vector for the particular non-zero value of the $(X_{j}-X^{\prime}_{j})$. Clearly, we have

[TABLE]

Furthermore, we use the expansion

[TABLE]

such that

[TABLE]

Now, condition II.(b) in proposition 1 reads

[TABLE]

Choosing one particular vector $\boldsymbol{y}$ one obtains a necessary condition for stability. Picking $\boldsymbol{y}=\boldsymbol{v}^{\prime}_{F}$ we find

[TABLE]

and expanding the right-hand side to leading order in $(X_{j}-X^{\prime}_{j})$ yields

[TABLE]

The matrix $\boldsymbol{X}$ is diagonal such that we can evaluate the expression further to obtain

[TABLE]

∎

Furthermore, we can derive two sufficient conditions for stability in terms of the algebraic connectivity. These results show that in the limit of high connectivity only the pure voltage dynamics determines the stability of the full system.

Corollary 6**.**

If the algebraic connectivity is positive $\lambda_{2}>0$ and for all nodes $j=1,\ldots,N$ we have

[TABLE]

where $\|\cdot\|_{2}$ is the operator $\ell_{2}$-norm, then an equilibrium point is stable.

Proof.

(1) If $\lambda_{2}>0$ this directly implies that $\boldsymbol{\Lambda}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ and condition I.(a) in proposition 1 is satisfied.

(2) Using Gershgorin’s circle theorem as in the proof of corollary 3 we find that the condition (38) implies that \boldsymbol{X}^{-1}-\boldsymbol{H}-\lambda_{2}^{-1}\|\boldsymbol{A}\|_{2}\|\boldsymbol{A}^{\top}\|_{2}\mbox{\mbox{1}!\mbox{l};} is positive definite on $\bar{\mathcal{D}}_{\perp}$. Noting that $\lambda_{2}^{-1}=\|\boldsymbol{\Lambda}^{+}\|_{2}$, this implies

[TABLE]

Hence, matrix $\boldsymbol{X}^{-1}-\boldsymbol{H}-\boldsymbol{A}\boldsymbol{\Lambda}^{+}\boldsymbol{A}^{\top}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$. Condition I.(b) in proposition 1 is satisfied and the equilibrium is stable. ∎

Corollary 7**.**

If the voltage subsystem is stable, i.e. $\boldsymbol{X}^{-1}-\boldsymbol{H}$ is positive definite, and the algebraic connectivity satisfies

[TABLE]

where $\|\cdot\|_{2}$ is the operator $\ell_{2}$-norm, then an equilibrium point is stable.

Proof.

(1) By assumption, we have $\boldsymbol{X}^{-1}-\boldsymbol{H}$ positive definite such that condition II.(a) is satisfied in proposition 1.

(2) Furthermore, the assumption (40) implies that

[TABLE]

such that the matrix $\boldsymbol{\Lambda}-\boldsymbol{A}^{\top}\boldsymbol{(}\boldsymbol{X}^{-1}-\boldsymbol{H})^{-1}\boldsymbol{A}$ is positive definite on $\bar{\mathcal{D}}_{\perp}$ and condition II.(b) in proposition 1 is also satisfied. Hence the equilibrium is stable. ∎

V Existence of Equilibria and Stability Map in Two Bus Systems

The analysis of the Jacobian reveals the stability of a given equilibrium and different routes to instability. We now investigate when a stable equilibrium exists and derive the global stability map of the grid. We focus on an elementary network of two coupled machines, one with positive effective power $P_{1}>0$ and one with negative effective power $P_{2}=-P_{1}$. All other machine parameters are assumed to be identical, $E_{2}^{f}=E_{1}^{f}>0$, $X_{2}-X^{\prime}_{2}=X_{1}-X^{\prime}_{1}>0$ , $B_{2,2}=B_{1,1}<0$ and $B_{2,1}=B_{1,2}>0$.

This example allows for a fully analytical solution of the equilibrium equations (5). Stability requires that phase difference satisfies $\delta_{1}-\delta_{2}\in(-\pi/2,\pi/2)$. Equation (5b) in (5) can then be solved for the phase difference

[TABLE]

Substracting (5c) for one machine from that for the other machine one can show that the voltages at both machines must be identical in a stable equilibrium,

[TABLE]

Using the result (41), we can now eliminate the phases from (5c). We are left with an equation containing only the state variable $E_{1}^{\circ}$:

[TABLE]

This is a fourth-order polynominal in $E_{1}^{\circ}$, which can be solved analytically leading to rather lengthy expressions. More importantly, the fundamental theorem of algebra tells us that there are exactly four solutions for a given set of system parameters. Discarding solutions where $E_{1,2}^{\circ}$ or $\delta_{1}^{\circ}-\delta_{2}^{\circ}$ are complex or $E_{1,2}^{\circ}<0$, which are physically not meaningful, we obtain a set of equilibria, which can be stable or unstable, though.

Solving equation (42) as a function of the system parameters provides a general stability map and bifurcation set. Here, we focus on the dependence on the machine parameters $P_{1}=-P_{2}$ and $X-X^{\prime}$, while keeping the transmission system parameters fixed. The stability map in Fig. 4 reveals three qualitatively different regimes. In region I the polynomial (42) has two real roots, one corresponding to a stable and one to an unstable equilibrium. Hence a stable operation of the grid is possible.

The phase difference decreases steadily with the increase of $X-X^{\prime}$, whereas the nodal voltages increase. Notably, the maximum transmitted power (the maximum power $P_{1}$ for which a stable solution exists, corresponding to the boundary of regions I and II) firstly decreases with $X-X^{\prime}$ and then increases for larger values of $X-X^{\prime}$ due to 2nd power dependence on the voltages on nodes. Near the border of the two regimes I-III the maximum transmitted power is effectively infinite, as the voltage diverges.

When the system is approaching the “phase border” between region I and II in the parameter space, two equilibria, stable and unstable, merge in an inverse saddle-node bifurcation, and the equilibrium is lost, accompanied with the birth of 2 complex roots of the polynomial (42). The border line can be obtained analytically from the determinant of the fourth-order polynomial (42) as it sign changes from $-$ to $+$ when passing through the line I-II. Therefore in region II all four roots of (42) are complex and no stable equilibrium exits. For $X-X^{\prime}=0$ the saddle-node bifurcation corresponds to a pure rotor-angle instability (cf. Fig. 1), for other values of $X-X^{\prime}$ it corresponds to a mixed instability, cf. Fig. 3.

In region III the polynomial (42) has two complex and two real roots with opposite signs. Only the positive real root corresponds to a physical solution, but it is always unstable, such that no stable operation is possible. By approaching the border I-III from below the coefficient of the fourth power in the polynomial (42) becomes infinitesimally small. Consequently the magnitude of one root goes to infinity on the border, thereby changing its sign from positive in region I to negative in region III (cf. Fig. 2).

In physical terms, the voltage magnitude increases with $X-X^{\prime}$ until it diverges on the border of the regions I-III. On the border one faces a pure voltage instability as illustrated in Fig. 2. No stable operation is possible if the difference of static and transient reactance $X-X^{\prime}$ exceeds a critical value (region III). Notably the divergence of the voltage magnitude implies that the phase difference $\delta_{1}^{\circ}-\delta_{2}^{\circ}\rightarrow 0$ on the I-III-border if the transmitted power is finite.

VI Conclusion and Outlook

This paper has investigated the third-order model of power system dynamics with a focus on the interplay of rotor-angle and voltage stability. We employed a linearization approach and derived necessary and sufficient conditions for local exponential stability and uncovered different routes to instability. For the simplified case of a two-bus system, the positions of the equilibria as well as a global stability map was obtained fully analytically.

Our paper provides several rigorous results which can deepen our understanding into the factors limiting power system stability. In particular, the Proposition 1 provides a decomposition of the Jacobian into the angle and voltage subspace by means of the Schur decomposition and thus allows to rigorously classify possible routes to instability. Stability of the angle subsystem requires that voltage phase angle differences remain small, see Corollaries 1 and 2. Stability of the voltage subsystem is threatened if the difference of the static and transient reactances $X-X^{\prime}$ becomes too large as expressed in corollary 3 and 4. Mixed instabilities can emerge due to the interplay of both subsystems.

Furthermore, our paper provides analytical insights into how the structure of a network determines its dynamics by linking stability properties to measures of connectivity. Future research can deepen the understanding and the applicability of our result. In particular, our approach can be extended towards models of increasing complexity such as the fourth order model. Analytic results can be compared to numerical simulations to test how tight the derived bounds are and to clarify the importance of line losses for the stability problem.

Acknowledgements.

All authors gratefully acknowledge support from the Helmholtz Association via the joint initiative “Energy System 2050 – A Contribution of the Research Field Energy”. KS, LRG, MM, DW also acknowledge support by the Helmholtz Association under grant no. VH-NG-1025 to DW and by the Federal Ministry of Education and Research (BMBF grant nos. 03SF0472). TF acknowledges further support from the Baden-Württemberg Stiftung in the Elite Programme for Postdocs. We thank Benjamin Schäfer, Moritz Thümler, Kathrin Schmietendorf, Oliver Kamps and Robin Delabays for stimulating discussions.

Appendix A Construction of the Schur decomposition

In the appendix we summarize some technical results which are used to proof the main results of the paper.

Lemma 3**.**

Consider the matrix $\boldsymbol{U}$ defined in (30). If $\boldsymbol{x}\in\tilde{\mathcal{D}}_{\perp}$, then also $\boldsymbol{y}=\boldsymbol{U}\boldsymbol{x}\in\tilde{\mathcal{D}}_{\perp}$.

Proof.

Let’s write $\boldsymbol{x}=(\boldsymbol{\xi}\boldsymbol{\epsilon})^{\top}$, then

[TABLE]

\mbox{\mbox{1}!\mbox{l};}\boldsymbol{\xi}=0 because $\boldsymbol{x}\in\tilde{\mathcal{D}}_{\perp}$. \mbox{\mbox{1}!\mbox{l};}\boldsymbol{\Lambda}^{+}\boldsymbol{A^{\top}}\boldsymbol{\epsilon}=0 because $\boldsymbol{\Lambda}^{+}$ maps onto $\bar{\mathcal{D}}_{\perp}$ as $\boldsymbol{\Lambda}$ maps onto $\bar{\mathcal{D}}_{\perp}$. Then (\mbox{\mbox{1}!\mbox{l};}\boldsymbol{\xi}-\boldsymbol{\Lambda}^{+}\boldsymbol{A^{\top}}\boldsymbol{\epsilon})=0 and therefore $\boldsymbol{y}\in\tilde{\mathcal{D}}_{\perp}$. Also $\boldsymbol{y}\neq 0$, as if it is equal to zero, then

[TABLE]

which is realized if and only if

[TABLE]

meaning $\boldsymbol{x}=0$, that contradicts with the condition $\boldsymbol{x}\neq 0$. ∎