homot.c

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

/* Homotopy Continuation Method: Main. */

#include "homot.h"

#ifdef ANSIPROTO
void AddLoad(void);
BOOLEAN CheckVIlimits(BOOLEAN flagVoltage,BOOLEAN flagCurrent);
#else
void AddLoad();
BOOLEAN CheckVIlimits();
#endif


/* ------- Global Variables ------ */
VALUETYPE *Dx,Dparam,param0,*x0,*x0p,Kh,Htol,AngTr,VoltageSum,VoltageSum0,
          DxiMax,VSFone,VSFinf,SF,TVI,ZeroDx,lambda_o,TotalPl,TotalQl,TotalPg,TotalQg;
INDEX RedNumEq,NewNumEq,CountSteps,NumSteps,TVIbus,TVIrank,VSFbus;
BOOLEAN flagReduceSystem,*DxZero,flagPrintTotalPl,flagPrintTotalQl,flagPrintTotalPg,flagPrintTotalQg;
AClist *Vlist=NULL,*Vlistp=NULL;
FILE *OutputHomot;
int field,SD0;
extern BOOLEAN flagBS,flagVloads;


/* --------------------------- DCsetup --------------------------------- */
#ifdef ANSIPROTO
BOOLEAN DCsetup(void)
#else
BOOLEAN DCsetup()
#endif
/* Setup inital DC control mode, i.e., rectifier -> alfa,Id (or P),
                                       inverter  -> gamma,Vd. */
{
  DCbusData *DCptrR,*DCptrI;
  BOOLEAN flag=FALSE,flagP=FALSE,flagT=FALSE;

  for (DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next) if(!strcmp(DCptrR->Type,"R")) {
    DCptrI=DCptrR->To;
    DCptrI->VdN=DCptrI->Vd;
    DCptrR->AlfaN=DCptrR->Alfa;
    if (strcmp(DCptrR->Cont1,"AL") && strcmp(DCptrR->Cont2,"AL")) flag=TRUE;
    else if (strcmp(DCptrR->Cont1,"AT") && strcmp(DCptrR->Cont2,"AT")) flag=TRUE;
    else if (!strpbrk(DCptrR->Cont1,"IP") && !strpbrk(DCptrR->Cont2,"IP")) flag=TRUE;
    else if (strcmp(DCptrI->Cont1,"GA") && strcmp(DCptrI->Cont2,"GA")) flag=TRUE;
         else if (strcmp(DCptrI->Cont1,"AT") && strcmp(DCptrI->Cont2,"AT")) flag=TRUE;
    else if (strcmp(DCptrI->Cont1,"VD") && strcmp(DCptrI->Cont2,"VD")) flag=TRUE;
    if (!strcmp(DCptrR->Cont1,"PA") ||!strcmp(DCptrR->Cont2,"PA") ||
        !strcmp(DCptrI->Cont1,"PA") ||!strcmp(DCptrI->Cont2,"PA")) flagP=TRUE;
    if (!strcmp(DCptrR->Cont1,"AL") ||!strcmp(DCptrR->Cont2,"AL")) strcpy(DCptrR->Cont1,"AL");
    else if(!strcmp(DCptrR->Cont1,"AT") ||!strcmp(DCptrR->Cont2,"AT")) {
      strcpy(DCptrR->Cont1,"AT");
      flagT=TRUE;
    }
    else strcpy(DCptrR->Cont1,"AL");
    if (flagP) strcpy(DCptrR->Cont2,"PA");
    else strcpy(DCptrR->Cont2,"ID");
    strcpy(DCptrI->Cont1,"GA");
    if (flagT) strcpy(DCptrI->Cont2,"AT");
    else strcpy(DCptrI->Cont2,"VD");
    if (DCptrR->Tap>DCptrR->TapMax) {
      flag=TRUE;
      DCptrR->TapMax=DCptrR->Tap;
      fCustomPrint(stderr,"***Warning: The program will change the maximum tap limit for rect. %s\n",DCptrR->Name);
    }
    if (DCptrR->Tap<DCptrR->TapMin) {
      flag=TRUE;
      DCptrR->TapMin=DCptrR->Tap;
      fCustomPrint(stderr,"***Warning: The program will change the minimum tap limit for rect. %s\n",DCptrR->Name);
    }
    if (DCptrI->Tap>DCptrI->TapMax) {
      flag=TRUE;
      DCptrI->TapMax=DCptrI->Tap;
      fCustomPrint(stderr,"***Warning: The program will change the maximum tap limit for inv. %s\n",DCptrI->Name);
    }
    if (DCptrI->Tap<DCptrI->TapMin) {
      flag=TRUE;
      DCptrI->TapMin=DCptrI->Tap;
      fCustomPrint(stderr,"***Warning: The program will change the minimum tap limit for inv. %s\n",DCptrI->Name);
    }
    if (DCptrI->Gamma!=DCptrI->GammaMin) {
      DCptrI->GammaMin=DCptrI->Gamma;
      fCustomPrint(stderr,"***Warning: The program will change the minimum gamma limit for inv. %s\n",DCptrI->Name);
    }
    fCustomPrint(stderr,"***Warning: The startup control mode for HVDC link from %s to %s\n",DCptrR->Name,DCptrI->Name);
    fCustomPrint(stderr,"            is: rectifier -> %s-%s\n",DCptrR->Cont1,DCptrR->Cont2);
    fCustomPrint(stderr,"                inverter  -> %s-%s\n",DCptrI->Cont1,DCptrI->Cont2);
  }
  return(flag);
}

/* -------------------- ChangeDCmode ----------------------- */
#ifdef ANSIPROTO
BOOLEAN ChangeDCmode(void)
#else
BOOLEAN ChangeDCmode()
#endif
/* Change DC control mode when DC variable at a limit. */
{
  DCbusData *DCptrR,*DCptrI;
  BOOLEAN flag=FALSE,flagp=FALSE;
  INDEX i;

  i=NacVar;
  for (DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next) if(!strcmp(DCptrR->Type,"R")) {
    DCptrI=DCptrR->To;
    if(!strcmp(DCptrR->Cont1,"AL") && strpbrk(DCptrR->Cont2,"IP") &&
       (DCptrR->Tap>=DCptrR->TapMax || DCptrR->Tap<=DCptrR->TapMin)) {
      flag=flagp=TRUE;
      strcpy(DCptrR->Cont1,"AT");
    }
    else if(!strcmp(DCptrR->Cont1,"AT") && strpbrk(DCptrR->Cont2,"IP") &&
            DCptrR->Alfa<=DCptrR->AlfaMin) {
                flag=flagp=TRUE;
      if (!strcmp(DCptrR->Cont2,"ID")) strcpy(DCptrI->Cont2,"ID");
      else strcpy(DCptrI->Cont2,"PA");
      strcpy(DCptrR->Cont1,"AL");
      strcpy(DCptrR->Cont2,"AT");
      strcpy(DCptrI->Cont1,"AT");
    }
    else if(!strcmp(DCptrR->Cont1,"AT") && strpbrk(DCptrR->Cont2,"IP") &&
            DCptrR->Alfa==DCptrR->AlfaN) {
      flag=flagp=TRUE;
      strcpy(DCptrR->Cont1,"AL");
    }
    else if(!strcmp(DCptrI->Cont1,"AT") && strpbrk(DCptrI->Cont2,"IP") &&
            DCptrI->Tap<=DCptrI->TapMin && DCptrI->Gamma<=DCptrI->GammaMin) {
      flag=flagp=TRUE;
      if (!strcmp(DCptrI->Cont2,"ID")) strcpy(DCptrR->Cont2,"ID");
      else strcpy(DCptrR->Cont2,"PA");
      strcpy(DCptrR->Cont1,"AT");
      strcpy(DCptrI->Cont1,"GA");
      strcpy(DCptrI->Cont2,"AT");
    }
    else if(!strcmp(DCptrI->Cont1,"AT") && strpbrk(DCptrI->Cont2,"IP") &&
                 DCptrI->Tap>DCptrI->TapMin && DCptrI->Gamma<=DCptrI->GammaMin) {
      flag=flagp=TRUE;
      if (!strcmp(DCptrI->Cont2,"ID")) strcpy(DCptrR->Cont2,"ID");
      else strcpy(DCptrR->Cont2,"PA");
      strcpy(DCptrR->Cont1,"AT");
      strcpy(DCptrI->Cont1,"GA");
      strcpy(DCptrI->Cont2,"VD");
    }
    else if(!strcmp(DCptrI->Cont1,"GA") && !strcmp(DCptrI->Cont2,"VD") &&
            (DCptrI->Tap<=DCptrI->TapMin || DCptrI->Tap>=DCptrI->TapMax)) {
      flag=flagp=TRUE;
      strcpy(DCptrI->Cont2,"AT");
    }
    else if(!strcmp(DCptrI->Cont1,"GA") && !strcmp(DCptrI->Cont2,"AT") &&
            DCptrI->Vd==DCptrI->VdN) {
      flag=flagp=TRUE;
      strcpy(DCptrI->Cont2,"VD");
    }
    if (flagH && flag) {
      if(strcmp(DCptrR->Cont1,"VD")&&strcmp(DCptrR->Cont2,"VD")) {
        i++; x0[i]=DCptrR->Vd;
      }
      if(strcmp(DCptrR->Cont1,"AT")&&strcmp(DCptrR->Cont2,"AT")) {
        i++; x0[i]=DCptrR->Tap*DCptrR->Ntrf;
      }
      if(strcmp(DCptrR->Cont1,"AL")&&strcmp(DCptrR->Cont2,"AL")) {
        i++; x0[i]=cos(DCptrR->Alfa);
      }
      if(strcmp(DCptrR->Cont1,"GA")&&strcmp(DCptrR->Cont2,"GA")) {
        i++; x0[i]=cos(DCptrR->Gamma);
       }
      i++; x0[i]=DCptrR->MVA;
      if(strcmp(DCptrR->Cont1,"PA")&&strcmp(DCptrR->Cont2,"PA")) {
        i++; x0[i]=DCptrR->P;
      }
      if(strcmp(DCptrR->Cont1,"QA")&&strcmp(DCptrR->Cont2,"QA")) {
        i++; x0[i]=DCptrR->Q;
      }
      if(strcmp(DCptrI->Cont1,"VD")&&strcmp(DCptrI->Cont2,"VD")) {
        i++; x0[i]=DCptrI->Vd;
      }
      if(strcmp(DCptrI->Cont1,"AT")&&strcmp(DCptrI->Cont2,"AT")) {
        i++; x0[i]=DCptrI->Tap*DCptrI->Ntrf;
      }
      if(strcmp(DCptrI->Cont1,"AL")&&strcmp(DCptrI->Cont2,"AL")) {
        i++; x0[i]=cos(DCptrI->Alfa);
      }
      if(strcmp(DCptrI->Cont1,"GA")&&strcmp(DCptrI->Cont2,"GA")) {
        i++; x0[i]=cos(DCptrI->Gamma);
       }
      i++; x0[i]=DCptrI->MVA;
      if(strcmp(DCptrI->Cont1,"PA")&&strcmp(DCptrI->Cont2,"PA")) {
        i++; x0[i]=DCptrI->P;
      }
      if(strcmp(DCptrI->Cont1,"QA")&&strcmp(DCptrI->Cont2,"QA")) {
        i++; x0[i]=DCptrI->Q;
      }
      i++; x0[i]=DCptrI->MVA;
      if(strcmp(DCptrR->Cont1,"ID")&&strcmp(DCptrR->Cont2,"ID")&&
         strcmp(DCptrI->Cont1,"ID")&&strcmp(DCptrI->Cont2,"ID")) {
        i++; x0[i]=DCptrR->Id;
      }
    }
    if (flagp) {
      flagp=FALSE;
      fCustomPrint(stderr,"***Warning: HVDC link from %s to %s will switch control mode\n",DCptrR->Name,DCptrI->Name);
      fCustomPrint(stderr,"            to: rectifier -> %s-%s\n",DCptrR->Cont1,DCptrR->Cont2);
      fCustomPrint(stderr,"                inverter  -> %s-%s\n",DCptrI->Cont1,DCptrI->Cont2);
    }
  }
  return(flag);
}

