潮流计算的牛顿—拉夫逊法

潮流计算的牛顿—拉夫逊法王亮2012109427电力系统潮流计算是研究电力系统稳态运行情况的一种基本电气计算。它的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。本文介绍了牛顿—拉夫逊法的基本原理以及在Matlab下用此方法做电力系统的潮流计算的初步知识。1、牛顿—拉夫逊法概要首先对一般的牛顿—拉夫逊法作一简单的说明。已知一个变量X函数为：0)(Xf（0）到此方程时，由适当的近似值)0(X出发，根据：,......)2,1()()()()()()1(nXfXfXXnnnn（1）反复进行计算，当)(nX满足适当的收敛条件就是上面方程的根。这样的方法就是牛顿—拉夫逊法。这一方法还可以做下面的解释，设第n次迭代得到的解与真值之差，即)(nX的误差为时，则：0)()(nXf（2）把)()(nXf在)(nX附近对用泰勒级数展开0......)(!2)()()()(2)()()(nnnnXfXfXfXf（3）上式省略去2以后部分0)()()()(nnXfXf（4）)(nX的误差可以近似由上式计算出来。)()()()(nnXfXf（5）比较两式，可以看出牛顿—拉夫逊法的休整量和)(nX的误差的一次项相等。用同样的方法考虑，给出n个变量的n个方程：0),,,(0),,,(0),,,(21212211nnnnXXXfXXXfXXXf（6）假定变量的初始值给定为(0)(0)(0)12nX,X,...X，并设(0)(0)(0)12nX,X,...X分别为各变量的修正值，则有(0)(0)(0)(0)(0)(0)1122nn(0)(0)(0)(0)(0)(0)1122nn(0)(0)(0)(0)(0)(0)112221nnnf(,,,)0f(,,XXXXXXXXXXXXXXX,)0...f(,,,)XXX0（7）将上式在初始值附近展开成泰勒级数，并略去修正量的二次及以上的高次项111112n000(0)(0)(0)(0)(0)(0)12n12n(0)(0)(0)(0)(0)(0)12n12n(0)(0)(0)(0)(222212n000nnn10)1210n022ffff(,,,)()0xxxffff(,,,)()0xxxffffXXXXX...XXXXXX...XX(,,,)(xxXXXX...(nnn00))0xX（8）将上式表示成矩阵形式，得11112(0)(0)(0)(0)12n1(0)(0)(0)(0)12nn00012222(0)(0212n000nnnn12n000)(0)12nnfff...xxxf(,,,)f(,,,)fff...xxxXXXXXXXX...XXXX...f(,,,)fff...xxx(0)（9）式中等号右边的矩阵nnxf都是对于(0)(0)(0)12n,,X,XX的值。这一矩阵称为雅可比（JACOBI）矩阵。按上述得到的修正向量(0)(0)(0)12nX,X,...X后，得到如下关系(1)(0)(0)111(1)(0)(0)222(1)(0)(0)nnnXXXXXX...XXX（10）这比(0)(0)(0)12n,,X,XX更接近真实值。这一步在收敛到希望的值以前重复进行，一般要反复计算满足1112121111,,,maxnnnnnnnnXXXXXX（11）为预先规定的小正数，1nnX是第n次迭代nX的近似值。二、极坐标形式的牛顿-拉夫逊法潮流计算1.节点的分类(1)PQ节点：给定Pi及Qi，求Ui及δi;(2)PV节点：给定Pi和Ui，求δi及Qi;(3)平衡节点：给定Ui和δi，求Pi和Qi；其中，P有功功率，Q无功功率，U节点电压，δ功率因数角设系统中有n个节点，其中m个PQ节点，n-m+1个PV节点。为了叙述方便，假定按照先PQ节点，再PV节点，再平衡节点的次序进行编号，即节点编号·种类个数1,2，…mPQ节点mm+1,m+2,...n-1PV节点n-m-1n平衡节点1图1.节点电压的分类2．节点功率方程为了简化，假定所有节点有含发电机和负荷。对各类节点列出电压用极坐标形式表示的功率方程（1）PQ节点的功率方程GiLiiijijijijjiGiLiiijijijijjiPPU(GcosBsin)QQU(GsinBcos)i1,2,..,m（12）式中，GiLiGiLiP,P,Q,Q都是给定量，而iiU,都待求。（2）PV节点的功率方程GiLiiijijijijjiGiLiiijijijijjiPPU(GcosBsin)QQU(GsinBcos)im1,...n1（13）式中，GiLiiP,P,U都是给定值，而GiLiiQ,Q,是待求量。（3）平衡节点GiLiiijijijijjiGiLiiijijijijjiPPU(GcosBsin)QQU(GsinBcos)in（14）式中，iiU,都是已知量，而GiLiGiLiP,P,Q,Q待求量将以上功率方程改写为iGiLiiijijijijjiiGiLiiijijijijjiP(x)PPU(GcosBsin)Q(x)QQU(GsinBcos)（15）从而可以得到修正方程(简写式)(k)(k)(k)(k)(k)(k)(k)(k)(k)(k)(k)PHNJQMLU'U'（16）其中T(k)(k)(k)(k)1n1T(k)(k)(k)(k)(k)(k)(k)1122mm2...U'U/UU/U...U/U（17）式（16）中，对于每个模块矩阵的对角线元素（ji）有iijijijijijijjiijjijijijijijjiijijijijijijjiijjijijijijijjPHUU(GsinBcos)PNUUU(GcosBsin)UQMUU(GcosBsin)QLUUU(GsinBcos)U(18)对于每个模块矩阵的对角线元素（ji）有iiiijijijijijjijjiiiiiijijijijijjiijiiiiijijijijijjiijiiiiiijijijijijjiijiPHUU(GsinBcos)PNUUU(GcosBsin)UQMUU(GcosBsin)QLUUU(GsinBcos)U(19)对修正方程（16）进行求解后，可以得出修正量(k)和(k)U'(k1)(k)(k)iii(k)(k1)(k)(k)iiii(k)i;i1,2,...n1UUUU;i1,2,...mU（20）3.初值的给定对于PQ节点电压有效值得初值，一般都给定为1，因为正常运行情况下各母线的电压都在额定值附近。至于PV节点和平衡节点，其电压的有效值都是给定不变的。4.牛顿法潮流计算的步骤及程序框图（1）用牛顿法计算潮流时，有以下的步骤：1.输入系统原始数据。包括节点分类和节点编号；2.形成节点导纳矩阵；3.给这各节点电压初始值(0)(0)U,；4.置迭代次数k=0；5.计算有功功率误差；6.判断误差是不是在范围内。如果是，则转向第11步，否则进行下一步；7.计算雅可比矩阵元素；8.解修正方程，得到'(k)(k)U,；9.计算各节点电压和相位的修正值，即新的初值；10.置k=k+1，返回第5步进行下一轮迭代；11.计算平衡节点注入有功功率和无功功率，计算各PV节点注入无功功率，计算元件两端的功率、电流、损耗，最后输出结果。（2）程序框图如下：启动输入原始数据形成节点导纳矩阵给定电压幅值和相位初值和迭代次数K=0按式（15）计算,求得雅克比矩阵各元素计算PV节点和平衡节点的注入功率，计算各支路的功率，电流，损耗求解修正方程，得求输出结果置k=k+1否是结束三．实例分析下图所示的简单电力系统中，系统中节点1、2为PQ节点，节点3为PV节点，节点4为平衡节点，已给定3.04.01js，2.03.02js，4.03P，02.13V，05.14V，04，网络各元件参数的标幺值如表2所示，给定电压的初始值如表2所示，收敛系数00001.0~123444V11jQP22jQP3V3P网络各元件参数的标幺值支路电阻电抗输电线路cy21变压器变比k1—20.030.090.02—1—30.020.050.02—2—30.040.08——2—40.00.05——3—40.030.07——各节点电压（初值）标幺值参数1.计算程序disp('电力系统极坐标下的牛顿-拉夫逊法潮流计算:');clearn=input('请输入结点数：n=');n1=input('请输入PV结点数：n1=');n2=input('请输入PQ结点数：n2=');节点i1234电压1.00+j0.01.0+j0.01.0+j0.01.05+j0.0isb=input('请输入平衡结点：isb=');pr=input('请输入精确度：pr=');K=input('请输入变比矩阵:K=');C=input('请输入支路阻抗矩阵：C=');y=input('请输入支路导纳矩阵：y=');U=input('请输入结点电压矩阵：U=');S=input('请输入各结点的功率：S=');Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros(1,n1+n2);form=1:nforR=1:nC(m,m)=C(m,m)+y(m,R);ifK(m,R)~=0C(m,m)=C(m,m)+1/(C(m,R)/(K(m,R)*(K(m,R)-1)));C(R,R)=C(R,R)+1/(C(m,R)/(1-K(m,R)));C(m,R)=C(m,R)/K(m,R);C(R,m)=C(m,R);endendendform=1:nforR=1:nifm~=RZ(m)=Z(m)+1/C(m,R);endendendform=1:nforR=1:nifm==RY(m,m)=C(m,m)+Z(m);elseY(m,R)=-1/C(m,R);endendenddisp('结点导纳矩阵:');disp(Y);disp('迭代中的雅克比矩阵:');G=real(Y);B=imag(Y);O=angle(U);U1=abs(U);k=0;PR=1;P=real(S);Q=imag(S);whilePRprform=1:n2UD(m)=U1(m);endform=1:n1+n2forR=1:nPT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endPT1(m)=sum(PT);PP(m)=P(m)-PT1(m);PP1(k+1,m)=PP(m);endform=1:n2forR=1:nQT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m