135(!
Transactions
on
Power
Systems,
Vol.
7,
No.
3,
August
1992
EFFICIENT METHODS FOR FINDING TRANSFER FUNCTION
ZEROS
OF
POWER
SYSTEMS
Nelson
Martins
Herminio
J.C.P. Pinto
Member IEEE
CEPEL
-
Centro de Pesquisas de Energia Eletrica
Caixa Postal 2754
20.001,
Rio
de Janeiro,
RJ,
Brazil
Leonard0
T.G. Lima
Member IEEE
AV. Pres. Vargas 54212202
20.000,
Rio
de Janeiro,
RJ,
Brazil
MARTE Engenharia
Abstract -This paper is probably the first to describe algorithms
suited to the efficient calculation of both proper and non-proper
transfer function zeros of linearized dynamic models for large
interconnected power systems. The paper also describes an
improvement to the well known AESOPS algorithm, formulating
it as an exact transfer function zero finding problem which is
efficiently solved by
a
Newton-Raphson iterative scheme. Large
power system results are presented in the paper.
w:
Power System Stability, Low Damped Oscillations,
Additional Feedback, Excitation Control, Transfer Function
Zeros, Large Scale Systems, Sparse Eigenanalysis.
I.
INTRODUCTION
The location of the zeros of the open-loop transfer function of
a feedback system is closely related to the ease or difficulty with
which the system is controlled. The movement of zeros following
system changes is
a
rather complex subject and little work has
been done in association with the power system problem [1,2].
The use of the
augmented system equations
(see Appendix)
for the smallsignal stability problem has already allowed the
efficient calculation of eigenvalues, eigenvectors, frequency
response plots, transfer function residues, participation factors and
step response plots [2,3,4,5,6,7,8,9] for large scale systems.
This paper comes in response to the need for efficient
algorithms for the calculation of transfer function zeros of large
power system dynamic models [4]. Newton-Raphson, Inverse
Iteration and Simultaneous Iteration algorithms [3,5,9], applied to
the
augmented system equation,
are described. With such
algorithms an optional facility may be added to comprehensive
packages for smallsignal stability analysis enabling engineers to
carry out controller design with extra valuable information.
The EPRI software for the analysis of smallsignal stability
of large scale power systems uses two alternative techniques to
compute eigenvalues which complement each other
[4]:
2. the Modified Arnoldi method, which has
a
sound math-
ematical basis and can find typically
up
to five
eigenvalues simultaneously.
A
recent paper [lo]
has
explained the AESOPS algorithm in
terms of traditional eigenvalue analysis. The results presented in
[lo]
for the 10-machine New En land system did not, however,
clearly attested the superiority
of
the algorithm proposed by its
authors over the original AESOPS in the eigensolution of large
practical power systems.
The power system eigenvalues were shown in (41 to be equal
to the zeros
of
a
special transfer function. This fact were not used
to advanta
e
in
[4]
due to the lack of an exact analytical
expression
81
such transfer function and of an adequate transfer
function zero finding method for large scale systems. These two
obstacles were obviated in the work reported in this paper, leading
to the improved
AESOPS
algorithm of sections V and VI.
All
the algorithms of this paper have been implemented
exploiting the
augmented system equations
sparse structure. The
notations adopted in the paper are defined
as
used.
.
II.
TRANSFER FUNCTION
ZEROS
Consider the dynamic system equations:
where A is
a
state
matrix of order
n,
is the state vector,
U
is
a
single input and
y
is
a
single output whichhave been specified.
The objective here is to find the zeros of the open-loop
transfer function
y(s)/u(s)
=
~t
(sI-A)-lb. From Root LOCUS
theory it
is
known that the closed-loop transfer function poles tend
to the open-loop transfer function zeros as the feedback loop gain
tends to infinity
[ll].
This concept was used to derive the basic
algorithm of this paper which is similar to that described in [12].
The closed-loop system will be defined here
as
having a
1.
the AESOPS algorithm which is
a
successful heuristical-
ly based algorithm and computes one eigenvalue at
a
control
signal
uproportion~
to
the
output
y:
time;
Papers
presented
at
the Seventeenth
PICA
Conference
at
the
Hyatt Regency Baltimore
Hotel,
Baltimore, Maryland, May
7
-
10
1991
Sponsored
by the
IEEE
Power
Engineerhg
Society
'
u(s)
=
Kds)
The poles of the closed-loop system are then the eigenvalues
of the state equation:
(3)
The eigenvalues of matrix
A,1
will coincide with the open-
loop transfer function zeros when the feedback ain
K
approaches
infinit
.
In
this case, matrix A,1 differs from by the introduc-
tion orvery large elements in the locations defined by the product
--
b
ct.
Matrix A,1
is
of the same order of the whole system, is real
and unsymmetric and its eigenvalues can be obtained by
a
stan-
dard QR routine
[13].
As
a
transfer function normally has less
zeros than poles, the
QR
eigensolution will contain extraneous
0885-8950/92$03.00
0
1992
EEE
1351
zeros which assume larger values
as
the feedback gain
K
is
increased. These extraneous zeros should theoretically
go
to
infinity with the feedback gain
K,
but this does not happen due to
rounding errors.
The closed-loop system poles, i.e., the eigenvalues of A,1 can
also
be
found by solving the generalized eigenvalue problem
Ag
4
=
A
B
4:
where
I
is the identity matrix,
Qt
is
a
row vector with all elements
equal to zero and
-
comprises both
x
and scalar input
U.
Note
that
as
the value
2K
tends to infinity, the matrix element
1/K
tends to zero. A
QZ
routine
[13]
for solving the generalized
eigenproblem of equation
(4)
directly deals with the case where the
matrix element
1/K
is identical to zero, and therefore the
extraneous zeros assume such large magnitudes that can easily be
identified and discarded. The solution of the generalized eigen-
value problem
of
equation
(4)
should therefore be preferred to the
method of
[12]
for finding
all
the transfer function zeros of
a
moderate size system.
III.
CALCULATION OF ZEROS
FOR
LARGE SCALE SYSTEMS
The use of
a
QZ
routine to solve for all the zeros of the
specified transfer function is
a
prohibitively expensive task in large
scale
s
stems. The only alternative in large system problems is to
solve
&
one zero at
a
time or for several zeros at a time located
around
a
fixed point which can be placed at will in various parts
of
the complex plane. Efficient algorithms can be developed to
exploit the sparse structure of the
augmented system equations
which are described in the Appendix.
The generalized eigenvalue problem of
(4
can be solved, one
zero
at a
time, by the inverse iteration algorit
h
m
[3],
whose basic
scheme is shown below for the case where the matrix element
1/K
is equal to zero:
a.
Solve for
Ek
+I:
b. Compute the vector Zk+lfor the next iteration:
Convergence occurs when the change in
at any iteration is
less
than some specified tolerance. In this algorithm the subscript
kis the iteration number,
q
the specified a proximate value of the
desired transfer function zero
zi
and
maJt&
is the element of
largest magnitude in this vector. The vector
_Zk
has arbitrary
initial value and corresponds to the zero direction vector at
convergence. After convergence, the factor
ll(2i-q)
will be
dominant in the element
maT(c(~~+~)
and the correct zero
zi
is given
by:
zi
=
q
+
1/ma.(Ek+J
Note that the transfer function zero and zero direction
vector of the matrix in equation
(5)
form
a
pair which is analogous
to the eigenvalueeigenvector pair of the state matrix A.
Equation
(5)
is now expressed in terms
of
the
augmented
system equations:
IV.
FINDING INVARIANT ZEROS
IN
THE
MULTI-INPUT-MULTI-OUTPUT CASE
The zero finding algorithms described in sections
I1
and
111
can readily be extended to the
multi-input-multi-output
case.
When
m
inputs and
m
outputs are simultaneously considered,
vectors
_b
and
tt
of equation
(4)
become matrices
B
and
C
of
appropriate dimensions. The invariant zeros
[15]
of a large scale
system matrix can be calculated by the inverse iteration and
simultaneous iteration algorithms. The transmission zeros of a
transfer function matrix are a subset of the system matrix
invariant zeros
[15].
A
brief result on
a
5-machine system
is
presented in section
VIII,
but further research is needed into this area.
V.
THE
AESOPS ALGORITHM FORMULATED
AS A ZERO FINDING PROBLEM
The AESOPS algorithm
(41
is
a
heuristically based one-at-a-
time eigenvalue method designed to compute the electromechan-
ical modes of oscillation for large power systems. The
AESOPS
algorithm is derived from the linearized equation of motion of
a
chosen generator, to which
a
complex frequency disturbance in the
mechanical torque is applied. At every iteration, a corrected value
for this complex frequency disturbance is applied until the system
becomes resonant. This iterative process is almost always
convergent and the converged complex frequency value corre
sponds to an electromechanical eigenvalue which is dominant at
the disturbed generator.
An interesting paper
[lo]
has suggested improvements to the
basic AESOPS algorithm, but lacked large scale system results to
substantiate its claims. In this section, an improved
AESOPS
algorithm is proposed which requires the calculation
of
the zeros of
a
specially tailored transfer function
[4].
Consider the block diagram of Figure
1
which describes the
torqueangle loop dynamics of the disturbed j-th generator in
a
large power system. The mechanical damping constant
D.
is here
assumed
to
be zero for brevity, but was fully considered in the
1352
computer
al
orithm implementation. The variables
A6(s
,
Ads),
subscript
j
to
relate them to the jth generator, but this subscript
was omitted for simplicity. The inertia constant of the jth
generator is denoted by
Hj.
ATel(s)
ant
ATm(s)
of this section should rigorously
3
1
have
a
The generalized eigenvalue problem described by equation
(10) cannot
be
adequately solved by the inverse iteration
algorithm since the matrix on the left part of the equation is
a
functional of the Laplace variable
s.
A
more convenient way to
solve this problem would be by using the Newton-Raphson
method,
as
described in the next section.
I
I
VI. A
NEWTON-RAPHSON SOLUTION SCHEME
FOR
THE
IMPROVED
AESOPS ALGORITHM
When
s
is
a
zero of the transfer function
ATm(s)/A6(s)
described in equation
(8)
it satisfies:
~~(sI-A,)
-lhs
+
4s)
=
o
(12)
Figure
1.
TorqueAngle Loop
of
Disturbed jth Generator
From the inspection
of
the Figure
1
one can write:
By choosing the mechanical torque and rotor angle as output
and input variables respectively, one gets:
Equation
(8)
can be expressed in the form:
ATm(s)
=
@x(s)
+
4s)
A6(s)
(9)
in which the output variable
ATm(s)
depends not only on the
vector
~(s)
but also
on
the system input
A6(s)
and its derivatives.
Let
A
be the
(mn)
state matrix of the global multimachine
power system. The
AESOPS
al
orithm requires the opening
of
the torqueangle loop of the distur%ed jth generator. The opening
of
this torqueangle loop implies making zero the
A6
and
Aw
states
of
the jth generator and letting the column of the
A6
state become
the input vector
5
to the system.
The zeros of
(9)
can therefore be found
by
solving the
generalized eigenvalue problem:
1
-
6
where
A'
is
a
matrix
of
order
(n-2)
due to the elimination of states
A6
and
Au
of the 5th generator. The vector
~(s)
used in this
section and the next is also of order
(n-2).
The term
4s)
is given
by:
d(s)
=
Cs
+
-+sa
2
H.
where
c6
is
a
real constant which depends on the system operating
point.
(11)
Transfer function zeros can be found, one at
a
time, through
use of
an
iterative algorithm such as Newton-Raphson. Solving
equation
(12)
is equivalent to solving:
(SI-A')
Z(S)
-
bs
=
0
-ct
x(s)
+
4s)
=
0
(13)
which is
a
non-linear system with
(n-1)
equations in
An-1)
unknowns. The unknowns are the Laplace variable
s
an the
vector
x
s)
which is of order
(n-2)
since the states
A6
and
Aw
of
the jth disturbed generator were removed.
The Newton-Raphson algorithm for solving
(13)
is given by:
a.
provide initial estimates
~0,
SO
b. calculate the vector of residues
f(Sk,&)
where the value of the input variable
6(s)
is set to unity.
c. Stop process if change in
f(sk,&)
is below the specified
tolerance.
increments
Axk
and
Ask:
d. Evaluate the Jacobian of (14) and solve for the new
=-
e.
Obtain
4.1
=
4
+
A&+i,
Sk+l
=
Sk
+
ASk+1
and return
to step "b".
For the solution of large scale problems, equations (14) and
(15) must be expressed in terms of the
augmented system equu-
tiom,
described in Appendix. Equation (15) is expressed below in
the desired form:
1353
-
~~lef(8)
--;c
where
J<
-
Jz’(J~)‘~J~’
=
A’.
The original AESOPS algorithm has the good characteristic
of converging to the dominant electromechanical modes of the
disturbed enerator in spite of bad initial values for
~0
and
so.
The
improved
XESOPS
algorithm described in this section also has the
same characteristics. This desired robustness was obtained by
using the augmented initial vector:
DZ(4
Static Compensator
SVC(s)
W.
EXTENDING THE AESOPS ALGORITHM CONCEPT
TO OTHER ACTIVE SYSTEM COMPONENTS
A
ain
one can note that the zeros of the transfer function
AV
ref
tk
s
/ABY(s)
of
(18)
are equal to the poles of the
closed
loop
system of Figure
2.
The remaining considerations are similar to
the material contained in section
V
of this paper. As the Newton-
Raphson method is very sensitive to the initial values given, there
is a need for
an
initialization vector in order to make this
algorithm converge to the desired dominant modes of the static
compensator.
MII.
RESULTS ON TRANSFER FUNCTION ZEROS
The open loop transfer function zeros of a plant are not
altered by the addition of a feedback controlrer. Consider the case
where the system has an unstable pair of poles and that feedback
stabilization is attempted through an input-output pair whose
transfer function exhibits an unstable pair of zeros in the
nei hborhood of the
poles
to be damped.
A
root locus branch
[ll]
wifexist between these neighboring pairs of poles and zeros, irre
spective
of
the feedback controller transfer function. Therefore, it
is not possible in practice to stabilize this system through this
feedback control loop.
The knowledge on the location of transfer function zeros
enables control engineers to carry out controller design more
effectively. The results presented in this section are intended to
show the potential of the algorithms developed and estimulate
power system control en ’neers to further investigate the practical
application of this extragcility.
1354
Cases Studied
No
Description
1
Base Case
2
Lower
Transfer
3
No AVR
in
Gl,GZ,Gs
4
PSS's
in
GlrG2,G3
5
PSS
in
GI
This system has
a
pair of unstable eigenvalues
A
=
+0.646 +j5.391
and any attempt to stabilize it throu h
excitation control on
G4
is bound to fail. Reference
[2]
shows
t%e
root locus of the critical ei envalues
as
the gain of
a
rotor speed-
derived stabilizer
at
the
84
generator is varied. The critical
electromechanical mode is seen to always remain unstable due to
the presence of
an
unstable pair of zeros in the
Ad(s)/A Vr4(~)
transfer function
z
=
+0.049
+j
5.902
The upperscript
4
exciter reference voltage. In the small signal stability area the
complex zeros and poles always occur in complex conjugate pairs.
A
complex conjugate pair (a
*j
b) is here typed as (a
+j
b) for
better readability
.
The 5-machine system has
28
ei envalues (poles and the
eigenroutine. The three extraneous zeros had magnitudes larger
than
105.
The critical zeros for different transfer functions are
presented in Table
1
and discussed in the following lines:
1,2.
The critical pair of zeros are identical for
Ad(s)/AVr4(s)
and
APt4(s)/A V,~(S),
where
Pt
is the
generator terminal power.
3.
The symbol
R4-8
denotes the apparent resistance of the
transmission line between buses
4
and
6.
This signal, for
this particular system, is worse then the two previous signals
since its critical zeros are more unstable.
4,5,6.
The critical pair of zeros for
Aw'(s)/A Vr'(S),
Aw2(s)/AVrZ(s), Aw~(s)/A Vr3(~)
transfer functions are
almost identical and very close to the unstable pair of
eigenvalues
(A
=
+0.646
+j
5.391).
This unstable pair of
eigenvalues is therefore not controllable from the excitation
systems of
GI, G2
and
G3.
7.
There is no troublesome pair of complex zeros in the
Ad(s)/A Vr7(s)
transfer function. Problem appears due to
a
real positive zero
(z
=
+7.012),
which informs in advance
of the detrimental action that
a
stabilizer located at the
synchronous motor
G7
would have
on
the system synchron-
izing torques.
8.
The possibility of stabilizing the system through the
function
Aw4
s
/A Vr'(S
is discarded due to the existance of
here would modulate the reference voltage of the
GI
exciter
and be derived from the rotor speed signal of
G4.
denotes
a
variable
o
\
generator
G4
and
Vr
the deviations in the
Ad(s)/A Vr4(s)
has
25
finite zeros whici were obtained
i
y the
QZ
a
higly unsta
v
e
pair
o
2
zeros. The stabilizer to be added
Crit. Poles Crit. Zeros
+.646+j5.391 +.049+j5.908
+.428+j5.610 -.023+j5.958
+.667+j5.315 -.242+j5.660
+.656+j5.380 -.562+j5.044
+.652+j5.386 -.427+j5.835
No
Transfer Function Critical Zeros
1
1
Considered
I
$0.049
+j
5.908
+0.049
+j
5.908
+0.249
+j
6.404
+0.655
+j
5.379
+0.650
+j
5.376
+0.654
+j
5.380
4.310
+j
5.748*
+0.899
+j
5.354
Table
1.
Critical Pair of Zeros for Different Transfer
Functions of the 5-Machine System
Note:
*
The function
Aw7(s)/A
Vr7(s)
has another critical zero of
value
z
=
+7.012
Table
2
shows the critical pair of poles for the 5-machine
system
together with the critical pair of zeros for the
Ad(s)/A Vr4(3)
transfer function. The various cases presented
are
descnbed below:
1.
Critical poles and zeros are presented for the base case
condition described in
[2].
2.
A
10
percent reduction on the power interchange between
G4
and
G,
machines causes the troublesome pair of zeros to
move slightly into the left-half plane, but the system
continues to present basically the same stabilization
problem.
3.
The critical pair of zeros become stable when automatic
excitation control is neglected on the generators
GI, G2
and
Gs.
A
sin le
PSS
at
G4
can now stabilize the system
through moiulation
of
the impedance loads at buses
1,
2,
3
and
5.
The maximum dampin achieved for the electre
mechanical eigenvalue is about
4k
since
it
will coincide with
z=
4.242
+j
5.660
for infinite gain
at
the
G4
stabilizer.
4.The power system stabilizers (PSS) in
GI, G2
and
Ga
practically do not alter the unstable eigenvalue pair but
have
a
strongly positive effect on the critical zeros:
z=
-0.562
+j
5.044.
The stabilizers at
GI, G2
and
Gt
are
therefore needed not for being able to damp the unstable
poles but for moving away the troublesome zeros.
5.
The presence of
a
stabilizer only in
G1
also has a highly
positive effect on
Aw4(s /A V,*(S),
since the critical zeros
become very well damped.
Result
No
5
of Table
2
informs in advance that
a
PSS
in
G4
could stabilize the system if another PSS
was
already present at
G1.
This result indicates that the transfer function matrix:
has well damped transmission zeros
(see
Section IV). This is
actually the
case,
since the least damped transmission zeros of the
matrix, calculated by
a
QZ routine, are
z
=
-1.839
+j
9.157
and
z=
-1.273
+j
6.635.
Results
on
the Brazilian Interconnected System
The power system analysed is
a
616
bus-50 generator model
of the South-Southeast Brazilian Interconnected System
(8
.
The
system stabilizers, are presented in Table
5
of the next section.
The reader should refer to
[8]
for additional information on this
system model.
Reference
[8)
described results showing that the inter-area
mode
(A1=
-0.0017
+j
3.511)
could be stabilized through a
properly tuned SVC located at the terminals of the Jacui
generating plant. The effectiveness of a SVC
at
this bus in
damping this interarea mode actually depends
on
whether its bus
voltage signal is local or remote.
There exists
a
zero
z1=
-0.03537
+j
3.535
in the transfer
function
A V(s)/A&l(s),
where
V1
denotes the voltage magnitude
of the Jacui generator bus. The proximity of zero
z1
caused pole
A1
to be invariably attracted to it as the SVC gain was raised.
least damped eigenvalues of this system, in the absence
o
t'
power