机电系统动力学仿真matlab课后答案刘白雁

二、仿真的分类
按模型分类
1、物理仿真：采用物理模型，有实物介入！ 具有效果逼真，精度高等优点，但造价高或耗时长， 大多在一些特殊场合下采用（如导弹、卫星一类飞 行器的动态仿真，发电站综合调度仿真与培训系统 等），具有实时性、在线的特点。 2、数学仿真：采用数学模型 在计算机上进行，具有非实时性、离线的特点，经 济、快速、实用。
《机电系统动态仿真——基于 MATLAB/Simulink》
刘白雁教授编著 机械工业出版社
2006，8
仿真软件的简介 一、仿真的发展
1、程序编程阶段： 所有问题（如：微分方程求解、 矩阵运算、绘图等）都是用高级算法语言（如C、 FORTRAN等）来编写。 2、程序软件包阶段： 出现了“应用子程序库”。 3、交互式语言阶段（仿真语言：仿真语言可用一 条指令实现某种功能，如“系统特征值的求解”， 使用人员不必考虑什么算法，以及如何实现等低 级问题。 4、模型化图形组态阶段：符合设计人员对基于模 型图形化的描述。
三、常见的几种仿真软件
PSPICE、ORCAD：通用的电子电路仿真软件， 适合于元件级仿真。 SYSTEM VIEW：系统级的电路动态仿真软 件 MATLAB：具有强大的数值计算能力，包含 各种工具箱，其程序不能脱离MATLAB环境 而运行，所以严格讲，MATLAB不是一种计 算机语言，而是一种高级的科学分析与计算软 件。 SIMULINK：是MATLAB附带的基于模型化 图形组态的动态仿真环境。
按计算机类型分类
1、模拟仿真：采用数学模型，在模拟计算机 上进行的实验研究。50年代 2、数字仿真：采用数学模型，在数字计算机 上借助于数值计算方法所进行的仿真实验。 60年代
3、混合仿真：结合了模拟仿真与数字仿真。 4、现代计算机仿真：采用先进的微型计算机，基于 专用的仿真软件、仿真语言来实现，其数值计算 功能强大，使用方便，易学。80年代以来

实验一MATLAB基本操作一、实验目的：①通过上机实验操作，使学生熟悉MATLAB实验环境，练习MATLAB命令、m文件，进行矩阵运算、图形绘制、数据处理。

②通过上机操作，使得学生掌握Matlab变量的定义和特殊变量的含义，理解矩阵运算和数组运算的定义和规则。

③通过上机操作，使得学生掌握数据和函数的可视化，以及二维曲线、三维曲线、三维曲面的各种绘图指令。

二、实验原理与说明Matlab是Matrix 和Laboratory两词的缩写，是美国Mathworks公司推出的用于科学计算和图形处理的可编程软件，经历了基于DOS版和Windows版两个发展阶段。

三、实验设备与仪器：PC电脑，Matlab7.0仿真软件四、实验内容、方法与步骤：数组运算与矩阵运算数组“除、乘方、转置”运算符前的“.”决不能省略，否则将按矩阵运算规则进行运算；执行数组与数组之间的运算时，参与运算的数组必须同维，运算所得的结果也与参与运算的数组同维。

A=[ 1 2 3; 4 5 6; 7 8 9];B=[-1 -2 -3；-4 -5 -6；-7 -8 -9];X=A.*BY=A*Bplot用于二维曲线绘图，若格式为plot（X，Y，’s’），其中X为列向量，Y是与X等行的矩阵时，以X为横坐标，按Y的列数绘制多条曲线；若X为矩阵，Y是向量时，以Y为纵坐标按X的列数（或行数）绘制多条曲线。

参考程序如下：t=(0:pi/100:pi)'y1=sin(t)*[-1 1];y2=sin(t).*sin(9*t);plot(t,y1, 'r:', t, y2, 'b-.')axis([0 pi, -1, 1])title('Drawn by Dong-yuan GE')程序运行界面如下：plot3用于三维曲线绘制，其使用格式与plot十分相似。

参考程序如下：t=0:0.02:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z,'b-', x,y,z,'o')程序运行界面如下：mesh与surf用于三维空间网线与曲面的绘制。

