nlf.cpp

计算电力潮流

#include "nlf.h"
#include <fstream>
#include <iostream>
#include <cstdlib>

nlf::nlf()
{
	if((in=fopen("in.txt","r"))==NULL)
    	puts("Data file open error!");
	if((out=fopen("out.txt","w"))==NULL)
		puts("Data file open error!");
}
nlf::~nlf()
{
	fclose(in);
	fclose(out);
}
void nlf::readdata()
{
	int i,q1,q2,f;
	PointData f1,f2,f3,f4,f5;
	//读取串连支路数据
	for(i=1;i<=number_branch;i++)
	{
	    fscanf(in,"%d %d %lf %lf %lf",&q1,&q2,&f1,&f2,&f3);
	    branch[i].i_node=q1;
		branch[i].j_node=q2;
		branch[i].r=f1;
		branch[i].x=f2;
		/*if(f3<0)
			branch[i].bk=f3;
		else
            branch[i].bk=f3/2;*/
		branch[i].bk=f3;
	}
	//读取接地支路数据
	for(i=1;i<=number_ground;i++)
	{
		fscanf(in,"%d %lf",&q1,&f1);
		ground[i].node=q1;
		ground[i].x=f1;
	}
   /*//读发电机数据
	for(i=1;i<=number_genrt;i++)
	{
		fscanf(in,"%d %lf %lf %lf ",&q1,&f1,&f2,&f3);
		genrt[i].node=q1;
		genrt[i].p=f1;
		genrt[i].q=f2;
		genrt[i].v=f3;
	}
	//读负荷数据
	for(i=1;i<=number_load;i++)
	{  
		fscanf(in,"%d %lf %lf",&q1,&f1,&f2);
		load[i].node=q1;
		load[i].p=f1;
		load[i].q=f2;
	}*/
    //读取所有节点信息
	int j=1,k=1,l=1;
	for(i=1;i<=N;i++)
	{
		fscanf(in,"%d %lf %lf %lf %lf %lf %d",&q1,&f1,&f2,&f3,&f4,&f5,&f);
		node[i].node=q1;
		node[i].Pg=f1;
		node[i].Qg=f2;
		node[i].Pl=f3;
		node[i].Ql=f4;
		node[i].v=f5;
		node[i].flag=f;

		if(node[i].flag==1)//提取负荷信息
		{
			load[k].node=node[i].node;
			load[k].p=-node[i].Pl;
			load[k].q=-node[i].Ql;
			k++;
		}
		else//提取发电机和平衡节点信息
		{
			genrt[j].node=node[i].node;
			genrt[j].p=node[i].Pg-node[i].Pl;
			genrt[j].q=node[i].Qg-node[i].Ql;
			genrt[j].v=node[i].v;
			j++;
		}
		if(node[i].flag==-1)
		{
			ppv[l]=node[i].node;
			l++;
		}

	}
    number_genrt=j-1;
    number_load=k-1;
	number_pv=l-1;

    /*//读取pv节点号
	for(i=1;i<=number_pv;i++)
	{
		fscanf(in,"%d",&q1);
		ppv[i]=q1; 
	}*/
	
}
void nlf::node_order()//求取原始节点的优化节点
{
    int i,i_node,j_node,node;
	int n,nn;
	int amount1,amount2;
    DATA_II  *data_ii;
	data_ii=(DATA_II *)calloc((N+1),sizeof(DATA_II));
	if(!data_ii)EXIT;

    for(i=1;i<=number_branch;i++)
	{   
		i_node=abs(branch[i].i_node);
		j_node=abs(branch[i].j_node);
        data_ii[i_node].node=i_node;
		data_ii[i_node].amount++;
		data_ii[j_node].node=j_node;
		data_ii[j_node].amount++;
	}
	
	for(n=1;n<=N-1;n++)
	{
        amount1=data_ii[n].amount;
		for(nn=n+1;nn<=N;nn++)
		{
			amount2=data_ii[nn].amount;
			if(amount1>amount2)
			{
				node=data_ii[n].node;
				data_ii[n].node=data_ii[nn].node;
				data_ii[nn].node=node;

				data_ii[n].amount=amount2;
				data_ii[nn].amount=amount1;

				amount1=amount2;
			}
		}
	}
	for( i=1;i<=N;i++)
		n_o[i]=data_ii[i].node;
	for( i=1;i<=N;i++)
	{
		node=n_o[i];
		o_n[node]=i;             //主要目的就是求出o_n[]（用于表示原始节点被优化后的节点号）
	}
	free(data_ii);
}

