Average Power Flow through a Synchronous Generator

Sam Chevalier

Guiding Question: When a synchronous generator experiences a forced oscillation, what

happens to the average power flow?

1 Synchronous Generator Model

We present the dynamic equations for a synchronous generator1 with the following characteristics.
• 2 pole
• salient-pole
• 3-phase
• wye-connected
• symmetric (balanced)
• 1 field winding and 3 damper windings on the rotor
• 3 identical stator windings displaced by exactly 120◦
The physical structure of this model is visualized in Figure 1. We define direct (d) and quadrature (q) axes
which are in line with and perpendicular to the axis of rotation of the rotor, respectively. We also arbitrarily
define δr which represents the angle between the fixed a phase axis of the stator and the direction of the q
axis.

𝑎𝑠

𝑞 axis

𝑐𝑠 ′ 𝑏𝑠 ′

𝑘𝑞2′ 𝛿𝑟
𝑘𝑑 ′
𝑘𝑞1′

𝑓𝑑 ′

𝑘𝑑 𝑎𝑠 axis

𝑓𝑑 𝑘𝑞2
𝑘𝑞1

𝑏𝑠 𝑐𝑠

𝑑 axis

𝑎𝑠 ′

Figure 1: 3-Phase Salient Pole Synchronous Generator (positive current flows into  and out of )

1 Equations and some figures adapted from Analysis of Electric Machinery by Paul Krause

1
Using Figure 1 as a guide, a simple circuit diagram may be employed to model the interior voltages and
currents. The system is modeled as though each rotor and stator coil acts as a transformer coil, capable of
transforming flux from the air gap into a voltage potential according to Faraday’s Law.

Rotor
Damper 1 Damper 2
+ 𝑣𝑘𝑞1 − + 𝑣𝑘𝑞2 −
Stator 𝑟𝑘𝑞1 𝑖𝑘𝑞1 𝑟𝑘𝑞2 𝑖𝑘𝑞2
𝑖𝑏𝑠 + 𝑟𝑠 𝑁𝑘𝑞1 𝑁𝑘𝑞2
𝑖𝑘𝑑
𝑖𝑐𝑠 + 𝑟𝑠
−
𝑖𝑎𝑠 + 𝑟𝑠
𝑁𝑠 𝑁𝑠 𝑣𝑘𝑑 𝑁𝑘𝑑 Damper 3
−
𝑁𝑠 +
𝑟𝑘𝑑
𝑖𝑓𝑑
−

𝑣𝑓𝑑 𝑁𝑓𝑑 Field Winding

+
𝑟𝑓𝑑

Figure 2: Circuit Model of the 3-Phase Synchronous Generator. Each of the 7 coils has an associated
resistance rx and coil turn count Nx (with a specified reactance).

Transformations
In order to simplify the circuit analysis, we make two key transformations. First, we reflect all rotor voltages,
currents, and impedances to the stator reference frame. In general, way may transform a voltage from a
Ni
secondary side of a transformer to the primary by noting that vi = N j
vj . Similar relationships exit for
current and impedance transformations. In a circuit diagram, reflected values all contain a prime (0 ) symbol.
This can be seen in Figure 3. The second transformation we employ is known as Park’s transformation. In
steady state, the rotor spins at synchronous frequency ωs . In order to over come the mathematical challenge
of a rotating rotor, we choose to analyze the stator’s circuit in the synchronously rotating reference frame of
the rotor. This eliminates the time varying inductance values of the system since all relationships are now
stationary. For example, once the q-axis, d-axis, and zero sequence voltages have been calculated, we can
convert them back into physically meaningful phase a, b, and c voltages (or currents, fluxes, etc.).
    
vas cos (δr )
sin (δr )
1 vqs
 vbs  =  cos δr − 2π 2π
3
sin δr − 3
1
  vds  (1)
2π 2π
vcs cos δr + 3 sin δr + 3 1 v0s

Transformed Circuit Model

Without further derivation, we present the equivalent circuit model for a 3-phase synchronous generator
which has undergone both transformations (transformer reflection and Park’s) described in the previous
section. Table 1 interprets the circuit diagram.

2
 Subscript Significance Superscript Significance
0
fd field winding 1 in d-axis transformed across transformer
kd damper 1 in d-axis r rotor reference frame
kq1 damper 1 in q-axis
kq2 damper 2 in q-axis
m magnetizing
l leakage

Table 1: Figure 3 Interpretation

′𝑟
Q Axis ′
+ 𝑣𝑘𝑞2
𝑟𝑘𝑞2 ′𝑟
𝑖𝑘𝑞2
𝐿′𝑙𝑘𝑞2
𝑟 λ𝑟𝑑𝑠 𝜔𝑟
𝑖𝑞𝑠 𝑟𝑠 𝐿𝑙𝑠
+ + − 𝐿′𝑙𝑘𝑞1
𝑟
𝑣𝑞𝑠 ′
𝐿𝑚𝑞 𝑟𝑘𝑞1
Zero Sequence
−
′𝑟 + ′𝑟
𝑟𝑠
𝑖𝑘𝑞1 𝑣𝑘𝑞1 +

𝑣0𝑠 𝐿𝑙𝑠
𝑖0𝑠
′𝑟 −
+ 𝑣𝑘𝑑
D Axis ′
𝑟𝑘𝑑 ′𝑟
𝑖𝑘𝑑
𝐿′𝑙𝑘𝑑
𝑟 λ𝑟𝑞𝑠 𝜔𝑟
𝑖𝑑𝑠 𝑟𝑠 𝐿𝑙𝑠
+ − + 𝐿′𝑙𝑓𝑑
𝑟
𝑣𝑑𝑠 ′
𝑟𝑓𝑑
𝐿𝑚𝑑
−
′𝑟 + ′𝑟
𝑖𝑓𝑑 𝑣𝑓𝑑

Figure 3: Transformed Circuit Model of the 3-Phase Synchronous Generator. λ is the flux linkage between
the rotor and stator coils. The product λω represents a voltage drop in the d-axis (where the field winding
is located) and a voltage increase in the q-axis. If the system is feeding balanced load, then v0s = 0 and the
zero sequence circuit can be omitted.

Although a very busy diagram, the circuit simplifies considerably in steady state since the damper winding
currents will be equal to 0. This reduced circuit is shown in Figure 4.

3
 Q Axis

𝑟 λ𝑟𝑑𝑠 𝜔𝑟
𝑖𝑞𝑠 𝑟𝑠 𝐿𝑙𝑠
+ + −
𝑟
𝑣𝑞𝑠
𝐿𝑚𝑞
−

D Axis

𝑟 λ𝑟𝑞𝑠 𝜔𝑟
𝑖𝑑𝑠 𝑟𝑠 𝐿𝑙𝑠
+ − + 𝐿′𝑙𝑓𝑑
𝑟
𝑣𝑑𝑠 ′
𝑟𝑓𝑑
𝐿𝑚𝑑
−
′𝑟 + ′𝑟
𝑖𝑓𝑑 𝑣𝑓𝑑

Figure 4: Steady State Circuit Model of the Transformed 3-Phase Synchronous Generator. Damper winding
effects have been removed.

Transformed Circuit Model State Equations

Using basic KVL, KCL and Newton’s Second Law, state equations can be formulated which describe the
dynamics of the full synchronous generator model given in Figure 3. Without proof, we note a simplified
expression for the electrical torque τe .
3
τe = (λds iqs − λqs ids ) (2)
2
This expression for electrical torque is used in constructing the system’s swing equation, where damping is
omitted. The variable ωr represents the true angular velocity of the rotor while δr represents the rotors
absolute angular position (not with respect to a synchronously rotating reference).

τm = J ω̇r + τe (3)
ωr = δ̇r (4)

Using the flux linkage variables (λ), we can write out the seven defined voltages in Figure 3 using basic
circuit theory. First we define the three stator voltages, where p denotes the derivative operator with respect
to time.

vqs = −rs iqs + ωr λds + pλqs (5)

vds = −rs ids − ωr λqs + pλds (6)
v0s = −rs i0s + pλ0s (7)

Next we define the four rotor voltages.

0 0 0 0
vkq1 = rkq1 ikq1 + pλkq1 (8)
0 0 0 0
vkq2 = rkq2 ikq2 + pλkq2 (9)
0 0 0 0
vf d = rf d if d + pλf d (10)
0 0 0 0
vkd = rkd ikd + pλkd (11)

The d-axis and q-axis current flows are given below.

4
 1
iqs = − (λqs − λmq ) (12)
Lls
1
ids =− (λds − λmd ) (13)
Lls
The flux linkage terms associated with each winding on the d-axis and q-axis are defined, again, using basic
circuit theory.
 0 0

λqs = −Lls iqs + Lmq −iqs + ikq1 + ikq2 (14)
 0 0

λds = −Lls ids + Lmd −ids + ikd + if d (15)
λ0s = −Lls i0s (16)
0 0 0
 0 0

λkq1 = Llkq1 ikq1 + Lmq −iqs + ikq1 + ikq2 (17)
0 0 0
 0 0

λkq2 = Llkq2 ikq2 + Lmq −iqs + ikq1 + ikq2 (18)
0 0 0
 0 0

λf d = Llf d if d + Lmd −ids + if d + ikd (19)
0 0 0
 0 0

λkd = Llkd ikd + Lmd −ids + if d + ikd (20)

It is also useful to define the flux linkages associated with the Lmq and Lmd inductors. We do so by balancing
the currents via KCL, but we write this current balance in terms of inductances and the defined flux linkages.

0 0 !
λqs λkq1 λkq2
λmq = Laq + 0 + 0 (21)
Lls Llkq1 Llkq2
0 0
!
λds λf d λkd
λmd = Lad + 0 + 0 (22)
Lls Llf d Llkd

In this case, it is clear that Laq and Lad are effective reactance parameters. These are defined as the parallel
combination of all reactances in each respective (q-axis and d-axis circuit).
!−1
1 1 1 1
Laq = + + 0 + 0 (23)
Lmq Lls Llkq1 Llkq2
!−1
1 1 1 1
Lad = + + 0 + 0 (24)
Lmq Lls Llf d Llkd

The flux linkage per second variables ψ can be used at the primary state variables, along with θr and ωr ,
which drive the dynamics of the system. Therefore, we must calculate their derivatives. This is accomplished
through Faraday’s Law and basic algebraic manipulation. The results are shown below.

5
 rs
λ̇rqs = vqs
r
− ωr λrds + λrmq − λrqs


(25)
Lls
rs
λ̇rds r r
= vds + ωr λqs + (λr − λrds ) (26)
Lls md
rs
λ̇0s = v0s − λ0s (27)
Lls
0
0
r 0
r
rkq1  r 0
r

λ̇kq1 = vkq1 + 0 λmq − λkq1 (28)
Llkq1
0
0
r 0
r
rkq2  r 0
r

λ̇kq2 = vkq2 + 0 λmq − λkq2 (29)
Llkq2
0
0 0 rf d  r 0

λ̇frd r
= vf d + 0 λmd − λfrd (30)
Llf d
0
r0 0
r rkd  r 0
r

λ̇kd = vkd + 0 λmd − λkd (31)
Llkd

At this point, the synchronous generator may be fully simulated using the state equations which have been
developed. There are a variety of methods for simulating this system, but the DC field voltage, the applied
mechanical torque, and the AC stator voltages are typically defined as external model inputs. Stator currents
are the external outputs. The procedure for simulating this system is outlined in Figure 5 (a similar procedure
is vaguely recommended in Analysis of Electric Machinery by Paul Krause)

Variable Key
𝑥ҧ = 𝑥𝑡𝑖−1
𝑥 = 𝑥𝑡𝑖
Input Output
𝑣𝑎𝑠
ҧ 𝑣𝑞𝑠
ҧ 𝑖𝑎𝑠
Forward Euler:
𝑣ҧ𝑏𝑠 𝑣ҧ𝑑𝑠 𝜆ሶ ҧ𝑞𝑑0
Park Compute Internal Compute New Inverse Park 𝑖𝑏𝑠
𝑣𝑐𝑠
ҧ Transformation 𝑣ҧ0𝑠 Flux Derivatives Flux States Transformation 𝑖𝑐𝑠
𝜆 = 𝜆ሶ ҧ + (∆𝑡)𝜆ҧ
𝑣𝑓𝑑
𝜆𝑞𝑑0
𝜏𝑒
Compute
𝑖𝑞𝑑0
𝜏𝑚 Swing Electrical Torque Compute Internal 𝑖𝑞𝑑0
Equation Current States 𝜔𝑟

