应用介绍

免费下载用于辐射状配电网系统的潮流计算的matlab代码，配电网潮流计算matlab代码。

clc;

clear all;

format short;

tic

m=load('loaddata12bus.m');

l=load('linedata12bus.m'); 

br=length(l);

no=length(m);

MVAb=100;

KVb=11;

Zb=(KVb^2)/MVAb;

% Per unit Values

for i=1:br

    R(i,1)=(l(i,4))/Zb;

    X(i,1)=(l(i,5))/Zb;

end

for i=1:no

    P(i,1)=((m(i,2))/(1000*MVAb));

    Q(i,1)=((m(i,3))/(1000*MVAb));

end

R

X

P

Q

C=zeros(br,no);

for i=1:br

    a=l(i,2);

    b=l(i,3);

    for j=1:no

        if a==j

            C(i,j)=-1;

        end

        if b==j

            C(i,j)=1;

        end

    end

end

C

e=1;

for i=1:no

    d=0;

    for j=1:br

        if C(j,i)==-1

            d=1;

        end

    end

    if d==0

        endnode(e,1)=i;

        e=e+1;

    end

end

endnode

h=length(endnode);

for j=1:h

    e=2;

    f=endnode(j,1);

   % while (f~=1)

   for s=1:no

     if (f~=1)

       k=1;  

       for i=1:br

           if ((C(i,f)==1)&&(k==1))

                f=i;

                k=2;

           end

       end

       k=1;

       for i=1:no

           if ((C(f,i)==-1)&&(k==1));

                f=i;

                g(j,e)=i;

                e=e+1;

                k=3;

           end            

       end

     end

   end

end

for i=1:h

    g(i,1)=endnode(i,1);

end

g;

w=length(g(1,:))

for i=1:h

    j=1;

    for k=1:no 

        for t=1:w

            if g(i,t)==k

                g(i,t)=g(i,j);

                g(i,j)=k;

                j=j+1;

             end

         end

    end

end

g;

for k=1:br

    e=1;

    for i=1:h

        for j=1:w-1

            if (g(i,j)==k) 

                if g(i,j+1)~=0

                    adjb(k,e)=g(i,j+1);            

                    e=e+1;

                else

                    adjb(k,1)=0;

                end

             end

        end

    end

end

adjb;

for i=1:br-1

    for j=h:-1:1

        for k=j:-1:2

            if adjb(i,j)==adjb(i,k-1)

                adjb(i,j)=0;

            end

        end

    end

end

adjb;

x=length(adjb(:,1));

ab=length(adjb(1,:));

for i=1:x

    for j=1:ab

        if adjb(i,j)==0 && j~=ab

            if adjb(i,j+1)~=0

                adjb(i,j)=adjb(i,j+1);

                adjb(i,j+1)=0;

            end

        end

        if adjb(i,j)~=0

            adjb(i,j)=adjb(i,j)-1;

        end

    end

end

adjb;

for i=1:x-1

    for j=1:ab

        adjcb(i,j)=adjb(i+1,j);

    end

end

b=length(adjcb);

% voltage current program

for i=1:no

    vb(i,1)=1;

end

for s=1:10

for i=1:no

    nlc(i,1)=conj(complex(P(i,1),Q(i,1)))/(vb(i,1));

end

nlc;

for i=1:br

    Ibr(i,1)=nlc(i+1,1);

end

Ibr;

xy=length(adjcb(1,:));

for i=br-1:-1:1

    for k=1:xy

        if adjcb(i,k)~=0

            u=adjcb(i,k);

            %Ibr(i,1)=nlc(i+1,1)+Ibr(k,1);

            Ibr(i,1)=Ibr(i,1)+Ibr(u,1);

        end

    end      

end

Ibr;

for i=2:no

      g=0;

      for a=1:b 

          if xy>1

            if adjcb(a,2)==i-1 

                u=adjcb(a,1);

                vb(i,1)=((vb(u,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));

                g=1;

            end

            if adjcb(a,3)==i-1 

                u=adjcb(a,1);

                vb(i,1)=((vb(u,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));

                g=1;

            end

          end

        end

        if g==0

            vb(i,1)=((vb(i-1,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));

        end

end

s=s+1;

end

nlc;

Ibr;

vb

vbp=[abs(vb) angle(vb)*180/pi]

toc;

for i=1:no

    va(i,2:3)=vbp(i,1:2);

end

for i=1:no

    va(i,1)=i;

end

va;

Ibrp=[abs(Ibr) angle(Ibr)*180/pi];

PL(1,1)=0;

QL(1,1)=0;

% losses

for f=1:br

    Pl(f,1)=(Ibrp(f,1)^2)*R(f,1);

    Ql(f,1)=X(f,1)*(Ibrp(f,1)^2);

    PL(1,1)=PL(1,1)+Pl(f,1);

    QL(1,1)=QL(1,1)+Ql(f,1);

end

Plosskw=(Pl)*100000

Qlosskw=(Ql)*100000

PL=(PL)*100000

QL=(QL)*100000

voltage = vbp(:,1)

angle = vbp(:,2)*(pi/180)