牛顿拉夫逊介绍(原理、计算方法、程序).doc

1、3 牛顿拉夫逊法概述3.1 牛顿-拉夫逊法基本原理电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。实际电力系统的潮流技术那主要采用牛顿-拉夫逊法。牛顿-拉夫逊法(简称牛顿法)在数

2、学上是求解非线性代数方程式的有效方法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程。即通常所称的逐次线性化过程。对于非线性代数方程组： 即 (3-1)在待求量x的某一个初始估计值附近，将上式展开成泰勒级数并略去二阶及以上的高阶项，得到如下的经线性化的方程组： (3-2)上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量 (3-3)将和相加，得到变量的第一次改进值。接着就从出发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求解的迭代格式为： (3-4) (3-5)上两式中：是函数对于变量x的一阶偏导数矩阵，即雅可比矩阵J；k为迭代次数。有上式可见，牛

3、顿法的核心便是反复形式并求解修正方程式。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快，具有平方收敛特性。牛顿潮流算法突出的优点是收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代45次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。牛顿法也具有良好的收敛可靠性，对于对以节点导纳矩阵为基础的高斯法呈病态的系统，牛顿法也能可靠收敛。牛顿法所需的内存量及每次迭代所需时间均较高斯法多。牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的节点上。对于正常运行的系统，各节点电压一般均在额定值附近，偏

4、移不会太大，并且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为平直电压)，如假定： 或 (3-6) 这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。解决这个问题的办法可以用高斯法迭代12次，以此迭代结果作为牛顿法的初值。也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。3.2 牛顿-拉夫逊法潮流求解过程以下讨论的是用直角坐标形式的牛顿拉夫逊法潮流的求解过程。当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量由于平衡节点的电压向量是给定的，因此待求

5、两共2(n-1)需要2(n-1)个方程式。事实上，除了平衡节点的功率方程式在迭代过程中没有约束作用以外，其余每个节点都可以列出两个方程式。对PQ节点来说，是给定的，因而可以写出 （3-7）对PV节点来说，给定量是，因此可以列出 （3-8）求解过程大致可以分为以下步骤：（1）形成节点导纳矩阵；（2）将各节点电压设初值U（3）将节点初值代入相关求式，求出修正方程式的常数项向量；（4）将节点电压初值代入求式，求出雅可比矩阵元素；（5）求解修正方程，求修正向量；（6）求取节点电压的新值；（7）检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步；（8）计算

6、支路功率分布，PV节点无功功率和平衡节点注入功率。以直角坐标系形式表示：迭代推算式采用直角坐标时,节点电压相量及复数导纳可表示为: (3-9)将以上二关系式代入上式中,展开并分开实部和虚部;假定系统中的第1,2,m号为PQ节点,第m+1,m+2,n-1为PV节点,根据节点性质的不同,得到如下迭代推算式: 对于PQ节点 (3-10)对于PV节点 (3-11)对于平衡节点平衡节点只设一个,电压为已知,不参见迭代,其电压为: (3-12)修正方程两组迭代式中包括2(n-1)个方程.选定电压初值及变量修正量符号之后代入,并将其按泰勒级数展开,略去二次方程及以后各项,得到修正方程如下: (3-13)其中

7、，； 。 (3-14)雅可比矩阵各元素的算式式(3-14)中, 雅可比矩阵中的各元素可通过对式(3-10)和(3-11)进行偏导而求得.当时, 雅可比矩阵中非对角元素为 (3-15)当时,雅可比矩阵中对角元素为: (3-16)由式(3-15)和(3-16)看出,雅可比矩阵的特点：矩阵中各元素是节点电压的函数,在迭代过程中,这些元素随着节点电压的变化而变化；导纳矩阵中的某些非对角元素为零时,雅可比矩阵中对应的元素也是为零.若,则必有；雅可比矩阵不是对称矩阵;雅可比矩阵各元素的表示如下： （3-17） （3-18） （3-19） （3-19） （3-20）3.3 牛顿拉夫逊法的程序框图输入原始数据

