pq_powerflow.c

电力系统潮流计算的程序

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

void yzb(int t, int* iu, double* u, double* di, int* nfd)

{
//****	本函数求因子表	****//

//参数1为标志(t=0 求B',t=1求B'')//
//参数2因子表上三角矩阵非零非对角元素的列足码
//参数3因子表上三角矩阵非零非对角元素的数值
//参数4因子表上三角矩阵对角元素
//参数5因子表上三角各行非零元素个数

	int i, j, k, i1, i2; 
	int jj, jj1, jj2, im, x, fd[NS];
	double ai, b[NS];
	nfd[1] = 1;
	for (i = 1; i <= n; i ++)
	{   
		//nobt[] 存放的是节点类型, 0: pq节点, -1: pv节点。                                                   
		if (((t != 1) || (nobt[i] != -1)) && i != mpj) // <---|
		{                                              //     |
			for (j = i + 1; j <= n; j ++)              //	  |
				b[j] = 0.0;                            //     | 
			b[i] = bii[i];                             //     |   
			if ((kk2 == 0) && (t == 1) && (nobt[i] != -1))// 存在(t == 1)的情况，不多余。   
				b[i] = b[i] + af[1] * ql[i] / v0[i] / v0[i];//af[1]
			i1 = ydz[i];
			i2 = ydz[i + 1] - 1;
			for (j = i1; j <= i2; j ++) 
			{
				k = iy[j]; 
				b[k] = yb[j];
			}
			b[mpj] = 0.0;
			if (t == 1)
			for (j = 1; j <= n; j ++)
				if (nobt[j] == -1) 
					b[j] = 0.0;
			i1 = i - 1;
			for (im = 1; im <= i1; im ++)
			{	  
	    		jj1 = nfd[im];
				jj2 = nfd[im + 1] - 1;
				for (jj = jj1; jj <= jj2; jj ++)
				{
					if(iu[jj] == i)
					{
		   				ai = u[jj] / di[im];
						for(k = jj; k <= jj2; k ++)
						{
							j = iu[k];
							b[j] = b[j] - ai * u[k];
						}
						break;
					}
				}
			}
			x = nfd[i];
			di[i] = 1.0 / b[i];
			ai = di[i];
			k = 0;
			i1 = i + 1;
			for (j = i1; j <= n; j ++)
			{
				if (fabs(b[j]) > 1.0e-15)
				{
					u[x] = b[j] * ai;
					iu[x] = j;
					k++;
					x++;
				}
			}	  
			fd[i] = k;
		} 
		else 
		{ 
			fd[i] = 0; 
			di[i] = 0.0;
		}
		nfd[i+1] = nfd[i] + fd[i];
	}
	fprintf(fp2, "\n********U*********");
	for (i = 1; i <= x; i ++)
	{ 
		if(i % 3 == 1) 
			fprintf(fp2, "\n");
		fprintf(fp2, "%10.5f%8i", u[i], iu[i]);
	}
	fprintf(fp2, "\n********DI********");
	printf1(di, n);
}

void printf1(double* 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 calc(int* iu, double* u, double* di, int* nfd, double* b)

{ 
//****	本函数利用因子表解线形方程组。(P417图F1-9)	****//
	
	double bi;
	int i, j, k, i1, i2;
	for (i = 1; i <= n; i ++) //前代过程。
    { 
		bi = b[i];
		i1 = nfd[i];
		i2 = nfd[i + 1];
		for (j = i1; j < i2; j ++) 
		{ 
			k = iu[j]; 
			b[k] = b[k] - bi * u[j]; 
		}
		b[i] = bi * di[i];
	}
	for (i = n; i >= 1; i --) //   回代过程。
    { 
		bi = b[i];
		i1 = nfd[i];
		i2 = nfd[i + 1] - 1;
		for (j = i2; j >= i1; j --) 
		{ 
			k = iu[j]; 
			bi = bi - b[k] * u[j];
		}
		b[i] = bi;
    }
}

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

{ 
//****	本函数形成节点导纳阵，一次形成，不分B'和B''	****//	
	int j1;
	double 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;       //   排除左、右节点号相等的情况。
		r = zr[l];
		x = zx[l];
		yk = zyk[l];	    //	zr[],zx[],zyk[]:支路三参数。		 
		zf = r * r + x * x;
		gij = r / zf;
		bij = -x / zf;
		if ((i1 > 0) && (j1 > 0))//	不是变压器支路。是一般支路。
		{   
			yg[ll] = yg[ll] - gij;
			yb[ll] = yb[ll] - bij;
			gii[i] = gii[i] + gij;
			bii[i] = bii[i] + bij + yk;
			gii[j] = gii[j] + gij;
			bii[j] = bii[j] + bij + yk; 
		}
		else	//	变压器支路。
		{ 
			if(j1 < 0) 
			{ 
				i=iabs(j1); 
				j=iabs(i1);
			}  
			//	若非标准变比在右(j)侧，则左、右互换，保证非标准变比在左侧。
			gii[j] = gii[j] + gij;
			bii[j] = bii[j] + bij;	//	标准变比侧。
			gii[i] = gii[i] + gij / yk / yk; 
			//	yk=zyk[]，对变压器支路指非标准变比(设在节点号为负的一侧)。
			bii[i] = bii[i] + bij / yk / yk;  
			//	非标准变比侧阻抗计算。
			yg[ll] = yg[ll] - gij / yk;
			yb[ll] = yb[ll] - bij / yk;
		}
		//◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆
		//注意此处处理双回线和导纳阵一不一样
		if((i1 != izl[l + 1]) || (j1 != jzl[l + 1])) 
			ll++;
	}
//	打印导纳矩阵。对角元实部为gii，虚部为bii，
//	非零非对角元实部为yb[]，虚部为yb[]，列足码为iy[]。
	fprintf(fp2, "*******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 printf2(double* aa, double* 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]);
    }
}

void jdgl(int kq0)   

{ 
//****	本函数计算节点功率(流程图见p421)	****//

	double ai, bi;      //ai,bi 是临时工作单元。
	double vi, vj, ci, sn, cs;
	int i, j, k, i1, i2;     
	for (i = 1; i <= n; i ++) 
		w[i]=0.0;// w[]存放节点功率,首先清零。
	for (i = 1; i <= n; i ++)
    { 
		vi = v[i];  // v[]存放节点电压幅值。
		ai = -bii[i]; //bii[]存放导纳阵对角元的虚部(gii[] + j*bii[])
		if (kq0 == 0) 
			ai = gii[i]; //gii[]存放导纳阵对角元的实部(gii[] + j*bii[])
		w[i] = w[i] + vi * vi * ai; 
		//对角元素的节点功率p = vi*vi*gii,q = vi*vi*bii
		if(i < n)  //导纳阵最后一行没有非对角元，故 i < n而不能i = n.