Power system transient stability simulation under uncertainty based on Taylor model arithmetic

Abstract

The Taylor model arithmetic is introduced to deal with uncertainty. The uncertainty of model parameters is described by Taylor models and each variable in functions is replaced with the Taylor model (TM). Thus, time domain simulation under uncertainty is transformed to the integration of TM-based differential equations. In this paper, the Taylor series method is employed to compute differential equations; moreover, power system time domain simulation under uncertainty based on Taylor model method is presented. This method allows a rigorous estimation of the influence of either form of uncertainty and only needs one simulation. It is computationally fast compared with the Monte Carlo method, which is another technique for uncertainty analysis. The proposed method has been tested on the 39-bus New England system. The test results illustrate the effectiveness and practical value of the approach by comparing with the results of Monte Carlo simulation and traditional time domain simulation.

Keywords

interval arithmetic / power systems / Taylor series expansion / Taylor model / time domain simulation / transient stability / uncertainty

Introduction

POWER system transient stability analysis is an important tool in power system research, which has a significant impact on the design and operation of both individual electrical utility companies and large interconnected power systems. There are two kinds of methods applied in transient stability studies, namely time domain simulation method and direct method. The time domain simulation method solves differential-algebraic equations (DAEs) to describe power system dynamics under a certain disturbance. From the time solutions of interested variables such as the rotor angles of certain generators it can find out if the system will preserve stability. The time domain simulation method has the advantage of adopting detailed models of electric components. The main methods for solving DAEs include implicit trapezoidal, modified Euler, Runge-Kutta and Taylor series methods.

In power system simulation, method selection and model selection are important; however, the obtaining of accurate model parameters also has a great influence on the accuracy of simulation results [1]. Therefore, it is very necessary to study the techniques for determining model parameters. If economic issues are ignored, it could help to obtain more accurate equipment parameters to install a large number of real-time measurement devices in power systems. Even so, the uncertainties of model parameters for power system simulation are inevitable. Such uncertainties arise mainly from two sources: one is from computational inaccuracy based on finite accuracy of computers; the other is from uncertainty of measurement for parameters of power system models. In the first case, the inaccuracies are small, at least initially. In the second case, however, inaccuracies can be large even from the beginning.

Earlier literature dealing with uncertainty in power systems are mostly about power flow and reliability assessment [2,3]. Only a few of them concern transient stability [4-6].

Recently, several methods concerning uncertainty in transient stability have been developed. Hiskens et al. have proposed a trajectory sensitivity approach to approximate the effect of uncertain parameter values on the outcome of time-step simulations [7]. The method essentially employs an augmented model that includes additional variables to represent the sensitivity of specified state variables to select parameters and initial conditions. Hockenberry et al. have explored a probabilistic collocation method (PCM) to enable the evaluation of uncertainty in power system simulations [8]. The PCM allows the uncertainty in transient behavior of power systems to be studied using only a handful of simulations, which relies on polynomial models of the relationship between the uncertain parameter in the system and the outputs of interest.

Widespread applications of interval arithmetic (IA) have appeared in recent years [9-15]. The uncertainty is an acute problem in power systems because of the complexity of power systems, so several interval arithmetic-based methods have been proposed for solving power system problems [12-15]. However, the naive use of IA is prone to pessimistic conclusions, which at times limits the industrial practicality of such methods. To decrease the overestimation of IA, the Taylor model arithmetic was proposed in 1996 by Martin Berz and his group. The Taylor model arithmetic is an extension of interval arithmetic with a comprehensive variety of applicable enclosure sets; details are discussed in Refs. [16-20].

In this paper, the Taylor model arithmetic is introduced to deal with the uncertainty problem in power system simulation. The uncertainty of model parameters is described by Taylor models (TMs) and each variable in functions is replaced with TM. Moreover, Taylor series method is adopted to compute differential equations [21]. The Taylor model method of power system time domain simulation under uncertainty is presented and realized.

Taylor model arithmetic

The Taylor model arithmetic is an extension of interval arithmetic with a comprehensive variety of applicable enclosure sets.

Definition 1 Taylor model. Let f:D\subset R^v\to R be a function that is (n+1) times continuously partially differentiable on an open set containing the domain D. Let x₀ be a point in D and P be the nth order Taylor polynomial of f around x₀. Let I be an interval such that

(1)f(x)\in P(x-x_{\mathrm{0}})+I\mathrm{for}\mathrm{all}x\in D.

Then we call the pair (P,I) an nth order Taylor model of f around x₀ on D.

Definition 2 Addition and multiplication of Taylor models. Let T_1,2=(P_1,2,I_1,2) be nth order Taylor models around x₀ over the domain D. We define

(2)T_{\mathrm{1}}+T_{\mathrm{2}}=(P_{\mathrm{1}}+P_{\mathrm{2}},I_{\mathrm{1}}+I_{\mathrm{2}}),

(3)T_{\mathrm{1}}{T}_{\mathrm{2}}=(P_{\mathrm{1.2}},I_{\mathrm{1.2}}),

where P_1.2 is the part of the polynomial P₁P₂ up to order n and I_1.2=B(P_e)+B(P₁)I₂+B(P₂)I₁+I₁I₂, where P_e is the part of the polynomial P₁P₂ of orders (n+1) to 2n, and B(P) denotes a bound of P on the domain D.

Definition 3 Intrinsic functions of Taylor models. Let T = (P, I) be a Taylor model of order n over the v-dimensional domain D = [a,b] around the point x₀. In the following, let f(x)∈P(x-x₀) +I be any function in the Taylor model, let c_f = f(x₀), and \overline{f} be defined by \overline{f}(x)=f(x)-c_f. Likewise, we define \overline{P} by \overline{P}(x-x_{\mathrm{0}})=P(x-x_{\mathrm{0}})-c_f, so that \left(\overline{P},I\right) is a Taylor model for \overline{f}. Some useful intrinsic functions which are needed in this paper are listed as follows.

1) Multiplicative inverse. Under the condition \forall x\in D, \mathrm{0}\notin B[P(x-x_{\mathrm{0}})+I], we write as follows:

(4){\displaystyle \frac{\mathrm{1}}{f(x)}}={\displaystyle \frac{\mathrm{1}}{c_f}}\left\{\mathrm{1}-{\displaystyle \frac{\overline{f}(x)}{c_f}}+{\displaystyle \frac{{\left(\overline{f}(x)\right)}^{\mathrm{2}}}{c_f^{\mathrm{2}}}}-\cdots +{(-\mathrm{1})}^k{\displaystyle \frac{{\left(\overline{f}(x)\right)}^k}{c_f^k}}\right\}+{(-\mathrm{1})}^{k+\mathrm{1}}{\displaystyle \frac{{\left(\overline{f}(x)\right)}^{k+\mathrm{1}}}{c_f^{k+\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+\theta \overline{f}(x)/c_f\right)}^{k+\mathrm{2}}}},

where 0<θ<1.

2) Sine.

(5)\sin f(x)=\sin c_f+\cos c_f\overline{f}(x)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}!}}\sin c_f{\left(\overline{f}(x)\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}!}}\cos c_f{\left(\overline{f}(x)\right)}^{\mathrm{3}}+\cdots +{\displaystyle \frac{\mathrm{1}}{(k+\mathrm{1})!}}{\left(\overline{f}(x)\right)}^{k+\mathrm{1}}J,

where

J=\{\begin{array}{l}-J_{\mathrm{0}},\mathrm{if}\mathrm{mod}(k,\mathrm{4})=\mathrm{1},\mathrm{2},\hfill \\ J_{\mathrm{0}},\mathrm{else},\hfill \end{array}

J_{\mathrm{0}}=\{\begin{array}{l}\cos \left(c_f+\theta \overline{f}(x)\right),\mathrm{if}k\mathrm{is}\mathrm{even,}\hfill \\ \sin \left(c_f+\theta \overline{f}(x)\right),\mathrm{else}.\hfill \end{array}

3) Cosine.

(6)\cos f(x)=\cos c_f-\sin c_f\overline{f}(x)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}!}}\cos c_f{\left(\overline{f}(x)\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}!}}\sin c_f{\left(\overline{f}(x)\right)}^{\mathrm{3}}+\cdots +{\displaystyle \frac{\mathrm{1}}{(k+\mathrm{1})!}}{\left(\overline{f}(x)\right)}^{k+\mathrm{1}}J,

where

J=\{\begin{array}{l}-J_{\mathrm{0}},\mathrm{if}\mathrm{mod}(k,\mathrm{4})=\mathrm{0},\mathrm{1},\hfill \\ J_{\mathrm{0}},\mathrm{else},\hfill \end{array}

J_{\mathrm{0}}=\{\begin{array}{l}\sin \left(c_f+\theta \overline{f}(x)\right),\mathrm{if}k\mathrm{is}\mathrm{even,}\hfill \\ \cos \left(c_f+\theta \overline{f}(x)\right),\mathrm{else}.\hfill \end{array}

In a word, we can obtain Taylor models of arbitrary functions by overloading the operations and basic math functions in our program.

Taylor model based transient stability simulation

For the sake of clarity, in this paper, TM variables are bracketed with ‘[ ]’ in order to distinguish them from fixed numbers.

Taylor model based differential equations

If some equation variables are Taylor models, the differential equations will be transformed to Taylor model based differential equations. The solving of differential equations needs a certain integral method. The higher order Taylor series method allows us to adopt a long integration step and thus shortens the integration time and decreases the number of arithmetic operations. For Taylor model arithmetic, the less number of arithmetic operations, the less overestimation the final result will preserve. Therefore, in this paper, Taylor series method is adopted to solve Taylor model based differential equations.

Taylor model based transient stability simulation for classical models

The center of inertia (COI) reference frame is adopted [13]. For simplicity and to highlight the characteristics of the proposed algorithm, the classical model of generators is adopted.

The center of inertia is defined by

(7)\{\begin{array}{l}[{\delta }_O]={\displaystyle \frac{\mathrm{1}}{[M_T]}}{\displaystyle \sum _{i=\mathrm{1}}^n[M_i][{\delta }_i],}\hfill \\ \left[{\stackrel{·}{\delta }}_O\right]={\displaystyle \frac{\mathrm{1}}{[M_T]}}{\displaystyle \sum _{i=\mathrm{1}}^n[M_i]\left[{\stackrel{·}{\delta }}_i\right],}\hfill \end{array}

where [δ_i] is the angle of generator i; [M_i] is the inertia constant of generator i; and [M_T]={\displaystyle \sum _{i=\mathrm{1}}^n[M_i]}.

The swing equations in COI reference frame are

(8)\{\begin{array}{l}[M_i]{\displaystyle \frac{\mathrm{d}\left[{\tilde{\omega }}_i\right]}{\mathrm{d}t}}=[P_{mi}]-[P_{ei}]-{\displaystyle \frac{[M_i]}{[M_T]}}[P_{\mathrm{COI}}],\hfill \\ {\displaystyle \frac{\mathrm{d}[{\theta }_i]}{\mathrm{d}t}}=\left[{\tilde{\omega }}_i\right],\hfill \end{array}

where [θ_i]=[δ_i]-[δ_O]; \left[{\tilde{\omega }}_i\right]\left[{\stackrel{·}{\delta }}_i\right]-\left[{\stackrel{·}{\delta }}_O\right]; [P_{\mathrm{COI}}]={\displaystyle \sum _{i=\mathrm{1}}^n([P_{mi}]-[P_{ei}])}; [P_mi] is the mechanical input of generator i; [P_ei] is the electrical output of generator i; [P_COI] is the COI accelerating power; [θ_i] and \left[{\tilde{\omega }}_i\right] are the rotor angle and speed of generator i in COI notation, respectively.

[P_ei] can be derived by

(9)\left[P_{ei}\right]=E_i^{\mathrm{2}}{G}_{ii}+{\displaystyle \sum _{j=\mathrm{1},j\ne i}^n\left(C_{ij}\sin [{\theta }_{ij}]+D_{ij}\cos [{\theta }_{ij}]\right),}

where C_ij = E_iE_jB_ij, D_ij = E_iE_jG_ij. B_ij and G_ij are susceptance and conductance between the internal generator nodes i and j, respectively. E_i and E_j are the electromotive force (EMF) of internal generators i and j, respectively.

The time domain simulation with Taylor series method under classical models needs arbitrary order derivatives of electrical power [P_ei], which can be obtained by

(10)[P_{ei}^{(m)}]={\displaystyle \sum _{j=\mathrm{1},j\ne i}^n(C_{ij}{\sin }^{(m)}[{\theta }_{ij}]+D_{ij}{\cos }^{(m)}[{\theta }_{ij}]).}

The time domain simulation process is actually the solution of Eqs. (8)-(10).

Implementation of arithmetic

A power system transient stability simulation based on Taylor model arithmetic is realized by C++ standard language. First, we design an interval class and Taylor model class, and define their various basic operations. Then the Taylor series method is used to implement transient stability simulation. Several programming skills are adopted to reduce calculation time significantly as follows:

1) Although we can easily design two functions for example DSIN() and DCOS() in terms of Eqs. (9) and (10) and then call them to calculate sin^(m)[θ_ij] and cos^(m)[θ_ij], this may involve many redundant computations. Taking the computation of sin^(m)[θ_ij] as an example, the 0th to (m–1)th derivatives of cos[θ_ij] will be used in the computation of sin^(m)[θ_ij] and 0th to mth derivatives of cos[θ_ij] will be used in the computation of sin^(m+1)[θ_ij]. If each derivative is recalculated and very higher Taylor series method is adopted, the computation time will be substantially increased. A better method is to save sin^(m)[θ_ij] in an array to avoid recalculating. For example, DSIN[i][j][m] is used to save sin^(m)[θ_ij] in the program.

2) Because:

(11){\sin }^{(m)}[{\theta }_{ij}]=-{\sin }^{(m)}[{\theta }_{ji}],

(12){\cos }^{(m)}[{\theta }_{ij}]={\cos }^{(m)}[{\theta }_{ji}].

We only need calculate the upper triangular matrixes of DSIN[ ][ ][ ] and DCOS[ ][ ][ ].

3) The problem of [θ_ij]^(m) is similar to sin^(m)[θ_ij], namely, [θ_ij]^(m) will be used many times in the computation, which can be saved in an array DDELT[i][j][m]. This method can reduce many subtraction operations, especially, at sufficiently high Taylor order.

By adopting the above three programming skills, it not only reduces computation time, thus significantly increasing computing speed, but also reduces the overestimation of the final result through reduction of computation..

Test result and analysis

The proposed method has been tested on the 10-unit 39-bus New England system with uncertain parameters in generator models. The system data are fully taken from Ref. [22]. It is assumed that a three-phase fault on bus 26 cleared at 5th cycle, the tripped line is 26-27. The one line diagram of the 10-unit 39-bus New England system is given in Fig. 1.

For the convenience of comparison, inertia constants and mechanical inputs of all generators are listed in Table 1.

The simulation results of the proposed method with a time step of 0.02 s are illustrated and compared to those of the Monte Carlo method and traditional time domain simulation method in some tables and figures.

Note that each figure listed below has five curves with some of which being overlapped probably. The top curve and bottom curve with the symbol ‘•’ show the bounds of the Taylor model result when a certain parameter is an interval. The middle curve with the symbol ‘+’ shows the result of traditional time domain simulation when the parameters are set the midpoints of given intervals. The other two curves with the symbol ‘*’ give the bounds of the Monte Carlo result when a certain parameter varies in an interval.

Since the rotor angle of the ith generator is most influenced by the uncertainty of its own parameters, we take the rotor angle of the ith generator as the example to plot result figures when one of its parameters is an interval.

1) The detailed simulation results with the above three methods at each integration step are listed in Table 2 when the mechanical input of G1 is an interval [0.50 P_m1, 1.50 P_m1], in accordance to which the resulting figures are shown in Fig. 2.

2) The resulting figures are shown in Fig. 3 when the mechanical input of G10 is an interval [0.50 P_m10, 1.50 P_m10].

3) The resulting figures are shown in Fig. 4 when the inertia constant of G1 is an interval [0.70 H₁,1.30 H₁].

4) The resulting figures are shown in Fig. 5 when the inertia constant of G1 is an interval [0.80 H₁,1.20 H₁].

5) The resulting figures are shown in Fig. 6 when the inertia constant of G8 is an interval [0.90 H₈,1.10 H₈].

6) The resulting figures are shown in Fig. 7 when the inertia constant of G8 is an interval [0.80 H₈,1.20 H₈].

From the above six resulting figures, we can obtain the following implications:

1) The result from Taylor model arithmetic when a certain parameter of the generator is an interval encloses the result of the traditional time domain simulation method when the parameter is set the midpoint of a given interval. At the same time, it also encloses the result calculated by Monte Carlo method. This demonstrates the validity of the proposed method.

2) It can be observed from Figs. 2 and 3 that the results of Taylor model arithmetic almost entirely overlap those of the Monte Carlo method. Plenty of tests show that in most cases the proposed method can obtain the simulation result as nearly precise as the Monte Carlo method can, and even in some cases the proposed method can obtain as precise result as the Monte Carlo method. However, judging the computation speed, the proposed method has absolute advantages because Monte Carlo simulation spends a lot of time in repeating simulations (e.g., 50000 times) while selecting random points in the appointed interval.

3) In Taylor model arithmetic, if inertia constants and mechanical inputs vary with the same proportion, the former may have a bigger influence than the latter on the overestimation of the final result. This is due to the fact that generator inertia constants appear more often in equations than mechanical inputs.

