数学建模论文-基于贪婪算法的公园内道路设计模型 37页

编程编程建模建模组别：7建模组员：公园道路设计问题组别：___________建模37 目录基于贪婪算法的公园内道路设计模型1摘要11问题重述22问题分析22.1问题1的分析32.2问题2的分析32.3问题3的分析33模型假设34符号说明35模型的建立与求解45.1问题1模型建立与求解45.1.1问题1模型的建立45.1.2问题1模型的求解65.2问题2模型的建立与求解85.2.1问题2模型的建立85.2.2问题2模型的求解95.2.3问题2模型的结果对比115.3问题3模型的建立与求解145.3.1问题3模型的建立145.3.2问题3模型的求解145.3.3问题3模型结果分析186模型评价196.1模型的优点206.2模型的缺点207模型推广208参考文献219附录2137 基于贪婪算法的公园内道路设计模型摘要本题通过公园内部规划设计道路，在已有的边界条件下，找出相应的最短路径最优解，可利用贪婪算法，通过局部优化逐步逼近最优解。对问题1，以给定的四个交叉点的情况下，寻找公园内的最短道路最优解，并满足约束条件。根据贪婪算法的基本原理，可以先用经典算法prim或kruskal对组成的点集构造最小生成树，再用floyd算法逐个筛选在最小生成树下的解，通过边和点集的不断更换，求得任意两点间的最短路径矩阵，进而求解最短道路的最优集。对问题2，考虑到交叉点的位置选取与交叉点数目对问题的双重影响，通过列举0,1,2,3四个情况下可能存在的交叉点进行建模。特别考虑到，任意两点间的最短路径要满足小于1.4两点间直线距离，可以利用椭圆的相关性质，在矩形区域相做交叉点的交点点集，简化问题2对交点位置的穷举过程。利用局部优化，在矩形区域内分别建立H型交叉点模型和双Y型交叉点模型，再通过交叉点数目的变化，求解得两种模型下的最优解，并进行对比，分析此两种模型下的最优程度，做出评判标准和相应交叉点坐标。对问题3，利用公园内增加矩形湖这一约束条件，对问题2中的不用交叉点模型进行验证，通过合理假设湖边道路存在并不计入总道路长，简化模型的操作，通过交替迭代法，建立在H型交叉点模型和双Y型交叉点模型基础上的局部最优解。再根据穷举法，对矩形区域内所有的点进行最优求解。并通过交替迭代法与穷举法的多次使用，逐步逼近全局最优解。对3个问题的综合分析后，贪婪算法下的最短路径最优解还可以在城市规划中，利用不同交叉点模型的局部优化程度的不同性，按功能和价值等对城市进行合理规划。关键词：贪婪算法　局部最优解　最小生成树　Floyd　prim　逐次逼近37 1问题重述现计划建一个形状为矩形公园，公园计划有若干个入口，需要建立一个模型去设计道路让任意两个入口相连，使总的道路长度和最小（公园四周不计入总长），前提要求是任意的两个入口之间的最短道路长不大于两点连线的1.4倍。设计对象可为如图1所示，长200米，宽100米的矩形公园。要求公园内新修的道路与四周的连接只能与8个路口相通，而不能连到四周的其它点，现要完成以下问题：1、假定公园内确定要使用4个道路交叉点为：A(50,75)，B(40,40)，C(120,40)，D(115,70）。通过建立模型并给出算法使公园内道路的总路程最短，画出道路设计，计算新修路的总路程。2、现在公园内可以任意修建道路，如何在满足条件下使总路程最少。建立模型并给出算法。给出道路交叉点的坐标，画出道路设计，计算新修路的总路程。3、若公园内有一条矩形的湖，新修的道路不能通过，但可以到达湖四周的边，示意图见图2，重复完成问题2 的任务。图1　公园入口坐标2问题分析本题是一个道路设计的最优化的问题，即是如何利用已有的边界，使公园内部新修路总长最小，同时满足以下两个约束条件：K1:任两个入口连通；K2:任两个入口的最短路径不超过其直线距离的1.4倍。2.1问题1的分析找出四个交叉点下的最小路问题，可以根据prim37 算法，可以求得公园内四个交叉点下的最小生成树，在此条件下利用Floyd算法计算出任意两点之间的最短路径，做出最短路径矩阵并进行校对。2.2问题2的分析在问题1的基础上，解决问题2的关键在于交叉点的数目与位置的分析与求解，利用本题中对最小路径的1.4倍约束，容易发现若以其中任意两点为焦点利用椭圆的几何性质就可以求得交叉点大致范围，同时在交叉点数目m为{1,2,3….}时，再次利用问题1求解此条件下的最小路径问题。2.3问题3的分析利用问题2中的不同交叉点模型下的最短路径最优解，主要分析此条件下的路径与湖的交叉问题，通过合理的假设令湖边的路长不计入总路长，利用穷举法，对最优路径逐步逼近得到最优解。3模型假设1、忽略道路宽度，交叉点的不影响道路长度；2、问题3中的沿湖道路长度不计入总道路长度；3、所有的道路均为直线，门与交叉点均为点参与计算；4、任意两点间最短路径不构成环路。4符号说明在模型构建过程中，为简化说明，将主要参与模型中计算的量用符号表示，如表1所示：表1　运算符号说明符号说明PiAimdijDDssijSxijX公园入口数（i=1,2,3….8）公园内的点数（i=1,2,3…n）公园内的交叉点数（m=n-8）从i点到j点的距离dij所构成的距离矩阵D矩阵的控制矩阵Ds=1.4D从i点到j点的最短路径距离从i点到j点的最短距离矩阵0~1决策变量xij所构成的决策变量矩阵37 根据上述说明，符号间还满足下列关系：其中，，同时其中xij满足0-1决策变量5模型的建立与求解5.1问题1模型建立与求解5.1.1问题1模型的建立问题一中已假定4个交叉点，故可构造公园内点集为：为尽可能地简化模型，则应使通过公园边界的道路长度尽可能地长，通过matlab计算任意两点间的边界距离和直线1.4倍距离（仅在1~8个点），并制成表格(仅制成下三角，下同)（见附录1）：表2　门口点与点间的边界距离P1P2P3P4P5P6P7P8P103014023024015513045P2011020027018516075P309022029527018537 P40130215240275P5085110195P6025110P7085P80利用matlab对Ds求解，可以得到表3的八个入口的1.4倍直线距离。表3　八个入口的1.4倍直线距离P1P2P3P4P5P6P7P8P1042196261.5197.9141.5140.644.8P20154221.3170.8141.5150.778.2P3089.6150.7224.1252.3226.7P40132.0241.3275.0282.1P50119154198.1P6035115.8P70105.9P80将上表与边界上的点间距和1.4倍直线距离进行对比可以发现，15，16，18,25，26,27,34,35,36,37这10个距离不满足小于两点间1.4倍直线距离的要求，则在以后的讨论中将视为重点讨论。考虑此问题的高计算度，通过对每个点和边权值在0-1决策变量下的所有情况，可以得到最短路方程：目标函数：约束条件：s.t.D-Ds<0在此基础上，A分别表示8个入口与4个交叉点，边表示两点之间的连线，将表1中的直线距离赋予相应的边的权值。并引用贪婪算法来求解问题1中的最小路径问题，下面简要介绍贪婪算法：37 根据贪婪算法的基本求解过程，问题划分为如下3步：1）忽略网中的环路存在，利用prim算法构造最小生成树（见附录2）；2）通过Floyd算法对任意两点间的最短路径求解，并构造相应的最短路径矩阵，同时检验最短矩阵是否满足控制条件K1，若不满足，继续搜索，若满足，执行下一步（见附录3）；3）进行边的删除和替换（允许环路的存在），使得任意点间最短路径满足不超过直线距离1.4倍的要求（即K2）。最后利用穷举法逐步逼近得到最优解。5.1.2问题1模型的求解将问题1中的交叉点与入口点构成的点集，按照贪婪算法的3个步骤进行运算，其中12个点集坐标如下：1）利用matlab，根据prim算法容易求解得到决策变量矩阵X中各个Xij对应的值，其中，0代表i点和j点之间没有连线，1表示i点和j点之间有连线，（0-1）矩阵X结果如下：利用0-1矩阵X可以得到在最大利用边界情况下的各点连线图：37 图3　最小生成树的最短路径图2）根据决策变量X所代表的连线情况，假定任意两点间最小路径不存在环路(即不考虑无连接的边界)，利用floyd算法求得最小生成树的最短路径矩阵S：利用表3对比，利用matlab做出最短路径矩阵S与控制条件矩阵Ds的0~1差值表，差值表中若只含有0的项（对角线上的0除外），表示Ds>S，则所求的最短路矩阵满足要求，是此情况下的最优解，否之则否，其中差值表如下：表4　最短路矩阵S与控制条件矩阵Ds的0~1差值表P1P2P3P4P5P6P7P8P101000101P20000111P301111037 P400010P50010P6010P700P80从上表中可以发现，0~1差值表中含有１的项有16项，其中有14项可以通过边界直接连接，这是因为模型建立在不含无连接边界的基础上，故将边界加入考虑可以得到上表中，只有15,25不满足条件。3）考虑公园边界可以节约道路总长，则可以利用穷举法，对不满足的点进行替换，容易知道此情况下的替换次数是有限的，从而接近满足K1条件的最优解。从图1中可以得到八条可替换路线：59,510,52,129,1210,122,1110,112，再次求解各个情况下的最短路径距离S，并判断是否满足条件K2从而，筛选出最优解。其中最优解如下图：图4　4个交叉点的最优道路设计图在此条件下，利用目标函数求解得，公园内道路总长最优解为：5.2问题2模型的建立与求解5.2.1问题2模型的建立37 在问题1模型的基础上，问题2允许公园内可修建任意道路，为简化模型，寻求最优解，通过控制交叉点数m（0,1,2,3,…），找到局部最优解。再利用交替迭代法，在matlab条件下，通过控制局部最优解，迭代接近全局最优解。5.2.2问题2模型的求解<1>　m=0根据表2，在公园内交叉点数为0时，有10个点间不能满足控制条件K2。所以，必须在公园内添加交叉点。<2>　m=1从题目的控制条件K2可以发现，若公园内只有一个交叉点，则任意两两选择入口，则分别以两个入口点做为椭圆的焦点，以两点间1.4倍焦距距离做半长轴画椭圆，则这两个椭圆必有一个交点（见附录4、5、6、7）。根据上述原理，分别选取P1、P7和P3、P5为两个椭圆的焦点，利用matlab可以绘制如下的椭圆交点图：图5　m=1椭圆交点图分析上图，发现两个椭圆并无交点，故公园内的交叉点数。<3>　m=2利用上述方法，对8个入口点分别做椭圆，为求解最短路径解，则椭圆须在矩形内部，将不在矩阵内部或是与一个边有两个交点的椭圆删去。由此可知，椭圆包围的最密集的地区出现交点的可能性最大，可以利用交替迭代法，在这个范围内逐步逼近最优解。椭圆交点图如下：37 图6　m=2椭圆交点图（a）H型算法：为简化模型，可预先设定P1与P8、P3与P4间存在直接连线。根据图6所示的大致区域分别在两个区域内选定两个点，记为Firststep：将两个点的搜索半径定为4，并选定两个区域内的一个角点，利用matlab的计算能力，穷举两个区域内的任意点，利用问题1中的floyd算法计算最短路径，与控制条件相比对，若满足则进行记录。Secondstep：在第一步的基础上，将搜索半径定为0.1，重新运行步骤一的过程，并再新结果重新记录。Thirdstep：得到在搜索半径4和搜索半径0.1条件下的结果进行比对，根据结果，考虑是否进行新一轮的穷举（见附录8）。（b）双Y型算法：在H型算法中，对两条道路进行设定，考虑到设定的道路对最优解的影响，通过减少设定路线从而逼近最优解（此外仅减少P3与P4间的连线）。根据H型算法中的基本思想，对双Y型算法中存在的最优解的点，同样利用穷举法在两个椭圆包含密集区域搜索（见附录9、10）。<4>　m=3考察<3>中m=2的两种情况，两个交叉点的相对位置偏离十分明显，考虑椭圆包含的密集程度，在双Y型算法的基础上，增设交叉点。令交叉点的坐标为并让F点与P5、H、I三个点相连接，通过固定H、I逐步逼近F的最优解，再通过固定P5、H逐步逼近I的最优解，以此类推，穷举出区域范围内满足条件的最短道路的最优解（见附录11、12）。37 5.2.3问题2模型的结果对比通过问题2对交叉点的数目与位置对最短路径距离的分析，可以得到在2个交叉点下，H型交叉点与双Y型交叉点的最优解，如下图7、图8：图7　H型交叉点的最优解图8　双Y型交叉点的最优解37 利用最短路径的目标函数，可以计算得到，H型算法下的两个交叉点的最优解为：H(66,55)　I(107,63)相应的最短路径为：其中，制作Ds与Sij的差值表，可得如表5所示，H型交叉点满足控制条件K2。表5　H型交叉点下的Ds与Sij的差值表P1P2P3P4P5P6P7P8P10.012.056.057.529.70.423.712.8P20.044.047.332.629.613.916.2P30.025.629.245.448.654.7P40.02.126.435.146.1P50.034.044.03.1P60.010.05.9P70.020.9P80.0同理，可以取得双Y型算法下的两个交叉点的最优解为：H(59.7,77.6)　I(173.1,43.6)相应的最短路径为：在交叉点数目为3的情况下，利用穷举法绘制出最优解，如图9所示：37 图9　m=3时的最短道路最优解同时，获得双Y型交叉点中Ds与Sij的差值表，如表6所示，双Y型交叉点满足控制条件K2。表6　双Y型交叉点的Ds与Sij的差值表P1P2P3P4P5P6P7P8P10.012.056.047.326.70.323.712.8P20.044.037.129.630.314.516.2P30.015.428.817.120.454.7P40.027.051.360.035.9P50.034.044.03.1P60.010.05.9P70.020.9P80.0可以得到三个交叉点的最优解为：H(62,70)　I(169,40)　F(117,87)相应的最短道路长度为：37 通过上述对比发现，在交叉点数目ｍ=2时，双Y型交叉点相较H型交叉点最短道路长度更短，故选取双Y型交叉点为交点数目为2时的最优解，而在交叉点数选为3时，与交叉点为2时的双Y型交叉点对应的最短道路长度之差，相距甚微，从而类推不必再求解四个交叉点的最短道路长度。最优解可以估算为m=2时的最短道路长度。即S=358.52985.3问题3模型的建立与求解5.3.1问题3模型的建立该问题在是建立在问题2的基础上加入了一个矩形湖（如图2所示）,湖的四个点坐标分别R1(140,70),R2(140,45),R3(165,45),R4(165,70)，并且湖边的道路长度不计入总长度，新建路径不能直接穿过。若原来道路通过湖则考虑如下两种情况：　　情况一：对于问题2中所求得的交叉点向矩形湖引垂线，并连接最近入口形成新的路径；情况二：重新确定新的范围，利用问题2的方法确定新的交叉点和路径；最后比较一、二两种情况，得出最优解。5.3.2问题3模型的求解（1）在问题2模型的基础上，首先利用双Y型算法中交叉点的位置的连线可以得到交叉点的位置图，如图10所示（见附录13）：37 图10　双Y型交叉点位置图情况一：根据模型的建立要求，可将模型分为两步讨论。Step1：假设I点移入湖内，分别对湖两边引垂线，易知此情况下的相应点的位置，向入口点引直线。Step2：在step1的基础上，使I点在湖内移动，利用穷举法逼近此条件下的最优解。根据上述两个步骤，可以在matlab中绘制得到此条件下的最短道路总长的最优解，如图11所示：37 图11　I点在湖内移动的最优解情况二：使Ｉ点不在湖内移动，分别在湖的领域内移动，同样利用情况一的解法，求得Ｉ点的最短道路长度最优解，如图12所示（见附录14、15）：图12　I点在湖边领域内的最短道路最优解（2）在（1）的基础上，考虑问题2中的3个交叉点在有湖情况下的最优解。根据（1）中的步骤，考虑3个交叉点间可以相互约束，共同逼近最优解的情况，分别固定H和I或湖边3个接点37 分别分析求解最优解，逼近情况可以用图13、图14、图15、图16表示：图13　在H型交叉点情况下对湖边最优交点解图14　固定湖边3个接点下对H和I的最优解37 图15　再次对湖边3个点的最优搜索图16　再次固定湖边3点对H和I的最优搜索5.3.3问题3模型结果分析从双Y型交叉点的最优解图11、图12，并计算相应的最短道路总长，可以推出图12为双Y型交叉点在有湖情况下的最优解，其相应的湖边37 接点和最优交叉点分别为：R1(140,70)，R3(165,45)，O（169.5,40.7）其最短道路总长为：S=323.716而在图13至图16中，利用matlab算得，在H型交叉点情况下，四步优化的最短道路总长分别为：S1=338.2392S2=330.0933S3=325.2299S4=322.6000可以推出在H型交叉点条件下，图16为最短道路总长最优解，相应的最优湖边接点分别为：A(140,61)，B(161,45)，C(165,50)利用matlab可以绘得，图16情况下的Ds与Sij的差值表，满足控制条件K2，如表7所示。表7　H型交叉点逼近最优解的Ds与Sij的差值表P1P2P3P4P5P6P7P8P10.012.016.031.526.00.310.712.8P20.044.021.428.929.713.93.3P30.00.662.851.254.441.7P40.02.126.435.157.2P50.034.024.03.1P60.010.05.9P70.020.9P80.0通过结果对比可以发现，双Y型交叉点并未和问题2中比H型交叉点的最优解更优，这主要与本题的假设，湖边道路不计入总道路长度有关，而两者的最后最优值相差并不是十分显著，故在公园内有湖的情况下，公园内的道路总长的最优解为：S=322.600006模型评价本题模型建立在任意两点间相连接，并在约束条件任意两点间最短路径小于两点间直线距离的1.4倍，公园内道路的最短路径问题。模型建立于贪婪算法基础上，利用编程的穷举法对区域内所满足的点进行一一验证。37 6.1模型的优点1、与一般的优化相比，本模型通过局部交替优化的方式，大大简化了优化过程。2、充分利用公园周围边界，使公园内道路最短。3、对交叉点在公园内的交叉模型和交叉点数目和列举，可以大概推算模型的最优解，简化了操作过程。6.2模型的缺点1、在matlab上的算法是对区域内的所有点进行穷举，加大了计算量和劳动量。2、在局部所求得的最优解，难以逼近或证明为全局最优解。7模型推广1、对比表5、表6容易发现，不同的交叉点模型对各个点间的最短道路长度的优化程度不同，如H型交叉点中的P4至P2、P3、P5、P7和P3至P7与双Y型对应的差值有较大的不同，但是可以得到相应的最优解，则在城市规划模型中，可以按要求对各区域职能的重要性进行分区。2、在问题3中，对模型做了一个假设，即湖边道路已建，并不记入总道路长度，故在要求湖边道路长度归入总长度时，对双Y型交叉点做再一次的穷举逼近，可以得到此条件下的最优解，如图17所示：图17　记入湖边道路的双Y型最优解37 8参考文献[1]姜启源,谢金星,叶俊.数学模型[M]（第三版）.北京:高等教育出版社,2003.8.[2]赵静,但琦.数学建模与数学实验（第二版）.北京:高等教育出版,2003.8.[3]薛翠平,张薇.最短路问题的灵敏度分析与最短路调整.燕山大学学报.2009年1月第33卷第一期.[4]黄辉,梁国宏,张生,何尚录.求解一类线性规划问题的原始贪婪算法和对偶贪婪算法及其相互关系.甘肃:兰州大学学报2007年2月第26卷第1期.[5]王英,刘天时.基于Kruskal算法的最短路径算法研究.重庆文理学院学报2009年12月第28卷第6期.[6]王志栋.城市道路网络设计的数学模型及其解法.大连:大连铁道学院学报,1997年3月第18卷第1期.[7]胡华,石世光.Floyd算法分析与演示系统设计.电脑学习2007年12月第6期.[8]伍少坤,黎夏,刘小平.基于城市扩张的动态选址模型.地理科学2008年6月第28卷第3期.[9]朱友勇.基于图像分割法的研究与应用.江南大学硕士学位论文.[10]陈丹娜,田毅.浅谈城市道路设计中的人性化设计.哈尔滨阿城区规划勘察设计院.9附录附录1两点间直线距离%array.m计算12个点间的直线距离array1=[200;500;1600;20050;120100;35100;10100;025;5075;4040;120,40;11570];%12个点的坐标，第一列表示横坐标，第二列表示纵坐标array2=zeros(12);%存放两点间距离的矩阵fori=1:12forj=i:12array2(i,j)=sqrt((array1(i,1)-array1(j,1))^2+(array1(i,2)-array1(j,2))^2);%计算距离array2(j,i)=array2(i,j);%矩阵式对阵矩阵endend附录2%Prim算法求最小生成树result=[];%存放最后结果的矩阵：第一行和第二行存放相邻边的顶点，第三行存放点间的距离37 p=1;tb=2:length(array2);whilelength(result)~=length(array2)-1temp=array2(p,tb);temp=temp(:);d=min(temp);[jb,kb]=find(array2(p,tb)==d);j=p(jb(1));k=tb(kb(1));result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[];end附录3%Folyd求最短路径n=12;dis=zeros(n);xij=zeros(n);forj=1:11dis(result(1,j),result(2,j))=result(3,j);dis(result(2,j),result(1,j))=dis(result(1,j),result(2,j));Xij(result(1,j),result(2,j))=1;Xij(result(2,j),result(1,j))=1;endfori=1:nforj=1:nifdis(i,j)==0dis(i,j)=inf;endifi==jdis(i,j)=0;endendendpath=zeros(n);fork=1:nfori=1:nforj=1:nifdis(i,j)>dis(i,k)+dis(k,j)dis(i,j)=dis(i,k)+dis(k,j);path(i,j)=k;endendendenddif=array2(1:8,1:8)*1.4;%点间距离的1.4倍37 dif>=dis(1:8,1:8)%判断不符合条件的顶点，若不符合则对应位置为0，反之，为1附录4%计算椭圆倾斜角，中心角array;k15=(array1(1,2)-array1(5,2))/(array1(1,1)-array1(5,1));ang15=atan(k15);x015=(array1(1,1)+array1(5,1))/2;y015=(array1(1,2)+array1(5,2))/2;ifang15<0ang15=pi+ang15;endk16=(array1(1,2)-array1(6,2))/(array1(1,1)-array1(6,1));ang16=atan(k16);x016=(array1(1,1)+array1(6,1))/2;y016=(array1(1,2)+array1(6,2))/2;ifang16<0ang16=pi+ang16;endk18=(array1(1,2)-array1(8,2))/(array1(1,1)-array1(8,1));ang18=atan(k18);x018=(array1(1,1)+array1(8,1))/2;y018=(array1(1,2)+array1(8,2))/2;ifang18<0ang18=pi+ang18;endk25=(array1(2,2)-array1(5,2))/(array1(2,1)-array1(5,1));ang25=atan(k25);x025=(array1(2,1)+array1(5,1))/2;y025=(array1(2,2)+array1(5,2))/2;ifang25<0ang25=pi+ang25;endk26=(array1(2,2)-array1(6,2))/(array1(2,1)-array1(6,1));ang26=atan(k26);x026=(array1(2,1)+array1(6,1))/2;y026=(array1(2,2)+array1(6,2))/2;ifang26<0ang26=pi+ang26;end37 k27=(array1(2,2)-array1(7,2))/(array1(2,1)-array1(7,1));ang27=atan(k27);x027=(array1(2,1)+array1(7,1))/2;y027=(array1(2,2)+array1(7,2))/2;ifang27<0ang27=pi+ang27;endk34=(array1(3,2)-array1(4,2))/(array1(3,1)-array1(4,1));ang34=atan(k34);x034=(array1(3,1)+array1(4,1))/2;y034=(array1(3,2)+array1(4,2))/2;ifang34<0ang34=pi+ang34;endk35=(array1(3,2)-array1(5,2))/(array1(3,1)-array1(5,1));ang35=atan(k35);x035=(array1(3,1)+array1(5,1))/2;y035=(array1(3,2)+array1(5,2))/2;ifang35<0ang35=pi+ang35;endk36=(array1(3,2)-array1(6,2))/(array1(3,1)-array1(6,1));ang36=atan(k36);x036=(array1(3,1)+array1(6,1))/2;y036=(array1(3,2)+array1(6,2))/2;ifang36<0ang36=pi+ang36;endk37=(array1(3,2)-array1(7,2))/(array1(3,1)-array1(7,1));ang37=atan(k37);x037=(array1(3,1)+array1(7,1))/2;y037=(array1(3,2)+array1(7,2))/2;ifang37<0ang37=pi+ang37;endk17=(array1(1,2)-array1(7,2))/(array1(1,1)-array1(7,1));ang17=pi+atan(k17);x17=(array1(1,1)+array1(7,1))/2;37 y17=(array1(1,2)+array1(7,2))/2;ang=[ang15ang16ang18ang25ang26ang27ang34ang35ang36ang37ang17];x0=[x015x016x018x025x026x027x034x035x036x037x17];y0=[y015y016y018y025y026y027y034y035y036y037y17];附录5%计算长半轴，短半轴c=[array2(1,5)array2(1,6)array2(1,8)array2(2,5)array2(2,6)array2(2,7)...array2(3,4)array2(3,5)array2(3,6)array2(3,7)array2(1,7)];c=c./2;a=c.*1.4;b=sqrt(a.^2-c.^2);附录6%画图函数%画椭圆functionellipse(x,y,a,b,ang)axis([02000100]);plot(y+a*sin(ang)*cos(0:pi/720:2*pi)+b*cos(ang)*sin(0:pi/720:2*pi),x+a*cos(ang)*cos(0:pi/720:2*pi)-b*sin(ang)*sin(0:pi/720:2*pi),"b");附录7画椭圆，在matlab命令窗口中输入代码：%画两个椭圆的程序ellipse(y0(11),x0(11),b(11),a(11),ang(11));holdon;ellipse(y0(8),x0(8),b(8),a(8),ang(8));holdoff;%画10个椭圆的程序fori=1:10ellipse(y0(i),x0(i),b(i),a(i),ang(i));holdon;end附录8%H形代码functionans=search1(x,y)min=inf;fori1=35:120forj1=1:99len1=figdis([i1x(2)],[j1y(2)])+figdis([i1x(6)],[j1y(6)]);if(len1+30<142)&(len1+25<150)fori2=i1:160forj2=1:99len2=figdis([i2x(5)],[j2y(5)])+figdis([i2x(3)],[j2y(3)]);%3105iflen2>151continue;37 end;len3=figdis([i2x(5)],[j2y(5)])+figdis([i1x(2)],[j1y(2)])+figdis([i1i2],[j1j2]);%29105len4=figdis([i1x(6)],[j1y(6)])+figdis([i2x(3)],[j2y(3)])+figdis([i1i2],[j1j2]);%31096iflen3>172continue;end;iflen4len1+110minl=len1+110;end;ifminl>222continue;end;ifminl+10>253continue;end;if110+len2>len3minl=len3;elseminl=len2+110;end;iflen1+85>minlminl=len1+85;end;ifminl>197continue;end;ifmin>len1+len2+figdis([i1i2],[j1j2])min=len1+len2+figdis([i1i2],[j1j2]);ani1=i1;ani2=i2;anj1=j1;anj2=j2;min1=len1;min2=len2;end;endendendendendans=[ani1anj1ani2anj2minmin1min2];附录9%双Y形算法%shuang_y.m双Y形程序short=inf;v12=array2(1,2);v23=array2(2,3);v56=array2(5,6);v67=array2(6,7);fori1=0:4:200forj1=100:-4:0fori2=i1:4:200forj2=0:4:100v2h=sqrt((i1-array1(2,1))^2+(j1-array1(2,2))^2);v6h=sqrt((i1-array1(6,1))^2+(j1-array1(6,2))^2);37 v5h=sqrt((i1-array1(5,1))^2+(j1-array1(5,2))^2);v3i=sqrt((i2-array1(3,1))^2+(j2-array1(3,2))^2);v4i=sqrt((i2-array1(4,1))^2+(j2-array1(4,2))^2);v5i=sqrt((i2-array1(5,1))^2+(j2-array1(5,2))^2);w=v2h+v5h+v6h+v3i+v4i+v5i;ifshort>=wifarray2(1,5)>=(v12+v2h+v5h)/1.4ifarray2(1,6)>=(v12+v2h+v6h)/1.4if(array2(2,5)>=(v2h+v5h)/1.4)|(array2(2,5)>=(v23+v3i+v5i)/1.4)ifarray2(2,6)>=(v2h+v6h)/1.4ifarray2(2,7)>=(v2h+v6h+v67)/1.4ifarray2(3,4)>=(v3i+v4i)/1.4if(array2(3,5)>=(v3i+v5i)/1.4)|array2(3,5)>=(v23+v2h+v5h)/1.4if(array2(3,6)>=(v3i+v5i+v56)/1.4)|(array2(3,6)>=(v23+v2h+v6h)/1.4)if(array2(3,7)>=(v3i+v5i+v56+v67)/1.4)|(array2(3,7)>=(v23+v2h+v6h+v67)/1.4)short=w;x1=i1;y1=j1;x2=i2;y2=j2;x=w;endendendendendendendendendendendendendend%计算总路程total_dis1=array2(1,8)+x;附录10%双Y形算法优化37 %shuang_y_youhua.m优化算法i01=x1-2;i02=x2-2;j01=y1-2;j02=y2-2;short=x;v12=array2(1,2);v23=array2(2,3);v56=array2(5,6);v67=array2(6,7);fori1=i01:0.1:i01+4forj1=j01+4:-0.1:j01fori2=i02:0.1:i02+4forj2=j02:0.1:j02+4v2h=sqrt((i1-array1(2,1))^2+(j1-array1(2,2))^2);v6h=sqrt((i1-array1(6,1))^2+(j1-array1(6,2))^2);v5h=sqrt((i1-array1(5,1))^2+(j1-array1(5,2))^2);v3i=sqrt((i2-array1(3,1))^2+(j2-array1(3,2))^2);v4i=sqrt((i2-array1(4,1))^2+(j2-array1(4,2))^2);v5i=sqrt((i2-array1(5,1))^2+(j2-array1(5,2))^2);w=v2h+v5h+v6h+v3i+v4i+v5i;ifshort>=wifarray2(1,5)>=(v12+v2h+v5h)/1.4ifarray2(1,6)>=(v12+v2h+v6h)/1.4if(array2(2,5)>=(v2h+v5h)/1.4)|(array2(2,5)>=(v23+v3i+v5i)/1.4)ifarray2(2,6)>=(v2h+v6h)/1.4ifarray2(2,7)>=(v2h+v6h+v67)/1.4ifarray2(3,4)>=(v3i+v4i)/1.4if(array2(3,5)>=(v3i+v5i)/1.4)|array2(3,5)>=(v23+v2h+v5h)/1.4if(array2(3,6)>=(v3i+v5i+v56)/1.4)|(array2(3,6)>=(v23+v2h+v6h)/1.4)if(array2(3,7)>=(v3i+v5i+v56+v67)/1.4)|(array2(3,7)>=(v23+v2h+v6h+v67)/1.4)short=w;x1=i1;y1=j1;x2=i2;y2=j2;x=w;endendendendendendendendend37 endendendendend%计算总路程total_dis2=array2(1,8)+x;附录11%三个点的程序%san_dian_youhua双Y形程序short=inf;v12=array2(1,2);v23=array2(2,3);v56=array2(5,6);v67=array2(6,7);fori1=0:10:200forj1=100:-10:0fori2=i1:10:200forj2=0:10:100fori3=0:10:200forj3=0:10:200v2h=sqrt((i1-array1(2,1))^2+(j1-array1(2,2))^2);v6h=sqrt((i1-array1(6,1))^2+(j1-array1(6,2))^2);vfh=sqrt((i1-i3)^2+(j1-j3)^2);v3i=sqrt((i2-array1(3,1))^2+(j2-array1(3,2))^2);v4i=sqrt((i2-array1(4,1))^2+(j2-array1(4,2))^2);vfi=sqrt((i2-i3)^2+(j2-j3)^2);v5f=sqrt((i3-array1(5,1))^2+(j3-array1(5,2))^2);w=v2h+v6h+vfh+v3i+v4i+vfi+v5f;ifshort>=wifarray2(1,5)>=(v12+v2h+vfh+v5f)/1.4ifarray2(1,6)>=(v12+v2h+v6h)/1.4if(array2(2,5)>=(v2h+vfh+v5f)/1.4)|(array2(2,5)>=(v23+v3i+vfi+v5f)/1.4)ifarray2(2,6)>=(v2h+v6h)/1.4ifarray2(2,7)>=(v2h+v6h+v67)/1.4ifarray2(3,4)>=(v3i+v4i)/1.4if(array2(3,5)>=(v3i+vfi+v5f)/1.4)|array2(3,5)>=(v23+v2h+vfh+v5f)/1.4if(array2(3,6)>=(v3i+vfi+v5f+v56)/1.4)|(array2(3,6)>=(v23+v2h+v6h)/1.4)|(array2(3,6)>=(v3i+vfi+vhf+v6h)/1.4)if(array2(3,7)>=(v3i+vfi+v5f37 +v56+v67)/1.4)|(array2(3,7)>=(v23+v2h+v6h+v67)/1.4)|(array2(3,7)>=(v3i+vfi+vhf+v6h+v67)/1.4)short=w;x1=i1;y1=j1;x2=i2;y2=j2;x3=i3;y3=j3;x=short;endendendendendendendendendendendendendendendend%计算总路程total_dis3=array2(1,8)+x;附录12%三个点中更加逼近的算法%san_dian_gengyouhua.m双Y形程序i01=x1-5;i02=x2-5;j01=y1-5;j02=y2-5;i03=x3-5;j03=y3-5;short=x;v12=array2(1,2);v23=array2(2,3);v56=array2(5,6);v67=array2(6,7);fori1=i01:i01+10forj1=j01:j01+1037 fori2=i1:i02+10forj2=j02:j02+10fori3=i03:i03+10forj3=j03:j03+10v2h=sqrt((i1-array1(2,1))^2+(j1-array1(2,2))^2);v6h=sqrt((i1-array1(6,1))^2+(j1-array1(6,2))^2);vfh=sqrt((i1-i3)^2+(j1-j3)^2);v3i=sqrt((i2-array1(3,1))^2+(j2-array1(3,2))^2);v4i=sqrt((i2-array1(4,1))^2+(j2-array1(4,2))^2);vfi=sqrt((i2-i3)^2+(j2-j3)^2);v5f=sqrt((i3-array1(5,1))^2+(j3-array1(5,2))^2);w=v2h+v6h+vfh+v3i+v4i+vfi+v5f;ifshort>=wifarray2(1,5)>=(v12+v2h+vfh+v5f)/1.4ifarray2(1,6)>=(v12+v2h+v6h)/1.4if(array2(2,5)>=(v2h+vfh+v5f)/1.4)|(array2(2,5)>=(v23+v3i+vfi+v5f)/1.4)ifarray2(2,6)>=(v2h+v6h)/1.4ifarray2(2,7)>=(v2h+v6h+v67)/1.4ifarray2(3,4)>=(v3i+v4i)/1.4if(array2(3,5)>=(v3i+vfi+v5f)/1.4)|array2(3,5)>=(v23+v2h+vfh+v5f)/1.4if(array2(3,6)>=(v3i+vfi+v5f+v56)/1.4)|(array2(3,6)>=(v23+v2h+v6h)/1.4)|(array2(3,6)>=(v3i+vfi+vhf+v6h)/1.4)if(array2(3,7)>=(v3i+vfi+v5f+v56+v67)/1.4)|(array2(3,7)>=(v23+v2h+v6h+v67)/1.4)|(array2(3,7)>=(v3i+vfi+vhf+v6h+v67)/1.4)short=w;x1=i1;y1=j1;x2=i2;y2=j2;x3=i3;y3=j3;37 x=short;endendendendendendendendendendendendendendendend%计算总路程total_dis4=array2(1,8)+x;附录13%第3、4两点到湖的距离之和最小求解%shortest34.m求到3、4点距离最短va=45;vb=35;short=va+vb;fori=160:0.01:200forj=0:0.01:50v3o=sqrt((i-array1(3,1))^2+(j-array1(3,2))^2);v4o=sqrt((i-array1(4,1))^2+(j-array1(4,2))^2);ifi<=165&j<=45vo=45-j;elseifj>=45&i>=165vo=i-165;elsevo=sqrt((i-165)^2+(j-45)^2);endendifv3o+v4o+vo89.6437continue;endlen31=110+figdis([x(2)x1],[y(2)y1])+figdis([x(6)x1],[y(6)y1]);len32=len1+figdis([x1x2],[y1y2])+figdis([x(6)x1],[y(6)y1]);len33=len1+figdis([x(5)x2],[y(5)y2])+85;len3=figmin([len31len32len33]);iflen3>224.1093continue;endlen3=figmin([len31+25len32+25len33+25]);iflen3>252.3886continue;endiflen33-85>150.7846continue;endlenans=32.0156+len32-len1+figdis([x2140],[y2j1])+figdis([x(5)x2],[y(5)y2])+figdis([x(2)x1],[y(2)y1])+figdis([x(3)i2],[045])+figdis([x(4)165],[y(4)j3]);ifmin>lenansmin=lenans;anj1=j1;anj3=j3;ani2=i2;endendendendans=[anj1ani2anj3min];附录15%搜索交叉点的位置37 functionans=backsort(x,y,lj1,lx2,lj3,ax1,ay1,ax2,ay2)min=inf;fori1=ax1-15:ax1+15forj1=ay1-15:ay1+15len1=figdis([i1x(2)],[j1y(2)])+figdis([i1x(6)],[j1y(6)]);if(len1+30<142)fori2=ax2-15:ax2+15forj2=ay2-15:ay2+15len2=figdis([x(3)lx2],[y(3)45])+lx2-140+lj1-45+figdis([i2140],[j2lj1]);iflen2+figdis([x(5)i2],[y(5)j2])>151continue;endlen3=figdis([i2x(5)],[j2y(5)])+figdis([i1x(2)],[j1y(2)])+figdis([i1i2],[j1j2]);len4=figdis([i1x(6)],[j1y(6)])+len2+figdis([i1i2],[j1j2]);iflen3>172continue;endiflen4len1+110minl=len1+110;endifminl>222continue;endiflen1+85>len3minl=len3;elseminl=len1+85;endifminl>197continue;endlenans=len1+len2+figdis([i1i2],[j1j2])+figdis([x(5)i2],[y(5)j2])-(lx2-140+lj1-45)+32.0156+figdis([x(4)165],[y(4)lj3]);%if(i1==ax1)&(j1==ay1)&(i2==ax2)&(j2==ay2)lenans%end;ifmin>lenansmin=lenans;ani1=i1;ani2=i2;anj1=j1;anj2=j2;endendendendendendans=[ani1anj1ani2anj2min];37 37