EQUAL AREA CRITERION

One of the methods for solving the swing

equation is by graphical integration, which is
called equal area criterion. This method is
applicable only for a single-machine-infinite-
bus or a two-machine system.

Assumptions
1. Pm does not change during the swing
2. Pe can be obtained from a steady state
solution of the system
3. Damping powers are negligible
4. Machines are represented by subtransient
reactances

Illustrative Example
For one-machine-infinite bus system, if a
sudden change of Pm occurs, find the
conditions for stability.

Machine Line Infinite bus

1
Solution
Pe

b A2
c e
Pm2
A1

Pm1
a

δ0 δs δm 90
o δmax δ

9 Initially, the operating point was a and the

rotor angle was δ = δ0 where Pm1 = Pe.
Hence, at this point Pm1 – Pe = Pa = 0
9 Finally, the operating point is b and the
rotor angle is δ = δs where Pm2 = Pe.
Hence, at this point Pm2 – Pe = Pa = 0
9 Pe cannot change instantaneously. Hence,
the difference Pm2 – Pe = Pa will start
accelerating the rotor. This will continue
until δ reaches and passes δs
9 As δ becomes greater than δs, Pe becomes
greater than Pm2 and the rotor starts
decelerating. This will continue until δ
reaches δm

2
9 δ starts decreasing and oscillating around
until it finally reaches the new steady state
δs
9 The excess energy stored in the inertia
during acceleration is shown as A1 = area
abc
9 The excess energy stored in the inertia
during acceleration is shown as A1 = area
abc. This energy is returned back during
deceleration as A2 = area bde
9 For net energy be zero A1 must equal A2

Mathematically
d 2δ
M 2 = Pa = Pm 2 − Pe
dt
Multiplying both sides by 2dδ/Mdt
d 2δ dδ 2 dδ
2 2 = Pa
dt dt M dt
The left-hand side can be written as
d ⎡⎛ dδ ⎞ 2 ⎤ 2 dδ
⎢⎜ ⎟ ⎥= Pa
dt ⎢⎣⎝ dt ⎠ ⎥⎦ M dt
Multiplying both sides by dt and integrating
2 δ
⎛ dδ ⎞ 2
⎜ ⎟ =
⎝ dt ⎠ M ∫δ P dδ
a
0

3
i.e.
12
dδ ⎛⎜ 2 ⎞
δ
=
dt ⎜⎝ M ∫δ Pa dδ ⎟⎟
0 ⎠
When the machine reaches the new steady
state δ is constant and dδ/dt = 0. Hence,
δm δm δs δm

∫δ P dδ = δ∫ ( P
a m2 − Pe ) dδ = ∫ ( Pm 2 − Pe ) dδ + ∫ ( Pm 2 − Pe )dδ = 0
δ0 δs
0 0

δs δm

∫δ ( P m2 − Pe ) dδ − ∫ ( Pe − Pm 2 ) dδ = A1 − A2 = 0
δs
0

Hence, the condition for the machine to reach

new steady state is that A1 = A2
Substituting Pe = Pmax sin δ
δs

A1 = Pm 2 (δ s − δ 0 ) − ∫ Pmax sin δ dδ
δ0
δm

A2 = ∫δ P max sin δ dδ − Pm 2 (δ m − δ s )
s

Note: δ0=sin-1(Pm1/Pmax) and δs=sin-1(Pm2/Pmax)

Evaluating the integrals and substituting
Pmax P
δm + cos δ m = δ 0 + max cos δ 0
Pm 2 Pm 2
which can be solved by trial and error for δm

4
Comments

9 Although δ is in degrees in the P-δ curve,

in all equations involved in the derivation
the angles are in radians
9 The maximum value of δm for stability is
δmax = 180˚ – δs.
9 If the operation were at the linear portion
of the P-δ curve, then for A1 = A2, δm = 2 δs
– δ0. Hence, the selection of the first trial
value of δm can be made wisely

5
 Example
T1 Line 1 T2
G Line 2
∞
a b

System data is as follows:

f=60 Hz, H = 3s, Pm=0.8 pu, EG=E∞=1.0 pu
+ −
G: x = x = x = j 0.25 pu , x = j 0.1 pu
'' 0

+ −
T1 & T2: x = x = x = j 0.1 pu
0

+ −
L1 & L2: x = x = j 0.4 pu , and x = j1.2 pu
0

(a) Obtain δm for a sudden loading of Pm = 1.3 pu

(b) Analyze the system stability if a three-phase fault
at point a occurs while the system operates at its
initial state. Let the fault clear itself at a time tc
which corresponds to δc = 70˚. Find also the
critical clearing angle δcc
(c) Repeat (b) for the fault at point b cleared at δc =
70˚ by removing the faulted line. Find also the
critical clearing angle δcc
(d) Consider a L-G fault at point a with the original
initial voltages, with Pm = 0.9 pu. Study the system
stability with the following switching operations:
1. The L-G fault of phase a occurs at t = 0
2. Circuit breakers at both ends of line 2 open at
δc1 = 60˚, stay open while δ advances 10˚, and
then reclose
3. Since the fault is not cleared, the circuit
breakers reopen after δ advances another 20˚,
and then remain open

6
Solution of Part (a)
The +ve sequence reactance diagram (symmetrical
fault) is as shown
0.1 0.4 0.1

0.25
a 0.4
E1 E 2

X 12 = 0 .6 5 p u
E1 E 2 1
Pe = sin δ = sin δ = 1.538 sin δ
X 12 0.65
⎛ Pm1 ⎞
⎟⎟ = sin −1 ⎛⎜
0.8 ⎞
δ 0 = sin −1 ⎜⎜ ⎟ = 31.33 = 0.547 rad
D

⎝ Pmax ⎠ ⎝ 1.538 ⎠
⎛ Pm 2 ⎞
⎟⎟ = sin −1 ⎛⎜
1.3 ⎞
δ s = sin −1 ⎜⎜ ⎟ = 57.67 = 1.007 rad
D

⎝ Pmax ⎠ ⎝ 1.538 ⎠
Pe

1.538

A2
1.3
A1

0.8

δ0 δs δm 90
o δmax δ

For stability A1 = A2
δs δm
Pm 2 (δ s − δ 0 ) − ∫ Pmax sin δ dδ = ∫δ P max sin δ dδ − Pm 2 (δ m − δ s )
δ0 s

7
 Pmax P
δm + cos δ m = δ 0 + max cos δ 0
Pm 2 Pm 2
δ m + 1.183 cos δ m = 1.558
Note: 2δs – δ0 < δm < 180˚– δs, i.e. 84˚ < δm < 122˚.
Solve by trial and error

δm δm + 1.183 cos δm
90 1.571
100 1.540
95 1.555
94 1.558

Hence, δm = 94˚.

Solution of Part (b)

During a 3-phase fault at point a, the electrical power
transfer during the fault is zero
Pe

1.538

A2

0.8

A1

δ0 δc δm δ

δ 0 = 31.33D = 0.547 rad , δ c = 70D = 1.222 rad

8
For stability A1 = A2
δm
Pm (δ c − δ 0 ) = ∫δ P max sin δ dδ − Pm (δ m − δ c )
c

0.8(1.222 − 0.547) = −1.538(cos δ m − cos 70D ) − 0.8(δ m − 1.222)

δ m + 1.9225 cos δ m = 1.204
Note: 2δc – δ0 < δm < 180˚– δ0, i.e. 110˚ < δm < 150˚.
Solve by trial and error

δm δm + 1.183 cos δm
120 1.133
110 1.262
115 1.195
114 1.207

Hence, δm ≈ 114˚.
To find the critical clearing angle δcc the power-angle
curve is shown as
Pe

1.538

A2

0.8

A1

δ0 δcc δmax
δ

Here, we have
δ max = 180D − δ 0 = 148.67D = 2.595 rad , cos δ max = −0.854

9
For stability A1 = A2
δ max
Pm (δ cc − δ 0 ) = ∫δ P max sin δ dδ − Pm (δ max − δ cc )
cc

0.8δ cc − 0.437 = −1.538( −0.854 − cos δ cc ) − 0.8( 2.595 − δ cc )

1.538 cos δ cc = 0.326
Hence, δcc ≈ 77.8˚.

Solution of Part (c)

During a 3-phase fault at point b, the reactance diagram
is as shown
0.1 0.1 0.45 0.2
0.4

0.25
b
0.2 0.2
E1 0.05 E2
E1 E2
⎯∆⎯
⎯
−Y
→
2.45

⎯Y⎯
⎯
−∆
→

Before the fault the power transfer is

Pebf = 1.538 sin δ
Here, during the fault the power transfer is
E1 E2 1
Pedf = sin δ = sin δ = 0.408 sin δ
X 12 2.45
After the fault with line 2 removed
X12 = 0.25 + 0.2 + 0.4 + 0.1 =0.85 pu
Hence, the power transfer is

