writejac.c

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

      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,TCSCptr->delta_t);
      sprintf(str,"Ftcsc%-d_%-d",k,++l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
  }
  for(k=0,STATCOMptr=dataPtr->STATCOMbus;STATCOMptr!=NULL;STATCOMptr=STATCOMptr->Next){
      k++; l=0;
      if(!strcmp(STATCOMptr->Cont,"PW") || !strcmp(STATCOMptr->Cont,"AL")){
        sprintf(str,"Istat%-d",k);i++;
        fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->I);
      } else {
	       sprintf(str,"Vrefc%-d",k);i++;
        fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->Vvar);
      }
      sprintf(str,"Vcont%-d_%-d",k,++l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"theta%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->theta);
      sprintf(str,"Cont%-d_%-d",k,++l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"Vdc%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->Vdc);
      sprintf(str,"Loss%-d_%-d",k,++l);  fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"k%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->k);
      sprintf(str,"ReS%-d_%-d",k,++l);   fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"alpha%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->alpha);
      sprintf(str,"ImS%-d_%-d",k,++l);   fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"Pstat%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->P);
      sprintf(str,"Pstat%-d_%-d",k,++l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      sprintf(str,"Qstat%-d",k);i++;
      fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,STATCOMptr->Q);
      sprintf(str,"Qstat%-d_%-d",k,++l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
  }
                             /* END FACTS */

  if (flagH) {
    fCustomPrint(OutFile,"%4d %8s %-11.5g\n",Jac->n1,"l",lambda);
    fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",Jac->n1,"dH",dF[Jac->n1]);
  }
  else if (flag) {
    N=NacVar+11*Ndc/2+3*Nsvc+NtcscVar+7*Nstatcom;      /* FACTS */
    for (i=N,ACptr=dataPtr->ACbus; ACptr!=NULL; ACptr=ACptr->Next){
      i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
      sprintf(str,"gP%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
      sprintf(str,"gQ%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      if(Acont && strpbrk(ACptr->Type,"A")){
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gPA%-d",ACptr->Area->N); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      }
      if (PQcont) for(ELptr=ACptr->Reg;ELptr!=NULL;ELptr=ELptr->Next) {
	       Eptr=ELptr->Eptr;
        if(strpbrk(Eptr->Type,"PQNM")) {
          if (Eptr->From==ACptr) {
            I=Eptr->From->Num;
            J=Eptr->To->Num;
          } else {
            J=Eptr->From->Num;
            I=Eptr->To->Num;
          }
          if(!strcmp(Eptr->Type,"PM")){
            i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
            sprintf(str,"gP%-d_%-d",I,J); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
          } else {
            i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
            sprintf(str,"gQ%-d_%-d",I,J); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
          }
        }
      }
      if (ACptr->Gen!=NULL) {
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gPg%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gQg%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gEq%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gEd%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gVd%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gVq%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gId%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gIq%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gVr%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gVi%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gIa%-d",ACptr->Num); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      }
    }
    for(k=0,DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next){
      if(!strcmp(DCptrR->Type,"R")){
        for (k++,l=1;l<=11;l++){
          i++; sprintf(str,"w%-d",i-N); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
          sprintf(str,"gdc%-d_%-d",k,l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
        }
      }
    }

                               /* FACTS */
    for(k=0,SVCptr=dataPtr->SVCbus;SVCptr!=NULL;SVCptr=SVCptr->Next){
      for (k++,l=1;l<=3;l++){
        i++; sprintf(str,"wsvc%-d_%-d",k,l); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gsvc%-d_%-d",k,l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      }
    }
    for(k=0,TCSCptr=dataPtr->TCSCbus;TCSCptr!=NULL;TCSCptr=TCSCptr->Next){
      for (k++,l=1;l<=7;l++){
        i++; sprintf(str,"wtcsc%-d_%-d",k,l); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gtcsc%-d_%-d",k,l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      }
    }
    for(k=0,STATCOMptr=dataPtr->STATCOMbus;STATCOMptr!=NULL;STATCOMptr=STATCOMptr->Next){
      for (k++,l=1;l<=7;l++){
        i++; sprintf(str,"wstat%-d_%-d",k,l); fCustomPrint(OutFile,"%4d %8s %-11.5g\n",i,str,x0[i-N]);
        sprintf(str,"gstat%-d_%-d",k,l); fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",i,str,dF[i]);
      }
    }
				                       /* END OF FACTS */

    fCustomPrint(OutFile,"%4d %8s %-11.5g\n",Jac->n1,"l",lambda);
    fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",Jac->n1,"gl",dF[Jac->n1]);
  }
  fCustomPrint(OutFile,"%4d %8s %-11.5g\n",0,"0",0.);
  fCustomPrint(OutFilep,"%4d %8s %-11.5g\n",0,"0",0.);
  fclose(OutFile);
  fclose(OutFilep);
}

/* --------------- WriteQmatrix ---------------------- */
#ifdef ANSIPROTO
void WriteQmatrix(INDEX count,VALUETYPE *vec)
#else
void WriteQmatrix(count,vec)
INDEX count;
VALUETYPE *vec;
#endif
{
  SparseMatrixElement *Jptr;
  ACbusData *ACptr;
  char Namebase[80],Name[80],MaxName[13];
  FILE *OutFile;
  INDEX i,j,I,J,m,MaxBus,MaxNum,Nvar;
  VALUETYPE val=0,valMax=-0.1;
  BOOLEAN flagTFbus=FALSE;

  strcpy(Namebase,NameParameter('0'));
  if(NullName(Namebase)) return;
  sprintf(Name,"%s%d.m",Namebase,count);
  OutFile=OpenOutput(Name);
  if (count==1) TFbus=IntegerParameter('1',0,1,9999);
  Nvar=NacVar+11*Ndc/2+3*Nsvc+NtcscVar+7*Nstatcom;     /* FACTS */
  for(i=1;i<=Nvar+1;i++) for(Jptr=Jac->RowHead[i];Jptr!=NULL;Jptr->Value=0,Jptr=Jptr->RowNext);
  ACFunJac(Jac,NULL,FALSE,TRUE,FALSE);
  DCFunJac(Jac,FALSE,TRUE);
  SVCFunJac(Jac,FALSE,TRUE);                  /*  FACTS  */
  TCSCFunJac(Jac,FALSE,TRUE);                 /*  FACTS  */
  STATCOMFunJac(Jac,FALSE,TRUE);              /*  FACTS  */
  fCustomPrint(OutFile,"lambda(%d)=%-10.6g;\n",count,lambda+lambda_o);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Full System Jacobian:\n","%%");
  fCustomPrint(OutFile,"N=%d;\nJ=zeros(N);\n",Nvar);
  for(i=1;i<=Nvar+1;i++) {
    for(Jptr=Jac->RowHead[i];Jptr!=NULL;Jptr=Jptr->RowNext) {
      I=Jptr->Row;
      if (OldRow->p[I]!=0) I=OldRow->p[I];
      J=Jptr->Col;
      if (OldCol->p[J]!=0) J=OldCol->p[J];
      if (I<=Nvar && J<=Nvar) fCustomPrint(OutFile,"J(%d,%d)=%-10.10g;\n",I,J,Jptr->Value);
    }
  }
  fCustomPrint(OutFile,"J=sparse(J);\n");
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Voltage (load) Buses:\n","%%");
  fCustomPrint(OutFile,"Vars=[(1:N)' zeros(N,1) zeros(N,1) ones(N,1)];\n");
  fCustomPrint(OutFile,"clear Buses Vbuses\n");
  for (m=0,ACptr=dataPtr->ACbus; ACptr!=NULL; ACptr=ACptr->Next){
    i=ACvar[ACptr->N];
    j=i+1;
    fCustomPrint(OutFile,"Vars(%d,2)=%d; Vars(%d,2)=%d;\n",i,-ACptr->Num,j,ACptr->Num);
    if (!strpbrk(ACptr->Type,"S")) fCustomPrint(OutFile,"Vars(%d,3)=1;\n",i);
    if ((ACptr->Cont!=NULL &&(QRcont || !strpbrk(ACptr->Type,"G"))) ||
       	(!Rcont && strpbrk(ACptr->Type,"T"))||
       	(!QRcont && strpbrk(ACptr->Type,"C"))){
      i++; m++;
      fCustomPrint(OutFile,"Vbuses(%d,1)=%d; Vbuses(%d,2)=%d;\n",m,i,m,ACptr->Num);
      fCustomPrint(OutFile,"Vars(%d,3)=1; Vars(%d,4)=0;\n",i,i);
      if (count==1) {
       	val=fabs(vec[i]);
       	if (TFbus==ACptr->Num && !flagTFbus) {
          flagTFbus=TRUE;
          TFnum=i;
          strcpy(TFname,ACptr->Name);
       	}
       	if (val>valMax) {
       	  valMax=val;
       	  MaxBus=ACptr->Num;
       	  MaxNum=i;
       	  strcpy(MaxName,ACptr->Name);
       	}
      }
    }
  }
  if (count==1 && !flagTFbus) {
    TFnum=MaxNum;
    TFbus=MaxBus;
    strcpy(TFname,MaxName);
  }
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Bus for test functions and tangent vector:\n","%%");
  fCustomPrint(OutFile,"l=%d;\n",TFbus);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Minimum singular value and real |e-value| for full J:\n","%%");
  fCustomPrint(OutFile,"[val,svecJ]=inviter(J'*J); svJ(%d)=sqrt(val);\n",count);
  fCustomPrint(OutFile,"[val,evecJ]=inviter(J); evJ(%d)=abs(val);\n",count);
  fCustomPrint(OutFile,"  %s Critical bus numbers and ranking of bus l:\n","%%");
  fCustomPrint(OutFile,"  [val,max_sv]=max(abs(svecJ));\n");
  fCustomPrint(OutFile,"  crsvJ(%d,1)=Vars(max_sv,2);\n",count);
  fCustomPrint(OutFile,"  [val,max_ev]=max(abs(evecJ));\n");
  fCustomPrint(OutFile,"  crevJ(%d,1)=Vars(max_ev,2);\n",count);
  fCustomPrint(OutFile,"  for i=1:N,\n");
  fCustomPrint(OutFile,"    if (Vars(i,2)==l), lnumJ=i; break; end\n");
  fCustomPrint(OutFile,"  end\n");
  fCustomPrint(OutFile,"  crsvJ(%d,2)=rankbus(svecJ,lnumJ);\n",count);
  fCustomPrint(OutFile,"  crevJ(%d,2)=rankbus(evecJ,lnumJ);\n",count);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Compute J_PF matrix:\n","%%");
  fCustomPrint(OutFile,"P0=zeros(N);\n");
  fCustomPrint(OutFile,"j=0;\n");
  fCustomPrint(OutFile,"for i=1:N,\n");
  fCustomPrint(OutFile,"  if(Vars(i,3)>0),\n");
  fCustomPrint(OutFile,"    j=j+1;\n");
  fCustomPrint(OutFile,"    Buses(j)=Vars(i,2);\n");
  fCustomPrint(OutFile,"    if (Vars(i,2)==l), lnumPF=j; end\n");
  fCustomPrint(OutFile,"    P0(Vars(i,1),j)=1;\n");
  fCustomPrint(OutFile,"  end\n");
  fCustomPrint(OutFile,"end\n");
  fCustomPrint(OutFile,"M=j;\n");
  fCustomPrint(OutFile,"for i=1:N,\n");
  fCustomPrint(OutFile,"  if(Vars(i,3)==0), j=j+1; P0(Vars(i,1),j)=1; end\n");
  fCustomPrint(OutFile,"end\n");
  fCustomPrint(OutFile,"P0=sparse(P0); Jp=P0'*J*P0;\n");
  fCustomPrint(OutFile,"J_PF=Jp(1:M,1:M);\n");
  fCustomPrint(OutFile,"%s Minimum singular value and real |e-value| for standard P.F. J_PF:\n","%%");
  fCustomPrint(OutFile,"[val,svecPF]=inviter(J_PF'*J_PF); svPF(%d)=sqrt(val);\n",count);
  fCustomPrint(OutFile,"[val,evecPF]=inviter(J_PF); evPF(%d)=abs(val);\n",count);
  fCustomPrint(OutFile,"  %s Critical bus numbers and ranking of bus l:\n","%%");
  fCustomPrint(OutFile,"  [val,max_sv]=max(abs(svecPF));\n");
  fCustomPrint(OutFile,"  crsvPF(%d,1)=Buses(max_sv);\n",count);
  fCustomPrint(OutFile,"  [val,max_ev]=max(abs(evecPF));\n");
  fCustomPrint(OutFile,"  crevPF(%d,1)=Buses(max_ev);\n",count);
  fCustomPrint(OutFile,"  crsvPF(%d,2)=rankbus(svecPF,lnumPF);\n",count);
  fCustomPrint(OutFile,"  crevPF(%d,2)=rankbus(evecPF,lnumPF);\n",count);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Compute J_QV matrix:\n","%%");
  fCustomPrint(OutFile,"P1=zeros(N);\n");
  fCustomPrint(OutFile,"Mp=%d; j=Mp;\n",m);
  fCustomPrint(OutFile,"for i=1:Mp,\n");
  fCustomPrint(OutFile,"  if (Vbuses(i,2)==l), lnumQV=i; end\n");
  fCustomPrint(OutFile,"  P1(Vbuses(i,1),i)=1;\n");
  fCustomPrint(OutFile,"end;\n");
  fCustomPrint(OutFile,"for i=1:N,\n");
  fCustomPrint(OutFile,"  if(Vars(i,4)>0), j=j+1; P1(Vars(i,1),j)=1; end\n");
  fCustomPrint(OutFile,"end\n");
  fCustomPrint(OutFile,"P1=sparse(P1); Jp=P1'*J*P1;\n");
  fCustomPrint(OutFile,"A1=Jp(1:Mp,1:Mp);  B1=Jp(1:Mp,Mp+1:N); C1=Jp(Mp+1:N,1:Mp); D1=Jp(Mp+1:N,Mp+1:N);\n");
  fCustomPrint(OutFile,"J_QV=A1-B1*(D1\\C1);\n");
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Minimum singular value and real |e-value| for J_QV:\n","%%");
  fCustomPrint(OutFile,"[val,svecQV]=inviter(J_QV'*J_QV); svQV(%d)=sqrt(val);\n",count);
  fCustomPrint(OutFile,"[val,evecQV]=inviter(J_QV,-0.001); evQV(%d)=abs(val);\n",count);
  fCustomPrint(OutFile,"  %s Critical V bus numbers and ranking of bus l:\n","%%");
  fCustomPrint(OutFile,"  [val,max_sv]=max(abs(svecQV));\n");
  fCustomPrint(OutFile,"  crsvQV(%d,1)=Vbuses(max_sv,2);\n",count);
  fCustomPrint(OutFile,"  [val,max_ev]=max(abs(evecQV));\n");
  fCustomPrint(OutFile,"  crevQV(%d,1)=Vbuses(max_ev,2);\n",count);
  fCustomPrint(OutFile,"  crsvQV(%d,2)=rankbus(svecQV,lnumQV);\n",count);
  fCustomPrint(OutFile,"  crevQV(%d,2)=rankbus(evecQV,lnumQV);\n",count);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Reduced determinant detD_ll, for bus l=%d, '%12s':\n","%%",TFbus,TFname);
  fCustomPrint(OutFile,"P2=speye(N);\n");
  fCustomPrint(OutFile,"v=P2(1,:); P2(1,:)=P2(%d,:); P2(%d,:)=v; \n",TFnum-1,TFnum-1);
  fCustomPrint(OutFile,"v=P2(2,:); P2(2,:)=P2(%d,:); P2(%d,:)=v; \n",TFnum,TFnum);
  fCustomPrint(OutFile,"Jp=P2'*J*P2;\n");
  fCustomPrint(OutFile,"A2=Jp(1:2,1:2);  B2=Jp(1:2,3:N); C2=Jp(3:N,1:2); D2=Jp(3:N,3:N);\n");
  fCustomPrint(OutFile,"detD_ll(%d)=det(A2-B2*(D2\\C2));\n",count);
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"%s Test Function t_ll, for bus l=%d, '%12s':\n","%%",TFbus,TFname);
  fCustomPrint(OutFile,"el=zeros(N,1); el(%d)=1; el=sparse(el);\n",TFnum);
  fCustomPrint(OutFile,"J_ll=(speye(N)-el*el')*J+el*el';\n");
  fCustomPrint(OutFile,"t_ll(%d)=el'*J*(J_ll\\el);\n",count);
  fclose(OutFile);
}