Accepted by 2017 IEEE Power & Energy Society General Meeting, Chicago, IL, July 16-20, 2017

Semi-Analytical Fault-on Trajectory Simulation and Its

Application in Direct Methods
Eric Abreut1, Bin Wang2, Student Member, IEEE, Kai Sun2, Senior Member, IEEE
1 2
Florida International University, Miami, FL University of Tennessee, Knoxville, TN
[email protected] {bwang, kaisun}@utk.edu

Abstract— This paper proposes a time series based semi- This method focuses on creating a numerical solution of the
analytical solution (SAS) which can approximately represent the nonlinear differential algebraic equations and takes a
solution of dynamical power system equations within a short moderate amount of time to complete. Although it is often
period of time. All system parameters and initial conditions are
represented by symbolic variables such that different system
accurate, it does not provide much information about the
operating conditions and contingencies can be completely stability margin of the system.
considered in the evaluation of any power system. The proposed The second method is known as the energy function
SAS is accurate at the initial point while its error increases over method or in other words the direct method [2]-[5]. This
time. It is observed that using the SAS for approximating the method focuses on evaluating a predefined energy function
fault-on trajectory is fairly accurate. Since faults in a power for a power system and works almost in real time. This
system are cleared after a short period of time, i.e 5 or 6 cycles
(.08 to .1 seconds), then the proposed SAS can be integrated with
allows for a screening of the stability of a power system to be
any direct methods to access the transient stability of the system. quick and somewhat accurate. Although it is not as accurate
A decoupling based direct method is used in this paper along as the time-domain method, it does determine a stability
with the proposed SAS for demonstration purposes. Case studies margin with a low error index and can be quite reliable in
performed on the IEEE 9-bus system and the simplified WECC many cases. Conducting continuous TSA’s requires fast
179-bus system show the accuracy of the proposed SAS and its processing and real-time calculations to determine the
application in transient stability analysis using direct methods.
stability of the system during a disturbance. The benefit of
Index Terms—Semi-analytical solution (SAS), time series, using the direct method is that its real-time speed allows for
simulation, direct method, linear decoupling, transient stability, online applications of TSA that are as fast as real-time
stability margin, transient energy function. monitoring. This will help monitor a power system using an
online program without any delays in the determination of the
system’s stability.
I. INTRODUCTION The basic idea behind a traditional direct method is to
Power systems must work under stable conditions in order compare the initial energy and the critical energy of the
to be able to function securely. Whenever a fault occurs, the system based on a pre-defined energy function. The critical
system’s stability begins to fluctuate. If the stability reaches a energy of the system defined at a specific point on its
point of no return, then the system completely loses its boundary, also known as controlling unstable equilibrium
stability and blackout may take place [1]. In order to take point (CUEP), is a point that the system cannot exceed or else
action and prevent the system from losing its stability, a it will completely lose its stability and never return to the
power system transient stability assessment (TSA) is used to post-fault stable equilibrium point. If the critical energy of the
guarantee both security and stability when operating the system is greater than its initial energy, then the system is
power system. stable and the stability margin can be determined by their
When monitoring a power system, a constantly running difference. Otherwise, the system is unstable. Still, there is
TSA must be conducted to determine whether or not the one step in a direct method that requires time-domain
system is going to remain stable or cause problems under simulation, i.e. the simulation of the fault-on trajectory.
credible contingencies in the near future. Therefore, TSA is Instead, if the semi-analytical solution (SAS) is acquired,
an extremely important routine in safety and security of the then it can be used for simulating the behaviors of a
power system and can be conducted using two different dynamical system [6]-[9]. This enables a possibility for
methods that measure the stability of a power system. replacing the time-domain simulation for fault-on period
The first method is known as the time-domain method. found in the direct method and also increasing the speed.

This work was supported in part by NSF CURENT Engineering Research

Center and NSF grant ECCS-1610025.
 This paper proposes an SAS, i.e. an explicit function in ⎧b10 − a11 = 0 ⎧b20 − a21 = 0 ⎧bn 0 − an1 = 0