10
 E1 E2 1
Peaf = sin δ = sin δ = 1.176 sin δ
X 12 0.85
The power-angle characteristics before, during, and
after the fault are shown as
Pe

Pe before fault
1.538

1.176 Pe after fault

A2
0.8
A1

0.408 Pe during fault

δ0 δc δm δmax
δ

δ 0 = 31.33D = 0.547 rad , δ c = 70D = 1.222 rad

⎛ 0.8 ⎞
δ max = 180D − sin −1 ⎜ ⎟ = 137.1 = 2.393 rad
D

⎝ 1.176 ⎠
For stability A1 = A2
δc δm
Pm (δ c − δ 0 ) − ∫ 0.408 sin δ dδ = ∫ 1.176 sin δ dδ − Pm (δ m − δ c )
δ0 δc

0.331 = 1.380 − 0.8δ m − 1.176 cos δ m

δ m + 1.47 cos δ m = 1.311
Solve by trial and error

δm δm + 1.183 cos δm
120 1.159
130 1.324

11
 135 1.317
137 1.316

Hence, δm > δmax, and the system is unstable.

To find the critical clearing angle replace δc by δcc and

δm by δmax, then for stability A1 = A2
δ cc δ max

Pm (δ cc − δ 0 ) − ∫ 0.408 sin δ dδ = ∫δ 1.176 sin δ dδ − P (δ

m max − δ cc )
δ0 cc

− Pmδ 0 + 0.408(cos δ cc − cos δ 0 ) = −1.176(cos δ max − cos δ cc ) − Pmδ max

Substituting and solving for δcc yields δcc = 69.7˚

Solution of Part (d)

Before the fault the power transfer is
⎛ 0.9 ⎞
Pebf = 1.538 sin δ and δ 0 = sin −1 ⎜ ⎟ = 35.8 = 0.625 rad
D

⎝ 1.538 ⎠
During the L-G fault, the sequence networks are as
shown

12
 0.1 0.4 0.1

a
0.25
0.4
E1 E 2
+ve

0.1 0.4 0.1

0.25
0.4 -ve

0.1 1.2 0.1

0.1
1.2 zero

0.35 × 0.3 0 .2 × 0 .7
Z− = j = j 0.162 Z0 = j = j 0.156
0.65 0 .9
Z − + Z 0 = j 0.317
0.35 0.3 0.981

E1 0.317 E2
⎯Y⎯
⎯
−∆
→

Here, X12 = 0.918 pu

E1 E2 1
Pedf = sin δ = sin δ = 1.02 sin δ
X 12 0.918
After the fault with line 2 removed
X12 = 0.25 + 0.2 + 0.4 + 0.1 =0.85 pu
Hence, the power transfer is
E1 E2 1
Peaf = sin δ = sin δ = 1.176 sin δ
X 12 0.85
The power-angle characteristics for the switching
sequences are shown as
13
 Pe

Pe before fault
1.538

1.176 Pe after fault

1.02
Pe during fault
0.9

δ0 δc1 δc2 δ δm δmax

c3 δ

For, δ c1 = 60D = 1.047 rad , δ c 2 = 70 = 1.222 rad , δ c 3 = 90 = 1.571 rad ,

D D

we have
δ c1
A1 = Pm (δ c1 − δ 0 ) − ∫ 1.02 sin δ dδ = 0.0625
δ0
δc2 δ c3
A2 = ∫ 1.176 sin δ dδ + ∫ 1.02 sin δ dδ
δ c1 δc2
δm
+ ∫ 1.176 sin δ dδ −Pm (δ m − δ c1 ) = 1.477 − 1.176 cos δ m − 0.9δ m
δ c3

For stability A1 = A2, i.e.

0.0625 = 1.477 − 0.9δ m − 1.176 cos δ m
δ m + 1.307 cos δ m = 1.571
Trial and error solution gives δm ≈ 90.05˚.
⎛ 0.9 ⎞
δ max = 180D − sin −1 ⎜ ⎟ = 130.1
D

⎝ 1.176 ⎠
Hence, δm < δmax, and the system is stable.

14
 NUMERICAL SOLUTIONS

Many techniques exist for solving non-linear

ordinary differential equations (ODE’s). The
widely used techniques in power system
stability studies are trapezoidal method and
Runge-Kutta method. For the purpose of this
discussion we shall assume that equations are
of the form
•

