電力潮流計算C源程序/**************************FLOW.C*********************************/

/*******************************************************************
* 這里提供的是電力系統潮流計算機解法的五個子程序,采用的方法是 *
* Newton_Raphson法. *
* 程序中所用的變量說明如下: *
* N:網絡節點總數. M:網絡的PQ節點數. *
* L:網絡的支路總數. N0:雅可比矩陣的行數. *
* N1:N0+1 K:打印開關.K=1,則打印;否則,不打印.*
* K1:子程序PLSC中判斷輸入電壓的形式.K1=1,則為極座標形式.否則*
* 為直角坐標形式. *
* D:有功及無功功率誤差的最大值. *
* G(I,J):Ybus的電導元素(實部). *
* B(I,J):Ybus的電納元素(虛部). *
* G1(I) :第I支路的串聯電導. B1(I):第I支路的串聯電納. *
* C1(I) :第I支路的pie型對稱接地電納. *
* C(I,J):第I節點J支路不對稱接地電納. *
* CO(I) :第I節點的接地電納. *
* S1(I) :第I節點的起始節點號. E1(I):第I節點的終止節點號. *
* P(I) :第I節點的注入有功功率. Q(I):第I節點的注入無功功率.*
* P0(I) :第I節點有功功率誤差. Q0(I):第I節點無功功率誤差. *
* V0(I) :第I節點(PV節點)的電壓誤差(平方誤差). *
* V(I) :第I節點的電壓誤差幅值. *
* E(I) :第I節點的電壓的實部. F(I):第I節點的電壓的虛部. *
* JM(I,J):Jacoby矩陣的第I行J列元素. *
* A(I,J):修正方程的增廣矩陣,三角化矩陣的第I行J列元素,運算結 *
* 束后A矩陣的最后一列存放修正的解. *
* P1(I) :第I支路由S1(I)節點注入的有功功率. *
* Q1(I) :第I支路由S1(I)節點注入的無功功率. *
* P2(I) :第I支路由E1(I)節點注入的有功功率. *
* Q2(I) :第I支路由E1(I)節點注入的無功功率. *
* P3(I) :第I支路的有功功率損耗. *
* Q3(I) :第I支路的無功功率損耗. *
* ANGLE(I):第I節點電壓的角度. *
*******************************************************************/
#include <math.h>
#include <stdio.h>

#define f1(i) (i-1) 
/* 把習慣的一階矩陣的下標轉化為C語言數組下標*/

#define f2(i,j,n) ((i-1)*(n)+j-1)
/* 把習慣的二階矩陣的下標轉化為C語言數組下標*/

/****************************************************
* 本子程序根據所給的支路導納及有關信息,形成結點 *
* 導納矩陣,如打印參數K=1,則輸出電導矩陣G和電納矩B *
****************************************************/
void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\
float *c,float *co,int k,int *s1,int *e1)
  {
    extern FILE *file4;
    FILE *fp;
    int i,j,io,i0;
    int pos1,pos2;
    int st,en;
    if(file4==NULL)
     {
       fp=stdout;
     }
    else
     {
       fp=file4; /* 輸出到文件 */
     } 

    /* 初始化矩陣G,B */
    for(i=1;i<=n;i++) 
     {
       for(j=1;j<=n;j++)
         {
           pos2=f2(i,j,n);
           g[pos2]=0;b[pos2]=0;
         }
     }

    /* 計算支路導納 */
    for(i=1;i<=l;i++)
     {
       /* 計算對角元 */ 
       pos1=f1(i);
       st=s1[pos1];en=e1[pos1];
       pos2=f2(st,st,n);
       g[pos2]+=g1[pos1];
       b[pos2]+=b1[pos1]+c1[pos1];
       pos2=f2(en,en,n);
       g[pos2]+=g1[pos1];
       b[pos2]+=b1[pos1]+c1[pos1];

       /* 計算非對角元 */ 
       pos2=f2(st,en,n);
       g[pos2]-=g1[pos1];
       b[pos2]-=b1[pos1];
       g[f2(en,st,n)]=g[f2(st,en,n)];
       b[f2(en,st,n)]=b[f2(st,en,n)];
     }

    /* 計算接地支路導納 */
    for(i=1;i<=n;i++)
     {
      /*　對稱部分　*/
      b[f2(i,i,n)]+=co[f1(i)];

      /* 非對稱部分　*/　
      for(j=1;j<=l;j++)
       {
        b[f2(i,i,n)]+=c[f2(i,j,l)];
       }
     }
    if(k!=1)
      {
        return; /* 如果K不為 1,則返回;否則,打印導納矩陣 */
      }
    fprintf(fp,"\n BUS ADMITTANCE MATRIX Y(BUS):");
    fprintf(fp,"\n ******************* ARRAY G ********************");
    for(io=1;io<=n;io+=5)
      {
       i0=(io+4)>n?n:(io+4);
       fprintf(fp,"\n");
       for(j=io;j<=i0;j++)
        {
         fprintf(fp,"%13d",j);
        }
       for(i=1;i<=n;i++)
        {
         fprintf(fp,"\n%2d",i);
         for(j=io;j<=i0;j++)
          {
           fprintf(fp,"%13.6f",g[f2(i,j,n)]);
          }
        }
       fprintf(fp,"\n");
      }

    fprintf(fp,"\n ******************* ARRAY B ********************");
    for(io=1;io<=n;io+=5)
     {
       i0=(io+4)>n?n:(io+4);
        fprintf(fp,"\n");
       for(j=io;j<=i0;j++)
        {
         fprintf(fp,"%13d",j);
        }
       for(i=1;i<=n;i++)
        {
          fprintf(fp,"\n%2d",i);
          for(j=io;j<=i0;j++)
          {
            fprintf(fp,"%13.6f",b[f2(i,j,n)]);
          }
        }
       fprintf(fp,"\n");
     }
    fprintf(fp,"\n************************************************");
  }

/*******************************************
* 本子程序根據所給的功率及電壓等數據 *
* 求出功率及電壓誤差量,并返回最大有功功率 *
* 以用于與給定誤差比較.如打印參數K=1,則輸 *
* 出P0,Q0(對PQ結點),V0(對PV結點). *
* 對應書上P185式(4-86)(4-87) *
*******************************************/

void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\
int n,float *e,float *f,int k,float *g,float *b,float *dd)
{
extern FILE *file4;
FILE *fp;
int i,j,l;
int pos1,pos2;
float a1,a2,d1,d;
if(file4==NULL)
{
fp=stdout; /* 輸出到屏幕 */
}
else
{
fp=file4; /* 輸出到文件　*/
}
l=n-1;
if(k==1)
{
fprintf(fp,"\n CHANGE OF P0,V**2,P0(I),Q0(I),V0(I) ");
fprintf(fp,"\n I P0(I) Q0(I)");
}
for(i=1;i<=l;i++)
{
a1=0;a2=0;
pos1=f1(i);
for(j=1;j<=n;j++)
{ 
/* a1,a2對應課本p185式(4-86)中括號內的式子 */
pos2=f2(i,j,n);
a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];
a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];
}
/* 計算式(4-86)(4-87)中的deltaPi　*/
p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;
if(i <= m)
{ /* 計算PQ結點中的deltaQi　*/
q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;
}
else
{ /* 計算PV結點中的deltaVi平方　*/
v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];
}

/* 輸出結果 */
if(k==1)
{
if(i<m)
{
fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);
}
else if(i==m)
{
fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);
fprintf(fp,"\n I P0(I) V0(I)");
}
else
{
fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]);
}
}
}

/* 找到deltaP和deltaQ中的最小者，作為收斂指標, 存在dd中 */
d=0;
for(i=1;i<=l;i++)
{
pos1=f1(i);
d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1];
if(d<d1)
{
d=d1;
}
if(i<=m)
{
d1=q0[pos1]>0?q0[pos1]:-q0[pos1];
if(d<d1)
{
d=d1;
}
}
} 
(*dd)=d;
}

/***************************************************
* 本子程序根據節點導納及電壓求Jacoby矩陣,用于求*
* 電壓修正量,如打印參數K=1,則輸出Jacoby矩陣. *
* 對應于課本P186式(4-89)(4-90) *
***************************************************/

void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k)
{
extern FILE *file4;
FILE *fp;
int i,j,i1,io,i0,ns;
int pos1,pos2;
if(file4==NULL)
{
fp=stdout;
}
else
{
fp=file4;
}

/* 初始化矩陣jm */
for(i=1;i<=n0;i++)
{
for(j=1;j<=n0;j++)
{
jm[f2(i,j,n0)]=0;
}
}

ns=n-1; /* 去掉一個平衡結點 */

/* 計算式(4-89)(4-90) */
for(i=1;i<=ns;i++)
{
/* 計算式(4-90) */
for(i1=1;i1<=n;i1++)
{
/* pos1是式(4-90)中的j */
pos1=f1(i1);

/* pos2是式(4-90)中的ij */
pos2=f2(i,i1,n);

if(i<=m) /* i是PQ結點 */
{
/* 計算式(4-90)中的Jii等式右側第一部分 */
jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1];

/* 計算式(4-90)中的Lii等式右側第一部分 */
jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];
}

/* 計算式(4-90)中的Hii等式右側第一部分 */
jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];

/* 計算式(4-90)中的Nii等式右側第一部分 */
jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1];
}

/* pos2是式(4-90)中的ii */
pos2=f2(i,i,n);

/* pos1是式(4-90)中的i */
pos1=f1(i);

if(i<=m) /* i是PQ結點 */　
{ 
/* 計算式(4-90)中的Jii */
jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];

/* 計算式(4-90)中的Lii */
jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[pos1];
}

/* 計算式(4-90)中的Hii */
jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];

/* 計算式(4-90)中的Jii */
jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];

if(i>m) /* PV結點 */
{
/* 計算式(4-90)中的Rii */
jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];

/* 計算式(4-90)中的Sii */
jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];
}

/* 計算式(4-89) */
for(j=1;j<=ns;j++)
{
if(j!=i)
{
/* pos1是式(4-89)中的i */
pos1=f1(i);

/* pos2是式(4-89)中的ij */
pos2=f2(i,j,n);

/* 計算式(4-89)中的Nij */
jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1];

/* 計算式(4-89)中的Hij */
jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1];

if(i<=m) /* i是PQ結點 */
{
/* 計算式(4-89)中的Lij (=-Hij) */
jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];

/* 計算式(4-89)中的Jij (=Nij) */
jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];
}
else /* i是PV結點 */
{
/* 計算式(4-89)中的Rij (=0) */
jm[f2(2*i-1,2*j-1,n0)]=0;

/* 計算式(4-89)中的Sij (=0) */
jm[f2(2*i-1,2*j,n0)]=0;
}
}
}
}
if(k!=1)
{
return;
}

/* 輸出Jacoby矩陣 */
fprintf(fp,"\n J MATRIX(C)");
for(io=1;io<=n0;io+=5)
{
i1=(io+4)>n0?n0:(io+4);
fprintf(fp,"\n");
for(j=io;j<=i1;j++)
{
fprintf(fp,"%10d",j);
}
for(i=1;i<=n0;i++)
{
fprintf(fp,"\n%2d",i);
for(j=io;j<=i1;j++)
{
fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);
}
}
}
fprintf(fp,"\n");
}
/**********************************************
* 本子程序用選列主元素的高斯消元法求解組 *
* 性方程組求各結點電壓修正量,如打印參數K=1,則*
* 輸出增廣矩陣變換中的上三角及電壓修正量.如果*
* 無唯一解,則給出信息,并停止程序運行. *
**********************************************/
void sevc ( float a[], int n0, int k, int n1)
{
extern FILE *file4;
FILE *fp;
int i,j,l,n2,n3,n4,i0,io,j1,i1;
float t0,t,c;
if(file4==NULL) fp=stdout;
else fp=file4;
for(i=1;i<=n0;i++)
{
l=i;
for(j=i;j<=n0;j++)
{
if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) )
{
l=j; /* 找到這行中的最大元 */
}
}
if(l!=i)
{ /* 行列交換 */
for (j=i;j<=n1;j++)
{
t=a[f2(i,j,n1)]; 
a[f2(i,j,n1)]=a[f2(l,j,n1)]; 
a[f2(l,j,n1)]=t;
}
}
if (fabs(a[f2(i,i,n1)]-0)<1e-10)
{ /* 對角元近似于0, 無解 */
printf("\nNo Solution\n"); 
exit (1);
}

t0=a[f2(i,i,n1)];
for(j=i;j<=n1;j++)
{
/* 除對角元 */
a[f2(i,j,n1)]/=t0;
}
if(i==n0)
{ /* 最后一行，不用消元 */
continue;
}

