飞行器系统仿真.docx
《飞行器系统仿真.docx》由会员分享,可在线阅读,更多相关《飞行器系统仿真.docx(31页珍藏版)》请在冰豆网上搜索。
飞行器系统仿真
《飞行器系统仿真与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)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]=mdlInitializeSizessizes=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;ylabe'lr(');
subplot(3,1,2);
plot(t,x(:
2));grid;ylabe'lq(');
subplot(3,1,3);
plot(t,x(:
3));grid;ylabe'ls(igma)';
[tx]=sim('CAD1_S');
hSFun=figure();
subplot(3,1,1);
plot(t,x(:
1));grid;ylabe'lr(');
subplot(3,1,2);
plot(t,x(:
2));grid;ylabe'lq(');
subplot(3,1,3);
plot(t,x(:
3));grid;ylabe'ls(igma)';
模型仿真结果:
S函数仿真结果:
题目2:
给出动态方程x(1x2)xtx1;x(0)1,x(0)0
⑴用MATLAB语言编写S-函数;
(2)用命令行gear/adams法对
(1)进行仿真,显示曲线x(t=O:
100);
⑶建立方框图,用RK45仿真50秒显示曲线
答:
(1)用MATLAB语言编写S—函数
function[sys,x0,str,ts]=CAD2_sfun(t,x,u,flotherwise
ag)
error('unhandledflag=,'num2str(flag))
switchflag
end
case0
function[sys,x0,str,ts]=mdlInitializeSizes
sizes=simsizes;
[sys,x0,str,ts]=mdlInitializeSizes;
=2;
case1
=0;
sys=mdlDerivatives(t,x,u);
=2;
case3
=0;
sys=mdlOutputs(t,x,u);
=1;
case{2,4,9}sys=[];
dx
(1)=x
(2);
sys=simsizes(sizes);
ts=[00];
functionsys=mdlOutputs(t,x,u)
sys=x;
functionsys=mdlDerivatives(t,x,u)
(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;ylabe'lx(');
subplot(2,1,2);
plot(t,x(:
2));grid;ylabe'ld(x/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;ylabe'lx(');
subplot(2,1,2);
plot(t,x(:
2));grid;ylabe'ld(x/dt');
ode15s(Gear)仿真结果:
ode113(Adams)仿真结果:
(3)建立方框图,用RK45仿真50秒,显示曲线方框图模型:
仿真结果:
问题3:
质量一弹簧系统,质量M,恢复系数K,阻力系数C,主动力P,动
力学方程为Mx(Cx2Mg)sign(x)KxPM=1kg,K=4kg/s2,C=100kg/m,
2o
g=s,?
=;;
(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','StepRespons,e'u'nit','normalized,''Position',[,,,]);step(ss0);grid;
hImpulse=figure();
set(hImpulse,'NumberTitle','off','Name','Impulse
Response,'unit','normalized','Position',[,,,]);
impulse(ss0);grid;
命令窗口输出结果:
A=
+008
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','Integratedbytheode23solve,.r.'.
'Units','Normalized','Position',[]);
subplot(2,1,1);
plot(t,x(:
1));grid;ylabe'lx(');
subplot(2,1,2);
plot(t,x(:
2));grid;ylabe'ld(x/dt');
仿真结果:
问题4:
给出一个系统,要求生成一个新Simulink模块,实现其功能
(1)Mask功能
(2)s-函数
答:
实现所需功能的S函数
function[sys,x0,str,ts]=[sys,x0,str,ts]=mdlInitializeSizes;
CAD01_04sfun_kernel(t,x,u,flag,ul,ur,yl,y
case3,
r)
sys=mdlOutputs(t,x,u,ul,ur,yl,yr);
switchflag,
case9,
case0,
sys=[];
function[sys,x0,str,ts]=mdlInitializeSizesif(u>=ur+yr)
sizes=simsizes;
=0;
=0;
=1;
=1;
=1;
=1;
sys=simsizes(sizes);
x0=[];
str=[];
ts=[00];
在Simulink中将调用S
y=yr;
elseif(u<=ul+yl)
y=yl;
elseif(u>ul+yl)&&(u
elseif(uur)y=u-ur;
else
y=0;
end
sys=y;
参数传递及初始化
用户界面:
测试结果
问题5:
已知系统A=[01;-1-2],B=[10;01],C=[10;01],D=[00;00],求系
统的状态空间方程(linmod),并分析系统的稳定性,练习仿真参数设置
答:
对模型进行线性化并分析稳定性
clear;clc;
[ABCD]=linmod('CAD5')
ssO=ss(A,B,C,D);
hpz=figure();
set(hpz,'NumberTitle','off,'Name',卩ole-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=[O101]的系统动态234]
平衡点
答:
系统框图模型系统的平衡点分析程序
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:
setthenumberofinputparametersto1*/
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:
themagnifyingratioisacquiredfromtheinputparameter*/
}
}
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设置如下的任务控制序列(MS6
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
ParameterSetType:
Keplerian
ecc:
w:
0deg
inc:
0deg
TA:
0deg
ParameterSetType:
Spherical
Decl:
0deg
Azimuth:
90deg
OtherEllipticOrbitParameters:
Ecc.Anom:
0deg
MeanAnom:
0deg
LongPeri:
0deg
Arg.Lat:
0deg
TrueLong:
0deg
VertFPA:
90deg
Vel.RA:
90deg
Vel.Decl:
0deg
TimePastPeriapsis:
0sec
GeodeticParameters:
GeocentricParameters:
SpacecraftConfiguration:
DragArea:
1e-006kmA2
SRPArea:
1e-006kmA2
DryMass:
500kg
FuelMass:
500kg
TotalMass:
1000kg
Area/MassRatio:
1e-009kmA2/kg
Pa
TankPressure:
5000
Cr:
Cd:
MCSSegmentType:
Propagate
Name:
Propagate
UserComment:
-PropagationDescription-
Propagatormodelused:
Earth_Point_Mass(Simplenumericaltwobody)
StoppingConditionInformation(GregorianDate[JulianDate]):
PropagationStatistics:
Numberofsteps:
581
Averagestepsize:
sec
Largeststepsize:
sec
Smalleststepsize:
sec
SatelliteStateatEndofSegment:
UTCGregorianDate:
2Jun200300:
00:
UTCJulianDate: