Small Signal Analysis of Dynamic Systems

Contents
Small Signal Analysis of Dynamic Systems ................................................................................................. 3
Lineaerisation of Non-linear systems for small signal analysis ................................................................ 5
Eigenvalues and modal matrices ............................................................................................................... 9
Free motion of a dynamic system ........................................................................................................... 11
Example E.1 ............................................................................................................................................ 14
Participation Factor Analysis .................................................................................................................. 19
Small signal stability analysis of multi-machine system ........................................................................ 20
Introduction to Coherency .......................................................................................................................... 25
Network Transformation ......................................................................................................................... 26
Elimination of Nodes .......................................................................................................................... 26
Aggregation of Nodes ......................................................................................................................... 28
Aggregation of Generator Units .......................................................................................................... 32
Coherency ............................................................................................................................................... 33
Coherency Recognition ........................................................................................................................... 34
Example E.2. Single Machine Model ..................................................................................................... 39
Example E.3. Small Signal Model of IEEE 9 bus system ...................................................................... 41
Appendix 1 – Axis transformation of example E.1................................................................................. 46
Appendix 2 – MATLAB Program to calculate the modes of the IEEE 9 bus system ............................ 49
Appendix 3 – Example of three mass system ......................................................................................... 52
Small Signal Analysis of Dynamic Systems
A dynamic system like the power system can be represented by 𝑛𝑛 ordinary
differential equations such as;
𝑥𝑥̇ 1 = 𝑓𝑓1 (𝑥𝑥1 , 𝑥𝑥2 , ⋯ , 𝑥𝑥𝑛𝑛 ; 𝑢𝑢1 , 𝑢𝑢2 , ⋯ , 𝑢𝑢𝑟𝑟 ; 𝑡𝑡)
𝑥𝑥̇ 2 = 𝑓𝑓2 (𝑥𝑥1 , 𝑥𝑥2 , ⋯ , 𝑥𝑥𝑛𝑛 ; 𝑢𝑢1 , 𝑢𝑢2 , ⋯ , 𝑢𝑢𝑟𝑟 ; 𝑡𝑡)
(1)
⋮
𝑥𝑥̇ 𝑛𝑛 = 𝑓𝑓𝑛𝑛 (𝑥𝑥1 , 𝑥𝑥2 , ⋯ , 𝑥𝑥𝑛𝑛 ; 𝑢𝑢1 , 𝑢𝑢2 , ⋯ , 𝑢𝑢𝑟𝑟 ; 𝑡𝑡)
Where
𝑥𝑥𝑖𝑖 ith state variable
𝑟𝑟 Number of inputs
𝑛𝑛 Order of the system
The same system can be represented by matrix form as;
𝒙𝒙̇ = 𝒇𝒇(𝒙𝒙, 𝒖𝒖, 𝑡𝑡)
Where
𝑥𝑥1 𝑢𝑢1 𝑓𝑓1
𝑥𝑥2 𝑢𝑢2 𝑓𝑓
𝒙𝒙 = � ⋮ � 𝒖𝒖 = � ⋮ � 𝒇𝒇 = � 2 �
⋮
𝑥𝑥𝑛𝑛 𝑢𝑢𝑟𝑟 𝑓𝑓𝑛𝑛
If the derivatives of the state variables are not explicit functions of t then, the vector
equation can be expressed as;
𝒙𝒙̇ = 𝒇𝒇(𝒙𝒙, 𝒖𝒖) (2)
These state variables create the output of the system such as;
𝒚𝒚 = 𝒈𝒈(𝒙𝒙, 𝒖𝒖) (3)
Where
𝑦𝑦1 𝑔𝑔1
𝑦𝑦2 𝑔𝑔2
𝒚𝒚 = � ⋮ � 𝒈𝒈 = � ⋮ � (4)
𝑦𝑦𝑚𝑚 𝑔𝑔𝑚𝑚
We call this representation of a system as the state representation. Here it may be
noted that the rate of change of each state is given by the state itself and the inputs.
If the present state of the system and the present inputs are known using the set of
equations (2) and (3), we can calculate the rate of change of each state. If we
consider a sufficiently small time step so that the rate of change of states within this
time step is constant, then the state after this time can be calculated using the present
states and the rate of change of the present states. This way the total system becomes
predictable. But unless the set of equations in equation (2) and (3) are linear
differential equations, this set of equations cannot be solved analytically.
If the system is required to be solved analytically, it is required to linearise this
system.
Lineaerisation of Non-linear systems for small signal analysis

As a solution to this problem, we can linearize this system around an operating point
using the Taylor Series. The Taylor series tells if 𝑦𝑦 is a function of multiple variables
𝑥𝑥1 , 𝑥𝑥2 , ⋯ 𝑥𝑥𝑛𝑛 , then;
𝑦𝑦 = 𝑓𝑓(𝑥𝑥1 , 𝑥𝑥2 , ⋯ 𝑥𝑥𝑛𝑛 )
If each 𝑥𝑥𝑖𝑖 is given a small change ∆𝑥𝑥𝑖𝑖 from its original position of 𝑥𝑥𝑖𝑖,0 , then if 𝑦𝑦
experience a change of ∆𝑦𝑦
𝑦𝑦 = 𝑦𝑦0 + ∆𝑦𝑦

𝑦𝑦0 = 𝑓𝑓�𝑥𝑥1,0 , 𝑥𝑥2,0 , ⋯ 𝑥𝑥𝑛𝑛,0 �

And
𝑦𝑦 = 𝑓𝑓�𝑥𝑥1,0 + ∆𝑥𝑥1 , 𝑥𝑥2,0 + ∆𝑥𝑥2 , ⋯ 𝑥𝑥𝑛𝑛,0 + ∆𝑥𝑥𝑛𝑛 �

Taylor series states

𝑛𝑛
𝜕𝜕𝜕𝜕
𝑦𝑦 = 𝑦𝑦0 + � � ∆𝑥𝑥
𝜕𝜕𝑥𝑥𝑖𝑖 0 𝑖𝑖
𝑖𝑖=1

Or
𝑛𝑛
𝜕𝜕𝜕𝜕
∆𝑦𝑦 = � � ∆𝑥𝑥
𝜕𝜕𝑥𝑥𝑖𝑖 0 𝑖𝑖
𝑖𝑖=1

It is seen this using Taylor series it is possible to linearize a non-linear function

around an operating point.
 Example

Say 𝑦𝑦 = 𝑥𝑥12 + 𝑥𝑥1 𝑥𝑥2 + 𝑥𝑥22 . Let’s try to linearize this function around (𝑥𝑥1 , 𝑥𝑥2 ) =
(2,3).

𝜕𝜕𝜕𝜕
� = 2𝑥𝑥1 + 𝑥𝑥2 |(2,3) = 7
𝜕𝜕𝑥𝑥1 (2,3)

𝜕𝜕𝜕𝜕
� = 𝑥𝑥1 + 2𝑥𝑥2 |(2,3) = 8
𝜕𝜕𝑥𝑥2 (2,3)

Now the small signal liner equation around the point (𝑥𝑥1 , 𝑥𝑥2 ) = (2,3) can be
given as;

∆𝑦𝑦 = 7∆𝑥𝑥1 + 8∆𝑥𝑥2

We can confirm this by applying (𝑥𝑥1 , 𝑥𝑥2 ) = (2.01,3.05) where 𝑦𝑦 = 2.012 +

2.01 × 3.05 + 3.052 = 19.4731

Now we can calculate using the linearized equation as, 𝑦𝑦0 = 22 + 2 × 3 +

32 = 19, ∆𝑦𝑦 = 7 × 0.01 + 8 × .05 = 0.47 gives in 𝑦𝑦 = 19.47.

Now if we consider that the system is operating around a point with prefix 0, then;
𝒙𝒙̇ 𝟎𝟎 = 𝒇𝒇(𝒙𝒙𝟎𝟎 , 𝒖𝒖𝟎𝟎 ) (5)
𝒚𝒚𝟎𝟎 = 𝒈𝒈(𝒙𝒙𝟎𝟎 , 𝒖𝒖𝟎𝟎 ) (6)
Now we give a small perturbation to this state
𝒙𝒙 = 𝒙𝒙𝟎𝟎 + ∆𝒙𝒙 (7)
𝒖𝒖 = 𝒖𝒖𝟎𝟎 + ∆𝒖𝒖 (8)

𝒙𝒙̇ = 𝒇𝒇(𝒙𝒙, 𝒖𝒖)

𝒙𝒙̇ = 𝒇𝒇(𝒙𝒙𝟎𝟎 + ∆𝒙𝒙, 𝒖𝒖𝟎𝟎 + ∆𝒖𝒖) (9)

Now by the linearity of the system around the operating point we can write
𝒙𝒙̇ = 𝒙𝒙̇ 𝟎𝟎 + ∆𝒙𝒙̇
𝒙𝒙̇ = 𝒇𝒇(𝒙𝒙𝟎𝟎 , 𝒖𝒖𝟎𝟎 ) + ∆𝒙𝒙̇ (10)
The individual equations of this vector equation can be expressed as follows using
the Taylor series
𝑥𝑥𝚤𝚤̇ = 𝑥𝑥𝑖𝑖0 + ∆𝑥𝑥𝑖𝑖
𝑥𝑥𝚤𝚤̇ = 𝒇𝒇𝒊𝒊 [(𝒙𝒙𝟎𝟎 + ∆𝒙𝒙), (𝒖𝒖𝟎𝟎 + ∆𝒖𝒖)]
𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖
= 𝒇𝒇𝒊𝒊 (𝒙𝒙𝟎𝟎 , 𝒖𝒖𝟎𝟎 ) + ∆𝑥𝑥1 + ∆𝑥𝑥2 + ⋯ + ∆𝑥𝑥
𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥2 𝜕𝜕𝑥𝑥𝑛𝑛 𝑛𝑛
𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖
+ ∆𝑢𝑢1 + ∆𝑢𝑢2 + ⋯ + ∆𝑢𝑢
𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢2 𝜕𝜕𝑢𝑢𝑛𝑛 𝑟𝑟 (11)
Gives
𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖 𝜕𝜕𝑓𝑓𝑖𝑖
∆𝑥𝑥̇ 𝑖𝑖 = ∆𝑥𝑥1 + ∆𝑥𝑥2 + ⋯ + ∆𝑥𝑥𝑛𝑛 + ∆𝑢𝑢1 + ∆𝑢𝑢
𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥2 𝜕𝜕𝑥𝑥𝑛𝑛 𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢2 2
(12)
𝜕𝜕𝑓𝑓𝑖𝑖
+ ⋯+ ∆𝑢𝑢
𝜕𝜕𝑢𝑢𝑛𝑛 𝑟𝑟
Similarly
𝜕𝜕𝑔𝑔𝑖𝑖 𝜕𝜕𝑔𝑔𝑖𝑖 𝜕𝜕𝑔𝑔𝑖𝑖 𝜕𝜕𝑔𝑔𝑖𝑖 𝜕𝜕𝑔𝑔𝑖𝑖
∆𝑦𝑦𝑖𝑖 = ∆𝑥𝑥1 + ∆𝑥𝑥2 + ⋯ + ∆𝑥𝑥𝑛𝑛 + ∆𝑢𝑢1 + ∆𝑢𝑢
𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥2 𝜕𝜕𝑥𝑥𝑛𝑛 𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢2 2
(13)
𝜕𝜕𝑔𝑔𝑖𝑖
+ ⋯+ ∆𝑢𝑢
𝜕𝜕𝑢𝑢𝑛𝑛 𝑟𝑟
Now this can be expressed as a matrix equation
∆𝒙𝒙̇ 𝒏𝒏×1 = 𝑨𝑨𝑛𝑛×𝑛𝑛 ∆𝒙𝒙𝑛𝑛×1 + 𝑩𝑩𝑛𝑛×𝑟𝑟 ∆𝒖𝒖𝑟𝑟×1 (14)
∆𝒚𝒚𝑚𝑚×1 = 𝑪𝑪𝑚𝑚×𝑛𝑛 ∆𝒙𝒙𝑛𝑛×1 + 𝑫𝑫𝑚𝑚×𝑟𝑟 ∆𝒖𝒖𝑟𝑟×1 (15)

Where
𝜕𝜕𝑓𝑓1 𝜕𝜕𝑓𝑓1 𝜕𝜕𝑓𝑓1 𝜕𝜕𝑓𝑓1
⎡ ⋯ ⎤ ⎡ ⋯ ⎤
⎢𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥𝑛𝑛 ⎥ ⎢𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢𝑟𝑟 ⎥
𝐴𝐴𝑛𝑛×𝑛𝑛 =⎢ ⋮ ⋱ ⋮ ⎥ 𝐵𝐵𝑛𝑛×𝑟𝑟 =⎢ ⋮ ⋱ ⋮ ⎥
⎢ 𝜕𝜕𝑓𝑓𝑛𝑛 ⋯
𝜕𝜕𝑓𝑓𝑛𝑛 ⎥ ⎢ 𝜕𝜕𝑓𝑓𝑛𝑛 ⋯
𝜕𝜕𝑓𝑓𝑛𝑛 ⎥
⎣𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥𝑛𝑛 ⎦ ⎣𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢𝑟𝑟 ⎦
 𝜕𝜕𝑔𝑔1 𝜕𝜕𝑔𝑔1 𝜕𝜕𝑔𝑔1 𝜕𝜕𝑔𝑔1
⎡ ⋯ ⎤ ⎡ ⋯ ⎤
⎢ 𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥𝑛𝑛 ⎥ ⎢ 𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢𝑟𝑟 ⎥
𝐶𝐶𝑚𝑚×𝑛𝑛 =⎢ ⋮ ⋱ ⋮ ⎥ 𝐷𝐷𝑚𝑚×𝑟𝑟 =⎢ ⋮ ⋱ ⋮ ⎥
⎢𝜕𝜕𝑔𝑔𝑚𝑚 ⋯
𝜕𝜕𝑔𝑔𝑚𝑚 ⎥ ⎢𝜕𝜕𝑔𝑔𝑚𝑚 ⋯
𝜕𝜕𝑔𝑔𝑚𝑚 ⎥
⎣ 𝜕𝜕𝑥𝑥1 𝜕𝜕𝑥𝑥𝑛𝑛 ⎦ ⎣ 𝜕𝜕𝑢𝑢1 𝜕𝜕𝑢𝑢𝑛𝑛𝑛𝑛 ⎦

Laplace transform of equation (14) and (15) results in;

𝑠𝑠∆𝒙𝒙(𝑠𝑠) − ∆𝒙𝒙(0) = 𝑨𝑨 ∙ ∆𝒙𝒙(𝑠𝑠) + 𝑩𝑩 ∙ ∆𝒖𝒖(𝑠𝑠) (16)
∆𝒚𝒚(𝑠𝑠) = 𝐶𝐶 ∙ ∆𝒙𝒙(𝑠𝑠) + 𝑫𝑫 ∙ 𝒖𝒖(𝑠𝑠) (17)
Equation 16 can be expressed as;
(𝑠𝑠𝑰𝑰 − 𝑨𝑨)∆𝒙𝒙(𝑠𝑠) = ∆𝒙𝒙(0) + 𝑩𝑩∆𝒖𝒖(𝑠𝑠)
∆𝒙𝒙(𝑠𝑠) = (𝑠𝑠𝑰𝑰 − 𝑨𝑨)−1 [∆𝒙𝒙(0) + 𝑩𝑩∆𝒖𝒖(𝑠𝑠)]
adj(𝑠𝑠𝑰𝑰 − 𝑨𝑨) (18)
∆𝒙𝒙(𝑠𝑠) = [∆𝒙𝒙(0) + 𝑩𝑩∆𝒖𝒖(𝑠𝑠)]
det(𝑠𝑠𝑰𝑰 − 𝑨𝑨)

