数学建模——合理开挖土地问题(附MATLAB源程序)

1、本页只是说明，论文从第二页开始，下载后请删除本页即可:论文内容：关于合理开挖土地问题的数学建模竞赛论文（含Matlab源程序）特别申明：本论文版权归百度文库账号dxzsk同学所有，仅限个人下载学习使用，其他人不得转载分享，侵权必究。以下是本论文原始题目：合理开挖土地问题：A市是一个山区城市,向山要地是A市发展的一个必然的选择,但是如何在一片山地之中选择合适的方位与开挖深度,从而使总的土石方量最小,就是一个十分有意义的课题.A市某工厂为了在一片长度为1500米,宽度为900米的山地之中,开挖出一个800米600米平坦连续的长方形地块作为工厂的厂房地基,前期已经在本块土地上测量出长、宽每隔30米的

2、网格的对应网格点的海拔高度（详细数据见附件）.请你考虑以下几个问题：问题（1）：用附件中的数据画出工厂的这片土地的三维图形与等高线图；问题（2）：从什么地方,什么海拔高度平整一块800米600米的连片土地能使总的土石方量最小?问题（3）：如果允许平整出来的土地为二层的台阶状地块,要求各地块的长、宽不少于60米,又将从什么地方、什么海拔高度分别开挖,能使总的土石方量最小?提示：在平整土地的过程中,有些地方是要挖山的,但有些地方是要填土的,假设填土的每立方米所需的费用为挖山的每立方米土石方所需费用的1/3.2013*大学金水节第五届研究生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的

3、竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3.

4、指导教师或指导教师组负责人 (打印并签名)： 日期： 2013 年 11 月 01 日赛区评阅编号（由赛区组委会评阅前进行编号）：2013年*大学金水节第五届研究生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）： 关于合理开挖土地问题的数学模型求解摘要随着山区城市的发展，向山要地开发成为了一种必然的选择。为了更好的确定平整地的起始点，做到节约成本的目的，我们根据几个决定性因素利用网格优化建立了相关的模型。主要讨论如何在一片山地之中选

5、择合适的方位与开挖深度，从而使总的土石方量最小，使得在实际工程中达到最可靠，最节约成本。采用MATLAB软件模拟功能，画出了该工厂的这片土地的三维图形与等高线图，并对山地地形作出分析。采用二重积分的集合意义，即二重积分是曲顶柱体的体积，利用这个思想，进行合理假设，对曲顶柱体进行拆分，分成多个细小的小区顶柱体，再对这些小曲顶柱体的体积求和，并考虑挖土，和填土的费用代价，求出土石方量最小时的海拔高度。在不同的底面投影位置下，分别求出对应的费用消耗，求出最小费用，从而决定从什么地方，什么海拔高度平整一块固定大小的长方形连片土地，能使总的土石方量最小，使工程施工经济上达到最优。关键词： 合理开挖土地

6、最优 土石方量最小一、问题提出及其分析1. 提出问题A市是一个山区城市，向山要地是A市发展的一个必然的选择，但是如何在一片山地之中选择合适的方位与开挖深度，从而使总的土石方量最小，就是一个十分有意义的课题.A市某工厂为了在一片长度为1500米，宽度为900米的山地之中，开挖出一个800米600米平坦连续的长方形地块作为工厂的厂房地基，前期已经在本块土地上测量出长、宽每隔30米的网格的对应网格点的海拔高度（详细数据见附件）。请你考虑以下几个问题：问题一：用附件中的数据画出工厂的这片土地的三维图形与等高线图；问题二：从什么地方，什么海拔高度平整一块800米600米的连片土地能使总的土石方量最小？2

7、. 问题分析1) 问题一：我们可以利用MATLAB软件直接调用surf函数和contour函数分别 画出这片土地的三维图形和等高线图。2) 问题二：考虑要平整一块800米600米的连片土地，且要使土石方量最小， 我们需要确定这块800米600米的连片土地的底面投影区域，以及海拔高度，使得其土石方量最小。 连片土地底面投影区域的确定：分析一：对于指定底面是800米600米大小区域内，取海拔的高度做方差，如果方差最小，说明海拔都比较接近平均值，就是说明这片区域比较平坦，挖土或填土的土石方量近似最小。我们可以通过这个思路解决问题，即为求出满足海拔方差值最小的底面位置的问题。分析二：让底面投影区域在长

