# CHAPTER 2.10: POWER SYSTEM STABILITY¶

## Example 2.10.1, Page number 270¶

In [1]:
import math

#Variable declaration
Z = 0.1        #Impedance of transmission line(p.u)
M = 0.3        #Stability margin
X = 1.0        #Constant(p.u)

#Calculation
sin_delta_0 = 1-M                               #Sin(δ_0)
delta_0 = math.asin(sin_delta_0)*180/math.pi    #δ_0(°)
P_0 = X/Z*sin_delta_0                           #Magnitude of P_0(p.u)

#Result
print('Operating power angle, δ_0 = %.2f° ' %delta_0)
print('P_0 = %.2f p.u' %P_0)

Operating power angle, δ_0 = 44.43°
P_0 = 7.00 p.u


## Example 2.10.2, Page number 270¶

In [1]:
#Variable declaration
x_s = 0.85     #Reactance(p.u)
x_T1 = 0.157   #Reactance(p.u)
x_T2 = 0.157   #Reactance(p.u)
x_l1 = 0.35    #Reactance(p.u)
x_l2 = 0.35    #Reactance(p.u)
E = 1.50       #Sending end voltage(p.u)
P_0 = 1.0      #Stable power output(p.u)

#Calculation
x = x_s+x_T1+x_T2+(x_l1/2)         #Total reactance(p.u)
P_max = E*V_L/x                    #Maximum power limit(p.u)
M = (P_max-P_0)/P_max*100          #Steady state stability margin(%)
V_Lmin = P_0*x/E                   #Minimum value of V_L(p.u)
E_min = P_0*x/V_L                  #Minimum value of E(p.u)

#Result
print('Minimum value of |E|, |E_min| = %.3f p.u' %E_min)
print('Minimum value of |V_L|, |V_Lmin| = %.3f p.u' %V_Lmin)
print('Maximum power limit, P_0 = %.2f p.u' %P_max)
print('Steady state stability margin, M = %.1f percent' %M)

Minimum value of |E|, |E_min| = 1.339 p.u
Minimum value of |V_L|, |V_Lmin| = 0.893 p.u
Maximum power limit, P_0 = 1.12 p.u
Steady state stability margin, M = 10.7 percent


## Example 2.10.3, Page number 270-271¶

In [1]:
#Variable declaration
E_1 = 1.25     #Sending end voltage(p.u)
x_d = 1.0      #Reactance(p.u)
x_T1 = 0.2     #Reactance(p.u)
x_l1 = 1.0     #Reactance(p.u)
x_l2 = 1.0     #Reactance(p.u)
x_T2 = 0.2     #Reactance(p.u)
E_2 = 1.0      #Receiving end voltage(p.u)
x_L = 1.0      #Shunt inductor reactance(p.u)
x_C = 1.0      #Shunt capacitor reactance(p.u)

#Calculation
#Case(a)
Z_1_a = x_d+x_T1+(x_l1/2.0)              #Reactance(p.u)
Z_2_a = x_T2+x_d                         #Reactance(p.u)
Z_3_a = x_L                              #Reactance(p.u)
Z_a = Z_1_a+Z_2_a+(Z_1_a*Z_2_a/Z_3_a)    #Transfer reactance(p.u)
P_max_1 = E_1*E_2/Z_a                    #Maximum power transfer if shunt inductor is connected at bus 2(p.u)
#Case(b)
Z_1_b = x_d+x_T1+(x_l1/2.0)              #Reactance(p.u)
Z_2_b = x_T2+x_d                         #Reactance(p.u)
Z_3_b = -x_C                             #Reactance(p.u)
Z_b = Z_1_b+Z_2_b+(Z_1_b*Z_2_b/Z_3_b)    #Transfer reactance(p.u)
P_max_2 = E_1*E_2/Z_b                    #Maximum power transfer if shunt capacitor is connected at bus 2(p.u)

#Result
print('Case(a): Maximum power transfer if shunt inductor is connected at bus 2, P_max1 = %.3f p.u' %P_max_1)
print('Case(b): Maximum power transfer if shunt capacitor is connected at bus 2, P_max2 = %.2f p.u' %P_max_2)

Case(a): Maximum power transfer if shunt inductor is connected at bus 2, P_max1 = 0.253 p.u
Case(b): Maximum power transfer if shunt capacitor is connected at bus 2, P_max2 = 1.45 p.u


## Example 2.10.4, Page number 271¶

In [1]:
import math

#Variable declaration
V = 400.0          #Voltage(kV)
L = 220.0          #Line length(km)
P = 0.58           #Initial real power transfer(p.u)
PF = 0.85          #Lagging power factor
V_L = 1.00         #Load bus voltage(p.u)
x_d = 0.460        #Reactance(p.u)
x_T1 = 0.200       #Reactance(p.u)
x_T2 = 0.15        #Reactance(p.u)
x_line = 0.7       #Reactance(p.u)

#Calculation
x = x_d+x_T1+x_T2+(x_line/2)                  #Net reactance(p.u)
phi = math.acos(PF)*180/math.pi               #Φ(°)
Q = P*math.tan(phi*math.pi/180)               #Reactive power(p.u)
E = ((V_L+(Q*x/V_L))**2+(P*x/V_L)**2)**0.5    #Excitation voltage of generator(p.u)
P_max = E*V_L/x                               #Maximum power transfer(p.u)
M = (P_max-P)/P_max*100                       #Steady state stability margin(%)

#Result
print('Maximum power transfer, P_max = %.2f p.u' %P_max)
print('Stability margin, M = %.f percent' %M)

Maximum power transfer, P_max = 1.35 p.u
Stability margin, M = 57 percent


## Example 2.10.5, Page number 271-272¶

In [1]:
import math

#Variable declaration
V_A = 1.0          #Voltage at bus A(p.u)
Z_AB = 1j*0.5      #Impedance(p.u)
S_DA = 1.0         #p.u
S_DB = 1.0         #p.u
V_B = 1.0          #Voltage at bus B(p.u)

#Calculation
#Case(i) & (ii)
X = abs(Z_AB)                                              #Reactance(p.u)
sin_delta = 1.0*X/(V_A*V_B)                                #Sin δ
delta = math.asin(sin_delta)*180/math.pi                   #δ(°)
V_2 = V_B
V_1 = V_A
Q_gB = (V_2**2/X)-(V_2*V_1*math.cos(delta*math.pi/180)/X)
#Case(iii)
V_2_3 = 1/2.0**0.5                                         #Solving quadratic equation from textbook
delta_3 = math.acos(V_2_3)*180/math.pi                     #δ(°)

#Result
print('Case(i)  : Q_gB = %.3f' %Q_gB)
print('Case(ii) : Phase angle of V_B, δ = %.f° ' %delta)
print('Case(iii): If Q_gB is equal to zero then amount of power transmitted is, V_2 = %.3f∠%.f° ' %(V_2_3,delta_3))

Case(i)  : Q_gB = 0.268
Case(ii) : Phase angle of V_B, δ = 30°
Case(iii): If Q_gB is equal to zero then amount of power transmitted is, V_2 = 0.707∠45°


## Example 2.10.6, Page number 272¶

In [1]:
import math
import cmath

#Variable declaration
A = 0.98*cmath.exp(1j*0.3*math.pi/180)       #Constant
B = 82.5*cmath.exp(1j*76.0*math.pi/180)      #Constant(ohm)
C = 0.0005*cmath.exp(1j*90.0*math.pi/180)    #Constant(mho)
D = A                                        #Constant
V_S = 110.0                                  #Sending end voltage(kV)
V_R = 110.0                                  #Receiving end voltage(kV)

#Calculation
alpha = cmath.phase(A)*180/math.pi                                                  #α(°)
beta = cmath.phase(B)*180/math.pi                                                   #β(°)
P_max = (V_S*V_R/abs(B))-(abs(A)*V_R**2/abs(B)*math.cos((beta-alpha)*math.pi/180))  #Maximum power transfer(MW)
B_new = abs(B)*math.sin(beta*math.pi/180)                                           #Constant(ohm)
beta_new = 90.0                                                                     #β(°)
P_max_new = (V_S*V_R/B_new)-(V_R**2/B_new*math.cos(beta_new*math.pi/180))           #Maximum power transfer(MW)

#Result
print('Steady state stability limit, P_max = %.1f MW' %P_max)
print('Steady state stability limit if shunt admittance is zero & series resistance neglected, P_max = %.2f MW' %P_max_new)
print('\nNOTE: Changes in the obtained answer from that of textbook is due to precision')

Steady state stability limit, P_max = 111.2 MW
Steady state stability limit if shunt admittance is zero & series resistance neglected, P_max = 151.16 MW

NOTE: Changes in the obtained answer from that of textbook is due to precision


## Example 2.10.8, Page number 273-275¶

In [1]:
import math
%matplotlib inline
from matplotlib.pyplot import plot,title,xlabel,ylabel,text,axis,grid,show

#Variable declaration
V = 33.0*10**3        #Line voltage(V)
R = 6.0               #Resistance per phase(ohm)
X = 15.0              #Reactance per phase(ohm)

