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

#define M 20 /*節(jié)點數(shù)、支路數(shù)極限值*/ 
#define N 5 /*迭代次數(shù)極限值*/ 

int n,m,dd=0,pq=0,pv=0,tt,qy; 
float eps; /*節(jié)點數(shù)、支路數(shù)、對地支路數(shù)、pq節(jié)點數(shù)、pv節(jié)點數(shù)、控制打印、互聯(lián)網(wǎng)數(shù)、精度*/ 

struct jiedian /*節(jié)點數(shù)據(jù)結構體*/ 
{int s; /*節(jié)點類型：1-PQ節(jié)點；2-PV節(jié)點；3-平衡節(jié)點*/ 
float p,q,e,f,v; /*節(jié)點的有功、無功、電壓實部、虛部、電壓輻值*/ 
}jiedian[M]; /*如引用節(jié)點1的無功，則為jiedian[1].q */ 

struct zhilu /*支路數(shù)據(jù)結構體*/ 
{int p1,p2,s; /*支路兩端節(jié)點號，支路類型：1-普通支路；2-變壓器支路；3-對地支路*/ 
float r,x,b,kt; /*支路的電阻、電抗、導納、變壓器非標準變比*/ 
}zhilu[M]; 

struct hulianwang /*互聯(lián)網(wǎng)結構體*/ 
{int num,pv; /*互聯(lián)區(qū)域號,指定的pv節(jié)點號*/ 
float p,eps2; /*規(guī)定的有功功率及其允許誤差*/ 
int count; /*每個互聯(lián)區(qū)域包括的節(jié)點數(shù)*/ 
int a[M]; /*包括的節(jié)點*/ 
} hulianwang[M]; 


static float G[M][M],B[M][M],G1[M][M],B1[M][M]; /*節(jié)點導納陣矩陣*/ 
static float ykb[2*M][2*M]; /*節(jié)點導納陣矩陣*/ 
static float yinzi[2*M][2*M]; /*因子表*/ 
static float P[M][M],Q[M][M]; /*潮流計算結果*/ 

FILE *fp1,*fp2; /*文件輸入、輸出指針*/ 


void input() /*從文件in.txt 輸入線路基本參數(shù)、節(jié)點數(shù)據(jù)、支路數(shù)據(jù)*/ 
{ 
int i,j,h; /*循環(huán)變量、節(jié)點（支路）類型*/ 

/**********************************打開電網(wǎng)參數(shù)文件**********************/ 
fp1=fopen("in.txt","r"); 
if (fp1==NULL) 
{ 
printf("Can not open file in.txt !\n"); 
exit(0); 
} 

/***********讀取節(jié)點數(shù),支路數(shù),互聯(lián)網(wǎng)數(shù),PQ節(jié)點數(shù)PV節(jié)點數(shù)和精度***********/ 
fscanf(fp1,"%d,%d,%d,%f\n",&n,&m,&qy,&eps); 

/************************讀取節(jié)點信息*************************************/ 
for(i=1;i<=n;i++) 
{ 
fscanf(fp1,"%d",&h); 
if(h==1) /*類型h=1是PQ節(jié)點*/ 
{ 
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節(jié)點*/ 
{ 
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是平衡節(jié)點*/ 
{ 
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++; 
} 
} 

/*******************************讀互聯(lián)網(wǎng)信息*************************************/ 
if(qy!=0) 
for(i=1;i<=qy;i++) /*輸入互聯(lián)網(wǎng)狀況*/ 
{ 
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**************************** 原始數(shù)據(jù)如下 ***********************************\n\n "); 
fprintf(fp2," 節(jié)點數(shù):%2d 支路數(shù):%2d 對地支路數(shù):%2d PQ節(jié)點數(shù):%2d PV節(jié)點數(shù):%2d 精度1:%f \n",n,m,dd,pq,pv,eps); 

fprintf(fp2,"\n------------------------------------------------------------------------------\n\n "); 
for(i=1;i<=pq;i++) 
fprintf(fp2," 節(jié)點%2d PQ節(jié)點 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," 節(jié)點%2d PV節(jié)點 P[%d]=%f V[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].v); 
fprintf(fp2," 節(jié)點%2d 平衡節(jié)點\n",i); 
fprintf(fp2,"\n-------------------------------------------------------------------------------\n\n "); 
for(i=1;i<=m;i++) 
{ 
if(zhilu[i].s==1) 
fprintf(fp2," 支路%2d：普通支路 相關節(jié)點:%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：變壓器支路 相關節(jié)點:%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： 對地支路 相關節(jié)點:%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," 互聯(lián)網(wǎng)區(qū)域:%2d 指定pv節(jié)點:%2d 規(guī)定有功功率:%8.5f 允許誤差:%8.5f 包括的節(jié)點: ",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() /*利用節(jié)點數(shù)據(jù)、支路數(shù)據(jù)形成新編號*/ 
{ 
int i,k,j,jd[N]; /*jd[N]記錄節(jié)點*/ 
struct jiedian tem; /*中間臨時節(jié)點所連支路數(shù)*/ 

for(i=1;i<=n;i++) /*對節(jié)點連接支路數(shù)賦0*/ 
{ 
jd[i]=0; 
} 

for(i=1;i<=n;i++) /*考慮每個節(jié)點*/ 
{ 
for(j=1;j<=m;j++) /*考慮每條支路*/ 
{ 
if((zhilu[j].p1==i||zhilu[j].p2==i)&&(zhilu[j].s!=3)) /*如果是這條支路的節(jié)點且不是接地支路*/ 
jd[i]++; 

} 
} 

for(i=1;i<pq;i++) /*對pq節(jié)點按所連支路的個數(shù)進行重新編號*/ 
{ 
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節(jié)點按所連支路的個數(shù)進行重新編號*/ 
{ 
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**************************** 節(jié)點優(yōu)化結果如下 *********************************\n\n "); 

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


void form_y() /*利用支路數(shù)據(jù)形成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++) /*節(jié)點導納矩陣的主對角線上的導納*/ 
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++) /*節(jié)點導納矩陣非主對角線上的導納*/ 
{ 
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****************************** 節(jié)點導納矩陣為 **********************************\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() /*利用節(jié)點數(shù)據(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節(jié)點子陣*/ 
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]); 
} 

}