Numerical Methods Coupled with Richardson Extrapolation for Computation of Transient Power Systems

The numerical solution of transient stability problems is a key element for electrical power system operation. The classical model for multi-machine systems is defined as a set of non-linear differential equations for the rotor speed and the generator angle for each electrical machine, this mathematical model is usually known as the swing equations. This paper presents how to use direct Richardson extrapolation of several orders for the numerical solution of the swing equations and compares it with other commonly used implicit and explicit solvers such as Runge-Kutta, trapezoidal, Shampine and Radau methods. A numerical study on a simple three machine system is used to illustrate the performance and implementation of algebraic Richardson extrapolation coupled to several solution methods. Normally, 1 Universidad Pontificia Bolivariana, whady.florez@upb.edu.co, http://orcid.org/0000-0003-3977-0371 , Medellín, Colombia. 2 Universidad Pontificia Bolivariana, jorgew.gonzalez@gmail.com, Medellín, Colombia. 3 Universidad Pontificia Bolivariana, Ahill@upb.edu.co, Medellín, Colombia. 4 Universidad Pontificia Bolivariana, gabriel.lopez@upb.edu.co, Medellín, Colombia. 5 Universidad Pontificia Bolivariana, diegolopezjimenez277@gmail.com, http://orcid.org/0000-0002-4213-6825 , Medellín, Colombia. Universidad EAFIT 65| Numerical Methods Coupled with Richardson Extrapolation for Computation of Transient Power Systems the order of accuracy of any numerical solution can be increased when Richardson Extrapolation is used. A numerical example is provided for an electrical grid consisting of three machines and nine buses undergoing a disturbance. It is shown that in this case Richardson extrapolation effectively increases the order of accuracy of the explicit methods making them competitive with the implicit methods.


Introduction
The transient stability analysis in electrical grids is important to evaluate the effects of disturbances in rotors in power systems operation, security, and maintenance.Several researchers have applied numerical integration methods for solving the classical swing equations [1], [2], [3] to obtain the behavior of the power angles and the rotational speeds in the time domain.Classical numerical methods both implicit and explicit [4], [5], [6], [7] can be employed to find accurate results because numerical integration of the non linear differential equations is performed.The transient stability analysis (TSA) is very important for the design of electric power systems because it allows to determine the system's ability to withstand large disturbances and to transition to a normal operating condition.These disturbances can be faults such as: short and open circuits on a transmission line, loss of a generator, loss of a load, gain of load or loss of a portion of transmission network [8], [9], [10].For multi-machine stability analysis, the condition of the system and grid configuration before and after the disturbance, must be known.Therefore, two preliminary stages are needed before the transient study: the pre-fault conditions for the system which are calculated by a load flow balance, then the representation of the grid before the fault is modified to take into account the conditions during and after the failure [9].Regular numerical simulations are also important during planning stages to gain knowledge of the system.Yet, even a well designed system may face transient instability issues [11].The application and assesment of innovative numerical methods for solving transient stability situations is currently of interest because it allows the development of high accuracy computational packages and it may also have additional advantages for research and academical purposes.The solution of TSA problems involves two main stages, the network analysis and the dynamic model solution [2] .The network analysis includes the power load solution for determining the initial bus angles and voltages and the computation of the reduced admittance matrix.The dynamic model solves the transient differential equations for each machine in the system yielding the values of the angular velocities and the angles.These variables are important in studies of transient stability enhancement control [12], [13], [14].The set of differential equations, which can be highly nonlinear depending on the machine model used, are solved by a numerical method, for each time step.
In this paper the combination of the Richardson extrapolation [15] with different numerical methods is studied to determine its advantages respect to tradditional methods, such as the trapezoidal integration commonly used in power system dynamic equations [16].Initially, in this work the classical swing equations for rotor dynamics are solved with several numerical methods both explicit and implicit such as euler, trapezoidal (semi-implicit), Runge kutta, Shampine [4] and Radau [17].Then, the above methods are combined with Richardson extrapolation to evaluate and compare the accuracy and the order improvement of the resulting combined schemes.