/* ----------------- Direction ---------------------------- */
#ifdef ANSIPROTO
int Direction(SparseMatrix *Mptr,VALUETYPE *vec,BOOLEAN flag)
#else
int Direction(Mptr,vec,flag)
SparseMatrix *Mptr;
VALUETYPE *vec;
BOOLEAN flag;
#endif
/* Find derivative -> dx/dlambda. */
{
  SparseMatrixElement *Jptr;
  AreaData *Aptr;
  ACbusData *ACptr;
  INDEX i,N;
  int m=0,PgMax,PgMaxH;

  N=Mptr->n1;
  if(!flag) {
    if(ExistParameter('d')) fCustomPrint(stderr,"Direction Factor Jacobian.\n");
    for(i=1;i<=N;i++) for(vec[i]=dx[i]=0,Jptr=Mptr->RowHead[i];Jptr!=NULL;Jptr->Value=0,Jptr=Jptr->RowNext);
    Aptr=ACFunJac(Mptr,&PgMax,FALSE,TRUE,FALSE);
    if(DCFunJac(Mptr,FALSE,TRUE)) return(-1);
    SVCFunJac(Mptr,FALSE,TRUE);                    /* FACTS */
    TCSCFunJac(Mptr,FALSE,TRUE);                   /* FACTS */
    STATCOMFunJac(Mptr,FALSE,TRUE);                /* FACTS */
    if (flagH) {
      PgMaxH=HFunJac(FALSE,TRUE,Aptr,dx);
      if(NewCol->p[N]==0) Jptr=Mptr->ColHead[N];
      else Jptr=Mptr->ColHead[NewCol->p[N]];
      while(Jptr!=NULL){
        i=OldRow->p[Jptr->Row];
        if(i==0) i=Jptr->Row;
        if (i<N) vec[i]= -Jptr->Value;
        Jptr=Jptr->ColNext;
      }
    }
    if (PgMax<0 && (!flagH || (flagH && PgMaxH<0)) ) {
      if (Aptr!=NULL) {
        fCustomPrint(stderr,"\nError: Area %d %s does not have any spinning reserves.\n",Aptr->N,Aptr->Name);
        fCustomPrint(stderr,"       Increase the maximum P generation in this area, otherwise\n");
      } else {
        fCustomPrint(stderr,"\nError: The system does not have any spinning reserves.\n");
        fCustomPrint(stderr,"       Increase the maximum P generation in this system, otherwise\n");
      }
      fCustomPrint(stderr,"       the Jacobian matrix becomes singular.\n");
      fCustomPrint(stderr,"Loading Factor -> %-6.3lf\n",lambda);
      WriteSolution(0,TrueParamStr(2),"Pg Max. Problems:");
      stopExecute(1);
    }
    if(NewCol->p[N]!=0) m=factor(Mptr);
    else m=WARNINGEXIT;
  }
  if (m==WARNINGEXIT || flag) {
    if(ExistParameter('d')) fCustomPrint(stderr,"Direction Order and Factor Jacobian.\n");
    DeleteJac(Mptr,NewRow,NewCol,OldRow,OldCol);
    Aptr=ACFunJac(Mptr,&PgMax,FALSE,TRUE,FALSE);
    if(DCFunJac(Mptr,FALSE,TRUE)) return(-1);
    SVCFunJac(Mptr,FALSE,TRUE);                    /* FACTS */
    TCSCFunJac(Mptr,FALSE,TRUE);                   /* FACTS */
    STATCOMFunJac(Mptr,FALSE,TRUE);                /* FACTS */
    if (flagH) {
      for (i=1; i<=N; vec[i]=dx[i]=0, i++);
      PgMaxH=HFunJac(FALSE,TRUE,Aptr,dx);
      for(Jptr=Mptr->ColHead[N];Jptr!=NULL;Jptr=Jptr->ColNext) {
        i=Jptr->Row;
        if (i<N) vec[i]= -Jptr->Value;
      }
    }
    if (PgMax<0 && (!flagH || (flagH && PgMaxH<0)) ) {
      if (Aptr!=NULL) {
        fCustomPrint(stderr,"\nError: Area %d %s does not have any spinning reserves.\n",Aptr->N,Aptr->Name);
        fCustomPrint(stderr,"       Increase the maximum P generation in this area, otherwise\n");
      } else {
        fCustomPrint(stderr,"\nError: The system does not have any spinning reserves.\n");
        fCustomPrint(stderr,"       Increase the maximum P generation in this system, otherwise\n");
      }
      fCustomPrint(stderr,"       the Jacobian matrix becomes singular.\n");
      fCustomPrint(stderr,"Loading Factor -> %-6.3lf\n",lambda);
      WriteSolution(0,TrueParamStr(2),"Pg Max. Problems:");
      stopExecute(1);
    }
    SortRowsColumns(Mptr);
    if(factorns(Mptr,alpha,RowPartition,ColPartition,NewRow,NewCol,OldRow,OldCol)){
      fCustomPrint(stderr,"*** Singular Jacobian (possible voltage collapse, contol or limit problems).\n");
      fCustomPrint(stderr,"    Try changing the load levels, controls or limits, or use the -F option.\n");
      fCustomPrint(stderr,"Loading Factor -> %-6.3lf\n",lambda);
      WriteSolution(0,TrueParamStr(2),"Singular Jacobian:");
      stopExecute(1);
    }
    SortRowsColumns(Mptr);
  }
  if (!flagH) {
    for (i=1; i<=N; vec[i]=0, i++);
    for (ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next) {
      i=ACvar[ACptr->N];
      vec[i]= ACptr->Pnl*pow(ACptr->V,ACptr->a)+ACptr->Pzl*ACptr->V*ACptr->V;
      vec[i+1]=ACptr->Qnl*pow(ACptr->V,ACptr->b)+ACptr->Qzl*ACptr->V*ACptr->V;
    }
  }
  repsolp(Mptr,vec,OldRow,NewCol);
  if (m==WARNINGEXIT) SD0=0;
  return(0);
}

