电力系统稳态分析-牛顿拉夫逊法

1、0 引言潮流是配电网络分析的基础，用于电网调度、运行分析、操作模拟和设计规划，同时也是电压优化和网络接线变化所要参考的内容。潮流计算通过数值仿真的方法把电力系统的详细运行情况呈现给工作人员，从而便于研究系统在给定条件下的稳态运行特点。随着市场经济的发展，经济利益是企业十分看重的，而线损却是现阶段阻碍企业提高效益的一大因素。及时、准确的潮流计算结果，可以给出配电网的潮流分布、理论线损及其在网络中的分布，从而为配电网的安全经济运行提供参考。从数学的角度来看，牛顿-拉夫逊法能有效进行非线性代数方程组的计算且具有二次收敛的特点，具有收敛快、精度高的特点，在输电网中得到广泛应用。随着现代计算机技术的发展

2、，利用编程和相关软件，可以更好、更快地实现配电网功能，本文就是结合牛顿-拉夫逊法的基本原理，利用C+程序进行潮流计算，计算结果表明该方法具有良好的收敛性、可靠性及正确性。1 牛顿-拉夫逊法基本介绍1.1 潮流方程对于N个节点的电力网络（地作为参考节点不包括在内），如果网络结构和元件参数已知，则网络方程可表示为： （1-1）式中，为N*N阶节点导纳矩阵；为N*1维节点电压列向量；为N*1维节点注入电流列向量。如果不计网络元件的非线性，也不考虑移相变压器，则为对称矩阵。电力系统计算中，给定的运行变量是节点注入功率，而不是节点注入电流，这两者之间有如下关系： （1-2）式中，为节点的注入复功率，是N

3、*1维列矢量；为的共轭；是由节点电压的共轭组成的N*N阶对角线矩阵。由（1-1）和（1-2），可得：上式就是潮流方程的复数形式，是N维的非线性复数代数方程组。将其展开，有： j=1,2，.，N （1-3）式中， 表示所有和相连的节点，包括。将节点电压用极坐标表示，即令，代入式（1-3）中则有：故有： i=1,2,N （1-4）式（1-4）是用极坐标表示的潮流方程。而节点功率误差： (1-5) (1-6)式中：，为节点给定的有功功率及无功功率。 1.2 牛顿-拉夫逊法基本原理1.2.1 牛拉法的一般描述牛拉法是把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，即非线性问题通过线性化逐

4、步近似，这就是牛拉法的核心。下面以非线性方程式的求解过程来进行说明。设电力网络的节点功率方程一般形式如下： （1-7）式中，为节点注入功率给定值；为对应的物理量和节点电压之间的函数表达式；为节点电压。写成功率偏差的形式： （1-8）应用牛拉法求解如下。在给定的初值处将式（1-8）作一阶泰勒展开：定义为潮流方程的雅克比矩阵，为在处的值，则有：用修正就得到的新值。如果用k表示迭代次数，写成一般的表达式，有： （1-9）对于潮流收敛的情况，应比更接近于解点。收敛条件为： 由简单迭代法收敛性分析的结论知，越接近解点，牛顿-拉夫逊法收敛越快，它具有二阶收敛速度。由图1.1可以直观地了解牛拉法的步骤：图1

5、.1 牛顿-拉夫逊法的几何解释1.2.2 极坐标的牛顿-拉夫逊法在极坐标中，有如下的形式： （1-10）共2n-r个方程，状态变量为：共2n-r个待求量。r 个PV节点的电压幅值给定，不需求解。潮流雅克比矩阵的维数是（2n-r）*(2n-r)，结构如下：上式右侧的对电压幅值的偏导数项中的电压幅值的阶数减少了1，为使雅克比矩阵的各部分子矩阵具有一致的形式，在实际计算中，常将该项乘以电压幅值，并选取作为待求的修正量，则雅克比矩阵可写成： （1-11）将式（1-10）和（1-11）代入式（1-9）的修正方程即可求得x的修正量，用它修正x直到为止。将式（1-11）用下式表示：其中每个字块的计算公式如下

