飞行器系统仿真.doc

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

2、=CAD1_sfun(t,x,u,flag)switch flagcase 0sys,x0,str,ts=mdlInitializeSizes;case 1sys = mdlDerivatives(t,x,u);case 3sys = mdlOutputs(t,x,u);case 2,4,9sys = ;otherwiseerror(unhandled flag=,num2str(flag)endfunction sys,x0,str,ts=mdlInitializeSizessizes=simsizes;sizes.NumContStates=3;sizes.NumDiscStates=0;

3、sizes.NumOutputs=3;sizes.NumInputs=0;sizes.DirFeedthrough=1;sizes.NumSampleTimes=1;sys=simsizes(sizes);str=;x0=5000,pi/3,pi/6;ts=0 0;functionsys=mdlDerivatives(t,x,u)vm=0.9*340;v=1.2*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=m

4、dlOutputs(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; 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(); s

5、ubplot(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)(,1)0(;1)1( xxtxx(1) 用 MATLAB 语言编写 S函数;(2) 用命令行 gear/adams 法对(1)进行仿真，显示曲线x(t=0：100);(3) 建立方框图,用 RK45 仿真 50 秒,显示曲线答：（1）

6、用 MATLAB 语言编写 S函数functionsys,x0,str,ts=CAD2_sfun(t,x,u,flag)switch flagcase 0sys,x0,str,ts=mdlInitializeSizes;case 1sys=mdlDerivatives(t,x,u);case 3sys=mdlOutputs(t,x,u);case 2,4,9sys=;otherwiseerror(unhandled flag=,num2str(flag)endfunction sys,x0,str,ts=mdlInitializeSizessizes=simsizes;sizes.NumCon

7、tStates=2;sizes.NumDiscStates=0;sizes.NumOutputs=2;sizes.NumInputs=0;sizes.DirFeedthrough=1;sizes.NumSampleTimes=1;sys=simsizes(sizes);str=;x0=1,0;ts=0 0;function sys=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 数值积分函数进行仿真，系统微分方

8、程 CAD2odefun.m：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 解算器入口 CAD2ode_call.m：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, Pos

9、ition, 0.1 0.4 0.4 0.4); 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, 0.55 0.4

10、 0.4 0.4); 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 秒,显示曲线方框图模型 CAD2.mdl：仿真结果：问题 3：质量弹簧系统，质量 M，恢复系数 K，阻力系数 C，主动力 P，动力学方程为 M=1kg, xCgsinxP ()(2K=4kg/s2, C=100kg/m, g=9.8m/ s2, =0

11、.1;（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 = si

12、ze(A); isControllable = (rank(Co, eps) - row); hStep = figure(); set(hStep, NumberTitle, off, Name, Step Response,unit,normalized,Position ,0.05,0.5,0.4,0.4); step(ss0);grid; hImpulse = figure(); set(hImpulse, NumberTitle, off, Name, Impulse Response,unit,normalized,Position ,0.5,0.5,0.4,0.4); impul

13、se(ss0);grid; 命令窗口输出结果：A =1.0e+008 *0 0.0000-0.0000 -1.3286B =01C =1 00 1D =00The 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 = 9.8; miu = 0.1; dx(1) = x(2); dx(2) = (sin(t)-K*x

14、(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:0.05:10, 0 0, options); hode23 = figure(); set(hode23, NumberTitle, off, Name, Integrated by the ode23 solver,.Units, Normalized, Position, 0.25 0.4 0.5 0.5); subp

15、lot(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-函数 y uylyuul urK=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=mdlI

16、nitializeSizes; case 3, sys=mdlOutputs(t,x,u,ul,ur,yl,yr); case 9, sys=; endfunction sys,x0,str,ts=mdlInitializeSizessizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 0; sizes.NumOutputs = 1; sizes.NumInputs = 1; sizes.DirFeedthrough= 1; sizes.NumSampleTimes = 1; sys = simsizes(sizes)

17、; x0 = ; str = ; ts = 0 0; function sys=mdlOutputs(t,x,u,ul,ur,yl,yr) if (u = ur + yr) y = yr; elseif (u ul + yl) 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)，并分析系统的稳定性，练习仿真参数设置答：对模型进行线性化并分析稳

18、定性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:row subdet(i) = det(P(1:i,1:i); endsubdet系统零极点图：存在正实部的极点，系统不稳定。问题 6：系统的动力学方程为 dx / dt

19、= 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 = 1 0; 0 1, 求：(1)系统动态平衡点(2)x(0)=1 1 1 1, ix=1 2 3 4,dx=0 1 0 1,idx=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;

20、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;ans=9x=-0.0909;-0.9091;-0.0909;-0.0909;u=2.9091;-1.9091;y=3.3636;-0.1818;ans=41问题 7：自学文件 C 与 M -s 函数模板和示例文件答：Simulink 中的示例文件实现了将输入信号放大为 2 倍输出的功能，自学时对示例程序进行改进，使之可以指定信号放大的倍数。语言 S

21、函数源代码#define S_FUNCTION_NAME CAD02_07sfun /*Modified: change the function name*/ #define S_FUNCTION_LEVEL 2 #include “simstruc.h“ static void mdlInitializeSizes(SimStruct *S) ssSetNumSFcnParams(S, 1);/*Revised: set the number of input parameters to 1*/ if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCou

22、nt(S) return; if (!ssSetNumInputPorts(S, 1) return; ssSetInputPortWidth(S, 0, DYNAMICALLY_SIZED); ssSetInputPortDirectFeedThrough(S, 0, 1); if (!ssSetNumOutputPorts(S,1) return; ssSetOutputPortWidth(S, 0, DYNAMICALLY_SIZED); ssSetNumSampleTimes(S, 1); ssSetOptions(S, SS_OPTION_EXCEPTION_FREE_CODE);

23、static void mdlInitializeSampleTimes(SimStruct *S) ssSetSampleTime(S, 0, INHERITED_SAMPLE_TIME); ssSetOffsetTime(S, 0, 0.0); static void mdlOutputs(SimStruct *S, int_Ttid) int_T i; InputRealPtrsTypeuPtrs = ssGetInputPortRealSignalPtrs(S,0); real_T *y = ssGetOutputPortRealSignal(S,0); int_T width = s

24、sGetOutputPortWidth(S,0); constmxArray *pmxRatio = ssGetSFcnParam(S,0);/*Revised: get the pointer to the parameter in the type of mxArray*/ constreal_T *pRatio = mxGetPr(pmxRatio);/*Revised: get the pointer to the parameter in the type of real_T*/ for (i=0; i0varargout: = sf( varargin: );elsesf(vara

25、rgin:);endStateflow是有限状态机（finite state machine）的图形工具，它可以用于解决复杂的逻辑问题，用户可以通过图形化工具实现在不同状态之间的转换。Stateflow可以直接嵌入到Simulink仿真模型中，并且在仿真的初始化阶段，SIMULINK会把Stateflow绘制的逻辑图形通过编译程序转换成C语言，使二者有机地结合在一起。Stateflow可以在SIMULINK Extra模块库中找到。Stateflow的仿真原理是有限状态机（finite state machine）理论，有限状态机是指系统含有可数的状态，在相应的状态事件发生时，系统会从当前状态

26、转移到与之对应的状态。在有限状态机中实现状态的转移是有一定条件的，同时相互转换的状态都会有状态转移事件，这样就构成了状态转移图。在SIMULINK的仿真窗口中，允许用户建立有限个状态以及状态转移的条件与事件，从而绘制出有限状态机系统，这样就可以实现对系统的仿真。Stateflow的仿真框图一般都会嵌入到Simulink仿真模型中，同时实现状态转移的条件或是事件即可以取自Stateflow仿真框图，也可以来自Simulink仿真模型。第二部分综合（60）2.J2HPOP3 采用 STK 轨道机动模块实现任意两个圆轨道间的霍曼转移答：利用轨道转移模块 Astrogator，结合轨道机动打靶功能 T

27、arget 实现霍曼转移。对卫星的 Astrogator 设置如下的任务控制序列 (MSC)。09 May 2013 15:17:05Satellite-Satellite1Astrogator Mission Control Sequence Summary -*-* MCS Segment Type: InitialState Name: Initial State User Comment: -Initial State Description- - Satellite State at End of Segment: - UTC Gregorian Date: 1 Jun 2003 1

28、2:00:00.00 UTC Julian Date: 2452792 Julian Ephemeris Date: 2452792.00074287 Time past epoch: 0 sec (Epoch in UTC Gregorian Date: 1 Jun 2003 12:00:00.00) State Vector in Coordinate System: Earth Centered Mean J2000 Parameter Set Type: Cartesian X: 6700.0000000000000000 km Vx: 0.0000000000000000 km/se

29、c Y: 0.0000000000000000 km Vy: 7.7131448355214580 km/sec Z: 0.0000000000000000 km Vz: 0.0000000000000000 km/sec Parameter Set Type: Keplerian sma: 6700.0000000000000000 km RAAN: 0 deg ecc: 0.0000000000000000 w: 0 deg inc: 0 deg TA: 0 deg Parameter Set Type: Spherical Right Asc: 0 deg Horiz. FPA: -3.

30、509791787161889e-015 deg Decl: 0 deg Azimuth: 90 deg |R|: 6700.0000000000000000 km |V|: 7.7131448355214580 km/sec Other Elliptic Orbit Parameters : Ecc. Anom: 0 deg Mean Anom: 0 deg Long Peri: 0 deg Arg. Lat: 0 deg True Long: 0 deg Vert FPA: 90 deg Ang. Mom: 51678.07039799377 km2/sec p: 6699.9999999

31、999991000 km C3: -59.49260325373134 km2/sec2 Energy: -29.74630162686567 km2/sec2 Vel. RA: 90 deg Vel. Decl: 0 deg Rad. Peri: 6700.0000000000000000 km Vel. Peri: 7.7131448355214580 km/sec Rad. Apo: 6700.0000000000000000 km Vel. Apo: 7.713144835521458 km/sec Mean Mot.: 0.06595979788783687 deg/sec Peri

32、od: 5457.869968191409 sec Period: 90.9644994698568 min Period: 1.51607499116428 hr Period: 0.06316979129851168 day Time Past Periapsis: 0 sec Time Past Ascending Node: -5457.869968191409 sec Beta Angle (Orbit plane to Sun): 22.0290337460786 deg Mean Sidereal Greenwich Hour Angle: 69.5612456281104 de

33、g Geodetic Parameters: Latitude: 0.01728125535160604 deg Longitude: -69.51750367179487 deg Altitude: 321.8630019297578400 km Geocentric Parameters: Latitude: 0.01717112559927386 deg Longitude: -69.51750367179487 deg Spacecraft Configuration: Drag Area: 1e-006 km2 SRP Area: 1e-006 km2 Dry Mass: 500 k

34、g Fuel Mass: 500 kg Total Mass: 1000 kg Area/Mass Ratio: 1e-009 km2/kg Tank Pressure: 5000 Pa Fuel Density: 999999999999.9999 kg/km3 Cr: 2.000000 Cd: 2.000000 *-* MCS Segment Type: Propagate Name: Propagate User Comment: -Propagation Description- Propagator model used: Earth_Point_Mass (Simple numer

35、ical twobody) Stopping Condition Information (Gregorian Date Julian Date): Propagation Statistics: Number of steps: 581 Average step size: 74.4868 sec Largest step size: 76.812 sec Smallest step size: 12.392 sec - Satellite State at End of Segment: - UTC Gregorian Date: 2 Jun 2003 00:00:00.00 UTC Ju

36、lian Date: 2452792.5 Julian Ephemeris Date: 2452792.50074287 Time past epoch: 43200 sec (Epoch in UTC Gregorian Date: 1 Jun 2003 12:00:00.00) State Vector in Coordinate System: Earth Centered Mean J2000 Parameter Set Type: Cartesian X: 5770.7341874767671000 km Vx: 3.9189766414393801 km/sec Y: -3404.

37、2072407123160000 km Vy: 6.6433594918105143 km/sec Z: 0.0000000000000000 km Vz: 0.0000000000000000 km/sec Parameter Set Type: Keplerian sma: 6700.0000000170658000 km RAAN: 0 deg ecc: 0.0000000000000336 w: 0 deg inc: 0 deg TA: 329.4632687491111 deg Parameter Set Type: Spherical Right Asc: 329.46326874

38、91111 deg Horiz. FPA: 9.565580324710149e-013 deg Decl: 0 deg Azimuth: 90 deg |R|: 6700.0000000172613000 km |V|: 7.7131448355114101 km/sec Other Elliptic Orbit Parameters : Ecc. Anom: 329.4632687491121 deg Mean Anom: 329.4632687491131 deg Long Peri: 0 deg Arg. Lat: 329.4632687491111 deg True Long: 32

39、9.4632687491111 deg Vert FPA: 89.99999999999905 deg Ang. Mom: 51678.07039805959 km2/sec p: 6700.0000000170658000 km C3: -59.49260325357981 km2/sec2 Energy: -29.7463016267899 km2/sec2 Vel. RA: 59.46326874911019 deg Vel. Decl: 0 deg Rad. Peri: 6700.0000000168411000 km Vel. Peri: 7.7131448355118941 km/sec Rad. Apo: 6700.0000000172904000