flow.cpp

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

	{
		// Adjacency.SetElement(i,i,1);
		Adjacency[i][i]=true;	
	}
	
	for(i=0;i<BranchNum;i++)	//设置关联矩阵
	{
		ACBranch * TheLine=pBranchData+i;
		if(TheLine->NodeI>=0&&TheLine->NodeJ>=0&&TheLine->BranchType!=CAP)  //去掉接地电容的影响
		{
			ni=OptimisedNode(TheLine->NodeI-1);
			nj=OptimisedNode(TheLine->NodeJ-1);
			// Adjacency.SetElement(ni,nj,1);
			// Adjacency.SetElement(nj,ni,1);
			Adjacency[ni][nj]=true;
			Adjacency[nj][ni]=true;
		}
	}
	
	
	delete [] outlet;
	delete [] NodeUsed;
	delete [] RelTable;
	
	/*
	fprintf(outfile,"节点优化编号如下：\n    原编号＝>优化编号  优化编号＝>原编号");
	for(i=0;i<NodeNum;i++)
	fprintf(outfile,"\n          %d => %d           %d => %d",i+1,OptimisedNodeTable[i]+1,i+1,OriginalNodeTable[i]+1);
	//fclose(outfile);
	*/
}
//--------------------------------------------------------------------------------
//------2004.10.29，针对原节点编号太耗时的缺点，写此改进版
void Flow::NodeSerialOptimize1()	//节点编号优化
{
	int (*A)[MAXOUTLET];	//声明关联矩阵
	int * outlet=new int[NodeNum];	//节点出线表
	bool * NodeUsed=new bool[NodeNum];	//节点处理完毕标志
	int * RelTable=new int[NodeNum];	//节点相连节点表
	int i,j,k,l;
	int ni,nj;
	
	
	A=new int[NodeNum][MAXOUTLET];
	for(i=0;i<NodeNum;i++)
	{
		NodeUsed[i]=false;
		outlet[i]=0;
		for(j=0;j<MAXOUTLET;j++)
			A[i][j]=0;
	}
	
	for(i=0;i<BranchNum;i++)	//设置关联矩阵
	{
		ACBranch * TheLine=pBranchData+i;
		ni=TheLine->NodeI-1;   
		nj=TheLine->NodeJ-1;
		if(ni>=0&&nj>=0&&TheLine->BranchType!=CAP)  //去掉接地电容的影响
		{
			A[ni][outlet[ni]]=nj;  	outlet[ni]++;
			A[nj][outlet[nj]]=ni;	outlet[nj]++;
		}
	}
	//使平衡节点编号为lastone
	OptimisedNodeTable[SlackNodeOrder-1]=NodeNum-1;
	OriginalNodeTable[NodeNum-1]=SlackNodeOrder-1;
	NodeUsed[SlackNodeOrder-1]=true;
	for(int count=0;count<NodeNum-1;count++)
	{
		//	if(count==SlackNodeOrder-1) continue;  
		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 mcount=0;
		int total=outlet[no];
		if(total==1) continue;		//只与一条支路相连，消去后不增加新支路
		for(j=0;j<total;j++)//与节点相连的节点数目
		{
			if(NodeUsed[A[no][j]]) continue;
			RelTable[mcount]=A[no][j];
			for(k=0;k<outlet[RelTable[mcount]];k++)
			{
				if(A[RelTable[mcount]][k]==no)  //除去已消去的支路
				{
					for(l=k;l<outlet[RelTable[mcount]]-1;l++)
						A[RelTable[mcount]][l]=A[RelTable[mcount]][l+1];
					A[RelTable[mcount]][l]=0;
					outlet[RelTable[mcount]]--;
					break;
				}
			}
			mcount++;
		}
		//添加生成的支路
		for(j=0;j<mcount-1;j++)
			for(k=j+1;k<mcount;k++)
			{
				//对于支路j->k
				for(l=0;l<outlet[RelTable[j]];l++)
					if(A[RelTable[j]][l]==RelTable[k]) break;
					if(l>=outlet[RelTable[j]])		//没找到
					{	//加到最后							
						A[RelTable[j]][outlet[RelTable[j]]]=RelTable[k];
						outlet[RelTable[j]]++;
					}
					
					//对于支路K->J			
					for(l=0;l<outlet[RelTable[k]];l++)
						if(A[RelTable[k]][l]==RelTable[j]) break;
						if(l>=outlet[RelTable[k]])		//没找到
						{	//加到最后
							A[RelTable[k]][outlet[RelTable[k]]]=RelTable[j];
							outlet[RelTable[k]]++;
						}
						
			}
			
	}
	
	//形成邻接矩阵
	// int ni,nj;
    for(i=0;i<NodeNum;i++)
	{
		// Adjacency.SetElement(i,i,1);
		Adjacency[i][i]=true;	
	}
	
	for(i=0;i<BranchNum;i++)	//设置关联矩阵
	{
		ACBranch * TheLine=pBranchData+i;
		if(TheLine->NodeI>=0&&TheLine->NodeJ>=0&&TheLine->BranchType!=CAP)  //去掉接地电容的影响
		{
			ni=OptimisedNode(TheLine->NodeI-1);
			nj=OptimisedNode(TheLine->NodeJ-1);
			//			Adjacency.SetElement(ni,nj,1);
			//			Adjacency.SetElement(nj,ni,1);
			Adjacency[ni][nj]=true;
			Adjacency[nj][ni]=true;
		}
	}	
		
	delete [] outlet;
	delete [] NodeUsed;
	delete [] RelTable;

	delete [] A;
	
	outfile<<std::endl<<"节点优化编号如下："<<std::endl<<"    原编号＝>优化编号  优化编号＝>原编号"<<std::endl;
	for(i=0;i<NodeNum;i++)
		outfile<<"         "<<i+1<<" => " <<OptimisedNodeTable[i]+1<<"          "<<i+1<<" => " <<OriginalNodeTable[i]+1<<std::endl;	
}

