w

'2"

'-).~4.-

\‘n ".\.. .

23.33;- ..

'4 r

m.“ -

O
A

‘0”?!

4'1r.

z

k “‘5'; u )u .

 

“Easy“, .

:r,
‘0-

"

o—‘V

...!. ‘
!‘("‘-‘.'f
' a Y:

 

 

u

, A. ‘ ‘ -, - pry-fined:
.'a.-v"l‘l'.:c -' 4" .v" . . ' V . '. 34;. "n‘-.‘- .
m.- ‘ . ~ , . , . ". - “gggw‘ffi .3— s»
., J 2‘ ~ . . ‘ .

1. "I
. .4 :532 rug
“"t‘s'i‘wfil; ,. ‘ ‘ .
' .(V*r'l‘.‘:7!;§ _-
9.3 "PM” ‘ I
Q...“ 3233392”;
"A.§.‘.W,:. J...

v 3"ng ..

1-‘

.
\v. ;
". '2

35:33.: ,
a -

L “J,
. ; 49:7,.
‘L’g;:k"t::p

. ,

’ 5
‘ ‘1 1'
'u ,

\1
“o .

s
“

‘ A: ,‘;. ,.-
( .1 iii-1255:"
' 5.: A I 4 ._
w»;- . ,

r 7-6 14.»:

A\'

"'3 .~ " “.1
K‘i’?’ .3"
' .3- v

I”;

”If n

, ”tr!

..

«a 2,;-
. .

i‘
thg

.3:
«“47
“u‘ ,‘ ;

£23425??- “5.45;?
{’5‘ - WNW? J-‘fi-"v‘fi’fi'.
' " 4.. no» r -;u

-' Hm “fiflfln
2M.- V- A. ‘
rfé‘éawr:
fair I,
\ e ’ m . u
.1“ 1;.

- 5- £4422.” .
212:1. I--”T'<" saw” :4“
'5‘- " a" '1 .
' J

.,..._. _
“.1.“ A:
,. ,

,.. .. “1..., ”
Z "’ . *‘3 " 4&3? .t: :
~< vine ~\l ‘. 4.. -,..,,_v(1r......-
. ~.*::=‘..7:...M., mic.“ -_ ..‘:“-....3 a
1M- ‘ -, VF”. .0409'1'333AWV. sn—pru: .L.£'___ .‘
um. m -- J ”bf”, z ...,.. ”g, L..-..~ V?
.1 y ‘9 1
.
J

'-

.'
Viz n J ‘1.

,
H . —;-4w'l
. "per J #53:: u
' .0- "! ‘
" {r 5‘ Jfiwflx
‘ R" 3- -">' ""L3-.--'.'.:'.'.w'.
. _ ,
. .. _ "‘r'L-Zu" .. '25
.
1'

W on-
c - a f‘
4"}:ch

N
"”47 -. .1;

IV raga-«3m . .rv

.u'.' 'v P1 ~

T!
Lfilés¢£¢ fun.“
":1. - my.

91".:

. n-
ZWJUMa")?
my..-
D‘
'3"
. .u
u. E
NI“ r‘
’3' ” 13": '

’3" ' , :r
’ ' dawn'nxxrr .7: ‘ .
Wait”! .1 rv~ vii!

 

 

williiiggiiilgimimuyul

éf‘ii'f3.;»

This is to certify that the

dissertation entitled
DETERMINING THE STEADY STATE SOLUTIONS OF

NONLINEAR MODELS OF POWER SYSTEMS
Homotopy Methods and Computer Implementation

presented by

SHIXIONG GUO

has been accepted towards fulfillment
of the requirements for

PH 0 D 0 degree in ELECTRICAL
ENGINEERING

1i Major professor A

Date 2 NOVQMLQI’ )4?0

 

 

 

MSU is an Affirmative Action/Equal Oppo rrrrr ' ty IMII!HUU'I 0- 12771

 

FT a

EJBRARY l
““1“!“ State
valves-guy j

 

 

L

PLACE IN RETURN BOX to remove We checkout from your record.
TO AVOID FINES return on or before due due.

F___________.____————-————————————1'l
DATE DUE DATE DUE DATE DUE

i i I
LT g

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative ActlorVEqueI Opportunity Institution
omens-9.1

 

DETERMINING THE STEADY STATE SOLUTIONS OF
NONLINEAR MODELS OF POWER SYSTEMS
Homotopy Methods and Computer Implementation

By
Shixiong Guo

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1990

(44%- 5%

ABSTRACT

DETERMINING THE STEADY STATE SOLUTIONS OF
NONLINEAR MODELS OF POWER SYSTEMS
Homotopy Methods and Computer Implementation

By

Shixiong Guo

Determining the steady state solutions (equilibria) of models of interconnected
power systems, known as the load flow problem, has become increasingly important,
and is presenting challenging problems, facing theoretical as well as applied research-
ers. The load flow problem has continuously received the attention of researchers due
to its essential role in the planning and operation of power systems. It is the core prob-
lem in the studies focusing on the stability and bifurcations of (models) of power sys-

tems.

In this thesis, we use powerful analytical tools and modern techniques from alge-
braic geometry to infer the number of the steady state (equilibrium) solutions of non-
linear models of power systems. We develop the theorems and the methods to predict
and determine the steady state (equilibrium) solutions for various levels of detailed
models of power systems. Sufficient conditions are provided which guarantee the pre-
cise number of solutions to the full-fledged load flow. The sufficient conditions are
cast in terms of properties of the physical admittance matrix of the power grid. Con-
sequently, these sufficient conditions are. placed on the structure or the topology of the
given power network. When the sufficient conditions are not satisfied, we develop the
cluster method to provide a "tighter" upper bound on the load flow solutions for spe-

cial power grid structures. Consequently, our results lower the upper bound for the

Shixt'ong Guo

many practical power grids which are normally sparsely connected

We also present the special homotopy method, due to Li et al, to reduce the com-
putational complexity, and to "guarantee" the finding of all possible solutions of the
load flow equations of power systems. We then develop the imbedding-based and the
homotopy-based heuristic methods to simplify the computations in finding the solutions
of the so—called deficient systems (with particular interest to power systems). More-
over, the methods render procedures which are directly implementable on digital serial
and parallel processors.

We specialize some of our results to prototype models and numerical examples of
power systems to illustrate as well as demonstrate the procedures capabilities. The
algorithmic techniques are then implemented to obtain the steady state (equilibrium)

solutions of various models of power systems.

To

my wife, Xuemei; my son, Zheng; and my Motherland.

iv

ACKNOWLEDGMENTS

I would like to give sincere thanks to my academic advisor, Dr. Fathi M.A.
Salam, for his continuous advice, encouragement, and support throughout my graduate
studies and my stay at Michigan State University.

For their interest and helpful suggestions, I also wish to thank the other members
of my committee: Dr. Hassan K. Khalil, Professor of Electrical Engineering; Dr.
Robert A. Schlueter, Professor of Electrical Engineering; and Dr. David Yen, Profes-

sor of Mathematics.

From the bottom of my heart, I would like to thank my wife, Xuemei, and my
son, Zheng, for their love, support, and understanding. I also owe a special thanks to

my parents for encouraging me to pursue higher education.
Finally, I wish to express my thanks to Dr. Lionel M. Ni, Professor of Computer
Science; and my colleague, Dr. X. Sun; for their helpful discussions in the Power Sys-

tem Research Group.

TABLE OF CONTENTS

LIST OF TABLES ................................................................................ ix

LIST OF FIGURES .............................................................................. x

NOMENCLATURE .............................................................................. xi

Chapter 1. INTRODUCTION ....................................................... 1

1.1. The Load Flow Problem: basic issues ..................................... 2

1.2. Scanning the Literature ............................................................. 3

1.3. The Contributions of This Thesis ............................................. 8
Chapter 2. (NONLINEAR) MODELS OF

POWER SYSTEMS ..................................................... 10

2.1. The Classical Model ................................................................. 10

2.1.1. The polar-coordinate form .............................................. 13

2.1.2. The rectangular-coordinate form .................................... 14

2.1.3. The zero-conductance form ............................................ 15

2.1.4. The decoupled load flow equations ................................ 16

2.1.5. The complex form ........................................................... 17

2.2. The Model with Internal and Terminal

Buses of a Generator ................................................................ 18

2.3. The Model Augmented by the Excitation System ................... 23
Chapter 3. ON THE NUMBER OF (EQUILIBRIUM) STEADY

STATE SOLUTIONS OF POWER SYSTEMS ........... 30

vi

3.1. A Deficient System and Its

Associated Homogeneous System ........................................... 31
3.2. Nonsingular Zeros on Projective Spaces ................................. 33
3.3. Application to Quadratic Polynomial Systems ........................ 35
3.4. The Classical Model of Power Systems .................................. 39
3.4.1. The all-PQ—bus load flow equations ............................... 40

3.4.1.1. Example: a fully-connected

3-bus power system structure ............................. 42
3.4.2. The all-PV-bus load flow equations ............................... 46
3.4.3. General load flow equations of power systems ............. 52

3.5. The Model with Internal and Terminal

Buses of a Generator ................................................................ 55

3.6. The Model Augmented by the Excitation System ................... 57
3.7. A Structure of Special Power Systems .................................... 59
3.7.1. The cluster method .......................................................... 59
3.7.2. Examples of power system networks ............................. 62
3.7.2.1. A not-fully-connected 3-bus network ................. 62

3.7.2.2. A 7-bus network ................................................. 62

3.7.2.3. The model with internal and terminal

buses of a generator ............................................ 65

3.8. Summary ................................................................................... 65
Chapter 4. THE SPECIAL HOMOTOPY METHOD .................. 67
4.1. The Basic Homotopy Method .................................................. 68
4.1.1. Applications to power systems ....................................... 71
4.1.1.1. A 3-bus numerical example ................................ 72

4.1.1.2. A 5-bus power system network .......................... 74

vii

4.2. The Special Homotopy Method ............................................... 77
4.2.1. The classical model ......................................................... 77
4.2.2. The model with internal and terminal

buses of a generator .................. ' ...................................... 87

4.2.3. The model augmented by the excitation system ............ 90

4.3. Summary ................................................................................... 94

Chapter 5. THE IMBEDDING-BASED METHOD SOLVING

FOR ROOTS OF DEFICIENT SYSTEMS ................ 95

5.1. The Basic Irnbedding Method .................................................. 96
5.2. The Imbedding-based Method:

Practical Heuristic Approach ................................................... 99

5.3. Numerical Examples of the Load Flow Equations .................. 101

5.3.1. A not-fully-connected 3-bus

power system network .................................................... 101

5.3.2. A 7-bus power system network ...................................... 104

5.4. Summary ................................................................................... 107
Chapter 6. THE HOMOTOPY-BASED METHOD TO

DEFICIENT SYSTEMS ............................................... 108

6.1. The Practical Heuristic Approach ............................................ 108

6.2. Numerical Examples of Power Systems .................................. 113

6.3. Summary ................................................................................... 118

Chapter 7. CONCLUSIONS AND SUGGESTIONS .................... 119

7.1. Conclusions ............................................................................... 119

7.2. Suggestions ................................................................................ 120

Appendix The Proof Of Theorem 3.4.3 ......................................... 121

BIBLIOGRAPHY ............................................................................ 129

viii

LIST OF TABLES

Table 4.1. The solutions of the 3—bus example (16 initial points) ............ 75
Table 4.2. The solutions of the 5-bus network (256 initial points) '. .......... 78
Table 4.3. The solutions of the 3-bus example (6 initial points) .............. 82
Table 4.4. The solutions of the 5 -bus network (70 initial points) ............. 83
Table 4.5. The solutions of the 7-bus network (924 initial points) ........... 86
Table 4.6. The solutions ofthe 4-bus example (20 initial points) 89
Table 4.7. The solutions of the 4-bus network

with an excitau’on system (70 initial points) ........................... 93
Table 5.1. The solutions of the not-fully-connected 3-bus example

by the imbedding-based method (4 initial poinm) ................... 105
Table 5.2. The solutions of the 7-bus network by the

imbedding-based method (288 initial points) .......................... 106
Table 6.1. The solutions of the 4-bus example by the

homotopy-based method (12 initial points) ............................. 115
Table 6.2. The solutions of the 7-bus network by the

homotopy-based method (288 initial points) ........................... 117

Figure 2.1.
Figure 2.2.

Figure 2.3.
Figure 3.1.
Figure 3.2.
Figure 3.3.
Figure 3.4.
Figure 4.1.
Figure 4.2.
Figure 4.3.
Figure 4.4.
Figure 4.5.
Figure 4.6.
Figure 5.1.

LIST OF FIGURES

The main (components of a power system. ............................. 12
The main components of a power system including

the direct axis synchronous reactance. .................................. 21
Block diagram of an excitation system. ................................. 24
The 3-bus power system structure. ........................................ 43
The partitioning of a network. ............................................... 60
The not-fully-connected 3-bus structure. ................................ 63
The partitioning of the 7-bus network. .................................. 64
The homotopy curves. ............................................................ 70
The fully-connected 3-bus example. ....................................... 73
The 5-bus network. ................................................................. 76
The 7-bus network. ................................................................. 85
The 4-bus example. ................................................................ 88
The 4 -bus network with an excitation system. ....................... 91
The not-fitlly-connected 3-bus example. ................................. 102

Nomenclature

moment of inertia.

damping coefficient.

angle at generator buses.

angle at buses excluding the slack bus,
angle at load buses.

mechanical power,

active power injection at generator buses,
active power injection at load buses,
reactive power injection at load buses,
complex power injection at buses,

direct axis transient open-circuit time constant .
direct axis synchronous reactance,

the exciter output voltage,

quadrature axis magnitude of voltage behind
transient reactance.

quadrature axis magnitude of voltage behind

synchronous reactance,

CI"I

the flux linkage of the field winding.
voltage magnitude at buses.

the reference voltage,

a voltage enor,

the measured terminal voltage,

the stabilizer output voltage,

the amplifier output limit voltage,
the saturation function,

the amplifier gain,

the amplifier time constant.

the time constant measured through a potential

transformer, rectified and filtered.
the stabilizer gain,

the stabilizer time constant,

the exciter gain.

the exciter time constant,

the homotopy function with homotopy parameter 0951,

the ”initial" starting system in the homotopy
function, i.e., H (.,O):=S (.),

the target system in the homotopy
function, i.e., H(.,1):=T"(.),

the projective space.

xii

Nomenclature

Chapter 1
INTRODUCTION

Determining the steady state constant solutions (i.e., equilibria) for any dynamic
system is the most basic and fundamental problem in the quest for analyzing and
investigating a system’s behavior. Once the equilibria are determined, one may calcu-
late the linearization of the vector field about a chosen equilibrium point, i.e. the Jaco-
bian matrix, and then calculate the eigenvalues of the Jacobian matrix. By Lyapunov’s
direct method, the eigenvalues generically infer the local stability or instability of the
chosen equilibrium point. These types of calculations have generally been programmed
onto digital computers. Alternatively, one may use the Lyapunov function techniques
within a neighborhood of the equilibrium point to analytically infer the local stability
or instability of the equilibrium point.

For many systems, an equilibrium point represents the operating point of that sys-
tem if it is (asymptotically) stable. It may otherwise belong to the boundary of the
basin of the stability (or attraction) for a certain desired operating point; in this case,
the equilibrium point may characterize a part of the boundary within its vicinity. (The
described view was specialized for the particular case of power systems in [7,8,10].)

The question is then: how to determine the equilibria?

Various numerical procedures have been used over the years with limited degree

of success. Examples of these numerical procedures are Newton-Raphson, Gauss,

Seidel, Poncel Fletcher etc.. All these procedures, however, suffer flour a "poor"
choice of initial conditions. Variations of Newton-Raphson procedures are common
and are believed to be frequent in use. Indeed instability or non-convergence of these
procedures have been reported sporadically in the literature. Yet, for a lack of better

alternatives, the use of these procedures have continued to dominate.

Recently, a study [11] by Thorp and Naqavi has articulated a form of instability
associated with the use of the Newton-Raphson method in determining equilibria for
the polar-coordinate models of the swing equations of power systems. The conclusion
of [11] illustrates via numerical computations of some power system examples that
each isolated equilibrium point attracts basins of initial guesses. However, the basins of
attraction for some equilibrium points can be very irregular, and/or extremely small.
The assertion in [11] is that fractal structure are present in their power system exam-

ples.

It would be a welcome relief therefore to develop numerical procedures or algo-
rithms that are independent from the choice of initial conditions of the algorithms and
would succeed despite a "poor" choice of the initial conditions. It would even be more
welcome if such numerical procedures are ggaranteed via mathematical foundations to

successfully find all (or some of) solutions.

1.1. The Load Flow Problem: basic issues

The load flow (or power flow) problem is the calculation of line loading given the
generation and demand levels for the normal balanced three-phase steady-state operat-
ing conditions of an electric power system ([1]-[5]). Basically, the problem can be
intuitively described as follows: given the forecast (real and reactive) load demands,
the secure operation of power systems entails determining the required power genera-
tion so that a known set of inequality constraints are satisfied. The load flow equa-

tions are a system of nonlinear models which relate the real and reactive power, or the

real power and voltage magnitudes at each node or bus in an electrical network operat-
ing in steady state. In general, load flow calculations are performed routinely for
power system planning, and in connection with system operation and control. The data
obtained from the load flow are used for the study of normal operating mode, con-
tingency analysis, outage security assessment, as well as optimal dispatching, stability
and bifurcations.

Although the full load flow equations have extensively been used for a long time,

there remain a number of very basic open questions.

(a) What is the least upper bound on the number of (complex) solutions of the full
fledged polynomial load flow?
(b) The more important but difficult problem is how to solve and obtain all (or some

of) the equilibria of the nonlinear models of power systems.

(c) What are the number of real system solutions of the load flow equations for a

given N-node power system?

These challenging questions have fascinated and puzzled mathematicians and engineers
for many years. Many researchers have devoted large part of their professional lives to
such problems. This thesis will address and provide answers to the first two

significant questions.

1.2. Scanning the Literature

Because of its difficulty, few theoretical investigations of the load flow equations
have been reported and with rather limited practical results ([12]-[17], [28], [30], [31]).
For instance, it is not yet possible to infer the number of solutions of the full-fledged

PV- and PQ-bus load flow equations.

The possible existence of multiple (stable) solutions to the load flow equations
has been realized both analytically as well as via simulation of realistic power net-
works ( see, e.g., [55], [57], [58], [74]). Related works using computational methods

have appeared in [47], [52], and [56] though they emphasize the security aspects of
power systems. Works that emphasize the optimal load flow computations are available
in [42], [53], [61], [62], [68], and [71]. Works that emphasize the stability and bifur-
cations of power systems have appeared in [6]-[10], [16], [17], [26], [29], [33], [41],
[49], [59], [60], [65], and [70]. Some theoretical analysis of various simple models of
the load flow have been initiated by Galiana ([50], [51]). Analysis on the real power
of the load flow equations began with the work of Tavora and Smith [31]. Recent
investigations made by Arapostathis, Sastry and Varaiya [27] have utilized tools from
bifurcation theory. The most recent work on the real power flow equations have been
reported by Baillieul and Bymes ([15]-[17]). As for the full fledged power flow equa-
tions, the work by Wu and Kumagai, namely [7 2], [73], represents the only analytical
treatment. Yet, none of these theoretical works determines upper bounds on the
number of solutions for a general N—node (PQ and PV buses) power system. None
addresses the dependence of the number of solutions on the network structure in addi-
tion to the number of nodes. More importantly, none of them addresses guaranteeing
the finding of all possible solutions of the load flow equations. The reason stems from
the difficulty of the problem analytically even though the mathematical tools employed
in the works ([15]-[17], [27], [46]) are sophisticated.

In the early work on multiple equilibria of swing equations, i.e. on multiple solu-
tions of the load flow equations restricted to be lossless PV buses, for an N -node
(excluding the slack bus) M power system model, Prabhakara et a1 [12] have
stated that the number of equilibria of the swing equations for an N -node (excluding

 

the slack bus) lossless power system model is 2”. The upper bound 2N for an N-node
power system (excluding the slack bus) has been rationalized as follows. Consider N
decoupled nodes connected only to the reference or infinity node. If one assumes that
each node connected to an infinity node produces 2 possible solutions, then the N
decoupled nodes will give rise to 2N possible solutions. It is also assumed that the

coupling would not increase the possible number of solutions. This is the best rational

for an otherwise a completely heuristic justification.

In the following extension on the number of solutions of the load flow equations,
Tamura et al ([13], [14]) have concluded that the number of solutions of the load flow
equations for an N-node power system (excluding the slack bus) is 2”. Tamura et al
have developed both basic and simplified algorithms to compute multiple solutions to
the load flow equations. These algorithms are based on the Newton-Raphson method
(with some optimal multiplier which is assumed to prevent divergence and occurrence
of oscillation) to compute the multiple solutions to the load flow equations. Even so,
there is no guarantee that these algorithms would find all possible solutions to the load
flow problem since the algorithms are based on the incorrect conclusion that the upper
bound on the number of solutions of the full load flow equations for an N-node power
system (excluding the slack bus) is 2N .

The subsequent works on the number of solutions of the load flow equations by
Baillieul and Bymes ([15]-[17]) have used some powerful results from Algebraic
Geometry, the Morse Theory and the Intersection Theory to study the load flow equa-

tions of the lossless PV-bus power system models. These theories were specialized to

 

the lossless PV-bus models of power systems cast as a special quadratic system of
polynomials. Their work has extensively exploited the particular features of the loss-
less PV-bus model. They have finally concluded that the number of complex solutions
of the load flow equations for an N-node (excluding the slack bus) power system (with
all buses (i.e. nodes) restricted to be PV buses) is bounded above by

N» = [if].
The most recent works by Li, Sauer, and Yorke ([18], [19]) have used different

tools fi'om Intersection Theory, and the Homotopy Continuation methods to develop

Homotopies to find all the solutions of general classes of polynomial systems. They

have developed some theorems that can be used to determine the upper bound on the
number of solutions for general classes of (generic) mlygomial systems. Some of
their applied results have stated that the number of solutions of the load flow equations
for an N-node (excluding the slack bus) power system (with all buses (i.e. nodes) res-
tricted to be PQ buses) is bounded above by ‘
N. = [ill
It appears that the fi'amework of Li, Sauer, and Yorke is more general and is sys-
tematic. It only applies to general (generic) polynomial systems however. It should be
emphasized that neither work is applicable to a general power system which includes
both PV and PQ buses. It would be of interest, however, to prove that this bound
applies to general PV- and PQ-bus power system. Moreover, for more detailed models
that include the excitation system, e.g., one has yet to develop representations that
describe the overall model in terms of systems of polynomials. Only through the gen-
eral polynomial systems approach of [18], [19] as well as our newly developed
representations have we been able to extend, then apply, the Homotopy method tech-

niques to general power system models.

Because of the nonlinearity of the models of power systems, numerical pro-
cedures have continued to be the only possible avenue for obtaining all (real) solutions
of the load flow problem [25]. However, traditional methods, such as the Newton
method and its variations, are not capable of solving for all the roots. The popular
IMSL package, for instance, can only find one root based on the MINPACK imple-
mentation of MJD. Powell’s hybrid algorithm and P. Wolfe’s secant method [54].
The Newton-Raphson method is also known to fail in many case studies. It is well
known that N ewton-Raphson method can be unpredictable when the initial guess is
poor, and it breaks down when the Jacobian at any stage of the iteration becomes

singular. Thence this method does not always converge to a solution, and if it does,

there is no guarantee that it would find all possible solutions to the load flow problem.
The authoritative review on the subject of computational load flow equations is still
perhaps [25]. Various subsequent works have employed variations of the Newton

method and have used parallel processing concepts ([48], [66], [67], [69]).

Solving for all the roots of any system of polynomial equations has been almost
impossible until the advent of the globally convergent probability-one basic homotopy
method (subsequently, referred to as the homotopy method) ([20]-[24], [35]-[39], [43]-
[45]). The method is globally convergent in the sense that it will cOnverge to solu-

tions of the problem from an arbitrary set of starting initial points ([2]-[5]).

Homotopy continuation methods have been proven to be superior to the quasi-
Newton methods ([76]). A recent survey of this method can be found in [77]. Here
we briefly describe the homotopy method as applied to solving systems of polynomial
equations and underline some of its properties. First, the numerical computation of the
homotopy method can be systematically implemented in parallel processors. Second,
the homotopy method is globally convergent, i.e., one may choose any initial guesses
and the homotopy method is guaranteed to converge to all solutions with probability

0116.

Since each homotopy curve only depends on the "initial" starting point, we can
trace the homotopy curves separately. This makes it possible to exploit the inherent
parallelism in the (polynomial) load flow model to take advantage of massively parallel
computers ([45], [63]). Within the last five years, various types of parallel machines
have been commercially produced [64]. Large-scale scientific computing is one of the
major application domains which demands the huge computing power of parallel
machines. However, the move into parallel territory requires new conceptual strategies
in formulating a problem, and new algorithms to shape the problem for parallel com-
puters. This implies that a brute-force approach to solve a large-scale problem, such

as the load flow equations of power systems, on parallel computers does not render the

problem tractable.

1.3. The Contributions of This Thesis

In this thesis, we use some powerful analytical tools and modern techniques from

algebraic geometry and the homotopy continuation methods (with their algorithmic

implementation onto computer) to investigate and calculate the steady state solutions of

various nonlinear models of power systems. In particular,

(A) We extend the model of the (equilibrium) steady-state equations in the appropri-

(B)

(C)

ate (complex) polynomial representations in chapter 2. This representation form is

convenient for the homotopy approach we pursue.

We develop theorems in chapter 3 to investigate and predict the number of solu-

tions of the full-fledgg (equilibrium) steady-state equations for various levels of

 

detailed models of power systems. Sufficient conditions are provided which
guarantee the precise number of solutions to the load flow. The sufficient condi-
tions are cast in terms of properties of the physical admittance matrix of the
power grid. Consequently, these sufficient conditions are placed on the structure
or the topology of the given power network. When the sufficient conditions are
not satisfied, we describe the cluster method to provide a "tighter" upper bound
on the load-flow solutions for special power grid structures. The upper bound,
therefore, depends on the topology or the structure of the power grid in addition
to the number of nodes. Consequently, our results reduce the upper bound for the
many practical power grids which are normally sparsely connected.

We present the special homotopy method in chapter 4 to reduce the computa-
tional complexity, and to guarantee the finding of all possible solutions of the
load flow equations of the polynomial power systems with probability one. The

amount of computation depends on the size of the system.

(D) We develop the "imbedding-based" method in chapter 5 to simplify the computer

CE)

computations of finding the solutions of the load flow equations, and to make the
algorithm capable of handling relatively large-sized power system networks. This
method, largely heuristic, draws from experience with the Homotopy procedures
and from properties of power systems.

We develop the "homotopy-based" method in chapter 6 to further reduce the
computational complexity in the finding of the geometrically isolated roots of

deficient systems, such as power systems.

The final conclusions and suggestions are given in chapter 7.

Chapter 2

(NONLINEAR) MODELS OF
POWER SYSTEMS

The mathematical model of representations of a power system should be chosen
or developed to accommodate, and be accommodated by, the techniques and the com-
putational facilities available. In the following, we present various representations of
power systems that can be used as per convenience. From the view point of the
homotopy-based or the imbedding-based methods, the representations that are in the

form of systems of polynomials are of interest.

2.1. The Classical Model

The well-known classical model is formulated by assuming that the flux linkage,
hp,- , of the field winding is constant, and a voltage regulator holds the magnitude of its
terminal voltage fixed by automatically varying the generator field excitation. We
have listed the nomenclature at the beginning of this thesis (see page xi). The model
of a power system consists of three main components (see Figure 2.1): generators,
loads, and a transmission network that connects generators and loads. Let the model
of the power system consist of N+1 buses or nodes. Let the generator buses be sub-
scripted from 1 to N8, and let the load buses be subscripted from Ng-I-l to N. We
choose to subscript the slack bus by N +1. We denote the node or bus complex admit-
tance matrix by [Y] where its ki-th component is yh- = G,“- +13,“- (j = 4:1). The

10

ll

term G,“- is the conductance of the line connecting buses k and i, and B,“- is the sus-
ceptance of the line connecting buses k and i. Each load is represented as a constant
real (Pd ) and reactive (Q4 ) power demand. Therefore, we refer to a load bus as a PQ

bus.

In the general case, a change in input, or load, or structure, or a sequence of such
changes, causes dynamic motion. If we ignore the system direct axis synchronous
reactances xd , then the motion of the i-th generator with round rotor is govemed by

the following swing equations

0=é,-to,, i=1,2,---,N8, (2.1.1a)
o = Mid),- +D,-m,- +F,-G(e,V)-R,-"', (2.1.1b)
t =1, 2 , , ”3’
with constraints
0 = Ef(e,V)—R,L, k =Ng+1,---,N, (2.1.1c)
o = Tf(e,V)-Q,{-, k =Ng+1,°--,N. (2.1.ld)

The load flow equations or the equilibrium equations are obtained by setting the

derivative terms to zero, that is,

0 = FiG(9,V)-P{", i=1,2,~‘,N8, (2.1.2a)
o = F,[-(e,V)-P,{-, It =N8+1,--°,N, (2.1.2b)
o = Tf(e,V)-Q,{-, lt =Ng+l,--°,N. (2.1.20)

The power flow, or load flow, problem has been formulated based on sinusoidal
Steady state nodal analysis of circuit theory. Load flow calculations are performed in

12

Ear? 833 a mo 3:259:50 Ea:— ofi. . fl .N Ouzwgnm

 

 

@

~+wz

983-0%
mews

 

 

#838 Z

eemmflfimefih

 

2:. Scam

Jo

 

._+Z

c.3353

 

 

 

O
QZO 388050
A

@

 

 

13

power system planning, operation and control. They are increasingly and extensively
being used to solve very large systems, to solve multiple cases for purposes such as
outage security assessment and voltage collapse, and as part of more involved calcula-

tions such as optimization, stability, and bifurcations of (models of) power systems.
There are three types of buses in a power system network:
(i) PQ bus: a bus where the real and reactive powers are specified.
(ii) PV bus: a bus where the real power and the voltage amplitude are specified.

(iii) A slack bus: a fictitious concept whereby one of the generator buses has only
its complex voltage specified. One purpose of this bus is to guarantee that the
total power injection into the network equals the total power consumed. Conse-

quently, ensuring the existence of steady state solutions (i.e., equilibria).

It is conventional to model loads as PQ-buses, one generator as a slack bus, and
the rest of the generators as PV buses. All phase angles are measured relative to the
angle of the slack bus which is set to be 0”“ = 0; furthermore the slack bus voltage is
set to be VN+1 = l per unit. In the following, we derive various forms of the power
system (equilibrium) steady-state models that are convenient representations for the

analytical development treated later.

2.1.1. The polar-coordinate form
In the polar coordinate representation, the complex voltage at the k-th bus can be
expressed as:
V,‘ eje‘ = V,‘ c059,, +j Vk sine,“
where V,‘ represents its voltage amplitude and 9* represents its phase angle. Let 9,“-
denote 6k - 9,- (lskjSNH).

At the k-th generator node, we have [40]

14

O = EiileV‘IGb'COSGh' + Basins“) - PF, (2.1.33)
k = 1, 2, ,N8.

At the k-th load node, we obtain

0 = giflvkvgoucoseb stab-sine”) - Pf, (2.1.3b)
o = Zf‘VkVi(Gb-sin65 - Bficoseu) — of, (2.1.3c)
k =Ng+1,---,N.

2.1.2. The rectangular-coordinate form
In the rectangular coordinate representation, we have
Vk cosGk-t-j Vk sinek =:)?k +j Yk, k =1,2,--°,N.

For convenience, we use the notation Pk to denote the real power injected at the
k-th bus including both the generator bus and the load bus. similarly, we use Q,‘ for

the "imaginary" power at the k—th bus. Therefore, the load flow equations can be as:

o = 1;;1[a,,(i,ri, +r,r,)+o,,ot,.r, 4,9,1] 4,, (2.1.4.)

o = 23.9343, to = 1,2, - - - ,Ng, (2.1.4b)
0 = £I1[Gfl(iifk 'thi) 'Beozltii + YJD] - th (2.1.46)

k=Ng+l,---,N,

15

01'
o = xksz (Gh-X- -BhY)+YkZ,-’!+1(35X;+thi) -P,,, (2.1.4a’)
k=l,2, ,N,
o = fluff—v.2, k =1,2,---,Ng, (2.1.4b’)

o = ykzN+l(G,,x— —E,,,f)— —X,);N+1(B,,,Ii,+G,,-1?,-) -Q,, (2.1.4c’)

k=Ng+1,-~-,N,

where 23,, represents the real part of the complex voltage E,‘ , while 1?; represents its
imaginary part. The equations (2.1.4b) and (2.1.4b’) represent the constraints on the

voltage amplitude of the PV-buses.

2.1.3. The zero-conductance form

In this form, the transmission line conductances are assumed negligible. Conse-

quently, we set Gh- = 0 in (2.1.4) to obtain the simplified model as follows:

0 = N+1[Eh(x 19,- 42,179] -P,,, k = 1,2, - - - ,N, (2.1.5a)
o = 123+ 1"?- v}, It = 1,2, - - . ,N,, (2.1.511)
o = N+1[Eh(x,,1f + 17,,1’ )] +Q,, k= N +1, ,N, (2.1.5c)
01'
o = nzN; lahx- -x,z,.N+IE,,,.1",. -P,,, It =1,2,---,N, (2.1.5a’)
_ “2 A2 2 _ . . . 9
0 - Xk +Yk -Vk9 k-1,2, ’Ng’ (2.1.5b)

o = x,z,N+IE,,-x +Y‘,2,.N+IE,,1?,- +Q,,, k =Ng+1 , . - - ,N. (2.1.5c’)

16

2.1.4. The decoupled load flow equations

In addition to assuming that Ga = 0, it is often assumed that the phase angle
difference between nodes (or buses) are small. Consequently, one may write
sine,“- = 6k,- and cosOb- = 1. These two assumptionsnow reduce the polar representa-

tion (2.1.3) to the so called decoupled load flow equations, namely

0 = v, Ntlvmb- 9,, -P,,, It =1,2,---,N, (2.1.6a)
o = V, [LN;IV,E,,, 1+Q,, It =Ng+l,---,N. (2.1.6b)

We remark that the simplified models, namely, Models (2.1.5) and (2.1.6), can be
used to obtain approximate solutions for the full-fledged model. The motivation for

doing so is three-fold:

(a) The solutions of the approximate models can be computed efficiently (through
parallel processing which will be tailored for each of the simplified models).

(b) A simplified model, with its solutions, may be used as the "initial configuration"

solution- set required in a given homotopy method.

(c) The (parallel processing) techniques developed for the simplified models are
important in their own right, since the simplified models are employed in various
applications ranging from planning to (transient) stability.

The next model is a new model that was introduced in [3, 4]. The model uses the
complex space directly and hence one does not need to go though the necessary steps
of complexifying the space to permit (complex) analysis and/or facilitate the computa-
tional procedures for solving the roots of systems of equations. We call this form "the

complex form".

17

2.1.5. The complex form

Let E, denote the complex voltage of bus k. It is of the form
E ,, = V], cosek + j V,, sine,“ where V,, represents the amplitude of - the voltage at bus
k. Let 5,, := Pk + ij denote the complex injected power at node k, then the injected
complex power balance equation can be expressed as [E t] [Y] E - S ‘ = 0, where [E]
= diaglEl, - - -,E,,,],1~: = [E], ~--,E~]T, s = [3,, - - . , SN11', and the superscript
* denotes the complex conjugate.

We will give the 2N equations that govern the load flow equations in each of the
following cases [3,4]:

(a) The load flow equations of the PQ-bus Network:
El liflyzE; - Sk = O, k =1, 2 , ' ‘ ' , N, (2.1.73)
E;2N;1y,,,E, —S,,‘=0, k=1,2,°-',N. (2.1.711)

(b) The load flow equations of the PV-bus network:

E,2N;1y;E,+E;z,N;y1,,E, -2R,=o, lt=1,2,---,N, (2.1.8a)
E,E,,‘—V,,1=o, ‘ k=1,2,---,N. (2.1.8b)

(c) The load flow equations for the general power network:

E,z,.N;1y,;-E, + E,,‘2,N;1y,,,E,- -2R,, = 0,1 = 1, 2 , - - - ,N,, (2.1.911)
E,E,,’-V,,1=o, It = 1,2 , - - - ,N,, (2.1.911)
'E,2;N;1y,‘,E,-‘ - S, = o, It = N,+1 , ~ . - ,N, (2.1.9c)

51:21:11th1- S; = O, k =Ng+l , ‘ ' ' , N. (2.1.9d)

18

2.2. The Model with Internal and Terminal Buses of a Generator

This model includes one circuit for the field Winding of the round-rotor machine.
Let the i-th (lSl' Sn) generator’s terminal voltage be denoted by V,eje‘, and its internal
stator voltage be Egg/'5‘, or its internal stator-based voltage be Eq,ej& [40]. Let the
1.111 (n +2SkSn-I-m-I-l) load-bus voltage be denoted by V,e1'°'. The motion of the 1-111

generator with round rotor is governed by the following swing equations [40].

o=3,-to,, i=1,2,-~,n, (2.2.111)
0 = Miéi—Pim't'Dimi‘l'Pf, i=1,2,°°°,n, (2.2.1b)
o = T,’,,,.E,,',+E,, -EF,., i=1,2,~-,n, (2.2.1c)

with the algebraic constraints

 

0 = P,“ -g,,(6,E, ,9 ,v,<1>), (2.2.111)
o = (2,;1 - h,,(8,E ,e,v,<t>), (2.2.111)
k=1,---,n,n+2,~--,n+m+1,
where
P1‘ = ft(6 .E,.e,v,¢) = E,,v,-s1:5,-e,) (2.2.211)

bdt'eq' Vi 5111(51‘ ‘91 )1 (bot = l/xd, )1

i=l,2,-°°,n.

19

At the i-th (ISi Sn) terminal generator bus, we have

C
I

_ 81‘(51eq61‘/1¢)

b4; Eqi V; 8111(6,’ '8‘) +
Zféllgij Vi Vj 005(61' 'ej) + bij V; Vj sin(6,--9j )] +

"+111“ [g ,7, V,- V,, cos(6,- 4111) + but Vi VI: 3111(91' '4’" )1 ’

=n+2

and

o = h,(5,E ,9,V,<b)
= 114,11, 1v, - E4,- cos(e, -5,. )1 +
21311 [81; V1 Vj 8111(91' '91) - bij Vi Vj 008(91 '9 j )1 +
$311513", V,- V,, sin(9,- —<1>,,) - b,,, V,- V,, cos(e,- —<t>,, )1.

At the k-th (n +2SkSn+m+1) load bus, we obtain

Pkd 81(5.Eq.9.V.¢)
= ”:11 [81d Vlt Vi 005(411 ‘91) + be Vt Vi sin(<b,, '91 )1 +

an+l[gu Vk VI COS(¢k 41) + by Vt VI sin(¢k 41 )1,

=n+2

and

de = hk(8’E 969V9¢)

(2.2.2b)

(2.2.2c)

(2.2.2d)

20

= 231mb- Vk Vi Sin(¢k-ei) - bk. Vk Vi cos(¢k “6" )] +

mpg, V, V, sin(<1>,-<b,) - 11,, V,V,cos(¢,-¢,)].

Let the complex variables E,, E,- and E, be defined by
E,- := E,,,-cosb,- +j Eq,sin5,-, (j = \CI)
,- := VicosG, + j V,sin9,,

i=1,2,°°°,n,

and

Bk 2: VkCOS¢k +j VkSIn‘Dk,

k=n+2,~°,n+m+l.

(2.2.2c)

The static equations of the model (see the depicted equivalent schematic diagram

in Figure 2.2) in polynomial representations are obtained by taking all time-derivatives

to be zero. Therefore

0 = P,-’"—P-‘

1

51(5). - ET”; -

=p,_g

1 - ~ .
351‘ (E.- — 5.11.1. on é —,o,,.).

i=1,2,°~-,n,

(2.2.3a)

21

.025an 39850:? mg 886 05

wficfiofi Samba, 530m a mo ficocanoo 5.2: 2; .N.N own—wrm

 

 

933-3
mused

 

 

 

N+c

 

36+:

#8332
:owmmmamcfih

 

as seam

 

T—1:}—1I- coo T—[j—ln. flu

@

TIM

lemw
C

$33.?6
£98050

8% H
. ~

 

 

 

 

22

and

o = EEE- -E,.~2,-, (£4,- =12”), (2.2.3b)

At the i-th (151' Sn) generator terminal, one obtains
o = g;(E,E’) (2.2.3c)
1 ~ «It —¢ e
= 35KB} ‘55))2: +
l u ~ —
"2'53 (5: " Etb’di +
—(Ei 2:: lqu, 4' Bi EE‘ilygEl) +
"'05; 33351?th + Ei Emilyikik)’
and
o = may) (2.2.3a)
1 . u -t e
= ZEi<Ei 'Ei))’di "

l u . -
2—jEi (E,- - Ei))’di +

"2—(85 n=+1y£;E[ “'E 21:; lqul)+

23
1 . .
"27(Ei "$151)" tkEk " Eiflm21yfiEk)°
At the k-th (n +2SkSn+m+1) load bus, one has
0 = Pf - gk(E, E‘) ' (2.2.3e)

= " -(Ek25"=il)’a5i. + Ek2in=+IIYHEi ) ‘

-(5221m2 lYuEz + 51:21:32 lYuEz ):

and
o = Q: - ME, 5‘) (2.2.30
= Q4 _ —(Ek2¢n=+lykiE: _ Ezzzflthi ) __

72(5): milYuEt " Ekzmzlmiz).

2.3. The Model Augmented by the Excitation System

In the classical model of the swing equations, we simplify the model by assuming
the flux linkage to be constant. This simplified model can be used to obtain approxi-
mate solutions for the full model. In some cases, experience indicates the simplification
of the model to be reasonable and convenient. It is generally suggested, in transient
stability studies, that generators close to the fault should be modeled in greater detail.
The complexity of a generator model is completed by adding the flux decay dynamics

in the rotor windings and the excitation system dynamics [75] (see Figure 2.3).

The action of the excitation system of the i-th (151' Sn )generator can be described

by the following differential equations

.889? coufimoxo an we Eaewflu x85 .m.No.5wmn~

 

 

86.56 9793 #86 Sad nod. no 97 9—... on-n~ NdéQo 86¢

 

 

 

 

 

 

 

 

 

 

 

KM «HR. gnaw mm MM NH §M> N<SM> <V~ Act .MH.

 

 

 

88m? 88385 05 we 83 Boa?

 

24

 

a:
_ 5
8258302

 

v

 

 

 

 

 

 

m
LEI
V. gm
bazssv.

 

 

 

 

 

83350 maxim 35>

 

 

 

 

 

 

 

 

 

 

 

 

w 5+3. + , _a
_> 30 4am g I g> :3

g
> Hogans.

 

 

 

 

 

 

 

 

 

 

 

 

 

mm

 

 

 

 

25

 

 

TRUI = -U1+V,-, (2.3.1.1)
TFU3 = —U3+ K2,” - an“: ”(5) , (2.3.1b)
TA VR = KA(V,¢f — U1 - U3) - VR, (2.3.1c)
TEE“ = VR - (SE +KE)EF,.. (2.3.1d)

The steady state equations of (2.3.1) are obtained by setting all derivative terms

to zero. Thus

 

 

o = Ul — Vi, _ (2.3.2a)

o = U3 - Kg," + KFEF‘g: “(5), (2.3.2b)

0 = KA(V,¢f - U1 - U3) — VR, (2.3.2c)

o = VR - (SE 4 KE)E,.-,-, (2.3.2d)
i=l,2,°-',n.

It can be shown, after elimination of the variables U 1, U 3 and VR from (2.3.2),
that the terminal voltage magnitude V; and the field voltage E p,- of the i-th generator

are constrained by the following algebraic equation

0 = (SE 4' KE)EFi +KAV“ -KAVIef’ (2.3.3)

Solving for BF,- from (2.3.3), we have

26

KA
5;; = KB Vuf - KB Vi, (where KB = m), (2.3.4)

Plugging (2.3.4) into (2.2.221), (2.2.2b), and (2.2.2c), we obtain

0 = Pi" "' bdiKB [an Vi sin(5,- ‘6‘) "' ViZSiI'Ksi -9,- )1, (2.153)

At the i-th (lSi Sn) terminal generator bus, (2.2.2b) yields

0 = g,(5 134,9 ,v,¢) (2.3.5b)

tat-KB [an V.sin<e.-—6.) - V.2sin<e.-8.-)1 +

21:313.,- v, V}. cos(e,.—ej) + by. v,- VJ. sin(9,-—8j )] +

3:331 [git Vi Vk (308(9“ _¢k) + bit Vi Vk Sin(9" 4k )] ,

and (2.2.2c) yields

0 = h,(5,£ ,e,v,¢) (2.3.5c)

12,,- (V,2 — K, an V,cos(e,-—5,-) + K, V}cos(e,--6,. )] +
2;:11Uij Vi Vj sin(9,- ’ej) - bi; Vi V} 005(9.‘ '9} )1 +
2137121 [8.]: Vi Vk SW9: 4’2) ‘ bit Vi Vk 005(9i ‘4’]; )1.

At the k-th (n +2Sk Sn +m +1) load bus, (2.2.2d) yields

27

Pkd = gk(89E :99V:¢)
= $313,, V,‘ V,cos(¢,,-e,-) + b,‘,- v,I V,- sin(<l>,,-e,- )] +
£12.32 1 [81d Vk V1 005(¢k 4’1) + bu V1: V1 Sin(¢k "4’1 )1.

and (2.2.2c) yields

de

hk (5 ’Eq 99 9V9¢)

”:11 [81d Vt Vi Sin(¢l ‘9: ) - bh' Vk Vi °°S(¢k 'ei )] +

21";"13 1 [81d Vk Vt Sin(¢k ‘01) - bu VA: VI 008“”); ‘4’: )1.

Let the complex variables 5,, if, and E, be defined by

E‘- Z= Vicos(5,-—6,-) +j Vi sin(8,--9,- ),

A

Bi 1: ViCOSGi +j Vi Sinai,

and

it I: VkCOSd’k +j Vksindh,

k'=n+2,-°-,n+m+l.

Then the static equations become

1 - —e
0 Pi". - iydiKB (Vref _ Viin - Ei )’

(2.3.Sd)

(2.3.56)

(2.3.6a)

28

associated with

0 = [fig—viz, i=1,2,--

and

0=éiéi."vi29 i=1.2."°

At the terminal generator buses, one has

1

or",

0 = —YdiKB (Vi - Vrefxii — it") +

2

l A A A A
3(1‘3: 315551. +5: ZiIYiIEI) 4'

1 e e an ..
3(5} "331%;51: + El 33.55%».

and

l

0 = —.-)’di['2Vi2 + KB (Vref ' ViXE'i + E: )1 +

21

1 . . . .
37(Ei2l'31yi'75; ‘ 35.21'Lilyi'151) +

l A ~ 0 5.
ZIP(5: ”$25352 ’ El XI:

"$21M 21: ).

(2.3.6b)

(2.3.6c)

(2.3.6d)

(2.3.6c)

29

At the load buses, one obtains

0 = ‘ “'(Ekzi'i‘il YIdEi +52 n-HMdEi ) "

'2—(Ekzm21YuEl + EkamiilYuEI)’

k=n+2,---,n+m+1,

and

0 = '51;(Ek251+1)’551 ”EkZi'SIYkiEI)‘

217(Ek2m21m51 - 51:27:32 mil).

k=n+2,-~,n+m+l.

(2.3.60

(2.3.6g)

Chapter 3

ON THE NUMBER OF (EQUILIBRIUM) STEADY STATE
SOLUTIONS OF POWER SYSTEMS

The fundamental theorem of algebraic geometry states that the number of the iso-
lated solutions of n polynomial equations in n complex unknowns is bounded above by

the total degree

d = Hin=1di ,
where d,- is the degree (i.e. the power of the highest ordered term) of the i-th polyno-

mial. This result is referred to as the Bézout theorem, and the total degree d is often
referred to as the Bézout number.

The load flow of a power system, expressed in the polynomial form (e. g. the
rectangular-coordinate form (2.1.4) or the complex form (2. 1.9)), comprise 2N polyno-
mial equations each with degree d,- = 2. Consequently, the number of solutions is
bounded by their Bézout number N

d = Hin=1di = 22”-

To get a feel for these numbers, we choose N = 10 to get 220 = 1,048,576.

However, it will be shown that the load flow models of power systems belong to
the class of deficient polynomial systems ([18], [19]) and thus the number of their

30

31

finite solutions are bounded by [21:7] . For N = 10, the bound on the number of solu-

tions is 8"] = 184,756. The difference is 220 - [3'8] = 863,820.

The basic homotopy methods ([31-[4], [35]-[38]) use 22” initial points to trace

22” homotopy curves. When applied to a deficient system such as a power system,

however, at most [213]] homotopy curves reach their finite solutions while the rest of

the homotopy curves grow unbounded as t approaches 1. From a computational point
of view, tracing the non-finite solution curves takes much more time than tracing the
finite solution curves. This would consequently represent wasted computational time,

and in some cases would cause serious difficulties in numerical computations.

In this chapter, we describe deficient systems and show that the (polynomial)
models of power systems are deficient. For the various models of power systems, we
also determine an upper bound on the number of solutions and conditions under which

this upper bound can be reached.

3.1. A Deficient System and Its Associated Homogeneous System

Suppose we are solving a system of n polynomial equations in n unknowns,

namely

F(x1.-°-.xn)=[f1(xl.---.x,.).---.f,.(x1.°--.x,.)]T=0. (3.1.1)

Let the degree of the i-th polynomial equation f,- be (1;, where i = 1, 2 , - - - , n.

Then we associate to F a homogeneous polynomial system, namely,

. .. .. T
F(xo,x1,---,x,,): = [f1(xo,x1,--°,x,,),---,f,,(xo,x1,-°-,x,,)] (3.1.2)

1‘
d d
= [xo'fro‘r/xo . ' ' ' . Xn/xo) . ‘ ‘ ° . xo'fno‘r/xo . ' ' ° . xn/xo)] °

32

According to the Bézout theory, the number of the isolated solutions, counting multi-
plicity, of a homogeneous polynomial system E (in the projective space CP") is
expected to be dlx . . .x (1,. However, the associated homogeneous system E and the

original polynomial system F are related by setting x0 = 1, that is,

.. .. ., T
F(19x19°°'9xu): = [f1(19x19'°°9xn)9°'°’fn(19x19°°.sxn)] (3°1°3)

= [f1(xl’°°°rxn)a°"9fn(xl9°H’xfl)]T = :F(xl"”’xn)°

It has been shown that the solutions of the associated homogeneous system 13' = O
in CF" consist of both the finite solutions (for x0 at O) and the non-finite solutions (for
x0 = O). The latter solutions are called "the solutions at infinity" [19]. For an example,
we consider the following 2 polynomial equations F (x 1, 1:7) in 2 variables [19],

from. 12) x12 -xr+X2-2 = 0. (3.1.4)

1112-4114-312-4 = 0.

f2(x1.x7,)

Its associated homogeneous system E (x0, 1:1, 17) is given by

It? —x0x1+xox§-2x& = 0, (3.1.5)

Nxo 11: 1?)

f2<x0’ x1,X2) x1x2-4x0x1+3x0x2-4x& = 0

The Bézout theorem actually states that the associated homogeneous system
E(xo, x1, x2) of (3.1.5) has d, x d, = 2 x 2 = 4 isolated solutions in CPZ. These solu-
tions include both the solutions for x0 at 0 (the finite solutions) and the solutions for
x0 = 0 (i.e. the solutions at infinity). Therefore, the number of the original system
solutions in C 2 (the two dimensional complex space) equals the number of the solu-

tions of the associated homogeneous system minus the number of its solutions at

33

infinity.

It can be checked by hand calculations that the associated homogeneous system
(3.1.5) has exactly 4 solutions, one of which is non-finite in terms of x0 = O (i.e. the
solution at infinity), which is basically the linear space 2 := [(x1, x2). x1 = O}, and the
remaining three are finite in terms of xo¢0 (these are the solutions of (3.1.4)).
Therefore, the number of the original system solutions in C 2 (the two dimensional
complex space) is 3 (: = 4.1).

In order to predict the number of solutions of (3.1.4) a priori, one must determine
the "solutions at infinity" of the associated homogeneous system (3.1.5). The relation-

ship between the original system (3.1.4) and the associated homogeneous system

(3.1.5) is clearly that ifxo = 1, then Feta, x1, x7) = F(x1,x2).

3.2 Nonsingular Zeros on Projective Spaces

In geometric terms for general polynomial systems [18], if Z is the common zero
set in the n-dimensional complex projective space CP" of the associated homogeneous

system E , then Z is the disjoint union
{(1,11,’ ' ° 9xn)lfi(x0vxl 9” ' ,X~)=0} U [(0.11 o” °9xn)lfi(x09xl 9” ° 9110:01-

Since E(1,x1,-~,x,,)= F(xl , - - - .19.). Z is then the disjoint union of the zeros
of F in C" and the zeros in the "hyperplane at infinity” for x0 = 0. Here the n-
dimensional projective space is denoted by

CP” :=[x=(xo,x1,-°-,x,,)lx,- e C,x¢(0,0,-~,O)]/" .

Herex" ysignifiesx=cyforsomenonzeroc e C.

Theorem ([18], [32]). Leth , - . - ,f, be polynomials in x=(x1 , - . - , x") of degrees

d1 , - - - , d" , respectively, and suppose that the common zero set of the associated

34

homogeneous system, say Z(f~1 , - ° - , f1.) in CF" is a disjoint union of nonsingular

points y, , - - - ,y, and nonsingular linear spaces Z1 , - - - , Z,. Then

dlx -- - xdn = r + 2213,], (3.2.1)

where for a linear space 2 = CP‘ , the equivalence [Z] is given by the coefficient of t‘

in the power series [18]

(1+:)“" ni;1(l+d,-t). ' (3.2.2)
As an illustration of this theorem, consider the case when d1 = . . . = d" = 2,
then
(1+2t)"(1+t)“" = (1+I+t)"(1+t)"" (3.2.3)
_ . . n
_ 2.."=0(1+r)"“t‘ [Jud-0‘“
= 2?=O(l+t)"iti[?].
_ n
Thence [Z] - 29:41.] .

Definition. We say apointy 6 2(71 , . - - ,f,) is nonsingular, if

 

rank 30%.”an

‘a(xo,...,x..)<y) = codimchtwwfm"). (32-4)

where codim denotes complex codimension with codimy (Z , CP") = n - e , and
= dimyl signifies complex dimension of Z in CF". A variety 2 is nonsingular if

each point y of the variety 2 is nonsingular [18].

35

3.3. Application to Quadratic Polynomial Systems

Consider the general quadratic polynomial system in 2N equations with 2N vari-

ables, namely

f1 = L11L12 = 0. (33.1)

fzw = £2N.1[:2N.2 = 0.

where
E51=En(x1,° °',xN)=a,-lx1+.. .+awa “fa"
isalinearforminthevariablesxl,---,x~ fori =1,2,--°,2N,and
Lj2=Lj2CYr."°.YN)=bj1Y1+---+bjwyN+Bj
is a linear form in the variables y1,---,yN forj=1,---,2N. Its associated

homogeneous system is given by

it = 1:111:12 = 0. (3.3.2)

fay = int/dim; = 0.

where

Li1=L;1(x1,- --,xN)=a,-1x1+...+awa + aixo,

36

and

in =£j2(yl 9 ° ' ' .yN) = bji)’1+-- ~+ ijYN + Bjxo.
j=1,2,---,2N.
Calculate the Jacobian matrix J(z) defined as

361,...,fm)

3(x0,x1,~-,xN,y1.'".yN).

 

With x0 = 0, one obtains

 

 

1(2, x0 = 0) = (3.3.3)
(11er + 311m (111er . .. Out/Liz b11141 . .. 1’er1
0/3142: + 521m 0th22 - . - (1an berzt - - . b2NL21
L(12IszItI.2 + BZNLMJ 02N.1L2N.2 ' ' ' 02N.NL2N.2 banmd ' ' ‘ b2NJVL2NJJ
where

Lil = Li1(x1, ' ’ ' ,xN)=a;1x1+ . . . +awx~,

and

Ln = Lj2(yl"'°vYN)=bjlyl+---+ij)’Ns
I=1.2."'.2N.

The solutions at infinity of (3.3.2) are the solutions of the following equations

obtained by ignoring all lower order degree terms of (3.3.1). Therefore the solutions at

37

infinity of (3.3.2) are the solutions of

f1 = («11le = 0"“:sz = L2N,lev.2 = 0. (3.3.4)

Let the coefficient matrices A and B be defined by

 

 

all 012 . . . am
021 022 . . . (12"
A= III III III III . (3.3.5)
50ml aZN.2 ' ' ' “ZNNJMN
and
bu bu . . . hm
hm hm . . . b2”
8: (3.3.6)
Pam b2N,2 ' ' ’ szwszva

 

 

If all the N x N minors of A and B are nonsingular, the solution sets at infinity of

(3.3.2), say Z 1 and Z 2, consist of the two linear subspaces

Z1={(0,x1,---,x,,,,y1,---,yN)e'Cl”1W |x1=...=xN=0]=CPN“,
and
Zz=[(0,x1,°--,xN,y1,~--,yN)e CPZNIy1=...=yN=0}=CPN"l.

The codim(Z, C102") would be N+1 (=2N-(N-1)). Plugging the solutions at
infinity Z1 and 22 into the Jacobian matrix (3.3.3) and removing all zero columns from

1(2), one obtains

f<zt> =

and

f(2,) =

38

allez
022L22

a1L 12
021122

(111le
‘12an

 

‘a’zNLszz am.1L2N,2 (12,”ng

bllLll
b21L21

blZLll
b22L21

511411
[3sz1

 

Ll3I.IlIL21\I.l b2N,1L2N.l b2N2L2N.l

Now define the extended matrices A and B as

>l

and

bi

p

111
0'2

 

’51
Ba

 

all
021

02111.1

bll
b21

012
022

02111.2

1’12
bzz

_BZN b2N,l b2N.2

Then the following theorem holds.

 

 

bZNNJ

 

 

aan
‘1an
I I I , (3.3.7)
“mwbmzj
me" ‘
b2NL21
I . ' (3.3.8)
bZNJVL2NJJ
9 (3.3.9)
2Nx(N+1)
(3.3.10)
7N><(N+1)

Theorem 3.3. The system (3.3.1) has exactly r isolated (complex) solutions if all

 

N x N minors of A of (3.3.5) and B of (3.3.6), and all (N+1) x (N+1) minors of the

extended matrix A of (3.3.9) and B of (3.3.10) are nonsingular, where r is given by

39
r= 22” - [Zr] " [22] =2”, ' 2115183 ' §N+1[21N]=[21y]’

and [2] denotes the equivalence of Z. '

m. Since 1(21) and J (22) are of the same form, we only need to prove
rankcJ(Z1) = codim(z,, CPZN) = (N+1). Then the same conclusion on J(z,) would
apply to J (Z 2). From the first assumption of the nonsingularity of all N x N minors
of the matrix B, there are at most (N - 1) linear forms L12 to be zero for
y $(0,0,° --,0). SupposethatthefirstiS(N—1)linearformsLnJc‘:1,2,---,i
are zero (If not, we can always reorder the polynomial equations (3.3.1) so that such a

condition is satisfied.) Therefore, the rank of I (Z 1) can be determined by checking

the rank of the following matrix
“in ai+l.1 ai+1.2 . . . “my
“n2 ai+2.1 ai+2.2 - . - “my
K= III III III III III . (3.3.11)
_°'2N 02m 021m ‘ ’ ‘ (’2ij (m—i)x(N+l)

 

 

From the assumption of nonsingularity of all (N+1) x(N+1) of A, the rank of
the matrix A- equals N+l. Consequently, Z 1 is nonsingular. Similarly, we can draw
the same conclusion that Z 2 is nonsingular. By the Intersection Theory, system (3.3.1)

has total r isolated solutions, where r is

r= 22” 4211-1221 =22” - 21:18")- 2342,”) = [213’]
(QED)

3.4. The Classical Model of Power Systems

Intersection theory implies that if there are no "solutions at infinity", the number

of solutions of the original system is exactly equal to the total degree d. And if "the

40

solutions at infinity", which consist of the intersection of all linear subspaces Z,-
(1Si St) are nonsingular, the original system has exactly r isolated solutions where r

is given in (3.2.1).

In the case of the load flow equations of power systems, it will be shown that
"the solutions at infinity" of the load flow equations exist because the load flow equa-
tions are deficient. Not all of the solutions at infinity of the load flow equations are
nonsingular, however, because a power system is not always fully connected (Note: we

say the system is fully connected if each bus is connected directly to all other buses.)

3.4.1. The all-PQ-bus load flow equations

In the PQ-bus networks, the injected complex power at each bus is specified. The
maximum number of solutions of such systems can be attained if all the conditions of
the following theorem are satisfied. We point out on the outset that the sufficient con-
ditions are generically satisfied, i.e., they are satisfied for all but a measure-zero set of

admittance values.

Theorem 3.4.1. There are exactly [ZNN] isolated complex solutions of the load flow

equations (2.1.7) for an all-PQ—bus N—node power system (excluding the slack bus), if
all k X]: (k =1,2,-- -,N) minors of the extended complex admittance matrix II

are nonsingular. Where the extended complex admittance matrix I7 is defined by

Y11 YI2 ... I’m YI.N+1
Mr )‘22 ~- Y2N Yum

Y= III III III III III. (3.4.1.1)
bYNl YN2 ' ' ' YNN YNJV+1J

 

 

Proof. Since all It x It (It = l, 2 , - - - ,N) minors of 1" are nonsingular, the solution
set at infinity is Z1 U Z2, where

41

z,:= {(0,E,,-~,EN,E; ,~-,E,;',)Io=E, ...=E~] =CPN-l,
and
z,:= ((0,E,,~-,EN,E; ,---,E,;)Io=EI =. ..=E,;]=CPN"1.
with codim(Z, CPZN) = 2N-(N-1) = N+1.
Calculating the matrices J(Z1) and 1(22) with respect to (£0,131 , ~ - - , E13)
according to (3.2.4) from the associated homogeneous System of (2.1.7) and removing

all zero columns from 1(21) and I (Z 2), one respectively obtains

0 a1 0 0
0 0 02 0
J‘(z,)= ° 0 0 II a” . , (3.4.1.2)

YMMEE YilEi MES YmEl
Y21v+152 Y2152 Y2252 -- Y2N52

 

 

. t . .
LYNNHEN YNlEN YanN " YNNENmegt/H)

and

YENHEl YilEl yz251 YZNEl
Y2Iv+152 Y2l52 Y2252 Y2N52

 

 

. E ”17 151v 3’1; 251v .. YXINEN
j = YNN+1 N .. . . .
(Z7) 0 b1 0 O . (3 4 1 3)
0 0 b2 0
L 0 o o .. bN )2NxtN+l)

where at: = ZjN-lyszj.’ and bk: = ZfilykjEj for k =1, 2 , ‘ ’ ' , N.

42

Consider the non-singularity of i (Z 7) first. Suppose that the first i rows (e. g.
E;=0, OSiS(N-1)) become zero. It can be shown that there are at most p=N—I‘-1 rows
of 112,) with bj=0 (OSjsP). In the worst case, suppose that bj (OSjSN-i-l) become
zero. Also assume that the load flow equations (2.1.7) are reordered in such a manner
that the first 1 rows and the last p rows of i(zz) are zero. We obtain the reduced
matrix by removing the first i rows and the last p rows from I(Z7_). From the assump-

tions of the theorem, the following reduced matrix

a a o e 1
rYi+1JV+lEi+l yt‘+l,lEi+1 , , , yt‘+l,i+lEi+l , , , yi+l.NEi+l
YNJVHEN YNJEN - - - YNJHEN - . - YIvNEN
0 b1 . . . 0 o o . 0 (3°4‘l°4)
_ 0 0 ' ' ’ bi+l ' ' ' 0 j (N+1)x(N+l)

 

 

is nonsingular. Therefore, white] (2,) equals N+1. Consequently, Z 2 is nonsingular.
Since I (Z 1) and I (Z 2) are of the same form, we can follow analogous procedures to
prove that Z, is nonsingular. By the Intersection Theorem (see page 32 of this thesis),

the number of isolated solutions is
[2"] were)

Thus, the proof of the theorem is completed. (Q.E.D)

3.4.1.1. Example: a fully-connected 3-bus power system structure

We now apply theorem 3.4.1 to investigate the number of solutions of the load
flow equations for an all-PQ-bus fully-connected 3-node power system as shown in

Figure 3.1. The system’s extended complex admittance matrix I7 is given by

43

.2335 82m? $33 39m 25. A .m 0.53%

 

 

 

25 82m

 

@

 

NCO

 

 

2%

a.»

 

2.»

vanm

 

 

l7 = )’11 Y12 Y13
3’21 Y22 Y23 '

The load flow equations for a general all-PQ-bus 3-node power system can be
expressed in the complex polynomial form (2.1.7) as

£10115; + thE + y'{3) - s1 = o, (3.4.1.5a)
1520315; + ygzzi; + y;3) - s2 = o, . (3.4.1.5b)
E10115, +y12152+y13)—s'{ = 0, (3.4.1.5c)
£30,115, + me, + y23) - s; = o. (3.4.1.5d)

Its associated homogeneous system is given by

510115; ”1‘25; + y13Eo) - 51133 = o, (3.4.1.6a)
520315; + figs; + y5380) - 5253 = o, (3.4.1.6b)
£10,113, + an52 + y13E0) - 5153 = 0, (3.4.1.6c)
5502,15, + m5: + MED) - 3353 = o. (3.4.1.6d)

The solutions at infinity, by definition, are the solutions of the following equations

5,0115; +y§2I5;) = o, (3.4.1.7a)
520315; ”3253) = o, (3.4.1.7b)
510,151+ ylez) = 0, (3.4.1.7c)

45

E202151+Y2252) = 0.

Define the matrices A and B by

r a

1 O
O 1

 

)’11 In,

 

3’21 Y22J

 

 

(3.4.1.7d)

(3.4.1.8)

It can be seen that if all k Xk (k = 1, 2) minors of l? are nonsingular, then all

2 x 2 minors of matrices A and B are nonsingular. Therefore, all the solutions at

infinity consist of the two linear subspaces

21 = ((0.51, 152; E}, 15;): o = 15,: 52},

