/*下面是數據輸入格式說明：
為了減少不必要的運算，本潮流程序中將嚴格規定數據的錄入格式.
具體要求如下：
1.數據共分六大組：即基本數據組、支路數據組、補償電容數據組、并聯電抗數據組、發電機節點數據組、負荷節點數據組。

2.在基本數據組內按照以下格式錄入以下數據：
  節點個數，支路條數，對地補償電容支路數，并聯電抗數、發電機節點個數，負荷節點個數，系統平均電壓，計算精度

3.在支路數據組內按照以下格式錄入以下數據：
  對于線路支路：  支路的第一個節點號，支路的第二個節點號，支路的電阻，支路的電抗，支路的導納
  對于變壓器支路：變壓器的第一個節點號，變壓器的第二個節點號，變壓器的電阻，變壓器的電抗，變壓器的非標準變比
注意：（1）支路兩端的節點號應把小號排在前面，大號排在后面；
  （2）各支路按其小節點號的順序排列；
  （3）對于線路支路，導納應為Y,而不是計算中的Y/2;
  （4）對于變壓器支路，變壓器非標準變比所在側的節點號應為負號。
另外，串聯電容也應該在系統中被認為是一條支路，且串聯電容的參數應當換算成為相應的電抗值，注意：此值為負。
特別提醒：本程序不支持同節點之間的并聯線路，需要先行歸化為一條線路，然后求解。

4.在補償電容數據組內按照以下格式錄入以下數據：
  所在節點號、導納值

5.在并聯電抗數據組內按照以下格式錄入以下數據：
  所在節點號、電抗值

6.在發電機節點數據組內按照以下格式錄入以下數據：
  發電機所在節點號，發電機所發有功，發電機所發無功，發電機節點的電壓
注意：（1）對于PV節點，無功應為發電機的無功上限；
  （2）對于PV節點，電壓為節點所需維持的電壓，應為負值；
  （3）對于PQ節點，電壓為計算過程中的節點的電壓初值。

7.在負荷節點數據組內按照以下格式錄入以下數據：
  負荷所在節點號，負荷有功，負荷無功，負荷節點的電壓
注意：（1）負荷節點一般為PQ節點，其電壓為計算過程中的節點的電壓初值；
  （2）負荷的有功和無功均應寫成注入節點功率的形式，即全部為負值；
  （3）特別的，系統的平衡節點應該被排在最后，即作為第N個節點，同時要求該節點為負荷節點；如果該點無負荷，則該點的負荷功率填零。

至于對所求解系統的描述性文字可以放在該文件的末尾，任意長度。
另外，本程序可以自動形成輸出文件，其文件名為輸入文件名后加"out"，且輸出文件名的后綴確定為".dat"。例如，對于程序已給的5節點系統，輸入文件名"ps-5"，后綴為".dat"，輸出文件名為"ps-5out"，后綴為".dat"。為了保證名稱的前后統一，因此建議：將輸入文件的名字寫成如下格式：****.dat;并且輸入文件的文件名的長度應不大于5（不計后綴）。











*/


/**********************************************************************/
/*****************PQ分解法潮流程序-BB方案,允許迭代次數2000次******************************/
/*********
        作者：李常剛，山東大學電氣工程學院02級3班，spring-_lee@163.com
        源于畢業設計題目，2006年4月3日
**********/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

  struct Branch_Type
  {
   int i,j;
   double R,X,YK;
  }; 
  struct Compensation_Type
  {
   int i;
   double Y;
  };
  struct Reactance_Type
  {
   int i;
   double X;
  };
  struct Generator_Type
  {
   int i;
   double P,Q;
   double V;
  };
  struct Load_Type
  {
   int i;
   double P,Q;
   double V;
  };
  struct PVNode_Type
  {
   int i;
   double V;
  };
  struct Yii_Type
  {
   double G,B;
  };
  struct Yij_Type
  {
   double G,B;
   int j;
  };
  struct NodalPow
  {
   double P,Q;
  };
  struct NodalVol
  {
    double V,theta;
  };
  struct GeneratorPower
  {
     double P,Q;
  };
  struct U_Type
  {
     double value;
     int j;
  };

  void Datain(int Nb,int Nl,int Nc,int Nr,int Ng,FILE *fp,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance,struct Load_Type *Load,struct Generator_Type *Generator,struct PVNode_Type *PVNode,int *Npv);/*從文件中讀入數據*/
  void AdmittanceMatrix(int N,int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,int *NYseq1,int *NYsum1);/*導納矩陣形成函數*/
  void AdmittanceMatrixAdd(int Nb,int Nc,int Nr,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance);/*導納矩陣追加接地支路函數*/
  void Factorial(int flag,int N,int Npv,struct PVNode_Type *PVNode,int *NUsum,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq,double *D,struct U_Type *U);/*因子表形成函數*/
  void NodePower(int flag,int N,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq);/*節點功率計算函數*/
  void Iteration(int flag,struct Generator_Type *Generator,struct Load_Type *Load,struct PVNode_Type *PVNode,struct NodalVol *NodeVol,struct NodalPow *NodePow,struct GeneratorPower *GenPower,int N,double *DI,double *MaxError,int *ErrNode);/*迭代函數*/
  void FormulaSolution(int flag,struct U_Type *U,double *D,int *NUsum,double *DI,int N,struct NodalVol *NodeVol,double V0);
  void NodeDataOutput(FILE *fp,struct NodalVol *NodeVol,struct Generator_Type *Generator,int N,struct GeneratorPower *GenPower,struct NodalPow *NodePow,struct Load_Type *Load,int Nl);
  void BranchDataOutput(FILE *fp,int Nb,int Nc,int Nr,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance,struct NodalVol *NodeVol);









main()
{
  int N;
  int Nb;
  int Nc;
  int Nr;
  int Ng;
  int Nl;
  double V0;
  double epsilon;
  struct Branch_Type *Branch;
  struct Compensation_Type *Compensation;
  struct Reactance_Type *Reactance;
  struct Generator_Type *Generator;
  struct Load_Type *Load;
  struct PVNode_Type *PVNode;
  int Npv=0;/*標記PV節點總個數*/
  struct Yii_Type *Yii,*Yiil;
  struct Yij_Type *Yij,*Yijl;
  int *NYseq,*NYsum;
  struct NodalPow *NodePow;
  struct NodalVol *NodeVol;
  struct GeneratorPower *GenPower;
  double *DI1,*DI2;
  double MaxError=0.0,MaxErrorTemp;
  int ErrNode,ErrNodeTemp;
  int Kp=1,Kq=1;


  int n,k,count;
  time_t now;
 



  int i;/*用以標記當前處理行號*/





 

  
  struct U_Type *U1,*U2;
  double *D1,*D2;
  int *NUsum1,*NUsum2;
  

  char FILENAME[20];
  FILE *fp;



 /* if((fp=fopen("README.txt","r"))==NULL)
  {
      printf("部分文件缺失，請確保文件的完整性!!\n");
      printf("按 ENTER 鍵返回！\n");
      scanf("%c",&ch);
      exit(1);
  }
*/

  /***********打開系統數據文件**********************/
  printf("Please Enter The Filename of The System:");
  gets(FILENAME);
  if((fp=fopen(FILENAME,"r"))==NULL)
  {
    printf("Cannot Find The File:%s\n",FILENAME);
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }



 /******************讀系統數據信息****************/
  fscanf(fp,"%d,%d,%d,%d,%d,%d,%lf,%lf",&N,&Nb,&Nc,&Nr,&Ng,&Nl,&V0,&epsilon);


 /******************給支路、發電機、負荷和PV節點分配內存******************/
  Branch=malloc((Nb+1)*sizeof(struct Branch_Type));
  if(Branch==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
  Compensation=malloc((Nc+1)*sizeof(struct Compensation_Type));
  if(Compensation==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
  Reactance=malloc((Nr+1)*sizeof(struct Reactance_Type));
  if(Reactance==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
  Generator=malloc((Ng+1)*sizeof(struct Generator_Type));
  if(Generator==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
  Load=malloc((Nl+1)*sizeof(struct Load_Type));
  if(Load==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
  PVNode=malloc(N*sizeof(struct PVNode_Type));
  if(PVNode==NULL)
  {
    printf("Cannot Get Memory for That Many Values!");
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
/**********內存分配完成***************/





  Datain(Nb,Nl,Nc,Nr,Ng,fp,Branch,Compensation,Reactance,Load,Generator,PVNode,&Npv);/*從文件中讀入數據*/




  for(n=0;1;n++)/*打開輸出文件*/
  {
    if(FILENAME[n]=='.')
    {
        FILENAME[n]='\0';
        strcat(FILENAME,"out.dat");
        break;
    }
  }
  if((fp=fopen(FILENAME,"w"))==NULL)/*打開數據輸出文件*/
  {
    printf("Cannot Find The File:%s\n",FILENAME);
    printf("Press ENTER to Escape!");
    getchar();
    exit(1);
  }
/************為導納矩陣分配內存,并形成導納矩陣*******************/

  Yii=malloc((N+1)*sizeof(struct Yii_Type));
  Yiil=malloc((N+1)*sizeof(struct Yii_Type));
  Yij=malloc((N+1)*sizeof(struct Yij_Type));//???
  Yijl=malloc((N+1)*sizeof(struct Yij_Type));
  NYseq=malloc((N+1)*sizeof(int));
  NYsum=malloc((N+1)*sizeof(int));


  AdmittanceMatrix(N,Nb,Yii,Yiil,Yij,Yijl,Branch,NYseq,NYsum);/*導納矩陣形成函數*/



/***********分配及導納矩陣的形成完成*******************/


/******************下面確認導納矩陣的正確性（僅用于調試）*********************

  fprintf(fp,"The Admittance Matrix is:\n");
  printf("The Admittance Matrix is:\n");

  for(n=1;n<=N;n++)
  {
    fprintf(fp,"Yii[%d]=%lf+j%lf\n",n,Yii[n].G,Yii[n].B);
    fprintf(fp,"Yiil[%d]=%lf+j%lf\n",n,Yiil[n].G,Yiil[n].B);


    printf("Yii[%d]=%lf+j%lf\n",n,Yii[n].G,Yii[n].B);
    printf("Yiil[%d]=%lf+j%lf\n",n,Yiil[n].G,Yiil[n].B);
  }
  for(n=1;n<=N-1;n++)
  {
    for(i=NYseq[n];i<=NYseq[n+1]-1;i++)
    {
        fprintf(fp,"Yij[%d][%d]=%lf+j%lf\n",n,Yij[i].j,Yij[i].G,Yij[i].B);
        fprintf(fp,"Yijl[%d][%d]=%lf+j%lf\n",n,Yijl[i].j,Yijl[i].G,Yijl[i].B);

        printf("Yij[%d][%d]=%lf+j%lf\n",n,Yij[i].j,Yij[i].G,Yij[i].B);
        printf("Yijl[%d][%d]=%lf+j%lf\n",n,Yijl[i].j,Yijl[i].G,Yijl[i].B);
    }
  }

*****************確認導納矩陣的正確性完成（僅用于調試）*********************/



/*****************下面形成因子表,采用BX方案(在B'中不忽略電阻，在B"中忽略電阻)*****************/
  U1=malloc((N-1)*(N-2)/2*sizeof(struct U_Type));
  U2=malloc((N-1)*(N-2)/2*sizeof(struct U_Type));
  D1=malloc(N*sizeof(double));
  D2=malloc(N*sizeof(double));
  NUsum1=malloc(N*sizeof(int));
  NUsum2=malloc(N*sizeof(int));

  Factorial(1,N,Npv,PVNode,NUsum1,Yii,Yij,NYseq,D1,U1);/*因子表形成函數*/
/**************************

  count=0;
  for(n=1;n<=N-1;n++)
  {
      fprintf(fp,"D[%d]=%lf;\n",n,D1[n]);

     }
  for(i=1;i<=N-1;i++)
  {
      
      for(n=1;n<=NUsum1[i];n++)
      {
          fprintf(fp,"U[%d][%d]=%lf\n",i,U1[count+n].j,U1[count+n].value);

      }
      count+=NUsum1[i];
  }


    count=0;

**************************************/
  AdmittanceMatrixAdd(Nb,Nc,Nr,Yii,Yiil,Yij,Yijl,Branch,Compensation,Reactance);/*導納矩陣追加接地支路函數*/
  Factorial(2,N,Npv,PVNode,NUsum2,Yii,Yij,NYseq,D2,U2);/*因子表形成函數*/


/********************************

  for(n=1;n<=N-1;n++)
  {
      fprintf(fp,"D[%d]=%lf;\n",n,D2[n]);
 
  }
  for(i=1;i<=N-1;i++)
  {
      
      for(n=1;n<=NUsum2[i];n++)
      {
         fprintf(fp,"U[%d][%d]=%lf\n",i,U2[count+n].j,U2[count+n].value);
      }
      count+=NUsum2[i];
  }



    
 
*************************************/



/****************下面利用所求得的因子表進行迭代求解******************/

   

  DI1=malloc(N*sizeof(double));
  DI2=malloc(N*sizeof(double));
  NodePow=malloc((N+1)*sizeof(struct NodalPow));
  NodeVol=malloc((N+1)*sizeof(struct NodalVol));


  /****先送電壓初值******/
  for(i=1;i<=N;i++)
  {
    NodeVol[i].V=V0;
    NodeVol[i].theta=0.0;
  }
  /*for(i=1;i<=Ng;i++)
  {
      n=Generator[i].i;
      NodeVol[n].V=Generator[i].V;
      NodeVol[n].theta=0.0;
  }
   for(i=1;i<=Nl;i++)
  {
      n=Load[i].i;
      NodeVol[n].V=Load[i].V;
      NodeVol[n].theta=0.0;
  }*/

  for(n=1;n<=Npv;n++)
  {
    i=PVNode[n].i;
    NodeVol[i].V=PVNode[n].V;
  }