《矿井布局问题论文_B09050228(徐力).docx》由会员分享，可在线阅读，更多相关《矿井布局问题论文_B09050228(徐力).docx（14页珍藏版）》请在三一办公上搜索。
1、钻井布局问题的数学模型摘要勘探部门在某地区找矿时，首先进行初步勘探，取几个位置钻井，取得地质资料；然后进行系统勘探，进行纵横等距的撒网式钻井。显然如果能尽可能多的在系统勘探时利用初步勘探的钻井资料，就能有效的节约费用。 在不考虑网格方向的情况下，本文首先给出了两个结论，即网格的位置由节点唯一确定，与原始矿井节点以及单位方格内的矿井映射点有相同的性质。这样就将问题等效为在单位方格内确定网格的一个节点。要解决这个问题，首先我们提出运用一般的搜索法对网格节点在单位方格内进行遍历（模型一）。通过对遍历算法进行有效的优化，大量减少了搜索的次数，进而初步计算得到了原井位最多有4个可被利用，并给出了方格节点
2、的坐标为：考虑到搜索算法的复杂度，我们给出了模型二，即在单位方格内通过确定每个矿井映射节点被利用时节点的区域，来找出方格内被这些区域覆盖次数最高的部分，显然如果将节点放在这部分内，将会有最多的点被利用，从而也就确定了节点的位置范围。运用MATLAB进行计算与判别，得到最多有4个可被利用，并求出了网格节点坐标具体的范围：当网格方向可以改变时，我们建立了模型三。考虑到判别条件是欧氏距离，可以将原题简化为一个圆形进行覆盖，圆的半径为，再用类比利用模型二进行判断，那么就能相应的找到最优规划。模型三首先进行了误差分析，根据假设的误差使用夹逼法则，然后，为了减小搜索范围，我们证明了Xi=ai+-s,Yi=
3、bi+-t时，最多有6个矿井可被利用。对于第三问的判定算法，我们仍然根据模型三，建立假设模型四。构造出两个极端情况，此时所有矿井均可被利用。具体算法的见问题三分析步骤。 最后我们对模型四的一个假设进行了检验。虽然这个假设严格的说并不成立，但通过我们用蒙特卡罗方法进行多次模拟，发现假设成立的概率极高。综上，我们可以先用模型四进行计算，再对结果进行检验，对极少数不成立的，可以综合特殊情况进行考虑。关键词 勘探 矿井 遍历算法 蒙特卡罗 数值分析 误差分析 假设检验一 问题重述勘探部门在某地区找矿。初步勘探时期已零散地在若干位置上钻井，取得了地质资料。进入系统勘探时期后，要在一个区域内按纵横等距的网
4、格点来布置井位，进行“撒网式”全面钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井的地质资料，不必打这口新井。因此，应该尽量利用旧井，少打新井，以节约钻探费用。比如，钻一口新井的费用为500万元，利用旧井资料的费用为10万元，则利用一口旧井就节约费用490万元。 设平面上有n个点Pi，其坐标为(ai, bi),i=1,2,n，表示已有的n个井位。新布置的井位是一个正方形网格N的所有结点（所谓“正方形网格”是指每个格子都是正方形的网格；结点是指纵线和横线的交叉点）。假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。整个网格是可以在平面上任意移
5、动的。若一个已知点Pi与某个网格结点Xi的距离不超过给定误差（0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。为进行辅助决策，勘探部门要求我们研究如下问题： 1)假定网格的横向和纵向是固定的（比如东西向和南北向），并规定两点间的距离为其横向距离（横坐标之差绝对值）及纵向距离（纵坐标之差绝对值）的最大值。在平面上平行移动网格N，使可利用的旧井数尽可能大。试提供数值计算方法，并对下面的数值例子用计算机进行计算。 2)在欧氏距离的误差意义下，考虑网格的横向和纵向不固定（可以旋转）的情形，给出算法及计算结果。 3)如果有n口旧井，给出判定这些井均可利用的条件和算法（你可以任意选
6、定一种距离）。 n=12个点的坐标如下表所示： i123456789101112ai0.50 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50 bi2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80 二 符号说明Pi(ai,bi) 原始勘测中任意一点Pi(ai，,bi) 将旧井移至单位正方形里面时候的新井Pif可利用的最大钻井数n 正方形的候选位置的数目 旧井离新表格节点的距离Q（s,t） 网格的一个节点ui 判断旧井节点是否在表格范围内的标志位Xi ，Yi a,b等
7、原始节点需要转化至单位正方形内部时候需要移动的整数大小。 三 模型假设1. 初步勘探时，所取的点比较分散，不存在两个点的资料被系统勘探时同一个点利用的情况。2. 系统勘探所取的点数比初步勘探时取的点多。3. 初步勘探所取的点在系统勘探的勘探范围之内。4. 假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。整个网格是可以在平面上任意移动的。5. 若一个已知点Pi与某个网格结点Xi的距离不超过给定误差（=0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。四 问题的分析及说明考虑到勘探部门在某地区找矿时是分为两步进行的，即：第一步，初步勘探，取得地质资料。第二步，
8、系统勘探，进行“撒网式”全面钻探。那么，为了节约钻探费用，我们自然希望充分利用第一步的数据，来减少勘探次数。对实际问题进行分析，可以认为初步勘探时钻井的位置可对应与二维坐标系中的点，成为初始点。系统勘探时的全面钻探可认为是二维坐标系中横坐标与纵坐标间距相等的点为节点所构成的网格。这样，问题的实质就转化为如何定位网格，使尽可能多的初始点位于网格节点的误差范围内。为了更好的表述问题，下面给出两条结论：结论一：假定网格的横向和纵向是固定的，那么只要确定其中一个节点，就可确定整个网格。证明：很显然，当网格节点中的一个点的坐标确定时，其他点的坐标可由如下公式确定：结论二：对于坐标系里的每一个点，可以定义
9、其映射点，映射点与原来的点对网格有相同的位置关系。证明：意为对向下取整。由于网格横向和纵向的单位都是1，所以当点位于网格中一个节点的误差范围内时，则其映射点必然位于的误差范围内。反之如果点不在任何一个节点的误差范围内，那么映射点也不在任何一个节点的误差范围内。五 模型的建立模型一 遍历法考虑已经给定点（矿井）以及网格的横向与纵向固定，那么问题就变为如何确定出网格的最优平移位置。由结论一，定出网格的位置只要确定出一个节点，同时，由于网格以1为单位，那么在单位方格内，必然有且只有一个节点，这样我们在四个点所构成的方块（单位方格）内对节点进行搜索，就能将网格的全部可能的平移情况进行遍历。考虑到题目中
10、所给数据精确到0.01，那么可取0.01为两个坐标的步长进行搜索。对于搜索到的每一个坐标，作为网格的节点，来确定整个网格，再计算到底有多少初始点能进入网格节点的误差范围，最后进行比较，选出容纳初始点最多的网格作为解决方案。搜索法从理论上是可行的，然而这种方法的计算量往往比较大。如上例，要搜索的节点个数为（即循环次数）：10010010000（次）而每一次循环都要对12个点依次进行判断，这必然十分复杂，当实际的数据量比较大，数据精度要求比较高时，这种直接搜索的方法计算起来就十分耗时，甚至是不可行的。为了减少循环次数，必须想办法对算法进行优化。这里我们考虑对于每一个初始点，由结论二，首先在单位方格
11、内找到该点的映射点，这样就将所有点都映射到了单位方格内。然后进行搜索，在映射点附近搜索网格节点，使映射点位于网格节点的误差范围内，再对节点构成的网格进行分析、比较，选出其中的最优。经过这样的过程，循环的次数为：1111121452（次）这样就大大减少了循环次数，优化了搜索算法，使之耗时短，使方案变得可行。 运用优化后的搜索算法，对所给数据进行计算，得到最多可以有4口井不必打而可以利用初步勘探的资料。可以利用的初始点为： (1.41,3.50) (3.37,3.51) (3.40,5.50) (8.38,4.50)网格坐标为：模型二 框图分析法考虑到遍历节点法毕竟是一种复杂度极高的算法，当数据量
12、大的时候就会有严重缺陷，因此这里给出一种简单易行而又十分精确的算法：方框图分析法。 首先，根据定理二，确定出初始点在单位方格内的坐标（即映射点），变换前后点的坐标如下：表1：坐标映射表n123456789101112初始点ai0.50 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50 bi2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80 变换点ai0.50 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50
13、bi0.00 0.50 0.50 0.51 0.50 0.00 0.24 0.10 0.01 0.50 0.41 0.80 使用matlab软件将变换前的点画在一坐标系下，如图1：点图：1同时将映射之后的点画于一个坐标系下面得到如图2：映射点图：2将所有的旧井节点平移至单位正方形可以得到Q=（x,y）|0=x=1,0=yPiai,bi ai=ai-ai, bi=bi-bi在上述变换后，问题1大致等价于用一个变长为2的正方形去覆盖尽可能多的Pi,正方形的中心就是网格的一个节点所在位置。如图-3所示：观察可得，正方形有n2个候选位置。若正方形左边经过Pi，右边经过Pj，则aj【ai, ai+2*】
14、,bj【bi-2*,bi】代入数据得：aj【ai, ai+0.1】,bj【bi-0.1,bi】（1）根据上表1可知：ai=【0.5,0.41,0.00,0.37,0.40,0.72,0.72,0.43,0.57,0.38,0.98,0.50】 bi=【0.00,0.50,0.50,0.51,0.50,0.00,0.24,0.10,0.01,0.50,0.41,0.80】;将以上数据代入（1）式得：aj0.5,0.6,0.41,0.51,0.00,0.10,0.37,0.47,0.40,0.50,0.72,0.82,0.72,0.82,0.43,0.53,0.57,0.67,0.38,0.48,
15、0.98,1.08,0.50,0.60; bj-0.1,0.00,0.40,0.50,0.40,0.50,0.41,0.51,0.40,0.50,-0.10,0.00,0.14,0.24,0.00,0.10, -0.09,0.01,0.40,0.50,0.31,0.41,0.70,0.80;经过将上述数据使用matlab数学软件分析可以得到，满足条件的旧井数目为4，方格节点的纵横坐标为（0.42,0.46），亦即（0.42+xi,0.46+yi）, xi，yiZ模型三 动态分析法网格平行移动，对网格上面的任一点Pi，当网格移动整数个单位时，Pi相对于最近的网格节点的距离不变， Pi在网格上平移
16、整数个单位时，其相对于最近的网格节点的距离不变。则现在把所有旧井节点平移至正方形：Q=(x,y)|0=x=1,0=y=1;设网格的一个节点为(s,t)，0=s,t设Pi:ui=1否则ui=0；则以上可以归纳为：（2）等价于以某一角度为步长旋转网格，一个半径为的圆覆盖更多的Pi，圆心是网格的一个节点位置。旋转模式：如下图可旋转的网格图4：图4Pi可利用当且仅当（3）最多可以利用旧井位为max1nui,其中，ui=10，。代入数据，即(s+xi-ai)2+t+yi-bin=0.0025xi,yi 化简该算式得：(-s-xi+ai)2+-t-yi+bin a=0.5 1.41 3.00 3.37 3
17、.40 4.72 5.43 7.57 8.38 8.98 9.50; b=2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.80; a=0.5 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50; plot(a,b,r) plot(a,b,*) hold on grid mkdir chapter_1999 cd chapter_1999 a=0.5 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50; b=0.00
18、 0.50 0.50 0.51 0.50 0.00 0.24 0.10 0.01 0.50 0.41 0.80; plot(a,b,*) grid on plot(a,b,*) grid on a=0.5 0.41 0.00 0.37 0.40 0.72 0.72 0.43 0.57 0.38 0.98 0.50; b=0.00 0.50 0.50 0.51 0.50 0.00 0.24 0.10 0.01 0.50 0.41 0.80; c=a=0.37&a d=b=0.41&b avg_x=(0.37+0.47)/2avg_x = 0.4200 avg_y=(0.51+0.41)/2avg
19、_y = 0.4600第二问：C语言代码：#include#include#define kesi 0.05int main()double a12=0.5,1.41,3.00,3.37,3.40,4.72,5.43,7.57,8.38,8.98,9.50;double b12=2.00,3.50,1.50,3.51,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.80;double c12,d12,u_112,u_212,s,t;int i=0,flag=0,k=0,j=0,p=0,x10;while(i12)ci=fabs(ai-floor(ai);di=fabs(bi-floor(bi);i+;for(i=0;i12;i+)printf(n);printf(c%d=%.2lf ,i,ci);printf(n);printf(d%d=%.2lf ,i,di);while(k12)for(j=0;j12;j+)u_1k=pow(ck-cj,2)+pow(dk-dj,2);u_2k=sqrt(u_1k);printf(nu_2%d=%lf ,j=%d n,k,u_2k,j);if(u_2k=0.1&u_2k!=0)flag+=1;/printf(k=%d,j=%dn,k,j);k+;printf(nflag=%dn,flag);return 0;
