原程序代码.txt

这个是电力系统的很经典的潮流算法

x1=fopen('中途1.txt');                                           
y1=fgetl(x1);
fclose(x1);
z1=str2num(y1);
n=z1(1);           
nl=z1(2);
O=z1(3);
pr=z1(4);
x2=fopen('中途2.txt');
y2=fgetl(x2);
fclose(x2);
B1=str2num(y2);                                                        
x3=fopen('中途3.txt');
y3=fgetl(x3);
fclose(x3);
B2=str2num(y3);
Y=zeros(n);                                                             
e=zeros(1,n);                                                           
f=zeros(1,n);                                                           
V=zeros(1,n);   
O=zeros(1,n);
S1=zeros(nl);
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));  
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;                   
    Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;                                
end
%以上求导纳矩阵               

G=real(Y);
B=imag(Y);
for i=1:n
    e(i)=real(B2(i,3));                                                  
    f(i)=imag(B2(i,3));
    V(i)=B2(i,4);                                                        
end
for i=1:n
    S(i)=B2(i,1)-B2(i,2);                                                
    B(i,i)=B(i,i)+B2(i,5);                                              
end
P=real(S);Q=imag(S);                                                      
 IT2=1; N0=2*n; N=N0+1; a=0;
while IT2~=0
    IT2=0; a=a+1;
    for i=1:n
        if i~=O                                                         
            C(i)=0;
            D(i)=0;
            for j1=1:n
                C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);
                D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);
            end
            P1=C(i)*e(i)+f(i)*D(i);
            Q1=f(i)*C(i)-D(i)*e(i);         
 %以上求P',Q'
            V2=e(i)^2+f(i)^2;   
            if B2(i,6)~=3                                                
                DP=P(i)-P1;            
                DQ=Q(i)-Q1;
                for j1=1:n
                    if j1~=O & j1~=i                                    
                        X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
                        X2=B(i,j1)*e(i)-G(i,j1)*f(i);
                        X3=X2;
                        X4=-X1;                                          
                        p=2*i-1;
                        q=2*j1-1;
                        J(p,q)=X3;
                        J(p,N)=DQ;                                        
                        m=p+1;
                        J(m,q)=X1;
                        J(m,N)=DP;
                        q=q+1;
                        J(p,q)=X4;
                        J(m,q)=X2;
                    elseif j1==i&j1~=O
                        X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
                        X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                        X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                        X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
                        p=2*i-1;
                        q=2*j1-1;
                        J(p,q)=X3;
                        J(p,N)=DQ;
                        m=p+1;
                        J(m,q)=X1;
                        J(m,N)=DP;
                        q=q+1;
                        J(p,q)=X4;
                        J(m,q)=X2;
                    end
                end
            else
                DP=P(i)-P1;
                DV=V(i)^2-V2;
                for j1=1:n
                    if j1~=O & j1~=i
                        X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
                        X2=B(i,j1)*e(i)-G(i,j1)*f(i);
                        X5=0;
                        X6=0;
                        p=2*i-1;
                        q=2*j1-1;
                        J(p,q)=X5;
                        J(p,N)=DV;
                        m=p+1;
                        J(m,q)=X1;
                        J(m,N)=DP;
                        q=q+1;
                        J(p,q)=X6;
                        J(m,q)=X2;
                    elseif j1==i&j1~=O
                        X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
                        X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                        X5=-2*e(i);
                        X6=-2*f(i);
                        p=2*i-1;
                        q=2*j1-1;
                        J(p,q)=X5;
                        J(p,N)=DV;
                        m=p+1;
                        J(m,q)=X1;
                        J(m,N)=DP;
                        q=q+1;
                        J(p,q)=X6;
                        J(m,q)=X2; 
                    end
                end
            end
        end
    end
    %以上求雅克比矩阵
    for k=3:N0
        k1=k+1; N1=N;
        for k2=k1:N1
            J(k,k2)=J(k,k2)./J(k,k);
        end
        J(k,k)=1;
        if k~=3
            k4=k-1;
            for k3=3:k4
                for k2=k1:N1
                    J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
                end
                J(k3,k)=0;
            end
            if k==N0,break; 
         end
            for k3=k1:N0
                for k2=k1:N1
                    J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
                end
                J(k3,k)=0;
            end
            else
                for k3=k1:N0
                    for k2=k1:N1
                        J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
                    end
                    J(k3,k)=0;
                end
            end
        end
        for k=3:2:N0-1
            L=(k+1)./2;
            e(L)=e(L)-J(k,N);
            k1=k+1;
            f(L)=f(L)-J(k1,N);
        end
        for k=3:N0
            DET=abs(J(k,N));     
            if DET>=pr
                IT2=IT2+1;
            end
        end
    for k=1:n
        V(k)=sqrt(e(k)^2+f(k)^2);   
        Z(k)=atan(f(k)./e(k))*180./pi;  
        E(k)=e(k)+f(k)*j;                                          
    end
    disp('求导纳矩阵为');
    disp(Y);
    disp('各节点的实际电压标么值为（节点号从小到大）');
    disp(E)
    disp('各节点的电压大小为（节点号从小到大）');
    disp(V);
    disp('各节点的电压相角为（节点号从小到大）');
    disp(Z);
    for p=1:n
        C(p)=0;
        for q=1:n
            C(p)=C(p)+conj(Y(p,q))*conj(E(q)); 
        end
        S(p)=E(p)*C(p);       
    end
    disp('各节点的功率为');
    disp(S);
    for i=1:nl               
        if B1(i,6)==0
            p=B1(i,1);
            q=B1(i,2);
        else p=B1(i,2);
             q=B1(i,1);
         end
         Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));                                                                 
     end
    for i=1:nl               
        if B1(i,6)==0
            p=B1(i,1);
            q=B1(i,2);
        else p=B1(i,2);
             q=B1(i,1);
         end
         Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
     end
    disp('各支路的功率损耗为');
      for i=1:nl               
        if B1(i,6)==0
            p=B1(i,1);
            q=B1(i,2);
        else p=B1(i,2);
             q=B1(i,1);
         end
         DS(i)=Si(p,q)+Sj(q,p);
         disp (DS(i));
     end
 end