Multirate and multi-stage holomorphic embedding for dynamic simulation of systems of differential-algebraic equations

In this paper, we present a novel multirate (MR) extension to the multi-stage holomorphic embedding (MSHE) method and use it to dynamically solve systems of differential-algebraic equations. As a test of real-world applicability, we analyze the algorithm for typical power system cases. When applied to the IEEE 34 bus and IEEE 9 bus benchmark system, the multirate approach significantly reduces the number of recalculations of most variables as compared to the existing MSHE algorithm, thus offering the chance of decreasing simulation time while providing simple means of error control. This integration of a multirate ansatz into the well-established MSHE solver potentially opens up systems for dynamic simulation that were up to now out of numerical reach.


Introduction
Dynamic simulation of power systems consists of solving initial-value problems of systems of differential algebraic equations (DAE).While conventional numerical approaches (e.g.Newton-Raphson or Gauss-Seidel) often struggle with slow convergence and also do not even guarantee convergence to a physical solution [1], the Holomorphic Embedding (HE) mechanism introduced by Trias in 2012 is able to solve both issues and offers theoretical guarantees on convergence, robustness and efficiency [2,3].It was soon extensively applied to loadflow problems [4][5][6][7] and adapted to dynamic simulation of power systems [8,9] where it is a topic of ongoing research [10].Its basic idea to approximate the variables in the system by complex power series can additionally be expanded by the Padé approximation [11,12].
When simulating large power systems, the occurrence of a wide range of time scales usually leads to fast varying variables dictating the simulation's step size and speed.In order to reduce the resulting computational load from solving variables on slow time scales, it is desired to logically separate different scales and solve each variable only with the accuracy required to capture its own dynamic behavior.Multirate algorithms address this phenomenon by effectively decoupling variables and updating them at different rates.They have been an active area of research for more than 4 decades and continue to be of interest in many fields [13][14][15][16][17][18][19][20].
This publication presents a new multirate approach in combination with the multi-stage HE methodology introduced in [9].Due to its capability of approximating the analytical solutions in the form of power series, HE is uniquely suited for the proposed multirate method.
We demonstrate that the usage of our multirate approach together with the application of Padé approximation significantly reduces the necessary amount of recalculations compared to the approach in [9].This shows the algorithm's potential of greatly reducing computational load and speeding up dynamic simulation considerably whilst taking advantage of the robustness and convergence properties of the HE formulation.
The paper is structured as follows.First, section 2 provides the description of the novel multirate technique.In sections 2.1 and 2.2 the basic HE approach and Padé approximations are revisited briefly.Then, we describe the simulated test cases, from simple toy problems to the IEEE 9 bus test case in section 3. Finally, section 4 presents our results, indicating a consistently lower number of recalculations of the HE coefficients than with the non-multirate version of the algorithm, and section 5 provides summarizing thoughts and an outlook into future research.

Methods and tools
In this section, we shortly describe the holomorphic embedding (HE) method for solving differential-algebraic equations (DAEs) in combination with Padé approximations.For more in-depths explanations, refer to References [9,21,22].Afterwards, our extension of this integration scheme with a multirate approach is explained in more detail.

Holomorphic embedding
We follow Reference [9] for the description of the holomorphic embedding process.Throughout this section, we will use the following definitions. .
The functions x and y are considered to be holomorphic in a neighborhood A of α 0 .We also assume that f and g are polynomial.Now, a typical system of DAEs is given for Additionally, we have some initial condition x(α 0 ) = x 0 , y(α 0 ) = y 0 .The solutions for x(α) and y(α) on the whole subset A can be approximated by truncated power series of order Î  N , designated x ˜and y˜: Expanding around α 0 , we can calculate the coefficients x[k] and y[k] using the following scheme.First, rewrite the system of DAEs (1) according to table 1.
In the general case, the transformed system of equations is given by where F and G denote the equations arising from applying the transformation rules in table 1 to f and g.Given the initial condition x 0 and y 0 , we can iteratively solve (2) for all coefficients up to arbitrary order N.The resulting truncated power series for x ˜and y˜have a finite convergence radius outside of which they become invalid (in the sense that they do not fulfill the DAE anymore).The full solution is therefore obtained by combining multiple locally valid solutions into a multi-stage solution.
Since we are interested in the dynamic simulation of power systems, assume that we restrict the domain A w.l.o.g. to a given time interval t t , f 0 [ ]on the domain , so Î  t .We define the vector-valued residual error Δ (t) of the system of DAEs as and the vector-valued relative residual error ò as Note that both definitions do not compute the deviation of the individual approximations from their true solutions, but instead the deviation of each full equation from its true value zero.It becomes clear that this error measure entails an inherent dependence of the algorithm on the internal structure of the DAE.Both Δ and ò can be chosen as measure for the DAE imbalance, the choice should be based on the values of the DAEs solutions.Although ò does not directly correspond to the relative deviation between the approximation and the true solution of a variable, it aims to take the involved orders of magnitude into account to some extent by dividing through the derivative ¢ x .However, the particular choice of this denominator is rather free, any other summand in f (x(t), y(t)) could be equally suitable.During our tests, no specific choice neither between Δ and ò, nor between different denominators in equation (4), showed significant advantages.Now we can find the value t max on our integration path, where either ||Δ|| ∞ or ||ò|| ∞ first exceeds a predefined error-tolerance  max .The value t max marks the transition from one solution to another.Evaluating the truncated power series of the current stage at t max yields initial conditions that can be used to obtain a new approximation for the next stage.In this manner the system of DAEs can be solved on the full simulation interval t t A description of relevant parameters can be found in table 2.

Padé approximation
In order to improve the radius of convergence of each truncated power series in x ˜and y˜, a Padé approximation (PA) can be applied to each approximation [11,12].In general, a Padé approximant to a scalar function f, denoted as [L/M] f , is a ratio of two polynomials p and q with p having degree L and q having degree M that fulfills the following condition: If f is given as a truncated power series (as it is the case for the approximations in x ˜and y˜), equating the coefficients will lead to a matrix equation whose solution yields the coefficients of p and q: In this paper, the algorithm described in [23] is employed to calculate the Padé coefficients due to its capability of producing numerically stable results via singular value decomposition [21,22].

Multirate approach
In [9], when  max is exceeded at a transition point, all variables are recalculated by solving the equations for the next stage.But depending on the DAE and its solution, it is possible that certain approximations in the system have much larger radii of convergence than others [24,25].One truncated power series could be perfectly valid for the next three stages while rapidly changing variables force the simulation to recalculate them nevertheless.In order to take advantage of this phenomenon, we adapted the algorithm to recycle solutions from the previous stage that are still valid within the error bounds for the next stage instead of calculating a complete new set of approximations with the HE mechanism.This approach results in different recalculation rates for the approximations of the variables in the system.The implementation works as follows: 1. Find the convergence radius of each equation via binary search over the whole simulation interval t t The search is terminated when the next step length is lower than a predefined parameter t step and the left boundary is taken as convergence radius.Step size in search to determine t max .δ Minimal improvement required to trigger the multirate shift.η Parameter for absolute multirate improvement condition.One possibility to determine δ. ξ Parameter for relative multirate improvement condition.The other possibility to determine δ.
2. Identify offending equations that do not exceed the minimal convergence radius t max by a certain amount δ (see also figure 1(a)).
3. Tag variables that do not occur in the offending equations as valid variables.
4. Transform coefficients of valid variables from previous stage into coefficients of power series around t max by applying a coefficient shift.

5.
Resolve offending equations for next stage using the known power series of the valid variables around t max .
In order to be superior to a non-multirate simulation, steps 2-4 should be faster than solving the whole system of DAEs.For simulations without Padé approximants this can be achieved through a coefficient shift of the power series [26]: be the k-th order coefficient for an expansion around t 0 .Coefficients for an expansion around t 1 are then obtained by a simple matrix multiplication: Coefficients of Padé approximants cannot be transformed as easily.The original power series that was obtained by the solver is not necessarily valid at transition point t max .This behaviour is in fact desired as it justifies the application of the Padé approximation with the goal of increasing the convergence radius of the approximations.Therefore one cannot do such coefficient shift to recalculate the corresponding Padé approximant, but rather needs to expand the Padé approximant into a power series around expansion point t max .Due to the fact that Padé approximants are simple rational functions this can be done in a straightforward manner.We use the algorithm outlined in [27] which exploits synthetic division of polynomials.

Multirate improvement condition
The value of the minimal improvement δ, depicted in figure 1(a), that is required to trigger the multirate algorithm can be determined either by absolute improvement

d x =
or as relative improvement compared to the minimal HE step This introduces the parameters ξ and η which inherently have a significant impact on the multirate algorithm.It should be noted that they cannot be used simultaneously.All variables in offending equations are resolved by the HE mechanism at the transition point t max , effectively resetting the convergence radius of these offending equations to t max .The complete multirate HE algorithm is illustrated in figure 2.An overview of the relevant parameters is given in table 2.

Application on example systems
The successful application of the HE algorithm to dynamic simulation of power systems and the benefit in comparison to other numerical solvers has been shown [9].This paper focuses on the performance of our novel multirate HE approach with respect to the classical HE method as a proof of concept.In this section, several example DAEs are solved using different combinations of the methods outlined above.Both small simple toy models and real-world IEEE benchmark systems are chosen as test cases in order to ensure that the algorithm produces consistent results independent of system size.

Toy models
As first simple tests, two small and simple systems are considered.Many DAEs of higher order and varying complexity can be formulated as linear first order systems by substitution of variables which is why those systems are highly relevant for the assessment of the proposed method.Therefore the first example consists of two weakly coupled linear first order systems with different time constants.
This example is denoted as 'toy model 1'.Note that A 11 only possesses purely real eigenvalues whereas the eigenvalues of A 22 have both real and imaginary parts, i.e. the subsystems have exponential and damped oscillatory solutions, respectively.Another class of ODE systems with a wide area of application in various fields is the harmonic oscillator.We formulate the equations for a system of 4 classical oscillators that are coupled through springs and driven by external forces on both ends (see figure 1(b)).Since one of the main purposes of the multirate extension is to correctly recognize distinct time scales in the system as well as mixed coupling strengths between the equations This system is denoted as 'toy model 2'.

Power system test cases
In order to evaluate and quantify the performance of the algorithm when applied to real-world systems, the linearized IEEE 34 node test feeder is simulated with the outlined methods [28][29][30].where x are the independent state variables (here capacitor voltages and coil currents) of the linearized system, u are the input variables (in this case the voltages of the three phase generator) and y are the output variables (in this case voltage or current measurements).Note that the solutions of all variables are sinusoidal due to the generator input voltages being modelled as a sine.
In addition to this steady-state system we also tested a dynamical real-world power system.To this end the IEEE 9 bus test system is simulated with a generator-side fault, namely a 10% reduction of the mechanical input to generator 1 from t = 1 until t = 2, inducing the typical electromagnetic transient dynamics on distinct time scales.This type of simulation is typical in stability analysis.The equations modelling the dynamic behaviour of a multi-machine system are taken from [31] and brought into suitable symbolic form as described in [33].Following [33], the power flow data, transmission lines and dynamic machine data are found in [31] and the machine sub-transient and turbine parameters are estimated from [34].The MATLAB code pertaining to [33] was published and was used to generate the equation system solved in this paper.After some variable substitutions, the system features 60 variables and equations.A single line diagram of the IEEE 9 bus test case is depicted in figure 3(b).

Results
The number of recalculations of each variable (denoted n c ), i.e. the amount of transition points where the power series coefficients need to be recalculated via the HE algorithm in order to obtain an approximation for that variable, is chosen as the primary measure of performance for the context of this publication.A coefficient shift performed in the multirate algorithm is thus not considered as a recalculation, since the simple matrix multiplication necessary to perform it requires negligible effort compared to the computation of the coefficients via the HE solver.No focus was laid on optimization and therefore run times are not discussed in this publication.This paper serves as a proof of concept for the multirate approach rather than a demonstration of superior performance to other well-established solvers.However, all numerical solutions are cross-checked and validated against MATLAB solvers.
Four simulation methods are compared: Usage of power series approximations with and without usage of Padé approximants, respectively with or without the novel multirate approach.Simulations utilizing only the power series representation are denoted by 'Normal' whereas simulations working with Padé approximations are referred to by 'Padé'.When using the novel multirate approach, the two methods are referred to as 'Normal MR' and 'Padeé MR', respectively.
The parameters for the simulation of the test systems are listed in table 3.

Toy models
We start by evaluating the simple toy models introduced in section 3.1.The simulation results are presented in figure 4.
In the case of toy model 1, the multirate algorithm leads to less recalculations and recognition of different timescales when paired with Padé approximation.Without it, all recalculations happen simultaneously for all variables, indicating a poor choice of η-see also section 4.2 directly below, where we explore the influence of η on the multirate capabilities of the algorithm.In this case, increasing η such that the single-rate solution is recovered is the least-effort option to prevent performance degradation.For toy model 2, the multirate algorithm is able to detect the structural independence of the two sets of oscillators and successfully delivers a lower number of average recalculations.We can also see that the positions x 1 , x 2 , x 3 of the 3 oscillators coupled strongly to the left wall (with lower oscillation frequency ω 1 ) need to be recalculated less often than the positions x 4 , x 5 , x 6 of the oscillators coupled strongly to the right wall (with higher oscillation frequency ω 2 ).The weak coupling in the middle effectively couples the different time scales and leads to an increase of needed recalculations for both adjacent oscillators.3.For (b), both versions of the multirate method recognize distinct time scales and produce consistent results with lower or equal average n C .For (a), the multirate method without Padé approximation leads to a higher number of recalculations while also using the same rate across state variables.This may be a sign of badly chosen η for this specific case; as a simple example, it is possible to choose η so big that the single-rate solution is recovered.The simulation results for the IEEE 34 node test case described in section 3.2 are depicted in figure 5 (system 1) and figure 6 (system 2) respectively.For system 1 it is found that both Padé versions of the algorithm need persistently less recalculations than their Normal counterparts.Also, although single variables exhibit a significantly increased amount of recalculations, both multirate methods again decrease average n C in comparison to their non-multirate counterparts by dragging down the baseline.The majority of variables benefits from multirate and only a few variables need smaller time steps and more recalculations.A closer look at the six peaking variables shows that all of them represent variables of the line between bus 808 and 810 (namely, input and output voltages as well as coil currents), indicating the physical location for the reason in the region of the first branch of IEEE 34.Another important observation is that the differential state-space equations of these variables are closely linked: 5 of the 6 contain exclusively peaking variables.This demonstrates that the form of the equation system-which variable occurs with which coefficient in which equation-inherently has a strong influence on the multirate algorithm.Regarding system 2 both Padé versions show the same desired behavior of decreasing average n C .Nevertheless it can be seen that usage of MR increases the total number of recalculations when used in 'Normal' mode.This could again originate from a bad choice of η and shows that the multirate improvement parameters need further inspection.We investigate them by means of test case 'system 1'.  ).Using Padé approximations, the multirate approach is able to reduce or leave invariant the number of necessary recalculations for every variable.In the 'Normal' case, recalculations of the same variables as in the Padé case can be reduced, but only at the cost of increasing recalculations for the other variables, leading to a higher average n c .

Parameter scan η
As indicated before, the parameter η has a large impact on multirate performance.If chosen too small, very small convergence radii can lead to "oscillating" violations of the error tolerance: If at least one equation is valid only slightly longer than t max (relative to t step ), this leads to an almost immediate recalculation directly after the previous one.Then the coefficient shifts performed at the transition point can result in the original (un-shifted) coefficients being invalid soon after that, hence triggering another recalculation.If too many such coefficient shifts are executed consecutively, they can become numerically unstable and lead to repeated recalculations although the original approximation would have still been valid without the repeated coefficient shifts.This effect can drastically decrease the effective step length and thereby increase the number of recalculations for multiple variables.In contrast, if η is chosen too large, it will prevent the trigger of multirate events altogether and therefore reduce a multirate simulation to a single-rate simulation-which may be desired in certain scenarios (e.g. the results for toy model 1 above).Figure 7 shows the number of recalculations for each state variable in the small IEEE 34 test case (system 1) for different values of η.It is evident how sensitive the simulation reacts to the choice of η.Simulations with low values of η exhibit the behavior seen above, namely a low baseline with high peaks.This characteristic shape fades into lower peaks and elevated plateaus as η is gradually increased, leading to a high average n C , and eventually evolves into a non-multirate simulation when η = 3.0.

