潮流.c

基于c语言的电力系统潮流计算程序（牛顿拉夫逊算法）

#include <stdlib.h>
#include <math.h> 
#include <stdio.h> 
#define N 6
#define ep 1e-5


/*矩阵求逆*/

int brinv(double a[], int n) 
{ int *is,*js,i,j,k,l,u,v; 
double d,p; 
is=malloc(n*sizeof(int)); 
js=malloc(n*sizeof(int)); 
for (k=0; k<=n-1; k++) 
{ d=0.0; 
for (i=k; i<=n-1; i++) 
for (j=k; j<=n-1; j++) 
{ l=i*n+j; p=fabs(a[l]); 
if (p>d) { d=p; is[k]=i; js[k]=j;} 
} 
if (d+1.0==1.0) 
{ free(is); free(js); printf("err**not inv\n"); 
return(0); 
} 
if (is[k]!=k) 
for (j=0; j<=n-1; j++) 
{ u=k*n+j; v=is[k]*n+j; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
if (js[k]!=k) 
for (i=0; i<=n-1; i++) 
{ u=i*n+k; v=i*n+js[k]; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
l=k*n+k; 
a[l]=1.0/a[l]; 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ u=k*n+j; a[u]=a[u]*a[l];} 
for (i=0; i<=n-1; i++) 
if (i!=k) 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ u=i*n+j; 
a[u]=a[u]-a[i*n+k]*a[k*n+j]; 
} 
for (i=0; i<=n-1; i++) 
if (i!=k) 
{ u=i*n+k; a[u]=-a[u]*a[l];} 
} 
for (k=n-1; k>=0; k--) 
{ if (js[k]!=k) 
for (j=0; j<=n-1; j++) 
{ u=k*n+j; v=js[k]*n+j; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
if (is[k]!=k) 
for (i=0; i<=n-1; i++) 
{ u=i*n+k; v=i*n+is[k]; 
p=a[u]; a[u]=a[v]; a[v]=p; 
} 
} 
free(is); free(js); 
return(1); 
} 



/*矩阵相乘*/
void matrixmul(double g[N][N],double h[N],double xx[N])


 {  int i,j; 
   xx[0]=0;xx[1]=0;xx[2]=0;xx[3]=0;xx[4]=0;xx[5]=0;
   for(i=0;i<=N-1;i++)
     for(j=0;j<=N-1;j++)
        xx[i]+=g[i][j]*h[j];
    
 }


/* fe矩阵 */

void funfe(double fe[N],double x[N])
 {fe[0]=fe[0]+x[0];    fe[1]=fe[1]+x[1];
  fe[2]=fe[2]+x[2];    fe[3]=fe[3]+x[3];
  fe[4]=fe[4]+x[4];    fe[5]=fe[5]+x[5];
 }


/* pq矩阵  */

void matrixpq(double pq[N],double G[4][4],double B[4][4],double e[4],double f[4])
{int i,j;
 double pp,qq;
 double p[3]={-0.3,-0.55,0.5}; 
 double q[2]={-0.18,-0.13};
 double u=1.1;

for(i=0;i<3;i++)
 { pp=0;qq=0;
   for(j=0;j<4;j++)
    { pp=pp+e[i]*(G[i][j]*e[j]-B[i][j]*f[j])+f[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
      if (i<2)
        qq=qq+f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
    }
   pq[2*i]=p[i]-pp;
   if (i<2)
        pq[2*i+1]=q[i]-qq;
   else
        pq[5]=u*u-(e[2]*e[2]+f[2]*f[2]);
   
  }

  
}


/*  e,f矩阵  */

void matrixef(double e[4],double f[4],double fe[N])
 {  e[0]=fe[1];    f[0]=fe[0];
    e[1]=fe[3];    f[1]=fe[2];
    e[2]=fe[5];    f[2]=fe[4];
    e[3]=1.05;     f[3]=0;
 }

/*  雅克比矩阵  */
void matrixjcb(double G[4][4],double B[4][4],double e[4],double f[4],double jcbjz[N][N])
 {int i,j,k;
  double bb;double aa;
  double H[3][3],M[3][3],J[2][3],L[2][3];
  for(i=0;i<3;i++)
    { for(j=0;j<3;j++)
      { if (i!=j)
         {H[i][j]=G[i][j]*f[i]-B[i][j]*e[i];
          M[i][j]=G[i][j]*e[i]+B[i][j]*f[i];
          if (i!=2)
             {J[i][j]=-M[i][j];
              L[i][j]=H[i][j];
             }
         }
        else 
         { bb=0;aa=0;
           for(k=0;k<4;k++)
            
               {bb=bb+G[i][k]*f[k]+B[i][k]*e[k];
                aa=aa+G[i][k]*e[k]-B[i][j]*f[k];
               }
             
             
          H[i][i]=G[i][i]*f[i]-B[i][i]*e[i]+bb;  
          M[i][i]=G[i][i]*e[i]+B[i][i]*f[i]+aa;
          if (i!=2)
            {J[i][i]=(-G[i][i])*e[i]-B[i][i]*f[i]+aa; 
             L[i][i]=(-B[i][i])*e[i]+G[i][i]*f[i]-bb;
            }
         }
       } 
    }    
 jcbjz[0][0]=H[0][0]; jcbjz[0][1]=M[0][0]; jcbjz[0][2]=H[0][1];
 jcbjz[0][3]=M[0][1]; jcbjz[0][4]=H[0][2]; jcbjz[0][5]=M[0][2];

 jcbjz[1][0]=J[0][0]; jcbjz[1][1]=L[0][0]; jcbjz[1][2]=J[0][1];
 jcbjz[1][3]=L[0][1]; jcbjz[1][4]=J[0][2]; jcbjz[1][5]=L[0][2];

 jcbjz[2][0]=H[1][0]; jcbjz[2][1]=M[1][0]; jcbjz[2][2]=H[1][1];
 jcbjz[2][3]=M[1][1]; jcbjz[2][4]=H[1][2]; jcbjz[2][5]=M[1][2];

 jcbjz[3][0]=J[1][0]; jcbjz[3][1]=L[1][0]; jcbjz[3][2]=J[1][1];
 jcbjz[3][3]=L[1][1]; jcbjz[3][4]=J[1][2]; jcbjz[3][5]=L[1][2];
 
 jcbjz[4][0]=H[2][0]; jcbjz[4][1]=M[2][0]; jcbjz[4][2]=H[2][1];
 jcbjz[4][3]=M[2][1]; jcbjz[4][4]=H[2][2]; jcbjz[4][5]=M[2][2];

 jcbjz[5][0]=0;       jcbjz[5][1]=0;       jcbjz[5][2]=0;
 jcbjz[5][3]=0;       jcbjz[5][4]=2*f[3];  jcbjz[5][5]=2*e[3];    



}




/*主函数*/

int main()
{ int i,j,h,num;
  double fe[N]={0,1,0,1,0,1};
  double jcbjz[N][N],b[N][N],pq[N],x[N],e[4],f[4];
  double G[4][4]={{1.0421,-0.5882,0,-0.4539},
                  {-0.5822,1.069,0,-0.4808},
                  {0,0,0,0},
                  {-0.4539,-0.4808,0,0.9346}};
 double B[4][4]={{-8.2429,2.3529,3.6666,1.8911},
                  {2.3529,-4.7274,0,2.4038},
                    {3.6666,0,-3.3333,0},
                        {1.8911,2.4038,0,-4.2616}};
   num=0;


do{
  h=0;
  matrixef(e,f,fe);  
  matrixjcb(G,B,e,f,jcbjz);

/*  for(i=0;i<6;i++)
  {for(j=0;j<6;j++)
     printf("%6.4f ",jcbjz[i][j]);
   printf("\n");
  }
*/
  matrixpq(pq,G,B,e,f); 
 for (i=0; i<=N-1; i++)
   for (j=0; j<=N-1; j++)
     b[i][j]=jcbjz[i][j];
i=0;
i=brinv(b,N);
if (i!=0)
 { num++;
  printf("\n%d",num);

printf("\n MAX pq :");
for(i=0;i<6;i++)
 printf("%6.5e ",pq[i]);
matrixmul(b,pq,x);
  for(i=0;i<=N-1;i++)
    if (fabs(x[i])<ep)
     h++;
  funfe(fe,x);
 printf("\nMAT X-fe IS:");
  for(i=0;i<=N-1;i++)
    printf("%6.5f ",x[i]);
 printf("\nMAT fe IS:");
  for(i=0;i<=N-1;i++)
    printf("%6.5f ",fe[i]);
 
  }
 
}
while(h!=6);
printf("\nnum is:  %d\n",num);
for(i=0;i<3;i++)
  printf("i=%d,Ui=%f,Ai=%f\n",i+1,sqrt(e[i]*e[i]+f[i]*f[i]),atan(f[i]/e[i])*180/3.14159);  
   
 getch();
}