and

22 = ((0,151, 52, 5;, 5;): o = E; = E; 1.

Since both codt°m(z,, CP“) and codim(Z2, CP“) are equal to 3 (= 4-1), 21 and

22 are nonsingular if rankcJ (Z ) equals 3. The singularity of 21 and 22 can be

checked by the following steps. First, taking the partial derivative with respect to

(E0, E1 , - - . , 5;) of (3.4.1.6), one obtains

30'1""..74)

 

r223551
Y2352
>335;
grab“;

 

2:15; +Yi252

o .
YllEl
YmEz

e '5 -
3(EO.Elv"'sz) 0'0

0
yitEi +2.2in
YnEl
Y2252

22151
YnEz
qul +le£2
0

22251
22252
0

 

321514?sz

(3.4. 1.9)

Plugging the solutions at infinity 21 and Z 2 into (3.4.1.9), one respectively

46

 

 

 

 

obtains
Wtw'nft)
J z ) = , (3.4.1.10a)
( ‘ a<Eo,Et.-~.Eal WM?"
0 yItEI +YIZE; . .0 . . o o
= o. o. YZIEI+XZIEZ o o
YlsEl )‘1151 )1251 8 3 I
hsEz Y2152 Y2252 ”)6
and
ad“ , . . . ,f4) I
12 = _ . . 3.4.1.10b
( 7) 3(50. El , . . . , £2) (50'31'52'0) ( )
YisEl 0 0 YilEl )3251
= )’2352 0 0 Y2152 Y2252
0 0 0 YllEl+Y1252 0
0 0 0 0 Y2151+Y2252 4rd

