牛顿拉夫逊迭代法极坐标潮流计算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;
intn;
n=2*pq+pv;
f
升级会员 