用C语言进行复杂电力系统的潮流计算
用C语言进行复杂电力系统的潮流计算
杨汉明
电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤，是指在一定的系统结构和运行条件下，确定系统运行状态的计算，也即是对各母线（节点）电压，各元件（支路）传输电线或功率的计算。通过计算出的节点电压和功率分布用以检查系统各元件是否过负荷，各点电压是否合理，以及功率损耗等。
即使对于一个简单的电力系统，潮流计算也不是一件简单就可以完成的事，其运算量很大，因此如果对于一个大的、复杂的电网来说的话，由于其节点多，分支杂，其计算量可想而知，人工对其计算也更是难上加难了。特别是在现实生活中，遇到一个电力系统不会像我们期望的那样可以知道它的首端电压和首端功率或者是末端电压和末端功率，而是只知道它的首端电压和末端功率，更是使计算变的头疼万分。为了使计算变的简单，我们就可以利用计算机，用C语言编程来实现牛顿-拉夫逊（Newton-Raphson）迭代法，最终实现对电力系统潮流的计算。
用牛顿-拉夫逊迭代法进行电力系统潮流计算的相关概念
1.节点导纳矩阵
如图所示的电力网络,将节点i和j的电压用i和j表示,它们之间的支路导纳表示为yij,那么有基尔霍夫电流定律可知注入接点I的电流i(设流入节点的电流为正)等于离开节点I的电流之和,因此有
j
I
(1-1)
(1-2)
如令 则可将(1-2)改写为:
I=1,2,…,n. (1-3)
上式也可以写为: I=YU (1-4)
其中Y为节点导纳矩阵，也称为稀疏的对称矩阵，它是n×n阶方阵。对角元Yii称为自导纳，它等于与该节点I直接相连的所有支路导纳总和；非对角元Yij（i≠j）称为互导纳或转移导纳，它等于连结节点I，j支路导纳的负数，且有Yij=Yji，当节点I，j之间没有支路直接相连时，Yij=Yji=0。
电力系统的分析计算中，往往要作不同的运行方式下的潮流计算，如果系统发生变化，如投切一条线路或一台变压器，由于改变了一条支路的状态或参数只影响该支路两端节点的自导纳和他们之间的互导纳，因而对每一种运行方式不必重新形成导纳矩阵，只需对原有导纳矩阵作相应的修改即可。
2．潮流计算的功率方程
在实际的电力系统中，已知的条件往往不是节点的注入电流而是负荷和发电机功率，且这些功率一般不随节点的电压变化而变化，而节点的电流则是随电压的变化而变化的，因此在已知节点导纳矩阵的情况下，必须用已知的节点功率来替代未知的节点注入电流，才能求出节点电压，每一个节点的注入功率方程式为：
I
（1-5）
（1-6）
(1-7)
节点注入电流用功率和电压表示为：
（1-8）
功率方程可以表示为：
（1-9）
3．节点分类
对于有n个节点的电力网络，可以列出n个功率方程，由图可知一个节点有四个变量：注入有功功率Pi，注入无功功率Qi，节点电压幅值Ui和相角。n个节点有4n个变量,但只有2n个关系式,所以为了使潮流有确定解,必须给定其中2n个变量。根据给定节点变量的不同，可以有以下三种类型的节点：
PQ节点：给定注入功率Pi，Qi，即已知PGi，PLi，QGi，QLi，待求Ui，δi。例如：降压变电所母线（负荷节点），固定出力的发电厂母线。
PV节点：给定了注入有功功率Pi（PGi，PLi），Ui和QLi，待求QGi（Qi），δi。例如：有一定无功电源的降压变电所母线，有一定储备的发电厂母线。
平衡节点：给定了Ui，δi和PLi ，QLi，待求PGi ，QGi，即Pi，Qi，用来平衡全电网的功率，通常在一个独立的电力系统中只设一个平衡节点。
4．牛顿-拉夫逊迭代法
牛顿-拉夫逊迭代法将解非线性方程组的过程转化为反复求与之相对应的线性方程的求解过程。
对于一个n维非线性方程组： n=1,2,3,…,n
假定其初值为x1(0),x2(0),…,xn(0),也即其近似解,它与真值之间的误差为也即各变量与真解之间的修正量。
将这n个方程式都在初值的附近展开成Taylor级数且忽略二次项及高次项,则可得修正方程
, I=1,2,…n. (1-10)
将修正方程写成矩阵形式:
(1-11)
其中令J= ,称之为雅可比(Jacobi)方阵。
它的第I行，第j列交点的元素为第I个函数对第j个变量xj的偏导数在点（x1(0),x2(0),…,xn(0)）的值，所以方程组是线性方程，可用于求出，从而得到新的近似解，
（1-12） 于是得到一般迭代式：
（1-13）
于是得到近似解：
（1-14）
迭代一直进行到Max{|yi-fi(x1(0),x2(0),…,xn(0))|}<ε或Max{|Δxi(k)|}<ε为止。
用牛顿-拉夫逊迭代法进行潮流计算
设网络中除参考节点外有n个节点，其中1个平衡节点（并令第n个节点为平衡节点），m个PQ节点（为第1~m个节点），有n-m-1个PV 节点。运用牛顿-拉夫逊法直接求解功率方程（1-9），并将Yij=Gij+jBij及代入得：
将实部和虚部分开得：
(2-1)
(2-2)
此外由于系统中还有PV节点，所以还应补充一组方程：
(2-3)
在式（2-3）中，ei、fi分别为迭代过程中求得的节点电压的实部和虚部，Pi为PQ节点和PV节点的注入有功功率，Qi为PQ节点的注入无功功率，Ui为PV节点的电压大小。由（2-1），（2-2），（2-3）三式所组成的方程组一共有2（n-1）个独立方程，其中，式（2-1）类型的有(n-1)个,包括除平衡节点以外所有节点的有功功率Pi的表达式;式(3-2)类型的有(m-1)个,包括所有PQ节点无功功率Qi的表达式,式(2-3)类型的有(n-m)个,包括所有PV节点电压Ui2的表达式.平衡节点s的功率和电压之所以不包括在这方程组内,是由于平衡节点的注入功率不可能事先给定,从而不可能列出Ps,Qs的表达式,而平衡节点的电压则不必求取.
于是建立修正方程式如下:
(2-4)
式中的
(2-5)
(2-6) (2-7)
当j≠I时雅可比方阵的各个元素分别为：
； ；
； ；
当j=I时雅可比方阵的各个元素分别为：
; ;
; ;
;
其中：
用C语言编程计算潮流的流程图
否
是
否
是
用编程方法求解实际问题
如图所示的一个电力网络，
1 3 4
1.25-j3.75 10.00-j30.00
1.667-j5.000
1.667-j5.00
2 5
如何来求解它的潮流分布呢？
已知： 为定值，其余四个节点都是PQ节点，且给定的注入功率分别为：
程序清单如下:
#include
#include
float divRe(b1,b2,b3,b4)
float b1,b2,b3,b4;
{
float a1r;
a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4);
return(a1r);
}
float divIm(b1,b2,b3,b4)
float b1,b2,b3,b4;
{
float a1i;
a1i=(b2*b3-b1*b4)/(b3*b3+b4*b4);
return(a1i);
}
float mulRe(b1,b2,b3,b4)
float b1,b2,b3,b4;
{
float a2r;
a2r=b1*b3-b2*b4;
return(a2r);
}
float mulIm(b1,b2,b3,b4)
float b1,b2,b3,b4;
{
float a2i;
a2i=b2*b3+b1*b4;
return(a2i);
}
float Max(float a[],int n)
{int i;
float max;
for(i=0;i
if(a[i]>a[i+1])
{max=a[i];a[i]=a[i+1];a[i+1]=max;}
return(max);
}
main()
{
int i,j,k,n,km;
float eps,sumpi1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,I1r,I1i;
float pi0[5],qi0[5],detpi[5],detqi[5],Iir0[5],Iii0[5],J0[8][8],detsi[8],detui[8],
u[8][8],l[8][8],y[8],ui1[8],H[4][4],N[4][4],J[4][4],L[4][4],ei1[5],fi1[5];
static float ybr[5][5]={{6.250,-5.000,-1.250,0,0},{-5.000,10.834,-1.667,-1.667,-2.500},
{-1.250,-1.667,12.917,-10.000,0},{0,-1.667,-10.000,12.917,-1.250},
{0,-2.500,0,-1.250,3.750}};
static float ybi[5][5]={{-18.750,15.000,3.750,0,0},{15.000,-32.500,5.000,5.000,7.500},
{3.750,5.000,-38.750,30.000,0},{0,5.000,30.000,-38.750,3.750},
{0,7.500,0,3.750,-11.250}};
float ei0[5]={1.06,1.0,1.0,1.0,1.0};
float fi0[5]={0,0,0,0,0};
float pi[5]={0,0.2,-0.45,-0.4,-0.6};
float qi[5]={0,0.2,-0.15,-0.05,-0.1};
k=0;
km=6;
eps=0.00001;
do{
k+=1;
printf("Now start...\n");
printf("The %d times\n",k);
for(i=1;i<5;i++)
printf("pi[%d]=%f",i,pi[i]);
sumpi2=0;
sumqi2=0;
for(i=1;i<5;i++)
{for(j=0;j<5;j++)
{
sumpi1=(ei0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])+fi0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j]));
sumpi2+=sumpi1;
}
pi0[i]=sumpi2;
printf("pi0[%d]=%f\t",i,pi0[i]);
sumpi2=0;
}
for(i=1;i<5;i++)
{for(j=0;j<5;j++)
{
sumqi1=(fi0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])-ei0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j]));
sumqi2+=sumqi1;
}
qi0[i]=sumqi2;
printf("qi0[%d]=%f\t",i,qi0[i]);
sumqi2=0;
}
for(i=1;i<5;i++)
{
detpi[i]=pi[i]-pi0[i];
detqi[i]=qi[i]-qi0[i];
printf("detpi[%d]=%f\t",i,detpi[i]);
printf("detqi[%d]=%f\t",i,detqi[i]);
}
for(i=1;i<5;i++)
{Iir0[i]=divRe(pi0[i],-qi0[i],ei0[i],-fi0[i]);
Iii0[i]=divIm(pi0[i],-qi0[i],ei0[i],-fi0[i]);
printf("Iir0[%d]=%f\t",i,Iir0[i]);
printf("Iii0[%d]=%f\t",i,Iii0[i]);
}
for(i=0;i<4;i++)
{for(j=0;j<4;j++)
if(i==j) {
H[i][j]=-ybi[i+1][j+1]*ei0[i+1]+ybr[i+1][j+1]*fi0[i+1]+Iii0[i+1];
N[i][j]=ybr[i+1][j+1]*ei0[i+1]+ybi[i+1][j+1]*fi0[i+1]+Iir0[i+1];
J[i][j]=-ybr[i+1][j+1]*ei0[i+1]-ybi[i+1][i+1]*fi0[i+1]+Iir0[i+1];
L[i][j]=-ybi[i+1][j+1]*ei0[i+1]+ybr[i+1][j+1]*fi0[i+1]-Iii0[i+1];
}
else {
H[i][j]=ybr[i+1][j+1]*fi0[i+1]-ybi[i+1][j+1]*ei0[i+1];
N[i][j]=ybr[i+1][j+1]*ei0[i+1]+ybi[i+1][j+1]*fi0[i+1];
J[i][j]=-ybi[i+1][j+1]*fi0[i+1]-ybr[i+1][j+1]*ei0[i+1];
L[i][j]=ybr[i+1][j+1]*fi0[i+1]-ybi[i+1][j+1]*ei0[i+1];
}
}
for(i=0;i<8;i++)
for(j=0;j<8;j++){
if(i%2==0&&j%2==0) J0[i][j]=H[i/2][j/2];
else if(i%2==0&&j%2!=0) J0[i][j]=N[i/2][(j-1)/2];
else if(i%2!=0&&j%2==0) J0[i][j]=J[(i-1)/2][j/2];
else J0[i][j]=L[i/2][(j-1)/2];
}
for(i=0;i<8;i++)
for(j=0;j<8;j++)
printf("J0[%d][%d]=%f\t",i,j,J0[i][j]);
for(i=0;i<8;i++)
{if(i%2==0) detsi[i]=detpi[(i+2)/2];
else detsi[i]=detqi[(i+1)/2];
printf("detsi[%d]=%f\t",i,detsi[i]);
}
for(i=0;i<8;i++) u[i][i]=1.000;
for(n=0;n<8;n++)
{for(i=n;i<8;i++)
{l[i][n]=J0[i][n];
for(j=0;j<=n-1;j++)
l[i][n]-=(l[i][j]*u[j][n]);
}
for(j=n+1;j<8;j++)
{u[n][j]=J0[n][j];
for(i=0;i<=n-1;i++)
u[n][j]-=(l[n][i]*u[i][j]);
u[n][j]/=l[n][n];
}
}
for(i=0;i<8;i++)
{y[i]=detsi[i];
for(j=0;j<=i-1;j++)
y[i]-=(l[i][j]*y[j]);
y[i]/=l[i][i];
}
for(i=7;i>=0;i--)
{detui[i]=y[i];
for(j=i+1;j
detui[i]-=(u[i][j]*detui[j]);
}
for(i=0;i<8;i++)
printf("detui[%d]=%f\t",i,detui[i]);
for(i=0;i<8;i++)
{
{if(i%2==0) ui1[i]=detui[i]+fi0[i/2+1];
else ui1[i]=detui[i]+ei0[(i+1)/2];
}
printf("ui1[%d]=%f\t",i,ui1[i]);
}
for(i=1;i<5;i++)
{ei1[i]=ui1[2*i-1];
fi1[i]=ui1[2*i-2];
}
for(i=1;i<5;i++)
{printf("ei1[%d]=%f\t",i,ei1[i]);
printf("fi1[%d]=%f\t",i,fi1[i]);
}
max=Max(detui,8);
printf("max=%f",max);
for(i=1;i<5;i++)
{ei0[i]=ei1[i];
fi0[i]=fi1[i];
}
for(i=1;i<5;i++)
{pi[i]=detpi[i]+pi0[i];
qi[i]=detqi[i]+qi0[i];
}
}while(max>eps&&k
printf("All do %d times\n",k);
sumir=0;
sumii=0;
for(i=0;i<5;i++)
{
I1r=mulRe(ybr[0][i],-ybi[0][i],ei0[i],-fi0[i]);
I1i=mulIm(ybr[0][i],-ybi[0][i],ei0[i],-fi0[i]);
sumir+=I1r;
sumii+=I1i;
}
pi[0]=mulRe(ei0[0],fi0[0],sumir,sumii);
qi[0]=mulIm(ei0[0],fi0[0],sumir,sumii);
printf("S1=%f+j%f\n",pi[0],qi[0]);
ei1[0]=ei0[0];
fi1[0]=fi0[0];
for(i=0;i<5;i++)
printf("u%d=%f<%f\n",i+1,sqrt(ei1[i]*ei1[i]+fi1[i]*fi1[i]),atan(fi1[i]/ei1[i])*180/3.14159);
}
最终求出平衡节点的功率和线路的功率以及各个节点的电压为：
利用这中方法进行潮流的计算，其收敛的速度是显而易见的，如上面这到题，在经过三次迭代以后已经达到我们要求的水平，但必须建立在有一个比较优的初值的基础上。
用C语言进行复杂电力系统的潮流计算.doc
