## Influence of Load Models on Equilibria, Stability and Algebraic Manifolds of Power System Differential-Algebraic System

Dan Wu and Bin Wang

Abstract -Load models have a great impact on voltage behaviors as well as power system transient dynamics. Extensive work has been done on this topic, proposing appropriate load models and capturing better load behaviors during transient. This paper presents a comprehensive study to investigate the geometric and topological changes induced by different load models for the traditional power system differential-algebraic equations. Specifically, we attempt to reveal the deformation of equilibria, stability regions, and algebraic manifolds during a continuous evolution of load model. Several findings are presented in the paper, some of which countering traditional recognitions and intuitions. A major discovery is that the load model with a large proportion of constant impedance and a small proportion of constant power exhibits much more complex features than the load model with the reversed proportions of impedance and power. The increase of complexity is thoroughly recorded and investigated by the changes of geometric properties and mutations of topological invariants in the sense of equilibria, stability regions, and algebraic manifolds for the DAE system. However, most of the changes seem to occur on unstable components of algebraic manifolds or near the singular boundary surfaces, suggesting a limited variation of dynamical behaviors on the stable component.

Index Terms -Load model, differential-algebraic equations, equilibrium, stability region, algebraic manifold

## I. INTRODUCTION

Load behaviors have a profound impact on power system dynamical performance. Inaccurate load models can introduce significant errors in stability analysis [1], fail to capture important phenomena, e.g. fault-induced delayed voltage recovery [2], [3], and even lead to power systems operated in modes where collapse and separation may occur [4].

Power system engineers and researchers have been continuously improving the accuracy of load models in the past several decades [5]-[8]. These efforts along with advanced parameter estimation techniques significantly reduce the gap between the actual measured load behaviors and the numerically simulated behaviors based on the model. Introducing subtler load models can potentially better approximate the measured load behaviors. With today's computing capabilities, there seems to be no reason to avoid an accurate load model, if possible, in transient stability simulations. The state-of-the-art load model is the composite load model developed by the Western Electricity Coordinating Council, which includes static components, four motors and an elec- tronic component [7]. Further improvements are ongoing in the field [8].

Although multiple load models are available, most of them favor numerical simulations in a qualitative way. When applied to the analytical transient stability analysis such as direct methods, they would face several difficulties in either derivation or analysis, for example, formulating an appropriate transient energy function [9] and handling of differentialalgebraic equations (DAEs) instead of ordinary differential equations (ODEs) [10]. This is the reason that most analytical methods were developed on classical power system model whose loads are represented by constant impedance [11], [12]. Further generalizations are possible but not always easy [9], [13]. It has been reported that static load models with properly identified parameters are adequate for transient stability analysis, since the transient stability is mainly about the real power behavior of the load while the static load models can capture the real power behavior with a fairly acceptable accuracy [14].

An earlier attempt in [15] studied the influence of load models on a power system DAE model using a simple system with a single machine, an infinite bus and one load. To study more complex behaviors, higher dimensional state spaces need to be investigated. Thus, this paper studies 5-Bus, 9Bus, and 14-Bus test cases in the center-of-inertia framework, and presents extensive investigations on equilibria, stability regions and algebraic manifolds of differential-algebraic equations. Several new interesting phenomena are revealed, visualized and analyzed when loads are represented by a combination of constant power and constant impedance. At a fixed loading level, when loads transition gradually from the constant power model to the constant impedance model,

- 1) the number of power flow solutions drastically increases, but not monotonically;
- 2) there is always a single stable equilibrium point (SEP) on the stable component of the algebraic manifold;
- 3) neither a type-1 unstable equilibrium point (UEP) on the stable component of the algebraic manifold should admit low voltage at only one bus, nor a solution with low voltage at only one bus should be type-1 UEP;
- 4) along with the gradual disappearance of singular surface, the stability region expands during the load transition, determined initially by the singularity bound-

ary and eventually by the stable manifolds of type-1 UEP(s).

- 5) the algebraic manifold is enlarged during the load transition, connected to its modulo 2 π replica, and eventually occupies the entire angle subspace;
- 6) the topological invariants of the algebraic manifold mutate, from a sphere to a quotient space of a few tori;
- 7) the change of topology for the algebraic manifold only occurs on the unstable components.

The rest of the paper is organized as follows: Section II introduces the power system differential-algebraic equations for transient stability analysis. Section III, IV and V respectively present numerical investigations of equilibria, stability region and algebraic manifold of the IEEE 9-bus system considering loads gradually transitioning from constant power model to constant impedance model. Section VI draws conclusions and envisions the future work.

## II. POWER SYSTEM DIFFERENTIAL-ALGEBRAIC MODEL FOR TRANSIENT STABILITY ANALYSIS

This section presents the modeling of power systems used for this study which can consider different load models. We first introduce the power system model represented by DAEs using the center-of-inertia (COI) framework. Then, we model the load as a combination of constant power load and constant impedance. With an introduced parameter to linearly combine these two load models, we design a continuous transition from one to the other. To identify multiple equilibrium points for these DAEs, an equivalent power flow problem in the COI framework is formualted. Finally, we show how to determine the type of the identified equilibrium points.

## A. Classical Model in COI Framework

Consider a general N -bus power grid with N g generator buses, including PV and slack buses, and N d PQ buses. Appearently, N = N g + N d . We adopt the classical DAE power injection model [16] for the transient stability analysis throughout this paper. A general form is presented by

