综合模型.ppt

1、第七章 综合模型第一节 截断切割,一、问题的提出 这是1997年全国大学生数学建模竞赛的B题，问题如下： 某些工业加工部门，经常需要从一个长方体中加工出一个尺寸,已知，位置预定的长方体（设长方体的对应表面是平面），通常采用截断切割的加工方式，这里“截断切割”是指将物体沿某个切割平面分成两部分，因此，一般情况下，须经过6次截断切割，分别截去原长方体的前后、左右、上下6个方向多余的部分。,设水平切割单位面积的费用是垂直切割单位面积费用的r倍，且当先后两次垂直切割的平面不平行时，因调整刀具需额外费用e，显然，若截去各个多余小块的先后顺序不同，则加工费用不同，试设计一种确定最优加工次序的方法，使加工费

2、用最少。用以下,实例验证你的方法：待加工长方体与成品长方体的长、宽、高分别为10；14.5；19和3；2；4二者左面，前面，底面之间的距离分别为6；7；9（单位cm），垂直切割费用为1元/cm2，r和e的数据有以下4组：a）r=1；e=0； b）r=1.5；e=0；c）r=8；e=0; d）r=1.5；2e15,二、 模型假设及符号说明（1）假设水平工作台接触的长方体底面是事先指定的；（2）假设第一次切割前调整刀具的费用不计；（3）假设加工前后，两长方体对应的表面平行；（4）假设垂直切割费用为1元/cm2，水平切割费用为r元/cm2；（5）假设调整一,次刀具的费用为e；（6）假设总的加工费用为

3、C；（7）假设待加工长方体的长、宽、高分别为a、b、c，为常数，成品长方体的长宽、高分别为A、B、C，也是常数；（8）假设待加工长方体与成品长方体的左面、前面、后面、下面、上面之间的距离为常数：,三、 问题分析及数学模型 根据所给的实际问题，从一个长方体加工出一个尺寸位置预定的长方体，通常情况下，需要经过6次截断切割，如果我们假定这6个切割面分别位于左、右、前、后、下、下面，那么可将它们相应地编号为1、2、3、4、5、6。这样这6个数字的一个全排列代表着一种切割顺序。例如，排列1 2 3 4 5 6代表的的切割顺序为,“左前后右下上”。我们将这6个数字的所有全排列记为集合S，即： ，因此问题转

4、化为求总的切割费用C在S上的最小值问题。 总的切割费用应包含二部分，一部分为加权的切割面积之和，另一部分是刀具调整费用之和，,因此，数学模型为： （7-1）其中： 表示 进行切割后，加工面 的切割面积 ：为相应切割面的权，由题意有： ：为调整刀具的次数。,四、 模型的分析及计算当e=0时（7-1）式成为： 由于集合S为有限集，只有 6！=720个元素，代表着720种切割方式，r给定切割顺序固定则很容易计算出加工费用。比较出所有不同切割方式下的费用，便可得到最小费用下的切割方式。下面给出MATLAB计算程序。,（1）建立可行的加工顺序表； （2）用穷举法求出720种切割方式下的费用； （3）求最

5、小费用及其对应的切割方式。在MATLAB编辑器中建立如下文件。%截断切割ch711%文件名：ch711.ma0=10 14.5 19;a1=3 2 4;d1=6 7 9;r=1;,d2=a0-a1-d1;d=d1,d2; d=d(1,4,2,5,3,6);p=0; for i=1:6 for j=1:6,if(j-i)=0, for k=1:6,if(k-i)*(k-j)=0 for l=1:6,if(l-i)*(l-j)*(l-k)=0, for m=1:6,if(m-i)*(m-j)*(m-k)*(m-l)=0 for n=1:6, if(n-i)*(n-j)*(n-k)*,(n-l)*(

6、n-m)=0, p=p+1;x(p,:)=i j k l m n; end end end end end end end end end,end endf=1 1 2 2 3 3;for p=1:720 o=x(p,:);cost=0;a=a0; for i=1:6 j=o(i);aa=a;aa(f(j)=; if f(j)=3 cost=cost+r*aa(1)*aa(2); else cost=cost+aa(1)*aa(2);,end a(f(j)=a(f(j)-d(j); end c(p)=cost;endminc=min(c),find(c=minc);minx=x(ans,:)执

7、行后输出最小切割费用mink =374,最优切割顺序minx = 5 3 1 6 4 2 5 3 6 1 4 2即当r=1,e=0时最优切割顺序有两条，它们分别是“下前左上后右”和“下前上左后右”，最小切割费用374元。当e0时，模型为（1）式，即,由于垂直切割的方向有两个，因此至少调整一次刀具，而垂直切割共进行四次，故调整刀具的次数至多3次于是额外费用有e，2e， 3e，三种可能。所以我们把全体切割顺序按刀具的调整费用分类，共分三类，同类的刀具调整费用是相等的，故可先分别求出在e=0时每一类的最小费用及加工顺序。再将每一类的最小切割费用,加上相应的刀具调整费用，便得到总 的加工费用。从而比较

8、各类加工顺序 求出整体最优的加工顺序。下面给出 MATBLE计算程序。 （1）计算出三类可行的加工顺序及相应的切割费用表； （2）对每一个最小费用的切割顺序，与刀具调整费用结合考虑； （3）求出总的费用最小时的切割顺序,及调整刀具次数。 %截断切割ch712 %文件名：ch712.mfunction minc,minx1,minx2,minx3= cutorde(a0,a1,d1,r)minc=inf,inf,inf;minx1=;minx2=;minx3=;k1=0;k2=0;k3=0;v1=1 2 3 4 5 6;,%建立可行的加工顺序表 x1 ,x2, x3及 相应的切割费用表for i

9、1=1:6,o1=v1(i1);v2=v1;v2(i1)=;for i2=1:5,o2=v2(i2);v3=v2;v3(i2)=;for i3=1:4,o3=v3(i3);v4=v3;v4(i3)=; for i4=1:3,o4=v4(i4);v5=v4;v5(i4)=;for i5=1:2,o5=v5(i5);o6=v5(3-i5); x=o1,o2 o3 o4 o5 o6; c=cost(x,a0,a1,d1,r);,z=adjustnum(x); switch z case 1,k1=k1+1;x1(k1,:) =x;c1(k1)=c; case 2,k2=k2+1;x2(k2,:) =

10、x;c2(k2)=c; case 3,k3=k3+1;x3(k3,:) =x;c3(k3)=c; end,end end end endendminc=min(c1) min(c2) min(c3);find(c1=minc(1);minx1=x1(ans,:);find(c2=minc(2);minx2=x2(ans,:);find(c3=minc(3);minx3=x3(ans,:);,%求切割费用的子函数function c=cost(x,a0,a1,d1,r)c=0;d2=a0-a1-d1;a=a0;for p=1:6 switch x(p) case 1,c=c+a(2)*a(3);

11、a(1)=a(1)-d1(1); case 2,c=c+a(2)*a(3);a(1)=a(1)-d2(1); case 3,c=c+a(1)*a(3);a(2)=a(2)-d1(2); case 4,c=c+a(1)*a(3);a(2)=a(2)-d2(2);,case 5,c=c+r*a(1)*a(2);a(3)=a(3)-d1(3); case 6,c=c+r*a(1)*a(2);a(3)=a(3)-d2(3); endend%求调整刀具次数的子函数function z=adjustnum(x)z=-1;v0=0;for p=1:6if x(p)5 if x(p)3,v=1; else v

12、=2; end if(v0-v)=0 z=z+1;v0=v; end end end在窗口中输入,a0=10,14.5,19;a1=3,2,4;d1=6,7,9;r=1.5;minc, minx1,minx2,minx3 =cutorde(a0,a1,d1,r)执行后输出minc = 442.5000 456.5000 437.5000minx1 = 3 5 4 1 6 2minx2 =,1 3 5 4 6 2minx3 = 3 1 5 4 6 2 3 5 1 4 6 2每一类的最小费用为：调整一次刀具： ，切割顺序为3 5 4 1 6 2；调整二次刀具： ，切割顺序为1 3 5 4 6 2。

13、调整三次刀具： ，切割顺,序为 3 1 5 4 6 2 3 5 1 4 6 2 在坐标系中，画出三条曲线l1，l2，l3，求出它们的交点，如图7-1。由图可得全体加工顺序中最小费用为,图71,当 时，最小费用为437.5+3e。对应的加工顺序为3、1、5、4、6、2和3、5、1、4、6、2。 当 时，最小费用为442.5+e。对应的加工顺序为3、5、4、1、6、2。,五、问题的进一步讨论 在(1)式中，|S|=6！=720为不 同的切割方式的总数，由于相邻二个平行面的次序交换并不影响总费 用C，因此，由容斥原理可计算出 |S|可减少至426，同时考虑到切割面k到相应的边界面的距离，称为 切割厚

14、度。不失一般性，假定,a1a2，b1b2，c1c2 （7-2） 当（2）式不成立时，假定即可，其它的作类似处理；在（7-2）式成立时，只考虑S的一个子集，即只考虑1在2前、3在4前、5在6前的那些排列，这个子集中含有元素720/23=90。（1）e=0时的优化方法,这时的优化准则为：将 排列成不增序列，相应的指标序列即为切割面的最优序列。在证明中，若假设r=1，否则在垂直方向作一个因子为1/r的变换，其切割面积和的最小化相当于原来的加权面积和最小化，因此在证明,中便不妨假设 同时考虑到e=0的情形，于是（7-1）式成立。 （7-3）其中,下面性质被称为相邻交换原则：设 ，此时假定第j个切割厚度

15、大于或等于第k个切割厚度，将相邻加工面k与j交换后得 则对（7-3）式给出的函数，必有 （7-4） 这个原则也就说明了上述的优化准则了。这主要是因为，从任何一个最优解出发，如果第j个厚度最大，第j,个切割面通过交换总可度调到第一个位置，由于（7-4）式保持为最优解,对于切割厚度次大者，也可以这样做 相邻交换原则的简单证明： 设相应于切割方式 中加工面j，k的切割面积为 ；设相应于切割方式 中加工面k，j的切割面积为 ，由几何知识可得 （7-5）,其中 为切割厚度。 事实上，上式二者均等于在加工方式x（或 ）下加工面k与j被切割之前与切割之后的体积差，由（7-5）式移项得：由于 ，因此有 即,而

16、 与 所含有的各项相同，从而证明了结论（7-4） （2）e0时的优化方法 用表达式（7-1）建立截断切割的数学模型之后，通常可采用分支定界法来求解。每一个分支相当于前面n个面的切割方式已经确定的情形，这时，后面n个面的切割面还可以用多种方式排列，相应费用的下界可以由,二部分和组成，一部分为加权切割面的面积之和，它可以由e=0时的优化准则得到，另一部分为刀具的调整费用的下界，它为0或e。 借助于e=0时的解来求e0时的解，是一种较为简单的方法。当e=0时的最优解 的调整刀具次数q（x）为1时，x使（7-1）式中的C取最小值，当x的刀具调整次数q （x）,为2时，则只能断言，x使（7-1）式在所有

17、刀具调整次数大于1的切割方式中达到最小，为求（7-1）式的整体最小，还应考虑刀具调整次数q（x）为1的切割方式，由前面所述，不妨设（7-2）式成立，这时，可只考虑1在2前、3在4前、5在6前的那些切割方式，在刀具调整次数为1时，在那些切割方式中，加工面1、2、3、4的,的排列方式只可能是1、2、3、4与3、4、1、2，然后将加工面5与6按照相邻交换原则插入到适当的位置，可以得到几个候选的切割方式，它们相当于一定意义下的局部解。这些局部解与上述x组成候选解，比较所有的候选解，便得到e0时的最优切割方式。 由于按排列1、2、3、4的次序,至多有2个单调下降段，可知加工面5与6按相邻单调下降方式插入

18、时，至多可得到3个局部解，对排列3、4、1、2就是这样。所以q（x）=2时，候选解总数至多为3+3+1=7，当x的刀具调整次数q（x）为3时，除了补充考虑刀具调整次数为1的切割方式外，还必须补充考虑刀具调整次数为2的切割方式，它们可由加工面5与6,按照相邻交换原则插到1、3、4、2与3、1、2、4。它们各具有2个单调下降段（二者的插入方式各不超过3），或者二者之一具有3个单调下降段，而另一个则必定只有一个单调下降段（前者的插入方式不超过6，后者的插入方式为1，总数不超过7），所以候选解的总数不超过7+6+1=14，这样就很容易求出e0时的最优切割顺序。,第二节 零件的参数设计 一、 实际问题

19、这是1997年全国大学生数学建模竞赛的A题，问题如下： 一件产品由若干零件组装而成，标志产品性能的某个参数取决于这些零件的参数。零件参数包括标定值和容差两部分。进行成批生产时，标定,值表示一批零件该参数的平均值，容差则给出了参数偏离其标定值的容许范围。若将零件参数视为随机变量，则标定值代表期望值，在生产部门无特殊要求时，容差通常规定为均方差的3倍。,进行零件参数设计，就是确定其标定值和容差。这时要考虑两方面因素：一是当零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大；二是零件容差的大小决定了其制造成本，容差设计得越小，成本越高试通过如下的具体问题给出一般

20、的零零件参数设计方法。粒子分离器某参,数（记作y）由7个零件的参数（记作x1，x2，x7）决定，经验公式为 y的目标值（记作y0）为1.50。当y偏,离y00.1时，产品为次品，质量损失为1000（元）；当y偏离y00.3时产品为废品，损失为9000（元）。零件参数的标定值有一定的容许变化范围；容差分为A、B、C三个等级，用与标定值的相对值表示，A等为1%B等为5%，C等为10%。7个零件参数标定值的容许范围及不同容差等级零件的成本（元）如表71（符号/表示无此等级零件）。,表71,现进行成批生产，每批产量1000个。在原设计中，7个零件参数的标定值为：x1=0.1,x2=0.3,x3=0.1

21、,x4=0.1,x5=1.5,x6=16,x7=0.75；容差均取最便宜的等级。 请你综合考虑y偏离y0造成的损失和零件成本，重新设计零件参数（包括标定值和容差），并与原设计比较，总费用降低了多少。,二、 模型的假设及符号说明 （1）模型的假设 假设各零件的参数为随机变量，且是相互独立的，它们都服从以零件标定值为均值，容差的三分之一为均方差的正态分布。 假设生产1000个批量产品时，每个零件参数只按一种容差等级进行生产，7个零件可以按不同的容差等级进行组合生产。,考虑到题目所给条件，产品参数y偏离y00.1时，产品为次品，质量损失为1000元，y偏离y00.3时，产品为废品，损失为9000元，

22、可假设产品参数y偏离目标值y0造成的单件产品质量损失函数L（y）与（y-y0）2成正比即 （7-6）易见,（2）符号说明 第i个零件的参数，它是随机变量； 第i个零件的标定值，即xi的期望值； 第i个零件的均方差； 第i个零件的容差，即 第i个零件的相对容差，即；,y0：产品质量参数的目标值； y：产品质量参数，它是随机变量；y：y的均方差；ai：标定值xi0的取值下限（已知的）bi：标定值xi0的取值上限（已知的）x0=（x10，x20，x70）标定值向量t=（t1，t2，t7）：相对容差向量,三、问题的分析及数学模型 （1）问题的分析 问题的目标函数是总费用函数最小，而总费用是由产品参数偏

23、离目标值引起的质量损失费用和产品的成本费用两部分组成。一般来说，标定值设计不合理或容差设计得太大，会使产品参数远离目标值，造成质量损失；而容差设计得太小，又会增加零,零件制造成本，综合考虑产品质量 损失费和成本费是问题的关键。 由题目所给条件知，单件产品的零件成本取决于相对容差等级ti，第I 种零件的成本记作ci（ti），于是单件产品的成本，即7个零件的总成本为 （7-7）产品参数y由x1，x2，x7决定，,记作y=f（x1，x2，x3，x4，x5，x6，x7），而xi是随机变量，且xiN（xi0，i2），所以产品参数y也是随机变量。进行成批生产时，平均每件产品的质量损失费用应该用损失函数L（

24、y）的期望来度量，它取决于零件参数标定值x0和容差t，记为：,由（7-6）式为了得到 的简单表达式，在 处对 作Taylor展开，并略去二阶及二阶以上项，有,其中 于是所以,通过以上分析讨论可知，成批生产时平均每件产品的总费用可表示为 （2）数学模型 问题的目标函数是总费用函数Z（x0，t）达到最小，且约束条件是标,定值x0落在给定的容许范围内，所以该问题的数学模型为： （7-8）,四、模型的分析及求解 1. 模型的分析 注意到在模型（7-8）中，决策变量x0是连续的，t是离散型的，而且t的取值共有23332=108种，对每种固定的t值，求解如下一系列子问题： （7-9）,得到最优解 和最优值

25、 ，然后对108个t比较 （7-10）得到全局最优解 和最优值 。 子问题（7-9）是非线性规划问题，它的计算量是很大的，为了减少计算子问题的次数，在程序设计时，每解出一个子问题，由（7-10）式算出当前的最优值 ，对待求解的子问题，先判断 是否成立，若成立，则该子问题不必再解。,2. 模型的求解 按照上面对模型（7-8）的分析，下面设计在MATLAB中计算模型（7-8）的程序。为了方便起见，记,分别建立函数文件v.m，u.m，fy.mfunction v=fun(x)v=1-0.36*x(2)0.56*x(4)(-0.56);function u=fun(x) u=1-2.62*v(x)1.

26、5*(x(4)/x(2)1.16;function fy=fun(x)fy=174.42*(x (1)/x(5)*(x(3)/(-x(1)+x(2)0.85*sqrt(u(x)/(x(6)*x(7);,分别建立函数u(x)对x2,x4的偏导数函数文件uu2.m,uu4.mfunction uu2=fun(x)uu2=0.792288*v(x)0.5*x(2)(-1.6)*x(4)0.6+3.0392*v(x)1.5*x(2)(- 2.16)*x(4)1.16; function uu4=fun(x)uu4=(-0.792288)*v(x)0.5 *x(2)0.6*x(4)(-0.4)-3.03

27、92*v(x)1.5*x(2)(- 1.16)*x(4)0.16;,分别建立函数fy(x)对x1,x2,x7的偏导数函数文件y11.m,y21.m,y71.mfunction y11=fun(x)y11=fy(x)*(1/x(1)+0.85/(x(2)-x(1);function y21=fun(x)y21=fy(x)*(-0.85)/（x(2)-x(1)+0.5*(1/u(x)*uu2(x);function y31=fun(x)y31=fy(x)*(0.85/x(3);,function y41=fun(x)y41=fy(x)*0.5*(1/u(x)*uu4(x);function y51

28、=fun(x)y51=fy(x)*(-1/x(5);function y61=fun(x)y61=fy(x)*(-0.5/x(6);function y71=fun(x)y71=fy(x)*(-0.5/x(7);,建立子问题（7-9）的目标函数文件ch721.mfunctionch721.g=ch721(x,p2,p3,p4,p6,p7)ch81=100000*(fy(x)-1.5)2+(100000/9)*(y11(x)*x(1)*0.05)2+(y21(x)*x(2)*p2)2+(y31(x)*x(3)*p3)2+(y41(x)*x(4)*p4)2+(y51(x)*x(5)*0.1)2+(

29、y61(x)*x(6)*p6)2+(y71(x)*x(7)*p7)2;g=-x(1);,用MATLAB语言设计计算模型（7-8）的程序，可存为M文件ch722.mtic %启动计时器p2=0.1;p3=0.1;p4=0.1;p6=0.1;p7=0.05; %原设计方案的相对容差x0=0.1,0.3,0.1,0.1,1.5,16,0.75; %原设计方案的标定y0=fy(x0); %原设计方案产品的参数值,Q0=ch81(x0,p2,p3,p4,p6,p7); %原设计方案的质量损失费Z0=Q0+200; %原设计方案的总费用Zmin=Z0; %原设计方案的总费用作为优化初值t1(2)=0.05

30、; %对各等级零件相对容差赋值t2(2)=0.05;t2(3)=0.1;t3(1)=0.01;t3(2)=0.05;t3(3)=0.1;,t4(1)=0.01;t4(2)=0.05;t4(3)=0.1;t5(3)=0.1;t6(1)=0.01;t6(2)=0.05;t6(3)=0.1;t7(1)=0.01;t7(2)=0.05;c2(2)=50;c2(3)=20;c3(1)=200;c3(2)=50;c3(3)=20;c4(1)=500;c4(2)=100;c4(3)=50;c6(1)=100;c6(2)=25;c6(3)=10;c7(1)=100;c7(2)=25;,for i2=2:3fo

31、r i3=1:3for i4=1:3for i6=1:3for i7=1:2C=75+c2(i2)+ c3(i3)+ c4(i4)+ c6(i6)+ c7(i7); %计算单件产品的成本p2=t2(i2); p3=t3(i3); p4=t4(i4); p6=t6(i6); p7=t7(i7);,Zminif C=Zminx0=0.1,0.3,0.1,0.1,1.5,16,0.75; %优化的初始值v1b=0.075,0.225,0.075,0.075,1.125,12,0.5625 %变量x的下界vub=0.125,0.375,0.125,0.125,1.875,20,0.935 %变量x的上

32、界options= ; %优化中的参数均采用缺省值,jacob= ; %将空阵赋给雅可比阵x,options=constr(ch81,x0,options,vlb,vub,jacob,p2,p3,p4,p6,p7); %计算非线性规化问题（7）Z=C+cs97(x,p2,p3,p4,p6,p7); %当t值给定时的最小总费用if Z=ZminZmin=Z; %当前的最小总费用j2=i2;j3=i3;j4=i4;j6=i6;j7=i7;Xmin=x;,Cmin=C;Qmin=ch81(x,p2,p3,p4,p6,p7);y=fy(x)endendendendendendend,T=0.05,t2

33、(j2),t3(j3),t4(j4),0.1,t6(j6),t7(j7) %输出最优设计方案的相对容差Xmin %输出最优设计方案的标定值Cmin %输出最优设计方案的成本费Qmin %输出最优设计方案的质量损失费Zmin %输出最优设计方案的总费用,y %输出最优设计方案产品的参数值 Q0 %输出原设计方案的质量损失费 Z0 %输出原设计方案的总费用 y0 %输出原设计方案产品的参数值 toc %关闭计时器，并显示 程序运行时间,注意：以上建立的一系列m文件一定要存盘。然后在MATLAB下键入ch722可得数值结果为：T=0.0500 0.0500 0.0500 0.1000 0.1000

34、0.0500 0.0500Xmin=0.0750 0.3750 0.1249 0.1250 0.1264 15.0567 0.7740Cmin= 275 Qmin= 465.0495,Zmin= 740.0495y= 1.4969Q0=6.2946e+003Z0=6.4946e+003y0= 1.7256elapsed_time= 134.0100,第三节 钻井布局,一、问题的提出 这是1999年全国大学生数学建模 竞赛的B题，问题如下： 已知勘探部门在某地区找矿，初步勘 探时已在若干位置上钻井，取得了地质资 料，进行统勘探时期后，要在一个区域内 按横纵等距的网格来布置井位，进行“撤 网式”全

35、面钻探。由于钻一口井的费用很 高，因此应尽量利用旧井，少打新井。,设平面上有n个点 ，其坐标为 表示已有的n个井位。新布置的井位是一个正方 形网格N的所有结点（指每个格子都是正方形网 格，结点是纵与横线的交叉点）。假定每个格 子边长（井位的纵横间距）都是非曲直单位 。 整个网格是可以在平面上任意移动。若一个已 知点 与某个网格点 的距离不超过给定误差 e（=0.05单位），则认为 处的旧井资料可以 利用，不必在结点 处打新井了。 为进行辅助决策，勘探部门要求我们研究,如下问题。 （1）假定网格的横向和纵向是固 定的（比如东西和南北向），并规定 两点间的距离为其横向距离（横坐标 之差绝对值）及纵

36、向距离（纵坐标之 差绝对值）得最大值。在平面上平行 移动网格N，使可利用的旧井数尽可 能大。试提供数值计算方法，并对下 面的数值例子（见下表）用计算机进,行计算。 （2）在欧氏距离的误差意义下，考虑 网格的横向不固定（可以旋转）的情 形，给出算法及计算的结果。 （3）如果有n口井，给出判定这些并 均可以利用的条件和算法（你可以任意 选定一种距离）,数值例子 n=12个点的坐标如下表,二、模型假设及符号说明 （1）假设n口井都可以进行生产， 没 有废弃的井。 （2）假设在一定的范围内，任何地 点钻新井都可以产油。 （3）假设钻新井的费用只与井口数 有关，与其它条件无关。 （4）假设在某直角坐标系

37、原点O最近 的结点坐标为,（6）假设 表示 与网格结点 的距离误差。 （7）假设 表示 的取整。 （8）假设网格的旋转角度为 三、 问题分析及数学模型 1. 网格横向和纵向均固定的情况 对于给定的直角坐标系 ，已知点 的坐标为 而在网格N中， 由假设已知，离坐标原点o最近的结点为,，于是 且对于网格N 中的任 一结点可表示为 ，其中 都为整 数 ，那么问题转化为对于给定的网格N 的 设计参数 如何计算可利用旧井数目问题。 已知点 与结点 的距离误差是沿坐标 轴方向的，由题意可知 与 的距离应小 于或等于，即要求正方形邻域。 中存在结点 时，旧井 才可以利,用，否则不能利用，若设 表示旧 井 可

38、以利用，则此时 ，否则 上述过程可以用模型表示为： 其中INT表示向 方向取整。 若设可利用井数为 则 问题归结为当 时求,的极大值问题。将 的取值离散化，对 满足（7-11）式的 进行计数，可取得 （7-12）式的值。 2. 对于网格横向不固定的情况 首先我们考虑用欧氏距离表示误差 而网格N不旋转的情况，由题意可 知，当且仅当园形领域 内存在结点， 时，结点 才可以 利用，此时仍记 否则 而,此时的圆形邻域含于问题 1 的正方形邻 域之中，而唯一可能含于正方形邻域中 的结点是 因此，旧井 才可以利用的条件为 此时， 所以旧井可以利用的个数为,3. 对于网格可以旋转的情况 原坐标系 旋转一角度

39、 得到新 坐标系为 其中 轴与 轴的夹角为 。由坐标旋转公式，可得点 在 新坐标系下的坐标为： 此处坐标原点o不一定是网格N的结点。 假设在坐标系 中，网格N中距,离原点最近的结点为 则 于是网格N 的设计参数就是 对 的 取值范围 离散化，对每个 的值， 运用公式（7-13）计算出各点 的新坐 标，然后再用上面的方法，进行搜索求 解。 4. 对于 n口井均可以利用的情况 不妨假设指定的井就是 则 符合题意的网格N ，应满足,其中 及 满足（7-13）式，z为整 数集合。 而对于（7-14）式有解的条件为： （7-15）,（7-16） 由于它们的形式一样，所以这里我们 只研究其中一组，如（7-

40、15）。对于 （7-16），由于 故而： ；其中 当 时，有 。 因此，整数变量 至多有两个可,能值 及 。 参数s的取值范围，可从下面的命 题中看出： 命题：当且仅当如下二者之一成立时 ，不等式组（7-15）有解 。 （） （） 证明： 若（）成立，则可取,若（）成立，则可取 反之，设（7-15）有解 及s。若 ，则可知 ，事实上 ，若 ，则自然有 ；否则 （注意）由 ；可知 ； 即 ， 因此,从而条件（）成立，同理，若 ， 则可断言 ，从而条件（）成立 。基于上述命题，可得判定不等式组 （7-13）-（7-14）是否有解的算法。 基本步骤是： 过程 ；对 计算：,若 或 输出“有解”，否则

41、输出 “无解”。 算法： 对 的取值范围 离散化： 令j=0。 (2)对 ，按公式（3）算出 的坐标 。 3）执行过程 ，若无解，,则转(5)：若有解则记录： （4）执行过程 ，若无解，则 转(5) ：若有解则记录： 转(6)。 （5）若j=m，则终止（问题无解）； 否则令j=j+1，转(2)。,（6）输出 终止（问题无解）。 四、模型分析及计算 由于问题（7-11）中要只要求沿着 横向和纵向移动网格，因此可按步长 为0.025（或以更小的步长）单位进行 移动，由于给定的以知井位均位于第 一条限，所以在平移时把网格想负方 向移动即可。 问题（7-12）中要求网格的横向不 固，可以旋转搜索，此时

42、可把直角坐,标系转化为极坐标系，点的位置有极 角和极径给出，利用 对角度 以 步长进行搜索，从 而转化为问题（7-11）进行求解。 下面给出问题（7-11），（7- 12） 的MATLAB计算程序。,%钻井布局ch73 %文件名：ch73.m dianxy=0.5,1.41,3,3.37,3.4,4.72,5.43,7.57,8.38,8.98,9.8; 2,3.5,1.5,3.51,5.5,2,6.24,4.1,2.01,4.5,3.41,0.8 wce=0.05; jnum,jno,ydian=zjing(dianxy,wce) jnum,jno,ydian,jiao=zjing2(dia

43、nxy,wce),functionjnum,jno,ydian,jiao=zjing2(dianxy,wce) jnum=0;ttdiannum=size(dianxy,2); diany=dianxy(2,:); %开始搜索 for i=0:wce/2:1 dianxy(1,:)=dianxy(1,:)+I; dianxy(2,:) for j=0:wce/2:1,mjnum=0;mjno=; dianxy(2,:)=dianxy(2,:)+j gedianxy=tound (dianxy); distance=dianxy-gedianxe; for k=1:ttdiannum if ab

44、s(distance(1,k),=wce,end end if mjnumjnum jnum=mjnum; jno=mjno; ydian=-I,-j; end end end,functionjnum,jno,ydian,jiao=zjing2(dianxy,wce) jnum=0;ttdiannum=size(dianxy,2); %直角坐标化为极坐标 dianjb=zeros(2,ttdiannum); newdianxy=zeros(2,ttdiannum); dianjb(1,:)=sqrt(sum(dianxy.2); dianjb(2,:)=atan(dianxy(2,:)./d

45、ianxy(1,:); %开始搜索 for h=0:pi/360:pi/2,%逆时针旋转h度角 fujiao=dianjb(2,:)-h； newdianxy(1,:)=dianjb(1,:).*cos(fujiao); newdianxy(2,:)=dianjb(1,:).*cos(fujiao); %在旋转后的直角坐标中搜索 newdianxy=newdianxy(2.:); for i=0:wce/2:1 newdianxy(1,:)=newdianxy(1,:)+i; %x zhou come back newdianxy(2,:)=newdiany;,for j=0:wce/2:1

46、mjnum=0;mjno=; newdianxy(2.;)=nwedianxy(2,:+j); gedianxy=round(newdianxy); distance=newdianxy-gedianxy; for k=1:ttdiannum If distance(1,k)2+distance(2,k)2= wce2; mjnum=mjnum+1 mjno=mjno,k;,end end if mjnumjnum jnum=mjnum; jno=mjno; ydian=-I,-j; jiao=h*18./pi; end end end,end 执行后输出 问题一 可利用井数 jnum=4 井

47、号：2 4 5 10 移动值为-0.3500，-0.1500 问题二,可利用井数 jnum=6 井号：jn0=1 6 7 8 9 11 移动值为-0.1000 0.2000 旋转角度为,第四节 钢管订购和运输一、问题的提出 这是2000年全国大学生数学建模竞赛的B题，问题如下： 要建设一条 的输送天然气的主管道，如图82所示。经筛选后可以生产这种主管道钢管的钢有 。,图中粗线表示铁路，单细线表示公路，双细线表示要铺设的管道（假设沿管道或者原来有公路，或者建有施工公路），圆圈表示火车站，每段铁路，公路和管道旁的阿拉伯数字表示里程（单位km）。 为方便计，1km主管道钢管称为1单位钢管。,图72,

48、一个钢厂如果承担制造这种钢管，至少需要生产500个单位。钢厂 在指定期限内能生产该钢管的最大数量为 个单位，钢管出厂销价1单位钢管为 万元，如下表：,1单位钢管的铁路运价如下表：,1000km以上每增加1至100km运价增加5万元。 公路运输费用为1单位钢管每公里0.1万元（不足整公里部分按整公里计算）。 钢管可由铁路、公路运往铺设地点（不只是运到点 ，而是管道全线）。 （1）请制定一个主管道钢管的订购和运输计划，使总费用最小（给出总费用）。,（2）请就1）的模型分析：哪个钢厂钢管的销价的变化对购运计划和总费用影响最大，哪个钢厂钢管的产量的上限的变化对购运计划和总费用的影响最大，并给出相应的数

49、字结果。 （3）如果要铺设的管道不是一条线，而是一个树形图，铁路、公路和管道构成网络，请就这种更一般的情形给出一种解决办法，并对图73按（1）要求给出模型和结果。,图73,二、模型的假设及符号说明 （1）模型的假设 钢管运输时，不计换车、转站的时间和费用。 不计装卸费用。 不计运输时由于运输工具出现故障等意外事故引起工期延误造成损失。（2）符号说明 ：一个单位钢管从 到 的最小运价（含购价）；,：从 到 的运 量； ： 得到的钢管向 的左边（西边）铺设的钢管数量； ：得到的钢管向 的 右边（东边）铺设的钢管数量； ： 到 的长度； ：钢厂 生产钢管的最大数量（已知的）；,三、问题的分析及数学模

50、型 钢管的订购和运输方案是直接影响工程费用的主要原因，因此，选取费用最短的路线运送货物，合理的订购计划是决定该工程费用的重要因素，首先利用图论的方法，来确定从钢管生产厂家到施工结点的费用最小路线，然后建立工程费用的优化模型，从中优化出最佳购运方案。,首先利用E.W.Dijkstra的最短路算法计算出一个单位钢管从到的最小费用系数阵，MATLAB程序如下：% 钢管订购和运输ch741% 文件名：ch741.m% ch741 函数是将距离转化为费用function f=ch741(ll)rr=ll;,if rr=0 real1=0; elseif rr0 real1=20;elseif rr300

51、.5 elseif rr450.5 elseif rr500.5 ,% 钢管订购和运输ch742% 最小费用矩阵 c% 文件名：ch742.m% 输入原始数据 A= ; for i=1:39 for j=1:39 A(i,j)=inf; endend for j=1:39 A(j,j)=0;end,A(1,2)=104; A(2,3)=301; A(3,4)=750; A(4,5)=606;A(5,6)=194; A(6,7)=205; A(7,8)=201; A(8,9)=680;A(9,10)=480; A(10,11)=300; A(11,12)=220;A(12,13)=210; A(

52、13,14)=420; A(14,15)=500;A(15,22)=20; A(15,39)=20; A(14,38)=30;A(14,21)=110; A(13,37)=62; A(12,35)=10;,A(11,34)=10; A(10,32)=70; A(9,31)=42; A(8,30)=12; A(7,16)=31; A(7,29)=10;A(6,28)=5; A(5,27)=10; A(4,26)=600; A(3,24)=2; A(2,23)=3; A(23,25)=450; A(24,25)=80; A(25,26)=1150;A(26,30)=1100; A(17,30)=1

53、200; A(30,16)=202;A(16,29)=20; A(29,28)=195; A(28,27)=306;,A(30,31)=720; A(31,18)=690; A(31,32)=520;A(32,33)=170; A(33,19)=690; A(33,36)=160; A(33,34)=88; A(34,20)=462; A(33,36)=160;A(36,35)=70; A(36,37)=320; A(37,38)=160;A(38,21)=70; A(38,39)=290; A(39,22)=30;A(2,1)=104; A(3,2)=301; A(4,3)=750;,A(5

54、,4)=606; A(6,5)=194; A(7,6)=205;A(8,7)=201; A(9,8)=680; A(10,9)=480;A(11,10)=300; A(12,11)=220; A(13,12)=210;(14,13)=420; A(15,14)=500; A(22,15)=20;A(39,15)=20; A(38,14)=30; A(21,14)=110;A(37,13)=62; A(35,12)=10; A(34,11)=10; A(32,10)=70; A(31,9)=42;,A(30,8)=12;A(16,7)=31; A(29,7)=10; A(28,6)=5;A(27

55、,5)=10; A(26,4)=600; A(24,3)=2;A(23,2)=3;A(25,23)=450; A(25,24)=80; A(26,25)=1150;A(30,26)=1100; A(30,17)=1200;A(16,30)=202;A(29,16)=20; A(28,29)=195; A(27,28)=306;A(31,30)=720; A(18,31)=690; A(32,31)=520;,A(33,32)=170; A(19,33)=690; A(36,33)=160;A(34,33)=88; A(20,34)=462; A(36,33)=160;A(35,36)=70; A(37,36)=320; A(38,37)=160;A(21,38)=70; A(39,38)=290; A(22,39)=30;A;,%将火车道长转化为费用 H=;HH=; % 先计算任意两火车站之间的最短距离 for sw=16:39 ; % H为距离矩阵，为火车费用矩阵 L=; % 为临时距离矩阵 for i=1:40 h(i)=0; end for i=16:39,L(i)=A(sw,i); end T=L(16);s=16; % T 为最小数的暂时空间 for i=17:39 if TL(i); T=L(i); % 找最小数 s