clear all;
close all;
clc;
%network data matrix: [buses|lines|transformers|shunt elements|base MVA|PV|PQ]
ND=[6 5 2 2 100 1 4];
%transmission lines data
TL=[1 1 6 0.123 0.518 0 55;
2 1 4 0.080 0.370 0 65;
3 4 6 0.087 0.407 0 30;
4 5 2 0.282 0.640 0 55;
5 2 3 0.723 1.050 0 40];
%Transformer Data
TR=[1 6 5 0 0.300 1 30;
2 4 3 0 0.133 1 55];
%Bus Data
Slack=[1 00 200 40 0 0 100 -50 1.02];
PV=[2 50 100 20 0 0 50 -25 1.02];
PQ=[3 55 13 1;
4 0 0 1;
5 30 18 1;
6 50 5 1];
%Shunt Element Data
SE=[1 4 2;
2 6 2.5];
n=ND(:,1);
l=ND(:,2);
t=ND(:,3);
s=ND(:,4);
bmva=ND(1,5);
pv=ND(:,6);
pq=ND(:,7);
HLCS=complex(0,TL(:,6));
TLZ=complex(TL(:,4),TL(:,5));
TRZ=complex(TR(:,4),TR(:,5));
y=zeros(n);
z=zeros(n);
for i=1:n
for j=1:n
if (i==j)
for k=1:l
if(TL(k,2)==i||TL(k,3)==i)
y(i,j)=y(i,j)+1/TLZ(k,1)+HLCS(k,1);
end
end
for k=1:t
if(TR(k,2)==i||TR(k,3)==i)
y(i,j)=y(i,j)+1/TRZ(k,1);
end
end
else
for k=1:l
if((TL(k,2)==i&&TL(k,3)==j)||(TL(k,2)==j&&TL(k,3)==i))
y(i,j)=-1/TLZ(k,1);
end
end
for k=1:t
if((TR(k,2)==i&&TR(k,3)==j)||(TR(k,2)==j&&TR(k,3)==i))
y(i,j)=-1/TRZ(k,1);
end
end
end
end
end
display('Admittance Matrix')
disp(y)
acc=input('accleration factor(1 to 1.6)=');
conv=input('convergence factor= ');
k=1;
%type allotment
for i=1:n
type(1)=1;
for j=1:pv
if PV(j,1)==i
type(i)=2;
end
end
for j=1:pq
if PQ(j,1)==i
type(i)=3;
end
end
end
for i=1:n %type check
type(i);
end
for i=1:n
switch type(i)
case 1
V(i,k)=Slack(i,9);
del(i,k)=0;
P(i,k)=0;
Q(i,k)=0;
continue;
case 2
for j=1:pv
if PV(j,1)==i
P(i,k)=PV(j,2)/bmva;
V(i,k)=PV(j,9);
qmax(i)=PV(j,7)/bmva;
qmin(i)=PV(j,8)/bmva;
del(i,k)=0;
Q(i,k)=0;
end
end
continue;
case 3
for j=1:pq
if PQ(j,1)==i
Qgen(i)=0;
for a=1:s
if SE(a,2)==i
Qgen(i)=SE(a,3);
end
end
P(i,k)=-PQ(j,2)/bmva;
Q(i,k)=(Qgen(i)-PQ(j,3))/bmva;
V(i,k)=1;
del(i,k)=0;
end
end
continue;
end
end
for k=1:35
for i=1:n
switch type(i)
case 1
V(i,k+1)=V(i,k);
del(i,k+1)=del(i,k);
P(i,k+1)=P(i,k);
Q(i,k+1)=Q(i,k);
continue;
case 2
AA=0;
for j=1:i-1
AA=AA+y(i,j)*V(j,k+1);
end
for j=i:n
AA=AA+y(i,j)*V(j,k);
end
Q(i,k+1)=-imag(conj(V(i,k))*AA);
if (Q(i,k+1)>qmin(i) && Q(i,k+1)<qmax(i))
AA=0;
for j=1:i-1
AA=AA+y(i,j)*V(j,k+1);
end
for j=i+1:n
AA=AA+y(i,j)*V(j,k);
end
Vtemp=((P(i,k)/conj(V(i,k))-1j*Q(i,k)/conj(V(i,k)))-AA)/y(i,i);
del(i,k+1)=angle(Vtemp);
Vmag=abs(V(i,k));
V(i,k+1)=complex(Vmag*cos(del(i,k+1)),Vmag*sin(del(i,k+1)));
else
if (Q(i,k+1)<qmin(i))
Q(i,k+1)=qmin(i);
elseif (Q(i,k+1)>qmax(i))
Q(i,k+1)=qmax(i);
end
type(i)=3;
V(i,k+1)=1;
del(i,k+1)=0;
end
P(i,k+1)=P(i,k);
continue;
case 3
AA=0;
for j=1:i-1
AA=AA+y(i,j)*V(j,k+1);
end
for j=i+1:n
AA=AA+y(i,j)*V(j,k);
end
V(i,k+1)=((P(i,k)/conj(V(i,k))-1j*Q(i,k)/conj(V(i,k)))-AA)/y(i,i);
del(i,k+1)=angle(V(i,k+1));
V(i,k+1)=V(i,k)+acc*(V(i,k+1)-V(i,k));
delV(i)=abs(abs(V(i,k+1))-abs(V(i,k)));
P(i,k+1)=P(i,k);
Q(i,k+1)=Q(i,k);
continue;
end
end
count=0;
for i=1:n % convergence check
if type(i)==3 && delV(i)<=conv
count=count+1;
end
end
tni=k;
if count>=pq %exit iterations based on conv check
break;
end
end
for i=1:n
if type(i)==1
AA=0;
for j=1:n
AA=AA+y(i,j)*V(j,k+1);
end
end
P(i,k+1)=real(conj(V(i,k+1))*AA);
Q(i,k+1)=-imag(conj(V(i,k+1))*AA);
end
display('Gauss Seidel Load Flow Solution');
display('Total no. of iterations ='); disp(tni);
for i=1:n
disp(i);
OP(i,:)= [abs(V(i,k+1)) del(i,k+1)*180*7/22 P(i,k+1)*bmva Q(i,k+1)*bmva];
disp(OP(i,:))
end