8、度为1500米，宽度为900米内枚举，此时计算出对应的土石方量最小的是，在这些所有的最小值中取得最小值的那块投影区域即为所求。我们采用此种方法求解。 连片土地海拔高度的确定：在平整土地的过程中，凸出的地方是要挖土的，凹下的地方是可以填土的，这样挖下来的土可以填充到凹下去的部分。通过计算体积的方式计算挖土和填土大小，求出最小值对应的海拔高度，即可确定平整连片土地的具体高度。二、问题的假设1、在平整土地的过程中，有些地方是要挖山的，但有些地方是要填土的，假设填土的每立方米所需的费用为挖山的每立方米土石方所需费用的1/3。2、假设除了挖土和填土以外，在平整土地的过程中其他作业（如登高）产生的费用都与

9、800米600米的连片土地所处的位置方向和海拔高度均无关。这样我们只需考虑挖土和填土的土石方量及费用，据此来考虑连片土地所处的位置方向和海拔高度。3、假设计算体积的过程中，分割的小曲顶柱体不能达到无穷小，而产生的误差忽略不计。实际计算中我们取一个很小的步长去划分，使其划分尽可能的小，产生的误差忽略不计。三、 符号说明符号说明L平整连片土地块的长W平整连片土地块的宽D平整连片土地块的底面所在闭区域平整连片土地块底面划分的子区域子区域的面积平整土地所需的挖土量平整土地所需的填土量C平整土地挖土和填土所需要的总费用四、模型的建立与求解4.1. 问题一求解1利用MATLAB软件的三维绘图功能，画出工厂

10、这片土地的三维图形如下：图表4.1.1山地的三维图形2，这片土地的等高线图如下：图表4.1.2山地的三维等高线图图表4.1.3山地的二维等高线图4.2. 问题二模型建立与求解4.2.1. 平整块的海拔高度的确定通过土石方量费用最小的原则，确定平整土地海拔的开挖高度。平整土地快所对应的的底面区域记为D，显然D是一个800米600米的矩形区域。设顶部海拔高度是一个非负连续函数，。假设 D的位置已经确定，下面我们来讨论在什么样的海拔高度下，土石方量费用最小。考虑底面区域D对应的山体是一个曲顶柱体，把区域D分成个小区域，为了简化运算，我们分成个等大的区域个等大的区域，每个区域都是边长为的正方形。再分别

11、以这些小区域的边界为准线，作母线平行于轴的柱面，这些柱面将原曲顶柱体分成个细的曲顶柱体。每个区域上任取一点。当很大时，所分的区域很小，边长的值很小，每个区域对应的海拔高度z可以近似为一个相同的值，记所有小区域的最小海拔高度为，最大海拔高度为设在平整块连片土地建在海拔高度为处。当小区域所对应的海拔高度时，说明这块区域需要挖土，挖土量为，这里表示的面积。时，说明这块区域需要填土，填土量为，总的挖土量，其中为满足条件的小区域的个数。总的填土量， ，其中为满足条件的小区域的个数。当，挖土量大于填土量，这样挖出来的土足够填充凹下的山地。总的土石方量费用为：当，挖土量小于填土量，这样挖出来的土不足以填充凹

12、下的山地。需要从别的地方挖土来填充。总的土石方量费用为：因此如果让平整块的底面投影区域取遍所有可能的值，再根据上述思想求解此底面投影区域对应土石方量费用最小的开挖海拔高度。按照最小的对应的底部和开挖海拔高度来开挖，所得平整的连片土地能使总的土石方量最小。4.2.2. 平整块的底面投影区域的确定将平整块区域投影到底面（XOY平面），确定开挖方向，即要确定平整块的底面位置。显然底面是800米600米的矩形，假设矩形的四个顶点分别为A、B、C、D。山地是长度为1500米，宽度为900米的矩形区域，因此平整块底面投影矩形ABCD可以在1500900的区域范围内任意移动，四个顶点不能超出这个范围我们要取

13、遍所有可能的位置，必须枚举所有的可能值。设矩形ABCD的边长为AB=W=600，AC=L=800，且各点坐标位置也已经确定，设各顶点的坐标分别为A（，）、B（，）、C（，）、D（，）。设AB边的中点为E，CD边的中点为F，那么唯一确定了线段EF。OxyADBCEF由中点公式，中点E，F的坐标分别可以表示为E（，），F（，），且EF的边长为EF=AC=800,因此只要矩形ABCD的位置和大小确定，那么EF的大小和位置就唯一确定。反过来如果EF的位置和大小确定了（E,F分别为短边的中点），那么矩形ABCD的位置和大小也唯一确定。下面我们分别说明矩形四个顶点坐标，以及四个边的直线方程计算方法。设EF

