Method of average voltages in integration step: theory and application

A numerical one-step method of average voltage values in integration step, which allows to algebraize the differential equations for electric circuits of electromechanical systems with semiconductor converters, was described in the paper. The method is proposed for the creation of digital power circuit real-time models of electrical systems in hybrid models that combine digital and physical components. Digital models in hybrid systems must continuously function over a long time period in conjunction with physical objects. The indicated method provides stability of mathematical models, and is also effective in the sense of the amount of calculations. The article describes the application of the method of average voltage values in the integration step to create a hybrid model of power generation system for a ship's power plant with diesel generator. The digital component of the hybrid model, in this case, is power conversion system that includes a synchronous generator driven by a diesel engine, while the physical part is real excitation controller of generator with automatic voltage regulator. The research results obtained on the hybrid model and their comparison with the results of the physical experiment at the laboratory plant are presented.


Introduction
Systems containing electric circuits are described by differential and algebraic equations, which are usually nonlinear.Their solution requires the application of numerical methods that carry out the algebraization of differential equations in integration step.The algebraization of equations for electric circuits by explicit methods is obviously complicated by the error of calculating currents, as a result of which, after multiple calculations, there is a deviation from the first Kirchhoff's law.The algebraization of these equations by implicit methods is complicated by the error of calculating voltages, as a result of which, after multiple calculations, there is a deviation from the second Kirchhoff's law.These deviations lead to instability of mathematical models, especially those that are used in conjunction with real objects, and in which the integration step is very small, and the calculation time is very long.Authors faced this problem, during the development of hybrid models.Such problems are mentioned in [1][2][3][4][5][6][7][8][9][10].In particular, it is noted that most of the electrotechnical objects are stiff systems [1], for which the use of explicit methods is not feasible, since it may cause numerical instability [2,3,6,7].Reducing the step of numerical integration for stiff systems, if explicit methods are used, slows down the calculation and complicates the real-time implementation of models that interact with real objects [4].Expediency to use implicit methods for stiff systems is mentioned in [5][6][7][8][9][10][11][12], in particular the trapezoidal method, which is A-stable.At the same time, it is noted in [5] for power systems, that in the case of numerical integration step, increasing the trapezoid method gives a worse result in comparison with the explicit (forward) Euler method.Also, the numerical oscillations of frequency at 1/(2Δt) (Δt-the numerical integration step), which are relevant for trapezoid method in case of stiff systems and which in some cases but not always can be eliminated by decreasing the step value, are noted in [1].The disadvantage of using the trapezoidal method that appears in numerical oscillations in case of high-frequency current changes in circuits with inductance or their breaks (for example, in case of triggering switches) is indicated in [6,7].In addition, implicit methods are known to require additional computations at the integration step, which slows down the calculation.To increase the accuracy of numerical integration, it is required to decrease the step.However, long-term calculations with a small step lead to accumulation of roundoff errors [8] and significant increase of calculation time.To eliminate this contradiction, methods of integrating with a variable step or procedures for changing the integration method during the calculation are proposed in [10,11].However, we note that the use of such methods for hybrid models is complicated, since the step in a real-time model is determined by performance of the processor.In addition, in case of the analysis of electric circuits with semiconductor gates, the step change is necessary to find the moment of gate turning-off when the gate current passes through zero.
The peculiarities of using the numerical integration methods for hybrid systems are discussed in [13,14].In particular, it is noted the need for numerical integration with the increased step, due to the limited processor performance, as a result of which explicit methods give significant errors.Also, it was noted that using of methods which are absolutely stable for numerical systems was not always stable in case of hybrid systems that contain an analog component (stability in this case depends on the step of numerical integration).According to the analysis carried out in [13], the advantage of using trapezoidal methods and backward differentiation formulas for hybrid systems was claimed, which gives a better result than implicit Euler method, as an example.
These problems are overcame by the authors using the method of average voltage values in the integration step, which ensures the stability of mathematical models in the indicated sense, and is also effective in terms of the amount of calculations.This method appeared as a result of the need to avoid problems in hybrid models of power generation units in which the generator excitation system is a real object, and the power scheme "turbogenerator-block transformer-power line" is a mathematical model.The basis of this method is first described in [15].The method allows algebraizing differential equations of an electromechanical system with semiconductor converters in the integration step.
The authors for a long time participated in the projects concerned with simulation of electromechanical systems in the field of electric power and transport [16] in the former USSR, Ukraine, Russia, and Poland, in which they used onestep methods of the Runge-Kutta, Adams, and trapezes.To calculate the transients in these projects, the use of mentioned methods was effective.In the real-time models, their operation was limited to 15-20 min, as a result of numerical disturbance related to the algebraization of the equations by the mentioned methods.Thus, using the explicit method, errors arise in the equations that correspond to the first Kirchhoff law, and in case of implicit method-in equations which correspond to the second Kirchhoff law.The problem was solved using the method of average voltage values in the integration step of the second order, which has the features of explicit and implicit methods and simultaneously eliminates mentioned drawbacks of explicit and implicit methods which take place in real-time models of electromechanical systems electric circuits.
Note, that the trapezoidal method, which is widely used in mathematical models, is equivalent in sense of accuracy to the first order method of medium voltage in the integration step, assuming that electrical circuits not include L, C elements connected in series and constraints are linear.In other cases, trapezoidal method is less accurate.