6、：对角元素： （1-12）非对角元素： （1-13）2 牛顿法潮流计算步骤2.1 程序流程图在了解了牛拉法的原理之后，明确程序编写思路，如图2.1、2.2所示。 其中图2.1中的“计算电压幅值和角度”步骤较多，单独用图2.2表示出来。图2.1 牛顿法计算潮流的程序框图 图2.2 电压幅值和角度求解步骤框图当不符合收敛的条件“amontk>1”时,即认为计算不收敛。具体程序见附录。2.2 计算步骤下面讨论的是极坐标形式的牛顿法求解过程，大致分为以下几个步骤： 形成节点导纳矩阵； 给各节点电压设初值（）；根据式（1-12）、（1-13）生成雅克比矩阵（H、N、M、L）； 将节点电压初值代入式

7、(1-5)、式(1-6)，求出修正方程式的常数项向量； 求解修正方程，得到电压幅值和角度；判断是否收敛，若收敛，计算平衡节点和线路功率； 输出结果，并结束。3 算例3.1 系统模型本文以图3.1所示电力网络为例，调用基于牛顿-拉夫逊法的C+程序。图3.1 系统模型其中节点4设为平衡节点，电压标幺值为1.05，计算误差为0.000001。3.2 输入与输出将图3.1所示模型的相关数据放在data.dat文件中图3.2 输入节点和支路数据对各个数字含义的解释如下：网络模型有四个节点，四条支路，编号见图3.1。第一个零下面三行数为支路参数，分别表示三条支路的起始和终止节点编号，后面的为电阻、电抗和电

8、纳，电导均为0，例如：1 2 0.1 0.4 0.01528。第二个零下面的为变压器支路，各数字意义同支路参数。接下去三行均为节点参数，分别表示注入有功功率和无功功率。调用text.cpp文件，得到运行结果，见图3.3和图3.4。图3.3 运行结果1图3.4 运行结果23.3 结果分析将上述仿真结果整理为表格3.1、3.2，其中“+”表示节点i输出功率给节点j，“-”表示节点j输出功率给i(纵向为i，横向为j)。表3.1节点有功功率输入与输出节点号123410+0.245981-0.5-0.0465632-0.2431600-0.31294930.500040.04821430.3196710

9、0表3.2节点无功功率输入与输出节点号123410-0.014708-0.029001-0.13618720.011050500-0.1403630.09701600040.104640.16025500根据表格计算：节点1有功功率：0+0.245981-0.5-0.046563=-0.300582 无功功率：0-0.014708-0.029001-0.136187=-0.179896节点2有功功率：-0.24316+0+0-0.312949=-0.556109无功功率：0.0110505+0+0-0.14036=-0.1293095节点3有功功率：0.5+0+0+0=0.5无功功率：0.09

10、7016+0+0+0=0.097016节点4有功功率：0.0482143+0.319671+0+0=0.3678853无功功率：0.10464+0.160255+0+0=0.264895根据已知条件，两个PQ节点的注入有功、无功分别为：P1=0.3，Q1=0.18；P2=0.55，Q2=0.13潮流计算误差：可见，误差均在允许范围内。线路损耗： 3.4 结论通过上面的分析与计算，验证了程序的正确性。由于编写过程的不足，线路损耗没能直接计算出来，而是需要手算，比较遗憾。程序在运行过程中，需要区分三种不同的节点，这由子程序保证实现。相比于快速分解法，牛拉法程序较为复杂，但更精确一点，潮流误差较小。

11、4 总结本文基于牛顿-拉夫逊潮流算法的基本原理，利用C+编程计算了一个4节点简单电力网络的潮流，并验证了运行成果，误差在允许范围之内。因为牛拉法计算过程中要不断生成新的雅各布矩阵，所以相对来说占用内存较多，但收敛速度快，这在程序运行过程中可以体现出来。本文程序并不是特别实用，因为真正的电力网络不可能只有几个节点，而且各种电力设备的情况也会复杂很多，因此程序会变得非常大，占用极大内存。但是，我还是通过这次练习，进一步巩固了书本上的理论知识，了解了实际操作的过程步骤。最后，感谢杨伟老师的悉心指导！参考文献1 朱红，赵琦，王庆宝. C+程序设计教程M. 北京：清华大学出版社，2010.2张伯明，陈寿

12、孙，严正.高等电力网络分析M.北京：清华大学出版社，2007.3吴明波. 牛顿-拉夫逊法在潮流计算中的应用J. 内蒙古科技与经济，2011, 21：111-112，115. 4明日科技. VisualC+从入门到精通M. 北京：清华大学出版社，2012.5 顾洁，陈章潮，徐蓓. 一种新的配电网潮流算法-改进牛顿拉夫逊法J.华东电力，2000，5:10-12.6 朱文强. 牛顿-拉夫逊法在配电网中的应用J. 水利科技，2004，3:55-56,58.7 刘明波，谢敏，赵维新. 大电网最优潮流计算M. 北京：科学出版社，2010.31附录/*/#include"iostream.h&qu

13、ot;#include"math.h"#include<fstream>#include<string>#include<cassert>#define nodeNumber readParameter6using namespace std;/class powerFlowCalculation /导纳矩阵的计算类private: int readDataAmount; /存放元件参数的读取个数int balanceNodeindex; /平衡结点号double balanceNodeVoltage; /平衡节点电压幅值double b

14、alanceNodePowerP; /平衡节点功率double balanceNodePowerQ;double calculationAccuracy; /计算精度double readParameter7; /宏定义double *conductance;double *susceptance ; double *admittanceAmplitude;double *admittanceAngle;double *lineData;int PVNodeNumber; /PV节点数double *jacobiMatrix; /雅克比矩阵double *PQData; /矩阵，存放PQ节点的

15、数据 double *PVData; /矩阵，存放PV节点的数据double *voltageAmplitude; /电压幅值double *voltageAngle; /电压角度double *constantVector; /常数向量double *lineConsumeG;double *lineConsumeB; ifstream instream;public:/下标的转换int converIndex(int i,int j) int serial;serial=(i-1)*nodeNumber+j;return serial;/线路的计算void countLineBranch(

16、double para,double*G,double*B)double GIJ,BIJ;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4); nij=converIndex(para1,para2);nji=converIndex(para2,para1);nii=converIndex(para1,para1);njj=converIndex(para2,para2);Gnji=-GIJ;Gnij=-GIJ;Bnji=-BIJ;Bnij=-BIJ;Gnii+

17、=GIJ;Gnjj+=GIJ;Bnii+= (BIJ+para5); Bnjj+= (BIJ+para5); /变压器的计算，规定para1为理想变压器侧，para2为归算阻抗侧void countTranformer(double para,double*G,double*B)double GIJ,BIJ,b_i,b_j;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4);b_i=BIJ*(1-para5)/para5/para5; b_j=BIJ*(para

18、5-1)/para5;nij=converIndex(para1,para2);nji=converIndex(para2,para1); nii=converIndex(para1,para1);njj=converIndex(para2,para2); Gnij=-GIJ/para5;Bnij=-BIJ/para5;Gnji=-GIJ/para5;Bnji=-BIJ/para5;Gnii+= GIJ/para5;Bnii+= (BIJ/para5+b_i);Gnjj+= GIJ/para5;Bnjj+= (BIJ/para5+b_j); /对地支路的计算void countGroundBr

19、anch(double para,double*G,double*B) int nii; nii=converIndex(para1,para1); Gnii+=para2;Bnii+=para3;/导纳矩阵的显示void display() cout<<endl<<"导纳矩阵为："<<endl;int nij; for(int i=1;i<=nodeNumber;i+) for(int j=1;j<=nodeNumber;j+) nij=converIndex(i,j); if(susceptance nij<0) c

20、out<<conductancenij<<susceptance nij<<"j"<<" " else cout<<conductancenij<<"+"<<susceptance nij<<"j"<<" " cout<<endl;cout<<endl;/显示任意秩矩阵void displayMatrix(double *matr,int ranki,int ran

21、kj) int nij; for(int i=1;i<=ranki;i+) for(int j=1;j<=rankj;j+) nij=(i-1)*rankj+j; cout<<matrnij<<" " cout<<endl;/cout<<endl;/构造函数,打开文件powerFlowCalculation(string inputfilename) instream.open(inputfilename.data();assert(instream.is_open();/析构函数，关闭文件 powerFlowCa

22、lculation() instream.close();/初始化，根据节点数开辟数组void initi() nodeNumber=readParameter0; /读取节点数 PVNodeNumber=0; conductance=new doublenodeNumber*nodeNumber; susceptance =new doublenodeNumber*nodeNumber; admittanceAmplitude=new doublenodeNumber*nodeNumber; admittanceAngle=new doublenodeNumber*nodeNumber; P

23、QData=new double(nodeNumber-1)*6; /开辟PQDataPVData=new double(nodeNumber-1)*5; /开辟PVDatavoltageAmplitude=new doublenodeNumber;voltageAngle=new doublenodeNumber;jacobiMatrix=new doublenodeNumber*nodeNumber*4;constantVector=new doublenodeNumber*2;lineData=new doublenodeNumber*(nodeNumber-1)/2*5;lineCon

24、sumeG=new doublenodeNumber*nodeNumber;lineConsumeB=new doublenodeNumber*nodeNumber;for(int a=1;a<=nodeNumber*(nodeNumber-1)/2*5;a+)lineDataa=0.0;for(int i=1;i<=nodeNumber;i+) /电压向量的初始化voltageAmplitudei=1.0;voltageAnglei=0.0; for(int j=1;j<=(nodeNumber*nodeNumber);j+) /导纳初始化 conductancej=0.0

25、;susceptance j=0.0;admittanceAmplitudej=0.0;admittanceAnglej=0.0;lineConsumeGj=0.0;lineConsumeBj=0.0; for(int k=0;k<=nodeNumber*nodeNumber*4;k+) /雅克比矩阵初始化 jacobiMatrixk=0.0; for(int n=1;n<=nodeNumber*2;n+) /常数向量初始化 constantVectorn=0.0; /计算导纳矩阵主程序void countAdmittance() int lineDataSerial=1; for

26、(int facility=0;facility<4;facility+) switch(facility) /不同元件参数的读取个数 case 0: /节点数，平衡节点，平衡节点电压，计算精度 readDataAmount=4; break; case 1: /线路 readDataAmount=6; break; case 2: /变压器readDataAmount=6; break;case 3: /对地支路readDataAmount=4; break;/end switch for(;) /参数的读取和操作 instream>>readParameter0; if

27、(readParameter0<0.5) break; /序号为0时退出循环for(int i=1;i<readDataAmount;i+) instream>>readParameteri; if(facility=1) lineDatalineDataSerial=readParameteri; lineDataSerial+; if(facility=2) lineDatalineDataSerial=readParameteri; if(i=5)lineDatalineDataSerial=0.0;lineDataSerial-;lineDatalineData

28、Serial*=readParameteri;lineDataSerial-;lineDatalineDataSerial*=readParameteri; lineDataSerial+=2; /end for if(i=5) lineDataSerial+; /end for if(facility=2) /end for loop i switch(facility) /不同元件的计算 case 0:initi(); balanceNodeindex=(int)readParameter1;balanceNodeVoltage=readParameter2;calculationAccu

29、racy=readParameter3;cout<<"节点数:"<<nodeNumber<<endl<<"平衡节点:"<<readParameter1<<endl<<"平衡节点电压:"<<readParameter2<<endl<<"计算精度:"<<readParameter3<<endl; break;case 1:countLineBranch(readParamete

30、r,conductance,susceptance );break;case 2:countTranformer(readParameter,conductance,susceptance );break;case 3:countGroundBranch(readParameter,conductance,susceptance );break; /end switch /end for loop /end for (facility)/end countAdmittance()/*雅克比矩阵*/double vectorAngle(double a,double b)if(fabs(a)&l

31、t;0.0000001 && fabs(b)<0.0000001) return 0.0; double alpha; alpha=atan(fabs(b)/fabs(a);if(a<=0 && b>=0)alpha=3.1415926-alpha;else if(a<=0 && b<=0)alpha=-3.1415926+alpha;else if(a>=0 && b<=0)alpha=-alpha; return alpha;/将直角坐标的导纳矩阵转换为极坐标形式void conver

32、AdmittanceMatrix() int n; for(int i=1;i<=nodeNumber;i+)for(int j=1;j<=nodeNumber;j+) n=converIndex(i,j);admittanceAmplituden=sqrt(conductancen*conductancen+susceptancen*susceptancen);admittanceAnglen=vectorAngle(conductancen,susceptancen);/读取PQ节点和PV节点的数据void readNodePowerData() for(int j=1;j&l

33、t;=(nodeNumber-1)*6;j+) /读取PQ节点功率数据存入PQDatainstream>>PQDataj; double dustbin; /将PQData与PVData之间的0去掉instream>>dustbin; for(int i=1;i<=(nodeNumber-1)*5;i+) /读取PV节点功率数据存入PVDatainstream>>PVDatai;if(i-1)%5=0)if(PVDatai=0) break;PVNodeNumber+=1;/end for loop read PVData/电压幅值和角度赋初值void

