power flow2.cpp

电力系统的牛顿-拉夫逊潮流计算。有程序和PPT

// power flow2.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include"math.h"
#include"stdio.h"
#include"stdlib.h"
#include "Conio.h"
#define N 4    /*定义节点数*/
#define M 2    /*定义PQ节点数*/
#define V  N-M-1   /*定义PV节点数*/

#define N 4    /*定义节点数*/
#define M 2    /*定义PQ节点数*/
#define V  N-M-1   /*定义PV节点数*/
#define n 6
#define precision 1e-16 

int gauss(double a[][n+2],double x[n+1])
{int i,j,k,r;double c;
   for(k=1;k<=n-1;k++)
    {r=k;
     for(i=k;i<=n;i++)
      if(fabs(a[i][k])>fabs(a[r][k]))r=i;
      if(fabs(a[r][k])<precision)
      {printf("\n det A=0.failed!");exit(0);}
       if(r!=k)
       {for(j=k;j<=n+1;j++)
       {c=a[k][j];a[k][j]=a[r][j];a[r][j]=c;}
       }
        for(i=k+1;i<=n;i++)
         {c=a[i][k]/a[k][k];
        for(j=k+1;j<=n+1;j++)
        a[i][j]=a[i][j]-c*a[k][j];
        } 
       } 
      if(fabs(a[n][n])<precision)  
        {printf("\n det A=0.failed!");exit(0);}
      for(k=n;k>=1;k--)
        {x[k]=a[k][n+1];
         for(j=k+1;j<=n;j++)
         x[k]=x[k]-a[k][j]*x[j];
         x[k]=x[k]/a[k][k];
          }       
        return(1);
        }
        