#Calculation
V_S = V/3**0.5                                                              #Sending end phase voltage(V)
V_R = V/3**0.5                                                              #Receiving end phase voltage(V)
beta = math.atan(X/R)*180/math.pi                                           #β(°)
Z = (R**2+X**2)**0.5                                                        #Impedance(ohm)
delta_0 = 0.0                                                               #δ(°)
delta_1 = 30.0                                                              #δ(°)
delta_2 = 60.0                                                              #δ(°)
delta_3 = beta                                                              #δ(°)
delta_4 = 90.0                                                              #δ(°)
delta_5 = 120.0                                                             #δ(°)
delta_6 = (math.acos(R/Z)*180/math.pi)+beta                                 #δ(°)
plot(
[delta_0,delta_1,delta_2,delta_3,delta_4,delta_5,delta_6],
[P_0,P_1,P_2,P_3,P_4,P_5,P_6],
marker='o',color='b'
)
text(70,14.12,'P_max = 14.12 MW/phase(app.)',color='b')
axis([0,150,0,16])
title('Plot of Power angle diagram')
xlabel('Electrical degree')
ylabel('Power in MW/phase')
grid()
show()
P_max = V_R/Z**2*(V_S*Z-V_R*R)/10**6                                                  #Maximum power transmitted(MW/phase)
delta_equal = 0.0                                                                     #δ With no phase shift(°)
P_no_shift = (V_R/Z**2)*(V_S*Z*math.cos((delta_equal-beta)*math.pi/180)-V_R*R)/10**6  #Power transmitted with no phase shift(MW/phase)

#Result
print('Maximum power the line is capable of transmitting, P_max = %.2f MW/phase' %P_max)
print('With equal voltage at both ends power transmitted = %.f MW/phase' %abs(P_no_shift))

Maximum power the line is capable of transmitting, P_max = 14.12 MW/phase
With equal voltage at both ends power transmitted = 0 MW/phase


## Example 2.10.9, Page number 275¶

In [1]:
#Variable declaration
V = 132.0*10**3              #Sending end voltage(V)
Z_line = complex(4,6)        #Line impedance per phase(ohm)

#Calculation
V_S = V/3**0.5                                         #Sending end phase voltage(V)
V_R = V/3**0.5                                         #Receiving end phase voltage(V)
Z = abs(Z_line)                                        #Impedance(ohm)
R = Z_line.real                                        #Resistance per phase(ohm)
P_max_phase = ((V_S*V_R/Z)-(R*V_R**2/Z**2))/10**6      #Maximum steady state power that can be transmitted over the line(MW/phase)
P_max_total = 3.0*P_max_phase                          #Maximum steady state power that can be transmitted over the line(MW)

#Result
print('Maximum steady state power that can be transmitted over the line, P_max = %.f MW (total)' %P_max_total)

Maximum steady state power that can be transmitted over the line, P_max = 1076 MW (total)


## Example 2.10.10, Page number 275-276¶

In [1]:
import math

#Variable declaration
E_1 = 1.1      #Sending end voltage(p.u)
x_d1 = 1.0     #Reactance(p.u)
x_T1 = 0.1     #Reactance(p.u)
x_l1 = 0.4     #Reactance(p.u)
x_l2 = 0.4     #Reactance(p.u)
x_T2 = 0.1     #Reactance(p.u)
E_2 = 1.0      #Receiving end voltage(p.u)
x_d2 = 1.0     #Reactance(p.u)
x_L = 1.0      #Shunt inductor reactance(p.u)
x_C = 1.0      #Static capacitor reactance(p.u)
delta = 30.0   #δ(°)

#Calculation
#Case(a)
Z_1_a = x_d1+x_T1+(x_l1/2.0)                                           #Reactance(p.u)
X_1_a = 1j*Z_1_a
Z_2_a = x_T2+x_d2                                                      #Reactance(p.u)
X_2_a = 1j*Z_2_a
Z_3_a = -x_C                                                           #Reactance(p.u)
X_3_a = 1j*Z_3_a
X_a = X_1_a+X_2_a+(X_1_a*X_2_a/X_3_a)                                  #Transfer reactance(p.u)
P_max_a = E_1*E_2/abs(X_a)                   #Maximum steady state power if static capacitor is connected(p.u)
P_a = P_max_a*math.sin(delta*math.pi/180)                               #Value of P(p.u)
Q_a = (E_1*E_2/abs(X_a))*math.cos(delta*math.pi/180)-(E_2**2/abs(X_a))  #Value of Q(p.u)
#Case(b)
Z_1_b = x_d1+x_T1+(x_l1/2.0)                                            #Reactance(p.u)
X_1_b = 1j*Z_1_b
Z_2_b = x_T2+x_d2                                                       #Reactance(p.u)
X_2_b = 1j*Z_2_b
Z_3_b = x_L                                                             #Reactance(p.u)
X_3_b = 1j*Z_3_b
X_b = X_1_b+X_2_b+(X_1_b*X_2_b/X_3_b)                                   #Transfer reactance(p.u)
P_max_b = E_1*E_2/abs(X_b)                   #Maximum steady state power if static capacitor is replaced by an inductive reactor(p.u)
P_b = P_max_b*math.sin(delta*math.pi/180)                               #Value of P(p.u)
Q_b = (E_1*E_2/abs(X_b))*math.cos(delta*math.pi/180)-(E_2**2/abs(X_b))  #Value of Q(p.u)

