lizhi.cpp

基于异步发电机的配电网潮流计算程序,适合于分布式发电对配网的影响分析

#include "stdafx.h"
#include "lizhi.h"
#include "math.h"
#include <iomanip.h>
#include "fstream.h"

CDistributionFlow::CDistributionFlow()
{
	N=0;
	LL=0;
	Tree=0;
	Loop=0;
	NW=0;
	PFN=0.0;
	PFMAX=0.0;
	DONE=0;
	VB=0.0;
	SB=0.0;
}

CDistributionFlow::~CDistributionFlow()
{
	m_Branch.RemoveAll();
	m_Node.RemoveAll();
	m_IT.RemoveAll();
	m_IL.RemoveAll();
	m_TempIT.RemoveAll();
	m_TempIL.RemoveAll();
	m_FatherNode.RemoveAll();
	m_Wind.RemoveAll();
}

void CDistributionFlow::Run()
{
	ReadData();
	GetTree();
	Calculate();
}

void CDistributionFlow::ReadData()
{
	int i;

	ifstream infile;
	infile.open("data30.txt",ios::nocreate);

	if(infile.fail())
	{
		cout<<"Error,Can't open this file!"<<endl;
	}

	infile>>N>>LL>>NW>>VB>>SB>>PFN>>PFMAX;//节点数，支路数,风机节点数，基准电压，基准功率

	m_Node.SetSize(N+1);
	m_FatherNode.SetSize(N+1);
	m_Branch.SetSize(LL+1);
	m_IT.SetSize(LL+1);
	m_IL.SetSize(LL+1);
	m_TempIT.SetSize(LL+1);
	m_TempIL.SetSize(LL+1);
	m_Wind.SetSize(NW+1);

	for(i=1;i<=N;i++)
	{
		m_FatherNode[i]=0;
	}

	for(i=1;i<=LL;i++)
	{
		m_TempIT[i]=0;
		m_TempIL[i]=0;
		m_IT[i].IT1=0;
		m_IT[i].IT2=0;
		m_IL[i].IL1=0;
		m_IL[i].IL2=0;
	}

	for(i=1;i<=LL;i++)
	{
		infile>>m_Branch[i].LJ>>m_Branch[i].LK>>m_Branch[i].R>>m_Branch[i].X;
	}

		for(i=1;i<=N;i++)
	{
		infile>>m_Node[i].P>>m_Node[i].Q;
	}

	for(i=1;i<=N;i++)
	{
		infile>>m_Node[i].E>>m_Node[i].F>>m_Node[i].V
			>>m_Node[i].PVFLAG>>m_Node[i].ISPHNODE;
	}

	for(i=1;i<=NW;i++)
	{
		infile>>m_Wind[i].NUM>>m_Wind[i].PN>>m_Wind[i].VN>>m_Wind[i].PE>>m_Wind[i].R1>>m_Wind[i].R2>>m_Wind[i].X1
			>>m_Wind[i].X2>>m_Wind[i].Xm>>m_Wind[i].QNunit>>m_Wind[i].VCN;
	}
	
	for(i=1;i<=NW;i++)
	{
		m_Wind[i].CN=0;
		m_Wind[i].OLDCN=0;
		m_Wind[i].NMAX=0;
		m_Wind[i].s=0.0;
		m_Wind[i].PF=0.0;
		m_Wind[i].Qe=0.0;
		m_Wind[i].Qc=0.0;
		m_Wind[i].Xf=m_Wind[i].X1+m_Wind[i].X2;
	}
	
	for(i=1;i<=NW;i++)
	{
		m_Wind[i].s=m_Wind[i].R2*(pow(m_Wind[i].VN,2)-sqrt(pow(m_Wind[i].VN,4)-4*pow(m_Wind[i].Xf,2)*pow(m_Wind[i].PN,2)))
			/(2*m_Wind[i].PN*pow(m_Wind[i].Xf,2));
		m_Wind[i].Qe=(pow(m_Wind[i].R2,2)+m_Wind[i].Xf*(m_Wind[i].Xm+m_Wind[i].Xf)*pow(m_Wind[i].s,2))*m_Wind[i].PN
			/(m_Wind[i].s*m_Wind[i].R2*m_Wind[i].Xm);
		m_Wind[i].PF=m_Wind[i].PN/(sqrt(pow(m_Wind[i].PN,2)+pow(m_Wind[i].Qe,2)));
		m_Wind[i].Qc=m_Wind[i].PN*(sqrt(1/pow(m_Wind[i].PF,2)-1)-sqrt(1/pow(PFMAX,2)-1));
		m_Wind[i].NMAX=int(m_Wind[i].Qc/m_Wind[i].QNunit)+1;
	}
	infile.close();
}

void CDistributionFlow::GetTree()
{
	int i,j;
	int J,K;
	int TempNode1,TempNode2;
	int *NADJ=new int[N+1];
	int (*NearBranch)[MAX_BRANCH+1]=new int[N+1][MAX_BRANCH+1];
	
	for(i=1;i<=N;i++)
	{
		NADJ[i]=0;
		for(j=1;j<=MAX_BRANCH;j++)
		{
			NearBranch[i][j]=0;
		}
	}

	for(i=1;i<=LL;i++)
	{
		J=m_Branch[i].LJ;
		K=m_Branch[i].LK;

		NADJ[J]+=1;
		NearBranch[J][NADJ[J]]=K;
		NADJ[K]+=1;
		NearBranch[K][NADJ[K]]=J;
	}

	TempNode1=1;
	TempNode2=0;
	m_FatherNode[TempNode1]=-1;

	while((Tree+Loop)<LL)
	{
		if(NADJ[TempNode1]>0)
		{
			i=1;
			while(NearBranch[TempNode1][i]==0)
			{
				i++;
			}
			TempNode2=NearBranch[TempNode1][i];
			if(m_FatherNode[TempNode2]==0)
			{
				m_FatherNode[TempNode2]=TempNode1;
				m_IT[Tree+1].IT1=TempNode1;
				m_IT[Tree+1].IT2=TempNode2;
				NearBranch[TempNode1][i]=0;
				NADJ[TempNode1]-=1;
				i=1;
				while(NearBranch[TempNode2][i]!=TempNode1)
				{
					i++;
				}
				NearBranch[TempNode2][i]=0;
				NADJ[TempNode2]-=1;
				TempNode1=TempNode2;
				Tree+=1;
			}
			else
			{
				m_IL[Loop+1].IL1=TempNode1;
				m_IL[Loop+1].IL2=TempNode2;
				NearBranch[TempNode1][i]=0;
				NADJ[TempNode1]-=1;
				i=1;
				while(NearBranch[TempNode2][i]!=TempNode1)
				{
					i++;
				}
				NearBranch[TempNode2][i]=0;
				NADJ[TempNode2]-=1;
				Loop+=1;
			}
			
		}
		else
		{
			TempNode2=m_FatherNode[TempNode1];
			TempNode1=TempNode2;
		}		
	}

	for(i=1;i<=Tree;i++)
	{
		j=1;
		while((m_IT[i].IT1!=(m_Branch[j].LJ)||m_IT[i].IT2!=(m_Branch[j].LK))&&(m_IT[i].IT1!=(m_Branch[j].LK)||m_IT[i].IT2!=(m_Branch[j].LJ)))
		{
			j++;
		}
		m_TempIT[i]=j;
	}

	for(i=1;i<=Loop;i++)
	{
		j=1;
		while((m_IL[i].IL1!=(m_Branch[j].LJ)||m_IL[i].IL2!=(m_Branch[j].LK))&&(m_IL[i].IL1!=(m_Branch[j].LK)||m_IL[i].IL2!=(m_Branch[j].LJ)))
		{
			j++;
		}
		m_TempIL[i]=j;
	}

	delete [] NADJ;
	delete [] NearBranch;
}