int _tmain(int argc, _TCHAR* argv[])
{
 
  int k=0,i,j;   /*k为迭代次数,n为节点个数*/
  double temp1,temp2;
   double g[N][N]={{1.0421,-0.5882,0.0,-0.4539},{-0.5822,1.069,0.0,-0.4808},
  {0.0,0.0,0.0,0.0},{-0.4539,-0.4808,0.0,0.9346}},b[N][N]={{-8.2429,2.3529,3.6666,1.8911},
  {2.3529,-4.7274,0.0,2.4038},{3.6666,0.0,-3.3333,0.0},{1.8911,2.4038,0.0,-4.2616}};
  double e[N]={1.0,1.0,1.0,1.05},f[N]={0.0,0.0,0.0,0.0};
  double p[N]={-0.3,-0.55,0.5,0.0},q[N]={-0.18,-0.13,0.0,0.0};
  double u[N]={0.0,0.0,1.1,1.05};
  double dp[N]={0.0},dq[N]={0.0};
  double du2[N]={0.0};
  double jac[2*(N-1)][2*(N-1)];
  double h[N-1][N-1],nn[N-1][N-1],jj[M][N-1],l[M][N-1],r[N-1][N-1],s[N-1][N-1];
  double dpq[2*(N-1)],dfe[2*(N-1)];
  double fe[2*(N-1)]={0.0,1.0,0.0,1.0,0.0,1.0};
  int re,det;

  double aa[n][n+1],a[n+1][n+2],x[n+1];

  do
  { for(i=0;i<N-1;i++)                  /*计算节点功率,电压不平衡量*/
     {  temp1=0;
        temp2=0;
        for(j=0;j<N;j++)
        {
           temp1+=(g[i][j]*e[j]-b[i][j]*f[j]);
           temp2+=(g[i][j]*f[j]+b[i][j]*e[j]);

        }
        dp[i]=p[i]-e[i]*temp1-f[i]*temp2;
      }
    for(i=0;i<N-2;i++)
      {  temp1=0;
         temp2=0;
         for(j=0;j<N;j++)
         {
            temp1+=(g[i][j]*e[j]-b[i][j]*f[j]);
            temp2+=(g[i][j]*f[j]+b[i][j]*e[j]);
         }
         dq[i]=q[i]-f[i]*temp1+e[i]*temp2;
       }

    for(i=M;i<N-1;i++)
       {  du2[i]=u[i]*u[i]-(e[i]*e[i]+f[i]*f[i]);
        }

    for(i=0;i<N-1;i++)                //显示dp
      printf("dp(%d)=%f\n",i+1,dp[i]);
	for(i=0;i<M;i++)                //显示dq
      printf("dq(%d)=%f\n",i+1,dq[i]);
    for(i=M;i<(M+V);i++)                //显示du2
      printf("du2(%d)=%f\n",i+1,du2[i]);

  
	   for(i=0;i<N-1;i++)                 /*计算雅克比矩阵各元素*/
       { temp1=0;
         temp2=0;
	       for(j=0;j<N;j++)
             {temp1+=(g[i][j]*f[j]+b[i][j]*e[j]);
		      temp2+=(g[i][j]*e[j]-b[i][j]*f[j]);
            }
		   
		   
		   for(j=0;j<N-1;j++)            /*计算H,N元素*/
          {if(j!=i)
           { h[i][j]=(-1)*b[i][j]*e[i]+g[i][j]*f[i];
             nn[i][j]=g[i][j]*e[i]+b[i][j]*f[i];

           }

          if(j==i)
           { h[i][i]=(-1)*b[i][i]*e[i]+g[i][i]*f[i]+temp1;
             nn[i][i]=g[i][i]*e[i]+b[i][i]*f[i]+temp2;

           }
          }
        }

     

     for(i=0;i<M;i++)                    /*计算J,L元素*/
        { temp1=0;
          temp2=0;
	       for(j=0;j<N;j++)
             {temp1+=(g[i][j]*f[j]+b[i][j]*e[j]);
		      temp2+=(g[i][j]*e[j]-b[i][j]*f[j]);
            }

			for(j=0;j<N-1;j++)
          {if(j!=i)
            { jj[i][j]=(-1)*nn[i][j];
               l[i][j]=h[i][j];
            
             }
            if(j==i)
             { jj[i][i]=(-1)*g[i][i]*e[i]-b[i][i]*f[i]+temp2;
                l[i][i]=(-1)*b[i][i]*e[i]+g[i][i]*f[i]-temp1;
              
             }
           }
         }
	 for(i=M;i<N-1;i++)                   //计算R,S元素
	 {for(j=0;j<N-1;j++)
      {if(j!=i)
            { 
             r[i][j]=0;
             s[i][j]=0;
             }
            if(j==i)
             { 
                r[i][i]=2*f[i];
                s[i][i]=2*e[i];
             }
           }

	 }

       for(i=0;i<N-1;i++)             /*将H,N元素放入雅克比矩阵*/
         { for(j=0;j<N-1;j++)
           {jac[2*i][2*j]=h[i][j];
            jac[2*i][2*j+1]=nn[i][j];
           }
         }
        for(i=0;i<M;i++)              /*将J,L元素放入雅克比矩阵*/
         { for(j=0;j<N-1;j++)
            { jac[2*i+1][2*j]=jj[i][j];
              jac[2*i+1][2*j+1]=l[i][j];
             }
          }
        for(i=M;i<N-1;i++)              /*将R,S元素放入雅克比矩阵*/
          { for(j=0;j<N-1;j++)
             { jac[2*M+1][2*j]=r[i][j];
               jac[2*M+1][2*j+1]=s[i][j];
             }
          }
     
        for(i=0;i<N-1;i++)            /*得到功率,电压不平衡向量*/
          { dpq[2*i]=dp[i];

           }
        for(i=0;i<M;i++)
          { dpq[2*i+1]=dq[i];

           }
        for(i=0;i<V;i++)
           {  dpq[2*M+1+i]=du2[M];

           }

          for(i=0;i<n;i++)
		  {for(j=0;j<n;j++)
		    aa[i][j]=jac[i][j];
		   
		     aa[i][j]=dpq[i];
		   } 	
		for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)
			a[i][j]=aa[i-1][j-1];
		
		
		
		
		
		det=gauss(a,x);/*用高斯消元法解修正方程*/

        
        for(i=0;i<2*(N-1);i++)      /*得到df与de的解*/
         dfe[i]=x[i+1];

        k=k+1;                      /*迭代次数加1*/
        if(k>=10)                   /*迭代次数过多,出错*/
         { printf("error");
           return(0);
         }

        for(i=0;i<2*(N-1);i++)      /*得到下一轮迭代f,e的初值*/
         { fe[i]+=dfe[i];
          }
        for(i=0;i<(N-1);i++)
          {  f[i]=fe[2*i];
             e[i]=fe[2*i+1];
          }
        for(i=0;i<2*(N-1);i++)
          dfe[i]=fabs(dfe[i]);

        for(i=0;i<2*(N-1);i++)     /*判断收敛条件*/
          { if(dfe[i]<=0.00001)
            re=0;
            else {re=1; break;}
          }

}
        while(re);

        for(i=0;i<N;i++)       /*显示节点电压*/
         {printf("e(%d)=%f\n",i+1,e[i]);
          printf("f(%d)=%f\n",i+1,f[i]);

         }

    


	return 0;
}