3117040428.cpp

课程设计时做得潮流程序

// 潮流计算.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include<fstream>
#include<stdio.h>
#include<iostream>
#include<stdlib.h>
#include <complex>
#include<math.h>
#include "NEquation.h"
#define nodenumber 9   //总接点数
#define PVshu 2        //PV接点数
#define PHshu 1        //平衡接点数
#define lineshu 9      //线路数
#define Eaf 0.000000000001    //结束循环值
using namespace std;
void main(int argc, char* argv[])
{	
   	   
//************************************************************************************************************
   ifstream infile("nodedate.txt",ios::in|ios::binary);     //导入节点参数
     if(!infile)
	 { cout<<"can not open"<<endl;
	   exit(1);
	 }
     int i=0,j=0;
	int nodebianhao[nodenumber+1];
	int nodetype[nodenumber+1];
	float Q[nodenumber+1];
	float V[nodenumber+1];
	float P[nodenumber+1];   
	float D[nodenumber+1];            //节点电压角度
    for(i=1;i<=nodenumber;i++)
	   infile>>nodebianhao[i]  >>nodetype[i]  >>P[i]    >>Q[i]   >>V[i]  >>D[i];
    cout <<"节点参数" <<endl;
	for(i=1;i<=nodenumber;i++)
		cout<<nodebianhao[i] <<'\t'<<nodetype[i] <<'\t'<<P[i] <<'\t'<<Q[i] <<'\t'<<V[i]<<'\t'<<D[i]<<endl;	
          cout<<"******************************************************"<<endl;
		  infile.close();
 
  
cout<<"平衡结点数"<<PHshu<<'\t'<<"PV结点数"<<PVshu<<'\t'<<endl;

 //**************************************************************************************************************
     ifstream infile2("linedate.txt",ios::in|ios::binary);//导入线路参数 
	 if(!infile2)
	 { cout<<"can not open"<<endl;
	   exit(1);
	 }
     int k;
	 float R=0,X=0,B=0;
	 i=0;
	 j=0;
	 complex<float>y[nodenumber+1][nodenumber+1];//导纳矩阵
	 
	 for(i=0;i<=nodenumber;i++)    //初始化
		 for(j=0;j<=nodenumber;j++)
			 y[i][j]=0;
	
	 for(k=0;k<lineshu;k++)
	    {
		 infile2>>i >>j >>R >>X  >>B;
	      complex<float> fu1(R,X);
		  complex<float> fu2(0,B);
		  complex<float> fu(1,0);
          y[i][j]=-fu/fu1;
	      y[j][i]=y[i][j];
	      y[i][i]=y[i][i]+fu/fu1+fu2;
	      y[j][j]=y[j][j]+fu/fu1+fu2;
	
        }
		infile2.close();
	cout<<"**********导纳矩阵y[][]**********"<<endl;
	for(i=1;i<=nodenumber;i++)
	{ 
		for(j=1;j<=nodenumber;j++)
		cout<<y[i][j];
	   	cout<<endl;
	}
 ofstream onfile("导纳矩阵.txt",ios::out);
    for(i=1;i<=nodenumber;i++)
	{ 
	   for(j=1;j<=nodenumber;j++)
         onfile<<y[i][j]<<'\t' ;
            onfile<<endl;
	}
   
   
   onfile.close();
 cout<<"***********************************"<<endl;
		
//*******************************************************************求导纳模和角度
	float Y[nodenumber+1][nodenumber+1],Yd[nodenumber+1][nodenumber+1];//&d-导纳角
       for(i=1;i<=nodenumber;i++)   
		 for(j=1;j<=nodenumber;j++)
		 { 
			 Y[i][j]=abs(y[i][j]);
		     Yd[i][j]=arg(y[i][j]);
		 }
cout<<"************导纳矩阵幅值***********************"<<endl;	  
    for(i=1;i<=nodenumber;i++)
	{	 
		for(j=1;j<=nodenumber;j++)
		  cout<<Y[i][j] <<'\t';
		 cout<<endl;
	}
ofstream onfile1("导纳模.txt",ios::out);
   for(i=1;i<=nodenumber;i++)
	{	
	     for(j=1;j<=nodenumber;j++)
		  onfile1<<Y[i][j] <<'\t';
		 onfile1<<endl;
	}
onfile1.close();
	
cout<<"***********导纳矩阵相角************************"<<endl;
	for(i=1;i<=nodenumber;i++)
	{	 
		for(j=1;j<=nodenumber;j++)
		   cout<<Yd[i][j] <<'\t';
		cout<<endl;
	}
ofstream onfile2("导纳相角.txt",ios::out);
   	for(i=1;i<=nodenumber;i++)
	{
		for(j=1;j<=nodenumber;j++)
		  onfile2<<Yd[i][j] <<'\t';
	     onfile2<<endl;
	}
	onfile2.close();

	 cout<<"***********************************"<<endl;
//*****************************************************************************计算不平衡量
    float PP[nodenumber+1]={0},QQ[nodenumber+1]={0},DP[nodenumber+1]={0},DQ[nodenumber+1]={0};
	float DD[nodenumber+1]={0};
	float DV[nodenumber+1]={0};
	 int a;
 	 int n; 
	 int Z;
	 ofstream ofx("J1234.txt") ;    //打印分块雅
	 ofstream onfile3("jacobian.txt",ios::out);    //输出jacobian
	 ofstream ob3("修正量.txt",ios::out);//打印每次不平衡量
	 ofstream ob5("电压相角.txt",ios::out);//打印每次电压相角.
     ofstream ob6("不平衡量.txt",ios::out);//
	  ofstream ob8("功率.txt",ios::out);//
for(Z=1;Z<20;Z++)   //总循环  循环次数最多为20
 {	
	for(i=1;i<=nodenumber+1;i++)
	{
		PP[i]=0;
	    QQ[i]=0;
	    DP[i]=0;
	    DQ[i]=0;
        DD[i]=0;
        DV[i]=0;
	}
	for(k=1;k<=nodenumber;k++)
	{ 
		for(n=1;n<=nodenumber;n++)
		{
		    PP[k]=PP[k]+V[k]*Y[k][n]*V[n]*cos(D[k]-D[n]-Yd[k][n]);
	        QQ[k]=QQ[k]+V[k]*Y[k][n]*V[n]*sin(D[k]-D[n]-Yd[k][n]);
		}
	}
	 ob8<<"************第"<<Z<<"次迭代**********";
	  ob8<<"\n**PP<<'\t'<<QQ**"<<endl;
    for(k=1;k<=nodenumber;k++)
	{ 
			ob8<<PP[k]<<'\t'<<P[k]<<'\t'<<QQ[k]<<'\t'<<Q[k]<<endl;
	
	}
   
 //*********************************************************************定义不平横量
	 for(k=1;k<=nodenumber;k++)//从DP[1]到DP[9]
	    DP[k]=P[k]-PP[k];
	
	 for(k=1;k<=nodenumber;k++)//从DQ[1]到DQ[9]
	    DQ[k]=Q[k]-QQ[k];
   ob6<<"************第"<<Z<<"次迭代**********";
   ob6<<"\n**DP<<'\t'<<DQ**"<<endl;
        for(k=1;k<=nodenumber;k++)
			ob6<<DP[k]<<'\t'<<DQ[k]<<endl;
//*************************************************************************定义修正量
    
	     
//**********************************************************************求jacobian
     float J1[nodenumber-PHshu][nodenumber-PHshu]={0};
	 float J2[nodenumber-PHshu][nodenumber-PVshu-PHshu]={0};
	 float J3[nodenumber-PVshu-PHshu][nodenumber-PHshu]={0};
	 float J4[nodenumber-PVshu-PHshu][nodenumber-PVshu-PHshu]={0}; 
	 float J[nodenumber-PHshu+nodenumber-PHshu-PVshu][nodenumber-PHshu+nodenumber-PHshu-PVshu]={0};
	for(i=PHshu+1;i<=nodenumber;i++)      //求J1  8X8
	  {   
		for(j=PHshu+1;j<=nodenumber;j++)
		{
			if(i!=j)	  
				J1[i-PHshu-1][j-PHshu-1]=V[i]*Y[i][j]*V[j]*sin(D[i]-D[j]-Yd[i][j]);
			else 
			{	
				for(a=1;a<=nodenumber;a++)  //
				   J1[i-PHshu-1][j-PHshu-1]-=V[i]*Y[i][a]*V[a]*sin(D[i]-D[a]-Yd[i][a]);
			    J1[i-PHshu-1][j-PHshu-1]+=V[i]*Y[i][j]*V[j]*sin(D[i]-D[j]-Yd[i][j]);
			}
		}
	  }
 ofx<<"\n\n*************J1*********************"<<endl;
    for(i=0;i<nodenumber-PHshu;i++)      //J1 
	{  
		for(j=0;j<nodenumber-PHshu;j++)
		   ofx << J1[i][j] <<"\t" ;
	       ofx << endl ;
	}

	  for(i=PHshu+1;i<=nodenumber;i++)      //求J2   8X6
	  {   
		  for(j=PHshu+PVshu+1;j<=nodenumber;j++)
		  {	if(i!=j)	  
				J2[i-PHshu-1][j-PHshu-PVshu-1]=V[i]*Y[i][j]*cos(D[i]-D[j]-Yd[i][j]);
			else 
			{	
				for(a=1;a<=nodenumber;a++)
				   J2[i-PHshu-1][j-PHshu-PVshu-1]+=Y[i][a]*V[a]*cos(D[i]-D[a]-Yd[i][a]);
			    J2[i-PHshu-1][j-PHshu-PVshu-1]+=V[i]*Y[i][j]*cos(Yd[i][j]);
			}
		  }
	  }

ofx<<"\n***********J2************************"<<endl;
     for(i=0;i<nodenumber-PHshu;i++)      //J2 
	 { 
	   for(j=0;j<nodenumber-PHshu-PVshu;j++)
		 ofx <<J2[i][j]<<'\t';
       ofx << endl ;
	 }

    for(i=PHshu+PVshu+1;i<=nodenumber;i++)      //求J3 6X8
	{   
		for(j=PHshu+1;j<=nodenumber;j++)
	
		{	if(i!=j) 
				 J3[i-PHshu-PVshu-1][j-PHshu-1]-=V[i]*Y[i][j]*V[j]*cos(D[i]-D[j]-Yd[i][j]);
			 else 
			 { 
				 for(a=1;a<=nodenumber;a++)
				     J3[i-PHshu-PVshu-1][j-PHshu-1]+=V[i]*Y[i][a]*V[a]*cos(D[i]-D[a]-Yd[i][a]);
				  J3[i-PHshu-PVshu-1][j-PHshu-1]-=V[i]*Y[i][j]*V[i]*cos(D[i]-D[j]-Yd[i][j]);
		
			 }
		}
	}

	  ofx <<"\n\n*************J3**********************"<<endl;
  for(i=0;i<nodenumber-PHshu-PVshu;i++)      //J3   
	{  
	    for(j=0;j<nodenumber-PHshu;j++)
           ofx<<J3[i][j]<<'\t';
		   ofx<<endl;
	}

  for(i=PHshu+PVshu+1;i<=nodenumber;i++)      //求J4   6X6
  {   
	     for(j=PHshu+PVshu+1;j<=nodenumber;j++)
		 {	if(i!=j)	  
				J4[i-PHshu-PVshu-1][j-PHshu-PVshu-1]=V[i]*Y[i][j]*sin(D[i]-D[j]-Yd[i][j]);
			else 
			{	for(a=1;a<=nodenumber;a++)
				  J4[i-PHshu-PVshu-1][j-PHshu-PVshu-1]+=Y[i][a]*V[a]*sin(D[i]-D[a]-Yd[i][a]);
				J4[i-PHshu-PVshu-1][j-PHshu-PVshu-1]-=V[i]*Y[i][j]*sin(Yd[i][j]);			
			}
		 }
  }
ofx<<"\n\n***************J4********************"<<endl;
    for(i=0;i<nodenumber-PHshu-PVshu;i++)      //J4
	{  
		for(j=0;j<nodenumber-PHshu-PVshu;j++)
           ofx<<J4[i][j]<<'\t';
		  ofx<<endl;
	}

//**************************************************************合并
	
	    for(i=0;i<nodenumber-PHshu+nodenumber-PHshu-PVshu;i++)
		 for(j=0;j<nodenumber-PHshu+nodenumber-PHshu-PVshu;j++)
		 {   if((i<nodenumber-PHshu)&&(j<nodenumber-PHshu))
				     J[i][j]=J1[i][j];
		     if((i<nodenumber-PHshu)&&(j>=nodenumber-PHshu))
		             J[i][j]=J2[i][j-nodenumber+PHshu]; 
             if((i>=nodenumber-PHshu)&&(j<nodenumber-PHshu))                                       
                     J[i][j]=J3[i-nodenumber+PHshu][j];
			 if((i>=nodenumber-PHshu)&&(j>=nodenumber-PHshu))
			         J[i][j]=J4[i-nodenumber+PHshu][j-nodenumber+PHshu];		
		 }
cout<<"*********************jacobian******************************"<<endl;
   cout <<"********第"<<Z<<"次迭代******\n";
	for(i=0;i<nodenumber-PHshu+nodenumber-PHshu-PVshu;i++)
	{  
			for(j=0;j<nodenumber-PHshu+nodenumber-PHshu-PVshu;j++)
		           cout<<J[i][j]<<'\t';
			cout<<endl;
	}

   onfile3<<"********第"<<Z<<"次迭代******"<<endl;
  for(i=0;i<nodenumber-PHshu+nodenumber-PHshu-PVshu;i++)
	{   
	        for(j=0;j<nodenumber-PHshu+nodenumber-PHshu-PVshu;j++)
		           
		           onfile3<<J[i][j]<<'\t';
		           onfile3<<endl;
	}


//***************************************************************************解方程
	 
NEquation ob;
ob.SetSize(nodenumber-PHshu+nodenumber-PHshu-PVshu);
	for(i=0;i<nodenumber-PHshu+nodenumber-PHshu-PVshu; i++)	
	{
			if(i<nodenumber-PHshu)
			  ob.Value(i)= DP[i+PHshu+1];
		    if(i>=nodenumber-PHshu)
              ob.Value(i)= DQ[i-nodenumber+PHshu+PHshu+PVshu+1];
		for(j=0;j<nodenumber-PHshu+nodenumber-PHshu-PVshu; j++)
			ob.Data(i,j) = J[i][j];   
	}
ob.Run();
	for(i=0;i<nodenumber-PHshu+nodenumber-PHshu-PVshu; i++)
	{   
		if(i<nodenumber-PHshu)
	        DD[i+PHshu+1]=ob.Value(i);
	    if(i>=nodenumber-PHshu)
            DV[i-nodenumber+PHshu+PHshu+PVshu+1]=ob.Value(i);
	}



        ob3<<"********第"<<Z<<"次迭代******\n";     //打印
		ob3<<"*****DV '\t' DD*****\n";
	for(i=1;i<=nodenumber;i++)
		{      
               ob3<<DV[i]<<'\t'<<DD[i]<<endl;
		}
//***************************************************************************修正量
    for(i=1;i<=nodenumber;i++)      
	{	V[i]+=DV[i];
        D[i]+=DD[i];
	}
//****************************************************************************打印
    ob5<<"********第"<<Z<<"次迭代******\n";
	ob5<<"*****电压 '\t' 相角*****\n";
	for(i=1;i<=nodenumber;i++)  
	{
		
               ob5<<V[i]<<'\t'<<D[i]<<endl;
	}

	  cout<<"********第"<<Z<<"次迭代**电压差值和相角的差值****\n";
	for(i=1;i<=nodenumber;i++)                                      //输出每次修正量 
	 {
		
               cout<<DV[i]<<'\t'<<DD[i]<<endl;
	 }
     cout<<"********第"<<Z<<"次迭代**电压和相角****\n";
	for(i=1;i<=nodenumber;i++)                                   //输出每次电压相角
	 {
		 
               cout<<V[i]<<'\t'<<D[i]<<endl;

	 }
	cout<<"********第"<<Z<<"次迭代****功率PQ差值**\n";
	for(i=1;i<=nodenumber;i++)                                     //输出每次功率PQ差值
	 {
		 
               cout<<DP[i]<<'\t'<<DQ[i]<<endl;

	 }
	cout<<"********第"<<Z<<"次迭代****功率PQ**\n";
	for(i=1;i<=nodenumber;i++)                                       //输出每次功率
	 {
		 
               cout<<PP[i]<<'\t'<<QQ[i]<<endl;

	 }
//*******************************************************************判断是否退出循环
   int r=0;
   for(i=1;i<=nodenumber;i++)  
   { 
	   if(fabs(DV[i])<Eaf)
		   r++;
       if(fabs(DD[i])<Eaf)
		  r++;
   }
   if(r==2*nodenumber)  //判断是否退出循环
      break;
}
 ofx.close();
 ob3.close();
 ob5.close();
 onfile3.close();
 ob6.close();
 ob8.close();
 //************************循环结束计算功率
 
    for(n=1;n<=nodenumber;n++)
	{  P[1]+=V[1]*Y[1][n]*V[n]*cos(D[1]-D[n]-Yd[1][n]);	
	   Q[1]+=V[1]*Y[1][n]*V[n]*sin(D[1]-D[n]-Yd[1][n]);	
	}
	for(k=PHshu+1;k<=PVshu+PHshu+1;k++)
	{ 
		for(n=1;n<=nodenumber;n++)
	       Q[k]=Q[k]+V[k]*Y[k][n]*V[n]*sin(D[k]-D[n]-Yd[k][n]);	
	}
  
 ofstream ob7("节点功率.txt",ios::out);
     for(k=1;k<=nodenumber;k++)
	 { 
	   ob7<< P[k]<<'\t'<<Q[k]<<endl;  
	        	
	 }
 
 ob7.close();
 //*******************************************************************计算线路损耗
ifstream inxl("linedate.txt",ios::in|ios::binary);//导入节点参数
     if(!infile)
	 { cout<<"can not open"<<endl;
	   exit(1);
	 }
 
     complex<float>yy[nodenumber+1][nodenumber+1]={0};//线路
	 complex<float>ss[nodenumber+1][nodenumber+1]={0};//Sij损耗
	 complex<float>sDs[nodenumber+1][nodenumber+1]={0};//DSij
     complex<float>VFC[nodenumber+1]={0};
	 complex<float>VF[nodenumber+1]={0};
	 complex<float>szong;
	  for(i=1;i<=nodenumber;i++)
	{ 
        VFC[i]=complex<float>(V[i]*cos(D[i]),-V[i]*sin(D[i]));    //电压的共轭
		VF[i]=complex<float>(V[i]*cos(D[i]),V[i]*sin(D[i]));      //复数形式的电压
	}   
	 for(i=0;i<=nodenumber;i++)    //初始化
		 for(j=0;j<=nodenumber;j++)
			 yy[i][j]=0;
	 for(k=0;k<lineshu;k++)
	 {
		 inxl >>i >>j >>R >>X  >>B;
	      
		   yy[i][j]=complex<float>(R/(R*R+X*X),X/(R*R+X*X));
	       yy[j][i]=y[i][j];
	       yy[i][0]=complex<float>(0,B);
	       yy[j][0]=complex<float>(0,B);
		   ss[i][j]=VF[i]*(VFC[i]*yy[i][0]+(VFC[i]-VFC[j])*yy[i][j]);  //接点输出功率
		   ss[j][i]=VF[j]*(VFC[j]*yy[j][0]+(VFC[j]-VFC[i])*yy[j][i]);  //接点输出功率
	 }
  inxl.close();

   for(i=1;i<=nodenumber;i++)
	{	 
		for(j=1;j<=nodenumber;j++)
		  if(i<j)	
		  {
			  sDs[i][j]=ss[i][j]+ss[j][i];    //线路损耗
			  sDs[j][i]=sDs[i][j];
		  }
   }
  ofstream ob9("线路损耗.txt",ios::out);//*************打印线路损耗
     ob9<<"\n*******线路损耗*******\n";
	  for(i=1;i<=nodenumber;i++)
	  {	 
		for(j=1;j<=nodenumber;j++)
		  if(i<j)	
           szong+=sDs[i][j];
	  }
     for(i=1;i<=nodenumber;i++)
	{	 
		for(j=1;j<=nodenumber;j++)
		  if(i<j)
			  ob9<<"***第"<<i<<"到第"<<j<<"条支路"<<ss[i][j]<<'\t'<<"***第"<<j<<"到第"<<i<<"条支路"<<ss[j][i]<<endl;
		     
			  ob9<<"***第"<<i<<j<<"条支路"<<sDs[i][j]<<endl;

	 }
	  ob9<<"***总****"<<szong<<endl;
	  ob9.close();


}