pq_powerflow.c

电力系统潮流计算的程序

//////////////////////////////////////////////////////////////////////
//							PQ分解法潮流	        			    //
//文件输入格式：节点总数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		//NS4、NS必须大于2*zls。
#define ZS		3000		//最大支路数
#define ZS2		ZS * 2
#define DKS		200			//最大电抗器数
#define N2		ZS * 4
#define N3		ZS * 8 + NS * 4  

FILE *fp1, *fp2;
char inname[12], outname[12];
//	fp1输入数据文件指针	fp2输出文件指针
//	inname[]输入数据文件名	outname[]输出数据文件名
int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;
//	节点总数n(包括联络节点)   支路数(回路数)zls  节点数nb(发电机和负荷)  
//	接地电抗数mdk	精度eps   平衡节点号mpj  
//	节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)  
//	最大迭代次数it1	 最低电压或最大功率误差节点号dsd
//	负荷静特性标志(=0考虑负荷静特性)
//	支路数(双回线算一条支路)
int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];
//	izl[],jzl[],idk[]:分别存放左、右节点号和电抗器节点号。
//	yds[]存放各行非零非对角元素的个数。
//	ydz[i]是第 i 行第一个非零非对角元素的首地址，
//	即在所有非零非对角元素中的次序号
//  iy[]存放列足码。
int nnew[NS4], old[NS], nob[NS], nobt[NS];
//	nnew[],old[]存放的是新、旧节点号。
//	nnew[i]中为i对应的新号
//	nob[]存放的是节点号。nobt[]存放的是节点类型, 0: pq节点, -1: pv节点。
double eps, dsm, vmin, dph, dqh, af[3];
//	eps迭代收敛精度，dsm最大功率误差
//	vmin:系统最低电压值。dph,dqh:系统有、无功损耗。
//	af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。
double v00;	
//  v00: 系统平均电压	ci,cj分别作为节点i,j的电压相角的临时存储单元。
double zr[ZS], zx[ZS], zyk[ZS], dkk[DKS], gii[NS], bii[NS], yg[ZS2], yb[ZS2];
double pg[NS], qg[NS], pl[NS], ql[NS], v0[NS], v[NS], va[NS];
//	支路电阻zr[]	支路电抗zx[]	输电线路充电容纳zyk[](y0/2)
//	接地电抗dkk[]	对角元实部gii[]	对角元虚部
//	非对角元实部yg[]	非对角元虚部yb[]
//	pg[],qg[],pl[],ql[]:发电机，负荷功率实、虚部
//  v[]是电压幅值,va[]是电压相角。
double w[NS2], kg[3], b[NS2];
int newsort[NS4];
//	newsort[i]存放i对应的老号

void  initial();
void pqflow();
void  out();
void dataio();
void bnsopt();
void zlsort(int* nnew);
void printo();
void printy();
void y2();
void ya0();
void yzb();
void jdgl(int kq0);
void bbhl(int kq0);
void calc();
int iabs(int a);
void branch_output();
void newval(double* aa);
void printc();
void iswap();
void swap();
void printf2(double* aa, double* bb, int n);
void calc(int* iu, double* u, double* di, int* nfd, double* b);
void printi(int* aa, int n);
void printf1(double* aa, int n);
int find(int k, int a[], int* z);
void yzb(int t, int* iu, double* u, double* di, int* nfd);
int isgn(int a, int b);
void yy1();
void y3();
void newtoold();


int main(void)

{
	initial();	//初始化
	pqflow();	//pq潮流计算
	out();		//输出节点和支路数据
	return 1;
}

int isgn(int a, int b)

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

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

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 printo()

{
//****	输出af[]、v00和节点排序后的支路、节点和接地电抗数据	****//

	int i;
	fprintf(fp2, "\n    ******AF AND V0 ******\n");
	fprintf(fp2, "\n  %7.3f%7.3f%7.3f\n", af[0], af[1], v00);
	printc('-', 78);
	fprintf(fp2, "\n\n    *******ZLB*******\n");
	for (i = 1; i <= zls; i ++)
    { 
		fprintf(fp2, "\n");
		fprintf(fp2, "%8d%8d%8d%8d", izl[i], jzl[i], old[abs(izl[i])], old[abs(jzl[i])]);
		fprintf(fp2, "%9.4f%9.4f%9.4f", zr[i], zx[i], zyk[i]);
    }
	printc('-', 78);
	fprintf(fp2, "\n\n*******BUS*******\n");
	for (i = 1; i <= nb; i ++)
    { 
		fprintf(fp2, "\n");
		fprintf(fp2, "%8d%8d%8d", nob[i],old[nob[i]], nobt[i]);
		fprintf(fp2, "%9.4f%9.4f%9.4f%9.4f%9.4f", pg[i], qg[i], pl[i], ql[i], v0[i]);
    }
	printc('-', 78);
	fprintf(fp2,"\n\n******DKK******\n");
	for (i = 1; i <= mdk; i ++)
    { 
		fprintf(fp2, "\n");
		fprintf(fp2, "%8d%8d%7.4f", idk[i], old[idk[i]], dkk[i]);
    }
}
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, "%lf %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, "%lf %lf %lf ", &zr[i], &zx[i], &zyk[i]);
    }
	for (i = 1; i <= nb; i ++)
	{   
		fscanf(fp1, "%d %d", &nob[i], &nobt[i]);
		fscanf(fp1, "%lf %lf %lf %lf %lf", &pg[i], &qg[i], &pl[i], 
			&ql[i], &v0[i]);
	}
	for (i = 1; i <= mdk; i ++) 
	{ 
		fscanf(fp1, "%d %lf", &idk[i], &dkk[i]); 
	}
	fclose(fp1);
}

void pqflow()	 

{ 
//****	PQ分解法计算潮流,程序框图见P164图3-16(从第 7 步起)	****//
	
	int kq0, iu1[N2], nfd1[NS], iu2[N2], nfd2[NS];
	int i, t;  
	double u1[N2], u2[N2], di1[NS], di2[NS];
	yy1();    
	yzb(0, iu1, u1, di1, nfd1); //form the B matrix of P-0 iteration
	y2();
	yzb(1, iu2, u2, di2, nfd2); //form the B matrix of Q-V iteration
	t = 0;     
	kq0 = 0;
	kg[0] = kg[1] = 1;
	do
	{
		jdgl(kq0); // calculating the power
		bbhl(kq0); // find out the maxi
		if (kq0 == 0)
			printf("P: %d\t%d\t%f\n", t, dsd, dsm);
		else
			printf("Q: %d\t%d\t%f\n", t, dsd, dsm);

		if (fabs(dsm) > eps)
		{ 
			kg[kq0]=1;
			if (kq0 == 0) 
				calc(iu1, u1, di1, nfd1, b);
			if (kq0 == 1) 
				calc(iu2, u2, di2, nfd2, b);
			for (i = 1; i <= n; i ++)
			{ 
				if(kq0 == 0 ) 
					va[i] = va[i] - b[i] / v00;
				else          
					v[i] = v[i] - b[i];
			}
		}
		else
			kg[kq0] = 0;
		if(kq0 == 0)   
			kq0 = 1;
		else         
		{ 
			kq0 = 0; 
			t ++;
		}
		if(t > it1) 
			break;
	}while((fabs(dsm) > eps) || (kg[kq0] != 0));	
	fprintf(fp2, "\n%s%d", "times = ", t);
}

void out()             

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

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

void newval(double* aa)

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

	int i, k1;
	for (i = 1; i <= n; i ++)  
		b[i] = 0.0;
	for (i = 1; i <= nb; i ++)