華北電力大學(xué)電力系統(tǒng)暫態(tài)上機(jī)程序設(shè)計(jì)完整版.doc

clc;clear; %-輸入已知條件-%bus_Num1=3; %節(jié)點(diǎn)數(shù)bus_Num2=5; %包括發(fā)電機(jī)節(jié)點(diǎn)的節(jié)點(diǎn)數(shù)branch_Num1=3; %線路數(shù)branch_Num2=5; %包括發(fā)電機(jī)支路的支路數(shù)branch1_No1=1,1,2;branch1_No2=2,3,3;branch2_No1=1,1,2,1,2;branch2_No2=2,3,3,4,5; %每條支路始末端編號(hào) %輸入支路各序阻抗，z1_branch表示支路正序阻抗，z2_branch表示支路負(fù)序阻抗，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; %-第一步：不考慮發(fā)電機(jī)節(jié)點(diǎn)計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣-%節(jié)點(diǎn)導(dǎo)納矩陣，Y1表示不計(jì)發(fā)電機(jī)節(jié)點(diǎn)的正序網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣，Y2表示不計(jì)發(fā)電機(jī)節(jié)點(diǎn)的負(fù)序網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣，Y0表示不計(jì)發(fā)電機(jī)節(jié)點(diǎn)的零序網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣，Y1=zeros(bus_Num1);%全為0矩陣（3*3）Y1(1,1)=1/(j*0.15);Y1(2,2)=1/(j*0.075);%請(qǐng)同學(xué)們求正序節(jié)點(diǎn)導(dǎo)納矩陣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;%負(fù)序等于正序 Y0=zeros(bus_Num1);Y0(1,1)=1/(j*0.0525);Y0(2,2)=1/(j*0.0263);%請(qǐng)同學(xué)們求零序節(jié)點(diǎn)導(dǎo)納矩陣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 %-第二步：考慮發(fā)電機(jī)節(jié)點(diǎn)計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣-%節(jié)點(diǎn)導(dǎo)納矩陣，YY1表示計(jì)及發(fā)電機(jī)節(jié)點(diǎn)的正序網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣，YY2表示計(jì)及發(fā)電機(jī)節(jié)點(diǎn)的負(fù)序網(wǎng)絡(luò)節(jié)點(diǎn)導(dǎo)納陣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);%請(qǐng)同學(xué)們求正序節(jié)點(diǎn)導(dǎo)納矩陣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;%負(fù)序等于正序 %-第三步：計(jì)算節(jié)點(diǎn)阻抗矩陣-%Z1=inv(Y1);Z2=inv(Y2);Z0=inv(Y0);clc;clear; %數(shù)據(jù)來(lái)源于教材電力系統(tǒng)暫態(tài)分析P77例（3-4），P143例（5-7） Matrix %計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣、節(jié)點(diǎn)阻抗矩陣，形成全局變量 Fault_Node=input(輸入短路點(diǎn)編號(hào) ;n Fault_Node=); Fault_Type=input(輸入短路類型 ;n(1)Fault_Type=0為三相短路;n(2)Fault_Type=1為a相接地短路;n(3)Fault_Type=2為a相經(jīng)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為對(duì)稱分量法的變換矩陣，見(jiàn)P87公式（4-4） a2 a 1 a a2 1; %-第一步：計(jì)算短路點(diǎn)的序電流，相電流-%計(jì)根據(jù)故障類型選擇不同的計(jì)算公式，計(jì)算故障點(diǎn)各序電流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) %有效值 %-第二步：計(jì)算各個(gè)節(jié)點(diǎn)的序電壓，相電壓-%故障分量 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=各節(jié)點(diǎn)相電壓;strUabc=T*Uall_node1 Uall_node2 Uall_node0. %相量Uabc_effective=abs(Uabc) %有效值 %-第三步：計(jì)算各個(gè)支路的序電流，相電流-%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); %或者不加符號(hào)，而是除以z0_branch(k),像第五步那樣子 mn(1)=m; mn(2)=n; mn %顯示支路號(hào) str= 該支路的序電流 Ibranch1(k) Ibranch2(k) Ibranch0(k) str= 該支路的相電流 Ibranchabc=T*Ibranch1(k) Ibranch2(k) Ibranch0(k). %相量 Ibranchabc_effective=abs(Ibranchabc) end %-第四步：計(jì)算發(fā)電機(jī)節(jié)點(diǎn)的序電壓，相電壓-%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表示未經(jīng)修正的電壓值 %發(fā)電機(jī)節(jié)點(diǎn)的序電壓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 發(fā)電機(jī)節(jié)點(diǎn)相電壓 for m=(bus_Num1+1):bus_Num2 m Ugabc=T*Ugall_node1(m) Ugall_node2(m) 0. %Uabc表示p節(jié)點(diǎn)的a,b,c三相電壓向量 Ugabc_effective=abs(Ugabc) %Uabc表示p節(jié)點(diǎn)的a,b,c三相電壓有效值end %-第五步：計(jì)算發(fā)電機(jī)支路的序電流，相電流-% 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表示未經(jīng)修正的電流值，即14支路和2-5支路的電流值（這兩個(gè)支路有零序電流），不是發(fā)電機(jī)支路的電流值 if m=1 z=0.0525*j; elseif m=2 z=0.0263*j; end mn str=編號(hào)對(duì)應(yīng)支路電流 Iggbranch0(k)=(Uall_node0(m)-0)/z Iggbranchabc=T*Iggbranch1(k) Iggbranch2(k) Iggbranch0(k). str=發(fā)電機(jī)端口電流 Igbranch1(k)=Iggbranch1(k)*(j*0.5+sqrt(3)/2) ; %發(fā)電機(jī)支路電流序分量 Igbranch2(k)=Iggbranch2(k)*(-j*0.5+sqrt(3)/2); Igbranchabc=T*Igbranch1(k) Igbranch2(k) 0. %發(fā)電機(jī)支路電流相量 end clear;clc; U=1; %系統(tǒng)電壓Eq0=1.972; %空載電動(dòng)勢(shì)UG0=1.21; %機(jī)端電壓Xe=0.504; %線路電抗值Xd=0.982; %同步電抗Xdd=0.344; %暫態(tài)電抗Xd_all=1.486; %系統(tǒng)電抗(包含同步電抗)Xdd_all=0.848; %系統(tǒng)電抗(包含暫態(tài)電抗)Tj=10; %慣性時(shí)間常數(shù)Td0=10; %勵(lì)磁繞組時(shí)間常數(shù) Ke=0.1;Kemax=0;deltamax=0;Pmax=0; while(Ke5.7) for delta=49:110 %請(qǐng)同學(xué)們求各個(gè)參數(shù) 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); %請(qǐng)同學(xué)們求P186公式7-14狀態(tài)方程各個(gè)參數(shù) 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; %系統(tǒng)額定頻率Tj=8.47; %歸算后的發(fā)電機(jī)慣性時(shí)間常數(shù)PT=1; %正常運(yùn)行時(shí)發(fā)電機(jī)向無(wú)窮大系統(tǒng)傳輸?shù)挠泄β蔖2M=0.48; %故障存在時(shí)發(fā)電機(jī)的最大功率P3M=1.38; %故障切除后發(fā)電機(jī)的最大功率 %下面是利用改進(jìn)歐拉法進(jìn)行逐段計(jì)算%需要注意故障切除前后電磁功率有躍變h=0.001; %設(shè)置步長(zhǎng)0.01sDuration=20; %設(shè)置計(jì)算時(shí)段長(zhǎng)度2s CutTime=input(輸入故障切除時(shí)間n); %故障發(fā)生時(shí)的功角變化過(guò)程Delta(1)=33.92; %初始功角Omega(1)=1; %初始轉(zhuǎn)速t(1)=0;for i=1:round(CutTime/h) d(i)=(Omega(i)-1)*360*f;%第一時(shí)段開(kāi)始時(shí)d(i)和o(i)的