电力系统计算报告.doc

电力系记录算报告

院 （系） 电气工程及自动化
授 课 教 师 胡林献
姓 名 张远实
学 号 13S106052

P-Q分解潮流法
简述P-Q分解潮流法基本原理、计算过程、计算框图。用C语言编制P-Q分解法潮流程序，并用电科院6机22节点系统加以验证。规定采用稀疏技术、因子表技术和节点优化技术，并考虑负荷静态特性。
P-Q分解潮流法基本原理：
P-Q分解法即是基于采用极坐标形式表达牛顿法，其依照电力系统实际运营状态线路参数R/X普通很小状况，对求解修正量修正方程系数矩阵加以简化，使其变为常数阵(即所谓等斜率)，且 P、Q迭代解耦。这样可减少每次迭代计算时间，提高计算速度，又不影响最后成果，因而是普通选用一种办法。但在低电压配电网中，当线路 R/X比值很大时，也许浮现不收敛状况，此时应考虑更换其他办法。
计算过程：
1、 形成有功迭代和无功迭代简化雅克比矩阵 和
2、 给定PQ节点初值和各节点电压相角初值
3、 作有功迭代，计算，解修正方程，得各节点电压相角修正值。
4、 作无功迭代，计算，解修正方程，得各节点电压幅值修正量。
5、 返回第三步，继续迭代到满足规定为止。
计算框图流程：

Y
Y
输入原始数据
形成矩阵 及 并进行三角分解
设PQ节点电压初值，各节点电压相角初值
制迭代计数k=0

用公式计算不平衡功率，计算

置
解修正方程求

置

用公式计算不平衡功率，计算

置

解修正方程求

置

计算平衡节点功率及所有线路功率
输出
Y
Y

算例描述
用电科院6机22节点算例进行验证，详细参数和网络拓扑图如下所示

表1 线路参数
支路号
首末端节点号
支路电阻
支路电抗
对地电纳/2
1
7-8
0.0106
0.0740
0.0
2
7-9
0.0147
0.1040
0.0
3
8-9
0.0034
0.0131
0.0
4
9-22
0.0559
0.2180
0.1954
5
12-13
0.00245
0.0255
1.395
6
14-19
0.0034
0.0200
0.0
7
16-19
0.0578
0.2180
0.1887
8
16-20
0.0163
0.0662
0.2353
9
16-21
0.0374
0.1780
0.164
10
16-18
0.0033
0.0333
0.0
11
19-21
0.0114
0.0370
0.0
12
20-22
0.0214
0.0859
0.3008
13
21-22
0.0150
0.0607
0.2198
14
8-22
0.0537
0.1900
0.1653
15
11-12
0.0033
0.0343
1.08797

表2 变压器支路数据
支路号
首末端节点号
电阻
电抗
变比
1
7-1
0.0
0.0150
1.050
2
9-2
0.0
0.0217
1.075
3
22-3
0.0
0.0124
1.100
4
19-4
0.0
0.0640
1.025
5
18-5
0.0
0.0375
1.050
6
17-6
0.0
0.0337
1.000
7
10-9
0.0
-0.002
1.000
8
11-10
0.0
0.0180
1.000
9
15-12
0.0
0.0180
1.000
10
17-13
0.0
0.0100
1.000
11
15-14
0.0
-0.002
1.000
12
16-17
0.0
0.0010
1.027

表3 并联电容器数据
节点号
电纳
12
-1.3665
16
0.50176
11
-1.3665
12
-1.3665

表4 母线功率数据
母线号
发电机输出有功
发电机输出无功
负荷有功
负荷无功
1
5.9631
1.7355
0.00
0.00
2
6
3.2
0.00
0.00
3
3.1
4.6
0.00
0.00
4
1.6
0.7
0.00
0.00
5
4.3
3.34
0.00
0.00
6
-0.01
1.0
0.00
0.00
7
0.00
0.00
0.00
0.00
8
0.00
0.00
2.87
1.44

9
0.00
0.00
3.76
2.21
10
0.00
0.00
0.00
0.00
11
0.00
0.00
0.00
0.00
12
0.00
0.00
0.00
0.00
13
0.00
0.00
0.00
0.00
14
0.00
0.00
0.00
0.00
15
0.00
0.00
0.00
0.00
16
0.00
0.00
5.0
2.9
17
0.00
0.00
0.00
0.00
18
0.00
0.00
4.3
2.6
19
0.00
0.00
0.864
0662
20
0.00
0.00
0.719
0.474
21
0.00
0.00
0.7
0.5
22
0.00
0.00
2.265
1.69