$$\dot { x } = f ( x, y ) \quad \quad \quad ( 1 \math$$

$$0 = g ( x, y ) \quad \quad \quad ( 1 b )$$

where x ∈ R 2 N g is the differential state vector; y ∈ R 2 N is the algebraic state vector; f : R 2 N → R 2 N g ; g : R 2 N → R 2 N .

For each generator bus 1 , a pair of differential equations captures the angular dynamics of the machine.

$$\begin{array} {$$

$$\dot { \delta } _ { i } & = \omega _ { i } - \$$

where subscript i = 1 2 , , . . . , N g is the index of generator bus; ω i is the rotor angular velocity; δ i is the rotor angle;

1 The generator bus in this dynamical model is the generator internal bus. The generator terminal bus is, thereby, a PQ bus.

ω s is the constant synchronous speed; H i is the inertia constant; P m,i is the constant mechanical power injected into the generator; P e,i is the electrical power delivered from the generator which is defined shortly below; D i is the damping coefficient.

Each generator bus also induces two algebraic equations.

$$\dots \Phi ^ { \dagger } \Phi ^ { \dagger } \$$

where G i,k is the ( i, k ) 's entry of bus conductance matrix; B i,k is the ( i, k ) 's entry of bus susceptance matrix; V m,i is the constant bus voltage magnitude.

We substitute (3a) into (2a) to eliminate P e,i , reducing the DAE system to R 2 N + N g . One can also substitute (3b) in all other equations to further reduce the DAE system. However, later in the paper we will need these voltage equations to search for other equilibria. Hence, we leave them explicit in our DAE system.

For each PQ bus, only two power balance equations need to be specified.

$$\L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L
\L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \l } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \l}} \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \Sigma } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L
& \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \ L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \
& \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \Sigma \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \Dag \Lambda } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L\L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \SQ } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \\ \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \ \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \| \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \R \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \N } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L. & \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \. & \LL } \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L \L$$

where subscript j = N g +1 , N g +2 , . . . , N is the index of PQ bus; P d,j and Q d,j are respectively the active and reactive power loads.

Therefore, the overall DAE system includes (2a), (2b), (3b), (4a), and (4b).

In this paper, we adopt the COI angle framework . Define 2 the COI angular velocity and angle.

$$\omega _ { C O I } \coloneqq \frac { 1 } { M } \sum _ { \substack { i = 1 \\ N } } ^ { N _ { g } } M _ { i } \omega _ { i } \quad \quad ( 5 a )$$

$$\delta _ { C O I } \coloneqq \frac { 1 } { M } \sum _ { i = 1 } ^ { \dots } M _ { i } \delta _ { i } \quad \quad ( 5 b )$$

where M i = 2 H /ω i s and M = ∑ M i Let's consider

.

$$\hat { \omega } _ { i } \coloneqq \omega _ { i } - \omega _ { C O I } \text{ \quad \quad } \text{(6a)}$$

$$\hat { \delta } _ { i } \coloneqq \delta _ { i } - \delta _ { C O I } \text{ \quad \ \ } \text{(6b)}$$

$$\hat { \delta } _ { j } \coloneqq \delta _ { j } - \delta _ { C O I } \text{ \quad \ \ } \text{(6c)}$$

2 Other angle reference choices include a particular infinite bus or an arbitrary PV bus.

Substitute (6) into (2), (3) and (4) we get

$$\dot { \hat { \omega } } _ { i } = & \dot { \omega } _ { i } - \frac { 1 } { M } \sum _ { n = 1 } ^ { N _ { g } } M _ { n } \dot { \omega } _ { n }$$

