pq.c

9节点PQ潮流计算程序

#include<stdio.h>
#include<math.h>
#define N 9         /*总结点数*/
#define M 6         /*PQ结点数*/
#define K 9         /*线路数*/
#define eps 1e-4
void guass(int n,int m,float c[],float b[][N],float x[])         /*高斯函数*/
{
float a[N][N],y[N];
int i,j,k;

for(i=0;i<m;i++)
{
        y[i]=c[i];
        for(j=0;j<m;j++) a[i][j]=b[i+n][j+n];
}

for(k=0;k<m-1;k++)
        for(i=k+1;i<m;i++)
        {
         for(j=k+1;j<m;j++) a[i][j]=a[i][j]-a[k][j]*a[i][k]/a[k][k];
         y[i]=y[i]-y[k]*a[i][k]/a[k][k];
        }

for(i=m-1;i>=0;i--)
{
        for(j=i+1;j<m;j++) y[i]-=a[i][j]*x[j];
        x[i]=y[i]/a[i][i];
}
}

struct line
{int Lindex;
        int Headnode;
        int Endnode;
        float R;
        float X;
        float b;
}line[K];
struct line *t;
main()
{
        float r[N]={0.0};         
        float u[N]={1.04,1.025,1.025,1.0,1.0,1.0,1.0,1.0,1.0};   
        float p[N]={1,1.63,0.85,0.0,0.0,0.0,-1.25,-0.9,-1.0};    
        float q[N]={1,1,1,0.0,0.0,0.0,-0.5,-0.3,-0.35};      
        float g[N][N]={0.0},b[N][N]={0.0};       
        float h[N][N]={0.0};
        float B[N][N]={0.0};
        float temp;
        float H[K][6];
        float lr,lx,lb1,lg,lb;
        int i,j;
        int ku=0,kr=0,kp=1,kq=1;
        void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]);         /*函数

申明*/

        FILE *fp;
        fp=fopen("data.txt","r");
        for(i=0;i<K;i++)
        { 
          for(j=0;j<6;j++)
           { fscanf(fp,"%f,",&temp);
             H[i][j]=temp;
           }
         }
        fclose(fp);

        for(i=0;i<K;i++)
          {
            line[i].Lindex=(int)H[i][0];line[i].Headnode=(int)H[i][1];line[i].Endnode=(int)H[i][2];
            line[i].R=H[i][3];line[i].X=H[i][4];line[i].b=H[i][5];}
        for(t=line;t<line+K;t++)
          {i=t->Headnode-1;
           j=t->Endnode-1;
           lr=t->R;
           lx=t->X;
           lb1=t->b;   
           lg=lr/(lr*lr+lx*lx);   
           lb=-lx/(lr*lr+lx*lx);   
           g[i][i]+=lg;
           g[j][j]+=lg;   
           b[i][i]+=lb+lb1;
           b[j][j]+=lb+lb1;    
           h[i][j]=h[j][i]=-lb1;   
           g[i][j]-=lg;
           g[j][i]-=lg;    
           b[i][j]-=lb;
           b[j][i]-=lb;    
           }
        getch();
        printf("\n=====jie dian dao na ju zhen=====\n");
        for(i=0;i<N;i++)
          {
            for(j=0;j<N;j++)
             {
               B[i][j]=b[i][j];
               printf("%8f,",B[i][j]); 
              }
           printf("\n");
           }

        printf("\n");
        getch();
        printf("\n=====gei ding chu zhi=====\n");
        for(i=0;i<N;i++)
          printf("u[%d]=%8f        p[%d]=%8f        q[%d]=%8f\n",i+1,u[i],i+1,p[i],i+1,q[i]);       
        printf("\n=====die dai qiu jie=====\n");              
        while(kp==1)
         {
           float ip,iq,max;
           float dp[N],dq[N],dr[N];
           float dpu[N],dqu[N];
           float y1[N-1],y2[M];
           float x1[N-1],x2[M];

           for(i=1;i<N;i++) /*算dp对应B',N-1维,除去平衡结点*/
            {
              ip=0;
              for(j=0;j<N;j++)   
              ip=ip+u[j]*(g[i][j]*cos(r[i]-r[j])+b[i][j]*sin(r[i]-r[j]));
              dp[i]=p[i]-u[i]*ip;
              dpu[i]=dp[i]/u[i];
              printf("dp[%d]=%8f\n",i+1,dp[i]);
              printf("\n");
              }
           getch();
           max=fabs(dpu[1]);
           for(i=1;i<N;i++)
             {
              if(fabs(dpu[i])>max)
              max=fabs(dpu[i]);
             }
         if (max>=eps)
         {
             for(i=0;i<N-1;i++)
             y1[i]=-dpu[i+1];        /*起值不同,为了对应,故加一*/
    
             guass(1,N-1,y1,B,x1);
             for (i=1;i<N;i++)
              {
               dr[i]=x1[i-1]/u[i];
                r[i]=r[i]+dr[i];
               }
             printf("\n===== di %d ci die dai hou dian ya        xiang jiao chu zhi =====\n",kr+1);
             for(i=1;i<N;i++)
             printf("r[%d]=%8f\n",i+1,r[i]);
             getch();
             printf("\n\n");
             kr=kr+1;kq=1;
          top:
            for(i=N-M;i<N;i++) /*算dq对应B",仅M维,除去平衡结点和PV结点*/
             {
              iq=0;
              for(j=0;j<N;j++)
                iq=iq+u[j]*(g[i][j]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r[j]));
              dq[i]=q[i]-u[i]*iq;
              dqu[i]=dq[i]/u[i];
              printf("dq[%d]=%8f\n",i+1,dq[i]);
              printf("\n");
              }
            max=fabs(dqu[N-M]);
            for (i=N-M;i<N;i++)
             {
               if(fabs(dqu[i])>max)
               max=fabs(dqu[i]);
              }
            if(max>=eps)
              {
               for(i=0;i<M;i++)
               y2[i]=-dqu[i+N-M];        /*同上,对应加N-M*/
  
               guass(N-M,M,y2,B,x2);
               for(i=N-M;i<N;i++)
                 u[i]=u[i]+x2[i-(N-M)];
               printf("\n=====di %d ci die dai dian ya chu zhi=====\n",ku+1);
               for(i=N-M;i<N;i++)
                 printf("u[%d]=%8f\n",i+1,u[i]);
               printf("\n\n");
               ku=ku+1;kp=1;
               }

             else
               {
                kq=0;
                if(kp==0)
                val(u,g,b,r,ku,kr,h);
                }

            }
      
            else
              {
               kp=0;
               if(kq==0)
                 val(u,g,b,r,ku,kr,h);
               else
                 goto top;
               }

           }

}
void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N])
         {
           float ps=0,pv1=0,pv2=0;
           float qs=0,qv1=0,qv2=0;
           float p[N][N]={0};
           float q[N][N]={0};
           float s[N][N];
           float dp[N][N]={0};
           float dq[N][N]={0};
           float ds[N][N];
           float dSp=0,dSq=0;
           int i,j;
           FILE *fp1;
           printf("\n=====ping heng jie dian gong lv =====\n");
           getch();
           for(i=0;i<N;i++)
           {
              ps=ps+u[0]*u[i]*(g[0][i]*cos(r[0]-r[i])+b[0][i]*sin(r[0]-r[i]));
              qs=qs+u[0]*u[i]*(g[0][i]*sin(r[0]-r[i])-b[0][i]*cos(r[0]-r[i]));
           }
           printf("sp=%8f+j(%8f)\n",ps,qs);
           printf("\n=====PV jie dian gong lv=====\n");
           getch();
           for(i=0;i<N;i++)
           {
              pv1=pv1+u[1]*u[i]*(g[1][i]*cos(r[1]-r[i])+b[1][i]*sin(r[1]-r[i]));
              qv1=qv1+u[1]*u[i]*(g[1][i]*sin(r[1]-r[i])-b[1][i]*cos(r[1]-r[i]));
           }
           printf("sv1=%8f+j(%8f)\n",pv1,qv1);
           for(i=0;i<N;i++)
           {
              pv2=pv2+u[2]*u[i]*(g[2][i]*cos(r[2]-r[i])+b[2][i]*sin(r[2]-r[i]));
              qv2=qv2+u[2]*u[i]*(g[2][i]*sin(r[2]-r[i])-b[2][i]*cos(r[2]-r[i]));
           }
           printf("sv2=%8f+j(%8f)\n",pv2,qv2);
           for(i=0;i<N;i++)
             for(j=0;j<N;j++)
             {
               p[i][j]=u[i]*u[i]*(-g[i][j])+u[i]*u[j]*(g[i][j]*cos(r[i]-r[j])+b[i][j]*sin(r[i]-r[j]));
               q[i][j]=u[i]*u[i]*(-h[i][j]+b[i][j])+u[i]*u[j]*(g[i][j]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r[j]));
               dp[i][j]=u[i]*u[i]*(-g[i][j])+u[i]*u[j]*(g[i][j]*cos(r[i]-r[j])+b[i][j]*sin(r[i]-r[j]))
           +u[j]*u[j]*(-g[j][i])+u[j]*u[i]*(g[j][i]*cos(r[j]-r[i])+b[j][i]*sin(r[j]-r[i]));
               dq[i][j]=u[i]*u[i]*(-h[i][j]+b[i][j])+u[i]*u[j]*(g[i][j]*sin(r[i]-r[j])-b[i][j]*cos(r[i]-r[j]))
           +u[j]*u[j]*(-h[j][i]+b[j][i])+u[j]*u[i]*(g[j][i]*sin(r[j]-r[i])-b[j][i]*cos(r[j]-r[i]));
             }
           printf("\n======die dai hou dian ya yu xiang jiao zhi ======\n");
           getch();
           for(i=0;i<N;i++)
           printf("u[%d]=%8f         r[%d]=%8f\n",i+1,u[i],i+1,r[i]);
           printf("\n=====xian lu gong lv=======\n");
           for(i=0;i<N;i++)
           { 
             getch();
             { for(j=0;j<N;j++)
         { printf("s[%d][%d]=%8f+j(%8f)\n",i+1,j+1,p[i][j],q[i][j]);
                  printf("\n");
                 }
               printf("\n");
             }
          }
          printf("\n");
          printf("\n=====xian lu sun hao=====\n");
          for(i=0;i<N;i++)
          { 
            getch();
            { for(j=i+1;j<N;j++)
         { printf("ds[%d][%d]=%8f+j%8f\n",i+1,j+1,dp[i][j],dq[i][j]);
                printf("\n");
                }
               printf("\n");
              }
          }
          printf("\n");
          printf("\n=====wang luo zong sun hao=====\n");
          getch();
          for(i=0;i<N;i++)
           {
             for(j=i+1;j<N;j++)
             { dSp+=dp[i][j];
               dSq+=dq[i][j];
               }
            }
          printf("dS=%8f+j(%8f)\n",dSp,dSq);  
          printf("\n======die dai ci shu========\n");
          printf("ku=%d\n",ku);
          printf("kr=%d\n",kr);
    
          printf("\n=======shu ju bao cun=====\n");
          fp1=fopen("jieguo.txt","w+");
          {
            fprintf(fp1,"xian lu cao liu:\n");
            for(i=0;i<N;i++)
             { 
               for(j=0;j<N;j++)
         fprintf(fp1,"s[%d][%d]=%8f+j%8f\n",i+1,j+1,p[i][j],q[i][j]);
               }
            fprintf(fp1,"\n");
            fprintf(fp1,"zhi lu sun hao:\n");
            for(i=0;i<N;i++)
             { 
               for(j=i+1;j<N;j++)
         fprintf(fp1,"ds[%d][%d]=%8f+j%8f\n",i+1,j+1,dp[i][j],dq[i][j]);
               }
            fprintf(fp1,"\n");
            fprintf(fp1,"wang luo zong sun hao:\n");
            fprintf(fp1,"dS=%8f+j(%8f)\n",dSp,dSq);
            fprintf(fp1,"\n");
            fprintf(fp1,"PV jie dian gong lv:\n");
            fprintf(fp1,"sv1=%8f+j(%8f)\n",pv1,qv1);
            fprintf(fp1,"sv2=%8f+j(%8f)\n",pv2,qv2);
            fprintf(fp1,"\n");
            fprintf(fp1,"ping heng jie dian gong lv:\n");
            fprintf(fp1,"sp=%8f+j(%8f)\n",ps,qs);
            fprintf(fp1,"\n");
            fprintf(fp1,"jie dian dian ya yu xiang jiao:\n");
            for(i=0;i<N;i++)
            fprintf(fp1,"u[%d]=%8f         r[%d]=%8f\n",i+1,u[i],i+1,r[i]);
          }
          printf("\n===========THE END==============\n");
          getch();

}