《数学建模-土地平整问题(共25页).doc》由会员分享，可在线阅读，更多相关《数学建模-土地平整问题(共25页).doc（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上土地平整问题摘要随着社会的快速发展，对于山区的土地整理显得格外重要。为了更好地合理利用土地，选择合适的土地整理地点，达到节约成本的目的。我们会建立模型，根据几个重要的制约因素，使得土石方量最小，而且成本最低，在实际的操作过程中更接近于科学实施。针对问题一：利用软件MATLAB中的绘图功能，对数据进行处理，得到这片土地的三维图形与等高线图。针对问题二：采用二重积分的意义，即曲顶柱体的体积。利用分割、取值、求和、取极限的思想，将山地分割成无数个小曲顶柱体，然后对这些小曲顶柱体进行求和，并结合土石方量最小，挖填土时费用使用的最少限制条件，从而求得在什么海拔高度进行平整。对于
2、平面位置的确定，利用枚举法，总费用最少，求得在什么位置，什么海拔高度平整，已达到最优的标准。针对问题三：类比于问题二，采用相同的方法，在不同的海拔高度上平整，已满足二层的台阶状地块的要求，达到挖填土石方量最小的目标。最后，对所建立的模型和求解方法的优缺点给出了客观的评价，并指出了改进的方法。关键词：土石方量最小；二重积分；枚举法；一、问题重述与分析1.1问题的重述 十堰市是一个山区城市，向山要地是十堰市发展的一个必然的选择，但是如何在一片山地之中选择合适的方位与开挖深度，从而使总的土石方量最小，就是一个有意义的命题了！某工厂为了在一片长度为1500米，宽度为1000米的山地之中，开挖出一个80
3、0米600米平坦连续的长方形地块作为工厂的厂房地基，前期已经在本块土地上测量出长、宽每隔30米的网格的对应网格点的海拔高度。问题： 用附件（1）中的数据画出工厂的这片土地的三维图形与等高线图； 从什么地方，什么海拔高度平整一块800米600米的连片土地能使总的土石方量最小？ 如果允许平整出来的土地为二层的台阶状地块，要求各地块的长、宽不少于60米，又将从什么地方、什么海拔高度分别开挖，能使总的土石方量最小？1.2问题的分析1.2.1问题一 我们可以利用MATLAB软件直接调用surf函数和contour函数分别画出这片土地的三维图形和等高线图。 1.2.2问题二 考虑要平整一块800米600米
4、的连片土地，且要使土石方量最小，我们需要确定这块800米600米的连片土地的底面投影区域，以及海拔高度，使得其土石方量最小。让底面投影区域在长度为1500米，宽度为900米内枚举，此时计算出对应的土石方量最小的是，在这些所有的最小值中取得最小值的那块投影区域即为所求。我们采用此种方法求解。在平整土地的过程中，凸出的地方是要挖土的，凹下的地方是可以填土的，这样挖下来的土可以填充到凹下去的部分。通过计算体积的方式计算挖土和填土大小，求出最小值对应的海拔高度，即可确定整连片土地的具体高度。二、基本假设1、在平整土地的过程中，有些地方是要挖山的，但有些地方是要填土的，假设填土的每立方米所需的费用为挖山
5、的每立方米土石方所需费用的1/3。2、假设除了挖土和填土以外，在平整土地的过程中其他作业（如登高）产生的费用都与800米600米的连片土地所处的位置方向和海拔高度均无关。这样我们只需考虑挖土和填土的土石方量及费用，据此来考虑连片土地所处的位置方向和海拔高度。3、假设计算体积的过程中，分割的小曲顶柱体不能达到无穷小，而产生的误差忽略不计。实际计算中我们取一个很小的步长去划分，使其划分尽可能的小，产生的误差忽略不计。三、符号约定符号说明平整连片土地块的长平整连片土地块的宽平整连片土地块的底面所在闭区域平整连片土地块底面划分的子区域子区域的面积平整土地所需的挖土量平整土地所需的填土量平整土地挖土和填
6、土所需要的总费用四、模型的建立4.1. 问题一求解1 利用MATLAB7.0软件的三维绘图功能，画出工厂这片土地的三维图形如图4.1.1所示：图4.1.1山地的三维图形2、 这片土地的三维等高线图如图4.1.2所示，二维等高线图4.1.3所示：图4.1.2山地的三维等高线图图4.1.3山地的二维等高线图4.2. 问题二模型建立4.2.1. 平整块的海拔高度的确定通过土石方量费用最小的原则，确定平整土地海拔的开挖高度。平整土地快所对应的的底面区域记为，显然是一个800米600米的矩形区域。设顶部海拔高度是一个非负连续函数，。假设 的位置已经确定，下面我们来讨论在什么样的海拔高度下，土石方量费用最
7、小。考虑底面区域对应的山体是一个曲顶柱体，把区域分成个小区域，为了简化运算，我们分成个等大的区域，每个区域都是边长为的正方形。再分别以这些小区域的边界为准线，作母线平行于轴的柱面，这些柱面将原曲顶柱体分成个细的曲顶柱体。每个区域上任取一点。当很大时，所分的区域很小，边长的值很小，每个区域对应的海拔高度z可以近似为一个相同的值，记所有小区域的最小海拔高度为，最大海拔高度为 （1） （2）设在平整块连片土地建在海拔高度为处。当小区域所对应的海拔高度时，说明这块区域需要挖土，挖土量为， （3）这里表示的面积。时，说明这块区域需要填土，填土量为 （4）总的挖土量， （5）其中为满足条件的小区域的个数。
8、总的填土量， （6）其中为满足条件的小区域的个数。当，挖土量大于填土量，这样挖出来的土足够填充凹下的山地。总的土石方量费用为： （7）当，挖土量小于填土量，这样挖出来的土不足以填充凹下的山地。需要从别的地方挖土来填充。总的土石方量费用为： （8）因此 （9）如果让平整块的底面投影区域取遍所有可能的值，再根据上述思想求解此底面投影区域对应土石方量费用最小的开挖海拔高度。按照最小的对应的底部和开挖海拔高度来开挖，所得平整的连片土地能使总的土石方量最小。4.2.2. 平整块的底面区域的确定将平整块区域投影到底面（平面），确定开挖方向，即要确定平整块的底面位置。显然底面是800米600米的矩形，假设矩
9、形的四个顶点分别为。山地是长度为1500米，宽度为900米的矩形区域，因此平整块底面投影矩形可以在1500900的区域范围内任意移动，四个顶点不能超出这个范围我们要取遍所有可能的位置，必须枚举所有的可能值。设矩形的边长为=600，=800，且各点坐标位置也已经确定，设各顶点的坐标分别为（，）、（，）、（，）、（，）。设边的中点为，边的中点为，那么唯一确定了线段。OxyADBCEF由中点公式，中点，的坐标分别可以表示为（，），（，），且的边长为=800,因此只要矩形的位置和大小确定，那么的大小和位置就唯一确定。反过来如果的位置和大小确定了（分别为短边的中点），那么矩形的位置和大小也唯一确定。下面
10、我们分别说明矩形四个顶点坐标，以及四个边的直线方程计算方法。 设=800，两点的坐标分别为（，）、（，），为了使计算不重复，保证的唯一性，我们不妨设 ，且当时, ,若EF所在直线的斜率存在，设为。当时，斜率存在，分以下四种情况讨论：(1)当,即平行于轴时，此时矩形如下图所示：xyOADBCEF 此时，四条边对应的直线方程分别为：，：，：，：， 四个顶点的坐标分别为：（，）, （，），（，）, （，） 矩形的充分必要条件是：(2)当时，xyOADBCEF 此时，四条边对应的直线方程分别为: ,即： : ,即：: 任意两条直线的交点即为顶点，、四个顶点的坐标分别由以下方程组解出：矩形的充分必要条件
11、是：(3)当，即垂直于轴时xyO 此时，四条边对应的直线方程分别为：，：，：，：， 四个顶点的坐标分别为：（，）, （，），（，）, （，） 矩形的充分必要条件是：(4)当，C 此时，四条边对应的直线方程分别为: ,即： : ,即：: 任意两天直线的交点即为顶点，、四个顶点的坐标分别由以下方程组解出：矩形的充分必要条件是：如果让在平面内取遍所有可能的值，再根据上述的思想求解可以确定平整块的底面投影区域的位置，再根据4.2.1的介绍求解海拔位置。为了让能取遍所有可能值，先在山地海底平面上任意取一点作为点，再根据点确定点的位置。这样只需考虑点取遍所有值的可能情况，根据EF的长度，确定点可能的有效位
12、置。平面区域内任意取一点， ，以为圆心，=800为半径作圆，显然F可能取得圆周上的所有点。Oxy为了避免重复，令 ，且当时, ,这样的取值为上图所示的圆弧上。这样线段就可以确定下来。4.3问题三模型的建立 我们可以继续使用问题二中的模型，但要对问题二中的模型结合问题进行修改。我们依然要确定平整块的海拔高度，而对平整块的底面区域的面积的求解我们依然采用问题二中的模型。 由于平整出来的土地为二层的台阶状地块，所以我们可以想象成两个曲柱体，这两个曲顶柱体所对应的高度不同，所以我们需要求出两个海拔高度。记所有小区域的最小海拔高度为，最大海拔高度为 （10） （11）设在平整块连片土地建在海拔高度为处，
13、为了求解方便，我们不妨设，若。当小区域所对应的海拔高度时，说明这块区域需要挖土，挖土量为。 （12）这里表示的面积。时，说明这块区域需要填土，填土量为。 （13）总的挖土量， （14）其中为满足条件的小区域的个数。总的填土量， （15）其中为满足条件的小区域的个数。当，挖土量大于填土量，这样挖出来的土足够填充凹下的山地。总的土石方量费用为： （16）当，挖土量小于填土量，这样挖出来的土不足以填充凹下的山地。需要从别的地方挖土来填充。总的土石方量费用为： （17）因此 （18）五、模型的求解5.1问题二的求解根据4.2.1和4.2.2介绍的方法，编写MATLAB程序求解。值得说明的是，在计算过程
14、中，我们根据已有的检测数据，用了插值的方法，多维插值拟合计算某一点的海拔高度，如线性插值，三次样条插值等，计算思想在数值分析中有说明，在这里不做细说，直接调用MATLAB中的多维插值函数interp2()求解即可。根据算法的基本思想，作出程序流程图见附录1，MATLAB程序见附录2.程序的运性结果为，海拔高度，即在以这样的构成的矩形区域为底面，在海拔米处平整一块800米600米的连片土地，所需要的土石方量费用最低。5.2问题三的求解程序的运性结果为，海拔高度，海拔高度，即在以这样的构成的矩形区域为底面，在海拔米处和在以构成的矩形区域为底面，在海拔65.76处平整一块土地所需要的土石方量费用最低
15、。六、模型的评价与推广6.1模型的优点本题中的模型使用了matlab软件进行了拟合，具有较高的可信性。同时模型中对于问题分了多种情况进行讨论，具有很高的参考价值。6.2模型的缺点本题中的模型求解过程比较繁琐，运算量较大，运算效率相对较低。6.3模型的推广 本题中的模型可以应用于其它土地平整问题中，只需将模型进行修改便可以推广。参考文献1陈纪修，於崇华，金路，数学分析第二版，M，北京：高等教育出版社，1999.2李庆扬，王能超，易大义，数值分析第五版，M，北京：清华大学出版社，2008.3韩明，张积林，李林，林杰等，数学建模案例，M，北京：同济大学出版社，2012.4侯进军，肖艳清，谭敏，高明柯
16、，数学建模方法与应用，M，南京：东南大学出版社，2012.5宋世德，郭满才，数学实验，M，北京：高等教育出版社，2002.6姜启源等，数学模型，M，北京：高等教育出版社，2003.7盛聚等，概率论与数理统计，M，北京：高等教育出版社，1983.8薛定宇，陈阳泉，高等应用数学问题的MATLAB求解，M，北京：清华大学出版社，2004。9何正风，MATLAB在数学方面的应用，M，北京：清华大学出版社，2012.%脚本文件drawpicture.m%用附件中的数据画出工厂的这片土地的三维图形与等高线x=0:30:900;y=0:30:1500;X,Y=meshgrid(x,y);Z=xlsread(
17、data,1,b3:az33);surf(X,Y,Z); %三维图xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);title(山地三维图);figure,contour(X,Y,Z); %二维等高线title(二维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);figure,contour3(Z); %三维等高线title(三维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);%函数文件getypos.m%功能：根据已知E点，及F点的x坐标，取得满足EF=800且yfye所对应的F点的坐标fun
18、ction yf=getypos(xe,ye,xf)%输入:E点的坐标(xe,ye),F点的x坐标xf%输出:F的y坐标yfyf=sqrt(800*800-(xe-xf)2) + ye;%函数文件isvalidpoint.m%功能：判断P是否是有效的点function ret = isvalidpoint( p )% p是向量，长度为2，分别存放x,y轴的坐标if (p(1) 0) & (p(1) 0) & (p(2) 0 A1 = 1 k;k -1; %AB和AD的直线方程系数 b1 = pe(1)+k*pe(2),(-sqrt(1+k*k)*width/2-pe(2)+k*pe(1); %
19、AB和AD的直线方程对应的常数项 pa = A1b1; A2 = 1 k; k (-1); %AB和BC的直线方程系数 b2 = (pe(1)+k*pe(2),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和BC的直线方程对应的常数项 pb = A2b2; A3 = 1 k; k -1; %CD和BC的直线方程系数 b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直线方程对应的常数项 pc = A3b3; A4 = 1 k;k -1; %CD和AD的直线方程系数 b4 = pf(1)+k*pf
20、(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直线方程对应的常数项 pd = A4b4; else if k 0 A1 = 1 k;k -1; %AB和BC的直线方程系数 b1 = pe(1)+k*pe(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和BC的直线方程对应的常数项 pb = A1b1; A2 = 1 k; k (-1); %AB和AD的直线方程系数 b2 = (pe(1)+k*pe(2),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和AD的直线方程对应的常数项 p
21、a = A2b2; A3 = 1 k; k -1; %CD和AD的直线方程系数 b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直线方程对应的常数项 pd = A3b3; A4 = 1 k;k -1; %CD和BC的直线方程系数 b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直线方程对应的常数项 pc = A4b4; end end endend A2 = 1 k; k (-1); %AB和AD的直线方程系数 b2 = (pe(1)+k*pe(2
22、),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和AD的直线方程对应的常数项 pa = A2b2; A3 = 1 k; k -1; %CD和AD的直线方程系数 b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直线方程对应的常数项 pd = A3b3; A4 = 1 k;k -1; %CD和BC的直线方程系数 b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直线方程对应的常数项 pc = A4b4; end end e
23、ndend%函数文件isinregion.m%功能：判断点pt是否在矩形区域内function ret = isinregion( pt , pe, pf)% pe, pf分别表示E,F点的坐标% pt表示待判断的点，是否在矩形区域内ret = 1;% 0表示不再区域内，1表示不再区域内long = 800;width = 600;x0 = pt(1);y0 = pt(2);if pe(2) = pf(2) %EF平行于x轴 if (x0 pf(1) | (y0 (pf(2) + width/2) ret = 0; endelse if x0 = pf(1) %EF平行于x轴 if (x0 (
24、pf(1) - long/2) | (y0 pf(2) ) ret = 0; end else k = (pf(2) - pe(2)/(pf(1) - pf(1); if k 0 ab = x0 + k*y0 - pe(1) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); ad = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); bc = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0)
25、 | (cd 0) ret =0; end if (ab 0) | (bc 0) | (cd 0) ret =0; end else if k 0 ab = x0 + k*y0 - pe(1) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); bc = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); ad = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0
26、; end if (ab 0) | (ad 0) | (cd 0) ret =0; end end end endendend ad = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); bc = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (bc 0) | (cd 0) ret =0; end else if k 0 ab = x0 + k*y0 - pe(1) -
27、 k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); bc = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); ad = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (ad 0) | (cd 0) ret =0; end end end endendend%函数文件getmaxminpoint.m%功能：取得矩形区域内的点对
28、应的横纵坐标的最大值和最小值function xmin, xmax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd )xs=pa(1), pb(1), pc(1), pd(1);ys=pa(2), pb(2), pc(2), pd(2);xmin = min(xs);xmax = max(xs);ymin = min(ys);ymax = max(ys);end%脚本文件main.m%计算并输出最小石方量对应的底面矩形坐标和海拔高度x=0:30:900;y=0:30:1500;long = 800;%平整块矩形区域的长width = 600;%平整块矩
29、形区域的宽X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);%disp(Z);dp=30;%步长dr=5;%分区域用dh=0.2;%海拔高度步长direct=zeros(4,2);%开挖方向，存放底面投影的四个顶点seaheight=0;%开挖海拔高度minmincost=0;%最小费用代价pexi=0;% E的横坐标while pexi 900 % 循环E的横坐标 peyi=0;% E的纵坐标 fprintf(pexi); disp(pexi); while peyi 1500 % 循环E的纵坐标 % E(pexi, peyi) pfxj=max(pex
30、i-long,0);%F的横坐标初始值 fprintf(peyi); disp(peyi); maxpfxj=min(pexi + long),900);%小于900范围内取值 while (pfxj maxpfxj) % 循环F的横坐标，穷举F的所有可能情况，由距离确定F的纵坐标 pfyj=getypos(pexi, peyi, pfxj);%F的纵向坐标 % F(pfxj, pfyj) pe = pexi peyi;%E点的坐标 ret=isvalidpoint(pe);%F判断是否是有效范围内的点 if ret = 0 pfxj = pfxj + dp; continue; end %取
31、得ABCD四点坐标，并判断是否超出总的范围 pf = pfxj, pfyj;%F点的坐标 pa, pb, pc, pd = getpoints( pe, pf );%取得ABCD四点坐标 %至少有一个点超出范围 if (isvalidpoint(pa) = 0) | (isvalidpoint(pb) = 0) | (isvalidpoint(pc) = 0) | (isvalidpoint(pd) = 0) pfxj = pfxj + dp; continue; end % 固定矩形区域内求土石方量费用最小 % mincost = 0;%最小代价 height = 0;%最小代价，对应的海拔
32、高度 MP=zeros(ceil(long+1)*(width+1)/dp/dp), 3);%分小区域数long*width/d/d，记录小区域的任意一点 xmin, xmax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd );%找到区域内的最高点和最低点的横纵坐标 px = xmin; index=1; while px xmax %fprintf(px); %disp(px); py = ymin; while py ymax %判断点P(px,py)是否在矩形区域内 在此写一个函数去判断.预留 pt = px, py; ret = isinre
33、gion( pt , pe, pf); if ret=0 py = py + dr; continue; end %如果该点在矩形区域内，记录该点.为了节约时间，以空间换时间。 MP(index, 1) = px; MP(index, 2) = py; MP(index, 3) = interp2(x, y, Z, px, py);%多维插值，取值 index = index + 1; py = py + dr; end px = px + dr; end %如果在区域内在此区域内求出体积最大值 mh = MP(:,3); minheight = min(mh); maxheight = max(mh); zh = minheight; while zh maxheight % 循环所有高度（从最低点的高度到最高点的高度） S = 0;% 记录挖土量 T = 0;% 记录填土量 index=1; len = size(mh, 1); while (index 0 S = S + temp*dr*dr;% 挖土量 else