潮流c.docx
《潮流c.docx》由会员分享，可在线阅读，更多相关《潮流c.docx（19页珍藏版）》请在冰豆网上搜索。
潮流c
潮流c
//潮流计算上机程序
//设只有一个平衡节点
//平行支路与接地支路只出现一种情况
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#definemax100//最大迭代次数
#defineM20//最多节点数
inti,j,kk,h;//定义循环控制变量
intk=1;//迭代次数
intn=0;//定义节点个数
intm=0;//定义支路个数
intpq=0;//定义PQ节点个数
intpv=0;//定义PV节点个数
floatY[M][M][2];//定义节点导纳矩阵
doubleeps=0.0001;//定义精度
floatD[M];//定义修正量矩阵
structJD//定义节点结构体
{
intnum,s;//s为节点类型
floatp,q,e,f,v,dp,dq,de,df,dv;//de,df,dp,dq,dv为各不平衡量
};
structZL//定义支路结构体
{
intnum,n1,n2;//n1,n2为与该支路相关的节点
floatr,x,zlp1,zlq1,zlp2,zlq2,dzlp,dzlq;//zlp1,zlq1,zlp2,zlq2,dzlp,dzlq分别为支路流通
的功率和支路功率损耗
};
/*structDN
{
intnum,n1,n2;
floatg,b;
};
structDNdn[M];*/
structJDjd[M];//定义节点数组
structZLzl[M];//定义支路数组
floatyakebi[2*（M-1）][2*（M-1）+1];//雅克比矩阵定义+1为为后面求解准备
FILE*fp1,*fp2;//定义文件指针
voidRead_Data（）
{
fp1=fopen（"input.txt","r"）;
if（fp1==NULL）
第1/13页
{
printf（"请先建立input.txt文件!
\n"）;
exit（0）;
}
fscanf（fp1,"%d,%d,%d,%d,%lf\n",&n,&m,&pq,&pv,&eps）;//从文件中获得节点个数、
支路个数、PQ节点个数、PV节点个数及精度
intii=0,jj=pq;
for（i=0;i<n;i++）
{
intnum1=0,s1=0;
fscanf（fp1,"%d,%d,",&num1,&s1）;
if（s1==1）
{
fscanf（fp1,"%f,%f\n",&jd[ii].p,&jd[ii].q）;//输入pq节点信息：
节点编号，节点
类型1，P+jq
jd[ii].num=num1;//把pq节点放在一起
jd[ii].s=s1;
ii++;
}
elseif（s1==2）
{
fscanf（fp1,"%f,%f\n",&jd[jj].p,&jd[jj].v）;//输入pv节点信息：
节点编号，节点
类型2，p，v
jd[jj].num=num1;//把pv节点放在一起
jd[jj].s=s1;
jj++;
}
elseif（s1==3）
{
fscanf（fp1,"%f,%f\n",&jd[n-1].e,&jd[n-1].f）;//输入平衡节点信息：
节点编号，节
点类型3，e+jf
jd[n-1].num=num1;//平衡节点放在最后
jd[n-1].s=s1;
}
}
for（i=0;i<m;i++）//从文件中获得支路阻抗和相关节点信息：
支
路编号，支路相关节点n1，n2，支路阻抗r+jx
fscanf（fp1,"%d,%d,%d,%f,%f\n",&zl[i].num,&zl[i].n1,&zl[i].n2,&zl[i].r,&zl[i].x）;//对对地
支路约定n1为所接节点，n2=0
fclose（fp1）;
fp2=fopen（"output.txt","w"）;
fprintf（fp2,"\n潮流上机\n****电力0803********\n"）;
第2/13页
fprintf（fp2,"\n原始数据\n"）;
fprintf（fp2,"\n***********************************************************************************\n"）;
fprintf（fp2,"节点数:
%d支路数:
%dPQ节点数:
%dPV节点数:
%d精度:
%f\n",n,m,pq,pv,eps）;
fprintf（fp2,"\n===================================================================================\n"）;
for（i=0;i<n;i++）
{
if（jd[i].s==1）fprintf（fp2,"节点%dPQ节点P=%fQ=%f\n",jd[i].num,jd[i].p,jd[i].q）;
if
（jd[i].s==2）fprintf（fp2,"节点%dPV节点P=%fV=%f\n",jd[i].num,jd[i].p,jd[i].v）;
if（jd[i].s==3）fprintf（fp2,"节点%d平衡节点e=%ff=%f\n",jd[i].num,jd[i].e,jd[i].f）;
}
fprintf（fp2,"\n===================================================================================\n"）;
for（i=0;i<m;i++）
fprintf（fp2,"支路%d相关节点:
%d%dR=%fX=%f\n",zl[i].num,zl[i].n1,zl[i].n2,zl[i].r,zl[i].x）;
fprintf（fp2,"\n==================================================================================\n"）;
}
//形成节点导纳矩阵考虑对地支路及并行支路
voidForm_Y（）
{
for（i=0;i<m;i++）//求支路的导纳
{
floatrr=zl[i].r;
floatxx=zl[i].x;
zl[i].r=zl[i].r/（rr*rr+xx*xx）;//把支路阻抗替代为支路导纳
zl[i].x=-zl[i].x/（rr*rr+xx*xx）;//把支路阻抗替代为支路导纳
}
for（i=0;i<M;i++）//节点导纳矩阵初始化防止因节点间没有支路而数组默认值对结果的影响
for（j=0;j<M;j++）
{
Y[i][j][0]=0.0;
第3/13页
Y[i][j][1]=0.0;
}
for（i=0;i<m;i++）//形成节点导纳矩阵
{
intn3=zl[i].n1-1;
intn4=zl[i].n2-1;
if（zl[i].n1==0）//对地支路节点导纳矩阵的形成
{
Y[n4][n4][0]=Y[n4][n4][0]+zl[i].r;//对角线元素加上对地导纳实部Y[n4][n4][1]=Y[n4][n4][1]+zl[i].x;//对角线元素加上对地导纳虚部}elseif（zl[i].n2==0）//对地支路节点导纳矩阵的形成{Y[n3][n3][0]=Y[n3][n3][0]+zl[i].r;//对角线元素加上对地导纳实部Y[n3][n3][1]=Y[n3][n3][1]+zl[i].x;//对角线元素加上对地导纳虚部}
else//非对地支路节点导纳矩阵的形成
{
Y[n3][n3][0]=Y[n3][n3][0]+zl[i].r;//对角线元素实部
Y[n3][n3][1]=Y[n3][n3][1]+zl[i].x;//对角线元素虚部
Y[n4][n4][0]=Y[n4][n4][0]+zl[i].r;//对角线元素实部
Y[n4][n4][1]=Y[n4][n4][1]+zl[i].x;//对角线元素虚部
Y[n3][n4][0]=-Y[n3][n4][0]-zl[i].r;//非对角线元素实部
Y[n3][n4][1]=-Y[n3][n4][1]-zl[i].x;//非对角线元素虚部
Y[n4][n3][0]=-Y[n4][n3][0]-zl[i].r;//非对角线元素实部
Y[n4][n3][1]=-Y[n4][n3][1]-zl[i].x;//非对角线元素虚部
}
}
fprintf（fp2,"\n计算过程及结果\n"）;
fprintf（fp2,"\n***********************************************************************************\n"）;
fprintf（fp2,"\n节点导纳矩阵为：
\n"）;
for（i=0;i<n;i++）
{
for（j=0;j<n;j++）
{
fprintf（fp2,"%10f+j%-10.6f",Y[i][j][0],Y[i][j][1]）;
}
fprintf（fp2,"\n"）;
}
fprintf（fp2,"\n==================================================================================\n"）;
第4/13页
}
voidfuchuzhi（）//赋电压迭代初值
{
for（i=0;i<n;i++）
{
if（jd[i].s!
=3）
{
jd[i].e=1;
jd[i].f=0;
}
}
}
//计算不平衡功率、电压
voidCalculate_Unbalanced_Value（）
{
fprintf（fp2,"\n不平衡量为:
\n\n"）;
for（i=0;i<n;i++）
{
if（jd[i].s==1）//判断为pq节点
{
floatsum1=0,sum2=0;
for（j=0;j<n;j++）
{
sum1+=（jd[i].e*（Y[i][j][0]*jd[j].e-Y[i][j][1]*jd[j].f）+jd[i].f*（Y[i][j][0]*jd[j].f+Y[i][j][1]*jd[j].e））;
sum2+=（jd[i].f*（Y[i][j][0]*jd[j
].e-Y[i][j][1]*jd[j].f）-jd[i].e*（Y[i][j][0]*jd[j].f+Y[i][j][1]*jd[j].e））;}
jd[i].dp=jd[i].p-sum1;//计算有功不平衡量
jd[i].dq=jd[i].q-sum2;//计算武功不平衡量
fprintf（fp2,"节点%d的不平衡量为:
dp=%-10fdq=%-10f\n",jd[i].num,jd[i].dp,jd[i].dq）;
}
elseif（jd[i].s==2）//判断为pv节点
{
floatsum1=0;
for（j=0;j<n;j++）
{
sum1+=（jd[i].e*（Y[i][j][0]*jd[j].e-Y[i][j][1]*jd[j].f）+jd[i].f*（Y[i][j][0]*jd[j].f+Y[i][j][1]*jd[j].e））;}jd[i].dp=jd[i].p-sum1;//计算有功不平衡量
jd[i].dv=（jd[i].v*jd[i].v-（jd[i].e*jd[i].e+jd[i].f*jd[i].f））;//计算电压不平衡量
fprintf（fp2,"节点%d的不平衡量为:
dp=%-10fdu^2=%-10f\n",jd[i].num,jd[i].dp,jd[i].dv）;
}
}
fprintf（fp2,"\n==================================================================================\n"）;
}
//形成雅克比矩阵
voidForm_Jacobi_Matric（）
{
intii=0;//临时变量
for（i=0;i<n;i++）
{
intjj=0;
if（jd[i].s==3）//判断为平衡节点
{
continue;//结束本次循环进入下次循环
}
elseif（jd[i].s==1）//判断为PQ节点
{
for（j=0;j<n;j++）
{
if（jd[j].s==3）//判断为平衡节点{continue;}elseif（i==j）//PQ节点i=j时{floatsum1=0,sum2=0;for（kk=0;kk<n;kk++）{sum1+=Y[i][kk][0]*jd[kk].f+Y[i][kk][1]*jd[kk].e;sum2+=Y[i][kk][0]*jd[kk].e-Y[i][kk][1]*jd[kk].f;}yakebi[ii][jj]=-sum1-Y[i][i][0]*jd[i].f+Y[i][i][1]*jd[i].e;//计算Hyakebi[ii][jj+1]=-sum2-Y[i][i][0]*jd[i].e-Y[i][i][1]*jd[i].f;//计算Nyakebi[ii+1][jj]=-sum2+Y[i][i][0]*jd[i].e+Y[i][i][1]*jd[i].f;//计算Jyakebi[ii+1][jj+1]=sum1-Y[i][i][0]*jd[i].f+Y[i][i][1]*jd[i].e;//计算Ljj+=2;
第6/13页
}}else//PQ节点i!
=j时{yakebi[ii][jj]=Y[i][j][1]*jd[i].e-Y[i][j][0]*jd[i].f;//计算Hyakebi[ii][jj+1]=-Y[i][j][0]*jd[i].e-Y[i][j][1]*jd[i].f;//计算Nyakebi[ii+1][jj]=-yakebi[ii][jj+1];//计算Jyakebi[ii+1][jj+1]=yakebi[ii][jj];//计算Ljj+=2;}}ii+=2;elseif（jd[i].s==2）//判断为PV节点{for（j=0;j<n;j++）{if（jd[j].s==3）//判断为平衡节点{continue;}elseif（i==j）//PV节点i=j时{floatsum1=0,sum2=0;for（kk=0;kk<n;kk++）{sum1+=Y[i][kk][0]*jd[kk].f+Y[i][kk][1]*jd[kk].e;sum2+=Y[i][kk][0]*jd[kk].e-Y[i][kk][1]*jd[kk].f;}yakebi[ii][jj]=-sum1-Y[i][i][0]*jd[i].f+Y[i][i][1]*jd[i].e;//计算Hyakebi[ii][jj+1]=-sum2-Y[i][i][0]*jd[i].e-Y[i][i][1]*jd[i].f;//计算Nyakebi[ii+1][jj]=-2*jd[i].f;//计算Ryakebi[ii+1][jj+1]=-2*jd[i].e;//计算Sjj+=2;}else//PV节点i!
=j时{yakebi[ii][jj]=Y[i][j][1]*jd[i].e-Y[i][j][0]*jd[i].f;//计算Hyakebi[ii][jj+1]=-Y[i][j][0]*jd[i].e-Y[i][j][1]*jd[i].f;//计算Nyakebi[ii+1][jj]=0;//计算Ryakebi[ii+1][jj+1]=0;//计算Sjj+=2;}}ii+=2;
第7/13页
}
}
//输出雅克比矩阵
fprintf（fp2,"\n雅克比矩阵为:
\n"）;
for（i=0;i<2*（n-1）;i++）
{
fprintf（fp2,"\n"）;
for（j=0;j<2*（n-1）;j++）
{
fprintf（fp2,"%10f",yakebi[i][j]）;
}
}
fprintf（fp2,"\n==================================================================================\n"）;
}
//用列主元法解修正方程
voidSlove_Eq
uaions（）
{
j=0;
for（i=0;i<n-1;i++）//把不平衡量加入雅克比矩阵形成增广矩阵
{
if（jd[i].s==1）
{
yakebi[j][2*（n-1）]=jd[i].dp;//把pq节点不平衡量加入
yakebi[j+1][2*（n-1）]=jd[i].dq;//把pq节点不平衡量加入
j+=2;
}
elseif（jd[i].s==2）
{
yakebi[j][2*（n-1）]=jd[i].dp;//把pv节点不平衡量加入
yakebi[j+1][2*（n-1）]=jd[i].dv;//把pv节点不平衡量加入
j+=2;
}
}
doublerr=0;//临时变量
for（i=0;i<2*（n-1）;i++）//修正方程消元
{
rr=fabs（yakebi[i][i]）;//选主元
kk=i;
for（j=i+1;j<2*（n-1）;j++）//选主元
{
if（fabs（yakebi[j][i]）>rr）
第8/13页
{
kk=j;
rr=fabs（yakebi[j][i]）;
}
}
if（rr==0）
{
printf（"修正方程奇异无解!
!
"）;
exit（0）;
}
elseif（kk!
=i）//判断是否需要换行
{
}for（j=i;j<=2*（n-1）;j++）//选主元交换行{floatxx=0;xx=yakebi[i][j];yakebi[i][j]=yakebi[kk][j];yakebi[kk][j]=xx;}}floatnum3=yakebi[i][i];//临时变量for（j=i;j<=2*（n-1）;j++）{yakebi[i][j]=yakebi[i][j]/num3;//规格化}for（j=i+1;j<2*（n-1）;j++）{floatnum4=yakebi[j][i];//临时变量for（kk=i;kk<=2*（n-1）;kk++）{yakebi[j][kk]=yakebi[j][kk]-yakebi[i][kk]*num4;//消元}}for（i=2*（n-1）-1;i>=0;i--）//回代过程{rr=0;for（j=i+1;j<2*（n-1）;j++）{rr+=D[j]*yakebi[i][j];}D[i]=yakebi[i][2*（n-1）]-rr;}
第9/13页
j=0;
for（i=0;i<n-1;i++）//求得各节点的修正电压量
{
jd[i].df=D[j];//赋值回各节点的修正电压量
jd[i].de=D[j+1];//赋值回各节点的修正电压量
jd[i].f-=jd[i].df;//求得各节点的修正后电压量
jd[i].e-=jd[i].de;//求得各节点的修正后电压量
j+=2;
}
h=0;
for（i=0;i<2*（n-1）;i++）
{
if（fabs（D[i]）>eps）
{
h=1;
break;
}
}
//输出修正量
fprintf（fp2,"\n电压修正量为:
\n"）;
for（i=0;i<n-1;i++）
{
fprintf（fp2,"节点%d:
de=%-10fdf=%-10f\n",jd[i].num,jd[i].de,jd[i].df）;}
fprintf（fp2,"\n==================================================================================\n"）;
//输出节点电压值
fprintf（fp2,"\n节点电压值为:
\n"）;
for（i=0;i<n-1;i++）
{
fprintf（fp2,"节点%d:
e=%-10ff=%-10f\n",jd[i].num,jd[i].e,jd[i].f）;}
fprintf（fp2,"\n==================================================================================\n"）;
}
//潮流结果计算及输出
voidPowerflow_Results（）
{
fprintf（fp2,"潮流计算结果"）;
fprintf（fp2,"\n\n***********************************************************************************\n\n"）;
fprintf（fp2,"\
升级会员 