系统仿真上机大作业

1、西安电子科技大学 系统仿真上机作业 指导老师：屈胜利 一 计算机辅助系统分析： 系统如下图所示其中为单位阶跃，为非线性器件 要求：1 当，时，用MATLAB画出开环Bode图，求出、。由其估计出、clear num=0 0 0 40 40; den1=conv(1 0,10 1); den2=conv(0.625 1,0.025 1); den=conv(den1,den2) bode(num,den)grid on;xlabel(Frequency rad/s,fontsize,10)title(Bode Diagram G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.0

2、25s+1) mag,phase,w=bode(num,den)；Gm,Pm,Wcg,Wcp=margin(mag,phase,w)所得参数为：den = 0.1563,6.5156 ,10.6500 ,1.0000 ,0 增益裕量Gm = 4.3168 相位裕量Pm = 10.0158 穿越频率Wcg = 5.1598 增益为0的频率Wcp =2.3975 所以可知：=5.1598 =10.0158计算、程序 sys=tf(num,den) sys=feedback(sys,1); y,t=step(sys); ytr=find(y=1); rise_time=t(ytr(1) ymax,t

3、p=max(y); peak_time=t(tp); max_overshoot=ymax-1 s=length(t); while y(s)0.98&y(s)=1); rise_time=t(ytr(1) ymax,tp=max(y); peak_time=t(tp); max_overshoot=ymax-1 s=length(t); while y(s)0.98&y(s)1.02 s=s-1; end settling_time=t(s)结果：rise_time =0.7772 max_overshoot = 0.7438 settling_time = 16.9691 可知：=0.77

4、72s =16.9691s =74.38%3，：仿真之，并由仿真结果求出、答：建立系统模型用自适应变步长方法（ode45）仿真可得：观察上图可知：超调量约为80%、上升时间Tr约为0.5S；调整时间Ts约为18S。定步长RK-2法H=0.08,仿真结果发散h=0.05，仿真收敛H=0.02，收敛H=0.001，收敛，但仿真精度明显降低综上可知：仿真的步长必须选择合理才能更准确的得出仿真结果，步长太大，截断误差大，甚至导致仿真失败，而步长太小，使舍入误差逐渐积累，也导致误差增大。所以步长的合理选择和仿真精度有很大的关系。4令图1的分别为：分别仿真，并计算、当GN为饱和特性时，仿真结果如下： tr

5、=1.37s ts=15.7s、 %=58%。当GN为死区特性时，仿真结果如下： tr=0.49s，ts=34.2s、 %=100%。，在之后，反馈点之前加上，仿真之，并计算、仿真图形： tr=0.43s，ts=17.3s、 %=72%。 对 3 和4 中、的r t 、s t 、 %比较，并解释差异的原因。比较自适应变步长法和定步长RK-2法的仿真结果得出，RK-2的ts较长，但两种方法的 tr、 %相差较小。原因：自适应变步长法中步长hk的大小与y的变化率有关，当系统趋于稳定时，y变化较慢，步长hk变大，而定步长时hk为一个较小值，故定步长RK-2法的ts较大。而在仿真前期,y变化较剧烈，两

6、种方法的步长hk都比较小，故其 tr、 %相差不大。若非线性环节加在G0(s)之前，死区特性仅作用于误差e(s)上，若非线性环节加在G0(s)之后，则死区特性作用于e(s)G0(s)上，在该系统中，G0（S）增益较大，所以后者对原系统的动态性能影响较大，从而 tr、ts、 %与原系统相差较多；二、病态系统（stiff）仿真（simulink）:单位阶跃:1 用自适应变步长法（RK45）仿真之解：当1=102，2=10-2时。仿真结果如下：2.用定步长四阶龙格库塔法仿真，并试着搜索收敛的步长h 的范围；若找不到h，将1 增大，2减小，用定步长四阶龙格库塔法仿真，寻找h。解：当1=102，2=10

7、-2时。收敛的步长h的范围是：h=0.027。当h=0.027时的仿真结果如下：3.用病态仿真算法仿真之以上三问，均打印出仿真曲线，计算暂态响应，并比较讨论之。答：Y（s）=G(s)R(s)=1s(10s+1)(0.1s+1)=1s-10099s+0.1+199s+10 y(t)=-10099e-0.1t+199e-10t+1。当t=10s时，可算出y(t)=0.39,用自适应变步长法（RK45）仿真得到的y(t)=0.164,用定步长四阶龙格库塔法仿真得到的y(t)50，系统严重病态，而RK4-5和RK4若取低步长则仿真速度极慢，若不如此则误差较大。三计算机辅助控制器设计：+-ABB要求：开

8、环,且，1 开关处于A时，系统性能满足上述要求否？仿真波形：Bode图为：由仿真波形可以知道，开关接A时，不能满足指标。2.开关处于B时，计算机辅助设计，使系统性能满足上述要求（1）基于频率法的串联超前校正主程序：function=question3_2(k)num=k*1;den=conv(1,0,1,1);G=tf(num,den);kc=1;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,5

9、);hold onstep(Gx_c,b,5);grid onfigure(2)bode(G,r)hold on bode(G*kc*Gc,b)grid onGm,Pm,Wcg,Wcp=margin(G*kc*Gc);r=PmWc=Wcpy,t=step(Gx_c);ytr=find(y=1);rise_time=t(ytr(1);tr=rise_timeymax,tp=max(y); peak_time=t(tp); max_overshoot=ymax-1; max_over=max_overshoot s=length(t);while y(s)0.98&y(s) t(2)ans = 1

10、.000000000000000e-03 s(2)ans = 0.140000000000000 theta_theory(2)ans = -1.520837232808838theta_observation(2)ans = -7.075789556117257e-05 theta_output(2)ans =-7.573799703386787e-05(2) 试根据表格中的数据使用MATLAB完成以下问题l 根据A、B两列数据确定飞机飞行时的理论运行轨迹和飞机的飞行速度；l 根据A、B、C三列数据确定当飞机到达M点时观测站O到M点的距离。编写脚本文件plane_rate.mplot(t,s

11、);grid on;%作飞机位移关于时间t的图象，可看出飞机是匀速运动i=2:143001;rate(i)=s(i)./t(i);%求运行轨迹中各点的速率，可看出运动速率恒为140m/si=1:143001;abs_theta_theory(i) = abs(theta_theory(i);%取观察站测得角度的绝对值min_theta_theory k = min(abs_theta_theory);%当飞机到达M点时，观察站测得角度应为最小角，故求得此最小角度时%对应的时间t(k)distance_OM = s(k)*cot(abs_theta_theory(1)%当时间为t(k)时，对应的

12、位移为s(k)，tan(theta)=PM/OM%其中PM=s(k),theta=theta_theory的最初的值，由此可求得OM的距离运行结果： planeratedistance_OM = 5.000030000000002e+02故飞机飞行速率为140m/s,OM间距离为500m。（3）设计标准二阶伺服系统G(s) 解：1确定合适的值，其单位阶跃响应的性能指标满足:系统的超调量%介于4.5%8.0%之间因为%=e-1-2*100%，使%介于4.5%8.0%之间则可取=0.68，此时%= %5.4280251339488，满足要求。 2.使用MATLAB仿真确定参数n,使得：当把C列数据

13、theta_theory作为该二阶伺服系统的输入时，系统的输出尽可能的接近E列所给出的theta_output。求出系统传递函数参数。用simulink设计标准二阶伺服系统，其中w的值待定设计的框图如下：再编写脚本文件，如下：for n=1:20;%取不同的w值时，求输出的theta_output sim(erjiesifuxitong);%运行simulink中的框图文件， %获得输出的的p_theta_output值，为一个长度为143001的列向量 for j=1:143001; d_value(j)=abs(p_theta_output(j)-theta_output(j); %求每组

14、输出的p_theta_output与理论值的差 end abs_dif(n)=sum(d_value);%取误差和，存到向量ads_dif中，再依次比较向量中 %的值，当取最小值时说明p_theta_output与理论值最接近，此时 %的w值即为标准w值end当w取120时，得到的abs_dif各值为： abs_difabs_dif = 1.0e+03 * Columns 1 through 4 6.019400510911621 2.739540882427722 1.644088863073810 1.100263457041429 Columns 5 through 8 0.773241

15、725337181 0.549005514111091 0.384051198965172 0.256977609508468 Columns 9 through 12 0.156862800155772 0.079109612117110 0.036344028042700 0.064234872892925 Columns 13 through 16 0.109981064663135 0.149590657196720 0.183731144417703 0.213275684364825 Columns 17 through 20 0.238998060855472 0.2615401

16、30309759 0.281422897899948 0.299065677805618由数据可看出，当w=11时，取得最小的差值，再取w=1012之间的数，求得其中最小差值。得到的simulink框图为：3.将C列theta_theory作为设计好的标准二阶伺服系统的输入信号，计算输出相对于输入的相对误差；作出输出随时间变化的曲线，以及相对误差随时间变化的曲线。作输出变化曲线和相对误差曲线求相对误差的脚本文件为：sim(erjiesifuxitong);%运行二阶伺服系统for i=1:143001;abs_dif(i)=abs(p_theta_output(i)-theta_theory(

17、i);%求输出p_theta_output的绝对误差relative_error(i)=abs_dif(i)/abs(theta_theory(i);%求输出的相对误差End作输出变化曲线为：从图中可见系统输出的曲线与给出的理论曲线几乎重合，可见选择的系统参数比较准确。当作相对误差曲线时，发现为如下图所示：接近于冲激函数。因为时间轴选用的过宽，无法看出细节，故用axis(0,2,0,2)指令截取相应的曲线，就能看清楚了：a) 时间轴为02时plot(t,relative_error);axis(0,2,0,1)b) 时间轴为270时，图形几乎无变化，即相对误差近似为0.c) 时间轴为7073时

18、，plot(t,relative_error);axis(70,73,0,40)由图可看出相对误差变化图在7073之间有一个猛增的过程d) 时间轴为73143时，plot(t,relative_error);axis(73,143,0,1)分析图可看出，在70s73s之间，相对误差有一个显著的冲激式上升过程，分析原因：相对误差的求取公式为Poutput-thetatheorythetatheory*100%而theta_theory的值在70s73s时近似取到了零，而theta_theory在分母上，故此时相对误差值会变得极大。当theta_theory取值远离零值时，由于绝对误差接近于0，故

19、相对误差也接近于0。而在0s2S区间可看到一个明显的振荡过程，分析原因：因为系统为二阶系统，输入的起始阶段响应会有超调量，但当经过一段时间系统稳定后，就不会再发生振荡。4.在上述设计好的伺服系统中加入合适的干扰信号和适当的非线性环节，使得当把C列数据theta_theory作为输入信号时系统的输出尽可能接近D列数据(t)的观测值theta_observation)。试：1) 确定干扰信号和非线性环节的相关参数，并画出系统simulink模型。l 作出theta_theory的时间曲线和theta_observation的时间曲线如图，蓝色为theta_theory的时间曲线，红色为theta_