Since all the k x It (1932) minors of l7 are nonsingular, both rankcJ (Z 1) and
rankcJ(Zg) equal 3 (:= N+1). Therefore, we conclude that system (3.4.1.5) has 6 [: = 2‘

- [Z d - [Z 2] = [3]] isolated complex solutions.

3.4.2. The all-PV-bus load flow equations

PV buses are mostly generation buses at which the injected active power is
specified and held fixed by turbine settings. A voltage regulator holds the magnitude of
the voltage fixed at a PV bus by automatically varying the generator field excitation.
This variation anticipated to cause the generated reactive power to vary in such a way
as to bring the terminal voltage magnitude to the specified value. The load flow equa-
tions in this case are described by (2.1.8). The number of solutions are specified by the

following theorem.

47

Theorem 3.4.2. There are exactly [213,] isolated complex solutions of the load flow

equations (2.1.8) for an all-PV-bus N-node power system (excluding the slack bus), if

the following conditions hold:

(a) all kl x k, (k, =1 , - - . ,N) minors of the complex admittance matrix Y are non-

singular,
(b) all k; x k; minors of the extended matrix I7 are nonsingular,
wherek2= 1 , - - - ,N/2 foreven N, andk2=1 , - . - ,(N+1)/2foroddN, and

where the complex admittance matrix Y and the extended admittance matrix I7 are

defined, respectively, by

 

 

r a
In In Y1N
221 Y22 Yul
Y =
_)’Nl YN2 yNNJ
and
r 1
yll Y12 ... le YlN+l
hr 222 YzN )‘2.N+l
Ll’m YNz ' ' ' YNN YNN+1J

 

 

