IEEE算例完全版

%单机无穷大系统动态过程 用隐式梯形积分法求解编程 tic clear all clc %初值参数 %%%同步发电机参数 Wref=1; %参考转速 W(1)=0.98; %转子角速度，简称转速 TJ(1)=40.28; %同步电机的转子机械惯性时间常数 Xd(1)=0.1460; %d轴同步电抗 Xd1(1)=0.0608; %d轴暂态电抗 Xq(1)=0.0969; %q轴同步电抗 Td0(1)=8.96; %d轴开路暂态时间常数 DX(1)=0.1; %阻尼系数 %%%输电线及变压器参数 XT=0.0576; %变压器阻抗 RL1=0.01; %一回路输电线路电阻 XL1=0.085; %一回路输电线路电抗 Y1=0.088i; %一回路输电线路对地导纳 RL2=0.017; %二回路输电线路电阻 XL2=0.092; %二回路输电线路电抗 Y2=0.079i; %二回路输电线路对地导纳 ZXL=XT*i+((RL1+XL1*i+2/(Y1))*(RL2+XL2*i+2/(Y2)))/(RL1+XL1*i+2/(Y1)+RL2+XL2*i+2/(Y2)); %双回输电线路的pi型等值阻抗 RL=real(ZXL); %等值线路的电阻 XL=imag(ZXL); %等值线路的电抗 %%%原动机及调速器传递函数参数 Ka=20; %转速测量比较放大系数 Kr=3/5; %硬反馈放大系数 Kb=6; %软反馈放大系数 Ts=4; %接力器时间常数 Tb=2.5; %软反馈时间常数 Tw=0.5; %水锤时间常数 %%%励磁系统传递函数参数 TR=0.001; %测量惯性环节时间常数 KA=20; %惯性放大环节放大倍数（！注意励磁的a大写，与调速的a小写区别） TA=0.1; %惯性放大环节时间常数 TF=0.02; %励磁电压负反馈惯性微分环节时间常数 KF=0.03; %励磁电压负反馈惯性微分环节放大倍数 TL=0.5; %直流励磁机传递函数，计及饱和作用的惯性环节时间常数 KL=1; %直流励磁机传递函数，计及饱和作用的惯性环节放大倍数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% u(1)=1.05; %%发电机端电压幅值 J(1)=0.00; %%电压相角 U(1)=u(1)*cos(J(1))+u(1)*sin(J(1))*1i;%机端电压向量 Vx(1)=real(U(1)); %机端电压实部 Vy(1)=imag(U(1)); %机端电压虚部 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% V0(1)=1.0; %%末端电压（无穷大母线处的电压，一般假设为恒值） Ix(1)=(Vx(1)*RL+Vy(1)*XL-V0(1)*RL)/(RL^2+XL^2); Iy(1)=(V0(1)*XL-Vx(1)*XL+Vy(1)*RL)/(RL^2+XL^2); I(1)=Ix(1)+Iy(1)*i; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% ED01=U(1)+I(1)*Xq(1)*i; %%求虚构电势 D(1)=angle(ED01); %%求发电机的转子角度 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% Eq1(1)=real(U(1))*cos(D(1))+imag(U(1))*sin(D(1))+Xd1(1)*(real(I(1))*sin(D(1))-imag(I(1))*cos(D(1))); %%%暂态电势 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% mid1=real(I(1))*cos(D(1))+imag(I(1))*sin(D(1)); %mid1=Iq q轴电流 mid2=real(I(1))*sin(D(1))-imag(I(1))*cos(D(1)); %mid2=Id d轴电流 Pe(1)=mid1*Eq1(1)-(Xd1(1)-Xq(1))*mid2*mid1; %电磁功率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% DW(1)=W(1); %转速 DD(1)=D(1); %功角 Pm(1)=Pe(1); %%电磁功率的初值 DPT(1)=Pm(1)+DX(1)*W(1); %%机械功率的初值 DU(1)=DPT(1); %导水翼开度 DUN(1)=DPT(1); %导水翼参考开度 DJ1(1)=0; %总反馈量 DEq1(1)=Eq1(1); %暂态电势 DV1(1)=sqrt(Vx(1)^2+Vy(1)^2);%电压中间变量v1 DV2(1)=0; %电压中间变量v2 DEfd(1)=Eq1(1)+(Xd(1)-Xd1(1))*(real(I(1))*sin(D(1))-imag(I(1))*cos(D(1)));%励磁电压 DV3(1)=(2*KL-1)*DEfd(1); %电压中间变量v3 VF0(1)=1.025; %%DV3(1)/KA+DV1(1);参考输入电压 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% h=0.02; %时间步长 DT=h; %时间步长 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%不变系数部分 aj(1)=h/(2*TJ(1)+DX(1)*h);%v adet=pi*50*h; %v a1=h/(2*Ts+Kr*h);a2=(2*Ts-Kr*h)/(2*Ts+Kr*h);%v a3=(1-Kb*h/(2*Ts)-h/(2*Tb))/(1+Kb*h/(2*Ts)+h/(2*Tb));a4=(Kb*h/(2*Ts))/(1+Kb*h/(2*Ts)+h/(2*Tb));%v a5=(Tw-h)/(Tw+h);a6=h/(Tw+h);a7=(h/Ts)/(1+h/Tw);%v aq(1)=h/(2*Td0(1)+h); aR=h/(2*TR+h);aF=(2*TF-h)/(2*TF+h);aF1=KF*h/(2*TF+h)/TL;aa=h/(2*TA+h);aL=h/(2*TL+(2*KL-1)*h); %aL1=2*TL/(2*TL+(2*KL-1)*h); %%%%%%%%%%%%% TotalTime=400; R_Step=TotalTime/DT; for k = 1:R_Step %%%%%%%%%%%%% 变系数部分 ii=1 ; WN(ii)=DW(ii)+aj(ii)*(DPT(ii) - Pe(ii) -2*DX(ii)*DW(ii)); DN(ii)=DD(ii)+adet*(DW(ii) - 2*Wref); UN(ii)=a2*DU(ii) + a1*(-Ka*(DW(ii) - 2*Wref) - DJ1(ii) +2*Kr*DUN(ii)); J1N(ii)=a3*DJ1(ii) + a4*(-Ka*(DW(ii) - 2*Wref) - Kr*(DU(ii) - 2*DUN(ii))); PTN(ii)=a5*DPT(ii) +a6*DU(ii)+a7*(Ka*(DW(ii) - 2*Wref) + DJ1(ii) + Kr*(DU(ii) - 2*DUN(ii))); Dmid1=(Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii)))*(Xd(ii) - Xd1(ii)); Eq1N(ii)=DEq1(ii)+aq(ii)*(DEfd(ii) - (2*DEq1(ii) + Dmid1)); V1N(ii)=DV1(ii)+aR*(sqrt(Vx(ii)^2 +Vy(ii)^2) -2*DV1(ii)); V2N(ii)=aF*(DV2(ii)) +aF1*(DV3(ii)) - aF1*(2*KL - 1)*DEfd(ii); V3N(ii)=DV3(ii)+aa*(KA*(2*VF0(ii) - DV1(ii) - DV2(ii)) -2*DV3(ii)); EfdN(ii)=DEfd(ii) +aL*(DV3(ii) - 2*(2*KL-1)*DEfd(ii)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% for k1=1:10 %%%%%%%%%%%%%%% 求取偏差方程 ii=1; jj=1; R(jj)=DW(ii)-aj(ii)*(DPT(ii) - Pe(ii)) - WN(ii); jj=jj+1; R(jj)=DD(ii)-adet*DW(ii)-DN(ii); jj=jj+1; R(jj)=DU(ii)+a1*Ka*DW(ii)+a1*DJ1(ii)-UN(ii); jj=jj+1; R(jj)=DJ1(ii)+a4*Ka*DW(ii)+a4*Kr*DU(ii)-J1N(ii); jj=jj+1; R(jj)=DPT(ii)-(a6+a7*Kr)*DU(ii)-a7*Ka*DW(ii)-a7*DJ1(ii) - PTN(ii); jj=jj+1; R(jj)=DEq1(ii)-aq(ii)*(DEfd(ii) -(Ix(ii)*sin(DD(ii)) -Iy(ii)*cos(DD(ii)))*(Xd(ii) -Xd1(ii))) -Eq1N(ii); jj=jj+1; R(jj)=DV1(ii)-aR*sqrt(Vx(ii)^2 + Vy(ii)^2)-V1N(ii); jj=jj+1; R(jj)=DV2(ii)-aF1*DV3(ii)+aF1*(2*KL - 1)*DEfd(ii) - V2N(ii); jj=jj+1; R(jj)=DV3(ii)+aa*KA*(DV1(ii)+DV2(ii)) - V3N(ii); jj=jj+1; R(jj)=DEfd(ii)-aL*DV3(ii) - EfdN(ii); jj=jj+1; Dmid1=(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii))); Dmid2=(Xd1(ii)-Xq(ii))*(Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii)))*Dmid1; R(jj)=Pe(ii)-Dmid1*DEq1(ii)+Dmid2; jj=jj+1; Dmid1=Vx(ii)*cos(DD(ii))+Vy(ii)*sin(DD(ii)); Dmid2=Xd1(ii)*(Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii))); R(jj)=DEq1(ii) - Dmid1 -Dmid2; jj=jj+1; R(jj)=Vx(ii)*sin(DD(ii)) - Vy(ii)*cos(DD(ii)) - Xq(ii)*(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii))); jj=jj+1; R(jj)=Vx(ii)-Ix(ii)*RL+Iy(ii)*XL-V0(ii); jj=jj+1; R(jj)=Vy(ii)-Ix(ii)*XL-Iy(ii)*RL; %%%%%%%%%%%%% if max(abs(R)) > 1.0e-10 %%%%%%%%%%%%% 求取雅克比矩阵 ii=1; JKB(1,1)=1; JKB(1,5)=-aj(ii); JKB(1,11)=aj(ii); %%%%%%%%%%% JKB(2,1)=-adet; JKB(2,2)=1; %%%%%%%%%%% JKB(3,1)=a1*Ka; JKB(3,3)=1; JKB(3,4)=a1; %%%%%%%%%%% JKB(4,1)=a4*Ka; JKB(4,3)=a4*Kr; JKB(4,4)=1; %%%%%%%%%%% JKB(5,1)=-a7*Ka; JKB(5,3)=-(a6+a7*Kr); JKB(5,4)=-a7; JKB(5,5)=1; %%%%%%%%%%% JKB(6,2)=aq(ii)*(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii)))*(Xd(ii)-Xd1(ii)); JKB(6,6)=1; JKB(6,10)=-aq(ii); JKB(6,12)=aq(ii)*sin(DD(ii))*(Xd(ii)-Xd1(ii)); JKB(6,13)=-aq(ii)*cos(DD(ii))*(Xd(ii)-Xd1(ii)); %%%%%%%%%%% JKB(7,7)=1; JKB(7,14)=-aR*Vx(ii)/sqrt(Vx(ii)^2 + Vy(ii)^2); JKB(7,15)=-aR*Vy(ii)/sqrt(Vx(ii)^2 + Vy(ii)^2); %%%%%%%%%%% JKB(8,8)=1; JKB(8,9)=-aF1; JKB(8,10)=aF1*(2*KL - 1); %%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% JKB(9,7)=aa*KA; JKB(9,8)=aa*KA; JKB(9,9)=1; %%%%%%%%%%% JKB(10,9)=-aL; JKB(10,10)=1; %%%%%%%%%%% Dmid1=(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii)))^2 - (Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii)))^2; JKB(11,2)=-DEq1(ii)*(-Ix(ii)*sin(DD(ii)) + Iy(ii)*cos(DD(ii))) + (Xd1(ii) - Xq(ii))*Dmid1; JKB(11,6)=-(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii))); JKB(11,11)=1; Dmid1=sin(DD(ii))*(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii))) + cos(DD(ii))*(Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii))); JKB(11,12)=-DEq1(ii)*cos(DD(ii)) + (Xd1(ii) - Xq(ii))*Dmid1; Dmid1=-cos(DD(ii))*(Ix(ii)*cos(DD(ii)) + Iy(ii)*sin(DD(ii))) + sin(DD(ii))*(Ix(ii)*sin(DD(ii)) - Iy(ii)*cos(DD(ii))); JKB(11,13)=-DEq1(ii)*sin(DD(ii)) + (Xd1(ii) - Xq(ii))*Dmid1; %%%%%%%%%%% JKB(12,2)=Vx(ii)*sin(DD(ii)) - Vy(ii)*cos(DD(ii)) - Xd1(ii)*(Ix(ii)*cos