题号：2编程计算潮流（电力系统分析）

最新问答

• Brita阿菜

  #include "stdio.h"
  #include "math.h"
  
  #define M 20 /*节点数、支路数极限值*/
  #define N 5 /*迭代次数极限值*/
  
  int n,m,dd=0,pq=0,pv=0,tt,qy;
  float eps; /*节点数、支路数、对地支路数、pq节点数、pv节点数、控制打印、互联网数、精度*/
  
  struct jiedian /*节点结构体*/
  {int s; /*节点类型：1-PQ节点；2-PV节点；3-平衡节点*/
  float p,q,e,f,v; /*节点的有功、无功、电压实部、虚部、电压辐值*/
  }jiedian[M]; /*如引用节点1的无功，则为jiedian[1].q */
  
  struct zhilu /*支路结构体*/
  {int p1,p2,s; /*支路两端节点号，支路类型：1-普通支路；2-变压器支路；3-对地支路*/
  float r,x,b,kt; /*支路的电阻、电抗、导纳、变压器非标准变比*/
  }zhilu[M];
  
  struct hulianwang /*互联网结构体*/
  {int num,pv; /*互联区域号,指定的pv节点号*/
  float p,eps2; /*规定的有功功率及其允许误差*/
  int count; /*每个互联区域包括的节点数*/
  int a[M]; /*包括的节点*/
  } hulianwang[M];
  
  static float G[M][M],B[M][M],G1[M][M],B1[M][M]; /*节点导纳阵矩阵*/
  static float ykb[2*M][2*M]; /*节点导纳阵矩阵*/
  static float yinzi[2*M][2*M]; /*因子表*/
  static float P[M][M],Q[M][M]; /*潮流计算结果*/
  
  FILE *fp1,*fp2; /*文件输入、输出指针*/
  
  void input() /*从文件in.txt 输入线路基本参数、节点、支路*/
  {
  int i,j,h; /*循环变量、节点（支路）类型*/
  
  /**********************************打开电网参数文件**********************/
  fp1=fopen("in.txt","r");
  if (fp1==NULL)
  {
  printf("Can not open file in.txt !\n");
  exit(0);
  }
  
  /***********读取节点数,支路数,互联网数,PQ节点数PV节点数和精度***********/
  fscanf(fp1,"%d,%d,%d,%f\n",&n,&m,&qy,&eps);
  
  /************************读取节点信息*************************************/
  for(i=1;i<=n;i++)
  {
  fscanf(fp1,"%d",&h);
  if(h==1) /*类型h=1是PQ节点*/
  {
  pq++;
  fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].q,&jiedian[i].e,&jiedian[i].f);
  jiedian[i].s=1;
  jiedian[i].v=sqrt(jiedian[i].e*jiedian[i].e+jiedian[i].f*jiedian[i].f);
  }
  else if(h==2) /*类型h=2是pv节点*/
  {
  pv++;
  fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].v,&jiedian[i].e,&jiedian[i].f);
  jiedian[i].s=2;
  jiedian[i].q=0;
  }
  else /*类型h=3是平衡节点*/
  {
  jiedian[i].p=0;
  jiedian[i].q=0;
  jiedian[i].e=1;
  jiedian[i].f=0;
  jiedian[i].v=1;
  jiedian[i].s=3;
  }
  }
  
  /*******************************读取支路信息*************************************/
  for(i=1;i<=m;i++)
  {
  fscanf(fp1,"%d",&h);
  if(h==1) /*类型h=1是普通支路*/
  {
  fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].b);
  zhilu[i].kt=1;
  zhilu[i].s=1;
  }
  if(h==2) /*类型h=2是变压器支路*/
  {
  fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].kt);
  zhilu[i].s=2;
  }
  if(h==3) /*类型h=3是接地支路*/
  {
  fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,&zhilu[i].b);
  zhilu[i].kt=1;
  zhilu[i].s=3;
  dd++;
  }
  }
  
  /*******************************读互联网信息*************************************/
  if(qy!=0)
  for(i=1;i<=qy;i++) /*输入互联网状况*/
  {
  fscanf(fp1,"%d,%d,%f,%f,%d",&hulianwang[i].num,&hulianwang[i].pv,&hulianwang[i].p,&hulianwang[i].eps2,&hulianwang[i].count);
  for(j=1;j<=hulianwang[i].count;j++)
  fscanf(fp1,",%d",&(hulianwang[i].a[j]));
  }
  
  fclose(fp1);
  
  /**********************************打开输出结果文件********************************/
  fp2=fopen("out.txt","w");
  if(fp2==NULL)
  {
  printf("Can not open file!\n");
  exit(0);
  }
  
  fprintf(fp2,"\n**************************** 原始如下 ***********************************\n\n ");
  fprintf(fp2," 节点数:%2d 支路数:%2d 对地支路数:%2d PQ节点数:%2d PV节点数:%2d 精度1:%f \n",n,m,dd,pq,pv,eps);
  
  fprintf(fp2,"\n------------------------------------------------------------------------------\n\n ");
  for(i=1;i<=pq;i++)
  fprintf(fp2," 节点%2d PQ节点 P[%d]=%f Q[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].q);
  
  for(i=pq+1;i<=pq+pv;i++)
  fprintf(fp2," 节点%2d PV节点 P[%d]=%f V[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].v);
  fprintf(fp2," 节点%2d 平衡节点\n",i);
  fprintf(fp2,"\n-------------------------------------------------------------------------------\n\n ");
  for(i=1;i<=m;i++)
  {
  if(zhilu[i].s==1)
  fprintf(fp2," 支路%2d：普通支路 相关节点:%2d,%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].b);
  else if(zhilu[i].s==2)
  fprintf(fp2," 支路%2d：变压器支路 相关节点:%2d,%2d R=%f X=%f Kt=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].kt);
  else
  fprintf(fp2," 支路%2d： 对地支路 相关节点:%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].r,zhilu[i].x,zhilu[i].b);
  
  }
  for(i=1;i<=qy;i++)
  {
  fprintf(fp2," 互联网区域:%2d 指定pv节点:%2d 规定有功功率:%8.5f 允许误差:%8.5f 包括的节点: ",hulianwang[i].num,hulianwang[i].pv,hulianwang[i].p,hulianwang[i].eps2);
  for (j=1;j<=hulianwang[i].count;j++)
  fprintf(fp2,"%2d ",hulianwang[i].a[j]);fprintf(fp2,"\n");
  }
  
  }
  
  void youhua() /*利用节点、支路形成新编号*/
  {
  int i,k,j,jd[N]; /*jd[N]记录节点*/
  struct jiedian tem; /*中间临时节点所连支路数*/
  
  for(i=1;i<=n;i++) /*对节点连接支路数赋0*/
  {
  jd[i]=0;
  }
  
  for(i=1;i<=n;i++) /*考虑每个节点*/
  {
  for(j=1;j<=m;j++) /*考虑每条支路*/
  {
  if((zhilu[j].p1==i||zhilu[j].p2==i)&&(zhilu[j].s!=3)) /*如果是这条支路的节点且不是接地支路*/
  jd[i]++;
  
  }
  }
  
  for(i=1;i {
  for(j=i+1;j<=pq;j++) /*从小到大排序算法*/
  {
  if(jd[i]>jd[j])
  {
  tem=jiedian[i];
  jiedian[i]=jiedian[j];
  jiedian[j]=tem;
  
  for(k=1;k<=m;k++) /*更新支路信息*/
  {
  if(zhilu[k].p1==j)
  zhilu[k].p1=i;
  else if(zhilu[k].p2==j)
  zhilu[k].p2=i;
  else if(zhilu[k].p1==i)
  zhilu[k].p1=j;
  else if(zhilu[k].p2==i)
  zhilu[k].p2=j;
  }
  
  }
  }
  }
  
  for(i=pq+1;i {
  for(j=i+1;j<=pq+pv;j++) /*从小到大排序算法*/
  {
  if(jd[i]>jd[j])
  {
  tem=jiedian[i];
  jiedian[i]=jiedian[j];
  jiedian[j]=tem;
  
  for(k=1;k<=m;k++) /*更新支路信息*/
  {
  if(zhilu[k].p1==j)
  zhilu[k].p1=i;
  else if(zhilu[k].p2==j)
  zhilu[k].p2=i;
  else if(zhilu[k].p1==i)
  zhilu[k].p1=j;
  else if(zhilu[k].p2==i)
  zhilu[k].p2=j;
  }
  
  }
  }
  }
  
  fprintf(fp2,"\n**************************** 节点优化结果如下 *********************************\n\n ");
  
  for(i=1;i<=n;i++)
  {
  if(jiedian[i].s==1)
  fprintf(fp2," 节点%2d PQ节点 P[%d]=%f Q[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].q);
  if(jiedian[i].s==2)
  fprintf(fp2," 节点%2d PV节点 P[%d]=%f V[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].v);
  if(jiedian[i].s==3)
  fprintf(fp2," 节点%2d 平衡节点 \n",i);
  }
  }
  
  void form_y() /*利用支路形成Y，注意对地支路*/
  {
  int i,j,k;
  float S;
  
  for(i=1;i<=n;i++)
  for(j=1;j<=n;j++)
  G[i][j]=B[i][j]=0;
  
  for(i=1;i<=m;i++) /*节点导纳矩阵的主对角线上的导纳*/
  for(j=1;j<=n;j++)
  if(zhilu[i].p1==j||zhilu[i].p2==j)
  {
  S=zhilu[i].r*zhilu[i].r+zhilu[i].x*zhilu[i].x;
  if(S==0) continue;
  G[j][j]+=zhilu[i].r/S;
  B[j][j]+=-zhilu[i].x/S;
  
  if(zhilu[i].s==1) /*如果是普通支路*/
  {
  B[j][j]+=zhilu[i].b/2;
  }
  else if(zhilu[i].s==2) /*如果是普通变压器支路*/
  {
  if(zhilu[i].p1==j)
  {
  G[j][j]+=(zhilu[i].r/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);
  B[j][j]+=(-zhilu[i].x/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);
  }
  else if(zhilu[i].p2==j)
  {
  G[j][j]+=(zhilu[i].r/S*(zhilu[i].kt-1))/zhilu[i].kt;
  B[j][j]+=(-zhilu[i].x/S*(zhilu[i].kt-1))/zhilu[i].kt;
  }
  }
  else if(zhilu[i].s==3) /*如果是对地支路*/
  {
  B[j][j]+=zhilu[i].b;
  }
  }
  
  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;
  if(S==0) continue;
  G[i][j]+=-zhilu[k].r/S;
  B[i][j]+=zhilu[k].x/S;
  if(zhilu[k].kt!=1.0)
  {
  G[i][j]/=zhilu[k].kt;
  B[i][j]/=zhilu[k].kt;
  }
  G[j][i]=G[i][j];
  B[j][i]=B[i][j];
  }
  
  fprintf(fp2,"\n\n****************************** 节点导纳矩阵为 **********************************\n");
  for(i=1;i<=n;i++)
  {
  fprintf(fp2,"\n ");
  for(j=1;j<=n;j++)
  fprintf(fp2,"%8.5f+j%8.5f ",G[i][j],B[i][j]);
  }
  }
  
  void form_j() /*利用节点和Y形成J*/
  
  {
  float ei,fi,a=0,b=0;
  int i1,j1,k1,i,j,k;
  
  for(i=1;i<=2*(pq+pv)+1;i++)
  for(j=1;j<=2*(pq+pv)+1;j++)
  ykb[i][j]=0;
  
  for(i=1;i<=pq;i++)
  for(j=1;j { i1=i;
  j1=j;
  ei=jiedian[i].e;
  fi=jiedian[i].f;
  
  if(i!=j) /*求i!=j时的H、N、J、L*/
  { ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /* H */
  ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /* N */
  ykb[2*i][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /* J */
  ykb[2*i][2*j]=-G[i1][j1]*ei-B[i1][j1]*fi; /* L */
  }
  
  else /*求i=j时的H、N、J、K*/
  { a=0;b=0;
  for(k=1;k<=n;k++)
  if(k!=i)
  { k1=k;
  a=a+G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;
  b=b+G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;
  }
  ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/
  ykb[2*i][2*j]=-2*B[i1][i1]*jiedian[i].f+a; /*L*/
  ykb[2*i-1][2*j]= 2*G[i1][i1]*jiedian[i].f+b; /*N*/
  ykb[2*i][2*j-1]=-2*B[i1][i1]*jiedian[i].e-b; /*J*/
  }
  }
  
  for(i=pq+1;i<=pq+pv;i++) /*形成pv节点子阵*/
  for(j=1;j { i1=i;
  j1=j;
  ei=jiedian[i].e;
  fi=jiedian[i].f;
  
  if(i!=j) /*求i!=j时的H、N*/
  {
  ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /*H*/
  ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /*N*/
  }
  
  else /*求i=j时的H、N、R、S*/
  { a=0;b=0;
  for(k=1;k<=n;k++)
  if(k!=i)
  { k1=k;
  a+=G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;
  b+=G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;
  }
  ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/
  ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].f+b; /*N*/
  ykb[2*i][2*j-1]=2*ei; /*R*/
  ykb[2*i][2*j]=2*fi; /*S*/
  }
  }
  
  fprintf(fp2,"\n\n\n******************************* 雅可比矩阵*************************************** \n");
  
  for(i=1;i<=2*(pq+pv);i++)
  {
  fprintf(fp2,"\n ");
  for(j=1;j<=2*(pq+pv);j++)
  fprintf(fp2,"%8.5f ",ykb[i][j]);
  }
  }
  
  void form_yz() /*利用J形成因子表*/
  {
  int i,j,k;
  float a[2*M][2*M],L[2*M][2*M],R[2*M][2*M];
  float x,y,z;
  
  for(i=1;i<=2*(pq+pv);i++)
  for(j=1;j<=2*(pq+pv);j++)
  { a[i][j]=0;R[i][j]=0; L[i][j]=0;}
  for(i=1;i<=2*(pq+pv);i++)
  for(j=1;j<=2*(pq+pv);j++)
  a[i][j]=ykb[i][j];
  
  for(i=1;i<=2*(pq+pv);i++)
  {
  for(j=1;j<=2*(pq+pv);j++)
  {
  if(i>j)
  {
  R[i][j]=0;
  x=0;
  for(k=1;k<=j-1;k++)
  {
  x=x+L[i][k]*R[k][j];
  }
  L[i][j]=a[i][j]-x;
  }
  if(i==j)
  {
  R[i][j]=1;
  y=0;
  for(k=1;k<=i-1;k++)
  {
  y=y+L[i][k]*R[k][i];
  }
  L[i][i]=a[i][j]-y;
  }
  if(i {
  L[i][j]=0;
  z=0;
  for(k=1;k<=i-1;k++)
  {
  z=z+L[i][k]*R[k][j];
  }
  R[i][j]=(a[i][j]-z)/L[i][i];
  }
  }
  }
  
  for(i=1;i<=2*(pq+pv);i++)
  {
  for(j=1;j<=2*(pq+pv);j++)
  {
  if(i yinzi[i][j]=R[i][j];
  if(i==j)
  yinzi[i][j]=1/L[i][j];
  if(i>j)
  yinzi[i][j]=L[i][j];
  }
  }
  
  fprintf(fp2,"\n\n\n******************************** 因子表 ***************************************** \n");
  
  for(i=1;i<=2*(pq+pv);i++)
  {
  fprintf(fp2,"\n ");
  for(j=1;j<=2*(pq+pv);j++)
  fprintf(fp2,"%8.5f ",yinzi[i][j]);
  }
  
  }

  浏览 325赞 140时间 2023-09-26
• 黄二小要奋斗

  #include "stdio.h"
  #include "math.h"

  浏览 439赞 115时间 2023-08-31
• 爱吃爱疯

  最讨厌精度计算了 T_T 稍微乱一点就错

  浏览 188赞 51时间 2023-07-07