支路追加法求阻抗阵.cpp

电网计算时

#include "math.h"
#include "stdio.h"
#define L 20                                   /*数组大小*/
//复数求倒子函数
void compute_reciprocal( float *z0_real,float *z0_imag,float *zr_real,float *zr_imag)
{
	float t1 ;
	t1=*z0_real**z0_real+*z0_imag**z0_imag;
    *zr_real=*z0_real/t1;
    *zr_imag=-1**z0_imag/t1;
}
//复数矩阵求逆的子函数
void reciprocal_complex_matrix(int n,float c0_real[L][L],float c0_imag[L][L],float cr_real[L][L],float cr_imag[L][L])
{
	int i,j,k,m,p,q;
	float t;
	float M[L][L]={0};//增广矩阵
    for(i=0;i<n;i++)  
	{
		for(j=0;j<n;j++)
			M[i][j]=c0_real[i][j];
		for(j=n;j<2*n;j++)
			M[i][j]=-1*c0_imag[i][j-n];
		M[i][i+2*n]=1;
	}
	for(i=n;i<2*n;i++)
	{
		for(j=0;j<n;j++)
			M[i][j]=c0_imag[i-n][j];
		for(j=n;j<2*n;j++)
			M[i][j]=c0_real[i-n][j-n];
	}
	//对M做2n次列主元Gauss消去
	for(j=0;j<2*n;j++)
	{
        for(i=j,m=j;i<2*n;i++)
			if(float(fabs(float(M[i][j])))>float(fabs(float(M[m][j]))))
				m=i;                         //找到主元所在行
		for(k=0;k<3*n;k++)
		{
			t=M[j][k];                       //主元所在行与对角元所在行交换
            M[j][k]=M[m][k];
		    M[m][k]=t;
		}
	    for(i=j+1;i<2*n;i++)
			for(k=3*n-1;k>=j;k--)
				M[i][k]=M[i][k]-M[j][k]*M[i][j]/M[j][j];//化为上三角阵
	}
	for(j=2*n-1;j>0;j--)
		for(i=0;i<j;i++)
			for(k=3*n-1;k>=0;k--)
				M[i][k]=M[i][k]-M[j][k]*M[i][j]/M[j][j];//化为单位阵
	for(i=0;i<2*n;i++)
		for(k=3*n-1;k>=0;k--)
			M[i][k]=M[i][k]/M[i][i];         //对角元化为1
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			cr_real[i][j]=M[i][j+2*n];
			cr_imag[i][j]=M[i+n][j+2*n];
		}
}

void main()               /*主函数*/
{
	int i,j,k,JF,JT,IFM[L],ITO[L];                  /*临时变量*/
    float R=0.0,X=0.0,G,B,r,x,Zll_real,Zll_imag,Zl_real[L]={0},Zl_imag[L]={0},K_real[L][L]={0},K_imag[L][L]={0};//,K0[L][L]={0},B1[L]={0},B2[L]={0},B0[L]={0};
    int Nnode;   /*节点数*/
    int Nbr;    /*支路数*/
    float Y0_real[L][L]={0},Y0_imag[L][L]={0};         /*原网的导纳矩阵*/
    float Z0_real[L][L]={0},Z0_imag[L][L]={0};         /*原网的阻抗矩阵*/
    int node0;   //追加树支的初始节点号
    int node_new;//追加树支的新增节点号
    float Radd_n;//追加树支的电阻
    float Xadd_n;//追加树支的电抗
    int node_i,node_j;//追加连支所在的两个节点号
    float Radd_br;    //追加连支的电阻
    float Xadd_br;    //追加连支的电抗
    float Znew_real[L][L]={0},Znew_imag[L][L]={0};         /*新网的阻抗矩阵*/
    FILE *fp;
/*读入网络原始数据,形成网络的节点导纳矩阵*/
    if((fp=fopen("原始数据输入.txt","r"))==NULL)
    {
		printf("Can not open the file.\n");
        return;
    }
    fscanf(fp,"节点数Nnode=%d\n支路数Nbr=%d\n",&Nnode,&Nbr);
    fscanf(fp,"电网络的参数(依次为支路编号，支路一端的节点号，支路另一端的节点号，支路电阻，支路电抗，其中0号节点为接地点）\n");
	for(j=0;j<Nbr;j++)
	{
		fscanf(fp,"%d\t%d\t%d\t%f\t%f\n",&i,&JF,&JT,&r,&x);
		IFM[i]=JF;
        ITO[i]=JT;
        R=r;
		X=x;
        compute_reciprocal(&R,&X,&G,&B);
        Y0_real[JF][JF]=Y0_real[JF][JF]+G;
        Y0_imag[JF][JF]=Y0_imag[JF][JF]+B;
        Y0_real[JT][JT]=Y0_real[JT][JT]+G;
        Y0_imag[JT][JT]=Y0_imag[JT][JT]+B;
        Y0_real[JF][JT]=Y0_real[JF][JT]-G;
        Y0_imag[JF][JT]=Y0_imag[JF][JT]-B;
        Y0_real[JT][JF]=Y0_real[JT][JF]-G;
        Y0_imag[JT][JF]=Y0_imag[JT][JF]-B;
	}
    fscanf(fp,"追加树支的初始节点号node0=%d,新增节点号node_new=%d\n",&node0,&node_new);
    fscanf(fp,"追加树支的阻抗：Radd_n=%f,Xadd_n=%f\n",&Radd_n,&Xadd_n);
    fscanf(fp,"追加连支所在的两个节点号node_i=%d,node_j=%d\n",&node_i,&node_j);
    fscanf(fp,"追加连支的阻抗：Radd_br=%f,Xadd_br=%f\n",&Radd_br,&Xadd_br);
    fclose(fp);
	for(i=0;i<Nnode;i++)
		for(j=0;j<Nnode;j++)
		{
			Y0_real[i][j]=Y0_real[i+1][j+1];
			Y0_imag[i][j]=Y0_imag[i+1][j+1];
		}
//导纳阵求逆得到阻抗阵
	reciprocal_complex_matrix(Nnode,Y0_real,Y0_imag,Z0_real,Z0_imag);
//追加树支后的阻抗阵
    for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
		{
			Znew_real[i][j]=Z0_real[i][j];
			Znew_imag[i][j]=Z0_imag[i][j];
		} 
        Znew_real[i][node_new-1]=Z0_real[i][node0-1];
		Znew_imag[i][node_new-1]=Z0_imag[i][node0-1];
        Znew_real[node_new-1][i]=Z0_real[node0-1][i];
		Znew_imag[node_new-1][i]=Z0_imag[node0-1][i];
	}
    Znew_real[node_new-1][node_new-1]=Z0_real[node0-1][node0-1]+Radd_n;
    Znew_imag[node_new-1][node_new-1]=Z0_imag[node0-1][node0-1]+Xadd_n;
