flow.cpp

该程序是一个含直流输电系统的潮流计算程序 程序运行需要包含Boost库文件,需要在VS2005环境以上(VC6.0不完全符合C++标准,可能有些变量作用域问题,无法通过调试) 潮流程序计算主要参考

#if defined(_MSC_VER) && _MSC_VER <= 1200
#pragma warning ( disable: 4503 4786 )
#endif

#ifdef NDEBUG
#pragma warning ( disable: 4267 )
#endif

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>

#ifndef NDEBUG
#include <cstdlib>
#endif

#include <boost/foreach.hpp>
#include <boost/integer_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_utility.hpp>
#include <boost/io/ios_state.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>

#ifndef NDEBUG
#include <boost/numeric/ublas/io.hpp>
#endif

#if !( defined(_MSC_VER) && _MSC_VER <= 1200 )
#include <boost/cast.hpp>
using boost::numeric_cast;
#endif

#include "Flow.h"

// 2009.3.29 Ignore the comments in the input data file.
template < class IStream >
void IgnoreComments( IStream& in )
{
  typedef typename IStream::char_type Char;
  static const Char beg_tag('/');
  static const Char end_tag('\n');

  if( in.eof() ) return;

  while ( IStream::sentry(in) && in.peek() == beg_tag ) {
    in.ignore(boost::integer_traits<int>::const_max,end_tag);
  }

  return;
}

// 2009.4.4 Count the injected edges resulted from the elimination of a vertex.
template < class Vertex, class Graph >
std::size_t CountInjectedEdges( Vertex const& vertex, Graph const& graph )
{
  typedef boost::graph_traits<Graph>::adjacency_iterator AdjVertexIter;
  typedef boost::graph_traits<Graph>::edges_size_type EdgeSizeType;

  EdgeSizeType inj_edges = 0;

  AdjVertexIter adj_ver_beg, adj_ver_end;
  boost::tie(adj_ver_beg,adj_ver_end) = boost::adjacent_vertices(vertex,graph);
  for( ; adj_ver_beg != adj_ver_end; ++adj_ver_beg ) {
    for( AdjVertexIter adj_ver_next = adj_ver_beg; ++adj_ver_next != adj_ver_end; ) {
      if( !boost::edge(*adj_ver_beg,*adj_ver_next,graph).second ) {
        ++inj_edges;    
      }
    }
  }

  return inj_edges;
}

// 2009.4.4 Eliminate the just labeled vertex and add injected edges into the graph.
template < class Vertex, class Graph >
void EliminateVertex( Vertex const& vertex, Graph& graph )
{
  typedef boost::graph_traits<Graph>::adjacency_iterator AdjVertexIter;

  AdjVertexIter adj_ver_beg, adj_ver_end;
  boost::tie(adj_ver_beg,adj_ver_end) = boost::adjacent_vertices(vertex,graph);
  for( ; adj_ver_beg != adj_ver_end; ++adj_ver_beg ) {
    for( AdjVertexIter adj_ver_next = adj_ver_beg; ++adj_ver_next != adj_ver_end; ) {
      if( !boost::edge(*adj_ver_beg,*adj_ver_next,graph).second ) {
        boost::add_edge(*adj_ver_beg,*adj_ver_next,graph);   
      }
    }
  }

  boost::clear_vertex(vertex,graph);

  return;
}

//---构造函数1-----------------------------------------------------------------------------
Flow::Flow()
: dc_correction_perm(0), 
  Y(), X(), B() 
{
	MaxIterNum=20;
	eps=0.00001;
	SumOfPV=0;
	bPFconverge=false;
	bBackupedJacobian=false;
	bInitialized=false;
}
//----构造函数2----------------------------------------------------------------------------
Flow::Flow(int MNodeNum,int MBranchNum,int MMaxIterNum,double MEps)
: dc_correction_perm(0), 
  Y(MNodeNum,MNodeNum),X(MNodeNum,MNodeNum)
{
	NodeNum=MNodeNum;
	BranchNum=MBranchNum;
	MaxIterNum=MMaxIterNum;
	eps=MEps;
	
	bInitialized=ACInitialize();
	
}
//--------------------------------------------------------------------------------
Flow::~Flow()
{
	if(bInitialized)
	{
		delete [] pNodeData;
		delete [] pBranchData;
		delete [] OptimisedNodeTable;
		delete [] OriginalNodeTable;
		delete [] Phase;
		delete [] U;
		delete [] pBranchLoss;
		
		delete [] Pij;
		delete [] Pji;
		delete [] Qij;
		delete [] Qji;
		delete [] Iij;
		delete [] Iji;
		delete [] NodalPowerP;
		delete [] NodalPowerQ;
//		delete [] ShuntLoad;
	}
	if(bBackupedJacobian) delete Jacobian;
}

