钻井布局问题的数学模型讨论.docx

1、钻井布局问题的数学模型讨论钻井布局问题的数学模型讨论 论文摘要： 本文利用数学规划模型，从旧井的可利用条件出发，对钻井的网格移动方案进行修整得出结论，当网格仅可平移时，共有四口旧井可被利用，分别为P2（1.4，3.50)、P4(3.37,3.51)、P5(3.40,5.50)、P10(8.30,4.50)。当网格可以旋转移动时，共有6口井可被利用，分别为：P1(0.5,2.0)、P6（4.72，2.00）、P7（4.72，6.24）、P8（5.43，4.10）、P9（7.57，2.01）、P11（8.98，3.41） 本文在最后对各种建井方案进行了进一步分析与讨论，将模型进一步完整化。 对第三

2、个问题进行了较深入的探讨与研究，建立了三种模型，对三种方案进行比较。 关键字：网格，平移，旋转，可利用旧井数，对称性,欧氏距离1 问题重述：勘探部门在某地区找矿。初步勘探时期已零散地在若干位置上钻井，取得了地质资料。进入系统勘探时期后，要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面钻探。由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井的地质资料，不必打这口新井。因此，应该尽量利用旧井，少打新井，以节约钻探费用。假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。整个网格是可以在平面上任意移动的。若一个已知点Pi与某个网格结点Xi的距

3、离不超过给定误差（=0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。 现提供数值计算方法：1) 假定网格的横向和纵向是固定的（比如东西向和南北向），并规定两点间的距离为其横向距离（横坐标之差绝对值）及纵向距离（纵坐标之差绝对值）的最大值。在平面上平行移动网格N，使可利用的旧井数尽可能大。试提供数值计算方法，并对下面的数值例子用计算机进行计算。 2) 在欧氏距离的误差意义下，考虑网格的横向和纵向不固定（可以旋转）的情形，给出算法及计算结果。 3) 如果有n口旧井，给出判定这些井均可利用的条件和算法（你可以任意选定一种距离）。数值例子n=12个点的坐标如下表所示i 1 2

4、3 4 5 6 7 8 9 10 11 12ai 0.50 1.41 3.00 3.37 3.40 4.72 4.72 5.43 7.57 8.38 8.98 9.50bi 2.00 3.50 1.50 3.51 5.50 2.00 6.24 4.10 2.01 4.50 3.41 0.802 问题的分析对问题一建立模型一，对问题二建立模型二 图1 题目已知中所给各点分布情况2.1 问题一网格的横向与纵向为固定的，所以对网格的移动为线性平移，故可以将网格看作先进行横向平移，在进行纵向平移。效果与直接平移相同，利用将平移分成两步来进行，利用决策变量及约束条件对决策目标进行图 2 问题一平移方案示

5、例2.2 问题二网格的横向和纵向不固定（可以旋转）的情况下进行与问题一相同的求解与讨论。将其分为两步，第一步以原点为轴，以某角度旋转网格，再沿着网格的自然方向进行平移，与第一问中的操作相同，此两步加起来与进行一步旋转平移效果相同。图3-1 旋转后的网格平移方案示例（第一步操作后）图3-2 问题二网格平移方案示例（第二步操作后）2.3 设平行移动的距离为xi,yi(xi为横向位移大小，yi为纵向位移大小)，设横、纵向位移的小数部分分别为xi,yi。由题意知，旧井的可利用与否仅与其距离最近的网点的距离有关，而与该网点是哪一个无关。故横向平移xi 的效果与横向平移xi的效果相同。故在此取横向平移距离

6、xi=0，1）。同理，横向平移距离yi =0，1）。在对第二个问题的解决中，设网格旋转角度为，分析知，当其转过/2角度后，再转，以前面进行的旋转效果相同，故在此假设的范围为0，/2）。2.4 修正函数 设修正函数，由于判断点的坐标与其相邻最近网格点的距离 若x 0.5，x自减1直到其小于0.5 若x-0.5，x自加1直到其大于-0.53 问题假设：每口井在坐标系中可被看为一个点。取每个格子的边长为一个单位。假定网格的横向和纵向是固定的（比如东西向和南北向）。平移网格时，非连续的平移，为以某一小单位渐进平移。最后取得一个确定的平移网格方案。（此处可改进，见模型二）设在任意点钻井均不存在自然障碍4

7、 符号说明：表格1-1 问题1中的符号说明点的初始横坐标点的初始纵坐标平移后的横坐标平移后的纵坐标横坐标平移的距离纵坐标平移的距离平移后可利用的旧井数aibiaicdxdyni=1,2,12表格2-1 问题2中的符号说明点初始横坐标点初始纵坐标旋转角度旋转后横坐标旋转后纵坐标横坐标平移纵坐标平移移动后横坐标移动后纵坐标移后可用旧井数aibiaicdxdya1ib1ini=1,2,12各井的坐标为pi（i=1，2n）移动后各点坐标为pi5 问题一模型的建立：目标函数 max s.t.0dx1;0dy1;问题解决方案： 将移动范围确定后，用穷举尝试法利用计算机编程实现对网格自然方向平移的模拟，存下

8、每种网格移动方案中的旧井利用数， 5.1 模型求解：由约束条件及决策目标建立C+源代码程序求解，得到横坐标平移数：x=0.58纵坐标平移数：y=0.45可利用旧井数：n=4表格1-2 不同步长下利用旧井最多的方案 步长移动距离 0.10.010.0010.0001横向移动距离dx0.60.580.580.580.58纵向移动距离dy0.50.450.450.450.45可用旧井数n444 44由于初始数据给定坐标的精确度为小数点后两位，故取移动的步长为0.01，即可得到与连续变换相同的答案。5.2 讨论：表格1-3 各点移动前后坐标及旧井利用数旧井标号P1P2P3P4P5P6P7P8P9P10

9、P11P12横坐标（移动前）0.501.413.003.373.404.724.725.437.578.388.989.50纵坐标（移动前）2.003.501.503.515.502.006.244.102.014.503.410.80横坐标（移动后）1.081.993.583.953.985.305.306.018.158.969.5610.08纵坐标（移动后）2.453.951.953.965.952.456.694.552.464.953.861.25移动后可用否否可否可可否否否否可否否由上表知，可利用的四口旧井分别为 P2、P4、P5、P10号。5.3 平移后与平移前的坐标比较图4 平

10、移前后的坐标示意图注：112为平移前各点坐标，112为平移后各点坐标。6 问题二模型而的建立：目标函数 N= MAX s.t.0dx1;0dy1;(其中：a1i=cos(tan ()-) +dxb1i=sin(tan ()-) +dy)问题解决方案： 将移动范围与旋转角度范围确定后，利用计算机编程用穷举尝试法先对网格进行旋转，再实现对网格自然方向平移的模拟，记录下每种网格移动方案中的旧井利用数,取其旧井利用数的最大的方案。6.1模型求解：表格2-2 不同步长下利用旧井最多的方案 的步长移动距离 10.10.01旋转的角度00.80.780.78横向移动距离dx0.580.250.210.21纵

11、向移动距离dy0.4500.890.89可用旧井数n4566由约束条件及决策目标编写C+程序求解（程序源代码见附录二），得到旋转角度=0.78 rad横坐标平移数：dx=0.21纵坐标平移数：dy=0.89最大可利用旧井数：n=6由分析知，当旋转角度的步长 为0.01时即可达到最大利用旧井的效果，而当步长进一步缩小时，会加大计算机的预算负荷，使运算时间大大增加，但并不能使方案更优化，所以在实际编程时，可直接取用0.01rad的步长，即可得到预期结果。6.2 移动以后各点情况表2-3 旋转后各点坐标及旧井可利用情况旧井标号P1P2P3P4P5P6P7P8P9P10P11P12横坐标（移动前）0.

12、501.413.003.373.404.724.725.437.578.388.989.50纵坐标（移动前）2.003.501.503.515.502.006.244.102.014.503.410.80横坐标（移动后）1.973.673.405.076.504.977.956.957.019.338.997.52纵坐标（移动后）1.962.39-0.151.012.41-1.012.01-0.01-3.00-1.80-3.00-5.22移动后可用否可否否否否可可可可否可否共有6口旧井可被利用，分别为P1、P6、P7、P8、P9、P11号旧井网格旋转方案为 ： 网格以原点为轴旋转弧度为0.78

