powerflow.cpp

电力系统牛顿拉夫逊法的潮流程序

    
  GenPower=malloc((N+1)*sizeof(struct GeneratorPower)); /*为发电机所发功率数组分配内存*/  
  fprintf(fp,"\t系统潮流计算结果：\n(1)迭代过程纪录：\n迭代次数\t    \t有功迭代\t\t无功迭代\n\t\t   ΔPmax\tP-Node\t ΔQmax\t\tQ-Node\n");
  for(k=0;1;k++)
  {
    fprintf(fp,"   %2d:\t",k);
    if(Kp==1)
    {
        NodePower(1,N,NodeVol,NodePow,Yii,Yij,NYseq);/*节点功率计算函数*/
        Iteration(1,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI1,&MaxError,&ErrNode);/*迭代*/
        fprintf(fp,"\t%10.7lf\t %d\t",MaxError,ErrNode);
        if(MaxError>=epsilon)
            FormulaSolution(1,U1,D1,NUsum1,DI1,N,NodeVol,V0);/*线性方程求解*/
        else
            Kp=0;
    }
    else
        fprintf(fp,"\t\t\t\t");
    if(Kq==1)
    {
        NodePower(2,N,NodeVol,NodePow,Yii,Yij,NYseq);/*节点功率计算函数*/
        Iteration(2,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI2,&MaxError,&ErrNode);/*迭代*/
        fprintf(fp,"%10.7lf\t %d\n",MaxError,ErrNode);
        if(MaxError>=epsilon)
            FormulaSolution(2,U2,D2,NUsum2,DI2,N,NodeVol,V0);/*线性方程求解*/
        else
            Kq=0;

    }
    else
        fprintf(fp,"\n");
    if(Kp==0&&Kq==0)
        break;
    if(k>2000)
    {
        fprintf(fp,"\n迭代次数超过2000次，系统不收敛!\n");
        fprintf(fp,"本系统为%d个节点、%d条支路系统，含有%d个补偿电容和%d个并联电抗\n",N,Nb,Nc,Nr);
        now=time(NULL);
        fprintf(fp,"\t\t处理时间： %s\n",ctime(&now));
        /*fclose(fp);
        exit(1);*/
        goto DATA;
    }

    
    printf("COMPUTING......\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
  }


  fprintf(fp,"\n总计迭代次数: %d 次!\n\n",k+1);
  fprintf(fp,"(2)潮流计算结果(节点电压)：\n Node\t\tV\t\tθ\t\tP\t\tQ\n");
  NodeDataOutput(fp,NodeVol,Generator,N,GenPower,NodePow,Load,Nl);
  fprintf(fp,"\n(3)潮流计算结果(支路功率):\n    Branch\t\tP\t\tQ\t\tQc\t\tQr\n");
  BranchDataOutput(fp,Nb,Nc,Nr,Branch,Compensation,Reactance,NodeVol);
  fprintf(fp,"\n本系统为%d个节点、%d条支路系统，含有%d个补偿电容和%d个并联电抗\n",N,Nb,Nc,Nr);
  now=time(NULL);
  fprintf(fp,"\n\t\t\t处理时间： %s\n",ctime(&now));

DATA:
  fprintf(fp,"\n\n附原始数据：\n");
  fprintf(fp,"节点数=%d 个,支路数=%d 条,补偿电容=%d 个,并联电抗=%d 个，发电机=%d 个,负荷=%d 个,系统平均电压=%lf,误差=%lf\n",N,Nb,Nc,Nr,Ng,Nl,V0,epsilon);
  for(n=1;n<=Nb;n++)
  {
    if((Branch[n].i>0)&&(Branch[n].j>0))
    {
        fprintf(fp,"节点%3d和节点%3d之间的电力线参数:",Branch[n].i,Branch[n].j);
        fprintf(fp,"R+jX=%e+j%e   Y0=%e\n",Branch[n].R,Branch[n].X,Branch[n].YK);
    }
    else
    {
        fprintf(fp,"节点%3d和节点%3d之间的变压器参数:",abs(Branch[n].i),abs(Branch[n].j));
        fprintf(fp,"R+jX=%e+j%e   K=%lf,",Branch[n].R,Branch[n].X,fabs(Branch[n].YK));
    if(Branch[n].i<0)
        i=abs(Branch[n].i);
    else
        i=abs(Branch[n].j);
    fprintf(fp,"非标准变比位于节点%3d侧\n",i);
    }
  }
  for(n=1;n<=Nc;n++)
  {
    fprintf(fp,"位于节点%3d的补偿电容:",Compensation[n].i);
    fprintf(fp,"Y=j%e\n",Compensation[n].Y);
  }
  for(n=1;n<=Nr;n++)
  {
    fprintf(fp,"位于节点%3d的并联电抗:",Reactance[n].i);
    fprintf(fp,"X=j%e\n",Reactance[n].X);
  }
  for(n=1;n<=Nl;n++)
  {
    fprintf(fp,"节点%3d的负荷  :",Load[n].i);
    fprintf(fp,"P+jQ=%lf-j%lf   V=%lf \n",Load[n].P,fabs(Load[n].Q),Load[n].V);
  }
  for(n=1;n<=Ng;n++)
  {
    fprintf(fp,"节点%3d的发电机:",Generator[n].i);
    fprintf(fp,"P+jQ=%lf+j%lf   V=%lf \n",Generator[n].P,Generator[n].Q,Generator[n].V);
  }
  for(n=1;n<Npv;n++)
  {
    if(n==1)
    fprintf(fp,"PV节点有%3d个,分别为：",Npv-1);
    fprintf(fp,"节点%3d,",PVNode[n].i); 
  }

  fprintf(fp,"\n\n\n");

  fclose(fp);
  free(Branch);
  free(Compensation);
  free(Reactance);
  free(Generator);
  free(Load);
  free(PVNode);
  free(Yii);
  free(Yiil);
  free(Yij);
  free(Yijl);
  free(NYseq);
  free(U1);
  free(U2);
  free(D1);
  free(D2);
  free(NUsum1);
  free(NUsum2);
  free(DI1);
  free(DI2);
  free(NodePow);
  free(NodeVol);
  free(GenPower);
}



void Datain(int Nb,int Nl,int Nc,int Nr,int Ng,FILE *fp,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance,struct Load_Type *Load,struct Generator_Type *Generator,struct PVNode_Type *PVNode,int *Npv)
{
  int n;
  /**********下面从数据文件中读取支路、发电机、负荷的数据***********/

  for(n=1;n<=Nb;n++)
  {
    fscanf(fp,"%d,%d,%lf,%lf,%lf",&(Branch[n].i),&(Branch[n].j),&(Branch[n].R),&(Branch[n].X),&(Branch[n].YK));
  }
  for(n=1;n<=Nc;n++)
  {
    fscanf(fp,"%d,%lf",&(Compensation[n].i),&(Compensation[n].Y));
  }
  for(n=1;n<=Nr;n++)
  {
    fscanf(fp,"%d,%lf",&(Reactance[n].i),&(Reactance[n].X));
  }
  for(n=1;n<=Ng;n++)
  {
    fscanf(fp,"%d,%lf,%lf,%lf",&(Generator[n].i),&(Generator[n].P),&(Generator[n].Q),&(Generator[n].V));
    if((Generator[n].V)<0)
    {
        (*Npv)++;
        PVNode[(*Npv)].i=Generator[n].i;
        PVNode[(*Npv)].V=-(Generator[n].V);
    }
  }
  for(n=1;n<=Nl;n++)
  {
    fscanf(fp,"%d,%lf,%lf,%lf",&(Load[n].i),&(Load[n].P),&(Load[n].Q),&(Load[n].V));
    if((Load[n].V)<0)
    {
        (*Npv)++;
        PVNode[(*Npv)].i=Load[n].i;
        PVNode[(*Npv)].V=-(Load[n].V);
    }
  }



  fclose(fp);/*关闭数据输入流*/
/******************读取支路、发电机、负荷的数据完成，并关闭系统数据文件********************/



}







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,int *NYsum1)
{
  int i,n;


  int *NYsum=malloc((N+1)*sizeof(int));
  int *NYseq=malloc((N+1)*sizeof(int));
  for(i=1;i<=N;i++)
  {
    Yii[i].G=0.0;
    Yii[i].B=0.0;
    Yiil[i].G=0.0;
    Yiil[i].B=0.0;
    NYsum[i]=0;
  }

/************下面形成不考虑接地支路时的导纳矩阵*******************/
  for(n=1;n<=Nb;n++)
  {
    int i=abs(Branch[n].i);
    int j=abs(Branch[n].j);
    double R=(Branch[n].R);
    double X=(Branch[n].X);
    double YK=(Branch[n].YK);

    double Zmag2=R*R+X*X;
    double Gij=R/Zmag2;
    double Bij=-X/Zmag2;

    double b_ij=-1.0/X;

    if((Branch[n].i<0)||(Branch[n].j<0))
    {
        Yij[n].G=-Gij/YK;
        Yij[n].B=-Bij/YK;
        Yijl[n].G=0.0;
        Yijl[n].B=-b_ij/YK;
    }
    else
    {
        Yij[n].G=-Gij;
        Yij[n].B=-Bij;
        Yijl[n].G=0.0;
        Yijl[n].B=-b_ij;

        /*这几句在最后一次执行时，执行ps.dat时NYsum[1]、[2]变为0，NYsum[3]=-21845,NYsum[4]=-21845,但是[1]、[2]、[3]均应为2，[4]应为0;在执行ps-5.dat时则没有这问题*//*现在通过在此子程序内定义新的NYsum和NYseq已经得到了解决*/

    }
    Yij[n].j=j;
    Yijl[n].j=j;

    if((Branch[n].i<0)||(Branch[n].j<0))
    {
        Yii[i].G=Yii[i].G+Gij/YK;
        Yii[i].B=Yii[i].B+Bij/YK;
        Yii[j].G=Yii[j].G+Gij/YK;
        Yii[j].B=Yii[j].B+Bij/YK;

        Yiil[i].B=Yiil[i].B+b_ij/YK;
        Yiil[j].B=Yiil[j].B+b_ij/YK;
    }
    else
    {
        Yii[i].G=Yii[i].G+Gij;
        Yii[i].B=Yii[i].B+Bij;
        Yii[j].G=Yii[j].G+Gij;
        Yii[j].B=Yii[j].B+Bij;

        Yiil[i].B=Yiil[i].B+b_ij;
        Yiil[j].B=Yiil[j].B+b_ij;
    }
    NYsum[i]+=1;

  }

  NYseq[1]=1;

/****************这里给各行的首元素定位********************/

  for(i=1;i<=N-1;i++)
  {
    NYseq[i+1]=NYseq[i]+NYsum[i];
  }

  for(n=1;n<=N;n++)
 {
  NYseq1[n]=NYseq[n];
 }




}


void AdmittanceMatrixAdd(int Nb,int Nc,int Nr,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance)
{
  int n;
  /*char ch;*/




/*******************下面在导纳矩阵中追加接地支路*********************/

  for(n=1;n<=Nb;n++)
  {
    int i=Branch[n].i;
    int j=Branch[n].j;
    double YK=Branch[n].YK;

    if(!(i<0||j<0))
    {
        double Bij=YK/2.0;
        double b_ij=YK/2.0;

        Yii[i].B=Yii[i].B+Bij;
        Yii[j].B=Yii[j].B+Bij;

        Yiil[i].B=Yiil[i].B+b_ij;
        Yiil[j].B=Yiil[j].B+b_ij;
    }
    else
    {
        double Gij;
        double Bij;
        double b_ij;
        if(i<0)
        {
            i=abs(i);
            Gij=Yij[n].G;
            Bij=Yij[n].B;

            b_ij=Yijl[n].B;

            Yii[i].G=Yii[i].G+(1.0-1.0/YK)*Gij;
            Yii[i].B=Yii[i].B+(1.0-1.0/YK)*Bij;
            Yiil[i].B=Yiil[i].B+(1.0-1.0/YK)*b_ij;


            Yii[j].G=Yii[j].G+(1.0-YK)*Gij;
            Yii[j].B=Yii[j].B+(1.0-YK)*Bij;
            Yiil[j].B=Yiil[j].B+(1.0-YK)*b_ij;
        }
        else
        {
            j=abs(j);
            Gij=Yij[n].G;
            Bij=Yij[n].B;

            b_ij=Yijl[n].B;

            Yii[j].G=Yii[j].G+(1.0-1.0/YK)*Gij;
            Yii[j].B=Yii[j].B+(1.0-1.0/YK)*Bij;
            Yiil[j].B=Yiil[j].B+(1.0-1.0/YK)*b_ij;


            Yii[i].G=Yii[i].G+(1.0-YK)*Gij;
            Yii[i].B=Yii[i].B+(1.0-YK)*Bij;
            Yiil[i].B=Yiil[i].B+(1.0-YK)*b_ij;
        }
    }
  }
/********************追加接地支路完成**********************/
/********************追加对地补偿电容***********************/
  for(n=1;n<=Nc;n++)
  {
    int i=Compensation[n].i;
    double Bij=Compensation[n].Y;
    double b_ij=Compensation[n].Y;

    Yii[i].B=Yii[i].B+Bij;
    
    Yiil[i].B=Yiil[i].B+b_ij;
   }
/********************追加对地并联电感***********************/
  for(n=1;n<=Nr;n++)
  {
    int i=Reactance[n].i;
    double Bij=-1/Reactance[n].X;
    double b_ij=-1/Reactance[n].X;

    Yii[i].B=Yii[i].B+Bij;
    
    Yiil[i].B=Yiil[i].B+b_ij;
    
  }



}



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)/*因子表形成函数*/
{
  int n_pv,i_pv,j,n_u,i_above;
  int i,count;
  double *B=malloc((N+1)*sizeof(double));
  double Btemp;
  
  n_pv=1;
  i_pv=PVNode[1].i;
  
  for(i=1;i<N;i++)
  {
    if((flag==2)&&(i==i_pv))
    {
        n_pv++;
        i_pv=PVNode[n_pv].i;
        NUsum[i]=0;
        D[i]=0.0;
        continue;
    }
    else
    {
        for(count=i+1;count<N;count++)
            B[count]=0.0;
        
        B[i]=Yii[i].B;
        
        for(count=NYseq[i];count<NYseq[i+1];count++)
        {
            j=Yij[count].j;
            B[j]=Yij[count].B;
        }
        
        if(flag==2)
        {
            for(count=1;count<=Npv;count++)
            {
                j=PVNode[count].i;
                B[j]=0.0;
            }
            
        }

        
        n_u=1;
        i_above=1;
        
        while(i_above<=(i-1))
        {
            count=1;
            
            while(count<=NUsum[i_above])
            {
                if(U[n_u].j==i)
                {   
                    
                    Btemp=U[n_u].value/D[i_above];