直角坐标系下牛顿法潮流计算(共27页).doc

《直角坐标系下牛顿法潮流计算(共27页).doc》由会员分享，可在线阅读，更多相关《直角坐标系下牛顿法潮流计算(共27页).doc（27页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上直角坐标系下牛顿法潮流计算1电力系统潮流计算潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。2节点导纳矩阵的形成在图1（a）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b）所示。将接于节点1

2、和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c）的等值网络，其中和分别称为节点1和4的注入电流源。 图1 电力系统及其网络以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下： （2-1）上述方程组经过整理可以写成 （2-2）式中，；。一般的，对于有个独立节点的网络，可以列写个节点方程 （2-3）也可以用矩阵写成 （2-4）或缩写为 （2-5）矩阵称为节点导纳矩阵。它的对角线元素称为节点的自导纳，其值等于接于节点的所有支路导纳之和。非对角线元素称为节点、 间的互导纳，它等于直接接于节点、间的支路导纳的负值。若节点、间不存在

3、直接支路，则有。由此可知节点导纳矩阵是一个稀疏的对称矩阵。3牛顿-拉夫逊法潮流计算3.1牛顿-拉夫逊法的基本原理牛顿拉夫逊法（NewtonRaphson法）是求解非线性方程代数方程组的有效迭代计算方法。在牛顿拉夫逊法的每一次迭代过程中，对非线性方程通过线性化处理逐步近似。下面以单变量加以说明。设有单变量非线性方程 (3-1)求解此方程时。先给出解的近似值它与真解的误差为，则将满足方程，即 (3-2)将(3-8)式左边的函数在附近展成泰勒级数，于是便得 （3-3）式中,分别为函数在处的一阶导数，.，n阶导数。如果差值很小，3-9式右端的二次及以上阶次的各项均可略去。于是，3-9便简化为 0 (3

4、-4)这是对于变量的修正量的现行方程式，亦称修正方程式。解此方程可得修正量 (3-5)用所求的去修正近似解，变得 (3-6)由于3-10是略去高次项的简化式，因此所解出的修正量也只是近似值。修正后的近似解同真解仍然有误差。但是，这样的迭代计算可以反复进行下去，迭代计算的通式是 (3-7)迭代过程的收敛判据为 (3-8)或 (3-9)式中，为预先给定的小正数。 这种解法的几何意义可以从图31得到说明。函数yf(x)为图中的曲线。f(x)0的解相当于曲线与x轴的交点。如果第k次迭代中得到，则过点作一切线，此切线同x轴的交点便确定了下一个近似值。由此可见，牛顿拉夫逊法实质上就是切线法，是一种逐步线性

5、化的方法。 应用牛顿法求解多变量非线性方程组3-1时，假定已给出各变量的初值，. ，令，. 分别为各变量的修正量，使其满足方程3-2即 (3-10)将上式中的n个多元函数在初始值附近分别展成泰勒级数,并略去含有，二次及以上阶次的各项，便得. (3-11)方程式3-17也可以写成矩阵形式 (3-12)方程式3-18是对于修正量， 的线性方程组,称为牛顿法的修正方程式.利用高斯消去法或三角分解法可以解出修正量，。然后对初始近似值进行修正 (i=1,2,.,n) (3-13)如此反复迭代，在进行k1次迭代时，从求解修正方程式 (3-14)得到修正量，并对各变量进行修正 (i=1,2,n) (3-15

6、)式3-20和3-21也可以缩写为 (3-16)和 (3-17) 式中的X和分别是由n个变量和修正量组成的n维列向量；F(X)是由n个多元函数组成的n维列项量；J是n阶方阵，称为雅可比矩阵，它的第i、j个元素是第n个函数对第j个变量的偏导数；上角标(k)表示阵的每一个元素都在点处取值。迭代过程一直到满足收敛判据 (3-18)或 (3-19)为止。和为预先给定的小正数。 将牛顿拉夫逊法用于潮流计算，要求将潮流方程写成形如方程式3-1的形式。由于节点电压可以采用不同的坐标系表示，牛顿拉夫逊法潮流计算也将相应的采用不同的计算公式。3.2节点电压用直角坐标表示是的牛顿-拉夫逊法潮流计算采用直角坐标时，

7、节点电压可表示为导纳矩阵元素则表示为将上述表示式代入的右端，展开并分出实部和虚部，便得（11-45）假定系统中的第1，2，3，m号节点为PQ节点，第i个节点的给定功率设为和，对对该节点可列写方程 （i=1,2，m） （11-46）假定系统中的第m+1，m+2，n-1号节点为PV节点，则对其中每一个节点可以列写方程 （i=m+1，m+2，n-1） （11-47）第n号节点为平衡点，其电压是给定的，故不参加迭代。式（11-46）和式（11-47）总共包含了2（n-1）个方程，待求的变量有也是2（n-1）个。我们还可看到，方程（11-46）和式（11-47）已经具备了方程组（11-34）的形式。因此

8、，不难写出如下的修正方程式 （11-48）式中上述方程中雅克比矩阵的各元素，可以对式（11-46）和式（11-47）求偏导数获得。当时 （11-49）当时 （11-50）修正方程式（11-48）还可以写成分块矩阵的形式 （11-51）式中，和都是二维列向量；是介方阵。对于PQ节点 （11-52）对于PV节点 （11-53）从表达式（1-7）（1-11）可以看到，雅克比矩阵有以下特点：（1） 雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。（2） 雅克比矩阵的子块中的元素的表达式只用到导纳矩阵中的对应元素。若，则必有。因此，式（1-9）式中分块形式的雅克比矩阵同节点导纳矩

9、阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧。（3） 无论在式（1-6）或式（1-9）中雅克比矩阵的元素或子块都不具有对称性。用牛顿-拉夫逊法计算潮流的流程框图示于图11-27。首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值和，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：（1） 按上一次迭代计算出的节点电压值和，利用式（11-46）和式（11-47）计算各类节点的不平衡量、和。（2） 按条件（11-43）校验收敛，即 如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续

10、计算。（3）利用式（11-49）和式（11-50）计算雅克比矩阵的各元素。（4）解修正方程式（11-48）求节点电压的修正量和。（5）修正各节点的电压（6） 迭代计数加1，返回第一步继续迭代过程。迭代结束后，还要计算出平衡节点的功率和网络中的功率分布。输电线路功率的计算公式如下：3.3牛顿-拉夫逊法潮流计算程序框图否是输出输入原始数据计算平衡节点功率及全部线路功率k=k+1解修正方程式，求计算雅可比矩阵各元素是否收敛计算及k=0给定节点电压初值形成节点导纳矩阵图2 牛顿-拉夫逊法潮流计算程序框图4 C语言实现直角坐标表示的牛顿拉夫逊的计算4.1计算节点导纳矩阵该电力系统的等值电路如图所示图3电

11、力系统等值电路根据节点方程和等值电路图，计算节点导纳矩阵的C语言子程序如下。G11=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r* z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12

12、r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120;G33=z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*

