pq.m

基于P-Q分解法的潮流计算适合初学者学习程序简单易懂

%本程序的功能是利用PQ分解法进行潮流计算
clear
clc
n=input('请输入节点数：n=');
nl=input('请输入支路数：nl=');
isb=input('请输入平衡母线节点号：isb=');
pr=input('请输入误差精度：pr=');
disp('请输入由支路参数形成的矩阵：B1=');
%输入B1
[fname, pname]= uigetfile('*.xls'); 
file=[pname fname]; 
x=xlsread(file);
B1=[x(:,1),x(:,2),(x(:,3)+x(:,4)*i),x(:,5)*i,x(:,6),x(:,7)];
disp('请输入由节点参数形成的矩阵：B2=');
%输入B2
[fname, pname]= uigetfile('*.xls'); 
file=[pname fname]; 
y=xlsread(file);
B2=[(y(:,1)+y(:,2)*i),(y(:,3)+y(:,4)*i),y(:,5),y(:,6),y(:,7),y(:,8)];
disp('请输入由节点号及其对地阻抗形成的矩阵：X=');
%输入X
[fname, pname]= uigetfile('*.xls'); 
file=[pname fname]; 
z=xlsread(file);
X=z;
na=input('请输入PQ节点数na=');
Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);
Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);
for i=1:n
    if X(i,2)~=0;
        p=X(i,1);
        Y(p,p)=1./X(i,2);
    end
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
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
    YI(p,q)=YI(p,q)-1./B1(i,3);
    Y(q,p)=Y(p,q);
    YI(q,p)=YI(p,q);
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
    YI(q,q)=YI(q,q)+1./B1(i,3);
    Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
    YI(p,p)=YI(p,p)+1./B1(i,3);
end %求导纳矩阵

G=real(Y);B=imag(YI);BI=imag(Y);
for i=1:n
    S(i)=B2(i,1)-B2(i,2);
    BI(i,i)=BI(i,i)+B2(i,5);
end
P=real(S);Q=imag(S);
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
    if B2(i,6)==2
        V(i)=sqrt(e(i)^2+f(i)^2);
        O(i)=atan(f(i)./e(i));
    end
end
for i=2:n
    if i==n
        B(i,i)=1./B(i,i);
    else IC1=i+1;
        for j1=IC1:n
            B(i,j1)=B(i,j1)./B(i,i);
        end
        B(i,i)=1./B(i,i);
        for k=i+1:n
            for j1=i+1:n
                B(k,j1)=B(k,j1)-B(k,i)*B(i,j1);
            end
        end
    end
end
p=0;q=0;
for i=1:n
    if B2(i,6)==2
        p=p+1;k=0;
        for j1=1:n
            if B2(j1,6)==2
                k=k+1;
                A(p,k)=BI(i,j1);
            end
        end
    end
end
for i=1:na
    if i==na
        A(i,i)=1./A(i,i);
    else k=i+1;
        for j1=k:na
            A(i,j1)=A(i,j1)./A(i,i);
        end
        A(i,i)=1./A(i,i);
        for k=i+1:na
            for j1=i+1:na
                A(k,j1)=A(k,j1)-A(k,i)*A(i,j1);
            end
        end
    end
end
ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;
while ICT2~=0|ICT3~=0
    ICT2=0;ICT3=0;
    for i=1:n
        if i~=isb
            C(i)=0;
            for k=1:n
                C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));
            end
            DP1(i)=P(i)-V(i)*C(i);
            DP(i)=DP1(i)./V(i);
            DET=abs(DP1(i));
            if DET>=pr
                ICT2=ICT2+1;
            end
        end
        end
    Np(K)=ICT2;
    if ICT2~=0
        for i=2:n
            DP(i)=B(i,i)*DP(i);
            if i~=n
                IC1=i+1;
                for k=IC1:n
                    DP(k)=DP(k)-B(k,i)*DP(i);
                end
            else 
                for LZ=3:i
                    L=i+3-LZ;
                    IC4=L-1;
                    for MZ=2:IC4
                    I=IC4+2-MZ;
                    DP(I)=DP(I)-B(I,L)*DP(L);
                    end
                end
            end
        end
        for i=2:n
            O(i)=O(i)-DP(i);
        end
        kq=1;L=0;
        for i=1:n
            if B2(i,6)==2
                C(i)=0;L=L+1;
                for k=1:n
                    C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));
                end
                DQ1(i)=Q(i)-V(i)*C(i);
                DQ(L)=DQ1(i)./V(i);
                DET=abs(DQ1(i));
                if DET >=pr
                    ICT3=ICT3+1;
                end
            end
        end
    else kp=0;
        if kq~=0;
            L=0;
            for i=1:n
                if B2(i,6)==2
                    C(i)=0;L=L+1;
                    for k=1:n
                        C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));
                    end
                    DQ1(i)=Q(i)-V(i)*C(i);
                    DQ(L)=DQ1(i)./V(i);
                    DET=abs(DQ1(i));
                end
            end
        end
    end
    Nq(K)=ICT3;
    if ICT3~=0
        L=0;
        for i=1:na
            DQ(i)=A(i,i)*DQ(i);
            if i==na
                for LZ=2:i
                    L=i+2-LZ;
                    IC4=L-1;
                    for MZ=1:IC4
                        I=IC4+1-MZ;
                        DQ(I)=DQ(I)-A(I,L)*DQ(L);
                    end
                end
            else
                IC1=i+1;
                for k=IC1:na
                    DQ(k)=DQ(k)-A(k,i)*DQ(i);
                end
            end
        end
        L=0;
        for i=1:n
            if B2(i,6)==2
                L=L+1;
                V(i)=V(i)-DQ(L);
            end
        end
        kp=1;
        K=K+1;
    else 
        kq=0;
        if kp~=0
           K=K+1;
        end
    end
    for i=1:n
        Dy(K-1,i)=V(i);
    end
end
disp('迭代次数');
disp(K);
disp('每次没有达到精度要求的有功功率个数为');
disp(Np);
disp('每次没有达到精度要求的无功功率个数为');
disp(Nq);
for k=1:n
    E(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;
    O(k)=O(k)*180./pi;
end
disp('各节点的电压标幺值E为（节点号从小到大排）：');
disp(E);
disp('各节点的电压V大小（节点号从小到大排）为：');
disp(V);
disp('各节点的电压相角O（节点号从小到大排）为：');
disp(O);
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('各节点的功率S（节点号从小到大排）为：');
disp(S);
disp('各条支路的首端功率Si（顺序同您输入B1时一样）为：');
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))));
    disp(Si(p,q));
end
disp('各条支路的末端功率Sj（顺序同您输入B1时一样）为：');
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))));
    disp(Sj(q,p));
end
disp('各条支路的功率损耗DS（顺序同您输入B1时一样）为：');
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
for i=1:K
    Cs(i)=i;
    for j=1:n
        Dy(K,j)=Dy(K-1,j);
    end
end
disp('以下是每次迭代后各节点的电压值（如图所示）');
plot(Cs,Dy),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');

save C:\mydata K E V O S Si Sj DS Y