1.cpp

生成导纳矩阵

//此程序可用于求复数矩阵的逆.

#include <iostream>
#include <iomanip>
#include < fstream >
#include < cstring >
using namespace std;

#include "1.h"

//////定义全局变量
int N;
int Nb;
int Ng;
int Nl;

NodeDATA *Node;
BranchDATA *Branch ;
Yii_type * Yii;
Yii_type * Y1ii;
Yij_type * Yij;
Yij_type * Y1ij;
int *NYsum;
int * NYseq;
complex **c;


void print(complex z2)
{
 //if (z2.imag>=0)
  cout<<setw(12)<<z2.real<<'+'<<'i'<<z2.imag;
 /*else
  cout<<setw(12)<<z2.real<<'-'<<'i'<<z2.imag;*/
}


//主程序入口函数main
void main()
{  
    std :: cout << "Beginning the flowing calculation:" << std :: endl
			<< "please input the file which have the information (*.dat):" << std :: endl;
    char ch[10];
	std :: cin >> ch;
	std :: cout << "get the file name :"<< ch << std ::endl  ;
	readfile ( ch );//读原始数据文件

    cout<<"电路的节点数N："<<N<<endl;
    BuildY1();
	OutputY();
    
	cout<<"接下来，对上面生成的导纳矩阵求逆!!"<<endl;
   
    int i,j,k,n,p,q,h;
	n=N;
    complex **a,**b,temp,temp2;
   
    a=new complex*[n+1];
    b=new complex*[n+1]; //
   
    for (i=1;i<=n;i++)
    {
      a[i]=new complex[n+1];
      b[i]=new complex[n+1];
    } 
   
   //单位矩阵
   for (i=1;i<=n;i++)
   {
     for (j=1;j<=n;j++)
	 {   
         if (i==j)
		 {
           b[i][j].real=1;
           b[i][j].imag=0;
		 }
         else
		 {
           b[i][j].real=0;
           b[i][j].imag=0;
		 }
	 }
   }
   
  cout<<"上面得到的完整的导纳矩阵为"<<endl;
  for (i=1;i<=n;i++)
   {
        for (j=1;j<=n;j++)
		{ 
		  a[i][j].real=c[i][j].real;
		  a[i][j].imag=c[i][j].imag;
		  print(a[i][j]);
		}
            
     cout<<endl;
     cout<<endl;
   }
   
   
   
   for (i=1;i<=n;i++)
   {
       for (k=i;k<=n;k++)//此语句用于判断该矩阵是否有逆
	   {
           if (a[k][i].real!=0 || a[k][i].imag!=0)
		   {
              h=k;//h记住每列中从对角线起第一个非零的元素所在行数
              k=n+1;
		   }
	   }
      
	   if (k==n+1)
	  {
        cout<<"您输入的矩阵没有逆!!!"<<endl;
        i=n+1;
	  }
      
	  else//重点
	  {
          if (h!=i)
		  {
              for (j=1;j<=n;j++)
			  {
                a[i][j]=a[i][j]+a[h][j];
                b[i][j]=b[i][j]+b[h][j];
			  }
		  }
       temp=a[i][i];
       for (j=1;j<=n;j++)
	   {
          a[i][j]=a[i][j]/temp;//把对角上的数单位化
          b[i][j]=b[i][j]/temp;
	   }
       for (p=1;p<=n;p++)
	   {
           if (i!=p)//此行不用动
		   {
             temp2=a[p][i];
             for (q=1;q<=n;q++)
			 {
                 a[p][q]=a[p][q]-temp2*a[i][q];
                 b[p][q]=b[p][q]-temp2*b[i][q];
			 }
		   }//if
	   }
	  }//else
   }//for
  
   
   if (k!=n+1)
   {
         cout<<"所求导纳矩阵的逆矩阵为:"<<endl;
        for (i=1;i<=n;i++)
		{
          for (j=1;j<=n;j++)
		     print(b[i][j]);
		  
       cout<<endl;
       cout<<endl;
		}
   }//if
   
    
   delete []a;
    delete []b;
   	delete [] Yii;
	delete [] Y1ii;
	delete [] Yij;
	delete [] Y1ij;
	delete [] NYsum ;
	delete [] NYseq ;
	delete [] Node;
	delete [] Branch;
}//main

