1 Differential algebraic equation (DAE)

A simple description of DAE is ordinary differential equation with algebraic equation. Different ways to look at DAE can tell us how DAE can be used. A system of DAE can be formed by a model of a system (ODE), which we want to investigate, have some algebraic. constrains,

\begin{aligned}x′ & = & f(x, y, t) \\ y & = & g(x, t) \end{aligned}

. This type of systems is called Constrained Variational problems. Transient stability problem is of this type. DAE can also come from model reduction and singular perturbations. Another view of the network equation can classify Transient stability problem as model reduction. Model reduction is a system of ODE of which some variables contribute a very small factor to the system. When we are interested in the first order approximation of its behavior, we eliminate the some of the variables.

\begin{aligned}x′ & = & f(x, y, t, ϵ) \\ ϵy′ & = & g(x, y, t, ϵ) \end{aligned}

, where ϵ → 0. The resulted system is a DAE system. The difference between DAE and ODE is that DAE can be more difficult to solve. The degree of difficulty can be indicated by the index of the DAE system. A system of ODE is classify as index 0. Next, we have index 1 DAEs. These DAE can subdivided into two classes, semi-explicit

\begin{aligned}x′ & = & f(x, y, t) \\ 0 & = & g(x, y, t) & ∃\frac{∂g}{∂y} \end{aligned} ,

and full-implicit system,

\begin{aligned}0 & = & f(x, x′, y, t) \\ 0 & = & g(x, y, t) & ∃\frac{∂g}{∂y} \end{aligned}

. There are higher index problem; however, we will not introduce them in this report. The underlining model of electrical transient analysis is based on the simple RLC circuit model. The transmittion line, transformer and capator is modeled by a perfect components, such as resistor, capticler and inductor. The infinite voltage and current sources is replace by the classical generator model and the 2-axis generator model. A classical generator model is a perfect voltage source with bookkeeping for the rotor angle added. The 2-axies generator model is more detail model of a generator. As a result, the rotor angle keep track by a 2-axies generator model is more accurate. The two type of generators contribute ODEs to the model’s equations. The network is contribute the algebraic. system of equation. The model’s equations are

\begin{aligned}x′ & = & f(x, V, I) \\ I & = & Y (x)V \end{aligned}

. As a result, electrical transient problem consists of solving semi-explicit index one differential algebraic equations.

1.1 Integration methods

In PTI software, modified Euler method is used. There are two types of modified Euler methods, explicit and implicit. Explicit modified Euler method is not stiff accurate. In other words, the explicit modified Euler method requires small step size to maintain the desire accuracy when solving stiff systems [9]. A stiff system is generally referred to a system, which some of the variables varies slowly and some variables varies very fast. In our case, the generator’s rotor angle varies slowly (governed by the ordinary differential equations), while the voltage of the network varies instantaneously (governed by the algebraic equations). Since the explicit modified Euler method is not stiff accurate, the numerical solution of the DAE system do not automatically satisfy the algebraic part of the DAE, ie, I_{n}−Y (x_{n})V ≠ 0. The formulation of implicit modified Euler method is the same as the trapezoidal rule. Implicit modified Euler method, on the other hand, have a very different properties. In ETMSP, electrical transient analysis problem is currently solved by applying trapezoidal rule to the system [7]. The system of equation becomes,

[ \begin{matrix}F \\ G \end{matrix} ] = [ \begin{matrix}x_{n+1} − x_{n} − \frac{h}{2}(f(x_{n+1}, V_{n+1}) + f(x_{n}, V_{n})) \\ Y V_{n+1} − I(x_{n+1}, V_{n+1}) \end{matrix} ] .  (1)

(2)

Trapezoidal rule is a stuff accurate, order 2 method [8]. A stuff accurate method implies that the step size used on solving the problem have reduced impact on the accuracy of the numerical result. The order 2 accurate means the order convergence is proportional to h^{2}. In addition, the trapezoidal rule is A-stable, so that the behavior of the system model will not be affected by the analytical behavior of the numerical method. Since trapezoidal rule is stuff accurate, the algebraic part of the DAE system is satisfied [2]. Consequently, the numerical result of this method is more accurate and the method allows to use large step size. Furthermore, new numerical model can be developed without any concern on the analytical behavior of the numerical DAE solver’s method. As a result, ETMSP is a good research tools on developing electrical component’s models. The implicit non-linear algebraic equations 1 are solved by Newton method. After applying Newton method to the system, the system becomes

[ \begin{matrix}F \\ G \end{matrix} ] = [ \begin{matrix}A & B \\ C & D \end{matrix} ] [ \begin{matrix}△x \\ △V \end{matrix} ] ,  (3)

(4)

where A is the Jacobian matrix for the generator equations, B is the Jacobian matrix for the coupling between generators, C is the Jacobian matrix for the coupling between generators and D is the Jacobian matrix for the network equations. The system then further decoupled into two systems,

x_{n+1}^{k+1} = x_{n+1}^{k} − A^{−1}(F_{n+1}^{k} + B(V_{n+1}^{k+1} − V_{n+1}^{k}) (5)
(Y + D − CA^{−1}B)(V_{n+1}^{k+1} − V_{n+1}^{k}) = −G_{n+1}^{k} + CA^{−1}F_{n+1}^{k}  (6)

(7)

The first system consists of all the generator equations. The sparcity pattern of A is block diagonal. Each diagonal block represents ODE of the generator’s model. As a result, the generator equation is mutually independent of other generator. The nonzero values of the second matrix contributed mostly by the network equations, Y; however, additional nonzero values may be added by D − CA^{−1}B, which is introduced by the decoupling process. The advantages of this decoupling reduce the over all size of the matrices that need to be factorized. This reduction can be show by the following calculation. Let a be the dimension of the Jacobian matrix of the generator equations, A and d be the dimension of the Jacobina matrix of the network equations, D.

\begin{aligned}(a + d)^{3} & = & (a^{2} + 2ad + d^{2})(a + d) \\ & = & a^{3} + 3a^{2}d + 3ad^{2} + d^{3} \end{aligned} This computation assumes that the matrices, A and

D, are dense. The reduction is discounted by using sparse matrix algorithms. However, this reduction can still be significant for sequential processing. The trade off for this processing method is the destruction of the sparcity pattern offered by the system.

1.2 DAE solver

As we see from the previous section, the trapezoidal rule is a excellent tool for developing models for new components. However, once the component is analyzed and determined stable. Trapezoidal rule becomes the limiting factor of the computation speed. In this section, we will introduce a new DAE solver package to replace trapezoidal rule. DASSL [11] [12] is specified designed for solving index one DAEs. The features of this packages are variable convergence order, variable step size and clean interface to linear sparse solver. DASSL [2], the DAE solver package, uses implicit backward difference formula. BDF provides formula for convergence from order 1 to 5. The analytical stability of order 1 and 2 BDF formula are A-stable. In fact, order 1 BDF is Euler method and order 2 BDF is modify Euler method. However, the gain in computational speed are expected form using the higher order of convergence formula. The solver package have automatic control for shifting from the lower formula to higher order formula, when the function exhibit stable behavior for higher order formula. At higher convergence order, the error of the result is reduced. Hence, larger step size can be used to compute the result with the same accuracy. However, the trade-off is that the higher order formula is not A-stable. The stability of

these higher formula is A(α)-stable [8]. \begin{matrix}order & 1 & 2 & 3 & 4 & 5 \\ α & 90 & 90 & 86.03 & 73.35 & 51.84 \end{matrix} As a re-

sult, using DASSL for model development requires in-depth knowledge of DASSL and the working knowledge of the model. However, if the model’s behavior is known, then DASSL can be used to provide fast simulations. DASSL [2] package also have a automatic variable step size control. Although ETMSP only allows fixed step size, its users usually use different step size for different phase of the simulation. For example, the user uses larger step size before fault, smaller step size just after the fault and larger step size after the disturbance is expected to be settled fault. In DASSL, the selection of the step size is based on the acceptable error, which is defined by the user, and the behavior of the function. If the history of the function indicates the function is “slow” varying, then the step size is increase until either the error is larger than specified or unstable result is computed. As a result, we expected further reduction of the runtime of a simulation. Furthermore, the DASSL package have a well defined interface to linear system solver [10]. This interface allows DASSL to call any parallel sparse matrix solver with minimum code changes. Currently, DASSL can use two major types of linear solver, direct solver and semi-iterative solver.

2 Sparse Solver

2.1 ETMSP

ETMSP uses a inhouse developed sparse matrix factorization package. In this section, we briefly introduce facilities of this sparse matrix package. The sparse package provides five major functions [1],

• Matrix reordering

• Symbolic factorization

• LU factorization

• Forward and Back solve

• Interface and Supporting functions

The matrix reordering strategy consists of static minimum degree ordering. The matrix is reordered once before the first factorization. The reordering is redone only, when there is significant changes to the matrix structure. Minimum degree ordering [6] is referred to placing the bus with minimum number of off diagonal elements before the bus other buses. This step is repeated for the n − 1 × n − 1 submatrix, where n is the dimension of the original matrix. The symbolic factorization can be used to discover the fill-ins patterns of the matrix during the factorization with specific ordering. Symbolic factorization is computed in a similar way to a factorizing a matrix; however, no numerical computation is done. As a result, symbolic factorization is much cheaper than a full factorization. The symbolic factorization function in the ETMSP sparse solver packages is not used. The LU factorization function of this package offers three functionalities. First, the LU factorization perform a full LU factorization on a symmetric structured matrix with ordering given by the matrix reordering function. Second, the function can partially un-factorization the matrix. This function is intend to provide a way to reduce the cost of factorization, when the changes of the matrix occur at the later half of the matrix. Third, partial factorize the previously partial un-factorization matrix. However, ETMSP only use the function for full LU factorization. Forward and Back solve functions is used in the very dishonest newton function. The interface function provides interface services for ETMSP to access the sparse package. Additional function to store and to restore the solver’s state is provided.

3 Semi-iterative solver

In the course of solving DAEs, DASSL [11] [12] uses the Newton’s method to solve systems of non-linear equations. Each Newton iteration requires to solve a system of linear equations, Ax = b. In the original DASSL [11], LU factorization is employed to solve the system of linear equations. However, the recent advances on Krylov subspace iterative linear solvers, for example, CGS, QMR and GMRES, has prompted the DASSL’s authors to incorporate GMRES iterative solver into DASSL [12]. In [4], [3] and [5], the authors show DASSL and Krylov subspace iterative solver combination, the sequential software package is renamed to DASPK, offers reduced memory storage for large problems. Experience shows that DASSL

A = [ \begin{matrix}3 & 1 & 0 & 0 \\ 1 & 3 & 1 & 0 \\ 0 & 1 & 3 & 1 \\ 0 & 0 & 1 & 3 \end{matrix} ]  C = [ \begin{matrix}\frac{1}{3} & 0 & 0 & 0 \\ 0 & \frac{1}{3} & 0 & 0 \\ 0 & 0 & \frac{1}{3} & 0 \\ 0 & 0 & 0 & \frac{1}{3} \end{matrix} ]

Λ_{A} = [ \begin{matrix}3.6180 & 0 & 0 & 0 \\ 0 & 4.6180 & 0 & 0 \\ 0 & 0 & 2.3820 & 0 \\ 0 & 0 & 0 & 1.3820 \end{matrix} ]

Λ_{CA} = [ \begin{matrix}1.2060 & 0 & 0 & 0 \\ 0 & 1.5393 & 0 & 0 \\ 0 & 0 & 0.7940 & 0 \\ 0 & 0 & 0 & 0.4607 \end{matrix} ]

Figure 1: The eigenvalues of A is clustered after the application of the precondition matrix C.

with semi-iterative solver computes the solution faster than DASSL with LU factorization. The test program is solving heat equations, so the timing of the test program has no bearing on predicting the DASSL’s performance on solving the electrical power problem. However, the test does demonstrate a possible way for exploration of fast on-line transient stability program. Currently, Concurrent DASPK is capable to use three iterative solvers, CGS, TFQMR and GMRES. When computing a Newton iteration, DASPK does not compute the Jacobian matrix (J). Instead, an approximation vector of Jy is generated by using a difference quotient [5],

Jy = F′(x)y ≈ \frac{F (x + σ) − F (x)}{σ}, σ ∈ ℜ.

Using this approximation, DASSL has shown to converge theoretically and numerically, when DASSL is applied to a system of ODEs. Although no difficulties are expected from using DASPK on the electrical power problem, a closer study of the interaction of DASPK with the electrical power DAEs is needed. In order to further speed up the convergence of the semi-iterative solving Ax = b, a preconditioner is used. The goal of preconditioner is to generate a precondition matrix (C) as close to the inverse of Jacobian matrix as possible with acceptable amount of computations. The basic function of the precondition matrix can be illustrate by the equation,

CJx = Cb.

If C = J^{−1}, then the system of linear equations is solved. However, J^{−1} is very difficult to generate, so C is a good approximation of J^{−1} at best. In this situation, the precondition matrix serves as a mechanism to cluster the eigenvalues of J as shown in figure 1. Although C ≈ J^{−1} does not solve the system, it maps the system of linear equations to a better behaved system in the point of view of the Krylov subspace methods. The underline process of the Krylov subspace methods is to find a function (f) such that f(0) = 1 and the roots of f are the eigenvalues of the matrix J. The function is found by increasing the index of the current function by one for each iteration as shown in figure 2. Another important property of iterative methods is that the large eigenvalues are found with relative ease. When a iterative solver is terminated at the ith iteration and i is less Figure 2: The minimal polynomial plot shows a polynomial with the eigenvalues as roots. The degree i polynomial shows the shape of the current polynomial at i iteration of the linear solver.

than dim J, the user accepts the function, f_{i}, as a “good” approximation to the function, f_{n}. Consequently, the clustered eigenvalues allows a lower order function to approximate the higher order function better, which translates into fast convergence of the linear solver.

References

[1] Extended transient-midterm stability program: Version 3.0.

[2] K.E. Brenan, S.L. Campbel, and L.R. Petzold. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. Elevier Science, 1989.

[3] Peter N. Brown. A local convergence theory for combined inexact-newton/finite-difference projection methods. SLAM Journal on Numerical Analysis, 1987.

[4] Peter N. Brown and Alan C. Hindmarsh. Matrix-free methods for stiff systems of ode’s. SLAM Journal on Numerical Analysis, 1986.

[5] Peter N. Brown and Alan C. Hindmarsh. Reduced storage matrix methods in stiff ode systems. Technical report, Numerical mathematics Group, Lawrence Livermore Natrional Laboratory, 1989.

[6] I. S. Duff, A. M. Erisman, and J. K. Reid. Direct Methods for sparse Matrices. Oxford Science, 1990.

[7] Electric Power Research Institute. Extended Transient-Midterm Stability Program: Version 3.0 Volumn 4: programmer’s Manual, Part I.

[8] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II. Springer-Verlag, 1980.

[9] J. H. Hubbard and B. H. West. Differential Equation: Part I. Springer-Verlag, 1990.

[10] Alvin Leung, Anthony Skjellum, and Geoffrey Fox. Concurrent dassl: A secondgeneration dae solver library. Scalable High Performance Computing Conference, Oct 1993.

[11] Linda R. Petzold, Peter N. Brown, and Alan C. Hindmarsh. ddassl.for - differential/algebraic solver, 1987.

[12] Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and Clement W. Ulrich. ddaspk.for - differential/algebraic solver, Sept 1990.