homotloa.cpp

用于电力系统潮流计算 c++程序

/* Homotopy Continuation Method: Load variables */

#include "homot.h"

/* ------- Global Variables ------ */
extern VALUETYPE *Dx,Dparam,param0,*x0,*x0p,Kh,Htol,SD0,AngTr,
                 VoltageSum,VoltageSum0,DxiMax,VSF,SF,ZeroDx,Tol;
extern INDEX NewNumEq,CountSteps,NumSteps;
extern AClist *Vlist,*Vlistp;
extern BOOLEAN flagReducedContinuation,flagReduceSystem,*DxZero,flagBS;


/* ------------------ LoadX0 ----------------------------- */
#ifdef ANSIPROTO
VALUETYPE LoadX0(BOOLEAN FlagLoadX0,BOOLEAN FlagUpdateVar,BOOLEAN FlagMakeDxZero)
#else
VALUETYPE LoadX0(FlagLoadX0,FlagUpdateVar,FlagMakeDxZero)
BOOLEAN FlagLoadX0,FlagUpdateVar,FlagMakeDxZero;
#endif
/* Load x0 vector and update AC/DC variables for power flow. */
{
  ACbusData *ACptr,*ACptrp,*Ptr,*BEptr;
  AClist *ALptr;
  DCbusData *DCptrR,*DCptrI,*DCptr;
  SVCbusData *SVCptr;                    /* FACTS */
  TCSCbusData *TCSCptr;                  /* FACTS */
  STATCOMbusData *STATCOMptr;            /* FACTS */
  ElementData *Eptr;
  ElementList *ELptr;
  VALUETYPE DPg,val,valp,vals,count,consp=0,Pg,Pmax,PgMax,Q,Qm;
  INDEX i,j;
  BOOLEAN Recover=TRUE;
  char Qmax[5],Qmin[5];

  if (ExistParameter('G')) Recover=FALSE;
  NewNumEq=Jac->n1 - 1 ;
  if (FlagLoadX0) VoltageSum0=Tol;
  if (FlagUpdateVar) VoltageSum=0;
  if (!Bl) {
    if (FlagLoadX0) param0=lambda;
    if (FlagUpdateVar) {
      lambda=param0+Dparam;
      if (lambda<=0) {
        val=fabs(lambda/Dparam);
        if(ExistParameter('d')) fprintf(stderr,"lambdaMin  %lf\n",val);
        lambda=0;}
      else val=0;
      if (val>consp) consp=val;
    }
  }
  for(ALptr=dataPtr->KGbus;ALptr!=NULL;ALptr=ALptr->Next) {
    ACptr=ALptr->AC;
    if (strpbrk(ACptr->Type,"S")) {
      i=ACvar[ACptr->N];
      if (FlagLoadX0)  {
        x0[i]=ACptr->Kg;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Kg=x0[i]+Dx[i];
      BEptr=ACptr;
    }
    else if(Acont){
      i=ACvar[ACptr->N]+2;
      if (FlagLoadX0) {
        x0[i]=ACptr->Kg;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Kg=x0[i]+Dx[i];
    }
  }
  for (Ptr=NULL,valp=0,count=0,ACptr=dataPtr->ACbus; ACptr!=NULL; ACptr=ACptr->Next){
    if (FlagUpdateVar) {
      if (ACptr->Area!=NULL) BEptr=ACptr->Area->Slack;
      DPg=ACptr->DPg;
      if (flagBS) {
        if (ACptr!=BEptr) Pg=ACptr->Pg+lambda*DPg;
        else              Pg=ACptr->Pg+BEptr->Kg;
      }
      else  Pg=ACptr->Pg+BEptr->Kg*DPg;
      Pmax=ACptr->Smax*ACptr->Smax-ACptr->Qg*ACptr->Qg;
      if (Pmax>0) Pmax=sqrt(Pmax);
      else        Pmax=99999999.;
      if (!flagSmax && Pmax<ACptr->PgMax) PgMax=Pmax;
      else if (!flagPgMax)                PgMax=ACptr->PgMax;
      else                                PgMax=99999999.;
      if (Pg>PgMax) {
        Pg=PgMax;
        DPg=0;
      }
      ACptr->PG=Pg;
      ACptr->DPG=DPg;
      ACptr->Pmax=PgMax;
      if (ACptr->Qmax==ACptr->Max) strcpy(Qmax,"Qmax");
      else                         strcpy(Qmax,"Smax");
      if (ACptr->Qmin==ACptr->Min) strcpy(Qmin,"Qmin");
      else                         strcpy(Qmin,"Smin");
    }
    if (!strcmp(ACptr->Type,"B")||!strcmp(ACptr->Type,"BA")) {
      i=ACvar[ACptr->N];
      if (FlagLoadX0) {
        x0[i]=ACptr->Ang;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Ang=x0[i]+Dx[i];
      i++;
      if (FlagLoadX0) {
        x0[i]=ACptr->V;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
        if (InList(ACptr,Vlistp)) VoltageSum0+=ACptr->V;
      }
      if (FlagUpdateVar) {
        ACptr->V=x0[i]+Dx[i];
        if (ACptr->V<=0){
          count++; val=fabs((ACptr->V-0.00001)/Dx[i]);
          if(ExistParameter('d')) fprintf(stderr,"%s %d %s Vzero %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
          ACptr->V=0.00001;}
        else val=0;
        if (val>consp) consp=val;
        if (InList(ACptr,Vlistp)) VoltageSum+=ACptr->V;
      }
      if(!Bl && FlagLoadX0 && FlagUpdateVar &&  FlagMakeDxZero && fabs(Dx[i])>valp) {
        valp=fabs(Dx[i]*param0/x0[i]);
        Ptr=ACptr;
      }
    }
    else if(strpbrk(ACptr->Type,"L")){
      i=ACvar[ACptr->N];
      if (FlagLoadX0) {
        x0[i]=ACptr->Ang;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Ang=x0[i]+Dx[i];
      i++; if (FlagLoadX0) x0[i]=lambda;
      if (FlagUpdateVar) {
        lambda=x0[i]+Dx[i];
        if (lambda<=0) {
          val=fabs(lambda/Dx[i]);
          if(ExistParameter('d')) fprintf(stderr,"lambdaMin  %lf\n",val);
          lambda=0;}
        else val=0;
        if (val>consp) consp=val;
      }
      if (FlagLoadX0) {
        param0=ACptr->V;
        if (InList(ACptr,Vlistp)) VoltageSum0+=ACptr->V;
      }
      if (FlagUpdateVar) {
        ACptr->V=param0+Dparam;
        if (ACptr->V<=0){
          count++; val=fabs((ACptr->V-0.00001)/Dx[i]);
          if(ExistParameter('d')) fprintf(stderr,"%s %d %s Vzero %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
          ACptr->V=0.00001;}
        else val=0;
        if (val>consp) consp=val;
        if (InList(ACptr,Vlistp)) VoltageSum+=ACptr->V;
        if(fabs(Dparam)>param0) {
          val=fabs(1-param0/fabs(Dparam)/8.0);
          if(ExistParameter('d')) fprintf(stderr,"%s %d %s lambdaMin %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);}
        else val=0;
        if (val>consp) consp=val;
      }
      if(FlagLoadX0 && FlagUpdateVar &&  FlagMakeDxZero && fabs(Dx[i])>valp) {
        valp=fabs(Dx[i]*param0/x0[i]);
        Ptr=ACptr;
      }
    }
    else if(strpbrk(ACptr->Type,"C")){
      i=ACvar[ACptr->N];
      if (FlagLoadX0) {
        x0[i]=ACptr->Ang;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Ang=x0[i]+Dx[i];
      i++;
      if (!QRcont) {
        if (FlagLoadX0) {
          x0[i]=ACptr->V;
          if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
          if (InList(ACptr,Vlistp)) VoltageSum0+=ACptr->V;
        }
        if (FlagUpdateVar) {
          ACptr->V=x0[i]+Dx[i];
          if (ACptr->V<=0){
            count++; val=fabs((ACptr->V-0.00001)/Dx[i]);
            if(ExistParameter('d')) fprintf(stderr,"%s %d %s Vzero %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
            ACptr->V=0.00001;}
          else val=0;
          if (val>consp) consp=val;
          if (InList(ACptr,Vlistp)) VoltageSum+=ACptr->V;
        }
        if(FlagLoadX0 && FlagUpdateVar &&  FlagMakeDxZero && fabs(Dx[i])>valp) {
          valp=fabs(Dx[i]*param0/x0[i]);
          Ptr=ACptr;
        }
      } else if (ACptr->Kbg<1) {
        ACptrp=ACptr->ContBus->AC;
        if (FlagLoadX0) {
          x0[i]=ACptr->V;
          if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
        }
        if (FlagUpdateVar) {
          ACptr->V=x0[i]+Dx[i];
          if (ACptr->V<=0){
            count++; val=fabs((ACptr->V-0.00001)/Dx[i]);
            if(ExistParameter('d')) fprintf(stderr,"%s %d %s Vzero %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
            ACptr->V=0.00001;}
          else if (strpbrk(ACptrp->cont,"Q")) {
            if (Recover && ACptrp->Qg>=ACptrp->Max && ACptr->V>=ACptr->VCont){
              count++; val=fabs((ACptr->V-ACptr->VCont)/Dx[i]);
              if(ExistParameter('d')) fprintf(stderr,"%s %d %s Recover%s %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,Qmax,val);
              ACptr->V=ACptr->VCont;}
            else if (Recover && ACptrp->Qg<=ACptrp->Min && ACptr->V<=ACptr->VCont){
              count++; val=fabs((ACptr->V-ACptr->VCont)/Dx[i]);
              if(ExistParameter('d')) fprintf(stderr,"%s %d %s Recover%s %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,Qmin,val);
              ACptr->V=ACptr->VCont;}
            else val=0;
          }
          else if (strpbrk(ACptrp->cont,"E")) {
            if (Recover && ACptrp->Gen->Eq>=ACptrp->Gen->EqMax && ACptr->V>=ACptr->VCont){
              count++; val=fabs((ACptr->V-ACptr->VCont)/Dx[i]);
              if(ExistParameter('d')) fprintf(stderr,"%s %d %s RecoverEqMax %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
              ACptr->V=ACptr->VCont;}
            else if (Recover && ACptrp->Gen->Eq<=ACptrp->Gen->EqMin && ACptr->V<=ACptr->VCont){
              count++; val=fabs((ACptr->V-ACptr->VCont)/Dx[i]);
              if(ExistParameter('d')) fprintf(stderr,"%s %d %s RecoverEqMin%lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
              ACptr->V=ACptr->VCont;}
            else val=0;
          }
          else if (strpbrk(ACptrp->cont,"I")) {
            if (Recover && ACptrp->Gen->Ia>=ACptrp->Gen->IaMax && ACptr->V>=ACptr->VCont){
              count++; val=fabs((ACptr->V-ACptr->VCont)/Dx[i]);
              if(ExistParameter('d')) fprintf(stderr,"%s %d %s RecoverIaMax %lf\n",ACptr->Type,ACptr->Num,ACptr->Name,val);
              ACptr->V=ACptr->VCont;}
            else val=0;
          }
          else val=0;
          if (val>consp) consp=val;
        }
      } else {
        if (FlagLoadX0) {
          x0[i]=ACptr->Qr;
          if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
        }
        if (FlagUpdateVar) ACptr->Qr=x0[i]+Dx[i];
      }
    }
    else if(strpbrk(ACptr->Type,"T")){
      i=ACvar[ACptr->N];
      if (FlagLoadX0) {
        x0[i]=ACptr->Ang;
        if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
      }
      if (FlagUpdateVar) ACptr->Ang=x0[i]+Dx[i];
      i++;
      if (Rcont) {
        for(ELptr=ACptr->Reg;ELptr!=NULL;ELptr=ELptr->Next){
          Eptr=ELptr->Eptr;
          if(!strcmp(Eptr->Type,"R")) {
            if (FlagLoadX0) {
              x0[i]=Eptr->Tap;
              if (FlagMakeDxZero && (DxZero!=NULL && (DxZero[i] || (flagReduceSystem && x0[i]!=0 && CountSteps>=NumSteps && fabs(Dx[i]/x0[i])<ZeroDx)))) { Dx[i]=0;  DxZero[i]=TRUE; NewNumEq--; }
            }
            if (FlagUpdateVar) {
              Eptr->Tap=x0[i]+Dx[i];
              if (Tlim && Eptr->Tap<=1/Eptr->Tmax+0.00001){
                count++; val=fabs((Eptr->Tap-1/Eptr->Tmax)/Dx[i]);
                if(ExistParameter('d')) fprintf(stderr,"%s %d %s %d %s Tmax %lf\n",Eptr->Type,Eptr->From->Num,Eptr->From->Name,