西电系统仿真大作业

----------系统仿真大作业指导老师：屈胜利班级：学号：姓名：系统仿真上机作业一、计算机辅助系统分析：系统如下图所示：其中r是单位阶跃，GN是非线性器件，G0(s)=)1025.0)(1625.0)(110(s1ssssK）（1．当GN=1、K=40时，用MATLAB画出开环Bode图，求出ωc、θB。由其估计出tr、ts、σ%。解：matlab程序如下：运行结果如下：由上图知：den=[0.1563,6.5156,10.6500,1.0000,0]增益裕量Gm=4.3168相位裕量Pm=10.0158穿越频率Wcg=5.1598增益为0的频率Wcp=2.3975所以得到：ωc=5.1598θB=10.0158仿真波形如下：求解tr、ts、σ%：由bode图：可将系统的开环传递函数化成G(s)=）（11040ss其相应的闭环传递函数为GB(s)=10s^2+s+40，此二阶系统可求出相应的=0.025，wn=2。所以可以估计出：18.67smatlab程序：运行结果为：所以得到：tr=0.7772sts=17.0985sσ%=81.68%经检验：由bode图估算出tr,σ%,ts的结果与正确值的差距不大.2．当GN1，K40时，用MATLAB画出根轨迹图，并求出K40时的闭环极点；由其估计出tr、ts、σ%。解：matlab程序为：运行结果：求系统闭环零极点matlab程序：运行结果：仿真图：由上图：闭环系统的所有极点为-40，-0.2274+2.4146j，-0.2274-2.4146j，-1.083，闭环系统的零点为-1。根据零极点对消原理，可将-1.083与-1对消掉，取主导极点，传递函数为G(s)=)4146.22274.0)(4146.22274.0(s40jsj=882.54548.02^s5.882s由上述G(s)可求得：tr=0.69，ts=15.39s，σ%=74.38%Matlab程序如下：运行结果：（rt）（%=74.38%）（st）与估算结果差距不大。3．当GN=1，K=40时，仿真之，并由仿真结果求出tr、ts、σ%。①用自适应变步长法②用定步长RK-2法，并观察步长与仿真结果的收敛于发散关系，以及仿真精度。解：G0(s))1025.0)(1625.0)(110(s1ssssK）（时：)40)(.61)(10.(s1256ssss）（①采用自适应变步长法：simulink仿真模块连接如下：仿真波形：计算求得：rt=0.69，st=17.14s，%=81.72%Matlab程序：ytr=find(simout(:,1)=1)；rise_time=tout(ytr(1))[ymax,tp]=max(y);peak_time=tout(tp);max_overshoot=ymax-1s=length(tout);whiley(s)0.98&y(s)1.02s=s-1;endsettling_time=tout(s)运行结果：rise_time=0.7630s（rt）max_overshoot=0.8169（%=81.69%）settling_time=17.1290s（st）与估算结果相差不大。②用定步长RK-2法：h取0.01：仿真图为：rt=0.69s，st=17.14s，%=81.72%h取0.02：仿真图：rt=0.70s，st=17.14s，%=81.69%h取0.03：仿真图：rt=0.69s，st=17.13s，%=81.73%h取0.04：rt=0.72s，st=17.16s，%=81.74%h取0.05：rt=0.71s，st=17.15s，%=81.70%h取0.052：综上可知：使用RK—2定步长法时，当h取0.05以上时，将会出现发散，因此，增大步长会使系统稳定性下降。在h小于等于0.05时，对应的系统的性能如下所示：4．令图中的K=40。①GN分别为：。分别仿真之，并由仿真结果求出tr、ts、σ%。解:1）加饱和非线性特性曲线：simulink仿真模块连接如下：步长h上升时间tr调整时间ts超调量σ%0.010.6917.140.81720.020.7017.140.81690.030.6917.130.81660.040.7217.160.81640.050.7117.150.8170仿真波形如下：计算求得：rt=0.91，st=16.36s，%=72.76%2）加死区非线性特性曲线：simulink仿真模块连接如下：仿真波形：计算求得：rt=0.7424，st=29.35s，%=80.67%②GN=1，在G0(s)之后，反馈点之前加。仿真之，并由仿真结果求出tr、ts、σ%。解：simulink仿真模块连接如下：仿真波形：计算如下：rt=0.7270s，st=17.27s，%=89.42%。③对3和4中①、②的tr、ts、σ%比较，并解释差异的原因。答：3中①﹑②仿真方法不同，从而精度不同最后导致性能差异。原因：自适应变步长法中步长hk的大小与y的变化率有关，当系统趋于稳定时，y变化较慢，hk变大；hk为较小值时，定步长RK-2法的ts较大。而在仿真前期，y变化较剧烈，两种方法的hk都比较小，故其tr，超调量相差不大。4中由下表知：①中加入死区非线性与饱和非线性相比，加入饱和非线性的系统性能更好，②的指标小于①，①中延迟大，二者不同说明在系统不同位置加入相同的非线性环节，对系统的影响不同。rtst%①o.9116.3672.76%①0.742429.3580.67②0.727017.27s89.42%二、病态系统（stiff）仿真（simulink）1．用自适应变步长法（RK45）仿真之解：选取τ1=1000，τ2=0.001（1）simulink连接图如下：仿真波形：2．用定步长四阶龙格库塔法仿真，并试着搜索收敛的步长h的范围；若找不到h，将τ1减小，τ2增大，用定步长四阶龙格库塔法仿真，寻找h。（2）解：simulink设置为定步长RK-4仿真法仿真波形：h=0.279h=0.28h=0.3结论：由上图分析可知，当h=0.279时，系统仿真收敛，当h0.279时,仿真将发散。3．用病态仿真算法仿真之。解:设置simulink为病态仿真分析比较：当设置1=1000，2=0.001时：系统输出1G()(101)(0.11)sss则求出y(t)=0.1101-10099199ttee由上面最后一张图:当t=6.2时，求出输出y(t)=0.3985用自适应变步长法RK45仿真得到的y(t)=0.3994用定步长RK--4（h=0.279）仿真得到的y(t)=0.3981用病态仿真算法仿真得到的y(t)=0.3987通过比较发现用病态仿真算法得到的暂态响应与理论计算值最接近，而定步长RK-4其次，但步长过大时，系统稳定性降低，自适应变步长法精度最低。并且该仿真系统本身就是一个病态系统maxmin10050s，所以用适应它的病态仿真算法最精确。三．计算机辅助控制器设计：要求：开环θB≥45o,ωc≥4.2,且tr≤0.4s，ts≤1.5s，σ%≤25%。1．开关处于A时，系统性能满足上述要求否？解：（1）开关出于A时，simulink连接图如下：仿真结果：系统bode图：Matlab程序为：程序运行结果：bode图为：求得：增益裕量gm=400.07，相位裕量：pm=17.9930穿越频率wcg=200.0277，截止频率（0dB）wcp=3.0837求tr,ts,%的matlab源程序：运行结果：上升时间tr=0.6447s，峰值时间tp=1.0541s%=0.6045s，调整时间ts=7.2762s所以不符合要求。2．开关处于B时，计算机辅助设计Gc(s)，使系统性能满足上述要求。解:将开关处于B，参照《自动控控原理》（刘丁著），需采用频率法超前矫正，使系统满足指标要求Matlab程序为：num=[10];den=[110];G=tf(num,den);kc=1.05;yPm=45+10;Gc=plsj(G,kc,yPm)%求超前矫正环节Gy_c=feedback(G,1)%求矫正前系统闭环传递函数Gx_c=feedback(G*kc*Gc,1)%求矫正后系统闭环传递函数figure(1)step(Gy_c,'r',10);%画校正前系统阶跃相应曲线holdonstep(Gx_c,'b-.',10)%画校正后系统阶跃相应曲线figure(2)bode(G,'r')%画校正前系统波特图holdonbode(G*Gc*kc,'b')%画校正前系统波特图子函数：functionGc=plsj(G,kc,yPm)G=tf(G);[mag,pha,w]=bode(G*kc);Mag=20*log10(mag);[Gm,Pm.Wcg,Wcp]=margin(G*kc);phi=(yPm-getfield(Pm,'Wcg'))*pi/180;alpha=(1+sin(phi))/(1-sin(phi));Mn=-10*log10(alpha);Wcgn=spline(Mag,w,Mn);T=1/(Wcgn*sqrt(alpha));Tz=alpha*T;Gc=tf([Tz1],[T1]);End运行结果得超前矫正环节为：校正后系统simulink模块连接图为：仿真波形为：较粗的那条为校正后的波形，其性能指标：rt=0.2604s，st=0.3712s，%=22.9534%满足系统指标要求。四.观测站(O点)为测得的某航班数据，当飞机到达某位置P时开始对飞机的相关数据进行记录，在位置P处时间t记为0，当飞机到观测站的距离达到最短时飞机所处的位置记为M点，飞机在t时刻所处的位置与观测站O点的连线到直线OM的夹角记为θ(t)。附表一给出了观测站某次记录数据：A列为记录时间t，其间隔为0.001s；B列为飞机飞行过程中t时刻飞机相对于P点的距离S(t)；C列给出了飞机t时刻角度θ(t)的理论值theta_theory（参考输入），D列给出了观测站实际上观测到的飞机在t时刻时角度θ(t)的观测值theta_observation；E列给出了当把C列数据theta_theory作为某标准二阶伺服系统G(s)的输入信号时该标准二阶伺服系统的输出值theta_output。作业要求：（1）将附表一中的数据导入matlab工作空间，使各列数据都能作为变量使用。解：（2）试根据表格中的数据使用MATLAB完成以下问题1.根据A、B两列数据确定飞机飞行时的理论运行轨迹和飞机的飞行速度；2.根据A、B、C三列数据确定当飞机到达M点时观测站O到M点的距离。解：1.matlab程序;飞行曲线：求解观测站O点到M点的距离。Matlab代码如下：运行结果：先拟合飞机位移与实践的关系，求解出飞机的飞行速度为140m/s，在根据三角形OMP求解出OM变得长度。其中data1（71430）对应飞机飞行到M的时刻。求解得OM的长度为500.0030m（3）设计标准二阶伺服系统G(s)ωn2,要求：S22ζωnωn21．确定合适的ζ值，其单位阶跃响应的性能指标满足:系统的超调量σ%介于4.5%~8.0%之间解：matlab程序：dete=[];foryipsin=0:0.1:1dt=100*exp(-yipsin*pi/(sqrt(1-yipsin^2)));dete=[dete,dt]endyipsin=0:0.1:1;plot(yipsin,dete,'-b')gridonxlabel('ζ取值');ylabel('σ取值%');仿真波形：由上图可确定取ζ为0.65时，可满足σ%介于4.5%~8.0%之间。验证，当=0.65=6.8%时，，满足要求2.使用MATLAB仿真确定参数ωn,使得：当把C列数theta_theory作为该二阶伺服系