$$& \top \, \varkappa _ { j, k } \subset \overline { \mu } ( \upsilon _ { j } - \upsilon _ { k } \downarrow ) & \quad \\ 0 = & Q _ { d, j } + V _ { j } \sum _ { k = 1 } ^ { N } V _ { k } ( G _ { j, k } \sin ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) & \\ & - B _ { j, k } \cos ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) ) & ( \\ ( 7 ) \text{ is the DAE system for which we will investig} \\. &. &. &. &. & \Delta \cdots &.$$

$$\begin{array} { l l } \text{Substitute} ( 6 ) \text{ into} ( 2 ), ( 3 ) \text{ and} ( 4 ) \text{ we get} & & \text{voltag} \\ \hat { \omega } _ { i } = & \dot { \omega } _ { i } - \frac { 1 } { M } \sum _ { n = 1 } ^ { N _ { g } } M _ { n } \dot { \omega } _ { n } & ( 7 a ) & \\ \hat { \delta } _ { i } = & \dot { \omega } _ { i } & ( 7 b ) & \text{$\text{$\text{$the$col}$} \\ V _ { i } = & V _ { m, i } & ( 7 c ) & C. \ E q \\ 0 = & P _ { d, j } + V _ { j } \sum _ { k = 1 } ^ { N } V _ { k } ( G _ { j, k } \cos ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) & \text{$\text{$\text{$to$} }$} \\ & + B _ { j, k } \sin ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) ) & ( 7 d ) & \text{$\text{$implic$ further} \\ 0 = & Q _ { d, j } + V _ { j } \sum _ { k = 1 } ^ { N } V _ { k } ( G _ { j, k } \sin ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) & \text{$\text{$to the$} \\ & - B _ { j, k } \cos ( \hat { \delta } _ { j } - \hat { \delta } _ { k } ) ) & ( 7 e ) & ( 2 b ) \, d \\ \text{$\text{$\text{$clng}$} } \\ \text{$\text{$\text{$qt.} } ( 7 ) \text{ is the DAE system for which we will investigate} & \text{to sol} \\ \text{$\text{$rec transient} $dynamic$ Note that $\mathbb{ }M_{\hat{ }e.}$}. = 0 \text{ canoecting} & \text{angles} \end{array}$$

Eqt. (7) is the DAE system for which we will investigate the transient dynamics. Note that ∑ M ω i ˆ i = 0 , suggesting that (7) still has one degree of degeneracy which comes from (7b). We will deal with this issue shortly below when using an equivalent power flow problem to find the equilibria.

## B. Load Modeling

Equation. (4) presents the power balance relation at each PQ bus 3 . A typical load model is the constant power model, assuming that power injections P d,j and Q d,j at each PQ bus are constants. This model is very useful for static voltage stability analysis and induces the traditional power flow problem. During the transient, however, voltages at PQ buses fluctuate, which can alter the power consumption from their designated values. To better capture the change of power consumption, the 'ZIP' model is formulated in the following way.

$$P _ { d, j } = P _ { 0, j } + I _ { p, j } V _ { j } + G _ { d, j } V _ { j } ^ { 2 } \text{ \quad \ \ } ( 8 \text{a} )$$

$$Q _ { d, j } = Q _ { 0, j } + I _ { d, j } V _ { j } + B _ { d, j } \overset { \cdot } { V } _ { j } ^ { 2 } \text{ \quad \ \ } ( 8 b )$$

where P 0 ,j and Q 0 ,j are the constant active and reactive power; I p,j and I d,j are the constant active and reactive current; G d,j is the load conductance; B d,j is the load susceptance.

In this paper, we ignore the constant current part in (8), i.e. I p,j = I q,j = 0 , to acquire a uni-directional change of load model from the constant power to the constant impedance. Consider

$$P _ { d, j } + j Q _ { d, j } = \alpha ( P _ { 0, j } + j Q _ { 0, j } ) + ( 1 - \alpha ) V _ { j } ^ { 2 } / Z _ { d, j } \ \ ( 9 ) \ \stackrel { \mathfrak { s } \, \text{cuda} } { D \ D e t. }$$

where α ∈ [0 , 1] ; Z d,j is the constant impedance.

We can track the change of transient dynamics by continuously changing α from 1 to 0 , which gradually convert the constant power model to the constant impedance model. P 0 ,j + jQ 0 ,j is fixed to be the designated load power consumption. When solving the traditional power flow high

3 We choose the flow convention that injecting power is the positive direction.

voltage solution V 0 ,j associated with P 0 ,j + jQ 0 ,j , we obtain the corresponding impedance by

$$Z _ { d, j } = V _ { 0, j } ^ { 2 } / ( P _ { 0, j } + j Q _ { 0, j } ) \quad \quad ( 1 0 )$$

This setting ensures that the high voltage solution is unchanged during the change of load model.

## C. Equivalent Power Flow Problem for COI Framework

To evaluate the transient dynamics of the system (5), we need to obtain its equilibrium points. Setting ˙ ˆ δ i to zero in (7b) implies that the relative angular velocity ˆ ω i vanishes, which further suggests that the true angular velocity ω i converges to the COI angular velocity ω COI that does not necessarily comply with the synchronous speed ω s . In this case, ˙ δ i in (2b) does not vanish. Hence, all the generator angles keep changing, while their relative differences stabilize. Therefore, to solve the relative angle differences we subtract all COI angles from the first COI angle. Specifically, we define

$$\tilde { \omega } _ { i } \coloneqq \hat { \omega } _ { i } - \hat { \omega } _ { 1 } = \omega _ { i } - \omega _ { 1 } \text{ \quad \ \ } ( 1 1 a )$$

$$\tilde { \delta } _ { i } \coloneqq \hat { \delta } _ { i } - \hat { \delta } _ { 1 } = \delta _ { i } - \delta _ { 1 } \text{ \quad \quad \ \ } ( 1 1 b )$$

$$\tilde { \delta } _ { j } \coloneqq \hat { \delta } _ { j } - \hat { \delta } _ { 1 } = \delta _ { j } - \delta _ { 1 } \text{ \quad \quad \ \ } ( 1 1 c )$$

Assuming that D i = D j for every pair i and j , the algebraic equations we are going to solve with respect to ˜ δ i are

$$\frac { \ e } { \ w } ^ { \ e } \quad 0 = & \frac { \omega _ { s } } { 2 H _ { i } } ( P _ { m, i } - P _ { e, i } ) - \frac { \omega _ { s } } { 2 H _ { 1 } } ( P _ { m, 1 } - P _ { e, 1 } ) \quad \ ( 1 2 a ) \\ \ s s \quad \ \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \ \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w \quad \ w$$

$$\overset { \varepsilon } { \text{ir} } \quad 0 = & V _ { m, i } - V _ { i } \quad \text{(12b)}$$

$$\begin{smallmatrix} 0 = & V _ { m, i } ^ { \dagger } - V _ { i } & & ( 1 2 b ) \\ \mathfrak { r } \\ 0 = & P _ { d, j } + V _ { j } \sum _ { k = 1 } ^ { N } V _ { k } ( G _ { j, k } \cos ( \tilde { \delta } _ { j } - \tilde { \delta } _ { k } ) \\ & + B _ { j, k } \sin ( \tilde { \delta } _ { j } - \tilde { \delta } _ { k } ) ) & & ( 1 2 c ) \\ \mathfrak { i } ) \\ 0 = & Q _ { d, j } + V _ { j } \sum _ { k = 1 } ^ { N } V _ { k } ( G _ { j, k } \sin ( \tilde { \delta } _ { j } - \tilde { \delta } _ { k } ) \\ \mathfrak { e } \\ & - B _ { j, k } \cos ( \tilde { \delta } _ { j } - \tilde { \delta } _ { k } ) ) & & ( 1 2 d ) \\ \mathfrak { e } \\ d \\ \quad \text{Equation. (12) is the equivalent power flow problem that } \\ \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \delta \end{smallmatrix}$$

Equation. (12) is the equivalent power flow problem that provides the equilibrium points to the DAE system (7). It is similar to the traditional power flow problem except that every active power balance equation at generator bus subtracts from the active power balance equation at a particular generator bus.

## D. Determining Stability and Type of Equilibrium

Suppose ( x , y 0 0 ) is an equilibrium point to the DAE system (1). Let's consider the Jacobian matrix J x for the dynamic states x . Suppose ∂g ∂y | ( x ,y 0 0 ) is nonsingular, then

$$J _ { x } \coloneqq \left [ \frac { \partial f } { \partial x } - \frac { \partial f } { \partial y } \Big ( \frac { \partial g } { \partial y } \Big ) ^ { - 1 } \frac { \partial g } { \partial x } \right ] _ { ( x _ { 0 }, y _ { 0 } ) } \quad ( 1 3 )$$

Gathering the eigenvalues of J x in Λ = { λ , . . . , λ 1 2 N g } . The point ( x , y 0 0 ) is said to be typek if there exists exactly k entries of Λ whose real parts are positive. If J x is hyperbolic and type0 , then we say ( x , y 0 0 ) is an SEP.

Note that (7) has one degree of degeneracy, the Jacobian matrix J x of dynamic states will always have a zero eigenvalue. Ignoring this zero eigenvalue, the type and stability of an equilibrium point can still be determined by the rest of the eigenvalues.

## E. Model Modification and Computational Procedure

As discussed in the above subsection, the types of equilibria are evaluated from the differential part of the DAE system (7). Therefore, the static power flow problem is not sufficient to determine the equilibria and their types. In this paper, we add generator classical model to each PV bus of the power flow problem, set the generator internal voltage bus as the new PV bus, and revise the generator terminal bus to PQ bus. The modified system is then a DAE system with more nodes and components. For example, the 9-bus system turns out to be a 12-bus system after modification since it includes three generator internal buses. We still call this system the 9-bus system since readers would be more familiar with the standard power flow system. The tested systems presented in this paper are all selected from the standard power flow systems, namely, 5-bus case [17], [18], 9-bus case, and 14bus case [19], with modifications on the dynamical parts. The dynamical parameters are selected from the Appendices of [20], and are included in Appendix I of this paper.

The computational procedure is given below.

- 1) Choose a static power flow problem. Solve the high voltage solution by a standard solver.
- 2) Add generator internal bus to each generator bus of the static power flow problem. Compute the internal bus voltage based on the previously solved high voltage solution. Revise the generator terminal bus to PQ bus.
- 3) Compute the equivalent load impedance at the high voltage solution. Choose the load model with given proportions of constant power and constant impedance.
- 4) Based on the new DAE model, solve the corresponding equivalent power flow problem (12) for many equilibrium points.
- 5) Evaluate the types and properties of these equilibrium points by (13).