Power system model
A multi-machine system has several synchronous generators which can be modelled by two differential nonlinear equations.This model is called the classical model or swing equations represented as follows: where n is the number of rotor machines, δ i is the rotor angle of the i − th machine, ω i is the angular speed, ω b is the base velocity of the system, M i is the moment of inertia, D i is the damping coefficient and P mi and P ei are the mechanical and electrical power respectively.The rotor angle is constant when the machines are in steady state when ω i = ω b that is also the velocity initial condition in the pre-fault state.The initial conditions of generator rotor angles and emf (E) are also obtained from the system prefault conditions [18].After an electrical disturbance the rotor angle changes from a constant value to another after the transient stage, assuming that stability is maintained.
The electrical power P e,i of machine i is calculated from where Y ij are the entries of a reduced admittance matrix [19], E i is the i-th voltage behind the transient reactance, and α ij are the admittance phase angles.Note that equations ( 1) and ( 2) can be seen either as a coupled algebraic-differential system or as a nonlinear differential system when equation ( 2) is subtituted in (1).The presented formulation allows

|68
Ingeniería y Ciencia reducing the number of unknown variables, because bus voltages as well as the equations of current injections at network buses are not needed [18], [3].During the stability study, the assumptions made are [16], [20], [21]: Mechanical power input is constant, damping is neglected, constant voltage behind the transient reactance model of the synchronous machine is valid, the mechanical rotor angle is assumed to be equal to the angle of the voltage behind the transient reactance, and the loads are represented by passive impedances.
An electrical network with n nodes from the standpoint of the generator terminals, has an admittance matrix Y defined by: To calculate the above matrix several steps are usually needed.First, the equivalent impedances or admittances are connected between the load buses and the reference node.Also, calculation of the fault impedances is added, and the admittance matrix is modified for each switching condition.Secondly, equivalent admittances are obtained for each impedance element.Finally, elements of the Y matrix are obtained as follows: Y ii is the sum of all the admittances connected to node i, and Y ij is the negative of the admittance between node i and j.After the admittance matrix has been obtained for the original network, it can be reduced by eliminating all but those corresponding to generators.This procedure known as Kron reduction [22] can be achieved by row-column operations considering that all the nodes have net zero current except for the internal generator ones.Thus, a reduced system is obtained as shown below where Now the system in ( 4) is partitioned to get ing.cienc., vol.13, no.26, pp.65-89, julio-diciembre.2017.
where n indicates the generator nodes and r is used for the remaining nodes.Now, expanding (6) and eliminating V r , the result is The reduced matrix (Y nn −Y nr Y −1 rr Y rn ) only contains the coefficients that correspond to the generator nodes.The network simplified equation represented in ( 7) is a convenient analytical technique that can be used when the loads are treated as constant impedances, this elimination procedure can be applied only to those nodes with zero current.The Kron reduction is the most important step prior to the transient solution.Also, it is a standard tool to obtain stationary and dynamically-equivalent reduced models for power flow studies, or in the reduction of differential-algebraic systems to lower dimensional dynamic models of power networks [22].

