数学模型作业-20150601-最终定稿

1、武汉理工大学数学模型作业题一：箱子的摆放策略某省内知名企业生产的产品用形状为长方体的箱子包装，使用叉车将这些箱子从生产车间运输至仓库。这些箱子叠放在叉车的正方形底板上，如下图所示，叉车置放箱子的底板是一个边长为1.1米的正方形。箱子的规格是统一的（所有箱子的长方形底面的尺寸相同）。通常在一次运输中，箱子像下图中这样横着放，或者竖着放。下图所示的便是一种可行的摆放方法，但不一定是最优的。现在这家企业需要你们帮助建立一个通用的优化模型，使得给定长方形箱子的长和宽之后，利用这个模型就能算出该如何摆放箱子（不需考虑箱子的高度，即只考虑摆放一层箱子），才能使得一次摆放的箱子数量最多。问题：如果不允许箱子

2、超出叉车底板（如上图所示情形），也不允许箱子相互重叠，建立一个优化模型，考虑如何摆放这些箱子，才能使摆放的箱子数量最多? 利用你们构建的模型，分别计算出对于下表中型号1、型号2和型号3的箱子，最多可以摆放多少个？该如何摆放？如果你们能画出摆放示意图，那么将有助于这家企业更快地理解你们的方法。长（米）宽（米）箱子型号10.30.24箱子型号20.60.4箱子型号30.30.21 问题分析根据题目的要求，实质上讨论的是建立优化模型解决同一规格货物的箱子在装载时对于二维平面的利用率的问题即属于二维平面优化问题，因为题目中明确要求不需考虑箱子的高度并且只考虑摆放一层箱子。在优化的过程中，要求箱子在堆放

3、时不允许超出叉车底板，也不允许出现箱子相互重叠的情况。对于这问题拟采用启发式的算法，通过循环嵌套的方法由外之内优化，直到内层的剩余矩形平面无法摆放小矩形为止。通过MATLAB编程，最后把优化的结果输出到Excel中。二模型假设1.箱子在堆放时互不叠压，且认为在堆放的过程中不发生挤压变形。2.箱子可以放在底板的任意位置，箱子的摆放必须与叉车四边平行或正交而不能斜放更不允许斜放在底板的角落避免箱子发生变形。3.考虑实际的装箱操作，箱子在堆放时先沿边缘放置，再底板中心填放的原则。4.模型在建立的过程中不考虑箱子的总重量是否在叉车的承重范围内，即只需考虑如何最大利用叉车底板的面积。5.箱子在堆放的过程

4、中不超出叉车的底板。三符号说明表1：符号说明符号解释说明L叉车底板的长W叉车底板的宽a箱子的长b箱子的宽i下L边上放置a的个数j下L边上放置b的个数n右W边上放置a的个数m右W边上放置b的个数k上L边上放置a的个数t上L边上放置b的个数p左W边上放置a的个数q左W边上放置b的个数Sum1每一次循环时返回的小矩形的个数Sum优化模型最终得到的小矩形的个数4 模型的建立与求解4.1模型的建立对于该问题的研究，常用的方法有启发式算法，贪婪算法，智能算法。参看文献1，学者陈端兵、黄文奇提出了一种贪婪算法，基于占角，占穴的概念对该问题进行了研究，该算法认为在每一次放置矩形时都要求是最优的，要求是合法的占

5、角动作，并且在此基础上要求穴度和贴边数是最大的，这实际上满足的是局部最优，但不一定就是整体最优。参看文献2，学者隋树林采用了一种启发式的算法对该问题进行了研究，本次采用的优化模型在文献2的基础上进行了改进和完善考虑了，非条件下的堆箱操作，同时也考虑了底板的边能被小矩形的边长整除情况下的堆箱操作。设在(,本次题目要求L=W，这属于特殊情况，采用该优化模型优化不失其一般性)的矩形内放入的小矩形块，要想获得尽量多的小矩形块就必须得利用a和b的各种组合方式，使得在L和W方向上的利用率尽可能高，这就要求在优化时必须在L和W边上同时对a，b进行组合优化。优化示意图如图1所示：图1：优化示意图如图1所示，设

