pq.txt

pq分解法潮流计算 用于高压网配电的潮流计算

#include "stdio.h" 
#include "math.h" 

#define M 20 /*节点数、支路数极限值*/ 
#define N 5 /*迭代次数极限值*/ 

int n,m,dd=0,pq=0,pv=0,tt,qy; 
float eps; /*节点数、支路数、对地支路数、pq节点数、pv节点数、控制打印、互联网数、精度*/ 

struct jiedian /*节点数据结构体*/ 
{int s; /*节点类型：1-PQ节点；2-PV节点；3-平衡节点*/ 
float p,q,e,f,v; /*节点的有功、无功、电压实部、虚部、电压辐值*/ 
}jiedian[M]; /*如引用节点1的无功，则为jiedian[1].q */ 

struct zhilu /*支路数据结构体*/ 
{int p1,p2,s; /*支路两端节点号，支路类型：1-普通支路；2-变压器支路；3-对地支路*/ 
float r,x,b,kt; /*支路的电阻、电抗、导纳、变压器非标准变比*/ 
}zhilu[M]; 

struct hulianwang /*互联网结构体*/ 
{int num,pv; /*互联区域号,指定的pv节点号*/ 
float p,eps2; /*规定的有功功率及其允许误差*/ 
int count; /*每个互联区域包括的节点数*/ 
int a[M]; /*包括的节点*/ 
} hulianwang[M]; 


static float G[M][M],B[M][M],G1[M][M],B1[M][M]; /*节点导纳阵矩阵*/ 
static float ykb[2*M][2*M]; /*节点导纳阵矩阵*/ 
static float yinzi[2*M][2*M]; /*因子表*/ 
static float P[M][M],Q[M][M]; /*潮流计算结果*/ 

FILE *fp1,*fp2; /*文件输入、输出指针*/ 


void input() /*从文件in.txt 输入线路基本参数、节点数据、支路数据*/ 
{ 
int i,j,h; /*循环变量、节点（支路）类型*/ 

/**********************************打开电网参数文件**********************/ 
fp1=fopen("in.txt","r"); 
if (fp1==NULL) 
{ 
printf("Can not open file in.txt !\n"); 
exit(0); 
} 

/***********读取节点数,支路数,互联网数,PQ节点数PV节点数和精度***********/ 
fscanf(fp1,"%d,%d,%d,%f\n",&n,&m,&qy,&eps); 

/************************读取节点信息*************************************/ 
for(i=1;i<=n;i++) 
{ 
fscanf(fp1,"%d",&h); 
if(h==1) /*类型h=1是PQ节点*/ 
{ 
pq++; 
fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].q,&jiedian[i].e,&jiedian[i].f); 
jiedian[i].s=1; 
jiedian[i].v=sqrt(jiedian[i].e*jiedian[i].e+jiedian[i].f*jiedian[i].f); 
} 
else if(h==2) /*类型h=2是pv节点*/ 
{ 
pv++; 
fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].v,&jiedian[i].e,&jiedian[i].f); 
jiedian[i].s=2; 
jiedian[i].q=0; 
} 
else /*类型h=3是平衡节点*/ 
{ 
jiedian[i].p=0; 
jiedian[i].q=0; 
jiedian[i].e=1; 
jiedian[i].f=0; 
jiedian[i].v=1; 
jiedian[i].s=3; 
} 
} 

/*******************************读取支路信息*************************************/ 
for(i=1;i<=m;i++) 
{ 
fscanf(fp1,"%d",&h); 
if(h==1) /*类型h=1是普通支路*/ 
{ 
fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].b); 
zhilu[i].kt=1; 
zhilu[i].s=1; 
} 
if(h==2) /*类型h=2是变压器支路*/ 
{ 
fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].kt); 
zhilu[i].s=2; 
} 
if(h==3) /*类型h=3是接地支路*/ 
{ 
fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].b); 
zhilu[i].kt=1; 
zhilu[i].s=3; 
dd++; 
} 
} 

/*******************************读互联网信息*************************************/ 
if(qy!=0) 
for(i=1;i<=qy;i++) /*输入互联网状况*/ 
{ 
fscanf(fp1,"%d,%d,%f,%f,%d",&hulianwang[i].num,&hulianwang[i].pv,&hulianwang[i].p,&hulianwang[i].eps2,&hulianwang[i].count); 
for(j=1;j<=hulianwang[i].count;j++) 
fscanf(fp1,",%d",&(hulianwang[i].a[j])); 
} 

fclose(fp1); 

/**********************************打开输出结果文件********************************/ 
fp2=fopen("out.txt","w"); 
if(fp2==NULL) 
{ 
printf("Can not open file!\n"); 
exit(0); 
} 

fprintf(fp2,"\n**************************** 原始数据如下 ***********************************\n\n "); 
fprintf(fp2," 节点数:%2d 支路数:%2d 对地支路数:%2d PQ节点数:%2d PV节点数:%2d 精度1:%f \n",n,m,dd,pq,pv,eps); 

fprintf(fp2,"\n------------------------------------------------------------------------------\n\n "); 
for(i=1;i<=pq;i++) 
fprintf(fp2," 节点%2d PQ节点 P[%d]=%f Q[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].q); 

for(i=pq+1;i<=pq+pv;i++) 
fprintf(fp2," 节点%2d PV节点 P[%d]=%f V[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].v); 
fprintf(fp2," 节点%2d 平衡节点\n",i); 
fprintf(fp2,"\n-------------------------------------------------------------------------------\n\n "); 
for(i=1;i<=m;i++) 
{ 
if(zhilu[i].s==1) 
fprintf(fp2," 支路%2d：普通支路 相关节点:%2d,%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].b); 
else if(zhilu[i].s==2) 
fprintf(fp2," 支路%2d：变压器支路 相关节点:%2d,%2d R=%f X=%f Kt=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].kt); 
else 
fprintf(fp2," 支路%2d： 对地支路 相关节点:%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].r,zhilu[i].x,zhilu[i].b); 

} 
for(i=1;i<=qy;i++) 
{ 
fprintf(fp2," 互联网区域:%2d 指定pv节点:%2d 规定有功功率:%8.5f 允许误差:%8.5f 包括的节点: ",hulianwang[i].num,hulianwang[i].pv,hulianwang[i].p,hulianwang[i].eps2); 
for (j=1;j<=hulianwang[i].count;j++) 
fprintf(fp2,"%2d ",hulianwang[i].a[j]);fprintf(fp2,"\n"); 
} 

} 

void youhua() /*利用节点数据、支路数据形成新编号*/ 
{ 
int i,k,j,jd[N]; /*jd[N]记录节点*/ 
struct jiedian tem; /*中间临时节点所连支路数*/ 

for(i=1;i<=n;i++) /*对节点连接支路数赋0*/ 
{ 
jd[i]=0; 
} 

for(i=1;i<=n;i++) /*考虑每个节点*/ 
{ 
for(j=1;j<=m;j++) /*考虑每条支路*/ 
{ 
if((zhilu[j].p1==i||zhilu[j].p2==i)&&(zhilu[j].s!=3)) /*如果是这条支路的节点且不是接地支路*/ 
jd[i]++; 

} 
} 

for(i=1;i<pq;i++) /*对pq节点按所连支路的个数进行重新编号*/ 
{ 
for(j=i+1;j<=pq;j++) /*从小到大排序算法*/ 
{ 
if(jd[i]>jd[j]) 
{ 
tem=jiedian[i]; 
jiedian[i]=jiedian[j]; 
jiedian[j]=tem; 


for(k=1;k<=m;k++) /*更新支路信息*/ 
{ 
if(zhilu[k].p1==j) 
zhilu[k].p1=i; 
else if(zhilu[k].p2==j) 
zhilu[k].p2=i; 
else if(zhilu[k].p1==i) 
zhilu[k].p1=j; 
else if(zhilu[k].p2==i) 
zhilu[k].p2=j; 
} 

} 
} 
} 

for(i=pq+1;i<pq+pv;i++) /*对pv节点按所连支路的个数进行重新编号*/ 
{ 
for(j=i+1;j<=pq+pv;j++) /*从小到大排序算法*/ 
{ 
if(jd[i]>jd[j]) 
{ 
tem=jiedian[i]; 
jiedian[i]=jiedian[j]; 
jiedian[j]=tem; 


for(k=1;k<=m;k++) /*更新支路信息*/ 
{ 
if(zhilu[k].p1==j) 
zhilu[k].p1=i; 
else if(zhilu[k].p2==j) 
zhilu[k].p2=i; 
else if(zhilu[k].p1==i) 
zhilu[k].p1=j; 
else if(zhilu[k].p2==i) 
zhilu[k].p2=j; 
} 

} 
} 
} 

fprintf(fp2,"\n**************************** 节点优化结果如下 *********************************\n\n "); 

for(i=1;i<=n;i++) 
{ 
if(jiedian[i].s==1) 
fprintf(fp2," 节点%2d PQ节点 P[%d]=%f Q[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].q); 
if(jiedian[i].s==2) 
fprintf(fp2," 节点%2d PV节点 P[%d]=%f V[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].v); 
if(jiedian[i].s==3) 
fprintf(fp2," 节点%2d 平衡节点 \n",i); 
} 
} 


void form_y() /*利用支路数据形成Y，注意对地支路*/ 
{ 
int i,j,k; 
float S; 

for(i=1;i<=n;i++) 
for(j=1;j<=n;j++) 
G[i][j]=B[i][j]=0; 

for(i=1;i<=m;i++) /*节点导纳矩阵的主对角线上的导纳*/ 
for(j=1;j<=n;j++) 
if(zhilu[i].p1==j||zhilu[i].p2==j) 
{ 
S=zhilu[i].r*zhilu[i].r+zhilu[i].x*zhilu[i].x; 
if(S==0) continue; 
G[j][j]+=zhilu[i].r/S; 
B[j][j]+=-zhilu[i].x/S; 

if(zhilu[i].s==1) /*如果是普通支路*/ 
{ 
B[j][j]+=zhilu[i].b/2; 
} 
else if(zhilu[i].s==2) /*如果是普通变压器支路*/ 
{ 
if(zhilu[i].p1==j) 
{ 
G[j][j]+=(zhilu[i].r/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt); 
B[j][j]+=(-zhilu[i].x/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt); 
} 
else if(zhilu[i].p2==j) 
{ 
G[j][j]+=(zhilu[i].r/S*(zhilu[i].kt-1))/zhilu[i].kt; 
B[j][j]+=(-zhilu[i].x/S*(zhilu[i].kt-1))/zhilu[i].kt; 
} 
} 
else if(zhilu[i].s==3) /*如果是对地支路*/ 
{ 
B[j][j]+=zhilu[i].b; 
} 
} 


for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上的导纳*/ 
{ 
i=zhilu[k].p1; 
j=zhilu[k].p2; 
S=zhilu[k].r*zhilu[k].r+zhilu[k].x*zhilu[k].x; 
if(S==0) continue; 
G[i][j]+=-zhilu[k].r/S; 
B[i][j]+=zhilu[k].x/S; 
if(zhilu[k].kt!=1.0) 
{ 
G[i][j]/=zhilu[k].kt; 
B[i][j]/=zhilu[k].kt; 
} 
G[j][i]=G[i][j]; 
B[j][i]=B[i][j]; 
} 


fprintf(fp2,"\n\n****************************** 节点导纳矩阵为 **********************************\n"); 
for(i=1;i<=n;i++) 
{ 
fprintf(fp2,"\n "); 
for(j=1;j<=n;j++) 
fprintf(fp2,"%8.5f+j%8.5f ",G[i][j],B[i][j]); 
} 
} 

void form_j() /*利用节点数据和Y形成J*/ 

{ 
float ei,fi,a=0,b=0; 
int i1,j1,k1,i,j,k; 

for(i=1;i<=2*(pq+pv)+1;i++) 
for(j=1;j<=2*(pq+pv)+1;j++) 
ykb[i][j]=0; 

for(i=1;i<=pq;i++) 
for(j=1;j<n;j++) 
{ i1=i; 
j1=j; 
ei=jiedian[i].e; 
fi=jiedian[i].f; 

if(i!=j) /*求i!=j时的H、N、J、L*/ 
{ ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /* H */ 
ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /* N */ 
ykb[2*i][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /* J */ 
ykb[2*i][2*j]=-G[i1][j1]*ei-B[i1][j1]*fi; /* L */ 
} 


else /*求i=j时的H、N、J、K*/ 
{ a=0;b=0; 
for(k=1;k<=n;k++) 
if(k!=i) 
{ k1=k; 
a=a+G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f; 
b=b+G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e; 
} 
ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ 
ykb[2*i][2*j]=-2*B[i1][i1]*jiedian[i].f+a; /*L*/ 
ykb[2*i-1][2*j]= 2*G[i1][i1]*jiedian[i].f+b; /*N*/ 
ykb[2*i][2*j-1]=-2*B[i1][i1]*jiedian[i].e-b; /*J*/ 
} 
} 


for(i=pq+1;i<=pq+pv;i++) /*形成pv节点子阵*/ 
for(j=1;j<n;j++) 
{ i1=i; 
j1=j; 
ei=jiedian[i].e; 
fi=jiedian[i].f; 

if(i!=j) /*求i!=j时的H、N*/ 
{ 
ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /*H*/ 
ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /*N*/ 
} 


else /*求i=j时的H、N、R、S*/ 
{ a=0;b=0; 
for(k=1;k<=n;k++) 
if(k!=i) 
{ k1=k; 
a+=G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f; 
b+=G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e; 
} 
ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ 
ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].f+b; /*N*/ 
ykb[2*i][2*j-1]=2*ei; /*R*/ 
ykb[2*i][2*j]=2*fi; /*S*/ 
} 
} 


fprintf(fp2,"\n\n\n******************************* 雅可比矩阵*************************************** \n"); 

for(i=1;i<=2*(pq+pv);i++) 
{ 
fprintf(fp2,"\n "); 
for(j=1;j<=2*(pq+pv);j++) 
fprintf(fp2,"%8.5f ",ykb[i][j]); 
} 
} 


void form_yz() /*利用J形成因子表*/ 
{ 
int i,j,k; 
float a[2*M][2*M],L[2*M][2*M],R[2*M][2*M]; 
float x,y,z; 

for(i=1;i<=2*(pq+pv);i++) 
for(j=1;j<=2*(pq+pv);j++) 
{ a[i][j]=0;R[i][j]=0; L[i][j]=0;} 
for(i=1;i<=2*(pq+pv);i++) 
for(j=1;j<=2*(pq+pv);j++) 
a[i][j]=ykb[i][j]; 

for(i=1;i<=2*(pq+pv);i++) 
{ 
for(j=1;j<=2*(pq+pv);j++) 
{ 
if(i>j) 
{ 
R[i][j]=0; 
x=0; 
for(k=1;k<=j-1;k++) 
{ 
x=x+L[i][k]*R[k][j]; 
} 
L[i][j]=a[i][j]-x; 
} 
if(i==j) 
{ 
R[i][j]=1; 
y=0; 
for(k=1;k<=i-1;k++) 
{ 
y=y+L[i][k]*R[k][i]; 
} 
L[i][i]=a[i][j]-y; 
} 
if(i<j) 
{ 
L[i][j]=0; 
z=0; 
for(k=1;k<=i-1;k++) 
{ 
z=z+L[i][k]*R[k][j]; 
} 
R[i][j]=(a[i][j]-z)/L[i][i]; 
} 
} 
} 

for(i=1;i<=2*(pq+pv);i++) 
{ 
for(j=1;j<=2*(pq+pv);j++) 
{ 
if(i<j) 
yinzi[i][j]=R[i][j]; 
if(i==j) 
yinzi[i][j]=1/L[i][j]; 
if(i>j) 
yinzi[i][j]=L[i][j]; 
} 
} 

fprintf(fp2,"\n\n\n******************************** 因子表 ***************************************** \n"); 

for(i=1;i<=2*(pq+pv);i++) 
{ 
fprintf(fp2,"\n "); 
for(j=1;j<=2*(pq+pv);j++) 
fprintf(fp2,"%8.5f ",yinzi[i][j]); 
} 

}