UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Transient stability analysis of multimachine power systems by catastrophe theory Mihirig, Ali Mohamed 1987

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1988_A1 M53.pdf [ 4.8MB ]
Metadata
JSON: 831-1.0098151.json
JSON-LD: 831-1.0098151-ld.json
RDF/XML (Pretty): 831-1.0098151-rdf.xml
RDF/JSON: 831-1.0098151-rdf.json
Turtle: 831-1.0098151-turtle.txt
N-Triples: 831-1.0098151-rdf-ntriples.txt
Original Record: 831-1.0098151-source.json
Full Text
831-1.0098151-fulltext.txt
Citation
831-1.0098151.ris

Full Text

TRANSIENT STABILITY ANALYSIS OF MULTIMACHINE POWER SYSTEMS BY CATASTROPHE THEORY by A l l MOHAMED MIHIRIG B.Sc, University of A l Fateh, 1978 M.A.Sc, University of B r i t i s h Columbia, 1984 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1987 ©Ali Mohamed M i h i r i g , 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or .by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date /4d\rcl\ <S y DE-6G/81) ABSTRACT Transient s t a b i l i t y analysis i s an Important part of power system planning and operation. For large power systems, such analysis i s very demanding i n computation time. On-line transient s t a b i l i t y assessment w i l l be necessary for secure and r e l i a b l e operation of power systems i n the near future because systems are operated close to t h e i r maximum l i m i t s . In the l a s t two decades, a vast amount of research work has been done i n the area of fast transient s t a b i l i t y assessment by d i r e c t methods. The major d i f f i c u l t i e s associated with d i r e c t methods are the l i m i t a t i o n s i n the power system model, determination of transient s t a b i l i t y regions and adaptation to changes i n operating conditions. In t h i s thesis catastrophe theory i s used to determine the transient s t a b i l i t y regions. Taylor series expansion i s used to f i n d the energy balance equation i n terms of c l e a r i n g time and system transient parameters. The energy function i s then put i n the form of a catastrophe manifold from which the b i f u r c a t i o n set i s extracted. The b i f u r c a t i o n set represents the transient s t a b i l i t y region i n terms of the power system transient parameters bounded by the transient s t a b i l i t y l i m i t s . The transient s t a b i l i t y regions determined are v a l i d for any changes i n loading conditions and f a u l t l o c a t i o n . The transient s t a b i l i t y problem i s dealt with i n the two dimensions of transient s t a b i l i t y l i m i t s and c r i t i c a l c l e a r i n g times. Transient s t a b i l i t y l i m i t s are given by the b i f u r c a t i o n set and the c r i t i c a l c l e a r i n g times are calculated from the catastrophe manifold equation. The method achieves a breakthrough i n the modelling problem because the e f f e c t s of ex c i t e r response, f l u x decay and - i i -systems damping can a l l be included i n the transient s t a b i l i t y analysis. Numerical examples of one-machine i n f i n i t e - b u s and multi-machine power systems show very good agreement with the time solution i n the p r a c t i c a l range of f i r s t swing s t a b i l i t y a n a l y sis. The method presented f u l f i l l s a l l requirements for on-line assessment of transient s t a b i l i t y of power systems. - i i i -TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS i v LIST OF TABLES v i i LIST OF FIGURES v i i i ACKNOWLEDGEMENT x i CHAPTER 1 INTRODUCTION 1 CHAPTER 2 REVIEW OF THE STABILITY PROBLEM AND ITS SOLUTION 5 2.1 Introduction 5 2.2 The S t a b i l i t y Problem 5 2.3 Basic Power System S t a b i l i t y Model 7 2.3.1 Synchronous Generator 8 2.3.2 E x c i t e r and Governor Control Systems 10 2.3.3 Transmission System and Loads 11 2.4 Solution of the Transient S t a b i l i t y Problem 12 2.4.1 Numerical Integration 12 2.4.2 Lyapunov's Direct Method 15 2.4.3 Pattern Recognition Method 18 2.5 Discussion of E x i s t i n g Methods 22 CHAPTER 3 APPLICATION OF CATASTROPHE THEORY TO TRANSIENT STABILITY ANALYSIS 26 3.1 Introduction ' 26 - i v -Page 3.2 The Catastrophe Theory 27 3.2.1 Catastrophe Theory and B i f u r c a t i o n Analysis 28 3.3 Ap p l i c a t i o n to the Steady State S t a b i l i t y Problem 29 3.3.1 C y l i n d r i c a l - R o t o r I n f i n i t e Bus Power System 31 3.3.2 Numerical Results 38 3.4 A p p l i c a t i o n to the Transient S t a b i l i t y Problem 38 3.4.1 Single Machine Infinite-bus Power System. 41 3.4.2 Results 48 3.4.3 Equivalent Two Machine Power System 52 3.5 Summary 54 CHAPTER 4 TRANSIENT STABILITY REGIONS OF MULTI-MACHINE POWER SYSTEMS 55 4.1 Introduction 55 4.2 Dynamic Equivalent of the C r i t i c a l Machines 56 4.2.1 The Transient S t a b i l i t y Regions 61 4.3 Taylor Expansion of the Accelerating Power 65 4.3.1 Transient S t a b i l i t y Regions Using Taylor Expansion 70 4.4 I d e n t i f i c a t i o n of the C r i t i c a l Machines 72 4.5 Numerical Examples 74 4.5.1 The Three-Machine System 75 4.5.2 CIGRE 7-Machine Test System 84 - v -Page 4.6 Discussion of Results 84 CHAPTER 5 INCLUSION OF DAMPING, FLUX DECAY AND EXCITATION RESPONSE 90 5.1 L i m i t a t i o n of the C l a s s i c a l Model 90 5.2 Damping 91 5.3 E x c i t a t i o n Response and Flux Decay 95 CHAPTER 6 CONCLUSION 107 REFERENCES 112 APPENDIX A 117 APPENDIX B 122 APPENDIX C 126 APPENDIX D 130 APPENDIX E 131 - v i -LIST OF TABLES Page TABLE 3.1 The 7-elementary catastrophes 30 TABLE 4.1 Network parameters of the three-machine power system. 75 TABLE 4.2 Generator data and i n i t i a l conditions of the three-machine power system 77 TABLE 4.3 C r i t i c a l c l e a r i n g time by the time sol u t i o n and the proposed methods for the 3-machine power system 80 TABLE 4.4 Data of the CIGRE 7-machine system 85 TABLE 4.5 C r i t i c a l c l e a r i n g time by the time sol u t i o n and the proposed methods for the 7-machine power system 86 - v i i -LIST OF FIGURES Page Figure 2.1 Pattern recognition c l a s s i f i c a t i o n procedure 19 Fiugre 2.2 One-machine i n i f i n i t e - b u s power system 24 Figure 2-3 C h a r a c t e r i s t i c s of KE and TM with load v a r i a t i o n 24 Figure 2.4 The decision functions for three d i f f e r e n t f a u l t locations 24 Figure 3.1 One-machine i n f i n i t e - b u s power system 32 Figure 3.2 The power angle curve ••• 32 Figure 3.3 The b i f u r c a t i o n set of the cusp catastrophe which represents the steady state s t a b i l i t y c r i t e r i a 37 Figure 3.4 The energy function s t a b i l i t y c r i t e r i a 39 Figure 3.5 One machine i n f i n i t e - b u s power system 47 Figure 3.6 The transient s t a b i l i t y l i m i t s given by the swallowtail catstrophe 47 Figure 3.7 Transient s t a b i l i t y region for a l l possible f a u l t locations and loading conditions 49 Figure 3.8 The transient s t a b i l i t y region i n terms of s t a b i l i t y l i m i t s and c r i t i c a l c l e a r i n g times 50 Figure 3.9 The equivalent two machine power system 53 Figure 4.1 The three-machine power system (WSCC) 76 Figure 4.2 The accuracy of Taylor series during the f a u l t on period using only the second order term 78 - v i i i -Page Figure A.3 The ac c e l e r a t i o n power for d i f f e r e n t f a u l t l o c a t i o n s . 79 Figure 4.4 The transient s t a b i l i t y region of the three-machine power system using the dynamic equivalent method 81 Figure 4.5 The transient s t a b i l i t y region for the three-machine system by the Raylor series method 82 Figure 4.6 The CIGRE test system 83 Figure 4.7 The transient s t a b i l i t y region for the CIGRE test system using the dynamic equivalent method 87 Figure 4.8 The transient s t a b i l i t y region for the CIGRE test system using the Taylor series method 88 Figure 5.1 E f f e c t of damping i n c l u s i o n on the accelerating power during fault-on period 96 Figure 5.2 E f f e c t of damping i n c l u s i o n on the accelerating power during fault-on period 96 Figure 5.3 The phasor diagram of the generator i n t e r n a l q uantities for the transient state 98 Figure 5.4 The t y p i c a l e x c i t a t i o n voltage-time response 100 Figure 5.5 The e f f e c t of f l u x decay and e x c i t a t i o n response on the accelerating power during the transient period... 106 - i x -Page Figure A . l Cross section of the swallowtail b i f u r c a t i o n set for U > 0 120 Figure A.2 Cross section of the swallowtail b i f u r c a t i o n set for U < 0 120 Figure A.3 The 3-D b i f u r c a t i o n set of the swallowtail catastrophe 120 Figure B.l The cusp manifold and i t s b i f u r c a t i o n set 123 Figure B.2 The cusp p o t e n t i a l function V(X) at d i f f e r e n t values of control variables 125 x -ACKNOWLEDGEMENT I would l i k e to thank Dr. M.D. Wvong for his continual encouragement and invaluable guidance throughout the course of this research. I would l i k e to express my sincere appreciation and indebteness to Dr. H.W. Dommel for his encouragement and f i n a n c i a l assistance. I am also indebted to the B.C. Hydro and Power Authority for the f i n a n c i a l support and to the E l e c t r i c a l Engineering Department for the computing support and the teaching assistantships provided. I wish to thank my wife, Najwa, for her encouragement and patience throughout my graduate program. Thanks are extended to G a i l Schmidt for her pleasant attit u d e and s k i l l e d typing of the manuscript. - x i -1 CHAPTER 1 INTRODUCTION The s t a b i l i t y of a power system implies that a l l i t s generators remain i n synchronism through normal and abnormal operating conditions. Transient s t a b i l i t y a r i s e s when a large disturbance such as a loss of generation, load or transmission l i n e s takes place i n the power system. The question of whether the power system w i l l s e t t l e to a new stable operating state or not i s known as the transient s t a b i l i t y problem. The ever increasing demand for e l e c t r i c a l energy requires ever larger interconnected power systems and a maximum generation. This raises great concern about the security of power systems when subjected to large d i s t u r -bances. Transient s t a b i l i t y , therefore, becomes an increasingly important consideration i n system planning and operation. Extensive s t a b i l i t y studies are needed i n order to ensure system security before a planning or operating d e c i s i o n i s made. Each contingency for each disturbance considered requires a large number of s t a b i l i t y studies to determine the c r i t i c a l c l e a r i n g time or system s t a b i l i t y l i m i t s . A t y p i c a l transient s t a b i l i t y study consists of obtaining the time s o l u t i o n to the power system d i f f e r e n t i a l and algebraic equations s t a r t i n g with the system conditions p r i o r to the transient. The power system equa-tions should include a l l s i g n i f i c a n t parameters that influence s t a b i l i t y such as generator controls, s t a b i l i t y controls and protective devices. The desired objectives of a transient s t a b i l i t y analysis are: 2 i . The s t a b i l i t y of the power system. Is i t stable or not, to what degree i s i t stable, and how far i s i t from the s t a b i l i t y l i m i t s ? i i . Time responses of generator v a r i a b l e s , bus voltages, currents, and active and reactive power. i i i . System quantities that a f f e c t the performance of protective devices. The above objectives are key issues i n power system design, planning and operation to ensure system s t a b i l i t y for d i f f e r e n t prescribed disturbances. The time s o l u t i o n method of s t a b i l i t y analysis, although i t i s very r e l i a b l e , accurate, and suitable for d i f f e r e n t modelling orders, has the following disadvantages; 1. The method involves numerical i n t e g r a t i o n of a large number of d i f f e r -e n t i a l equations for each disturbance considered. A large number of r e p e t i t i v e simulations i s required for each case to determine either the s t a b i l i t y l i m i t s or the c r i t i c a l c l e a r i n g time. This procedure Is very time consuming i n the system planning stage where a large number of cases need to be considered. 2. In system operations, there are si t u a t i o n s where f a s t s o l u t i o n i s need-ed to make operational decisions. These si t u a t i o n s could be d i f f e r e n t from those previously considered during planning. Since the time s o l u -t i o n method i s slow, the system operator has either to overreact to ensure system security or to make decisions that may put s t a b i l i t y of the system at r i s k . 3. The power system operating conditions change during the' course of the day and year, while s t a b i l i t y studies are done o f f - l i n e for ce r t a i n 3 severe cases. This leads to improper decisions i n some cases and, therefore, may increase operational expenditures. From the above, i t i s clear that the power industry greatly needs an a l t e r n a t i v e method to solve the transient s t a b i l i t y problem. The alt e r n a -t i v e method cannot e n t i r e l y replace the time sol u t i o n , but i t should reduce the number of simulations for each case and hence save a great deal of computation time. A new method i s also needed for system operations where the system operator has only minutes or hours to make an important opera-t i o n a l decision. The best s o l u t i o n , of course, i s to have an on-line method that deals with system operations on a real-time b a s i s . The desired method for fast analysis of transient s t a b i l i t y should s a t i s f y the following important requirements. 1. Provide fast and r e l i a b l e answers to the transient s t a b i l i t y problem when a s p e c i f i e d disturbance i s given. The answer must p a r t i c u l a r l y indicate whether the system i s stable or not. 2. Include s u f f i c i e n t modelling options so that the answer provided i s i n the range of the actual system response. 3. Provide the necessary information about the degree of s t a b i l i t y and system se c u r i t y so that the operator can make proper decisions to ensure system s e c u r i t y . 4. If th i s fast method i s to be used for real-time system operations, i t should be adaptable to changes i n operating conditions, d i f f e r e n t disturbances, and s t a b i l i t y c o n t r o l s . Although extensive research has been conducted i n t h i s area, l i t t l e of the previous requirements have been achieved so f a r . A great deal of 4 research i s s t i l l needed to ei t h e r improve the e x i s t i n g methods or propose new methods i n order to f u l f i l l a l l these requirements. This research i s motivated by the challenging problem of developing a s a t i s f a c t o r y fast d i r e c t method for the assessment of transient s t a b i l i t y . The main objectives of this thesis are: 1. To develop a f a s t , r e l i a b l e d i r e c t method for the assessment of trans-ient s t a b i l i t y of power systems. 2. To improve the power system model used by including a l l necessary options to get as accurate a system response as possible. 3. To predict the comprehensive transient s t a b i l i t y region with security boundaries for possible on-line assessment of transient s t a b i l i t y . A review of the l i t e r a t u r e on the transient s t a b i l i t y problem and presently a v a i l a b l e solutions are given i n Chapter 2. In Chapter 3 a new d i r e c t method for s t a b i l i t y assessment i s introduced with preliminary a p p l i -cations on a simple one-machine i n f i n i t e - b u s power system. In Chapter 4, the technique i s applied to multi-machine power systems using dynamic equivalents and Taylor series expansions. Transient s t a b i l i t y regions are al s o given by using the catastrophe theory. Chapter 5 presents the cap-a b i l i t y of the new method to include e x c i t a t i o n , f l u x decay and system damping i n the power system model. Chapter 6 summarizes the conclusions and achievements of t h i s project. 5 CHAPTER 2 REVIEW OF STABILITY PROBLEM AND ITS SOLUTION 2*1 Introduction The s t a b i l i t y problem of power systems has been given new importance since the famous blackout i n Northeast U.S.A. i n 1965. Considerable research e f f o r t has gone into the s t a b i l i t y i n v e s t i g a t i o n both for o f f - l i n e and on-line purposes. I t started a chain reaction a f f e c t i n g planning, operation and control procedures of e l e c t r i c power systems. Since then several studies have been conducted and new concepts and di r e c t i o n s have been suggested to prevent i n s t a b i l i t y and ensure security and r e l i a b i l i t y of power systems. In t h i s chapter, a b r i e f review of the s t a b i l i t y problem i s presented, followed by a l i t e r a t u r e review of the d i f f e r e n t methods used and suggested to analyze i t . 2.2 The S t a b i l i t y Problem A stable power system implies that a l l i t s interconnected generators are operating i n synchronism with the network and with each other. Problems a r i s e when the generators o s c i l l a t e because of disturbances that occur from transmission f a u l t s on switching operations. There are two types of s t a b i l i t y problem i n power systems. F i r s t l y , the steady state s t a b i l i t y problem which refers to the s t a b i l i t y of power systems when a small disturbance occurs i n the power systems such as a gradual change i n loads, manual or automatic changes i n e x c i t a t i o n , 6 i r r e g u l a r i t i e s i n prime-mover input. Obviously these small disturbances cannot cause loss of synchronism unless the system i s operating at, or very near, i t s steady state s t a b i l i t y l i m i t . This l i m i t i s the greatest power that can be transmitted on a s p e c i f i e d c i r c u i t , under c e r t a i n operating conditions i n the steady state, without loss of synchronism. The analysis of steady state s t a b i l i t y requires the so l u t i o n of power flow equations and swing equations over a period of a few miniutes. Governor and exciters should also be included i n the steady state s t a b i l i t y a n a l y s i s . Secondly, there i s the transient s t a b i l i t y problem which arises from a large disturbance i n the power system such as a sudden loss of generators or loads, switching operations, or f a u l t s with subsequent c i r c u i t i s o l a t i o n . Such large disturbances create a power unbalance between supply and demand i n the system. This unbalance takes place at the generator shafts and causes the rotors to o s c i l l a t e u n t i l a new steady state operating point i s reached; or u n t i l the rotors continue to o s c i l l a t e and deviate from each other and f i n a l l y some generators w i l l lose synchronism. Loss of synchronism must be prevented or c o n t r o l l e d because i t has a d i s t u r b i n g e f f e c t on voltages, frequency and power, and i t may cause serious damage to generators, which are the most expensive elements i n power systems [1 ] . The generators which tend to lose synchronism should be tripped, i . e . disconnected from the system before any serious damage occurs, and subse-quently brought back to synchronism. While t h i s can be done r e a d i l y with gas and water turbine generators, steam turbine generators require many hours to r e b u i l d steam and the operator has to shed load to compensate for the loss of generators. Loss of synchronism may also cause some protective 7 relays to operate f a l s e l y and t r i p the c i r c u i t breakers of unfaulted l i n e s . In such cases the problem becomes very complicated and may r e s u l t i n more generators l o s i n g synchronism. 2.3 Basic Power System S t a b i l i t y Model The equation of motion of n generator rotors i n a power system of n-machines i s given by , i = 1,2, . M.6. + D. 6 = p - P = P i i ' i i i e i a i n where e i = 2 ( G E E . cos 6 + B, . E. E. s i n 6 ) j=l i j i J i j i j i j i j ' (2.1) (2.2) M. 10 e i a i G ,B i j * i j 6 i j V E j i n t e r n a l rotor angle of machine 2. to 1^ = i n e r t i a constant, moment of i n e r t i a , angular speed, damping c o e f f i c i e n t , mechanical power input, e l e c t r i c a l power output, ac c e l e r a t i n g power. r e a l and imaginary parts of reduced nodal admittance matrix. 6 - 6 . i j i n t e r n a l voltages of machine i , j . Under steady s t a t e c o n d i t i o n s , P .=0, and 6. i s constant. When a a i i {• d i s t u r b a n c e o c c u r s , P . becomes d i f f e r e n t from zero and equation (2.1) a i d e s c r i b e s the behaviour of 6 w i t h time. For machine i to be stable, 6 8 must s e t t l e to a constant value again and P . must return to zero. For a power system of n machines, n swing equations have to be solved i n order to decide whether the power system i s stable or not when a s p e c i f i c disturbance occurs.. The s o l u t i o n of the swing equations obtained depends upon the model of the power system elements such as generators (machine, exciter and governor), transmission l i n e s and loads. Furthermore, the choice of the power system model w i l l depend upon the s t a b i l i t y study to be c a r r i e d out and the period of analysis [2]. 2.3.1 Synchronous Generator The most s i m p l i f i e d model of the sychronous generator i s the so-c a l l e d c l a s s i c a l model. In t h i s model the generator i s described by a constant voltage behind transient impedance. This model i s acceptable for the f i r s t swing transient s t a b i l i t y analysis which has a period of one second or less [2]. The i n t e r n a l voltage i s constant i n magnitude only. This representation neglects the e f f e c t of saliency and assumes constant f l u x linkages and a small change i n speed. The voltage equation i s given by |E| / 6 - Vfc + r a I t + j x- l t (2.3) where | E | : magnitude of voltage behind transient impedance. 6: angle of |E| ( v a r i a b l e ) . V t terminal voltage. r armature resistance. a transient reactance. I terminal current. t 9 When saliency and changes i n f i e l d f l u x linkages are taken into account, i t i s no longer possible to represent the generator completely by an equivalent c i r c u i t . The r e s u l t i n g equations have time-dependent c o e f f i c i e n t s unless they are transformed by the use of Park's transformation into the so-called d i r e c t - and quadrature-axis variables [2]. The model then involves both d i f f e r e n t i a l and phasor-algebraic equations [3]. The d i f f e r e n t i a l equation i s : i dE , - 1 = J - ( K - . - O (2.4) dt ' f d 1  Ado i » where E^ i s the v o l t a g e p r o p o r t i o n a l to f i e l d f l u x l i n k a g e , T^ o i s the d i r e c t - a x i s transient open-circuit time constant, E ^ i s the quadrature-axis f i e l d v o l t a g e and E^ . i s the v o l t a g e p r o p o r t i o n a l to f i e l d current. The phasor equations of the above quantities are: E* = E - j ( x - x ' ) l , (2.5) q q q d d I t a t d d <}q E = V + r I + jx I (2.7) q t a t q t where E^ i s the voltage behind quadrature-axis reactance x^, 1^ , 1^  are the components of I f c along the d i r e c t - and quadrature-axes, r e s p e c t i v e l y . Here the machine angle 6 i s the angle of E^. The above model i s desirable i n the transient s t a b i l i t y analysis for the cases when a longer period of analysis i s needed when a high-speed voltage-regulator i s considered [2]. 10 2.3.2 E x c i t e r and Governor Control Systems The e x c i t e r i s a device that supplies f l u x to the synchronous generator and d i r e c t l y controls the output voltage of the generator. The e x c i t e r control system provides the proper f i e l d voltage to maintain a desired system voltage during transient and steady state operations. There are many types of e x c i t e r control systems i n use i n power systems [2]. The basic components of an e x c i t e r control system are the regulator, amplifier and e x c i t e r . The regulator measures the actual regulated voltage and determines the voltage deviation. The deviation s i g n a l i s then amplified to provide the s i g n a l required to c o n t r o l the e x c i t e r f i e l d current that changes the exciter output voltage and hence r e s u l t s i n a new e x c i t a t i o n l e v e l for the generator. Reference [3] gives the d e t a i l e d d i f f e r e n t i a l equations of the e x c i t e r control system. The e x c i t e r c o n t r o l system must be Included i n the transient s t a b i l i t y studies i f the period of analysis i s longer than one second. Modern fast h i g h - c e i l i n g e x c i t e r s operate even within a period of one second and i f such an e x c i t e r i s considered then the e x c i t e r c o n t r o l system must be included i n the generator model. The rotor speed of the synchronous generator i s c o n t r o l l e d by the governor control system by varying the mechanical power input. The governor c o n t r o l system can be taken into account i n the s t a b i l i t y studies by solving i t s d i f f e r e n t i a l equations simultaneously with the r e s t of the system equations. I t i s necessary to include the governor e f f e c t i f the period of analysis Is longer than one second [3]. 11 2.3.3 Transmission System and Loads Transmission l i n e s are usually represented by nodal type equations. Current equations are written f or each node i n the following manner: n 1. = E Y E (2.8) j-1 J J where E. are nodal voltages, Y. . are the c o e f f i c i e n t s of the i ' * 1 row of the J i j network admittance matrix and 1^ i s the input current into node i due to a generator or load. The power input into the network at node i i s given by P e i - V E i O (2-9) Loads are usually represented as constant complex power, constant impedance or constant current at constant power fa c t o r . The equation for constant complex power i s S± = P i + j Q ± = E± li = constant (2.10) for constant impedance i s E i Z. = y~ = constant (2-11) i and for constant current at constant power factor i s Transmission l i n e s and loads a f f e c t d i r e c t l y the s t a b i l i t y l i m i t s [4]. The types of load models mentioned above are not adequate models i f the load i s of a dynamic type. I t i s very important to use more advanced load models 12 for dynamic loads e s p e c i a l l y when high voltage v a r i a t i o n s are expected [5,6]. 2.4 Solution of the Transient S t a b i l i t y Problem Transient s t a b i l i t y refers to the amount of power that can be stably transmitted when the power system i s subjected to a large disturbance such as, f a u l t s , loss of l i n e s , loss of generators or loads and sudden changes i n the t i e - l i n e flow. The main objective of a transient s t a b i l i t y study i s to assess the possible loss of synchronism due to such large disturbances. Transient s t a b i l i t y studies enable power engineers to set up protective devices and to design proper s t a b i l i t y controls that ensure no loss of synchronism. The growth of large interconnected power systems demands c a r e f u l transient s t a b i l i t y studies because some large disturbances could have catastrophic r e s u l t s . On the other hand, rapid growth makes i t extremely d i f f i c u l t to carry out such c a r e f u l and d e t a i l e d s t a b i l i t y studies. An e a s i l y reached conclusion i s that an exact answer i s d i f f i c u l t to obtain [7]. There are three main approaches to solving the transient s t a b i l i t y problem. These approaches w i l l be discussed i n the following sections. 2.4.1. Numerical Integration Methods These methods solve the power system d i f f e r e n t i a l equations (swing equations) during and a f t e r the transient period. From the response of 13 the swing curves of the generators i n the power system an experienced power engineer can assess the s t a b i l i t y of the power system. The numerical i n t e g r a t i o n methods can be divided into two basic types: 1. One-step methods: these methods require only information concerning the previous time point to ca l c u l a t e the values at the next time point. Consider a set of non-linear d i f f e r e n t i a l equations ^ • = f ( y , t ) (2.13) where y represents a set of n variables and f a vector of n functions. The one-step method calculates at point 1 s t a r t i n g from point 0 using a s t r a i g h t l i n e defined by the d e r i v a t i v e at point 0 yi = yo + d r l o h ( 2 a 4 ) where h i s the time i n t e r v a l , t ^ - t ^ . The new point can be used to c a l c u l a t e the next one, and the i t e r a t i o n continues u n t i l the sol u t i o n i s obtained. These methods have high accuracy when the time i n t e r v a l Is small, but the accuracy decreases r a p i d l y as the time i n t e r v a l increases. Euler's method and the step by step method given i n [1] are t y p i c a l one-step methods. 2. Multi-step methods: these methods require for each new point, information about two or more preceding points. The predictor-corrector method [8] calculate a point y as follows 14 *n+l * }n \ V k + bk f ^ tn-k ) ( 2 ' 1 5 ) k=0 k=-l where appropriate values for the a and b constants are chosen for the pre d i c t o r and corrector formulas. f ( y , t ) indicates the derivative with respect to time. Equation (2.15) i s the same as equation (2.14) for the one-step method with a =b =1, and a l l other constants equal to zero. r n n ' n The number series approach [3] Is another type of multi-step method where the generator angles versus time are obtained using the following i t e r a t i v e equation 6 + - 6 + i ^ L + i ^ P _ (2.16) n+2 n w a n+1 Equation (2.16) i s derived from the swing equation by evaluating 6(t) i n t e g r a t i o n by the trapezoidal r u l e . Other methods have been applied to t h i s problem such as Runge-Kutta method [10], an improved Gauss-Seidel approach [11], and the phase-plane techniques [12]. The main advantage of multistep methods i s that long time i n t e r v a l s can be used with high accuracy and s t a b i l i t y . In general, the numerical i n t e g r a t i o n methods are widely accepted and used because of the fact that these methods provide det a i l e d r e s u l t s and are capable of including d e t a i l e d models of power system elements. The main disadvantage of the numerical int e g r a t i o n methods for transient s t a b i l i t y studies are the formidable computation time required, the need for human i n t e r p r e t a t i o n of the swing curves f o r assessing the transient s t a b i l i t y and the fact that the method cannot be used for on-line assessments of s t a b i l i t y . 15 In spite of the high speed computation of modern computers and several refinements and developments of the technique [13-15] the d i g i t a l a n a lysis of even a moderate sized power system i s s t i l l too time consuming. Analog and hybrid simulation [16] have been suggested. In these schemes, a l l the d i f f e r e n t i a l and a l g e b r i c equations are solved i n analog form and the rest of the work ( i n i t i a l i z a t i o n , switching and l o g i c operations) i s taken care of by the d i g i t a l computer. But t h i s approach i s s t i l l too slow for on-line a p p l i c a t i o n s . 2.4.2 Lyapunov's Direct Method To improve on the slow computational speed of numerical int e g r a t i o n methods, Lyapunov's s t a b i l i t y theory has been applied to the transient s t a b i l i t y problem. Let the system be represented by a set of d i f f e r e n t i a l equations dX. ^JL = f ± ( X 1 § X 2, ... X n, t ) , i = l , . . . n . (2.17) The o r i g i n of the s t a t e v e c t o r X i s c o n s i d e r e d to be s t a b l e , i . e . f i ( 0 , 0 ... 0,t) =0, i = l , . . . n . When the set of i n i t i a l conditions at the s t a r t of the disturbance are a l l known, Lyapunov's d i r e c t method w i l l determine, without f i n d i n g the actual t r a j e c t o r y of the state v a r i a b l e s , i f the systm w i l l asymptotically reach a stable point. Therefore, much integration time can be saved as s t a b i l i t y can be determined d i r e c t l y from i n i t i a l conditions.;. The method assumes a s c a l a r , continuous Lyapunov function V(X,t) exi s t s i n the neighborhood of the o r i g i n such that, 16 1. V(0,t) - 0 2. V(X,t) > 0 X e R, X * 0 3. f f f i d l <0 X e R dt where R i s a region around the stable point X=0 and i s c a l l e d the region of s t a b i l i t y . I f the above conditions are v a l i d for V(X,t), i t can be shown that for a disturbance X e R (X * 0) then X e R as t -*• fl0. The main d i f f i c u l t y here i s that the determination of the boundary of the s t a b i l i t y region i s very time consuming. The a p p l i c a t i o n of t h i s method to the transient s t a b i l i t y problem of power systems can be summarized i n the following steps [17-19]. 1. I n i t i a l conditions are determined by a load flow study. 2. The admittance matrix of the power system i s reduced to the number of g e n e r a t i n g buses f o r d u r i n g - d i s t u r b a n c e and p o s t - d i s t u r b a n c e conditions. 3. The steady state equilibrium point of the post-disturbance system i s found by solving the equation P., - P = 0 , i = 1, ... n i e i where P^ and P g^ are the mechanical input power and the e l e c t r i c a l output power, re s p e c t i v e l y . The same equation has also to be solved to f i n d the closest unstable equilibrium point. Such a point helps to define the region of s t a b i l i t y around the stable state point. 4. The Lypanov's function V(X,t) i s found and c a l c u l a t e d . Thus the s t a b i l i t y regions for V=constant can be found. 17 5. Integration of the d i f f e r e n t i a l equations for the during-disturbance period. At each new point, the function V i s calculated with post-disturbance conditions to determine whether the system i s stable when the disturbance i s cleared. The f i r s t instance of i n s t a b i l i t y defines the c r i t i c a l c l e a r i n g time for the disturbance. Numerous methods have been suggested to derive a Lyapunov function with the desired properties [20,21]. I t should be noted here that d i f f e r e n t V functions may y i e l d d i f f e r e n t answers i n the sense that they comprise d i f f e r e n t subregions of the s t a b i l i t y region around the stable point [15]. The transient energy function method [22] i s the most promising approach among the d i r e c t methods. Remarkable progess has been made i n recent years [23,24]. However, the method i s s t i l l not suitable for on-line a p p l i c a t i o n because of some p r a c t i c a l and computational d i f f i c u l t i e s [25]. The advantages which Lyapunov's d i r e c t method has over numerical i n t e g r a t i o n methods are the following: 1. This method i s computationally f a s t e r . Here the c r i t i c a l c l e a r i n g time can be obtained i n a single i n t e g r a t i o n while the simulation technique requires r e p e t i t i v e Integration for d i f f e r e n t assumed c l e a r i n g times. 2. For a given c l e a r i n g time, t h i s method w i l l indicate whether the system i s stable or not. This s o l u t i o n i s d i r e c t and i t does not require the i n t e g r a t i o n of the system d i f f e r e n t i a l equations beyond the c l e a r i n g time. 3. I t i s shown that the value of the Lyapunov function and the c r i t i c a l c l e a r i n g time are related [26]. Hence, this method can y i e l d s t a b i l i t y margins or s t a b i l i t y indices without making actual c a l c u l a t i o n s . Looked 18 at from another angle, t h i s method can indicate whether a c l e a r i n g time i s adequate or not. Despite the above a t t r a c t i v e features, Lyapunov's d i r e c t method has not received general acceptance from the u t i l i t i e s . The main reasons for t h i s are the following p r a c t i c a l d i f f i c u l t i e s associated with t h i s method: 1. Approximate mathematical model of the system; many s i m p l i f y i n g assumptions have been made i n a r r i v i n g at a mathematical model of the power system to make i t suitable for i t s analysis by Lyapunov's d i r e c t method. The c r i t i c s of t h i s method question the v a l i d i t y of the s i m p l i f i c a t i o n s which are necessary f o r construction of a proper Lyapunov function. Moreoever, the Lyapunov functions are not unique and several such functions can be constructed by using d i f f e r e n t dynamic models of the power system. 2. The method i s conservative and, although some stable point are r e a d i l y i d e n t i f i e d , others are inconclusive. Also, the method cannot predict i n s t a b i l i t y , therefore, i t may produce too many f a l s e alarms. 3. Determination of the region of s t a b i l i t y ; from a p r a c t i c a l a p p l i c a t i o n point of view, t h i s has been the most consuming part of th i s method. In power systems, t h i s requires determination of the unstable equilibrium points closest to the post f a u l t stable equilibrium points, and involves the s o l u t i o n of n-nonlinear algebraic equations. 2.4.3 Pattern Recognition Method * The object of th i s method i s to determine a function S(X) [27] such that: 19 > 0 for stable X S(X) - { < 0 for unstable X S(X) i s c a l l e d the decision function and X i s the state variable , where (X — , X£ ••• X^). The c l a s s i f i c a t i o n procedure i n the pattern recognition method i s shown i n F i g . (2.1) and summarized i n the following steps: Data_* Tra ining set Feature extraction Training procedure Decision Figure (2.1) 1. Training Set: In order to obtain S(X), a t r a i n i n g set of various operating condi-tions must be a v a i l a b l e . Each operating condition or state i s s p e c i f i e d by variables such as i n j e c t i o n powers, loads, generation powers, load flows, voltage magnitudes and rotor angles. Each v a r i a b l e describing a p a r t i c u l a r X condition constitutes a component of the pattern vector X = ( X ^ Xj ••• x n ) * I d e a l l y , every c o n c e i v a b l e p a t t e r n should be included i n the t r a i n i n g set so that i t covers the whole spectrum of the system operating conditions. Each pattern of the t r a i n i n g set must be c l a s s i f i e d whether i t i s stable or unstable. 20 2. Feature extraction: The number of variables obtainable from the t r a i n i n g set can be very large. In pattern recognition, i t i s not desirable nor i s i t necessary to use a l l the av a i l a b l e variables to obtain the c l a s s i f i e r function. Therefore, we have to choose the most Important variables that contain the necessary information about s t a b i l i t y . This can be done by p r a c t i -c a l experience with the power system or by using s t a t i s t i c a l measures. The following function F [28] provides an easy and straightforward measure of the information i n each v a r i a b l e . s where m s = mean of the v a r i a b l e i n the stable c l a s s . m u = mean of the v a r i a b l e i n the unstable c l a s s . cr s = variance of the v a r i a b l e i n the stable c l a s s . a u = variance of the v a r i a b l e i n the unstable c l a s s . Feature extraction begins with the computations of the F values for a l l the components of the pattern vectors In the t r a i n i n g set. The vari a b l e with the largest F i s selected as the f i r s t feature. Then the redundant information i s removed by c a l c u l a t i n g the c o r r e l a t i o n factor between the f i r s t feature and the rest of the var i a b l e s , excluding variables with high c o r r e l a t i o n factor with the f i r s t feature. The procedure of c a l c u l a t i n g F with the rest of the features Is repeated and se l e c t i o n continues [27]. 21 3. Training Procedure: Having selected the s i g n i f i c a n t features, the next step i s to obtain the decision function S(X). This function can be f i r s t order such as S(X) = W + W. X, + ... + W X o 1 1 m m The weighting c o e f f i c i e n t s (WQ, ... Wffi) are determined such that S(X) > 0 i f X i s stable and S(X) < 0 i f X i s unstable There are many t r a i n i n g procedures a v a i l a b l e such as the least squares method [29,30], l i n e a r programming [30,31] and an optimal search algorithm [32]. Another possible c l a s s i f i e r function i s of second order or quadratic form [27] given below m m m S(X) = W + Z W, Z. + Z Z W,. Z, Z, ° 1-1 1 1 1=1 j-1 i j 1 i A f t e r the weights are determined, t h e i r performance are checked by using them to c l a s s i f y the patterns i n the t r a i n i n g set. This method has been improved [33-34], and i t i s found that non-para-metric approaches are more r e l i a b l e than the parametric approaches mentioned above. The non-parametric approaches r e l y on experience to select the proper variables to be used i n s t a b i l i t y c l a s s i f i c a t i o n [35]. The main advantages of pattern recognition are: 1. I t i s suitable for on-line assessment of transient s t a b i l i t y i f a r e l i -able decision function i s determined. 22 2. The method i s independent of the model of the power system that i s used. The main disadvantages of pattern recognition are: 1. O f f - l i n e computations to obtain the decision function are excessive. 2. The decision function S(X) may very well be s e n s i t i v e to system configuration so that a d i f f e r e n t S(X) i s needed for any change i n the system configuration. 3. Correct c l a s s i f i c a t i o n of 100% cannot be achieved with large power systems because i t involves unlimited o f f - l i n e computations. 2.5 Discussion of E x i s t i n g Methods The growth of large interconnected power systems dic t a t e s the necessity for on-line transient s t a b i l i t y assessment to ensure the security of the power system and maintain r e l i a b l e service. The^ main problem of large power systems i s that i t i s impossible to simulate exactly for on-line purposes. Most d i r e c t methods provide global s o l u t i o n to the s t a b i l i t y problem with high approximation, except the method using the energy-type Lyapunov function [24] which y i e l d s s o l u t i o n for the i n d i v i d u a l machine. But with t h i s method the problem Is how to define the s t a b i l i t y region for an i n d i v i d u a l machine and that only a s i m p l i f i e d power system model i s used. One may expect better r e s u l t s from the pattern recognition method because of the fact that i t i s independent of the model used t o simulate the power system. This method can be used for an on-line assessment of s t a b i l i t y even for i n d i v i d u a l machines, but i t cannot be generalized because 23 each decision function i s good only for the system trained f o r . The motivation for t h i s work i s to search for a method which i s su i t a b l e for an on-line assessment of transient s t a b i l i t y that can be used for every i n d i v i d u a l generator i n the power system. I f such a method can be found then the problem of transient s t a b i l i t y can be dealt with on-line. The most suitable method (from among the e x i s t i n g methods) for such an a p p l i c a t i o n i s pattern recognition because of the fac t that i t i s independent of the systems equations which couple the generators together. 2.4.1 Ap p l i c a t i o n of Pattern Recognition to a Simple Power System A recent technique of pattern recognition [35] i s applied to a single machine-Infinite bus system shown In F i g . 2.2. Two indices are used as state v a r i a b l e s , the k i n e t i c energy (K.E) and transmission margin (TM), where K.E = - MW2 2 and TM = P - P max c M = i n e r t i a constant of the machine W = speed p = maximum power for pos t - f a u l t condition max P = l i n e power at the instant a f a u l t i s cleared c K.E and TM have opposite c h a r a c t e r i s t i c s with loading conditions as shown i n F i g . 2.3. The decision function S(X) i n terms of K.E. and TM can be written i n the form j . 16 H - 3.0 | E | « 1.25 pu 2d • j 0.1.6 pu j.52 j.52 T, Z T, j.16 0 Fig. 2.2 One machine infinite-bus power system KE(%) 1 0 0 90 80 70 60 50 TM(%) 110 105 100 70 80 90 100 (%) L o a d Fig. 2.3 Characteristics of KE and TM with load variation T M KE Fig. 2.4 The decision function for three different fault locations 25 S(.X) = W z + <t> x o S(X) > 0 : stable S(X) < 0 : unstable where X = (K.E, TM), W i s c o e f f i c i e n t vector <t> = constant o The technique i s applied to the system shown i n F i g . 2.2 for three d i f f e r e n t f a u l t locations indicated. Three d i f f e r e n t decision functions were determined for the three f a u l t locations as shown i n F i g . 2.4. For f a u l t l o c a t i o n (1), S i ( X ) = -TM + 20.71 K.E + 2.44. For f a u l t l o c a t i o n (2), S 2 ( X ) = -TM + 31.51 K.E 4 3.09 For f a u l t l o c a t i o n (3), S 3 ( X ) = -TM - 681.8 K.E + 56.1. In conclusion, to represent a l l possible f a u l t locations we have to b u i l d up a surface that consists of a l l possible f a u l t l o c a t i o n s . This surface should contain a large number of l i n e s i n order to cover the whole region of s t a b i l i t y . Therefore, the decision function S (X) should represent the equilibrium surface of the s t a b i l i t y region. This conclusion brought up the idea of using catastrophe theory to define the s t a b i l i t y region. 26 CHAPTER 3 APPLICATION OF CATASTROPHE THEORY TO POWER SYSTEM STABILITY A preliminary a p p l i c a t i o n of the catastrophe theory, a mathematical method, to the s t a b i l i t y assessment of power systems i s introduced i n t h i s chapter. The theory i s reviewed i n Section 3.2 followed by applications to both steady state and transient s t a b i l i t y of a simple one-machine i n f i n i t e - b u s power system. This method provides a very good t o o l to v i s u a l i z e steady state and transient s t a b i l i t y regions. 3.1 Introduction The sustained increasing demand for e l e c t r i c a l power requires larger interconnected power systems and operation at or near to f u l l capacity. Therefore, transient s t a b i l i t y of power systems becomes a major factor i n planning and day-to-day operations and there i s a need for f a s t on-line s o l u t i o n of transient s t a b i l i t y to predict any possible loss of synchronism and to take the necessary measures to restore s t a b i l i t y . Lyapunov's d i r e c t method and pattern recognition have been introduced f o r f a s t assessment of transient s t a b i l i t y and eventually to implement these methods for on-line a p p l i c a t i o n s . Catastrophe theory has been applied to the study of s t a b i l i t y of various dynamic systems [36] and i n recent years to the steady state s t a b i l i t y problem of power systems by Sallam and Dineley [37]. An a t t r a c t i v e feature of catastrophe theory i s that the s t a b i l i t y regions are 27 defined i n terms of the system parameters bounded by the l i n e s of s t a b i l i t y l i m i t s . The a p p l i c a t i o n of catastrophe theory i s extended to the transient s t a b i l i t y problem. A well defined transient s t a b i l i t y region i n terms of system parameters i s obtained; these regions are suitable for fast on-line a p p l i c a t i o n s . 3.2 Catastrophe Theory I t i s a natural phenomenon that sudden changes can occur as a r e s u l t of smooth or gradual changes. Examples might include the breakdown of an i n s u l a t o r when voltage i s b u i l t up gradually, the collapse of a bridge by gradual load increases and the loss of synchronism of generators i n a power system when subject to smooth changes i n operating conditions. The term "catastrophe" i s used for such sudden changes that are caused by smooth a l t e r a t i o n s . C atastrophe theory was o r i g i n a l l y p resented by the French mathematician Rene Thorn and published i n h i s book "Str u c t u r a l S t a b i l i t y and Morphogenesis" [38]. Thorn used d i f f e r e n t i a l topology to explain sudden changes i n morphogenesis. The idea then attracted the attention of many s c i e n t i s t s and has been applied to a v a r i e t y of sciences [36]. Catastrophe theory i s a theory that explores the region of sudden changes i n dynamic systems and deals with the properties of d i s c o n t i n u i t i e s d i r e c t l y . I t has been defined as the study from a q u a l i t a t i v e point of view of the ways the solutions to d i f f e r e n t i a l equations may change [39]. 28 Consider a system whose behaviour i s usually smooth but which sometimes (or i n some places) exhibits d i s c o n t i n u i t i e s . Suppose the system has n state variables and c o n t r o l l e d by m independent variab l e s , and suppose a smooth p o t e n t i a l function to describe the system dynamics e x i s t s . Given such a system, what catastrophe theory t e l l s us Is the following: The number of q u a l i t a t i v e d i f f e r e n t configurations of d i s c o n t i n u i t i e s that can occur depends not on the number of state v a r i a b l e s , which may be very large, but on the number of control v a r i a b l e s , which i s generally small. P a r t i c u l a r l y , i f the number of control variables i s not greater than four, then there are only seven d i s t i n c t types of catastrophes known as the seven elementary catastrophes, and i n none of these are more than two state variables involved [40]. 3.2.1. Catastrophe Theory and B i f u r c a t i o n Analysis Consider a continuous p o t e n t i a l function V(X,C) which represents the system behaviour, where X are the state variables and C are the control v a r i a b l e s . The p o t e n t i a l function V(V,C) can be mapped i n terms of i t s c o n t r o l variables C to define the continuous region. Let the p o t e n t i a l function be represented as V(X,C) : M X C * R (3.1) where M, C are m a n i f o l d s In the s t a t e space R° and the control space R T r e s p e c t i v e l y . Now we define the catastrophe manifold M as the equilibrium surface t h a t r e p r e s e n t s a l l c r i t i c a l p o i n t s of V(X,C). I t i s the subset R Q X R C defined by v x V C ( X ) = o < 3 - 2 > 29 where V C ( X ) = V(X,C) and V x i s p a r t i a l d e r i v a t i v e w i t h r e s p e c t to X. Equation (3.2) i s the set of a l l c r i t i c a l points of the function V(X,C). Next we fi n d the s i n g u l a r i t y set, S, which i s the subset of M that consists of a l l degenerate c r i t i c a l points of V. These are the points at which v x v c(x) = 0 and V x V C(X) = 0 (3.3) The s i n g u l a r i t y set, S, i s then projected down onto the control space R r by eliminating the state variables X using (3.3) and (3.2), to obtain the b i f u r c a t i o n set, B. The b i f u r c a t i o n set provides a proj e c t i o n of the s t a b i l i t y region of the function V(X,C), i . e . i t contains a l l non-degenerate c r i t i c a l points of the function V bounded by the degenerate c r i t i c a l point at which the system e x h i b i t s sudden changes when i t i s subject to small changes. The seven elementary catastrophes of r<4 are l i s t e d i n Table (3.1). The geometric analysis of the catastrophes that are used i n t h i s thesis are presented i n d e t a i l i n Appendix ( A l ) . A s i m p l i f i e d analysis of the seven elementary catastrophes Is given i n Reference [40]. 3.3 Applications to the Steady State S t a b i l i t y Problem The a p p l i c a t i o n of catastrophe theory to the steady state s t a b i l i t y of power systems was f i r s t introduced i n Reference [37]. However, that a p p l i c a t i o n was l i m i t e d only to s a l i e n t - p o l e type synchronous generators. When the same procedure was applied to c y l i n d r i c a l - r o t o r generators, Catastrophe Control Space Dimension State Space Dimension Function Catastrophe Manifold Fold 1 1 1/3 x 3-ax 2 x -a Cusp 2 1 x^-ax-1/2 bx 2 3 . x -a-bx Swallowtail 3 1 1/5 x 5-ax-l/2 b x 2 - l / 3 c x 3 4 , 2 x -a-bx-cx B u t t e r f l y 4 1 1/6 x 6-ax-l/2 bx 2-l/3 cx 3-dx 4 5 2 3 x -a-bx-cx -dx Hyberbolit 3 2 3 3 x +y +ax+by+cxy 2 3x +a+cy 3y2+b+cx E l i p t i c 3 2 3 2 2 2 x -xy +ax+by+cx +cy 2 2 3x -y +a+2cx -2xy+b+2cy Parabolic 4 2 2 4 2 2 x y+y +ax+by+cx +dy 2xy+a+2cx x2+4y3+b+2dy Table 3.1 The Seven-Elementary Catastrophes 31 u n j u s t i f i e d answers were obtained [ 3 7 ] . In t h i s section, a new procedure for applying catastrophe theory to the steady state s t a b i l i t y of c y l i n d r i c a l rotor generators i s presented using Taylor series expansion around the i n i t i a l steady state operating point. This section i s a complement to the work done by Sallam and Dineley and i t introduces a general procedure which i s applicable to a l l types of generators. Furthermore, the damping and governor e f f e c t s can also be included using t h i s method. 3 . 3 . 1 C y l i n d r i c a l - r o t o r I n f i n i t e - b u s Power System Consider the one-machine i n f i n i t e - b u s power system of F i g . 3 . 1 . The steady state output power P g i s given by P = P s i n 6 ( 3 . 4 ) e m where P = m X , d S -is--the machine i n t e r n a l voltage, V i s the i n f i n i t e - b u s voltage, i s the mchine reactance and 6 i s the rotor angle (angle of i n t e r n a l voltage E). The power angle curve i s given i n F i g . 3 . 2 . Suppose the input power to the machine i s increased smoothly from P Q to P^. The machine w i l l o s c i l l a t e between 6Q and 6 2 a s shown In F i g . 3 . 2 . I f the new operating angle i s higher than the maximum power l i m i t , the machine w i l l lose synchronism. The energy of o s c i l l a t i o n s between 6 and 6^ i s given by 62 f P d6 * ( 3 . 5 ) 6 a o 32 jo.72 g i n f i n t P bus E = 1.71 V = 1. / o . o 6 = 36.5 o jo.72 Fig. 3.1 60 6 1 62 Angle Fig. 3.2 The power-angle curve 33 P = P - P s i n 6 a i m (3.6) P i s the acc e l e r a t i n g power. We now expand e q u a t i o n (3.6) by Taylor series around 6^  at t=0, P a becomes P ( 2 ) t 2 P (t) = P ( 0 ) + P ( 1 > t + J—- + ... a a a 2 j (3.7) where (m) _ a , m d t o (3.8) Therefore P ( 0 ) - P - P s i n 6 a i m i (1) _ d P d t P cos 6 6 = 0 m o (2) d P (1) d t P s i n 6 6 m o P cos 6 6 m o (3.9) = -P cos 6 6 m 0 34 / Q \ D P (3) _ a (2) d t = P cos 6 6 J + p s i n 6 ( 2 6 ) 6 m o m o ,(4) _ a  a " dt + P s i n 6 6 6 - p cos 6 6 = 0 m o m o - - P cos 6 N , 6 ' * + 3P s i n 6. 6' m 0 m 0 0 where 6 = 0 at t = 0 Since the change i n the power input of the machine i s smooth and small, the Taylor seri e s expansion determines the exact region of o s c i l l a t i o n , from 6 ^ to 6 ^ . The ac c e l e r a t i n g power i s then given by p ( 2 ) t 2 p ( 2 ) t 4 p ( 6 ) t 6 P (t) = P ( 0 ) + a + - 5 + ± a v ' a (3.10) 2! 4! 6! Let 6 = 6 Q + \ Y t 2 , ( 6 - 6 Q ) ^ y t 2 , 6 ' - \ Y t 2 where Y ,(0) M Put 6 - 6 Q = 6 ' = ^  y t and substitute i n equation (3.10) to obtain P<°>6. P ( 4 ) ,(6) P ( 6 ' ) - p ( 0 ) + ± : + 'J- bl + ± 6 ,3 (3.11) 6Y 18 Y 6 ' represents the change i n rotor angle corresponding to changes i n power i . e . a t t = 0 , 6 ' = 0 E v a l u a t i n g e q u a t i o n (3.5) to obtain the energy function form 6 ^ to 35 6^, we obtain the energy function, p ( l ) 6 , 2 p ( 2 ) 6 , 3 p ( 3 ) 6 , 4 V = P ( } 6' + .1 + _? + — (3.12) a 2 3! 4! Equation (3.12) represents a catastrophe cusp equation, see Table (3.1). The f i r s t d e rivative of the energy function (3.12) that represents the catastrophe manifold (equilibrium surface) of a l l operating points of the power system i s given by V „ v = P ( ° > + P ( 1 ) 6' + P <2>iZ + ! f (3.13) 6' a a a 2! 3! To put equation (3.12) i n the cusp catastrophe manifold, the second order term must be eliminated. Let 6' = X - a P <2> and a = a (3.14) p (3) a Substitute equation (3.14) i n (3.13) to obtain 2 p (2) p O) V V = (P ( 0 ) + K ) + (P ( 1 ) - — ) X + JL X 3 - 0 X v a a' a , 2P ( 3 ) 6 a p ( 2 ) 2 = X 3 + (P ( 1 ) - X + -J— ( P ^ 0 ) + K a) - 0 p (3) a 2 p (3) p (3) a a a 36 X 3 + a X + b = 0 (3.15) where p ( 2 ) 2 a = _ i _ ( P CD -If ) p O) 3 2 p (3) a a b - _1_ (p a<°> + K „ ) p (3) 8 a P (2) P O) a3 " * 2! 3! Equation (3.15) i s exactly the cusp catastrophe manifold (see Appendix B). The steady state s t a b i l i t y region i s calculated from equation (3.15) by projecting i t down onto the cont r o l space (a,b), where: V 2 V = 3 X 2 + a = 0 (3.16) A. From equation (3.16) and equation (3.15) the state variable X can now be eliminated to get the b i f u r c a t i o n set B, 4 a 3 + 27 b 2 = 0 (3.17) Equation (3.17) represents the steady state s t a b i l i t y l i m i t s of the power system i n terms of the control parameters a and b as shown i n F i g . 3.3. 37 F i g . 3.3 The b i f u r c a t i o n set of the cusp catastrophe which represents the steady state s t a b i l i t y region. 38 3.3.2 Numerical Results Consider the one-machine i n f i n i t e - b u s system of F i g . 3.1, given the following data H = 3 s , V r o = 1. /0.0 pu E = 1.71 pu , 6 q = 36.5 pu X, = 1.05 pu , X_ = 0.36 pu d t Following the procedure i n the previous section, the sol u t i o n of equation (3.15) w i l l give the new operating angle f o r any change i n the input power. The steady state s t a b i l i t y region i s shown In F i g . 3.3. For any loading condition, the lo c a t i o n of the control parameters (a,b) on the b i f u r c a t i o n set of F i g . 3.3 w i l l determine the steady state s t a b i l i t y of the system. The r e s u l t s obtained i n F i g . 3.3 agree with the power equation (3.4). I t should be clear that the change i n the loading condition must be small and within the determining region of the Taylor series order used. However, a h i g h e r order T a y l o r s e r i e s can be used with a s u i t a b l e h i g h e r catastrophe. This method can be extended further to Include the speed governor e f f e c t . 3.4 A p p l i c a t i o n to the Transient S t a b i l i t y Problem The a p p l i c a t i o n of the catastrophe theory to the steady state s t a b i l i t y problem i s f a i r l y straightforward, as shown i n the previous section. In the case of transient s t a b i l i t y , the s i t u a t i o n i s d i f f e r e n t because there are two switchings ( d i s c o n t i n u i t i e s ) during the transient 39 > > . 3.4 The energy function stability c r i t e r i a a. stable b. c r i t i c a l l y stable c« unstable 40 period, one at f a u l t occurrence and the other at f a u l t clearance. Before we attempt to apply the catastrophe theory we need to f i n d a continuous function that represents the system behaviour during the transient period. We need also to define the degenerate and non-degenerate c r i t i c a l points i n terms of transient s t a b i l i t y . Consider again the one-machine i n f i n i t e - b u s system of F i g . 3.1. If a f a u l t occurs on one of the l i n e s near the machine bus, the rotor w i l l accelerate and gain k i n e t i c energy. I f the f a u l t i s cleared at the c r i t i c a l c l e a r i n g time, the k i n e t i c energy generated by the f a u l t w i l l be absorbed by the system and the gained energy at the end of the transient period w i l l be exactly zero; the system i s considered to be c r i t i c a l l y stable. A t y p i c a l energy curve for a f a u l t cleared at d i f f e r e n t c l e a r i n g times i s shown i n F i g . 3.4. In terms of catastrophe theory, curve (b) of F i g . 3.4 can be considered as the energy equilibrium surface or catastrophe manifold at which the k i n e t i c energy equals the p o t e n t i a l energy of the system. A l l non-degenerate c r i t i c a l points l i e on the energy equilibrium surface which corresponds to c r i t i c a l c l e a r i n g times. The degenerate c r i t i c a l points which are defined by the b i f u r c a t i o n set, are the points which correspond to the transient s t a b i l i t y l i m i t s of the power system at which any small disturbance of the power system w i l l drive the system unstable regardless of the c l e a r i n g time. In summary, since the energy function during the transient period i s continuous and represents the power system behaviour, the 'energy balance equation w i l l represent the equilibrium surface that decides the c r i t i c a l 41 c l e a r i n g time and the transient s t a b i l i t y l i m i t s define the degenerate c r i t i c a l points. 3»4.i. Single-Machine Infinite-bus Power System [42] Consider the power system of F i g . 3.5 which consists of one machine, an i n f i n i t e bus and two transmission l i n e s ab and dc. Representing the synchronous machine by a constant voltage source behind a reactance (the c l a s s i c a l model [2]), the swing equation representing the system behaviour i s given by .2, M — — = P. - P - P (3.18) , i e a d t z where P = P s i n 6' Is the e l e c t r i c a l power output, e max M = i n e r t i a constant of the machine. . > P^ = mechanical input power (constant during t r a n s i e n t ) . P^ = acc e l e r a t i n g power. 6 = rotor angle. Consider a three-phase f a u l t on l i n e ab cleared at the c r i t i c a l c l e a r i n g t i m e t . The a c c e l e r a t i n g power P w i l l e x h i b i t two ° c a d i s c o n t i n u i t i e s , one when the f a u l t occurs and the other when the f a u l t i s c l e a r e d . M u l t i p l y equation (3.18) by 6 and integrate with respect to time using the po s t - f a u l t network conditions to obtain 1 M 6 2 = P cos 6 + P, 6 - P cos 6 - P 6 (3.19) 2 c m c l e m m i m ? where 42 6^ = c r i t i c a l c l e a r i n g angle. 6^ = w£ = speed at c r i t i c a l c l e a r i n g . = maximum power of po s t - f a u l t network. 6^ = unstable equilibrium angle (maximum angle). The L.H.S. of equation (3.19) represents the k i n e t i c energy (K.E) generated during the fault-on period and the R.H.S. represents the po t e n t i a l energy (P.E) of the post - f a u l t network. I f the f a u l t i s cleared at the c r i t i c a l c l e a r i n g time, then K.E = P.E or K.E - P.E = 0 (3.20) and f o r s t a b i l i t y K.E < P.E (3.21) Thus equation (3.19) represents the equilibrium surface for the transient s t a b i l i t y or i n terms of catastrophe theory i t represents the catastrophe manifold, N, which i s the gradient of the energy i n t e g r a l function V (6^) . N = V. V(6 ) = 1 M 6 2 - P cos 6 - P, 6 + P cos 6 + P, 6 = 0 o c ' _ c m c i c m m i m c 2 (3.22) We define also the s i n g u l a r i t y set, S, as the set of transient s t a b i l i t y l i m i t s of the power system as V 2 V(6 ) = 0 (3.23) o c c Using Taylor series expansion to approximate 6 and 6 as a function of time c c (which has been reported to be a very good approximation for the f i r s t swing transient s t a b i l i t y analysis [43]) we get 43 6 = W = Yt and 6 = 6 +1 yt 1 (3.24) c c c c < > 2 c where Y = acceleration at the instant of f a u l t occurrence Y = I [P. - P (t ) ] (3.25) M Let X = I Y t 2 (3.26) 2 C and K = p 6 + p cos 6 m m m m S u b s t i t u t i n g f o r 6^ and 6^ by equations (3.24) and (3.26), replacing the cosine term (cos 6^) by i t s expansion, then equation (3.22) becomes (6 + X ) 2 (6 + X ) 4 N = V R V(6 ) = MYX - P [ l - — - + — ° ] - P f 6 + X ) + K = 0 6 c c m a . 21 4! 1 ° (3.27) Truncating the cosine expansion to the fourth order we get (6 + X ) 2 (6 + X ) 4 N = MYX - P +—2 P - — P - P , ( 6 + X ) + K = 0 (3.28) ' ma 2 ! ma ^, m a i o ' Since equation (3.28) Is a function of c l e a r i n g time which i s usually small, thus the fourth-order truncation would be enough to represent the equilibrium surface N. 44 N = - J l i x 4 - J ^ 6 x 3 + [—2-) P X 2 + (MY + 6 P - J± 6 3 - P J X 24 6 ° 4 o ma 6 o i ' + ( J 5 i 6 2 - 6 4 - P - P fi + K) = 0 2 0 24 0 1 1 1 3 P P Let A = , B = _S. 6 24 6 ° 2-6 2 C = ( °_) P , D = (MY + 6 P - 6 J - P ) v y m a oma t o i 4 6 and E = J U i 6 2 - 6 4 - P - P 6 + K) 2 o 24 ° 1 1 1 3 N = -AX4 - BX3 + CX2 + DX + E = 0 (3.29) Equation (3.29) is structurally stable and in the form of the swallowtail catastrophe (Appendix A) except we need to eliminate the cubic term. Let X = y - a N = -A(y-a) 4 - B(y-a) 3 + C(y-a) 2 + D(y-a) + E = 0 = -Ay4 + (4A<x - B)y 3 + (C-6Aa2 + 3Ba)y2 + (D - 2ca - 3Ba2 + 4Aa3)y - Aa 4 + Bot3 + Cot2 - Dot + E = 0 I B then a = 4 A Substitute for ot and get 45 8 A 2 B 8 .2 A 4 A 16 A2 256 A 3 divide by -A 4 rC 3 B 2 i 2 rD 1 CB 1 B^ \ N = y - I- . Jy - I- Jy A 8 .2 A 2 .2 8.3 Let u = v = 2 2 4 _ E + D * _ - £ _ _ - _1_ 5_ = 0 1 4 A 2 1 6 A 3 2 5 6 A 4 - £ + !•-) A 8 A 2 - ( E - I c ! - - ! ^ ) (3.30) A 2 . 2 8.3 A A 2 4 E . D B C B 3 B and w = - — + — 4 A 2 1 6 A 3 2 5 6 A 4 then N becomes N = y 4 + uy 4 + vy + w = 0 (3.31) Equation (3.31) i s the manifold of the swallowtail catastrophe 46 1 5 1 3 1 2 V(y) = ± y + - uy J + ± vy^ + wy (3.32) 5 3 2 We need to f i n d the s i n g u l a r i t y set, S, which i s the subset of N that consists of a l l the degenerate c r i t i c a l points of V. These are the points at which V V(y) = M = 0 y and V y 2 V(y) = 0 (3.33) V y 2 V(y) = 4y 3 + 2uy + u = 0 (3.34) I t i s i n t e r e s t i n g to note that the control variable (U) i n t h i s case i s constant and negative. 3 B 2 We have u = - (— + ) A 8.2 A P P i. » 1 1 1 3 T> ma c where A = , B = o 24 6 ° and C = P (-ma 2 - 6 2 thus, u = - (12 - 66 2 + 66 2 ) = - 12 o o This r e s u l t reduces the b i f u r c a t i o n manifold from three dimensions to only two dimensions i n u and W. We can plot the b i f u r c a t i o n set of the F i g . 3.6 The transient s t a b i l i t y l i m i t s given by the swallowtail catastrophe for the power system of F i g . 3 .5 48 power system using equations (3.31) and (3.34) (see Appendix A for d e t a i l s ) . The boundaries of the b i f u r c a t i o n set of F i g . 3.6 represent the degenerate transient s t a b i l i t y l i m i t s of the power system. I t should be noted that for a generator the control variables u and W are p o s i t i v e and for a motor they are negative. The region of transient s t a b i l i t y of the power system has to be i n s i d e the p o s i t i v e side of F i g . 3.6 and the following conditions have to be met: 1 2 Since X = _ Yt > 0 2 c where t i s the c r i t i c a l c l e a r i n g time and X = y - a then y - a > 0 (3.35) or y > a (3.36) I B c we have a = = o 4 A hence, y > 6^  (3.37) For the s t a b i l i t y region bounded by equation (3.31), (3.34) and the c o n s t r a i n t s y > 6 Q equation (3.31) has four c r i t i c a l points, two maxima and two minima. The c r i t i c a l c l e a r i n g time i s represented by the f i r s t p o s i t i v e c r i t i c a l point of y > 6 q. 3.4.2 Results The transient s t a b i l i t y region of the system shown i n F i g . 3.5 i s defined by equations (3.31), (3.34) and (3.37). Changing the f a u l t 49 50 F i g . 3.8 The transient s t a b i l i t y region i n terms of s t a b i l i t y l i m i t s and c r i t i c a l c l e a r i n g times. 51 l o c a t i o n , and f o r each f a u l t the loading condition (P^) i s changed from .1 to 1.2 pu. The shaded areas of Figure 3.7 represent the region of transient s t a b i l i t y for a l l possible f a u l t locations and loading conditions. A more comprehensive view of the s t a b i l i t y region can be v i s u a l i z e d by p l o t t i n g the region i n three dimensions where y i s the s i n g u l a r i t y (y > a > 0) representing the c r i t i c a l c l e a r i n g time T c - / 2 ( y " C ) (3.38) Y Figure 3.8 shows the transient s t a b i l i t y region i n three dimensions v. W and t . The t h i r d dimension t i s c a l c u l a t e d from the catastrophe ' c c r manifold and p l o t t e d over the shaded arc of Figure 3.7. This q u a l i t a t i v e view of the s t a b i l i t y region cannot be achieved by other d i r e c t methods of transient s t a b i l i t y , because these other methods use state variables to define the region of s t a b i l i t y which depends on operating conditions, f a u l t locations and time. But by using the catastrophe theory, the s t a b i l i t y region i s defined by two c o n t r o l v a r i a b l e s , u and W, only and includes a l l possible f a u l t locations and loading conditions. The control variables u and W can be written as a function of operating conditions and f a u l t l o c a t i o n s by using equations (3.30) and s u b s t i t u t i n g for A, B, C, D and E to get u = 2L (My + P j + 8(6 3 - 6 ) (3-39) p rma w = 24(6 2 + til 6 + * L 6 - JL + 1) o p o p o p ma ma ma In f a c t the r e g i o n of s t a b i l i t y defined by u, W and T £ i s fixed and 52 the machine i s stable when i t i s operating Inside the region. These r e s u l t s are I d e n t i c a l to that of the time s o l u t i o n . 3.4.3 Equivalent Two Machine Power System [44] Consider the power system of F i g . 3.9 which consists of three power plants A, B and C. A three-phase f a u l t i s applied on l i n e 5-6 near bus 6. I t i s found that machines B and C are swinging coherently against machine A. In t h i s case B and C can be combined together to form a large machine B+C, and the system can be reduced to a two-machine power system. The swing equation of machine A against the equivalent machine (B+C) can be written i n the form ,2c M - — = P - (P + P s i n (6 - a)) (3.40) 7 I c m d t Z The same procedure of the previous section can be applied to c a l c u l a t e the equilibrium surface from the energy balance equation [44]. S i m i l a r transient s t a b i l i t y regions of F i g s . 3.7 and 3.8 are also obtained i n t h i s case. The c r i t i c a l c l e a r i n g time calculated from the catastrophe manifold equation for t h i s case i s .348 seconds and that calculated by time solution i s .35 seconds. The method shows very close agreement with the numerical i n t e g r a t i o n method. When either the f a u l t l o c a t i o n or the loading condi-tions are changed, the machines may respond to the disturbance i n d i f f e r e n t coherency. This dictates d i f f e r e n t equivalence c a l c u l a t i o n s for d i f f e r e n t f a u l t locations and loading conditions. Therefore, a faster method F i g . 3.9 The equivalent two machine power system 5 4 i s needed to assess the response of the system to d i f f e r e n t contingencies. This problem w i l l be addressed i n the next chapter. 3 . 5 Summary The preliminary applications of catastrophe theory to the transient s t a b i l i t y problems of simple power systems have shown three a t t r a c t i v e advantages: 1. The transient s t a b i l i t y region i s well defined i n terms of the control v a r i a b l e s (system parameters) regardless of the state v a r i a b l e s . The s t a b i l i t y regions are bounded by the s t a b i l i t y l i m i t s which provide very good insi g h t to the security boundaries. Adequate s t a b i l i t y controls can be also designed and applied according to the transient s t a b i l i t y l i m i t s . 2. Extremely few computations are needed to define the s t a b i l i t y region, each transient case can be e a s i l y calculated i n terms of the system parameters without r e p e t i t i v e i t e r a t i o n s . 3 . The method i s a serious candidate to be used for on-line assessment of transient s t a b i l i t y , s e c u r i t y analysis and for the a p p l i c a t i o n of transient s t a b i l i t y controls such as dynamic breaking, fast valving, etc. 55 CHAPTER 4 TRANSIENT STABILITY REGIONS OF MULTIMACHINE POWER SYSTEMS 4.1 Introduction The transient s t a b i l i t y problem of multi-machine power systems i s much more complicated than the simple power systems analyzed i n Chapter 3. The analysis i n t h i s case involves every machine i n the power system without equivalencing any machine. The swing equation of each machine depends upon the movements of the other machines coupled to i t . There have been some major d i f f i c u l t i e s associated with the a p p l i c a t i o n of d i r e c t methods to the transient s t a b i l i t y problem of multi-machine power systems. Although noticeable progress has been reported i n dealing with these d i f f i c u l t i e s [45-48], the challenge of fast on-line transient s t a b i l i t y assessment i s s t i l l not overcome. These challenging problems are [49]: i . Power system model: applications of d i r e c t methods are l i m i t e d to the c l a s s i c a l model of power systems. In t h i s model, generators are represented by constant voltage magnitude behind transient reactance, loads are represented by constant impedances, transfer conductances are usually neglected or approximated, although they are s i g n i f i c a n t network elements [50]. i i . Although f a s t e xciters respond within the f i r s t swing period (.1 second), none of the e x i s t i n g d i r e c t methods consider e x c i t a t i o n e f f e c t s . 56 i i i . In p r a c t i c e , power engineers are usually not interested i n the c r i t i c a l c l e a r i n g time (which i s the main goal of d i r e c t methods) but rather i n the amount of power that can be delivered without r i s k i n g system security for s p e c i f i e d c l e a r i n g times. i v . Network reduction and c a l c u l a t i o n s of stable and unstable equilibrium points for large power systems are time consuming. For unstable equilibrium points, e x i s t i n g algorithms cannot always be r e l i e d upon to converge to the r i g h t s o l u t i o n . Two methods w i l l be introduced i n t h i s chapter to provide an energy function suitable for the a p p l i c a t i o n of catastrophe theory, namely the c r i t i c a l machines dynamic equivalent method and a method using Taylor series expansion of the accelerating power during the transient period. Both methods need the i d e n t i f i c a t i o n of the c r i t i c a l machines for each disturbance considerd. The transient s t a b i l i t y regions are then defined i n terms of the system parameters which enable power engineers to define the security regions and apply preventative c o n t r o l methods. The d i f f i c u l t i e s mentioned above w i l l be dealt with i n the proposed methods. Transfer conductance w i l l be f u l l y Included; damping and e x c i t a t i o n response w i l l be also included and discussed i n Chapter 5. 4.2 Dynamic Equivalent of the C r i t i c a l Machines The equation of motion of machine i i n a multi-machinfe power system, using c l a s s i c a l model representation i s given by 57 mi - P e i (4.1) n where P e i (4.2) = e l e c t r i c a l power output of machines Pmi = mechanical power input = i n t e r t i a constant 6 = rotor angle = speed deviation = i n t e r n a l voltage g^j = transfer conductance b . = transfer susceptance When a f a u l t occurs i n a large power system only a few machines ttct.iv?.i-.y reBpon.se to the f a u l t and tend to lose synchronism. These machines are known as the c r i t i c a l machines [45,47]. Therefore, i t i s enough to study the behaviour of the c r i t i c a l machines with respect to the rest of the power system i n order to evaluate the transient s t a b i l i t y of the system for a s p e c i f i c f a u l t . Consider that machine k i s the c r i t i c a l machine(s) for a s p e c i f i c disturbance. This machine i s considered to be o s c i l l a t i n g against the rest of the power system which i s not s i g n i f i c a n t l y affected by the disturbance and w i l l be considered moving as one machine n Let M, 0 (4-3) 58 1 n  U M Q i * k 1 1 Mg and 6 Q a r e , r e s p e c t i v e l y , the i n e r t i a constant and the angle of the center of angle of the power system excluding the c r i t i c a l machine. Let e k = 6 k - 6 0 - 6, - i _ I ( M l ) K M Q i * k 1 1 also k * M Q i*k 1 9 . - 6 . - 1 - S (M 6 ) (4.5) k * M Q i * k 1 1 We substitute equation (4.1) into (4.5) to obtain the swing equation of machine k with respect to the center of angle (COA). Pmk " Pek 1 n \~ JV^-r- * ( P m i - P e i ) <4'6> \ M o i , t k V k " Pmk ' Pek " ^ X <Pmi " P e i > <4'7> M Q i* k S u b s t i t u t i n g f o r P ^ and P g^ from equation (4.2) and separating 9^ from the rest of the system we get the swing equation of the c r i t i c a l machine k against the rest of the power system. This i s explained i n the following steps: > 59 F i r s t we l e t <= E± g ±^ and C^-^y^ Equation (4.7) becomes "kV Pmk " Dkk - X < D i j C O S °kj + Scj S i n V M, n n - JL Z (P - Z (D. cos 0 + C, , s i n 9 ) j M Q i * k m i J - l l j i j l j l j Separate the term B^ k from the l a s t term and obtain Vk - Pmk " Dkk " ^ X ^ Pmi " j k < D i j C O S 0 i j + C i j S i n G i j ) ) K n + — ^ (D cos 9 + C s i n 9 ) M Q i*k ± k i k l k l k - z^ ( D k j cos e k j + s i n e k j ) Let P = P - D,, - ^ Z (P . - S (D cos 9 + C s i n 9 )) (4.8) k mk kk M i # k mi i j I j l j i j then W = P k ~ ^ b S l n °k " 3 C O S \ ] ( 4 * 9 ) where KL n n a = _ S (D„ cos 9 + C., s i n 9 ) - Z (C s i n 9 - D cos 9 ) M n i*k l k 1 ± k 1 j*k J J 6 0 n M, n b = £ (D • s i n 9 + C cos 0 ) - _J± Z (D., s i n 9 - C,. cos 9 ) j*k k 3 2 M Q i * k l f c i i k i Equation ( 4 . 9 ) can be written i n a more convenient form \ \ - pk - T k s i n <\ - V ( 4 a o ) where - 1 a a = tan _ k b . 2 _,_ v 2 x l / 2 and T, = (a + b ) k Equation (4 .10) i s a simple form representing the motion of the c r i t i c a l machine for a c e r t a i n disturbance. Since we assumed that the rest of the system i s not responding to the disturbance, i t i s reasonable to use the pre-disturbance angles © q to ca l c u l a t e the parameters P^, and c t ^ . The stable and unstable equilibrium points of equation (4 .10) can be e a s i l y computed by solving equation (4 .11) for 9^ P k " T k S i n ( Q k " V = ° ( 4 - 1 D and the unstable equilibrium point (UEP) Is 9 U = n - ef (4 .12) k k We note here that we have two sets of the parameters P^, T r and a^; one set for the fault-on network and another for the post-fault network. 61 4.2.1 The Transient S t a b i l i t y Region During the transient period, an exchange of energy takes place between the rotor of the c r i t i c a l machine and the p o s t - f a u l t network. The k i n e t i c energy generated by the accelerating power during the fault-on period must be f u l l y absorbed by the p o s t - f a u l t network i n order to maintain s t a b i l i t y . • C O M u l t i p l y i n g e q u a t i o n (4.10) by and integra t i n g between 0^ and with respect to time for the fault-on network parameters we obtain the k i n e t i c energy generated by the f a u l t • 2 k - e = \ \Ql = p k ( e k • e k ) - T k [ c o s ( e k - 4) - c o s ( e k • ° k ) ] (4.13) f f f C where P^, T^ and oc^  are the system parameters for fault-on network and 0^ i s the c l e a r i n g angle. The p o t e n t i a l energy of the p o s t - f a u l t network i s derived i n the same manner by integra t i n g equation (4.10) between and 0^ using the p o s t - f a u l t network parameters, we obtain • 2 - I M ^ = ?l (9j - e£) - TP[cos(e^ - o^P) - cos(Q^ - of)] (4.14) Note that the speed deviation at the unstable equilibrium point i s zero (e£ = o ) . 6 2 Equation ( 4 . 1 4 ) represents the energy function of the c r i t i c a l machine during the transient period. The L.H.S. represents the f a u l t k i n e t i c energy and the R.H.S. represents the p o t e n t i a l energy for the post - f a u l t network. The energy balance equation for c r i t i c a l c l e a r i n g becomes: • 2 I M - pP e£ - T£ cos (Qp - o£) + k U = 0 . . ( 4 . 1 5 ) where k U = P £ 9^  + T£ cos (GJJ - c£) r again we represent 9^  by Taylor ser i e s expansion and ek = V c where Y, = — [ P , - P^Ct-. ] k w 1 k ek v OH " 1 . 2 L e t X = ~ V c By r e p l a c i n g cos (9^ - ot^) i n equation (4.15) with the cosine ser i e s expansion up to the fourth order, we obtain M Y,x - pP(9° + X) - T J[1 - 1 L_+ 1 L_ + k U = 0 2 ! 4! Let P - 9° - aP then we get the catastrophe manifold equation 63 k + i p x 3 + i ( i - 5l) x 2 + (MY, - pP) 2 + Tl B - i f33 24 6 2 2 k k k 6 )x + ( k u - T P - p P e ° + ! k p 2 - ^ p 4 ) = 0 (4.16) 2 24 24 M u l t i p l y equation (4.16) by - to give: T P f 2 T X 4 - 4BX 3 - 12(1 - *L)X2 - I 4 (MY, - pP + B - A p 3 )X 2 T p k k k 6 k ~ ( k U " T k - P k 9 k + ~ P2 " ~ - 0 T p 2 24 Ak (4.17) We need to eliminate the t h i r d order term i n order to have equation (4.17) i n the swallowtail catastrophe manifold form Let A 3 = - 4 B A 2 = - 12(1 2 A x = - 2± (MYk - P j + TP B - ± B 3 ) 4 64 A Q - - 1* ( k U - - 0P + Tl B 2 - Tl B 4) T P 2 24 Xk and x = y - u (4.18) A3 u = _ (4.19) 4 We obtain the swallowtail catastrophe manifold 4 2 y +uy + u y + w = 0 (4.20) 3 2 where u = (A„ - — A ) 8 » . (A, - ^ I t f L ) 1 2 8 2 / A A 1 A 3 . A 2 A 3 3 . 4. w = (A - + - A- ) U 4 16 256 The b i f u r c a t i o n set B can then be defined by 4y 3 + 2uy + u = 0 ... (4.21) The transient s t a b i l i t y region i s formed i n the shape of the swallowtail b i f u r c a t i o n set (see Appendix A) bounded by the transient s t a b i l i t y l i m i t s . The region i s defined i n terms of the system parameters which can be e a s i l y c a l c u l a t e d . In summary, the procedure of c a l c u l a t i n g the transient s t a b i l i t y regions for multi-machine power systems using the critical-machines dynamic-equivalent method i s outlined i n the following steps: 6 5 1. Define the c r i t i c a l machine(s) for the f a u l t considered. 2. Calculate the parameters Pfc, T f c and to form the dynamic equivalent of the c r i t i c a l machine(s). 3. Calculate the catastrophe manifold parameters u, u and w from which the c r i t i c a l c l e a r i n g time and the degree of s t a b i l i t y are determined. I t i s important to note that the transient s t a b i l i t y region using t h i s method is the same shape regardless of f a u l t l o c a t i o n and loading condition. This makes the method a good candidate for on-line assessment of the transient s t a b i l i t y and the degree of security of power systems. 4.3 Taylor Expansion of the Accelerating Power The chaotic nature of the swing equations of multi-machine power systems i s a major d i f f i c u l t y i n solving the transient problem and i n t r y i n g to include more de t a i l e d system models In d i r e c t methods. Rec a l l i n g the swing equation of machine i i n an n-machine power system from section (4.2): V i - Pmi < D i j C O S 6 i j + C i j S l n 6 i j > = P f l i ( i = 1, ... n) (4.22) where P ^ i s the accelerating power of machine i . In the conventional time s o l u t i o n , equations (4.22) are solved by n u m e r i c a l i n t e g r a t i o n for the transient period ( t n -*• t ), from the instant 66 of f a u l t occurrence to the unstable equilibrium point. I t i s already known that each po s t - f a u l t network has a fi x e d energy-absorbing capacity, c a l l e d the c r i t i c a l energy [51]. The c r i t i c a l energy can be calculated o f f - l i n e for any post-fault condition. This energy i s the numerical i n t e g r a t i o n of the p o s t - f a u l t accelerating power from the stable equilibrium point to the unstable equilibrium point. We suggest using Taylor series expansion to f i n d the energy function as a function of time by expanding the accelerating power equation around the instant of f a u l t occurrence at t=0. This provides the advantage of elimin a t i n g the cal c u l a t i o n s of the stabel and unstable equilibrium points on-line and s i m p l i f i e s the procedure of defining the c r i t i c a l machines. The energy function of machine i i n a multi-mchine power system during the transient period i s given by: t t t c f c u V i = f P a l ( t ) d t + f P a i ( t ) d t + ' P a i ( t ) d t ( 4 ' 2 3 ) where P ^ i s fault-on a c c e l e r a t i n g power. P f ( t ) i s po s t - f a u l t accelerating power, a i t i s the c l e a r i n g time. t ,t i s the times at SEP and UEP. s u The f i r s t two terms of the R.H.S. of equation (4.23) are the k i n e t i c energy generated during fault-on, and the t h i r d term i s the c r i t i c a l energy of the post - f a u l t network. The portion of the k i n e t i c energy needed to move 6 7 the rotor from 6r t to 6 of the po s t - f a u l t network does not contribute to the 0 s r i n s t a b i l i t y [51]. Therefore, i t should be subtracted from the energy f u n c t i o n i n order to write the energy balance equation; t h i s portion i s t s / P n J ( t ) dt (4.24) 0 a l Subtracting (4.24) from (4.23), and assuming that machine i Is the c r i t i c a l machine and the f a u l t i s cleared at c r i t i c a l c l e a r i n g , we obtain the energy balance equation, t c t u V. = / (pf ( t ) - p j ( t ) ) d t + / PJJ(t) dt = 0 1 Q a a t a 8 (4.25) Eq u a t i o n (4.25) i s evaluated by replacing P (t) and P P ( t ) with t h e i r 3. 3. Taylor series expressions, i n the form: P (t) = P ( 0 ) + i - P ( 1 ) t + i - P ( 2 > t 2 + (4.26) a a u a 2 , a where P ( m ) = 1 a d t m (4.27) t=0 and the accelerating power of machine i n P - P - E (D . cos 6 + C.. s i n 6 ) a i mi i j i j i j i j We note that at t=0 = 0 68 and Therefore 6^ = 6^ ( I n i t i a l operating condition) P , ( 0 ) = P - P . (t = 0) a i mi e i = P mi n j = \ < D i j c ° s 6 i j + c i j s i n (4.28) (1) _ d a i dt 3 1 = Z (D.. s i n 5 ° - C.. cos 6.1)6.. t - 0 J-l i j i j i J i J i j Since i j = 0 t-0 Therefore P f 1 * = 0 a i (4.29) *£2i - -a l d t a i n t=0 ^ ( D ± j cos 6 ± ° + C ± j s i n 6 ^ ) 6 ^ + £ ( D i j s i n 6±i Since 6 = 0 n P ! 2 ) = S (D s i n 6 ±° a i J-l j C i j C O S 6 i j ) 6 i j n c 0 c 0 C i j C ° S 6 i j 6 i j (4.30) where 6 . ^ = p ( i n i t i a l a cceleration) I j a l I The derivations of the rest of the Taylor s e r i e s c o e f f i c i e n t s are given i n Appendix C. I t i s found that a l l odd order c o e f f i c i e n t s are zeros; 69 t h i s reduces the number of Taylor series c o e f f i c i e n t s to one h a l f . Therefore, Pal = P a 0 ) + L p i 2 > < 2 + M 4 ) ' 4 + - M - 'D d i <* 2! 4! The accelerating power c o e f f i c i e n t s have to be calculated for the system during the f a u l t and for the post f a u l t condition. Usually i n large power systems the c r i t i c a l c l e a r i n g time i s very small ( t y p i c a l l y < .5 second). For t h i s period we w i l l use up to the s i x t h order of the acc e l e r a t i n g power expansion for both fault-on and post-fault conditions. Thus equation (4.25) becomes; C c £ D P ( 2 ) f - P ( 2 ) P P ( 4 ) f - P ( 4 ) P V . - J ( < P l 0 ) - P < 0 ) ) + 0 ^ _ ) t 2 • £ t 4 1 0 • a 2! 4! p ( 6 ) f _ p ( 6 ) P + (_f 1 ) t 6 ) d t + / P P ( t ) dt - 0 (4.32) 6! t a s The l a s t term i s the c r i t i c a l energy of each s p e c i f i c post-fault condition. The c r i t i c a l energy i s evaluated o f f - l i n e for each post-fault condition by the trapezoidal rule as follows t v c r - J U P j ( t ) dt s V (1+1) = V (k) +1 (P (t)(k+l) - P A ( t ) ( k ) ) - (t(k+l) - t ( k ) ) 70 4.3.1 Transient Stability Regions Using Taylor Expansion The Taylor series expansion of the accelerating power provides flexible choice over the seven elementary catastrophes of Table (3.1). The choice of a specific catastrophe type depends on the highest order of the Taylor series that w i l l give the best required results. It is needless to add that the lower order catastrophes are easier to visualize. We start with the resultant Taylor series expansion of equation (4.32). P( 2 > p ( 4> . p (t) = p (o) + — t z + - i _ t* + ... a a 2! 4! The rotor angle can also be represented by the form 6 = 6 + I Y t 2 I Oi 2 T i which provides good agreement with the time solution for t < .5 second [52]. * c 1 2 We let x = i " °0i = - Y i fc 2 substitute for t in the accelerating power series p(2) p(4) P (x) = P ( 0 ) + _!_ X + _ ! _ X 2 + ... (4.33) a a Y 2 Y i 6y±z The energy balance equation can be evaluated with respect to the angle advances x x < c V = J P (x)dx + V =0 i Q a 71 P ( 3 ) p ( 4 ) 3 = P ( 0 ) x + — x 2 + _ i _ J L + . . . + V = 0 ( 4 . 3 4 ) 3 2Y 2 C r 18 The choice of the determinant order decides the catastrophe manifold type as follows X 2 = f o l d V 3 X = cusp 4 = swallowtail X 5 = b u t t e r f l y A Here we consider the cusp catastrophe manifold given i n Appendix B. The cusp catastrophe manifold from equation ( 4 . 2 4 ) i s given by p<4> P < 2 > v = a x3 + a x2 + p ( 0 ) x + V = 0 1 8 Y 2 " i 3 » T T P « ' 2 1 8 Y 2 P < 0 ) 1 8 r 2 = xJ + — _ — x z + a X + IV = 0 ( 4 . 3 5 ) P ( 4 > p ( A ) P ( 4 ) ^ a a a = X 3 + ^ X 2 + A^X + AQ = 0 ( 4 . 3 6 ) In order to put equation ( 4 . 3 6 ) i n the standard cusp manifold we need to 2 eliminate the second order term (X ). We l e t X = y - B A P = — We get 72 3 A2 2 2 A 1 A 2 y J + (A - — — ) y + (A + — — - -i_£) = 0 1 3 27 3 (4.37) or y + u y + u = 0 (4.38) where u u Equation (4.38) i s i n the standard form of the cusp catastrophe. The b i f u r c a t i o n set i s defined by equation (4.38) and the set of degenerate c r i t i c a l points. The b i f u r c a t i o n set given i n Appendix B consists of a l l stable points i n terms of the system parameters u and u. 4.4 I d e n t i f i c a t i o n of the C r i t i c a l Machines Both methods presented i n Sections (4.2) and (4.3) r e l y on the a c c u r a t e i d e n t i f i c a t i o n of the c r i t i c a l machines f o r a s p e c i f i e d disturbance. Correct i d e n t i f i c a t i o n has been achieved by c a l c u l a t i n g the unstable e q u i l i b r i u m points for a l l machines i n the power systems; the machine which 3y 2 + u - 0 (4.29) 73 has the highest unstable equilibrium point i s i d e n t i f i e d as the c r i t i c a l machine [47], Although this method provides correct i d e n t i f i c a t i o n , i t s main drawback i s the c a l c u l a t i o n of the UEP's, which i s time consuming, and i n some cases, wrong answers may be obtained [48]. Another method i s to use the ac c e l e r a t i o n at the instant of f a u l t occurrence as a f i r s t i d e n t i f i c a t i o n and then c a l c u l a t e the c r i t i c a l c l e a r i n g time for each machine. The machine with the lowest c r i t i c a l c l e a r i n g time i s the c r i t i c a l machine [45]. For small and medium s i z e power systems, the c a l c u l a t i o n of UEP's i s not d i f f i c u l t and usually the correct answer i s achieved i n a reasonable time, but for large power systems the r i g h t answer i s not guaranteed. In t h i s t h e s i s , the c r i t i c a l machines are i d e n t i f i e d as follows: I. Calculate the i n i t i a l a c c e l e r a t i o n for each machine as follows \ " ^ t Pmi " P e i ( t > M i where Pg^Ctg^ *"s t ^ i e e ^ e c t r i c a i power output d u r i n g f a u l t at the i n s t a n t of f a u l t occurrence. i i . The machines which have high and p o s i t i v e i n i t i a l accelerations are i n j e c t i n g k i n e t i c energy to the system; therefore, they a l l contribute to the system I n s t a b i l i t y . These machines are combined to form a c r i t i c a l group. i i i . The c r i t i c a l energies of the c r i t i c a l machines are also added together to form the global energy-absorbing capacity of the c r i t i c a l group. 74 In most cases c r i t i c a l machines can be I d e n t i f i e d e a s i l y by the f i r s t a c c e l e r a t i o n , e s p e c i a l l y when the f a u l t i s close to one of the generator terminals. The procedure of combining a group of c r i t i c a l machines has to be c a r r i e d out only when a f a u l t occurs at non-generator buses or far from generator buses. 4.5 Numerical Examples In t h i s section two examples are presented to demonstrate the v a l i d i t y and advantages of the a p p l i c a t i o n of catastrophe theory to transient s t a b i l i t y assessment of power B y s t e m s . Transient s t a b i l i t y regions i n terms of the system parameters are given for each example. Three and seven machine power systems are used, three-phase short c i r c u i t s are considered at d i f f e r e n t locations and a comparison between the presented methods and the step-by-step time s o l u t i o n i s given for each example. Each s h o r t - c i r c u i t case considered i s evaluated by the following steps. i . I d e n tify the c r i t i c a l machine or machines as explained i n Section 4.4. i i . Dynamic equivalent method; c a l c u l a t e the dynamic equivalent parameters of equation (4.10) for the c r i t i c a l machine(s) (as given i n Section 4.2). The b i f u r c a t i o n set parameters and c r i t i c a l c l e a r i n g time are then calc u l a t e d . i i i . Taylor series method; c a l c u l a t e the Taylor series c o e f f i c i e n t s of the accel e r a t i n g power for the c r i t i c a l machine (as given i n Section 4.3). The b i f u r c a t i o n set parameters and the c r i t i c a l c l e a r i n g time are then calculated and the case i s located on the transient s t a b i l i t y region. 75 4.5.1 The Three-Machine System This system with nine buses, three machines and three loads, i s widely referred to i n the l i t e r a t u r e as the Western Systems Coordinating C o u n c i l (WSCC) test system. A s i n g l e - l i n e diagram of the system i s shown i n F i g . 4.1. Transmission l i n e parameters and loads are given i n per unit on a 100-MVA base i n Table 4.1. Generator data and i n i t i a l operating conditions are given i n Table 4.2. Bus No. Admittance G (pu) B Generators 1 1-4 0 - 8.446 2 2-7 0 - 5.485 3 3-9 0 - 4.168 Transmission 4-5 1.365 -11.604 ... • Lines 4-6 1.942 -10.511 5-7 1.188 - 5.975 6-9 1.282 - 5.588 7-8 1.617 -13.698 8-9 1.155 - 9.784 Shunt Admittances Load A 5-0 1.261 - 0.263 Load B 6-0 0.878 - 0.035 Loda C 8-0 0.969 - 0.160 4-0 0.167 7-0 0.227 9-0 0.283 Table 4.1 Network parameters of the three-machine system 77 Generator data I n i t i a l Conditions Gen. No. H P E d mo (MW/MVA) (pu) (pu) (pu) (degree) 1 23.64 .0608 2.269 1.096 6.95 2 6.40 .1198 1.6 1.102 13.49 3 3.01 .1813 1.0 1.024 8.21 Table 4.2 Generator data and i n i t i a l conditions of the three-machine system Three-phase short c i r c u i t s are considered at d i f f e r e n t l o c a t i o n s . The transient s t a b i l i t y i s evaluated for each f a u l t by the step-by-step time s o l u t i o n , dynamic equivalent and Taylor series methods. A comparison between the three methods i s given i n Table 4.3 i n terms of the c r i t i c a l c l e a r i n g time. Both methods show very good agreement with the time s o l u t i o n . F i g . 4.2 shows the accuracy of the Taylor series method using terms up to the second order term for the angle during the fault-on period. F i g . 4.3 shows a comparison between the Taylor series approximation of the a c c e l e r a t i n g power and the time s o l u t i o n for d i f f e r e n t f a u l t locations during the transient period. The transient s t a b i l i t y region using the dynamic equivalent method i s shown i n F i g . 4.4. A l l stable cases are shown inside the region i n terms of the system parameters. F i g . 4.5 show the transient s t a b i l i t y region using the Taylor series method. Stable points are marked inside the region. 250.0 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.5 time i n seconds F i g . 4.2 Accuracy of Taylor s e r i e s during the fault-on period using only the second order term. time s o l u t i o n Taylor s e r i e s method C.C.T. c r i t i c a l c l e a r i n g time 0.3 time i n seconds - A f f e r e n t f a u l t locations F i e . 4 .3 The ac c e l e r a t i n g power for d i f f e r e n t ra time s o l u t i o n Taylor se r i e s Fault at C r i t i c a l Clearing Time [second) Bus # Time Solution (step-by-step) Dynamic equivalent method Taylor Series method 4 .29-.30 .27 .29 5 .22-.23 .20 .22 6 .51-.|2 .46 .49 7 .11-.12 .1 .11 8 .28-.29 .26 .28 9 .40-.41 .38 .41 Table 4.3 C r i t i c a l c l e a r i n g time by the time so l u t i o n and the proposed methods for the 3-machine power system. 8 1 F i g . 4.4 The transient s t a b i l i t y region of the 3-machine power system using the dynamic equivalent method. Stable cases are marked (x) ins i d e the region. F i g . 4.5 The transient s t a b i l i t y region of the 3-machine system using Taylor s e r i e s method with stable cases marked (x) inside the region. 83 84 4.5.2 CIGRE 7-machine Test System The CIGRE 225 KV test system i s shown i n F i g . 4.6. I t has 10 buses and 13 l i n e s . Buses 1 through 7 are generating buses while loads are located at buses 2, 4, 6, 7, 8, 9 and 10. Generator, load and transmission l i n e data are given i n Table 4.4. The base values used are 225 KV and 100 MVA. Three-phase f a u l t s are applied and the transient s t a b i l i t y i s evaluated for each f a u l t . The c r i t i c a l c l e a r i n g time i s calculated for each case using the three methods. Table 4.5 gives a comparison between the time s o l u t i o n and the two methods presented i n t h i s chapter i n terms of the c r i t i c a l c l e a r i n g time. The transient s t a b l i t y regions are shown i n Figures 4.7 and 4.8. Both methods show very good agreement with the time s o l u t i o n plus well defined transient s t a b i l i t y regions i n terms of the system parameters v a l i d for a l l loading conditions and f a u l t l o c a t i o n s . 4.6 Discussion of Results The r e s u l t s of the two examples presented i n the previous section have shown the f e a s i b i l i t y of c a t a s t r o p h e theory a p p l i c a t i o n to multi-machine power systems. The accuracy of the methods i s very good when compared with the step-by-step time s o l u t i o n which i s used as the benchmark for accuracy. The methods proposed are fast and give the c r i t i c a l c l e a r i n g time through a single computation instead of using the time so l u t i o n method r e p e t i t i v e l y and human i n t e r p r e t a t i o n to evaluate a s i m i l a r case. 85 GENERATORS Bus p base X M pm E 6° (MVA) (%)( ) (MW s /rad) (MW) (p.u) (°) 1 100 7.4 6.02 217 1.106 7.9 2 100 11.8 4.11 120 1.156 -0.2 3 100 6.2 7.59 256 1.098 6.5 4 100 4.9 9.54 300 1.110 3.9 5 100 7.4 6.02 230 1.118 7.0 6 100 7.1 6.77 160 1.039 3.6 7 100 8.7 5.68 174 1.054 7.9 LOADS Bus p Q Bus P Q (MW) (MVar) (MW) (MVar) 2 200 120 8 100 50 4 650 405 9 230 140 6 80 30 10 90 45 7 90 40 LINES Bus R X uC/2 (ohm) (ohm) (MS) 1 - 3 5 24.5 200 1 - 4 5 24.5 100 2 - 3 22.8 62.6 200 2 - 10 8.3 32.3 300 3 - 4 6 39.5 300 3 - 9 5.8 28 200 4 - 5 2 10 200 4 - 6 3.8 10 1200 4 ~ 9 24.7 97 200 4 - 10 8.3 33 300 6 - 8 9-5 31.8 200 7 - 8 6 39.5 300 8 - 9 24.7 97 200 ( ) These values include the transformer's reactances and are expressed on a 100 MVA base. Table 4.4 Data for the CIGRE 7-machine system (taken from Ref. [20]) Fault at C r i t i c a l Clearing Time [second) Rus # Time Solution (step-by-step) Dynamic equivalent method Taylor Series method 1 .35-.36 .34 .36 2 .41-.42 .40 .37 3 .39-.40 .38 .39 4 .50-.51 .52 .48 5 .35-.36 .34 .35 6 .52-.53 .51 .51 7 .33-.34 .33 .33 8 .44-.45 .42 .45 Table 4.5 C r i t i c a l c l e a r i n g time by the time so l u t i o n and the proposed methods for the 7-machine power system. 87 88 89 The Importance of the transient s t a b i l i t y regions provided by the b i f u r c a t i o n set of the catastrophe manifold i s not only i n terms of speed and accuracy. I t provides another important dimension to the transient s t a b i l i t y problem-the s t a b i l i t y l i m i t s . The method checks f o r the v i o l a t i o n s of the transient s t a b i l i t y l i m i t s of the post-fault network i n terms of the system parameters. Of course, i f the l i m i t s are exceeded the system i s unstable. This important feature i s not possible with e x i s t i n g d i r e c t methods, i . e . , the e x i s t i n g d i r e c t methods cannot give any solu t i o n when the s t a b i l i t y l i m i t s of the post-fault network are v i o l a t e d . The accuracy of the Taylor series method becomes a problem when t i s higher than .5 second [52]. To obtain good accuracy beyond t h i s l i m i t higher order terms have to be included i n the computations. This w i l l complicate the s t a b i l i t y region and slow down the c a l c u l a t i o n procedure. In pr a c t i c e , however, t h i s problem i s very rare. Faults i n large power systems are usually tripped i n a few cycles, t y p i c a l l y .1 second. Interest i s i n the delivery of maximum power at low c l e a r i n g time without r i s k i n g the power system s e c u r i t y . This p r a c t i c a l consideration can e a s i l y be handled by the proposed methods without loss of accuracy. The o f f - l i n e c a l c u l a t i o n s of the c r i t i c a l energy for each case considered can be improved by c a l c u l a t i n g the closest unstable equilibrium point for the c r i t i c a l machine. This reduces the computation time and s i m p l i f i e s the procedure of f i n d i n g the c r i t i c a l energy [45]. 9 0 CHAPTER 5 INCLUSION OF DAMPING, FLUX DECAY AND EXCITATION RESPONSE The most challenging problem in the application of direct methods to transient stability of power systems is the validity of the simplified model used to represent power systems. In this Chapter a significant improvement in the modelling problem is presented. The new model proposed includes damping, f i e l d flux decay and excitation response. 5.1 Limitation of the Classical Model Power u t i l i t i e s are hesitant to accept the direct methods of transient stability assessment mainly because they raise doubts about the validity of the classical model. In this model i t is assumed that the flux linking the main f i e l d winding remains constant during the transient. This may be true only i f the exciter does not respond during the f i r s t swing period (1 second or less). Modern excitation systems can reach f u l l response within .1 second, so that the classical model assumption is not valid for such exciters. In fact, during the last decade trends in the design of power system components have resulted in more reliance on fast-response and high ceiling voltage exciters [53]. Although fast-response exciters are more desirable and widely used in modern power systems, none of the existing direct methods of transient s t a b i l i t y assessment have considered the excitation effect. The reason is 91 that Lyapunov's function becomes so complicated that the d i r e c t methods lose the merit of being fast methods. In c l a s s i c a l model assumptions, damping e f f e c t s are usually neglected although i t a f f e c t s d i r e c t l y the accelerating power. In some cases the damping power has appreciable value and so a f f e c t s the f i r s t swing s t a b i l i t y [54]. 5.2 Damping In power systems there are sources of p o s i t i v e damping which tend to damp out o s c i l l a t i o n s r e s u l t i n g from disturbances. This damping i s due to the c h a r a c t e r i s t i c s of the mechanical system, generator and loads. The mechanical system damping r e s u l t s from the increase i n shaft torque with the decrease i n speed. The per unit damping torque c o e f f i c i e n t i s defined as the negative of the per unit change i n torque for each per unit change i n speed [54]. The generator e l e c t r i c a l damping i s associated with the currents that are induced i n the rotor windings due to a disturbance or o s c i l l a t i o n s . The damping torque due to the f i e l d winding i s usually small because of the r e l a t i v e l y long time constant. The damper winding component of the damping torque i s quite appreciable and i t usually a f f e c t s accelerating power during disturbances. C y l i n d r i c a l rotor generators may develop appreciable torque due to the eddy currents which are induced by asynchronous o p e r l t i o n , such as s l i p p i n g poles or o s c i l l a t i o n s . 92 Loads are also a source of damping i n power systems. Induction and synchronous motors with t h e i r mechanical shaft loads develop damping torque during disturbances. The damping power (D6) i s usually added to the i n e r t i a l power i n the swing equation. The damping c o e f f i c i e n t D includes the various damping power components, both mechanical and e l e c t r i c a l . The damping c o e f f i c i e n t s are usually i n the range of 1-3 pu. This represents mechanical damping, generator and load damping. Larger values are also reported i n the l i t e r a t u r e [2]. The swing equation of machine i In an n-machine power system i s given by M. 6 + D 6 = P - P . . .. i i i mi e i i = l , . . . n (5.1) where D i s the damping c o e f f i c i e n t . The accelerating power becomes P a i - Pmi " P e i " D * i ( 5 ' 2 ) The same procedure of S e c t i o n 4.3, i s used to expand P &^ i n t o a Taylor seri e s to get the energy function during the fault-on period. R e c a l l i n g e q u a t i o n s (4.26) and (4.27) P 4 ( t ) = P ( 0 ) + ! - p ( 1 ) t + ! - p ( 2 > t 2 + . . . p(m) = a i a , m dt t=0 we have 93 P = P , - Z (D, . cos 6 + C. . sin 6 ) - D& ai mi j = 1 i j i j i j i j i ( 5 . 3 ) We note that at t = 0 , 6^ = 0 . Therefore ai ( 0 ) = P_, " E (D*4 cos 6 ° + c4 . sin 6 ° ) - D(0) " m i j=l i J ' i j " i j ' i j ' ( 5 . 4 ) (1) _ d P a i ai dt t=0 1 1 0 - - E (-D4, sin 6,, + C,, cos 6 v 6 i j i j - i j i j i j D)6f> = - D 6|°) ( 5 . 5 ) (0) where 6 (0 ) _ ai M. ( 5 . 6 ) (2 ) _ ^ a i ai dt' t=0 n 0. ,? 0 ,2 - = ( " D i J C ° S ^ J ' C i J S ± n ' i ? ( 6 i > ^ ( - D ± j sin 6 ± ° + cos 6 ± ° ) 6^ - 0*6^ - E (D sin 6 ° - C cos 6 ° ) 6 ° - D*6^0 » ( 5 . 7 ) 94 P ( 1 ) where 6 0 = — (5.8) M i D e r i v a t i o n of p a ^ ^ a n& p a i ^ a r e 8 i v e n i n Appendix D. We note here that a l l Taylor series c o e f f i c i e n t s e x i s t when damping i s taken into account. P, = P $ 0 > + P < 1 > t + ! ? l ^ t 2 + ... (5.8) a i a i a i 2J The energy equilibrium surface i s then formed using equation (4.23) V i = ' Q ( P a i " P a i > d t + I P a ? d t + V c r = 0 <5'9> This becomes = A 4 t 4 + ^ t 3 + A ^ 2 + A^t + A Q + k g + V £ r - 0 (5.10) whe re A • - i - (P f ( n ) - P P ( n ) ) n . a n t c K = / P P dt - constant s J Q a i Here i s t r u n c a t e d a f t e r the 4th order term. We eliminate the 3rd order term to obtain the standard swallowtail catastrophe manifold (Appendix A) as given i n Section 3.4.1 to get X 4 + uX 2 + vX + w - 0 (5.11) The transient s t a b i l i t y region i s defined through of the b i f u r c a t i o n set i n 'y terms of the system parameters u, v and w which includes the damping term. 95 The e f f e c t of damping on the accelerating power during the transient period i s shown i n F i g . 5.1 for a short c i r c u i t on the power system of F i g . 3.5. Another example i s given i n F i g . 5.2 for a short c i r c u i t at bus 7 of the three-machine power system given i n F i g . 4.1. I t i s clear that the damping tends to reduce the f a u l t k i n e t i c energy and r e s u l t s i n more stable systems. 5.3 E x c i t a t i o n Response and Flux Decay When a f a u l t occurs near a generator bus, transient currents i n the armature c i r c u i t induce other currents i n the rotor c i r c u i t which c a r r i e s the f i e l d current. The f l u x l i n k i n g the armature c i r c u i t w i l l decay according to the e f f e c t i v e time constant of the f i e l d c i r c u i t . This time constant i s i n order of several seconds at no load, and one second or higher under load. The f l u x decay decreases the generator i n t e r n a l voltage and hence reduces the s t a b i l i t y l i m i t s . Both steady-state and transient s t a b i l i t y can be improved by e x c i t a t i o n control systems [55]. Fast speed of response and high c e i l i n g voltage exciters can p a r t i c u l a r l y improve transient s t a b i l i t y . With the help of f a s t transient f o r c i n g of e x c i t a t i o n and the boost of i n t e r n a l machine f l u x , the output power of the generator can be increased during the transient period. This reduces the accelerating power and r e s u l t s i n improved transient performance. Fast exciters also ensure that subsequent swings are smaller than the f i r s t swing. This i s important it or modern low i n e r t i a generators and weakly damped power systems where transient i n s t a b i l i t y a f t e r the f i r s t swing may occur. 96 0i8_. ___ . o.o H 1 1 1 1 1 1 ' 1 — ^ 1 ' 1 • — 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 time in seconds Fig. 5.1 0.8 T I 0.0-i 1 1 " 1 " 1 < 1 1 1 • 1 r J — 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 time in seconds Fig. 5.2 Effect of damping inclusion on the accelerating power during the fault-on period for two test cases. 9 7 The vector diagram for the transient state of a c y l i n d r i c a l rotor generator i s given i n F i g . 5.3. We assume that the i n t e r n a l voltage lags behind the q-axis by a constant angle <t>^  a l l the time. The equations of motion of generator i i n an n-machine power system with f i e l d f l u x decay and e x c i t a t i o n response are expressed as follows n M i 6 i = Pmi " j = \ < D i j C 0 S 6 i j + C i j 6 l n 6 i J > and dE' T'j < — — " EfA< ~ E % - ~ x ^ ) 1*4 (5.12) doi ^ f d i q i d i d l ' d i where T^ Q^ : d-axis transient open c i r c u i t time constant E q i : v o ^ - t a 8 e a l ° n 8 the q-axis related to E^ by angle 4>^ . E f d i : e x c i t a t * - o n voltage applied to f i e l d winding. X,.,X' : d-axis synchronous and transient reactance, d i d i I,. : d-axis current, d i The terms of the swing equation were previously defined. From the vector diagram of F i g . 5.3, and since <t>^  i s constant E . = E cos <t> q i i i Equation (5.12) becomes • d E i Efdi Tdoi — - — r " \ - c xdi - V d i . ( 5 a 3 ) dt cos 4> Equation (5.13) can be modified i n terms of the power equation and assuming that 4>^  i s small we obtain [56]: 98 F i g . 5.3 The vector diagram of the generator i n t e r n a l quantities for the transient state. 99 T d o i — - — - - v ( x d i - x d i > , V g i J E J c o s 6 i J + b i J E J s i n 6 i J ) ( 5 - 1 4 ) dt c o s ^ j=l J J J J J The f i r s t term of the R.H.S. of equation (5.14) represents the exciter response while the second and t h i r d terms represent the drop i n the i n t e r n a l v o l t a g e due to the f l u x decay. I t i s clear that i f E ^ ^ i s fa s t and high t h e n the drop of i n t e r n a l v o l t a g e (E^) due to the f l u x decay can be eliminated. A t y p i c a l e x c i t a t i o n voltage-time response curve i s shown i n F i g . 5.4. OA i s the generator rated load f i e l d voltage before disturbance occurs. The s t r a i g h t l i n e AC i s drawn such that the area ACD i s equal to the area ABD enclosed by the actual response. The exc i t e r response i s given by the rate of Increase or decrease of the ex c i t e r voltage, i . e . , the slope CE *"* AO determined by i n F i g . 5.4. The time, i n seconds, for the ex c i t e r OE voltage to reach 95 percent of c e i l i n g voltage i s known as exc i t e r voltage response time [57]. An e x c i t e r with a voltage response time of 0.1 second or less i s taken to be a fa s t high i n i t i a l response e x c i t e r . The e x c i t e r voltage response i s given by - t / T E... = E + (E„ - E )(1 - e e ) (5.15) r d i o c o where E i s the i n i t i a l e x c i t e r voltage o E maximum c e i l i n g voltage c T e x c i t a t i o n system time constant Substituting by E... i n equation (5.14) we get r d i E x c i t e r Voltage A time i n seconds 5.* The t J P i d - c i t a t i o n , o l t a g e - t i * e response 101 , d E, E - (E -E )e dt cos <t>. n " E i " ( X d i " X d i ) \ i s i l E l C O S 6 i j + b i j E j S l n 6 i j } ( 5 * 1 6 ) j The function F ( t ) w i l l now be expanded i n a Taylor series to evaluate equation (5.16) for the transient period (t < .5s). F ( t ) i s expanded around t=0 as follows F ( t ) - F ( 0 ) + F a > + ^ 2 ) t 2 + _ t=0 2! (5.17) where ,(n) _ d u F ( t ) dt t=0 iTherefore E n F ( 0 ) . _o_ _ £ i ( 0 ) _ ( x ^ _ + b ^ E j B l n f i ^ ) cos <t>, (5.18) F ( D . d F(t) dt t=0 102 -tli 1 <V Eo ) e ° d E i n T cos 4^ dt M j=l J J J + b i j E j cos 6 ^ ) 6 ^ at t = 0 i j dt T 1 doi Thus . (E -E ) (0) F ( 1 ) = l (5.19) \ C ° S * 1 T d o i p ( 2 ) ^ F ( 3 ) ^ f ( n ) ^ e v a l u a t e d i n the same way. Equation (5.16) i s evaluated using (5.17), and we get E. = E(o) + F ( 0 ) t + I — + I — (5.20) 1 2 6 where E(o) i s the i n i t i a l generator i n t e r n a l voltage at t=to-. Equation (5.20) represents the generator i n t e r n a l voltage during the transient period including e x c i t e r response and f l u x decay. Truncation of E^ depends on the l e n g t h of the transient period. As we explained i n the previous chapter the Taylor series method can give good r e s u l t s up to a period of .5 seconds. In p r a c t i c e , f a u l t s are usually cleared f a s t e r than the .5 s period (a t y p i c a l c l e a r i n g time i n large power ^  systems i s .1 second). 103 Equation ( 5 . 2 0 ) can be used with the method presented i n Section 4 . 3 to include e x c i t a t i o n response and f l u x decay i n the assessment of transient s t a b i l i t y . The swing equation of the c r i t i c a l machine i i n an n-machine system i s given by n V i = Pmi~ f x E i E j < S i j c o s 6 i j + b i j s i n 6±j> ( 5 > 2 1 )  = P a i From equation ( 5 . 2 0 ) l e t the i n t e r n a l voltage of the c r i t i c a l machine be ( 0 ) F ( 1 > t 2 F ( 2 ) t 3 E = E(o) + F t U ; t + _ — + £ _ ( 5 . 2 2 ) 1 2 6 S u b s t i t u t e ( 5 . 2 2 ) i n ( 5 . 2 1 ) and expand P i n a Taylor series as i n Section 4 . 3 to obtain P ( 2 ) t 2 P 4 • P ( ° > + P < 1 > t + a ... a i a a 21 where P ( n ) = *A d t n t = ° The accelerating power c o e f f i c i e n t s are evaluated and given as follows P a 0 ) - P m i ~ ^ E I ( 0 ) E j ( g ± j cos 6t° + b ± j s i n 6±°) ( 5 . 2 3 ) p ( D = - I F < ° > E . (g., cos 6 ° + b., s i n 6 ° ) ( 5 . 2 4 ) a j = 1 3 i j i j i j IJ 104 P a 2> - - ^ F U > E j (8lJ cos 6^ + b ± J s i n 6^) + ^ E ± ( 0 ) E j ( g ± j s i n 6±° - b ± j cos 6 ^ ) ^ ° (5-25) The d e r i v a t i o n of these c o e f f i c i e n t s and P ^ 3 \ P^4^ are given i n Appendix 3 3 E. The energy equilibrium surface (catastrophe manifold) can now be written down i n terms of the accelerating power c o e f f i c i e n t s as i n equation x (5.9) V i = 0' ( P a i " P a ? ) d t + ' Q P a ? d t + V c r = 0 < 5" 2 6> Equation (5.26) i s evaluated and truncated a f t e r the 4th order term to get the swallowtail catastrophe manifold. This truncation i s v a l i d for t < .5 second., Thus equation (5.26) becomes V. = A . t 4 + a 0 t 2 + A.t + A_ + K + V = 0 (5.27) I H I 1 0 s cr where A = — (P f ( n ) - P P ( n ) ) n i a i a i n! t c K = / P P dt = constant s J Q ao V = c r i t i c a l energy cr 105 Again we eliminate the t h i r d order term of equation (5.27) to get the standard swallowtail catastrophe (Appendix A) i n the form X 4 + uX 2 + vX + w = 0 (5.28) The e f f e c t of f l u x decay and e x c i t e r response on the accelerating power during the transient period i s shown i n F i g . 5.4 and F i g . 5.5 using examples shown i n F i g . 3.5 and F i g . 4.1. I t i s apparent that the area under the accelerating power i s reduced. This area represents the k i n e t i c energy generated by the f a u l t . Therefore, the degree of s t a b i l i t y i s increased and the c r i t i c a l c l e a r i n g time i s also increased. This means that when fast e x c i t e r s are used the system can generate more power for the same cl e a r i n g time. 0.8-0.6-0.4 0.2-0.0 c l a s s i c a l model proposed method 0.00 0.05 —I ' 1 0.10 0.15 time i n seconds F i g . 5.5 0.20 0.25 0.8-0.6-0.4-0.2 0.0 c l a s s i c a l model proposed method 1 0.00 — I — 0.05 0.20 0.10 0.15 time i n seconds F i g . 5.6 The effect of flux decay and excitation response on the accelerating power during the transient period for two test cases. 0.25 107 CHAPTER 6 CONCLUSIONS In t h i s t h e s i s , the transient energy function was e f f i c i e n t l y used with catastrophe theory to define comprehensive transient s t a b i l i t y regions for power systems. E x p l i c i t boundaries of the transient s t a b i l i t y regions were i d e n t i f i e d by the b i f u r c a t i o n set of the catastrophe manifolds. The method i s made p r a c t i c a l by using the concept of c r i t i c a l machines since for a s p e c i f i e d f a u l t , only a few machines show s i g n i f i c a n t o s c i l l a t i o n s during the transient period. The thesis suggests two methods to solve the transient s t a b i l i t y problem of multi-machine power systems. In the f i r s t method, the c r i t i c a l machine(s) i s i d e n t i f i e d and a two-machine dynamic equivalent for the power system i s formed. The c r i t i c a l machine i s singled out and a center of angle group i s formed for the rest of the system. The energy balance equation i s derived from the equation of motion of the c r i t i c a l machine against the center of angle. At c r i t i c a l c l e a r i n g , the energy balance equation forms the equilibirum surface of the catastrophe manifold from which the transient s t a b i l i t y region i s derived by the b i f u r c a t i o n technique. This s t a b i l i t y region i s v a l i d for d i f f e r e n t loading conditions and f a u l t l o c a t i o n s . This method has the advantage that the c r i t i c a l energy i s calculated d i r e c t l y by using the closest unstable equilibrium point. However, this method can only be used with the c l a s s i c a l model. * In the second method, the c r i t i c a l machine i s i d e n t i f i e d and a Taylor s e r i e s expansion for i t s a c c e l e r a t i n g power for a s p e c i f i e d f a u l t i s 108 constructed. The energy balance equation i s then derived by integrating the ac c e l e r a t i n g power. The c r i t i c a l energy i s calculated o f f - l i n e by numerical i n t e g r a t i o n for the f a u l t considered. Again the energy balance equation at c r i t i c a l c l e a r i n g represents the catastrophe manifold from which the trans-ient s t a b i l i t y region i s determined. This region i s bounded by the s t a b i l i -ty l i m i t s and i s v a l i d for d i f f e r e n t loading conditions and f a u l t locations. This method i s l i m i t e d to a period of .5 second for best accuracy. However, t h i s i s adequate for most large power systems since t y p i c a l c l e a r i n g times are i n the range of a few cycles or .1 second. Model Improvements were made possible using t h i s method. E x c i t a t i o n response, f l u x decay, and damping, were included i n the s t a b i l i t y analysis by the Taylor series expansion method. The r e s u l t s obtained by the proposed methods are i n good agreement with those obtained by the time s o l u t i o n method. A large number of simula-tions are presently needed i n system planning i n order to determine the c r i t i c a l c l e a r i n g time or s t a b i l i t y l i m i t s . The proposed methods give d i r e c t l y the c r i t i c a l c l e a r i n g time and s t a b i l i t y l i m i t s with good accuracy and less computation. Therefore, they can be used to greatly reduce the large number of time simulations i n the system planning stage. In system operations, the proposed methods provide f a s t solutions to the transient s t a b i l i t y problem with d e f i n i t e s t a b i l i t y boundaries so that corrective a c t i o n can be taken to prevent i n s t a b i l i t y . The major contributions of t h i s research are the following: 1. For the f i r s t time, the catastrophe theory i s applied to the transient s t a b i l i t y problem. Comprehensive transient s t a b i l i t y regions are 109 calculated i n terms of the power system parameters. These regions are v a l i d for any loading condition and f a u l t l o c a t i o n . A l l e x i s t i n g d i r e c t methods, on the other hand, require new computation for any chf.nge i n operating conditions. This i s an important consideration i n implementing d i r e c t methods for real-time assessment of transient s t a b i l i t y and system s e c u r i t y . Another unique advantage of the methods presented i s the i n c l u s i o n of the e x c i t a t i o n response during the transient period i n the s t a b i l i t y a n a l y s i s . A l l e x i s t i n g d i r e c t methods are l i m i t e d to the c l a s s i c a l model [47] which neglects the e x c i t e r response. With the i n c l u s i o n of damping, fl u x decay and e x c i t a t i o n response i n the presented methods, a l l factors that d i r e c t l y e f f e c t transient s t a b i l i t y are taken into account. This r e s u l t makes the d i r e c t methods more r e a l i s t i c . The methods presented predict any v i o l a t i o n of s t a b i l i t y l i m i t s and, therefore transient i n s t a b i l i t y . This r e s u l t i s very important; i t helps to ensure the security of power systems for any disturbance considered. I t also enables power system planners to design proper s t a b i l i t y controls to prevent system i n s t a b i l i t y . E x i s t i n g d i r e c t methods either neglect or approximate the e f f e c t of the transfe r conductances. However, i n some cases transfer conductances can have an appreciable e f f e c t on system performance [50] . In this research, the transfer conductances are f u l l y represented. The a p p l i c a t i o n of catastrophe theory i s extended to the steady state s t a b i l i t y problem of c y l i n d r i c a l - r o t o r generators. The method 110 presented (Chapter 3) i s also adaptable to include the e f f e c t of governors on steady state s t a b i l i t y . 6. The achievements of t h i s research have made the on-line security assessment more f e a s i b l e than ever before. The important requirements for real-time applications are speed, accuracy and information. The methods presented need minimal c a l c u l a t i o n s to locate the system para-meters with respect to the s t a b i l i t y region. Good agreement with the time so l u t i o n method were obtained (Chapter 4). The l o c a t i o n of the system parameters, for each case, on the transient s t a b i l i t y region gives enough information on the degree of s t a b i l i t y of the power system. The outcome of t h i s thesis gives new motivation to the a p p l i c a t i o n of f a s t d i r e c t methods to the transient s t a b i l i t y problem of power systems. New areas of research need to be explored i n order to reach the ultimate goal of the d i r e c t methods which i s on-line assessment of transient s t a b i l i -ty. The future research should include the following: 1. The c a l c u l a t i o n of the c r i t i c a l energy of the p o s t - f a u l t network i s s t i l l time consuming and needs to be done o f f - l i n e . With some approxi-mations, the dynamic equivalent method can give fast answers by calcu-l a t i n g the closest unstable equilibrium points. However, a fast and exact method i s s t i l l needed i n order to complete the o v e r a l l on-line approach. 2. S t a b i l i t y controls such as f a s t valving, braking r e s i s t o r s , single pole switchings, series capacitors and generator trippings are usually applied i n p r a c t i c e to restore transient s t a b i l i t y of power systems. I l l The i n c l u s i o n of these controls i n the proposed d i r e c t methods i s of great in t e r e s t to power u t i l i t i e s . This research considered only single disturbances. M u l t i p l e d i s t u r -bances should also be considered i n future research. More can also be done i n the area of steady state s t a b i l i t y . The presented method (Chapter 3) i s adaptable to include speed governor co n t r o l systems. This Is an i n t e r e s t i n g area to investigate and an on-line steady state s t a b i l i t y assessment approach can be developed using the catastrophe theory. 112 REFERENCES [1] E.W. Kimbark, "Power System S t a b i l i t y " , V o l . I, John Wiley & Sons Inc., New York, 1948. [2] P.M. Anderson and A.A. Fouad, "Power System Control and S t a b i l i t y " , V o l . I, Iowa State University Press, Ames, Iowa, 1977. [3] G.W. Stagg and A.H. El-Abiad, "Computer Methods i n Power System Analysis", McGraw-Hill, 1968. [4] C. Concordia, "Representation of Loads", IEEE PES Winter Meeting, paper 75 CH0970-4-PWR, 1975. [5] IEEE Committee Report, "Dynamic Load Dynamics Simulation E f f e c t s and Determination of Load Constants", Paper T73-089-0, IEEE Trans. PAS, Vol. 92, March/April 1973, pp. 600-609. [6] F. John Meyer and K.Y. Lee, "Improved Dynamic Load Model for Power System S t a b i l i t y Studies", IEEE PES, Winter Meeting 1982, 82 WM 126-1. [7] A.M. Sasson, "The Power System S t a b i l i t y Problems and i t s Solution Methods", IEEE PES Summer Meeting, Los Angeles, C a l i f . , July 1970. [8] W.D. Humpage and B. Stott, "Predictor-Corrector Methods of Numerical Integration i n D i g i t a l Computer Analysis of Power System Transient S t a b i l i t y " , Proc. IEE, V o l . 112, pp. 1557-1565, 1965. [9] N.D. Rao and H.N. Rao, "Solution of Transient S t a b i l i t y Problems through the Number Series Approach", Proc. IEE, Vol. I l l , pp. 775-788, 1964. [10] D.L. Johnson, J.B. Ward, "The Solution of Power System S t a b i l i t y Problems by Means of D i g i t a l Computers", Trans. AIEE-PAS, V o l . 75, pp. 1321-1329, 1957. [11] S.N. Talukdar, " I t e r a t i v e Multi-Step Methods for Transient S t a b i l i t y Studies", IEEE Paper No. 70TP 48-PWR, Winter Power Meeting, New York, 1970. [12] N.D. Rao and H.N. Rao, "Phase-Plane Techniques for the Solution of Transient S t a b i l i t y Problems", Proc. IEE, Vol. 110, 1963, pp. 1451-1461. [13] K.N. Stanton and S.N. Tulukdar, "New Integration 'Algorithms for Transient S t a b i l i t y Studies", IEEE Trans., V o l . PAS-89, May/June 1970, pp. 785-991. [14] H.W. Dommel and N. Sato, "Fast Transient S t a b i l i t y Solutions", IEEE Trans., Vol. PAS-91, July/August 1972, pp. 1643-1650. 113 [15] K. Prabhashanker and W. J a n i s c h e w s y j , " D i g i t a l S i m u l a t i o n o f Multi-Machine Power Systems for S t a b i l i t y Studies", IEEE Trans., V o l. PAS-87, Jan. 1968, pp. 73-80. [16] L.H. Michaels, "The AC/Hybrid Power System Simulator and i t s Role i n System Security", Proc. of the Power Industry Computer Applications Conference, Boston, Mass., May 1971, pp. 128-136. [17] A.H. El-Abiad and K. Naguppun, "Transient S t a b i l i t y Regions of Multi-Machine Power Systems", IEEE Trans., Vol. PAS-85, pp. 169-179. [18] J.L. Willems and J.C. Willems, "The Appl i c a t i o n of Lyapunov Methods to the Computations of Transient S t a b i l i t y Regions for Multi-Machine Power Systems", IEEE Trans. PAS., v o l . 89, No. 5/6, pp. 795-801, May/June 1970. [19] F.S. Prabhakara and A.H. El-Abiad, "A S i m p l i f i e d Determination of Transient S t a b i l i t y Regions for Lyapunov Methods", IEEE Trans. PAS, Vol.-94, No. 2, March/April 1975, pp. 672-689. [20] M.A. Pai, "Power System S t a b i l i t y " , V o l. 3, North Holland Publishing Co. [21] M. Ribbens-Pavella and F.J. Evans, "Direct Methods for Studying Dynamics of Large-Scale E l e c t r i c Power Systems - A Survey", Automatica, Vol. 21, No. 1, pp. 1-21, 1985. [22] T. Athay, R. Podmore and S. Vivmani, "A P r a c t i c a l Method f o r Direct Analysis of Transient S t a b i l i t y " , IEEE Trans., PAS-98, March/April 1979, pp. 573-584. [23] A.A. Fouad and S.E. Stanton, "Transient S t a b i l i t y of Multi-Machine Power System Part I and I I " , IEEE Trans. PAS-100, pp. 3808-3424, Aug. 1981. [24] M. Ribbens-Pavella, Th. Nan Cutsem, R. Ohifoui and B. Toumi, "Energy Type Lyapunov-Like Direct C r i t e r i a for Rapid Transient S t a b i l i t y A n a l y s i s " , Proceedings of International Symposium on Power Systems S t a b i l i t y , Ames, Iowa, May 1985. [25] A.A. Fouad, "Application of Transient Energy Functions to P r a c t i c a l Power System Problems", IEEE Speical P u b l i c a t i o n on Rapid Analysis of Transient S t a b i l i t y , Power Winter Meeting, 1987. [26] A. Valveer, et a l . , "Relative Transient S t a b i l i t y Power System Margin Index", IEEE PES Summer Meeting and EHV/UHV Conference, Paper C73 474-4, Vancouver, Canada, July 1973. [27] C K . Pang, F.S. Prabhakara, A.H. El-Abiad and A.J. Koivo, "Security Evaluation i n Power Systems Using Pattern Recognition", IEEE PES Winter Meeting, New York, Feb. 1973. 114 [28] K.S. Fu, "Sequential Method i n Pattern Recognition and Machine Learning", Academic Press, 1968. [29] W. Meisel, "Computer Oriented Approaches to Pattern Recognition", Academic Press, 1972. [30] K.S. Fu and J.M. Mendel, "Adaptive Learning and Pattern Recognition", Academic Press, 1970. [31] G. Hadley, "Linear Programming", Addison Wesley, 1962. [32] P.H. Mengert, "Solution of Linear I n e q u a l i t i e s " , IEEE Trans, on Computers, Vol. C-19, pp. 124-131, 1970. [33] C L . Gupta, "Transient Security Assessment of Power Systems by Pattern Recognition and Lyapunov's Direct Method", Ph.D. th e s i s , Purdue University, West Lafayette, Indiana, 1976. [34] H. Hakkimmashhadi, "Fast Transient S t a b i l i t y Assessment", Ph.D. thesis , Purdue University, West Lafayette, Indiana, Aug. 1982. [35] S. Yamashira, T. Koike, A.H. El-Abiod, "Fast Transient Security Assessment and Enhancement Using Pattern Recognition", Proc. of 8th PSCC, Aug. 1984. [36] T. Poston and I. Stewart, "Catastrophe Theory and I t s Applications", Pitman Publishing Co., London, 1979. [37] A.A. Sallam and J.L. Dineley, "Catastrophe theory as a Tool for Determining Synchronous Power System Dynamic S t a b i l i t y " , IEEE Trans. PAS, Vol. 102, pp. 622-630, March 1983. [38] R. Thorn, "St r u c t u r a l S t a b i l i t y and Morphogenesis", Benjamin-Addison Wesley, New York, 1975. [39] I. Stewart, "Elementary Catastrophe Theory", IEEE Trans, on C i r c u i t s and Systems, V o l . CAS-30, pp. 578-586, August 1983. [40] P.T. Saunders, "An Introduction to Catastrophe Theory", Cambridge University Press, 1980. [41] R.F. Chu, Discussion to Reference [37], IEEE Trans. PAS, Vol.-102, pp. 629-630, March 1983. [42] M.D. Wvong and A.M. M i h i r i g , "Catastrophe Theory Applied to Transient S t a b i l i t y Assessment of Power Systems", Proc. IEE, Vol. 133, Pt. C, No. 6, Sept. 1986. [43] M. Ribbens-Pavella, Thn Van Cutsem, R. Odifaoui and B. Toumi, "Energy Type Liapuniv-Like Direct C r i t e r i a for Rapid Transient S t a b i l i t y A nalysis", P r o c of the International Symposium on Power System S t a b i l i t y , Ames, Iowa, May 1985. 115 [44] A.M. M i h i r i g and M.D. Wvong, "Application of Catastrophe Theory to Transient S t a b i l i t y Analysis of Power Systems", Proceedings of LASTED International Conference on High Technology i n the Power Industry, Bosman, Montana, August 1986. [45] Y. Xue, Th. Van. Cutsem and Ribbens-Pauella, "A Simple Direct Method for Fast Transient S t a b i l i t y Assessment of Large Power Systems", IEEE PES Winter Meeting, New Orleans, Louisiana, February 1-6, 1987. [46] R. Rahimi and G. Schaffer, "Power System Transient S t a b i l i u t y Indexes for On-Line Analysis of 'Worst Case' dynamic Contingencies", IEEE PES Summer Meeting, Paper 86 SM332-1, Mexico Cit y , Mexico, July 20-25, 1986. [47] A.N. Michel, A.A. Fouad and V. V i t t a l , "Power System Transient S t a b i l i t y Using Individual Machine Energy Functions", IEEE Trans., CAS-30, May 1983, pp. 266-276. [48] N. Kakimoto, Y. Ohnogi, H. Matsuda and H. Shibuya, "Transient S t a b i l i t y Analysis of Large-Scale Power System by Lyapunov's Direct Method", IEEE/PES, 83SM 403-3, July, 1983. [49] A.A. Fouad, "Application of Transient Energy functions to P r a c t i c a l Power System Problems", IEEE Special P u b l i c a t i o n on Rapid Analysis of Transient S t a b i l i t y , No. 87 TH0169-3-PWR. [50] V.R. Sastry, " V a l i d i t y of Nelgecting Transfer Conductances i n Transient S t a b i l i t y Studies", P r o c IEE, Vol. 120, No. 12, December 197?, [51] A.A. Fouad and S.E. Stanton, "Transient S t a b i l i t y of A Multi-Machine Power System", Part I and I I , IEEE Trans., Vol. PAS-100, pp. 3408-3424, July 1981. [52] S. Furuya and S. Iwamoto, "Fast Transient S t a b i l i t y Solution Using Taylor Expansion and Energy Function", EE i n Japan, Vol. 105, No. 3, 1985. [53] C. Concordia and P.G. Brown, " E f f e c t s of Trends i n Large Steam Turbine Generator Parameters on Power System S t a b i l i t y " , IEEE Trans. PAS-90, pp. 993-1000, 1971. [54] S.B. Crary, "Power System S t a b i l i t y " , V o l. 2, John Wiley & Sons, Inc., New York, 1947. [55] IEEE Work Group on Special S t a b i l i t y Controls PSEC, ?k Bibliography on the App l i c a t i o n of Discrete Supplementary Controls to Improve Power System S t a b i l i t y " , IEEE Trans., Vol. PWRS-2, No. 2, May 1987. [56] N. Kakimoto, Y. Ohsawa and M. Hayashi, "Transient S t a b i l i t y Analysis of Multi-machine Power Systems with F i e l d Flux Decays Via Lyapunov's Direct Method", IEEE Trans., Vol. PAS-99, No. 5, Sept/Oct. 1980. 116 [57] M.S. Sarma, "Synchronous Machines (Their theory, S t a b i l i t y and E x c i t a t i o n Systems)", Gordon and Breach Science Publishers, New York. 117 APPENDIX A The swallowtail catastrophe: The p o t e n t i a l function i s V(X) = X + uX + vx + wX The equilibrium surface M i s the hyper surface 5X + 3uX + 2uX + w = 0 (A.l) and the s i n g u l a r i t y set i the subset of M for which the equation 20X + 6ux + 2v = 0 (A.2) The b i f u r c a t i o n set, B, i s a three-dimensional surface i n the control space u, u and w. Since we are concerned only with the q u a l i t a t i v e behaviour of the system, and therefore want p r i m a r i l y to be able to plot the b i f u r c a t i o n set B. Let C be a plane u = constant i n C (the control space), B w i l l be a u u surve i n C, and i f we can sketch t h i s curve for a l l values of u we can b u i l d up the complete surface B. Equation (A.2) implies that u i s an odd function of X, and that i s together with ( A . l ) , implies that w i s an even function X. Hence w i s an even f u n c t i o n of u, and so f o r any u the curve B u i s symmetric about the w-axis. Nest we d i f f e r e n t i a t e (A.l) and (A.2) obtaining dX dw (A.3) 118 and — = - (30X 2 + 3u) (A.4) dX the rest of the term i n (A.3) having vanished i n account of (A.2). We now have to consider the cases u > 0 and u < 0 separately. I f U i s p o s i t i v e , then — cannot v a n i s h . Hence u i s a s t r i c t l y dX monotone function of X and the equation ^ = - 2 (A.5) d u i s v a l i d everywhere. Moreover, equation (A.2) implies that X u < 0 with e q u a l i t y only when X = u = 0, at which point w also vanishes. I t f o l l o w s that B y i s smooth, that w i s large when |x| i s large, and dw t h a t the s i g n a l of — i s the same as that of u, v a n i s h i n g only at the d u o r i g i n . This enables us to draw F i g . ( A . l ) . I f u i s n e g a t i v e , then — vanishes for two r e a l values of X, ±/ ——• dX 10 dw Consequently, — vanishes f o r three values of X, these two together with dX X = 0 as b e f o r e , and i t f o l l o w s that B u has a c r i t i c a l point at X = 0 and cusps at the other two points. To determine the type of the c r i t i c a l point, we notice that Equation (A.2) implies that for |x| < / the product Xu cannot be negative. Since 10 X and u also vanish together, i t follows that i f u i s small and pos i t i v e so i s X, and — i s then n e g a t i v e . T h i s together w i t h the f a c t that B u i s d u 119 symmetric about the w-axis, establishes that the c r i t i c a l point i s a r e l a t i v e maximum. F i n a l l y , we note that i f u = 0 then either X = 0 or X = ±/ -1^ .. We 10 have just seen that X = 0 corresponds to a maximum at the o r i g i n , and su b s t i t u t i n g into equation (A.l) we f i n d that both the other roots give 9 u 2 w = . Hence B has a point of s e l f - i n t e r s e c t i o n on the p o s i t i v e w-axis. 20 u We then check that |x| large implies that both | u| and w are also large and then, using the values of the parameter X to t e l l us that the order of the p o i n t s we have found i s : s e l f - i n t e r s e c i t o n , cusp, maximum, cusp s e l f - i n t e r s e c t i o n , we can draw F i g . (A.2). And since the equation of the l i n e of points of s e l f - I n t e r s e c t i o n i s the parabola 2 I 20 v - * L . u « 0 We can put the curves B u together to form the surface B shown i n F i g . (A.3). The o r i g i n of the name "swallowtail" i s now apparent. To f i n d the form of the p o t e n t i a l i n each of the three regions into which B divides C, i t i s s u f f i c i e n t to consider points for which u = 0 and u < 0 then the so l u t i o n of Equation (A.l) i s X 2 = — (- 3u ± /(9u 2 - 20w)) 10 There are three cases: 2 9u (a) w > E q u a t i o n ( A . l ) has no r e a l r o o t s and V has no c r i t i c a l 20 120 Fig.(A.2 ) Fig.(A.1 ) Cross section of the b i f u r c a t i o n set for u > 0 and u < 0 r>f the swallowtail The b i f u r c a t i o n set of tne * 121 points. 9 u 2 2 0 < w < Because /(9u - 20w) i s r e a l and less than the r e a l and 20 2 p o s i t i v e - 3 u , both solutions for X are r e a l and p o s i t i v e and V has four c r i t i c a l points, two maxima and two minima. 2 w < 0 B o t h s o l u t i o n s f o r X a r e r e a l but one i s n e g a t i v e . Consequently V has only two c r i t i c a l points, one minimum and one maximum. 122 APPENDIX B The cusp catastrophe: The p o t e n t i a l function i s V(X) = X 4 + uX 2 + uX (B.l) so the equilibrium surface i s a three-dimensional space i n x, u and u given by 4X 3 + 2uX + u = 0 (B.2) and the s i n g u l a r i t y set i s the subset of the equilibrium surface such that the de r i v a t i v e of (B.2) i s also equal to zero. I t i s given by 12X 2 + 2u = 0 (B.3) We f i n d the b i f u r c a t i o n set by eliminating the state variable X from (B.2) and (B.3), we obtain 8U"3 + 27u 2 - 0 (B.4) Equation (B.4) i s the pr o j e c t i o n of the three-dimensional manifold of equation (B.2) onto the control space (u-u). The cusp manifold and the b i f u r c a t i o n set i s shown i n F i g . ( B . l ) . Equation (B.2) has three r e a l roots within the b i f u r c a t i o n set region, or when 8u 3 + 27 u 2 < 0 But when 8u 3 + 27u 2 > 0 123 F i g . B . l The cusp manifold and i t s b i f u r c a t i o n set. 124 there i s only one r e a l root. Figure (B.2) shows the b i f u r c a t i o n set of u-u plane i n which the fuctions V(X) i s sketched for d i f f e r e n t values of the parameters u and u. „ , . ».* Tne cusp p o t e n t i a l f - c t i o n VCX) at d i f f e r e n t values of the c o n t r o l v a r i a b l e s . 126 APPENDIX C Taylor series expansion of the accelerating power: The accelerating power of machine I i n an n-machine power system i s given by (Equation 4.22) P = P , - Z (D. . cos 6 + c. . s i n 6 ) a i mi ± - 1 i j i j i j i j ' ( C I ) We expand P ^ i n the neighborhood of t=0 to obtain, P , t ) - P<°> + 1- P ^ t + 1- P< 3>t 2 + ... + l- p ( n ) t 3.x 3 - . 3 T 3 . 3 1! 2! m! m ( C 2 ) where ,(m) _ a dt m t-0 We note that at t=0, w (speed deviation) = 0 and 6. = 6 ( i n i t i a l angle) therefore n P a 0 ) - V - £ <»1] " s 6 i J + «1J 6 l n 6i3> (C.3) 3 dt ^ 0 0, - Z (D s i n 6 - c . cos 6 ) t-0 J - l 3 J J J Since therefore P t=0 (1) _ u. . = 0 = 0 ( C 4 ) To s i m p l i f y the der i v a t i v e equations we l e t 127 n 0 0 A = Z ( D ^ cos 6 j L j + c ± j s i n 6^) (C.5) and n B = 2 ( D l j s i n 6 ±° - c ± j cos 6^) (C.6) the derivatives of A and B with respect to time are dt 1 J dt J (1) * (2) for convenience l e t 6^ y = 6^ , 6^ = 6 ... and so on therefore P<°> = P , - A a mi dt a ± j P<3) . d_ (1) = d_ (1) dt 3 dt « « A 6.!^ + B 6.(.2) = B c (2) "IJ I j " ° 6 ± ^ ' (C.7) P ( 3 ) . i . P ( 2 ) . 2 A 6 ( D 6 (2) _ B 6 ( D 3 + A 6 CD 6 (2) + B 6 (3) a . a i j i j i j i j I j IJ dt (C8) = 0 128 p CD P (1) where 6< 3> - - = 0 M i M j P<4> = 1_ P<3> = B 6< 4> + A 6 < 1 > 6 < 3 > a d t a i j i j i j + 2A 6±<2> - B 6 1J 1> 6±<2> + 2 A S ^ 1 ^ , ^ - 2B 6 ^ ) 6±<2> " 3 - \ f \ ? - A \ ^ P ( 4> - B 6 /*> + 3 A 6 <2) 2 ( C 9 ) a i j i j P<5> - *_ P<4> - B 6/5)+ 5A 6//>6//> + 1 0 A ^ "> B 6//> 26//> a d a i j i j i j i j i j i j i j - 15 B 61<j2>26i<j1> - 10 A 6 1< 1>\< 2> + B 6 ^ - 0 (C.10) P<«> = 1- P<5> = B 6 < 6 > + 6 A 6//>6//>+ 15 A 6< 2>6 <*>- 15 B 6 / ^ V / * a d a i j i j i j i j i j i j i j + 10 A 6±<3> - 60 B * ± ™ * £ h ™ - 20 A 6 ^ 6 ±< 3>- 15 B 64<2>" 45 A 61<l>6152>2 + 15 B 6 ™ " * ™ + A 6 ™ ' B 6 ^ + 15 A ^ 2 )6 ±5 4> - 15 B 6 ±< 2> 3 ( C l l ) 129 Note that 6 ^ ) = 6^3> - 6 ± \ 5 ) = 6±<7> ... = 0 130 APPENDIX D Derivation of p a ^ ^ an<* p a ^ ^ including damping: P < 3) - ±_ P <2> a i d t a i From equation (C.8) of Appendix C (3) P } ' (without damping) = 0 a l Therefore , ™ - - » » « > C D Also p {«> - i - P »> a l d t a i From equation (C.4) we obtain P <4> - B 6 <*> + 3 A 6 C2>2 - D 6 (5) a i i j i j i j Since 6 1 J 5 ) = 0 then P <4> - B 6 <*> + 3 A a/2)' a i i j i j - ^ ( D ^ - m 6 ±° - c i j C o s * ± ] ) * £ K 3 , A 0. <• (2)' + c l j S i n 6 ± j) 6 ± J (D.2) 131 APPENDIX E Taylor series c o e f f i c i e n t s of the accelerating power including e x c i t a t i o n response. The accelerating power of machine i i s P a i - Pmi " ^  V j <«lj C O S 5 i j + b i j S i n 6ij> ( E a ) and E = E(0) + F ( 0 ) t + l i - i t 2 + ... (E.2) 1 2 R e c a l l i n g equation (5.16) d E i E c ~ ( E c " V " t / T e Td 0 ; — - - F ^ - - — e - Er< v- xdi> i dt cos 4> S E j ( g ± J cos 6 ± J + b t j s i n 6 ± j) (E.3) the c o e f f i c i e n t s F ( 0 ) , F ( 1 ) are g i v e n i n S e c t i o n (5.3), F ( 2 \ F ( 3 ) are derived from (E.3) and given as follows. (E - E ) (1) n 0 F ( 2 ) = J±__c 0 J _ + _ s s i n 5 0 x 2 " S * i W ' J = 1 j j j e b c s 6 ° , CE.4) 132 ,(3) _ 1 ( E c " V ,(2) •z 3 cos • e "doi (E.5) The der i v a t i o n of the accelerating power c o e f f i c i e n t s are as follows Let A ± i = Ej ( g ^ cos 6^ + b ± j s i n 6 ±°) Bij • Ej ( 8U 8 i n 6ij " bij C°S \ ^ (B.6) then Let d A dt d B M - - B 6 ° ij ij 0 dt I 0 = 6 (1) g (1) = 6 (2) 6ij \i ' 6ij 6ij ' *• etc. Substituting equations (E.6), and (E.2) to ( E . l ) to obtain the Taylor series c o e f f i c i e n t s of P ^ we get P<°> - , , - I E<°> A , , a mi i i j (E.7) ,(D _ d Pai dt = - S ( F ^ A . , " E(0) B 6 t=0 j=l J 2 1 J i F ( 0 > A j=l ij (E.8) 4 t i \ d p < dt t=0 ^ I 6 f ( 1 > B i j 6 i f } + 3 E<°> A i j 6 i f > 2 133 >(2) _ d P a i dt n E (E(0) B 4 4 6 ^ 2 ) - F ( 1 ) A,,) t=0 j=l i j i j i j ' (B.9) ,(3) . ^ a i dt" t=0 n (E.10) + E(O) B ± j 6 ±< 4 ) - F ( 3 ) A ± j ] ( E . l l ) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0098151/manifest

Comment

Related Items