The poles of ∆𝑥𝑥(𝑠𝑠) are given by det(𝑠𝑠𝑰𝑰 − 𝑨𝑨) = 0. As this is a 𝑛𝑛 × 𝑛𝑛 matrix equation
there will be 𝑛𝑛 values which satisfy it.
Eigenvalues and modal matrices

Eigen matrix, left modal matrix and the right modal matrix of any 𝑛𝑛 × 𝑛𝑛 matrix 𝐴𝐴
can be defined as follows.
If
𝑨𝑨𝑛𝑛×𝑛𝑛 𝝓𝝓𝑛𝑛×1 = 𝜆𝜆𝝓𝝓𝑛𝑛×1 (19)

solutions of 𝜆𝜆 is called a Eigen values of matrix 𝐴𝐴. If we solve the above equation,
we can get
(𝑨𝑨 − 𝜆𝜆𝑰𝑰)𝝓𝝓 = 𝟎𝟎 (20)

In order to have non trivial solutions for 𝝓𝝓, det(𝑨𝑨 − 𝜆𝜆𝑰𝑰) = 0. As this results in an
equation of nth order for 𝜆𝜆, generally there are 𝑛𝑛 solutions for 𝜆𝜆 and there are 𝑛𝑛
column vectors for 𝝓𝝓.
i.e.
𝑨𝑨𝝓𝝓𝑖𝑖 = 𝜆𝜆𝑖𝑖 𝝓𝝓𝑖𝑖 where 𝑖𝑖 = 1 ⋯ 𝑛𝑛 (21)
𝜙𝜙1𝑖𝑖
𝜙𝜙
𝝓𝝓𝑖𝑖 = � 2𝑖𝑖 � (22)
⋮
𝜙𝜙𝑛𝑛𝑛𝑛

The right modal matrix 𝝓𝝓 is defined as

𝝓𝝓 = [𝝓𝝓1 𝝓𝝓2 𝝓𝝓3 ⋯ 𝝓𝝓𝑖𝑖 ⋯ 𝝓𝝓𝑛𝑛 ] (23)
𝜙𝜙11 𝜙𝜙12 𝜙𝜙13 ⋯ 𝜙𝜙1𝑖𝑖 ⋯ 𝜙𝜙1𝑛𝑛
⎡𝜙𝜙 𝜙𝜙22 𝜙𝜙23 ⋯ 𝜙𝜙2𝑖𝑖 ⋯ 𝜙𝜙2𝑛𝑛 ⎤
⎢ 21 ⎥
𝝓𝝓 = ⎢𝜙𝜙31 𝜙𝜙32 𝜙𝜙33 ⋯ 𝜙𝜙3𝑖𝑖 ⋯ 𝜙𝜙3𝑛𝑛 ⎥ (24)
⎢ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥
⎣𝜙𝜙𝑛𝑛1 𝜙𝜙𝑛𝑛2 𝜙𝜙𝑛𝑛3 ⋯ 𝜙𝜙𝑛𝑛𝑛𝑛 ⋯ 𝜙𝜙𝑛𝑛𝑛𝑛 ⎦

Similarly the left modal matrix can be defined as;

 𝝍𝝍𝑖𝑖 𝑨𝑨 = 𝜆𝜆𝑖𝑖 𝝍𝝍𝑖𝑖 where 𝑖𝑖 = 1 ⋯ 𝑛𝑛 (25)
𝝍𝝍𝑖𝑖 = [𝜓𝜓𝑖𝑖1 𝜓𝜓𝑖𝑖2 ⋯ 𝜓𝜓𝑖𝑖𝑖𝑖 ] (26)
𝝍𝝍 = [𝝍𝝍1𝑇𝑇 𝝍𝝍𝑇𝑇2 ⋯ 𝝍𝝍𝑇𝑇𝑖𝑖 ⋯ 𝝍𝝍𝑇𝑇𝑛𝑛 ]𝑇𝑇 (27)

𝜓𝜓11 𝜓𝜓12 𝜓𝜓13 ⋯ 𝜓𝜓1𝑛𝑛

⎡𝜓𝜓 𝜓𝜓22 𝜓𝜓23 ⋯ 𝜓𝜓2𝑛𝑛 ⎤
⎢ 21 ⎥
⎢𝜓𝜓31 𝜓𝜓32 𝜓𝜓33 ⋯ 𝜓𝜓3𝑛𝑛 ⎥
𝝍𝝍 = ⎢ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥ (28)
⎢ 𝜓𝜓𝑖𝑖1 𝜓𝜓𝑖𝑖2 𝜓𝜓𝑖𝑖3 ⋯ 𝜓𝜓𝑖𝑖𝑖𝑖 ⎥
⎢ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥
⎣𝜓𝜓𝑛𝑛1 𝜓𝜓𝑛𝑛2 𝜓𝜓𝑛𝑛3 ⋯ 𝜓𝜓𝑛𝑛𝑛𝑛 ⎦

The diagonal matrix of Eigen values if defined as;

𝜆𝜆1 0 0 ⋯ 0
⎡ ⎤
0 𝜆𝜆2 0 ⋯ 0
⎢ ⎥
𝚲𝚲 = ⎢ 0 0 𝜆𝜆3 ⋯ 0⎥ (29)
⎢⋮ ⋮ ⋮ ⋱ ⋮⎥
⎣0 0 0 ⋯ 𝜆𝜆𝑛𝑛 ⎦

The equation system 𝑨𝑨𝝓𝝓𝑖𝑖 = 𝜆𝜆𝑖𝑖 𝝓𝝓𝑖𝑖 can be re-written as;

𝑨𝑨𝑨𝑨 = 𝝓𝝓𝚲𝚲 (30.1)
−𝟏𝟏 −𝟏𝟏 −𝟏𝟏
Which can be written as 𝑨𝑨 = 𝝓𝝓𝚲𝚲𝝓𝝓 and 𝝓𝝓 𝑨𝑨 = 𝚲𝚲𝝓𝝓 by post and pre
−𝟏𝟏
multiplication by 𝝓𝝓 and 𝝓𝝓 consecutively.
Similarly, the equation system 𝝍𝝍𝑖𝑖 𝑨𝑨 = 𝜆𝜆𝑖𝑖 𝝍𝝍𝑖𝑖 can be written as
𝝍𝝍𝑨𝑨 = 𝚲𝚲𝝍𝝍 (30.2)
By the comparison of equation (30.1) and (30.2) we can derive
𝝍𝝍𝝍𝝍 = 𝐈𝐈 (31)
𝝍𝝍 = 𝝓𝝓−1 (32)
And by equation (30.1) we can further write,
𝝓𝝓−1 𝑨𝑨𝑨𝑨 = 𝚲𝚲 (33)
Free motion of a dynamic system

The free motion of a dynamic equation occurs when the input ∆𝑢𝑢 = 0 in equation
(14). At this state the state equation of the system given in equation (14) becomes;
∆𝒙𝒙̇ = 𝑨𝑨∆𝒙𝒙 (34)
In this equation, the rate of change of a state Δ𝑥𝑥𝚤𝚤̇ is a linear combination of all the
states as,
𝑥𝑥̇ 𝑖𝑖 = 𝑎𝑎𝑖𝑖1 Δ𝑥𝑥1 + 𝑎𝑎𝑖𝑖2 Δ𝑥𝑥2 + 𝑎𝑎𝑖𝑖3 Δ𝑥𝑥3 + ⋯ + 𝑎𝑎𝑖𝑖𝑖𝑖 Δ𝑥𝑥𝑛𝑛 (35)
It is not possible to derive a simple solution and it is not possible to identify the
contributions from each state variable to the final performance in this coupled form
of state representation.
Let us transform the physical state vector to a new state vector as,
Δ𝒙𝒙 = 𝝓𝝓𝝓𝝓 (36)
where 𝝓𝝓 is the right modal matrix of 𝑨𝑨.
Now
Δ𝒙𝒙̇ = 𝝓𝝓𝒛𝒛̇ (37)
Applying this to equation (34)
𝝓𝝓𝒛𝒛̇ = 𝑨𝑨𝑨𝑨𝑨𝑨 (38)
𝒛𝒛̇ = 𝝓𝝓−1 𝑨𝑨𝑨𝑨𝑨𝑨 (39)
𝒛𝒛̇ = 𝚲𝚲𝒛𝒛 (40)

It is seen that by this transformation we have achieved a set of decoupled variables

so that
𝑧𝑧̇1 𝜆𝜆1 0 ⋯ 0 ⋯ 0 𝑧𝑧1
⎡ 𝑧𝑧̇ ⎤ ⎡ 0 𝜆𝜆 ⋯ 0 ⋯ 0 ⎤ ⎡ 𝑧𝑧 ⎤
2
⎢ 2⎥ ⎢ 2
⎥⎢ ⋮ ⎥
⎢⋮⎥=⎢⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⎥
(41)
⎢ 𝑧𝑧̇𝑖𝑖 ⎥ ⎢ 0 0 ⋯ 𝜆𝜆𝑖𝑖 ⋯ 0 ⎥ ⎢ 𝑧𝑧𝑖𝑖 ⎥
⎢⋮⎥ ⎢⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣𝑧𝑧̇𝑛𝑛 ⎦ ⎣ 0 0 ⋯ 0 ⋯ 𝜆𝜆𝑛𝑛 ⎦ ⎣𝑧𝑧𝑛𝑛 ⎦
Has decoupled solutions 𝑧𝑧̇𝑖𝑖 = 𝜆𝜆𝑖𝑖 𝑧𝑧𝑖𝑖 where 𝑖𝑖 = 1,2, ⋯ 𝑛𝑛 (42)
Results in 𝑧𝑧𝑖𝑖 (𝑡𝑡) = 𝑧𝑧𝑖𝑖 (0)𝑒𝑒 𝑖𝑖 where 𝑖𝑖 = 1,2, ⋯ 𝑛𝑛
𝜆𝜆 𝑡𝑡 (43)
But from equation (36)
Δ𝒙𝒙 = 𝝓𝝓𝝓𝝓
∆𝑥𝑥𝑖𝑖 = 𝜙𝜙𝑖𝑖1 𝑧𝑧1 + 𝜙𝜙𝑖𝑖2 𝑧𝑧2 ⋯ + 𝜙𝜙𝑖𝑖𝑖𝑖 𝑧𝑧𝑛𝑛 (44)
∆𝑥𝑥𝑖𝑖 (𝑡𝑡) = 𝜙𝜙𝑖𝑖1 𝑧𝑧1 (𝑡𝑡) + 𝜙𝜙𝑖𝑖2 𝑧𝑧2 (𝑡𝑡) ⋯ + 𝜙𝜙𝑖𝑖𝑖𝑖 𝑧𝑧𝑛𝑛 (𝑡𝑡)
From equation (43)
∆𝑥𝑥𝑖𝑖 (𝑡𝑡) = 𝜙𝜙𝑖𝑖1 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝜙𝜙𝑖𝑖2 𝑧𝑧2 (0)𝑒𝑒 𝜆𝜆2𝑡𝑡 ⋯ + 𝜙𝜙𝑖𝑖𝑖𝑖 𝑧𝑧𝑛𝑛 (0)𝑒𝑒 𝜆𝜆𝑛𝑛𝑡𝑡 (45)

The response of each state variable to the disturbance is a combination of the

performance of all the Eigen values. The equation (45) tells us how individual Eigen
values, 𝜆𝜆1 , 𝜆𝜆2 , ⋯ , 𝜆𝜆𝑛𝑛 create the final response, ∆𝑥𝑥𝑖𝑖 (𝑡𝑡). It is seen that the individual
Eigen values connect to the final response of the system by the left Eigen matrix
elements, 𝜙𝜙𝑖𝑖1 , 𝜙𝜙𝑖𝑖2 , ⋯ , 𝜙𝜙𝑖𝑖𝑖𝑖 and the initial values of the transformed variables
𝑧𝑧1 (0), 𝑧𝑧2 (0), ⋯ , 𝑧𝑧𝑛𝑛 (0).
This means that the output of a system is generated by the stimulation of individual
modes. In this simple understanding, we can assume that Eigen modes are the
various mutually independent hidden oscillatory modes of the power system. Once
the power system is undergoing one mode of oscillation, all the state variables are
oscillating. Equation (45) tells us as linkage between each mode of oscillation with
each state variable.
Now it is required to find how the initial values of the transformed variable are
created by the initial values of the state variables. This can be achieved through the
right Eigen matrix.
Similarly from equation (36)
Δ𝒙𝒙 = 𝝓𝝓𝝓𝝓
𝒛𝒛 = 𝝓𝝓−1 Δ𝒙𝒙
𝒛𝒛 = 𝝍𝝍Δ𝒙𝒙 (46)
𝑧𝑧𝑖𝑖 (0) = 𝜓𝜓𝑖𝑖1 Δ𝑥𝑥1 (0) + 𝜓𝜓𝑖𝑖2 Δ𝑥𝑥2 (0) ⋯ + 𝜓𝜓𝑖𝑖𝑖𝑖 Δ𝑥𝑥𝑛𝑛 (0)
Let us call 𝑧𝑧𝑖𝑖 (0) = 𝐶𝐶𝑖𝑖
Then
∆𝑥𝑥𝑖𝑖 (𝑡𝑡) = 𝜙𝜙𝑖𝑖1 𝐶𝐶1 𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝜙𝜙𝑖𝑖2 𝐶𝐶2 𝑒𝑒 𝜆𝜆2𝑡𝑡 ⋯ + 𝜙𝜙𝑖𝑖𝑖𝑖 𝐶𝐶𝑛𝑛 𝑒𝑒 𝜆𝜆𝑛𝑛𝑡𝑡 (47)
 where 𝑧𝑧𝑗𝑗 (0) = 𝐶𝐶𝑗𝑗 = 𝜓𝜓𝑗𝑗1 Δ𝑥𝑥1 (0) + 𝜓𝜓𝑗𝑗2 Δ𝑥𝑥2 (0) ⋯ + 𝜓𝜓𝑗𝑗𝑗𝑗 Δ𝑥𝑥𝑛𝑛 (0) (48)

Now let’s assume that only Eigen mode 1 is exited so other than 𝑧𝑧1 (0), all the other
initial conditions are null, i.e. 𝑧𝑧2 (0) = 0, ⋯ 𝑧𝑧𝑖𝑖 (0) = 0, ⋯ 𝑧𝑧𝑛𝑛 (0) = 0 Then from
equation (45) we see that the first column of right Eigen matrix tells us how each
state is excited by Eigen mode 1.
∆𝑥𝑥1 (𝑡𝑡) = 𝜙𝜙11 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡
∆𝑥𝑥2 (𝑡𝑡) = 𝜙𝜙21 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡
∆𝑥𝑥3 (𝑡𝑡) = 𝜙𝜙31 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡
⋮
∆𝑥𝑥𝑛𝑛 (𝑡𝑡) = 𝜙𝜙𝑛𝑛1 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡
How can we excite only mode 1?
The answer is given in equation (48). The left Eigen matrix lets us how the initial
states excite each Eigen mode. According to equation 48, initial conditions of all the
state variables, Δ𝑥𝑥1 (0), Δ𝑥𝑥2 (0) ⋯ Δ𝑥𝑥𝑛𝑛 (0) contribute to the creation of 𝑧𝑧1 (0)
𝑧𝑧1 (0) = 𝐶𝐶1 = 𝜓𝜓11 Δ𝑥𝑥1 (0) + 𝜓𝜓12 Δ𝑥𝑥2 (0) ⋯ + 𝜓𝜓1𝑛𝑛 Δ𝑥𝑥𝑛𝑛 (0)
Therefore, it is seen that the time response of a system in free motion can be
described in terms of the Eigen values and the left and right Eigen vectors. The
excitation of each Eigen mode is determined by the linear combination of initial
condition of the state.
Example E.1

For example, let us take a two dimensional case of a second order system given by
∆𝑥𝑥̇ = 𝐴𝐴∆𝑥𝑥 where
−2.09 0.36
𝐴𝐴 = � �
−0.27 −0.91
The initial conditions of this free motion is given by;
1.5
𝑥𝑥(0) = � �
2
We can derive left and right Eigen matrices and the Eigen matrix as shown in the
below figure using MATLAB.
Therefore, as per the MATLAB calculations,
−0.97 −0.32
𝜙𝜙 = � �
−0.24 −0.95
−2 0
Λ=� �
0 −1
−1.12 0.38
𝜓𝜓 𝑇𝑇 = � �
0.29 −1.15
Using equation (47) we can calculate the how each eigen mode contributes to the
response of each state variable as;
𝑥𝑥1 (𝑡𝑡) = 𝐶𝐶1 𝜙𝜙11 𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝐶𝐶2 𝜙𝜙12 𝑒𝑒 𝜆𝜆2𝑡𝑡 (E1.1)

𝜙𝜙11 = −0.97 𝜙𝜙12 = −0.32

𝜆𝜆1 = −2 𝜆𝜆2 = −1

𝑥𝑥1 (𝑡𝑡) = 𝐶𝐶1 (−0.97)𝑒𝑒 −2𝑡𝑡 + 𝐶𝐶2 (−0.31)𝑒𝑒 −𝑡𝑡

Similarly using equation (47)
𝑥𝑥2 (𝑡𝑡) = 𝐶𝐶1 𝜙𝜙21 𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝐶𝐶2 𝜙𝜙22 𝑒𝑒 𝜆𝜆2𝑡𝑡 (E1.2)

𝑥𝑥2 (𝑡𝑡) = 𝐶𝐶1 (−0.24)𝑒𝑒 −2𝑡𝑡 + 𝐶𝐶2 (−0.95)𝑒𝑒 −𝑡𝑡

And
−1.12 0.38
𝜓𝜓 𝑇𝑇 = � �
0.29 −1.15
Using equation (48), now we can write how the initial states of the state variables
stimulate each eigen axis using the above right eigen matrix as;
𝐶𝐶1 = 𝜓𝜓11 Δ𝑥𝑥1 (0) + 𝜓𝜓12 Δ𝑥𝑥2 (0)
𝐶𝐶2 = 𝜓𝜓21 Δ𝑥𝑥1 (0) + 𝜓𝜓22 Δ𝑥𝑥2 (0)
𝐶𝐶1 = −1.12 × 1.5 + 0.38 × 2 = −0.99
𝐶𝐶2 = 0.29 × 1.5 − 1.15 × 2 = −1.87
Now these values can be applied to equations (E1.1) and (E1.2)
𝑥𝑥1 (𝑡𝑡) = (−0.99)(−0.97)𝑒𝑒 −2𝑡𝑡 + (−1.87)(−0.31)𝑒𝑒 −𝑡𝑡 (E1.3)
= 0.91𝑒𝑒 −2𝑡𝑡 + 0.58𝑒𝑒 −𝑡𝑡
𝑥𝑥2 (𝑡𝑡) = (−0.99)(−0.24)𝑒𝑒 −2𝑡𝑡 + (−1.87)(−0.95)𝑒𝑒 −𝑡𝑡 (E1.4)
= 0.23𝑒𝑒 −2𝑡𝑡 + 1.78𝑒𝑒 −𝑡𝑡

The two state variables can be plotted on a Cartesian space as shown below.

x2
Eigen Axis 2
2.5

C2Φ21,C2Φ22
≡ (0.58,1.78) A ≡ x1(0),x2(0)≡(1.5,2)
2 0

C 0.1

0.2
1.5
0.3
Movement of the second
eigen component 0.4
0.5
0.6
1
0.7

1
C1Φ11,C1Φ12
≡ (0.91,0.23)
0.5 Eigen Axis 1

2
B
3 Movement of the first eigen component
0 7654
0 0.5 1 1.5 2 2.5 x1

The point A represents the initial condition of the system. The two Eigen
components of the initial position is shown as point B and point C. The time
response of the system as given in equation (E1.3 and E1.4) shows that the first
Eigen component decays faster than the second Eigen component. Therefore, it is
seen that the time trajectory of the two states is asymptotic to second Eigen axis.
These trajectories for various initial values as shown in below figure shows that the
two states converges to zero as the two Eigen values are negative and they are
asymptotic to the axis corresponding to the lower Eigen value.
The axis transformation of related to this example is calculated in Appendix 1.
2.5

1.5

0.5

0
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
-0.5

-1

-1.5

-2

-2.5

The behaviour of each Eigen value in the total performance is determined by

individual values of the corresponding right Eigen vector. The Eigen values occur
either as a single or in conjugate pair. The mode shape of each mode means the
Argon diagram representation of the corresponding element of the Eigen vector.
For example, consider a state space model shown below.
𝑥𝑥 ̇ 3 8 𝑥𝑥1
� 1� = � �� �
𝑥𝑥̇ 2 −9 −5 𝑥𝑥2
Let us assume that the state 𝑥𝑥1 corresponds to Generator 1 and state 𝑥𝑥2 corresponds
to Generator 2.
Now the right Eigen matrix of the system can be calculated as;
−0.32 − 0.61i −0.32 + 0.61i 0.69∡118.1° 0.69∡ − 118.1°
𝜙𝜙 = � �=� �
0.73 0.73 0.73 0.73

−1 + 7.48i 0
Λ=� �
0 −1 − 7.48i

im
im
Φ21=-0-0.3234+j0.6050
G1

Φ12=0.7276 Φ22=0.7276

G2 Re G2 Re

G1
Φ11=-0-0.3234-j0.6050

λ1=-1+7.48i λ2=-1-7.48i

The system has an oscillatory behaviour in reference to the two conjugate Eigen
7.48 1
value pairs with a frequency of = 1.19Hz. The damping ratio is 𝜁𝜁 = =
2𝜋𝜋 √12 +7.482
0.018. The modal shapes shows that when the mode 1 is exited the two generators
swing 118 degrees out of phase.
Participation Factor Analysis

The excitation of one state variable of the system initially will excite the modal states
(transformed states) as described by equation (46) as given below.
𝒛𝒛 = 𝝍𝝍Δ𝒙𝒙
Let’s say only the initial condition pertaining to kth element for the state variables is
present, or only the kth state variable is disturbed, then;
0
⎡ 0 ⎤
⎢ ⎥
⎢ ⋮ ⎥
∆𝒙𝒙(𝟎𝟎) =
⎢∆𝑥𝑥𝑘𝑘 (0)⎥
⎢ ⋮ ⎥
⎣ 0 ⎦
Then
𝜓𝜓1𝑘𝑘 ∆𝑥𝑥𝑘𝑘 (0)
⎡ ⎤
𝜓𝜓2𝑘𝑘 ∆𝑥𝑥𝑘𝑘 (0)
⎢ ⎥
⎢ ⋮ ⎥
𝒛𝒛 =
𝜓𝜓 ∆𝑥𝑥
⎢ 𝑗𝑗𝑗𝑗 𝑘𝑘 (0) ⎥
⎢ ⋮ ⎥
⎣𝜓𝜓𝑛𝑛𝑛𝑛 ∆𝑥𝑥𝑘𝑘 (0)⎦
But from Equation (45),
∆𝑥𝑥𝑖𝑖 (𝑡𝑡) = 𝜙𝜙𝑖𝑖1 𝑧𝑧1 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝜙𝜙𝑖𝑖2 𝑧𝑧2 (0)𝑒𝑒 𝜆𝜆2𝑡𝑡 ⋯ + 𝜙𝜙𝑖𝑖𝑖𝑖 𝑧𝑧𝑛𝑛 (0)𝑒𝑒 𝜆𝜆𝑛𝑛𝑡𝑡
The effect of the same state variable ∆𝑥𝑥𝑘𝑘 (𝑡𝑡) can be given by;

∆𝑥𝑥𝑘𝑘 (𝑡𝑡) = 𝜙𝜙𝑘𝑘1 𝜓𝜓1𝑘𝑘 ∆𝑥𝑥𝑘𝑘 (0)𝑒𝑒 𝜆𝜆1𝑡𝑡 + 𝜙𝜙𝑘𝑘2 𝜓𝜓2𝑘𝑘 ∆𝑥𝑥𝑘𝑘 (0)𝑒𝑒 𝜆𝜆2𝑡𝑡 + ⋯ + 𝜙𝜙𝑘𝑘𝑘𝑘 𝜓𝜓𝑗𝑗𝑗𝑗 ∆𝑥𝑥𝑘𝑘 (0)𝑒𝑒 𝜆𝜆𝑗𝑗𝑡𝑡
+ ⋯ + 𝜙𝜙𝑘𝑘𝑘𝑘 𝜓𝜓𝑛𝑛𝑛𝑛 ∆𝑥𝑥𝑘𝑘 (0)𝑒𝑒 𝜆𝜆𝑛𝑛𝑡𝑡
∆𝑥𝑥𝑘𝑘 (𝑡𝑡) is oscillating in multiple modes. The jth mode oscillation is governed by the
coefficient 𝜙𝜙𝑘𝑘𝑘𝑘 𝜓𝜓𝑗𝑗𝑗𝑗 , .i.e. the participation of the initial condition excitation of the kth
state variable in the jth mode is given by 𝜙𝜙𝑖𝑖𝑖𝑖 𝜓𝜓𝑗𝑗𝑗𝑗 . This is called the participation
factor of kth state in jth mode.
𝑃𝑃𝑘𝑘𝑘𝑘 = 𝜙𝜙𝑘𝑘𝑘𝑘 𝜓𝜓𝑗𝑗𝑗𝑗
Small signal stability analysis of multi-machine system

Let us assume there are 𝑛𝑛 number of generators connected to a transmission network.

The swing equation of the 𝑖𝑖 𝑡𝑡ℎ generator can be written in standard notation as,
𝑑𝑑Δ𝜔𝜔𝑖𝑖
Δ𝑃𝑃𝑚𝑚𝑚𝑚 − Δ𝑃𝑃𝑒𝑒𝑒𝑒 = 𝑀𝑀𝑖𝑖 + 𝐷𝐷𝑖𝑖 Δ𝜔𝜔𝑖𝑖 (49)
𝑑𝑑𝑑𝑑

Assuming the governor response of the machine is much slower than the power
angle response, we can assume Δ𝑃𝑃𝑚𝑚𝑚𝑚 = 0. Then,
𝑑𝑑Δ𝜔𝜔𝑖𝑖 −Δ𝑃𝑃𝑒𝑒𝑒𝑒 𝐷𝐷𝑖𝑖
= − Δ𝜔𝜔𝑖𝑖 (50)
𝑑𝑑𝑑𝑑 𝑀𝑀𝑖𝑖 𝑀𝑀𝑖𝑖

And,
𝑑𝑑Δ𝛿𝛿𝑖𝑖
= Δ𝜔𝜔𝑖𝑖 (51)
𝑑𝑑𝑑𝑑

The transmission network can be modelled by the standard network equations as,
𝐼𝐼 ̅ 𝑌𝑌� ⋯ 𝑌𝑌�1𝑖𝑖 ⋯ 𝑌𝑌�1𝑛𝑛 𝑉𝑉�1
⎡ 1 ⎤ ⎡ 11 ⎤⎡ ⎤
⎢ ⋮̅ ⎥ ⎢ �⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎢ 𝐼𝐼𝑖𝑖 ⎥ = ⎢ 𝑌𝑌𝑖𝑖1 ⋯ 𝑌𝑌�𝑖𝑖𝑖𝑖 ⋯ 𝑌𝑌�𝑖𝑖𝑖𝑖 ⎥ ⎢ 𝑉𝑉�𝑖𝑖 ⎥
⎢⋮⎥ ⎢ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣𝐼𝐼𝑛𝑛̅ ⎦ ⎣𝑌𝑌�𝑛𝑛1 ⋯ 𝑌𝑌�𝑛𝑛𝑛𝑛 ⋯ 𝑌𝑌�𝑛𝑛𝑛𝑛 ⎦ ⎣𝑉𝑉�𝑛𝑛 ⎦
Current injection at any node can be written as
𝑛𝑛

𝐼𝐼𝑖𝑖̅ = 𝑌𝑌�𝑖𝑖𝑖𝑖 𝑉𝑉�𝑖𝑖 + � 𝑌𝑌�𝑖𝑖𝑖𝑖 𝑉𝑉�𝑗𝑗

𝑗𝑗=1;𝑗𝑗≠𝑖𝑖

Now we can write the voltages in their polar form as 𝑉𝑉�𝑖𝑖 = 𝑉𝑉𝑖𝑖 ∡𝛿𝛿𝑖𝑖 = 𝑉𝑉𝑖𝑖 𝑒𝑒 𝑗𝑗𝛿𝛿𝑖𝑖 but the
admittance is maintained in its Cartesian form as 𝑌𝑌�𝑖𝑖𝑖𝑖 = 𝐺𝐺𝑖𝑖𝑖𝑖 + 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 .
𝑛𝑛

𝐼𝐼𝑖𝑖̅ = (𝐺𝐺𝑖𝑖𝑖𝑖 + 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 )𝑉𝑉𝑖𝑖 𝑒𝑒 𝑗𝑗𝛿𝛿𝑖𝑖 + � �𝐺𝐺𝑖𝑖𝑖𝑖 + 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 � 𝑉𝑉𝑗𝑗 𝑒𝑒 𝑗𝑗𝛿𝛿𝑗𝑗 (52)
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖
Complex power injected in any node 𝑖𝑖 can be defined as,
𝑆𝑆𝑖𝑖̅ = 𝑃𝑃𝑖𝑖 + 𝑗𝑗𝑄𝑄𝑖𝑖
𝑆𝑆𝑖𝑖̅ = 𝑉𝑉�𝑖𝑖 𝐼𝐼𝑖𝑖∗̅
𝑛𝑛

= 𝑉𝑉𝑖𝑖 𝑒𝑒 𝑗𝑗𝛿𝛿𝑖𝑖 �(𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 )𝑉𝑉𝑖𝑖 𝑒𝑒 −𝑗𝑗𝛿𝛿𝑖𝑖 + � �𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 � 𝑉𝑉𝑗𝑗 𝑒𝑒 −𝑗𝑗𝛿𝛿𝑗𝑗 �
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖
𝑛𝑛

= (𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 )𝑉𝑉𝑖𝑖2 + 𝑉𝑉𝑖𝑖 � �𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 � 𝑉𝑉𝑗𝑗 𝑒𝑒 𝑗𝑗�𝛿𝛿𝑖𝑖 −𝛿𝛿𝑗𝑗�

𝑗𝑗=1;𝑗𝑗≠𝑖𝑖
𝑛𝑛

= (𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 )𝑉𝑉𝑖𝑖2 + � 𝑉𝑉𝑖𝑖 𝑉𝑉𝑗𝑗 �𝐺𝐺𝑖𝑖𝑖𝑖 − 𝑗𝑗𝐵𝐵𝑖𝑖𝑖𝑖 � �𝑐𝑐𝑐𝑐𝑐𝑐�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 � + 𝑗𝑗𝑗𝑗𝑗𝑗𝑗𝑗�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 ��
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖

We can separate the real and imaginary parts of the equation to derive expressions
for real and reactive power.
𝑛𝑛

𝑃𝑃𝑖𝑖 = 𝐺𝐺𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖2 + � 𝑉𝑉𝑖𝑖 𝑉𝑉𝑗𝑗 �𝐺𝐺𝑖𝑖𝑖𝑖 𝑐𝑐𝑐𝑐𝑐𝑐�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 � + 𝐵𝐵𝑖𝑖𝑖𝑖 𝑠𝑠𝑠𝑠𝑠𝑠�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 �� (53)
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖
𝑛𝑛

𝑄𝑄𝑖𝑖 = −𝐵𝐵𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖2 + � 𝑉𝑉𝑖𝑖 𝑉𝑉𝑗𝑗 �𝐺𝐺𝑖𝑖𝑖𝑖 𝑠𝑠𝑠𝑠𝑠𝑠�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 � − 𝐵𝐵𝑖𝑖𝑖𝑖 𝑐𝑐𝑐𝑐𝑐𝑐�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 �� (54)
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖

This is called the hybrid network equation.

It is observed that active and reactive power flows are non-linear functions of power
angles as 𝑃𝑃 = 𝑃𝑃(𝑉𝑉, 𝛿𝛿 ) and 𝑄𝑄 = 𝑄𝑄 (𝑉𝑉, 𝛿𝛿 ). We can linearize these equations around
an operating point to derive the change in active and reactive power injection in
every node due to marginal changes of power angle and voltage as,
 𝜕𝜕𝑃𝑃1 𝜕𝜕𝑃𝑃1 𝜕𝜕𝑃𝑃1 𝜕𝜕𝑃𝑃1
⎡ ⋯ ⋮ ⋯ ⎤
⎢ 𝜕𝜕𝛿𝛿1 𝜕𝜕𝛿𝛿𝑛𝑛 𝜕𝜕𝑉𝑉1 𝜕𝜕𝑉𝑉𝑛𝑛 ⎥
∆𝑃𝑃1 ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥ ∆𝛿𝛿1
⎡ ⋮ ⎤ ⎢ 𝜕𝜕𝑃𝑃𝑛𝑛 𝜕𝜕𝑃𝑃𝑛𝑛 𝜕𝜕𝑃𝑃𝑛𝑛 𝜕𝜕𝑃𝑃𝑛𝑛 ⎥ ⎡ ⋮ ⎤
⎢ ⋯ ⋮ ⋯
⎢ ∆𝑃𝑃 ⎥ ⎢ ⎢ ⎥
⎢ 𝑛𝑛 ⎥ 𝜕𝜕𝛿𝛿1 𝜕𝜕𝛿𝛿𝑛𝑛 𝜕𝜕𝑉𝑉1 𝜕𝜕𝑉𝑉𝑛𝑛 ⎥ ⎢∆𝛿𝛿𝑛𝑛 ⎥
⎢ ⋯ ⎥⎢ ⋯ ⎥
⎢ ⋯ ⎥= ⋯ ⋯ ⋯ ⋮ ⋯ ⋯
⎢ ∆𝑄𝑄1 ⎥ ⎢ 𝜕𝜕𝑄𝑄1 ⋯
𝜕𝜕𝑄𝑄1
⋮
𝜕𝜕𝑄𝑄1
⋯
𝜕𝜕𝑄𝑄1 ⎥ ⎢ ∆𝑉𝑉1 ⎥
⎢ ⋮ ⎥ ⎢ 𝜕𝜕𝛿𝛿1 𝜕𝜕𝛿𝛿𝑛𝑛 𝜕𝜕𝑉𝑉1
⎥
𝜕𝜕𝑉𝑉𝑛𝑛 ⎥ ⎢ ⋮ ⎥
⎢
⎣∆𝑄𝑄𝑛𝑛 ⎦ ⎢ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥ ⎣ ∆𝑉𝑉𝑛𝑛 ⎦
⎢𝜕𝜕𝑄𝑄𝑛𝑛 ⋯
𝜕𝜕𝑄𝑄𝑛𝑛
⋮
𝜕𝜕𝑄𝑄𝑛𝑛
⋯
𝜕𝜕𝑄𝑄𝑛𝑛 ⎥
⎣ 𝜕𝜕𝛿𝛿1 𝜕𝜕𝛿𝛿𝑛𝑛 𝜕𝜕𝑉𝑉1 𝜕𝜕𝑉𝑉𝑛𝑛 ⎦
This matrix can be partitioned into four sub matrices as
∆𝑷𝑷 𝑯𝑯 𝑴𝑴 ∆𝜹𝜹
� �=� �� �
∆𝑸𝑸 𝑵𝑵 𝑲𝑲 ∆𝑽𝑽
The elements of each matrix can now be derived by partial differentiation of equation
(53) and (54). For example, the elements of matrix 𝑯𝑯 can be derived as,
For 𝑖𝑖 ≠ 𝑗𝑗 𝜕𝜕𝑃𝑃𝑖𝑖
𝐻𝐻𝑖𝑖𝑖𝑖 = = 𝑉𝑉𝑖𝑖 𝑉𝑉𝑗𝑗 �𝐺𝐺𝑖𝑖𝑖𝑖 𝑠𝑠𝑠𝑠𝑠𝑠�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 � − 𝐵𝐵𝑖𝑖𝑖𝑖 𝑐𝑐𝑐𝑐𝑐𝑐�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 �� (55)
𝜕𝜕𝛿𝛿𝑗𝑗
𝑛𝑛
For 𝑖𝑖 = 𝑗𝑗 𝜕𝜕𝑃𝑃𝑖𝑖
𝐻𝐻𝑖𝑖𝑖𝑖 = = � 𝑉𝑉𝑖𝑖 𝑉𝑉𝑗𝑗 �−𝐺𝐺𝑖𝑖𝑖𝑖 𝑠𝑠𝑠𝑠𝑠𝑠�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 �
𝜕𝜕𝛿𝛿𝑖𝑖 (56)
𝑗𝑗=1;𝑗𝑗≠𝑖𝑖
+ 𝐵𝐵𝑖𝑖𝑖𝑖 𝑐𝑐𝑐𝑐𝑐𝑐�𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 ��

Similarly, we can derive the elements of M, N and K matrices as well.

Anyhow, it is interesting to note that in the H matrix the diagonal comprises of the
negative summation of all the off diagonal elements of a particular raw as, from
equation (55) and (56)
𝑛𝑛

� 𝐻𝐻𝑖𝑖𝑖𝑖 = −𝐻𝐻𝑖𝑖𝑖𝑖
𝑗𝑗=1,𝑗𝑗≠𝑖𝑖

Or we can express this as the summation of the raw elements of H matrix is zero.
 𝑛𝑛 𝑛𝑛
(57)
𝐻𝐻𝑖𝑖𝑖𝑖 + � 𝐻𝐻𝑖𝑖𝑖𝑖 = � 𝐻𝐻𝑖𝑖𝑖𝑖 = 0
𝑗𝑗=1,𝑗𝑗≠𝑖𝑖 𝑗𝑗=1

The transmission network can be defined incorporating the generator impedances

and the loads into the admittance matrix so that finally all the generator nodes will
appear as current injections, as shown in the below figure.

I1 L1 I1
I1
G1 G1

I2 L2 I2
G2 I2 G2

Ii Transmission Li Ii
Equivalent
G3 Ii Network G3
Network
-- --
-- --
In -- --
Ln In
Gn Gn

In such definition the, if the generator voltages are considered constant, then we can
assume ∆𝑽𝑽 = 𝟎𝟎 or the ∆𝑽𝑽 submatrix is a null matrix. This assumption will simplify
our power flow equation as,
∆𝑷𝑷 = 𝑯𝑯∆𝜹𝜹
This can be expressed as,
𝑛𝑛

∆𝑃𝑃𝑖𝑖 = � 𝐻𝐻𝑖𝑖𝑖𝑖 ∆𝛿𝛿𝑗𝑗 (58)

𝑗𝑗=1

Applying equation (58) to equation (50),

𝑛𝑛
𝑑𝑑Δ𝜔𝜔𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐷𝐷𝑖𝑖
=� ∆𝛿𝛿𝑗𝑗 − Δ𝜔𝜔𝑖𝑖 (60)
𝑑𝑑𝑑𝑑 𝑀𝑀𝑖𝑖 𝑀𝑀𝑖𝑖
𝑗𝑗=1

Now equation (51) and (60) can be expressed in matrix form as,
 0 0 0 0 ⋮ 1 0 ⋯ 0
⎡ 0 0 0 0 ⋮ 0 1 ⋯ 0 ⎤
∆𝛿𝛿̇ ⎢ ⎥ ∆𝛿𝛿1
⎡ 1⎤ ⎢ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮
⎥ ⎡ ∆𝛿𝛿2 ⎤
⎢ ∆𝛿𝛿2̇ ⎥ ⎢ 0 0 0 0 ⋮ 0 0 ⋯ 1 ⎥⎢ ⋯ ⎥
⎢ ⋯ ⎥ ⎢ − − − − ⋮ − − − − ⎥⎢ ⎥
⎢ ∆𝛿𝛿𝑛𝑛̇ ⎥ ⎢ 𝐻𝐻11 𝐻𝐻12 𝐻𝐻1𝑛𝑛 𝐷𝐷1 ⎥ ⎢ ∆𝛿𝛿𝑛𝑛 ⎥
⎢ ⋯ ⎥= − − ⋯ − ⋮ − 0 ⋯ 0 (61)
⎢ ⎥ ⎢ 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 ⎥⎢ ⋯ ⎥
⎢ ∆𝜔𝜔̇ 1 ⎥ ⎢ 𝐻𝐻21 𝐻𝐻22 𝐻𝐻2𝑛𝑛 𝐷𝐷2 ⎥ ⎢∆𝜔𝜔1 ⎥
⎢∆𝜔𝜔̇ 2 ⎥ ⎢− 𝑀𝑀2 − ⋯ − ⋮ 0 − ⋯ 0 ⎥ ⎢∆𝜔𝜔2 ⎥
𝑀𝑀2 𝑀𝑀2 𝑀𝑀2
⎢ ⋮ ⎥ ⎢ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥
⎥⎢ ⋮ ⎥
⎣∆𝜔𝜔̇ 𝑛𝑛 ⎦ ⎢ 𝐻𝐻𝑛𝑛1 𝐻𝐻𝑛𝑛2 𝐻𝐻𝑛𝑛𝑛𝑛 𝐷𝐷𝑛𝑛 ⎣∆𝜔𝜔𝑛𝑛 ⎦
⎢− − ⋯ − ⋮ 0 0 ⋯ − ⎥
⎣ 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 ⎦
Introduction to Coherency
A typical power system is so large that a power system analysis program may take
long time and many resources to completely model it. Therefore, usually a part of
the system called internal subsystem is modelled in detail while the remainder of the
system called external subsystem is represented by an equivalent system. The
division of the network to internal and external subsystem is graphically shown in
figure below. The division of the network to external and internal subsystem is done
based on the assumption that the disturbance occurs only inside the internal
subsystem. The replacement of the external subsystem by a reduced equivalent
therefore shall not affect the overall outcome of the analysis. The nodes between
the external subsystem and the internal subsystem are called the border nodes.
{L}

{B}

Internal External
{G}
Subsystem Subsystem

The equivalent of the external subsystem is created using measurements or by

reduction methods. Basically, the model reduction techniques widely used in power
system analysis are physical reduction and topological reduction. The replacement
of parallel operating generators and lines by physical equivalents is done in physical
reduction. The network equivalent models created by eliminating and aggregating
network nodes to reduce the external network is used in the topological reduction.
In the elimination of nodes the respective nodes completely vanish from the system
and in the aggregation of nodes, such nodes are replaced by one equivalent node.
The equivalent of the external subsystem created through the physical and
topological reduction results in a model comprising of standard elements such as
equivalent generators, equivalent lines and equivalent loads enabling the external
subsystem to be modelled using standard software.
Network Transformation

Elimination of Nodes
Let’s consider a network of 𝑚𝑚 + 𝑛𝑛 nodes where 𝑛𝑛 number of nodes is eliminated to
create a reduced network. The elimination of the node set {𝐸𝐸} from the total network
of nodes set {𝑅𝑅 ∪ 𝐸𝐸} results the node set {𝑅𝑅} remaining.
Therefore, |𝑅𝑅 ∪ 𝐸𝐸 | = 𝑚𝑚 + 𝑛𝑛, where |𝑅𝑅| = 𝑚𝑚 and |𝐸𝐸 | = 𝑛𝑛.
The elimination process shall ensure that the nodal voltages and currents of the
remaining node set {𝑅𝑅} does not change due to the elimination process.
{R} {E} {R}

Unreduced IR-ΔIR Reduced

IR
Network Network

m n m

The network equations in the split form can be written for the complete network as
𝐼𝐼 𝑌𝑌 𝑌𝑌𝑅𝑅𝑅𝑅 𝑉𝑉𝑅𝑅 (62)
� 𝑅𝑅 � = � 𝑅𝑅𝑅𝑅 �� �
𝐼𝐼𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝐸𝐸

Which gives
𝐼𝐼𝑅𝑅 = 𝑌𝑌𝑅𝑅𝑅𝑅 𝑉𝑉𝑅𝑅 + 𝑌𝑌𝑅𝑅𝑅𝑅 𝑉𝑉𝐸𝐸 (63)
𝐼𝐼𝐸𝐸 = 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝑅𝑅 + 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝐸𝐸 (64)
−1 −1
𝑉𝑉𝐸𝐸 = 𝑌𝑌𝐸𝐸𝐸𝐸 𝐼𝐼𝐸𝐸 − 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝑅𝑅
−1 −1
𝐼𝐼𝑅𝑅 = 𝑌𝑌𝑅𝑅𝑅𝑅 𝑉𝑉𝑅𝑅 + 𝑌𝑌𝑅𝑅𝑅𝑅 (𝑌𝑌𝐸𝐸𝐸𝐸 𝐼𝐼𝐸𝐸 − 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝑅𝑅 )
−1
𝐼𝐼𝑅𝑅 = (𝑌𝑌𝑅𝑅𝑅𝑅 − 𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 −1
𝑌𝑌𝐸𝐸𝐸𝐸 )𝑉𝑉𝑅𝑅 + 𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 𝐼𝐼𝐸𝐸 (65)
Results in

𝐼𝐼𝑅𝑅 −1
𝑌𝑌𝑅𝑅𝑅𝑅 − 𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 −1
𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 𝑉𝑉𝑅𝑅 (66)
� �=� −1 � � 𝐼𝐼 �
𝑉𝑉𝐸𝐸 −1
−𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 𝐸𝐸
This can be written as
𝐼𝐼 𝑌𝑌𝑅𝑅 𝐾𝐾𝐼𝐼 𝑉𝑉𝑅𝑅 (67)
� 𝑅𝑅 � = � −1 � � 𝐼𝐼 �
𝑉𝑉𝐸𝐸 𝐾𝐾𝑉𝑉 𝑌𝑌𝐸𝐸𝐸𝐸 𝐸𝐸

