pointac.cpp

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

/* Point of Collapse:  AC Hessian. */

#include "pointl.h"


/* ------------------ ACFunHes ----------------------------- */
#ifdef ANSIPROTO
void ACFunHes(BOOLEAN flagF,BOOLEAN flagJ)
#else
void ACFunHes(flagF,flagJ)
BOOLEAN flagF,flagJ;
#endif
/* Construct the AC part of the PoC Jacobian. */
{
  SparseMatrixElement *Jptr;
  ACbusData *ACptr,*To,*From,*BEptr,*BSptr;
  AClist *ALptr;
  ElementList *ELptr;
  ElementData *Eptr;
  VALUETYPE Vi,Vj,di,dj,gij,bij,gsij,bsij,gji,bji;
  VALUETYPE v1,dv1,v2,dv2,v3,dv3,v4,dv4,v5,v6,dv6,v7,v8,dv8,v9,v10,dv10,v11;
  VALUETYPE Ddi,Dvi,Ddip,Dvip,Ddj,Dvj,Dval,a,b,Dlambda,DPg;
  VALUETYPE Ra,Xd,Xq,dg,Vd,Vq,Vr,Vim,Iq,Id,Ir,Iim,Ia;
  INDEX N,i,j,k,l;

  N=NacVar+11*Ndc/2+3*Nsvc+NtcscVar+7*Nstatcom;     /*  FACTS  */
  for(v1=0,i=1+N;i<Jac->n1;i++) {
    if (fabs(x0[i-N])>v1) {
      k=i-N;
      v1=fabs(x0[k]);
      if(x0[k]>0) v2=1.0;
      else v2= -1.0;
    }
    if (flagJ) JacElement(Jac,Jac->n1,i,0.);
    if (flagF) dF[i]=0;
  }
  if (flagF) dF[Jac->n1]=v1-1;
  if (flagJ) JacElement(Jac,Jac->n1,k+N,v2);
  if (flagJ) for(i=1;i<=Jac->n1;i++) {
    Jptr=Jac->ColHead[i];
    while(Jptr!=NULL) {
      j=Jptr->Col;
      if (OldCol->p[j]) j=OldCol->p[j];
      l=Jptr->Row;
      if (OldRow->p[l]) l=OldRow->p[l];
      if (flagJ && j<=N && l<=N) JacElement(Jac,j+N,l+N,Jptr->Value);
      Jptr=Jptr->ColNext;
    }
  }
  for(ALptr=dataPtr->KGbus;ALptr!=NULL;ALptr=ALptr->Next) {
    BSptr=ALptr->AC;
    if(strpbrk(ALptr->AC->Type,"S"))  break;
  }
  for(ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next){
    if (ACptr->Area!=NULL) BSptr=ACptr->Area->Slack;
    i=ACvar[ACptr->N];
    Vi=ACptr->V;
    di=ACptr->Ang;
    a=ACptr->a;
    b=ACptr->b;
    From=ACptr;
    Ddi=Dvi=Ddip=Dvip=0;
    /*if (flagF) dF[Jac->n1]+=-ACptr->DPl*x0[i]-ACptr->DQl*x0[i+1];*/
    for(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(Eptr->From==From) {
        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)*Eptr->Tap*Eptr->Tap-gij;
        bsij=(Eptr->B1+Eptr->B)*Eptr->Tap*Eptr->Tap-bij;
        gji=(Eptr->G*cos(Eptr->Ang)+Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        bji=(-Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
      } 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;
        gji=(Eptr->G*cos(Eptr->Ang)-Eptr->B*sin(Eptr->Ang))*Eptr->Tap;
        bji=(Eptr->G*sin(Eptr->Ang)+Eptr->B*cos(Eptr->Ang))*Eptr->Tap;
      }
      v1=sin(di-dj)*(gij*x0[i]+gji*x0[j]);
      dv1=cos(di-dj)*(gij*x0[i]+gji*x0[j]);
      v2=cos(di-dj)*(bij*x0[i]-bji*x0[j]);
      dv2= -sin(di-dj)*(bij*x0[i]-bji*x0[j]);
      v3=cos(di-dj)*(gij*x0[i+1]-gji*x0[j+1]);
      dv3= -sin(di-dj)*(gij*x0[i+1]-gji*x0[j+1]);
      v4=sin(di-dj)*(bij*x0[i+1]+bji*x0[j+1]);
      dv4=cos(di-dj)*(bij*x0[i+1]+bji*x0[j+1]);
      v5= -2*((gij+gsij)*x0[i]-(bij+bsij)*x0[i+1]);
      if (Acont && From->Area!=To->Area) {
        BEptr=From->Area->Slack;
        if(!strpbrk(BEptr->Type,"S")) {
          k=ACvar[BEptr->N]+2;
          if (Eptr->Meter==From) {
            v6= -(gij*sin(di-dj)-bij*cos(di-dj))*x0[k];
            dv6= -(gij*cos(di-dj)+bij*sin(di-dj))*x0[k];
            v7= -2*(gij+gsij)*x0[k];
          } else {
            v6= -(gji*sin(dj-di)-bji*cos(dj-di))*x0[k];
            dv6=(gji*cos(dj-di)+bji*sin(dj-di))*x0[k];
            v7=0;
          }
        } else v6=dv6=v7=0;
        BEptr=To->Area->Slack;
        if(!strpbrk(BEptr->Type,"S")) {
          k=ACvar[BEptr->N]+2;
          if (Eptr->Meter==From) {
            v8=(gij*sin(di-dj)-bij*cos(di-dj))*x0[k];
            dv8=(gij*cos(di-dj)+bij*sin(di-dj))*x0[k];
            v9=2*(gij+gsij)*x0[k];
          } else {
            v8=(gji*sin(dj-di)-bji*cos(dj-di))*x0[k];
            dv8= -(gji*cos(dj-di)+bji*sin(dj-di))*x0[k];
            v9=0;
          }
        } else v8=dv8=v9=0;
      } else v6=dv6=v7=v8=dv8=v9=0;
      if (PQcont && strpbrk(Eptr->Type,"PM")){
        k=ACvar[Eptr->Cont->N]+1+Eptr->Cont->Ncont-Eptr->Ncont;
        if (Acont && strpbrk(Eptr->Cont->Type,"A")) k++;
        if (Eptr->Cont==From) {
          v10= -(gij*sin(di-dj)-bij*cos(di-dj))*x0[k];
          dv10= -(gij*cos(di-dj)+bij*sin(di-dj))*x0[k];
          v11= -2*(gij+gsij)*x0[k];
        } else {
          v10=(gji*sin(dj-di)-bji*cos(dj-di))*x0[k];
          dv10= -(gji*cos(dj-di)+bji*sin(dj-di))*x0[k];
          v11=0;
        }

      } else if (PQcont && strpbrk(Eptr->Type,"QN")){
        k=ACvar[Eptr->Cont->N]+1+Eptr->Cont->Ncont-Eptr->Ncont;
        if (Acont && strpbrk(Eptr->Cont->Type,"A")) k++;
        if (Eptr->Cont==From) {
          v10=(gij*cos(di-dj)+bij*sin(di-dj))*x0[k];
          dv10= -(gij*sin(di-dj)-bij*cos(di-dj))*x0[k];
          v11=2*(bij+bsij)*x0[k];
        } else {
          v10= -(gji*cos(dj-di)+bji*sin(dj-di))*x0[k];
          dv10= -(gji*sin(dj-di)-bji*cos(dj-di))*x0[k];
          v11=0;
        }
      } else v10=dv10=v11=0;
      if(flagF){
        if (!strpbrk(From->Type,"S")) dF[i+N]+=Vi*Vj*(-v1+v2+v3+v4+v6+v8+v10);
        if (From->Cont!=NULL)         dF[i+N+1]+=Vj*(dv1-dv2-dv3-dv4-dv6-dv8-dv10)+Vi*(v5+v7+v9+v11);
      }
      if(flagJ){
        if (!strpbrk(From->Type,"S")){
          Ddj=Vi*Vj*(dv1-dv2-dv3-dv4-dv6-dv8-dv10);
          Ddi=Ddi-Ddj;
          Dvj=Vi*(-v1+v2+v3+v4+v6+v8+v10);
          Dvi=Dvi+Vj*(-v1+v2+v3+v4+v6+v8+v10);
          if (!strpbrk(To->Type,"S"))  JacElement(Jac,i+N,j,Ddj);
          if (To->Cont!=NULL)  JacElement(Jac,i+N,j+1,Dvj);
          else JacElement(Jac,i+N,j+1,0.);
          if (PQcont && strpbrk(Eptr->Type,"PQMN")){
            k=ACvar[Eptr->Cont->N]+1+Eptr->Cont->Ncont-Eptr->Ncont;
            if (Acont && strpbrk(Eptr->Cont->Type,"A")) k++;
            if(!strcmp(Eptr->Type,"RQ"))
              Dval=Vi*Vj/Eptr->Tap*(-v1+v2+v3+v4+v6+v8+v10);
            else if(!strcmp(Eptr->Type,"RP"))
              Dval=Vi*Vj*(-dv2+dv1-dv4-dv3-dv6-dv8-dv10);
            else Dval=0;
            JacElement(Jac,i+N,k,Dval);
          }
          else if(Rcont && Eptr->Cont!=NULL) {
            k=ACvar[Eptr->Cont->N]+1;
            if(!strcmp(Eptr->Type,"R"))
              Dval=Vi*Vj/Eptr->Tap*(-v1+v2+v3+v4+v6+v8+v10);
            else Dval=0;
            JacElement(Jac,i+N,k,Dval);
          }
        }
        if (From->Cont!=NULL) {
          Ddj= -Vj*(-v1+v2+v3+v4+v6+v8+v10);
          Ddip=Ddip-Ddj;
          Dvj=dv1-dv2-dv3-dv4-dv6-dv8-dv10;
          Dvip=Dvip+v5+v7+v9+v11;
          if (!strpbrk(To->Type,"S"))  JacElement(Jac,i+1+N,j,Ddj);
          if (To->Cont!=NULL) JacElement(Jac,i+1+N,j+1,Dvj);
          else JacElement(Jac,i+1+N,j+1,0.);
          if (PQcont && strpbrk(Eptr->Type,"PQMN")){
            k=ACvar[Eptr->Cont->N]+1+Eptr->Cont->Ncont-Eptr->Ncont;
            if (Acont && strpbrk(Eptr->Cont->Type,"A")) k++;
            if(!strcmp(Eptr->Type,"RQ")) {
              Dval=Vj/Eptr->Tap*(dv1-dv2-dv3-dv4-dv6-dv8-dv10);
              if(Eptr->From==From) Dval=Dval+2*Vi/Eptr->Tap*(v5+v7+v9+v11);}
            else if(!strcmp(Eptr->Type,"RP"))
              Dval=Vj*(-v2+v1-v4-v3-v6-v8-v10);
            else Dval=0;
            JacElement(Jac,i+1+N,k,Dval);
          } else if(Rcont && Eptr->Cont!=NULL) {
            k=ACvar[Eptr->Cont->N]+1;
            if(!strcmp(Eptr->Type,"R")) {
              Dval=Vj/Eptr->Tap*(dv1-dv2-dv3-dv4-dv6-dv8-dv10);
              if(Eptr->From==From) Dval=Dval+2*Vi/Eptr->Tap*(v5+v7+v9+v11);}
            else Dval=0;
            JacElement(Jac,i+1+N,k,Dval);
          }
        } else {
          JacElement(Jac,i+1+N,j,0.);
          JacElement(Jac,i+1+N,j+1,0.);
          if (PQcont && strpbrk(Eptr->Type,"PQMN")){
            k=ACvar[Eptr->Cont->N]+1+Eptr->Cont->Ncont-Eptr->Ncont;
            if (Acont && strpbrk(Eptr->Cont->Type,"A")) k++;
            JacElement(Jac,i+1+N,k,0.);
          } else if(Rcont && Eptr->Cont!=NULL) {
            k=ACvar[Eptr->Cont->N]+1;
            JacElement(Jac,i+1+N,k,0.);
          }
        }
      }
    }
    if (flagF) {
      DPg=ACptr->DPG;
      if(DPg) {
        k=ACvar[BSptr->N];
        if (strpbrk(BSptr->Type,"S")) dF[k+N]+=DPg*x0[i];
        else if(Acont)                dF[k+N+2]+=DPg*x0[i];
      }
      if (ACptr->Cont!=NULL)
        dF[i+N+1]+=-2*Vi*(ACptr->G*x0[i]-ACptr->B*x0[i+1])
                   -2*Vi*(ACptr->Pz+ACptr->Pzl*lambda)*x0[i]
                   -a*pow(Vi,a-1)*(ACptr->Pn+ACptr->Pnl*lambda)*x0[i]
                   -2*Vi*(ACptr->Qz+ACptr->Qzl*lambda)*x0[i+1]
                   -b*pow(Vi,b-1)*(ACptr->Qn+ACptr->Qnl*lambda)*x0[i+1];
      if(QRcont && strpbrk(ACptr->Type,"G")) {
        j=ACvar[ACptr->Cont->N];
        if (strpbrk(ACptr->cont,"V")) dF[j+N+1]+=ACptr->Kbg*x0[i+1];
        else if (ACptr->Gen!=NULL) {
          j=ACptr->Gen->Nvar;
          if (strpbrk(ACptr->cont,"E"))   dF[j+N+1]+=x0[i+1];
          else                            dF[j+N+11]+=x0[i+1];
        }
      }
      else if(strpbrk(ACptr->Type,"V,Q,S,Z") || (!QRcont && strpbrk(ACptr->Type,"G")) ) {
        if (strpbrk(ACptr->cont,"V")) dF[i+N+1]+=x0[i+1];
        else if (ACptr->Gen!=NULL) {
          j=ACptr->Gen->Nvar;
          if (strpbrk(ACptr->cont,"E"))   dF[j+N+1]+=x0[i+1];
          else                            dF[j+N+11]+=x0[i+1];
        }
      }
    }
    if (flagJ){
      if (!strpbrk(From->Type,"S")) JacElement(Jac,i+N,i,Ddi);
      if (From->Cont!=NULL) {
        Dvip=Dvip-2*(ACptr->G*x0[i]-ACptr->B*x0[i+1])
             -2*(ACptr->Pz+ACptr->Pzl*lambda)*x0[i]
             -a*(a-1)*pow(Vi,a-2)*(ACptr->Pn+ACptr->Pnl*lambda)*x0[i]
             -2*(ACptr->Qz+ACptr->Qzl*lambda)*x0[i+1]
             -b*(b-1)*pow(Vi,b-2)*(ACptr->Qn+ACptr->Qnl*lambda)*x0[i+1];
        JacElement(Jac,i+N,i+1,Dvi);
        JacElement(Jac,i+1+N,i,Ddip);
        JacElement(Jac,i+1+N,i+1,Dvip);
        Dlambda=-(2*Vi*ACptr->Pzl+a*pow(Vi,a-1)*ACptr->Pnl)*x0[i]
                -(2*Vi*ACptr->Qzl+b*pow(Vi,b-1)*ACptr->Qnl)*x0[i+1];
        if (Dlambda) {
          JacElement(Jac,i+1+N,Jac->n1,Dlambda);