自来水管道问题

1、自来水管道连接规划模型 摘要在我们的日常生活中，有这样一个问题，那就是需要通过自来水管道将自来水运输至各个用户处，本文将着重分析讨论自来水管道连接规划问题，即在自来水管道铺设过程中在绕开障碍物的前提下的最优路径问题，使自来水管道将各个供水点用最短路径连接。采用合理的方法排除障碍区域中的点，是对管道连接的效率、能耗、可行性起到决定性作用，是一个非常实际的问题。本文采用面积分析的方法，提供一种解决位于障碍区域中点的判定方法，在二维坐标系上标定各点，障碍区域用由阴影覆盖的凸多边形表出，通过对点坐标之间的向量运算判定各点是否位于阴影区域，在剔除障碍区中的点后，采用Kruskal避圈算法生成最优路径，对

2、于通过障碍区域的线段，采用将其权值设定为（Inf）的处理方法，最终通过Matlab2009a编程、绘图，给出管道最优连接方案，解决本问题。最终通过Matlab2009a编程实现。Kruskal避圈算法是Kruskal在1956年提出的最小生成树算法，它的思路很容易理解。Kruskal算法每次选择n-1条边，所使用的贪婪准则是：从剩下的边中选择一条不会产生回路的具有最小权重的边加入以选择的边的集合中。最后我们对模型的可行性，合理性，科学性进行了详细分析，对模型的优劣点进行阐述，得出了对模型的评价以及推广。 关键词：管道连接 面积法 障碍点筛选 有效线段的筛选 矩阵的运算 Kruskal算法 权值

3、 最小生成树目录摘要2一.问题重述4二问题分析7三模型假设83.1 基本假设83.2符号和变量的说明8四模型建立8五模型求解105.1 筛选有效用户105.2有效线段的筛选115.3利用Kruskal算法求最小生成树11六模型检验与评价12七.模型推广12八参考书目12九附录13附录一：分工13附录二：问题重述附录13附录三：Matlab代码1（做图）16附录四：求解最小生成树17附录五：计算生成树长度40附录六：Matlab代码277一. 问题重述自来水是人们日常生活中不可缺少的生活要素，然而自来水管网的组建却有很多问题需要解决。一般来说，我们假设管网中任意两个用户之间存在直线段相连，但是在

4、连接过程中，有些区域是必须绕开的，这些必须绕开的区域我们称为障碍区域。表1给出了若干个可能的用户的地址的横纵坐标，可能的用户的含义是：如果用户的地址不在障碍区域内，那么该用户就是需要使用自来水的用户（即有效用户），否则如果用户的地址在障碍区域内，那么该用户就是无效用户（即不要将该用户连接在网络中）。表2-表5是分别是4个障碍区域必须要覆盖的点的坐标，而对应障碍区域就是覆盖这些要覆盖的点的最小凸集。(1)请您判定表1中那些用户为有效用户。(2)请设计一个算法将有效用户连接起来，并且连接的距离总和最小。表1若干个可能的用户的地址的横纵坐标可能的用户的序号可能的用户横坐标可能的用户纵坐标1.0000

5、95.012958.27922.000023.113942.34963.000060.684351.55124.000048.598233.39515.000089.129943.29076.000076.209722.59507.000045.646857.98078.00001.850476.03659.000082.140752.982310.000044.470364.052611.000061.543220.906912.000079.193737.981813.000092.181378.332914.000073.820768.084615.000017.626646.10951

6、6.000040.570656.782917.000093.547079.421118.000091.69045.918319.000041.027060.286920.000089.36505.026921.00005.789141.537522.000035.286830.499923.000081.316687.436724.00000.98611.500925.000013.889176.795026.000020.276597.084527.000019.872299.008328.000060.379278.886229.000027.218843.865930.000019.88

