牛顿拉夫逊极坐标形式的潮流上机c源程序Word文件下载.docx
《牛顿拉夫逊极坐标形式的潮流上机c源程序Word文件下载.docx》由会员分享，可在线阅读，更多相关《牛顿拉夫逊极坐标形式的潮流上机c源程序Word文件下载.docx（22页珍藏版）》请在冰豆网上搜索。
/*节点有功、无功功率，功率模值，电压模值,阻抗角
牛顿--拉夫逊中功率不平衡量、电压修正量*/
}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节点*/
{
%lf,%lf,%lf\n"
jd[i].q=-1.567;
if（h==3）/*类型h=3是平衡节点*/
%lf,%lf\n"
}
=m;
i++）/*输入支路阻抗*/
%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华北电力大学电力1104谭程凯201105030119\n"
**********原始数据*********\n"
================================================================================\n"
节点数:
%d支路数:
%dPQ节点数:
%dPV节点数:
%d精度:
%f\n"
n,m,pq,pv,eps）;
------------------------------------------------------------------------------\n"
=pq;
i++）
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;
=pq+pv;
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）;
平衡节点:
节点%de[%d]=%ff[%d]=%f\n"
jd[n].num,jd[n].num,jd[n].U,jd[n].num,jd[n].zkj）;
-------------------------------------------------------------------------------\n"
支路%d:
相关节点:
%d,%d非标准变比:
kx=%fR=%fX=%f\n"
i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x）;
==============================================================================\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（）
i++）/*节点导纳矩阵元素附初值*/
for（j=1;
j<
j++）
G[i][j]=B[i][j]=0;
for（i=1;
i++）/*节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中（考虑非标准变比）*/
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）;
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）
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<
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;
G[i][j]-=c2/zl[k].kx;
B[i][j]-=d2/zl[k].kx;
/*--------------------------输出节点导纳矩阵------------------------------*/
\n\n*********潮流计算过程*********\n"
\n==============================================================================\n"
\n节点导纳矩阵为:
"
\n"
k++）
fprintf（fp2,"
%f"
G[i][k]）;
if（B[i][k]>
=0）
fprintf（fp2,"
+j"
%f"
B[i][k]）;
else
B[i][k]=-B[i][k];
-j"
B[i][k]=-B[i][k];
\n------------------------------------------------------------------------------\n"
/*-------------------------------牛顿——拉夫逊-------------------------------*/
voidCalculate_Unbalanced_Para（）
if（jd[i].ty==1）/*计算PQ节点不平衡量*/
t=jd[i].num;
cc[t]=dd[t]=0;
for（j=1;
{
for（a=1;
a<
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节点不平衡量*/
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].q=jd[i].U*dd[t];
i++）/*形成不平衡量矩阵D[M]*/
{
D[2*i-1]=jd[i].dp;
D[2*i]=jd[i].dq;
=n-1;
D[pq+i]=jd[i].dp;
\n不平衡量为:
/*输出不平衡量*/
t=jd[i].num;
\ndp[%d]=%f"
i,D[2*t-1]）;
\ndq[%d]=%f"
i,D[2*t]）;
i,D[pq+t]）;
voidForm_Jacobi_Matric（）/*形成雅克比矩阵*/
{
i++）/*形成pq节点子阵*/
n;
{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;
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））;
ykb[2*i][pq+j]=-Ui*Uj*（G[i1][j1]*cos（zi-zj）+B[i1][j1]*sin（zi-zj））;
i++）/*形成pv节点子阵*/
for（j=1;
inti1=jd[i].num;
if（j<
=pq时的H、N*/
ykb[pq+i][2*j-1]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;
ykb[pq+i][2*j]=Ui*Uj*（G[i1][j1]*cos（zi-zj）+B[i1][j1]*sin（zi-zj））;
else/*求j>
pq时的H*/
{
=j时的H*/
ykb[pq+i][pq+j]=Ui*Uj*（G[i1][j1]*sin（zi-zj）-B[i1][j1]*cos（zi-zj））;
else/*求i=j时的H*/
aa[i]=0;
intk1=jd[k].num;
H=0;
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）））;
/*------------------------------输出雅克比矩阵--------------------------------*/
\n\n雅克比矩阵为:
"
=（2*pq+pv）;
i++）
=2*pq+pv;
%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};
//定义
升级会员 