开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 行业应用 > 能源行业(电力石油煤炭) > 查看源码
pq_decouple.C
资源名称：PQ_Power_Flow_Pro.rar　[点击查看]
上传用户：wxl722722
上传日期：2022-07-11
资源大小：10k
文件大小：21k
源码类别：

能源行业(电力石油煤炭)

开发平台：

C/C++

pq_decouple.C：源码内容
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  6. void Datain(int Nb,int Nl,int Ng,FILE *fp,struct Branch_Type *Branch,struct Load_Type *Load,struct Generator_Type *Generator,struct PVNode_Type *PVNode,int *Npv);  /*从初始数据文件中读入数据*/
  7. void AdmittanceMatrix(int N,int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,int *NYseq1);  /*形成导纳矩阵函数*/
  8. void AdmittanceMatrixAdd(int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch);  /*形成导纳矩阵追加接地支路函数*/
  9. void Factorial(int flag,int N,int Npv,struct PVNode_Type *PVNode,int *NUsum,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq,double *D,struct U_Type *U);  /*形成因子表函数*/
  10. void NodePower(int flag,int N,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq);  /*节点功率计算函数*/
  11. void Iteration(int flag,struct Generator_Type *Generator,struct Load_Type *Load,struct PVNode_Type *PVNode,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct GeneratorPower *GenPower,int N,double *DI,double *MaxError,int *ErrNode);  /*迭代计算函数*/
  12. void FormulaSolution(int flag,struct U_Type *U,double *D,int *NUsum,double *DI,int N,struct NodalVol *NodeVol,double V0);  /*线性方程组求解函数*/
  13. void NodeDataOutput(FILE *fp,struct NodalVol *NodeVol,struct Generator_Type *Generator,int N,struct GeneratorPower *GenPower,struct NodalPow *NodePow,struct Load_Type *Load,int Nl);  /*节点数据输出函数*/
  14. void BranchDataOutput(FILE *fp,int Nb,struct Branch_Type *Branch,struct NodalVol *NodeVol);  /*支路数据输出函数*/
  16.   struct Branch_Type
  17.   {
  18.    int i,j;
  19.    double R,X,YK;
  20.   }; 
  21.   
  22.   struct Generator_Type
  23.   {
  24.    int i;
  25.    double P,Q;
  26.    double V;
  27.   };
  28.   
  29.   struct Load_Type
  30.   {
  31.    int i;
  32.    double P,Q;
  33.    double V;
  34.   };
  35.   
  36.   struct PVNode_Type
  37.   {
  38.    int i;
  39.    double V;
  40.   };
  41.   
  42.   struct Yii_Type
  43.   {
  44.    double G,B;
  45.   };
  46.   
  47.   struct Yij_Type
  48.   {
  49.    double G,B;
  50.    int j;
  51.   };
  52.   
  53.   struct NodalPow
  54.   {
  55.    double P,Q;
  56.   };
  57.   
  58.   struct NodalVol
  59.   {
  60.     double V,theta;
  61.   };
  62.   
  63.   struct GeneratorPower
  64.   {
  65.      double P,Q;
  66.   };
  67.   
  68.   struct U_Type
  69.   {
  70.      double value;
  71.      int j;
  72.   };
  73.   
  75. main()
  76. {
  77.   int N;                             /*系统节点总数*/
  78.   int Nb;                            /*系统中支路数*/
  79.   int Ng;                            /*发电机节点总数,不包括平衡机*/
  80.   int Nl;                            /*负荷节点总数*/
  81.   double V0;                         /*系统标称电压*/
  82.   double epsilon;                    /*迭代收敛所要求的精确度*/
  83.   struct Branch_Type *Branch;        /*支路数据数组*/
  84.   struct Generator_Type *Generator;  /*发电机数据数组*/
  85.   struct Load_Type *Load;            /*负荷数据数组*/
  86.   struct PVNode_Type *PVNode;        /*PV节点数据数组*/
  87.   int Npv=0;                         /*PV节点个数*/                             
  88.   struct Yii_Type *Yii,*Yiil;        /*导纳对角矩阵数组*/
  89.   struct Yij_Type *Yij,*Yijl;        /*导纳非对角矩阵数组*/
  90.   int *NYseq;                        /*矩阵非对角元素记录数组*/
  91.   struct NodalPow *NodePow;          /*节点功率数组*/
  92.   struct NodalVol *NodeVol;          /*节点电压数组*/
  93.   struct GeneratorPower *GenPower;   /*发电机功率数组*/
  94.   double *DI1,*DI2;                  /*存放误差项的数组*/
  95.   double MaxError=0.0;               /*每次迭代过程中最大误差项*/
  96.   int ErrNode;                       /*最大误差项所在节点*/
  97.   int Kp=1,Kq=1;                     /*有功迭代和无功迭代的标记*/
  98.   int n,k,count;                     /*K为迭代次数*/
  99.   int i;                          
  100.   struct U_Type *U1,*U2;             /*存放因子表非对角元素的数组*/
  101.   double *D1,*D2;                    /*存放因子表对角元素的数组*/
  102.   int *NUsum1,*NUsum2;               /*矩阵非对角元素记录数组*/ 
  104.   char FILENAME[20];
  105.   FILE *fp;                          
  107.   /****************打开初始数据原文件*****************/
  108.   printf("Please Enter The Filename of The System:");
  109.   gets(FILENAME);
  110.   if((fp=fopen(FILENAME,"r"))==NULL)
  111.   {
  112.     printf("Cannot Find The File:%sn",FILENAME);
  113.     printf("Press ENTER to Escape!");
  114.     exit(0);
  115.   }
  117.   /********************读系统基本数据信息*********************/
  118.   fscanf(fp,"%d,%d,%d,%d,%lf,%lf",&N,&Nb,&Ng,&Nl,&V0,&epsilon);
  120.   /*************给支路、发电机、负荷、PV节点分配内存**********/
  121.   Branch=malloc((Nb+1)*sizeof(struct Branch_Type));
  122.   Generator=malloc((Ng+1)*sizeof(struct Generator_Type));
  123.   Load=malloc((Nl+1)*sizeof(struct Load_Type));
  124.   PVNode=malloc(N*sizeof(struct PVNode_Type));
  126.   Datain(Nb,Nl,Ng,fp,Branch,Load,Generator,PVNode,&Npv);  /*从初始文件中读入数据*/
  128.   for(n=0;1;n++)  /*建立计算结果输出文件*/
  129.   {
  130.     if(FILENAME[n]=='.')
  131.     {
  132.         FILENAME[n]='';
  133.         strcat(FILENAME,"out.dat");
  134.         break;
  135.     }
  136.   }
  137.   if((fp=fopen(FILENAME,"w"))==NULL)  /*打开计算结果输出文件*/
  138.   {
  139.     printf("Cannot Find The File:%sn",FILENAME);
  140.     printf("Press ENTER to Escape!");
  141.     exit(0);
  142.   }
  143.   
  144.   /***********为导纳矩阵分配内存并形成导纳矩阵***************/
  145.   Yii=malloc((N+1)*sizeof(struct Yii_Type));
  146.   Yiil=malloc((N+1)*sizeof(struct Yii_Type));
  147.   Yij=malloc((N+1)*sizeof(struct Yij_Type));
  148.   Yijl=malloc((N+1)*sizeof(struct Yij_Type));
  149.   NYseq=malloc((N+1)*sizeof(int));
  150.   
  151.   AdmittanceMatrix(N,Nb,Yii,Yiil,Yij,Yijl,Branch,NYseq);  /*导纳矩阵形成函数*/
  153.   /*****形成因子表,采用BX方案(B'中忽略充电电容和非标准变比，B"中忽略电阻)*****/
  154.   U1=malloc((N-1)*(N-2)/2*sizeof(struct U_Type));
  155.   U2=malloc((N-1)*(N-2)/2*sizeof(struct U_Type));
  156.   D1=malloc(N*sizeof(double));
  157.   D2=malloc(N*sizeof(double));
  158.   NUsum1=malloc(N*sizeof(int));
  159.   NUsum2=malloc(N*sizeof(int));
  161.   Factorial(1,N,Npv,PVNode,NUsum1,Yii,Yij,NYseq,D1,U1);  /*形成因子表B'*/
  163.   AdmittanceMatrixAdd(Nb,Yii,Yiil,Yij,Yijl,Branch);  /*导纳矩阵追加接地支路函数*/
  165.   Factorial(2,N,Npv,PVNode,NUsum2,Yiil,Yijl,NYseq,D2,U2);  /*形成因子表B"*/
  167.   /****************下面利用所求得的因子表进行迭代求解******************/
  168.   DI1=malloc(N*sizeof(double));
  169.   DI2=malloc(N*sizeof(double));
  170.   NodePow=malloc((N+1)*sizeof(struct NodalPow));
  171.   NodeVol=malloc((N+1)*sizeof(struct NodalVol));
  173.   /********************先送电压初值***********************/
  174.   for(i=1;i<=N;i++)
  175.   {
  176.     NodeVol[i].V=V0;
  177.     NodeVol[i].theta=0.0;
  178.   }
  179.   for(n=1;n<=Npv;n++)
  180.   {
  181.     i=PVNode[n].i;
  182.     NodeVol[i].V=PVNode[n].V;
  183.   }
  184.     
  185.   GenPower=malloc((N+1)*sizeof(struct GeneratorPower));  
  186.    
  187.   fprintf(fp,"ttt系统潮流计算结果 n(1)迭代过程纪录：n迭代次数t    t有功迭代tt无功迭代ntt   ΔPmaxtP-Nodet ΔQmaxttQ-Noden");
  188.   
  189.   for(k=0;1;k++)
  190.   {
  191.     fprintf(fp,"  %2d:t",k);
  192.     if(Kp==1)
  193.     {
  194.         NodePower(1,N,NodeVol,NodePow,Yii,Yij,NYseq);  /*节点功率计算函数*/
  195.         Iteration(1,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI1,&MaxError,&ErrNode);  /*迭代计算函数*/
  196.         fprintf(fp,"t%10.7lft %dt",MaxError,ErrNode);
  197.         if(MaxError>=epsilon)
  198.             FormulaSolution(1,U1,D1,NUsum1,DI1,N,NodeVol,V0);  /*线性方程组求解函数*/
  199.         else
  200.             Kp=0;
  201.     }
  202.     else
  203.         fprintf(fp,"tttt");
  204.     if(Kq==1)
  205.     {
  206.         NodePower(2,N,NodeVol,NodePow,Yii,Yij,NYseq);  /*节点功率计算函数*/
  207.         Iteration(2,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI2,&MaxError,&ErrNode);  /*迭代计算函数*/
  208.         fprintf(fp,"%10.7lft %dn",MaxError,ErrNode);
  209.         if(MaxError>=epsilon)
  210.             FormulaSolution(2,U2,D2,NUsum2,DI2,N,NodeVol,V0);  /*线性方程组求解函数*/
  211.         else
  212.             Kq=0;
  213.     }
  214.     else
  215.         fprintf(fp,"n");
  216.     if(Kp==0&&Kq==0)
  217.         break;
  218.     if(k>1000)
  219.     {
  220.         fprintf(fp,"n迭代次数超过1000次，系统不收敛!n");
  221.     }
  222.     printf("...");
  223.   }
  225.   fprintf(fp,"n总迭代次数为: %d 次!nn",k+1);
  226.   fprintf(fp,"(2)潮流计算结果(节点电压)：n NodettVttθttPttQn");
  227.   NodeDataOutput(fp,NodeVol,Generator,N,GenPower,NodePow,Load,Nl);
  228.   fprintf(fp,"n(3)潮流计算结果(支路功率):n    BranchttPttQn");
  229.   BranchDataOutput(fp,Nb,Branch,NodeVol);
  230.   fprintf(fp,"nnn");
  232.   fclose(fp);
  234.   free(Branch);
  235.   free(Generator);
  236.   free(Load);
  237.   free(PVNode);
  238.   free(Yii);
  239.   free(Yiil);
  240.   free(Yij);
  241.   free(Yijl);
  242.   free(NYseq);
  243.   free(U1);
  244.   free(U2);
  245.   free(D1);
  246.   free(D2);
  247.   free(NUsum1);
  248.   free(NUsum2);
  249.   free(DI1);
  250.   free(DI2);
  251.   free(NodePow);
  252.   free(NodeVol);
  253.   free(GenPower);
  254. }
  257. void Datain(int Nb,int Nl,int Ng,FILE *fp,struct Branch_Type *Branch,struct Load_Type *Load,struct Generator_Type *Generator,struct PVNode_Type *PVNode,int *Npv)
  258. {
  259.   int n;
  260.   /*********************从初始数据文件中读取支路、发电机、负荷数据************************/
  261.   for(n=1;n<=Nb;n++)
  262.   {
  263.     fscanf(fp,"%d,%d,%lf,%lf,%lf",&(Branch[n].i),&(Branch[n].j),&(Branch[n].R),&(Branch[n].X),&(Branch[n].YK));
  264.   }
  265.   for(n=1;n<=Ng;n++)
  266.   {
  267.     fscanf(fp,"%d,%lf,%lf,%lf",&(Generator[n].i),&(Generator[n].P),&(Generator[n].Q),&(Generator[n].V));
  268.     if((Generator[n].V)<0)
  269.     {
  270.         (*Npv)++;
  271.         PVNode[(*Npv)].i=Generator[n].i;
  272.         PVNode[(*Npv)].V=-(Generator[n].V);
  273.     }
  274.   }
  275.   for(n=1;n<=Nl;n++)
  276.   {
  277.     fscanf(fp,"%d,%lf,%lf,%lf",&(Load[n].i),&(Load[n].P),&(Load[n].Q),&(Load[n].V));
  278.     if((Load[n].V)<0)
  279.     {
  280.         (*Npv)++;
  281.         PVNode[(*Npv)].i=Load[n].i;
  282.         PVNode[(*Npv)].V=-(Load[n].V);
  283.     }
  284.   }
  285.   fclose(fp);            
  286. }
  289. void AdmittanceMatrix(int N,int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,int *NYseq1)
  290. {
  291.   /*****************不考虑接地支路时，形成节点导纳矩阵******************/
  292.   int i,n;
  293.   int *NYsum=malloc((N+1)*sizeof(int));
  294.   int *NYseq=malloc((N+1)*sizeof(int));
  295.   for(i=1;i<=N;i++)
  296.   {
  297.     Yii[i].G=0.0;
  298.     Yii[i].B=0.0;
  299.     Yiil[i].G=0.0;
  300.     Yiil[i].B=0.0;
  301.     NYsum[i]=0;
  302.   }
  303.   for(n=1;n<=Nb;n++)
  304.   {
  305.     int i=abs(Branch[n].i);
  306.     int j=abs(Branch[n].j);
  307.     double R=(Branch[n].R);
  308.     double X=(Branch[n].X);
  309.     double YK=(Branch[n].YK);
  311.     double Zmag2=R*R+X*X;
  312.     double Gij=R/Zmag2;
  313.     double Bij=-X/Zmag2;
  315.     double b_ij=-1.0/X;
  317.     if((Branch[n].i<0)||(Branch[n].j<0))
  318.     {
  319.         Yij[n].G=-Gij/YK;
  320.         Yij[n].B=-Bij/YK;
  321.         Yijl[n].G=0.0;
  322.         Yijl[n].B=-b_ij/YK;
  323.     }
  324.     else
  325.     {
  326.         Yij[n].G=-Gij;
  327.         Yij[n].B=-Bij;
  328.         Yijl[n].G=0.0;
  329.         Yijl[n].B=-b_ij;
  330.     }
  331.     Yij[n].j=j;
  332.     Yijl[n].j=j;
  334.     if((Branch[n].i<0)||(Branch[n].j<0))
  335.     {
  336.         Yii[i].G=Yii[i].G+Gij/YK;
  337.         Yii[i].B=Yii[i].B+Bij/YK;
  338.         Yii[j].G=Yii[j].G+Gij/YK;
  339.         Yii[j].B=Yii[j].B+Bij/YK;
  341.         Yiil[i].B=Yiil[i].B+b_ij/YK;
  342.         Yiil[j].B=Yiil[j].B+b_ij/YK;
  343.     }
  344.     else
  345.     {
  346.         Yii[i].G=Yii[i].G+Gij;
  347.         Yii[i].B=Yii[i].B+Bij;
  348.         Yii[j].G=Yii[j].G+Gij;
  349.         Yii[j].B=Yii[j].B+Bij;
  351.         Yiil[i].B=Yiil[i].B+b_ij;
  352.         Yiil[j].B=Yiil[j].B+b_ij;
  353.     }
  354.     NYsum[i]+=1;
  355.   }
  357.   NYseq[1]=1;
  358.   for(i=1;i<=N-1;i++)
  359.   {
  360.     NYseq[i+1]=NYseq[i]+NYsum[i];
  361.   }
  362.   for(n=1;n<=N;n++)
  363.   {
  364.     NYseq1[n]=NYseq[n];
  365.   }
  367. }
  369. void AdmittanceMatrixAdd(int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch)
  370. {
  371.   int n;
  372.   /***********************在导纳矩阵中追加接地支路************************/
  373.   for(n=1;n<=Nb;n++)
  374.   {
  375.     int i=Branch[n].i;
  376.     int j=Branch[n].j;
  377.     double YK=Branch[n].YK;
  379.     if(!(i<0||j<0))
  380.     {
  381.         double Bij=YK/2.0;
  382.         double b_ij=YK/2.0;
  384.         Yii[i].B=Yii[i].B+Bij;
  385.         Yii[j].B=Yii[j].B+Bij;
  387.         Yiil[i].B=Yiil[i].B+b_ij;
  388.         Yiil[j].B=Yiil[j].B+b_ij;
  389.     }
  390.     else
  391.     {
  392.         double Gij;
  393.         double Bij;
  394.         double b_ij;
  395.         if(i<0)
  396.         {
  397.             i=abs(i);
  398.             Gij=Yij[n].G;
  399.             Bij=Yij[n].B;
  401.             b_ij=Yijl[n].B;
  403.             Yii[i].G=Yii[i].G+(1.0-1.0/YK)*Gij;
  404.             Yii[i].B=Yii[i].B+(1.0-1.0/YK)*Bij;
  405.             Yiil[i].B=Yiil[i].B+(1.0-1.0/YK)*b_ij;
  408.             Yii[j].G=Yii[j].G+(1.0-YK)*Gij;
  409.             Yii[j].B=Yii[j].B+(1.0-YK)*Bij;
  410.             Yiil[j].B=Yiil[j].B+(1.0-YK)*b_ij;
  411.         }
  412.         else
  413.         {
  414.             j=abs(j);
  415.             Gij=Yij[n].G;
  416.             Bij=Yij[n].B;
  418.             b_ij=Yijl[n].B;
  420.             Yii[j].G=Yii[j].G+(1.0-1.0/YK)*Gij;
  421.             Yii[j].B=Yii[j].B+(1.0-1.0/YK)*Bij;
  422.             Yiil[j].B=Yiil[j].B+(1.0-1.0/YK)*b_ij;
  425.             Yii[i].G=Yii[i].G+(1.0-YK)*Gij;
  426.             Yii[i].B=Yii[i].B+(1.0-YK)*Bij;
  427.             Yiil[i].B=Yiil[i].B+(1.0-YK)*b_ij;
  428.         }
  429.     }
  430.   }
  432. }
  434. void Factorial(int flag,int N,int Npv,struct PVNode_Type *PVNode,int *NUsum,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq,double *D,struct U_Type *U)
  435. {
  437.   /**************************因子表形成函数***************************/
  438.   int n_pv,i_pv,j,n_u,i_above;
  439.   int i,count;
  440.   double *B=malloc((N+1)*sizeof(double));
  441.   double Btemp;
  442.   n_pv=1;
  443.   i_pv=PVNode[1].i;
  444.   
  445.   for(i=1;i<N;i++)
  446.   {
  447.     if((flag==2)&&(i==i_pv))
  448.     {
  449.         n_pv++;
  450.         i_pv=PVNode[n_pv].i;
  451.         NUsum[i]=0;
  452.         D[i]=0.0;
  453.         continue;
  454.     }
  455.     else
  456.     {
  457.         for(count=i+1;count<N;count++)
  458.             B[count]=0.0;
  459.         
  460.         B[i]=Yii[i].B;
  461.         
  462.         for(count=NYseq[i];count<NYseq[i+1];count++)
  463.         {
  464.             j=Yij[count].j;
  465.             B[j]=Yij[count].B;
  466.         }
  467.         
  468.         if(flag==2)
  469.         {
  470.             for(count=1;count<=Npv;count++)
  471.             {
  472.                 j=PVNode[count].i;
  473.                 B[j]=0.0;
  474.             }
  475.             
  476.         }
  478.         
  479.         n_u=1;
  480.         i_above=1;
  481.         
  482.         while(i_above<=(i-1))
  483.         {
  484.             count=1;
  485.             
  486.             while(count<=NUsum[i_above])
  487.             {
  488.                 if(U[n_u].j==i)
  489.                 {   
  490.                     
  491.                     Btemp=U[n_u].value/D[i_above];
  492.                     while(count<=NUsum[i_above])
  493.                     {
  494.                         j=U[n_u].j;
  495.                         B[j]=B[j]-Btemp*U[n_u].value;
  496.                         count++;
  497.                         n_u++;
  498.                     }
  499.                     break;
  500.                 }
  501.                 count++;
  502.                 n_u++;
  503.             }
  504.             i_above++;
  505.         }
  506.         
  507.         Btemp=1.0/B[i];        
  508.         D[i]=Btemp;
  510.         count=0;
  512.         for(j=i+1;j<N;j++)
  513.         {
  514.             if(B[j]!=0.0)
  515.             {   
  516.                 U[n_u].value=B[j]*Btemp;
  517.                 U[n_u].j=j;
  518.                 count++;
  519.                 n_u++;
  520.             }
  521.         }
  523.         NUsum[i]=count;
  524.     }
  525.   }
  526.  free(B);
  527. }
  529. void NodePower(int flag,int N,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq)
  530. {
  531. /*********************计算各节点功率函数**********************/
  532.   double A,B,Vi;
  533.   int i,n,j;
  534.   double VV,theta;
  535.   for(i=1;i<=N;i++)
  536.   {
  537.     if(flag==1)
  538.         NodePow[i].P=0.0;
  539.     else
  540.         NodePow[i].Q=0.0;
  541.   }
  543.   for(i=1;i<=N;i++)
  544.   {
  545.     Vi=NodeVol[i].V;
  547.     if(flag==1)
  548.     {
  549.         A=Yii[i].G;
  550.     }
  551.     else
  552.     {
  553.         A=-Yii[i].B;
  554.     }
  555.         
  556.     if(flag==1)
  557.         NodePow[i].P+=Vi*Vi*A;
  558.     else
  559.         {NodePow[i].Q+=Vi*Vi*A;}
  560.         
  561.     if(i==N)
  562.     {
  563.         break;
  564.     }
  565.     else
  566.     {
  567.         for(n=NYseq[i];n<=NYseq[i+1]-1;n++)
  568.         {
  569.             if(flag==1)
  570.             {
  571.                 A=Yij[n].G;
  572.                 B=Yij[n].B;
  573.             }
  574.             else
  575.             {
  576.                 A=-Yij[n].B;
  577.                 B=Yij[n].G;
  578.             }
  579.             j=Yij[n].j;
  581.             VV=Vi*NodeVol[j].V;
  582.             theta=NodeVol[i].theta-NodeVol[j].theta;
  583.             A=A*VV*cos(theta);
  584.             B=B*VV*sin(theta);
  585.             if(flag==1)
  586.             {
  587.                 NodePow[i].P+=(A+B);
  588.                 NodePow[j].P+=(A-B);
  589.             }
  590.             else
  591.             {
  592.                 NodePow[i].Q+=(A+B);
  593.                 NodePow[j].Q+=(A-B);
  594.             }
  595.         }
  596.     }
  597.   }
  598. }
  600. void Iteration(int flag,struct Generator_Type *Generator,struct Load_Type *Load,struct PVNode_Type *PVNode,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct GeneratorPower *GenPower,int N,double *DI,double *MaxError,int *ErrNode)
  601. {
  602.   /*****************迭代计算函数*****************/
  603.   int i=1,n_g=1,n_l=1,n_pv=1,i_g=Generator[1].i,i_l=Load[1].i,i_pv=PVNode[1].i;
  604.   double Vi,Wi,Wtemp;
  605.   (*MaxError)=0.0;
  606.     
  607.  do
  608.   {
  609.     Vi=NodeVol[i].V;
  610.     if(i==i_l)      
  611.     {
  612.         if(flag==1)
  613.         {
  614.             Wi=Load[n_l].P;
  615.         }
  616.         else
  617.         {
  618.             Wi=Load[n_l].Q;
  619.         }
  620.         n_l+=1;
  621.         i_l=Load[n_l].i;
  622.     }
  623.     else
  624.     {
  625.         Wi=0.0;
  626.     }
  627.     Wtemp=Wi;
  628.     if(flag==1)
  629.         Wi=Wi-NodePow[i].P;
  630.     else
  631.         Wi=Wi-NodePow[i].Q;
  633.     if(i==i_g)
  634.     {
  635.         if(flag==1)
  636.         {
  637.             NodePow[i].P=Wtemp;
  638.             GenPower[i_g].P=-Wi;
  639.         }
  640.         else
  641.         {
  642.             NodePow[i].Q=Wtemp;
  643.             GenPower[i_g].Q=-Wi;
  644.         }
  646.         if(flag==1)
  647.         {
  648.             Wi+=Generator[n_g].P;
  649.         }
  650.         else
  651.         {
  652.             Wi+=Generator[n_g].Q;
  653.         }
  654.         n_g+=1;
  655.         i_g=Generator[n_g].i;
  656.     }
  658.     if(i==N)
  659.     {
  660.         break;                 
  661.     }
  662.     else
  663.     {
  664.         if(flag==2&&i==i_pv)
  665.         {
  666.             n_pv+=1;
  667.             i_pv=PVNode[n_pv].i;
  668.             DI[i]=0.0;
  669.         }
  670.         else
  671.         {
  672.             if(fabs(Wi)>(*MaxError))
  673.             {
  674.                 (*MaxError)=fabs(Wi);
  675.                 (*ErrNode)=i;
  676.             }
  677.                 
  678.             DI[i]=Wi/Vi;
  679.         }
  680.     }
  681.     i+=1;
  682.   }while(1);
  683. }
  685. void FormulaSolution(int flag,struct U_Type *U,double *D,int *NUsum,double *DI,int N,struct NodalVol *NodeVol,double V0)
  686. {
  687.   /******************进行线性方程组的求解********************/
  688.   int n_u;
  689.   int i,count;
  690.   int j;
  691.   double DItemp,Dtheta,DV;
  692.   n_u=1;
  693.     
  694.   for(i=1;i<=N-1;i++)
  695.   {
  696.     DItemp=DI[i];
  697.     for(count=1;count<=NUsum[i];count++)
  698.     {
  699.         j=U[n_u].j;
  700.         DI[j]=DI[j]-DItemp*U[n_u].value;
  701.         n_u++;
  702.     }
  703.     DI[i]=DItemp*D[i];
  704.   }
  706.   for(i=N-1;i>=1;i--)
  707.   {
  708.     DItemp=DI[i];
  709.     for(count=1;count<=NUsum[i];count++)
  710.     {
  711.         n_u-=1;
  712.         j=U[n_u].j;
  713.         DItemp=DItemp-DI[j]*U[n_u].value;
  714.     }
  715.     DI[i]=DItemp;
  716.   }
  718.   for(i=1;i<=N-1;i++)
  719.   {
  720.     if(flag==1)
  721.     {
  722.         Dtheta=DI[i]/V0;
  723.         NodeVol[i].theta=NodeVol[i].theta-Dtheta;
  724.     }
  725.     else
  726.     {
  727.         DV=DI[i];
  728.         NodeVol[i].V-=DV;
  729.     }
  730.   }
  731. }
  733. void NodeDataOutput(FILE *fp,struct NodalVol *NodeVol,struct Generator_Type *Generator,int N,struct GeneratorPower *GenPower,struct NodalPow *NodePow,struct Load_Type *Load,int Nl)
  734. {
  735.   /***************节点数据输出**********************/
  736.   double Vmin=NodeVol[1].V;
  737.   double V,theta,P,Q;
  738.   int i_g=Generator[1].i;
  739.   int VminNode=1;
  740.   int n_g=1;
  741.   int i;
  742.     
  743.   for(i=1;i<=N;i++)
  744.   {
  745.     theta=NodeVol[i].theta/3.14159*180;
  746.     V=NodeVol[i].V;
  747.     if(V<Vmin)
  748.     {
  749.         Vmin=V;
  750.         VminNode=i;
  751.     }
  752.     if(i==i_g)
  753.     {
  754.         P=GenPower[i].P;
  755.         Q=GenPower[i].Q;
  756.         n_g+=1;
  757.         i_g=Generator[n_g].i;
  758.     }
  759.     else
  760.     {
  761.         P=0.0;
  762.         Q=0.0;
  763.     }
  764.     if(i!=N)
  765.         fprintf(fp,"  %dt   %10.7lft   %10.7lft  %10.7lft  %10.7lfn",i,V,theta,P,Q);
  766.     else
  767.         fprintf(fp,"  %dt   %10.7lft   %10.7lft  %10.7lft  %10.7lfn",i,V,theta,NodePow[i].P-Load[Nl].P,NodePow[i].Q-Load[Nl].Q);
  769.   }
  770.     fprintf(fp,"系统最低电压=%10.7lf,节点=%dn",Vmin,VminNode);
  771. }
  773. void BranchDataOutput(FILE *fp,int Nb,struct Branch_Type *Branch,struct NodalVol *NodeVol)
  774. {
  775.   /***************支路数据输出***************/
  776.   double PLoss=0.0,QLoss=0.0;
  777.   int n;
  778.   int i,j;
  779.   double R,X,YK,Y,theta,Ei,Ej,Fi,Fj,Vi,Vj,DE,DF;
  780.   double Zmag2,Ir,Ii;
  781.   double Pij,Qij,Pji,Qji;
  782.   for(n=1;n<=Nb;n++)
  783.   {
  784.     i=abs(Branch[n].i);
  785.     j=abs(Branch[n].j);
  786.     R=Branch[n].R;
  787.     X=Branch[n].X;
  788.     YK=Branch[n].YK;
  790.     Vi=NodeVol[i].V;
  791.     theta=NodeVol[i].theta;
  792.     Ei=Vi*cos(theta);
  793.     Fi=Vi*sin(theta);
  795.     Vj=NodeVol[j].V;
  796.     theta=NodeVol[j].theta;
  797.     Ej=Vj*cos(theta);
  798.     Fj=Vj*sin(theta);
  799.             
  800.     if(Branch[n].i<0||Branch[n].j<0)
  801.     {
  802.         if(Branch[n].i<0)
  803.         {
  804.             Ei=Ei/YK;
  805.             Fi=Fi/YK;
  806.         }
  807.         else
  808.         {
  809.             Ej=Ej/YK;
  810.             Fj=Fj/YK;
  811.         }
  812.         YK=0.0;
  813.     }
  814.             
  815.     DE=Ei-Ej;
  816.     DF=Fi-Fj;
  817.     Zmag2=R*R+X*X;
  818.     Ir=(DE*R+DF*X)/Zmag2;
  819.     Ii=(DF*R-DE*X)/Zmag2;
  820.             
  821.     Pij=Ir*Ei+Ii*Fi;
  822.     Qij=Ir*Fi-Ii*Ei;
  823.             
  824.     Pji=-(Ir*Ej+Ii*Fj);
  825.     Qji=-(Ir*Fj-Ii*Ej);
  826.             
  827.     Qij-=(Vi*Vi*YK/2.0);
  828.     Qji-=(Vj*Vj*YK/2.0);
  829.             
  830.     PLoss=PLoss+Pij+Pji;
  831.     QLoss=QLoss+Qij+Qji;
  832.             
  833.     fprintf(fp,"  %3d->%3dt   %10.7lft   %10.7lfn  %3d->%3dt   %10.7lft   %10.7lfn",i,j,Pij,Qij,j,i,Pji,Qji);
  834.   }
  836.   fprintf(fp,"     损耗t   %10.7lft   %10.7lfn",PLoss,QLoss);
  837. }