牛顿拉夫逊迭代法极坐标潮流计算C语言程序综述.docx

牛顿拉夫逊迭代法极坐标潮流计算C语言程序综述.docx

《牛顿拉夫逊迭代法极坐标潮流计算C语言程序综述.docx》由会员分享，可在线阅读，更多相关《牛顿拉夫逊迭代法极坐标潮流计算C语言程序综述.docx（26页珍藏版）》请在冰豆网上搜索。

牛顿拉夫逊迭代法极坐标潮流计算C语言程序综述

/*利用牛顿-拉夫逊迭代法（极坐标形式），计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。

所有参数应归算至标幺值下。

/*可计算最大节点数为100，可计算PQ,PV,平衡节点*/

/*可计算非标准变比和平行支路*/

#include

#include

#include

#defineM100/*最大矩阵阶数*/

#defineNl100/*迭代次数*/

inti,j,k,a,b,c;/*循环控制变量*/

intt,l;

doubleP,Q,H,J;/*中间变量*/

intn,/*节点数*/

m,/*支路数*/

pq,/*PQ节点数*/

pv;/*PV节点数*/

doubleeps;/*迭代精度*/

doubleaa[M],bb[M],cc[M],dd[M],max,rr,tt;/*中间变量*/

doublemo,c1,d1,c2,d2;/*复数运算函数的返回值*/

doubleG[M][M],B[M][M],Y[M][M];/*节点导纳矩阵中的实部、虚部及其模方值*/

doubleykb[M][M],D[M],d[M],dU[M];/*雅克比矩阵、不平衡量矩阵*/

structjd/*节点结构体*/

{intnum,ty;/*num为节点号，ty为节点类型*/

doublep,q,S,U,zkj,dp,dq,du,dj;/*节点有功、无功功率，功率模值，电压模值,阻抗角

牛顿--拉夫逊中功率不平衡量、电压修正量*/

}jd[M];

structzl/*支路结构体*/

{intnumb;/*numb为支路号*/

intp1,p2;/*支路的两个节点*/

doublekx;/*非标准变比*/

doubler,x;/*支路的电阻与电抗*/

}zl[M];

FILE*fp1,*fp2;

voiddata（）/*读取数据*/

{

inth,number;

fp1=fopen（"input.txt","r"）;

fscanf（fp1,"%d,%d,%d,%d,%lf\n",&n,&m,&pq,&pv,&eps）;/*输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/

for（i=1;i<=n;i++）/*输入节点编号、类型、输入功率和电压初值*/

{

fscanf（fp1,"%d,%d",&number,&h）;

if（h==1）/*类型h=1是PQ节点*/

{

fscanf（fp1,"%lf,%lf,%lf,%lf\n",&jd[i].p,&jd[i].q,&jd[i].U,&jd[i].zkj）;

jd[i].num=number;

jd[i].ty=h;

}

if（h==2）/*类型h=2是pv节点*/

{

fscanf（fp1,",%lf,%lf,%lf\n",&jd[i].p,&jd[i].U,&jd[i].zkj）;

jd[i].num=number;

jd[i].ty=h;

jd[i].q=-1.567;

}

if（h==3）/*类型h=3是平衡节点*/

{

fscanf（fp1,",%lf,%lf\n",&jd[i].U,&jd[i].zkj）;

jd[i].num=number;

jd[i].ty=h;

}

}

for（i=1;i<=m;i++）/*输入支路阻抗*/

fscanf（fp1,"%d,%lf,%d,%d,%lf,%lf\n",&zl[i].numb,&zl[i].kx,&zl[i].p1,&zl[i].p2,&zl[i].r,&zl[i].x）;

fclose（fp1）;

if（（fp2=fopen（"output.txt","w"））==NULL）

{

printf（"cannotopenfile!

\n"）;

exit（0）;

}

fprintf（fp2,"电力系统潮流计算\n"）;

fprintf（fp2,"**********原始数据*********\n"）;

fprintf（fp2,"================================================================================\n"）;

fprintf（fp2,"节点数:

%d支路数:

%dPQ节点数:

%dPV节点数:

%d精度:

%f\n",

n,m,pq,pv,eps）;

fprintf（fp2,"------------------------------------------------------------------------------\n"）;

for（i=1;i<=pq;i++）

fprintf（fp2,"PQ节点:

节点%dP[%d]=%fQ[%d]=%f\n",

jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].q）;

for（i=pq+1;i<=pq+pv;i++）

fprintf（fp2,"PV节点:

节点%dP[%d]=%fU[%d]=%f初值Q[%d]=%f\n",

jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].U,jd[i].num,jd[i].q）;

fprintf（fp2,"平衡节点:

节点%de[%d]=%ff[%d]=%f\n",

jd[n].num,jd[n].num,jd[n].U,jd[n].num,jd[n].zkj）;

fprintf（fp2,"-------------------------------------------------------------------------------\n"）;

for（i=1;i<=m;i++）

fprintf（fp2,"支路%d:

相关节点:

%d,%d非标准变比:

kx=%fR=%fX=%f\n",

i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x）;

fprintf（fp2,"==============================================================================\n"）;

}

/*------------------------------------复数运算函数--------------------------------------*/

doublemozhi（doublea0,doubleb0）/*复数求模值函数*/

{mo=sqrt（a0*a0+b0*b0）;

returnmo;

}

doubleji（doublea1,doubleb1,doublea2,doubleb2）/*复数求积函数a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/

{a1=a1*cos（b1）;

b1=a1*tan（b1）;

c1=a1*a2-b1*b2;

d1=a1*b2+a2*b1;

returnc1;

returnd1;

}

doubleshang（doublea3,doubleb3,doublea4,doubleb4）/*复数除法求商函数*/

{c2=（a3*a4+b3*b4）/（a4*a4+b4*b4）;

d2=（a4*b3-a3*b4）/（a4*a4+b4*b4）;

returnc2;

returnd2;

}

/*--------------------------------计算节点导纳矩阵----------------------------------*/

voidForm_Y（）

{

for（i=1;i<=n;i++）/*节点导纳矩阵元素附初值*/

for（j=1;j<=n;j++）

G[i][j]=B[i][j]=0;

for（i=1;i<=n;i++）/*节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中（考虑非标准变比）*/

for（j=1;j<=m;j++）

if（zl[j].p1==i）

{

if（zl[j].kx==1）

{

mozhi（zl[j].r,zl[j].x）;

if（mo==0）continue;

shang（1,0,zl[j].r,zl[j].x）;

G[i][i]+=c2;

B[i][i]+=d2;

}

else

{

mozhi（zl[j].r,zl[j].x）;

if（mo==0）continue;

shang（1,0,zl[j].r,zl[j].x）;

G[i][i]+=c2/zl[j].kx+c2*（1-zl[j].kx）/（zl[j].kx*zl[j].kx）;

B[i][i]+=d2/zl[j].kx+d2*（1-zl[j].kx）/（zl[j].kx*zl[j].kx）;

}

}

elseif（zl[j].p2==i）

{

if（zl[j].kx==1）

{

mozhi（zl[j].r,zl[j].x）;

if（mo==0）continue;

shang（1,0,zl[j].r,zl[j].x）;

G[i][i]+=c2;

B[i][i]+=d2;

}

else

{

mozhi（zl[j].r,zl[j].x）;

if（mo==0）continue;

shang（1,0,zl[j].r,zl[j].x）;

G[i][i]+=c2/zl[j].kx+c2*（zl[j].kx-1）/zl[j].kx;

B[i][i]+=d2/zl[j].kx+d2*（zl[j].kx-1）/zl[j].kx;

}

}

for（k=1;k<=m;k++）/*节点导纳矩阵非主对角线上（考虑非标准变比）的元素*/

if（zl[k].kx==1）

{

i=zl[k].p1;

j=zl[k].p2;

mozhi（zl[k].r,zl[k].x）;

if（mo==0）continue;

shang（1,0,zl[k].r,zl[k].x）;

G[i][j]-=c2;

B[i][j]-=d2;

G[j][i]=G[i][j];

B[j][i]=B[i][j];

}

else

{

i=zl[k].p1;

j=zl[k].p2;

mozhi（zl[k].r,zl[k].x）;

if（mo==0）continue;

shang（1,0,zl[k].r,zl[k].x）;

G[i][j]-=c2/zl[k].kx;

B[i][j]-=d2/zl[k].kx;

G[j][i]=G[i][j];

B[j][i]=B[i][j];

}

/*--------------------------输出节点导纳矩阵------------------------------*/

fprintf（fp2,"\n\n*********潮流计算过程*********\n"）;

fprintf（fp2,"\n==============================================================================\n"）;

fprintf（fp2,"\n节点导纳矩阵为:

"）;

for（i=1;i<=n;i++）

{

fprintf（fp2,"\n"）;

for（k=1;k<=n;k++）

{

fprintf（fp2,"%f",G[i][k]）;

if（B[i][k]>=0）

{

fprintf（fp2,"+j"）;

fprintf（fp2,"%f",B[i][k]）;

}

else

{

B[i][k]=-B[i][k];

fprintf（fp2,"-j"）;

fprintf（fp2,"%f",B[i][k]）;

B[i][k]=-B[i][k];

}

}

}

fprintf（fp2,"\n------------------------------------------------------------------------------\n"）;

}

/*-------------------------------牛顿——拉夫逊-------------------------------*/

voidCalculate_Unbalanced_Para（）

{

for（i=1;i<=n;i++）

{

if（jd[i].ty==1）/*计算PQ节点不平衡量*/

{

t=jd[i].num;

cc[t]=dd[t]=0;

for（j=1;j<=n;j++）

{

for（a=1;a<=n;a++）

{

if（jd[a].num==j）

break;

}

P=Q=0;

P=jd[a].U*（G[t][j]*cos（jd[i].zkj-jd[a].zkj）+B[t][j]*sin（jd[i].zkj-jd[a].zkj））;

Q=jd[a].U*（G[t][j]*sin（jd[i].zkj-jd[a].zkj）-B[t][j]*cos（jd[i].zkj-jd[a].zkj））;

cc[t]+=P;

dd[t]+=Q;

}

jd[i].dp=jd[i].p-jd[i].U*cc[t];

jd[i].dq=jd[i].q-jd[i].U*dd[t];

}

if（jd[i].ty==2）/*计算PV节点不平衡量*/

{

t=jd[i].num;

cc[t]=dd[t]=0;

for（j=1;j<=n;j++）

{

for（a=1;a<=n;a++）

{

if（jd[a].num==j）

break;

}

P=Q=0;

P=jd[a].U*（G[t][j]*cos（jd[i].zkj-jd[a].zkj）+B[t][j]*sin（jd[i].zkj-jd[a].zkj））;

Q=jd[a].U*（G[t][j]*sin（jd[i].zkj-jd[a].zkj）-B[t][j]*cos（jd[i].zkj-jd[a].zkj））;

cc[t]+=P;

dd[t]+=Q;

}

jd[i].dp=jd[i].p-jd[i].U*cc[t];

jd[i].q=jd[i].U*dd[t];

}

}

for（i=1;i<=pq;i++）/*形成不平衡量矩阵D[M]*/

{

D[2*i-1]=jd[i].dp;

D[2*i]=jd[i].dq;

}

for（i=pq+1;i<=n-1;i++）

{

D[pq+i]=jd[i].dp;

}

fprintf（fp2,"\n不平衡量为:

"）;/*输出不平衡量*/

for（i=1;i<=pq;i++）

{

t=jd[i].num;

fprintf（fp2,"\ndp[%d]=%f",i,D[2*t-1]）;

fprintf（fp2,"\ndq[%d]=%f",i,D[2*t]）;

}

for（i=pq+1;i<=n-1;i++）

{

t=jd[i].num;

fprintf（fp2,"\ndp[%d]=%f",i,D[pq+t]）;

}

}

voidForm_Jacobi_Matric（）/*形成雅克比矩阵*/

