电力系统潮流计算完整c语言程序(含网损计算的最终版)
// flow.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include
#include
#include "NEquation.h"
#include "math.h"
#include "config.h"
void test()
{
NEquation ob1;
ob1.SetSize(2);
ob1.Data(0,0)=1;
ob1.Data(0,1)=2;
ob1.Data(1,0)=2;
ob1.Data(1,1)=1;
ob1.Value(0)=5;
ob1.Value(1)=4;
ob1.Run();
printf("x1=%f\n",ob1.Value(0));
printf("x2=%f\n",ob1.Value(1));
}
void GetData()//Read the data
{
FILE *fp;
int i;
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\data.txt","r"))==NULL)
{
printf("Can not open the file named 'data.txt' \n");
return;
}
for(i=0;i<=Bus_Num-1;i++)
{
fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].Voltage,&gBus[i].Phase,
&gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type);
}
for(i=0;i<=Line_Num-1;i++)
{
fscanf(fp,"%d,%d,%d,%f,%f,%f,%f",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J,
&gLine[i].R,&gLine[i].X,&gLine[i].B,&gLine[i].k);
}
fclose(fp);
}
void GetYMatrix()
{
int i,j;
FILE *fp;
// calculate the Y matrix
int l,bal_num;
float r,x,d1,g,b;
for(i=0;i
if(gBus[i].Type==3)
bal_num=i+1;
}
for(i=0;i
{
if(gLine[i].No_I==bal_num)
gLine[i].No_I=1;
else if(gLine[i].No_I==1)
gLine[i].No_I=bal_num;
}
{
if(gLine[i].No_J==bal_num)
gLine[i].No_J=1;
else if(gLine[i].No_J==1)
gLine[i].No_J=bal_num;
}
}
for(i=0;i
gY_G[i][j]=0.0;
gY_B[i][j]=0.0;
}
for(l=0;l
i=gLine[l].No_I-1;
j=gLine[l].No_J-1;
r=gLine[l].R;
x=gLine[l].X;
d1=r*r+x*x;
g=r/d1;
b=-x/d1;
if(gLine[l].k==0)
{
gY_G[i][i]=gY_G[i][i]+g;
gY_G[j][j]=gY_G[j][j]+g;
gY_B[i][i]=gY_B[i][i]+b+gLine[l].B;
gY_B[j][j]=gY_B[j][j]+b+gLine[l].B;
gY_G[i][j]=gY_G[i][j]-g;
gY_G[j][i]=gY_G[j][i]-g;
gY_B[i][j]=gY_B[i][j]-b;
gY_B[j][i]=gY_B[j][i]-b;
}
else
{
gY_G[i][i]=gY_G[i][i]+g/gLine[l].k+g*(1-gLine[l].k)/gLine[l].k/gLine[l].k;
gY_G[j][j]=gY_G[j][j]+g/gLine[l].k+g*(gLine[l].k-1)/gLine[l].k;
gY_B[i][i]=gY_B[i][i]+b/gLine[l].k+b*(1-gLine[l].k)/gLine[l].k/gLine[l].k+gLine[l].B;
gY_B[j][j]=gY_B[j][j]+b/gLine[l].k+b*(gLine[l].k-1)/gLine[l].k+gLine[l].B;
gY_G[i][j]=gY_G[i][j]-g;
gY_G[j][i]=gY_G[j][i]-g;
gY_B[i][j]=gY_B[i][j]-b;
gY_B[j][i]=gY_B[j][i]-b;
}
}
// output the Y matrix
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\ymatrix.txt","w"))==NULL)
{
printf("Can not open the file named 'ymatrix.txt' \n");
return ;
}
fprintf(fp,"---Y Matrix---\n");
for(i=0;i<=Bus_Num-1;i++)
{
for(j=0;j<=Bus_Num-1;j++)
{
fprintf(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);
}
}
fclose(fp);
}
void SetInitial()
{
int i;
for(i=0;i<=Bus_Num-1;i++)
{
gf[i]=gBus[i].Volt
age*sin(gBus[i].Phase);
ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);
}
}
void GetUnbalance()
{
int i,j,n;
float a;
FILE *fp;
for(i=0;i
if(gBus[i].Type==3)
{
a=gBus[i].Voltage;
gBus[i].Voltage=gBus[0].Voltage;
gBus[0].Voltage=a;
a=gBus[i].Phase;
gBus[i].Phase=gBus[0].Phase;
gBus[0].Phase=a;
a=gBus[i].GenP;
gBus[i].GenP=gBus[0].GenP;
gBus[0].GenP=a;
a=gBus[i].GenQ;
gBus[i].GenQ=gBus[0].GenQ;
gBus[0].GenQ=a;
a=gBus[i].LoadP;
gBus[i].LoadP=gBus[0].LoadP;
gBus[0].LoadP=a;
a=gBus[i].LoadQ;
gBus[i].LoadQ=gBus[0].LoadQ;
gBus[0].LoadQ=a;
a=gBus[i].Type;
gBus[i].Type=gBus[0].Type;
gBus[0].Type=a;
}
}
for(i=0,n=1;i
gDelta_P[i]=gBus[n].GenP-gBus[n].LoadP;
if(gBus[n].Type==1)
gDelta_Q[i]=gBus[n].GenQ-gBus[n].LoadQ;
else
gDelta_Q[i]=gBus[n].Voltage*gBus[n].Voltage-(ge[n]*ge[n]+gf[n]*gf[n]);
for(j=0;j
gDelta_P[i]=gDelta_P[i]-ge[n]*(gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j])-gf[n]*(gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j]);
if(gBus[n].Type==1)
gDelta_Q[i]=gDelta_Q[i]-gf[n]*(gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j])+ge[n]*(gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j]);
}
}
for(i=0;i
gDelta_PQ[2*i]=gDelta_P[i];
gDelta_PQ[2*i+1]=gDelta_Q[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\unbalance.txt","w"))==NULL)
{
printf("Can not open the file named 'unbalance.txt' \n");
return ;
}
fprintf(fp,"---Unbalance---\n");
for(i=0;i<=2*Bus_Num-3;i++)
{
fprintf(fp,"Unbalance[%d]=%10.5f\n",i+1,gDelta_PQ[i]);
}
fclose(fp);
}
void GetJaccobi()
{
FILE *fp;
int i,j,n;
float ia[Bus_Num-1],ib[Bus_Num-1];
for(i=0,n=1;i
ia[i]=0;
ib[i]=0;
for(j=0;j
ia[i]=ia[i]+gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j];
ib[i]=ib[i]+gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j];
}
}
for(i=0,n=1;i
for(j=0;j
if(i!=j)
{
gJaccobi[2*i][2*j]=-gY_B[n][j+1]*ge[n]+gY_G[n][j+1]*gf[n];
gJaccobi[2*i][2*j+1]=gY_G[n][j+1]*ge[n]+gY_B[n][j+1]*gf[n];
if(gBus[n].Type==2)
{
gJaccobi[2*i+1][2*j]=0;
gJaccobi[2*i+1][2*j+1]=0;
}
else
{
gJaccobi[2*i+1][2*j]=-gJaccobi[2*i][2*j+1];
gJaccobi[2*i+1][2*j+1]=gJaccobi[2*i][2*j];
}
}
else
{
gJaccobi[2*i][2*j]=-gY_B[n][j+1]*ge[n]+gY_G[n][j+1]*gf[n]+ib[i];
gJaccobi[2*i][2*j+1]=gY_G[n][j+1]*ge[n]+gY_B[n][j+1]*gf[n]+ia[i];
if(gBus[n].Type==2)
{
gJaccobi[2*i+1][2*j]=2*gf[n];
gJaccobi[2*i+1][2*j+1]=2*ge[n];
}
else
{
gJaccobi[2*i+1][2*j]=-gY_G[n][j+1]*ge[n]-gY_B[n][j+1]*gf[n]+ia[i];
gJaccobi[2*i+1][2*j+1]=-gY_B[n][j+1]*ge[n]+gY_G[n][j+1]*gf[n]-ib[i];
}
}
}
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\jaccobi.txt","w"))==NULL)
{
printf("Can not open the file named 'jaccobi
.txt' \n");
return ;
}
fprintf(fp,"---Jaccobi Matrix---\n");
for(i=0;i<=2*Bus_Num-3;i++)
{
for(j=0;j<=2*Bus_Num-3;j++)
{
fprintf(fp,"jaccobi(%d,%d)=%10.5f\n",i+1,j+1,gJaccobi[i][j]);
}
}
fclose(fp);
}
void GetRevised()
{
int i,j;
FILE *fp;
NEquation ob1;
ob1.SetSize(2*(Bus_Num-1));
for(i=0;i<2*(Bus_Num-1);i++)
for(j=0;j<2*(Bus_Num-1);j++)
ob1.Data(i,j)=gJaccobi[i][j];
for(i=0;i<2*(Bus_Num-1);i++)
ob1.Value(i)=gDelta_PQ[i];
ob1.Run();
for(i=0;i
gDelta_f[i]=ob1.Value(2*i);
gDelta_e[i]=ob1.Value(2*i+1);
gDelta_fe[2*i]=gDelta_f[i];
gDelta_fe[2*i+1]=gDelta_e[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\revised.txt","w"))==NULL)
{
printf("Can not open the file named 'revised.txt' \n");
return;
}
fprintf(fp,"---Revised---\n");
for(i=0;i<=2*Bus_Num-3;i++)
{
fprintf(fp,"revised[%d]=%10.5f\n",i+1,gDelta_fe[i]);
}
fclose(fp);
}
void GetNewValue()
{
int i;
FILE *fp;
for(i=0;i
gf[i+1]=gf[i+1]+gDelta_f[i];
ge[i+1]=ge[i+1]+gDelta_e[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\newvalue.txt","w"))==NULL)
{
printf("Can not open the file named 'newvalue.txt' \n");
return ;
}
fprintf(fp,"---New Value---\n");
for(i=0;i
fprintf(fp,"f(%d)=%10.5f,e(%d)=%10.5f\n",i+1,gf[i],i+1,ge[i]);
}
fclose(fp);
}
int main(int argc, char* argv[])
{
int i,Count_Num;
float maxValue;
test();
GetData();
GetYMatrix();
SetInitial();
for(Count_Num=0;Count_Num<=100;Count_Num++)
{
GetUnbalance();
GetJaccobi();
GetRevised();
GetNewValue();
maxValue=fabs(gDelta_fe[0]);
for(i=1;i<=2*(Bus_Num-1)-1;i++)
{
if(maxValue
maxValue=fabs(gDelta_fe[i]);
}
}
if(maxValue
break;
}
}
printf("%d\n",Count_Num);
for(i=0;i<=Bus_Num-1;i++)
{
printf("%10.5f\n",sqrt(ge[i]*ge[i]+gf[i]*gf[i]));
}
printf("\n***平衡节点功率***\n");
int j;
FILE *fp;
float gPs,gQs;
gPs=0;
gQs=0;
for(j=0;j<=Bus_Num-1;j++)
{ gPs=gPs+ge[0]*(gY_G[0][j]*ge[j]-gY_B[0][j]*gf[j])+gf[0]*(gY_G[0][j]*gf[j]+gY_B[0][j]*ge[j]);
gQs=gQs+gf[0]*(gY_G[0][j]*ge[j]-gY_B[0][j]*gf[j])-ge[0]*(gY_G[0][j]*gf[j]+gY_B[0][j]*ge[j]);
}
printf("平衡节点功率为%10.5f+j%10.5f\n",gPs,gQs);
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\网络损耗.txt","w"))==NULL)
{
printf("Can not open the file named '网络损耗.txt' \n"); //由于是在主程序,故下一行不用用到return
}
fprintf(fp,"***平衡节点功率***\n");
fprintf(fp,"平衡节点功率为%10.5f+j%10.5f\n",gPs,gQs);
//计算PV节点的无功功率
printf("\n***PV节点的注入无功***\n");
fprintf(fp,"\n***PV节点的注入无功***\n");
float Q_PV;
for(i=0;i
if(gBus[i].Type==2)
{
Q_PV=0;
for(j=0;j
Q_PV+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[j]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
}
printf("第%d个节点的注入无功是：%10.5f\n",i+1,Q_PV);
fprintf(fp,"第%d个节点的注入无功是：%10.5f\n",i+1,Q_PV);
}
}
//计算线路的损耗
printf("\n******线路网损******\n");
fprintf(fp,"\n***线路网损***\n");
int First_plot,End_plot;
float First_B,End_B;
float r,x,d,g,b;
float Branch_G,Branch_B;
float A,B,C,D;
float P_First,Q_First,P_End,Q_End;
float P_Loss,Q_Loss;
for(i=0;i
First_plot=gLine[i].No_I-1; //线路一侧连接的节点号减1
End_plot=gLine[i].No_J-1;
r=gLine[i].R;
x=gLine[i].X;
d=r*r+x*x;
g=r/d;
b=-x/d;
if(gLine[i].k==0) //无变压器的支路,其中部分式子已经省略零项//
{
First_B=gLine[i].B;
End_B=gLine[i].B;
Branch_G=g;
Branch_B=b;
A=-gf[First_plot]*First_B+(ge[First_plot]-ge[End_plot])*Branch_G-Branch_B*(gf[First_plot]-gf[End_plot]);
B=ge[First_plot]*First_B+(ge[First_plot]-ge[End_plot])*Branch_B+Branch_G*(gf[First_plot]-gf[End_plot]);
P_First=ge[First_plot]*A+gf[First_plot]*B;
Q_First=gf[First_plot]*A-ge[First_plot]*B;
C=-gf[End_plot]*End_B+(-ge[First_plot]+ge[End_plot])*Branch_G-Branch_B*(-gf[First_plot]+gf[End_plot]);
D=ge[End_plot]*End_B+(-ge[First_plot]+ge[End_plot])*Branch_B+Branch_G*(-gf[First_plot]+gf[End_plot]);
P_End=ge[End_plot]*C+gf[End_plot]*D;
Q_End=gf[End_plot]*C-ge[End_plot]*D;
P_Loss=P_First+P_End;
Q_Loss=Q_First+Q_End;
printf("线路%d%d上首端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_I,P_First,Q_First);
printf("线路%d%d上末端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_J,P_End,Q_End);
printf("线路%d%d上的功率损耗:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,P_Loss,Q_Loss);
printf("\n");
fprintf(fp,"线路%d%d上首端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_I,P_First,Q_First);
fprintf(fp,"线路%d%d上末端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_J,P_End,Q_End);
fprintf(fp,"线路%d%d上的功率损耗:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,P_Loss,Q_Loss);
fprintf(fp,"\n");
}
else //含有变压器的支路//
{
First_B=b*(1-gLine[i].k)/(gLine[i].k*gLine[i].k);
End_B=b*(gLine[i].k-1)/gLine[i].k;
Branch_B=b/gLine[i].k;
A=-gf[First_plot]*First_B-Branch_B*(gf[First_plot]-gf[End_plot]);
B=ge[First_plot]*First_B+(ge[First_plot]-ge[End_plot])*Branch_B;
P_First=ge[First_plot]*A+gf[First_plot]*B;
Q_First=gf[First_plot]*A-ge[First_plot]*B;
C=-gf[End_plot]*End_B-Branch_B*(-gf[First_plot]+gf[End_plot]);
D=ge[End_plot]*End_B+(-ge[Fi
rst_plot]+ge[End_plot])*Branch_B;
P_End=ge[End_plot]*C+gf[End_plot]*D;
Q_End=gf[End_plot]*C-ge[End_plot]*D;
P_Loss=P_First+P_End;
Q_Loss=Q_First+Q_End;
printf("线路%d%d上首端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_I,P_First,Q_First);
printf("线路%d%d上末端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_J,P_End,Q_End);
printf("线路%d%d上的功率损耗:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,P_Loss,Q_Loss);
printf("\n");
fprintf(fp,"线路%d%d上首端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_I,P_First,Q_First);
fprintf(fp,"线路%d%d上末端%d的线路功率:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,gLine[i].No_J,P_End,Q_End);
fprintf(fp,"线路%d%d上的功率损耗:P=%10.5f,Q=%10.5f\n",gLine[i].No_I,gLine[i].No_J,P_Loss,Q_Loss);
fprintf(fp,"\n");
}
}
while(true)
{
}
return 0;
}