电力系统设计潮流计算C程序.docx
《电力系统设计潮流计算C程序.docx》由会员分享，可在线阅读，更多相关《电力系统设计潮流计算C程序.docx（13页珍藏版）》请在冰豆网上搜索。
电力系统设计潮流计算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
升级会员 