Where
−1
𝑌𝑌𝑅𝑅 = 𝑌𝑌𝑅𝑅𝑅𝑅 − 𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 (68)
−1
𝐾𝐾𝐼𝐼 = 𝑌𝑌𝑅𝑅𝑅𝑅 𝑌𝑌𝐸𝐸𝐸𝐸 (69)
−1
𝐾𝐾𝑉𝑉 = −𝑌𝑌𝐸𝐸𝐸𝐸 𝑌𝑌𝐸𝐸𝐸𝐸 (70)
From the first row of the matrix
[𝐼𝐼𝑅𝑅 ]𝑚𝑚×1 = [𝑌𝑌𝑅𝑅 ]𝑚𝑚×𝑚𝑚 [𝑉𝑉𝑅𝑅 ]𝑚𝑚×1 + [𝐾𝐾𝐼𝐼 ]𝑚𝑚×𝑛𝑛 [𝐼𝐼𝐸𝐸 ]𝑛𝑛×1 (71)
If 𝐾𝐾𝐼𝐼 𝐼𝐼𝐸𝐸 is defined as Δ𝐼𝐼𝑅𝑅
𝐼𝐼𝑅𝑅 = 𝑌𝑌𝑅𝑅 𝑉𝑉𝑅𝑅 + Δ𝐼𝐼𝑅𝑅 (72)
Any network is sufficiently described by the nodal currents and nodal voltages.
Therefore, the impedance matrix 𝑌𝑌𝑅𝑅 uniquely describes the retained nodes and the
links interlinking them. The matrix [𝐾𝐾𝐼𝐼 ]𝑚𝑚×𝑛𝑛 gives the currents referred from the
eliminated nodes {𝐸𝐸 } to the retained nodes {𝑅𝑅}. Each equivalent current is a
combination of the eliminated currents. However until such time Δ𝐼𝐼𝑅𝑅 is a non-zero
matrix, the elimination is not fully complete
The second row of the matrix can be expanded as
−1 ]
[𝑉𝑉𝐸𝐸 ]𝑛𝑛×1 = [𝐾𝐾𝑉𝑉 ]𝑛𝑛×𝑚𝑚 [𝑉𝑉𝑅𝑅 ]𝑚𝑚×1 + [𝑌𝑌𝐸𝐸𝐸𝐸 𝑛𝑛×𝑛𝑛 [𝐼𝐼𝐸𝐸 ]𝑛𝑛×1 (73)

Now, if the nodal power injection at each node in the elimination subset {𝐸𝐸} is
𝑆𝑆 ∗
represented by an equivalent shunt admittance 𝑌𝑌𝐸𝐸𝐸𝐸 = 𝑖𝑖� 2 with the appropriate sign
𝑉𝑉𝑖𝑖
to denote whether the power at the node is injected or absorbed.
 {R} {E} {R} {E}

Unreduced Unreduced
IR IR
Network Network
-YE1
-YE2
m n m n
-YEn

With such augmentation the diagonal elements of the 𝑌𝑌𝐸𝐸𝐸𝐸 can be added with these
admittances as,
𝑌𝑌𝐸𝐸𝐸𝐸11 − 𝑌𝑌𝐸𝐸1 𝑌𝑌𝐸𝐸𝐸𝐸12 ⋯ 𝑌𝑌𝐸𝐸𝐸𝐸𝐸𝐸1 (74)
𝑌𝑌𝐸𝐸𝐸𝐸21 𝑌𝑌𝐸𝐸𝐸𝐸22 − 𝑌𝑌𝐸𝐸2 ⋯ 𝑌𝑌𝐸𝐸𝐸𝐸𝐸𝐸2
𝑌𝑌𝐸𝐸𝐸𝐸,𝑛𝑛𝑛𝑛𝑛𝑛 =� �
⋮ ⋮ ⋱ ⋮
𝑌𝑌𝐸𝐸𝐸𝐸𝐸𝐸1 𝑌𝑌𝐸𝐸𝐸𝐸𝐸𝐸2 ⋯ 𝑌𝑌𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸 − 𝑌𝑌𝐸𝐸𝐸𝐸

so that the reduced model does not have equivalent currents, as, Δ𝐼𝐼𝑅𝑅 = 0.
Therefore,
𝐼𝐼𝑅𝑅 = 𝑌𝑌𝑅𝑅 𝑉𝑉𝑅𝑅 (75)

Aggregation of Nodes
Let’s assume that there are [R,A] cluster of generator nodes out of which the cluster
[A] are coherent and can be aggregated to a node a. However, in order to aggregate
the generators, it is essential to aggregate the generation nodes of the network in to
a single generator as shown in figure below.
{R} {A} {R}

Yfa a
Unreduced IR Reduced
Ia
IR Ia
Network Yif Network
Va Va

This aggregation is subjected to two constraints.

a. It shall not change the currents and voltages, 𝑰𝑰𝑹𝑹 and 𝑽𝑽𝑹𝑹
b. It shall not change the sum of the aggregation of real and reactive power
injection to the network, and therefore, 𝑆𝑆𝑎𝑎 = ∑𝑖𝑖∈[𝐴𝐴] 𝑆𝑆𝑖𝑖
The admittance matrix of the unreduced network can be written as,
𝐼𝐼1 𝑌𝑌11 ⋯ 𝑌𝑌1𝑟𝑟 𝑌𝑌1,𝑟𝑟+1 ⋯ 𝑌𝑌1,𝑟𝑟+𝑎𝑎 𝑉𝑉 (76)
⎡ ⋮ ⎤ ⎡ ⋮ ⎤⎡ 1 ⎤
⋱ ⋮ ⋮ ⋱ ⋮ ⋮
⎢ ⎥ ⎢ 𝑌𝑌 ⋯ 𝑌𝑌𝑟𝑟𝑟𝑟 𝑌𝑌𝑟𝑟,𝑟𝑟+1 ⋯
⎥
𝑌𝑌𝑟𝑟,𝑟𝑟+𝑎𝑎 ⎢ 𝑉𝑉𝑟𝑟 ⎥
⎢ 𝐼𝐼𝑟𝑟 ⎥ = ⎢ 𝑟𝑟1 ⎥⎢ ⎥
⎢𝐼𝐼𝑟𝑟+1 ⎥ ⎢𝑌𝑌𝑟𝑟+1,1 ⋯ 𝑌𝑌𝑟𝑟+1,𝑟𝑟 𝑌𝑌𝑟𝑟+1,𝑟𝑟+1 ⋯ 𝑌𝑌𝑟𝑟+1,𝑟𝑟+𝑎𝑎 ⎥ ⎢𝑉𝑉𝑟𝑟+1 ⎥
⎢ ⋮ ⎥ ⎢ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣𝐼𝐼𝑟𝑟+𝑎𝑎 ⎦ ⎣𝑌𝑌𝑟𝑟+𝑎𝑎,1 ⋯ 𝑌𝑌𝑟𝑟+𝑎𝑎,𝑟𝑟 𝑌𝑌𝑟𝑟+𝑎𝑎,𝑟𝑟+1 ⋯ 𝑌𝑌𝑟𝑟+𝑎𝑎,𝑟𝑟+𝑎𝑎 ⎦ ⎣𝑉𝑉𝑟𝑟+𝑎𝑎 ⎦
If this is written as an equation of sub matrices corresponding to the retained nodes,
[R] and aggregate nodes, [A], then,
𝑰𝑰 𝒀𝒀 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑹𝑹 (77)
� 𝑹𝑹 � = � 𝑹𝑹𝑹𝑹 �� �
𝑰𝑰𝑨𝑨 𝒀𝒀𝑨𝑨𝑨𝑨 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨
Now the admittance matrix of the reduced network can be written as,
𝐼𝐼1 𝑌𝑌11 ⋯ 𝑌𝑌1𝑟𝑟 𝑌𝑌1,𝑎𝑎 𝑉𝑉1 (78)
⋮ ⋮ ⋱ ⋮ ⋮ ⋮
� � = � 𝑌𝑌 ⋯ 𝑌𝑌𝑟𝑟𝑟𝑟 𝑌𝑌𝑟𝑟,𝑎𝑎 𝑉𝑉𝑟𝑟 �
� �
𝐼𝐼𝑟𝑟 𝑟𝑟1
𝐼𝐼𝑎𝑎 𝑌𝑌𝑎𝑎,1 ⋯ 𝑌𝑌𝑎𝑎,𝑟𝑟 𝑌𝑌𝑎𝑎,𝑎𝑎 𝑉𝑉𝑎𝑎
Which can be written as
𝑰𝑰 𝒀𝒀 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑹𝑹 (79)
� 𝑹𝑹 � = � 𝑹𝑹𝑹𝑹 �� �
𝐼𝐼𝑎𝑎 𝒀𝒀𝒂𝒂𝒂𝒂 𝑌𝑌𝑎𝑎𝑎𝑎 𝑉𝑉𝑎𝑎
Where 𝒀𝒀𝑹𝑹𝑹𝑹 is a column matrix and 𝒀𝒀𝑨𝑨𝑨𝑨 is a row matrix. 𝐼𝐼𝑎𝑎 and 𝑉𝑉𝑎𝑎 and values.
Now we can equate the two matrix equations as
𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑹𝑹 + 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑨𝑨 = 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑹𝑹 + 𝒀𝒀𝑹𝑹𝑹𝑹 𝑉𝑉𝑎𝑎 (80)
In order for this condition to be satisfied
𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑨𝑨 = 𝒀𝒀𝑹𝑹𝑹𝑹 𝑉𝑉𝑎𝑎 (81)
or
𝒀𝒀𝑹𝑹𝑹𝑹 = 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑨𝑨 𝑉𝑉𝑎𝑎−1 (82)
If we consider 𝑽𝑽𝑨𝑨 𝑉𝑉𝑎𝑎−1 as 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕
Then we can write
 𝒀𝒀𝑹𝑹𝑹𝑹 = 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (83)

Now the second condition can be elaborated as,

𝑽𝑽𝑻𝑻𝑨𝑨 𝑰𝑰𝑨𝑨∗ = 𝑉𝑉𝑎𝑎 𝐼𝐼𝑎𝑎∗ (84)
Now we can equate the second row of the two network equations,
𝑰𝑰𝑨𝑨 = 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑹𝑹 + 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨 (85)

𝐼𝐼𝑎𝑎 = 𝒀𝒀𝒂𝒂𝒂𝒂 𝑽𝑽𝑹𝑹 + 𝑌𝑌𝑎𝑎𝑎𝑎 𝑉𝑉𝑎𝑎 (86)

Now we can substitute in power balance equation
𝑽𝑽𝑻𝑻𝑨𝑨 (𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑹𝑹 + 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨 )∗ = 𝑉𝑉𝑎𝑎 (𝒀𝒀𝒂𝒂𝒂𝒂 𝑽𝑽𝑹𝑹 + 𝑌𝑌𝑎𝑎𝑎𝑎 𝑉𝑉𝑎𝑎 )∗ (87)

𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽∗𝑹𝑹 + 𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽∗𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝒀𝒀∗𝒂𝒂𝒂𝒂 𝑽𝑽∗𝑹𝑹 + 𝑉𝑉𝑎𝑎 𝑌𝑌𝑎𝑎𝑎𝑎
∗
𝑉𝑉𝑎𝑎∗
This can only be satisfied if
𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽∗𝑹𝑹 = 𝑉𝑉𝑎𝑎 𝒀𝒀∗𝒂𝒂𝒂𝒂 𝑽𝑽∗𝑹𝑹 (88)
i.e.
𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝒀𝒀∗𝒂𝒂𝒂𝒂 (89)

and
𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨∗ = 𝑉𝑉𝑎𝑎 𝑌𝑌𝑎𝑎𝑎𝑎
∗
𝑉𝑉𝑎𝑎∗ (90)
as
𝑽𝑽𝑨𝑨 𝑉𝑉𝑎𝑎−1 = 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (91)

𝑽𝑽𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (92)

𝑽𝑽𝑻𝑻𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (93)

Substituting in above equation

𝑉𝑉𝑎𝑎 𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀∗𝑨𝑨𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝒀𝒀∗𝒂𝒂𝒂𝒂 (94)

𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀∗𝑨𝑨𝑨𝑨 = 𝒀𝒀∗𝒂𝒂𝒂𝒂 (95)

 𝒀𝒀𝒂𝒂𝒂𝒂 = 𝑽𝑽∗𝑻𝑻
𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀𝑨𝑨𝑨𝑨 (96)

Using the second equation

𝑽𝑽𝑻𝑻𝑨𝑨 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨∗ = 𝑉𝑉𝑎𝑎 𝑌𝑌𝑎𝑎𝑎𝑎
∗
𝑉𝑉𝑎𝑎∗ (97)

Substituting 𝑽𝑽𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕

𝑽𝑽𝑻𝑻𝑨𝑨 = 𝑉𝑉𝑎𝑎 𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (98)

and
𝑽𝑽𝑨𝑨∗ = 𝑽𝑽∗𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝑉𝑉𝑎𝑎∗ (99)

Gives
𝑉𝑉𝑎𝑎 𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽∗𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝑉𝑉𝑎𝑎∗ = 𝑉𝑉𝑎𝑎 𝑌𝑌𝑎𝑎𝑎𝑎
∗
𝑉𝑉𝑎𝑎∗ (100)

results in
∗
𝑌𝑌𝑎𝑎𝑎𝑎 = 𝑽𝑽𝑻𝑻𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀∗𝑨𝑨𝑨𝑨 𝑽𝑽∗𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (101)

∗𝑻𝑻
𝑌𝑌𝑎𝑎𝑎𝑎 = 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 (102)

Now we compile the matrix equation

𝑰𝑰 𝒀𝒀 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑹𝑹 (103)
� 𝑹𝑹 � = � 𝑹𝑹𝑹𝑹 �� �
𝐼𝐼𝑎𝑎 𝒀𝒀𝒂𝒂𝒂𝒂 𝑌𝑌𝑎𝑎𝑎𝑎 𝑉𝑉𝑎𝑎
as

𝑰𝑰 𝒀𝒀𝑹𝑹𝑹𝑹 𝒀𝒀𝑹𝑹𝑹𝑹 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝑽𝑽 (104)

� 𝑹𝑹 � = � ∗𝑻𝑻 ∗𝑻𝑻 � � 𝑹𝑹 �
𝐼𝐼𝑎𝑎 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝒀𝒀𝑨𝑨𝑨𝑨 𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 𝑉𝑉𝑎𝑎

where
𝑽𝑽𝑨𝑨,𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕 = 𝑽𝑽𝑨𝑨 𝑉𝑉𝑎𝑎−1 (105)

If the classical transient model is considered and the voltage magnitude of the system
is considered to be maintained at the transient and the difference of voltage occurs
only in the angle, then the angle of 𝑉𝑉𝑎𝑎 , 𝛿𝛿𝑎𝑎 can be derived as,
 𝑀𝑀𝑖𝑖 𝛿𝛿𝑖𝑖′ (106)
𝛿𝛿𝑎𝑎 = �
𝑀𝑀𝑖𝑖
𝑖𝑖∈[𝐴𝐴]
whereas the magnitude of the voltages are considered as initial value of the voltage
�𝑉𝑉𝑎𝑎,𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 � .

Aggregation of Generator Units

In order to aggregate the network nodes it is essential to aggregate the generators
fonnection to the subset of nodes {𝐴𝐴} which are subjected to aggregation. It is
understood that all the coherent generators oscillates together having a steady angle
deviation. Therefor all such generators can be assumed to be connected on the same
shaft as shown in the figure below.
Pm1 Pm2 Pmn Pma

M1 M2 Mn Ma

Pe2 Pen Pea

Pe1

Then we can aggregate these generators in the mechanical context as

𝑀𝑀𝑎𝑎 = � 𝑀𝑀𝑖𝑖 (110)

𝑖𝑖∈{𝐴𝐴}

and

𝑃𝑃𝑚𝑚𝑚𝑚 = � 𝑃𝑃𝑚𝑚,𝑖𝑖 (111)

𝑖𝑖∈{𝐴𝐴}

𝑃𝑃𝑒𝑒𝑒𝑒 = � 𝑃𝑃𝑒𝑒,𝑖𝑖 (112)

𝑖𝑖∈{𝐴𝐴}
Where 𝑀𝑀𝑖𝑖 is the inertial constant, 𝑃𝑃𝑚𝑚,𝑖𝑖 the mechanical power input and 𝑃𝑃𝑒𝑒,𝑖𝑖 the
mechanical power output of the ith generator. Same parameters of the aggregated
equivalent generator are 𝑀𝑀𝑎𝑎 , 𝑃𝑃𝑚𝑚,𝑎𝑎 , and 𝑃𝑃𝑒𝑒,𝑎𝑎 .
Coherency

It is now clear that the aggregation of nodes can only be done if the transformation
𝑽𝑽𝑨𝑨
ratio remains constant during the transient for the subset of nodes identified for
𝑉𝑉𝑎𝑎
the transient. Therefore, at any time 𝑡𝑡
For 𝑖𝑖 ∈ [𝐴𝐴]
𝑉𝑉𝑖𝑖 (𝑡𝑡) 𝑉𝑉�𝑖𝑖 (107)
= = 𝑉𝑉𝐴𝐴,𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡,𝑖𝑖 = 𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶
𝑉𝑉�𝑎𝑎 (𝑡𝑡) 𝑉𝑉�𝑎𝑎

The circumflex denotes the initial stable state for which the reduced network model
is constructed.
For any node 𝑖𝑖, 𝑗𝑗 ∈ [𝐴𝐴] this condition can further be elaborated as,
�𝚤𝚤 (𝑡𝑡) 𝑉𝑉𝑖𝑖 (𝑡𝑡) �
𝑉𝑉 � 𝑉𝑉�𝑖𝑖 � � (108)
= 𝑒𝑒 𝑗𝑗[𝛿𝛿𝚤𝚤 (𝑡𝑡)−𝛿𝛿𝚥𝚥(𝑡𝑡)] = 𝑒𝑒 𝑗𝑗[𝛿𝛿𝚤𝚤 −𝛿𝛿𝚥𝚥] = 𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶
𝑉𝑉�𝑗𝑗 (𝑡𝑡) 𝑉𝑉𝑗𝑗 (𝑡𝑡) 𝑉𝑉�𝑗𝑗
Any two nodes which satisfies the above condition can be considered as coherent.
Usually load nodes are not coherent. Only the nodes which are far away from the
fault remain with constant voltage magnitude and angle. However, it is not
phenomenal to find generators which swing together. Those generator nodes qualify
for the above nodal aggregation and therefore can be represented by a single
generator connected to the aggregated node.
In the classical generator model, the transient voltage of the generator bus is
represented by a constant magnitude voltage with the power flow angle as 𝐸𝐸𝑖𝑖 ∠𝛿𝛿𝑖𝑖 . If
the pre-fault volage of the generator bus is denoted with the circumflex as 𝐸𝐸�𝑖𝑖 ∠𝛿𝛿̂𝑖𝑖
Where for all 𝑖𝑖, 𝑗𝑗 ∈ [𝐴𝐴]

𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑗𝑗 = 𝛿𝛿̂𝑖𝑖 − 𝛿𝛿̂𝑗𝑗 (109)

This can be graphically represented for three generators 𝑖𝑖, 𝑗𝑗, 𝑘𝑘 as shown in figure
below. Generators 𝑖𝑖 and 𝑗𝑗 are coherent and generator 𝑘𝑘 is not coherent with them.

δ
δi

δj

δk

Variation of the rotor angle of tow coherent one generator

which is not coherent with the two.

Coherency Recognition

However, there lies a greater problem as to how the coherency be identified without
running the complete model to find the generators oscillate in phase. If the detailed
dynamic simulation can be run, then the requirement of model reduction through
aggregating the coherent generators is not required. There are generalised methods
available in ascertaining coherency without going through the complete dynamic
analysis. The detailed research works carried out in this area is known as ‘coherency
recognition’ which does not have time for a discussion in this course. Even though
the detailed analysis of the methodology of coherency recognition is quite deep, it
may be possible to analyse the coherency when one of the nodes 𝑘𝑘 in retained node
subset,{𝑅𝑅} is disturbed.
All the load nodes of the retained node subset {𝑅𝑅} are eliminated through modelling
the loads as impedances, and all unconnected nodes are eliminated by above
discussed reduction methods, to create the generator node subset {𝐵𝐵}.
 {R} {R}

{B} {B}

Network
Reduced for
Network Coherency
Check

A disturbance of the voltage of generator node 𝑘𝑘 ∈ {𝐵𝐵} influences the generators in

the Retained Subsystem through the equivalent transfer admittances of the network.
1
2

Network for Coherency

i
{B} Hik

Check
Mi
j
Mj Hjk
k

If the branch impedances are considered totally reactive (Transfer Conductance,

𝐺𝐺𝑖𝑖𝑖𝑖 = 0). Then the active power produced by the generators can be expressed as,

𝑃𝑃𝑖𝑖 = �𝐸𝐸𝑖𝑖2 �𝐺𝐺𝑖𝑖𝑖𝑖 + � 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 sin 𝛿𝛿𝑖𝑖𝑖𝑖 (113)

𝑘𝑘∈{𝐵𝐵}

Where
𝛿𝛿𝑖𝑖𝑖𝑖 = 𝛿𝛿𝑖𝑖 − 𝛿𝛿𝑘𝑘 (114)
Let’s say that the voltage disturbance occurs at a border node 𝑘𝑘 from the initial angle
of 𝛿𝛿𝑘𝑘0 to a value 𝛿𝛿𝑘𝑘 = 𝛿𝛿𝑘𝑘0 + Δ𝛿𝛿𝑘𝑘 and the voltage at other nodes are constant, this
change of angle will result in a change in the generation at node 𝑖𝑖 as
Δ𝑃𝑃𝑖𝑖 = 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 [sin(𝛿𝛿𝑖𝑖𝑖𝑖 + Δ𝛿𝛿𝑘𝑘 ) − sin 𝛿𝛿𝑖𝑖𝑖𝑖 ] (115)
Δ𝑃𝑃𝑖𝑖 = 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 [sin 𝛿𝛿𝑖𝑖𝑖𝑖 cos Δ𝛿𝛿𝑘𝑘 + cos 𝛿𝛿𝑖𝑖𝑖𝑖 sin Δ𝛿𝛿𝑘𝑘 − sin 𝛿𝛿𝑖𝑖𝑖𝑖 ] (116)
As 𝛿𝛿𝑘𝑘 is small
Δ𝑃𝑃𝑖𝑖 ≃ 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 [sin 𝛿𝛿𝑖𝑖𝑖𝑖 ∙ 1 + cos 𝛿𝛿𝑖𝑖𝑖𝑖 ∙ Δ𝛿𝛿𝑘𝑘 ] − sin 𝛿𝛿𝑖𝑖𝑖𝑖 ] (117)
Δ𝑃𝑃𝑖𝑖 ≃ 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 cos 𝛿𝛿𝑖𝑖𝑖𝑖 Δ𝛿𝛿𝑘𝑘 (118)
We can call 𝐸𝐸𝑖𝑖 𝐸𝐸𝑘𝑘 𝐵𝐵𝑖𝑖𝑖𝑖 cos 𝛿𝛿𝑖𝑖𝑖𝑖 = 𝐻𝐻𝑖𝑖𝑖𝑖 which is the synchronizing power constant
between bust 𝑖𝑖 and 𝑘𝑘.
This power flow will cause a rotor acceleration of
𝑑𝑑𝜔𝜔𝑖𝑖 (119)
Δ𝑃𝑃𝑖𝑖 (Δδ𝑘𝑘 ) = 𝑀𝑀𝑖𝑖
𝑑𝑑𝑑𝑑