控制系统仿真第六章课后题作业6.1在图6.1中，已知单位负反馈系统被控对象的传递函数为)1001.0)(11.0()(++=s s s K s G 试编写matlab 程序，设计系统的超前矫正器Gc(s),要求：1）在斜坡信号r （t)=2t 作用下，系统的稳态误差ess<=0.002;2)校正后系统的相位裕度Pm 范围为：45~55；3）绘制系统校正后的bode 图和阶跃响应曲线。

程序：>> s=tf('s');>> G=1000/(s*(0.1*s+1)*(0.001*s+1));>> margin(G) % 绘制校正前的bode 图>> figure(2)>> sys=feedback(G,1);>> step(sys) %绘制校正前的单位阶跃响应曲线>> [Gm,Pm]=margin(G); %该句值计算bode图的增益裕量Gm和相位裕量Pm >> [mag,phase,w]=bode(G); %该句只计算bode图上多个频率点w对应的幅值和相位>> QWPm=50; %取矫正后的相位为50>> FIm=QWPm-Pm+5;>> FIm=FIm*pi/180;>> alfa=(1-sin(FIm))/(1+sin(FIm));>> adb=20*log10(mag);>> am=10*log10(alfa);>> wc=spline(adb,w,am);>> T=1/(wc*sqrt(alfa));>> alfat=alfa*T;>> Gc=tf([T 1],[alfat 1]) %校正器的传递函数Transfer function:0.01794 s + 1-------------0.00179 s + 1>> figure(3)>> margin(Gc*G) % 系统矫正后的bode图>> figure(4)>> step(feedback(Gc*G,1)) % 校正后的单位阶跃响应曲线作业6.2在图6.1中，已知单位负反馈系统被控对象的传递函数为)102.0)(11.0()(++=s s s K s G 试编写matlab 程序，设计系统的滞后校正器Gc （s ）,要求：1) 在斜坡信号r(t)=t 作用下，系统的稳态误差ess 01.0≤；2) 校正后系统的相位裕度Pm 范围为：40~50；3) 绘制系统矫正前后的bode 图和阶跃响应曲线。