Richardson extrapolation combined with the numerical solution methods
The numerical solution of the power system equations can be improved by means of Richardson extrapolation which is an efficient technique for enhancing the accuracy of integration schemes [15].To explain how Richardson extrapolation can be applied to the solution of the swing equations, let consider A(h) the numerical solution of the system (1) which depends on the step size h, thus where a 0 is the solution that is being sought as the time step h approaches 0, i.e.A(h) → a 0 , and O(h p+n )is the Landau symbol representing the truncation error terms of order p + n an higher.Substituting different step sizes h,h/2,h/4,etc.into (8), the following system of equations can be assembled, For example, if only the first two terms of the series (8) are taken, the numerical solution method of order p would be written as and the following system is obtained solving the system of equations yields, where it is clear that the new approximation a 0 , being of order p + 1, will be more accurate than both A(h) and A(h/2) [15].Equation (12) shows that Richardson extrapolation can be used to accelerate the convergence of a numerical method.If a numerical method for solving a differential equation like Runge-Kutta, trapezoidal, etc. converges as O(h p ), the error can be reduced by a factor of h to produce a technique that converges as O(h p+1 ), and this can be used iteratively to produce other approximations that converge at any chosen rate.Hence, Richardson extrapolation can be used to obtain a higher order approximation for the dynamic equations ( 1) by increasing the order of any numerical method used to solved them.Note that the approximation a 0 given in equation ( 12) involves two numerical solutions, one with the total time step A(h), and another one with two shorter steps of size h/2 used to obtain the second numerical approximation A(h/2) at the end of the interval.Higher order Richardson extrapolation formulae can be deduced from equations ( 8) and ( 9) if more terms of the series (8) are retained.For example, if three terms of the series are used, the following system is obtained ing.cienc., vol.13, no.26, pp.65-89, julio-diciembre.2017.
which leads to the new approximation Similarly, if four terms of ( 8) are employed, a system of four equations is obtained whose solution gives the following higher order approximation It can be seen that to obtain the higher order approximations defined in equations ( 15)-( 16), the problem in each time interval must be solved several times but using different time step sizes.The first solution A(h) in one step with length h, the solution A(h/2) obtained with two smaller steps of size h/2, A(h/4) with four h/4 size steps, A(h/8) with eight h size steps and so on.Although this procedure increases the computational cost, it also improves the order of accuracy.Thus, if for example Runge Kutta method of order p = 2 were used to obtain the approximations A(h) and A(h/2), equation ( 12) would give a higher order solution of order p + 1 = 3.
Likewise, if the mentioned Runge Kutta method were used to obtain the approximate solutions A(h), A(h/2), A(h/4) and A(h/8), equation ( 16) would define a new approximation of order p + 3 = 5 higher than the original Runge Kutta used.Some stability theorems about the combination of different solution methods with Richardson extrapolation are discussed in [15].A short discussion of the solution methods used in the present work that can be combined with Richardson extrapolation is presented in the following section in order to understand the obtained results.

Overview of numerical solution methods for electrical power systems
There are several methods for the approximate solution of the swing equations such as: Euler's method, semi-implicit Euler method (trapezoidal), collocation, Runge Kutta, shampine, Radeau, among many more.Some of these methods are classified as implicit or explicit, in fact some methods

|72
Ingeniería y Ciencia like Runge-Kutta can be stated as either of these categories.Furthermore, the solution schemes may be classified according to the kind of problem i.e. there are special techniques for non-linear ordinary differential equations (ODE), and other for algebraic-differential systems.All of the mentioned methods are essentially iterative or sequential procedures that allow to calculate the power angle δ and the rotor velocity ω as a function of time, and to determine whether or not the system is stable.

Euler method
It is the most basic explicit method for numerical solution of ordinary differential equations and it can also be regarded as the simplest first order Runge-Kutta method.Most solution methods require the ODE to be written in the form for the case of the stability equations (1) the function vector is defined as and the Euler method defines the forward step where h is the time step size.The Euler method can be numerically unstable depending on the step size, especially for stiff equations.The implicit variant of the Euler method is defined as and it can clearly be observed that in this case the equations are non-linear and they need to be solved by Newton method or other similar iterative techniques to obtain the value of y n+1 .

Semi-implicit methods
The simplest semi implicit method is the trapezoidal integration rule [23] that is based on Euler method by averaging the derivative values, note that this method requeries non-linear iteration to find y n+1 value.
A semi-implicit method for the solution of the differential-algebraic equations (DAE) describing transient stability models, is proposed by Milano [3].DAE systems are formulated as g(t, y) = 0 As stated in equations ( 1) and ( 2), the function g(t, y) defines the electrical power P ei .This formulation coupled with an implicit integration scheme has two advantages when compared to the conventional explicit formulation: it reduces the computational burden and increases the sparsity of the Jacobian matrix of the system.

Shampine implicit method
An implicit method for the solution of sets of first order differential equations was developed by Shampine [4].In this method the problem must be formulated as The mentioned method is based on the Lagrangian form of polynomial interpolation and it is applied to improve the MATLAB built-in solvers (ode15i) [24].Furthermore, Shampine algorithm is usually considered as embedded method that produces an estimate of the local truncation error and as a result allows to control the error with adaptive stepsize.

Runge-Kutta methods
There are several families of Runge-Kutta methods (RK), both explicit and implicit and they may also have different orders.For example, the fourth order explicit RK methods is [23] The explicit RK methods can be generalized as The coefficients can be arranged in the following tableau form [25], The explicit RK methods have a bounded stability region [26] specially when stiff or discontinous problems are concerned.The inestability issues motivate the development of implicit RK methods [26] which are usually expressed as, Equations ( 27)-(28) are non linear and need to be solved to find the values of the k i coefficients.

