powerflow.cpp

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

                    while(count<=NUsum[i_above])
                    {
                        j=U[n_u].j;
                        B[j]=B[j]-Btemp*U[n_u].value;
                        count++;
                        n_u++;
                    }
                    break;
                }
                count++;
                n_u++;
            }
            i_above++;
        }
        
        Btemp=1.0/B[i];
        
        D[i]=Btemp;
        count=0;
        for(j=i+1;j<N;j++)
        {
            if(B[j]!=0.0)
            {   
                U[n_u].value=B[j]*Btemp;
                U[n_u].j=j;
                count++;
                n_u++;
            }
        }
        NUsum[i]=count;
    }
  }
 free(B);
}










void NodePower(int flag,int N,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq)
{
    /*计算节点功率*/

  double A,B,Vi;
  int i,n,j;
  double VV,theta;
  for(i=1;i<=N;i++)
  {
    if(flag==1)
        NodePow[i].P=0.0;
    else
        NodePow[i].Q=0.0;
  }
    
    

  for(i=1;i<=N;i++)
  {
    Vi=NodeVol[i].V;

    if(flag==1)
    {
        A=Yii[i].G;
    }
    else
    {
        A=-Yii[i].B;
    }
        
    if(flag==1)
        NodePow[i].P+=Vi*Vi*A;
    else
        {NodePow[i].Q+=Vi*Vi*A;}
        
    if(i==N)
    {
        break;
    }
    else
    {
        
        for(n=NYseq[i];n<=NYseq[i+1]-1;n++)
        {
            if(flag==1)
            {
                A=Yij[n].G;
                B=Yij[n].B;
            }
            else
            {
                A=-Yij[n].B;
                B=Yij[n].G;
            }
            j=Yij[n].j;
            VV=Vi*NodeVol[j].V;
            theta=NodeVol[i].theta-NodeVol[j].theta;
            A=A*VV*cos(theta);
            B=B*VV*sin(theta);
            if(flag==1)
            {
                NodePow[i].P+=(A+B);
                NodePow[j].P+=(A-B);
            }
            else
            {
                NodePow[i].Q+=(A+B);
                NodePow[j].Q+=(A-B);
            }
        }
    }
    
  }
    
}







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)
{
    /*迭代函数*/
   

  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;
  double Vi,Wi,Wtemp;
    

  (*MaxError)=0.0;

    

    
  do
  {
    Vi=NodeVol[i].V;

    if(i==i_l)      
    {
        if(flag==1)
        {
            Wi=Load[n_l].P;
        }
        else
        {
            Wi=Load[n_l].Q;
        }
        n_l+=1;
        i_l=Load[n_l].i;
    }
    else
    {
        Wi=0.0;
    }
    Wtemp=Wi;
    if(flag==1)
        Wi=Wi-NodePow[i].P;
    else
        Wi=Wi-NodePow[i].Q;


    if(i==i_g)
    {
        if(flag==1)
        {
            NodePow[i].P=Wtemp;
            GenPower[i_g].P=-Wi;
        }
        else
        {
            NodePow[i].Q=Wtemp;
            GenPower[i_g].Q=-Wi;
        }


        if(flag==1)
        {
            Wi+=Generator[n_g].P;
        }
        else
        {
            Wi+=Generator[n_g].Q;
        }
        n_g+=1;
        i_g=Generator[n_g].i;
    }
        

    if(i==N)
    {
        break;/*退出循环，进行MaxError与epsilon的比较*/
    }
    else
    {
        if(flag==2&&i==i_pv)
        {
            n_pv+=1;
            i_pv=PVNode[n_pv].i;
            DI[i]=0.0;
        }
        else
        {
            if(fabs(Wi)>(*MaxError))
            {
                (*MaxError)=fabs(Wi);
                (*ErrNode)=i;
            }
                
            DI[i]=Wi/Vi;
        }
        
    }i+=1;

  }while(1);
    
    
}





