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)
#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)
V_L = 1.0 #Load 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)
#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)
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)
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))
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')
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 #δ(°)
P_0 = (V_R/Z**2)*(V_S*Z*math.cos((delta_0-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_1 = 30.0 #δ(°)
P_1 = (V_R/Z**2)*(V_S*Z*math.cos((delta_1-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_2 = 60.0 #δ(°)
P_2 = (V_R/Z**2)*(V_S*Z*math.cos((delta_2-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_3 = beta #δ(°)
P_3 = (V_R/Z**2)*(V_S*Z*math.cos((delta_3-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_4 = 90.0 #δ(°)
P_4 = (V_R/Z**2)*(V_S*Z*math.cos((delta_4-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_5 = 120.0 #δ(°)
P_5 = (V_R/Z**2)*(V_S*Z*math.cos((delta_5-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
delta_6 = (math.acos(R/Z)*180/math.pi)+beta #δ(°)
P_6 = (V_R/Z**2)*(V_S*Z*math.cos((delta_6-beta)*math.pi/180)-V_R*R)/10**6 #Power received(MW/phase)
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))
#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)
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)
#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)
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)
acceleration_1 = acceleration*math.pi/180.0 #Acceleration(rad/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')
import math
#Variable declaration
f = 50.0 #Frequency(Hz)
P = 4.0 #Number of poles
alpha = 200.0 #Acceleration(°/sec^2)
alpha_rad = 3.49 #Acceleration(rad/sec^2)
n = 10.0 #Number of cycle
#Calculation
t = 1/f*n #Time(sec)
delta_rel = ((alpha_rad*2)**0.5*0.5)**2 #Relation of change in rotor angle with time(rad)
delta = delta_rel*t**2 #Change in torque angle(rad)
delta_deg = delta*180/math.pi #Change in torque angle in that period(°)
rpm_rad = (alpha_rad*2*delta)**0.5 #r.p.m(rad/sec)
rpm = rpm_rad*60.0/(math.pi*P) #r.p.m
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)
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)
#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)
#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)
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)
#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)
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)
M = G*H/(math.pi*f) #Angular momentum(MJ-sec/rad)
acceleration = P_a/M #Accelerating power(rad/sec^2)
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)
#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)
#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)
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)
w = 2.0*math.pi*N_s/60 #ω(rad/sec)
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)
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')
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 = math.asin(P_m/P_max) #δ_0(radians)
delta_0_degree = delta_0*180/math.pi #δ_0(°)
delta_m = math.pi-math.asin(P_m/(gamma_2*P_max)) #δ_1(radians)
delta_m_degree = delta_m*180/math.pi #δ_1(°)
delta_c = math.acos((P_m/P_max*(delta_m-delta_0)+gamma_2*math.cos(delta_m)-gamma_1*math.cos(delta_0))/(gamma_2-gamma_1)) #Clearing angle(radians)
delta_c_degree = delta_c*180/math.pi #Clearing angle(°)
#Result
print('Critical clearing angle, δ_c = %.2f° ' %delta_c_degree)
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 = math.asin(P_m/P_max) #δ_0(radians)
delta_0_degree = delta_0*180/math.pi #δ_0(°)
delta_f = math.pi-math.asin(P_m/(gamma_2_P_max)) #δ_1(radians)
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 = math.acos(1.0/(gamma_2-gamma_1)*((delta_f-delta_0)*math.sin(delta_0)+(gamma_2*math.cos(delta_f)-gamma_1*math.cos(delta_0)))) #Clearing angle(radians)
delta_c_degree = delta_c*180/math.pi #Clearing angle(°)
#Result
print('Critical angle, δ_c = %.2f° ' %delta_c_degree)
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 = math.asin(sin_delta_0) #δ_0(radians)
delta_0_degree = delta_0*180/math.pi #δ_0(°)
gamma_1 = 1.0/x #γ_1
delta_m = math.pi-math.asin(sin_delta_0/(gamma_2)) #δ_m(radians)
delta_m_degree = delta_m*180/math.pi #δ_m(°)
delta_c = math.acos(1.0/(gamma_2-gamma_1)*((delta_m-delta_0)*math.sin(delta_0)+(gamma_2*math.cos(delta_m)-gamma_1*math.cos(delta_0)))) #Clearing angle(radians)
delta_c_degree = delta_c*180/math.pi #Clearing angle(°)
#Result
print('Critical clearing angle, δ_c = %.f° ' %delta_c_degree)
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))