𝛿𝑟ҧ
Forward Euler:
𝛿𝑟 𝛿𝑟
𝜔
ഥ𝑟 Rotor Angle
𝛿𝑟 = 𝛿𝑟ҧ + (∆𝑡)𝜔
ഥ𝑟

Figure 5: Simulation Procedure

2 Modeling the Interconnection of 2 Synchronous Generators:

Network Dynamics Excluded
We wish to model the interconnection of 2 synchronous generators in order to further trace how oscillations
propagate through the system. We do so by considering the diagram shown in Figure 6.

6
 (A)
Stator 1 3 Lines Stator 2
𝑟𝑠 𝑖𝑏𝑠1 𝑖𝑏𝑠2 + 𝑟𝑠

𝑟𝑠 𝑖𝑐𝑠1 𝑖𝑐𝑠2 + 𝑟𝑠
𝑟𝑠 𝑖𝑎𝑠1 𝑖𝑎𝑠2 + 𝑟𝑠
𝑁𝑠 𝑁𝑠 𝑁𝑠 𝑁𝑠
− −
𝑁𝑠 𝑁𝑠

(B)
Stator 1 Equivalent Line Stator 2

𝑟𝑠 𝑖𝑏𝑠1 V1 𝑒
𝑗𝜃1 V2 𝑒 𝑗𝜃2 𝑖 𝑟𝑠
𝑏𝑠2 +

𝑟𝑠 𝑖𝑐𝑠1 𝑖𝑐𝑠2 + 𝑟𝑠
𝑟𝑠 𝑖𝑎𝑠1 𝑖𝑎𝑠2 + 𝑟𝑠
𝑁𝑠 𝑁𝑠 𝑁𝑠 𝑁𝑠
− −
𝑁𝑠 𝑁𝑠

Figure 6: Double Synchronous Generator Interconnection. Panel (A) shows the full interconnection while
panel (B) shows the simplified connection.

We model the dynamics of each generator individually, and we use a set of algebraic relationships to couple
these dynamics. We first re-state the simplified flux dynamics of a 2 pole machine. Damper winding voltages
are excluded since they always have a value of 0, and we assume a balanced system, so we eliminate the 0
sequence terms entirely.
rs
λ̇qs = −ωr λds + (λmq − λqs ) + vqs
Lls
rs
λ̇ds = ωr λqs + (λmd − λds ) + vds
Lls
0
0 rkq1  0

λ̇kq1 = 0 λmq − λkq1
Llkq1
0
0 rkq2  0

λ̇kq2 = 0 λmq − λkq2
Llkq2
0
0 r 
fd 0
 0
λ̇f d = 0 λrmd − λfrd + vfrd
Llf d
0
0 rkd  0

λ̇kd = 0 λmd − λkd
Llkd
τm − τe
ω̇ =
J
˙δ = ω − w
r r s
=ω

We treat the exciter voltage vf d and the generator’s mechanical torque τm as independent inputs. The flux

7
variables λmq and λmd are defined through the following algebraic relationships.
0 0 !
λqs λkq1 λkq2
λmq = Laq + 0 + 0
Lls Llkq1 Llkq2
0 0
!
λds λf d λkd
λmd = Lad + 0 + 0
Lls Llf d Llkd

The unknown voltage inputs (vqs and vds ) may be solved for in terms of flux state variables and known
inductance terms.

vqs = −rs iqs + ωr λds + λ̇qs

vds = −rs ids − ωr λqs + λ̇ds

We constrain these differential equations with the following set of algebraic equations.
3
P = (ids vds + iqs vqs )
2
3
Q = (ids vqs − iqs vds )
2
3
τe = (λds iqs − λqs ids )
2
The unknown voltage inputs (iqs and ids ) may be solved for in terms of flux state variables and known
inductance terms.
1
iqs = − (λqs − λmq )
Lls
1
ids =− (λds − λmd )
Lls
One last highly important relationship exists between the generator and the network buses. We define the bus
voltage phasor V ejθ as the bus voltage at the node interfacing the generator with the system. Additionally,
δ r is the angle of the rotor with respect to some arbitrary synchronously rotating reference frame.
√

vds = 2V sin δ r − θ
√

vqs = 2V cos δ r − θ

Finally, we may employ the AC power flow equations to approximate the power transfer between the genera-
tors. Since phasor analysis is applied to the transmission line power flow, the line reactance is not considered
to be dependent on the system frequency. Assuming small, low frequency power oscillations, this assumptions
should have almost no effect on the analysis.

P = Vi2 Gii + Vi Vj [Gij cos(θij ) + Bij sin(θij )]

3
(ids vds + iqs vqs ) = Vi2 Gii + Vi Vj [Gij cos(θij ) + Bij sin(θij )]
2
Q = −Vi2 Bii + Vi Vj [Gij sin(θij ) − Bij cos(θij )]
3
(ids vqs − iqs vds ) = −Vi2 Bii + Vi Vj [Gij sin(θij ) − Bij cos(θij )]
2
In simulating this system, we assume both generators are identical. The algebraic and differential equations

8
for each generator are summarized below.

0 = Vi2 Gii + Vi Vj [Gij cos(θij ) + Bij sin(θij )] − 23 (ids vds + iqs vqs )


0 = −Vi2 Bii + Vi Vj [Gij sin(θij ) − Bij cos(θij )] − 32 (ids vqs − iqs vds )




3
2 (λds iqs − λqs i
ds ) − τe



 0=√




 0 = √ 2V sin δ r − θ
− vds
 0 = 2V cos δr − θ − vqs


Algebraic Equations = 0 = − L1ls (λds − λmd ) − ids
1
 0 = − Lls(λqs − λ0mq ) − iqs 



0

λf d
 λkd
 λds
 0=L


 ad L ls
+ 0
Llf d
+ 0
Llkd
− λmd

  0 0

 0 = Laq λLqs + λ0kq1 + λ0kq2 − λmq



ls Llkq1 Llkq2
 rs

 λ̇qs = −ωr λds + Lls (λmq − λqs ) + vqs
rs
λ̇ds = ωr λqs + Lls (λmd − λds ) + vds




 0  
 0 rkq1 0
λ̇ = λ − λ



 kq1 0
Llkq1 mq kq1

0

  
 0 rkq2 0
 λ̇kq2 = 0

Llkq2
λ mq − λ kq2
Differential Equations = 0  
0 rf d r 0
r 0
r
λ̇ = λ − λ f d + vf d

 0
 fd md

 Llf d
 0  
 λ̇0 = r0kd λmd − λ0




 kd Llkd kd
 τm −τe


 ω̇ = J
 ˙

δ r = ωr − ws = ω

3 Modeling the Interconnection of 2 Synchronous Generators and

an Infinite Bus: Network Dynamics Included
In order to model the interconnection of 2 synchronous generators and an infinite bus, we consider the
network diagram depicted by figure 7. In this model, network dynamics are NOT excluded, so we no longer
are able to employ the steady state power flow equations to model the transfer of power.

𝑟∞ 𝐿∞
𝑣∞𝑎 𝑖∞ 𝑎
𝑣∞𝑏 𝑖 ∞𝑏
𝑣∞𝑐 𝑖∞ 𝑐

𝑟𝑙 𝐿𝑙
𝑟𝑠 𝑖1𝑎 𝑖2 𝑎 + 𝑟𝑠

𝑟𝑠 𝑖1𝑏 𝑖2 𝑏 + 𝑟𝑠
𝑟𝑠 𝑖1𝑐 𝑖2𝑐 + 𝑟𝑠
𝑁𝑠 𝑁𝑠 𝑁𝑠 𝑁𝑠
− −
𝑁𝑠 𝑁𝑠

Stator 1 Stator 2

Figure 7: Double Synchronous Generator & Infinite Bus System

9
Analysis of a Multi-Generator System with qd0 Network Transfor-
mation
In order to overcome the challenges of analyzing a network in the abc sequence while analyzing each machine
in qd0 sequence (with each sequence in its own reference frame too), we choose to transform the network
such that we mathematically eliminate all abc terminal voltages and currents. This approach is given by
Figure 8.

Analyze Network in qd0 sequence Analyze Generator in

of synchronous reference rotor reference frame

𝑒 𝑒 𝑟
𝑖𝑎𝑏𝑐 𝑖𝑞𝑑0 𝑖𝑞𝑑0 𝑖𝑞𝑑0
(𝐾𝑠 )−1 Network (𝐾 𝑟 )−1 Generator
𝑒 𝑒 𝑟
𝑣𝑎𝑏𝑐 𝑣𝑞𝑑0 𝑣𝑞𝑑0 𝑣𝑞𝑑0

Inverse Park Transform Transform from rotor to

to for post analysis synchronous reference frame

Figure 8: Reference Frame Transformations

To perform this analysis, we first consider the balanced 3 phase network depicted in Figure 9 where each
voltage and current variable represents an instantaneous time domain value.

𝑣1𝑎 𝑟𝑙 𝐿𝑙 𝑣2𝑎

𝑖𝑎
𝑣1𝑏 𝑟𝑙 𝐿𝑙 𝑣2𝑏

𝑖𝑏
𝑣1𝑐 𝑟𝑙 𝐿𝑙 𝑣2𝑐

𝑖𝑐

Figure 9: Reference Frame Transformations

We write out the network voltage equations.

dia
v2a = v1a − ia rl − Ll
dt
dib
v2b = v1b − ib rl − Ll
dt
dic
v2c = v1c − ic rl − Ll
dt
We then apply a Park’s Transformation Ks (θe ) to the system equations, where θe is chosen to be the
instantaneous angle of the synchronously rotating reference frame.

v qd0 = Ks v abc
cos θe − 2π cos θe + 2π



  cos (θe )
2  3
3

Ks = sin (θe ) sin θe − 2π
3 sin θe + 2π
3

3 1 1 1
2 2 2

We choose to apply this transformation to the entire network, where r̂l and L̂l are diagonal resistance and

10
inductance matrices of the network, respectively.
diabc
Ks v 2abc = Ks v 1abc − Ks r̂l iabc − Ks L̂l
dt
d
v 2qd0 = v 1qd0 − Ks r̂l Ks−1 Ks iabc − Ks L̂l Ks−1 Ks iabc



dt
d
v 2qd0 = v 1qd0 − Ks r̂l Ks−1 iqd0 − Ks L̂l Ks−1 iqd0



dt
If the network is balanced, than we have that r̂l = rl I3 , where I3 is the 3 by 3 identity matrix. The same is
true for the inductance matrix.
Ks r̂l Ks−1 = rl Ks I3 Ks−1
= r̂l
Ks L̂l Ks−1 = Ll Ks I3 Ks−1
= L̂l
We apply these identities, along with the product differentiation rule, to the network voltage expression.
   
d −1 d
v 2qd0 = v 1qd0 − r̂l iqd0 − Ks L̂l Ks iqd0 + Ks−1 iqd0
dt dt
   
d −1 −1 d
= v 1qd0 − r̂l iqd0 − L̂l Ks Ks iqd0 + Ks Ks iqd0
dt dt
  
d d
= v 1qd0 − r̂l iqd0 − L̂l Ks Ks−1 iqd0 + iqd0
dt dt
Without showing the algebra (see Krause), we make the following statement.
 
0 1 0
d
Ks Ks−1 = ωe  −1 0 0 
dt
0 0 0
 >
Therefore, we arrive at the following simplified expression, where idq0 = id −iq 0 .
 
d
v 2qd0 = v 1qd0 − r̂l iqd0 − L̂l ωe idq0 + iqd0
dt
We can write this out explicitly for the voltage of each sequence. Since we assume a balanced system, we
neglect the 0 sequence voltages and currents.
diq
v2q = v1q − rl iq − ωe Ll id − Ll (32)
dt
did
v2d = v1d − rl id + ωe Ll iq − Ll (33)
dt
Using (32) and (33) as a guide, we may write the network equations for the system shown in 7. We first
consider the connection between the infinite bus and generator 1.
di∞q
v1q = v∞q − r∞ i∞q − ωe L∞ i∞d − L∞
dt
di∞d
v1d = v∞d − r∞ i∞d + ωe L∞ i∞q − L∞
dt
Since infinite bus currents are not known explicitly, we write them in terms of i1 and i2 since we know that
i1 + i2 + i∞ = 0. Therefore, i∞ = −i1 − i2 .
 

di1q di2q
v1q = v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + L∞ + (34)
dt dt
 

di1d di2d
v1d = v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q + L∞ + (35)
dt dt