6、在下L边上放置i个a边，则在该边上可放置的b的个数为（在MATLAB里面表示取整语句，在本文中均采用这样的方式表示取整）；同样设在右W边上放置n个a边，那么在该边上可放置的b边的个数为。依次原理，可利用上边的的个数确定该边上的个数，可利用左边的个数确定该边的个数，上面的排样的约束条件可表示为：（1）（2）（3）（4）（5）（6）（7）（8）（9）（10）其中条件（1）-（4）是保证边和边的利用率，即要保证边上所留的间隙不能再放入边（即较小边），其中条件（5）-（8）是保证在排样时不能出现小矩形之间重叠的情况，条件（9）-（10）分别说明了和，和之间的关系。但是上面的排法可能会形成中间没有矩形的

7、情况，如图2所示：图2：中间没有矩形的二维排样图对于中间没有矩形的情况同样可以采用该优化方法继续优化，相当于再一次调用该优化程序，形成循环嵌套，最终可得到修正后的排样图，如图3所示。图3：修正后的二维排样图4.2优化模型程序实现对于模型的求解采用循环嵌套的方式求解，图4是该模型的主程序的流程图。图4：主程序流程图二维优化函数在文献2的基础上做了改进，当且可以采用文献2的思想，但是当不满足：且该条件时，此时要对模型进行适当的修改，这时可以认为，基于此条件再进行循环嵌套最终得到这种情况的优化结果,出现这种情况即意味着二维优化即将完成，图5是二维优化函数（ewyh.m自定义的函数）的算法流程图。图5

8、：二维优化函数流程图下面将以一段伪代码说明上述的二维优化函数的算法：Step1：给sum0赋初值，sum0=0;计算底板的边长和小矩形的边相除的余数：，。判断和的值是否满足：是，执行Step2，否，执行Step3.Step2：是否满足：且；是，执行step2.1，否，判断是否满足：且；是，执行step2.2，否执行step2.3.Step2.1：；把此时的i，j，n，m，k，t，p，q的值分别赋给，x1，x2，x3，x4，x5，x6，x7，x8执行step4。Step2.2：；把此时的i，j，n，m，k，t，p，q的值分别赋给，x1，x2，x3，x4，x5，x6，x7，x8执行step4。St

9、ep2.3：for q=0:1:计算左W边剩余长度可放入a边的个数；for i=0:1:计算下L边的剩余长度可放入b边的个数；for n= :1:计算右W边剩余长度可放入b边的个数m；for :1:计算上L边剩余长度可放入b边的个数t；计算已经放置的小矩形的个数：；计算L1和W1的值：；（为了保 不出现矩形块重叠的情况）是否满足：的条件，是，执行Step2.4.否，跳出循环执行step4。Step2.4:判断sum2是否满足：的条件,是把sum2的值付给sum0即；，把此时的i，j，n，m，k，t，p，q的值分别赋给，x1，x2，x3，x4，x5，x6，x7，x8。否，跳出循环，记录当前的值x

10、1，x2，x3，x4，x5，x6，x7，x8，sum0执行step4。Step3：出现这种情情况即意味着模型的优化即将结束。for i=0:1:计算剩余长度堆放b的个数，在这种情况下，m，k，t，p的值均为0。for n=0:1:for q=0:1:计算已经放置的小矩形个数：；Step3.1：判断sum2是否满足：的条件,是把sum2的值付给sum0即；，把此时的i，j，n，m，k，t，p，q的值分别赋给，x1，x2，x3，x4，x5，x6，x7，x8。否，跳出循环，记录当前的值x1，x2，x3，x4，x5，x6，x7，x8，sum0执行step4。Step4：然后再对x1，x2，x3，x4，

11、x5，x6，x7，x8的值进行调整和修改，因为如果x1=0则x8=0，x2=0则x3=0，x4=0则x5=0，x6=0则x7=0，最后把计算的结果返回给主程序，程序结束。下面是一个实际的算例，当a=119mm，b=100mm，L=3000mm，W=1000mm优化的结果如下i=5;1，j=24;0，n=5;0，m=4;0，k=21;0，t=5;0，p=5;0，q=4;1，sum1=249;1，sum=250最终的优化示意图如图6所示：图6：算例二维排样图4.3模型求解针对本次作业的要求，有3种类型的箱子，下面将分别对三种类型的箱子进行优化，并得到对应的二维排样图。（1） 第一种箱子：，。调用该

