TIME-STACKING METHOD FOR DYNAMIC SIMULATIONS

Abstract

A time-stacking method is disclosed. The time-stacking method simulates dynamics so as to obtain time-domain trajectories which correspond to a disturbance or event. The method includes providing a model for the system. The model includes differential equations and algebraic equations. The method also includes solving the differential equations and the algebraic equations over a predetermined number of time steps simultaneously using an implicit integration scheme.

Description

TECHNICAL FIELD

This invention relates to dynamic simulations. More specifically, this invention relates to a time-stacking method of simulating dynamics of a system by solving multiple time steps simultaneously using an implicit integration scheme.

BACKGROUND

Dynamic simulation plays a significant role in secure operations and planning of systems. In power networks, dynamic simulations solve a large set of differential algebraic equations for thousands of times that describe the electro-mechanical interaction of generators as well as power electronic devices and their controllers. The solutions of these equations determine the time-series trajectory of a power grid when subject to disturbances such as short-circuit faults, generator tripping, line switching or other severe contingencies.

The dynamics of a power system can be represented by a set of first-order differential equations (1a) and algebraic equations (1b):

$\begin{array}{cc}\{\begin{array}{c}\stackrel{.}{x}=f\left(x,y\right)\\ 0=g\left(x,y\right)\left(1b\right)\end{array}& \left(1a\right)\end{array}$

Where x represents a vector of state variables describing generator or load dynamics (e.g., generator rotor angle, rotor speed, sub-transient voltages or flux) and y represents a vector of algebraic variables (e.g., voltages from all or selected buses in the network and the current injections from generators or loads, which can be expressed as magnitudes and phase angles, or real and imaginary parts).

A major obstacle of introducing dynamic simulation into online security assessment is the computationally intensive time-domain solution of numerous differential and algebraic equations involved in modeling the system dynamics and network. The traditional method for dynamic simulation is serial in nature and can only solve the network dynamics one step after another. This process is very slow for a larger-scale power network model, such as the Western or Eastern Interconnection in North America.

A common serial dynamic simulation method is the partitioned or alternating solution. The partitioned-solution method is predominately used in industrial systems or commercial simulation software programs. In this method, the differential equation set of (1a) is discretized and then solved for x_γ (x at time step γ). The algebraic equation set (1b) is then solved for y_γ.

The partitioned-solution scheme includes an explicit integration method and an implicit integration for solving the equations, as shown in FIG. 1.

A. Partitioned Method with Explicit Integration

Using an explicit integration method, equation 1(a) is converted into an algebraic equation set

x _γ =F(x_γ−1,y_γ−1).  (2)

which can be used to solve for x_γ using the previous step values. The value of x_γ is used to solve (3) for y_γ.

0=g(x_γ,y_γ).  (3)

B. Partitioned Method with Implicit Integration

Using an implicit integration method, equation (1a) is converted into an algebraic equation set

x _γ =F(x_γ,y_γ).  (4)

The initial guess of x_γ can be obtained based on x_γ−1 and y_γ−1 using any explicit integration formula. The value of y_γ is then used in equation (4) to obtain x_γ, which is used in equation (3) to solve y_γ. This process will be repeated until the solution is converged.

SUMMARY

The present invention describes methods of solving multiple time steps of system dynamics simultaneously so that the total computation time may be significantly reduced using parallel computing techniques. The method uses a single set of equations that models multiple time intervals simultaneously, which includes applying high-performance computing to solve the equations simultaneously in parallel.

In one embodiment of the present invention, a time-stacking method is disclosed. The time-stacking method simulates dynamics of a system to obtain time-domain trajectories of the system that correspond to a disturbance or event. The method includes a model for testing the system. The model includes differential equations and algebraic equations. The method also includes solving the differential equations and the algebraic equations over a predetermined number of time steps simultaneously using an implicit integration scheme. The method may also include applying a numerical solver to compute a final result of the equations. A supercomputer or high-performance computer may be applied for solving the equations simultaneously in parallel.

The implicit integration scheme includes, but is not limited to, at least one of the following: Backward Euler, Trapezoidal, Runge-Kutta, Newton, and Dishonest Newton.

The differential equations may represent vectors of state variables, and the algebraic equations may represent vectors of algebraic variables.

The vectors of state variables include, but are not limited to, generator and its controller dynamics. The generator dynamics may include rotor angle and generator speed.

The vectors of algebraic variables include, but are not limited to, voltages from all or selected buses in the network and the current injections from generators or loads, which can be expressed as magnitudes and phase angles, or real and imaginary parts.

In one embodiment, the time steps comprise intervals of at least 1 microsecond. The time steps of intervals may vary over a range of at least 1 microsecond so as to obtain a calculation that corresponds to the disturbance or event.

In another embodiment of the present invention, a method of simulating power system dynamics of a system so as to obtain time-domain trajectories which correspond to a disturbance or event is disclosed. The method includes providing a model for the system. The model includes differential equations and algebraic equations. The method further includes solving the differential equations and the algebraic equations over a predetermined number of time steps simultaneously using an implicit integration scheme. The method also includes applying a numerical solver to compute a final result of the equations, and applying a supercomputer or high-performance computer for solving the equations simultaneously.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart of a prior art partitioned-solution method using implicit or explicit methods, for simulating dynamics of a system.

FIG. 2 is a flow chart of a time-stacking method of simulating dynamics of a system so as to obtain time-domain trajectories which correspond to a disturbance or event, in accordance with one embodiment of the present invention.

FIG. 3A shows a Butcher tableau for an explicit Runge-Kutta integration method.

FIG. 3B shows a Butcher tableau for an implicit Runge-Kutta integration method.

FIG. 4 shows a Butcher tableau for a classical explicit Runge-Kutta integration method.

FIGS. 5A-5C shows stability domains for (A) Forward Euler, (B) Backward Euler, and (C) Trapezoidal integration methods.

FIG. 6 is a table showing the variables used in the code and the corresponding notation for a model for the system with a modified Euler integration method.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention relates to time-stacking methods of solving multiple time steps of system dynamic performance simultaneously by using an implicit integration scheme. The present invention solves a single set of differential and algebraic equations simultaneously over a predetermined number of time steps. The parallel time-stacking methods solve multiple time steps simultaneously. As a result, the time-stacking methods of the present invention significantly speed up simulations of system performance.

In the simultaneous time-stacking method of the present invention, as shown in FIG. 2, equation (1a) is converted into a set of algebraic equations, which are lumped together with equation (1b) to form a single larger algebraic equations set. Variables X_{γ and} y_γ are solved simultaneously over a predetermined number of times steps using an implicit integration scheme. As an example, writing equation (4) in the form of

0=H(x_γ,y_γ),  (5)

and then combining equation (5) and equation (3), a single set of algebraic equations can be obtained. As one example, Newton's method can be used to solve the algebraic equation set, requiring the construction and solution at each iteration of the Jacobian matrix

$\begin{array}{cc}J=\left[\begin{array}{cc}{J}_{\mathrm{Fx}}& J_{\mathrm{Fy}}\\ J_{\mathrm{gx}}& J_{\mathrm{gy}}\end{array}\right].& \left(6\right)\end{array}$

No matter which scheme is used, those non-linear algebraic equations need to be solved using, for example, Gauss or Newton type procedures. Dishonest Newton or very dishonest Newton methods are often used in industrial software to increase computation speed. When Newton-like method is used, linear solvers are included as well.

Mathematic Formulation of Multi-Machine System with Classical Model

When generators are represented using the classical model, the equations of motion for generator i in per unit are:

$\begin{array}{cc}\{\begin{array}{c}{\stackrel{.}{\delta }}_{\left(i\right)}={\omega }_b{\omega }_{\left(i\right)}\\ {\stackrel{.}{\omega }}_{\left(i\right)}=\frac{1}{2H_{\left(i\right)}}\left(P_{m\left(i\right)}-P_{e\left(i\right)}-D_{\left(i\right)}{\omega }_{\left(i\right)}\right)\left(7b\right)\end{array}& \left(7a\right)\end{array}$

where ω_(i) is the per unit speed deviation for generator i, δ_(i) is the angular position of the rotor of generator i in electrical radians with respect to a synchronously rotating reference, H_(i) is the inertia constant of generator i using system base VA², D_(i) is the damping factor or coefficient of generator i in pu torque/pu speed deviation, P_m(i) is the mechanical power input of generator i, and ω_b is the base rotor electrical speed in radians per second. In this case, state variables in (1) only include rotor speed and angle, i.e.,

x=[δ ₍₁₎,ω₍₁₎, . . . ,δ_(i),ω_(i), . . . ,δ_(n),ω_(n)]^T,  (8)

where n is the number of generators.

Using the classical model to study the transient stability for a multi-machine system, it is often assumed:

1) Mechanical power input (P_m(i)) is constant.

2) The mechanical rotor angle of a machine coincides with the angle of the voltage behind the transient reactance.

3) The network is assumed to be in the sinusoidal stead state.

4) Load is represented by the passive impedances, i.e., the dynamics of load is ignored.

Therefore, let Y″ (with a dimension equal to m by m) denote the nodal admittance matrix of an m-bus system comprised of n generator buses and m−n load buses. One can add machine internal buses and include load impedance into the admittance matrix, resulting in an extended Y-bus—Y′ (with a dimension equal to m+n by m+n). The network equation becomes:

$\begin{array}{cc}\left[\begin{array}{c}{I}_E\\ 0\end{array}\right]=Y^{\prime }\left[\begin{array}{c}E\\ V\end{array}\right].& \left(9\right)\end{array}$

The network equation can be reduced to

I _E =YE,  (10)

where Y=Y_nn−Y_nm Y⁻¹_mmY_mn. Numerical evaluation of direct solvers for large sparse, symmetric linear equations could be used to solve Y⁻¹_mmY_mn. For each generator internal voltage bus i, the injection current in system reference can be expressed as

$\begin{array}{cc}{\stackrel{\_}{I}}_{\left(i\right)}^{\mathrm{sys}}\stackrel{\Delta }{=}I_{\mathrm{re}\left(i\right)}+jI_{\mathrm{im}\left(i\right)}=\sum _{k=1}^n{\stackrel{\_}{Y}}_{\left(\mathrm{ik}\right)}{\stackrel{\_}{E}}_{\left(k\right)}& \left(11\right)\end{array}$

where a bar above the notation represents a complex number, Ī^sys_(i) is the injection current of bus i in system reference, and Ē(k) is the stator internal voltage of bus k in system reference. Let Ē(k)=E_re(k)+jE_im(k) and Y_(ik)=G_(ik)+jB_(ik), resolving (11) into real and imaginary part yields,

$\begin{array}{cc}{I}_{\mathrm{re}\left(i\right)}=\sum _{k=1}^n\left[E_{\mathrm{re}\left(k\right)}G_{\left(\mathrm{ik}\right)}-E_{\mathrm{im}\left(k\right)}B_{\left(\mathrm{ik}\right)}\right],& \left(12a\right)\\ I_{\mathrm{im}\left(i\right)}=\sum _{k=1}^n\left[E_{\mathrm{re}\left(k\right)}B_{\left(\mathrm{ik}\right)}-E_{\mathrm{im}\left(k\right)}G_{\left(\mathrm{ik}\right)}\right].& \left(12b\right)\end{array}$

The transformation from machine dq0 reference frame to system reference frame is:

$\begin{array}{cc}\left[\begin{array}{c}\mathrm{re}\\ \mathrm{im}\end{array}\right]=\left[\begin{array}{cc}\sin \delta & \cos \delta \\ -\cos \delta & \sin \delta \end{array}\right]\left[\begin{array}{c}d\\ q\end{array}\right],& \left(13\right)\end{array}$

In classic model, since E′_d=0, applying the above transformation to I_re, I_im, E_re, and E_im in equation (12) and after manipulations, the following equations are obtained in machine reference frame:

$\begin{array}{cc}0=\sin {\delta }_{\left(i\right)}I_{d\left(i\right)}+\cos {\delta }_{\left(i\right)}I_{q\left(i\right)}-\sum _{k=1}^n\left[\begin{array}{c}\left(\cos {\delta }_{\left(k\right)}E_{q\left(k\right)}^{\prime }\right)G_{\left(\mathrm{ik}\right)}-\\ \left(\sin {\delta }_{\left(k\right)}E_{q\left(k\right)}^{\prime }\right)B_{\left(\mathrm{ik}\right)}\end{array}\right],& \left(14a\right)\\ 0=-\cos {\delta }_{\left(i\right)}I_{d\left(i\right)}+\sin {\delta }_{\left(i\right)}I_{q\left(i\right)}-\sum _{k=1}^n\left[\begin{array}{c}\left(\cos {\delta }_{\left(k\right)}E_{q\left(k\right)}^{\prime }\right)B_{\left(\mathrm{ik}\right)}+\\ \left(\sin {\delta }_{\left(k\right)}E_{q\left(k\right)}^{\prime }\right)G_{\left(\mathrm{ik}\right)}\end{array}\right].& \left(14b\right)\end{array}$

The active power at air gap can be expressed as:

P _e(i) =E′ _q(i) I _q(i).  (15)

Replacing P_e(i) in equation (7) by equation (15) yields:

$\begin{array}{cc}\{\begin{array}{c}{\stackrel{.}{\delta }}_{\left(i\right)}={\omega }_b{\omega }_{\left(i\right)}\\ {\stackrel{.}{\omega }}_{\left(i\right)}=\frac{1}{2H_{\left(i\right)}}\left(P_{m\left(i\right)}-E_{q\left(i\right)}^{\prime }I_{q\left(i\right)}-D_{\left(i\right)}{\omega }_{\left(i\right)}\right)\left(16b\right)\end{array}& \left(16a\right)\end{array}$

Combining equation (16) and equation (14) results in a set of differential-algebraic equations (DAE) in the same form as equation (1), where:

x=[δ ₍₁₎,ω₍₁₎, . . . ,δ_(i),ω_(i), . . . ,δ_(n),ω_(n)]^T,  (17)

y=[I _d(1) ,I _q(1) , . . . ,I _d(i) ,I _q(i) , . . . ,I _d(n) ,I _q(n)]^T,  (18)

and the other notations represent parameters.

Numerical Integration Using Trapezoidal Method

The trapezoidal integration method is used in many time-domain simulation software and can be used to solve the classical model multi-machine system described above.

For time step γ, applying Trapezoidal rule to equation (16) and then combining the resulted algebraic equations with equation (14), yields:

$\hspace{1em}\begin{array}{cc}\{\begin{array}{c}0={\delta }_{\left(i,\gamma \right)}-{\delta }_{\left(i,\gamma -1\right)}-\frac{h}{2}\left({\omega }_b{\omega }_{\left(i,\gamma -1\right)}+{\omega }_b{\omega }_{\left(i,\gamma \right)}\right)\\ 0={\omega }_{\left(i,\gamma \right)}-{\omega }_{\left(i,\gamma -1\right)}-\frac{h}{2}\left[\begin{array}{c}\frac{1}{2H_{\left(i\right)}}\left(\begin{array}{c}{P}_{m\left(i,\gamma -1\right)}-E_{q\left(i\right)}^{\prime }I_{q\left(i,\gamma -1\right)}-\\ D_{\left(i\right)}{\omega }_{\left(i,\gamma -1\right)}\end{array}\right)+\\ \frac{1}{2H_{\left(i\right)}}\left(P_{m\left(i,\gamma \right)}-E_{q\left(i\right)}^{\prime }I_{q\left(i,\gamma \right)}-D_{\left(i\right)}{\omega }_{\left(i,\gamma \right)}\right)\end{array}\right]\left(19b\right)\\ 0=\sin {\delta }_{\left(i,\gamma \right)}I_{d\left(i,\gamma \right)}+\cos {\delta }_{\left(i,\gamma \right)}I_{q\left(i,\gamma \right)}-\sum _{k=1}^n\left[\begin{array}{c}\left(\cos {\delta }_{\left(k,\gamma \right)}E_{q\left(k\right)}^{\prime }\right)G_{\left(i,k\right)}-\\ \left(\sin {\delta }_{\left(k,\gamma \right)}E_{q\left(k\right)}^{\prime }\right)B_{\left(i,k\right)}\end{array}\right]\left(19c\right)\\ 0=-\cos {\delta }_{\left(i,\gamma \right)}I_{d\left(i,\gamma \right)}+\sin {\delta }_{\left(i,\gamma \right)}I_{q\left(i,\gamma \right)}-\sum _{k=1}^n\left[\begin{array}{c}\left(\cos {\delta }_{\left(k,\gamma \right)}E_{q\left(k\right)}^{\prime }\right)B_{\left(i,k\right)}+\\ \left(\sin {\delta }_{\left(k,\gamma \right)}E_{q\left(k\right)}^{\prime }\right)G_{\left(i,k\right)}\end{array}\right]\left(19d\right)\end{array}& \left(19a\right)\end{array}$

where h is the integration step length. The first two are corresponding to equations in F in equation (4); the last two are corresponding to equations in g in equation (3). Combining these equations results in:

$\begin{array}{cc}0=H\left(x_{\gamma },y_{\gamma }\right)=\left[\begin{array}{c}{F}_1\left(x_{\gamma },y_{\gamma }\right)\\ F_2\left(x_{\gamma },y_{\gamma }\right)\\ \dots \\ g_1\left(x_{\gamma },y_{\gamma }\right)\\ g_2\left(x_{\gamma },y_{\gamma }\right)\\ \dots \end{array}\right],\mathrm{where}x_{\gamma }={\left[{\delta }_{\left(1,\gamma \right)},{\omega }_{\left(1,\gamma \right)},\dots ,{\delta }_{\left(i,\gamma \right)},{\omega }_{\left(i,\gamma \right)},\dots ,{\delta }_{\left(n,\gamma \right)},{\omega }_{\left(n,\gamma \right)}\right]}^T\mathrm{and}y_{\gamma }=\left[I_{d\left(1,\gamma \right)},I_{q\left(1,\gamma \right)},\dots ,{I_d}_{\left(i,\gamma \right)},{I_q}_{\left(i,\gamma \right)},\dots ,{I_d}_{\left(n,\gamma \right)},{I_q}_{\left(n,\gamma \right)}\right].& \left(20\right)\end{array}$

Both Newton and Gauss methods can be used to solve the equations.

1) Newton: In order to solve the equation using Newton method, one needs to find the Jacobian Matrix first. The Jacobian matrix of equation (6), restated here for convenience, is:

$\begin{array}{cc}J=\left[\begin{array}{cc}{J}_{\mathrm{Fx}}& J_{\mathrm{Fy}}\\ J_{\mathrm{gx}}& J_{\mathrm{gy}}\end{array}\right].& \left(21\right)\end{array}$

Each sub-matrix is calculated below:

$\begin{array}{cc}{J}_{\mathrm{Fx}}=\mathrm{diag}\left(J_{\mathrm{Fx}}^i\right)\mathrm{and}J_{\mathrm{Fx}}^i=\left[\begin{array}{cc}1& -h{\omega }_b/2\\ 0& 1+{\mathrm{hD}}_{\left(i\right)}/\left(4H_{\left(i\right)}\right)\end{array}\right].& \left(22\right)\\ J_{\mathrm{Fy}}=\mathrm{diag}\left(J_{\mathrm{Fy}}^i\right)\mathrm{and}J_{\mathrm{Fy}}^i=\left[\begin{array}{cc}0& 0\\ 0& {\mathrm{hE}}_{q\left(i\right)}^{\prime }/\left(4H_{\left(i\right)}\right)\end{array}\right].& \left(23\right)\\ J_{\mathrm{gx}}J_{\mathrm{gx}}=\left[\begin{array}{cccc}{J}_{\mathrm{gx}\left(11\right)}& J_{\mathrm{gx}\left(12\right)}& \dots & J_{\mathrm{gx}\left(1n\right)}\\ J_{\mathrm{gx}\left(21\right)}& J_{\mathrm{gx}\left(22\right)}& \dots & J_{\mathrm{gx}\left(2n\right)}\\ ⋮& ⋮& \ddots & ⋮\\ J_{\mathrm{gx}\left(n1\right)}& J_{\mathrm{gx}\left(n2\right)}& \dots & J_{\mathrm{gx}\left(\mathrm{nn}\right)}\end{array}\right]& \left(24\right)\\ J_{\mathrm{gx}\left(\mathrm{ii}\right)}=\left[\begin{array}{cc}\begin{array}{c}{I}_{d\left(i,\gamma \right)}\cos {\delta }_{\left(i,\gamma \right)}-I_{q\left(i,\gamma \right)}\sin {\delta }_{\left(i,\gamma \right)}+\\ G_{\left(\mathrm{ii}\right)}E_{q\left(i\right)}^{\prime }\sin {\delta }_{\left(i,\gamma \right)}+B_{\left(i,i\right)}E_{q\left(i\right)}^{\prime }\cos {\delta }_{\left(i,\gamma \right)}\end{array}& 0\\ \begin{array}{c}{I}_{d\left(i,\gamma \right)}\sin {\delta }_{\left(i,\gamma \right)}+I_{q\left(i,\gamma \right)}\cos {\delta }_{\left(i,\gamma \right)}+\\ B_{\left(\mathrm{ii}\right)}E_{q\left(i\right)}^{\prime }\sin {\delta }_{\left(i,\gamma \right)}-G_{\left(i,i\right)}E_{q\left(i\right)}^{\prime }\cos {\delta }_{\left(i,\gamma \right)}\end{array}& 0\end{array}\right].& \left(25\right)\\ J_{\mathrm{gx}\left(\mathrm{ik}\right)}=\left[\begin{array}{cc}{G}_{\left(i,k\right)}E_{q\left(k\right)}^{\prime }\sin {\delta }_{\left(k,\gamma \right)}+B_{\left(i,k\right)}E_{q\left(k\right)}^{\prime }\cos {\delta }_{\left(k,\gamma \right)}& 0\\ B_{\left(i,k\right)}E_{q\left(k\right)}^{\prime }\sin {\delta }_{\left(k,\gamma \right)}-G_{\left(i,k\right)}E_{q\left(k\right)}^{\prime }\cos {\delta }_{\left(k,\gamma \right)}& 0\end{array}\right].& \left(26\right)\\ J_{\mathrm{gy}}=\mathrm{diag}\left(J_{\mathrm{gy}}^i\right)\mathrm{and}J_{\mathrm{gy}}^i=\left[\begin{array}{cc}\sin {\delta }_{\left(i,\gamma \right)}& \cos {\delta }_{\left(i,\gamma \right)}\\ -\cos {\delta }_{\left(i,\gamma \right)}& \sin {\delta }_{\left(i,\gamma \right)}\end{array}\right].& \left(27\right)\end{array}$

Similarly, one can find the analytical expression of Jacobian matrix when higher order machine model and dynamic load model are used. Among the reviewed the references, no one suggests that there exists any difficulty to find the analytical expression of Jacobian matrix. The method to find Jacobian matrix in general case (rather than classic machine model and constant impedance) is described in equation [6]. In equations [26] and [27], the analytical expression of the Jacobian matrix is provided when transient machine model, exciter model, and dynamic load model are used.

2) Gauss: Rearranging equation (19) in the form:

X _γ =G(X_γ),  (28)

the equation can be directly solved using Gauss or Gauss-Seidel method.

Time Stacking Method—Parallel in Time

In this method, in one embodiment, the equations in (20) are combined to form a large set of algebraic equations:

$\begin{array}{cc}\{\begin{array}{c}0=H\left(x_{\gamma },y_{\gamma }\right)\\ 0=H^{\prime }\left(x_{\gamma },y_{\gamma },x_{\gamma +1},y_{\gamma +1}\right)\\ \dots \\ 0=H^{\prime }\left(x_{\gamma +p-2},y_{\gamma +p-2},x_{\gamma +p-1},y_{\gamma +p-1}\right)\end{array}& \left(29\right)\end{array}$

where p represents the number of time steps in parallel, and H and H′ are comprised of equations similar to (19). In H, the previous time step variables are known and the only unknown are the current time step variables. In H′, variables of two future adjacent time steps are unknown. The variables to be determined are:

[(x_γ,y_γ),(x_γ+1,y_γ+1), . . . ,(x_γ+p−1,y_γ+p−1)].

This equation set can be solved similarly as a conventional one-step approach. The corresponding Jacobian matrix is

$\begin{array}{cc}\Lambda =\left[\begin{array}{ccccc}{J}_{\gamma }& 0& \dots & 0& 0\\ J_{\left(\gamma +1,\gamma \right)}& J_{\gamma +1}& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & J_{\left(\gamma +p-1,\gamma +p-2\right)}& J_{\gamma +p-1}\end{array}\right]& \left(30\right)\end{array}$

where the diagonal sub-matrix is shown in equation (21). In the subscript of the off-diagonal matrix, the first number represents the equation set index, and the second number represents variable set index.

$\begin{array}{cc}{J}_{\left(j+1,i\right)}=\left[\begin{array}{cc}{J}_{F_{j+1}x_j}& J_{F_{j+1}y_j}\\ 0& 0\end{array}\right],& \left(31\right)\end{array}$

where j=γ, γ+1, . . . , γ+p−2. Each sub-matrix is calculated below:

$\begin{array}{cc}J_{F_{j+1}x_j}=\mathrm{diag}\left(J_{F_{j+1}x_j}^i\right)\mathrm{and}J_{F_{j+1}x_j}^i=\left[\begin{array}{cc}-1& -h{\omega }_b/2\\ 0& -1+{\mathrm{hD}}_{\left(i\right)}/\left(4H_{\left(i\right)}\right)\end{array}\right].& \left(32\right)\\ J_{F_{j+1}y_j}=\mathrm{diag}\left(J_{F_{j+1}y_j}^i\right)\mathrm{and}J_{F_{j+1}y_j}^i=\left[\begin{array}{cc}0& 0\\ 0& h\text{?}/\left(4H_{\left(i\right)}\right)\end{array}\right].\text{?}\text{indicates text missing or illegible when filed}& \left(33\right)\end{array}$

Numerical integration methods can be categorized by two attributes: 1) single-step or multi-step; 2) explicit or implicit.

1) Single-step methods only utilize one previous point and its derivative to determine the value at current step. On the other hand, multi-step methods express the value at current step as a function of the previous points and their derivative values. Runge-Kutta (RK) methods are the most famous family of single-step methods, while Adams family is the best known multi-step methods.

2) Explicit methods evaluate the value of current step explicitly as a function of values at previous steps, while implicit methods expresses the value at current step as a function of values at both current and previous steps.

Runge-Kutta and Adams (both explicit and implicit) methods will be briefly described below.

A. Runge-Kutta Methods

For an ordinary differential equation {dot over (x)}=f(x), RK methods have the form:

$\begin{array}{cc}{x}_{n+1}=x_n+\sum _{i=1}^sb_iK_{i,}& \left(34\right)\end{array}$

The Butcher tableaus for explicit and implicit RK methods are shown in FIGS. 3A and 3B, respectively.

1) Examples of Explicit RK Methods:

• • Forward Euler is the simplest explicit RK method with s=1 and b1=1. Its formula is x_n+1=x_n+hf(x_n).
  • One of the family of Runge-Kutta methods is so commonly used that it is often referred to as “RK4”, “classical Runge-Kutta”, or simply as “the Runge-Kutta method”. Its Butcher tableau is shown in FIG. 4.

2) Examples of Implicit RK Methods:

• • The simplest example of an implicit RK method is the backward Euler method.
  • Trapezoidal method is the second order implicit RK method.

B. Adams Methods

Adams methods have the form:

$\begin{array}{cc}\sum _{k=0}^s{\alpha }_kx_{n+k}=h\sum _{k=0}^s{\beta }_kf\left(x_{n+k}\right).& \left(35\right)\end{array}$

For Adams-Bashforth (explicit Adams methods), β=0; for Adams-Moulton (implicit Adams methods), β≠0.

Multi-step methods attempt to gain efficiency by keeping and using the information from previous steps rather than discarding it. For example, the formula of two-step Adams-Bashforth method is

$\begin{array}{cc}\begin{array}{c}{x}_{n+1}=x_n+\frac{h}{2}\left[3f\left(x_n\right)-f\left(x_{n-1}\right)\right]\\ =x_n+\mathrm{hf}\left(x_n\right)+\frac{h^2}{2}\underset{\underset{f^{\prime }\left(x_n\right)+O\left(h\right)}{}}{\frac{f\left(x_n\right)-f\left(x_{n-1}\right)}{h}}\\ =x_n+\mathrm{hf}\left(x_n\right)+\frac{h^2}{2}f^{\prime }\left(x_n\right)+O\left(h^3\right).\end{array}& \left(36\right)\end{array}$

Therefore, two-step Adams-Bashforth is equivalent to considering up to the second derivative term in the Taylor series. It should be noted that idea of multi-step integration methods is not approximating the first-order derivative at current step by combination of first-order derivative at previous steps. Instead, the method utilizes the information at previous steps to obtain a higher order approximation of Taylor series. One disadvantage of multi-step methods is the integration process has to be restarted whenever a discontinuity occurs.

Forward Euler, backward Euler, and Trapezoidal methods are degenerate cases of Adams methods, which only utilize one previous point and its derivate to determine the value at a current step.

Numerical Instability and Hyper-Stability

A numerical integration method is chosen so that:

