电力网潮流电压计算例题与matlab程序

电力网潮流电压计算例题与MATLAB程序

编写佘名寰

本文介绍了两个电力网潮流计算例题。一例为5个节点5条支路；另一例有6个节点7条支路，有PQ负荷节点也有PV发电机节点，变压器支路标么变比不为1。本文给出了完整的计算程序和计算结果。程序包括计算网络节点导纳矩阵和潮流电压两部分，程序中引用的M函数见参考文献1。例题选自研究生教材，比较典型，计算结果与教材核对基本一致。本文可供电力专业师生和基层技术人员潮流计算时参考。

【例1.3】如图1-4所示，五个节点三条母线简单系统，系统参数均用标么值（pu.）表示。节点4、5通过变压器与节点2、3相连，变比为1.05:1，节点1、2、3为PQ负荷节点，节点4为PV发电机节点，节点5为参考节点，节点4、5电压均为1.05，图中箭头表示各支路的正方向，各个节点功率如表2.1所示,计算各节点电压、发电机无功功率和系统潮流分布：

图1-4[例1.3]网络接线图

【例1.3】程序说明：

本例题与参考文献2的【例2.3.1】不同之处在于有两个变比不为1的变压器支路，并且引入了PV发电机节点。由于网络结构与例图1-1完全相同，网络节点导纳矩阵已在文献1之1.3节求出不再赘述。潮流计算程序参考【例2.3.1】

节点电压符号矩阵：一共五个节点，独立节点数N1=4，PV节点数N2=1

节点电压幅值符号矩阵为 u=sym('[u1,u2,u3,u4,u5]');

节点电压相角符号矩阵为 delt=sym('[d1,d2,d3,d4,d5]');

节点功率给定值为：三个PQ节点，一个PV 节点

p=[-1.6,-2,-3.7,5];q=[-0.8,-1.0,-1.3];

【例1.3】源程序NU114N.m

%******* NU114N.m , example 2-3-2, fig 2-2 ********

% %The following Program for load flow calculation is based on MATLAB2007 clear

% bus 1,2,3 is PQ bus,bus 4 is PV bus,bus 5 is slack bus

global Np Nb bt p m k

% Np is number of node point,Nb is number of braches,

% p the transformer turns ratio,with off-nominal tap-setting,fig

2-4 ,PSCLF p6,

% bt is two-dimensional array ,line one is point 1,line two is point 2, % m is number of node point, k is number of transformers

Np=5;Nb=5;

bt=[2,3;4,5];p=[1/1.05,1/1.05];m=5;k=2;

% nstart--the start point of branches ,nend -- the end point,

% mm -- network incidence matrix

nstart=[2,2,1,4,5];

nend=[1,3,3,2,3];

mm=ffm(nstart,nend);

% zb1,the series impedances of transmission line

% yb1,the series admittances of transmission line

zb1r=[0.04,0.08,0.1,0.0,0.0];

zb1i=[0.25,0.30,0.35,0.015,0.03];

zb1=zb1r+zb1i*j;

yb1=zb1.^(-1);