Radau methods
These are fully implicit methods where the coefficient matrix (28) may have any structure, and they are usually appropiate for stiff systems.Radau methods are A-stable [5] but expensive to implement like all the implicit methods that need some sort of Newton iteration [17], [6], [5].Some alternatives have been proposed to reduce the computational cost that are usually based on approximate Jacobian matrix and dimensionality reduction [17].

Numerical results
In this section the proposed Richardson method of several orders equations ( 13),( 15) and ( 16), will be used to obtain the numerical solution for the IEEE 3-machine 9-bus system under a three-phase fault.The base MVA

|76
Ingeniería y Ciencia and system frequency are considered to be 100 MVA and 60 Hz, respectively.
Numerical experiments using different solution methods combined with Richardson extrapolation are presented to compare their performance.In the analyzed power network, Generator 1 is a hydro electric plant, while 2 and 3 are steam generators.Tables 1-2 show the load flow data needed for the solution, the electric power computation and the initial conditions.The simulated disturbance initiating the transient stage occurs near bus 7 at the end of line 5-7.The fault is cleared in 5 cycles (0.083 s) by opening line 5-7.obtain the initial bus voltages, the initial machine input and the power output.
• Compute the internal machine voltages • Compute the prefault system admittance matrix Y pre and its Kron reduction as in ( 7) (Table 3) • At the initial time t = 0 set the faulted bus voltages to zero • Compute the faulted system admittance matrix Y f ault by setting to zero the rows and columns of the initial matrix Y pre that correspond to the disturbed line.Also, the fault impedance is added as requiered and the admittance matrix is determined for each switching condition.Calculate the Kron reduction for Y f ault .
• Compute the post-fault system admittance matrix Y post and its kron reduction.Solve the differential equations by using different solvers alone and combined with Richardson extrapolation to compare the results.

|78
Ingeniería y Ciencia with both explicit and implicit solution methods to improve their order of accuracy.In Figure 1 the rotor angle variation with time is shown for each machine.This figure compares the results obtained by Euler explicit method with constant step width and by the Euler method improved by Richardson extrapolation with the reference solution given by DIGSILENT [27] and shown by [9].It is noticeable that the explicit Euler method with error order O(h 2 ) is improved by Richardson equation ( 16) with order O(h 5 ) according to what was explained in section 3, and the same trend can be observed for the angular velocity depicted in Figure 2  A plot of the angle differences of the machines in the system, δ 2 −δ 1 and δ 3 − δ 1 is shown in Figure 3, where it can be seen again that Richardson method considerably improves an explicit method such as Euler.Note that the solution is carried out in two swings to show that relative rotor angles swing together with respect to the simulation time and therefore the system is stable.A suitable criterion for stability is that if the rotor angle differences reach maximum values and then decrease, the system is stable [21].However, if any of the angle differences increases indefinitely the system is considered unstable because at least one machine would be out of synchronism.
where δ ref,i ,ω ref,i represent the reference solutions already mentioned.Note that the mean errors stay fairly uniform after the fault has been cleared.
The comparison of the absolute errors for different numerical solutions is shown in Figure 6 and 7 for the angle and velocity respectively.The total absolute errors for the three generators are defined as Angle error= Angular velocity error= where ndata is the total number of time values or steps in the numerical solution.Figure 4 and 5 present the comparison of some explicit and implicit methods, methods of different orders and the effect of Richardson extrapolation applied to some of these methods.Low order explicit methods such as Euler, which has order 1 i.e.O(h 2 ), yield higher errors than other explicit methods of higher order or the implicit methods.Note that the overall error is lower when the solution is calculated by methods of higher order, or when implicit or semi-implicit methods are used, or when Richardson extrapolation is combined with any method to improve its accuracy.Most of the methods presented in Figures 6 and 7 have constant time step except the Shampine method which includes adaptive control of error with variable step size and therefore its associated error is very small.Additionally, when Richardson extrapolation is combined with an integration method the results are noticeably improved and thus Richardson extrapolation can make explicit methods comparable with implicit and multi-step