//--------------------------------------------------------------------------------
//读入数据
//--------------------------------------------------------------------------------
//通过节点序号来设置节点内容
//参数: 序号,节点编号,节点类型,有功/相角,无功/电压
//--------------------------------------------------------------------------------
bool Flow::ReadData( char* in, char* out )
{
  infile.open(in/*,ios::nocreate*/); // 2008.3.27 Nocreate is not a standard stream flag.	
	if( infile.fail() ) {
		std::cerr << "Opening input file failed!" << std::endl;
		return false;
	}

  // Read in the controlling parameters of load flow calculation.
  IgnoreComments(infile);
  infile >> MaxIterNum;
  infile >> eps;
  infile >> BaseS;

  // Read in the AC system data.
  if ( !ReadACData() ) {
    infile.close();
    return false;
  }

  // Read in the DC system data.
  if ( !ReadDCData() ) {
    infile.close();
    return false;
  }

	infile.close();

	outfile.open(out);
	if( outfile.fail() ) {
		std::cerr << "Open output file failed!" << std::endl;
		return false;
	}

	outfile << setiosflags(ios::right) << setprecision(6);
	outfile << "交流系统节点数    " << setw(10) << NodeNum << '\n'
          << "交流系统支路数    " << setw(10) << BranchNum << '\n'
          << "直流系统节点数    " << setw(10) << dc_bus_count << '\n'
          << "直流系统支路数    " << setw(10) << dc_branch_count << '\n'
          << "计算收敛精度      " << setw(10) << eps << '\n'
          << "功率基值(MVA/MW)  " << setw(10) << BaseS << '\n';

  bool bSuccess = PreHandData();

	return bSuccess;	   
}

bool Flow::ReadACData()
{
  // Read in AC system information.
  IgnoreComments(infile);
	infile>>NodeNum;
	infile>>BranchNum;

	if ( !ACInitialize() )
		return false;		
	
  { // Read in AC bus data.
    for( int i = 0; i < NodeNum; ++i ) {
      IgnoreComments(infile);
      infile >> pNodeData[i].Num;
      infile >> pNodeData[i].Name;
      infile >> pNodeData[i].AreaNum;
      infile >> pNodeData[i].LossZoneNum;
      infile >> pNodeData[i].NodeType;
      infile >> pNodeData[i].Fv;
      infile >> pNodeData[i].Fa;
      infile >> pNodeData[i].Lp;
      infile >> pNodeData[i].Lq;
      infile >> pNodeData[i].Gp;
      infile >> pNodeData[i].Gq;
      infile >> pNodeData[i].BaseVolt;
      infile >> pNodeData[i].DesireVolt;
      infile >> pNodeData[i].MaxVarOrV;
      infile >> pNodeData[i].MinVarOrV;
      infile >> pNodeData[i].ShuntG;
      infile >> pNodeData[i].ShuntB;
      infile >> pNodeData[i].RemoteBusNum;
    }
  }

  { // Read in AC branch data.
    for( int i = 0; i < BranchNum; ++i ) {
      IgnoreComments(infile);
      infile >> pBranchData[i].NodeI;
      infile >> pBranchData[i].NodeJ;
      infile >> pBranchData[i].LfArea;
      infile >> pBranchData[i].LossZone;
      infile >> pBranchData[i].Circuit;
      infile >> pBranchData[i].BranchType;
      infile >> pBranchData[i].R;
      infile >> pBranchData[i].X;
      infile >> pBranchData[i].B;
      infile >> pBranchData[i].RateNo1;
      infile >> pBranchData[i].RateNo2;
      infile >> pBranchData[i].RateNo3;
      infile >> pBranchData[i].ControlBusNum;
      infile >> pBranchData[i].ControlBusSide;
      infile >> pBranchData[i].TRatio;
      infile >> pBranchData[i].TAngle;
      infile >> pBranchData[i].MaxTap;
      infile >> pBranchData[i].MinTap;
      infile >> pBranchData[i].StepSize;
      infile >> pBranchData[i].MaxLimit;
      infile >> pBranchData[i].MinLimit;
    }
  }

  return true;
}

