牛顿-拉夫逊迭代法电力网潮流计算方法与程序

1、牛顿-拉夫逊迭代法电力网潮流计算方法与程序编写 佘名寰 牛顿-拉夫逊迭代法在电力网潮流计算中因其收敛性较好获得广泛运用，该算法的难点是需反复计算功率方程中雅可比矩阵各个元素表达式。本文简叙了牛顿-拉夫逊迭代法潮流计算的基本公式，通过例题介绍了用牛顿-拉夫逊法计算电力网潮流电压的MATLAB程序。程序采用MATLAB语言的符号矩阵简化了雅可比矩阵系数的计算。本文可供电力系统电气技术人员和大专院校电力类专业师生参考。2牛顿-拉夫逊迭代法潮流计算2.1 极坐标型式节点功率方程式由电源注入或从负载流出节点的电流统称节点电流，由节点电压和节点电流所求得的节点注入功率为：Sti=Vi IiIi=YijVj

2、 (I,j=1, 2, .n)Pti=Vi 1nVjGijcosij+Bijsinij Qti=Vi 1nVjGijsinij-Bijcosij (i,j=1, 2.n) （2-01）式中 Pti, QtI节点注入的有功和无功功率 Vi，Vj节点 i, j电压幅值ij=i-j 节点 i, j电压的相角差Gij, Bij 节点导纳矩阵的元素，YIJ=GIJ+jBIJ节点功率平衡关系为：Pgi-PlI= PtiQgi-Qli= QtiPi=Pgi-PlI-Pti=0Qi=Qgi-Qli-Qti=0 （2-02）Pgi, Qgi节点i发电机输入有功和无功功率Pli , Qli节点i负荷有功和无功功率

3、Pi，Qi-节点i不平衡功率不平衡功率的微分 d(Pi), d(Qi) dpi=-(ptiv1v1+ ptiv2v2+pti11+pti22.) （i=1,2,n） dqi=-(qtiv1v1+ qv2v2+qti11+qti22.) （2-03）对于n个节点系统可得如下矩阵形式修正方程式；P1P2:PnQ1Q2:Qn= P11P21:Pn1Q11Q21:Qn1 P12P22:Pn2Q12Q22:Qn2 P1nP2n:PnnQ1nQ2n:Qnn P1v1P2v1:Pnv1Q1v1Q2v1:Qnv1 P1v2P2v2:Pnv2Q1v2Q2v2:Qnv2 P1vnP2vn:PnvnQ1vnQ2vn

4、:Qnvn 12:nv1v2:vn （2-04）式中偏微分矩阵为雅可比矩阵。若系统总共有n 个节点，其中r个为PV节点，去掉平衡节点则电压角度变量有n-1-个，电压幅值v变量有n-1-r个，方程式共有2n-r-2个。2.2 牛顿-拉夫逊迭代法解算过程； 指定各节点电压的初值；v1,v2,vn; 1,2,n.一般设V=1,=0；将所给定的电压值代入式（2-01），（2-02）计算各节点不平衡功率；将电压初值代入求雅可比矩阵系数；解方程式（2-04），求节点电压幅值和相角的修正值 vi,i;用所求得的电压修正值修正各节点电压： v(1)=v(0)+v(0) (1)=(0)+(0)用新的节点电压重复

5、以上五个步骤计算直到计算精度合符要求。 牛顿-拉夫逊迭代法计算的难点是求雅可比矩阵系数，本文采用MATLAB语言符号矩阵运算jacobian(f,x)命令求解函数f(x)对x的雅可比矩阵。2.3牛顿-拉夫逊迭代法计算程序 我们通过例题说明牛顿-拉夫逊迭代法计算潮流程序。【例2.3.1】 网络接线如图2-1所示。支路阻抗分别为Z12=j0.1,Z23=j0.1,Z31=j0.1;三条支路两端的对地电纳皆是j0.01。给定的注入功率分别为：S1=-1.7192-j0.7346;S2=0.6661。节点1是负荷节点即PQ节点；节点2给定的电压大小为V2=1.05,其属于PV节点；节点3是平衡节点，其

6、电压保持定值，大小为V3=1.0。Z13Z12Z23YC3YC1YC1V3=1.0312YC2P1=1.7192Q1=0.7346P2=0.6661V2=1.05Z21=j0.1Z23=j0.1Z13=j0.1YC1=J0.02YC2=J0.02YC3=J0.02 图2-1 例2.3.1网络接线图2.3.1.1 程序参变量说明 Y=G+j*B 节点导纳矩阵，G 为电导，B 为电纳； P 节点有功功率给定值，正为发电机注入，负为负荷功率； q 节点无功功率给定值，正为发电机注入，负为负荷功率； u 节点电压幅值的符号矩阵； delt节点电压相角的符号矩阵； pt 节点有功功率符号表达式； qt

7、节点无功功率符号表达式； pp节点有功功率不平衡值符号表达式； qq节点无功功率不平衡值符号表达式； uu 节点电压幅值数值矩阵； dd 节点电压相角数值矩阵； PP 节点不平衡功率的数值矩阵 N1 网络独立节点总数； N2 网络PV 节点总数； Sm 节点功率矩阵； Smn 支路功率矩阵； J1, J2, J 节点不平衡功率雅可比符号矩阵. JJ 节点不平衡功率雅可比数值矩阵2.3.1.2 源程序源程序如下：%* NUE78.m example 2-3-1 fig 2-1 * %The following Program for load flow calculation is based

