CHAPTER – 11: Power System Stability 11.1
Assume that three machines whose ratings and inertia constants are respectively
given by G1 , H 1 , G2 , H 2 , and G3 , H 3 are operating in synchronism, that is G1 H 1 d
2
δ 1
G s π f dt 2 G2 H 2 d
2
δ 2 2
G s π f
dt
G3 H 3 d
2
δ 3
G s π f dt 2
= (Pi1 − Pu1 )
G1 Gs
= (Pi 2 − Pu 2 ) = (Pi 3 − Pu 3 )
G2 Gs
G3 Gs
Since the machines are swinging coherently,δ 1 = δ 2 = δ 3 = δ . Thus ,
⎛ G1 H 1 G2 H 2 G3 H 3 ⎞ 1 d 2δ G1 ⎜⎜ ⎟⎟ P P P P P P + + = + + − + + ( ) ( ) 1 2 3 1 1 1 i i i u u u Gs G s ⎠ π f dt 2 Gs ⎝ G s H eq d 2δ π f
2
dt
= Pieq ( pu ) − Pueq ( pu )
⎛ G H G H G H ⎞ where H eq = ⎜⎜ 1 1 + 2 2 + 3 3 ⎟⎟ Gs G s ⎠ ⎝ G s ⎛ G H G H G H ⎞ H eq = ⎜⎜ 1 1 + 2 2 + 3 3 ⎟⎟ Gs G s ⎠ ⎝ G s » 200*4/1000+500*3/1000+750*5/1000 = 6.05 MVJ/MVA M=
11.2
GH
180 × f
= 1000*6.05/(180*50) = 0.6722 MJ-sec/elect-deg
Synchronous speed ω = 2π ×
2 × 50 2
rads/sec = 314.16 rads/sec
KE = 1/2*75000*(314.1593)^2*10^(-6) = 3.7011e+003 MJ H = (3.7011e+03)/250 = 14.8044 MJ/MVA M = 250*14.8044/(180*50) = 0.4112 MJ-sec/elec. deg 11.3
» M=175/(180*50) = 0.0194 MJ-sec/rad At the time of fault, power input is given by » Pi=50*0.8 = 40 MW Power output during fault » Pu=0.4*40 = 16 MW Accelerating power is as follows » Acclpower=Pi-Pu = 24 MW Acceleration is given by 2
» Accl=Acclpower/M = 1.2343e+003 elec deg/sec 11.4
Pre-fault current » Ig=Pi*10^6/(sqrt(3)*11*10^3*0.8) = 874.7731 A Ig=Ig*(0.8-i*0.6) = 6.9982e+002 -5.2486e+002i A Generator voltage per phase » Vgpp=(11*10^3/sqrt(3)+i*Ig*2)/(10^3) = 7.4006 + 1.3996i = 7.5318 kV Line to line generator voltage » Vgll=abs(Vgpp)*sqrt(3) = 13.0454 kV Pre-fault rotor angle 0
» delta=asin(Pi/(Vgll*11))*180/pi = 16.1854
2
» Accl=Accl*pi/180 = 21.5423 elec. rads/sec » t=5/50 = 0.1000 sec
Rotor velocity at the time of occurrence of fault » ddeltadt=sqrt(2*Accl*(delta*pi/180)) = 3.4887 elec. rads/sec = 66.6290 rpm Rotor velocity at the end of acceleration period =
120 × 50 2
+ 66.6290
= 1566.6290 rpm Rotor angle at the end of the acceleration period 0
» deltan=Accl*t^2+delta*pi/180 = 0. 0.4979 rads = 28.5283 11.5 The swing equations of the two machines are H 1 d
2
δ 1
= P1i − P1u
2
π f
dt
H 2 d π f
2
δ 2
= P2i − P2u
2
dt
⎛ P1i − P1u ⎞ ⎛ P2i − P2u ⎞ ⎟ ⎟⎟ = − π f ⎜ π f ⎜ ⎜ ⎟ ⎜ 2 2 H H dt dt 1 2 ⎝ ⎠ ⎝ ⎠ ⎛ H P − H 1 P2i H 2 P1u − H 1 P2u ⎞ ⎟⎟ = π f ⎜⎜ 2 1i − H H H H 1 2 1 2 ⎝ ⎠ 2 1 H 1 H 2 ⎛ d (δ 1 − δ 2 ) ⎞ H 2 P1i − H 1 P2i H 2 P1u − H 1 P2u ⎜ ⎟= − 2 ⎜ ⎟ H H H 1 + H 2 + π f H 1 + H 2 ⎝ dt 1 2 ⎠ d
2
1 π f
δ 1
−
H eq
d
d
2
2
δ 2
δ 12
where H eq =
H eq =
Peqi =
H 1 H 2 H 1 + H 2
H 1 H 2 H 1 + H 2
=
H 1 + H 2
, δ 12 = (δ 1 − δ 2 ), Peqi =
4.5 × 6.0 10.5
H 2 P1i − H 1 P2i
Max power =
11.6
= Peqi − Pequ
2
dt
=
H 1 + H 2
= 2.5714 MJ / MVA
6.0 × 1.25 − 4.5 × 1.1 10.5
1.2 × 1.15
(0.15 + 0.4 + 0.15)
= 0.2429 pu
= 1.9714 pu
Input data: » Pl=0.8;E=1.0;pf=0.85;H=4.0; 0
» phi=-acos(pf) = -0.5548 rads. = -31.7883 Load current:
H 2 P1i − H 1 P2i
, andPequ =
H 2 P1u − H 1 P2u H 1 + H 2
» I=Pl/(E*pf) = 0.9412 pu = 0.8000 - 0.4958i pu Pre-fault transfer reactance between the internal generator bus an d infinite bus » Xl=0.25+0.15 = 0.4000 pu Generator voltage behind transient reactance » Edash=E+i*I*Xl = 1.1983 + 0.3200i pu = 1.2403 pu Power angle 0
» delta=angle(Edash) = 0.2610 rads. = 14.9514
During fault transfer reactance between generator internal bus and infinite bus » Xlf=1.15; Power output during fault » Poutf=(abs(Edash)*E/Xlf)*sin(delta) = 0.2783 pu Accelerating power at the time of fault » Paccl=0.8-Poutf = 0.5217 pu Inertia constant » M=H/(pi*50) = 0.0255 MJ-sec/elec. rad. Acceleration at the time of fault 2
» Accl=Paccl/M = 20.4886 elec. rads./sec 2
d
δ 2
dt
=
1.2403 * 1.0 ⎛ ⎞ sin δ ⎟ = 31.3725 − 42.2941sin δ ⎜ 0.8 − 0.0255 ⎝ 1.15 ⎠ 1
' 11.7 From the previous problem E =1.2403 pu, M = 0.0255 MJ-sec/elec. rad.,
δ0 =
0
0.2610 rads. = 14.9514 . Power input Pi 0 = 0.8 pu Fault duration t = 10 cycles = 0.2 secs. Rotor angle δ at the time of fault clearance is given by δ =
Pi 0
t + δ 0 = 0.8885 rads. 2
2 M
Accelerating area A1 = Pi0 (δ – δ0) = 0.5020 1.2403 × 1.0
Decelerating area A2 =
0.4
[cos(δ ) − cos δ m ] − 0.8(δ m − δ )
= 3.1007[cos(δ ) − cos δ m ] − 0.8(δ m − δ ) For stability, A1 = A2 , that is, 3.1007[cos(δ ) − cos δ m ] − 0.8(δ m − δ ) = 0.5020 f (δm) = 3.1007 cos δm + 0.8 δm – 2.1641 df (δ m ) d δ m
= 3.1007 sin δ m
» tolr=0.0001; » [deltam] = P1207(tolr); initial estimate of deltam60*pi/180 0
deltam = 66.5220
0
Maximum swing of the rotor angle is δm = 66.5220 . Since δm is less than ( π – δ0) the system will remain stable. (ii) Assume that the critical clearing angle is δc. Accelerating area A1 = Pi0 (δc – δ0) = 0.8 (δc – 0.2610 ) = 0.8 δc – 0.2088 Decelerating area A2 =
δ m
∫
δ c
3.1007 sin δ d δ − Pi 0 (δ m − δ c )
A2 = 3.1007[cos δ c + cos δ 0 ] − 0.8(π − δ 0 − δ c )
= 3.1007 cos δc + 0.8 δc +0.6912 For stability A1 = A2, that is, 0.8 δc – 0.2088 = 3.1007 cos δc + 0.8 δc +0.6912 δ c
⎛ (0.2088 + 0.6912 ) ⎞ 0 = cos −1 ⎜ − ⎟ = 1.8653 rads. = 106.8733 3.1007 ⎝ ⎠
t c =
11.8
2 M (δ c − δ 0 ) Pi 0
=
2 × 0.0255(1.8653 − 0.2610 ) 0.8
= 0.3198 secs. = 16 cycles
Input data: » Pi0=0.8;Einfint=1.0;Xlprefault=0.4;Xlfault=1.05;Xlpfault=0.55;pf=1.0; Load current » I=0.8/(Einfint*pf) = 0.8000 pu Voltage behind generator transient reactance » Edash=Einfint+i*Xlprefault*I = 1.0000 + 0.3200i pu = 1.05 pu Initial power angle 0
» delta0=angle(Edash) = 0.3097 rads. = 17.7447 Pre fault maximum power » PmA=Edash*Einfint/Xlprefault = 2.6249 pu Post fault maximum power » PmB=Edash*Einfint/Xlpfault = 1.9090 pu During fault maximum power » PmC=Edash*Einfint/Xlfault = 1.0000 pu Maximum power angle
0
» deltam=(180-(asin(Pi0/PmB)*180/pi)) = 155.2243 = 2.7092 rads
Cosine of critical clearing power angle » cosdeltac=(Pi0*(deltam-delta0)+PmB*cos(deltam)-PmC*cos(delta0))/(PmBPmC) = -0.8427 Critical clearing power angle 0
» deltac=acos(cosdeltac)*180/pi = 147.4307 Critical clearing time computation:
» tc=sqrt(2*M*(147.4307-17.7447)/0.8) = 0.3796 secs. = 18.9800 cycles. Plotting the power output versus power angle » delta=linspace(0,pi,500); » P1=2.6249*sin(delta); » P2=1.0*sin(delta); » P3=1.9090*sin(delta); » plot(delta*180/pi,P1,'k-') » grid » hold on » plot(delta*180/pi,P2,'k-.') » plot(delta*180/pi,P3,'k--') » x=linspace(0,180,500); » y=0.8; » plot(x,y,'k-') » xlabel('Power angle in degrees') » ylabel('Power output in pu')
» legend('Pre-fault power output','During fault power output','Post fault power output')
3 Pre-fault power output During fault power output Post fault power output
2.5
2
u p n i t u p t u 1.5 o r e w o P
1
0.5
0
11.9
0
20
40
60 80 100 120 Power angle in degrees
140
160
180
Input data: » Pi0=1.0;Edash=1.20;Einfinit=1.0;Xlprefault=0.6;Xlfault=1.0267; » tf=15/50;H=3.5; (i) Rotor angle prior to the fault 0
» delta0=asin(Pi0*Xlprefault/(Edash*Einfinit)) = 0.5236 rads. = 30.0000
(ii) Computation of generator output, accelerating power, a nd rotor acceleration » Poutgen=(Edash*Einfinit/Xlfault)*sin(delta0) = 0.5844 pu » Paccl=Pi0-Poutgen = 0.4156 pu » M=H/(pi*50) = 0.0223 MJ-sec./elec. rad. 2
» Accl=Paccl/M = 18.6522 rads./sec
Computation of rotor angle and decelerating power at the instance of fault clearance:
0
» delta=delta0+(Paccl/(2*M))*tf^2 = 1.3629 rads. = 78.0912 Generator output at the instance of fault clearance » Poutgenf=(Edash*Einfinit/Xlfault)*sin(delta) = 1.1436 pu
Genrator accelerating power at the instance of fault clearance » Paccl=Pi0-Poutgenf = -0.1436 pu
11.10 f unct i on ydot = p1210( t , y) ; ydot =t ^3- y^2;
Input to the command window » tspan=[0 1]; » y0=0.5; » [t,y]=ode23('p1210',tspan,y0) Output t=
y=
0
0.5000
0.1000
0.4762
0.2000
0.4549
0.3000
0.4367
0.4000
0.4226
0.5000
0.4144
0.6000
0.4141
0.7000
0.4242
0.8000
0.4477
0.9000
0.4875
1.0000
0.5469
Plot of y vs i » plot(t,y,'k--'); » hold on » grid » xlabel('t') » ylabel('y') » title('Solution of differential eq. of problem 12.10') » hold off »
Solution of differential eq. of problem 12.10 0.56 0.54 0.52 0.5 y 0.48
0.46 0.44 0.42 0.4
0
0.1
0.2
0.3
0.4
0.5 t
0.6
0.7
0.8
0.9
1
11.11 The second order differential swing equation is written as two first order
differential equations, that is, .
Z 1 = δ and Z 2 = δ
⎡ Z 1 ⎤ ⎡δ ⎤ Z = ⎢ ⎥ = ⎢ . ⎥, ⎣ Z 2 ⎦ ⎢⎣δ ⎥⎦
⎡ . ⎤ ⎡ . ⎤ . ⎢ Z 1 ⎥ ⎡ Z 2 ⎤ ⎢ Z 2 ⎥ Z = ⎢ . ⎥ = ⎢ . ⎥ = ⎢ ⎥ ⎢ Z 2 ⎥ ⎢⎣ Z 2 ⎥⎦ ⎣⎢ Pa ⎦⎥ ⎣ ⎦
Function ‘p1211’ comput es t he st at e vari abl es f unct i on Zdot = p1211( t , Z) ; % I nput i s t i me t and Z = [ Z( 1) ; Z( 2) ] % Out put i s Zdot i f t >= 0. 25; Pa=0. 75- 1. 5*si n( Z( 1) ) ; el se Pa=0. 75- 0. 30*si n( Z( 1) ) ; end Zdot=[ Z( 2) ; Pa] ;
The f ol l owi ng sc r i pt f i l e execut es
» tspan=[0 2.0];
ode23
% Sets the time span
» Z0=[30*pi/180;0]; % Assigns the initial rotor angle in radians » [t,Z]=ode23('P1211',tspan,Z0)
% Executes ode23
Output: .
t vs δ and t vs
δ is
obtained and plotted
» x=Z(:,1);y=Z(:,2);
% Sets first and second columns of Z equal to x and y respectively
» subplot(2,1,1) » plot(t,x*180/pi,'k--') » grid » xlabel('Time t in sec.') » ylabel('Rotor angle in radians') » title('Rotor Angle Versus Time') » subplot(2,1,2) » plot(t,y,'k--') » grid » xlabel('Time t in sec.') » ylabel('Rotor velocity in radians/sec.') » title('Rotor Velocity Versus Time') »
Rotor Angle Versus Time 38 s n a 36 i d a r n 34 i e l g 32 n a r o 30 t o R
28
0
0.2
0.4
0.6
. c e 0.15 s / s n 0.1 a i d a r 0.05 n i y t 0 i c o l e-0.05 v r o t o -0.1 R 0
0.2
0.4
0.6
0.8 1 1.2 1.4 Time t in sec. Rotor Angle Versus Time
0.8 1 1.2 Time t in sec.
1.4
1.6
1.8
2
1.6
1.8
2
11.12 For derivation please see 7.3.1 Power flow equations in [Y bus] frame of reference.
The power flow equation is given by P = E 1 Y 11 cos θ 11 + E 1 E 2Y 12 cos(θ 12 − δ 12 ) where δ12 = δ1 – δ2 2
(I)
Eq. (I) modifies as follows for the problem case P = E 1 Y 12 cos θ 12 + E 1 E 2Y 12 cos(θ 12 − δ 12 ) 2
For maximum power transfer;
(θ 12 − δ 12 ) = 0. Hence, δ12 = θ12. θ12 of a transmission line
0
having resistance is < 90 . If Y 12 =
1 Z
and cos θ 12 =
R Z
. Assuming Z = R + j X .
The maximum steady limit is written as P=−
E 12 R Z 2
+
E 1 E 2 Z
(II)
[Note: In deriving (I) in Sec.7.3.1, θ12 was substituted by (180 – θ12)] 0
When resistance is neglected, θ11 = 90 . Therefore power transferred is given by P =
E 1 E 2 X
Pmax =
sin δ 12
E 1 E 2 X
0
when δ12 = 90
(III)
Comparison of (II) and (III) shows that the steady limit is more when resistance of the line is neglected. 11.14 Computation of reactance between generator terminal and infinite bus
» X12=0.4/2+0.4 = 0.6000 pu Pf angle: 0
» phi=-acos(0.85) = -0.5548 rads. = -31.7883 Load current
» I=0.75/(1.0*0.85) = 0.8824 pu = 0.7500 - 0.4648i pu Voltage behind generator reactance » E2=1.0+i*X12*I = 1.2789 + 0.4500i pu Magnitude of voltage behind generator reactance » Emag2=abs(E2) = 1.3557 pu Power angle δ0 prior to the disturbance » delta0 = angle(E2) = 0.3383 rads. = 19.38540 Synchronizing coefficient » Ps=(Emag2*1.0/X12)*cos(delta0) = 2.1315 pu Coefficient of damping power » Pd=0.15;
» M=8/(pi*50) = 0.0509 MJ-sec/elec. rad Damping factor » Pd/M = 2.9452 » Ps/M = 41.8514 Natural frequency of oscillations
Ps M
= 41. 8514 = 6.4693 rads/sec.
The equation representing the motion of the rotor following a disturbance is, thus, written as ..
.
Δ δ + 2.9452Δ δ + 41.8514 = 0 The second order differential equation is written as two first order differential equations for employing MATLAB to plot variation of rotor angle δ and rotor velocity ω against time.
⎡ x1 ⎤ ⎡Δδ ⎤ x = ⎢ ⎥ = ⎢ . ⎥ ⎣ x 2 ⎦ ⎢⎣Δ δ ⎥⎦ The two first order differential equations become .
x1 = x2 .
x 2 = −2.9452 x2 − 41.8514 x1
The MATLAB function ‘p1214’ computes the state variables f unct i on xdot = p1214( t , x) ; xdot=[ x( 2) ; - 2. 9452*x( 2) - 41. 8514*x( 1) ] ;
The following script file executes ode23 » tspan=[0 5.0]; » x0=[12.5*pi/180;0]; » [t,x]=ode23('P1214',tspan,x0)
» x1=x(:,1);y1=x(:,2); » plot(t,x1*180/pi,'k-') » grid » hold on » plot(t,y1,'k-.') » xlabel('Time t in secs.') » ylabel('Rotor angle in degrees') » ylabel('Rotor velocity in rads./sec') » ylabel('Rotor angle in degrees & velocity in rads./sec') » legend('Rotor angle delta','Rotor velocity') » title('Plot of rotor angle and velocity versus time') » Plot of rotor angle and velocity versus time 14 Rotor angle delta Rotor velocity
12
c e s / . s 10 d a r n 8 i y t i c o 6 l e v & 4 s e e r g 2 e d n i e 0 l g n a r -2 o t o R
-4 -6
0
0.5
1
1.5
2 2.5 3 Time t in secs.
3.5
4
4.5
5