20、observation的时间曲线。编写M文件，使每选定一种非线性环节或添加一种干扰都能得到对应的曲线，再将对应曲线与给定的观察曲线作比较，即能找到对应的参数。Thetaobservation.m%plot(t,theta_theory);hold on;%作理论曲线plot(t,theta_observation,r);hold on;%作观察所得曲线sim(erjiesifuxitong);%运行仿真框图plot(t,p_theta_output);%作仿真结果曲线hold on;首先在前向通道中加一个死区非线性环节，得到如下图所示波形：此时输出曲线与观察曲线相差较大，故不正确。再将死区非线

21、性放到回路外：得到曲线：输出曲线在中间停滞区与观察曲线近似，为得到死区宽度，将题设所给的理论曲线与观察曲线打印出来如图：从上图中可看出，死区宽度为1，故将死区的起始与截止参数改为如下所示：此时能得到与观察曲线非线性弯曲部分重合的曲线：此即找到非线性环节。l 为找到干扰环节，我首先将观察曲线与输出曲线作对比，如下图：为求出合适的干扰信号使蓝色曲线接近红色曲线，先将两曲线作差(i=1:143001;diff(i)=p_theta_output(i)-theta_observation(i);plot(t,diff);)得到如下曲线：由图可看出，差值diff的时间曲线近似为一个正弦波形，故应输入一个

22、正弦形状的干扰信号。为求出干扰信号的具体差值，我再对diff进行一次线性拟合，编写函数文件如下：function fun_sin = xianxingnihe(x,t)for i=1:143001 fun_sin(i)=x(1)*sin(x(2)*t(i)+x(3)*pi);end再输入指令行x(0)=0.1 0.1 0.1x,resnorm=lsqcurvefit(xianxingnihe,x0,t,diff),即对diff曲线进行输入为时间t的线性拟合，拟合结果如下：x =-0.049815057234210 0.100054669976175 -0.006617837055889resn

23、orm = 1.820876780757486再作出拟合所得曲线与真实的diff曲线作比较：plot(t,diff);hold on;plot(t,xianxingnihe(x,t),r);得到图为：由图可知，拟合结果近似接近于原曲线，因为原曲线幅值较小，故将此拟合结果作为干扰信号输入。输入系统框图如图：其中正弦波参数采用拟合曲线得到的参数，如下输出此框图的结果为：由图可看出，两条曲线近似重合，故此系统即为标准的输入系统。2) 画出系统的输出随时间变化的曲线。plot(t,theta_observation);hold on;%作理论观察曲线sim(erjiesifuxitong);%运行仿真

24、框图plot(t,p_theta_output,r);%作仿真结果曲线，为红色hold on;3) 画出输出相对于输入theta_theory的相对误差随时间变化的曲线，以及输出相对于D列数据(t)的观测值theta_observatuon)的相对误差曲线。输出相对于输入theta_theory的相对误差随时间变化的曲线：i=1:143001;diff_a(i)=abs(p_theta_output(i)-theta_theory(i)./abs(theta_theory(i);plot(t,diff_a);hold on;axis(0 143 0 1);l 输出相对于D列数据(t)的观测值theta_observatuon)的相对误差曲线i=1:143001;diff_a(i)=abs(p_theta_output(i)-theta_observation(i)./abs(theta_observation(i);plot(t,diff_a);hold on;axis(0 143 0 0.1);五. 幕布系统由电机驱动电路、电机及幕布三部分组成。通过电机驱动电路给电机输入电信号，电机轴转动带动幕布