五一数学建模A题不确定性下的最短路径问题CUMT赖增强

1、2 0 15 年暑期数学建模培训第一次模拟五一数学建模A题不确定性下的最短路径问题CUMT赖增强文档编制序号：KK8UY-LL9IO69-TTO6M3-MTOL89-FTT688承诺书我们仔细阅读了数学建模联赛的竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、 电子邮件、网上咨询等）与本队以外的任何人（包括指导教师）研究、 讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的，如果引用别人的成果 或其它公开的资料（包括网上查到的资料），必须按照规定的参考文献 的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如

2、 有违反竞赛规则的行为，我们愿意承担由此引起的一切后果。我们授权数学建模联赛赛组委会，可将我们的论文以任何形式进行公 开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正 式发表等）。我们参赛选择的题号为（从A/B/C中选择一项填写）：我们的参赛报名号为：参赛组别（研究生或本科或专科）： wt所属学校（请填写完整的全名）中国矿业大学南湖校区参赛队员（打印并签名）：1.赖增强2 .兰卫旗3 .李康杰口期：2015 年 8 月 11 B获奖证书邮寄地址：中国矿业大学南湖校区桃4B5032邮政编码: 22nl6收件人姓名：赖增强 联系电话：183612308562015年暑期数学建模培训第

3、一次模拟编号专用页竞赛评阅编号（由竞赛评委会评阅前进行编号）： 评阅记录评 阅 人评分备注裁剪线裁剪线裁剪线竞赛评阅编号（由竞赛评委会评阅前进行编号）：参赛队伍的参赛号码：1（请各参赛队提前填写好）：不确定条件下的最优路径问题摘要本文针对如何在复杂的交通环境下寻找一条可靠、快速、安全的最优 路径的问题，考虑到交通堵塞、恶劣天气、路途成本等不确定因素对司 机路径选择的影响，建立多个不确定条件下的最优路径模型。对于问题一，我们在各个路段所用时间服从正态分布N(，”)的 基础上，建立了在不确定条件下求最短路的NP模型，给每个路段设定一 个预留到达的时间t,为了尽可能准确的到达目的地，选取95%的概率

4、， 满足PTt?95%,那么最优路径的定义就是预留时间最小的那个路径， 将其转换为标准的正态分布，通过标准的正态分布得到了在不确定性条 件下车辆从起点到终点预留时间的数学表达式：廿 十76。计算得对 应的t (绕城)=34. 645min, t (市区)=54. 675min,那么最优路径为绕 城快速路。对于问题二，在第一问定义的基础上进一步引入Bool系数B(a,k),在 搜集得到的具体的交通网络中，建立了一个从起点到终点路径为 EL B(a,k)T：nk的正态分布，通过求最小预留时间 小)之田叫+ 0)-小,谈力，得出最优路径的算法。其中E,泮”= Xa=i B (% k)E Tlink