Proof. The proof is essentially the same as the proof of Theorem 3.4.1. We first give
the associated homogeneous system of (2.1.8) by

0 = EIQjN-lytjgj. + y;,N+lEO) + (3.4.2.13)

El (2,411} 115} 4' YlN+lEo) '- ”1502'

48

0 = Eunice-E; +)’2,N+lEo) +
52 (2,5021% + Y2N+lEo) " 2P253-
0 =. EN (2,010de + YNNHEO) *-
EN(ZjY-1YNjEj 4‘ WWI/+150) " 2PNEti .
0 = 5,5; - V353, (3.4.2.1b)
0 = 525; - V2253.

o = ENE); — vfieg.

It can be shown that if assumption (a) of Theorem 3.4.2 holds, then the solution
sets at infinity of (3.4.2.1) are the union of the two disjoint linear subspaces Z 1 and 22,

where these two sets are of the following forms:

z,= {(0,E,,~-,EN,E; ,~-,E,;)|0=El=...=EN} =0“,
and

22: {(0,E,,-~,EN,E; ,-~,E,;)Io=z~:; =...=1~:,;} =CPN".

It can be seen also that the codim(z, CPz”) equals N+1. In the following, we
will show that the rankcJ (Z 1) and rankcJ(Z1) are bOth equal to N+1. It then leads to
[18] the nonsingularity of the solution sets at infinity, namely, 2, and 22. The polyno-
mial equations of the load flow of the all-PV-bus power systems may appear more

49

complicated if one compares them with the load flow equations (2.1.7) of the all-PQ-
bus power systems. However, the proof of the nonsingularity of the zero sets at
infinity, namely, Z, and 22 can be carried out by following the same procedures as in

the proof of Theorem 3.4.1.

Similarly, calculating the J(z) at £0 = 0 according to (3.2.4) and removing all zero

columns from J (2,), one obtains

J‘(z,) Izod, = (3.4.2.2)
r . O ‘ . ‘
YlN+lEt 01 411151 21251 . YlNEl
Y2N+152 Y2152 02 + Y2252 Y2N52
YNN+15N YNlEN YNzEN 0N + YNNEN
o E; o .. o '
.9. 0 ’55 " .9.
L O 0 0 " EN )2Nx(N+l)

 

 

where at: = zjilyL-Ef, k =1, 2 , - - - , N as before.

Because of symmetry, it can be seen that I (Z 1) and I(Zz) are of the same form.
Therefore, we only need to check I (Z 1), then by similarity the same conclusion can be
applied to I (Z 7).

For convenience, we separate I (Z 1) into two parts. The first N rows of I (Z 1) are
called part 1, and the last N rows of I (Z 1) are called part 2. Now consider the solu-
tions at infinity (21I15,;=o,ls,:1 #0 for 1:1: 1 , - - -,t', k2=i+1,--,,N
OSi <N }. From assumption (a) of Theorem 3.4.2, it follows that there are at most (N-
H) rows of part 1 in which at become zero. Because of E; = 0 and 5,; #0

(k1=1,---,i, k2=i+1,---,N), all the rows containing only 5,; areremoved

50

from I(Zl), and all the rows and the columns that share the element E,"2 in part 2 are

removed from the matrix I (Z 1). Call the remaining matrix .7 (Z 1).

 

 

0 at 0 o
0 0 02 0
f(2,) = o _ o . o .. II a,- . (3.4.2.3)
Yi+lN+15i+l _ Yi+l.l£i+l yi'+l.ZEi+l " yi+uE,-'+1
E‘ E' E' '-
k)’N.I\I+1 N YNl N YN,2 N yNJEN JNx(i+1)

Therefore, we can simply check the rank of .7. It can be shown that the nonsingularity
of 21 is verified by showing the rank of J’ to be 1+1 (:= N+l-(N-i)).

The matrix .7121) can further be simplified by removing all zero rows from it.

There are two cases to be considered.
(1) i < N - i.
For convenient, we assume that the first i rows of I- become zeros because of (ti

= 0 (lSjSi) (this is the worst case), these zero rows are removed from the matrix .7.

The new matrix is denoted by I1: = I1(Z 1)

r \
Yi+lN+lEi:I-l yi'+l.lEl:<I-l Yi+l.25t:+l .. Yi+lr5i:+l
yi+2N+lEi+2 Yi+2.1Et+2 Yi+225i+2 .. Yi+2dEl+2
frat) = 222 III III " III . (3.4.2.4)

 

 

E. E. E. I .E.
‘ YNN+1 N YN.1 N YN,2 N YNn N J(N-i)x(t'+l)
Clearly, the number of rows are greater than or at least equal to the number of

columns of I1. By assumption (b) of the theorem, we conclude that the rankcj 1 equals

i+1.

51

(2) i 2 N-i.

From assumption (a) of the theorem, there are in this case at most (N-i-l) rows
with a; = 0. For simplicity, we suppose that the first (N-i-l) rows (this is the worst
case) become zero, and are removed from .7. The remaining matrix is denoted by .7 : =

i(z,)

~ 0 0 II I (I). II at
J(Z = a e e e . 3.4.2.5)
l) Yi+1N+lEi+l yi‘+l.lEi+l -- yt‘+lJV-iEi+l " yi+1,iEi+l (

 

 

L YNJV+1EN YNlEN " YNJV-iEN " YNJEN , (i+l)x(t'+l)
Since a, (N -i SkSi ) are nonzero, we remove all the rows and the columns that

share of from I. The remaining (N-i) x (N-i) matrix is denoted by .72. Therefore

Yl+1N+1Ei:+l yt'+1.15i:+1 Yi-o-leiE-l , , , yi+lJV-i-lEi:+l
yi'+2,N+lEi+2 yt'+2.15t'+2 Yi+225i+2 hamlet-15:4
j2(Zt) = III III III III III . (3.4.2.6)
E‘ E‘ E‘ ._ E‘
L)’N.IlI+l N YN,1 N YN,2 N INN-s 1 N .(N-i)x(N—i)

 

 

Obviously, from assumption (b) of the theorem, rankcjz equals N-i. Therefore,
rankcI(Z,) equals N+1 (i.e. (N -i )+(N -i )+(i-(N—i -1))=N+1). This infers that Z, is
nonsingular. The same conclusion can be obtained for Z 2 by taking exactly the same
procedures as above. From the Intersection Theory, system (2.1.8) has total r isolated

solutions in total. Where r is given by

restattzsomsetn- amiable].

52

This completes the proof of the theorem. (QED)

3.4.3. General load flow equations of power systems

In general, a power system consists of both PQ and PV buses. The load flow
(2.1.9) is governed by both the PQ-bus and the PV-bus load flow equations combined.
In the following, we give explicit sufficient conditions on the (complex) admittance
matrix of the load flow network which would ensure that the load flow has a precise
(complex) number of solutions. We point out on the outset that the sufficient condi-
tions are generically satisfied, i.e., they are satisfied for all but a measure-zero set of
admittance values. From a practical viewpoint, however, we may consider that the
"physi " admittance matrix for a fully-connected network satisfies the sufficient con-
ditions naturally. We can follow almost the same procedures to analyze the number of
solutions for such a system. First we rewrite the general load flow equations (2.1.9) as

follows:

At the PV buses:

0 = E,(zj”;;ly;jEJ-‘) + Enzjfllyuej) — 2P,, (3.4.3.1a)

0 = E~,(2j’iily§,.j5;) + 513, (Ej'iilYNPIEI) ‘ 2PM:

o = 5,3; -V,2, (3.4.3.1b)

o = EN‘EX,‘ -V,&".

At the PQ buses:

53

0 = EN,+l(sz=ll)’N,+l.jEjI)“SN,+1t (3-4-3-10)

EN(EjN=+ 'l lyNjEIj -SN9

C
II

o = 5,;‘+,(z;"=~;lyN‘,,JEj)-s;,‘,,, (3.4.3.1d)

0 = E~< firm-E t>-St‘t.

There are 2N complex variables in 2N polynomial equations. The first 2”: equa-
tions are the PV-bus load flow equations, and the last 2(N —Ng) equations are the PQ-
bus load flow equations. The following theorem holds.

 

Theorem 3.4.3. There are exactly [219’] isolated complex solutions of the road flow

equations (3.4.3.1) for an N-node PV- and PQ-bus power system (excluding the slack
bus), if the following conditions hold:

(a) all It, x It, (k, = 1 , - - - ,N) minors of the complex admittance matrix Y are non-
singular,

(b) all [:2 x Ikz minors of the extended matrix 17 are nonsingular,

where 1:2: 1 , - - - , [N -N8/2] forevenNg, and k2: 1 , . --, [N -(Ng-1)/2} for
odng.

The m. of this theorem is given in the appendix. Note that the proof of Theorem
3.4.1 and 3.4.2 immediately follow from Theorem 3.4.3.

While Theorems 3.4.1, 3.4.2, and 3.4.3 conclude the exact number of solutions
when the stated assumptions hold, they say nothing when those assumptions are not

satisfied. The following theorem gives an upper bound on the number of solutions of a

general power 51.3mm

Theorem 3.4.-l. The number of isolated solutions of the load flow equations for an N-

 

r. A.
node power sysmm (excluding the slack bus) is bounded above by RH .

Proof. Let IKE ) denore the load flow equations of an N-node power system exclud-
ing the slack bus, where E isdefincd as E = (E, .' . . .Ey.EI . ‘ ' ~ .E;]. And let
T denore its associated homogeneous system. Define a system SKI:1 ) as follows:

P \

(XilaliEI + al.%'+l)x(£I:1bliEi. + bl..\'+1)
(Zilah’Ei + 02;\’+1)X(2,‘I:1b215: + b231,)

(At
A
("111
V
II
F
'N
b
t»
l J
v

......

N - N o
(25:102VJIE1I + 02v.v.1)xlzr=1b2.v.i£l + blow“!

 

 

where the complex constants (1,} and 1),]- (lSim, 15} SN +1) are chosen randomly.
The random choice ensures that S (E ) has [21y] distinct solutions.
Let S denote the associated homogeneous system to S . Let Z(S ) denote the solu-

tion set of the associated homogeneous system S. Similarly, let Z('I‘) denote the solu-

tion set of the associated homogeneous system T. Clearly,

(1) the points of Z(S) lying at infinity are nonsingular as a consequence to the ran-
dom choice of the coefficients (0,7) and (1),,- ), and

(2) every point of Z(S) at infinity is also a point of 2(7).

Now define a homotopy function, H (E , t),

H103“. t) . ..
H(E, t) = -;- =(1-t)cS(E)+tT(E), (3.4.3.3)
H,.(E. t)

By Theorem 2.2 in [19], for almost all choices of the complex coefficients (0U), (bit)

55

(lSiSZN , lSjSN +1) and c, every one of the smooth and disjoint [21y] solution paths

of H(E, t) = 0, parameterized by t e [0.1), beginning at the zeros of S, will either
converge to an isolated solution of T (E) = 0, or grows unbounded to a "zero at
infinity", as t approaches 1. For each isolated solution of I'(E ) = 0, there must be a
solution path of H (E, t) = 0 converging to it. This implies that the system T(E ) has

at most [213,] isolated solutions. Thus, the theorem is proved.

3.5. The Model with Internal and Terminal Buses of a Generator

For convenient, we renumber the generator buses and load buses in the following
way: the i-th generator’s internal stator-based voltage is E,- (ISi Sn ), and its terminal
voltage is denoted by E, (n +1Sl $2n ), and the k-th load-bus voltage is denoted by E,
(2n+1$lt$2n+m ). The slack bus is taken as [2n+m+l]. In this arrangement, the static
polynomial equations (2.2.3) become:

0 = Pin-Pie (3.5.13)

1 e
= P," - 315.15? ‘ Er:+i))’t'.n+i +

ENE: "’ Badman]: (.Yi.n+i =Ydt).

At the generator terminals, (2.2.3c) yields

56
o = g,(E, E‘) (3.5.1b)
= -§-[E,(E," "' EII-rtb’III—n + 51751 " El-n)yl,l-n] +
$451 EJSHYIIEI 4' El. 2"‘I""Il)’til'3l)

'art-l

l=n+1,n+2,~--,2n,

and (2.2.3d) yields

0 = ms, 5‘) (3.5.1c)

1

