=Paper=
{{Paper
|id=Vol-1839/MIT2016-p38
|storemode=property
|title= Specification and analysis of transients in electrical power systems using the methodology of hybrid systems
|pdfUrl=https://ceur-ws.org/Vol-1839/MIT2016-p38.pdf
|volume=Vol-1839
|authors=Yury Shornikov,Dmitry Dostovalov,Maria Nasyrova
}}
== Specification and analysis of transients in electrical power systems using the methodology of hybrid systems==
https://ceur-ws.org/Vol-1839/MIT2016-p38.pdf
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Specification and Analysis of Transients in
Electrical Power Systems Using the
Methodology of Hybrid Systems
Yury Shornikov1,2 , Dmitry Dostovalov1,2 , and Maria Nasyrova1
1
Novosibirsk State Technical University,
20, Prospekt K. Marksa, 630073 Novosibirsk, Russia
http://nstu.ru
2
Design Technological Institute of Digital Techniques, SB RAS
6, Academica Rzhanova St., 630090 Novosibirsk, Russia
http://kti.nsc.ru
shornikov@inbox.ru,d.dostovalov@corp.nstu.ru,maria_myssak@mail.ru
Abstract. This paper discusses discrete-continuous (hybrid) systems
and corresponding simulation tools. Modern hybrid systems (HS) for-
malism can be effectively used in problem-oriented environments of com-
puter analysis. One of the many HS applications is the study of tran-
sients 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 associ-
ated 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 pro-
cesses in EPS. The solver of ISMA uses the library of classical and original
numerical methods intended for solving systems of differential-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, nu-
merical analysis, library of numerical methods, event detection, principle
circuits.
1 Introduction
Hybrid systems (HS) theory is a modern and versatile apparatus for mathe-
matical 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
445
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
continuous modes [1], [2], [3], [4]. Each mode is given by a set of differential-
algebraic equations with the following constraints:
𝑦 ′ = 𝑓 (𝑥, 𝑦, 𝑡), 𝑥 = 𝜙(𝑥, 𝑦, 𝑡),
𝑝𝑟 : 𝑔(𝑦, 𝑡) < 0, (1)
𝑡 ∈ [𝑡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 tradition-
ally 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 elec-
trical machines. Mode parameters of other elements belong to the synchronously
rotating coordinate system. The equations for the currents and voltages are writ-
ten 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 cir-
cuit 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 ser-
vices and techniques for computational experiments. Leading native (RastrWin,
ANARES) and foreign (EUROSTAG, DIgSILENT PowerFactory, PSS○E) R soft-
ware 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 anal-
ysis problems of transient processes in power systems requires the development
446
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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 or-
ganization 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
Fig. 1. Architecture of ISMA.
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 ki-
netics (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.
447
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
3 Test Model of EPS
A test circuit of institute ”Energosetproject” [9] is given as an example of solv-
ing 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 electrome-
chanical transients is built. Synchronous machines are described by Park-Gorev
equations in normal form [7], [8], [9]. Generator 𝐺1 is a powerful hydroelectric
Fig. 2. Scheme of EPS.
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 syn-
chronous machine 𝐺𝑖 , 𝑖 = 1, . . . , 6 have the differential-algebraic form (2) and
(3).
𝑑𝛹𝑖
= −𝑈𝑖 − 𝜔𝑖 · 𝛾 · 𝛹𝑖 ,
𝑑𝑡
𝑑𝛹𝑓 𝑖 𝑟𝑓 𝑖
= 𝐸𝑞𝑒𝑖 · − 𝑟𝑓 𝑖 · 𝑖𝑓 𝑖 , (2)
𝑑𝑡 𝑥𝑎𝑑𝑖
𝑑𝜔𝑖 1
= · (𝑀𝑇 𝑖 − 𝑀𝑖 ),
𝑑𝑡 𝑇𝐽𝑖
448
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
1 1
𝛹𝑑𝑖 = · 𝑥𝑑𝑖 · 𝑖𝑑𝑖 + · 𝑥𝑎𝑑𝑖 · 𝑖𝑓 𝑖 ,
𝜔𝑛𝑜𝑚 𝜔𝑛𝑜𝑚
1
𝛹𝑞𝑖 = · 𝑥𝑞𝑖 · 𝑖𝑞𝑖 ,
𝜔𝑛𝑜𝑚 (3)
1 1
𝛹𝑓 𝑖 = · 𝑥𝑓 𝑖 · 𝑖𝑓 𝑖 + · 𝑥𝑎𝑑𝑖 · 𝑖𝑑𝑖 ,
𝜔𝑛𝑜𝑚 𝜔𝑛𝑜𝑚
0 = −𝑀𝑖 + 𝑖𝑞𝑖 · 𝛹𝑑𝑖 − 𝑖𝑑𝑖 · 𝛹𝑞𝑖 ,
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
0 1
frequency, 𝜔𝑛𝑜𝑚 is a rated frequency, 𝛾 = , 𝛹𝑓 𝑖 is an interlinkage of
−1 0
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).
𝑑𝛿𝑖1
= 𝜔𝑖 − 𝜔1 , 𝑖 = 2, . . . , 6. (4)
𝑑𝑡
The equation for the block transformer of the generator 𝐺1 has of the form:
1 𝑑𝐼1 𝑈1 𝑈8 𝜔1
· = − − · 𝛾 · 𝐼1 , (5)
𝜔𝑛𝑜𝑚 𝑑𝑡 𝑥𝑚1 𝑥𝑚1 𝜔𝑛𝑜𝑚
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 :
1 𝑑𝐼10 𝑈9 𝑈8 𝑟𝑝8−9 · 𝑈8 𝜔1
· = − − − · 𝛾 · 𝐼10 , (6)
𝜔𝑛𝑜𝑚 𝑑𝑡 𝑥𝑝8−9 𝑥𝑝8−9 𝑥𝑝8−9 𝜔𝑛𝑜𝑚
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).
1 𝑑𝑈8 𝜔1
· = 𝑥𝑐8 · 𝐼20 − · 𝛾 · 𝑈8 . (7)
𝜔𝑛𝑜𝑚 𝑑𝑡 𝜔𝑛𝑜𝑚
449
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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).
1 𝑑𝐼15 𝑈7 𝑟𝑙7 𝜔1
· = − · 𝐼15 − · 𝛾 · 𝐼15 . (8)
𝜔𝑛𝑜𝑚 𝑑𝑡 𝑥𝑙7 𝑥𝑙7 𝜔𝑛𝑜𝑚
The complete mathematical model of the EPS contains equations of excita-
tion systems and speed regulators of generators as well as the current balance
equations at the nodes and coordinate transformation. Coordinate transforma-
tion (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 ) · 𝐼𝑖 ,
(9)
𝑈𝑖(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)
(|𝛼𝑖 | − 𝑍𝐻𝑖 ) · 𝑠𝑖𝑔𝑛(𝛼𝑖 ), 𝑓 𝑜𝑟|𝛼𝑖 | > 𝑍𝐻𝑖 ,
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
450
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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)
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 = 𝑦𝑛 +
ℎ𝑛+1 𝜙𝑛 , 𝑛 = 0, 1, 2, . . . . As a result we obtain the dependence of the predicted
integration step ℎ𝑛+1 .
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
∑︁
𝑦𝑛+1 = 𝑦𝑛 + 𝑝𝑖 𝑘𝑖 , 𝑘𝑖 = ℎ𝑛 𝑓 (𝑦𝑛 + 𝛽𝑖𝑗 𝑘𝑗 ),
𝑖=1 𝑗=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 im-
plementation 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
451
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Fig. 3. Scheme of a power supply of neighborhood.
Table 1. Library of Numerical Methods
Method (𝑝, 𝑚) Description
DISPF (5, 6) Stability control, systems of medium and low stiffness
RADAU5 (3, 3) Stiff systems
DISPF1 RADAU Adaptive method DISPF in combination with RADAU5 with
stiffness control, essentially stiff systems
DP78ST (8, 13) Stability control, variable order and step, systems of medium
stiffness and high precision
RKF78ST (7, 13) Stability control, variable order and step, systems of medium
stiffness and high precision, based on Runge-Kutta-Feldberg
method[17]
RK2ST (2, 2), Explicit methods with stability control for analysis of non-stiff
RK3ST (2, 3) systems
DISPS1 Algorithm of variable order with adaptive stability region
MK22 (2, 2), Freezing of Jacobean matrix, stiff systems
MK21 (2, 2)
MK11F Algorithm of analysis of implicit problems
452
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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 Runge-
Kutta 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 Runge-
Kutta 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 al-
ternating 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 inde-
pendent 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 numeri-
cal methods. API is a public interface of the computing module consisted of
public classes and interfaces used by other components to interact with the im-
plemented 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)
𝜕𝑦 𝜕𝑡
453
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Theorem. The choice of the step according to the formula
𝑔𝑛
ℎ𝑛+1 = (𝛾 − 1) , 𝛾 ∈ (0, 1), (13)
𝜕𝑔𝑛 𝜕𝑔𝑛
· 𝜙𝑛 +
𝜕𝑦 𝜕𝑡
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, . . . . Con-
verting recurrently this expression we get 𝑔𝑛+1 = 𝛾 𝑛+1 𝑔0 . Given that 𝛾 < 1,
then 𝑔𝑛 → 0 takes place when 𝑛 → ∞. In addition, condition 𝛾 > 0 implies that
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
area 𝑔(𝑦𝑛 , 𝑡𝑛 ) = 0, which completes the proof.
6.1 Control of Event Function in the Integration Algorithm
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.
454
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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.
Fig. 4. The moment of generator 𝐺4 turbine.
The results of six-machine EPS simulation correspond with those in the orig-
inal 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 sys-
tem 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 con-
tributes 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 sci-
entific research area for many years. Modern formalisms and methods for analysis
of complex systems can be effectively used by domain specialists only if domain-
specific 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 prepa-
ration for computation including choice of the effective step size and numerical
method, start and control of the computational experiment process. Expanding
455
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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 com-
ponents 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 fre-
quently 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 special-
ist 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-DE-
EPPKA2-CBHE-JP, Innovative teaching and learning strategies in open mod-
elling 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 Stiff Hybrid Systems:
Monograph. Publishing House of Novosibirsk State Technical University, Novosi-
birsk (2012)
456
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
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 Ex-
isting 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” Soft-
ware. Certificate of State Registration of Computer Programs # 2015617235. Fed-
eral 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 Stiff Systems. Nauka, Novosibirsk (1997)
17. Fehlberg E.: Classical fifth-, sixth-, seventh- and eighth order Runge-Kutta formu-
las 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 differential
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 stiff hybrid systems with one-
sided 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 differential equations II. Stiff and
differential-algebraic problems. 1996, Berlin: Springer.
22. Shornikov, Yu.V., Novikov, A.E., Novikov, E.A. Approximation of the Jacobi ma-
trix in (2,2)-method for solving stiff 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
differential equations. Computer, 1963, 5, 329330.
457