Small-Signal Stability Analysis of Implicit Integration Methods for Power Systems with Delays

—The paper focuses on the accuracy and stability of implicit numerical methods when utilized for the Time Domain Integration (TDI) of power systems with inclusion of time delays. In particular, a small-disturbance analysis technique is proposed to evaluate the numerical distortion that TDI methods induce to the dynamic modes of power systems modeled as Delay Differential Algebraic Equations (DDAEs). The case study illus- trates the features of the proposed technique through simulations conducted using the IEEE 14-bus test system, and considering three examples of implicit integration methods, namely Backward Euler Method (BEM), Implicit Trapezoidal Method (ITM), and 2-stage Radau IIA.


I. INTRODUCTION A. Motivation
Power system modeling and stability analysis in the presence of time delays has attracted increasing attention in recent years, mainly due to the potential negative impacts that measurement and communication latency can have on the stability of automatic regulation loops, e.g. in wide area damping controllers [1]- [7]. To date, the most successful technique to evaluate the effect of time delays on the dynamic behavior of a power system following a disturbance is to carry out a TDI of the system's equations through a suitable numerical scheme. However, it is known that delays can worsen the accuracy and even make unstable otherwise very robust numerical schemes, such as the Implicit Trapezoidal Method (ITM) [8]. It is thus crucial to have reliable tools able to evaluate the accuracy of the numerical schemes employed for the TDI of dynamic power system models with inclusion of delays. This paper proposes a novel framework to address precisely this issue.

B. Literature Review
The conventional dynamic power system model is formulated as a set of non-linear Differential Algebraic Equations (DAEs) [9], [10]. These equations are known to be stiff and their TDI is typically conducted through an implicit numerical scheme, in order to guarantee that numerical stability is maintained for every integration time step size. Including time delays to the power system model changes the nature of its equations, leading to a set of non-linear Delay Differential Algebraic Equations (DDAEs) [2], for which many of the known results about the stability properties of TDI methods for DAEs lose their validity.
The simplest approach to numerically integrate a set of DDAEs is by modifying standard implicit methods so that they account also for the delayed terms. In this regard, the modifications required by the ITM to integrate power system models with time delays were discussed in [2], [11]. On the other hand, theoretical results from numerical analysis indicate that such approach may not be adequate and suggest that accurate integration of Delay Differential Equations (DDEs) requires the application of special methods developed to this aim, e.g. see [8], [12]. In this regard, we note that most insights on the stability characterization of a TDI method for time-delay systems are derived based on the behavior of the method when applied to a scalar, linear test DDE. The major limitation of such approach is that results provide only rough and qualitative information, since they are typically not generalizable for systems of DDEs. As a matter of fact, it has been proven that no A-stable natural Runge-Kutta (RK) method is asymptotically stable on the whole class of asymptotically stable linear systems of DDEs [8], [13].
The stability and precision of implicit numerical methods employed for the TDI of DDAE power system models is a problem that has not been systematically discussed in the literature. Regarding the precision, the standard tool used by most solvers is truncation error analysis, which, although it provides good estimates of the deviation between exact and numerically computed trajectories, it yet cannot capture the ability of a method to prevent the exponential growth of truncation errors. The goal of this paper is to provide a novel approach to assess, in a unified way, the stability and accuracy of numerical methods employed for the TDI of DDAE power system models.

C. Contributions
The contributions of the paper are twofold: • A novel framework, based on Small-Signal Stability Analysis (SSSA), is proposed to evaluate the numerical approximation introduced by TDI methods, when applied for the integration of power systems modeled as DDAEs.
• A discussion on how the presence of delays impacts the spurious oscillations and/or overdamping introduced to the dynamic modes of power systems by the Backward Euler Method (BEM), the ITM, and the 2-stage Radau IIA method. To the best of our knowledge, this is the first study that aims to systematically study and observe this numerical effect.

D. Organization
The remainder of the paper is organized as follows. Section II outlines the modeling and stability analysis of power systems impacted by time delays. Section III describes the proposed approach to evaluate the accuracy and stability of TDI methods. Section IV presents a case study based on the IEEE 14-bus benchmark system that showcases important features of the proposed approach. Finally, conclusions are drawn and future work is outlined in Section V.

