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

一、作业要求

编写程序计算图1所示算例系统的潮流及三相短路电流。 潮流计算：方法不限，计算系统的节点电压和相角。

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

二、电路图及参数

7289354f(3)61 图1 3机9节点系统 表1 9节点系统支路参数

支路 1~4 2~7 3~9 4~5 4~6 5~7 6~9 7~8 8~9 R（.） 0 0 0 X（.） B/2(TK) （.）

表2 9节点系统发电机参数

发电机编号 1 2 3 节点类型 平衡 PV PV PG（.）

VG（.） （.） xd E（.） 表3 9节点系统负荷参数

节点编号 4 5 6 7 8 9 节点类型 PQ PQ PQ PQ PQ PQ

三、计算步骤

(1) 进行系统正常运行状态的潮流计算，求得Ui(0) (2) 形成不含发电机和负荷的节点导纳矩阵YN；

Pi（.） 0 0 1 0 Qi（.） 0 0 0 )和导纳yi(1/jxdi)的并联组合；节点负荷用恒(3) 将发电机表示为电流源Ii(Ei/jxdi阻抗的接地支路表示；形成包括所有发电机支路和负荷支路的节点导纳矩阵Y，即在YN中的发电机节点和负荷节点的自导纳上分别增加发电机导纳yi和负荷导纳yLD,i（yLD,iSLD,iPLDijQLDi）； Vi2Vi21*(4) 利用ZY，计算节点阻抗矩阵，从而得到阻抗矩阵中的第f列；

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

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

(7) 利用公式(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 +

电压幅值：

电压相角：

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 + 近似计算结果： 短路电流： 模值： 相角：

节点电压模值：

五、程序流程图

开 始 输入电网及节点注入数据 计算节点导纳矩阵参数 设置节点电压初值x[0] 设置k←0及最大迭代次数Kmax 计算f（x）及雅克比矩[k]Max|fi（x）|Kmax? Y

六、程序及输入文件

文件： Sheet1 9节点系统支路参数 首节点编号i 1 2 3 4 4 5 6 7 8 Sheet2

9节点系统发电机及负荷参数 节点编号 1 2 3 4 5 6 7 8 9 节点类型 平衡 PV PV PQ PQ PQ PQ PQ PQ P（.） 0 0 -1 0 Q（.） 0 0 0 U（.） Xd（.） E（.） 尾节点编号j 4 7 9 5 6 7 9 8 9 R（.） 0 0 0 X（.） B/2（.） 0 0 0 文件： 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);

Z=Y^(-1);

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) ' ' 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);

num2str(I)]; for i=1:n

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

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