/* 消元 */
j1=i+1;
for(i1=j1;i1<=n0;i1++)
{
c=a[f2(i1,i,n1)];
for(j=i;j<=n1;j++)
{
a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;
}
}
}

if(k==1)
{ /* 輸出上三角矩陣 */
fprintf(fp,"\nTrianglar Angmentex Matrix ");
for(io=1;io<=n1;io+=5)
{
i0=(io+4)>n1?n1:(io+4);
fprintf(fp,"\n");
fprintf(fp," ");
for(i=io;i<=i0;i++)
{
fprintf(fp,"%12d",i);
}
for(i=1;i<=n0;i++)
{
fprintf(fp,"\n");
fprintf(fp,"%2d",i);
for(j=io;j<=i0;j++)
{
fprintf(fp,"%15.6f", a[f2(i,j,n1)]);
}
}
}
}

/* 回代求方程解 */
n2=n1-2;
for(i=1;i<=n2;i++)
{
n3=n1-i;
for(i1=n3;i1<=n0;i1++)
{
n4=n0-i;
a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];
}
}

if(k!=1)
{
return;
}

/* 輸出電壓修正值 */
fprintf(fp,"\nVoltage correction E(i), F(i) :");
for(io=1;io<=n0;io+=4)
{ 
i1=(io+1)/2;
i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);
fprintf(fp,"\n");
for(j=i1;j<=i0;j++)
{
fprintf(fp,"%16d%16d",j,j);
}
i1 = 2*i0;
fprintf(fp,"\n");
for(i=io;i<=i1;i++)
{
fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);
}
}
}


#define Pi 3.1415927/180
void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\
int e1[],int s1[],float g1[],float b1[],float c1[],float c[],\
float co[],float p1[],float q1[],float p2[],float q2[],float p3[],\
float q3[],float p[],float q[],float v[],float angle[],int k1)
{
extern FILE *file4;
FILE *fp;
float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;
int i,i1,j,m1,ns,pos1,pos2,km;
ns=n-1;
if(file4==NULL)
{
fp=stdout;
}
else
{
fp=file4;
}

fprintf(fp,"\nTHE RESULT ARE:");
if(k1==1)
{
for(i=0;i<n;i++)
{
angle[i]*=Pi;
e[i]=v[i]*cos(angle[i]);
f[i]=v[i]*sin(angle[i]);
}
}
t1=0.0;t2=0.0;
for(i=1;i<=n;i++)
{
pos1=f1(i);pos2=f2(n,i,n);
t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];
t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];
}
pos1=f1(n);
p[pos1]=t1*e[pos1];
q[pos1]=-t2*e[pos1];
m1=m+1;
for(i1=m1;i1<=ns;i1++)
{
t1=0;t2=0;
for(i=1;i<=n;i++)
{
pos1=f1(i);pos2=f2(i1,i,n);
t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];
t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];
}
pos1=f1(i1);
q[pos1]=f[pos1]*t1-e[pos1]*t2;
}
for(i=0;i<n; i++)
{
cm=co[i];
if(cm!=0)
{
q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;
}
}
fprintf(fp,"\nBUS DATA");
fprintf(fp,"\nBUS VOLTAGE ANGLE(DEGS.) BUS P BUS Q");
for(i=0;i<n;i++)
{
v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);
x=e[i];
y=f[i];
z=y/x;
angle[i]=atan(z);
angle[i]/=Pi;
fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[i],p[i],q[i]);
}
fprintf(fp,"\n LINE FLOW ");
for(i=1;i<=l;i++)
{
pos1=f1(i);
st=s1[pos1];
en=e1[pos1];
x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];
x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];
y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];
y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];
p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];
q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];
p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];
q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];
for(j=1;j<=n;j++)
{
cm=c[f2(j,i,l)];
if(cm!=0.0)
{
km=1;
if(en==j)
{
km=2;
}
if(km==1)
{
q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;
}
else
{
q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;
}
}
}
p3[pos1]=p1[pos1]+p2[pos1] ;
q3[pos1]=q1[pos1]+q2[pos1] ;
fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%13.6e%13.6e",\
i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1],\
e1[pos1],s1[pos1],p2[pos1],q2[pos1]);
}
}