四川灾区巡视最优化路线的设计.doc

四川灾区巡视最优化路线的设计摘 要去年五月十二日四川遭遇了数十年未见的巨大地震，对当地的各方面都造成了难以计数的破坏，对这种情况，我们除了对死者表示哀悼外，还可以借助我们自己掌握的一点微薄之力，使用数学工具建立一些模型，让灾情的巡视路线设计的更合理，为有关部门救灾提供建议，为此我们按照题意设计了一组和多组巡视灾区的简略模型和详细模型，在模型中使用了直接汉密尔顿法，破圈法，模拟退火法，主观划分法，基于基点的划分法等。并通过均衡值法验证结论，求出了相对合理的结果。对于问题一，由于无法判断能否形成汉密尔顿圈，并且点数较少，我们就假设存在汉密尔顿圈的情况下，使用逻辑方法判断去求汉密尔顿圈，之后使用模拟退火法验证结论，结果证实了判断。对于问题二、三，我们主要采用基于基点的划分法来解决问题，结果均衡性都不错，这也说明了这一方法的可行性。关键词：汉密尔顿圈 模拟退火法 基于基点的划分法一 问题重述2008年5月12日四川发生大地震，下图为四川地震一截图。成都领导决定带领有关部门负责人到受灾地区各县（市）巡视，巡视路线指导从成都所在地出发，走遍各受灾县（市），又回到成都政府所在地的路线。路上合理考虑阻抗。问题1、假定巡视人员在各县停留时间t=1.5小时，市停留时间T=2.5小时，汽车行驶一般公路（国道）速度v=60公里/小时，行驶高速速度V=89公里/小时。若分一组巡视地图上的方框地区，如何设计巡视线路，最短时间为多少？若要在18小时内完成巡视，至少应分几组？请给出这种分组下你认为最佳的巡视路线。问题2、在上述关于t、v、T 和V的假定下，若分三组（路）巡视方框地区，试设计总路程最短且各组在巡视时间上尽可能均衡的巡视路线。问题3、考虑更为具体的受灾地区，如附表1所示（与问题1和问题2比较需要巡视的地区数增多），分2组考虑最佳巡视路线。 问题4、进一步考虑通用模型。已知有M个受灾地点，给定各点之间的距离和路况，在给定t、v、T 和V的假定下，若分2组（路）巡视受灾地点，试设计巡视时间最短的巡视路线。同学有兴趣的话，可以进一步考虑若分N组（路）巡视受灾地点，试设计巡视时间最短的巡视路线。附表1：四川主要受灾地区汶川县北川县绵竹市都江堰市广元市青川县成 都什邡市安 县平武县彭州市茂 县江油市理 县雅 安眉 山巴 中南 充遂 宁乐 山甘 孜广 安泸 州凉 山自 贡资 阳内 江说明：两地间的距离可以通过INTERNET网查询，路况可做一个大致估计，路上的行驶不需要很精确，但要大致合理。模型的假设1） 把各个县市看做质点。2） 所查资料真实可信。3） 巡视路线可以重复。4） 无任何意外发生巡视中断。5） 对于某些要经过多次的县市,只停留一次。6） 对于巡视时已巡视过的市县路过时不必停留。7） 对于工作人员休息问题，可认为巡视人员在路上休息。符号的说明1）图中黑线表示一般的公路和国道。2）图中加粗的红线表示高速公路。3）图中加粗的蓝线表示区域分界线。 3）各质点用表示。4）图中的数字表示两地间的行驶时间。5）：市停留时间；：县停留时间。6）各组巡视时间为，表示均衡度。（）问题一问题的分析 问题一要求总巡视时间最短，对于整个区域，可以采用破圈法求出一个最小生成树，再适当添加路线（所花时间尽量少），能使巡视路线构成一个回路，由此可求出一个走遍所有的县市的路线所花时间为29.69（小时），对于分两组巡视整个区域，要求巡视路线尽量均衡，可求出每组路线巡视时间，再分析它们的均衡度。 模型的建立与求解 各县市可经过的线路图如下 方法一 由于最后方案为汉密尔顿图，且点的数目不多，有几个点的度数为2，只有一条必经之路，如果遇到度数大于2的点，选择去掉多余的最终形成一个各点度数为2的接近汉密尔顿图的圈。所以下面我们列举几种可能方案首先编号：设县市为，成都为,德阳，绵阳，北川，茂县，汶川，映秀，都江堰，绵竹 ，什邡 ，彭州 。一号路线：成都都江堰,若走绵竹，无论走哪条路都会造成点的重复，走彭州，什邡一线也不行，所以只有一条线路可成，为成都都江堰映秀汶川茂县北川绵阳德阳绵竹什邡彭州成都二号路线：成都，如果选择彭州，则必须经什邡绵竹德阳绵阳北川茂县汶川映秀都江堰，最后回到成都。三号路线：成都德阳绵阳北川茂县汶川映秀都江堰绵竹什邡彭州成都 设巡视人员在各县停留时间=1.5小时，市停留时间=2.5小时，比较两个方案一号方案=0.74+2.5+0.76+1.5+0.97+1.5+0.66+1.5+1.56+1.5+1.14+2.5+0.64+2.5+0.62+2.5+0.38 +2.5+0.53+2.5+0.51=29.51h 二号方案=0.51+0.53+038+0.62+064+1.14+1.56+0.66+0.97+0.76+0.74+6+15=29.51h 三号方案 =0.69+2.5+0.64+2.5+1.14+1.5+1.56+1.5+0.66+1.5+0.97+1.5+0.76+2.5+1.61+2.5+0.38+2.5+0.53+2.5+0.51=30.45h 由于,所以选择一号方案。路线为:成都都江堰映秀汶川茂县北川绵阳德阳绵竹什邡彭州成都，最短时间为29.51小时。 方法二 根据运筹学中的破圈法，先选择成都都江堰绵竹市德阳市成都，由图中数据知都江堰绵竹市的行驶时间最长，故舍去，舍去后如下图所示：再选择北川县绵阳市德阳市绵竹市北川县，由上图中数据知，应舍去北川县至绵竹市的行驶时间为1.16的那条后如下图所示： 选择绵竹市德阳市什邡市绵竹市环形线路，舍去最长行驶时间路线：绵竹市德阳市，如下图所示：然后选择线路：成都彭州市什邡市德阳市成都，可看出成都德阳市路线所花时间最长，故舍去后如下图所示： 此时整个环路就剩下最小的一个圈，若舍去花时间0.74的路线：都江堰成都，巡视回来时就要走路线：都江堰彭州市成都，所花时间为0.6+0.51=1.110.74，故舍去不满足所花时间最短，应舍去路线：都江堰彭州市，如下图所示： 最后还剩最大的圈：成都都江堰映秀镇汶川县茂县北川县绵阳市德阳市什邡市彭州市成都，按破圈法应该舍去茂县北川县的路线。若舍去茂县北川县，巡视小组就得走来回两次路线，所花时间更长。此处，我们应灵活运用破圈法，故不舍去。综上所述，最终最短路线为：成都都江堰映秀镇汶川县茂县北川县绵阳市德阳市什邡市绵竹市什邡市彭州市成都（反向也可）。 将巡视各市县路段所花时间相加，有： =0.74+0.76+0.97+0.66+1.56+1.14+0.64+0.42+0.38+0.38+0.53+0.51=8.69而且巡视队员在各县市有停留，从所个图中可看出共有六个市，四个县，停留时间为：62.5+41.5=21（小时），得最短时间为8.69+21=29.69（小时）。方法三由于形成汉密尔顿圈的充分条件很难给出，只能采用其他方法去逼近最佳解，如果经过的点较多，计算量将极速增长。同时为了印证前面方法的正确性，我们借鉴了较为新颖的模拟退火法。模拟退火算法的主要思想 就函数最小值问题来说，模拟退火的主要思想是：在搜索区间（二维平面中）随机游走（即随机选择点），再以Metropolis抽样准则，使随机游走逐渐收敛于局部最优解。而温度即是Metropolis算法中的一个重要控制参数，可以认为这个参数的大小控制了随时过程向局部或全局最优解移动的快慢。 冷却参数表、领域结构和新解产生器、接受准则和随机数产生器（即Metropolis算法）一起构成算法的三大支柱。 l 重点抽样与Metroplis算法： Metropolis是一种有效的重点抽样法，其算法为：系统从能量一个状态变化到另一个状态时，相应的能量从E1变化到E2，概率为p = exp - (E2- E1)/kT 。如果E2 E1，系统接收此状态，否则，以一个随机的概率接收此或丢弃此状态。这种经常一定次数的迭代，系统会逐渐趋于一引稳定的分布状态。 重点抽样时，新状态下如果向下则接受（局部最优），若向上（全局搜索），以一定机率接受。模拟退火方法从某个初始解出发,经过大量解的变换后,可以求得给定控制参数值时组合优化问题的相对最优解。然后减小控制参数T的值，重复执行Metropolis算法，就可以在控制参数T趋于零时，最终求得组合优化问题的整体最优解。控制参数的值必须缓慢衰减。 其中温度是一个Metropolis的重要控制参数，模拟退火可视为递减控制参数什时Metroplis算法的迭代。开始T值大，可能接受较差的恶化解，随着T的减小，只能接受较好的恶化解，最后在T趋于0时，就不再接受任何恶化解了。 在无限高温时，系统立即均匀分布，接受所有提出的变换。T的衰减越小，T到达终点的时间越长；但可使马可夫链越小，到达准平衡分布的时间越短， 参数的选择： 我们称调整模拟退火法的一系列重要参数为冷却进度表。它控制参数T的初值及其衰减函数，对应的MARKOV链长度和停止条件，非常重要。 一个冷却进度表应当规定下述参数： 1 控制参数t的初值t0; 2. 控制参数t的衰减函数； 3.马尔可夫链的长度Lk。（即每一次随机游走过程，要迭代多少次，才能趋于一个准平衡分布，即一个局部收敛解位置） 4 结束条件的选择 有效的冷却进度表判据： 一算法的收敛：主要取决于衰减函数和马可夫链的长度及停止准则的选择 二算法的实验性能：最终解的质量和CPU的时间 参数的选择： 一）控制参数初值T0的选取 一般要求初始值t0的值要充分大，即一开始即处于高温状态，且Metropolis的接收率约为1。 二）衰减函数的选取 衰减函数用于控制温度的退火速度，一个常用的函数为：T(n + 1) = K*T(n)，其中K是一个非常接近于1的常数。 三）马可夫链长度L的选取 原则是：在衰减参数T的衰减函数已选定的前提下，L应选得在控制参数的每一取值上都能恢复准平衡。 四）终止条件 有很多种终止条件的选择，各种不同的条件对算法的性能和解的质量有很大影响，我们选择一个最优解与最新的一个最优解的之差小于某个容差，即可停止此次马尔可夫链的迭代。 步1： 设定初始温度T，给定一个初始的巡视路线。 步2 ：步3 -8循环K次 步3：步 4-7循环M次 步4：随机选择路线的一段 步5：随机确定将选定的路线反转或移动，即两种调整方式：反转、移动。 步6：计算代价D，即调整前后的总路程的长度之差 步7：按照如下规则确定是否做调整： 如果D0,则按照EXP(-D/T)的概率进行调整 步8:T*0.9-T,降温 运行结果为：成都德阳绵阳北川茂县汶川映秀都江堰绵竹什邡彭州成都，总巡视时间为29.51，这样就和方法一中的方案一结果一致要想在18小时内完成巡视线路，而以上三种方法算出时间均大于18，显然分一组是不可能的。接着考虑分组的情况主观性分析：若假定分为两组巡视回路，首先将整个区域分成两部分，分区域的依据是：各组在各自区域内所经过的路线与所花时间约占总路线和总时间的二分之一。区域划分图如下图所示：第一组路线：成都都江堰映秀镇汶川茂县北川茂县汶川映秀镇都江堰成都，所需的巡视时间为： =(0.74+0.76+0.97+0.66+1.56)*2+8.5=17.88（小时）第二组路线：成都彭州市什邡市绵竹市德阳绵阳德阳成都，所需的巡视时间为： = 0.51+0.53+0.38+0.62+0.64+1.33+12.5=16.51（小时） 可见 ，显然时间是够的，从而分成两组是可行的。 定义均衡度为 注：均衡度越小，分配的线路越合理。根据文献及相关资料，当均衡度小于25%时，我们认为模型是合理的。代入得：基于基点的划分方法 根据题目中所述和前种方法的估计，我们的工作是给出集合的一个分组，目标函数是： 总路线最短 两组尽量均衡 但是，要给出这样一个优化问题的最优解是十分困难的。可以想象，当两组的总时间比较均匀时，两组的总时间不会离最优解太远。我们对这个问题的解法是首先根据一定的原则给出初始的划分，当然，我们希望这个初始分组离最优解不太远，然后计算各组的路线长度，再根据一定的原则进行调整，使各组尽量均匀，并使路线总和也比较理想。从巡视路线图中,先确定基点，即离目的地最远的点，由于每组都是要经过成都的，所以取成都为一个基点，其它的基点取与成都距离最远的点，接着将离各基点范围距离最近的点归于相应基点范围内，各基点轮流取至范围内点，直至点分完为止，如有少于基点数的点余下，则视其与各基点的距离远近决定归属，最后将成都这一基点范围内的点按距离及实际情况分给其余基点范围，成都为各基点的公共点。从上面方法我们也假定划分为两组，那么取三个基点，其中一个为成都，其他的基点为汶川，绵阳，分属的点为成都，彭州，都江堰，什邡；映秀，汶川，茂县；绵阳，德阳，绵竹，北川三组，然后将成都这组的点调整到其他组，结果为成都，彭州，都江堰，映秀，汶川，茂县；成都，什邡，绵阳，德阳，绵竹，北川两部分。如下图：第一组路线：成都都江堰映秀镇汶川茂县汶川映秀镇都江堰彭州市成都，所需的巡视时间为： =0.66+0.97+0.76)*2+0.74+0.6+0.51+4.5+5= 16.13（小时）第二组路线：成都德阳绵阳北川绵竹什邡市德阳成都，所需的巡视时间为： =0.69*2+0.64+1.14+1.16+0.38+0.42+10+1.5=16.62（小时）可见 ，显然时间是够的，从而分成两组是可行的。 两组的均衡度为：（二）问题二 问题的分析与求解 采用与前一小题相同的方法。然后在本题中，要求划分三组，那么取四个基点，其中一个为成都，其他的基点为汶川，什邡，绵阳，分属的点为成都，彭州，都江堰；映秀，汶川，茂县；北川，绵阳，德阳；绵竹，什邡四组，然后将成都这组的点调整到其他组，结果是成都德阳绵阳北川；成都彭州什邡绵竹；成都都江堰映秀镇汶川茂县三部分。 就可将区域分为，如下图：可看出每组均沿原路返回第一组巡视的时间： (0.74+0.76+0.97+0.66)*2+4.5+2.5 = 13.26（小时）第二组巡视的时间： (0.69+0.64+1.14)*2+5+1.5=11.44 （小时）第三组巡视的时间： (0.51+0.53+0.38)*2+7.5 =10.34 （小时）此三组的均衡度为：（3） 问题三问题的分析与求解在原有的问题上扩展所需要巡视市县的数目，根据附件1可做出各市县之间的线路图，如下：根据题意可知该题与第一小题的后半小题极为相似，可继续使用基于基点的划分方法，取基点为成都，茂县，广元。逐次调整，可将区域分为,如下图： 第一组线路：成都眉山乐山雅安成都都江堰汶川茂县理县茂县北川平武青川北川绵竹什邡彭州成都，所需的巡视时间为： =47.74 (小时)第二组路线：成都自贡泸州内江资阳遂宁广安南充巴中广元安县江油安县成都，所需巡视时间为： =50.58（小时）此两组的均衡度为： （四）模型的评价与改进本模型运用图论中最佳推销员路径问题的相关结论，建立了关于该类问题的严格的优化模型，但是可以说求得这类问题的精确最优解是十分困难的。因此本模型为了解决这类问题采用了近似的算法，得到的结果比较理想。2 . 模型主观性相对较强，求出的均衡度比较符合要求。而且由于信息不够完善，本文在建模中考虑的因素不多，比如在巡视过程中的阻抗问题未加以考虑。本文中的模型适用范围相对较窄，对市县较多等情况下的通用模型解决力度不够，耗时较多。对此可以采用更先进的算法或原理，比如遗传算法，神经网络算法，对问题加以解决。3. 附录：#include #include #define TFACTR 0.9#define MAX 10000.0#define MBIG 1000000000#define MSEED 161803398#define MZ 0#define FAC (1.0/MBIG)float weight1212;int iorder11;void anneal(int ncity) int irbit(unsigned long *iseed); int metrop(float de,float t); float ran(long *idum); float reverse_len(int ncity,int n); void reverse(int ncity,int n); float transe_len(int ncity,int n); void transe(int ncity,int n); int ans,nover,nlimit,i1,i2; int i,j,nsucc,nn,idec; unsigned long k; static int n7; long idum; unsigned long iseed; float path,de,t,T; nlimit=100*ncity;/*max succ try numbers */ path =0.0; for(i=1;incity;i+) /* primitive path length */ i1=iorderi; i2=iorderi+1; path+=weighti1i2; i1=iorderncity; i2=iorder1; path+=weighti1i2; idum=-1; iseed=111; T=2.8; for(j=1;j=30;j+) nsucc=0; for(k=1;k=n1) +n2; nn=1+(n1-n2+ncity-1)%ncity); while(nn=nlimit) break; printf(n %s %10.6f %s %10.6f j=%d,T=,T,path len=,path,j); printf(successful moves: %6dn,nsucc); for(i=1;i=11;i+) printf(%5d,iorderi); T*=TFACTR; if(nsucc=0) return; float reverse_len(int ncity,int n) float de; int xx5,yy5; int j,ii; n3=1+(n1+ncity-2)%ncity); n4=1+(n2%ncity); for(j=1;j=4;j+) xxj=iordernj; de=-weightxx1xx3; de-=weightxx2xx4; de+=weightxx1xx4; de+=weightxx2xx3; return de;void reverse (int ncity,int n) int nn,j,k,l,itmp; nn=(1+(n2-n1+ncity)%ncity)/2; for(j=1;j=nn;j+) k=1+(n1+j-2)%ncity); l=1+(n2-j+ncity)%ncity); itmp=iorderk; iorderk=iorderl; iorderl=itmp; float transe_len(int ncity,int n) float /*xx7,yy7,*/ de; int j,ii,xx7; n4=1+(n3%ncity); n5=1+(n1+ncity-2)%ncity); n6=1+(n2%ncity); for(j=1;j=6;j+) /*ii=iordernj; xxj=xii; yyj=yii;*/ xxj=iordernj; de=-weightxx2xx6; de-=weightxx1xx5; de-=weightxx3xx4; de+=weightxx1xx3; de+=weightxx2xx4; de+=weightxx5xx6; return de;void transe(int ncity,int n) int m1,m2,m3,nn,j,jj,jorder100; m1=1+(n2-n1+ncity) % ncity); /* n1 to n2 city numbers*/ m2=1+(n5-n4+ncity) % ncity); /*n4 to n5 */ m3=1+(n3-n6+ncity) % ncity); /*n3. n6 .*/ nn=1; for(j=1;j0) for(j=1;j0) for(j=1;j=m3;j+) jj=1+(j+n6-2)%ncity); jordernn+=iorderjj; for(j=1;j=ncity;j+) iorderj=jorderj;int metrop(float de,float t) float ran(long *idum); static long gljdum=1; return de0.0 |ran(&gljdum)exp(-de/t);float ran(long *idum) static int inext,inextp; static long ma56; static int iff=0; long mj,mk; int i,ii,k; if(*idum0|iff=0) iff=1; mj=MSEED-(*idum0?-*idum:*idum); mj%=MBIG; ma55=mj; mk=1; for(i=1;i54;i+) ii=(21*i)%55; maii=mk; mk=mj-mk; if(mkMZ) mk+=MBIG; mj=maii; for(k=1;k=4;k+) for(i=1;i=55;i+) ma