A. DDAE Model
The dynamic power system model is conventionally described by a set of DAEs, as follows [9]: where x = x(t) ∈ R n and y = y(t) ∈ R m denote the state and algebraic variables, respectively; f : R n+m → R n , g : R n+m → R m , are non-linear functions; and 0 m,1 denotes the m×1 zero matrix. Discrete variables in (1) are represented implicitly, i.e., a discontinuous change in the system gives rise to a jump from (1) to a new set of equations of the same form. Inclusion of delays in the right-hand side of (1) changes the DAEs into a set of DDAEs of retarded type, as follows [2]: where x d ∈ R n d , y d ∈ R m d , are the delayed or retarded state and algebraic variables, respectively. In compact form, (2) can be rewritten as follows: , and:

B. Small-Signal Stability Analysis
The proposed approach to assess the accuracy and stability of numerical methods for the TDI of power systems with delays, which will be described in Section III, is based on SSSA. Thus, we first introduce the reader to the SSSA of power systems with inclusion of delays. For simplicity, let assume that the delay system includes a single constant delay τ > 0. Then, we have that x d = x(t − τ ) and (3) becomes: Linearization of (5) around an equilibrium point gives: where A 0 and A 1 are the Jacobian matrices associated to the delay-free and delayed variables of the system, respectively; andx indicates the deviation of x from the equilibrium. If the system includes multiple, say ν, delays, τ k > 0, k = 1, 2, . . . , ν, the last expression is generalized as: The eigenvalues of system (7) are determined from the solution of the corresponding characteristic equation [14]: where s ∈ C. Then, (7) is stable if and only if all finite eigenvalues have negative real parts. The presence in (8) of the exponential function implies the existence of infinitely many roots [14]. In this paper, this problem is solved with the technique proposed in [15]. That is, (7) is transformed to an equivalent system of Partial Differential Equations (PDEs) of infinite size. Then, the PDE system is reduced to a finite dimensional problem through Chebyshev's spectral discretization. Chebyshev's discretization technique has been used for the eigenvalue analysis of delayed power system models, e.g. in systems with constant and stochastic delays affecting the stability of control loops [2], [6], [16] and has proved to show very good accuracy if a proper number of interpolation nodes is selected. A sparse version of the same technique has been recently proposed in [17].

III. PROPOSED APPROACH
A TDI method for power systems with time delays is a discrete-time approximation employed to solve system (3) for a defined time period and set of initial conditions. In this section, we describe the proposed approach to evaluate the amount of approximation introduced by implicit TDI methods to the representation of the dynamic modes of system (3).

A. Single Delay
Assume for simplicity that the integration time step h is constant, and that the system includes a single constant delay, which is an integer multiple of the integration time step, i.e. τ = ch, where c ∈ N. We provide the following definition: Definition 1. In an implicit form, a TDI method applied to system (3) can be described by a discrete-time system, as follows: where r = n + m; η : R 4r → R r is a vector of non-linear functions; and x t : The discrete-time system (9) covers the family of RK methods, i.e. the most important family of methods for the integration of DDEs. In general, RK methods have been shown to be preferable to linear multi-step methods for the integration of DDEs, since they are characterized by simplicity of implementation, low computational complexity and high accuracy, e.g. see [8] and relevant references therein. Linearization of (9) gives: The linearized method can be rewritten as follows: where C i are, in general, matrix functions of h, E, A 0 , A 1 . We arrive at the following proposition.
Proposition 1. The stability of (11) can be assessed by studying the stability of the following linear discrete-time system: where: with C † having dimensions r × cr.
The proof of Proposition 1 is provided in the Appendix. Then, the stability of (12) can be studied by finding its eigenvalues, which are the roots of the characteristic equation: whereẑ ∈ C. In particular, (12) is asymptotically stable if and only if all finite eigenvalues have magnitude less than one.
Proposition 2. The stability of (16) can be assessed by studying the stability of a linear discrete-time system in the form of (12), where: with C † having dimensions r × c ν r. The proof of Proposition 2 is provided in the Appendix.