11
We next consider the connection between bus 1 and bus 2.
di2q
v2q = v1q + rl i2q + ωe Ll i2d + Ll (36)
dt
di2d
v2d = v1d + rl i2d − ωe Ll i2q + Ll (37)
dt
Therefore, equations (34), (35), (36), and (37) define the network equations which are necessary to simulate
the system. We note that in order to determine v∞q and v∞d , we are required to apply a Park’s Transfor-
mation on the infinite bus voltages. Since, by definition, the infinite bus voltage frequency is constant for
all time, θe = ωe t.
   
v∞q v∞ a
 v∞d  = Ks  v∞b 
v∞ 0 v∞c
  cos (θe ) cos θe − 2π cos θe + 2π


 
v∞a
2  3
3

= sin (θe ) sin θe − 2π
3 sin θe + 2π
3
  v∞ b 
3 1 1 1
2 2 2
v∞c

To relate the network variables to the generator variables, we must define the transformation matrix r K e
which transforms voltages and currents in the rotor reference frame to voltages and currents in the syn-
chronous reference frame (or vice versa, depending on the need). We define δr as the rotor angle for a given
generator.

v eqd0 =r K e v rqd0
 
cos(θe − δr ) − sin(θe − δr ) 0
r e
K =  sin(θe − δr ) cos(θe − δr ) 0 
0 0 1

Moving forward, we denote v eqd0 = v qd0 with no superscript. These are the voltages (or currents) in the
synchronous reference frame. Alternatively, we denote v rqd0 as the voltages (or currents) in the rotor reference
frame. Finally, we may specify the differential equations and the accompanying algebraic constraints which
are necessary to simulate each generator along with the network interconnections.

12
 0 = [r K e ] irqd0 − iqd0


0 = [r K e ] v rqd0 − v qd0




0 = − L1ls (λds − λmd ) − irds





0 = − L1ls (λqs − λmq )
− irqs



Generator Algebraic Equations = 0 = 23 λ r r
ds iqs − λqs ids − τe (38)
 0 0
λf d
 λkd
λds
− λmd




 0 = Lad Lls + 0
Llf d
+ 0
Llkd

  0 0

λkq1 λkq2
 λqs
− λmq


 0 = Laq Lls + 0
Llkq1
+ 0
Llkq2

λ̇qs = −ωr λds + Lrls

 r


s
(λmq − λqs ) + vqs
rs r
λ̇ds = ωr λqs + Lls (λmd − λds ) + vds




 0  
r

 0 0


 λ̇kq1 = L0kq1 λmq − λkq1
lkq1


 0  
 0 rkq2 0
λ̇kq2 = L0 λmq − λkq2


Generator Differential Equations = 0
lkq2
  (39)
0 rf d 0 0
λmd − λf d + vf d



 λ̇f d = L0
 lf d
0

  
 0 r 0


 λ̇kd = L0kd λmd − λkd
lkd

ω̇r = τmJ−τe





δ̇r = ωr

 

di1q di2q
v 1 = v ∞ + r∞ i 1 + i2 + ωe L ∞ (i1 + i2 ) + L ∞ +
q q q q d d
 dt dt 




di di
v1d = v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q + L∞ dt1d + dt2d

Network Equations =
di2
v = v1q + rl i2q + ωe Ll i2d + Ll dtq

 2q


 di
v2d = v1d + rl i2d − ωe Ll i2q + Ll dt2d
(40)
0 0 0 0
r r
Specifically, for each generator, we have 10 state variables: λqs , λds , λkq1 , λkq2 , λf d , λkd , vqs , vds , δr , and
ωr . Therefore, each generator adds ten additional state variables. The network equations, in fact, do not
introduce any new state variables: they simply represent a set of constraints on the aforementioned state
variables.

Steady State Power Flow Solution

Before attempting to simulate this entire system, we first choose to write out the solution for the system in
steady state. To do so, we treat the infinite bus as if it was a swing bus (where RMS voltage magnitude and
phase angle are both specified), and we treat both generators like PV buses, where RMS voltage magnitude
and real power are both specified (we could alternatively specify reactive power instead of voltage magnitude).
In the time domain, we note that the active power injections have the following definitions.

P1 = p1 (t)
= i1a v1a + i1b v1b + i1c v1c
P2 = p2 (t)
= i2a v2a + i2b v2b + i2c v2c
P∞ = p3 (t)
= i∞a v∞a + i∞b v∞b + i∞c v∞c

Since this system is operating in steady state, we may use phasor analysis on each phase of the generator

13
where V = Vejδ .
 ∗  ∗  ∗
V 1a − V 2a V1a − V∞a V1b − V2b V1b − V∞b V1c − V2c V1c − V∞c
P1 + jQ1 = V1a + + V1b + + V1c +
rl + jωs ll r∞ + jωs l∞ rl + jωs ll r∞ + jωs l∞ rl + jωs ll r∞ + jωs l∞
 ∗  ∗  ∗
V 2a − V 1a V 2b − V 1b V2c − V1c
P2 + jQ2 = V2a + V 2b + V2c
rl + jωs ll rl + jωs ll rl + jωs ll
 ∗  ∗  ∗
V∞a − V1a V∞b − V1b V∞c − V1c
P∞ + jQ∞ = V∞a + V∞b + V∞c
r∞ + jωs l∞ r∞ + jωs l∞ r∞ + jωs l∞

If we separate these equations into their real and imaginary parts, we obtain 6 equations with 6 unknowns:
δs1 , δs2 , Q1 , Q2 , P∞ , and Q∞ (we arbitrarily specify δsi = 0). Therefore, an iterative method may be used
to solve for these variables. This is the first step used to determine the initial conditions for this system.
Once the unknown values have been solved for, the terminal abc voltages and currents can be solved for. For
example, the phase a voltages on generator 1 is given below.
√
v1a (t) = 2 |V1a | cos(2πωs t + δs1 )

To determine the phase a current on generator 1, we solve the current phasor expression and then extrapolate
to the time domain.
V1a − V2a V1a − V∞a
I 1a = +
rl + jωs ll r∞ + jωs l∞
√
i1a (t) = 2 |I1a | cos(2πωs t + ∠I1a )

Initializing Other State Variables

Once power flow is solved, we initialize the other state variables by setting all derivatives equal to 0 (except
for the rotor angles, since their derivatives are equal to synchronous frequency) and solving the associated
0
system of equations. We are particularly interested in solving for τm , the input mechanical torque and vf d ,
the field voltage (transformed across the machine windings). Both of these parameters represent explicit
control settings which are adjusted by generator operators.

Generator Power and Energy Analysis

To consider how a forced oscillation affects the power and energy oscillations in the system, we first balance
the powers in the stationary circuits given by figure (10).

′ ′
+ 𝑣𝑘𝑞2 + 𝑣𝑘𝑑
Q Axis ′
𝑟𝑘𝑞2 D Axis ′
𝑟𝑘𝑑 ′
′ 𝑖𝑘𝑑
𝑖𝑘𝑞2 𝐿′𝑙𝑘𝑑
𝐿′𝑙𝑘𝑞2
λ𝑑𝑠 𝜔𝑟 𝑟 λ𝑞𝑠 𝜔𝑟
𝑟
𝑖𝑞𝑠 𝑟𝑠 𝐿𝑙𝑠 𝑖𝑑𝑠 𝑟𝑠 𝐿𝑙𝑠
+ + − 𝐿′𝑙𝑘𝑞1 + − + 𝐿′𝑙𝑓𝑑
𝑟 𝑟
𝑣𝑞𝑠 ′
𝑟𝑘𝑞1 𝑣𝑑𝑠 ′
𝑟𝑓𝑑
𝐿𝑚𝑞 𝐿𝑚𝑑
− −
′ + ′ ′ + ′
𝑖𝑘𝑞1 𝑣𝑘𝑞1 𝑖𝑓𝑑 𝑣𝑓𝑑

Figure 10: Transformed Circuit Model of the 3-Phase Synchronous Generator. The zero sequence circuit is
omitted. All variables are given in the rotor’s local reference frame, but the r superscript is only given to
the terminal voltages and currents.

14
According to Tellegen’s theorem, the instantaneous powers in each circuit must sum to 0. We now balance
the powers in each circuit, beginning with the Q-axis. Employing the passive sign convention (PSC), power
consumption is positive while power generation is negative.

 0 2 0  0 2 0
2
Resistive Dissipation ⇒ikq1 rkq1 + ikq2 rkq2 + irqs rs
0  0 d 0  0  0 d 0 
Inductor Power Injection ⇒ ikq1 Llkq1 i + ikq2 Llkq2 i +
dt kq1 dt kq2
   
d r d
irqs Lls


iqs + (imq ) Lmq imq
dt dt
Power Tranfer ⇒ − λds ωr
Terminal Power Output ⇒ irqs vqs
r

Terminal Power Input ⇒ 0

The D-axis circuit power may be balanced similarly, but this circuit contains a power input term (from the
field excitation).

 0 2 0  0 2 0
2
Resistive Dissipation ⇒ikd rkd + if d rf d + (irds ) rs
0  0 d 0  0  0 d 0 
Inductor Power Injection ⇒ ikd Llkd i + if d Llf d i +
dt kd dt f d
   
d r d
(irds ) Lls ids + (imd ) Lmd imd
dt dt
Power Tranfer ⇒ λqs ωr
Terminal Power Output ⇒ irds vds
r

0 0
Terminal Power Input ⇒ − if d vf d

Full Machine Power Balance

We next consider the power balance of the entire machine. To do so, we first consider the conservative forces
present in a single generator.
1 2
Rotational Kinetic Energy of the Turbine ⇒ δ̇ J
2 r
1
Magentic Potential Energy of the Coils ⇒ q̇> L(δr )q̇
2
where L(δr ) is given by (68). We can therefore state the Lagrangian.

L (xj , ẋj , qk q̇) = T ∗ − V + Wm

∗
− We
1 2 1 >
= δ̇r J − 0 + q̇ L(δr )q̇ − 0
2 2
For the generalized force, we determine Lagrange’s equation:
   
d ∂L ∂L
τm = −
dt ∂ δ˙r ∂δr
1
= δ¨r J − q̇> L0 (δr )q̇ (41)
2

15
And for the generalized voltages, we determine Lagrange’s equation. The voltages in (42) represent the
voltage drop across the inductors in Figure 2.
 
va + ia ra
 vb + ib rb 
 
 vc + ic rc 
L0 (δr )δ˙r q̇ + L(δr )q̈ = 
 
 vf d − if d rf d (42)


 −ikd rkd 
 
 −ikq1 rkq1 
−ikq2 rkq2

We may calculate the energy supplied by the shaft by noting that, in general, energy is equal to the
product of torque and angular displacement: E = τ (4θ).

Eshaft = τm (δr − δr0 )

where δr is the rotor angle as a function of time and δr0 is the initial angular position of the rotor. To
write this expression more carefully, we must consider the situation where τm is a function of time instead of
constant. In this case, we integrate the product of the shaft power and the derivative of the angular velocity.
Z tf
Eshaft = τm δ̇r dt
t0
d
We are further interested in the flow of energy, or the power, produced by the shaft, where Pshaft = dt Eshaft .

Pshaft = τm δ̇r
 
1
= δ¨r J − q̇> L0 (δr )q̇ δ̇r
2

The power provided by the turbine shaft represents one power input. The second power input is provided
by the product of the excitation field current and the excitation field voltage.
3 0 0
Pexc = v i
2 fd fd
= vf d i f d

Therefore, the system has the following total input power.

Pin = Pexc + Pshaft

The input power must be equal to the sum of four factors: 1) output power, 2) resistive dissipation, 3)
turbine inertia power, and 4) inductive power. The output power is defined as the product of the terminal
voltage and current (the following identity is shown in the Appendix).

Pout = v1a i1a + v1b i1b + v1c i1c

3 r r
vq1 iq1 + vdr1 ird1


=
2
Resistive dissipation consumes energy on both the rotor and stator side of the generator.

Pdiss = q̇> Rq̇

As the turbine accelerates and decelerates, it injects positive and negative power into the system. The
rotational kinetic energy of the turbine has the following value.
1 2
Eturb = Jω
2 r

16
The derivative of energy yields and expression for turbine power. Clearly, a positive power value corresponds
to a power consumption.
d
Pturb = Eturb
dt
= Jωr ω̇r

Finally, we must consider the instantaneous power consumed by the inductors. We note that the energy
stored in the inductors has the following expression.
1 >
EL =
q̇ L(δr )q̇
2
The derivative of this expression will compute the power injection of the inductors.
d
PL = EL
dt
1 >
q̇ L(δr )q̈ + L̇(δr )q̇ + q̈> L(δr )q̇


=
2
1 >
q̇ L(δr )q̈ + q̇> L̇(δr )q̇ + q̈> L(δr )q̇

=
2
If we use Park’s Transform to eliminate the abc variables, the inductor matrix becomes a static matrix, as
is shown in equation (69) in the appendix. We use Faraday’s law to determine the power drop across these
inductors.
  >
3 d
P̂L = λ̂ q̂˙
2 dt
  >
3 d ˙
= L̂(δr )q̂ q̂˙
2 dt
  >
3 ˙
= L̂(δr )¨q̂ + L̂(δr )q̂˙ q̂˙
2

Since L̂(δr ) is a static inductance matrix, its time derivative is 0. Therefore, the expression for inductor
power simplifies.
 
3 ¨>
P̂L = q̂ L̂(δr )> q̂˙
2

We may now state the power balance for an entire generator system, where PL = P̂L .

Pin = Pout + Pdiss + Pturb + PL

This may be stated mathematically.

3 1 >
(vq1 iq1 + vd1 id1 ) + q̇> Rq̇ + Jωr ω̇r + q̇ L(δr )q̈ + q̇> L̇(δr )q̇ + q̈> L(δr )q̇

vf d if d + τm ωr =
2 2
We note, though, that ω̇r J − 12 q̇> L0 (δr )q̇ = τm . We make this substitution.
 
1 > 0 3 1 >
vf d if d + ω̇r J − q̇ L (δr )q̇ ωr = (vq1 iq1 + vd1 id1 ) + q̇> Rq̇ + Jωr ω̇r + q̇ L(δr )q̈ + q̇> L̇(δr )q̇ + q̈> L(δr )q̇

2 2 2
1 3 1 >
vf d if d − q̇> L0 (δr )q̇ωr = (vq1 iq1 + vd1 id1 ) + q̇> Rq̇ + q̇ L(δr )q̈ + q̇> L̇(δr )q̇ + q̈> L(δr )q̇

2 2 2

17
Setting up the ODEs for Linearization
In order to perform linearization and eigenvalue analysis, we present a slightly alternative simulation approach
where we entirely eliminate the terminal qd axis voltages. This in turn eliminates the algebraic variables
present in the DAE system given by equations (38-40). The reason for doing so it as follows. In order to build
the Jacobian of the linearized the system, we wish to deal with a set of ODEs. Currently, our expressions
for q and d axis flux derivatives contain voltage terms. For example:
rs r
λ̇qs = −ωr λds + (λmq − λqs ) + vqs
Lls
In our system, since infinite bus voltage is fixed, we can express all unknown voltages in terms of currents
and the known infinite bus voltage. Deriving an expression for voltage in the rotor reference frame, such as
r
vqs , becomes a slightly arduous task since we must incorporate reference frame transformations in order to
eliminate all of the algebraic constraints. In particular, we consider the first two equations of (38), but we
entirely ignore the 0 sequence portion of the transformation.
    r 
iq cos(θe − δr ) − sin(θe − δr ) iq
=
id sin(θe − δr ) cos(θe − δr ) ird
    r 
vq cos(θe − δr ) − sin(θe − δr ) vq
=
vd sin(θe − δr ) cos(θe − δr ) vdr
We treat vqr and vdr as independent state variables, so we may write vq and vd explicitly as functions of vqr
and vdr . We write out these functions for generators 1 and 2.
v1q = cos(θe − δr1 )v1rq − sin(θe − δr1 )v1rd (43)
v 1d = sin(θe − δr1 )v1rq + cos(θe − δr1 )v1rd (44)
v2q = cos(θe − δr2 )v2rq − sin(θe − δr2 )v2rd (45)
v 2d = sin(θe − δr2 )v2rq + cos(θe − δr2 )v2rd (46)
We write out a similar set of equations for the terminal generator currents.
i1q = cos(θe − δr1 )ir1q − sin(θe − δr1 )ir1d (47)
i1d = sin(θe − δr1 )ir1q + cos(θe − δr1 )ir1d (48)
i 2q = cos(θe − δr2 )ir2q − sin(θe − δr2 )ir2d (49)
i2d = sin(θe − δr2 )ir2q + cos(θe − δr2 )ir2d (50)
Unfortunately, irq and ird are functions of flux state variables, so iq and id must be written in term of the flux
state variables (moving forward, we drop the subscript s which was used on terminal voltages and currents
r
in the previous section, such that vqs = vqr , for example). For a generic generator, we have the following
definitions for the terminal currents in the rotor reference frame.
0 0
!!
r 1 λds λf d λkd
id = − λds − Lad + 0 + 0
Lls Lls Llf d Llkd
0 0 !!
r 1 λqs λkq1 λkq2
iq = − λqs − Laq + 0 + 0
Lls Lls Llkq1 Llkq2

We now turn our attention to the network equations which are originally given by (40). We first substitute
in the updated expressions for voltages (43-46).

  

di1 di2


 cos(θe − δr1 )v1rq − sin(θe − δr1 )v1rd = v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + L∞ dtq + dtq
  

di di
sin(θe − δr1 )v1rq + cos(θe − δr1 )v1rd = v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q + L∞ dt1d + dt2d

Network Equations =
 cos(θ − δ )v r − sin(θ − δ )v r di2q


 e r2 2q e r2 2d = cos(θe − δr1 )v1rq − sin(θe − δr1 )v1rd + rl i2q + ωe Ll i2d + Ll dt
di
sin(θe − δr2 )v2q + cos(θe − δr2 )v2rd
r
= sin(θe − δr1 )v1rq + cos(θe − δr1 )v1rd + rl i2d − ωe Ll i2q + Ll dt2d



18
Next, we choose to write each of these expression in terms of a terminal voltage variable.
 di di2q

1q
v∞q +r∞ (i1q +i2q )+ωe L∞ (i1d +i2d )+L∞ dt + dt +sin(θe −δr1 )v1r
d
v1rq = cos(θe −δr1 )
di1 di2

v∞d +r∞ (i1d +i2d )−ωe L∞ (i1q +i2q )+L∞ d
dt + dt
d −sin(θe −δr1 )v1rq
r
v1d = cos(θe −δr ) (51)
1
di2q
cos(θe −δr1 )v1rq −sin(θe −δr1 )v1r +rl i2q +ωe Ll i2d +Ll +sin(θe −δr2 )v2r
v2rq = d
cos(θe −δr2 )
dt d

di2
sin(θe −δr1 )v1rq +cos(θe −δr1 )v1r +rl i2d −ωe Ll i2q +Ll dtd −sin(θe −δr2 )v2rq
v2rd = d
cos(θe −δr2 )

We are able to eliminate all of the voltage terms on the right hand side of each line of (51) via substitution.
 di di2q
 h  di di2
i
1q




1
v∞q + r∞ i1q + i2q + ωe L∞ i1 + i2 + L∞ + + tan(θe − δr1 ) v ∞ + r∞ i1 + i2 − ωe L∞ i1q + i2q + L∞ d + d
d d dt dt d d d dt dt
r
v1 =
q


cos(θe − δr1 ) 1 + tan2 (θe − δr1 )
(52)
 di di2
 h  di di2q
i
1q


1d d



v ∞ + r∞ i1 + i2 − ωe L∞ i1q + i2q + L∞ + − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ i1 + i2 + L∞ +
d d d dt dt d d dt dt
r
v1 =

d
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )
(53)

di2q
h di2
i
cos(θe − δr1 )v r − sin(θe − δr1 )v r + rl i2q + ωe Ll i2 + Ll + tan(θe − δr2 ) sin(θe − δr1 )v r + cos(θe − δr1 )v r + rl i2 − ωe Ll i2q + Ll d
1q 1d d dt 1q 1d d dt
r
v2 =
q


cos(θe − δr2 ) 1 + tan2 (θe − δr2 )
(54)

di2
h di2q
i
sin(θe − δr1 )v r + cos(θe − δr1 )v r + rl i2 − ωe Ll i2q + Ll d − tan(θ − δ
e r2 ) cos(θe − δr1 )v r − sin(θe − δr1 )v r + rl i2q + ωe Ll i2 + Ll
1q 1d d dt 1q 1d d dt
r
v2 =

d
cos(θe − δr2 ) 1 + tan2 (θe − δr2 )
(55)

At this point, these voltage expressions may be substituted into the equations for λ̇qs and λ̇ds in (??) for
both generators 1 and 2. In order to build the Jacobian for the system though, we must determine values for
the current derivatives in terms of flux values. This is because these current derivative themselves contain
and depend on terminal voltage variables. We consider these individually, beginning with i̇2q and i̇2d since
these are the only two derivatives which show up in (54) and (55).
d d d
i2 = cos(θe − δr2 )ir2q − sin(θe − δr2 )ir2d
dt q dt dt
= cos(θe − δr2 )i̇r2q − ir2q sin(θe − δr2 ) θ̇e − δ̇r2 − sin(θe − δr2 )i̇r2d − cos(θe − δr2 )ir2d θ̇e − δ̇r2
   
 
1
= cos(θe − δr2 )i̇r2q −
 
(λ2mq − λ2qs ) sin(θe − δr2 ) θ̇e − δ̇r2 (56)
Lls
 
1
− sin(θe − δr2 )i̇r2d − cos(θe − δr2 )
 
(λ2md − λ2ds ) θ̇e − δ̇r2 (57)
Lls
d d d
i 2d = sin(θe − δr2 )ir2q + cos(θe − δr2 )ir2d
dt dt dt
= sin(θe − δr2 )i̇r2q + cos(θe − δr2 ) θ̇e − δ̇r2 ir2q + cos(θe − δr2 )i̇r2d − sin(θe − δr2 ) θ̇e − δ̇r2 ir2d
   
 
r
  1
= sin(θe − δr2 )i̇2q + cos(θe − δr2 ) θ̇e − δ̇r2 (λ2mq − λ2qs ) (58)
Lls
 
 1
+ cos(θe − δr2 )i̇r2d − sin(θe − δr2 ) θ̇e − δ̇r2

(λ2md − λ2ds ) (59)
Lls
Next, we define terminal current derivatives in terms of flux variables.
 
1
ir2d = (λ2md − λ2ds )
Lls
 
1
ir2q = (λ2mq − λ2qs )
Lls

19
We take the derivative and expand.
0 0
!
λ̇2f d
 
Lad λ̇ Lad 1
i̇r2d = 0 + 2kd
0 + λ̇2ds −
Lls Llf d Llkd L2ls Lls
0 0 !
λ̇2kq1 λ̇2kq2
 
Laq Laq 1
i̇r2q = 0 + 0 + λ̇2qs −
Lls Llkq1 Llkq2 L2ls Lls

We now substitute in values for the flux derivatives.

 0     0   
rf d 0 0 rkd 0
0 λ2md − λ2f d + v2f d 0 λ2md − λ2kd 
Lad  L
 lf d
Llkd
i̇r2d = 0 + 0

Lls  Llf d Llkd 

  
rs Lad 1
+ ω2r λ2qs + (λ2md − λ2ds ) + v2rd −
Lls L2ls Lls
 0     0   
rkq1 0 rkq2 0
0 λ 2mq − λ2kq1 0 λ 2mq − λ2kq2
Laq  L
 lkq1
Llkq2
i̇r2q

= 0 + 0

Lls  Llkq1 Llkq2 

  
rs r Laq 1
+ −ω2r λ2ds + (λ2mq − λ2qs ) + v2q −
Lls L2ls Lls

Finally, we pull out the terminal voltage terms and we introduce two new variables, Q∗2d and Q∗2q , in order
to simplify the notation.
 0     0   
rf d 0 0 rkd 0
0 λ2md − λ2f d + v 2f d 0 λ2md − λ 2kd
Lad  L
 lf d
Llkd
i̇r2d =

0 + 0

Lls  Llf d Llkd 

    
rs Lad 1 r Lad 1
+ ω2r λ2qs + (λ2md − λ2ds ) − + v2d −
Lls L2ls Lls L2ls Lls
 
Lad 1
= Q∗2d + v2rd −
L2ls Lls
 0    r 0   
rkq1 0
kq2
0
0 λ2mq − λ2kq1 0 λ2mq − λ2kq2 
Laq  L
 lkq1
Llkq2
i̇r2q = 0 + 0

Lls  Llkq1 Llkq2 

    
rs Laq 1 r Laq 1
+ −ω2r λ2ds + (λ2mq − λ2qs ) − + v 2q −
Lls L2ls Lls L2ls Lls
 
Laq 1
= Q∗2q + v2rq −
L2ls Lls

These expressions can be placed back into (56) and (58). Performing all of these substitutions will be
extremely laborious, so we simply choose to pull out the voltage terms and leave the substitution work to

20
MATLAB.
    
