1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
//Chapter 12
//Example 12.11
//page 488
//To plot swing curves for fault cleared at 0.275s and 0.08s of a multimachine system
clear;clc;
xd1=%i*0.067;xd2=%i*0.1;
//primitive admittances of the lines
y45=1/(0.018+%i*0.11); B45=%i*0.113;
y51=1/(0.004+%i*0.0235); B51=%i*0.098;
y41=1/(0.007+%i*0.04); B41=%i*0.041;
z24=(%i*0.022);
z35=(%i*0.04);
//Bus data and prefault load-flow values in PU
V1=1.0; P1=-3.8083; Q1=-0.2799; Pl1=0; Ql2=0;
V2=1.0194+%i*0.1475; P2=3.25; Q2=0.6986; Pl2=0; Ql2=0;
V3=1.0121+%i*0.1271; P3=2.10; Q3=0.3110; Pl3=0; Ql3=0;
V4=1.0146+%i*0.0767; P4=0; Q4=1.0; Pl4=1.0; Ql4=0.44;
V5=1.0102+%i*0.0439; P5=0; Q5=0; Pl5=0.5; Ql5=0.16;
// To find voltage behind transient reactances before the occurance of fault
I2=(P2-%i*Q2)/(V2');
E2=V2+I2*xd1;
E1=V1;
I3=(P3-%i*Q3)/(V3');
E3=V3+I3*xd2;
//converting loads into their admittances
Yl4=(Pl4-%i*Ql4)/(V4*V4');
Yl5=(Pl5-%i*Ql5)/(V5*V5');
//forming augmented Bus admittance matrix before the occurance of fault
Y11=y41+y51;Y12=0;Y13=0;Y14=-y41;Y15=-y51;
Y21=Y12;Y22=1/(xd1+z24);Y23=0;Y24=-(1/(z24+xd1));Y25=0;
Y31=0;Y32=0;Y33=1/(z35+xd2);Y34=0;Y35=-1/(z35+xd2);
Y41=Y14;Y42=Y24;Y43=Y34;Y44=y41+Yl4+y45+B45+B41-Y24;Y45=-y45;
Y51=Y15;Y52=Y25;Y53=Y35;Y54=Y45;Y55=Yl5+y45+y51+B45+B51-Y35;
Ybus=[Y11 Y12 Y13 Y14 Y15;
Y21 Y22 Y23 Y24 Y25;
Y31 Y32 Y33 Y34 Y35;
Y41 Y42 Y43 Y44 Y45;
Y51 Y52 Y53 Y54 Y55];
printf('\n Augmented prefault bus admittance matrix (in PU) is given by\n\n Ybus=\n');
disp(Ybus);
///////to find the Ybus during fault
Ybus_1=Ybus([1:3,5],[1:3,5]);
n=4
for k=1:n-1
for j=1:n-1
Ybus_during_fault(k,j)=Ybus_1(k,j)-(Ybus_1(k,n)*Ybus_1(n,j))/Ybus_1(n,n);
end
end
printf('\n\n\n Bus admittance matrix during fault (in PU) is given by\n\n Ybus_during_fault=\n');
disp(Ybus_during_fault);
//to find Ybus after the fault has been cleared
Y45=0;Y54=0;Y44=Y44-y45-B45;Y55=Y55-y45-B45;
Ybus_2=[Y11 Y12 Y13 Y14 Y15;
Y21 Y22 Y23 Y24 Y25;
Y31 Y32 Y33 Y34 Y35;
Y41 Y42 Y43 Y44 Y45;
Y51 Y52 Y53 Y54 Y55];
//eliminating node 5 from Ybus_2
n=5
for k=1:n-1
for j=1:n-1
Ybus_3(k,j)=Ybus_2(k,j)-(Ybus_2(k,n)*Ybus_2(n,j))/Ybus_2(n,n);
end
end
//eliminating node 4 to get post fault Ybus
n=4
for k=1:n-1
for j=1:n-1
Ybus_post_fault(k,j)=Ybus_3(k,j)-(Ybus_3(k,n)*Ybus_3(n,j))/Ybus_3(n,n);
end
end
printf('\n\n\n Bus admittance matrix postfault (in PU) is given by\n\n Ybus_post_fault=\n');
disp(Ybus_post_fault);
printf('\n\n\n');
//During fault power angle equation
delta3=0:0.1:180;
Pe2f=0;
Pe3f=(abs(E3'))^2*real(Ybus_during_fault(3,3))+abs(E1')*abs(E3')*abs(Ybus_during_fault(3,1))*cosd(delta3-atand(imag(Ybus_during_fault(1,3))/real(Ybus_during_fault(1,3))));
//Postfault power angle equations
delta2=0:0.1:180;
Pe2pf=(abs(E2'))^2*real(Ybus_post_fault(2,2))+abs(E1')*abs(E2')*abs(Ybus_post_fault(2,1))*cosd(delta2-atand(imag(Ybus_post_fault(1,2))/real(Ybus_post_fault(1,2))));
Pe3pf=(abs(E3'))^2*real(Ybus_post_fault(3,3))+abs(E1')*abs(E3')*abs(Ybus_post_fault(3,1))*cosd(delta3-atand(imag(Ybus_post_fault(1,3))/real(Ybus_post_fault(1,3))));
//mechanical inputs which are assumed to be constant are given by
Pm2=max(real(E2*I2'));
Pm3=max(real(E3*I3'));
//xdot function defining the swing equations of each of the machines
function xdot=mac2(t,x,tc)
xdot(1)=x(2);
if t>tc then
xdot(2)=180*50*(Pm2-(0.6012+8.365*sind(x(1)-1.662)))/12;//swing equation after clearing the fault
else
xdot(2)=180*50*(Pm2-(0))/12; //swing equation before clearing the fault
end
endfunction
function xdot=mac3(t,x,tc)
xdot(1)=x(2);
if t>tc then
xdot(2)=180*50*(Pm3-(0.1823+6.5282*sind(x(1)-0.8466)))/9;//swing equation after clearing the fault
else
xdot(2)=180*50*(Pm3-(0.1561+5.531*sind(x(1)-0.755)))/9; //swing equation before clearing the fault
end
endfunction
//to find the solution of swing equation to draw the swing curves
//to draw the swing curves for machines 2 and 3 for example12.11 for clearing at 0.275 sec
subplot(2,1,1)
x_1_0=[19.354398,0]';t0=0; T=0:0.01:1;T=T';
x_2_0=[18.2459,0]';tc=0.275;
sol1=ode(x_1_0,t0,T,mac2);
sol2=ode(x_2_0,t0,T,mac3);
plot(T(1:20),sol1(1,1:20)',T,sol2(1,:)');
set(gca(),"grid",[1 1]);
legend('Machine 2','Machine 3',[,1]);
title('Swing Curves for machines 2 and 3 of Example 12.11 for a clearing at '+string(tc)+' s');
xstring(0.55,59,'Machine 1 is reference (infinte bus)');
xlabel('Time (in seconds)----->');
ylabel('Torque Angle (delta,deg)----->');
//to draw the swing curves for machines 2 and 3 for example12.11 for clearing at 0.08 sec
subplot(2,1,2)
x_1_0=[19.354398,0]';t0=0; T=0:0.01:1;T=T';
x_2_0=[18.2459,0]';tc=0.08;
sol1=ode(x_1_0,t0,T,mac2);
sol2=ode(x_2_0,t0,T,mac3);
plot(T,sol1(1,:)',T,sol2(1,:)');
set(gca(),"grid",[1 1]);
legend('Machine 2','Machine 3',[,4]);
title('Swing Curves for machines 2 and 3 of Example 12.11 for a clearing at '+string(tc)+' s');
xstring(0.44,43,'Machine 1 is reference (infinte bus)');
xlabel('Time (in seconds)----->');
ylabel('Torque Angle (delta,deg)----->');
f=get("current_figure");
f.figure_position=[0,15];
f.figure_size=[565,1000];
|
