P-Q潮流计算

＃i nclude <cmath>    
＃i nclude <iostream> 
#define N 10
 //#define n 5      /*备用*/
using namespace std;

int main()
{
 int n;        //统计节点个数
 int i,j,k;      //循环变量
 int r;          //PV节点个数
    //int r=0;
    int m;          //PQ节点个数
 double R[N][N];/*R[5][5]={{0,0.02,0.08,0,0},{0.02,0,0.06,0.06,0.04},{0.08,0.06,0,0.01,0},
 {0,0.06,0.01,0,0.08},{0,0.04,0,0.08,0}};*/                //节点间的电阻
 double X[N][N];/*X[5][5]={{0,0.06,0.24,0,0},{0.06,0,0.18,0.18,0.12},{0.24,0.18,0,0.03,0},
 {0,0.18,0.03,0,0.24},{0,0.12,0,0.24,0}};*/    //节点间的电抗
 double g1[N][N]={0};    //节点间的电导
 double b1[N][N]={0};    //节点间的电纳
 double G[N][N]={0};/*={{6.25,-5.0,-1.25,0,0},{-5.0,10.834,-1.667,-1.667,-2.5},
                     {-1.25,-1.667,12.917,-10.0,0},{0,-1.667,-10.0,12.917,-1.25},
      {0,-2.5,0,-1.25,3.75}};
                      导纳矩阵电导*/ 
 double B[N][N]={0};/*={{-18.75,15.0,3.75,0,0,},{15.0,-32.5,5.0,5.0,7.5},
                        {3.75,5.0,-38.75,30.0,0},{0,5.0,30.0,-38.75,3.75},
         {0,7.50,0,3.75,-11.25}};
                        导纳矩阵电纳*/

 double P[N]={0};//={0.2,-0.45,-0.4,-0.6}; //节点有功初值，平衡节点除外
 double Q[N]={0};//={0.2,-0.15,-0.05,-0.1}; //节点无功初值，平衡节点除外
 //double P[5]={0,0.2,-0.45,-0.4,-0.6}; //节点有功初值，平衡节点除外
 //double Q[5]={0,0.2,-0.15,-0.05,-0.1}; //节点无功初值，平衡节点除外
 
 double B1[N-1][N-1];     //系数矩阵1 
 double B2[N-1][N-1];     //系数矩阵2
 double A1[N-1][N-1];     //系数矩阵1的逆矩阵
 double A2[N-1][N-1];     //系数矩阵2的逆矩阵 
 
 double a[N][2*N];        //求逆矩阵时用的矩阵
 double w=0,s=0,l,m1;     //求逆矩阵时用的变量
 
 double U[N];             //节点电压幅值
 double rad[N];           //节点电压幅角
 //double U[N]={1.06,1,1,1,1};
 //double rad[N]={0,0,0,0,0};

 double cp[N];            //有功不平衡量
 double cq[N];            //无功不平衡量
 

    cout << "请输入网络节点数：";
 cin  >>  n ;
 cout << "请输入PV节点个数(r<n):" << "r=" ;
 cin  >> r;
 m=n-r-1;
 cout << "PQ节点个数(n-r-1)为："   << "m=" << m << endl;
 cout << "请输入节点间电阻：\n";  /*************************************/
 for( i=0 ; i<n ; i++ )      /* 0      0.02    0.08    0     0    */
 {                                  /* 0.02   0       0.06    0.06  0.04 */
  for( j=0 ; j<n ; j++ )         /* 0.08   0.06    0       0.01  0    */
  {                              /* 0      0.06    0.01    0     0.08 */
   cin >> R[i][j];            /* 0      0.04    0       0.08  0    */
  }                         /*************************************/
 }
 
 cout << "\n请输入节点间电抗：\n";  /*************************************/
 for( i=0 ; i<n ; i++ )            /* 0       0.06   0.24    0     0    */  
 {                                  /* 0.06    0      0.18    0.18  0.12 */
  for( j=0 ; j<n ; j++ )         /* 0.24    0.18   0       0.03  0    */
  {                              /* 0       0.18   0.03    0     0.24 */
   cin >> X[i][j];            /* 0       0.12   0       0.24  0    */
  }                         /*************************************/
 }
/*********************************************************************************************
接着计算节点间的导纳值
*********************************************************************************************/
 for( i=0 ; i<n ; i++ )
 {
  for( j=0 ; j<n ; j++ )
  {
   if( (R[i][j]!=0)||(X[i][j]!=0) )
   {
    g1[i][j]=R[i][j]/(R[i][j]*R[i][j]+X[i][j]*X[i][j]);
    b1[i][j]=-X[i][j]/(R[i][j]*R[i][j]+X[i][j]*X[i][j]);
   }
   else
   {
    g1[i][j]=0;
    b1[i][j]=0;    
   }
  }
 }
/*********************************************************************************************
下面形成导纳矩阵
*********************************************************************************************/
 for( i=0 ; i<n ; i++ )
 {
  for( j=0 ; j<n ; j++ )
  {
   if(j!=i)  G[i][j]=-g1[i][j];
   else for( k=0 ; k<n ; k++ )
   {
    G[i][j]+=g1[i][k];
   }
  }
 }

 for( i=0 ; i<n ; i++ )
 {
  for( j=0 ; j<n ; j++ )
  {
   if(j!=i)  B[i][j]=-b1[i][j];
   else for( k=0 ; k<n ; k++ )
   {
    B[i][j]+=b1[i][k];
   }
  }
 }
 /*for( i=0 ; i<n ; i++ )
 {
  cout << endl;
  for( j=0 ; j<n ; j++ )
  {
   cout << B[i][j] << " ";
  }
 }*/
/*********************************************************************************************
下面形成系数矩阵B1，B2
*********************************************************************************************/
 for( i=0 ; i<n-1 ; i++ )
 {
  for( j=0 ; j<n-1 ; j++ )
  {
   B1[i][j]=B[i+1][j+1];
  }
 }
 for( i=0 ; i<n-r-1 ; i++ )
 {
  for( j=0 ; j<n-r-1 ; j++ )
  {
   B2[i][j]=B[i+1][j+1];
  }
 }

 /*for( i=0 ; i<n-1 ; i++ )                         /*输出系数矩阵B1
 {
  cout << endl;
  for( j=0 ; j<n-1 ; j++ )
  {
   cout << B1[i][j] << "   " ;
  }
 }                                
 for( i=0 ; i<n-r-1 ; i++ )                        /*输出系数矩阵B2
 {
  cout << endl;
  for( j=0 ; j<n-r-1 ; j++ )
  {
   cout << B2[i][j] << "   ";
  }
 }*/
/*********************************************************************************************
下面形成系数矩阵B1逆矩阵A1
*********************************************************************************************/
 int n1=n-1;
 for (i=0;i<n1;i++)
 {
  for (j=0;j<n1;j++)
        {
   a[i][j]=B1[i][j];
  }
 }
 for (i=0;i<n1;i++)
 {
  for (j=n1;j<2*n1;j++)
  {
   if((i+n1)==j)
             a[i][j]=1;
   else
    a[i][j]=0;
  }
 }
 for (i=0;i<n1;i++)
   {
   w=a[i][i];  
   for(j=i;j<2*n1;j++)
 {
  a[i][j]=a[i][j]/w;
    }
    for(k=i+1;k<n1;k++)
    {
           s=a[k][i];
     for(j=i;j<2*n1;j++)
     {
              a[k][j]=a[k][j]-s*a[i][j];
     }
    }
   }
 for(i=n1-1;i>=0;i--)
  {
   l=a[i][i];
   for(j=i;j<2*n1;j++)
   {
    a[i][j]=a[i][j]/l;
   }
    for(k=i-1;k>=0;k--)
    {
     m1=a[k][i];
           for(j=i;j<2*n1;j++)
     {
      a[k][j]=a[k][j]-m1*a[i][j];
     }
    }
  }
 for(i=0;i<n1;i++)
 {
  for(j=n1;j<2*n1;j++)
  {
    A1[i][j-n1]=a[i][j] ;
  }
 }
 /*for(i=0;i<n1;i++)
 {
  cout << endl;
  for(j=0; j<n1 ; j++ )
  {
   co