pq.m

这是一个对电力系统潮流计算用matlab语言来实现pq法的例子。

%PQ分解迭代法
%按何仰赞书p46页例题编程
clear;

 N=input('请输入节点数:n=');
 m=input('请输入PV节点数:m=');
 nl=input('请输入支路数:nl=');
 ph=input('请输入平衡母线节点号ph=');
 jd=input('请输入误差精度jd=');
 % B=input('请输入由支路参数形成的矩阵:B=');
    
A=[-0.30-0.18i    1.0      0   2;             %1号PQ节点
    -0.55-0.13i   1.0      0   2;             %2号PQ节点
      0.5         1.1     1.1  3;              %3号PV节点
      0           1.05     0  1;]             %4号平衡节点 
  % A=input('请输入各节点参数形成的矩阵:A=');


Y=zeros(N);e=zeros(1,N);f=zeros(1,N);V=zeros(1,N);

for i=1:nl
     p=B(i,1);q=B(i,2);
    Y(p,q)=Y(p,q)-1./(B(i,3)*B(i,5));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(B(i,3)*B(i,5)^2)+B(i,4);
    Y(p,p)=Y(p,p)+1./B(i,3)+B(i,4);
end  %形成节点导纳矩阵          
disp('节点导纳矩阵Y为:');
disp(Y);

S=[-0.30-j*0.18        -0.55-j*0.13          0.5       0       ];%假定PV节点的Q初值为0.05
V=[1                   1                     1.1            1.05    ];%假定PQ节点的V初值为1
Vs=[1                  1                     1.1            1.05    ];
alpha=[0               0                     0              0];

G=real(Y);
B=imag(Y);
Ps=real(S);
Qs=imag(S);

n=N-m-1;%PQ节点数目
P=zeros(N-1);
Q=zeros(N-1);
DeltVP=zeros(N-1,1);
DeltVQ=zeros(n,1);
DeltP=zeros(N-1,1);
DeltQ=zeros(n,1);
B1=zeros(N-1);   % B'
B2=zeros(m);     % B''

B1=B(1:(N-1),1:(N-1));
B2=B(1:n,1:n);    
kp=1;
kq=1;
for k=0:100
    %计算PQ,PV节点的P误差    
    kp=1;
    MAX1=0;
    for i=1:(N-1)
        sum1=0;
        for j=1:N
            tmp1=V(j)*( G(i,j)*cos(alpha(i)-alpha(j)) + B(i,j)*sin(alpha(i)-alpha(j)) );
            sum1=sum1+tmp1;
        end
        P(i)=V(i)*sum1;
        DeltP(i)=Ps(i)-P(i);%DeltP      
        
        if(MAX1<abs(DeltP(i))) MAX1=abs(DeltP(i));
        end
        DeltP(i)=DeltP(i)/V(i);  % 12-68'left
    end
    
    %判断是否收敛,e=0.0001
    e=1E-5;   
    if(MAX1<e)  
        kp=0;
        if(kq==0) break;
        end
    else
        DeltVP=-inv(B1)*DeltP;
        for i=1:(N-1)
            alpha(i)=alpha(i)+DeltVP(i)/V(i);
        end
    end   
    
    %计算PQ节点的Q误差
    kq=1;
    MAX2=0;
    for i=1:n
        sum1=0;
        for j=1:N
            tmp1=V(j)*( G(i,j)*sin(alpha(i)-alpha(j)) - B(i,j)*cos(alpha(i)-alpha(j)) );
            sum1=sum1+tmp1;
        end
        Q(i)=V(i)*sum1;
        DeltQ(i)=Qs(i)-Q(i);%DeltQ
        
        if(MAX2<abs(DeltQ(i))) MAX2=abs(DeltQ(i));
        end
        DeltQ(i)=DeltQ(i)/V(i);
    end
    
    %判断是否收敛,e=0.0001
    e=1E-5;   
    if(MAX2<e)  
        kq=0;
        if(kp==0) break;
        end
    else
        DeltVQ=-inv(B2)*DeltQ;
        for i=1:n
            V(i)=V(i)+DeltVQ(i);
        end
    end 
end
k
MAX1
MAX2
V
alpha*180/pi