12、优化程序得到：i=2，j=2，n=2，m=2，k=2，t=2，p=2，q=2，sum=16，根据以上优化结果可得到排样图如图7所示。图7：型号1箱子二维排样图（2） 第二种箱子：，。调用该优化程序得到：i=1，j=1，n=1，m=1，k=1，p=1，t=1，q=1，sum=4，根据以上优化结果可得到排样图如图8所示。图8：型号2箱子二维排样图（3） 第三种箱子：，。调用该优化程序得到：i=1，j=4，n=1，m=4，k=3，t=1，p=3，q=1，sum=20，根据以上优化结果可得到排样图如图9所示。图9：型号3箱子二维排样图5 结果分析根据上面的优化结果，可得到每一种型号的箱子在装箱后，对于

13、底板的面积利用率如表2所示：表2：结果分析型号排放个数利用率（%）1（a=0.3m，b=0.2m）1695.212（a=0.6m，b=0.4m）479.343（a=0.3m，b=0.2m）2099.17从表2可看出，对于型号1和型号3的底板的利用率都可以达到95%以上，理论上在底板上能放箱子的个数：，对于型号1的箱子，对于型号2的箱子，对于型号3的箱子可见型号1和型号3的优化的结果均达到了理论计算值（这也是为什么在二维优化函数中（ewyh.m）计算的值sum2要满足条件：），但上述理论值计算是基于不考虑实际的装箱操作，如果箱子可分解对于型号2的箱子当然可以装载5个，可见在考虑实际装箱和实际情况

14、型号2的箱子的个数已经达到最优解。6 模型评价表3：本文算法优化结果和文献2比较（长度单位：mm）底板尺寸箱子尺寸文献2计算结果利用率（%）本文计算结果利用率（%）1201,23360,3512896.0612896.061201,23360,509298.639298.631201,23370,3013198.3113198.311201,23380,506998.636998.631800,900590,297997.40997.40500,40053,3710098.059997.07表3是和文献2的算例结果做比较，从表中可看出，本文的优化结果和文献2的优化结果能很好的吻合，并且使得平面

15、的利用率达到95%左右。表4：本文算法优化结果和文献3比较（长度单位：mm）底板尺寸箱子尺寸文献3计算结果利用率（%）本文计算结果利用率（%）100,608,514798.0014898.67200,608,529899.33300100500,40052,389896.829997.813000,23360,3532196.4332196.433000,1000200,1808096.008096.003000,1000240,1309497.769497.763000,1000120,10024899.202501003000,1000467,2302693.082693.08表4是和文献

16、3的算例结果做比较，从表中可看出，本文的优化结果整体较文献3更优，并且使得平面的利用率达到95%左右，从表4可见对于第二种类型和第七种类型的箱子，本文的算法可使得利用率达到100%，可见本文算法的优越性，优化模型所用的MATLAB程序如附录所示。6.1模型的优点（1） 该方法计算简便，可操作性强。（2） 是解决目前集装箱利用率低问题的一种有效的方法。（3） 这种方法同样可应用于材料的剪裁问题（针对同一规格坯料）。（4） 该方法是基于整体最优的目的编程实现的这也是为什么在程序里会有这一条件，在循环嵌套的过程中每循环一次就记录该循环下的小矩形个数并和上一次循环下的小矩形个数做比较并把较大值记录下来

17、。6.2模型的缺点（1） 程序中若出现多个可行解，优化结果只能提供一种优化方案。（2） 模型的假设过于理想化。（3） 没法证明该方法是最优，但是算法的基本思想是基于整体最优。7 参考文献1 陈端兵，黄文奇.求解矩形paking问题的贪心算法 J.计算机工程，2007,33（4）：2 隋树林，邵巍，高自友.同一尺寸货物三维装箱问题的一种启发式算法 J.信息与控制，2005,34（4）：3 薛莲.同一规格货物集装箱在装载问题研究及其在物流行业的应用 D.天津：天津大学计算机科学与技术学院，2008.4 农健恒，崔耀东.同尺寸物品装箱的动态规划算法 J.计算机应用与软件，2014,31（7）：5 陈

