Syllabus 8EE5 COMPUTER BASED POWER SYSTEM LAB 1. Fault analysis and verify the results using MATLAB or any available software for the cases: (i) LG Fault (ii) LLG Fault (iii) LL Fault and (iv) 3 -Phase Fault 2. Load flow analysis for a given system (for 3 to 6 bus) using (i) Gauss Seidal (ii) Newton Raphson (iii) Fast Decoupled Method and verify results using MATLAB or any available software 3. Study of voltage security analysis
4. Study of economic load dispatch problem with different methods. 5. Study of transient stability analysis using MATLAB/ETAP Software.
EXPERIMENT NO. 1
Object :Fault analysis (for 3 to 6 bus) and verify the results using MATLAB or any available software for the cases: (i) LG Fault (ii) LL Fault (iii) LLG Fault and (iv) 3-Phase Fault
Software Tool: MATLAB Theory: (i)
Single-Line-to-Ground Fault
Let a 1LG fault has occurred at node k of a network. The faulted segment is then as shown in Fig.1 where it is assumed that phase-a has touched the ground through an impedance Zf. Since the system is unloaded before the occurrence of the fault we have (1)
Fig. 1 Representation of 1LG fault Also the phase-a voltage at the fault point is given by (2) From (1) we can write
(3) Solving (3) we get
(4) This implies that the three sequence currents are in series for the 1LG fault. Let us denote the zero, positive and negative sequence Thevenin impedance at the faulted point as Z kk0 ,Z kk1 and Z kk2respectively. Also since the Thevenin voltage at the faulted phase is Vfwe get three sequence circuits that are similar to the ones shown in Fig. 1 We can then write
(5) Then from (4) and (5) we can write
(6)
(7) The Thevenin equivalent of the sequence network is shown in Fig. 8.3.
Fig. 2 Thevenin equivalent of a 1LG fault.
(ii)
Line-to-Line Fault
The faulted segment for an L-L fault is shown in Fig. 3 where it is assumed that the fault has occurred at node k of the network. In this the phases b and c got shorted through the impedance Zf. Since the system is unloaded before the occurrence of the fault we have (8)
Fig. 3Representation of L-L fault. (9)
(10)
(11) Therefore no zero sequence current is injected into the network at bus k and hence the zero sequence remains a dead network for an L-L fault. The positive and negative sequence currents are negative of each other. Now from Fig. 3 we get the following expression for the voltage at the faulted point (12)
(13) Moreover since I fa0 = Ifb0 = 0 and Ifa1 = - I fb2, we can write (14)
Therefore combining (12) - (14) we get
(15)
Equations (12) and (15) indicate that the positive and negative sequence networks are in parallel.
(16) The sequence network is then as shown in Fig.4 . From this network we get
Fig. 4 Thevenin equivalent of an LL fault
(iii)
Double- Line -to Ground Fault
The faulted segment for a 2LG fault is shown in Fig. 5 where it is assumed that the fault has occurred at node k of the network. In this the phases b and c got shorted through the impedance Zfto the ground. Since the system is unloaded before the occurrence of the fault
(17) we have the same condition as (8) for the phase-a current. Therefore
Fig. 5 Representation of 2LG fault.
(18) Also voltages of phases b and c are given by
(19)
Therefore
(20)
We thus get the following two equations from (19)
(21)
Substituting (18) and (20) in (21) and rearranging we get
(22)
(23)
Also since I fa= 0 we have
(24) The Thevenin equivalent circuit for 2LG fault is shown in Fig. 6. From this figure we get
(25) (26) The zero and negative sequence currents can be obtained using the current divider principle as
Fig. 6 Thevenin equivalent of a 2LG fault.
Code:
% Thismatlab code for fault calculation of 4 bus system % L-L fault L-G fault L-L-G fault % input quantity is bus voltage and % Zero positive and negative sequence components of % all calculation in per unit alpha=cos(2*pi/3)+i*sin(2*pi/3); z1=input('enter the value of positive sequence impedance ') z2=input('enter the value of negative sequence impedance ') z0=input('enter the value of zero sequence impedance ') zf=input('input the ground impedance if system is solidly grounded zf=0') v=input('Enter the value of thevenin in voltage in P.u. '); k=input('enter the type fault of fault if L-L press 1 fault L-G fault press 2 L-L-G fault press 3 for press 4 for L-L-L'); if (k==1) i_LL_fault=3*v/(z1+z2+z0+3*zf); elseif k==2 i_LG_fault=v/(z1+z2+zf); i_LG_fault=(alpha^2-alpha)*v/(z1+z2+zf) else if k==3 disp('For double line to ground fault current'); zg=input('the value of zg ground impedance') ia1=v/((z1+zf+(z2+zf)*(z0+zf+3*zg))/(z2+z0+2*zf+3*zg)) ia2=-ia1*(z0+z2+2*zf+3*zg)/(z2+z0+2*zf+3*zg) ia0=ia1+ia2
else phase_fault=v/z1 end end Output: enter the value of positive sequence impedance .3 z1 = 0.3000 enter the value of negative sequence impedance .3 z2 = 0.3000 enter the value of zero sequence impedance .1 z0 = 0.1000 input the ground impedance if system is solidly grounded zf=00 zf = 0 Enter the value of thevenin in voltage in P.u. 1 enter the type fault of fault if L-L press 1 fault L-G fault press 2 L-L-G fault press 3 for press 4 for L-L-L4 phase_fault = 3.3333
Experiment No. 2 Object :Load flow analysis for a given system (for 3 to 6 bus) using (i) Gauss Seidal (ii) Newton Raphson (iii) Fast Decoupled Method and verify results using MATLAB or any available software
Software Tool: MATLAB (I) GAUSS SIEDEL METHOD Exercise: (a) Using the Gauss-Sedel method determine the phasor values of the voltage at the load buses 2 and 3 (P-Q) accurat to four decimal places. (b) Find the slack bus real and reactive power (c) Determine the line flow and line losses. Construct a power flow diagram showing the direction of line flow
Code: y12=10-j*20; y13=10-j*30; y23=16-j*32; V1=1.05+j*0; iter =0; S2=-2.566-j*1.102; S3=-1.386-j*.452;
V2=1+j*0; V3=1+j*0; for I=1:10; iter=iter+1; V2 = (conj(S2)/conj(V2)+y12*V1+y23*V3)/(y12+y23); V3 = (conj(S3)/conj(V3)+y13*V1+y23*V2)/(y13+y23); disp([iter, V2, V3]) end V2= .98-j*.06; V3= 1-j*.05; I12=y12*(V1-V2); I21=-I12; I13=y13*(V1-V3); I31=-I13; I23=y23*(V2-V3); I32=-I23; S12=V1*conj(I12); S21=V2*conj(I21); S13=V1*conj(I13); S31=V3*conj(I31); S23=V2*conj(I23); S32=V3*conj(I32); I1221=[I12,I21] I1331=[I13,I31] I2332=[I23,I32] S1221=[S12, S21 (S12+S13) S12+S21] S1331=[S13, S31 (S31+S32) S13+S31] S2332=[S23, S32 (S23+S21) S23+S32]
Output: 1.0000
0.9825 - 0.0310i 1.0011 - 0.0353i
2.0000
0.9816 - 0.0520i 1.0008 - 0.0459i
3.0000
0.9808 - 0.0578i 1.0004 - 0.0488i
4.0000
0.9803 - 0.0594i 1.0002 - 0.0497i
5.0000
0.9801 - 0.0598i 1.0001 - 0.0499i
6.0000
0.9801 - 0.0599i 1.0000 - 0.0500i
7.0000
0.9800 - 0.0600i 1.0000 - 0.0500i
8.0000
0.9800 - 0.0600i 1.0000 - 0.0500i
9.0000
0.9800 - 0.0600i 1.0000 - 0.0500i
10.0000
0.9800 - 0.0600i 1.0000 - 0.0500i
I1221 = 1.9000 - 0.8000i -1.9000 + 0.8000i I1331 = 2.0000 - 1.0000i -2.0000 + 1.0000i I2332 = -0.6400 + 0.4800i 0.6400 - 0.4800i S1221 = Columns 1 through 3 1.9950 + 0.8400i -1.9100 - 0.6700i 4.0950 + 1.8900i Column 4 0.0850 + 0.1700i
S1331 = Columns 1 through 3 2.1000 + 1.0500i -2.0500 - 0.9000i -1.3860 - 0.4520i Column 4 0.0500 + 0.1500i S2332 = Columns 1 through 3 -0.6560 - 0.4320i 0.6640 + 0.4480i -2.5660 - 1.1020i Column 4 0.0080 + 0.0160i
(i)
Newton -Raphson Method :
Exercise: Figure shows the one line diagram of a simple three bus system with generators at buses 1 and 3. The magnitude of voltage at bus 1 is adjusted to 1.05 pu. Voltage magnitude at bus 3 is fixed at 1.04 pu with a real power generation of 200 MW. A load consisting of 400MW and 250 Mvar is taken from bus 2. Line impedance are marked in per unit on a base of 100 MVA and the line charging susceptance are neglected. Obtain the power flow solution by the NW methods.
Code : iter=0; pwracur = 0.00025; % Power accuracy DC = 10;
% Set the maximum power residual to a high value
while max(abs(DC)) > pwracur iter = iter +1 P=[V(2)*V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1))+V(2)^2*Y(2,2)*cos(t(2,2))+ ... V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3)); V(3)*V(1)*Y(3,1)*cos(t(3,1)-d(3)+d(1))+V(3)^2*Y(3,3)*cos(t(3,3))+ ... V(3)*V(2)*Y(3,2)*cos(t(3,2)-d(3)+d(2))];
Q= -V(2)*V(1)*Y(2,1)*sin(t(2,1)-d(2)+d(1))-V(2)^2*Y(2,2)*sin(t(2,2))- ... V(2)*V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3)); J(1,1)=V(2)*V(1)*Y(2,1)*sin(t(2,1)-d(2)+d(1))+... V(2)*V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3)); J(1,2)=-V(2)*V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3)); J(1,3)=V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1))+2*V(2)*Y(2,2)*cos(t(2,2))+... V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3)); J(2,1)=-V(3)*V(2)*Y(3,2)*sin(t(3,2)-d(3)+d(2)); J(2,2)=V(3)*V(1)*Y(3,1)*sin(t(3,1)-d(3)+d(1))+... V(3)*V(2)*Y(3,2)*sin(t(3,2)-d(3)+d(2)); J(2,3)=V(3)*Y(2,3)*cos(t(3,2)-d(3)+d(2)); J(3,1)=V(2)*V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1))+... V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3)); J(3,2)=-V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3)); J(3,3)=-V(1)*Y(2,1)*sin(t(2,1)-d(2)+d(1))-2*V(2)*Y(2,2)*sin(t(2,2))-... V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3)); DP = Ps - P; DQ = Qs - Q; DC = [DP; DQ] J DX = J\DC d(2) =d(2)+DX(1); d(3)=d(3) +DX(2);
V(2)= V(2)+DX(3); V, d, delta =180/pi*d; end P1= V(1)^2*Y(1,1)*cos(t(1,1))+V(1)*V(2)*Y(1,2)*cos(t(1,2)-d(1)+d(2))+... V(1)*V(3)*Y(1,3)*cos(t(1,3)-d(1)+d(3)) Q1=-V(1)^2*Y(1,1)*sin(t(1,1))-V(1)*V(2)*Y(1,2)*sin(t(1,2)-d(1)+d(2))-... V(1)*V(3)*Y(1,3)*sin(t(1,3)-d(1)+d(3)) Q3=-V(3)*V(1)*Y(3,1)*sin(t(3,1)-d(3)+d(1))-V(3)*V(2)*Y(3,2)*... sin(t(3,2)-d(3)+d(2))-V(3)^2*Y(3,3)*sin(t(3,3)) ter = 1 DC = -2.8600 1.4384 -0.2200 J= 54.2800 -33.2800 24.8600 -33.2800 66.0400 -16.6400 -27.1400 16.6400 49.7200 DX = -0.0453 -0.0077 -0.0265
V= 1.0500 0.9735 1.0400 d= 0 -0.0453 -0.0077 iter = 2 DC = -0.0992 0.0217 -0.0509 J= 51.7247 -31.7656 21.3026 -32.9816 65.6564 -15.3791 -28.5386 17.4028 48.1036 DX = -0.0018 -0.0010 -0.0018
V= 1.0500 0.9717 1.0400 d= 0 -0.0471 -0.0087 iter = 3 DC = 1.0e-003 * -0.2166 0.0382 -0.1430 J= 51.5967 -31.6939 21.1474 -32.9339 65.5976 -15.3516 -28.5482 17.3969 47.9549 DX = 1.0e-005 * -0.3856 -0.2386
-0.4412 V= 1.0500 0.9717 1.0400 d= 0 -0.0471 -0.0087 P1 = 2.1842 Q1 = 1.4085 Q3 = 1.4618 iter = 1 DC = -2.8600 1.4384 -0.2200 J= 54.2800 -33.2800 24.8600 -33.2800 66.0400 -16.6400
-27.1400 16.6400 49.7200 DX = -0.0453 -0.0077 -0.0265 V= 1.0500 0.9735 1.0400 d= 0 -0.0453 -0.0077 iter = 2
DC = -0.0992 0.0217 -0.0509 J= 51.7247 -31.7656 21.3026 -32.9816 65.6564 -15.3791
-28.5386 17.4028 48.1036 DX = -0.0018 -0.0010 -0.0018 V= 1.0500 0.9717 1.0400 d= 0 -0.0471 -0.0087 iter = 3
DC = 1.0e-003 * -0.2166 0.0382 -0.1430 J= 51.5967 -31.6939 21.1474
-32.9339 65.5976 -15.3516 -28.5482 17.3969 47.9549 DX = 1.0e-005 * -0.3856 -0.2386 -0.4412 V= 1.0500 0.9717 1.0400 d= 0 -0.0471 -0.0087
P1 = 2.1842 Q1 = 1.4085 Q3 = 1.4618
(iii) Fast Decoupled Method Exercise : Figure shows the one line diagram of a simple three bus system with generators at buses 1 and 3. The magnitude of voltage at bus 1 is adjusted to 1.05 pu. Voltage magnitude at bus 3 is fixed at 1.04 pu with a real power generation of 200 MW. A load consisting of 400MW and 250 Mvar is taken from bus 2. Line impedance are marked in per unit on a base of 100 MVA and the line charging susceptance are neglected. Obtain the power flow solution by the NW methods.
Code: Code: % Example 6.12 Fast decoupled method V1= 1.05; V2 = 1.0; V3 = 1.04; d1 = 0; d2 = 0; d3=0; Ps2=-4; Ps3 =2.0; Qs2= -2.5;
YB = [ 20-j*50 -10+j*20 -10+j*20 26-j*52 -10+j*30 -16+j*32
-10+j*30
-16+j*32 26-j*62];
Y = abs(YB); t=angle(YB); B =[-52 32; 32 -62] Binv = inv(B) iter=0; pwracur = 0.0003; % Power accuracy DC = 10;
% Set the max of power mismatch to a high value
while max(abs(DC)) > pwracur iter = iter +1; P2= V2*V1*Y(2,1)*cos(t(2,1)-d2+d1)+V2^2*Y(2,2)*cos(t(2,2))+ ... V2*V3*Y(2,3)*cos(t(2,3)-d2+d3); P3= V3*V1*Y(3,1)*cos(t(3,1)-d3+d1)+V3^2*Y(3,3)*cos(t(3,3))+ ... V3*V2*Y(3,2)*cos(t(3,2)-d3+d2); Q2=-V2*V1*Y(2,1)*sin(t(2,1)-d2+d1)-V2^2*Y(2,2)*sin(t(2,2))- ... V2*V3*Y(2,3)*sin(t(2,3)-d2+d3); DP2 = Ps2 - P2; DP2V = DP2/V2; DP3 = Ps3 - P3; DP3V = DP3/V3; DQ2 = Qs2 - Q2; DQ2V = DQ2/V2; DC =[DP2; DP3; DQ2]; Dd = -Binv*[DP2V;DP3V]; DV = -1/B(1,1)*DQ2V;
d2 =d2+Dd(1); d3 =d3+Dd(2); V2= V2+DV; angle2 =180/pi*d2; angle3 =180/pi*d3; R = [iter d2 d3 V2 DP2 DP3 DQ2]; disp(R) end Q3=-V3*V1*Y(3,1)*sin(t(3,1)-d3+d1)-V3^2*Y(3,3)*sin(t(3,3))- ... V3*V2*Y(3,2)*sin(t(3,2)-d3+d2); P1= V1^2*Y(1,1)*cos(t(1,1))+V1*V2*Y(1,2)*cos(t(1,2)-d1+d2)+ ... V1*V3*Y(1,3)*cos(t(1,3)-d1+d3); Q1=-V1^2*Y(1,1)*sin(t(1,1))-V1*V2*Y(1,2)*sin(t(1,2)-d1+d2)- ... V1*V3*Y(1,3)*sin(t(1,3)-d1+d3); S1=P1+j*Q1 Q3 B= -52
32
32 -62 Binv = -0.0282 -0.0145 -0.0145 -0.0236
Columns 1 through 6 1.0000 -0.0605 -0.0089
0.9958 -2.8600
1.4384
Column 7 -0.2200 Columns 1 through 6 2.0000 -0.0565 -0.0080
0.9653
0.1759 -0.0710
0.9657
0.6403 -0.4570
Column 7 -1.5790 Columns 1 through 6 3.0000 -0.0442 -0.0087 Column 7 0.0219 Columns 1 through 6 4.0000 -0.0448 -0.0090
0.9730 -0.0214
0.0012
0.9731 -0.1534
0.1129
Column 7 0.3652 Columns 1 through 6 5.0000 -0.0477 -0.0087 Column 7 0.0067 Columns 1 through 6
6.0000 -0.0476 -0.0086
0.9714
0.0005
0.0026
0.9713
0.0360 -0.0262
0.9717
0.0009 -0.0014
Column 7 -0.0861 Columns 1 through 6 7.0000 -0.0469 -0.0087 Column 7 -0.0041 Columns 1 through 6 8.0000 -0.0469 -0.0087 Column 7 0.0201 Columns 1 through 6 9.0000 -0.0471 -0.0087
0.9718 -0.0084
0.0061
0.9717 -0.0005
0.0005
Column 7 0.0016 Columns 1 through 6 10.0000 -0.0471 -0.0087 Column 7 -0.0047 Columns 1 through 6 11.0000 -0.0471 -0.0087 Column 7 -0.0005
0.9717
0.0020 -0.0014
Columns 1 through 6 12.0000 -0.0471 -0.0087
0.9717
0.0002 -0.0002
Column 7 0.0011 Columns 1 through 6 13.0000 -0.0471 -0.0087
0.9717 -0.0005
0.0003
0.9717 -0.0001
0.0000
Column 7 0.0002 Columns 1 through 6 14.0000 -0.0471 -0.0087 Column 7 -0.0003 S1 = 2.1843 + 1.4085i Q3 = 1.4617
Experiment No. 3
Object :Study of economic load dispatch problem with different methods. Software Tool:MATLAB Theory: Economic Distribution of Loads between the Units of a Plant To determine the economic distribution of a load amongst the different units of a plant, the variable operating costs of each unit must be expressed in terms of its power output. The fuel cost is the main cost in a thermal or nuclear unit. Then the fuel cost must be expressed in terms of the power output. Other costs, such as the operation and maintenance costs, can also be expressed in terms of the power output. Fixed costs, such as the capital cost, depreciation etc., are not included in the fuel cost. The fuel requirement of each generator is given in terms of the Rupees/hour. Let us define the input cost of an unit- i ,fi in Rs./h and the power output of the unit as Pi. Then the input cost can
Rs./h be expressed in terms of the power output as
(1)
The operating cost given by the above quadratic equation is obtained by approximating the power in MW versus the cost in Rupees curve. The incremental operating cost of each unit is
(2) Rs./MWh then computed as Let us now assume that only two units having different incremental costs supply a load. There will be a reduction in cost if some amount of load is transferred from the unit with higher incremental cost to the unit with lower incremental cost. In this fashion, the load is transferred from the less efficient unit to the more efficient unit thereby reducing the total operation cost. The load transfer will continue till the incremental costs of both the units are same. This will be optimum point of operation for both the units. The above principle can be extended to plants with a total of N number of units. The total fuel cost will then be the summation of the individual fuel cost fi,i = 1, ... , N of each unit, i.e.,
(3)
Let us denote that the total power that the plant is required to supply by PT, such that where P1, ... , PNare the power supplied by the N different units. The objective is minimize fTfor a given PT. This can be achieved when the total difference
(5) dfTbecomes zero, i.e.,
(6) Now since the power supplied is assumed to be constant we have
(7) Multiplying (6) by λ and subtracting from (5) we get
(8) The equality in (7) is satisfied when each individual term given in brackets is zero. This gives us Also the partial derivative becomes a full derivative since only the term fiof fTvaries with Pi, i = 1, ...,N . We then have
(9)
Generating Limits It is not always necessary that all the units of a plant are available to share a load. Some of the units may be taken off due to scheduled maintenance. Also it is not necessary that the less efficient units are switched off during off peak hours. There is a certain amount of shut down and start up costs associated with shutting down a unit during the off peak hours and servicing it back on-line during the peak hours. To complicate the problem further, it may take about eight hours or more to restore the boiler of a unit and synchronizing the unit with the bus. To meet the sudden change in the power demand, it may therefore be necessary to keep more units than it necessary to meet the load demand during that time. This safety margin in generation is called spinning reserve . The optimal load dispatch problem must then incorporate this startup and shut down cost for without endangering the system security. (10) The power generation limit of each unit is then given by the inequality constraints The maximum limit Pmaxis the upper limit of power generation capacity of each unit. On the other hand, the lower limit Pminpertains to the thermal consideration of operating a boiler in a thermal or nuclear generating station. An operational unit must produce a minimum amount of power such that the boiler thermal components are stabilized at the minimum design operating temperature Economic Sharing of Loads between Different Plants So far we have considered the economic operation of a single plant in which we have discussed how a particular amount of load is shared between the different units of a plant. In this problem we did not have to consider the transmission line losses and assumed that the losses were a part of the load supplied. However if now consider how a load is distributed between the different plants that are joined by transmission lines, then the line losses have to be explicitly included in the economic dispatch problem. In this section we shall discuss this problem.
(11) When the transmission losses are included in the economic dispatch problem, we can modify wherePLOSSis the total line loss. Since PTis assumed to be constant, we have (12 ) In the above equation dPLOSSincludes the power loss due to every generator, i.e.
(13)
Also minimum generation cost implies dfT = 0 as given in (5). Multiplying both (12) and (13) by
(14) λ and combining we get
Adding (14) with (5) we obtain
(16)
The above equation satisfies when Again since
from (16) we get (18) whereLiis called the penalty factor of load- i and is given by
Consider an area with N number of units. The power generated are defined by the vector (19) Then the transmission losses are expressed in general as whereB is a symmetric matrix given by
The elements Bijof the matrix B are called the loss coefficients. These coefficients are not constant but vary with plant loading. However for the simplified calculation of the penalty factor Lithese coefficients are often assumed to be constant. When the incremental cost equations are linear, we can use analytical equations to find out the economic settings. However in practice, the incremental costs are given by nonlinear equations that may even contain nonlinearities. In that case iterative solutions are required to find the optimal generator settings. Code: %Economical Load dispatch assuming 0 transmission losses %input quantity is no of unit %this code is limitation is that no of unit less then one n=input('enter the no of unit'); kr=input('input total power') %The form of IFC is %ifc1=dC1/dP1=x1*p1+C1; %ifc2=dC2/dP2=x2*p2+C2; for t=1:n disp('for unit=') t x(t)=input('enter the value of xi'); c(t)=input('enter the value of ci'); end if n==2
a=[1 1 ;x(1) -x(2)] p=inv(a)*[kr;c(1)-c(2)] else if n==3 a=[1 1 1; x(1) -x(2) 0; 0 x(2) -x(3)] p=inv(a)*[kr;c(1)-c(2);c(3)-c(2)] else n==4 a=[1 1 1 1; x(1) -x(2) 0 0; 0 x(2) -x(3) 0; 0 0 x(3) -x(4)] p=inv(a)*[kr;c(1)-c(2);c(3)-c(2);c(4)-c(3)] end end Output : enter the no of unit3 input total power1000 kr = 1000 for unit= t= 1 enter the value of xi.3 enter the value of ci.4 for unit= t=2
enter the value of xi.2 enter the value of ci.3 for unit= t= 3 enter the value of xi.4 enter the value of ci.3 a= 1.0000
1.0000
1.0000
0.3000 -0.2000 0
0.2000 -0.4000
p= 307.9231 461.3846 230.6923
0
EXPERIMENT NO. 4 Object :Study of transient stability analysis using MATLAB/ETAP Software. Software Tool :Matlab Theory: Equal Area Criterion The real power transmitted over a lossless line is given by
(1) Now consider the situation in which the synchronous machine is operating in steady state delivering a power Psequal to Pmwhen there is a fault occurs in the system. Opening up of the circuit breakers in the faulted section subsequently clears the fault. The circuit breakers take about 5/6 cycles to open and the subsequent post-fault transient last for another few cycles. The input power, on the other hand, is supplied by a prime mover that is usually driven by a steam turbine. The time constant of the turbine mass system is of the order of few seconds, while the electrical system time constant is in milliseconds. Therefore, for all practical purpose, the mechanical power is remains constant during this period when the electrical transients occur. The transient stability study therefore concentrates on the ability of the power system to recover from the fault and deliver the constant power Pmwith a possible new load angle δ. Consider the power angle curve shown in Fig. Suppose the system is operating in the steady state delivering a power of Pmat an angle of δ0 when due to malfunction of the line, circuit breakers open reducing the real power transferred to zero. Since Pmremains constant, the accelerating power Pabecomes equal to Pm. The difference in the power gives rise to the rate of change of stored kinetic energy in the rotor masses. Thus the rotor will accelerate under the constant influence of non-zero accelerating power and hence the load angle will increase. Now suppose the circuit breaker re-closes at an angleδc.The power will then revert back to the normal operating curve. At that point, the electrical power will be more than the mechanical power and the accelerating power will be negative. This will cause the machine decelerate. However, due to the inertia of the rotor masses, the load angle will still keep on increasing. The increase in this angle may eventually stop and the rotor may start decelerating, otherwise the system will lose synchronism. Note that
(2)
Fig.1 Power-angle curve for equal area criterion.
(3) Hence multiplying both sides of above equation by
and rearranging
(4) Multiplying both sides of the above equation by dtand then integrating between two arbitrary angles δ0 Now suppose the generator is at rest at δ0. We then have dδ/ dt= 0. Once a fault occurs, the machine starts accelerating. Once the fault is cleared, the machine keeps on accelerating before it reaches its peak at δc, at which point we again have dδ / dt= 0. Thus the area of accelerating is given from as
(5) In a similar way, we can define the area of deceleration. In Fig., the area of acceleration is given by A1 while the area of deceleration is given by A2.
(6)
Exercise:1 A loss-free generator supplies 0.5 p.u. to an infinite bus , the steady state limit of the system being 1 p.u. Determine whether the generator will remain in synchronism if the prime mover input is abruptly increased by 0.3 p.u.Write a MATLAB code to find the transient stability using equal area criterion.
Code: p_m=input('enter the mechanical output in per unit') x=linspace(0,pi,1000); x1=asin(p_m/4) y=p_e*sin(x); plot(x,y) hold on y=x*.000000001+p_m; x_crit=pi-x1 plot(x,y) hold off %for checking stability % using equal area criterion y=quad(@myfun,0,x1); acc_area=p_m*(x1-0)-y y1=quad(@myfun,x1,x_crit); de_acc_area=y1-p_m*(x_crit-x1)
if (acc_area
Output: enter the mechanical output in per unit =3 p_m = 3 x1 = 0.8481 x_crit = 2.2935 y= 0 0.4596 0.9131 1.6458 2.3168 2.6761 3.0000 y= 0.2302 0.6875 y= 1.2855 1.9906 y= 2.5006 2.8428 acc_area = 1.1899 y= 3.0000 3.4584 3.7839 4.0000 3.7839 3.4584 3.0000 y= 3.2448 3.6387 y= 3.9456 3.9456
y= 3.8780 3.9864 y= 3.9864 3.8780 y= 3.6387 3.2448 de_acc_area = 0.9551 system is un_stable
