华北电力大学 电力系统稳态潮流上机计算程序结果

// loadflow.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#define NODE_TOTAL_NUM 9
struct NodeData
{
unsigned int Index; //node index
unsigned int Type; //node type: PQ:1,PV:2,balance point:0
float FirstInput; //PQ or PV:active power,balance point:V
float SecondInput; //PQ:reactive power,PV:V,balance point:angle
};
struct NodeData gNodeData[NODE_TOTAL_NUM]; // define a global variable of node


#define NODE_DATA_FILENAME "E:\\loadflow\\node.txt"
#define LINE_TOTAL_NUM 9


struct LineData
{
unsigned int Index; //line index
unsigned int Type; //line type: Line:1,tranformer:2,line connected to ground:3
unsigned int Status; //line status: on:1,off:0
unsigned int SrcNode; //source node index of this line
unsigned int DestNode;//destination node index of this line
float Resistance; //resistance of the line
float Reactance; //reactance of the line
float Conductance; //Line: conductance B/2,transformer:change ratio
};
struct LineData gLineData[LINE_TOTAL_NUM];
#define LINE_DATA_FILENAME "E:\\loadflow\\line.txt"


void SolveEquation(unsigned int Dimension,float FactorMatrix[14][14], float ConstVector[14]) //求解方程的函数 Ax=B
{
unsigned int i,j,k;

//消元过程
for(i=0;i{
//规格化过程:(每次消元前,先规格化,以便以下各行消元时,消元系数直接取待消
//列元素的值即可,也便于回代过程,而运算量并不增加)
for( j = i+1; j < Dimension; j++ )
{
FactorMatrix[i][j] = FactorMatrix[i][j] / FactorMatrix[i][i];
}
ConstVector[i] = ConstVector[i]/FactorMatrix[i][i];


for( j = i+1; j < Dimension; j++ ) //消去第i列(从i+1行到最后一行)
{
if( FactorMatrix[j][i] != 0 ) //如果第j行第i列元素本就是0,则不需本列对应的消元过程
{
for( k = i + 1; k < Dimension; k++ ) //当FactorMatrix[i][k]=0,a[j][k]值不变,可省去运算
if( FactorMatrix[i][k] != 0 )
FactorMatrix[j][k] = FactorMatrix[j][k] - FactorMatrix[j][i] * FactorMatrix[i][k];
}
//常数项的消元
ConstVector[j] = ConstVector[j] - FactorMatrix[j][i]* ConstVector[i];
}
}

//回代过程
for( i = Dimension-1; i > 0; i-- ) //Dimension-2:最后一个变量可直接获得,从n-1个变量求起
{
for(j = Dimension-1; j > i-1 ; j-- )
{
if( FactorMatrix[i-1][j] != 0 )
ConstVector[i-1] = ConstVector[i-1] - FactorMatrix[i-1][j] * ConstVector[j]; //a[i][k]=0时可以不算
}
}
return;
}


int main(int argc, char* argv[])
{
printf("Hello World!\n");
//第一阶段 输入节点和线路数据

FILE* fp; //file structure
int i;
fp=fopen(NODE_DATA_FILENAME,"rb"); //open and read the node data输入节点数据
printf("\n节点数据\n");
for(i=0;i{
fscanf(fp,"%d %d %f %f",&gNodeData[i].Index,
&gNodeData[i].Type,
&gNode

Data[i].FirstInput,
&gNodeData[i].SecondInput);
printf("%d %d %f %f\n",gNodeData[i].Index,
gNodeData[i].Type,
gNodeData[i].FirstInput,
gNodeData[i].SecondInput);
}
fclose(fp);



fp=fopen(LINE_DATA_FILENAME,"r"); //open and read the line data线路数据
printf("\n线路数据\n");
for(i=0;i{
fscanf(fp,"%u%u%u%u%u%f%f%f", &gLineData[i].Index,
&gLineData[i].Type,
&gLineData[i].Status,
&gLineData[i].SrcNode,
&gLineData[i].DestNode,
&gLineData[i].Resistance,
&gLineData[i].Reactance,
&gLineData[i].Conductance);
printf("%u %u %u %u %u %f %f %f\n",gLineData[i].Index,
gLineData[i].Type,
gLineData[i].Status,
gLineData[i].SrcNode,
gLineData[i].DestNode,
gLineData[i].Resistance,
gLineData[i].Reactance,
gLineData[i].Conductance);

}
fclose(fp);

//第二阶段 计算线路导纳矩阵
struct daona //结构体存放导纳元素中的电导和电纳
{
float diandao;//电导
float dianna;//电纳

};
struct daona D[NODE_TOTAL_NUM][NODE_TOTAL_NUM];//定义一个二维数组存放导纳元素并对定义的数组赋0
int j;
for (i=0;ifor (j=0;j{
D[i][j].dianna=0;
D[i][j].diandao=0;
}

for (i=0;i{
float tab;
float lineg,lineb;
int x,y;
x=gLineData[i].SrcNode-1; //第i条线路首尾节点与 数组行列 的对应关系
y=gLineData[i].DestNode-1;
tab=pow(gLineData[i].Resistance,2)+pow(gLineData[i].Reactance,2);//第i条线路电阻与电抗平方之和
lineg=gLineData[i].Resistance/tab; //第i条线路电导
lineb=(-1)*gLineData[i].Reactance/tab; //第i条线路电纳


if(gLineData[i].Type==1)//对线路为导线的类型
{
D[x][y].dianna+=(-1)*lineb;//首尾节点之间互导纳中的电纳
D[x][y].diandao+=(-1)*lineg;//首尾节点之间互导纳中的电导

D[x][x].diandao+=lineg;//首节点自导纳中电导
D[x][x].dianna+=lineb+gLineData[i].Conductance;//首节点自导纳中电纳

D[y][y].dianna+=lineb+gLineData[i].Conductance;//尾节点自导纳中电纳

D[y][y].diandao+=lineg;//尾节点自导纳中电导

D[y][x].dianna=D[x][y].dianna;//首节点与尾节点之间导纳与尾节点与首节点之间导纳相等
D[y][x].diandao=D[x][y].diandao;
}


if(gLineData[i].Type==2)//对线路为变压器的类型 对变压器进行π型等值 计算方法与导线型相同
{
D[x][y].dianna+=(-1)*lineb/gLineData[i].Conductance;
D[x][y].diandao+=(-1)*lineg/gLineData[i].Conductance;

D[x][x].diandao+=lineg/gLineData[i].Conductance+(1-gLineData[i].Conductance)*lineg/(gLineData[i].Conductance*gLineData[i].Conductance);
D[x][

x].dianna+=lineb/gLineData[i].Conductance+(1-gLineData[i].Conductance)*lineb/(gLineData[i].Conductance*gLineData[i].Conductance);

D[y][y].dianna+=lineb/gLineData[i].Conductance+(gLineData[i].Conductance-1)*lineb/gLineData[i].Conductance;
D[y][y].diandao+=lineg/gLineData[i].Conductance+(gLineData[i].Conductance-1)*lineg/gLineData[i].Conductance;

D[y][x].dianna=D[x][y].dianna;
D[y][x].diandao=D[x][y].diandao;
}
}
printf("\n导纳矩阵中电导分量\n");

for (i=0;i<9;i++)
for(j=0;j<9;j++)
{
printf("%.3f ",D[i][j].diandao);//输出导纳矩阵中电导分量
if(j==8)
{printf("\n");}
}

printf("\n导纳矩阵中电纳分量\n");
for (i=0;i<9;i++)
for(j=0;j<9;j++)
{
printf("%.3f ",D[i][j].dianna);//输出导纳矩阵中电纳分量
if(j==8)
{printf("\n");}
}


//第三阶段 对节点电压赋初值

struct V
{
double U;
double O;
};
struct V dianya[NODE_TOTAL_NUM];//将节点电压用极坐标表示，并赋予初值
for (i=0;i{
switch (gNodeData[i].Type)
{
case 1:
dianya[i].U=1.0;//对PQ节点电压幅值相位都未知
dianya[i].O=0.0;
break;
case 2:
dianya[i].U=gNodeData[i].SecondInput;//对PV节点电压幅值已知，相位未知
dianya[i].O=0.0;
break;
case 0:
dianya[i].U=gNodeData[i].FirstInput;//平衡节点电压幅值相位都已知
dianya[i].O=gNodeData[i].SecondInput;
break;

}
}
printf(" 电压相位\n");
for(i=0;i{
printf("%f\n%f\n",dianya[i].U,dianya[i].O);
}
//第四阶段 计算节点PQ的不平衡量
int k;
for(k=0;k<20;k++)
{
printf("第%d次叠代\n",k+1);



float BPH[14];
int k=0,l=0;
for (i=0;i<9;i++)//对PQ节点都有有功P的不平衡量，计算得P的不平衡量并依次放在数组中
{

if(gNodeData[i].Type==1)
{
float N=0;
for(j=0;j<9;j++)
{
N+=dianya[i].U*dianya[j].U*(D[i][j].diandao*cos(dianya[i].O-dianya[j].O)+D[i][j].dianna*sin(dianya[i].O-dianya[j].O));


}

BPH[l]=gNodeData[i].FirstInput-N;
l++;

}
if(gNodeData[i].Type==2)
{
float N=0;
for(j=0;j<9;j++)
{
N+=dianya[i].U*dianya[j].U*(D[i][j].diandao*cos(dianya[i].O-dianya[j].O)+D[i][j].dianna*sin(dianya[i].O-dianya[j].O));

}

BPH[l]=gNodeData[i].FirstInput-N;
l++;

}
}



for(i=0;i<9;i++)//对PV节点只有无功Q的不平衡量，计算得Q的不平衡量并依次放在数组中
{
float M=0;
if(gNodeData[i].Type==1)
{
for(j=0;j<9;j++)
M+=dianya[i].U*dianya[j].U*(D[i][j].diandao*sin(dianya[i].O-dianya[j].O)-D[i][j].dianna*cos(dianya[i].O-dianya[j].O));
BPH[l]=gNodeData[i].SecondInput-M;
l++;
}

}


printf("\nPQ不平衡量\n");

for(i=0;i<14;i++)
{
printf("%f\n",BPH[i]);
}

float det=0;
for(i=0;i<8;i++)//对循环阶段，计算最大不平衡量，选择是否跳出循环迭代
{
if(fabs(BPH[i])>det)
det=fabs(BPH[i]);
}

if(det<0.000001)
break;





//第

五阶段 计算雅可比矩阵
float H[8][8],N[8][6],J[6][8],L[6][6];//将雅克比矩阵中HNJL元素分离，并分别存放在四个数组中

for(i=0;i<8;i++)//对HNJL数组各元素进行赋0
for(j=0;j<8;j++)
H[i][j]=0;

for(i=0;i<8;i++)
for(j=0;j<6;j++)
N[i][j]=0;

for(i=0;i<6;i++)
for(j=0;j<8;j++)
J[i][j]=0;

for(i=0;i<6;i++)
for(j=0;j<6;j++)
L[i][j]=0;


int m=0,n=0,s=0,z=0;


for(i=0;i<9;i++)//计算雅可比矩阵中H元素，除平衡节点，每个节点与所有节点（除平衡节点）对应，可得H元素
{
if(gNodeData[i].Type!=0)
{
for(n=0;n<8;n++)
{

if(i!=n)
{
H[i][n]=dianya[i].U*dianya[n].U*(D[i][n].diandao*sin(dianya[i].O-dianya[n].O)-D[i][n].dianna*cos(dianya[i].O-dianya[n].O));
}
else
{
for(s=0;s<9;s++)
{
if(s!=i)
{
H[i][n]+=(-1)*dianya[i].U*dianya[s].U*(D[i][s].diandao*sin(dianya[i].O-dianya[s].O)-D[i][s].dianna*cos(dianya[i].O-dianya[s].O));
}
}
}
}
}
}


for(i=0;i<9;i++)//计算雅可比矩阵中N元素 除平衡节点，每个节点与PQ节点对应，可得N元素
{
if(gNodeData[i].Type==1||gNodeData[i].Type==2)
{
n=0;
for(z=0;z<9;z++)
if(gNodeData[z].Type==1)
{
if(i!=z)
{
N[i][n]=dianya[i].U*dianya[z].U*(D[i][z].diandao*cos(dianya[i].O-dianya[z].O)+D[i][z].dianna*sin(dianya[i].O-dianya[z].O));
}
if(i==z)
{
for(s=0;s<9;s++)
{
if(s!=i)
{
N[i][n]+=dianya[i].U*dianya[s].U*(D[i][s].diandao*cos(dianya[i].O-dianya[s].O)+D[i][s].dianna*sin(dianya[i].O-dianya[s].O));
}
}
N[i][n]+=2*dianya[i].U*dianya[i].U*D[i][i].diandao;
}
n++;
}
}
}


m=0;

for(i=0;i<9;i++)//计算雅可比矩阵中J元素 PQ节点与所有节点（除平衡节点）对应，可得J元素
{
if(gNodeData[i].Type==1)
{
for(n=0;n<8;n++)
{
if(i!=n)
{
J[m][n]=(-1)*dianya[i].U*dianya[n].U*(D[i][n].diandao*cos(dianya[i].O-dianya[n].O)+D[i][n].dianna*sin(dianya[i].O-dianya[n].O));
}
if(i==n)
{
for(s=0;s<9;s++)
{
if(s!=i)
{
J[m][n]+=dianya[i].U*dianya[s].U*(D[i][s].diandao*cos(dianya[i].O-dianya[s].O)+D[i][s].dianna*sin(dianya[i].O-dianya[s].O));
}
}
}
}
m++;
}
}


m=0;
z=0;

for(i=0;i<9;i++)//计算雅可比矩阵中L元素 PQ节点与PQ节点对应 可得L元素
{
if(gNodeData[i].Type==1)
{z=0;
for(n=0;n<9;n++)
{
if(gNodeData[n].Type==1)
{
if(i!=n)
{
L[m][z]=dianya[i].U*dianya[n].U*(D[i][n].diandao*sin(dianya[i].O-dianya[n].O)-D[i][n].dianna*cos(dianya[i].O-dianya[n].O));
}
else
{
for(s=0;s<9;s++)
{
if(s!=i)
{
L[m][z]+=dianya[i].U*dianya[s].U*(D[i][s].diandao*sin(dianya[i].O-dianya[s].O)-D[i][s].dianna*cos(dianya[i].O-dianya[s].O));
}
}
L[m][z]+=(-2)*dianya[i].U*dianya[i].U*D[i][i].dianna;
}
z++;

}
}

m++;
}
}




printf("\n H阵\n");//输出H阵


for(i=0;i<8;i++)
{
for(j=0;j<8;j++)
printf("%.3f ",H[i][j]);
printf("\n");
}
printf("\n N阵\n");//输出N阵

for(i=0;i<8;i++)
{
for(j=0;j<6;j++)
printf("%.3f ",N[i][j]);
printf("\n");
}
printf("\n J阵\n");//输出J阵

for(i=0;i<6;i++)
{
for(j=0;j<8;j++)
printf("%.3f ",J[i][j]);
printf("\n");
}

printf("\n L阵\n");//输出L阵



for(i=0;i<6;i++)
{
for(j=0;j<6;j++)
printf("%.3f ",L[i][j]);
printf("\n");
}




float ykb[14][14];//根据求得的H，N，J，L子阵合成一个二维数组，形成雅克比矩阵
for(i=0;i<14;i++)
for(j=0;j<14;j++)
{
if(i<8&&j<8)
ykb[i][j]=H[i][j];
if(i<8&&j>=8)
ykb[i][j]=N[i][j-8];
if(i>=8&&j<8)
ykb[i][j]=J[i-8][j];
if(i>=8&&j>=8)
ykb[i][j]=L[i-8][j-8];
}


printf("\n 雅克比矩阵为\n");

for(i=0;i<14;i++)
{
for(j=0;j<14;j++)
printf("%.3f ",ykb[i][j]);
printf("\n");
}


// 第六阶段 解方程，返回电压的不平衡量

SolveEquation(14,ykb,BPH);
for(i=0;i<8;i++)
dianya[i].O+=BPH[i];//根据返回的电压不平衡量，修正节点电压相位

for(j=0;j<9;j++)
{
if(gNodeData[j].Type==1)
{
dianya[j].U+=BPH[i]*dianya[j].U;//根据返回的电压不平衡量，修正节点电压幅值
i++;
}

}
printf(" 电压修正结果\n");
for(i=0;i{
printf("%f\n%f\n",dianya[i].U,dianya[i].O);//计算得到新的电压幅值相位后，带回第四阶段开始 继续计算
}


}
//第七阶段 迭代循环 得到较准确的电压幅值和相位



//第八阶段 循环迭代之后 计算各节点注入功率
struct IN
{
float P;
float Q;
};
struct IN ZHRU[9];
for(i=0;i<9;i++)
{
ZHRU[i].P=0;
ZHRU[i].Q=0;
}
for(i=0;i<9;i++)
{
if(gNodeData[i].Type==0)//对平衡节点进行计算，得到注入P，Q
for(j=0;j<9;j++)
{
ZHRU[i].P+=dianya[i].U*dianya[j].U*(D[i][j].diandao*cos(dianya[i].O-dianya[j].O)+D[i][j].dianna*sin(dianya[i].O-dianya[j].O));
ZHRU[i].Q+=dianya[i].U*dianya[j].U*(D[i][j].diandao*sin(dianya[i].O-dianya[j].O)-D[i][j].dianna*cos(dianya[i].O-dianya[j].O));
}
if(gNodeData[i].Type==1)//对PQ节点，PQ都已知
{
ZHRU[i].P=gNodeData[i].FirstInput;
ZHRU[i].Q=gNodeData[i].SecondInput;
}
if(gNodeData[i].Type==2)//对PV节点，P已知，Q未知，计算节点注入无功功率Q
{
ZHRU[i].P=gNodeData[i].FirstInput;
for(j=0;j<9;j++)
ZHRU[i].Q+=dianya[i].U*dianya[j].U*(D[i][j].diandao*sin(dianya[i].O-dianya[j].O)-D[i][j].dianna*cos(dianya[i].O-dianya[j].O));

}

}
for(i=0;i<9;i++)
{
printf("\n %d节点注入有功%f 无功%f\n",i+1,ZHRU[i].P,ZHRU[i].Q);
}
printf("\n");
//第九阶段 计算各线路损耗
struct IN SUNHAO[9];
for(i=0;i<9;i++)
{
SUNHAO[i].P=0;
SUNHAO[i].Q=0;
}

int s,e;
for(i=0;i<9;i++)
{
s=gLineData[i].SrcNode-1;//首节点
e=gLineData[i].DestNode-1;//

尾节点
if(gLineData[i].Type==2)
{
SUNHAO[i].P=(dianya[s].U*dianya[s].U-2*dianya[s].U*dianya[e].U*cos(dianya[s].O-dianya[e].O)+dianya[e].U*dianya[e].U)*D[s][e].diandao*(-1);
SUNHAO[i].Q=(-1)*(dianya[s].U*dianya[s].U-2*dianya[s].U*dianya[e].U*cos(dianya[s].O-dianya[e].O)+dianya[e].U*dianya[e].U)*D[s][e].dianna*(-1);
}
if(gLineData[i].Type==1)
{
SUNHAO[i].P=(dianya[s].U*dianya[s].U-2*dianya[s].U*dianya[e].U*cos(dianya[s].O-dianya[e].O)+dianya[e].U*dianya[e].U)*D[s][e].diandao*(-1);
SUNHAO[i].Q=(-1)*(dianya[s].U*dianya[s].U-2*dianya[s].U*dianya[e].U*cos(dianya[s].O-dianya[e].O)+dianya[e].U*dianya[e].U)*D[s][e].dianna*(-1)-dianya[s].U*dianya[s].U*gLineData[i].Conductance-dianya[e].U*dianya[e].U*gLineData[i].Conductance;

}
}
for(i=0;i<9;i++)
{
printf("\n %d线路损耗有功%f 无功%f\n",i+1,SUNHAO[i].P,SUNHAO[i].Q);
}
float ZP=0,ZQ=0;
for(i=0;i<9;i++)
{
ZP+=SUNHAO[i].P;
ZQ+=SUNHAO[i].Q;
}
printf("线路总网损\n有功 %f\n无功 %f",ZP,ZQ);

return 0;
}