void Flow::NodeSerialOptimize2()
{
  typedef boost::adjacency_list<boost::setS,boost::vecS,boost::undirectedS> Graph;

  typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;

  typedef boost::graph_traits<Graph>::vertex_iterator VertexIter;
  typedef boost::graph_traits<Graph>::adjacency_iterator AdjVertexIter;

  typedef boost::graph_traits<Graph>::edge_iterator EdgeIter;
  typedef boost::graph_traits<Graph>::edges_size_type EdgeSizeType;

  typedef std::ptrdiff_t DiffType;

  typedef std::vector<bool> BooleanCont;
  typedef BooleanCont::iterator BooleanIter;

  // Step 1. Build a graph that models the network.
  Graph graph(NodeNum);
  {
    for ( int i = 0; i < BranchNum; ++i ) {
      ACBranch const& branch = pBranchData[i];

      int from = branch.NodeI - 1;
      int to = branch.NodeJ - 1;

      // Skip grounded branches.
      if ( from >= 0 && to >= 0 && branch.BranchType != CAP ) {
        boost::add_edge(from,to,graph);
      }
    }
  }

  // Step 2. Prepare auxiliary data.
  BooleanCont labeled(NodeNum,false);

  // Step 3. Label the slack bus and then eliminate it from the graph.
  {
    Vertex slack = SlackNodeOrder - 1;

    OptimisedNodeTable[slack] = NodeNum - 1;
    OriginalNodeTable[NodeNum-1] = 
#if defined(_MSC_VER) && _MSC_VER <=1200
      slack;
#else
      numeric_cast<int>(slack);
#endif
    labeled[slack] = true;

    EliminateVertex(slack,graph);
  }

  // Step 4. Execute dynamic ordering on the remaining vertices.
  {
    VertexIter ver_beg, ver_end;
    boost::tie(ver_beg,ver_end) = boost::vertices(graph);

    BooleanIter const labeled_beg = labeled.begin();
    BooleanIter const labeled_end = labeled.end();

    for( int id = 0; id < NodeNum - 1; ++id ) {
      // Locate the first unlabeled vertex.
      BooleanIter labeled_iter = std::find(labeled_beg,labeled_end,false);

      // Get the first unlabeled vertex and assume that it has the fewest injected edges.
      DiffType ver_tar_off_set = std::distance(labeled_beg,labeled_iter);

      VertexIter ver_iter = ver_beg;
      std::advance(ver_iter,
#if defined(_MSC_VER) && _MSC_VER <=1200
        ver_tar_off_set
#else
        numeric_cast<std::iterator_traits<VertexIter>::difference_type>(ver_tar_off_set)
#endif
      );
      Vertex ver_tar = *ver_iter;

      EdgeSizeType ver_tar_inj_edges = CountInjectedEdges(ver_tar,graph);

      // Search the unlabeled vertex with fewest injected edges after its elimination.
      for ( ver_iter = ver_beg; ver_iter != ver_end; ++ver_iter ) {
        DiffType ver_cur_off_set = 
#if defined(_MSC_VER) && _MSC_VER <= 1200
          std::distance(ver_beg,ver_iter);  
#else
          numeric_cast<DiffType>(std::distance(ver_beg,ver_iter));
#endif
        // Skip the labeled vertex.
        if ( labeled[ver_cur_off_set] ) continue;

        // Skip the same vertex.
        if ( *ver_iter == ver_tar ) continue;

        // Check the count of current vertex's inject edges.
        Vertex ver_cur = *ver_iter;
        EdgeSizeType ver_cur_inj_edges = CountInjectedEdges(*ver_iter,graph);

        if ( ver_cur_inj_edges < ver_tar_inj_edges ) {
          ver_tar_off_set = ver_cur_off_set;
          ver_tar = ver_cur;
          ver_tar_inj_edges = ver_cur_inj_edges;
        }                
      }

      // Order the target vertex.
      OptimisedNodeTable[ver_tar] = id;
      OriginalNodeTable[id] = 
#if defined(_MSC_VER) && _MSC_VER <=1200
        ver_tar;
#else
        numeric_cast<int>(ver_tar);
#endif

      labeled[ver_tar_off_set] = true;

      // Eliminate the just labeled vertex and add injected edges into the graph.
      EliminateVertex(ver_tar, graph);
    }  
  }

  { // Construct the adjacent matrix.
    { // Set diagonal elements.
      for ( int i = 0; i < NodeNum; ++i ) {
        // Adjacency.SetElement(i,i,1);
        Adjacency[i][i]=true;	
      }
    }

    { // Set off-diagonal elements.
      for( int i = 0; i < BranchNum; ++i ) {
        ACBranch const& branch = pBranchData[i];
        if( branch.NodeI >= 0 && branch.NodeJ >= 0 && branch.BranchType != CAP ) {
          int ni = OptimisedNode(branch.NodeI-1);
          int nj = OptimisedNode(branch.NodeJ-1);
          // Adjacency.SetElement(ni,nj,1);
          // Adjacency.SetElement(nj,ni,1);
          Adjacency[ni][nj]=true;
          Adjacency[nj][ni]=true;
        }
      }	
    }
  }

  { // Output node order maps.
    boost::io::ios_flags_saver flag_saver(outfile);

    outfile << '\n' 
            << "节点优化编号如下：" << '\n' 
            << std::setw(24) << "原始编号 -> 优化编号" 
            << std::setw(24) << "优化编号 -> 原始编号"
            << '\n';
    for( int i = 0; i < NodeNum; ++i )
      outfile << std::setw(12) << std::right << i + 1 << " -> " 
              << std::setw(8) << std::left << OptimisedNodeTable[i] + 1 
              << std::setw(12) << std::right << i + 1 << " -> " 
              << std::setw(8) << std::left << OriginalNodeTable[i] + 1 
              << '\n';
    
    flag_saver.restore();
  }

  return;
}

