APPLICATION of CATASTROPHE THEORY to VOLTAGE STABILITY ANALYSIS of POWER SYSTEMS Thorhallur Hjartarson B.A.Sc. University of Iceland, 1988 A THESIS SUBMITTED IN PARTIAL FUMLLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 1990 © Thorhallur Hjartarson, 1990 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 of u^ l€c4ric&^ EL^ Qi*\eeP* The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract In this thesis catastrophe theory is applied to the voltage stability problem in power systems. A general model for predicting voltage stability from the system conditions is presented and then applied to both a simple 2-bus explanatory power system and to a larger more realistic power system. The model is based on the swallowtail catastrophe which with its three control variables is able to determine the voltage stability of the system. The model is derived direcdy from the systems equations. The voltage stability of the system at each specified system bus is determined by comparing the values of the swallowtail catastrophe control variables with those of the unique region of voltage stability. The control variables are calculated from the system operating conditions. If the control variables specify a point inside the stability region, the system is voltage stable; otherwise it is voltage unstable. ii Table of Contents Abstract ii List of Figures v List of Tables vi Acknowledgments vii 1 Introduction 1 1.1 Literature Review on Voltage Stability 3 1.2 Literature Review on Applications of Catastrophe Theory 7 2 A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. 9 2.1 Introduction 9 2.2 Developing the Model 9 2.3 Testing the Model 14 3 A General Catastrophe Model for Voltage Stability in Power Systems 19 3.1 General System Equations 19 3.2 Solving the System 20 3.3 General Swallowtail Catastrophe Modelling 26 4 Examination of the General Model on a Test System 29 4.1 Test System 29 4.2 Initial System Conditions 29 4.3 Voltage Variations 31 4.4 Load Variations 33 4.5 Load and Voltage Variations 37 5 Conclusions 39 Appendix A The Catastrophe Theory and the Swallowtail Catastrophe 41 A.1 Introduction 41 A.2 The Elementary Catastrophes 42 iii A3 The Swallowtail Catastrophe 44 Appendix B Program Code: General 6-Bus Case 48 References 69 iv List of Figures 1.1 Steady-State Stability Region 7 1.2 Transient Stability Region 8 2.3 Simple Power System 9 2.4 Bifurcation Set For The Swallowtail Catastrophe 14 2.5 Voltage Versus MVA At Load Bus (pf = 0.9285) 15 2.6 The Voltage Stability Manifold of the Swallowtail Catastrophe 16 2.7 Manifold foru < 0 16 2.8 Manifold for u > 0 17 2.9 The Travel of the Voltage Curve through the Catastrophe 18 2.10 Voltage Curve near Critical Point 18 3.11 Non-Radial Power System 22 3.12 Bifurcation Set For The Swallowtail Catastrophe 28 4.13 A 6-Bus Power System 30 A.14 The Tail of the Swallow 46 A.15 Manifold for s = 0 47 A.16 Manifold for s > 0 47 v List of Tables 2.1 The Travel of the Voltage Curve through the Catastrophe 17 4.2 Load Flow Results for 6-Bus Power System 30 43 Control Variables for Initial Conditions 31 4.4 Control Variables with 82.6% Initial Voltage Values 32 4.5 Control Variables with 58% of Voltage Values at Bus 1 33 4.6 Control Variables with 52% of Voltage Values at Bus 2 33 4.7 Control Variables with 29% of Voltage Values at Bus 3 34 4.8 Control Variables with 72% of Voltage Values at Bus 4 34 4.9 Control Variables with 70% of Voltage Values at Bus 5 35 4.10 Control Variables with 5.5 times the Initial Reactive Loads 35 4.11 Control Variables with 240 times the Initial Reactive Load at Bus 1 36 4.12 Control Variables with 465 times the Initial Reactive Load (small) at Bus 2 36 4.13 Control Variables with 59 times the Initial Reactive Load at Bus 3 36 4.14 Control Variables with 11.3 times the Initial Reactive Load at Bus 4 37 4.15 Control Variables with 8.5 times the Initial Reactive Load at Bus 5. . . 37 4.16 Effects of both Load and Voltage on Voltage Stability 38 vi Acknowledgments I would like to thank my supervisor Dr. M.D. Wvong for his invaluable insight and assistance with the work involved in this thesis. I would also like to thank all friends and family who have been supportive through the ordeal of my educations era. Special thanks to friends Ivar and Iris and of course their new bom daughter Alma, for constantly feeding me real home made meals and providing me with good company during lonely winter nights. Also special thanks to all the friends who provided many motivating discussions and interchanging of views over the odd beer. vii Chapter 1: Introduction Chapter 1 Introduction The blackout problem of electric power systems has traditionally been associated with the steady state and transient stability problems. Steady state and transient stability are the phenomena involved in connection with the loss of a major portion of a grid due to the inability of certain generators to maintain synchronism in the face of small and large disturbances respectively. These types of instability are, generally speaking, well understood today. System stability is being preserved to a greater extent than ever before by the advent of faster and more effective stabilizers, and more reliable protection systems. In recent years a category of instability, usually termed voltage instability or collapse, and associated with the inability of a power system to maintain bus voltage magnitudes, has been responsible for several major blackouts world-wide. Earlier stability problems were concerned with the relationship between active power and phase angles of generators. The static voltage stability problem is, on the other hand, concerned with the relationship between reactive power and voltage magnitudes of generators. As power systems become heavily loaded, there is a possibility that the power systems might suffer from a cascading voltage collapse due to lack of reactive power. Also there is the dynamic voltage stability problem, which involves frequency as a parameter as well. In this thesis, only the static voltage stability problem is considered since the dynamic case is more complex and the static problem is currently being researched in greater detail. As frequency is a critical parameter in the balance between real (MW) generation and real (MW) load throughout the power system, so transmission voltage levels reflect the balance between the supply and demand of reactive power. While frequency is uniform throughout the power system, voltage levels can vary markedly across a transmission network, which is designed to operate at a particular voltage level. As a result, it is generally accepted that the voltage stability problem, 1 Chapter J: Introduction which is associated with the inability of a power system to maintain bus voltage magnitudes, is due to a deficit of reactive power at certain buses in the network. The actual process of collapse may therefore be triggered by some form of disturbance, resulting in significant changes in the reactive power balance in the system. The operating environment of many present-day power systems substantially increases the vulnerability of the system to reactive deficit problems and therefore difficulties in maintaining system voltage profiles. Several factors have contributed to this situation. There is increasing difficulty in obtaining power plant sites in the vicinity of major power consumers. Also, the exploitation of hydro power resources has proceeded spectacularly to a point where remote, large generation plants have been developed. As a result, electrical power is often transported through high capacity lines over long distances from generators to consumer. Furthermore, the strengthening of transmission networks has been curtailed in general by high costs, and in particular cases by the difficulty of acquiring right-of-way. This has resulted in increased loading and exploitation of the older circuits thus resulting in increasing voltage stability problems. Other factors include the relative decrease in the reactive power outputs of generating units, and shifts in power flow patterns associated with changing fuel costs and generator availability. Numerous approaches to predict voltage collapse have been suggested (see following section). A fast alternative method is to apply catastrophe theory to the reactive power equations where an unique solution set for the stability exists. If the solution according to the catastrophe theory for any particular operating condition falls within this unique set then the system is voltage stable for that particular condition. In this thesis, the swallowtail catastrophe is applied to the voltage stability problem. In Chapter 2 a well known simple voltage stability example is modeled and examined. A more general model for any interconnected power system is derived in Chapter 3 which is then, in Chapter 4, applied 2 Chapter 1: Introduction to a larger more realistic example. 1.1 Literature Review on Voltage Stability One of the first to address the problem of voltage collapse was W.R.Lachs [1]. There the phenomena is explained and an example is given which shows how a voltage collapse can occur. The importance of reactive compensation is discussed in detail since the reactive compensation on a EHV system must provide for both an overall, and a regional, balance and be able to withstand any feasible reactive disturbance on the system. In [2] a more recent review of the problem was done. Work about voltage stability conditions, proximity indicators, control strategies and planning network reinforcement is discussed. Some of the points made follow here. Some work has been reported on defining and establishing voltage stability criteria, i.e., criteria that may be used to determine whether or not an operating condition is stable from a voltage stability viewpoint It is suggested that an operating condition is stable from the voltage viewpoint if every load bus voltage increases when a source voltage increases or when a shunt capacitor is switched in at a load bus. Transformer taps are a major contributing factor in system voltage collapse and voltage instability is characterized in association with the slow tap-changing transformer dynamics. Stability conditions are derived in terms of allowable transformer taps settings using eigenvalue analysis. One criterion for voltage stability of a given operating condition states that for an operating point, voltage stability is ascertained when at mat operating point, an elementary increase of reactive demand is met by a finite increase in reactive power generation. Analytical computation of an index, which is defined on the basis of this criterion, for the simple two-bus system has been presented. 3 Chapter 1: Introduction As noted in [3], a region-wise framework is presented which accounts for several of the mechanisms to predict voltage collapse and can serve as a basis for comparing the effectiveness of performance indices, to predict, on-line, voltage collapse problems in power systems. The basis of the framework is the voltage stability region, which accounts for both static and dynamic mechanisms of voltage collapse. Based on this region, static and dynamic performance indices are denned to predict the static mechanisms of voltage collapse in the input or injection space and the dynamic mechanism in the post-contingency state space, respectively. Sekine et al. also studied the dynamic phenomena of voltage collapse [4] by the method of dynamic simulation using induction motor models. From the viewpoint of dynamic phenomena, the voltage collapse starts locally at the weakest node and spreads out to the other weak nodes. In [5] the stability limit problem as a static divergence or bifurcation characterized by the disappearance of an equilibrium point was studied. Beyond this point, solutions to the load flow equations cannot be obtained. In [6] a suggestion on the use of the minimum singular value of the Jacobian matrix of the load flow equations as a security index was made, and static control strategies based on this index were derived. Mori et al. presented a method [7] for estimating critical points on static voltage instability in electric power systems. The difficulty of evaluating critical points in y (specified value) space although it is easy to find a critical point in x (voltage) space with the conventional method is clarified. This is because there exists a set of specified values that provides a singular point in voltage space. Emphasis is put on evaluating the singular points in y space that is closest to operating conditions. A nonlinear programming technique is utilized to evaluate critical points in voltage and specified value space. In [8] a voltage stability index based on the feasibility of solution to the load flow equations at each bus was developed. 4 Chapter I: Introduction In [9] a security measure to indicate vulnerability to voltage collapse based on an energy function for system models that include voltage variation and reactive loads is denned. The system dynamic model, the energy function and the security measure are first motivated in a simple radial system. Application of the new measure and its computational aspects are then examined in a standard 30-bus example (New England System). The new measure captures nonlinear effects such as var limits on generators that can influence the vulnerability to collapse. The behavior of the index with respect to network load increases is nearly linear over a wide range of load variation, facilitating prediction of the onset of collapse. One common drawback of all these methods is that the operating constraints on system equipment (i.e., MW and MVAR limits of system generating units) are not taken into consideration. Production capabilities of generating units are important considerations, more so since voltage collapse is considered to be a reactive power problem. One suggestion is the use of repeated load flow computations, as power injections are increased, to determine the voltage stability limit. Having determined the limit, the margin to collapse is then available. However, besides being computationally very demanding, this approach may be inadequate due to the unreliable behavior of the Newton-Raphson method of load flow analysis in the vicinity of the voltage stability limit. This behavior is linked to the singularity of the Jacobian at the voltage stability limit, and the existence of close multiple load flow solutions around that limit. Another suggestion is the use of a combination of load flow analysis and sensitivity parameters in order to reduce the computation time and circumvent the numerical ill-conditioning known to occur as the voltage stability limit is approached. Another possible short-cut to the repeated load flow calculations is a quadratic extrapolation method which Borresmans et. al. proposed. These methods are inherently approximate. The third suggestion is an approximate method to evaluate the condition at the voltage stability 5 Chapter I: Introduction limit In [10] the problem of voltage collapse is approached by simulation of a power system using a slightly modified transient stability program. A sufficiently complex system with 39 buses and 10 generators is used in simulations. The system is stressed by progressively increasing the system load through a multiplier k. A very small change of k (order 1%) is used as collapse-inducing disturbance. Total system voltage collapse was observed after the disturbance. Significance of these findings and directions for future research are presented and applicability of real-time control discussed. Several voltage collapses have had a period of slowly decreasing voltage followed by an accelerating collapse in voltage. In [11] this type of voltage collapse based on a centre manifold voltage collapse model is analyzed. The essence of this model is that the system dynamics after bifurcation are captured by the centre manifold trajectory and it is a computable model that allows prediction of voltage collapse. Both physical explanations and computational considerations of this model are presented. The use of static and dynamic models to explain voltage collapse are clarified. Voltage collapse dynamics are demonstrated on a simple power system model. In [12] a method of determining the voltage stability limit of a general multimachine power system was presented. In the method, the search for the voltage stabihty limit is formulated as an optimization problem of maximizing the system total MVA load. With this formulation, difficulties related to singularity of the load flow equations Jacobian matrix, and convergence of the load flow solution around the voltage stability limit, are avoided. Voltage dependence of the MVA loads may be taken into account in determining the voltage stability limit. The method also accommodates device constraints or limitations in system controls (e.g. generator VAR limits and limits on transformer tap settings). A security margin is denned which serves as a measure of the security of the system as far as voltage collapse is concerned. 6 Chapter 1: Introduction 12 Literature Review on Applications of Catastrophe Theory In [13] M.D. Wvong and A.M. Mining introduce the application of catastrophe theory to transient stability assessment of power systems. A derivation is shown for a one-machine infinite-bus power system for which a swallowtail catastrophe manifold is derived from the swing equation of the system. From this a stability region is defined which is valid for changing load conditions and fault locations so that the stability of the system can be easily assessed. Figure 1.1 shows the steady-state stability region while Figure 1.2 shows the transient stability regioa This method of transient stability assessment has several advantages: • The regions of stability obtained by the method are well defined in terms of the control variables and the critical clearing time regardless of the state variables. • The computation required to define the stability are few and can be done in a very short time. • The method may be used for an on-line assessment of transient and steady-state (dynamic) stability. 7 Chapter I: Introduction Figure 1.2: Transient Stability Region In [14] A.A. Sallam did similar studies on steady state stability assessment in power systems with similar results. Catastrophe theory has also been applied to other fields in science with success. Some of the work done is described below: In [15] catastrophe theory is used to classify caustics and their associated traveltime diagrams. These traveltime diagrams are multivalued functions of the position coordinate when caustics are present. In [16] an experimental and theoretical investigation is presented for the forced vibration of a one-degree-of-freedom system with a non-linear restoring force. It is shown that the characteristics of these systems can be described by the cusp catastrophe model. And in [17] and [18] the rainbow effect in ion channeling is analyzed by the use of the catastrophe theory. 8 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. Chapter 2 A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. 2.1 Introduction In this chapter we are going to look at the application of catastrophe theory to voltage stability problems in a simple power system. This will give us insight into how the catastrophe model is developed. In later chapters we will do the same for a more realistic interconnected multi-bus power system. The simple 2-bus power system is shown in Figure 2.3. The problem is formulated with regard to voltage collapse at the load bus (bus 2). V1P1.Q1 V2P2.Q2 Source Load Figure 23: Simple Power System 22 Developing the Model A catastrophe model for voltage stability is derived from the following system equations: Miori + (ViV2/Xi)sin(*1 - S2) = Pi (2.1a) -(ViV2/Xi)cos(tf2 - Si) + (l/Xi)V? = Q i (2.1b) (V1V2/X1)sin(^2 - Si) = P 2 (2.1c) 9 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. (ViVj/XOcosto - *i) - ((1/Xi) - B)V?. = Q 2 (2.1d) Since we are only dealing with static voltage stability our angles do not change, hence 6j = 0. In our case bus 2 (the load bus) will be looked at from a voltage stability viewpoint. Therefore equation (2.Id) will be the equation that the catastrophe model will be built on. It will be the changes in V2 diat will be approximated. But since Vi is also present in equation (2.Id), we need to solve for it in terms of V2. This can be done from equation (2.1b) as follows: V 2 - V , V 2 cos (<52 -6^- Q J X J = 0 (2.2) This is a quadratic equation with the solution: Vi = \ (v2 cos(62 -60± v / (V 2 2 cos2 (^ 2 -^) )+4QiX 1 ) (2.3) Substitution of equation (2.3) into equation (2.1d) will give [(Vl/(2X0)cos2(*2-^) +(V 2 / (2X 1 ) )cos(« 2 - ^)^/(V 2 2 cos2(^-^ 1 )-r4Q 1 X 1 )] -((1/X,) - B)V22 = Q 2 This can be rearranged as: + V , ( C ° S '>>) ^ V | co.' (h - «,) + 4Q.X, where Q 2 is a potential function of the variable V2. The catastrophe manifold Ms is the set M s = {V a eR: / s (V 3 )=0} (2.4) (2.5) (2.6) 10 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. where fs is the derivative of our potential function with respect to the state variable V2: / S(V 2) = 2V,[(C 0 S'<*-'•>) - ( J - _ B)] y ^ « r f « i - W + 4Q1X1 + 2(C°S(2^ \Zi\, where ZR is the load impedance and ZL is the line impedance. Figure 2.5 shows the variation of VR, the load bus voltage, against MVA demand SR, at constant power factor. Point A (where VR = V^LT) represents the critical system state. The upper segment (VR > V^,{) is considered the stable operating region, i.e. < 1. It can be seen that in the stable region, increasing the sending end voltage increases the receiving end voltage whereas, in the unstable region, increasing the sending end voltage actually reduces the receiving end voltage. The catastrophe model was compared with the above analysis in the following manner. Two voltage points were taken from the graph in Figure 2.5, where one point was in the voltage stable area and the other in the unstable area. By feeding the same data used in Figure 2.5 into the catastrophe 14 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. Figure 2.5: Voltage Versus MVA At Load Bus (pf = 0.9285) model it could be determined whether the model agreed with the analysis. All the stable points were found to lie above the catastrophe manifold shown highlighted in Figure 2.6, while all the unstable points were found to lie beneath it The critical point A in the voltage curve shown in Figure 2.5 was found to be close to S = 1.576397 p.u. The turning point in the curve according to the Catastrophe model was almost the same or close to 5 = 1.576292 p.u. This is the point at which our system shows stability as suggested by Figure 2.6. For the swallowtail catastrophe shown above there are basically three different cases, with only one being stable, that could result from our analysis. These are, above the manifold as shown in Figures 2.7 and 2.8, inside the tail as shown in Figure 2.7, and outside the manifold as shown in Figures 2.7 and 2.8. 15 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. V Figure 2.6: The Voltage Stability Manifold of the Swallowtail Catastrophe W Outside Manifold Figure 2.7: Manifold for u < 0 It was found that almost all the points of the voltage curve, stable or unstable, resulted in a control variable u > 0. Only the points close to the critical point of the curve resulted in a control variable u < 0. As the voltage curve went from unstable to stable (the critical point) it travelled from 'outside the manifold' to 'inside the tail', the two unstable regions, to 'above the manifold', the stable 16 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. W i ^ ^ b o v e 1. Outside Manifold Figure 2.8: Manifold for u > 0 Regions of the Catastrophe Intervals of S (p.u.) Intervals of V? (p.u.) Outside Manifold, u > 0 (Unstable) 0-1.569 0-0.598 Outside Manifold, u < 0 (Unstable) 1.569-1.573 0.598-0.609 Inside Tail (Unstable) 1.573-1.576 0.609-0.630 Above Manifold, u < 0 (Stable) 1.576-1.573 0.630-0.658 Above Manifold, u> 0 (Stable) 1.573-0 0.658-1.05 Table 2.1: The Travel of the Voltage Curve through the Catastrophe region. Table 2.1 shows the points at which the voltage curve entered each region of the catastrophe. Similarly this can be shown as Figure 2.9 illustrates where the points shown are referred to Figure 2.10, which is a close up of the voltage curve of Figure 2.5 near the critical point. It should be noted that since the catastrophe is actually three-dimensional, the curve in Figure 2.9 is not exact It only serves for explanatory purposes and is not drawn to scale. 17 Chapter 2: A Catastrophe Model for Voltage Stability in a Simple 2-Bus System. W A Figure 2.9: The Travel of the Voltage Curve through the Catastrophe Figure 2.10: Voltage Curve near Critical Point 18 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems Chapter 3 A General Catastrophe Model for Voltage Stability in Power Systems In the previous chapter, voltage stability was studied for a simple 2-bus power system using a swallowtail catastrophe model. To apply the method to larger, more realistic, power systems it is useful to develop a general model for voltage stability of any power system. 3.1 General System Equations A general swallowtail catastrophe model can be developed for an n-bus power system directly from the systems equations. The most generalized way of stating the power flow equations of a power system [19] is n Si = Pi+jQi = ^{YsViVj + M^+jD^} j=i D = E { l V ' l l V iK G « " JBijJetoe-* + MjS] + jDj*j} j= i = E { l v i l l v i l(G S - JBiiXcos (Si - 6j) + j sin (Si -Ss)] + MjS] + j D ^ D = E { l v i l l v i lK G » c o s (ft ~ ft) + By sin (4 - Sj)) -rj(Gij sin (Si - Sj) - By cos (S{ - Sj))] + MjS] +}DjSi} (3.22) Since we are only dealing with static voltage stability our angles do not change, hence Sj Now by separating the real and imaginary terms: Pi = E l v i l l v j l[G i i c o s (ft - ft) + B « s i n (4 - ft D (3.23) Qi = E l v i l l v j l [G i i s i n(ft - ft) - B « c o s (ft - ft 19 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems or since we have P i = E irSi^cos (ft - ft) - x * s i n (ft - ft)] j=i Kij + A « ^ = E i rSi t^ s i »(ft - ft) + xy cos (ft - 6^)] (3.25) Again since static voltage stability is being modelled, all phase angles remain constant. It is only in dynamic voltage stability that changes in phase angles occur. The following constants can therefore be denned: Rij cos (Si — 6}) — Xij sin (Si - ft)\ <8 = u y 7 (3.26) , _ / Rij sin (ft - ft) + Xij cos (Sj - 6$) * = I Ri + Xi , Voiiagc stability is synonymous with reactive power stability where voltage is the dynamic variable (likewise, the phase angle is the variable in real power stability studies). For this reason only the reactive power equations are of interest in this study. Therefore our set of equations is: Q i = £ I V M I V ^ (3.27) j=i 3.2 Solving the System We now define any bus k as the one to be examined from a voltage stability viewpoint Its equation on which the catastrophe model will be built is a single case of the set of equations above or D Q k = ElVkHVjIVkj (3.28) 20 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems where Qk is a potential function of the variables VJt,Vj (j = l,2,..,n). But these variables are related through the other system equations. Therefore by solving the system equations the variables Vj (j = 1,2, ..,n) could be stated in terms of only V*. Essentially what we have is a set of n functions of the form: Qi = fi(Vi,Vj) (3.29) which can be rearranged as V i = g i(V J) (3.30) if Qi is considered to be constant. However, Q^ is not considered a constant since it is our potential function. Al l the reactive powers effectively vary, but for buses other than the bus under study the variation is small enough to be neglected. Thus by going through the system and solving at each node i for the corresponding voltage Vi and substituting into the next equation and so on, our potential function can be written in terms of V% only. For a strictly radial system this can be done in a straightforward manner and does not need further details. On the other hand for a non-radial system (loop or circular connections) problems will arise because of the nonlinearities of the equations. Figure 3.11 shows an example of this. Here we have four buses which are all interconnected, and have two connections to the rest of the system at buses u and z. To be able to solve the whole system we therefore need to be able to solve implicitly for the voltages at buses u and z. 2 1 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems bus U other buses Figure 3.11: Non-Radial Power System The functions at the four buses upc,y,z will be Qx = fx(vx,vu,vy,vz) Qy = f y ( V y , V u , V x , V z ) (3.31) Qu = fu(V u? V x , V y , V z , Vj) Qz = fz(V z , V u , V x» Vy,Vj ) where K and V, are voltages at other buses connected to buses u and z. This can be rewritten as V ^ g ^ V v . V , ) Vy = g y ( V u , V „ V x ) (3.32) V\l = gu(V x , Vy , V z , Vj) • V z = g x ( V u , V x , V y , V j ) This cannot be solved explicidy for Vu or Vz and therefore voltages V{ and Vj can't be solved either and so on. 22 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems Note, however, if VX has an explicit solution in terms of V u , i.e., V x = gx(Vu) (3.33) then Vu can be solved as V u=g u(V i,V x) = g u(V i,g x(V u)) = hu(V0 (3.34) where bus i precedes bus u in sequence. The same then goes for Vj etc. Therefore we need to write VU and VZ in the same form as shown in equation (3.33), that is with an explicit solution of only one variable. Since we are dealing with small changes in voltage and since these changes are going to be smaller the further away we are from the specific bus of study they can be dismissed as shown below. When solving for Vu, then VX and VY can be written as where VXQ, Vvo and V2o are the initial voltage values at buses x,y and z. Now we rewrite the equation Vx«gxi(V u,V ¥ O,V l 0) = h x(V u) (3.35) Vy «gyl(V u,V x 0,V z 0) = MVu) for Vu as V u = g„(V x ,V y ,V„Vi)wg u (g x i (V u ,V y o ,V x o) ,g > - i (V u ,V x 0 ,V ,o) ,V l 0 ,Vi ) (3.36) or V u = hu(Vi,Vxo,Vyo,V zo)=nu(Vi) (3.37) since Vro, Vyo and Vzo are constants. Similarly when solving for VZ, VX or V^ can be written as V x»6« 2 ( V „ V y o , V u o ) = hx(V,) (3.38) V y «gy 2(V z ,V xo,V uo) = h y(V z). 23 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems Again we rewrite the equation for VZ as V, = gz(V u,V x,V y,Vj) * gz(V uo,gx2(V z,V uo,V yo),g y 2(V u,V uo,V x 0),V j) (3.39) or V, = h,(Vj, V u 0, Vxo.Vyo) = h,(Vj). (3.40) Each of the solutions of equation (3.30) will be a quadratic solution of equation (3.29) which can be rewritten as Qi = Vi J2 VjV^ + v? £ V2ij j=lj#i j=lj# where Vj = VJQ when systems are inteconnected as described above. E VjVio Vf + V i J - ^ Qi E V>2ij E V^ij = 0 which will give the solution: , E Vi^ ia 2 D E 2« + 4Q; D \ 2- va / E V^ij N \ i=ij# / i=to#i This will, however, become unsolvable as we go through the system since the degree of the solution will increase by one at each step. But this problem can be overcome by approximating all voltages at buses far away1 from the specific bus as (3.41) (3.42) (3.43) Vi = V i 0 + AVi (3.44) 1 far away = all buses but the ones next to the specific bus. 24 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems by disregarding higher orders of AVj. Now equation (3.41) can be solved simply for AV, or even better it can be solved in the same way for Vi = V;o + A V ; . Qi^CVio + A V O J2 VjVis + (Vio + A V ; ) 2 E V>2ij = (Vi0 + AVi) £ VjVia + (Vi20 + 2Vi0AVi) E V*i (Vio + AVO E VjViij +(2V i 0(V i 0 + AVi)-V20) E ^ (3.45) = (V i 0 + A V i ) E V J ^ S + 2 V i o E ^ -v 2 0 E ^ Therefore, Qi + v2o E V>2ij Vi = V i 0 + A V i = —s J ~ 1 J ? " „ (3.46) E VjViij + 2V i 0 E ^ However, the solutions of Vj which have previously been solved in terms of V,- could have been substituted into equation (3.46) leaving only one unknown voltage variable. We shall call this variable, Vj. Now equation (3.46) can be written as Vi = r -^-v (3.47) li + mjVi where fc;,/,,m, are constants. This solution of Vj is then substituted into the remaining system equations to solve for V/ in the same way. A l l buses, except those near bus k, are solved this way. Since the voltage changes are larger the closer to the specific bus (the voltage collapse bus) more accurate voltage solutions are desired. This time the precise solutions of equation (3.43) are used for solving these buses next to the specific bus. The format of that equation will change somewhat though since all known Vj solutions of equation 25 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems (3.47) have been substituted in. The solution will be of the following form: V i = -^(ai + AVO + ^ ' + A V O ' + l i = ~ ( a i + A V k ) + \ JtfVl + 2aiAV k + a? + 7 i (3.48) where £ is the specific bus. The potential equation at bus it is similar to equation (3.41) or, D II Qk = Vk £ VjVikj + Vl Y, ^2kj (3.49) where all the Vj's have previously been solved so that only one variable is left in the equation, namely Vk. 33 General Swallowtail Catastrophe Modelling Since our potential function Qk is of the same form as previously shown for a 2-bus example, the steps involved in developing the catastrophe model are basically the same from here on and there is no reason for going through them again in detail. As said previously, the catastrophe manifold Ms is the set Ms = {Vk 6 R : f(Vk) = 0} (3.50) where fs is the derivative of our potential function Qf with respect to the state variable Vk By applying Taylor's approximations of functions for the changes in voltage, the voltage at any time can be written as: V 2 = V 2 0 + AV 2 + A V | A V | A V J 2 + 6 + 24 (3.51) or with x = A V 2 as _ _ . . X X V 2 = V 2 „ 4 - x + _ + - + - (3.52) 26 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems After approximating the square roots with a 4th order Taylor series and collecting like terms we get: a4X4 + aax3 + a2x2 + aix + ao = 0 (3.53) where ao = V k 0 52 TJVikj + V k 2 0 £ V2kj D D a i = E TjV'ikj+2Vko 52 ^2kj j=lj#k j=lj*k ^ 2 = 2 E T/Vikj + (V k 0 + 1) 52 V>2kj (3.54) j=lj#k j=lj/k *> = \ E T J W ^ + I) E ^ J=U#k v ' j=lj#k £r!(! Tj j? the r-th constant in the corresponding taylor series. Now by letting, x = y + a , where, a = our equation can be written as: y4 + uy2 + vy + w = 0 (3.55) where u,v and w are our control variables: If 3 a?A u = — I a2 ) a4 \ 8 a4 / If a2 , 1 ai! \ _ v = — a i - — — + - ^ | (3.56) &4\ 2 an 8 a| / _ ao aj an a2 a2 3 a4 W ~ a 7 ~ T a | l 6 a i ~ 2 5 6 a i 27 Chapter 3: A General Catastrophe Model for Voltage Stability in Power Systems This is precisely the same swallowtail catastrophe model for voltage stability as derived earlier for our simple 2-bus example. The bifurcation set for this catastrophe is shown here again in Figure 3.12. 28 Chapter 4: Examination of the General Model on a Test System Chapter 4 Examination or the General Model on a Test System 4.1 Test System Even though the catastrophe model had proven to work for a simple 2-bus system, studies needed to be done on larger systems for realistic purposes. This realistic system was needed to have the following four characteristics and conditions for our purposes: • Non-radial system. • Losses accounted for. • Load flow results available. • Generator and load buses. The system chosen [20] was a 6-bus single loop system as shown in Figure 4.13. This system has line-losses parameters available, a load flow has been computed and it has two generator buses and four load buses. Table 4.2 shows the load flow results which are needed for the catastrophe analysis. Bus 6 of the system is a slack bus which means that the voltage remains constant and the demand of power required would always be met. Thus bus 6 was excluded from our cases of voltage stability studies (No voltage change, no voltage collapse). All the other five buses were of course examined though as the following sections show. This system was then formulated according to the general model described in Chapter 3. 4.2 Initial System Conditions First the system was analyzed using the initial load flow conditions given in Table 4.2. Table 4.3 shows the results from the catastrophe program at each of the five buses. The stability and the unstability regions are identical to the ones shown and discussed in Chapter 2 since we are still dealing with the swallowtail catastrophe. 29 Chapter 4: Examination of the Genera! Model on a Test System 0.08+J0.37 p.u. 0.123+j0.518p.u. t=0.974 •J34.1 p.u. _ l jO.133 p.u. t=0.909 0.723+j 1.050 p.u. jO.3 p.u. 0.282+J0.64 p.u. -J28.5 p.u. Figure 4.13: A 6-Bus Power System Bus Number Voltage Magnitude, V (p.u.) Voltage Angle, 6 (°) Reactive Power, Q (p.u.) 1. 0.932061 -9.841571 0.218091 2. 1.003495 -12.770030 0.016298 3. 1.103648 -3.389008 0.370004 4. 0.924158 -12.282888 0.179988 5. 0.922237 -12.182555 0.153299 6. 1.05 0 0.497890 Table 4.2: Load Flow Results for 6-Bus Power System As expected the system is voltage stable at all buses for the initial conditions but of interest 30 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9888 7.9888 11.7420 Stable 2. 1.1586 2.5105 3.9941 Stable 3. 5.3547 7.1235 17.1505 Stable 4. 1.0936 3.3229 3.1817 Stable 5. 5.1193 6.9659 7.1436 Stable Table 4.3: Control Variables for Initial Conditions is to find out which bus is the least or most stable and so on. In later sections we will vary some of the state variables and then see clearly which buses are more sensitive to voltage collapse. But it is possible to predict from the above control variables which buses would be most susceptible to voltage instability although such a prediction might not be very accurate. This kind of prediction would be based on the proximity to the catastrophe manifold. Each of low control variable u, high v and low w by themselves or especially in combination could indicate this proximity. By looking at Table 4.3 it can be seen that both buses 2 and 4 have a relatively low control variable u. Similarly buses 1 and 3 have a relatively high control variable v and buses 2 and 4 have a low control variable w. From this it could be concluded that again buses 2 and 4 would be more vulnerable to voltage collapse than buses 1, 3 and 5. This would make sense since both these buses (2 and 4) are the load buses furthest away from the two generators. 43 Voltage Variations As previously stated the two state variables for static voltage collapse are voltage and reactive power. Lower voltage and/or higher reactive loads would make a power system more vulnerable to voltage instability. In this section the effects of strictly voltage variations on voltage collapse are examined using the catastrophe model. Voltage values were lowered until the model predicted a voltage collapse 31 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9887 7.9887 7.6945 Stable 2. 0.7667 1.8239 3.9941 Stable 3. 5.3014 6.9908 12.3260 Stable 4. -661.50 6940.07 -19200.91 Unstable 5. 5.2232 7.1251 9.6979 Stable Table 4.4: Control Variables with 82.6% Initial Voltage Values. at one of the buses. However since a new load flow would effectively be needed for each new voltage profile the exact voltage value at which the collapse occurs may not be accurate. But the purpose of this test is not to find the exact voltage collapse points for the power system in Figure 4.13 as much as it is to prove that the catstrophe model works for voltage collapse predictions. The first test was done by varying all the voltage profiles equally. Voltage collapse was found to occur at bus 4 at 82.6% of the value of the initial voltages as shown in Table 4.4. The other buses are still stable at this point but the whole system would collapse though due to the instability at bus 4. Individual voltage variations are also of interest although there are of course infinite combinations of the different bus voltages. The tests done here involve varying only one bus voltage value while the other voltages remain constant at the initial values. Since bus 4 proved to be the weakest one in the uniform voltage variation test above it was of particular interest for this examination and as can be seen from the following tables it still has the highest bus voltage for which it goes unstable (the weakest bus). Studies done on bus 3, a generator bus, found that voltage collapse occurs at a very low voltage 32 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9888 7.9888 2.3813 Unstable 2. 1.1586 2.5105 3.9941 Stable 3. 5.4484 ' 7.2298 19.5364 Stable 4. 0.9080 2.8032 2.7308 Stable 5. 5.1121 6.9491 6.8755 Stable Table 4.5: Control Variables with 58% of Voltage Values at Bus 1. Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9886 7.9886 11.7386 Stable 2. 0.2643 2.7717 1.4320 Unstable 3. 5.2575 7.0099 15.5127 Stable 4. 0.8394 2.6126 2.5673 Stable 5. 5.1188 6.9471 6.7781 Stable Table 4.6: Control Variables with 52% of Voltage Values at Bus 2. (29%) and the model cannot be used for some of the other buses at this voltage (buses 1 and 5). This is not at all unreasonable since there is no way the system would support these voltage conditions. Also voltage collapse rarely occurs at generator buses because of greater voltage control. 4.4 Load Variations Now we will examine the effects of reactive load variations with the voltage constant at the initial values. This testing was done in a similar way as before, only in this case by increasing the loads 33 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. - - - -2. 1.1381 2.4769 3.9709 Stable 3. 5.0389 6.8004 2.0884 Unstable 4. 1.0709 3.2673 3.1439 Stable 5. - - - -Table 4.7: Control Variables with 29% of Voltage Values at Bus 3. Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9887 7.9887 10.6383 Stable 2. 1.1536 2.5235 4.0283 Stable 3. 5.3788 7.1301 17.0745 Stable 4. 28.7936 131.5285 131.0456 Unstable 5. 5.1146 6.9520 6.9110 Stable Table 4.8: Control Variables with 72% of Voltage Values at Bus 4 . until voltage collapse occurred at one of the buses. The first test involved increasing all the bus loads uniformly until voltage instability was reached. This happened at 5.5 times the initial loads as shown in Table 4.10. Again it was bus 4 which failed first. The tests for the individual changes of loads at each bus with the other loads constant at the initial values provided some interesting results as the following tables show. The model does not 34 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9888 7.9888 12.0828 Stable 2. 1.1018 2.3576 3.8180 Stable *> J). 5.4118 7.1429 17.0061 Stable 4. 1.0936 3.3229 3.1817 Stable 5. 5.0856 6.9298 2.1775 Unstable Table 4.9: Control Variables with 70% of Voltage Values at Bus 5 . Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9893 7.9893 12.2171 Stable 2. 3.9361 8.6089 10.4715 Stable 3. 5.6998 7.4385 17.3363 Stable 4. 15.3072 12.5569 2.3516 Unstable 5. 5.1128 6.9538 6.9585 Stable Table 4.10: Control Variables with 53 times the Initial Reactive Loads. provide a stability census at the load change bus simply because the reactive power at that bus is not a variable in the catastrophe model. This is because the potential function of the catastrophe model is the reactive power equation at the bus in question. This is fully acceptable since for the value of the function to change then the variables of the function must change as well. As before only the initial load flow data was used, so the data is not accurate, but this is sufficient to show that the catastrophe model is valid. This inaccurate data is probably the reason for the high 35 Chapter 4: Examination of the General Model on a Test System Bus Control Control Variable, Control Stability Census Number Variable, u V Variable, w 2. 5.6703 12.0631 5.4013 Unstable 3. 5.6220 7.4188 31.6700 Stable 4. 0.7265 2.4308 2.5738 Stable 5. 5.1047 6.9267 6.5458 Stable Table 4.11: Control Variables with 240 times the Initial Reactive Load at Bus 1. Bus Control Control Variable, Control Stability Census Number Variable, u V Variable, w 1. 5.7234 8.2060 9.6336 Stable 3. -5.9235 -14.7004 -16.2437 Unstable 4. 0.9568 2.8519 2.6492 Stable 5. 5.1190 6.9724 7.2764 Stable Table 4.12: Control Variables with 465 times the Initial Reactive Load (small) at Bus 2. Bus Control Control Variable, Control Stability Census Number Variable, u V Variable, w 1. 5.9901 7.9901 14.5094 Stable 2. 1.7603 3.1642 1.5576 Stable 4. 1.2349 3.2494 1.4174 Unstable 5. 5.1954 7.0921 10.6266 Stable Table 4.13: Control Variables with 59 times the Initial Reactive Load at Bus 3. 36 Chapter 4: Examination of the General Model on a Test System Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9887 7.9886 11.5493 Stable 2. 0.9214 1.9520 3.3876 Stable 3. - - - -5. 2.6522 .10.5990 -0.4502 Unstable Table 4.14: Control Variables with 11.3 times the Initial Reactive Load at Bus 4. Bus Number Control Variable, u Control Variable, V Control Variable, w Stability Census 1. 5.9888 7.9888 11.6117 Stable 2. 1.1765 2.5387 4.0117 Stable 3. 5.2966 7.1363 17.7766 Stable 4. 16.0800 6.8469 -0.2222 Unstable Table 4.15: Control Variables with 8.5 times the Initial Reactive Load at Bus 5. load before bus 1 collapses and why in one case a solution is not found for bus 3, a generator bus (Table 4.14). 4.5 Load and Voltage Variations Now since the effects of both the state variables have been individually looked at it would be interesting to examine how variations of both of them simultaneously would affect the voltage stability and the catastrophe model. It is an obvious conclusion from the proceeding sections that lower voltage and higher reactive loads bring a system closer to voltage instability. Table 4.16 shows how the two state variables can bring the system quickly to collapse if both are changing. The higher the load becomes the higher the critical voltage point (the voltage value for which the system 37 Chapter 4: Examination of the General Model on a Test System Reactive Load Factor Critical Voltage Factor for Collapse (Bus 4 in All Cases) 1.0 0.826 2.0 0.823 3.0 0.873 4.0 0.929 5.0 0.979 5.5 1.0 6.0 1.022 7.0 1.06 Table 4.16: Effects of both Load and Voltage on Voltage Stability. becomes unstable). It should be noted that the factors in Table 4.16 are uniform over all buses and collapse occurs at bus 4 (the weakest bus) in all cases. 38 Chapter 5: Conclusions Chapter 5 Conclusions It has been shown i n this thesis that catastrophe theory can be used to predict voltage collapse i n power systems. A general swallowtail catastrophe model has been developed and applied to a realistic power system. F r o m this model the voltage stability status for any given system condition can be predicted. The model is based on finding the control variables of a swallowtail catastrophe based on the system equations for a specific bus and determining i f they l ie within a unique stability region i n the control space. This is determined for one bus at a time where the single state variable of the swallowtail catastrophe is the voltage at the particular bus. Prior to this al l other bus voltages are solved i n terms of this voltage. However because of the nonlinearity and the possibility of interconnections i n the system some approximations need to be done to obtain this single state variable. In some cases the approximations involve using initial voltage values for some of the bus voltages far away from the particular bus. But since during voltage collapse, voltage changes obviously take place, and although those changes might be minor at buses far away from the point o f collapse the effects o f these approximations need to be studied further. The model was first tested on a simple 2-bus power system which had some previous voltage stability results available. The results from the catastrophe model were compared with the known results and when proved to be accurate the model was extended to a more realistic interconnected power system. F o r a six-bus test system the model predicted voltage collapse at high loads and/or l o w voltages as would be expected. N o general methods are available to predict voltage collapse quickly and easily i n this way, so the use o f catastrophe theory is novel and very promising. 39 Chapter 5: Conclusions These results are very encouraging for further work using catastrophe theory for voltage stability studies especially as a possible online stability tool. More research lies ahead before this will be realized because the problem of voltage stability is very complex. Studies need to be done on very large power systems where voltage stability data exists. One such system is the so-called New England power system which has been designated as a voltage stability study system. Also it is of interest to try to formulate dynamic voltage stability as well, but that would include the power angles as a state variable and therefore increase the complexity of the system. This would bring us closer to a general voltage and power stability solution for power systems. 40 Appendix A The Catastrophe Theory and the Swallowtail Catastrophe A.1 Introduction Catastrophe theory [21] grows where algebra, calculus, and topology meet each other, and is concerned with the study of real-valued functions of several real variables. As a part of mathematics, catastrophe theory is a theory about singularities. When applied to scientific problems, therefore, it deals with the properties of discontinuities directly, without reference to any specific underlying mechanism. The theory attempts to study how the qualitative nature of the solutions of equations depends on the parameters that appear in the equations. Elementary catastrophe theory is the study of how the equilibria ^j{CQ) of V($j\ Ca) change as the control parameters Ca change. Catastrophe theory shows that the number of qualitatively different configurations of discontinuities that can occur depends upon the number of control variables, which are generally few, and not upon the number of state variables, which may generally be many. If wc have a family of functions V : F x C ^ R (A.56) where J 1 is a manifold, usually Rn, and C is another manifold, usually Rr. Rn is the state space and RT is the control space. The state space has dimension n and the control space has dimension r. A manifold is a term to indicate a high-order surface (hypersurface), e.g., a 1-dimensional (lst-order) manifold is a curve, a 2nd-order manifold is a contour, etc. The catastrophe manifold, N, is the subset Rn x RT denned by V,V c(x) = 0 (A.57) where Vc(x) = V(x,c) and x is the state variable. This is the set of all critical points of all the potentials Vc in the family V and V x is the partial derivative with respect to x. 41 N o w we find the singularity set, S, which is the subset of the manifold, N, that consists of all singular points o f V. These are the points at which V xV c(x) = 0 (A.58) and V 2V c(x) = 0 (A.59) The singularity set 5 is then projected down onto the control space RT to obtain the bifurcation set B. The bifurcation set is the image o f the catastrophe manifold N i n the control space C. The bifurcation set B provides the projection o f the stability region of all possible stable points o f V i n terms o f the control variables, which usually represent the system parameters. A.2 The Elementary Catastrophes F o r systems with less than five control variables, there are seven distinct types o f catastrophes [21] and no more than two state variables are involved in any of these. These catastrophes are called the elementary catastrophes. • r = l , rt=l. The fold catastrophe. V F ( u , x ) = U X + i x 3 (A.60) • r=2, n=l. The cusp catastrophe. V c (u ,x) = u i x + ^u 2 x 2 + ix 4 (A.61) • r=3, n = l . The swallowtail catastrophe. V s (u ,x) = m x + ^u 2 x 2 + ^ x 3 + Ix* (A.62) I o o 42 r=3, n=2. The hyperbolic umbilic catastrophe. V D ( U , X ) = UiXi + U 2 X 2 + U3X1X2+X1 + Xj (A.63) . r=3, n=2. The elliptic umbilic catastrophe. VE(u,x) = uixi + u2x2 + un(x2 + x2)-t-Xj - 3xix2 (A.64) • r=4, n=l. The butterfly catastrophe. VD(u,x) = uix + l U 2 x 2 + \x* + ±x4 + ^ x6 (A.65) z o 4 0 • r=4, «=2. The parabolic umbilic catastrophe. Vn(u,x) = uixj + u2x2 -(- U3X 2 + U4X 2+x 2x 2 + x2 (A.66) The set in Rn x RT where the differential of Vu is zero turns out to be a differential manifold for each one of the elementary catastrophes. It is sometimes called the catastrophe manifold, and in other contexts the equilibrium set. This set is then projected onto RT. The collection of the critical points of the projection is of particular interest. It is sometimes called the catastrophe set. The image of the catastrophe set under the projection is the set of all critical values of the projection. It is sometimes called the bifurcation set. In this thesis, catastrophe theory has been applied to voltage stability analysis in power systems for the purpose of finding stable and unstable regions directly from the system equations. The elementary catastrophe chosen for the study was the swallowtail. Some reasoning was behind this choice. Firstly, our system equations only have one state variable in terms of voltage stability, namely 43 voltage. Therefore the hyperbolic, elliptic and the parabolic umbilics were eliminated since they all have two state variables. This left us with a choice of anything from a one-dimensional to a four-dimensional control space. The fold was immediately thought to be too simple to produce accurate enough results. For similar reasons the swallowtail prevailed over the cusp, i.e., without getting to complicated it would be more accurate. Finally, the butterfly was considered to be complicated, simply for the reason that the control space would be 4-dimensional and hard to visualize. Also some previous studies on power systems had been done using the swallowtail catastrophe [13],[14]. A3 The Swallowtail Catastrophe The derivative of the potential function Vs(u,x) = uix + ^u2x2 + ^u3xn + ^x5 (A.67) z o o with respect to the state variable x is /5(u,x) = = 11!+ u2x + uax2 + x4 (A.68) The catastrophe manifold Ms is the set Ms = {(u,x) € R n x R : /5(u,x) = 0} (A.69) The catastrophe set Cs is the set of all points in Ms for which d e t ef^x) _ d e t , + 2 ^ + ^ = u 2 + 2 u ^ x + _ 0 ( A > 7 0 ) ox v ' Cs is therefore the solution set of the system ui + u2x + U.-JX2 + x4 = 0 (A.71) u2 + 2u3x + 4xn = 0 Eliminating the second term in the first equation yields an equivalent system of equations: ui - unx2 - 3x4 = 0 (A.72) u 2 + 2u3x + 4x3 = 0 44 From this we obtain C s = {(st2 + 3t4, -2st - 4t3, s, t) : (s, t) e R2} (A.73) The bifurcation set Bs is the projection of Cs onto it*3: B s = {(st2 + 3t4, -2st - 4t3, s, t) : (s, t) G R2} (A.74) We can have a closer look at Bs by studying its cross sections with the planes E? x {s}. The projection onto B? of the cross section is the image of the real line under the mapping As : R —• R2,A8(t) = (st2 + 3t4, -2st - 4t3) (A.75) that is, BB n (R2 x s) = A6(R) x {s} (A.76) Since X"(—t) = (X\(t), — A|(t)) , the image X"(R) is symmetric in theui axis. \"(R) is a curve in R?. Its slope can be computed from (Xs)'(t) = 2(s + 6t 2 )(t , - 1 ) . We obtain t^O (A 77) The curve turns left. The direction of movement along A*(.ft) is given by m . W i - 5 f c * ! ! l , 11 (A7S1 m~wm'^mw(t'-l) 45 Figure A.14: The Tail of the Swallow If s < 0, s + 6t2 changes sign from positive to negative at t\ = --^/(-s/6) and again from negative to positive at t2 = y/(-s/6). At all other points 6 is continuous. We have l i m S(t) = t - > t r / ( t i , - l ) = - l i m « ( t ) S 2 4s / - S l i m *(t) = . 1 (A.79) (t 2 ,-l)= - lim S(t) t - * t ^ s2 4s f-s The points A s ( t i ) and A s ( t 2 ) are cusp points. The tail of the swallow can now be seen in Figure A.M. If 5 = 0, s + 6t2 = 6f2, hence 6(t) is continuous, except for t = 0. Even there it has a limit: lim^(t) = (0,-1). The curve is shown in Figure A.15. If s > 0, s + 6t2 > 6t2, hence £(<) is continuous. X3 is an embedding. The curve is a one-dimensional submanifold. It is shown in Figure A. 16. 46 s = o Figure A.15: Manifold for s * * ( R ) _ _ _ _ S>0 Figure A.16: Manifold for s 47 Appendix B Program Code: General 6-Bus Case The swallowtail catastrophe model for voltage stability was programmed using the C programming language [22], The program was mostly built on function prototyping were each parameter and variables were calculated in their own function. Data was read from a file and then the program found the control variables and determined voltage stability from them. The program is a general program for the 6-bus example where any or al l of the buses could be examined for voltage stability. A general program for any power system would need to use arrays or pointers to greater extent and subprograms would be needed for the system solving prior to the catastrophe manipulation. A n example of the datafile read i n to the program follows the program listing. / * G E N E R A L S W A L L O W T A I L C A T A S T R O P H E P R O G R A M F O R A N A L Y S I N G V O L T A G E S T A B I L I T Y F O R A S P E C I F I C 6 - B U S E X A M P L E * / • i n c l u d e < m a t h . h > • i n c l u d e < s t d i o . h > / * G l o b a l V a r i a b l e s D e c l a r a t i o n s * / d o u b l e X[6] , R [ 6] , X c l , Xc5, Q [6] , V O [ 6] , t l , t5, d [ 6] , a l [ 6] , b e [ 6] , o m [ 6] , k s i [ 6] , j [ 7 ] , k [ 7 ] , l [ 7 ] , m [ 7 ] ; d o u b l e X c l , X c 5 , t l , t 5 ; i n t c h o s e b u s , c b u s , b u s , n b u s , c a s ; / * M a i n P r o g r a m * / m a i n () < / * F u n c t i o n D e c l a r a t i o n s * / 48 double Vt () ,Vtlow() ,T() , Tlow (), controlU () , controlV (), controlW (), paramAO () ,paramAl(),paramA2(),paramA3(),paramA4() , constantCl () , constantC2 (), constantC3 () ,constantC4(),constantC5(),constantC6() ,S0<),S1(),S2(),S3(),S4(),T0(),T1(),T2() ,T3<),T4(),P0(),P1(),P2(),P3(),P4(),R0() , R l ( ) , R 2 ( ) , R 3 ( ) , R 4 ( ) , p l O ( ) , p l l ( ) , p l 2 ( ) , p l 3 O , p l 4 ( ) , r l 0 < > , r l l ( ) , r l 2 ( ) , r l 3 ( ) , r l 4 0,k4(),14 0,ni4(),k3<),13(),m3 0 , j2 0,k2 0,12 0 f m 2 ( ) , k l O , l i t ) , ml () , alpha5 {) ,beta5 O , omega5 (), ksiS O ,alpha4 O,beta4 O,omega4 0 ,ksi4() ,alpha3(),beta3(),omega3(), ksi3{) ,alpha2 0 ,beta2 0 , omega2 0 , ksi2 () , a l p h a l ( ) , b e t a l ( ) , o m e g a l 0 , k s i l 0 ; /* V a r i a b l e s */ double contv.x; i n t i , c ; /* Function f o r Input of Data C a l l e d */ dataf (); /* Function f o r Choosing Voltage Collapse Bus */ buschoice(); f o r (c = 1; c <= cas; C + + ) ( i f (chosebus == 6) { cbus = c; ) e l s e { cbus = chosebus; i /* I n i t i a l i z i n g Arrays */ m[0] = 0.0;m[l] = 0.0;m[2] = 0.0;m[3] = 1[0] = 0.0;1[1] = 0.0;1[2] •= 0.0;1[3] = k[0] = 0.0;k[l] = 0.0;k[2] = 0.0;k[3] = j[0] = 0 . 0 ; j [ l ] = 0.0;j[2] = 0.0;j[3] = k[0] = -V0[0];k[6] •= -V0[0]; a l [ l ] = alpha! ();al[2) •= alpha2 () ; a l [3] al[ 5 ] = alpha5{);al[0]=0.0; b e [ l ] - b e t a l ().-be [2J = beta2 O ;be [3] = be[5] = beta5();be[0]=0.0; om[l] = omegal () ; om[2] «= omega2 () ;om[3] om[5] = omega5();om[0]=0. 0; k s i [ l ] = k s i l ( ) ; k s i [ 2 ] = k s i 2 ( ) ; k s i [ 3 ] k s i [ 5 ] = k s i 5 ( ) ; k s i [ 0 ] = 0.0; /* Control f o r Bus Constants C a l c u l a t i o n s */ f o r ( i = 5; i >= cbus + 1; i ~ ) < bus = i ; i f (bus == 5) { i f (bus == cbus + 1) { k[bus] = k4(); /* Loop for Each Bus */ /* Case of A l l Buses Examined */ /* Case of Single Bus Examined */ 0.0;m[4] = 0.0;m[5] •= 0.0;m[6]=0.0 0.0;1[4] = 0.0;1[5] = 0.0;1[6]=0.0 0.0;k[4] = 0.0;k[5J = 0.0;k[6]=0.0 0.0;j[4] = 0.0;j[5] = 0.0;j[6)=0.0 = alpha3 () ; a l [4] •= alpha4 0 ; beta3();be[4] = beta-4 O ; = omega3 () ;om[4] = omega4(); = ksi3 () ; k s i [4] = ksi4 () ; 49 e l s e ) 1 [ b u s ] m f b u s ] ) f k [ b u s ] 1 [ b u s ] m [ b u s ] } e l s e IK); k l Oi 1 1 0 ; m l ( ) ; { n b u s = b u s I f ( b u s == ( k [ b u s ] 1 [ b u s ] m [ b u s ] ) e l s e ( j [ b u s ] k [ b u s ] 1 [ b u s ] m [ b u s ] ) + 1 ; c b u s + 1 ) k 3 ( ) , 1 3 0 ; m3 O ; j 2 ( ) k 2 ( ) 1 2 0 m 2 ( ) C o n s t a n t s o f T y p e F o u r C o n s t a n t s o f T y p e O n e C o n s t a n t s o f T y p e T h r e e C o n s t a n t s o f T y p e T w o f o r ( 1 = 1 ; i <= c b u s - 1 ; i + + ) 1 ) b u s = i ; i f ( b u s • { i f ( b u s == c b u s - 1) el se k [ b u s ] 1 [ b u s ] m [ b u s ] ) < k [ b u s ] 1 [ b u s ] m [ b u s ] > e l s e e l s e j [ b u s ) k [ b u s ] 1 [ b u s ] m [ b u s ] ) n b u s = b u s i f ( b u s == ( k [ b u s ] 1 [ b u s ] m [ b u s ] ) k 4 () 14 0 m l 0 k l O; 1 1 0 ; m l ( ) ; - 1 ; c b u s -= k 3 ( ) ; = 1 3 0 , = n > 3 ( ) , j 2 ( ) ; k 2 ( ) ; 1 2 0 ; m2 () ; 1) C o n s t a n t s o f T y p e F o u r C o n s t a n t s o f T y p e O n e C o n s t a n t s o f T y p e T h r e e C o n s t a n t s o f T y p e T w o ) T y p e O n e B u s b e s i d e S l a c k B u s T y p e T w o G e n e r a l B u s T y p e T h r e e B u s b e s i d e V o l t a g e S t u d y B u s 50 Type Four Bus beside Voltage Study Bus and beside Slack Bus /* Output of Results */ p r i n t f ("\n") ; p r i n t f ( " C o l l a p s e at Bus Number %i\n",cbus); p r i n t f ( " C o n t r o l V a r i a b l e U: %4. 4f", controlU ( ) ) ; p r i n t f ( " ,V: %4 . 4f", controlV () ) ; p r i n t f ( " ,W: % 4 . 4 f c o n t r o l W ( ) ) ; p r i n t f ("\n") ; /* Checking f o r P o s i t i o n i n g With Respect to Manifolds */ i f (fabs(Vt()) > f a b s ( c o n t r o l V ( ) ) ) . { i f (controlUO > 0.0) < p r i n t f ( " S y s t e m Point Above Manifold, U > 0"); ) el s e ( i f (controlWO > 0.0 ii controlWO < 5. 0/4 . 0*pow (controlU () , 2. < pr i n t f ( " S y s t e m Point Inside Manifold"); } e l s e ( i f (controlWO > 5.0/4 . 0*pow (controlU () , 2. 0) ) < prin t f ( " S y s t e m Point Above Manifold, U < 0"); } e l s e { i f (fabs(Vtlowl) ) < fabs(controlV() ) ) { printf ( " S y s t e m Point Inside Manifold"); ) e l s e { p r i n t f ( " S y s t e m Point Below Manifold"); ) )}}) e l s e < p r i n t f ( " S y s t e m Point Outside Manifold"); ) p r i n t f ( " \ n \ n " ) ; ) ) /* Function f o r Input of Data */ dataf () { I* V a r i a b l e Declarations */ FILE *Y, * fopen () ; char junk[8], t ; double mV[6],mQ[6], t e s t ; f l o a t temp; i n t i , a ; i f ( (Y « fopen ("input.dat", V ) | == NULL) /* Opening D a t a F i l e */ 51 p r i n t f ( " E r r o r , Data F i l e NOT Found!"); ) e l s e < fo r ( i = 0; i <= 5; ++i) < fscanf (Y,"%s %s %f",junk, junk,stemp); d [ i l = (double) temp; ) f o r ( i = 0; i <= 5; ++i) { fscanf(Y,"%s %s %f",junk,junk,Stemp); V C[i] = (double) temp; } f o r ( i = 0; i <= 5; ++i) ( f s c a n f ( Y , " % s %s %f",junk,junk,stemp); Q[i] = (double) temp; ) f o r ( i = 0; i <= 5; ++i) ( fscanf(Y,"%s %s %f",junk,junk,Stemp); R[i] = (double) temp; ) f o r ( i = 0; i <= 5; ++i) ( fscanf(Y,"%s %s %f",junk,junk,Stemp); X [ i ] = (double) temp; ) fscanf(Y,"%s %s %f",junk,junk,Stemp); X c l = (double) temp; fs c a n f (Y,"%s %s %f",junk,junk,Stemp); Xc5 = (double) temp; fsca n f ( Y , " % s %s %f",junk,junk,Stemp); t l = (double) temp; fsca n f ( Y , " % s %s %f",junk,junk,stemp); t5 = (double) temp; f o r ( i = 0; i <= 5; ++i) ( f s c a n f ( Y , " % s %s %f",junk,junk,Stemp); m V [ i ) = (double) temp; V0[i) = V0[i]*mV[i]; ) f o r ( i = 0; i <= 5; ++i) < fscanf(Y,"%s %s %f",junk,junk,Stemp); mQ[i] = (double) temp; Q[i] = Q[i]*mQ[i]; } f c l o s e (Y) ; ); a = 1; ret u r n (a) ; } /* Function f o r Choosing Voltage Collapse Bus buschoice () < i n t a; pri n t f ( " W h i c h Bus (1-5), 6 i f a l l buses:"); scanf("%i",Schosebus); i f (chosebus •== 6) < cas = 5; /* A l l Buses */ ) e l s e /* Phase Angles */ /* Voltages */ /* Reactive Loads */ /* Line Resistances */ /* Line Reactances */ /* Capacitor */ /* Reactances. */ /* Transformer */ /* Ratio. */ /* Voltage Factors /* Load Factors */ 52 < cas =1; /* Only One Bus */ } a = 1; r e t u r n (a) ; ) /* Function f o r F i n d i n g Catastrophe Point Vt */ /* For the Manifold Comparising, a Swallowtail Manifold Point */ double Vt () { double v,pow(); v = 2.0*controlU()*T() + A . 0*pow (T () , 3. 0) ; re t u r n (v) ; ) /* Function f o r F i n d i n g the Catastrophe Point Vtlow */ /* For the Manifold Comparising, a Swallowtail Manifold Point */ double Vtlow() { double v, pow () ; v = 2.0*controlU()*Tlow() + A.0*pow(Tlow(),3.0); r e t u r n (v) ; ) /* Function f o r F i n d i n g the Catastrophe V a r i a b l e t */ /* To Find the Swallowtail Manifold Point Vt */ double TO ( double t, sqrt (), pow (); i f (controlWO < 0.0 ss controlUO > 0.0) ( t = 0.0; ) e l se ( t = sqrt (-controlUO/6.0 + sqrt (pow (controlU () , 2 . 0)/36. 0 + controlW () /3. 0) ) ; ) r e t u r n (t) ; ) /* Function f o r F i n d i n g the Catastrophe V a r i a b l e tlow */ /* To Find the Swallowtail Manifold Point Vtlow */ double Tlowl) ( double t, sqrt () ,pow () ; i f (controlWO < 0.0 it controlUO > 0.0) { t = 0.0; ) e l s e { 53 t - sqr t ( - c o n t r o l U ( ) /6 .0 - sqr t (pow (controlU (), 2. 0 ) /36 . 0 + c o n t r o l W ( ) / 3 . 0 ) ) ; } r e t u r n (t) ; ) / * Func t ion for F i n d i n g C o n t r o l V a r i a b l e U * / double contro lU( ) { double u; u = 1.0/paramA4 () * (paramA2 O - 3.0 4 paramA3 () *paramA3 () / (8.0*paramA4 () ) ) ; r e t u r n (u); ) / * F u n c t i o n f o r F i n d i n g C o n t r o l V a r i a b l e V * / double c o n t r o l V O double v; v = 1.0/paramA4 ()* (paramAl <) - paramA2 () *paramA3 () / (2. 0*paramA4 ()) + paramA3 () *paramA3 () *paramA3 () / (8.0*paramA4 () *paramA4 () ) )S 0 ; r e t u r n (v); ) / * F u n c t i o n for F i n d i n g C o n t r o l V a r i a b l e W * / double contro lWO ( double w; w = paramAO ()/paramA4 () - paramAl () *paramA3 () / (4 . 0*paramA4 () * paramA4()) + paramA2()*paramA3()*paramA3()/(16.0*paramA4()* paramA4()*paramA4()) - 3.0*paramA3()*paramA3()*paramA3()* paramA3 () / (256. 0*paramA4 () *paramA4 () *paramA4 () *paramA4 () ) S () ; r e t u r n (w) ; ) / * F u n c t i o n f or F i n d i n g Catastrophe Parameter Ao * / double paramAO () ( double aO; aO = c o n s t a n t C l f ) + constantC2 () *V0 [cbus] + 0.5* (POO + ROO) + constantC3() *V0[cbus]*V0[cbus]*S0() + constantC40*V0[cbus]*S00 + constantC5() *V0[cbus]*V0[cbus]*T0() + cons tantC6( )*V0[cbus ]*T0O; r e t u r n (aO) ; ) / * F u n c t i o n f o r F i n d i n g Catastrophe Parameter A l * / double paramAl0 ( double a l ; a l = constantC2() + 0.5*(P1() + R l () ) + const ant C3 () * (VO [cbus] *V0 [cbus ] *S1 () + 2*V0 [cbus] *S0 () ) + constantC4 ()* (V0[cbus]*Sl () + SO () ) + constantC5 () * (VO [cbus] *V0 [cbus] *T1 () + 2*V0 [cbus] *T0 () ) + constantC6()* (V0[cbus]*Tl () + T O O ) ; r e t u r n ( a l ) ; ) / * F u n c t i o n for F i n d i n g Catastrophe Parameter A2 * / double paramA2() 54 double a2; a2 = constantC2 ()/2.0 + 0.5MP2O + R2{)) + constantC3 {) * (VO [cbus] *V0[cbus]*S2() + 2*V0[cbus]*Sl () + (V0[cbus] + 1.0)*S0()) + constantC4 () * (VO [cbus] *S2 () + S l ( ) + SOO/2.0) + constantC5 () * (VO [cbus] *V0 [cbus] *T2() + 2*V0[cbus]*Tl () + (VOtcbus] + 1.0)*T0()) + constantC6()*(V0[cbus]*T2() + T l () + TOO/2.0); ret u r n (a2) ; ) /* Function f o r F i n d i n g Catastrophe Parameter A3 */ double paramA3 () < double a3; a3 •= constantC2 O/6.0 + 0.5* (P3() + R3()) + constantC3 () * (VO [cbus] *V0[cbus]*S3() + 2*V0[cbus]*S2() + (V0[cbus] + 1.0)*S1() + (V0[cbus]/3.0 + 1.0)*S0O) + constantC4()*(V0[cbus]*S3() + S2 0 + S1O/2.0 + SOO/6.0) + constantC5()*(V0 [cbus] *V0 [cbus] *T3 0 + 2*V0 [cbus] *T2 () + (V0[cbus] + 1.0)*T1() + (VOtcbus] /3.0 + 1.0)*T0O) + constantC6()*(V0[cbus]*T3() + T2 0 + T1O/2.0 + TOO/6.0); return(a3); 1 /* Function f o r F i n d i n g Catastrophe Parameter A4 */ double paramA4() { double a4; a4 = constantC20/24.0 + O.S*(P4() + R4 () ) + constantC3 () * (VO [cbus] *V0[cbus]*S4() + 2*V0[cbus]*S3() + (V0[cbus] + 1.0)*S2() + (V0(cbus]/3.0 + 1.0)*S1() + ((VOtcbus] + 1. 0) /12 . 0) *S0 () ) + constantC4 () * (VOtcbus) *S4 () + S3 0 + S2O/2.0 + S1O/6.0 + S0O/24.0) + constantC5 () * (VO [cbus ] *V0 [cbus] *T4 () + 2*V0[cbus]*T3() + (VOtcbus] + 1.0)*T2() + (VO [cbus]/3.0 + 1.0)«T1() + ((V0[cbus] + 7.0)/12.0)*T0O) 4 constantC6 () * (VO [cbus] *T4 () + T3() + T2O/2.0 + T1O/6.0 + TOO/24.0); ret u r n (a4) ; ) /* Function f o r F i n d i n g Constant CI */ /* A Constant For C o l l e c t i o n of Like 'Terms */ double constantCl() ( double CI; CI = -1.0/2.0*(k[cbus-l)*al[cbus] + k[cbus+1]*ksi[cbus]); re t u r n (CI); ) /* Function f o r F i n d i n g Constant C2 */ /* A Constant For C o l l e c t i o n of Like Terms */ double constantC2() { double C2; C2 *= 2.0*om[cbus] + 2.0*be(cbus] - 1 [cbus-1] * a l [cbus] - 1 [cbus+l] * k s i [cbus] r e t u r n (C2) ; ) /* Function f o r F i n d i n g Constant C3 */ 55 / * A Constant For C o l l e c t i o n of L i k e Terms * / double constantC3() { double C3; C3 = ( 1 . 0 / 2 . 0 ) * l [ c b u s - l ] * l [ c b u s - 1 ] ; r e t u r n (C3) ; ) / * F u n c t i o n f o r F i n d i n g Constant C4 * / / * A Constant For C o l l e c t i o n of L i k e Terms * / double constantC4 () < double C4; C4 = ( 1 . 0 / 2 . 0 ) * 1 [ c b u s - 1 ] * k [ c b u s - l ] ; r e t u r n (C4); ) / * F u n c t i o n f o r F i n d i n g Constant C5 * / / * A Constant For C o l l e c t i o n of L i k e Terms * / double constantC5() ( double C5; C5 = ( 1 . 0 / 2 . 0 ) * l [ c b u s + l ] * l [ c b u s + U ; r e t u r n (C5) ; ) / * F u n c t i o n f o r F i n d i n g Constant C6 * / / * A Constant For C o l l e c t i o n of L i k e Terms * / double constantC6() { double C6; C6 = ( 1 . 0 / 2 . 0 ) * k [ c b u s + l ] * l [ c b u s + l ] ; return (C6) ; j / * F u n c t i o n f o r F i n d i n g Constant SO * / / * Constant from T a y l o r S e r i e s For Square Root * / double S0() < double SO; 50 = 1 . 0 / s q r t ( p l O ( ) ) ; r e t u r n ( S O ) ; ) / * F u n c t i o n f o r F i n d i n g Constant SI * / / * Constant from T a y l o r S e r i e s For Square Root * / double SI () { double SI ; 51 = - ( 1 . 0 / 2 . 0 ) * p l l ( ) / p o w ( p l 0 ( ) , 1 . 5 ) ; r e t u r n (SI) ; ) / * F u n c t i o n for F i n d i n g Constant S2 * / / * Constant from T a y l o r S e r i e s For Square Root * / 56 double S2() < double S2; 52 •= (1.0/ (2.0*pow(pl0 ( ) , 0 . 5 ) ) ) * (pow(pll () ,2.0) / (4.0*pow(pl0 0,2.0)) - p l 2 0 / p l O ( > ) ; return(S2); } /* Function f o r F i n d i n g Constant S3 */ /* Constant from Taylor Series For Square Root */ double S3 () ( double S3; 53 = (1.0/ (2.0*pow (plO () ,0.5) ) ) * (3.0*pll ()*pl2 0 /(2.0*pow (plO () ,2.0) ) - p l 3 ( ) / p l 0 ( ) - 5.0*pow(pllO,3.0)/(8.0*pow(plO(),3.0))); return(S3); > /* Function f o r F i n d i n g Constant S4 */ /* Constant from Taylor Series For Square Root */ double SI () < double S4,pow(); 54 = (1.0/ (2.0*pow(pl0 () , 0.5) ) ) * (3.0*pll () *pl3 () / (2 . 0*pow (plO (), 2 . 0) ) - p l 4 ( ) / p l 0 ( ) + 3.0*pow(pl2 () ,2.0) / (4.0*pow(pl0 () ,2.0) ) - 15.0*pow(pll () ,2.0) *pl2 0 / (8.0*pow(pl0 (), 3.0)) + 35.0*pow ( p l l () , 4.0) / (64.0*pow(pl0() , 4.0) ) ) ; return (S4) ; ) /* Function f o r F i n d i n g Constant TO */ /* Constant from Taylor Series For Square Root */ double TOO i double TO; TO = 1.0/sqrt (rlO () ) ; retu r n (TO); ) /* Function f o r F i n d i n g Constant TI */ /* Constant from Taylor Series For Square Root */ double T I 0 ( double TI; TI •= - ( 1 . 0 / 2 . 0 ) * r l l ( ) / p o w ( r l 0 ( ) , 1 . 5 ) ; r e t u r n ( T I ) ; } /* Function f o r F i n d i n g Constant T2 */ /* Constant from Taylor Series For Square Root */ double T2() ( double T2; T2 = (1.0/ (2.0*pow (rlO () , 0.5) ) ) * (pow ( r l l 0 , 2. 0) / (4 . 0*pow (rlO () ,2.0) ) - r l 2 0 / r l O ( ) ) ; r e t u r n (T2); ) 57 /* Function f o r Fi n d i n g Constant T3 */ /* Constant from Taylor Series For Square Root */ double T3 () ( double T3; T3 = (1.0/ (2.0*pow (rlO () , 0.5) ) ) * <3.0*rll () * r l 2 () / (2. 0*pow (rlO () , 2 . 0) ) - r l 3 ( ) / r l 0 ( ) - 5.0*pow(rll () , 3.0) / (8.0*pow(rl0 () ,3.0) ) ) ; retu r n (T3); } /* Function f o r F i n d i n g Constant T4 . */ /* Constant from Taylor Series For Square Root */ double T4 () ( double T4,pow () ; T4 = (1.0/ <2.0*pow (rlO () , 0.5) ) ) * (3. 0 * r l l () * r l 3 () / (2. 0*pow (rlO <) , 2. 0) ) - r l 4 ( ) / r l 0 ( ) + 3.0*pow(rl2 () , 2.0) / (4.0*pow (rlO () ,2.0) ) - 15.0*pow(rll (),2.0)*rl2() / (8. 0*pow (rlO (), 3.0) ) + 35.0*pow(rll 0,4.0) / (64 . 0*pow (rlO () , 4 . 0) ) ) ; retu r n (T4) ; ) /* Function f o r F i n d i n g Constant P0 */ /* Constant from Taylor Series For Square Root */ double P0() < double P0,sqrt 0; P0 = sqrt (plO () ) ; retu r n (P0) ; } /* Function f o r Fin d i n g Constant PI */ / * C o n s t a n t f r o m Taylor Series For Square Root */ double PIO < double P I , s q r t ( ) ; PI = ( 1 . 0 / 2 . 0 ) * p l l ( ) / s q r t ( p l 0 ( ) ) ; r e t u r n (PI) ; ) /* Function f o r F i n d i n g Constant P2 */ /* Constant from Taylor Series For Square Root */ double P2 () ( double P2,pow(); P2 = (1.0/2.0)*pow(pl0(),0.5)*(pl2()/pl0() - powtpll (),2.0)/ (4.0*pow(pl0 () , 2.0) ) ) ; re t u r n (P2) ; } /* Function f o r F i n d i n g Constant P3 */ /* Constant from Taylor Series For Square Root */ double P3 () { double P3,pow(); P3 = (1.0/2.0) *pow(pl0 (), 0.5) * (pl3 O/pl0 () - p l l () *pl2 () / (2. 0*pow (pi 0 () , 2. 0) ) + pow ( p l l (),3.0) / (8. 0*pow(pl0 (),3.0) ) ) ; 58 r e t u r n (P3) ; ) /* Function f o r F i n d i n g Constant P3 */ /* Constant from Taylor Series For Square Root */ double P4() { double P4, pow () ; P4 = (1.0/2.0) *pow(plO (), 0.5) * (pl4 0/plO () - p l l0 * p l 3 ( ) / ( 2 . 0 * p o w ( p i 00 , 2 - pow (pl2 () , 2. 0) / (4 . 0*pow (plO () , 2. 0) ) + 3.0*pow ( p l l () ,2.0) *pl2 () / (8.0*pow(pl0 () , 3.0) ) - 5.0*pow(pll () , 4.0) / (64.0*pow (plO () , 4 . 0) ) ) ; retu r n (P 4); ) /* Function for F i n d i n g Constant R0 */ /* Constant from Taylor Series For Square Root */ double R0 () < double R0,pow(); R0 = pow(rl0 () , 0.5) ; retu r n (R0) ; ) /* Function f o r Find i n g Constant Rl */ /* Constant from Taylor Series For Square Root */ double Rl () ( double Rl,pow(); Rl = (1.0/2.0) * r l l ()/pow(rl0 () , 0.5) ; r e t u r n ( R l ) ; ) / * F u n c t i o n f o r Fi n d i n g Constant R2 */ /* Constant from Taylor Series For Square Root */ double R2() ( double R2, pow () ; R2 = (1.0/2.0) *pow(rl0 () , 0.5) * (rl2 O / r l O () - pow(rll (),2.0)/ (4.0*pow (rlO () , 2. 0) ) ) ; return (R2) ; } /* ' Function f o r F i n d i n g Constant R3 */ /* Constant from Taylor Series For Square Root */ double R3 () ( double R3,pow () , sqrt () ; R3 = (1.0/2.0)*sqrt (rlOO )* ( r l 3 ()/rlO () - r l l 0 * r l 2 O / (2 . 0*pow (rlO () , 2 . 0) + pow ( r l l () ,3.0) / (8.0*pow(rl0 () ,3.0) )) ; retu r n (R3) ; ) /* Function f o r F i n d i n g Constant R3 */ /* Constant from Taylor Series For Square Root */ double R4 () f 59 double R4,pow () , sqrt () ; R4 = (1.0/2.0) *s q r t (rlD() )* (rl4 ()/rlO () - rl1()*rl3()/(2.0*pow(rlO<),2.0)) - pow(rl2 () ,2.0) /(4.0*pow(rl0 () ,2.0) ) + 3.0*pow(rll0,2.0)*rl2()/(8.0*pow(rlO(), 3.0)) - 5.0*pow ( r l l () , 4.0) / (64.0*pow (rlO () , 4.0) ) ); return(R4); ) /* Function f o r F i n d i n g Constant plO */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double plO () < double p0,pow(); pO = pow(1[cbus+1],2.0)*pow(V0[cbus] ,2.0) + 2.0*k[cbus+l]*l[cbus+l]*V0[cbus] + m[cbus+l] + pow(k[cbus+1],2.0); return(pO); ) /* Function f o r F i n d i n g Constant p l l */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double p l l () ( double p i , pow (); p i = 2.0*(pow(l[cbus+l],2.0)*V0[cbus] + k[cbus+1]*1[cbus+1]) ; return (pi) ; ) /* Function f o r F i n d i n g Constant pl2 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double pl2() < double p2,pow(); p2 = pow(1[cbus+1],2.0)*(V0[cbus] + 1.0) + k[cbus+1]*1[cbus+1]; return (p2) ,-) /* Function f o r F i n d i n g Constant pl3 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double pl3() { double p3,pow () ; p3 = pow(l[cbus+l],2.0)*(V0[cbus]/3.0 + 1.0) + k[cbus+1]*1[cbus+1]/3.0; ret u r n (p3) ; ) /* Function f o r F i n d i n g Constant pl4 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double pi4 () { double p4,pow(); p4 = pow(1[cbus+1],2.0)*(V0[cbus] + 7.0)/12.0 + k[cbus+1]*1[cbus+1]/12. 0; return(p4); } /* Function f o r F i n d i n g Constant plO */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ 60 double rlO() < double rO,pow(); rO = pow(l [cbus-1] ,2.0)*pow(V0[cbus],2.0) + 2.0*k[cbus-1]*1[cbus-1]*V0[cbus] + m[cbus-l] + pow(k[cbus-1],2.0); r e t u r n (rO); } /* Function for F i n d i n g Constant r l l V /* Parameter of C o l l e c t e d Terms Inside Square Root */ double r l l () i double rl,pow(); r l = 2.0*(pow(l[cbus-l],2.0)*V0[cbus] + lc [cbus-1 J *1 [cbus-1] ); re t u r n ( r l ) ; ) /* Function f o r Fin d i n g Constant r l 2 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double r l 2 () { double r2,pow (); r2 = pow(l[cbus-l],2.0)*(V0[cbus] + 1.0) + k [cbus-1 ] *1 [cbus-1 ] ; r e t u r n ( r 2 ) ; ) /* Function f o r F i n d i n g Constant r l 3 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double r l 3 ( ) { double r3,pow(); r3 = pow(l[cbus-l],2.0)*(V0[cbus]/3.0 + 1.0) + k[cbus-1]*1[cbus-1]/3.0; r e t u r n (r3); ) /* Function f o r F i n d i n g Constant r l 4 */ /* Parameter of C o l l e c t e d Terms Inside Square Root */ double r l 4 () i 1 double r4,pow(); r4 = pow(1[cbus-1],2.0)*(V0[cbus] + 7.0J/12.0 + k[cbus-1]*1[cbus-1]/12.0; r e t u r n (r 4) ; ) /* Constants for Type 4 */ double k4 () I double k4; k4 = al[bus]*V0[0]/(be[bus] + om[bus]); r e t u r n (k4) ; ) double 14 0 ( double 14; 14 = k s i [bus] / (befbus] + om[bus]); 61 r e t u r n (14) ; ) double m4() { double m4; m4 = 4.0*Q[bus]/(be[bus] + omfbus]); return(m4); ) ,'* Constants f o r Type 3 */ double k3 () < double k3; k3 = (al[bus]*j[nbus]/(l[nbus] + m[nbus]*V0[bus]))/(be[bus] + om[bus]); re t u r n (k3) ; ) double 13 0 < double 13; 13 = ksi[bus]/(be[bus] + om[bus]); return(13); } double m3() { double m3; m3 = (4.0*(Q[bus] - a l [bus] *k [nbus] / (1 [nbus] + m[nbus] *V0 [bus] ) ) ) / (be [bus + om[bu s]) ; return(m3); ) /* Constants for Type 2 */ double j2() ( double j2; j2 = k s i [bus] *m[nbus]*VO [bus] *V0 [bush-return (j2) ; ) double k2 () { double k2; k2 = C[bus]*1[nbus] + (2.0*m[nbus)*V0[bus] + 1[nbus])*pow(VO[bus],2.0) *(be[bus] + om[bus]); re t u r n (k2) ; ) double 12() { double 13; 13 - VO[bus]*(be[bus] + om[bus])*(3.0*m[nbus]*V0[bus] + 2.0*l[nbus)) + al[bus]*k[nbus] - m[nbus] *Q[bus] ; re t u r n (13) ; ) 62 d o u b l e m2 () { d o u b l e m2; m2 = 2 . O ' k s i [ b u s j * m [ n b u s ] * V 0 [ b u s ] + k s i [ b u s ] * 1 [ n b u s ] ; r e t u r n ( m 2 ) ; > / * C o n s t a n t s f o r T y p e 1 * / d o u b l e k l () { d o u b l e k l ; k l = Q [ b u s ] + (om[bus) + b e [ b u s ] ) * p o w ( V O [ b u s ] , 2 . 0 ) ; r e t u r n ( k l ) ; ) d o u b l e 11 ( ) { d o u b l e 1 1 ; 11 = a l [bus]*V0[0] + 2 . 0 * ( o m [ b u s ] + be [bus] ) *V0 [bus ] ; r e t u r n (11) ; ) d o u b l e m l ( ) { d o u b l e m l ; ml = k s i [bus] ; r e t u r n ( m l ) ; ) / * C o n s t a n t s f r o m R e a c t i v e Power E q u a t i o n a t B u s 5 * / d o u b l e a l p h a 5 () { d r v j V l e a l p h ; a l p h = ( X [ 5 ] * c o s ( d [ 5 ) - d [ 0 ] ) + R [5 ] * s i n (d [ 5 ] - d [0] ) ) / (pow (R [5] , 2 . 0) + p o w ( X [ 5 J ,2.0) ) ; r e t u r n ( a l p h ) ; ) d o u b l e b e t a 5 ( ) ( d o u b l e b e t ; b e t = - X [ 5 ] / ( p o w ( R [ 5 ] , 2 . 0 ) + pow (X [5] , 2 . 0) ) - 1 . 0 / X c 5 ; r e t u r n ( b e t ) ; ) d o u b l e o m e g a 5 0 ( d o u b l e omeg; omeg = 1. 0/(pow ( t 5 , 2 . 0) *X [ 4 ] ) ; r e t u r n ( o m e g ) ; 1 d o u b l e k s i 5 () ( d o u b l e a l p h ; a l p h = - ( c o s ( d [ 5 ] - d [ 4 ] ) ) / ( 2 . 0 ' t 5 * X [ 4 ) ) ; r e t u r n ( a l p h ) ; ) 63 /* Constants from Reactive Power Equation at Bus 4 */ double alpha4() ( double alph; alph =