8、on MATLAB7.0clearu=sym(u1,u2,u3);delt=sym(d1,d2,d3);G=zeros(3);B=-19.98,10,10;10,-19.98,10;10,10,-19.98;Y=G+j*B;p(1)=-1.7192;q(1)=-0.7346;p(2)=0.6661;k=0;precision=1;N1=2;%the N1 is the amount of the PQ and PV busN2=1;%the N2 is the amount of the PV bus for m=1:N1 for n=1:N1+1 pt(n)=u(m)*u(n)*(G(m,n

9、)*cos(delt(m)-delt(n)+B(m,n)*sin(delt(m)-delt(n); end pp(m)=p(m)-sum(pt); end for m=1:N1-N2 for n=1:N1+1 qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n)-B(m,n)*cos(delt(m)-delt(n); end qq(m)=q(m)-sum(qt); end J1=jacobian(pp,d1,d2,u1); J2=jacobian(qq,d1,d2,u1); J=vertcat(J1,J2); uu=1.0,1.05,1.0;dd=0,0,0;

10、 while precision0.00001 u1=uu(1);u2=uu(2);u3=uu(3); d1=dd(1);d2=dd(2);d3=dd(3); for m=1:N1 PP(m)=eval(pp(m); end for m=1:N1-N2 PP(N1+m)=eval(qq(m); end JJ=eval(J); du=-inv(JJ)*PP;precision=max(abs(du); for n=1:N1 dd(n)=dd(n)+du(n); endfor n=1:N1-N2 uu(n)=uu(n)+du(N1+n); end k=k+1;endk-1,dd,uu%*% the

11、 following program is used to calculate the Sm and Smnfor n=1:N1+1 U(n)=uu(n)*(cos(dd(n)+j*sin(dd(n);end Um=conj(U); I=Y*Um;Sm=diag(Um)*conj(I)for m=1:N1+1 for n=1:N1+1 Smn(m,n)=U(m)*(conj(U(m)-conj(U(n)*conj(-Y(m,n); endend2.3.1.3 程序说明节点导纳矩阵：Y11=Y22=Y33=0+j(-10-10+0.02)=0-j19.98 Y12=Y23=Y13=0+j10.0

12、节点电压符号矩阵：一共三个节点，独立节点数N1=2，PV节点数N2=1 设节点电压幅值符号矩阵为 u=sym(u1,u2,u3); 节点电压相角符号矩阵为 delt=sym(d1,d2,d3); p(1)，q(1)，p(2)为节点功率给定值，节点2为PV节点，q(2)待求； pt(n),pp(m)相关语句是按公式(2-01),(2-02)建立Pi符号表达式；P1=pp(1)=-2149/1250-10*u1*u2*sin(d1-d2)-10*u1*u3*sin(d1-d3);P2=pp(2)= 6661/10000+10*u1*u2*sin(d1-d2)-10*u2*u3*sin(d2-d3)

13、 qt(n),qq(m)相关语句是建立Qi符号表达式，只针对PQ节点；Q1=qq(1)= -3673/5000-999/50*u12+10*u1*u2*cos(d1-d2)+10*u1*u3*cos(d1-d3) J1，J2，J 语句是建立公式（2-03）中的雅可比符号矩阵；J = -10*u1*u2*cos(d1-d2)-10*u1*u3*cos(d1-d3), 10*u1*u2*cos(d1-d2), -10*u2*sin(d1-d2)-10*u3*sin(d1-d3) 10*u1*u2*cos(d1-d2), -10*u1*u2*cos(d1-d2)-10*u2*u3*cos(d2-d3

14、), 10*u2*sin(d1-d2) -10*u1*u2*sin(d1-d2)-10*u1*u3*sin(d1-d3), 10*u1*u2*sin(d1-d2), -999/25*u1+10*u2*cos(d1-d2)+10*u3*cos(d1-d3) while 循环语句是将所假设的节点电压幅值和相角代入符号矩阵，求出节点不平衡功率和雅可比数值矩阵 PP, JJ ，解方程（2-04）,求出节点电压幅值和相角的修正值vi,i,反复迭代，直到精度合符要求。 Um, Smn 语句计算节点注入功率和支路功率。2.3.1.4 程序计算结果：各次迭代节点电压幅值(uu)和相角(dd)的变化： K uu

15、1 uu2 uu3 dd1 dd2 dd3 1 0.9890 1.0500 1.0000 -0.0909 -0.0137 0 2 0.9850 1.0500 1.0000 -0.0923 -0.0138 0 3 0.9849 1.0500 1.0000 -0.0923 -0.0138 0迭代次数k-1=3节点电压相角（弧度）dd =节点 1 2 3 -0.0923 -0.0138 0节点电压幅值（标么值）uu =节点 1 2 30.9849 1.0500 1.0000节点功率矩阵（正号表示功率注入）Sm =节点 1 -1.7192 - 0.7346i 2 0.6661 + 1.2189i 3 1.0531 - 0.3265i支路功率（正号表示功率方向与支路方向一致）Smn =节点 1 2 3 1 0 -0.8111 - 0.6089i -0.9081