Specification and Analysis of Transients in Electrical Power Systems Using the Methodology of Hybrid Systems

Abstract

This paper discusses discrete-continuous (hybrid) systems and corresponding simulation tools. Modern hybrid systems (HS) formalism can be efectively used in problem-oriented environments of computer analysis. One of the many HS applications is the study of transients in electrical power systems (EPS). A special module is developed in ISMA (translated from Russian ”Instrumental Facilities of Machine Analysis”) simulation environment to support the research of transient processes by Park-Gorev's equations. Discrete behavior of EPS associated with nonlinear characteristics of generator speed regulators. Also, the EPS operating mode can be changed upon the occurrence of certain events: switching, short circuit, breakage of power lines, etc. Therefore HS methodology is adequate for description and study of transient processes in EPS. The solver of ISMA uses the library of classical and original numerical methods intended for solving systems of diferential-algebraic equations with discontinuities. The original algorithm of correct event detection is developed for processing gaps, which is an integral part of numerical analysis.

Keywords

• Hybrid systems
• modal behavior
• transient processes
• numerical analysis
• library of numerical methods
• event detection
• principle circuits

Introduction

Hybrid systems (HS) theory is a modern and versatile apparatus for mathematical description of the complex dynamic processes in systems with different physical nature (mechanical, electrical, chemical, biological, etc.). Such systems are characterized by points of discontinuity in the first derivative of the phase variables. HS behavior can be conveniently described as sequential changes of continuous modes [ 1 ], [ 2 ], [ 3 ], [ 4 ]. Each mode is given by a set of differentialalgebraic equations with the following constraints: (1)

), : (, ) < 0, ∈ [ 0, ], ( 0) = 0, ( 0) = 0, where ∈ , ∈ , ∈ , : × × → , : × → . A scalar function (, ) is called event function or guard [ 1 ], [ 4 ], [ 5 ], [ 6 ]. The inequality (, ) < 0 means that the phase trajectory in the current state should not cross the border (, ) = 0. Therefore HS state is determined by a predicate . The system is in the appropriate state when = . Events occurring in violation of this condition and leading to a transition to a different state without crossing the border, called the one-sided [ 4 ], [ 6 ]. Such events are of particular practical interest. For example, due to changes in the electrical power systems (EPS) configuration, the operating mode cannot be determined at the time the event occurred.

Models of EPS based on the Park-Gorev equations [ 7 ], [ 8 ], [ 9 ] are traditionally used to describe the electromagnetic and electromechanical processes when studying synchronous operation of generators and solving many other problems in the analysis of electric power systems. Generator mode parameters are written in the rotating coordinate system − associated with respective rotors of electrical machines. Mode parameters of other elements belong to the synchronously rotating coordinate system. The equations for the currents and voltages are written by the laws of Ohm’s and Kirchhoff’s for each of the axes of the coordinate system in accordance with the circuit topology. Thus, equations of electric circuit and its elements completely match the class (1). Therefore, the developed mathematical and instrumental tools for HS simulation can be unified to the electricity problems.

New formalism and methodology for the analysis of complex dynamic systems should be implemented in a problem-oriented environment with plenty of services and techniques for computational experiments. Leading native (RastrWin, ANARES) and foreign (EUROSTAG, DIgSILENT PowerFactory, PSS ○RE) software systems for the calculation of steady-state and transient processes in EPS implement the traditional models and methods of analysis [ 9 ], [ 10 ], [ 11 ]. Experts in the power industry almost never use modern methodology of hybrid systems. Therefore, the task of developing custom tools with object-oriented interface and input language, a new formalism and original interpretation mechanisms is new and topical. 2

Specification Languages in ISMA

Software for instrumental analysis of HS in ISMA[ 4 ], [ 12 ], [ 13 ] is unified to the problems of different nature: the study of simple dynamic processes, automatic control, chemical kinetics, electrical engineering. Unification to computer analysis problems of transient processes in power systems requires the development of tools for specification of EPS program models. Fig. 1 shows the architecture of instrumental environment ISMA. Designed architecture allows customizing the environment to a new application with minimal modifications in the organization of interaction of available modules and libraries with object-oriented graphics editor and model interpreter. Modules discussed in this paper are grayed out. The tools provide five different input languages for computer analysis using the methodology of hybrid systems: multi-purpose symbolic (LISMA PDE)[ 14 ] and graphic (Harel charts) languages, as well as thematic language of block diagrams of automatic control systems (ACS) and the equations of chemical kinetics (LISMA+). The language of principal EPS circuits (LISMA EPS) is also thematic. In the proposed architecture, subject-oriented interfaces interact with the computing core of the system through a universal internal representation of hybrid models. This ensures the continuity of the developed software for new applications with its own characteristics. 3

Test Model of EPS

A test circuit of institute ”Energosetproject” [ 9 ] is given as an example of solving the electricity problems. Schematic diagram of the closed energy system with two voltage levels and six synchronous machines of different types and powers is presented in Fig. 2. A mathematical model for the calculation of electromechanical transients is built. Synchronous machines are described by Park-Gorev equations in normal form [ 7 ], [ 8 ], [ 9 ]. Generator 1 is a powerful hydroelectric plant. Generators 2 and 3 simulate a thermal power plant with a small power which aggregates operate on power line with different rated voltage 500 and 220 kV. Generators 4 and 5 simulate a thermal power plant with a high power, which aggregates also operate on power lines of different voltage. 6 helps to simulate synchronous compensators mounted to the hub substation.

Here are the equations of main network elements. The equations of synchronous machine , = 1, . . . , 6 have the differential-algebraic form (2) and (3).

= − − · · , = = 1 ·

− · , · ( − ), (2) = = = 1 · 8 = 8 · 20 − 1 · · 8. where = ( , ) are projections of interlinkage of the stator windings on the axis and , = ( , ) are projections of the stator windings voltages, = ( , )

are projections of the stator current, is rotor rotation frequency, is a rated frequency, = , is an interlinkage of excitation winding, is a electromotive force, and are synchronous inductances of longitudinal and transverse axis, is inductive reactance of the excitation winding, is resistance of the field winding, and are inductances of longitudinal and transverse stator reaction, is a field winding current, is inertia constant,