18、德良，陈治亚.三维装箱问题的智能启发式算法 J.中南林业科技大学学报，2009,29（3）：附录：MATLAB编程程序主程序：clear allclca=input('请输入箱子的长度值a=');%长度和宽度一定不要混淆必须按照提示正确输入b=input('请输入箱子的宽度值b=');L=input('请输入底板的长度L=');W=input('请输入底板的宽度W=');r=0;sum=0;while L>=a&&W>=b%用于判断剩余矩形能否再放入小矩形块 r=r+1; i(r),j(r),n(r),

19、m(r),k(r),t(r),p(r),q(r),sum1(r)=ewyh(a,b,L,W); L=L-(j(r)+t(r)*b; W=W-(m(r)+q(r)*b; if L<W Z1=W;Z2=L;L=Z1;W=Z2; end sum=sum+sum1(r);%计算能放入的矩形的总个数endi=i'j=j'n=n'm=m'k=k't=t'p=p'q=q'sum1=sum1's=i,j,n,m,k,t,p,q,sum1;xlswrite('优化结果.xls',s,'A2:I3');%

20、把优化以后的结果输出到excel中，指定输出的范围为A2到I3供工人查看,表头是在excel里面手动输入.自定义函数：function x1,x2,x3,x4,x5,x6,x7,x8,sum0 = ewyh( a,b,L,W )sum0=0;y1=rem(L,a);y2=rem(L,b);y3=rem(W,a);y4=rem(W,b);%先确定L,W能不能被小矩形的边a,b整除，通过取余的方法来判断，这样能保证底板的利用率尽可能高.if L>=a+b&&W>=a+b if y1=0&&y4=0 i=fix(L/a);q=fix(W/b);j=0;n=

21、0;m=0;k=0;t=0;p=0; x1=i; x2=j; x3=n; x4=m; x5=k; x6=t; x7=p; x8=q; sum0=q*i+j*n+m*k+t*p; elseif y2=0&&y3=0 j=fix(L/b);n=fix(W/a);i=0;m=0;k=0;t=0;p=0;q=0; x1=i; x2=j; x3=n; x4=m; x5=k; x6=t; x7=p; x8=q; sum0=q*i+j*n+m*k+t*p; else for q=0:1:fix(W/b)%下面的循环程序是对i,j,n,m,k,t,p,q进行遍历然后根据条件跳出循环. p=fi

22、x(W-q*b)/a); for i=0:1:fix(L/a) j=fix(L-i*a)/b); for n=fix(q*b/a):1:fix(W/a) m=fix(W-n*a)/b); for k=fix(j*b/a):1:fix(L/a) t=fix(L-a*k)/b); sum2=q*i+j*n+m*k+t*p; W1=(m+q)*b; L1=(t+j)*b; if W1<W&&L1<L if sum2>sum0&&sum2<=fix(L*W/a/b)%用于判断目前计算的总和是否达到最大值，所以每循环一次就把sum1的值和上一次循环

23、计算的sum0做比较. sum0=sum2;%把较大值赋值给sum0. x1=i;%最后跳出循环时把当前的满足条件的i，j，n，m，k，t，p，q的值分别赋给x1,x2,x3,x4,x5,x6,x7,x8. x2=j; x3=n; x4=m; x5=k; x6=t; x7=p; x8=q; end end end end end end endelse for i=0:1:fix(L/a) j=fix(L-i*a)/b); m=0;%在这种情况下m、k、t、p均为0. k=0; t=0; p=0; for n=0:1:fix(W/a) for q=0:1:fix(W/b) sum2=i*q+j