13、rad的角度后，再将网格沿着新的自然方向，横移0.21个单位后，再纵向移动0.89个单位。6.3 移动后各点与移动前的位置示意图图5 移动前后坐标示意图注：1*12*为旋转后的点，1-12-为移动后的点。放大的字体为可利用的旧井所在点7 对模型一和模型二的讨论：显然，方案一为方案二的一种特殊情况。在方案一的基础上，添加新的参变量，将模型进一步完善。新变量的引入，使网格变化方案更加的灵活多样，有了更多的方案选择，自然可以找到更适合的方案，使旧井利用率得到提升（如图1） 第一问的旧井利用情况 第二问中的旧井利用情况 图6 两种情况的旧井利用率情况比较（注：图中黑色部分为可利用旧井数占总旧井数的百分

14、比） 旧井利用率得到提升，自然费用就相应的减少，见图2，两种网格移动方案的减少费用比较情况。图7 两种情况下的费用比较注： 1为第一问中的旧井利用方案的费用节省情况，2为第二问中的旧井利用方案的费用节省情况。显然，网格的移动方案的限制越少，旧井的利用率便会在某限定范围内相应的得到提高。8 对模型的修正：8.1 对模型一的修整由于非连续的移动网格，所以，通过计算机模拟计算得到的为旧井利用数最高的方案之一。而实际上的网格移动方案应该是在一个范围内。通过计算机编程，得到：dx = 0.58，0.63 dy = 0.45，0.54网格按此范围内的方案移动均可达到最大旧井利用数 n=4（C+程序源代码见

15、附录4）8.2 对模型二的修整由于非连续的移动网格，所以，通过计算机模拟计算得到的为旧井利用数最高的方案之一。而实际上的网格移动方案应该是在一个范围内。通过计算机编程，得到：dx = 0.21，0.27dy = 0.89，0.97=0.78，0.79网格按此范围内的方案移动均可达到最大旧井利用数 n=6（C+程序源代码见附录5）9 问题三的分析与建模题目要求给出判定的条件和算法，而不同的算法得出的结论必定不同。在此，我们沿用欧式距离的误差意义用于判断旧井是否可用。给出如下三种模型：9.1 模型3-1利用第二个问题的模型二，将已知各N个点的坐标代入到前面的模型中，求其利用旧井最大数，当n=N时，

16、即为当前的N个旧井均可以被利用。显然，此模型略显简单，可进一步改进。9.21 模型3-2：引理：对于任意整数N，存在任意两个正整数N1,N2（零不包括在内），使N= 则N1,N2是该等式的唯一解。取N点坐标为（x,y）由坐标转换公式得以适当的角度旋转坐标轴后的N点的坐标为(xcos-ysin,ycos+xsin)必要条件：任意该两点的距离d满足 （1）当该条件满足时，我们只能判断该两点满足条件，还需判断其它所有点是否满足该条件。取任一点X，以步长0.01从0到 /2找寻，当xcos-ysin与接近的整数相差不超过，此时由于（1）公式ycos+xsin也与接近的整数相差不超过，此时将该代入个点，

17、当出现xcos-ysin与接近的整数相差超过，则继续寻找，若每个都满足ycos+xsin也与接近的整数相差不超过，则这些点都能被利用。9.2.2 模型3-2的讨论： 由于第三问用了筛选，即当不满足条件时就更换，与第二问相比，运算量将大大减少。且该充分条件的判定的加入会筛掉许多数值，比第二个模型更完善。9.3 模型3-3：（该模型将规划的网格图形改变）经过分析，有如下三种对称图形（包括初始方案）9.3.1 方案一：形状规则以概率论方面来说，可利用井的个数n=S*n0，其中S为单个点周围的可行解域（本题以欧氏空间计算），n0为单位面积需打井的个数（为常数），所以改变形状不能优化。9.3.2 方案二

18、：形状不规则我们可以根据实际给点的情况拉长图形的一边，使点位于图形交点附近，为度量图形形变的严重程度，再每单元占单位面积的假设下，我们以边长改变x后该单位变化的面积S作为度量，S值越小代表图形形变程度越小，S值最小的图形为最佳方案。假设一：单位面积打一口井，所以一个单元占一个面积。假设二：该点放入该点属于的单元计算。假设三：该偏离点在该单元的一条边上（为了解题方便，可以这样设，用在边上的点算出的形变度大的话，在中间的点的形变度也一定大）1.正方形网络：当其中一个边发生变化T时 单位边长正方形变化面积为：S=1T=T2.三角形网络：当三角形的一个边发生变化x时： 单位边长三角形面积 S=单位面积