void CDistributionFlow::Calculate()
{
	int i,j,m,mm,J,K,Iteration,TempNode1;
	double SumE,SumF,SumR,SumX,PLoop,QLoop,AAA,BBB;
	double wucha,jingdu,TempE,TempF,TempV,xx,xx1;
	double *PNode=new double[N+1];
	double *QNode=new double[N+1];
	double *PLL=new double[Tree+1];
	double *QLL=new double[Tree+1];
	double *ELL=new double[Tree+1];
	double *FLL=new double[Tree+1];
	double *Error=new double[N+1];
	
	for(i=1;i<=Tree;i++)
	{
		PLL[i]=0;
		QLL[i]=0;
		ELL[i]=0;
		FLL[i]=0;
		Error[i]=0;
	}

	for(i=1;i<=NW;i++)
	{
		mm=m_Wind[i].NUM;

		m_Wind[i].s=m_Wind[i].R2*(pow(m_Node[mm].V*VB,2)-sqrt(pow(m_Node[mm].V*VB,4)-4*pow(m_Wind[i].Xf,2)*pow(m_Wind[i].PE,2)))
			/(2*m_Wind[i].PE*pow(m_Wind[i].Xf,2));
		m_Wind[i].Qe=(pow(m_Wind[i].R2,2)+m_Wind[i].Xf*(m_Wind[i].Xm+m_Wind[i].Xf)*pow(m_Wind[i].s,2))*m_Wind[i].PE
			/(m_Wind[i].s*m_Wind[i].R2*m_Wind[i].Xm);
		
		m_Wind[i].Qc=m_Wind[i].CN*m_Wind[i].QNunit*(pow(m_Node[mm].V*VB,2)/pow(m_Wind[i].VCN,2));
		m_Node[mm].P=-m_Wind[i].PE/SB;
		m_Node[mm].Q=(-m_Wind[i].Qc+m_Wind[i].Qe)/SB;
	}

	wucha=0.005;
	jingdu=0.00001;
	Iteration=0;
	DONE=0;
	while(DONE!=1)//
	{
		//电流回推，由树干的顶端反向递推各树支的潮流，先不计连支
		for(i=2;i<=N;i++)
		{
			xx=pow(m_Node[i].E,2)+pow(m_Node[i].F,2);
			PNode[i]=(m_Node[i].P*m_Node[i].E+m_Node[i].Q*m_Node[i].F)/xx;
			QNode[i]=(m_Node[i].P*m_Node[i].F-m_Node[i].Q*m_Node[i].E)/xx;
		}

		for(i=Tree;i>=1;i--)
		{
			J=m_IT[i].IT1;
			K=m_IT[i].IT2;

			PLL[i]=PNode[K];
			QLL[i]=QNode[K];

			PNode[J]+=PLL[i];
			QNode[J]+=QLL[i];
		}
		
		//计算连支，修正潮流
		for(i=1;i<=Loop;i++)
		{
			J=m_IL[i].IL1;
			K=m_IL[i].IL2;
			TempNode1=J;
			SumE=0.0;
			SumF=0.0;
			SumR=m_Branch[m_TempIL[i]].R;
			SumX=m_Branch[m_TempIL[i]].X;
			PLoop=0.0;
			QLoop=0.0;
			//求连支电流PLoop,QLoop
			while(TempNode1!=K)
			{
				j=1;
				while(m_IT[j].IT2!=TempNode1)
				{
					j++;
				}
				m=m_TempIT[j];
				SumE-=PLL[j]*m_Branch[m].R-QLL[j]*m_Branch[m].X;
				SumF-=PLL[j]*m_Branch[m].X+QLL[j]*m_Branch[m].R;
				SumR+=m_Branch[m].R;
				SumX+=m_Branch[m].X;
				TempNode1=m_FatherNode[TempNode1];
			}
			xx1=pow(SumE,2)+pow(SumF,2);
			PLoop=(SumE*SumR+SumF*SumX)/xx1;
			QLoop=(-SumE*SumX+SumF*SumR)/xx1;

			//对支路电流修正
			TempNode1=J;
			while(TempNode1!=K)
			{
				j=1;
				while(m_IT[j].IT2!=TempNode1)
				{
					j++;
				}
				PLL[j]+=PLoop;
				QLL[j]+=QLoop;
				TempNode1=m_FatherNode[TempNode1];
			}
		}

		//电压前推
		wucha=0.0;
		for(i=1;i<=Tree;i++)
		{
			J=m_IT[i].IT1;
			K=m_IT[i].IT2;
			m=m_TempIT[i];
			ELL[i]=PLL[i]*m_Branch[m].R-QLL[i]*m_Branch[m].X;
			FLL[i]=PLL[i]*m_Branch[m].X+QLL[i]*m_Branch[m].R;
			TempE=m_Node[K].E;
			TempF=m_Node[K].F;
			TempV=m_Node[K].V;
			m_Node[K].E=m_Node[J].E-ELL[i];
			m_Node[K].F=m_Node[J].F-FLL[i];
			m_Node[K].V=sqrt(pow(m_Node[K].E,2)+pow(m_Node[K].F,2));
			
			Error[K]=fabs(TempE-m_Node[K].E);
			if(fabs(TempF-m_Node[K].F)>Error[K])
			{
				Error[K]=fabs(TempF-m_Node[K].F);
			}
			if(fabs(TempV-m_Node[K].V)>Error[K])
			{
				Error[K]=fabs(TempV-m_Node[K].V);
			}
			if(Error[K]>wucha)
			{
				wucha=Error[K];
			}
		}
		if(wucha>jingdu)
		{
			DONE=0;
		}
		else
		{
			DONE=1;
		}

		for(i=1;i<=NW;i++)
		{
			m_Wind[i].OLDCN=m_Wind[i].CN;
		}

		for(i=1;i<=NW;i++)
		{
			mm=m_Wind[i].NUM;
			m_Wind[i].s=m_Wind[i].R2*(pow(m_Node[mm].V*VB,2)-sqrt(pow(m_Node[mm].V*VB,4)-4*pow(m_Wind[i].Xf,2)*pow(m_Wind[i].PE,2)))
				/(2*m_Wind[i].PE*pow(m_Wind[i].Xf,2));
			m_Wind[i].Qe=(pow(m_Wind[i].R2,2)+m_Wind[i].Xf*(m_Wind[i].Xm+m_Wind[i].Xf)*pow(m_Wind[i].s,2))*m_Wind[i].PE
				/(m_Wind[i].s*m_Wind[i].R2*m_Wind[i].Xm);
			m_Wind[i].Qc=m_Wind[i].CN*m_Wind[i].QNunit*(pow(m_Node[mm].V*VB,2)/pow(m_Wind[i].VCN,2));
			m_Node[mm].P=-m_Wind[i].PE/SB;
			m_Node[mm].Q=(-m_Wind[i].Qc+m_Wind[i].Qe)/SB;
			m_Wind[i].PF=fabs(m_Node[mm].P)/(sqrt(pow(fabs(m_Node[mm].P),2)+pow(m_Node[mm].Q,2)));
			while(m_Wind[i].PF>PFN)
			{
				m_Wind[i].CN=m_Wind[i].CN-1;
				m_Wind[i].Qc=m_Wind[i].CN*m_Wind[i].QNunit*(pow(m_Node[mm].V*VB,2)/pow(m_Wind[i].VCN,2));
				m_Node[mm].Q=(-m_Wind[i].Qc+m_Wind[i].Qe)/SB;
				m_Wind[i].PF=fabs(m_Node[mm].P)/(sqrt(pow(fabs(m_Node[mm].P),2)+pow(m_Node[mm].Q,2)));
			}
			while(m_Wind[i].PF<=PFN)
			{
				BBB=m_Wind[i].CN;
				m_Wind[i].CN=m_Wind[i].CN+1;
				m_Wind[i].Qc=m_Wind[i].CN*m_Wind[i].QNunit*(pow(m_Node[mm].V*VB,2)/pow(m_Wind[i].VCN,2));
				m_Node[mm].Q=(-m_Wind[i].Qc+m_Wind[i].Qe)/SB;
				m_Wind[i].PF=fabs(m_Node[mm].P)/(sqrt(pow(fabs(m_Node[mm].P),2)+pow(m_Node[mm].Q,2)));
				AAA=m_Wind[i].PF;
			}
		}

		for(i=1;i<=NW;i++)
		{
			if(m_Wind[i].CN>m_Wind[i].NMAX)
			{
				ofstream outfile;
				outfile.open("error.txt");
				outfile<<"要求补偿容量超过限制，出错！"<<endl;
				outfile.close();
				return;
			}
		}

		for(i=1;i<=NW;i++)//由于不满足功率因数要求，电容器投切组数发生变化，置收敛标志为0
		{
			if(m_Wind[i].CN!=m_Wind[i].OLDCN)
			{
				DONE=0;
			}
		}
		
		Iteration++;
	}

	ofstream outfile;
	outfile.open("result100.txt");
	outfile<<"迭代次数：  "<<Iteration<<endl;
	outfile<<"最终误差：  "<<wucha<<endl;
	outfile<<"节点号"<<"   "<<"风电机输出功率"<<"   "<<"投入并联电容台数"<<endl;
	for(i=1;i<=NW;i++)
	{
		outfile<<m_Wind[i].NUM<<"   "<<m_Wind[i].PE<<"   "<<m_Wind[i].CN<<endl;
	}
	outfile<<endl;
	outfile<<"节点号"<<"   "<<"节点电压实部"<<"   "<<"节点电压虚部"<<"   "<<"节点电压幅值"<<endl;
	for(i=1;i<=N;i++)
	{
		outfile<<i<<"   "<<m_Node[i].E<<"   "<<m_Node[i].F<<"   "<<m_Node[i].V<<endl;
	}
	outfile.close();

	delete [] PNode;
	delete [] QNode;
	delete [] PLL;
	delete [] QLL;
	delete [] ELL;
	delete [] FLL;
	delete [] Error;
}