Solution of Swing Equation

• 1. Step by Step Method: Home Work

• 2. Euler Method: Home Work
• 3. Modified Euler Method:
• The initial power angle is given by
𝑃𝑚
• 𝛿0 = sin−1
𝑃𝑚𝑎𝑥
• The rotor is running at synchronous speed and change in speed ∆𝜔=0
• The swing eqn.can be written in two First order differential eqn.
𝑑𝛿
• = ∆𝜔,
𝑑𝑡
𝑑∆𝜔 𝜋𝑓 𝜋𝑓
• = 𝑃𝑚 − 𝑃𝑒 = 𝑃 , where 𝑃𝑎 = Accelerating power
𝑑𝑡 𝐻 𝐻 𝑎
• These eqns. can be solved by predicting and subsequently correcting variables in each iteration.
• It is a two-step procedure for solving each of the solution state. At any state 𝑖 + 1, the following eqns.
need to be solved. The value at the end of the step (𝑡1 = 𝑡0 + ∆𝑡)
𝑝 𝑑𝛿 𝑝 𝑑∆𝜔
• 𝛿𝑖+1 = 𝛿𝑖 + . ∆𝑡 and ∆𝜔𝑖+1 = ∆𝜔𝑖 + . ∆𝑡
𝑑𝑡 ∆𝜔𝑖 𝑑𝑡 𝛿𝑖
𝑝 𝑝
• Using the predicted value (represented as subscript 𝑝 ), 𝛿𝑖+1 and∆𝜔𝑖+1 , we get for derivatives at the
end of interval
𝑑𝛿 𝑝 𝑑∆𝜔 𝜋𝑓
• 𝑑𝑡 ∆𝜔𝑝
=∆𝜔𝑖+1 ; 𝑑𝑡 𝛿 𝑝
= 𝐻 𝑃𝑚 − 𝑃𝑒 𝑝
𝑖+1 𝑖+1 𝛿𝑖+1
• Then the average value of two derivatives is used to find the corrected value (represented as subscript
𝑐 ),
𝑑𝛿 𝑑𝛿
+
𝑑𝑡 ∆𝜔 𝑑𝑡 ∆𝜔𝑝
𝑐 𝑖
• 𝛿𝑖+1 = 𝛿𝑖 + 𝑖+1
x ∆𝑡 and
2
𝑑∆𝜔 𝑑∆𝜔
+
𝑑𝑡 𝛿 𝑑𝑡 𝛿𝑝
𝑐 𝑖
• ∆𝜔𝑖+1 = ∆𝜔𝑖 + 2
𝑖+1
x ∆𝑡
𝑐 𝑐
• For the next step of operation, 𝛿𝑖+1 = 𝛿𝑖+1 and ∆𝜔𝑖+1 = ∆𝜔𝑖+1
• The process continues and in each step, it is possible to get the value of 𝛿𝑖𝑐 against the corresponding
state.
4. Runge-Kutta Method(2nd Order): home work
5. Runge-Kutta Method(4th Order)
• The swing equation can be written as
𝐻 𝑑2 𝛿
• = 𝑃𝑚 − 𝑃𝑚𝑎𝑥 𝑠𝑖𝑛𝛿 = 𝑃𝑚 − 𝑃𝑒 𝑃𝑚 input power assumed to be constant
𝜋𝑓𝑠 𝑑𝑡 2
• The above equation can be split into 2 single order equation,
𝑑𝛿 𝑑𝛿
• = ∆𝜔 = 0 − 𝜔0 = 𝜔 − 2𝜋𝑓0 = 𝑓 𝜔
𝑑𝑡 𝑑𝑡
𝑑∆𝜔 𝜋𝑓
• = 𝑃𝑚 − 𝑃𝑒 = 𝑓 𝛿
𝑑𝑡 𝐻
• We are interested in getting 𝛿 and 𝜔 as function of time t.
• Let us assume that 𝑥 = 𝑡, 𝑦 = 𝛿, 𝑧 = 𝜔
𝑑𝛿 𝑑𝑦
• Hence, =𝑓 𝜔 = =𝑓 𝑥, 𝑦, 𝑧 (2a)
𝑑𝑡 𝑑𝑥
𝑑∆𝜔 𝑑𝑧
• And, =𝑓 𝛿 = =𝑔 𝑥, 𝑦, 𝑧 (2b)
𝑑𝑡 𝑑𝑥
• The fist step in the transient stability study is the load flow equation to obtain system condition prior to
disturbance.
• Each bus must be characterized by 𝑃, 𝑄, 𝑉, 𝛿 𝑎𝑛𝑑 any 2 of them must be known.
• Let us assume that the 𝑖 − 𝑡ℎ machine is represented by a constant voltage source behind transient
′ ′
reactance 𝑥𝑑′ . 𝐸𝑖(0) = 𝐸𝑡𝑖 + 𝑟𝑎𝑖 + 𝑗𝑥𝑑𝑖 𝐼𝑡𝑖 (3)
where, 𝐸𝑡𝑖 = terminal voltage of 𝑖 − 𝑡ℎ machine .
𝑟𝑎𝑖 = armature resistance of 𝑖 − 𝑡ℎ machine.
𝐼𝑡𝑖 = current through 𝑖 − 𝑡ℎ machine.
′ ′ ′
Also 𝐸𝑖(0) = 𝑒𝑖(0) + 𝑗𝑓𝑖(0) in rectangular form. (4)
• The initial internal voltage angle(power angle) is given by
′
−1 𝑓𝑖(0)
• 𝛿𝑖(0) = 𝑡𝑎𝑛 ′ (5)
𝑒𝑖(0)