34、 vectorAssignment() /PQ节点的值已经在初始化中赋值了，这里给平衡节点和PV节点赋值int n;voltageAmplitudebalanceNodeindex=balanceNodeVoltage;voltageAnglebalanceNodeindex=0.0;for(int i=1;i<=PVNodeNumber;i+)n=PVDataconverIndex(i,2);voltageAmplituden=PVDataconverIndex(i,3);voltageAnglen=0.0;/生成L矩阵double creatL(int i,int j)double

35、sum=0.0;double alpha=0.0;double beta=0.0;if(i=j)alpha=admittanceAngleconverIndex(i,i);for(int k=1;k<=nodeNumber;k+)beta=-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*sin(beta);/end for loop ksum*=voltageAmplitudei;sum=voltageAmplitudei*voltageAmplitud

36、ei*admittanceAmplitudeconverIndex(i,i)*sin(alpha);sum=-sum;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*sin(beta); return sum;/生成N矩阵double creatN(int i,int j)double sum=0.0;double alpha=0.0;double bet

37、a=0.0;if(i=j)alpha=admittanceAngleconverIndex(i,i);for(int k=1;k<=nodeNumber;k+)beta=-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta);/end for loop k sum*=voltageAmplitudei;sum+=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex

38、(i,i)*cos(alpha); sum=-sum;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=-voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*cos(beta);return sum;/生成J矩阵double creatJ(int i,int j)double sum=0.0;double beta;if(i=j)for(int k=1;k<=nodeNumber;k+)if(k=i)

39、continue;beta=voltageAnglei-voltageAnglek-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta);/end for loop ksum*=voltageAmplitudei;sum=-sum;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitude

40、converIndex(i,j)*voltageAmplitudej*cos(beta);/end for ifreturn sum;/生成H矩阵double creatH(int i,int j)double sum=0.0;double beta;if(i=j)for(int k=1;k<=nodeNumber;k+)if(k=i) continue;beta=voltageAnglei-voltageAnglek-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmpli

41、tudek*sin(beta);sum*=voltageAmplitudei;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=-1*voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*sin(beta);/end for ifreturn sum;/生成雅各比矩阵void creatJacobiMatrix() int serial;int serialH;int serialN;int serialJ;i

42、nt serialL;for(int i=1;i<=2*(nodeNumber-1);i=i+2)for(int j=1;j<=2*(nodeNumber-1);j=j+2)serialH=(i-1)*nodeNumber*2+j;serialN=(i-1)*nodeNumber*2+j+1;serialJ=(i-1)*nodeNumber*2+j+nodeNumber*2; serialL=(i-1)*nodeNumber*2+j+nodeNumber*2+1;jacobiMatrixserialH=creatH(i+1)/2,(j+1)/2);jacobiMatrixseria

43、lN=creatN(i+1)/2,(j+1)/2);jacobiMatrixserialJ=creatJ(i+1)/2,(j+1)/2);jacobiMatrixserialL=creatL(i+1)/2,(j+1)/2);/end for loop j/end for loop ifor(int n=(nodeNumber-PVNodeNumber)*2;n<=nodeNumber*2;n=n+2)for(int m=1;m<=nodeNumber*2;m+)if(n=m)serial=(m-1)*nodeNumber*2+n;jacobiMatrixserial=1.0;els

44、eserial=(m-1)*nodeNumber*2+n;jacobiMatrixserial=0.0;serial=(n-1)*nodeNumber*2+m;jacobiMatrixserial=0.0;/end for if/end for loop m/end for loop nserial=nodeNumber*nodeNumber*4-nodeNumber*2-1;jacobiMatrixserial=1.0;/ rank为6的下标转换int converSerial(int i,int j,int rank=6) int serial;serial=(i-1)*rank+j;re

45、turn serial;/生成常数向量void creatConstantVector() int i; /i为P、Q的下标double sum=0.0;double beta=0.0;for(int n=1;n<=(nodeNumber-1)*2;n+) /n为向量的下标i=(n+1)/2; /i为节点号 sum=0.0;if(n+1)%2=1) /求Q点for(int j=1;j<=nodeNumber;j+)beta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum+=voltageAmplitudei*voltageAmplitudej*admittanceAmplitudeconverIndex(i,j)*sin(beta);/end for loop jif(n<=(nodeNumber-PVNodeNumber-1)*2)constantVectorn=-(PQDataconverSerial(i,4)-PQDataconverSerial(i,6)-sum);elseconstantVectorn=0.0;else /求P点for(int jj=1;jj<