牛顿法程序.c

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

		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 bbhl(int kq0) 

{ 
//****	本函数计算各节点的功率误差，求最大功率误差dsm	****//
//****	和常数项b[i]。(程序框图见P423)	****//	

	int i;
	float vi, vj;
	float pl0, pg0;
	dsm = 0.0;   // dsm 即为最大功率误差。
	for (i = 1; i <= n; i ++)
	{
		vi = v[i]; //  v[]存放节点电压幅值。
		vj = v0[i];
		// v0[]存放节点初始电压幅值。v0[]存放的是最后一个节点数据。
		// vi[]
		// nopt[] = -1: pv节点，对于pv节点，即为该点应维持的电压值。
		// nopt[] =  0: pq节点，对于pq节点，即为系统平均电压值。
		// vj 此时表示的是节点正常电压的会定值。
		if (kq0 == 0) 
		{
			pl0 = pl[i];
			pg0 = pg[i]; 
		}   //  负荷p,发电机p
		else 
        {
			pl0 = ql[i]; 
			pg0 = qg[i]; 
		}   //  负荷q,发电机q
		if (kk2 == 0) 
			pl0 = pl0 * ((vi - vj) * af[kq0] / vj + 1.0);   
		//  考虑负荷静特性。(公式F1-39,P423)
		if (nobt[i] == -1 && kq0 == 1) 
			qg[i] = w[i] - pl0;      //  pv节点
		if (i == mpj && kq0 == 0) 
			pg[i] = w[i] - pl0;      //  平衡节点
		b[i] = pg0 + pl0 - w[i];
		//pv节点(nobt[] = -1)和平衡节点(mpj)不参与
		//求最大功率误差和常数项的运算
		if (((kq0 != 1) || (nobt[i] != -1)) && (i != mpj))
		{
			if (fabs(b[i]) > fabs(dsm)) 
			{ 
				dsm = b[i]; 
				dsd = i;
			}    // dsm 即为最大功率误差,dsd存放其对应的节点号。
			b[i] = b[i] / vi;      //  计算修正方程式的常数项。
		}
		else 
			b[i]=0.0;
		// pv节点(nobt[] = -1)和平衡节点(mpj)不参与
		// 求最大功率误差和常数项的运算。
	}
}

void nod()   //   

{
//****	输出节点数据和最小电压幅值、相角(角度)及其节点号	****//
//****	(程序框图见p426 F1-16)	****//

	float vi, ci;
	int i, j, oldnumber;
	printc('+', 72);
	fprintf(fp2, "\n%5s%8s%10s%11s%11s%11s%11s\n", "I", 
		"V", "CA", "PL", "QL", "PG", "QG");
	vmin = v[1];
	dsd = 1;
	for (i = 1; i <= n; i ++)
	{
		j = nnew[i];
		oldnumber = newsort[i];//转化为相应旧号
		ci = va[j] * 180.0/3.1416; //弧度转化为角度。
		vi = v[j];
		if (vi < vmin) 
		{ 
			vmin = vi; 
			dsd = j;
		}    //  vmin即为最小电压,dsd存放其对应的新节点号。
		fprintf(fp2, "\n%5d%11.5f%12.6f", oldnumber, vi, ci);
		fprintf(fp2, "%11.5f%11.5f%11.5f%11.5f", pl[j], ql[j],
			pg[j], qg[j]);
	}
	printc('-', 72);
}


void branch()    

{ 
//****	本函数输出支路数据。(程序框图见p428 F1-17)	****//

	int ii, jj;
	float r, x, yk, zf, vi, vj, ci, cj; 
	int i, j, l;
	float de, df, ei, ej, fi, fj, fii, fir, pij, pji, qij, qji;
	dph = 0.0;      //   统计系统有功网损。
	dqh = 0.0;      //   统计系统无功网损。
	fprintf(fp2, "\n%5s%5s%10s%12s%12s%12s\n ", "I", "J", "PIJ",
		"QIJ", "PJI", "QJI");
	for (i =1; i <= mdk; i ++)
    { 
		j = idk[i];
		dkk[i] = v[j] * v[j] / dkk[i];
    }
	for (l = 1; l <= zls; l ++)
    { 
		ii = iabs(izl[l]);    //  izl[]: 支路左节点号。
		jj = iabs(jzl[l]);    //  jzl[]: 支路右节点号。
		i=nnew[ii];
		j=nnew[jj];         //   转换为新节点号。
		ii = newsort[ii];
		jj = newsort[jj];//转化为相应旧号
		r = zr[l];
		x = zx[l];
		yk = zyk[l];
		vi = v[i];      //    v[]: 电压幅值。
		ci = va[i];     //   va[]: 电压相角。
		vj = v[j];
		cj = va[j];
		//支路左、右节点电压值由极坐标转换为直角坐标
		ei = vi * cos(ci); 
		fi = vi * sin(ci);   
		// ei: 支路左节点电压实部，fi: 支路左节点电压虚部。
		ej = vj * cos(cj); 
		fj = vj * sin(cj);    
		// ej: 支路右节点电压实部，fj: 支路右节点电压虚部。
		if ((izl[l] < 0) || (jzl[l] < 0)) //  变压器支路。
		{
			if (izl[l] < 0) 
			{ 
				ei = ei / yk; 
				fi = fi / yk;
			}     //  yk=zyk[l]
			else           
			{
				ej = ej / yk; 
				fj = fj / yk;
			}
			yk = 0.0 ;
		}
		de = ei - ej;
		df = fi - fj;
		zf = r * r + x * x;
		fii = (de * r + df * x) / zf;
		fir = (df * r - de * x) / zf;
		pij = fii * ei + fir * fi;
		qij = fii * fi - fir * ei;
		pji = -fii * ej - fir * fj;
		qji = -fii * fj + fir * ej;
		qij = qij - vi * vi * yk;
		qji = qji - vj * vj * yk;
		dph = dph + pij + pji;
		dqh = dqh + qij + qji;
		fprintf(fp2, "\n%5d%5d%12.5f%12.5f%12.5f%12.5f", ii, jj, pij,
			qij, pji, qji);
	}
	fprintf(fp2, "\n\n%8s%19s%18s%13s", "DPH", "DQH", "VMIN", "DSD");
	fprintf(fp2, "\n %10.5f%18.5f%18.5f%12d", dph, dqh, vmin, newsort[old[dsd]]);
	fprintf(fp2, "\n     \n%22sTHE END OF OUTPUT", " ");
}

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", izl[i], 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", 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%7.4f", idk[i], dkk[i]);
    }
}

void iswap(int* m, int* n)

{
//****	本函数交换m和n值	****//
	
	int p;
	p = *m; 
	*m = *n; 
	*n = p;
}

void swap(float* m, float* n)

{ 
//****	本函数交换m和n值	****//

	float p;
	p = *m; 
	*m = *n; 
	*n = p;
}

void yy0()
{
//****	形成导纳阵满阵	****//
	
	int i2, k1, k2, k3, k4, ip, nzls1, nzls2, jy[ZS2];
	int z;
	for (i =  1; i <= n; i ++)
		yds[i] = 0;
	for (i = 1; i <= n; i ++)
    { 
		i1 = ydz[i];
		i2 = ydz[i + 1] - 1;
		for (j = i1; j <= i2; j ++) 
			jy[j] = i;//存放行足码，即对应下三角阵的列足码
    }
	nzls1 = nzls + 1;
	nzls2 = nzls + nzls;
	for (i = nzls1; i <= nzls2; i ++)
    { 
		i1 = i - nzls;//i1从1开始，i从nzls + 1开始
		iy[i] = jy[i1];
		jy[i] = iy[i1];
		yg[i] = yg[i1];
		yb[i] = yb[i1];
    }//iy[]中存放满阵的列足码，jy[]中存放满阵的行足码
	
	fprintf(fp2, "\nzzzzzzzzzzzzzzzzzzzzzzzz\n");
	for (z = 1; z <= nzls2; z ++)
	{
		if ((z % 5) == 0)
			fprintf(fp2, "\n");
		fprintf(fp2, "i = %d %d\t", z, iy[z]);
		fprintf(fp2, "j = %d %d\t", z, jy[z]);
	}

	n1 = nzls2 - 1;
	//冒泡排序法，先按行后按列的顺序
	for (i = 1; i <= n1; i ++)
    { 
		i1 = i + 1;
		ip = i;
		k1 = jy[i];
		k3 = iy[i];
		for (j = i1; j <= nzls2; j ++)
		{
			k2 = jy[j];
			k4 = iy[j];
			if (k2 < k1 || (k2 == k1 && k4 < k3)) 
			{ 
				ip = j; 
				k1 = k2; 
				k3 = k4;
			}
		}
		if (i != ip)
		{ 
			iswap(&iy[i], &iy[ip]);
			iswap(&jy[i], &jy[ip]);
			swap(&yg[i], &yg[ip]);
			swap(&yb[i], &yb[ip]);
		}
    }
	//
	for (j = 1; j <= nzls2; j ++) 
	{ 
		i = jy[j]; 
		yds[i] ++;//统计每行的非零非对角元素个数
	}
	ydz[1] = 1;
	for (i = 1; i <= n; i ++)
		ydz[i + 1] = ydz[i] + yds[i];
	//ydz[]每行第一个非零非对角元素在非零元素中的顺序
}