d Laq 1 1
i2q = cos(θe − δr2 ) Q∗2q + v2rq
 
2 − − (λ2mq − λ 2qs ) sin(θe − δr2 ) θ̇e − δ̇r2
dt Lls Lls Lls
    
Lad 1 1
− sin(θe − δr2 ) Q∗2d + v2rd
 
2 − − cos(θ e − δ r2 ) (λ 2md − λ 2ds ) θ̇e − δ̇r2
Lls Lls Lls
    
d Laq 1  1
i2d = sin(θe − δr2 ) Q∗2q + v2rq

− + cos(θ e − δ r2 ) θ̇ e − δ̇ r2 (λ 2mq − λ 2qs )
dt L2ls Lls Lls
    
∗ r Lad 1   1
+ cos(θe − δr2 ) Q2d + v2d − − sin(θe − δr2 ) θ̇e − δ̇r2 (λ2md − λ2ds )
L2ls Lls Lls

We again pull out the voltage terms and we introduce two more variables, Q∗∗ ∗∗
2q and Q2d in order to simplify
notation.
d h i  1 
∗
(λ2mq − λ2qs ) sin(θe − δr2 ) θ̇e − δ̇r2 − sin(θe − δr2 ) Q∗2d
   
i2q = cos(θe − δr2 ) Q2q −
dt Lls
 
1  
− cos(θe − δr2 ) (λ2md − λ2ds ) θ̇e − δ̇r2
Lls
   
r Laq 1 r Lad 1
+ v2q cos(θe − δr2 ) − − v2d sin(θe − δr2 ) −
L2ls Lls L2ls Lls
   
∗∗ r Laq 1 r Lad 1
=Q2q + v2q cos(θe − δr2 ) − − v2d sin(θe − δr2 ) −
L2ls Lls L2ls Lls
 
d h i 1
i2d = sin(θe − δr2 ) Q∗2q + cos(θe − δr2 ) θ̇e − δ̇r2 (λ2mq − λ2qs ) + cos(θe − δr2 ) Q∗2d
   
dt Lls
 
  1
− sin(θe − δr2 ) θ̇e − δ̇r2 (λ2md − λ2ds )
Lls
   
r Laq 1 r Lad 1
+ v2q sin(θe − δr2 ) − + v2d cos(θe − δr2 ) −
L2ls Lls L2ls Lls
   
∗∗ r L aq 1 r Lad 1
=Q2d + v2q sin(θe − δr2 ) − + v2d cos(θe − δr2 ) −
L2ls Lls L2ls Lls

We now return to our expressions for terminal voltage in (54) and (55) and we pull out the current derivative
terms.
r r
 r r

cos(θe − δr1 )v1q − sin(θe − δr1 )v1 + rl i2q + ωe Ll i2d + tan(θe − δr2 ) sin(θe − δr1 )v1q + cos(θe − δr1 )v1 + rl i2d − ωe Ll i2q
v2rq = d d
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
h di i
Ll tan(θe − δr2 )
h di i
Ll 2q 2d
+ +
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) dt cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) dt
 
sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q − tan(θe − δr2 ) cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d
v2rd = d d
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
h di i
Ll tan(θe − δr2 )
h di i
Ll 2d 2q
+ −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) dt cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) dt

21
We then substitute in our newly derived current derivative expressions.
 
cos(θe − δr1 )v1q − sin(θe − δr1 )v1 + rl i2q + ωe Ll i2d + tan(θe − δr2 ) sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q
r r
v2rq = d d
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
    