5、, VarTath = Za=i P 3, k) VarT产，fflVarTjath 的 根式不具有线性可加性。不能用经典的dijkstra算法求解。对此采用基 于双目标规划的思路，利用第K短路径算法，分别对ETrh,运用matlab编程，找出各自前十条最短路径。之后在其并集 中找出最优路径：Vl-*V3-*V4-*V8o由此建立了求最短路的NPK模型。 最后从时间的渐进性态上分析模型的复杂性和收敛性。对于问题三，我们只考虑各路段空间上的相关性，并用概率论中的 协方差来表示这种耦合关系，建立了 NPK模型。得出可靠时间的数学表 达式廿 ETrh-1( P S(a,k)VarTlink + =2

6、cov(a- l,a)；求 解得最优路径：Vl->V3->V5->V8o关键词：NPK模型，K-短路，预留时间，正态分布，渐进性态。一.问题重述目前，交通拥挤和事故正越来越严重的困扰着城市交通，如何在复 杂的交通环境寻找一条可靠、快速、安全的最优路径，已经成为所有驾 驶员的共识。传统的最优路径就是行驶时间最短的路径，这是基于理想 交通情况下分析的，而事实上，在现实生活中，受到很多不确定性因素 的影响，例如：交通工具、恶劣天气、突发事件，导致车辆的行驶时间 均存在不确定性。为了减少交通阻塞所耽误的时间，尽可能快的到达目 的地，解决不确定性性条件下的最优路径问题，现依次提出以下问

7、题：（一）针对一般的交通网络，假设已知每条路段行驶时间的均值和 标准差，请建立相关的数学模型，定量地分析车辆行驶时间的不确定 性，然后给出在不确定性条件下车辆从起点到终点的最优路径的定义和 数学表达式。并将此模型应用到图1的例子中会选择哪条道路。（二）根据第一问的定义，假设已知每条路段行驶时间的均值和标准差，设计算法搜索最优路径，并将该算法应用到具体的交通网络中,用计算结果验证算法的有效性。从理论上分析算法的收敛性、复杂性等 性质。（三）建立数学模型，描述下游路段发生交通拥堵使车辆减速或者 排队，导致上游路段发生拥堵的交通路段之间驾驶时间的相关性，并将 此相关性应用到第一问和第二问的最优路径搜

8、索问题中，并设计算法解 决考虑相关性的最优路径搜索问题，给出算例验证算法的有效性。从理 论上分析算法的收敛性、复杂性等问题。二.基本假设符号说明3. 1基本假设3.1. 2假设仅相邻两条路段之间具有相关性。3.2. 3假设任意两条相邻路段组成的协方差矩阵为一个随机的正半定矩阵。3. 2符号说明3.2. 1ylink两相邻节点之间的路段3.2.2j'path从起点到终点的路径3.2.3Bool系数。当路径a通过此路段k时为3.2.4VarT对随机变量T求方差运算3.2.5标准正太分布累计概率函数的逆函数3.2.6cov(tl, t2)随机向量tl, t2协方差四.模型的建立与求解4. 1

9、问题一分析本问题是对于题设的交通网络，已知每条路段行驶时间的均值和标准差6条件下。定量分析车辆行驶时间的不确定性，以及给出在不 确定性条件下最优路径的定义和数学表达式。假设各个路段所用时间服从正态分布N（八”）模型，利用MATLAB 模拟（附录一）两条路径的正态分布图：）给每个路段设定一个预留到达的时间t,为了尽可能准确的到达目的地，选取95%的概率，满足PTWt295机 那么最优路径的定义就是 预留时间最小的那个路径。建立与求解已知到达目的地所用时间和预留时间满足:PTWt295%,将其转 换为标准的正态分布：P?< 字295乐 得到。 （字）295%, 二幺2。0），（其中P为95%

10、）即可解得每条路段的tN +l. 645 5,取廿+l. 645 5。我们将上诉模型称之为不确定条件下的 NP模型。应用已知的数学表达式，将图1所给的数据： 1=33,61=1； 2=30,6 2=15；带入公式计算出：t （绕城）=34. 645min, t （市区）=54. 675min,最优路径为绕城快速路。4.2问题二根据第一问中的最优路径定义和相关数学表达式，对于一般交通网络，可以结合K短路径算法建立NPK模型，最后从时间的渐进性态上分析其复杂性和收敛性。与求解对于一般交通网络，为了方便设计算法找到最优路径，我们可以将其转化为图论问题。将八个城市看做节点构成一个节点集：V=VI, V

11、2, V3, V4, V5, V6, V7, V8各个城市之间的道路看做边集：E=V1-*V2, V1-V3, Vl-*4, V1-V5, V1-V6, V1-V7,V2-V3, V1-*V4, V2-V7, V3-V4, V3-V5, V4-V5,V4-V8, V5-V6, V4-V8, V6-*V7, V6-V8, V7-*V8八个城市之间交通网络数据图图在第一问定义的基础上，针对每条路段引入Bool系数B(a,k),当该 路段被选择时为1,否则为0。那么从起点到终点的路径可表示为 EL B(a,k)T可知其?从正态分布。通过求该路径的最小预留时间 t心二E或叫+旷1 (p) JvarTj

12、ath，得出最优路径。对于产ET广町+t( p)JvarT泮7 其中：ETjath = Ea=i P(a,k) ET：nkVarT£ath= XL P (a, k)VarT产由于VarTh的根式不具有线性可加性。不能用经典的dijkstra 算法求解。对此我们基于双目标规划的思路，分别将ET*h, VarT泮?作为边权，运用matlab编程(附录二)，求出前十条最短路 径。表表根据上述两表数据，运用公式：t(Bin)=ETrh-1(P)JvarTath求出并集中的前16条最优路径的预留时间(表为了直观进行对比，将上表用excel制得如下柱状图：(图)由图可知最优路径为Vl-*V3-&

13、gt;V4-*V8o(图)收敛性分析：两个多项式时间算法之和还为多项式时间算法，其复杂 性比列举的低。当问题的规模趋向无穷大时，时间复杂度的数量级将表 现为渐进性态。即当路径K趋于无穷时，该模型一定收敛。4. 3问题三分析在问题三中，要求进一步考虑各路段之间行驶时间的相关性。我们用 概率论中的协方差来表示这种耦合关系，并建立推广的NPK模型。与求解在问题二中我们已得出最短预留时间的数学表达式：仁产田叫+k1()鼠口泮叮为方便模型的建立与求解。在此我们假设仅相邻两条路段之间具有 相关性。根据协方差的性质Var (T1+T2)=Var (Tl)+Var (T2) +2cov (Tl, T2);可以

14、得出t= E懊叫+8 39皿丁产+比=2皿Q一1闻称t为可靠时间。以下图为例：,为从A到B的一条路径。可靠时间t = E匹叫+t ( p ) JvarTath,其中ET泮阿二eT1+ET2, JvarTath=VVarTl + VarT2 + covTl + T2o由于问题二中，我们没有给出任意两条路段其时间随机向量（tl, t2） 的密度函数。无法具体求出协方差covtl+t2。对此我们假设任意两条 相邻路段组成的协方差矩阵为一个随机的正半定矩阵。在matlab中随机 函数rand （）基础上得出如下协差阵（附录三）：表对问题二找出的预留时间最小的前16条路分别求其可靠时间（表）：表运用ex

15、cel制图直观比较可靠时间大小（表）：表故得出最优路径为第三条V1->V3->V5->V8.（图）表五.模型评价本论文针对在不确定性条件下求解最优路径的问题，建立了以求最小 预留时间t为目标的NP模型，并对问题一给出了合理的解答。在此基础 上运用双目标规划的思想,结合求k短路径的方法，在没有考虑非相邻路 段间相关性基础上，针对更一般的问题建立了推广的NPK模型。复杂性 低，并随k增大具有较强收敛性。但在求均值与方差最短的并路径时， 没有设计出相应的算法，且本文只针对k较大时有效。六.参考文献1袁东，肖广冰.详解matlab快速入门与应用.北京：电子工业出版社,2011,73-

16、80.2邵虎，林兴强,孟强,谭美琳.基于出行问题可靠性的交通分配流问题.管理科学学报.2009. 12（5） : 1-4.3百度百科，多目标规划，2015-8-11.4桂云丽.变分不等式的算法研究.西安电子科技大学硕士论文.2010:1-8.附录一.为二维正太分布图function Y=funl(x);Y=(1/(2*pi*l)*exp(-(x-33)-2/(2*1*1);Y=(l/(2*pi*15) *exp(-(x-30) *2/(2*15*15);subplot (1,2,1);x y=meshgrid(25:0. 1:40);z = l/(2*pi*l).*exp(-(x-33). *

17、2/(2*1*1);h= mesh (x, y, z);set(h, ' edgecolor' , ' none' , ' facecolor' , ' interp，);subplot(1, 2, 2);x y=meshgrid(-50:0. 1:100);z = l/(2*pi*15).*exp(-(x-30). 2/(2*15*15);h= mesh (x, y, z);set(h, ' edgecolor' , none' , ' facecolor' , ' interp?);附录

18、二%第K短路算法function shortestPaths, totalCosts=kShortestPath(netCostMatrix, source, destination, k_paths)if source > size(netCostMatrix, 1) destination > size (netCostMatrix, 1)warning(, The source or destination node are not part of netCostMatrix,);shortestPaths=;totalCosts=;elsek=l;path cost = d

19、ijkstra(netCostMatrix, source, destination);%P is a cell array that holds all the paths found so far:if isempty(path)shortestPaths= ;totalCosts=;elsepath_number = 1;Ppath_number, Ij = path; Ppath_number, 2 = cost;current_P = path_number;%X is a cell array of a subset of P (used by Yen's algorith

20、m below):size_X=l;Xsize_X = path_number; path; cost;%S path_number x 1S(path_number) = path(1); adeviation vertex is the first node initially% K = 1 is the shortestpath returned by dijkstra():shortestPathsk= pathtotalCosts(k)= cost;while (k < k_paths &&size X = 0 )%remove P from Xfor i=l：

21、length(X)if Xi1 = current_Psize_X = size_X - 1;X(i) = 口%delete cellbreak;endendP_ 二 Pcurrent_P, 1) ; %P_ is current P, just to make is easier for the notations%Find w in (P_, w) in set S, w was the dev vertex used to found P_w = S (current-?);for i = 1: length(P_)if w = P_(i)w_index_in_path = i;ende

22、ndfor index_dev_vertex= w_index_in_path： length(P_) 一1 %index_dev_vertex is index in P_ of deviation vertex temp_netCostMatrix = netCostMatrix;%Remove vertices in P before index_dev_vertex and there incident edgesfor i = l： index_dev_vertex-lV = P_ ；temp_netCostMatrix (v, O=inf;temp_netCostMatrix(:,

23、 v)=inf;end%remove incident edge of v if v is in shortestPaths (K) U P with similar sub_path to P_.SP_sameSubPath=;index =1;SP-sameSubPathindex=P_；for 1 = 1： length(shortestPaths)if length(shortestPathsi) = index_dev_vertexif P_(l：index_dev_vertex)=shortestPathsi(1：index_dev_vertex)index = index+1;S

24、P_sameSubPath(index=shortestPaths i;endendendv_ = P_(index_dev_vertex);or j = 1： length(SP_sameSubPath)next 二 SP_sameSubPathj(index_dev_vertex+l);temp_netCostMatrix(v_, next)=inf;end%get the cost of the sub path before deviation vertex vsub_P = P_(l：index_dev_vertex);cost sub P=0;for i = 1: length(s

25、ub_P)-1cost_sub_P = cost_sub_P + netCostMatrix (sub_P (i), sub_P(i+D); end%call dijkstra between deviation vertex to destination nodedev_p c= dijkstra (temp_netCostMatrix,P_(index_dev_vertex), destination);if isempty(dev_p)path_number = path_number + 1;Ppath_number, 1 = sub_P(l:end-1) dev_p ; %conca

26、tenatesub path- to -vertex -to- destinationP path_number, 2)二cost_sub_P + c ;S(path-number) = P_(index_dev_vertex);size_X = size_X + 1;Xsize_X = path_number;Ppath-number, 1 ;Ppath_number, 2) );else%warning (' k=/d, isempty (p) true! n , k);endend%Step necessary otherwise if k is bigger than numb

27、er of possible paths%the last results will get repeated !if size_X > 0shortestXCost= X13; %cost of pathshortestX= X11;%ref number of pathfor i = 2 : size_Xif Xi3 < shortestXCostshortestX= Xi 1;shortestXCost= Xi3;endendcurrent P = shortestX;k = k+1;shortestPaths kJ 二 P current_P, 1;totalCosts(k

28、)= P current_P, 2;else%k = k+1;endendendendfunction shortestPath, totalCost = dijkstra(netCostMatrix,S, d)n = size (netCostMatrix, 1);for i = 1:n% initialize the farthest node to be itself;farthestPrevHop(i) = i; % used to compute the RTS/CTS range;farthestNextHop (i)= i;end% all the nodes are un-vi

29、sited;visited(1：n)= false;distance(1:n) = inf;% it stores the shortest distancebetween each node and the source node;parent(1:n)= 0;distance (s)= 0;for i = 1:(n-1),temp = 口;for h = 1：n,if visited(h) % in the tree;temp=temp distance(h);elsetemp=temp inf;endend;t, u = min (temp) ; % it starts from nod

30、e with the shortest distance to the source;visited(u)= true; % mark it as visited;for v = 1:n, % for each neighbors of node u;if ( ( netCostMatrix(u, v) + distance(u) < distance(v) distance (v) = distance(u) + netCostMatrix(u, v);% updatethe shortest distance when a shorter shortestPath is found;

31、 parent(v) = u; % update its parent;end;end;end;shortestPath =口;if parent (d) = 0% if there is a shortestPath!t = d；shortestPath =d;while t = sp = parent (t);shortestPath =p shortestPath;if netCostMatrix (t, farthestPrevHop (t) < netCostMatrix (t, p) farthestPrevHop(t)= p;end;if netCostMatrix(p,

32、farthestNextHop(p) < netCostMatrix(p, t) farthestNextHop(p)= t;end;t = p;end;end;totalCost = distance(d);%return;function = TestKShortestPath(case_number)switch 2case 1netCostMatrix = inf 10 20 inf inf 40 80 inf;10 inf 20 30 inf inf 70 inf;20 20 inf 10 50 inf inf inf;inf 30 10 inf 30 inf inf 60;7

33、0 inf 50 30 inf 30 inf 40;40 inf inf inf 30 inf 20 60;80 70 inf inf inf 20 inf 40;inf inf inf 60 40 60 40 infJ;source=l;destination=8;k = 10;case 2netCostMatrix =inf 4 25 inf 36 225 25 inf;4 inf 100 25 inf inf 225 inf;25 100 inf 49 4 inf inf inf;inf 25 49 inf 144 inf inf 100;36 inf 4 144 inf 225 inf 4;225 inf inf inf 225 inf 64 100;25 225 inf inf inf 64 inf 25;inf inf inf 100 4 100 25 inf;source=l;destination=8;k = 10;otherwiseerror (/ The only case options available are 1, 2 o