降落伞优化选择的整数线性规划模型

降落伞优化选择的整数线性规划模型摘要本文讨论了降落伞合理选择使费用最低的问题。

通过对问题的分析，最大化载重量，最小化选购降落伞费用。

以牛顿定律建立微分模型，以空投物资重量2000千克，每种降落伞最大载重量为约束条件建立整数线性规划模型。

通过分步优化，最后以整数规划来解决这一问题。

首先，找出数据之间的关系，运用物理学和整数线性规划建立模型，并运用MATLABR软件描点作图进行数据拟合的方法,得出载重为300kg，半径为3米的降落伞从500米高空下降时的运动曲线，发现降落伞后期趋于做匀速直线运动.当降落伞作匀速直线运动时，求出空气阻力系数为2.959,落地速度为17.5794.在求出每种降落伞最大载重量，并通过隔离载重物体并进行受力分析，求出相应半径降落伞绳索长度,进而算出每种半径的降落伞的绳索费。

最后，根据每种降落伞的总成本关系把问题转化为整数线性规划问题，用LINGO解得到要购买半径为3m的降落伞数量为6把时总费用最少，总费用为4932元。

本文主要研究了降落伞优化选择问题。

主要优点是：本文通过建立优化选择的整数线性规划模型求解，思路清晰，并大量运用计算机运算使计算误差减少，最终使得降落伞的选择最优;另一方面,本文所建的模型简单合理,具有较强的推广意义。

主要缺点：在建立模型时，忽略了降落伞在实际应用中，会受到天气、风等一些自然因素的影响，使得模型与实际有些误差；本模型未考虑降落伞打开时间，将其假设成在下降时伞就已经打开；虽然大量运用计算机运算，但其中还是有不可避免的误差。

关键词: 数据拟合；单目标优化；微分方程；整数线性规划.一、问题的提出：为向灾区空投救灾物资共2000kg，需选购一些降落伞。

已知空投高度为500m，要求降落伞落地时的速度不能超过20m/s。

降落伞面为半径r的半球面，用每根长l共16根绳索连接着载重m，示意图如图1。

图1每个降落伞的价格由3部分组成。

伞面价格由半径r决定（见表1）；绳索每米为4元，其他费用200元。

降落伞在降落过程中受到空气的阻力，可以认为与降落伞速度和伞面积的乘积成正比.为了确定阻力的大小，用半径3m、载重300kg的降落伞从500m 高度作降落试验，测得各时刻的高度（表2）。

试确定降落伞的选购方案，即共需多少个，每个伞的半径多大（在表1中选择），在满足空投要求的条件下，使费用最低。

表二二、 模型的基本假设1 降落伞下降过程中只受重力及空气阻力的作用,其他因素忽略;2 降落伞的质量和绳索质量忽略;3 每个降落伞载的物重都不会超过降落过程中的最大载重;4 阻力系数是常数，与其他因素无关;5 降落伞的落地速度不会超过20m/s;6 救灾物资2000kg 可以任意分配;7 假设空投物资的瞬时伞已打开。

三、 符号的说明k阻力系数 r 伞半径f空气阻力r m 半径为r 的降落伞的最大载重 r s 半径为r 的降落伞的伞面面积 ()t H t 时刻降落伞的下降高度 ()t v t 时刻降落伞的下降速度 r n 购买半径为r 的降落伞数目 1C 伞面费 2C 绳索费 3C 固定费用l 降落伞每根绳索的长度g重力加速度,2/8.9sm g =iN 表示第i 种伞的个数,i=1,2,3,4,5四、 问题分析题目要求满足空投物资要求、费用最低等方面，制定一个降落伞合理选择的方案。

问题要求使总费用最小，而总费用受伞面费用、绳索费用和固定费用三部分影响，伞面费用和绳索费用受半径r 影响，而且半径r 的大小还影了伞的最大载重量。

另一方面，由于限制了最大速度，速度受空气阻力和载重量的影响，空气阻力与阻力系数k 有关，综上所述，我们需要解决的问题如下：1 通过对降落伞下降过程的受力分析，求出阻力系数k 值；2 确定每种伞的最大载重量m;3 求出每种伞的绳索费用；4 求出每种伞每把的总费用；5 根据要求建立模型并求解。

