yangshuo.c

故障分析 用c语言编程的 利用pq pv节点进行的潮流计算

#include <stdio.h>
 #include <math.h>
 #define  M 50           /*矩阵阶数*/
 #define  N 50        /*迭代次数*/


 int i,j,k,l,z;              /* 循环变量 */
 int n,                      /* 节点数 */
     m,                      /* 支路数 */
	 dd,                     /*对地支路数*/
     pq,                     /* PQ节点数 */
     pv;                   /* PV节点数 */
   int byqn;           /* 变压器数*/ 
 float eps,	                 /* 精度 */
       max;
 char filename[10];

 static float G[M][M],B[M][M],
	      D[M],ykb[M][M];

 struct jiedian                    /* 节点结构体 */
      {  int num,s;                /* num为节点号，s为节点类型*/
	     float p,q,e,f,v;
      } jiedian[M];

 struct zhilu                      /* 支路结构体 */
      {  int num;
	 int p1,p2;                /*支路的两个节点*/
	 float r,x;                /*支路的电阻与电抗*/
      } zhilu[M];
struct byq             /*变压器结构体*/
      { int num;
        int p1,p2;
		float r,x,k;
		float s;
       } byq[M];
  FILE *fp;

 void data1()                 /* 读取数据 */
 { int h,numb;
   printf("               潮  流  上  机  计  算            \n");
   printf("        -----------------------------------------\n");
   printf("                                                 \n");
   printf("                    华北电力大学                 \n");
   printf("                      张炜                       \n");
   printf("                                                 \n");
   printf("              班级：  电力045                    \n");
   printf("              学号：200401010526                 \n");
   printf("              Email:[e-mail]           \n");          /*一些基本的信息*/
   printf("                                                 \n");
   printf("        ---------------------------------------  \n");
 printf("Please putin the in filename!\n");
 scanf("%s",&filename);
   if((fp=fopen(filename,"r"))==NULL)
     { printf(" can not open the input file,please check it!\n");
       exit(0);
     }
   fscanf(fp,"%d,%d,%d,%d,%d,%d,%f\n",&n,&m,&dd,&byqn,&pq,&pv,&eps);
   j=1;k=pq+1;


   for(i=1;i<=n;i++)                                        /*输入节点类型的输入功率和节电电压初值*/
      {
	fscanf(fp,"%d,%d",&numb,&h);
	if(h==1)                                                /*类型h=1是PQ节点*/
	  { fscanf(fp,",%f,%f,%f,%f\n",&jiedian[j].p,&jiedian[j].q,
				       &jiedian[j].e,&jiedian[j].f);
	    jiedian[j].num=numb;
	    jiedian[j].s=h;
	    j++;
	   }
	if(h==2)                                                 /*类型h=2是pv节点*/
	  { fscanf(fp,",%f,%f,%f,%f\n",&jiedian[k].p,&jiedian[k].v,
				       &jiedian[k].e,&jiedian[k].f);
	    jiedian[k].num=numb;
	    jiedian[k].s=h;
	    k++;
	   }
	if(h==3)                                                 /*类型h=3是平衡节点*/
	  { fscanf(fp,",%f,%f\n",&jiedian[n].e,&jiedian[n].f);
	    jiedian[n].num=numb;
	    jiedian[n].s=h;
	   }
       }
   for(i=1;i<=m;i++)                                         /*输入支路阻抗*/
       fscanf(fp,"%d,%d,%d,%f,%f\n",&zhilu[i].num,  &zhilu[i].p1,&zhilu[i].p2,   &zhilu[i].r,  &zhilu[i].x);
   for(i=1;i<=byqn;i++)
	   fscanf(fp,"%d,%d,%d,%f,%f,%f,%f\n",&byq[i].num,&byq[i].p1,&byq[i].p2,&byq[i].r,&byq[i].x,&byq[i].k,&byq[i].s);
    /*变压器编号num，所连节点p1，p2，阻抗r，x，变比k，容量s*/
   fclose(fp);
   printf("Please putin the out filename!\n");
        scanf("%s",&filename);
   if((fp=fopen(filename,"w"))==NULL)
     { printf("  can not write the output file,please check !\n");
       exit(0);
       }
  
  fprintf(fp,"\n            电力系统潮流计算上机  电力045  张炜   200401010526\n");
  fprintf(fp,"\n                   ********** 原始数据 **********\n");
  fprintf(fp,"  ===================================================================\n");
  fprintf(fp,"  节点数:%2d  支路数:%2d  对地支路数:%2d  变压器数:%2d  PQ节点数:%2d  PV节点数:%2d   精度:%f\n",
								 n,m,dd,byqn,pq,pv,eps);
  fprintf(fp,"  -------------------------------------------------------------------\n");
  for(i=1;i<=pq;i++)
  fprintf(fp,"   PQ节点      节点%2d       P[%d]=%f    Q[%d]=%f\n",
	       jiedian[i].num,jiedian[i].num,jiedian[i].p,jiedian[i].num,jiedian[i].q);
  for(i=pq+1;i<=pq+pv;i++)
  fprintf(fp,"   PV节点      节点%2d       P[%d]=%f     V[%d]=%f\n",
	       jiedian[i].num,jiedian[i].num,jiedian[i].p,jiedian[i].num,jiedian[i].v);
  fprintf(fp,"  平衡节点     节点%2d       e[%d]=%f     f[%d]=%f\n",
	       jiedian[n].num,jiedian[n].num,jiedian[n].e,jiedian[n].num,jiedian[n].f);
  fprintf(fp,"  -------------------------------------------------------------------\n");
  for(i=1;i<=m;i++)
  fprintf(fp,"    支路%2d     相关节点:%2d,%2d   R=%f  X=%f\n",
		      i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x);
  for(i=1;i<=byqn;i++)
	  fprintf(fp,"变压器编号%2d    相关节点:%2d,%2d   R=%f  X=%f  变比K=%f\n",byq[i].num,
	  byq[i].p1,byq[i].p2,byq[i].r,byq[i].x,byq[i].k);
  fprintf(fp,"  ===================================================================\n");

 }


 void form_y()                    /* 形成节点导纳矩阵 */
 {
  float S,g,b,g1,b1,g2,b2;int bn;
  for(i=1;i<=m;i++)               /*节点导纳矩阵的主对角线上的导纳*/
  { for(j=1;j<=n;j++)/*m个支路，n个节点*/
	 if((zhilu[i].p1==j)||(zhilu[i].p2==j))/*.p1.p2第i个支路的两个节点号*/
	   { S=zhilu[i].r*zhilu[i].r+zhilu[i].x*zhilu[i].x;/*  [1/(r+jx)]=(r-jx)/(r*r+x*x) */
	     G[j][j]+=zhilu[i].r/S;
	     B[j][j]+=-zhilu[i].x/S;
	   for(bn=1;bn<=byqn;bn++)
	   {if(byq[bn].p1==j)
		{ S=byq[bn].r*byq[bn].r+byq[bn].x*byq[bn].x;
		 g=byq[bn].r/S;
		 b=-byq[bn].x/S;
		 g2=g/byq[bn].k;
		 b2=b/byq[bn].k; 
		 g1=g*((1-byq[bn].k)/(byq[bn].k*byq[bn].k));
		 b1=b*((1-byq[bn].k)/(byq[bn].k*byq[bn].k));
		 G[j][j]=G[j][j]+g2+g1;
		 B[j][j]=B[j][j]+b2+b1;}
		 if(byq[bn].p2==j)
		 {S=byq[bn].r*byq[bn].r+byq[bn].x*byq[bn].x;
		 g=byq[bn].r/S;
		 b=-byq[bn].x/S;
		 g2=g/byq[bn].k;
		 b2=b/byq[bn].k;
		 g1=g*(byq[bn].k-1)/byq[bn].k;
		 b1=b*(byq[bn].k-1)/byq[bn].k;
		 G[j][j]=G[j][j]+g2+g1;
		 B[j][j]=B[j][j]+b2+b1;}
	 }
	 }
  }

  for(j=1;j<=m;j++)              /*对地导纳*/
  if((zhilu[j].p1==0)||(zhilu[j].p2==0))
    {S=zhilu[j].r*zhilu[j].r+zhilu[j].x*zhilu[j].x;
     if(S==0) {G[zhilu[j].p1][zhilu[j].p2]=B[zhilu[j].p1][zhilu[j].p2]=G[zhilu[j].p2][zhilu[j].p1]=B[zhilu[j].p2][zhilu[j].p1]=0;continue;}
        G[zhilu[j].p1][zhilu[j].p2]+=zhilu[i].r/S;
        B[zhilu[j].p1][zhilu[j].p2]+=-zhilu[i].x/S;
        G[zhilu[j].p2][zhilu[j].p1]=G[zhilu[j].p1][zhilu[j].p2];
        B[zhilu[j].p2][zhilu[j].p1]=B[zhilu[j].p1][zhilu[j].p2];
    }

  for(k=1;k<=m;k++)              /*节点导纳矩阵非主对角线上的导纳*/
      { i=zhilu[k].p1;
	j=zhilu[k].p2;
	S=zhilu[k].r*zhilu[k].r+zhilu[k].x*zhilu[k].x;
	G[i][j]+=-zhilu[k].r/S;
	B[i][j]+=zhilu[k].x/S;
	G[j][i]=G[i][j];
	B[j][i]=B[i][j];
	}
   for(bn=1;bn<=byqn;bn++)
      {i=byq[bn].p1;j=byq[bn].p2;
       S=byq[bn].r*byq[bn].r+byq[bn].x*byq[bn].x;
        G[i][j]+=-(byq[bn].r/S)/byq[bn].k;
        B[i][j]+=(byq[bn].x/S)/byq[bn].k;
        G[j][i]=G[i][j];
        B[j][i]=B[i][j];}
		
  /*输出节点导纳矩阵*/
  fprintf(fp,"\n\n             ********* 计算结果 *********\n");
  fprintf(fp,"\n  节点导纳矩阵为:");
  for(i=1;i<=n;i++)
     { fprintf(fp,"\n     ");
       for(j=1;j<=n;j++)
	   fprintf(fp,"%8.5f+j%8.5f   ",G[i][j],B[i][j]);
  }
    fprintf(fp,"\n");
   }


 
 void form_ykb()                      /* 形成雅可比矩阵 */
 { float ei,ej,fi,fj,a=0,b=0;
   int i1,j1,k1;
   float dP[M],dQ[M],dU2[M];
   
   
   for(i=1;i<=2*(pq+pv);i++)           /*初始化矩阵*/
      { dP[i]=0;
	dQ[i]=0;
	dU2[i]=0;
	D[i]=0;
	for(j=1;j<=2*(pq+pv)+1;j++)
	   ykb[i][j]=0;
      }


    for(i=1;i<=pq;i++)                   /*PQ节点函数残量*/
      {  ei=jiedian[i].e;
	 fi=jiedian[i].f;
	 for(j=1;j<=n;j++)
	  { i1=jiedian[i].num;
	    j1=jiedian[j].num;
	    ej=jiedian[j].e;
	    fj=jiedian[j].f;
	    dP[i]+=ei*(G[i1][j1]*ej-B[i1][j1]*fj)
		  +fi*(G[i1][j1]*fj+B[i1][j1]*ej);
	    dQ[i]+=fi*(G[i1][j1]*ej-B[i1][j1]*fj)
		  -ei*(G[i1][j1]*fj+B[i1][j1]*ej);
	    }
	 dP[i]=jiedian[i].p-dP[i];
	 dQ[i]=jiedian[i].q-dQ[i];
	}


   for(i=pq+1;i<=pq+pv;i++)                /*PV节点函数残量*/
	{  ei=jiedian[i].e;
	   fi=jiedian[i].f;
	   for(j=1;j<=n;j++)
	    { i1=jiedian[i].num;
	      j1=jiedian[j].num;
	      ej=jiedian[j].e;
	      fj=jiedian[j].f;
	      dP[i]+=ei*(G[i1][j1]*ej-B[i1][j1]*fj)
		    +fi*(G[i1][j1]*fj+B[i1][j1]*ej);
	     }
	   dP[i]=jiedian[i].p-dP[i];
	   dU2[i]=jiedian[i].v*jiedian[i].v-ei*ei-fi*fi;
	}


   for(i=1;i<=pq;i++)                                   /*形成函数残量矩阵D[M]*/
      {  D[2*i-1]=dP[i];
	 D[2*i]=dQ[i]; }/*pq节点部分，奇数位为dP,偶数位为dQ*/
   for(i=pq+1;i<=pq+pv;i++)
      {  D[2*i-1]=dP[i];
	 D[2*i]=dU2[i];  }/*pv节点部分，奇数位为dP,偶数位为dV2*/

                              /*形成pq节点子阵*/
   for(i=1;i<=pq;i++)                                 
      for(j=1;j<n;j++)
	 {  i1=jiedian[i].num;
	    j1=jiedian[j].num;
	    ei=jiedian[i].e;
	    ej=jiedian[j].e;
	    fi=jiedian[i].f;
	    fj=jiedian[j].f;
	 
		if(i!=j)                                     /*求i!=j时的H、N、J、L*/
	    {  ykb[2*i-1][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi;      /* H */
	       ykb[2*i-1][2*j]=G[i1][j1]*ei+B[i1][j1]*fi;         /* N */
	       ykb[2*i][2*j-1]=-G[i1][j1]*ei-B[i1][j1]*fi;        /* J */
	       ykb[2*i][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi;          /* L */
	     }