## III. EQUILIBRIA WITH DIFFERENT LOAD MODELS

This section applies an efficient method [21], [22] to solve (12) for multiple equilibra of the DAE model (7). We gradually change the load model from constant power to constant impedance and observe the change of equilibra.

## A. Number of Equilibria

The numbers of equilibria at different load models are summarized in Fig. 1. Note that we apply (9) as our load model for different proportions of constant power and impedance associated with α . One observation is that a small proportion of power with a large proportion of impedance load induces much more equilibria than a large proportion

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-7084d603-bb54-463b-8720-a6dca4db8059.jpg)

(a)

5-Bus Case

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-1407e970-8c6a-4813-9b12-0e63636cae18.jpg)

(b)

9-Bus Case

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-c798ef18-2ace-4fcc-a661-8d8b298a112d.jpg)

Fig. 1: Numbers of Equilibria with Different Load Models

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-2073f8c8-2941-4620-b5bf-d2251635fbf3.jpg)

(a)

Number of Equilibria Change

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-a611dab8-a3d1-48e1-bac4-8b10741b2a66.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-5783b4e6-0685-4a77-a4e0-8a0cd1bee476.jpg)

3

(b) Solution Pair Disappearance

(c) Solution Pair Appearance

Fig. 2: Non-monotone Change of Solutions in 14-Bus Case

of power with a small proportion of impedance load. This may be counter-intuitive because the constant impedance load is considered simpler than the constant power load. But as the constant power is gradually replaced by the constant impedance, it is equivalent to reducing the node power injection while increasing the shunt impedance. When the node power injection declines, it actually can intersect with more PV (QV) curves (surfaces), thus resulting in more equilibria.

The second observation is that, as the load model approaches the constant impedance model, the number of equilibria increases much faster. This phenomenon may also be explained by the reason that many PV (QV) curves (surfaces) occur at a low node power injection level. It suggests that a small change of load model can result in many equilibria appearing or disappearing at a small constant power injection level. For example, in Fig. 1(c), when the constant impedance increases from 80% to 90% , the number of equilibria doubles. It may raise a reasonable concern about the validity of transient stability analysis based on some heuristically or conceptually assigned load models. However, this concern can be partially resolved by the observations in the next subsection when we further investigate the types of equilibria on the stable component of algebraic manifold.

Fig. 1(b) and Fig. 2(a) show that the number of equilibria does not necessarily increase monotonically as the proportion of constant power decreases. In Fig. 1(b), the number of equilibira declines around the pure constant impedance load model. In Fig. 2(a), the number of equilibia temporarily reduces at 18 81% . impedance proportion, and returns as the impedance proportion increases to 23% . Fig. 2(b) depicts the process of this temporary equilibria reduction. As the impedance proportion increases from 18% to 18 80612% . , a pair of equilibria collides on the real plane, and becomes a complex-valued solution pair. This collision happens between a type-1 equilibrium and a type-2 equilibrium. Fig. 2(c), on the other hand, shows the acquisition of another equilibria pair which is not the same colliding equilibria pair in Fig. 2(b). As the impedance proportion increases from 22 98% . to 23 01% . , a pair of complex-valued solutions collides and becomes a real-valued equilibria pair from which one is type-0 and the other is type-1. At a particular proportion value between 23 01% . and 23 05% . , the type-1 equilibrium turns into type-0.

