1、牛顿拉夫逊潮流计算 张海强 学号： 2017200385 西南交通大学电气工程学院 Email: 摘要 : 牛顿拉夫逊法收敛性好，迭代次数少 ,在电力系统等领域得到较好的应用 。 本文主要介绍了 牛顿拉夫逊法对 电力系统 进行 潮流计算 的基本步骤 ， 并利用 MATLAB 仿真软件对简单电路的潮流进行了运算。 关键词： 牛顿 拉夫逊 算法 ； 潮流 计算 ； MATLAB. 1 引言 电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。在电力系统规划的设计和现有电力系
2、统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性。可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。所以潮流计算是研究电力系统的一种很重要和基础的计算。 电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运 行方式，后者则用于正在运行系统的经常监视及实时控制。 利用电子数字计算机进行电力系统潮流计算从 50年代中期就已经开始。在这 20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为下面几点： （ 1）计算方法的可靠性或收敛性； （
3、2）对计算机内存量的要求； （ 3）计算速度； （ 4）计算的方便性和灵活性。 电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解法都离不开迭代。因此，对潮流计算方法，首先要求它能可靠地收敛 ，并给出正确答案。由于电力系统结构及参数的一些特点，并且随着电力系统不断扩大，潮流计算的方程式阶数也越来越高，对这样的方程式并不是任何数学方法都能保证给出正确答案的。这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。 牛顿拉夫逊算法是数学中求解非线性方程的典型方法，能快速求出其他方法求不出或者难以求出的解 ,牛顿拉夫逊潮流计算方法收敛性好，迭代次数少，在潮流计算方法中得到广泛
4、的应用。所以本文将从牛顿法的基本原理、潮流计算方程、实例仿真三个部分对牛顿法潮流计算进行介绍。 2 牛顿法基本原理 牛顿 -拉夫逊法是解非线性方程式的有效方法。牛顿拉夫逊法潮流计算是目前最为广泛、效果最好的一种潮流计算方法。这种把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，即逐次线性化过程，这就是牛顿法的核心。我们以如下非线性方程式的求解过程为例来说明： （ 1） 0)( xf设 为该方程式的初值。而真正解 x在它的近旁： （ 2） 式中： 为初始值 的修正量。如果求得 ，则由式（ 3-2）就可以得到真正解 x。 为此将式 （ 3） 按泰勒级数展开 （ 4） 当我们选择的初始值
5、比较好，即 很小时，式（ 3-4）中包含的 和更高阶次项可以略去不计。因此，式（ 3-4）可以简化为 （ 5） 这是对于变量 的形式方程式，用它可以求出修正量 。 由于式（ 5）是式（ 4）的简化结果，所以由式（ 5）解出 后，还不能得到方程式（ 1）的真正解。实际 上，用 对 修正后得到的 ： （ 6） 现在如果再以作为初值 ，解式（ 3-5就能得到更趋近真正解的 ： （ 7） 这样反复下去，就构成了不断求解非线性方程式的逐次线性化过程。第 t次迭代时的参数方程为 （ 8） （ 9） 上式左端可以看成是近似解 引起的误差，当 时，就满足了原方程式（ 1），因而 就成为该方程的解。式中 是函数
6、 在 点的一次导数，也就是曲线在 点的斜率，如图 1所示，修正量 则是由 点的切线与横轴的交点来确定，由图 1可以直观的看出牛顿法的求解过程。 )0(x )0()0( xxx )0(x )0( )0(x0)( )0()0( xxf0! )()()1(!2 )()()()()( )0()0()()(2)0()0()0()0()0()0()0( nxxfxxfxxfxfxxf nnn)0(x 2)0( )( x0)()( )0()0()0( xxfxf)0(x )0(x)0()0(x )1(x)0()0()1( xxx )1(x 2(x)1()1()2( xxx 0)()( )()()( ttt
7、xxfxf)()()( )( ttt xxfxf )tx 0)( )( txf)(tx )( )(txf 0)( xf )(t)(t )(tt图 1 牛顿 -拉夫逊法几何解释 现在把牛顿法推广到多变量非线性方程组的情况。设有变量 的非线性联立方程组 ： （ 10） 给定各变量初值 ，假设 为其修正量，并使其满足 （ 11） 对以上 n个方程式分别按泰勒级数展开，当忽略 所组成的二次项和高次项时，可以得到 0)( )1( txf)( xfy )( )( txf)1( tx )(tx)1( tx)( txXYXnxxx 21,0),(0),(0),(21212211nnnnxxxfxxxfxxxf
8、)0()0(2)0(1 , nxxx )0()0(2)0(1 , nxxx 0),(0),(0),()0()0()0(2)0(2)0(1)0(1)0()0()0(2)0(2)0(1)0(12)0()0()0(2)0(2)0(1)0(11nnnnnnnxxxxxxfxxxxxxfxxxxxxf)0()0(2)0(1 , nxxx （ 12） 式 中 ： 为函数 对 自 变 量 的 偏 导 数 在 点（ ）处的值。把上式写成矩阵形式： （ 13） 这是变量 的线性方程组，称为牛顿法的修正方程，通过它可以解出 ，并可以进一步求得 （ 14） 式中 向真正解逼近了一步，如果再以它们作为初值重复解式（
9、13）修正方程式，等到更接近真解的 ，如此迭代下去，并按式（ 14）进行修正，直到满足收敛要求为止并停止迭代计算，这就构成了牛顿法的迭代过程。 一般第 t次迭代式的修正方程为 （ 15） 上式可以简写为 0),(0),(0),()0(0)0(202)0(101)0()0(2)0(1)0(02)0(2022)0(1012)0()0(2)0(12)0(01)0(2021)0(1011)0()0(2)0(11nnnnnnnnnnnnnxxfxxfxxfxxxfxxfxxfxxfxxxfxxfxxfxxfxxxf0iixf),( 21 ni xxxf j )0()0(2)0(1 , nxxx )0()
10、0(2)0(1002010202201201021011)0()0(2)0(1)0()0(2)0(12)0()0(2)0(11),(),(),(nnnnnnnnnnnxxxxfxfxfxfxfxfxfxfxfxxxfxxxfxxxf)0()0(2)0(1 , nxxx )0()0(2)0(1 , nxxx )0()0()1()0(2)0(2)1(2)0(1)0(1)1(1nnn xxxxxxxxx)1()1(2)1(1 , nxxx )2()2(2)2(1 , nxxx )()(2)(1212221212111)()(2)(1)()(2)(12)()(2)(11),(),(),(tntttnn
11、tntntntttntttnttntntttnttxxxxfxfxfxfxfxfxfxfxfxxxfxxxfxxxf(16) 其中 ， 其中的 为第 t次迭代时的雅克比矩阵； 同理可以得到第 t次迭代时的修正量： （ 17） 同样，也可以写出类似（ 14）的算式 （ 18） 这样反复交替的解式（ 16）及式（ 18）就可以使 逐步趋近方程式的真正解。当满足人为收敛条件时，即 或 （ 19） 迭代结束，式中 为预先给定的小正数。 3 牛顿法潮流计算方程 3.1 节点功率方程 电力系统的负荷习惯用功率表示，对于有 n个节点的电力系统，系统中各节点注入电流与注入功率以标幺值表示的关系为 i=1， 2
12、， ， n （ 20） 式中 表示其共轭复数。将此关系式代入节点电压方程的通式，可得到以节点注入功率表示的节点电压方程 : (21) 上述的方程式，通常称为功率方程。根据方程中的节点电压向量表示的不同，可以得到不)()()( )( ttt XJXF ),(),(),()()()(2)(1)()(2)(12)()(2)(11)(tnttntntttntttxxxfxxxfxxxfXFtnntntntntttntttxfxfxfxfxfxfxfxfxfJ212221212111)()(tJ)()(2)(1)(tntttxxxX)()()1( ttt XXX )1(tX 1)()(2)(1 ,(m
13、ax tntti xxxf 2)(max tix21,*)( iiii UjQPUSI nj jijiii UYUjQP 1)( 同形式的功率方程。 若节点电压向量以直角坐标表示，即以复数平面上实轴与虚轴上的投影表示可写成 （ 22） 其共轭值为 （ 23） 导纳表示为 （ 24） 把这两关系式代回式（ 21）的功率方程中，展开后再将功率方程的实部和虚部分别写成有功、无功功率分离的节点方功率方程： （ 25） 式中 ： i=1， 2， ， n为各节点的编号 。 若节点电压以极坐标表示，则 或写成 （ 26） 将其同导纳的复数表达式一起代入式（ 21）的功率方程，进整理可以得到 （ 27） 式中
14、： i与 j节点电压的相角差 。 由式（ 25）和（ 27）给出的功率方程表示方法避免了复数运算，因此，在潮流计算中普遍采用 。 3.2 修正方程 采用牛顿法计算潮流时，需要对功率方程进行修改。下面将根据在不同坐标内的修改进行讨论： （ 1）在直角坐标系内时，由 PQ 节点功率方程（ 25）可知：节点 i 的注入功率是各点电压的函数，设节点的电压已知，代入式（ 25），可以求出节点 i的有功及无功功率 ，它们与给定的 PQ 节点的注入功率 的差值应满足以下方程 iii jfeU iii jfeU ijijij jBGY njnjjijjijiiijjijiinjnjjijjijiiijjiji
15、ieBfGefBeGfQeBfGffBeGeP1 11 1)()()()(ijii eUU iiiii jUUU si nc o s )c o ss i n()s i nc o s(11njijijijijjiinjijijijijjiiBGUUQBGUUPjiij ii QP,isisQP,（ 28） 对于 PV 节点，已知节点的注入有功功率及节点电压大小，记作 ，其节点的有功功率应满方程： （ 29） 对于平衡节点，因为其电压给定，故不需要迭代求解。 通过以上分析可见，式（ 28）和式（ 29）共 2（ n-1）个方程，待求量共 2（ n-1）个。将上述 2（ n-1）个方程按泰勒级数展开
16、，并略去修正量的高次方项后得到修正方程如下： （ 30） 其中雅克比矩阵的各元素可以对式（ 28）和式（ 29）求偏导数获得。 对于非对角元素（ ）有 0)()(0)()(1 11 1njnjjijjijiiijjijiisiisinjnjjijjijiiijjijiisiisieBfGefBeGfQQQQeBfGffBeGePPPPisisUP,)1, . . . . . . . ,2,1(0)(0)()(22221 1 nmmifeUUeBfGffBeGePPPPiiisinjnjjijjijiiijjijiisiisi112211 , nn fefefe UJW Tnnmmmm UPUP
17、QPQPW 2 112 1111 Tnnmmmm fefefefeU 111111 12112112112121121121121111111111111111211211211212121121121111111111111111111111111111111111111111111111111111111nnnnmnmnmnmnnnnnnnmnmnmnmnnnnmnmmmmmmmmmmmnmnmmmmmmmmmmmnmnmmmmmmmmmmmnmnmmmmmmmmmmmnnmmmmnnmmmmfUeUfUeUfUeUfUeUfPePfPePfPePfPePfUeUfUeUfUeUfUeUf
18、PePfPePfPePfPePfQeQfQeQfQeQfQeQfPePfPePfPePfPePfQeQfQeQfQeQfQeQfPePfPePfPePfPePJji（ 31） 对于对角元素（ 有 （ 32） 由上述表达式可以看到，雅克比矩阵具有以下特点： 1） 各元素是各节点电压的函数，迭代过程中每迭代一次各节点电压都要变化，因而各元素每次也变化； 2） 雅克比矩阵不具有对称性； 3） 互导纳 ，与之对应的非对角元素亦为零，此外因非对角元素，故雅克比矩阵是稀疏矩。 当在极坐标系内时，由功率方程（ 27）可知节点 i的注入功率是各节点电压幅值和相角的函数。代入式（ 27）可以求出节点 i的有功功
19、率和无功功率，它们与给定的 PQ节点的注入功率 的差值满足下面方程： （ 33） 式中： i与 j节点电压的相角差 。 在有 n个节点的系统中，假定第 号节点为 PQ节点，第 m+1n-1号节点为 PV节点，第 n号节点为平衡节点。 和 是给定的， PV节点的电压幅值 也是给定的，因此，只剩下 n-1个节点的电压相角 和 m个节点的电压幅值 是未0)(22jijiiijiijjijiiijiijjijifUeUfGeBeQfPfBeGfQeP)jiiiiiiinjiiiiiijijjijiinjiiiiiijijjijiinjiiiiiijijjijiinjiiiiiijijjijiiffUe
20、eUfBeGfBeGfQfGeBeBfGeQeBfGeBfGfPfBeGfBeGeP22)()()()(2211110ijY 022 j ij i fUeU isisQP,0)coss i n(0)s i ncos(11njijijijijjiiinjijijijijjiiiBGUUQQBGUUPPjiij m1nVn11 nm VV 121 , n mVVV , 21知量。由式（ 33）可知一共包含了 n-1+m 方程式，正好同未知量的数目相等，而直角坐标形式的方程少了 n-1-m个。 由方程（ 33）可以写出修正方程 (34) 式中 （ 35） 其中： H 是 阶方阵，其元素为 ； N 是
21、 阶矩阵，其元素为 ； K是 阶矩阵，其元素为 ； L是 阶矩阵，其元素为 。 对式（ 33）求偏导数，可得雅克比矩阵元素的表达式如下： 非对角元素（ ） （ 36） 对角元素（ ） （ 37） 4 基于 MATLAB的实例仿真 4.1 潮流计算原始资料参数 图 2是一个简单的三相系统的单线图，节点 1、 3连接发电机，节点 1的电压幅值为 1.05 VVLK NQPD 12mDmnmnVVVVVVVVQQQQPPPP2121121211212;)1()1( nn jiij P mn )1(iijij VPVN )1( nmjiij QK mmjijij VQVL ji)co ss i n()
22、s i nco s()s i nco s()co ss i n(ijijijijjiijijijijijjiijijijijijjiijijijijijjiijBGVVLBGVVKBGVVNBGVVHjiiiiiiiiiiiiiiiiiiiiiiijQBVLPGVKPGNQBVH2222（标幺值），节点 3的电压幅值为 1.04（标幺值），发电机有功功率为 200MW。节点 2连接一负荷，负荷吸收有功功率为 400MW，无功功率为 250Mvar。线路阻抗如图所示，忽略线路充电电纳，基准功率为 100MVA。 图 2 简单三相系统 负荷以及发电机阻抗转换成标幺值为： 2 ( 4 0 0 2 5
23、 0 ) 4 . 0 2 . 5 p u100s c h jSj 3 200 2.0 pu100schP 节点 1作为参考节点 (平衡节点 )，节点 2、 3的电压初始估计值为： (0)2 1.0 0.0Vj、 (0)3 1.04 0.0Vj 4.2 牛顿法计算潮流步骤 根据以上对牛顿法进行电力系统潮流计算原理的介绍可得牛顿法计算潮流步骤如下： 给这各节点电压初始值 )0()0( ,fe ； 将以上电压初始值代入公式，求修正方程的常数项向量 )0(2)0()0( )(, VQP ； 将电压初始值在带入公式，求出修正方程中系数矩阵的各元素。 解修正方程式 )0()0( , fe ； 修正各节点电
24、压 )0()0()1( eee ， )0()0()1( fff ； 将 )1(e ， )1(f 在带入方程式，求出 )1(2)1()1( )(, VQP ； 检验是否收敛，即 )()( ,m ax kiki QP 如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回（ 2）进行下次迭代计算，直到收敛为止。 牛顿法潮流计算计算程序图如图 3所示。 图 3 潮流计算程序流程图 根据程序流程图可编制潮流计算程序，本文主要以以节点导纳矩阵计算程序为例进行介绍，具体代码如下： functionY=ybus1(zdata) nl=zdata(:,1);nr=zdat
25、a(:,2);R=zdata(:,3);X=zdata(:,4); nbr=length(zdata(:,1);nbus=max(max(nl),max(nr); Z=R+j*X; %支路阻抗 y=ones(nbr,1)./Z; %支路导纳 Y=zeros(nbus,nbus); %将 Y初始化为 0 For k=1:nbr; %费对角元素的数值 if nl(k)0 Y(nr(k),nl(k)=Y(nl(k),nr(k); end 输入原始数据 形成节点导纳阵 给定电压初值 e（ 0） 、 f（ 0） k=0 根据公式计算 用公式计算雅克比矩阵各元素 高斯法解修正方程，求 修正节点电压 是 计
26、算平衡节点功率及全部线路功率 输出 k=k+1 否 end for n=1:nbus %对角元素数值 for k=1:nbr if nl(k)=n|nr(k)=n Y(n,n)=Y(n,n)+y(k); else ,end end end 函数 Y=ybus1（ zdata ） 用来形成节点导纳矩阵， zdata是线路数据输入变量，包括四项，前两项是节点编号，后两项是线路电阻和电抗，均以标幺值表示，函数返回节点导纳矩阵。对于电力系统程序来说，节点导纳矩阵算法是非常简单和基础的，因此这里详细讲解，要求读者了解并掌握这种算法。编程时先将线路阻抗转换为导纳，然后将矩阵 Y初始化为零，执行第一次循环时
27、，查找线路数据，输入非对角元素，最后，在嵌套循环中寻找与节点有关的元素，对角元素随之形成。 则导纳矩阵为： 2 0 5 0 1 0 2 0 1 0 3 01 0 2 0 2 6 5 2 1 6 3 21 0 3 0 1 6 3 2 2 6 6 2busj j jY j j jj j j 潮流计算具体程序如本文附件所示，求得结果： 3 1.4617 puQ ， 1 2.1842 puP ， 1 1.4085 puQ 。 参考文献 1 王晓茹 , 高仕斌 . 电力系统分析 M. 高等教育出版社 , 2011. 2 张智星 . MATLAB程序设计与应用 M. 清华大学出版社 , 2002. 3 阮
28、沈勇 . MATLAB程序设计 M. 电子工业出版社 , 2004. 4 牛辉 , 郭志忠 . 电流注入模型的电力系统潮流计算 J. 电网技术 , 1998, 22(11):39-41. 5 杨镜非 , 李维 . 电力系统潮流计算组件的实现 J. 电网技术 , 2000, 24(4):57-59. 6 薛振宇 , 房大中 . 基于双向迭代的交直流互联电力系统潮流计算 J. 电力系统自动化 , 2013, 37(5):61-67. 7 罗日成 , 李卫国 . 电力系统潮流计算的符号分析方法 J. 电网技术 , 2005, 29(10):25-29. 附件：设计程序和运行结果 1 设计程序 V =
29、 1.05; 1.0; 1.04; d = 0; 0; 0; Ps=-4; 2.0; Qs= -2.5; YB = 20-j*50 -10+j*20 -10+j*30 -10+j*20 26-j*52 -16+j*32 -10+j*30 -16+j*32 26-j*62; Y= abs(YB); t = angle(YB); iter=0; pwracur = 0.00025; % Power accuracy DC = 10; % Set the maximum power residual to a high value while max(abs(DC) pwracur iter = i
30、ter +1 P=V(2)*V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1)+V(2)2*Y(2,2)*cos(t(2,2)+ . V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3); V(3)*V(1)*Y(3,1)*cos(t(3,1)-d(3)+d(1)+V(3)2*Y(3,3)*cos(t(3,3)+ . V(3)*V(2)*Y(3,2)*cos(t(3,2)-d(3)+d(2); Q= -V(2)*V(1)*Y(2,1)*sin(t(2,1)-d(2)+d(1)-V(2)2*Y(2,2)*sin(t(2,2)- . V(2)*V(3)*Y(2,
31、3)*sin(t(2,3)-d(2)+d(3); J(1,1)=V(2)*V(1)*Y(2,1)*sin(t(2,1)-d(2)+d(1)+. V(2)*V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3); J(1,2)=-V(2)*V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3); J(1,3)=V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1)+2*V(2)*Y(2,2)*cos(t(2,2)+. V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3); J(2,1)=-V(3)*V(2)*Y(3,2)*sin(t(3,2)-d(3)+
32、d(2); J(2,2)=V(3)*V(1)*Y(3,1)*sin(t(3,1)-d(3)+d(1)+. V(3)*V(2)*Y(3,2)*sin(t(3,2)-d(3)+d(2); J(2,3)=V(3)*Y(2,3)*cos(t(3,2)-d(3)+d(2); J(3,1)=V(2)*V(1)*Y(2,1)*cos(t(2,1)-d(2)+d(1)+. V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3); J(3,2)=-V(2)*V(3)*Y(2,3)*cos(t(2,3)-d(2)+d(3); J(3,3)=-V(1)*Y(2,1)*sin(t(2,1)-d(2
33、)+d(1)-2*V(2)*Y(2,2)*sin(t(2,2)-. V(3)*Y(2,3)*sin(t(2,3)-d(2)+d(3); DP = Ps - P; DQ = Qs - Q; DC = DP; DQ J DX = JDC d(2) =d(2)+DX(1); d(3)=d(3) +DX(2); V(2)= V(2)+DX(3); V, d, delta =180/pi*d; end P1= V(1)2*Y(1,1)*cos(t(1,1)+V(1)*V(2)*Y(1,2)*cos(t(1,2)-d(1)+d(2)+. V(1)*V(3)*Y(1,3)*cos(t(1,3)-d(1)+d
34、(3) Q1=-V(1)2*Y(1,1)*sin(t(1,1)-V(1)*V(2)*Y(1,2)*sin(t(1,2)-d(1)+d(2)-. V(1)*V(3)*Y(1,3)*sin(t(1,3)-d(1)+d(3) Q3=-V(3)*V(1)*Y(3,1)*sin(t(3,1)-d(3)+d(1)-V(3)*V(2)*Y(3,2)*. sin(t(3,2)-d(3)+d(2)-V(3)2*Y(3,3)*sin(t(3,3) 2 运行结果 Chp6ex10 iter = 1 DC = -2.8600 1.4384 -0.2200 J = 54.2800 -33.2800 24.8600 -3
35、3.2800 66.0400 -16.6400 -27.1400 16.6400 49.7200 DX = -0.0453 -0.0077 -0.0265 V = 1.0500 0.9735 1.0400 d = 0 -0.0453 -0.0077 iter = 2 DC = -0.0992 0.0217 -0.0509 J = 51.7247 -31.7656 21.3026 -32.9816 65.6564 -15.3791 -28.5386 17.4028 48.1036 DX = -0.0018 -0.0010 -0.0018 V = 1.0500 0.9717 1.0400 d = 0 -0.0471 -0.0087 iter = 3 DC = 1.0e-03 * -0.2166 0.0382 -0.1430 J = 51.5967 -31.6939 21.1474 -32.9339 65.5976 -15.3516 -28.5482 17.3969 47.9549 DX = 1.0e-05 * -0.3856 -0.2386 -0.4412 V = 1.0500 0.9717 1.0400 d = 0 -0.0471 -0.0087 P1 = 2.1842 Q1 = 1.4085 Q3 = 1.4618
