电力系统分析潮流计算程序.docx
《电力系统分析潮流计算程序.docx》由会员分享,可在线阅读,更多相关《电力系统分析潮流计算程序.docx(19页珍藏版)》请在冰豆网上搜索。
电力系统分析潮流计算程序
#include
#include
#include
#definePI3.14159
//节点参数结构体
structNodeType
{
intN;//节点号
intType;//节点类型
doublee;//电压幅值
doublef;//电压相角
doublePd;//负荷有功
doubleQd;//负荷无功
doublePs;//出力有功
doubleQs;//出力无功
doubleBc;//并联电容的电抗值
};
//支路参数结构体
structBranchType
{
intNbr;//支路号
intNl;//首节点
intNr;//末节点
doubleR;//支路电阻
doubleX;//支路电抗
doubleBn;//对地电抗
doubleKt;//支路变比
};
//******************************************全局变量声明***************************************
intn;//节点数
intnPQ;//PQ节点数
intnPV;//PV节点数
intnbr;//支路数
intng;//发电机台数
intMark=0;//标记支路参数是否已经转换
double**G;//导纳矩阵G部分
double**B;//导纳矩阵B部分
double*dS;//功率不平衡量
double*mid1,*mid2;//求功率不平衡量时的中间变量
double*Us;//电压初值
doubleerror=1;//误差值
doubleiteration=0.000001;//误差精度
double**Jacob;//雅克比矩阵
double**invJac;//雅克比矩阵的逆
double*dfe;//节点电压修正值
structNodeType*Node;//读入时的节点参数结构体
structBranchType*Branch;//读入时的支路参数结构体
//*********************************************主程序******************************************
voidmain()
{
voidLoadData();
voidFormY();//形成导纳矩阵
voidDeltaS();//求功率不平衡量
voidFormJacob();//形成雅克比矩阵
voidInvJac();//求雅克比矩阵的逆
voidUpdateU();//修正电压值
voidCalculatePQ();
voidPrint1(double*,int);
voidPrint2(double**,int,int);
intkk;//迭代次数
LoadData();
FormY();
printf("iteration=%lf\n",iteration);
kk=0;
DeltaS();
while(error>iteration&&kk<50)
{
FormJacob();
UpdateU();
DeltaS();
kk++;
}
printf("迭代次数为%4d\n",kk);
CalculatePQ();
printf("error=%e\n",error);
}
//*********************************************************************************************
voidLoadData()
{
inti,j;
inttN,tType;
doublete,tf,tPd,tQd,tPs,tQs,tBc;//用于重新排列节点信息的临时变量
FILE*fp;//文件指针
charfilename[50]={""};
printf("请输入数据文件名:
");
scanf("%s",filename);
if((fp=fopen(filename,"r"))==NULL)
{
printf("cannotopenthefile:
data.txt\n");
return;
}
fscanf(fp,"%d",&n);
printf("节点个数为:
%d\n",n);
//为节点参数申请空间
Node=(structNodeType*)malloc(sizeof(structNodeType)*n);
//读取节点参数
printf("调整前的节点参数为:
\n");
for(i=0;ifscanf(fp,"%d%d%lf%lf%lf%lf%lf%lf%lf",&Node[i].N,&Node[i].Type,&Node[i].e,&Node[i].f,&Node[i].Pd,&Node[i].Qd,&Node[i].Ps,&Node[i].Qs,&Node[i].Bc);
//计算PQ节点和PV节点的个数
for(i=0;i{
if(Node[i].Type==1)
nPQ++;
elseif(Node[i].Type==2)
nPV++;
}
printf("PQ节点个数:
%d\n",nPQ);
printf("PV节点个数:
%d\n",nPV);
//重新排列节点参数(冒泡法)
for(j=0;jfor(i=0;i{
if(Node[i].Type>Node[i+1].Type)
{
tN=Node[i].N;Node[i].N=Node[i+1].N;Node[i+1].N=tN;
tType=Node[i].Type;Node[i].Type=Node[i+1].Type;Node[i+1].Type=tType;
te=Node[i].e;Node[i].e=Node[i+1].e;Node[i+1].e=te;
tf=Node[i].f;Node[i].f=Node[i+1].f;Node[i+1].f=tf;
tPd=Node[i].Pd;Node[i].Pd=Node[i+1].Pd;Node[i+1].Pd=tPd;
tQd=Node[i].Qd;Node[i].Qd=Node[i+1].Qd;Node[i+1].Qd=tQd;
tPs=Node[i].Ps;Node[i].Ps=Node[i+1].Ps,Node[i+1].Ps=tPs;
tQs=Node[i].Qs;Node[i].Qs=Node[i+1].Qs;Node[i+1].Qs=tQs;
tBc=Node[i].Bc;Node[i].Bc=Node[i+1].Bc;Node[i+1].Bc=tBc;
}
}
//为电压初值申请空间
Us=(double*)malloc(sizeof(double)*(n-1));
for(i=0;iUs[i]=Node[i].e;
//读取支路参数
fscanf(fp,"%d",&nbr);
printf("支路个数为:
%d\n",nbr);
//为支路参数申请空间
Branch=(structBranchType*)malloc(sizeof(structBranchType)*nbr);//读入的支路参数结构体
//读入支路参数
for(i=0;ifscanf(fp,"%d%d%d%lf%lf%lf%lf",&Branch[i].Nbr,&Branch[i].Nl,&Branch[i].Nr,&Branch[i].R,&Branch[i].X,&Branch[i].Bn,&Branch[i].Kt);
//支路节点号参数调整
for(i=0;i{
Mark=0;
for(j=0;j{
if(Branch[i].Nl==Node[j].N&&Mark==0)
{
Branch[i].Nl=j+1;
Mark=1;
}
}
}
for(i=0;i{
Mark=0;
for(j=0;j{
if(Branch[i].Nr==Node[j].N&&Mark==0)
{
Branch[i].Nr=j+1;
Mark=1;
}
}
}
fclose(fp);
}
//***********************************************************************************************************************
//*********************************************形成导纳矩阵**************************************************************
voidFormY()
{
inti,j;
doubleZ2;//存储Z^2=R^2+X^2
G=(double**)malloc(sizeof(double*)*n);//为G申请空间
B=(double**)malloc(sizeof(double*)*n);//为B申请空间
for(i=0;i{
G[i]=(double*)malloc(sizeof(double)*n);
B[i]=(double*)malloc(sizeof(double)*n);
}
//初始化G、B
for(i=0;ifor(j=0;j{
G[i][j]=0;
B[i][j]=0;
}
//计算非对角元素
for(i=0;i{
Z2=Branch[i].R*Branch[i].R+Branch[i].X*Branch[i].X;
if(Branch[i].Kt==0)//非变压器支路
{
G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2;
B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2;
G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1];
B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1];
}
else//变压器支路
{
G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2/Branch[i].Kt;
B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2/Branch[i].Kt;
G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1];
B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1];
}
}
//计算对角元素
for(i=0;ifor(j=0;j{
Z2=Branch[j].R*Branch[j].R+Branch[j].X*Branch[j].X;
if(Branch[j].Kt==0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i))//非变压器支路
{
G[i][i]=G[i][i]+Branch[j].R/Z2;
B[i][i]=B[i][i]-Branch[j].X/Z2;
}
elseif(Branch[j].Kt!
=0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i))//变压器支路
{
G[i][i]=G[i][i]+Branch[j].R/Z2/Branch[j].Kt;
B[i][i]=B[i][i]-Branch[j].X/Z2/Branch[j].Kt;
}
}
//将对地电纳加入到对角元素中
for(i=0;i{
if(Branch[i].Kt==0)//非变压器支路
{
B[Branch[i].Nl-1][Branch[i].Nl-1]+=Branch[i].Bn;
B[Branch[i].Nr-1][Branch[i].Nr-1]+=Branch[i].Bn;
}
else//变压器支路
{
B[Branch[i].Nl-1][Branch[i].Nl-1]-=(1-Branch[i].Kt)/Branch[i].Kt/Branch[i].Kt/Branch[i].X;
B[Branch[i].Nr-1][Branch[i].Nr-1]-=(Branch[i].Kt-1)/Branch[i].Kt/Branch[i].X;
}
}
//将并联电容加入到对角元素中
for(i=0;iB[i][i]=B[i][i]+Node[i].Bc;
}
//*************************************************************************************************
//*****************************************求deltaP,deltaQ*****************************************
voidDeltaS()//计算功率不平衡量
{
inti,j;
//为中间变量申请空间
mid1=(double*)malloc(sizeof(double)*n);
mid2=(double*)malloc(sizeof(double)*n);
//为功率不平衡量申请空间
dS=(double*)malloc(sizeof(double)*2*(n-1));
//求功率不平衡量
for(i=0;i{
//初始化中间变量
mid1[i]=0;
mid2[i]=0;
for(j=0;j{
mid1[i]=mid1[i]+G[i][j]*Node[j].e-B[i][j]*Node[j].f;
mid2[i]=mid2[i]+G[i][j]*Node[j].f+B[i][j]*Node[j].e;
}
dS[2*i]=Node[i].Ps-Node[i].Pd-(Node[i].e*mid1[i]+Node[i].f*mid2[i]);
if(idS[2*i+1]=Node[i].Qs-Node[i].Qd-(Node[i].f*mid1[i]-Node[i].e*mid2[i]);
else
dS[2*i+1]=Us[i]*Us[i]-(Node[i].e*Node[i].e+Node[i].f*Node[i].f);
}
error=0;
for(i=0;i<2*(n-1);i++)
{
if(dS[i]<0&&error<-dS[i])
error=-dS[i];
elseif(dS[i]>0&&errorerror=dS[i];
}
}
//*************************************************************************************************
//*********************************************雅克比矩阵******************************************
voidFormJacob()
{
inti,j;
//为雅克比行列式申请空间
Jacob=(double**)malloc(sizeof(double*)*2*(n-1));
for(i=0;i<2*(n-1);i++)
Jacob[i]=(double*)malloc(sizeof(double)*2*(n-1));
//初始化雅克比行列式
for(i=0;i<2*(n-1);i++)
for(j=0;j<2*(n-1);j++)
Jacob[i][j]=0;
for(j=0;j{
//求H,N
for(i=0;i{
if(i!
=j)
{
Jacob[2*i][2*j]=B[i][j]*Node[i].e-G[i][j]*Node[i].f;
Jacob[2*i][2*j+1]=-G[i][j]*Node[i].e-B[i][j]*Node[i].f;
}
else
{
Jacob[2*i][2*i]=B[i][i]*Node[i].e-G[i][i]*Node[i].f-mid2[i];
Jacob[2*i][2*i+1]=-G[i][i]*Node[i].e-B[i][i]*Node[i].f-mid1[i];
}
}
//求J,L
for(i=0;i{
if(i!
=j)
{
Jacob[2*i+1][2*j]=G[i][j]*Node[i].e+B[i][j]*Node[i].f;
Jacob[2*i+1][2*j+1]=B[i][j]*Node[i].e-G[i][j]*Node[i].f;
}
else
{
Jacob[2*i+1][2*i]=G[i][i]*Node[i].e+B[i][i]*Node[i].f-mid1[i];
Jacob[2*i+1][2*i+1]=B[i][i]*Node[i].e-G[i][i]*Node[i].f+mid2[i];
}
}
//求R,S
for(i=nPQ;i{
if(i==j)
{
Jacob[2*i+1][2*i]=-2*Node[i].f;
Jacob[2*i+1][2*i+1]=-2*Node[i].e;
}
}
}
}
//*************************************************************************************************
//********************************************雅克比矩阵求逆***************************************
voidInvJac()
{
inti,j,k;
doubletemp;//中间变量
//为雅克比矩阵的逆申请空间
invJac=(double**)malloc(sizeof(double*)*2*(n-1));
for(i=0;i<2*(n-1);i++)
invJac[i]=(double*)malloc(sizeof(double)*2*(n-1));
//求逆
for(i=0;i<2*(n-1);i++)
for(j=0;j<2*(n-1);j++)
{
if(i!
=j)
invJac[i][j]=0;
else
invJac[i][j]=1;
}
for(i=0;i<2*(n-1);i++)
{
for(j=0;j<2*(n-1);j++)
{
if(i!
=j)
{
temp=Jacob[j][i]/Jacob[i][i];
for(k=0;k<2*(n-1);k++)
{
Jacob[j][k]-=Jacob[i][k]*temp;
invJac[j][k]-=invJac[i][k]*temp;
}
}
}
}
for(i=0;i<2*(n-1);i++)
if(Jacob[i][i]!
=1)
{
temp=Jacob[i][i];
for(j=0;j<2*(n-1);j++)
invJac[i][j]=invJac[i][j]/temp;
}
}
//*************************************************************************************************
//*********************************************电压修正********************************************
voidUpdateU()
{
voidInvJac();//求雅克比矩阵的逆
inti,j;
dfe=(double*)malloc(sizeof(double)2*(n-1));
InvJac();
for(i=0;i<2*(n-1);i++)
{
dfe[i]=0;
for(j=0;j<2*(n-1);j++)
dfe[i]-=invJac[i][j]*dS[j];
}
for(i=0;i{
Node[i].e+=dfe[2*i+1];
Node[i].f+=dfe[2*i];
}
}
voidCalculatePQ()
{
inti,j;
inttN,tType;
doublete,tf,tPd,tQd,tPs,tQs,tBc;//用于重新排列节点信息的临时变量
//计算平衡节点功率
mid1[n-1]=0;
mid2[n-1]=0;
for(j=0;j{
mid1[n-1]=mid1[n-1]+G[n-1][j]*Node[j].e-B[n-1][j]*Node[j].f;
mid2[n-1]=mid2[n-1]+G[n-1][j]*Node[j].f+B[n-1][j]*Node[j].e;
}
Node[n-1].Ps=Node[n-1].e*mid1[n-1];
Node[n-1].Qs=-Node[n-1].e*mid2[n-1];
//计算PV节点的Q