📄 jacobac.cpp
字号:
/* 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.);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -