GSRD.rar_flow

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