C. Eigenvalue Mapping and Validity of SSSA
In both single-delay and multiple-delay cases, the eigenvalues of system (12) represent, in the Z-plane, the numerically distorted by the TDI method dynamic modes of system (3). Let z k be the eigenvalue of (12) that approximates the k-th mode of the DDAE power system model. The latter is represented by the eigenvalue s k = α+ȷβ, which is a root of the characteristic equation (7). Then, the actual and distorted modes s k ,ẑ k , can become directly comparable by mapping the one to the domain of the other. Mappingẑ k from the Z to the S plane, we have: where log(·) denotes the complex logarithm. Then, the numerical distortion caused to the k-th mode by the TDI method is: while the distortion caused to the damping of this mode is: where ζ k = −α/(α 2 + β 2 ). Positive values of d ζ,k in (21) indicate that the mode is overdamped, whereas negative values indicate that the mode is underdamped. Equations (20) and (21) are based on SSSA and thus they are in principle valid only at a steady state solution of the system. Around the steady state solution, these equations are accurate measures of the distortion of the system's modes given h, or vice versa, can be employed to determine the value of h required to achieve a required level of precision. These measures are also good estimates of the approximation introduced by TDI methods under varying operating conditions, owing first, to that the structure of the dynamic modes and the stiffness of a power system model are features that do not alter dramatically by changing the operating point, and second, that TDI methods maintain certain qualitative properties, e.g. numerical stability properties, when applied to different systems of the same class. This suggests that, for a given network, the proposed analysis does not need to be repeated often and, possibly, can be done only once. Similar considerations that support the statement above are given, e.g., in [18]- [20].

D. Examples
In this section we consider three examples of implicit RK methods applied for the TDI of (3), namely the BEM, the ITM, and the 2-stage Radau IIA. The same methods are then utilized in the simulations of Section IV. The goal here is to show how the linearized version of each method can be formulated so that Proposition 2 can be readily applied. For generality, we assume that the system includes ν delays. The single delay case can be retrieved by substituting ν = 1.
Backward Euler Method: The Butcher tableau [21] of the BEM is: Applied for the TDI of (3), the method takes the form: Linearization of (23) gives: or equivalently: The last system is in the form of (11), where: Note that including the delayed Jacobians A k in (24) is necessary to capture the delay effects on the system. This should not be confused with the structure of the Jacobian used in the iterative solution of each point of the integration, e.g. in Newton's method, where delays can be viewed as known constants and thus omitted, see [2]. Implicit Trapezoidal Method: The Butcher tableau of the ITM is: 0 0 0 1 0.5 0.5 0.5 0.5 Applying the method for the TDI of (3) and linearizing, gives: System (28) can be written in the form of (11), where: 2-Stage Radau IIA: The Butcher tableau of the 2-stage Radau IIA method is: 1/3 5/12 −1/12 Applying the method for the TDI of (3) and linearizing, gives: Solving (31) for u t−h+h/3 and substituting in (32) yields: where M = (E − 5h 12 A 0 ) −1 . System (34) can be rewritten in the form of (11), where:

