homotpro.c

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

      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    else if(strpbrk(ACptr->Type,"Z")) {
      sprintf(str,"d%-d",ACptr->Num); i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Qz%-d",ACptr->Num); i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    else if(strpbrk(ACptr->Type,"S")){
      sprintf(str,"kg%-d",ACptr->Num); i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Qg%-d",ACptr->Num); i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    if(Acont && strpbrk(ACptr->Type,"A")){
      sprintf(str,"kg%-d",ACptr->Num); i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    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,"RP")){
           sprintf(str,"a%-d_%-d",I,J); i++;
           fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
         } else if(strpbrk(Eptr->Type,"PM")){
           sprintf(str,"P%-d_%-d",I,J); i++;
           fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
         } else if(!strcmp(Eptr->Type,"RQ")){
           sprintf(str,"1/t%-d_%-d",I,J); i++;
           fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
         } else {
           sprintf(str,"Q%-d_%-d",I,J); i++;
           fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
         }
      }
    }
    if (ACptr->Gen!=NULL) {
      i=ACptr->Gen->Nvar;
      if (!strpbrk(ACptr->cont,"E")) sprintf(str,"Eq%-d",ACptr->Num);
      else                           sprintf(str,"Qg%-d",ACptr->Num);
      i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"dg%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Vr%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Vi%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Ir%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Ii%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Vq%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Vd%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Iq%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      sprintf(str,"Id%-d",ACptr->Num); i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      if (!strpbrk(ACptr->cont,"I")) sprintf(str,"Ia%-d",ACptr->Num);
      else                           sprintf(str,"Qg%-d",ACptr->Num);
      i++; fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
  }
  for(k=0,DCptrR=dataPtr->DCbus;DCptrR!=NULL;DCptrR=DCptrR->Next){
    DCptrI=DCptrR->To;
    if(!strcmp(DCptrR->Type,"R")){
      k++;
      for(j=1;j<=2;j++){
        if(j==1) { DCptr=DCptrR; strcpy(type,"r"); }
        else { DCptr=DCptrI; strcpy(type,"i"); }
        if(strcmp(DCptr->Cont1,"VD")&&strcmp(DCptr->Cont2,"VD")) {
          sprintf(str,"Vd%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
        if(strcmp(DCptr->Cont1,"AT")&&strcmp(DCptr->Cont2,"AT")) {
          sprintf(str,"t%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
        if(strcmp(DCptr->Cont1,"AL")&&strcmp(DCptr->Cont2,"AL")) {
          sprintf(str,"al%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
        if(strcmp(DCptr->Cont1,"GA")&&strcmp(DCptr->Cont2,"GA")) {
          sprintf(str,"ga%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
        sprintf(str,"S%1s%-d",type,k); i++;
        fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        if(strcmp(DCptr->Cont1,"PA")&&strcmp(DCptr->Cont2,"PA")) {
          sprintf(str,"P%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
        if(strcmp(DCptr->Cont1,"QA")&&strcmp(DCptr->Cont2,"QA")) {
          sprintf(str,"Q%1s%-d",type,k); i++;
          fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
        }
      }
      if(strcmp(DCptrR->Cont1,"ID")&&strcmp(DCptrR->Cont2,"ID")&&
         strcmp(DCptrI->Cont1,"ID")&&strcmp(DCptrI->Cont2,"ID")) {
        sprintf(str,"Id%-d",k); i++;
        fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
      }
    }
  }
                             /* FACTS */
  for(k=0,SVCptr=dataPtr->SVCbus;SVCptr!=NULL;SVCptr=SVCptr->Next){
    k++;
    sprintf(str,"Qsvc%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Bv%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    if(!strcmp(SVCptr->Cont,"AL")){
      sprintf(str,"alpha%-d",k);i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    else {
      sprintf(str,"Vrefc%-d",k);i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
  }

  for(k=0,TCSCptr=dataPtr->TCSCbus;TCSCptr!=NULL;TCSCptr=TCSCptr->Next){
    k++;
    sprintf(str,"Ptcsc%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Qtcsck%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Qtcscm%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Be%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"alpha%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Itcsc%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"delta%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
  }

  for(k=0,STATCOMptr=dataPtr->STATCOMbus;STATCOMptr!=NULL;STATCOMptr=STATCOMptr->Next){
    k++;
    if(!strcmp(STATCOMptr->Cont,"PW") || !strcmp(STATCOMptr->Cont,"AL")){
      sprintf(str,"Istat%-d",k);i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    } else {
      sprintf(str,"Vrefc%-d",k);i++;
      fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    }
    sprintf(str,"theta%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Vdc%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"k%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"alpha%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Pstat%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
    sprintf(str,"Qstat%-d",k);i++;
    fCustomPrint(Out,"%4d %8s %-11.5g\n",i,str,vector[i]/Max);
  }
                             /* END FACTS */

  if (Option!='Y') {
    if (!Bl) fCustomPrint(Out,"%4d %8s %-11.5g\n",++i,"l",Dparam);
    else {
      sprintf(str,"V%-d",BlPtr->Num);
      fCustomPrint(Out,"%4d %8s %-11.5g\n",++i,str,Dparam);
    }
  }
  fCustomPrint(Out,"0 0 0.0\n");
  fclose(Out);
}


/* --------------------------- IndicesMatlab --------------------------------- */
#ifdef ANSIPROTO
void IndicesMatlab(INDEX count)
#else
void IndicesMatlab(count)
INDEX count;
#endif
/* Print plotting and other Matlab commands needed by the -0 option (VS indices) */
{
  char Namebase[80],Name[80];
  FILE *OutFile;
  INDEX i;
  

  strcpy(Namebase,NameParameter('0'));
  if(NullName(Namebase)) return;
  sprintf(Name,"%s.m",Namebase);
  OutFile=OpenOutput(Name);
  fCustomPrint(OutFile,"clear lambda evJ svJ evPF svPF evQV svQV detD_ll t_ll\n");
  fCustomPrint(OutFile,"clear crsvJ crevJ crsvPF crevPF crsvQV crevQV\n");
  fCustomPrint(OutFile,"warning off\n");
  for(i=1; i<=count; i++) fCustomPrint(OutFile,"%s%d\n",Namebase,i);
  fCustomPrint(OutFile,"figure; plot(lambda,evJ); \n");
  fCustomPrint(OutFile,"title('Full matrix |e-val.| index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"figure; plot(lambda,svJ); \n");
  fCustomPrint(OutFile,"title('Full matrix sing. val. index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"disp(' '); disp('Critical bus numbers and bus l rank for J indices are stored in crevJ and crsvJ')\n");
  fCustomPrint(OutFile,"figure; plot(lambda,evPF); \n");
  fCustomPrint(OutFile,"title('Power flow matrix |e-val.| index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"figure; plot(lambda,svPF); \n");
  fCustomPrint(OutFile,"title('Power flow matrix sing. val. index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"disp(' '); disp('Critical bus numbers and bus l rank for J_PV indices are stored in crevPF and crsvPF')\n");
  fCustomPrint(OutFile,"figure; plot(lambda,evQV,lambda,svQV,'-.'); \n");
  fCustomPrint(OutFile,"title('J_{QV} matrix |e-val.| and sing. val. indices');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"legend('e-v','s.v.');\n");
  fCustomPrint(OutFile,"disp(' '); disp('Critical bus numbers and bus l rank the J_QV indices are stored in crevQV and crsvQV')\n");
  fCustomPrint(OutFile,"figure; plot(lambda,detD_ll); \n");
  fCustomPrint(OutFile,"title('Reduced det. index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"figure; plot(lambda,-t_ll); \n");
  fCustomPrint(OutFile,"title('Test func. index');\n");
  fCustomPrint(OutFile,"xlabel('lambda [p.u.]');\n");
  fCustomPrint(OutFile,"disp(' '); disp('Bus l used for the red. det. and test func. indices is: %d %s')\n",TFbus,TFname);
  fCustomPrint(OutFile,"warning on\n\n");
  fclose(OutFile);
  /* --------------- Create 'rankbus.m' file needed for Matlab computations --------------- */
  OutFile=OpenOutput("rankbus.m");
  fCustomPrint(OutFile,"function num_rank=rankbus(vec,num)\n");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s Rank entry 'num' in abs(vec).\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"N=length(vec);\n");
  fCustomPrint(OutFile,"[val,I]=sort(abs(vec));\n");
  fCustomPrint(OutFile,"I=flipud(I);\n");
  fCustomPrint(OutFile,"for i=1:N,\n");
  fCustomPrint(OutFile,"  if (I(i)==num), num_rank=i; break; end\n");
  fCustomPrint(OutFile,"end\n");
  fclose(OutFile);
  /* ----------------  Create 'inviter.m' file needed for Matlab computations --------------- */
  OutFile=OpenOutput("inviter.m");
  fCustomPrint(OutFile,"function [e_val,v]=inviter(A,e_o,tol,iter,warn)\n");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s Inverse iteration method to compute a real e-value.\n","%%");
  fCustomPrint(OutFile,"%s Designed for sparse matrices, but works with full matrices.\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s      [v,e_val]=inviter(A,e_o,tol,iter,warn)\n","%%");
  fCustomPrint(OutFile,"%s \n","%%");
  fCustomPrint(OutFile,"%s Input:  A    -> NxN matrix\n","%%");
  fCustomPrint(OutFile,"%s         e_o  -> Optional eigenvalue guess (default 0)\n","%%");
  fCustomPrint(OutFile,"%s          tol  -> Optional convergence tolerance (default 1e-4)\n","%%");
  fCustomPrint(OutFile,"%s         iter -> Optional maximum number of iterations (default 30)\n","%%");
  fCustomPrint(OutFile,"%s         warn -> Use 0 to cancel display of no convergence warning message\n","%%");
  fCustomPrint(OutFile,"%s                 (default 1, i.e., display warning)\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s OutFileput: e_val -> Eigenvalue\n","%%");
  fCustomPrint(OutFile,"%s         v     -> Eigenvector\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"%s Copyright (c) Claudio Canizares, Shu Zhang, 1996, 2006.\n","%%");
  fCustomPrint(OutFile,"%s University of Waterloo, Waterloo, Canada\n","%%");
  fCustomPrint(OutFile,"%s\n","%%");
  fCustomPrint(OutFile,"if nargin<2, e_o=0; sign=1;\n");
  fCustomPrint(OutFile,"else\n");
  fCustomPrint(OutFile,"  if e_o<0, sign=-1;  else, sign=1;  end\n");
  fCustomPrint(OutFile,"end\n");
  fCustomPrint(OutFile,"if nargin<3, tol=1e-4; else, tol=abs(tol); end\n");
  fCustomPrint(OutFile,"if nargin<4, iter=30; else, iter=abs(iter); end\n");
  fCustomPrint(OutFile,"if nargin<5, warn=1; end\n");
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"N=size(A);\n");
  fCustomPrint(OutFile,"count=1; \n");
  fCustomPrint(OutFile,"conv=0;\n");
  fCustomPrint(OutFile,"J=sparse(A-e_o*eye(N(1)));\n");
  fCustomPrint(OutFile,"[L,U,P]=lu(J);\n");
  fCustomPrint(OutFile,"while (count<=2)\n");
  fCustomPrint(OutFile,"  v=ones(N(1),1);\n");
  fCustomPrint(OutFile,"  for I=1:iter,\n");
  fCustomPrint(OutFile,"   e_val=sign*1/norm(v,1);\n");
  fCustomPrint(OutFile,"   v=v*e_val;\n");
  fCustomPrint(OutFile,"   z=L\\(P*v);\n");
  fCustomPrint(OutFile,"   vp=U\\z;\n");
  fCustomPrint(OutFile,"   if(norm((1/e_val)*v-vp,1)<tol) break; end,\n");
  fCustomPrint(OutFile,"   v=vp;\n");
  fCustomPrint(OutFile,"  end\n");
  fCustomPrint(OutFile,"  if (I<iter)\n");
  fCustomPrint(OutFile,"    conv=1;\n");
  fCustomPrint(OutFile,"    break;\n");
  fCustomPrint(OutFile,"  else\n");
  fCustomPrint(OutFile,"    sign=-sign;\n");
  fCustomPrint(OutFile,"  end\n");
  fCustomPrint(OutFile,"  count=count+1;\n");
  fCustomPrint(OutFile,"end\n");
  fCustomPrint(OutFile,"e_val=e_o+sign*1/norm(vp,1);\n");
  fCustomPrint(OutFile,"v=e_val*vp;\n");
  fCustomPrint(OutFile,"\n");
  fCustomPrint(OutFile,"if (conv==0 & warn==1)\n");
  fCustomPrint(OutFile,"  disp(' ')\n");
  fCustomPrint(OutFile,"  disp('Warning: Inverse iteration method failed to converge')\n");
  fCustomPrint(OutFile,"  str=sprintf('         for tol=%s6.4e, iter=%sd.',tol,iter);\n","%%","%%");
  fCustomPrint(OutFile,"  disp(str)\n");
  fCustomPrint(OutFile,"  disp(' ')\n");
  fCustomPrint(OutFile,"end\n");
  fclose(OutFile);
}