8、启动形成导纳矩阵 给定电压初值、置对于PQ节点，按式3-10计算、对于PU节点，按式3-11计算、是否用式3-15,3-16计算雅克比矩阵各元素解修正方程式，求用，修正节点电压以按系统的潮流分布计算平衡节点功率及线路功率输出以图3-1牛顿拉夫逊法的程序框图4 潮流仿真程序4.1 Matlab简介目前电子计算机已广泛应用于电力系统的分析计算，潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求：（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信性（4）适应性好，亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强（5）简单。 MATLA

9、B是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。另外，MATLAB提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信

10、号处理（SIGNAL PROCESSING）、控制系统（CONTROL SYSTEMS）、神经网络（NEURAL NETWORKS）、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。MATLAB设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便性和灵活性有着直接的关系。原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。4.2 矩阵的运算矩阵是MATLAB数据存储的基本单元，而矩阵的运算是MATLAB语言的核心，在MATLAB

11、语言系统中几乎一切运算均是以对矩阵的操作为基础的。矩阵的基本数学运算包括矩阵的四则运算、与常数的运算、逆运算、行列式运算、秩运算、特征值运算等基本函数运算，这里进行简单介绍。（1）四则运算矩阵的加、减、乘运算符分别为“+，*” ，用法与数字运算几乎相同，但计算时要满足其数学要求 在MATLAB中矩阵的除法有两种形式：左除“”和右除“/”。在传统的MATLAB算法中，右除是先计算矩阵的逆再相乘，而左除则不需要计算逆矩阵直接进行除运算。通常右除要快一点，但左除可避免被除矩阵的奇异性所带来的麻烦。在MATLAB6中两者的区别不太大。（2）与常数的运算 常数与矩阵的运算即是同该矩阵的每一元素进行运算。

12、但需注意进行数除时，常数通常只能做除数。（3）基本函数运算矩阵的函数运算是矩阵运算中最实用的部分，常用的主要有以下几个：det(a) 求矩阵a的行列式eig(a) 求矩阵a的特征值inv(a)或a (-1) 求矩阵a的逆矩阵rank(a) 求矩阵a的秩trace(a) 求矩阵a的迹（对角线元素之和）我们在进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数学运算，为有所区别，我们称之为数组运算。（4）基本数学运算数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别，数组的乘除法是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.”。前面讲过

13、常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之间的除法运算没有任何限制。另外，矩阵的数组运算中还有幂运算（运算符为 . ）、指数运算（exp）、对数运算（log）、和开方运算（sqrt）等。有了“对应元素”的规定，数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。（5）逻辑关系运算 逻辑运算是MATLAB中数组运算所特有的一种运算形式，也是几乎所有的高级语言普遍适用的一种运算。 4.3 牛顿拉夫逊法潮流计算程序本程序的功能是用 牛顿-拉夫逊法进行潮流计算，详见附录。结束语在电力系统的潮流计算算法中牛顿一拉夫逊法是得

14、到电力系统研究人员认可的算法之一，在本文中我们采用牛顿一拉夫逊法，主要是同时考虑到计算的准确和程序的运行速度。没有采用经常用的高斯迭代法，而是采用了传统的逆阵方法，是考虑到用MATLAB实现高斯迭代将会通过很多的循环迭代才能实现，而逆阵可以直接通过命令来求解，这必然可以大大节省时间。对于不能求逆的矩阵我们通过在电力系统中至少有一条支路上有接地支路来实现其求逆。对于P-Q分解法不能使用于有些R/X比较大的电力系统的缺点，可以通过在电力系统中并联补偿法或虚构节点来得以解决。通过实例计算分析，取得了比较满意的效果。基于MATLAB的电力系统潮流计算使计算机在计算、分析、研究复杂的电力系统潮流分布问题

15、上又前进了一步。不管采用什么算法，所有的潮流计算都是基于矩阵的迭代运算。而MATLAB语言正是以处理矩阵见长，实践证明，MATLAB语言在电力系统潮流计算仿真研究中的应用是可行的，而且由于其强大的矩阵处理功能，完全可以应用于电力系统的其它分析计算中；用MATLAB语言编程效率高，程序调试十分方便，可大大缩减软件开发周期，如果像控制界一样开发出电力系统自己的专用工具箱，将系统分析用的一些基本计算以函数的形式直接调用，那么更高层次的系统软件也可以很容易地实现。参考文献1 于永源，杨绮雯电力系统分析(第二版)北京:中国电力出版社,20042 刘从爱电力工程基础济南：山东科学技术出版社，19973 电

16、气工程师手册第二版编辑委员会电气工程师手册北京：机械工业出版社，20004 何仰赞，温增银电力系统分析.武汉：华中科技大学出版社，20025 吴际舜,侯志捡. 电力系统潮流计算的计算机方法. 上海：上海交通大学出版社, 2000 6 陈珩. 电力系统稳态分析（第三版）. 北京：中国电力出版社, 20077 孙丽华电力工程基础北京：机械工业出版社：20068 李维波. MATLAB在电气工程中的应用.北京 ：中国电力出版社，20079 陈恳直角坐标牛顿一拉夫逊法潮流计算新解法J电力系统及其自动化学报，1999，10 张铮MATLAB程序设计与实例应用M北京：中国铁道出版社，2O0311 张志涌,

17、徐彦琴. MATLAB教程-基于6.x版本.北京：北京航空航天大学出版社,200112 A. Abur, F. Magnago, Y. Lu, Educational Toolboxes for Power System Analysis, IEEE Computer Applications in Power,2000,13(4):31-3513 F. Milano, An Open Source Power System Analysis Toolbox, IEEE Trans. on Power Systems, 2005, 20(3) :1199 -1206.14 R.D. Zimme

18、rman, C.E. Murrillo-Snches, D. Gan, MatPower, A Matlab Power System Simulation Package. Version 3.1b2.Users Manual, Power Systems Enginnering Research Center, Cornell University, Ithaca, NY,Online Available: http:/www.pserc.cornell.edu/matpower15 Power System Analysis Toolbox, Documentation for PSAT

19、 version 2.0.0 82,February7,2007,OnlineAvailable:http:/www.power.uwaterloo.ca/fmilano/psat.htm16 Voltage Stability Toolbox (VST) - PC/Unix Version, Online Available:http:/www.pagesdrexel.edu/hgk22/VST/voltage_stability_toolbox.htm17 J. D. Weber, Implementation of a Newton-Based Optimal Power Flow in

20、to a Power System Simulation Environment, MS Thesis, University of Illinois, Urbana, 199718 Gh. CrQinR, Y-H. Song, Gh. GrigoraT, Optimal Operation and Planning of Power Systems, Casa de EditurR VENUS, IaTi, 2003致 谢值此本科论文完成之际，首先要感谢我的导师唐萃老师。唐老师从一开始的论文方向的选定，到最后的整篇文论的完成，都非常耐心的对我进行指导。给我提供了大量数据资料和建议，告诉我应该

21、注意的细节问题，细心的给我指出错误，修改论文。唐老师诲人不倦的工作作风，一丝不苟的工作态度，严肃认真的治学风格给我留下深刻的影响，值得我永远学习。在此，谨向导师唐老师致以崇高的敬意和衷心的感谢！学部领导孙亲锡主任、电气工程及其自动化系主任刘沛、辅导员张甸老师等为我提供了良好的研究条件，谨向各位表示诚挚的敬意和谢忱。感谢我的同学郑鑫这学期来和我一起完成毕业论文的过程中给予我的帮助。附录Matlab仿真程序n=input(请输入节点数：n=);n1=input(请输入支路数：n1=);isb=input(请输入平衡母线节点号：isb=);pr=input(请输入误差精度：pr=);B1=input