14、=800，E、F两点的坐标分别为E（，）、F（，），为了使计算不重复，保证EF的唯一性，我们不妨设 ，且当时, ,若EF所在直线的斜率存在，设为。当时，斜率存在，分以下四种情况讨论：(1)当,即EF平行于y轴时此时矩形如下图所示：xyOADBCEF此时，四条边对应的直线方程分别为AB：，BC：，CD：，AD：，四个顶点的坐标分别为：A（，）, B（，），C（，）, D（，）矩形ABCD的充分必要条件是：(2)当时，xyOADBCEF此时，四条边对应的直线方程分别为AB: ,即： CD: ,即：AD：BC: 任意两条直线的交点即为顶点，A、B、C、D四个顶点的坐标分别由以下方程组解出：矩形ABC

15、D的充分必要条件是：(3)当，即EF垂直于x轴时xyOADBCEF此时，四条边对应的直线方程分别为AB：，BC：，CD：，AD：，四个顶点的坐标分别为：A（，）, B（，），C（，）, D（，）矩形ABCD的充分必要条件是：(4)当，xyODABCEF此时，四条边对应的直线方程分别为AB: ,即： CD: ,即：BC：AD: 任意两天直线的交点即为顶点，A、B、C、D四个顶点的坐标分别由以下方程组解出：矩形ABCD的充分必要条件是：如果让EF在平面内取遍所有可能的值，再根据上述的思想求解可以确定平整块的底面投影区域的位置，再根据4.2.1的介绍求解海拔位置。为了让EF能取遍所有可能值，先在山地

16、海底平面上任意取一点作为E点，再根据E点确定F点的位置。这样只需考虑E点取遍所有值的可能情况，根据EF的长度，确定F点可能的有效位置。平面区域内任意取一点E， ，以E为圆心，R=L=800为半径作圆，显然F可能取得圆周上的所有点。OxyEFnF1F3F2为了避免重复，令 ，且当时, ,这样F的取值为上图所示的圆弧上。这样线段EF就可以确定下来。4.2.3. 编程求解具体开挖位置根据4.2.1和4.2.2介绍的方法，编写MATLAB程序求解。值得说明的是，在计算过程中，我们根据已有的检测数据，用了插值的方法，多维插值拟合计算某一点的海拔高度，如线性插值，三次样条插值等，计算思想在数值分析中有说明

17、，在这里不做细说，直接调用MATLAB中的多维插值函数interp2()求解即可。根据算法的基本思想，作出程序流程图见附录1，MATLAB程序见附录2.程序的运性结果为，A，B，C，C，海拔高度，即在以这样的ABCD构成的矩形区域为底面，在海拔米处平整一块800米600米的连片土地，所所需要的土石方量费用最低。五、模型的评价与推广模型最大的优点在于对原始数据进行拟合时，使用MATLAB软件描述山地的三维图形和等高线图形，在此基础上可以得到具有很高的拟合度和适度性。对模型的进一步讨论便可以得到一系列可靠而实用的信息，所得结论与客观事实很好地吻合，从而进一步说明模型是合理的。在解决合理开挖土地的问

18、题中，通过土石方量费用最小的原则，确定平整土地海拔的开挖高度。模型的缺点在于在建立模型时，由于要提高精度，必须减小循环迭代步长，从而使循环次数增大，增大了程序的运算量，运行数据量比较大，会消耗大量的运算时间，所以造成程序运行过程中结果出来比较慢，算法效率比较低。本模型虽然讨论的是以通过土石方量费用最小的原则来确定平整土地海拔的开挖高度，但是可以推广到各种合理开挖土地的问题中，用类似的方法进行建模求解。解决该问题时，我们首先利用MATLAB软件的三维绘图功能，画出工厂这片土地的三维图形，再确定平整块开挖的方向。本文所建立的模型可以推广到其他各种模型。参考文献1 陈纪修，於崇华，金路，数学分析 第

19、二版，M，北京：高等教育出版社，1999. 2 李庆扬，王能超，易大义，数值分析 第五版，M，北京：清华大学出版社，2008. 3 韩明，张积林，李林，林杰等，数学建模案例，M，北京：同济大学出版社，2012. 4 侯进军，肖艳清，谭敏，高明柯，数学建模方法与应用，M，南京：东南大学出版社，2012. 5 宋世德，郭满才，数学实验，M，北京：高等教育出版社，2002. 6 姜启源等，数学模型，M，北京：高等教育出版社，2003.7 盛聚等，概率论与数理统计，M，北京：高等教育出版社，1983.8 薛定宇，陈阳泉，高等应用数学问题的MATLAB求解，M，北京：清华大学出版社，2004。9 何正风

