flow080605.cpp

电力系统潮流

// flow.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream.h>
#include <fstream.h>
#include "NEquation.h"
#include "math.h"
#include "config.h"

void test()
{
	NEquation ob1;
	ob1.SetSize(2);
	ob1.Data(0,0)=1;
	ob1.Data(0,1)=1;
	ob1.Data(1,0)=2;
	ob1.Data(1,1)=4;

	ob1.Value(0)=12;
	ob1.Value(1)=12;
	ob1.Run();

	printf("%5.3f,%5.3f\n",ob1.Value(0),ob1.Value(1));
}


//获取节点和线路数据
void GetData()
{
	FILE *fp;
	int i;

	if((fp=fopen("D:\\潮流计算080603\\flow\\data\\data.txt","r"))==NULL)
	{ 
		printf("Can not open the file named 'data.txt' \n");
		return;
	}

	for(i=0;i<=Bus_Num-1;i++)
	{
		fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].Voltage,&gBus[i].Phase,
			&gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type);
	}

	for(i=0;i<=Line_Num-1;i++)
	{
		fscanf(fp,"%d,%d,%d,%f,%f,%f",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J,
			&gLine[i].R,&gLine[i].X,&gLine[i].B);	
	}
	
/*	for(i=0;i<=Bus_Num-1;i++)
	{
		fscanf(fp,"%d,%d,%f,%f",&gShunt[i].Num,&gShunt[i].NumI,&gShunt[i].G,&gShunt[i].B);
	}
*/
	fclose(fp);
}
//计算节点导纳矩阵
void GetYMatrix()
{
	int i,j,bus1,bus2;
	float r,x,d,g,b;
	FILE *fp;

	for(i=0;i<=Bus_Num-1;i++)
	{
		for(j=0;j<=Bus_Num-1;j++)
		{
			gY_G[i][j]=0;
			gY_B[i][j]=0;
		}
	}

	for(i=0; i<=Line_Num-1; i++)
	{
		bus1=gLine[i].No_I-1;
		bus2=gLine[i].No_J-1;
		r=gLine[i].R;
		x=gLine[i].X;
		d=r*r+x*x;
		if(d>1e-6)
		{
			g=r/d;
			b=-x/d;
		}
		else
		{
			g=0;
			b=0;
		}
	
		gY_G[bus1][bus1]+=g;
		gY_G[bus2][bus2]+=g;
		gY_G[bus1][bus2]-=g;
		gY_G[bus2][bus1]-=g;
		gY_B[bus1][bus1]+=b+gLine[i].B/2;
		gY_B[bus2][bus2]+=b+gLine[i].B/2;
		gY_B[bus1][bus2]-=b;
		gY_B[bus2][bus1]-=b;
	}

	// output the Y matrix 
	if((fp=fopen("D:\\潮流计算080603\\flow\\data\\ymatrix.txt","w"))==NULL)
	{
		printf("Can not open the file named 'ymatrix.txt' \n");
		return ;
	}

	fprintf(fp,"---Y Matrix---\n");
	printf("			the Y matrix is calcalated below \n");
	for(i=0;i<=Bus_Num-1;i++)
	{
		int k;
		for(j=0;j<=Bus_Num-1;j++)
		{
			fprintf(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);
		
				
			printf("%7.3f,%7.3f\t",gY_G[i][j],gY_B[i][j]);
			for(k=0;k<=Bus_Num-1;k++)
			{
				if ((k+1)%5) break;
				else printf("\n");
				
			}
		}
	
	}

				printf("\n");
				printf("\n");
	fclose(fp);	
	
}
//设电压初值
void SetInitial()
{
	int i;

	for(i=0;i<=Bus_Num-1;i++)
	{
		gf[i]=gBus[i].Voltage*sin(gBus[i].Phase);
		ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);
	}
}
//求取功率不平衡量
void GetUnbalance()
{	
	float Bus_P[Bus_Num-1],Bus_Q[Bus_Num-1];
	int i,j;

	for (i=1;i<=Bus_Num-1;i++)
	{
		Bus_P[i-1]=0;
		Bus_Q[i-1]=0;
		for (j=0;j<=Bus_Num-1;j++)
		{
			Bus_P[i-1]+=ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
			Bus_Q[i-1]+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
		
		}
	}
	printf("	the unbalanced value of Power is calcalated below\n");
	for (i=1;i<=Bus_Num-1;i++)
	{
		
		gDelta_P[i-1]=gBus[i].GenP-gBus[i].LoadP-Bus_P[i-1];
		gDelta_Q[i-1]=gBus[i].GenQ-gBus[i].LoadQ-Bus_Q[i-1];
	
		printf("%f	 %f\n",gDelta_P[i-1],gDelta_Q[i-1]);
	}
		printf("\n");
		printf("\n");
}



