PQ分解法计算潮流

一、PQ分解法的原理P-Q分解法是牛顿-拉夫逊法潮流计算的一种简化方法。P-Q分解法利用了电力系统的一些特有的运行特性，对牛顿-拉夫逊法做了简化，以改进和提高计算速度。的基本思想是根据电力系统实际运行特点：通常网络上的电抗远大于电阻，则系统母线电压幅值的微小变化对用功功率的改变影响很小。同样，母线电压相角的的改变对无功功率的影响较小。因此，节点功率方程在用极坐标形式表示时。它的修正方程式可简化为：00PHQLUU将P、Q分开来迭代计算，因此大大地减少了计算工作量。但是H、L在迭代过程中仍将不断变化，而且又都是不对称矩阵。对牛顿法的进一步简化。为把上式中的系数矩阵简化成迭代过程中不变的对称矩阵。在一般情况下线路两端的电压相角ij是不大的，因此可以认为：cos1sinijijijijGB=2iiiiQUB=考虑到上述关系，可以得到：ijiijjijiijjHUBULUBU节点的功率增量为：11(cossin)(sincos)niisijijijijijjniisijijijijijjPPUUGBQQUUGBP-Q分解法的特点：以一个n-1阶和一个n-m-1阶线性方程组代替原有的2n-m-1阶线性方程组；修正方程的系数矩阵B’和B”为对称常数矩阵，且在迭代过程中保持不变；P-Q分解法具有线性收敛特性，与牛顿-拉夫逊法相比，当收敛到同样的精度时需要的迭代次数较多。二、程序说明1.数据说明Branch1.txt:支路参数矩阵第1列为支路的首端编号；第2列为支路的末端编号（首端编号小于末端编号）；第3列为之路的阻抗；第4为支路的对地容抗；第5列为支路的变比；第6列为折算到那一侧的标志Branch2.txt:节点参数矩阵第1列为节点所接发电机的功率；第2列为节点负荷的功率；第3列为节点电压的初始值；第4列为PV节点的电压V给定值；第5列为节点所接的无功补偿设备的容量；第6列为节点分类标号igl，其中igl=1为平衡节点，igl=2为PQ节点，igl=3为PV节点。n为节点数、sb为平很母线节点号（固定为1）、pr为误差精度。2.运行将上述三个文档及主程序文件放在同一个路径下，程序运行时，按照相应的提示输入如下数据：n为节点数、sb为平很母线节点号、pr为误差精度，运行，得出结果。三、流程图开始输入原始数据形成矩阵B’和B”设置PQ节点电压初值，各节点电压相角的初值置迭代计数n=0计算不平衡有功功率max{}?kipP是置Kp=0Kq=0？是解修正方程，并修正后电压相角置Kq=1计算不平衡无功功率无功收敛？否解修正方程，求得修正后电压置Kp=1，k=k+1计算平衡节点功率及全部线路功率输出是否是置Kq=0Kp=0？是否Kp=1，Kq=1四、源程序及运行结果x=size(B2);n=x(1);m=length(find(B2(:,1)==1));y=size(B1);G=zeros(n);B=zeros(n);fori=1:y(1)QS=B1(i,1);%起始节点ZZ=B1(i,2);%终止节点R=B1(i,3);X=B1(i,4);BB=B1(i,5);K=1/B1(i,6);%变比RX2=R^2+X^2;G(QS,ZZ)=-K*R/RX2;G(ZZ,QS)=-K*R/RX2;B(QS,ZZ)=K*X/RX2;B(ZZ,QS)=K*X/RX2;G(QS,QS)=G(QS,QS)+R/RX2;B(QS,QS)=B(QS,QS)-X/RX2;B(QS,QS)=B(QS,QS)+BB/2;G(ZZ,ZZ)=G(ZZ,ZZ)+K^2*R/RX2;B(ZZ,ZZ)=B(ZZ,ZZ)-K^2*X/RX2;B(ZZ,ZZ)=B(ZZ,ZZ)+BB/2;end;fori=1:nB(i,i)=B(i,i)+B2(i,8);end;Delta_P=zeros(n,1);Delta_Q=zeros(n,1);Delta_V_Ample=zeros(n,1);Delta_V_Angle=zeros(n,1);PP=find(B2(:,1)3);%PQ/PV节点QQ=find(B2(:,1)==1);%PQ节点%V_Ample=B2(:,2);%V_Angle=B2(:,3)*pi/180;V_Ample=B2(:,2);V_Ample(QQ)=1.0;V_Angle=zeros(n,1);Pis=B2(:,4)/100;Qis=B2(:,5)/100;Pld=B2(:,6)/100;Qld=B2(:,7)/100;Pis(PP)=Pis(PP)-Pld(PP);Qis(QQ)=Qis(QQ)-Qld(QQ);kp=0;kq=0;ep=1e-7;%误差k=0;while(~(kp*kq))k=k+1;%----------------------------------------------------------------------Delta_V_Angle_ij=V_Angle*ones(1,n)-(V_Angle*ones(1,n))';%计算相角差矩阵S=G.*cos(Delta_V_Angle_ij)+B.*sin(Delta_V_Angle_ij);%Sij=Gij*cos(delta_ij)+Bij*sin(delta_ij)P=V_Ample.*(S*V_Ample);%P为计算出来的母线功率Delta_P(PP)=Pis(PP)-P(PP);%有功不平衡量%[Pis(PP),P(PP),Pis(PP)-P(PP),(Pis(PP)-P(PP))./Pis(PP)*100]ifmax(abs(Delta_P))ep%有功不平衡量小于允许范围就做标记kp=1;end;ifkp==0%有功不平衡量大于允许范围就作修正Delta_V_Angle(PP)=-(B(PP,PP)\(Delta_P(PP)./V_Ample(PP)))./V_Ample(PP);%解修正方程V_Angle(PP)=V_Angle(PP)+Delta_V_Angle(PP);%修正相角kq=0;%这是必须的，认为无功不平衡量还是不满足end;%----------------------------------------------------------------------Delta_V_Angle_ij=V_Angle*ones(1,n)-(V_Angle*ones(1,n))';%计算相角差矩阵T=G.*sin(Delta_V_Angle_ij)-B.*cos(Delta_V_Angle_ij);%Tij=Gij*sin(delta_ij)-Bij*cos(delta_ij)Q=V_Ample.*(T*V_Ample);%Q为计算出来的母线功率Delta_Q(QQ)=Qis(QQ)-Q(QQ);%无功不平衡量%[Qis(QQ),Q(QQ),Qis(QQ)-Q(QQ),(Qis(QQ)-Q(QQ))./Qis(QQ)*100]ifmax(abs(Delta_Q))ep%无功不平衡量小于允许范围就做标记kq=1;end;ifkq==0%无功不平衡量大于允许范围就作修正Delta_V_Ample(QQ)=-B(QQ,QQ)\(Delta_Q(QQ)./V_Ample(QQ));%解修正方程V_Ample(QQ)=V_Ample(QQ)+Delta_V_Ample(QQ);%修正幅值kp=0;%这是必须的，认为有功不平衡量还是不满足end;%----------------------------------------------------------------------[Delta_P(PP)',Delta_Q(QQ)';V_Ample(QQ)',V_Angle(PP)'*180/pi];%每次迭代结果end;[V_Ample';B2(:,2)';V_Ample'-B2(:,2)';V_Angle'*180/pi;B2(:,3)';V_Angle'*180/pi-B2(:,3)']';P=P+Pld;Q=Q+Qld;V=[V_Ample,V_Angle*180/pi,P,Q]%母线全部信息，幅值、相角、发电有功、发电无功Sij_flow=zeros(y(1),6);%y(1)为支路总数,Sij_flow为所有支路上流动的功率Sij_loss=zeros(y(1),4);%y(1)为支路总数,Sij_loss为所有支路上损耗的功率fori=1:y(1)ii=B1(i,1);jj=B1(i,2);BB=B1(i,5);kk=B1(i,6);Vi=V_Ample(ii);Vj=V_Ample(jj);ag=Delta_V_Angle_ij(ii,jj);%ag为i,j节点相角差Sij_flow(i,1)=ii;Sij_flow(i,2)=jj;Sij_flow(i,3)=-G(ii,jj)*(Vj^2-Vi*Vj*cos(ag))+Vi*Vj*B(ii,jj)*sin(-ag);Sij_flow(i,4)=B(ii,jj)*(Vj^2-Vi*Vj*cos(ag))+Vi*Vj*G(ii,jj)*sin(-ag)-Vj^2*BB/2;Sij_flow(i,5)=G(ii,jj)*(Vi^2-Vi*Vj*cos(ag))-Vi*Vj*B(ii,jj)*sin(ag);Sij_flow(i,6)=-B(ii,jj)*(Vi^2-Vi*Vj*cos(ag))-Vi*Vj*G(ii,jj)*sin(ag)+Vi^2*BB/2;Sij_loss(i,1)=ii;Sij_loss(i,2)=jj;Sij_loss(i,3)=Sij_flow(i,3)-Sij_flow(i,5);Sij_loss(i,4)=Sij_flow(i,4)-Sij_flow(i,6);end;S=[Sij_flow,Sij_loss(:,3:4)]%全部线路的功率PLoss_Sum_Line=sum(S(:,7));%合计每条线路上的有功损耗求网损（S阵说明：第一、二列为节点序号，第三四列和第五六列分别为支路的首端和末端功率，第七列为该支路有功损耗，第八列为该支路无功损耗）PLoss_Sum_GeLd=sum(P)-sum(Pld)%总的母线注入功率和总的负荷功率之差求网损fid=fopen('潮流计算结果.txt','wt');fprintf(fid,'节点结果:');fprintf(fid,'\n');fprintf(fid,'电压幅值相角发电有功发电无功');fprintf(fid,'\n');fprintf(fid,'%f\t%f\t%f\t%f\n',V');fprintf(fid,'线路功率结果:');fprintf(fid,'\n');fprintf(fid,'节点i节点jPijQijPjiQjiLine_P_LossLine_Q_Loss');fprintf(fid,'\n');fprintf(fid,'%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n',S');fprintf(fid,'网损PLoss_Sum_Geld=%f',PLoss_Sum_GeLd);fclose(fid)运行结果节点结果:电压幅值相角发电有功发电无功1.0600000.0000002.323958-0.1688591.045000-4.9809750.4000000.4239941.010000-12.7176000.0000000.2339231.018626-10.323542-0.000000-0.0000001.020255-8.7836450.0000000.0000001.070000-14.235716-0.0000000.1217511.061947-13.361635-0.000000-0.0000001.090000-13.3616350.0000000.1