/* --------------------------- ChangeParam --------------------------------- */
#ifdef ANSIPROTO
BOOLEAN ChangeParam(void)
#else
BOOLEAN ChangeParam()
#endif
/* Change parameter for Homotopy method. */
{
  INDEX i;

  if (DxiMax>1) {
    Dparam=0;
    for (i=1;i<Jac->n1;i++) Dx[i]=0;
    LoadX0(FALSE,TRUE,TRUE);
    if (BlPtr->N==Bl) {
      strcpy(BlPtr->Type,"B");
      if(BlPtr->Area!=NULL && BlPtr->Area->Slack==BlPtr) strcat(BlPtr->Type,"A");
      Bl=0;
      BlPtr->Cont=BlPtr;
      if(ExistParameter('d')) fCustomPrint(stderr,"Change param. to lambda.\n");
      return(TRUE);
    } else {
      strcpy(BlPtr->Type,"BL");
      if(BlPtr->Area!=NULL && BlPtr->Area->Slack==BlPtr) strcat(BlPtr->Type,"A");
      Bl=BlPtr->N;
      BlPtr->Cont=NULL;
      if(ExistParameter('d')) fCustomPrint(stderr,"Change param. to V at %s %d %s.\n",BlPtr->Type,BlPtr->Num,BlPtr->Name);
      return(TRUE);
    }
  } else return(FALSE);
}

/* ------------------------- VoltageSensFactor --------------------------- */
#ifdef ANSIPROTO
void VoltageSensFactor(VALUETYPE *vec,BOOLEAN first)
#else
void VoltageSensFactor(vec,first)
VALUETYPE *vec;
BOOLEAN first;
#endif
/* Calculate Voltage Sensitivity Factor and
   Tangent Vector Index and rank for a bus */
{
  ACbusData *ACptr;
  VALUETYPE val;
  INDEX i,l,N=0,MaxBus=0;
  ACranklist  *RankList,*Ptr,*PtrMax,*PrevPtr,*PrevPtrMax;


  if (first) TVIbus=IntegerParameter('f',0,1,9999);
  RankList=NULL;
  for(VSFinf=VSFone=TVI=0,ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next){
    if (ACptr->Cont!=NULL){
      i=ACvar[ACptr->N];
      val=fabs(vec[i+1]);
      Ptr=RankList;
#ifdef WINDOWS
      RankList= new ACranklist;
#else
      RankList=(ACranklist *) malloc(sizeof(ACranklist));
      if (RankList==NULL) {ErrorHalt("Insufficient memory to allocate rank data (option -fnum)."); stopExecute(ERROREXIT);}
#endif
      RankList->AC=ACptr;
      RankList->val=val;
      RankList->Next=Ptr;
      N++;
      if (val>VSFone) {
        VSFone=val;
        VSFbus=MaxBus=ACptr->Num;
      }
      VSFinf=VSFinf+val;
      if (TVIbus==ACptr->Num) TVI=val;