jacobac.c

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

/* AC Function and Jacobian. */

#include "jacob.h"

/* ------------------ ACFunJac ----------------------------- */
#ifdef ANSIPROTO
AreaData *ACFunJac(SparseMatrix *Mptr,int *val,BOOLEAN flagF,BOOLEAN flagJ,BOOLEAN flagP)
#else
AreaData *ACFunJac(Mptr,val,flagF,flagJ,flagP)
SparseMatrix *Mptr;
int *val;
BOOLEAN flagF,flagJ,flagP;
#endif
/* Construct the AC part of the Jacobian. */
{
  ACbusData *ACptr,*To,*From,*BEptr;
  AClist *ALptr;
  AreaData *Aptr;
  ElementList *ELptr;
  ElementData *Eptr;
  VALUETYPE Pi,Qi,Vi,Vj,di,dj,SPij,dPijd,dPiid,dQijd,dQiid,val1;
  VALUETYPE gij,bij,gsij,bsij,SQij,dPijv,dPiiv,dQijv,dQiiv,val2;
  VALUETYPE Ra,Xd,Xq,Qg,Eq,dg,Vr,Vim,Ir,Iim,Vq,Vd,Iq,Id,Ia;
  VALUETYPE Pl,Ql,Pg,DPg,SPg=0;
  INDEX i,j,k,l;

  if (val!=NULL) *val=0;
  for(ALptr=dataPtr->KGbus;ALptr!=NULL;ALptr=ALptr->Next) {
    BEptr=ALptr->AC;
    if(strpbrk(ALptr->AC->Type,"S"))  break;
  }
  if (Acont) for(Aptr=dataPtr->Area;Aptr!=NULL;Aptr->SPg=0,Aptr=Aptr->Next);
  if (Bl) l=ACvar[Bl]+1;
  else if (flagH)  l=Mptr->n1;
  else l=0;
  for (ACptr=dataPtr->ACbus; ACptr!=NULL; ACptr=ACptr->Next){
    if (ACptr->Area!=NULL) BEptr=ACptr->Area->Slack;
    i=ACvar[ACptr->N];

  /* ------------------- Power Mismatches ------------------- */
    Vi=ACptr->V;
    di=ACptr->Ang;
    if (flagP) {Pl=ACptr->Pl; Ql=ACptr->Ql;}
    else {
      Pl=(ACptr->Pn+lambda*ACptr->Pnl)*pow(Vi,ACptr->a)+
         (ACptr->Pz+lambda*ACptr->Pzl)*Vi*Vi;
      Ql=(ACptr->Qn+lambda*ACptr->Qnl)*pow(Vi,ACptr->b)+
         (ACptr->Qz+lambda*ACptr->Qzl)*Vi*Vi;
    }
    Pg=ACptr->PG;
    DPg=ACptr->DPG;
    if (Acont && Narea>1) ACptr->Area->Slack->Area->SPg+=DPg;
    else SPg+=DPg;
    Pi=Pg-Pl-Vi*Vi*ACptr->G;
    Qg=ACptr->Qg;
    Qi=Qg-Ql+Vi*Vi*ACptr->B;
    dPiid=dQiid=0;
    dPiiv=dQiiv=0;
    SPij=SQij=0;
    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==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)*Eptr->Tap*Eptr->Tap-gij;
        bsij=(Eptr->B1+Eptr->B)*Eptr->Tap*Eptr->Tap-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;
      }
      dPijd=Vi*Vj*(gij*sin(di-dj)-bij*cos(di-dj));
      dPiid=dPiid-dPijd;
      dPijv=Vi*(gij*cos(di-dj)+bij*sin(di-dj));
      dPiiv=dPiiv+Vj/Vi*dPijv-2*Vi*(gij+gsij);
      dQijd= -Vj*dPijv;
      dQiid=dQiid-dQijd;
      dQijv=dPijd/Vj;
      dQiiv=dQiiv+Vj/Vi*dQijv+2*Vi*(bij+bsij);
      SPij=SPij+Vi*Vi*(gij+gsij)+dQijd;
      SQij=SQij-Vi*Vi*(bij+bsij)-dPijd;
      if(flagJ){
        /* df/ddelta_j */
        if (!strpbrk(To->Type,"S")) {
          JacElement(Mptr,i,j,dPijd);
          JacElement(Mptr,i+1,j,dQijd);
        }
        /* df/dV_j */
        if (To->Cont!=NULL) {
          JacElement(Mptr,i,j+1,dPijv);
          JacElement(Mptr,i+1,j+1,dQijv);
        } else if (flagH && strpbrk(To->Type,"L")) {
          JacElement(Mptr,i,Mptr->n1,dPijv);
          JacElement(Mptr,i+1,Mptr->n1,dQijv);
        } else {
          JacElement(Mptr,i,j+1,0.);
          JacElement(Mptr,i+1,j+1,0.);
        }
        /* df/dControl_ij */
        if (PQcont && strpbrk(Eptr->Type,"PQMN")){
          if (Acont && strpbrk(Eptr->Cont->Type,"A")) k=1; else k=0;
          j=ACvar[Eptr->Cont->N]+1+k+Eptr->Cont->Ncont-Eptr->Ncont;
          if(!strcmp(Eptr->Type,"RQ")) {
            val1= -dQijd/Eptr->Tap;
            val2=dPijd/Eptr->Tap;
            if (Eptr->From==ACptr){
              val1=val1-2*Vi*Vi*(Eptr->G1+Eptr->G)*Eptr->Tap;
              val2=val2+2*Vi*Vi*(Eptr->B1+Eptr->B)*Eptr->Tap;
            }
          }
          else if(!strcmp(Eptr->Type,"RP")) {
            val1=dPijd;
            val2=dQijd;
            if (Eptr->From!=ACptr) {
              val1= -val1;
              val2= -val2;
            }
          } else val1=val2=0;
          JacElement(Mptr,i,j,val1);
          JacElement(Mptr,i+1,j,val2);
        }
        else if(Rcont && Eptr->Cont!=NULL) {
          j=ACvar[Eptr->Cont->N]+1;
          if(!strcmp(Eptr->Type,"R")) {
            val1= -dQijd/Eptr->Tap;
            val2=dPijd/Eptr->Tap;
            if (Eptr->From==ACptr){
              val1=val1-2*Vi*Vi*(Eptr->G1+Eptr->G)*Eptr->Tap;
              val2=val2+2*Vi*Vi*(Eptr->B1+Eptr->B)*Eptr->Tap;
            }
          } else val1=val2=0;
          JacElement(Mptr,i,j,val1);
          JacElement(Mptr,i+1,j,val2);
        }
      }
    }
    if (flagF) {
      dF[i]=Pi-SPij;
      dF[i+1]=Qi-SQij;
    }
    if (flagJ) {
      dPiiv=dPiiv-2*Vi*ACptr->G;
      dQiiv=dQiiv+2*Vi*ACptr->B;
      if (!flagP) {
        dPiiv=dPiiv-2*(ACptr->Pz+ACptr->Pzl*lambda)*Vi-
              ACptr->a*(ACptr->Pn+ACptr->Pnl*lambda)*pow(Vi,ACptr->a-1.0);
        dQiiv=dQiiv-2*(ACptr->Qz+ACptr->Qzl*lambda)*Vi-
              ACptr->b*(ACptr->Qn+ACptr->Qnl*lambda)*pow(Vi,ACptr->b-1.0);
      }
      /* df/dKg */
      j=ACvar[BEptr->N];
      if(DPg) {
        if (strpbrk(BEptr->Type,"S")) { JacElement(Mptr,i,j,DPg);}
        else if(Acont) { JacElement(Mptr,i,j+2,DPg);}
      }
      /* df/dlambda */
      if (l) {
        if (ACptr->Pzl || ACptr->Pnl) JacElement(Mptr,i,l,-ACptr->Pzl*Vi*Vi-ACptr->Pnl*pow(Vi,ACptr->a));
        if (ACptr->Qzl || ACptr->Qnl) JacElement(Mptr,i+1,l,-ACptr->Qzl*Vi*Vi-ACptr->Qnl*pow(Vi,ACptr->b));
      }
      /* df/ddelta_i */
      if (!strpbrk(ACptr->Type,"S")) {
        JacElement(Mptr,i,i,dPiid);
        JacElement(Mptr,i+1,i,dQiid);
      }
      /* df/dV_i */
      if (ACptr->Cont!=NULL) {
        JacElement(Mptr,i,i+1,dPiiv);
        JacElement(Mptr,i+1,i+1,dQiiv);
      } else if (flagH && strpbrk(ACptr->Type,"L")) {
        JacElement(Mptr,i,Mptr->n1,dPiiv);
        JacElement(Mptr,i+1,Mptr->n1,dQiiv);
      } else {
        JacElement(Mptr,i,i+1,0.);
        JacElement(Mptr,i+1,i+1,0.);
      }
      /*  df/dQg   */
      if(strpbrk(ACptr->Type,"V,Q,S,G,Z")) {
        if(QRcont && strpbrk(ACptr->Type,"G")) {
          j=ACvar[ACptr->Cont->N];
          val1=ACptr->Kbg;
        } else {
          j=i;
          val1=1.;
        }
        if (strpbrk(ACptr->cont,"V")) JacElement(Mptr,i+1,j+1,val1);
        else                          JacElement(Mptr,i+1,j+1,0.);
        if (ACptr->Gen!=NULL) {  /*Gen. Model BEGIN */
          j=ACptr->Gen->Nvar;
          if (strpbrk(ACptr->cont,"I")) JacElement(Mptr,i+1,j+11,1.);
          else                          JacElement(Mptr,i+1,j+11,0.);
          if (strpbrk(ACptr->cont,"E")) JacElement(Mptr,i+1,j+1,1.);
          else                          JacElement(Mptr,i+1,j+1,0.);
        } /* Gen. Model  END */
      }
    }

  /* ---------------------- Area Control --------------------- */
    if (Acont && strpbrk(ACptr->Type,"A")){
      i=i+2;
      SPij=0;
      for (ELptr=BEptr->Area->Elem; ELptr!=NULL; ELptr=ELptr->Next){
        Eptr=ELptr->Eptr;
        From=Eptr->From;
        To=Eptr->To;
        if (Eptr->From!=Eptr->Meter) {
          From=Eptr->To;
          To=Eptr->From;
        }
        val2=1.0;
        if (Eptr->Meter->Area!=ACptr->Area) val2= -1.0;
        Vi=From->V;
        di=From->Ang;
        Vj=To->V;
        dj=To->Ang;
        if(Eptr->From==Eptr->Meter) {
          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;
        } 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;
        }
        dPiid= val2*(-Vi*Vj*(gij*sin(di-dj)-bij*cos(di-dj)));
        dPiiv= val2*(-2*Vi*(gij+gsij)+Vj*(gij*cos(di-dj)+bij*sin(di-dj)));
        dPijd= -dPiid;
        dPijv=val2*(Vi*(gij*cos(di-dj)+bij*sin(di-dj)));
        SPij=SPij+val2*Vi*Vi*(gij+gsij)-Vj*dPijv;
        if(flagJ){
          j=ACvar[From->N];
          if (!strpbrk(From->Type,"S")) {  JacElement(Mptr,i,j,dPiid);}
          if (From->Cont!=NULL) JacElement(Mptr,i,j+1,dPiiv);
          else if (flagH && strpbrk(From->Type,"L")) JacElement(Mptr,i,Mptr->n1,dPiiv);
          else JacElement(Mptr,i,j+1,0.);
          j=ACvar[To->N];
          if (!strpbrk(To->Type,"S")) { JacElement(Mptr,i,j,dPijd);}
          if (To->Cont!=NULL) JacElement(Mptr,i,j+1,dPijv);
          else if (flagH && strpbrk(To->Type,"L")) JacElement(Mptr,i,Mptr->n1,dPijv);
          else JacElement(Mptr,i,j+1,0.);