牛顿法程序.c

一个用牛顿法编的潮流程序

//////////////////////////////////////////////////////////////////////
//							牛顿极坐标法潮流	        			//
//文件输入格式：节点总数n(包括联络节点),支路数zls					//
//节点数(发电机和负荷)nb,接地电抗数mdk,迭代精度eps					//
//考虑负荷静特性标志kk2(0考虑),优化标志(0不优化)					//
//最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[]	//
//支路容纳zyk[],节点号nob[]及标志nobt[](0-PQ,-1-PV)					//
//发电机和负荷有功、无功pg[],qg[],pl[],ql[]							//
//电压v0[](pv节点输入实际值，PQ节点任输入一值)						//
//电抗节点号idk[],电抗值dkk[]										//
//////////////////////////////////////////////////////////////////////
#include <math.h>
#include <stdio.h>
#define NS 2000
#define NS2 NS*2
#define NS4 1000
#define ZS 3000
#define ZS2 ZS*2
#define DKS 200
#define N2 ZS*4
#define N3 ZS*8 + NS*4

FILE *fp1, *fp2;
int i, j, k, l, i1, i2, i3, j2, j3, k1, n1, ll;
int n, zls, nzls, nb, kk2, mdk, mpj, ifop, bnsopton, it1, t, dsd;
int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];
int nnew[NS4], old[NS], nob[NS], nobt[NS];
float eps, dsm, vmin, dph, dqh, af[3];
float r, x, yk, zf, gij, bij, v00, vi, vj, vij, ci, cj, cij, sn, cs;
float zr[ZS], zx[ZS], zyk[ZS], dkk[DKS], gii[NS], bii[NS], yg[ZS2], yb[ZS2];
float pg[NS], qg[NS], pl[NS], ql[NS], v0[NS], v[NS], va[NS], ve[NS], vf[NS];
float w[NS2], kg[3], b[NS2];
char inname[12], outname[12];
int newsort[NS4];
//	newsort[i]存放i对应的老号

void  initial();
void newton();
void  out();

void dataio();
void bnsopt();
void zlsort();
void nod();
void out();
void printo();
void printy();
void y2();
void ya0();
void yy0();
void bbhl();
void branch();
void newval();
void printc();
void iswap();
void swap();
void solution(float* bb, int* nbbs, int* ibb, float* ww);
void jcb(float* bb, int* nbbs, int* ibb);
void oldtonew();
int find(int k, int a[], int* z);
int iabs(int a);
int isgn(int a, int b);
void printf1(float* aa, int n);


void main()

{
	initial();
	newton();
	out();// out the node and branch data
}

int isgn(int a, int b)

{
//****	本函数功能返回值为a的绝对值b的符号	****//
//参数1提供值，参数2提供符号//

	if (b < 0)
		if (a > 0)
		a = -a;
	return a;
}

void dataio()  
{ 
//****	系统数据初始化	****//

	int i;

	af[0] = 0.6;  
	af[1] = 2.0;//af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。   
	v00 = 1.000;//系统平均电压 

	printf("\nplease input the name of data file\n");
	scanf("%s", inname);
	fp1 = fopen(inname, "r");
	printf("\nplease output the name of data file\n");
	scanf("%s", outname);
	fp2 = fopen(outname, "w");
	fscanf(fp1, "%d %d %d %d", &n, &zls, &nb, &mdk);    
	// the number of node ,branches, node
	fscanf(fp1, "%f %d %d %d %d", &eps, &kk2, &mpj, 
		&bnsopton, &it1); 
	//precision, swing node,sort the node,iteration numbers
	for (i = 1; i <= zls; i ++)
    { 
		fscanf(fp1, "%d %d", &izl[i], &jzl[i]);
		fscanf(fp1, "%f %f %f ", &zr[i], &zx[i], &zyk[i]);
    }
	for (i = 1; i <= nb; i ++)
	{   
		fscanf(fp1, "%d %d", &nob[i], &nobt[i]);
		fscanf(fp1, "%f %f %f %f %f", &pg[i], &qg[i], &pl[i], 
			&ql[i], &v0[i]);
	}
	for (i = 1; i <= mdk; i ++) 
	{ 
		fscanf(fp1, "%d %f", &idk[i], &dkk[i]); 
	}
	fclose(fp1);
}

int find(int k, int a[], int* z)

{
//****	本函数查找a[]中是否有fabs(k)有则返回0，无则返回1	****//
//参数1为待查找量，参数2待搜索数组，参数3返回k在a[]中的次序号//

	int i;
	for (i = 1; i <= n; i ++)
		if (iabs(k) == a[i])
		{
			*z = i;
			return 1;
		}
	return 0;
}

void oldtonew()

{
//****	本函数将输入数据中的节点号变成从1开始的连续节点号	****//

	int i, j, k, ii1, ii2, zls2, k1, k2, k3, k4, ip;
	zls2 = zls + zls;

	for (i = 1; i <= zls2; i ++)
		newsort[i] = 0;
	ii1 = 0;

	for (i = 1; i <= zls; i ++)
	{
		k = izl[i];
		if (!find(k, newsort, &ii2))
		{	
			ii1 ++;
			newsort[ii1] = iabs(k);
		}
		k = jzl[i];
		if (!find(k, newsort, &ii2))
		{	
			ii1 ++;
			newsort[ii1] = iabs(k);
		}
	}
	for (i = 1; i <= ii1-1; i ++)
	{
		for (j = i+1; j <= ii1; j ++)
		{
			if (newsort[i] > newsort[j])
			{
				k = newsort[i];
				newsort[i] = newsort[j];
				newsort[j] = k;
			}
		}
	}
	for (i = 1; i <= zls; i ++)
	{
		k = izl[i];
		if (find(k, newsort, &ii2))
		{	
			izl[i] = isgn(ii2, k);
		}
		else 
			printf("error!");
		k = jzl[i];
		if (find(k, newsort, &ii2))
		{	
			jzl[i] = isgn(ii2, k);
		}
		else
			printf("error!");
		printf("izl[%d] = %d, jzl[%d] = %d\n", i, izl[i], i, jzl[i]);
	}
	for (i = 1; i <= nb; i ++)
	{
		for (j = 1; j <= n; j ++)
			if (nob[i] == newsort[j])
			{
				nob[i] = j;
				break;
			}
		printf("nob[%d] = %d\n", i, nob[i]);
	}
	for (j = 1; j <= n; j ++)
	{
		if (mpj == newsort[j])
		{
			mpj = j;
			break;
		}
	}
	//电抗器节点号转变
	for (j = 1; j <= mdk; j ++)
	{
		for (i = 1; i <= n; i ++)
		{
			if (idk[j] == newsort[i])
			{
				idk[j] = i;
				break;
			}
		}
	}
}
void initial()

{ 
//****	本函数进行初始化工作	****//

	int i, k1;
	dataio();//输入原始数据
	oldtonew();//转化为新号
	if (bnsopton == 0)    //节点不优化，新节点号即为老节点号。
		for (i = 1; i <= n; i ++) 
		{
			old[i] = i;  
			nnew[i] = i;
		}   
	else  
		bnsopt();	 //节点优化	
	mpj = nnew[mpj];//mpj:平衡节点
	zlsort(nnew);	// sort the r,x and b
	for (i = 1; i <= mdk; i ++) 
	{ 
		k1 = idk[i]; 
		idk[i] = nnew[k1];
	}
	for (i = 1; i <= n; i ++)   
	{ 
		v[i] = v00; 
		va[i] = 0.0;
	}  //  所有节点的电压幅值初值都为1.000(v00)，电压相角初值都为0 。
	  // exchange the node before and after sort
	for (i = 1; i <= n; i ++)     
		yds[i] = 0;      // the immediate
	for (i = 1; i <= nb; i ++)  
	{ 
		k1 = nnew[nob[i]]; 
		yds[k1] = nobt[i]; 
	}
	for (i = 1; i <= n; i ++)     
		nobt[i] = yds[i];
	newval(pg);
	newval(qg);
	newval(pl);
	newval(ql);
	newval(v0);
	for (i = 1; i <= n; i ++)  //  nobt[] is type of node
		if (nobt[i] == -1) 
			v[i] = v0[i];      //   nob[] is serials numbe
	//nobt[] = -1: pv节点，v0[]存放的是最后一个节点数据，
	//对于pv节点，即为该点应维持的电压值。
	//nobt[] =  0: pq节点，v0[]存放的是最后一个节点数据，
	//对于pq节点，即为系统平均电压值。
	printo();
	//输出af[]、v00和节点排序后的支路、节点和
	//接地电抗数据(仅仅查看中间结果)
	ya0();//获得yds[]、ydz[]、列足码iy[]。( P407 )
}



void out()             

{ 
//****	本函数输出节点和支路数据	****//

	zlsort(old);    // recover the data if sorted
	nod();  // node data
	branch();   //branch data
	printc('-', 78);
	printc('*', 78);
	fprintf(fp2, "\n");
}

int iabs(int a)

{
//****	本函数求绝对值	****//

	a = (a >= 0) ? a : -a;
	return a;
}

void printf1(float* aa, int n)

{
//****	本函数输出aa[i],i=1-n	****//
	int i;
	for (i = 1; i <= n; i ++)
    { 
		if(i % 5 == 1) 
			fprintf(fp2, "\n");
		fprintf(fp2, "%9.5f", aa[i]);
    }
	fprintf(fp2, "\n\n");
}
void printf2(float* aa, float* bb, int n)

{
//****	本函数打印aa[i],bb[i]i=1-n	****//
	int i;
	for (i = 1; i <= n; i ++)
    { 
		if (i % 2 == 1) 
			fprintf(fp2, "\n");
		fprintf(fp2, "%5d%10.4f%10.4f", i, aa[i], bb[i]);
    }
}


// print the former n items of  colletion aa(int)
void printi(int* aa, int n)

{
//****	本函数输出aa[1]-aa[n]的值	****//
	int i;
	for(i = 1; i <= n; i ++)
    { 
		if(i % 10 == 1) 
			fprintf(fp2, "\n");
		fprintf(fp2, "%5d", aa[i]);
    }
}

void printc(char aa, int n)

{ 
//****	本函数输出n个aa字符	****//
	
	int i;
	fprintf(fp2, "\n");
	for (i = 1; i <= n; i ++) 
		fprintf(fp2, "%c", aa);
}

 /* print the matrix of Y, the diagnal item in G(real) and B(virtual),
 the undiagonal item in yg,yb and iy*/
void printy()

{ 
	fprintf(fp2, "\n*******GII,BII********");
	printf2(gii, bii, n);
	fprintf(fp2, "\n*******YYYYY********");
	for (i = 1; i <= nzls; i ++)
    { 
		if (i % 2 == 1) 
			fprintf(fp2,"\n");
      fprintf(fp2, "%10.4f%10.4f%8d", yg[i], yb[i], iy[i]);
    }
}

void newval(float* aa)

{ 
//****	本函数将旧号换成新号	****//

	int i, k1;
	for (i = 1; i <= n; i ++)  
		b[i] = 0.0;
	for (i = 1; i <= nb; i ++) 
	{ 
		k1 = nnew[nob[i]]; 
		b[k1] = aa[i]; 
	}
	for (i = 1; i <= n; i ++) 
		aa[i] = b[i];
}

void zlsort(int* nnew)	////	。

{
//****	本函数进行支路数据排序	****//
//小节点号放左边，大节点号放右边//
//左右皆按从小到大顺序排列//

	int ip, k1, k2, k3, k4;
	int i, j;	
	for (i = 1; i <= zls; i ++)
    { 
		k3 = izl[i];
		k4 = jzl[i];
		k1 = iabs(k3);
		k2 = iabs(k4);	    //	 原节点号。
		izl[i] = isgn(nnew[k1], k3);  // 新节点号。
		jzl[i] = isgn(nnew[k2], k4);
		k3 = izl[i];
		k4 = jzl[i];
		k1 = iabs(k3);
		k2 = iabs(k4);
		if (k1 > k2) 
		{ 
			izl[i] = k4; 
			jzl[i] = k3;
		}
    }
	for (i = 1; i <= zls - 1; i ++)
    { 
		ip = i;
		k1 = iabs(izl[i]);
		k3 = iabs(jzl[i]);
		for (j = i + 1; j <= zls; j ++)
		{ 
			k2 = iabs(izl[j]);
			k4 = iabs(jzl[j]);
			if(k2 < k1 || (k2 == k1 && k4 < k3)) 
			{ 
				ip = j; 
				k1 = k2; 
				k3 = k4;
			}
		}
		if(i != ip)
		{ 
			iswap(&izl[i], &izl[ip]);
			iswap(&jzl[i], &jzl[ip]);
			swap(&zr[i], &zr[ip]);
			swap(&zx[i], &zx[ip]);
			swap(&zyk[i], &zyk[ip]);
		}
    }
}

void bnsopt()
{
//****	节点优化	****//
	
	int ii1, ii2, zls2, nomax;
	int i, j, l, k1, k; 
	int temp;
	zls2 = zls + zls;
	for (i = 1; i <= zls2; i ++) 
		old[i] = nnew[i] = 0;//先清零。由此可知：NS4、NS必须大于2*zls。
	for (i = 1; i <= zls; i ++) 
	{
		old[i] = iabs(izl[i]); 
		old[i + zls] = iabs(jzl[i]);
	}
	//变压器节点号由正变负，old[]前zls个为左节点号，后zls个为右节点号。
	for (i = 1; i <= zls2; i ++)//	冒泡法排序。
    { 
		k1 = i + 1;
		for (j = k1; j <= zls2; j ++)
		if (old[i] > old[j]) 
			iswap(&old[i], &old[j]);
		//交换整数old[i]、old[j]。小节点号排在支路左侧。
    }
  
	nomax = old[zls2];//nomax 即是最大节点号。Iee30.dat ---- 30
	l = 1;
	for (i = 1; i <= n; i ++)
    { 
		ii1 = old[l];
		old[i] = ii1;
		for (j = l; j <= zls2; j ++)
		{ 
			ii2 = old[j];
			if (ii1 != ii2) 
			{
				l = j; 
				break;
			}
			nnew[i] ++;
		}
	}
	for (i = 1; i <= n - 1; i ++)
    { 
		for (j = i + 1; j <= n; j ++)
			if (nnew[j] <= nnew[i])
		if ((nnew[j] != nnew[i]) || (old[j] < old[i]))
		{  
			iswap(&old[i], &old[j]);
			iswap(&nnew[i], &nnew[j]);
		}
	}
	for (i = 1; i <= nomax; i ++) 
		nnew[i] = 0;
	for (i = 1; i <= n; i ++) 
	{
		j = old[i];  
		nnew[j] = i;
	}
}

void ya0()   //( P407 )

{
//****	本函数获得yds[]、ydz[]、列足码iy[]	****//

	int i, j, l, ll;
	for (i = 1; i <= n; i ++) 
		yds[i]=0;   //  yds[]存放各行非零非对角元素的个数。
	ll = 1;
	for (l = 1; l <= zls; l ++)
    { 
		i = iabs(izl[l]);
		j = iabs(jzl[l]);
		if (i == j) 
			continue;
		iy[ll] = j;
		if ((i != iabs(izl[l + 1])) || (j != iabs(jzl[l + 1]))) 
		{
			ll ++; 
			//ll统计总支路数(双回线算一条支路)
			yds[i] ++;
		}  
    }
	nzls = ll - 1; //总支路数(双回线算一条支路)
	ydz[1] = 1;
	for (i = 1; i <= n - 1; i ++) 
		ydz[i + 1] = yds[i] + ydz[i];   //由yds[]得ydz[]。
//  ydz[i]是第 i 行第一个非零非对角元素的首地址，
//	即在所有非零非对角元素中的顺序号。
	fprintf(fp2, "\n*******YDZ********"); 
	printi(ydz, n);
	fprintf(fp2, "\n*******YDS********");
	printi(yds, n);
	fprintf(fp2, "\n\n");
}

void y2()		

{ 
//****	本函数形成节点导纳阵，一次形成，不分B'和B''	****//	
	int j1;
	float r, x, yk, zf, gij, bij;
	int i, j, i1, l, ll;     
	for (i = 1; i <= n  ; i ++) 
	{ 
		gii[i] = 0.0;  
		bii[i] = 0.0;
	}//  导纳阵对角元(与节点一一对应)先清零。
	for (i = 1; i <= mdk; i ++) 
	{ 
		j = idk[i];    
		bii[j] = -1.0 / dkk[i];
	}//计算接地支路导纳，只影响导纳阵对角元(自导纳)。
	for (i = 1; i <= zls; i ++) 
	{ 
		yg[i] = 0.0;   
		yb[i] = 0.0;
	} //导纳阵非零非对角元(与支路一一对应)先清零。
	ll = 1;
	for (l = 1; l <= zls; l ++)
	{   
		i1 = izl[l];		//	支路左节点号。
		j1 = jzl[l];		//	支路右节点号。
		i = iabs(i1);		//	变压器支路有一节点号为负值。
		j = iabs(j1);
		if (i == j) 
			continue;       //   排除左、右节点号相等的情况。