adc.c

本程序为c语言环境下的电力系统潮流计算程序 文件输入输出

		if(l!=i)
		{	/* 行交换 */
			for (j=1;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"); 
			
		}

		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 F(i), E(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,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起始节点：%3d终止节点：%3d起始节点注入线路功率:%10.6f+%10.6fj线路损耗：%10.6f+%10.6fj终止节点注入线路功率：%10.6f+%10.6fj",\
				i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1],\
				p2[pos1],q2[pos1]);
	}
}


/*形成增广矩阵子程序*/
void makeA(float *jm,float *a,float *p0,float *q0,float *v0,int m,int n0)
{
	int i,j1;
	int k;
	int n1;
	n1=n0+1;
	for(i=1,k=1;i<=n0-1;i+=2,k++)
	{
		for(j1=1;j1<=n0;j1++)
			a[f2(i,j1,n1)]=jm[f2(i,j1,n0)];
		a[f2(i,n1,n1)]=p0[f1(k)];
		for(j1=1;j1<=n0;j1++)
			a[f2(i+1,j1,n1)]=jm[f2(i+1,j1,n0)];
		if(k<=m)	
		    a[f2(i+1,n1,n1)]=q0[f1(k)];
		else
			a[f2(i+1,n1,n1)]=v0[f1(k)];
		

	}

}

/*从增广矩阵中提取解*/
void getsolution(float *a,float *e0,float *f0,int n0)
{
	int i,j,k,n1;
	n1=n0+1;
	for(i=1,k=1;i<=n0;i=i+2,k++)
	{
		f0[f1(k)]=a[f2(i,n1,n1)];
		e0[f1(k)]=a[f2(i+1,n1,n1)];
	}

	
}
float findmax(float *e0,float *f0,int n)/*寻找最大值*/
{
	float max;
	int i,pos1;
	max=fabs(e0[0]);
	for(i=1;i<=n;i++)
	{
		pos1=f1(i);
		if(fabs(e0[pos1])>max)max=fabs(e0[pos1]);
		if(fabs(f0[pos1])>max)max=fabs(f0[pos1]);
	}
	return max;
}




FILE *file4;
void main()
{
	/*声明变量*/
	int N=0,M=0,L=0;
	int N1,N0=0;
	int K,K1;
	int timesmax=20;
	int times=0;
	int i,j,pos1,pos2,pos;
	float D;
	float ep=1e-5;
	float tempC;
	float *G,*B,*G1,*B1,*C1,*C,*CO,*P,*Q,*P0,*Q0,*V0,*V,*E,*F,*JM,*A,*P1,*Q1,*P2,*Q2,*P3,*Q3,*ANGLE;
	float *E0,*F0;
	int *S1,*E1;
	char infile[]="e:\\0411\\3118041127\\in.txt";
	char outfile[]="e:\\0411\\3118041127\\out.txt";
	FILE *infp;
	pos=0;
   /*Read the input data*/
      if((infp=fopen(infile,"r"))==NULL)
        { printf("Can not open the file named 'in.txt' \n");
         
         }
   fscanf(infp,"%d,%d,%d",&N,&M,&L);
   N0=2*N-2;
   N1=N0+1;
	/*分配内存*/
	G=(float *)malloc(N*N*sizeof(float));
	B=(float *)malloc(N*N*sizeof(float));
	G1=(float *)malloc(L*sizeof(float));
	B1=(float *)malloc(L*sizeof(float));
	C1=(float *)malloc(L*sizeof(float));
	C=(float *)malloc(N*L*sizeof(float));
	CO=(float *)malloc(N*sizeof(float));
	S1=(int *)malloc(L*sizeof(int));
	E1=(int *)malloc(L*sizeof(int));
	P=(float *)malloc(N*sizeof(float));
	Q=(float *)malloc(N*sizeof(float));
	P0=(float *)malloc(M*sizeof(float));
	Q0=(float *)malloc(M*sizeof(float));
	V0=(float *)malloc(N*sizeof(float));
	V=(float *)malloc(N*sizeof(float));
	E=(float *)malloc(N*sizeof(float));
	F=(float *)malloc(N*sizeof(float));
	E0=(float *)malloc(N*sizeof(float));
	F0=(float *)malloc(N*sizeof(float));
	JM=(float *)malloc((2*N-2)*(2*N-2)*sizeof(float));
	A=(float *)malloc((2*N-2)*(2*N-1)*sizeof(float));
	P1=(float *)malloc(L*sizeof(float));
	Q1=(float *)malloc(L*sizeof(float));
	P2=(float *)malloc(L*sizeof(float));
	Q2=(float *)malloc(L*sizeof(float));
	P3=(float *)malloc(L*sizeof(float));
	Q3=(float *)malloc(L*sizeof(float));
	ANGLE=(float *)malloc(N*sizeof(float));
	/*分配结束*/
	/*初始化输入参数矩阵*/
	   for(i=1;i<=N;i++)
 	 	  {
 	 	  pos1=f1(i);
		  CO[pos1]=P[pos1]=Q[pos1]=E[pos1]=F[pos1]=V[pos1]=0.0;
		  for(j=1;j<=L;j++)
		    {pos2=f2(i,j,L);
		     C[pos2]=0;
 	 	    }
 	 	  }
	  for(i=1;i<=L;i++)
 	 	 {
 	 	  pos1=f1(i);
		  G1[pos1]=B1[pos1]=C1[pos1]=0.0;
		  S1[pos1]=E1[pos1]=0;
 	 	 }
  /*读入n个节点的数据*/
    for(i=1;i<=N;i++)
    {
    pos1=f1(i);
	if(i<=M)
		fscanf(infp,"%f,%f,%f,%d,%f",P+pos1,Q+pos1,CO+pos1,&pos,&tempC);
	else if(i>M&&i<N)
		fscanf(infp,"%f,%f,%f,%d,%f",P+pos1,V+pos1,CO+pos1,&pos,&tempC);
	else if(i=N)
		fscanf(infp,"%f,%f,%f,%d,%f",E+pos1,F+pos1,CO+pos1,&pos,&tempC);
	if(pos!=0)
		C[f2(i,pos,L)]=tempC;
    }
			
   /*读入l条支路的数据*/
    for(i=1;i<=L;i++)
    {
    pos1=f1(i);
    fscanf(infp,"%d,%d,%f,%f,%f",S1+pos1,E1+pos1,G1+pos1,B1+pos1,C1+pos1);
	}

    fclose(infp);
	
	 if((file4=fopen(outfile,"a"))==NULL)
        printf("Can not open the file named 'out.txt' \n");
	 
	/*生成节点导纳矩阵*/
	ybus(N,L,M,G,B,G1,B1,C1,C,CO,1,S1,E1);
	/*设置PQ,PV结点初始电压*/
	for(i=1;i<=N-1;i++)
	{  
		pos1=f1(i);
		if(i<=M)
		{
			E[pos1]=1.0;
			F[pos1]=0.0;
		}
		else  /*PV节点*/
		{
			E[pos1]=V[pos1];
			F[pos1]=0.0;
		}
	}

	/*初始化功率电压误差矩阵*/
	for(i=1;i<=N-1;i++)
		{
			pos1=f1(i);
			P0[pos1]=Q0[pos1]=V0[pos1]=0.0;
			E0[pos1]=F0[pos1]=0.0;

	}
	do
	{
		times++; /*记录循环次数*/
		/*修整电压值*/
		for(i=1;i<=N-1;i++)
		{
			pos1=f1(i);
			E[pos1]+=E0[pos1];
			F[pos1]+=F0[pos1];
			
		
		}		
	/*计算功率电压差值*/
		dpqc(P,Q,P0,Q0,V,V0,M,N,E,F,1,G,B,&D);
	/*计算jm矩阵*/
		jmcc(M,N,N0,E,F,G,B,JM,1);
		/*形成增广矩阵*/
		makeA(JM,A,P0,Q0,V0,M,N0);
		/*解方程组*/
		sevc(A,N0,1,N1);
		/*提取电压修整量*/
		getsolution(A,E0,F0,N0);
		
		/*寻找修整量的最大值*/
		D=findmax(E0,F0,N-1);
		fprintf(file4,"\nTIMES=%d,D=%f\n",times,D);
	}
	while(times<timesmax&&D>ep);
	/*输出结果*/
	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,0);

		fclose(file4);
}