系统仿真上机大作业

1、系统仿真上机大作业 姓名：吴江平 学号：04103043班级：041031老师：屈胜利 一、计算机辅助系统分析：yu+-+re：单位阶跃：非线性器件1 ，时，用MATLAB画出开环Bode图，求出、。由其估计出、M程序：%*%功能：画输入为单位阶跃，开环传递函数G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系统的Bode图及求增益裕量Gm、相位裕量Pm、相位穿越-180°时的交界频率Wcg和增益为0时的频率Wcp%*num=40 40;den=conv(conv(conv(1 0,10 1),0.625 1),0.025 1);w=logspac

2、e(-3,3,100); mag phase w=bode(num,den,w); %返回系统频率特性的幅值和相角向量和对应的频率向量magdB=20*log10(mag); %幅值转换成dB值subplot(2,1,1)semilogx(w,magdB,'k-','LineWidth',2) %画幅频特性Bode图grid onxlabel('Frequecy(rad/s)','FontSize',14)ylabel('Gain(dB)','FontSize',14)h1=title('B

3、ode Diagram of G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)')set(h1,'FontWeigh','bold','FontSize',12)subplot(2,1,2)semilogx(w,phase,'k-','LineWidth',2) %画相频特性Bode图grid onxlabel('Frequecy(rad/s)','FontSize',14)ylabel('Phase(deg)','F

4、ontSize',14)Gm Pm Wcg Wcp=margin(mag,phase,w)%*%功能：求输入为单位阶跃，开环传递函数G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系统的上升时间，峰值时间，超调量和调整时间%*sys=tf(num,den); %建立系统的开环传递函数模型sys=feedback(sys,1); %建立系统的闭环传递函数模型y,t=step(sys); %求出该系统的单位阶跃响应ytr=find(y>=1);rise_time=t(ytr(1) %计算上升时间ymax,tp=max(y);peak_time=t(

5、tp) %计算峰值时间max_overshoot=ymax-1 %计算超调量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %计算调整时间figureh2=plot(t,y,'k',t,ones(length(t),1),'k-.') %绘制响应曲线set(h2(1),'LineWidth',2)xlabel('Time(sec)','FontSize',14)ylabel('Reponse',

6、'FontSize',14)h2=title('Plot of Unit-step Reponse of G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)');set(h2,'FontWeigh','bold','FontSize',12)结果：Gm Pm Wcg Wcp=4.3221 10.0223 5.1653 2.3981rise_time peak_timel max_overshoot setting_time=0.7772 1.2953 0.8168 17.2280所以

7、 =Wcp=2.3981，b=Pm=10.0223 =0.7772，=17.2280，=0.8168截图：2，时，用MATLAB画出根轨迹图，并求出时的闭环极点；由其估计出、M程序：%*%功能：画输入为单位阶跃，开环传递函数G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系统的根轨迹，返回闭环零极点%*num=40 40;den=conv(conv(conv(1 0,10 1),0.625 1),0.025 1);sys=tf(num,den); rlocus(sys)sys=feedback(sys,1); p,z=pzmap(sys) %返回零极点结果：

8、p = -40.1616 -0.2274 + 2.4146i -0.2274 - 2.4146i -1.0837 =0.7772，=17.2280，=0.8168截图： 3、，时，仿真之，并由仿真结果求出，自适应变步长法M程序：sim('f01031'); %进行模型仿真plot(ScopeData(:,1),ScopeData(:,2) %ScopeData为示波器输出到工作空间的矩阵，其中第一列为仿真时间，第二列为仿真波形grid onh=title('自适应变步长ode45法');set(h,'FontWeigh','bold

9、9;,'Color','red,'FontSize',20)结果： =0.685，=17.14，=0.8173截图：定步长二阶RK法截图：当步长大于或等于0.050时发散，否则收敛。4令图1的分别为：分别仿真之，并计算、（从仿真结果）结果： =0.895，=15.14，=0.7142结果： =0.746，=19.60，=0.8193，在之后，反馈点之前加上，仿真之，并计算、结果： =0.721，=14.94，=0.8973对3；4、的、比较，并解释差异的原因 自适应变步长比用定步长法tr,ts,小，因为自适应变步长可自动调整步长。受非线性环节的影响，tr

10、,ts,都相对变大。二、病态系统（stiff）仿真（simulink）:单位阶跃1 用自适应变步长法（RK45）仿真之2 用定步长四阶RK法仿真，并试着搜索收敛的步长h的范围步长0.1时发散0.01发散0.0020.00294稳定。3 用病态仿真之以上三问，均打印出仿真曲线，计算暂态响应，并比较讨论之1.>,(=0.421，=0.981，=0.941)2.>,(=0.421，=1.012，=0.952)3.>,(=0.426，=0.981，=0.948)三计算机辅助控制器设计：+-ABB要求：开环,且，1 开关处于A时，系统性能满足上述要求否？M程序：num=10;den=1

11、 1 0;sys=tf(num,den); sys=feedback(sys,1); y,t=step(sys); ytr=find(y>=1);rise_time=t(ytr(1) %计算上升时间ymax,tp=max(y);max_overshoot=ymax-1 %计算超调量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %计算调整时间bode(num,den)grid onmag phase w=bode(num,den); %返回系统频率特性的幅值和相角向量和对应的频率向量

12、Gm Pm Wcg Wcp=margin(mag,phase,w) %求增益裕量Gm、相位裕量Pm、相位穿越-180°时的交界频率Wcg和增益为0时的频率Wcp结果：Gm Pm Wcg Wcp= 4.0007e+003 18.0184 200.0277 3.0832rise_time max_overshoot setting_time= 0.5961 0.6042 7.3516所以 =Wcp=3.0832，b=Pm=18.0184 =0.5962，=7.3516，=0.6042由此可以看出，开关处于A时，系统性能显然不能满足以下要求：开环,且，截图：2 开关处于B时，计算机辅助设计

13、，使系统性能满足上述要求。基于频率法的串联超前校正M程序：主程序：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);hold onstep(Gx_c,'b',5);grid onfigure(2)bode(G,'r&

14、#39;)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)<1.02 s=s-1;endsettling_time=t(

15、s+1);ts=settling_time 子程序：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/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(Tz,1,T,1);结果：Transfe

16、r function:0.4536 s + 1-0.1126 s + 1 Transfer function: 10-s2 + s + 10 Transfer function: 4.536 s + 10-0.1126 s3 + 1.113 s2 + 5.536 s + 10R = 49.7706Wc = 4.4248tr = 0.3982s（上升时间）max_over = 0.2235（超调量）ts =1.2279s（调整时间）截图：校正前后阶跃响应曲线：校正前后bode图：注释：蓝色线表示校正后，红色线表示校正前上机作业题四仿真在位移x处加阶跃干扰的倒立摆非线性模型（已建好，但仿真不成功）

17、，可以修改各种仿真参数，如算法、误差限、步长、时间等，但是不能修改模型结构，要求仿真成功，说明修改过程，解释物理意义，并说明如何验证结果不稳定是系统本身不稳定还是仿真不当造成的？备注：以下是倒立摆系统建模过程。倒立摆系统如图1所示，它由在导轨上运动的小车和与小车相连的摆杆组成。图 1 倒立摆系统我们假设：M：小车的质量 0.5kgm：摆杆的质量 0.2kgb：小车与基座之间的摩擦系数 0.1N/m/sl：小车与摆杆连接点到摆杆重心之间的距离 0.3mI：摆杆的转动惯量 0.006F：作用在小车上的力 x：小车的位置座标（水平）：摆杆与竖直方向的夹角 将倒立摆系统分为两个独立的部分：小车和摆杆，分别如图2和图3所示。图 2 小车图 3 摆杆由小车水平方向的受力情况可以得到下面等式： Error! No bookmark name given.Error! No bookmark name given.由摆杆水平方向的受力情况可以得到下面等式： 将代入中得：Error! No b