FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 33, N o 1, March 2020, pp. 119-131 https://doi.org/10.2298/FUEE2001119M © 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND IMPROVING TIME INTEGRATION SCHEME FOR FET ANALYSIS OF POWER SYSTEM ANGLE STABILITY  Marin Mandić, Ivica Jurić-Grgić, Nedjeljka Grulović-Plavljanić University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Split, Croatia Abstract. This paper presents improved algorithm for numerical analysis of power system angle stability achieved by improvement of the time integration when forming a local system of equations for power system finite elements (FE). Previously developed local system of equations of power system angle stability has been obtained using the generalized trapezoidal rule (ϑ - method). Improvement of accuracy was obtained by using Heun's method. Numerical solutions obtained using Heun’s method and using the generalized trapezoidal rule are compared to Electromagnetic Transients Program (EMTP). It has been shown that Heun’s method yields the results with much higher accuracy comparing to results obtained by generalized trapezoidal rule. Key words: Heun’s method, finite element technique, angle stability, time domain analysis. 1. INTRODUCTION The assessment of angle or transient stability of an electrical system plays a fundamental role, as it allows contingencies classification and provides indications for the design and planning of the power systems. Transient stability or also known as large-disturbance rotor angle stability concerns the ability of the system to withstand critical disturbances such as three-phase short circuits. As it is reported in [1-7], it is very difficult to analyze simultaneously the slow electromechanical and fast electromagnetic phenomena. This paper presents improved algorithm for numerical analysis of power system angle stability. Complex numerical analysis is performed simultaneously by both the time- domain and the frequency-domain. The developed numerical model is based on a finite element technique (FET) procedure and time-varying phasors. In contrast to existing transient stability programs (TSP) that are limited to the first swing stability analysis, the developed numerical model can analyze large-disturbance rotor angle stability. The basis of developed numerical model for analysis of power system angle stability is application of FET to the power system so that the considered system is divided into smaller parts, Received August 28, 2019; received in revised form October 2, 2019 Corresponding author: Marin Mandić University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Split, Croatia (E-mail: mmandi00@fesb.hr)  120 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ which are treated as separated finite elements (FE) with a certain number of local nodes during entire time-domain simulation. Fundamentals of FET analysis of power system angle stability have been described in [8] where the generalized trapezoidal rule (ϑ - method) is used for temporal integration when forming a FE equations system. The use of the generalized trapezoidal rule for the time integration when forming a FE equations system can sometimes cause numerical oscillations or numerical diffusion of results. Numerical oscillations and/or numerical damping of results are the consequences of the choice of the numerical integration parameter (ϑ). In principle, its choice is based on numerical experiments. The influence of ϑ on stability of one-step time integration methods for the first order initial value problems is discussed in [9-10], where the stability behaviour is investigated by examining the free response of the FE model. It was found that methods for which 0.5 1  are unconditionally stable. The aim of this article is to improve the algorithm for numerical analysis of power system angle stability by improving the time integration scheme when forming FE equations system. Instead of using the generalized trapezoidal rule when forming a FE matrix, the accuracy of a numerical solution can be improved using one of the variants of the Runge- Kutta second-order method also known as Heun's method [11] for the time integration. Sometimes Heun‟s method is also referred as the improved Euler‟s method. Heun‟s method is suitable because, as a single-step method, it connects variables only at the beginning and at the end of the time interval and thus preserves the simplicity of numerical procedure. It is a single-step predictor-corrector method that improves the estimation of the slope for the time interval using an average value of slopes at the beginning and at the end of the time interval. Desire for calculation speed must be balanced with the appropriate accuracy and stability required. We could get higher accuracy by multi-stage Runge-Kutta second-order method but it leads to time-consuming problem. Theoretically, one can develop infinite number of time schemes and their combinations, but the idea was to improve accuracy of FET model computation for power system angle stability analysis with one-step method that is not time-consuming. Numerical solutions obtained using Heun‟s method and using the generalized trapezoidal rule for time integration parameter ϑ = 0.66 are compared to Electromagnetic Transients Program (EMTP). It has been shown that Heun‟s method yields the results with much higher accuracy comparing to results obtained by generalized trapezoidal rule. 2. FET ANALYSIS - HEUN‟S METHOD Fundamentals of FET analysis of power system angle stability have been described in [8] where the generalized trapezoidal rule (also known ϑ - method) is used for temporal integration when forming a FE equations system. In the following text, it will be shown how to derive a local system of equations for power system finite elements by using Heun‟s method for temporal integration instead of the generalized trapezoidal rule. Numerical model of turbine governor and excitation system have been included into synchronous generator FE model where numerical models have been improved using Heun's method [11] for the time integration. Finite elements of utility equivalent source, three-phase transformers and transmission lines are defined by the system of algebraic equations as it is already described in [8]. Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 121 2.1. FE Synchronous Generator – Heun’s method According to [12-13], salient pole synchronous generator subtransient phasor model is defined by the following system of algebraic and differential equations where time- domain swing equation (6) has been also included: q q q d dE U R I X I      (1) d d d q qE U R I X I      (2) o( ) o do f d d d dE T E X X I E dt          (3) ( ) q do o d d d q dE T E X X I E dt            (4) d d ( ) qo q q q dE T X X I E dt        (5) = P - P i m e d M dt   (6) = = ω - ω n d dt   (7) 1 1sin 2 d qo e d d q X XE U U P = δ + sin2δ X 2 X X             (8) where E'o, E"q, E"d are transient and subtransient generator voltages. Variable Ef represents generator field voltage while variables Uq, Ud, Iq and Id are axial components of the terminal voltage U and stator current I . In order to perform the time integration according the Heun's method, it is necessary to put relations (3-7) into following form: 1 ( , , ) o f d o dE f E I E dt   (9) 2 ( , , ) q o d q dE f E I E dt    (10) 3 ( , ) d q d dE f I E dt   (11) 4 ( ) ( , ) m e d f P P dt   (12) 5 ( ) d f dt    (13) where: 122 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ 1 ( ) ( , , ) f d d d o f d o do do do E X X I E f E I E T T T           (14) 2 ( ) ( , , ) qo d d d o d q do do do EE X X I f E I E T T T             (15) 3 ( ) ( , ) q q q d q d qo qo X X I E f I E T T         (16) 4 ( , ) m e m e i i P P f P P M M   (17) 5 ( )f     (18) According to Heun's method, the state of variables at the end of the time interval is defined by the following corrector equations: 1 1 ( , , ) ( , , ) 2 2 p o o f d o f d o t t E E f E I E f E I E             (19) 2 2 ( , , ) ( , , ) 2 2 p q q o d q o d q t t E E f E I E f E I E               (20) 3 3 ( , ) ( , ) 2 2 p d d q d q d t t E E f I E f I E            (21) 4 4 ( , ) ( , ) 2 2 m e m e t t f P P f P P           (22) 5 5 ( ) ( ) 2 2 t t f f            (23) According to (14-18), the predicted slopes 1 ( , , ) p f d o f E I E    , 2 ( , , ) p o d q f E I E     , 3 ( , ) p q d f I E   , 4 ( , )m ef P P   , 5 ( )f   at the end of the time interval are: 1 ( ) ( , , ) p p f d d d o f d o do do do E X X I E f E I E T T T                 (24) 2 ( ) ( , , ) p p qo d d d o d q do do do EE X X I f E I E T T T                 (25) 3 ( ) ( , ) p p q q q d q d qo qo X X I E f I E T T              (26) 4 ( , ) m e m e i i P P f P P M M       (27) 5 ( )f        (28) Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 123 where p o E  , p q E  , p d E  are predicted transient and subtransient voltage values at the end of the time interval defined by the following predictor equations: 1 ( , , ) p o o f d o E E t f E I E      (29) 2 ( , , ) p q q o d q E E t f E I E       (30) 3 ( , ) p d d q d E E t f I E      (31) where 1 ( , , )f d of E I E , 2 ( , , )o d qf E I E , 3 ( , )q df I E are defined by equations (14-16). The time integration of equations (3-7) using the generalized trapezoidal rule (  - method) can be found in [8]. With application of the Heun‟s method [11] on Equations (3-7), which includes combining (19-31) and separating the variables at the end of time interval to the left-hand side and the variables at the beginning of time interval to the right-hand side, the following system of algebraic equations is obtained: 2 2 2 2 2 2 ( ) 2 2 ( ) 2 2 2 2 ( ) ( ) ( ) ( ) 2 ( ) 2 ( ) 2 ( ) f d d d o o do do fo d d d o do do do do f d d d o do do do E t t I X X E E T T E tE t t I X X E t T T T T E t t I X X E t T T T                                                      (32) 22 2 2 2 2 ( ) 2 2 ( ) 2 2 2 2 ( )( ) ( ) ( ) 2 ( ) 2 ( ) 2 ( ) o d d d q q do do q qo d d d do do do do qo d d d do do do E t t I X X E E T T E t E tE t t I X X T T T T E tE t t I X X T T T                                                        (33) 2 2 2 2 ( ) 2 2 ( ) ( ) 2 2 2 ( ) ( ) ( ) 2 ( ) q q q d d d qo qo q q q d d qo qo qo q q q qo t I X X E t E E T T t I X X E t E t T T T t I X X T                                        (34) 2 2 2 2 m e m e i i i i P t P t P t P t M M M M                  (35) 2 2 t t        (36) 124 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ The variables in Equations (32-36) marked by “+” denote variables at the end of the time interval, while variables without mark denote variables at the beginning of the time interval. According to [8], synchronous generator can be defined as FE with three local nodes where FE local system has been defined by the system of algebraic Equation (8) and Equations (32-37). 1 1 { } [ ] { } [ ] { } g e g e o I Z Z E         (37) Dynamic phasors of local matrix as well as dynamic phasors of vectors in Equation (37) are defined as follows: 0 0 [ ] 0 0 0 0 e e e e Z Z Z Z             1 1 12 / 3 2 / 3 01 02 03 0 0 0 { } { } { } j j jT T o E E E E E e E e E e                  1 2 3 { } { } T g      1 2 3 { } { } T g I I I I where: [ cos ( )] j e d q d Z R j X e X X          2 2 = (E ) +(E ) o d q E   2 2 I= ( ) +( ) d q I I  =  +  - angle between stator current phasor I and „q‟ axis - power angle 2.2. Model of Hydraulic Digital Turbine Governor – Heun’s Method According to [14], numerical model of digital turbine governor is defined by the following system of algebraic and differential equations: s ref r RP = P Δω P  (38) s1 sa sb scP = P + P +P (39) sa p sP k P  (40) 1g P q    (41) _v open v v_closeq q q  (42) min maxq q q  (43) sc s sc d d dP dP P T = k dt dt   (44) Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 125 v g v g p s1 dq T q T T = P dt    (45) v dq = q dt (46) sb i s dP =k P dt  (47) 2 2 g g r r dP dq P T =δ T dt dt       (48) ( ) m m 11 w 11 23 11 13 21 w dP dq P a T = a δ a a a a T dt dt           (49) where:  - rotor angular frequency n - nominal angular frequency Pm, Pe - mechanical and electrical power kp - PID controller proportional gain ki - PID controller integral gain kd - PID controller derivative gain Td - PID controller derivative time constant Pref - referent value of mechanical power Tp - pilot time constant Tr - transient droop time constant Tw - water inertia time constant Tg - gate time constant qv_open - gate max opening speed qv_close - gate min closing speed qmax - maximal gate position qmin- minimal gate position  - permanent droop coefficient a11, a13, a21 - turbine coefficients a23 - turbine gain The time integration of above equations using the generalized trapezoidal rule (  - method) can be found in [8]. Using the same procedure shown in chapter 2.1. for time integration of Equations (44 - 49), that is application of Heun‟s method, the following system of algebraic equations is obtained: 2 2 2 2 )(2 )( )(2 )( 2 222 d sc d snd d sc d snd d sc d snd scsc T tP T tPk T tP T Pkt T tP T Pkt PP                      (50)       snsnss P t P t PP 22 (51) 126 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ 2 2 2 2 111 )(2 )( 2 )( 2222 p v pg s p v pg s p v pg s vv T tq TT tP T tq TT Pt T tq TT Pt qq                     (52)       vv q t q t qq 22 (53)       sisisbsb Pk t Pk t PP 22 (54) 2 2 2 2 2 2 22 2 )( 2 )( 2 222 r g r v r g v r g vgg T tP T tq T Pt q t T tP q t PP                       (55) 23 23 11 13 21 11 11 11 23 23 11 13 21 11 11 11 2 2 2 23 23 11 13 21 2 2 2 11 11 ( ) 2 2 2 ( ) 2 2 2 ( ) ( ) ( ) ( ) 2 ( ) 2 2 v m m m w w v m w w v m w w q t a t q a a a a P t P P a T a a T q a t t a a a a q P t a T a a T q a t t q a a a a P t a T a T                                                           2 11 ( ) w a T  (56) The numerical model of turbine governor is defined by the system of algebraic Equations (38-43) and Equations (50-56). 2.3. Model of Excitation System – Heun’s Method According to [15], numerical model of excitation system is defined by the following system of algebraic and differential equations: 0f C REF T f a E v v v v k    (57) in maxRm R Rv v v  (58) 1T r T t r r dv k v v dt T T     (59) 1 aR R C a a kdv v v dt T T     (60) 1f e f R e e dE k u v dt T T     (61) f f f f f dv dE v T k dt dt    (62) where: Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 127 ka, kr, kf, ke - regulator amplification factors Ta, Tr, Tf, Te - time constants vT, vf, vR , vC - regulator signals Ef - field voltage Ef 0 - initial field voltage vt - phase voltage vREF - referent phase voltage The time integration of above equations using the generalized trapezoidal rule (ϑ - method) can be found in [8]. Using the same procedure shown in chapter 2.1 for time integration of Equations (59 - 62), that is application of Heun‟s method, the following system of algebraic equations is obtained: 2 2 2 2 2 )( 2 )( 2222 r T r tr r T r tr r T r tr TT T tv T tvk T tv T vkt T tv T vkt vv                     (63) 2 2 2 2 2 )( 2 )( 2222 a R a ca a R a ca a R a ca RR T tv T tvk T tv T vkt T tv T vkt vv                     (64) 2 2 2 2 2 )( 2 )( 2222 e f e Re e f e Re e f e Re ff T tE T tvk T tE T vkt T tE T vkt EE                     (65) 2 2 2 2 2 )( 2 )( 1 2 2 1 22 1 f f f R e e f e f f f f R e e f e f f f f R e e f e f ff T tv T tv T k E T k T vt T v T k E T kt T tv T v T k E T kt vv                                        (66) The numerical model of excitation system is defined by the system of algebraic Equations (57-58) and Equations (63-66). 3. TEST EXAMPLE In order to verify a new numerical procedure, angle stability of the power system, shown in Figure 1, has been analysed. At the beginning of simulation, the generator 1 and generator 4 are in operating mode with P=108 MW (0.9 pu) and Q =0 Mvar (0 pu) while generator 2 and generator 3 are in operating mode with P=84 MW (0.7 pu) and Q =52.3 Mvar (0.436 pu). Global system of equations is obtained by assembling local system of equations of each FE according to the standard assembling procedure [9-10]. At the moment t=0.8 s, the three-phase fault with clearing fault time of 200 ms has been set at the busbar „D‟ causing synchronous generators rotor angle (δ) oscillations and variation of mechanical powers (Pm) due to turbine governor control. 128 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ The generators data are: 120 MVA n S  , 108 MW n P  , 52.3 MvarnQ  , 14.4 kVnU  , 4811 A n I  , 690 A fo I  , 1192 AfnI  . Available transient and subtransient generators time constants and reactances are: 1.01dx  pu, 0.666qx  pu, 0.287qx  pu, 0.421dx  pu, 0.268dx  pu, 0.198x  pu, 0.00236r  pu, 3.4 sdT  , 0 7.5 sdT  , 0.0554 sdT , 0.0876 sqT  , 8.63 s mG T  . The turbine governors time constants, permanent and temporary droop coefficients, maximal and minimal gate opening speed, maximal and minimal gate position, amplification factors and turbine coefficients are: Tp = 0.02 s, Tr = 4.5 s, Tw = 2.9 s, Tg = 0.5 s, _ 0.15 v open q  , v_close 0.08q  max 0.8q  , min 0.1q  , 0.04  , δ = 0.378, kp = 15, ki = 0.8, kd = 1.5, Td = 5 s, a11= 0.5, a13 = 1, a21 = 1.5, a23 = 1. The excitation systems time constants and amplification factors are: 0.02 s r T  , 0.001 s a T  , 0.1 s e T  , 0.1 sfT  , 1rk  , 50ak  , 1ek  , 0.001fk  . Fig. 1 Single line diagram of the electric power system Step-up transformers as well as transmission lines reactance‟s have been taken 0.1 pu. To demonstrate that Heun‟s method is more accurate than ϑ - method, numerical solutions obtained using Heun‟s method and using the generalized trapezoidal rule for time integration parameter ϑ = 0.66 are compared to Electromagnetic Transients Program (EMTP) as it shown in Figures (2-5). The below results, showing rotor angles (δ) and mechanical powers (Pm) of generators 1 and 2, clearly demonstrate that Heun‟s method yields the results with much higher accuracy comparing to results obtained by generalized trapezoidal rule. Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 129 Fig. 2 Generator „1‟ rotor angle (δ) during angle stability simulation Fig. 3 Generator „2‟ rotor angle (δ) during angle stability simulation Since different values of ϑ yield different time-stepping schemes, all these schemes vary in accuracy. As it can be seen from Figures (2-5), numerical solutions obtained by Heun‟s method yields highly accurate results while the generalized trapezoidal rule with 0.66  causes numerical damping of results. Furthermore, the generalized trapezoidal rule with 0.5  leads to numerical oscillations. Numerical oscillations and/or numerical 130 M. MANDIĆ, I. JURIĆ-GRGIĆ, N. GRULOVIĆ-PLAVLJANIĆ damping of results are the consequences of the choice of the numerical integration parameter (ϑ). In principle, its choice is based on numerical experiments. The influence of ϑ on stability of one-step time integration methods for the first order initial value problems is discussed in [9-10], where the stability behaviour is investigated by examining the free response of the FE model. It was found that methods for which 0.5 1  are unconditionally stable. Fig. 4 Generator „1‟ mechanical power (Pm) during angle stability simulation Fig. 5 Generator „2‟ mechanical power (Pm) during angle stability simulation Improving Time Integration Scheme for FET Analysis of Power System Angle Stability 131 4. CONCLUSION In this paper an improved time integration scheme for numerical analysis of power system angle stability using FET was presented. Improvement of accuracy was achieved by using Heun's method for the time integration scheme when forming a FE equations system. The proposed numerical method has been tested and compared to results obtained by ϑ - method and EMTP program. It has been shown that Heun‟s method for time integration gives results with higher accuracy than the generalized trapezoidal rule (ϑ – method) previously used. In some cases where the generalized trapezoidal rule is used, the choice of the numerical integration parameter (ϑ) generates numerical oscillations and/or numerical damping of results. Using the proposed approach, these unwanted consequences are avoided. REFERENCES [1] Z. de Souza, L. Bil, "Unified computational tool for transient and long‐term stability studies", IET Gener Transm Distrib, vol. 3, no.2, pp. 173‐181, 2009. [2] S. Henschel, Analysis of electromagnetic und electromechanical power system transients with dynamic phasors, Ph.D. dissertation, University of British Colombia, 1999. [3] V. Venkatasubramanian, "Tools for dynamic analysis of the general large power system using time- varying phasors", International Journal of Electrical Power and Energy Systems, vol. 16, no. 6, pp. 365- 376, 1994. [4] T. Odun-Ayo, M.L. Crow, "An analysis of power system transient stability using stochastic energy functions", International Transactions on Electrical Energy Systems, vol. 23, no. 2, pp. 151-165, 2013. [5] R. Adapa, J. Reeve, "A new approach to dynamic analysis of AC networks incorporating detailed modeling of DC systems. II. Application to interaction of DC and weak AC systems", IEEE Transactions on Power Delivery, vol. 3, no. 4, pp. 2012-2019, 1988. [6] M.R. Zarate, C.T. Van, M. Federico, J. Conejo Antonio, "Securing transient stability using time-domain simulations within an optimal power flow", IEEE Trans Power System., vol. 25, no. 1, pp. 243‐253, 2010. [7] B.Y. Bagde, B.S. Umre, K.R. Dhenuvakonda, "An efficient transient stability‐constrained optimal power flow using biogeography‐based algorithm", International Transactions on Electrical Energy Systems, vol. 28, no. 1, e2467, 2018. [8] I . Jurić-Grgić, N . Grulović-Plavljanić, M. Dabro, "An Analysis of Power System Transient Stability Using Finite Element Technique", International Transactions on Electrical Energy Systems, vol. 29, no. 1, e2647. 2019 [9] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, McGraw-Hill: Vol 1, London, UK, 1989. [10] O.C. Zienkiewicz, K. Morgan, Finite Elements and Approximation, John Willey & Sons: New York, USA, 1983. [11] S.C. Chapra, R.C. Canale, Numerical Methods for Engineers (7th Edition), McGraw-Hill: New York, 2015. [12] J. Arrillaga, C.P. Arnold, B.J. Harker, Computer Modelling of Electrical Power Systems. John Wiley and Sons: New York, USA, 1983. [13] J. Arrillaga, C.P. Arnold, Computer Analysis of Power Systems. John Wiley and Sons: New York, USA, 1990. [14] M. Dabro, I. Jurić-Grgić, R. Lucić, "EMTP-RV Model of Hydraulic Digital Governor", International Review on Modelling and Simulations, vol. 4, no. 6, pp. 1-5, 2011. [15] IEEE Standard 421.5-2016: Recommended Practice for Excitation System Models for Power System Stability Studies, 2016.