X = F ( X (t ), t ) → n first order DE
X (t 0 ) = X 0 → initial conditions
where X is n-dimensional.

1. Trapezoidal Method
•

f n +1
x(t )

fn

t
t t+h
From
•

X = F ( X (t ), t )

15
t +h • t +h

∫ X (t )dt = ∫ F ( X (t ), t )dt
t t

X (t + h) − X (t ) =
h
2
[
F ( X (t + h), t ) + F ( X (t ), t ) ]
Where area under curve is approximated by a
trapezoid. In general
X n+1
h
(
= X n + F n+1 + F n
2
)
Comments
1. Trapezoidal method is very popular for
stability studies due to its numerical
stability
2. Trapezoidal method has 2nd order accuracy

2. Runge-Kutta Method

Several forms of the Runge-Kutta method

exist. The commonly used one is the fourth-
order method. In this method, evaluate
(1)
K = F ( X n , tn )
( 2) h (1) h
K = F ( X n + K , tn + )
2 2
( 3) h ( 2) h
K = F ( X n + K , tn + )
2 2

16
 ( 4) ( 3)
K = F ( X n + hK , t n + h)
Then,
X n +1 (
h (1) ( 2) ( 3)
= X n + K + 2K + 2K + K
6
( 4)
)
Comments
1. Runge-Kutta method is very popular for
stability studies due to its excellent starting
properties
2. Four evaluations of the function per step is
required. This is the major disadvantage of
the method for high-dimensional power
applications

Application of The Trapezoidal Method

The trapezoidal rule can be applied to the

differential equations of the machine as
follows.
Recall
X n+1 = X n +
h
2
(
F n+1 + F n )
The swing equation is usually expressed in
two first order forms as

17
dω ω s
= (Pm − Pe (t ) )
dt 2 H
dδ
= ω − ωs
dt
Let
⎛ f1 ⎞
F = ⎜⎜ ⎟⎟
⎝ f2 ⎠
Then
ωs
f1, n = (Pm − Pe (n) )
2H
f 2, n = ωn − ωs
ωs
f1, n +1 = (Pm − Pe (n + 1) )
2H
f 2 , n +1 = ω n + 1 − ω s
Applying the trapezoidal rule,
h ωs
ω n +1 = ω n + (2 Pm − Pe (n) − Pe (n + 1) )
2 2H
and
h
δ n +1 = δ n + (ω n +1 + ω n − 2ω s )
2
Or by combining,
ωs h2
δ n +1 = δ n + h(ω n − ω s ) + (2 Pm − Pe (n) − Pe (n + 1) )
8H

18
Generally, the two resulting algebraic
equations can be written as
ωsh
ω n +1 = (Pm − Pe (n + 1) ) + α
4H
and
h
δ n +1 = (ω n +1 − 2ω s ) + β
2
where α and β are defined as
ωsh h
α = ωn + (Pm − Pe (n) ) β = δ n + ωn
4H 2
Observe: α and β are determined from
previous step

Example

T Line 1

G Line 2
∞
a

System data is as follows:

H = 5s, Pm=0.8 pu, V∞ = 1.0∠0 pu
D
f=60 Hz,
G: xd = j 0.2 pu and Vt = 1.05 pu
'

T: x = j 0.1 pu L1 & L2: x = j 0.4 pu

(a) Determine the swing equation of the
machine rotor

19
(b) If a three-phase fault at t = 0 with fault
impedance of Zf = j 0.1 pu to neutral occurs
at point a, develop the swing equation and
find the initial rotor acceleration.
(c) If the fault is cleared in 9 cycles by tripping
the faulted line, determine the swing
equation.
(d) Using the trapezoidal method with a time
step of 0.03 s, solve the derived swing
equations and study the system stability.

Solution of Part (a)

The reactance diagram is as shown
I 0.1 0.4

0.2
a 0.4

E∠δ Vt V∞

The transfer reactance between the Vt and the

infinite bus is 0.1 + 0.4/2 = 0.3 pu. Assume α
is the angle of Vt. Hence,
VtV∞ 1.05 × 1
Pe = sin α = sin α = 0.8
X 0.3
This gives α = sin-1 (0.8/3.5) = 13.21˚.

