pflow.cpp

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

#ifdef WINDOWS
      delete Jptr;
#else
      free(Jptr);
#endif
      Jptr=Jptrp;
    }
    Mptr->ColHead[k]=Mptr->RowHead[k]=NULL;
    P1Row->p[k]=P1Col->p[k]=0;
    P2Row->p[k]=P2Col->p[k]=0;
  }
}

/* --------------------------- Pflow --------------------------------- */
#ifdef ANSIPROTO
int Pflow(int iter,BOOLEAN flagDCLimits,BOOLEAN flagLimits,BOOLEAN flagFirst)
#else
int Pflow(iter,flagDCLimits,flagLimits,flagFirst)
int iter;
BOOLEAN flagDCLimits,flagLimits,flagFirst;
#endif
/* Power Flow Routine. */
{
  SparseMatrixElement *Jptr;
  AreaData *Aptr;
  int i,m=0,PgMax,PgMaxH;
  INDEX j,k,N,N1,N2,N3,MaxIterRun;
  VALUETYPE MaxdFi1,cons,val;
  BOOLEAN flag=FALSE,flagp=FALSE,flags=FALSE;


  N3=iter;
  if (flagH) MaxIterRun=MaxIter+iter-1;
  else MaxIterRun=MaxIter;
  RowPer=NewRow; ColPer=NewCol;
  NewtonRapson:
  ACFunJac(Jac,&PgMax,TRUE,FALSE,flagFirst);
  if(DCFunJac(Jac,TRUE,FALSE)) return(-iter);
  SVCFunJac(Jac,TRUE,FALSE);                  /*  FACTS  */
  TCSCFunJac(Jac,TRUE,FALSE);                 /*  FACTS  */
  STATCOMFunJac(Jac,TRUE,FALSE);              /*  FACTS  */
  if (flagH) HFunJac(TRUE,FALSE,NULL,Dx);
  else if (iter!=1 && flagPoC) {
    ACFunHes(TRUE,FALSE);
    if(DCFunHes(TRUE,FALSE)) return(-iter);
    SVCFunHes(TRUE,FALSE);                   /* FACTS  */
    TCSCFunHes(TRUE,FALSE);                  /* FACTS  */
    STATCOMFunHes(TRUE,FALSE);               /* FACTS  */
  }
  MaxdFi=Norm(dF,Jac->n1,&N1);
  fprintf(stderr,"\nIteration: %2d  ",N3);
  PrintMismatch(MaxdFi,0,N1);
  for (i=N3;i<=MaxIterRun;i++){
    N=Jac->n1;
    /* if (ExistParameter('d')) fprintf(stderr,"i %d   MaxdFi %lf   Tol %lf\n",i,MaxdFi,Tol); */
    if ((MaxdFi>Tol) || flagDCLimits) {
      if ((i>1 || flagR || flagL) && !flagDCLimits && (flagLimits || (i==N3))) {
        if(ExistParameter('d')) fprintf(stderr,"Factor Jacobian.\n");
        for(k=1; k<=N; k++) for(Jptr=Jac->RowHead[k];Jptr!=NULL;Jptr->Value=0,Jptr=Jptr->RowNext);
        Aptr=ACFunJac(Jac,&PgMax,TRUE,TRUE,flagFirst);
        if(DCFunJac(Jac,TRUE,TRUE)) return(-iter);
        SVCFunJac(Jac,TRUE,TRUE);                    /* FACTS */
        TCSCFunJac(Jac,TRUE,TRUE);                   /* FACTS */
        STATCOMFunJac(Jac,TRUE,TRUE);                /* FACTS */
        if (flagH) PgMaxH=HFunJac(TRUE,TRUE,Aptr,Dx);
        else if (iter!=1 && flagPoC) {
           ACFunHes(TRUE,TRUE);
           if(DCFunHes(TRUE,TRUE)) return(-iter);
           SVCFunHes(TRUE,FALSE);                   /* FACTS  */
           TCSCFunHes(TRUE,FALSE);                  /* FACTS  */
           STATCOMFunHes(TRUE,FALSE);               /* FACTS  */
        }
        if (PgMax<0 && (!flagH || (flagH && PgMaxH<0)) ) {
           if (Aptr!=NULL) {
              fprintf(stderr,"\nError: Area %d %s does not have any spinning reserves.\n",Aptr->N,Aptr->Name);
              fprintf(stderr,"       Increase the maximum P generation in this area, otherwise\n");
           } else {
              fprintf(stderr,"\nError: The system does not have any spinning reserves.\n");
              fprintf(stderr,"       Increase the maximum P generation in this system, otherwise\n");
           }
           fprintf(stderr,"       the Jacobian matrix becomes singular.\n");
           if (flagFirst) InitializeLoad();
           WriteSolution(--i,TrueParamStr(2),"Pg Max. Problems:");
           exit(1);
        }
        m=factor(Jac);
      }
      if ((i==1 && !flagR && !flagL) || m==WARNINGEXIT || flagDCLimits) {
        if(ExistParameter('d')) fprintf(stderr,"Order and factor Jacobian.\n");
        flagDCLimits=FALSE;
        if (i>1 || flagR || flagL) DeleteJac(Jac,NewRow,NewCol,OldRow,OldCol);
        Aptr=ACFunJac(Jac,&PgMax,TRUE,TRUE,flagFirst);
        if(DCFunJac(Jac,TRUE,TRUE)) return(-iter);
        SVCFunJac(Jac,TRUE,TRUE);                    /* FACTS */
        TCSCFunJac(Jac,TRUE,TRUE);                   /* FACTS */
        STATCOMFunJac(Jac,TRUE,TRUE);                /* FACTS */
        if (flagH) PgMaxH=HFunJac(TRUE,TRUE,Aptr,Dx);
        else if (iter!=1 && flagPoC) {
           ACFunHes(TRUE,TRUE);
           if(DCFunHes(TRUE,TRUE)) return(-iter);
           SVCFunHes(TRUE,FALSE);                   /* FACTS  */
           TCSCFunHes(TRUE,FALSE);                  /* FACTS  */
           STATCOMFunHes(TRUE,FALSE);               /* FACTS  */
        }
        if (PgMax<0 && (!flagH || (flagH && PgMaxH<0)) ) {
           if (Aptr!=NULL) {
              fprintf(stderr,"\nError: Area %d %s does not have any spinning reserves.\n",Aptr->N,Aptr->Name);
              fprintf(stderr,"       Increase the maximum P generation in this area, otherwise\n");
           } else {
              fprintf(stderr,"\nError: The system does not have any spinning reserves.\n");
              fprintf(stderr,"       Increase the maximum P generation in this system, otherwise\n");
           }
           fprintf(stderr,"       the Jacobian matrix becomes singular.\n");
           if (flagFirst) InitializeLoad();
           WriteSolution(--i,TrueParamStr(2),"Pg Max. Problems:");
           exit(1);
        }
        SortRowsColumns(Jac);
        if(factorns(Jac,alpha,RowPartition,ColPartition,NewRow,NewCol,OldRow,OldCol)){
           fprintf(stderr,"*** Singular Jacobian (possible voltage collapse, contol or limit problems).\n");
           fprintf(stderr,"    Try changing the load levels, controls or limits, or use the -F option.\n");
           if (flagFirst) InitializeLoad();
           WriteSolution(--i,TrueParamStr(2),"Singular Jacobian:");
           exit(1);
        }
        SortRowsColumns(Jac);
      }
      fprintf(stderr,"Iteration: %2d  ",i);
      for(j=1;j<=N;j++) dx[j]=dF[j];
      repsolp(Jac,dx,OldRow,NewCol);
      if (m==WARNINGEXIT) SD0=0;
      N2=10; k=j=0;
      while(j<=N2){
        cons= j;
        if(j==0) cons= -1; else cons=1.0/pow(2.0,cons);
        UpdateACvar(cons,j,TRUE,!ExistParameter('G'));
        UpdateDCvar(cons,j,!flag);
        UpdateSVCvar(cons,j);                        /* FACTS */
        UpdateTCSCvar(cons,j);                       /* FACTS */
        UpdateSTATCOMvar(cons,j);                    /* FACTS */
        if (iter!=1 && flagPoC) UpdateEvector(cons);
        ACFunJac(Jac,&PgMax,TRUE,FALSE,flagFirst);
        flags=DCFunJac(Jac,TRUE,FALSE);
        SVCFunJac(Jac,TRUE,FALSE);                   /* FACTS */
        TCSCFunJac(Jac,TRUE,FALSE);                  /* FACTS */
        STATCOMFunJac(Jac,TRUE,FALSE);               /* FACTS */
        if(!flags) {
            if (flagH) HFunJac(TRUE,FALSE,NULL,Dx);
            else if (iter!=1 && flagPoC) {
              ACFunHes(TRUE,FALSE);
              flags=DCFunHes(TRUE,FALSE);
              SVCFunHes(TRUE,FALSE);                   /* FACTS  */
              TCSCFunHes(TRUE,FALSE);                  /* FACTS  */
              STATCOMFunHes(TRUE,FALSE);               /* FACTS  */
           }
        }
        if (!flags) {
           val=Norm(dF,N,&N1);
           PrintMismatch(val,k++,N1);
           if (j==10 && val>2*MaxdFi) N2=20;
           else if (MaxdFi>val) {
              MaxdFi1=MaxdFi;
              MaxdFi=val;
              break;
           }
        }
        j++;
      }
      /* if (ExistParameter('d')) fprintf(stderr,"j %d   MaxdFi %lf   Tol %lf   MaxdFi1 %lf   tol %lf\n",j,MaxdFi,Tol,MaxdFi1,tol); */
      if ((MaxdFi>Tol) && (j>N2 || (fabs(MaxdFi1-MaxdFi)/MaxdFi)<=tol)) {
        if (flagL) return(-(++i));
        if (!flagR) {
          flagp=CheckRlimits();
          /*  Apply Q limits after convergence
          if (!ExistParameter('G') && !flagp) flagp=CheckQlimits();
          */
          if (!flagp) flagp=CheckVlimits();
          if (!flagp) flagp=CheckQlimits();
          if (!flagp) flagp=flag=CheckDClimits();
          if (!flagp) flagp=flagDCLimits=ChangeSVCmode();       /* FACTS */
          if (!flagp) flagp=flagDCLimits=ChangeTCSCmode();      /* FACTS */
          if (!flagp) flagp=flagDCLimits=ChangeSTATCOMmode();  /* FACTS */
        } else {
          flagp=CheckRlimits();
          if (!flagp) flagp=CheckVlimits(); else CheckVlimits();
          if (!flagp) flagp=CheckQlimits(); else CheckQlimits();
          if (!flagp) flagp=flagDCLimits=ChangeDCmode(); else flagDCLimits=ChangeDCmode();
          if (!flagp) flagp=flagDCLimits=ChangeSVCmode(); else flagDCLimits=ChangeSVCmode();          /* FACTS */
          if (!flagp) flagp=flagDCLimits=ChangeTCSCmode(); else flagDCLimits=ChangeTCSCmode();        /* FACTS */
          if (!flagp) flagp=flagDCLimits=ChangeSTATCOMmode(); else flagDCLimits=ChangeSTATCOMmode();  /* FACTS */
          flagL=TRUE;
       }
       if(!flagp) {
          if (flagR) return(-(++i));
          MaxdFi=val;
          fprintf(stderr,"\n *** The case diverges (possible voltage collapse or AC/DC/FACTS control\n");
          fprintf(stderr,"     problems).  Try changing the load levels or AC/DC/FACTS controls, or\n");
          fprintf(stderr,"     use the -F option or decrease the tolerance between two\n");
          fprintf(stderr,"     consecutive iterations with the -t option.\n");
          if (flagFirst) InitializeLoad();
          WriteSolution(i,TrueParamStr(2),"Divergence:");
          exit(1);
        }
      }
    } else break;
  }
  /* if (ExistParameter('d')) fprintf(stderr,"i %d   MaxdFi %lf   Tol %lf\n",i,MaxdFi,Tol); */
  if (i>MaxIterRun && MaxdFi>Tol) {
     if (flagR) return(-i);
     fprintf(stderr,"\n *** The case has not been solved (possible voltage collapse, AC/DC/FACTS control\n");
     fprintf(stderr,"     problems, or too few iterations).  Try running the case using the -F\n");
     fprintf(stderr,"     option, or change load levels, AC/DC/FACTS controls, or increase the maximum\n");
     fprintf(stderr,"     number of iterations with the -M option.\n");
     if (flagFirst) InitializeLoad();
     WriteSolution(--i,TrueParamStr(2),"Unsolved case:");
     exit(1);
  }
  /*    Apply Q limits after convergence
  N3=i;
  if (ExistParameter('G') && !flagR) {
     flagp=CheckQlimits();
     if (flagp) goto NewtonRapson;
  }
  */
  if (flagFirst) InitializeLoad();
  return(i);
}