Genesis of the method of average voltage values in integration step for the algebraization of the equations for an electric circle
The equations for a branch of an electric circuit containing time dependencies of: electromotive force e(t) source, voltage on resistance Ri(t), induced voltage d dt , voltage on capacitance u C (t), and to which voltage u(t) is applied as a result of the influence of the electric circuit for this branch, are as follows: We accept that in the integration step, Δt, R and C are constant.Flux linkages ψ is a function of the current in a branch and of other factors taken into consideration by the vector of variables ⃗ q.By integration of the first Eq.( 1), the following equation is obtained: where U, U R , U C , U L , and E are the average values of voltages and electromotive force in the integration step Δt: (1) Representing the integration for the voltages U R and U C by a Taylor's series, the following equations are obtained: Equation (3) after the transformation is written as: where m is number of members of a Taylor series taken into consideration, which corresponds to the order of the method of average voltages in the integration step.Thus, when m = 1-is the first order method, when m = 2-is the second-order method, etc.When m = 1, the current of the branch in a step is a straight line, and the voltage on the capacitance is a parabola.When m = 2, the current of the branch in a step is a parabola, and the voltage on the capacitance is described by a third-order polynomial.

The equation of an electric circle in the method of average voltage values in integration step
The first Eq.( 4) corresponds to the one shown in Fig. 1 which is the equivalent circuit, for the analysis of which the methods of the DC circuit theory can be applied.This scheme allows to determine the current of the branch at the end of the integration step.This branch is described by the following equations: (3) where R S is the step resistance of the branch: E S -is the step electromotive force (e.m.f.), which is determined by the initial conditions for the integration step: i 1 , ψ 1 , and ⃗ q 1 are the variable values at the end of integra- tion step.
The electrical circuit is described by Eq. ( 5) for each of the branches, as well as by the algebraic equations based on the first Kirchhoff's law for currents in independent nodes at the end of the integration step.These equations form a system of electric equilibrium equations.
In case of the presence of electrical machines in the system, the equations system is supplemented by magnetic equilibrium equations, as well as by mechanic equilibrium equations, which allow to calculate the flux linkages of the branches based on the currents of the branches and their mutual placement.
In general, the system of equations is nonlinear and is solved by the Newton method.The algorithm for solving by the Newton method will be explained in a simplified example for one branch in which U and ⃗ q 1 are known.To do this, we shall designate: The increments of current and flux linkages on each iteration are determined as: ,where i 0 1 , 0 1 are zeroth approximations for the current and flux linkage: = −L  -is the dynamic inductance of the electric branch as a function of current i 1 , Therefore: The values of current and flux linkage at the first approximation are written as: -unit of measurement is [Ω]; Δt ,

Fig. 1 E quivalent circuit of the branch
The next approximation is implemented according to the algorithm described for the previous approximation.If the flux linkage is a linear function of current, the algorithm is reduced to a single zeroth iteration.

Description of the scheme
The described method was repeatedly used by the authors to create hybrid models that combine physical and computer (digital) parts.One example was the model of power generation system of the ship's power plant with a diesel generator.
The digital component of the hybrid model (Fig. 2) is power conversion system that contains the synchronous generator G, driven by a diesel engine D with a speed control system, as well as an autonomous load consisting of the synchronous machine G2 and RL branches.(7)

The calculation of voltages and derivatives of currents
The calculation scheme of the power conversation system implemented by the digital model is shown in Fig. 3.This scheme is used for calculation of voltages and current derivative.The scheme consists of six branches that simulate the phases A, B, and C of the stator winding, the excitation winding, and the damping system (the damping system is modeled by two short-circuited RL branches on the axes d and q) of synchronous generator.The inductance and resistance of load are taken into consideration in the inductances and resistances of the stator windings.In addition, the branches of stator windings contain capacitors which take into consideration the load capacity.
The authors use the second-order method of averages voltages in integration step (m = 2).This method requires the calculation of currents derivative in the branches.
The presented scheme is described by the following vector equation according to the second Kirchhoff's law: where  The relation between the voltages on the branches and the potentials of the independent nodes v is described by the following expression: where = 1 1 1 0 0 0 -is the incidence matrix.
Equations for the derivatives of currents according to Kirchhoff's first law is as follows: From Eqs. ( 8) and ( 9), the derivatives of currents are obtained: After substitution ( 11) into (10), the following equation is obtained: The potential of the independent node v is determined from Eq. ( 12), and then, the voltage on the branches is determined from Eq. ( 9) and the derivatives of the currents are determined from Eq. ( 11).