void newton()

{
//****	本函数为牛顿法极坐标潮流	****//

	float dsm1, ww[NS2], bb[N3];
	int nbbs[NS2], ibb[N3];
	y2();
	yy0();
	t = 0;
	do
	{
		jcb(bb, nbbs, ibb);
		bbhl(0);
		dsm1 = dsm;
		for (i = 1; i <= n; i ++) 
			ww[i] = b[i] * v[i];//ww[i],deltaP,deltaQ
		for (i = 1; i <= n; i ++)
		{
			w[i] = w[i + n];
		}
		bbhl(1);
		if (fabs(dsm) > fabs(dsm1)) 
			dsm1 = dsm;
		printf("%d\t%f\n", t, dsm1);
		for (i = 1; i <= n; i ++) 
		{
			ww[i + n] = b[i] * v[i];//ww[i],deltaQ，和修正没有除电压对应
			//所以这里没有乘电压
		}
		solution(bb, nbbs, ibb, ww);
		for (i = 1; i <= n; i ++)
		{
			va[i] = va[i] - b[i];
			v[i] = v[i] - b[n + i] * v[i];
		}
		t ++;
		fprintf(fp2, "\n%d  %d  %f", t, dsd, dsm1);
	} while(fabs(dsm1) > eps && t <= it1);
}

void solution(float* bb, int* nbbs, int* ibb, float* ww)

{ 
//****	本函数求解方程组	****//

	int im, i0, jj, jj1, jj2, nn, x, nfd[NS], iu[N3];
	float bi, u[N3];
	double dd;
	x = 1;
	nfd[1] = 1;
	nn = n * 2;
	dd = 1.0;
	for (i = 1; i <= nn; i ++)
    {
		for (j = 1; j <= nn; j ++) 
			b[j] = 0.0;
		i1 = nbbs[i];
		i2 = nbbs[i + 1] - 1;
		for (j = i1; j <= i2; j ++) 
		{ 
			i0 = ibb[j]; 
			b[i0] = bb[j];
		}
		i1 = i - 1;
		for (im = 1; im <= i1; im ++)
		{ 
			if (fabs(b[im]) < 1.0e-15) 
				continue;
			jj1 = nfd[im];
			jj2 = nfd[im + 1] - 1;
			for (jj = jj1; jj <= jj2; jj ++) 
            {
				j  = iu[jj];
				b[j] = b[j] - b[im] * u[jj];
            }
			ww[i] = ww[i] - b[im] * ww[im];
        }
		for (j = i + 1; j <= nn; j ++) 
        { 
			if (fabs(b[j]) > 1.0e-15)//不等于0 
			{ 
				u[x] = b[j] / b[i];
				iu[x] = j;
				x ++;
            }
		}
		if (fabs(b[i]) > 1.0e-15) //不等于0
		{
			dd = dd * b[i];
			ww[i] = ww[i] / b[i];
		}
		nfd[i + 1] = x;
	}
	fprintf(fp2, "\ndd=%e", dd);
	for (i = 1; i <= nn; i ++) 
		b[i] = ww[i];
	for (i = nn; i >= 1; i --) 
    { 
		jj1 = nfd[i];
		jj2 = nfd[i + 1] - 1;
		bi = b[i];
		for (jj = jj2; jj >= jj1; jj --) 
        { 
			j = iu[jj];
			bi = bi - b[j] * u[jj];
        }
		b[i] = bi;
    }
}

void jcb(float* bb, int* nbbs, int* ibb)

{
//****	本函数形成雅可比矩阵	****//
//****	见《电力系统》(南京工学院p144)	****//
	
	float bb1, bb2, ah[NS], an[NS], aj[NS], al[NS];
	int i20, i2, i0, nbb[NS2];
	i1 = 0;
	i2 = nzls * 4 + n * 2 - 2;
	for (i = 1; i <= i2; i ++) 
	{ 
		ibb[i] = 0; 
		bb[i] = 0.0;
	}
	for (i = 1; i <= n + n + 1; i ++) 
	{ 
		nbb[i] = 0; 
		nbbs[i] = 0;
	}
	for (i = 1; i <= n; i ++)
    {
		w[i] = w[i + n] = 0.0;
		for (j = 1; j <= n; j ++) 
			ah[j] = an[j] = aj[j] = al[j] = 0.0;
		j3 = ydz[i];
		j2 = ydz[i + 1] - 1;
		for (j = j3; j <= j2; j ++)
		{	
			i0 = iy[j];
			vi = v[i];
			vj = v[i0];
			vij = vi * vj;
			ci = va[i];
			cj = va[i0];
			cij = ci - cj;
			sn = sin(cij);
			cs = cos(cij);
			bb1 = vij * (yg[j] * cs + yb[j] * sn);
			bb2 = vij * (yg[j] * sn - yb[j] * cs);
			w[i] = w[i] + bb1;//累加H
			w[i + n] = w[i + n] + bb2;//累加J
			if (i0 != mpj) 
				ah[i0] = -bb2;
			if (nobt[i0] != -1) 
				an[i0] = -bb1;
			if (i0 != mpj) 
				aj[i0] = bb1;
			if (nobt[i0] != -1) 
				al[i0] = -bb2;
		}
		ah[i] = w[i + n];//Hii
		aj[i] = -w[i];//Jii
		w[i] = w[i] + vi * vi * gii[i];
		w[i + n] = w[i + n] - vi * vi * bii[i];
		if (i == mpj) 
			continue;
		an[i] = -vi * vi * gii[i] - w[i];//Nii
		al[i] = vi * vi * bii[i] - w[i + n];//Lii
		for (j = 1; j <= n; j ++)
		{ 
			if (fabs(ah[j]) > 1.0e-10)//!=0
			{ 
				nbb[i] ++;//nbb[i]第i行的非零元素个数 
				i1 ++ ;//统计非零元素个数
				bb[i1] = ah[j];//非零元素值
				ibb[i1] = j;//列足码
			}
		}
		for (j = 1; j <= n; j ++ )
		{ 
			if (fabs(an[j]) > 1.0e-10)//!=0
			{ 
				nbb[i] ++;//nbb[i]第i行的非零元素个数 
				i1 ++;//统计非零元素个数
				bb[i1] = an[j];//非零元素值
				ibb[i1] = j + n;//列足码
			}
		}
		if (nobt[i] == -1) 
			continue;
		for (j = 1; j <= n; j ++)
		{ 
			if (fabs(aj[j]) > 1.0e-10)//!=0
			{
				nbb[i + n] ++;//nbb[i]第i行的非零aj个数
				i2 ++;//统计非零元素个数？
				bb[i2] = aj[j];//非零元素值
				ibb[i2] = j;//列足码
			}
		}
		for (j = 1; j <= n; j ++)
		{ 
			if (fabs(al[j]) > 1.0e-10)//nbb[i]第i行的非零aj个数
			{
				nbb[i + n] ++; //nbb[i]第i行的非零al个数
				i2 ++;//统计非零元素个数？
				bb[i2] = al[j];//非零元素值
				ibb[i2] = j + n;//列足码
			}
		}
    }

	nbbs[1] = 1;
	for (i = 1; i <= n + n + 1; i ++) 
		nbbs[i + 1] = nbbs[i] + nbb[i];
	//nbbs[]每行第一个非零元素在非零元素中的顺序
	i20 = nzls * 4 + n * 2 - 1;
	for (i = i20; i <= i2; i ++)
    { 
		i3 = i1 + i - i20 + 1;
		bb[i3] = bb[i];//非零元素值
		ibb[i3] = ibb[i];//列足码
    }
	for (i = i3 + 1; i <= i2; i ++) 
		bb[i] = 0.0;
	fprintf(fp2,"\n********BB********");
	printf1(bb, i3);
	printi(ibb, i3);
	printi(nbb, 70);
	printi(nbbs, 70);
}