20
Therefore, Vt = 1.05∠13.21D
Vt − V∞ 1.05∠13.21D − 1.0∠0D
I= = = 0.8 + j 0.074 = 0.803∠ − 5.29D
jX j 0.3
E ' = Vt + jIxd' = 1.05∠13.21D + j 0.2 × 0.803∠ − 5.29 D = 1.111∠21.09 D pu
The total reactance = 0.2 + 0.1 + 0.4/2= 0.5 pu
Hence,
1.111 × 1
Pe = sin δ = 2.222 sin δ
0 .5
Then, the swing equation is given by
d 2δ ωs 377
= ( Pm − Pe ) = (0.8 − 2.222 sin δ )
dt 2 2 H 10
i.e.
d 2δ
2
= 37.7(0.8 − 2.222 sin δ )
dt

Solution of Part (b)

During a 3-phase fault at point b, the reactance
diagram is as shown

0.1 0.3 0.2

0.4

0.2 a
0.4 0.1
0.1

⎯
⎯→

21
 1.1

⎯Y⎯
⎯
−∆
→

1.111 × 1
Pe = sin δ = 1.010 sin δ
1 .1
Then, the swing equation is given by
d 2δ
2
= 37.7(0.8 − 1.010 sin δ )
dt
At the start of the fault δ0 = 21.09˚. Substitute
d 2δ
dt 2
= 37 . 7 (
0 . 8 − 1 . 010 sin 21 . 09)D
= 16 . 46 rad/s 2

Solution of Part (c)

The reactance diagram is as shown
0.1 0.4

0.2
a

The total reactance = 0.2 + 0.1 + 0.4 = 0.7 pu

Hence,
1.111 × 1
Pe = sin δ = 1.587 sin δ
0 .7

22
Then, the swing equation is given by
d 2δ
2
= 37.7(0.8 − 1.587 sin δ )
dt

For the above 3 parts we can write

⎧37.7(0.8 − 2.222 sin δ ) t<0
d δ ⎪
2

2
= ⎨ 37.7(0.8 − 1.010 sin δ ) 0 ≤ t < 0.15
dt ⎪ 37.7(0.8 − 1.587 sin δ )
⎩ t ≥ 0.15

Solution of Part (d)

The above swing equations can be written as:

Before fault
• •
ω = 37.7(0.8 − 2.222 sin δ ) δ = ω − 377
δ0 = 21.09˚=0.368 rad ω0 = ωs = 377
During fault
• •
ω = 37.7(0.8 − 1.010 sin δ ) δ = ω − 377
After fault
• •
ω = 37.7(0.8 − 1.587 sin δ ) δ = ω − 377

Let h = ∆t = 0.03 s. Then applying the

trapezoidal rule
h
β n = δ n + ωn
2
23
 ωsh
α n = ωn + (Pm − Pe (n) )
4H
ωh
ω n +1 = s (Pm − Pe (n + 1) ) + α n
4H
h
δ n +1 = (ω n +1 − 2ω s ) + β n
2
Substituting,
β n = δ n + 0.015ω n
α n = ω n + 0.5655(0.8 − Pe ( n) )
ω n +1 = 0.5655(0.8 − Pe ( n + 1) ) + α n
δ n +1 = 0.015 (ω n +1 − 754 ) + β n
Step-by-Step Solution
Initial condition
δ0 = 21.09˚=0.368 rad ω0 = ωs = 377
Pe (0) = 1.010 sin δ 0 = 0.3633
β 0 = δ 0 + 0.015ω 0 = 6.022
α 0 = ω 0 + 0.5655(0.8 − Pe (0) ) = 377.247

n = 0 ( t = 0.03 s)
Assume δ does not change much (∆t is small)
Pe (1) ≅ Pe (0) = 0.3633
Then,
ω1 = 0.5655(0.8 − Pe (1) ) + α 0 = 377.4939
δ 1 = 0.015 (ω1 − 754 ) + β 0 = 0.3744

24
n = 1 ( t = 0.06 s)
Pe ( 2) ≅ Pe (1) = 1.010 sin δ 1 = 0.36937
β1 = δ1 + 0.015ω1 = 6.0368
α1 = ω1 + 0.5655(0.8 − Pe (1) ) = 377.4374
Then,
ω 2 = 0.5655(0.8 − Pe ( 2) ) + α1 = 377.6809
δ 2 = 0.015 (ω 2 − 754 ) + β1 = 0.3920

For n ≥ 2 ( t ≥ 0.09 s), repeat as in previous

steps. Note: At t ≥ 0.15 s, obtain Pe form
Pe = 1.587 sin δ

25