## B. Types of Equilibria

After enumerating equilibrium points for different load models, we evaluate their types displaying in Fig. 3. A DAE system (7) can be regarded as a confined dynamical system whose dynamical flows are restricted on the algebraic manifold defined by the algebraic equations of the DAE system. The type of an equilibrium point is determined by its differential part in (13). Therefore, multiple SEPs can exist on different components of the algebraic manifold. However, even some equilibria are classified as stable (type-0) in the analysis, the physical system may not be able to work at them. The algebraic part of a DAE system is usually derived from reducing fast dynamics, or other simplifications. Hence, an SEP in the DAE model may not be 'stable' in the physical sense. At this point, equilibria on the stable component of the algebraic manifold are of most interest because they are usually consistent with the types in the real situation.

One observation from Fig. 3 is that every tested case in each load model only exhibits one SEP (which is the high-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-ffda20b6-7bf3-4f13-b14d-e0f3b581b6a9.jpg)

(a) 5-Bus Case

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-93736b30-155a-48ce-87b4-ca2caf5593e0.jpg)

(b) 9-Bus Case

(c) 14-Bus Case

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-52efb20e-2f4d-4b7e-836f-cf810bd3cdb6.jpg)

(d) 14-Bus Case on Stable Manifold

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-7a1b84dc-3f3f-45eb-9228-19c15f8bf947.jpg)

Fig. 3: Equilibria Types for Different Load Models

voltage solution) on the stable component of the algebraic manifold. Ref. [23] discussed the uniqueness of load flow solution for radial distribution network, however, it is still an open question whether the SEP is unique or not on the stable component of algebraic manifold for a general power grid DAE system. Moreover, although the total number of equilibria can change drastically as shown in Fig. 1, the number of equilibria on the stable component of algebraic manifold remains small for all tested cases in this paper. This observation may partially relieve the concern raised in the previous subsection: the change of load model can substantially alter the number of equilibria, but only very limited equilibria are modified on the stable component of algebraic manifold. More investigations are needed.

A more interesting finding is that the equilibria on the stable component of algebraic manifold become complicated with different types when the proportion of constant power reduces to a particular value, which is usually less than 20% . Fig. 3(b) and Fig. 3(d) present the increase of complexity for the equilibria. This phenomenon is somehow counterintuitive because the constant power load is thought to be more complicated than the constant impedance load. Our analysis shows the opposite, at least in the sense that a small proportion of constant power with a large proportion of constant impedance will result in more complex dynamical behaviors.

TABLE I: Selected Equilibria from 9-Bus Case

| Bus           | V 1    | V 2    | V 3    | V 4    | V 5    |
|---------------|--------|--------|--------|--------|--------|
| Equilibrium 1 | 0.3519 | 0.7600 | 0.7565 | 0.0696 | 0.1584 |
| Equilibrium 2 | 0.5889 | 0.8282 | 0.7490 | 0.4042 | 0.0037 |
| Bus           | V 6    | V 7    | V 8    | V 9    | Type   |
| Equilibrium 1 | 0.6060 | 0.5874 | 0.5939 | 0.1488 | 1      |
| Equilibrium 2 | 0.5948 | 0.6505 | 0.7081 | 0.4883 | 0      |

Another interesting discovery may also counter people's intuition: neither a type-1 equilibrium point (on the stable component of the algebraic manifold) should admit only one low bus voltage, nor a solution with only one low bus voltage should be type-1. For example, in Table I, we present two selected equilibria from the 9-bus case at 97% proportion of impedance load. The first equilibrium is type-1 on the stable component of algebraic manifold. The voltage magnitudes on bus 4, 5 and 9 are all below 0 2 . p.u. While the second equilibrium is type-0 on the unstable component of algebraic manifold. It only has one voltage magnitude at bus 5 below 0 2 . p.u. Other tested cases also exhibit a similar phenomenon. This observation suggests that an initial guess with only one low voltage magnitude cannot necessarily lead the Newton's method 4 to find a type-1 equilibrium point. It also suggests that an equilibrium with only one low voltage magnitude can also be another type.

## IV. STABILITY REGIONS WITH DIFFERENT LOAD MODELS

Section III shows that power system DAEs can have multiple SEPs on different components of the algebraic manifold. This section will investigate the load model influence on the stability region (or the stability boundary which is the boundary surface of region of attraction) with respect to the high-voltage SEP. Unlike the power system ODE dynamical models whose stability boundaries are composed of the stable manifolds of certain type-1 UEPs, the unstable manifolds of which converge to the SEP, the stability boundaries of power system DAE dynamical models consist of two components: the first component is the same as the ODE model, while the second one is a set of points whose trajectories reach singular surfaces [15], [24]-[26].

4 Any method that requires a starting point may fail to identify a type-1 equilibrium point from the mentioned initial guess.

Fig. 4 shows the cross-sections of the stability regions on the zero-speed plane 5 for the 9-bus case with different load models. These intersections can be regarded as particular slices of the corresponding stability regions of the power system DAE models in the δ -ω space. For simplicity, we refer this intersection to 'stability region' in the following discussions. In Fig. 4, green solid circles are equilibria on the stable component, red curves represent stability boundaries identified by Algorithm 1, and blue stars represent the projection of singular surfaces on the relative angle plane.

As shown in Fig. 4(a)-4(c), when the proportion of impedance load increases from 0% to 80% , the stability region of the 9-bus case expands gradually. During this process, there is only one equilibrium, which is stable, and the stability boundary is exactly the same as the singularity boundary.

When the proportion of impedance load reaches 90% , the singularity boundaries break up into multiple segments in the plot, and a type-1 UEP appears around ( δ 21 , δ 31 ) = (2 4 . , 0 2) . , whose stable manifolds contribute to a small portion of the stability boundaries around this UEP, as shown in Fig. 4(d) and 4(e). Note that the stability boundary near ( -3 , -1) does not overlap with any singularity boundaries because the stability region shown in the plot is just an angle slice that does not necessarily comply with the angle slice of the singularity boundary. The dash dot line in Fig. 4(e) represents an unstable trajectory initialized at ( δ 21 , δ 31 ) = ( -2 08 . , 0 19) . , i.e. a point outside stability region but inside the singularity region. As a comparison, the black trajectory initialized at ( δ 21 , δ 31 ) = ( -2 48 . , 0 07) . , i.e. a point in the stability region but close to the boundary, first crosses the SEP and travels near the type-1 UEP. But it converges to the SEP eventually.

When the load is represented by a pure impedance, the singularity boundary completely disappears and the stability boundary is fully determined by stable manifolds of the two type-1 UEPs, as shown in Fig. 4(f).

Another interesting finding is that except for the case with 100% impedance load model, the majority of the stability boundary is determined by the singular surface. When the proportion of constant power reduces from 100% to less than 5% , the singular surface expands in the relative angle subspace and disappears eventually. In the meantime, the expanded area tends to include more UEPs, which contribute to a portion of the stability boundary.

## V. ALGEBRAIC MANIFOLDS WITH DIFFERENT LOAD MODELS

Our analysis started with equilibria in Section III, and extended to stability regions on the stable component of algebraic manifold in Section IV. Now in this section we directly investigate the algebraic manifold to acquire a global view of how the load model influences the DAE system. The

(rad)

31

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-08a43f73-4213-483d-a888-fdc53dec00b6.jpg)

(a)

0%

Impedance

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-8bd64583-f81e-4ebc-9e56-aa2d3795e6e7.jpg)

(b)

40%

Impedance

(rad)

31

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-82ce0368-3df3-4fa0-9b43-b5b040daffc8.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-75ecf1da-e76f-4c00-896e-8e44791a970c.jpg)

Fig. 4: 9-Bus Case Stability Boundary

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-7a697cda-a366-41df-917c-5dd7202d2920.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-3b87ab49-4e37-401f-aee8-024ead718ed4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-8ab7937e-b03d-4f07-8591-ec6ca391eaa1.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-a7afe1ec-1bf9-439a-a26e-3c0f06ca21d5.jpg)

(c)

80%

Impedance

(rad)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-7c5f58d4-bbaf-4ef0-a829-fb9e2ad78821.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-9c49027b-ad7b-4a38-8dc9-e2b947e31e4b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-e475cd2f-dfb2-4b18-80a9-ed1757d88a98.jpg)

Fig. 5: 9-Bus Case Algebraic Manifolds

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-e94e4dbc-61c6-44a2-a1a9-3d5a91c43c9d.jpg)

Algorithm 1 Numerical identification of the stability boundary

- 1: In the δ 21 -δ 31 plane, set up M unit vectors with angles respectively taking 0 , 360 /M , ..., 360( M -1) /M , say n 1 , ..., n M . M = 180 is used in this paper.
- 2: Set r to be a small step, e.g. 0.1 used in this paper. Along the direction determined by each of the unit vectors, say n j , let ( δ 210 , δ 310 ) = r n j and conduct steps 3-6 below.
- 3: Use ( ˜ δ 1 , ω ˜ 1 , δ ˜ 2 , ω ˜ 2 , δ ˜ 3 , ω ˜ ) = (0 3 , 0 , δ 210 , 0 , δ 310 , 0) as the initial state to numerically solve DAEs in (7) over a period of time (5 seconds in this paper) for δ 21 ( ) t and δ 31 ( ) t .
- 4: If r &lt; /epsilon1 ( /epsilon1 = 0 001 . in this paper) record ( δ 210 , δ 310 ) as an estimate of the stability boundary in direction n j and then go to step 2 for next direction.
- 5: If max ( δ 21 ) -min ( δ 21 ) &gt; 2 π or max ( δ 31 ) -min ( δ 31 ) &gt; 2 π , let r = r/ 2 , ( δ 210 , δ 310 ) = ( δ 210 , δ 310 ) -r n j and go to step 3.
- 6: Let ( δ 210 , δ 310 ) = ( δ 210 , δ 310 ) + r n j and go to step 3.

algebraic manifold is defined by the algebraic set of (7c), (7d), and (7e). We only focus on the 9-bus case because it is complicated enough to demonstrate interesting phenomena, while not computationally extensive.

Fig. 5 depicts algebraic manifolds under different load models. Specifically, the X -axis represents the generator angle ˜ δ 2 which equals to the angle difference δ 21 := δ 2 -δ 1 as described by (11b). The Y -axis similarly represents ˜ δ 3 = δ 31 . The Z -axis is chosen to be the voltage magnitude on bus9 . The algebraic manifolds are shown in pink surfaces on which the blue curves represent the singularity boundary surfaces. The yellow diamonds are the equilibria that we analyzed in Section III.

From Fig. 5(a) to Fig. 5(f) one can observe that the algebraic manifold enlarges as the proportion of impedance load increases. This is consistent with the finding that constant power model induces smaller stability region than the constant impedance model [27]. Since the X and Y axes represent angles, they should repeat themselves by any 2 nπ . When the impedance model dominates the load model, the algebraic manifold connects to its 2 nπ replica, as shown in Fig. 5(f). It suggests that the DAE system dynamics may stabilize after several angular cycles of 2 π .

A new finding is that a large proportion of impedance in the load really complicates the geometry of the algebraic manifold, and substantially alters the topology of the manifold as well. For example, with constant power load model, the algebraic manifold depicted in Fig. 5(a) is homeomorphic to a sphere. When the impedance load increases to 20% in Fig. 5(b), a small depression area appears. It induces negative curvatures and creates its own singularity boundary which is shown as the small blue circle in the plot. When the impedance proportion increases further to 40% , the small depression area becomes a concave bubble inside the outer

