飞行器系统仿真.docx
《飞行器系统仿真.docx》由会员分享，可在线阅读，更多相关《飞行器系统仿真.docx（28页珍藏版）》请在冰豆网上搜索。
飞行器系统仿真
《飞行器系统仿真与CAD》学习报告
第一部分仿真（40）
题目1：
给定导弹相对于目标的运动学方程组为
r（0）=5km,q（0）=60deg,（0）=30deg,V=,Vm=,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）
switchflag
case0
[sys,x0,str,ts]=mdlInitializeSizes;
case1
sys=mdlDerivatives（t,x,u）;
case3
sys=mdlOutputs（t,x,u）;
case{2,4,9}
sys=[];
otherwise
error（'unhandledflag=',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=[00];
functionsys=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;
functionsys=mdlOutputs（t,x,u）
sys=x;
调用S函数的模型框图
（3）框图仿真结果：
S函数仿真结果：
（4）命令输入
clear;clc
[tx]=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'）;
[tx]=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：
给出动态方程
;
（1）用MATLAB语言编写S—函数;
（2）用命令行gear/adams法对
（1）进行仿真，显示曲线x（t=0：
100）;
（3）建立方框图,用RK45仿真50秒,显示曲线
答：
（1）用MATLAB语言编写S—函数
function[sys,x0,str,ts]=CAD2_sfun（t,x,u,flag）
switchflag
case0
[sys,x0,str,ts]=mdlInitializeSizes;
case1
sys=mdlDerivatives（t,x,u）;
case3
sys=mdlOutputs（t,x,u）;
case{2,4,9}
sys=[];
otherwise
error（'unhandledflag=',num2str（flag））
end
function[sys,x0,str,ts]=mdlInitializeSizes
sizes=simsizes;
=2;
=0;
=2;
=0;
=1;
=1;
sys=simsizes（sizes）;
str=[];
x0=[1,0];
ts=[00];
functionsys=mdlDerivatives（t,x,u）
dx
（1）=x
（2）;
dx
（2）=1-t*x
（1）-（1-x
（1）^2）*x
（2）;
sys=dx;
functionsys=mdlOutputs（t,x,u）
sys=x;
（2）直接调用ode数值积分函数进行仿真，系统微分方程：
functiondx=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;
[tx]=ode15s（@CAD01_02odefun,0:
100,[10]）;
hGear=figure（）;
set（hGear,'NumberTitle','off','Name','IntegratedbytheGearalgorithm','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'）;
[tx]=ode113（@CAD01_02odefun,0:
100,[10]）;
hAdams=figure（）;
set（hAdams,'NumberTitle','off','Name','IntegratedbytheAdamsalgorithm','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=1kg,K=4kg/s2,C=100kg/m,g=s2,=;;
（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）;
[rowcol]=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*
0
B=
0
1
C=
10
01
D=
0
0
Thesystemiscontrolled
（3）阶跃响应：
脉冲响应:
（4）对原模型进行仿真，P=sin（t）（使用Simulink）
仿真结果：
（5）对原模型进行仿真，P=sin（t）（使用ode23）
系统微分方程：
functiondx=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-55e-5]）;
[tx]=ode23（@CAD3odefun,0:
:
10,[00],options）;
hode23=figure（）;
set（hode23,'NumberTitle','off','Name','Integratedbytheode23solver',...
'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,ul,ur,yl,yr）
switchflag,
case0,
[sys,x0,str,ts]=mdlInitializeSizes;
case3,
sys=mdlOutputs（t,x,u,ul,ur,yl,yr）;
case9,
sys=[];
end
function[sys,x0,str,ts]=mdlInitializeSizes
sizes=simsizes;
=0;
=0;
=1;
=1;
=1;
=1;
sys=simsizes（sizes）;
x0=[];
str=[];
ts=[00];
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）&&（uy=u-ul;
elseif（uur）
y=u-ur;
else
y=0;
end
sys=y;
在Simulink中将调用S函数的模块进行封装
参数传递及初始化
用户界面：
测试结果
问问题5：
已知系统A=[01;-1-2],B=[10;01],C=[10;01],D=[00;00],求系统的状态空间方程（linmod），并分析系统的稳定性，练习仿真参数设置
答：
对模型进行线性化并分析稳定性
clear;clc;
[ABCD]=linmod（'CAD5'）
ss0=ss（A,B,C,D）;
hpz=figure（）;
set（hpz,'NumberTitle','off','Name','Pole-zeromapofthelinmodsystem'）;
pzmap（ss0）;
sgrid;
[rowcol]=size（A）;
P=lyap（A,eye（row））;
fori=1:
row
subdet（i）=det（P（1:
i,1:
i））;
end
subdet
系统零极点图：
存在正实部的极点，系统不稳定。
问题6：
系统的动力学方程为dx/dt=Ax+Bu,y=Cx+Du,A=[0100;0010;0001;-1-2-3-4],B=[12;34;23;45],C=[1122;2354];D=[10;01],求：
（1）系统动态平衡点
（2）x（0）=[1111]’,ix=[1234]’,dx=[0101]’,idx=[1234]’,的系统动态平衡点
答：
系统框图模型
系统的平衡点分析
程序
clear;clc;
[x,u,y,dx,options]=trim（'CAD6'）;
options（10）
x0=[1111];ix=[1234];dx=[0101];idx=[1234];
[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函数源代码
#defineS_FUNCTION_NAMECAD02_07sfun
/***Modified:
changethefunctionname*/
#defineS_FUNCTION_LEVEL2
#include""
staticvoidmdlInitializeSizes（SimStruct*S）
{
ssSetNumSFcnParams（S,1）;/***Revised:
setthenumberofinputparameters
to1*/
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）;
}
staticvoidmdlInitializeSampleTimes（SimStruct*S）
{
ssSetSampleTime（S,0,INHERITED_SAMPLE_TIME）;
ssSetOffsetTime（S,0,;
}
staticvoidmdlOutputs（SimStruct*S,int_Ttid）
{
int_Ti;
InputRealPtrsTypeuPtrs=ssGetInputPortRealSignalPtrs（S,0）;
real_T*y=ssGetOutputPortRealSignal（S,0）;
int_Twidth=ssGetOutputPortWidth（S,0）;
constmxArray*pmxRatio=ssGetSFcnParam（S,0）;/***Revised:
getthe
pointertotheparameterinthetypeofmxArray*/
constreal_T*pRatio=mxGetPr（pmxRatio）;/***Revised:
getthepointerto
theparameterinthetypeofreal_T*/
for（i=0;i*y++=（*pRatio）*（*uPtrs[i]）;/***Revised:
themagnifyingratiois
acquiredfromtheinputparameter*/
}
}
staticvoidmdlTerminate（SimStruct*S）{}
#ifdefMATLAB_MEX_FILE
#include""
#else
#include""
#endif
封装及用户界面：
问题8自学文件Stateflow示例文件
function[varargout]=stateflow（varargin）
%STATEFLOWOpensSIMULINKandcallssfnewwhenappropriate.
%Copyright1995-2002TheMathWorks,Inc.
ifnargout>0
[varargout{:
}]=sf（varargin{:
}）;
else
sf（varargin{:
}）;
end
Stateflow是有限状态机（finitestatemachine）的图形工具，它可以用于解决复杂的逻辑问题，用户可以通过图形化工具实现在不同状态之间的转换。
Stateflow可以直接嵌入到Simulink仿真模型中，并且在仿真的初始化阶段，SIMULINK会把Stateflow绘制的逻辑图形通过编译程序转换成C语言，使二者有机地结合在一起。
Stateflow可以在SIMULINKExtra模块库中找到。
Stateflow的仿真原理是有限状态机（finitestatemachine）理论，有限状态机是指系统含有可数的状态，在相应的状态事件发生时，系统会从当前状态转移到与之对应的状态。
在有限状态机中实现状态的转移是有一定条件的，同时相互转换的状态都会有状态转移事件，这样就构成了状态转移图。
在SIMULINK的仿真窗口中，允许用户建立有限个状态以及状态转移的条件与事件，从而绘制出有限状态机系统，这样就可以实现对系统的仿真。
Stateflow的仿真框图一般都会嵌入到Simulink仿真模型中，同时实现状态转移的条件或是事件即可以取自Stateflow仿真框图，也可以来自Simulink仿真模型。
第二部分综合（60）
2.
J2
HPOP
3采用STK轨道机动模块实现任意两个圆轨道间的霍曼转移
答：
利用轨道转移模块Astrogator，结合轨道机动打靶功能Target实现霍曼转移。
对卫星的Astrogator设置如下的任务控制序列（MSC）。
09May201315:
17:
05
Satellite-Satellite1
AstrogatorMissionControlSequenceSummary
---------------------------------------------------------------------------------------------------
***--------------------------------------------------------------------------***
MCSSegmentType:
InitialState
Name:
InitialState
UserComment:
-InitialStateDescription-
----------------------------------
SatelliteStateatEndofSegment:
----------------------------------
UTCGregorianDate:
1Jun200312:
00:
UTCJulianDate:
2452792
JulianEphemerisDate:
Timepastepoch:
0sec（EpochinUTCGregorianDate:
1Jun200312:
00:
StateVectorinCoordinateSystem:
EarthCenteredMeanJ2000
ParameterSetType:
Cartesian
X:
kmVx:
km/sec
Z:
kmVz:
km/sec
ParameterSetType:
Keplerian
sma:
kmRAAN:
0deg
ecc:
w:
0deg
inc:
0degTA:
0deg
ParameterSetType:
Spherical
Decl:
0degAzimuth:
90deg
OtherEllipticOrbitParameters:
Ecc.Anom:
0degMeanAnom:
0deg
LongPeri:
0degArg.Lat:
0deg
TrueLong:
0degVertFPA:
90deg
Vel.RA:
90degVel.Decl:
0deg
TimePastPeriapsis:
0sec
GeodeticParameters:
GeocentricParameters:
SpacecraftConfiguration:
DragArea:
1e-006km^2
SRPArea:
1e-006km^2
DryMass:
500kg
FuelMass:
500kg
升级会员 