飞行器系统仿真.docx

1、飞行器系统仿真飞行器系统仿真与 CAD学习报告第一部分仿真(40)题目1给定导弹相对于目标的运动学方程组为r(0) = 5km, q(0) = 60deg, ?(0) = 30deg,V = , V , 1Ma = 340m/s, k = 2(1)建立系统的方框图模型；(2)用MATLAB语言编写S函数(3)用窗口菜单对(1), (2进行仿真，动态显示结果(4)用命令行对(1), (2进行仿真，以图形显示结果答:(2)用MATLAB语言编写S函数function sys,x0,str,ts=CAD1_sfun(t,x,u,flag) switch flagcase0sys,x0,str,ts=

2、mdlInitializeSizes;case1sys = mdlDerivatives(t,x,u);case3sys = mdlOutputs(t,x,u);case2,4,9sys = ;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;functionsys=mdlDerivatives(t,x,

3、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;clct x = sim(CAD1); hSimulink = figure(); subplot(3, 1, 1);plot(t,x(:,1); grid; ylabelr();subplot(3, 1

4、, 2);plot(t,x(:,2); grid; ylabelq( );subplot(3, 1, 3);plot(t,x(:,3); grid; ylabels(igma);t x = sim(CAD1_S);hSFun = figure();subplot(3, 1, 1);plot(t,x(:,1); grid; ylabelr();subplot(3, 1, 2);plot(t,x(:,2); grid; ylabelq( );subplot(3, 1, 3);plot(t,x(:,3); grid; ylabels(igma);模型仿真结果：S函数仿真结果：题目 2：给出动态方程

5、x (1 x 2 ) x tx 1; x (0) 1, x (0) 0 用MATLAB语言编写S-函数;(2)用命令行gear/adams法对(1)进行仿真，显示曲线 x(t=O: 100); 建立方框图，用RK45仿真50秒显示曲线答:(1)用MATLAB语言编写S函数functionsys,x0,str,ts=CAD2_sfun(t,x,u,fl otherwiseag)error(unhandled flag=,num2str(flag)switch flagendcase0function sys,x0,str,ts=mdlInitializeSizessizes=simsizes;s

6、ys,x0,str,ts=mdlInitializeSizes;=2;case1=0;sys=mdlDerivatives(t,x,u);=2;case3=0;sys=mdlOutputs(t,x,u);=1;case2,4,9 sys=;dx(1)=x(2);sys=simsizes(sizes);ts=0 0;function sys=mdlOutputs(t,x,u)sys=x;function sys=mdlDerivatives(t,x,u)(2)直接调用 ode 数值积分函数进行仿真，系统微分方程：function dx = CAD01_02odefun(t, x)dx(1) =

7、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; ylabelx();subplot(2, 1, 2);plo

8、t(t, x(:,2); grid; ylabeld( x/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; ylabelx();subplot(2, 1, 2);plot(t, x(:,2); grid; ylabeld( x/dt );o

9、de15s(Gear)仿真结果：ode113(Adams)仿真结果：( 3)建立方框图 ,用 RK45 仿真 50 秒,显示曲线 方框图模型：仿真结果:问题3:质量一弹簧系统，质量M，恢复系数K，阻力系数C,主动力P,动力学方程为 Mx (Cx2 Mg)sign(x) Kx PM=1kg, K=4kg/s 2, C=100kg/m,2 og= s, ?=；(1)在原点处用linmod线性化，求线性系统的A,B,C,D;(2)对线性模型，判断能控性；(3)对线性模型，求阶跃、脉冲响应曲线；(4)对原模型进行仿真，P=sin(t)(使用Simulink);(5)对原模型进行仿真，P=sin(t)(

10、使用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, Step Respons,eunit,normalized,Positi

11、on ,); step(ss0);grid;hImpulse = figure();set(hImpulse,NumberTitle, off, Name, ImpulseResponse,unit,normalized,Position ,);impulse(ss0);grid;命令窗口输出结果：A =+008B =01C =1 00 1D =00The system is controlled3）阶跃响应：脉冲响应 :( 4)对原模型进行仿真， P=sin(t) (使用 Simulink) 仿真结果：(5)对原模型进行仿真， P=sin(t) (使用 ode23) 系统微分方程：funct

12、ion 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,

13、 Integrated by the ode23 solve,.r.Units, Normalized, Position, );subplot(2, 1, 1);plot(t, x(:,1); grid; ylabelx();subplot(2, 1, 2);plot(t, x(:,2); grid; ylabeld( x/dt );仿真结果：问题4：给出一个系统，要求生成一个新Simulink模块，实现其功能(1)Mask功 能(2)s-函数答：实现所需功能的S函数function sys,x0,str,ts = sys,x0,str,ts=mdlInitializeSizes;CAD01

14、_04sfun_kernel(t,x,u,flag,ul,ur,yl,ycase3,r)sys=mdlOutputs(t,x,u,ul,ur,yl,yr);switch flag,case9,case0,sys=;function sys,x0,str,ts=mdlInitializeSizes if (u = ur + yr)sizes = simsizes;= 0;= 0;= 1;= 1;= 1;= 1;sys = simsizes(sizes);x0 = ;str = ;ts = 0 0;在 Simulink 中将调用 Sy = yr;elseif (u ul + yl) & (u ul

15、) y = u - ul;elseif (u ur) y = u - ur;elsey = 0;endsys = y;参数传递及初始化用户界面：测试结果问题 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 = linm od(CAD5)ssO = ss(A, B, C, D);hpz = figure();set(hp z, NumberTitle, off, Name,卩

16、ole-zero map of the linmod systempzmap(ss0);sgrid;row col = size(A);P = lyap(A, eye(row);for i = 1:rowsubdet(i) = det(P(1:i,1:i);endsubdet系统零极点图：存在正实部的极点，系统不稳定。问题 6：系统的动力学方程为 dx / dt = Ax + Bu, y = Cx + Du, A = 0 1 0 0; 0 0 1 0;0 0 0 1; -1 -2 -3 -4, B = 1 2 ; 3 4; 2 3; 4 5, C = 1 1 2 2; 2 3 5 4; D =

17、 1 0; 0 1求,：(1) 系统动态平衡点(2) x(0)=1 1 1 1 , ix=1 2 3 4 ,dx=O 1 0 1 的系统动态 2 3 4平衡点答：系统框图模型 系统的平衡点分析 程序clear; clc;x, u, y, dx, options = trim(CAD6 );options(10)x0=1 1 1 1; ix=1 2 3 4;dx=0 1 0 1;idx=1 2 3 4;x, u, y, dx, options = trim(CAD01_06, x0, , ,ix, , , dx, idx);options(10)运行结果x=0;0;0;0;u=0;0;y=0;0

18、;ans=9x=;u=;y=;ans=41问题 7：自学文件 C 与 M -s 函数模板和示例文件答：Simulink中的示例文件实现了将输入信号放大为 2倍输出的功能，自学时对示例程序进行改进，使之可以指定信号放大的倍数。语言S函数源代码#define S_FUNCTION_NAME CAD02_07sfun/*Modified: change the function name*/#define S_FUNCTION_LEVEL 2#include static void mdlInitializeSizes(SimStruct *S)ssSetNumSFcnParams(S, 1);/*

19、Revised: set the number of input parameters to 1*/if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S) return;if (!ssSetNumInputPorts(S, 1) return;ssSetInputPortWidth(S, 0, DYNAMICALLY_SIZED); ssSetInputPortDirectFeedThrough(S, 0, 1);if (!ssSetNumOutputPorts(S,1) return; ssSetOutputPortWidth(S, 0, D

20、YNAMICALLY_SIZED);ssSetNumSampleTimes(S, 1);ssSetOptions(S, SS_OPTION_EXCEPTION_FREE_CODE);static void mdlInitializeSampleTimes(SimStruct *S)ssSetSampleTime(S, 0, INHERITED_SAMPLE_TIME); ssSetOffsetTime(S, 0, ;static void mdlOutputs(SimStruct *S, int_Ttid)int_T i;InputRealPtrsTypeuPtrs = ssGetInputP

21、ortRealSignalPtrs(S,0);real_T *y = ssGetOutputPortRealSignal(S,0);int_T width = ssGetOutputPortWidth(S,0);constmxArray *pmxRatio = ssGetSFcnParam(S,0);/*Revised: get thepointer to the parameter in the type of mxArray*/constreal_T *pRatio = mxGetPr(pmxRatio);/*Revised: get the pointer tothe parameter

22、 in the type of real_T*/for (i=0; i0varargout: = sf( varargin: );elsesf(varargin:);endStateflow是有限状态机(finite state machine的图形工具，它可以用于解决 复杂的逻辑问题，用户可以通过图形化工具实现在不同状态之间的转换。Stateflow/可以直接嵌入到Simulink仿真模型中，并且在仿真的初始化阶段，SIMULINK会把Stateflow绘制的逻辑图形通过编译程序转换成 C语言，使二者有机 地结合在一起。Stateflow可以在SIMULINK Extra模块库中找到。Sta

23、teflow的仿真原理是有限状态机(finite state machine理论，有限状态机 是指系统含有可数的状态， 在相应的状态事件发生时， 系统会从当前状态转移到 与之对应的状态。 在有限状态机中实现状态的转移是有一定条件的， 同时相互转 换的状态都会有状态转移事件，这样就构成了状态转移图。在 SIMULINK 的仿真 窗口中，允许用户建立有限个状态以及状态转移的条件与事件， 从而绘制出有限状态机系统，这样就可以实现对系统的仿真。Stateflow的仿真框图一般都会嵌入 到Simuli nk仿真模型中，同时实现状态转移的条件或是事件即可以取自 Stateflow仿真框图，也可以来自Sim

24、ulink仿真模型。第二部分综合（ 60）2.J2HPOP3 采用 STK 轨道机动模块实现任意两个圆轨道间的霍曼转移答：利用轨道转移模块 Astrogator,结合轨道机动打靶功能Target实现霍曼 转移。对卫星的Astrogator设置如下的任务控制序列（MS609 May 2013 15:17:05Satellite-Satellite1Astrogator Mission Control Sequence SummaryMCS Segment Type: InitialStateName: Initial StateUser Comment: -Initial State Descr

25、iption-Satellite State at End of Segment:UTC Gregorian Date: 1 Jun 2003 12:00:UTC Julian Date: 2452792Julian Ephemeris Date:Time past epoch: 0 sec (Epoch in UTC Gregorian Date: 1 Jun 2003 12:00:State Vector in Coordinate System: Earth Centered Mean J2000Parameter Set Type: CartesianParameter Set Typ

26、e: Keplerianecc:w:0 deginc:0 degTA:0 degParameter Set Type: SphericalDecl:0 degAzimuth:90 degOther Elliptic Orbit Parameters :Ecc. Anom:0 degMean Anom:0 degLong Peri:0 degArg. Lat:0 degTrue Long:0 degVert FPA:90 degVel. RA:90 degVel. Decl:0 degTime Past Periapsis:0 secGeodetic Parameters:Geocentric

27、Parameters:Spacecraft Configuration:Drag Area: 1e-006 kmA2SRP Area: 1e-006 kmA2Dry Mass: 500 kgFuel Mass: 500 kgTotal Mass: 1000 kgArea/Mass Ratio: 1e-009 kmA2/kgPaTank Pressure: 5000Cr:Cd:MCS Segment Type: PropagateName: PropagateUser Comment: -Propagation Description-Propagator model used: Earth_Point_Mass(Simple numerical twobody)Stopping Condition Information (Gregorian Date Julian Date):Propagation Statistics:Number of steps: 581Average step size: secLargest step size: secSmallest step size: secSatellite State at End of Segment:UTC Gregorian Date: 2 Jun 2003 00:00:UTC Julian Date: