Chapter 14 : Power System Stability

Example 14.1, Page No 380

In [1]:
import math
#initialisation of variables
WR2 = 5.82
Smach = 1333
n = 1800

#Calculations
ft_lb = 746 / 550.0
w = 2 * math.pi * n / 60
H = ft_lb * WR2 * w**2 / (2 * 32.2 * Smach)

#Results
print('The inertia constant in MJ/MVA is = %.2f' %H)
print('Converting H to a 100-MVA system base,units in MJ/MVA = %.2f' %(H * Smach / 100))

The inertia constant in MJ/MVA is = 3.27
Converting H to a 100-MVA system base,units in MJ/MVA = 43.55

Example 14.2, Page No 181

In [2]:
import math
#initialisation of variables
P1 = 500.0
pf1 = 0.85
V1 = 20.0
n1 = 3600
P2 = 1333.0
pf2 = 0.9
V2 = 22.0
n2 = 1800.0
Pbase = 100
H1 = 4.8
H2 = 3.27

#Calculations
KE = H1 * P1 + H2 * P2

#Results
print('The total kinetic energy of rotation of the two machines in MJ is = %.2f' %KE)
print('The inertia constant for the equivalent machine on 100-MVA base in MJ/MVA is = %.2f' %(KE/Pbase))

The total kinetic energy of rotation of the two machines in MJ is = 6758.91
The inertia constant for the equivalent machine on 100-MVA base in MJ/MVA is = 67.59

Example 14.3, Page No 386

In [3]:
import math
#initialisation of variables
Pm =1
Vt = 1.0
V_ib = 1.0
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4

#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm*X/(Vt*V_ib))*180/math.pi
Vt1 = Vt*complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1

#Results
print('The terminal voltage after considering the angle is,in per unit= {0:.2f}+{1:.2f}i'.format(Vt1.real, Vt1.imag))
print('The output current from the generator in per unit is= {0:.2f}+{1:.2f}i'.format(I.real, I.imag))
print('The transient internal voltage in per unit= {0:.2f}+{1:.2f}i'.format(E1.real, E1.imag))
print('Power angle equation is')
print("\n Pe = %.2f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pmax)

The terminal voltage after considering the angle is,in per unit= 0.95+0.30i
The output current from the generator in per unit is= -0.15+1.00i
The transient internal voltage in per unit= 0.92+0.50i
Power angle equation is

 Pe = 2.10 * sin(delta) 
 where delta is the machine rotor angle wrt to the infinite bus

Example 14.4 Page No 388

In [4]:
import math
#initialisation of variables
H = 5.0
Pm =1.0
Vt = 1.0
V_ib = 1.0
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt*complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
y10 = complex(3.33)
y32 = complex(2.5)
y30 = complex(5)
y20 = complex(5)
Ybus =  [[0,0,0],[0,0,0],[0,0,0]]

#Calculations
Ybus[0][0] = -y10
Ybus[0][1] =0
Ybus[0][2] = y10
Ybus[1][0] = Ybus[0][1]
Ybus[1][1] = -(y32 + y30)
Ybus[1][2] = y32
Ybus[2][0] = Ybus[0][2]
Ybus[2][1] = Ybus[1][2]
Ybus[2][2] = -(y10 + y32 + y30)
print('Ybus formed by inspection is')
print(Ybus)

Ybus_new = [[0,0],[0,0]]
for c in range(0, 2):
    for d in range(0, 2):
        Ybus_new[c][d] = Ybus[c][d] - ((Ybus[c][2]*Ybus[2][d]) / Ybus[2][2])
        
print('Ybus formed after elimination of Bus 3')
print(Ybus_new)
Pmax = abs(E1) * V_ib * abs(Ybus_new[0][1])
delta = 28.44
Pa = Pm - Pmax * math.sin(delta * math.pi / 180)
b = 180 * Pa / H
#Results
print('The power abgle equation is')
print(" Pe = %.3f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pmax)
print('The swing equation is')
print("(%.2f/180f) * d(delta)^2/dt^2 = %.2f - %.2fsin(delta) \n" %(H,Pm,Pmax))
print(" Intial Accelerating power is %.3f per unit \n" %(Pa))
print(" Initial acceleration is %.2f*f \n where f is the system frequency" %(b))

Ybus formed by inspection is
[[(-3.33-0j), 0, (3.33+0j)], [0, (-7.5-0j), (2.5+0j)], [(3.33+0j), (2.5+0j), (-10.83-0j)]]
Ybus formed after elimination of Bus 3
[[(-2.306094182825485+0j), (0.7686980609418281+0j)], [(0.7686980609418281+0j), (-6.922899353647276+0j)]]
The power abgle equation is
 Pe = 0.807 * sin(delta) 
 where delta is the machine rotor angle wrt to the infinite bus
The swing equation is
(5.00/180f) * d(delta)^2/dt^2 = 1.00 - 0.81sin(delta) 

 Intial Accelerating power is 0.616 per unit 

 Initial acceleration is 22.16*f 
 where f is the system frequency

Example 14.5, Page No 389

In [5]:
import math
#initialisation of variables
H = 5
Pm =1
Vt = 1
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4

#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
y12 = 1 / complex((X1_g + X1_t + X1_l1))
Y12 = -y12
Pe = abs(E1) * V_ib * abs(Y12)

#Results
print('The power abgle equation is')
print(" Pe = %.3f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pe)
print('The swing equation is')
print(" (%.2f/180f) * d(delta)^2/dt^2 = %.2f - %.2fsin(delta) \n" %(H,Pm,Pe))

The power abgle equation is
 Pe = 1.500 * sin(delta) 
 where delta is the machine rotor angle wrt to the infinite bus
The swing equation is
 (5.00/180f) * d(delta)^2/dt^2 = 1.00 - 1.50sin(delta)

Example 14.6, Page No 392

In [6]:
import math
#initialisation of variables
delta = 28.44
H = 5
ws = 377
Pm =1
Vt = 1.0
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4

#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1
Sp = Pmax * math.cos(delta * math.pi / 180)
wn = math.sqrt(ws * Sp / (2 * H))
fn = wn / (2 * math.pi)
T = 1 / fn

#Results
print("The angular frequency of oscillation is %.3f elec rad/s " %wn)
print("The corresponding frquency of oscillation is %.2f Hz " %fn)
print("The period of oscillation is %.3f s" %T)

The angular frequency of oscillation is 8.343 elec rad/s 
The corresponding frquency of oscillation is 1.33 Hz 
The period of oscillation is 0.753 s

Example 14.7 Page No 392

In [7]:
import math
#initialisation of variables
delta = 28.44
H = 5
ws = 377
Pm =1
Vt = 1
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4

#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1
delta_rad = delta * math.pi / 180
delta_cr = math.acos((math.pi - 2 * delta_rad) * math.sin(delta_rad) - math.cos(delta_rad))
t_cr = math.sqrt(4 * H * (delta_cr - delta_rad) / (ws * Pm))

#Results
print("Critical clearing angle = %.3f elec rad " %(delta_cr))
print("Critical clearing angle for the system = %.3f s" %(t_cr))

Critical clearing angle = 1.426 elec rad 
Critical clearing angle for the system = 0.222 s