表5 无功可调母线数据
母线号
电压幅值（标幺值）
无功下限值
无功上限制
1
1.0
-5
10
3
1.0
-5
5
6
1.0
-5
6

表6 发电机参数
母线号
暂态电抗
转子惯性时间常数
阻尼系数
发电机有功出力下限
发电机有功出力上限
1
0.0150
140.82
0.00
3
16.5
2
0.0382
30.00
0.00
1.2
6.6
3
0.0396
79.50
0.00
1.5
8.25
4
0.1210
15.68
0.00
0.4
2.2
5
0.0480
39.20
0.00
1.02
5.61
6
0.1976
2.62
0.00
0.2
1.1

图1 电科院6机22节点系统图
计算成果
I V CA PL QL PG QG
1 1.00000 0.00000 0.00000 0.00000 5.96312 1.73549
2 0.97384 -11.311893 0.00000 0.00000 600000 3.0
3 1.00000 -27.459057 0.00000 0.00000 3.10000 3.14731
4 1.02190 -25.168610 0.00000 0.00000 1.60000 0.70000
5 1.04392 -28.205656 0.00000 0.00000 4.30000 3.34000
6 1.00000 -37.566566 0.00000 0.00000 -0.01000 0.91664
7 1.02697 -5.247170 0.00000 0.00000 0.00000 0.00000
8 0.96847 19.788813 -2.87000 -1.44000 0.00000 0.00000
9 0.98081 -19.738201 -3.76000 -2.21000 0.00000 0.00000
10 0.97985 -19.304344 0.00000 0.00000 0.00000 0.00000
11 0.99050 -23.173771 0.00000 0.00000 0.00000 0.00000
12 0.99359 -30.499034 0.00000 0.00000 0.00000 0.00000
13 0.98198 -35.528090 0.00000 0.00000 0.00000 0.00000
14 1.00020 -30.695861 0.00000 0.00000 0.00000 0.00000
15 1.00103 -30.720282 0.00000 0.00000 0.00000 0.00000
16 0.99308 -27.750976 -5.00000 -2.90000 0.00000 0.00000
17 0.96911 -37.546641 0.00000 0.00000 0.00000 0.00000
18 0.98360 -37.696751 -4.30000 -2.60000 0.00000 0.00000
19 1.00776 -31.018334 -0.86400 -0.66200 0.00000 0.00000
20 1.01475 -35.674534 -0.71900 -0.47400 0.00000 0.00000
21 1.101546 -32.104404 -0.70000 -0.50000 0.00000 0.00000
22 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.0 -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.38055 -2.04951 -3.59392 -1.18599
12 13 3.38055 -1.10216 -3.35200 -1.32306
12 15 0.21337 -0.40993 -0.21337 041383
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

源程序
#include "math.h"
#include "stdio.h"

#define NS //最大节点数
#define NS2 NS * 2
#define NS4 1000 //NS4、NS必要不不大于2*zls。
#define ZS 3000 //最大支路数
#define ZS2 ZS * 2
#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(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 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 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 kint 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] = %d\n"，i，izl[i]，i，jzl[i]);
}
for (i = 1；i <= nb；i ++)
{
for (j = 1；j <= nj ++)
if (nob[i] == newsort[j])
{
nob[i] = j;
break;
}
printf("nob[%d] = %d\n"，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;
}
}
}
//hu -1-8
printc('-'，58);
fprintf(fp2，"\n\n *******newsort*******\n");
for (i = 1；i <= n；i ++)
{
fprintf(fp2，"%8d%8d\n"，i，newsort[i]);
}
fprintf(fp2，"\n *******newsort*******\n");

// hu -1-8
}

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.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%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，"\n\n*******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,"\n\n******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 file\n");
scanf("%s"，inname);
fp1 = fopen(inname，"r");
printf("\nplease output the name of data file\n");
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分解法计算潮流 ****//

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：%d\t%d\t%f\n"，t，dsd，dsm);
else
printf("Q：%d\t%d\t%f\n"，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] / v00;
else
v[i] = v[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 <= NS；i ++) //hu -1-9
{
u[i] = 0.0;
} //hu -1-9
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)状况，