31715) (EII ‘ EII-nD’III-n ‘ ENE: ' El-n )le-n] +
l e a

-27(Et2€'.*£i*‘yt-Et " EIIXZIanYIi Ei)

I =n+1,n+2, - - - ,2n.

At the load bus, (2.2.3c) yields

0' = Pf - g,(£, E‘) (3.5.1d)
= P: " %(Ek2izgn+£+ly;i5i. 4' 5: EntrifimEi).
It =2n+1, . - - ,2n+m,
and (2.2.30 yields

0 = Q: - Inns, 5‘) (3.5.1c)

57
l
= Q: ' ‘2)th 25.54.5125; ‘ 5; ngIth‘Ei)’

k = 2n+1, - - - , 2n+m.

It can be seen that system (3.5.1) is of the form of (2.1.9) in 4n+2m polynomial

equations with (4n+2m) unknowns. From Theorem 3.4.4, the number of complex

solutions of the polynomial equations (3.5. 1) is bounded above by [3:33;] .

3.6. The Model Augmented by the Excitation System

For system (2.3.6) which has 5n+2m polynomial equations in (5n+2m) unknowns,
the upper bound on the number of complex solutions can be determined by the follow-

ing theorem.

Theorem 3.6. Let n denote the number of generators excluding a slack bus. Let m
denote the number of load buses. The load flow equations (2.3.6) including both flux
decay in the field winding of the rotor and the excitation system have at most

30+"! [5" 't'ZM]
l

'=2n m isolated (complex) solutions.

Proof. The theorem can be proved by the homotopy method. As the proof of
Theorem 3.4.4 in the previous section. Let T(E ) denote the static equations (2.3.6).
We define

E =[V1.“nVntitw'wEWEItwid-

Let 7' denote its associated homogeneous system. Let 2(7' ) denote the set of zeros of

the associated homogeneous system 7‘ . .

58

Construct the "initial" starting system as

[2,"=,(a,,V, + “innit + 01.2“: it)
+ 21'1““ 1,3n+jEIj) + a1]
Xlzt'"=1(briVi 4' b1,n+i E? + b1.2n+t' E?)
+ Efi'l(bl.3n+jE;) 4‘ Bl]

3(5) = ...... , (3.6.1)

[Eli-IWStthJVi + “5n+2m.ntiE-i + “5n+2m.2n+iét’)
M
+ 2j=l(aSn+2m.3n+jEj l: asn+2ml .
xl2t°"=1(bsn+2m.t Vi + b5n+?dn.n+iEi + b5n+2m.2n+iEi)
~ e
+ 2fi1(b5n+2m.3n+jEj) + BSIH-ZM]

 

 

t J

where all the complex constants (age), (by- ), ((1,), and ([3,) (1Si55n+2m, lSjS3n+m)

are chosen randomly. It can be shown that S (E) has 3“" [5" 72m]

'=2n m isolated (com-

plex) solutions.

Let Z(S) denote the set of solutions of the associated homogeneous system S = O.

Analogous to the steps in the proof of theorem 3.4.4, it can be shown that
(1) the sets of solutions at infinity of S , consisting of the two subsets
21 = {(0,E)e CP5"+2~ IO=V,0=E, , . - -,i:,,,) = crew-1
and
z2 = [(0,E)e CP5"+7”' Io=v,o=§{ ,- - 313;. } = CP’Mm‘l
are nonsingular due to randomly choosing all the complex constants (0,7), (b,,- ), (on),
and (Bi) (1SiSSn-t-2m, lSjS3n+In), and
(2) every point of the sets Z, and 22 is also a point of 2(7).

Therefore, by Theorem 2.2 in [19], we conclude that by the random choice of the

complex constants (0,,- ), (by ), ((1,), ([3,), and c, every one of the smooth and disjoint

2,3221%, [5" TM] solution paths of

59

H(E, t) = (1-t)cS(E)+ tr‘uz) = 0
emanating from the zeros of S , will either converge to an isolated solutions of T(E) =
0, or grows unbounded to a "zero at infinity", as t approaches 1. For each isolated
solution of T‘(E) = 0, there must be a solution path of I-I(E, t) = o converging to it.

Thus, the proof of the theorem is completed.

3.7. A Structure of Special Power Systems

The conditions given in Theorems 3.4.1, 3.4.2, and 3.4.3 can not be satisfied for a

system which is not fully connected. Therefore, the number of solutions of the load

flow equations for such a system may be much smaller than the number [213,] . In the

following, we will discuss special power system structures in which there is only one

node connecting two subsystems (see Figure 3.2 for example).

3.7.1. The cluster method

Theoretically speaking, an N+1-node original system can be separated into sub-
systems, and the solutions of the load flow equations of the original system are
bounded by the product of the upper bound on the number of solutions of the load
flow equations of the subsystems. As an illustration, the power system network given
in Figure 3.2(a) is separated into two parts called subsystem 1 (Figure 3.2(b)) and
subsystem 2 (Figure 3.2(c)) by re-choosing the common bus as a new reference bus-
One may determine the number of the solutions of the load flow equations of subsys-
tem 1 and subsystem 2, if possible. Then the product of the number of solutions of
the load flow equations of subsystem 1 and subsystem 2 will determine a reduced
upper bound on the number of the solutions of the load flow equations of the original
system (Figure 3.2(a)). Thus the total number of solutions of the load flow equations
of the original system is bounded above by the upper bound on the number of solu—

tions of the load flow equations of subsystem 1 times the upper bound on the number

£8.50: a mo weiouga 2F . Nd ounwwm

 

 

 

 

a Esmxmnsm .on

 

.2 53:38 .on

 

 

.8393 Refine 2n. .3

25 gum—m

 

 

 

 

 

 

61

of solutions of the load flow equations of subsystem 2. In fact, the following theorem

holds.

Theorem 3.7. The number of complex solutions of the load flow equations of the ori—

 

ginal system shown schematically in Figure 3.2(a) is bounded above by the upper
bound on the number of solutions of the load flow equations of subsystem 1 (Flgure
3.2(b)) times the upper bound on the number of solutions of the load flow equations of

subsystem 2 (Figure 3.2(c)) with the common node as a new reference bus.

m. The validity of the theorem is obvious. Assume that subsystem 1 is solvable.
For each known solution, say EC , we solve for the solutions of subsystem 2 by re-
choosing the common node as a reference bus. Then, all the solutions of the original
system are included in the combination of all the solutions of subsystem 1 with all the
solutions of subsystem 2. This implies that the number of solutions of the load flow
equations of the original system is bounded above by the upper bound on the number
of solutions of the load flow equations of subsystem 1 "times" the upper bound on the
number of solutions of the load flow equations of subsystem 2. Thus, the proof of the

theorem is completed.

The importance of Theorem 3.7 is that we can partition the original system into
as many subsystems as possible if the original system consists of subsystems which
can be separated by one (shared) node. Then the total number of solutions of the load
flow equations of the original system is bounded above by the product of the upper
bounds on the number of solutions of the load flow equations of all the subsystems.
The results are appealing from practical view points since power networks are sparsely
connected, i.e. each bus is connected to relatively few (neighboring) nodes. Note that
one may extend the applicability to multiple sub-systems whereby each two subsys-

tems are connected via only one node.

62

3.7.2. Examples of power system networks

Here are three examples to be used to illustrate the procedures and to demonstrate

their capabilities.

3.7.2.1. A not-fully-connected 3-bus network

A simple 3-bus network is shown in Figure 3.3(a). It can be seen that there are
only two lines interconnected in serious between the buses. Since two lines are shar-
ing a common bus, we can separated it into two subsystems as shown in Figure 3.3(b)
and 3.3(c). By Theorem 3.4.4, the upper bound on the number of solutions for each
subsystem is 2. Therefore, the number of solutions of the load flow equations for the

system shown in Figure 3.3(a) is bounded above by 4.

3.7.2.2. A 7-bus network

The second example of a 7-bus network is given in Figure 3.4. Since it is am

completely interconnected, by Theorem 3.4.4 we only can say that the number of com-

plex solutions of the load flow equations is bounded above by [162] = 924.

We separate the system into three parts motivated by Theorem 3.7 by taking the
shared bus as a reference bus for both subsystem 2 and subsystem 3. Now we can
easily calculate the upper bound on the number of solutions of the load flow equations
for each of them. Subsystem 1 is a 4-node system, by Theorem 3.4.4 it has at most 20
solutions. Subsystem 2 is a 3-node system in which the upper bound on the number of
solutions is 6. The last subsystem contributes at most 2 solutions. Therefore, the total
number of solutions of the load flow equations of the 7-bus network shown in Figure
3.4(a) is bounded above by 20 x 6 x 2 = 240.

63

930:9. 2.9m eouooeeootéado: 2:. .m.m uuswmm

 

 

 

.N 83¢.»QO .on

8:232 IV

®

 

 

._ 82:35 3

®
._|

 

was mom—m awn

lilo

 

 

 

 

 

 

 

.5393 REE—o 2F .3

@ @

Im.
a as

 

«x

 

 

 

 

3838: 35-5 05 we memeefita 9F ..v.m ohflwfinm

 

 

 

 

.m duodenum A3

0|.le

3552

 

.N Eamamnsm .on

 

 

 

 

 

3:29—9—

 

 

 

 

4 83%35 .Ae

 

 

9.0...

 

 

.803? 353.5 05. .3

cw a.»

 

 

 

 

 

 

65

3.7.2.3. The model with internal and terminal buses of a generator

Since the internal bus and the terminal bus at each generator bus of the model
given in Figure 2.2 are only connected by single line with the direct axis synchronous
reactance, thus we can separate the original system shown in Figure 2.2 into the n+1
subsystems in which each of the first n subsystems is a 2-bus network which has sin-
gle line connecting by the direct axis synchronous reactance yd,- (151‘ Sn ), and the
(n+1)th subsystem is an N+1-node (N=n+m) system with the slack bus as n+1.

Clearly, each of the first 11 subsystems contributes at most 2 complex solutions, and
subsystem n+1 has at most [2:] . Therefore, from Theorem 3.7, the number of com-

plex solutions of the polynomial equations (3.5.1) is bounded above by

2"[22‘13’] 4427.132")-

3.8. Summary

We use powerful results from Algebraic Geometry, and the Intersection Theory to
show that a power system of 2N (polynomial) load-flow equations in 2N complex unk-
nowns has fewer than, and in many cases only a small fraction of, the total degree
(i.e., the Bézout number) of solutions. In fact, the number of solutions of the load-

flow equations for a given general N-node power system excluding the slack bus is

bounded above by or equal to N” = [213,] . We also give sufficient conditions in which

this bound can be reached. Loosely speaking, these sufficient conditions mean that the
power network is fully connected (i.e., every node is connected to all other nodes).

Many practical power networks, however, are not fully connected, and in fact
connected only via a common node or a branch. We have also stated a theorem that is
the foundation of a cluster method to determine a "tighter" upper bound on the number

of solutions for power systems with special su'uctures or topologies.

The theorems and the method developed in this paper are applicable to load flow

66

models that are augmented by the excitation system as well; and thus can be used in
investigating the voltage collapse problem due to the load flow/or the excitation sys-

tem.

Chapter 4

THE SPECIAL HOMOTOPY METHOD

We have presented models for the steady state equations of power systems in
polynomial form, and developed the theorems and the cluster method to predict the
number of solutions of the static equations of power systems. In light of the multipli-
city of solutions, one is interested in obtaining all possible solutions to the steady state
equations. Practically speaking, only stable solutions are of interest. However, all or at
least the so-called type-one (or index-one) solutions are desired in the study of (tran-

sient) stability ([6]-[10], [33]).

Traditional methods are not capable of solving for all the steady state solutions.
The most common and acceptable numerical method is the Newton method and its
variants. The Newton-er method is also known to fail in the finding of the all possi-
ble solutions to the load flow problems, see. [25], for instance. A good initial guess is
very helpful to the success of this (traditional) method.

Solving all solutions of a power system of polynomial equations has been almost
impossible until the advent of the globally convergent probability-one basic homotopy
method ([19], [43]). The basic homotopy method has been proven to be superior to
Newton-like (traditional) method in systematical implementation of the numerical com-
putations both in the sequential and parallel machines to take advantage of massively
parallel computers. The basic homotopy method is also globally convergent in the

sense that one can arbitrarily choose any set of initial guesses, this method can

67

68

converge to all possible solutions from the starting initial guesses with probability one.
Here we briefly describe the homotopy method as applied to solving systems of poly-

nomial equations and underline some of its properties.

4.1. The Basic Homotopy Method
Solving a target system T (Z) of n polynomial equations in n complex variables

may be represented as

t \

 

 

7:1(21’ 22 . ‘ ' ' , Zn)

A T ’ ’ . . . ’

7(2): ”(2‘ 2i. 2") =0. (4.1.1)
7al:n(zlv 22 9 . . . 9 z")
L J

where Z is an n-dimensional complex vector [21 , ° - - , n]. Let the degree of 7",- be
di for 1 Si S n. It can be seen that if n = 1, then there are (11 complex roots includ-
ing multiplicity roots. In general, (4.1.1) has at most d = dlxdzx. . .xdn isolated
roots.

Consider an "initial" starting system S (Z) of polynomial equations,

.. all?! "br
S(Z) = --- = 0, (4.1.2)
d
082,,“ ’bn
where 01, 02 , ' - - , a, and b1, b2 , - - ° , b, are 2n complex constants. Each polyno-

mial, ang‘ - b; = 0, has exaCtly d,- distinct complex roots and can be easily obtained.

It can be shown that S (Z) has d = Hgld; distinct complex roots.
Define a homotopy function, H (Z , t),

H1(Zv ‘)

H(Z,t)= =(1—t)§(Z)+tT(Z), (4.1.3)
H..(Z.t)

69

where t is a real parameter varying from 0 to 1. Clearly, we have

ML 0) =s“(2).

and
H(Z, 1) = T‘(Z).

In [24], it is shown that for almost every choice of the constants a,- and b,- , for
1 51' s n, T‘(Z) = 0 can be solved by tracing the solution curves of the following
equation H(Z, t) = O.

Pictorially, we may think of having d known distinct roots at t = 0 (based on
(4.1.2)). As t gradually increases from O to 1, the value of each root will be changed
according to (4.1.3) and a curve based on the changing values is formed. When t = l,

the final value is one of the desired roots. These curves are called homotopy curves.

To solve a target system of polynomials using the homotopy method, a numerical
approach is used to trace all the homotopy curves. The following pr0perties of the

homotopy curve delineated in [43] are essential in tracing the curve.

(P1) The solution set of H (Z, t) = 0 has (1 disjoint smooth homotopy curves for

t 6 [0,1). In other words, from each root of the "initial" starting system

H(z,0)=§(2)=o,

there emanates a smooth curve.
(P2) The homotopy curve is always forward, i.e., as t increases, the length of the
curve increases correspondingly.

(P3) The homotopy curve does not stop or cease in midway, i.e. it is well-defined over

t 6 [0,1).

70

.3250 390:5: 23. 4.? oaswmnm

 

A u H a 25 833

28080: a 888 8a 305 :05 .882 a BE AN; 3 5
.3328 5 mos r502. 023 3088.5 3

65:28 .5095 m_ PEG .3985: 5
.Ai_.NV=vouANK we muons—om 05 93 858 can 2F A3
.Suefiiv o u S m Co 80:28 2: oz 358 53 2:. A:

 

 

 

 

 

”diatom
M l K) \I) o
850m
use a 823 380:5: 853
. can u
3.12.5: 8 £35
bani cm 823 @8080: A93 (
1 N

 

 

ANY? + ANV mod n 9.5m

 

 

71

(P4) If T(Z) has exactly p roots including multiplicity, then there are exactly p end

points of these curves at t = 1.

In the degenerate case, however, the number of r00ts is less than d. In this case,

the rest of the curves will go to infinity as t approaches 1 (see Figure 4.1).

It may seem that the homotopy method is inefficient for solving large scale poly-
nomial systems. However, the robustness, stability, and accuracy properties of the
homotopy method make it the only and best choice for solving polynomial and certain

nonlinear systems.

4.1.1. Applications to power systems

Two power systems examples are presented here to demonstrate that all possible
solutions of the load flow equations can be computed simultaneously using the basic

homotopy method.

We now describe the application of the basic homotopy method. Let the "initial"

starting system S (E ) be constructed as

 

 

01512 ’ b1
$‘(”) - “Nag-b” (414a)
“AH-18:2 "' bN+l . . O
_ “with?2 " b21v .
Obviously, (4.1.43) has d = 2?” distinct complex roots.
Define the homotopy function
H(E, t) = (1-t) 3(3) + t 1““(2), (4.1.4b)

where E = [El , - ° ~ , EMEI , - - - , Em is the (complex) voltage of the given load

flow equations in the complex form of the classical model. Consequently, to solve for

72

all the solutions of the load flow equations flE) = O, a numerical approach is used to

trace all 22” homotopy curves. The following properties of the homotopy curves are

essential in successful curve-uacing. For "almost all" (a, b) e C 2"’xC 2” with a =

(01,---,am)andb=(b1,---,b;w),

(a) the solution set of H (E , t) = 0 consists of smooth l-d manifolds parameterized
by t e [0, 1), and

(b) each isolated solution of T(E ) = 0 is reached by a finite path emanating from a

solution of S (E ) = 0.

4.1.1.1. A 3-bus numerical example

For a 3-node example showed in Figure 4.2, the load flow equations T(E ) for a
3-node system is of the 4 polynomial equations (bus 3 as slack bus) in 4 complex vari-
ables. In the complex form, I” (E) is expressed by

-

51:23:10553) + 5:23:105 El) - 2pm

,~ k=12

1(5) = . (4.1.5)

5115; — V,2

E E“ - V2
2 2 2

L .1

 

 

Consider the "initial" starting system S (E ) as

 

 

01312-191

A ~ - 02522 -b2 .-

S( )— 033:2 -b3 —0, (4.1.6)
52,522 -b4J

where a1, 02, a3, a4 and b1, b2, b3, b4. are complex constants which are chosen at

random.

73

.2955 39m 38258.35 25. .N..—u ouzwwm

 

 

 

o no. a
as coca

   

c4n~>

60568.5 5 8:: Exams—mag Ho «33 ed M 8m

 

 

o. 5.5.0

 

 

3 u _>
3. u an

 

 

74

Obviously, S (E) has d = 24 = 16 distinct complex roots. Each quadratic polyno-

mial of the form, a, 13,2 - b,- = o, contributes 2 distinct easily obtained complex roots.

Defining the homotopy function

Htié. t) = (I—t)S“(E) + tT‘tE).

we have executed our computation by tracing 16 homotopy curves. For each of the 24
solutions of S (E ) = 0 acting as an initial point, we use the Newton-Raphson iterative
method to compute the solution in each time interval. And thus tracing a homotopy
curve from t = 0 to t = l, we obtain 6 finite solutions. Table 4.1 lists the 6 (finite and

isolated) solutions of T(E ) = O for a 3-bus example shown in Figure 4.2.

4.1.1.2. A S-bus power system network

For a S-bus network depicted in Figure 4.3, the "initial" starting system S (E ) is

 

 

given by
alElz “b1
§(') — 0‘5} 4’4 (417)
- 055:2 ‘bs ’ . ° °
use“? - b3]
where a = (a;) and b = (b5). i = l, . . . ,8, are randomly chosen complex constants.

Clearly, s“ (E) has d = 28 = 256 distinct complex roots (each quadratic polynomial of
the form, aiEiz -b,- = 0, conuibuting 2 distinct easily obtained complex roots).
The load flow equations 1°(E) of a S-bus network, expressed in the complex

form, is the following 8 polynomial equations in 8 complex variables (bus 5 as a slack
bus)

 

 

 

 

 

 

75

 

 

 

 

oooooo6m+ooooooé 866666636666; Naooo6fiooo¢oo6 c
oooooo6m+ooooooé Nooooo6rwwaaoo6. Noooo6m+mcaaoo6 . m
6666666186666; 3866.9..38664 wooooo6m+wm~oooAt v
66666661666666; :ooow6mtm566m6- 3663613336. m
86666.6?66884 8836336036- 383.2-366606. N
oooooo.9+oooooo.~ 966669663666- mooooo6m+vooooQT _ ‘

m 3m . . N Sm _ 3m mfifififi

 

FE, u c mow—38> x3688 2E.

 

 

 

@509 RES 6: 29:88 39m 2: we ESE—Em 05. .So 3an

 

76

3.838: 25% 2F .m..v 0.53%

 

 

 

636966
663.6366
636.6166
6nn66_+66
63661.66

 

 

a5 2.8.» 9
852:3

 

068 gm

 

 

 

 

 

 

 
     
    

 

 

 

 

 

 

 

A83 <>2 666 a no x833: 25-0 560m
.Em A3585 5 85. £32555 .8 £3 3 u _ > @ 3 .... so
. N :v . H :—
gong: 6.2.73.9 9 fl
suprae €57,856
€24,856 sagas
gorse _ _ 3.38.9 _ o moo."
moduao w . u w ®
tons; ,6 @ 2o 80 ®

 

as:

 

3o "8; 8W4 @

582

 

 

 

 

77

l'3'1231-11’111‘5+ 5125=1Y11Er2PG1

. ~ 54:15:11’41'51' +E4ES=1Y4151 '21’04
T(E) = 5151 —5V12 . (4.1.8)

5221442151-5221=1Y2151-21Q02

54:15:041'51' E4Z§=1Y41°E121QD4

 

 

To solve the polynomial equations (4.1.8) using the hom0topy method,

H(E‘, 1) = (1-1)$"(E) + 11"(E),

a numerical approach is used to trace all 256 homotopy curves. The 54 homotopy
curves converge to the finite solutions in which the 10 solutions are the system solu-

tions of f0?) = 0 listed in the Table 4.2.

4.2. The Special Homotopy Method

It has been shown that the load flow equations of a power system are deficiem.

That is, the number of solutions can turn out to be a lot smaller than the Bézout

number. When using the "basic homotopy method" to solve for the solutions to the

polynomial load flow equations of power system, we need to trace all the total degree

d (the Bézout number) homotopy curves in ‘order to obtain all possible isolated solu-

tions. This would represent wasted computational time. We will present the special
homotopy method to reduce the number of homotopy curves to be traced.

4.2.1. The classical model

Since the number of solutions to the load flow equations described by (2.1.7),

(2. 1.8), and (2. 1.9) is at most [219'] , by the "basic homotopy method" we need to trace

all 22" homotopy curves. There are at most [219,] homOtopy curves which converge to

solutions, and the rest of the homotopy curves (at least ZM-[zlflvb will go to

78

 

 

 

 

 

 

 

 

 

 

 

0889915889 — mowwmodmbmmw 8.: magmeémgmho NN SSdfiwmowEd 3889-32.35 A:
888618089 fl ESSdfianSd. 233.21%an ~61 833.2-33861 3 Kemdflmwcwat? a
888.2%. _ awcm 5.919.836 cawvoédwmhocd 833.918 _ 58.9 985913336 w
888618998; Smwcndmdgawvd. wwogfiom1wawe561 S 28.9% Gang? SoSwdTwowcamd1 5
8886389664 guahadmémwmnod 21mm — ficmévmevfic gm gmodméoga—od 88891350ch 0
888.9%; SvSndfiggqq 83 3.9153 5.6. waowm.c_1ua§m.c1 gnmcdfihvovnbd. m
Socoodmaéooooo. _ 388.212 38.6 9.58.918386 $32 635336 Somwmdm- Gummmd V
888.9%; 333.9133? .c cagmodvn 19.86 hwocwcdfinfighfic avokudfio— 336 m
oooooodtgoocé §b§.o_+wa_§.c1 cvSvmd: Sand? whammmdfinwhgfi .c. n "N — wfiovnnuvncd. N
cocooodioooooo. — boobsdvnwbogd 38 8.91% 33.6 SwahodTa—cvwad mammodvbnmaad g
m mflm V mam M mflm N mflm g mam

 

 

 

 

 

 

21.7 snow—38> 5388 05.

 

 

 

@58 3:5 e3 e852. ”BA .5 do 88.28 2:. .3. 2.3

 

79

infinity [3].
It has been shown, by the proof of Theorem 3.4.4, that the special homotopy

method only needs to trace [21y] homotopy curves, and all possible solutions of the

load flow equations can be obtained by following [w] roots of the "initial" starting

system S“ (E‘ ) which is given by

 

 

(2.501151 4‘ Gum)
’ j E” ‘ X(Z}:1b11'51" + bum)
.1(-) (2.1102151 4' “LA/+1)
52(8) ><(Z1'A-'.=1b21'51" + b2.N+1)
$03) = "' = --- , (4.2.1.1)
.. I? . . .
152M )1 (215102111151: awn/+1)
_><(2.£1b21v,1'51' + b2NJV+1)J

 

 

where the complex matrices a = (01,-) and b = (by) (lSiSZN, lSjSN-H) are chosen
randomly. Since the "initial" starting system S (E) is trivial, we can easily calculate
the solutions by setting the first half of the N 85’s to be zero to solve for El, 52 , . ..
, EN and setting the second half of the N Sj’s (iati) to be zero to solve for E},
E; , - - - , 5);. Then all the solutions of S (171') = 0 are obtained by the combinations of
all the values of (El, 52 , ~ - - , EN) with 1111 the values 01105;, E; ,---,E,;). It can
be shown that S (E) = O has exactly [219’] roots. Thus, the number of homotopy

curves that need to be traced is [213,].

Let TU?) be the steady state polynomial equations of power systems excluding

the model with excitation system. We define the homotopy function

IKE, 1) = (1-t)c§(i§) + (fat). (4.2.1.2)

80

To solve for the polynomial equations (4.2.1.2), a numerical approach is used to
trace the homotopy curves.

Let s,- = t,- — t;_1 be stepping distance from :54 to the next sampling point ti. Ini-
tially, we have to. The value of E (to), which defines the initial point of the curve, is a
root of 8(5) = 0. Now we want to solve 5'01) for 11 = to + s1. In general, given
§(t,-_1), we have to solve for 5(1)) for t,- = :54 + 3,. For si sufficiently small, the root
of the homotopy function at IQ is near the root of the homotopy function at ti. Thus,
to solve for the root at ti, we choose the root at time :14 as initial guess. The collec-
tion of the roots traces a homotopy curve emanating from a root of S (E‘ ) = O. The
Newton-Raphson iterative method fits well in executing the "local" computation at
every 1;.

In a serial computer, the stepping distance is usually chosen to be a constant (i.e.,
s,- = s for all i). If the value of s is large, two undersirable facts could happen. One is
that the Newton-Raphson method may not converge; the other is the possibility of
curve-merging. Consequently, some roots may be missed. Therefore, it is common to
begin with a small or a conservative choice of a fixed stepping distance. Either
increasing the number of sampling points due to a small s or increasing the possibility

of repetitive computation due to a large s, more computation time is wasted.

In order to speedup the computation, one may wish to consider a greater stepping
distance. To expedite the finding of a homotopy curve, we use dynamic stepping
interval. The largest allowable value of s is 0.2. When divergence or curve merging is
detected, the stepping size decreases automatically until convergence is achieved or
curve-merging is prevented.

We have executed our computation using the approaches discussed. For each of
the 22” solutions of S (E ) = 0 acting as an initial point, we use the Newton-Raphson

iterative method to compute the solution in the each time interval. And thus tracing a

81

homotopy curve from t = 0 to 1' =1.
Specializing the deficient polynomials for a 3-bus example shown in Figure 4.2,

we choose 3 (E ) as follows:

(01151 + “1252 4’ “13)(b115l+ b1252+ b13)
(112,15 1 + 1122132 + 112901215; + one; + b”) 4 2 1 3)
(031151+ 03252 + 033)(b315l + b32552 + baa) ' ( ° ' .
5(04151 4' 04252 + a43)(b415l + b11252 + 1’43)‘

l

U!)
A
v

II

 

 

The matrices a = (ay) and b = (by) (151' S4, lSjSB) are the randomly chosen
complex constants.

Since S (E ) of (4.2.1.3) has totally 6 isolated solutions, we need to trace 6 homo-
topy curves according to (4.2.1.2) with t increasing from t=0 to t=l. All the 6 solu-
tions are obtained by tracing the 6 homotopy curves beginning at the solutions of S (I?)
= 0. These 6 solutions are listed in Table 4.3.

We redo the 5-bus network shown in Figure 4.3. Let S (E ) be constructed as

(21:1“11'31' + 015)
x<zi4=1bliEi. 4' b 15)
(22:10in + 025)
x(z,;,b2,-i,’ + 122,) (4 2 1 4)

U»

A

V
II

(2:14:10 8.1' 51' + 085)
“254:? 8.1' 51" + has)

 

 

It has been shown that for almost all choices of the complex constants (ay ), (by)
(151' $8, 1.955), and c each isolated solution of f(i) = o in (4.1.3) is obtained by

tracing a finite path of

mi, 1) = (1-1)c.§(i)+ 11"(i)

82

 

 

 

 

 

 

 

oooooo.o_+oooooo._ ._ Sooodfigmoood mooooodvowoaaod o
oooocodm+oooooog 9886133361 338.2133an m
oooooodiooooooé neocodTwNmood 3886133361 v
88861888; oamecw.o_1c:8m.o1 Vomwowdfivogafi? n
8886180084 “Seowdiamooomd. Soowflofihvooono. N
88861888; choooodzwmoood. vmoooodToowoooA- _

m 3m m mam _ mam fififiwfl

 

 

 

 

 

AD, .1. a momma—3 x0388 2E.

 

 

353 3:5 8 29:88 £51m 05 me 98:28 23. .mé 035p.

 

83

 

 

 

 

 

 

 

 

 

 

 

88.9.5868; 8:89-88 86 3859-8886 2 52.9-v8wg6 588948536 A:
688618864 $389-$886- ~mcm89$vmmm _6- 8869-3386- wmgew9-wm6m36- a
88.984 8889-8886 3839-82.86 8839-63566 _8~m.9-3mona6 w
88.9.5868. — 6588.98836- 858 _ 9-88:6. 8589-8886. 6589-8886- 5
886.98; K836986686 3829-3836 5389-66886 68689-3886 0
88.92.6868.“ 88398886- 35319-3836- 8889-89.86. 8839-863156- m
88.9.7868; 8889-36366 $589-N8N86 58899-8336 26809-8386 .v
8686.9...8884 3389-8856 8889-6386 388.9-Svhfi6 949489-3886 m
8898.— 858.9+8N86- 35309-3886- Nanwmm9-8526- 855.915.5886- N
88.98. _ 53589-30866 $88.96 ~ 886 8889-68366 8889-96386 _
m 35 #4 25 m 25 N 25 _ mam

 

 

 

 

 

 

AB." snow—33> onEOo 2F

 

 

 

A353 395 63 9838: man-m 05 we $8638 23. .vé 03mg.

 

84

emanating from a solution of f (E) = 0. The 54 finite paths converge to the complex
solutions, and the 10 computed system solutions- of the S-bus network are listed in
Table 4.4. Observe the slight numerical difference in the values of the solutions in

Table 4.4 as compared to the values of the solutions in Table 4.2.

From the computational point of view, the basic homotopy method needs to trace
22" homotopy curves to obtain all possible solutions of the load flow equations for an
N mode power systems excluding the slack bus. The amount of computational effort is
exponentially increasing with the size of a system. Consequently, the basic homotopy

method is only limited to a smaller-sized system.

Since the special homotopy method reduces the number of homotopy curves from
22" to [210’] , it makes possible 'to solve for the load flow solutions of a 7 -bus network
shown in Figure 4.4.

The load flow equations for a 7-bus network with bus 7 as a slack bus expressed

in the complex form are given by
r 7 e t ‘
ElEi=IYIiEi ‘ 51

~ ~ 5 627-1? a3; - 5 6
T(E ) = - . (4.2.1.5)
EiELth-‘Ei ' 5;

 

 

_E;27=1)'6i5t ' 5;,