• Machine current prior to disturbance are calculated from

𝑃𝑡𝑖 −𝑗𝑄𝑡𝑖
• 𝐼𝑡𝑖 = for 𝑖 = 1,2,3, … … . . , 𝑚, (6)
𝐸𝑡𝑖
• where, 𝑃𝑡𝑖 𝑎𝑛𝑑 𝑄𝑡𝑖 are real and reactive power generation of 𝑖 − 𝑡ℎ machine and 𝐸𝑡𝑖 is the terminal
voltage of 𝑖 − 𝑡ℎ machine taken as ref. phasor.
• 𝑃𝑚𝑖(0) initial mechanical input power= electrical air gap power 𝑃𝑒𝑖 prior to
disturbance and can be obtained from
• 𝑃𝑒𝑖 = 𝑃𝑡𝑖 + 𝐼𝑟𝑖 2 𝑟𝑎𝑖 (7)
• For steady state, 𝛿0 = 𝛿𝑖(0) . For the transient analysis(when disturbance occurs
due to fault),
• The network admittance matrix must be now changed to include equivalent
circuit of machines and changes in loads.
• Each element representing a machine is a branch of a new bus and each element
representing a load is a link to ground.
• According to Runge-kutta derivation,it is required to solve 2 simultaneous
differential eqns. (2)
• General Expression of a Runge-kutta 4th order numerical technique
1
• 𝑦1 = 𝑦0 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4
6
1
• 𝑧1 = 𝑧0 + 𝑙1 + 2𝑙2 + 2𝑙3 + 𝑙4 where,
6
• 𝑘1 = 𝑓 𝑥0 , 𝑦0 , 𝑧0 ℎ, 𝑙1 = 𝑔 𝑥0 , 𝑦0 , 𝑧0 ℎ
ℎ 𝑘1 𝑙1 ℎ 𝑘1 𝑙1
• 𝑘2 = 𝑓 𝑥0 + , 𝑥0 + , 𝑧0 + ℎ, 𝑙2 = 𝑔 𝑥0 + , 𝑥0 + , 𝑧0 + ℎ
2 2 2 2 2 2
ℎ 𝑘2 𝑙2 ℎ 𝑘2 𝑙2
• 𝑘3 = 𝑓 𝑥0 + , 𝑥0 + , 𝑧0 + ℎ, 𝑙3 = 𝑔 𝑥0 + , 𝑥0 + , 𝑧0 + ℎ,
2 2 2 2 2 2
• 𝑘4 = 𝑓 𝑥0 + ℎ , 𝑥0 + 𝑘3 , 𝑧0 + 𝑙3 ℎ, 𝑙4 = 𝑔 𝑥0 + ℎ , 𝑥0 + 𝑘3 , 𝑧0 + 𝑙3 ℎ,

• According to Runge-kutta 4th order numerical technique,

1
• ∆𝛿𝑖(𝑡+∆𝑡) = 𝑘1𝑖 + 2𝑘2𝑖 + 2𝑘3𝑖 + 𝑘4𝑖 (8-a)
6
1
• And ∆𝜔𝑖(𝑡+∆𝑡) = 𝑙1𝑖 + 2𝑙2𝑖 + 2𝑙3𝑖 + 𝑙4𝑖 (8-b)
6
• 𝑘 ′ 𝑠 𝑎𝑛𝑑 𝑙 ′ 𝑠 are the changes in 𝛿𝑖 and 𝜔𝑖 respectively obtained using derivatives evaluated at
predetermined points.
• STEP 1: The initial estimates of changes are obtained from
 𝑘1𝑖 = 𝜔𝑖(𝑡) − 2𝜋𝑓 ∆𝑡
𝜋𝑓
 𝑙1𝑖 = ( ) 𝑃𝑚𝑖 − 𝑃𝑒𝑖(𝑡) ∆𝑡
𝐻𝑖
 At 𝑡 = 0, when disturbance just occurs, 𝜔𝑖(𝑡) =𝜔𝑖(0)= synchronous angular velocity(2𝜋𝑓)
 𝑃𝑒𝑖(𝑡) = elect. Power calculated at 𝛿𝑖(𝑡) = 𝛿𝑖(0)
 STEP 2: The second set of estimates of changes in 𝛿𝑖 and 𝜔𝑖 are obtained from