The calculation of the currents by the method of averages voltages in integration step
The calculation scheme of the power conversion system, taking into consideration the statements of the method of average voltage in integration step, is shown in Fig. 4.
If given (5), the vector equation for this scheme is as follows: where the matrix of step resistance: in which R A , R B , and R C are the total resistances of branches consisting of the stator winding and generator load; R f is the resistance of the excitation winding; R D and R Q are the resistances of branches represented the synchronous machine damping system; vector of the step e.m.f: the vector of flux linkages (value at the end of the integration step): the vector of flux linkages (value at the beginning of the integration step): the vector of currents (value at the end of the integration step): the vector of the average at the integration step values of voltage: Fig. 3 Calculation scheme of the power conversation system for the calculation of voltages and derivatives of currents Equation ( 13) is supplemented by the equation for the currents according to Kirchhoff's first law: where = 1 1 1 0 0 0 -is the incident matrix.
The flux linkages in Eq. ( 13) are a function of the currents and the rotation angle γ: The relation between the average in the integration step values of the branch voltages and the potentials of independent nodes is described by the equations: The system of Eqs. ( 13), ( 16)-( 18) is a mathematical model of the presented system.This system of equations is nonlinear and Newton's method is used to solve it.
Denote for initial approximation: where: The vector equation for sequential approximations is as follows: where The increments of currents and average on the step value of voltage (potential of an independent node) on each iteration are determined as: V = 0. Thus, Eq. ( 22) is written as: From the system of Eq. ( 23), the following equation is obtained: where The algorithm of mathematical modeling is cyclic and consists of the following actions: Fig. 4 Calculation scheme of the power conversation system for the current calculation Set initial conditions (values of variables at the beginning of the step).Define the parameters of the scheme in the step: S , ⃗ E S .Set the zeroth approximations for the currents ⃗ i 0 1 and the average potential V 0 (for example, it can be equal to the values at the beginning of the step).
Otherwise, proceed to the 3.3 and 3.4.From Eq. ( 23) and ( 24), the increments ΔV and Δ ⃗ i 1 are calculated.The variables at the end of the iteration Proceeding to 3.1 with the changing of initial approximations.
After determining of the currents, the voltages on the capacitors are determined as:

The modeling results
Figures 5, 6, 7, 8, 9 shows the results obtained on the developed hybrid model, compared with the results of a physical experiment on a laboratory plant where the synchronous generator of type GCf84a/4, 27 kVA was used.The research was carried out for the regime of initial excitation of the generator without a load (Figs. 5, 6), for the load decrementing (Fig. 7), and for the load incrementing (Fig. 8) regimes, as well as for the steady-state regime (Fig. 9).The presented oscillograms show a high convergence of the results obtained on the hybrid model and on the physical plant.In particular, the maximum deviation is near 9%, and the average deviation is near 6%.
The other example of using average voltages in the integration step is described in [17].In this case, the developed nonlinear mathematical model is used for synthesis of fuzzy-control system for a power generation unit.

Conclusions
The application of the method of average voltages in the integration step to algebraize the differential equations describing the electrical circuits ensures the stability of mathematical models in their computer applications.
The combination of the method of average voltages in the integration step with the Newton method for solving systems of nonlinear algebraic equations allows to effectively simulate processes in nonlinear electromechanical systems.
In case of using the method in real-time models, the discretization step is determined by two factors: the rate of change of variables and the sampling of the signals in the physical part of the hybrid model.Within the physical signal sample, the step for the second-order method determines the time range in which the currents are described by the second-order polynomial and the voltages on the capacitances are described by the third-order polynomial.
The presented example is one of the numerous projects developed by the authors in the field of modeling of electromechanical systems.

Compliance with ethical standards
The physical part of the hybrid model is the real generator's excitation controller of type Digital AVR DECS 200 with an automatic voltage regulator.The signals of voltage U g and stator's current I g of generator are emitted from the computer model to the physical excitation controller.The output control signal of the pulse-width converter, which forms the pulse excitation voltage of generator, is sent to the computer model from the excitation controller.For the transmission of signals between physical and digital parts, the galvanic separation devices were used; these are not shown in the diagram.
is the vector of the currents in the branches; ⃗ e tr =  d ⃗ i dt -is the vector of the transformation e.m.f., -is the matrix of the dynamic inductances of the branches; ⃗ e r = p d d ⃗ i-is the vector of the rotation e.m.f., L-is the matrix of the static inductances of the branches, ω-is the angular speed, p-is the number of pairs of the poles, γ-is the rotation angle; ⃗ u C = u CA , u CB , u CC , 0, 0, 0 T -is the vector of the voltages on the capacitors; ⃗ e = 0, 0, 0, e f , 0, 0 T -is the vector of the e.m.f.applied to the branches.(

Fig. 2
Fig. 2 Functional scheme of the hybrid model of the power generation system of the ship's power plant

Fig. 5 Fig. 6 Fig. 7 Fig. 8
Fig. 5 Generator's output voltage (instantaneous values) during the initial excitation: a computer simulation; b experiment (the measurement sampling corresponds to 150,000 points per 6 s)