yb=diag(conj(yb1'));

y=mm*yb*(mm)';

% yb0 ,the shunt admittances of transmission line

yb0i=[0.25,0.25,0.0,0.0,0.0];

yb0=0+yb0i*j;

y0=diag(conj(yb0'));

yg=mm*y0*(mm)';

yn=diag(diag(yg));

yy=yn-yg;

y=y+yn;

% yn0,the node-admittance of shunt capacitor

%yn0i=[];

%yn0=0+yn0i*j;

%ync=diag(conj(yn0'));

%y=y+ync;

[Y,YY]=fdt1(y,yy);

G=real(Y);B=imag(Y);

% YY is the shund admittances of the line and transformers at each end, % Y is node-admittance matrix for network

%******************************************************************** *

u=sym('[u1,u2,u3,u4,u5]');delt=sym('[d1,d2,d3,d4,d5]');

% u -- node voltage magnitude,delt -- angle

p=[-1.6,-2,-3.7,5];q=[-0.8,-1.0,-1.3];

k=0;precision=1;

N1=4;%the N1 is the amount of the PQ and PV bus

N2=1;%the N2 is the amount of the PV bus

for m=1:N1

for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt( n)));

end

pp(m)=p(m)-sum(pt);

end

for m=1:N1-N2

for n=1:N1+1

qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt( n)));

end

qq(m)=q(m)-sum(qt);

end

J1=jacobian(pp,['d1,d2,d3,d4,u1,u2,u3']);

J2=jacobian(qq,['d1,d2,d3,d4,u1,u2,u3']);

J=vertcat(J1,J2);

uu=[1.0,1.0,1.0,1.05,1.05];dd=[0,0,0,0,0];

while precision>0.00001

u1=uu(1);u2=uu(2);u3=uu(3);u4=uu(4);u5=uu(5);

d1=dd(1);d2=dd(2);d3=dd(3);d4=dd(4);d5=dd(5);

for m=1:N1

PP(m)=eval(pp(m));

end

for m=1:N1-N2

PP(N1+m)=eval(qq(m));

end

JJ=eval(J);

du=-inv(JJ)*PP';precision=max(abs(du));

for n=1:N1

dd(n)=dd(n)+du(n);

end

for n=1:N1-N2

uu(n)=uu(n)+du(N1+n);

end

k=k+1;

end

k-1,dd*180/pi,uu

%***********************************************

% the following program is used to calculate the Sm and Smn

for n=1:N1+1

U(n)=uu(n)*(cos(dd(n))+j*sin(dd(n)));

end

Um=conj(U');

I=Y*Um;

Sm=diag(Um)*conj(I)

for m=1:N1+1

for n=1:N1+1

Smn(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n))+U(m)*conj(U(m)*YY (m,n));

end

end

Smn

【例1.3】潮流计算结果：

迭代次数

ans =

4

节点1至5电压相角（度）

ans =

-4.7785 17.8535 -4.2819 21.8433 0

节点1至5电压幅值（pu）

uu =

0.8622 1.0779 1.0364 1.0500 1.0500

节点1至5功率（pu）

Sm =

-1.6000 - 0.8000i

-2.0000 - 1.0000i

-3.7000 - 1.3000i

5.0000 + 1.8131i

2.5794 + 2.2994i

支路功率（行列编号均为节点1至5）

Smn =

0 -1.4662 - 0.4091i -0.1338 - 0.3909i 0 0

1.5845 + 0.6726i 0 1.4155 - 0.2443i -5.0000 - 1.4282i 0

0.1568 + 0.4713i -1.2774 + 0.2032i 0 0 -2.5794 - 1.9745i

0 5.0000 + 1.8131i 0 0 0

0 0 2.5794 + 2.2994i 0 0

【例1.4】如图1-5所示，六个节点网络。节点2、1间和节点4、3间通过变压器相连，变比分别为0.909:1和0.976:1，节点1、2、3、4为PQ负荷节点，节点5为PV发电机节点，节点电压为1.10;节点6为参考节点,节点电压为1.05。输电线路和变压器数据见表2.2,各个节点功率如表2.3所示,计算各节点电压、发电机无功功率和系统潮流分布：

图1-5[例1.4]网络接线图

表1-2 输电线路和变压器数据汇编

线路首端节点编号线路末端

节点编号

充电功率换

算导纳BC/2

线路阻抗换

算电导G

线路阻抗换

算电纳B

变压器变比

6 2 0.0070 0.2791 -1.2910

6 4 0.0090 0.2170 -0.9137

5 1 0.0 0.2224 -0.3230

5 3 0.0 0.2383 -0.6542

2 1 0.0 0.0 -3.7594 0.909 2 4 0.0076 0.2771 -1.1625

4 3 0.0 0.0 -1.6667 0.976

表1-3 各节点功率

【例1.4】程序说明：

本例题计算程序可分为形成节点导纳矩阵和计算节点电压潮流两部分。

计算节点导纳矩阵参考文献1第1.3节。这里有6个节点7条支路，节点编号原则是先PQ节点，后PV节点，平衡节点编号最大，支路序号如图所示。各支路的首末端节点编号用数组nstart和nend表示，各支路导纳用数组yb1r,yb1i表示，各支路两端充电电容导纳为数组yb0i:

Np=6;Nb=7;

nstart=[6,6,2,3,5,5,1];nend=[4,2,4,4,3,1,2];

yb1r=[0.2170,0.2791,0.2771,0.0,0.2383,0.2224,0.0];

yb1i=[-0.9137,-1.2910,-1.1625,-1.6667,-0.6542,-0.3230,-3.7594];

yb0i=[0.0099,0.0070,0.0076,0.0,0.0,0.0,0.0];

两台变压器跨接在节点2-1和4-3之间，参见图1-2，i节点对应2、4，j节点对应1、3，变比为1：1/0.909和1:1/0.976,程序中用数组bt[ ]和p[ ]表示，m=6,k=2表明用6个节点两条变压器支路：

bt=[2,4;1,3];p=[1/0.909,1/0.976];m=6;k=2;

计算潮流电压参照例题[1.3]。这里独立节点数N1=5，PV节点数N2=1，节点电压初始给定值幅值为数组uu,相角为dd，节点有功功率为p,无功功率为q:

uu=[1.000,1.000,1.000,1.000,1.100,1.050];dd=[0,0,0,0,0,0];

p=[-.275,0,-0.150,-0.250,0.250];q=[-0.065,0,-0.090,-0.02

【例1.4】源程序NU117N.m

%*******NU117N.M example 2.3.3, fig 2-3 ************************** % %The following Program for load flow calculation is based on MATLAB R2007 clear

% bus 1,2,3,4 is PQ bus,bus 5 is PV bus,bus 6 is slack bus

u=sym('[u1,u2,u3,u4,u5,u6]');delt=sym('[d1,d2,d3,d4,d5,d6]');

global Np Nb bt p m k

Np=6;Nb=7;

bt=[2,4;1,3];p=[1/0.909,1/0.976];m=6;k=2;

nstart=[6,6,2,3,5,5,1];nend=[4,2,4,4,3,1,2];

mm=ffm(nstart,nend);

% yb1,the series admittance of transmission line

yb1r=[0.2170,0.2791,0.2771,0.0,0.2383,0.2224,0.0];

yb1i=[-0.9137,-1.2910,-1.1625,-1.6667,-0.6542,-0.3230,-3.7594];

yb1=yb1r+yb1i*j;

%yb1=zb1.^(-1);

yb=diag(conj(yb1'));

y=mm*yb*(mm)';

% yb0 ,the shunt admittances of transmission line

yb0i=[0.0099,0.0070,0.0076,0.0,0.0,0.0,0.0];

yb0=0+yb0i*j;

y0=diag(conj(yb0'));

yg=mm*y0*(mm)';

yn=diag(diag(yg));

yy=yn-yg;

y=y+yn;

[Y,YY]=fdt1(y,yy);

%******************************************************************** ******

G=real(Y);B=imag(Y);

p=[-.275,0,-0.150,-0.250,0.250];q=[-0.065,0,-0.090,-0.025];

k=0;precision=1;

N1=5;%the N1 is the amount of the PQ and PV bus

N2=1;%the N2 is the amount of the PV bus

for m=1:N1

for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt( n)));

end

pp(m)=p(m)-sum(pt);

end

for m=1:N1-N2

for n=1:N1+1

qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt( n)));

end

qq(m)=q(m)-sum(qt);

end

J1=jacobian(pp,['d1,d2,d3,d4,d5,u1,u2,u3,u4']);

J2=jacobian(qq,['d1,d2,d3,d4,d5,u1,u2,u3,u4']);

J=vertcat(J1,J2);

uu=[1.000,1.000,1.000,1.000,1.100,1.050];dd=[0,0,0,0,0,0];

while precision>0.0000001

u1=uu(1);u2=uu(2);u3=uu(3);u4=uu(4);u5=uu(5);u6=uu(6);

d1=dd(1);d2=dd(2);d3=dd(3);d4=dd(4);d5=dd(5);d6=dd(6);

for m=1:N1

PP(m)=eval(pp(m));

end

for m=1:N1-N2

PP(N1+m)=eval(qq(m));

end

JJ=eval(J);

du=-inv(JJ)*PP';precision=max(abs(PP));

for n=1:N1

dd(n)=dd(n)+du(n);

end

for n=1:N1-N2

uu(n)=uu(n)+du(N1+n);

end

k=k+1;

end

k-1,dd*180/pi,uu

%***********************************************

% the following program is used to calculate the Sm and Smn

for n=1:N1+1

U(n)=uu(n)*(cos(dd(n))+j*sin(dd(n)));

end

Um=conj(U');

I=Y*Um;

Sm=diag(Um)*conj(I)

for m=1:N1+1

for n=1:N1+1

Smn(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n))+U(m)*conj(U(m)*YY (m,n));

end

end

Smn

【例1.4】程序计算结果：

ans =

4

节点电压相角（度）

ans =

-12.6740 -9.7623 -12.3814 -12.2420 -2.7721 0

节点电压幅值（pu）

uu =

1.0001 0.9296 0.9220 0.9204 1.1000 1.0500

节点功率（pu）

Sm =

-0.2750 - 0.0650i

0.0000 - 0.0000i

-0.1500 - 0.0900i

-0.2500 - 0.0250i

0.2500 + 0.0958i

0.4744 + 0.2166i

支路功率矩阵（6行6列均对应节点编号）

Smn =

0 -0.1953 - 0.0801i 0 0 -0.0797 + 0.0151i 0

0.1953 + 0.0919i 0 0 0.0456 - 0.0059i 0 -0.2410 - 0.0861i 0 0 0 -0.0035 - 0.0323i -0.1465 - 0.0577i 0

0 -0.0452 - 0.0052i 0.0035 + 0.0330i 0 0 -0.2084 - 0.0529i 0.0892 - 0.0013i 0 0.1608 + 0.0971i 0 0 0

0 0.2529 + 0.1275i 0 0.2215 + 0.0891i 0 0

参考文献：

1电力网节点导纳矩阵计算例题与程序著者佘名寰2016年百度文库

2牛顿-拉夫逊迭代法电力网潮流计算方法与程序著者佘名寰2016年百度文库