CHAPTER 15: POWER SYSTEM STABILITY¶

Example 15.1, Page number 561¶

In [1]:
G = 50.0            #Rating of machine(MVA)
f = 50.0            #Frequency of turbo generator(Hz)
V = 11.0            #Voltage rating of machine(kV)
H = 9.0             #Cycle corresponding to 180 ms
P_0 = 40.0          #Pre-fault output power(MW)
delta_0 = 20.0      #Rotor angle at instant of fault(degree)

P_0_close = 0                              #Output power at instant of reclosing(MW)
P_a = P_0 - P_0_close                      #Net accelerating power(MW)
delta_sqr = P_a*180*f/(G*H)                #double derivative(elect.degrees/sec^2)

from scipy.integrate import quad

def integrand1(t):                         #Integrates the double derivative to 800*t
return delta_sqr
a, err = quad(integrand1, 0, 180*10**-3)   #Rotor velocity(electrical degrees/sec)

def integrand2(t):                         #Integrates the double derivative to 400*t^2
return delta_sqr*t
b, err = quad(integrand2, 0, 180*10**-3)
delta = delta_0 + b                        #Rotor angle(electrical degrees)

print('Rotor angle at the instant of reclosure = %.2f° electrical degrees' %delta)
print('Rotor velocity at the instant of reclosure = %.1f electrical degrees/sec' %a)

Rotor angle at the instant of reclosure = 32.96° electrical degrees
Rotor velocity at the instant of reclosure = 144.0 electrical degrees/sec


Example 15.2, Page number 571¶

In [1]:
V = 1.0         #Infinite bus voltage(p.u)
E = 1.0         #e.m.f of finite generator behind transient reactance(p.u)
X_T = 0.8       #Transfer reactance(p.u)
P_i = 0.5       #Input power(p.u)
P_i_d = 0.8     #p.u
P_0 = 0.5       #Output power(p.u)
P = 0.5         #Power(p.u)

import math
P_m = E*V/X_T                       #Amplitude of power angle curve(p.u)
delta_0 = math.asin(P_i/P_m)        #Radians
delta = math.asin(P_i_d/P_m)        #Radians
delta_m = math.pi-delta             #Radians
A_acc = P_i_d*(delta-delta_0)-P_m*(math.cos(delta_0)-math.cos(delta))   #Possible area of accceleration
A_dec = P_m*(math.cos(delta)-math.cos(delta_m))-P_i_d*(delta_m-delta)   #Possible area of deceleration

if(A_acc < A_dec):
print('System is stable')
stability = A_dec/A_acc
print('Margin of stability = %.2f' %stability)
else:
print('System is not stable')

System is stable
Margin of stability = 12.59


Example 15.3, Page number 572-575¶

In [1]:
x = 0.25        #Transient reactance(p.u)
E = 1.0         #e.m.f of finite generator behind transient reactance(p.u)
x_T = 0.1       #Reactance of transformer(p.u)
x_L = 0.4       #Reactance of one line(p.u)
P_i = 0.25      #Pre-fault power(p.u)

import math
X_T = x+x_T+(x_L/2)                 #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T                      #Amplitude of power angle curve at pre-fault state(p.u)
X_T1 = 1.45                         #Transfer reactance b/w finite generator & infinite bus at faulty state(p.u).Refer texbook problem for figure
P_m1 = E**2/X_T1                    #Amplitude of power angle curve at faulty state(p.u)
r1 = X_T/X_T1
delta_0 = math.asin(P_i/P_m)        #Radians
delta_1 = math.asin(P_i/(r1*P_m))   #Radians
delta_m = math.pi - delta_1         #Radians
from scipy.integrate import quad

def integrand1(delta):
return r1*P_m*math.sin(delta)
a, err = quad(integrand1, delta_0, delta_1)
A_acc = P_i*(delta_1-delta_0) - a

def integrand2(delta):
return r1*P_m*math.sin(delta)
b, err = quad(integrand2, delta_1, delta_m)
A_dec = b - P_i*(delta_m-delta_1)
limit = 0.5648                      #Obtained by iterations.Refer textbook.Here assigned directly.

