多机系统短路故障后时域仿真.zip

function [f1,f2,o,pt0]=calculate(gen,tranxt,line,VD0,SLD0,S20,S30) global E1 E2 E3 Z11 Z12 Z13 a11 a12 a13 Z22 Z21 Z23 a22 a21 a23 Z33 Z31 Z32 a33 a31 a32 global z11 z12 z13 a_11 a_12 a_13 z22 z21 z23 a_22 a_21 a_23 z33 z31 z32 a_33 a_31 a_32 %% 运行初态计算 j=sqrt(-1); X1=gen(1,2)+tranxt(1); X1B=X1+line(1,1); X2=gen(2,2)+tranxt(2); XBD=tranxt(4)+line(2,1); X3=gen(3,2)+tranxt(3); YLD=conj(SLD0)/(VD0^2); ZLD=1/YLD; E3=((VD0+imag(S30)*X3/VD0)^2+(real(S30)*X3/VD0)^2)^0.5; %E3=(()^0.5+()^0.5)^0.5 o_3=atan2d(real(S30)*X3/VD0,VD0+imag(S30)*X3/VD0); %3号发电机功角初始值 SBD2=SLD0-S30; %SBD'' VB=((VD0+imag(SBD2)*XBD/VD0)^2+(real(SBD2)*XBD/VD0)^2)^0.5; oB=atan2d(real(SBD2)*XBD/VD0,VD0+imag(SBD2)*XBD/VD0); SBD1=SBD2+j*XBD*((real(SBD2))^2+(imag(SBD2))^2)/VD0^2; S1=SBD1-S20; pt1=real(S1); %1号发电机有功功率初始值 E1=((VB+imag(S1)*X1B/VB)^2+(real(S1)*X1B/VB)^2)^0.5; % 利用S1,UB,X1B可以计算E1 o1B=atan2d(real(S1)*X1B/VB,VB+imag(S1)*X1B/VB); o_1=oB+o1B; %1号发电机功角初始值 E2=((VB+imag(S20)*X2/VB)^2+(real(S20)*X2/VB)^2)^0.5; % 利用S2,UB,X2可以计算E2 o2B=atan2d(real(S20)*X2/VB,VB+imag(S20)*X2/VB); o_2=oB+o2B; %2号发电机功角初始值 pt2=real(S20); pt3=real(S30); o=[o_1,o_2,o_3]; %发电机功角初始值 pt0=[pt1,pt2,pt3]; %发电机有功功率初始值 %% 故障初态计算 % (a)短路时各短路点负荷输入阻抗的计算 XLD2=0.35*abs(ZLD); XD=(gen(3,3)+tranxt(3))*XLD2/(gen(3,3)+tranxt(3)+XLD2); XB=(XBD+XD)*(gen(2,3)+tranxt(2))/(XBD+XD+gen(2,3)+tranxt(2)); Xff2=(line(1,1)+XB)*(gen(1,3)+tranxt(1))/(line(1,1)+XB+gen(1,3)+tranxt(1)); % 在零序网络中，短路点零序输入阻抗 XB0=(line(2,2)+tranxt(4)+tranxt(3))*tranxt(2)/(line(2,2)+tranxt(4)+tranxt(3)+tranxt(2)); Xff0=(line(1,2)+XB0)*tranxt(1)/(line(1,2)+XB0+tranxt(1)); X0=Xff2*Xff0/(Xff2+Xff0); % (b)短路时系统输入阻抗和转移阻抗的计算 e3=0;i3=1; VD=j*X3*i3; ID=YLD*VD; IBD=i3+ID; vb=VD+j*XBD*IBD; e2=0; I2=(e2-vb)/(j*X2); IAB=-I2+IBD; va=vb+j*line(1,1)*IAB; IA=va/(j*X0); I1=IA+IAB; e1=va+j*X1*I1; Z31=e1/i3; a31=90-(atan2d(imag(Z31),real(Z31))); Z13=Z31;a13=a31; Z21=e1/(-I2); a21=90-(atan2d(imag(Z21),real(Z21))); Z12=Z21;a12=a21; Z11=e1/I1; a11=90-(atan2d(imag(Z11),real(Z11))); % 令E1=0,E2~=0,则从B点左边部分网络的总阻抗为 Xb=line(1,1)+X1*X0/(X1+X0); iAB=-vb/(j*Xb); i2=-iAB+IBD; e2=vb+j*X2*i2; Z32=e2/i3; a32=90-(atan2d(imag(Z32),real(Z32))); Z23=Z32;a23=a32; Z22=e2/i2; a22=90-(atan2d(imag(Z22),real(Z22))); % 现在通过网络化简求取输入阻抗Z33，节点D左边部分网络总阻抗为 Xd=XBD+X2*Xb/(X2+Xb); Z33=j*X3+(j*Xd)*ZLD/(j*Xd+ZLD); a33=90-(atan2d(imag(Z33),real(Z33))); %% 故障线路切除后系统状态计算 e3=0;e2=0;i3=1; I1=IAB; e1=vb+j*(X1+2*line(1,1))*I1; z31=e1/i3; a_31=90-(atan2d(imag(z31),real(z31))); z13=z31;a_13=a_31; z21=e1/(-I2); a_21=90-(atan2d(imag(z21),real(z21))); z12=z21;a_12=a_21; z11=e1/I1; a_11=90-(atan2d(imag(z11),real(z11))); % 令E1=0,E2~=0,以计算z32和z22 I1=-vb/(j*(X1+2*line(1,1))); I2=IBD-I1; e2=vb+j*X2*I2; z32=e2/i3; a_32=90-(atan2d(imag(z32),real(z32))); z23=z32;a_23=a_32; z22=e2/I2; a_22=90-(atan2d(imag(z22),real(z22))); % 现在通过网络化简求取输入阻抗z33，节点D左边部分网络总阻抗为 Xd=XBD+X2*(X1+2*line(1,1))/(X1+X2+2*line(1,1)); z33=j*X3+(j*Xd)*ZLD/(j*Xd+ZLD); a_33=90-(atan2d(imag(z33),real(z33))); f1.p1=@P1; f1.p2=@P2; f1.p3=@P3; f2.p1=@P_1; f2.p2=@P_2; f2.p3=@P_3; f1=struct2cell(f1); f2=struct2cell(f2); %把结构体数组转换成元胞数组 end %% 故障期间的功率特性分析 function p1=P1(o12,o13) global E1 E2 E3 Z11 Z12 Z13 a11 a12 a13 p1=E1^2*sind(a11)/(abs(Z11))+E1*E2*sind(o12-a12)/(abs(Z12))+E1*E3*sind(o13-a13)/(abs(Z13)); end function p2=P2(o21,o23) global E1 E2 E3 Z22 Z21 Z23 a22 a21 a23 p2=E2^2*sind(a22)/(abs(Z22))+E2*E1*sind(o21-a21)/(abs(Z21))+E2*E3*sind(o23-a23)/(abs(Z23)); end function p3=P3(o31,o32) global E1 E2 E3 Z33 Z31 Z32 a33 a31 a32 p3=E3^2*sind(a33)/(abs(Z33))+E3*E1*sind(o31-a31)/(abs(Z31))+E3*E2*sind(o32-a32)/(abs(Z32)); end %% 故障切除后的功率特性分析 function p_1=P_1(o12,o13) global E1 E2 E3 z11 z12 z13 a_11 a_12 a_13 p_1=E1^2*sind(a_11)/(abs(z11))+E1*E2*sind(o12-a_12)/(abs(z12))+E1*E3*sind(o13-a_13)/(abs(z13)); end function p_2=P_2(o21,o23) global E1 E2 E3 z22 z21 z23 a_22 a_21 a_23 p_2=E2^2*sind(a_22)/(abs(z22))+E2*E1*sind(o21-a_21)/(abs(z21))+E2*E3*sind(o23-a_23)/(abs(z23)); end function p_3=P_3(o31,o32) global E1 E2 E3 z33 z31 z32 a_33 a_31 a_32 p_3=E3^2*sind(a_33)/(abs(z33))+E3*E1*sind(o31-a_31)/(abs(z31))+E3*E2*sind(o32-a_32)/(abs(z32)); end