Question

Answer part D

The single line diagram of a power network is shown in the figure. Bus#1 is a slack bus. The scheduled powers for bus#2 and bus#3 are given. The impedances shown in the figure are all in per-unit considering a power base of 100 MVA. 30 400 MW 320 MVAr Slack V-140 j0.0125 jo.0s 00 MW 270 MVAr A. Use the Gauss Seidel technique to determine voltages at bus#2 & bus#3. (Start with an initial guess 140 for both buses). [Only TWO iterations] D. Add a resistance of R-001 PU, R=0.02 PU and R-0.008 PU in series to lines 1-2, 2-3 and 1-3 respectively. Run a load flow in MATLAB. Compare your results with what you got in part C

Answer 1

% data to be entered
% Line From-bus To-bus R X B Trafo. tap ratio   
line_data= [ 1 1 2 0.01 0.033 0 1   
2 1 3 0.02 0.0125 0 1
3 2 3 0.008 0.05 0 1 ];

% bus no. V Pg Qg Pd Qd Bus Type Angle V_limit
bus_data= [ 1 1 0 0 0 0 1 0 1
2 1 0 0 4 3.2 3 0 1   
3 1 0 0 3 2.7 3 0 1   
];
% Base MVA is 100.
%Bus Type: 1- Slack bus, 2- PV bus, 3- PQ bus
%Pd,Qd: Active and Reactive power consumed by PQ bus (in p.u)
%Pg,Qg: Active and Reactive power generated by PV bus (in p.u)
l=max(line_data(:,1));
n=max(bus_data(:,1));
s=0;
for i=1:n
if((bus_data(i,7)==2)||(bus_data(i,7)==1))
s=s+1;
end
end

% ybus formation
ybus=zeros(n);

for i=1:1:l
  
from_bus=line_data(i,2);
to_bus=line_data(i,3);
b=complex(0,line_data(i,6));
r=line_data(i,4);
x=line_data(i,5);
z=complex(r,x);
y=1./z;
a=line_data(i,7);
if a~=1
  
ybus(from_bus,from_bus)=ybus(from_bus,from_bus)+ y*(1-a)+ y*a;
ybus(to_bus,to_bus)=ybus(to_bus,to_bus)+ y*(a-1)*a+y*a;
ybus(from_bus,to_bus)=ybus(from_bus,to_bus)- y*a;
ybus(to_bus,from_bus)=ybus(to_bus,from_bus)- y*a;
else
  
ybus(from_bus,from_bus)=ybus(from_bus,from_bus)+ y+ b/2;
ybus(to_bus,to_bus)=ybus(to_bus,to_bus)+ y+ b/2;
ybus(from_bus,to_bus)=ybus(from_bus,to_bus)- y;
ybus(to_bus,from_bus)=ybus(to_bus,from_bus)- y;
end
end

% gauss siedal code
V_complex=zeros(n,1);
P=zeros(n,1);
Q=zeros(n,1);
del_V=1;
itr=1;
Z=zeros(l,1);
I=zeros(l,1);
for a=1:n
V_complex(a)= complex(bus_data(a,2),0);
end
for p=2:n
P(p)=(bus_data(p,3)-bus_data(p,5));
Q(p)=(bus_data(p,4)-bus_data(p,6));
P(p)=P(p)/100;
Q(p)=Q(p)/100;
end
while itr<100 && del_V>0.00000001
itr=itr+1;
for i=2:n
t=0;
for k=1:n
if(i~=k)
t=t+(ybus(i,k)*V_complex(k));
end
end
temp1=V_complex(i);
V_complex(i)=((1/ybus(i,i))*(((conj(complex(P(i),Q(i))))/conj(V_complex(i)))-t));
temp2=V_complex(i);
del_V=temp2-temp1;
end
end

%Calculating current values
for i=1:l
f=line_data(i,2);
g=line_data(i,3);
I(i)=-ybus(f,g)*(V_complex(f)-V_complex(g));
end