Parameter scan ξ
The parameter ξ has a similar impact on the multirate method as η, since they both are measures for the improvement condition.An analogous study about the influence of ξ shown in figure 8 thus features similar characteristics as in figure 7.But it is important to note that in principle ξ offers a more intuitive and individual control over the simulation, as it can be directly related to the relevant time scales in the system.If these are known a priori, ξ can be utilized to suppress any unwanted multirate behavior by setting it accordingly.Additionally, if the time scales are not known beforehand, the behavior of the multirate algorithm can allow conclusions about the system.All in all this analysis shows that while η and ξ need to be chosen carefully, the overall simulation performance can potentially benefit greatly from our multirate approach.

General results
The results for the IEEE 9 bus test case are depicted in figure 9.The Padé methods are again able to reduce the average number of recalculations while in Normal mode the average n C is drastically increased.This behavior seems to be consistent over multiple system formulations, illustrating once more the advantage of the Padé version over the Normal version.

Dynamics
Furthermore, it is of great interest to investigate the dynamic behavior of the outlined methods.In order to understand how the fault and the subsequent dynamics reflect in the number of recalculations, a timedependent moving average over a specified time window is computed.The results for selected variables are shown in figure 10.It can be observed that the Padé versions n c emulate the solution dynamics: Following the fault, a fast oscillation is accompanied by increased number of recalculations, which decreases as the solution stabilizes again.In contrast, all variables exhibit the same behavior when simulated with the Normal MR algorithm, namely a center peak that has no obvious origin in the solution dynamics.Once again, a different choice of the multirate improvement condition could at least recover the single-rate solution.
Figure 10(d) points out another important result: Although the variable ¢ E p,3 has the value zero over the whole simulation interval, its approximation is not valid over the same interval which may contradict intuition.Highly active equations containing ¢ E p,3 force it to get recalculated although it is itself constant.This demonstrates nicely the impact of the underlying structure of the equation system on the multirate performance that can only be reduced to a certain degree, e.g. by analytical simplification of the system.

Conclusion and outlook
This paper presents the successful integration of our novel multirate extension into the multi-stage holomorphic embedding formulation and the application to dynamic simulation of simple toy systems as well as the realworld IEEE 34 bus and IEEE 9 bus systems.The number of recalculations for each variable is selected as performance measure to establish a proof of concept for the method.A direct comparison between simulations using only power series and those also using Padé approximations and between the respective multirate and non-multirate versions is provided, verifying the advantage of Padé approximants over pure power series in all  E p,3 is recalculated although its value is constant.This illustrates the influence of the equation system structure.
cases.When applied to real-world systems, the new multirate approach is able to reliably decrease the average number of recalculations, proving its potential to reduce the overall computational load and simulation time during dynamic simulation of systems.Performance optimizations that are necessary for quantitative tests of this potential will be shown in a future publication.
In the future, the presented multirate variant of the multi-stage HE technique may present further advantages than simply computation speed.Due to the calculation of convergence radii inherent to the method, statements about the stability of the simulated system may be inferable, possibly even attributable to individual equations thanks to the multirate nature of the approach.Furthermore, the order of the method may be variably increased or decreased based on the convergence radius, to dynamically increase or decrease not only the timestep, but the approximation order itself, on a per-equation basis.Additionally it would be interesting to investigate further the influence of the equation structure, error measure and possibly equation-bound improvement parameters on the multirate performance.

