Tunable Difference Equations for the Time-Domain Simulation of Power System Operating States by JESUS C A L V I N O - F R A G A Ingeniero Electronico, Universidad Simon Bolivar, Venezuela, 1987 M . S c , The University o f British Columbia, 1999 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF THE R E Q U I R E M E N T S FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department o f Electrical and Computer Engineering) W e accept this thesis as conforming to the required standard t T H E UNIVERSITY OF BRITISH C O L U M B I A July 2003 © Jesus Calvino-Fraga, 2003 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the requirements f o r an a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e h e a d o f my d e p a r t m e n t o r by h i s o r h e r r e p r e s e n t a t i v e s . It i s understood that copying o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Jesus Department o f Electrical The U n i v e r s i t y o f B r i t i s h V a n c o u v e r , Canada Date J u l y , 2.A 2003 Calvino-Fraga and Computer E n g i n e e r i n g Columbia Abstract This report presents the development o f a new family o f difference equations for the numerical solution o f systems o f differential equations arising in electric power circuits. The main advantage o f these difference equations is their adjustability, which allows for the tuning o f the resulting formulas while using large time steps. The tuned formulas have been used successfully in the time-domain simulation o f 50/60-Hz solutions o f electric power circuits for linear and non-linear problems. The linear steady state and the non-linear power flow problems can be solved, with no error either in magnitude or in phase, at a comparable computational cost to traditional techniques such as the well-known phasor analysis and Newton-Raphson iterations. Additionally, the techniques developed do not use complex numbers, resulting in a faster matrix factorization for linear systems, or a smaller Jacobian matrix for non-linear systems. The main advantage ofthe proposed techniques is that due to the time-domain nature o f the simulation, where the transition from one time solution to the next represents a small transition between states, the technique is able to closely follow the trajectory o f the system even for ill-conditioned power flow problems where traditional Newton-Raphson iterations fail to converge, such as when the system buses are beyond or near the point o f voltage collapse. With the proposed novel formulas, the computational cost of generating P V curves across all busses o f the power system network is greatly reduced, and the resulting curves clearly show the point o f voltage collapse as well as any unstable voltage/power combination thereafter. The ability to quickly obtain P V curves makes the proposed techniques particularly well suited for online power system monitoring o f operating 11 states. A l s o , since these new difference equations introduce very little error for frequencies lower than the tuned frequency, they are also well suited for the solution o f dynamic power flow problems such as those resulting from loads variation, transformers tap change, and outages. 111 Table of Contents Abstract ii Table o f Contents iv List o f Tables vii List o f Figures viii Acknowledgments xiii C H A P T E R 1: Introduction and Overview 1 1.1 Background and Motivation .- 1.2 Electric Power Circuit Analysis Techniques 1 4 1.2.1 Nodal Analysis 1.2.2 Modified Nodal Analysis 16 1.2.3 Phasor Analysis 18 1.2.4 Power flow Analysis 19 1.3 5 Outline o f the Thesis C H A P T E R 2: 22 Adjustable Methods 23 2.1 Introduction 23 2.2 Derivation o f the Adjustable Methods 25 2.3 Example Method Derivation 26 2.4 Magnitude and Phase Distortion 28 2.5 Stability 30 C H A P T E R 3: 3.1 Tuned Time-Domain Linear A C Steady State Solutions Introduction 34 34 iv 3.2 Tuning the Adj ustable Methods 36 3.3 Circuit Simulation Example 42 3.4 Reducing the Number o f Circuit Nodes 43 3.5 Estimating Magnitude and Phase 46 3.6 Tuned Time-Domain Steady State Simulations 47 3.7 Computational Cost o f the Tuned Time-Domain Methods 53 C H A P T E R 4: Solving the Static Power flow Problem 56 4.1 Introduction 56 4.2 Modeling the Power System Elements 57 4.2.1 Modeling Transmission Lines 57 4.2.2 Modeling Transformers 58 4.2.3 Modeling Shunts 59 4.2.4 Modeling P Q Loads 59 4.2.5 Modeling P V Sources 63 4.3 Tuned Time-Domain Static Power flow Test Cases 69 4.3.1 I E E E 14-Bus test system 69 4.3.2 C I G R E 32-Bus test system 72 4.3.3 I E E E 57-Bus test system 74 4.4 Computational Cost o f the Tuned Time-Domain Power flow 77 4.5 Generation o f P V Curves 78 4.5.1 Effect o f the Current Increment Size 81 4.5.2 Comparison with Published P V Curves 82 4.5.3 Computing the P V Curves for all System P Q Buses 84 v 4.5.4 A Note on Load Voltage Stability Indices C H A P T E R 5: Dynamic Power flow Using Tuned Time-Domain Formulas 86 88 5.1 Introduction 88 5.2 Dynamic Modelling of U L T C Transformers 90 5.3 Dynamic Modelling of Loads 95 5.4 Simulating Outages 99 5.5 Example of Dynamic Power flow Solution C H A P T E R 6: Conclusions and Recommendations for Future Work 101 106 6.1 Summary and Conclusions 106 6.2 Recommendations for Future Work 109 Bibliography Ill vi List of Tables Table 3.1 Tuned values o f k and co for the difference equations obtained with polynomials o f orders 4 to 9 37 Table 3.2 The tuned methods obtained with interpolating polynomials o f order 4 to 9 are A(a)-stable. For the method obtained with a seventh order interpolating polynomial is almost A-stable 40 Table 3.3 Steady state voltages o f the circuit i n Figure 3.9. The D C values are obtained when a 10 V dc voltage source is connected i n series with V I 50 Table 3.4 Steady state solution times for phasor analysis and a seventh order tuned timedomain formula 55 Table 4.1 Node voltages for the system o f Figure 4.13, solved using Newton-Raphson iterations and tuned time-domain power flow 72 Table 4.2 Static power flow voltages solution for the circuit o f Figure 4.14 and Figure 4.15. 74 Table 4.3 Power flow solution times using the tuned time-domain and Newton-Raphson methods 78 Table 5.1 Some typical coefficients for the exponential load model o f equation (5.2) 95 Table 5.2 Step values over time for the load connected at bus 5 i n the system o f Figure 4.11. .- 96 vn List of Figures Figure 1.1 Simple circuit for Nodal Analysis example 6 Figure 1.2 Equivalent circuit o f an inductor when a difference equation is used to model it. 8 Figure 1.3 Circuit with a capacitor and an inductor (top), and its equivalent discrete representation (below) 9 Figure 1.4 Magnitude distortion and phase error introduced by the backward Euler formula. 10 Figure 1.5 L o c i o f the roots for the B D F formulas o f order 1 to 6 13 Figure 1.6 Circuit used in the simulations to test the stability o f different numerical methods. 14 Figure 1.7 For the circuit o f Figure 1.6 the product hA (the two dots, h=\ms) is inside the area in which the difference equation is unstable '. 14 Figure 1.8 Voltage across the capacitor o f Figure 1.6. After a number o f time steps the simulation becomes unstable since the product hA is not in the area o f stability (Figure 1.7) 15 Figure 1.9 One way o f obtaining a stable response is to decrease (or increase) h so the product hA is in the area o f stability 15 Figure 1.10 The voltage i n the capacitor o f Figure 1.6 after decreasing the value o f h to 0.5 ms is stable, since the product hA is in the area o f stability as shown in Figure 1.9 16 Figure 1.11 Simple circuit for Modified Nodal Analysis ( M N A ) example 17 Figure 1.12 Simple circuit for Phasor Analysis example 19 viii Figure 2.1 For a number o f discrete points (diamonds) a fitting polynomial can be derived (dashed line). The derivative at the last discrete point can be approximated using different points, resulting in slopes o f lines SO (k=0.0) to S3 (k=1.0) 24 Figure 2.2 Magnitude distortion and phase error introduced by equation (2.17) for the indicated values o f k 29 Figure 2.3 Stability zone for equation (2.17) for the indicated values o f k 30 Figure 2.4 Circuit used to test equation (2.21) 31 Figure 2.5 Magnitude distortion and phase error introduced by equation (2.21). The dashed line corresponds to the trapezoidal rule, included here for comparison purposes 32 Figure 2.6 Voltage in phase ' a ' o f B u s _ l when the trapezoidal rule is used to simulate the circuit o f Figure 2.4. Notice the numerical oscillations when C B 1 is tripped at 60 ms. 33 Figure 2.7 Voltage in phase ' a ' o f B u s _ l when equation (2.21) is used to simulate the circuit of Figure 2.4 33 Figure 3.1 Magnitude distortion and phase error introduced by the adjustable method o f order 5 for &=0.0 to 0.4 i n 0.1 increments 35 Figure 3.2 Magnitude distortion and phase error o f tuned difference equations derived with polynomials o f orders 4 to 9 and the values o f k in Table 3.1 Figure 3.3 L o c i o f the roots for the tuned methods o f orders four to nine 38 39 Figure 3.4 This Maple 6 script can be use to computed the coefficients o f the adjustable difference equations, as well as to find their tuned values for interpolating polynomials of order > 4 41 Figure 3.5 Circuit used to test equation (3.9) 42 ix Figure 3.6 The solid line is the voltage in the inductor o f Figure 3.5 obtained using phasor analysis. The dots are the values obtained when simulating the same circuit using the tuned difference equation (3.9) Figure 3.7 Discrete equivalents for R L and R C circuits 43 45 Figure 3.8 Program flow chart to compute the steady state with the tuned formulas. The simulation is stopped once the steady state is found, or a maximum number o f time steps is reached 48 Figure 3.9 Circuit from which the steady state is computed using equation (3.10). Inductances are in m H and capacitances are in [iF 50 Figure 3.10 Magnitude and phase change over time for the voltage at node ' D ' o f the circuit of Figure 3.9 for A C excitation using the tuned formula o f order 7 51 Figure 3.11 Magnitude and phase change over time for the voltage at node ' D ' o f the circuit of«Figure 3.9 using the trapezoidal rule Figure 4.1 Equivalent jt circuit o f transmission line 52 58 Figure 4.2 Typical two winding transformer model used in power flow studies (top) and its equivalent circuit (bottom). The values o f Y l , Y 2 , and Y 3 are given by equation (4.1). 59 Figure 4.3 Damped numerical oscillations caused by a non-tuned formula derived from a seventh order interpolating polynomial. Similar behaviour is observed with the tuned formulas 61 Figure 4.4 Circuits used to test the P Q load model. The circuit on the right yields a solution, whereas the circuit in the left is i l l conditioned and the voltage at node ' B ' collapses.. 61 Figure 4.5 Magnitude and phase changes over time for the circuits o f Figure 4.4 62 Figure 4.6 PI controller used to adjust the phase angle o f the P V sources 65 Figure 4.7 Circuit used to test the P V sources model 65 Figure 4.8 Voltages magnitude and phase for the circuit o f Figure 4.7 when Kp=0.11 and Ki=0.011 66 Figure 4.9 Voltages magnitude and phase for the circuit o f Figure 4.7 when K p = l .0 and Ki=0.1 67 Figure 4.10 Voltages magnitude and phase for the circuit o f Figure 4.7 when Kp=0.01 and Ki=0.001 68 Figure 4.11 I E E E 14-bus test system 70 Figure 4.12 Voltages magnitude and phase change over time for the I E E E 14-bus test system solved using tuned time-domain power flow 71 Figure 4.13 32-bus test system. The resulting voltages magnitudes and phase for this system are presented in Table 4.1 73 Figure 4.14 I E E E 57-bus test system, sheet 1 o f 2 75 Figure 4.15 I E E E 57-bus test system, sheet 2 o f 2 76 Figure 4.16 Load P V curve at bus A for the circuit o f Figure 4.4. The solid line corresponds to the values obtained from the tuned time-domain simulation. The dashed line corresponds to a fifth order interpolating polynomial obtained from the solid line. The diamonds are the solutions obtained using Newton-Raphson iterations 80 Figure 4.17 Effect o f the Al size in the generated P V curves. A small value o f Al results in a more accurate and detailed P V curve, but at a higher computational cost. The plot shows the curves generated with a Al o f 0.5, 0.05, and 0.005 xi 81 Figure 4.18 Load P V curve for bus N207 o f the system o f Figure 4.13. The dashed line is the curve reported in [37]. The solid line was obtained using the tuned time-domain power flow described in this report 83 Figure 4.19 A l l load P V curves for the 14-bus system o f Figure 4.11 85 Figure 4.20 The load P V curve for bus 38 o f the I E E E 57-bus system showing the point at which the voltage at bus 31 collapses first (diamond) 87 Figure 5.1 Basic connexion o f a U L T C transformer 91 Figure 5.2 Effect o f the U L T C transformers on the static power flow solution o f the system of Figure 4.11 93 Figure 5.3 Effect o f the U L T C transformers in the P V curves o f buses 4 and 5 for the system of Figure 4.11 94 Figure 5.4 Magnitude and phase for the voltages at buses 4 and 5 for the system o f Figure 4.11 when the load at bus 5 changes with the values given in Table 5.2 97 Figure 5.5 Magnitude and phase for the voltages at buses 4 and 5 for the system o f Figure 4.11 when the load at bus 5 changes linearly with time 98 Figure 5.6 Magnitude and phase for the voltages at buses 4 and 5 for the system of Figure 4.11 when the line connecting both buses is removed 100 Figure 5.7 Dynamic voltage power flow solution at buses N201, N202, and N204 for the O G R E 32-bus system o f Figure 4.13 using pure P Q loads 104 Figure 5.8 Dynamic voltage power flow solution at buses N201, N202, and N204 for the C I G R E 32-bus system o f Figure 4.13 using voltage dependant loads xii 105 Acknowledgments I owe a tremendous debt o f gratitude to Dr. Jose R. Marti, my Ph. D . thesis supervisor. His constant guidance, advice, help, encouragement, and patience as well as financial support, have been invaluable during my Ph. D . research. I want to thank my friends Dr. Roberto Rosales and Dr. Luis Linares for the help, advice, and encouragement they gave me in all these years at U B C . A l s o , I want to thank Dr. Ray Burge and Peter Dickson for helping me when I needed it the most. Finally I want to thank my senseis, Dr. George Iwama, Dr. Darin Bennet, and Claudio Lerner not only for all their teachings but also for their example o f excellence and discipline. Xlll CHAPTER 1: Introduction and Overview 1.1 Background and Motivation Numerical methods for the solution o f systems o f differential equations are at the heart o f most transient analysis electric circuit simulators. Popular computer programs such as the E M T P , S P I C E , as well as many o f their predecessors and derivatives use such methods. For some, such as the E M T P , only the well known Backward Euler and Trapezoidal one step methods are used [1;2]. Multistep methods have been also used to simulate electric circuits. For instance, the well-known Adams-Bashforth and Adams-Moulton methods are sometimes used, although they lack much o f the stability required for a general-purpose electric circuit simulator [3]. It was in 1952 when the multistep Backward Differentiation Formulae ( B D F ) methods were first described [4], and in 1963 when Dahlquist [5] described the properties o f such multistep methods. Dahlquist demonstrated that such methods could be derived using Taylor series approximations from which a truncation error is computed. Moreover, he studied the methods stability properties; in particular, a very desirable property called A stability. This property is defined as the ability o f a given method to converge to a stable solution when it is used to solve a system o f first order differential equations, where the real part o f the eigenvalues o f the system are all negative. Since most passive electric circuit have negative real part eigenvalues, A-stability is a very desirable property to have in the numerical method used to simulate such circuits. 1 B D F methods came to prominence in 1971 when Gear [6] used them in his generalpurpose O D E solver . A t the time it was known that such methods were A-stable for orders 1 one and two, and A(ct)-stable for orders three to six [7]. The A(ct)-stability property has to do with the ability o f a given multistep method to converge to a stable solution not for all real part negative eigenvalues but for just a range which often comprises eigenvalues far apart from each other. This seemed to be an appropriate solution since most electric circuits show this eigenvalues range. Systems with this characteristic are often referred as 'stiff systems' [8]. They comprise not only electric circuits, but many other mechanical and chemical systems as well [9]. B D F methods are the most widely used methods for the solution o f stiff systems o f differential equations [10]. Unfortunately this class o f methods become unstable for orders larger than six [11; 12]. Many techniques have been suggested to increase the stability and order o f B D F methods: [12;13;14;15;16;17] just to name a few. Although very effective, in general, these approaches require very small or adjustable time steps that impose an important burden in the performance o f a general-purpose circuit simulator. Simpler numerical methods often result in faster simulators. One such method that stands over the rest is the trapezoidal method, whose properties are well described by Dahlquist in [5]. The trapezoidal method has also drawbacks, in particular, the oscillatory behaviour o f this method when discontinuities are came across is well known [2; 18; 19]. One of the first solutions to this problem consisted o f the addition o f external damping [20], which is functionally equivalent to a class o f adjustability or tuning presented by many other authors[3;21]. Since the addition o f damping results i n an effective modification o f the Reference [6] includes the FORTRAN II source code for the general-purpose ODE solver DIFSUB that became very popular in the 1970's. That is the reason why BDF methods are also known as Gear's Formulas. 1 2 problem being solved, techniques that are more convenient were eventually developed. Gear and Osterby in [22] explored the alternative o f order reduction to get rid o f undesirable numerical oscillations. It was not until the Critical Damping Algorithm ( C D A ) was introduced by Marti and L i n [2] that an efficient, simple, and elegant method to suppress the numerical oscillations o f the trapezoidal rule was obtained. Nevertheless, the problem remains when the solution o f the system has to be obtained in real-time, and the time step reduction required by C D A is not possible [23]. Another important problem often found in electric circuit simulations is the so-called oscillatory problem. In many systems o f ordinary differential equations, such as those resulting from electric circuits, complex conjugate eigenvalues or forcing periodical functions result in oscillatory responses. The use o f phasor analysis, described in many basic circuit analysis textbooks, is a well-established technique for the solution o f the linear oscillatory problem in electric circuits. Unfortunately, the presence o f non-linear elements in the system being simulated renders this technique unusable and alternative methods, such as Newton-Raphson iterations are used [24],[25]. Even for the simpler non-linear power flow problem involving linear network elements but non-linear loads, phasor analysis is not appropriate, and again iterative techniques such as Newton-Raphson iterations, or one o f its variations, are used [26]. This thesis describes a family o f new methods used for the simulation o f electric circuits. These methods have very desirable properties such as A-stability, high accuracy, and low computational burden. In particular, two kinds o f methods are described. The first one, derived from a third-order interpolating polynomial, shows stability properties similar to that o f the Backward Euler method, and the accuracy o f the Trapezoidal method. 3 This method is useful for the real time simulation o f circuits with discontinuities, but it is mainly presented in this report as a simple example o f the general method's derivation. The other group o f methods are derived from interpolating polynomials o f fourth or greater order, tuned at a particular frequency. These are fixed time-step methods, which produce an exact solution for linear and non-linear A C circuits, both in magnitude and phase for the frequency they are tuned to. Since almost all electric power systems work at one fixed frequency, the potential use o f these methods for the solution o f such circuits is evident. In this thesis, tuned time-domain methods are used to accurately solve linear steadystate circuits as well as non-linear static and dynamic power flow electric power circuits. For the simulation o f linear steady-state circuits, the methods are faster than the well-known phasor analysis for circuits with a large number o f nodes. When used to solve the static power flow problem, the cost o f the tuned time-domain methods is similar to that o f NewtonRaphson iterations. But unlike Newton-Raphson iterations, the tuned time-domain methods produce useful results for ill-conditioned systems, such as those undergoing voltage collapse. The ability o f the tuned time-domain methods to solve electric power systems at or beyond the point of voltage collapse makes them ideal for the economical generation o f P V curves. A l s o , since the methods are still accurate enough for frequencies much lower than the tuned frequency, they can also be used to simulate the dynamics o f power flow electric circuits. 1.2 Electric Power Circuit Analysis Techniques In this section, four electric power circuit analysis techniques, relevant to this research are briefly described. The first two, nodal analysis and modified nodal analysis are used to implement the algorithms and techniques described in this research. The other two, 4 phasor analysis and power flow analysis are the techniques against which this research is validated and compared. 1.2.1 Nodal Analysis O f the available techniques for the computer simulation o f electric circuits, nodal analysis is the easiest to implement [27]. A variation o f nodal analysis, the so called modified nodal analysis or M N A [28], was used on this research to implement a circuit simulator to develop and test the new simulation techniques. Nodal analysis is based on K i r c h h o f f s current law, which states that the algebraic sum o f currents entering a node is zero. Nodal analysis consist on finding the admittance matrix Y and the vector o f independent current sources / so that YV = I, (1.1) where the dimension o f matrix Y is equal to the number o f nodes in circuit minus one (for the reference node), and solving for the voltages V. Matrix Y and vector / can be built by inspection following these three simple rules: 1) The diagonal elements o f the Y matrix come from the sum o f all admittances at the corresponding node. 2) The non-diagonal elements o f the Y matrix come from the negative o f the admittance connecting the corresponding nodes. 3) The elements o f the independent current sources / results from the sum o f all current sources coming into the corresponding node. 5 1 2 g3 g2 gi I2 Figure 1.1 Simple circuit for Nodal Analysis example. For instance, consider the circuit o f Figure 1.1. Notice first that it has two nodes (excluding the ground reference node), therefore the Y matrix w i l l have two rows and two The sum o f admittances at node 1 is gl+g3. columns. g2+g3. Similarly, for node 2 the sum is The negative o f the admittance connecting nodes 1 and 2 is -g3. The sum o f all current sources entering node 1 is II, and the sum o f all current sources entering node 2 is 12. Therefore -g3 g2 + g3 -g3 VI 11 V2 12 (1.2) The system o f equations (1.2) can be generalized to any number n o f nodes as Y Y 1 1 ll Y 21 1 Y„i ln l 12 Y 22 l Y„2 ••• ••• 'V, Y 2n 1 (1.3) Y m or for one particular node i as (1.4) Nodal analysis can be applied to dynamic elements such as capacitors and inductors as well. This is achieved by discretizing the continuous time element using a difference 6 equation. The application o f a difference equation in the simulation o f electric circuits is based in the substitution o f the element being simulated by an equivalent one which w i l l perform the same task as the differential equation that represents such element [1]. Let us, for example, perform this substitution for the voltage across an inductor, given by V (t) = L L^P-. (1.5) at The idea is to replace the differential part o f the equation by a difference formula. Let us choose, for simplicity, the Backward Euler formula di (t) i(t)-i(t-h) L dt (1.6) h where h is a time step, small enough so the approximation does not introduce too much error. Replacing (1.6) back into (1.5) (1.7) h and getting i(t) out results in i(t)«jV (t) + i(t-h) L . (1.8) Equation (1.8) is usually represented by the equivalent circuit shown in Figure 1.2, in which the discrete equivalent o f the inductor is a resistor with a current source in parallel, whose value depends on the previous value o f the current on the inductor: h = i(t - h) L _h 7 (1.9) •-o Figure 1.2 Equivalent circuit o f an inductor when a difference equation is used to model it. The procedure described above is also valid for capacitors, resulting i n a similar equivalent circuit. Figure 1.3 shows a circuit with a capacitor and an inductor as well as its equivalent discretized representation. Applying nodal analysis to the discretized circuit results in the system o f equations gl + gL -gL -gl gL + gC I\-hL [V2_ hL + hC + I2 (1.10) where hL and hC change from time step to time step accordingly to discretization formula used. 8 Figure 1.3 Circuit with a capacitor and an inductor (top), and its equivalent discrete representation (below). A s mentioned before, the circuit o f Figure 1.2 is just an approximation. H o w good this approximation is can be determined either by finding the transfer function o f the equivalent circuit and comparing it with that o f the real circuit as proposed i n [2], or plotting the frequency response o f the Z-transform o f the difference equation [29] . 2 Then, for an inductor o f one Henry, discretized using the backward Euler formula, it can be shown [2] that 1 X{Z): z (1.11) hz-l is its Z-transform. N o w , replacing z with e^ , the magnitude and phase of X(z) can be plotted wh as shown i n Figure 1.4. These plots represent both the magnitude and phase errors for the discretized inductor for the range o f frequencies D C to l/(2h) . Page 494, equation 12-83, but it is easier to followfromequations 12-86 and 12-87 in page 495. 9 2.0-r1.8-1.6 4- cc ^ 0.8-0.6 - - \ | j < f j | ; \ 0.4-- j \ i \ \ j { j \ 0.2-- \ | j \ I ] | I \ 0.0-I 0.00 1 1 1 1 1 j 1 1 1 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 Frequency in per unit of (1/h) -30--45--60--75- - 0.00 Frequency in per unit of (1/h) gure 1.4 Magnitude distortion and phase error introduced by the backward Euler formula. 10 To analyze the stability properties o f a given rule many authors ([6], [29], [27], [3]) refer to the initial-value problem v'(t) = kv(t), v(0) = l , (1.12) v(0 = e " , 0-13) which has the known solution where k is the time constant or eigenvalue o f equation (1.12). This time constant is, i n general, a complex number whose real part is negative: Re(l)<0. For larger systems, it is possible to apply a modal transformation i n such a way that a set o f independent equations with the same elementary form can be obtained. In that case, k w i l l represent the eigenvalues of the state matrix o f the system and v(t) w i l l be a vector. For example, consider the following formula, which is the B D F second order formula, and analyze its stability properties: 3v(t)- 4v(t-h) + v(t-2h)~ h(2v '(t)). (1.14) If we use the initial value problem o f equation (1.12) and replace it i n equation (1.14), we get 3v(t) - 4v(t - h ) + v(t-2h) = h(2Av(t)). (1.15) Following the procedure outlined in [27] we want the product hk h A = 3 (t)-4 (t-h V V ) + v( -2h)_ t 2v(t) N o w , lets evaluate the function using equation (1.13) for all possible eigenvalues jco, and see what values o f hk lay inside it: hM«>)= 3 2 Backward differentiation formulas. 11 e M • (1.17) Simplifying hXM= 3 - 4 -J E + W 2 + P- ^ 2 e . (1.18) Regrouping we get the loci o f the roots o f equation (1.14) 3e W<») = 2jUi -4e ju> Z +1 2g • (1.19) If we evaluate equation (1.19) for 0<ax2jt, and plot the result in the complex plane, we get the locus o f hX (Figure 1.5). N o w , by plotting the actual values hX o f the simulated system in the complex plane (the eigenvalues multiplied by the time step), it is possible to determine i f the rule is stable or not for such system. For some formulas the area o f stability is inside the closed curve, for example, all Adams-Bashforth formulas, and all AdamsMoulton formulas o f order larger than two ; for others, such as Backward Euler, Trapezoidal , 4 5 and all B D F formulas, the area o f stability is outside the curve [29]. If the area o f stability lays inside the closed curve the method is explicit, otherwise the method is implicit. Observe from Figure 1.5 that the first and second order B D F 6 are stable for all eigenvalues X with negative real part, despite the value o f the time step h. This property, once more, is known as A-stability [5]. A-stability is important because for most practical electric circuits the real part o f the eigenvalues is negative. To demonstrate the importance o f A-stability, consider the circuit o f Figure 1.6, whose eigenvalues are the roots o f 4 The second order Adams-Moulton is also the Trapezoidal rule. The stability zone of the Trapezoidal rule is the complete complex left plane. It can be pictured as a circle of infinite ratio centered in the right hand side of the complex plane touching the origin. The first order B D F is also the Backward Euler formula. 5 6 12 -8 —4 0 4 8 12 Re(z) 16 20 24 28 32 Figure 1.5 L o c i o f the roots for the B D F formulas o f order 1 to 6. (Lfifi; )S 2 +(L + C^Rj )S + (R + R ) = 0. 1 2 (1.20) First, the circuit o f Figure 1.6 is simulated using the fourth order B D F formula and A=1.0 ms. The area o f stability for the B D F fourth order formula, and the values hX for the given circuit are shown in Figure 1.7, which in this case lie inside the curve, that is, in the area in which the rule is not stable. The effect on the time-domain simulation is shown in Figure 1.8 where the voltage across the capacitor is plotted. noticeable unstable. Notice that after about 0.04 s the output is If the value o f h is set to 0.5 ms, then the product hX lies outside the curve (Figure 1.9), and the output o f the time-domain simulation is numerically stable (Figure 1.10). 13 Figure 1.6 Circuit used i n the simulations to test the stability o f different numerical methods. Figure 1.7 For the circuit o f Figure 1.6 the product hA (the two dots, h=lms) is inside the area in which the difference equation is unstable. 14 Figure 1.8 Voltage across the capacitor o f Figure 1.6. After a number o f time steps the simulation becomes unstable since the product hA is not in the area o f stability (Figure 1.7). 12.0-1 •n -12.0 u . 1 ! : 1 : : : , : i i i i I I I i I I -9.6 -7.2 -4.8 -2.4 0.0 Re(z) 2.4 4.8 7.2 9.6 12.0 Figure 1.9 One way o f obtaining a stable response is to decrease (or increase) h so the product hA is in the area o f stability. 15 6.01 4.83.62.4> 1 - " 2 S"o.oCO 1 -1.2-2.4-3.6-4.8- 1 1 I — i — i — i — -6.00.0000 0.0099 0.0199 0.0298 0.0398 0.0498 0.0597 0.0697 0.0796 0.0895 0.0995 Time (s) Figure 1.10 The voltage in the capacitor o f Figure 1.6 after decreasing the value o f h to 0.5 ms is stable, since the product hA. is in the area o f stability as shown in Figure 1.9. 1.2.2 Modified Nodal Analysis Standard nodal analysis, as described i n the previous section, can only deal with independent current sources. Although it is possible to add independent voltage sources to the nodal analysis formulation by means o f Norton equivalent circuits, modified nodal analysis or M N A [28] is often simpler and easier to implement. Additionally, it solves for the current coming out o f the voltage sources, a handy figure to have when power computations are needed. M N A basically consists on the aggregation o f extra equations and variables to the standard nodal analysis formulations so to impose new circuital conditions, as for example, an applied voltage between two nodes. Consider for instance the circuit o f Figure 1.11, which is the same circuit o f Figure 1.1 with a voltage source o f value E between nodes 1 and 2. Assuming the current from this additional voltage source moves from the positive terminal to the negative terminal, it w i l l 16 add to the total current o f node 2 and subtract to the total current o f node 1 as well as to set the voltage difference between nodes 1 and 2 to a given value: V1-V2 =E (1.21) where IE+ and Is- are the currents flowing into the positive and negative terminals o f the voltage source. The modified system o f equations is therefore gl + g3 -g3 g -g3 1 VI 2 + g3 -1 V2 -1 0 h 1 11- - 12 (1.22) E Comparing the system o f equations (1.22) and (1.2) we observe that: a) The matrix dimension has been incremented by one, despite the number o f nodes remaining constant. b) The current through the voltage source is one o f the unknowns, along with the node voltages. Not only independent voltage sources, but also dependent current and voltage sources, ideal transformers, and operational amplifiers can be represented using M N A [27]. But since they are not relevant for this report the reader is referred to the aforementioned reference. E (+ ') 93 g2 gi Figure 1.11 Simple circuit for Modified Nodal Analysis ( M N A ) example. 17 1.2.3 Phasor Analysis In many electrical circuits, such as those arising in electric power and communications, the most frequent type o f source excitation is sinusoidal. Phasor analysis is similar to nodal analysis, but using complex numbers instead. This is possible because phasor representations for sources, inductors, and capacitors can be easily obtained. For instance, consider the top circuit o f Figure 1.12, for which the following transformations are performed I\*cos(wt + Ph\)^I\\Ph\ 12 * cos(wt + Phi) -> I\\Ph2 (1.23) 1 COL C->ytoC,y resulting in the bottom circuit o f Figure 1.12. The voltages o f the resulting circuit can be solved using nodal analysis: J coZ, J (X)L f c o C - _L A V2_ ~ n\Ph\~ n\Phi (1.24) coZ, Similarly to the admittance matrix and the vector o f independent sources, the vector o f voltages resulting from the solution o f (1.24) is made up o f complex numbers as well. M N A techniques can be applied to phasor analysis too, as described in the preceding section, allowing for the easy incorporation o f independent voltage sources. 18 -j/wL II LRhl I2l£h2 gi jwc r • . — * — • T - Figure 1.12 Simple circuit for Phasor Analysis example. 1.2.4 Power flow Analysis "Power flow" or "load flow" analysis as it is some times called, is a fundamental tool for power system operation and planning. According to Bergen & Vittal [26], power flow is the most common o f the electric power system computer simulations performed. It is somehow similar to phasor analysis in the sense that the magnitude and phase o f the system buses voltages have to be computed, but instead o f using currents as the independent variables, it uses powers. To describe how power flow analysis is performed, let us start with the definition o f complex power S,=VX, (1.25) The complex conjugate of the current is used because of the way the complex power S is defined. The real part of the complex power is called active power or P and is given by P=\V\\I\cos(9), where 9 is the angle between the voltage and the current; the number resulting from cos(9) is known as the power factor. The imaginary part of S is called reactive power or Q and is given by Q=\ V\ \I\sin(9). If V=\and I=\I\ef ^ the complex power is S=Vl'=\ V\ \I\e/' e^ =\V\\I\(cos(9)+) sin(9))=P+jQ. f e) 19 where Sj, V , and t are phasor quantities. A s i n the nodal analysis formulation equation (1.4), the net current into a bus can be computed as ( V n I = =z^; J U=i (i-26) 4=1 Replacing equation (1.26) into equation (1.25) results i n S^V^YX (1.27) k=\ To simplify the computations, the following substitutions are performed to equation (1.27) V, = \KW A di k = ^-e ,Y =G +jB , k lk lk (1.28) lk resulting i n S^fyWV^^-jB,,). (1.29) A p p l y i n g Euler's substitution to the exponential part, results i n S, =P,+jQi =Z|^||^|(cos(9 ) + 7sin(e ))(G -y5 ). l t f t l k l t (1.30) Finally, separating real and imaginary parts yields the equations to be solved for each o f the system buses > =Z l < I \ * I P F v c o s ^ - + * ( < - *» • B sin B 6 a = ZI^IKI(GftSin(e -e )-^cos(e -e )) < t J t 0 - !) 3 (1.32) In typical power flow problems the active and reactive power on some buses is specified, therefore those are known as P Q buses, for which both equation (1.31) and equation (1.32) must be used to compute the bus voltage magnitude and phase. For other buses, the magnitude o f the voltage and active power are specified, the so called P V buses, 20 for which only the active power equation (1.31) is needed to compute the bus voltage angle. Since the total flow o f power is not known a priori, it is customary to set one system bus as free to provide any P or Q. For this swing or slack bus, both the voltage magnitude and phase are assumed, and neither equation (1.31) nor equation (1.32) are required. Notice from equations (1.31) and (1.32) that the resulting system o f equations is nonlinear due to the products o f voltages. The preferred method for the solution o f this nonlinear system o f equations is Newton-Raphson iterations. It consists o f the linearization of the system to be solved using the first term o f the Taylor's series expansion. Before proceeding with such linearization, equations (1.31) and (1.32) can be rewritten as f (v ,...,v ,e ,...,Qj = f ^ Pi l n l j l , k= / (F a 1 > ..:,K ,e ,...,ej = B 1 (1-33) 2l^||^|(G sin(G -e )-5 cos(e -et))-a a t / tt < k=\ which after linearization with the first term o f the Taylor's series becomes, ¥PI dfpi ae, - ae„ d\v\ dfpi dfpi s\v m A6 7 fp d dfpn n ¥ P „ PI df Pn ae, "" ae„ d\v.\ "' a i r fp n fi d ¥QI_ Q ¥QI_ ¥QI_ ae, "' ae„ d\v\ "' d\vm JQ, - (1.34) /QI A IK. I fcQm ae„ a ^ | - d\vm where n is the sum o f the number o f PQ and PV buses for which the phase angle 9 o f the voltages need to be computed; and m is the number o f PQ buses, for which only the voltages magnitude need to be computed. This results in a Jacobian matrix o f dimension (n+m). The 21 system o f equations (1.34) is solved by assuming an initial value for the buses voltage magnitude and phase with which to compute the values o f the functions in equation (1.33) and the Jacobian matrix o f the system o f equations (1.34). Then the changes on the buses voltage magnitude and phase are computed from (1.34), and the process is repeated until the changes in the buses voltage magnitude and phase reach a given tolerance. 1.3 Outline ofthe Thesis In Chapter 1 we reviewed some o f the required background information relevant for this thesis. In Chapter 2 the derivation o f adjustable difference equations is presented as well as their application to time-domain simulations. The tuning o f the adjustable methods is presented in Chapter 3, with their application i n the solution o f forced sinusoidal circuits in the time-domain and comparison to phasor analysis. Chapter 4 deals with the application o f the tuned time-domain methodology described in Chapter 3 to the solution o f the power flow problem and its comparison to Newton-Raphson iterations. A l s o in Chapter 4, some unique properties o f the tuned time-domain methodology are used to assess the voltage stability o f electric power systems. simulated. In Chapter 5, the dynamics o f the power flow are modelled and Finally, Chapter 6 presents the conclusions o f this work, as well as some recommendations for future work. 22 CHAPTER 2: Adjustable Methods 2.1 Introduction Linear multistep methods are usually denoted with the following difference equation[3] j=0 j=0 ^ -> ' A l y'^myhtZO.f^fft.y,) are the method's coefficients, h is a constant time step, and q is usually the where Oj and order o f the method. Equation (2.1) can represent, for example, the Trapezoidal, Backward Euler, Adams-Moulton, Adams-Bashforth, and B D F methods by selecting adequate values for Oj and fy. In fact, the values o f ccj and can be chosen so that equation (2.1) meets any desired criteria. Generally the criteria to choose the method's coefficients CCJ and intends o to maximize the order o f the method p so that C = C =... = Cp=0,Cp+1*0 0 (2.2) 1 where Co to C +i are given by p j=0 (2.3) J- -if<*j — —ij'~% + i!fj J J ! (i-l)'.U i ,i = l,2,... j In this chapter a technique for choosing the values o f Oj and is described, so that the coefficients render a method with some stability and accuracy properties that are useful for the simulation o f electric circuits. Before going to the mathematical derivation o f the adjustable methods, a graphical explanation o f the technique can be useful. Reference [3], starting at page 126 has a more complete and detailed explanation. 23 Consider, for example, the six discrete points represented by the diamonds in Figure 2.1. A polynomial o f order five can be derived from those six points (using for example the Matlab function polyfit), from which the dashed curve o f Figure 2.1 is obtained. N o w , let us estimate the derivative at the last discrete point. If the previous discrete point is used, the derivative results in the slope represented by line SO. This is equivalent to the well known Backward Euler method. If the derivative o f the fifth order polynomial represented by the dashed line is used, the results is the slope represented by line S3, which is the same as the fifth order B D F . Intermediate points from the derived polynomial can be used as well, as depicted by the slopes of lines SI and S2. In this chapter, the formulation for methods resulting i n slopes similar to those o f S1 and S2 is derived, and their stability and accuracy properties described. 53 § ' 1 \ 2 \ \ : -V \ \ \ 1 ...» :.. * \ ; •. / \ / / .. ; K : / ^ \ - \ \ \ \ ^ - \ \ sa _ ^ - o \ ~~~~;~*f- . . . A V 1 1 1 1 1 2 3 4 ! Figure 2.1 For a number o f discrete points (diamonds) a fitting polynomial can be derived (dashed line). The derivative at the last discrete point can be approximated using different points, resulting in slopes o f lines SO (k=0.0) to S3 (k=1.0). 24 2.2 Derivation of the Adjustable Given a number o f Methods n+1 values v(t), v(t-h), v(t-nh), sampled at points t, (t-h), (t-nh) we can assume there is a polynomial o f degree n v(t)=f a t m u m (2.4) m=0 which produces the exact values at the selected sampled points. The coefficients o f such polynomial can be found by solving the system o f equations m=0 vO-hJ^aJt-h)" m=0 (2.5) v(t-nh)=^ajt-nh)" m=0 The solution o f (2.5) produces coefficients that are a function o f the sampled values v(t), v(t-h), v(t-nh), and the time step h such that o = f(y(t), v(t-h),...,v(t-nh),h) a =f(y(t),v(t-h),...,v(t-nh),h) l (2.6) „ = f( (0, (t -h),...,v{t- nh), h), a v v from which the interpolating polynomial o f equation (2.4) can be completed. If the derivative of equation (2.4) is taken, the well-known backward difference formula o f order equal to the order o f the interpolating polynomial is obtained: ^2>o.<" m=l 25 (2.7) On the other hand, from the basic definition o f the derivative, dv(t) dt = v'(0 = lim v(t)-v{t-h) (2.8) h and replacing a variable h by kh where h is fixed and H s a real variable, v'(0 = lim v(t)-v(t-kh) (2.9) kh The delay term v(t-kh) in (2.9) can be represented by the interpolating polynomial o f equation (2.4). Furthermore, it is possible, by selecting small, positive, and different from zero values for k and h to get a close approximate value o f v '(t) v(t)-v(t-kh) v(0 = (2.10) kh B y setting k to one, the backward Euler formula is obtained. A l s o , by setting k to zero the B D F formula whose degree corresponds to the used interpolating polynomial is obtained. Therefore, this technique permits the creation o f a customized formula for the solution o f differential equations which accuracy, stability, and frequency response can be adjusted by choosing a value o f k between zero and one. 2.3 Example Method Derivation The above-described methodology can be better illustrated with an example. Consider the third order polynomial v(t) =a t + a t + a t+a 3 (2.11) 2 3 2 1 to approximate the value o f v'(t) as a function o f k. i To find the coefficients of this polynomial four values are required, v(t), v(t-h), v(t-2h),and v(t-3h): 26 v(t) =a t + a t + a,t+a 3 2 3 2 0 v(t-h)=a (t-h) + a (t-h) + a,(t-h)+a 3 2 3 2 0 v(t-2h)=a (t-2h) + a (t-2h) + a,(t-2h)+a 3 2 3 2 0 v(t-3h)=a (t-3h) + a (t-3h) + a,(t-3h)+a 3 2 3 2 0 Solving this system o f equations, the coefficients a$, a.2, aj, and ao are given by _ -v(t-3h) + 3v(t-2h)-3v(t-h) + v(t) 6h _ -v(t -3h) + 4v(t -2h)~ 5v(t -h) + 2v(t) , • 3 (2.13) 2 h '~ a _ -2v(t-3h) + 9v(t-2h)6h a =v(t) 18v(t-h) + llv(t) 0 What is needed now is v(t-kh). This term comes from the interpolating polynomial o f equation (2.4) and the calculated values o f as, a2, aj, and ao, as follows v(t-kh) = --(-v(t-3h)+3v(t-2h)-3v(t-h)+v(t))k 3 6 +L(- ( -3h)+4v(t-2h)-5v(t-h)+2v(t))k • 2 v t (- ) 2 14 --(-2v(t -3h)+9v(t -2h)- J8v(t -h)+v(t ))k+v(t) 6 What is now left is to evaluate the approximate derivative using equation(2.10): v'(t) = (-(-v(t -3h) + 3v(t -2h)- 3v(t -h) + v(t ))k 6 3 -L(- ( -3h) v + 4v(t-2h)-5v(t-h) t +-(-2v(t-3h) 6 + 2v(t))k 2 + 9v(t-2h)-18v(t-h) - ( ) 2-15 + v(t))k)/kh If terms are collected from equation (2.15) the result is v'(t) = ((-k -k+^MO 6 o 2 (-k -2k+-)v(t-2h) 2 2 2 + (~k 2 2 + (--k 6 2 2 2 +lk-3)v(t-h)+ +-k--)v(t-3h))/h 3 Equation (2.16) is the difference equation as a function o f k. It is more convenient, especially for higher order interpolating polynomials, to write equation (2.16) as 27 (2.16) (2.17) where n is the order o f the used polynomial, in this case n—3 and k = (±k -k o 2 0 2 1 + ^), o 2 k =(--k +-k--). (2.18) 2 3 3 6 2 3 Notice that when k=\ equation (2.17) becomes v(0 v(t)-v(t-h) (2.19) h which is the Backward Euler formula, similar to equation (1.6). O n the other hand, when k=0, the third order B D F is obtained as v'(t)~{(^v(0-3v(t-h)+?-v(t-2h)-U ho 2 3 (2.20) 2.4 Magnitude and Phase Distortion Figure 2.2 shows the magnitude and phase distortion introduced by equation (2.17), computed as described i n section 1.2.1. A s expected, setting k to one results i n the same distortion as the Backward Euler formula shown i n Figure 1.4. Similarly, setting k to zero results i n the distortion introduced by the third order B D F . A n y intermediate values w i l l result i n a method with different distortion characteristics. In Chapter 3 o f this work, this ability to select the distortion characteristic o f the method by adjusting k is used to find methods tuned for a particular frequency, so that the resulting formula produces exact results at that frequency. 28 -10+ 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Frequency in per unit of (1/h) 0.40 0.45 0.50 ;ure 2.2 Magnitude distortion and phase error introduced by equation (2.17) for the indicated values o f k. 29 2.5 Stability Although the main focus o f this thesis is in the tuned difference equations o f next chapter and their applications, it is interesting to see how it is possible to adjust a difference equation so it becomes A-stable. For example, Figure 2.3 shows how the stability zone o f equation (2.17) changes with the value o f k computed as described in section 1.2.1. For instance, it is possible by carefully selecting the value o f k to obtain a formula which is A stable. It so happens that this the case for equation (2.17) when k is set to approximately 0.14, resulting in the difference equation v'(t) ~-{1.6966v(t)-2.6598v(t-h h -3 -2 -1 0 1 )+1.2298v(t-2h)-0.2666v(t-3h)). 2 Re(z) 3 4 5 6 7 Figure 2.3 Stability zone for equation (2.17) for the indicated values o f k. 30 (2.21) It is interesting to see the distortion introduced by equation (2.21). The solid line o f Figure 2.5 corresponds to the magnitude and phase characteristics o f equation (2.21) while the dashed line corresponds to the trapezoidal rule. The main characteristics o f the proposed method are A-stability and low error at low frequencies. The properties o f the method are very similar o f those o f the trapezoidal rule up to one tenth o f the maximum frequency 1/h. For instance, when the circuit o f Figure 2.4 is simulated using the trapezoidal rule, the output voltage o f phase ' a ' o f B u s l shows the numerical oscillations o f Figure 2.6, after circuit breaker C B 1 is tripped . 9 O n the other hand, when the circuit is simulated with equation (2.21), the output voltage at phase ' a ' o f B u s _ l is free o f such numerical oscillations, as shown in Figure 2.7. Figure 2.4 Circuit used to test equation (2.21) In Microtran, the UBC version of the EMTP, numerical oscillations are suppressed using the critical damping algorithm (CDA) as described in [2]. 31 -10+ 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Frequency in per unit of (1/h) 0.40 0.45 0.50 gure 2.5 Magnitude distortion and phase error introduced by equation (2.21). The dashed line corresponds to the trapezoidal rule, included here for comparison purposes. 32 _g I 0 1 0.01 1 0.02 I I I I I I I I 0.03 0.04 0.05 Time s 0.06 0.07 0.08 0.09 0.1 Figure 2.6 Voltage in phase ' a ' of B u s _ l when the trapezoidal rule is used to simulate the circuit of Figure 2.4. Notice the numerical oscillations when C B 1 is tripped at 60 ms. Time s Figure 2.7 Voltage in phase ' a ' of B u s _ l when equation (2.21) is used to simulate the circuit of Figure 2.4. 33 CHAPTER 3: Tuned Time-Domain Linear AC Steady State Solutions 3.1 Introduction Two possible techniques for the solution o f linear A C steady state circuits are phasor analysis and time-domain simulations using small time-steps. A n alternative to these techniques, discussed in this thesis, is the adjustable methods described in the previous chapter, which can be tuned-up to produce exact results at one frequency. To illustrate this finding, consider for example the magnitude and phase distortion introduced by the adjustable method o f order five o f Figure 3.1 plotted for k=0 to 0.4 in 0.1 increments. Notice from Figure 3.1 that the magnitude plot, for any o f the displayed values o f k, starts at one, raises slightly, and then decreases, crossing one again. Similarly, for the phase plot, the phase distortion starts at zero, but before reaching the final value o f 90°, it crosses zero twice. Furthermore, notice that for k=0.4 10 the magnitude plot crosses one and the phase plot crosses zero for almost the same value co o f the normalized frequency ! Not only 11 they are exact, but also the time step o f these methods is considerably larger than that o f traditional non-tuned time-domain methods. In this chapter, the derivation o f tuned methods is presented, as well as their application for the simulation o f one-frequency linear steady state electric power circuits. Later, in Chapter 4, tuned methods are used to solve the nonlinear static power flow problem, and in Chapter 5 to solve the dynamic power flow problem. 11 A more accurate valuefromTable 3.1 is £=0.4029346295063755. The normalizedfrequencyis the horizontal axis of Figure 3.1. It is defined as the product ay=f h. 34 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Frequency in per unit of 0.35 (1/h) 0.40 0.45 0.50 Figure 3.1 Magnitude distortion and phase error introduced by the adjustable method of order 5 for &=0.0 to 0.4 in 0.1 increments. 35 3.2 Tuning the Adjustable For Methods polynomials o f orders equal or greater than four, we can compute at least one value o f k and one value o f the normalized angular frequency co for which the magnitude distortion is one and the phase error is zero, for co different from zero. Since the normalized frequency is directly proportional to the time-step h (note 11, page 34), once co is found h is known. The values o f k and co can be found from the z-domain transfer function o f the corresponding difference equation. For B D F - l i k e methods, the transfer function o f the difference equation has the form (3.1) -jrm> I*. m=0 where n is the order o f the interpolating polynomial and k m are the difference equation coefficients as a function o f k, similar to those o f equation (2.18). A l s o , since e -jmo _ ^ j_j. j ( j^ cos mco s n ^ 2) m0d equation (3.1) becomes H(v>) = - — ^K (cos(m(x))m . (3.3) j • sin(mu))) m=0 Assuming that at one particular co there is no distortion, that is H(CO)=1,WQ have ^ K (cos(mco) m - j • sin(mco)) = y'co. (3.4) m=0 Separating real and imaginary parts, a system o f two equations and two unknowns is formed, 36 X m K C O S ( m = 0 = ™ ) g,(k,(X>) m=0 (3.5) n co + ^ K sin( wco ) = 0 = g (k,(x)) m 2 m=0 This system o f equations can be solved using Newton-Raphson iterations: dk "AA:" dto dg Aco 2 dk As (3.6) 2. dto _ where ^iL Y^s.cos(mo) dk to dk dk to dk = (3.7) = -Y km • sin( ma) OtO =0 M = / + ^fc m• cos( mcoj m Solving the system o f equations (3.5) for difference equations obtained with polynomials o f orders four to nine results in the values o f Table 3.1. n k CO 4 5 6 '7 8 9 0.1995183215439153 0.4029346295063755 0.5577432820878876 0.6761264770927172 0.7672635252838457 0.8369482736463172 0.9648838221480567 1.5328796706951471 1.9006879375994667 2.1526576264068528 2.3321297134871948 2.4636521899340823 Table 3.1 Tuned values o f k and OJ for the difference equations obtained with polynomials o f orders 4 to 9. 37 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Frequency in per unit of (1/h) 0.40 0.45 0.50 Figure 3.2 Magnitude distortion and phase error of tuned difference equations derived with polynomials of orders 4 to 9 and the values of k i n Table. 3.1. 38 39 From the tuned value o f co, the value o f the time-step h can be computed using h=—, (Q =2nf , <*>, s s (3.8) where f is the frequency at which the difference equation is required to produce exact results. s Figure 3.2 shows the magnitude and phase characteristics o f the tuned difference equations obtained with the values o f k in Table 3.1 for polynomials o f orders four to nine, while Figure 3.2 shows the stability regions for such methods, both obtained as described in section 1.2.1. From Figure 3.3 it can be observed that all the tuned methods are reasonably stable, but the tuned difference equation obtained with a seventh order interpolating polynomial is almost A-stable 12 [5] as shown in Table 3.2. The difference equation for this curve is v'(t) = (1.338v(t)-2.068v(t-h)+1.518v(t-2h)-1.441v(t-3h) +1.007v(t-4h)-0.465v(t-5h)+0.126v(t-6h)-0.015v(t-7h))/h. The coefficients o f equation (3.9) have been rounded to three decimal digits. More accurate coefficients or coefficients for different order interpolating polynomials can be obtained with the Maple 6 script o f Figure 3.4. n a(deg) 84.65845 86.96453 89.07581 89.98968 89.38110 87.13993 4 5 6 7 8 9 Table 3.2 The tuned methods obtained with interpolating polynomials o f order 4 to 9 are A(a)-stable. For the method obtained with a seventh order interpolating polynomial is almost A-stable. a for an A-stable method is 90°. 40 # ctune.m: F i n d the tuned c o e f f i c i e n t s o f the d i f f e r e n c e # e q u a t i o n u s i n g a p o l y n o m i a l o f o r d e r n. restart; Digits:=16: n:=7: # t h i s i s the o r d e r o f the p o l y n o m i a l to u s e . v:=array(0..n): a:=array(0..n): equ:=array(0..n): k:=array(0..n): # Create a p o l y n o m i a l o f the d e s i r e d o r d e r i : = ' i ' : i p o l y : = s u m ( ( a [ n - i ] * t ( n - i ) ) , i=0..n); A # Create the system o f . e q u a t i o n s f o r i from 0 to n do t:=-i*h; equ[i]:=v[i]=ipoly; end do: # S o l v e the system o f e q u a t i o n s assign(solve({seq(equ[i],i=0..n)},{seq(a[i],i=0..n)})): # C o n s t r u c t the d i f f e r e n c e e q u a t i o n t:=-K*h: diffeq:=simplify((v[0]-ipoly)/K*h): # From t h e d i f f e r e n c e e q u a t i o n get the formulas f o r i from 0 to n do v[i]:=0: end do: f o r i from 0 to n do v[i]:=1; k[i] :=sort(simplify(diffeq/h)); v[i]:=0; end do; f o r k[0..n] # D i s p l a y k[0..n] as a f u n c t i o n o f K. f o r i from 0 t o n do k[i]:=k[i]; end do; # Compute the tuned v a l u e s o f K and w. i : = ' i ' : eqR:=sum(k[i]*cos(i*wO), i=0..n)=0: i : = ' i ' : eql:=wO+sum(k[i]*sin(i*wO), i=0..n)=0: a s s i g n ( f s o l v e ( { e q R , e q l } , {K, wO}, {K=0.15..1.0, w0=0.9..evalf(Pi)})); # D i s p l a y the r e s u l t s K:=K; wO:=wO; f o r i from 0 to n do k[i]:=k[i]; end do; Figure 3.4 This Maple 6 script can be use to computed the coefficients o f the adjustable difference equations, as well as to find their tuned values for interpolating polynomials o f order > 4. 41 3.3 Circuit Simulation Example Before using equation (3.9) to simulate a circuit, the value o f h, the time-step, must be computed. Suppose a circuit with a steady state solution at 60 H z has to be solved. Using equation (3.8) for n=7 results in , 2.1526576264068528 h= ., = 5.7101017004519526ms c n i n i n 1 n n n A C i n e . (3.10) 2%60 The circuit o f Figure 3.5 was simulated using difference equation (3.9) and the time step given by equation (3.10). The simulation was started using A C steady state initial conditions computed with phasor analysis. The voltage across inductor L I is presented in Figure 3.6. The solid line represents the voltage obtained using phasor analysis, while the scattered dots are the voltages obtained after a time domain simulation performed with equation (3.9). The dots match perfectly with the corresponding values obtained using standard phasor analysis. Figure 3.5 Circuit used to test equation (3.9). 42 0.00 0.02 0.04 0.06 0.08 0.10 Time (s) 0.12 0.14 0.16 0.18 0.20 Figure 3.6 The solid line is the voltage i n the inductor o f Figure 3.5 obtained using phasor analysis. The dots are the values obtained when simulating the same circuit using the tuned difference equation (3.9). 3.4 Reducing the Number of Circuit Nodes It w i l l be shown later in this chapter, that the cost o f the tuned time-domain simulation is proportional to the cube o f the number o f nodes i n the circuit to be solved. Therefore, it is important to reduce the number o f nodes in the circuit. The most important reduction is achieved when series passive elements are discretized together so that the internal nodes connecting them are unnecessary. In particular, the series connection o f resistors and inductors ( R L ) as well as resistors and capacitors ( R C ) can be reduced. This is also necessary for the impedances represented with complex numbers found in phasor analysis, which need to be converted either to a R L or R C equivalent circuit for time-domain 43 simulations. For instance, a complex impedance o f the form R+jX is converted to a R L series circuit, where L=X/(2itf). If the imaginary part o f the complex impedance is negative, that is R-jX, then it is converted to a R C series circuit, where C=l/(2jtpQ. The R L and R C circuits' discrete representation is similar to that o f Figure 1.2, and depicted in Figure 3.7. For the R L series connection, the voltage o f the series is given by v(t) = Ri(t) + L ^ - . dt (3.11) Replacing the derivative with equation (2.17) results in v(t) = Ri(t) + ^ £ k i(t - mh) . m " (3.12) m=0 Collecting the first term o f the current i(t) results in the formula we are looking for v(t)= („ Lk R + —t n i(t) + -Y k i(t-mh), d n m (3.13) m=l from which the values o f the equivalent circuit are SRL Lk, R +- 0 h (3.14) h =-^-f k i(t-mh) RL J m n m=l For the tuned time-domain steady state solution o f a given circuit, the coefficients ko to k„ o f equation (3.14) would be computed using the values o f k o f Figure 3.1. 44 hRL + gRL hRC + gRC Figure 3.7 Discrete equivalents for R L and R C circuits. Similarly for the R C circuit the voltage across the series combination is given by: vft) = Rift) + — lift )dt. C (3.15) J Taking the derivative o f equation (3.15) results in dt dt (3.16) C Replacing the derivatives with equation (2.17) results in 1 " " m=0 R " 1 n C m=0 (3.17) Collecting the vft) and ift) terms from equation (3.17) yields f k vft) + Y k fvft-mh)-Rift-mh)) 0 J = m k R 0 + h — HO, (3.18) J m=l From which the discrete equivalent circuit values are: SRC ~~ KLR 0 + C (3.19) =f-Y k (vft-mh)-Ri(t-mh)) Rc h J 45 m 3.5 Estimating Magnitude and Phase To compute the A C steady state o f a circuit using any o f the tuned difference equations such as equation (3.9), the magnitude and phase o f each node in the circuit must be calculated using the discrete points resulting from the transient simulation. Since equation (3.9) is not only exact for the tuned frequency but also for D C , it is possible, by using three previously computed points to match a function o f the form v(t) = V + dc 42V ms cos(2xft + Q) (3.20) by solving the system o f equations v(t ) = V + V cos(cdt + Q) dc v(t-h) P = V + V cos(oi(t-h) dc + B) P , (3.21) [v(t -2h) = V + V cos(u(t -2h) + Q) dc P where v(t), v(t-h), and v(t-2h) are the simulation results at times t, t-h, and t-2h, and (3.22) u> = 2xf A l s o , by making the following substitutions in the systems o f equations (3.21) cos( <x>(t-h) + B) = cos( oit+Q) cos( u>h) + sin( oit+Q) sin( 00/2) cos( <j)(t-2h) + Q) = cos( cat + Q) cos( co2h) + sin( oit+Q) sin( o)2h), V =v(t)- (3.23) V cos((ot + Q) dc p we obtain in the reduced system o f equations \v( t-h) \v(t-2h) = v(t)-V p cos( oit + Q) + V (cos( oit + Q) cos( wh) + sin( cor + 6 j sin( (ah) p = v(t)-VpCos((tit + Q) + Vp(cos((X)t + Q)cos(o)2h) + sin((Dt+Q Collecting terms in (3.24) results in \v(t) - v(t - h) = Vp cos((ot + Q)(l + cos((x)h)) + V sin((at + Q)sin(u>h) p \v(t)-v(t-2h) = V cos( (x)t+Q)(l+cos( w2h)) + V sin( cat + Q) sin( (a2h) P P 46 (3.24) Applying this substitutions to (3.25) Nj = V sin((tit + Q),N =V cos((s}t p 2 + Q), p (3.26) yields the linear system o f equations v(t)-v(t-h) = N sin((x)h) + N (l+cos((x>h)) 1 vft )-v(t-2h) 2 = N, sin((x)2h ) + N fl 2 (3.27) + cos(o)2h )) Solving the system o f equation (3.27) for Nj and N2 allows for the computation o f V , 8, and p V , as dc Q= V dc atan(^-)-i£>t =v(t)- (3.28) V cos fat + Q) P 3.6 Tuned Time-Domain Steady State Simulations The algorithm for the simulation o f an electric circuit using the tuned difference equations is similar to the one used i n the Electromagnetics Transient Program E M T P described in [30]. The flow chart o f Figure 3.8 shows the program sequence, which includes most o f the steps presented in section 1.2.1. The main difference, besides the use o f the tuned formulas, is the computation o f the A C magnitudes and phases for the node voltages as described in section 3.5. It is the change o f the A C magnitudes and phases after each time step what is used to terminate the simulation. If the desired tolerance is not achieved in a preset number o f steps, the simulation is terminated and an error message is displayed. A n electric power simulation program, which follows the flow chart o f Figure 3.8, was developed using the C programming language. The program is not only capable o f solving the linear steady state o f a circuit, but also the static and dynamic power flow as well. A l s o , 47 the discretization techniques used in the program are selectable, so it was possible to test the tuned difference equations generated with different interpolating polynomials. Begin Update Sources Read Case Solve Node Voltages Check Case Compute Currents Compute IMag. & Pha 1 Report Errors Build Y Matrix Steady\Yes State' 1 Report Results Abort LU Factor Y Matrix Yes Store Results Finish Figure 3.8 Program flow chart to compute the steady state with the tuned formulas. The simulation is stopped once the steady state is found, or a maximum number o f time steps is reached. 48 A n y o f the tuned difference equations obtained with the values o f Table 3.1 can be used to compute the steady stafe o f linear A C circuits. Since a simulation performed with any o f these methods is a transient simulation, a number o f time steps have to be executed in order for the transients to damp. The main advantage these tunable methods have over other discretization methods is their considerably larger time step. For instance, to solve the 60 H z steady state problem with the trapezoidal rule with an error o f 1%, a time step o f 0.5 ms is required. The same problem can be solved using equation (3.9) with a time step o f about 5.71 ms with no error. That is, the trapezoidal rule would require at least ten times more simulation steps than equation (3.9), even though it is less accurate for a steady state solution. In practice this difference is even larger, since equation (3.9) w i l l introduce additional damping for any frequency different from the tuned frequency as shown in Figure 3.10. A s an example, consider the circuit o f Figure 3.9, for which the A C steady state is computed. When this circuit is simulated using equation (3.9), the resulting node voltages, using the fitting function (3.20) on the discrete points obtained, are shown in Table 3.3. The steady state was obtained in 278 time steps using as stop criteria a change i n 6 o f less than 10" degrees. 5 These results are identical to those obtained with the steady state option o f Microtran ( U B C version o f the E M T P ) . Notice from Figure 3.2 that the tuned methods not only produce exact results for the tuned frequency co, but also for D C . For instance, when a direct current voltage source o f 10 V is connected i n series with V I , the D C values o f Table 3.3 are obtained. Therefore, the methods are capable o f computing the D C steady state i n addition to the A C steady state. See equation (3.10). 49 Figure 3.9 Circuit from which the steady state is computed using equation (3.10). Inductances are in m H and capacitances are i n uJF. Node V d c (V) V r m s (V) #(deg) A 10.0000 7.071068 0.0000 B 10.0000 11.47991 -16.9305 C 9.99700 -17.1180 11.46359 D 9.99700 24.92909 -31.3557 Table 3.3 Steady state voltages o f the circuit in Figure 3.9. The D C values are obtained when a 10 V dc voltage source is connected i n series with V I . Figure 3.10 shows how the magnitude and phase o f the voltage phasor at node ' D ' o f the circuit o f Figure 3.9 change over time when a 60 H z source is applied and the seventh order tuned method is used (equation (3.9)). After an initial transient, the voltage phasor quickly converges to the steady state solution. O n the other hand, Figure 3.11 shows the magnitude and phase for the voltage at the same node when the circuit is simulated using the trapezoidal rule with a time step o f 166.67 us. Since the trapezoidal rule does not introduce any damping by itself, all the circuit natural frequencies remain in the simulation outputs for a long time. After 10000 time steps the simulation using the trapezoidal rule was aborted, as indicated by the program flowchart o f Figure 3.8, since the simple technique described in section 3.5 to compute the magnitude and phase o f voltages fails to produce an acceptable result. If the trapezoidal rule were to be used, some kind o f additional filtering would be required instead, as suggested in [31]. 50 0.8 Time (s) 1 Figure 3.10 Magnitude and phase change over time for the voltage at node ' D ' o f the circuit of Figure 3.9 for A C excitation using the tuned formula o f order 7. 51 15 10 52 3.7 Computational Cost ofthe Tuned Time-Domain Methods The time required to compute the steady state of a linear circuit using the tuned formulas may be smaller than the time required when using phasor analysis. For instance the number of numerical operations required to L U factor a n x n real matrix is given approximately by [32]: #mult = -n 3 + 3 0(n ) 2 # adds / subs = l-n + 0(n ) 3 2 (3.29) #divs- — n + 0(n) 2 2 and to solve the system of equations using back substitution: BS = 2n -n~2n . (3.30) 2 Therefore if the tuned formulas are used to get the steady state, and it takes q steps to damp all non-tuned harmonic components, the total number of operations is given by Sol(real) » | n + 2qn • 3 (3.31) 2 On the other hand, if the circuit is solved using phasor analysis, and we consider that each complex addition or subtraction takes two real additions or subtractions, each complex multiplication takes four real multiplications and two real additions, and each complex division takes six real multiplications plus two real divisions and three real additions, the number of operations would be: #mult = -n 3 3 + 0(n ) 2 # adds/subs = -n + 0(n ), 3 #divs = n + 0(n) 3 2 53 2 (3.32) or ex)~—n 3 (3.33) The breakpoint at which the tuned time-domain and phasor solutions start to take approximately the same number o f operations is given by: n =>q~n. (3.34) That is, when the number o f step evaluations is smaller than the number o f nodes i n the circuit, the time required to obtain the solution using the tuned formulas w i l l be less than the time required solving the same circuit using phasor analysis. To verify this, two programs were written in C/C++ to solve two standard I E E E cases (with 118 and 300 nodes), one using phasor analysis, the other using the tuned formula derived from a seventh order polynomial (3.9) as described in section 3.6. Both programs were compiled with Microsoft Visual C++ 6.0, with the same speed optimizations, and using the same algorithm to solve the system o f equations: L U decomposition with back substitution. The stop criterion for the proposed tuned formula was a change o f less than 0.0001 degrees in the final node angles. The results are shown in Table 3.4 for five different processor types and speeds. Notice that the tuned method solution is faster in all cases since the required accuracy was obtained with a number o f evaluation steps q smaller than the number o f nodes, as suggested by equation (3.34). Please notice that matrix sparsity techniques were not used, but the above comparison remains valid even i f a different solution technique is applied to both the phasor analysis and tuned time-domain methods. 54 CPU type 118 nodes Phasor Tuned (^=78) 0.93s 0.6s 300 nodes Phasor Tuned (4=150) 9.34s 6.49s Ppro 200MHz PII 0.34s 0.21s 4.026s 400MHz Celeron 0.45s 4.797s 0.38s 1GHz Athlon 0.22s 0.15s •2.553s 900MHz PHI 0.16s 2.293s 0.21s 866MHz Table 3.4 Steady state solution times for phasor analysis and domain formula. 55 2.673s 3.655s 1.512s 1.852s a seventh order tuned time- CHAPTER 4: Solving the Static Power flow Problem 4.1 Introduction Power flow analysis is the fundamental tool for many electric power system studies. It is the primary tool used for contingency analysis and clearance evaluation. It is also required for transient and dynamic stability studies, as well as to generate the base cases for security analysis and transfer capability studies. The most frequently used technique for power flow analysis consists in the solution o f a system o f non-linear equations using iterative methods such as the Newton-Raphson method, described i n section 1.2.4, or one o f its derivatives, such as the decoupled power flow [26;33]. In this chapter, the difference equations used i n the tuned time-domain solution o f the steady state problem are used to solve the non-linear static power flow problem. Before proceeding with the solution o f such non-linear power flow problems, suitable models for the elements commonly found in this kind o f analysis, such as P Q loads and P V developed. 1 4 sources, are Then the new tuned time-domain power flow approach is used to solve three standard tests systems, and the results are compared with those obtained using NewtonRaphson iterations. Finally, the developed tuned time-domain power flow is used to obtain detailed P V curves . 15 As they are connected to the PQ and PV buses described in section 1.2.4. A PV curve shows the variation of a bus voltage as a function of the reactive power, at either that bus or any combination of buses in the power system. They are useful to asses the stability ofthe power system. 15 56 4.2 Modeling the Power System Elements In order to solve the non-linear power flow problem, adequate circuital models for P Q loads, P V sources, transmission lines, and transformers are required. For transmission lines and transformers the standard models available for iterative power flow solution are used. P Q loads and P V sources are modelled using current and voltage sources respectively. The swing generator is modelled i n time domain simulations by using an independent voltage source. Other power flow elements such as phase shifters and loads different from the P Q load can be easily modelled as well. However, since they rarely appear i n standard power flow cases, they w i l l not be described in this chapter. 4.2.1 Modeling Transmission Lines Transmission lines are modelled in the tuned time-domain power flow simulator in a similar way as they are treated i n the traditional iterative methods. To this effect the equivalent jr. circuit o f Figure 4.1 is used. The circuit is composed o f four elements: a) The series resistance R, representing ohmic losses and skin effect. b) The shunt admittance G , representing leakage currents and corona effect. This value is so small that it is normally neglected in power flow analysis. For example, it is not even found in the I E E E common format for exchange o f solved power flow data [34]. If used, it is represented as two admittances o f half value at both ends o f the transmission line. c) Series reactance X , representing internal and external magnetic flux linkages o f the transmission line conductors. For the tuned time-domain power flow, the series reactance X and series resistance R are represented as one discrete R L element as outlined in section 3.4. 57 d) Shunt susceptance B , representing the voltage charge between conductors and ground. It is modelled as two shunt capacitances o f half value at both ends o f the transmission line. R X Figure 4.1 Equivalent Jt circuit o f transmission line. 4.2.2 Modeling Transformers Typically two-winding transformers i n power flow studies are represented by the top circuit o f Figure 4.2. This model consists o f an equivalent branch impedance and an ideal transformer with the desired turns ratio ' a ' . To simplify its inclusion in the admittance matrix required in nodal analysis, the equivalent circuit at the bottom o f Figure 4.2 is normally used [35]. The value o f the equivalent admittances are given by (4.1) Y2 = Y3 = A s in the case o f the transmission line, the resulting complex admittances o f Figure 4.2 are represented by their equivalent R L circuit. For the static power flow solution, a fixed-tap transformer model is used. Later, in section 5.2, the under-load tap-change transformer ( U L T C ) model is described, since it is needed for dynamic power flow simulations. 58 Figure 4.2 Typical two winding transformer model used i n power flow studies (top) and its equivalent circuit (bottom). The values o f Y l , Y 2 , and Y 3 are given by equation (4.1). 4.2.3 Modeling Shunts Shunts are admittances connected in parallel to the system buses. capacitive or inductive, and are modelled as a capacitor or inductor. I E E E common data format for power flow cases [34], They are either Accordingly to the capacitors are represented as a positive shunt reactance value +B, while inductors are represented as a negative reactance value - B . 4.2.4 Modeling P Q Loads One way o f modeling P Q loads in time-domain simulations consists o f an adjustable sinusoidal current source whose magnitude and phase changes over time, until the specified complex power is obtained. This is accomplished by using the solved voltage at the node where the load is connected. For instance, from the definition of complex power: p+Qj = Vxl'. 59 (4.2) If magnitude and phase are used instead, (4.3) O n time-domain simulations the peak values o f currents and voltages are used, therefore (4.4) from which the magnitude and phase o f the current source used to represent the load can be computed as (4.5) (3 = a - 9 A very important factor to take into consideration for this load model to work is the numerical transient that occurs every time there is a change in the circuit. The characteristics of this transient depends both on the discretization technique used and how much the circuit changes [2]. For the formulas proposed in this report the numerical transient damps completely after the number o f time steps, counting from the last circuit change, is equal to the order o f the method used. This numerical transient is illustrated i n Figure 4.3 for a nontuned formula derived from a seventh order interpolating polynomial. For example, i f the tuned formula obtained with a seventh order interpolating polynomial is used, then the current source should be adjusted every seventh time step. To test the P Q model, the simple circuits o f Figure 4.4, obtained from [33] are used. The circuit on the left produces a stable result as shown in Figure 4.5, while the circuit on the right, with the load values suggested in [33] , is i l l conditioned and the voltage collapses, as 16 it is also shown in Figure 4.5. 16 p. 1010. Q was slightly incremented from 0.15 to 0.152 to accelerate voltage collapse. 60 12.5 12 - 0.0245 0.025 0.0255 Time s 0.026 0.0265 Figure 4.3 Damped numerical oscillations caused by a non-tuned formula derived from a seventh order interpolating polynomial. Similar behaviour is observed with the tuned formulas. V1 V2 Phase=0.0 Phase=0.0 V: 0.855373 Ph:-13.5219 Voltage=1.0 V: 0.482312 Ph:-55.2707 Voltage=1.0 X1 0.25 X3 0.25 X4 X2 r P=0.8 Q=0.4 1.843 =0.152 Figure 4.4 Circuits used to test the P Q load model. The circuit on the right yields a solution, whereas the circuit in the left is ill conditioned and the voltage at node ' B ' collapses. 61 1.3r A B 1.2k 1.1 E 0.81- < 0.7 0.61 0.5 0.4 L 10 Time (s) 12 14 16 18 20 ou 20 10 0 '-10 1. -20 -30 -40 il 1 r — VC 1 -50 -60 i- 10 Time (s) 12 14 16 18 20 F i g u r e 4.5 M a g n i t u d e and phase changes o v e r t i m e for the circuits o f F i g u r e 4.4. 62 4.2.5 Modeling PV Sources P V sources are modelled with a sinusoidal voltage source in the tuned time-domain simulator. determined. In this case, the magnitude o f the voltage is known and the phase has to be Since the current coming out o f the sinusoidal voltage source is needed to compute the active power, the modified nodal analysis ( M N A ) model o f such sources is used [27]. The use o f the M N A voltage source augments the size o f the circuit Y matrix by one row and column as explained in section 1.2.2. The steady state is obtained by changing the phase o f the voltage source until the desired active power comes out o f the voltage source. From formula (4.4) 9 = a-(3 (4.6) and the actual active power coming from the voltage source is given by P c uai=0.5*V *I cos(a-$) a t p (4.7) p In equation (4.7), a is slowly changed so that Pactual gradually becomes equal to the required active power P. A proportional-integral controller (PI), as shown in Figure 4.6, is well suited for this purpose. Typical values o f the proportional constant Kp range from 0.001 to 0.5, while the integrative constant Ki, using five terms in the summation, can go from 0.0001 to 0.05. The larger the values o f Kp and Ki the faster the convergence to the steady state, but the risk o f the controller going unstable increases. For most cases a Kp of 0.01 and a Ki o f 0.001 seems to work satisfactorily. Consider for example the circuit o f Figure 4.7 taken from [33]. For Kp=0.\\ and £7=0.011 the steady state is obtained with a tolerance o f 10" i n 332 steps, while the 5 magnitude and phases converge as depicted in Figure 4.8. If Kp=\ and Ki=0A the number o f steps required to achieve steady state is 3328 time steps, and the magnitude and phase 63 waveforms are those o f Figure 4.9. Similarly for Kp=0.01 and £7=0.001 it takes 2124 time steps to achieve the steady state, with the magnitude and phase changing as shown in Figure 4.10. A s observed from Figure 4.9, large Kp and Ki constant values results in faster transition near the final steady state, but the response is somehow noisy. O n the other hand, small constant values result in a much smoother response as shown i n Figure 4.10, but it takes a longer time to achieve the steady state. It is therefore o f crucial importance to select appropriate values for these constants. A rule o f thumb consists in using the total generated active power and the total amount o f load active power to compute the PI controller constants as Kp = 10\Pg~Pl\, Ki = (4.8) Kp 10 where Pg is the total specified active power generated at the P V buses, and PI is the total specified load active power both at the P V and P Q buses. A l s o , notice that the constants are selected to accelerate convergence not to accurately model the dynamic behaviour o f the P V generators. Capacitor banks and other reactive power sources are also often modelled as P V generators. In this case the active power is set to zero and the capacitor bank is expected to maintain the desired bus voltage providing any necessary reactive power in the process. Reactive power limits can be taken into account by monitoring the reactive power generated by the P V source. If the limit is reached, the phase o f the voltage source is computed using a = arctan(Qij2«.) + (3, 64 (4.9) and the voltage source magnitude is computed using y QumU ( 4 1 Q ) 0.5*I sin(a-$) p p Equations (4.9) and (4.10) are functionally equivalent to the P Q load described i n section 4.2.4. Figure 4.6 PI controller used to adjust the phase angle o f the P V sources. Bus-3 Bus-2 I P=0.9 V=1.0 Bus-1 Phase=0 X'=0.5 Voltage=0.9 X=0.15 ho- V: 1.000 Ph: 28.37 X'=0.93 -oH V: 0.944 Ph: 20.15 Figure 4.7 Circuit used to test the P V sources model. 65 />. '^v., BUS-1 BUS-2 l l l l 0 0.2 0.4 0.6 0.8 1 Time (s) 1.2 1.4 1.6 1.8 2 1 20 f—' J- s '. r r r J 1 1 °- 10 r-1 r 1 BUS-1 > BUS-2 J . .J... 5 jj J. 1 r \r 0 1 0 0.2 0.4 0.6 0.8 1 Time (s) 1.2 1.4 1.6 1.8 2 Figure 4.8 Voltages magnitude and phase for the circuit o f Figure 4.7 when Kp=0.11 and Ki=0.011 66 j 111111.. t h Wit \ BUS-1 BUS-2 0.85 0 2 4 6 8 10 Time (s) 12 14 16 18 20 30 25 ... i. v , „ 20 J [ IJ 11 11 ll BUS-1 BUS-2 1 1 I II II II Jl. 0 0 2 4 6 8 10 Time (s) 12 14 16 18 20 ;ure 4.9 Voltages magnitude and phase for the circuit o f Figure 4.7 when Kp=1.0 and Ki=0.1 67 | 0.95 E < BUS-1 BUS-2 0.9 0.85 10 12 14 Time (s) Time (s) Figure 4.10 Voltages magnitude and phase for the circuit of Figure 4.7 when Kp=0.01 and Ki=0.001 68 4.3 Tuned Time-Domain Static Power flow Test Cases To validate the functionality and accuracy o f the tuned formulas and techniques described, the models o f section 4.2 were added to the simulator program described in section 3.6 and applied to three test cases: the standard I E E E 14-bus and 57-bus test systems, 17 and the C I G R E 32-bus test system. The tests cases were solved both using tuned time- domain power flow and Newton-Raphson iteration power flow, and the voltages magnitude and phase resulting from both methods compared. For the tuned time-domain power flow the seventh order tuned formula was used for all the cases. The tolerance for all the cases and methods was set 10" (see section 3.6). 5 4.3.1 IEEE 14-Bus test system The standard I E E E 14-bus test system o f Figure 4.11 [36] consists o f 14 buses; five generator buses o f which one is the swing generator and the other four are P V sources; eleven P Q loads; fifteen transmission lines; five transformers, o f which three are off-nominal tap transformers; and one shunt reactance. K p and K i for all the P V generators were set using equation (4.8). The final voltages magnitude and phase are shown in the circuit o f Figure 4.11 in the box close to each bus, while the same voltages magnitude and phases are shown in Figure 4.12 as they change over time. The results obtained with the tuned timedomain program are identical to those obtained using the Newton-Raphson iterations program. Conseil International des Grands Reseaux Electriques. 69 ni JTJ-DW — Q.tr_ •2 c £ £ « S 0 > |2 S'cocnScooS Figure 4.11 I E E E 14-bus test system. 70 -5 -10 -15 -20 0.5 1.5 Time (s) 2.5 Figure 4.12 Voltages magnitude and phase change over time for the I E E E 14-bus test system solved using tuned time-domain power flow. 71 4.3.2 CIGRE 32-Bus test system The static power flow o f the 32-bus system shown i n Figure 4.13 [37] was also obtained using the tuned time-domain static power flow. It consists of: the swing generator, four P V buses, 19 P Q loads, 25 lines, 15 tap transformers, and 12 shunts. Table 4.1 shows the buses magnitudes and phases obtained both using Newton-Raphson i n three iterations and the tuned time-domain simulation, obtained i n 525 time steps. Except for the last digit i n some o f the buse phases, the results obtained with both methods are identical. Bus N201 Newton-Raphson Tuned Time-Domain IVI Phase IVI Phase 1.0426 -9.6963 1.0426 -9.6963 -11.3695 N202 1.0236 -11.3695 1.0236 N203 1.0247 -9.9861 1.0247 -9.9861 N204 1.0134 -11.3138 1.0134 -11.3138 N205 1.0332 -9.7571 1.0332 -9.7571 N206 1.0256 -10.1259 1.0256 -10.1259 N207 1.0254 -11.6906 1.0254 -11.6906 N101 1.0547 -5.7613 1.0547 -5.7613 N102 1.0286 -7.1895 1.0286 -7.1895 N103 1.0193 -7.0523 1.0193 -7.0523 N104 1.0296 -5.8108 1.0296 -5.8108 N105 1.0552 -5.8521 1.0552 -5.8521 N106 1.0305 -5.9568 1.0305 -5.9568 N107 1.0303 -7.5213 1.0303 -7.5213 Ml 1.0000 4.1614 1.0000 4.1614 M2 1.0000 1.5591 1.0000 1.5592 Nl 1.0502 -0.1696 1.0502 -0.1696 N10 1.0840 2.5467 1.0840 2.5468 Nil 1.0881 0.8021 1.0881 0.8021 N13 1.0593 -0.3546 1.0593 -0.3546 N14 1.0814 0.8045 1.0814 0.8046 N16 1.0400 -1.6802 1.0400 -1.6802 N2 1.0485 -0.3843 1.0485 -0.3843 N3 1.0398 -1.6628 1.0398 -1.6628 N4 1.0442 -0.9636 1.0442 -0.9635 N5 1.0400 -1.4283 1.0400 -1.4283 N6 1.0443 -1.5498 1.0443 -1.5498 N7 1.0303 -2.0440 1.0303 -2.0440 N8 1.0552 -0.6752 1.0552 -0.6751 N9 1.0554 T 0.6351 1.0554 -0.6351 N12 1.0977 -1.3056 1.0977 -1.3055 N15 1.0920 0.0000 1.0920 0.0000 Table 4.1 Node voltages for the system o f Figure 4.13, solved using Newton-Raphson iterations and tuned time-domain power flow. Reference [37] includes this case in the common data format described in [34], 72 Figure 4.13 32-bus test system. The resulting voltages magnitudes and phase for this system are presented in Table 4.1. 73 4.3.3 IEEE 57-Bus test system Table 4.2 shows the voltages magnitude and phase for the buses o f the standard I E E E 57-bus system shown in Figure 4.14 and Figure 4.15 [36]. The system consists of: the swing generator, six P V generators, 42 P Q loads, 63 lines, 17 tap transformers, and 3 shunts. The results obtained both with Newton-Raphson i n three iterations and the tuned time-domain static power flow are identical. The time-domain solution was obtained in 528 time steps. Bus BUS_1 Newton-Raphson IVI Phase BUS 2 BUS_3 BUS_4 1.0400 1.0100 0.9850 0.9808 BUS 5 BUS_6 0.9765 0.9800 BUS_7 0.9842 -7.6014 BUS_8 1.0050 -4.4779 BUS_9 BUS_10 0.9800 0.9862 -9.5847 -11.4497 1.0050 0.9800 0.9862 BUS_11 0.9740 -10.1932 0.9740 BUS 12 BUS_13 BUS_14 1.0150 -10.4712 0.9789 0.9702 -9.8035 -9.3503 1.0150 0.9789 0.9702 0.9880 1.0134 -7.1902 0.9880 -9.3503 -7.1902 -8.8589 1.0134 -8.8589 -5.3959 1.0175 -5.3959 BUS_18 1.0175 1.0007 -11.7296 1.0007 BUS_19 0.9702 -13.2265 0.9702 BUS 20 BUS_21 0.9638 -13.4443 -12.9290 BUS_15 BUS_16 BUS_17 1.0085 0.0000 -1.1882 -5.9881 -7.3374 -8.5464 Tuned Time-Domain IVI Phase 1.0400 0.0000 -8.6741 1.0100 0.9850 0.9808 0.9765 0.9800 0.9842 Bus Newton-Raphson IVI Phase BUS_30 0.9627 -1.1882 -5.9881 -7.3374 BUS 31 BUS_32 BUS_33 0.9359 0.9499 0.9476 -8.5464 -8.6741 BUS 34 BUS_35 0.9592 0.9662 -18.5520 -14.1490 -13.9062 -7.6014 -4.4779 BUS_36 BUS_37 0.9758 0.9849 -13.6348 -13.4459 -9.5847 BUS_38 BUS_39 -12.7346 -13.4910 -13.6582 -18.7196 -19.3838 -18.5123 Tuned Time-Domain IVI Phase 0.9627 -18.7196 0.9359 0.9499 0.9476 0.9592 0.9662 0.9758 -19.3838 -18.5123 -18.5520 -14.1490 -13.9062 0.9849 -13.6348 -13.4459 1.0128 0.9828 -12.7346 -13.4910 0.9728 0.9962 -13.6582 0.9665 -15.5328 -11.3544 BUS_40 1.0128 0.9828 0.9728 -10.4712 BUS_41 0.9962 -14.0767 -9.8035 BUS 42 BUS_43 0.9665 1.0096 -15.5328 -11.3544 BUS 44 1.0168 -11.8565 1.0096 1.0168 BUS_45 1.0360 -9.2701 1.0360 -9.2701 1.0598 -11.1161 1.0598 -11.7296 BUS_46 BUS_47 BUS_48 -12.5116 -12.6107 1.0333 -13.2265 1.0333 1.0274 -11.1161 -12.5116 1.0274 -12.6107 0.9638 1.0085 -13.4443 BUS_49 1.0362 -12.9361 1.0362 -12.9290 BUS_50 1.0233 -13.4127 1.0233 -12.9361 -13.4127 -12.5334 -11.4497 -10.1932 -14.0767 -11.8565 BUS_22 1.0097 -12.8743 1.0097 -12.8743 BUS_51 1.0523 -12.5334 1.0523 BUS_23 1.0083 -12.9396 1.0083 -12.9396 BUS_52 0.9804 -11.4976 0.9804 -11.4976 BUS_24 0.9992 0.9992 -13.2921 0.9825 0.9825 -18.1732 BUS_53 BUS_54 0.9709 0.9963 -12.2526 -11.7097 0.9709 BUS_25 BUS_26 -13.2921 -18.1732 0.9963 -12.2526 -11.7097 0.9588 -12.9813 0.9588 -12.9813 BUS_55 1.0308 -10.8011 1.0308 -10.8011 BUS 27 0.9815 -11.5136 0.9815 -11.5136 BUS 56 0.9684 -16.0651 0.9684 -16.0651 BUS_28 0.9967 -10.4816 0.9967 -10.4816 BUS_57 0.9648 -16.5837 0.9648 -16.5837 BUS_29 1.0102 -9.7718 1.0102 -9.7718 Table 4.2 Static power flow voltages solution for the circuit o f Figure 4.14 and Figure 4.15. 74 OO i- o do II » ta c3 « m •2 I— ^icc £ ' a -o 05 •- Q.IT c o m w x Q.p<i! c o a S o §>c °>CO W S t 0 O 5 f OJ O CM o •<- o odd T-O do CLO o o o a: x co Figure 4.14 I E E E 57-bus test system, sheet 1 o f 2. 75 Figure 4.15 I E E E 57-bus test system, sheet 2 o f 2. 76 4.4 Computational Cost of the Tuned Time-Domain Power flow A t this point, it is possible to compare the computational cost o f the tuned timedomain simulation against the standard Newton-Raphson iterations power flow. As mentioned in section 1.2.4, the size o f the Jacobian matrix when using Newton-Raphson iterations is given by NR = 2*PQ + PV, (4.11) where PQ is the number o f P Q buses, and PV is the number o f P V buses. Assuming that it takes IT iterations to reach the desired tolerance, the number o f floating point operations required to solve the system using Newton-Raphson iterations is SolNR ~IT*^*(NR) . 3 (4.12) O n the other hand, the size o f the Y matrix for the time domain simulation using the tuned formulas and M N A voltage sources for the P V elements and swing generator is given by TD = PQ + 2*(PV + 1J (4.13) while the number o f floating point operation required to reach the steady state is given by SolTTD~^(TD) 3 + 2q(TDf, (4.14) where q is the number o f time steps required to reach the steady state. Therefore, the tuned time-domain simulation could result in a considerable smaller matrix to be L U decomposed i f the number o f P Q buses is much larger than the number o f P V buses. Unfortunately, it is just the opposite i f the P V buses outnumber the P Q buses. However, even i f the number o f P Q buses is larger than the number o f P V buses, the tuned time-domain may still result in a larger number o f operations i n some cases. For example, for the system o f Figure 4.13, which has 27 P Q buses and 4 P V buses, NR=5&, while TD=2>1. The solution using Newton- 77 Raphson converges in IT—3 iterations. The time domain solution reaches the steady state in q=535 steps. Therefore the approximate number o f floating point operations for each method is: SolNR *>3*-*(58f = 390224 3 SolTTD = ^(37/ . (4.15) +2 * 535* (37/ = 1465743 Therefore, the tuned time-domain solution, for the system o f Figure 4.13, takes four times more floating-point operations. This is often the case for small cases, as presented in Table 4.3. Fortunately, the number o f time steps required by the tuned method depends on the time constants o f the system not on its size and the situation would change radically for a larger system. For instance, i f P^J^OOO and PF=1000, and assuming IT=3 and g=1000 steps, results on aSo/Affi=6.859xl0 and Sb/Y7Y)=1.130xl0 . 12 12 Case TTD NR IEEE-14 0.02 s O.001 s CIGRE-32 0.03 s 0.03 s IEEE-57 0.08 s 0.15 s Table 4.3 Power flow solution times using the tuned time-domain and NewtonRaphson methods. 4.5 Generation of PV Curves Despite its computational cost compared with Newton-Raphson iterations, one clear advantage o f the tuned time-domain over iterative methods is the way the voltages change from step to step, even when the voltage collapses. Unlike Newton-Raphson iterations which just fail to converge for an ill-conditioned system, the tuned time-domain simulation w i l l show which voltages become unstable and collapse, as exemplified i n Figure 4.5. One way of assessing the voltage stability of a power system is by means o f a curve showing the 78 relation between the active power P and the bus voltage V . This P V curve shows the distance between the current operating point and the point at which the system collapses. Figure 4.5 suggests that it is possible to get the load P V curve profile o f a system node beyond the point o f voltage collapse, not only accurately but also economically. Figure 4.16 shows the load P V curve at node ' A ' for the left circuit o f Figure 4.4. The current source used to model the P Q load at node ' A ' was slowly incremented, maintaining the phase angle constant, which results also in a constant power factor, and the voltage and active power were recorded after each increment. That is, instead o f using equation (4.5) to compute the injected current at the P Q bus, the following formula was used: ' ' ' ' " . (3 = a - 6 = + (4,6) where Al is chosen to obtain a P V curve with the accuracy and detail required. Notice that by using this methodology there is no stable solution beyond the nose o f the P V curve and the simulated power system w i l l eventually collapse, but since the solution follows a trajectory towards the system instability point the states before reaching this point are displayed in the P V curve. For the P V curve o f Figure 4.16, Al was set to 0.005. Since the changes are very small from one solution point to the next, the steady state is reached quickly and with little error. The curve was obtained after 6608 time-step evaluations and consists o f 944 points using the tuned discretization formula obtained from a seventh order interpolating polynomial (equation (3.9)). Notice that the curve clearly shows the point o f voltage collapse and points thereafter. 79 1.2 P=10.000, V=0.987 lh~ 0.8 h 5? 0.6 > 0.4 h 0.2 h 0 0 50 100 150 PL(MW) Figure 4.16 Load P V curve at bus A for the circuit o f Figure 4.4. The solid line corresponds to the values obtained from the tuned time-domain simulation. The dashed line corresponds to a fifth order interpolating polynomial obtained from the solid line. The diamonds are the solutions obtained using Newton-Raphson iterations. For comparison purposes, some solutions obtained with Newton-Raphson iterations are also plotted in Figure 4.16 (the diamonds), showing a perfect match with the results obtained with tuned time-domain power flow. When the curve is very detailed, as that o f Figure 4.16, the point o f voltage collapse can be found directly from the tabular data resulting from the tuned time-domain simulation. A n alternative method consist in fitting the tabular data using a given order polynomial, derivate it, and find the tip point o f the P V curve. The dashed line of Figure 4.16 comes from a fifth order polynomial obtained with the Matlab function 'polyfit'. Due to the large number of points used, as well as the availability o f points around the tip of the P V curve, the fitting is excellent, resulting in a very accurate P V curve peak computation. 80 4.5.1 Effect of the Current Increment Size A very small increment value Al in equation (4.16) results in a very detailed and accurate load P V curve similar to the one shown in Figure 4.16. That kind o f accuracy is unlikely to be needed, therefore it would be wise to select an appropriate value o f Al such that the resulting P V curve is still accurate enough, but with fewer points, thus resulting i n faster computation times. Figure 4.17 shows the load P V curve at node ' A ' for the left circuit o f Figure 4.4, computed using three different values o f Al: 0.005, 0.05, and 0.5. The curves obtained with Al =0.005 and 0.05 are very similar, while the curve obtained with Al =0.5, which only has 10 points, is noticeably different. A s a general rule o f thumb, curves with 100 to 200 points seem to be accurate and detailed enough for engineering purposes. — 0.005 0.05 0.5 0.2 h 0-1 0 1 1 50 100 PL(MW) 1 150 Figure 4.17 Effect ofthe Al size in the generated P V curves. A small value o f Al results in a more accurate and detailed P V curve, but at a higher computational cost. The plot shows the curves generated with a Al o f 0.5, 0.05, and 0.005. 81 4.5.2 Comparison with Published PV Curves Figure 4.18 shows the load P V curve o f node N207 for the 32-bus system o f Figure 4.13. Dashed line (1) shows the P V curve as reported in [37] which was obtained using the 19 technique described in [38], while the solid line (2) shows the P V curve as obtained using the tuned time-domain technique with AI=0.1, resulting i n a curve with 88 points. Although similar in form, there is a visible difference on the resulting P V peaks, where the results o f the tuned time-domain curve are more conservative. The difference results from the way the curve reported in [37] is obtained. According to [38], the load for which the P V curve is required, is replaced by an equivalent admittance o f the form: P - iO s=- ~ylT , Y L L (4-17) s and the system is solved by what looks like standard Newton-Raphson iterations. The resulting admittance is then incremented by a small value and the system solved once again. The procedure is repeated a number o f times, until there are enough points to draw the curve. This technique results i n P V curves that do not keep the power factor constant, and therefore the difference with the results obtained using the tuned time-domain technique. Furthermore, the technique described in [38] is computationally more expensive than the tuned timedomain power flow technique describe in this report, since a number o f Newton-Raphson iteration solutions would have to be computed to obtain a detailed curve. This document is available in PDF formatfromhttp://thunderbox.uwaterloo.ca/~claudio/vswg/Chapter4.pdf. The figure in reference is 4-3.24(a). To obtain the figurefromthe document, it was printed to a file using a postscript driver. The resulting file was edited, and the table of values from which the figure is plot was exported to a Matlab filefromwhich the original figure was reconstructed. 19 82 1 2 1h \ \ \ N . . 0.8 h \ \ \ \ ] I 11 // '/ / > / / / / 0.4 h 0.2 h Ql 1 I I I I I 0 100 200 300 PL(MW) 400 500 600 Figure 4.18 Load P V curve for bus N207 o f the system o f Figure 4.13. The dashed line is the curve reported i n [37]. The solid line was obtained using the tuned time-domain power flow described in this report. 83 4.5.3 Computing the PV Curves for all System P Q B u s e s The technique described at the beginning o f section 4.5 to obtain the load P V curves can be applied to obtain all system load P V curves. To achieve this objective, the following changes were made to the tune time-domain program described in section 3.6: 1. A flag was added to each P Q load to indicate that the P V curve is required. If a load is sharing a bus with a P V generator, then the P V curve is not obtainable since the active power w i l l remain constant despite the load. 2. Once the static flow is completed, all the sources values are stored. This includes the independent P Q current P V voltage sources, as well as the historic current sources o f inductors, capacitors, R L , and R C series connections. 3. The P V curve is computed for each flagged load at a time. Before doing so, the static power flow steady state stored as in point 2 is restored. Once a number o f points for the P V curve are obtained, the resulting table is saved to disk. Experimental results show that the cost o f each P V curve generation, which requires the computation o f about 100 points, is approximately equal to twice the computational cost of the equivalent tuned time-domain static power flow method. For example, consider Figure 4.19 which shows all the P V curves for the system o f Figure 4.11. The static power flow for this system was obtained in 0.02 seconds while the eight P V curves were obtained 20 in 0.35 seconds. Similarly for the 32-bus system o f Figure 4.13, the static power flow was obtained in 0.04 seconds, and its 12 load P V curves were obtained i n 0.68 seconds. Finally, for the system o f Figure 4.14 and Figure 4.15, the static power flow was obtained in .0.08 seconds, and its 35 P V curves were obtained in 4.6 seconds. 2 0 Pentium Celleron 1 GHz, running Microsoft Windows XP. 84 C:/Eevs/steady/ieee14.net (BUS-4) C:/Eevs/steady/ieee14.net (BUS-5) P=47.800, V=1.018 P=7.600, V= 1.020 P=982.84C , V=0.655 j> fl QQ7 pi=u.yy / nf • nfftQ7Q pi—w.a/y P=s776:95f; V»0.593 + - „,--' 10 400 600 800 PL(MW) C:/Eevs/s1eady/ieee14.net (BUS-9) 1000 0 P=29.500, V=1.056 100 150 200 250 PL(MW) C:/Eevs/steady/ieee14.net (BUS-11) 600 800 PL(MW) C:/Eevs/steady/ieee14.net (BUS-10) nf pi— U.O*t 1 300 100 150 PL(MW) C:/Eevs/steady/ieee14.net (BUS-12) P-6.100, V-1;055 P=190.106, V=0.566 * " P<-°- 9 6 7 : P=i 84.698 V 0.57? X --' ' 50 100 150 PL(MW) C:/Eevs/steady/ieee14.net (BUS-13) 50 100 150 PL(MW) C:/Eevs/steady/ieee14.net (BUS-14) P=13.500, V=1.050 P=14.900, V=1.036 pl=0.9 50 P=171.096, V= 1571 * 50 P=3.500, V.1.057 _ 400 P=9.000, V.i.051 &0.6 > 50 200 pf=0.948 P=2f 2.393, V=3.563 X 100 150 PL(MW) 200 250 300 0 50 PL(MW) 100 Figure 4.19 A l l load P V curves for the 14-bus system o f Figure 4.11. 85 4.5.4 A Note on Load Voltage Stability Indices Voltage stability indices are scalar numbers used to asses the proximity to voltage collapse o f a power system. M a n y indices have been proposed as described in [37], but o f particular interest for this research is the P V curve based index proposed i n [38] which assesses the proximity o f the current operating point to the tip o f the curve: P P Lmgn- r — M -P A X "Q N p (A ' ] Q\ maxn where P n maX is the peak value o f the P V curve at node n and P no is the current operating active power at node n. Unfortunately, this index, as defined, is not taking into account the voltage collapse at different buses i n the system. This condition occurs i n some o f the nodes o f the system o f Figure 4.14 and Figure 4.15. Consider, for example, the P V curve at bus 38 in Figure 4.20, where a voltage collapse is detected first at bus 31. Equation (4.18) should be modified, therefore, to show this behaviour as follows: P -P p stablen Lmgn- r p no > (A 1 Q\ V+- > ly stablen where Pstablen is the last value o f the P V curve obtained before a voltage collapse in any o f the system buses. 86 ure 4.20 The load P V curve for bus 38 o f the I E E E 57-bus system showing the point at which the voltage at bus 31 collapses first (diamond). 87 CHAPTER 5: Dynamic Power flow Using Tuned TimeDomain Formulas 5.1 Introduction The selection o f the correct power system analysis tools for the study o f some particular system phenomena greatly affects the results obtained, as well as the time required to obtain such results. Traditional power flow methods are often used because o f their efficiency and simplicity, but they are usually inadequate for the simulation o f long time dynamic phenomena that often arises i n power systems [37;39]. O n the other hand, timedomain methods are better suited for this kind o f simulations since the transition from one time solution to the next represents the dynamics o f the simulated system. Furthermore, time domain methods are free o f numerical problems around the point o f system instabilities which allows them to closely follow the trajectory o f the system solution even for i l l conditioned power flow problems where traditional Newton-Raphson iterations fail to converge. Unfortunately, the numerical integration o f the differential-algebraic equations as well as the handling o f non-linear elements often posses a high computational cost that makes traditional time-domain techniques very slow. In the previous chapter, the tuned time-domain formulas were successfully used to solve the non-linear static power flow problem. Fortunately, the time-step required by this new technique is considerably larger than that o f traditional time-domain simulations (e.g. E M T P [30]), and due to the exactness the method achieves for the tuned-in frequency it constitutes a viable alternative for the solution o f the static power flow problem. Since the results are obtained using a time-domain simulation, it is possible to model and simulate the 88 dynamic behaviour o f the power system elements involved in the power flow problem, such as time and voltage dependencies. Before applying the tuned time-domain formulas to the solution o f the dynamic loadflow problem, the dynamic system time-scale must be well thought-out. That is, the time constants o f the system under study should be considered, so that the system, once discretized with the tuned time-domain formulas, does not have much introduced error in the simulation process. Consider for example the magnitude and frequency response o f the tuned formulas o f Figure 3.2, from which four main observations can be derived: 1) N o distortion is introduced by any o f the formulas at zero frequency ( D C ) . 2) N o distortion is introduced by any ofthe formulas at some tuned frequency. 3) The distortion introduced at low frequencies near D C is very small. Furthermore, it is not the magnitude error, which is almost flat, but the phase shift which introduces the most distortion. 4) The distortion introduced around the tuned frequency is also very small. Traditional time-domain simulations which use a fixed time-step, such as the E M T P , require a time step small enough to accurately represent the smallest time constant o f the simulated system. This time step depends on the discretization technique used, for example, for the trapezoidal rule it is often considered that a time step o f one tenth o f the maximum time constant is enough to simulate the circuit with a maximum distortion o f 3% [2]. This criterion can also be applied to the tuned difference equations described in this report. Consider for example Figure 3.2, from which it can be observed that for the formula derived from a seventh order interpolating polynomial the phase shift is less that 3% at about a 89 normalized frequency o f 0.025. I f a tuned frequency o f 60 H z is used, the resulting time step, using the values o f Table 3.1 and equation (3.8), would be approximately 5.71 ms, which corresponds to a smallest accurately discretized system time constant o f 0.025*5.71 ms=0.23 s , enough to represent many dynamic power system behaviours [37] [33]. In this chapter, simple dynamic models for the tuned time-domain dynamic power flow simulation are described. They include under-load tap-changing transformers ( U L T C ) , as well as time varying and voltage dependent P Q load. A l s o , the method used to handle outages in the tuned time-domain simulation program is described. Finally a system example is presented, which uses the developed models and techniques, with its corresponding time domain solution. 5.2 Dynamic Modelling of ULTC Transformers Under-load tap-changing transformers ( U L T C ) are used to control the amount o f active and reactive power delivered to the loads as well as to regulate the voltage at the buses they are connected to [33] [26]. U L T C transformers are usually made with 8, 16, or 32 steps, with a regulation range o f ± 1 0 % o f the rated tap value [40], but it is not uncommon to find transformers with a range o f regulation o f ± 7 . 5 % and ± 5 % . Figure 5.1 shows the basic connection o f an U L T C transformer. The tap changing contacts are placed near the center o f the transformer to keep the coil physically and electrically symmetric i n order to avoid electrical and mechanical stress. In addition to the If a formula derivedfroma fourth order polynomial were to be used, the smaller accurately discretized time constant would be 0.15*2.55 ms=382.5 |as! Unfortunately this formula is less stable than the one derived from a seventh order polynomial as shown in Figure 3.3. 21 90 step change size, the desired voltage range, the number o f taps, the delay for the first tap movement and delay time o f any tap change thereafter is usually specified [33]. The inclusion o f the dynamic model o f U L T C transformers in the tuned time-domain program o f section 3.6 is similar to that described in section 4.2.2. The equivalent circuit o f Figure 4.2 is adjusted according to the selected tap value, and is then included in the 7 matrix of the system o f equations (1.1). The Y matrix has to be L U decomposed every time the U L T C transformer position changes, before solving the system for the node voltages and voltage source currents. Figure 5.1 Basic connexion o f a U L T C transformer. Unlike standard Newton-Raphson iterations which produce a real number for the tap position that may or may not be one o f the actual available physical values, the tuned timedomain method described in this report selects the closest tap position available to satisfy the transformer operation conditions. For example, consider Figure 5.2 which shows all the buses magnitudes and phases after solving the static power flow o f the system o f Figure 4.11. The tap changing transformers are required to maintain the voltage at buses 4 and 5 to 1.02 and 1.027 V (PU) respectively. In Figure 5.2 the action o f the tap transformers is clearly 91 visible as a step change on the system buses voltages magnitude and phase, until the desired bus voltages are obtained. When the load P V curves are obtained using the model o f the U L T C ' s instead o f fixed-tap transformers some interesting results are obtained. For instance, Figure 5.3 shows the P V curves for nodes 4 and 5 o f the system o f Figure 4.11, both with fixed-tap and U L T C transformers. Although the initial effect o f the U L T C is to increase the bus voltage as the load at such bus is increased, shortly after the effect reverses, resulting in a less stable system, as shown by the peak values o f the P V curves in Figure 5.3. This undesirable effect on the stability o f the system caused by the U L T C is also commented on by some authors, for 22 example Kundur in [33] , who made the observation for a similar case. Page 986: "The net effect of each tap movement of transformer T6 is to reduce bus 11 voltage rather than increase it." 92 1.2 1.15 Time (s) Figure 5.2 Effect of the U L T C transformers on the static power flow solution of the system o f Figure 4.11. 93 BUS-4 1.2 I --- Fixed-tap ULTC 0.8 \ \ \ \ 1 £ 0.6 / j / 0.4h 0.2h 200 400 600 800 1000 PL(MW) BUS-5 1.2 Fixed-tap ULTC --- 0.8 \ \ \ 0,0.6 > \ / / I y 0.4 0.2 0 0 100 200 300 400 PL(MW) 500 600 700 800 Figure 5.3 Effect o f the U L T C transformers in the P V curves o f buses 4 and 5 for the system of Figure 4 . 1 1 . 94 5.3 Dynamic Modelling of Loads The simplest way o f modelling the dynamic behaviour o f loads consists in the adjustment o f its complex power value over time: S(t) = P(t) + Q(t)j. (5.1) The implementation o f equation (5.1) on the tuned time-domain power flow program described i n section 3.6 is straight forward, since the P Q loads are modelled as current sources. A l s o , since any changes on the current sources are only reflected in the vector o f independent sources o f the system o f equations (1.1), the computational cost associated with such changes is very small. Furthermore, the variation o f loads with voltage as described i n [41] can be easily modelled as well. Assuming the system frequency to be constant, the load variation with voltage can be represented using an exponential model, p=p (yr 0 v - ' (5 2) where Po and Qo are the nominal power values at nominal voltage Vo, and pv and qv are the coefficients selected according to the load type as shown in Table 5.1. Load Type pv qv Incandescent lights Fluorescent lights Heaters Induction Motors, full load Induction Motors, half load 1.6 1.2 2.0 0.1 0.2 0 3.0 0 0.6 1.6 Table 5.1 Some typical coefficients for the exponential load model o f equation (5.2). 95 The effect o f a dynamic load model i n a simulation can be illustrated when the load connected at bus 5 i n the system o f Figure 4.11 is step-incremented by the values given i n Table 5.2. After 10 seconds (when the steady state has been reached), the value o f the load is increased every five seconds keeping the power factor constant and with the U L T C transformers activated so as to compensate for the voltage drop at the controlled buses. Figure 5.4 shows the voltages magnitude and phase for buses 4 and 5. Figure 5.4 clearly shows the effect o f the load variation on the buses voltages, as well as the compensating effect o f the U L T C transformers. Notice that the final load value is near the point o f voltage collapse, as indicated by the nose peak o f Figure 5.3. Time (s) P Q 10.0 2.5 0.5263 15.0 3.5 0.7368 20.0 4.5 0.9474 25.0 5.5 1.1579 30.0 6.5 1.3684 35.0 7.5 1.5789 ' Table 5.2 Step values over time for the load connected at bus 5 in the system o f Figure 4.11. The load can be dynamically adjusted using any function. One function that is found in the literature is a lineal change over time o f the load values [37;38]. Figure 5.5 shows the voltages magnitude and phase at buses 4 and 5 for the system o f Figure 4.11 when the load at bus 5 is incremented linearly using the equation P m - P final initial ^final P(t) = m(t-t ) inMal where P/j«af=7.5, Wia/=10.0, tf i=35.0, ina (5.3) ^initial + P initial and Pmrnai is the resulting value o f the static power flow solution. Similarly, Q(t) was set to keep the power factor constant. 96 1.05 BUS-5 BUS-4 0.95 0.9 \ a 0.85 0.8 0.75 0.7 0.65 10 20 30 40 50 Time (s) 60 10 BUS-5 BUS-4 -10 -20 -30 -40 -50 10 20 30 40 50 60 Time (s) Figure 5.4 Magnitude and phase for the voltages at buses 4 and 5 for the system of Figure 4.11 when the load at bus 5 changes with the values given in Table 5.2. 97 1.05 30 Time (s) Figure 5.5 Magnitude and phase for the voltages at buses 4 and 5 for the system of Figure 4.11 when the load at bus 5 changes linearly with time. 98 5.4 Simulating Outages Outages such as the loss o f a transmission line or generator can be simulated using a simple switch model. Switches can be modelled using a resistance whose value is very high when the switch is open and very low when the switch is closed. In any case, the matrix o f admittances Y o f the system o f equations (1.1) changes whenever the state o f the switch changes, which results in a new L U decomposition o f the corresponding system o f equations, and the costs associated with it (section 4.4). Resistor values o f 1 0 Q for an opened switch and 1 0 " Q for a closed switch, which 1 0 10 are a few orders o f magnitude different from those normally encountered in power flow cases, seem to work satisfactorily. The operation o f the switch is considered instantaneous, i n the sense that the current across the switch is not expected to cross zero before toggling its state. This behaviour differs from the E M T P implementation o f switches, since i n the tuned timedomain simulation the time steps are much larger, and the point o f zero current crossing probably w i l l be far from any o f the discrete time points anyhow. Figure 5.6 shows the voltages magnitude and phase for buses 4 and 5 o f the system o f Figure 4.11. The fixed load at bus 5 was incremented to 0.76+0.16j so to emphasize the effect o f the outage: after 40 seconds the lines connecting buses 4 and 5 were removed by means o f two switches at both ends o f the line. The resulting transient is clearly visible, as well as the effect o f the U L T C transformers that managed to bring the buses voltage back to their desired values just after a few tap changes. 99 1.05 r~ l BUS-5 BUS-4 10 0 20 30 40 50 60 Time (s) I BUS-5 BUS-4 a) - 5 I \ 0 10 20 30 40 50 60 Time (s) Figure 5.6 Magnitude and phase for the voltages at buses 4 and 5 for the system o f Figure 4.11 when the line connecting both buses is removed. 100 5.5 Example of Dynamic Power flow Solution Linearly changing loads, U L T C transformers, and outage simulation mechanism are required to simulate the dynamic power flow behaviour o f the O G R E 32-bus system o f Figure 4.13, as presented in [37]. According to [37], the C I G R E 32-bus system was originally created to study the voltage collapse o f the northern part o f Belgium in 1982 and reported i n [42]. The series o f events that drove the system to collapse are described as follows: 1. A t t=30 s, the loads at buses N201 through N207 started to increase steadily at a rate o f 30% in 7200 s. A l l other loads remain constant. Although the U L T C transformers in the actual system were controlled by the operator, the tap o f the transformers connected to loads N201 to N207 was allowed to change automatically in the simulation. 2. The transmission line connecting buses N 3 and N 1 6 was tripped off at t=5000 s. 3. A t t=7230 s the loads at buses N201 to N 2 0 7 stop increasing. 4. A t 7400 s, generator M 2 trips off. For the tuned time-domain dynamic power flow simulation, the U L T C transformers are assumed to have 32 steps with a range o f 10%, with a first tap movement delay o f 30 s and 30 s thereafter. It is also assumed that the linearly increasing P Q loads keep their power factor constant as described in section 5.3. For the outage modelling o f the line between nodes N 3 and N 1 6 , switches were placed at both ends o f the line which trip off at the required outage time o f 5000 s. The tuned time-domain dynamic power flow simulation was run for a total simulation time o f 7500 seconds. For a time step o f approximately 5.71 ms, corresponding to the tuned 101 formula derived from a seventh order polynomial, the number o f time steps required to complete the simulation is around one million. In order to save disk space and speed-up plotting, the output o f the simulation was saved every 700th time step, or around every fourth second. The simulation was completed in about 68 seconds i n a 1 G H z Pentium Celeron running Windows X P . Figure 5.7 shows the voltages magnitude and phase for buses N 2 0 1 , N202, and N 2 0 4 that result from the tuned time-domain simulation. A s the loads increase over time, the U L T C transformers adjust their tap position so as to keep the voltage at the load buses more or less constant. A t t=5000 s the line connecting busses N 3 and N 1 6 trips off, which results in the large transients shown in the plots. The system recovers from this outage, but the U L T C transformers quickly reach their tap limit at t=5165 s. From that point on, the buses voltage magnitude decrease steadily until the voltage quickly collapses at t=7160 s when the buses voltage become erratic. The results o f Figure 5.7 are very similar to those presented in [37] but the predicted time for the collapse is different since the aforementioned reference states that the voltage collapse is due to the tripping o f generator M 2 at t=7400 s. This is caused by the way loads were represented i n reference [37]; "constant impedances with internal taps (no limits) to avoid power flow convergence problems", while i n the tuned time-domain dynamic power flow simulation they were represented as pure P Q loads. To obtain better results, the loads can be represented as suggested in [42], the original source for the C I G R E 32-bus case, as a mix o f fixed and voltage dependent loads given by Page 4-46, Figure 4.3-25. 102 P = P (0.6 + 0.4V ) 2 0 (5-4) Q = Q (0.5 + 0.5V ) 2 0 where the quadratic part results from equation (5.2), with pv=qv=2 and V =l. 0 Figure 5.8 shows the magnitude and phase for the voltages at buses N 2 0 1 , N202, and N204 when using equation (5.4). Since the loads do not increase as fast as before due to their voltage dependence, the system remains stable after the loads stop increasing at t=7230 s. It is now the loss o f the generator at bus M l at t=7400 s which precipitates the voltage collapse as shown in the figure. 103 i TWyy <luulJ- uu «UI]l) u UUU|J U ""UUIJU ' i °UU1J i '•. M n o i Uu ul|- "Hull u K j 0 1000 2000 3000 4000 . ., ul 5000 : N202 \ . [ N204 6000 7000 Time (s) -5 I -10 N201 N202 -15 N204 -20 -25 0 1000 2000 1 1 3000 4000 5000 6000 7000 Time (s) Figure 5.7 Dynamic voltage power flow solution at buses N 2 0 1 , N202, and N204 for the O G R E 32-bus system o f Figure 4.13 using pure P Q loads. 104 1.1 1.08 1.06 1.04 N201 N202 1.02 > N204 CD | 1 Q. E < 0.98 0.96 0.94 0.92 h 0.9 1000 2000 3000 4000 Time (s) 5000 6000 7000 -10 N201 N202 -15 N204 to CO .c Q. -20 -25 -30 1000 2000 3000 4000 Time (s) 5000 6000 7000 ure 5.8 Dynamic voltage power flow solution at buses N 2 0 1 , N202, and N204 for the O G R E 32-bus system o f Figure 4.13 using voltage dependant loads. 105 CHAPTER 6: Conclusions and Recommendations for Future Work. 6.1 Summary and Conclusions In Chapter 1 the description o f some o f the basic techniques used in the simulation o f electric power circuits was presented. The use o f nodal analysis and M N A for the formulation o f the circuit equations was presented, and the use o f difference equations to convert the resulting system o f differential equations into a set o f algebraic equations was described. A l s o in Chapter 1 important properties o f the difference equations, such as A - stability and distortion were presented. Finally, phasor analysis for the solution o f linear A C circuits and Newton-Raphson iterations for the solution o f non-linear power flow problems were described. The development o f adjustable difference equations, suitable for electric circuit simulators, was presented in Chapter 2. These new difference equations were derived from interpolating polynomials o f a given order, and their coefficients were written as a function of an adjusting parameter k. B y setting this parameter k to zero, the resulting difference equation becomes the B D F o f the same order o f the used interpolating polynomial. O n the other hand, i f the adjusting parameter is set to one, the well known Backward Euler formula is obtained. B y setting the parameter k to values between zero and one, the adjustable difference equations can be attuned to meet some stability and accuracy criteria. example, the adjustable For difference equation resulting from a third order interpolating polynomial can be made A-stable with a A: o f about 0.14, and the resulting formula w i l l 106 exhibit accuracy and stability properties very similar to those o f the well known trapezoidal formula. When interpolating polynomials o f order higher than four are used, the resulting adjustable difference equations can be tuned to a particular frequency, so that there is no distortion at such frequency. This results i n difference equations that produce exact results for linear circuits working at one particular frequency, such as for example those arising i n electric power circuits. Chapter 3 describes the derivation o f these tuned formulas. The resulting tuned formulas use fixed time steps, considerably larger than those required by traditional formulas such as the Trapezoidal and Backward Euler rules for the particular case of simulating o f one frequency circuits and result in faster and more accurate simulations. Additionally, these tuned formulas show very good stability properties, in particular the tuned formula derived from a seventh order interpolating polynomial, which is almost A stable. A l s o i n Chapter 3, linear electric power circuits are simulated using the tuned formulas in the time-domain, and the results compared with the solution o f phasor analysis, where a perfect match is obtained. The computational cost o f the tuned time-domain simulation is compared against phasor analysis, and it was concluded that the linear tuned time-domain simulation w i l l result in faster solutions when the number o f nodes i n the circuit being simulated is larger than the number o f steps required for the time-domain transients to damp. In Chapter 4 the solution o f the static power flow problem using the tuned timedomain methods developed in Chapter 3 was achieved. Before doing so, suitable models for elements commonly found in static power flow cases were developed. 107 These include transmission lines, fixed tap transformers, shunts, and in particular P Q loads and P V sources. P Q loads were modelled using a sinusoidal current source whose magnitude and phase constantly adapt so as to match the required complex power load. Similarly, P V sources were modelled with a sinusoidal voltage source for which the magnitude is known, and the phase is adjusted using a simple PI controller to match the required active power. A l s o i n Chapter 4, three standard static power flow problems were solved using the tuned time-domain simulator. These include the I E E E 14-bus, the O G R E 32-bus, and the I E E E 57-bus cases, where the results match accurately those obtained using standard Newton-Raphson iterations. The computational cost o f the tuned time-domain static power flow is then compared to standard Newton-Raphson iterations. It was determined that the tuned time-domain static power flow would be faster reaching a solution than the standard Newton-Raphson iterations for systems with many more P Q buses than P V buses. Finally, in Chapter 4, the stability properties o f the tuned time-domain static power flow, which is capable o f producing usable results on and around the point o f voltage collapse, is used to generate P V curves. A n adjustable current source is used to represent a varying P Q load with constant power factor. The resulting voltage, for a given active power, is then obtained from which a load P V curve can be plotted. The resulting P V curves are validated against the results obtained using Newton-Raphson iterations and compared with published load P V curves where the results show good agreement. In Chapter 5 the dynamic power flow solution using the tuned time-domain formulas is obtained. This is possible because the tuned time-domain formulas are not only exact for the tuned frequency but also for low frequencies near D C . Models for the ULTC transformers, time varying loads and switches used to simulate outages were developed. 108 Finally the dynamic power flow simulation o f the O G R E 32-bus system is presented and the results compared with published data, where once again good agreement is obtained with the results obtained using the tuned time-domain method. 6.2 Recommendations for Future Work The results obtained in this research, due to their fundamental nature, in particular the tuned difference equations for the solution o f differential equations, may have many applications On the solution o f many other oscillatory problems found in nature. Concentrating into the scope o f this report, this research can be further expanded by combining it with the research o f other members and former members o f the power group at UBC: 1. Marti in [43;44] extended the concepts o f K r o n ' s diakoptics [45] to multi-area Thevenin equivalents ( M A T E ) i n which each subsystem can be solved with a different solution technique while maintaining a simultaneous solution o f the entire network. Linares i n [46] developed highly efficient computer algorithms to implement M A T E in real-time simulations. Hollman in [47] mapped M A T E algorithms to a PC-cluster architecture with multiple processing nodes. Using M A T E ' S methodology in which a large circuit is torn into smaller more manageable circuits, the solution o f static and dynamic power flow problems using the tuned time-domain method can be greatly speed up. 2. Louie in [48] researched the modeling o f electric power circuit loads. His research can be used to greatly improve the results obtained with the tuned time-domain method by replacing the simple load models presented in this research with more accurate models. 109 3. Moreira also in [44], studied the simulation o f electric circuits using latency. In his research an electric circuit is split into sections with remarkable different time constants, where each section is solved independently with a different time-step and latter combined to obtain the correct overall solution. B y using this solution methodology the tuned time-domain method can be used to simulate the slow changing part o f the circuit, while E M T P - l i k e methods would be used to solve the fast changing part o f the circuit. 4. Finally, L u k i c in [31] has demonstrated the feasibility o f using time-domain simulations for the solution o f many traditional power flow problems. developed tuned formulas are also directly applicable to her research. 110 The Bibliography [1] Dommel, H . W . , "Digital computer solution o f electromagnetic transients in single and multiphase networks," IEEE Transactions on Power Apparatus and Systems, v o l . pas-88, no. 4, pp. 388-399, 1969. [2] Marti, J. R. and L i n , J., "Suppression o f numerical oscillations in the E M T P , " IEEE Transactions on Power Systems., vol. 4 pp. 739-745, June 1989. [3] Ascher, U . M . and Petzold, L . R., Computer methods for ordinary differential equations and differential-algebraic equations Philadelphia: Society for Industrial and Applied Mathematics, 1998. [4] Curttis, C . F. and Hirschfelder, J. O., "Integration o f stiff equations," U.S. National Academy of Science, v o l . 38 pp. 235-243, 1952. Proceedings [5] Dahlquist G . G , " A special stability problem for linear multistep methods," BIT, v o l . 3 pp. 27-43, 1963. [6] Gear, C . W . , Numerical initial value problems in ordinary differential Englewood Cliffs, N . J . : Prentice-Hall, 1971. [7] Wildund, O., " A note on unconditionally stable linear multistep methods," BIT, v o l . 7 pp. 65-70, 1967. [8] Spijker, M . N . , "Stiffness in numerical initial-value problems," Journal Of Computational And Applied Mathematics, v o l . 72, no. 2, pp. 393-406, A u g . 1996. equations [9] Verwer, J. G . , B l o m , J. G . , V a n Loon, M . , and Spee, E . J., " A comparison o f stiff ode solvers for atmospheric chemistry problems," Atmospheric Environment, v o l . 30, no. l , p p . 49-58, Jan. 1996. [10] Hairer, E . and Wanner, G . , "On the instability ofthe B D F formulas," SIAM on Numerical Analysis., no. 20., pp. 1206-1209, 1983. [11] Fornberg, B . , "Calculation o f Weights i n Finite Difference Formulas," SIAM Review, vol. 40, no. 3, pp. 685-691, Sept. 1998. [12] Fredebeul, C , " A - B D F : A Generalization o f the Backward Differentiation Formulae," SIAM Journal on Numerical Analysis, v o l . 35, no. 5, pp. 1917-1938, Oct. 1999. [13] Hosseini, S. M . and Hojjati, G . , "Matrix free M E B D F method for the solution o f stiff systems o f O D E s , " Mathematical And Computer Modeling, v o l . 29, no. 4, pp. 67-77, Feb. 1999. Ill Journal [14] W u , X . Y . , " A sixth-order A-stable explicit one-step method for stiff systems," Computers & Mathematics With Applications, v o l . 35, no. 9, pp. 59-64, M a y 1998. [15] Bettis, D . G . , "Stabilization o f finite difference methods o f numerical integration," Celestial Mechanics, v o l . 2, no. 3, pp. 282, 1970. [16] Shampine, L . F . and Reichelt, M . W . , "The M A T L A B O D E suite," SIAM Journal on Scientific Computing, v o l . 18, no. 1, pp. 1, 1997. [17] Trnka, O., Hartman, M . , and Svoboda, K . , " A n alternative semi-implicit Euler method for the integration o f highly stiff nonlinear differential equations," Computers & Chemical Engineering, v o l . 21, no. 3, pp. 277, 1997. [18] Cremascoli, P. and Gubian, P. Mathematical behavior o f integration methods with non continuous first derivative devices i n circuit simulation. 2, 1225-1228. 1993. Proceedings o f 36th Midwest Symposium on Circuits and Systems. [19] Gubian, P. and Zanella, M . Stability properties o f integration methods i n S P I C E transient analysis. 2701-2704. 1991. I E E E International Symposium on Circuits and Systems. [20] Alvarado, F. L . , Lasseter, R . H . , and Sanchez, J. J., "Testing o f trapezoidal integration with damping for the solution o f power transient problems," IEEE Transactions on Power Apparatus and Systems., no. 102., pp. 3783-3790, Dec. 1983. [21] Smith, J. M . , Mathematical modeling and digital simulation for engineers and scientists N e w York: Wiley, 1977. [22] Gear, C . W . and Osterby, O., "Solving ordinary differential equations with discontinuities," ACM Transactions on Mathematical Software., v o l . 10, no. 1, pp. 23-44, Mar. 1984. [23] Calvino-Fraga, J. Implementation o f a Real Time Simulator for Protective Relay Testing. 6-15-1999. The University o f British Columbia. Master Thesis Report. [24] Aprille, T. J., Jr. and Trick, T. N . , "Steady-state analysis o f nonlinear circuits with periodic inputs," Proceedings of the IEEE, v o l . 60, no. 1, pp. 108, 1972. [25] Perkins, B . K . , Marti, J. R., and Dommel, H . W . , "Nonlinear elements i n the E M T P : steady-state initialization," IEEE Transactions on Power Systems, v o l . 10, no. 2, pp. 593, 1995. [26] Bergen, A . R. and Vittal, V . , Power systems analysis, 2nd ed. ed. Upper Saddle River, N J : Prentice H a l l , 2000. [27] V l a c h , J. and Singhal, K . , Computer methods for circuit analysis and design N e w York: V a n Nostrand Reinhold, 1983. 112 [28] Chung-Wen, H . , Ruehli, E . , and Brennan, P. A . , "The modified nodal approach to network analysis," IEEE Transactions on Circuits and Systems, v o l . C A S - 2 2 , no. 6, pp. 504-509, June 1975. [29] Chua, L . O. and L i n , P., Computer Aided Analysis of Electronic circuits: Algorithms & computational techniques. Englewood Cliffs, N . J . : Prentice-Hall, Inc., 1975. [30] Dommel, H . W . , EMTP Theory Book, Second ed. Vancouver, B C : Microtran Power System Analysis Corporation, 1996. [31] L u k i c , M . , Feasibility of security assessment in power systems using full time-domain solutions in the EMTP Vancouver: University o f British Columbia. Master Thesis Report, 2000. [32] Burden, R. L . , Reynolds, A . C , and Faires, J. D . , Numerical analysis, 2d ed. Boston, Mass.: Prindle, Weber & Schmidt, 1981. [33] Kundur, P., Power system stability and control N e w Y o r k : M c G r a w - H i l l , 1994. [34] "Common format for exchange o f solved power flow data," IEEE Transactions Power Apparatus and Systems, vol. P A S - 9 2 , no. 6, pp. 1916-1925, N o v . 2000. [35] Grainger, J. J. and Stevenson, W . D . , Power system analysis N e w York: M c G r a w H i l l , 1994. [36] Freris, L . L . and Sasson, A . M . , "Investigation o f the load-flow problem," Proc.IEE, vol. 115, no. 10, pp. 1459-1470, A u g . 1968. [37] Voltage Stability Assessment: Concepts, Practices and Tools. Canizares, C . Editor. I E E E / P E S Power System Stability Subcommittee Special Publication. 2002. on [38] Nagao, T., Tanaka, K . , and Takenaka, K . , "Development o f static and simulation, programs for voltage stability studies o f bulk power system," IEEE Transactions on Power Systems, v o l . 12, no. 1, pp. 273-281, Feb. 1997. [39] V a n Cutsem, T., Jacquemart, Y . , Marquet, J. N . , and Pruvot, P., " A comprehensive analysis o f mid-term voltage stability," IEEE Transactions on Power Systems, v o l . 10, no. 3, pp. 1173-1182, A u g . 1995. [40] Westinghouse, E . C , Electrical transmission and distribution reference book, [4th completely rewritten ed.] East Pittsburgh, 1950. [41] Berg, G . J., "Power-system load representation," Proceedings of the Institution Electrical Engineers, v o l . 120, no. 3, pp. 344-348, Mar. 1973. [42] Indices Predicting Voltage Collapse Including Dynamic Phenomena. Hatziargyriou, N . D . and V a n Cutsem, T. technical report T F 38-02-11, C I G R E . 1994. 113 of [43] Marti, J. R. Multi-area thevenin equivalent. A circuit concept. Internal memo to H . W . D o m m e l and L.R.Linares.Department o f Electrical Engineering, The University of British Columbia. 1993. [44] Marti, J. R., Linares, L . R., Hollman, J. A . , and Moreira, F. A . O V N I : Integrated software/hardware solution for real-time simulation o f large power systems. Conference Proceedings o f the 14th power system computation conference P S C C 0 2 . 2002. 6-24-2002. [45] K r o n , G . , "Tensorial analysis o f integrated transmission systems, part III - the primitive division," AIEE Transactions, vol. 71 pp. 814-821, 1952. [46] Linares-Rojas, L . R., OVNI (Object Virtual Network Integrator) Vancouver: University o f British Columbia. Ph. D . Thesis Report, 2000. [47] Marti, J. R., Hollman, J. A . , and Calvino-Fraga, J. Implementation o f a real-time distributed network simulator with PC-cluster. 223-227. 2000. Trois-Rivieres, Que., Canada, I E E E Comput. Soc. Proceedings International Conference on Parallel Computing i n Electrical Engineering. P A R E L E C 2000. [48] Louie, K . W., Aggregation of voltage and frequency dependent electrical loads Vancouver: University o f British Columbia. Ph. D . Thesis Report, 1999. 114
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Tunable difference equations for the time-domain simulation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Tunable difference equations for the time-domain simulation of power system operating states Calviño-Fraga, Jesús 2003
pdf
Page Metadata
Item Metadata
Title | Tunable difference equations for the time-domain simulation of power system operating states |
Creator |
Calviño-Fraga, Jesús |
Date Issued | 2003 |
Description | This report presents the development of a new family of difference equations for the numerical solution of systems of differential equations arising in electric power circuits. The main advantage of these difference equations is their adjustability, which allows for the tuning of the resulting formulas while using large time steps. The tuned formulas have been used successfully in the time-domain simulation of 50/60-Hz solutions of electric power circuits for linear and non-linear problems. The linear steady state and the non-linear power flow problems can be solved, with no error either in magnitude or in phase, at a comparable computational cost to traditional techniques such as the well-known phasor analysis and Newton-Raphson iterations. Additionally, the techniques developed do not use complex numbers, resulting in a faster matrix factorization for linear systems, or a smaller Jacobian matrix for non-linear systems. The main advantage of the proposed techniques is that due to the time-domain nature of the simulation, where the transition from one time solution to the next represents a small transition between states, the technique is able to closely follow the trajectory of the system even for ill-conditioned power flow problems where traditional Newton-Raphson iterations fail to converge, such as when the system buses are beyond or near the point of voltage collapse. With the proposed novel formulas, the computational cost of generating PV curves across all busses of the power system network is greatly reduced, and the resulting curves clearly show the point of voltage collapse as well as any unstable voltage/power combination thereafter. The ability to quickly obtain PV curves makes the proposed techniques particularly well suited for online power system monitoring of operating states. Also, since these new difference equations introduce very little error for frequencies lower than the tuned frequency, they are also well suited for the solution of dynamic power flow problems such as those resulting from loads variation, transformers tap change, and outages. |
Extent | 4700737 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-17 |
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.0091165 |
URI | http://hdl.handle.net/2429/15027 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-860019.pdf [ 4.48MB ]
- Metadata
- JSON: 831-1.0091165.json
- JSON-LD: 831-1.0091165-ld.json
- RDF/XML (Pretty): 831-1.0091165-rdf.xml
- RDF/JSON: 831-1.0091165-rdf.json
- Turtle: 831-1.0091165-turtle.txt
- N-Triples: 831-1.0091165-rdf-ntriples.txt
- Original Record: 831-1.0091165-source.json
- Full Text
- 831-1.0091165-fulltext.txt
- Citation
- 831-1.0091165.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-0091165/manifest