哈尔滨工业大学电力系统计算课程报告

电力系统计算课程报告哈尔滨工业大学电力系统计算课程报告基于P-Q法潮流计算分析一、P-Q分解法基本原理1、 基本思想把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。2、 算法根据牛顿法是从数学上严格推导的，并没有考虑电力系统这个具体对象的特点。在高压电力系统中有功功率潮流主要与各节点电压向量的角度有关，无功功率潮流则主要受各节点电压向量幅值的影响。P-Q法是对牛顿法基于电力系统特点的简化算法。第一步简化，有功无功分解开来进行迭代。第二步简化，把系数矩阵简化为在迭代过程中不变的对称矩阵。二、P-Q分解法计算过程1、 基本计算公式修正方程式：功率误差方程式：2、 程序原理框图三、程序源代码本程序使用C+语言实现，应用面向对象设计思想，运用了多文件编译技术，并且，程序中使用了静态节点优化算法、稀疏矩阵技术、因子表分解等算法，并考虑了负荷静特性。详细内容如附录所示四、算例描述本算法用电科院6机22节点算例进行验证，详细参数和网络拓扑图如下：1、 算例参数2、 网络拓扑图五、计算结果I V CA PL QL PG QG1 1.00000 0.000000 0.00000 0.00000 5.96312 1.735492 0.97384 -11.311893 0.00000 0.00000 6.00000 3.200003 1.00000 -27.459057 0.00000 0.00000 3.10000 3.147314 1.02190 -25.168610 0.00000 0.00000 1.60000 0.700005 1.04392 -28.205656 0.00000 0.00000 4.30000 3.340006 1.00000 -37.566566 0.00000 0.00000 -0.01000 0.916647 1.02697 -5.247170 0.00000 0.00000 0.00000 0.000008 0.96847 -19.788813 -2.87000 -1.44000 0.00000 0.000009 0.98081 -19.738201 -3.76000 -2.21000 0.00000 0.0000010 0.97985 -19.304344 0.00000 0.00000 0.00000 0.0000011 0.99050 -23.173771 0.00000 0.00000 0.00000 0.0000012 0.99359 -30.499034 0.00000 0.00000 0.00000 0.0000013 0.98198 -35.528090 0.00000 0.00000 0.00000 0.0000014 1.00020 -30.695861 0.00000 0.00000 0.00000 0.0000015 1.00103 -30.720282 0.00000 0.00000 0.00000 0.0000016 0.99308 -37.750976 -5.00000 -2.90000 0.00000 0.0000017 0.96911 -37.546641 0.00000 0.00000 0.00000 0.0000018 0.98360 -37.696751 -4.30000 -2.60000 0.00000 0.0000019 1.00776 -31.018334 -0.86400 -0.66200 0.00000 0.0000020 1.01475 -35.674534 -0.71900 -0.47400 0.00000 0.0000021 1.01546 -32.104404 -0.70000 -0.50000 0.00000 0.0000022 1.05792 -29.749725 -2.26500 -1.69000 0.00000 0.00000- I J PIJ QIJ PJI QJI 1 7 5.96312 1.73549 -5.96312 -1.15693 2 9 6.00000 3.20000 -6.00000 -2.14195 3 22 3.10000 3.14731 -3.10000 -2.90531 4 19 1.60000 0.70000 -1.60000 -0.51308 5 18 4.30000 3.34000 -4.30000 -2.31986 6 17 -0.01000 0.91664 0.01000 -0.88833 7 8 3.48120 0.74380 -3.35384 0.14533 7 9 2.48191 0.41314 -2.39368 0.21112 8 9 -0.28186 -0.83922 0.28470 0.85016 8 22 0.76571 -0.74611 -0.71214 0.59561 9 10 3.63860 -0.48476 -3.63860 0.45674 9 22 0.71038 -0.64457 -0.66894 0.39951 10 11 3.63860 -0.45674 -3.63860 0.70886 11 12 3.63860 -2.04951 -3.59392 -1.18599 12 13 3.38055 -1.10216 -3.35200 -1.32306 12 15 0.21337 -0.40993 -0.21337 0.41383 13 17 3.35200 1.32306 -3.35200 -1.18839 14 15 -0.21337 0.41339 0.21337 -0.41383 14 19 0.21337 -0.41339 -0.21264 0.41772 16 17 -3.34200 -2.06023 3.34200 2.07671 16 18 0.00027 0.28284 -0.00000 -0.28014 16 19 -0.51158 -0.08567 0.52751 -0.23199 16 20 -0.59323 -0.39926 0.59959 -0.04958 16 21 -0.55346 -0.14284 0.56509 -0.13266 19 21 0.42113 -0.33465 -0.41788 0.34519 20 22 -1.31859 -0.42442 1.35499 -0.07584 21 22 -0.84721 -0.71254 0.86108 0.29604六、附录- 22 -#include "math.h"#include "stdio.h"#define NS2000/最大节点数#define NS2NS * 2#define NS41000 /NS4、NS要大于2*zls。#define ZS3000/最大支路数#define ZS2ZS * 2#define DKS200/最大电抗器数#define N2ZS * 4#define N3ZS * 8 + NS * 4 FILE *fp1, *fp2;char inname12, outname12;/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 izlZS, jzlZS, idkDKS, ydsNS, ydzNS, iyZS2;/izl,jzl,idk:分别存放左、右节点号和电抗器节点号。/yds存放各行非零非对角元素的个数。/ydzi是第 i 行第一个非零非对角元素的首地址，/即在所有非零非对角元素中的次序号/ iy存放列足码。int nnewNS4, oldNS, nobNS, nobtNS;/nnew,old存放的是新、旧节点号。/nnewi中为i对应的新号/nob存放的是节点号。nobt存放的是节点类型, 0: pq节点, -1: pv节点。double eps, dsm, vmin, dph, dqh, af3;/eps迭代收敛精度，dsm最大功率误差/vmin:系统最低电压值。dph,dqh:系统有、无功损耗。/af0和af1分别是负荷有功功率、无功功率静态特性系数。double v00;/ v00: 系统平均电压ci,cj分别作为节点i,j的电压相角的临时存储单元。double zrZS, zxZS, zykZS, dkkDKS, giiNS, biiNS, ygZS2, ybZS2;double pgNS, qgNS, plNS, qlNS, v0NS, vNS, vaNS;/支路电阻zr支路电抗zx输电线路充电容纳zyk(y0/2)/接地电抗dkk对角元实部gii对角元虚部/非对角元实部yg非对角元虚部yb/pg,qg,pl,ql:发电机，负荷功率实、虚部/ v是电压幅值,va是电压相角。double wNS2, kg3, bNS2;int newsortNS4;/newsorti存放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) = ai)*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 +)newsorti = 0;ii1 = 0;for (i = 1; i <= zls; i +)k = izli;if (!find(k, newsort, &ii2)ii1 +;newsortii1 = iabs(k);k = jzli;if (!find(k, newsort, &ii2)ii1 +;newsortii1 = iabs(k);for (i = 1; i <= ii1-1; i +)for (j = i+1; j <= ii1; j +)if (newsorti > newsortj)k = newsorti;newsorti = newsortj;newsortj = k;for (i = 1; i <= zls; i +)k = izli;if (find(k, newsort, &ii2)izli = isgn(ii2, k);else printf("error!");k = jzli;if (find(k, newsort, &ii2)jzli = isgn(ii2, k);elseprintf("error!");printf("izl%d = %d, jzl%d = %dn", i, izli, i, jzli);for (i = 1; i <= nb; i +)for (j = 1; j <= n; j +)if (nobi = newsortj)nobi = j;break;printf("nob%d = %dn", i, nobi);for (j = 1; j <= n; j +)if (mpj = newsortj)mpj = j;break;/电抗器节点号转变for (j = 1; j <= mdk; j +)for (i = 1; i <= n; i +)if (idkj = newsorti)idkj = i;break;void initial() /*本函数进行初始化工作*/int i, k1;dataio();/输入原始数据oldtonew();/转化为新号if (bnsopton = 0) /节点不优化，新节点号即为老节点号。for (i = 1; i <= n; i +) oldi = i; nnewi = i; else bnsopt(); /节点优化mpj = nnewmpj;/mpj:平衡节点zlsort(nnew);/ sort the r,x and bfor (i = 1; i <= mdk; i +) k1 = idki; idki = nnewk1;for (i = 1; i <= n; i +) vi = v00; vai = 0.0; / 所有节点的电压幅值初值都为1.000(v00)，电压相角初值都为0 。 / exchange the node before and after sortfor (i = 1; i <= n; i +) ydsi = 0; / the immediatefor (i = 1; i <= nb; i +) k1 = nnewnobi; ydsk1 = nobti; for (i = 1; i <= n; i +) nobti = ydsi;newval(pg);newval(qg);newval(pl);newval(ql);newval(v0);for (i = 1; i <= n; i +) / nobt is type of nodeif (nobti = -1) vi = v0i; / 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", af0, af1, v00);printc('-', 78);fprintf(fp2, "nn *ZLB*n");for (i = 1; i <= zls; i +) fprintf(fp2, "n");fprintf(fp2, "%8d%8d%8d%8d", izli, jzli, oldabs(izli), oldabs(jzli);fprintf(fp2, "%9.4f%9.4f%9.4f", zri, zxi, zyki); printc('-', 78);fprintf(fp2, "nn*BUS*n");for (i = 1; i <= nb; i +) fprintf(fp2, "n");fprintf(fp2, "%8d%8d%8d", nobi,oldnobi, nobti);fprintf(fp2, "%9.4f%9.4f%9.4f%9.4f%9.4f", pgi, qgi, pli, qli, v0i); printc('-', 78);fprintf(fp2,"nn*DKK*n");for (i = 1; i <= mdk; i +) fprintf(fp2, "n");fprintf(fp2, "%8d%8d%7.4f", idki, oldidki, dkki); void dataio() /*系统数据初始化*/int i;af0 = 0.6; af1 = 2.0;/af0和af1分别是负荷有功功率、无功功率静态特性系数。 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, nodefscanf(fp1, "%lf %d %d %d %d", &eps, &kk2, &mpj, &bnsopton, &it1); /precision, swing node,sort the node,iteration numbersfor (i = 1; i <= zls; i +) fscanf(fp1, "%d %d", &izli, &jzli);fscanf(fp1, "%lf %lf %lf ", &zri, &zxi, &zyki); for (i = 1; i <= nb; i +) fscanf(fp1, "%d %d", &nobi, &nobti);fscanf(fp1, "%lf %lf %lf %lf %lf", &pgi, &qgi, &pli, &qli, &v0i);for (i = 1; i <= mdk; i +) fscanf(fp1, "%d %lf", &idki, &dkki); fclose(fp1);void pqflow() /*PQ分解法计算潮流,程序框图见P164图3-16(从第 7 步起)*/int kq0, iu1N2, nfd1NS, iu2N2, nfd2NS;int i, t; double u1N2, u2N2, di1NS, di2NS;yy1(); yzb(0, iu1, u1, di1, nfd1); /form the B matrix of P-0 iterationy2();yzb(1, iu2, u2, di2, nfd2); /form the B matrix of Q-V iterationt = 0; kq0 = 0;kg0 = kg1 = 1;dojdgl(kq0); / calculating the powerbbhl(kq0); / find out the maxiif (kq0 = 0)printf("P: %dt%dt%fn", t, dsd, dsm);elseprintf("Q: %dt%dt%fn", t, dsd, dsm);if (fabs(dsm) > eps) kgkq0=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 ) vai = vai - bi / v00;else vi = vi - bi;elsekgkq0 = 0;if(kq0 = 0) kq0 = 1;else kq0 = 0; t +;if(t > it1) break;while(fabs(dsm) > eps) | (kgkq0 != 0);fprintf(fp2, "n%s%d", "times = ", t);void out() /*本函数输出节点和支路数据*/zlsort(old); / recover the data if sorted/newtoold();node_output(); / node databranch_output(); /branch dataprintc('-', 78);printc('*', 78);fprintf(fp2, "n");void newval(double* aa) /*本函数将旧号换成新号*/int i, k1;for (i = 1; i <= n; i +) bi = 0.0;for (i = 1; i <= nb; i +) k1 = nnewnobi; bk1 = aai; for (i = 1; i <= n; i +) aai = bi;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, fdNS;double ai, bNS;nfd1 = 1;for (i = 1; i <= n; i +) /nobt 存放的是节点类型, 0: pq节点, -1: pv节点。 if (t != 1) | (nobti != -1) && i != mpj) / <-| / |for (j = i + 1; j <= n; j +) / |bj = 0.0; / | bi = biii; / | if (kk2 = 0) && (t = 1) && (nobti != -1)/ 存在(t = 1)的情况，不多余。 bi = bi + af1 * qli / v0i / v0i;/af1i1 = ydzi;i2 = ydzi + 1 - 1;for (j = i1; j <= i2; j +) k = iyj; bk = ybj;bmpj = 0.0;if (t = 1)for (j = 1; j <= n; j +)if (nobtj = -1) bj = 0.0;i1 = i - 1;for (im = 1; im <= i1; im +) jj1 = nfdim;jj2 = nfdim + 1 - 1;for (jj = jj1; jj <= jj2; jj +)if(iujj = i) ai = ujj / diim;for(k = jj; k <= jj2; k +)j = iuk;bj = bj - ai * uk;break;x = nfdi;dii = 1.0 / bi;ai = dii;k = 0;i1 = i + 1;for (j = i1; j <= n; j +)if (fabs(bj) > 1.0e-15)ux = bj * ai;iux = j;k+;x+; fdi = k; else fdi = 0; dii = 0.0;nfdi+1 = nfdi + fdi;fprintf(fp2, "n*U*");for (i = 1; i <= x; i +) if(i % 3 = 1) fprintf(fp2, "n");fprintf(fp2, "%10.5f%8i", ui, iui);fprintf(fp2, "n*DI*");printf1(di, n);void printf1(double* aa, int n)/*本函数输出aai,i=1-n*/int i;for (i = 1; i <= n; i +) if(i % 5 = 1) fprintf(fp2, "n");fprintf(fp2, "%9.5f", aai); 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 = bi;i1 = nfdi;i2 = nfdi + 1;for (j = i1; j < i2; j +) k = iuj; bk = bk - bi * uj; bi = bi * dii;for (i = n; i >= 1; i -) / 回代过程。 bi = bi;i1 = nfdi;i2 = nfdi + 1 - 1;for (j = i2; j >= i1; j -) k = iuj; bi = bi - bk * uj;bi = bi; void zlsort(int* nnew)/。/*本函数进行支路数据排序*/小节点号放左边，大节点号放右边/左右皆按从小到大顺序排列/int ip, k1, k2, k3, k4;int i, j;for (i = 1; i <= zls; i +) k3 = izli;k4 = jzli;k1 = iabs(k3);k2 = iabs(k4); / 原节点号。izli = isgn(nnewk1, k3); / 新节点号。jzli = isgn(nnewk2, k4);k3 = izli;k4 = jzli;k1 = iabs(k3);k2 = iabs(k4);if (k1 > k2) izli = k4; jzli = k3; for (i = 1; i <= zls - 1; i +) ip = i;k1 = iabs(izli);k3 = iabs(jzli);for (j = i + 1; j <= zls; j +) k2 = iabs(izlj);k4 = iabs(jzlj);if(k2 < k1 | (k2 = k1 && k4 < k3) ip = j; k1 = k2; k3 = k4;if(i != ip) iswap(&izli, &izlip);iswap(&jzli, &jzlip);swap(&zri, &zrip);swap(&zxi, &zxip);swap(&zyki, &zykip); void bnsopt()/*节点优化*/int ii1, ii2, zls2, nomax;int i, j, l, k1, k; int temp;zls2 = zls + zls;for (i = 1; i <= zls2; i +) oldi = nnewi = 0;/先清零。由此可知：NS4、NS必须大于2*zls。for (i = 1; i <= zls; i +) oldi = iabs(izli); oldi + zls = iabs(jzli);/变压器节点号由正变负，old前zls个为左节点号，后zls个为右节点号。for (i = 1; i <= zls2; i +)/冒泡法排序。 k1 = i + 1;for (j = k1; j <= zls2; j +)if (oldi > oldj) iswap(&oldi, &oldj);/交换整数oldi、oldj。小节点号排在支路左侧。 nomax = oldzls2;/nomax 即是最大节点号。Iee30.dat - 30l = 1;for (i = 1; i <= n; i +) ii1 = oldl;oldi = ii1;for (j = l; j <= zls2; j +) ii2 = oldj;if (ii1 != ii2) l = j; break;nnewi +;for (i = 1; i <= n - 1; i +) for (j = i + 1; j <= n; j +)if (nnewj <= nnewi)if (nnewj != nnewi) | (oldj < oldi) iswap(&oldi, &oldj);iswap(&nnewi, &nnewj);for (i = 1; i <= nomax; i +) nnewi = 0;for (i = 1; i <= n; i +) j = oldi; nnewj = i;void ya0() /( P407 )/*本函数获得yds、ydz、列足码iy*/int i, j, l, ll;for (i = 1; i <= n; i +) ydsi=0; / yds存放各行非零非对角元素的个数。ll = 1;for (l = 1; l <= zls; l +) i = iabs(izll);j = iabs(jzll);if (i = j) continue;iyll = j;if(i != iabs(izll + 1) | (j != iabs(jzll + 1) ll +; /ll统计总支路数(双回线算一条支路)ydsi +; nzls = ll - 1; /总支路数(双回线算一条支路)ydz1 = 1;for (i = 1; i <= n - 1; i +) ydzi + 1 = ydsi + ydzi; /由yds得ydz。/ ydzi是第 i 行第一个非零非对角元素的首地址，/即在所有非零非对角元素中的顺序号。fprintf(fp2, "n*YDZ*"); printi(ydz, n);fprintf(fp2, "n*YDS*");printi(yds, n);fprintf(fp2, "nn");void printi(int* aa, int n)/*本函数输出aa1-aan的值*/int i;for(i = 1; i <= n; i +) if(i % 10 = 1) fprintf(fp2, "n");fprintf(fp2, "%5d", aai); 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 +) giii = 0.0; biii = 0.0;/ 导纳阵对角元(与节点一一对应)先清零。for (i = 1; i <= mdk; i +) j = idki; biij = -1.0 / dkki;/计算接地支路导纳，只影响导纳阵对角元(自导纳)。for (i = 1; i <= zls; i +) ygi = 0.0; ybi = 0.0; /导纳阵非零非对角元(与支路一一对应)先清零。ll = 1;for (l = 1; l <= zls; l +) i1 = izll;/支路左节点号。j1 = jzll;/支路右节点号。i = iabs(i1);/变压器支路有一节点号为负值。j = iabs(j1);if (i = j) continue; / 排除左、右节点号相等的情况。r = zrl;x = zxl;yk = zykl; /zr,zx,zyk:支路三参数。 zf = r * r + x * x;gij = r / zf;bij = -x / zf;if (i1 > 0) && (j1 > 0)/不是变压器支路。是一般支路。 ygll = ygll - gij;ybll = ybll - bij;giii = giii + gij;biii = biii + bij + yk;giij = giij + gij;biij = biij + bij + yk; else/变压器支路。 if(j1 < 0) i=iabs(j1); j=iabs(i1); /若非标准变比在右(j)侧，则左、右互换，保证非标准变比在左侧。giij = giij + gij;bi