牛顿拉夫逊迭代法极坐标潮流计算C语言程序.docx
- 文档编号:3983312
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:27
- 大小:21.79KB
牛顿拉夫逊迭代法极坐标潮流计算C语言程序.docx
《牛顿拉夫逊迭代法极坐标潮流计算C语言程序.docx》由会员分享,可在线阅读,更多相关《牛顿拉夫逊迭代法极坐标潮流计算C语言程序.docx(27页珍藏版)》请在冰豆网上搜索。
牛顿拉夫逊迭代法极坐标潮流计算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,
);
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);
fprintf(fp2,\
==============================================================================\n);
fprintf(fp2,\
节点导纳矩阵为:
);
;.
.
for(i=1;i<=n;i++)
{
fprintf(fp2,\
);
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);
}
/*-------------------------------牛顿——拉夫逊-------------------------------*/
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,\
不平衡量为:
);/*输出不平衡量*/
for(i=1;i<=pq;i++)
{
t=jd[i].num;
fprintf(fp2,\
dp[%d]=%f,i,D[2*t-1]);
fprintf(fp2,\
dq[%d]=%f,i,D[2*t]);
}
;.
.
for(i=pq+1;i<=n-1;i++)
{
t=jd[i].num;
fprintf(fp2,\
dp[%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雅克比矩阵为: ); for(i=1;i<=(2*pq+pv);i++) { fprintf(fp2,\ ); 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};//定义数组
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 牛顿 拉夫逊 迭代法 坐标 潮流 计算 语言 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)