IV. CASE STUDY
This section presents simulation results based on the IEEE 14-bus system. The system consists of 14 buses, 5 synchronous machines equipped with Automatic Voltage Regulators (AVRs), 12 transmission lines, 4 transformers, and 12 loads. Moreover, the machine connected to bus 1 is equipped with a Power System Stabilizer (PSS). The DAE system model has in total 52 state and 92 algebraic variables. The full static and dynamic data are taken from [10] and the gain of the AVR of the machine at bus 1 is reduced by 2.5 times compared to these data, a modification that secures for the system an adequate delay margin to facilitate a full comparison among the examined TDI methods for both small and large delay and time step sizes. Simulations in this section are executed using Dome [22] and eigenvalues are computed with LAPACK [23].
Let us consider first the system without any delay. The eigenvalue analysis shows that the system is stable and that the most poorly damped mode is the complex pair −0.976 ± ȷ5.970. Hereafter we will refer to this mode as Mode 1. Figure 1 shows, for different integration time step sizes, how the dynamic modes of the DAE power system model are approximated by the three TDI methods discussed in Section III-D. The plots are drawn by first solving for each method equation (8) and then mapping the computed roots to the S-plane according to (19), so that they are directly comparable to the eigenvalues of the DAE system, i.e. the roots of (8), where, for the delay-free case, one has A k = 0 r,r .
Results indicate that the BEM overdamps the system's dynamics (see Fig. 1a), which is as expected, since the method is known to be hyperstable, see also the relevant discussion in [20]. Moreover, the ITM, while being very accurate for time steps in the order of 10 −2 s, it introduces underdamped oscillations when employed with large time steps (see Fig. 1b). Finally, under the same time step, the most accurate among the examined methods is the 2-stage Radau IIA method. We now introduce a delay. With this aim, we assume that the input signal of the PSS connected to the machine at bus 1 is impacted by a constant delay τ . By varying the delay magnitude and repeating the eigenvalue analysis we find that the delay margin of the system is 0.39 s. In fact, for τ = 0.4 s, the system has one unstable mode 0.007 ± ȷ7.498, which is thus the most critical for the delay-dependent stability of the system. Hereafter we will refer to this mode as Mode 2.
Note that, in the delay-free system, Mode 2 is represented by the complex pair −2.719 ± ȷ8.898. The eigenvalues of the DDAE system are determined considering Chebyshev's discretization technique (see also Section II-B) with 8 nodes, which is a choice that provides a good compromise between precision and computational burden.     Fig 2b. The magnitude of the overdamping per se in this case is small, yet, interestingly, the response of the method is in the opposite direction from what one would expect, that would be, similar to the delay-free case, all modes to be more or less underdamped.
In both cases illustrated in Fig. 2, the BEM overdamps the system's oscillations. This property of the BEM is well known. On the other hand, what is of interest in this work and not well known is how the numerical approximation is impacted due to the presence of the delay. To this aim, we carry out a number of simulations considering a varying delay magnitude and a constant time step, and we compute for each scenario and method the numerical and damping distortion according to (20) and (21). Results are presented for Mode 1 and Mode 2.    Tables I and II. In this case, the Radau IIA method is the most accurate among the three methods for all delays considered. Furthermore, for large delays, the overdamping introduced by the BEM decreases, which is also consistent with Figs. 3b and 4b. Finally, as the delay increases, the underdamping introduced by the ITM to Mode 1 increases, whereas for Mode 2, it decreases, so that for delays larger than 0.2 s the mode appears overdamped. This is in line with the discussion of Fig. 2 provided above.
We consider that the delay is fixed at 0.2 s, and we examine how the distortion of the two modes varies as the time step h is increased. The results are presented in Fig. 5, and show that, as expected, the magnitude of d s for the two modes for all three methods increases with the increase of h. The different response of the ITM in capturing the damping of the two modes is confirmed also in this figure. Furthermore, Fig. 5f indicates that the Radau IIA may induce underdamping or overdamping depending on the step size (see Mode 1). A relevant remark is that, assuming only constant delays that are integer multiples of h limits our ability to draw plots like Fig. 5 with full accuracy. The assumption is made in this work for the sake of simplicity, e.g. we have no need for explicit calls on interpolation in the definition of matrices F and G. The task of extending the formulation for delays of any value, including time-varying and stochastic delays, is left as a task for future work.

V. CONCLUSIONS
The paper presents a framework to evaluate the numerical approximation that TDI methods introduce when employed for the numerical integration of power system models impacted by time delays. The proposed framework is based on SSSA and is formulated in a general way that covers many methods, including the most important family of numerical methods for the integration of DDEs, i.e. implicit RK methods. We will 22nd Power Systems Computation Conference dedicate future work to evaluate the computational burden of our approach when applied to large-scale systems, as well as to extend it for time-varying and stochastic delays, plus for delays that are not integer multiples of the time step.