#Result
print('Case(a): Maximum steady state power if static capacitor is connected, P_max = %.3f p.u' %P_max_a)
print('         Value of P = %.3f p.u' %P_a)
print('         Value of Q = %.3f p.u' %Q_a)
print('Case(b): Maximum steady state power if static capacitor is replaced by an inductive reactor, P_max = %.3f p.u' %P_max_b)
print('         Value of P = %.3f p.u' %P_b)
print('         Value of Q = %.4f p.u' %Q_b)

Case(a): Maximum steady state power if static capacitor is connected, P_max = 1.134 p.u
Value of P = 0.567 p.u
Value of Q = -0.049 p.u
Case(b): Maximum steady state power if static capacitor is replaced by an inductive reactor, P_max = 0.287 p.u
Value of P = 0.144 p.u
Value of Q = -0.0124 p.u


## Example 2.10.11, Page number 303¶

In [1]:
#Variable declaration
f = 50.0        #Frequency(Hz)
G = 100.0       #Rating of generator(MVA)
H = 5.0         #Inertia constant(MJ/MVA)
P_a = 20.0      #Acceleration power(MVA)

#Calculation
GH = G*H                #Energy stored in rotor at synchronous speed(MJ)
M = GH/(180*f)          #Angular momentum
acceleration = P_a/M    #Acceleration(°/sec^2)

#Result
print('Kinetic energy stored in the rotor at synchronous speed, GH = %.f MJ' %GH)
print('Acceleration = %.f°/sec^2' %acceleration)

Kinetic energy stored in the rotor at synchronous speed, GH = 500 MJ
Acceleration = 360°/sec^2


## Example 2.10.12, Page number 303-304¶

In [1]:
import math

#Variable declaration
f = 50.0          #Frequency(Hz)
P = 4.0           #Number of poles
G = 20.0          #Rating of generator(MVA)
H = 9.0           #Inertia constant(kWsec/MVA)
P_m = 26800.0     #Rotational loss(hp)
P_e = 16000.0     #Electric power developed(kW)

#Calculation
GH = G*H                                       #Energy stored in rotor at synchronous speed(MJ)
P_m_kW = P_m*0.746                             #Rotational loss(kW)
P_a = P_m_kW-P_e                               #Acceleration power(kW)
P_a1 = P_a/1000.0                              #Acceleration power(MW)
M = GH/(180*f)                                 #Angular momentum
acceleration = P_a1/M                          #Acceleration(°/sec^2)

#Result
print('Kinetic energy stored in the rotor at synchronous speed, GH = %.f MJ' %GH)
print('Acceleration = %.f°/sec^2 = %.2f rad/sec^2' %(acceleration,acceleration_1))
print('\nNOTE: ERROR: H = 9 kW-sec/MVA, not 9 kW-sec/kVA as mentioned in the textbook statement')

Kinetic energy stored in the rotor at synchronous speed, GH = 180 MJ
Acceleration = 200°/sec^2 = 3.48 rad/sec^2

NOTE: ERROR: H = 9 kW-sec/MVA, not 9 kW-sec/kVA as mentioned in the textbook statement


## Example 2.10.13, Page number 304¶

In [1]:
import math

#Variable declaration
f = 50.0           #Frequency(Hz)
P = 4.0            #Number of poles
alpha = 200.0      #Acceleration(°/sec^2)
n = 10.0           #Number of cycle

#Calculation
t = 1/f*n                                #Time(sec)
delta = delta_rel*t**2                   #Change in torque angle(rad)
delta_deg = delta*180/math.pi            #Change in torque angle in that period(°)
speed_rotor = (120*f/P)+rpm              #Rotor speed at the end of 10 cycles(r.p.m)

#Result
print('Change in torque angle in that period, δ = %.4f rad = %.f elect degree' %(delta,delta_deg))
print('Rotor speed at the end of 10 cycles = %.2f r.p.m' %speed_rotor)

Change in torque angle in that period, δ = 0.0698 rad = 4 elect degree
Rotor speed at the end of 10 cycles = 1503.33 r.p.m


## Example 2.10.14, Page number 304¶

In [1]:
import math

#Variable declaration
Power = 20.0*10**3       #Rating of generator(kVA)
PF = 0.8                 #Lagging power factor
fault = 0.5              #Reduction in output under fault
P = 4.0                  #Number of poles
f = 50.0                 #Frequency(Hz)

