华北电力大学电力系统暂态上机程序设计完整版.doc

clc;clear; %-输入已知条件-%bus_Num1=3; %节点数bus_Num2=5; %包括发电机节点的节点数branch_Num1=3; %线路数branch_Num2=5; %包括发电机支路的支路数branch1_No1=1,1,2;branch1_No2=2,3,3;branch2_No1=1,1,2,1,2;branch2_No2=2,3,3,4,5; %每条支路始末端编号 %输入支路各序阻抗，z1_branch表示支路正序阻抗，z2_branch表示支路负序阻抗，z0_branch表示支路零序阻抗z1_branch(1)=j*0.1;z1_branch(2)=j*0.1;z1_branch(3)=j*0.1; z2_branch=z1_branch; z0_branch(1)=j*0.2;z0_branch(2)=j*0.2;z0_branch(3)=j*0.2; %-第一步：不考虑发电机节点计算节点导纳矩阵-%节点导纳矩阵，Y1表示不计发电机节点的正序网络节点导纳阵，Y2表示不计发电机节点的负序网络节点导纳阵，Y0表示不计发电机节点的零序网络节点导纳阵，Y1=zeros(bus_Num1);%全为0矩阵（3*3）Y1(1,1)=1/(j*0.15);Y1(2,2)=1/(j*0.075);%请同学们求正序节点导纳矩阵for m=1:branch_Num1Y1(branch1_No1(m),branch1_No1(m)=Y1(branch1_No1(m),branch1_No1(m)+1/z1_branch(m);Y1(branch1_No1(m),branch1_No2(m)=Y1(branch1_No1(m),branch1_No2(m)-1/z1_branch(m);Y1(branch1_No2(m),branch1_No1(m)=Y1(branch1_No2(m),branch1_No1(m)-1/z1_branch(m);Y1(branch1_No2(m),branch1_No2(m)=Y1(branch1_No2(m),branch1_No2(m)+1/z1_branch(m);end Y2=Y1;%负序等于正序 Y0=zeros(bus_Num1);Y0(1,1)=1/(j*0.0525);Y0(2,2)=1/(j*0.0263);%请同学们求零序节点导纳矩阵for m=1:branch_Num1Y0(branch1_No1(m),branch1_No1(m)=Y0(branch1_No1(m),branch1_No1(m)+1/z0_branch(m);Y0(branch1_No1(m),branch1_No2(m)=Y0(branch1_No1(m),branch1_No2(m)-1/z0_branch(m);Y0(branch1_No2(m),branch1_No1(m)=Y0(branch1_No2(m),branch1_No1(m)-1/z0_branch(m);Y0(branch1_No2(m),branch1_No2(m)=Y0(branch1_No2(m),branch1_No2(m)+1/z0_branch(m);end %-第二步：考虑发电机节点计算节点导纳矩阵-%节点导纳矩阵，YY1表示计及发电机节点的正序网络节点导纳阵，YY2表示计及发电机节点的负序网络节点导纳阵z1_branch(4)=j*0.05;z1_branch(5)=j*0.025; z2_branch=z1_branch; YY1=zeros(bus_Num1);YY1(4,4)=1/(j*0.1);YY1(5,5)=1/(j*0.05);%请同学们求正序节点导纳矩阵for m=1:branch_Num2YY1(branch2_No1(m),branch2_No1(m)=YY1(branch2_No1(m),branch2_No1(m)+1/z1_branch(m);YY1(branch2_No1(m),branch2_No2(m)=YY1(branch2_No1(m),branch2_No2(m)-1/z1_branch(m);YY1(branch2_No2(m),branch2_No1(m)=YY1(branch2_No2(m),branch2_No1(m)-1/z1_branch(m);YY1(branch2_No2(m),branch2_No2(m)=YY1(branch2_No2(m),branch2_No2(m)+1/z1_branch(m);end YY2=YY1;%负序等于正序 %-第三步：计算节点阻抗矩阵-%Z1=inv(Y1);Z2=inv(Y2);Z0=inv(Y0);clc;clear; %数据来源于教材电力系统暂态分析P77例（3-4），P143例（5-7） Matrix %计算节点导纳矩阵、节点阻抗矩阵，形成全局变量 Fault_Node=input(输入短路点编号 ;n Fault_Node=); Fault_Type=input(输入短路类型 ;n(1)Fault_Type=0为三相短路;n(2)Fault_Type=1为a相接地短路;n(3)Fault_Type=2为a相经10欧姆电阻接地短路n(4)Fault_Type=3为bc两相相间短路n(5)Fault_Type=4为bc两相短路接地nFault_Type=) a=-0.5+j*sqrt(3)/2;T=1 1 1 % T为对称分量法的变换矩阵，见P87公式（4-4） a2 a 1 a a2 1; %-第一步：计算短路点的序电流，相电流-%计根据故障类型选择不同的计算公式，计算故障点各序电流if Fault_Type=0 I_Fault1=1/Z1(Fault_Node,Fault_Node); I_Fault2=0; I_Fault0=0; elseif Fault_Type=1 I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node); I_Fault2=I_Fault1; I_Fault0=I_Fault1; elseif Fault_Type=2 I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node)+30*60/1152); I_Fault2=I_Fault1; I_Fault0=I_Fault1; elseif Fault_Type=3 I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node); I_Fault2=-I_Fault1; I_Fault0=0; elseif Fault_Type=4 I_Fault1=1/(Z1(Fault_Node,Fault_Node)+Z2(Fault_Node,Fault_Node)*Z0(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node); I_Fault2=-I_Fault1*Z0(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node); I_Fault0=-I_Fault1*Z2(Fault_Node,Fault_Node)/(Z2(Fault_Node,Fault_Node)+Z0(Fault_Node,Fault_Node); end; str=短路电流;strI_Fault1I_Fault2I_Fault0 %故障电流序分量Iabc=T*I_Fault1 I_Fault2 I_Fault0. %相分量相量Iabc_effective=abs(Iabc) %有效值 %-第二步：计算各个节点的序电压，相电压-%故障分量 Ifault_node1=zeros(bus_Num1,1);Ifault_node2=zeros(bus_Num1,1);Ifault_node0=zeros(bus_Num1,1);for m=1:bus_Num1 if m=Fault_Node Ifault_node1(m,1)=-1*I_Fault1; Ifault_node2(m,1)=-1*I_Fault2; Ifault_node0(m,1)=-1*I_Fault0; else Ifault_node1(m,1)=0; Ifault_node2(m,1)=0; Ifault_node0(m,1)=0; endend Ufault_node1=zeros(bus_Num1,1);Ufault_node2=zeros(bus_Num1,1);Ufault_node0=zeros(bus_Num1,1);Ufault_node1=Y1Ifault_node1;Ufault_node2=Y2Ifault_node2;Ufault_node0=Y0Ifault_node0; %正常分量 Unormal_node1=ones(bus_Num1,1);Unormal_node2=zeros(bus_Num1,1);Unormal_node0=zeros(bus_Num1,1); %故障后电压值=故障分量+正常分量Uall_node1=Ufault_node1+Unormal_node1Uall_node2=Ufault_node2+Unormal_node2Uall_node0=Ufault_node0+Unormal_node0 %序分量str=各节点相电压;strUabc=T*Uall_node1 Uall_node2 Uall_node0. %相量Uabc_effective=abs(Uabc) %有效值 %-第三步：计算各个支路的序电流，相电流-%for k=1:branch_Num1 m=branch1_No1(k); n=branch1_No2(k); Ibranch1(k)=-(Uall_node1(m)-Uall_node1(n)*Y1(m,n); Ibranch2(k)=-(Uall_node2(m)-Uall_node2(n)*Y1(m,n); Ibranch0(k)=-(Uall_node0(m)-Uall_node0(n)*Y0(m,n); %或者不加符号，而是除以z0_branch(k),像第五步那样子 mn(1)=m; mn(2)=n; mn %显示支路号 str= 该支路的序电流 Ibranch1(k) Ibranch2(k) Ibranch0(k) str= 该支路的相电流 Ibranchabc=T*Ibranch1(k) Ibranch2(k) Ibranch0(k). %相量 Ibranchabc_effective=abs(Ibranchabc) end %-第四步：计算发电机节点的序电压，相电压-%Igfault_node1=zeros(bus_Num2,1);Igfault_node2=zeros(bus_Num2,1);Igfault_node0=zeros(bus_Num2,1);for m=1:bus_Num2 if m=Fault_Node Igfault_node1(m,1)=-1*I_Fault1; Igfault_node2(m,1)=-1*I_Fault2; Igfault_node0(m,1)=-1*I_Fault0; else Igfault_node1(m,1)=0; Igfault_node2(m,1)=0; Igfault_node0(m,1)=0; endend Ugfault_node1=zeros(bus_Num2,1);Ugfault_node2=zeros(bus_Num2,1);Ugfault_node1=YY1Igfault_node1;Ugfault_node2=YY2Igfault_node2; %正常分量 Ugnormal_node1=ones(bus_Num2,1);Ugnormal_node2=zeros(bus_Num2,1); %故障后电压值=故障分量+正常分量Uggall_node1=Ugfault_node1+Ugnormal_node1;Uggall_node2=Ugfault_node2+Ugnormal_node2; %Uggall_node表示未经修正的电压值 %发电机节点的序电压for m=(bus_Num1+1):bus_Num2 m Ugall_node1(m)=Uggall_node1(m)*(j*0.5+sqrt(3)/2) Ugall_node2(m)=Uggall_node2(m)*(-j*0.5+sqrt(3)/2)end 发电机节点相电压 for m=(bus_Num1+1):bus_Num2 m Ugabc=T*Ugall_node1(m) Ugall_node2(m) 0. %Uabc表示p节点的a,b,c三相电压向量 Ugabc_effective=abs(Ugabc) %Uabc表示p节点的a,b,c三相电压有效值end %-第五步：计算发电机支路的序电流，相电流-% for k=(branch_Num1+1):branch_Num2 m=branch2_No1(k); n=branch2_No2(k); mn(1)=m; mn(2)=n; Iggbranch1(k)= (Uall_node1(m)-Uggall_node1(n)/z1_branch(k) Iggbranch2(k)= (Uall_node2(m)-Uggall_node2(n)/z2_branch(k) %Iggbranch表示未经修正的电流值，即14支路和2-5支路的电流值（这两个支路有零序电流），不是发电机支路的电流值 if m=1 z=0.0525*j; elseif m=2 z=0.0263*j; end mn str=编号对应支路电流 Iggbranch0(k)=(Uall_node0(m)-0)/z Iggbranchabc=T*Iggbranch1(k) Iggbranch2(k) Iggbranch0(k). str=发电机端口电流 Igbranch1(k)=Iggbranch1(k)*(j*0.5+sqrt(3)/2) ; %发电机支路电流序分量 Igbranch2(k)=Iggbranch2(k)*(-j*0.5+sqrt(3)/2); Igbranchabc=T*Igbranch1(k) Igbranch2(k) 0. %发电机支路电流相量 end clear;clc; U=1; %系统电压Eq0=1.972; %空载电动势UG0=1.21; %机端电压Xe=0.504; %线路电抗值Xd=0.982; %同步电抗Xdd=0.344; %暂态电抗Xd_all=1.486; %系统电抗(包含同步电抗)Xdd_all=0.848; %系统电抗(包含暂态电抗)Tj=10; %惯性时间常数Td0=10; %励磁绕组时间常数 Ke=0.1;Kemax=0;deltamax=0;Pmax=0; while(Ke5.7) for delta=49:110 %请同学们求各个参数 d=delta*pi/180; a=Ke2*Xe2/Xd_all2-1; b=2*(Ke2*Xd*Xe*U*cos(d)/Xd_all2+(Eq0+Ke*UG0); c=Ke2*Xd2*U2/Xd_all2-(Eq0+Ke*UG0)2; Eq=-(b-sqrt(b2-4*a*c)/(2*a); Eqq=Eq*Xdd_all/Xd_all+(Xd_all-Xdd_all)*cos(d)*U/Xd_all; UGd=U*Xd*sin(d)/Xd_all; UGq=Eqq*Xe/Xdd_all+Xdd*U*cos(d)/Xdd_all; UG=sqrt(UGd2+UGq2); K1=Eqq*U*cos(d)/Xdd_all+U2*(Xdd_all-Xd_all)*cos(2*d)/(Xdd_all*Xd_all); K2=U*sin(d)/Xdd_all; K3=Xdd_all/Xd_all; K4=(Xd_all-Xdd_all)*U*sin(d)/Xdd_all; K5=UGd*Xd*cos(d)/(UG*Xd_all)-UGq*Xdd*sin(d)/(UG*Xdd_all); K6=UGq*(Xd_all-Xd)/(UG*Xdd_all); A=zeros(3); %请同学们求P186公式7-14状态方程各个参数 A(1,2)=2*pi*50; A(2,1)=-K1/Tj; A(2,3)=-K2/Tj; A(3,1)=-(K4+Ke*K5)/Td0; A(3,3)=-(1/K3+Ke*K6)/Td0; V=eig(A); V1=real(V(1); V2=real(V(2); V3=real(V(3); if (V10)&(V20)&(V30) Ke_result=Ke; delta_result=delta; P_result=Eq*U*sin(delta*pi/180)/Xd_all; if(Pmax=P_result) Pmax=P_result; Kemax=Ke_result; deltamax=delta_result; end; else s=sprintf(Ke=%f, delta=%d, P=%fn,Ke_result,delta_result,P_result); disp(s); break; end; end; Ke=Ke+0.1; end;s=sprintf(Kemax=%f, deltamax=%d, Pmax=%fn,Kemax,deltamax,Pmax);disp(s); clear;clc; f=50; %系统额定频率Tj=8.47; %归算后的发电机惯性时间常数PT=1; %正常运行时发电机向无穷大系统传输的有功功率P2M=0.48; %故障存在时发电机的最大功率P3M=1.38; %故障切除后发电机的最大功率 %下面是利用改进欧拉法进行逐段计算%需要注意故障切除前后电磁功率有跃变h=0.001; %设置步长0.01sDuration=20; %设置计算时段长度2s CutTime=input(输入故障切除时间n); %故障发生时的功角变化过程Delta(1)=33.92; %初始功角Omega(1)=1; %初始转速t(1)=0;for i=1:round(CutTime/h) d(i)=(Omega(i)-1)*360*f;%第一时段开始时d(i)和o(i)的