//将结果保存到"运行结果.txt"中
   if((fp=fopen("运行结果.txt","w+"))==NULL)
   {
	   printf("Can not open the file.\n");
       return;
   }
   fprintf(fp,"\n原网的阻抗阵:\n",node0);
   for(i=0;i<Nnode;i++)
   {
	   for(j=0;j<Nnode;j++)
		   fprintf(fp,"%f+j%f\t",Z0_real[i][j],Z0_imag[i][j]);
	   fprintf(fp,"\n");
   }
   fprintf(fp,"\n在节点%d追加一条树支后的阻抗阵:\n",node0);
   for(i=0;i<Nnode+1;i++)
   {
	   for(j=0;j<Nnode+1;j++)
		   fprintf(fp,"%f+j%f\t",Znew_real[i][j],Znew_imag[i][j]);
	   fprintf(fp,"\n");
   }
//追加连支后的阻抗阵
	Zll_real=Z0_real[node_i-1][node_i-1]+Z0_real[node_j-1][node_j-1]-2*Z0_real[node_i-1][node_j-1]+Radd_br;
	Zll_imag=Z0_imag[node_i-1][node_i-1]+Z0_imag[node_j-1][node_j-1]-2*Z0_imag[node_i-1][node_j-1]+Xadd_br;
    for(k=0;k<Nnode;k++)
	{
		Zl_real[k]=Z0_real[k][node_i-1]-Z0_real[k][node_j-1];
		Zl_imag[k]=Z0_imag[k][node_i-1]-Z0_imag[k][node_j-1];
	}
	for(i=0;i<Nnode;i++)
		for(j=0;j<Nnode;j++)
			K_real[i][j]=Zl_real[i]*Zl_real[j];
	for(i=0;i<Nnode;i++)
		for(j=0;j<Nnode;j++)
			K_imag[i][j]=Zl_imag[i]*Zl_imag[j];
	for(i=0;i<Nnode;i++)
		for(j=0;j<Nnode;j++)
		{
			Znew_real[i][j]=Z0_real[i][j]-K_real[i][j]/Zll_real;
			Znew_imag[i][j]=Z0_imag[i][j]-K_imag[i][j]/Zll_real;
		}
//将结果保存到"运行结果.txt"中
   fprintf(fp,"\n在节点%d和节点%d之间追加一条连支后的阻抗阵:\n",node_i,node_j);
   for(i=0;i<Nnode;i++)
   {
	   for(j=0;j<Nnode;j++)
		   fprintf(fp,"%f+j%f\t",Znew_real[i][j],Znew_imag[i][j]);
	   fprintf(fp,"\n");
   }
   fclose(fp);
}
//输出导纳阵
/*	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Y0_real[i][j]);
		printf("\n");
	}
    printf("\n");
	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Y0_imag[i][j]);
		printf("\n");
	}
	printf("\n");
//输出阻抗阵
	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Z0_real[i][j]);
		printf("\n");
	}
    printf("\n");
	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Z0_imag[i][j]);
		printf("\n");
	}
	printf("\n");

//输出追加连支后的阻抗阵
    printf("输出追加连支后的阻抗阵\n");
	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Znew_real[i][j]);
		printf("\n");
	}
    printf("\n");
	for(i=0;i<Nnode;i++)
	{
		for(j=0;j<Nnode;j++)
			printf("%f\t",Znew_imag[i][j]);
		printf("\n");
	}*/