电力系统计算报告

1、外南虎万学电力系统计算报告院 （系）电气工程及自动化授课教师 胡林献姓名 张远实学号 学S106052P-Q 分解潮流法简述 P-Q 分解潮流法的基本原理、计算过程、计算框图。用 C 语言编制 P-Q 分解法潮流程序，并用电科院 6 机 22 节点系统加以验证。要求采用稀疏技术、 因子表技术和节点优化技术， 并考虑负荷静态特性。P-Q 分解潮流法的基本原理：P-Q 分解法即是基于采用极坐标形式表示的牛顿法， 其根据电力系统实际运行状态的线路参数 R/X 通常很小的情况，对求解修正量的修正方程系数矩阵加以简化，使其变为常数阵 ( 即所谓的等斜率) ，且P、 Q 迭代解耦。这样可减少每次迭代的计算

2、时间，提高计算速度，又不影响最终结果，因此是通常选用的一种方法。但在低电压配电网中，当线路R/X 比值很大时，可能出现不收敛情况，此时应考虑更换其它方法。计算过程：1、 形成有功迭代和无功迭代的简化雅克比矩阵B' 和 B ''2、 给定 PQ 节点的初值和各节点电压相角初值n3、 作有功迭代PiPisPiPisVij1Vj(Gijcos ijBij sinij ) ， 计算Pi(k)/Vi(k)， 解修正方程i(k 1)i(k)i(k) ， 得各节点电压相角的修正值。n4、 作无功迭代QiQisQiQisViVj(Gij sin ij Bij cos ij) ， 计算Q

3、i(k)/Vi(k) ， 解修正方程 Vi(k 1) Vi (k) Vi (k) ， 得各节点电压幅值的修正量。5、 返回第三步，继续迭代到满足要求为止。编辑版 word计算框图流程：编辑版 word输入原始数据*_ ' _ "形成矩阵B及B并进行三角分解设PQ节点电压初值，各节点电压相角初值制迭代计数k=0Kp1,Kq1用公式计算不平衡功率 p(k),计算 p(k)/Vi(k)max p(k)解修正方程求i(k)(k 1)(k)(k)置Kp 0Kq 0?置KQ 1用公式计算不平衡功率Qi(k),计算Qi(k) /Vi(k)解修正方程求Vi(k)max Q(k)Y置Kq 0V

4、i(k1)V的Vi(k)Kp 0?置Kp 1计算平衡节点功率及全部线路输出算例描述用电科院 6 机 22 节点算例进行验证，详细参数和网络拓扑图如下所示表 1 线路参数支路号首末端节点号支路电阻支路电抗对地电纳 /217-80.01060.07400.027-90.01470.10400.038-90.00340.01310.049-220.05590.21800.1954512-130.002450.02551.395614-190.00340.02000.0716-190.05780.21800.1887816-200.01630.06620.2353916-210.03740.17800

5、.1641016-180.00330.03330.01119-210.01140.03700.01220-220.02140.08590.30081321-220.01500.06070.2198148-220.05370.19000.16531511-120.00330.03431.08797表 2 变压器支路数据支路号首末端节点号电阻电抗变比17-10.00.01501.05029-20.00.02171.075322-30.00.01241.100419-40.00.06401.025518-50.00.03751.050617-60.00.03371.000710-90.0-0.002

6、1.000811-100.00.01801.000915-120.00.01801.0001017-130.00.01001.0001115-140.0-0.0021.0001216-170.00.00101.027表 3 并联电容器数据节点号电纳12-1.3665160.5017611-1.366512-1.3665表 4 母线功率数据母线发电机输出有功发电机输出无功负荷有功负荷无功15.96311.73550.000.00263.20.000.0000.0000.0054.33.340.000.006-0.011.00.000.0070.000.00

7、0.000.0080.000.002.871.4490.000.003.762.21100.000.000.000.00110.000.000.000.00120.000.000.000.00130.000.000.000.00140.000.000.000.00150.000.000.000.00160.000.005.02.9170.000.000.000.00180.000.004.32.6190.000.000.8640.662200.000.000.7190.474210.000.000.70.5220.000.002.2651.69表 5 无功可调母线数据母线电压幅值（标幺值）无功

8、下限值无功上限制11.0-51031.0-5561.0-56表 6 发电机参数母线暂态电抗转子惯性时阻尼系数发电机有功发电机有功间常数出力下限出力上限10.0150140.820.00316.520.038230.000.001.26.630.039679.500.001.58.2540.121015.680.0048039.200.001.025.6160.19762.620.000.21.1编辑版 word图1电科院6机22节点系统图计算结果I V CA PL QL PG QG1 1.00000 0.00000 0.00000 0.00000 5.96312 1.735

9、492 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

10、.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.5280

11、90 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 -27.750976 -5.00000 -2.90000 0.00000 0.0000017 0.96911 -37.546641 0.00000 0.00000 0.00000 0.00000编辑版word18 0.98360 -37.696751 -4.30000 -2.60000 0.000

12、00 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.101546 -32.104404 -0.70000 -0.50000 0.00000 0.0000022 1.05792 -29.749725 -2.26500 -1.69000 0.00000 0.00000I J PIJ QIJ PJI QJI1 7 5.96312 1.73549 -5.96312 -1.156932 9 6.00000 3

13、.20000 -6.00000 -2.141953 22 3.10000 3.14731 -3.10000 -2.905314 19 1.60000 0.70000 -1.60000 -0.513085 18 4.30000 3.34000 -4.30000 -2.319866 17 -0.01000 0.91664 0.01000 -0.888337 8 3.48120 0.74380 -3.35384 0.145337 9 2.48191 0.41314 -2.39368 0.211128 9 -0.28186 -0.83922 0.28470 0.850168 22 0.76571 -0

14、.74611 -0.71214 0.595619 10 3.63860 -0.48476 -3.63860 0.456749 22 0.71038 -0.64457 -0.66894 0.3995110 11 3.63860 -0.45674 -3.63860 0.7088611 12 3.38055 -2.04951 -3.59392 -1.1859912 13 3.38055 -1.10216 -3.35200 -1.3230612 15 0.21337 -0.40993 -0.21337 0.4138313 17 3.35200 1.32306 -3.35200 -1.1883914 1

15、5 -0.21337 0.41339 0.21337 -0.4138314 19 0.21337 -0.41339 -0.21264 0.4177216 17 -3.34200 -2.06023 3.34200 2.0767116 18 0.00027 0.28284 -0.00000 -0.2801416 19 -0.51158 -0.08567 0.52751 -0.2319916 20 -0.59323 -0.39926 0.59959 -0.0495816 21 -0.55346 -0.14284 0.56509 -0.1326619 21 0.42113 -0.33465 -0.41

16、788 0.3451920 22 -1.31859 -0.42442 1.35499 -0.0758421 22 -0.84721 -0.71254 0.86108 0.29604源程序#include "math.h"#include "stdio.h"#define NS2000/ 最大节点数#define NS2NS * 2#define NS41000/NS4、 NS 必须大于 2*zls。#define ZS#define ZS2#define DKS#define N2#define N33000/ 最大支路数ZS * 2200/ 最大电抗器

17、数ZS * 4ZS * 8 + NS * 4FILE *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 最低电压或最

18、大功率误差节点号 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:

19、 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 支

20、路电抗 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();vo

21、id 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(char aa, int n);/void printc();void iswap(int* m, int* n);/void iswap();/void swap();void swap(double* m, double* n);void printf2(double* aa, double* bb, int

22、 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();void node_output(); /int

23、 pq_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 =

24、 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)111 +;ne

25、wsortii1 = iabs(k);k = jzli;if (!find(k, newsort, &ii2)112 +;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

26、(ii2, k);elseprintf("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

27、 = %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;/hu 2013-1-8printc('-', 58);fprintf(fp2, "nn*newsort*n");for (i = 1; i <= n; i +)fprintf(fp2, &

28、quot;%8d%8dn", i, newsorti);fprintf(fp2, "n*newsort*n");/ hu 2013-1-8void initial()/* 本函数进行初始化工作*/int i, k1;dataio();/ 输入原始数据oldtonew();/ 转化为新号if (bnsopton = 0)/ 节点不优化，新节点号即为老节点号。for (i = 1; i <= n; i +)oldi = i;nnewi = i;elsebnsopt(); / 节点优化mpj = nnewmpj;/mpj: 平衡节点zlsort(nnew); /

29、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; /所有节点的电压巾I值初值都为1.000(v00),电压相角初值都为/ 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;

30、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 节点，即为系统平均电压值

31、。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&q

32、uot;);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, &

33、quot;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",

34、 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&quo

35、t;);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,s

36、ort 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

37、", &pgi, &qgi, &pli, &qli, &v0i);for (i = 1; i <= mdk; i +)fscanf(fp1, "%d %lf", &idki, &dkki);fclose(fp1);void pqflow()/* PQ 分解法计算潮流*/int kq0, iu1N2, nfd1NS, iu2N2, nfd2NS;int i, t;double u1N2, u2N2, di1NS, di2NS;yy1();yzb(0, iu1, u1, di1, nfd1); /form

38、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;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)

39、;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;elsevi = vi - bi;elsekgkq0 = 0;if(kq0 = 0)kq0 = 1;elsekq0 = 0;t +;if(t > it1)break;while(fabs(dsm) > eps) | (kgkq0 != 0);f

40、printf(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)/* 本函数将

41、旧号换成新号*/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因子表上三角矩阵非零非对角元素的数值/

42、 参数4因子表上三角矩阵对角元素/ 参数5因子表上三角各行非零元素个数int i, j, k, i1, i2;编辑版 wordint jj, jj1, jj2, im, x, fdNS;double ai, bNS;nfd1 = 1;for (i = 1; i <= NS; i +)ui = 0.0;for (i = 1; i <= n; i +)/hu 2013-1-9/hu 2013-1-9/nobt 存放的是节点类型, 0: pq 节点 , -1: pv 节点。if (t != 1) | (nobti != -1) && i != mpj) / <-|/

43、for (j = i + 1; j <= n; j +)bj = 0.0;bi = biii;if (kk2 = 0) && (t = 1) && (nobti != -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

44、;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;/存在 (t = 1) 的情况，不多余。break;x = nfdi;dii = 1.0 / bi;ai = dii;k = 0;i1 = i + 1;for (j = i1; j <= n; j +)if (fabs(bj) &g

45、t; 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, "%9.5f%7d", ui, iui);/fprintf(fp2, "n*xxxx*”);/ fprintf(fp2, "%8d%10.5f%8d"

46、,x,u38,fd22);fprintf(fp2, "n*fd*");printi(fd, n);fprintf(fp2, "n*nfd*”);printi(nfd, n);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"

47、;, 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;11 = nfdi;12 = nfdi + 1;for (j = i1; j < i2; j +)k = iuj;bk = bk - bi * uj;bi = bi * dii;f

48、or (i = n; i >= 1; i -) / 回代过程。bi = bi;11 = nfdi;12 = 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

49、);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);

50、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);/变压器节点号由正变负，01d口前zls个为左节点号，后zls个为右节点号。for (i = 1; i <= zls2; i +)/ 冒泡法排序。k1 = i + 1;f0r (j = k1; j <= z1s2; j +)if (01di &