xhtxt.m

潮流计算

clear
a=textread('IEEE4.txt');
t=cputime;
for n=1:a(1,1)
V(1,n)=1;S(1,n)=0;
end
V(a(1,2))=a(1,3);
n=0;h=2;j=0;k=0;m=0;d=0;branchi=[];
while n<5
    h=h+1;
    if a(h,1)==0
        n=n+1;
    end   
    if n<2
        j=j+1; linei(j)=a(h,2);linej(j)=a(h,3);liner(j)=a(h,4);linex(j)=a(h,5);lineb(j)=a(h,6);
    end
    if ((n==2)&(a(h+1,1)~=0))
        k=k+1;branchi(k)=a(h+1,2);branchg(k)=a(h+1,3);branchb(k)=a(h+1,4);
    end 
    if  ((n==3)&(a(h+1,1)~=0)) 
        m=m+1;powi(m)=a(h+1,2);Pgb(powi(m))=a(h+1,3);Qgb(powi(m))=a(h+1,4);Pdi(powi(m))=a(h+1,5);Qdi(powi(m))=a(h+1,6);
    end
    if ((n==4)&(a(h+1,1)~=0))
        d=d+1;pvi(d)=a(h+1,2);pvv(d)=a(h+1,3);pvql(d)=a(h+1,4);pvqt=a(h+1,5);V(a(h+1,2))=a(h+1,3);
    end
end   
 n=0;h=3;
while a(h,1)~=0
    n=n+1;h=h+1;
end
G=zeros(a(1,1));B=G;
for k=1:n;
h=linei(k);j=linej(k);r=liner(k);x=linex(k);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);G(h,j)=-GIJ;G(j,h)=-GIJ;B(j,h)=-BIJ;B(h,j)=-BIJ;G(h,h)=G(h,h)+GIJ;
b=lineb(k);B(h,h)=B(h,h)+BIJ+b;G(j,j)=G(j,j)+GIJ;B(j,j)=B(j,j)+BIJ+b;
end
for k=(n+2):numel(linei)
h=linei(k);j=linej(k);r=liner(k);x=linex(k);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);k0=lineb(k);G(h,j)=-GIJ/k0;G(j,h)=-GIJ/k0;
B(h,j)=-BIJ/k0;B(j,h)=-BIJ/k0;G(h,h)=G(h,h)+GIJ/k0/k0;B(h,h)=B(h,h)+BIJ/k0/k0;G(j,j)=G(j,j)+GIJ;B(j,j)=B(j,j)+BIJ+b;
end       
for k=1:numel(branchi)
    h=branchi(k);G(h,h)=G(h,h)+branchg(k);B(h,h)=B(h,h)+branchb(k);
end
Y0=G+i*B;Y=abs(Y0);C=angle(Y0);kmax=20;k=0;P=Pgb-Pdi;Q=Qgb-Qdi;P(a(1,2))=0;Q(a(1,2))=0;
for k=1:kmax
    m=0;n=0;Sij=zeros(a(1,1),a(1,1));AP=zeros(1,a(1,1));AQ=zeros(1,a(1,1));J=zeros(2*a(1,1),2*a(1,1));APQ=zeros(2*a(1,1),1);
    for m=1:a(1,1)
        for n=1:a(1,1)
            Sij(m,n)=S(m)-S(n)-C(m,n);                
            AP(m)=AP(m)-V(m)*Y(m,n)*V(n)*cos(Sij(m,n));AQ(m)=AQ(m)-V(m)*Y(m,n)*V(n)*sin(Sij(m,n));              
         end               
        if m~=a(1,2)
          AP(m)=AP(m)+P(m);AQ(m)=AQ(m)+Q(m);
        else
            AP(m)=0;AQ(m)=0;
        end
        for h=1:numel(pvi)
            AQ(pvi(h))=0;
        end
    end     
    if (max(AP)<a(1,4))&(max(AQ)<a(1,4))
        break;
    end
    for m=1:a(1,1)
        for n=1:a(1,1)
            if m~=n                
               J(2*m-1,2*n-1)=-V(m)*Y(m,n)*V(n)*sin(Sij(m,n));J(2*m-1,2*n)=-V(m)*Y(m,n)*V(n)*cos(Sij(m,n));J(2*m,2*n-1)=-J(2*m-1,2*n);
               J(2*m,2*n)=J(2*m-1,2*n-1);
            else
                for h=1:a(1,1)
                  J(2*m-1,2*m-1)=J(2*m-1,2*m-1)+V(m)*Y(m,h)*V(h)*sin(Sij(m,h));J(2*m-1,2*n)=J(2*m-1,2*n)-V(m)*Y(m,h)*V(h)*cos(Sij(m,h));
                  J(2*m,2*n-1)=J(2*m,2*n-1)-V(m)*Y(m,h)*V(h)*cos(Sij(m,h));J(2*m,2*n)=J(2*m,2*n)-V(m)*Y(m,h)*V(h)*sin(Sij(m,h));
                end
                  J(2*m-1,2*m-1)=J(2*m-1,2*m-1)-V(m)*Y(m,m)*V(m)*sin(Sij(m,m));J(2*m-1,2*n)=J(2*m-1,2*n)-V(m)*V(m)*Y(m,m)*cos(C(m,m));
                  J(2*m,2*n-1)=J(2*m,2*n-1)+V(m)*Y(m,m)*V(m)*cos(Sij(m,m));J(2*m,2*n)=J(2*m,2*n)+V(m)*V(m)*Y(m,m)*sin(C(m,m));
            end
        end
    end
    for n=1:a(1,1)
        APQ(2*n-1,1)=AP(n);APQ(2*n,1)=AQ(n);J(2*a(1,2)-1,2*n-1)=0;J(2*a(1,2)-1,2*n)=0;J(2*a(1,2),2*n-1)=0;J(2*a(1,2),2*n)=0;
        J(2*n-1,2*a(1,2)-1)=0;J(2*n,2*a(1,2)-1)=0;J(2*n-1,2*a(1,2))=0;J(2*n,2*a(1,2))=0;
    end
    J(2*a(1,2)-1,2*a(1,2)-1)=1;J(2*a(1,2),2*a(1,2))=1;
    for n=1:numel(pvi)
        m=pvi(n);
        for h=1:2*a(1,1)
            J(2*m,h)=0;J(h,2*m)=0;
        end
        J(2*m,2*m)=1;
    end 
    ASV=-J\APQ;
    for n=1:a(1,1)
        S(n)=S(n)+ASV(2*n-1);V(n)=V(n)+V(n)*ASV(2*n);
    end
end
cputime-t