void FormulaSolution(int flag,struct U_Type *U,double *D,int *NUsum,double *DI,int N,struct NodalVol *NodeVol,double V0)
{

    /*迭代完成后，进行线性方程组的求解*/
  int n_u;
  int i,count;
  int j;
  double DItemp,Dtheta,DV;
  n_u=1;
    
  for(i=1;i<=N-1;i++)
  {
    DItemp=DI[i];
    for(count=1;count<=NUsum[i];count++)
    {
        j=U[n_u].j;
        DI[j]=DI[j]-DItemp*U[n_u].value;
        n_u++;
    }
    DI[i]=DItemp*D[i];
  }
  for(i=N-1;i>=1;i--)
  {
    DItemp=DI[i];
    for(count=1;count<=NUsum[i];count++)
    {
        n_u-=1;
        j=U[n_u].j;
        DItemp=DItemp-DI[j]*U[n_u].value;
    }
    DI[i]=DItemp;
  }




  for(i=1;i<=N-1;i++)
  {
    if(flag==1)
    {
        Dtheta=DI[i]/V0;
        NodeVol[i].theta=NodeVol[i].theta-Dtheta;
    }
    else
    {
        DV=DI[i];
        NodeVol[i].V-=DV;
    }
  }

 

}



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)
{
    /*节点数据输出*/
  double Vmin=NodeVol[1].V;
  double V,theta,P,Q,temp;
  int i_g=Generator[1].i;
  int VminNode=1;
  int n_g=1;
  int i;
    
  for(i=1;i<=N;i++)
  {
    theta=NodeVol[i].theta/3.14159*180;
    if(theta>180||theta<-180)
    {
        theta=(int)theta%360+modf(theta,&temp);
        if(theta>180)
            theta-=360;
        else
        if(theta<-180)
            theta+=360;
    }
    V=NodeVol[i].V;
    if(V<Vmin)
    {
        Vmin=V;
        VminNode=i;
    }
    else
    {
        ;
    }
    if(i==i_g)
    {
        P=GenPower[i].P;
        Q=GenPower[i].Q;
        n_g+=1;
        i_g=Generator[n_g].i;
    }
    else
    {
        P=0.0;
        Q=0.0;
    }
    if(i!=N)
        fprintf(fp,"  %d\t   %10.7lf\t   %10.7lf\t  %10.7lf\t  %10.7lf\n",i,V,theta,P,Q);
    else
        fprintf(fp,"  %d\t   %10.7lf\t   %10.7lf\t  %10.7lf\t  %10.7lf\n",i,V,theta,NodePow[i].P-Load[Nl].P,NodePow[i].Q-Load[Nl].Q);

  }
    fprintf(fp,"系统最低电压=%10.7lf,节点=%d\n",Vmin,VminNode);
}



void BranchDataOutput(FILE *fp,int Nb,int Nc,int Nr,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance,struct NodalVol *NodeVol)
{
  /*支路数据输出*/
  double PLoss=0.0,QLoss=0.0;
  int n;
  int i,j;
  double R,X,YK,Y,theta,Ei,Ej,Fi,Fj,Vi,Vj,DE,DF;
  double Zmag2,Ir,Ii;
  double Pij,Qij,Pji,Qji;
  for(n=1;n<=Nb;n++)
  {
    i=abs(Branch[n].i);
    j=abs(Branch[n].j);
    R=Branch[n].R;
    X=Branch[n].X;
    YK=Branch[n].YK;

    Vi=NodeVol[i].V;
    theta=NodeVol[i].theta;
    Ei=Vi*cos(theta);
    Fi=Vi*sin(theta);

    Vj=NodeVol[j].V;
    theta=NodeVol[j].theta;
    Ej=Vj*cos(theta);
    Fj=Vj*sin(theta);
            
    if(Branch[n].i<0||Branch[n].j<0)
    {
        if(Branch[n].i<0)
        {
            Ei=Ei/YK;
            Fi=Fi/YK;
        }
        else
        {
            Ej=Ej/YK;
            Fj=Fj/YK;
        }
        YK=0.0;
    }
    else
    {
        ;
    }
            
    DE=Ei-Ej;
    DF=Fi-Fj;
    Zmag2=R*R+X*X;
    Ir=(DE*R+DF*X)/Zmag2;
    Ii=(DF*R-DE*X)/Zmag2;
            
    Pij=Ir*Ei+Ii*Fi;
    Qij=Ir*Fi-Ii*Ei;
            
    Pji=-(Ir*Ej+Ii*Fj);
    Qji=-(Ir*Fj-Ii*Ej);
            
    Qij-=(Vi*Vi*YK/2.0);
    Qji-=(Vj*Vj*YK/2.0);
            
    PLoss=PLoss+Pij+Pji;
    QLoss=QLoss+Qij+Qji;
            
    fprintf(fp,"  %3d->%3d\t   %10.7lf\t   %10.7lf\n  %3d->%3d\t   %10.7lf\t   %10.7lf\n",i,j,Pij,Qij,j,i,Pji,Qji);
  }
  for(n=1;n<=Nc;n++)
  {
    i=Compensation[n].i;
    Y=Compensation[n].Y;

    Vi=NodeVol[i].V;

    Qij=-(Vi*Vi*Y); 

    QLoss=QLoss+Qij;

    fprintf(fp,"  %3d->0  \t\t\t\t\t   %10.7lf\n",i,Qij);
  }
  for(n=1;n<=Nr;n++)
  {
    i=Reactance[n].i;
    X=Reactance[n].X;

    Vi=NodeVol[i].V;

    Qij=Vi*Vi/X;   

    QLoss=QLoss+Qij;

    fprintf(fp,"  %3d->0  \t\t\t\t\t\t\t   %10.7lf\n",i,Qij);
  }


  fprintf(fp,"     损耗\t   %10.7lf\t   %10.7lf\n",PLoss,QLoss);
         
}