22、(请输入由支路参数形成的矩阵：B1=);B2=input(请输入各节点参数形成的矩阵：B2=);X=input(请输入由节点参数形成的矩阵：X=);Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=seros(1,n);O=zeros(1,n);S1=zeros(n1);for i=1:n if X(i,2)=0; p=X(i,1); Y(p,p)=1./X(i,2); end end for i=1:n1 if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Y(p,q)=Y(p,q)-1.

23、/(B1(i,3)*B1(i,5); Y(p,q)=Y(p,q); Y(p,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)2)+B1(i,4)./2; Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; end %求导纳矩阵 G=real(Y);B=imag(Y); for i=1:n e(i)=real(B2(i,3); f(i)=imag(B2(i,3); V(i)=B2(i,4); endfor i=1:n S(i)=B2(i,1)-B2(i,2); B(i,i)=B(i,i)+B2(i,5); endP=rea(S);Q=imag(S);ICT1=0;

24、IT2=1;NO=2*n;N=NO+1;a=0;while IT2=0 IT2=0;a=a+1; for i=1:n; if i=isb C(i)=0; D(i)=0; for j1=1:n C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); end P1=C(i)*e(i)+f(i)*D(i); Q1=f(i)*C(i)-D(i)*e(i); %求P,Q V2=e(i)2+f(i)2; if B2(i,6)=3 DP=P(i)-P1; DQ=Q(i)-Q1; for j1=1:n if j

25、1=isb&j1=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X3=X2; X4=-X1; p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; end endelse DP=P(i)-P1; DV=V(i)2-V2; for j1=1:n if j1=isb&j1=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i)

26、; X5=0; X6=0; p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2; elseif j1=i&j1=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X5=-2*e(i); X6=-2*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,

27、q)=X2; end end end end end %求雅可比矩阵for k=3:N0 k1=k+1;N1=N; for k2=k1:N1 J(k,k2)=J(k,k2)./J(k,k);endJ(k,k)=1;if k=3 k4=k-1; for k3=3:k4 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end end for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end endendfor k=3:2

28、:N0-1 L=(k+1)./2; e(L)=e(L)-J(k,N); k1=k+1; f(L)=f(L)-J(k1,N);end for k=3:N0 DET=abs(J(k,N); if DET=pr IT2=IT2+1 endend ICT2(a)=IT2 ICT1=ICT1+1;for k=1:n dy(k)=sqrt(e(k)2+f(k)2);endfor i=1:n Dy(k)=sqrt(e(k)2+f(k)2); endfor i=1:n Dy(ICT1,i)=dy(i); end end %用高斯消去法解“w=-J*V”disp(迭代次数);disp(ICT1);disp(没有

29、达到精度要求的个数);disp(ICT2);for k=1:n V(k)=sqrt(e(k)2+f(k)2); O(k)=atan(f(k)./e(k)*180./pi;endE=e+f*j;disp(各节点的实际电压标么值E为（节点号从小到大的排列）：);disp(E);disp(各节点的电压大小V为（节点号从小到大的排列）：);disp(V);disp(各节点的电压相角O为（节点号从小到大的排列）：);disp(O);for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); end S(p)=E(p)*C(p);end dis

30、p(各节点的功率S为（节点号从小到大排列）：）； disp(S); disp(各条支路的首端功率Si为（顺序同您输入B1时一样）：）； for i=1:n1 if B1（i,6）=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); endSi(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5)-conj(E(q)*conj(1./(B1(i,3)*B1(i,5); disp(Si(p.q);enddisp(各条支路的末端功率Sj为（顺序同您的输入B1时一样）：）；for i=1:n1 i

31、f B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(xonj(E(q)./B1(i,5)-conj(E(p)*xonj(1./(B1(i,3)*B1(i,5);disp(Sj(q,p);end disp(各条支路的功率损耗DS为（顺序同您输入B1时一样）：）；for i=1:n1 if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end DS(i)=Si(p,q)+Sj(q,p); disp(DS(i); end for i=1:ICT1 Cs(i)=i; enddisp(以下是每次迭代后各节点的电压值（如图所示）；plot(Cs,Dy),xlabel(迭代次数），ylabel(电压),title(电压迭代次数曲线);