time, for the fault-on system which, when evaluated at the ⎪b − 2a = 0 ⎪b − 2a = 0 ⎪
fault clearing time, can give the system states required for
⎪ 11 12 , ⎪ 21 22 ,…, ⎪bn1 − 2an 2 = 0 (5)
⎨ ⎨ ⎨
calculating the initial energy. Thus, time-domain simulation ⎪ M ⎪M ⎪M
for the fault-on period is avoided and simulation-free direct ⎪⎩b1m = 0 ⎪⎩b2 m = 0 ⎪⎩bnm = 0
methods can be achieved.
There are different direct methods differing from each These (m+1)×n equations are usually enough for solving
other in conservativeness and reliability. For convenience, the the (m+1)×n unknown coefficients aij. After obtaining aij,
decoupling based direct method [5] is adopted in this substituting them back into (2) will give the explicit semi-
research. This method uses linear decoupling transformation analytical solution for the ODEs in (1).
and constructs a fictitious single machine infinite bus (SMIB) Note that (i) the coefficients in the SAS only depend on
power system for each oscillation mode. The transient energy the initial conditions x0 and the function f; (ii) semi-analytical
function method is applied to each of the SMIB systems and solution is accurate at the initial point and deviates from the
then the stability margins are calculated. true solution over time; (iii) the error is small, or ignorable,
only within a short period, e.g. 0.1s-0.2s for power systems
II. TIME SERIES BASED SEMI-ANALYTICAL SOLUTION which will be investigated in the case study section.
This section introduces the derivation of the time series
III. APPLICATION FOR SIMULATING POWER SYSTEM
based SAS for dynamical systems governed by ordinary FAULT-ON TRAJECTORY
differential equations (ODEs).
Given a dynamical system in (1). This subsection introduces the application of the time
series based SAS to power system fault-on trajectory
simulation. Consider an N-machine power system ODEs in
x& = f ( x) (1)
(6), where all machines are represented by classical model and
where x is an n×1 vector containing the state variables and f is all loads are constant impedances.
a linear/nonlinear vector field.
Assume (2) to be an approximate of the solution to (1) Di & ωs
δ&&i + δi + ( Pmi − Pei ) = 0 (6)
with initial conditions x0, where m is a pre-specified order and 2H i 2Hi
aij are unknowns to be determined.
m
Pei = Ei2Gi + ∑ (C ij sin(δ i − δ j ) + Dij cos(δ i − δ j ) ) (7)
⎛ x1 (t ) ⎞ ⎛ a10 + a11t + L + a1mt m ⎞ j =1, j ≠ i

