📄 pointdc.cpp
字号:
/* Point of Collapse: DC Hessian. */
#include "pointl.h"
/* ------------------ DCFunHes ----------------------------- */
#ifdef ANSIPROTO
BOOLEAN DCFunHes(BOOLEAN flagF,BOOLEAN flagJ)
#else
BOOLEAN DCFunHes(flagF,flagJ)
BOOLEAN flagF,flagJ;
#endif
/* Construct the DC part of the PoC Jacobian. */
{
INDEX i,j,k,l,kp,lp,m,n,N,DCvar[16];
DCbusData *DCptrR,*DCptrI;
ACbusData *BEptr;
VALUETYPE Sa1,Sa2;
VALUETYPE Vdr,Vr,ar,cosar,cosgr,Xcr,Sr,Pr,Dr,Id;
VALUETYPE Vdi,Vi,ai,cosai,cosgi,Xci,Si,Pi,Di,Rd;
BOOLEAN dVr=FALSE,dVi=FALSE;
i=NacVar; N=NacVar+11*Ndc/2+3*Nsvc+NtcscVar+7*Nstatcom; /* FACTS */
for(DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next){
DCptrI=DCptrR->To;
if (!strcmp(DCptrR->Type,"R")){
Id=DCptrR->Id; Rd=DCptrR->Rd;
Vdr=DCptrR->Vd; Vdi=DCptrI->Vd;
k=ACvar[DCptrR->AC->N]; l=ACvar[DCptrI->AC->N];
kp=ACvar[DCptrR->AC->N]+1; lp=ACvar[DCptrI->AC->N]+1;
if (DCptrR->AC->Cont!=NULL) dVr=TRUE;
if (DCptrI->AC->Cont!=NULL) dVi=TRUE;
Vr=DCptrR->AC->V; Vi=DCptrI->AC->V;
ar=DCptrR->Tap*DCptrR->Ntrf; ai=DCptrI->Tap*DCptrI->Ntrf;
Xcr=DCptrR->Xc; Xci=DCptrI->Xc;
cosar=cos(DCptrR->Alfa); cosai=cos(DCptrI->Alfa);
cosgr=cos(DCptrR->Gamma); cosgi=cos(DCptrI->Gamma);
Pr=DCptrR->P; Pi=DCptrI->P;
Sr=DCptrR->MVA; Si=DCptrI->MVA;
Dr=Sr*Sr-Pr*Pr; Di=Si*Si-Pi*Pi;
if (Dr<=0 || Di<=0) return(TRUE);
if (Acont && DCptrR->Meter!=NULL) {
if(DCptrR->Meter==DCptrR) { Sa1= -1; Sa2=1;}
else { Sa1=1; Sa2= -1;}
BEptr=DCptrR->Area->Slack;
if (!strpbrk(BEptr->Type,"S")) m=ACvar[BEptr->N]+2;
else m=0;
BEptr=DCptrI->Area->Slack;
if (!strpbrk(BEptr->Type,"S")) n=ACvar[BEptr->N]+2;
else n=0;
} else m=n=0;
if (flagF) {
j=i+N;
if (dVr) dF[kp+N]=dF[kp+N]-K1*ar*cosar*x0[i+1]-K1*ar*Id*x0[i+2]
-ar*(cosar+cosgr)*x0[i+9];
if (strcmp(DCptrR->Cont1,"VD") && strcmp(DCptrR->Cont2,"VD"))
dF[++j]=x0[i+1]+Id*x0[i+3]+x0[i+11];
if (strcmp(DCptrR->Cont1,"AT") && strcmp(DCptrR->Cont2,"AT"))
dF[++j]= -K1*Vr*cosar*x0[i+1]-K1*Vr*Id*x0[i+2]-Vr*(cosar+cosgr)*x0[i+9];
if (strcmp(DCptrR->Cont1,"AL") && strcmp(DCptrR->Cont2,"AL"))
dF[++j]= -K1*ar*Vr*x0[i+1]-ar*Vr*x0[i+9];
if (strcmp(DCptrR->Cont1,"GA") && strcmp(DCptrR->Cont2,"GA"))
dF[++j]= -ar*Vr*x0[i+9];
dF[++j]=x0[i+2]+Sr/sqrt(Dr)*x0[i+4];
if (strcmp(DCptrR->Cont1,"PA") && strcmp(DCptrR->Cont2,"PA")) {
dF[++j]=x0[i+3]-Pr/sqrt(Dr)*x0[i+4]+x0[k];
if (m!=0 && DCptrR->Meter==DCptrR) dF[j]=dF[j]-Sa1*x0[m];
if (n!=0 && DCptrR->Meter==DCptrR) dF[j]=dF[j]-Sa2*x0[n];
}
if (strcmp(DCptrR->Cont1,"QA") && strcmp(DCptrR->Cont2,"QA"))
dF[++j]=x0[i+4]+x0[kp];
if (dVi) dF[lp+N]=dF[lp+N]-K1*ai*cosgi*x0[i+5]-K1*ai*Id*x0[i+6]
-ai*(cosai+cosgi)*x0[i+10];
if (strcmp(DCptrI->Cont1,"VD") && strcmp(DCptrI->Cont2,"VD"))
dF[++j]=x0[i+5]-Id*x0[i+7]-x0[i+11];
if (strcmp(DCptrI->Cont1,"AT") && strcmp(DCptrI->Cont2,"AT"))
dF[++j]= -K1*Vi*cosgi*x0[i+5]-K1*Vi*Id*x0[i+6]-Vi*(cosai+cosgi)*x0[i+10];
if (strcmp(DCptrI->Cont1,"AL") && strcmp(DCptrI->Cont2,"AL"))
dF[++j]= -ai*Vi*x0[i+10];
if (strcmp(DCptrI->Cont1,"GA") && strcmp(DCptrI->Cont2,"GA"))
dF[++j]= -K1*ai*Vi*x0[i+5]-ai*Vi*x0[i+10];
dF[++j]=x0[i+6]+Si/sqrt(Di)*x0[i+8];
if (strcmp(DCptrI->Cont1,"PA") && strcmp(DCptrI->Cont2,"PA")) {
dF[++j]=x0[i+7]-Pi/sqrt(Di)*x0[i+8]+x0[l];
if (m!=0 && DCptrI->Meter==DCptrI) dF[j]=dF[j]-Sa1*x0[m];
if (n!=0 && DCptrI->Meter==DCptrI) dF[j]=dF[j]-Sa2*x0[n];
}
if (strcmp(DCptrI->Cont1,"QA") && strcmp(DCptrI->Cont2,"QA"))
dF[++j]=x0[i+8]+x0[lp];
if (strcmp(DCptrR->Cont1,"ID") && strcmp(DCptrR->Cont2,"ID") &&
strcmp(DCptrI->Cont1,"ID") && strcmp(DCptrI->Cont2,"ID"))
dF[++j]=K2*Xcr*x0[i+1]-K1*ar*Vr*x0[i+2]+Vdr*x0[i+3]+K2*Xci*x0[i+5]
-K1*ai*Vi*x0[i+6]-Vdi*x0[i+7]+sqrt(2.0)*(Xcr*x0[i+9]+Xci*x0[i+10])
-Rd*x0[i+11];
}
if (flagJ) {
j=i;
if (strcmp(DCptrR->Cont1,"VD") && strcmp(DCptrR->Cont2,"VD")) DCvar[1]= ++j;
else DCvar[1]=0;
if (strcmp(DCptrR->Cont1,"AT") && strcmp(DCptrR->Cont2,"AT")) DCvar[2]= ++j;
else DCvar[2]=0;
if (strcmp(DCptrR->Cont1,"AL") && strcmp(DCptrR->Cont2,"AL")) DCvar[3]= ++j;
else DCvar[3]=0;
if (strcmp(DCptrR->Cont1,"GA") && strcmp(DCptrR->Cont2,"GA")) DCvar[4]= ++j;
else DCvar[4]=0;
DCvar[5]= ++j;
if (strcmp(DCptrR->Cont1,"PA") && strcmp(DCptrR->Cont2,"PA")) DCvar[6]= ++j;
else DCvar[6]=0;
if (strcmp(DCptrR->Cont1,"QA") && strcmp(DCptrR->Cont2,"QA")) DCvar[7]= ++j;
else DCvar[7]=0;
if (strcmp(DCptrI->Cont1,"VD") && strcmp(DCptrI->Cont2,"VD")) DCvar[8]= ++j;
else DCvar[8]=0;
if (strcmp(DCptrI->Cont1,"AT") && strcmp(DCptrI->Cont2,"AT")) DCvar[9]= ++j;
else DCvar[9]=0;
if (strcmp(DCptrI->Cont1,"AL") && strcmp(DCptrI->Cont2,"AL")) DCvar[10]= ++j;
else DCvar[10]=0;
if (strcmp(DCptrI->Cont1,"GA") && strcmp(DCptrI->Cont2,"GA")) DCvar[11]= ++j;
else DCvar[11]=0;
DCvar[12]= ++j;
if (strcmp(DCptrI->Cont1,"PA") && strcmp(DCptrI->Cont2,"PA")) DCvar[13]= ++j;
else DCvar[13]=0;
if (strcmp(DCptrI->Cont1,"QA") && strcmp(DCptrI->Cont2,"QA")) DCvar[14]= ++j;
else DCvar[14]=0;
if (strcmp(DCptrR->Cont1,"ID") && strcmp(DCptrR->Cont2,"ID") &&
strcmp(DCptrI->Cont1,"ID") && strcmp(DCptrI->Cont2,"ID")) DCvar[15]= ++j;
else DCvar[15]=0;
if (dVr){
if(DCvar[2]) JacElement(Jac,kp+N,DCvar[2],-K1*cosar*x0[i+1]-K1*Id*x0[i+2]-(cosar+cosgr)*x0[i+9]);
if(DCvar[3]) JacElement(Jac,kp+N,DCvar[3],-K1*ar*x0[i+1]-ar*x0[i+9]);
if(DCvar[4]) JacElement(Jac,kp+N,DCvar[4],-ar*x0[i+9]);
if(DCvar[15]) JacElement(Jac,kp+N,DCvar[15],-K1*ar*x0[i+2]);
} else {
if(DCvar[2]) JacElement(Jac,kp+N,DCvar[2],0.);
if(DCvar[3]) JacElement(Jac,kp+N,DCvar[3],0.);
if(DCvar[4]) JacElement(Jac,kp+N,DCvar[4],0.);
if(DCvar[15]) JacElement(Jac,kp+N,DCvar[15],0.);
}
if (strcmp(DCptrR->Cont1,"VD") && strcmp(DCptrR->Cont2,"VD")) {
if(DCvar[15]) JacElement(Jac,DCvar[1]+N,DCvar[15],x0[i+3]);
}
if (strcmp(DCptrR->Cont1,"AT") && strcmp(DCptrR->Cont2,"AT")){
if(dVr) JacElement(Jac,DCvar[2]+N,kp,-K1*cosar*x0[1+1]-K1*Id*x0[i+2]-(cosar*cosgr)*x0[i+9]);
else JacElement(Jac,DCvar[2]+N,kp,0.);
if(DCvar[3]) JacElement(Jac,DCvar[2]+N,DCvar[3],-K1*Vr*x0[i+1]-Vr*x0[i+9]);
if(DCvar[4]) JacElement(Jac,DCvar[2]+N,DCvar[4],-Vr*x0[i+9]);
if(DCvar[15]) JacElement(Jac,DCvar[2]+N,DCvar[15],-K1*Vr*x0[i+2]);
}
if (strcmp(DCptrR->Cont1,"AL") && strcmp(DCptrR->Cont2,"AL")){
if(dVr) JacElement(Jac,DCvar[3]+N,kp,-K1*ar*x0[1+1]-ar*x0[i+9]);
else JacElement(Jac,DCvar[3]+N,kp,0.);
if(DCvar[2]) JacElement(Jac,DCvar[3]+N,DCvar[2],-K1*Vr*x0[1+1]-Vr*x0[i+9]);
}
if (strcmp(DCptrR->Cont1,"GA") && strcmp(DCptrR->Cont2,"GA")){
if(dVr) JacElement(Jac,DCvar[4]+N,kp,-ar*x0[i+9]);
else JacElement(Jac,DCvar[4]+N,kp,0.);
if(DCvar[2]) JacElement(Jac,DCvar[4]+N,DCvar[2],-Vr*x0[i+9]);
}
JacElement(Jac,DCvar[5]+N,DCvar[5],-Pr*Pr*x0[i+4]/(Dr*sqrt(Dr)));
if(DCvar[6]) JacElement(Jac,DCvar[5]+N,DCvar[6],Sr*Pr*x0[i+4]/(Dr*sqrt(Dr)));
if (strcmp(DCptrR->Cont1,"PA") && strcmp(DCptrR->Cont2,"PA")){
JacElement(Jac,DCvar[6]+N,DCvar[5],Sr*Pr*x0[i+4]/(Dr*sqrt(Dr)));
JacElement(Jac,DCvar[6]+N,DCvar[6],-Sr*Sr*x0[i+4]/(Dr*sqrt(Dr)));
}
if (dVi){
if(DCvar[9]) JacElement(Jac,lp+N,DCvar[9],-K1*cosgi*x0[i+5]-K1*Id*x0[i+6]-(cosai+cosgi)*x0[i+10]);
if(DCvar[10]) JacElement(Jac,lp+N,DCvar[10],-ai*x0[i+10]);
if(DCvar[11]) JacElement(Jac,lp+N,DCvar[11],-K1*ai*x0[i+5]-ai*x0[i+10]);
if(DCvar[15]) JacElement(Jac,lp+N,DCvar[15],-K1*ai*x0[i+6]);
} else {
if(DCvar[9]) JacElement(Jac,lp+N,DCvar[9],0.);
if(DCvar[10]) JacElement(Jac,lp+N,DCvar[10],0.);
if(DCvar[11]) JacElement(Jac,lp+N,DCvar[11],0.);
if(DCvar[15]) JacElement(Jac,lp+N,DCvar[15],0.);
}
if (strcmp(DCptrI->Cont1,"VD") && strcmp(DCptrI->Cont2,"VD")) {
if(DCvar[15]) JacElement(Jac,DCvar[8]+N,DCvar[15],-x0[i+7]);
}
if (strcmp(DCptrI->Cont1,"AT") && strcmp(DCptrI->Cont2,"AT")){
if(dVi) JacElement(Jac,DCvar[9]+N,lp,-K1*cosgi*x0[i+5]-K1*Id*x0[i+6]-(cosai+cosgi)*x0[i+10]);
else JacElement(Jac,DCvar[9]+N,lp,0.);
if(DCvar[10]) JacElement(Jac,DCvar[9]+N,DCvar[10],-Vi*x0[i+10]);
if(DCvar[11]) JacElement(Jac,DCvar[9]+N,DCvar[11],-K1*Vi*x0[i+5]-Vi*x0[i+10]);
if(DCvar[15]) JacElement(Jac,DCvar[9]+N,DCvar[15],-K1*Vi*x0[i+6]);
}
if (strcmp(DCptrI->Cont1,"AL") && strcmp(DCptrI->Cont2,"AL")){
if(dVi) JacElement(Jac,DCvar[10]+N,lp,-ai*x0[i+10]);
else JacElement(Jac,DCvar[10]+N,lp,0.);
if(DCvar[9]) JacElement(Jac,DCvar[10]+N,DCvar[9],-Vi*x0[i+10]);
}
if (strcmp(DCptrI->Cont1,"GA") && strcmp(DCptrI->Cont2,"GA")){
if(dVi) JacElement(Jac,DCvar[11]+N,lp,-K1*ai*x0[i+5]-ai*x0[i+10]);
else JacElement(Jac,DCvar[11]+N,lp,0.);
if(DCvar[9]) JacElement(Jac,DCvar[11]+N,DCvar[9],-K1*Vi*x0[i+5]-Vi*x0[i+10]);
}
JacElement(Jac,DCvar[12]+N,DCvar[12],-Pi*Pi*x0[i+8]/(Di*sqrt(Di)));
if(DCvar[13]) JacElement(Jac,DCvar[12]+N,DCvar[13],Si*Pi*x0[i+8]/(Di*sqrt(Di)));
if (strcmp(DCptrI->Cont1,"PA") && strcmp(DCptrI->Cont2,"PA")){
JacElement(Jac,DCvar[13]+N,DCvar[12],Si*Pi*x0[i+8]/(Di*sqrt(Di)));
JacElement(Jac,DCvar[13]+N,DCvar[13],-Si*Si*x0[i+8]/(Di*sqrt(Di)));
}
if (strcmp(DCptrR->Cont1,"ID") && strcmp(DCptrR->Cont2,"ID") &&
strcmp(DCptrI->Cont1,"ID") && strcmp(DCptrI->Cont2,"ID")){
if (dVr) JacElement(Jac,DCvar[15]+N,kp,-K1*ar*x0[i+2]);
else JacElement(Jac,DCvar[15]+N,kp,0.);
if(DCvar[1]) JacElement(Jac,DCvar[15]+N,DCvar[1],x0[i+3]);
if(DCvar[2]) JacElement(Jac,DCvar[15]+N,DCvar[2],-K1*Vr*x0[i+2]);
if (dVi) JacElement(Jac,DCvar[15]+N,lp,-K1*ai*x0[i+6]);
else JacElement(Jac,DCvar[15]+N,lp,0.);
if(DCvar[8]) JacElement(Jac,DCvar[15]+N,DCvar[8],-x0[i+7]);
if(DCvar[9]) JacElement(Jac,DCvar[15]+N,DCvar[9],-K1*Vi*x0[i+6]);
}
}
i=i+11;
}
}
return(FALSE);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -