钻井布局设计(数学建模论文)

1、课 程 设计（论文）任务书专业班级：  应用数学0801 学生姓名：    指导教师（签名）：      一、课程设计（论文）题目钻井布局模型二、本次课程设计（论文）应达到的目的课程设计（论文）是实践教学环节的重要组成部分，其目的是通过课程设计（论文）加深学生对本课程基本知识的理解，提高综合运用知识的能力；掌握本课程的主要内容、工程设计或撰写论文的步骤和方法；提高制图能力，学会应用有关设计资料进行设计计算和理论分析的方法，以提高学生独立分析问题、解决问题的能力，逐步增强实际工程训练。三、本次课程设计（论文）任务的主要内容和要

2、求（包括原始数据、技术参数、设计要求等）              勘探部门在某地区找矿。初步勘探时期已零散地在若干位置上钻井，取得了地质资料。进入系统勘探时期后，要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井的地质资料，不必打这口新井。因此，应该尽量利用旧井，少打新井，以节约钻探费用。 四、应收集的资料及主要参考文献：    

3、60;       收集资料：网格遍历求解、穷举算法、MATLAB编程 等参考文献：1周承高, 廖园.优化方法及应用程序设计M. 北京: 中国铁道出版社, 1989.2 张圣勤编著.MATLAB7.0实用教程M.北京：机械工业出版社，20093刘来福, 曾文艺.数学模型与数学建模M. 北京:  北京师范大学出版社, 1997, 8.4 谢兆鸿编著.数学建模技术M.北京：中国水利水电出版社，2003五、审核批准意见教研室主任（签字）        

4、    目录摘要 31 问题重述 42 模型建设 53 相关变量及其说明 54 模型的建立 54.1 问题1） 54.2 问题2） 64.3 钻井布局算法实现 74.4 对问题3）的讨论 105 参考文献 126 致谢 137 附录 14钻井布局模型摘要勘探部门在某地区找矿。要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井，不必打新井。因此应该尽量利用旧井，少打新井，以节约钻探费用。本文根据钻井布局的实际情况，对钻进布局问题做了深入的研究。力求找到变化中的不变量。

5、讨论了在一个区域内按纵横等距的点来布置井位的问题，通过移动网格使旧井利用最多，采用网格遍历求解的方法得到很好的计算精度。最后得到问题1)的解为：I=2,4,5,10时的旧井可以利用；经分析可知问题1）是问题2）的特殊形式，故问题2)可采用类似的方法，求得的解为：I=1, 6, 7, 8, 9, 11的旧井可以利用。关键词：网格遍历法，结点，移动，覆盖1 问题重述勘探部门在某地区找矿。要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井，不必打新井。因此应该尽量利用旧井，少打新井，以节约钻探费用。设

6、平面上有n个点pi，其坐标为(ai,bi)，I=1，2，n，表示以有的n个井位。新布置的井位是一个正方形网格N的所有结点（所谓“正方形网格”是指每个格子都是正方形的网格；结点是指纵线和横线的交叉点）。假定每个格子的边长都是1单位。整个网格是可以在平面上任意移动的。若一个已知点pi与某个网格结点Xi的距离不超过给定误差=(0.05单位)，则认为pi处的旧井资料可用，不必在结点Xi处打新井。为进行辅助决策，勘探部门要求我们研究以下问题：1） 假定网格的横向和纵向是固定的（比如东西南北），并规定两点间的距离为其横向距离（横坐标之差绝对值）及纵向距离（纵坐标之差绝对值）的最大值。在平面上平行移动网格N

7、，使可利用的旧井数尽可能大。试提供数值计算方法，并对下面的数值例子用计算机进行计算。2） 在欧氏距离的误差意义下，考虑网格的横向和纵向不固定（可以旋转）的情形，给出算法及计算结果。3） 如果有n口旧井，给出判定这些井均可利用的条件和算法（你可以任意选定一种距离）。数值例子 n=12个点的坐标如下表所示：i123456789101112ai0.501.413.003.373.404.724.725.437.578.388.989.50bi2.003.501.503.515.502.006.244.102.014.503.410.802 模型假设1） 正方形网格N可以覆盖n个井位；2) 每个格子的

8、边长为1；3) 旧井位点pi(ai,bi)初始坐标位于第一象限，网格N可以任意移动；4) 把旧井位点看作质点；初始位置网格的某条横线和纵线与坐标轴重合，并记交点Q与原点O重合3 相关变量及其说明1) pi(ai,bi) I=1,2,n：旧井位点的坐标；2) Dpixi：pi到结点xi的距离；3) ：误差；4) (x1i,y1i),(x2i,y2i)：pi所在方格对角线顶点坐标；（x1i<x2i,y1i<y2i）;5) B：网格每次移动的步长 C：每次旋转的角度；6) m,n：横向与纵向移动的次数（横向向右，纵向向上）；7) s：旋转的次数（逆时针）；8) Q(p,q)：变化后Q点的

9、坐标；9) 如果第pi个旧井可利用；则Ti=1否则Ti=0；10) Nj=：第j个位置可利用旧井的总数；如果平移，则j=1,2,mxn;如果旋转，则j=1,2,m x n x s;4 模型的建立4.1 问题1）1 计算移动后pi所在方格对角线顶点坐标,移动后，原对角线顶点 坐标变为（x1i + m.B , y1i + n.B），（x2i + m.B , y2i + n.B）（其中:m.B,n.B<=1）； 如果（x1i+mB>=ai&&y1i+nB>=bi），则(x1i +mB-1,y1i+nB)è(x2i,y2i).即将新坐标作为pi所在方格的(x

10、2i,y2i)同理：如果（x1i+mB>ai&&y1i+nB<bi）则(x1i +mB-1,y1i+nB）è(x1i,y1i).如果（x1i+mB<ai&&y1i+nB>bi）则(x1i +mB,y1i+nB-1）è(x1i,y1i).如果（x1i+mB <ai&&y1i+nB<bi）则（x1i +mB,y1i+nB）è(x1i,y1i).2. 求Dpixi.令d1=|x1i-ai|,d2=|ai-x2i|,d3=|y1i-bi|,d4=|bi-y2i|.则Dpixi=mc(d1,

11、d2,d3,d4) mc(di)指di中次小的.3. 判断误差并确定pi是否可用 如果Dpixi<=,Ti=1如果Dpixi>,Ti=04. 计算NjNj= j=1,2,.,m×n5. N=max(Nj),同时记下m,n值及Q(mB,nB).根据上述步骤，利用计算机计算（见附录）可得最后结果为：可利用点为第2, 4, 5,10 点4.2 问题2）同一坐标系下，将网格N旋转相当于将pi反方向旋转每次旋转的角度c=/r,r=max()所以 c 并且 s*c<=/21. 执行问题1）中1后，以逆时针旋转，旋转后点pi的坐标变为（ai',bi'）其中ai&#

12、39;=*cos(arctan(bi/ai)+s*c),bi'=*sin(arctan(bi/ai)+s*c)2. 求旋转后pi所在方格的对角线顶点坐标如果（ai'>=x1i&&bi'<=y2i）,则（x1i,y1i）è(x1i,y1i)如果（ai'<x1i&&bi'<y2i）,则（x1i-1,y1i）è(x1i,y1i)如果（ai'<x1i&&bi'>y2i）,则（x1i-1,y1i+1）è(x1i,y1i)3. 求 Dpi

13、xid1'=,d2'=d3'=,d4'=则Dpixi=min(d1',d2',d3',d4')4. 执行问题1）的3，4，5根据以上步骤，过程中取B=c=0.05，用计算机编程求得可利用点为6个，即第1, 6, 7, 8, 9, 11 点.4.3 钻井布局算法实现：我们研究的问题是：已知旧井的位置，确定可以重用旧井的最多的个数。假设有一网格可以在井分布的区域移动，这样只要旧井与网格的某个节点距离<=0.05那么认为这口旧井是可以重用的，具体计算机实现时我们考虑到网格的每个小方格是全等的，这样我们在移动网格时只需整体上下移动1

14、单位即可，另外我们是以一个参考点Q(开始在原点)在一单位的小方格中移动来处理的。移动步长我们设定为0.01，即共要移动100*100次。在第二问，考虑旋转的情况下，便于计算，我们移动的是原来的旧井而不是方格，这样便于处理，同时在第一问的基础上每移动一次，我们都要对原来的旧井点坐标逆时针旋转90度（我们步长暂设定为1度，根据计算机实现情况，当步长为1时运行时间已达到15分钟左右），计算可以重用的旧井的个数。到了第三问，我们还是考虑到网格的全等对称性，考虑将所有的旧井移动到近量小的范围，我们用逆向思维考虑：将网格放在适当的位置，如果旧井都可以重用的话，一定可以将这些井的横纵坐标加上m,n单位（n,

15、m可以小于零）将所有的旧井移动到某个节点的周围，并且新位置井位与这个节点的距离<=0.05。这样只要将原来的节点移动到点（1，1）周围，做一个尽量小且能包围所有新位置井位的圆，如果圆的直径<=2*0.05那么这些旧井都可以重用。计算机实现流程图如下：井点坐标数据值显示井点坐标旧井点数数据值输入程序入口处理结果数据处理选择是否需要考虑旋转可进行相应处理纪录未移动前pi所在方格的定点坐标纪录可重用井点变量加1,JNk1:=JNk1+1用MD与0.05比较，若MD<0.05认为此旧井可以重用函数MinD具体实现过程count具体实现求出pi到四个顶点的最短距离MD求出网格移动后pi

16、I所在方格的定点坐标x1,y1,x2,y2纪录s,t的值s,t循环10000 返回循环s,t继续移动网格的位置结果输出：没移动网格一次时，可重用旧井的个数及相应网格的位置（Q点的坐标(s*0.01,t*0.01)）以上是在没有考虑旋转时的流程图，对于考虑旋转时的情况，只需改变一下count过程和MinD函数即可，还需要增加一层循用来考虑旋转90度的情况，解决解决方案与不考虑旋转时大同小异。 用程序解决问题一，我们将运行结果存到文件jieguo.xls中经分析知在我们的精确度范围内最多可以重用井的个数是4个，相应网格的位置是：祥见jieguo.xls。在算法优化的基础上我们还可以对这两个点的周围

17、区域在进行循环查找看看是否有更多的旧井可以重用。 最后我们分析一下我们的算法复杂度，这对于优化算法，提高程序的执行效率是非常重要的。解决问题一时总的运行次数大约是n*10000次，若考虑旋转时运行次数更是庞大的惊人：n*900000，因此程序执行时间是比较长的，大约需要8s左右，分析输出的数据知，很多结果是没用的或者不需要程序循环，因此，我们可以这样优化程序：开始的时候我们的循环步值适当大些，这样可以粗略的找到合适的网格位置，我们再在这个位置上进行细化最终得到最终最优结果，这样处理我们的程序循环次数会减少到至少0.1倍，同样在考虑旋转的情况下 我们也作这样的处理以便优化程序、提高运行效率、减少

18、运行时间。4.4 对问题 3）的讨论1 根据题意，n口旧井均可利用的条件是由n口旧井的位置决定的，也就是说，n口旧井在满足一定位置关系的条件下都是可以重用的。我们要计算的就是这种关系。2. 如果旧井pi(ai,bi)可利用，则pi'(ai+m,bi+n)也可利用。由此，可把所有的旧井位点移到某个点的周围。令d=max(pi pj)（问题1）中的距离）。 pi pj:表示移动后旧井位点任意两点的距离，如果d<=2 则满足条件，即n口旧井均可利用。3. 对具体坐标取值得讨论。我们考虑将所有的旧井点移动到点（1，1）周围。现假设旧井点坐标的小数部分是(xi,yi)，现在分4种情况对旧井

19、点坐标移动到点（1，1）的周围后新坐标的值进行讨论。(1)、xi>=0.5并且yi>=0.5 那么新坐标是(xi,yi)è(xi,yi)；(2)、xi<0.5并且yi>=0.5 那么新坐标是(xi+1,yi)è(xi,yi)；(3)、xi>=0.5并且yi<0.5 那么新坐标是(xi,yi+1)è(xi,yi)；(4)、xi<0.5并且yi<0.5 那么新坐标是(xi+1,yi+1)è(xi,yi)；4. 接下来要求pipj的最大距离d（问题1）中的距离）。令U1=min(xi),U2=max(xi);V1

20、=min(yi),V2=max(yi);（I=1,2,n）则d=max(U2-U1,V2-V1)由此可得所有旧井均可利用的条件：d<=2。5 参考文献1周承高, 廖园.优化方法及应用程序设计M. 北京: 中国铁道出版社, 1989.2 3刘来福, 曾文艺.数学模型与数学建模M. 北京:  北京师范大学出版社, 1997, 8.4 谢兆鸿编著.数学建模技术M.北京：中国水利水电出版社，20036 致谢本篇论文的创作，离不开老师耐心而正确的指导，离不开各类参考资料给予的启发和帮助。在此，我表示郑重的感谢。同时对在本篇论文创作过程中给予过我无私帮助的学长和各位同学表示诚挚的感谢。另外

21、，由于时间紧迫和水平有限，本篇论文难免有所错误和疏漏，请您批评指正，不吝赐教。7 附录：7.1 问题1)的源程序：（1）、自定义过程及函数 function MinD(k:integer): real;var secondD:array1.4 of real; i,j:integer; L:real;/做交换中间变量begin/求的Pi点与所在方框四个点即距离 SecondD1:=abs(x1k-Pk1);SecondD2:=abs(y1k-Pk2); SecondD3:=abs(x2k-pk1);SecondD4:=abs(y2k-Pk2); for i:=1 to 3 do for j:=

22、i+1 to 4 do begin if SecondDi<SecondDj then begin L:=SecondDi; SecondDi:=SecondDj; SecondDj:=L; end; end; result:=SecondD3;end;function MinDK(k:integer): real;var secondD:array1.4 of real;/存储四个点的即距离 i,j:integer; L:real;begin /求的Pi点与所在方框四个点即距离SecondD1:=sqrt(Pk1-x1k)*(Pk1-x1k)+(Pk2-y1k)*(Pk2-y1k);S

23、econdD2:=sqrt(Pk1-x1k)*(Pk1-x1k)+(Pk2-y2k)*(Pk2-y2k);SecondD3:=sqrt(Pk1-x2k)*(Pk1-x2k)+(Pk2-y1k)*(Pk2-y1k);SecondD4:=sqrt(Pk1-x2k)*(Pk1-x2k)+(Pk2-y2k)*(Pk2-y2k);/对距离进行排序，同时返回最小值 for i:=1 to 3 do for j:= 1 to 4 do begin if SecondDi<SecondDj then begin l:=SecondDi; SecondDi:=SecondDj; SecondDj:=L;

24、end; end; result:=SecondD3;end;procedure TForm1.count(a:real;b:real;x:real;y:real;k:integer);varboo1,boo2:boolean;begin if x>a then boo1:=true else boo1:=false; if y>b then Boo2:=true else boo2:=false; / 求出坐标 if boo1 and boo2 then begin x2k:=x1k; y2k:=y1k; x1k:=x1k-1;y1k:=y1k-1; end; if boo1 a

25、nd not boo2 then begin x2k:=x1k-1; y2k:=y1k+1; end; if not boo1 and boo2 then begin x2k:=x1k+1; y2k:=y1k-1; end;end;procedure TForm1.countK(a:real;b:real;x:real;y:real;k:integer);varboo1,boo2:boolean;begin boo1:=a<x;boo2:=b>y; / 求出坐标 if boo1 and boo2 then begin end; if not boo1 and boo2 then b

26、egin x2k:=x1k+1; y2k:=y1k+1; end; if not boo1 and not boo2 then begin x2k:=x1k+1; y2k:=y1k-1; end; if boo1 and not boo2 then begin x2k:=x1k-1; y2k:=y1k-1; end;end;（2）、变量定义部分 var Form1: TForm1; P,P1:array of ar;/存储旧井的坐标 JN:array1.10000 of ar1; /存储结果 x1,y1,x2,y2:array of real; /存储某个P点所在的方格对角线的俩顶点的坐标 J

27、W:array1.10000 of string; /不考虑旋转的情况结果 i:integer; n:integer; /旧井的个数 s,t:integer;/s,t网格移动部数（3）、程序实现部分 var k:integer; MD:real;/存储最短距离 mess:string;begin for s:=0 to 99 do for t:=0 to 99 do begin for k:=1 to n do begin x1k:=int(Pk1); y1k:=int(Pk2); x2k:=x1k+1+0.01*s; y2k:=y1k+1+0.01*t; x1k:=x1k+0.01*s; y

28、1k:=y1k+0.01*t; count(Pk1,Pk2,x1k,y1k,k);/求得i所在的方格的对郊县的两顶点的坐标 MD:=MinD(k); if MD<0.05 then begin JNs*100+t+11:=JNs*100+t+11+1;/如果最小距离小与0.05 那么可用井的个数JNs*t+11加1； JWs*100+t+1:=JWs*100+t+1+' '+inttostr(k); end; end; JNs*100+t+12:=s; JNs*100+t+13:=t; mess:=floattostr(JNs*100+t+11)+'个'+char(9)+floattostr(JNs*100+t+12*0.01)+char(9)+floattostr(JNs*100+t+13*0.01)+char(9)+JWs*100+t+1; end;7.2 问题2)的源程序： var k,o:integer;/k,循环n次计算n各界点；o 循环90次 MDk:real;/存储最短距离 mess:string; f_x,f_y,bi:real;beginfor k:=1 to n do for o:=1 to 2 do P1ko:=Pko; /初始化P1k 用于初始化 for s:=0 to 99 do