⎜ ⎟ ⎜ ⎟
x (t ) a + a t + L + a2 mt m ⎟ where i ∈{1,2,…,m}, δi, Pmi, Pei, Ei, Hi and Di respectively
x(t ) = ⎜ 2 ⎟ ≈ ⎜ 20 21 (2)
⎜ M ⎟ ⎜ M ⎟ represent the absolute rotor angle, mechanical power,
⎜ ⎟ ⎜⎜ m⎟
⎟ electrical power, field voltage, the inertia constant and
⎝ xn (t ) ⎠ ⎝ an 0 + an1t + L + anmt ⎠ damping constant of machine i, and Gi, Cij, and Dij represent
network parameters including loads; ωs represents the
Then, the derivative of x is calculated and shown in (3).
synchronous angular frequency, i.e. 50Hz or 60Hz.
Substitute (2) into the vector field f and calculate the Taylor
Following the procedures introduced in section II, the final
series in t up to the order m and obtain (4), where each b is a
equations set corresponding to (5) has such a characteristic
function of a.
(note that n=2N): the equation set consisting of the k-th row of
each equation set in (5), as shown in (8) where a=0 when k=m,
⎛ x&1 (t ) ⎞ ⎛ a11 + 2a12t + L + ma1mt m −1 ⎞ can always be solved as a linear problem with the solution of
⎜& ⎟ ⎜ ⎟
⎜ x2 (t ) ⎟ ⎜ a21 + 2a22t + L + ma2 mt m −1 ⎟ (3)
(8) for k-1, k-2, …, 2 and 1. Such characteristic guarantees the
x& (t ) = = solvability of the problem in (5) for power systems. In
⎜ M ⎟ ⎜ M ⎟
⎜ ⎟ ⎜⎜ ⎟
m −1 ⎟
addition, this characteristic also allows the symbolic SAS.
⎝ x&n (t ) ⎠ ⎝ an1 + 2an 2t + L + manmt ⎠ ⎧b1, k +1 − ka1k = 0
⎪b
⎛ a10 + a11t + L + a1mt m ⎞ ⎛ b10 + b11t + L + b1mt m ⎞ ⎪ 2, k +1 − ka2 k = 0 (8)
⎜ ⎟ ⎜ ⎟ ⎨
⎜ a20 + a21t + L + a2 m t m ⎟ ⎜ b20 + b21t + L + b2 mt m ⎟ (4) ⎪M
f ( x) = f ≈ ⎪⎩bn ,k +1 − kank = 0
⎜ M ⎟ ⎜ ⎟
⎜⎜ m⎟
⎟ ⎜⎜ m⎟
⎟
⎝ an 0 + an1t + L + anmt ⎠ ⎝ bn 0 + bn1t + L + bnm t ⎠ The following will use the derived symbolic SAS for
simulating the fault-on dynamics of power systems. Note that
Substitute (3) and (4) into (1), equate the terms on both
for power systems under different fault conditions, the only
sides by term in t and obtain (m+1)×n equations as shown in
difference is found in the parameters Gi, Cij and Dij in (7). This
(5).
does not impact the derived symbolic solution at all. Thus, for
different faults the only thing we need to do is to use their
(
⎧a12 = D1δ&10 / ωs + G1 E12 + C12 sin δ120 + C13 sin δ130
⎪
corresponding parameters for Gi, Cij and Dij, while everything ⎪ + D12 cos δ120 + D13 cos δ130 − Pm1 ) ⋅ ωs 4 H1
else remains the same. ⎪
Since a short-circuit event, or usually called fault, in a
⎪
(
⎪a13 = 2a12 D1δ&10 / ωs + C12δ120 cos δ120 + C13δ130 cos δ130
power system is always cleared by the protection system no
later than 5-6 cycles, i.e. 0.08s-0.1s, after its occurrence, the
⎪ − D δ 0 cos δ 0 − D δ 0 cos δ 0 ) ⋅ ω 12 H
⎪ 12 12 12 13 13 13 s 1
SAS should be valid within that period and can be directly
used for determining any point on the fault-on trajectory. With ⎪ 22
(
⎪a = D δ& / ω + G E + C sin δ + C sin δ 0
2 2
0
s 2 2
2
21
0
21 23 23

the last system state and the decoupling based direct method, ⎪ + D cos δ 0 + D cos δ 0 − P ) ⋅ ω 4 H
the transient stability of the system under this fault can be ⎪ 21 21 23 23 m2 s 2
⎨
evaluated. Such application will be demonstrated in the next
section.
( &
⎪a23 = 2a22 D2δ 2 / ωs + C21δ 21 cos δ 21 + C23δ 23 cos δ 23
0 0 0 0 0

⎪
⎪ − D21δ 21 cos δ 21 − D23δ 23 cos δ 23 ) ⋅ ωs 12 H 2
0 0 0 0

IV. CASE STUDIES ⎪

A. Tests on IEEE 9-bus power system
(
⎪a32 = D3δ&3 / ωs + G3 E3 + C31 sin δ 31 + C32 sin δ 32
0 2 0 0

⎪
⎪ + D31 cos δ 310 + D32 cos δ 320 − Pm3 ) ⋅ ωs 4 H 3
⎪
(
⎪a33 = 2a32 D3δ&30 / ωs + C31δ 310 cos δ 310 + C32δ 320 cos δ 320
⎪
⎪ − D31δ 310 cos δ 310 − D32δ 320 cos δ 320 ) ⋅ ωs 12 H 3 (10)
⎩
where δij0 = δi0 - δj0.

Fig. 1. IEEE 9-bus power system

By the proposed approach, the 3rd-order time series based

SAS of the IEEE 9-bus system [10] is shown in (9), where ai0
and ai1 are initial conditions δi0 and δ&i 0 of generator i,
respectively, and ai2 and ai3 are derived as shown in (10).

In the first test, a three-phase fault is added at bus 7 and

