西电系统仿真大作业

1、 -系统仿真大作业 指导老师： 屈胜利 班 级： 学 号： 姓 名： 系统仿真上机作业一、计算机辅助系统分析：系统如下图所示：其中 r 是单位阶跃，GN 是非线性器件，G0 (s)= 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

2、w=5.1598 b qB =10.0158仿真波形如下：求解tr 、 ts 、 %：由bode图：可将系统的开环传递函数化成G(s)=+，其相应的闭环传递函数为GB(s)=10s2+s+40，此二阶系统可求出相应的x=0.025，wn =2。所以可以估计出： 18.67smatlab程序：运行结果为： 所以得到：tr=0.7772s ts=17.0985s % =81.68% 经检验：由bode图估算出tr, % ,ts的结果与正确值的差距不大.2 当 GN = 1， K = 40 时，用 MATLAB 画出根轨迹图，并求出 K = 40 时的闭环极点；由其估计出 tr 、 ts 、 % 。

3、解：matlab程序为：运行结果： 求系统闭环零极点matlab程序：运行结果： 仿真图：由上图：闭环系统的所有极点为-40，-0.2274+ 2.4146j，-0.2274-2.4146j， -1.083，闭环系统的零点为-1。根据零极点对消原理，可将-1.083与-1对消掉，取主导极点，传递函数为G(s)= = 由上述G(s)可求得： tr=0.69， ts=15.39s， %=74.38%Matlab程序如下： 运行结果： （） （=74.38%） （）与估算结果差距不大。3当 GN = 1， K = 40 时，仿真之，并由仿真结果求出 tr 、 ts 、 % 。用自适应变步长法用定步长

4、 RK-2 法，并观察步长与仿真结果的收敛于发散关系，以及仿真精度。解：G0 (s) = K=40时： =采用自适应变步长法：simulink仿真模块连接如下：仿真波形：计算求得： =0.69，=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);while y(s)0.98&y(s)1.02 s=s-1;end settling_time=tout(s)运行结果： rise

5、_time = 0.7630s （） max_overshoot =0.8169 （=81.69%） settling_time =17.1290 s （）与估算结果相差不大。用定步长 RK-2 法：h取0.01：仿真图为：=0.69s，=17.14s，=81.72%h取0.02：仿真图：=0.70s，=17.14s，=81.69%h取0.03：仿真图：=0.69s，=17.13s，=81.73%h取0.04：=0.72s，=17.16s，=81.74%h取0.05：=0.71s，=17.15s，=81.70%h取0.052：综上可知：使用RK2定步长法时，当h取0.05以上时，将会出现发散，

6、因此，增大步长会使系统稳定性下降。在h小于等于0.05时，对应的系统的性能如下所示：步长h上升时间tr调整时间ts 超调量%0.010.6917.140.81720.020.7017.14 0.81690.030.6917.130.81660.040.7217.160.81640.050.7117.150.81704令图中的 K = 40 。 GN 分别为：。分别仿真之，并由仿真结果求出 tr 、 ts 、 % 。解:1）加饱和非线性特性曲线： simulink仿真模块连接如下：仿真波形如下：计算求得： =0.91，=16.36s，=72.76%2）加死区非线性特性曲线：simulink仿真模

7、块连接如下：仿真波形：计算求得： =0.7424，=29.35s，=80.67% GN = 1，在 G0 (s) 之后，反馈点之前加。仿真之，并由仿真结果求出 tr 、 ts 、 % 。解：simulink仿真模块连接如下：仿真波形：计算如下：=0.7270s， =17.27s， =89.42%。对 3 和 4 中、的 tr 、 ts 、 % 比较，并解释差异的原因。答：3中仿真方法不同，从而精度不同最后导致性能差异。原因： 自适应变步长法中步长hk的大小与y的变化率有关，当系统趋于稳定时，y变化较慢，hk变大；hk为较小值时，定步长RK-2法的ts较大。而在仿真前期，y变化较剧烈，两种方法的

8、hk都比较小，故其tr，超调量相差不大。4中由下表知：中加入死区非线性与饱和非线性相比，加入饱和非线性的系统性能更好，的指标小于，中延迟大，二者不同说明在系统不同位置加入相同的非线性环节，对系统的影响不同。o.9116.3672.76%0.742429.3580.670.727017.27s89.42%二、病态系统（stiff）仿真（simulink）1用自适应变步长法（RK45）仿真之解：选取1 =1000，2=0.001（1） simulink连接图如下：仿真波形：2 用定步长四阶龙格库塔法仿真，并试着搜索收敛的步长 h 的范围；若找不到h，将1 减小， 2 增大，用定步长四阶龙格库塔法仿

9、真，寻找 h。（2） 解：simulink设置为定步长RK-4仿真法仿真波形：h=0.279h=0.28h=0.3结论：由上图分析可知，当h0.279时, 仿真将发散。3 用病态仿真算法仿真之。解: 设置simulink为病态仿真分析比较：当设置=1000，=0.001时：系统输出则求出y(t)=由上面最后一张图:当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通过比较发现用病态仿真算法得到的暂态响应与理论计算值最接近，而定步长R

10、K-4其次，但步长过大时，系统稳定性降低，自适应变步长法精度最低。并且该仿真系统本身就是一个病态系统 ，所以用适应它的病态仿真算法最精确。三计算机辅助控制器设计：要求：开环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,

11、 % 的matlab源程序： 运行结果： 上升时间tr =0.6447s ， 峰值时间tp =1.0541s %=0.6045s ，调整时间ts =7.2762s所以不符合要求。2 开关处于 B 时，计算机辅助设计 Gc (s) ，使系统性能满足上述要求。解:将开关处于B，参照自动控控原理（刘丁著），需采用频率法超前矫正，使系统满足指标要求Matlab程序为：num=10;den=1 1 0;G=tf(num,den);kc=1.05;yPm=45+10;Gc=plsj(G,kc,yPm) %求超前矫正环节Gy_c=feedback(G,1) %求矫正前系统闭环传递函数Gx_c=feedbac

12、k(G*kc*Gc,1) %求矫正后系统闭环传递函数 figure(1)step(Gy_c,r,10); %画校正前系统阶跃相应曲线hold onstep(Gx_c,b-.,10) %画校正后系统阶跃相应曲线figure(2)bode(G,r) %画校正前系统波特图hold onbode(G*Gc*kc,b) %画校正前系统波特图子函数：function Gc=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/1

13、80;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(Tz 1,T 1);End运行结果得超前矫正环节为：校正后系统simulink模块连接图为：仿真波形为：较粗的那条为校正后的波形，其性能指标：=0.2604s，=0.3712s，=22.9534%满足系统指标要求。四 .观测站(O 点)为测得的某航班数据，当飞机到达某位置 P 时开始对飞机的相关数据进行记录，在位置 P 处时间 t 记为 0，当飞机到观测站的距离达到最短

14、时飞机所处的位置记为 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。作

15、业要求：（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

16、（3） 设计标准二阶伺服系统 G(s) =n2,要求：S2 + 2n + n21 确定合适的值，其单位阶跃响应的性能指标满足:系统的超调量%介于 4.5%8.0%之间 解：matlab程序： dete=;for yipsin=0:0.1:1 dt=100*exp(-yipsin*pi/(sqrt(1-yipsin2); dete=dete,dtendyipsin=0:0.1:1;plot(yipsin,dete,-b)grid on xlabel(取值); ylabel(取值%);仿真波形： 由上图可确定取为0.65时，可满足%介于4.5%8.0%之间。验证，当，满足要求2. 使用 MATLA

17、B 仿真确定参数n,使得：当把 C 列数theta_theory 作为该二阶伺服系统的输入时，系统的输出尽可能的接近 E 列所给出的 theta_output。求出系统传递函数参数。解：simulink仿真模块连接为：仿真波形：当wn=1时：Wn=2时Wn=5由上图：wn=5时，两条曲线基本吻合，只是在起始阶段存在误差，故wn大于5时，只需观察两条曲线在起始阶段的差异即可。Wn=7Wn=9Wn=10Wn=12由上图可知，当wn大于9时，两条曲线的差异开始增大，故当wn取9时，系统的输出更接近theta_output。=0.65，wn=9传递函数为：G(s)=n / S2 + 2n + n2=9

18、/（s2+11.7+81）3. 将 C 列 theta_theory 作为设计好的标准二阶伺服系统的输入信号，计算输出相对于输入的相对误差；作出输出随时间变化的曲线，以及相对误差随时间变化的曲线。解：用matlab作图：。不会。（4） 在上述设计好的伺服系统中加入合适的干扰信号和适当的非线性环节，使得当把 C列数据 theta_theory 作为输入信号时系统的输出尽可能接近 D 列数据(t)的观测值theta_observatuon)。试：1. 确定干扰信号和非线性环节的相关参数，并画出系统 simulink 模型。2. 画出系统的输出随时间变化的曲线。3. 画出输出相对于输入 theta_theory 的相对误差随时间变化的曲线，以及输出相对于 D 列数据(t)的观测值 theta_observatuon)的相对误差曲线。解：1非线性环节相关参