飞行器系统仿真
《飞行器系统仿真与CAD》学习报告
第一部分仿真(40)
题目1：给定导弹相对于目标的运动学方程组为
r(0) = 5km, q(0) = 60deg, (0) = 30deg,V = , V
= , 1Ma = 340m/s, k = 2
m
（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 1
sys = mdlDerivatives(t,x,u);
case 3
sys = mdlOutputs(t,x,u);
case {2,4,9}
sys = [];
otherwise
error('unhandled flag=',num2str(flag))
end
function [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(2===+-+?
???x 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 flag
case 0
[sys,x0,str,ts]=mdlInitializeSiz es;
case 1
sys=mdlDerivatives(t,x,u);
case 3
sys=mdlOutputs(t,x,u);
case {2,4,9}sys=[];
otherwise
error('unhandled
flag=',num2str(flag))
end
function
[sys,x0,str,ts]=mdlInitializeSiz es
sizes=simsizes;
=2;
=0;
=2;
=0;
=1;
=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数值积分函数进行仿真，系统微分方程：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 ????
+++=()()2μM=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', 'Step
Response','unit','normalized','Position' ,[,,,]);
step(ss0);grid;
hImpulse = figure();
set(hImpulse, 'NumberTitle', 'off', 'Name', 'Impulse Response','unit','normalized','Position' ,[,,,]); impulse(ss0);grid;
命令窗口输出结果：
A =
+008 *
B =
1
C =
1 0
0 1
D =
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-函数
答：实现所需功能的S函数
function [sys,x0,str,ts] =
CAD01_04sfun_kernel(t,x,u,flag,u l,ur,yl,yr) switch flag, case 0,
[sys,x0,str,ts]=mdlInitializeSiz es;
case 3,
sys=mdlOutputs(t,x,u,ul,ur,yl,yr );
case 9,
sys=[];
end
function
[sys,x0,str,ts]=mdlInitializeSiz es
sizes = simsizes;
= 0;
= 0;
= 1;
= 1; = 1;
= 1;
sys = simsizes(sizes);
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)
y = yl;
elseif (u >ul + yl) && (u elseif (u y = u - ur; else y = 0; end sys = 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:row subdet(i) = det(P(1:i,1:i)); end subdet 系统零极点图： 存在正实部的极点，系统不稳定。 问题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 = [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];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=9 x=[;;;];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);/***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, DYNAMICALLY_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 = ssGetInputPortRealSignalPtrs(S,0); real_T *y = ssGetOutputPortRealSignal(S,0); int_T width = ssGetOutputPortWidth(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; i *y++ = (*pRatio) *(*uPtrs[i]);/***Revised: the magnifying ratio is acquired from the input parameter*/