• • if the continuous system is stable, the discrete time system is also stable.
  • if the continuous time system is unstable, the discrete time system is also unstable.

However, an integration method may create numerical instability, i.e., the continuous time system is stable, but the discrete time system is unstable. On the other hand, an integration method may also create hyper-stability, i.e., the continuous time system is unstable, but the discrete time is stable. Forward Euler, backward Euler, and Trapezoidal methods will be used to explain the idea of the above concepts.

A. Forward Euler

Let H(s) and H(z) denote the s-domain (corresponding to continuous time) and z-domain (corresponding to discrete time) transfer functions, respectively. Let s_p=σ_p+jω_p denote a pole of H(s). The continuous time system is stable if σ_p<0, for all the poles of H(s). The discrete system is stable when |z_p|<1 for all the poles of H(z). The z-domain transfer function using forward Euler integration can be obtained by

$\begin{array}{cc}H\left(z\right)=H\left(s\right){|}_{s=\frac{x-1}{h}}.& \left(37\right)\end{array}$

The poles of H(z) can be obtained by z_p=1+hs_p.

When the continuous system is unstable, i.e., there exists at least one pole with σ_p>0, then:

|z_p|=√{square root over ((1+hσ_p)²+(hω_p)²)}{square root over ((1+hσ_p)²+(hω_p)²)}≦1, ∀h≦0.  (38)

Therefore the discrete system is unstable, which means forward Euler does not create hyper-stability.

When the continuous system is stable, i.e., σ_p<0 for all the poles H(s), then:

$\begin{array}{cc}\begin{array}{c}z_p=\sqrt{{\left(1+h{\sigma }_p\right)}^2+{\left(h{\omega }_p\right)}^2}\\ =\sqrt{1+h\left({\sigma }_p^2+{\omega }_p^2\right)\left(h+\frac{2{\sigma }_p}{{\sigma }_p^2+{\omega }_p^2}\right)}\end{array}>1,\mathrm{when}h>\frac{-2{\sigma }_p}{{\sigma }_p^2+{\omega }_p^2}\left(i.e.,1+{\mathrm{hs}}_p>1\right),& \left(39\right)\end{array}$

which means when integration step length h is not small enough, forward Euler will cause numerical instability.

Defining R(hs) R as stability function of integration method, then the stability domain is

S={μεC,|R(μ)|≦1}.  (40)

The stability function of forward Euler is

R(μ)=1+μ.  (41)

Its stability domain is shown in FIG. 5A.

B. Backward Euler

The stability function of backward Euler is

$\begin{array}{cc}R\left(\mu \right)=\frac{1}{1-\mu }.& \left(42\right)\end{array}$

Its stability domain is shown in FIG. 5B.

When the continuous system is stable, i.e., σ_p<0 for all the poles of H(s), μ falls in the stability region for any integration time step length. Therefore, the discrete time system is stable, which means backward Euler method does not create numerical instability.

When the continuous system is unstable, i.e., there exists at least one pole with σ_p>0, the corresponding μ_p=hs_p falls in the stability region when h>2σ_p/(σ²_p+w²_p) (i.e., |h s_p−1|>1), which means the discrete time system is stable. Therefore, backward Euler may create hyper-stability.

C. Trapezoidal

The stability function of forward Euler is

$\begin{array}{cc}R\left(\mu \right)=\frac{1+\mu /2}{1-\mu /2}.& \left(43\right)\end{array}$

Its stability domain is shown in FIG. 5C. As can be seen, the stability property of the discrete time system is exactly the same as the continuous time system for all the time step length, i.e., the Trapezoidal method avoids both numerical instability and hyper-stability.

Classical Model with Modified Euler Method

The variables used in the code and the corresponding notation herein are listed in the table of FIG. 6. The equations are provided based on which existing code is developed.

A. Equations Translated Back from Existing Code

For each generator at time step 1 (main):

$\begin{array}{ccc}& {\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}\leftarrow \frac{S_{\mathrm{MVA}}}{{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}}& \\ & \stackrel{\_}{V}=E_{T\left(i\right)}{}^{j{\theta }_{\left(i\right)}}& \\ I=\sqrt{P_{e\left(i,1\right)}^2+Q_{e\left(i\right)}^2}/E_{T\left(i\right)}\times & \phi ={\tan }^{-1}\left(\frac{Q_{e\left(i\right)}}{P_{e\left(i,1\right)}}\right)& \stackrel{\_}{I}=I{}^{j\left({\theta }_{\left(i\right)}-\phi \right)}\\ {\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& & \\ \stackrel{\_}{E}=\stackrel{\_}{V}+\hspace{0.17em}{\mathrm{jX}}_{d\left(i\right)}^{\prime }\stackrel{\_}{I}& E_{\mathrm{im}}=\mathrm{imag}\left(\stackrel{\_}{E}\right)& E_{\mathrm{re}}=\mathrm{real}\left(\stackrel{\_}{E}\right)\\ {\delta }_{\left(i,1\right)}={\tan }^{-1}\left(\frac{\mathrm{imag}\left(\stackrel{\_}{E}\right)}{\mathrm{real}\left(\stackrel{\_}{E}\right)}\right)& {\omega }_{\left(i,1\right)}^p=1& \\ & ={\mathrm{j}}^{-{\mathrm{j\delta }}_{\left(i,1\right)}}\stackrel{\Delta }{=}{}^{j\left(90°-{\delta }_{\left(i,1\right)}\right)}& \\ {\stackrel{\_}{E}}^{\prime }=\stackrel{\_}{E}& E_{d\left(i,1\right)}^{\prime }=\mathrm{real}\left(E^{\prime }\right)& E_{q\left(i,1\right)}^{\prime }=\mathrm{imag}\left(E^{\prime }\right)\\ \stackrel{\_}{I}\leftarrow \stackrel{\_}{I}& I_{\mathrm{dg}\left(i\right)}=\mathrm{real}\left(I\right)& I_{\mathrm{qg}\left(i\right)}=\mathrm{imag}\left(I\right)\\ & I_{d\left(i\right)}=\mathrm{real}\left(I\right)/{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& I_{q\left(i\right)}=\mathrm{imag}\left(I\right)/\\ & & {\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}\\ \stackrel{\_}{V}\leftarrow \stackrel{\_}{V}& E_{d\left(i\right)}=\mathrm{real}\left(\stackrel{\_}{V}\right)& E_{q\left(i\right)}=\mathrm{imag}\left(\stackrel{\_}{V}\right)\\ & V_{\mathrm{ex}\left(i\right)}=E_{q\left(i,1\right)}^{\prime }& \\ & P_{m\left(i,1\right)}=P_{e\left(i,1\right)}{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& \end{array}$

For each time step γ, for each generator i:

P _m(i,γ+1) =P _m(i,γ)

mac_em1 (machine to system reference frame):

E _re(i,γ)=sin δ_(i,γ)E′_d(i,γ)+cos δ_(i,γ)E′_q(i,γ)real[ E′_(i,γ)/R_(i,γ)]

E _im(i,γ)=−cos δ_(i,γ)E′_d(i,γ)+sin δ_(i,γ)E′_q(i,γ)imag[ E′_(i,γ)/R_(i,γ)]

After terminal voltage Ē(_i,γ) is obtained for all the buses, for each generator i:

i_simu_innerloop2 (system reference frame):

${\stackrel{\_}{I}}_{\left(i,\gamma \right)}=\sum _k{\stackrel{\_}{Y}}_{\left(\mathrm{ik}\right)}{\stackrel{\_}{E}}_{\left(k,\gamma \right)}$ $I_{\mathrm{re}\left(i,\gamma \right)}=\mathrm{real}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}\right]$ $I_{\mathrm{im}\left(i,\gamma \right)}=\mathrm{imag}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}\right]$