4) By comparing Figs. 4 and 7, we can observe that in Taylor model arithmetic, when varying with the same proportion, the resulting interval from the inertia constant of G1 is sharper than that of G8 with the inertia constant of G1 being larger than that of G8. This is due to the fact that generator inertia constants are often to be used as divisors in calculation.

Conclusions

Taylor model arithmetic is an effective method to deal with uncertainty problems by decreasing the overestimation of interval arithmetic.

Power system transient stability simulation under uncertainty based on Taylor model arithmetic ensures the simulation result contains all possible solutions when the parameters vary in the given domain.

The proposed method has been validated by comparing with the Monte Carlo simulation method. The proposed method is computationally superior to Monte Carlo simulation method although the result of the former is not too much conservative than that of the latter.

References

[1] Kosterev D N, Tylor C W, Mittelstadt W A. Model validation for the August 10, 1996 WSCC system outage. IEEE Transactions on Power Systems, 1999, 14(3): 967–979

[2] Dopazo J F, Klitin O A, Sasson A M. Stochastic load flows. IEEE Transactions on Power Apparatus and Systems, 1975, 94(2): 299–309

[3] Allan R N, Billinton R, Breipohl A M, Grigg C H. Bibliography on the application of probability methods in power system reliability evaluation 1992-1996. IEEE Transactions on Power Systems, 1999, 14(1): 51–57

[4] Timko K J, Bose A, Anderson P M. Monte Carlo simulation of power system stability. IEEE Transactions on Power Apparatus and Systems, 1983, PAS-102(10): 3453–3459

[5] Anderson P M, Bose A. A probabilistic approach to power system stability analysis. IEEE Transactions on Power Apparatus and Systems, 1983, PAS-102(8): 2430–2439

[6] Wu F F, Tsai Y K. Probabilistic dynamic security assessment of power systems I: basic model. IEEE Transactions on Circuits and Systems, 1983, 30(3): 148–159

[7] Hiskens I A, Pai M A, Nguyen T B. Bounding uncertainty in power system dynamic simulations. In: Proceedings of IEEE Power Engineering Society Winter Meeting, Singapore. 2000, 2: 1533–1537

[8] Hockenberry J R, Lesieutre B C. Evaluation of uncertainty in dynamic simulations of power system models: the probabilistic collocation method. IEEE Transactions on Power Systems, 2004, 19(3): 1483–1491

[9] Moore R E. Interval arithmetic and automatic error analysis in digital computing. Dissertation for the Doctoral Degree. Stanford: Stanford University, 1962

[10] Alefeld G, Mayer G. Interval analysis: theory and applications. Journal of Computational and Applied Mathematics, 2000, 121(1-2): 421–464

[11] Hickey T, Ju Q, Van Emden M H. Interval arithmetic: from principles to implementation. Journal of the ACM, 2001, 48(5): 1038–1068

[12] Wang Z, Alvarado F L. Interval arithmetic in power flow analysis. IEEE Transactions on Power Systems, 1992, 7(3): 1341–1349

[13] Das B. Radial distribution system power flow using interval arithmetic. International Journal of Electrical Power & Energy Systems, 2002, 24(10): 827–836

[14] Sarić A T, Stanković A M. An application of interval analysis and optimization to electric energy markets. IEEE Transactions on Power Systems, 2006, 21(2): 515–523

[15] Chaturvedi A, Prasad K, Ranjan R. Use of interval arithmetic to incorporate the uncertainty of load demand for radial distribution system analysis. IEEE Transactions on Power Delivery, 2006, 21(2): 1019–1021

[16] Berz M, Hoffstatter G. Computation and application of Taylor polynomials with interval remainder bounds. Reliable Computing, 1998, 4(1): 83–97

[17] Makino K. Rigorous Analysis of nonlinear motion in particle accelerators. Dissertation for the Doctoral Degree. East Lansing: Michigan State University, 1998

[18] Neumaier A. Taylor forms–use and limits. Reliable Computing, 2003, 9(1): 43–79

[19] Revol N, Makino K, Berz M. Taylor models and floating-point arithmetic: proof that arithmetic operations are validated in COSY. Journal of Logic and Algebraic Programming, 2005, 64(1): 135–154

[20] Makino K, Berz M. Taylor models and other validated functional inclusion methods. International Journal of Pure and Applied Mathematics, 2003, 4(4): 379–456

[21] Barrio R. Performance of the Taylor series method for ODEs/DAEs. Applied Mathematics and Computation, 2005, 163(2): 525–545

[22] Pai M A. Energy Function Analysis for Power System Stability. Norwell, MA: Kluwer Academic Publishers, 1989

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg