NUMERICAL INSTABILITIES IN POWER SYSTEM TRANSIENT SIMULATIONS By Antonio Emilio Angueth de Araujo B. Sc. (Electrical Engineering) Universidade Federal de Minas Gerais, Brazil M. Sc. (Electrical Engineering) Universidade Federal de Minas Gerais: Brazil A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June 1993 © Antonio Emilio Angueth de Araujo, 1993 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. Electrical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: `3 9 3 Abstract In power system simulations, control systems act upon the electric power system in a variety of ways. There is, therefore, a necessity of simulating the control system equations together with the power system equations. This is the case not only in stability studies, but also in the simulation of electromagnetic transient phenomena. Because of the different structure of the two sets of equations, computational efficiency can be improved if both sets of equations are solved separately, with the result from one set of equations showing up in the other set of equations one or more time steps later. This procedure, however efficient, can lead to numerical instability. In the simulation of the control system itself, the modelling of nonlinearities in closed loops poses some problems as well. The most efficient way of handling a nonlinear block is to open the closed loop and use, instead of the current value of the variable, a predicted value of the feedback path. This can also lead to numerical instabilities. In this thesis a simultaneous solution of the control and power system equations is proposed. For the cases of nonlinearities in closed loops, an iterative solution is suggested. The implementation of the method is carried out in the Electromagnetic Transients Program (EMTP). Several comparisons between the proposed and other methods are performed. In cases in which the analytical solution is not available, the comparison is made against the results of the TAGS package ("Transient Analysis of Control Systems") which ii is integrated with the EMTP version of DCG/EPRI ("Development Coordination Group/Electric Power Research Institute"). It is shown through simple, but typical, case studies that the proposed method is not only accurate and efficient, but also numerically stable, thus representing a significant improvement to the ability of the EMTP to reliably simulate electromagnetic transients in power systems when control systems are involved. iii Table of Contents ii Abstract^ List of Tables^ viii List of Tables^ viii List of Figures^ ix List of Figures^ xii 1 Introduction^ 1 2 Literature Review^ 4 2.1 Introduction ^4 2.2 EMTP-TACS Approach ^4 2.3 TACS Solution Method ^6 2.3.1 Transfer Functions - Linear Blocks ^6 2.3.2 Nonlinear Function Blocks ^8 2.3.3 Transfer Function with Limiters ^9 2.3.4 TACS-EMTP Interface ^ 11 2.4 Time Delays in the TACS-EMTP Approach ^ 12 2.5 Final Remarks ^ 14 iv 3 Time Delays and the Stability and Accuracy of the Simulation^15 3.1 Introduction ^15 3.2 Numerical Instability ^16 3.2.1 DC Circuit ^16 3.2.2 Arc modelling ^18 3.2.3 Sinusoidal Oscillator ^20 3.3 Stable Inaccurate Results 3.3.1 Internal Time Delays ^21 ^21 3.3.2 Converter Simulation - External Time Delays ^ 22 4 Simultaneous Solution of Power and Control System Equations ^25 4.1 Introduction ^25 4.2 Basic Considerations ^26 4.2.1 Elimination of the External Time Delay ^ 26 4.2.2 Elimination of the Internal Time Delays ^ 26 4.3 Actual Implementation ^ 28 4.3.1 Interface Between the Written Program and the EMTP . ^28 4.3.2 Capabilities of the Program ^33 5 Comparison Between Simultaneous and Non-simultaneous Approach 41 5.1 Introduction 5.2 Control System Case Studies ^ 5.2.1 Second-Order Linear System ^ ^41 42 42 5.2.2 Static and Dynamic Limiters in Open Loops ^ 42 5.2.3 FORTRAN Function in a Closed Loop ^ 44 5.2.4 Special Element Block ^ 5.3 Simulations of Control and Power Systems ^ 44 44 5.3.1 Arc Modelling ^ 44 5.3.2 Single-Phase Variable Load ^ 48 5.3.3 Three-Phase Variable Load ^ 53 5.3.4 Simple HVDC System Model ^ 64 ^70 5.4 Converter Simulations 5.4.1 Single-Phase Thyristor Rectifier Bridge ^72 5.4.2 Three-Phase Thyristor Inverter ^72 5.4.3 Remarks on Converter Simulation with the EMTP ^ 72 5.5 Final Remarks ^ 76 6 Conclusions^ 81 Bibliography^ 83 A Interpreter Basic Principles^ 86 A.1 Compilation of Arithmetic Expressions ^ 87 A.2 From Infix to Polish Notation ^ 89 A.3 Evaluation of a Reverse Polish expression ^ 100 B SS and EMTP Input Files^ 113 B.1 Control System Case Studies - SS Input Files Only ^ 113 B.1.1 Second-Order linear system (Section 5 2 1)^ 113 B.1.2 Static and Dynamic Limiters (Section 5.2.2) ^ 114 B.1.3 FORTRAN Function in a Closed Loop (Section 5.2.3) ^ 116 B.1.4 Special Element Blocks (Section 5.2.4) ^ 117 vi B.2 Control and Power System Case Studies - SS and EMTP Input Files 119 B.2.1 Arc Modelling (Section 5.3.1) ^ 119 B.2.2 Single-Phase Variable Load (Section 5.3.2) ^ 120 B.2.3 Step Load Variation (Section 5.3.3) ^ 122 B.2.4 Pulse Load Variation (Section 5.3.3) ^ 125 B.2.5 Load Variation with Thyristor Controlled Reactor (Section 5.3.3)128 B.2.6 HVDC System Simple Model (Section 5.3.4) ^ 135 ^ 138 B.2.8 Three-Phase Thyristor Inverter (Section 5.4) ^ 141 B.2.7 Single-Phase Thyristor Rectifier Bridge (Section 5.4) vii List of Tables A.1 Simple Precedence Function ^ 89 A.2 Transformation of a — b x c + d/5 into Reverse Polish ^ 90 A.3 Input and Stack Precedence Functions ^ viii 91 List of Figures 1.1 Interaction between power system and control system ^ 2 2.1 Typical diagram of a control system ^ 5 2.2 Transfer function block 7 2.3 Control system with linear transfer function blocks 2.4 Nonlinear function blocks 2.5 Static and dynamic limiters 2.6 Time delay between TACS and EMTP ^ 11 2.7 Possible simultaneous solution with TACS-EMTP ^ 13 3.1 Numerical instability in a DC circuit 16 3.2 Current and Voltage in the circuit of Figure 3.1 3.3 Arc modelling with EMTP-TACS ^ 18 3.4 z-transform domain analysis ^ 19 3.5 Numerical instability in the arc modelling 19 3.6 Block diagram and exact solution of^(t)^x(t) = 0 [1] ^ 20 3.7 Instability due to the internal time delay ^ 21 3.8 Inaccuracy due to the internal time delay [1] 22 3.9 Three-phase thyristor inverter current - TACS inaccurate result [2] 24 4.1 Partition of the augmented matrix 27 4.2 Compensation method in the EMTP ^ ^ ^ 9 ^ ^ ix 8 ^ ^ ^ ^ 10 17 29 4.3 Subroutine CONNEC and its interface with the EMTP ^ 30 4.4 Control actions upon the power system ^ 31 4.5 Matrix formation inside CONNEC ^ 32 ^36 4.6 Controlled integrator block 5.1 Second-order linear system: (a) Diagram; (b) Solution for a unit step ^43 input. 5.2 Static limiter: (a) Diagram; (b) Simulation results. ^45 5.3 Dynamic limiter: (a) Diagram; (b) Simulation results ^46 5.4 Nonlinearity in a closed loop: (a) Diagram; (b) Simulation results [1] ^47 5.5 Special element blocks - Diagram and simulation results (source = cos wt, w = 377 rad/s and time delay = 2 ms). ^ 47 5.6 Arc modelling - superposition of the analytical and the SS simulation ^48 results 5.7 R-L load representation; c = R/L ^ 5.8 Active power variation ^ 5.9 Block diagram - Single-phase variable load [1] 49 50 ^51 5.10 Simulation Results: (a) TACS solution [1]; (b) SS results ^52 5.11 Sample power system [1]. ^ 54 5.12 Diagram of the Thevenin reduced system [1]. ^ 54 5.13 Load modelling - step variation [1] ^56 5.14 (a) Load variation - step variation, (b) TACS simulation result [1], (c) SS simulation result. 5.15 Load modelling - pulse variation [1]. ^57 ^58 5.16 (a) Load variation - pulse variation, (b) TACS simulation result [1], ^59 (c) SS simulation result. 5.17 Thyristor controlled reactor: (a) Circuit diagram; (b) Circuit model; (c) TCR s-domain model ^ 5.18 Basic TCR voltage controller - 61 ic/L modelling [1] 5.19 TCR and its control system - a simple model ^62 ^62 5.20 SS block diagram. ^63 5.21 Load variation. ^64 5.22 Load variation without TCR: (a) TACS simulation [1]; (b) SS simulation. 65 5.23 Load variation with TCR: (a) TACS simulation [1]; (b) SS simulation. 66 5.24 Power system with a DC line [1]. ^ 67 ^68 5.25 Id and a modelling [1] 5.26 EMTP-SS simulation diagram [1]. ^ 68 5.27 (a) Voltage tracking system; (b) DC system equations; (c) AC current ^69 injection [1] 5.28 IDC, ALPHA and IACF comparisons; (a) TACS simulation [1]; (b) SS ^70 simulation. 5.29 Voltage at bus 1 (BUS1) comparison; (a) TACS simulation [1]; (b) SS ^71 simulation. ^73 5.30 Single-phase thyristor rectifier. 5.31 Single-phase rectifier simulations - voltage vd: (a) TACS results; (b) ^74 SS results 5.32 Closing and opening of EMTP switches. ^ 75 5.33 Backtracking procedure. ^77 5.34 Inverter simulation: instability due to thyristor misfiring ^78 xi Acknowledgements To my wife Angela, for her having put aside her own career and interests to kindly accompany me and, with love, make my undertaking possible. To my children Thiago and Henrique, who made my life a little more complicated but full of love and fun. To Dr. H.W. Dommel and Dr. J.R. Marti for their invaluable help, direction and advice. To Dr. Nelson Santiago for his help during the period of my thesis writing. To my colleagues for making my stay at UBC memorable. Finally, to CNPq (Brazilian Research Council), to Universidade Federal de Minas Gerais and to The University of British Columbia for their financial support. xii Chapter 1 Introduction During the 1960's and 1970's the modelling and simulation of power systems transient phenomena reached a high degree of sophistication. Numerical techniques, as well as models of various power system components, were developed for digital transient simulation software packages. With the increased use of HVDC transmission systems, and the wide application of other power electronic devices in power systems, it became necessary to simulate the control systems acting upon the power system in addition to the power system itself. This is very important in cases where the dynamic interactions between the control and the power system are of major concern. In general, the control system acts upon the power system in three different ways: • issuing commands to open or close switches and thyristors; • controlling voltage sources; or • controlling current sources. The decisions made by the control system are based on the information which it receives from the power system, such as node or branch voltages, branch currents or state of the switches. Figure 1.1 depicts the general situation. 1 Chapter 1. Introduction^ Power System 2 switch status voltage current '-7 switch command voltage current Control System Figure 1.1: Interaction between power system and control system The power system component models are relatively easy to incorporate into nodal network equation representations. Control system components are not as easy to model. They do not fit into the nodal network representation, and the associate matrices are not symmetric. The easiest and most efficient way of handling these differences, is to simulate the power system equations separately from the simulation of the control system equations. Both sets of equations are solved separately, and then interfaced with each other through the exchange of information about the values of the parameters which they share. This introduces time delays between both solutions and makes the solution procedure non-simultaneous. The simulation of the control system itself may require some special treatment as well. Whenever a nonlinearity exists in a closed loop in the control system, there is an efficient and non-simultaneous procedure to overcome this difficulty - opening the loop and using the past value as the actual value of the output of the feedback path. Chapter 1. Introduction^ 3 Non-simultaneity is the basis for the efficiency of the solution method outlined above. However, this procedure can lead to two distinct and important problems inaccuracy and numerical instability. These two problems are addressed in this thesis. A general simultaneous solution method of the control and power system equations is proposed. Comparisons between the proposed and other methods, showing its accuracy and stability, are performed. The thesis is divided into six chapters as outlined below: • Chapter One - Introduction • Chapter Two - Literature Review ^ An overview of the relevant technical literature is presented. • Chapter Three - Effects of Time Delays on the Stability and Accuracy of the Simulation ^ The main objective of this chapter is to identify the most important effects of time delays (non-simultaneity) on the stability and accuracy of the solution. • Chapter Four - Simultaneous Solution of Power and Control System Equations ^ The proposed method is described and the implementation of the method is discussed. • Chapter Five - Comparison between simultaneous and non-simultaneous approaches ^ The objective here is to show the main advantages and limitations of the simultaneous approach when compared with the other approaches. This is done through several case studies. • Chapter Six - Conclusions Chapter 2 Literature Review 2.1 Introduction Several current digital computer programs for power system transient analysis are known as "Electromagnetic Transients Program". Most of them are variations of the basic EMTP [3]. See reference [4] for a complete description of the various versions. Control system simulation capabilities were originally added to the BPA EMTP (Bonneville Power Administration EMTP). Although other programs not directly related to the EMTP exist which are able to simulate control and power systems, this chapter will cover only the EMTP-related literature, because of the wide-spread use of the EMTP compared to other programs. 2.2 EMTP-TACS Approach One of the first works on the simulation of control and power systems using general purpose programs was done by L. Dube and H. Dommel [5]. Their technique was incorporated into a digital computer program called TACS (Transient Analysis of Control Systems), which became part of the EMTP. TACS is still used today in the DCG/EPRI EMTP. The main objective at that time was the simulation of HVDC converter controls. 4 Chapter 2. Literature Review^ From internal source 5 To EMTP G(s) Kl/s Vcos(180- P) From EMTP Sensor Figure 2.1: Typical diagram of a control system Some important features of TACS are: • its code is independent of the EMTP, and interfaces with the EMTP with one time step delay; • it can accept arbitrary interconnections of a set of control system building blocks; • it solves the control system equations by implicit integration with the trapezoidal rule. Figure 2.1 shows a typical diagram of a control system to be simulated by the EMTP-TACS approach, where the control receives some information from the network, performs some actions and then gives information back to the network. Although there have been many additions to the TACS set of blocks over the years [6], the basic blocks used to simulate any control system are: • transfer functions — s-domain rational polynomial functions which are the most basic building blocks in a control system; • sources — time dependent sources, such as DC, pulse, step, etc; ^ 6 Chapter 2. Literature Review^ • Fortran expressions — a general function given by a Fortran statement which describes the relationships between the inputs and the output of a block; • devices — preprogrammed subsystems of control systems, such as relay operated switches, sample and track circuits, etc. 2.3 TACS Solution Method 2.3.1 Transfer Functions — Linear Blocks A general transfer function block is shown in Figure 2.2. It describes a relationship between the output X and the inputs Ui as a rational function G(s) = N(s)/D(s) with a gain K and m < n. Such an n-th order transfer function is transformed into a system of first-order differential equations which are integrated with the trapezoidal rule. After eliminating internally created state variables for second and higher order derivatives, a single equation relating the output to the inputs is formed: cx(t) = K • d • u(t) -+ hist(t — At) , (2.1) where u is the sum of all input variables. The coefficients c and d are calculated using the recursive formula: (-2) i {^(i 2 \i+i d . +1 i^ (At) ' • • • + (71 ( rdn}, ) At where ei ) is the binomial coefficient, and n^2 • co = E() 2 D i i . 0 At where c = co . The formulas for d are found replacing D by N in the above formulas. Chapter 2. Literature Review^ 7 Ul 0 U U2 +N s+.+Ns 1 K D + D s + + Dn s O l X n U3 Figure 2.2: Transfer function block The calculation of the history term hist(t — At) is not simple and the reader is referred to Chapter 4 and reference [7] for details. After discretizing all transfer functions, it is possible to build a matrix equation (A)(X) = (hist). To clarify this point, assume the system of Figure 2.3 is to be simulated [7]. Applying Eq. 2.1 to the blocks of this figure, the related matrix equation will be Ca —Kbd b —K c d, 0 Ka da 0 c, 0 Ka d a 0 0 Cd 0 Cb — 0 K add Kada 0 0 0 0 Kbdb / xi x2 X3 0 x4 0 U1 \ (hist a histb hist histd (2.2) U2 If u i and u 2 are known sources at each time step, the unknowns are x l , x 2 , x 3 and x 4 . The previous equation can be written as follows: (A x )(x) = (hist) — (A xu )(u),^ (2.3) Chapter 2. Literature Review^ 8 -U2 Ul X1 a X4 b d Figure 2.3: Control system with linear transfer function blocks where Ca Ka d a Ka d a 0 —Kbdb 0 0 Cb —Kc d, 0 c, 0 0 0 cd —Kddd ( A xe =-and ASU ( —Ka da 0 0^Kb db 0^0 0^0 The Equations 2.3 can now be solved for the unknowns at each time step. Triangular factorization and ordering techniques with a sparsity storage scheme are used in the solution. The factorization is performed once, and at each time step there is one forward solution and one back substitution [5]. 2.3.2 Nonlinear Function Blocks Nonlinear function blocks are not included in the previous simultaneous solution of the whole system. Instead, they are solved sequentially as soon as their input signals become available at some stage in the back substitution mentioned above. 9 Chapter 2. Literature Review^ System —4 Section A N.L. System Section System Section B N.L. System Section Open Loop^Closed Loop System Section N.L. Delay System Section 1-- Decoupling Figure 2.4: Nonlinear function blocks In an open loop configuration, the nonlinearity does not pose any problem since its input does not depend on its output. Problems arise in a closed loop configuration, where the sequential solution forces the introduction of a time-delay to decouple the input from the output of the nonlinear block (see Figure 2.4). The introduction of this time delay occurs whenever a nonlinear block exists inside a closed loop configuration. Many of the Fortran Functions and Devices blocks are nonlinear and are likely to cause the introduction of such time delays, which the program user may not even be aware of. 2.3.3 Transfer Function with Limiters The general transfer function block has an additional feature that enables the user to use two types of limiters. The function of these limiters is to limit the output of the block if a certain condition is reached. The output of the static limiter shown in Figure 2.5 is Chapter 2. Literature Review^ 10 x G(s) Xrain Dynamic Limiter Static Limiter Figure 2.5: Static and dynamic limiters x = Ku xmin Xmas if Xmin < Ku < Xmas, if Ku < xmin, if Ku> Xmas. Each of the three equations above is a linear algebraic equation of the form of Eq. 2.1, with c = d = 1, hist = 0 inside the limits, and c = 1, d = 0, hist = (X mas or x rnin ) at the limits. Contrary to static limiters, the limiting action of dynamic limiters changes the dynamic behaviour of the transfer function block. These limiters are applied only for first-order transfer functions [7]. The equations for the dynamic limiter are (see Figure 2.5) I x +T"t = Ku if X sa m < X < X mas , X = X m i a^if X < Xmin and (Ku — x) < 0, x = X mas^if x> X mas and (Ku — x) > 0. The dynamic limiter equations can also be put in the format of Eq. 2.1. Inside the limits, the coefficients are: d = K; 2T c^(1+ —A7 ); hist = Ku(t —^(1 2t — At — )x(t — At). 11 Chapter 2. Literature Review^ EMTP TACS t-2At^t-/.t t-20t t-:At Figure 2.6: Time delay between TACS and EMTP Outside the limits the coefficients are the same for both limiters. Although both types of limiters fit into Eq. 2.1, when they occur in closed loops their input will depend on their output. This situation is similar to the case of nonlinear blocks where the introduction of time delays is unavoidable. This is discussed in Section 2.4. 2.3.4 TACS-EMTP Interface As mentioned before, TACS solves the set of control system equations independently from the EMTP and has to communicate with it to exchange information on how the control action is supposed to influence the power system. Figure 1.1 in Chapter 1 outlines the interface process between both programs. It is clear from this figure that for TACS to issue any control signal to the power system it must have received information from the EMTP previously. Figure 2.6 shows how the process takes place as a function of time. TACS solves the control system equations from (t — 20t) to (t — At). Then it gives the values of the variables to the EMTP. The EMTP solves the power system equations for (t — At) to (t). At t, the EMTP gives to TACS the values of the variables Chapter 2. Literature Review ^ 12 in order for TACS to advance the solution from (t - At) to (t). Therefore, the EMTP always uses information from TACS delayed by one time step. 2.4 Time Delays in the TACS-EMTP Approach To summarize what has been discussed so far, there are two types of time delays in the TACS-EMTP approach: • Time delay between TACS and the EMTP - external time delays; • Time delays generated by TACS itself due to nonlinear blocks or transfer function blocks with limiters - internal time delays. In 1983/84 Ma Ren-ming made some revisions to the original TACS code with respect to the order in which the blocks of the control system are solved [7, 8]. One of the objectives of this revision was to decrease the possibility of the introduction of internal time delays in a variety of situations. To illustrate the process, consider Figure 2.7. It consists of 10 transfer function blocks with three of them having limits. The first four blocks form an independent set of equations, as do the last four blocks. To simultanously solve the system comprising the first four blocks, the output of F2 must be calculated first. This means that its equation (Eq. 2.1) must be the last equation of equations 2.3 (the first one to be solved during the backsubstitution - see Section 2.3.1). The limits constraints are applied as soon as the output is available. With the output of F2 known, the output of F3 can be calculated and its equation must be solved next in the backsubstitution process. Then F4 and Fl are solved in this order. F5 and F6 are solved next and in the last set of four blocks the order is: F10 -+F7 -*F8 Chapter 2. Literature Review^ HI Fl F2 13 HI F5 F6 HI LOW F4 LOW F3 LOW Figure 2.7: Possible simultaneous solution with TACS-EMTP The ordering used in this case avoids the introduction of time delays in the solution procedure, but the ordering approach is far from general. Any ordering scheme fails if there are two limiters in one single closed loop. In addition, an ordering algorithm is needed to detect the loops with limiters or nonlinear blocks. Although TACS-EMTP has a built-in ordering routine, a recent improvement [9] has thoroughly dealt with this subject. In terms of what types of problems time delays can cause in the simulation of control and power systems, Lima [10] showed in 1985 that a numerical instability can arise in the simulation of an arc of a circuit breaker, due to the time delay between TACS and EMTP. A more recent work [11] describes numerical instability in simulations of a lead acid battery (DC circuit) due to the same time delay. Examples of numerical oscillations and inaccuracies due to internal time delays are described in [1]. In 1986/87 a study of several problems concerning TACS, after about ten years of existence, was carried out at The University of Wisconsin [12]. This work classifies, among others, the problem of time delays within TACS and between TACS and EMTP as a major TACS problem. It also contains examples of numerically bounded Chapter 2. Literature Review^ 14 oscillations caused by this problem. 2.5 Final Remarks The major work related to the elimination of time delays is reference [8]. This work reduces many of the internal time delays, but not completely. The problem of the time delay between TACS and EMTP has not been considered either by that work or by any other work. The importance of the present thesis lies in the fact that it addresses the problem of eliminating both types of time delays. How important this is will become clear in the next chapter, where the effects of the time delays in transient simulations will be considered. Several cases of inaccuracy, instabilities and bounded oscillations will be analyzed. Chapter 3 Time Delays and the Stability and Accuracy of the Simulation "...if the feedback is delayed too long, oscillations will happen. A simple example that illustrates feedback instability is the common home shower. Typically the shower begins with the water too cold, and the user turns up the hot water to get the temperature he wants. If the adjustment is too strong (he turns the knob too far), he will soon find that the shower is too hot, whereupon he rapidly turns back to cold and soon finds it is too cold. If the reactions are too strong, or alternately the total system (pipes, valve, and human) is too slow, there will result a "hunting" that grows more and more violent as time goes on." R.W. Hamming° ° Numerical Methods for Scientists and Engineers — Dover Publications — New York — 1973. 3.1 Introduction The objective of this chapter is to illustrate, with examples, the numerical instability and stable inaccurate results that arise from the introduction of time delays in transient simulations. The examples used are purposely simple to show the concepts involved. However, some of them have practical importance in specific applications as well. 15 Chapter 3. Time Delays and the Stability and Accuracy of the Simulation^16 R= I t= 0 TACS b EMTP Figure 3.1: Numerical instability in a DC circuit 3.2 Numerical Instability 3.2.1 DC Circuit Figure 3.1 shows a circuit to be simulated by TACS-EMTP [11]. The switch and resistor R2 are simulated inside the EMTP and the voltage source V and resistor R 1 are simulated inside TACS. TACS gives to the EMTP the voltage Va b and the EMTP gives to TACS the current I, as indicated in the figure. TACS is, therefore, working as a current controlled voltage source. At any time step n, I(n) = Vab(n) R2 where Va b(n) is the latest value TACS gives to the EMTP. As TACS lags one time step behind the EMTP, Kb(n) = 1 — R i I(n — 1). Figure 3.2 shows a plot of I versus t for the first 10 time steps, using the two previous formulas for different values of the ratio K = R I /R2 (0.5, 1 and 1.5). From the values of I(n) and Va t, given in the formulas above one can get the Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^17 40 K=1.5 e z . 0.75 2° C' 112 2, 0.5 -20 0.25 0 0 0 2^4^6 -40 ^ 0^2^4^6^8^10 10 Time (s) Time (s) Figure 3.2: Current and Voltage in the circuit of Figure 3.1 following difference equation: I(n) = 1^R 1 T(^ ^1 ). ) — R2 °-1 ' Applying the z-transform technique [13] = 1^R1 — Z-1/(2). R2 R2 — Solving for /(z) gives /(z) ^/R2 z z Ri/R2 The pole of /(z) is z = —R 1 /R 2 . The stability of the solutions /(z) is guaranteed if the poles zi satisfy jzil < 1. Therefore for /(z) to be stable, R 1 < R2. When R l = R2) 121 = 1 and there is a bounded oscillation around the correct value. For R2 = 2R1 (K = 1.5) the solution is clearly unstable (see Figure 3.2). Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^18 I(t) = (I/R,2 )V(t-At) EMTP I TACS V 1 / R2 (a) (b) Figure 3.3: Arc modelling with EMTP-TACS 3.2.2 Arc modelling Consider the circuit of Figure 3.3a where the resistor R2 is simulated by a current source inside TACS [10]. The simulation arrangement is shown in Figure 3.3b. This circuit can be used to model an arc in a circuit breaker by allowing R2 to vary according to the arc characteristic. The EMTP-TACS program discretizes the differential equation of the system using the trapezoidal rule of integration. If the system Laplace transform algebraic equation is available, it is possible to directly obtain the z-transform difference equation by replacing s in the Laplace domain by (2/At) (z=-1 z ) [14]. The circuit of Figure 3.3b is, therefore, transformed to the circuit of Figure 3.4. Carrying out the stability analysis for the current I one obtains I(z) = 1 v1 /Zi j_ (ZR2)(Z1 1 +Z2) Z1 Z2 The poles of I(z) are the roots of 1+ (zR 2 )(Z1 + Z2 ) Zi Z2 Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^19 Ze--(2/At)[(z 1)/(z+1)]L - Figure 3.4: z-transform domain analysis 100 ".7e 50 -50 -100 - -150 ^ 0 1^2^3^4^5^6^7 Time (Rs) Figure 3.5: Numerical instability in the arc modelling With At.— 0 these roots are: z = 1 and z = Ri/R2. For R 1 > R2 the system is unstable even with At^0. Figure 3.5 shows one unstable simulation for L 0.3H, R 1 = 10051, R2 = 10Q and At = 10ms. Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^20 20 40 60 80 100 Time (s) Figure 3.6: Block diagram and exact solution of (t) x(t) = 0 [1] 3.2.3 Sinusoidal Oscillator The examples of instabilities that have been shown so far are caused by the time delay between the power system solution and the control system solution. The example of this section illustrates the effect of time delays inside the control system solution procedure [1]. Consider the block diagram of Figure 3.6 that represents the differential equation (t) x(t) = 0^x(0) = 0, i(0) = 1. The exact solution of this equation is also shown in Figure 3.6. The introduction of one time delay in the feedback path of the block diagram causes the system to become unstable (see Figure 3.7). Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^21 20 40 60 80 100 Time (s) Figure 3.7: Instability due to the internal time delay 3.3 Stable Inaccurate Results 3.3.1 Internal Time Delays As mentioned before, according to the EMTP-TACS approach, only pure transfer function blocks will be solved simultaneously (see Chapter 2). The introduction of any other type of control system block in closed loops will introduce time delays in the solution procedure. Figure 3.8 shows an example of the effect of the introduction of a Fortran function block in a closed loop control system [1]. The function used in the example is the x 1 function which would not have had any influence on the system output if the simulation had correctly been performed. The output of the system with the TACS-introduced time delay jumps to 1.0 instantaneously l . In this case, not only the 1 The latest information from R.H. Lasseter of July 15, 1993 indicates that the wrong answer is not only caused by the time delay, but by a coding error as well. The time delay by itself will still produce incorrect answers, in the form of damped oscillations around the correct solution. Chapter 3. Time Delays and the Stability and Accuracy of the Simulation^22 OUT 0 — TACS solution Correct solution 20^40^60^80 ^ 100 Time (ms) Figure 3.8: Inaccuracy due to the internal time delay [1] time delay, but also the ordering scheme [1], makes TACS work inaccurately. Note that the system is stable throughout the simulation. 3.3.2 Converter Simulation — External Time Delays A recent paper [2] shows another case of the effect of time delays in the simulation of control and power system equations. A simulation of a three-phase thyristor inverter is described. The comparison between the EMTP-TACS and the simultaneous solution is shown Figure 3.9. The three-phase inverter described in the paper is the one simulated in [15] using the EMTP-TACS approach and is also shown in Figure 3.9. To show how sensitive Chapter 3. Time Delays and the Stability and Accuracy of the Simulation ^23 an inverter is to the variation of its firing angle, consider the following values: VLL^480V at 60Hz; L l = 0.2mH; L2 = 1.0mH; Ld = 16mH; E = 600V; Firing angle a = 150°. The simulation is performed with a discretizing time step At = 50ps. A time delay of one At means a time delay of 1.08° electric degrees at 60 Hz. The inverter current Id as a function of a is given by [16] Id^= 1.35VLL cos a + E ^ 3wL3/71- where L, = L1 + L2. With the values shown above one has: 648.23 cos(150)+600 Id = 0.432 = 89.39A for a = 150°, = 0.432^76.49A for a = 151°. 648.23 cos(151)+600 The difference between the two values above is approximately 14.4%. From Figure 3.9 the value of Id which TACS calculates at t = 80ms is 50A. Chapter 3. Time Delays and the Stability and Accuracy of the Simulation^24 Ld 80 "? 60 t 2 40 ii 0 20 90 ^ 100 Time (ms) Figure 3.9: Three-phase thyristor inverter current - TACS inaccurate result [2] Chapter 4 Simultaneous Solution of Power and Control System Equations "There are no good, general methods for solving systems of more than one nonlinear equation." W.H. Press, et al." "...so Newton's method finds its chief use after we have got ourselves close to a root." Forman S. Actonb a Numerical Recipes — The Art of Scientific Computing, Second Edition — Cambridge University Press — 1992. b Nurnerical Methods That Work — Mathematical Association of America — Washington, D.C. — 1990. 4.1 Introduction The general concept of the simultaneous solution of the power system and the control system sets of equations is developed in this chapter. It involves two distinct procedures, each one of them dealing with the elimination of one type of time delay. The implementation of the concepts is carried out in one of the many versions of the EMTP. The main concerns during the implementation were: to keep the EMTP code modifications to a minimum, and to write a code capable of simulating significant case studies while keeping the code to a reasonable size. The code was written as a development code for testing ideas. At this stage, no attempt towards optimizing parts of the code was done. 25 Chapter 4. Simultaneous Solution of Power and Control System Equations ^26 4.2 Basic Considerations 4.2.1 Elimination of the External Time Delay The basic idea of the simultaneous solution is to assemble both sets of equations (from the power and from the control systems) in one single matrix equation. By doing this, the time delay between the solution of both sets of equation is eliminated. The assembling procedure would produce an augmented matrix which is not symmetric in general (see Chapter 2, Eq. 2.2). This fact must be accounted for in a full implementation of this method in the EMTP code, which, currently, handles only symmetric matrices. The present EMTP code works with symmetric matrices only, while TACS must work with non-symmetric ones. In the proposed scheme, the augmented matrix consists of the elements of both matrices put together, plus the connection elements between the power and control systems (these elements are hidden by the nonsimultaneous approach currently in use). The number of these elements depends on the number of elements in the power system controlled by the control system. This does not include switch opening or closing commands. A reasonable assumption is that this number is small. Therefore, the nonsymmetric augmented matrix would require more memory for storage (compared to the sum of memory required by the two separate matrices, one symmetric and the other nonsymmetric), as well as a small amount of extra computational time to handle the connection elements. 4.2.2 Elimination of the Internal Time Delays The existence of internal time delays is due to nonlinearities in closed loops inside the control system. Therefore, to eliminate them, iterations must be allowed. Chapter 4. Simultaneous Solution of Power and Control System Equations ^27 Y ^Y Linear Equations Y ^Y Nonlinear Equations 11 21 12 22 Figure 4.1: Partition of the augmented matrix Although the two opening statements of this chapter sound discouraging, the chosen method for solving nonlinear equations was the Newton-Raphson procedure. This method has fast convergence close to the root. The efficiency of the entire implementation will depend on this fast convergence characteristic. As for the closeness to the root, except at moments of discontinuity of the solution (switching operations), the previous value of a variable will constitute a very good guess for its present root (value). The application of the Newton-Raphson method to the whole augmented matrix would be very inefficient. Therefore, the procedure used is similar to the one described in [17]. The iterative method is confined only to those elements which contain the nonlinearities . Figure 4.1 shows the partition of the augmented matrix, which is necessary to isolate the nonlinear elements. Note that the nonlinear elements are not only those from the control system, but also those from the power system. By the same token, the linear part of the matrix contains linear elements from both systems. The reduction of the matrix to the nonlinear element matrix is illustrated in what Chapter 4. Simultaneous Solution of Power and Control System Equations ^28 follows. The complete matrix equation is [[ 1;12 11 1 m[2Y(11 ,] t )1 ) [(e 12 (( tt )) 1 ) — ( [[ kik2 1)^(4.4) where the vector [ei(t)] not only contains the node voltages in the power system, but also includes the control system block outputs (see Chapter 2, Eq. 2.2). The vector [ki] includes the history terms from both systems; and x1 = ( x2 x n, is the vector of nonlinear variables which [Y22(Y, 1)] depends on. The submatrix [Y22(Y,t)} contains the information about the nonlinear elements and the submatrix [Y21 ] contains the information on how those elements are interconnected to the other elements. Therefore, Eq. 4.4 can be reduced to the subset 2: ({1'22(Y, t )} - [Y21][Yii] -1 [ Y12])[e2(t) ] = [ k2] - [Y211[Y111 -1 [k1]. This is the equation that is to be solved by the Newton-Raphson method. 4.3 Actual Implementation 4.3.1 Interface Between the Written Program and the EMTP One of the methods used by the EMTP to handle nonlinearities is the compensation method. This method excludes the nonlinear branch from the network and replaces it by a current source (see Chapter 12 of [7]). The value of the current will depend on the Thevenin equivalent network (seen from the nonlinear element nodes) and on the nonlinear characteristic itself. Therefore, the equation of the network and of the Chapter 4. Simultaneous Solution of Power and Control System Equations ^29 v km Network Equation / Element Equation Linear Network Solution km i km Figure 4.2: Compensation method in the EMTP nonlinear element must be solved simultaneously. Fig 4.2 illustrates the situation. The way of describing the nonlinear characteristic is, normally, point by point in the EMTP input file. The MicroTran version of the EMTP [18] makes this compensation method approach accessible through an interface which does not require any code changes in the EMTP proper. The EMTP calculates at each time step the Thevenin equivalent circuit of the network seen from whatever nodes are specified in the input file. It then calls (at the same time step) a subroutine "CONNEC" to which it passes the values of the Thevenin impedance and voltage (open circuit voltage) through the argument list. This subroutine must return the value of the currents through the nonlinear elements (Figure 4.3). Therefore, CONNEC enables the user to write a piece of code to describe his/her own nonlinear elements and to interface these elements with the EMTP, without having to touch the EMTP code itself. The important feature here is that, although the nonlinear elements are not included directly in the EMTP matrix, Chapter 4. Simultaneous Solution of Power and Control System Equations^30 EMTP Z'Thev ' open (Linear Network) >^ CONNEC Current Figure 4.3: Subroutine CONNEC and its interface with the EMTP they are solved simultaneously with the network equations. In the simultaneous solution method implementation, full advantage of the CONNEC facility was taken. All code additions (except for the part of program that reads the control system input file and parts of the interpreter - discussed later on) were implemented inside CONNEC. Any controlled power system element must be considered as a nonlinear element to be included inside CONNEC together with the control system equations in this approach. Two consequences of working inside CONNEC must be pointed out: • all the nonlinear elements to be modelled, regardless where they come from (control or power system), must be included in CONNEC, since the network must be linear for the compensation method to be applied; and • there are actually two matrices being built (one inside the EMTP and the other one inside CONNEC). Chapter 4. Simultaneous Solution of Power and Control System Equations ^31 Z They .-..././..",- Control System T Vopen Vcontrol 'input ^ ^ •^ (a) Vinput Control System 'control ^• (b) Figure 4.4: Control actions upon the power system As mentioned before (Chapter 1), there are three ways through which the control system may act upon the network: switch commands (put aside for a while); controlled voltage source; and controlled current source. The two controlled source situations are depicted in Figure 4.4. Figure 4.4a shows the case of a current controlled voltage source. In this situation, the equation to be included in the control system equations is Vcontrol — ZThev X 'input = Vopen, Chapter 4. Simultaneous Solution of Power and Control System Equations ^32 Figure 4.5: Matrix formation inside CONNEC while in the case of Figure 4.4b (voltage controlled current source), the equation is Vinput + ZThev X 'control = Vopen• This is done for each branch which is acted upon by the control system. A typical matrix formation inside CONNEC is shown in Figure 4.5. The switch commands are handled separately. CONNEC controls the firing time of thyristors and requires the EMTP to actually fire them. The firing action appears as changes in the Thevenin equivalent circuit of the network. Chapter 4. Simultaneous Solution of Power and Control System Equations ^33 4.3.2 Capabilities of the Program In what follows, the program capabilities are outlined. The elements which the program can model are similar to those modelled by TACS. The input format for elements similar to TACS elements use similar formats to facilitate the translation of TACS input files. The present implementation can handle most (but not all) of the TACS blocks and is sufficient to allow fairly complex cases to be simulated. Linear Elements Linear elements are transfer functions blocks, which model any n th -order transfer function of the form \ Output^No +^+ + Nm S m G°) — Input K Do + D 1 s + + Dnsn Figure 2.2 of Chapter 2 shows this type of block. For transient solutions, the Laplace operator s can be replaced by the differential operator Tic transforming the equation above into the following differential equation: dx^dnxit, , du^, dmu ) ^= ix wf o^4- • • • + v m D o x^+ • • • -4- .u n dt^dtn^dt^dtm This equation can now be decomposed into a system of n first-order differential equations, as follows: x1 = u1 = dx^dx1 dx,,--1 x2 =^• • • x n = ^ dt'^dt '^dt du^dui^dum_i U2 -=^Um = ^ dt^dt^dt The differential equation becomes an algebraic equation Do x + D i +.. . + D rix , = K (Nou^+... + Armum). Chapter 4. Simultaneous Solution of Power and Control System Equations ^34 The trapezoidal rule of integration can now be used to numerically integrate the first - order differential equations . The internal variables x; can then be expressed as 2 xt = t xi_ 1 — {x1(t — Ot) + Qt xi_1(t —fit)} u7 = Qt 2lj_1 — {u,(t — Ot) + Qt u^_1(t — Ot)}, for i=1. .. nand j = 1.. .m. Expressing x n in terms of x n _ 1 and then x n _ 2 in terms of xn_1 and so on, only x is left at the end of this process. If the same procedure is used for u, the following output -input relationship is found: cx(t) = K d • u(t) + hist(t — Ot). The history terms are recursively calculated using the following formulation: hist 1 = Kd i u(t) — c i x(t) — hist l (t — Ot) + hist 2 (t — zt) histt = Kd u(t) — c;x(t) — hist1(t — Ot) + histi +l (t — Lit) ; hist n = Kdn u(t) — c n x(t) where hist = hist 1 . The coefficients c and d are calculated once using the recursive formula: 2 ),+1d t+1 + ... + ) n(2 i Ot^i^Ot^i Ot )ndn}, c, = ca_ 1 + (-2)'{ z (2 )'D1 + z + 1 ( ( ^ Chapter 4. Simultaneous Solution of Power and Control System Equations ^35 where (a.) is the binomial coefficient, and n^2 • co = E()'Di i.0 At where c = co . The formulas for d are found replacing D by N in the above formulas. The current implementation accepts up to 5 (five) different signals as inputs of this block. Therefore, the input in the equation above can be Input = +u1 u 2 u 3^u5. The transfer function block also accepts the two types of limiters described in Chapter 2 (Figure 2.5) - static and dynamic limiters. FORTRAN Blocks The FORTRAN blocks accept fairly general FORTRAN expressions as the relationship between their output and input. One can have, for example, a Xp = 3.6 • e -1-35.4 - 4.822 • e -217 • xP + 1.2217. The program accepts the following operators and functions: algebraic operators^: +, - , *, /, ** relational operators^: .EQ., .NE., .LT., .LE., .GE., .GT. logical operators^: .0R., .AND., .NOT., .EQV., .NEQV. FORTRANfunctions^:SIN, COS, TAN, COTAN, SINH, COSH, TANH, ASIN, ACOS, ATAN, EXP, LOG, LOG10, SQRT, ABS Appendix A gives a complete description of the FORTRAN block capabilities. Chapter 4. Simultaneous Solution of Power and Control System Equations ^36 Inputs (up to 5) Cntrl> 0 Cntrl ^ Named_Reset Cntrl< 0 Do • OUT(t) = (Do + L • K • INPUT(t)-1- hist(t — At) where hist(t) = K • INPUT(t)— (Do — nt • D 1 ) • OUT(t) OUT(t) = named_reset(t) Figure 4.6: Controlled integrator block Special Elements There are four blocks, called special element blocks, that play a special role described bellow. Time delay^This block delays the input by - T S. RMS block^This block calculates the RMS value of the input in the previous 360° electrical degrees. Controlled integrators^This block is a normal integrator that accepts two controlling signals. Figure 4.6 illustrate how the element works. Point to point functions The general output-input relationship of this block - - is given point by point (two-column format) in the input file. Sources Several types of sources have been implemented in the program. They are listed below (cf. [6]): Chapter 4. Simultaneous Solution of Power and Control System Equations ^37 • Sources with User-defined parameters Source Type Sinusoidal Pulse Ramp Level Signal User-defined Parameters amplitude and frequency start time, amplitude, period and width start time, period and amplitude amplitude • Built-in sources Internal Name TIMEX ISTEP DELTAT ZERO MINUS1 PLUS1 INFNTY PI Output simulation time in seconds number of the current simulation time step size of the time step (TIMEX=ISTEP x DELTAT) 0.0 —1. 1. -boo (very large number) 7r The Interpreter The basic function of the interpreter is to compile a FORTRAN-like expression written in an input file - describing the input-output relationship of a FORTRAN block - in a way that the program can calculate the value of the expression given the values of its variables. This is done in a sequence of subroutines shown in Appendix A. The basic principles behind the construction of such an interpreter are also discussed in that appendix. The process of taking the expression from the input file and transforming it into a vector of symbols to be evaluated from the values of the constituent variables, is divided into two parts. The first part consists of taking the expression and building the vector of symbols. Each FORTRAN block gives rise to a specific vector of symbols. Chapter 4. Simultaneous Solution of Power and Control System Equations ^38 This is done outside the time loop. The second part of the process occurs inside the time loop and is the evaluation of the expressions according to the values which their variables take. Initialization The network steady-state initialization is performed by the EMTP. An ideal situation would be a simultaneous steady-state solution of the network and the control system. Although a feasible alternative [19], this has not been implemented in the EMTP-TACS program. For the purpose of this thesis, only a DC initialization was implemented. This initialization is simpler and yet does not limit too much the range of the program applications. The DC initialization process assumes that the input and output signals of the control system blocks are DC quantities in steady state. In this case, the output-input relationship of a linear (transfer function) block is Xpc = , No ---„ o UDC• This equation fits into the same format of Eq. 2.1, Chapter 2 with c = 1, d = No /Do , and hist = 0 and a system of equations is formed and solved for the unknowns xpc. Integrators pose a problem for this procedure because D o = 0 for this type of linear block. Therefore, their steady-state output must be supplied by the user. Since the initialization procedure implemented at this time is not iterative, it cannot calculate the output of nonlinear blocks. The user must supply those values to the program for the time being. Chapter 4. Simultaneous Solution of Power and Control System Equations ^39 The Newton Method • The "Dishonest" Newton Method There is a provision in the program to speed up the Newton-Raphson method with the so-called "dishonest" Newton method. Assume one has a vector of nonlinear functions Ft involving variables x i (i = I, 2, ... , n) expressed as fol, lows: x2,..., x 7i ) = 0^i = 1,2,...,n. If x denotes the variable vector and F denotes the function vector, at iteration h Jx(h) • (5x = — F x (h) where aFi ax e . The solution vector is then modified at each iteration x(h+1)(h)^r = X + OX. The "dishonest" Newton method calculates the Jacobian matrix J only once (instead of at each time step) and this is done outside the iterative loop. Therefore, the steps, after calculating the Jacobian matrix, are: J^= —F X( h ) and x (h+1) = x(h) Sx. The application of this method is controlled by the user, who can disable the normal Newton method from the input file. Chapter 4. Simultaneous Solution of Power and Control System Equations ^40 • The calculation of the Jacobian matrix The Jacobian matrix is numerically calculated by the program from the values of the functions. The term Ji; of the matrix is calculated as follows: aF, Fi (x i ,...,x i +Ax i ,...,x n )—Fi (x i ,...,x i —Ax i ,...,x,,) ax^ 2 • Ax e If the initial increment Ax i is not enough to provide small increments in the function Fi, it is automatically divided by half and the Jacobian matrix is calculated again. Memory Limitation of the Current Implementation The compilation of the current PC-implementation is performed by a compiler that addresses only 640 kB of memory. Although this poses a restriction on the number of elements in the EMTP and in the number of blocks in the control system simulation program, the implementation is capable of simulating cases with a considerable degree of complexity. Chapter 5 Comparison Between Simultaneous and Non-simultaneous Approach 5.1 Introduction This chapter presents several cases simulated with the simultaneous solution program (SS) described in the previous chapter. It also compares results with TACS simulations, where available. The case studies give an overall view of the behaviour of the implementation of SS, exploring key issues such as accuracy, efficiency and stability. To focus the attention on these key issues, the cases presented here were chosen to be simple, without any additional unnecessary complexity. They can be considered, however, prototype situations that can occur in more complex cases where the key issues are difficult to identify. The first group of cases consists of simulations where only the control system is involved. Static and dynamic limiters, as well as nonlinear functions in closed loops, are simulated. The second group contains cases where the control system and the network work together. The cases illustrate the strength of the simultaneous solution approach. 41 Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 42 For the reader's reference, Appendix B shows the input files for SS and the EMTP (when available) of all cases discussed in this chapter. 5.2 Control System Case Studies This set of cases shows the behaviour of SS in simulations involving control systems alone. The EMTP part of it is disabled and only the program described in the previous chapter is active. Although the cases are simple, they are useful to show how the program performs in general. Note that the case of Section 5.2.3 shows an important advantage of SS over TACS. 5.2.1 Second-Order Linear System Figure 5.1a shows the diagram of the following second-order transfer function: w2 \ Output H(s) = Input —^ s 2 gw,i s + co n2 This system has an analytical solution for a step input [20]: Output(t) = 1 — 1 — • e - wn t • sin(wOt + 0) where /13 = — 4. 2 and 0 = tan' Og. Figure 5.1b shows the SS simulation for the case where w = 4.0 rad/s, = 0.125 and unit step input. The analytical result is superimposed to, and is indistinguishable from the SS result. 5.2.2 Static and Dynamic Limiters in Open Loops Two simple cases are considered in this section. They are taken from Chapter 5 of reference [1]. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 43 S Output ^> Input (a) 1.6 \f\i \ 1.2 0.8 -------- 0.4 o o 2^4^6^8^10 Time (s) (b) Figure 5.1: Second-order linear system: (a) Diagram; (b) Solution for a unit step input. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 44 Figure 5.2a shows the diagram of a static limiter system. Figure 5.2b shows its solutions with and without the static limiter. Figure 5.3 illustrates a similar situation for a dynamic limiter. 5.2.3 FORTRAN Function in a Closed Loop As mentioned earlier, a nonlinear function in a closed loop can lead TACS to produce wrong results. Figure 3.8 of Chapter 3 is reproduced in Figure 5.4, with the SS solution added. 5.2.4 Special Element Block To illustrate the application of special element blocks, Figure 5.5 shows a system uses many of them. Some of these blocks will again be used in more realistic cases in subsequent sections. 5.3 Simulations of Control and Power Systems In this section, both the SS and the EMTP are active. The cases show the behaviour of both programs working together, and demonstrate the advantages of the simultaneous solution as the cases are presented. For these cases, Appendix B shows not only the SS input files, but also the related EMTP files. 5.3.1 Arc Modelling Section 3.2.2 of Chapter 3 mentions a case of severe instability due to EMTP-TACS interaction. This case is, therefore, a good test for the new program. Figure 5.6 shows the superposition of the analytical and the SS simulation results for the system of Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 45 u ^ 2.0 X 1000 s -2.0 Static Limiter u(t), x(t) 6^ 4 (a) x(t) - without limiter x(t) - with limiter 2 0 -2 -4 -6 0 0.04^0.08 0.12 Time (s) (b) Figure 5.2: Static limiter: (a) Diagram; (b) Simulation results. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 46 u> 1000 S / -0.5 Dynamic Limiter u(t), x(t) 6 (a) 4 2 0 -2 -4 -6 0 0.04^0.08 0.12 Time (s) (b) Figure 5.3: Dynamic limiter: (a) Diagram; (b) Simulation results. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 47 4 3 +5 0 OUT 2 — TACS solution 1SS and Correct solutions 0^ 0^20^40^60^80^100 Time (ms) Figure 5.4: Nonlinearity in a closed loop: (a) Diagram; (b) Simulation results N. var3 ^var4 ^var2 y—x2 *1 K=1 var5 var4 Time Delay var5^ var6 RMS ---> var6 varl 0.02 ^ time (s) 0.04 Figure 5.5: Special element blocks - Diagram and simulation results (source = cos wt, w = 377 rad/s and time delay = 2 ms). Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 48 0 0^0.1 ^ 0.2 Time (s) Analytical solution: i(t) 1 ^ri Ril1R2^e R, t, - I' Figure 5.6: Arc modelling - superposition of the analytical and the SS simulation results. Figure 3.4 with the same parameters used in Figure 3.5. Both results are practically identical. 5.3.2 Single-Phase Variable Load Figure 5.7a shows an R-L load, of which the active power P is variable (timedependent), and the ratio P/Q is kept constant [1}. Such an R-L circuit has a Laplace domain representation as shown in Figure 5.7b if zero initial conditions are assumed. The block diagram representation of the circuit is shown in Figure 5.7c. The assumed variation of the active power P is shown in Figure 5.8 and the system parameters are as follows: R= 3 SI (P/Q = 3 - for L. 2.653 mH) Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 49 (a) >^ 1 (b) 1 s -1- c Source Fortran Block s-Block (c) Figure 5.7: R-L load representation; c= R/L. I(s) > Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 50 15 P (W) 1.5 100 150 200 Time (s) 50 Figure 5.8: Active power variation. E m = 3.1623 V It can be shown [1] that the value of L is dependent on the value of P as follows: 1 2wP ,P Q L =^kQ+ TD ) where Q is the reactive power and P/Q = R/wL. If P/Q is kept constant, L will depend only on P. Therefore, the variation of power shown in Figure 5.8 can be accomplished by varying 1/L from 377 (1.5 W) to 3770 (15 W). The block diagram for this variable load problem is shown in Figure 5.9. The variable load is modelled as a control system inside SS and the sinusoidal source (e(t) = sin wt) is modelled inside the EMTP. PINST is the instantaneous power and PAVER is the average power in the circuit'. 1 The instantaneous power of any electrical component is p(t) = v(t) • i(t), where v(t) is the voltage across the component and i(t) is the current through it. The average power over one period of time P is 1 f P =^t p(t)dt. If the voltage and current waveforms have half wave symmetry, P can be calculated using T = I/2f = 0.00833 s for 60 Hz. The values of P can be found continuously using the rolling Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 51 Signal from EMTP to SS K=1 SCR / IOUT SCR*(CONS +PULSE) Sinusoidal Source CONS s + 1131 - -) * PINST PULSE^IOUT 377 Controlled Current Source ^> 3393 120 \_. Sources Delay 0.00833 PAVER Network v SS Figure 5.9: Block diagram - Single-phase variable load [1]. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 52 -15 40^80^120 0 160 ^ 200 Time (s) 20 15 10 5 0 -:.2"N,_." -5 -10 — -15 -- (a) — ..--- iout 0 40^80^120 160^200 Time (s) (b) Figure 5.10: Simulation Results: (a) TACS solution [1]; (b) SS results. The simulation results obtained with TACS and SS are shown in Figure 5.10 for the current and the average power (PAVER). Note that both programs produce practically identical results. average power P(t) [1, Chapter 2] as follows: 1 P(t) = 7 f^1 f—T 0 p(t)dt — T 0^p(t)dt. ,- ., Note that the second integral has the value of the first integral delayed by T seconds. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 53 5.3.3 Three-Phase Variable Load The next three case studies are related to the 230 kV power system of Figure 5.11 [1]. The Thevenin equivalent circuit of the subsystem inside the dashed line is calculated as seen from bus 7. The reduced system is shown in Figure 5.12 along with the variable load connected to bus 13. The Thevenin equivalent values are (in sequence values): Zero Sequence: Ro = 0.13 SI L o = 27.71 mH Positive Sequence: = 0.06 C2 L i = 39.99 mH. The sequence values of the transmission lines are: Zero Sequence: Ro = 0.03167 SZ/km L o = 3.222 mH/km Co = 0.00787 pF/km Positive Sequence: R l = 0.0243 Il/km = 0.9238 mH/km Cl = 0.0126 pF/km. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 54 r 13.8 kV 4 It T 50 MVA 5 0.9 pf 500/v1i/A 10% 2 . 200 MVA X"=12% 200 MV 10 X=10% mi 400 MVA1 -IN 42 200MVA 200 ..- 10%iC MVA 500 MVA X"=1 0.9pf . 10 n 100MVA 3 15mi I 90 mi 12 230kV 60 mi 1 00MVI 10% ^J 180 mi 200 MVA 50 MVA 50 M A 0.9 pf X"=12% X =10% Figure 5.11: Sample power system [1]. Bus7A They A ^BuslA B111112A Three Phase Bus1B Busl2B They B Thevenin Bus7B Equivalent ^ Bus1C^Bus12C T ev C Sequence Bus7C Values 144.4 km three phase 24.14 km three phase transmission lime^transmission lime Transformer Bus 13 A IL13 Sources Variable Load 1 Figure 5.12: Diagram of the Thevenin reduced system [1]. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 55 The line lengths are: from bus 1 to bus 7 — / 1 _ 7 = 144.4 km from bus 1 to bus 12 — 4-12 = 24.14 km. The transformer between bus 12 and bus 13 is represented by its series inductance: LT = 70.16 mH seen from the 230 kV side. As in the single-phase variable load case, it is assumed that the ratio Active Power^P = 3.042 Reative Power Q wL will remain constant during the load variation.. In this case [1, Chapter 2] 1^P3 — 41.63 • V2: where P3 is the three-phase power in MW and Vpu is the per unit value when K ase is 230000/0 . . The variable load, as in the previous case, is simulated inside either TACS or SS. Step Load Variation The load modelling diagram for this case is shown in Figure 5.13. The load variation is shown in Figure 5.14a. There are four three-phase power levels with the following corresponding values of 1/L: = 0.1795 for 9.5 MW; Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 56 IL13 Network as in Figure 5.12 Network SS Bus13A Bus13B Sources K=1 K=1 0.1795 2.3715 CONS1 Busl3A* :CONS1+CONS24 CONS3+CONS4) OONS2 Bus13C Sources K=1 2.8650 1.2293 CC NS3 C 3NS4 Bus13B* Bus13C* CONS1+CONS2+ CONS3+CONS4) ;CONS1+CONS2+ CONS3+CONS4) INT A 1 s+1146 INT_B 1 s+1146 PINSTA CBA INT_C 1 s+1146 120 Delay 8.3ms PA Figure 5.13: Load modelling - step variation [1]. 1 1 1 = 2.551 for 114 MW; = 5.416 for 190 MW; = 6.645 for 206 MW. The voltage of phase A on bus 1 is shown in Figures 5.14b and 5.14c for TACS and SS simulations. The good agreement between the results is readily noticeable. ^ ^ ^ Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 57 206 190 114 9.5 40 80^120 160^200 Time (ms) (a) 300 >^200 100 a) 0)^0 cu 4-.^-100 -200 -300 0 40^80^120 160 200 Time (ms) (b) 300 5 200 ..^100 a) co^0 al .....^ -100 0 >^-200 -300 0 40^80^120 160 ^ 200 Time (ms) (c) Figure 5.14: (a) Load variation — step variation, (b) TAGS simulation result [1], (c) SS simulation result. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 58 IL13 Network as in Figure 5.12 Network SS Bus 13A Bus 13B Busl3C K=1 K=1 Source CBA K=1 Sourc 0.1795 CONS1 5.2364 CONS2 • Busl3A* Busl3B* (CONS1+CONS2) (CONS1+CONS2) INT—A 1 s+1146 PINSTA Busl3C* CONS1+CONS2) INT_B 1 s+1146 1 s+1146 x 120 8 Delay 8.3ms PA Figure 5.15: Load modelling - pulse variation [1]. Pulse Load Variation The load modelling diagram used to simulate this case is shown in Figure 5.15. Figure 5.16 shows the load variation and the voltage at phase A on bus 1 for both TACS and SS simulations, respectively. Both results are practically identical. Load Variation with Thyristor Controlled Reactor A thyristor controlled reactor (TCR) is the basic building block of static Var compensators (SVC). It is usually connected in parallel with a fixed capacitor. The TCR Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 59 186.2 9.8 0 ^ 40^80^120^160^200 Time (ms) (a) 300 ^ 200 a) 100 -100 — O >^-200 -300 ^ 0 40^80^120^160 Time (ms) (b) 5- 200 300 200 100 rn ca 0 -100 — 0 >^-200 - -300 ^ 0 40^80^120 Time (ms) (c) 160 200 Figure 5.16: (a) Load variation — pulse variation, (b) TACS simulation result [1], (c) SS simulation result. ^ Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 60 consists of an inductor in series with an anti-parallel thyristor (two oppositely poled thyristors), which conduct on half-cycles of the supply frequency as indicated in Figure 5.17a. The RMS value of the fundamental component of the current through inductor is Il^— sin , = ^ ), 7r wL^ where o is the conduction interval. This current has the same value as that of a a — sin °normal inductor with L n „,,„ / = L/K where IC = ^ (cf. Figure 5.17b). ir If all current harmonics are neglected, the inductance of Figure 5.17b is the fun- damental frequency model for the TCR of Figure 5.17a, which leads to the block diagram of Figure 5.17c. The TCR control system is represented, in this model, by the parameter IC. A basic voltage controller for the TCR is shown in Figure 5.18 [1]. A voltage coming from the power system (VRMS pu ) is compared with the unity source (UNITY). The difference is integrated with a gain and multiplied by 1/L. If VRMS p .—UNITY> 0, a positive error signal is provided to the integrator. Epu, which plays the role of IC in this model, and BLgam (K/L in Figure 5.17c) will increase. The TCR and its control system is shown in Figure 5.19. This constitutes a simple TCR model which ignores harmonics. Consider again the power system of Figure 5.12. Capacitors are added to bus 1 (7 pF) to provide the necessary dynamic range of capacitive reactive power [1]. A TCR is also added to bus 1 with L = 0.8 H. The block diagram for the SS simulation is shown in Figure 5.20. In the first simulation, the TCR is disconnected from the network. The variable load is assumed to have the variation of Figure 5.21. The TACS and SS results are shown in Figure 5.22. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 61 ma Current fundamental component: / 1 = lirma VW L/K (a) ( a -sin a ) . 7r I Mirn a — sin a IC = ^ 7r I^i(t)>^I I ^I I <^ v(t) ^> I (b) V(s) I(s) S x/L (c) Figure 5.17: Thyristor controlled reactor: (a) Circuit diagram; (b) Circuit model; (c) TCR s-domain model Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 62 Figure 5.18: Basic TCR voltage controller - K/L modelling [1]. V (from the network) Vrmspu /^ E P UNITY ^ 1 L 0 BLgam I (to the network) Figure 5.19: TCR and its control system - a simple model. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 63 Network SS Bus1A^Bus1B^Bus1C K=1 K=1 EMTP1^EMTP2 ^VRMS^ RMS 1 VRMSpu I EMTP3 300 132.7876 ^PLUS1 Source Epu CD CD Figure 5.20: SS block diagram. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 64 P (MW) 209 9.5 50 100 150 200 Time (s) Figure 5.21: Load variation. With the same load variation, the TCR is connected to the network and the simulations result are shown in Figure 5.23. Note the excellent agreement between TACS and SS simulations. 5.3.4 Simple HVDC System Model Assume that the following modifications are made in the power system of Figure 5.11: • capacitors at bus 1 and 10 are removed; • generation at bus 3 is increased from 200 to 500 MW; • the double circuit lines between bus 2 and 9 are removed; and • load at bus 10 is increased from 300 to 500 MVA. To make the system operational under these circumstances, the AC transmission line between bus 1 and 2 is now replaced by a DC line carrying 400 MW. The reactive power consumption of the two converters is assumed to be compensated locally. Figure 5.24 shows the system under study [1]. Three events will be simulated: • initial DC current ramp-up (from 0% at time zero up to 100% at 50 ms); ^ Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 65 270 Vrms 180 a) cu O 90 0 -90 -180 -270 ^ 0 40^80^120 160 200 Time (ms) (a) f7 111 ^270 ^ j^180^► A^kAA^A 90 ^ a) 0)^o as 4-^-90 -180 -270 0 40^80^120^160^200 Time (ms) (b) Figure 5.22: Load variation without TCR: (a) TACS simulation [1]; (b) SS simulation. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 66 5' G) cc; 0 270 ^ Vrms / 180k90 0 -90 — -180 -270 40^80^120^160^200 0 Time (ms) (a) 270 > .. a) 0) co .... O > f ll I T T T Vrms 180 90 0 -90 — -180 — -270 0^ i I 40^80^120 A A ^ A 160 Time (ms) (b) Figure 5.23: Load variation with TCR: (a) TACS simulation [l]; (b) SS simulation. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 67 1 50 MVA^500MVA 0.9 pf 500MVA 3 200 MVA 200 MVA SOOMVA 120 mi 400 MW^DC line 1 50 00 10 mi^ MVA 0.9pf 500 MVA 6 As 11 100MVA 90 mi 00 MVA 60 mi 9 15 mi 200M VA 0.95 pf 12 11 100MVA 300 MVA 0.9 pf 50 MVA 50 MVA Figure 5.24: Power system with a DC line [1]. • DC line to ground fault at 100 ms, followed by • a complete DC load rejection at 150 ms. To simulate these events in terms of the converter current (Id ) and firing angle (a), the diagram of Figure 5.25 is used. The AC system of Figure 5.24 is reduced to its single-phase equivalent, and the HVDC system is simulated using the current injection method [1]. Figure 5.26 illustrates the simulation diagram. In addition to the modelling for Id and a (Figure 5.25), there are three more subsystems to be simulated. They are (Figure 5.27): the voltage tracking system, the DC system itself, and the AC current injection system. The reader is referred to reference [1] for a complete description of this system. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 68 Magnitude = 1.6 Turns on at t = 50 ms Source Ramp up to 1.6 in 50 ms Source IDC TIMEX 107.79* (X111 ^-286*XP ^)1 X1*(X.GE.0) ^>3.6*6 135"3:4.822*e 17."+1.2217 0.1 ALPHA Source Magnitude = 0.31416 Figure 5.25: Id and a modelling [1]. V( t) IDC IACF SS ALPHA Figure 5.26: EMTP SS simulation diagram [1]. - Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 69 from Network cos cot sin cot K=1 Vm= (CPC1+SPS1) Vm CV CP cos cot + SP sin cot Vm eV — Cl* sin o)t -S1 cos cot Vm ^> cV 8V (a) IDC ^ 1 Ilpk — a=0.8933 2*Xc*IDC coe(DT) cos(AF) ^ 0'43 *Vm we coe(AF) + coe(DT) 2 (b)^ cV^eV +elan:DC term = cV*ces, + eV*sin, IACF IACF = - Ilpk*term (c) Figure 5.27: (a) Voltage tracking system; (b) DC system equations; (c) AC current injection [1]. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 70 ALPHA IACF 40^80^120 160 200 Time (ms) (a) ALPHA 160^200 Time (ms) (b) Figure 5.28: IDC, ALPHA and IACF comparisons; (a) TACS simulation [1]; (b) SS simulation. Four parameters from the TACS solution were compared against the SS solution: currents IDC and IACF; firing angle ALPHA; and voltage at bus 1 (BUS1). Figures 5.28 and 5.29 show the results of both solutions. 5.4 Converter Simulations Converter simulations were briefly discussed in Chapter 3, Section 3.3.2. This section presents additional results concerning the thyristor rectifier and inverter simulations. Comparisons between TACS and SS results are also provided. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 71 300 100 0) co 0 -100 -300 40^80^120^160 200 Time (ms) (a) 300 100 a) 0) co^-100 0 -300 0^40^80^120 160 200 Time (ms) (b) Figure 5.29: Voltage at bus 1 (BUS1) comparison; (a) TACS simulation [1]; (b) SS simulation. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 72 5.4.1 Single Phase Thyristor Rectifier Bridge - A single-phase thyristor rectifier is shown in Figure 5-30a. The main functions of its controller with its block diagram appears in Figure 5-30b. The TACS [16] and SS simulation results for this rectifier are shown in Figure 5-31. Although the curve of id is the same in both simulations, the curve of v d shows a difference. The TACS solution presents pronounced peaks at the times of thyristor firings (see Figure 5.31). It was possible to reproduce these peaks with the SS program by introducing, artificially, one time delay between the program and the EMTP. This shows that the origin of the peaks is the external time delay between TACS and EMTP. 5.4.2 Three Phase Thyristor Inverter - In the case of a three-phase thyristor inverter, the effect of the external time delay between TACS and EMTP appears on the DC side current (id). Figure 3.9 shows the inverter diagram and the TACS and SS simulation results. As stated before, there is a difference of 14% between the two currents. If, again, one artificially introduces a time delay between SS and the EMTP, it is possible to reproduce the TACS results. 5.4.3 Remarks on Converter Simulation with the EMTP Although not apparent from the discussion above, there is an extra non-simultaneous phenomenon in effect in the thyristor converter simulations. This effect is a characteristic of the EMTP and is present both in the SS and TACS results. Thyristors are modelled by the EMTP as controlled switches. Their closing times are controlled by the equations describing the firing angle control. The opening action Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 73 From EMTP Source 1.0 K=1 (a) h To EMTP To EMTP rsA 120 S Source 0.25 • RAMP Comparator / :RAMP .GE. VCONTL, • PLS1-2 Delay Block 1/2 cycle VCONTL VSA RAMP Firing Signals ■11 PLS1-2 PLS3-4 (b) t Figure 5.30: Single-phase thyristor rectifier. PLS34 •^ ^ Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 74 270 ^d 200 130 - a) ^ho 60-^ vd Cd^-10 - 0 430 ^ -150 ^ 35^ 36 37 Time (s) (a) >^ 200 r- 1 30 Pili^ 60 CD b.0^-10 Cd 17( 4dri. 80 " - - — -150 ^ 35 36 37 Time (s) (b) Figure 5.31: Single-phase rectifier simulations — voltage vd: (a) TACS results; (b) SS results. Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 75 i ?witch opens SWITCH current forced to zero in the next step T FIRE 1 closed here Figure 5.32: Closing and opening of EMTP switches. takes place as soon as the current goes through zero, from a positive to a negative value. Since the EMTP uses a fixed time step to integrate the differential equations of the system, the closing and opening actions will generally not occur at the specified instants of time [7]. Instead, there is a delay between the desired time and the real one. Figure 5-32 illustrates the situation for a given firing angle. The closing action takes place at the time step nearest to TFIRE (some EMTP versions may close at the time step following TFIRE [7]). The opening action takes place after the current has gone through zero, which is detected by the current sign change. As presented in [21] this problem can be eliminated by the following procedure: 1. interpolation to obtain the variables at the point of discontinuity (switch opening or closing); and 2. integration of the system equations after the discontinuity over two half time steps (At/2) using the backward Euler method, to properly re-initialize the Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 76 history terms of all the system components [14]. The interpolation is performed to catch the exact moment of the discontinuity. The use of the backward Euler method is necessary to re-initialize the variables after the switch operation. Reference [7] suggests that the integration over just one half time step after the interpolation is enough to re-initialize all the variables. This is true only for isolated lumped elements, when one can assume that either a constant voltage (across capacitances) or constant current (through inductances) will occur during the one half time backward Euler step. One of those two conditions is essential for re-initializations of the variables at the instant of discontinuity. For elements connected either in series or in parallel, two backward Euler integration over two half time steps are needed to provide the right re-initialization for the variables of the system. This method is illustrated in Figure 5.33. The effectiveness of the backtracking procedure can be seen in Figure 5.34, where the inverter of Figure 3.9 is simulated with the normal EMTP (EMTP in the figure) and with a modified version of the EMTP that includes the backtracking method (bktrack in the figure). Because the normal EMTP misfires the thyristors, its solution shows a short-circuit of the inverter that does not exist, as demonstrated by the modified EMTP solution. The time step used in the simulations was 200ps. 5.5 Final Remarks As mentioned in the introduction to this chapter, the case studies presented here were intended to show the performance of a new, more accurate and stable approach of the simultaneous solution, when compared to the current non-simultaneous option. In general, when TACS performs well, so does the proposed method in terms Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 77 i( t) interpolation iLt? ........ ---> t t o +A t t - At 1 : et^et 2^2 1 Vab 1 a ........, 1 \..-----^---.-^ 1 t CLOSE^I 41, t + At a^ ■ CLOSE t - At^t At^At 2^2 backward Euler steps Figure 5.33: Backtracking procedure. b Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 78 1000 a) c3) 500 0 C 0 a) L.. '5 EMTP bktrack LiartA ^ vi -500 vd -1000 0.16 0.185 0.21 0.235 0.26 Time (s) Figure 5.34: Inverter simulation: instability due to thyristor misfiring Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 79 of accuracy. In those cases where there is an instability in the TACS solution, the proposed method outperforms TACS, showing no instability in its solution. In cases where an analytical solution exists, it is possible to observe how accurate the SS approach performs. The efficiency of the SS method is related to the application of the NewtonRaphson method, which is used whenever there is a nonlinearity in a closed loop. Nonlinearities outside closed loops pose no problem since they do not require any iterative solution method. As mentioned in Chapter 4, for slowly converging cases, there is an option of the "dishonest" Newton method, which can speed up the solution somewhat. Although most of the cases simulated by the author, and the great majority of the cases shown in this thesis, do not require more than three iterations, on average, at each time step, some cases may require more than three iterations [2]. This situation is caused by the fact that the convergence of the Newton-Raphson method can be guaranteed only if one is close enough to a root. Since almost any kind of equation can be simulated with SS, one cannot guarantee either fast convergence, nor convergence per si, in the general case. As mentioned before, at each time step, the guess for the value of the root is normally its previous value. This may be sufficient in the majority of situations. In cases of poor convergence, one may try to improve the method of prediction of the present value of the root. The author has performed a theoretical analysis on the stability and accuracy of various prediction method [22] and has implemented the linear prediction method in the current version of the SS program. The application of the linear prediction method is controlled by the user. The case shown in Figure 5.24, for example, required 8 iterations to converge when the normal guess method is used (previous value) and the "dishonest" Newton Chapter 5. Comparison Between Simultaneous and Non-simultaneous Approach 80 method is applied. If the linear prediction and the full Newton method is applied the number of iterations decreases to 4. Timing is not shown in this thesis because the code developed is not a production code. In a production code, simulation time would increase because of iterations, but that is better than having unreliable answers. Chapter 6 Conclusions The main contributions of this thesis are: • it proposes and describes a method for solving the network equations and the control system equations simultaneously; • it describes the implementation of the method within the framework of the EMTP, which is one of the most widely used general purpose programs for the calculation of transient phenomena in electric power system; • it presents a set of case studies showing the main characteristics of the present implementation and how well it performs when compared with other existing implementations. The proposed method of simultaneous solution of the control and power systems equations has proved to be, not only accurate, but also stable. A number of example circuits were simulated to illustrate the inaccuracy and instability problems that can occur when artificial time delays are introduced in the simulation (as in TACS). The situations typified in the case studies are likely to occur in more complex simulations where the effect of these time delays can be much harder to isolate. On the programming side, the main objective was to write enough code to demonstrate the feasibility of the proposal. Code optimization was not considered during 81 Chapter 6. Conclusions^ 82 the development of the ideas. In this sense, more work has to be done if this code is to become production-grade code. It may be appropriate to indicate some of the improvements which the current implementation might need to be transformed into a useful general-purpose code. The first decision is related to the question whether the control system should be simulated inside a subroutine (as it is in the current implementation). A better solution would be to merge the control and network equations into one single system of equations. This would require a major code modification in the EMTP. The new code would have to be able, among other things, to handle non-symmetric matrices. Sparsity techniques would still be applied. One necessary modification in the current SS code, in the case of a major code revision, is the implementation of a procedure similar to the CDA, to damp numerical oscillations originating from within the control system. The current CDA procedure is capable of damping only numerical oscillations caused by events occurring inside the network. Although the implementation has performed well in terms of accuracy and stability, some cases may pose some problems in terms of efficiency, requiring an excessive number of iterations. To better guess the value of the variables and improve the convergence of the Newton method, prediction methods should be further studied. Another possibility for improvements is the interpreter of functions outlined in Appendix A. This interpreter can be extensively improved and, perhaps, also be used as a front-end interpreter which can allow digital control systems to be read in from the EMTP input file. Bibliography [1] R.H. Lasseter. EMTP Workbook IV. Electric Power Research Institute, Palo Alto - California, 1989. [2] A.E.A. Araujo, H.W. Dommel, and J.R. Marti. Simultaneous Solution of Power and Control Systems Equations. Paper 93 WM 241-0-PWRS, IEEE PES Winter Power Meeting, Columbus, Ohio, 1993. [3] H.W. Dommel. Digital Computer Solutions of Electromagnetic Transients in Single- and Multiphase Networks. Trans. Power App. and Syst., PAS-88:388399, 1969. [4] W.F. Long. Power Systems Analysis Using Electromagnetic Transients Programs. IEEE/PES Transmission and Distribution Conference, September 1991. [5] L. Dube and H. W. Dommel. Simulation of Control Systems in an Electromagnetic Transients Program with TA CS. In Proceedings of IEEE PICA Conf., pages 266-271, 1977. [6] Bonneville Power Administration, Portland-Oregon. Electromagnetic Transients Program Rule Book, April 1982. [7] H. W. Dommel. Electromagnetic Transients Program Reference Manual (EMTP Theory Book). Department of Electrical Engineering - The University of British Columbia, Vancouver - Canada, 1986. [8] Ma Ren-ming. The Challenge of Better EMTP-TACS Variable Ordering. EMTP Newsletter, 4(4):1-6, August 1984. [9] D.G. Chapman, A.J. Helse, F.L. Alvarado, and N.J. Balu. Ordering and Initialization Algorithms for Control Systems Models. IEEE Trans. Power Delivery, 3(3):1189-, July 1988. [10] J.A. Lima. Numerical Instability due to EMTP-TACS Inter-relation. EMTP Newsletter, 5(1):21-33, 1985. 83 84 Bibliography^ [11] M. Ceraolo. A Case of Numerical Oscillations caused by the TACS-EMTP Interface. EMTP Newsletter, 5(3):4-10, September 1992. [12] R.H. Lasseter. Tacs improvements - classification of problems - final report. Technical report, Department of Electrical and Computer Engineering - University of Wisconsin, Madison, June 1987. [13] Robert A. Gabel and Richard A. Roberts. Signals and Linear Systems. John Wileys & Sons, Inc., New York - USA, 1987. [14] J.R. Marti and J. Lin. Suppression of Numerical Oscillations in the EMTP. IEEE Trans. on Power Syst., 5(2):394 402, May 1990. - [15] N. Mohan. Computer Exercises for Power Electronics Education. Department of Electrical Engineering-University of Minnesota, Minneapolis, 1990. [16] N. Mohan, T.M. Undeland, and W.P. Robbins. Power Electronics Converters, Applications and Design. John Wiley & Sons, New York - USA, 1989. - [17] H.W. Dommel. Nonlinear and Time Varying Elements in Digital Simulation of Electromagnetic Transients. Trans. Power App. and Syst., PAS-90:2561-2567, November-December 1971. - [18] MicroTran Power System Analysis Corporation, Vancouver-Canada. Transients Analysis Program Reference Manual, 1991. [19] X.J. Cheng and C. Hatziadoniu. A Gradient-Based Algorithm for the State Initialization of Control Systems. IEEE Trans. on Power Syst., 6(4):1349-1355, November 1991. [20] Richard C. Dorf. Modern Control System. Addison-Wesley Publishing Company, Inc, New York - USA, 1989. [21] A.E.A. Araujo, H.W. Dommel, and J.R. Marti. Converter Simulations with the EMTP: Simultaneous Solution and Backtracking Technique. Paper APT 28618 - Accepted for publication, Joint International Power Conference (IEEE NTUA) - Athens Power Tech, September 5-8, 1993, Athens, Greece. [22] A.E.A. Araujo, H.W. Dommel, and J.R. Marti. Numerical Instabilities in Power System Transients Simulations. In Proceedings of the IASTED Conference, Power Systems and Engineering, pages 176-180, 1992. Bibliography^ 85 [23] Jean-Paul Trembley and Paul G. Sorenson. The theory and Practice of Compiler Writing. McGraw-Hill Book Company, New York — USA, 1985. [24] S. Bhattacharya. Personal communication to Dr. H. W. Dommel. 1991. Appendix A Interpreter Basic Principles This appendix is based on references [23] and [24]. Its objective is to broadly describe the principles used to write the interpreter which reads, from the input file, FORTRAN-like mathematical expressions and calculates their values, given the values of their variables, at each time step. The project of the interpreter begins with the definitions about the syntax, the operators, the functions, types of variables and function arguments it will accept. These definitions for the present interpreter are as follows: algebraic operators^: +, - , *, / , ** relational operators^: .EQ . , . NE . , . LT . , . LE . , . GE . , . GT . logical operators^: .0R., .AND., .NOT., .EQV., .NEQV. FORTRAN functions^:SIN, COS, TAN, COTAN, SINH, COSH, TANH, ASIN, ACOS, ATAN, EXP, LOG, LOG10, SQRT, ABS variables^: any 6 character alphanumeric name arguments^: any variable or any numerical arguments, either integer or real parenthesis nesting^: accepted In addition, all the rules for constructing standard FORTRAN expressions apply. 86 A.1 Compilation of Arithmetic Expressions The way an arithmetic expression is written in the input file is the normal way to write them, for example: (a + b) (c d). This type of expression is called infix expression. This format is easier to understand but it is not efficient to evaluate. A complex expression in the infix format normally requires parenthesization. The expression below is known as a fully parenthesized expression ((w ^ (x / y))^(z x 5)). ^ 3^1^2 2 The integers below the algebraic operators specify the parenthetical level of each operator. To evaluate this expression, the subexpression which contains the operator with the highest parenthetical level is calculated first. If two or more operators have the same parenthetical level, they are evaluated from left to right. The process continues with the subexpressions at the next highest levels until all the operations have been performed. The repeated scanning of an infix expression can be avoided if it is first transformed to a parenthesis free format such as (operand) (operand)(operator) ^(suffix) or (operator) (operand) (operand) ^(prefix) instead of the infix form (operand) (operator) (operand) ^(prefix). The prefix and suffix forms of the notation are also called prefix Polish and reverse Polish, respectively. 87 The evaluation of a Polish expression is much easier than the related infix expression. The infix expression (x x ((y + z) — 5)) has the following reverse polish equivalent form xyz + 5 — x. A reverse Polish expression can be evaluated by applying the following rules repeatedly until all operators have been processed: 1. Find the leftmost operator in the expression. 2. Select the two operands immediately preceding the operator found in step 1. 3. Perform the indicated operation. 4. Replace the operator and its operands with the result. The previous reverse Polish expression is evaluated as follows (applying the four rules above): a -------, First evaluation x yz + 5 — x ---3 (y + z). # „....-....., Second evaluation x a5 — x —4 (y + z) — 5. 'Y e....—"—., Final evaluation x/9 x -- 7 = x x ((y + z) — 5)). Comparing this evaluation process with the evaluation of the infix expression shown previously, one can see how much easier it is to calculate Polish expressions. 88 Table A.1: Simple Precedence Function Symbol Precedence — x, / Variables and nonnegative integers f 1 2 3 0 A.2 From Infix to Polish Notation The conversion of an infix expression without parenthesis to a reverse Polish expression requires an operator precedence function and the use of a stack. As seen above, in both forms of the Polish format (suffix and prefix) the variables and constants appear in the same order as in the infix format. The operators, however, are reordered and this is the reason for the use of a stack. The output of the conversion is put in an output string. For the simple example discussed in this section, assume the existence of the four arithmetic operators, variables and nonnegative integers. Table A.1 shows the precedence function for this case. The special symbol # has the smallest precedence value in the table and is used to initialize the stack and to ensure that the stack is always nonempty. The stack is initialized with the symbol #. The algorithm to be described is mainly based on the comparison of the precedent function value of the current input and the precedent function value of the top element of the stack. If the precedent value of the input symbol is greater than that of the symbol on top of the stack, then the 89 Table A.2: Transformation of a— b x c+ d/5 into Reverse Polish Current input Contents of stack ^Reverse Polish symbol^(rightmost symbol is expression top of stack)^(output string) a^#a a b^#—b^a x^#— x^ab c^#— x c^ab #+^abc x — d^#+d^abc x — / #+/^abc x —d 5^#+/5^abc x —d tt^ abc x —d5 + current input symbol is pushed onto the stack and the next input symbol is scanned. If, however, the precedence value of the current input symbol is less than or equal to that of the stack top symbol, this top element is written in the output string. After this, the same current input symbol is compared to the new element on the stack. The process stops when all the input symbols have been analyzed. This algorithm is illustrated in Table A.2, for the infix expression a — b x c+ d/5. To handle more complex parenthesized expressions two precedent functions are needed: the input precedence function f and the stack precedence function g. The right and left parenthesis perform the same function as the special symbol # in the simple algorithm previously described. Table A.3 shows the two functions for the implementation of this thesis. The complete algorithm to translate infix to reverse Polish is: 90 Table A.3: Input and Stack Precedence Functions Input precedence Stack precedence function, f^function, g 21^—1 (^ variables and constants^19^20 intrinsic function (COS, SIN, etc.) 17^18 ** (exponential)^16^15 *, /^ 13^14 +, -^ 11^12 .EQ.,.NE.,.GE.,.GT.,.LT.,.LE. ^9^10 .NOT.^ 7^8 .AND.^ 5^6 3^4 .OR.^ .EQV.,NEQV.^ 1^2 0^not inserted in stack —1^not inserted in stack )^ Symbol^ 1. Initialize TOP NEXT 4- 1 4- STACK(TOP) OUTBUFFER 4- " 2. Scan the infix expression Repeat through step 5 while there are unprocessed tokens (symbols) 3. Get next input token and set right parenthesis flag 91 LAST 4 NEXT - NEXT 4-- GET_NEXT_TOKEN (INFIX) RIGHTJ'AREN 4-- false 4. Put in the stack if token higher precedence than NEXT Repeat while f (NEXT)< g(STACK (TOP) ) TEMP STACK (TOP) - TOP 4 - TOP - 1 if f (NEXT) < g(TEMP) then OUTBUFFER 4 OUTBUFFER//TEMP - end if if TOP=0 then write(Invalid expression') end if 5. Push NEXT to STACK but not `)' or ' , '. The function of the delimiter 'Y is to separate arguments of an intrinsic function. Therefore this delimiter must not close the left parenthesis. Its action is similar to a right parenthesis and a comma should not be inserted in the stack at any time. if not RIGHT_PAREN and NEXT comma then TOP 4 - TOP + 1 STACK (TOP) 4-- NEXT end if 6. Finished? or error 92 if TOP=0 then exit the program el se write('1nvalid expression) end if The FORTRAN code for the precedence function f and g is shown right below. Next it is shown the FORTRAN code for subroutine POLISH, which transforms a infix to a reverse Polish expression. The subroutine EVAL (also shown below) takes the output of POLISH (OUTPUT) and transforms it into a character vector that has, in each position, a different token of the reverse Polish expression. INTEGER FUNCTION PFUNC(X) CHARACTER*6 X C^This Integer Function gives the input precedence function f C^Logical operators IF(X .EQ. '.EQV. ' .OR. X .EQ. '.NEQV.') THEN PFUNC = 1 ELSE IF(X .EQ. '.OR. ^') THEN PFUNC = 3 ELSE IF(X .EQ. '.AND. ') THEN PFUNC = 5 ELSE IF(X .EQ. '.NOT. ') THEN PFUNC = 7 C^Relational operators ELSE IF(X .EQ. '.EQ.^' .OR. X .EQ. '.NE.^' .0R. 1^X .EQ. '.LE.^' .OR. X .EQ. '.LT.^' .OR. 2^X .EQ. '.GE.^' .OR. X .EQ. '.GT.^') THEN PFUNC = 9 C^Arithmetic operators ELSE IF(X .EQ. '+ ^' .OR. X .EQ.^') THEN PFUNC = 11 ELSE IF(X .EQ. '* ^' .OR. X .EQ. '/^') THEN 93 PFUNC = 13 C^Exponential ') THEN ELSE IF(X .EQ.^'** PFUNC = 16 C^Intrinsic functions ' .OR. X ELSE IF(X .EQ. 'COS ' .OR. X 1^X .EQ. 'SIN ' .OR. X 2^X .EQ. 'TAN 3^X .EQ. 'ATAN2 ' .OR. X 4^X .EQ. 'COSH ' .OR. X PFUNC = 17 ELSE IF(X .EQ. 'DIM ' .OR. X .EQ. 'SIGN ' .OR. X 1^X 2^X .EQ. 'SQRT ' .OR. X ' .OR. X 3^X .EQ. 'EXP PFUNC = 17 C^Parenthesis ELSE IF(X .EQ. '(^') THEN PFUNC = 21 ELSE IF(X .EQ. ')^') THEN PFUNC = -1 C^Delimiter ELSE IF(X .EQ.^') THEN PFUNC = 0 C^Variables and Constants ELSE PFUNC = 19 ENDIF RETURN END .EQ. .EQ. .EQ. .EQ. .EQ. 'ACOS^' 'ASIN^' 'ATAN^' 'TANH^' 'SINH^') .OR. .OR. .0R. .OR. THEN .EQ. .EQ. .EQ. .EQ. 'AMOD^' 'ALOG^' 'ALOG10' 'ABS^' .OR. .OR. .OR. ) THEN INTEGER FUNCTION GFUNC(X) CHARACTER*6 X C^This Integer Function gives the stack precedence function g C^Logical operators IF(X .EQ. '.EQV.' .OR. X .EQ. '.NEQV.') THEN GFUNC = 2 ELSE IF(X .EQ. '.OR.') THEN GFUNC = 4 94 ^ ELSE IF(X .EQ. '.AND.') THEN GFUNC = 6 ELSE IF(X .EQ. '.NOT.') THEN GFUNC = 8 C^Relational operators ^ELSE IF(X^.EQ.^'.EQ.' .OR.^X^.EQ.^'.NE.' .OR. 1^X^.EQ.^'.LE.' .OR.^X^.EQ.^'.LT.' .OR. 2^X^.EQ.^'.GE.' .OR.^X^.EQ.^'.GT.') THEN GFUNC = 10 C Arithmetic operators ') THEN '^.OR.^X^.EQ. ELSE IF(X^.EQ.^'+ GFUNC = 12 ') THEN '^.OR.^X^.EQ.^'/ ELSE IF(X^.EQ.^'* GFUNC = 14 Exponential C ') THEN ELSE IF(X .EQ. '** GFUNC = 15 Intrinsic functions C ' .OR. X .EQ. 'ACOS ' .OR. ELSE IF(X^.EQ. 'COS ' .OR. X .EQ. 'ASIN ' .OR. 1^X^.EQ. 'SIN ' .OR. X .EQ. 'ATAN ' .OR. 2^X^.EQ. 'TAN 3^X^.EQ. 'ATAN2 ' .OR. X .EQ. 'TANH ' .OR. 4^X^.EQ. 'COSH ' .OR. X .EQ. 'SINH ') THEN GFUNC = 18 ' .OR. X .EQ. 'AMOD ' .OR. ELSE IF(X .EQ. 'DIM 1^X^.EQ. 'SIGN ' .OR. X .EQ. 'ALOG ' .OR. 2^X^.EQ. 'SQRT ' .OR. X .EQ. 'ALOG10' .OR. ' .OR. X .EQ. 'ABS^' ) THEN 3^X^.EQ. 'EXP GFUNC = 18 C^Parenthesis ELSE IF(X .EQ. '(^') THEN GFUNC = -1 ELSE IF(X .EQ. ')^') THEN GFUNC = -1 C^Delimiter ELSE IF(X .EQ.^') THEN GFUNC = 0 C^Variables and Constants ELSE 95 ^ GFUNC = 20 ENDIF RETURN END SUBROUTINE POLISH(INPUT,OUTPUT) C^Transforms a expression from the infix to reverse Polish form. In C^the transformation, it uses two functions: input precedence C^function PFUNC; and stack precedence function GFUNC. C^INPUT has its tokens separated by colons (delimiter). For example, C^COS;(;VAR1;+;VAR2;) for the infix expression COS(VAR1+VAR2) C^OUTPUT also has its tokens separated by colons. For example, for C^the same infix expression, OUTPUT would be VAR1;VAR2;+;COS. CHARACTER*6 COMMA CHARACTER*400 INPUT, OUTPUT, OUTBUF INTEGER^TOP, PFUNC, GFUNC, IBUF, POINTR INTEGER^SIX, BUFSIZ, MAXSTK LOGICAL^RPAREN PARAMETER^(SIX = 6, BUFSIZ = 400) PARAMETER^(MAXSTK = 100, COMMA = CHARACTER*6 STKTOP, STACK(MAXSTK), NEXT, BLANK6 EXTERNAL^PFUNC, GFUNC DATA^BLANK6/' C^Initialize variables DO 10 I = 1, BUFSIZ OUTBUF(I:I) = BLANK6(1:1) 10^CONTINUE DO 15 I = 1, MAXSTK STACK(I)^= BLANK6 15^CONTINUE IBUF = 1 TOP = 1 POINTR = 1 STACK(TOP) = '( C^Get next token 100 CONTINUE C^Get the next token from the input string and put in the C^variable NEXT. If no more token is available ISTAT is set C^to ZERO otherwise ISTAT is one. 96 NEXT = BLANK6 J = 0 ISTAT = 1 DO 50 I = POINTR, BUFSIZ IF(INPUT(I:I) .EQ.^GOTO 60 IF(INPUT(I:I) .NE. BLANK6(1:1)) THEN J = J + 1 IF(J .GT. SIX) THEN WRITE(*,*) 'TOKEN IS GREATER THAN 6 CHARACTERS' STOP 'INTERNAL ERROR' ENDIF NEXT(J:J) = INPUT(I:I) ENDIF 50^CONTINUE C^Check for J ? IF(J .GT. 0) THEN WRITE(*,*) 'LAST TOKEN DOES NOT HAVE A ENDING DELIMITER' STOP 'ERROR' ENDIF ISTAT = 0 60^CONTINUE C^Reset the pointer POINTR = I + 1 IF(ISTAT .EQ. 0) THEN NEXT = ') ENDIF C^Set right parenthesis flag RPAREN = .FALSE. C^Start STKTOP = STACK(TOP) 70^IF(PFUNC(NEXT) .LE. GFUNC(STKTOP)) THEN TOP = TOP - 1 IF(PFUNC(NEXT) .LT. GFUNC(STKTOP)) THEN CALL TOBUFF(IBUF,OUTBUF,STKTOP) ELSE RPAREN = .TRUE. GOTO 150 ENDIF IF(TOP .EQ. 0) THEN 97 WRITE(*,*) 'INVALID EXPRESSION' STOP 'ERROR' ENDIF STKTOP = STACK(TOP) GOTO 70 ENDIF 150 CONTINUE C^Push NEXT to STACK but not ")" or "," C^The function of the delimiter "," is to separate arguments of an C^intrinsic function. Therefore, this delimiter must not close the C^left parenthesis. Its action is similar to a right parenthesis C^and a comma should not be inserted in the stack at any time. IF((.NOT. RPAREN) .AND. NEXT.NE . COMMA) THEN TOP = TOP + 1 IF(TOP .GT. MAXSTK) THEN WRITE(*,*) 'STACK OVERFLOW' STOP 'ERROR' ENDIF STACK(TOP) = NEXT ENDIF C^Continue till ISTAT = 0 IF(ISTAT .NE. 0) THEN GOTO 100 ENDIF C^Finished? IF(TOP .EQ. 0) THEN OUTPUT = OUTBUF RETURN ELSE WRITE(*,*) 'INVALID EXPRESSION' STOP 'ERROR' ENDIF END SUBROUTINE TOBUFF(IBUF,BUFF,TOKEN) C^Routine TOBUFF writes the token to the output buffer OUTBUF. C^It also trims each token and inserts a token delimiter ';' C^after each token. PARAMETER^(SIX = 6, BUFSIZ = 400) 98 CHARACTER*400 BUFF CHARACTER*6^TOKEN, TRIMTK C^Trims the incoming token J = 0 DO 10 I = 1, SIX IF(TOKEN(I:I) .NE. ") THEN J = J + 1 TRIMTK(J:J) = TOKEN(I:I) ENDIF 10^CONTINUE IF(IBUF .LE. BUFSIZ) THEN IF(J .GT. 0) THEN BUFF(IBUF:80) = TRIMTK(1:J) // ';' IBUF = IBUF + J + 1 ENDIF ENDIF RETURN END SUBROUTINE EVAL(OUTPUT,N,ISTACK) C^Creates a stack ISTACK (a vector) that contains, in each position, C^one token of OUTPUT. ISTACK has a proper form to be evaluated by C^subroutine CALCU. IMPLICIT^REAL*8 (A-H O-Z), INTEGER*4 (I-N) PARAMETER^(NSTKSIZ=400,NVAR=30) CHARACTER*400 OUTPUT CHARACTER*6 ISTACK(NVAR) N = 0 K = 0 DO 10 I = 1, NVAR ISTACK(I) = ' 10^CONTINUE DO 50 I = 1, NVAR J = 0 100^CONTINUE K = K + 1 IF(K .LE. 400) THEN IF(OUTPUT(K:K) .EQ. ';') GOTO 50 IF(OUTPUT(K:K) .NE. ") THEN 2 99 J = J + 1 IF(J .LE. 6) THEN ISTACK(I)(J:J) = OUTPUT(K:K) ELSE GOTO 50 ENDIF ENDIF GOTO 100 ENDIF GOTO 200 50^CONTINUE 200 CONTINUE DO 400 I = 1, NVAR IF(ISTACK(I) .NE.^') N = N + 1 400 CONTINUE RETURN END A.3 Evaluation of a Reverse Polish expression To evaluate the reverse Polish expression the following algorithm is used: 1. Initialize TOP 4- 0 VALUE 4- 0 2. Scan the reverse Polish expression Repeat through step 5 while there are unprocessed tokens (symbols) 3. Get next input token 1 00 NEXT < — GET_NEXT_TOKEN (REVERSE) 4. Check if it is an operand (control system variable or real number) if NEXT = operand then TOP 4-- TOP + 1 STACK (TOP) = Value_of_(NEXT) GOTO step 3 end if 5. Find what type of operator? Apply proper operations 6. Finished? or error if TOP=1 then VALUE = STACK (TOP) exit the program else write(Invalid expression') end if The general algorithm to read the infix expression from the control system input file and transforms it into a reverse Polish expression to be calculated at each time step is described below. 1. Outside the time loop. (a) Read infix expression and take all the blank characters out of it. 101 (b) Separate the tokens by colons (delimiters). (c) Call subroutine POLISH (d) Call subroutine EVAL 2. Inside the time loop — Call, at each time step, the subroutine CALCU (FORTRAN code shown below) to evaluate the expression (following the algorithm above). CALCU has the ability to recognize control system variable from real numbers that could appear in the FORTRAN expressions. C^Interpreter SUBROUTINE CALCU(ISTACK,N,VALUE,II) IMPLICIT^REAL*8 (A-H 2 O-Z), INTEGER*4 (I-N) PARAMETER^(NSTKSIZ=400,NVAR=30,NFOR=18) CHARACTER*6 ISTACK(NSTKSIZ),IVAR CHARACTER*6 CHECK1, CHECK2 DIMENSION^STACK(NSTKSIZ) EXTERNAL^RGHTJST,CONVERTION INCLUDE^VAR1 C^Initialize ITOP = 0 C^Start traversing the input stack (N tokens) DO 10 I = 1, N C^Check for a operator IF(ISTACK(I) .EQ. '+^') THEN IF(ITOP .LT. 2) THEN WRITE(*,*) 'ERROR IN + OPERATION' WRITE(*,20) II STOP 'ERROR' ELSE VALUE = STACK(ITOP) + STACK(ITOP-1) 102 ITOP = ITOP - 1 STACK(ITOP) = VALUE ENDIF ELSE IF(ISTACK(I) .EQ. '-^') THEN IF(ITOP .GT. 1) THEN VALUE = STACK(ITOP-1) - STACK(ITOP) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE IF (ITOP .EQ. 1) THEN STACK(ITOP) = -STACK(ITOP) ELSE WRITE(*,*) 'ERROR IN - OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '*^') THEN IF(ITOP .GT. 1) THEN VALUE = STACK(ITOP-1) * STACK(ITOP) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN * OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '/^') THEN IF(ITOP .GT. 1) THEN IF(STACK(ITOP) .EQ. O.DO) THEN WRITE(*,*) 'DIVISION BY ZERO;' WRITE(*,20) II CALL ERROR('Error') STOP ENDIF VALUE = STACK(ITOP-1) / STACK(ITOP) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN / OPERATION' WRITE(*,20) II 103 ELSE ELSE ELSE ELSE STOP 'ERROR' ENDIF IF(ISTACK(I) . EQ. '** ^') THEN IF(ITOP .GT. 1) THEN IF(STACK(ITOP-1) .LT. 0.D0) THEN WRITE(*,*) 'EXPONENTIATION OF NEGATIVE NUMBER' WRITE(*,20) II CALL ERROR('Error') STOP ENDIF VALUE = STACK(ITOP-1) ** STACK(ITOP) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN ** OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'COS^') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DCOS(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN COSINE OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'SIN^') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DSIN(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN SINE OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'TAN^') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DTAN(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN TAN OPERATION' 104 ELSE ELSE ELSE ELSE WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'COSH ') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DCOSH(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN COSH OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'SQRT ') THEN IF(ITOP .GE. 1) THEN IF(STACK(ITOP) .LT. O.DO) THEN WRITE(*,*) 'SQUARE ROOT OF NEGATIVE NUMBER' WRITE(*,20) II CALL ERROR('Error') STOP ENDIF STACK(ITOP) = DSQRT(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN SQRT OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'EXP^') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DEXP(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN EXP OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'ACOS ') THEN IF(ITOP .GE. 1) THEN IF(STACK(ITOP).LT.-1.DO.OR.STACK(ITOP).GT.1.D0) THEN WRITE(*,*) 'ACOS ARGUMENT EITHER >1 OR <-1;' WRITE(*,20) II CALL ERROR('Error') 105 STOP ENDIF STACK(ITOP) = DACOS(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ACOS OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. 'ASIN ') THEN IF(ITOP .GE. 1) THEN IF(STACK(ITOP) .LT. O.DO) THEN WRITE(*,*) 'ASIN ARGUMENT EITHER >1 OR <1;' WRITE(*,20) II CALL ERROR('Error') STOP ENDIF STACK(ITOP) = DASIN(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ASIN OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. 'ATAN ') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DATAN(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ATAN OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. 'TANH ') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DTANH(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN TANH OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. 'SINH ') THEN 106 ELSE ELSE ELSE ELSE IF(ITOP .GE. 1) THEN STACK(ITOP) = DSINH(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN SINH OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'ALOG ') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DLOG(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ALOG OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'ALOG10') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DLOG10(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ALOG10 OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. 'ABS ^') THEN IF(ITOP .GE. 1) THEN STACK(ITOP) = DABS(STACK(ITOP)) ELSE WRITE(*,*) 'ERROR IN ABS OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF IF(ISTACK(I) .EQ. '.EQV. ') THEN IF(ITOP .GT. 1) THEN K = 1 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE 107 ELSE WRITE(*,*) 'ERROR IN .EQV. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.NEQV.') THEN IF(ITOP .GT. 1) THEN K = 2 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .NEQV. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.OR. ') THEN IF(ITOP .GT. 1) THEN K = 3 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .OR. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.AND. ') THEN IF(ITOP .GT. 1) THEN K = 4 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE 108 ELSE WRITE(*,*) 'ERROR IN .AND. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.NOT. ') THEN IF(ITOP .GE. 1) THEN K = 5 A = STACK(ITOP) CALL CONVERTION1(A,VALUE) STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .NOT. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.EQ. ') THEN IF(ITOP .GT. 1) THEN K = 6 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .EQ. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.NE. ') THEN IF(ITOP .GT. 1) THEN K = 7 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .NE. OPERATION' 109 WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.LE. ') THEN IF(ITOP .GT. 1) THEN K = 8 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .LE. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.LT. ') THEN IF(ITOP .GT. 1) THEN K = 9 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .LT. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.GE. ') THEN IF(ITOP .GT. 1) THEN K = 10 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .GE. OPERATION' 110 WRITE(*,20) II STOP 'ERROR' ENDIF ELSE IF(ISTACK(I) .EQ. '.GT. ') THEN IF(ITOP .GT. 1) THEN K = 11 A = STACK(ITOP-1) B = STACK(ITOP) CALL CONVERTION(A,B,VALUE,K) ITOP = ITOP - 1 STACK(ITOP) = VALUE ELSE WRITE(*,*) 'ERROR IN .GT. OPERATION' WRITE(*,20) II STOP 'ERROR' ENDIF ELSE C^It must be an operand, therefore place the operand on the C^top of the stack ITOP = ITOP + 1 READ(ISTACK(I),*,IOSTAT=NERROR) STACK(ITOP) IF(NERROR .GT. 0) THEN READ(ISTACK(I),'(1A6)') CHECK1 CALL RGHTJST(CHECK1) C^Check if the operand is a control system variable or a real number. C^If it is a variable, its value is in T_IVAR. DO 9 KK = 1, LL(II) CHECK2 = IVAR(II,KK) IF(CHECK1 .EQ. CHECK2) STACK(ITOP) = T_IVAR(II,KK) 9^CONTINUE ENDIF ENDIF 10^CONTINUE C^Check if ITOP = 1 IF( ITOP .EQ. 1) THEN 111 20 30 40 VALUE = STACK(1) ELSE WRITE(*,30) ITOP WRITE(*,20) II WRITE(*,40) N ENDIF FORMAT(' FORTRAN EXPRESSION it',I3) FORMAT(' ERROR IN THE END OF STACK - ITOP =',I3) FORMAT(' NTOKENS = ',I3) RETURN END Appendix B SS and EMTP Input Files B.1 Control System Case Studies — SS Input Files Only B.1.1 Second-Order linear system (Section 5.2.1) C Second-order linear system -- step unit response C C^s-block^1/s ^ 1 VAR2 + VAR1 - LOOP2 1.D0 1.D0 1.D0 C C^s-block -- 1/s ^ 1^OUT - LOOP1 + VAR2 1.D0 1.D0 1.D0 C C^k-block -- *16 0 VAR1 +SOURCE^ 16.D0 C C^k-block^auxiliar block 0 LOOP1 + OUT^ 1.D0 C C^k-block -- *16 0 LOOP2 + OUT^ 16.D0 C C^source block -- unit step function ^ 21SOURCE^1.D0 O.DO C C^fortran expression block 113 C C^output variables 33 OUT VAR1 C C^connection variables C C^initial conditions 77 VAR2^O.DO 77 OUT^O.DO C C^end of file B.1.2 Static and Dynamic Limiters (Section 5.2.2) C^Static limiter C C^k-block^auxiliar block INP + SRC23 + AUX C C^k-block^auxiliar block AUX + SRC21 C C^s-block^(10**3)/s 1 VAR1 + INP 1.D0 1.D0 C C^k-block -- limiter block OUT + VAR1^ 1.D0 -2.D0 2.D0 C C^source blocks -- input square wave 21 SRC21^-1.D0^ 23 SRC23^2.D0^20.D-3^10.D-3^ C C^fortran expression block C 114 O.DO O.DO C^output variables 33 INP OUT C C^connection variables C C^initial conditions 77 OUT^-4S.D-1 C C^end of file -- static limiter C^Dynamic limiter C C^k-block^auxiliar block INP + SRC23 + AUX^ 1.D0 C C^k-block^auxiliar block AUX + SRC21^ 1.D0 C C^s-block^(10**3)/s with limits (-.6,1.5) 1 OUT + INP 1.D0 1.D0 1.D3-5.D-1 1.5D0 C C^source blocks -- input square wave 21 SRC21^-1.D0^ 23 SRC23^2.D0^20.D-3^10.D-3^ C C^fortran expression block C C^output variables 33 INP OUT C C^connection variables 115 O.DO 0.D0 C C^initial conditions 77^OUT^-6.D-1 C C^end of file B.1.3 FORTRAN Function in a Closed Loop (Section 5.2.3) C FORTRAN function in a closed loop -- Work Book IV C C^s-block 1 VAR1 -FDBACK +SOURCE^ 6.D0^1.D0 1.D0 C C^k-block 0 VAR2 + VAR1 C C^k-block 0 OUT + VAR2 C C^k-block OFDBACK + OUT C C^source block 22SOURCE^1.D0^-60.D0^90.0 C C^Fortran expression block C C^output variables 33 OUT VAR2FDBACK C C^connection variables C 116 1.D0 1.D0 1.D0-2.D-3 2.D-3 25.D-2 ^ C^initial conditions C C^end of file B.1.4 Special Element Blocks (Section 5.2.4) C RMS, NONLINEAR FUNCTION, TIME DELAY AND FORTRAN EXPRESSION BLOCKS C &CASE37.OUT^(Output file name) &CASE37.AUX^(Auxiliar file name - shows how SS reads the input file) &YN^(Apply fast Newton Method)(Save number of iterations) C^k-block 0 VAR1 + SCR^ 1.D0 C C C^Time delay block 66 VAR5 VAR4^2.D-3 C C^source block 22^SCR^1.D0^60.D0 C C^Fortran expression block 99 VAR4 =VAR2*1.D0 C C^point-to-point nonlinear function --> y = x**2 88 VAR2^VAR1 C C^x^ y C -2.D0^4.D0 -1.05D0^1.1025D0 -1.D0^1.D0 -9.5D-01^9.03D-01 -9.0D-01^8.10D-01 -8.5D-01^7.22D-01 -8.0D-01^6.40D-01 -7.5D-01^5.62D-01 117 -7.0D-01 -6.5D-01 -6.0D-01 -5.5D-01 -5.0D-01 -4.5D-01 -4.0D-01 -3.5D-01 -3.0D-01 -2.5D-01 -2.0D-01 -1.5D-01 -1.0D-01 -5.0D-02 0.D0 6.00D-02 1.00D-01 1.50D-01 2.00D-01 2.50D-01 3.00D-01 3.60D-01 4.00D-01 4.50D-01 5.00D-01 5.60D-01 6.00D-01 6.50D-01 7.00D-01 7.50D-01 8.00D-01 8.50D-01 9.00D-01 9.50D-01 1.D0 1.05D0 1.1D0 2.D0 99999.D0 4.90D-01 4.22D-01 3.60D-01 3.02D-01 2.50D-01 2.02D-01 1.60D-01 1.22D-01 9.00D-02 6.25D-02 4.00D-02 2.25D-02 1.00D-02 2.50D-03 0.D0 2.50D-03 1.00D-02 2.25D-02 4.00D-02 6.25D-02 9.00D-02 1.22D-01 1.60D-01 2.02D-01 2.50D-01 3.02D-01 3.60D-01 4.23D-01 4.90D-01 5.63D-01 6.40D-01 7.23D-01 8.10D-01 9.03D-01 1.D0 1.1025D0 1.21D0 4.D0 C C^RMS block ^ 60.D0 60 VAR6 VAR5 C 118 C^output variables 33 VAR1 VAR4 VAR5 VAR6 C C^initial conditions ^ 77 VAR2 1.D0 ^ 77 VAR4 1.D0 C C^end of file B.2 Control and Power System Case Studies SS and EMTP Input Files B.2.1 Arc Modelling (Section 5.3.1) EMTP Input File &LIMA.INP^ * SS input file ^ Case identification card Resistor as a voltage contolled current source ^-1 ^ Time card .0001^.2 ^ Lumped RLC branch 1^2^ 300. 2^ 100. ^ 92^2 2.^0. 9999999. $>>>>>>>End of level 1: Linear and nonlinear elements<<<<<<<<< $>>>>>>>End of level 2: Switches and piecewise linear elements * ^ Voltage or current sources 11^1 1^1 $>>>>>>>End of level 3: Sources<<««<<<««<<<<<««««<<<<< 1^**** All voltages will be printed **** $>>>>>>>Level 5: End of data case<<<<<<<<<<<<<<<<<<<<<<< 119 11 SS Input File C^Arc modelling C Input file name &LIMA.OUT^ Auxiliar file name &LIMA.AUX^ &YN^(Apply fast Newton Method?),(Save number of iterations?) C C^k-block 0 ITACS + VEMTP^ 1.D-1 C C^output variables 33 ITACS C C^connection variables -- variables shared by SS and EMTP 44 ITACS VEMTP^0.D0 C C^initial conditions ^ 77 ITACS 0.D0 C C^end of file B.2.2 Single-Phase Variable Load (Section 5.3.2) EMTP Input File C^Single phase variable load C C &CASE17.0UT &CASE17.AUX &YN C C C^s-block 1 ITACS +^INT^ 1131.^1.D0 1.D0 C C 1.D0 s -block 120 1.D0 1^PB + PINST^ 1.D0 120.D0 C C^k-block 0 AUX1 + SRC21 C C^k-block 0 AUX2 + VEMTP C C^k-block 0 AUX5 + SCR23 C C^k-block 0 PAVER +^PB - PDEL C 377.D0 1.D0 1.D0 1.D0 C^time-delay block 66 PDEL^PB^8.33D-3 C C^source block 21 SRC21^1.D0 ^ ^ 23 SCR23^3393.D0 1.D-1 5.D-2^ C C^Fortran expression block 99 INT = AUX2*(AUX1+AUX5) 99 PINST = AUX2*ITACS C C^output variables 33 ITACS PAVER C C^connection variables 44 ITACS VEMTP^31623D-4 C C^initial conditions ^ 77^INT 0.D0 ^ 77 PINST 0.D0 ^ 77^PB O.DO C 121 5.D-2 C^end of file SS Input File &CASE17.INP ^ Case identification card Resistor as a voltage contolled current source Time card ^ 0.00005^0.2^ 1 * Lumped RLC branch 1^2^1.D-19 92^2 2.^O. 9999999. $>>>>>»End of level 1: Linear and nonlinear elements<<<<<<<<< $>>>>>>>End of level 2: Switches and piecewise linear elements ^ ^ Voltage or current sources 14^1 1^3.1623^60. $>>>>>>>End of level 3: Sources «««««««««««««««< 1^**** All voltages will be printed **** $>>>>>>>Level 5: End of data case<<<<<««««««««« B.2.3 Step Load Variation (Section 5.3.3) EMTP Input File &CASE22.INP C^.^ .^. Case identification card THREE-PHASE VARIABLE LOAD-STEP VARIATION C LOAD REPRESENTED INSIDE SS (NODE BUS13A) C C^ . Time card 50.E-6^.2^1 C C^.^ .^. Symmetric pi circuit 1 THEVA BUS7A^83.E-3^35. 2 THEVB BUS7B^23.E-3^-5.^83.E-3^35. 3 THEVC BUS7C^23.E-3^-5.^23.E-3^-5. C C^.^ . Lumped RLC branch 122 83.E-3^35. ^ BUS1A BUS1B BUS1C BUS12ABUS13A BUS12BBUS13B BUS12CBUS13C 92 IL13A 2.^0. 9999999. ^ 92 IL13B IL13A ^ 92 IL13C IL13A 4.011 4.011 4.011 70.16 70.16 70.16 C C^. .^. Constant-parameters line -1 BUS7A BUS1A .03167 3.222.00787 144.4 -2 BUS7B BUS1B .0243 .9238^.0126 144.4 -3 BUS7C BUS1C -1 BUS1ABUS12A .03167 3.222.00787 24.14 -2 BUS1BBUS12B .0243 .9238^.0126 24.14 -3 BUS1CBUS12C $ = = End of level 1: Linear and nonlinear elements ^ BUS13A IL13A^0.^99999. BUS13B IL13B^0.^99999. BUS13C IL13C^0.^99999. $ = = = End of level 2: Switches and piecewise linear elements C C^.^ .^. Voltage or current sources 14 THEVA 1 60. 187.79 14 THEVB 1 187.79 60. -120. 14 THEVC 1 187.79 60. 120. $ = = = End of level 3: Sources -1. -1 . -1. C C^ . Voltage-output nodes BUS1A $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case ^ SS Input File C^CASE22.INP C C^Phase A C C^s-block ^ 1ITACSA + INT_A 11468.D-1^1.D0 1.D0 123 ^ 1.D0 C C^k-block ^ O AUX2A +VEMTPA 1.D0 C C^Phase B C C^s-block ^ 1ITACSB + INT_B 11468.D-1^1.D0 1.D0 1.D0 C C^k-block ^ O AUX2B +VEMTPB 1.D0 C C^Phase C C C^s-block ^ 1ITACSC + INT_C 11468.D-1^1.D0 1.D0 1.D0 C C^k-block ^ O AUX2C +VEMTPC 1.D0 C C Auxiliar Blocks C C^k-block O CONS1 +SRC211^ 1.D0 C C^k-block O CONS2 +SRC212^ 1.D0 C C^k-block O CONS3 +SRC213^ 1.D0 C C^k-block O CONS4 +SRC214^ 1.D0 C C^k-block O AUX2 + CONS1 + CONS2 + CONS3 + CONS4 ^ C C^source block 124 1.D0 21SRC211^1795.D-4 21SRC212 23715.D-4 21SRC213 28650.D-4 21SRC214 12293.D-4 5.D-2 1.D-1 15.D-2 C C^fortran expression block 99 INT_A =AUX2A*AUX2 99 INT_B =AUX2B*AUX2 99 INT_C =AUX2C*AUX2 C C^output variables 33 CONS4 C C^connection variables 44ITACSAVEMTPA^0.D0 44ITACSBVEMTPB^O.DO 44ITACSCVEMTPC^0.D0 C C^initial conditions ^ 77 INT_A O.DO ^ 77 INT_B O.DO ^ 77 INT_C O.DO C C^end of file B.2.4 Pulse Load Variation (Section 5.3.3) EMTP Input File kCASE23.INP C^.^ .^. Case identification card THREE-PHASE VARIABLE LOAD-PULSATING LOAD VARIATION C LOAD REPRESENTED INSIDE SS (NODE BUS13A) C C^.^ . Time card 50.E-'6^.2^1 1 C C^.^ . Symmetric pi circuit 1 THEVA BUSTA^0.083334.563 125 ^ ^ ^ 2 THEVB BUS7B 0.083334.563 0.0233-5.427 ^ ^ ^ 0.083334.563 3 THEVC BUS7C 0.0233-5.427 0.0233-5.427 C C^.^ . Lumped RLC branch ^ BUS1A 4.011 ^ BUS1B 4.011 ^ BUS1C 4.011 ^ BUS12ABUS13A 70.16 ^ BUS12BBUS13B 70.16 ^ BUS12CBUS13C 70.16 92 IL13A 2.^0. 9999999. ^ 92 IL13B IL13A ^ 92 IL13C IL13A C C^.^ .^. Constant-parameters line .03167 3.222.00787 144.4 -1 BUS7A BUS1A -2 BUS7B BUS1B .0243 .9238^.0126 144.4 -3 BUS7C BUS1C -1 BUS1ABUS12A .03167 3.222.00787 24.14 -2 BUS1BBUS12B .0243 .9238^.0126 24.14 -3 BUS1CBUS12C $ = = End of level 1: Linear and nonlinear elements ^ BUS13A IL13A^0. 99999. BUS13B IL13B^0. 99999. BUS13C IL13C^0. 99999. $ = = = End of level 2: Switches and piecewise linear elements C C^. .^. Voltage or current sources 14 THEVA 1^187.79 60. 14 THEVB 1^187.79 60.^-120. 14 THEVC 1^187.79^60.^120. $ = = = End of level 3: Sources ^ C C^.^ . Voltage-output nodes BUS1A $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case ^ SS Input File C^CASE23.INP C C^Phase A 126 ^ ^ C C^s-block ^ 1ITACSA + INT_A 1.D0 11468.D-1^1.D0 1.D0 C C^k-block ^ 0 AUX2A +VEMTPA 1.D0 C C^k-block ^ O AUX3A +ITACSA 1.D0 C C^Phase B C C^s-block ^ 1ITACSB + INT_B 1.D0 11468.D-1^1.D0 1.D0 C C^k-block ^ O AUX2B +VEMTPB 1.D0 C C^k-block ^ O AUX3B +ITACSB 1.D0 C C^Phase C C C^s-block ^ 1ITACSC + INT_C 1 .DO 11468.D-1^1.D0 1.D0 C C^k-block ^ O AUX2C +VEMTPC 1.D0 C C^Auxiliary blocks C C^k-block ^ O AUX3C +ITACSC 1.D0 C C^k-block ^ O CONS1 + SRC21 1.D0 C C^k-block 127 ^ 0 CONS2 + SRC23^ C C^k-block ^ 0 AUX2 + CONS1 + CONS2 C C^source block 21 SRC21^1795.D-4 ^ ^ 5.D-2^ 23 SRC23^524.D-2 1.D-1 C 1.D0 1.D0 5.D-2 C^fortran expression block 99 INT_A =AUX2A*AUX2 99 INT_B =AUX2B*AUX2 99 INT_C =AUX2C*AUX2 C C^output variables 33 AUX2B C C^connection variables 44ITACSAVEMTPA^0.D0 44ITACSBVEMTPB^0.D0 44ITACSCVEMTPC^0.D0 C C^initial conditions ^ 77 INT_A 0.D0 ^ 77 INT_B O.DO ^ 77 INT_C O.DO C C^end of file B.2.5 Load Variation with Thyristor Controlled Reactor (Section 5.3.3) Without TCR EMTP Input File &CASE24.INP C 128 C^.^ .^. Case identification card THREE-PHASE VARIABLE LOAD WITHOUT TCR C LOAD REPRESENTED INSIDE Ss (NODE BUS13A) C . Time card C^.^ ^ 50.E-6^.2^1 1 C . Symmetric pi circuit C^.^ 1 THEVA BUS7A^0.083334.563 2 THEVB BUS7B^0.0233-5.427^0.083334.563 3 THEVC BUS7C^0.0233-5.427^0.0233-6.427^0.083334.563 C . Lumped RLC branch C^.^ 7. BUS1A^ BUS1B^ 7. 7. BUS1C^ BUS12ABUS13A^ 70.16 BUS12BBUS13B^ 70.16 BUS12CBUS13C^ 70.16 92 IL13A^ 2.^0. 9999999. 92 IL13B^IL13A^ 92 IL13C^IL13A^ 11 2 3 C C^.^ .^. Constant-parameters line -1 BUS7A BUS1A^.03167 3.222.00787 144.4 -2 BUS7B BUS1B^.0243 .9238 .0126 144.4 -3 BUS7C BUS1C -1 BUS1ABUS12A^.03167 3.222.00787 24.14 -2 BUS1BBUS12B^.0243 .9238 .0126 24.14 -3 BUS1CBUS12C $ = = End of level 1: Linear and nonlinear elements ^ BUS13A IL13A^O.^99999. BUS13B IL13B^O.^99999. BUS13C IL13C^O.^99999. $ = = = End of level 2: Switches and piecewise linear elements C . Voltage or current sources C^.^ 14 THEVA 1^187.79^60.^ 14 THEVB 1^187.79^60.^-120.^ 14 THEVC 1^187.79^60.^120.^ $ = = = End of level 3: Sources ^ C 129 -1. -1. -1. . Voltage-output nodes C^.^ BUS1A $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case ^ SS Input File C^CASE24.INP C &CASE24.0UT &CASE24.AUX &YN C C C^Phase A C C^s-block 1ITACSA + INT_A 11468.D-1^1.D0 1.D0 1.D0 C C^k-block O CONS1 +SRC211 1.D0 C C^k-block O CONS2 +SRC212 1.D0 C C^k-block 0 CONS3 +SRC213 1.D0 C C^k-block 0 CONS4 +SRC214 1.D0 C C^k-block O AUX2A +VEMTPA 1.D0 C C^k-block O AUX2 + CONS1 + CONS2 + CONS3 + CONS4 1.D0 C C^Phase B C C^s-block 1ITACSB + INT_B 11468.D-1^1.D0 1.D0 1.D0 130 C C^k-block ^ 0 AUX2B +VEMTPB 1.D0 C C^Phase C C C^s-block ^ 1ITACSC + INT_C 11468.D-1^1.D0 1.D0 1.D0 C C^k-block ^ 0 AUX2C +VEMTPC 1.D0 C C^source block 21SRC211^1795.D-4 21SRC212^O.DO^ 21SRC213^64658.D-4^ 21SRC214^O.DO^ 5.D-2 1.D-1 15.D-2 C C^fortran expression block 99 INT_A =AUX2A*AUX2 99 INT_B =AUX2B*AUX2 99 INT_C =AUX2C*AUX2 C C^RMS meter block 60 VRMS AUX2A^60.DO C C^output variables 33 VRMS C C^connection variables 44ITACSAVEMTPA^0.D0 44ITACSBVEMTPB^O.DO 44ITACSCVEMTPC^0.D0 C C^initial conditions 77 INT_A^ O.DO 131 ^ 77 INT_B 0.D0 ^ 77 INT_C 0.DO C C^end of file With TCR EMTP Input File &CASE25.INP C C^.^ .^. Case identification card THREE-PHASE VARIABLE LOAD WITH TCR AT BUS 1 C LOAD REPRESENTED INSIDE SS (NODE BUS13A) C C^. 50.E-6 .2 . Time card 1 1 C C^.^.^. 1 THEVA BUS7A 2 THEVB BUS7B 3 THEVC BUS7C . Symmetric pi circuit 0.083334.663 0.0233-5.427^0.083334.563 0.0233-5.427^0.0233-5.427 0.083334.563 C C^. BUS1A BUS1B BUS1C BUS12ABUS13A BUS12BBUS13B BUS12CBUS13C IL13A IL13B IL13C IX13A IX13B IX13C 92^TCRA 2.^O. 9999999. 92^TCRB TCRA 92^TCRC TCRA . Lumped RLC branch 7. 7. 7. 70.16 70.16 70.16 6389.^5571. 6389.^5571. 6389.^5571. 177.36154.65 177.36154.65 177.36154.65 1 2 3 C C^. .^. Constant-parameters line -1 BUS7A BUS1A^.03167 3.222.00787 144.4 132 -2 BUS7B BUS1B .0243 .9238^.0126 144.4 -3 BUS7C BUS1C -1 BUS1ABUS12A .03167 3.222.00787 24.14 -2 BUS1BBUS12B .0243 .9238^.0126 24.14 -3 BUS1CBUS12C $ = = End of level 1: Linear and nonlinear elements BUS13A IL13A 0. .1 .0001 BUS13B IL13B 0. .1 .0001 BUS13C IL13C 0. .1 .0001 .0001 BUS13A IX13A .1 999999. BUS13B IX13B .1 999999. .0001 .0001 BUS13C IX13C .1 999999. BUS1A^TCRA 0. 99999. BUS1B^TCRB 0. 99999. BUS1C^TCRC 0. 99999. $ = = =^End of level 2: Switches and piecewise linear elements ^ C C^. 14 THEVA 14 THEVB 14 THEVC $ = = = . Voltage or current sources 1^187.79 60. 1^187.79 60. -120. 1^187.79 60. 120. End of level 3: Sources -1. -1. -1. C . Voltage-output nodes BUS1A $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case ^ C SS Input File C^CASE25.INP C &CASE25.0UT &CASE25.AUX &YN C C C^s-block 300/s 1 EPU +VRMSPU - UNITY 1.D0 1.D0 300.D0 C C C^k-block ^ 0 BLGAM + EPU 125.D-2 0.D0125D-2 133 C C 7630823661.D-12 OVRMSPU + VRMS C C^k-block O UNITY + SRC21 C C^k-block O EMTP1 +VEMTP1 C C^k-block O EMTP2 +VEMTP2 C C^k-block 0 EMTP3 +VEMTP3 C C^Phase A - Busi C^s-block 1/s 1ITACS1 + BGAMA 1.D0 1.D0 C C^Phase B - Busi C^s-block^1/s 1ITACS2 + BGAMB 1.D0 1.D0 C C^Phase C - Bus1 C^s-block 1/s 1ITACS3 + BGAMC 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 C C^source block 21 SRC21^1.D0 C C^fortran expression block 99 BGAMA =EMTP1*BLGAM 99 BGAMB =EMTP2*BLGAM 99 BGAMC =EMTP3*BLGAM C 134 ^ C^RMS meter block 60 VRMS EMTP1^60.DO C C^output variables 33VRMSPU C C^connection variables 44ITACS1VEMTP1^0.D0 44ITACS2VEMTP2^0.D0 44ITACS3VEMTP3^0.D0 C C^initial conditions 77 BGAMA^ 0.D0 77 BGAMB^ 0.D0 77 BGAMC^ 0.D0 77 EPU^ 0.D0 77ITACS1^ 0.D0 77ITACS2^ 0.D0 77ITACS3^ 0.D0 C C^end of file B.2.6 HVDC System Simple Model (Section 5.3.4) EMTP Input File &CASE40.INP * * Case identification card ^ CSTP alone - Time delay block - 2nd case ^ Time card 5.D-5^.2^2^1^ ^ Lumped RLC branch ^ SCR BUS1 1.1572 44.58 ^ BUS1 9.66 92 BUS1 2.^0. 135 1 9999999. $>>>>>>>End of level 1: Linear and nonlinear elements<<<<<<<<< $>>>>>>>End of level 2: Switches and piecewise linear elements ^ Voltage or current sources 14^SCR 1^187.794^60.^67.43^ -1.0 $>>>>>>>End of level 3: Sources«<<<<<<<<<<<<<<<<<<<<<<<<<<<<< BUS1 $>>>>>>>Level 5: End of data case<<<<<<<««<<<<«««<< SS Input File C CASE40.INP C BASIC HVDC MODELS - FIG. 7-10 - WORKBOOK IV - CHAP 7 - PAGE 7-12 C Input file name &CASE40_1.0UT^ &CASE40_1.AUX^ Auxiliar file name &YN^(Apply fast Newton Method?),(Save number of iterations?) C C C^k-block O IDC +SCR211 + AUX3 + DCF^ 1.D0 C C^k-block O X + AUX2 +SCR212^ 1.D0 C C^k-block O ALPHA + CONT +SCR215^ 1.D0 C C^k-block ^ 0 AUX3 + SCR24 1.D0 C C^k-block ^ 1.D0 O AUX4 +SCR216 C C^k-block ^ O BUS1 + VEMTP 1.D0 C C^k-block ^ O AUX1 +^PI 1.D0 C C^k-block ^ XC +SCR213 1.D0 O C C^k-block 136 ^ O 1.D0 A +SCR214^ C C^s-block^120/s 1^C1P + PRODC^ 1.D0 1.D0 120.D0 C C^s-block^120/s 1 S1P + PRODS 1.D0 1.D0 120.D0 C C^k-block O C1 +^C1P -^C1D^ 1.D0 C C^k-block S1 +^S1P -^S1D^ O 1.D0 C C^k-block O AUX2 + TIMEX^ 1.D0 C C^Time delay block 66^C1D^ClP^1667.D-5 66^S1D S1P^1667.D-5 C C^source block 21SCR211^1.6 21SCR212^-0.1 21SCR213^9.08 21SCR214^0.8933 21SCR215^0.31416 21SCR216^107.79 24 SCR24^1.6 25 TIMEX 25^PI 50.D-3^150.D-3 50.D-3^ C C^Fortran expression blocks 99 DCF^=1.D6*AUX4*(ABS(XP)**3.15)*EXP((0.0-286.)*XP) 99 CONT =3.6*EXP((0.0-1.0)*135.*XP)-4.822*EXP((0.0-217.)*XP)+1.222 99^XP =(ABS(X))*(X.GE.0) 99^at =120*AUX1*AUX2 137 50.D-3 99 PRODC =BUS1*COS(wt) 99 PRODS =BUS1*SIN(wt) 99^VM =SQRT(C1*C1+S1*S1) 99^CV =(C1*COS(wt)+S1*SIN(wt))/VM 99^SV =(C1*SIN(wt)-S1*COS(wt))/VM C C^DC system equation C 99 PHI^=ACOS(COS(ALPHA)-(XC*IDC)/(A*SQRT(3)*VM)) C C^AC current injection C 99 IACF =(CV*COS(PHI)+SV*SIN(PHI))*A*2.0*SQRT(3)*IDC/AUX1 C C^output variables 33^C1^CV IDC ALPHA C C^connection variable 44 IACF VEMTP^77.614 C C^initial conditions 77^C1P^ 1.D0 77^S1P^ O.DO 77 DCF^ 0.D0 77 CONT^ O.DO 77 PRODC^77.614 77 PRODS^ O.DO 77^VM^ 3.6 77^CV^ 1.D0 77^SV^ 0.D0 77 PHI^.31434D0 77 IACF^ 0.D0 77^XP^ 0.D0 77^wt^ O.DO C C^end of file B.2.7 Single-Phase Thyristor Rectifier Bridge (Section 5.4) EMTP Input File 138 ^ &MOHAN5.INP C^.^ .^. Case identification card 1-PHASE THYRISTOR RECTIFIER C . Time card C^.^ 50.0E-6^50.0E-3 ^ .5E-3 C . Lumped RLC branch C^.^ VSA VMA^ 0.2 VMA^VA^ 1.0 POSP^POS^ 16.0 POS NEG^2.0 POSP NEG^1.0E9 POSP^VA^200.0^1.0 POSP^POSP^VA VA NEG POSP^VA NEG POSP^VA POSA POSP^0.01 POSG POSP^0.01 ^ 92 VSA 4. 0. 9999999. 92 POS VSA $ = = End of level 1: Linear and nonlinear elements 1 2 C C^.^ . Time-controlled switch -1^VA POSA^ 0.002^60. -1^NEG^ 0.002^60. -1^POSG^ 0.002^60. -1^NEG^VA^ 0.002^60. $ = = = End of level 2: Switches and piecewise linear elements C C^.^ .^. Voltage or current sources 14^VSA^169.7^60.0^-90.0 $ = = = End of level 3: Sources C^.^.^.^.^.^.^. Voltage-output nodes VSA POSP $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case SS Input File C^MOHAN5.INP C 139 1 2 3 4 Input file name &MOHAN5.0UT^ Auxiliar file name &MOHANS.AUX^ &YY^(Apply fast Newton Method?),(Save number of iterations?) C C^k-block 0 VSA + VEMTP^ OVCONTL +SRC211^ CO AUX1 + AUX^ OPLS1-2 + AUX^ 1.D0 1.D0 1.D0 1.D0 C C^Time delay blocks 66PLS3-4PLS1-2 0.00833333 C66 AUX2 AUX1^50.D-6 C C^controlled integrator block ^ 58 RAMP+SRC212 VSA C C^source blocks 21SRC211^0.25 21SRC212^1. C C^fortran expression 99 AUX = RAMP .GE. VCONTL C C^output variables 33 RAMPPLS1-2PLS3-4 C C^connection variables 44 ICSTP VEMTP^0.D0 55PLS1-2 55PLS1-2 55PLS3-4 S5PLS3-4 C C^end of file 140 ^ 120.^1. ^ B.2.8 Three-Phase Thyristor Inverter (Section 5.4) EMTP Input File &MOHAN8.INP . Case identification card C^.^ 3-PHASE THYRISTOR INVERTER C . Time card 50.0E-6 100.0E-3^ C^ .5E-3 C . Lumped RLC branch C^.^ VSA VMA^ 0.2 VSB VMB VSA VMA VSC VMC VSA VMA VMA^VA^ 1.0 VMB^VB VMA^VA VMC^VC VMA^VA POSP^POS^ 16.0^ POS^NEG^1.0E-6 NEG^VA^1.0E5 NEG^VB^1.0E5 NEG^VC^1.0E5 POSP NEG^1.0E9^ VSA VSC^1.0E9^ 92 VSA^ 4. 4. 9999999. 92 VSC VSA $ = = End of level 1: Linear and nonlinear elements 1 2 2 1 2 C . Time-controlled switch C^.^ -1^VA POSP^ 2.E-3^60. -1^NEG^VC^ 2.E-3^60. -1^VB POSP^ 2.E-3^60. -1^NEG^VA^ 2.E-3^60. -1^VC POSP^ 2.E-3^60. -1^NEG^VB^ 2.E-3^60. $ = = = End of level 2: Switches and piecewise linear elements C .^. Voltage or current sources 14^VSA^391.9^60.0^-60.0 14^VSB^391.9^60.0^-180.0 14^VSC^391.9^60.0^-300.0 C^ 141 1 2 3 4 5 6 ^ 11^POS-1 -600.0E6 11^NEG-1^600.0E6 $ = = = End of level 3: Sources ^ .^. Voltage-output nodes C^ VSA $ = = = End of level 4: User-defined voltage output $ = = = Level 5: End of data case ^ SS Input File C^MOHAN8.INP C C MOHAN8_1.PLO - Simultaneous solution C MOHAN8_2.PLO - One time-step delay C MOHAN8_3.PLO - Two time-step delay C Input file name &MOHAN8.0UT^ Auxiliar file name &MOHAN8.AUX^ &YN^(Apply fast Newton Method?),(Save number of iterations?) C C^k-block 0^VSA +VEMTPA 0^VSC +VEMTPC 0^VSAC +^VSA OVCONTL +SRC211 0^AUX1 +^AUX 0^AUX3 +^AUX2 0 PULS1 +^AUX4 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 1.D0 VSC C C^Time delay 66 PULS2 PULS1 66 PULS3 PULS2 66 PULS4 PULS3 66 PULS5 PULS4 66 PULSE PULS5 66^DCMP^AUX1 66^AUX4^AUX3 blocks 0.00277778 0.00277778 0.00277778 0.00277778 0.00277778 250.D-6 100.D-6 C C^controlled integrator block ^ 58 RAMP+SRC212 VSAC C C^source blocks 142 ^ 120.^1. 21SRC211 .833333333 21SRC212^1. C C^fortran expression 99 AUX = RAMP .GE. VCONTL 99 AUX2 = (.NOT. DCMP) .AND. AUX1 C99 GATE1 = PULS1 .OR. PULS2 C99 GATE2 = PULS2 .OR. PULS3 C99 GATE3 = PULS3 .OR. PULS4 C99 GATE4 = PULS4 .OR. PULS5 C99 GATES = PULS5 .OR. PULS6 C99 GATE6 = PULS6 .0R. PULS1 C C^output variables 33 PULS1 PULS2 PULS3 PULS4 PULS5 PULS6 RAMP C C^connection variables 44ICSTAPVEMTPA^0.D0 44ICSTCPVEMTPC^0.D0 55 PULS1 55 PULS2 55 PULS3 55 PULS4 55 PULS5 55 PULS6 C C^initial conditions ^ 77 AUX O.DO ^ 77 AUX2 O.DO C C^end of file 143
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical instabilities in power system transient simulations
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical instabilities in power system transient simulations Araújo, António E. A. 1993
pdf
Page Metadata
Item Metadata
Title | Numerical instabilities in power system transient simulations |
Creator |
Araújo, António E. A. |
Date Issued | 1993 |
Description | In power system simulations, control systems act upon the electric power system in a variety of ways. There is, therefore, a necessity of simulating the control system equations together with the power system equations. This is the case not only instability studies, but also in the simulation of electromagnetic transient phenomena. Because of the different structure of the two sets of equations, computational efficiency can be improved if both sets of equations are solved separately, with the result from one set of equations showing up in the other set of equations one or more time steps later. This procedure, however efficient, can lead to numerical instability. In the simulation of the control system itself, the modelling of nonlinearities in closed loops poses some problems as well. The most efficient way of handling a nonlinear block is to open the closed loop and use, instead of the current value of the variable, a predicted value of the feedback path. This can also lead to numerical instabilities. In this thesis a simultaneous solution of the control and power system equations is proposed. For the cases of nonlinearities in closed loops, an iterative solution is suggested. The implementation of the method is carried out in the Electromagnetic Transients Program (EMTP). Several comparisons between the proposed and other methods are performed. Incases in which the analytical solution is not available, the comparison is made against the results of the TAGS package ("Transient Analysis of Control Systems") which is integrated with the EMTP version of DCG/EPRI ("Development Coordination Group/Electric Power Research Institute"). It is shown through simple, but typical, case studies that the proposed method is not only accurate and efficient, but also numerically stable, thus representing a significant improvement to the ability of the EMTP to reliably simulate electromagnetic transients in power systems when control systems are involved. |
Extent | 4636282 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-09-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065142 |
URI | http://hdl.handle.net/2429/2121 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1993-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1993_fall_phd_araujo_antonio.pdf [ 4.42MB ]
- Metadata
- JSON: 831-1.0065142.json
- JSON-LD: 831-1.0065142-ld.json
- RDF/XML (Pretty): 831-1.0065142-rdf.xml
- RDF/JSON: 831-1.0065142-rdf.json
- Turtle: 831-1.0065142-turtle.txt
- N-Triples: 831-1.0065142-rdf-ntriples.txt
- Original Record: 831-1.0065142-source.json
- Full Text
- 831-1.0065142-fulltext.txt
- Citation
- 831-1.0065142.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065142/manifest