mac_em2 (system to machine reference frame):

$I_{d\left(i\right)}=\sin {\delta }_{\left(i,\gamma \right)}I_{\mathrm{re}\left(i,\gamma \right)}-\cos {\delta }_{\left(i,\gamma \right)}I_{\mathrm{im}\left(i,\gamma \right)}\stackrel{\Delta }{=}\mathrm{real}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}{}_{\left(i,\gamma \right)}\right]I_{q\left(i\right)}=\cos {\delta }_{\left(i,\gamma \right)}I_{\mathrm{re}\left(i,\gamma \right)}+\sin {\delta }_{\left(i,\gamma \right)}I_{\mathrm{im}\left(i,\gamma \right)}\stackrel{\Delta }{=}\mathrm{imag}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}{}_{\left(i,\gamma \right)}\right]\begin{array}{ccc}{I}_{\mathrm{dg}\left(i\right)}=I_{d\left(i\right)}{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& I_{\mathrm{qg}\left(i\right)}=I_{q\left(i\right)}{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& \\ E_{d\left(i\right)}=E_{d\left(i,\gamma \right)}^{\prime }+& E_{q\left(i\right)}=E_{q\left(i,\gamma \right)}^{\prime }-& E_{T\left(i\right)}=\sqrt{E_{d\left(i\right)}^2+E_{q\left(i\right)}^2}\\ X_{d\left(i\right)}^{\prime }I_{\mathrm{qg}\left(i\right)}& X_{d\left(i\right)}^{\prime }I_{\mathrm{dg}\left(i\right)}& \\ P_{e\left(i,\gamma \right)}=E_{q\left(i\right)}I_{q\left(i\right)}+& Q_{e\left(i,\gamma \right)}=E_{q\left(i\right)}I_{d\left(i\right)}-& \\ E_{d\left(i\right)}I_{d\left(i\right)}& E_{d\left(i\right)}I_{q\left(i\right)}& \\ & E_{d\left(i,\gamma +1\right)}^{\prime }=E_{d\left(i,\gamma \right)}^{\prime }& E_{q\left(i,\gamma +1\right)}^{\prime }=E_{q\left(i,\gamma \right)}^{\prime }\\ & {\stackrel{.}{\delta }}_{\left(i,\gamma \right)}={\omega }_b\left[{\omega }_{\left(i,\gamma \right)}^p-1\right]& \end{array}{\stackrel{.}{\omega }}_{\left(i,\gamma \right)}^p=\frac{1}{2H_{\left(i\right)}}\left[P_{m\left(i,\gamma \right)}-P_{e\left(i,\gamma \right)}{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}-D_{\left(i\right)}\left({\omega }_{\left(i,\gamma \right)}^p-1\right)\right]S_{e\left(i\right)}=\left[E_{d\left(i\right)}+jE_{q\left(i\right)}\right]\left[I_{d\left(i\right)}-jI_{q\left(i\right)}\right])\stackrel{\Delta }{=}\left(\stackrel{\_}{V}\right){\left(\stackrel{\_}{I}\right)}^{*}$ $\begin{array}{c}\stackrel{\_}{V}=E_d+jE_q\\ =\left(E_d^{\prime }+X_d^{\prime }I_{\mathrm{qg}}\right)+j\left(E_q^{\prime }-X_d^{\prime }I_{\mathrm{dg}}\right)\\ =\left(E_d^{\prime }+jE_q^{\prime }\right)-jX_d^{\prime }\left(I_{\mathrm{qg}}+jI_{\mathrm{dg}}\right)\\ ={\stackrel{\_}{E}}^{\prime }-jX_d^{\prime }{\stackrel{\_}{I}}_g\end{array}$

B. A “Clean” but Equivalent Version of Equations

For each generator at time step 1 (main):

$\begin{array}{ccc}& {\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}\leftarrow \frac{S_{\mathrm{MVA}}}{{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}}& \\ & D_{\left(i\right)}\leftarrow \frac{D_{\left(i\right)}}{{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}}& \\ & H_{\left(i\right)}\leftarrow \frac{H_{\left(i\right)}}{{\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}}& \\ & \stackrel{\_}{V}=E_{T\left(i\right)}{}^{j{\theta }_{\left(i\right)}}& \\ I=\sqrt{P_{e\left(i,1\right)}^2+Q_{e\left(i\right)}^2}/E_{T\left(i\right)}\times & \phi ={\tan }^{-1}\left(\frac{Q_{e\left(i\right)}}{P_{e\left(i,1\right)}}\right)& \stackrel{\_}{I}=I{}^{j\left({\theta }_{\left(i\right)}-\phi \right)}\\ {\mathrm{GEN}}_{\mathrm{MVA}\left(i\right)}& & \\ & \stackrel{\_}{E}=\stackrel{\_}{V}+jX_{d\left(i\right)}^{\prime }\stackrel{\_}{I}& \\ {\delta }_{\left(i,1\right)}={\tan }^{-1}\left(\frac{\mathrm{imag}\left(\stackrel{\_}{E}\right)}{\mathrm{real}\left(\stackrel{\_}{E}\right)}\right)& {\omega }_{\left(i,1\right)}=0& \\ & E_{q\left(i\right)}=E_{\left(i,1\right)}& \\ & P_{m\left(i,1\right)}=P_{e\left(i,1\right)}& \end{array}$

For each time step γ, for each generator i:

P _m(i,γ+1) =P _m(i,γ)

mac_em1 (machine to system reference frame):

E _re(i,γ)=cos δ_(i,γ)E′_q(i)

E _im(i,γ)=sin δ_(i,γ)E′_q(i)

After terminal voltage Ē(_i,γ) is obtained for all the buses, for each generator i:

i_simu_innerloop2 (system reference frame):

${\stackrel{\_}{I}}_{\left(i,\gamma \right)}=\sum _k{\stackrel{\_}{Y}}_{\left(\mathrm{ik}\right)}{\stackrel{\_}{E}}_{\left(k,\gamma \right)}$ $I_{\mathrm{re}\left(i,\gamma \right)}=\mathrm{real}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}\right]$ $I_{\mathrm{im}\left(i,\gamma \right)}=\mathrm{imag}\left[{\stackrel{\_}{I}}_{\left(i,\gamma \right)}\right]$

mac_em2 (system to machine reference frame):

$I_{q\left(i\right)}=\cos {\delta }_{\left(i,\gamma \right)}I_{\mathrm{re}\left(i,\gamma \right)}+\sin {\delta }_{\left(i,\gamma \right)}I_{\mathrm{im}\left(i,\gamma \right)}$ $P_{e\left(i,\gamma \right)}=E_{q\left(i\right)}^{\prime }I_{q\left(i\right)}$ ${\stackrel{.}{\delta }}_{\left(i,\gamma \right)}={\omega }_b{\omega }_{\left(i,\gamma \right)}$ ${\stackrel{.}{\omega }}_{\left(i,\gamma \right)}=\frac{1}{2H_{\left(i\right)}}\left[P_{m\left(i,\gamma \right)}-P_{e\left(i,\gamma \right)}-D_{\left(i\right)}{\omega }_{\left(i,\gamma \right)}\right]$ ${\omega }_{\left(i,\gamma \right)}^p={\omega }_{\left(i,\gamma \right)}+1$

In compliance with the statute, embodiments of the invention have been described in language more or less specific as to structural and methodical features. It is to be understood, however, that the entire invention is not limited to the specific features and/or embodiments shown and/or described, since the disclosed embodiments comprise forms of putting the invention into effect.

Claims

1. A time-stacking method of simulating dynamics of a system so as to obtain time-domain trajectories which correspond to a disturbance or event, comprising: a. providing a model for the system, wherein the model includes differential equations and algebraic equations; andb. solving the differential equations and the algebraic equations over a predetermined number of time steps simultaneously using an implicit integration scheme.

2. The method of claim 1 further comprising applying a numerical solver to compute a final result of the equations.

3. The method of claim 2 further comprising applying a supercomputer or high-performance computer for solving the equations simultaneously in parallel.

4. The method of claim 1 wherein the implicit integration scheme comprises at least one of the following: Backward Euler, Trapezoidal, Runge-Kutta, Newton, and Dishonest Newton.

5. The method of claim 1 wherein the differential equations represent vectors of state variables and the algebraic equations represent vectors of algebraic variables.

6. The method of claim 5 wherein the vectors of state variables include generator dynamics.

7. The method of claim 6 wherein the generator dynamics include rotor angle and generator speed.

8. The method of claim 5 wherein the vectors of algebraic variables include voltages from selected buses in a network and current injections from generators or loads, expressed as magnitudes and phase angles, or real and imaginary parts.

9. The method of claim 1 wherein the time steps comprise intervals of at least 1 microsecond.

10. The method of claim 9 wherein the time steps of intervals vary over a range of at least 1 microsecond so as to obtain a calculation that corresponds to the disturbance or event.

11. A method of simulating power system dynamics of a system so as to obtain time-domain trajectories which correspond to a disturbance or event, comprising: a. providing a model for the system, wherein the model includes differential equations and algebraic equations;b. solving the differential equations and the algebraic equations over a predetermined number of time steps simultaneously using an implicit integration scheme;c. applying a numerical solver to compute a final result of the equations; andd. applying a supercomputer or high-performance computer for solving the equations simultaneously.

12. The method of claim 11 wherein the implicit integration scheme comprises at least one of the following: Backward Euler, Trapezoidal, Runge-Kutta, Newton, and Dishonest Newton.

13. The method of claim 11 wherein the differential equations represent vectors of state variables and the algebraic equations represent vectors of algebraic variables.

14. The method of claim 13 wherein the vectors of state variables include generator dynamics.

15. The method of claim 14 wherein the generator dynamics include rotor angle and generator speed.

16. The method of claim 13 wherein the vectors of algebraic variables include voltages from selected buses in a network and current injections from generators or loads, expressed as magnitudes and phase angles, or real and imaginary parts.

17. The method of claim 11 wherein the time steps comprise intervals of at least 1 microsecond.

18. The method of claim 17 wherein the time steps of intervals vary over a range of at least 1 microsecond so as to obtain a calculation that corresponds to the disturbance or event.