7、1449.831131.00001.527421.396332.000074.678664.349233.000044.509632.003634.000093.181596.009935.000046.599472.663236.000041.864941.195337.000084.622174.456638.000052.515226.794739.000020.264743.992440.000067.213793.338041.000083.811868.333242.00001.964021.256043.000068.127783.923844.000037.948162.878

8、545.000083.179613.377346.000050.281320.713347.000070.947160.719948.000042.889262.988849.000030.461737.047750.000018.965457.514851.000019.343145.142552.000068.22234.389553.000030.27642.718554.000054.167431.268555.000015.08731.286356.000069.789838.396757.000037.837368.311658.000086.00129.284259.000085

9、.36553.533860.000059.356361.239561.000049.655260.854062.000089.97691.576063.000082.16291.635564.000064.491019.007565.000081.797458.691866.000066.02285.758167.000034.197136.756868.000028.972663.145169.000034.119471.763470.000053.407969.266971.000072.71138.407972.000030.929045.435573.000083.849644.182

10、874.000056.807235.325075.000037.041415.360676.000070.274067.564577.000054.657169.921378.000044.488072.750979.000069.456747.838480.000062.131055.484281.000079.482112.104782.000095.684345.075483.000052.259071.588384.000088.014289.284285.000017.295627.310286.000097.974725.476987.000027.144786.560388.00

11、0025.232923.235089.000087.574280.487290.000073.730690.839891.000013.651923.189492.00001.175723.931393.000089.38984.975494.000019.91387.838495.000029.872364.081596.000066.144319.088797.000028.440984.386998.000046.922417.390099.00006.478117.0793100.000098.833599.4295表2障碍区域1必须要覆盖的点的坐标顶点序号顶点的横坐标顶点的纵坐标13

12、.2060 12.9166217.457119.337734.7576 20表3障碍区域2必须要覆盖的点的坐标顶点序号顶点的横坐标顶点的纵坐标150 30253.746548.4490346.922257.1195433.320739.8050543.112356.3187表4障碍区域3必须要覆盖的点的坐标顶点序号顶点的横坐标顶点的纵坐标154.698270253.746590346.922280表5障碍区域4必须要覆盖的点的坐标顶点序号顶点的横坐标顶点的纵坐标190752809537080二问题分析建立模型要达到的目的就是节省管道，即在满足每个有效用户用水的情况下，使得铺设的管道最短。因此，

13、自来水的管道问题可以看做是一个最优化问题，目标函数是求铺设的管道最短。由实际可知不是每两个用户之间都可以用直线相连，必须绕开一些障碍物也就是所谓的障碍区，所以我们应该首先要解决的就是找出这些障碍区域，然后再判断所给出的点是否位于障碍区域内，这样就筛选出了有效用户。接下来就是要把剩下的点用直线连接起来，通过障碍区域的线段视为无效线段把其剔除，筛选出有效线段。最后就是计算出这些有效线段的总和。三模型假设3.1 基本假设1. 假设任意两个用户之间均可用直线连接；2. 文中给出所有点的坐标值准确无误；3. 障碍区域就是障碍顶点围成的凸多边形区域；4. 有效用户都能通过自来水管道获得自来水供应；5. 要

14、保证在任意两点间线段不过障碍区的情况下，求解连接形成的最短路径；3.2符号和变量的说明表6 论文符号说明符号含义X记录100个用户点的坐标信息A障碍区1的各顶点坐标信息B障碍区2的各顶点坐标信息C障碍区3的各顶点坐标信息D障碍区4的各顶点坐标信息SIGN记录各用户点是否在障碍区，若在对应位置记为1；若不在，则对应位置记为0INSIGN记录在障碍区的用户点的序号n记录保留用户点的个数NUM记录任意两用户点之间可用线段连接起来且不过障碍区的线段DIS记录不在障碍区各用户点之间可用不过障碍区线段连接的线段的长度EE记录生成的最小生成树的各点及各线段信息sum表示产生的最小生成树中所有管道的总长四模型

15、建立 5.1.问题一的模型建立问题一是判断这100个点中哪些点属于有效点，即有效用户。首先利用matlab做出这一百个点的相应位置的图，其代码见附录三做出此图，分析可知：要求出哪些用户为有效用户，可用面积法对其进行筛选。这样就先得根据障碍区域的顶点坐标求出每个障碍区域的面积，然后求出各用户点与各障碍区域任意两个顶点所围成的三角形面积之和，比较面积，若两面积相等，则该点在障碍区域内，视为无效点，即无效用户，否则用户点不在障碍区域内，为有效用户。根据障碍区的顶点坐标，可做出相应的图形，代码见附录三，图如下：五模型求解5.1 筛选有效用户用面积法确定是否为有效点。面积法的原理：确定各障碍区的面积以及

16、用户点与各障碍区任意两个定点构成的三角形的面积之和，比较上面两个面积，若相等，则该用户点在障碍区内为无效用户，否则，用户点不在障碍区内为有效用户。运用向量的方法求解障碍区面积S若障碍区是三角形，对应各顶点坐标分别为(x1,y1)，(x2,y2), (x3,y3)。则a=(x2-x1,y2-y1),b=(x3-x1,y3-y1)。由于三角形面积S=|a|*|b|*sin<a,b>/2,向量a,b外积的模长|a×b|=|a|*|b|*sin<a,b>则有S=|a×b|/2;若障碍区为五边形，对应点为(x1,y1)，(x2,y2), (x3,y3), (x

17、4,y4)，(x5,y5)。则划分成三个三角形，各三角形的顶点分别为(x1,y1)，(x2,y2), (x3,y3)；(x3,y3), (x4,y4)，(x5,y5)；(x1,y1)，(x3,y3), (x5,y5)。再用求三角形面积的方法求解即可。筛选完毕的结果如下：INSIGN = 4 23 36 99n = 96 所以在障碍区的点的序号分别为：4 23 36 99。 无效用户的信息为：(4.0000,48.5982,33.3951);(23.0000,81.3166,87.4367); (36.0000,41.8649,41.1953); (99.0000,6.4781,17.0793)

18、;有效用户的个数是：96。5.2有效线段的筛选 已筛选出有效用户，就要求出有效用户之间以最短的线段线段相连，但是这些线段必须是有效线段，若两用户之间以线段相连了，但是这条线段通过了障碍区域，此时，这条线段就是无效线段。此时需要筛选出有效线段，首先要求出任意两个有效用户之间的直线与过各障碍区域任意两个顶点之间的直线的交点坐标，然后用向量法判断该交点是否在两用户的线段上和障碍区顶点为端点的线段上，若在，则为无效线段，否则为有效线段。 运用矩阵的方法求解两直线之间的交点坐标 如果任意两个有效用户点的坐标分别为A、B，同一障碍区任意两个顶点坐标为M、N。则由解线性方程组的方法有，运用Matlab求解该

19、线性方程组=A。5.2.2运用向量法判断线段是否为有效线段若求得的交点坐标为P(x,y)，则通过向量关系PM=PN,可以求的。若0,则该线段为有效线段；若<0,则要考虑向量关系PA=PB,若0，则该线段为有效线段，否则，该线段为无效线段，生成的矩阵见附录四，在m矩阵中存储。5.3利用Kruskal算法求最小生成树筛选出有效线段之后就要求出这些有效点之间怎样用线段连起来，可能连接的方式有很多种，但为了更好的节省资源，我们就要求出最短的路线，因此就涉及到最优化问题。考虑在连接96个用户的自来水管道问题还要考虑所有的管道之和最短，就要用到最优化问题，利用Kruskal算法设计matlab程序进

20、行生成最小生成树，再将各边连接起来求出最短路径长度，画出连接图形。一共有96个有效的用户，因此我们只考虑这96个点，上面已经筛选出了有效线段，这些点只有连接和不连接，我们可以先求出邻接矩阵，再根据有效点进行修改，将无效线段的权值设为（Inf），有效线段的权值设为其距离，这样就可以得到一个新的邻接矩阵DIS，见附录四。然后用冒泡法对所有有效线段进行排序，用Kruskal算法进行最小生成树的计算，代码见附录五，最后求得的最小生成树如下图，求得的结果为sum =653.0196 六模型检验与评价 修改最小生成树任意一条边的连接，得出生成树的总和大于前面求出的sum,所以，求出的最小生成树为最佳方案。

21、此生成树已经是最优的，并且不经过障碍区，因此，此模型就可以应用于具体的实际生活中。 七.模型推广该模型可以应用到修建连接若干个地点之间的公路长度和最短的问题中，也可以应用到若某地发生事故，警察以最最快的速度到达出事地点。 八参考书目1 汪晓银，周保平，数学建模与数学实验，北京：科学出版社，20102 薛定宇，陈阳泉，高等应用数学问题的Matlab求解，清华大学出版社，20083 赵静，但琦等，数学建模与数学实验，高等教育出版社，20084 贺超英，MATLAB应用于实验教程，电子工业出版社，20105 李国朝，MATLAB基础及应用，北京大学出版社，20116 马莉，MATLAB数学实验与建模

22、，清华大学出版社，2010 >> hold on; for i=1:100 x=A(i,2); y=A(i,3); plot(x,y,'o')end附录三：Matlab代码1（做图）A=1 3.2060 12.9166;2 17.4571 19.3377;3 4.7576 20;B=1 50 30;2 53.7465 48.4490;3 46.9222 57.1195;5 43.1123 56.3187;4 33.3207 39.8050;C=1 54.6982 70;2 53.7465 90;3 46.9222 80;D=1 90 75; 2 80 95; 3 7