|82
Ingeniería y Ciencia methods.In short, there are several alternatives to lower the numerical error: To use shorter integration step, to use adaptive step error control, to use high order methods either explicit or implicit and to improve the methods by coupling them with Richardson extrapolation.Also, the previous alternatives could be used combined among them to obtain highly accurate solutions.Usually, in most numerical methods as the step size becomes smaller, the less error is obtained.This behavior is observed in Table (4) where the fourth order Runke Kutta method combined with Richardson extrapolation, and the Shampine method were used with different step sizes to approximate the solution of equation (1).For the Runge Kutta method the step size value was kept constant on the entire integration interval, while due to its adaptive nature, in the Shampine method the step value h refers to the maximum initial step size allowed.The order of the Runge Kutta method was increased by combining it with Richardson extrapolation according to formulae ( 13),( 15) and ( 16), hence theoretically the combined methods would have errors of order O(h 6 ),O(h 7 ) and O(h 8 ) respectively.For example, it can be seen, for an explicit method like Runge Kutta of fourth order (RK4) combined with Richardson of error order O(h p+1 ), the mean error goes from 1.41 × 10 −2 to 3.40 × 10 −3 for step sizes h = 1 × 10 −2 and h = 1 × 10 −4 respectively.Hence, the reduction of the step size produced almost 4 times less error.Similar error reduction can be accomplished if the Richardson order increases from O(h p+1 ) to O(h p+3 ) with a constant step of h = 1 × 10 −2 , the error value 1.41 × 10 −2 lowers to 2.84 × 10 −3 which means a reduction almost by a factor of 5. Similar results can be observed in the other rows and columns of Table 4 that illustrates how, from an accuracy standpoint, an explicit method improved by Richardson extrapolation can be made competitive with other more sophisticated implicit methods such as Shampine.

Method
The error reduction effect associated with Richardson method can be confirmed in Figure 8.It can be seen that when the order of the Richardson extrapolation is higher the error becomes smaller regardeless of the integration method being explicit or implicit.However, this improvement in accuracy comes at the expense of computational time as presented in Table 5, which shows that the cpu time becomes higher either when the order of the integration method is high or the order of Richardson is increased.The lower order methods need few operations, therefore their CPU time is lower.In contrast, implicit methods take longer not only because they need more operations to advance one step, but also because they require Newton iterations to account for their inherent non-linearity.Note that when the Richardson order is increased, the CPU time is doubled.However, it can be seen that a lower order explicit method combined with Richardson becomes comparable in time and accuracy with a higher order implicit method.For example, an Euler method with Richardson extrapolation of error order O(h p+3 ) takes 0.23 seconds, while the trapezoidal spends 1.38 seconds to achieve less accuracy than the improved Euler method as was shown in Figure 6 and 7.

Conclusions
This paper presented the numerical solution of a transient stability problem for the IEEE 3-machine 9-bus system using different numerical methods.The integration methods are compared, and the effect of combining them with Richardson extrapolation of several orders is analyzed.The obtained numerical results are in agreement with a reference solution calculated by software Digsilent and also presented in several literature references.The proposed Richardson extrapolation schemes of several orders are presented based on an algebraic approach rather than in the conventional recursion method.It was found that Richardson approach really improves the explicit methods making them competitive with the implicit ones.Also, although Richardson method increases the computational time when increasing the

|86
Ingeniería y Ciencia order of the approximation, the solution time of the explicit methods improved by Richardson stays well below their implicit counterparts, while the accuracy of the solution improves considerably.

Figure 6 :Figure 7 :
Figure 6: Overall angle error comparison for different methods

Figure 8 :
Figure 8: Effect of increasing the order of Richardson extrapolation.

Table 1 :
Load flow for IEEE 3-machine 9-bus system

Table 2 :
Generator internal voltages and their initial angles The steps for the numerical experiments are the following:• After modeling the network, a load flow analysis is performed to ing.cienc., vol.13, no.26, pp.65-89, julio-diciembre.2017.

Table 3 :
Reduced admittance matrices numerical solution of the transient problem that describe the power network dynamics during a disturbance.The main purpose of the program is to show how Richardson extrapolation can be used together .

Table 4 :
Effect of the step size on several solution methods

Table 5 :
CPU time for the numerical solution by different methods