//计算雅可比矩阵
void GetJaccobi()
{
	float JH[2*(Bus_Num-1)][2*(Bus_Num-1)];
	float JN[2*(Bus_Num-1)][2*(Bus_Num-1)];
	float JJ[2*(Bus_Num-1)][2*(Bus_Num-1)];
	float JL[2*(Bus_Num-1)][2*(Bus_Num-1)];
	float Ja[Bus_Num-1][Bus_Num-1];
	float Jb[Bus_Num-1][Bus_Num-1];
	int i,j,m=0;

	for(i=0;i<2*(Bus_Num-1);i++)
	{
		for(j=0;j<2*(Bus_Num-1);j++)
		{
			JH[i][j]=0;
			JN[i][j]=0;
			JJ[i][j]=0;
			JL[i][j]=0;
			Ja[i][j]=0;
			Jb[i][j]=0;
		}
	}
	printf ("		the Jaccobi matrix is calculated below\n");
	for(i=0;i<2*(Bus_Num-1);i++)
	{
		for(j=0;j<2*(Bus_Num-1);j++)
		{
			if (i!=j)
			{
			JH[i][j]=-gY_B[i][j]*ge[i]+gY_G[i][j]*gf[i];
			JN[i][j]=gY_G[i][j]*ge[i]+gY_B[i][j]*gf[i];
			JJ[i][j]=-JN[i][j];
			JL[i][j]=JH[i][j];
			Ja[i][j]+=gY_G[i][i]*ge[i]-gY_B[i][i]*gf[i]+gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j];
			Jb[i][j]+=gY_B[i][i]*ge[i]+gY_G[i][i]*gf[i]+gY_B[i][j]*ge[j]+gY_G[i][j]*gf[j];
										
					printf("JMN(%d,%d):	%7.3f   %7.3f\n",i+1,j+1,JH[i][j],JN[i][j]);
					printf("JJL(%d,%d):	%7.3f   %7.3f\n",i+1,j+1,JJ[i][j],JL[i][j]);
				
				
			}
			else 
			{
			Ja[i][j]+=gY_G[i][i]*ge[i]-gY_B[i][i]*gf[i]+gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j];
			Jb[i][j]+=gY_B[i][i]*ge[i]+gY_G[i][i]*gf[i]+gY_B[i][j]*ge[j]+gY_G[i][j]*gf[j];
			JH[i][j]=-gY_B[i][j]*ge[i]+gY_G[i][j]*gf[i]+Jb[i][j];
			JN[i][j]=gY_G[i][j]*ge[i]+gY_B[i][j]*gf[i]+Ja[i][j];
			JJ[i][j]=-JN[i][j]+Ja[i][j];
			JL[i][j]=JH[i][j]-Jb[i][j];

					printf("JMN(%d,%d):	%7.3f   %7.3f\n",i+1,j+1,JH[i][j],JN[i][j]);
					printf("JJL(%d,%d):	%7.3f   %7.3f\n",i+1,j+1,JJ[i][j],JL[i][j]);
					
			
			}

			
		}
	}


}
/*
//计算电压不平衡量
void GetRevised()
{		
	int i,j;
	

	text();
	ob1.SetSize(1);
	ob1.Data(0,0)=1;

}
//计算电压新值
void GetNewValue()
{
	int i;
	
	for(i=0;i<=Bus_Num-2;i++)
	{
		gf[i]=gf[i]+gDelta_f[i];
		ge[i]=ge[i]+gDelta_e[i];
	}
}

*/
int main(int argc, char* argv[])
{
	int i,Count_Num;
	float minValue;

//	test();
	
	GetData();
	GetYMatrix();
	SetInitial();
	
	GetUnbalance();
	GetJaccobi();

/*	for(Count_Num=0;Count_Num<=100;Count_Num++)
	{
		GetUnbalance();
		GetJaccobi();
		GetRevised();
		GetNewValue();

		minValue=gDelta_fe[0];
		for(i=0;i<=2*(Bus_Num-1)-1;i++)
		{
			if(minValue>gDelta_fe[i])
			{
				minValue=gDelta_fe[i];
			}
		}
		if(minValue<Precision)
		{
			break;
		}
	}*/


//	printf("%d",Count_Num);
	return 0;
}