电力潮流计算.txt

这是电力系统中常常会用到的潮流的分析和计算的一个程序,读者可以编译生成软件试用

// 潮流计算程序.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"


#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define N 9       /*总结点数*/
#define M 6       /*PQ结点数*/
#define K 9       /*线路数*/
#define eps 1e-4
#define double float
void guass(int n,int m,double c[],double b[][N],double x[])       /*高斯函数*/
{
double 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;
      double R;
      double X;
      double b;
}line[K];
struct line *t;
void main()
{
      double r[N]={0.0};         
      double u[N]={1.04,1.025,1.025,1.0,1.0,1.0,1.0,1.0,1.0};   
      double p[N]={1,1.63,0.85,0.0,0.0,0.0,-1.25,-0.9,-1.0};    
      double q[N]={1,1,1,0.0,0.0,0.0,-0.5,-0.3,-0.35};      
      double g[N][N]={0.0},b[N][N]={0.0};       
      double h[N][N]={0.0};
      double B[N][N]={0.0};
      double temp;
      double H[N][6];
      double lr,lx,lb1,lg,lb;
      int i,j;
      int ku=0,kr=0,kp=1,kq=1;
      void val(double u[N],double g[N][N],double b[N][N],double r[N],int ku, int kr,double h[N][N]);       /*函数申明*/

      FILE *fp;
     if ((fp=fopen("data.txt","r"))==NULL)
      {
       printf("cannot open this file.\n");
       exit(0);
      }
      //fp=fopen("data.txt","r");
      for(i=0;i<N;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;    
         }
      getchar();
      printf("\n=====节点导纳矩阵=====\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");
      getchar();
      printf("\n=====给定初值=====\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=====迭代求解=====\n");              
      while(kp==1)
       {
        double ip,iq,max;
         double dp[N],dq[N],dr[N];
         double dpu[N],dqu[N];
         double y1[N-1],y2[M];
        double 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");
            }
         getchar();
         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===== 第 %d 次迭代后电压      相角初值 =====\n",kr+1);
           for(i=1;i<N;i++)
           printf("r[%d]=%8f\n",i+1,r[i]);
           getchar();
           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=====第 %d 次迭代电压初值=====\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(double u[N],double g[N][N],double b[N][N],double r[N],int ku, int kr,double h[N][N])
       {
         double ps=0,pv1=0,pv2=0;
        double qs=0,qv1=0,qv2=0;
         double p[N][N]={0};
         double q[N][N]={0};
         double s[N][N]={0};
        double dp[N][N]={0};
         double dq[N][N]={0};
         double ds[N][N]={0};
         double dSp=0,dSq=0;
         int i,j;
         FILE *fp1;
         printf("\n=====平衡节点功率 =====\n");
         getchar();
         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 节点功率=====\n");
         getchar();
         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======迭代后电压和相交值 ======\n");
         getchar();
         for(i=0;i<N;i++)
         printf("u[%d]=%8f       r[%d]=%8f\n",i+1,u[i],i+1,r[i]);
         printf("\n=====线路功率=======\n");
         for(i=0;i<N;i++)
         { 
           getchar();
           { 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=====线路损耗=====\n");
        for(i=0;i<N;i++)
        { 
          getchar();
          { 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=====网络中损耗=====\n");
        getchar();
        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======迭代次数========\n");
        printf("ku=%d\n",ku);
        printf("kr=%d\n",kr);
    
        printf("\n=======数据保存=====\n");
        fp1=fopen("jieguo.txt","w+");
        {
          fprintf(fp1,"线路潮流:\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,"支路损耗:\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,"网络中损耗:\n");
          fprintf(fp1,"dS=%8f+j(%8f)\n",dSp,dSq);
          fprintf(fp1,"\n");
          fprintf(fp1,"PV 节点功率:\n");
          fprintf(fp1,"sv1=%8f+j(%8f)\n",pv1,qv1);
          fprintf(fp1,"sv2=%8f+j(%8f)\n",pv2,qv2);
          fprintf(fp1,"\n");
          fprintf(fp1,"平衡节点功率:\n");
          fprintf(fp1,"sp=%8f+j(%8f)\n",ps,qs);
          fprintf(fp1,"\n");
          fprintf(fp1,"节点电压和相交值:\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===========结束==============\n");
        getchar();

}