电力系统分析潮流计算程序.docx

电力系统分析潮流计算程序.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;i

fscanf（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;j

for（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;i

Us[i]=Node[i].e;

//读取支路参数

fscanf（fp,"%d",&nbr）;

printf（"支路个数为：

%d\n",nbr）;

//为支路参数申请空间

Branch=（structBranchType*）malloc（sizeof（structBranchType）*nbr）;//读入的支路参数结构体

//读入支路参数

for（i=0;i

fscanf（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;i

for（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;i

for（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;i

B[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（i

dS[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&&error

error=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