///生成导纳矩阵的子程序////
/*******************************************/
void readfile ( char filename[10] )
{	
	////读原始数据文件
	
	ifstream file ( filename );	
	///如果文件读取失败，输出提示
	if (file.fail())
	{
		std ::cout << "Cannot open file " << filename << std :: endl;
		return;
	}
	//读入节点号，支路个数，发电机节点个数，负荷节点个数
	file >> N >> Nb >> Ng >>Nl;
	/*
	std ::cout <<"节点号，支路个数，发电机节点个数，负荷节点个数分别为："<< std ::endl 
			<< N << "	" << Nb  << "	" 
			<< Ng  << "	" << Nl << "\n" ;
	*/
	//生成动态数组，用来存放节点的信息
	Node = new NodeDATA[N];
	Branch = new BranchDATA [Nb];

	
	//读取支路信息
	for (int  m = 0 ;m < Nb ;m++ )
	{
		file >> Branch [m].i >> Branch [m].j >> Branch [m] .R 
			>> Branch [m] .X >>Branch [m] .Yk >>Branch [m] .type;
	}

	return ;
}
/*******************************************/
void BuildY1 ()
{     
	////建立导纳矩阵

	Yii = new Yii_type[N];
	Y1ii = new Yii_type[N];
	Yij = new Yij_type[Nb];
	Y1ij = new Yij_type[Nb];
	NYsum = new int [N];
	NYseq = new int [N];
	//////////////////////////
	//////学生编写程序部分
	//////节点导纳矩阵的形成
	///////
	//////////////////////////
	for ( int i =0; i < N; i++)////i=0;
	{
		Yii[i].G = 0.0;
		Yii[i].B = 0.0;
		Yij[i].G = 0.0;
		Yij[i].B = 0.0;
		NYsum[i] = 0;
	}									//清零

	for(int k=0;k<Nb;k++)				
	{
		int i=abs(Branch[k].i);			//节点号
        int j=abs(Branch[k].j);
        double R=Branch[k].R;			//电阻
		double X=Branch[k].X;			//电抗
		double YK=Branch[k].Yk;			//非标准变比
		double Zmag2=R*R+X*X;
		double Gij=R/Zmag2;				//导纳实部
		double Bij=-X/Zmag2;			//导纳虚部
        			
		if ((Branch[k].i<0) ||(Branch[k].j<0))
		{	Yij[k].G=-Gij/YK;
		    Yij[k].B=-Bij/YK;
        }								//变压器
		else
		{
			Yij[k].G=-Gij;
		    Yij[k].B=-Bij; 
		}								//线路
		Yij[k].j=j;
		
										
       if ((Branch[k].i<0)||(Branch[k].j<0))
		{
			
		   if(Branch[k].i<0)            //变压器在支路的i端
		   {
			    i=abs(i);               
				Yii[i-1].G=Yii[i-1].G+Gij/YK/YK;
				Yii[i-1].B=Yii[i-1].B+Bij/YK/YK;
			    Yii[j-1].G=R/Zmag2;
		    	Yii[j-1].B=-X/Zmag2;
		   }
		   else							//变压器在j端
		   {
			    j=abs(j);              
				Yii[j-1].G=Yii[j-1].G+Gij/YK/YK;
				Yii[j-1].B=Yii[j-1].B+Bij/YK/YK;
				Yii[i-1].G=Yii[i-1].G+Gij;
				Yii[i-1].B=Yii[i-1].B+Bij;
			
		   }
           
		}
		else
		{	
			
			Yii[i-1].G=Yii[i-1].G+Gij;
			Yii[i-1].B=Yii[i-1].B+Bij+YK/2;
            Yii[j-1].G=Yii[j-1].G+Gij;
		    Yii[j-1].B=Yii[j-1].B+Bij+YK/2;
			   
		}
		NYsum[i-1]=NYsum[i-1]+1;				//不接地支路个数，累计非对角元素非零的个数
	}
	NYseq[0]=1;								//每一行非零元素在导纳矩阵的序号
	for( i=0;i<N-1;i++)
	{
		NYseq[i+1]=NYseq[i]+NYsum[i];
	}
	return ;
}

/*******************************************/
void OutputY ()
{  ////输出节点导纳矩阵
	c=new complex*[N+1];
	for (int m=1;m<=N;m++)
     c[m]=new complex[N+1];
     
   
	std :: cout <<"节点导纳矩阵的上三角矩阵如下："<<"\n";
	for (int i = 0; i < N ; i++)
	{	
		for (int j = i ;j < N;j++)
			{
				int k = 0;
				if (i == j)
				{
					std ::cout << Yii[i].G<<"+i"<< Yii[i].B << "	";
					c[i+1][j+1].real=Yii[i].G;
					c[i+1][j+1].imag=Yii[i].B;
				}
				else
				{	
					for (int t = NYseq[i];t < NYseq[i+1]; t++)
					{
						if(j == (Yij[t-1].j-1))
						{
							k = -1;
							std ::cout << Yij [t-1].G<<"+i"<< Yij[t-1].B<< "	";
                               	c[i+1][j+1].real=Yij[t-1].G;
                                c[j+1][i+1].real=Yij[t-1].G;
						        c[i+1][j+1].imag=Yij[t-1].B;
								c[j+1][i+1].imag=Yij[t-1].B;
						}
					}
					if ( k == 0)
					{
						std ::cout << "0"<<"+i"<<"0"<< "	";
					}
				}
			}
		std ::cout <<"\n";
	}
	return ;
}