readepri.cpp

用于潮流计算的程序请让我下载我需要的吧感谢了啊

/* Read AC data in WSCC/EPRI/BPA/ANAREDE format
   (EPRI's P.F. V 5.0 input data formats). */

#include "readdata.h"

/* ------- Global Variables ------ */
VALUETYPE G,B,G1,B1,G2,B2,Tap,Ang;
extern FILE *InputDataFile;


/* ---------------- Multiply 2 Complex Numbers ----------------------------- */
#ifdef ANSIPROTO
void Multiply(VALUETYPE *a,VALUETYPE *b,VALUETYPE c,VALUETYPE d)
#else
void Multiply(a,b,c,d)
VALUETYPE *a,*b,c,d;
#endif
{
  VALUETYPE r,i;

  r=(*a)*c - (*b)*d;
  i=(*b)*c + (*a)*d;
  *a=r; *b=i;
}

/* ---------------- Divide 2 Complex Numbers ----------------------------- */
#ifdef ANSIPROTO
void Divide(VALUETYPE *a,VALUETYPE *b,VALUETYPE c,VALUETYPE d)
#else
void Divide(a,b,c,d)
VALUETYPE *a,*b,c,d;
#endif
{
  VALUETYPE D;

  D=c*c+d*d;
  if (!D) { c=d=0; D=1;}
  Multiply(a,b,c/D,-d/D);
}


/* ---------------- AddSection --------------------------- */
#ifdef ANSIPROTO
BOOLEAN AddSection(ACbusData *From,ACbusData *To,char *Line,char *Ckt,INDEX Sec)
#else
BOOLEAN AddSection(From,To,Line,Ckt,Sec)
ACbusData *From,*To;
char *Line,*Ckt;
INDEX Sec;
#endif
{
  ElementList *ELptr;
  ElementData *Eptr;
  VALUETYPE Ge,Be,Ge1,Be1,Ge2,Be2,Gs,Bs;

  for(ELptr=From->Elem;ELptr!=NULL;ELptr=ELptr->Next){
    Eptr=ELptr->Eptr;
    if (((From==Eptr->From&&To==Eptr->To)||(From==Eptr->To&&To==Eptr->From))
        &&!strcmp(Ckt,Eptr->Ckt)&&strcmp(Eptr->Zone,"")){
      if (!strncmp(Line,"T",1)&&!strpbrk(Eptr->Type,"RT")) strcpy(Eptr->Type,"T");
      if(Sec<Eptr->Sec) {
        Ge=Eptr->G;   Eptr->G=G/(Eptr->Tap*Eptr->Tap);   G=Ge;
        Be=Eptr->B;   Eptr->B=B/(Eptr->Tap*Eptr->Tap);   B=Be;
        Ge1=Eptr->G1; Eptr->G1=G1/(Eptr->Tap*Eptr->Tap); G1=Ge1;
        Be1=Eptr->B1; Eptr->B1=B1/(Eptr->Tap*Eptr->Tap); B1=Be1;
        Ge2=Eptr->G2; Eptr->G2=G2/(Eptr->Tap*Eptr->Tap); G2=Ge2;
        Be2=Eptr->B2; Eptr->B2=B2/(Eptr->Tap*Eptr->Tap); B2=Be2;
      } else {
        Eptr->G=Eptr->G/(Tap*Tap);
        Eptr->B=Eptr->B/(Tap*Tap);
        Eptr->G1=Eptr->G1/(Tap*Tap);
        Eptr->B1=Eptr->B1/(Tap*Tap);
        Eptr->G2=Eptr->G2/(Tap*Tap);
        Eptr->B2=Eptr->B2/(Tap*Tap);
      }
      if ((fabs(Eptr->G)<1e-7 && fabs(Eptr->B)<1e-7)||(fabs(G)<1e-7 && fabs(B)<1e-7)) {
         fCustomPrint(stderr,"***Warning: A section of element %d %s - %d %s\n",
                  Eptr->From->Num,Eptr->From->Name,Eptr->To->Num,Eptr->To->Name);
         fCustomPrint(stderr,"            has a zero impedance.  Check the section data for this element.\n");
      }
      if (!strncmp(Line,"T",1)&&!strpbrk(Eptr->Type,"RT")) strcpy(Eptr->Type,"T");
      Eptr->Tap*=Tap;
      Eptr->Ang+=Ang;
      Eptr->Sec=Sec;
      Gs=Eptr->G+G+Eptr->G2+G1;
      Bs=Eptr->B+B+Eptr->B2+B1;
      Ge=Ge1=Eptr->G;  Be=Be1=Eptr->B;
      Ge2=G;  Be2=B;
      Multiply(&Ge,&Be,G,B);
      Divide(&Ge,&Be,Gs,Bs);
      Multiply(&Ge1,&Be1,Eptr->G2+G1,Eptr->B2+B1);
      Divide(&Ge1,&Be1,Gs,Bs);
      Ge1+=Eptr->G1;  Be1+=Eptr->B1;
      Multiply(&Ge2,&Be2,Eptr->G2+G1,Eptr->B2+B1);
      Divide(&Ge2,&Be2,Gs,Bs);
      Ge2+=G2;       Be2+=B2;
      Eptr->G=Ge;    Eptr->B=Be;
      Eptr->G1=Ge1;  Eptr->B1=Be1;
      Eptr->G2=Ge2;  Eptr->B2=Be2;
      return(TRUE);
    }
  }
  return(FALSE);
}


/* ---------------- ReadWSCC ----------------------------- */
#ifdef ANSIPROTO
void ReadWSCC()
#else
void ReadWSCC()
#endif
/* Read Bus and Element data in WSCC format. */
{
  ACbusData *ACptr,*ACptrp,*ACptrs;
  ElementData *Eptr;
  AreaData *Aptr;
  char Line[BUFLEN],Name[31],str[6];
  VALUETYPE KV,KVp,KVs,R,X,Bx,Taps,Tap1,Tap2,Imax;
  int i,j,k,s;
  INDEX Sec;
  BOOLEAN flag=FALSE,flagPrint=TRUE;


  for(i=0;i<=2;strcpy(dataPtr->Title[i],"\n"),i++);
  Sn=100.0;
  RealParameter('$',&Sn,1.0,100000000.0);
  ACptr=NULL;
  for(;;){ /* Reading Loop */
    if (fgets(Line,BUFLEN,InputDataFile)==NULL){
      ErrorHalt("Missing END or Anarede's FIM card.");
      break;
    }
    LineNum++;
    if (!strncmp(Line,"C",1) || !strncmp(Line,"BAS",3)) continue;

  /* ------------------------ Title ----------------------------- */
    else if (!strncmp(Line,"HDG",3)||!strncmp(Line,"TITU",4)) {
      i=0;
      for(;;){ /* Title Reading Loop */
        if (fgets(Line,BUFLEN,InputDataFile)==NULL){
          ErrorHalt("Missing BAS or Anarede's D title cards.");
          break;
        }
        LineNum++;
        if (strncmp(Line,"BAS",3)&&strncmp(Line,"D",1)) {
          strcpy(dataPtr->Title[i],Line);
          i++;
          if (i>2) break;
        }
        else break;
      }
    }

  /* --------------- AC bus data -------------------------------- */
    else if (!strncmp(Line,"B  ",3) || !strncmp(Line,"BQ ",3) ||
        !strncmp(Line,"BV ",3) || !strncmp(Line,"BE ",3) ||
        !strncmp(Line,"BC ",3) || !strncmp(Line,"BG ",3) ||
        !strncmp(Line,"BT ",3) || !strncmp(Line,"BX ",3) ||
        !strncmp(Line,"BA ",3) || !strncmp(Line,"BF ",3) ||
        !strncmp(Line,"BJ ",3) || !strncmp(Line,"BK ",3) ||
        !strncmp(Line,"BL ",3) || !strncmp(Line,"BS ",3)){
      GetStr(Line,7,12,12,Name);
      KV=GetValue(Line,15,4,0);
      ACptr=ACbusInList(0,Name,KV,Nac,1);
      if (ACptr->V>0) {
         fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
         ErrorHalt("The AC bus was previously defined (check B cards).");
      }
      if (ACptr->N==0) { Nac++; ACptr->Num=ACptr->N=Nac;}
      GetStr(Line,19,2,2,ACptr->Zone);
      GetStr(Line,4,3,3,ACptr->Owner);
      ACptr->V=GetValue(Line,58,4,3);
      if (ACptr->V<=0) ACptr->V=1;
      ACptr->VCont=ACptr->V;
      ACptr->Pl=GetValue(Line,21,5,0)/Sn;
      ACptr->Ql=GetValue(Line,26,5,0)/Sn;
      ACptr->G=GetValue(Line,31,4,0)/Sn;
      ACptr->B=GetValue(Line,35,4,0)/Sn;
      ACptr->Pg=GetValue(Line,43,5,0)/Sn;
      ACptr->PgMax=GetValue(Line,39,4,0)/Sn;
      if (ACptr->PgMax==0) ACptr->PgMax=99999999.;
      if (ACptr->PgMax<ACptr->Pg) {
         ACptr->PgMax=99999999.;
         fCustomPrint(stderr,"***Warning: Bus %d %s has its maximum generating power PgMax < Pg.\n",ACptr->N,ACptr->Name);
         fCustomPrint(stderr,"            PgMax will be given value of 99999999.\n");
      }
      ACptr->Pmax=ACptr->PgMax;
      ACptr->Qg=GetValue(Line,48,5,0)/Sn;
      if ((!strncmp(Line,"B  ",3)||!strncmp(Line,"BJ ",3)||!strncmp(Line,"BF ",3))
          && strcmp(ACptr->Type,"BC")) {
         ACptr->Cont=ACptr;
         if (!strncmp(Line,"BJ ",3)){
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            fCustomPrint(stderr,"***Warning: BJ bus has been transformed to a B bus.\n");
         }
         if (!strncmp(Line,"BF ",3)){
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            fCustomPrint(stderr,"***Warning: BF bus has been transformed to a B bus.\n");
         }
      }
      else if (!strncmp(Line,"BC ",3)) strcpy(ACptr->Type,"BC");
      else if (!strncmp(Line,"BT ",3)) strcpy(ACptr->Type,"BT");
      else if (!strncmp(Line,"BE ",3) || !strncmp(Line,"BK ",3)){
         Nvolt++;
         ACptr->Qg=0;
         strcpy(ACptr->Type,"BQ");
         strcpy(ACptr->cont,"V");
         ACptr->Qmax=99999999.;
         ACptr->Qmin= -99999999.;
         if (!strncmp(Line,"BK ",3)){
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            fCustomPrint(stderr,"***Warning: BK bus has been transformed to a BE bus.\n");
         }
      }
      if (!strncmp(Line,"BQ ",3) || !strncmp(Line,"BK ",3)){
         Nvolt++;
         ACptr->Qg=0;
         strcpy(ACptr->Type,"BQ");
         strcpy(ACptr->cont,"V");
         ACptr->Qmax=GetValue(Line,48,5,0)/Sn;
         ACptr->Qmin=GetValue(Line,53,5,0)/Sn;
         if (ACptr->Qmax<=ACptr->Qmin) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("AC bus Q limits are wrong: Qmin >= Qmax.");
         }
         if (!strncmp(Line,"BL ",3)){
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            fCustomPrint(stderr,"***Warning: BL bus has been transformed to a BQ bus.\n");
         }
      }
      else if (!strncmp(Line,"BG ",3)){
         Nvolt++;
         ACptr->Qg=0;
         strcpy(ACptr->Type,"BG");
         strcpy(ACptr->cont,"V");
         ACptr->Qmax=GetValue(Line,48,5,0)/Sn;
         ACptr->Qmin=GetValue(Line,53,5,0)/Sn;
         if (ACptr->Qmax<=ACptr->Qmin) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("AC bus Q limits are wrong: Qmin >= Qmax.");
         }
         GetStr(Line,66,12,12,Name);
         KV=GetValue(Line,74,4,0);
         if (!strcmp(Name,"            ")) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("The controlled bus has not been defined.");
         }
         ACptrp=ACbusInList(0,Name,KV,Nac,1);
         if (ACptrp->N==0) { Nac++; ACptrp->Num=ACptrp->N=Nac;}
         if (!strcmp(ACptr->Name,"            ")) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("The controlled bus has not been defined.");
         }
         if(!strcmp(ACptrp->Type,"B")) strcpy(ACptrp->Type,"BC");
         else if(strcmp(ACptrp->Type,"BC")) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("The voltage controlled bus is not a PQ bus.");
         }
         ACptr->Cont=ACptrp;
         ACptrp->Cont=NULL;
         ACptrp->Kbg++;
         ACptr->Kbg=GetValue(Line,78,3,0)/100;
         if (flag2Vcontrol) {
           ACptrp->Kbg1=ACptrp->Kbg1+ACptr->Qmax;
           ACptrp->Kbg2=ACptrp->Kbg2+ACptr->Qmin;
         }
         if (ACptr->Kbg<=0) ACptr->Kbg=1.;
      }
      else if (!strncmp(Line,"BV ",3)||!strncmp(Line,"BA ",3)) {
         strcpy(ACptr->Type,"BV");
         strcpy(ACptr->cont,"Q");
         ACptr->VCont=ACptr->Qg;
         ACptr->Vmax=GetValue(Line,58,4,3);
         if (ACptr->Vmax<=0) ACptr->Vmax=10.;
         ACptr->Vmin=GetValue(Line,62,4,3);
         if (ACptr->Vmin<=0) ACptr->Vmin=0.001;
         if (ACptr->Vmax<=ACptr->Vmin) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("AC bus V limits are wrong: Vmin >= Vmax.");
         }
         ACptr->Qmax=99999999.;
         ACptr->Qmin= -99999999.;
         ACptr->Cont=ACptr;
         if (!strncmp(Line,"BA ",3)){
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            fCustomPrint(stderr,"***Warning: BA bus has been transformed to a BV bus.\n");
         }
      }
      else if (!strncmp(Line,"BX ",3)) {
         ACptr->Vmax=GetValue(Line,58,4,3);
         if (ACptr->Vmax<=0) ACptr->Vmax=10.;
         ACptr->Vmin=GetValue(Line,62,4,3);
         if (ACptr->Vmin<=0) ACptr->Vmin=0.001;
         if (ACptr->Vmax<=ACptr->Vmin) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("AC bus V limits are wrong: Vmin >= Vmax.");
         }
         ACptr->Cont=ACptr;
      }
      else if (!strncmp(Line,"BS ",3)) {
         strcat(ACptr->Type,"S");
         strcpy(ACptr->cont,"V");
         Nslack++;
         ACptr->Qg=0;
         AngSlack=0;
         ACptr->Qmax=GetValue(Line,48,5,0)/Sn;
         ACptr->Qmin=GetValue(Line,53,5,0)/Sn;
         if (ACptr->Qmax<=ACptr->Qmin) {
            fCustomPrint(stderr,"Input Line-> %d\n%s",LineNum,Line);
            ErrorHalt("AC bus Q limits are wrong: Qmin >= Qmax.");
         }
      }