电力系统潮流计算的C语言实现

第一篇：电力系统潮流计算的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 #define ZS #define ZS2 1000 3000 ZS * 2

//NS4、NS必须大于2*zls。//最大支路数

// // // // // // // #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] = %dn”, 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] = %dn”, 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.3fn”, af[0], af[1], v00);printc('-', 78);fprintf(fp2, “nn

*******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, “nn*******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,“nn******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 filen”);scanf(“%s”, inname);fp1 = fopen(inname, “r”);printf(“nplease output the name of data filen”);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: %dt%dt%fn”, t, dsd, dsm);else

printf(“Q: %dt%dt%fn”, 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];} } 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 ++)

{

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(im = 1;im <= i1;im ++){

jj1 = nfd[im];

jj2 = nfd[im + 1]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, “nn”);}

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]1;

for(j = i2;j >= i1;j--)

{

k = iu[j];

bi = bi1;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]);

}

} }

void bnsopt(){ iswap(&jzl[i], &jzl[ip]);swap(&zr[i], &zr[ip]);swap(&zx[i], &zx[ip]);swap(&zyk[i], &zyk[ip]);//**** 节点优化 ****//

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 <= n1;//总支路数(双回线算一条支路)ydz[1] = 1;for(i = 1;i <= ngij;yb[ll] = yb[ll]gij / yk;yb[ll] = yb[ll]1;

//即为第 i 行的非零元素个数。

for(k = i1;k <= i2;k ++)//对第 i 行的所有非零元素进行操作。{

if(kq0!= 0)

{

ai =-yb[k];

bi = yg[k];}//yb[]存放导纳阵非对角元的虚部(yg[] + j*yb[])else

{

ai = yg[k];

bi = yb[k];}//yg[]存放导纳阵非对角元的实部(yg[] + j*yb[])j = iy[k];//iy[]存放的是列足码。

vj = vi * v[j];// v[]存放节点电压幅值。

ci = va[i]bi;//非对角元素 j 的功率

} }

void bbhl(int kq0)

{

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

int i;double vi, vj;double 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 *((vipl0;

// pv节点 if(i == mpj && kq0 == 0)

pg[i] = w[i]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)不参与

// 求最大功率误差和常数项的运算。

} }

node_output()

//

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

double vi, ci;int i, j, oldnumber;printc('+', 72);fprintf(fp2, “n%5s%8s%10s%11s%11s%11s%11sn”, “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_output()

{

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

int ii, jj;double r, x, yk, zf, vi, vj, ci, cj;

int i, j, l;double 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%12sn ”, “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 = eifj;

zf = r * r + x * x;

fii =(de * r + df * x)/ zf;

fir =(df * rfir * ei;

pji =-fii * ejvi * vi * yk;

qji = qji1;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]);

}

} for(i = 1;i <= nb;i ++){

k = nob[i];

nob[i] = newsort[k];

printf(“nob[%d] = %dn”, i, nob[i]);} for(i = 1;i <= nb-1;i ++){

for(j = i+1;j <= nb;j ++)

{

if(nob[i] > nob[j])

{

k = nob[i];

nob[i] = nob[j];

nob[j] = k;

}

} } printf(“nob[%d] = %dn”, i, nob[i]);}

void yy1()

{

//**** 本函数形成节点导纳阵(不包括接地支路)****// 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 <= 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;bij =-1/x;yg[ll] = yg[ll]bij;gii[i] = gii[i] + gij;bii[i] = bii[i] + bij;gii[j] = gii[j] + gij;bii[j] = bii[j] + bij;

if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1])))

ll++;} // 打印导纳矩阵。对角元实部为gii，虚部为bii，// 非零非对角元实部为yb[]，虚部为yb[]，列足码为iy[]。

fprintf(fp2, “*******GII(1),BII(1)********n”);printf2(gii,bii,n);}

void y3()

{

//**** 本函数形成节点导纳阵，追加接地支路 ****//

int j1;double r, x, yk, zf, gij, bij;int i, j, i1, l, ll, kk = 0;

for(i = 1;i <= mdk;i ++)

{

j = idk[i];

bii[j] = bii[j]gij / yk;

yb[ll] = yb[ll]-bij / yk;

}

if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1])))

{

ll++;

kk = 0;

}

else

kk = 1;} // 打印导纳矩阵。对角元实部为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]);

} }

第二篇：电力系统通用潮流计算C语言程序

#include #include #include #include using namespace std;

//节点号

类型

负荷有功

负荷无功母线数据（类型 1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS { int busno;int type;float Pd;float Qd;};

//发电机数据 节点号

有功发电

电压幅值 struct Generator { int busno;float Pg;float Vg;};

//支路信息 节点I 节点J R X B/2 k struct Line { int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2// void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){ double dP=0,dQ=0,dV=0;int i,j;for(i=0;i

double A=0,B=0;

for(j=0;j

A+=YG[i][j]*e[j]-YB[i][j]*f[j];

B+=YG[i][j]*f[j]+YB[i][j]*e[j];

}

dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i];

dP=P[i]-e[i]*A-f[i]*B;

W[2*i]=dP;

dQ=Q[i]-f[i]*A+e[i]*B;

if(type[i]==1)

W[2*i+1]=dQ;

else W[2*i+1]=dV;} } //Jacobi矩阵// void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){

int i,j;

for(i=0;i

for(j=0;j

if(i!=j){

if(type[i]==1){

Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]);

Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];

Ja[2*i+1][2*j]=Ja[2*i][2*j+1];

Ja[2*i+1][2*j+1]=-Ja[2*i][2*j];

}

else {

Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i];

Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i];

Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0;

}

}

else {

double a[50]={0},b[50]={0};

for(int k=0;k

a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]);

b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]);

Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i];

Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];

if(type[i]==1){

Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i];

Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i];

}

else {

Ja[2*i+1][2*j]=-2*e[i];

Ja[2*i+1][2*j+1]=-2*f[i];

}

}

}

}

} }

//高斯消元法解方程组函数// void gauss(double a[][101],int n)

{

int i,j,k;double c;

for(k=0;k

c=a[k][k];

for(j=k;j<=n;j++)a[k][j]/=c;

for(i=k+1;i

c=a[i][k];

for(j=k;j<=n;j++)a[i][j]-=c*a[k][j];

}

}

a[n-1][n]/=a[n-1][n-1];

for(k=n-2;k>=0;k--)

for(j=k+1;j

void main(){ ifstream fin;int N=0,GS=0,LD=0,ZLs=0;//节点数

发电机数

负荷数 BUS

*B;

Generator

*G;Line

*L;

//从文本中读入原始数据到数组中//

fin.open(“C:data.txt”);if(!fin)

{

cout<<“输入数据文件不存在！”<

getchar();

}

int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50];

int i,j,l;

支路数//

for(i=0;;i++){

fin>>m1[i];

if(m1[i]==0)break;

fin>>m2[i]>>m3[i]>>m4[i];

N++;} B

=new BUS[N];for(i=0;i

B[i].busno=m1[i];

B[i].type=m2[i];

B[i].Pd=m3[i];

B[i].Qd=m4[i];}

for(i=0;;i++){

fin>>m1[i];

if(m1[i]==0)break;

fin>>m4[i]>>m3[i];

GS++;} G

=new Generator[GS];for(i=0;i

G[i].busno=m1[i];

G[i].Pg=m4[i];

G[i].Vg=m3[i];}

for(i=0;;i++){

fin>>m1[i];

if(m1[i]==0)break;

fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i];

ZLs++;}

L

=new Line[ZLs];for(i=0;i

L[i].busi=m1[i];

L[i].busj=m2[i];

L[i].R=m3[i];

L[i].X=m4[i];

} L[i].B=m5[i];L[i].k=m6[i];

LD=N-GS;fin.close();

//节点导纳矩阵形成// double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i

YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;} } for(l=0;l

} for(i=0;i

for(j=i;j

K[i][j]=K[j][i];K[j][i]=1;

}

for(j=0;j

if(i!=j){

YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]);

YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]);

}

}

} //修正后// for(l=0;l

// 求A=e+f// double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i

C[2*i]=1;} else C[2*i]=V[i];

}

double W[100]={0},Ja[100][101]={0};

//调用Jacobi函数和高斯函数//

for(int t=1;t<10;t++){

for(i=0;i<2*N-2;i++){

e[i]=C[2*i];

f[i]=C[2*i+1];

}

fun1(YG,YB,e,f,type,N,W,P,Q,V);

double it=fabs(W[0]);

for(i=1;i<2*N-2;i++){

if(it

it=fabs(W[i]);j=i;

}

}

//中间迭代过程//

cout<

cout<

if(it<0.00001)break;

Jacobi(YG,YB,e,f,type,N,Ja);

for(i=0;i<2*N-2;i++){

Ja[i][2*N-2]=W[i];

}

//高斯消元法解方程//

gauss(Ja,2*N-2);

for(i=0;i<2*N-2;i++){

D[i]=-Ja[i][2*(N-1)];

C[i]+=D[i];

}

}

//平衡节点//

for(i=0;i

double a=0,b=0;

for(int j=0;j

a+=(YG[i][j]*e[j]-YB[i][j]*f[j]);

b+=(YB[i][j]*e[j]+YG[i][j]*f[j]);

}

P[i]=e[i]*a+f[i]*b;

Q[i]=f[i]*a-e[i]*b;

}

//支路//

double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k

i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B);

}

//全网数据// int high=1,low=1;

double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]);for(k=0;kVh){

Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;} if(sqrt(e[k]*e[k]+f[k]*f[k])

Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;} }

//输出数据到文件databak.txt//

ofstream fout;

fout.open(“C:databak.txt”);

fout<<“节点”<

fout<

for(i=0;i

fout<

}

for(j=0;j

i=G[j].busno-1;

fout<}

fout<<“支路 ”<

fout<

for(k=0;k

i=L[k].busi-1;j=L[k].busj-1;

fout<}

fout<<“全网数据”<

fout<

fout<fout.close();}

第三篇：电力系统潮流计算

南 京 理 工 大 学

《电力系统稳态分析》

课程报告

姓名

XX

学 号： 5*** 自动化学院 电气工程

基于牛顿-拉夫逊法的潮流计算例题编程报学院(系): 专

业: 题

目: 任课教师 硕士导师 告

杨伟 XX

2015年6月10号

基于牛顿-拉夫逊法的潮流计算例题编程报告

摘要：电力系统潮流计算的目的在于：确定电力系统的运行方式、检查系统中各元件是否过压或者过载、为电力系统继电保护的整定提供依据、为电力系统的稳定计算提供初值、为电力系统规划和经济运行提供分析的基础。潮流计算的计算机算法包含高斯—赛德尔迭代法、牛顿-拉夫逊法和P—Q分解法等，其中牛拉法计算原理较简单、计算过程也不复杂，而且由于人们引入泰勒级数和非线性代数方程等在算法里从而进一步提高了算法的收敛性和计算速度。同时基于MATLAB的计算机算法以双精度类型进行数据的存储和运算, 数据精确度高，能进行潮流计算中的各种矩阵运算，使得传统潮流计算方法更加优化。

一 研究内容

通过一道例题来认真分析牛顿-拉夫逊法的原理和方法（采用极坐标形式的牛拉法），同时掌握潮流计算计算机算法的相关知识，能看懂并初步使用MATLAB软件进行编程，培养自己电力系统潮流计算机算法编程能力。

例题如下：用牛顿-拉夫逊法计算下图所示系统的潮流分布，其中系统中5为平衡节点，节点5电压保持U=1.05为定值，其他四个节点分别为PQ节点，给定的注入功率如图所示。计算精度要求各节点电压修正量不大于10-6。

二 牛顿-拉夫逊法潮流计算 1 基本原理

牛顿法是取近似解x(k)之后，在这个基础上，找到比x(k)更接近的方程的根，一步步地迭代，找到尽可能接近方程根的近似根。牛顿迭代法其最大优点是在方程f(x)=0的单根附近时误差将呈平方减少，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点的电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成节点电压新的初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。2 基本步骤和设计流程图

形成了雅克比矩阵并建立了修正方程式，运用牛顿-拉夫逊法计算潮流的核心问题已经解决，已有可能列出基本计算步骤并编制流程图。由课本总结基本步骤如下：

1)形成节点导纳矩阵Y；

2)设各节点电压的初值，如果是直角坐标的话设电压的实部e和虚部f；如果是极坐标的话则设电压的幅值U和相角a；

3)将各个节点电压的初值代入公式求修正方程中的不平衡量以及修正方程的系数矩阵的雅克比矩阵；

4)解修正方程式，求各节点电压的变化量，即修正量； 5)计算各个节点电压的新值，即修正后的值；

6)利用新值从第（3）步开始进入下一次迭代，直至达到精度退出循环； 7)计算平衡节点的功率和线路功率，输出最后计算结果； ① 公式推导

② 流程图

三

matlab编程代码

clear;

% 如图所示1，2，3，4为PQ节点,5为平衡节点

y=0;

% 输入原始数据，求节点导纳矩阵

y(1,2)=1/(0.07+0.21j);

y(4,5)=0;y(1,3)=1/(0.06+0.18j);

y(1,4)=1/(0.05+0.10j);

y(1,5)=1/(0.04+0.12j);

y(2,3)=1/(0.05+0.10j);

y(2,5)=1/(0.08+0.24j);

y(3,4)=1/(0.06+0.18j);

for i=1:5

for j=i:5

y(j,i)=y(i,j);

end

end

Y=0;

% 求节点导纳矩阵中互导纳

for i=1:5

for j=1:5

if i~=j

Y(i,j)=-y(i,j);

end

end

end

% 求节点导纳矩阵中自导纳

for i=1:5

Y(i,i)=sum(y(i,:));

end

Y

% Y为导纳矩阵

G=real(Y);

B=imag(Y);% 输入原始节点的给定注入功率

S(1)=0.3+0.3j;

S(2)=-0.5-0.15j;

S(3)=-0.6-0.25j;

S(4)=-0.7-0.2j;

S(5)=0;

P=real(S);

Q=imag(S);

% 赋初值，U为节点电压的幅值，a为节点电压的相位角

U=ones(1,5);

U(5)=1.05;

a=zeros(1,5);

x1=ones(8,1);

x2=ones(8,1);

k=0;

while max(x2)>1e-6

for i=1:4

for j=1:4

H(i,j)=0;

N(i,j)=0;

M(i,j)=0;

L(i,j)=0;

oP(i)=0;

oQ(i)=0;

end

end

% 求有功、无功功率不平衡量

for i=1:4

for j=1:5

oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

end

oP(i)=oP(i)+P(i);

oQ(i)=oQ(i)+Q(i);

end

x2=[oP,oQ]';

% x2为不平衡量列向量

% 求雅克比矩阵

% 当i~=j时，求H,N,M,L

for i=1:4

for j=1:4

if i~=j

H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,j)=H(i,j);

M(i,j)=-N(i,j);

end

end

end

% 当i=j时，求H,N,M,L

for i=1:4

for j=1:5

if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))

end

end

N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);

L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);

end

J=[H,N;M,L]

% J为雅克比矩阵

x1=-((inv(J))*x2);

% x1为所求△x的列向量

% 求节点电压新值，准备下一次迭代

for i=1:4

oa(i)=x1(i);

oU(i)=x1(i+4)*U(i);

end

for i=1:4

a(i)=a(i)+oa(i);

U(i)=U(i)+oU(i);

end

k=k+1;

end

k,U,a

% 求节点注入功率

i=5;

for j=1:5

P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);

Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);

end

S(5)=P(5)+Q(5)*sqrt(-1);

S

% 求节点注入电流

I=Y*U'

四

运行结果

节点导纳矩阵

经过五次迭代后的雅克比矩阵

迭代次数以及节点电压的幅值和相角（弧度数）

节点注入功率和电流

五 结果分析

在这次学习和实际操作过程里：首先，对电力系统分析中潮流计算的部分特别是潮流计算的计算机算法中的牛顿-拉夫逊法进行深入的研读，弄明白了其原理、计算过程、公式推导以及设计流程。牛顿-拉夫逊法是求解非线性方程的迭代过程，其计算公式为FJX，式中J为所求函数的雅可比矩阵；X为需要求的修正值；F为不平衡的列向量。利用x(*)=x(k+1)+X(k+1)进行多次迭代，通过迭代判据得到所需要的精度值即准确值x(*)。六 结论

通过这个任务，自己在matlab编程，潮流计算，word文档的编辑功能等方面均有提高，但也暴漏出一些问题：理论知识储备不足，对matlab的性能和特点还不能有一个全面的把握，对word软件也不是很熟练，相信通过以后的学习能弥补这些不足，达到一个新的层次。

第四篇：电力系统潮流计算发展史

电力系统潮流计算发展史

对潮流计算的要求可以归纳为下面几点：

（1）算法的可靠性或收敛性（2）计算速度和内存占用量（3）计算的方便性和灵活性

电力系统潮流计算属于稳态分析范畴，不涉及系统元件的动态特性和过渡过程。因此其数学模型不包含微分方程，是一组高阶非线性方程。非线性代数方程组的解法离不开迭代，因此，潮流计算方法首先要求它是能可靠的收敛，并给出正确答案。随着电力系统规模的不断扩大，潮流问题的方程式阶数越来越高，目前已达到几千阶甚至上万阶，对这样规模的方程式并不是采用任何数学方法都能保证给出正确答案的。这种情况促使电力系统的研究人员不断寻求新的更可靠的计算方法。

在用数字计算机求解电力系统潮流问题的开始阶段，人们普遍采用以节点导纳矩阵为基础的高斯-赛德尔迭代法（一下简称导纳法）。这个方法的原理比较简单，要求的数字计算机的内存量也比较小，适应当时的电子数字计算机制作水平和电力系统理论水平，于是电力系统计算人员转向以阻抗矩阵为主的逐次代入法（以下简称阻抗法）。

20世纪60年代初，数字计算机已经发展到第二代，计算机的内存和计算速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗矩阵是满矩阵，阻抗法要求计算机储存表征系统接线和参数的阻抗矩阵。这就需要较大的内存量。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行计算，因此，每次迭代的计算量很大。

阻抗法改善了电力系统潮流计算问题的收敛性，解决了导纳法无法解决的一些系统的潮流计算，在当时获得了广泛的应用，曾为我国电力系统设计、运行和研究作出了很大的贡献。但是，阻抗法的主要缺点就是占用计算机的内存很大，每次迭代的计算量很大。当系统不断扩大时，这些缺点就更加突出。为了克服阻抗法在内存和速度方面的缺点，后来发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需存储各个地区系统的阻抗矩阵及它们之间的联络线的阻抗，这样不仅大幅度的节省了内存容量，同时也提高了节省速度。

克服阻抗法缺点的另一途径是采用牛顿-拉夫逊法（以下简称牛顿法）。牛顿法是数学中求解非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳矩阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流程序的计算效率。自从20世纪60年代中期采用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为直到目前仍被广泛采用的方法。

在牛顿法的基础上，根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造，得到了P-Q分解法。P-Q分解法在计算速度方面有显著的提高，迅速得到了推广。

牛顿法的特点是将非线性方程线性化。20世纪70年代后期，有人提出采用更精确的模型，即将泰勒级数的高阶项也包括进来，希望以此提高算法的性能，这便产生了保留非线性的潮流算法。另外，为了解决病态潮流计算，出现了将潮流计算表示为一个无约束非线性规划问题的模型，即非线性规划潮流算法。

近20多年来，潮流算法的研究仍然非常活跃，但是大多数研究都是围绕改进牛顿法和P-Q分解法进行的。此外，随着人工智能理论的发展，遗传算法、人工神经网络、模糊算法也逐渐被引入潮流计算。但是，到目前为止这些新的模型和算法还不能取代牛顿法和P-Q分解法的地位。由于电力系统规模的不断扩大，对计算速度的要求不断提高，计算机的并行计算技术也将在潮流计算中得到广泛的应用，成为重要的研究领域。

第五篇：电力系统潮流计算程序设计

电力系统潮流计算程序设计

姓名：韦应顺

学号：2011021052 电力工程学院

牛顿—拉夫逊潮流计算方法具有能够将非线性方程线性化的特点，而使用MATLAB语言是由于MATLAB语言的数学逻辑强，易编译。

【】【】1.MATLAB程序12

Function tisco %这是一个电力系统潮流计算的程序 n=input（‘n请输入节点数：n=’）； m=input(‘请输入支路数：m=’);ph=input(‘n请输入平衡母线的节点号：ph=’)； B1=input（‘n请输入支路信号：B1=’）;%它以矩阵形式存贮支路的情况，每行存贮一条支路 %第一列存贮支路的一个端点 %第二列存贮支路的另一个端点 %第三列存贮支路阻抗

%第四列存贮支路的对地导纳

%第五列存贮变压器的变比，注意支路为1 %第六列存贮支路的序号

B2=input(‘n请输入节点信息：B2=’)； %第一列为电源侧的功率 %第二列为负荷侧的功率 %第三列为该点的电压值

%第四列为该点的类型：1为PQ，2为PV节点，3为平衡节点 A=input（‘n请输入节点号及对地阻抗：A=’）； ip=input（‘n请输入修正值：ip=’）； %ip为修正值);Y=zeros（n）；

Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）；

Y(p,q)=Y(p,q);f=zeros（1，n）；

no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2；

End for i=1：n

G=real(Y);if A（i2）=0

B=imag(Y);p=A(i1)；

Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while

max(abs(De))>ip&max(abs(Dfi)>ip

t=t+1;

e=e+De；

f=f+Df;

[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)；

J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)；

[De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v);

jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m

a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3));

b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2;

c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%该子程序是用来求取Df for i=1:n

If

i=ph

C(i)=0;

D(i)=0;

For j=i:n

C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end

P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If

B2(i4)=2 p=2*i-1;

Df(p)=P(i)-P1;p=p+1;else p=2*i-1;

Df(p)=P(i)-P1;p=p+1;

Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no=[];end

function [De,Df]=hxf(J,Df,ph,n,no)%该子函数是为求取De Df DX=JDf;DX1=DX;

x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0;

For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else

DX=[DX1,0,0];End k=0;

[x,y]=size(DX);For i=1:2:x K=k+1;

Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i);

X2=G(i,j)*e(i)+B(i,j)*f(i);%该子程序是用来求取jacci矩阵

for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if

J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End

2实例验证 【例题】设有一系统网络结线见图1，各支路阻抗和各节点功率均已以标幺值标示于图1中，其中节点2连接的是发电厂，设节点1电压保持U1＝1.06定值，试计算其中的潮流分布，请输入节点数：n=5 请输入支路数：m=7 请输入平衡母线的节点号：ph=l 请输入支路信息：

BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 请输入节点信息：

B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 请输入节点号及对地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 请输入修正值：ip=0．000 0l

参考文献

[1]陈珩．电力系统稳定分析[M]．北京：中国电力出版社，2002：139—187．

[2]郑阿奇．MATLAB实用教程[M]．北京：电子工业出版社，2005：1-243．

[3] 束洪春，孙士云，等．云电送粤交商流混联系统全过 程动态电压研究[J】．中国电力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn锄ic Voltage stabil时of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4．

[4] 朱新立，汤涌，等．大电网安全分析的全过程动态仿 真技术[J】．电网技术，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn锄ic simulation for the stabilhy a眦lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28．

[5]Roytelm锄I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．发电机励磁系统数学模型及参 数对电网动态稳定性分析结果影响的研究[J】．继电 器，2007，35(21)：22-27．

SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g锄e翰to璐baScd ∞de诅iled excitation system models柚d parameterS t0 power铲id dyn锄ic stabil时【J】．Relay，2007，35(2 1)：22-27．

[7] 方思立，朱方．快速励磁系统对系统稳定的影响[J】．中 国电机工程学报，1986，6(1)：20．28．

FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞

excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28．

[8] 刘取．电力系统稳定性及发电机励磁控制[M】．北京： 中国电力出版社，2007．

LIU Qu．Power system S诅bility锄d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro诅ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996．

[10] 陈利芳，陈天禄．浅谈自并励励磁系统在大容量机组 中的应用【J】．继电器，2007，35(1)：8l培4． CHEN Li-f抽岛CHEN Tian—lIL Application of 辩l仁exci组tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84．

[11] 方思立，刘增煌，孟庆和．大型汽轮发电机自并励励 磁系统的应用条件【J]．中国电力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu锄g，MENG Qin争hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63．

[12]梁小冰，黄方能．利用EMTDC进行长持续时间过程 的仿真研究【J】．电网技术，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan争眦ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57．

[13]王卉，陈楷，彭哲，等．数字仿真技术在电力系统中 的应用及常用的几种数字仿真工具【J】．继电器，2004，32(21)：7l一75．

wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眦hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75．

[14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI时studies【s】．