(a) 40% Impedance

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-0d2e3395-537f-4b0d-aac3-6497eca634c7.jpg)

Fig. 6: 9-Bus Case Algebraic Manifolds: Inner Bubbles

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01111v1-figure-192039a3-c5c9-4543-aa63-cde9434039a2.jpg)

manifold and introduces more equilibria on it, shown in Fig. 5(c) (global view) and Fig. 6(a) (local view). As we increase the load impedance to 60% , two inner bubbles reach the other side of the outer surface, punctuating the manifold to create two tunnels shown in Fig. 5(d) (global view) and Fig. 6(b) (local view). Now the manifold is no longer homeomorphic to a sphere, but homeomorphic to two tori glued together. Increasing the impedance proportion further enlarges the tunnels and creates more complicated structures as shown in Fig. 5(e) and 5(f).

In this particular example, however, both the geometric and the topological changes occur on the unstable components of algebraic manifold. Whether the stable component can exhibit a similar change requires further investigations.

## VI. CONCLUSION

Load models are recognized being influential on voltage behaviors and transient dynamics. Numerous works have been done on this topic to build appropriate load models and to better capture static and dynamical behaviors.

In this paper, we specifically focused on the geometric and topological changes of the solution sets for a power system DAE model with different load models. A few findings were recorded and analyzed when the load model continuously transitions from the constant power model to the constant impedance model. One major counter-intuitive discovery is that a load model with a large proportion of constant impedance and a small proportion of constant power introduces much more complex geometric and topological structures than a load model with a small proportion of constant impedance and a large proportion of constant power. Specifically, the number of equilibria increases dramatically when the constant power is largely replaced by the constant impedance. However, this increase of equilibrium quantity is not necessarily monotone. We also emphasized that the appearance of type-1 UEP is not necessarily associated to a single node experiencing ultra-low voltage magnitude. The stability region, as expected, starts from a small region determined by the singular surface, expands when the proportion of constant impedance increases, and is finally constrained by the stable manifolds of type-1 UEPs. In this case, the algebraic manifold is connected to all its modulo 2 π replica. The fundamental group of the algebraic manifold is, however, more complicated when the load is dominated by the impedance.

These findings certainly provide deep insights to the fundamental geometric and topological structures of power system dynamics associated with different load models. Future investigations will include, not limited to, how to identify changes of equilibria only on the stable component of algebraic manifold; what changes of other load models can exert on power system dynamical models; how to identify the appropriate distance from the SEP to the nearest stability boundary point in the power system DAE model.

## ACKNOWLEDGEMENT

The authors gratefully acknowledge Prof. Bernard Lesieutre, Prof. Chris Demarco, Prof. Franz-Eric Wolter, and Prof. Marija Ilic for many helpful discussions. We also thank Prof. Eytan Modiano and Prof. Le Xie for their great supports.

## REFERENCES