//-------------在形成导纳矩阵之前就要设置-------------------------------------
void Flow::SwitchPallelCompenTech(bool bUseIt)
{
	if(!bParallelCompenTechSwitched)
	{
		bUseParallelCompenTech=bUseIt;
		bParallelCompenTechSwitched=true;
	}
}
//------统计是否有小阻抗和大R/X比的支路
//------需要对支路扫描三次，第一次统计R和X的均值
//------第二次找小阻抗和大R/X比的支路，如果有则作标记
//------第三次正式记录这些支路
//------并计算替代的正常支路和并补支路的G，B
//
void Flow::StatisAbnormalBranches()
{
	int i;
	ACBranch * pb;
	double AvergR,AvergX,AvergZ2;
	double R,X,Z2;
	double SumR,SumX;
	double RXratio;
	double Threshold=5;  //门槛,临界比例，这是认为选择的
	
	nAbnormalBranchNum=0;
	if(!bUseParallelCompenTech) return;
	SumR=0.0;SumX=0.0;
	for(i=0;i<BranchNum;i++)	
	{
		pb=pBranchData+i;
		if(pb->BranchType!=LN && pb->BranchType!=TS) continue;
		SumR+=fabs(pb->R);
		SumX+=fabs(pb->X);
	}
	AvergR=SumR/BranchNum;
	AvergX=SumX/BranchNum;
	RXratio=AvergR/AvergX;
	AvergZ2=AvergR*AvergR+AvergX*AvergX;
	
	for(i=0;i<BranchNum;i++)	
	{
		pb=pBranchData+i;
		if(pb->BranchType!=LN && pb->BranchType!=TS) continue;
		R=pb->R;
		X=pb->X;
		if(fabs(R/X)>Threshold)	//大比值
		{	pb->bExcludeRX=true;  nAbnormalBranchNum++;}
		else	//小阻抗
		{
			Z2=(R*R+X*X)*Threshold*Threshold*100;
			if(Z2<AvergZ2) 
			{pb->bExcludeRX=true;  nAbnormalBranchNum++;}
		}
		
	}
	//
	int Counter=0;
	Complex Yij;
	double Gij,Bij;
	ParalCompensBranches= new ParallelCompensBranch[nAbnormalBranchNum];
	for(i=0;i<BranchNum;i++)	
	{
		pb=pBranchData+i;
		if(pb->bExcludeRX)
		{
			R=pb->R; X=pb->X;
			Yij=Complex(1.0,0)/Complex(R,X);
			Gij=Yij.Real; Bij=Yij.Imag;
			ParalCompensBranches[Counter].nBranchOrder=i;
			ParalCompensBranches[Counter].NodeI=pb->NodeI;
			ParalCompensBranches[Counter].NodeJ=pb->NodeJ;
			ParalCompensBranches[Counter].BranchType=pb->BranchType;
			if(pb->BranchType==LN)
			{ParalCompensBranches[Counter].B_2=pb->B_2;
			}
			else
			{ParalCompensBranches[Counter].K=pb->K;
			}
			ParalCompensBranches[Counter].Gn=1/AvergX;  
			ParalCompensBranches[Counter].Bn=-1/AvergX;
			ParalCompensBranches[Counter].Gc=Gij-1/AvergX; //注意这里的推导，令Gc=0
			ParalCompensBranches[Counter].Bc=Bij+1/AvergX;
			Counter++;
		}
	}
	
}