24、*n+m*k+t*p; if sum2>=sum0&&sum2<=fix(L*W/a/b) sum0=sum2; x1=i; x2=j; x3=n; x4=m; x5=k; x6=t; x7=p; x8=q; end end end endendif x1=0|x8=0%下面的4个if语句是对x1,x2,x3,x4,x5,x6,x7,x8的值进行调整 x1=0;x8=0;endif x2=0|x3=0 x2=0;x3=0;endif x4=0|x5=0 x4=0;x5=0;endif x6=0|x7=0 x6=0;x7=0;end作业二：投资效益排序下表是我国198

25、42000 年宏观投资的一些数据，试利用主成分分析对投资效益进行分析和排序。表1 19842000年宏观投资效益主要指标年份投资效果系数（无时滞）投资效果系数（时滞一年）全社会固定资产交付使用率建设项目投产率基建房屋竣工率19840.710.490.410.510.4619850.400.490.440.570.5019860.550.560.480.530.4919870.620.930.380.530.4719880.450.420.410.540.4719890.360.370.460.540.4819900.550.680.420.540.4619910.620.900.380.560

26、.4619920.610.990.330.570.4319930.710.930.350.660.4419940.590.690.360.570.4819950.410.470.400.540.4819960.260.290.430.570.4819970.140.160.430.550.4719980.120.130.450.590.5419990.220.250.440.580.5220000.710.490.410.510.461 问题分析根据题目的要求，采用主成分分析对投资效益进行分析和排序，从而达到用较少的变量去解释原有资料中的大部分变异，将许多相关性很高的变量转化成彼此不相关的变量

27、，即所谓主成分，并用以解释资料的综合性指标。试图在力保数据信息丢失最少的原则下，对多变量的截面数据进行最佳综合简化，也就是说，对高维变量空间进行降维处理并计算综合得分最后根据综合得分进行排序和分析。2 主成分分析及排序本案例运用主成分分析方法综合评价我国1984-2000年宏观投资的一些数据，并对投资效益进行排序。（1） 对原始数据进行标准化处理假设进行主成分分析的指标变量有m个：，共有n个评价对象，第i个评价对象的第j个指标的取值为。将个项指标转换成标准化指标：其中，即为第j个指标的样本均值和样本标准差。对应地，称，（）为标准化指标变量。在处理过程先把原始数据写入sj.txt中并存在MATL

28、AB运行的当前文件夹中（即Current Folder里面），在本次案例分析中指标变量有5个即，有17个评价对象即。记在本案例中则有：通过数据标准化处理可得到：（2） 计算相关系数矩阵R相关系数矩阵，其中，式中是第个指标与第个指标的相关系数。通过MATLAB计算可得到相关系数矩阵为：（3） 计算特征值和特征向量计算相关系数矩阵R的特征值并按从大到小进行排序，及对应的特征向量，其中，由特征向量组成m个新的指标变量。其中是第一主成分，是第二主成分，是第m主成分。针对本次作业，通过MATLAB计算可得到特征值和相应的特征向量为：；。对应的特征向量构成的矩阵为:（4） 选择（）个主成分，计算综合评价值

29、计算特征值的信息贡献率和累积贡献率。为主成分的信息贡献率；为主成分的累积贡献率，理论上当时，则可选择前个指标变量作为个主成分，代替原有的的m个指标变量，从而可对个主成分进行综合分析。计算综合得分其中为第个主成分的信息贡献率，根据综合得分值就可以进行评价并对投资效益进行排序。本次案例的计算在后续的（5）中给出计算的结果。（5） 综合评价及投资效益排序表1：主成分分析结果序号特征根贡献率（%）累积贡献率（%）13.134362.686662.6921.168323.367086.0530.35027.003693.0640.22584.1516297.5750.12132.4266100从表1可以

30、看出，前3个特征根的累积贡献率就达到90%以上（超过了85%），主成分分析的效果很好，下面暂取前3个主成分进行综合评价，前3个特征根对应的特征向量如表2所示，但是累积贡献率超过85%作为选取的原则，这其实一种不严格的规定，近年关于这方面的研究和讨论很多，有很多比较成熟的方法，下面将讨论特征值因子的筛选方法。单纯的考虑累积贡献率有时是不够的需要做后续的判断，即还需考虑选择的主成分对原始变量的贡献值，这个问题可用相关系数的平方来表示，如果选取的主成分为，.，则选取的主成分对原始变量的贡献值为：这里表示与的相关系数1。特征向量 标准化变量第一特征向量0.49050.5254-0.48710.0671

31、-0.4916第二特征向量-0.29340.0490-0.28120.89810.1606表2：标准化变量的前3个特征值对应的特征向量第三特征向量0.51090.43370.37140.14770.6255利用MATLAB编程计算得到，。没出现的情况，可见选择的主成分都包含了原有变量的信息，并且累积贡献率达到了93.06%，能很好地解释原有资料中的信息。由此可得到三个主成分分别为：从主成分的系数可以看出，第一主成分主要反映了，投资效果系数（无时滞），投资效果系数（时滞一年）、全社会固定资产交付使用率、基建房屋竣工率4个指标的信息，第二主成分主要反映了建设项目投产率，第三主成分主要反映了，投资效

32、果系数（无时滞）、投资效果系数（时滞一年）、基建房屋竣工率3个指标的信息，但是3个主成分共同反映了该5个指标的信息。把每年对应的原始5个指标的标准化数据带入3个主成分的表达式，便可得到每一年的3个主成分的值如表3所示。表3：每一年对应的3个主成分的值年份第一主成分（）第二主成分（）第三主成分（）19840.6988-1.6351-0.04151985-1.05570.38200.61961986-0.9266-1.18351.085419871.5373-0.64520.46621988-0.2090-0.4387-0.46601989-1.3160-0.6019-0.084319900.60

33、31-0.67050.063419911.72060.05160.309719923.04300.5044-0.693919932.92162.55280.263719941.13410.57220.22481995-0.2757-0.2405-0.34791996-1.29960.5036-0.62071997-1.70180.0903-1.45541998-3.26301.41280.43561999-2.31390.98090.282820000.6988-1.6351-0.0415分别以三个主成分的贡献率为权重，构建出主成分综合评价模型：，把每一年的三个主成分的值带入该式，便可得到每一

34、年投资效益的综合评价值及排序结果，如表4所示。表4：排名及综合得分年份199319921991199419871990名次123456综合得分2.44641.97861.11230.86040.84560.2258年份198420001995198819851996名次789101112综合得分0.05310.0531-0.2534-0.2662-0.5292-0.7405年份19861989199719991998名次1314151617综合得分-0.7789-0.9715-1.1476-1.2015-1.68483 结果分析每一年宏观投资的效益存在较大的差异，1993、1992、1991、

35、1994、1987年的投资效益较好，主要表现在投资效果系数（无时滞），投资效果系数（时滞一年）等方面。1990、1984、2000年的投资效益也较好。1995、1988、1985、1996、1986年的投资效益处于中等。1989、1997、1999、1998的投资效益最差在今后的宏观投资时一定要避免出现这种情况。参看文献1和文献4，在计算的过程中对pcacov（R）语句求得的特征向量U做了修改，修改了其正负号，即让每个特征向量乘以该特征向量各元素和的符号函数值（和大于0函数值为1，和小于0函数值为-1），使得特征向量的分量和大于0。在交互式选择主成分的个数时，程序中不仅考虑了累积贡献率大于85

36、%还考虑了选取的主成分对于原始变量的贡献值，剔除了选择的主成分不包含原有变量的可能性，从而使得后续的分析和排序更加准确。同时也达到了用较少的变量去解释原始资料中的大部分变异，从而将许多相关性很高的变量转化成彼此相互独立或不相关的变量即主成分。在根据累积贡献率交互式选择主成分的个数时，考虑了选取的主成分对于原始变量的贡献值在用MATLAB实现时首先是求矩阵的相关系数矩阵（这时选择的是三个主成分），通过MATLAB得到的相关系数矩阵为：记，理论上有，由于MATLAB计算的误差问题值不为零，但可以看到它们的数量级均达到了-16次方近似为0，可以接受这样的计算误差，在第二节中（5）涉及的主成分对于原始变量的贡献值的计算有：图1是MATLAB计算的程序流程图：图1：主成分分析程序流程图四.参考文献1司守奎.数学建模算法与应用 M.北京：国防工业出版社，2011:2许淑娜，李长坡.