📄 flow.cpp
字号:
// flow.cpp : Defines the entry point for the console application.
//sys. of 4 nodes, with one PV node,designed by maqiyan,guided by professor qizheng
//库函数说明
#include "stdafx.h"
#include <iostream.h>
#include <fstream.h>
#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)=1;
ob1.Data(1,0)=2;
ob1.Data(1,1)=4;
ob1.Value(0)=12;
ob1.Value(1)=12;
ob1.Run();
printf(" the test value is: %5.3f\t%5.3f\n",ob1.Value(0),ob1.Value(1));
}
//获取节点,变压器，并联支路和线路数据
void GetData()
{
FILE *fp;
int i;
if((fp=fopen("D:\\潮流计算080607\\MQY flow\\PVflow\\data\\data2PV.txt","r"))==NULL)
{
printf("Can not open the file named 'data2PV.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",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J,
&gLine[i].R,&gLine[i].X,&gLine[i].B);
}
for(i=0;i<=Bus_Num-1;i++)
{
fscanf(fp,"%d,%d,%f,%f",&gShunt[i].Num,&gShunt[i].NumI,&gShunt[i].G,&gShunt[i].B);
}
for(i=0;i<TR_Num;i++)
{
fscanf(fp,"%d,%d,%d,%f,%f,%f",
&TR1[i].Num,&TR1[i].No_I,&TR1[i].No_J,&TR1[i].k,&TR1[i].R,&TR1[i].X);
}
fclose(fp);
}
//计算节点导纳矩阵
void GetYMatrix()
{
int i,j,bus1,bus2;
float r,x,d,g,b,g1,g2,b1,b2; //g1,g2,b1，b2用于变压器PI型模型对地导纳之路计算
float r1=0,r2=0,x1=0,x2=0;
FILE *fp;
for(i=0;i<=Bus_Num-1;i++)
{
for(j=0;j<=Bus_Num-1;j++)
{
gY_G[i][j]=0;
gY_B[i][j]=0;
}
}
for(i=0; i<=Line_Num-1; i++)
{
bus1=gLine[i].No_I-1;
bus2=gLine[i].No_J-1;
r=gLine[i].R;
x=gLine[i].X;
d=r*r+x*x;
g1=0;
g2=0;
b1=0;
b2=0;
//变压器PI型模型数据处理(节点之间连接用阻抗，对地用导纳)
for(j=0;j<TR_Num;j++)
{
if ((gLine[i].No_I==TR1[j].No_I)&&(gLine[i].No_J==TR1[j].No_J))
{
r/=TR1[j].k;
x/=TR1[j].k;
d=r*r+x*x;
r1=TR1[j].R/(TR1[j].k*TR1[j].k-TR1[j].k);
r2=TR1[j].R/(1-TR1[j].k);
x1=TR1[j].X/(TR1[j].k*TR1[j].k-TR1[j].k);
x2=TR1[j].X/(1-TR1[j].k);
g1=r1/(r1*r1+x1*x1);
g2=r2/(r2*r2+x2*x2);
b1=-x1/(r1*r1+x1*x1);
b2=-x2/(r2*r2+x2*x2);
}
}
if(d>1e-6)
{
g=r/d;
b=-x/d;
}
else
{
g=0;
b=0;
}
gY_G[bus1][bus1]+=g+g1;
gY_G[bus2][bus2]+=g+g2;
gY_G[bus1][bus2]-=g;
gY_G[bus2][bus1]-=g;
gY_B[bus1][bus1]+=b+gLine[i].B/2+b1;
gY_B[bus2][bus2]+=b+gLine[i].B/2+b2;
gY_B[bus1][bus2]-=b;
gY_B[bus2][bus1]-=b;
}
//将对地并联支路导纳写入节点导纳矩阵
for(i=0; i<=Bus_Num-1;i++)
{
for(j=0; j<=Bus_Num-1;j++)
{
if (i==j)
{
gY_G[i][j]+=gShunt[i].G;
gY_B[i][j]+=gShunt[i].B;
}
}
}
// 输出节点导纳距阵
if((fp=fopen("D:\\潮流计算080607\\MQY flow\\PVflow\\data\\ymatrix2.txt","w"))==NULL)
{
printf("Can not open the file named 'ymatrix2.txt' \n");
return ;
}
fprintf(fp,"---Y Matrix---\n");
printf(" the Y matrix is calcalated below \n");
for(i=0;i<=Bus_Num-1;i++)
{
int k;
for(j=0,k=0;j<=Bus_Num-1;j++,k++)
{
fprintf(fp,"Y(%d,%d)=(%7.3f,%7.3f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);
printf("%7.3f,%7.3f\t",gY_G[i][j],gY_B[i][j]);
if ((k+1)%(Bus_Num)) continue;
else printf("\n");
}
}
printf("\n");
printf("\n");
fclose(fp);
}
//设电压初值
void SetInitial()
{
int i;
for(i=0;i<=Bus_Num-1;i++)
{
gf[i]=gBus[i].Voltage*sin(gBus[i].Phase);
ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);
}
}
//求取功率不平衡量
void GetUnbalance()
{
float Bus_P[Bus_Num-1],Bus_Q[Bus_Num-1]; //定义节点有功，无功数组
int i,j,k;
for (i=0;i<Bus_Num-1;i++)
{
Bus_P[i]=0;
Bus_Q[i]=0;
UU[i]=0;
for (j=0;j<Bus_Num-1;j++)
{
if (gBus[j].Type==0&&gBus[i].Type==0)
{
Bus_P[i]+=ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])
+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
Bus_Q[i]+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])
-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
}
if (gBus[j].Type==1&&gBus[i].Type==1)
{
UU[i]=gf[i]*gf[i]+ge[i]*ge[i];
}
}
}
printf(" the unbalanced value of Power and voltage is calcalated below\n");
for (i=0;i<Bus_Num-1;i++)
{
if (gBus[i].Type==0)
{
gDelta_P[i]=gBus[i].GenP-gBus[i].LoadP-Bus_P[i];
gDelta_Q[i]=gBus[i].GenQ-gBus[i].LoadQ-Bus_Q[i];
printf("gDelta_P[%d]=%f gDelta_Q[%d]=%f\n",i+1,gDelta_P[i],i+1,gDelta_Q[i]);
}
if (gBus[i].Type==1)
{
gDelta_UU[i]=UU[i]-(gf[i]*gf[i]+ge[i]*ge[i]);
gDelta_P[i]=gBus[i].GenP-gBus[i].LoadP-Bus_P[i];
printf("gDelta_P[%d]=%f gDelta_UU[%d]=%f\n",i+1,gDelta_P[i],i+1,gDelta_UU[i]);
}
printf("\n");
printf("\n");
}
for (i=0,k=0;i<Bus_Num-1;i++)
{
k=i;
if (gBus[i].Type==0)
{
gDelta_PQ[2*k]=gDelta_P[i];
gDelta_PQ[2*k+1]=gDelta_Q[i];
}
if (gBus[i].Type==1)
{
gDelta_PQ[2*k]=gDelta_P[i];
gDelta_PQ[2*k+1]=gDelta_UU[i];
}
}
printf(" the gDelta_PQ Matrix is calculated below:\n");
for (i=0,k=0;i<2*(Bus_Num-1);i++)
{
printf("gDelta_PQ[%d]=%7.3f\n",i+1,gDelta_PQ[i]);
}
printf("\n");
printf("\n");
}
//计算雅可比矩阵用自定义函数
//自定义函数1
float add1 (int i,int j)
{
float add=0;
for (j=0;j<=Bus_Num-1;j++)
{
if(j!=i)
{
add+=gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j];
}
}
return add;
}
//自定义函数2
float add2 (int i,int j)
{
float add=0;
for (j=0;j<=Bus_Num-1;j++)
{
if(j!=i)
{
add+=gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j];
}
}
return add;
}
//计算雅可比矩阵
void GetJaccobi()
{
float JH[Bus_Num-1][Bus_Num-1];
float JN[Bus_Num-1][Bus_Num-1];
float JJ[Bus_Num-1][Bus_Num-1];
float JL[Bus_Num-1][Bus_Num-1];
float JR[Bus_Num-1][Bus_Num-1];
float JS[Bus_Num-1][Bus_Num-1];
int i,j,m,n,k; //m,n用于建立雅可比矩阵，k用于雅可比矩阵输出
for(i=0;i<Bus_Num-1;i++)
{
for(j=0;j<Bus_Num-1;j++)
{
JH[i][j]=0;
JN[i][j]=0;
JJ[i][j]=0;
JL[i][j]=0;
JR[i][j]=0;
JS[i][j]=0;
}
}
for(i=0;i<Bus_Num-1;i++)
{
for(j=0;j<Bus_Num-1;j++)
{
if ((i!=j)&&(gBus[j].Type==0))
{
JH[i][j]=-gY_B[i][j]*ge[i]+gY_G[i][j]*gf[i];
JN[i][j]=gY_G[i][j]*ge[i]+gY_B[i][j]*gf[i];
JJ[i][j]=-JN[i][j];
JL[i][j]=JH[i][j];
m=2*i;
n=2*j;
gJaccobi[m][n]=JH[i][j];
gJaccobi[m][n+1]=JN[i][j];
gJaccobi[m+1][n]=JJ[i][j];
gJaccobi[m+1][n+1]=JL[i][j];
}
if ((i==j)&&(gBus[j].Type==0))
{
JH[i][j]=2*gY_G[i][j]*gf[i]+add1(i,j); //调用自定义函数1
JN[i][j]=2*gY_G[i][j]*ge[i]+add2(i,j); //调用自定义函数2
JJ[i][j]=-2*gY_B[i][j]*gf[i]+add2(i,j); //调用自定义函数2
JL[i][j]=-2*gY_B[i][j]*ge[i]-add1(i,j); //调用自定义函数1
m=2*i;
n=2*i;
gJaccobi[m][n]=JH[i][j];
gJaccobi[m][n+1]=JN[i][j];
gJaccobi[m+1][n]=JJ[i][j];
gJaccobi[m+1][n+1]=JL[i][j];
}
if ((i!=j)&&(gBus[j].Type==1))
{
JH[i][j]=-gY_B[i][j]*ge[i]+gY_G[i][j]*gf[i];
JN[i][j]=gY_G[i][j]*ge[i]+gY_B[i][j]*gf[i];
JR[i][j]=0;
JS[i][j]=0;
m=2*i;
n=2*j;
gJaccobi[m][n]=JH[i][j];
gJaccobi[m][n+1]=JN[i][j];
gJaccobi[m+1][n]=JR[i][j];
gJaccobi[m+1][n+1]=JS[i][j];
}
if ((i==j)&&(gBus[j].Type==1))
{
JH[i][j]=2*gY_G[i][j]*gf[i]+add1(i,j); //调用自定义函数1
JN[i][j]=2*gY_G[i][j]*ge[i]+add2(i,j); //调用自定义函数2
JR[i][j]=2*gf[j];
JS[i][j]=2*ge[j];
m=2*i;
n=2*i;
gJaccobi[m][n]=JH[i][j];
gJaccobi[m][n+1]=JN[i][j];
gJaccobi[m+1][n]=JR[i][j];
gJaccobi[m+1][n+1]=JS[i][j];
}
}
}
printf(" the Jaccobi Matrix is calculated below:\n");
for (i=0,k=0;i<2*(Bus_Num-1);i++)
{
for (j=0;j<2*(Bus_Num-1);j++)
{
printf("%7.3f\t",gJaccobi[i][j]);
k++;
(k%6==0)?printf("\n"):printf("");
}
}
printf("\n");
printf("\n");
}
//计算电压不平衡量
void GetRevised()
{
int i=0,j=0,k=0;
NEquation ob2;
ob2.SetSize(2*(Bus_Num-1));
for (i=0;i<2*(Bus_Num-1);i++)
{
for (j=0;j<2*(Bus_Num-1);j++)
{
ob2.Data(i,j)=gJaccobi[i][j];
}
}
/* //显示方程系数矩阵
printf(" the OBJ of the equation is calculated below:\n");
for (i=0,k=0;i<2*(Bus_Num-1);i++)
{
for (j=0;j<2*(Bus_Num-1);j++)
{
printf("ob2.Data(%d,%d)=%7.3f\t",i+1,j+1,ob2.Data(i,j));
k++;
(k%6==0)?printf("\n"):printf("");
}
}
*/
printf("\n");
printf("\n");
for (i=0;i<2*(Bus_Num-1);i++)
{
ob2.Value(i)=gDelta_PQ[i];
printf("ob2.Value(%d)=%7.3f\n",i+1,ob2.Value(i));
}
ob2.Run();
/*
printf("\n");
printf(" the root of the equation is calculated below:\n");
for (i=0;i<2*(Bus_Num-1);i++)
{
printf("ob2.Value(%d)=%7.3f\n",i+1,ob2.Value(i));
}
*/
for (i=0;i<2*(Bus_Num-1);i++)
{
gDelta_fe[i]=ob2.Value(i);
}
printf(" the gDelta_fe Matrix is calculated below:\n");
for (i=0,k=0;i<2*(Bus_Num-1);i++)
{
printf("gDelta_fe[%d]=%7.3f\n",i+1,gDelta_fe[i]);
}
printf("\n");
printf("\n");
for (i=0,k=0;i<2*(Bus_Num-1);i++,i++)
{
k=i/2;
gDelta_f[k]=gDelta_fe[i];
gDelta_e[k]=gDelta_fe[i+1];
}
printf(" the gDelta_e && gDelta_f Matrix is calculated below\n");
for (i=0;i<Bus_Num-1;i++)
{
printf("gDelta_e[%d]=%7.3f\tgDelta_f[%d]=%7.3f\n",
i+1,gDelta_e[i],i+1,gDelta_f[i]);
}
printf("\n");
printf("\n");
}
//计算电压新值
void GetNewValue()
{
int i;
for(i=0;i<Bus_Num-1;i++)
{
gf[i]+=gDelta_f[i];
ge[i]+=gDelta_e[i];
}
}
//主函数
int main(int argc, char* argv[])
{
int i,j,Count_Num; //定义循环变量和循环次数变量
float maxValue; //定义电压偏差最小值变量
FILE *fp;
// test(); //调用测试方程函数
GetData(); //调用读取数据函数
GetYMatrix(); //调用节点导纳矩阵函数
SetInitial(); //调用电压初始化函数
//调试程序用
GetUnbalance();
GetJaccobi();
// test();
GetRevised();
GetNewValue();
/*
for(Count_Num=0;Count_Num<=20;Count_Num++) //迭代循环
{
GetUnbalance(); //调用功率不平衡量函数
GetJaccobi(); //调用雅可比函数
GetRevised(); //调用电压修正函数
GetNewValue(); //调用电压新值函数
maxValue=gDelta_fe[0]; //电压偏差最大值赋初值
for(i=0;i<=2*(Bus_Num-1)-1;i++)
{
if(maxValue<gDelta_fe[i])
{
maxValue=gDelta_fe[i];
}
}
if(maxValue<Precision)
{
break;
}
}
*/
printf("Count_Num=%d\n",Count_Num);
printf(" the final voltage is calculated below:\n");
if((fp=fopen("D:\\潮流计算080607\\MQY flow\\PVflow\\data\\Bus Voltage.txt","w"))==NULL) //确定输出文件路径
{
printf("Can not open the file named 'Bus Voltage.txt' \n");
}
for(i=0;i<=Bus_Num-1;i++) //计算最终电压
{
printf("U[%d]=%f\tU_angle[%d]=%f\n",i+1,sqrt(ge[i]*ge[i]+gf[i]*gf[i]),
i+1,atan(gf[i]/ge[i])*(180/3.1415926));
fprintf(fp,"U[%d]=%f\tU_angle[%d]=%f\n",i+1,sqrt(ge[i]*ge[i]+gf[i]*gf[i]),
i+1,atan(gf[i]/ge[i])*(180/3.1415926));
}
fclose(fp);
//计算平衡节点功率
float Ps,Qs; //定义平衡节点功率变量
Ps=0; //给平衡节点功率赋初值
Qs=0;
for(i=0;i<Bus_Num;i++) //计算平衡节点功率
{
for(j=0;j<Bus_Num;j++)
{
if(gBus[i].Type==2) //平衡节点条件判断
{
Ps+=gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j];
Qs+=-gY_B[i][j]*ge[j]-gY_G[i][j]*gf[j];
}
}
if(gBus[i].Type==2)
{
Ps=Ps*ge[i]-Qs*gf[i];
Qs=Ps*gf[i]+Qs*ge[i];
}
}
//屏幕输出平衡节点功率
printf(" the power of the balance node is calculated below:\n");
printf("Ps=%f\t\tQs=%f\n",Ps,Qs);
//计算线路功率
float LP1[Bus_Num][Bus_Num],LQ1[Bus_Num][Bus_Num]; //定义线路功率数组
float LP2[Bus_Num][Bus_Num],LQ2[Bus_Num][Bus_Num];
for(i=0;i<Bus_Num;i++)
{
for(j=0;j<Bus_Num;j++)
{ //给线路功率数组赋初值
LP1[i][j]=0;
LQ1[i][j]=0;
LP2[i][j]=0;
LQ2[i][j]=0;
}
}
for(i=0;i<Bus_Num;i++) //计算线路功率数组
{
for(j=0;j<Bus_Num;j++)
{
LP1[i][j]=-ge[i]*(gY_G[i][j]*(ge[i]-ge[j])+gY_B[i][j]*(gf[j]-gf[i]))+
gf[i]*(gY_B[i][j]*(ge[j]-ge[i])+gY_G[i][j]*(gf[j]-gf[i]));
LQ1[i][j]=-ge[i]*(gY_B[i][j]*(ge[j]-ge[i])+gY_G[i][j]*(gf[j]-gf[i]))-
gf[i]*(gY_G[i][j]*(ge[i]-ge[j])+gY_B[i][j]*(gf[j]-gf[i]));
LP2[j][i]=-ge[j]*(gY_G[j][i]*(ge[j]-ge[i])+gY_B[j][i]*(gf[i]-gf[j]))+
gf[j]*(gY_B[j][i]*(ge[i]-ge[j])+gY_G[j][i]*(gf[i]-gf[j]));
LQ2[j][i]=-ge[j]*(gY_B[j][i]*(ge[i]-ge[j])+gY_G[j][i]*(gf[i]-gf[j]))-
gf[j]*(gY_G[j][i]*(ge[j]-ge[i])+gY_B[j][i]*(gf[i]-gf[j]));
}
}
//屏幕输出线路功率数组
printf("\n\n the power of the line is calculated below:\n\n");
for(i=0;i<Bus_Num;i++)
{
for(j=0;j<Bus_Num;j++)
{
printf("LP1[%d][%d]=%8.5f\tLQ1[%d][%d]=%8.5f\n",i+1,j+1,LP1[i][j],i+1,j+1,LQ1[i][j]);
printf("LP2[%d][%d]=%8.5f\tLQ2[%d][%d]=%8.5f\n",i+1,j+1,LP2[i][j],i+1,j+1,LQ2[i][j]);
}
}
float gDelta_LP[Line_Num],gDelta_LQ[Line_Num]; //定义线路功率损耗数组
float gDelta_LPT=0,gDelta_LQT=0; //定义线路总功率损耗
for(i=0;i<Line_Num;i++) //计算线路损耗
{
gDelta_LP[i]=LP1[gLine[i].No_I-1][gLine[i].No_J-1]+LP2[gLine[i].No_J-1][gLine[i].No_I-1];
gDelta_LQ[i]=LQ1[gLine[i].No_I-1][gLine[i].No_J-1]+LQ2[gLine[i].No_J-1][gLine[i].No_I-1];
}
if((fp=fopen("D:\\潮流计算080607\\MQY flow\\PVflow\\data\\line loss.txt","w"))==NULL) //确定输出文件路径
{
printf("Can not open the file named 'line loss.txt' \n");
}
//屏幕显示线路功率损耗
printf("\n\n the loss of the line is calculated below:\n");
for(i=0;i<Line_Num;i++)
{
printf("gDelta_LP[%d]=%8.5f\t\tgDelta_LQ[%d]=%8.5f\n",
i+1,gDelta_LP[i],i+1,gDelta_LQ[i]);
fprintf(fp,"gDelta_LP[%d]=%8.5f\t\tgDelta_LQ[%d]=%8.5f\n",
i+1,gDelta_LP[i],i+1,gDelta_LQ[i]);
}
fprintf(fp,"Ps=%f\tQs=%f\n",Ps,Qs);
fclose(fp);
for(i=0;i<Line_Num;i++)
{
gDelta_LPT+=gDelta_LP[i];
gDelta_LQT+=gDelta_LQ[i];
}
//屏幕显示线路总功率损耗
printf("\n\n the total loss of the line is calculated below:\n");
printf("%f+j%f\n",gDelta_LPT,gDelta_LQT);
return 0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -