flow.c

C语言编写的潮流计算程序

		    jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[pos1];
		}

		/* 计算式(4-90)中的Hii */
		jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];

		/* 计算式(4-90)中的Jii */
		jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];
		
		if(i>m) /* PV结点 */
		{
			/* 计算式(4-90)中的Rii */
			jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];

			/* 计算式(4-90)中的Sii */
		    jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];
		}

		/* 计算式(4-89) */
		for(j=1;j<=ns;j++)
		{
			if(j!=i)
			{
				/* pos1是式(4-89)中的i */
				pos1=f1(i);
				
				/* pos2是式(4-89)中的ij */
				pos2=f2(i,j,n);

				/* 计算式(4-89)中的Nij */
				jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1];

				/* 计算式(4-89)中的Hij */
				jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1];

				if(i<=m) /* i是PQ结点 */
				{
					/* 计算式(4-89)中的Lij (=-Hij) */
					jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];

					/* 计算式(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)]);
		}
	}
}


#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;
	FILE *fp;
	float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;
	int i,i1,j,m1,ns,pos1,pos2,km;
	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%17d%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]);
	}
}