is turbine moment, is electromagnetic torque. The mutual angles of the generator rotors are defined relative to the ︂( generator 1 by the formula (4).

The equation for the block transformer of the generator 1 has of the form: 1 = − 1, = 2, . . . , 6. where 1 is inductive resistance of the transformer.

The equations for generators 2–6 are written likewise considering numbering and using voltages (1) and currents (1), = 2, . . . , 6, given to coordinate system of the base machine 1. The equations of the autotransformers 1– 3 are similar to (5).

The equations for the longitudinal elements of the power line 3: where 8−9 are inductive resistance of a branch, 8−9 is active resistance.

Inductive conductivities of power line are equated to conductivities on the start and on the end node. Thus, the node 8 gets the equations (7). For the other power lines program models are written in form (6) and (7). (3) (4) (5) (7)

Loads are given by active-inductive elements and reactor is given only by inductance. For example, the load at the node 7 is given by equation (8).

The complete mathematical model of the EPS contains equations of excitation systems and speed regulators of generators as well as the current balance equations at the nodes and coordinate transformation. Coordinate transformation (9) is necessary for communication of machines currents and voltages at the nodes connecting generators to the coordinate system of the base machine 1. (1) = (1 · cos 1 + · sin 1 ) · , (1) = (1 · cos 1 + · sin 1 ) · , = 2, . . . , 6.

Synchronous machines 1, 4, 5 and 6 has installed automatic excitation regulators (AER) of strong action. Generators 2 and 3 are equipped with AER of proportional action. As models of turbine regulators applied the ones used in the software package MUSTANG [ 9 ]. Speed controller model is described by no more than two differential equations. Speed regulators of turbine of generators 2– 5 have a deadband = ︃{ 0,

| | < , (| | − ) · ( ), | | > , = 2, . . . , 6.

(10) (8) (9)

Here is signal to input of the speed control system (in relative units), is deadband, is movement of the clutch of the centrifugal pendulum. All values are given in relative units.

General mathematical model of the analyzed EPS contains 279 nonlinear differential-algebraic equations. The hybrid behavior of the system under study is due to the equations (10). 4

Graphical Specification

The program model of the test circuit [ 9 ] contains 173 lines in LISMA PDE. This form of presentation is useful when checking the correctness of the mathematical notation and experimenting with various mathematical models of components and the network as a whole. However, for the specialist the great importance has the ability to quickly change the scheme, adding and removing elements and connections, edit the properties of the elements. Editing a text model in this case becomes very time-consuming and increases the likelihood of errors. In such cases, the representation of the problem as a circuit diagram of the EPS is preferred. To do this a graphical editor of principal schemes is developed in ISMA. The interface of the editor is shown in Figures 2 and 3. Fig. 3 shows a power supply circuit of the neighborhood in Novosibirsk [ 15 ]. Presentation of computer models in the form of concepts of EPS is more compact and familiar to user. Tools converting graphical models to universal internal representation of HS hide from the user the process of constructing a mathematical model and allow to quickly getting the results of a computational experiment in the treated and convenient way. 5

Library of Numerical Methods in ISMA

This section is devoted to the integration algorithms of variable order and step based on explicit methods of Runge-Kutta type. The algorithms are applied to numerically solve the Cauchy problem for ODE systems of the following form: ′ = ( ), ( 0) = 0, 0 ≤ ≤ . (11) ℎ +1 , integration step ℎ +1.

Consideration of autonomous problem does not reduce the generality because non-autonomous problem always can be cast to autonomous by introducing an additional variable. Particular attention should be paid to the choice of the integration method. Fully implicit methods cannot be used because they require the calculation of ( ) at a potentially dangerous area, where the model is not defined. Therefore here we will use explicit methods with solution: +1 = + = 0, 1, 2, . . . . As a result we obtain the dependence of the predicted Considering that explicit methods are known by poor stability this paper examines integration methods with accuracy and stability control. Generally accuracy and stability control are used to limit the size of the integration step. As a result projected step ℎ +1 is calculated as follows.

The choice of the next integration step size is based on the proved theorem [ 16 ] and can be written as follows:

ℎ +1 = max[ℎ , min(ℎ , ℎ )], where ℎ and ℎ

are step sizes obtained as a result of accuracy control and stability control respectively. This formula allows to stabilize the step behavior in the area of solution establishing where stability plays a decisive role. Because the presence of this area severely limits the use of explicit methods for solving stiff problems. Suppose that for numerical solution of problem (11) the following implicit methods of Runge-Kutta type is used: +1 = + , = ℎ ( +

), ︁∑ where and are real -dimensional vector-functions, ℎ is an integration step, are the method stages, and

are numerical coefficients.

Peculiarities of numerical analysis are defined by the configuration and implementation of the solver in the scheme interpreter. Solver is configured to numerical analysis not only of smooth dynamical systems but also systems with ordinary discontinuity and stiff systems [ 4 ]. For the analysis of the stiff modes new original -phasic methods of -order (Table 1), developed by the authors, are included in the solver library.

DISPF(5, 6), DP78ST(8, 13), RKF78ST(7, 13), RK2ST(2, 2), RK3ST(2, 3), and DISPS1 are integration algorithms based on explicit methods of the RungeKutta type [ 4 ][ 16 ][ 20 ]. All of these algorithms have accuracy and stability control and are applied to solving multi-mode problems. RADAU5(3, 3), MK22(2, 2), MK21 (2, 2) are based on fully-implicit schemes which are aimed at solving single-mode systems. The first one is the implicit fifth order three stage RungeKutta method, the second and third schemes belong the (m, k)-class of numerical methods [ 4 ][ 21 ][ 22 ]. DISPF1 RADAU is the variable structure algorithm of alternating order and step based on explicit Runge-Kutta numerical formulas of the first, second, and fifth order and the implicit Runge-Kutta type method of the fifth order. This algorithm has accuracy and stability control and can be applied to single-mode as well as multi-mode systems [ 4 ]. Finally, MK11F is the algorithm based on the L-stable Rosenbrock method aimed at solving implicit problems [ 19 ][ 23 ].

Libraries of standard blocks and numerical methods are implemented as independent application modules that are loaded at run time. This approach allows to allocate in the application programming interface (API) a set of functions and classes required for the implementation of element libraries and numerical methods. API is a public interface of the computing module consisted of public classes and interfaces used by other components to interact with the implemented solvers and to create new ones. API classes describes the subject area and declares the type of systems and problems recognized by the solver. Using the API any user with basic knowledge of object-oriented programming able to develop and built in the system new typical elements and numerical methods without recompiling the entire system. 6

Event Detection in Hybrid Systems

The correct analysis of hybrid models is significantly depends on the accuracy of detection [ 6 ][ 18 ] of the change of the local states of the HS. Therefore, the numerical analysis is necessary to control not only the accuracy and stability of the calculation, but also the dynamics of the event-function. The degree of approximation by the time the event occurred is defined by the behavior of event driven function.

Analyze the behavior of the event function (, ). Let the method of the form +1 = + ℎ , where function is calculated in point , is used for calculations. Then the event-function (, ) at point +1 has a form +1 = ( + ℎ , + ℎ ). Decomposing the +1 in a Taylor series and taking into account the linearity of +1, we obtain the dependence of +1 of the projected step ℎ : +1 = + ℎ ( · + ).

(12) Theorem. The choice of the step according to the formula

· + provides the event-dynamics behavior as a stable linear system, the solution of which is asymptotically approaching to the surface (, ) = 0.

Proof. Substituting (13) in (12), we have +1 = , = 0, 1, 2, . . . . Converting recurrently this expression we get +1 = +1 0. Given that < 1, then → 0 takes place when

→ ∞. In addition, condition > function does not change sign. Therefore, when 0 < 0, < 0 will be valid for all . Then the guard condition will never cross the potentially dangerous 0 implies that area ( , ) = 0, which completes the proof. We complete the implicit problem’s integration algorithm by the algorithm of the step control that takes into account the event function dynamics. Let the solution at the point is calculated with the step ℎ . In addition, the new accuracy step ℎ +1 is computed by the formula (13). Then the approximate solution at the point +1 is calculated as follows:

Step 1. Calculate the functions

Step 2. Calculate where = .

= ( , ), = ( , ) , = ( , ) .

′ = · +

, Step 3. If ′ < 0, then ℎ +1 = ℎ +1 and go to the Step 6.

Step 4. Calculate the new ”Event” step ℎ +1 by the formula ℎ +1 = ( − 1) .

′ Step 5. Calculate the new step ℎ +1 by the formula ℎ +1 = min(ℎ +1, ℎ +1). Step 6. Go to the next integration step.

In the Step 3, unlike the previously presented algorithm [ 4 ], we determine the direction of event-function change. Near the boundary regime denominator (13) will be positive, and away from the boundary (, ) = 0 it becomes negative. Then, defining the direction of event-function change, we do not impose any further restrictions on the integration step if the event-function is removed from the state boundary. 7

Simulation Results

At time t=0 s a far electrical load decreases by 10% initiating an electromagnetic transient. Numerical experiment results are shown in Fig. 4. The calculations are performed by RK3ST algorithm with initial step 0.000001 s. Results from [ 9 ] obtained by implicit Euler method with integration step 0.00001 s are shown in the same figure.

The results of six-machine EPS simulation correspond with those in the original source [ 9 ] and do not contradict the theoretical concepts that confirms the correctness of the method used. Discrepancy between the results is based not only on different numerical methods used by [ 9 ] and the authors. In [ 9 ] the system is considered as a continuous one, so the problem of correct hybrid system event detection arisen from regulators’ deadzone is not dealt with [ 10 ]. It contributes to the accumulation of a global calculation error. Using a correct event detection algorithm let obtain better simulation results, as shown by the authors in [ 20 ]. It should be noted that the sensitivity to chosen formalisms (continuous, hybrid) varies from one phase variable to another. 8

Conclusion

Computer-aided analysis of discrete-continuous systems has been an actual scientific research area for many years. Modern formalisms and methods for analysis of complex systems can be effectively used by domain specialists only if domainspecific software tools have been developed. Such software frees end-users from routine translation from a mathematical model to its program implementation and helps to carry out a computer experiment. This software solves problems of the optimal representation of a model in the computer memory, a model preparation for computation including choice of the effective step size and numerical method, start and control of the computational experiment process. Expanding a modeling and simulation environment must demand of minimal changes of existing components as well as development of model specification tools.

For a description and analysis of transients in power systems and their components the use of the methodology of hybrid systems surrounded by the tools of computer analysis is proposed. Approaches to specification of EPS in ISMA instrumental environment are presented. Textual specification is fully consistent with the mathematical notation and can be used for experiments with different models of network elements. Graphical specification convenient for specialists in the field of electricity when the network configuration of the circuit is frequently changing. It should be noted that ISMA has tools for translation from graphical to textual model. Thus, it is possible to verify the program models. The main advantage is the dramatic increase in the efficiency of research when a specialist using an instrumental service deals exclusively with the analysis of the results of their design decisions. While in the traditional practice energy specialist to get the results requires additional expertise in the field of computational mathematics and programming.

The architecture of the solver for the continuous behavior of hybrid systems is proposed. The library of numerical methods can be easily extended by new methods of the Runge-Kutta type as well as other one-step integration methods for ODE systems. The API also provides the mechanism to add implementations of algorithms that deal with another type of systems. The new original method of switching point’s localization is proposed. The algorithm easily complements the existing numerical solvers based on explicit and semi-explicit schemes.

The presented results of the test problem (six-machine EPS) calculation are obtained using the considered approaches and methods. Thus, the correctness of theoretical assumptions, mathematical and algorithmic software is constructively proved.

Acknowledgments. This work was supported by the grant of the Russian Foundation for Basic Research (RFBR grant 17-07-01513) and by the grant of the Education, Audiovisual and Culture Agence (EU), programme ERASMUS + Capacity building in higher education, project 573751-EPP-1-2016-1-DEEPPKA2-CBHE-JP, Innovative teaching and learning strategies in open modelling and simulation environment for student-centered engineering education.

References

1. Harel , D. , Pnueli , A. , Schmidt , J.P. , Sherman , R.: On the Formal Semantics of Statecharts . In: Proc. 2nd IEEE Symp. On Logic in Computer Science , pp. 54 - 64 . IEEE Press, New York ( 1987 )

2. Maler , O. , Manna , Z. , Pnueli , A. : From Timed to Hybrid systems . In: Real-Time: Theory in Practice, LNCS , vol. 600 , pp. 447 - 484 . Springer-Verlag ( 1992 )

3. Senichenkov , Yu .B.: Numerical Simulation of Hybrid Systems . Publishing House of Politechnical University, St. Petersburg ( 2004 )

4. Novikov , E.A. , Shornikov , Yu.V.: Computer Simulation of Stif Hybrid Systems: Monograph. Publishing House of Novosibirsk State Technical University, Novosibirsk ( 2012 )

5. Lee , E.A. , Zhenq , H.: Operational Semantics of Hybrid Systems . In: Proc. of Hybrid Systems: Computational and Control (HSCC) , vol. 3414 , pp. 25 - 53 . Zurich ( 2005 )

6. Esposito , J.M. , Kumar , V.: A State Event Detection Algorithm for Numerically Simulating Hybrid Systems with Model Singularities . In: ACM Transactions on Modeling and Computer Simulation , vol. 17 , issue 1, pp. 1 - 22 ( 2007 )

7. Venikov , V.A. : Transient Electromechanical Processes in Electrical Systems . Vysshaja Shkola, Moscow ( 1985 )

8. Lukashov , Je.S., Kalyuzhnyj , A.H. , Lizalek , N.N. : Long-Term Transients in Power Systems . Nauka, Novosibirsk ( 1985 )

9. Fomina , T. Yu .: Development of the Algorithm for Calculating the Transients in Complex Regulated EPS: Dissertation of the Candidate of Technical Sciences (PhD) . Moscow ( 2014 )

10. System Research in the Energy Sector . In: Proceedings of Young Scientists ESI SB RAS , vol. 44 . ESI SB RAS , Irkutsk ( 2014 )

11. Nepsha , F.S. , Otdel'nova, G.V. , Savinkina , O.A. : Compare the Functionality of Existing Software to Calculate and Analyze of Electric Modes . In: Vestnik of Kuzbass State Technical University, vol. 2 ( 96 ), pp. 116 - 118 ( 2013 )

12. Shornikov , Yu .V., Bessonov , A.V. : The Core Components of ”ISMA 2015” Software . Certificate of State Registration of Computer Programs # 2015617235 . Federal Service For Intellectual Property (Rospatent), Moscow ( 2015 )

13. Shornikov , Yu .V., Bessonov , A.V. : A Unified Approach to Computer Simulation of Hybrid Systems . In: Information Technology of Modeling and Control , vol. 3 ( 93 ), pp. 289 - 298 ( 2015 )

14. Shornikov , Yu .V., Bessonov , A.V. : The Component for Specification of Hybrid Systems in LISMA PDE Language . Certificate of State Registration of Computer Programs # 2015617191 . Federal Service For Intellectual Property (Rospatent), Moscow ( 2015 )

15. Skurihina , K.A. , Arestova , A.Yu. , Armeev , D.V. : A Study of Dynamic Properties of MicroGrid in Parallel Operation with the Power System . In: Siberian Journal of Science , vol. 1 ( 15 ), pp. 93 - 102 ( 2015 )

16. Novikov , E.A. : Explicit Methods for Stif Systems . Nauka, Novosibirsk ( 1997 )

17. Fehlberg E. : Classical fifth-, sixth-, seventh- and eighth order Runge-Kutta formulas with step size contro . In: Computing , vol. 4 , pp. 93 - 100 ( 1969 )

18. Esposito , J.M. , Kumar , V. , Pappas , G.J.: Accurate event detection for simulating hybrid systems . In: Hybrid Systems: Computation and Control (HSCC) , Volume LNCS 2034 , pp. 204 - 217 ( 2001 )

19. Novikov , E.A. , Yumatova , L.A. Some methods for solving of ordinary diferential equations unresolved with respect to derivative . Reports of the USSR Academy of Sciences , 1987 , 295 ( 4 ), 809 - 812 .

20. Dostovalov , D.N. , Shornikov , Yu.V. Simulation of stif hybrid systems with onesided events in the ISMA environment . Computer simulation 2012: Proceedings of the international workshop, St. Petersburg , 20 - 21 September 2012 , 2012 , 34 - 41 .

21. Hairer , E. , Wanner , G. Solving ordinary diferential equations II. Stif and diferential-algebraic problems . 1996 , Berlin: Springer.

22. Shornikov , Yu .V., Novikov , A.E. , Novikov , E.A. Approximation of the Jacobi matrix in (2,2)-method for solving stif systems . Proceedings of the Russian Higher School of Academy of Sciences of Russian Federation , 2008 , 1 ( 10 ), 31 - 44 .

23. Rosenbrock , H.H. Some general implicit processes for the numerical solution of diferential equations . Computer , 1963 , 5 , 329330 .