第四节PQ分解法潮流计算

129第四节PQ分解法潮流计算一、PQ分解法的基本方程式60年代以来N—R法曾经是潮流计算中应用比较普遍的方法，但随着网络规模的扩大（从计算几十个节点增加到几百个甚至上千个节点）以及计算机从离线计算向在线计算的发展，N—R法在内存需要量及计算速度方面越来越不适应要求。70年代中期出现的快速分解法比较成功的解决了上述问题，使潮流计算在N—R法的基础上向前迈进了一大步，成为取代N—R法的算法之一。快速分解法（又称P—Q分解法）是从简化牛顿法极坐标形式计算潮流程序的基础上提出来的。它的基本思想是根据电力系统实际运行特点：通常网络上的电抗远大于电阻值，则系统母线电压副值的微小变化V对母线有功功率的改变P影响很小。同样，母线电压相角的少许改变，也不会引起母线无功功率的明显改变Q。因此，节点功率方程在用极坐标形式表示时，它的修正方程式可简化为：VVLHQP/00（4—19）这就是把2（n—1）阶的线性方程组变成了两个n—1阶的线性方程组，将P和Q分开来进行迭代计算，因而大大地减少了计算工作量。但是，H，L在迭代过程中仍然在不断的变化，而且又都是不对称的矩阵。对牛顿法的进一步简化（也是最关键的一步），即把（4—19）中的系数矩阵简化为在迭代过程中不变的对称矩阵。在一般情况下，线路两端电压的相角ij是不大的（不超过10○～20○）。因此，可以认为：ijijijijBGsin1cos（4—20）此外，与系统各节点无功功率相应的导纳BLDi远远小于该节点自导纳的虚部，即130iiiiLDiBVQB2因而iiiiBVQ2（4—21）考虑到以上关系，式（4—19）的系数矩阵中的各元素可表示为：ijjiijBVVH（i,j=1,2,………，n-1）（4—22）ijjiijBVVL（i,j=1,2,……………，m）(4—23)而系数矩阵H和L则可以分别写成：11,1122,1111,1111,222222121211,1121211111nnnnnnnnnnnnVBVVBVVBVVBVVBVVBVVBVVBVVBVH=1211,12,11,11,222211,11211121nnnnnnnnVVVBBBBBBBBBVVV=11DDBVV（4—24）mmmmmmmmmmmmmVBVVBVVBVVBVVBVVBVVBVVBVVBVL22122222212121121211111=mmmmmmmmVVVBBBBBBBBBVVV2121222211121121=22''DDVBV（4—25）将（4—24）和（4—25）式代入（4—19）中，得到13111DDVBVPVBVQD''2用11DV和12DV分别左乘以上两式便得：111'DDVBPV（4—26）VBQVD''12（4—27）这就是简化了的修正方程式，它们也可展开写成：1122111,12,11,11,222211,11211112211nnnnnnnnnnVVVBBBBBBBBBVPVPVP（4—28）mmmmmmmmmVVVBBBBBBBBBVQVQVQ212122221112112211（4—29）在这两个修正方程式中系数矩阵元素就是系统导纳矩阵的虚部，因而系数矩阵是对称矩阵，且在迭代过程中保持不变。这就大大减少了计算工作量。用极坐标表示的节点功率增量为：njijijijijjiisinjijijijijjiisiBGVVQQBGVVPP110)cossin(0)sincos(（4—30）式（4—28）、（4—29）和（4—30）构成了P—Q分解法迭代过程的基本方程式。132二、计算步骤和程序框图（1）给定各节点电压的初始值)0()0(,iiV；（2）代入式（4—30）计算各节点有功功率iP，并求出iiVQ/；（3）解修正方程式（4—28），得出各节点电压相角修正量i；（4）修正各节点电压的相角i)()()1(kikiki（5）式（4—30）求得各节点无功功率误差iQ，并求出iiVQ/（6）求解修正方程式（4—29），得出各节点电压幅值的修正量iV；（7）修正各节点电压的幅值iV，)()()1(kikikiVVV（8）回（2）进行迭代，直到各节点功率误差iP及iQ都满足收敛条件P—Q分解法程序框图：一、程序清单及打印结果是置KQ=0KP=0?是是Max?}{)(QkiQ否否解修正方程（4—28）求)(ki用公式（4—30）计算不平衡功率)(kiP，计算)()(/kikiVP输入原始数据形成矩阵B’和B’’并进行三角分解设PQ节点电压初值，各节点电压相角初值置迭代计数k=0Kp=1,kQ=1?}max{')(pkiP)()()1(kikiki置1QK用公式（4—30）计算不平衡功率)(kiQ，计算)()(/kikiVQ置Kp=0KQ=0?否是133%本程序的功能是用P-Q分解法进行潮流计算n=input('请输入节点数:n=');nl=input('请输入支路数:nl=');isb=input('请输入平衡母线节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入由支路参数形成的矩阵:B1=');B2=input('请输入由节点参数形成的矩阵:B2=');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);fori=1:nlifB1(i,6)==0p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);end134Y(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);fori=1:nS(i)=B2(i,1)-B2(i,2);BI(i,i)=BI(i,i)+B2(i,5);endP=real(S);Q=imag(S);fori=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfori=1:nifB2(i,6)==2V(i)=sqrt(e(i)^2+f(i)^2);O(i)=atan(f(i)./e(i));endendfori=2:nifi==nB(i,i)=1./B(i,i);135elseIC1=i+1;forj1=IC1:nB(i,j1)=B(i,j1)./B(i,i);endB(i,i)=1./B(i,i);fork=i+1:nforj1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendendp=0;q=0;fori=1:nifB2(i,6)==2p=p+1;k=0;forj1=1:nifB2(j1,6)==2k=k+1;A(p,k)=BI(i,j1);endendendendfori=1:naifi==naA(i,i)=1./A(i,i);elsek=i+1;forj1=k:na136A(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);fork=i+1:naforj1=i+1:naA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendendICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;whileICT2~=0|ICT3~=0ICT2=0;ICT3=0;fori=1:nifi~=isbC(i)=0;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));endDP1(i)=P(i)-V(i)*C(i);DP(i)=DP1(i)./V(i);DET=abs(DP1(i));ifDET=prICT2=ICT2+1;endendendNp(K)=ICT2;ifICT2~=0137fori=2:nDP(i)=B(i,i)*DP(i);ifi~=nIC1=i+1;fork=IC1:nDP(k)=DP(k)-B(k,i)*DP(i);endelseforLZ=3:iL=i+3-LZ;IC4=L-1;forMZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);endendendendfori=2:nO(i)=O(i)-DP(i);endkq=1;L=0;fori=1:nifB2(i,6)==2C(i)=0;L=L+1;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);138DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));ifDET=prICT3=ICT3+1;endendendelsekp=0;ifkq~=0;L=0;fori=1:nifB2(i,6)==2C(i)=0;L=L+1;fork=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));endendendendNq(K)=ICT3;ifICT3~=0L=0;fori=1:naDQ(i)=A(i,i)*DQ(i);ifi==na139forLZ=2:iL=i+2-LZ;IC4=L-1;forMZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIC1=i+1;fork=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;fori=1:nifB2(i,6)==2L=L+1;V(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;ifkp~=0K=K+1;end140endenddisp('迭代次数');disp(K);disp('每次没有达到精度要求的有功功率个数为');disp(Np);disp('每次没有达到精度要求的无功功率个数为');disp(Nq);fork=1:nE(k)=V(k)*cos(O(k))+V