newton_1.asv

可进行电力系统多节点系统的优化潮流计算

%%%%%%%%%%%%%%%%%power flow %%%%%%%%%
%%%%%%%%%newton%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%newton（补偿前）%%%%%%%%%%%%%%%
%Lfybus;
function (energy,Vm,t) = newton(Qg,QL,tap);
%%%%%%%%%%%%%%%%%%
j=sqrt(-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
    RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
    GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nbus = length(busdata(:,1));
Ym=abs(Ybus); 
t = angle(Ybus);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    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); 
%%%%%%%%%%对节点10，24进行补偿%%%%%%%%%%%%
  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
j2=0;NNSlack=0;NNPQ=0;
for k=1:nbus
    if kb(k) ~= 1
        j2=j2+1; NNSlack(j2)=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;
maxerror = 1; 
converge=1;
iter = 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; 
            Q(ii)=QK;  
        end
        if kb(ii) == 2
           if Qmax(ii) ~= 0
               Qgc = Q(ii)*basemva + Qd(ii) - Qsh(ii);
               if  iter <= 20                 % 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
deltad;
abs(V)'
busdata(:,3)=Vm';
busdata(:,4)=deltad';
Pgt = sum(Pg);
Qgt = sum(Qg);
Pdt = sum(Pd);
Qdt = sum(Qd);
toc