Figure 1 .
Figure 1.(a) Improvement condition imposed on multirate algorithm to identify offending equations.The equations convergence radii are marked by short black lines.At each transition all variables occurring in offending equations are marked invalid and re-solved via HE, effectively resetting their convergence radii to a max .(b)Drawing of the coupled oscillators modeled by toy model 2. The system is driven by external forces from both sides and coupled through springs with different constants.

Figure 2 .
Figure 2. Multirate HE algorithm.If the Padé approximation is not applied, that step is skipped, proceeding directly with determination of a max .If multirate is switched off, the multirate improvement condition and the coefficient shift are skipped.

Figure 3 .
Figure 3. Overview of the power system test cases.(a) IEEE 34 node test feeder and areas of the real-world examples.System 1 corresponds to a subregion, System 2 to the complete test feeder.Adapted from [28].(b) Single line diagram of the IEEE 9 bus test case.[33].

Figure 4 .
Figure 4. Number of recalculations n C of the state variables for simple toy models: (a) coupled linear first order ODE system; (b) coupled oscillator system.Simulation parameters are listed in table3.For (b), both versions of the multirate method recognize distinct time scales and produce consistent results with lower or equal average n C .For (a), the multirate method without Padé approximation leads to a higher number of recalculations while also using the same rate across state variables.This may be a sign of badly chosen η for this specific case; as a simple example, it is possible to choose η so big that the single-rate solution is recovered.

Figure 5 .
Figure 5. Number of recalculations n C of the state variables for the IEEE 34 test feeder system 1.Both multirate versions recognize distinct timescales and produce consistent results with lower average n C than their single-rate counterparts.

Figure 6 .
Figure 6.Number of recalculations n c of the state variables for the full IEEE 34 test feeder ('system 2').Using Padé approximations, the multirate approach is able to reduce or leave invariant the number of necessary recalculations for every variable.In the 'Normal' case, recalculations of the same variables as in the Padé case can be reduced, but only at the cost of increasing recalculations for the other variables, leading to a higher average n c .

Figure 7 .
Figure 7. Number of recalculations n C of the state variables in system 1 test case with varying parameter η when using (a) Normal simulation mode and (b) Padé simulation mode.All other parameters are kept constant as in table3.Except for the dashed blacked line which corresponds to a non-multirate simulation, multirate is always switched on.The average number of recalculations is lower compared to non-multirate when η is chosen below 0.3 (Normal) or 0.4 (Padé).With η = 3.0 the simulation reduces to nonmultirate.

Figure 8 .
Figure 8. Number of recalculations n C of the state variables in system 1 with varying parameter ξ when using (a) Normal simulation mode and (b) Padé simulation mode.All other parameters are kept constant as in table 3. Except for the dashed blacked line which corresponds to a non-multirate simulation, multirate is always switched on.The average number of recalculations is lower compared to non-multirate only if η is 1 • 10 −6 (Normal) or below 2 • 10 −6 (Padé).Choosing 2 • 10 −5 reduces the simulation to non-multirate.

Figure 9 .
Figure 9. Number of recalculations for the IEEE 9 bus test case with fault.The Padé versions are able to reduce n C when MR is used, the 'Normal' algorithms show a significant increase in n C when using multirate.

Figure 10 .
Figure 10.Dynamic of the multirate algorithm presented in form of a moving average (blue and orange shaded areas for Normal and Padé respectively) of selected variables with the dynamic solution of each variable depicted by blue and red lines.(a) to (c) demonstrate the ability of the Padé methods to map the solution dynamics onto n C .(d) ¢E p,3 is recalculated although its value is constant.This illustrates the influence of the equation system structure.

Table 2 .
Overview of relevant parameters and variables for the multirate multi-stage HE algorithm.Truncation order of each variable's Padé approximation. max Maximum error for all equations.If an equation is correct within a margin of error  max , it is considered valid.t max Transition point from one HE solution to another.t step

Table 3 .
Parameters for simulation of example systems.The relative error ò is used as error measure in all cases except IEEE 9 which uses the absolute error.The time intervals for the IEEE 34 systems are chosen due to computational resources.