19、三角形边长X= 形变度S=0.33x3.六边形网络：当六边形的一个边发生变化时t时：单位边长六边形面积S=1单位面积六边形边长形变度由上可以得出结论，六边形网络为最佳规划。所以，在此模型中， 沿用欧式距离的误差意义，建井布局网格选用正六边形。参考文献：1. 姜启源 数学模型M 北京: 高等教育出版社, 2003, 82.冯杰 数学模型原理与案例 北京：科学出版社，2007.13.韩中庚 数学建模方法及其应用 北京：高等教育出版社 2005，6附录1： 第一个问题解决时所用的C+程序源代码如下：#include#includevoid change(double &x) for(;x0.5;)

20、x=x-1;void main() double a12=0.5,1.41,3.0,3.37,3.4,4.72,4.72,5.43,7.57,8.38,8.98,9.50; double b12=2.0,3.50,1.5,3.51,5.5,2.00,6.24,4.10,2.01,4.50,3.41,0.80; double x,y,x0,y0; double e,c; int i,n,m; m=0; for(x=0;x1;x+=0.001) for(y=0;y1;y+=0.001) n=0; for(i=0;i12;i+) e=x+ai; c=y+bi; change(e); change(c

21、); if(fabs(e) =0.05&fabs(c) m)x0=x;y0=y;m=n; coutx=x0endl; couty=y0endl; coutn=mendl;附录二：第二个问题求解时所用的C+程序源代码：#include#includevoid change(double &x)for(;x=0.5;) x=x-1; for(;x=-0.5;) x=x+1;void main() double a12=0.50,1.41,3.00,3.37,3.40,4.72,4.72,5.43,7.57,8.38,8.98,9.50; double b12=2.00,3.50,1.50,3.51

22、,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.80; double x,y,x0,y0,sta0; double e,c; double t,s; double sta; int i,n,m; m=0; for(x=0;x1;x+=0.001) for(y=0;y1;y+=0.001) for(sta=0;sta90;sta=0.1) n=0; for(i=0;i12;i+) e=x+ai; c=y+bi; s=sin(atan2(y,x)-sta); t=cos(atan2(y,x)-sta); change(s); change(t); if(sqrt(s*s

23、+t*t)m) x0=x;y0=y;sta0=sta;m=n; coutx=x0endl; couty=y0endl; coutsta=sta0endl; coutn=mendl;附录3 ： 问题二判断旧井利用情况 与 旋转网格后 各井的新坐标 程序#include#includevoid change(double &x) for(;x=0.5;) x=x-1; for(;x=-0.5;) x=x+1;void main() double a12=0.50,1.41,3.00,3.37,3.40,4.72,4.72,5.43,7.57,8.38,8.98,9.50; double b12=2

24、.00,3.50,1.50,3.51,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.80; double p,q; double sta=0.78; for(int i=0;i12;i+) p=sqrt(ai*ai+bi*bi)*cos(atan2(bi,ai)-sta)+0.21; q=sqrt(ai*ai+bi*bi)*sin(atan2(bi,ai)-sta)+0.89; coutaip ; coutbiqendl; change(p); change(q); if(sqrt(p*p+q*q)=0.05) coutendl-the ipoint is avai

25、lable!-endlendlendl; else coutendl; 附录4：修正后模型一的程序源代码#include#includevoid change(double &x) for(;x=0.5;) x=x-1;void main() double a12=0.5,1.41,3.0,3.37,3.4,4.72,4.72,5.43,7.57,8.38,8.98,9.50; double b12=2.0,3.50,1.5,3.51,5.5,2.00,6.24,4.10,2.01,4.50,3.41,0.80; double x,y,x0,y0,x1,y1; double e,c; int

26、i,n,m; m=0; for(x=0;x1;x+=0.01) for(y=0;y1;y+=0.01) n=0; for(i=0;i12;i+) e=x+ai; c=y+bi; change(e); change(c); if(fabs(e)=0.05&fabs(c)m)x0=x;y0=y;m=n;x1=x;y1=y; else if(n=m) if(xx0)x1=x; if(yy0)y1=y; coutx=x0x1endl; couty=y0y1endl; coutm=mendl;附录五 修正后的方案二的C+源程序#include#includevoid change(double &x) for(;x=0.5;) x=x-1; for(;x=-0.5;) x=x+1;void main() double a12=0.50,1.41,3.00,3.37,3.40,4.72,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.8