潮流计算的pq分解法用c语言实现的.docx
《潮流计算的pq分解法用c语言实现的.docx》由会员分享，可在线阅读，更多相关《潮流计算的pq分解法用c语言实现的.docx（47页珍藏版）》请在冰豆网上搜索。
潮流计算的pq分解法用c语言实现的
//////////////////////////////////////////////////////////////////////
//PQ分解法潮流//
//文件输入格式：
节点总数n（包括联络节点）,支路数zls//
//节点数（发电机和负荷）nb,接地电抗数mdk,迭代精度eps//
//考虑负荷静特性标志kk2（0考虑）,平衡节点号，优化标志（0不优化）//
//最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[]//
//支路容纳zyk[],节点号nob[]及标志nobt[]（0-PQ,-1-PV）//
//发电机和负荷有功、无功pg[],qg[],pl[],ql[]//
//电压v0[]（pv节点输入实际值，PQ节点任输入一值）//
//电抗节点号idk[],电抗值dkk[]//
//////////////////////////////////////////////////////////////////////
#include"math.h"
#include"stdio.h"
#defineNS2000//最大节点数
#defineNS2NS*2
#defineNS41000//NS4、NS必须大于2*zls。
#defineZS3000//最大支路数
#defineZS2ZS*2
#defineDKS200//最大电抗器数
#defineN2ZS*4
#defineN3ZS*8+NS*4
FILE*fp1,*fp2;
charinname[12],outname[12];
//fp1输入数据文件指针fp2输出文件指针
//inname[]输入数据文件名outname[]输出数据文件名
intn,zls,nb,mdk,mpj,bnsopton,it1,dsd,kk2,nzls;
//节点总数n（包括联络节点）支路数（回路数）zls节点数nb（发电机和负荷）
//接地电抗数mdk精度eps平衡节点号mpj
//节点优化（标志）bnsopton（=0节点不优化，!
=0节点优化）
//最大迭代次数it1最低电压或最大功率误差节点号dsd
//负荷静特性标志（=0考虑负荷静特性）
//支路数（双回线算一条支路）
intizl[ZS],jzl[ZS],idk[DKS],yds[NS],ydz[NS],iy[ZS2];
//izl[],jzl[],idk[]:
分别存放左、右节点号和电抗器节点号。
//yds[]存放各行非零非对角元素的个数。
//ydz[i]是第i行第一个非零非对角元素的首地址，
//即在所有非零非对角元素中的次序号
//iy[]存放列足码。
intnnew[NS4],old[NS],nob[NS],nobt[NS];
//nnew[],old[]存放的是新、旧节点号。
//nnew[i]中为i对应的新号
//nob[]存放的是节点号。
nobt[]存放的是节点类型,0:
pq节点,-1:
pv节点。
doubleeps,dsm,vmin,dph,dqh,af[3];
//eps迭代收敛精度，dsm最大功率误差
//vmin:
系统最低电压值。
dph,dqh:
系统有、无功损耗。
//af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。
doublev00;
//v00:
系统平均电压ci,cj分别作为节点i,j的电压相角的临时存储单元。
doublezr[ZS],zx[ZS],zyk[ZS],dkk[DKS],gii[NS],bii[NS],yg[ZS2],yb[ZS2];
doublepg[NS],qg[NS],pl[NS],ql[NS],v0[NS],v[NS],va[NS];
//支路电阻zr[]支路电抗zx[]输电线路充电容纳zyk[]（y0/2）
//接地电抗dkk[]对角元实部gii[]对角元虚部
//非对角元实部yg[]非对角元虚部yb[]
//pg[],qg[],pl[],ql[]:
发电机，负荷功率实、虚部
//v[]是电压幅值,va[]是电压相角。
doublew[NS2],kg[3],b[NS2];
intnewsort[NS4];
//newsort[i]存放i对应的老号
voidinitial（）;
voidpqflow（）;
voidout（）;
voiddataio（）;
voidbnsopt（）;
voidzlsort（int*nnew）;
voidprinto（）;
voidprinty（）;
voidy2（）;
voidya0（）;
voidyzb（）;
voidjdgl（intkq0）;
voidbbhl（intkq0）;
voidcalc（）;
intiabs（inta）;
voidbranch_output（）;
voidnewval（double*aa）;
voidprintc（）;
voidiswap（）;
voidswap（）;
voidprintf2（double*aa,double*bb,intn）;
voidcalc（int*iu,double*u,double*di,int*nfd,double*b）;
voidprinti（int*aa,intn）;
voidprintf1（double*aa,intn）;
intfind（intk,inta[],int*z）;
voidyzb（intt,int*iu,double*u,double*di,int*nfd）;
intisgn（inta,intb）;
voidyy1（）;
voidy3（）;
voidnewtoold（）;
intmain（void）
{
initial（）;//初始化
pqflow（）;//pq潮流计算
out（）;//输出节点和支路数据
return1;
}
intisgn（inta,intb）
{
//****本函数功能返回值为a的绝对值b的符号****//
//参数1提供值，参数2提供符号//
if（b<0）
if（a>0）
a=-a;
returna;
}
intfind（intk,inta[],int*z）
{
//****本函数查找a[]中是否有fabs（k）有则返回0，无则返回1****//
//参数1为待查找量，参数2待搜索数组，参数3返回k在a[]中的次序号//
inti;
for（i=1;i<=n;i++）
if（iabs（k）==a[i]）
{
*z=i;
return1;
}
return0;
}
voidoldtonew（）
{
//****本函数将输入数据中的节点号变成从1开始的连续节点号****//
inti,j,k,ii1,ii2,zls2,k1,k2,k3,k4,ip;
zls2=zls+zls;
for（i=1;i<=zls2;i++）
newsort[i]=0;
ii1=0;
for（i=1;i<=zls;i++）
{
k=izl[i];
if（!
find（k,newsort,&ii2））
{
ii1++;
newsort[ii1]=iabs（k）;
}
k=jzl[i];
if（!
find（k,newsort,&ii2））
{
ii1++;
newsort[ii1]=iabs（k）;
}
}
for（i=1;i<=ii1-1;i++）
{
for（j=i+1;j<=ii1;j++）
{
if（newsort[i]>newsort[j]）
{
k=newsort[i];
newsort[i]=newsort[j];
newsort[j]=k;
}
}
}
for（i=1;i<=zls;i++）
{
k=izl[i];
if（find（k,newsort,&ii2））
{
izl[i]=isgn（ii2,k）;
}
else
printf（"error!
"）;
k=jzl[i];
if（find（k,newsort,&ii2））
{
jzl[i]=isgn（ii2,k）;
}
else
printf（"error!
"）;
printf（"izl[%d]=%d,jzl[%d]=%d\n",i,izl[i],i,jzl[i]）;
}
for（i=1;i<=nb;i++）
{
for（j=1;j<=n;j++）
if（nob[i]==newsort[j]）
{
nob[i]=j;
break;
}
printf（"nob[%d]=%d\n",i,nob[i]）;
}
for（j=1;j<=n;j++）
{
if（mpj==newsort[j]）
{
mpj=j;
break;
}
}
//电抗器节点号转变
for（j=1;j<=mdk;j++）
{
for（i=1;i<=n;i++）
{
if（idk[j]==newsort[i]）
{
idk[j]=i;
break;
}
}
}
}
voidinitial（）
{
//****本函数进行初始化工作****//
inti,k1;
dataio（）;//输入原始数据
oldtonew（）;//转化为新号
if（bnsopton==0）//节点不优化，新节点号即为老节点号。
for（i=1;i<=n;i++）
{
old[i]=i;
nnew[i]=i;
}
else
bnsopt（）;//节点优化
mpj=nnew[mpj];//mpj:
平衡节点
zlsort（nnew）;//sortther,xandb
for（i=1;i<=mdk;i++）
{
k1=idk[i];
idk[i]=nnew[k1];
}
for（i=1;i<=n;i++）
{
v[i]=v00;
va[i]=0.0;
}//所有节点的电压幅值初值都为1.000（v00），电压相角初值都为0。
//exchangethenodebeforeandaftersort
for（i=1;i<=n;i++）
yds[i]=0;//theimmediate
for（i=1;i<=nb;i++）
{
k1=nnew[nob[i]];
yds[k1]=nobt[i];
}
for（i=1;i<=n;i++）
nobt[i]=yds[i];
newval（pg）;
newval（qg）;
newval（pl）;
newval（ql）;
newval（v0）;
for（i=1;i<=n;i++）//nobt[]istypeofnode
if（nobt[i]==-1）
v[i]=v0[i];//nob[]isserialsnumbe
//nobt[]=-1:
pv节点，v0[]存放的是最后一个节点数据，
//对于pv节点，即为该点应维持的电压值。
//nobt[]=0:
pq节点，v0[]存放的是最后一个节点数据，
//对于pq节点，即为系统平均电压值。
printo（）;
//输出af[]、v00和节点排序后的支路、节点和
//接地电抗数据（仅仅查看中间结果）
ya0（）;//获得yds[]、ydz[]、列足码iy[]。
（P407）
}
voidprinto（）
{
//****输出af[]、v00和节点排序后的支路、节点和接地电抗数据****//
inti;
fprintf（fp2,"\n******AFANDV0******\n"）;
fprintf（fp2,"\n%7.3f%7.3f%7.3f\n",af[0],af[1],v00）;
printc（'-',78）;
fprintf（fp2,"\n\n*******ZLB*******\n"）;
for（i=1;i<=zls;i++）
{
fprintf（fp2,"\n"）;
fprintf（fp2,"%8d%8d%8d%8d",izl[i],jzl[i],old[abs（izl[i]）],old[abs（jzl[i]）]）;
fprintf（fp2,"%9.4f%9.4f%9.4f",zr[i],zx[i],zyk[i]）;
}
printc（'-',78）;
fprintf（fp2,"\n\n*******BUS*******\n"）;
for（i=1;i<=nb;i++）
{
fprintf（fp2,"\n"）;
fprintf（fp2,"%8d%8d%8d",nob[i],old[nob[i]],nobt[i]）;
fprintf（fp2,"%9.4f%9.4f%9.4f%9.4f%9.4f",pg[i],qg[i],pl[i],ql[i],v0[i]）;
}
printc（'-',78）;
fprintf（fp2,"\n\n******DKK******\n"）;
for（i=1;i<=mdk;i++）
{
fprintf（fp2,"\n"）;
fprintf（fp2,"%8d%8d%7.4f",idk[i],old[idk[i]],dkk[i]）;
}
}
voiddataio（）
{
//****系统数据初始化****//
inti;
af[0]=0.6;
af[1]=2.0;//af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。
v00=1.000;//系统平均电压
printf（"\npleaseinputthenameofdatafile\n"）;
scanf（"%s",inname）;
fp1=fopen（inname,"r"）;
printf（"\npleaseoutputthenameofdatafile\n"）;
scanf（"%s",outname）;
fp2=fopen（outname,"w"）;
fscanf（fp1,"%d%d%d%d",&n,&zls,&nb,&mdk）;
//thenumberofnode,branches,node
fscanf（fp1,"%lf%d%d%d%d",&eps,&kk2,&mpj,
&bnsopton,&it1）;
//precision,swingnode,sortthenode,iterationnumbers
for（i=1;i<=zls;i++）
{
fscanf（fp1,"%d%d",&izl[i],&jzl[i]）;
fscanf（fp1,"%lf%lf%lf",&zr[i],&zx[i],&zyk[i]）;
}
for（i=1;i<=nb;i++）
{
fscanf（fp1,"%d%d",&nob[i],&nobt[i]）;
fscanf（fp1,"%lf%lf%lf%lf%lf",&pg[i],&qg[i],&pl[i],
&ql[i],&v0[i]）;
}
for（i=1;i<=mdk;i++）
{
fscanf（fp1,"%d%lf",&idk[i],&dkk[i]）;
}
fclose（fp1）;
}
voidpqflow（）
{
//****PQ分解法计算潮流,程序框图见P164图3-16（从第7步起）****//
intkq0,iu1[N2],nfd1[NS],iu2[N2],nfd2[NS];
inti,t;
doubleu1[N2],u2[N2],di1[NS],di2[NS];
yy1（）;
yzb（0,iu1,u1,di1,nfd1）;//formtheBmatrixofP-0iteration
y2（）;
yzb（1,iu2,u2,di2,nfd2）;//formtheBmatrixofQ-Viteration
t=0;
kq0=0;
kg[0]=kg[1]=1;
do
{
jdgl（kq0）;//calculatingthepower
bbhl（kq0）;//findoutthemaxi
if（kq0==0）
printf（"P:
%d\t%d\t%f\n",t,dsd,dsm）;
else
printf（"Q:
%d\t%d\t%f\n",t,dsd,dsm）;
if（fabs（dsm）>eps）
{
kg[kq0]=1;
if（kq0==0）
calc（iu1,u1,di1,nfd1,b）;
if（kq0==1）
calc（iu2,u2,di2,nfd2,b）;
for（i=1;i<=n;i++）
{
if（kq0==0）
va[i]=va[i]-b[i]/v00;
else
v[i]=v[i]-b[i];
}
}
else
kg[kq0]=0;
if（kq0==0）
kq0=1;
else
{
kq0=0;
t++;
}
if（t>it1）
break;
}while（（fabs（dsm）>eps）||（kg[kq0]!
=0））;
fprintf（fp2,"\n%s%d","times=",t）;
}
voidout（）
{
//****本函数输出节点和支路数据****//
zlsort（old）;//recoverthedataifsorted
//newtoold（）;
node_output（）;//nodedata
branch_output（）;//branchdata
printc（'-',78）;
printc（'*',78）;
fprintf（fp2,"\n"）;
}
voidnewval（double*aa）
{
//****本函数将旧号换成新号****//
inti,k1;
for（i=1;i<=n;i++）
b[i]=0.0;
for（i=1;i<=nb;i++）
{
k1=nnew[nob[i]];
b[k1]=aa[i];
}
for（i=1;i<=n;i++）
aa[i]=b[i];
}
voidyzb（intt,int*iu,double*u,double*di,int*nfd）
{
//****本函数求因子表****//
//参数1为标志（t=0求B',t=1求B''）//
//参数2因子表上三角矩阵非零非对角元素的列足码
//参数3因子表上三角矩阵非零非对角元素的数值
//参数4因子表上三角矩阵对角元素
//参数5因子表上三角各行非零元素个数
inti,j,k,i1,i2;
intjj,jj1,jj2,im,x,fd[NS];
doubleai,b[NS];
nfd[1]=1;
for（i=1;i<=n;i++）
{
//nobt[]存放的是节点类型,0:
pq节点,-1:
pv节点。
if（（（t!
=1）||（nobt[i]!
=-1））&&i!
=mpj）//<---|
{//|
for（j=i+1;j<=n;j++）//|
b[j]=0.0;//|
b[i]=bii[i];//|
if（（kk2==0）&&（t==1）&&（nobt[i]!
=-1））//存在（t==1）的情况，不多余。
b[i]=b[i]+af[1]*ql[i]/v0[i]/v0[i];//af[1]
i1=ydz[i];
i2=ydz[i+1]-1;
for（j=i1;j<=i2;j++）
{
k=iy[j];
b[k]=yb[j];
}
b[mpj]=0.0;
if（t==1）
for（j=1;
升级会员 