《飞行器系统仿真与CAD 》学习报告第一部分仿真(40)题目1：给定导弹相对于目标的运动学方程组为q k q V q V q r q V q V r m m ,sin )sin(),cos(cosr(0) = 5km, q(0) = 60deg, (0) = 30deg,V = , V m = , 1Ma = 340m/s, k = 2 （1） 建立系统的方框图模型； （2） 用MATLAB 语言编写S —函数 （3） 用窗口菜单对(1), (2)进行仿真，动态显示结果; （4）用命令行对(1), (2)进行仿真，以图形显示结果答： （1）（2）用MATLAB 语言编写S 函数function [sys,x0,str,ts]=CAD1_sfun(t,x,u,flag) switch flag case 0[sys,x0,str,ts]=mdlInitializeSizes; case 1sys = mdlDerivatives(t,x,u); case 3sys = mdlOutputs(t,x,u); case {2,4,9} sys = []; otherwiseerror('unhandled flag=',num2str(flag))endfunction [sys,x0,str,ts]=mdlInitializeSizes sizes=simsizes;=3;=0;=3;=0;=1;=1;sys=simsizes(sizes);str=[];x0=[5000,pi/3,pi/6];ts=[0 0];function sys=mdlDerivatives(t,x,u)vm=*340;v=*340;k=2;dx(1)=vm*cos(x(2))-v*cos(x(2)-x(3));dx(2)=(v*sin(x(2)-x(3))-vm*sin(x(2)))/x(1); dx(3)=k*dx(2);sys=dx;function sys=mdlOutputs(t,x,u)sys=x;调用S函数的模型框图(3)框图仿真结果：S函数仿真结果：（4）命令输入clear;clc[t x ] = sim('CAD1');hSimulink = figure();subplot(3, 1, 1);plot(t,x(:,1)); grid; ylabel('r'); subplot(3, 1, 2);plot(t,x(:,2)); grid; ylabel('q'); subplot(3, 1, 3);plot(t,x(:,3)); grid; ylabel('sigma'); [t x ] = sim('CAD1_S');hSFun = figure();subplot(3, 1, 1);plot(t,x(:,1)); grid; ylabel('r'); subplot(3, 1, 2);plot(t,x(:,2)); grid; ylabel('q'); subplot(3, 1, 3);plot(t,x(:,3)); grid; ylabel('sigma');模型仿真结果：S 函数仿真结果：题目2：给出动态方程0)0(,1)0(;1)1(2x x tx x x x ;(1) 用MATLAB语言编写S—函数;(2) 用命令行gear/adams法对(1)进行仿真，显示曲线x(t=0：100);(3) 建立方框图,用RK45仿真50秒,显示曲线答：（1）用MATLAB语言编写S—函数function[sys,x0,str,ts]=CAD2_sfu n(t,x,u,flag)switch flagcase 0[sys,x0,str,ts]=mdlInitializeSiz es;case 1sys=mdlDerivatives(t,x,u);case 3sys=mdlOutputs(t,x,u);case {2,4,9}sys=[];otherwiseerror('unhandledflag=',num2str(flag))endfunction[sys,x0,str,ts]=mdlInitializeSiz es sizes=simsizes;=2;=0;=2;=0;=1;=1;sys=simsizes(sizes);str=[];x0=[1,0];ts=[0 0];functionsys=mdlDerivatives(t,x,u)dx(1)=x(2);dx(2)=1-t*x(1)-(1-x(1)^2)*x(2); sys=dx;function sys=mdlOutputs(t,x,u) sys=x;（2）直接调用ode数值积分函数进行仿真，系统微分方程：function dx = CAD01_02odefun(t, x)dx(1) = x(2);dx(2) = 1-(1-x(1)*x(1))*x(2) - t*x(1);dx = dx';调用ode解算器入口：clear; clc;[t x] = ode15s(@CAD01_02odefun, 0:100, [1 0]);hGear = figure();set(hGear, 'NumberTitle', 'off', 'Name', 'Integrated by the Gear algorithm', 'Units', 'Normalized', 'Position', [ ]);subplot(2, 1, 1);plot(t, x(:,1)); grid; ylabel('x');subplot(2, 1, 2);plot(t, x(:,2)); grid; ylabel('dx/dt');[t x] = ode113(@CAD01_02odefun, 0:100, [1 0]);hAdams = figure();set(hAdams, 'NumberTitle', 'off', 'Name', 'Integrated by the Adams algorithm', 'Units', 'Normalized', 'Position', [ ]);subplot(2, 1, 1);plot(t, x(:,1)); grid; ylabel('x');subplot(2, 1, 2);plot(t, x(:,2)); grid; ylabel('dx/dt');ode15s(Gear) 仿真结果：ode113(Adams) 仿真结果：（3）建立方框图,用RK45仿真50秒,显示曲线方框图模型：仿真结果：问题3：质量—弹簧系统，质量M ，恢复系数K ，阻力系数C ，主动力P ，动力学方程为M x C x Mg sign x Kx P()()2M=1kg, K=4kg/s 2, C=100kg/m, g= s 2, =;;（1）在原点处用linmod 线性化，求线性系统的A,B,C,D; （2）对线性模型，判断能控性;（3）对线性模型，求阶跃、脉冲响应曲线;（4）对原模型进行仿真，P=sin(t) (使用Simulink); （5）对原模型进行仿真，P=sin(t) (使用ode23)答：（1）①线性化时需在模型中制定输入端、输出端（状态），如下图，状态选为位置和速度②linmod函数应用于该系统会出现奇异，故选用改进的linmod2 函数: clc;clear;[A,B,C,D]=linmod2('CAD3');ss0 = ss(A, B, C, D);Co = ctrb(ss0) ;[row col] = size(A);isControllable = ~(rank(Co, eps) - row);hStep = figure();set(hStep, 'NumberTitle', 'off', 'Name', 'StepResponse','unit','normalized','Position' ,[,,,]);step(ss0);grid;hImpulse = figure();set(hImpulse, 'NumberTitle', 'off', 'Name', 'ImpulseResponse','unit','normalized','Position' ,[,,,]);impulse(ss0);grid;命令窗口输出结果：A =+008 *B =1C =1 00 1D =The system is controlled （3）阶跃响应：脉冲响应:（4）对原模型进行仿真，P=sin(t) (使用Simulink)仿真结果：（5）对原模型进行仿真，P=sin(t) (使用ode23)系统微分方程：function dx = CAD3odefun(t, x)M = 1; K = 4; C = 100; g = ; miu = ;dx(1) = x(2);dx(2) = (sin(t)-K*x(1)-sign(x(2))*(C*x(2)*x(2)+miu*M*g))/M; dx = dx';仿真入口程：clc;clear;options = odeset('RelTol',1e-3,'AbsTol',[1e-5 5e-5]); [t x] = ode23(@CAD3odefun, 0::10, [0 0], options); hode23 = figure();set(hode23, 'NumberTitle', 'off', 'Name', 'Integrated by the ode23 solver',...'Units', 'Normalized', 'Position', [ ]); subplot(2, 1, 1);plot(t, x(:,1)); grid; ylabel('x'); subplot(2, 1, 2);plot(t, x(:,2)); grid; ylabel('dx/dt');仿真结果：问题4：给出一个系统, 要求生成一个新Simulink 模块，实现其功能(1)Mask 功能（2）s-函数yuy ly u u l u rK=1答：实现所需功能的S 函数 function [sys,x0,str,ts] =CAD01_04sfun_kernel(t,x,u,flag,ul,ur,yl ,yr)switch flag,case 0,[sys,x0,str,ts]=mdlInitializeSizes;case 3,sys=mdlOutputs(t,x,u,ul,ur,yl,yr);case 9, sys=[]; endfunction[sys,x0,str,ts]=mdlInitializeSizessizes = simsizes;= 0;= 0;= 1;= 1;= 1;= 1;sys = simsizes(sizes); x0 = [];str = [];ts = [0 0]; functionsys=mdlOutputs(t,x,u,ul,ur,yl,yr) if (u >= ur + yr)y = yr;elseif (u <= ul + yl)y = yl;elseif (u >ul + yl) && (u <ul)y = u - ul;elseif (u <ur + yr) && (u >ur)y = u - ur;elsey = 0;endsys = y;在Simulink中将调用S函数的模块进行封装参数传递及初始化用户界面：测试结果问问题5：已知系统A = [0 1; -1 -2], B = [1 0; 0 1], C = [1 0; 0 1], D = [0 0; 0 0], 求系统的状态空间方程(linmod)，并分析系统的稳定性，练习仿真参数设置答：对模型进行线性化并分析稳定性clear; clc;[A B C D] = linmod('CAD5')ss0 = ss(A, B, C, D);hpz = figure();set(hpz, 'NumberTitle', 'off', 'Name', 'Pole-zero map of the linmod system');pzmap(ss0);sgrid;[row col] = size(A);P = lyap(A, eye(row));for i = 1:rowsubdet(i) = det(P(1:i,1:i));endsubdet系统零极点图：存在正实部的极点，系统不稳定。