𝑑𝑑𝜔𝜔𝑖𝑖 Δ𝑃𝑃𝑖𝑖 (Δδ𝑘𝑘 ) (120)

=
𝑑𝑑𝑑𝑑 𝑀𝑀𝑖𝑖
𝑑𝑑𝜔𝜔𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 (121)
= Δ𝛿𝛿
𝑑𝑑𝑑𝑑 𝑀𝑀𝑖𝑖 𝑘𝑘
We can write the same equation for another generator 𝑗𝑗 as
𝑑𝑑𝜔𝜔𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 (122)
= Δ𝛿𝛿𝑘𝑘
𝑑𝑑𝑑𝑑 𝑀𝑀𝑗𝑗
Now if the two generators are electromechanically coherent then the rotor
acceleration of the two generators essentially needs to be same.
Therefore ideal coherency means
𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑗𝑗𝑗𝑗 (123)
=
𝑀𝑀𝑖𝑖 𝑀𝑀𝑗𝑗

𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑗𝑗𝑗𝑗 (124)

− = 𝜌𝜌
𝑀𝑀𝑖𝑖 𝑀𝑀𝑗𝑗

Where 𝜌𝜌 ≈ 0
Equation (61) states that
 0 0 0 0 ⋮ 1 0 ⋯ 0
⎡ 0 0 0 0 ⋮ 0 1 ⋯ 0 ⎤
∆𝛿𝛿̇ ⎢ ⎥ ∆𝛿𝛿1
⎡ 1⎤ ⎢ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮
⎥ ⎡ ∆𝛿𝛿2 ⎤
̇
⎢ ∆𝛿𝛿2 ⎥ ⎢ 0 0 0 0 ⋮ 0 0 ⋯ 1 ⎥⎢ ⋯ ⎥
⎢ ⋯ ⎥ ⎢ − − − − ⋮ − − − − ⎥⎢ ⎥
⎢ ∆𝛿𝛿𝑛𝑛̇ ⎥ ⎢ 𝐻𝐻11 𝐻𝐻12 𝐻𝐻1𝑛𝑛 𝐷𝐷1 ⎥ ⎢ ∆𝛿𝛿𝑛𝑛 ⎥
⎢ ⋯ ⎥= − − ⋯ − ⋮ − 0 ⋯ 0
⎢ ⎥ ⎢ 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 ⎥⎢ ⋯ ⎥
⎢ 1 ⎥ ⎢ 𝐻𝐻21
∆𝜔𝜔̇ 𝐻𝐻22 𝐻𝐻2𝑛𝑛 𝐷𝐷2 ⎥ ⎢∆𝜔𝜔1 ⎥
⎢∆𝜔𝜔̇ 2 ⎥ ⎢− 𝑀𝑀2 − ⋯ − ⋮ 0 − ⋯ 0 ⎥ ⎢∆𝜔𝜔2 ⎥
𝑀𝑀2 𝑀𝑀2 𝑀𝑀2
⎢ ⋮ ⎥ ⎢ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⎥
⎥⎢ ⋮ ⎥
⎣∆𝜔𝜔̇ 𝑛𝑛 ⎦ ⎢ 𝐻𝐻𝑛𝑛1 𝐻𝐻𝑛𝑛2 𝐻𝐻𝑛𝑛𝑛𝑛 𝐷𝐷𝑛𝑛 ⎥ ⎣∆𝜔𝜔𝑛𝑛 ⎦
⎢− − ⋯ − ⋮ 0 0 ⋯ −
⎣ 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 ⎦

From which we can extract the generator dynamic equations as,

⋱ 1 2 ⋯ 𝑖𝑖 ⋯ 𝑗𝑗 ⋯ 𝑘𝑘 ⋯ 𝑛𝑛 (125)
⎡ 𝐻𝐻1𝑗𝑗 ⎤
𝐻𝐻11 𝐻𝐻12 𝐻𝐻1𝑖𝑖 𝐻𝐻1𝑘𝑘 𝐻𝐻1𝑛𝑛
⎢1 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎥
⎢ 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 ⎥
⎢ 𝐻𝐻21 𝐻𝐻22 𝐻𝐻2𝑖𝑖 𝐻𝐻2𝑗𝑗 𝐻𝐻2𝑘𝑘 𝐻𝐻2𝑛𝑛 ⎥
∆𝜔𝜔̇ 1 2 − − ⋯ − ⋯ − ⋯ − ⋯ − ∆𝛿𝛿1
⎡∆𝜔𝜔̇ ⎤ ⎢ 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 ⎥ ⎡
∆𝛿𝛿
⎤
⎢ 2⎥ ⎢ ⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ 2⎥
⋮ ⋮
⎢ ⎥ ⎢ 𝐻𝐻𝑖𝑖1 𝐻𝐻𝑖𝑖2 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 ⎥ ⎢ ⎥
∆𝜔𝜔̇ 𝑖𝑖 ⎥ ⎢ Δ𝛿𝛿𝑖𝑖 ⎥
⎢ ⎥ ⎢ 𝑖𝑖 −
𝑀𝑀𝑖𝑖
−
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖 ⎥ ⋮
⋮ ⎢
⎢ ⎥= ⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ Δ𝛿𝛿 ⎥
∆𝜔𝜔̇ ⎢
⎢ 𝑗𝑗 ⎥ 𝐻𝐻𝑗𝑗1 𝐻𝐻𝑗𝑗2 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 ⎥ ⎢ 𝑗𝑗 ⎥
⎢
⎢ ⋮ ⎥ ⎢ 𝑗𝑗 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎢ ⋮ ⎥
𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 ⎥ Δ𝛿𝛿
⎢∆𝜔𝜔̇ 𝑘𝑘 ⎥ ⎢ ⎢ 𝑘𝑘 ⎥
⎢ ⋮ ⎥ ⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣∆𝜔𝜔̇ 𝑛𝑛 ⎦ ⎢𝑘𝑘 𝐻𝐻𝑘𝑘1 𝐻𝐻𝑘𝑘2 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 ⎥ ⎣∆𝛿𝛿 ⎦
− − ⋯ − ⋯ − ⋯ − ⋯ − ⎥ 𝑛𝑛
⎢ 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 ⎥
⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥
⎢𝑛𝑛 𝐻𝐻𝑛𝑛1 𝐻𝐻𝑛𝑛2 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 ⎥
− − ⋯ − ⋯ − ⋯ − ⋯ −
⎣ 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 ⎦

Now equation means that there are certain elements in the column 𝑘𝑘 of the dynamic
coefficient matrix in equation (126) are similar or nearly similar. That means for a
disturbance occurring in the generator 𝑘𝑘 results in a coherent disturbance in such
cluster of generators.
 ⋱ 1 2 ⋯ 𝑖𝑖 ⋯ 𝑗𝑗 ⋯ 𝑘𝑘 ⋯ 𝑛𝑛 (126)