Ll Laq 1 Lad 1
+ Q∗∗ r
2q + v2q cos(θe − δr2 ) − − v2rd sin(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls L2ls Lls
    
Ll tan(θe − δr2 ) Laq 1 Lad 1
+ Q∗∗ r
2d + v2q sin(θe − δr2 ) − + v2rd cos(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls L2ls Lls
 
sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q − tan(θe − δr2 ) cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d
v2rd = d d
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
    
Ll Laq 1 Lad 1
+ Q∗∗ r
2d + v2q sin(θe − δr2 ) − + v2rd cos(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls L2ls Lls
    
Ll tan(θe − δr2 ) Laq 1 Lad 1
− Q∗∗ r
2q + v2q cos(θe − δr2 ) − − v2rd sin(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls L2ls Lls

Once again, we pull out voltage terms.

    
L aq 1 Laq 1
Ll cos(θe − δr2 ) L2
− Lls
Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
v2rq 1 − ls
− ls 
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
 
cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d + tan(θe − δr2 ) sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q
d d
=
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
  
Ll Lad 1
+ Q∗∗ r
2q − v2d sin(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls
  
Ll tan(θe − δr2 ) Lad 1
+ Q∗∗ r
2d + v2d cos(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls
    
Lad 1 Lad 1
Ll cos(θe − δr2 ) L2
− Lls
Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
v2rd 1 − ls
− ls 
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
 
sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q − tan(θe − δr2 ) cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d
d d
=
cos(θe − δr2 ) (1 + tan2 (θe − δr2 ))
  
Ll Laq 1
+ Q∗∗ r
2d + v2q sin(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls
  
Ll tan(θe − δr2 ) Laq 1
− Q∗∗ r
2q + v2q cos(θe − δr2 ) −
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) L2ls Lls

We are able to perform on last set of substitutions and rearrangements in order to isolate the terminal

22
voltages.
 
cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d + tan(θe − δr2 ) sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q
v2rq = d
  d
 
Laq 1 Laq 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls
 
Ll Q∗∗ ∗∗
2q + tan(θe − δr2 )Q2 d
+    
Laq 1 Laq 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls
h    i
Lad 1 Lad 1
Ll − sin(θe − δr2 ) L2
− Lls
+ tan(θe − δr2 ) cos(θe − δr2 ) L2
− Lls  
+ ls
  ls
  v2rd
Laq 1 Laq 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls
 
sin(θe − δr1 )v1rq + cos(θe − δr1 )v1r + rl i2d − ωe Ll i2q − tan(θe − δr2 ) cos(θe − δr1 )v1rq − sin(θe − δr1 )v1r + rl i2q + ωe Ll i2d
v2rd = d
  d
 
Lad 1 Lad 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls
 
Ll Q∗∗ ∗∗
2 − tan(θe − δr2 )Q2q
d
+    
Lad 1 Lad 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls
h    i
Laq 1 Laq 1
Ll sin(θe − δr2 ) L2
− Lls
− tan(θe − δr2 ) cos(θe − δr2 ) L2
− Lls  
+ ls
  ls
  v2rq
Lad 1 Lad 1
cos(θe − δr2 ) (1 + tan2 (θe − δr2 )) − Ll cos(θe − δr2 ) L2
− Lls
− Ll tan(θe − δr2 ) sin(θe − δr2 ) L2
− Lls
ls ls

In order to simplify the expressions, we develop a new set of denominator coefficient variables.
   
2

Laq 1 Laq 1
D2q = cos(θe − δr2 ) 1 + tan (θe − δr2 ) − Ll cos(θe − δr2 ) − − Ll tan(θe − δr2 ) sin(θe − δr2 ) −
L2 Lls L2 Lls
 ls   ls 
2

Lad 1 Lad 1
D2d = cos(θe − δr2 ) 1 + tan (θe − δr2 ) − Ll cos(θe − δr2 ) − − Ll tan(θe − δr2 ) sin(θe − δr2 ) −
L2ls Lls L2ls Lls
We develop numerator coefficients as well, and we pull out all terminal voltage terms.
 
1 rl i2q + ωe Ll i2d + tan(θe − δr2 ) rl i2d − ωe Ll i2q
T2q =
D2q
 
cos(θ e − δ r1 ) + tan(θ e − δr2 ) sin(θe − δr1 )
Tvq v
2q 1q
r
= v1rq
D2q
 
tan(θe − δr2 ) cos(θe − δr1 ) − sin(θe − δr1 )
Tvd v
2q 1d
r
= v1rd
D2q
h i
Ll Q∗∗ 2q + tan(θe − δr2 )Q2d
∗∗

T22q =
D2q
h    i
Ll − sin(θe − δr2 ) LLad 2 − L
1
ls
+ tan(θe − δr2 ) cos(θe − δr2 ) LLad 2 − L
1
ls
T32q = ls ls

D2q
 
1 rl i2d − ωe Ll i2q − tan(θe − δr2 ) rl i2q + ωe Ll i2d
T2d =
D2d
 
vq r sin(θ e − δ r1 ) − tan(θ e − δr2 ) cos(θe − δr1 )
T2d v1q = v1rq
D2q
 
cos(θe − δr1 ) + tan(θe − δr2 ) sin(θe − δr1 )
Tvd v
2d 1d
r
= v1rd
D2q
h i
Ll Q∗∗ 2d − tan(θe − δr2 )Q2q
∗∗

T22d =
D2d
h    i
L 1 Laq 1
Ll sin(θe − δr2 ) Laq 2 − L
ls
− tan(θ e − δ r 2 ) cos(θ e − δr 2 ) 2
Lls
− Lls
T32d = ls

D2d

23
We substitute in these coefficients.
v2rq =T12q + Tvq r vd r
 2
3 r
   
2q v1q + T2q v1d + T2q + T2q v2d

v2rd =T12d + Tvq

     
2d
v1rq + Tvd r 2 3 r
2d v1d + T2d + T2d v2q

Next, we solve for generator 2 terminal voltages in terms of generator 1 terminal voltages.
vq vq
 r   r   1  r   r   
r
v2q =T2q + T2q v1q + T2q v1d + T2q + T2q T2d + T2d v1q + T2d v1d + T2d + T2d v2rq
1 vd 2 3 vd 2 3

v2rd =T12d + Tvq vq

          
2d
v1rq + Tvd r 2 3 1 r vd r 2 3 r
2d v1d + T2d + T2d T2q + T2q v1q + T2q v1d + T2q + T2q v2d

T12q + Tvq + T22q + T32q T12d + Tvq

r
 vd r
        
2q v1q + T2q v1 2d
v1rq + Tvd r
2d v1 + T22d
v2rq = d

d

1 − T32q T32d
 
T12q + T22q + T32q T12d + T22d Tvq 3 vq
2q + T2q T2d
  Tvd + T3 Tvd  
=
+
v1rq + 2q 3 2q 3 2d
v1rd
1 − T32q T32d 1 − T32q T32d 1 − T2q T2d
vq
Tvq
         
T12d + T2d v1rq vd
+ T2d v1r + T22d + T32d T12q + 2q v1rq + Tvd r
2q v1 + T22q
v2rd = d

d

1 − T32d T32q
 
T12d + T22d + T32d T12q + T22q Tvq + T32d Tvq
2q
  Tvd + T3 Tvd  
=
+
2d

v1rq + 2d 3 2d 3 2q
v1rd
1− T32d T32q 1− T32d T32q 1 − T2d T2q

We now move on to consider the terminal voltages at generator 1. We state with the current derivatives in
the synchronous reference frame.
d d d
i1q = cos(θe − δr1 )ir1q − sin(θe − δr1 )ir1d
dt dt dt
= cos(θe − δr1 )i̇r1q − ir1q sin(θe − δr1 ) θ̇e − δ̇r1 − sin(θe − δr1 )i̇r1d − cos(θe − δr1 )ir1d θ̇e − δ̇r1
   
 
r 1  
= cos(θe − δr1 )i̇1q − (λ1mq − λ1qs ) sin(θe − δr1 ) θ̇e − δ̇r1 (60)
Lls
 
1
− sin(θe − δr1 )i̇r1d − cos(θe − δr1 )
 
(λ1md − λ1ds ) θ̇e − δ̇r1 (61)
Lls
d d d
i1 = sin(θe − δr1 )ir1q + cos(θe − δr1 )ir1d
dt d dt dt
= sin(θe − δr1 )i̇r1q + cos(θe − δr1 ) θ̇e − δ̇r1 ir1q + cos(θe − δr1 )i̇r1d − sin(θe − δr1 ) θ̇e − δ̇r1 ir1d
   
 
 1
= sin(θe − δr1 )i̇r1q + cos(θe − δr1 ) θ̇e − δ̇r1

(λ1mq − λ1qs ) (62)
Lls
 
 1
+ cos(θe − δr1 )i̇r1d − sin(θe − δr1 ) θ̇e − δ̇r1

(λ1md − λ1ds ) (63)
Lls
Next, we define terminal current derivatives in terms of flux variables.
 
r 1
i 1d = (λ1md − λ1ds )
Lls
 
r 1
i1q = (λ1mq − λ1qs )
Lls
We take the derivative of each expression and expand.
0 0
!
Lad λ̇1f d
 
r λ̇1kd Lad 1
i̇1d = 0 + 0 + λ̇1ds −
Lls Llf d Llkd L2ls Lls
0 0 !
Laq λ̇1kq1 λ̇1kq2
 
Laq 1
i̇r1q = 0 + 0 + λ̇ 1qs −
Lls Llkq1 Llkq2 L2ls Lls

24
We now substitute in values for the flux derivatives.
 0     0   
rf d 0 0 rkd 0
0 λ 1md − λ 1f d + v 1f d 0 λ 1md − λ1kd
  
r Lad  Llf d Llkd  rs r Lad 1
i̇1d = +  + ω1r λ1qs + Lls (λ1md − λ1ds ) + v1d −
 0 0

Lls  Llf d Llkd L2ls Lls

 0  0
    
rkq1 0 rkq2 0
0 λ1mq − λ1kq1 0 λ1mq − λ1kq2   
Laq  Llkq1 Llkq2
 + −ω1r λ1ds + rs (λ1mq − λ1qs ) + v1r Laq 1
i̇r1q

=  0 + 0 −
Lls  Llkq1 Llkq2  Lls q
L2ls Lls

Finally, we pull out the terminal voltage terms.

 0     0   
rf d 0 0 rkd 0
0 λ1md − λ1f d + v1f d 0 λ1md − λ1kd 
Lad  L
 lf d
Llkd
i̇r1d = 0 + 0

Lls  Llf d Llkd 

    
rs Lad 1 r Lad 1
+ ω1r λ1qs + (λ1md − λ1ds ) − + v1d −
Lls L2ls Lls L2ls Lls
 
L ad 1
= Q∗1d + v1rd −
L2ls Lls
 0    r 0   
rkq1 0
kq2
0
0 λ 1mq − λ1kq1 0 λ 1mq − λ 1kq2
Laq  L
 lkq1
Llkq2
i̇r1q

= 0 + 0

Lls  Llkq1 Llkq2 

    
rs Laq 1 r Laq 1
+ −ω1r λ1ds + (λ1mq − λ1qs ) − + v 1q −
Lls L2ls Lls L2ls Lls
 
Laq 1
= Q∗1q + v1rq −
L2ls Lls

Finally, these expressions can be placed back into (60) and (62). Performing all of these substitutions will
be extremely laborious, so we simply choose to pull out the voltage terms and leave the substitution work
to MATLAB.
    
d ∗ r Laq 1 1  
i1q = cos(θe − δr1 ) Q1q + v1q 2 − − (λ1mq − λ1qs ) sin(θe − δr1 ) θ̇e − δ̇r1
dt Lls Lls Lls
    
∗ r L ad 1 1  
− sin(θe − δr1 ) Q1d + v1d 2 − − cos(θe − δr1 ) (λ1md − λ1ds ) θ̇e − δ̇r1
Lls Lls Lls
    
d ∗ r L aq 1   1
i1 = sin(θe − δr1 ) Q1q + v1q − + cos(θe − δr1 ) θ̇e − δ̇r1 (λ1mq − λ1qs )
dt d L2ls Lls Lls
    
Lad 1  1
+ cos(θe − δr1 ) Q∗1d + v1rd

− − sin(θ e − δ r ) θ̇ e − δ̇ r (λ1md − λ1ds )
L2ls Lls 1 1
Lls

25
We again pull out the voltage terms.

d h i  1 
∗
 
i1q = cos(θe − δr1 ) Q1q − (λ1mq − λ1qs ) sin(θe − δr1 ) θ̇e − δ̇r1
dt Lls
 
 ∗  1  
− sin(θe − δr1 ) Q1d − cos(θe − δr1 ) (λ1md − λ1ds ) θ̇e − δ̇r1
Lls
   
r L aq 1 r Lad 1
+ v1q cos(θe − δr1 ) − − v1d sin(θe − δr1 ) −
L2ls Lls L2ls Lls
   
∗∗ r Laq 1 r Lad 1
=Q1q + v1q cos(θe − δr1 ) − − v1d sin(θe − δr1 ) −
L2ls Lls L2ls Lls
 
d h i 1
i1 = sin(θe − δr1 ) Q∗1q + cos(θe − δr1 ) θ̇e − δ̇r1
 
(λ1mq − λ1qs )
dt d Lls
 
 ∗    1
+ cos(θe − δr1 ) Q1d − sin(θe − δr1 ) θ̇e − δ̇r1 (λ1md − λ1ds )
Lls
   
r Laq 1 r Lad 1
+ v1q sin(θe − δr1 ) − + v1d cos(θe − δr1 ) −
L2ls Lls L2ls Lls
   
∗∗ r Laq 1 r L ad 1
=Q1d + v1q sin(θe − δr1 ) − + v1d cos(θe − δr1 ) −
L2ls Lls L2ls Lls

We may substitute these current derivative values (and the previously derived ones) into the network voltage
equations.



r v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + tan(θe − δr1 ) v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q
v 1q = +
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


   
L∞ di1q L∞ di2q
+
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) cos(θe − δr1 ) 1 + tan2 (θe − δr1 )



dt dt
   
tan(θe − δr1 )L∞ di1d tan(θe − δr1 )L∞ di2d
+ +
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) cos(θe − δr1 ) 1 + tan2 (θe − δr1 )



dt dt



v∞ + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d )
v1rd = d +
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


   
L∞ di1d L∞ di2d
+
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) cos(θe − δr1 ) 1 + tan2 (θe − δr1 )



dt dt
   
tan(θe − δr1 )L∞ di1q tan(θe − δr1 )L∞ di2q
− −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) cos(θe − δr1 ) 1 + tan2 (θe − δr1 )



dt dt

26



v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + tan(θe − δr1 ) v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q
v1rq =
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


    
L∞ ∗∗ r Laq 1 r Lad 1
+
Q1q + v1q cos(θe − δr1 ) 2 − L − v1d sin(θe − δr1 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) L ls ls L2ls Lls
    
L∞ Laq 1 Lad 1
+
Q∗∗ r
2q + v2q cos(θe − δr2 ) 2 − − v2rd sin(θe − δr2 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) Lls Lls L2ls Lls
    
tan(θe − δr1 )L∞ Laq 1 Lad 1
+
Q∗∗ r
1d + v1q sin(θe − δr1 ) − + v1rd cos(θe − δr1 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls Lls L2ls Lls
    
tan(θe − δr1 )L∞ Laq 1 Lad 1
+
Q∗∗ r
2d + v2q sin(θe − δr2 ) 2 − L + v2rd cos(θe − δr2 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) L ls ls L2ls Lls



v∞ + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d )
v1rd = d
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


    
L∞ ∗∗ r Laq 1 r Lad 1
+
Q1d + v1q sin(θe − δr1 ) 2 − L + v1d cos(θe − δr1 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) L ls ls L2ls Lls
    
L∞ Laq 1 Lad 1
+
Q∗∗ r
2d + v2q sin(θe − δr2 ) 2 − + v2rd cos(θe − δr2 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) Lls Lls L2ls Lls
    
tan(θe − δr1 )L∞ Laq 1 Lad 1
−
Q∗∗ r
1q + v1q cos(θe − δr1 ) − − v1rd sin(θe − δr1 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls Lls L2ls Lls
    
tan(θe − δr1 )L∞ Laq 1 Lad 1
−
Q∗∗ r
2q + v2q cos(θe − δr2 ) 2 − L − v2rd sin(θe − δr2 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) L ls ls L2ls Lls

Now we pull out the generator 1 terminal voltage terms.

27
     
Laq 1 Laq
L∞ cos(θe − δr1 ) L2ls
− Lls tan(θe − δr1 )L∞ sin(θe − δr1 ) − L1ls
L2ls
v1rq 1 − 2

− 
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


cos(θe − δr1 ) 1 + tan (θe − δr1 )



v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + tan(θe − δr1 ) v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q
=
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


  
L∞ ∗∗ r Lad 1
+ Q 1q − v 1d sin(θ e − δ r1 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls


Lls
    
L∞ ∗∗ r Laq 1 r Lad 1
+ Q2q + v2q cos(θe − δr2 ) − − v2d sin(θe − δr2 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls L2ls


Lls Lls
  
tan(θe − δr1 )L∞ Lad 1
+ 2

Q∗∗ r
1d + v1d cos(θe − δr1 ) 2 −
cos(θe − δr1 ) 1 + tan (θe − δr1 ) Lls Lls
    
tan(θe − δr1 )L∞ ∗∗ r Laq 1 r Lad 1
+ Q 2d + v 2q sin(θ e − δ r2 ) − + v 2d cos(θ e − δ r2 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls L2ls


Lls Lls
    
L∞ cos(θe − δr1 ) LLad
2 − L
1
ls
sin(θe − δr1 ) tan(θe − δr1 )L∞ LLad 2 − L
1
ls
v1rd 1 − ls
− ls 
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) cos(θe − δr1 ) 1 + tan2 (θe − δr1 )



v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d )
=
cos(θe − δr1 ) 1 + tan2 (θe − δr1 )


  
L∞ ∗∗ r Laq 1
+ Q1d + v1q sin(θe − δr1 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls


Lls
    
L∞ ∗∗ r Laq 1 r Lad 1
+
Q2d + v2q sin(θe − δr2 ) 2 − + v2d cos(θe − δr2 ) −
2
cos(θe − δr1 ) 1 + tan (θe − δr1 ) Lls Lls L2ls Lls
  
tan(θe − δr1 )L∞ Laq 1
−
Q∗∗ r
1q + v1q cos(θe − δr1 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls Lls
    
tan(θe − δr1 )L∞ ∗∗ r Laq 1 r Lad 1
− Q2q + v2q cos(θe − δr2 ) − − v2d sin(θe − δr2 ) −
cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) L2ls L2ls


Lls Lls

28
     
Laq 1 Laq 1
v1rqcos(θe − δr1 ) 1 + tan2 (θe − δr1 ) − L∞ cos(θe − δr1 )


− − tan(θ e − δ r )L ∞ sin(θ e − δ r ) −
L2ls Lls 1 1
L2ls Lls



=v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + tan(θe − δr1 ) v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q
  
Lad 1
+ L∞ Q∗∗ 1q − v r
1d sin(θ e − δ r1 ) −
L2 Lls
  ls   
Laq 1 Lad 1
+ L∞ Q∗∗ 2q + v r
2q cos(θ e − δ r2 ) − − v r
2d sin(θ e − δ r2 ) −
L2ls Lls L2ls Lls
  
Lad 1
+ tan(θe − δr1 )L∞ Q∗∗ r
1d + v1d cos(θe − δr1 ) −
L2ls Lls
    
Laq 1 Lad 1
+ tan(θe − δr1 )L∞ Q∗∗ 2d + v r
2q sin(θ e − δ r2 ) − + v r
2d cos(θ e − δr2 ) −
L2ls Lls L2ls Lls
    
Lad 1 Lad 1
v1rd cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) − L∞ cos(θe − δr1 )


− − sin(θ e − δ r1 ) tan(θ e − δ r1 )L ∞ −
L2ls Lls L2ls Lls



=v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d )
  
∗∗ r Laq 1
+ L∞ Q1d + v1q sin(θe − δr1 ) −
L2 Lls
  ls   
∗∗ r Laq 1 r Lad 1
+ L∞ Q2d + v2q sin(θe − δr2 ) − + v2d cos(θe − δr2 ) −
L2ls Lls L2ls Lls
  
Laq 1
− tan(θe − δr1 )L∞ Q∗∗ r
1q + v1q cos(θe − δr1 ) −
L2ls Lls
    
∗∗ r Laq 1 r Lad 1
− tan(θe − δr1 )L∞ Q2q + v2q cos(θe − δr2 ) − − v2d sin(θe − δr2 ) −
L2ls Lls L2ls Lls

We now choose to define two denominator values.

   
2

Laq 1 Laq 1
D1q = cos(θe − δr1 ) 1 + tan (θe − δr1 ) − L∞ cos(θe − δr1 ) − − tan(θe − δr1 )L∞ sin(θe − δr1 ) −
L2ls Lls L2ls Lls
   
Lad 1 Lad 1
= cos(θe − δr1 ) 1 + tan2 (θe − δr1 ) − L∞ cos(θe − δr1 )


D1d − − sin(θ e − δ r )L∞ tan(θ e − δ r ) −
L2ls Lls 1 1
L2ls Lls

29
We use these denominator definitions to simplify the expressions.
"

 #
r v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d ) + tan(θe − δr1 ) v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q
v1q =
D1q
" #
L∞ Q1q + L∞ Q2q + tan(θe − δr1 )L∞ Q1d + tan(θe − δr1 )L∞ Q∗∗
∗∗ ∗∗ ∗∗
2d
+
D1q
    
tan(θe − δr1 )L∞ cos(θe − δr1 ) LLad2 − L
1
ls
− L ∞ sin(θ e − δ r 1
) Lad
L 2 − L
1
ls
 v1r
ls ls
 
+
D1q d

    
L 1 Laq 1
L∞ cos(θe − δr2 ) Laq2 − L
ls
+ tan(θ e − δ r 1
)L ∞ sin(θ e − δ r 2
) 2
Lls
− Lls
h i
+ ls  v2r
D1q q

    
tan(θe − δr1 )L∞ cos(θe − δr2 ) LLad2 − L
1
ls
− L ∞ sin(θ e − δ r 2
) Lad
2
Lls
− 1
Lls
 v2r
ls
 
+
D1q d

"

#
r v∞d + r∞ (i1d + i2d ) − ωe L∞ i1q + i2q − tan(θe − δr1 ) v∞q + r∞ i1q + i2q + ωe L∞ (i1d + i2d )
v 1d =
D1d
" #
∗∗ ∗∗ ∗∗
+L∞ Q1d + L∞ Q2d − Q1q tan(θe − δr1 )L∞ − tan(θe − δr1 )L∞ Q∗∗ 2q
+
D1d
    
L 1 L 1
L∞ sin(θe − δr1 ) Laq2 − L
ls
− tan(θe − δr1 )L∞ cos(θe − δr1 ) Laq 2 − L
ls
h i
+ ls ls  v1r
D1d q

    
L 1 Laq 1
L∞ sin(θe − δr2 ) Laq2 − L
ls
− tan(θ e − δ r 1
)L ∞ cos(θ e − δ r 2
) 2
Lls
− L ls
h i
+ ls  v2r
D1d q

    
tan(θe − δr1 )L∞ sin(θe − δr2 ) LLad2 − L
1
ls
+ L ∞ cos(θ e − δ r 2 ) Lad
2
Lls
− 1
Lls
 v2r
ls
 
+
D1d d

We assign variable names to the coefficients attached to the voltage variables (we employ alternative alpha-
betic variables for simplified notation).
h i
v1rq = Tq1 + Tq2 v1rd + Tq3 v2rq + Tq4 v2rd
   
h i
= a + b v1rd + c v2rq + d v2rd
   
h i h i
v1rd = Td1 + Td2 v1rq + Td3 v2rq + Td4 v2rd
 
h i h i
= e + f v1rq + g v2rq + h v2rd
 

At this point, we import the expressions for v2rq and v2rd (we again employ alternative alphabetic variables
for simplified notation).
 
r
T12q + T22q + T32q T12d + T22d Tvq 3 vq
2q + T2q T2d
 r  Tvd 3 vd 
2q + T2q T2d

v 2q =
3 3

+
3 3

v1q +
3 3

v1rd
1 − T2q T2d 1 − T2q T2d 1 − T2q T2d
   
=i + j v1rq +k v1rd
 
T12d + T22d + T32d T12q + T22q Tvq + T32d Tvq
2q
  Tvd + T3 Tvd  
v2rd =
+
2d

v1rq + 2d 3 2d 3 2q
v1rd
1 − T32d T32q 1 − T32d T32q 1 − T2d T2q
   
=l + m v1rq + n v1rd

30
With these four equations, we may solve for the terminal voltage values purely in terms of currents and other
state variables.
h h i  i h h i  i
v1rq = a + b v1rd + c i + j v1rq + k v1rd + d l + m v1rq + n v1rd
 
h i h h i  i h h i  i
v1rd = e + f v1rq + g i + j v1rq + k v1rd + h l + m v1rq + n v1rd

We group like terms to solve the expressions.

h i
v1rq = a + ci + dl + (cj + dm) v1rq + (b + ck + dn) v1rd
 
h i
v1rd = e + gi + hl + (f + gj + hm) v1rq + (gk + hn) v1rd
 

⇓
 
a + ci + dl + (b + ck + dn) v1rd
v1rq =
(1 − (cj + dm))
Now we are able to perform substitution.
"  #
r a + ci + dl + (b + ck + dn) v1rd
+ (gk + hn) v1rd
 
v1d = e + gi + hl + (f + gj + hm)
(1 − (cj + dm))
     
a + ci + dl b + ck + dn
+ (gk + hn) v1rd
 
= e + gi + hl + (f + gj + hm) + (f + gj + hm)
1 − cj − dm 1 − cj − dm
      
r b + ck + dn a + ci + dl
v1d 1 − (f + gj + hm) + (gk + hn) = e + gi + hl + (f + gj + hm)
1 − cj − dm 1 − cj − dm
 
a+ci+dl
e + gi + hl + (f + gj + hm) 1−cj−dm
v1rd =  h   i
b+ck+dn
1 − (f + gj + hm) 1−cj−dm + (gk + hn)

We may now solve for v1rq .

 a+ci+dl

e+gi+hl+(f +gj+hm)( 1−cj−dm )
a + ci + dl + (b + ck + dn) b+ck+dn
(1−[(f +gj+hm)( 1−cj−dm )+(gk+hn)])
v1rq =
(1 − cj − dm)

Generator 2 terminal voltages follow similarly.

" #
a + ci + dl + (b + ck + dn) e + gi + hl
 
v2rq =i + j  a+ci+dl+(b+ck+dn) 
(1 − (cj + dm))
 

1 − (f + gj + hm) (1−(cj+dm))
+ (gk + hn)
" #
e + gi + hl
+k   a+ci+dl+(b+ck+dn)  

1 − (f + gj + hm) (1−(cj+dm))
+ (gk + hn)
" #
a + ci + dl + (b + ck + dn) e + gi + hl
 
v2rd =l + m  a+ci+dl+(b+ck+dn) 
(1 − (cj + dm))
 

1 − (f + gj + hm) (1−(cj+dm))
+ (gk + hn)
" #
e + gi + hl
+n   a+ci+dl+(b+ck+dn)  

1 − (f + gj + hm) (1−(cj+dm))
+ (gk + hn)

31
Linearization and Eigenvalue Analysis
Assume we have the following set of ODEs.

ẋ = f (x, y)
ẏ = g(x, y)

We perform Taylor Series expansion about an equilibrium point (or a steady state power flow solution in
our context) x0 , y0 .
   
d d
ẋ ≈ f (x0 , y0 ) + f (x, y) ∆x + f (x, y) ∆y
dx x0 ,y0 dy x0 y0
   
d d
ẏ ≈ g(x0 , y0 ) + g(x, y) ∆x + g(x, y) ∆y
dx x0 ,y0 dy x0 ,y0

We subtract the operating point from both sides in order to determine how the derivatives change, and we
build the Jacobian matrix.
  " d d
#  
∆ẋ dx f (x, y) dy f (x, y) ∆x
≈ d d
∆ẏ dx g(x, y) dy g(x, y)
∆y
x0 ,y0

We perform a similar process with our multi-generator system. For our given set of system parameters, we
build the system Jacobian and we find the 16 eigenvalues (each machine requires 8 differential equations,
and the network differential equations have been written in terms of the flux variables of the generator).
Only two of these eigenvalues have oscillatory (imaginary) components.

λ9,10 = −0.0025 ± j0.0379

λ12,13 = −0.0003 ± j0.0162

These oscillatory modes correspond to 0.6032 Hz and 0.2578 Hz, respectively. Of course, λ12,13 is the more
dominant mode of the two since is has a larger real part (thus a smaller damping component). We next
perform the following experiment. We put a oscillation on the shaft of generator 1 such that the mechanical
power input oscillates at frequency fosc at magnitude ατm0 where α = 0.1 (10% of the nominal torque).

τm = τm0 + α · τm0 · sin (2πfosc t)

In the first plot (Figure 11), we trace the standard deviation of four variables: δ1 , δ2 , ω1 , and ω2 and we
compare their differences. This effectively represents the magnitude of the phase oscillations.

32
 0.7 0.05
</1
</2
0.6 <!1 0
<!2

0.5 -0.05
Standard Deviation

Standard Deviation
0.4 -0.1

0.3 -0.15

0.2 -0.2

0.1 -0.25 </1 ! </2

<!1 ! <!2
zero
0 -0.3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Forced Oscillation Frequncy (Hz) Forced Oscillation Frequncy (Hz)

Figure 11: The left panel shows the standard deviation (a measure of magnitude) of variablesδ1 , δ2 , ω1 , and
ω2 , and the right panel shows the difference in standard deviations (the difference in magnitudes) of the
phase angles and the frequencies of the generators. There are noticeable peaks at both oscillatory modes of
0.2578 Hz and 0.6032 Hz.

In the next plot, we trace the phase of the phase oscillations. When both angles are oscillating in phase, the
phase shift is close to 0, and when both angles are out of phase, the phase shift is close to π.

3.5
Phase Shift between Voltage Phases

2.5

1.5

0.5 ?/1 ! ?/2

:
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Forced Oscillation Frequncy (Hz)

Figure 12: The phase shift between the voltage phase angles is shown. Between modes, the phase shift flops
from 0 to π.

Linearization of a Single Machine

When oscillatory mechanical torque (or power) is applied to a synchronous machine, the integral of the
machine’s power output may tend positive or negative. We wish to determine the conditions for when this
power integral will tend positive or negative. To do so, we fix the three phase terminal voltage of the machine
and study its output. In doing so, we define the three phase voltages accordingly, where Vm is the peak
voltage magnitude of a given phase voltage.

33
 va (t) = Vm cos (ωs t)
 
2π
vb (t) = Vm cos ωs t −
3
 
2π
vc (t) = Vm cos ωs t +
3

To determine the generator’s qd axis voltages, we apply park’s transformation.

     
2 2π 2π
vqs = cos (δr ) va (t) + cos δr − vb (t) + cos δr + vc (t)
3 3 3
        
2 2π 2π 2π 2π
= cos (δr ) Vm cos (ωs t) + cos δr − Vm cos ωs t − + cos δr + Vm cos ωs t +
3 3 3 3 3
        
2 2π 2π 2π 2π
= Vm cos (δr ) cos (ωs t) + cos δr − cos ωs t − + cos δr + cos ωs t +
3 3 3 3 3
     
2 2π 2π
vds = sin (δr ) va (t) + sin δr − vb (t) + sin δr + vc (t)
3 3 3
        
2 2π 2π 2π 2π
= sin (δr ) Vm cos (ωs t) + sin δr − Vm cos ωs t − + sin δr + Vm cos ωs t +
3 3 3 3 3
        
2 2π 2π 2π 2π
= Vm sin (δr ) cos (ωs t) + sin δr − cos ωs t − + sin δr + cos ωs t +
3 3 3 3 3

In order to further simplify these identities, we may incorporate the following two identities.

1
cos(x) cos(y) = [cos (x − y) + cos (x + y)]
2
1
sin(x) cos(y) = [sin (x + y) + sin (x − y)]
2
We apply the first identity to the expression for vqs .

        
2 2π 2π 2π 2π
vqs = Vm cos (δr ) cos (ωs t) + cos δr − cos ωs t − + cos δr + cos ωs t +
3 3 3 3 3
    
2 4π 4π
= Vm cos (δr − ωs t) + cos (δr + ωs t) + cos (δr − ωs t) + cos ωs t + δr − + cos (δr − ωs t) + cos ωs t + δr +
6 3 3
= Vm cos (δr − ωs t)


= Vm cos δ r

We are able to make this last substitution since δ r ≡ δr − ωs t. We now perform similar analysis on the
expression for vds .

        
2 2π 2π 2π 2π
vds = Vm sin (δr ) cos (ωs t) + sin δr − cos ωs t − + sin δr + cos ωs t +
3 3 3 3 3
     
2 4π 4π
= Vm sin (δr + ωs t) + sin (δr − ωs t) + sin δr + ωs t − + sin (δr − ωs t) + sin δr + ωs t + + sin (δr − ωs t)
6 3 3
= Vm sin (δr − ωs t)


= Vm sin δ r

34
Using these simplified expressions for vqs and vds , we may state the differential equations which model the
system.
rs



 λ̇ qs = −ω r λ ds + L ls
(λ mq − λ qs ) + V m cos δ r
rs


 λ̇ds = ωr λ0 qs + Lls (λmd − λds ) + Vm sin δ r



  
 λ̇0 = r0kq1 λmq − λ0




 kq1 Llkq1 kq1
0

  
 λ̇kq2 = r0kq2 λmq − λ0kq2
0


Llkq2
0   (64)
0 r fd r 0
r 0
r
λ̇f d = L0 λmd − λf d + vf d




 lf d
 0  
0 r 0

λ̇kd = L0kd λmd − λkd




 lkd
ω̇ = τmJ−τe




 ˙

δ r = ωr − ωs
There are three nonlinear terms in this set of equations: f1 = ωr λds , f2 = ωr λqs , and f3 = τe . To linearize
this system, we apply the first order Taylor Series. We start with f1 and f2 . The 0 superscript indicates a
variable evaluated at its steady state value.

∂f1 ∂f1
f1 ≈ f1 + 4ωr + 4λds
ωr0 ,λ0ds ∂ωr ωr0 ,λ0ds ∂λds ωr0 ,λ0ds

≈ ωr0 λ0ds + λ0ds ωr − ωr0 + ωr0 λds − λ0ds

≈ λ0ds ωr + ωr0 λds − ωr0 λ0ds

∂f2 ∂f2
f2 ≈ f2 + 4ωr + 4λqs
ωr0 ,λ0qs ∂ωr ωr0 ,λ0qs ∂λds ωr0 ,λ0qs

≈ ωr0 λ0qs + λ0qs ωr − ωr0 + ωr0 λqs − λ0qs

≈ λ0qs ωr + ωr0 λqs − ωr0 λ0qs

The electrical torque expression can be expanded and then linearized

f3 = τe
3
= (λds iqs − λqs ids )
2 
3 λqs λds + λds λmq − λqs λmd − λds λqs
=
2 Lls
0 0 0
    0
 
λqs λkq1 λkq2 λf d λkd
λqs λds + Laq λds Lls + L0 + L0 − Lad λqs λLds + 0 + 0 − λds λqs 
3 lkq1 lkq2
ls Llf d Llkd
=  
2 Lls 

In this expression, there are 6 flux variables, so the Taylor Series Linearization must include derivatives.
Each of these derivatives are evaluated at a steady state equilibrium point.
∂f3 ∂f3 0 ∂f3 0 ∂f3 ∂f3 0 ∂f3 0
f3 ≈ f30 + 4λqs + 4λkq1 + 4λkq2 + 4λds + 4λf d + 4λkd
∂λqs ∂λds ∂ωr ∂λds ∂ωr ∂λds
We may now insert the approximations for f1 , f2 , and f3 into the set of differential equations in (64). Once
done, we simulate the full model and the linearized model in order to determine its effectiveness.

35
Linearization of a Single Machine: Damper Windings Removed
We wish to remove the 3 damper windings present in the generator and then perform a linearization. In
doing so, we end up with the following sets of magnetizing flux variables.

 
λqs
λmq = Laq (65)
Lls
0 !
λds λf d
λmd = Lad + 0 (66)
Lls Llf d

Our initial set of differential equations has the following form, where θv is the angle of the generator bus
voltage. When tied to an infinite bus, θv = 0 for all time.

λ̇qs = − (ωr + ωs ) λds + Lrls

(λmq − λqs ) + Vm cos δ r − θv
 s

λ̇ds = (ωr + ωs ) λqs + Lrls


s
(λmd − λds ) + Vm sin δ r − θv




 0  
0 rf d 0 0
λ̇f d = L0 λmd − λf d + vfrd (67)
 lf d
 ω̇ = τm −τe


 J
 ˙

δ =ω −ω
r r s



First, we assume that δ r − θv  1, so cos δ r − θv ≈ 1 and sin δ r − θv ≈ δ r − θv . Next, we apply a Taylor
Series Expansion to each element. We begin by writing the electrical torque in terms of flux state variables.

3
τe = (λds iqs − λqs ids )
2 
3 λds λmq − λqs λmd
=
2 Lls
  !!
3 Laq − Lad 0 Lad
= λqs λds − λqs λf d 0
2 L2ls Lls Llf d

Now we apply Taylor Series.

 
rs Laq rs
4λ̇qs ≈ − ω 0r + ωs 4λds − λ0ds 4ω r +
   
− 4λqs
L2ls Lls
  " #
 0   0  rs Lad rs rs Lad 0
4λ̇ds ≈ ω r + ωs 4λqs + λqs 4ω r + 2 − 4λds + 0 4λf d + [Vm ] 4δ r
Lls Lls Lls Llf d
" 0 # " 0 0 #
0 rf d Lad rf d Lad rf d 0
4λ̇f d ≈ 0 4λds + 0 0 − 0 4λf d
Llf d Lls Llf d Llf d Llf d
""   !#    " !# #
3 0 Laq − Lad 0
0 Lad 0 Laq − Lad Lad 0
4ω̇ r ≈ − λds − λf d 0 4λqs + λqs 4λds − λ0qs 0 4λf d
J2 L2ls Lls Llf d L2ls Lls Llf d

4δ˙ r ≈ 4ω r

We can transform these expressions into the form ẋ = Ax + Bu, where α is a scaling factor for the oscillating
input torque.

36
  h i 
rs Laq rs
 0   0 
L2ls
− Lls −ω r − ωs 0 −λds 0
     
4λ̇qs 4λq
h i
rs Lad rs rs Lad
  0   0  
 ω r + ωs L 2 − Lls 0 λqs [Vm ] 
 4λ̇ds  
 ls0
Lls Llf d  4λd
0 0
 0
      0

 4λ̇f d =
 
0
rf d Lad
0
rf d Lad
0 0 −
rf d
0 0 0

 4λf
L Lls Llf d Llf d Llf d
4ω̇ r 4ω
   lf d

         
4δ˙ r
  h  i
3 00
 3 0 Lad −Laq Lad −Laq
Lad Lad 4δ
3 0 3 0

 J2 λds + J2 λf d L L0 J2 λqs J2 λqs 0 0
 0

L2ls ls lf d L2ls Lls Llf d 
0 0 0 1 0

We may determine the linearized electrical power output.

3
Pe = (vqs iqs + vds ids )
2 !!
 
3
Laq λqs
Lad Lad 0 λds
= Vm cos δ r λqs − + Vm sin δ r λds + λ −
Lls Llf d f d Lls
0
2 L2ls Lls L2ls
  !!
3 Laq λqs Lad Lad 0 λds
= Vm λqs − + Vm δ r λds + λ −
Lls Llf d f d Lls
0
2 L2ls Lls L2ls

Appendix I - Terminal Active Power

The instantaneous power at the terminals of a generator may be written in terms of the qd sequence variables
in the synchronous reference frame through an algebraic relationship.
p1 =v1a i1a + v1b i1b + v1c i1c
  
= cos (ωs t) vqs1 + sin (ωs t) vds1 cos (ωs t) isq1 + sin (ωs t) isd1 +
2π 2π 2π s 2π s
h     ih     i
cos ωs t − vqs1 + sin ωs t − vds1 cos ωs t − iq1 + sin ωs t − id1 +
3 3 3 3
2π 2π 2π s 2π s
h     ih     i
cos ωs t + vqs1 + sin ωs t + vds1 cos ωs t + iq1 + sin ωs t + id1
3 3 3 3
 
= cos2 (ωs t) vqs1 isq1 + sin2 (ωs t) vds1 isd1 + sin (ωs t) vds1 cos (ωs t) isq1 + cos (ωs t) vqs1 sin (ωs t) isd1 +
2π 2π 2π 2π s 2π 2π s
h             i
cos2 ωs t − vqs1 isq1 + sin2 ωs t − vds1 isd1 + sin ωs t − vds1 cos ωs t − iq1 + cos ωs t − vqs1 sin ωs t − id1 +
3 3 3 3 3 3
2π 2π 2π 2π s 2π 2π s
h             i
cos2 ωs t + vqs1 isq1 + sin2 ωs t + vds1 isd1 + sin ωs t + vds1 cos ωs t + iq1 + cos ωs t + vqs1 sin ωs t + id1
3 3 3 3 3 3

We employ the following trig identity: cos2 (x) + cos2 x + 2π x − 2π

2


3 + cos 3 = 1.5. This holds for the sine
and the cosine terms. We simplify this expression by leveraging this identity and gathering like terms.

p1 =1.5 vqs1 isq1 + 1.5 vds1 isd1 +

   
        
s s 2π 2π 2π 2π
iq1 vd1 sin (ωs t) cos (ωs t) + sin ωs t − cos ωs t − + sin ωs t + cos ωs t + +
3 3 3 3
        
2π 2π 2π 2π
isd1 vqs1 cos (ωs t) sin (ωs t) + cos ωs t − sin ωs t − + cos ωs t + sin ωs t +
3 3 3 3
3 s s
= vq1 iq1 + vds1 isd1


2
3 r r
= vq1 iq1 + vdr1 ird1


2

37
Appendix II - Time Varying Inductance Matrix
The induction matrix for the salient pole generator is given by equation (68). Its component sub-matrices
follow.  
Ls Lsr
L(δr ) = (68)
L>sr Lr

− 12 LA − LB cos 2 δr − π3
− 21 LA − LB cos 2 δr + π3



Lls + LA − LB cos 2δr

Ls =  − 12 LA − LB cos 2 δr − π3
Lls + LA − LB cos 2 δr − 2π 3 − 12 LA − LB cos 2 (δr + π)

− 12 LA − LB cos 2 δr + π3 − 12 LA − LB cos 2 (δr + π) Lls + LA − LB cos 2 δr + 2π 3
 
Llf d + Lmf d Lf dkd 0 0
 Lf dkd L lkd + L mkd 0 0 
Lr =  
 0 0 Llkq1 + Lmkq1 Lkq1kq2 
0 0 Lkq1kq2 Llkq2 + Lmkq2
 
Lsf d sin (δr + 0)
Lskd sin (δr + 0)
Lskq1 cos (δr + 0)
Lskq2 cos (δr + 0)

Lsr =  Lsf d sin δr − 2π 2π
3
Lskd sin δr − 3
Lskq1 cos δr − 3

2π
Lskq2 cos δr − 2π3


2π 2π
Lsf d sin δr + 3 Lskd sin δr + 3 Lskq1 cos δr + 2π3 Lskq2 cos δr + 2π3

We use this matrix to determine the flux linkages associated with each of the generator’s seven coils.

λ = L(δr )q̇

The vector λis the magnetic flux linkage associated with each of the seven coils, and q̇ is the corresponding
current flow through each coil. All values are instantaneous values.
 >
q̇ = i = −ia −ib −ic ikq1 ikq2 if d ikd
 >
λ= λa λb λc λkq1 λkq2 λf d λkd

Appendix III - Static Inductance Matrix

If we apply a Park’s Transform, it can be shown (see Krause) that a constant, time invariant matrix can be
used to link the flux linkage and the current variables. All variables are given in the reference frame of the
rotor δr .

   
λqs −iqs

Lls + Lmq 0 Lmq Lmq 0 0
 λds  
0 Lls + Lmd 0 0 Lmd Lmd  −ids 
0 0
    
 λkq1  
Lmq 0 Llf d + Lmf d Lf dkd 0 0  ikq1 
=
 0
   0

λkq2 Lmq 0 Lf dkd Llkd + Lmkd 0 0 ikq2
  
    
0 0
0 Lmd 0 0 Llkq1 + Lmkq1 Lkq1kq2
    
 λf d   if d 
0 0
λkd 0 Lmd 0 0 Lkq1kq2 Llkq2 + Lmkq2 ikd
(69)
λ̂ = L̂(δr )q̂˙

˙
We define λ̂ and q̂.

h 0 0 0 0
i>
q̂˙ = î = −iqs −ids ikq1 ikq2 if d ikd
h 0 0 0 0
i>
λ̂ = λqs λds λkq1 λkq2 λf d λkd

And we denote the constant inductance matrix as L̂(δr ):

38
  
Lls + Lmq 0 Lmq Lmq 0 0

 0 Lls + Lmd 0 0 Lmd Lmd 

 Lmq 0 Llf d + Lmf d Lf dkd 0 0 
L̂(δr ) =   (70)

 Lmq 0 Lf dkd Llkd + Lmkd 0 0 

 0 Lmd 0 0 Llkq1 + Lmkq1 Lkq1kq2 
0 Lmd 0 0 Lkq1kq2 Llkq2 + Lmkq2

Appendix IV - Resistive Matrix

The diagonal matrix of resistances for the rotor and stator coils are given below.
 
ra 0 0 0 0 0 0
 0
 rb 0 0 0 0 0  
 0
 0 rc 0 0 0 0  
R=  0 0 0 rf d 0 0 0   (71)
 0
 0 0 0 rkd 0 0  
 0 0 0 0 0 rkq1 0 
0 0 0 0 0 0 rkq2

39