void nlf::re_order()
{ 
	int i,old_node,new_node;
	for(i=1;i<=number_branch;i++) //支路节点优化
	{
		int i_node,j_node;
		i_node=abs(branch[i].i_node);
		j_node=abs(branch[i].j_node);
		branch[i].i_node=o_n[i_node];
		branch[i].j_node=o_n[j_node];
	}
	for(i=1;i<=number_genrt;i++)//发电机机节点优化
	{
		old_node=genrt[i].node;
		new_node=o_n[old_node];
		genrt[i].node=new_node;
	}
	for(i=1;i<=number_load;i++)//负荷节点优化
	{
		old_node=load[i].node;
		new_node=o_n[old_node];
		load[i].node=new_node;
	}
	for(i=1;i<=number_ground;i++)//接地支路节点优化
	{
		old_node=ground[i].node;
		new_node=o_n[old_node];
		ground[i].node=new_node;
	}
	for(i=1;i<=number_pv;i++)//制造pv节点标志
	{
		old_node=ppv[i];
		new_node=o_n[old_node];
		ppv[i]=new_node;
		pv[new_node]=1;
	}
	NS=o_n[NS];             //平衡节点优化
}
void nlf::branch_rank()
{
	int i,node,i_node,j_node,ii_node,jj_node;
	int n,nn;
	for(i=1;i<=number_branch;i++) //此时支路两端的节点号不能改变，能改变的只能是两端的命名标号
	{
		node=branch[i].i_node;
		if(branch[i].bk<0)
			branch[i].j_node=-branch[i].j_node;//根据初始输入把变压器bk端的j_node节点标为负
		i_node=abs(branch[i].i_node);
		j_node=abs(branch[i].j_node);
		if(i_node>j_node)
		{
			branch[i].i_node=branch[i].j_node;
			branch[i].j_node=node;            
		}                 //始终把变压器变比端节点号标为负,并且i_node端放小节点，j_node端放大节点
	}   
	for(n=1;n<=number_branch-1;n++)
	{
		i_node=abs(branch[n].i_node);
		j_node=abs(branch[n].j_node);
        for(nn=n+1;nn<=number_branch;nn++)
		{
			ii_node=abs(branch[nn].i_node);
			jj_node=abs(branch[nn].j_node);
			if(i_node>ii_node||(i_node==ii_node&&j_node>jj_node))
			{
				PointData r,x,bk;
				node=branch[n].i_node;
				branch[n].i_node=branch[nn].i_node;
				branch[nn].i_node=node;

				node=branch[n].j_node;
				branch[n].j_node=branch[nn].j_node;
				branch[nn].j_node=node;

				r=branch[n].r;
				branch[n].r=branch[nn].r;
				branch[nn].r=r;

				x=branch[n].x;
				branch[n].x=branch[nn].x;
				branch[nn].x=x;

				bk=branch[n].bk;
				branch[n].bk=branch[nn].bk;
				branch[nn].bk=bk;

				i_node=ii_node;
				j_node=jj_node;
			}
		}
	}
	for(i=1;i<=number_branch;i++)
		branch[i].number=i;
}

void nlf::form_ymatrix()
{
    int i,i_node,j_node,*add;
	add=(int *)calloc((N+1),sizeof(int));
	if(!add)EXIT;
	PointData gij,bij,z;

	for(i=1;i<=number_branch;i++)//串联支路导纳阵
	{   
        i_node=abs(branch[i].i_node);
		j_node=abs(branch[i].j_node);
		z=branch[i].r*branch[i].r+branch[i].x*branch[i].x;
		gij=branch[i].r/z;
		bij=-branch[i].x/z;
		if(branch[i].bk<0)
		{ 
			PointData bk;
			bk=fabs(branch[i].bk);
			y_matrix[i].Gij=-gij/bk;  //互导纳
			y_matrix[i].Bij=-bij/bk;
			y_matrix[i].number=j_node;
            
			if(branch[i].j_node<0)    //自导纳
			{
				yy_matrix[j_node].Gii+=gij/(bk*bk);
				yy_matrix[j_node].Bii+=bij/(bk*bk);

				yy_matrix[i_node].Gii+=gij;
				yy_matrix[i_node].Bii+=bij;
			}
			else
			{
				yy_matrix[i_node].Gii+=gij/(bk*bk);
                yy_matrix[i_node].Bii+=bij/(bk*bk);

				yy_matrix[j_node].Gii+=gij;
				yy_matrix[j_node].Bii+=bij;
			}
			add[i_node]++;
			continue;
		}
        y_matrix[i].Gij=-gij;
		y_matrix[i].Bij=-bij;
		y_matrix[i].number=j_node;

		yy_matrix[i_node].Gii+=gij;
		yy_matrix[i_node].Bii+=bij+branch[i].bk;
		yy_matrix[j_node].Gii+=gij;
		yy_matrix[j_node].Bii+=bij+branch[i].bk;

		add[i_node]++;
	}
	for(i=1;i<=number_ground;i++)//并联支路导纳阵
	{
		int node;
		node=ground[i].node;
		yy_matrix[node].Bii+=ground[i].x;
	}
	yy_matrix[1].addres=1;
	for(i=1;i<=N;i++)
		yy_matrix[i+1].addres=yy_matrix[i].addres+add[i];
	free(add);
	//yy_matrix[i].address表示第i行第一个非零对角元素在总的非零对角元素中排列的序号

	using namespace std;     //输出到test.txt
	ofstream out_stream;
	out_stream.open("test.txt");
	if(out_stream.fail())
		puts("test file open error!\n");
	
	for(i=1;i<30;i++)
	{
		out_stream<<i<<" "<<yy_matrix[i].addres<<" "<<yy_matrix[i].Gii<<" "<<yy_matrix[i].Bii<<endl;
		out_stream<<i<<" "<<y_matrix[i].number<<" "<<y_matrix[i].Gij<<" "<<y_matrix[i].Bij<<endl;
	}
	out_stream.close();
}   

