zhxll.c

电力系统雅各比行列式计算,牛顿---拉夫逊进行潮流计算（c源程序）

				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)]);
		}
	}
}

/****************************************************
 *   本子程序计算线路功率,平衡节点功率,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,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]);
	}
}
void aaaaa(float *p0,float *q0,float *a,float *v0,int n,float *jm)
{int i,j;
 float m0;
 j=4*(n-1)*(n-1);
 for(i=0;i<j;i++)
 {a[i]=jm[i];}
 for(i=1;i<7;i++)
 { for(j=41;j>7*i-1;j--)
   {a[j]=a[j-1];}}
 a[6]=q0[0];
 a[13]=p0[0];
 a[20]=q0[1];
 a[27]=p0[1];
 a[34]=v0[2];
 a[41]=p0[2];}

void efc(float *e,float *f,float *a)
{ int i,j=0;
  float a0[6];
  for(i=0;i<6;i++)
   a0[i]=-a[(6+i*7)];
  for(i=0;i<3;i++)
   { e[i]+=a0[j];
     f[i]+=a0[j+1];
     j+=2;
   }
}


main()
{int n,m,l,k,z;
  float g[N][L],b[N][L];
  float g1[L]={0.588235,0,0.453858,0.480769},
	b1[L]={-2.352941,-3.666667,-1.891074,-2.403846},
	c1[L]={0.015280,0,0.019200,0.014130},
	c[N][L]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},
	co[N]={0.03448,0.02941,0,0.03333};
   int s1[L]={1,1,1,4},e1[L]={2,3,4,2};
  float p[N]={-0.300000,-0.550000,0.500000,0},
	q[N]={-0.180000,-0.130000,0,0},
	v[N]={0,0,1.100000,1.050000},
	e[N]={1,1,1.100000,1.050000},
	f[N]={0,0,0,0};
  float p0[N],q0[N],v0[N];
  float p1[L],q1[L],p2[L],q2[L],p3[L],q3[L],jm[2*N-2][2*N-2];
  float d,*D;
  int n0,n1,row,line,k1=0;
  float angle[N];
  float a[6][7];
    D=&d;
    n0=2*N-2;
    n1=n0+1;
    n=4;l=4;k=1;m=2;
    ybus(n,l,g,b,g1,b1,c1,c,co,k,s1,e1);
    scanf("%d",&z);
    dpqc(p,q,p0,q0,v,v0,m,n,e,f,k,g,b,D);
    jmcc(m,n,n0,e,f,g,b,jm,k);
printf("\n************************************************");
    scanf("%d",&z);
while(D[0]>0.00001)
  { k=0;
   jmcc(m,n,n0,e,f,g,b,jm,k);
   aaaaa(p0,q0,a,v0,n,jm);
  n1=n0+1;
  sevc(a,n0,k,n1);
  efc(e,f,a);
  dpqc(p,q,p0,q0,v,v0,m,n,e,f,k,g,b,D);
  }
  k=1;
  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);
printf("\n************************************************");
   scanf("%d",&z);
}