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

资源描述

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

《钻井布局问题的数学模型讨论.docx》由会员分享，可在线阅读，更多相关《钻井布局问题的数学模型讨论.docx（22页珍藏版）》请在冰豆网上搜索。

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

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

论文摘要：

本文利用数学规划模型，从旧井的可利用条件出发，对钻井的网格移动方案进行修整得出结论，当网格仅可平移时，共有四口旧井可被利用，分别为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）

本文在最后对各种建井方案进行了进一步分析与讨论，将模型进一步完整化。

对第三个问题进行了较深入的探讨与研究，建立了三种模型，对三种方案进行比较。

关键字：

网格，平移，旋转，可利用旧井数，对称性,欧氏距离

1问题重述：

勘探部门在某地区找矿。

初步勘探时期已零散地在若干位置上钻井，取得了地质资料。

进入系统勘探时期后，要在一个区域内按纵横等距的网格点来布置井位，进行“撒网式”全面钻探。

由于钻一口井的费用很高，如果新设计的井位与原有井位重合（或相当接近），便可利用旧井的地质资料，不必打这口新井。

因此，应该尽量利用旧井，少打新井，以节约钻探费用。

假定每个格子的边长（井位的纵横间距）都是1单位（比如100米）。

整个网格是可以在平面上任意移动的。

若一个已知点Pi与某个网格结点Xi的距离不超过给定误差ε（=0.05单位），则认为Pi处的旧井资料可以利用，不必在结点Xi处打新井。

现提供数值计算方法：

1）假定网格的横向和纵向是固定的（比如东西向和南北向），并规定两点间的距离为其横向距离（横坐标之差绝对值）及纵向距离（纵坐标之差绝对值）的最大值。

在平面上平行移动网格N，使可利用的旧井数尽可能大。

试提供数值计算方法，并对下面的数值例子用计算机进行计算。

2）在欧氏距离的误差意义下，考虑网格的横向和纵向不固定（可以旋转）的情形，给出算法及计算结果。

3）如果有n口旧井，给出判定这些井均可利用的条件和算法（你可以任意选定一种距离）。

数值例子n=12个点的坐标如下表所示

i123456789101112

ai0.501.413.003.373.404.724.725.437.578.388.989.50

bi2.003.501.503.515.502.006.244.102.014.503.410.80

2问题的分析

对问题一建立模型一，对问题二建立模型二

图1题目已知中所给各点分布情况

2.1问题一

网格的横向与纵向为固定的，所以对网格的移动为线性平移，故可以将网格看作先进行横向平移，在进行纵向平移。

效果与直接平移相同，利用将平移分成两步来进行，利用决策变量及约束条件对决策目标进行

图2问题一平移方案示例

2.2问题二

网格的横向和纵向不固定（可以旋转）的情况下进行与问题一相同的求解与讨论。

将其分为两步，第一步以原点为轴，以某角度旋转网格，再沿着网格的自然方向进行平移，与第一问中的操作相同，此两步加起来与进行一步旋转平移效果相同。

图3-1旋转后的网格平移方案示例（第一步操作后）

图3-2问题二网格平移方案示例（第二步操作后）

2.3

设平行移动的距离为xi,yi（xi为横向位移大小，yi为纵向位移大小），设横、纵向位移的小数部分分别为[xi],[yi]。

由题意知，旧井的可利用与否仅与其距离最近的网点的距离有关，而与该网点是哪一个无关。

故横向平移xi的效果与横向平移[xi]的效果相同。

故在此取横向平移距离[xi]=[0，1）。

同理，横向平移距离[yi]=[0，1）。

在对第二个问题的解决中，设网格旋转角度为θ，分析知，当其转过π/2角度后，再转，以前面进行的旋转效果相同，故在此假设θ的范围为[0，π/2）。

2.4修正函数

设修正函数，由于判断点的坐标与其相邻最近网格点的距离

若x≥0.5，x自减1直到其小于0.5

若x≤-0.5，x自加1直到其大于-0.5

3问题假设：

每口井在坐标系中可被看为一个点。

取每个格子的边长为一个单位。

假定网格的横向和纵向是固定的（比如东西向和南北向）。

平移网格时，非连续的平移，为以某一小单位渐进平移。

最后取得一个确定的平移网格方案。

（此处可改进，见模型二）

设在任意点钻井均不存在自然障碍

4符号说明：

表格1-1问题1中的符号说明

点的初始横坐标

点的初始纵坐标

平移后的横坐标

平移后的纵坐标

横坐标平移的距离

纵坐标平移的距离

平移后可利用的旧井数

ai’’

i={1,2,…,12}

表格2-1问题2中的符号说明

点初始横坐标

点初始纵坐标

旋转角度

旋转后横坐标

旋转后纵坐标

横坐标平移

纵坐标平移

移动后横坐标

移动后纵坐标

移后可用旧井数

ai’

a1i’

b1i

i={1,2,…,12}

各井的坐标为pi（i=1，2…n）

移动后各点坐标为pi’

5问题一模型的建立：

目标函数

max

s.t.

0≤dx<1;

0≤dy<1;

问题解决方案：

将移动范围确定后，用穷举尝试法利用计算机编程实现对网格自然方向平移的模拟，存下每种网格移动方案中的旧井利用数，

5.1模型求解：

由约束条件及决策目标建立C++源代码程序求解，得到

横坐标平移数：

x=0.58

纵坐标平移数：

y=0.45

可利用旧井数：

n=4

表格1-2不同步长下利用旧井最多的方案

步长

移动距离

0.1

0.01

0.001

0.0001

…

横向移动距离dx

0.6

0.58

纵向移动距离dy

0.5

0.45

可用旧井数n

由于初始数据给定坐标的精确度为小数点后两位，故取移动的步长为0.01，即可得到与连续变换相同的答案。

5.2讨论：

表格1-3各点移动前后坐标及旧井利用数

旧井标号

P10

P11

P12

横坐标（移动前）

0.50

1.41

3.00

3.37

3.40

4.72

5.43

7.57

8.38

8.98

9.50

纵坐标（移动前）

2.00

3.50

1.50

3.51

5.50

2.00

6.24

4.10

2.01

4.50

3.41

0.80

横坐标（移动后）

1.08

1.99

3.58

3.95

3.98

5.30

6.01

8.15

8.96

9.56

10.08

纵坐标（移动后）

2.45

3.95

1.95

3.96

5.95

2.45

6.69

4.55

2.46

4.95

3.86

1.25

移动后

可用否

否

可

否

可

否

可

否

由上表知，可利用的四口旧井分别为P2、P4、P5、P10号。

5.3平移后与平移前的坐标比较

图4平移前后的坐标示意图

注：

1~12为平移前各点坐标，1'~12’为平移后各点坐标。

6问题二模型而的建立：

目标函数

N=MAX

s.t.

0≤dx<1;

0≤dy<1;

（其中：

a1i’=

×cos（tan

（

）-θ）+dx’

b1i’=

×sin（tan

（

）-θ）+dy’）

问题解决方案：

将移动范围与旋转角度范围确定后，利用计算机编程用穷举尝试法先对网格进行旋转，再实现对网格自然方向平移的模拟，记录下每种网格移动方案中的旧井利用数,取其旧井利用数的最大的方案。

6.1模型求解：

表格2-2不同步长下利用旧井最多的方案

θ的步长

移动距离

0.1

0.01

…

旋转的角度θ

0.8

0.78

横向移动距离dx

0.58

0.25

0.21

纵向移动距离dy

0.45

0.89

可用旧井数n

由约束条件及决策目标编写C++程序求解（程序源代码见附录二），得到

旋转角度θ=0.78rad

横坐标平移数：

dx’=0.21

纵坐标平移数：

dy’=0.89

最大可利用旧井数：

n=6

由分析知，当旋转角度θ的步长为0.01时即可达到最大利用旧井的效果，而当步长进一步缩小时，会加大计算机的预算负荷，使运算时间大大增加，但并不能使方案更优化，所以在实际编程时，可直接取用0.01rad的步长，即可得到预期结果。

6.2移动以后各点情况

表2-3旋转后各点坐标及旧井可利用情况

旧井标号

P10

P11

P12

横坐标（移动前）

0.50

1.41

3.00

3.37

3.40

4.72

5.43

7.57

8.38

8.98

9.50

纵坐标（移动前）

2.00

3.50

1.50

3.51

5.50

2.00

6.24

4.10

2.01

4.50

3.41

0.80

横坐标（移动后）

1.97

3.67

3.40

5.07

6.50

4.97

7.95

6.95

7.01

9.33

8.99

7.52

纵坐标（移动后）

1.96

2.39

-0.15

1.01

2.41

-1.01

2.01

-0.01

-3.00

-1.80

-3.00

-5.22

移动后

可用否

可

否

可

否

可

否

共有6口旧井可被利用，分别为P1、P6、P7、P8、P9、P11号旧井

网格旋转方案为：

网格以原点为轴旋转弧度为0.78rad的角度后，再将网格沿着新的自然方向，横移0.21个单位后，再纵向移动0.89个单位。

6.3移动后各点与移动前的位置示意图

图5移动前后坐标示意图

注：

1*~12*为旋转后的点，1-~12-为移动后的点。

放大的字体为可利用的旧井所在点

7对模型一和模型二的讨论：

显然，方案一为方案二的一种特殊情况。

在方案一的基础上，添加新的参变量θ，将模型进一步完善。

新变量的引入，使网格变化方案更加的灵活多样，有了更多的方案选择，自然可以找到更适合的方案，使旧井利用率得到提升（如图1）

第一问的旧井利用情况第二问中的旧井利用情况

图6两种情况的旧井利用率情况比较

（注：

图中黑色部分为可利用旧井数占总旧井数的百分比）

旧井利用率得到提升，自然费用就相应的减少，见图2，两种网格移动方案的减少费用比较情况。

图7两种情况下的费用比较

注：

1为第一问中的旧井利用方案的费用节省情况，2为第二问中的旧井利用方案的费用节省情况。

显然，网格的移动方案的限制越少，旧井的利用率便会在某限定范围内相应的得到提高。

8对模型的修正：

8.1对模型一的修整

由于非连续的移动网格，所以，通过计算机模拟计算得到的为旧井利用数最高的方案之一。

而实际上的网格移动方案应该是在一个范围内。

通过计算机编程，得到：

dx=[0.58，0.63]

dy=[0.45，0.54]

网格按此范围内的方案移动均可达到最大旧井利用数n=4

（C++程序源代码见附录4）

8.2对模型二的修整

由于非连续的移动网格，所以，通过计算机模拟计算得到的为旧井利用数最高的方案之一。

而实际上的网格移动方案应该是在一个范围内。

通过计算机编程，得到：

dx=[0.21，0.27]

dy=[0.89，0.97]

θ=[0.78，0.79]

网格按此范围内的方案移动均可达到最大旧井利用数n=6

（C++程序源代码见附录5）

9问题三的分析与建模

题目要求给出判定的条件和算法，而不同的算法得出的结论必定不同。

在此，我们沿用欧式距离的误差意义用于判断旧井是否可用。

给出如下三种模型：

9.1模型3-1

利用第二个问题的模型二，将已知各N个点的坐标代入到前面的模型中，求其利用旧井最大数，当n=N时，即为当前的N个旧井均可以被利用。

显然，此模型略显简单，可进一步改进。

9.2．1模型3-2：

引理：

对于任意整数N，存在任意两个正整数N1,N2（零不包括在内），使N=

则N1,N2是该等式的唯一解。

取N点坐标为（x,y）

由坐标转换公式得以适当的角度θ旋转坐标轴后的N点的坐标为（xcosθ-y

sinθ,ycosθ+xsinθ）

必要条件：

任意该两点的距离d满足

（1）

当该条件满足时，我们只能判断该两点满足条件，还需判断其它所有点是否满足该条件。

取任一点X，以步长0.01从0到π/2找寻θ，当xcosθ-ysinθ与接近的整数相差不超过

，此时由于

（1）公式ycosθ+xsinθ也与接近的整数相差不超过

，此时将该θ代入个点，当出现xcosθ-ysinθ与接近的整数相差超过

，则继续寻找θ，若每个都满足ycosθ+xsinθ也与接近的整数相差不超过

，则这些点都能被利用。

9.2.2模型3-2的讨论：

由于第三问用了筛选，即当不满足条件时就更换θ，与第二问相比，运算量将大大减少。

且该充分条件的判定的加入会筛掉许多数值，比第二个模型更完善。

9.3模型3-3：

（该模型将规划的网格图形改变）

经过分析，有如下三种对称图形（包括初始方案）

9.3.1方案一：

形状规则

以概率论方面来说，可利用井的个数n=S*n0，其中S为单个点周围的可行解域（本题以欧氏空间计算），n0为单位面积需打井的个数（为常数），所以改变形状不能优化。

9.3.2方案二：

形状不规则

我们可以根据实际给点的情况拉长图形的一边，使点位于图形交点附近，为度量图形形变的严重程度，再每单元占单位面积的假设下，我们以边长改变△x后该单位变化的面积△S作为度量，△S值越小代表图形形变程度越小，△S值最小的图形为最佳方案。

假设一：

单位面积打一口井，所以一个单元占一个面积。

假设二：

该点放入该点属于的单元计算。

假设三：

该偏离点在该单元的一条边上（为了解题方便，可以这样设，用在边上的点算出的形变度大的话，在中间的点的形变度也一定大）

1.正方形网络：

当其中一个边发生变化△T时

单位边长正方形变化面积为：

△S=1

△T=△T

2.三角形网络：

当三角形的一个边发生变化Δx时：

单位边长三角形面积S=

单位面积三角形边长X=

∴形变度ΔS=

=0.33Δx

3.六边形网络：

当六边形的一个边发生变化时Δt时：

单位边长六边形面积S=1

单位面积六边形边长

形变度

由上可以得出结论，六边形网络为最佳规划。

所以，在此模型中，沿用欧式距离的误差意义，建井布局网格选用正六边形。

参考文献：

[1].姜启源数学模型[M]北京:

高等教育出版社,2003,8

[2].冯杰数学模型原理与案例北京：

科学出版社，2007.1

[3].韩中庚数学建模方法及其应用北京：

高等教育出版社2005，6

附录1：

第一个问题解决时所用的C++程序源代码如下：

#include

voidchange（double&x）

{

for（;x≥0.5;）

{x=x-1;}

}

voidmain（）

{

doublea[12]={0.5,1.41,3.0,3.37,3.4,4.72,4.72,5.43,7.57,8.38,8.98,9.50};

doubleb[12]={2.0,3.50,1.5,3.51,5.5,2.00,6.24,4.10,2.01,4.50,3.41,0.80};

doublex,y,x0,y0;

doublee,c;

inti,n,m;

m=0;

for（x=0;x<1;x+=0.001）

for（y=0;y<1;y+=0.001）

{n=0;

for（i=0;i<12;i++）

{e=x+ai;

c=y+bi;

change（e）;

change（c）;

if（fabs（e）<=0.05&&fabs（c）<=0.05）

n++;

}

if（n>m）{x0=x;y0=y;m=n;}

}

cout<<"x="<

cout<<"y="<

cout<<"n="<

}

附录二：

第二个问题求解时所用的C++程序源代码：

#include

voidchange（double&x）

{

for（;x>=0.5;）

{x=x-1;}

for（;x<=-0.5;）

{x=x+1;}

}

voidmain（）

{

doublea[12]={0.50,1.41,3.00,3.37,3.40,4.72,4.72,5.43,7.57,8.38,8.98,9.50};

doubleb[12]={2.00,3.50,1.50,3.51,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.80};

doublex,y,x0,y0,sta0;

doublee,c;

doublet,s;

doublesta;

inti,n,m;

m=0;

for（x=0;x<1;x+=0.001）

for（y=0;y<1;y+=0.001）

for（sta=0;sta<90;sta=0.1）

{

{n=0;

for（i=0;i<12;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+t*t）<0.05）

n++;

}

if（n>m）

{x0=x;y0=y;sta0=sta;m=n;}

}

cout<<"x="<

cout<<"y="<

cout<<"sta="<

cout<<"n="<

}

附录3：

问题二判断旧井利用情况与旋转网格后各井的新坐标程序

#include

voidchange（double&x）

{

for（;x>=0.5;）

{x=x-1;}

for（;x<=-0.5;）

{x=x+1;}

}

voidmain（）

{

doublea[12]={0.50,1.41,3.00,3.37,3.40,4.72,4.72,5.43,7.57,8.38,8.98,9.50};

doubleb[12]={2.00,3.50,1.50,3.51,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.80};

doublep,q;

doublesta=0.78;

for（inti=0;i<12;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;

cout<<"a["<"<

cout<<"b["<"<

change（p）;

change（q）;

if（sqrt（p*p+q*q）<=0.05）

cout<

-----"<

elsecout<

}

附录4：

修正后模型一的程序源代码

#include

voidchange（double&x）

{

for（;x>=0.5;）

{x=x-1;}

}

voidmain（）

{

doublea[12]={0.5,1.41,3.0,3.37,3.4,4.72,4.72,5.43,7.57,8.38,8.98,9.50};

doubleb[12]={2.0,3.50,1.5,3.51,5.5,2.00,6.24,4.10,2.01,4.50,3.41,0.80};

doublex,y,x0,y0,x1,y1;

doublee,c;

inti,n,m;

m=0;

for（x=0;x<1;x+=0.01）

for（y=0;y<1;y+=0.01）

{n=0;

for（i=0;i<12;i++）

{e=x+ai;

c=y+bi;

change（e）;

change（c）;

if（fabs（e）<=0.05&&fabs（c）<=0.05）

n++;

}

if（n>m）{x0=x;y0=y;m=n;x1=x;y1=y;}

elseif（n==m）

{

if（x>x0）x1=x;

if（y>y0）y1=y;

}

cout<<"x="<

cout<<"y="<

cout<<"m="<

附录五修正后的方案二的C++源程序

#include

voidchange（double&x）

{

for（;x>=0.5;）

{x=x-1;}

for（;x<=-0.5;）

{x=x+1;}

}

voidmain（）

{

doublea[12]={0.50,1.41,3.00,3.37,3.40,4.72,4.72,5.43,7.57,8.38,8.98,9.50};

doubleb[12]={2.00,3.50,1.50,3.51,5.50,2.00,6.24,4.10,2.01,4.50,3.41,0.8

展开阅读全文