电力系统潮流及短路电流计算程序

班级：

姓名：

学号：

作业要求

编写程序计算图1所示算例系统的潮流及三相短路电流。

潮流计算：方法不限，计算系统的节点电压和相角。

短路电流：4号母线发生金属性三相短路时（z f=0），分别按照精确算法和近似算法计算短路电流、系统中各节点电压以及网络中各支路的电流分布，并对两种情况下的计算结果进行比较。

电路图及参数

图1 3机9节点系统

表1 9节点系统支路参数

计算步骤

进行系统正常运行状态的潮流计算，求得(0)i U &

形成不含发电机和负荷的节点导纳矩阵Y N ；

将发电机表示为电流源i

I &(/i di E jx ''=&)和导纳i y (1/di jx '=)的并联组合；节点负荷用恒阻抗的接地支路表示；形成包括所有发电机支路和负荷支路的节点导纳矩阵Y ，即在Y N 中的发电机节点和负荷节点的自导纳上分别增加发电机导纳i y 和负荷导纳,LD i

y （*

,,22

LD i LDi LDi

LD i

i i S P jQ y V V -==）； 利用1

Z Y

-=，计算节点阻抗矩阵，从而得到阻抗矩阵中的第f 列；

利用公式(6-7)或(6-10)计算短路电流；

利用公式(6-8)或(6-11)计算系统中各节点电压；

利用公式(6-9)计算变压器支路的电流；对输电线路利用П型等值电路计算支路电流。

计算结果

节点导纳矩阵Yn ： Columns 1 through 5

0 0 0 0 + 0

0 0 0 0 0

0 0 0 0 0

0 + 0 0 +

0 0 0 +

0 0 0 + 0

0 0 + 0 0 +

0 0 0 0 0

0 0 0 + 0 0 Columns 6 through 9

0 0 0 0

0 0 + 0 0

0 0 0 0 +

+ 0 0 0

0 + 0 0

0 0 +

0 + 0

0 + +

+ 0 +

电压幅值：

电压相角：

节点有功：

节点无功：

修正后的节点导纳矩阵Y：

Columns 1 through 5

0 0 0 0 + 0

0 0 0 0 0

0 0 0 0 0

0 + 0 0 +

0 0 0 +

0 0 0 + 0

0 0 + 0 0 +

0 0 0 0 0

0 0 0 + 0 0 Columns 6 through 9

0 0 0 0

0 0 + 0 0

0 0 0 0 +

+ 0 0 0

0 + 0 0

0 0 +

0 + 0

0 + +

+ 0 +

节点阻抗矩阵Z的第4列：

+

+

+

+

+

+

+

+

+

精确计算结果：

短路电流：

模值：

相角：

节点电压模值：

支路电流：

i j Iij

1 4 7 9 5 +

4 6 +

5 7 +

6 9 +

7 8 +

8 9 +

近似计算结果：

短路电流：

模值：

相角：

节点电压模值：五、程序流程图

六、程序及输入文件

文件：

文件：

l=9;%支路数

n=9;%节点数

m=6;%PQ节点数

Yn=zeros(n);%初始化节点导纳矩阵Y DATA1=xlsread('',1);

%计算节点导纳矩阵Y

for k=1:l

i=DATA1(k,1);

j=DATA1(k,2);

R=DATA1(k,3);

X=DATA1(k,4);

B2=DATA1(k,5);

Yn(i,i)=Yn(i,i)+1i*B2+1/(R+1i*X);

Yn(j,j)=Yn(j,j)+1i*B2+1/(R+1i*X);

Yn(i,j)=Yn(i,j)-1/(R+1i*X);

Yn(j,i)=Yn(j,i)-1/(R+1i*X);

end

disp('节点导纳矩阵Yn：');

disp(Yn);

G=real(Yn);

B=imag(Yn);

DATA2=xlsread('',2);

P=zeros(1,n);

Q=zeros(1,n);

U=ones(1,n);

P(2:n)=DATA2(2:n,3);

Q(4:n)=DATA2(4:n,4);

U(1:3)=DATA2(1:3,5);

%设置节点电压初值

e(1)=DATA2(1,5);

e(2:n)=;

f(1:n)=;

%设置迭代次数

t=0;

tmax=10;

while t<=tmax

%计算f(x)

a(1:n)=;

c(1:n)=;

for i=2:n

for j=1:n

a(i)=a(i)+G(i,j)*e(j)-B(i,j)*f(j);

c(i)=c(i)+G(i,j)*f(j)+B(i,j)*e(j);

end

end

for i=2:n

deltaP(i)=P(i)-e(i)*a(i)-f(i)*c(i);

end

for j=4:n

deltaQ(j)=Q(j)-f(j)*a(j)+e(j)*c(j);

end

for k=2:3

deltaU2(k)=U(k)*U(k)-e(k)*e(k)-f(k)*f(k);

end

fx=[deltaP(2:n) deltaQ(4:n) deltaU2(2:3)]';

%计算雅克比矩阵J

for i=2:n

for j=2:n

if i~=j

H(i,j)=-(G(i,j)*e(i)+B(i,j)*f(i));

N(i,j)=B(i,j)*e(i)-G(i,j)*f(i);

else

H(i,j)=-a(i)-(G(i,i)*e(i)+B(i,i)*f(i));

N(i,j)=-c(i)+(B(i,i)*e(i)-G(i,i)*f(i));

end

end

end

for i=4:n

for j=2:n

if i~=j

M(i,j)=B(i,j)*e(i)-G(i,j)*f(i);

L(i,j)=G(i,j)*e(i)+B(i,j)*f(i);

else

M(i,j)=c(i)+(B(i,i)*e(i)-G(i,i)*f(i));

L(i,j)=-a(i)+(G(i,i)*e(i)+B(i,i)*f(i));

end

end

end

for i=2:3

for j=2:n

if i~=j

R(i,j)=0;

S(i,j)=0;

else

R(i,j)=-2*e(i);

S(i,j)=-2*f(i);

end

end

end

J=[H(2:n,2:n) N(2:n,2:n);M(4:n,2:n) L(4:n,2:n);R(2:3,2:n) S(2:3,2:n)];

if max(abs(fx))<

%输出结果

break;

else

%求解修正方程获得dx

dx=-J^(-1)*fx;

dx=dx';

e(2:n)=e(2:n)+dx(1:n-1);

f(2:n)=f(2:n)+dx(n:2*(n-1));

t=t+1;

end

end

if t>tmax

str='潮流计算不收敛';

disp(str);

else

a(1:n)=;

c(1:n)=;

for i=1:n

for j=1:n

a(i)=a(i)+G(i,j)*e(j)-B(i,j)*f(j);

c(i)=c(i)+G(i,j)*f(j)+B(i,j)*e(j);

end

end

for i=1:n

U(i)=e(i)+1i*f(i);

amp(i)=abs(U(i));

arg(i)=angle(U(i));

P(i)=e(i)*a(i)+f(i)*c(i);

Q(i)=f(i)*a(i)-e(i)*c(i);

end

disp('电压幅值：');

disp(amp);

disp('电压相角：');

disp(arg);

disp('节点有功：');

disp(P);

disp('节点无功：');

disp(Q);

end

%计算短路电流

f=4;

zf=;

%修正节点导纳矩阵

Xd=DATA2(1:3,6);

E=DATA2(1:3,7);

for i=1:3

Ii(i)=E(i)/(1i*Xd(i));

end

Y=Yn;

for i=1:3

Y(i,i)=Y(i,i)+1/(1i*Xd(i));

end

for j=4:n

Y(j,j)=Y(j,j)+(-P(j)+1i*Q(j))/(U(j)*U(j)); end

disp('修正后的节点导纳矩阵Y：');

disp(Y);

disp('节点阻抗矩阵Z的第4列：');

disp(Z(:,4));

%精确计算

disp('精确计算结果：');

U0=U;

If=U0(f)/(Z(f,f)+zf);

amp=abs(If);

arg=atand(imag(If)/real(If));

disp('短路电流：');

disp('模值：');

disp(amp);

disp('相角：');

disp(arg);

for i=1:n

U(i)=U0(i)-Z(i,f)*If;

amp=abs(U);

end

disp('节点电压模值：');

disp(amp);

disp('支路电流：');

str=['i ''j '' Iij'];

disp(str);

for k=1:l

i=DATA1(k,1);

j=DATA1(k,2);

r=DATA1(k,3);

x=DATA1(k,4);

z=r+1i*x;

I=(U(i)-U(j))/z;

str=[num2str(i) ' ' num2str(j) ' ' num2str(I)];

disp(str);

end

%近似计算

disp('近似计算结果：');

U0(1:n)=;

If=U0(f)/(Z(f,f)+zf);

amp=abs(If);

arg=atand(imag(If)/real(If));

disp('短路电流：');

disp('模值：');

disp(amp);

disp('相角：');

disp(arg);

U(i)=U0(i)-Z(i,f)*If;

amp=abs(U);

end

disp('节点电压模值：'); disp(amp);