where E; (151' $6) denotes the complex voltage at i-th bus, and St = P,- + jQ; denotes
the complex power injected into the i-th bus. The superscript * denotes the complex

conjugate.

85

.2833: man-h 2F 4:» Ounwwm

 

 

 

 

     
 

 

 

 

 

 

 

 

 

 

 

 

 

283 <>282 a 8 222
8 285852 2 3.2 5232.552 .8 use 8
2.5" o 58.5" 30 5.5" 0
255 u am 28 u an \H/ 225 u 82
28.255 2 m w m
55.255 5 \ 35.8... C 2.5.38.5
28.255 W
88.255 5
28.255 5
55.255 N . . 52.258.
88.255 _ E 5.58 a 58.5 ...-25.5
1.0 259» o.
85223 88 an R 5%.:
\
M «2.2.58.5 KWK 3.2.55.5 «WK 3.5.5.85 N.
mfifiw ”Ergo £55" 80

 

 

86

 

 

 

 

 

 

 

 

 

 

 

 

888255855 588225883 858525883 888518885 5 25
3535532555 328.258m85 855325535 535.52-385.55 c 25
8555555555 83525855 558.258.5555 2.585.25nfi85 m Sm
55225285535 2:2.25885 «5:22-885... 8558528885 v 25
5552-285; 5582258825 25522835 8R22§885 m 3m
”58.2-2555 5533255535 58_8.2-§9n5 888253585 m 25
2.8325558... 535258535 5855255558.? 23282555555 _ 3m
v m N _
muons—om no “38:: 2F

 

 

 

A E," 5892? 229:8 23.

 

 

@205 222 551833: 5:. 25 5o 8828 2: .3 222

 

87

Consider the "initial" system

‘

(26:10 li El 4" al7)
X<Es=lbliEi. + bu)
$650215: + 027)
X(26=1b255i' + 1’27) .

- - ° (4.2.1.6)

VJ)
A
R}:
V
II

(2.-21012.15: + ‘112.7)
.
~><(2;6=1b12.i5i 4’ blz7)J

 

 

Using the special homotopy method to solve for NE) = O, we follow the homo-
topy curves emanating from the roots of the trivial system § (5‘) to the desired solu-
tions of H(E‘, 1) = 0 according to (4.2.1.2). All the 288 complex solutions including 4
system solutions are obtained by tracing total 924 homotopy curves. Those 4 system

solutions are listed in the Table 4.5.

4.2.2. The model with internal and terminal buses of a generator

In this model, the steady state equations (3.5.1) include the effect of the direct
axis synchronous reactances which do not appear in the classical model. Let N be
defined, at this time, as N = 2n+m. Then to solve for the system (3.5.1), we can use
the same so?) of (4.2.1.1) as the "initial" starting system. It has been shown that for a
random choice of complex constants a = (05,-) and b = (by) (lSiSZN, ISjSN+1), all

possible solutions of (3.5.1) will be obtained by tracing [213,] homotopy curves accord-

ing to (4.2.1.2) which emanate from the solutions of 5 (E ) = 0. Here a 4—bus example
shown in Figure 4.5 is used to demonstrate the procedures. Where bus 1 is considered
as the internal bus of the generator, and bus 2 is considered as the terminal bus of the
generator. Bus 4 is taken as a slack bus, and bus 3 is treated as load bus. The admit-

tance between bus 1 and bus 2 signifies the direct axis synchronous admittance.

88

63:58 25% 23. .mé 0.53%

 

 

 

oVo:

      
  

     

2.210

805mg: 5 50:: 822855.: 05 .8 853
85- ... a:

 

2.5.2.85

£5.32

55 u a

 

 

5.: u _>
«5 u 5:

 

 

89

 

 

 

 

 

55.2.1595. _ whit ~cdmtnm¢cmod mamamodm+ ~©mm ~ m5 ¢ww ~8.0m+av5§. —t .V
508.2%. ~ wan 25.516?” wwbd ggmcdm+b©v¢ and wvevmcdmtwacamdt M
9595519595. ~ amnv ~55TNMM© ~55 Nvmv ~Odmt©bav~m5 hn¢mmcdm+¢efiv§d N
55.51955. g 5 235.95ng m5 ©¢m§5¢N bead N fimhvadm... ~ hm ”8.5 m
5:828 me
Q mam— m mam N mzm ~ mam 62:5: 05.

 

 

 

 

 

AB, u snow—33> 229:8 2.5.

 

 

 

35:: EB: omv £9538 $5-5 on: w: 28:38 05. 6.: 035,—.

 

90

The steady state equations of a 4-node network depicted in Figure 4.5 in the poly-

nomial form can be expressed as

f . . . 1
512az=101i5i ) + El£i2=101iEi) ' 2pm
515; - Vl

- ~ 5224-10'21'51'.)' 52

T E = g - a , 4.2.2.1

( ) EZZitlo'ZiEi)'52 ( )
5321;2031'51') " S3

_ 5325;203:190 - 53 j

where S,‘ = P, + )0, denote the complex injected power at node k (k =2,3).

 

 

Choose an "initial" starting system § (5' ) as
(2110155: 4’ 014)
x(£.'3=1b 135: + bu)
3
X(Xi3=1b2i Ea" 4' b24)

0))
A
[1'11
v
II

(4.2.2.2)

(23:106i5i + 064)
‘X(X:3=1bsi5i. 4' bu!

where the complex constants a = (03,-) and b = (by) (lSiSG, IS] 54) are chosen at ran-

 

 

dom.
Define the homotopy function as
mi, t) = (1—:)c§(£‘) + (THE). (4.2.2.3)
Finally, the 12 complex finite solutions are obtained at t=1 by tracing 20 solution

paths of H (E , t) from t = 0. We obtain the 4 computed system solutions listed in
Table 4.6.

4.2.3. The model augmented by the excitation system

The model, as described above, is completed by adding the excitation system in

it. Thus, the steady state polynomial representation 1° (5‘) includes additional two kinds

91

.8893 8:365 :a .23 x838: 25-: 2:20.? ouzwmnm

 

 

 

 

 

 

 

 

94 no.5- 53.0 on
u8> mv— mm <V~

83mm» 8:385 2: we as:

 

owe

 

 

A8525 5 8:: Sagas. 2: .8 .23

 

3:53

o.~:..c.a

3" No a

o5 u a 3.3.:

3 u _>
No u 5:

 

“$33

2.9" na
8.? u n:

 

 

92

of terms in which one term is contributed by the direct axis synchronous susceptances

in generators and the other term is created by the excitation systems.

It has been shown, in the proof of Theorem 3.6, that the solution paths of
H (E , t) = (l-t )c§(i') + tfli‘ ) beginning at the solutions of 5' (E ) = 0 will converge
to all possible solutions of THE" ) = 0. For a 4-bus example shown in Figure 4.6, let E 1
denote E 1. Let E 2 denote E‘ l, and let E 3 denote E." 1. Then the steady state equations

fli‘ ) including the excitation system are given by

p

y12K8(Vl - refXEl "'EI) - 2PGI

£15; - V12
15215; - V12
y21K8(Vref ‘ V0051" 5;) 4'
To?) = 5,22,03,12,.34- E; 2202,12,.) = 0, (4.2.3.1)

Y21[2V12 - Kan/.4 - Vow. + ED] +
5224.2038.‘> - Bantams.)
5325:203i5i.) - 5 3
( 5325;203:530 ‘5; J

 

 

where S 3 = P3 + jQ3 denote the complex injected power at node 3.

In fact, the polynomial equations fa?) of (4.2.3.1) represent the load flow equa-
tions of a 3-bus power system network by choosing bus 1 as a generator bus including
both excitation system and the the direct axis synchronous susceptance between the
terminals bus and internal bus of the generator. According to the proof of the

theorem, we can choose an starting system in the following:

f (23:01:15: + 014V1 + “15)1
X(253=1b1,-E; + b14V1 + bls)
(25110235,- + auV1+ 025)
42,111)“; + b24V1 + b1.)

0))
A
Rn
V
II

(4.2.3.2)

(2532107553 + 074V1 + 075)
f<23=1b7i 5: + b74V1 + b75>J

 

 

93

 

 

 

 

 

 

 

 

 

 

 

 

 

66666662666666; 856663. 59.66 avwnwm6m$hchmo¢ mmmeQOfimoSSA- w
66666661666666.“ N5686v3vmvo6 mnmwvm.6_-§mmoo6. anmoo.6_+a_m636 h
66666666+6666666 Eva—626669.66 ~w6w—m.6_-vman 36. 6n>n666¢~a§a6 0
66666661666666; ovmm 86723—66 Sawwc6finv336 mnhws.6_+mcw§6 m
66666661666666; mane—6.9666566 vasoo6m+av$a666 unwwm6.6_¢:m66._- .94
66666661686664 Svm 8.262366 macaqomdncnocé masqofinwvaoo. T m
6666666_+66686._ 3662.9-3226 6R6nm.6_-n6§666. 336667328. 7 N
6666666166686; 6868672386 66:666mi8m86 ww~¢36¥©$v86 N
am mm mm _m ”Magma.

 

 

NE, u c moms—3 53an on...

 

 

 

was: was as
88?? 85885 SW 5:» x838: £544 26 60 $8638 23. .N.. 64 03:6.

94

where all complex constants (0:1): and (by) (151' S7, ISj 55) are chosen randomly.

Similarly, we define the homotopy function

mi, t) = (1-:)cs‘(ié) + (IRE). (4.2.3.3)

Tracing the homotopy curves emanating from the roots of S“ (E ), we eventually

obtained 8 system solutions which are listed in Table 4.7.

4.3. Summary

Solving for all possible solutions of the load flow equations of power systems has
been impossible until the advent of the homotopy methods. In this chapter, we present
the special homotopy method that reduces the computational complexity (comparing
with the basic homotopy methods), and still guarantees in finding all the steady state
(equilibrium) solutions of various levels of detailed (nonlinear) models of power sys-

tems.

Chapter 5

THE IMBEDDING-BASED METHOD
SOLVING FOR ROOTS OF DEFICIENT SYSTEMS

The homotopy continuation methods are globally convergent, i.e., one may
choose any set of initial guesses, and by successively incrementing a parameter t, such
methods conve_rge to all solutions with probability one. However, the serious problem
of the current homotopy methods arise in their implementation. Even for a small-sized
polynomial system, say T (x ), the homotopy methods may require the tracing of many
useless homotopy curves that would not lead to solutions of the system 1“ (x ). The
number of the (useless) homotopy curves grows exponentially to the point where
efficient computation is rendered impossible even on very fast and/or parallel comput-
ers. As an example, we consider the deficient systems of the load flow equations of

power systems. The number of paths that need to be traced even by "the special
homotopy method" for a S-bus system [3] is 70 [ = [3]]. This number increases
rapidly to 924 [ = [162]] when the number of buses is 7 (i.e., a 7-bus system). This
means that we must follow all 924 paths in order to obtain the 4 system solutions.

From computational point of view, most of the computer computational effort in
tracing the homotopy curves is spent tracing the curves which do not converge! That
is, more computational effort is allocated for "useless" homotopy curves which would

not lead to a solution of the target system of polynomials T (x ). In contrast, the

95

96

computational time spent to trace convergent homotopy curves is relatively small. Con-
sequently, it is detrimental from a computational point of view to get rid of nonconver—

gent homotopy curves at any cost.

5.1. The Basic Imbedding Method

Solving a system of n polynomial equations in 11 variables may be expressed as

 

 

(A 1
7:105)
. T
T(x) = .29) = 0, (5.1.1)
Thor)
t J
whereeachfiJ =1,2,~-,n isapolynomialinx =(xl,---,x,,)ofdegreed,-.

The basic idea of the imbedding methods for computing a solution x' of (5.1.1)

consists of "embedding" system (5.1.1) into a family of maps

H(x, t) = 0, t s [0,1]; H: C" x [0,1] -) C" (5.1.2a)

defined by a homotopy or one-parameter embedding function H satisfying the two pro-
perties

H (x (O), O) 0 for a given x (O), (5.1.2b)

and

H(x(l), 1) = for") = o for all x‘. (5.1.2c)
Therefore, the homotopy curves are the solution paths that emanate from the given
vector x(0) to the unknown solution vector x(1) = x' .

Some well-known homotopy methods derived from the embedding function
(5.1.2) which guarantee finding all possible solutions with probability one are listed as

97
follows.

(1) The basic homotopy method

The number of homotopy curves that need to be traced equals the total degree d
(the Bézout number) '
d = Hit-14:.
where d; for 19' Sn denotes the degree of the i-th component T,- of the target system
T (x ). When applied to the load flow equations of power systems, we must trace 22”
homotopy curves in order to obtain all possible solutions of the load flow equations for

an N-node power system excluding the slack bus.

(2) The special homotopy method
The special homotopy method reduces the number of homotopy curves to be
traced from 22” to [2:] in finding all possible solutions of the load flow equations of

N-node power systems excluding the slack bus.

(3) The basic imbedding method [34]
Start with the system T(x) = 0 of n (nonlinear) polynomial equations
T1, T2 , - - -- , T, in 11 unknown variables x1,- --,x,, and choose an embedding

parameter I e C. The system T (x) can always be written in the following form.

for) = S(x)+F(x) 1": C" —->C" (5.1.3)
with
film)
S‘m = 5.297). (5.1.4)
fax.)

 

 

98

and

Ul(x19 ° ' . 9x";
f2(xlr°°'9xn)

F(x)= ..
quIV'WXn)

, (5.1.5)

 

 

where deg{S,-(X;)l is required to be greater than, or at least equal to, deg[f,-(x)},

i =1,2,- - - ,n. Moreover, thefunction Sbi = 1,2, - - - ,n hasonly simple zeros.

The problem of finding all the zeros of such a system has been solved by Drexler

[34] using a homotopy continuation method.

Giving a system T(x) = 0, one can choose an embedding parameter t, and con-

struct the homotopy continuation method in the following way:

H(x,t) = S(x)+tF(x), (5.1.6)

H: C" xC -)C";H(x, 1) = T‘(x).

It has been shown [34] that, for almost every curve W in the complex plane C, simple
solutions of S (x) = 0 generate (bounded, continuous, and differentiable) homotopy
curves that converge to all the solutions of the system H(x, 1) = S (x) + F(x) = 0 as t
approaches l+i0. These curves are parameterized by the complex embedding parameter
t e C, connecting 0 + 10 and 1 +i0. The justification of this method relies on relag
tively sophisticated ideas from algebraic geometry and transversality theory [34]. How-
ever, a drawback is that Sim) is at least of the same degree asfi(x); so for a heavily
deficient system, of which the number of solutions is only a small fraction of the total
degree d, only a few of the solution curves reach the system solutions and the rest of
the solutions, which do not converge to the system solutions, will go to infinity. This
would represent wasted computational time, and cause serious problems for numerical

computations.

99

5.2. The Imbedding-based Method: Practical Heuristic Approach

It is of great interest from a computational view-point to _re_dgg_e_ the number of
homotopy curves traced. It is desirable to trace a number of curves that is comparable
to the number of the actual (complex) solutions of a given polynomial system T (x ).
For deficient systems where their structures may reveal a "tight" bound on the number
of solutions, we suggest to use an initial polynomial system S (x) having simple solu-

tions comparable in number to the actual solutions of T (x) = 0. Such an drastically

 

reducing the number of homotopy curves traced.

In the following we describe the "imbedding-based" method that is based on our
experience with many simulations. Although this formulation may not have a theoreti-
cal base, we have found it to work satisfactorily in many simulations of power system
examples. Indeed, our approach renders the homotopy methods practically feasible in
the sense that it _redu_cg§ the number of traced homotopy curves to a number compar-
able to the number of actual solutions of a given (target) system T(x ). We begin by
describing the proposed imbedding based approach.

Let T (x): C " —) C " be a polynomial system. Let C[O,l) (excluding the point 1

+10) denote a set of complex curves from 0 +10 to 1 + i0. Let S‘ (x) be of the form

p

510:) syn—bl
fax) 32(1)“b2
S‘(x) = 222 = III . (5.2.1)

 

 

 

 

gnu) 573(1)-bn

b J

Assume that S (x) has only simple zeros. Consider the homotopy

H(x,t) = S(x)+tF(x), H: C" xC —> C"; (5.2.2)

100

T(x) = S(x) + F(x).

Then for random choices of complex coefficients b,- e C (ISISn) and t e C [0,1].

The zero set
{(x, t) e C" x C[O,1) lI-I(x, t) = 0]

consists of solution paths x1(t) , - - - , x" (t) emanating from the roots of H(x, O) =

S(x) at t = 0. The zeros of T(x) may be obtained by tracing solution paths to t = 1.

The boundness of solution paths xl(t) , - - - , x, (t) have been discussed by
Drexler in [34]. In the following we explain the elements of the results and their impli-

cations.

As in [34], we begin by using the elimination theory to compute a "Resultante"
function which is only a function of a single variable, e.g., x,‘ . The coefficients of such
function are finite degree polynomial of the (complex) imbedding parameter t. The
poles of the "Resultante" function occur when the leading (polynomial) coefficient of
the "Resultante" equals zero. Since the coefficient is a polynomial in the complex
parameter t, the coefficient becomes zero only at finite values of t. Consequently, for
almost all curves in the complex plane parameterized by t and connecting the points 0
+ 10 and 1 + i0, the "Resultante" function has "no poles". This follows since the zeros
of the "Resultante" function of each x,‘ is the same as the k-th "component" or projec-
tion of the zeros of the homotopy function H(x, t). The homotopy solution curves,
therefore, do not go to infinity (in the sense of some norm) for almost all curves traced

by the parameter t in the complex plane.

The approach actually says that if a strict bound on the number of complex solu-
tions of the target system T (x) is known, then all the solutions of the target system
T(x) may be obtained by tracing Nf (the number of complex zeros of F(x)) solution

paths, even with the assumption deg{S,-(x)} S deg(T,-(x)}.

101

The number of solutions of deficient systems, such as power systems, is bounded
by [2:], where N is the number of buses excluding the slack bus. This result is

obtained for the load flow equations of power systems including transmittances. Yet,
based on numerous simulations and applications of several homotopies to many power
system example models, this upper bound is usually too large. Indeed, the topology or
the figuration of a power system limits the actual number of solutions of this deficient
system. In the following, we apply the proposed approach to the same 3-bus example

without complete interconnection and the same 7-bus example of power systems.

5.3. Numerical Examples of the Load Flow Equations

The following two examples will be used to illustrate the computational pro-

cedures and to present a considerable saving in the computational effort.

5.3.1. A not-fully-connected 3-bus power system network

For a 3-bus example shown in Figure 5.1, the load flow equations in complex

form can be expressed as

q

31:23:101359 + 5:253:10551) - ”Gt
1‘ = 1.2

a.

S.

V
11

(5.3.1)
5151 " V12

1 £25; - v,2

 

 