【4.2】程序：num=[5,0];den=conv([1,1],conv([1,2],[1,3])); [numc,denc]=cloop(num,den);[z,p,k]=tf2zp(numc,denc);[A,B,C,D]=tf2ss(numc,denc);g_zp=zpk(z,p,k)g_tf=tf(numc,denc)g_ss=ss(A,B,C,D)运行结果：Zero/pole/gain:5 s----------------------------------(s+0.4432) (s^2 + 5.557s + 13.54)Transfer function:5 s----------------------s^3 + 6 s^2 + 16 s + 6a =x1 x2 x3x1 -6 -16 -6x2 1 0 0x3 0 1 0b =u1x1 1x2 0x3 0c =x1 x2 x3y1 0 5 0d =u1y1 0【4.3】程序：A=[0 0 0 -1;1 0 0 -2;0 1 0 -3;0 0 1 -4]; B=[0;0;0;1];C=[1 0 0 0];g_ss=ss(A,B,C,D)[num,den]=ss2tf(A,B,C,D);g_tf=tf(num,den)[z,p,k]=ss2zp(A,B,C,D);g_zpk=zpk(z,p,k)运行结果：a =x1 x2 x3 x4x1 0 0 0 -1x2 1 0 0 -2x3 0 1 0 -3x4 0 0 1 -4b =u1x1 0x2 0x3 0x4 1c =x1 x2 x3 x4y1 1 0 0 0d =u1y1 0Continuous-time model.Transfer function:-3.109e-015 s^3 - s^2 - 3.331e-015 s - 4.441e-016 -------------------------------------------------s^4 + 4 s^3 + 3 s^2 + 2 s + 1Zero/pole/gain:- s^2----------------------------------------------(s+0.6724) (s+3.234) (s^2 + 0.0936s + 0.4599)【5.1】（1）程序num=[0,10];den=conv([1,0],[1,7,17]); [numc,denc]=cloop(num,den,-1); G=tf(numc,denc)[y,t]=step(G);plot(t,y,'b-')C=dcgain(G);n=1;while y(n)<0.1*Cn=n+1;endm=1;while y(m)<0.9*Cm=m+1;endrisetime=t(m)-t(n)[Y,k]=max(y); percentovershoot=100*(Y-C)/Ci=length(t);while(y(i)>0.98*C)&(y(i)<1.02*C) i=i-1;endsettlingtime=t(i)运行结果：Transfer function:10-----------------------s^3 + 7 s^2 + 17 s + 10risetime =2.7312percentovershoot =-0.4399settlingtime =5.1372图：0123456700.10.20.30.40.50.60.70.80.91（2）程序k=[10,100,1000];t=linspace(1,20,200);num=1;den=conv([1,0],[1,7,17]);for j=1:3;s1=tf(num*k(j),den);sys=feedback(s1,1)y(:,j)=step(sys,t);endplot(t,y(:,1),'r',t,y(:,2),'b',t,y(:,3),'g')gtext('k=10');gtext('k=100');gtext('k=1000') 运行结果：Transfer function:10-----------------------s^3 + 7 s^2 + 17 s + 10Transfer function:100------------------------s^3 + 7 s^2 + 17 s + 100Transfer function:1000-------------------------s^3 + 7 s^2 + 17 s + 1000图：024681012141618200.20.40.60.811.21.41.61.8图：02468101214161820-3-2-1123422【6.1】程序：（1）num1=[1,1];den1=conv([1,0,0],conv([1,2],[1,4]));sys1=tf(num1,den1)rlocus(sys1)运行结果：-12-10-8-6-4-2024-8-6-4-202468Root LocusReal Axis I m a g i n a r y A x i s（2）num2=[1,1];den2=conv([1,0],conv([1,-1],[1,4,16]));sys2=tf(num2,den2)rlocus(sys2)运行结果：-10-8-6-4-2024-8-6-4-202468Root LocusReal Axis I m a g i n a r y A x i s（3）num3=[1,8];den3=conv([1,0,0],conv([1,3],conv([1,5],conv([1,7],[1,15])))); sys3=tf(num3,den3)rlocus(sys3)运行结果：-30-25-20-15-10-5051015-20-15-10-505101520Root LocusReal Axis I m a g i n a r y A x i s【6.3】程序：num=[1,2];den=conv([1,0],conv([1,4],conv([1,8],[1,2,5])));sys=tf(num,den)rlocus(sys)[k,poles]=rlocfind(sys)运行结果：Transfer function:s + 2---------------------------------------s^5 + 14 s^4 + 61 s^3 + 124 s^2 + 160 sSelect a point in the graphics windowselected_point =0.0296 + 2.2826i k =135.8815poles =-7.3248-5.41040.0145 + 2.3021i0.0145 - 2.3021i -1.2939图：-20-15-10-5051015-15-10-551015Root LocusReal Axis I m a g i n a r y A x i s【7.3】程序（1）画波特图num=[50];den=conv([1,0],conv([1,10],[3,1]));sys=tf(num,den)sys1=feedback(sys,1)bode(sys)grid图（1）-150-100-50050100M a g n i t u d e (d B)10-210-1100101102103-270-225-180-135-90P h a s e (d e g )Bode DiagramFrequency (rad/sec)程序（2）画奈奎斯特图num=[50];den=conv([1,0],conv([1,10],[3,1]));sys=tf(num,den)sys1=feedback(sys,1)nyquist(sys)grid图（2）-16-14-12-10-8-6-4-20-300-200-100100200300Nyquist DiagramReal Axis I m a g i n a r y A x i s程序（3）画零极点图num=[50];den=conv([1,0],conv([1,10],[3,1]));sys=tf(num,den)sys1=feedback(sys,1)pzmap(sys1)gird图（3）P ole-Zero MapReal Axis I m a g i n a r y A x i s -12-10-8-6-4-20-1.5-1-0.50.511.5程序（4）计算相角裕量和幅值裕量num=[50];den=conv([1,0],conv([1,10],[3,1]));sys=tf(num,den)sys1=feedback(sys,1)[gm,pm,wcg,wcp]=margin(sys)运行结果Transfer function:50---------------------3 s^3 + 31 s^2 + 10 sTransfer function:50--------------------------3 s^3 + 31 s^2 + 10 s + 50gm =2.0667pm =7.5615wcg =1.8257wcp =1.2645程序（5）绘制阶跃响应曲线num=[50];den=conv([1,0],conv([1,10],[3,1])); sys=tf(num,den)sys1=feedback(sys,1)step(sys1)图（5）00.20.40.60.811.21.41.61.82Step ResponseTime (sec)A m p l i t u d e【7.4】程序如下：num=[300];den=conv([1,0,0],conv([0.2,1],[0.02,1]));sys=tf(num,den)margin(sys)grid波特图如下：-150-100-50050100M a g n i t u d e (d B )10-1100101102103-360-315-270-225-180P h a s e (d e g )Bode DiagramGm = Inf , P m = -78 deg (at 11 rad/sec)Frequency (rad/sec)【9.3】程序：A=[-2 2 -1;0 -2 0;1 -4 0];B=[0;0;1];C=[1,0,0];D=0;M=ctrb(A,B)m=rank(M)if m==3;disp('系统可控')elsedisp('系统不可控')endN=obsv(A,C)n=rank(N)if n==3;disp('系统可观')elsedisp('系统不可观') endsys=ss(A,B,C,D) [num,den]=ss2tf(A,B,C,D) sys1=tf(num,den)[z,p,k]=ss2zp(A,B,C,D)运行结果：M =0 -1 20 0 01 0 -1m =2系统不可控N =1 0 0-2 2 -13 -4 2n =2系统不可观a =x1 x2 x3x1 -2 2 -1x2 0 -2 0x3 1 -4 0b =u1x1 0x2 0x3 1c =x1 x2 x3y1 1 0 0d =u1y1 0 Continuous-time model.0 0 -1 -2den =1 4 5 2Transfer function:-s - 2---------------------s^3 + 4 s^2 + 5 s + 2z =-2p =-1-1-2k = -1【10.1】（1）程序：A=[0,1,0,0;0,5,0,0;0,0,-7,0;0,0,0,-8]; B=[0;1;0;1];C=[1,2,3,4];D=zeros(1,1);G_ss=ss(A,B,C,D)运行结果：a =x1 x2 x3 x4x1 0 1 0 0x2 0 5 0 0x3 0 0 -7 0x4 0 0 0 -8u1x1 0x2 1x3 0x4 1c =x1 x2 x3 x4y1 1 2 3 4d =u1y1 0（2）：程序：[num1,den1]=ss2tf(A,B,C,D); p=roots(den1)i=0;for k=1:1:length(p)if real(p(k))>0i=i+1;endendif i>0disp('系统不稳定');elsedisp('系统稳定');end运行结果：p =5.0000-8.0000-7.0000系统不稳定（3）（4）程序：AA=[0,1,0;0,5,0;0,0,-8];BB=[0;1;1];P=[-1,-2,-8];K=acker(AA,BB,P);i=4;K(4)=0;Kpp=eig(A-B*K)sys1=tf(num1,den1);[y1,t]=step(sys1);plot(t,y1)hold onA_feedback=A-B*K;[num2,den2]=ss2tf(A_feedback,B,C,D); sys2=tf(num2,den2);[y2,t]=step(sys2);plot(t,y2,'r')gridgtext('反馈前')gtext('反馈后')运行结果：K =2 8 0 0pp =-8-2-1-7图形：01234560123456【13.1】程序：A=[0,1;0,0];B=[0;1];C=[1,0];D=zeros(1,1);G_ss=ss(A,B,C,D)M=ctrb(A,B);if rank(M)==2disp('系统完全能控'); elsedisp('系统不完全能控'); endS=[1,0];N=obsv(A,S);if rank(N)==2disp('(A,S)可观测'); elsedisp('(A,S)不可观测'); endR=1;Q=[1,0;0,0];[K,P,E]=Lqr(A,B,Q,R)A_new=A-B*K;G_new=ss(A_new,B,C,D);t=linspace(0,5,100)';y1=step(G_ss,t);y2=step(G_new,t);plot(t,y1,'r:',t,y2,'b-')gridgtext('反馈前')gtext('反馈后')运行结果：a =x1 x2x1 0 1x2 0 0b =u1x1 0x2 1c =x1 x2y1 1 0d =u1y1 0Continuous-time model. 系统完全能控(A,S)可观测K =1.0000 1.4142P =1.4142 1.00001.0000 1.4142E =-0.7071 + 0.7071i-0.7071 - 0.7071i图形：00.51 1.52 2.53 3.54 4.5502468101214。