if(A_acc < A_dec):
print('System is Stable')
stability = A_dec/A_acc
print('Margin of stability = %.2f' %stability)
else:
print('System is not stable')
print('Transient stability limit = %.4f p.u' %limit)
print('\nNOTE : ERROR : angle delta_0 = 7.9° = 0.13788 radian not 0.014 radian as in textbook')

System is Stable
Margin of stability = 38.31
Transient stability limit = 0.5648 p.u

NOTE : ERROR : angle delta_0 = 7.9° = 0.13788 radian not 0.014 radian as in textbook


Example 15.4, Page number 575-578¶

In [1]:
x = 0.25        #Transient reactance(p.u)
E = 1.0         #e.m.f of finite generator behind transient reactance(p.u)
x_T = 0.1       #Reactance of transformer(p.u)
x_L = 0.4       #Reactance of one line(p.u)
P_i = 0.7       #Pre-fault power(p.u)

import math
X_T = x+x_T+(x_L/2)                 #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T                      #Amplitude of power angle curve at pre-fault state(p.u)
X_T1 = 1.45                         #Transfer reactance b/w finite generator & infinite bus at faulty state(p.u).Refer texbook problem for figure
P_m1 = E**2/X_T1                    #Amplitude of power angle curve at faulty state(p.u)
r1 = X_T/X_T1
X_T2 = x+x_T+x_L                    #Transfer reactance for post fault state(p.u)
r2 = X_T/X_T2
P_m2 = r2*P_m
delta_0 = math.asin(P_i/P_m)        #Radians
delta_1 = math.asin(P_i/(r2*P_m))   #Radians
delta_m = math.pi - delta_1         #Radians
delta_c = 0.7                       #Specified value(radians)
from scipy.integrate import quad

def integrand1(delta):
return r1*P_m*math.sin(delta)
a, err = quad(integrand1, delta_0, delta_c)
A_acc = P_i*(delta_c-delta_0) - a

def integrand2(delta):
return r2*P_m*math.sin(delta)
b, err = quad(integrand2, delta_c, delta_m)
A_dec = b - P_i*(delta_m-delta_c)
cos_delta_cr = ((delta_m-delta_0)*math.sin(delta_0)-r1*math.cos(delta_0)+r2*math.cos(delta_m))/(r2-r1)
delta_cr = math.acos(cos_delta_cr)*180/math.pi

if(A_acc < A_dec):
print('System is Stable')
stability = A_dec/A_acc
print('Margin of stability , K = %.2f' %stability)
else:
print('System is not stable')
print('Critical clearing angle for a certain pre-fault power = %.2f°' %delta_cr)
print('Critical clearing time will be known from circuit-breaker specifications')

System is Stable
Margin of stability , K = 7.98
Critical clearing angle for a certain pre-fault power = 111.48°
Critical clearing time will be known from circuit-breaker specifications


Example 15.5, Page number 578-580¶

In [1]:
P_i = 0.75      #Pre-fault power(p.u)
f = 50.0        #Frequency(Hz)
H = 6.0         #Value of H for finite machine(sec)
x_G = 0.2       #Reactance of machine(p.u)
x_T = 0.1       #Reactance of transformer(p.u)
x_L = 0.4       #Reactance of line(p.u)
V = 1.0         #Voltage of infinite bus(p.u)
E = 1.0         #e.m.f of finite generator behind transient reactance(p.u)

import math
X_T = x_G+x_T+(x_L)                           #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T                                #Amplitude of power angle curve at pre-fault state(p.u)
delta_0 = math.asin(P_i/P_m)                  #Radians
delta_0a = delta_0*180/math.pi
delta_cr = math.acos((math.pi-2*delta_0)*math.sin(delta_0)-math.cos(delta_0))
delta_cra = delta_cr*180/math.pi
t_cr = ((delta_cra-delta_0a)*2*H/(180*f*P_i))**0.5

print('Critical clearing angle for circuit breaker at bus 1 = %.2f°' %delta_cra)
print('Time for circuit breaker at bus 1 ,t_cr = %.3f sec' %t_cr)

Critical clearing angle for circuit breaker at bus 1 = 77.42°
Time for circuit breaker at bus 1 ,t_cr = 0.285 sec