//-----建立导纳矩阵--------------------------
void Flow::BuildMatrixY()		//形成导纳矩阵
{
	int i;
	ACBranch * tb;
	int ni,nj;		//优化后的节点编号
	Complex tii,tij,tji,tjj,tc;

  Y.DeleteAllElements();  // 2008.05.21 Clear all non-zero elements.

	for(i=0;i<BranchNum;i++)
	{
		tb=pBranchData+i;
		if(tb->bExcludeRX) 
			continue;  //!!!!!附加的部分
		
		ni=OptimisedNode(tb->NodeI-1);
		if(tb->BranchType!=CAP) nj=OptimisedNode(tb->NodeJ-1);
		
		tii=GetY(ni,ni);	//获得导纳矩阵原来元素的值
		if(tb->BranchType!=CAP)
		{
			tij=GetY(ni,nj);
			tji=GetY(nj,ni);
			tjj=GetY(nj,nj);
		}
		
		switch(tb->BranchType)
		{
		case LN:
			tc=Complex(1.0,0)/Complex(tb->R,tb->X);
			tii=tii+tc;
			tii.Imag+=tb->B_2;
			tjj=tjj+tc;
			tjj.Imag+=tb->B_2;
			tij=tij-tc;
			tji=tji-tc;
			break;
		case TS:
			{   //I为非标准侧
				Complex y10,y12,y20;
				tc=Complex(1.0,0)/Complex(tb->R,tb->X);
				y12=tc/tb->K;
				y20=y12/(tb->K);
				y10=tc;
				tii=tii+y20;
				tjj=tjj+y10;
				tij=tij-y12;
				tji=tji-y12;
			}
			break;
		case CAP:	//接地电容
			{
				tii=tii+Complex(0,tb->B);
			}
			break;
		case TSP:	//移相器
			{
				Complex k;	//复数变比
				Complex y10,y12,y20;
				k=Complex(tb->K*cos(tb->TAngle),tb->K*sin(tb->TAngle));
				tc=Complex(tb->R,tb->X);	//Z
				y12=Complex(1.0,0)/tc/k;
				y10=(Complex(1,0)-k.Conj())/k.Conj()/k/tc;
				y20=(k-Complex(1,0))/k/tc;
				tii=tii+y20;
				tjj=tjj+y10;
				tij=tij-y12;
				tji=tji-y12;
			}
			break;
		default:
			assert(false);
		}
		SetY(ni,ni,tii);	//设置导纳矩阵的值
		if(tb->BranchType!=CAP){
			SetY(ni,nj,tij);
			SetY(nj,ni,tji);
			SetY(nj,nj,tjj);
		}
	}
	
	//建立导纳矩阵还需要加入各节点导纳
	for(i=0;i<NodeNum;i++)
	{
		if(fabs(pNodeData[i].ShuntG)>THRESHOLD || fabs(pNodeData[i].ShuntB)>THRESHOLD)
		{
			int n=pNodeData[i].Num;
			ni=OptimisedNode(n);
			tii=GetY(ni,ni);
			tii=tii+Complex(pNodeData[i].ShuntG,pNodeData[i].ShuntB);
			SetY(ni,ni,tii);