9 V. ' PARAMETER ESTIMATION APPLIED TO POWER SYSTEM MODELS SUBJECT TO RANDOM DISTURBANCES by HANS MULLER D i p l . El.-Ing. ETH Swiss Federal I n s t i t u t e of Technology, Zurich 1973 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE SCIENCE Department of E l e c t r i c a l Engineering We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June, 1979 (c) Hans Muller, 1979 In p r e s e n t i n g this thesis an a d v a n c e d d e g r e e a t the L i b r a r y I further for agree scholarly by h i s of shall this written in p a r t i a l the U n i v e r s i t y make it freely that permission fulfilment of of Columbia, British available for for extensive the requirements reference copying o f I agree and this thesis for It financial is understood gain shall that not copying or be a l l o w e d w i t h o u t Hans M u l l e r - ^ E l e c t r i c a l Engineer in, of The U n i v e r s i t y of British 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date 6 thesis 9/7/79 Columbia or publication permission. Department that study. p u r p o s e s may be g r a n t e d by t h e Head o f my D e p a r t m e n t representatives. for my ABSTRACT For control of an e l e c t r i c generator, i t i s desirable to approximate the remote part of the power system by a low order model, whose parameters can be estimated from measurements. The proposed model for the remote system consists of a large synchronous machine and a random varying load to account for the dynamic behaviour and for the small f l u c t u a t i o n s that are always present. The maximum l i k e l i h o o d algorithm i s a good method for parameter estimation of a noisy system. I t can even be used f o r estimation from measurements alone, without applied disturbance, which could be of great advantage i n a power system. The algorithm i s derived i n a form s u i t a b l e for e f f i c i e n t d i g i t a l processing of sampled measurements from a l i n e a r continuous system with physical parameters. The performance of the algorithm i s tested by computer simulation of a s i m p l i f i e d model for several cases of deterministic and/or stochastic input. I t i s demonstrated, that good estimates can be found from the output of a disturbed system when no i n t e n t i o n a l input i s applied. ii TABLE OF CONTENTS Abstract i i Acknowledgement v Nomenclature 1. 2. 3. vi INTRODUCTION 1 1.1 Dynamic Power System Models 1 1.2 Parametric Models for I d e n t i f i c a t i o n 2 1.3 Input Signal Requirements 3 1.4 Sources of Input Noise i n Power Systems 6 1.5 Methods for Parameter Estimation 6 THE MAXIMUM LIKELIHOOD ALGORITHM FOR PARAMETER ESTIMATION 9 2.1 State Space Description and Likelihood Function 9 2.2 The Continuous- Discrete Kalman F i l t e r 2.3 Steady State of the Kalman F i l t e r 14 2.4 Maximum Likelihood Estimates 17 2.5 Least Squares Algorithms 19 2.6 Optimization Methods 21 2.7 Correlated Input Noise 24 - 12 DYNAMIC POWER SYSTEM MODEL 26 3.1 Structure of the Model 26 3.2 Derivation of the State Equations 28 3.3 Dynamic Equations f o r one Machine 29 3.4 Load Bus Equations 35 3.5 Simplified Dynamic Model 36 3.6 Steady State Values 37 iii 4. 5. IMPLEMENTATION OF ALGORITHM 39 4.1 Calculation of Likelihood Values 39 4.2 Optimization Procedure 40 4.3 Computation E f f i c i e n c y 40 4.4 Simulation Programs 42 RESULTS 43 5.1 Parameters and Input Data 43 5.2 Output Data 44 5.3 Convergence of Parameter Estimates 50 6. CONCLUSIONS 58 7. BIBLIOGRAPHY 60 iv ACKNOWLEDGEMENT I would l i k e to thank my thesis supervisor, Dr. M. D. Wvong, as well as Dr. M. Davies f o r t h e i r guidance, advice and encouragement during the research work and w r i t i n g of t h i s t h e s i s . My thanks also to Ms. S. Y. Kim f o r typing the t h e s i s . My stay i n Canada was made possible by an award of the Canada Council. This grant and the f i n a n c i a l support of the University of B r i t i s h Columbia are g r a t e f u l l y acknowledged. v ' NOMENCLATURE General: small l e t t e r , eg. x or x column- vector T x row- vector c a p i t a l l e t t e r , eg. A matrix or s c a l a r quantity of power system T A E[x], x transpose of matrix A mean of random vector x T E[x.x ] covariance matrix of random vector x 6(t) Dirac delta 6,, kl Kronecker d e l t a L, V l i k e l i h o o d function (scalar) Time: t continuous time t^, k discrete time, t ^ = k.At State space description: x state vector y (model) output z (actually measured) output u control input vector w input noise vector v measurement noise vector a vector of unknown p h y s i c a l parameters A, B, C, H system matrices of continuous- F, G, H system matrices of d i s c r e t e - time model Q, R noise covariance vector matrices vi vector time model Kalman f i l t e r : x(-) , x(+) extrapolated and updated state estimate P(-) , P(+) covariance matrix of x ( - ) , x(+) respectively e innovations vector P e innovations covariance matrix K Kalman gains Power system quantities: a l l quantities except <5 are per unit quantities E, V voltage I current P, Q r e a l and reactive power X reactance X' transient reactance T torque or time constant T' transient time constant M machine i n e r t i a D macine damping 6 angle ( i n rad) between machine q- axis and synchronous frame of reference a angle ( i n rad) between machine q- axis and bus voltage Aw speed deviation from synchronous speed vii Indices: e electrical m mechanical d d- axis q q- axis i l o c a l generator i j external macine j t terminal of l o c a l generat L load FD voltage G speed governor bus regulator viii 1. I 1.1 INTRODUCTION Dynamic Power System Models There i s continued interest to improve dynamic and transient s t a b i l i t y l i m i t s i n e l e c t r i c power systems. new methods to achieve this goal. great promise [1]. power system. Modern control theory o f f e r s Linear optimal control designs show These techniques i n turn require a good model for the The interest here i s a model f o r a single generator to be controlled which i s connected to a large system. The c l a s s i c a l approach i s to model the generator as a machine connected to an i n f i n i t e bus [ 2 ] , This implies that the machine under consideration has no influence upon the dynamic behaviour of the rest of the power system. This assumption i s too r e s t r i c t i v e , e s p e c i a l l y for a large generator. On the other hand, each machine i n the entire interconnected system may be modelled i n d i v i d u a l l y . The r e s u l t i n g model i s then much too large and must be s i m p l i f i e d to become tractable. An intermediate approach i s to model the generator under consideration (the " i n t e r n a l " system) i n some d e t a i l and postulate a simple model for the remaining "external" system. Parameters of this simple model are then i d e n t i f i e d from measurements taken at the " i n t e r n a l " system generator. This method has two advantages: F i r s t , only the dynamics of the external system that interact strongly with the i n t e r n a l system heed be represented, and second, the parameters could be estimated by an on-line computer and the model therefore adapted to the actual state of the power network. 2. The proposed external system usually consists of one (or several) very large synchronous machines [3-7]. Yu, et a l . [ 3 ] used a small d i s t u r b - ing input s i g n a l and estimated the parameters from the system response. It was assumed that there were no disturbances i n the external system. In the following, additional random fluctuations i n the external system are considered. The problems treated are the modelling of these disturbances and estimating of parameters i n a noisy system. 1.2 Parametric Models for I d e n t i f i c a t i o n o _ Basic properties of i d e n t i f i c a t i o n problems are treated by Astrom and Eykhoff i n [8]. In the following i t i s assumed that a mathematical model for a power system i s derived from a p r i o r i knowledge and the i d e n t i f i c a t i o n problem i s reduced to the estimation of several unknown parameters. State- space models have been chosen for several reasons: they are e a s i l y obtained from d i f f e r e n t i a l equations, can describe multi-output systems e a s i l y , and form the basis for optimal control. (a) Deterministic models I f a l l inputs and outputs of the system can be measured accurately, a general l i n e a r state space model i s of the form: x(k+l) = F(a) . x(k) + G(a) . u(k) (1-1) y_(k) = H(a) . x (k) (1-2) a_ i s the vector of unknown parameters and s h a l l not depend on time. The discrete-time form of the state equation i s chosen because the input and output signals s h a l l be sampled and processed by a discrete 3. machine (computer). I f a c o n t i n u o u s - t i m e model x(t) = A(o).x(t) + B(a).u(t) i s derived (1-3) from d i f f e r e n t i a l e q u a t i o n s , i t can be e x a c t l y form (1-1) by i n t e g r a t i o n over one time s t e p , changes o n l y a t the d i s c r e t e time p o i n t s . converted provided that This into the i n p u t u_(t) i s not a severe restriction f o r a power system. (b) Stochastic The and models system may have i n p u t s ( d i s t u r b a n c e s ) which cannot be measured the measurements o f the o u t p u t s may be i n a c c u r a t e . A s t a t e space model f o r a system w i t h a d d i t i v e n o i s e i s : v(k) x ( k + l ) = F ( o ) . x ( k ) + G(a) . u(k) + M(a) . w(k) d-4) y_(k) d-5) = H(a). x(k) + v(k) and w(k) a r e random v a r i a b l e s of input and measurement n o i s e respectively. A d d i t i o n a l assumptions l i k e zero mean, independent sequences a r e u s u a l l y made. I t s h o u l d be noted t h a t variables t o o , i . e . , x.(k) and y_(k) the model s t a t e and output a r e now random are described by a p r o b a b i l i t y density function. 1.3 Input S i g n a l Requirements L e t x(0)=0_. Then i d e n t i f i c a t i o n i s o n l y requirements about the i n p u t condition be a r e met. possible F o r many e s t i m a t o r s , o i f some m i n i m a l a sufficient o f " p e r s i s t e n t e x c i t a t i o n " was o u t l i n e d by Astrom i n interpreted that the i n p u t [ 8 ] , I t may s i g n a l must have some amplitude f o r e v e r y frequency i n t h e band o f i n t e r e s t . Stochastic models have two types of input s i g n a l s , the control _u(k) which i s deterministic and the noise w(k) (a) Deterministic inputs These include disturbances Such disturbances which i s stochastic. applied for the sake of estimation. are a compromise between what i s desired for estimation purposes and what can be r e a l i z e d under physical and operational constraints. For estimation of power system dynamics input signals can be applied at the governor or the voltage regulator. Since any interference into normal opera- tion i s undesirable, these signals should be small and of short Yu, et a l . duration. [3] proposed a step change of mechanical torque for a short time, the simplest s i g n a l with "persistent e x c i t a t i o n " . (b) Noise inputs Observed changes i n the system under t e s t , which are neither to applied inputs nor to parameter changes may stochastic inputs w. These inputs may due be attributed to a d d i t i o n a l be considered a nuisance since they make a c o r r e l a t i o n of applied input and measured output much more d i f f i c u l t . But on the other hand, one might ask i f the noise alone disturbs the system enough so that i d e n t i f i c a t i o n i s possible from output measurements only. A white or colored Gaussian noise sequence for example would meet the requirement of persistent e x c i t a t i o n . The l e s s i s known about these noise inputs, the system i d e n t i f i c a t i o n gets more d i f f i c u l t , since not only system parameters the , but also parameters of the postulated noise process have to be estimated. But t h i s i s compensated by the advantage that the estimation process requires no interference of the normal operation of the system at a l l , an advantage that would count very high for power system identification. In fact, a few attempts have been made to estimate power system dynamics from normal operating data. Lindahl and Lyung [4] t r i e d to f i t several models by processing actual power system measurements and estimated parameters of a very simple 5th order model. They achieved not unreasonable estimates but expressed a desire to introduce a d d i t i o n a l i n t e n t i o n a l d i s t u r b ances. Price et a l . [7] tuned simulation programs to produce output signals that were very s i m i l a r to actual f i e l d measurements. They concluded that parameter estimation from on-line measurements i s possible, but that a choice of several model structures should be a v a i l a b l e , based upon knowledge from system planning data. (c) E f f e c t s of noise and feedback One has to be careful i n choosing an input s i g n a l that i t i s not the r e s u l t of a feedback. An example of what can happen i s given i n [8] f o r a transfer function model U G(z) 1—6 noise w H(z) An attempt to i d e n t i f y G(z) from measurements would y i e l d the estimate Y(z) _ 1 u(z) H(z) ' For a power system model the e x c i t a t i o n voltage and the mechanical torque can be treated as input signals only i f the feedback loops of the governor and voltage regulator are open. 6. 1.4 Sources of Input Noise i n Power Systems F i e l d measurements of t e r m i n a l q u a n t i t i e s on a g e n e r a t o r show f l u c t u a t i n g s i g n a l s even i f no i n t e n t i o n a l disturbance (e.g. i s applied. These f l u c t u a t i o n s come from changes i n the i n t e r n a l or e x t e r n a l system the model o f the e n t i r e system s h o u l d take them i n t o account. o n l y f l u c t u a t i o n s i n the e x t e r n a l system power demand w i l l be The s i m p l e s t way [7]) In t h i s and thesis, considered. t o account f o r these demand changes i s to model the a c t i v e and r e a c t i v e power demand as s t o c h a s t i c p r o c e s s e s w i t h the steady power as mean v a l u e . state Hence a l l l o a d dynamics s h o u l d be i n c l u d e d i n the one- machine model t h a t d e s c r i b e s the dynamics o f the e x t e r n a l system. P r i c e e t a l . [5-7] i n t r o d u c e d power demand as a s t o c h a s t i c i n p u t ) i n t h e i r model, whereas L i n d a h l } et al.[4] introduced a general noise vector w w i t h unknown c o v a r i a n c e . The e a s i e s t noise process w h i t e n o i s e and from an e s t i m a t i o n s t a n d p o i n t Gaussian t h i s r a t h e r r e s t r i c t i v e assumption i s made f o r most of following chapters. How a somewhat more r e a l i s t i c n o i s e sequence can i n t r o d u c e d i n the power system model i s shown i n s e c t i o n 1.5 is be 2-7. Methods f o r Parameter E s t i m a t i o n There a r e a v a r i e t y o f methods a v a i l a b l e , the problem i s to one which y i e l d s an u n b i a s e d , computing (a) the reliable, s t a b l e estimate with find reasonable effort. Parameter e s t i m a t i o n as s t a t e e s t i m a t i o n F o r s t a t e e s t i m a t i o n i n l i n e a r systems, the Kalman f i l t e r v e r y e f f i c i e n t and w e l l proven e s t i m a t o r . The parameters a_ i n eq. is a (1-4) 7. can be considered as system states which obey the equation cx(k+l) = tx(k) . The parameter estimation problem can therefore be viewed as estimation of the states of the augmented model [9]: ' x(k+l)' " F(a(k)) 0 " ' 2L( )' k * a(k+l) 0 I a(k) This method has two drawbacks: The augmented system becomes nonlinear even i f the o r i g i n a l system was l i n e a r and the order becomes higher too. It i s not too well suited f o r estimation of unknown parameters [10] , but i t should be considered for parameter tracking i f very good i n i t i a l estimates are available. (b) Least squares methods Least squares methods involve minimization quadratic function of some error. of a cost which i s a The error r e s u l t s from comparison of the model behaviour with the measurements on the r e a l system. The error i s normally so defined that i t i s l i n e a r i n the parameters. The cost function can be minimized i n one step. "equation error" methods y i e l d e f f i c i e n t algorithms mented as recursive estimators. Unfortunately, biased estimates when noise i s present These so c a l l e d and are e a s i l y imple- they can produce strongly i n the system [8]. In "output e r r o r " methods, the error i s the difference between measured output and model output. This error i s not l i n e a r i n the parameters and nonlinear techniques are applied to optimize method was used by the cost function. Yu, et a l . [3] for a deterministic system. This I t can be generalized for the case with noise and i s related to the maximum l i k e l i h o o d method (see section 2-5). (c) Maximum l i k e l i h o o d method The maximum l i k e l i h o o d method i s designed for stochastic system models where the outputs are random variables. model parameters i n a way that maximizes delivers the measured output. I t b a s i c a l l y adjusts the the p r o b a b i l i t y that the model I t i s reported to y i e l d stable, unbiased estimates [8] and has been shown useful for power system estimation [4-7, 11]. The maximum l i k e l i h o o d method works f o r systems with several outputs and also for systems with only noise inputs. Although i t i s numerically expensive, i t can be implemented on a modern large computer. 2. THE MAXIMUM LIKELIHOOD ALGORITHM FOR PARAMETER ESTIMATION In this chapter, the maximum l i k e l i h o o d algorithm for parameter estimation i s developed. optimal f i l t e r The major component of the estimator i s a l i n e a r (Kalman f i l t e r ) . The f i l t e r presented processing of sampled measurements from a continuous i s useful for d i g i t a l process. of a general white input noise, the model and f i l t e r should be In the case rearranged with the f i l t e r gain as unknown stochastic parameters to insure uniqueness. Alternatively, a n input-noise model of low order may be preferable. A comparison with a least squares algorithm i s made. for nonlinear function optimization are presented. Some methods The system model i s extended to include correlated noise input. 2.1 State Space Description and Likelihood Function Let the physical system be described by the following equations: x(t) y(t ) k = A(ct) . x(t) + B(a) = H(a) . u(t) + C(a) . ( t ) w . x(t ) +v(t ) k (2-1) (2-2) k with the vectors: x: state w: u: control input plant or input noise y: model output v: measurement noise a: unknown physical parameters The mixed continuous (2-1) and d i s c r e t e (2-2) form indicates that the model i s derived from d i f f e r e n t i a l equations but measurements are sampled and processed i n d i s c r e t e time, w i s a Gaussian white noise process with 10. zero mean and constant covariance Q, i . e . : E[w(t)] = 0 E[w(t). W (T)] T = Q . 6(t-r) (2-3) v i s a Gaussian white noise vector with zero mean and constant covariance R, i.e.: E[v(k)] = 0 E[v(k) . v ( * ) ] = R . 6 (2-4) T k A where the index k denotes time t. . k In addition to a , Q and R may be unknown too. With these assumptions for v and w, the p r o b a b i l i t y density function of the output pr(y(k)) can be calculated for a l l k as a function of a , Q and R from Eqs. (2-1, 2-2). The l i k e l i h o o d function i s the value of the j o i n t p r o b a b i l i t y density for the actually measured output sequence pd ' ,„ y(N),y(N-l),...,y(l) . . . , z ( l ) | a,Q,R) ments. (z(N) , z(N--l) ,. where z ( l ) ....z(N) denote the a v a i l a b l e output measure- The maximum l i k e l i h o o d method attempts to find those parameter t*. estimates a, Q, R for which the l i k e l i h o o d function i s maximized. Calculation of the l o g l i k e l i h o o d function [10] The j o i n t p r o b a b i l i t y density function i s computed r e c u r s i v e l y by application of Bayes' r u l e . With the abbreviation Y(k) = [ y ( k ) , y ( k - l ) , . . . , y ( l ) ] andZ(k) - [z(k),z(k-l),...,z(1)](2-5) pd Y(k) (Z(k)|o,Q,R) = pd (k) (z(k)|z(k-l);a,Q,R) . pd Y(k _ ( Z ( k - l ) | a , Q , R ) (2-6) 1) If L denotes the negative logarithm of the l i k e l i h o o d function, Eq. (2-6) becomes 11. -L(k;a,Q,R) = -L(k-l;a,Q,R) + In pd (k) (z(k)jZ(k-l);a,Q,R) (2-7) For L(0)=0 and with the recursion carried out N L(N;a,Q,R) = - £ In pd k=l y V (z(k)|Z(k-l);a,Q,R) (2-8) ; Minimizing L w i l l maximize the l i k e l i h o o d , since the logarithm i s a monotonic function. I t remains to evaluate the c o n d i t i o n a l p r o b a b i l i t y density of y(k). Because of the Gaussian assumption for the noises w and v, this density i s also a normal d i s t r i b u t i o n . The normal d i s t r i b u t i o n for a random v a r i a b l e i s P where d (z) y e=z-y = (2.T7. C T 2 ) . exp(-e / a ) 2 _ i (2-9) 2 2 , y=E[y], a = E[(y-y) ] 2 2 The corresponding multivariable normal d i s t r i b u t i o n for a random vector i s P d y ( k ) . P (z(k)|z(k-l)) -:L y =[(2.u) . m det P ( k | k - l ) ] . e x [ - \ e ( k ) _ i y T P (k|k-l).e(k)] where m = dimension (2-10) of the vector y e(k) = z(k) - y(k|k-l) y(k|k-l) = E[y(k)| (k-l);a,Q,R] Z P (k|k-1) = E [ ( y ( k ) - y ( k ) ) ( y ( k ) - y ( k ) ) T y |z (k-1) ;a,Q,R] The conditional mean and the covariance of y(k) are provided by an optimal l i n e a r f i l t e r (Kalman f i l t e r ) that processes the measurements z ( l ) ... z(k-l) for the s p e c i f i e d values of a, Q, R. 12. 2.2 The Continuous-Discrete Kalman F i l t e r The Kalman f i l t e r for the model (2-1, 2-2) i s a combination of the discrete and continuous f i l t e r s described by Gelb [9]. (a) Propagation i n the absence of measurements Between t, .. and ti, no measurements are available and the state k-1 Kestimate x(t) and i t s error covariance P(t) are simply extrapolated: x (t) = A . x(t) + B . u(t) (t _ <t<t ) k 1 (2-11) k P (t) = A . P(t) + P ( t ) . A + C.Q.C (t _ <t<t ) T T k (b) 1 k (2-12) Update with new measurement At time t k a new measurement, z ( t ) , k update state estimate and covariances. point and the indices t and t k + k i s available and i s used to There i s a discontinuity at this (or simply k- and k+) are used to indicate the quantities before and a f t e r the updating takes place. y(k-) = H- x(k-) output estimate (2-13) innovations (2-14) innovations covariance (2-15) x(k+) = x(k-)+K(k)-e(k) updated state estimate (2-16) P(k+) updated error covariance (2-17) e(k) = z(k) - y(k-) P (k) = H • P(k-)*H e T + R = [I - K(k)«H] • P(k-) K(k) = P(k-)-H »[H.P(k-).H +R]~ T T 1 Kalman gains (2-18) y(k-) and P (k) are the required mean and covariance y(k|k-l) and P (k|k-1) e i n the l i k e l i h o o d formula (2-10). v 13. Substitution of (2-10) into (2-8) y i e l d s N L(N;cc,Q,R) = j £ m-TT+ln det P^(t) + e ( k ) - P (k) .e(k) ] T e (2-19) k=l The c a l c u l a t i o n of the l i k e l i h o o d value can be done by a s l i g h t l y modified Kalman f i l t e r . Solution of the D i f f e r e n t i a l equations (2-11, 2-12) I t i s possible to integrate both x ( t ) and P(t) numerically, for example by the Runge-Kutta excessive. method, but the computational burden would be Since the model parameters are a l l time-invariant, Eq. (2-11) and (2-12) may be converted into difference equations. x(t) = A-x(t) + B-u(t) t k + 1 <t<t , x ( t + ) given k k x (2-11) the superposition i n t e g r a l i s x ^k > •<v k-i -* i*-i = t ) ( t t i o n matrix $(x) = exp(A'x). ) $ (t -x)« B-u(x)dx with the transi- + I f the input u i s allowed to change i t s value only at the time points t ^ , k=l ... N (see section 1-2), then the constant term B.u(x)=B.u(t _^) can be taken out of the i n t e g r a l . k With At = t ^ - t ^ (2-20) F(a) = *(At,a) = exp(A(a)- At) G(a) = f At (2-21) -1 exp[A-(At-x)].dx-B = A (a)•[F(a)-I]-B(a) (2-22) 0 the integrated version of (2-11) i s x(t-) = F(a).x(t+_ ) + G(a).u(tHh_ ) 1 1 (2-23) For the coveriance equation (2-12) the same approach as f o r the Matrix • R i c a t t i equation [9, p.136 ] i s taken. T For the equation P = A-P+P.A +CQC (2-12) the transformations X = P«y and u = -A «y define a l i n e a r system r -A o" T • A X X with the solution r y(At)' (F'V 0 y (0) x (0) X(At) where T <F-V 1 r-A exp 0- * At A C.Q.C (2-24) and for P P(t~) = T • F T + F • P(t+ .) k-1' (2-25) There are many methods to evaluate the matrix exponentials numerically, see section 4-3. 2.3 Steady State of the Kalman F i l t e r The c a l c u l a t i o n of one l i k e l i h o o d value L(N;a, Q,R) as developed i n the previous sections involves two steps T (1) For given values of a, Q, R and P(0) = E[x(0)-x (0)], the system matrices F(a), G(a), H(a) , the Kalman gains K ( t ) , k=l...N and the k innovations covariance P ( t j ) , k=l...N are calculated. g c (2)1 The measurements u(0) ... u(N-l); z ( l ) .... z(N) are processed by the Kalman f i l t e r to y i e l d the innovation sequence e(l)...e(N) and L. Since both the system dynamics F, G, H and the noise s t a t i s t i c s Q, R are assumed independent of time, the Kalman f i l t e r w i l l reach a s t a t i s t i c a l steady state, i . e . , the matrices P ( t ) , P ( t - ) , P (t£)» ^ ^ k ^ k constant matrices P(+), P ( - ) , P , K. e converge to The estimation algorithm requires r e l a t i v e l y long data sequences (> 100 p o i n t s ) , the time dependent covariance and gain matrices can therefore be replaced by the steady-state matrices, with the additional advantage, that dependence upon the unknown i n i t i a l covariance P(0) i s removed. The steady state value could be obtained by solving the matrix equation P(t^) = ^(^-i^* However, this system of equations i s d i f f i c u l t to solve. nonlinear algebraic It i s numerically more convenient to calculate the Kalman gains as shown i n the f i r s t step above with an a r b i trary matrix P(0) and use the values to which the gains converge. For the power system model considered i n section 3-5 choice P(0)=0, convergence to 5 d i g i t s was and the usually obtained i n 5 to 10 steps, an i n d i c a t i o n that s t a t i s t i c a l steady-state can safely be assumed. Summary of Equations for the Likelihood Function System model (process): x(t) = A(a) • x(t)•+ B(a).u(t) + C(a)-w(t) (2-1) z(k) = H ( a ) - x ( t ) + v(k) (2-2) E[w(t)-w (x)] (2-3) k T = Q • 6(t-x) E[v(k)-v (£)] = R • 6^ T (2-4) Kalman F i l t e r : x(k-) > = F-x(k-l+) + G.u(k-1+) (2-23) e(k) = z(k) - H-x(k-) (2-13/14) x(k+) = x(k-) + K-e(k) (2-16) l i k e l i h o o d function (constant term mir neglected) : N L(N; a, Q,R) = ^ {N-£n det P + £ e (k)-P . e(k)} k=l T 1 (2-19) Steady state matrices: K= l i m K(k) , k-x» P = l i m P„(k) from , e e K P(k-) = F-P(k-1+).F + T-F T T ; P(0+) = 0 (2-25) P (k) = H-P(k-)-H + R (2-15) K(k) = P(k-)-H . P (k) (2-18) P(k+) = [I-K(k)-H].P(k-) (2-17) T e T _ 1 e Remark: This form of the Kalman f i l t e r i s v a l i d also i n the case of perfect measurements are (R=0), although numerical d i f f i c u l t i e s may a r i s e when inverses calculated. Continuous-discrete r e l a t i o n s h i p : F = exp(At-A(a)) (F"V G = A .(F _ 1 -A 0 " = exp At • F (2-21, I)-B C.Q.C 2-22) 0-T (2-24) A By use of the i d e n t i f y T T x • A • y = trace (A* y • x ) the time independent matrix P £ can be taken out of the summation i n the l i k e l i h o o d formula. Hence L(N;a,Q,R) = N-£n det P 2.4 + trace(P 1 N £ e(k)-e (k)) k=l (2-27) Maximum Likelihood Estimates The existence of optimal parameter estimates i s now investigated. Uniqueness of optimal parameters The major drawback of state space models i s t h e i r non-uniqueness; many combinations of the matrices A,B,C,H account for the same input-output sequence. [11]. For multivariable systems even canonical forms are not unique For a p h y s i c a l l y derived model with a small set of parameters a, this need not be a b i g problem. parameters, Care should be exercised i n the choice of so that the transfer functions H- (sI-A) ^"B and H-(sI-A) ^"C are unique. o _ Both Astrom [11] and Kashyap [12] point out, that i n general the noise s t a t i s t i c s (Q and R) cannot be determined by the l i k e l i h o o d method. D i f f e r e n t noise covariances may give the same values for K and P e and there- o fore the same l i k e l i h o o d L. function of a, K and P . e Astrom suggests that L be minimized as a The Innovations Representation [8,11] The equations of Kalman f i l t e r and system model can be combined to the following model: x(t~) = F ( a ) . x ( t ~ _ ) + G ( a ) - u ( t _ ) + K ' - e C t ^ ) (2-28) z ( t ) = H(a).£(t ) + e ( t ) (2-29) 1 k k k k 1 where K'=F.K and e(t^) , k=l..N i s a sequence of independent vectors with zero mean and covariance P . These equations can be solved f o r the innovae tions e i f a, K', u and z are given. The l i k e l i h o o d function i s n 1 L(N;a,K;P )= ^ {N-ln det P + trace(P e 6 6 The minimization with respect to P T . £ e(k)-e (k))} N £ (2-30) k=l i s done a n a l y t i c a l l y . g derivatives of L to each element of P ?e opt = I -1 Setting the to zero gives the condition N I e(k)-e (k) k=l (2-31) T The maximum l i k e l i h o o d algorithm i s then restated as: N minimize V(N;cc,K*) = det £ e(k)-e (k) k=l (2-32) e(k) = z(k) - H(a).x(k) (2-33) x(k) = F(a)-x(k-l) + G(a)-u(k-l) + K'e(k-l) (2-34) x(0) = 0, e(0) = 0 This i s very s i m i l a r to a minimum variance estimator; the r e l a t i o n s h i p i s shown i n the next section. Kashyap [12] provides a proof of convergence and the necessary conditions for estimation from an innovations representation. However, he uses a d i f f e r e n t type of input-output description to avoid the ambiguities of the state-space model. The conclusions of [12] can be interpreted i n the following way: Under the assumptions of - zero mean, Gaussian input noise with a c e r t a i n covariance - zero mean, independent Gaussian measurement noise 19. - stable feedback system and l i n e a r optimal f i l t e r ( i . e . eigenvalues of both F and F-K'H inside unit c i r c l e ) the maximum l i k e l i h o o d estimates of both the parameters i n F and H and the f i l t e r gains K w i l l converge to the true values as the number of samples N goes to i n f i n i t y (NB: This includes the case with no deterministic inputs). Choise of reference model For the t h e o r e t i c a l reasons mentioned, the innovations description should be used for estimation of a general model. This was done by Lindahl [4] , who made no s p e c i a l assumptions about the manner i n which the input noise enters the system. On the other hand, i f there i s a simple physical description of the source of the input noise, then the noise covariance Q may contain only a few unknown e n t r i e s . The simple power system model of section 3-5 contains only one stochastic parameter (the load demand variance), but at least s i x Kalman gains. In t h i s case, the o r i g i n a l function L(a,Q,R) should be optimized, since i t involves fewer parameters. Price et a l . [ 5 ] also choose this approach for their model. 2.5 Least Squares Algorithms (a) Single output systems For a p h y s i c a l l y derived model with measurement noise the output error method i s used. J(a) The cost function i s N• I e(k/a) , = k=l e(k/a) = z(k)-y(k/a) (2-35) The output error or r e s i d u a l e i s the difference between the measured output z and the output y of a reference model. Nonlinear optimization methods are used to minimize J ( a ) . Convergence properties of the estimates depend on the reference model that provides the y(k). I f the reference model i s of the form x(k) = F(a)-x(k-l) + G(a)-u(k-l) (2-36) y(k) = H(a)-x(k) (2-37) which does not include any noise, but the measurements are created by a process disturbed by noise at the input (w(k)), then the residuals e(k) w i l l be a correlated sequence. In this case the least squares estimate i s biased [8]. For transfer function models with noise, the least squares method i s generalized by p r e - f i l t e r i n g the measurements u(k) and z(k) i n such a way that the residuals e(k) become uncorrelated. i s accomplished by the Kalman f i l t e r . For a state space model, t h i s I f the innovations model (2-33, 2-34) i s used as reference instead of (2-36,2-37), then the residuals are the uncorrelat innovations. In f a c t , f o r a single output system the l i k e l i h o o d function (2-30) is N V(a,k) = det I e(k)e (k) k=l X N = £ e^(k)=J(a,K) k=l (2-38) In the case of a single output, the maximum l i k e l i h o o d estimate and the least squares estimate (for the innovation model) are accordingly i n d e n t i c a l . In addition, i t i s also the minimum variance estimate, unbiased approximation of the variance of e. since J(ot,K) i s an 21. (b) Multiple outputs The least squares cost function f o r several outputs i s T J(a) = I e (k/a)-W.e(k/a) ; k=l N e(k/a)=z(k)-y(k/a) (2-39) with an a r b i t r a r y weighting matrix W. I f the e(k) are independent Gaussian vectors (from an innovations model) then i t can be shown that the optimal choice f o r W i s the inverse of the covariance matrix and f o r this choice the least squares estimate i s equal to the minimum variance estimate with the cost function J'(a) = trace N N £ e(k)<e (k) = 1 1 e (k) k=l k=l i (2-40) 2 1 this i s d i f f e r e n t from the l i k e l i h o o d function which involves the of the covariance of e. determinant Kashyap [12] points out that the l i k e l i h o o d function i s the correct c r i t e r i o n and that there are cases where estimation with Eq. (2-40) does not y i e l d the true parameters. 2.6 Optimization Methods With a minimum f o r the l i k e l i h o o d function guaranteed and the necessary equations to compute function values a v a i l a b l e , the problem of how to find the minimum value s t i l l remains. A great number of nonlinear optimization methods e x i s t , but i t i s not easy to find one that i s both r e l i a b l e and e f f i c i e n t . A l l the methods discussed here are described i n Adby and Dempster [13]; some of them are available at UBC see Manual UBC NLP [14]. as subroutines, (a) Search methods These methods attempt to find an optimum by using function evaluations only. Problems a r i s i n g from computation of derivatives are avoided, but a large number of function evaluations are needed. A simple method i s the nonlinear simplex technique, a more sophisticated one i s Powell' method of conjugate d i r e c t i o n s. (b) Gradient methods Gradient methods are based on a f i r s t order Taylor expansion f o r the Objective function T V(a +Ao) = V(a )+gV O )Aa, o o 3T g = -~da (2-41) They are normally superior to search methods f o r functions with continuous derivatives. Examples are steepest descent and conjugate gradient methods. A second order Taylor expansion involving the Hessian matrix V(a°+Aa) = V(a°)+g .Aa+ -|Aa . H •Aa T (2-42) T leads to Newton's method with the update formula Aa = -H . -1 g (2-43) The d i r e c t use Newton's method i s l i m i t e d , because the Hessian i s usually very expensive to evaluate and the method i s prone to divergence f o r i n i t i a l values not near the optimum. quasi-Newton But there are many approximations, known as methods, that seek to approximate theHessian d e f i n i t e matrix using only derivatives of f i r s t order. by a p o s i t i v e - For least squares problems the Gauss-Newton method i s often applied, with several possible modifications. Another group consists of the variable metric methods, the best known of which i s the DFP (Davidon-Fletcher-Powell) algorithm. Bard [15] compared several of the gradient methods f o r several s t a t i c parameter estimation problems and found that the Gauss-Newton methods usually performed best. (c) Evaluation of gradients For the innovations model (2-31, 2-32) and the l i k e l i h o o d V'(a,K) = j £n V(a,k) = j In det S; S = £ e(k)-e (k) T function (2-44) k 1 the derivative to the p h y s i c a l parameter <x i s given by 9 * ' ^ i - I « < ! f 7 • S" ) 1 i and the s e n s i t i v i t y equations " I ^ k S" ..*) (2-45) 1 i 3e/3a ± = -3H/3a • x(k) + H . 3x(k)/3a 3x/3a ± =3F/3a -x(k-l)+F-3x(k-l)/3a +3G/3a -u(k-l) i i (2-46) ± i (2-47) The s e n s i t i v i t y matrices 3F/3o^ , 3G/3a^ are f i n a l l y obtained from 3A/3a^, 3B/aa . ± For a physical power system model as developed i n Chapter 3, the a n a l y t i c a l derivatives of the system matrices w.r.t. the parameters become very complicated. A computationally simpler approach i s to perturb each parameter i n turn by a small amount Aa^ and evaluate the gradient from the approximations 3V/3a % [V(a+Aa ) - V(a)]/Aa ± ± (2-48) Each perturbation A a ^ must be small enough that the l i n e a r approximation is. true, but large enough that roundoff errors are kept small. A gradient evaluation requires as many function evaluations as there are unknown parameters. Choice of an optimization method The t h e o r e t i c a l advantage of the gradient methods seems questionable i n view of the computational e f f o r t required to compute the gradient. The main reason f o r this i s that the model derived by physical p r i n c i p l e s tends to have only a small number of parameters, which enter the model i n a complex way. The option to experiment with several methods on the actual model should be l e f t open. Nonlinear At UBC this i s e a s i l y possible by working with the Optimization Monitor [14], which o f f e r s a choice of optimization methods. From the limited experience that was gained by working with the model developed i n section 3-5, i f seems that Powell's method of conjugate directions performs best f o r this type of model. 2.7 Correlated Input Noise The model proposed i n section 2-1 may be inadequate f o r two reasons: - The input noise w i s not white but band limited and would better be modelled as a f i r s t order Gauss-Markov process: w(t) = A -w(t) + s ( t ) , s ( t ) white noise (2-49) The input noise w may influence the output d i r e c t l y , so that Eq. i(2-2) i s replaced by (2-50) y(t) = H(a).x(t)+D(a)-w(t)+v(t) The r e s u l t i n g complications f o r the estimator can be resolved by including w i n an augmented state vector [6]: x(t) " x(t) A(<x) B(a)' + _w(t) 0 w . w(t) "0 (2-51) • u(t) + s(t). . 0 x(t ) k y ( t ) = [H(a) D(a)] + k w(t ) v ( t k ) (2-52) k This augmented stated space model has the same stochastic properties as the o r i g i n a l one. The elements of the matrix A r? may be considered a d d i t i o n a l unknown parameters or may be assigned empirical values. 3. DYNAMIC POWER SYSTEM MODEL A two machine equivalent with random load demand i s used as a dynamic power system model. A general method to obtain a l i n e a r i z e d state space model from a system of mixed algebraic and d i f f e r e n t i a l equations i s presented and applied for the two-machine case. the model i s used for t e s t i n g the estimation 3.1 A s i m p l i f i e d version of algorithm. Structure of the Model The .following model i s used for parameter estimation: _rmm_ *di' q i x -TYYVW Vt Fig. 3.1 local generator external system - l o c a l generator ( i ) : This i s the generator identified. The generator for which a dynamic bus model s h a l l be i s described as a s a l i e n t pole machine of 3rd order (state variables 6 , u^, E ^ ) with an a d d i t i o n a l voltage regulator ± and governor, each of f i r s t order (state variables E,.,., T . ) . A l l para' f d i mx b meters are known. Measurements can be made at this l o c a l generator They include the terminal quantities V , P , Q t t t only. and the rotor speed OK. - external system: The t r a d i t i o n a l concept of an i n f i n i t e bus i s replaced by a simple external system consisting of: - a large external round-rotor synchronous machine (j) which the dynamics of the external system. represents A t h i r d order description includes the equations of motion and armature reaction (state v a r i ables 6 ., oj. , E' . ) J J qj unknown parameters: M. - t o t a l external 3 inertia D_. - t o t a l external damping X_. - steady-state reactance of machine j X' - transient reactance 3 T* .- transient open c i r c u i t time - a bus load P + jQ . Li constant No dynamics are included for t h i s load, but Li the load demand i s assumed to vary randomly i n the following P = P + AP , where P i s the steady-state load and AP T way: is a Gaussian white-noise process. - the network interconnection, which consists of the unknown transmission l i n e reactance X bus. t between the machine terminal and the f i c t i t i o u s Only algebraic equations are considered. These are obtained from phasor diagrams for synchronous speed. Depending on the size of the load P i s operating as a motor or a generator. i s always employed. n , Q n , the external machine The sign convention for a generator 28. 3.2 Derivation of the State Equations The following approach was chosen because i t allows both s i m p l i f i - cation of the r e s u l t i n g model as w e l l as generalization for more machines. The desired model i s l i n e a r , v a l i d for small transient deviations from steady state. I t i s obtained by a f i r s t nonlinear equations about the steady state. order approximation of the The steady state i s evaluated separately, i t depends upon the unknown parameters (see section 3-6). From the physical model, two types of equations are obtained: (1) state equations of one machine with respect to the load bus 3-3): (6, of (section They can be w r i t t e n i n the form x=f(x,u,r) where x are the states a), E^ , . . . ) , u are the control inputs and r are intermediate variables the load bus (<5 , V ). The l i n e a r i z e d one machine model developed i n section 3-3 i s very s i m i l a r to the one of De Mello and Concordia [2] with two differences: - the voltage V of the load bus i s v a r i a b l e Li - the machine angle <5^ i s d i f f e r e n t from the angle between machine q-axis and bus voltage (a.=S.-<5 ) X X Li The output equations are of the form y=h(x,r,v) where v stands f o r the measurement noise. The l i n e a r i z e d equations of a l l machines are combined and one machine angle can be eliminated. (2) Algebraic equations f o r the connecting network (load bus, section 3-4): Energy balance conditions at the load bus are expressed by equations g(x,r,w)=0 where w stands f o r the random loads P , T Q. T In order to obtain a single state space model, the equations are solved for the intermediate equations. In the nonlinear which can be solved by variables r and substituted i n t o the state form, t h i s constitutes a load-flow i t e r a t i o n only. very simple model of section 3-5 problem However, i n the l i n e a r i z e d case Ag(Ax,Ar,Aw)=0 can be solved d i r e c t l y by Gaussian elimination. 3.3 the s u b s t i t u t i o n For the i s done a n a l y t i c a l l y . Dynamic Equations f o r One Machine - equations of motion M • Aw + D • Aw = T I (> = ID m o e (3-1) - Aw (3-2) - T a l l quantities i n per u n i t , except 6[rad] and W =120TT rad/sec o J- o w 2 M= — i n e r t i a constant base D damping constant T mechanical torque T e > ^P =e V,• d Id, +q V q«I e l e c t r i c a l torque ^ m Aco = to -w rot o w o . . • «_ • speed deviation 6 angle between rotor axis and synchronous frame of - e l e c t r i c a l equations assumptions: - 3 windings: g=0 d,q,f - resistance i n stator winding n e g l i g i b l e l u x changes i n nde gand negligible - fspeed deviation l i g iq-axis ble reference from Kimbark [16], p. 73 *d = M f * f " d ' d T q y *f q = L f X q f f' f " T d E L o = R f I . M f ' d x q • I + * f with X, = a) • L J , X d o d ' q = w *L_ , o q ' n T, = do R,. f „ (o M ' - w .M Q q o f f q L f r f X' = X, - | - M " d d 2 ff f FD R f f the following equations are obtained: T do V V • ^q d = X q = E q q + E q • I = E FD ' h q - X,«I, = E' - X'-I, d d q d d <"> 3 3 (3-4) (3-5) v In the following a l i n e a r i z e d model with the load bus as reference i s developed. 31. Power transfer between two buses: Fig. 3.2 This s i m p l i f i e d phasor diagram of F i g . 3*2 i s used i n the following. Equations f o r both machines ( i ) and (j) are obtained through the s u b s t i t u tions : Machine ( i ) : E=E \ , V=V ( L , X =X^. + X d a=a.=6.-6 , I = I L Machine ( i ) : d d ± t , X =X .+ X q q t , , I =I q E=E . , V=V , X =X! , X =X. qj L ' d 3 ' q 3 1 T a=a.=6.-6 , I =1 ,. , I =1 . J J L d dj q qj T incoming power at the load bus V: S = P+jQ = V- I = <V +jV ) . ( I - j I ) d d (3-6) (3-7) Q - ~ d-Iq + V . I (3-8) V q d V,= V. sina = X .1 d q q V = V. cosa = E - X . I q d d 32. Substitution of I dJ and I yields: q V-E«sin a +, —1 V^ „ «sin2a . ' 1 XJ 2 q „ p _ 1_ 1 ? x Q _ = V-E- cos a n X x (3-9) d 1_ _ 1_ r •2 - ,T2 v^'Isin^a X d q X (3-10) d X d For small variations AV,AE,Aa: 9P 9V AP with . AV+ 3P 9E 3P ,AE+ — . .Aa da a (3-11) 1 9P _ V E - cosa 1 1 + V cos 2a 3a " XJ l q " di 2 X 3P 3E V* sina 3P 3V E« sing X, + V . s i n 2a X 1 ^ 1 V x q where V, E, a on the right hand side are steady state values. Equations for generator ( i ) - e l e c t r i c a l torque with the appropriate s u b s t i t u t i o n s : AT . £ ex 3P with K li ex 3a 2i 7 V .E -cosa L q l -xT.+x dx 3P K ei 3E'. qi (3-12) AP . = K -Aa.+ K 'AE' + K .-AV ex lx i 2x q /x L V t sina. i L' X'+X d t T 2 + V L • c o s 2 a i v v x -+ -t qx x X'+X„ dx t 33. ei 8 P K . E 71 " I v f and qi' i S i n a • - % ^ are equal to V + and " s i 2 «i I x q l + x t - x- x i + t j i n the paper of De Mello and Concordia [2] - i n t e r n a l voltage Eq.: X -X' AE . =AE' . + (X -X' ).AI,. = AE'. + r J - ~ T qx qx dx dx dx qx X +X (AE' . —A(V .cosai)) qx L ~ - or AE . = - r - • AE'. + K,..Aa.+ K .•AV qx qx 4x x 8x L (3-13) Q with = v dx t X..+X,. dx t 3i „ dx 4 i K = dx 3rrTx— ' v dx K 8i = „ i t dx dx X'.+X • dx t ~ s i n a C O S a i - terminal voltage V. 2 V t AV 2 = dt V V A V qt A V dt' dt V V t = fc + A ' W ' X . = X .AI . = J l • A(V-sina) dt qx qx X .+X L qx t v n AV 2 + n X =AE\ - X' .AI,. = , ,„ qt qx dx dx Y x + d i x t X* -AE', + , J. q di t 1 -A(V .cosa) X hence AV = K Aa .+ K, . . AE' . + K . . AV t 5x x 6x qx 9x L (3-14) with hf ¥ t • Xqx^ T t \ t • xdxp rt • V ^ ° i 34. V K T K - 6i . 9i X * di t q t V X' .4X t dt ai — ' = TrTT t qi . ' t s i n a , i + qt di "v" ' Y T T T t di ' t c o s a i - terminal power AP_ t = AP = AT . L ei (3-15) Again, the constants K„. through K,. are i d e n t i c a l with K„ through K, i n [2] Jl ox j b - voltage regulator T FD'^FDi + F D i " FD A E K < ref " V AV A ( 3 " 1 6 ) - governor and turbine T_ ' A T . + AT . = K_ . (A co -Ao .) G mi mi -G ref I AV • ref (3-17) and Aco . are input signals through which an i n t e n t i o n a l disturbance ref can be applied. Summary of equations for machine ( i ) A6 . =co 1 o . Aw . i M-J «Au). = T . - T . - D - A j j i ' I mi ei A X T - A T . = - A T .• + YL (W .-Aw.) G mi mi » ref i r T' ,-AE' = -AE dot T q FD- FDi AE = .+AE-. qi " FDi A E FDi + K FD ( A V ref- t> A V 35. where AT . = K,.•Aa. + K_.•AE'. + K_.•AV ei l i i 2x qi 7x L T AE . = r^=— • AE' . + K. .* Aa. + K ." AV qx. qx 4x x 8x L 0 T AV. = K . •Aa. + K,.•AE'. + K ..AV t 5x x 6x qx 9x L c rt Summary of equations for machine (j) A6. = w„ . Aw. M. • Aw. = -D. • Aco. - AT . 3 1 3 3 J e T' .AE'. = -AE . doj qj qj AT . = K * Aa. + K • AE' . + K_ . • AV ej lj J 2j qj 7j L AE . = ^ — qj AE' . + K..•Aa. + K .* AV qj 4j j 8j L a In the combined model for both machines, one machine angle can be eliminated, since only r e l a t i v e positions of the machine axes are of i n t e r e s t , A<5 = A<5_^-A6_., This reduces the order of the model by one. 3.4 Load Bus Equations Three more equations are required to eliminate AV^^\a^,Aa2 from the state equations. P TL Two r e l a t i o n s are from power continuity conditions: = P. . + P . Lj Lx T Q = Q .+ Q . L Lj Li and one i s an angular r e l a t i o n s h i p Ac$..-A<5 =Aa -Aa„ = A6 0 n T 36. hence K,.Aa. +K.,.Aa. +K-.AE . + K „ . A E ' . + (K-. ,+K_ . )-* AV lx x 1] ] 2x q i 2j qj 7x 7j L 1 T AQ = L A6 (3-18) = Aa^ - Ao^ The 3 equations can be solved for AV^Ao^jAa i s best done numerically. by Gaussian elimination. This The expressions are substituted i n the state equations of the machines, which i n e f f e c t w i l l modify the constants 3.5 .-K^. Simplified Dynamic Model Since the estimation algorithm posed some problems and i s expensive in computing time, a s i m p l i f i e d model was used: - the voltage change on the load bus AV i s assumed to be n e g l i g i b l e - the load bus i s assumed to keep the frequency of the external machine Aa ± ^ A6 Aa ^ 2 0 Under these assumptions, only the equations of motion of the machine (j) r e main relevant and i t s e l e c t r i c a l torque can be expressed as AT . = AP - AT . ej L ex T (3-19) State equations for s i m p l i f i e d model: (AOVJ-AW.) A6 = w o v (3-20) y 1 M. • Aw-4. = -D4 . A all +AT . -AT . i 1 mx ex (3-21) M.'AtD. = -D.-Au). -AT . 3 3 3 3 ej (3-22) 1 37, T . .AT . = -AT . + K^(Au> ,-AOJ,) G mi mi G ref (3-23) 1 T do' T FD AT , A E qi • A E FD " FDi = - A E A E FDi qi + K ( 3 FD ( A V (3 ref- V A . = K -A6 + K *AE'. ei l i 2i qi AE . = K..A6 + q AV. t AT A E i 4 1 K q 25) (3-27) i = K .A6 + K-.AE'. 5i 6i qi (3-28) C . = AP -AT . ei L ei T h i s i s a 6 t h o r d e r model ~ 2 4 ) (3-26) AE'. 3 i " (3-29) w i t h t h e 3 unknown parameters Mj , D j , X unknown c o v a r i a n c e o f t h e random i n p u t AP e s t i m a t i o n p r o c e d u r e a r e : speed ^ . Measurable q u a n t i t i e s , terminal voltage AV t > t and t h e f o r the t e r m i n a l power AP =AT . t 3.6 ei Steady State Values The steady state i s computed from the known system parameters (M^, ^ j i ' •••)» estimates of the unknown parameter state quantities V Q, P Q> Qj:0" t ^ e l n and the terminal steady d e x 0 i s omitted i n the following. Fig. 3.3 Steady state phasor diagram The equations (3-9, 3-10) applied for the terminal of machine ( i ) are: P. = V .E' . .sinijj. t qx l V .E .,..siniK t qdx x X . T di V .E ...cosiji. t qdi x X . Y + V o ^ . siniJj. cosil;. t x i X . qi V_ t X qi and V .V .cos(-* ) V .V .sin(-i|» ) t L t t L V - t t -P. = w i t h the s o l u t i o n s = arc tan *V r (3-29) a r c tan t qi p . x' E' . = „ . + V -cos*, qx V .sxnijj. t x t V .sinijij. ^ (3-30) x. (3-31) 1X . qi (3-32) t (3-33) X« dx > 4. IMPLEMENTATION OF ALGORITHM This chapter contains a short description of the computer implementation of the estimator and of some methods to reduce computing time. 4-1 Calculation of Likelihood Values The following flow-chart shows the major steps i n the c a l c u l a t i o n of one l i k e l i h o o d value: arguments a, Q, R + c a l c u l a t i o n of system matrices A(a) , B(a), C(a) equations of power system model, sections 3-5, 3-6 4conversion to discrete-time system A(a) ,B(a) C.Q.C Q = T-F method see section 4-3 F(a),G(a) T D 4c a l c u l a t i o n of steady-state Kalman gains K and (expected) innovations covariance P i t e r a t i v e equations, see section 2-3 e + processing of the measurement arrays u, z with the Kalman f i l t e r and summation of the innovations (sample) covariance S equations, see section 2-3 4l i k e l i h o o d value L(a, Q, R) The f i r s t of the steps above contains the physical model. following steps are independent of i t . for The The computation sequence i s simpler the case without input noise when there i s no Kalman f i l t e r . 40. 4-2 Optimization Procedure The Monitor for nonlinear optimization available at UBC [14] was employed to find parameter estimates. I t supports i n t e r a c t i v e processing and switching among d i f f e r e n t optimization methods. The following user routines were added: FUNCTION XDFUNC - C a l c u l a t i o n of a l i k e l i h o o d value as described i n section 4-1 SUBROUTINE INIT - an i n i t i a l i z a t i o n routine that must be c a l l e d f i r s t and asks for the necessary data. Most of the a v a i l a b l e algorithms i n [14] were t r i e d , but none of the methods based on gradient evaluation performed w e l l . A possible cause i s that the parameter steps taken f o r the difference approximation not suitable. of the d e r i v a t i v e were Good r e s u l t s were achieved by the conjugate d i r e c t i o n s algorithm (routine POWEL). 4-3 Computation E f f i c i e n c y The number of function evaluations i s high for any optimization algorithm. Since a function evaluation involves the s e t t i n g up and running of a l i n e a r f i l t e r , computation time i s long even for a small system. A considerable reduction i s achieved by using the following techniques. (a) F i l t e r equations and matrix operations The mathematical operations of the f i l t e r are simple but highly repetitive. Since m u l t i p l i c a t i o n s of a vector by a matrix are among the most frequent operations, a l l matrices are stored by rows, so that a. . becomes A(J,I) i n Fortran. by a single index only. In most subroutines, array elements are addressed A package of short matrix subroutines was written i n Fortran. order i n which the elements are stored i s chosen c a r e f u l l y . The These matrix routines run several times faster than t h e i r more general counterparts i n the system l i b r a r y . (b) Conversion of continuous to discrete time This can be done by numerical integration. I f x(0) i s chosen as the i - t h unit vector, then integration of the system x=Ax over one time step At w i l l y i e l d the i - t h column of the transmission matrix F. The number of integration steps required i s therefore equal to the order of the system. It was found that by using the Runge-Kutta method, up to 80% of the t o t a l computing time was spent on this integration. A much faster procedure i s obtained through approximation of the matrix exponential by i t s s e r i e s : At-2 F = exp(A-At) = I + A f A+ jf~ A+ 2 ... (4-1) In general this method may suffer from several numerical problems, such as roundoff errors and slow convergence. But for the model i n section 3-5, approximation by the f i r s t eight terms i s accurate to 6 d i g i t s and about ten times faster than numerical integration. The complete recursive algorithm for F, G and T i s obtained from the series expansion of Eqs. (2-21, 2-22, 2-24): F(n) = F(n-l) + • A(n) F(0)=I (4-2) G(n) = G(n-l) + • B(n) G(0)=0 (4-3) T(n) = T(n-l) + ^ r - • T(n) n! T(0)=0 (4-4) 42. A(n) = A(n-l) • A A(1)=A=A-At (4-5) B(n) = B(n-l) • B B(l)=B=B"At (4-6) T(n) = - T ( n - l ) - A (c) Achieved The one computing T +A(n-1)-f T(l)=f=C.Q-C .At T times e v a l u a t i o n of a l i k e l i h o o d v a l u e i n p u t , one output (4-7) and f o r a system o f 6th o r d e r , with 100 measurement p o i n t s r e q u i r e s the f o l l o w i n g approximate CPU-times ( f o r an Amdahl 470 computer): 24 msec f o r d e t e r m i n i s t i c case (no Kalman gain) 55 msec f o r the s t o c h a s t i c case A t y p i c a l e s t i m a t i o n run w i t h 6 i t e r a t i o n s of the c o n j u g a t e t i o n s a l g o r i t h m r e q u i r e s about 100 This i s s t i l l f u n c t i o n e v a l u a t i o n s o r 5 sec CPU-time. c o n s i d e r a b l e f o r such a s i m p l e model. F u r t h e r improvements a r e p o s s i b l e . ( e s p e c i a l l y the f i l t e r ) c o u l d p r o b a b l y precision. Most o f the c a l c u l a t i o n s be done w i t h s i n g l e i n s t e a d of double Of most importance would be a f a s t e r , r e l i a b l e method f o r com- p u t i n g the steady 4.4 s t a t e Kalman g a i n s . S i m u l a t i o n Programmes A l l d a t a used f o r the e s t i m a t i o n p r o c e s s was puter s i m u l a t i o n with the same model. same f o r e s t i m a t i o n and and direc- fulfills estimator simulation. obtained from a com- Many of the s u b r o u t i n e s were The s i m u l a t i o n data i s therefore i d e a l a l l assumptions of s t r u c t u r a l and n o i s e p r o p e r t i e s t h a t presupposes. the the 5 DATA AND RESULTS The simulation data and the performance of the maximum l i k e l i h o o d estimation method are shown f o r the following cases: Case 1 - deterministic input only Case 2 - combined deterministic and stochastic inputs Case 3 - stochastic input only Case 4 - stochastic input and measurement noise 5.1 Parameters and Input Data The following parameters were used i n simulation program: - sampling time step At=0.05 sec (a) Deterministic parameters ( a l l per unit) - i n t e r n a l generator: D.=10 M.=5 X K di = 1 FD = 5 T* =8 di X =0.6 X X • X .=0.3 d q K =10 ° G — external system (unknown parameters a) - transmission reactance: - unknown machine: M.=9 X =0.5 t D.=26 - steady state values: v to = 1 - 0 5 p to = 0 - 9 Q to = 0 - 3 V 0.5 (b) Deterministic input The control input i n case 1 and 2 i s a short pulse applied at the reference input of the governor/turbine loop. AW f(pu.) re 0.25 0 0.75 t(sec) 5s»- 2 -10r3 F i g . 5.1 (c) Deterministic Input Stochastic input A random disturbance AP^ was applied i n cases 2 to 4. noise with a variance Q=E[AP corresponds ]=10 4 was considered. Only white In d i s c r e t e time this to a random sequence with a standard deviation of a =/Q'At = m if Li _3 2.2'10 p.u. For cases 1 to 3, measurement noise i s neglected. _9 measurement noise with covariance R=10 simulated 5.2 In case 4, -3 (o =0.032-10 v - ) was added to the output Au)j_. Output Data The simulated outputs were AP s i s t s of a sequence of 100 points. t , AV t and Aco^. Each example con- Some examples for each case are shown in F i g . 5.2 to 5.9. Table 5.1 contains for a l l simulated sequences the mean values m, the standard deviation a, the maximum absolute amplitude and the expected |A| innovation standard deviation /p^T of a Kalman f i l t e r with the true parameters. Table 5-1: Simulation Results: Case 1 2 Signal 4 Fig. |A| -11.3 (23.9) 94.1 _ 5.2 AV, 0.19 (2.02) 7.66 - 5.3 Aco^ -0.43 (0.62) 1.84 - 5.4 t -8.78 (30.3) 85.7 same 5.5 AV 0.31 (4.10) 11.2 as - -0.33 (0.83) 2.46 case 3 5.6 2.53 21.8 42.7 3.33 5.7 0.388 - AP A p t t A(o^ 3 m AP t AV t 0.120 3.07 6.26 Aw^ 0.099 0.445 0.881 0.0268 5.8 Ao). 0.1011 0.447 0.884 0.0806 5.9 i A l l values x 10 per unit. 46. 49. 5.3 Convergence of parameter estimate Parameters were estimated from the simulated measurements f o r a l l four cases. In a l l examples the cost function was minimized by the conjugate directions algorithm (routine P0WEL[14]). Only single output sequences were considered. This i s s u f f i c i e n t i n the cases without measurement noise, because one output contains a l l available information about the system. Estimation from multiple output sequences was not attempted, because numerical problems were encountered i n the c a l c u l a t i o n of the Kalman gains. This should be further investigated for the case with measurement noise. Table 5-2 contains a summary of the obtained optimal estimates for a l l examples. F i g . 5-10 to 5-18 show how the estimates and are followed by comments for the i n d i v i d u a l cases. converged T a b l e 5-2 •H I p>4 to I cr o C •H ^2 W cd "4" m vO r~- oo m m mi m in m CM CO rH i I m Estimation Results m I rH | 00 ON ON ON ON ON ON ON ON m o o O rH rH | rH rH rH rH rH rH CN CO o vO CN + rH o o O • + • CM + U CM 4-> o • rH + cd CM O • + o T3 CD 60 U CD > • -d- X> H cd * o o 00 o rH * <a o m rH o o m vO m CO I VO + cd rH o o m m + S ON ON ON co in CN rH Cd e o m CN VD o m CN rH m C o in ON CM m rH vO rH + O CO CM + rH + r- CM rH O CO •H 00 oo o CN CM O CO o m en oo rH II co CO + + m o CO CO • 00 I i o + CM I O I CO I II a o i u 3 o o o 00 CM ON CM r- vO 00 ON ON oo oo 00 00 rH o CO CM vO oo 00 00 CN vO O rH CO m CN • II a o vO CM oo ON II CD CD O to CU C - H CO & m •rl 4J •—• U cd m in I m cd co o n e <u o =S= 4-1 - H vO o 00 • cd 4J CO •HQ) O CM CM CO CM <r 60 3 ft O 0) CO cd u P-i < PH < -H 3 < •H 3 < 4J < CO 4-1 < 3 < 3 < •H -H 3 < • rt • m m II II a 4-1 O m • c 3< a CO • •« m •H CO M o •—s <—s oo o o m o- CO o o o- m O N-^ + I rH II o o rH v-'t\ 1 o co • CO (X O CM o o rH II Pi o* m rH oo o I vo 1O •H 4-1 cd cj •K •H *a OJ CO •H T3 1 rH O a CM O CO u CD 4-1 CD 6 cd u cd Pi d) 3 M 4-1 CD Xi 4J • • S3 F i g . 5-10 Case 1 - deterministic output AP F i g . 5-11 Case 2 - mixed output AP input oc (p.u) F i g . 5-12 Case 2 - mixed input output AID. F i g . 5-13 Case 2 - mixed input output Aw F i g . 5-17 Case 4 - w i t h measurement n o i s e output Aw., F i g . 5-18 Case 4 - w i t h measurement n o i s e output Aw. 56. (a) Case 1 - deterministic In this s p e c i a l case there i s no Kalman f i l t e r . The estimator i s the ordinary least squares estimator with the cost function V = I e (k). k 2 Rapid convergence to the true parameters was 5-10) and Aw^ measurements. obtained for A P ( F i g . t The algorithm f a i l e d for the AV t measurements, because the i n i t i a l steps were too large. (b) Case 2 - deterministic and stochastic input The same cost function was of a Kalman f i l t e r . used, but now the e(k) are the innovations For n e g l i g i b l e measurement noise (R-K)), the Kalman gains depend on the parameters a only and the expected innovations covariance i s proportional to Q. Q i s therefore estimated P £ from the sample covariance P =V/N. e Speed of convergence i s comparable to the deterministic case, but the estimates have some bias ( F i g . 5-11, 5-12). The f i r s t step of the algorithm always seems to make M. and D. l a r g e r , even i f they are J too large (Fig. 5-13). already J I t seems preferable to s t a r t with estimates that are smaller than expected. (c) Case 3 - stochastic input only Estimation with the cost function V= £ e (k) did not bring s a t i s f a c t o r y 2 results. S l i g h t l y biased estimates of M^. and D_. were found i f X constant at the true value 0.5. t was held But this i s only one l o c a l minimum of V and lower minima can be found for quite d i f f e r e n t values o f M, , and X. fc 57. Good r e s u l t s were obtained with the o r i g i n a l l i k e l i h o o d cost function L=£n det P + £ e ( k ) / P , where P 2 e filter. e e i s the estimated covariance of the Kalman Estimates from a l l of the three outputs converged (Figs. 5-14 to 5-17), with a small bias for M. and 3 (d) t and a rather large bias for D.. J Case 4 - input and measurement noise In the f i r s t example ( F i g . 5-l7) Q and R are held constant at the true values and only a i s optimized. The r e s u l t s are s i m i l a r as i n case 3. In the second example ( F i g . 5-18) only R i s f i x e d , a and Q are estimated simultaneously. converge reasonably. After four i t e r a t i o n steps, the estimates seem to However, i n the next i t e r a t i o n step very d i f f e r e n t estimates were found, for s l i g h t l y lower function value. For t h i s short measurement sequence, the algorithm w i l l not converge to the true values. 58. 6 CONCLUSIONS The maximum l i k e l i h o o d method i s a general method f o r estimation of parameters i n a l i n e a r system disturbed by both input and measurement noise. The stochastic parameters can be modelled either as unknown gains of a l i n e a r f i l t e r or as parameters of a low order model derived from physical knowledge. The second approach i s given preference and a low order model for the remote part of a power system i s proposed that explains small f l u c t u a tions by varying load demand. Estimation r e s u l t s were obtained for a s i m p l i f i e d system from one output at a time. They show that the maximum l i k e l i h o o d method can i n fact produce good estimates for a heavily disturbed system. Parameter estimates from a system with noise input only are also possible. However, convergence problems may arise for short measurement sequences. False convergence to a wrong set of parameters i s possible, e s p e c i a l l y i f the measurement noise i s appreciable. This case would need further investigation with longer data sequences and a more r e a l i s t i c power system model. Estimation of the dynamics of a r e a l power system from measurements alone should be possible i f ... - accurate measurements of small fluctuations are a v a i l a b l e . - the dynamic behaviour and the source of fluctuations of the external system with respect to a generator can be represented with s u f f i c i e n t accuracy by a very low order model. To answer t h e s e q u e s t i o n s , a c t u a l measurements s h o u l d be i n v e s t i g a t e d . from a generator 60. 7 1. BIBLIOGRAPHY Y.N. Yu and H.A. Moussa, " O p t i m a l s t a b i l i z a t i o n o f a m u l t i - m a c h i n e system", IEEE T r a n s , on PAS, v o l . PAS-91, pp. 1174-82 (May 1972). 2. F.P. DeMello and C. C o n c o r d i a , "Concepts o f synchronous machine s t a b i l i t y as a f f e c t e d by e x c i t a t i o n c o n t r o l " , IEEE T r a n s , on PAS, v o l . PAS-88, pp. 316-329 ( A p r i l 3. Y.N. Yu, M.A. 1969). E l - S h a r k a w i and M.D. Wvong, " E s t i m a t i o n o f unknown large power system dynamics", IEEE T r a n s , on PAS, v o l . PAS-98, pp. 279-289 (January 1979). 4. S. L i n d a h l and L. L j u n g , " E s t i m a t i o n o f power g e n e r a t o r dynamics from normal o p e r a t i n g d a t a " , 3rd IFAC Symposium, The Hague (June 1972), P a r t I , pp. 367-374, Paper pp 2. 5. W. P r i c e , F.C. Schweppe, E.M. G u l a c h e n s k i and R.F. S i l v a , "Maximum l i k e l i h o o d i d e n t i f i c a t i o n o f power system dynamic e q u i v a l e n t s " , P r o c e e d i n g s 1974 IEEE Conference on D e c i s i o n and C o n t r o l , pp. 579-586. 6. R.D. M a s i e l l o and F.C. Schweppe, " M u l t i - m a c h i n e e x c i t a t i o n v i a system i d e n t i f i c a t i o n " , 454, 7. stabilization IEEE Trans, on PAS, v o l . PAS-94, pp. 444- (March 1975). W. P r i c e , D.N. Ewart, E.M. G u l a c h e n s k i and R.F. S i l v a , "Dynamic e q u i v a - l e n t s from o n - l i n e measurements", IEEE T r a n s , on PAS, v o l . PAS-94, pp. 1349-1357, ( J u l y 1975). o 8. K . J . Astrom and P. E y k h o f f , "System i d e n t i f i c a t i o n - a survey", Auto- m a t i c a , v o l . 7 (1971), pp. 123-161. 9. A. Gelb ( e d . ) , " A p p l i e d o p t i m a l e s t i m a t i o n " , The A n a l y t i c Corporation ( M I T - P r e s s ) , 1974, 74-1604. Sciences 10. F.C. Schweppe, "Uncertain dynamic systems", Prentice H a l l , 1973, 7211801. o _ 11. K.J. Astrom, "Modelling and i d e n t i f i c a t i o n of power system components" from E. Handschin (ed.), "Real time control of e l e c t r i c power systems" E l s e v i e r Publishing Company, 1972.' 12. R.L. Kashyap, "Maximum l i k e l i h o o d i d e n t i f i c a t i o n of stochastic l i n e a r systems", IEEE Trans, on AC, v o l . AC-15, pp. 25-34 (Feb. 1970). 13. P.R. Adby and M.A.H. Dempster, "Introduction to optimization methods", Chapman and H a l l mathematic s e r i e s , 1974, 74-4109. 14. M. Patterson, "UBC-NLP, Nonlinear Centre, Univ. of B.C., function optimization", Computing 1978. 15. Y. Bard, "Comparison of gradient methods for the s o l u t i o n of nonlinear parameter estimation problems", SIAM Journal Numer. Anal., v o l . 7, pp. 157-186, (March 1970). 16. E.W. Kimbark, "Power system s t a b i l i t y : Synchronous machines", Dover, 1968, 68-12937.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Parameter estimation applied to power system models...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Parameter estimation applied to power system models subject to random disturbances Muller, Hans 1979
pdf
Page Metadata
Item Metadata
Title | Parameter estimation applied to power system models subject to random disturbances |
Creator |
Muller, Hans |
Date Issued | 1979 |
Description | For control of an electric generator, it is desirable to approximate the remote part of the power system by a low order model, whose parameters can be estimated from measurements. The proposed model for the remote system consists of a large synchronous machine and a random varying load to account for the dynamic behaviour and for the small fluctuations that are always present. The maximum likelihood algorithm is a good method for parameter estimation of a noisy system. It can even be used for estimation from measurements alone, without applied disturbance, which could be of great advantage in a power system. The algorithm is derived in a form suitable for efficient digital processing of sampled measurements from a linear continuous system with physical parameters. The performance of the algorithm is tested by computer simulation of a simplified model for several cases of deterministic and/or stochastic input. It is demonstrated, that good estimates can be found from the output of a disturbed system when no intentional input is applied. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0094649 |
URI | http://hdl.handle.net/2429/21569 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Unknown |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1979_A7 M94.pdf [ 2.83MB ]
- Metadata
- JSON: 831-1.0094649.json
- JSON-LD: 831-1.0094649-ld.json
- RDF/XML (Pretty): 831-1.0094649-rdf.xml
- RDF/JSON: 831-1.0094649-rdf.json
- Turtle: 831-1.0094649-turtle.txt
- N-Triples: 831-1.0094649-rdf-ntriples.txt
- Original Record: 831-1.0094649-source.json
- Full Text
- 831-1.0094649-fulltext.txt
- Citation
- 831-1.0094649.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0094649/manifest