#Calculation
P_m = Power*PF           #Output power before fault(kW)
P_e = fault*P_m          #Output after fault(kW)
P_a = P_m-P_e            #Accelerating power(kW)
w_s = 4.0*math.pi*f/P    #Speed
T_a = P_a*10**3/w_s      #Accelerating torque at the time the fault occurs(N-m)

#Result
print('Accelerating torque at the time the fault occurs, T_a = %.2f N-m' %T_a)

Accelerating torque at the time the fault occurs, T_a = 50929.58 N-m


## Example 2.10.16, Page number 304-305¶

In [1]:
#Variable declaration
S = 1000.0          #Rating of generator(MVA)
N = 1500.0          #Speed of alternator(r.p.m)
WR_sq = 5.0*10**6   #WR^2(lb.ft^2)

#Calculation
H = 2.31*10**-10*WR_sq*N**2/S      #Inertia constant(MJ/MVA)
H_100 = H*1000.0/100               #Inertia constant on 100 MVA(MJ/MVA)

#Result
print('Value of inertia constant, H = %.1f MJ/MVA' %H)
print('Value of inertia constant in 100 MVA base, H = %.f MJ/MVA' %H_100)

Value of inertia constant, H = 2.6 MJ/MVA
Value of inertia constant in 100 MVA base, H = 26 MJ/MVA


## Example 2.10.17, Page number 305¶

In [1]:
#Variable declaration
MVA_1 = 500.0        #Rating of generator(MVA)
H_1 = 4.0            #Inertia constant(MJ/VA)
MVA_2 = 1000.0       #Rating of generator(MVA)
H_2 = 3.5            #Inertia constant(MJ/VA)
MVA = 100.0          #Base MVA

#Calculation
KE_T = H_1*MVA_1+H_2*MVA_2     #Total KE of the system(MJ)
H_total = KE_T/MVA             #Equivalent H for the two to common 100MVA base(MJ/MVA)

#Result
print('Equivalent H for the two to common 100 MVA base, H = %.f MJ/MVA' %H_total)

Equivalent H for the two to common 100 MVA base, H = 55 MJ/MVA


## Example 2.10.18, Page number 305¶

In [1]:
import math

#Variable declaration
MVA = 210.0           #Rating of generator(MVA)
P = 2.0               #Number of poles
f = 50.0              #Frequency(Hz)
MI = 60.0*10**3       #Moment of inertia(kg-mt^2)

#Calculation
N = 120.0*f/P                            #Speed(r.p.m)
KE = 1.0/2*MI*(2*math.pi*N/f)**2/10**6   #Energy stored in the rotor at rated speed(MJ)
H = KE/MVA                               #Inertia constant(MJ/MVA)
G = MVA
M = G*H/(180*f)                          #Angular momentum(MJ-sec/elect.degree)

#Result
print('Energy stored in the rotor at the rated speed, KE = %.2e MJ' %KE)
print('Value of inertia constant, H = %.2f MJ/MVA' %H)
print('Angular momentum, M = %.3f MJ-sec/elect.degree' %M)

Energy stored in the rotor at the rated speed, KE = 4.26e+03 MJ
Value of inertia constant, H = 20.30 MJ/MVA
Angular momentum, M = 0.474 MJ-sec/elect.degree


## Example 2.10.19, Page number 305¶

In [1]:
#Variable declaration
P_accl = 30.0         #Acceleration power(MVA)
M = 0.474             #Angular momentum(MJ-sec/elect.degree). From Example 10.18

#Calculation
acceleration = P_accl/M      #Acceleration of the rotor(elect.degree/sec^2)

#Result
print('Acceleration of the rotor = %.2f elect.degree/sec^2' %acceleration)

Acceleration of the rotor = 63.29 elect.degree/sec^2


## Example 2.10.20, Page number 305¶

In [1]:
import math

#Variable declaration
MVA = 50.0           #Rating of alternator(MVA)
P = 4.0              #Number of poles
f = 50.0             #Frequency(Hz)
KE = 150.0           #Kinetic energy stored in rotor(MJ)
P_m = 25.0           #Machine input(MW)
P_e = 22.5           #Developed power(MW)
n = 10.0             #Number of cycles

#Calculation
P_a = P_m-P_e            #Accelerating power(MW)
H = KE/MVA               #Inertia constant(MJ/MVA)
G = MVA
M_deg = G*H/(180*f)      #Angular momentum(MJ-sec/elect.degree)
t = 1/f*n                #Time(sec)
delta = 1.309*t**2       #Term in δ

#Result
print('Accelerating power = %.3f rad/sec^2' %acceleration)
print('New power angle after 10 cycles, δ = (%.3f + δ_0) rad' %delta)

Accelerating power = 2.618 rad/sec^2
New power angle after 10 cycles, δ = (0.052 + δ_0) rad


## Example 2.10.21, Page number 305-306¶

In [1]:
#Variable declaration
f = 50.0        #Frequency(Hz)
P = 4.0         #Number of poles
G = 20.0        #Rating of turbo-generator(MVA)
V = 13.2        #Voltage(kV)
H = 9.0         #Inertia constant(kW-sec/kVA)
P_s = 20.0      #Input power less rotational loss(MW)
P_e = 15.0      #Output power(MW)

#Calculation
KE = G*H                         #Kinetic energy stored(MJ)
M = G*H/(180*f)                  #Angular momentum(MJ-sec/elect.degree)
P_a = P_s-P_e                    #Accelerating power(MW)
alpha = P_a/M                    #Acceleration(elect.degree/sec^2)
alpha_deg = alpha/2.0            #Acceleration(degree/sec^2)
alpha_rpm = 60.0*alpha_deg/360   #Acceleration(rpm/sec)

#Result
print('Case(a): Kinetic energy stored by rotor at synchronous speed, GH = %.f MJ' %KE)
print('Case(b): Acceleration, α = %.f degree/sec^2' %alpha_deg)
print('         Acceleration, α = %.2f rpm/sec' %alpha_rpm)

Case(a): Kinetic energy stored by rotor at synchronous speed, GH = 180 MJ
Case(b): Acceleration, α = 125 degree/sec^2
Acceleration, α = 20.83 rpm/sec


## Example 2.10.22, Page number 306¶

In [1]:
#Variable declaration
f = 50.0        #Frequency(Hz)
P = 4.0         #Number of poles
G = 20.0        #Rating of turbo-generator(MVA)
V = 13.2        #Voltage(kV)
H = 9.0         #Inertia constant(kW-sec/kVA)
P_s = 20.0      #Input power less rotational loss(MW)
P_e = 15.0      #Output power(MW)
n = 10.0        #Number of cycles

#Calculation
KE = G*H                         #Kinetic energy stored(MJ)
M = G*H/(180*f)                  #Angular momentum(MJ-sec/elect.degree)
P_a = P_s-P_e                    #Accelerating power(MW)
alpha = P_a/M                    #Acceleration(elect.degree/sec^2)
alpha_deg = alpha/2.0            #Acceleration(degree/sec^2)
alpha_rpm = 60.0*alpha_deg/360   #Acceleration(rpm/sec)
t = 1.0/f*n                      #Time(sec)
delta = 1.0/2*alpha*t**2         #Change in torque angle(elect.degree)
N_s = 120*f/P                    #Synchronous speed(rpm)
speed = N_s+alpha_rpm*t          #Speed at the end of 10 cycles(rpm)

#Result
print('Change in torque angle in that period, δ = %.f elect degrees.' %delta)
print('Speed in rpm at the end of 10 cycles = %.2f rpm' %speed)

Change in torque angle in that period, δ = 5 elect degrees.
Speed in rpm at the end of 10 cycles = 1504.17 rpm


## Example 2.10.23, Page number 306¶

In [1]:
import math

#Variable declaration
G = 20.0        #Rating of turbo-generator(MVA)
PF = 0.75       #Lagging power factor
fault = 0.5     #Fault reduces output power
N_s = 1500.0    #Synchronous speed(rpm). From Example 10.22

#Calculation
P_prefault = PF*G         #Pre-fault output power(MW)
P_a = P_prefault*fault    #Post-fault output power(MW)
T_a = P_a*10**6/w         #Accelerating torque at the time of fault occurrence(N-m)

#Result
print('Accelerating torque at the time of fault occurrence, T_a = %.f N-m' %T_a)

Accelerating torque at the time of fault occurrence, T_a = 47746 N-m


## Example 2.10.24, Page number 306-307¶

In [1]:
import math
import cmath
from sympy import Symbol

#Variable declaration
x_d = 1j*0.2        #Transient reactance of generator(p.u)
P_e = 0.8           #Power delivered(p.u)
V_t = 1.05          #Terminal voltage(p.u)
H = 4.0             #Inertia constant(kW-sec/kVA)
x_t = 1j*0.1        #Transformer reactance(p.u)
x_l = 1j*0.4        #Transmission line reactance(p.u)
V = 1.0             #Infinite bus voltage(p.u)
f = 50.0            #Frequency(Hz)

