csflow.cpp

本程序采用牛顿拉夫逊进行电力网络的潮流分布计算

					/* 计算式(4-89)中的Jij (=Nij) */
					jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];
				}
				else	/* i是PV结点 */
				{
					/* 计算式(4-89)中的Rij (=0) */
					jm[f2(2*i-1,2*j-1,n0)]=0;

					/* 计算式(4-89)中的Sij (=0) */
					jm[f2(2*i-1,2*j,n0)]=0;
				}
			}
		}
	}
	if(k!=1)
	{
		return;
	}
	
	/* 输出Jacoby矩阵 */
	fprintf(fp,"\n J                 MATRIX(C)");
	for(io=1;io<=n0;io+=5)
	{
		i1=(io+4)>n0?n0:(io+4);
		fprintf(fp,"\n");
		for(j=io;j<=i1;j++)
		{
			fprintf(fp,"%10d",j);
		}
		for(i=1;i<=n0;i++)
		{
			fprintf(fp,"\n%2d",i);
			for(j=io;j<=i1;j++)
			{
				fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);
			}
		}
	}
	fprintf(fp,"\n");
}

/**********************************************
 *     本子程序用选列主元素的高斯消元法求解组 *
 * 性方程组求各结点电压修正量,如打印参数K=1,则*
 * 输出增广矩阵变换中的上三角及电压修正量.如果*
 * 无唯一解,则给出信息,并停止程序运行.        *
 **********************************************/
void sevc ( float a[], int n0, int k, int n1)
{
	extern FILE *file4;
	FILE *fp;
	int i,j,l,n2,n3,n4,i0,io,j1,i1;
	float t0,t,c;
	if(file4==NULL) fp=stdout;
	else fp=file4;
	for(i=1;i<=n0;i++)
	{
		l=i;
		for(j=i;j<=n0;j++)
		{
		   if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) )
		   {
			   l=j; /* 找到这列中的最大元 */
		   }
		}
		if(l!=i)
		{	/* 行交换 */
			for (j=i;j<=n1;j++)
			{
				t=a[f2(i,j,n1)]; 
				a[f2(i,j,n1)]=a[f2(l,j,n1)]; 
				a[f2(l,j,n1)]=t;
			}
		}
		if (fabs(a[f2(i,i,n1)]-0)<1e-10)
		{	/* 对角元近似于0, 无解 */
			printf("\nNo Solution\n"); 
		    //exit (1);
		}

		t0=a[f2(i,i,n1)];
		for(j=i;j<=n1;j++)
		{
			/* 除对角元 */
			a[f2(i,j,n1)]/=t0;
		}
		if(i==n0)
		{   /* 最后一行，不用消元 */
			continue;
		}

		/* 消元 */
		j1=i+1;
		for(i1=j1;i1<=n0;i1++)
		{
			c=a[f2(i1,i,n1)];
			for(j=i;j<=n1;j++)
			{
				a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;
			}
		}
	}

	if(k==1)
	{	/* 输出上三角矩阵 */
		fprintf(fp,"\nTrianglar Angmentex Matrix ");
		for(io=1;io<=n1;io+=5)
		{
			i0=(io+4)>n1?n1:(io+4);
			fprintf(fp,"\n");
			fprintf(fp,"       ");
			for(i=io;i<=i0;i++)
			{
				fprintf(fp,"%12d",i);
			}
			for(i=1;i<=n0;i++)
			{
				fprintf(fp,"\n");
				fprintf(fp,"%2d",i);
				for(j=io;j<=i0;j++)
				{
				   fprintf(fp,"%15.6f", a[f2(i,j,n1)]);
				}
			}
		}
	}

	/* 回代求方程解 */
	n2=n1-2;
	for(i=1;i<=n2;i++)
	{
		n3=n1-i;
		for(i1=n3;i1<=n0;i1++)
		{
			n4=n0-i;
			a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];
		}
	}
	
	if(k!=1)
	{
		return;
	}

	/* 输出电压修正值 */
	fprintf(fp,"\nVoltage correction E(i), F(i) :");
	for(io=1;io<=n0;io+=4)
	{ 
		i1=(io+1)/2;
		i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);
		fprintf(fp,"\n");
		for(j=i1;j<=i0;j++)
		{
		   fprintf(fp,"%16d%16d",j,j);
		}
		i1 = 2*i0;
		fprintf(fp,"\n");
		for(i=io;i<=i1;i++)
		{
		   fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);
		}
	}
}

/****************************************************
 *   本子程序计算线路功率,平衡节点功率,PV节点无功功 *
 * 率及线路的功率损耗并输出.如选择参数K1=1,则表示输 *
 * 入为极座标.                                      *
 ****************************************************/