⎡ 𝐻𝐻1𝑗𝑗 ⎤
𝐻𝐻11 𝐻𝐻12 𝐻𝐻1𝑖𝑖 𝐻𝐻1𝑘𝑘 𝐻𝐻1𝑛𝑛
⎢1 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎥
⎢ 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 ⎥
⎢ 𝐻𝐻21 𝐻𝐻22 𝐻𝐻2𝑖𝑖 𝐻𝐻2𝑗𝑗 𝐻𝐻2𝑘𝑘 𝐻𝐻2𝑛𝑛 ⎥
∆𝜔𝜔̇ 1 2 − − ⋯ − ⋯ − ⋯ − ⋯ − ∆𝛿𝛿1
⎡∆𝜔𝜔̇ ⎤ ⎢ 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 ⎥ ⎡
∆𝛿𝛿2
⎤
⎢ 2
⎥ ⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ ⎥
⋮ ⋮
⎢ ⎥ ⎢ 𝐻𝐻𝑖𝑖1 𝐻𝐻𝑖𝑖2 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 ⎥ ⎢ ⎥
∆𝜔𝜔̇ 𝑖𝑖 ⎥ Δ𝛿𝛿𝑖𝑖
⎢ ⎥ ⎢ 𝑖𝑖 −
𝑀𝑀𝑖𝑖
−
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖 ⎥ ⎢ ⋮ ⎥
⋮ ⎢
⎢ ⎥= ⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ Δ𝛿𝛿 ⎥
⎢
⎢ ∆𝜔𝜔̇ 𝑗𝑗 ⎥ 𝐻𝐻𝑗𝑗1 𝐻𝐻𝑗𝑗2 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 ⎥ ⎢ 𝑗𝑗 ⎥
⎢
⎢ ⋮ ⎥ ⎢ 𝑗𝑗 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎢ ⋮ ⎥
𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 ⎥ Δ𝛿𝛿
⎢∆𝜔𝜔̇ 𝑘𝑘 ⎥ ⎢ ⎢ 𝑘𝑘 ⎥
⎢ ⋮ ⎥ ⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣∆𝜔𝜔̇ 𝑛𝑛 ⎦ ⎢𝑘𝑘 𝐻𝐻𝑘𝑘1 𝐻𝐻𝑘𝑘2 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 ⎥ ⎣∆𝛿𝛿 ⎦
− − ⋯ − ⋯ − ⋯ − ⋯ − ⎥ 𝑛𝑛
⎢ 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 ⎥
⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥
⎢𝑛𝑛 𝐻𝐻𝑛𝑛1 𝐻𝐻𝑛𝑛2 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 ⎥
− − ⋯ − ⋯ − ⋯ − ⋯ −
⎣ 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 ⎦

Now if there is a cluster of generators which coherently oscillates for a disturbance

in any of the nodes, then as shown in equation (127) nearly equal values shall be in
all the columns of the dynamic coefficient matrix corresponding to the particular set
of generators.
⋱ 1 2 ⋯ 𝑖𝑖 ⋯ 𝑗𝑗 ⋯ 𝑘𝑘 ⋯ 𝑛𝑛 (127)
⎡ 𝐻𝐻1𝑗𝑗 ⎤
𝐻𝐻11 𝐻𝐻12 𝐻𝐻1𝑖𝑖 𝐻𝐻1𝑘𝑘 𝐻𝐻1𝑛𝑛
⎢1 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎥
⎢ 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 𝑀𝑀1 ⎥
⎢ 𝐻𝐻21 𝐻𝐻22 𝐻𝐻2𝑖𝑖 𝐻𝐻2𝑗𝑗 𝐻𝐻2𝑘𝑘 𝐻𝐻2𝑛𝑛 ⎥
∆𝜔𝜔̇ 1 2 − − ⋯ − ⋯ − ⋯ − ⋯ − ∆𝛿𝛿1
⎡∆𝜔𝜔̇ ⎤ ⎢ 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 𝑀𝑀2 ⎥ ⎡
∆𝛿𝛿
⎤
⎢ 2⎥ ⎢ ⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ 2⎥
⋮ ⋮
⎢ ⎥ ⎢ 𝐻𝐻𝑖𝑖1 𝐻𝐻𝑖𝑖2 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 𝐻𝐻𝑖𝑖𝑖𝑖 ⎥ ⎢ ⎥
∆𝜔𝜔̇ 𝑖𝑖 ⎥ ⎢ Δ𝛿𝛿𝑖𝑖 ⎥
⎢ ⎥ ⎢ 𝑖𝑖 −
𝑀𝑀𝑖𝑖
−
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖
⋯ −
𝑀𝑀𝑖𝑖 ⎥ ⋮
⋮ ⎢
⎢ ⎥= ⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥ ⎢ Δ𝛿𝛿 ⎥
∆𝜔𝜔̇ ⎢
⎢ 𝑗𝑗 ⎥ 𝐻𝐻𝑗𝑗1 𝐻𝐻𝑗𝑗2 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 𝐻𝐻𝑗𝑗𝑗𝑗 ⎥ ⎢ 𝑗𝑗 ⎥
⎢
⎢ ⋮ ⎥ ⎢ 𝑗𝑗 − − ⋯ − ⋯ − ⋯ − ⋯ − ⎢ ⋮ ⎥
𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 𝑀𝑀𝑗𝑗 ⎥ Δ𝛿𝛿
⎢∆𝜔𝜔̇ 𝑘𝑘 ⎥ ⎢ ⎢ 𝑘𝑘 ⎥
⎢ ⋮ ⎥ ⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⎥
⎣∆𝜔𝜔̇ 𝑛𝑛 ⎦ ⎢𝑘𝑘 𝐻𝐻𝑘𝑘1 𝐻𝐻𝑘𝑘2 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 𝐻𝐻𝑘𝑘𝑘𝑘 ⎥ ⎣∆𝛿𝛿 ⎦
− − ⋯ − ⋯ − ⋯ − ⋯ − ⎥ 𝑛𝑛
⎢ 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 𝑀𝑀𝑘𝑘 ⎥
⎢⋮ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⋱ ⋮ ⎥
⎢𝑛𝑛 𝐻𝐻𝑛𝑛1 𝐻𝐻𝑛𝑛2 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 𝐻𝐻𝑛𝑛𝑛𝑛 ⎥
− − ⋯ − ⋯ − ⋯ − ⋯ −
⎣ 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 𝑀𝑀𝑛𝑛 ⎦
Example E.2. Single Machine Model

∞
P+jQ
X
G

M, D
V1 V2

Single machine connected to the infinite bus-bar

Equation (61) can be written for a single machine model shown above as,
0 1
∆𝛿𝛿 ̇ 𝐷𝐷 � � ∆𝛿𝛿 �
� � = � 𝐻𝐻
∆𝜔𝜔̇ − − ∆𝜔𝜔
𝑀𝑀 𝑀𝑀
Where
𝑉𝑉1 𝑉𝑉2 𝑉𝑉1 𝑉𝑉2
𝐻𝐻 = cos(𝛿𝛿1 − 𝛿𝛿2 ) = cos 𝛿𝛿
𝑋𝑋 𝑋𝑋
The Eigen values of this equation is,

−𝐷𝐷 𝐷𝐷2 𝐻𝐻
𝜆𝜆1 , 𝜆𝜆2 = ±� 2−
2𝑀𝑀 4𝑀𝑀 𝑀𝑀

And thee response can be given as

𝐷𝐷 𝐷𝐷2 𝐻𝐻 𝐷𝐷2 𝐻𝐻
𝑒𝑒 𝜆𝜆1 𝑡𝑡 = 𝑒𝑒 −2𝑀𝑀𝑡𝑡 �cos �� 2 − 𝑡𝑡� + 𝑗𝑗 sin �� 2 − 𝑡𝑡��
4𝑀𝑀 𝑀𝑀 4𝑀𝑀 𝑀𝑀

𝐷𝐷 𝐷𝐷2 𝐻𝐻 𝐷𝐷2 𝐻𝐻
𝑒𝑒 𝜆𝜆2 𝑡𝑡 = 𝑒𝑒 −2𝑀𝑀𝑡𝑡 �cos �� 2 − 𝑡𝑡� − 𝑗𝑗 sin �� 2 − 𝑡𝑡��
4𝑀𝑀 𝑀𝑀 4𝑀𝑀 𝑀𝑀

There is a point of transition of the roots from imaginary to real when 𝛿𝛿 increases. This transition point can
be calculated as,

𝐷𝐷2
𝐻𝐻𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 =
4𝑀𝑀
And the power delivery angle at the transition can be calculated as,

𝑉𝑉1 𝑉𝑉2 𝐷𝐷2

𝐻𝐻𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 = cos 𝛿𝛿𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 =
𝑋𝑋 4𝑀𝑀
Gives

𝐷𝐷2 𝑋𝑋
𝛿𝛿𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 = cos −1 � �
4𝑀𝑀 ∙ 𝑉𝑉1 𝑉𝑉2

This can be graphically represented as,

P, KE

KE

δtransition δ

Transition of behaviour with the increase of power delivery angle

Example E.3. Small Signal Model of IEEE 9 bus system

Let’s consider the IEEE 9 bus model we used in the transient analysis. The network
is given as below in figure 1 and the pre-fault load flow data Is given in figure 2.

Figure 1-IEEE 9 Bus system (All impedances in 100MVA base)

 Figure 2-Pre-fault load flow data of IEEE (-bus system

The Generator parameters are given as;

And the internal voltages of the generators are given as;

𝐸𝐸1 = 1.0566∡2.2717°

𝐸𝐸2 = 1.0502∡19.7315°

𝐸𝐸3 = 1.0170∡13.1752°
The Ybus Matrix once constructed as given in MATLAB program in Appendix 2, can
be given as
𝑌𝑌𝑏𝑏𝑏𝑏𝑏𝑏
− 8.45i 0 0 8.45i 0 0 0 0 0
⎡ 0 −5.5i 0 0 0 0 5.49i 0 0 ⎤
⎢ ⎥
0 0 − 4.2i 0 0 0 0 0 4,2i
⎢ ⎥
⎢ 8.45i 0 0 3.3 − 30.4i −1.4 + 11.6i −1.9 + 10.5i 0 0 0 ⎥
=⎢ 0 0 0 1.4 + 11.6i 3.9 − 17.7i 0 −1.2 + 6.0i 0 0 ⎥
⎢ 0 0 0 −1.9 + 10i 0 4.1 − 16.0i 0 0 1.3 + 5.6i ⎥
⎢ 0 5.5i 0 0 −1.2 + 6.0i 0 2.8 − 24.9i −1.6 + 13.7i 0 ⎥
⎢ 0 0 0 0 0 0 −1.6 + 13.7i 3.7 − 23.7i −1.2 + 9.8i⎥
⎣ 0 0 4.2i 0 0 −1.3 + 5.6i 0 −1.2 + 9.8i 2.4 − 19.3i ⎦

The reduced Ybus Matrix and the state matrix of the system is as given below,
0.90 − 2.93𝑖𝑖 0.31 + 1.54𝑖𝑖 0.23 + 1.25𝑖𝑖
𝑌𝑌𝑏𝑏𝑏𝑏𝑏𝑏,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 = �0.31 + 1.54𝑖𝑖 0.43 − 2.71𝑖𝑖 0.22 + 1.10𝑖𝑖 �
0.23 + 1.25𝑖𝑖 0.22 + 1.10𝑖𝑖 0.28 − 2.36𝑖𝑖

−0.78 −937.59 0 511.83 0 425.76

⎡ 1 0 0 0 0 0 ⎤
⎢ ⎥
0 1836.44 −2.81 −2980.46 0 1144.02 ⎥
𝐴𝐴 = ⎢
⎢ 0 0 0 1 0 0 ⎥
⎢ 0 3251.83 0 1145.26 −5.98 −5686.94⎥
⎣ 0 0 0 0 1 0 ⎦

The Eigen matrix of the system is given as;

−0.8 + 10.8𝑖𝑖 0 0 0 0 0
⎡ 0 −0.8 − 10.8𝑖𝑖 0 0 0 0 ⎤
⎢ 0 0 −1.2 + 57.1𝑖𝑖 0 0 0
⎥
⎢ ⎥
⎢ 0 0 0 −1.2 − 57.1𝑖𝑖 0 0 ⎥
⎢ 0 0 0 0 −2.8 + 78.8𝑖𝑖 0 ⎥
⎣ 0 0 0 0 0 −2.8 − 78.8𝑖𝑖⎦

The mode shape oscillations corresponding to these three pairs of conjugate motes
is as shown in figure 3.
Figure 3-Mode Shape Oscillations
The Right Eigen Matrix of the System is as shown in table below and the mode
shapes of oscillations are as shown in figure below.
State Mode 1 (𝑓𝑓 = 1.72, Mode 2 (𝑓𝑓 = 9.09, 𝜉𝜉 = Mode 1 (𝑓𝑓 = 12.54,
Variable 𝜉𝜉 = −0.75) −1.23) 𝜉𝜉 = −2.80)
Real Imag. Real Imag. Real Imag.
𝜔𝜔1 0.627 0.000 -0.236 0.010 -0.047 0.001
𝛿𝛿1 -0.004 -0.058 0.000 0.004 0.000 0.001
𝜔𝜔2 0.598 -0.005 0.962 0.000 -0.309 0.020
𝛿𝛿2 -0.004 -0.055 0.000 -0.017 0.000 0.004
𝜔𝜔3 0.490 -0.005 0.138 0.001 0.950 0.000
𝛿𝛿3 -0.004 -0.045 0.000 -0.002 0.000 -0.012

Mode Shapes of the Machine δ

It is seen that in mode 1 oscillation all the machines oscillates in-phase. In the mode
2 oscillation, G1 and G2 machines oscillates 180 degrees out of phase around G3.
In mode 3 oscillation, G2 and G3 machines oscillation 180 degrees out of phase
around G1.
Appendix 1 – Axis transformation of example E.1

Let’s plot 𝑧𝑧1 and 𝑧𝑧2 axis of example E.1, in x-y plane.
In this example the right Eigen matrix of the system ∆𝑥𝑥̇ = 𝐴𝐴∆𝑥𝑥 is;
−0.97 −0.32
𝜙𝜙 = � �
−0.24 −0.95
Therefore, the system can be transformed in to the decoupled domain by;
𝒙𝒙 = 𝝓𝝓𝝓𝝓
𝑥𝑥
If we consider the coupled domain axis are x,y, then 𝒙𝒙 = �𝑦𝑦�, we can plot the
location of z1 and z2 axis by plotting a cluster of points as below.
We can plot the 𝑧𝑧1 axis by transforming the set of 𝑧𝑧 vectors
𝑧𝑧1 −5 −4 0 5
�𝑧𝑧 � = � � , � � , ⋯ � � ⋯ � �
2 0 0 0 0
to x-y domain by 𝒙𝒙 = 𝝓𝝓𝝓𝝓
𝑥𝑥 −0.97 −0.32 𝑧𝑧1
�𝑦𝑦� = � �� �
−0.24 −0.95 𝑧𝑧2
For example
𝑥𝑥 −0.97 −0.32 −5 (−0.97) × −5 + (−0.32) × 0 4.85
�𝑦𝑦� = � �� � = � �=� �
−0.24 −0.95 0 (−0.24) − 5 + (−0.95) × 0 1.20
We can write this as,
𝑥𝑥 −0.97 −0.32 𝑧𝑧1
�𝑦𝑦� = � �� �
−0.24 −0.95 0
This represents an axis such as
−0.24
𝑥𝑥 = −0.97𝑧𝑧1 and 𝑦𝑦 = −0.24𝑧𝑧1 which gives 𝑦𝑦 = 𝑥𝑥 = 0.2474𝑥𝑥
−0.97

Similarly, we can transform 𝑧𝑧2 axis to x-y plane by transforming

𝑧𝑧1 0 0 0 0
�𝑧𝑧 � = � � , � � , ⋯ � � ⋯ � �
2 −5 −4 0 5
Such as
 𝑥𝑥 −0.97 −0.32 0
�𝑦𝑦� = � �� �
−0.24 −0.95 𝑧𝑧2
−0.95
𝑥𝑥 = −0.32𝑧𝑧2 and 𝑦𝑦 = −0.95𝑧𝑧1 which gives 𝑦𝑦 = 𝑥𝑥 = 2.0688𝑥𝑥
−0.32

The calculation is shown below with the two z-axis resulted in by such calculation.
𝑥𝑥 𝑧𝑧1
Therefore, it shall be noted that �𝑦𝑦� axis system can be transformed to �𝑧𝑧 � axis
2
system and visa-versa by plotting in the co-ordinate system shown in the figure.
phi
-0.9701 -0.3162
-0.2425 -0.9487
x z
4.8505 1.2125 -5 0
3.8804 0.97 -4 0
2.9103 0.7275 -3 0
1.9402 0.485 -2 0
0.9701 0.2425 -1 0
0 0 0 0
-0.9701 -0.2425 1 0
-1.9402 -0.485 2 0
-2.9103 -0.7275 3 0
-3.8804 -0.97 4 0
-4.8505 -1.2125 5 0
1.581 4.7435 0 -5
1.2648 3.7948 0 -4
0.9486 2.8461 0 -3
0.6324 1.8974 0 -2
0.3162 0.9487 0 -1
0 0 0 0
-0.3162 -0.9487 0 1
-0.6324 -1.8974 0 2
-0.9486 -2.8461 0 3
-1.2648 -3.7948 0 4
-1.581 -4.7435 0 5
 -5

-4

-3

-2

-1 -5
-4
-3
-1 -2
0
1
2
3 1
5 4
2

5
Appendix 2 – MATLAB Program to calculate the modes of the IEEE 9 bus
system
% Small signal model of the IEEE 9 bus system

clear;

alpha=[0;9.3;4.7];
v=[1.040;1.025;1.025].*(cos(alpha/180*pi)+j*sin(alpha/180*pi));

%Generator impedences
xdd1=j*0.0608;
xdd2=j*0.1198;
xdd3=j*0.1813;
xdd=[xdd1;xdd2;xdd3];

P_Gen=[.716;1.63;0.85];
Q_Gen=[.27;.062;-.109];
E_Gen=[v+Q_Gen.*xdd./v+j*P_Gen.*xdd./v];
theta=angle(E_Gen);
v=abs(E_Gen).*(cos(theta)+j*sin(theta));

%transformer impedances
trf14=j*0.0576;
trf27=j*0.0625;
trf39=j*0.0586;

z14=xdd1+trf14;
z27=xdd2+trf27;
z39=xdd3+trf39;

y14=1/z14;
y27=1/z27;
y39=1/z39;

z45=.0100+j*0.0850;
z46=.0170+j*0.0920;
z57=0.0320+j*0.1610;
z69=0.0390+j*0.1700;
z78=.0085+j*0.0720;
z89=0.0119+j*0.1008;

y45=1/z45;
y46=1/z46;
y57=1/z57;
y69=1/z69;
y78=1/z78;
y89=1/z89;

%load calculation
VA=0.996*exp(j*deg2rad(-4));
PA=125/100;
QA=50.0/100;
YA=PA/VA^2-j*QA/VA^2;
VB=1.013*exp(j*deg2rad(-3.7));
PB=90/100;
QB=30/100;
YB=PB/VB^2-j*QB/VB^2;

VC=1.016*exp(j*deg2rad(0.7));
PC=100/100;
QC=35/100;
YC=PC/VC^2-j*QC/VC^2;

%Shut admittances of lines

B45=j*0.088;
B57=j*0.153;
B78=j*0.0745;
B89=j*0.1045;
B69=j*0.179;
B46=j*0.079;

y40=B45+B46;
y70=B57+B78;
y90=B89+B69;

y50=YA+B57+B45;
y60=YB+B46+B69;
y80=YC+B89+B78;

Ybus=[y14 0 0 -y14 0 0 0 0 0
0 y27 0 0 0 0 -y27 0 0
0 0 y39 0 0 0 0 0 -y39
-y14 0 0 y14+y45+y46+y40 -y45 -y46 0 0 0
0 0 0 -y45 y45+y57+y50 0 -y57 0 0
0 0 0 -y46 0 y46+y69+y60 0 0 -y69
0 -y27 0 0 -y57 0 y27+y57+y78+y70 -y78 0
0 0 0 0 0 0 -y78 y78+y89+y80 -y89
0 0 -y39 0 0 -y69 0 -y89 y39+y69+y89+y90];

ynn=Ybus(1:3,1:3);
ynr=Ybus(1:3,4:9);
yrn=Ybus(4:9,1:3);
yrr=Ybus(4:9,4:9);
Ybusnew=ynn-ynr*inv(yrr)*yrn;

G=real(Ybusnew);
B=imag(Ybusnew);
vv=abs(v);
H=[0 0 0
0 0 0
0 0 0];

theta=theta/180*pi;

for i=1:3
tot=0;
for k=1:3
if i==k
H(i,k)=0;
 else
H(i,k)=vv(i)*vv(k)*(G(i,k)*(sin(theta(i)-theta(k)))-
B(i,k)*(cos(theta(i)-theta(k))));
tot=tot+H(i,k);
end
end
H(i,i)=-tot;
end

M=[23;6.4;3.01]/(180*50);
D=[.002 .002 .002];

A=[-D(1)/M(1) -H(1,1)/M(1) 0 -H(1,2)/M(1) 0 -H(1,3)/M(1)

1 0 0 0 0 0
0 -H(2,1)/M(2) -D(2)/M(2) -H(2,2)/M(2) 0 -H(2,3)/M(2)
0 0 1 0 0 0
0 -H(3,1)/M(3) 0 -H(3,2)/M(2) -D(3)/M(3) -H(3,3)/M(3)
0 0 0 0 1 0]

% V is the right eigen matrix where A*V=V*lambda

% W is the left eigen matrix where W'*A=lambda*W'

[V,lambda,W]=eig(A);

%participation factor calculation

part=zeros(6);
W=W';

for i=1:6
tott=0;
for k=1:6
part(k,i)=abs(V(k,i))*abs(W(i,k));
tott=tott+part(k,i);
end
for kk=1:6
part(kk,i)=part(kk,i)/tott;
end
end
Appendix 3 – Example of three mass system

Consider three masses connected by three spring-damper systems. This system is having three degrees of
freedom. We can describe dynamics of this system by three equations for three masses as;

𝑚𝑚1 𝑥𝑥̈1 = −𝑘𝑘1 (𝑥𝑥1 − 𝑥𝑥2 ) − 𝑐𝑐1 (𝑣𝑣1 − 𝑣𝑣2 )

𝑚𝑚2 𝑥𝑥̈ 2 = 𝑘𝑘1 (𝑥𝑥1 − 𝑥𝑥2 ) − 𝑘𝑘2 (𝑥𝑥2 − 𝑥𝑥3 ) + 𝑐𝑐1 (𝑣𝑣1 − 𝑣𝑣2 ) − 𝑐𝑐2 (𝑣𝑣2 − 𝑣𝑣3 )

= 𝑘𝑘1 𝑥𝑥1 − (𝑘𝑘1 + 𝑘𝑘2 )𝑥𝑥2 + 𝑘𝑘2 𝑥𝑥3 + 𝑐𝑐1 𝑣𝑣1 − (𝑐𝑐1 + 𝑐𝑐2 )𝑣𝑣2 + 𝑐𝑐2 𝑣𝑣3

𝑚𝑚3 𝑥𝑥̈ 3 = 𝑘𝑘2 𝑥𝑥2 − 𝑘𝑘2 𝑥𝑥3 + 𝑐𝑐2 𝑣𝑣2 − 𝑐𝑐2 𝑣𝑣3

The state space interpretation of this system can be written as;

𝑥𝑥̇1 = 𝑣𝑣1
𝑘𝑘1 𝑘𝑘1 𝑐𝑐1 𝑐𝑐1
𝑣𝑣̇ 1 = − 𝑥𝑥1 + 𝑥𝑥2 − 𝑣𝑣1 + 𝑣𝑣
𝑚𝑚1 𝑚𝑚1 𝑚𝑚1 𝑚𝑚1 2

𝑥𝑥̇ 2 = 𝑣𝑣2

𝑘𝑘1 (𝑘𝑘1 + 𝑘𝑘2 ) 𝑘𝑘2 𝑐𝑐1 (𝑐𝑐1 + 𝑐𝑐2 ) 𝑐𝑐2

𝑣𝑣̇ 2 = 𝑥𝑥1 − 𝑥𝑥2 + 𝑥𝑥3 + 𝑣𝑣1 − 𝑣𝑣2 + 𝑣𝑣
𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 3

𝑥𝑥̇ 3 = 𝑣𝑣3
𝑘𝑘2 𝑘𝑘2 𝑐𝑐2 𝑐𝑐2
𝑣𝑣̇ 3 = 𝑥𝑥2 − 𝑥𝑥3 + 𝑣𝑣2 − 𝑣𝑣
𝑚𝑚3 𝑚𝑚3 𝑚𝑚3 𝑚𝑚3 3

And the state space representation of this system in free motion is given by,
 0 1 0 0 0 0
⎡ 𝑘𝑘1 𝑐𝑐1 𝑘𝑘1 𝑐𝑐1 ⎤
𝑥𝑥̇1 ⎢− − 0 0 ⎥ 𝑥𝑥1
⎡𝑣𝑣̇ ⎤ ⎢ 𝑚𝑚1 𝑚𝑚1 𝑚𝑚1 𝑚𝑚1 ⎥ ⎡𝑣𝑣 ⎤
⎢ 1⎥ ⎢ 0 0 0 1 0 0 ⎥ ⎢ 1⎥
𝑥𝑥
⎢𝑥𝑥̇ 2 ⎥ = ⎢ 𝑘𝑘1 𝑐𝑐1 (𝑘𝑘1 + 𝑘𝑘2 ) (𝑐𝑐1 + 𝑐𝑐2 ) 𝑘𝑘2 𝑐𝑐2 ⎥ ⎢ 2 ⎥
⎢𝑣𝑣̇ 2 ⎥ ⎢ 𝑚𝑚 − − 𝑣𝑣2
𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 𝑚𝑚2 ⎥ ⎢𝑥𝑥 ⎥
⎢𝑥𝑥̇ 3 ⎥ ⎢ 2 ⎥ ⎢ 3⎥
0 0 0 0 0 1 ⎥ ⎣𝑣𝑣 ⎦
⎣𝑣𝑣̇ 3 ⎦ ⎢ 𝑘𝑘2 𝑐𝑐2 𝑘𝑘2 𝑐𝑐2 3
⎢ 0 0 − − ⎥
⎣ 𝑚𝑚3 𝑚𝑚3 𝑚𝑚3 𝑚𝑚3 ⎦

If we consider 𝑚𝑚1 = 𝑚𝑚2 = 𝑚𝑚3 = 1, 𝑐𝑐1 = 𝑐𝑐2 = 𝑐𝑐3 = 0, and, 𝑘𝑘1 = 𝑘𝑘2 = 𝑘𝑘3 = 1, for simplicity, then the
state equation becomes,
𝑥𝑥̇1 0 1 0 0 0 0 𝑥𝑥1
⎡𝑣𝑣̇ ⎤ ⎡
1 −1 0 1 0 0 0⎤ ⎡𝑣𝑣1 ⎤
⎢ ⎥ ⎢ ⎥⎢ ⎥
⎢𝑥𝑥̇ 2 ⎥ = ⎢ 0 0 0 1 0 0⎥ ⎢𝑥𝑥2 ⎥
⎢𝑣𝑣̇ 2 ⎥ ⎢ 1 0 −2 0 1 0⎥ ⎢𝑣𝑣2 ⎥
⎢𝑥𝑥̇ 3 ⎥ ⎢ 0 0 0 0 0 1⎥ ⎢𝑥𝑥3 ⎥
⎣𝑣𝑣̇ 3 ⎦ ⎣ 0 0 1 0 −1 0⎦ ⎣𝑣𝑣3 ⎦

The right Eigen matrix, 𝑉𝑉, the left Eigen matrix, 𝑊𝑊, and the Eigen matrix, Λ, of this state matrix is as shown
below.
0 + 1.7321𝑖𝑖 0 0 0 0 0
⎡ 0 0 − 1.7321𝑖𝑖 0 0 0 0 ⎤
⎢ ⎥
0 0 0 + 1𝑖𝑖 0 0 0
Λ=⎢ ⎥
⎢ 0 0 0 0 − 1𝑖𝑖 0 0 ⎥
⎢ 0 0 0 0 −9.3823 × 10−9 0 ⎥
⎣ 0 0 1 0 0 −9.3823 × 10 −9 ⎦

0.2041i − 0.2041i 0.5 0.5 −0.5774 0.5774

⎡ −0.3536 −0.3536 0.5i − 0.5i 0 0 ⎤
⎢ ⎥
− 0.4082i 0.4082i 0 0 −0.5774 0.5774⎥
V=⎢
⎢ 0.7071 0.7071 0 0 0 0 ⎥
⎢ 0.2041i 0.2041i −0.5 −0.5 −0.5774 0.5774⎥
⎣ −0.3536 −0.3536 − 0.5i 0 0 0 ⎦

−0.3536 −0.2041i 0.7071 − 0.4082i −0.3536 0.2041i

⎡−0.3536 0.2041i 0.7071 0.4082i −0.3536 −0.2041i⎤
⎢ ⎥
0.5000 0.5000i 0 0 −0.5000 −0.5000i⎥
W=⎢
⎢ 0.5000 0.5000i 0 0 −0.5000 0.5000i ⎥
⎢ 0 0.57735 0 0.5774 0 0.57735 ⎥
⎣ 0 0.57735 0 0.5774 0 0.57735 ⎦
If the right Eigen vector values for displacement states, and velocity states are extracted separately, in order
to identify the displacements and velocities in the three modes of oscillation it is seen that the corresponding
right Eigen values are as shown below and the corresponding pure modes of oscillations are as shown in
the below figures.
 0.2041i 0.5 −0.5774
𝑉𝑉𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 = �− 0.4082i 0 −0.5774�
0.2041i −0.5 −0.5774
−0.3536 0.5i 0
𝑉𝑉𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 = � 0.7071 0 0�
−0.3536 − 0.5i 0
Any oscillation of these three mass system will comprises of a mixture of these three pure modes of
oscillation in various percentages. The actual composition of a particular oscillation is dependent on the
initial excitation which is described by the left Eigen matrix and the corresponding left Eigen vectors.
We can use practical values for mass, spring constant and the damping coefficient to find that the mode
shapes have different magnitudes and phase angles which describe the oscillation of the three masses when
exited in each mode.

The below MATLAB script and the corresponding mode shapes are good example for such.
m1=1.5;
m2=1;
m3=1.25;

k1=0.75;
k2=1;
k3=1.25;

c1=0.1;
c2=0.15;
c3=0.1;

A=[0 1 0 0 0 0
-k1/m1 -c1/m1 k1/m1 c1/m1 0 0
0 0 0 1 0 0
k1/m2 c1/m2 -(k1+k2)/m2 -(c1+c2)/m2 k2/m2 c2/m2
0 0 0 0 0 1
0 0 k2/m3 c2/m3 -k2/m3 -c2/m3]
[V,lambda,W]=eig(A)
V_r=real(V);
V_i=imag(V);

The mode shapes for the dominant two modes of oscillation are as shown below.