void nlf::initia_pqv()//初始化各节点的有功、无功和电压2
{
	int i;
	for(i=1;i<=N;i++)
	{
		dpqv[i].p=0;
		dpqv[i].q=0;
		dpqv[i].v2=0;
		dpqv[i].dp=0;
		dpqv[i].dq=0;
		dpqv[i].dv2=0;

		voltage[i].e=1;
		voltage[i].f=0;

	}

	for(i=1;i<=number_genrt;i++)//对pv节点和ns节点的初始电压进行重新设置
	{   
		int node;
        node=genrt[i].node;	
		if(pv[node]==1||node==NS)
		{
			voltage[node].e=genrt[i].v;
			voltage[node].f=0;
		}
	}
	for(i=1;i<=number_genrt;i++)
	{
		int node;
		node=genrt[i].node;
		dpqv[node].p+=genrt[i].p;
		dpqv[node].q+=genrt[i].q;
        if(pv[node]==1)
			dpqv[node].v2=genrt[i].v*genrt[i].v;
	}
	for(i=1;i<=number_load;i++)
	{
		int node;
		node=load[i].node;
		dpqv[node].p+=load[i].p;
		dpqv[node].q+=load[i].q;
	}
	/*using namespace std;     //输出到test.txt
	ofstream out_stream;
	out_stream.open("test.txt");
	if(out_stream.fail())
		puts("test file open error!\n");
	
	for(i=1;i<N+1;i++)
	{
		out_stream<<i<<" "<<dpqv[i].p<<" "<<dpqv[i].q<<" "<<dpqv[i].v2<<" "<<dpqv[i].dp<<" "<<dpqv[i].dq<<" "<<dpqv[i].dv2<<" "<<
			voltage[i].e<<" "<<voltage[i].f<<" "<<endl;
	}
	out_stream.close();*/
} 
void nlf::initia_dpqv()
{
	int i,j;
	for(i=1;i<=N;i++)
	{
		ab[i].aii=0;
		ab[i].bii=0;
	}
	for(i=1;i<=N;i++)//按节点的先后处理dpqv修正量
	{
		if(i==NS)continue;
		int n,end,number;
		for(j=1;j<i;j++)//求取aii和bii （1---i的节点累加）
		{
			for(end=yy_matrix[j+1].addres-1,n=yy_matrix[j].addres;n<=end;n++)
			{
				if(y_matrix[n].number==i)
				{
					ab[i].aii+=y_matrix[n].Gij*voltage[j].e-y_matrix[n].Bij*voltage[j].f;
					ab[i].bii+=y_matrix[n].Gij*voltage[j].f+y_matrix[n].Bij*voltage[j].e;
			       // break;考虑双回路	
				}
			}
		}
		ab[i].aii+=yy_matrix[i].Gii*voltage[i].e-yy_matrix[i].Bii*voltage[i].f;//（自身累加）
		ab[i].bii+=yy_matrix[i].Gii*voltage[i].f+yy_matrix[i].Bii*voltage[i].e;
		
		for(end=yy_matrix[i+1].addres-1,n=yy_matrix[i].addres;n<=end;n++)//（i---N的节点累加）
		{
            number=y_matrix[n].number;
			ab[i].aii+=y_matrix[n].Gij*voltage[number].e-y_matrix[n].Bij*voltage[number].f;
			ab[i].bii+=y_matrix[n].Gij*voltage[number].f+y_matrix[n].Bij*voltage[number].e;
		}
        dpqv[i].dp=dpqv[i].p-(ab[i].aii*voltage[i].e+ab[i].bii*voltage[i].f);//求出dp
		if(pv[i]==1)
			dpqv[i].dv2=dpqv[i].v2-(voltage[i].e*voltage[i].e+voltage[i].f*voltage[i].f);//求出dv2
		else
			dpqv[i].dq=dpqv[i].q-(ab[i].aii*voltage[i].f-ab[i].bii*voltage[i].e);//求出dq
	}
}
void nlf::gene_matrix()//按节点边形成边消去
{
    PointData *UP,*UQ;
	UP=(PointData *)calloc((2*N+1),sizeof(PointData));//申请两个工作数组
	if(!UP)EXIT;
    UQ=(PointData *)calloc((2*N+1),sizeof(PointData));
	if(!UQ)EXIT;