#define Pi 3.1415927/180
void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\
           int e1[],int s1[],float g1[],float b1[],float c1[],float c[],\
           float co[],float p1[],float q1[],float p2[],float q2[],float p3[],\
           float q3[],float p[],float q[],float v[],float angle[],int k1)
{
	extern FILE *file4;/**file6;*/
	FILE *fp;
	float t1,t2,cm,x,y,z,x1,x2,y1,y2;
	int i,i1,j,m1,ns,pos1,pos2,km,st,en;
	ns=n-1;
	if(file4==NULL)
	{
		fp=stdout;
	}
	else
	{
		fp=file4;
	}

	fprintf(fp,"\nTHE RESULT ARE:");
	if(k1==1)
	{
		for(i=0;i<n;i++)
		{
			angle[i]*=Pi;
			e[i]=v[i]*cos(angle[i]);
			f[i]=v[i]*sin(angle[i]);
		}
	}
	t1=0.0;t2=0.0;
	for(i=1;i<=n;i++)
	{
		pos1=f1(i);pos2=f2(n,i,n);
		t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];
		t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];
	}
	pos1=f1(n);
	p[pos1]=t1*e[pos1];
	q[pos1]=-t2*e[pos1];
	m1=m+1;
	for(i1=m1;i1<=ns;i1++)
	{
		t1=0;t2=0;
		for(i=1;i<=n;i++)
		{
			pos1=f1(i);pos2=f2(i1,i,n);
			t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];
			t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];
		}
		pos1=f1(i1);
		q[pos1]=f[pos1]*t1-e[pos1]*t2;
	}
	for(i=0;i<n; i++)
	{
		cm=co[i];
		if(cm!=0)
		{
			q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;
		}
	}
	fprintf(fp,"\nBUS DATA");
	fprintf(fp,"\nBUS     VOLTAGE      ANGLE(DEGS.)      BUS P          BUS Q");
	for(i=0;i<n;i++)
	{
		v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);
		x=e[i];
		y=f[i];
		z=y/x;
		angle[i]=atan(z);
		angle[i]/=Pi;
		fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[i],p[i],q[i]);
	}
	fprintf(fp,"\n LINE FLOW ");
	for(i=1;i<=l;i++)
	{
		pos1=f1(i);
		st=s1[pos1];
		en=e1[pos1];
		x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];
		x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];
		y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];
		y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];
		p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];
		q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];
		p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];
		q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];
		for(j=1;j<=n;j++)
		{
			cm=c[f2(j,i,l)];
			if(cm!=0.0)
			{
				km=1;
				if(en==j)
				{
					km=2;
				}
				if(km==1)
				{
					q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;
				}
				else
				{
					q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;
				}
			}
		}
		p3[pos1]=p1[pos1]+p2[pos1] ;
		q3[pos1]=q1[pos1]+q2[pos1] ;
		fprintf(fp,"\n%2d%8d%11d %13.6e %13.6e %13.6e %13.6e\n%10d%11d %13.6e %13.6e",\
				i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1],\
				e1[pos1],s1[pos1],p2[pos1],q2[pos1]);
	}
}	


#define L 5					//网络的支路总数
#define N 5					//节点数
#define M 3					//网络的PQ节点数
#define N0 2*(N-1)			//雅可比矩阵的行数
#define N1 N0+1
int _tmain(int argc, _TCHAR* argv[])
{
	int K=0;					//打印开关.K=1,则打印;否则,不打印
	int K1=0;					//子程序PLSC中判断输入电压的形式;K1=1,则为极座标形式.否则为直角坐标形式. 
	float D=1e-4;				//有功及无功功率误差的最大值（绝对值）.
	float *DD;					//误差量最大者
	float p=1.0;;
	DD=&p;
	float G[f2(N,N,N)+1];				//Ybus的电导元素(实部).  
	float B[f2(N,N,N)+1];				//Ybus的电纳元素(虚部).   
	float G1[L]={0.624024,0.754717,0.829876,0,0};	//第I支路的串联电导.      
	float B1[L]={-3.900156,-2.641509,-3.112033,-63.492063,-31.746032};	//第I支路的串联电纳

	//float C1[L]={0.25,0,0.25,3.023432,1.511716};			//C1(I)第I支路的pie型对称接地电纳.?????????
	float C1[L]={0.25,0,0.25,0,0};			//C1(I)第I支路的pie型对称接地电纳

	//float C[f2(N,L,L)+1]={0,0,0,0,0, 0,0,0,3.023432,0, 0,0,0,0,0,
	//				0,0,0,-6.198035, 0,0,0,0,0,-3.099017};	//C(I,J)第I节点J支路不对称接地电纳.  ???????
	float C[f2(N,L,L)+1]={0,0,0,0,0, 0,0,0,3.023432,0, 0,0,0,0,1.511716,
					0,0,0,-3.174603,0, 0,0,0,0,-1.587302};	//C(I,J)第I节点J支路不对称接地电纳

	float CO[N]={0,0,0,0,0};				//CO(I) :第I节点的接地电纳
	int S1[L]={1,1,2,2,3};					//第I支路的起始节点号
	int E1[L]={2,3,3,4,5};					//第I支路的终止节点号
	float P[N]={-1.6,-2,-3.7,5.0,0};				//第I节点的注入有功功率
	float Q[N]={-0.8,-1,-1.3,0,0};				//第I节点的注入无功功率
	float P0[N]={0};				//第I节点有功功率误差
	
	float V0[N]={0};				//第I节点(PV节点)的电压误差(平方误差)
	float V[N]={1,1,1,1.05,1.05};				//第I节点的电压幅值
	//float E[N]={1,1,1,1.05,1.05};				//第I节点的电压的实部
	float E[N]={1,1,1,1,1.05};				//第I节点的电压的实部
	float F[N]={0,0,0,0,0};				//第I节点的电压的虚部
	float JM[f2(N0,N0,N0)+1];			//Jacoby矩阵的第I行J列元素
	float A[f2(2*N-2,2*N-1,2*N-1)+1]={0};		//A(I,J)修正方程的增广矩阵,三角化矩阵的第I行J列元素,
										//运算结束后A矩阵的最后一列存放修正的解.  
	float P1[L];				//第I支路由S1(I)节点注入的有功功率
	float Q1[L];				//第I支路由S1(I)节点注入的无功功率
	float P2[L];				//第I支路由E1(I)节点注入的有功功率
	float Q2[L];				//第I支路由E1(I)节点注入的无功功率
	float P3[L];				//第I支路的有功功率损耗
	float Q3[L];				//第I支路的无功功率损耗
	float ANGLE[N];			//第I节点电压的角度
	float Q0[N]={0};			//第I节点有功功率误差
	int ii,jj;
	int cnt=0;				//迭代次数
	file4=fopen("cs_flow_1.txt","w");
	ybus(N,L,M,G,B,G1,B1,C1,C,CO,1,S1,E1);//计算电导和导纳参数
	dpqc(P,Q,P0,Q0,V,V0,M,N,E,F,1,G,B,DD);
	//while (arrayMax(N,Q0)>D || arrayMax(N,P0)>D){
	while (*(DD)>D){
		cnt++;
		jmcc(M,N,N0,E,F,G,B,JM,1);//算jacob矩阵
		//算增广的修正方程矩阵，最后一列放误差向量，与书上的方程顺序稍有不同
		for (ii=1;ii<=N0;ii++){
			for(jj=1;jj<=N0;jj++){
				A[f2(ii,jj,N1)]=JM[f2(ii,jj,N0)];
			}
		}
		for(ii=1;ii<=N0;ii++){
			if(ii<=2*M){
				if(ii%2==1)		A[f2(ii,2*N-1,2*N-1)]=Q0[f1((ii+1)/2)];
				else	A[f2(ii,2*N-1,2*N-1)]=P0[f1(ii/2)];
			}
			else{
				if(ii%2==1)		A[f2(ii,2*N-1,2*N-1)]=V0[f1((ii+1)/2)];
				else	A[f2(ii,2*N-1,2*N-1)]=P0[f1(ii/2)];
			}
		}
		sevc(A,N0,1,N1);//解线性方程组
		for(ii=1;ii<=N0;ii++){
			if(ii%2==1){
				E[f1((ii+1)/2)]-=A[f2(ii,2*N-1,2*N-1)];
			}
			else F[f1(ii/2)]-=A[f2(ii,2*N-1,2*N-1)];
		}//修正解
		dpqc(P,Q,P0,Q0,V,V0,M,N,E,F,1,G,B,DD);//重新计算误差
	}

	/*float Atest[f2(4,5,5)+1]={1,2,0,0,1, 0,1,0,0,2, 0,0,1,0,1, 2,0,0,1,1};
	sevc ( Atest, 4, 1, 5);*/
	plsc(  N,  L,  M,  G ,  B ,  E ,  F ,\
             E1 ,  S1 ,  G1 ,  B1 ,  C1 ,  C ,\
             CO ,  P1 ,  Q1 ,  P2 ,  Q2 ,  P3 ,\
             Q3 ,  P ,  Q ,  V ,  ANGLE , K1);
	//file4=fclose("cs_flow_1");
	extern FILE *file4;/**file6;*/
	/*FILE *fp;
	if(file4==NULL) fp=stdout;
	else fp=file4;
	fprintf(fp,"\n迭代总次数为:");
	fprintf(fp,"%d",cnt);*/
	return 0;
}