𝑙1𝑖
 𝑘2𝑖 = {𝜔𝑖(𝑡) + }− 2𝜋𝑓 ∆𝑡
2
𝜋𝑓
 𝑙2𝑖 = ( ) 𝑃𝑚𝑖 − 𝑃𝑒𝑖(1) ∆𝑡 for 𝑖 = 1,2,3, … … . . , 𝑚,
𝐻𝑖
𝑘1𝑖
 Thus, 𝑃𝑒𝑖(1)= elect. Power corresponding to power angle 𝛿𝑖(𝑡) =
2
 Before 𝑙2𝑖 can be calculated , new component for the voltage corresponding to machine buses
can be calculated from
′ ′ 𝑘1𝑖
 𝑒𝑖(1) = 𝐸𝑖(0) 𝑐𝑜𝑠 𝛿𝑖(𝑡) +
2
′ ′ 𝑘1𝑖
 𝑓𝑖(1) = 𝐸𝑖(0) 𝑠𝑖𝑛 𝛿𝑖(𝑡) +
2
 ′
𝑘1𝑖 𝐸𝑖(1) −𝐸𝑡𝑖(1) ′
• Knowing 𝛿𝑖(𝑡) + , 𝐼𝑡𝑖(1) = , 𝑤ℎ𝑒𝑟𝑒, 𝑍𝑖 = 𝑟𝑎𝑖 + 𝑗𝑥𝑑𝑖 ,
2 𝑍𝑖
𝐸𝑡𝑖 1 = terminal voltage of 𝑖−𝑡ℎ machine .
𝑘1𝑖
• When 𝛿𝑖 = 𝛿𝑖(𝑡) + which can be obtained from the new load flow study, we get
2
′ ∗
• 𝑃𝑒𝑖(1) =Real 𝐸𝑖(1) 𝐼𝑡𝑖(1) , Thus , 𝑙2𝑖 can be calculated.
• Step 3: The third set of estimates are obtained from
𝑙2𝑖
 𝑘3𝑖 = {𝜔𝑖(𝑡) + }− 2𝜋𝑓 ∆𝑡
2
𝜋𝑓
 𝑙3𝑖 = ( ) 𝑃𝑚𝑖 − 𝑃𝑒𝑖(2) ∆𝑡 for 𝑖 = 1,2,3, … … . . , 𝑚,
𝐻𝑖
𝑘2𝑖
• Where, 𝑃𝑒𝑖(2)=𝑓 𝛿 = 𝑓 𝛿𝑖(𝑡) + ∆𝑡 ,
2
• 𝑃𝑒𝑖(2) is obtained from 2nd solution of the network eqn. with internal voltage angles equal to
𝑘
𝛿𝑖(𝑡) + 2𝑖 and component of the voltage of machine buses equal to
2
′ ′ 𝑘2𝑖
 𝑒𝑖(2) = 𝐸𝑖(1) 𝑐𝑜𝑠 𝛿𝑖(𝑡) +
2
′ ′ 𝑘2𝑖
 𝑓𝑖(2) = 𝐸𝑖(1) 𝑠𝑖𝑛 𝛿𝑖(𝑡) +
2
• Step 4: The 4th set of estimates are obtained from
 𝑘4𝑖 = {𝜔𝑖(𝑡) +𝑙3𝑖 } − 2𝜋𝑓 ∆𝑡
𝜋𝑓
 𝑙4𝑖 = ( ) 𝑃𝑚𝑖 − 𝑃𝑒𝑖(3) ∆𝑡 for 𝑖 = 1,2,3, … … . . , 𝑚,
𝐻𝑖

• Where, 𝑃𝑒𝑖(3)is obtained from 3rd solution of the network eqn. with internal voltage angles equal
to 𝛿𝑖(𝑡) + 𝑘3𝑖 and component of the voltage of machine buses equal to
′ ′
 𝑒𝑖(3) = 𝐸𝑖(2) 𝑐𝑜𝑠 𝛿𝑖(𝑡) + 𝑘3𝑖
′ ′
 𝑓𝑖(3) = 𝐸𝑖(2) 𝑠𝑖𝑛 𝛿𝑖(𝑡) + 𝑘3𝑖
 Knowing different 𝑘 ′ 𝑠 𝑎𝑛𝑑 𝑙 ′ 𝑠 , find values of 𝛿 ′ 𝑠 𝑎𝑛𝑑 𝜔′ 𝑠 after an interval of ∆𝑡are computed.
 i.e. 𝛿𝑖(𝑡+∆𝑡) = 𝛿𝑖(𝑡) + ∆𝛿𝑖(𝑡+∆𝑡) and
 𝜔𝑖(𝑡+∆𝑡) = 𝜔𝑖(𝑡) + ∆𝜔𝑖(𝑡+∆𝑡)
 The numerical technique used for the solution of swing eqn are useful for computer calculation.