writejac.cpp

用于潮流计算的程序请让我下载我需要的吧感谢了啊

    }

                               /* 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);
}