数学建模降落伞的选择Word文档下载推荐.docx
《数学建模降落伞的选择Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数学建模降落伞的选择Word文档下载推荐.docx（22页珍藏版）》请在冰点文库上搜索。
半径r/m
2
2.5
3
3.5
4
费用/元
65
170
350
660
1000
表格1不同半径的降落伞伞面价格
t/s
6
9
12
15
18
21
24
27
30
x/m
500
470
425
372
317
264
215
160
108
55
1
表格2降落伞试验的时刻t与高度x的观测值
问题的分析
这是一个有约束的优化问题，目标函数是降落伞的总费用，为了实用上的方便，不妨只选一种规格（伞半径）的降落伞，于是总费用是降落伞的个数与每个降落伞价格的乘积，而你决策变量是降落伞数量（记作n）和每个伞的半径r。
虽然r和n都只能取有限个离散值，但是，对n和r的各种组合进行枚举计算，逐个验算是否满足约束条件，比较费用，是相当繁琐的，并且缺乏一般性。
我们宁可先将n和r看作连续变量，建立优化模型，求得最优解后，再按题目要求做适当调整。
约束条件主要是伞的落地速度不能超过20m/s，为表述这一条件需要建立并求解降落伞速度满足的微分方程，而方程中的重要参数——空气阻力系数——又要通过测量数据（表格2）作拟合得到。
显然，由于测量数据是时间和高度，所以需要找出速度和高度之间的关系。
确定费用函数的关键是找出伞面价格与伞半径的关系，它可以根据所给数据（表格1）用适当的函数来拟合，观察这些数据的散点图，用幂函数拟合比较合适。
建立降落伞下落的微分方程时，关键是对所受阻力的分析，显然，阻力随着降落速度和伞面积的增加而变大。
模型假设
（1）伞面价格c1与伞半径r的关系，用幂函数c1=arb（a,b为待定参数）按照表格1数据来拟合；
载重m位于球心下方面处，每根绳索的长度l=
r。
（2）降落伞在空中只受到向下的重力和向上的空气阻力的作用，阻力与降落速度和伞面积的乘积成正比，阻力系数用表格2数据作拟合；
降落伞初速为零。
符号说明
降落伞数目
所投物的质量
物体在
时刻的高度
时刻的速度
阻力系数
降落伞的面积
重力加速度
时间
降落伞的半径
降落伞的总费用
每个降落伞的三面价格
每个降落伞的绳索价格
其他费用
模型建立
（1）目标函数
n个降落伞的总费用，记作C。
每个降落伞的费用由伞面价格c1=arb,绳索价格
和其他费用c3=200组成，于是
（20）
（2）伞的速度和高度
记时刻t伞的速度为v（t）,高度为x（t）,空气阻力为kr2v,k是待定参数，按照牛顿第二定律，v（t）满足
（21）
其中，m=2000/n（一个降落伞的载重），g=9.8m/s2,方程（21）的解为
（22）
对速度函数积分，并注意到t=0时x=500，得到伞的高度x（t）为
（23）
（3）约束条件
降落伞落地速度不超过20m/s,即当（23）式的x=0时的解得的根t,代入（22）式后满足
，此外还有
的附加条件，
整个优化模型可记作
（24）
当参数a,b,k用所给数据拟合确定后，即可求解模型（24）得到n，r（实数值），然后再作适当调整。
模型求解
（1）
参数估计
a，b的估计：
先将然后对于表14.5数据用线性最小二乘法和MATLAB软件编程[1]
得到a=4.3039，b=3.9779
与数据的拟合效果见图14.5.
图表1由表14.5数据拟合参数a,b
下面取a=4.3b=4
K的估计：
用表14.6数据估计k，注意到作降落试验时n=1，m=300，r=3，于是（23）式应改为
（25）
有表14.6数据利用MATLAB软件作非线性最小二乘法拟合，编程[2]：
得到k=18.4583
与数据的拟合效果见图14.6.下面取k=18.5
图表2由表14.6数据拟合参数k
（2）优化模型求解
将参数估计得到的a=4.3，b=4，k=18.5代入优化模型（24），用MATLAB的优化工具箱求解编程[3]
得到x=6.00722.969527.0408
c=4.8245e+003
根据题目要求，将结果调整为n=6，r=3，验证落地速度是否不超过20m/s,为此，先由（23）求解非线性方程：
（26）
再将得到的t代入（22）式计算出落地速度，编程[4]：
得到t=27.4867v=19.6196
落地速度符合要求。
最后，按照（20）式计算总费用（其中c1用实际价格350元），得到
C=6*（350+90.5*3+200）=4920（元）。
添加问题
风速偏向与风速大小将影响降落伞下降过程中的拉直时间，从而影响降落伞落地时的速度大小，请分析出风向与风速对拉直时间的关系，确定适合利用降落伞空投物资的风速范围以及风速的偏向。
问题分析
降落伞拉直过程是降落伞开伞过程的一个重要阶段，早起拉直过程动力学模型，都是假设伞系统和气流方向一致，处于理想的直线状态。
但实际上降落伞在拉直过程中，会收到各种因素的影响，导致降落伞的拉直方向几乎不可能与气流速度方向一致。
当在拉直过程中受到风的影响时，拉直时间会随之改变，从而影响落地的速度，为了保证物资可以安全的抵达地面，我们需要建立风速与风向对降落伞拉直时间的关系模型，从中解出适合利用降落伞空投物资的风速与风向范围。
由于限制条件要求降落伞落地时的速度不能超过20m/s，以及上面已经求出的最优解的限制，我们可以得出在高度H处必须将降落伞打开，于是对拉直时间T做出了要求，即
。
所以我们需要先利用已知的数据，做出模型，在
的范围内，求解出风速与风向的范围
（1）在降落伞拉直打开的过程中，我们假设此时降落伞做的是自由落体运动，即此时不受阻力的影响。
（2）假设在伞面打开以后，降落伞的运动过程与上一问中的运动情况完全相同。
拉直时间
最低开伞高度
风向偏角
风速大小
约束条件
上一问题中，在不加拉直时间的情况下，我们求出费用最低的最优解，此时k=18.5，n=6，r=3，代入下式中：
（27）
假设落地速度v=20m/s，我们可以解出t=27.0860s,H-=491.6858m
又因为拉直过程做的是自由落体运动，有
（28）
解得，
（29）
（1）变量说明
在该模型中，我们只研究风速偏向，以及风速大小对拉直时间的影响，其他外界环境因素可先忽略。
在此，我们对风速偏向做出三维图解，如下图所示：
（2）数据分析
偏角
风速
0°
20°
40°
60°
80°
100°
120°
140°
160°
180°
0m/s
1.27
30m/s
1.37
1.35
1.29
1.22
1.16
1.15
1.18
1.19
1.20
1.21
60m/s
1.50
1.40
1.26
1.10
0.91
0.82
1.07
1.14
经过网络数据统计，我们得到表1-1、表1-2
表1-1拉直时间随风速偏角的变化
风速（m/s）
5
10
20
25
拉直时间（s）/0度顺风
1.28
1.3
1.32
1.34
拉直时间（s）/90度顺风
1.25
1.24
拉直时间（s）/180度逆风
1.23
35
40
45
50
60
1.4
1.43
1.45
1.47
1.49
1.5
1.11
1.08
1.05
1.01
0.98
0.95
表1-2拉直时间随风速大小的变化
利用MATLAB，导入数据得到图1-1、图1-2
图1-1
图1-2
（3）数据处理
风速偏向与拉直时间关系
由图1-1我们可以看出，风速0m/s与风速60m/s是两个边界图像，即0m/s~60m/s之间的风速值所对应的图像都在0m/s图像和60m/s图像之间，因此我们只需考虑0m/s和60m/s两种情况即可。
（a）风速大小为0m/s时
由散点图可知拟合函数为三角函数，我们可以假设风速偏向与拉直时间的关系为：
（30）
将数据代入matlab软件中进行拟合（编程[5]），可以得到：
利用matlab软件，进行误差分析,得到：
方差=0.02482
标准差=0.05955
根据方差和标准差分析，可知误差相对较小，因此我们可以得到风速偏向于拉直时间之间的函数关系为：
当风速为0m/s时
根据式（29），即
，我们排除偶然误差等因素，可以看出当风速大小为0m/s时，偏向角范围是0~180°
，即偏向角可以取任意值。
（b）风速大小为60m/s时
（31）
将数据代入matlab软件中进行拟合（编程[5]），可以得到：
利用matlab软件，进行误差分析,得到：
方差=0.001087
标准差=0.01246
当风速为60m/s时
，我们可以得到风速为60m/s时，风速偏向范围是:
26.8717°
≤θ≦180°
风速大小与拉直时间关系
由图1-2我们可以看出，拉直时间t与风速大小v成正比关系，有
t=av+b（30）
图中，0°
顺风呈递增图像，90°
侧风和180°
逆风都成递减图像，因此我们只需对0°
顺风的情况进行研究，若0°
顺风时满足约束条件的要求，则其他偏向角也可满足条件。
图1-3
利用上面所提供的数据，利用matlab，将数据进行拟合（编程[6]），如图1-3我们可以得到：
a=0.0041
b=1.2593
方差=0.0005863
标准差=0.0073
拟合图像如图1-3所示
即可得到拉直时间t与风速大小v之间的函数关系为：
t=0.0041v+1.2593（31）
，我们可以得到风速的大小范围是:
（4）总结说明
为了使物资安全抵达地面，根据最后计算结果，我们可以得到适合空投物资的环境为：
风速偏角范围：
风速大小范围：
附录：
编程[1]：
r=[22.533.54];
c1=[651703506601000];
lgc1=log（c1）;
lgr=log（r）;
A=polyfit（lgr,lgc1,1）;
b=A
（1）;
a=exp（A
（2））;
rr=2:
0.01:
4;
cc1=a*rr.^b;
plot（r,c1,'
+'
rr,cc1）,grid
编程[2]：
functionf=sanf（k,t）
f=500-980*t/3/k+98000*（1-exp（-3*k*t/100））/9/k^2;
t=[036912151821242730];
x=[500470425372317264215160108551];
k0=10;
k=lsqcurvefit（@sanf,k0,t,x）;
tt=0:
0.1:
30;
f=500-980*tt/3/k+98000*（1-exp（-3*k*tt/100））/9/k^2;
plot（t,x,'
tt,f）,grid
编程[3]：
functionc=sanc（x）%目标函数,
x
（1）=n,x
（2）=r
c=x
（1）*（4.3*x
（2）^4+90.5*x
（2）+200）;
function[cc,cceq]=sancc（x）%约束条件，x（3）=t
k=18.5;
p=19600/k/x
（1）/x
（2）^2;
q=1-exp（-k*x
（1）*x
（2）^2*x（3）/2000）;
cc=p*q-20;
%不等式约束
cceq=500-p*x（3）+p^2*q/9.8;
%等式约束
x0=[3,4,30];
v=[1,2,0];
u=[10,4,50];
[x,c]=fmincon（@sanc,x0,[],[],[],[],v,u,@sancc）
编程[4]
functionff=sany（t,k,n,r）
p=19600/k/n/r^2;
q=1-exp（-k*n*r^2*t/2000）;
ff=500-p*t+p^2*q/9.8;
n=6;
r=3;
t0=30;
t=fzero（@sany,t0,[],k,n,r）
v=p*q
编程[5]：
v=[0pi/92*pi/9pi/34*pi/9pi/25*pi/92*pi/37*pi/98*pi/9pi];
t=pianjiao（4,:
）;
t1=pianjiao（2,:
再利用cftool工具里的三角函数拟合进行拟合
编程[6]
v=fengsu（1,:
t=fengsu（2,:
plot（v,t）
polyfit（v,t,1）
holdon
t1=0.0041*v+1.2593;
plot（v,t1,'
r*-'
）
%xlable（'
风速偏角/（°
）'
%x轴注解
%ylable（'
拉直时间/s'
%y轴注解
title（'
拉直时间随风速的变化'
%图形标题
legend（'
原始数据'
'
拟合函数'
%图形注解
gridon;
%显示网格线
升级会员 