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

0引言
潮流是配电网络分析的基础，用于电网调度、运行分析、操作模拟和设计规划，同时也是电压优化和网络接线变化所要参考的内容。潮流计算通过数值仿真的方法把电力系统的详细运行情况呈现给工作人员，从而便于研究系统在给定条件下的稳态运行特点。随着市场经济的发展，经济利益是企业十分看重的，而线损却是现阶段阻碍企业提高效益的一大因素。及时、准确的潮流计算结果，可以给出配电网的潮流分布、理论线损及其在网络中的分布，从而为配电网的安全经济运行提供参考。从数学的角度来看，牛顿-拉夫逊法能有效进行非线性代数方程组的计算且具有二次收敛的特点，具有收敛快、精度高的特点，在输电网中得到广泛应用。随着现代计算机技术的发展，利用编程和相关软件，可以更好、更快地实现配电网功能，本文就是结合牛顿-拉夫逊法的基本原理，利用C++程序进行潮流计算，计算结果表明该方法具有良好的收敛性、可靠性及正确性。
1牛顿-拉夫逊法基本介绍
1.1潮流方程
对于N个节点的电力网络（地作为参考节点不包括在内），如果网络结构和元件参数已知，则网络方程可表示为：
（1-1）
式中，为N*N阶节点导纳矩阵；为N*1维节点电压列向量；为N*1维节点注入电流列向量。如果不计网络元件的非线性，也不考虑移相变压器，则为对称矩阵。
电力系统计算中，给定的运行变量是节点注入功率，而不是节点注入电流，这两者之间有如下关系：
（1-2）
式中，为节点的注入复功率，是N*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牛拉法的一般描述
牛拉法是把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，即非线性问题通过线性化逐步近似，这就是牛拉法的核心。下面以非线性方程式的求解过程来进行说明。
设电力网络的节点功率方程一般形式如下：
（1-7）
式中，为节点注入功率给定值；为对应的物理量和节点电压之间的函数表达式；为节点电压。写成功率偏差的形式：
（1-8）
应用牛拉法求解如下。在给定的初值处将式（1-8）作一阶泰勒展开：

定义为潮流方程的雅克比矩阵，为在处的值，则有：

用修正就得到的新值。如果用k表示迭代次数，写成一般的表达式，有：
（1-9）
对于潮流收敛的情况，应比更接近于解点。收敛条件为：
由简单迭代法收敛性分析的结论知，越接近解点，牛顿-拉夫逊法收敛越快，它具有二阶收敛速度。由图1.1可以直观地了解牛拉法的步骤：

图1.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）用下式表示：

其中每个字块的计算公式如下：
对角元素： （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）；
④ 将节点电压初值代入式(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。第一个零下面三行数为支路参数，分别表示三条支路的起始和终止节点编号，后面的为电阻、电抗和电纳，电导均为0，例如：1 2 0.1 0.4 0.01528。第二个零下面的为变压器支路，各数字意义同支路参数。接下去三行均为节点参数，分别表示注入有功功率和无功功率。
调用text.cpp文件，得到运行结果，见图3.3和图3.4。

图3.3 运行结果1

图3.4运行结果2
3.3结果分析
将上述仿真结果整理为表格3.1、3.2，其中“+”表示节点i输出功率给节点j，“-”表示节点j输出功率给i(纵向为i，横向为j)。
表3.1节点有功功率输入与输出
节点号
1
2
3
4
1
0
+0.245981
-0.5
-0.046563
2
-0.24316
0
0
-0.312949
3
0.5
0
0
0
4
0.0482143
0.319671
0
0
表3.2节点无功功率输入与输出
节点号
1
2
3
4
1
0
-0.014708
-0.029001
-0.136187
2
0.0110505
0
0
-0.14036
3
0.097016
0
0
0
4
0.10464
0.160255
0
0
根据表格计算：
节点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.097016+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结论
通过上面的分析与计算，验证了程序的正确性。由于编写过程的不足，线路损耗没能直接计算出来，而是需要手算，比较遗憾。程序在运行过程中，需要区分三种不同的节点，这由子程序保证实现。相比于快速分解法，牛拉法程序较为复杂，但更精确一点，潮流误差较小。
4总结
本文基于牛顿-拉夫逊潮流算法的基本原理，利用C++编程计算了一个4节点简单电力网络的潮流，并验证了运行成果，误差在允许范围之内。因为牛拉法计算过程中要不断生成新的雅各布矩阵，所以相对来说占用内存较多，但收敛速度快，这在程序运行过程中可以体现出来。本文程序并不是特别实用，因为真正的电力网络不可能只有几个节点，而且各种电力设备的情况也会复杂很多，因此程序会变得非常大，占用极大内存。但是，我还是通过这次练习，进一步巩固了书本上的理论知识，了解了实际操作的过程步骤。最后，感谢杨伟老师的悉心指导！

参考文献

[1] 朱红，赵琦，王庆宝. C++程序设计教程[M]. 北京：清华大学出版社，2010.
[2]张伯明，陈寿孙，严正.高等电力网络分析[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.
附录

//***************************************************//
#include"iostream.h"
#include"math.h"
#include<fstream>
#include<string>
#include<cassert>
#define nodeNumber readParameter[6]
using namespace std;
///////////////////////////////////////////
class powerFlowCalculation //导纳矩阵的计算类
{
private:
int readDataAmount; //存放元件参数的读取个数
int balanceNodeindex; //平衡结点号
double balanceNodeVoltage; //平衡节点电压幅值
double balanceNodePowerP; //平衡节点功率
double balanceNodePowerQ;
double calculationAccuracy; //计算精度
double readParameter[7]; //宏定义
double *conductance;
double *susceptance ;
double *admittanceAmplitude;
double *admittanceAngle;
double *lineData;
int PVNodeNumber; //PV节点数
double *jacobiMatrix; //雅克比矩阵
double *PQData; //矩阵，存放PQ节点的数据
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(double para[],double*G,double*B)
{
double GIJ,BIJ;
int nii,njj,nij,nji;
GIJ=para[3]/(para[3]*para[3]+para[4]*para[4]);
BIJ=-para[4]/(para[3]*para[3]+para[4]*para[4]);
nij=converIndex(para[1],para[2]);
nji=converIndex(para[2],para[1]);
nii=converIndex(para[1],para[1]);
njj=converIndex(para[2],para[2]);
G[nji]=-GIJ;
G[nij]=-GIJ;
B[nji]=-BIJ;
B[nij]=-BIJ;
G[nii]+=GIJ;
G[njj]+=GIJ;
B[nii]+= (BIJ+para[5]);
B[njj]+= (BIJ+para[5]);
}
//变压器的计算，规定para[1]为理想变压器侧，para[2]为归算阻抗侧
void countTranformer(double para[],double*G,double*B)
{
double GIJ,BIJ,b_i,b_j;
int nii,njj,nij,nji;
GIJ=para[3]/(para[3]*para[3]+para[4]*para[4]);
BIJ=-para[4]/(para[3]*para[3]+para[4]*para[4]);
b_i=BIJ*(1-para[5])/para[5]/para[5];
b_j=BIJ*(para[5]-1)/para[5];
nij=converIndex(para[1],para[2]);
nji=converIndex(para[2],para[1]);
nii=converIndex(para[1],para[1]);
njj=converIndex(para[2],para[2]);
G[nij]=-GIJ/para[5];
B[nij]=-BIJ/para[5];
G[nji]=-GIJ/para[5];
B[nji]=-BIJ/para[5];
G[nii]+= GIJ/para[5];
B[nii]+= (BIJ/para[5]+b_i);
G[njj]+= GIJ/para[5];
B[njj]+= (BIJ/para[5]+b_j);
}
//对地支路的计算
void countGroundBranch(double para[],double*G,double*B)
{
int nii;
nii=converIndex(para[1],para[1]);
G[nii]+=para[2];
B[nii]+=para[3];
}
//导纳矩阵的显示
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)
cout<<conductance[nij]<<susceptance [nij]<<"j"<<" ";
else
cout<<conductance[nij]<<"+"<<susceptance [nij]<<"j"<<" ";
}
cout<<endl;
}
cout<<endl;
}
//显示任意秩矩阵
void displayMatrix(double *matr,int ranki,int rankj)
{
int nij;
for(int i=1;i<=ranki;i++)
{
for(int j=1;j<=rankj;j++)
{
nij=(i-1)*rankj+j;
cout<<matr[nij]<<" ";
}
cout<<endl;
}
//cout<<endl;
}
//构造函数,打开文件
powerFlowCalculation(string input)
{
instream.open(input());
assert(instream.is_open());
}
//析构函数，关闭文件
~powerFlowCalculation()
{
instream.close();
}
//初始化，根据节点数开辟数组
void initi()
{
nodeNumber=readParameter[0]; //读取节点数
PVNodeNumber=0;
conductance=new double[nodeNumber*nodeNumber];
susceptance =new double[nodeNumber*nodeNumber];
admittanceAmplitude=new double[nodeNumber*nodeNumber];
admittanceAngle=new double[nodeNumber*nodeNumber];

PQData=new double[(nodeNumber-1)*6]; //开辟PQData
PVData=new double[(nodeNumber-1)*5]; //开辟PVData
voltageAmplitude=new double[nodeNumber];
voltageAngle=new double[nodeNumber];
jacobiMatrix=new double[nodeNumber*nodeNumber*4];
constantVector=new double[nodeNumber*2];
lineData=new double[nodeNumber*(nodeNumber-1)/2*5];
lineConsumeG=new double[nodeNumber*nodeNumber];
lineConsumeB=new double[nodeNumber*nodeNumber];
for(int a=1;a<=nodeNumber*(nodeNumber-1)/2*5;a++)
lineData[a]=0.0;
for(int i=1;i<=nodeNumber;i++) //电压向量的初始化
{
voltageAmplitude[i]=1.0;
voltageAngle[i]=0.0;
}
for(int j=1;j<=(nodeNumber*nodeNumber);j++) //导纳初始化
{
conductance[j]=0.0;
susceptance [j]=0.0;
admittanceAmplitude[j]=0.0;
admittanceAngle[j]=0.0;
lineConsumeG[j]=0.0;
lineConsumeB[j]=0.0;
}
for(int k=0;k<=nodeNumber*nodeNumber*4;k++) //雅克比矩阵初始化
jacobiMatrix[k]=0.0;
for(int n=1;n<=nodeNumber*2;n++) //常数向量初始化
{
constantVector[n]=0.0;
}
}
//计算导纳矩阵主程序
void countAdmittance()
{
int lineDataSerial=1;
for(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>>readParameter[0];
if (readParameter[0]<0.5) break; //序号为0时退出循环
for(int i=1;i<readDataAmount;i++)
{
instream>>readParameter[i];
if(facility==1)
{
lineData[lineDataSerial]=readParameter[i];
lineDataSerial++;
}
if(facility==2)
{
lineData[lineDataSerial]=readParameter[i];
if(i==5)
{
lineData[lineDataSerial]=0.0;
lineDataSerial--;
lineData[lineDataSerial]*=readParameter[i];
lineDataSerial--;
lineData[lineDataSerial]*=readParameter[i];
lineDataSerial+=2;
}//end for if(i==5)
lineDataSerial++;
}//end for if(facility==2)
}//end for loop i
switch(facility) //不同元件的计算
{
case 0:
initi();
balanceNodeindex=(int)readParameter[1];
balanceNodeVoltage=readParameter[2];
calculationAccuracy=readParameter[3];
cout<<"节点数:"<<nodeNumber<<endl
<<"平衡节点:"<<readParameter[1]<<endl
<<"平衡节点电压:"<<readParameter[2]<<endl
<<"计算精度:"<<readParameter[3]<<endl;
break;
case 1:
countLineBranch(readParameter,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)<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 converAdmittanceMatrix()
{
int n;
for(int i=1;i<=nodeNumber;i++)
{
for(int j=1;j<=nodeNumber;j++)
{
n=converIndex(i,j);
admittanceAmplitude[n]=sqrt(conductance[n]*conductance[n]+susceptance[n]*susceptance[n]);
admittanceAngle[n]=vectorAngle(conductance[n],susceptance[n]);
}
}
}
//读取PQ节点和PV节点的数据
void readNodePowerData()
{
for(int j=1;j<=(nodeNumber-1)*6;j++) //读取PQ节点功率数据存入PQData
instream>>PQData[j];
double dustbin; //将PQData与PVData之间的0去掉
instream>>dustbin;
for(int i=1;i<=(nodeNumber-1)*5;i++) //读取PV节点功率数据存入PVData
{
instream>>PVData[i];
if((i-1)%5==0)
{
if(PVData[i]==0) break;
PVNodeNumber+=1;
}
}//end for loop read PVData
}
//电压幅值和角度赋初值
void vectorAssignment() //PQ节点的值已经在初始化中赋值了，这里给平衡节点和PV节点赋值
{
int n;
voltageAmplitude[balanceNodeindex]=balanceNodeVoltage;
voltageAngle[balanceNodeindex]=0.0;
for(int i=1;i<=PVNodeNumber;i++)
{
n=PVData[converIndex(i,2)];
voltageAmplitude[n]=PVData[converIndex(i,3)];
voltageAngle[n]=0.0;
}
}
//生成L矩阵
double creatL(int i,int j)
{
double sum=0.0;
double alpha=0.0;
double beta=0.0;
if(i==j)
{
alpha=admittanceAngle[converIndex(i,i)];
for(int k=1;k<=nodeNumber;k++)
{
beta=-admittanceAngle[converIndex(i,k)];
sum+=admittanceAmplitude[converIndex(i,k)]*voltageAmplitude[k]*sin(beta);
}//end for loop k
sum*=voltageAmplitude[i];
sum=voltageAmplitude[i]*voltageAmplitude[i]*admittanceAmplitude[converIndex(i,i)]*sin(alpha);
sum=-sum;
}
else
{
beta=voltageAngle[i]-voltageAngle[j]-admittanceAngle[converIndex(i,j)];
sum=voltageAmplitude[i]*admittanceAmplitude[converIndex(i,j)]*voltageAmplitude[j]*sin(beta);
}
return sum;
}
//生成N矩阵
double creatN(int i,int j)
{
double sum=0.0;
double alpha=0.0;
double beta=0.0;
if(i==j)
{
alpha=admittanceAngle[converIndex(i,i)];
for(int k=1;k<=nodeNumber;k++)
{
beta=-admittanceAngle[converIndex(i,k)];
sum+=admittanceAmplitude[converIndex(i,k)]*voltageAmplitude[k]*cos(beta);
}//end for loop k
sum*=voltageAmplitude[i];
sum+=voltageAmplitude[i]*voltageAmplitude[i]*admittanceAmplitude[converIndex(i,i)]*cos(alpha);
sum=-sum;
}
else
{
beta=voltageAngle[i]-voltageAngle[j]-admittanceAngle[converIndex(i,j)];
sum=-voltageAmplitude[i]*admittanceAmplitude[converIndex(i,j)]*voltageAmplitude[j]*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) continue;
beta=voltageAngle[i]-voltageAngle[k]-admittanceAngle[converIndex(i,k)];
sum+=admittanceAmplitude[converIndex(i,k)]*voltageAmplitude[k]*cos(beta);
}//end for loop k
sum*=voltageAmplitude[i];
sum=-sum;
}
else
{
beta=voltageAngle[i]-voltageAngle[j]-admittanceAngle[converIndex(i,j)];
sum=voltageAmplitude[i]*admittanceAmplitude[converIndex(i,j)]*voltageAmplitude[j]*cos(beta);
}//end for if
return 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=voltageAngle[i]-voltageAngle[k]-admittanceAngle[converIndex(i,k)];
sum+=admittanceAmplitude[converIndex(i,k)]*voltageAmplitude[k]*sin(beta);
}
sum*=voltageAmplitude[i];
}
else
{
beta=voltageAngle[i]-voltageAngle[j]-admittanceAngle[converIndex(i,j)];
sum=-1*voltageAmplitude[i]*admittanceAmplitude[converIndex(i,j)]*voltageAmplitude[j]*sin(beta);
}//end for if
return sum;
}
//生成雅各比矩阵
void creatJacobiMatrix()
{
int serial;
int serialH;
int serialN;
int serialJ;
int 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)*node