系统仿真大作业飞机剖析.docx
《系统仿真大作业飞机剖析.docx》由会员分享，可在线阅读，更多相关《系统仿真大作业飞机剖析.docx（20页珍藏版）》请在冰豆网上搜索。
系统仿真大作业飞机剖析
作业要求
系统仿真作业：
观测站（O点）为测得的某航班数据，当飞机到达某位置P时开始对飞机的相关数据进行记录，在位置P处时间t记为0，当飞机到观测站的距离达到最短时飞机所处的位置记为M点，飞机在t时刻所处的位置与观测站O点的连线到直线OM的夹角记为θ（t）。
附表一给出了观测站某次记录数据：
A列为记录时间t，其间隔为0.001s；B列为飞机飞行过程中t时刻飞机相对于P点的距离S（t）；C列给出了飞机t时刻角度θ（t）的理论值theta_theory（参考输入），D列给出了观测站实际上观测到的飞机在t时刻时角度θ（t）的观测值theta_observation；E列给出了当把C列数据theta_theory作为某标准二阶伺服系统G（s）的输入信号时该标准二阶伺服系统的输出值theta_output。
作业要求：
（1）将附表一中的数据导入matlab工作空间，使各列数据都能作为变量使用。
（2）试根据表格中的数据使用MATLAB完成以下问题
1.根据A、B两列数据确定飞机飞行时的理论运行轨迹和飞机的飞行速度；
2.根据A、B、C三列数据确定当飞机到达M点时观测站O到M点的距离。
（3）设计标准二阶伺服系统G（s）
要求：
1．确定合适的ζ值，其单位阶跃响应的性能指标满足:
系统的超调量σ%介于4.5%~8.0%之间
2.使用MATLAB仿真确定参数ωn,使得：
当把C列数据theta_theory作为该二阶伺服系统的输入时，系统的输出尽可能的接近E列所给出的theta_output。
求出系统传递函数参数。
3.将C列theta_theory作为设计好的标准二阶伺服系统的输入信号，计算输出相对于输入的相对误差；作出输出随时间变化的曲线，以及相对误差随时间变化的曲线。
（4）在上述设计好的伺服系统中加入合适的干扰信号和适当的非线性环节，使得当把C列数据theta_theory作为输入信号时系统的输出尽可能接近D列数据（θ（t）的观测值theta_observatuon）。
试：
1.确定干扰信号和非线性环节的相关参数，并画出系统simulink模型。
2.画出系统的输出随时间变化的曲线。
3.画出输出相对于输入theta_theory的相对误差随时间变化的曲线，以及输出相对于D列数据（θ（t）的观测值theta_observatuon）的相对误差曲线。
一、完成内容
（1）将附表一中的数据导入matlab工作空间，使各列数据都能作为变量使用。
编写脚本readxlsx.m，读取各列数据的值：
%将附表一中的数据导入matlab工作空间，使各列数据都能作为变量使用
[typ,desc,fmt]=xlsfinfo（'finaltest.xlsx'）
%读取表格文件的信息
t=xlsread（'finaltest.xlsx',1,'A2:
A143002'）;
%读取excel文件中第一个工作表中A列的全部数据，用数组t（）保存数据
s=xlsread（'finaltest.xlsx',1,'B2:
B143002'）;
%读取excel文件中第一个工作表中B列的全部数据，用数组s（）保存数据
theta_theory=xlsread（'finaltest.xlsx',1,'C2:
C143002'）;
%读取excel文件中第一个工作表中C列的全部数据，用数组theta_theory（）保存数据
theta_observation=xlsread（'finaltest.xlsx',1,'D2:
D143002'）;
%读取excel文件中第一个工作表中D列的全部数据，用数组theta_observation（）保存数据
theta_output=xlsread（'finaltest.xlsx',1,'E2:
E143002'）;
%读取excel文件中第一个工作表中E列的全部数据，用数组theta_output（）保存数据
formatlong;
%设定显示格式为long，因为表格中数据有9位有效值
运行结果：
依次可取得各数据的值
>>t
（2）
ans=
1.000000000000000e-03
>>s
（2）
ans=
0.140000000000000
>>theta_theory
（2）
ans=
-1.520837232808838
>>theta_observation
（2）
ans=
-7.075789556117257e-05
>>theta_output
（2）
ans=
-7.573799703386787e-05
（2）试根据表格中的数据使用MATLAB完成以下问题
●根据A、B两列数据确定飞机飞行时的理论运行轨迹和飞机的飞行速度；
●根据A、B、C三列数据确定当飞机到达M点时观测站O到M点的距离。
编写脚本文件plane_rate.m
plot（t,s）;gridon;
%作飞机位移关于时间t的图象，可看出飞机是匀速运动
i=2:
143001;
rate（i）=s（i）./t（i）;
%求运行轨迹中各点的速率，可看出运动速率恒为140m/s
i=1:
143001;
abs_theta_theory（i）=abs（theta_theory（i））;
%取观察站测得角度的绝对值
[min_theta_theoryk]=min（abs_theta_theory）;
%当飞机到达M点时，观察站测得角度应为最小角，故求得此最小角度时
%对应的时间t（k）
distance_OM=s（k）*cot（abs_theta_theory
（1））
%当时间为t（k）时，对应的位移为s（k），tan（theta）=PM/OM
%其中PM=s（k）,theta=theta_theory的最初的值，由此可求得OM的距离
运行结果：
>>planerate
distance_OM=
5.000030000000002e+02
故飞机飞行速率为140m/s,OM间距离为500m。
（3）设计标准二阶伺服系统G（s）
解：
1．确定合适的ζ值，其单位阶跃响应的性能指标满足:
系统的超调量σ%介于4.5%~8.0%之间
因为σ%=
，使σ%介于4.5%~8.0%之间则可取ζ=0.68，此时σ%=%5.4280251339488，满足要求。
2.使用MATLAB仿真确定参数ωn,使得：
当把C列数据theta_theory作为该二阶伺服系统的输入时，系统的输出尽可能的接近E列所给出的theta_output。
求出系统传递函数参数。
用simulink设计标准二阶伺服系统，其中w的值待定设计的框图如下：
再编写脚本文件，如下：
forn=1:
20;%取不同的w值时，求输出的theta_output
sim（'erjiesifuxitong'）;%运行simulink中的框图文件，
%获得输出的的p_theta_output值，为一个长度为143001的列向量
forj=1:
143001;
d_value（j）=abs（p_theta_output（j）-theta_output（j））;
%求每组输出的p_theta_output与理论值的差
end
abs_dif（n）=sum（d_value）;%取误差和，存到向量ads_dif中，再依次比较向量中
%的值，当取最小值时说明p_theta_output与理论值最接近，此时
%的w值即为标准w值
end
当w取1～20时，得到的abs_dif各值为：
>>abs_dif
abs_dif=
1.0e+03*
Columns1through4
6.0194005109116212.7395408824277221.6440888630738101.100263457041429
Columns5through8
0.7732417253371810.5490055141110910.3840511989651720.256977609508468
Columns9through12
.0791********
Columns13through16
0.1099810646631350.1495906571967200.183********77030.213275684364825
Columns17through20
0.2389980608554720.2615401303097590.2814228978999480.299065677805618
由数据可看出，当w=11时，取得最小的差值，再取w=10~12之间的数，求得其中最小差值。
得到的simulink框图为：
3.将C列theta_theory作为设计好的标准二阶伺服系统的输入信号，计算输出相对于输入的相对误差；作出输出随时间变化的曲线，以及相对误差随时间变化的曲线。
作输出变化曲线和相对误差曲线
求相对误差的脚本文件为：
sim（'erjiesifuxitong'）;%运行二阶伺服系统
fori=1:
143001;
abs_dif（i）=abs（p_theta_output（i）-theta_theory（i））;
%求输出p_theta_output的绝对误差
relative_error（i）=abs_dif（i）/abs（theta_theory（i））;
%求输出的相对误差
End
作输出变化曲线为：
当作相对误差曲线时，发现为如下图所示：
接近于冲激函数。
因为时间轴选用的过宽，无法看出细节，故用axis（[0,2,0,2]）指令截取相应的曲线，就能看清楚了：
a）时间轴为0～2时
plot（t,relative_error）;axis（[0,2,0,1]）
b）时间轴为2～70时，
图形几乎无变化，即相对误差近似为0.
c）时间轴为70～73时，plot（t,relative_error）;axis（[70,73,0,40]）
由图可看出相对误差变化图在70～73之间有一个猛增的过程
d）时间轴为73～143时，plot（t,relative_error）;axis（[73,143,0,1]）
分析图可看出，在70s～73s之间，相对误差有一个显著的冲激式上升过程，分析原因：
相对误差的求取公式为
而theta_theory的值在70s～73s时近似取到了零，而theta_theory在分母上，故此时相对误差值会变得极大。
当theta_theory取值远离零值时，由于绝对误差接近于0，故相对误差也接近于0。
而在0s~2S区间可看到一个明显的振荡过程，分析原因：
因为系统为二阶系统，输入的起始阶段响应会有超调量，但当经过一段时间系统稳定后，就不会再发生振荡。
4.在上述设计好的伺服系统中加入合适的干扰信号和适当的非线性环节，使得当把C列数据theta_theory作为输入信号时系统的输出尽可能接近D列数据（θ（t）的观测值theta_observation）。
试：
1）确定干扰信号和非线性环节的相关参数，并画出系统simulink模型。
●作出theta_theory的时间曲线和theta_observation的时间曲线如图，蓝色为theta_theory的时间曲线，红色为theta_observation的时间曲线。
编写M文件，使每选定一种非线性环节或添加一种干扰都能得到对应的曲线，再将对应曲线与给定的观察曲线作比较，即能找到对应的参数。
Thetaobservation.m
%plot（t,theta_theory）;holdon;%作理论曲线
plot（t,theta_observation,'r'）;holdon;%作观察所得曲线
sim（'erjiesifuxitong'）;%运行仿真框图
plot（t,p_theta_output）;%作仿真结果曲线
holdon;
首先在前向通道中加一个死区非线性环节，得到如下图所示波形：
此时输出曲线与观察曲线相差较大，故不正确。
再将死区非线性放到回路外：
得到曲线：
输出曲线在中间停滞区与观察曲线近似，经过多次尝试，发现将死区的起始与截止参数改为如下所示：
此时能得到与观察曲线非线性弯曲部分重合的曲线：
此即找到非线性环节。
●为找到干扰环节，我作出了极大的努力。
首先将观察曲线与输出曲线作对比，如下图：
为求出合适的干扰信号使蓝色曲线接近红色曲线，先将两曲线作差（i=1:
143001;diff（i）=p_theta_output（i）-theta_observation（i）;plot（t,diff）;）得到如下曲线：
由图可看出，差值diff的时间曲线近似为一个正弦波形，故应输入一个正弦形状的干扰信号。
为求出干扰信号的具体差值，我再对diff进行一次线性拟合，编写函数文件如下：
functionfun_sin=xianxingnihe（x,t）
fori=1:
143001
fun_sin（i）=x
（1）+x
（2）*sin（4.5*pi/143*t（i）+x（3）*pi）;
end
再输入指令行
>>[x,resnorm]=lsqcurvefit（@xianxingnihe,x0,t,diff）,即对diff曲线进行输入为时间t的线性拟合，拟合结果如下：
x=
0.000135938260289
-0.049733591012013
.020*********
resnorm=
2.238947129380221
再作出拟合所得曲线与真实的diff曲线作比较：
plot（t,diff）;holdon;plot（t,xianxingnihe（x,t）,'r'）;
得到图为：
由图可知，拟合结果近似接近于原曲线，因为原曲线幅值较小，故将此拟合结果作为干扰信号输入。
输入系统框图如图：
其中正弦波参数参照拟合曲线采用，如下
输出此框图的结果为：
由图可看出，两条曲线近似重合，故此系统即为标准的输入系统。
2）画出系统的输出随时间变化的曲线。
sim（'erjiesifuxitong'）;%运行仿真框图
plot（t,p_theta_output）;%作仿真结果曲线
3）画出输出相对于输入theta_theory的相对误差随时间变化的曲线，以及输出相对于D列数据（θ（t）的观测值theta_observatuon）的相对误差曲线。
●输出相对于输入theta_theory的相对误差随时间变化的曲线：
i=1:
143001;
diff_a（i）=abs（p_theta_output（i）-theta_theory（i））./abs（theta_theory（i））;
plot（t,diff_a）;holdon;
axis（[014301]）;
●输出相对于D列数据（θ（t）的观测值theta_observatuon）的相对误差曲线
i=1:
143001;
diff_a（i）=abs（p_theta_output（i）-theta_observation（i））./abs（theta_observation（i））;
plot（t,diff_a）;holdon;
axis（[014300.1]）;
图形分析：
i.在曲线的起始阶段相对误差有极大的震荡过程，因为在此阶段系统还没有稳定，系统输出值具有较大变动。
ii.在曲线的70s～80s范围内有一个相对误差猛增的过程，分析原因应为此阶段theta_observation的值无限接近于0，故相除得到误差无限大。
iii.在曲线其他阶段，相对误差保持在0.01以下，故认为设计的系统满足题设要求。
升级会员 