tefacdc.cpp

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

/* AC/DC Transient Energy Functions. */

#include <stdlib.h>
#ifndef WINDOWS
#include <stdio.h>
#else
#include "pfwstdio.h"
#endif
#include <math.h>
#include <string.h>
#include "constant.h"
#include "param.h"
#include "sparse.h"
#include "pflow.h"

#ifdef ANSIPROTO
void TEFac(BOOLEAN flag);
void TEFdc(FILE *Out);
void MatlabV(FILE *Out);
void TEFMatlabFiles(void);
void Print(FILE *File,int spaces,int width,int decimals,VALUETYPE val);
#else
void TEFac();
void TEFdc();
void MatlabV();
void TEFMatlabFiles();
void Print();
#endif

/* ------- Global Variables ------ */
extern Data *dataPtr;
extern SparseMatrix *Jac;
extern INDEX Nac,Ngen,Ndc,NacVar,Narea;
extern INDEX *ACvar;
extern VALUETYPE *x0p,lambda,Sn;
extern AClist *Vlist;
extern BOOLEAN flagBS,flagPgMax;
VALUETYPE Vac;
extern int field;
extern BOOLEAN flagPrintTotalPl,flagPrintTotalQl,flagPrintTotalPg,flagPrintTotalQg;


/* ------------------ TEFac ----------------------------- */
#ifdef ANSIPROTO
void TEFac(BOOLEAN flag)
#else
void TEFac(flag)
BOOLEAN flag;
#endif
/* Caculate AC TEF. */
{
  ACbusData *ACptr,*BEptr,*To;
  ElementList *ELptr;
  ElementData *Eptr;
  VALUETYPE Vi,Vj,di,dj,Pi,Qi,Pl,Ql,val;
  VALUETYPE gii,bii,gij,bij,gsij,bsij,DPg,Pg;
  INDEX i,j;

  if (Narea<2){
    for(BEptr=dataPtr->ACbus;BEptr!=NULL;BEptr=BEptr->Next)
      if(strpbrk(BEptr->Type,"S")) break;
  }
  for(Vac=0,ACptr=dataPtr->ACbus; ACptr!=NULL; ACptr=ACptr->Next) {
    if(ACptr->Area!=NULL) BEptr=ACptr->Area->Slack;
    if (flagBS) {if (ACptr!=BEptr) DPg=0; else DPg=1;}
    else DPg=ACptr->DPg;
    i=ACvar[ACptr->N];
    Vi=ACptr->V;
    di=ACptr->Ang;
    if (flag) {
      x0p[i]=di;
      x0p[i+1]=Vi;
    }
    Pl = (ACptr->Pn+lambda*ACptr->Pnl)*pow(Vi,ACptr->a);
    Ql = (ACptr->Qn-lambda*ACptr->Qnl)*pow(Vi,ACptr->b);
    Pg=ACptr->Pg+DPg*BEptr->Kg;
    if (!flagPgMax && Pg>ACptr->PgMax) Pg=ACptr->PgMax;
    Pi = Pg-Pl;
    Qi = -Ql;
    Qi=Qi+ACptr->Qg;
    gii=ACptr->G+ACptr->Pz+lambda*ACptr->Pzl;
    bii=ACptr->B+ACptr->Qz+lambda*ACptr->Qzl;
    for(val=0,ELptr=ACptr->Elem; ELptr!=NULL; ELptr=ELptr->Next) {
      Eptr=ELptr->Eptr;
      if (Eptr->From==ACptr) To=Eptr->To;
      else To=Eptr->From;
      j=ACvar[To->N];
      Vj=To->V;
      dj=To->Ang;
      if (flag) {
        x0p[j]=dj;
        x0p[j+1]=Vj;
      }
      if(Eptr->From==ACptr) {
        gij=(Eptr->G*cos(Eptr->Ang)-Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        bij=(Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
        gsij=(Eptr->G1+Eptr->G)*pow(Eptr->Tap,2.0)-gij;
        bsij=(Eptr->B1+Eptr->B)*pow(Eptr->Tap,2.0)-bij;
      } else {
        gij=(Eptr->G*cos(Eptr->Ang)+Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        bij=(-Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
        gsij=Eptr->G+Eptr->G2-gij;
        bsij=Eptr->B+Eptr->B2-bij;
      }
      val=val+0.5*bij*Vi*Vj*cos(di-dj);
      val=val-gij*x0p[i+1]*x0p[j+1]*cos(x0p[i]-x0p[j])*di-gij*x0p[j+1]*sin(x0p[i]-x0p[j])*Vi;
      gii=gii+gij+gsij;
      bii=bii+bij+bsij;
    }
    val=val-0.5*bii*Vi*Vi-Pi*di;
    if (Vi>0) val=val-Qi*log(Vi);
    val=val+gii*x0p[i+1]*x0p[i+1]*di;
    Vac=Vac+val;
  }
}

/* ------------------ TEFdc ----------------------------- */
#ifdef ANSIPROTO
void TEFdc(FILE *Out)
#else
void TEFdc(Out)
FILE *Out;
#endif
/* Caculate DC TEF and coupling terms. */
{
  DCbusData *DCptrR,*DCptrI,*DCptr;
  VALUETYPE a,b,d,t,cosag,Id,Xc;
  INDEX j;

  for(DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next){
    DCptrI=DCptrR->To;
    if(!strcmp(DCptrR->Type,"R")){
      for (j=1;j<=2;j++) {
        if (j==1) {
          DCptr=DCptrR;
          cosag=cos(DCptr->Alfa);
        } else {
          DCptr=DCptrI;
          cosag=cos(DCptr->Gamma);
        }
        Id=DCptr->Id;
        t=DCptr->Tap*DCptr->Ntrf;
        Xc=DCptr->Xc;
        a= -Xc*Xc*Id*Id;
        b=2*sqrt(2.0)*t*Xc*Id*cosag;
        /* c=2*t*t*(1-cosag*cosag); */
        d=8*t*t*Xc*Xc*Id*Id;
        Print(Out,0,field,4,a);  fprintf(Out,"\t");
        Print(Out,0,field,4,b);  fprintf(Out,"\t");
        Print(Out,0,field,4,d);  fprintf(Out,"\t");
      }
    }
  }
}


/* --------------------------- MatlabV --------------------------------- */
#ifdef ANSIPROTO
void MatlabV(FILE *Out)
#else
void MatlabV(Out)
FILE *Out;
#endif
/* Write Matlab commands in output file for ploting profiles. */
{
  AClist *Lptr;
  char LineType[4][5];
  DCbusData *DCptr,*DCptrp;
  INDEX count,countp,i,j,k,l,js,ks;
  VALUETYPE k1,k2,ki,kp,P11,P12,P22,P23,P33,Beta,KV;

  strcpy(LineType[0],"'-'");
  strcpy(LineType[1],"'-.'");
  strcpy(LineType[2],"':'");
  strcpy(LineType[3],"'--'");
  fprintf(Out,"];\n");
  fprintf(Out,"%c Plot profiles:\n",'%');
  fprintf(Out,"%c Change the value of K to scale L.F.\nK=1;\n",'%');
  fprintf(Out,"figure;\nplot(");
  for(count=countp=0,Lptr=Vlist;Lptr!=NULL;Lptr=Lptr->Next) {
    count++;
    if (countp<4) {
      if  (!strcmp(Lptr->Type,"V")) {
        KV=Lptr->AC->KV;
        if (KV<=0) KV=1;
        if (countp==0) fprintf(Out,"K*x(:,1),%5.1f*x(:,%1d)",KV,count+1);
        else           fprintf(Out,",K*x(:,1),%5.1f*x(:,%1d)",KV,count+1);
        fprintf(Out,",%s", LineType[countp]);
      } else {
        if (countp==0) fprintf(Out,"K*x(:,1),x(:,%1d)",count+1);
        else           fprintf(Out,",K*x(:,1),x(:,%1d)",count+1);
        fprintf(Out,",%s", LineType[countp]);
      }
    }
    countp++;
  }
  if (countp<4 && flagPrintTotalPl) {
    if (countp==0) fprintf(Out,"K*x(:,1),x(:,%1d)",count+1);
    else           fprintf(Out,",K*x(:,1),x(:,%1d)",count+1);
    fprintf(Out,",%s", LineType[countp]);
    countp++;
  }
  if (countp<4 && flagPrintTotalQl) {
    if (countp==0) fprintf(Out,"K*x(:,1),x(:,%1d)",count+1);
    else           fprintf(Out,",K*x(:,1),x(:,%1d)",count+1);
    fprintf(Out,",%s", LineType[countp]);
    countp++;
  }
  if (countp<4 && flagPrintTotalPg) {
    if (countp==0) fprintf(Out,"K*x(:,1),x(:,%1d)",count+1);
    else           fprintf(Out,",K*x(:,1),x(:,%1d)",count+1);
    fprintf(Out,",%s", LineType[countp]);
    countp++;
  }
  if (countp<4 && flagPrintTotalQg) {
    if (countp==0) fprintf(Out,"K*x(:,1),x(:,%1d)",count+1);
    else           fprintf(Out,",K*x(:,1),x(:,%1d)",count+1);
    fprintf(Out,",%s", LineType[countp]);
    countp++;
  }
  fprintf(Out,");\n");
  fprintf(Out,"legend(");
  for(countp=0,Lptr=Vlist;Lptr!=NULL;Lptr=Lptr->Next) {
    if (countp<4) {
      if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"V")) {
        if (countp==0) fprintf(Out,"'kV_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'kV_{%s}'",Lptr->AC->Name);
      }
      else if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"D")) {
        if (countp==0) fprintf(Out,"'deg_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'deg_{%s}'",Lptr->AC->Name);
      }
      else if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"PL")) {
        if (countp==0) fprintf(Out,"'L MW_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'L MW_{%s}'",Lptr->AC->Name);
      }
      else if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"QL")) {
        if (countp==0) fprintf(Out,"'L MVar_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'L MVar_{%s}'",Lptr->AC->Name);
      }
      else if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"PG")) {
        if (countp==0) fprintf(Out,"'G MW_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'G MW_{%s}'",Lptr->AC->Name);
      }
      else if (Lptr->AC!=NULL && !strcmp(Lptr->Type,"QG")) {
        if (countp==0) fprintf(Out,"'G MVar_{%s}'",Lptr->AC->Name);
        else           fprintf(Out,",'G MVar_{%s}'",Lptr->AC->Name);
      }
      else if(Lptr->Area!=NULL) {
        if (countp==0) fprintf(Out,"'A MW_{%s}'",Lptr->Type,Lptr->Area->Name);
        else           fprintf(Out,",'A MW_{%s}'",Lptr->Type,Lptr->Area->Name);