五、 模型的建立与求解5.1 阻力系数k 的确定降落伞下降时对降落伞进行受力分析有f mg ma -=,由于开始时不同时刻的加速度是不同的,即()t a a =,设初始时()00=v ,所以由上式有()⎪⎩⎪⎨⎧=-=00v kvs mg dt dvm即 ()⎪⎩⎪⎨⎧=-=00v m kvs g dtdv解上述微分方程,过程见附录一, 解得()()1mt s k eks mgks mgt v --=由物体下降高度H 和时间的二次微分等于加速度建立方程得：()()222/()/0000(0)0d H dt m g ksv m d H d H H =-===用MATLAB 解微分方程得：（程序见附录二）又当kg m m r 300,3==降落伞从m 500高空下降时其运动轨迹可由计算机模拟出来,采用MATLAB 的数据拟合函数并画图如下：该程序见附录三222222()k s t mmgt m g m g H t eksk sk s-=+-可得到其运动曲线如图所示图2可以发现降落伞在后期的运动曲线几乎是线性的,所以可以把降落伞后期的运动看成是匀速直线运动.对降落伞进行受力分析,有mg =,而kvs f =其中s 为降落伞的伞面面积. 取22,3,8.9,300rs m r g kg m π====，估算出s m v k /17,9.2≈≈由()mt s k eksmg ksmg t v --=把π⨯⨯=====232,17,9.2,8.9,300s v k g m 代入上式，可以用MATLAB 作出速度与时间的图象， 如下图：图3可以发现降落伞s 9以后速度几乎不变，这说明降落伞后期是作匀速直线运动的，所以降落伞后期匀速运动的速度可以这样确定:9秒以后的数据用最小二乘法进行线性拟合,设δ++=b vt t H )(,其中δ符合正态分布,采用Matlab 的polyfit 函数 程序见附录四，17.5794/v m s=因为降落伞为半球面,所以2222/4r r s ππ==由此解得 959.2925794.178.9300=⨯⨯⨯==πvsmg k5.2 降落伞载重的确定因为在高度一定时，物体在下落过程所用时间是不变的，因此可以把()v t 关系式中的t 看作一个定值，所以：由()1式可得m v 关于的函数()mt s k eksmg ksmg m v --=因为降落伞落地时,()500=t H , 即500222222=-+-sk g m esk g m ksmgtmt s k降落伞的最大载重当且仅当v 达到最大,即s m v /200=时取得,由此我们可以证明如下命题：()mt s k eksmg ksmg m v --=是关于m 的增函数.证明: 见附录五;反之m 关于v 的函数()v m 也是增函数. 令 ,500m H =,/20s m v =则()4式变为()()2222222220500k s tmk s t mm g m g v t e ks ks s r v H m gt m g m g H t e ks k sk sπ--⎧=-⎪⎪⎪=⎪=⎨⎪=⎪⎪=+-⎪⎩[]2.959,9.8,22.533.54k g r ===把2=r ,5.2=r ,=r ,.3=r ,4=r 代入上式可以解得满足空投条件下的各种半径的降落伞的最大载重： 运算程序见附录六;2151m kg=2.5237m kg=3341m kg=3.5464m kg=4607m kg=5.3 半径为4,5.3,3,5.2,2=r m时降落伞的绳索费因为载重m 位于球心正下方球面处,所以载重到球心的距离等于降落伞球面的半径,由几何关系得到2=l r ,其示意图如右图l图4 解得半径为4,5.3,3,5.2,2=r 的长度分别为123452.83, 3.53, 4.24, 4.95, 5.66l m l m l m l m l m =====求得绳索费如下表:5.4 确定降落伞的选购方案要使总费用最小,则要取每种半径的总成本最小,则用LINGO 可计算得:()12345151*n 237*n 341*n 464*n 607*n 2000;3++++≥N10, N20, N36, N40, N50;===== 4932c = 计算程序见附录七六、 模型结果的说明由上用LINGO 解得N10, N20, N36, N40, N50;===== 4932c = 即要购买降落伞的方案为:七、 模型的评价与推广7.1 优点：本模型对降落伞的运动作简化,即在后期,把降落伞的运动看作为匀速直线运动,则其运动方程为一次线形函数,对其求导,即可求得运动速度.从而可知运动速度与空投高度,运动时间无关.所以只要满足空投条件就可以降低空投高度,以减少空投难度. 7.1 缺点：因为在未知阻力系数k 和速度v 的情况下，根据题中所给数据和一些假设，估算出k 和v 的值，虽然其与后来算出值误差不大，但在一定程度上还是对后面的数值计算产生影响. 7.3 推广：当降落伞的半径仍为2m,2.5m,3m,3.5m,4m 五种时,其它条件不变,现在救灾物资很多,超过kg 2000,要求确定选购方案,则只需将(3)式的右端改为其它数据,如7000,8000等,就可求出相应的选购方案及总费用.改进：由于本模型假设的是在物资抛落的瞬时伞已打开，而在实际情况中，物资抛落后应有一段自由落体运动。

在模型的改进时应考虑到这一点，以便让模型更切合实际。

八、 参考文献[1] 姜启源．数学模型．北京：高等教育出版社，2008 [2] 薛翼．数值分析与实验．北京：北京工业大学出版社，2007[3] 赵静,但琦.数学建模与数学实验.高等教育出版社, 2008附录附录一：>> clear>> v=dsolve('m*Dv+k*s*v-m*g=0','v(0)=0','t')v =g/k/s*m-exp(-k*s/m*t)*g/k/s*m>> pretty(v)k s texp(- ------------) g mg m m--- - -----------------------------k s k s >>附录二>>H=dsolve('D2H*m+k*s*DH=m*g','H(0)=0,DH(0)=0','t')H =g/k^2/s^2*m^2*exp(-k*s/m*t)+g/k/s*m*t-1/k^2/s^2*m^2*g >> pretty(H)2 k s tg m exp(- -----) 2m g m t m g----------------- + ----- - -----2 2 k s 2 2k s k s >>附录三表一程序：>> x=[0 3 6 9 12 15 18 21 24 27 30];>> y=[500 470 425 372 317 264 215 160 108 55 1];>> plot(x,y)表二程序：>> m=300;>> g=9.8;>> k=2.9;>> v=17;>> s=2*9*pi;>> x=[0 3 6 9 12 15 18 21 24 27 30];>> V=m*g/(k*s)-m*g/(k*s)*exp(-k*s*x/m);>> plot(x,V)附录四>> clear>> X=(9:3:30);>> H=[128 183 236 285 340 392 445 499];>> P=polyfit(X,H,1)P =17.5794 -29.2976附录五命题：()m t s k e ks mg ks mg m v --=是关于m 的增函数.证明: 因为 ()()03222<-=--=---m ts k m t s k m t s k e m ks gt dmm v d e m gt e ks g ks g dmm dv m t s k m t s k e m gte ks gks g dm m dv ----=∴)(为严格单调减函数又因为 0)(lim =∞→dm m dv mm t s k m t s k e mgte ks gks gdm m dv ----=∴)(0> 由单调性的判别法, ()m t s k e ks mgks mg m v --=为m 的严格增函数所以命题成立附录六在matlab 中建立一个名位myfun 的m 文件,如下:function F=myfun(x)r=2.5;g=9.8;k=2.9458;s=2*pi*r^2;F=[x(1)^2*g/(k^2*s^2)*exp(-k*s*x(2)/x(1))+x(1)*g*x(2)/(k*s)-x(1)^2*g/(k^2*s^2)-500;g*x(1)/(k*s)-g*x(1)/(k*s)*exp(-k*s*x(2)/x(1))-20];在matlab 的commend window 中输入以下命令:>> clear>> x0=[1;1];>> options=optimset('Display','iter');>> x=fsolve(@myfun,x0,options)附录七min =447*n1+596*n2+822*n3+1177*n4+1563*n5;151*n1+237*n2+341*n3+464*n4+607*n5>2000;n1>=0;n2>=0;n3>=0;n4>=0;n5>=0;@gin(n1);@gin(n2);@gin(n3);@gin(n4);@gin(n5);运行结果：Global optimal solution found at iteration: 43Objective value: 4932.000Variable Value Reduced CostN1 0.000000 447.0000N2 0.000000 596.0000N3 6.000000 822.0000N4 0.000000 1177.000N5 0.000000 1563.000Row Slack or Surplus Dual Price1 4932.000 -1.0000002 46.00000 0.0000003 0.000000 0.0000004 0.000000 0.0000005 6.000000 0.0000006 0.000000 0.0000007 0.000000 0.000000。