cleared after 5 cycles by tripping the line 5-7. The fault-on Fig. 2. Comparison of fault-on trajectories by SAS and RK4
trajectory by the proposed 3rd-order SAS and the 4th-order
With the last point from the fault-on trajectory, the
Runge-Kutta (RK4) method are compared in Fig. 2, where the
decoupling based direct method can be used to determine the
visible red curve is by the proposed SAS the blue curve
stability margin of the system. The second test investigates
(which cannot be seen under the red curve) is by RK4. Testing
different faults, where the fault duration for each case is fixed
all other faults will give similar phenomena. Thus, it can be
at 5 cycles. The result is shown in TABLE I, where the
clearly seen that the proposed SAS is fairly accurate during
stability margins for the two oscillatory modes are shown in
the fault-on period.
columns three and four and the smaller one is selected as the
stability margin of the case. Column five gives the ranking of
⎛ δ1 (t ) ⎞ ⎛ a10 + a11t + a12t + a13t ⎞
2 3
the severities of these faults by the decoupling based direct
⎜ ⎟
δ (t ) = ⎜⎜ δ 2 (t ) ⎟⎟ = ⎜ a20 + a21t + a22t 2 + a23t 3 ⎟ (9) method while the last two columns give the ranking by the
⎜ δ (t ) ⎟ ⎜ a + a t + a t 2 + a t 3 ⎟ CCT. It can be seen that the two rankings roughly align with
⎝ 2 ⎠ ⎝ 30 31 32 33 ⎠ each other.
 TABLE I STABILITY MARGIN AND RANKING TABLE less severe such that errors from SAS may not have significant
Tripped Fault Mode Mode 2 Ranking ΔVn Ranking CCT
Line Bus 1 by ΔVn by CCT /s
impacts on the analysis results.
5-7 7 2.09 4.0e+03 1 2.09 1 .179 In addition, it is found that a higher-order SAS will give
7-8 7 5.81 169.91 2 5.81 2 .195 better accuracy within the same time window than that at a
5-7 5 8.72 4.5e+03 3 8.72 8 .353 lower order. However, a higher-order SAS has a longer
6-9 9 8.99 145.52 4 8.99 3 .231 expression and will cost more time to evaluate.
7-8 8 13.02 80.50 5 13.02 5 .297
8-9 8 20.55 67.29 6 20.55 6 .324
8-9 9 34.97 21.11 7 21.11 4 .249
6-9 6 25.90 1.2e+03 8 25.90 9 .430
4-6 4 33.19 7.4e+03 9 33.19 7 .329
4-6 6 50.98 998.22 10 50.98 10 .493

In the third test, the fault durations of different faults are

increased to corresponding critical clearing times (CCTs) and
lead to marginally stable cases. The purpose of these tests are
to demonstrate the accuracy of the SAS over a longer time.
To quantity the error of the proposed SAS, the largest Fig. 3. Example one (0.231 second fault-on simulation)
angle difference between the trajectory by SAS and the
trajectory by RK4 is selected as an index.
TABLE II shows the errors of the SAS with fault durations
at 5/60 seconds and CCT, respectively, for each fault. It can be
seen that (i) the error is fairly small, e.g. 2e-3 deg, if only
taking a 5-cycle fault duration into consideration; (ii) the error
is small for faults with small CCTs while large for faults with
large CCTs. This is reasonable since the SAS will be less
accurate over a longer time window.

TABLE II STABILITY MARGIN ERROR INDEX

CCT Error Index in Degrees Error Index in Degrees
(fault cleared after 5 cycles) (Fault cleared at CCT)
.179 0.001203 0.20
.195 0.001203 0.34 Fig. 4. Example two (0.430 second fault-on simulation)
.231 0.000486 0.49
.249 0.000486 0.77 B. Tests on WECC 179-bus power system
.297 0.001325 4.35 The second test system was conducted on a simplified
.324 0.001325 7.21 WECC power system [11], which has 29 machines and 179
.329 0.000270 1.38 buses as shown in Fig. 7. This test aims at demonstrating the
.353 0.000721 6.69 accuracy and the applicability of the proposed SAS on large
.430 0.000474 19.94 power systems and its application in direct methods.
.493 0.000474 41.57 Following the same procedures in Section IV-A, once the
fault-on trajectory is obtained from the SAS and the stability
The following gives two examples from TABLE II. In margins are calculated for each case using the final states of
example one, a three-phase fault is added at bus 9 and cleared the system. The stability margins can be ranked into a table
after 0.231s by tripping line 6-9 and the fault-on trajectory is from lowest to the greatest. Different cases used in this study
shown in Fig.3. In the second example, a three-phase fault is were picked every few cases from a low CCT to a higher CCT
added at bus 6 and cleared after 0.430s by tripping line 6-9 from all possible N-1 contingencies. Since the fault duration
and the resulting trajectory is shown in Fig.4. In both Fig. 3 used in each of these cases was .05 seconds, then any CCT
and 4, the red curves represent the trajectories by the proposed lower than .05 would make the system unstable. Any of the
SAS and blue curves are by the RK4. cases in whose CCTs are greater than .05 should be stable.
Although the error increases over a longer time window, The results are shown in TABLE III.
that should not discourage the use of the proposed SAS since The results shown in TABLE III show that all unstable
the faults with small CCTs can be accurately handled and they cases have been captured by the proposed methods, which
are the ones that may cause stability problems which should demonstrates the applicability to large power systems. Fig. 6
be of more interests. For those faults with large CCTs, they are
shows the trajectories from 4th-order SAS and RK4, VI. ACKNOWLEDGEMENT
respectively, which are close to each other. The authors want to thank Mr. Nan Duan from the
department of EECS at the University of Tennessee,
Knoxville, for his insightful comments on this work.

REFERENCES
[1] P. Kundur et al., "Definition and classification of power system
stability IEEE/CIGRE joint task force on stability terms and
definitions," IEEE Trans. on Power Syst., vol. 19, no. 3, pp. 1387-1401,
Aug. 2004.
[2] Chiang, Hsiao-Dong. Direct Methods for Stability Analysis of Electric
Power Systems: Theoretical Foundation, BCU Methodologies, and
Applications. Hoboken, US: Wiley, 2010.
[3] P. Varaiya, F. F. Wu and Rong-Liang Chen, "Direct methods for
transient stability analysis of power systems: Recent results," in
Proceedings of the IEEE, vol. 73, no. 12, pp. 1703-1715, Dec. 1985.
[4] M. Pavella, D. Ernst and D. Ruiz-Vega, Transient Stability of Power
Systems, 2000, Kluwer.
[5] B. Wang, K. Sun and X. Su, "A decoupling based direct method for
power system transient stability analysis," IEEE PESGM, Denver, CO,
2015, pp. 1-5.
[6] R. Barrio, M. Rodriguez, A. Abad, F. Blesa, "Breaking the limits: the
Taylor series method," Applied Mathematics and Computation, vol.
217, no., pp. 7940-7954, 2011.
[7] N. Duan, K. Sun, "Power system simulation using the multi-stage
Adomian decomposition method," IEEE Transactions on Power
Systems, in press.
[8] N. Duan, K. Sun, "Application of the Adomian decomposition method
for semi-analytic solutions of power system differential algebraic
equations," PowerTech, Eindhoven, 2015
[9] S. Wang, G. Radman, G. Zheng and W. Gao, "Analytical trajectory
extrapolation for power systems," IEEE PES General Meeting, San
Diego, CA, 2011, pp. 1-7.
[10] B. Wang, K. Sun, "Power system differential-algebraic
equations," arXiv prepoint arXiv:1512.05185, 2015
[11] S. Maslennikov, B. Wang, Q. Zhang, et al, "A Test Cases Library for
Fig. 5. WECC 179-bus power system Methods Locating the Sources of Sustained Oscillations," IEEE PES
General Meeting, Boston, MA, July 17-21, 2016.

TABLE III LARGE SYSTEM STABILITY AND CCT RANKING TABLE

Tripped Fault Ranking by ΔVn Ranking by CCT
Line Bus ΔVn CCT /s
31-80 80 1 -1.00 3 .043
24-25 24 2 -0.98 1 .026
22-23 23 3 -0.79 2 .035
114-171 171 4 -0.42 4 .049
115-127 127 5 -0.01 5 .071
130-131 130 6 1.81 6 .095
108-133 108 7 5.71 15 .392
14-21 21 8 5.92 10 .225
Fig. 6. Example one (0.1 second fault-on simulation) 19-25 19 9 14.15 11 .259
83-172 172 10 15.98 7 .102
V. CONCLUSIONS 104-135 104 11 16.35 17 .511
This paper proposes a time series based semi-analytical 48-55 55 12 16.94 18 .584
solution (SAS) that can be used to simulate the fault-on 136-152 136 13 19.61 16 .428
trajectory of power systems. The resulting trajectories are 41-58 41 14 29.30 19 .757
compared to the 4th-order Runge-Kutta method and prove to 49-64 49 15 33.30 20 1.228
be fairly accurate. Then, the last point from the fault-on 69-72 69 16 33.54 8 .138
trajectory is fed into a decoupling based direct method to 82-87 82 17 69.04 9 .178
determine the stability of the system under the tested fault. 115-127 115 18 129.19 12 .286
Future work includes the applications of the SAS to other 111-173 173 19 176.95 14 .376
direct methods, e.g. BCU, and the post-fault simulation. 82-91 91 20 247.32 13 .315