bool Flow::ReadDCData()
{
  // Read in DC system information.
  IgnoreComments(infile);
  infile >> dc_bus_count;
  infile >> dc_branch_count;

  if ( !DCInitialize() ) {
    return false;
  }

  {  // Read DC bus data one by one.
    for ( int i = 0; i < dc_bus_count; ++i ) {
      IgnoreComments(infile);
      infile >> dc_buses[i];

      int id = dc_buses[i].GetID();
      id2index[id] = i;
      index2id[i] = id;
    }
  }

  { // Read DC branch data one by one.
    for ( int i = 0; i < dc_branch_count; ++i ) {
      IgnoreComments(infile);
      infile >> dc_branches[i];
    }
  }

  return true;
}

//得到节点和支路的数据后，开始初始化
bool Flow::ACInitialize()
{
	int bSuccess;
	
	pNodeData=new ACBus[NodeNum];			//建立节点数据指针数组
	pBranchData=new ACBranch[BranchNum];
	OptimisedNodeTable=new int[NodeNum];
	OriginalNodeTable = new int[NodeNum];
	Phase = new double[NodeNum];
	U = new double[NodeNum];
	pBranchLoss = new Complex[BranchNum];
	NodalPowerP=new double [NodeNum];
	NodalPowerQ=new double [NodeNum];
	
/*	//-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
	if(NodeNum==3141)
	{
		ShuntLoad=new CShuntLoad[NodeNum];
		ifstream ShuntFile;
		ShuntFile.open(".\\输入数据\\并联导纳负荷.dat",ios::nocreate);	
		if(ShuntFile.fail())
		{
			std::cerr<<"open ShuntLoad file error!"<<std::endl;
			return false;
		}
		for(int i=0; i<NodeNum ; i++)
		{
			ShuntFile>>ShuntLoad[i].i>>ShuntLoad[i].ShuntP>>ShuntLoad[i].ShuntQ;
			ShuntLoad[i].ShuntP/=BaseS;
			ShuntLoad[i].ShuntQ/=BaseS;
		}
		ShuntFile.close();
	}
*/	//-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
	
	SumOfPQ=NodeNum-1;
	SumOfPV=0;
	
	Pij=new double[BranchNum];  //支路有功功率
	Pji=new double[BranchNum];
	Qij=new double[BranchNum];  //支路无功功率
	Qji=new double[BranchNum];  //支路无功功率
	
	Iij=new Complex[BranchNum];  //支路电流
	Iji=new Complex[BranchNum];  //支路电流
	
	for(int i=0;i<NodeNum;i++)	//建立缺省的节点对照表
	{
		OptimisedNodeTable[i]=i;
		OriginalNodeTable[i]=i;
		Phase[i]=0;
		pNodeData[i].P=0;
		pNodeData[i].Q=0; 
	}	
	
	bSuccess=Y.Initialize(NodeNum,NodeNum);
	if(!bSuccess) return false;
	bSuccess=X.Initialize(NodeNum,NodeNum);
    if(!bSuccess) return false;
	bSuccess=B.Initialize(NodeNum,NodeNum);
    if(!bSuccess) return false;
	//	bSuccess=Adjacency.Initialize(NodeNum,NodeNum);
	Adjacency=new bool *[NodeNum];
	{
		for(int i=0;i<NodeNum;i++)
			Adjacency[i]=new bool[NodeNum];
	}
	{
		for(int i=0;i<NodeNum;i++)
			for(int j=0;j<NodeNum;j++)
				Adjacency[i][j]=false;
	}
		if(!bSuccess) return false;
		
		bPFconverge=false;
		bBackupedJacobian=false;
		bInitialized=true;
		
		//-----
		bParallelCompenTechSwitched=false;
		bUseParallelCompenTech=false;
		nAbnormalBranchNum=0;
		
		return true;
}

bool Flow::DCInitialize()
{
  dc_buses.assign(dc_bus_count,DCBus());
  dc_branches.assign(dc_branch_count,DCBranch());

  return true;
}

bool Flow::PreHandData()
{
	SumOfPV=0;
	SumOfLoadP=0.0;
	SumOfGenP=0.0;
	
	//统计PV节点，识别节点类型，找出平衡节点，计算节点净功率
	{
		for(int i=0;i<NodeNum;i++)
		{
			if(pNodeData[i].NodeType==2) SumOfPV++;
			else if(pNodeData[i].NodeType==1 || pNodeData[i].NodeType==0) pNodeData[i].NodeType=PQ; //统一设置为PQ节点
			else if(pNodeData[i].NodeType==3) SlackNodeOrder =pNodeData[i].Num;
			pNodeData[i].Num-=1;
			pNodeData[i].Gp/=BaseS;
			pNodeData[i].Gq/=BaseS;
			pNodeData[i].Lp/=BaseS;
			pNodeData[i].Lq/=BaseS;
			if(pNodeData[i].NodeType!=3) 
			{
				pNodeData[i].P=pNodeData[i].Gp-pNodeData[i].Lp;  //有功净功率
				SumOfGenP+=pNodeData[i].Gp;
			}
			pNodeData[i].Q=pNodeData[i].Gq-pNodeData[i].Lq;  //无功净功率
			//  if(pNodeData[i].NodeType!=PQ) 
			pNodeData[i].V=pNodeData[i].DesireVolt;
			
			SumOfLoadP+=pNodeData[i].Lp;
			
		}
	}

	{
		for ( int i = 0; i < BranchNum; ++i ) {
			pBranchData[i].B/=2;  //!!!!已经是B/2
			if(fabs(pBranchData[i].TRatio)>THRESHOLD) 
			{pBranchData[i].BranchType=TS;
			pBranchData[i].K=pBranchData[i].TRatio;
			}
		}
	}

	SumOfPQ=NodeNum-SumOfPV-1;
	
	return true;
}
//--------------------------------------------------------------------------------
void Flow::NodeSerialOptimize()	//节点编号优化
{
	SparseMatrix<int> A(NodeNum,NodeNum);	//声明关联矩阵
	int * outlet=new int[NodeNum];	        //节点出线表
	bool * NodeUsed=new bool[NodeNum];	    //节点处理完毕标志
	int * RelTable=new int[NodeNum];	    //节点相连节点表
	int i;
	for(i=0;i<NodeNum;i++)
	{
		NodeUsed[i]=false;
	}
	for(i=0;i<BranchNum;i++)	//设置关联矩阵
	{
		ACBranch * TheLine=pBranchData+i;
		if(TheLine->NodeI>=0&&TheLine->NodeJ>=0&&TheLine->BranchType!=CAP)  //去掉接地电容的影响
		{
			A.SetElement(TheLine->NodeI-1,TheLine->NodeJ-1,1);   
			A.SetElement(TheLine->NodeJ-1,TheLine->NodeI-1,1);
		}
	}
	for(int count=0;count<NodeNum;count++)
	{
		for(i=0;i<NodeNum;i++)	//统计出线数
		{
			if(NodeUsed[i]) continue;
			int sum=0;
			for(int j=0;j<NodeNum;j++)
			{
				if(NodeUsed[j]) continue;
				if(A.GetElement(i,j)==1)
					sum++;
			}
			outlet[i]=sum;
		}
		if(count!=NodeNum-1) outlet[SlackNodeOrder-1]=NodeNum;  //使平衡节点编号为lastone
		int min=MAXIMUM,no=0;
		for(i=0;i<NodeNum;i++)		//寻找出线最少节点
		{
			if(NodeUsed[i]) continue;
			if(outlet[i]<min)
			{
				min=outlet[i];
				no=i;
			}
		}
		
		NodeUsed[no]=true;	//找到出线最少节点
        OptimisedNodeTable[no]=count;
		OriginalNodeTable[count]=no;
		
		//计算消去节点影响
		int total=outlet[no];	//与节点相连的节点数目
		int mcount=0;
		for(int j=0;j<NodeNum;j++)
		{
			if(mcount==total) break;
			if(A.GetElement(no,j)==1)
			{
				RelTable[mcount]=j;
				mcount++;
			}
		}
		if(total==1) continue;
		else
			for(int j=0;j<total-1;j++)
				for(int k=j+1;k<total;k++)
					A.SetElement(RelTable[j],RelTable[k],1);	//添加生成的支路
	}
	
	//形成邻接矩阵
	int ni,nj;
    for(i=0;i<NodeNum;i++)