基于曲面拟合与二重积分的沙矿投资利益分析模型Word格式.docx

1、蓝色隧道公司（B.T.I）提出，只要该矿含沙量超过400万吨，并且表土体积小于沙子体积19%，则开采一定获利。为该公司判断是否应该对这项工程投资。问题三：B.T.I公司计算出每月可能开采10万吨沙子，并有可能实施。根据以下数据给出收益分析。一些有关的数据：表土平均密度为1350千克/立方米；沙子的平均密度1620千克/立方米；月贴现率1%；沙子的开采率最大为100000吨/月。庄园主的索赔清单：补偿对现有农作物及房屋的损坏，一次性需付200000镑；开采一吨沙子付给地区使用费0.5472镑；倾倒表土要占用土地，每倾倒一吨表土付0.121镑；沙子开采完后，该地区的修建和美化付1000000镑。B

2、.T.I公司开采沙子的费用：将表土挖出并运到附近每吨需0.3672镑；开采一吨沙子并将其运到市场需0.4592镑。B.T.I公司所的利润：沙子的销售单价2.96镑/吨。2 模型假设1.假设用于探查结果的各网点之间经过插值之后的距离相较于整个开采地区而言足够小，可视为连续的点；2.假设表土及沙子在整个开采区内的密度均匀，都接近平均密度；3.假设在将表土移至别处的过程中，能恰好将其移除干净且带走的沙子可以忽略不计；4.假设相邻网格点之间的距离相等，即都为50米；3符号说明表1 符号说明及解释符号解释与说明z1z2z3A1 地表面的散点数据沙子顶面的散点数据沙子底面的散点数据沙子的开挖费A2 表面土

3、的开挖费B1 沙子的运输费B2 表面土的处理费C1 沙子的销售单价D1 沙层的厚度D2 表面土的厚度E1 沙子的密度E2表面土的密度W沙土总重量W1 表面土重量W2沙子重量p 开采每吨沙土所获利润q月贴现率4 问题的分析对于问题一，表1所给数据有部分因为沼泽地无法取样而缺失。由于必须有该开采地的全部三层数据才能够实现绘制等高线，因此考虑估算出未能取样到的沼泽地部分的三层数据。运用题中表1所给数据中关于开采地中三层均能够取样到的部分，即B列到M列、1行到7行的数据，采取拟合的方法求出开采地各层的曲面方程，并据此估算出其余未能够取样的地区的地表高度、沙层顶部的高度、沙层底部的高度。进而将三层数据分

4、别抽离，在matlab中运用contour语句可实现绘出地表、沙层顶部、沙层低部的等高线。对于问题二，判断是否应该对该项工程投资。该矿含沙量超过400万吨及表土体积小于沙子体积19%即可进行投资。在问题一中已拟合出开采地各层的曲面方程，接下来运用二重积分求体积的方法求出表土体积和沙子体积，即可算出两者体积值的关系。求出沙子体积，并已知沙子的平均密度，接下来运用公式即可算出该矿含沙量。对于问题三，由于每月开采量的限制，和月贴现率的设定，越晚赚得的利益，贴现到现在则越少。根据这个事实我们容易得出，按照该地区不同位置沙子体积与表土体积的比例，沙子所占比例越高的地方应最优先开采。从而我们制定了施工原则

5、、设计了利益函数，对相关的算法设计其程序，不难得到各个月份的收益分析，最终确定呢整个项目投资可收回的利润。5 模型建立与求解根据问题的分析，针对所要解决的几个问题，建立相应的数学模型，并采取相应方法求解，最终得到所有问题在一定条件下的最优可行解。5.1问题一的模型建立与问题求解5.1.1 模型的建立数据的取舍：首先考虑到整个开采区有沼泽地，从而导致少量某些网格点无数据。针对这一问题，可以从已有数据出发，取出尽可能大范围的mn的关于自变量x、y的数据矩阵，然后用matlab将其描点作图。通过观察整个开采地网格点数据的分布，不难看出以第17行与第213列组成的数据矩阵范围最大。因此选取其为描点作图

6、的样本数据。数据的处理：由于各网格点之间的距离相差较大，为了在得到光滑的曲面时误差尽量的小。首先对每个网格进行插值：将单位长度为50米的网格线细分为50份，变成单位长度为1米的新的网格点距离。对于开采地平面的坐标建立及单位的确定，同时考虑到坐标建立的简洁性和确保结果的准确性。故采取以网格平面的左上角为坐标原点，取x=1,2,9为行，y=1,2,15为列。,坐标轴的单位长度之所以取1而非50，一方面为了数据处理的简单，另一方面并不影响由已知数据到未知数据的回归预测。接下来将所取数据按照地表、沙层顶部和沙层底部分为3组。对每一组数据分别用matlab中的SFtool工具箱进行拟合测试，得出符合条件

7、最好的拟合结果。根据得出的函数表达式，进行沼泽地网格点的缺失数据的预测。根据预测的结果，对原数据所形成的函数进行检验，根据检验结果得出该模型的可行性。5.1.2 模型求解以第一组地表数据为例，其余两组数据步骤完全相同。对全部数据进行取舍之后的结果如表2所示：表2 地表高度中心部分取样234567891011121322.422.522.62323.123.223.423.52423.822.722.823.324.224.122.923.623.923.724.522.322.222.1将上述坐标值运用matlab的SFtool工具箱，在Polynomial选项下将x与y均设置为最高次数5次，

8、以保证拟合的精度。然后即可得到拟合的三维图像以及对应的函数表达式。对应的等高线也可随即得到。将地表面的数据输入工具箱后得到如下结果：p00 + p10*x + p01*y + p20*x2 + p11*x*y + p02*y2 + p30*x3 + p21*x2*y + p12*x*y2 + p03*y3 + p40*x4 + p31*x3*y + p22*x2*y2 + p13*x*y3 + p04*y4 + p50*x5 + p41*x4*y + p32*x3*y2 + p23*x2*y3 + p14*x*y4 + p05*y5 p00 = 5.5 （-3.234, 14.23） p10

9、= 1.097 （-1.367, 3.562） p01 = 19.79 （10.32, 29.27） p20 = 0.1297 （-0.4886, 0.748） p11 = -1.402 （-2.396, -0.4082） p02 = -8.415 （-12.56, -4.269） p30 = -0.03693 （-0.1203, 0.04642） p21 = 0.08438 （-0.02146, 0.1902） p12 = 0.3875 （0.1553, 0.6196） p03 = 1.697 （0.8183, 2.575） p40 = 0.003434 （-0.002139, 0.0090

10、07） p31 = -0.004379 （-0.01124, 0.002483） p22 = -0.01212 （-0.02386, -0.0003806） p13 = -0.04361 （-0.07074, -0.01648） p04 = -0.1633 （-0.2534, -0.0733） p50 = -0.0001188 （-0.0002646, 2.703e-05） p41 = 0.0002014 （1.533e-06, 0.0004012） p32 = -0.0001728 （-0.0005017, 0.0001561） p23 = 0.001288 （0.0006834, 0.00

11、1892） p14 = 0.00111 （-0.0001654, 0.002385） p05 = 0.006215 （0.002639, 0.009792）Goodness of fit: SSE: 1.66 R-square: 0.9446 Adjusted R-square: 0.9271 RMSE: 0.1623 根据工具箱给出的结果，拟合的地面的三维曲线如图1所示。图1 地面的三维曲线地表曲面函数方程为： f（z1）= 5.5+1.097*x+ 19.79*y+0.1297*x2-1.402*x*y -8.415*y2-0.03693*x3+0.08438*x2*y +0.3875*x

12、*y2+ 1.697*y3+ 0.003434*x4- 0.004379*x3*y- 0.01212*x2*y2- 0.04361*x*y3-0.1633*y4 -0.0001188*x5 +0.0002014*x4*y -0.0001728*x3*y2+ 0.001288*x2*y3+ 0.00111*x*y4 + 0.006215*y5。而且在给出的参数相关度分析中，R-square: 0.9446，R的平方接近于一，表明观测值f（z1）与拟合值较近，表明拟合效果很好。 同样的方法，将另外两个面输入工具箱，可得其曲面的数据如下：拟合的沙层顶部的三维曲线如图2所示。图2 沙层顶部的三维曲线沙

13、层顶部曲线方程为：f（z2） = 4.203 + 3.395*x +15.05*y-0.1935*x2-2.702*x*y-5.333*y2 -0.02551*x3+0.2749*x2*y+0.4546*x*y2+1.024*y3+0.002965*x4-0.008088*x3*y-0.04028*x2*y2 -0.01842*x*y3 -0.106*y4 -9.831e-05*x5 + 0.0001355*x4*y + 0.0005345*x3*y2 + 0.001587*x2*y3 -0.0002715*x*y4 +0.00441*y5。R-square: 0.9165。R的平方接近于一，

14、表明观测值f（z1）与拟合值较近，表明拟合效果很好。拟合的沙层底部的三维曲线如图3所示。图3 沙层底部的三维曲线沙层低部曲线方程为：f（z3） = 21.21 -12.53*x + 17.77*y +2.245*x2+1.189*x*y-13.48*y2 -0.2575*x3 -0.004664*x2*y -0.006048*x*y2 + 3.484*y3+0.017*x4-0.01379*x3*y+0.01125 *x2*y2-0.0325*x*y3-0.3879*y4-0.0004371*x5+0.0004018*x4*y+0.0008946*x3*y2-0.001641*x2*y3 +

15、0.00344*x*y4 + 0.01559*y5。 0.9301。对于这三个曲面方程的拟合，其中部分参数的置信区间包含原点，根据统计回归的知识，这一些项的参数可能会不精确，或者这一些项应该不存在。但通过我们作图发现原方程可以非常好的表示取样散点的分布，而且减少部分参数的方程模拟效果并不理想，考虑到相关指数都逼近一，所以我们决定保留工具箱做出的方程原型。所得地表的等高线图如图4所示。图4 地表的等高线图所得沙层顶部的等高线图如图5所示。图5 沙层顶部的等高线图所得沙层底部的等高线图如图6所示：图6 沙层底部的等高线图5.2 问题二的模型建立与问题求解5.2.1 问题二模型的建立题中所给开采地的

16、地图为网格图，图中曲线围起来的面积为可取样到的有表土和沙子的面积。本模型近似认为在一个网格中，若可取样到的面积在该网格的面积占该网格1/2及以上，则认为表土和沙子占有全部该网格。若可取样到的面积在该网格的面积占该网格1/2以下，则认为该网格不存在表土和沙子。于是不存在表土和沙子的网格为：第一行的第一个、第二个、第七个、第八个，第二行的第一个，第三行的第一个，第四行的第一个，第六行的最后一个，第七行的最后一个，第八行的第五个、第六个、第七个、第十个、第十一个、第十二个、第十三个、第十四个、第十五个。其余网格均存在表土和沙子。图7 近似前的区域图 图8 近似后的区域图用XOY表示题中所给开采地的地

17、图中可取样到的有沙子和表土的平面部分。地下存在三层曲面，即地表、沙层顶部、沙层顶部。设地表曲面方程为，沙层顶部曲面方程为，沙层低部曲面方程为。考虑用二重积分求解体积问题，则表示地表到到平面（x，y）的体积部分，表示沙层顶部到到平面（x，y）的体积部分，表示沙层低部到平面（x，y）的体积部分，利用公式求出表土体积，利用公式求出沙子体积。根据近似之后的地区图形，我们现将全部区域进行二重积分，再减去划出的部分区域的积分。利用Matlab中int（f（z）,y1,y2）,x1,x2）命令计算积分，由此得到的值可以近似作为沙矿区域沙子层的体积和表土层的体积。5.2.2 问题二模型的求解在第一问的求解中我

18、们已经得到了三个曲面的方程，对于求各曲面到XOY平面形成的体积，利用matlab中int函数求二重积分得到。对地表面的积分结果： z1_volumn= 2.2235e+03。对沙子顶面的积分结果z2_volumn = 1.9058e+03。对沙子底面的积分结果z3_volumn = 219.0490。由于题中所给开采地的地图中一个网格的边长代表50m，所以最后经积分所求的体积结果应乘以2500。沙子体积4.216847072638891e+06。表土体积7.943015323319440e+05。经计算表土体积与沙子体积的比值为18.8%，小于19%，满足B.T.I公司的要求。已知表土的平均密

19、度：1350千克/立方米，沙子的平均密度：1620千克/立方米，利用公式算出：沙子质量：6.831292257675003e+06，即683万吨，大于400吨，满足B.T.I公司的要求。B.T.I公司要求的含沙量超过400万吨，表土体积小于沙子体积19%，经计算均满足，所以判断可以对这项工程投资。5.3问题三的模型建立与问题求解5.3.1 开采每吨沙土所获利润函数的构建我们选择只考察沙层厚度大于0. 1米的那些区域作为开采的限制范围. 由合理的假设,该区域没有特别复杂的地质构造，不需牵扯地质知识.可以认为经过了插值之后的网格点处的函数值就是以此点为中心的单位方块内沙层厚度的平均值，这是一种粗糙

20、的估计，但用于本问题易于计算. 计算出沙子的体积是4216800（ m3 ） ,重量是6830000 （ T ） ; 沙子上层的表面土重量是1070000 （ T ） . 由月开挖量10万吨的限制,所以完成该地区沙子的全部开采需68个月工期.但是显然整片土地并不值得完全开采. 原因在于部分区域开采卖出的收益小于开采所需的各项支出,尤其是沙矿边缘的区域根本不值得开采. 究竟哪些区域值得开采,开采施工的顺序如何设计？ 这需要我们建立一个利润函数来表示开采每吨沙土所获的利润大小。假设构建的利润函数h （x , y ） ,在该函数的指导下安排施工区域和顺序。由于 （i）开采一吨沙子付给地区使用费0.5

21、472镑。 （ii）开采一吨沙子并将其运到市场需0.4592镑。 （iii）沙子的销售单价：2.96镑/吨 所以其开采每吨沙子可获利1.9536镑。 （i） 倾倒表土要占用土地，每倾倒一吨表土付0.121镑。 （ii） 将表土挖出并运到附近每吨需0.3672镑。所以每挖掉一吨表土要支出0.4882镑。开采每吨沙土所获利润函数,即单位面积利润与单位面积沙土重量之比函数因为考虑到月开采量有10万吨的限制,我们只选出了该函数大于等于1的范围,我们的施工设计局限于该范围.可以看到上述三个函数的定义式里,只有D1、D2两个变量,其分别是沙层和表土层的厚度,其他都是已知参数. 而D1、D2的厚度可在前两问

22、模型中的曲面方程做差得到。5.3.2 安排施工顺序和区域首先需要确定施工顺序,可以证明,在开采总量确定,每月开采量固定10万吨的前提下, 考虑到贴现的因素,要得到最大的效益,每月的利润应追求在当前条件下尽可能达到最大,即施工按开采每吨沙土所获利润的大小顺序施工,也就是使得上面构建的利润函数h （ x , y ） 递减。对确定的几个划分区域,我们应该遵循利润从多到少的顺序来一块块进行施工.进一步地,我们可以推导出的开采沙子施工原则.开采沙子施工原则: 对所有未施工过的沙矿区域,应按开采每吨沙土所获利润由大到小来安排施工进程, 即使得上面构建利润函数h （x , y ） 递减。利用matlab编程

23、可求得,该函数在坐标（ 11, 3） 处取最大值38.535,意为此处附近动工开挖一吨沙土获利最大,为38.535元.初始工程设计原则上以此点为中心,在其周围按等值线层层展开进行。每次从开挖区域中选出利润最大的一个坐标点，以它为中心挖去十万吨沙子。之后将该区域抛出再进行下一次选择。但是这会带来开初头几个月内施工区域的支离破碎,我们在设计工程时应照顾施工方便而适当予以变动.，根据具体工作情况随机变动。由于月贴现率必须先经过贴现，求出其纯现值。因此，假设从工程开始的第n个月有p镑利润，其纯现值为：.其中i是月贴现率，利用等比级数公式，经过贴现之后的月利润和是：当p=187783.34，i=0.01

24、时，可以求得贴现后的利润总和9874746.5， 但是从中我们还应该减去农作物与房屋赔偿费20万，由于地区建设美化费知道工程结束才付，经贴现应该为镑。所以最终得到的总利润是9874746.5-200000-469447=9205300, 由于计算出的总营业额是14611587.3.可得其利润是63.3%， 所以这项工程是值得投资的。5.3.3 面向实际工作跟细致的分析本题采用的模型是建立在散点上找出的最大值，虽然可以近似的求解，但不够精确，针对本问题更好的方法是使用前两问得出的曲面方程，列出二重积分的方程，利用matlab优化工具箱求解出最大利益的积分上下限，多次重复直到收益不足支出时停止，因

25、为曲面方程是连续的函数，取积分制比散点图的累加值要精确很多。 但因为这种方法需要使用多元积分函数的优化，使得问题求解过于复杂，本文依然选用散点进行利润的计算。6 模型评价及改进6.1 模型的评价优点：该模型利用拟合的方法运用matlab中SFtool拟合工具箱求得地表、沙层顶部、沙层低部的曲面方程，并运用contour语句绘制出等高线，并利用int函数求解出二重积分进而求得沙子和表土的体积。该模型对实际问题有严格的理论推导，对数学软件matlab有很强的应用,很好的体现了数学理论和计算机编程在解决该问题中所起的作用。该模型通过高次幂函数，精确拟合出的曲线方程，并重新规划的开采地地图中有沙土的面

26、积，进而积分求解出开采地下各层体积，很好的体现了求解的严谨性。缺点：作为露天开采，本问题没有考虑复杂的力学结构，没有将地质学相关知识应用到本模型中。推广：该问题的处理方法可应用到其他资源的开采中，如煤炭等。6.2 模型的改进在拟合曲面方程时，某些参数的置信区间包括原点，但并没有将其舍弃，因为经拟合检验保留这些参数能更好地反映散点的趋势。但曲线方程的形式仍需要进一步改进。该模型默认以幂函数的函数形式拟合曲面方程，但还需验证其他形式的函数是否能够产生更好的拟合效果。参考文献1武亦文,张子方,赵勇. 沙子的开采J. 数学的实践与认识,2007,11:99-105.2胡庆婉. 使用MATLAB曲线拟合工具