newton.m

短路电流计算 短路电 流计算 短路电流计算

Vm=0; delta=0;nSlack=0;nPQ=0;maxerror = 1; converge=1;iter = 0;
nbus = length(busdata(:,1));Ym=abs(Ybus); t = angle(Ybus);
for k=1:nbus
    n=busdata(k,1); kb(n)=busdata(k,2); Vm(n)=busdata(k,3); delta(n)=busdata(k, 4);
    Pd(n)=busdata(k,5); Qd(n)=busdata(k,6); Pg(n)=busdata(k,7); Qg(n) = busdata(k,8);
    Qmin(n)=busdata(k, 9); Qmax(n)=busdata(k, 10); Qsh(n)=busdata(k, 11);
    if Vm(n) <= 0  Vm(n) = 1.0; 
    else delta(n) = pi/180*delta(n);
         P(n)=(Pg(n)-Pd(n))/basemva;Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;
         S(n) = P(n) + j*Q(n);
    end
end
jj=0;NNSlack=0;NNPQ=0;
for k=1:nbus
    if kb(k) ~= 1
        jj=jj+1; NNSlack(jj)=k;     
    else
        nSlack = nSlack+1;    %  nSlack:number of slack bus
    end
    if kb(k) == 0, nPQ = nPQ+1;NNPQ(nPQ)=k;end   %  nPQ:number of PQ bus
end
m=nbus-nSlack;n=m+nPQ;Jacobain(n,n)=0;
while maxerror >= accuracy & iter <= maxiter % Test for max. power mismatch
    iter =iter+1;
    Jacobain=0; 
    for ii=1:n
        for jj=1:n
            if ii<=m
                R=NNSlack(ii);
                PK=0;QK=0; J22=0;
                for l=1:nbus
                   if Ym(R,l)~=0
                       cita=t(R,l)- delta(R) + delta(l);
                       PK=PK+Vm(R)*Vm(l)*Ym(R,l)*cos(cita);
                       QK=QK-Vm(R)*Vm(l)*Ym(R,l)*sin(cita);
                       J22=J22+Vm(l)*Ym(R,l)*cos(cita);
                   end
                end
                if kb(R)==2,Q(R)=QK;end     % Swing Q of PV bus
                DC(ii)=P(R)-PK;              % delta P(ii)
                if jj<=m
                   C=NNSlack(jj);
                   if R~=C       %off diagonalelements of J1
                      Jacobain(ii,jj)=-Vm(R)*Vm(C)*Ym(R,C)*sin(t(R,C)- delta(R) + delta(C)); 
                   else
                      Jacobain(ii,jj)=QK-Vm(R)*Vm(R)*Ym(R,R)*sin(t(R,R));%diagonalelements of J1
                   end
                else  
                   C=NNPQ(jj-m);
                   if R~=C   %off diagonalelements of J2
                      Jacobain(ii,jj)=Vm(R)*Ym(R,C)*cos(t(R,C)- delta(R) + delta(C));
                   else
                     Jacobain(ii,jj)=J22+Vm(R)*Ym(R,R)*cos(t(R,R));  %diagonalelements of J2
                   end
                end 
            else
                R=NNPQ(ii-m);
                QK=0;J33=0;J44=0;
                for l=1:nbus
                    if Ym(R,l)~=0
                        cita=t(R,l)- delta(R) + delta(l);
                        QK=QK-Vm(R)*Vm(l)*Ym(R,l)*sin(cita);
                        J33=J33+Vm(R)*Vm(l)*Ym(R,l)*cos(cita);
                        J44=J44-Vm(l)*Ym(R,l)*sin(cita);
                    end
                end
                DC(ii)=Q(R)-QK;  % delta Q(ii)
                if jj<=m
                    C=NNSlack(jj);
                    if R~=C              %off diagonalelements of J3
                        Jacobain(ii,jj)=-Vm(R)*Vm(C)*Ym(R,C)*cos(t(R,C)- delta(R)+delta(C));
                    else
                        Jacobain(ii,jj)=J33-Vm(R)*Vm(R)*Ym(R,R)*cos(t(R,R));%diagonalelements of J3
                    end
                else
                    C=NNPQ(jj-m);
                    if R~=C            %off diagonalelements of J4
                        Jacobain(ii,jj)=-Vm(R)*Ym(R,C)*sin(t(R,C)- delta(R) + delta(C));
                    else
                        Jacobain(ii,jj)=J44-Vm(R)*Ym(R,R)*sin(t(R,R));   %diagonalelements of J4
                    end
                 end
            end
        end
    end
    for ii=1:nbus
        if kb(ii)==1
            PK=0;QK=0;
            for jj=1:nbus
                if Ym(ii,jj)~=0
                    cita=t(R,l)- delta(R) + delta(l);
                    PK=PK+Vm(ii)*Vm(jj)*Ym(ii,jj)*cos(cita);
                    QK=QK-Vm(ii)*Vm(jj)*Ym(ii,jj)*sin(cita);
                end
            end
            P(ii)=PK;  % Swing P of slack bus
            Q(ii)=QK;  % Swing Q of slack bus
        end
        if kb(ii) == 2
           if Qmax(ii) ~= 0
               Qgc = Q(ii)*basemva + Qd(ii) - Qsh(ii);
               if iter <= 7                  % Between the 2th & 6th iterations
                  if iter > 2                % the Mvar of generator buses are
                     if Qgc  < Qmin(ii),       % tested. If not within limits Vm(n)
                         Vm(ii) = Vm(ii) + 0.01;    % is changed in steps of 0.01 pu to
                     elseif Qgc  > Qmax(ii),   % bring the generator Mvar within
                         Vm(ii) = Vm(ii) - 0.01; % the specified limits.
                     end
                  end
               end
            end
        end
    end
    DX=Jacobain\DC';
    for ii=1:n
        if ii<=m
             delta(NNSlack(ii))=delta(NNSlack(ii))+DX(ii);
        else
             Vm(NNPQ(ii-m))=Vm(NNPQ(ii-m))+DX(ii);
        end
    end
    maxerror=max(abs(DC));
     if iter == maxiter & maxerror > accuracy 
   fprintf('\nWARNING: Iterative solution did not converged after ')
   fprintf('%g', iter), fprintf(' iterations.\n\n')
   fprintf('Press Enter to terminate the iterations and print the results \n')
   converge = 0; pause, else, end
end
%----------- end while ----------------------
if converge ~= 1
   tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 
   tech=('                   Power Flow Solution by Newton-Raphson Method');
end   
V = Vm.*cos(delta)+j*Vm.*sin(delta);
deltad=180/pi*delta;
for n = 1:nbus
     if kb(n) == 1
     S(n)= P(n)+j*Q(n);
     Pg(n) = P(n)*basemva + Pd(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);    
     elseif  kb(n) ==2
     S(n)=P(n)+j*Q(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);  
  end
end
busdata(:,3)=Vm'; busdata(:,4)=deltad';
Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);