#Calculation
x_12 = x_d+x_t+(x_l/2)                                     #Reactance b/w bus 1 & 2(p.u)
y_12 = 1/x_12                                              #Admittance b/w bus 1 & 2(p.u)
y_21 = y_12                                                #Admittance b/w bus 2 & 1(p.u)
y_10 = 0.0                                                 #Admittance b/w bus 1 & 0(p.u)
y_20 = 0.0                                                 #Admittance b/w bus 2 & 0(p.u)
Y_11 = y_12+y_10                                           #Admittance at bus 1(p.u)
Y_12 = -y_12                                               #Admittance b/w bus 1 & 2(p.u)
Y_21 = -y_12                                               #Admittance b/w bus 2 & 1(p.u)
Y_22 = y_21+y_20                                           #Admittance at bus 2(p.u)
x_32 = x_t+(x_l/2)                                         #Reactance b/w bus 3 & 1(p.u)
theta_t = math.asin(P_e*abs(x_32)/V_t)*180/math.pi         #Angle(°)
V_t1 = V_t*cmath.exp(1j*theta_t*math.pi/180)               #Terminal voltage(p.u)
I = (V_t1-V)/x_32                                          #Current(p.u)
E = V_t1+I*x_d                                             #Alternator voltage(p.u)
sin = Symbol('sin δ')
P_e1 = 2.0*abs(E)*sin                                      #Developed power(p.u)
P_m_P_e = P_e-P_e1
M = 2*H/(2*math.pi*f)                                      #Angular momentum
acc = (P_e-P_e1)*2*math.pi*f/(2*H)                         #Acceleration = α (rad/sec^2)

#Result
print('Swing equation is,  %.4f*α = ' %(M) + repr(P_m_P_e))
print('\nNOTE: Swing equation is simplified and represented here')
print('      ERROR: x_d = 0.2 p.u, not 0.1 p.u as mentioned in textbook statement')

Swing equation is,  0.0255*α = -2.22295398614037*sin δ + 0.8

NOTE: Swing equation is simplified and represented here
ERROR: x_d = 0.2 p.u, not 0.1 p.u as mentioned in textbook statement


## Example 2.10.26, Page number 308-309¶

In [1]:
import math

#Variable declaration
X_d = 0.25         #Transient reactance of generator(p.u)
X_t1 = 0.15        #Reactance of transformer(p.u)
X_t2 = 0.15        #Reactance of transformer(p.u)
X_t3 = 0.15        #Reactance of transformer(p.u)
X_t4 = 0.15        #Reactance of transformer(p.u)
X_l1 = 0.20        #Reactance of line(p.u)
X_l2 = 0.20        #Reactance of line(p.u)
X_tr = 0.15        #Reactance of transformer(p.u)
P_m = 1.0          #Power delivered(p.u)
E = 1.20           #Voltage behind transient reactance(p.u)
V = 1.0            #Infinite bus voltage(p.u)

#Calculation
X_14 = X_d+((X_t1+X_t2+X_l1)/2)+X_tr                       #Reactance before fault(p.u)
x_1_b = X_t1+X_t2+X_l1                                     #Reactance(p.u). From figure (b)
x_2_b = X_l2+X_t4                                          #Reactance(p.u). From figure (b)
x_1 = x_1_b*X_t3/(x_1_b+x_2_b+X_t3)                        #Reactance(p.u). From figure (c)
x_2 = x_1_b*x_2_b/(x_1_b+x_2_b+X_t3)                       #Reactance(p.u). From figure (c)
x_3 = X_t3*x_2_b/(x_1_b+x_2_b+X_t3)                        #Reactance(p.u). From figure (c)
X_14_fault = x_1+X_d+x_2+X_tr+((x_1+X_d)*(x_2+X_tr)/x_3)   #Reactance under fault(p.u)
X_14_after_fault = X_d+X_t1+X_l1+X_t2+X_tr                 #Reactance after fault is cleared(p.u)
P_max = V*E/X_14                                           #Maximum power transfer(p.u)
gamma_1 = (V*E/X_14_fault)/P_max                           #γ_1
gamma_2 = (V*E/X_14_after_fault)/P_max                     #γ_2
delta_0_degree = delta_0*180/math.pi                       #δ_0(°)
delta_m_degree = delta_m*180/math.pi                       #δ_1(°)
delta_c_degree = delta_c*180/math.pi                       #Clearing angle(°)

#Result
print('Critical clearing angle, δ_c = %.2f° ' %delta_c_degree)

Critical clearing angle, δ_c = 58.57°


## Example 2.10.27, Page number 309-310¶

In [1]:
import math

#Variable declaration
f = 50.0               #Frequency(Hz)
P_m = 1.0              #Power delivered(p.u)
P_max = 1.8            #Maximum power(p.u)
gamma_1_P_max = 0.4    #Reduced maximum power after fault(p.u)
gamma_2_P_max = 1.30   #Maximum power after fault clearance(p.u)

#Calculation
delta_0_degree = delta_0*180/math.pi                       #δ_0(°)
delta_f_degree = delta_f*180/math.pi                       #δ_1(°)
gamma_1 = gamma_1_P_max/P_max                              #γ_1
gamma_2 = gamma_2_P_max/P_max                              #γ_2
delta_c_degree = delta_c*180/math.pi                       #Clearing angle(°)

#Result
print('Critical angle, δ_c = %.2f° ' %delta_c_degree)

Critical angle, δ_c = 55.35°


## Example 2.10.28, Page number 310¶

In [1]:
import math

#Variable declaration
sin_delta_0 = 0.45     #Supplying percent of peak power capacity before fault
x = 4.0                #Reactance under fault increased
gamma_2 = 0.7          #Peak power delivered after fault clearance

#Calculation
delta_0_degree = delta_0*180/math.pi                    #δ_0(°)
gamma_1 = 1.0/x                                         #γ_1
delta_m_degree = delta_m*180/math.pi                    #δ_m(°)
delta_c_degree = delta_c*180/math.pi                    #Clearing angle(°)

#Result
print('Critical clearing angle, δ_c = %.f° ' %delta_c_degree)

Critical clearing angle, δ_c = 73°


## Example 2.10.30, Page number 310-311¶

In [1]:
import math

#Variable declaration
f = 60.0          #Frequency(Hz)
P = 6.0           #Number of poles
H = 4.0           #Inertia constant(p.u)
P_e = 1.0         #Power supplied by generator(p.u)
E = 1.2           #Internal voltage(p.u)
V = 1.0           #Infinite bus voltage(p.u)
X = 0.3           #Line reactance(p.u)
del_t = 0.05      #Δt = Interval step size(sec)

#Calculation
P_max = E*V/X                                #Maximum power(p.u)
delta_0 = math.asin(P_e/P_max)*180/math.pi   #δ_0(°)
G = P_e
M = G*H/(180*f)                              #Angular momentum(p.u)
P_a_0 = 1.0/2*(P_e-0)                        #(p.u)
alpha_0 = P_a_0/M                            #α_0(°/sec^2)
del_w_r_1 = alpha_0*del_t                    #Δω_r_1(°/sec)
w_r_1 = 0+del_w_r_1                          #ω_r_1(°/sec)
del_delta_1 = w_r_1*del_t                    #Δδ_1(°)
delta_1 = delta_0+del_delta_1                #δ_1(°)
P_a_1 = 1.0*(P_e-0)                          #(p.u)
alpha_1 = P_a_1/M                            #α_1(°/sec^2)
del_w_r_2 = alpha_1*del_t                    #Δω_r_2(°/sec)
w_r_2 = del_w_r_1+del_w_r_2                  #ω_r_2(°/sec)
del_delta_2 = w_r_2*del_t                    #Δδ_2(°)
delta_2 = delta_1+del_delta_2                #δ_2(°)
del_w_r_3 = del_w_r_2                        #Δω_r_3(°/sec)
w_r_3 = w_r_2+del_w_r_3                      #ω_r_3(°/sec)
del_delta_3 = w_r_3*del_t                    #Δδ_3(°)
delta_3 = delta_2+del_delta_3                #δ_3(°)
del_w_r_4 = del_w_r_2                        #Δω_r_4(°/sec)
w_r_4 = w_r_3+del_w_r_4                      #ω_r_4(°/sec)
del_delta_4 = w_r_4*del_t                    #Δδ_4(°)
delta_4 = delta_3+del_delta_4                #δ_4(°)
del_w_r_5 = del_w_r_2                        #Δω_r_5(°/sec)
w_r_5 = w_r_4+del_w_r_5                      #ω_r_5(°/sec)
del_delta_5 = w_r_5*del_t                    #Δδ_5(°)
delta_5 = delta_4+del_delta_5                #δ_5(°)

#Result
print('Power angle, δ_0 = %.2f° ' %delta_0)
print('Value of δ vs t are:')
print('_________________________')
print('  t(Sec)   :   δ(degree)')
print('_________________________')
print('   %.1f     :    %.2f°' %(0,delta_0))
print('   %.2f    :    %.2f°' %((del_t),delta_1))
print('   %.2f    :    %.2f°' %((del_t+del_t),delta_2))
print('   %.2f    :    %.2f°' %((del_t*3),delta_3))
print('   %.2f    :    %.2f°' %((del_t*4),delta_4))
print('   %.2f    :    %.2f°' %((del_t*5),delta_5))

Power angle, δ_0 = 14.48°
Value of δ vs t are:
_________________________
t(Sec)   :   δ(degree)
_________________________
0.0     :    14.48°
0.05    :    17.85°
0.10    :    27.98°
0.15    :    44.85°
0.20    :    68.48°
0.25    :    98.85°