电力系统潮流计算C程序.docx
- 文档编号:29622217
- 上传时间:2023-07-25
- 格式:DOCX
- 页数:24
- 大小:19.36KB
电力系统潮流计算C程序.docx
《电力系统潮流计算C程序.docx》由会员分享,可在线阅读,更多相关《电力系统潮流计算C程序.docx(24页珍藏版)》请在冰豆网上搜索。
电力系统潮流计算C程序
电力系统潮流计算
注:
这是一个基于N-R法的潮流计算通用程序,仅提供了子程序,需要做些处理才能成为一个可运行的计算程序!
此程序非我原创,仅与大家共享!
!
!
/*******************************************************************
* 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是*
* Newton_Raphson法. *
* 程序中所用的变量说明如下:
*
* N:
网络节点总数. M:
网络的PQ节点数. *
* L:
网络的支路总数. N0:
雅可比矩阵的行数. *
* N1:
N0+1 K:
打印开关.K=1,则打印;否则,不打印.*
* K1:
子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则*
* 为直角坐标形式. *
* D:
有功及无功功率误差的最大值. *
* G(I,J):
Ybus的电导元素(实部). *
* B(I,J):
Ybus的电纳元素(虚部). *
* G1(I):
第I支路的串联电导. B1(I):
第I支路的串联电纳. *
* C1(I):
第I支路的pie型对称接地电纳. *
* C(I,J):
第I节点J支路不对称接地电纳. *
* CO(I):
第I节点的接地电纳. *
* S1(I):
第I节点的起始节点号. E1(I):
第I节点的终止节点号.*
* P(I) :
第I节点的注入有功功率. Q(I):
第I节点的注入无功功率.*
* P0(I):
第I节点有功功率误差. Q0(I):
第I节点无功功率误差.*
* V0(I):
第I节点(PV节点)的电压误差(平方误差). *
* V(I) :
第I节点的电压误差幅值. *
* E(I) :
第I节点的电压的实部. F(I):
第I节点的电压的虚部. *
* JM(I,J):
Jacoby矩阵的第I行J列元素. *
* A(I,J):
修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结*
* 束后A矩阵的最后一列存放修正的解. *
* P1(I):
第I支路由S1(I)节点注入的有功功率. *
* Q1(I):
第I支路由S1(I)节点注入的无功功率. *
* P2(I):
第I支路由E1(I)节点注入的有功功率. *
* Q2(I):
第I支路由E1(I)节点注入的无功功率. *
* P3(I):
第I支路的有功功率损耗. *
* Q3(I):
第I支路的无功功率损耗. *
* ANGLE(I):
第I节点电压的角度. *
*******************************************************************/
#include
#include
#definef1(i)(i-1)
/*把习惯的一阶矩阵的下标转化为C语言数组下标*/
#definef2(i,j,n)((i-1)*(n)+j-1)
/*把习惯的二阶矩阵的下标转化为C语言数组下标*/
/****************************************************
* 本子程序根据所给的支路导纳及有关信息,形成结点*
*导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B *
****************************************************/
voidybus(intn,intl,intm,float*g,float*b,float*g1,float*b1,float*c1,\
float*c,float*co,intk,int*s1,int*e1)
{
externFILE*file4;
FILE*fp;
inti,j,io,i0;
intpos1,pos2;
intst,en;
if(file4==NULL)
{
fp=stdout;
}
else
{
fp=file4;/*输出到文件*/
}
/*初始化矩阵G,B*/
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
{
pos2=f2(i,j,n);
g[pos2]=0;b[pos2]=0;
}
}
/*计算支路导纳*/
for(i=1;i<=l;i++)
{
/*计算对角元*/
pos1=f1(i);
st=s1[pos1];en=e1[pos1];
pos2=f2(st,st,n);
g[pos2]+=g1[pos1];
b[pos2]+=b1[pos1]+c1[pos1];
pos2=f2(en,en,n);
g[pos2]+=g1[pos1];
b[pos2]+=b1[pos1]+c1[pos1];
/*计算非对角元*/
pos2=f2(st,en,n);
g[pos2]-=g1[pos1];
b[pos2]-=b1[pos1];
g[f2(en,st,n)]=g[f2(st,en,n)];
b[f2(en,st,n)]=b[f2(st,en,n)];
}
/*计算接地支路导纳*/
for(i=1;i<=n;i++)
{
/* 对称部分 */
b[f2(i,i,n)]+=co[f1(i)];
/*非对称部分 */
for(j=1;j<=l;j++)
{
b[f2(i,i,n)]+=c[f2(i,j,l)];
}
}
if(k!
=1)
{
return;/*如果K不为1,则返回;否则,打印导纳矩阵*/
}
fprintf(fp,"\n BUSADMITTANCEMATRIXY(BUS):
");
fprintf(fp,"\n*******************ARRAYG********************");
for(io=1;io<=n;io+=5)
{
i0=(io+4)>n?
n:
(io+4);
fprintf(fp,"\n");
for(j=io;j<=i0;j++)
{
fprintf(fp,"%13d",j);
}
for(i=1;i<=n;i++)
{
fprintf(fp,"\n%2d",i);
for(j=io;j<=i0;j++)
{
fprintf(fp,"%13.6f",g[f2(i,j,n)]);
}
}
fprintf(fp,"\n");
}
fprintf(fp,"\n*******************ARRAYB********************");
for(io=1;io<=n;io+=5)
{
i0=(io+4)>n?
n:
(io+4);
fprintf(fp,"\n");
for(j=io;j<=i0;j++)
{
fprintf(fp,"%13d",j);
}
for(i=1;i<=n;i++)
{
fprintf(fp,"\n%2d",i);
for(j=io;j<=i0;j++)
{
fprintf(fp,"%13.6f",b[f2(i,j,n)]);
}
}
fprintf(fp,"\n");
}
fprintf(fp,"\n************************************************");
}
/*******************************************
* 本子程序根据所给的功率及电压等数据 *
*求出功率及电压误差量,并返回最大有功功率*
*以用于与给定误差比较.如打印参数K=1,则输*
*出P0,Q0(对PQ结点),V0(对PV结点). *
*******************************************/
voiddpqc(float*p,float*q,float*p0,float*q0,float*v,float*v0,intm,\
intn,float*e,float*f,intk,float*g,float*b,float*dd)
{
externFILE*file4;
FILE*fp;
inti,j,l;
intpos1,pos2;
floata1,a2,d1,d;
if(file4==NULL)
{
fp=stdout;/*输出到屏幕*/
}
else
{
fp=file4; /*输出到文件 */
}
l=n-1;
if(k==1)
{
fprintf(fp,"\n CHANGEOFP0,V**2,P0(I),Q0(I),V0(I)");
fprintf(fp,"\n I P0(I) Q0(I)");
}
for(i=1;i<=l;i++)
{
a1=0;a2=0;
pos1=f1(i);
for(j=1;j<=n;j++)
{
/*a1,a2对应课本p185式(4-86)中括号内的式子*/
pos2=f2(i,j,n);
a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];
a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];
}
/*计算式(4-86)(4-87)中的deltaPi */
p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;
if(i<=m)
{/*计算PQ结点中的deltaQi */
q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;
}
else
{/*计算PV结点中的deltaVi平方 */
v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];
}
/*输出结果*/
if(k==1)
{
if(i { fprintf(fp,"\n%2d%15.6e%15.6e",i,p0[pos1],q0[pos1]); } elseif(i==m) { fprintf(fp,"\n%2d%15.6e%15.6e",i,p0[pos1],q0[pos1]); fprintf(fp,"\n I P0(I) V0(I)"); } else { fprintf(fp,"\n%2d%15.6e%15.6e",i,p0[pos1],v0[pos1]); } } } /*找到deltaP和deltaQ中的最小者,作为收敛指标,存在dd中*/ d=0; for(i=1;i<=l;i++) { pos1=f1(i); d1=p0[pos1]>0? p0[pos1]: -p0[pos1]; if(d { d=d1; } if(i<=m) { d1=q0[pos1]>0? q0[pos1]: -q0[pos1]; if(d { d=d1; } } } (*dd)=d; } /*************************************************** * 本子程序根据节点导纳及电压求Jacoby矩阵,用于求* * 电压修正量,如打印参数K=1,则输出Jacoby矩阵. * * 对应于课本P186式(4-89)(4-90) * ***************************************************/ voidjmcc(intm,intn,intn0,float*e,float*f,float*g,float*b,float*jm,intk) { externFILE*file4; FILE*fp; inti,j,i1,io,i0,ns; intpos1,pos2; if(file4==NULL) { fp=stdout; } else { fp=file4; } /*初始化矩阵jm*/ for(i=1;i<=n0;i++) { for(j=1;j<=n0;j++) { jm[f2(i,j,n0)]=0; } } ns=n-1;/*去掉一个平衡结点*/ /*计算式(4-89)(4-90)*/ for(i=1;i<=ns;i++) { /*计算式(4-90)*/ for(i1=1;i1<=n;i1++) { /*pos1是式(4-90)中的j*/ pos1=f1(i1); /*pos2是式(4-90)中的ij*/ pos2=f2(i,i1,n); if(i<=m)/*i是PQ结点*/ { /*计算式(4-90)中的Jii等式右侧第一部分*/ jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1]; /*计算式(4-90)中的Lii等式右侧第一部分*/ jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; } /*计算式(4-90)中的Hii等式右侧第一部分*/ jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; /*计算式(4-90)中的Nii等式右侧第一部分*/ jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1]; } /*pos2是式(4-90)中的ii*/ pos2=f2(i,i,n); /*pos1是式(4-90)中的i*/ pos1=f1(i); if(i<=m)/*i是PQ结点*/ { /*计算式(4-90)中的Jii*/ jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1]; /*计算式(4-90)中的Lii*/ jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[pos1]; } /*计算式(4-90)中的Hii*/ jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1]; /*计算式(4-90)中的Jii*/ jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1]; if(i>m)/*PV结点*/ { /*计算式(4-90)中的Rii*/ jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1]; /*计算式(4-90)中的Sii*/ jm[f2(2*i-1,2*i,n0)]=-2*f[pos1]; } /*计算式(4-89)*/ for(j=1;j<=ns;j++) { if(j! =i) { /*pos1是式(4-89)中的i*/ pos1=f1(i); /*pos2是式(4-89)中的ij*/ pos2=f2(i,j,n); /*计算式(4-89)中的Nij*/ jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1]; /*计算式(4-89)中的Hij*/ jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1]; if(i<=m)/*i是PQ结点*/ { /*计算式(4-89)中的Lij(=-Hij)*/ jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)]; /*计算式(4-89)中的Jij(=Nij)*/ jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)]; } else/*i是PV结点*/ { /*计算式(4-89)中的Rij(=0)*/ jm[f2(2*i-1,2*j-1,n0)]=0; /*计算式(4-89)中的Sij(=0)*/ jm[f2(2*i-1,2*j,n0)]=0; } } } } if(k! =1) { return; } /*输出Jacoby矩阵*/ fprintf(fp,"\nJ MATRIX(C)"); for(io=1;io<=n0;io+=5) { i1=(io+4)>n0? n0: (io+4); fprintf(fp,"\n"); for(j=io;j<=i1;j++) { fprintf(fp,"%10d",j); } for(i=1;i<=n0;i++) { fprintf(fp,"\n%2d",i); for(j=io;j<=i1;j++) { fprintf(fp,"%12.6f",jm[f2(i,j,n0)]); } } } fprintf(fp,"\
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 电力系统 潮流 计算 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)