20、，MATLAB在数学方面的应用，M，北京：清华大学出版社，2012.附录文件1 流程图问题二解决方案的程序流程图如下图所示：2. 求解过程中的程序（已经在MATLAB Version 7.10.0.499(R2010a)环境下调试通过）。%脚本文件drawpicture.m%用附件中的数据画出工厂的这片土地的三维图形与等高线x=0:30:900;y=0:30:1500;X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);surf(X,Y,Z); %三维图xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);title(山地三维图);fig

21、ure,contour(X,Y,Z); %二维等高线title(二维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);figure,contour3(Z); %三维等高线title(三维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);%函数文件getypos.m%功能：根据已知E点，及F点的x坐标，取得满足EF=800且yfye所对应的F点的坐标function yf=getypos(xe,ye,xf)%输入:E点的坐标(xe,ye),F点的x坐标xf%输出:F的y坐标yfyf=sqrt(800*800-(xe-xf)2

22、) + 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); %AB和AD的直线方程对应的常数项 pa = A1b1; A2 = 1 k; k (-1); %AB和BC的直线方程系数 b2 = (pe(1)+k*pe(2),(sqrt(

23、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(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直线方程对应的常数项 pd = A4b4; else if k 0 A1 =

24、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的直线方程对应的常数项 pa = A2b2; A3 = 1 k; k -1; %CD和AD的直线方程系数 b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*

25、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),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和AD的直线方程对应的常数项 pa = A2b2; A3 = 1 k; k -1; %CD和

26、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%函数文件isinregion.m%功能：判断点pt是否在矩形区域内function ret = isinregion( pt , pe, pf)% pe, pf

27、分别表示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 (pf(1) - long/2) | (y0 pf(2) ) ret = 0; end else k = (pf(2) - pe(2)/(pf(1) - pf(1); if

28、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) | (cd 0) ret =0; end if (ab 0) | (bc 0) | (cd 0) ret =0; end else if k 0 ab = x0 + k*

29、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; end if (ab 0) | (ad 0) | (cd 0) ret =0; end end end endendend ad = -k*x0 + y0 - sqrt

30、(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) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); bc = -k*x0 + y0 - sqrt(1 + k2)*width/2

31、 - 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%功能：取得矩形区域内的点对应的横纵坐标的最大值和最小值function xmin, xmax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd )xs=pa(

32、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;%平整块矩形区域的宽X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);%disp(Z);dp=30;%步长dr=5;%分区域用dh=0.2;%海

33、拔高度步长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(pexi-long,0);%F的横坐标初始值 fprintf(peyi); disp(peyi); maxpfxj=min(pexi + long),900);%小于900范围内

34、取值 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 %取得ABCD四点坐标，并判断是否超出总的范围 pf = pfxj, pfyj;%F点的坐标 pa, pb, pc, pd = getpoints( pe, pf );%取得A

35、BCD四点坐标 %至少有一个点超出范围 if (isvalidpoint(pa) = 0) | (isvalidpoint(pb) = 0) | (isvalidpoint(pc) = 0) | (isvalidpoint(pd) = 0) pfxj = pfxj + dp; continue; end % 固定矩形区域内求土石方量费用最小 % mincost = 0;%最小代价 height = 0;%最小代价，对应的海拔高度 MP=zeros(ceil(long+1)*(width+1)/dp/dp), 3);%分小区域数long*width/d/d，记录小区域的任意一点 xmin, xm

36、ax, 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 = isinregion( pt , pe, pf); if ret=0 py = py + dr; continue; end %如果该点在矩形区域内，记录该点.为了节约时间，以空间换时

37、间。 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

38、;% 记录填土量 index=1; len = size(mh, 1); while (index 0 S = S + temp*dr*dr;% 挖土量 else T = T - temp*dr*dr;% 填土量 end index = index + 1; end if (S T) cost = 4*T/3; else cost = T/3 + S; end if (cost mincost) | (mincost = 0) mincost = cost; height = zh; end zh = zh+dh; end % 固定矩形区域 % if (minmincost mincost) | (minmincost = 0) minmincost = mincost; seaheight = height; direct = pa, pb, pc, pd; end pfxj = pfxj + dp; end peyi = peyi + dp; end pexi = pexi + dp;endfprintf(the final direct);disp(direct);fprintf(the final seaheight);disp(seaheight); %至少有一个点