nl_lf.c

自己编写的牛顿-拉夫逊法潮流计算程序

	x=p->b; 
	b=r*r+x*x;
	r/=b; 
	x/=-b; 
	b=p->c; 
	yr[i][i]+=r;
	yi[i][i]+=x;
	r/=b; 
	x/=b; 
	yr[j][j]+=r/b; 
	yi[j][j]+=x/b;
	yr[i][j]+=-r; 
	yi[i][j]+=-x; 
	yr[j][i]+=-r; 
	yi[j][i]+=-x; 
	}
	num_bran=0;
	for(i=1;i<=num_node-1;i++)
	for(j=i+1;j<=num_node;j++)
	if((yr[i][j]!=0) || (yi[i][j]!=0)) num_bran++;
}
/*        牛顿法潮流        */
void ntlf()
{

/*        定义局部变量        */

/*        ju              放雅可比阵消取过程中非零上三角元素的列号            */
/*        iu              放雅可比阵消取过程中各行非零上三角元素在ju中的位置  */
/*        uu              放雅可比阵消取过程中非零上三角元素的数值            */
/*        xrg             放雅可比阵对应与每一节点的实部行                    */
/*        xig             放雅可比阵对应与每一节点的虚部行                    */
/*        pol             按序放各负荷节点在load中的位置                      */
/*        pog             按序放各发电机节点在gene中的位置                    */
/*        pos             计PV节点无功所在行及平衡机有功和无功所在行的位置    */

	int nd,ne,nb,fa,i,j,k,i1,i2,j1,j2,iter,no_error;
	int *pog,*pol,*pos,*iu,*ju;
	double vi,di,vii,p,q,vij,dij,a,b,c,d,e,f,error;
	double *xrg,*xig,*mis,*uu;

	/*        申请变量空间        */
	pog=(int *)calloc(num_node+1,sizeof(double));
	if(pog==NULL)
	{
		printf("No memory in pog!\n"); 
		fprintf(fpo,"No memory in pog!\n"); 
		exit(0);
	}
	pol=(int *)calloc(num_node+1,sizeof(double));
	if(pol==NULL)
	{ 
		printf("No memory in pol!\n");
		fprintf(fpo,"No memory in pol!\n");
		exit(0); 
	}
	nd=2*num_node+1;
	pos=(int *)calloc(nd,sizeof(double));
	if(pos==NULL)
	{
		printf("No memory in pos!\n"); 
		fprintf(fpo,"No memory in pos!\n"); 
		exit(0);
	}
	xrg=(double *)calloc(nd,sizeof(double));
	if(xrg==NULL)
	{  
		printf("No memory in xrg!\n"); 
		fprintf(fpo,"No memory in xrg!\n"); 
		exit(0);  
	}
	xig=(double *)calloc(nd,sizeof(double));
	if(xig==NULL)
	{
		printf("No memory in xig!\n");  
		fprintf(fpo,"No memory in xig!\n");  
		exit(0); 
	}
	mis=(double *)calloc(nd,sizeof(double));
	if(mis==NULL)
	{ 
		printf("No memory in mis!\n"); 
		fprintf(fpo,"No memory in mis!\n"); 
		exit(0); 
	}
	iu=(int *)calloc(nd,sizeof(int));
	if(iu==NULL)
	{ 
		printf("No memory in iu!\n");   
		fprintf(fpo,"No memory in iu!\n");   
		exit(0);
	}
	ne=60*num_bran; 
	nb=num_bran/2;
	ju=(int *)calloc(ne,sizeof(int));
	if(ju==NULL)
	{
		printf("No memory in ju!\n");
		fprintf(fpo,"No memory in ju!\n");
		exit(0); 
	}
	uu=(double *)calloc(ne,sizeof(double));
	if(uu==NULL) 
	{ 
		printf("No memory in uu!\n");  
		fprintf(fpo,"No memory in uu!\n");  
		exit(0); 
	}
	
	/*        负荷和发电机节点在load和gene中的位置,并置电压初值        */
	for(i=1;i<=num_load;i++) 
		pol[load[i].i]=i;
	for(i=1;i<=num_node;i++) // 置各节点电压初值 
	{ 
		va[i]=0.0e0;
		vm[i]=1.0e0; 
	}
	for(j=1;j<=num_gene;j++)
	{
		i=gene[j].i; 
		pog[i]=j; 
		k=gene[j].j;
		if(k<=0)
		{
			vm[i]=gene[j].c; // 置PV和平衡机节点电压初值 
			pos[i+i]=1; // 置PV节点和平衡机无功平衡所在行位置
			if(k==0) 
				pos[i+i-1]=1; // 置平衡机节点有功平衡所在行位置
		}
	}

	/*        迭代开始        */
	iter=1;
	for(;;)
	{
	loop:
	/*		按节点循环        */
		for(error=0.0e0,fa=1,i=1;i<=num_node;i++)
		{
			/*		初始化		*/
			i2=i+i; 
			i1=i2-1; 
			di=va[i];  
			vi=vm[i]; 
			vii=vi*vi;
			zero(xrg,1,nd);
			zero(xig,1,nd);
			p=0.0e0;
			q=0.0e0;

			/*		计算ij元素		*/
			for(j=1;j<=num_node;j++)
			{
				if(j==i) continue;
				if((yr[i][j]==0) & (yi[i][j]==0)) continue;
				vij=vi*vm[j]; 
				dij=di-va[j]; 
				j2=j+j;  
				j1=j2-1;
				a=vij*yr[i][j]; 
				b=vij*yi[i][j];  
				c=cos(dij); 
				d=sin(dij);
				e=a*c+b*d; 
				f=a*d-b*c;
				/*Pi方程存在时，计算Jacobian元素Hij，Nij*/
				if(pos[j1]==0)  
				{  
					xrg[j1]=f; 
					xig[j1]=-e;
				}
				/*Qi方程存在时，计算Jacobian元素Jij，Lij*/
				if(pos[j2]==0) 
				{ 
					xrg[j2]=e; 
					xig[j2]=f;
				}
				p+=e; // 叠加P和Q
				q+=f;
			}

			/*		计算ii元素		*/
			a=yr[i][i]*vii; 
			b=yi[i][i]*vii;
			xrg[i1]=-q; 
			xrg[i2]=p+2.0e0*a;
			xig[i1]=p; 
			xig[i2]=q-2.0e0*b;
			p+=a; 
			q-=b;
			/*		叠加不平衡功率		*/
			pq_node[i][0]=p; 
			mis[i1]=-p;
			pq_node[i][1]=q; 
			mis[i2]=-q;
			
			/*		计算负荷		*/
			j=pol[i];
			if(j>0)
			{
				p=load[j].a;  
				q=load[j].b;
				/*		叠加不平衡功率		*/
				pq_load[i][0]=p; 
				mis[i1]-=p;
				pq_load[i][1]=q;
				mis[i2]-=q;
			}

			/*		计算发电机导致的不平衡量		*/
			j=pog[i];
			if(j>0)
			{
				k=gene[j].j; 
				mis[i1]+=gene[j].a;
				if(k>0)
					mis[i2]+=gene[j].b;
				else
				{
					mis[i2]=0.0e0;
					if(k==0) 
						mis[i1]=0.0e0;
				}
			}

			/*		取最大不平衡量		*/
			a=fabs(mis[i1]); 
			b=fabs(mis[i2]);
			if(a>error) 
			{  
				error=a; 
				no_error=i;  
			}
			if(b>error) 
			{  
				error=b; 
				no_error=i; 
			}
			
			/*		逐行消去		*/
			if(pos[i1]==0)
				eliminate(i1,xrg,mis,nd,&fa,iu,ju,uu);
			else  
				/*		空行同样记录行位置		*/
				iu[i1]=fa;
			if(pos[i2]==0) 
				eliminate(i2,xig,mis,nd,&fa,iu,ju,uu);
			else  
				iu[i2]=fa;
			/*		检查数组体积,如果不够则增加数组单元数		*/
			if(fa>=ne)
			{
				free(ju); 
				free(uu);
				ne+=nb;
				ju=(int *)calloc(ne,sizeof(int));
				if(ju==NULL)
				{ 
					printf("No memory in ju!\n");
					fprintf(fpo,"No memory in ju!\n");
					exit(0);
				}
				uu=(double *)calloc(ne,sizeof(double));
				if(uu==NULL) 
				{  
					printf("No memory in uu!\n");  
					fprintf(fpo,"No memory in uu!\n");  
					exit(0); 
				}
				/*		重新计算		*/
				goto loop;
			} // loop结束
		}
		/*		求解电压		*/
		for(i=nd-2;i>0;i--)
		{
			for(k=iu[i+1],j=iu[i];j<k;j++) 
				mis[i]-=uu[j]*mis[ju[j]];
		}
		for(i=1;i<=num_node;i++)
		{
			i2=i+i; 
			i1=i2-1;
			va[i]+=mis[i1]; 
			vm[i]*=1.0e0+mis[i2];
		}

		/*		输出误差情况		*/
		printf("  迭代次数:%2d    最大迭代误差: %12.8f     发生在节点:%5d\n",iter,error,no_error);
		fprintf(fpo,"  迭代次数:%2d    最大迭代误差: %12.8f     发生在节点:%5d\n",iter,error,no_error);

		/*		迭加次数自加		*/
		iter++;

		/*		收敛判断		*/
		if(error<=error_max) 
		{ 
			conv=0; 
			break;
		}
		if((iter>iter_max)||(error>1.0e6))
		{
			conv=1;
			break; 
		}
	}
	free(pog);
	free(pol);
	free(pos);
	free(xrg);
	free(xig);
	free(mis);
	free(iu);
	free(ju);
	free(uu);
}

/*        雅可比矩阵按行消去        */
void eliminate(m,x,s,n,fa,iu,ju,uu)
int m,n,*fa,*iu,*ju;
double *x,*s,*uu;
{
	int i,j,k;
	double aa,bb;
	iu[m]=*fa;
	/*		对待消区的第m行之第1至第(m-1)列进行消去		*/
	for(i=1;i<m;i++)
	{
		aa=x[i]; 
		if(aa==0.0e0)
			continue;
		for(k=iu[i+1],j=iu[i];j<k;j++)
			x[ju[j]]-=aa*uu[j];/*xrg/xig消去*/
		s[m]-=aa*s[i];/*功率消去*/
	}
	aa=x[m];
	/*		记消去后第m行的非零元素		*/
	for(i=m+1;i<n;i++)
	{
		bb=x[i]; 
		if(bb==0.0e0) 
			continue;
		ju[*fa]=i; 
		uu[*fa]=bb/aa; 
		(*fa)++;
	}
	s[m]/=aa; // 功率同除
}

/*        数组 *a 从 “start” 到 “end” 清零        */
void zero(double *a,int start,int end)
{
	int i;
	for(i=start;i<end;i++)  
		a[i]=0.0e0;
}

/*        定义两维实数数组        */
double **TwoArrayAlloc_D(int m,int n)
{
	double *x,**y;
	int i;
	x=(double *)calloc(m*n,sizeof(double));
	if(x==NULL)  {  printf("\nno memory in xd");   exit(0);  }
	y=(double **)calloc(m,sizeof(double *));
	if(y==NULL)  {  printf("\nno memory in yd");   exit(0);  }
	for(i=0;i<m;i++) y[i]=&x[n*i];
	return(y);
}

/*        释放两维实数数组        */
void TwoArrayFree_D(double **x)
{
	free(x[0]);
	free(x);
}