- [1] W.-J. Lee, M.-S. Chen, and L. B. Williams, 'Load model for stability studies,' IEEE Transactions on Industry Applications , vol. IA-23, no. 1, pp. 159-165, 1987.
- [2] D. Kosterev, A. Meklin, J. Undrill, B. Lesieutre, W. Price, D. Chassin, R. Bravo, and S. Yang, 'Load modeling in power system studies: Wecc progress update,' in 2008 IEEE Power and Energy Society General Meeting -Conversion and Delivery of Electrical Energy in the 21st Century , pp. 1-8, July 2008.
- [3] D. Kosterev and J. Undrill, 'Load modeling in wecc.' DOE-NERC FIDVR Workshop, 2009.
- [4] CIGRE Task Force 38.02.05, 'Load modeling and dynamics.' Electra, 1990.
- [5] IEEE Task Force on Load Representation for Dynamic Performance, 'Load representation for dynamic performance analysis,' IEEE Trans. on Power Syst. , vol. 8, pp. 472-482, May 1993.
- [6] IEEE Task Force on Load Representation for Dynamic Performance, 'Standard load models for power flow and dynamic performance simulation,' IEEE Trans. on Power Syst. , vol. 10, pp. 1302-1313, Aug 1995.
- [7] A. Arif, Z. Wang, J. Wang, B. Mather, H. Bashualdo, and D. Zhao, 'Load modelinga review,' IEEE Transactions on Smart Grid , vol. 9, pp. 5986-5999, Nov 2018.
- [8] Q. Huang, R. Huang, B. J. Palmer, Y. Liu, S. Jin, R. Diao, Y. Chen, and Y. Zhang, 'A reference implementation of wecc composite load model in matlab and gridpack.' arXiv preprint 1708.00939.
- [9] Y. Min and L. Chen, 'A transient energy function for power systems including the induction motor model,' Science in China Series E: Technological Sciences , vol. 50, no. 5, pp. 575-584, 2007.
- [10] K. R. Padiyar, Structure Preserving Energy Functions in Power Systems: Theory and Applications . CRC Press, 2017.
- [11] M. A. Pai, Power System Stability: Analysis by the Direct Method of Lyapunov . North-Holland Publishing Company, 1981.
- [12] H.-d. Chiang, Direct Methods for Stability Analysis of Electric Power Systems: Theoretical Foundation, BCU Methodologies, and Applications . John Wiley &amp; Sons, 2011.
- [13] I. A. Hiskens and R. J. Davy, 'Lyapunov function analysis of power systems with dynamic loads,' in Proceedings of 35th IEEE Conference on Decision and Control , vol. 4, pp. 3870-3875 vol.4, Dec 1996.
- [14] Y. Li, H. . Chiang, B. . Choi, Y. . Chen, D. . Huang, and M. G. Lauby, 'Representative static load models for transient stability analysis: development and examination,' IET Generation, Transmission Distribution , vol. 1, pp. 422-431, May 2007.
- [15] H. Wu, R. Guo, Z. Han, and D. Gan, 'Geometrical structure of constraint manifold in power system differential-algebraic model,' in 2006 International Conference on Power System Technology , pp. 1-7, Oct 2006.
- [16] P. W. Sauer, M. A. Pai, and J. H. Chow, Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox . John Wiley &amp; Sons, 2017.
- [17] F. Salam, L. Ni, S. Guo, and X. Sun, 'Parallel processing for the load flow of power systems: the approach and applications,' in Decision and Control, 1989., Proceedings of the 28th IEEE Conference on , pp. 2173-2178, IEEE, 1989.
- [18] F. Salam, L. Ni, X. Sun, and S. Guo, 'Parallel processing for the steady state solutions of large-scale non-linear models of power systems,' in Circuits and Systems, 1989., IEEE International Symposium on , pp. 1851-1854, IEEE, 1989.
- [19] R. D. Christie, 'Power systems test case archive.' University of Washington, Seattle, 1999.
- [20] P. M. Anderson and A. A. Fouad, Power system control and stability . John Wiley &amp; Sons, 2008.
- [21] D. Wu and B. Wang, 'Holomorphic embedding based continuation method for identifying multiple power flow solutions,' IEEE Access , 2019.
- [22] B. Lesieutre and D. Wu, 'An efficient method to locate all the load flow solutions-revisited,' in 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton) , pp. 381-388, IEEE, 2015.
- [23] H.-D. Chiang and M. E. Baran, 'On the existence and uniqueness of load flow solution for radial distribution power networks,' IEEE Transactions on Circuits and Systems , vol. 37, no. 3, pp. 410-416, 1990.
- [24] V. Venkatasubramanian, H. Sch¨ attler, and J. Zaborszky, 'Stability regions for differential-algebraic systems,' in Systems, Models and Feedback: Theory and Applications , pp. 385-402, Springer, 1992.
- [25] V. Venkatasubramanian, H. Schattler, and J. Zaborsky, 'Dynamics of large constrained nonlinear systems-a taxonomy theory [power system stability],' Proceedings of the IEEE , vol. 83, no. 11, pp. 1530-1561, 1995.
- [26] H.-D. Chinag and L. F. C. Alberto, Stability Regions of Nonlinear Dynamical Systems: Theory, Estimation, and Applications . Cambridge University Press, 2015.
- [27] E. Sorrentino, O. Salazar, and D. Chvez, 'Effect of generator models and load models on the results of the transient stability analysis of a power system,' in 2009 44th International Universities Power Engineering Conference (UPEC) , pp. 1-5, Sep. 2009.

## APPENDIX A

## 9-BUS SYSTEM PARAMETERS

The dynamical parameters are selected from the Appendices of [20] and have been converted to the same base power at 100 MVA.

TABLE II: 5-Bus Case Dynamical Data

|   Bus | Rated Power   |   X ′ d |      r |   T ′ do |     H |
|-------|---------------|---------|--------|----------|-------|
|     1 | 100 MVA       |  0.22   | 0.0035 |      5.9 | 4.985 |
|     5 | 147.1 MVA     |  0.2033 | 0.0023 |      4.3 | 4.31  |

## TABLE III: 9-Bus Case Dynamical Data

|   Bus | Rated Power   |   X ′ d |      r |   T ′ do |       H |
|-------|---------------|---------|--------|----------|---------|
|     1 | 233 MVA       |  0.1391 | 0.0007 |     5.14 |  9.5996 |
|     2 | 270 MVA       |  0.0948 | 0.0006 |     4.8  | 11.151  |
|     3 | 270 MVA       |  0.0948 | 0.0006 |     4.8  | 11.151  |

## TABLE IV: 14-Bus Case Dynamical Data

|   Bus | Rated Power   |   X ′ d |      r |   T ′ do |      H |
|-------|---------------|---------|--------|----------|--------|
|     1 | 330 MVA       |  0.0961 | 0.0005 |     6    | 9.9    |
|     2 | 147.1 MVA     |  0.2033 | 0.0023 |     4.3  | 4.31   |
|     3 | 192 MVA       |  0.1641 | 0.0019 |     5    | 5.952  |
|     6 | 100 MVA       |  0.22   | 0.0035 |     5.9  | 4.985  |
|     8 | 100.1 MVA     |  0.3137 | 0.0049 |     6.55 | 3.1201 |
Dan Wu Bin Wang
Influence of Load Models on Equilibria, Stability and Algebraic Manifolds of Power System Differential-Algebraic System
Load models have a great impact on voltage behaviors as well as power system
transient dynamics. Extensive work has been done on this topic, proposing
appropriate load models and capturing better load behaviors during transient.
This paper presents a comprehensive study to investigate the geometric and
topological changes induced by different load models for the traditional power
system differential-algebraic equations. Specifically, we attempt to reveal the
deformation of equilibria, stability regions, and algebraic manifolds during a
continuous evolution of load model. Several findings are presented in the
paper, some of which countering traditional recognitions and intuitions. A
major discovery is that the load model with a large proportion of constant
impedance and a small proportion of constant power exhibits much more complex
features than the load model with the reversed proportions of impedance and
power. The increase of complexity is thoroughly recorded and investigated by
the changes of geometric properties and mutations of topological invariants in
the sense of equilibria, stability regions, and algebraic manifolds for the DAE
system. However, most of the changes seem to occur on unstable components of
algebraic manifolds or near the singular boundary surfaces, suggesting a
limited variation of dynamical behaviors on the stable component.