homot.c

电力系统分析计算 学习调试程序 UNIX / LINUX / CYGWIN 系统使用

    }
  }
  if (first && TVI==0) {
    TVIbus=MaxBus;
    TVI=VSFone;
  }
  /*    Rank buses based on Tangent Vector   */
  if (TVI!=0) for(l=1;l<=N;l++) {
    if (RankList==NULL || RankList->Next==NULL) break;
    for(val=-0.1,PrevPtr=Ptr=RankList;Ptr!=NULL;PrevPtr=Ptr,Ptr=Ptr->Next) {
      if (Ptr->val>val) {
        val=Ptr->val;
        PtrMax=Ptr;
        PrevPtrMax=PrevPtr;
      }
    }
    if(PtrMax!=NULL) {
      if (PrevPtrMax!=PtrMax) PrevPtrMax->Next=PtrMax->Next;
      else RankList=PtrMax->Next;
      if (ExistParameter('d')) {
        fCustomPrint(stderr,"Ranked List: ");
        fCustomPrint(stderr,"%4d %12s %-11.5g\n",PtrMax->AC->Num,PtrMax->AC->Name,PtrMax->val);
      }
      if (TVIbus==PtrMax->AC->Num) {TVIrank=l; break;}
#ifdef WINDOWS
      delete PtrMax;
#else
      free(PtrMax);
#endif
    }
  }
  for(Ptr=RankList;Ptr!=NULL;){
    PrevPtr=Ptr->Next;
#ifdef WINDOWS
    delete Ptr;
#else
    free(Ptr);
#endif
    Ptr=PrevPtr;
  }
}


/* --------------------------  CheckVIlimits ------------------------ */
#ifdef ANSIPROTO
BOOLEAN CheckVIlimits(BOOLEAN flagVoltage,BOOLEAN flagCurrent)
#else
BOOLEAN CheckVIlimits(flagVoltage,flagCurrent)
BOOLEAN flagVoltage,flagCurrent;
#endif
{
  ACbusData *ACptr;
  ElementData *Eptr;
  ElementList *ELptr;
  VALUETYPE Pij,Qij,Vi,di,Vj,dj,Iij,G,B,Gi,Bi,Pl,Ql,Gp,Bp,Gj,Bj,Pji,Qji,Iji;
  BOOLEAN FlagStopContinuation=FALSE;

  TotalPl=TotalQl=TotalPg=TotalQg=0;
  for (ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next) {
    TotalPg=TotalPg+ACptr->PG;
    Pl=(ACptr->Pn+lambda*ACptr->Pnl)*pow(ACptr->V,ACptr->a)+
       (ACptr->Pz+lambda*ACptr->Pzl)*ACptr->V*ACptr->V;
    ACptr->PL=Pl;
    TotalPl=TotalPl+Pl;
    Ql=(ACptr->Qn+lambda*ACptr->Qnl)*pow(ACptr->V,ACptr->b)+
       (ACptr->Qz+lambda*ACptr->Qzl)*ACptr->V*ACptr->V;
    ACptr->QL=Ql;
    TotalQl=TotalQl+Ql;
    TotalQg=TotalQg+ACptr->Qg;
    if(ACptr->CheckVlimits && ACptr->Vlmax>ACptr->Vlmin && flagVoltage) {
      if(ACptr->V<=ACptr->Vlmin) {
        fCustomPrint(stderr,"***Warning: Bus %d %s has reached or exceeded its minimum V limit.\n",ACptr->Num,ACptr->Name);
        ACptr->CheckVlimits=FALSE;
        FlagStopContinuation=TRUE;
      }
      else if(ACptr->V>=ACptr->Vlmax) {
        fCustomPrint(stderr,"***Warning: Bus %d %s has reached or exceeded its maximum V limit.\n",ACptr->Num,ACptr->Name);
        ACptr->CheckVlimits=FALSE;
        FlagStopContinuation=TRUE;
      }
    }
    if (flagCurrent) for (ELptr=ACptr->Elem;ELptr!=NULL;ELptr=ELptr->Next) {
      Eptr=ELptr->Eptr;
      if (Eptr->CheckIlimits && Eptr->Imax>0)  {
        Vi=Eptr->From->V;  di=Eptr->From->Ang;
        Vj=Eptr->To->V;    dj=Eptr->To->Ang;
        G=(Eptr->G*cos(Eptr->Ang)-Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        B=(Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
        Gi=(Eptr->G1+Eptr->G)*pow(Eptr->Tap,2.0)-G;
        Bi=(Eptr->B1+Eptr->B)*pow(Eptr->Tap,2.0)-B;
        Gp=(Eptr->G*cos(Eptr->Ang)+Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        Bp=(-Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
        Gj=Eptr->G+Eptr->G2-Gp;
        Bj=Eptr->B+Eptr->B2-Bp;
        Pij=Vi*Vi*(Gi+G)-Vi*Vj*(G*cos(di-dj)+B*sin(di-dj));
        Qij= -Vi*Vi*(Bi+B)-Vi*Vj*(G*sin(di-dj)-B*cos(di-dj));
        Iij=sqrt(Pij*Pij+Qij*Qij)/Vi;
        Pji=Vj*Vj*(Gj+Gp)-Vi*Vj*(Gp*cos(dj-di)+Bp*sin(dj-di));
        Qji= -Vj*Vj*(Bj+Bp)-Vi*Vj*(Gp*sin(dj-di)-Bp*cos(dj-di));
        Iji=sqrt(Pji*Pji+Qji*Qji)/Vj;
        if(Iij>=Eptr->Imax) {
          fCustomPrint(stderr,"***Warning: Element %d %s to %d %s has reached\n",Eptr->From->Num,Eptr->From->Name,Eptr->To->Num,Eptr->To->Name);
          fCustomPrint(stderr,"            or exceeded its maximum I limit at its first bus.\n");
          Eptr->CheckIlimits=FALSE;
          FlagStopContinuation=TRUE;
        }
        else if(Iji>=Eptr->Imax) {
          fCustomPrint(stderr,"***Warning: Element %d %s to %d %s has reached\n",Eptr->From->Num,Eptr->From->Name,Eptr->To->Num,Eptr->To->Name);
          fCustomPrint(stderr,"            or exceeded its maximum I limit at its second bus.\n");
          Eptr->CheckIlimits=FALSE;
          FlagStopContinuation=TRUE;
        }
      }
    }
  }
  return(FlagStopContinuation);
}


/* ------------------------- AddLoad --------------------------- */
#ifdef ANSIPROTO
void AddLoad(void)
#else
void AddLoad()
#endif
/* Add initial loading to load buses  */
{
  ACbusData *ACptr;

  for(ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next){
     ACptr->Pn=ACptr->Pn+lambda*ACptr->Pnl;
     ACptr->Pz=ACptr->Pz+lambda*ACptr->Pzl;
     ACptr->Qn=ACptr->Qn+lambda*ACptr->Qnl;
     ACptr->Qz=ACptr->Qz+lambda*ACptr->Qzl;
  }

}

/* --------------------------- Homot --------------------------------- */
#ifdef ANSIPROTO
int Homot(void)
#else
int Homot()
#endif
/* Build voltage profiles. */
{
  VALUETYPE NDx,cons,PreviousLambda,StoplambdaVal,FirstVoltageSumUp=Tol;
  int iter,iterp,CountCyclingSteps,CountVupSteps,NumHVsolutions,
      Steps,IndicesSteps,MaxSteps,i,N,n,CountStepsAfterVup=0;
  BOOLEAN flagLimits=FALSE,flagSTOP=FALSE,flagFACTSlimits=FALSE,flagt=FALSE,stopValueSet=FALSE,
          flagMakeVlist=TRUE,flagCountVupAfterFirstLimit=FALSE,flagD=FALSE,
          ForceflagCountVupAfterFirstLimit=TRUE,flagFirstStep=TRUE,
          flagCycling=FALSE,flagRlimits=FALSE,flagQlimits=FALSE,
          flagVlimits=FALSE,flagVIlimits=FALSE,flagDClimits=FALSE,
          flagVoltage=FALSE,flagCurrent=FALSE,
          flagSVClimits=FALSE,flagTCSClimits=FALSE,flagSTATCOMlimits=FALSE;  /* FACTS */

  BlPtr=NULL;
  Dparam=param0=0;
  AngTr=180/PI;
  Kh=1;
  RealParameter('k',&Kh,1e-4,1e4);
  Htol=0.000001;
  RealParameter('o',&Htol,1e-10,1e-1);
  if (ExistParameter('H')) OutputHomot=(FILE *) OpenOutput(NameParameter('H'));
  else OutputHomot=(FILE *) OpenOutput(NameParameter('c'));
  N=Jac->n1;
  RedNumEq=N-1;
  CountSteps=0;
  NumSteps=IntegerParameter('U',10,2,100);
  iter=Pflow(1,FALSE,TRUE,TRUE);
  if(iter<0) { fclose(OutputHomot); return(iter);}
  else fCustomPrint(stderr,"**** Initial Power Flow Solved ****    \n\n");
  lambda=SD0=0;
  RealParameter('L',&lambda,-1e6,1e6);
  lambda_o=lambda;
  NumHVsolutions=IntegerParameter('2',5,2,100);
  flagVoltage=ExistParameter('7');
  flagCurrent=ExistParameter('8');
  if (lambda!=0) {
    if ((ExistParameter('D') && (!NullName(NameParameter('D')))) || flagVloads) flagD=TRUE;
    iter=Pflow(iter,flagD,TRUE,FALSE);
    if (iter>0) {
      fCustomPrint(stderr,"**** Initial Lambda Case Solved ****    ");
      fCustomPrint(stderr,"Starting loading factor -> %-10.6lg\n\n",lambda);
      AddLoad();
    } else {fclose(OutputHomot); return(iter);}
    lambda=0;
  }
  flagt=DCsetup();
  field=IntegerParameter('O',6,6,10);
  MaxSteps=IntegerParameter('z',1000,1,9999);
  StoplambdaVal=lambda;
  if (ExistParameter('u')) {
    ZeroDx=0.001;
    flagReduceSystem=TRUE;
#ifdef WINDOWS
    DxZero= new BOOLEAN[N];
#else
    DxZero=(BOOLEAN *) calloc(N,sizeof(BOOLEAN));
    if (DxZero==NULL) {
      fclose(OutputHomot);
      ErrorHalt("Insufficient memory to allocate vector for system reduction in the Cont. Method.");
      stopExecute(ERROREXIT);
    }
#endif
    RealParameter('u',&ZeroDx,0.,0.2);
    if (Kh<1) ZeroDx=Kh*ZeroDx;
  } else { flagReduceSystem=FALSE; ZeroDx=0;  DxZero=NULL;}
  for (i=1; i<=N; i++) { if (DxZero!=NULL) DxZero[i]=FALSE;}
  for (flagL=flagR=TRUE,CountCyclingSteps=CountVupSteps=Steps=1;;) {
    Homotopy:
    if(ExistParameter('d')) fCustomPrint(stderr,"\nContinuation Step=%d\n",Steps);
    Dparam=0;
    PreviousLambda=lambda;
    /* Predictor Step */
    if(Direction(Jac,Dx,flagt)<0) {fclose(OutputHomot); return(-iter);}
    if (SD0==0) SD0=DetSign;
    for (NDx=SF=0,i=1;i<=N-1;i++) {
      if(fabs(Dx[i])>NDx) NDx=fabs(Dx[i]);
      SF=SF+fabs(Dx[i]);
    }
    if (flagFirstStep) {
      if (ExistParameter('u') && ZeroDx>0) flagReducedContinuation=TRUE;
      flagFirstStep=FALSE;
    }
    if (ExistParameter('f')) VoltageSensFactor(Dx,flagMakeVlist);
    if (flagMakeVlist) {
      MakeVlist(OutputHomot);
      CheckVIlimits(flagVoltage,flagCurrent);
      VoltProf(TRUE,OutputHomot);
      if (ExistParameter('0')) {IndicesSteps=Steps; WriteQmatrix(IndicesSteps,Dx);}
      flagMakeVlist=FALSE;             
    }
    Steps++;
    if (ExistParameter('f')) {
      if (Kh<0) {VSFone=VSFinf=SF=TVI=0; VSFbus=TVIrank=0;}
      Print(OutputHomot,0,6,4,VSFone);
      fCustomPrint(OutputHomot,"    %6d    ",VSFbus);
      Print(OutputHomot,0,6,4,VSFinf);
      fCustomPrint(OutputHomot,"    ");
      Print(OutputHomot,0,6,4,SF);
      fCustomPrint(OutputHomot,"    ");
      if (TVI!=0) Print(OutputHomot,0,6,4,1./TVI);
      else fCustomPrint(OutputHomot,"%6d",0);
      fCustomPrint(OutputHomot,"    %6d",TVIrank);
    }
    fCustomPrint(OutputHomot,"\n");
    if (NDx) {
      if(ExistParameter('d')) fCustomPrint(stderr,"Lambda->%lf  SD0=%2d  Det.Sign=%2d  Kh=%lf\n",lambda,SD0,DetSign,Kh);
      if (Bl) { cons= -1; Kh= -fabs(Kh); SD0=0;}
      else { if(SD0*DetSign<0) { SD0=DetSign; Kh= -Kh;} cons=Kh;}
      if(ExistParameter('d')) fCustomPrint(stderr,"%38s Kh=%lf\n","",Kh);
      Dparam=cons/NDx;
      for (i=1;i<=N-1;i++) Dx[i]=Dparam*Dx[i];
      if (flagSTOP) {
         if (ExistParameter('Z') && !NullName(NameParameter('Z'))) PrintDirection('Z',Dx,1.);
         break;
      }

      if (flagReducedContinuation) {
        /* Stop system reduction when V goes up, and start again
           after V goes down for 5 steps. */
        CountSteps++;
        if (VoltageSum>FirstVoltageSumUp  && Kh<0) flagReduceSystem=FALSE;
        if (!flagReduceSystem) {
          if (VoltageSum<=FirstVoltageSumUp && Kh<0) CountStepsAfterVup++;
          if (CountStepsAfterVup>=5) flagReduceSystem=TRUE;
        }
      }

      /* Update variables from predictor */
      cons=LoadX0(TRUE,TRUE,TRUE);
      if (!CountVupSteps) FirstVoltageSumUp=VoltageSum0;
      if (flagReducedContinuation && flagReduceSystem && CountSteps>=NumSteps) {
        CountSteps=0;
        if (NewNumEq<RedNumEq) { RedNumEq=NewNumEq; SD0=0; }
      }
      if(ExistParameter('d')) fCustomPrint(stderr,"cons=%lf Dparam=%lf\n",cons,Dparam);
      if(ExistParameter('d')) fCustomPrint(stderr,"Num.Steps=%d  Num.Eqs.System=%d   Num.Eqs.Reduced=%d\n",CountSteps,N-1,RedNumEq);
      if(!flagt) { if (ExistParameter('H')) flagt=ChangeParam();}
      else flagt=FALSE;
      if (flagt) goto Homotopy;

      /*  Enforce Limits */
      if (cons>0) {
        Dparam=fabs(1-cons+Htol*0.1)*Dparam;
        if(ExistParameter('d')) fCustomPrint(stderr,"(1-%lf)Dparam->%lf\n",cons,Dparam);
        for (i=1;i<=N-1;i++) Dx[i]=fabs(1-cons+Htol*0.1)*Dx[i];
        LoadX0(FALSE,TRUE,TRUE);
        if ((!ExistParameter('G') || !flagCountVupAfterFirstLimit || !CountVupSteps) && fabs(Dparam)<=Htol) {
          flagRlimits=CheckRlimits();
          flagVlimits=CheckVlimits();
          flagQlimits=CheckQlimits();
          flagDClimits=ChangeDCmode();
          flagSVClimits=ChangeSVCmode();          /* FACTS */
          flagTCSClimits=ChangeTCSCmode();        /* FACTS */
          flagSTATCOMlimits=ChangeSTATCOMmode();  /* FACTS */
          if (flagRlimits || flagVlimits || flagQlimits || flagDClimits ||
              flagSVClimits || flagTCSClimits || flagSTATCOMlimits)  flagLimits=TRUE;  /* FACTS */
          else                                                       flagLimits=FALSE;
          if (flagDClimits || flagSVClimits || flagTCSClimits || flagSTATCOMlimits) flagFACTSlimits=TRUE;   /* FACTS */
          else                                                                      flagFACTSlimits=FALSE;  /* FACTS */
          if (flagLimits) {SD0=0; if(ForceflagCountVupAfterFirstLimit) flagCountVupAfterFirstLimit=TRUE;}
        } else flagLimits=flagFACTSlimits=FALSE;   /* FACTS */
      } else flagLimits=flagFACTSlimits=FALSE;     /* FACTS */

      /* Corrector step */
      for (n=1;n<=10;n++) {
        if (n==10) flagL=flagR=FALSE;
        iterp=Pflow(iter,flagFACTSlimits,flagLimits,FALSE);
        if (flagL == FALSE) flagL=flagR=TRUE;
        /* if ((iterp>0 && lambda>-1e-4) || fabs(Dparam)<=1e-6)  break;*/
        if (iterp>0 && lambda>-1e-4)  break;
        Dparam=Dparam/2;
        for (i=1;i<=N-1;i++) Dx[i]=Dx[i]/2;
        LoadX0(FALSE,TRUE,TRUE);
        if(ExistParameter('d')) fCustomPrint(stderr,"Dparam/2^%d->%lf\n",n,Dparam);
      }

      /* Turning point computations */
      if (flagCountVupAfterFirstLimit && VoltageSum>FirstVoltageSumUp && Kh<0) CountVupSteps++;
      else CountVupSteps=0;
      if(ExistParameter('d')) fCustomPrint(stderr,"V summation->%lf     First V sum Up->%lf    Vup Steps->%d\n",VoltageSum,FirstVoltageSumUp,CountVupSteps);
      if (ExistParameter('G') && flagCountVupAfterFirstLimit &&
          (CountVupSteps>=NumHVsolutions || (CountCyclingSteps>0 && iterp<0))) {
        Kh= -Kh;
        ForceflagCountVupAfterFirstLimit=flagCountVupAfterFirstLimit=FALSE;
        if (iterp<0) iterp= -iterp;
      }

      /* Stop and output commands */
	  // when lamda hits its max, set the stop value
      if(lambda<PreviousLambda && !stopValueSet) {
        stopValueSet=TRUE;
        RealParameter('S',&StoplambdaVal,0.,1.);
        StoplambdaVal = PreviousLambda*StoplambdaVal;
        if (ExistParameter('E') && !NullName(NameParameter('E'))) {
          if(ExistParameter('d')) fCustomPrint(stderr,"Write right e-vector at lambda->%lf\n",PreviousLambda);
          PrintDirection('E',Dx,1.0);
        }
      }
	  // Sometimes there is a small decrease in lamda before it actually hits its max
	  // the following will reset the stopValueSet flag incase that happened
	  if(stopValueSet && lambda>PreviousLambda){
		  stopValueSet = FALSE;
	  }

      if ((MaxSteps>0 && Steps>=MaxSteps) || (lambda-StoplambdaVal)<=1e-5 ) flagSTOP=TRUE;
      if (iterp<0)  {
        iter= -iterp;
        fCustomPrint(stderr,"\n *** The case diverges (possible AC/DC control problems).\n");
        fCustomPrint(stderr,"     Try reducing the step size (-k option).\n");
        fCustomPrint(stderr,"Loading Factor -> %-6.3lf\n",lambda);
        WriteSolution(--iter,TrueParamStr(2),"Divergence:");
        stopExecute(1);
      } else iter=iterp;
      flagVIlimits=CheckVIlimits(flagVoltage,flagCurrent);
      VoltProf(FALSE,OutputHomot);
      if (flagVIlimits) break;
      if (ExistParameter('0') && Kh>0) {IndicesSteps=Steps; WriteQmatrix(IndicesSteps,Dx);}
      if (fabs(lambda-param0)<=Htol) CountCyclingSteps++;
      else CountCyclingSteps=0;
      if (CountCyclingSteps>=5 && !CountVupSteps){
        if (flagCycling) {
          fCustomPrint(stderr,"\n *** Cycling due to possible AC/DC/FACTS control problems.  Try using the debug (-d)\n");
          fCustomPrint(stderr,"     option to detect the limit that is creating the problem.\n");
          fCustomPrint(stderr,"Loading Factor -> %-6.3lf\n",lambda);
          WriteSolution(iter,TrueParamStr(2),"Divergence:");
          stopExecute(1);
        } 
		else Kh=-Kh;
        flagCycling=TRUE;
      }
    } else break;
  }
  if (ExistParameter('m')) MatlabV(OutputHomot);
  if(OutputHomot!=NULL)
	fclose(OutputHomot);
  if (ExistParameter('O')) TEFMatlabFiles();
  if (ExistParameter('0')) IndicesMatlab(IndicesSteps);
  return(iter);
}