自来水管道连接规划模型论文.doc

自来水管道连接规划模型自来水管道连接规划问题自来水管道连接规划模型（一）摘要：自来水是人们日常生活中不可缺少的生活要素，我们分析自来水管道连接最优问题，即在自来水管道铺设过程中在绕开障碍物的前提下的最短路径问题，使自来水管道将各个供水点用最短路径链接。根据对目标点的数据进行筛选与分析，先用面积法排除障碍区域，再对剩余点采用kruskal算法生成最优路的方案。初始给定的供水点中存在位于障碍区域中的点，需要采用合理的方法排除障碍区域中的点。本文将采用面积分析的方法，提供一种解决障碍区域判定的切实可行的方法，在二维坐标系上标定各点，障碍区域用由阴影覆盖的凸多边形表出，通过对点坐标之间的向量运算判定各点是否位于阴影区域，最终通过matlab编程实现。在确定并剔除障碍区中的点位后，采用kruskal算法生成最优路径，对于通过阴影区域的线段，将其权值设定为无穷大，最终通过编程、绘图，给出管道最优连接方案，解决本问题。最后我们对模型进行了整体评价，并提出改进之处。（二）关键词：管道连接 面积法 障碍点筛选 最短路kruskal算法 权值 最小生成树一. 问题重述自来水是人们日常生活中不可缺少的生活要素，然而自来水管网的组建却有很多问题需要解决。一般来说，我们假设管网中任意两个用户之间存在直线段相连，但是在连接过程中，有些区域是必须绕开的，这些必须绕开的区域我们称为障碍区域。表1给出了若干个可能的用户的地址的横纵坐标，可能的用户的含义是：如果用户的地址不在障碍区域内，那么该用户就是需要使用自来水的用户（即有效用户），否则如果用户的地址在障碍区域内，那么该用户就是无效用户（即不要将该用户连接在网络中）。表2-表5是分别是4个障碍区域必须要覆盖的点的坐标，而对应障碍区域就是覆盖这些要覆盖的点的最小凸集。(1) 请您判定表1中那些用户为有效用户。(2) 请设计算法筛选有效用户之间的有效线段。(3)请设计一个算法将有效用户用有效线段连接起来，并且连接的距离总和最小。表1（见附录一）表2障碍区域1必须要覆盖的点的坐标顶点序号顶点的横坐标顶点的纵坐标13.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二模型假设1， 用户之间以直线连接;2， 障碍区中的用户可以忽略，认为是无效用户；3， 以管道总距离最小为目的；4， 障碍区域是障碍顶点围成的凸多边形区域；5， 在非障碍区用户之间可确保用直线连接，且直线不通过障碍区域；三符号说明表6 论文符号说明符号表示对象a用户点的坐标b障碍区1的各顶点坐标c障碍区2的各顶点坐标d障碍区3的各顶点坐标e障碍区4的各顶点坐标sign记录各用户点是否在障碍区，若在对应位置记为1；若不在，则对应位置记为0outsign无效用户点的序号n有效用户点的个数num记录任意两用户点之间可用线段连接起来且不过障碍区的线段dis连接的长度m最小生成树的点以及连接的信息sum最小生成树管道的总长四问题分析先排除障碍区域。如果用户点位于凸边型障碍物之外，则为有效用户，否则为无效用户。再将任意两个有效用户连接，如果连接通过障碍区域之内，则为无效线段。最后通过有效点和有效连接生成最小生成树，并且连接有效用户点，画出连接路线图形，并计算生成树所成长度。根据对模型的合理假设，障碍区域即为已知若干障碍区顶点围成的凸多边形，故解决此问题的关键在于在已建立的二维坐标系中，寻找到一种合理的算法能够判定出点是否位于障碍区域中。通过直观判断，阴影区域的构成由表7给出：表7 障碍区域构成障碍区域编号构成1由3个无效用户坐标点围成的三角形2由5个无效用户坐标点围成的凸五边形3由3个无效用户坐标点围成的三角形4由3个无效用户坐标点围成的三角形运用面积法进行筛选点，对所有点进行筛选，找到并排除障碍区域中的无效用户，再把任意两个有效用户点之间用线段连接，运用向量法设计筛选线段的程序，筛选出所有不过障碍区的线段。最后设计程序，将所有有效用户点连接起来，并使管道总距离最小。这是一个典型的最小生成树问题，但相较以往最小生成树问题又有着其特别之处，以往的无障碍的情况下只需要使用弗洛伊德算法即可。但因为障碍区域的干扰，这使得坐标系并非是一个连通区域，该无法直接使用。这就需要我们在对问题进行合理假设的前提下，对已有算法进行改良。我们通过对穿过障碍区的线段赋权值为无穷大的方法，利用kruskal算法，生成最优路径。五模型的建立与求解： 5.1.问题一的模型建立和求解5.1.1运用向量的方法求解障碍区面积s若障碍区是三角形，对应各顶点坐标分别为(x1,y1)，(x2,y2), (x3,y3)。则a=(x2-x1,y2-y1),b=(x3-x1,y3-y1)。由于三角形面积s=|a|*|b|*sin/2,向量a,b外积的模长|ab|=|a|*|b|*sin;则有s=|ab|/2;若障碍区为五边形，对应点为(x1,y1)，(x2,y2), (x3,y3), (x4,y4)，(x5,y5)。则划分成三个三角形，各三角形的顶点分别为(x1,y1)，(x2,y2), (x3,y3)；(x3,y3), (x4,y4)，(x5,y5)；(x1,y1)，(x3,y3), (x5,y5)。再用求三角形面积的方法求解即可。5.1.2求用户点与任意两个同一障碍区的顶点构成三角形的面积之和s15.1.3判断有效用户 如果s=s1，则该用户在障碍区内，为无效用户。反之为有效用户。则筛选完毕的结果如下： 在障碍区的点的序号分别为：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);有效用户的个数是：96。 100个点是否在障碍区的情况如下图：5.2连接有效用户 求出过任意两个有效用户点的直线m与过各障碍区中任意两个顶点的直线l的交点坐标，再运用向量法判断该交点是否在以上述两有效用户点为端点上的线段m1和以上述障碍区顶点为端点的线段l1上，然后判断过该两个有效用户点的连接是否有效。5.2.1运用矩阵的方法求解两直线之间的交点坐标 如果任意两个有效用户点的坐标分别为a(x1,y1)、b(x2,y2)，同一障碍区任意两个顶点坐标为m(x3,y3)、n(x4,y4)。则两直线方程分别为：（1）（2）；则由解线性方程组的方法有，线性方程组的的系数矩阵为： ；在运用matlab求解该线性方程组时，不妨把 分别设为： 可以求得=a。5.2.2判断线段是否为有效线段若求得的交点坐标为p(x,y)，则通过向量关系pm=pn,可以求的。若0,则该线段为有效线段；若0,则要考虑向量关系pa=pb,若0，则该线段为有效线段，否则，该线段为无效线段。5.3问题三的模型建立与求解若要在n个用户之间连接自来水管道，由于每一个用户与其余n-1个用户之间都可能连接自来水管道。因此，在n个用户之间，最多可能连接n(n-1)/2条自来水管道然而，在连接n个用户之间的管道时，最少只需连接n-1条管道。也就是说只需要n-1条管道线路就可以把n个用户之间的自来水管道连通。现在要考虑在连接n个用户的自来水管道的同时要保证所有的管道长度之和最短，这就涉及了最优化的问题。利用kruskal算法思想设计matlab程序进行最小生成树所需边的筛选,并且设计算法将筛选出来的构成最小生成树的各边连接起来，求出最短路径长度，并画出连接图形。5.3.1利用kruskal算法思想求解最小生成树设计96个用户之间的带权图，并作出邻接矩阵dis,再根据求得的有效线段与无效线段对邻接矩阵进行修改，将邻接矩阵中对应无效线段的位置的值修改为inf，可以得到一个新的邻接矩阵dis。接下来，用冒泡排序法对所有有效线段长度按从小到大的顺序进行排序。这时，需要借助kruskal算法进行最小生成树的计算。然后把最小生成树对应边的线段长度、起点、终点信息记录在矩阵ee中。生成最小生成树时，从长度最短的边开始选取。首先不妨设一个196的标记向量l用于记录被选取的点的序号，初始状态向量l的各元素依次为各用户序号，在选取线段为边后，将对应两点的序号m与n取最小值，并将向量l中所有与m位置元素相等的元素位置及所有与n位置的元素相等的元素位置都赋值为该最小值，如此循环知道向量l中所有元素均相等时停止；同时可以设一向量r来依次记录被选点的序号，直到所有用户点被无重复地被记录。在按线段长度从小到大的顺序选择边时，设线段端点用户的序号为m与n。这时需要考虑如下4种情况：如果在向量r中m和n均没有被记录,则该线段可以被选为最小生成树的边，将对应线段的信息记录在矩阵ee中，同时在r中添加记录m和n的值，并按照上述步骤更新向量l。如果在向量r中m被记录而n没有被记录，则该线段可以被选为最小生成树的边，将对应线段的信息记录在矩阵m中，同时在r中添加记录n的值，并按照上述步骤更新向量l。如果在向量r中n被记录而m没有被记录，则该线段可以被选为最小生成树的边，将对应线段的信息记录在矩阵ee中，同时在r中添加记录m的值，并按照上述步骤更新向量l。如果在向量r中m和n均被记录,则需要借助向量l来判断是否该线段可以被选为最小生成树的边：a. 如果向量l中对应的m位置与n位置的元素值相等，则该线段不是最小生成树的边，直接跳过到下一步判断。b. 如果向量l中对应的m位置与n位置的元素值不相等，则该线段是最小生成树的边，将对应线段的信息记录在矩阵m中，同时只需要更新向量l。通过上述方法，即可产生最小生成树，其各边信息记录在矩阵m中。 5.3.2设计matlab程序求出最小生成树长度并将各边连接起来要计算最小生成树的长度，只需要借助for循环将m矩阵中记录长度相加即可。具体算发如下： sum=0;for i=1:p-1 sum=sum+m(1,i);endsum可以求得最小生成树的长度为：sum= 653.0196；最后，借助plot函数画出最小生成树的图形。算法如下： hold on; for i=1:100 x=a(i,2); y=a(i,3); plot(x,y,o)end for i=1:n-1 x1=al(m(2,i),2); y1=al(m(2,i),3); x2=al(m(3,i),2); y2=al(m(3,i),3); x=x1,x2; y=y1,y2; plot(x,y)endfor i=1:3 x1=b(i,2); y1=b(i,3); x2=b(mod(i,3)+1,2); y2=b(mod(i,3)+1,3); x=x1,x2; y=y1,y2; plot(x,y,m)endfor i=1:5 x1=c(i,2); y1=c(i,3); x2=c(mod(i,5)+1,2); y2=c(mod(i,5)+1,3); x=x1,x2; y=y1,y2; plot(x,y,m)endfor 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 for i=1:3 x1=e(i,2); y1=e(i,3); x2=e(mod(i,3)+1,2); y2=e(mod(i,3)+1,3); x=x1,x2; y=y1,y2; plot(x,y,m)end连接形成的最小生成树的图形如下图所示：由图可知，该最小生成数是合理的。 六模型检验 首先可以通过对所画最小生成树图形的观察，看是否有回路，由图易知图形中无回路，则通过修改最小生成树中任意边的连接，计算修改后的最小生成树的长度sum与sum进行比较。可得sumsum,则该模型所生成的最小生成树的长度最短，即运用该模型进行自来水管道的连接所需要的自来水管长度最短。十附录附录：若干个可能的用户的地址的横纵坐标可能的用户的序号可能的用户横坐标可能的用户纵坐标1.000095.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.109516.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.881449.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.878545.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.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.182874.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.000025.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附录三：解决该自来水管道连接问题的matlab程序：%问题一，进行有效点的筛选a=1.0000 95.0129 58.2792;2.0000 23.1139 42.3496;3.0000 60.6843 51.5512;4.0000 48.5982 33.3951;5.0000 89.1299 43.2907;6.0000 76.2097 22.5950;7.0000 45.6468 57.9807;8.0000 1.8504 76.0365;9.0000 82.1407 52.9823;10.0000 44.4703 64.0526;11.0000 61.5432 20.9069;12.0000 79.1937 37.9818;13.0000 92.1813 78.3329;14.0000 73.8207 68.0846;15.0000 17.6266 46.1095;16.0000 40.5706 56.7829;17.0000 93.5470 79.4211;18.0000 91.6904 5.9183;19.0000 41.0270 60.2869;20.0000 89.3650 5.0269;21.0000 5.7891 41.5375;22.0000 35.2868 30.4999;23.0000 81.3166 87.4367;24.0000 0.9861 1.5009;25.0000 13.8891 76.7950;26.0000 20.2765 97.0845;27.0000 19.8722 99.0083;28.0000 60.3792 78.8862;29.0000 27.2188 43.8659;30.0000 19.8814 49.8311;31.0000 1.5274 21.3963;32.0000 74.6786 64.3492;33.0000 44.5096 32.0036;34.0000 93.1815 96.0099;35.0000 46.5994 72.6632;36.0000 41.8649 41.1953;37.0000 84.6221 74.4566;38.0000 52.5152 26.7947;39.0000 20.2647 43.9924;40.0000 67.2137 93.3380;41.0000 83.8118 68.3332;42.0000 1.9640 21.2560;43.0000 68.1277 83.9238;44.0000 37.9481 62.8785;45.0000 83.1796 13.3773;46.0000 50.2813 20.7133;47.0000 70.9471 60.7199;48.0000 42.8892 62.9888;49.0000 30.4617 37.0477;50.0000 18.9654 57.5148;51.0000 19.3431 45.1425;52.0000 68.2223 4.3895;53.0000 30.2764 2.7185;54.0000 54.1674 31.2685;55.0000 15.0873 1.2863;56.0000 69.7898 38.3967;57.0000 37.8373 68.3116;58.0000 86.0012 9.2842;59.0000 85.3655 3.5338;60.0000 59.3563 61.2395;61.0000 49.6552 60.8540;62.0000 89.9769 1.5760;63.0000 82.1629 1.6355;64.0000 64.4910 19.0075;65.0000 81.7974 58.6918;66.0000 66.0228 5.7581;67.0000 34.1971 36.7568;68.0000 28.9726 63.1451;69.0000 34.1194 71.7634;70.0000 53.4079 69.2669;71.0000 72.7113 8.4079;72.0000 30.9290 45.4355;73.0000 83.8496 44.1828;74.0000 56.8072 35.3250;75.0000 37.0414 15.3606;76.0000 70.2740 67.5645;77.0000 54.6571 69.9213;78.0000 44.4880 72.7509;79.0000 69.4567 47.8384;80.0000 62.1310 55.4842;81.0000 79.4821 12.1047;82.0000 95.6843 45.0754;83.0000 52.2590 71.5883;84.0000 88.0142 89.2842;85.0000 17.2956 27.3102;86.0000 97.9747 25.4769;87.0000 27.1447 86.5603;88.0000 25.2329 23.2350;89.0000 87.5742 80.4872;90.0000 73.7306 90.8398;91.0000 13.6519 23.1894;92.0000 1.1757 23.9313;93.0000 89.3898 4.9754;94.0000 19.9138 7.8384;95.0000 29.8723 64.0815;96.0000 66.1443 19.0887;97.0000 28.4409 84.3869;98.0000 46.9224 17.3900;99.0000 6.4781 17.0793;100.0000 98.8335 99.4295;b=1 3.2060 12.9166;2 17.4571 19.3377;3 4.7576 20;c=1 50 30;2 53.7465 48.4490;3 46.9222 57.1195;5 43.1123 56.3187;4 33.3207 39.8050;d=1 54.6982 70;2 53.7465 90;3 46.9222 80;e=1 90 75; 2 80 95; 3 70 80;%准备用户坐标矩阵和障碍区坐标矩阵sign=zeros(1,100);%生成1行100列的记录矩阵，若点不在障碍区则对应位置记为1，否则记为0for n=1:100a1=b(2,2)-b(1,2),b(2,3)-b(1,3),0;a2=b(3,2)-b(2,2),b(3,3)-b(2,3),0;a3=b(1,2)-b(3,2),b(1,3)-b(3,3),0;s=norm(cross(a1,a2)/2;%该障碍区面积x=a(n,2);y=a(n,3);z1=x-b(1,2),y-b(1,3),0;z2=x-b(2,2),y-b(2,3),0;z3=x-b(3,2),y-b(3,3),0;s1=norm(cross(z1, a1)/2;s2=norm(cross(z2, a2)/2;s3=norm(cross(z3, a3)/2;s1=s1+s2+s3;%由用户点和任意对应两个该障碍区顶点构成的三角形面积之和if(s1=s)m1=0 ;else m1=1;end%第一个障碍区内的点判断b1=c(2,2)-c(1,2),c(2,3)-c(1,3),0;b2=c(3,2)-c(2,2),c(3,3)-c(2,3),0;b3=c(4,2)-c(3,2),c(4,3)-c(3,3),0;b4=c(5,2)-c(4,2),c(5,3)-c(4,3),0;b5=c(1,2)-c(5,2),c(1,3)-c(5,3),0;s1=norm(cross(b1,b2)/2;s2=norm(cross(b3,b4)/2;s3=norm(cross(b1+b2,b3+b4)/2;s=s1+s2+s3;%该障碍区面积x=a(n,2);y=a(n,3);z1=x-c(1,2),y-c(1,3),0;z2=x-c(2,2),y-c(2,3),0;z3=x-c(3,2),y-c(3,3),0;z4=x-c(4,2),y-c(4,3),0;z5=x-c(5,2),y-c(5,3),0;s1=norm(cross(z1,b1)/2;s2=norm(cross(z2,b2)/2;s3=norm(cross(z3,b3)/2;s4=norm(cross(z4,b4)/2;s5=norm(cross(z5,b5)/2;s1=s1+s2+s3+s4+s5;%由用户点和任意对应两个该障碍区顶点构成的三角形面积之和if(s1=s)m2=0;else m2=1;end%第二个障碍区内的点判断c1=d(2,2)-d(1,2),d(2,3)-d(1,3),0;c2=d(3,2)-d(2,2),d(3,3)-d(2,3),0;c3=d(1,2)-d(3,2),d(1,3)-d(3,3),0;s=norm(cross(c1,c2)/2;%该障碍区面积x=a(n,2);y=a(n,3);z1=x-d(1,2),y-d(1,3),0;z2=x-d(2,2),y-d(2,3),0;z3=x-d(3,2),y-d(3,3),0;s1=norm(cross(c1,z1)/2;s2=norm(cross(c2,z2)/2;s3=norm(cross(c3,z3)/2;s1=s1+s2+s3;%由用户点和任意对应两个该障碍区顶点构成的三角形面积之和if(s1=s)m3=0;else m3=1;end%第三个障碍区内的点判断d1=e(2,2)-e(1,2),e(2,3)-e(1,3),0;d2=e(3,2)-e(2,2),e(3,3)-e(2,3),0;d3=e(1,2)-e(3,2),e(1,3)-e(3,3),0;s=norm(cross(d1,d2)/2;%该障碍区面积x=a(n,2);y=a(n,3);z1=x-e(1,2),y-e(1,3),0;z2=x-e(2,2),y-e(2,3),0;z3=x-e(3,2),y-e(3,3),0;s1=norm(cross(d1,z1)/2;s2=norm(cross(d2,z2)/2;s3=norm(cross(d3,z3)/2;s1=s1+s2+s3;%由用户点和任意对应两个该障碍区顶点构成的三角形面积之和if(s1=s)m4=0 ;else m4=1;end%第四个障碍区内的点判断if(m1&m2&m3&m4)sign(n)=1;else sign(n)=0;endendsign ;%以上判断不在障碍区的点，其标号记录在sign矩阵中用1表示outsign=find(sign=0)%记录在障碍区的用户点p=0;for i=1:100 p=p+sign(i);endp%计算保留用户点的个数p%问题二，进行有效线段的筛选num=zeros(100,100);for i=1:100 for j=i+1:100 num(i,j)=1; endendb1=(b(1,3)-b(2,3)/(b(1,2)-b(2,2);bb1=(b(1,3)*b(2,2)-b(1,2)*b(2,3)/(b(1,2)-b(2,2);b2=(b(2,3)-b(3,3)/(b(2,2)-b(3,2);bb2=(b(2,3)*b(3,2)-b(2,2)*b(3,3)/(b(2,2)-b(3,2);b3=(b(3,3)-b(1,3)/(b(3,2)-b(1,2);bb3=(b(3,3)*b(1,2)-b(3,2)*b(1,3)/(b(3,2)-b(1,2);bb=b1 bb1;b2 bb2;b3 bb3;c1=(c(1,3)-c(2,3)/(c(1,2)-c(2,2);cc1=(c(1,3)*c(2,2)-c(1,2)*c(2,3)/(c(1,2)-c(2,2);c2=(c(2,3)-c(3,3)/(c(2,2)-c(3,2);cc2=(c(2,3)*c(3,2)-c(2,2)*c(3,3)/(c(2,2)-c(3,2);c3=(c(3,3)-c(4,3)/(c(3,2)-c(4,2);cc3=(c(3,3)*c(4,2)-c(3,2)*c(4,3)/(c(3,2)-c(4,2);c4=(c(4,3)-c(5,3)/(c(4,2)-c(5,2);cc4=(c(4,3)*c(5,2)-c(4,2)*c(5,3)/(c(4,2)-c(5,2);c5=(c(5,3)-c(1,3)/(c(5,2)-c(1,2);cc5=(c(5,3)*c(1,2)-c(5,2)*c(1,3)/(c(5,2)-c(1,2);cc=c1 cc1;c2 cc2;c3 cc3;c4 cc4;c5 cc5;d1=(d(1,3)-d(2,3)/(d(1,2)-d(2,2);dd1=(d(1,3)*d(2,2)-d(1,2)*d(2,3)/(d(1,2)-d(2,2);d2=(d(2,3)-d(3,3)/(d(2,2)-d(3,2);dd2=(d(2,3)*d(3,2)-d(2,2)*d(3,3)/(d(2,2)-d(3,2);d3=(d(3,3)-d(1,3)/(d(3,2)-d(1,2);dd3=(d(3,3)*d(1,2)-d(3,2)*d(1,3)/(d(3,2)-d(1,2);dd=d1 dd1;d2 dd2;d3 dd3;e1=(e(1,3)-e(2,3)/(e(1,2)-e(2,2);ee1=(e(1,3)*e(2,2)-e(1,2)*e(2,3)/(e(1,2)-e(2,2);e2=(e(2,3)-e(3,3)/(e(2,2)-e(3,2);ee2=(e(2,3)*e(3,2)-e(2,2)*e(3,3)/(e(2,2)-e(3,2);e3=(e(3,3)-e(1,3)/(e(3,2)-e(1,2);ee3=(e(3,3)*e(1,2)-e(3,2)*e(1,3)/(e(3,2)-e(1,2);ee=e1 ee1;e2 ee2;e3 ee3;for i=1:100 for j=i+1:100 x1=a(i,2); y1=a(i,3); x2=a(j,2); y2=a(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-b(k,2) y-b(k,3); n1=x-b(mod(k,3)+1,2) y-b(mod(k,3)+1,3); m=m1/n1; aa=x-x1,y-y1; bb=x-x2,y-y2; if(