flow手算帮手.c

电力系统潮流计算程序 对电力系统自动化的学生是必要的

 /*读取PQ节点的注入功率*/
 for(i=1;i<=npq;i++)
    {fscanf(fp,"%f",PE+i);
     fscanf(fp,"%f",QUE+i);
    }
 
 /*读取PV节点的P && V并得出U^2*/
 for(i=npq+1;i<=n-1;i++)
    {fscanf(fp,"%f",PE+i);
     fscanf(fp,"%f",e+i);
     QUE[i]=e[i]*e[i];
    }
 
 /*读取平衡节点的e && f*/
 fscanf(fp,"%f",e+n);
 fscanf(fp,"%f",f+n);
 

/*读取各支路的起点号b，终点号e，支路阻抗R && X*/
 for(i=1;i<=z;i++)
    {fscanf(fp,"%d",&zl[i].b);
     fscanf(fp,"%d",&zl[i].e);
     fscanf(fp,"%f",&zl[i].RG);
     fscanf(fp,"%f",&zl[i].XB);
    }

 /*读取计算精度E*/
 fscanf(fp,"%f",&E);
 
 /*读取最大迭代次数km*/
 fscanf(fp,"%d",&km);

  /*给PQ节点的e赋初值*/
 for(i=1;i<=npq;i++)
    e[i]=1;
 
 printf("                          欢迎使用潮流精灵2007!\n"); 
 fprintf(fp1,"                          欢迎使用潮流精灵2007!\n"); 


 /*求取节点导纳矩阵并输出*/
 fprintf(fp1,"求节点导纳矩阵：\n"); 
 dn(fp1);
 fprintf(fp1,"\n导纳矩阵中的各元素:Y11=......Y12=......Ynn=......\n\n"); 
 
 fprintf(fp1,"形成导纳矩阵YB:\nYB=\n");
 for(i=1;i<=n;i++)
   {for(j=1;j<=n;j++)
      fprintf(fp1,"   %f+j(%f)",G[i][j],B[i][j]);
    fprintf(fp1,"\n"); 
   }
 fprintf(fp1,"\n\n\n"); 

 
 
 
 /*该循环为迭代循环，直到达到规定精度或最大迭代次数为止*/
 do
{ /*输出分隔线*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  

 fprintf(fp1,"第%d次迭代：K=%d\n\n",K,K);

 fprintf(fp1,"1、计算各PQ、PV节点功率的不平衡量，及PV节点电压的不平衡量:\n\n");
 fprintf(fp1,"取：\n");
 for(i=1;i<=npq;i++)
    fprintf(fp1,"   U%d(%d)=e%d(%d)+f%d(%d)=%f+j%f\n",i,K,i,K,i,K,e[i],f[i]);
 for(i=npq+1;i<=n-1;i++)
    fprintf(fp1,"   U%d(%d)=e%d(%d)+f%d(%d)=%f+j%f\n",i,K,i,K,i,K,e[i],f[i]);
 fprintf(fp1,"   节点%d是平衡节点,保持U%d=e%d+f%d=%f+j%f为定值。\n\n",n,n,n,n,e[n],f[n]);
 
 
 
 /*计算PQ节点的P和Q*/
 fprintf(fp1,"a/计算各PQ、PV节点功率：\n");
 
 fprintf(fp1,"Pi(%d)=求和... Qi(%d)=求和... Ui(%d)=ei(%d)^2+fi(%d)^2\n\n",K,K,K,K,K);
 fprintf(fp1,"下面的符号别忘了带角标!\n\n");
 
 PQ(fp1);
 /*计算PV节点的P和V*/
 PV(fp1);


 /*计算PQ节点的DP,DQ*/
 fprintf(fp1,"于是：DPi=Pi-Pi(%d);DQi=Qi-Qi(%d);DUi^2=Ui^2-(ei(%d)^2+fi(%d)^2)\n\n",K,K,K,K);
 for(i=1;i<=npq;i++)
    {DP[i]=PE[i]-P[i];
     DQU[i]=QUE[i]-QU[i];
     fprintf(fp1,"DP%d=P%d-P%d(%d)=%f-(%f)=%f\n",i,i,i,K,PE[i],P[i],DP[i]);
     fprintf(fp1,"DQ%d=Q%d-Q%d(%d)=%f-(%f)=%f\n",i,i,i,K,QUE[i],QU[i],DQU[i]);
    }
 /*计算PV节点的DP,D(U^2)*/
 for(i=npq+1;i<=n-1;i++)
    {DP[i]=PE[i]-P[i];
     DQU[i]=QUE[i]-QU[i];
     fprintf(fp1,"DP%d=P%d-P%d(%d)=%f-(%f)=%f\n",i,i,i,K,PE[i],P[i],DP[i]);
     fprintf(fp1,"D(U%d)^2=(U%d)^2-(e%d(%d)^2+f%d(%d)^2)=%f-((%f)^2+(%f)^2)=%f\n",i,i,i,K,i,K,QUE[i],e[i],f[i]);
    }
 fprintf(fp1,"\n\n");

 fprintf(fp1,"2、计算雅可比矩阵中各元素:\n\n");


 /*计算PQ节点的aii和bii*/
 for(i=1;i<=npq;i++)
   {ai[i]=(P[i]*e[i]+QU[i]*f[i])/(e[i]*e[i]+f[i]*f[i]);
    bi[i]=(P[i]*f[i]-QU[i]*e[i])/(e[i]*e[i]+f[i]*f[i]);  
    fprintf(fp1,"I%d=(P%d-jQ%d)/U%d*=(%f-j%f)/(%f-j%f)=%f+j%f=a%d%d+jb%d%d\n",i,i,i,i,P[i],QU[i],e[i],f[i],ai[i],bi[i],i,i,i,i); 
   }
 fprintf(fp1,"\n\n");

/*计算PV节点的aii和bii*/
 for(i=npq+1;i<=n-1;i++)
   {fprintf(fp1,"Q%d=",i);
    QU[i]=0;
    for(j=1;j<=n;j++)
       {QU[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
        fprintf(fp1,"+%f*(%f*%f-%f*%f)-%f*(%f*%f+%f*%f)\n   ",f[i],G[i][j],e[j],B[i][j],f[j],e[i],G[i][j],f[j],B[i][j],e[j]);
       }
     fprintf(fp1,"=%f\n",QU[i]);
    ai[i]=(P[i]*e[i]+QU[i]*f[i])/(e[i]*e[i]+f[i]*f[i]);
    bi[i]=(P[i]*f[i]-QU[i]*e[i])/(e[i]*e[i]+f[i]*f[i]);  
    fprintf(fp1,"I%d=(P%d-jQ%d)/U%d*=(%f-j%f)/(%f-j%f)=%f+j%f=a%d%d+jb%d%d\n",i,i,i,i,P[i],QU[i],e[i],f[i],ai[i],bi[i],i,i,i,i); 
   }
 fprintf(fp1,"\n\n");


fprintf(fp1,"雅可比矩阵的各个元素分别为：\n\n");
 
/*求PQ节点形成的雅可比矩阵*/
 PQykb(fp1);


/*求PV节点形成的雅可比矩阵*/
 PVykb(fp1);

fprintf(fp1,"3、所以可以得到K=%d时的雅可比矩阵:\n",K);
 /*输出雅克比矩阵*/
fprintf(fp1,"J(%d)=\n",K);
for(i=1;i<=2*(n-1);i++)
   {for(j=1;j<=2*(n-1);j++)
      fprintf(fp1,"  %f",J[i][j]);
   fprintf(fp1,"\n");
   }
fprintf(fp1,"\n\n");

fprintf(fp1,"至此，可以建立修正方程组如下：....\n\n");


/*用列主元高斯消去法来接上面行成的线性方程组*/
for(i=1;i<=n-1;i++)
   {b1[2*i-1]=DP[i];
    b1[2*i]=DQU[i];
   }
lzy(2*(n-1),J,b1,fp1);
 
/*把解方程的结果赋给df与de并求最大不平衡量的绝对值dm*/
dm=0;
for(i=1;i<=n-1;i++)
   {df[i]=b1[2*i-1];
    de[i]=b1[2*i];
    if(fabs(df[i])>dm)
       dm=fabs(df[i]);
    if(fabs(de[i])>dm)
       dm=fabs(de[i]);
   }

fprintf(fp1,"解得：\n");
for(i=1;i<=n-1;i++)
   {fprintf(fp1,"   df%d=%f\n",i,df[i]);
    fprintf(fp1,"   de%d=%f\n",i,de[i]);
   }
fprintf(fp1,"\n\n");

fprintf(fp1,"因为：ei(%d)=ei(%d)+dei(%d); fi(%d)=fi(%d)+dfi(%d)\n所以：\n",K+1,K,K,K+1,K,K);
/*根据不平衡量求f与e的修正值并输出*/
for(i=1;i<=n-1;i++)
   {f[i]+=df[i];
    e[i]+=de[i];
   }

for(i=1;i<=n-1;i++)
   {fprintf(fp1,"   e%d(%d)=e%d(%d)+de%d(%d)=%f+%f=%f\n",i,K+1,i,K,i,K,e[i]-de[i],de[i],e[i]);
    fprintf(fp1,"   f%d(%d)=f%d(%d)+df%d(%d)=%f+%f=%f\n",i,K+1,i,K,i,K,f[i]-df[i],df[i],f[i]);
   }
fprintf(fp1,"\n\n");

K++;

}
while(K<km && dm>E);

 
 /*输出分隔线*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 
 /*输出各节点的电压*/
 fprintf(fp1,"各节点电压为:\n");
 for(i=1;i<=n;i++)
    fprintf(fp1,"   U%d=e%d+jf%d=%f+j(%f)\n",i,i,i,e[i],f[i]);


/*输出分隔线*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 /*计算平衡功率S~n*/
 fprintf(fp1,"平衡节点功率为:\n");
 i=n;
 PE[i]=0;
 for(j=1;j<=n;j++)
    PE[i]+=e[i]*(G[i][j]*e[j]-B[i][j]*f[j])+f[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
 QUE[i]=0;
 for(j=1;j<=n;j++)
    QUE[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
 fprintf(fp1,"S~%d=......\n",i);
 fprintf(fp1,"   =(%f+j%f)[",e[i],f[i]);
 for(j=1;j<=n;j++)
    fprintf(fp1,"(%f-j(%f))(%f-j(%f))",G[i][j],B[i][j],e[j],f[j]);
    
 fprintf(fp1,"]\n   =%f+j(%f)\n",PE[i],QUE[i]);




/*输出分隔线*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 
 /*计算线路功率*/
 xlgl(fp1);

 
 /*求网络总损耗*/
 fprintf(fp1,"\n法二:\n");
 for(i=npq+1;i<=n-1;i++)
    {fprintf(fp1,"Q%d=",i);
	 QUE[i]=0;
     for(j=1;j<=n;j++)
     {QUE[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
      fprintf(fp1,"+%f*(%f*%f-%f*%f)-%f*(%f*%f+%f*%f)\n   ",f[i],G[i][j],e[j],B[i][j],f[j],e[i],G[i][j],f[j],B[i][j],e[j]);
     }
	 fprintf(fp1,"=%f\n\n\n",QUE[i]);
    }
 SP=0;SQ=0;
 fprintf(fp1,"DDS~=");
 for(i=1;i<=n;i++)
    {SP+=PE[i];
     SQ+=QUE[i];
	 fprintf(fp1,"+(%f+j%f)",PE[i],QUE[i]);
    }
 fprintf(fp1,"\n    =%f+%f",SP,SQ);

 fprintf(fp1,"\n\n******************************************THE END**********************************************");
   
 fclose(fp);
 fclose(fp1);
 
 printf("\n\n\n\n           潮流计算已完成！请打开输出文件out.txt查看计算结果！\n\n\n\n\n\n                       ");
 
}