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

1牛顿-拉夫逊迭代法电力网潮流计算方法与程序编写佘名寰牛顿-拉夫逊迭代法在电力网潮流计算中因其收敛性较好获得广泛运用，该算法的难点是需反复计算功率方程中雅可比矩阵各个元素表达式。本文简叙了牛顿-拉夫逊迭代法潮流计算的基本公式，通过例题介绍了用牛顿-拉夫逊法计算电力网潮流电压的MATLAB程序。程序采用MATLAB语言的符号矩阵简化了雅可比矩阵系数的计算。本文可供电力系统电气技术人员和大专院校电力类专业师生参考。2．牛顿-拉夫逊迭代法潮流计算2.1极坐标型式节点功率方程式由电源注入或从负载流出节点的电流统称节点电流，由节点电压和节点电流所求得的节点注入功率为：Sti=V̇iÎiİi̇=∑YijV̇j(I,j=1,2,….n)Pti=Vi∑Vjn1(Gijcosθij+Bijsinθij)Qti=Vi∑Vjn1(Gijsinθij−Bijcosθij)(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=QtiΔPi=Pgi-PlI-Pti=0ΔQi=Qgi-Qli-Qti=0（2-02）Pgi,Qgi——节点i发电机输入有功和无功功率Pli,Qli——节点i负荷有功和无功功率ΔPi，ΔQi--节点i不平衡功率不平衡功率的微分d(ΔPi),d(ΔQi)d(∆pi)=−(∂pti∂v1∆v1+∂pti∂v2∆v2…+∂pti∂ϑ1∆ϑ1+∂pti∂ϑ2∆ϑ2….)（i=1,2,…,n）d(∆qi)=−(∂qti∂v1∆v1+∂q∂v2∆v2…+∂qti∂ϑ1∆ϑ1+∂qti∂ϑ2∆ϑ2….)（2-03）对于n个节点系统可得如下矩阵形式修正方程式；2[∆P1∆P2:∆Pn∆Q1∆Q2:∆Qn]=[∂∆P1∂ϑ1∂∆P2∂ϑ1:∂∆Pn∂ϑ1∂∆Q1∂ϑ1∂∆Q2∂ϑ1:∂∆Qn∂ϑ1∂∆P1∂ϑ2∂∆P2∂ϑ2:∂∆Pn∂ϑ2∂∆Q1∂ϑ2∂∆Q2∂ϑ2:∂∆Qn∂ϑ2⋯⋯⋯⋯⋯⋯⋯⋯∂∆P1∂ϑn∂∆P2∂ϑn:∂∆Pn∂ϑn∂∆Q1∂ϑn∂∆Q2∂ϑn:∂∆Qn∂ϑn∂∆P1∂v1∂∆P2∂v1:∂∆Pn∂v1∂∆Q1∂v1∂∆Q2∂v1:∂∆Qn∂v1∂∆P1∂v2∂∆P2∂v2:∂∆Pn∂v2∂∆Q1∂v2∂∆Q2∂v2:∂∆Qn∂v2⋯⋯⋯⋯⋯⋯⋯⋯∂∆P1∂vn∂∆P2∂vn:∂∆Pn∂vn∂∆Q1∂vn∂∆Q2∂vn:∂∆Qn∂vn][Δϑ1∆ϑ2:∆ϑn∆v1∆v2:∆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)⑹用新的节点电压重复以上五个步骤计算直到计算精度合符要求。牛顿-拉夫逊迭代法计算的难点是求雅可比矩阵系数，本文采用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是平衡节点，其电压保持定值，大小为V3=1.0。3图2-1[例2.3.1]网络接线图2.3.1.1程序参变量说明Y=G+j*B节点导纳矩阵，G为电导，B为电纳；P节点有功功率给定值，正为发电机注入，负为负荷功率；q节点无功功率给定值，正为发电机注入，负为负荷功率；u节点电压幅值的符号矩阵；delt节点电压相角的符号矩阵；pt节点有功功率符号表达式；qt节点无功功率符号表达式；pp节点有功功率不平衡值符号表达式；qq节点无功功率不平衡值符号表达式；uu节点电压幅值数值矩阵；dd节点电压相角数值矩阵；PP节点不平衡功率的数值矩阵N1网络独立节点总数；N2网络PV节点总数；Sm节点功率矩阵；Smn支路功率矩阵；J1,J2,J节点不平衡功率雅可比符号矩阵.JJ节点不平衡功率雅可比数值矩阵2.3.1.2源程序源程序如下：%*****NUE78.mexample2-3-1fig2-1*******************************%ThefollowingProgramforloadflowcalculationisbasedonMATLAB7.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;Z13Z12Z23YC3YC1YC1V3=1.0312YC2P1=1.7192Q1=0.7346P2=0.6661V2=1.05Z21=j0.1Z23=j0.1Z13=j0.1YC1=J0.02YC2=J0.02YC3=J0.024N1=2;%theN1istheamountofthePQandPVbusN2=1;%theN2istheamountofthePVbusform=1:N1forn=1:N1+1pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));endpp(m)=p(m)-sum(pt);endform=1:N1-N2forn=1:N1+1qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));endqq(m)=q(m)-sum(qt);endJ1=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];whileprecision0.00001u1=uu(1);u2=uu(2);u3=uu(3);d1=dd(1);d2=dd(2);d3=dd(3);form=1:N1PP(m)=eval(pp(m));endform=1:N1-N2PP(N1+m)=eval(qq(m));endJJ=eval(J);du=-inv(JJ)*PP';precision=max(abs(du));forn=1:N1dd(n)=dd(n)+du(n);endforn=1:N1-N2uu(n)=uu(n)+du(N1+n);endk=k+1;endk-1,dd,uu%***********************************************5%thefollowingprogramisusedtocalculatetheSmandSmnforn=1:N1+1U(n)=uu(n)*(cos(dd(n))+j*sin(dd(n)));endUm=conj(U');I=Y*Um;Sm=diag(Um)*conj(I)form=1:N1+1forn=1:N1+1Smn(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.98Y12=Y23=Y13=0+j10.0⑵节点电压符号矩阵：一共三个节点，独立节点数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)⑸qt(n),qq(m)相关语句是建立ΔQi符号表达式，只针对PQ节点；ΔQ1=qq(1)=-3673/5000-999/50*u1^2+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),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)的变化：6Kuu1uu2uu3dd1dd2dd310.98901.05001.0000-0.0909-0.0137020.98501.05001.0000-0.0923-0.0138030.98491.05001.0000-0.0923-0.01380迭代次数k-1=3节点电压相角（弧度）dd=节点123-0.0923-0.01380节点电压幅值（标么值）uu=节点1230.98491.05001.0000节点功率矩阵（正号表示功率注入）Sm=节点1-1.7192-0.7346i20.6661+1.2189i31.0531-0.3265i支路功率（正号表示功率方向与支路方向一致）Smn=节点12310-0.8111-0.6089i-0.9081-0.1063i20.8111+0.7149i0-0.1450+0.5260i30.9081+0.1925i0.1450-0.4990i0参考文献：①电气工程专业毕业设计指南电力系统分册，陈跃主编，中国水利水电出版社2003年②电力系统计算，西安交通大学清华大学等合编，水利电力出版社，1978年③MATLAB电力系统设计与分析吴天明等编箸国防工业出版社2007年④MATLAB基础与应用简明教程张平等编著北京航空航天大学出版社2001年