{

for（i=1;i<=pq;i++）/*形成pq节点子阵*/

for（j=1;j

{inti1=jd[i].num;

intj1=jd[j].num;

doubleUi=jd[i].U;

doubleUj=jd[j].U;

doublezi=jd[i].zkj;

doublezj=jd[j].zkj;

if（j<=pq）/*求j<=pq时的H、N、J、L*/

{

if（i!

=j）/*求i!

=j时的H、N、J、L*/

{

ykb[2*i-1][2*j-1]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;/*H*/

ykb[2*i-1][2*j]=Ui*Uj*（G[i1][j1]*cos（zi-zj）+B[i1][j1]*sin（zi-zj））;/*N*/

ykb[2*i][2*j-1]=-ykb[2*i-1][2*j];/*J*/

ykb[2*i][2*j]=ykb[2*i-1][2*j-1];/*L*/

}

else/*求i=j时的H、N、J、L*/

{

aa[i]=0;bb[i]=0;

for（k=1;k<=n;k++）

{

intk1=jd[k].num;

H=J=0;

H=jd[k].U*（G[i1][k1]*sin（jd[i].zkj-jd[k].zkj）-B[i1][k1]*cos（jd[i].zkj-jd[k].zkj））;

J=jd[k].U*（G[i1][k1]*cos（jd[i].zkj-jd[k].zkj）+B[i1][k1]*sin（jd[i].zkj-jd[k].zkj））;

aa[i]=aa[i]+H;

bb[i]=bb[i]+J;

}

ykb[2*i-1][2*j-1]=-Ui*（aa[i]-Ui*（G[i1][i1]*sin（jd[i].zkj-jd[i].zkj）-B[i1][i1]*cos（jd[i].zkj-jd[i].zkj）））;/*H*/

ykb[2*i][2*j-1]=Ui*（bb[i]-Ui*（G[i1][i1]*cos（jd[i].zkj-jd[i].zkj）+B[i1][i1]*sin（jd[i].zkj-jd[i].zkj）））;/*J*/

ykb[2*i-1][2*j]=ykb[2*i][2*j-1]+2*Ui*Ui*G[i1][i1];/*N*/

ykb[2*i][2*j]=-ykb[2*i-1][2*j-1]-2*Ui*Ui*B[i1][i1];/*L*/

}

}

else/*求j>pq时的H、J*/

{ykb[2*i-1][pq+j]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;/*H*/

ykb[2*i][pq+j]=-Ui*Uj*（G[i1][j1]*cos（zi-zj）+B[i1][j1]*sin（zi-zj））;/*J*/

}

}

for（i=pq+1;i<=n-1;i++）/*形成pv节点子阵*/

for（j=1;j

{

inti1=jd[i].num;

intj1=jd[j].num;

doubleUi=jd[i].U;

doubleUj=jd[j].U;

doublezi=jd[i].zkj;

doublezj=jd[j].zkj;

if（j<=pq）/*求j<=pq时的H、N*/

{

ykb[pq+i][2*j-1]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;/*H*/

ykb[pq+i][2*j]=Ui*Uj*（G[i1][j1]*cos（zi-zj）+B[i1][j1]*sin（zi-zj））;/*N*/

}

else/*求j>pq时的H*/

{

if（i!

=j）/*求i!

=j时的H*/

ykb[pq+i][pq+j]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;/*H*/

else/*求i=j时的H*/

{

aa[i]=0;

for（k=1;k<=n;k++）

{

intk1=jd[k].num;

H=0;

H=jd[k].U*（G[i1][k1]*sin（jd[i].zkj-jd[k].zkj）-B[i1][k1]*cos（jd[i].zkj-jd[k].zkj））;

aa[i]=aa[i]+H;

}

ykb[pq+i][pq+j]=-Ui*（aa[i]-Ui*（G[i1][i1]*sin（jd[i].zkj-jd[i].zkj）-B[i1][i1]*cos（jd[i].zkj-jd[i].zkj）））;/*H*/

}

}

}

/*------------------------------输出雅克比矩阵--------------------------------*/

fprintf（fp2,"\n\n雅克比矩阵为:

"）;

for（i=1;i<=（2*pq+pv）;i++）

{

fprintf（fp2,"\n"）;

for（j=1;j<=2*pq+pv;j++）

{

fprintf（fp2,"%f",ykb[i][j]）;

}

}

}

voidSolve_Equations（）/*求解修正方程组（LU分解法）*/

{

doublel[Nl][Nl]={0};//定义L矩阵

doubleu[Nl][Nl]={0};//定义U矩阵

doubley[Nl]={0};//定义数组Y

doublex[Nl]={0};//定义数组X

doublea[Nl][Nl]={0};//定义系数矩阵

doubleb[Nl]={0};//定义右端项

doublesum=0;

inti,j,k,s