13、z12r+z12m*z12m);G13=G31=-k*z13r/(z13r*z13r+z13m*z13m);B13=B31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i5;i+)for(j

14、=1;j5;j+)printf(%f+%fj,Gij,Bij);printf( );printf(n);4.2计算各节点不平衡量按上一次迭代计算出的节点电压值和（当k=0时即为给定值），利用式（1-2和式（1-3）计算各类节点的不平衡量、和。同样的可以算出其部分C语言程序如下：printf(迭代次数k1=%dn,k1);for (i=1;i5;i+)float a=0,b=0;for(j=1;j5;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+Bij*ej;Pi=Psi-(ei*a+fi*b);/计算有功功率的增量Qi=Qsi-(fi*a-ei*b);/计算无功功率的增量V32

15、=V3s*V3s-e3*e3;printf(有功功率增量P1=%f,P1); printf( ,);printf(有功功率增量P2=%f,P2); printf( ,);printf(有功功率增量P3=%f,P3);printf(无功功率增量Q1=%f,Q1); printf( ,);printf(无功功率增量Q2=%f,Q2); printf( ,);printf(电压增量V32=%f,V32);printf(n);计算后校验收敛条件 fabs(P2)?fabs(P1):fabs(P2);max=maxfabs(P3)?max:fabs(P3);max=maxfabs(Q1)?max:fab

16、s(Q1);max=maxfabs(Q2)?max:fabs(Q2);max=maxfabs(V32)?max:fabs(V32);printf(max=%fn,max);/*while (max0.00001)经计算，各节点的不平衡量都未满足收敛条件，于是继续计算。4.3计算雅克比矩阵 按公式（2-1）和（2-2）计算雅克比矩阵的各元素，方程中雅克比矩阵的各元素，可以对式（1-2）和式（1-3）求偏导数获得。当时，当时，计算雅可比矩阵的C程序如下：for(j=1;j4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G1m*em-B1m*fm;

17、d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;4.4解修正方程 修正方程式 式中 修正各节点电压后，迭代计数加1，然后返回继续迭代过程。4.5计算平衡节点功率迭代结束后，还要算出平衡点的功率

18、和网络中的功率分布。输电线路功率的计算公式如下。计算平衡点功率和网络中功率分布的程序如下。p12=-2*e1*f1*y120-(e1*(e1-e2)-f1*(f1-f2)*G12+(e1*(f1-f2)+f1*(e1-e2)*B12; q12=-(e1*e1-f1*f1)*y120+(e1*(e1-e2)+f1*(f1-f2)*B12+(e1*(f1-f2) +f1*(e1-e2)*G12;5程序调试结果程序调试结果如图所示图4程序调试结果图总 结本设计采用直角坐标形式的牛顿拉夫逊法作常规潮流计算。P-Q分解法利用了电力系统的一些特有的运行特点，对牛顿拉夫逊法作了简化，可提高计算速度，但较难理

19、解，牛顿拉夫逊法的雅克比矩阵在每次迭代过程中都有变化，需要重新形成和求解，这占用了计算的大部分时间，成为牛顿拉夫逊计算速度不能提高的主要因素，但收敛性好，物理概念也较为清晰。应用C语言进行编程层次清晰，便于按模块化方式组织，易于调试和维护，C语言的表现能力和处理能力极强，极大的提高了计算效率，同时也培养了我们快速的掌握应用程序的能力和解决实际问题的能力。通过本次设计，我们更好的理解和掌握电力系统潮流计算的原理和方法，也更加系统的理解C程序设计语言的编程技巧。今后的研究方向仍围绕着如何改善传统潮流计算方法以及提高计算机潮流计算的可靠性、收敛性以及速度和灵活性。参考文献1 何仰赞等.电力系统分析上

20、册M武汉：华中科技大学出版社.2 何仰赞等.电力系统分析下册M武汉：华中科技大学出版社.3 陈珩.电力系统稳态分析. 水利电力出版社,19944 熊静琪. 计算机控制技术. 电子工业出版社，ISBN 7- 5053-8706- 5.5张伯明,陈寿孙.高等电力网络分析M.北京：清华大学出版社,1996.附录主程序如下：#include #include #include #define N 2#define M 7main()float z12r,z12m,y120,z13r,z13m,k,z14r,z14m,y140,z24r,z24m,y240,G55,B55,J77;float e5=0,

21、1,1,1.1,1.05,f5=0,P5,Q5,Ps5=0,-0.30,-0.55,0.5,xe4,xf4;float Qs5=0,-0.18,-0.13,V3s=1.1,V4S=1.05;float V32,max,P4,Q4;float a1=0,b1=0;int i,j,n,s,k1=0;float LMM=0,UMM=0,sigma1,sigma2,bM,yM,xM;float p12,p13,p14,p21,p24,p31,p41,p42,q12,q13,q14,q21,q24,q31,q41,q42;printf(请输入z12的实部和虚部n); /输入电路中的阻抗scanf(%f%

22、f,&z12r,&z12m);printf(请输入z13的实部和虚部n);scanf(%f%f,&z13r,&z13m);printf(请输入z14的实部和虚部n);scanf(%f%f,&z14r,&z14m);printf(请输入z24的实部和虚部n);scanf(%f%f,&z24r,&z24m);printf(请输入y120的值n);scanf(%f,&y120);printf(请输入y140的值n);scanf(%f,&y140);printf(请输入y240的值n);scanf(%f,&y240);printf(请输入变比k的值n);scanf(%f,&k);/*计算导纳矩阵*G1

23、1=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z2

24、4m)+y240+y120;G33=z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*z12r+z12m*z12m);G13=G31=-k*z13r/(z13r*z13

25、r+z13m*z13m);B13=B31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i5;i+)for(j=1;j5;j+)printf(%f+%fj,Gij,Bij);printf( )

26、;printf(n);/形成节点导纳矩阵/*printf(n);/*/计算各节点不平衡量loop1:printf(迭代次数k1=%dn,k1);for (i=1;i5;i+)float a=0,b=0;for(j=1;jfabs(P2)?fabs(P1):fabs(P2);max=maxfabs(P3)?max:fabs(P3);max=maxfabs(Q1)?max:fabs(Q1);max=maxfabs(Q2)?max:fabs(Q2);max=maxfabs(V32)?max:fabs(V32);printf(max=%fn,max);/*while (max0.00001)/*形成雅

27、克比矩阵*for(j=1;j4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G

28、1j*f1;for(j=1;j4;j+)if(2=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G2m*em-B2m*fm; d+=G2m*fm+B2m*em;J2*N-1j*N-1=-c-G2j*e2-B2j*f2;J2*N-1j*N=-d+B2j*e2-G2j*f2;J2*Nj*N-1=d+B2j*e2-G2j*f2;J2*Nj*N=-c+G2j*e2+B2j*f2;elseJ2*N-1j*N-1=-G2j*e2-B2j*f2;J2*Nj*N=G2j*e2-B2j*f2;J2*N-1j*N=B2j*e2-G2j*f2;J2*Nj*N-1=B2j*e2-G2j

29、*f2;for(j=1;j4;j+)if(3=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G3m*em-B3m*fm; d+=G3m*fm+B3m*em;J3*N-1j*N-1=-c-G3j*e3-B3j*f3;J3*N-1j*N=-d+B3j*e3-G3j*f3;J3*Nj*N-1=-2*e3;J3*Nj*N=-2*f3;elseJ3*N-1j*N-1=-G3j*e3-B3j*f3; J3*N-1j*N=B3j*e3-G3j*f3; J3*Nj*N-1=0; J3*Nj*N=0;printf(雅克比矩阵是:n);for(i=1;i7;i+)for(j=1;j

30、7;j+)printf(%f,Jij);printf( );printf(n);/*计算修正方程*for(i=1;iM;i+)Lii=1;for(i=1;iM;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;nM;n+)for(j=n;jM;j+)sigma1=0;for(s=0;s=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma1;for(i=n;iM;i+)sigma2=0;for(s=0;s=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P2;b4=Q2;b5=P3;b6=V3

31、2;for(i=1;iM;i+)sigma1=0;for(n=1;n=1;i-)sigma2=0;for(n=i+1;nM;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3;xe3=-x5;xf1=-x2;xf2=-x4;xf3=-x6;printf(节点电压:n);for(i=1;i4;i+)ei+=xei; fi+=xfi;for(i=1;i4;i+)printf(e%d=,i); printf(%f,ei); printf( ,);for(i=1;i4;i+)printf(f%d=,i); printf(%f,fi); printf

32、( ,);printf(n);k1=k1+1;goto loop1;for(j=1;j5;j+)a1+=G4j*ej-B4j*fj; b1+=G4j*fj+B4j*ej;P4=e4*a1+f4*b1;Q4=f4*a1-e4*b1;printf(P4+Q4=%f+j%f,P4,Q4);printf(n);p12=-2*e1*f1*y120-(e1*(e1-e2)-f1*(f1-f2)*G12+(e1*(f1-f2)+f1*(e1-e2)*B12;q12=-(e1*e1-f1*f1)*y120+(e1*(e1-e2)+f1*(f1-f2)*B12+(e1*(f1-f2)+f1*(e1-e2)*G12;p13=(e1*e1-f1*f1)*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+2*e1*f1*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-(e1*(e1-e3)-f1*(f1-f3)*G13+(e1*(f1-f3)+f1*(e1-e3)*B13;q13=(e1*e1-f1*f1)*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-2*e1*f1*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+(e1*(e