Where I? =[E1 , - - - ,EN,EI , - - 35,311? is the complex conjugate of B.

We have proved, in chapter 3, that the system shown in Figure 5.1 has exaCtly 4
solutions. Therefore, we only need to . follow 4 solution curves by the imbedding-

based method.

Consider the polynomial equations S (E) in the following:

102

29588 39m caooacoo$==to= 2E; .m ouflwmnm

 

 

owe.“
3.. seam

  

. u >
Ageaé a 8:: sages. c .33 ed u 8m

 

o. 2.56

 

 

 

 

103

 

 

01512-171

. .- 02522 "bz

S( )= t , (5.3.2)
0351 -b3
b0452 -b4J

where a1, 02, a3, a4 and b1, b2, b3, b4 are complex constants which are chosen ran-
domly. Obviously, S (E ) has d = 22 = 4 distinct complex roots. Each of the first two
quadratic polynomial of the form 0,5,2 - b,- = 0 contributes 2 distinct easily computed
complex solutions.

Define a homotopy continuation function as

H(E, t) =(1-t)S(E)+tT(E): t e C[O,1]. (5.3.3)

To solve the polynomial equations (5.3.3), numerically, 4 solution curves need to be
traced. Let a complex parameter t be chosen as t = reje, where r is a real parame-
ter with r 6 [0,1], and 6 is a real parameter which varies from 90 to 0. Initially, we
have r=0 and 9 = 60 where 90 is randomly chosen. Let 5'1 be stepping distance
from r=0 to the next sampling point r1. Let t, be defined as

‘1 = (0+sl)ejeo(1'81) = 7181.61.

In general, given n-1, we have

ti = (ri-l + 51' )el(91-1 ‘ 903:).

Clearly, 9 will be zero when r reaches 1. So the complex function reje is a good
candidate for such a complex parameter t. In general, from the value E (tH), we need
to solve for E05) for r,- = rH 4» s,- and 9,- = 95-1 — 90s,. For 5',- sufficiently small,
the zero of the homotopy function at IN is closed to the zero of the homotopy func-
tion at t,- . The Newton-Raphson iterative method can be used to execute the "local"

computation at every ti.

104

Because of the quadratic convergence of the Newton—Raphson algorithm, there is
less motivation to use "acceleration factors" [78] in Newton-Raphson iteration at every
t,- . However, for the case of nonconvergent iteration, acceleration factor can be used to

render a divergent case as convergent. The technique is

[E](l+l) = [E](l) _ a J-1[H(E', Ml“)

at the l-th iteration. The quantity (1 is the acceleration factor and it generally lies in the
range 0.2541513. Acceleration factors below 1.00 actually slow convergence, and

acceleration factors above 1.00 will speed convergence.

In numerical computation, if the value of s is large, there is the possibility of
curve-merging. And thus some roots may be missed. Consequently, it is common to
begin with a small or a conservative choice of a fixed stepping distance. Either
increasing the number of sampling points due to a small 5 or increasing the possibility

of repetitive computation due to a large s, more computation time is wasted.

Table 5.1 lists the 4 solutions of the load flow equations of a 3-bus example

obtained as the 4 solution'curves of (5.3.3) evolved toward their values at t = 1.

5.3.2. A 7-bus power system network

For the 7-bus network depicted in Figure 4.4 of chapter 4, T (E ), expressed in
complex form, is the following 12 polynomial equations in 12 complex variables (bus
7 is taken as a slack bus) as given in (4.2.1.5).

Since the number of complex solutions of such a 7 -bus network is 288 which has
been shown by numerical simulation in Chapter 4, by the imbedding-based method, we
can follow 288 solution paths in the finding of the system solutions. Let S (E) be
defined as follows:

105

 

 

 

 

 

888.o_+888._ 888.9...888A 88881588.“ v
888818884 8.98.7 ~88dm+m884 m
888.2+888._ 8.2-88.7 ~88§+88QT m
88.9.7888; 888.98; ~88dm+3oaoadt fl

M £5 N mam g mam mflgfifl

 

 

 

 

 

A _- u swamps—Q, one—mace 2F

 

 

fleece 3:5 3 852: eoweeaaeeoeé 2: E
29:88 warm wouooccoo$=ééoa 05 80 28:38 25. Ad 2931

 

106

 

 

 

 

 

 

 

 

 

 

 

 

888.2588; 88861888; 9583:8884 8889.888; b 8m
”Emceefimweno 3832-538... 3389-393 8§q9¢8§5 c 3m
Raweoéafised Senadwnaefid amenaovfiwwfic commaeéflad m 3m
o$§.c_-wme:e.c RESeEBSé Ssnaovneged 25853333 v 25
getéégac $328883: ozefidveeflanc 3%253383 m 25
8239-3835 22853852 ~_Reo.o_.§§.o . ”Eaeeezsd m 3m
Sagerseosd 83353.32»... nflog.o_+§omo¢ osaodwtmmwo; fl 25
v m N _
80:28 .8 “36:: 0E.

 

 

 

A E," a momma; 8388 2.5.

 

 

358 3:5 was cocoa 8393835

o5 .3 #838: 25A. 2: mo £8328 2E. .N.m ofifl.

 

107

 

 

513 ‘b1
53-172
532 “b3
E42 -b4
552 "bs .
5‘03): 562-126 , (5.3.4)
5:2 'b7
15; -b8
15; -b9
.53 ’bIZJ

where all complex coefficients and parameters are randomly chosen. Since each of the
first two equations of (5.3.4) contributes 3 complex solutions, and each equation with
degree 2 contributes 2 complex solutions, thus the total number of the complex solu-
tions of S“ (1? )=0 is 32 x 25 x14 = 288.

Solving for T(E‘) = o by the imbedding-based method described above, we trace
all 288 solution paths of H(E‘, t) = O as t increases from O to 1. We obtain the 4 sys-
tem solutions which are listed in the Table 5.2.

5.4. Summary

In this chapter, we present an imbedding-based method which is shown to be
computationally efficient to calculate the system solutions of a power system. This
method can be used when an upper bound on the number of complex solutions of
deficient system is known. Consequently, one can choose an "initial" starting system
§ (x ), of which the number of solutions is near but greater than, or at least equal to,
the number of solutions of the target system To: ). We use a 3-bus and a 7 -bus power

systems to demonstrate the computational capabilities.

Chapter 6

THE HOMOTOP Y-BASED METHOD
TO DEFICIENT SYSTEMS '

In the previous chapter, we suggested the imbedding-based method in finding the
solutions of deficient systems. The imbedding-based method allows us to choose t only
in such a way that t is a curve in the plane of complex numbers connecting the points
0 + i0 and 1 + i0 continuously. Consequently, some serious problems may arise in
implementation of the algorithm and numerical computational efficiency. The follow-
ing homotopy-based method, as it will be shown, is to use real parameter t with
t 5 [0,1]. The proposed homotopy-based method still keeps the properties of the
imbedding-based method, but with a real parameter t.

6.1. The Practical Heuristic Approach

Solve a system of n polynomial equations in n unknowns

.. .. .. . T
T(x) = [71(x).T2(x).-~.T,(x)] = 0. (6.1.1)
where each i}, k=l,2,---,n is a polynomial in x =(x1,---,x,,)e C" of
degree (1;.
Let S" (x,b) be constructed as

108

109

f

3:1(11’b1)
52(x2,b7)
for»)

Xk, bk 5 C, (6.1.2)

 

 

S“.<x..b.>
.. 4,
With Sk(xk, bk) = Xk‘ - bk, 15k Sh.

Let R(x, a) be defined as

rr1(xl,al)‘
'2(12:02)
R(x,a)= jjj _ : xbakeC, (6.1.3)

 

 

{..(xn. an),
d?
with rk(xk, at) = akxk‘ , d“ = max{dk, (1%], IS]: Sn.
DefineH:C’I ><[0,1]><C" xC" —>C" by

Ph1(x, t, b, a)‘
”(1. t. b. a)

(6. 1.4)

 

 

Lh..<x. t. b. a)

(1-:)S‘(x, b) + (To) + t(1-r)R(x, a),

witha=(a1,---,a,,)e C"andb=(b1,---,b,,)e C".

Some Results:

Let NT denote the number of zeros of fix ), or an upper bound on the number of

110

zeros of Tat). Assume d, = mad“ 2 NT. For a random choices of (b, a), except

for a set of measure zero, in C", the zero set
114(0): {(x,!)e C" x [0.1) |H(x,t,a,b) = 0}

consists of disjoint analytic solution paths emanating from the roots of H(x, O) = S“ (x)
at t = O. The zeros of T (1:) may be obtained by tracing solution paths to t = 1.

The theoretical results that follow can be shown by following the same pro-
cedures delineated by Chow, Mallet-Paret, and Yorke in [43]. In the following, we will
exploit the steps in A.P. Morgan [36].

(A) O is a regular value of H(., ., b) restricted to C " x [0,1) for almost all choices of
b e C". Consequently, from the Implicit Function Theorem, it results in that
H "1(0) consists of smooth, disjoint paths.

(B) Let 1? denote the highest order homogeneous terms of H. O is a regular value of
H(., ., a) on [(x,t)lx at O] for almost all a e C". It turns out that the solution
paths of H(x,t) are bounded on t 6. (0,1).

(C) dt/ds at O on solution paths H‘1(O), where s denotes the arc length of a homotopy
curve, implies that no path can be homeomorphic to a circle. Therefore, each path
(x(s), t(s)) would eventually go to infinity or would lead to a root of H(., 1) when
t approaches 1.

The M of the nonsingularity of the Jacobian matrix 3H (J: , t )le of the homo-
topy function H (x , t) with respect to x at x, and t, needs the following Transversality
Theorem, where x, is a root ofH(x, t) at t,.

Transversality Theorem ([24], [35], [43]). Let V : R‘I , U c R'" be Open, and let F:

V x U -) R" be C' where r > maxIO, m-p}. H0 6 RF is a regular value ofF, then

for almost all a e V, except in a set of measure 0 in R", O is a regular value ofFa:

U -) R”, where Fa(x) =F(x, a).

111

In applying the Transversality Theorem above, we regard the complex space as a

real space of two dimensions. We now apply this theorem to the homotopy function H

defined in (6.1.4). Since hb‘(., ., b) with respect to b, (bT=[b1 , - - - . 6,1)
am. ., b)
h J» ., b) = —— (6.1.5)
5 db,
= -(1-t) a: O, (O<t <1,1$k Sn).
Therefore, the Jacobian matrix
M. .. b)‘
Hb(., ., b) = (6.1.6)
th. ('9 ., by

 

 

is of full rank equal to Zn. Hence, by the Transversality Theorem, 0 is a regular value
of H(., ., b) for almost all b e C".

The proof of statement (B) can be carried out as follows [24]. Given the domain

W=[(x1,---,x,,)e C" Ixi‘¢O,--°,x,-I¢O,
i1,i2,"°,ij€(l,2,°'°,n) and Xk=0

for k¢(i1,~-,ij),15j Sn].
We illustrate the proof of statement (B) for the simple case

W={(x1,°--,x,-)eC‘olxjat.0forj=l,~-~,i}.

112
H(x,r,a)=[5,(x,t,a),---,5,,(x,t,a)] (6.1.7)
with

H: C" x (0,1) x C" —> C" (61.8)

be comprised of the highest order homogeneous terms of H, that is

bk(x,t,ak) = (l—t)Sk(x)+tTk(x)+t(1-t)rk(x,ak), isk Sn,

where for 1 S k S n, S,, and 7", are the homogeneous part of St, and 7‘}, respec-

tively, consisting of all terms with degree d“. Now, since

85" - 1 a" — 1 d" 0 15'5“ 619
Ba,- ““"aa,. -:(—:)x, a: . 1 z. (..)

we conclude that the Jacobian of I? has real-rank 2i at solutions for almost all
a e C". Therefore, 0 is a regular value of H(., ., a) on W x (0,1) for almost all
a e C". From the lemma 2.3 of [43], H(xo, to) = 0 implies that x0 = 0. Hence, all
the solution paths emanating from the roots of S (x) are bounded by a constant k(to) >
O with any-to 5 [0,1). If not the case, then there exists to e [0.1) and a sequence
(1(3), t(s)) e C" x [0, to] where H(x(S), t(s)) = O and br(s)| -> oo. We may suppose

x6)

Ix (s )l
degree less than d“, and 71,, is homogeneous of degree d,‘ it follows that

->w e C" as t(s) —> re [0, to]. Since bk -b,‘ is a polynomial in x of

lx(s)|-d"[bk(x (s ), t(s )) - h,‘(x(s ), t(s ))] -> 0, (6.1.10a)

andso

113

5,.(fl :6 )) Ixts)l“'"i.(x<s). t(s)) (6.1.10b)

|x(s)|’

Ix (s)|-d"[bk(x (s), t(s )) - hk(x (S), t(s))] -) 0.

Hence, it (w, t) = 0. But 19 e 0, this contradicts lemma 2.3 [43].

Now let (x (s), t(s)) be a local parametrization of a solution path with respect to
the arc length s. We can assume [36] that dx/ds and dt/ds are not both zero at s for
which (x(s), t(s)) is defined. Then from H(x(s), t(s)) = O, we have

fifl£+§flfl

31 ds 3: ds = 0. (6.1.11)

Since elf/ax is nonsingular, dt/ds = 0 implies that dx/ds = 0. Thus, dt/ds is always
nonzero. This means that x can be parametrized by t such that x (t) consists of a path
when t increases from t=0 to t=1. The monotonicity of t implies that no path can be

homeomorphic to a circle.

6.2. Numerical Examples of Power Systems

We redo a 4-bus and a 7-bus examples to illustrate its procedures and to demon-
strate their capabilities.
For a 4-bus example shown in Figure 4.5, we have shown, fi'om chapter 3, that

the system (3.5.1) has at most

2" [11. )1

complex solutions. Therefore, we only need to trace

2"[2("+"’)] = “[3 = 12

n +m
solution paths.

For given the steady-state equations fa? ) in (4.2.2.1), we construct an "initial"

starting system as

114

 

 

E?-bl
522 ‘b2
. - E32 -b3
= .2.1
5(5) E; -b., . (6 )
52 “b5
5; "b6

where (6,.) (151' $6) are the randomly chosen complex constants. Obviously, S (E ) has
12 complex solutions because the first equation E? — bl = 0 contributes 3 solutions

and the each equation with degree 2 contributes 2 solutions.

Define a homotopy function H (E ,t) as follows:

H(E, t) = (H) S‘(E)-+ t T(E) +t(l-t)R(E). (6.2.2)

Where R (E) is given by

r w

alEl3
02522
a3E32.
04542 '
05552
.0653.

R (E) = (6.2.3)

 

 

We trace the 12 solution paths of (6.2.2) emanating from the roots of S (E) for a
random choice of complex coefficients (by) and (bi) (151' S6). We obtain the 4 system
solutions which are listed in Table 6.1. (We use acceleration factor in the range

0.29513. More details can be found in chapter 5.)

For a 7-bus example shown in Figure 4.4, the load flow equations are given in
(4.2.1.5). Construct M?) as

115

 

 

 

 

 

888.9%; vovnmodmtfioaomod nmmamc.o_+vwma—m.o finwfisdiocaadt V
9508610895. — wwvfigdtfimvwwbd bwowmodm+w hmzmd mowvmodmth ~N§dt m
9508619568. — ~vnv—o.cmtmmme ad va ~odmtnhmav>d obwmmo.c_+mmv§d N
$619595. g 82V ~o.omt wonvwod hamgdmshmgd w wmhvo.cm+w©wwg.o —

v 25 m 35 N 25 _ mam “emu—“86%.

 

 

 

 

 

AB. u gems—g 5388 BE.

 

 

 

3:89 325 NS @052: women-wmouoEo:

05 .3 29538 meme 2: mo meoufiom 2F 4.0 £an

 

116

 

 

,
E13-b1
E3 -b2
13,2 —b3
5} -b4
E} -b5
S‘(")= E} -b, , (6.2.4)
t
512 ”b7
E; -b8
E; -b9
E6'b12J

where the complex constant numbers (b;) (19' 512) are chosen at random. It can be

seen that the total number of the complex solutions of S (1? )=0 is 32 x 25 x14 = 288.

Let R (E ) be chosen as

015?
02823
03532
04542
R (E) = asE} . (6.2.5)
asEE
075:2

 

 

Since the system has 288 finite solutions, thus we only need to trace 288 solution
paths according to (6.2.2). By the homotopy-based method. We obtain the 4 system
solutions listed in Table 6.2. (Similarly, we also use acceleration factor in the range

0.2501513 to avoid divergent problem.) '

117

 

 

 

 

 

 

 

 

 

 

 

 

ooooocdioooooé oooooodiooooooé $616958; oooooodm+o8684 N. mam
§Nm¢.9-mvwwvm.o own Rodfiavwnecd Nghmodmtmwgd n: bvcdmtnvoobad O mam
mcwbvcdmtmghwcd NNEEdTZnVNfio hmwmmodmémaahhd figmadmtamuahad m mam
Own aflomtmeotvd chzedméovgvd Nanhfidmtvahmfiod hmwnacdmtgmad .V mam
8wm>ficmtng=d mwowmficmtmmamofld Emma—.2; guano awnnmmdmtvmsad m mzm
wNmNncdmtachvd cmcmncdmtebhwwnd §~8dfimw "nvnd mnomncdmthaanhod N mam
anemvmdmivovgd onmmw —.c_+vwao:.c unnwhadmivmmcocdt mthSdTfi “ego; ~ mam
v m N —
80:28 we Spas: 05. .

 

 

 

21.71. c momma—9, x2980 2:.

 

 

eased 3:5 36 8568 8.3-3932

65 .3 M838: 39$ 05 me 80:28 2F .me 2an

 

118

6.3. Summary

In this chapter, the homotopy-based method is developed and shown to be com-
putationally efficient compared to current homotopy methods. This method can also
be used when the number of solutions of a polynomial deficient system is known.
Consequently, one can choose an "initial" starting system S (x), of which the number
of solutions can be exactly equal to the actual number of solutions of the target system
f(x). We use a 4—bus and a 7-bus power system networks to demonstrate the compu-

tational capabilities and efficiency of the algorithm.

Chapter 7
CONCLUSIONS AND SUGGESTIONS

The fundamental theorem of algebraic geometry states that the number of isolated
complex solutions of 2N polynomials in 2N complex variables is bounded by the total
degree d = maid,- of polynomials. This is the statement of the Bézout theorem.
When using the "basic homotopy method" to solve for the solutions to 2N polynomi-
als, such as 2N load flow equations of power system, we need to trace all the total

degree d (the Bézout number) homotopy curves in finding all possible solutions.

7.1. Conclusions

In this thesis, we use powerful results from algebraic geometry and homotopy
methods to determine the number of (complex) solutions for various models of power
systems. We develop the theorems to predict the upper bound on the number of solu-
tions of the full-fledged steady-state equations for various levels of the detailed models
of power systems. We give sufficient conditions under which this bound can be
reached. We observe that the sufficient conditions are naturally satisfied for fully con-
nected power system with generic coefficients. We also develop a cluster method to

predict the number of solutions for not-fully-connected special power systems.

Basic homotopy method is computationally expensive for finding of solutions to
the (equilibrium) steady state equations of power systems. The special homotopy
method is applied to reduce the computational complexity and to guarantee finding all

119

120

the solutions to the load flow equations of power systems with probability one. We
develop the imbedding-based and the homotopy-based methods that pertain to pro-
cedures to simplify the computations in finding the solutions of power systems. More-

over, these procedures are directly implementable on digital sequential and parallel
processors.

7.2. Suggestions

An application of these theorems and the methods developed in this thesis yields
explicit quantitative information regarding the number of complex solutions for given
(polynomial) equations of power systems. However, for very deficient systems this
upper bound may still be too large. From the numerical results of the power system
examples, the number of the real system solutions of the load flow equations of power
systems depends on the values of the system parameters for a given structure, while
the number of the complex solutions is changed only with the system structure. An
important open question is: what is the upper bound on the number of the real system
solutions of the load flow equation for a given N-node power grid?

Solving all solutions of the load flow equations of power systems has been almost
impossible until the advent of the globally convergent homotopy method. However, the
current homotopy continuation methods are computationally expensive for finding the
solutions to very deficient systems. The amount of computational effort grows
exponentially to the point where efficient computation is rendered impossible even on
very fast and/or parallel computers. Consequently, the current homotopy continuation
methods become incapable for a larger-sized deficient system. An important question
facing mathematicians and scientific researchers is to develop a method with a theoreti-
cal proof to reduce the (computer) computational complexity in the finding of all (or
some of) the solutions of the so—called deficient systems. Even for a large-sized power

system, the methods can still compute all (or some of) the solutions.

APPENDIX

Appendix
The Proof of Theorem 3 .43

For the system (3.4.3.1), its associated homogeneous system is given by

0 81(Efi1y1j5j. 4’ YIN+150) + (A-1.1a)

E;(2jN=lylej + Y1.N+150) ' 21"1156:

0 = EN1(2.IN=1y1;:-J'E; + .YIG‘JVHEO) +

51;, (Eff-.QNrjEj + mph/+150) ' ”ml-5:6:

o '= ElE; - vag, (A-1.1b)

Ema}. - VglEg,

C
II

0 = EN,+1(£ji-IYI;,+1JE; 4' yN,+1.N+lEO) " SN,+156 . (A'l-lc)

o = EN (2,1665; + atheist) - Sheri.

121

122 Appendix

0 = 51;,“ (ZjN:1)’N,+1.jEj +)’N,+l.~+1Eo)' SN.+1E(% . (A'l-ld)

0 = Engine-E,- + moist) - 5,353.

The Bézout theorem says that the solutions of the associated homogeneous sys-
tem (A-1.1) include both the solutions for E 0 at O (i.e. the finite solutions) of the origi-
nal system (3.4.3.1) and the solutions for E o = O (i.e. the solutions at infinity). There-
fore, the number of the original system solutions equals the number of the solutions of

the associated homogeneous system minus the number of the solutions at infinity.

However, the solutions at infinity of (A-l. l) are the solutions with 50 = 0,

namely, the solutions of

0 51(2j’ityijED + EI(Z}‘.’..1yt,-E,-). (A-l.2a)

0 = EN,(2jA;1yN,JEj‘) 4’ 51;, (211'Y-=1yN,.jEj)’

o = E15}, - (A-1.2b)
. t

0 = EN‘EN‘ ,

o = E~,..<2,’-Y...y&,...,-E,-‘>. (A-ch)

o = Entzj’ilynjb‘h.

123 Appendix

0 = EN,+1(21N=IYN,+1JEJ')’ (A'l-Zd)

O = EN<£f=lYNj E j )-

From assumption (a) of the theorem, we can see that the solution set of (A-l.2)
consist of two linear subspaces, say, Z 1 and Z 2 which are the disjoint union of

z.=t<o.z,.-~.E~.z1.-~.s.;>lo=zt=...=s~l=0“.
and

22= {(0.El,---,EN,EI ,~-,E,;)I0=EI =...=E,;) =CPN-1.
The codimension of Z 1 and Z; in the complex projective space CPZN satisfy the
equality codim(Z1, cpl” ) = codim (22, CE") = N+1.

In the following, we will Show through tedious calculations that if all the condi-
tions of assumption (b) of the theorem hold, then the subspaces Z, and Z 2 are non-
singular. Consequently, by the fundamental theorem in (3.2.1), we conclude that the

number of isolated complex solutions of the polynomial system (3.4.3.1) (or (2.1.9)) is
given by

- 2N 2N 2N
"22N'[zll'[zd=2m'2£°l[i]' '25M[i]=[N]'
We now begin the lengthy calculations. We first reorder system (A-1.1) in the

following way for convenience.

0 51(2y”=1yI,-Ef+ yiJvnEo) + (A-1.3a)

51(ZjZ1YIjEj'i'YIN-l-1EO) ‘ 219153 =fl’

o = Beater/5.45; + new +

124 Appendix

5117, (if-ItlmeEj 4' YN,.N+150) ’ 2PM]; 2 = pr

0 = Efi‘+l(2jN=1yN‘+l.jEj + YN,+1JV+IEO) " SNfilEg = ill/#11 (A’1'3b)

0 = Bit/(Ejz-JYNjEj +YNJV+1EO)-SNEOZ =fN-

o = ElE; - vleg =f'N,,, (A-1.3c)

0 = £111,513, ‘VhiEti =f~+~,t

0 = EN,+1(Z,'N=1)’1;,+145; '1’ Minn/+150) ' SN,+156 = f N+N,+11 (A'1-3d)

0 = EN(Z,'N=1)'1;,'E; ‘1' YNNHEO) - 311156 =f2N-

The solution sets Z 1 and Z 2 at infinity are nonsingular, if

 

 

a " ’ . . . , "
J(Zl) = rankc (fl f2”) 1, (x) = codim,(Zl,CPm) = N+l
8(EOQS°H9ENsEl 9°°°vEN)
for every pointx e Z], and
a " ’ . . . , "
J(Zz) = rankc (fl f7”) (y) = codimy(Zz,CP7N) = N+l

