nrlf303p.m

This is the file load flow in power system

% THIS PROGRAM IS nrlf in polar coordinates
clear all;clc;close;
inputfile;
no_s=find(bus_spec(:,2)==0);
if(length(no_s)==1)
    line=[linedata(:,1) linedata(:,2)];n=length(bus_spec(:,1));
    [linedata(:,1),jj]=find((repmat(linedata(:,1),1,length(bus_spec(:,1)))==repmat(bus_spec(:,1)',length(linedata(:,1)),1))');
    [linedata(:,2),jj]=find((repmat(linedata(:,2),1,length(bus_spec(:,1)))==repmat(bus_spec(:,1)',length(linedata(:,2)),1))');
    pv=find(bus_spec(:,2)==1);pq=find(bus_spec(:,2)==2);
    lambda=input('Enter the load multiplier lambda value :');
    iter=1;eps=1e-4;%input('Enter the tolerance value:');
    choice=1;
    tic;
    Ybus=sparse(linedata(:,1),linedata(:,2),-1./(linedata(:,3).*conj(linedata(:,5))),n,n)+sparse(linedata(:,2),linedata(:,1),-1./(linedata(:,3).*linedata(:,5)),n,n)...
        +sparse(linedata(:,1),linedata(:,1),linedata(:,4)+1./(linedata(:,3).*abs(linedata(:,5)).^2),n,n)+sparse(linedata(:,2),linedata(:,2),linedata(:,4)+1./linedata(:,3),n,n)+sparse(1:n,1:n,bus_spec(:,10),n,n);
    mismatch=zeros(2*(n-1),1);Np=[];Nq=[];id1=zeros(length(pv),1);id2=id1;
    G=real(Ybus);B=imag(Ybus);Vn=bus_spec(:,4);
    Vmag=bus_spec(:,3);Ang=bus_spec(:,5)*pi/180;
    [E,F]=pol2cart(Ang,Vmag);V=E+j*F;
    if choice~=1
        Np=input('input Np :');Nq=input('input Nq :');end

    while (max(abs(mismatch))>=eps&iter<=300|iter==1)
        Pcalc=real(V.*conj(Ybus*V));Qcalc=imag(V.*conj(Ybus*V));
        if choice==1
            ap=(bus_spec(:,11)+bus_spec(:,12).*abs(V./Vn)+bus_spec(:,13).*abs(V./Vn).^2);
            bp=(bus_spec(:,12).*1./(abs(Vn))+2*bus_spec(:,13).*abs(V)./abs(Vn).^2);
            bq=bp;aq=ap;
        else
            ap=abs(V./Vn).^Np;aq=abs(V./Vn).^Nq;
            bp=Np*abs(V./Vn).^(Np-1).*abs(V);bq=Nq*abs(V./Vn).^(Nq-1).*abs(V);
        end
        Psp=real(bus_spec(:,6))*(1+lambda)-real(bus_spec(:,7))*(1+lambda).*ap;
        Qsp=imag(bus_spec(:,6))*(1+lambda)-imag(bus_spec(:,7))*(1+lambda).*aq;
        if iter>1
            Vcheck=repmat(V,1,length(pv))+sparse(pv,[1:length(pv)],bus_spec(pv,3).*(cos(angle(V(pv)))+j*sin(angle(V(pv))))-V(pv,1));
            Icheck=diag(Ybus(pv,:)*Vcheck);
            Qcheck=imag(nonzeros(sparse(pv,[1:length(pv)],bus_spec(pv,3).*(cos(angle(V(pv)))+j*sin(angle(V(pv)))))).*conj(Icheck))+imag(bus_spec(pv,7))*(1+lambda).*aq(pv);
            id1=Qcheck>=bus_spec(pv,8);id2=bus_spec(pv,9)>=Qcheck;
            Qsp(pv(id1))=bus_spec(pv(id1),8) -imag(bus_spec(pv(id1),7))*(1+lambda).*aq(pv(id1));
            Qsp(pv(id2))=bus_spec(pv(id2),9)-imag(bus_spec(pv(id2),7)).*(1+lambda).*aq(pv(id2));
        end
        pvpq=pv(id1|id2);no_pv=pv(~(id1|id2));Vmag(no_pv)=bus_spec(no_pv,3);EE=sparse(repmat(E,1,n));FF=sparse(repmat(F,1,n));
        J1=(G.*EE'-B.*FF').*FF-(G.*FF'+B.*EE').*EE;J2=(G.*EE'-B.*FF').*EE+(G.*FF'+B.*EE').*FF;
        J{iter}=[spdiags(-Qcalc-diag(B).*abs(V).^2,[0],J1) spdiags(Pcalc+diag(G).*abs(V).^2+real(bus_spec(:,7)).*(1+lambda).*bp,[0],J2);
            spdiags(Pcalc-diag(G).*abs(V).^2,[0],-J2) spdiags(Qcalc-diag(B).*abs(V).^2+imag(bus_spec(:,7)).*(1+lambda).*bq,[0],J1)];
        mismatch=[Psp([pq;pv])-Pcalc([pq;pv]);Qsp([pq;pvpq])-Qcalc([pq;pvpq])];
        inc=inv(J{iter}([pq;pv;pq+n;pvpq+n],[pq;pv;pq+n;pvpq+n]))*mismatch;
        Ang([pq;pv])=Ang([pq;pv])+inc([1:n-1]);Vmag([pq;pvpq])=Vmag([pq;pvpq])+inc([n:end]).*abs(V([pq;pvpq]));
        ang{iter}(:,:)=Ang;vmag{iter}(:,:)=Vmag;
        iter=iter+1;
        [E,F]=pol2cart(Ang,Vmag);V=E+j*F;
    end
    toc
    if iter>300
        disp('Divergence Encountered!!!.');
    else
        Sd=bus_spec(:,7);
        S=V.*conj(Ybus*V);SG=S+Sd;
        flow1=V(linedata(:,1)).*conj((V(linedata(:,1))-V(linedata(:,2)))./(linedata(:,3).*linedata(:,5))+V(linedata(:,1)).*(linedata(:,4)+(1-linedata(:,5))./(linedata(:,5).^2.*linedata(:,3))));
        flow2=V(linedata(:,2)).*conj((V(linedata(:,2))-V(linedata(:,1)))./(linedata(:,3).*linedata(:,5))+V(linedata(:,2)).*(linedata(:,4)+(linedata(:,5)-1)./(linedata(:,5).*linedata(:,3))));
        bus_spec=[bus_spec(:,1) bus_spec(:,2) bus_spec(:,3) bus_spec(:,4) bus_spec(:,5) real(bus_spec(:,6)) imag(bus_spec(:,6)) real(bus_spec(:,7))...
            imag(bus_spec(:,7)) bus_spec(:,8) bus_spec(:,9) imag(bus_spec(:,10)) bus_spec(:,11) bus_spec(:,12) bus_spec(:,13)];
        linedata=[line(:,1) line(:,2) real(linedata(:,3)) imag(linedata(:,3)) imag(linedata(:,4)) real(linedata(:,5)) imag(linedata(:,5))];
        printout('NRLF(POLAR) INPUT.txt','NRLF(polar form) program inputfile\n','load multipier (lambda):',lambda,...
            'busdata :','[Busno]      [status]     [|V|]       [Vn]        [theta]     [Pg]        [Qg]        [Pd]        [Qd]        [Qmax]      [Qmin]      [Shunt B]   [Cp]        [Cc]        [Ci] ',...
            bus_spec,'linedata :','[from bus]   [to bus]      [R]         [X]       [(B/2)]     [r_tap]     [i_tap]',linedata);
        printout('NRLF(POLAR) OUTPUT.txt','NRLF(polar form) program results\n','load multipier (lambda):',lambda,'total no. of iterations for convergence','',iter-1,'voltages after convergence :','   Vmag      Ang(deg)',[Vmag Ang*180/pi],...
            'slack bus :','',bus_spec(no_s,1),'pv buses after convergence :','',bus_spec(no_pv,1),'powers :','  busno       P(inj)      Q(inj)      P(gen)      Q(gen)',[bus_spec(:,1)  real(S)  imag(S) real(SG)  imag(SG)],...
            'power flows :','[From bus(i)][To bus(j)]   Pij         Qij         Pji         Qji',[line(:,1) line(:,2) real(flow1)  imag(flow1)  real(flow2)  imag(flow2)],'[total real power loss    total reactive power loss]','',[sum(real(S)) sum(imag(S))],...
            '[total real power generation    total reactive power generation ]','',[sum(real(SG)) sum(imag(SG))],'[total real power demand    total reactive power demand]','',[sum(real(Sd)) sum(imag(Sd))]);
    end
else
    disp('check your slack bus');
end
% open('NRLF(POLAR) OUTPUT.txt');