23、0 80for i=1:3 x1=A(i,2); y1=A(i,3); x2=A(mod(i,3)+1,2); y2=A(mod(i,3)+1,3); X=x1,x2; Y=y1,y2; plot(X,Y,'m')endfor i=1:5 x1=B(i,2); y1=B(i,3); x2=B(mod(i,5)+1,2); y2=B(mod(i,5)+1,3); X=x1,x2; Y=y1,y2; plot(X,Y,'m')endfor i=1:3 x1=C(i,2); y1=C(i,3); x2=C(mod(i,3)+1,2); y2=C(mod(i,3)+1,

24、3); X=x1,x2; Y=y1,y2; plot(X,Y,'m')end for i=1:3 x1=D(i,2); y1=D(i,3); x2=D(mod(i,3)+1,2); y2=D(mod(i,3)+1,3); X=x1,x2; Y=y1,y2; plot(X,Y,'m')end 附录六：Matlab代码2>> NUM=zeros(100,100);for i=1:100 for j=i+1:100 NUM(i,j)=1; endendb1=(A(1,3)-A(2,3)/(A(1,2)-A(2,2);bb1=(A(1,3)*A(2,2)-A

25、(1,2)*A(2,3)/(A(1,2)-A(2,2);b2=(A(2,3)-A(3,3)/(A(2,2)-A(3,2);bb2=(A(2,3)*A(3,2)-A(2,2)*A(3,3)/(A(2,2)-A(3,2);b3=(A(3,3)-A(1,3)/(A(3,2)-A(1,2);bb3=(A(3,3)*A(1,2)-A(3,2)*A(1,3)/(A(3,2)-A(1,2);BB=b1 bb1;b2 bb2;b3 bb3;c1=(B(1,3)-B(2,3)/(B(1,2)-B(2,2);cc1=(B(1,3)*B(2,2)-B(1,2)*B(2,3)/(B(1,2)-B(2,2);c2=(B

26、(2,3)-B(3,3)/(B(2,2)-B(3,2);cc2=(B(2,3)*B(3,2)-B(2,2)*B(3,3)/(B(2,2)-B(3,2);c3=(B(3,3)-B(4,3)/(B(3,2)-B(4,2);cc3=(B(3,3)*B(4,2)-B(3,2)*B(4,3)/(B(3,2)-B(4,2);c4=(B(4,3)-B(5,3)/(B(4,2)-B(5,2);cc4=(B(4,3)*B(5,2)-B(4,2)*B(5,3)/(B(4,2)-B(5,2);c5=(B(5,3)-B(1,3)/(B(5,2)-B(1,2);cc5=(B(5,3)*B(1,2)-B(5,2)*B(1

27、,3)/(B(5,2)-B(1,2);CC=c1 cc1;c2 cc2;c3 cc3;c4 cc4;c5 cc5;d1=(C(1,3)-C(2,3)/(C(1,2)-C(2,2);dd1=(C(1,3)*C(2,2)-C(1,2)*C(2,3)/(C(1,2)-C(2,2);d2=(C(2,3)-C(3,3)/(C(2,2)-C(3,2);dd2=(C(2,3)*C(3,2)-C(2,2)*C(3,3)/(C(2,2)-C(3,2);d3=(C(3,3)-C(1,3)/(C(3,2)-C(1,2);dd3=(C(3,3)*C(1,2)-C(3,2)*C(1,3)/(C(3,2)-C(1,2);

28、DD=d1 dd1;d2 dd2;d3 dd3;e1=(D(1,3)-D(2,3)/(D(1,2)-D(2,2);ee1=(D(1,3)*D(2,2)-D(1,2)*D(2,3)/(D(1,2)-D(2,2);e2=(D(2,3)-D(3,3)/(D(2,2)-D(3,2);ee2=(D(2,3)*D(3,2)-D(2,2)*D(3,3)/(D(2,2)-D(3,2);e3=(D(3,3)-D(1,3)/(D(3,2)-D(1,2);ee3=(D(3,3)*D(1,2)-D(3,2)*D(1,3)/(D(3,2)-D(1,2);EE=e1 ee1;e2 ee2;e3 ee3;for i=1:1

29、00 for j=i+1:100 x1=X(i,2); y1=X(i,3); x2=X(j,2); y2=X(j,3); a=(y1-y2)/(x1-x2); b=(y1*x2-x1*y2)/(x1-x2); for k=1:3 M=a -1;BB(k,1) -1; N=b 0;BB(k,2) 0; if(det(M)=0) P=MN; x=P(1,1); y=P(2,1); M1=x-A(k,2) y-A(k,3); N1=x-A(mod(k,3)+1,2) y-A(mod(k,3)+1,3); m=M1/N1; aa=x-x1,y-y1; bb=x-x2,y-y2; if(xor(aa(1

30、),bb(1)|xor(aa(1),bb(1)mm=0; else mm=aa/bb; end if(m<0) if(mm<0) NUM(i,j)=NUM(i,j)*0; else NUM(i,j)=NUM(i,j)*1; end else NUM(i,j)=NUM(i,j)*1; end end end for k=1:5 M=a -1;CC(k,1) -1; N=b 0;CC(k,2) 0; if(det(M)=0) P=MN; x=P(1,1); y=P(2,1); M1=x-C(k,2) y-A(k,3); N1=x-C(mod(k,5)+1,2) y-C(mod(k,5)

31、+1,3); m=M1/N1; aa=x-x1,y-y1; bb=x-x2,y-y2; if(xor(aa(1),bb(1)|xor(aa(1),bb(1)mm=0; else mm=aa/bb; end if(m<0) if(mm<0)NUM(i,j)=NUM(i,j)*0; else NUM(i,j)=NUM(i,j)*1; end else NUM(i,j)=NUM(i,j)*1; end end end for k=1:3 M=a -1;DD(k,1) -1; N=b 0;DD(k,2) 0; if(det(M)=0) P=MN; x=P(1,1); y=P(2,1); M

32、1=x-D(k,2) y-D(k,3); N1=x-D(mod(k,3)+1,2) y-D(mod(k,3)+1,3); m=M1/N1; aa=x-x1,y-y1; bb=x-x2,y-y2; if(xor(aa(1),bb(1)|xor(aa(1),bb(1)mm=0; else mm=aa/bb; end if(m<0) if(mm<0)NUM(i,j)=NUM(i,j)*0; else NUM(i,j)=NUM(i,j)*1; end else NUM(i,j)=NUM(i,j)*1; end end end for k=1:3 M=a -1;EE(k,1) -1; N=b

33、 0;EE(k,2) 0; if(det(M)=0) P=MN; x=P(1,1); y=P(2,1); M1=x-E(k,2) y-E(k,3); N1=x-E(mod(k,3)+1,2) y-E(mod(k,3)+1,3); m=M1/N1; aa=x-x1,y-y1; bb=x-x2,y-y2; if(xor(aa(1),bb(1)|xor(aa(1),bb(1)mm=0; else mm=aa/bb; end if(m<0) if(mm<0)NUM(i,j)=NUM(i,j)*0; else NUM(i,j)=NUM(i,j)*1; end else NUM(i,j)=NU

34、M(i,j)*1; end end end endendNUM=NUM+NUM'for i=1:100 if(SIGN(i)=0) for j=1:100 NUM(i,j)=0; NUM(j,i)=0; end endendNUM%问题三，最小生成树的计算AL=A(1:INSIGN(1)-1,:);A(INSIGN(1)+1:INSIGN(2)-1,:);. A(INSIGN(2)+1:INSIGN(3)-1,:);A(INSIGN(3)+1:INSIGN(4)-1,:);. A(INSIGN(4)+1:100,:);%AL记录不在障碍区的用户点的信息DIS=zeros(p,p);fo

35、r i=1:p for j=1:p AAA=AL(i,2)-AL(j,2),AL(i,3)-AL(j,3); DIS(i,j)=norm(AAA); endendfor i=1:p for j=1:p if(i=j) if(NUM(AL(i,1),AL(j,1)=0)DIS(i,j)=inf; end end endendDIS%记录不在障碍区各用户点之间可用不过障碍区线段连接的线段的长度%以下是kruscal算法求解最短路问题%图论最小生成树Kruskal避圈算法(使用时根据题目修改w和n)%W为邻接矩阵W=DIS;n=p;%有p个点k=1;for i=1:n-1for j=i+1:nif

36、(W(i,j)=inf) x(1,k)=W(i,j);%记录边 x(2,k)=i;%记录起点 x(3,k)=j;%记录终点 k=k+1;endendendk=k-1; %统计边数 k为边数%步骤一%冒泡法给边的大小排序for i=1:kfor j=i+1:k if x(1,i)>x(1,j) a=x(1,i);x(1,i)=x(1,j);x(1,j)=a; a=x(2,i);x(2,i)=x(2,j);x(2,j)=a; a=x(3,i);x(3,i)=x(3,j);x(3,j)=a; endendendR=zeros(1,n); %R用于记录已选的点的标号for i=1:n l(i)=

37、i;end%l 用于对点的筛选，防止形成回路%初始时选e1加入集合EEEE(1,1)=x(1,1); %EE矩阵的第一行记录最小生成树的边长EE(2,1)=x(2,1); %EE矩阵的第二行记录边的起点EE(3,1)=x(3,1); %EE矩阵的第三行记录边的终点a=min(l(x(2,1),l(x(3,1);l(x(2,1)=a;l(x(3,1)=a;%用于筛选点，防止形成回路b=1;%记录EE中边数R(1)=x(2,1);R(2)=x(3,1);c=2;for i=2:k%步骤四if b=n-1 %如果树中边数达到n-1break %算法终止end%步骤二O1=find(R=x(2,i);O2=find(R=x(3,i);size1=size(O1,2);size2=size(O2,2);if(size1=0&