8(5009...9EN9EI 9...9E&)
for every pointy e 22.

Calculating the Jacobian matrices I (Z 1) and J(Zz) of (A-l.3) and removing all

the zero columns, we respectively obtain

i(z1) =

and

WthC we have SCt 0;: = 211055;, bi: = EfilyijEj’ 65k:

(lsthsN).

€1.N+l “1+9"

eN‘.N+l €N,.l

eNJV+l 9111.1

0 E;
o

o 8
0 1')

 

r

flJV+l b1+f11

fN,.N+l f N,,l

o o
o o
0 El
6 t’)

 

f NJV+1 f A”

k

€N,+l.~+l eN‘+l,l ..

f N,+1.N+l vam ..

125

€1,111, C l.N‘+l

aNtflN: “N: e”: N1“

€N,+1.N, €N,+lJv,+1 ,,
chm, €N.N,+1
0 0
EN. a 0
0 Nit-1
i) 0
fl”. f 1N,+l
bN‘+}N‘ ,N. f”: ”a“
0 bN.+l
0 0
0 O
5,, o
fN,+1.N, f N,+lJv,+1 ..
fN..N, f NN‘H

81‘”

8N. .N
eN‘+l,N

eNN

0N

 

J

flJV

fN.,N

° C

2Nx(N+l)

 

J 2NX(N+1)

Appendix

(A- 1.4)

, (A-1.5)

= qu;9 311d f it: = ”25'

It is known, by the definition in 3.2, that Z, and Z; are nonsingular if rankcj(21)

and rm¢j(22) are both equal to N+l. It should be observed that i(Z,) and f(Z,_) have

the same form simply by reordering the rows of f(Zp). Therefore, f(Z,) and i(z2)

have the same rank. Hence, we only need to check the rank of f (Z 1). In i (Z 0, now

126 Appendix

we assume that the first E; = 0 (OSkISi, iSN-l) and the remaining E}: :0 (1:24,).
From assumption (a) of the theorem, at most (N -i-l) of the a,- terms are equal to zero.
There are now two cases to be discussed.

(A) i s N‘.

For convenience, we partition f(Zl) into two parts. The first N rows of i(Z1)' are
called Part 1, and the remaining N rows are called Part 2.

(a) i < N —t‘ .

We consider the worst case in which (1) the first i rows of Part 2 become zero
because of E; = 0 (09,50, and (2) the first i rows of Part 1 and the last [N-i-l-i)
rows of Part 2 become zero because 5; = 0 (0919') and at most (N-i-l) of the a,-
terms equal to zero. Note that the condition N -N‘ -(N -i-l-i) 2 0 is required, other-
wise, we take N—N1-(N-i-1-i) = 0. Remove all zero rows from NZ!) and remove all
rows and columns that share the element 5*; (kpatkl) from 1(21). Then remove all rows
and columns that have a common element a, (patj) in Part 2 from f (21). The final
reduced matrix is denoted by 31(21)

' 1
€i+uv+1 8.41.1 8.41.: 61,421+: 8.4m

€N,.N+1 319.1 Z: en": €N,,2t'+2 eN‘N
€N,+uv+1 ¢N,+1.1 -- €N,+1.t ¢N,+r,2z+2 ~- €N,+uv

],(z,) = (A-l.6)

 

 

e e .. e ' C ’ .. 3

Since the condition N-N‘-(N-i-1-i)2 0 implies i 2 (N,—1)/2, the matrix ],(z,) is
nonsingular from assumption (b) of the theorem. Therefore, we have

rankcflzl) = N+1 [ := (N-iHN,-i)+(N—N,-(N-i-l-i))].

127 Appendix

(b) i 2 N -i .

In this case, we still assume that the first i rows of Part 2 become zero since 5.2
= O (0919'). We also assume ai = O (OSjSN-i-l) (i.e. we are considering the worst
case). Similarly, we first remove all zero rows from i (Z 1). We then remove all rows

and columns that share the element Es (kpatkl) in Part 2 from j(Zl). Finally, we

remove all rows and columns that have a common element a, (paej) from f(Z,). The

remaining matrix is given by

' l

€i+r,N+1 6.41.1 cum-H

€N,.N+r 8N,.r I eN‘N-i-l

A-1.7
€N,+r,~+r eN‘HJ -- €N,+1.N-i—r ( )

Jztzo =

 

 

L e1~I.I~I+r €N.t " emu-r AN-iNN-i)
The condition i 2 N-i implies N-i s N I2. Therefore, the matrix 12(2 1) is non-

singular if condition (b) of the theorem holds. We conclude that rankcflzl) = N +1 [
:= (N-i)+(N-(N-i-l))].

(B) i > N‘.

In the following, we will follow the same procedures as above to show the non-
singularity of Z, and 22.

(a) i < N-i.

In the worst case, we assume that E; = O (OS/5131'). We also assume that a,- = 0

(OSjSN—i-l). Remove all zero rows from i (Z t). and remove all rows and columns
that have a common element a, (patj) from f (Z 1). The remaining matrix is denoted by

]3(Z r)

128 Appendix
€i+uv+1 ei+l,l 35+uv-i—1
€r+2JV+1 €i+2.1 ei+2JV-i-l
]3(z,) = (A-l.8)
Len/NH €N.1 eNN-i-ngvqygt/q)

 

 

From i > N‘, we have N-t' < N—N8 <N-N8/2. Therefore, the matrix .73(Z,) is

nonsingular if assumption (b) of the theorem holds. Consequently, we prove that

rankcj(Zl) equals N+1 [ := (N-i+N-(N-i-l))].

(b)i ZN-t'.

We still consider the worst case in which 5;, = O (OS/:19“), and a,- = 0

(OSjSN—i-l). Similarly, we remove all zero rows from 1(21). and remove all rows
and columns that have a common element a, (patj) from i (Z r)- The remaining matrix

is denoted by 74(21).

r w
€i+1.~+r 6.41.1 chum-1

ei+2JV+l et+2.r €r+2.N-t'-r

74(2 r) '-" (A'1-9)

 

 

\ eN,N+l €N.1 eNJV-i-l J (N—c'NN-i)

Obviously, the condition i 2 N -i implies N -i s N l2. We conclude that the
matrix 74(2 1) is nonsingular if assumption (b) of the theorem is satisfied. Therefore,

the rank of i(z,) equals N+l [ := (N-t'+N—(N-i—l))].

Since the matrices NZ!) and i (Z 2) are of the same form, we can follow the same
procedure to prove that rankcf(Zz) equals N+l. Consequently, we prove that Z 1 and 22

are nonsingular if assumption (b) of the theorem holds. This completes the proof.

BIBLIOGRAPHY

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

BIBLIOGRAPHY

S.X. Guo, and F.M.A. Salam, "Determining the Solutions of the Load Flow of
Power Systems: Theoretical Results and Computer Implementation", Proc. of
IEEE 29th Conference on Decision and Control (CDC), Honolulu, Hawaii, De-
cember 5-7, 1990.

S.X. Guo, and F.M.A. Salam, "The Imbedding-based Method For Finding the
Steady States of Models of Power Systems", Memorandum, MSU/EE/s, August
1990.

F.M.A. Salam, L. Ni, S. Guo, and X. Sun, "Parallel Processing for the Load
Flow of Power Systems: the Approach and Applications", Proc. of IEEE 28th
Conference on Decision and Control (CDC), pp. 2173-2178, Tampa, Florida,
December 15-17, 1989.

EMA. Salam, L. Ni, X. Sun, and S. Guo, "Parallel Processing for the Steady
State Solutions of Large-Scale Nonlinear Models of Power Systems", Proc. of
IEEE International Symposium on Circuits and Systems (ISCAS), pp. 1851-
1854, Portland, Oregon, CA, May 9-11, 1989.

L.M. Ni, F.M.A. Salam, T.H. Tzen, X. Sun, and S. Guo, "PowerCube- A
Software Package for Solving Load Flow Problems", The Proceedings of the
32nd Midwest Symposium on Circuits and Systems, PP. 14-16, Urbana, Illinois,
August 1989.

F.M.A. Salam, S. Bai, S. Guo, "Chaotic Dynamics Even in the Highly Damped
Swing Equations of Power Systems", Proc. of IEEE 27th Conference on Deci-
sion and Control (CDC), pp. 681-683, Austin, Texas, December 1988.

F.M.A. Salam, "Cmrent Issues in the Stability and Control of Interconnected
Power Systems", Proc. of IEEE 24th Conference on Decision and Control
(CDC), pp. 830-831, Ft. Lauderdale, Florida, December 1985.

F.M.A. Salam, "Asymptotic Stability and Estimating the Region of Attraction
oggtlge Swing Equations", Journal of Systems and Control Letters, Vol. 7, No. 3,
1 .

F.M. Salam, LE. Marsden, and P. Varaiya, "Arnold Difi‘usion in the Swing

Equations of a Power System", IEEE Trans. Circuits and Systems, Vol. CAS-
31, pp. 673-688, August 1984.

129

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[13]

[19]

[20]

[21]

[22]

[23]

130 Bibliography

EMA. Salam, "Power Systems Transient Stability: the Critical Clearing
Time", EBB 23rd Conference on Decision and Control (CDC). Pp. 179-184,
Las Vegas, NV., 1984.

J. S. Thorp and SA. Naqavi, "Load Flow Fractals", Proc. of EBB 28th Confer-
ence on Decision and Control (CDC), pp. 1-822 1827, Tampa, Florida, De-
cember 15-17,1989.

F.S. Prabhakara, and AH. Bl-Abiad, "A Simplified Determination of Transient
Stability Regions for Lyapunov Methods", EBB Trans. Power App. Syst. Vol.
PAS-94. PP. 672-689, 1975.

Y. Tamura, K. Iba, and S. Iwamoto, "A Method for Finding Multiple Load-
flow Solutions for General Power Systems", EBB PBS Winter Meeting A80
043-0, Feb. 1980.

Y. Tamura, Y. Nakanishi, and S. Iwamoto, "On the Multiple Solution Structure,
Singular Point and Existence Condition of the Multiple Load-flow Solutions",
EBB PBS Winter Meeting A80 044-8, Feb. 1980.

J. Baillieul and CI. Bymes, "Geometric Critical Point Analysis of Lossless
Power System Models", EBB Trans. Circuits & Syst., Vol. GAS-29, No. 11,
Nov. 1982.

J. Baillieul and CI. Bymes, "The Load Flow Equations for a 3-Node Electrical
Power System", Syst. & Control Letters 2, pp. 321-329, 1983.

J. Baillieul and CI. Bymes, "The Singularity Theory of the Load Flow Equa-
gions3§3r a 83-Node Electrical Power System", Syst. & Control Letters 2, pp.
30- , 19 3.

T.Y. Li, T. Sauer, J. Yorke, "Numerical Solution of a Class of Deficient Poly-
nomial Systems", SIAM J.Num. Anal. 24, pp. 435-451, 1987.

T.Y. Li, T. Sauer, J. Yorke, "The Random Product Homotopy and Deficient Po-
lynomial Systems", Numer. Math., 51, pp. 481-500, 1987.

T.Y. Li, T. Sauer, "Regularity Results for Solving Systems of Polynomials by
Homotopy Method", Numer. Math., 50, pp. 283-289, 1987.

T.Y. Li, T. Sauer, J. Yorke, "The Cheater’s Homotopy: an Efficient Procedure
for Solving Systems of Polynomial Equations", SIAM J. Num. Anal. Vol. 26,
No. 5, pp. 1241-1251, 1989.

T.Y. Li, T. Sauer, J. Yorke, "Numerically Determining Solutions of Systems of
Pglynomzial SI’istgrations", Bulletin of The American Mathematical Society, Vol.
1 , No. , l .

T.Y. Li, T. Sauer, "A Simple Homotopy for Solving Deficient Polynomial Sys-
tems", Japan J. Appl. Math., 6, pp. 409-419, 1989.

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]
[33]

[34]

[35]

[36]

[37]

[38]

[39]

131 Bibliography

T.Y. Li, "On Chow, Mallet-Paret, and Yorke homotopy for Solving Systems of
Polynomials", Bulletin of the Institute of Mathematics, Academica Sinica 11,
pp. 433-437, 1983.

B. Stott, "Review of Load Flow Calculation Methods", Proc. EBB, Vol. 62,
pp. 916-929, 1974.

A. Arapostathis, S. Sastry and P. Varaiya, "Bifurcation Analysis of the Load
Flow Equations", 19th EBB Conference on Decision and Control (CDC), pp.
641-644. 1980.

A. Arapostathis, S.S. Sastry and RP. Varaiya, "Global Analysis of Swing
Dynamics", EBB Trans. Circuits and Systems 29, pp. 673-679, ,1982.

A. Arapostathis, S.S. Sastry and P. Varaiya, "Analysis of power-flow equation",
Electrical Power and Energy Systems, Vol. 3, pp. 115-126, July 1981.

A. Arapostathis and P. Varaiya, "Behavior of Three-node Power Netwo ",
Electrical Power and Energy Systems, Vol. 5, pp. 22-30, January 1983.

R.K. Mehra, "Bifurcation and Catastrophes in Power System Stability Prob-
lems", in: Optimization Days, Montreal, Canada, May 1977.

CJ. Tavora and OJ .M. Smith, "Equilibrium Analysis of Power Systems", EBB
Trans. Power App. Syst. Vol. PAS-91, pp. 1131-1137, 1972.

W. Fulton, "Intersection Theory", Springer-Verlag, New York, 1984.

E. Abed and PP. Varaiya, "Nonlinear Oscillations in Power Systems", Electric
Power and Energy Systems, Vol. 6(1), pp. 37-43, 1984.

FJ. Drexler, "A Homotopy Method for the Calculation of All Zero-dimensional
Polynomial Ideals, in Continuation Meth ", I-I. Wacker, ed., Academic Press,
New York, pp. 69-93, 1978.

A.H. Wright, "Finding all Solutions to a System of Polynomial Equations",
Math. Comput. 44, pp. 125-133, 1985.

A.P. Morgan, "A Homotopy for Solving Polynomial Systems", Applied Math.
Comp. 18, pp. 87-92, 1986.

A. Morgan and A. Sommese, "Computing all Solutions to Polynomial Systems
Using Homotopy Continuation", Applied Mathematics and Computation, Vol.
24, pp. 115-138, 1987.

Walter Zulehner, "A Simple Homotopy Method for Determining all Isolated
Solutions to Polynomial Systems", Math. Comp. 50, pp. 167-177, 1988.

B. Allgower and K. Georg, "Simplicial and Continuation Methods for Approxi-
mating Fixed Points", SIAM Rev., 22, pp. 28-85, 1980.

[40]

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[501'

[51]

[52]

[53]

132 Bibliography

A. R. Bergen, "Power Systems Analysis", Bnglewood Cliffs, N. J. Proentice-
Hall, 1986.

AR. Bergen and DJ. Hill, "A Structure Preserving Model for Power Systems
Stability Analysis", EBB Trans. Power App. Syst., Vol. PAS-100, pp. 25-35,
January 1981.

J. Carpentier, "Optimal Power Flows", Electrical Power and Energy Systems,
Vol. 1, No. 1, April 1979.

Shui-Nee Chow, John Mallet-Paret, and James A. Yorke, "A Homotopy
Method for Locating all Zeros of a System of Polynomials", Functional
Differential Equations and Approximation of Fixed Points, Springer Lecture
Notes in Mathematics, No. 730, pp. 77-88, 1979.

Shui-Nee Chow, John Mallet-Paret, and James A. Yorke, "Finding Zeroes of
Maps: Homotopy Methods That are Constructive with Probability One", Math.
Comp., (32). Pp. 887-899, 1978.

SN. Chow, L.M. Ni and Y.Q. Shen, "A Parallel Homotopy Method for Solving
a System of Polynomial Equations", Proc. of the 3rd SIAM Int’l Conf. on
Parallel Processing for Scientific Computing, 1987.

P. Dersin and AH. Levis, "Feasibility Sets of Steady-state Loads in Electric
Power Netwo ", EBB Trans. Power App. Syst., Vol. PAS-101, pp. 60-70,
January 1982. -

T.E. Dy Liacco, "System Security: the Computer’s Role", EEE Spectrum,
Vol. 15, pp. 43-50, June 1978.

J. Fong, and C. Pottle, "Parallel Processing of Power System Analysis Problems
via Simple Parallel Microcomputer Structure", EBB Trans. Power App. Syst.,
Vol. PAS-97, pp. 1834-1841, September 1978.

A.A. Fouad, "Stability Theory-criteria for Transient Stability", Proc. Conference
on System Engineering for Power: Status and Prospects, Henniker, New
Hampshire, 1975.

PD. Galiana, "Analytic Properties of the Load Flow Problem", Proc. of Int.

. Symp. Circuits and Systems, Special Session on Power Systems, EBB Catalog

No. 77CH1188-2 CA8. Pp. 802-816, 1977.

PD. Galiana, "Analytical Investigation of the Power Flow Equations", Proc.
1983 American Control Conference, San Francisco, pp. 411-415, June 1983.

L.P. Hajdu, R. Podmore, "Security Enhancement for Power Systems", Proc. En-
gineering Foundations Conf. on Systems Engineering for Power, New England
College, Henniker, New Hampshire, August 17-22, 1975.

H.H. Happ, "Optimal Power Dispatch - a Comprehensive Survey", EEE
Trans. Power App. Syst., Vol. PAS-96, No. 3, pp. 841-854, May/June 1977.

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

[64]

[65]

[66]

[67]

[63]

133 Bibliography

IMSL Library, "Fortran Subroutines for Mathematics and Statistics", User’s
Manual, Edition 9.2, IMSL, November 1984.

B.K. Johnson, "Extraneous and False Load Flow Solutions", EEE Trans.
Power App. Syst., Vol. PAS-96, pp. 524-534, March/April 1977.

R.J. Kaye and FF. Wu, "Dynamic Security Regions of Power Systems", EBB
Trans. Circuits and Systems, Vol. CAS-29, pp. 612-623, September 1982.

A. Klos and A. Kemer, "Non-uniqueness of Load Flow Solution", Proc. PSCC
V. 3.1/8. Cambridge, July 1975. '

A.J. Korsak, "On the Question of Uniqueness of Stable Load Flow Solutions",
EBB Trans. Power App. Syst., Vol. PAS-91, pp. 1093-1100, May/June 1972.

N. Kopell and RB. Washbum, Jr., "Chaotic Motions in Two-degree-of-freedom
Swing Equations", EEE Trans. Circuits and Systems, Vol. GAS-29, pp. 738-
746, November 1982.

H.G. Kwatrry, A.K. Pasrija, and L.Y. Bahar, "Static Bifurcation in Electric
Power Networks: Loss of Steady-state Stability and Voltage Collapse", EBB
Transactions on Circuits and Systems, vol. CAS-33, pp. 981-991, 1986.

(3.8. Lauer, N.R. Sandell, D.P. Bertsekav, and T.A. Posbergh, "Solution of
Large-scale Optimal Unit Commitment Problems," EBB Trans. Power App.
Syst., Vol. PAS-101, pp. 79-86, January 1982.

GE Marks, "A Method of Combining High-speed Contingency Load Flow
Analysis with Stochastic Probability Methods to Calculate a Quantitative Meas-
ure of Overall Power System Reliability", EBB PBS Winter Meeting, New
York, Paper A78 pp. 228-229, February 1978.

L.M. Ni and CE Wu, "Design Trade-offs for Process Scheduling in Tightly
Coupled Multiprocessor Systems", Proc. of the 1985 Int’l Conf. on Parallel Pro-
cessing, pp. 63-70, August 1985.

 

R. Olson, "Parallel Processing in a Message-based Operating System", EBB
Software. pp. 39-49, July 1985.

MA. Pai, "Power System Stability Analysis by Direct Method of Lyapunov",
Amsterdam and New York: North Holland, 1981.

R. Pritchard and C. Pottle, "High-speed Power Flows Using Attached Scientific
(array) Processors", EEE Trans. Power App. Syst., Vol. PAS-101, pp. 249-253,
January 1982.

B. Stott, O. Alsac, "Fast Decoupled Load Flow", EEE Trans. Power App.
Syst., Vol. PAS-93. PP. 859-869, May 1974.

B. Stott, O. Alsac, and J .L. Marinho, "The Optimal Power Flow", Presented at
the EPRIISIAM International Conference on Electric Power Problems and the
Mathematical Challenge, pp. 18-20, Seattle, Washington, March 1980.

[69]

[70]

[71]

[72]

[73]

[74]

[75]

[76]

[77]

[78]

134 Bibliography

W.F. Tinney, C.B. Hart, "Power Flow Solution by Newton’s Method", EBB
Trans. Power Appl Syst., Vol. PAS-86, pp. 1449-1460, November 1967.

N. Tsolas, A. Arapostathis and P. Varaiya, "A Structure Preserving Energy
Function for Power System Transient Stability Analysis", EBB Trans. Circuits
and Systems, pp. 1041-1050, October 1985.

RF. Wu, G. Gross, J.F. Luini, and RM. Look, "A Two-stage Approach to
Solving Large-scale Optimal Power Flows", Proc. 11th PICA Conference, pp.
126-136, Clevelnad, May 1979. '

F.F. Wu, "Theoretical Study of the Convergence of the Fast Decoupled Load
Flow", EBB Trans. Power App. Syst., Vol. PAS-98, No. 1, pp. 268-275, Janu-
ary 1977.

EB. Wu and S. Kumagai, "Steady-state Security Regions of Power Systems",
EEE Trans. Circuits and Systems, Vol. GAS-29, pp. 703-711, November
1982.

RF. Wu and RD. Masiello, eds., "Special Issue on Computers in Power Sys-
tem Operations", Proc. of the EBB, December 1987.

EBB Committee, "Computer Representation of Excitation System", EEE
Trans. Power Appar. Syst. pp. 1460—1464, June 1968.

L.T. Watson and MR. Scott, "Solving Spline-collocation Approximations to
Nonlinear Two-point Boundary-value Problems by a Homotopy Method", Ap-
plied Mathematics and Computation, Vol. 24, pp. 333-357, 1987.

L.T. Watson, "Numerical Linear Algebra Aspects of Globally Convergent
Homotopy Meth ", SIAM Review, pp. 529-545, December 1986.

GT. Heydt, "Computer Analysis Methods for Power Systems", Macmillan Pub-
lishing Company, New York, 1986.

"Illllllllllfil‘llllllllll“