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

2015年暑期数学建模培训第一次模拟承诺书我们仔细阅读了数学建模联赛的竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与本队以外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其它公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们愿意承担由此引起的一切后果。我们授权数学建模联赛赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。 我们参赛选择的题号为（从A/B/C中选择一项填写）：我们的参赛报名号为：参赛组别（研究生或本科或专科）：本科所属学校（请填写完整的全名）中国矿业大学南湖校区参赛队员(打印并签名)：1.赖增强2.兰卫旗3.李康杰日期：2015年8月11日获奖证书邮寄地址：中国矿业大学南湖校区桃4B5032邮政编码：221116收件人姓名：赖增强联系电话/p>

2015年暑期数学建模培训第一次模拟编号专用页竞赛评阅编号（由竞赛评委会评阅前进行编号）：评阅记录评阅人评分备注裁剪线裁剪线裁剪线竞赛评阅编号（由竞赛评委会评阅前进行编号）：参赛队伍的参赛号码：1（请各参赛队提前填写好）：不确定条件下的最优路径问题摘要本文针对如何在复杂的交通环境下寻找一条可靠、快速、安全的最优路径的问题，考虑到交通堵塞、恶劣天气、路途成本等不确定因素对司机路径选择的影响，建立多个不确定条件下的最优路径模型。对于问题一，我们在各个路段所用时间服从正态分布N(μ，δ2)的基础上,建立了在不确定条件下求最短路的NP模型，给每个路段设定一个预留到达的时间t，为了尽可能准确的到达目的地，选取95%的概率，满足P{T≤t}≫95%,那么最优路径的定义就是预留时间最小的那个路径，将其转换为标准的正态分布，通过标准的正态分布得到了在不确定性条件下车辆从起点到终点预留时间的数学表达式：t=μ+Φ-1δ。计算得对应对于问题二，在第一问定义的基础上进一步引入Bool系数β(a,k),在搜集得到的具体的交通网络中，建立了一个从起点到终点路径为a=1nβa,kTalink的正态分布，通过求最小预留时间t(min)=E[Tkpath]+Φ-1Var[Tkpath],得出最优路径的算法。其中ETKpath=a=1nβa,kE[Talink]，Var[Tkpath]=a=1对于问题三，我们只考虑各路段空间上的相关性，并用概率论中的协方差来表示这种耦合关系，建立了NPK模型。得出可靠时间的数学表达式t=ETKpath+Φ关键词：NPK模型，K-短路，预留时间，正态分布，渐进性态。一.问题重述目前，交通拥挤和事故正越来越严重的困扰着城市交通,如何在复杂的交通环境寻找一条可靠、快速、安全的最优路径，已经成为所有驾驶员的共识。传统的最优路径就是行驶时间最短的路径，这是基于理想交通情况下分析的，而事实上，在现实生活中，受到很多不确定性因素的影响，例如：交通工具、恶劣天气、突发事件，导致车辆的行驶时间均存在不确定性。为了减少交通阻塞所耽误的时间，尽可能快的到达目的地，解决不确定性性条件下的最优路径问题，现依次提出以下问题：（一）针对一般的交通网络，假设已知每条路段行驶时间的均值和标准差，请建立相关的数学模型，定量地分析车辆行驶时间的不确定性，然后给出在不确定性条件下车辆从起点到终点的最优路径的定义和数学表达式。并将此模型应用到图1的例子中会选择哪条道路。起点起点终点市区道路均值30分钟，标准差15分钟图.1中国矿业大学徐州火车站绕城快速路均值33分钟，标准差1分钟（二）根据第一问的定义，假设已知每条路段行驶时间的均值和标准差，设计算法搜索最优路径，并将该算法应用到具体的交通网络中，用计算结果验证算法的有效性。从理论上分析算法的收敛性、复杂性等性质。（三）建立数学模型，描述下游路段发生交通拥堵使车辆减速或者排队，导致上游路段发生拥堵的交通路段之间驾驶时间的相关性，并将此相关性应用到第一问和第二问的最优路径搜索问题中，并设计算法解决考虑相关性的最优路径搜索问题，给出算例验证算法的有效性。从理论上分析算法的收敛性、复杂性等问题。二.基本假设符号说明3.1基本假设3.1.1各个路段所用时间服从正态分布。3.1.2假设仅相邻两条路段之间具有相关性。3.1.3假设任意两条相邻路段组成的协方差矩阵为一个随机的正半定矩阵。3.2符号说明3.2.1Tnlink3.2.2Tmpath3.2.3βa,kBool系数。当路径a通过此路段3.2.4Var[T]对随机变量T求方差运算3.2.5Φ-1（ρ）3.2.6cov(t1,t2)随机向量t1,t2协方差四.模型的建立与求解4.1问题一4.1.1问题分析本问题是对于题设的交通网络，已知每条路段行驶时间的均值μ和标准差δ条件下。定量分析车辆行驶时间的不确定性，以及给出在不确定性条件下最优路径的定义和数学表达式。假设各个路段所用时间服从正态分布N(μ，δ2)模型，利用MATLAB模拟（附录一）两条路径的正态分布图[1]：（图）图给每个路段设定一个预留到达的时间t，为了尽可能准确的到达目的地，选取95%的概率，满足P{T≤t}≥95%。那么最优路径的定义就是预留时间最小的那个路径。4.1.2模型建立与求解已知到达目的地所用时间和预留时间满足：P{T≤t}≥95%，将其转换为标准的正态分布：P{＜}≥95%，得到Φ0（）≥95%，≥Φ0-1（ρ），（其中ρ为95%）即可解得每条路段的t≥μ+1.645δ，取t=μ+1.645δ。我们将上诉模型称之为不确定条件下的NP模型。应用已知的数学表达式，将图1所给的数据：μ1=33，δ1=1；μ2=30，δ2=15；带入公式计算出：t（绕城）=34.645min，t（市区）=54.675min，最优路径为绕城快速路。4.2问题二4.2.1问题分析根据第一问中的最优路径定义和相关数学表达式，对于一般交通网络，可以结合K短路径算法建立NPK模型，最后从时间的渐进性态上分析其复杂性和收敛性。4.2.2模型建立与求解对于一般交通网络，为了方便设计算法找到最优路径，搜集了八个城市之间路段的时间均值和时间标准差，列表如下(表)：表我们可以将其转化为图论问题。将八个城市看做节点构成一个节点集：V={V1,V2,V3,V4,V5,V6,V7,V8}各个城市之间的道路看做边集：E={V1→V2，V1→V3，V1→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}则可将八个城市之间的交通网络图看做一个无向赋权图G（图.）每条路为图中的边。八个城市之间交通网络数据图图.在第一问定义的基础上，针对每条路段引入Bool系数βa,k,当该路段被选择时为1，否则为0。那么从起点到终点的路径可表示为a=1nβa,kΤkpath，可知其服从正态分布。通过求该路径的最小预留时间t(min)=E[Tk对于t（min）=E[Tkpath]+E[Tkpath]=a=1nβ由于Var[Tkpath]的根式不具有线性可加性。不能用经典的dijkstra算法求解。对此我们基于双目标规划的思路[3]，分别将E[Tkpath]，表表根据上述两表数据，运用公式:t（min）=E[Tkpath]+求出并集中的前16条最优路径的预留时间（表）表为了直观进行对比，将上表用excel制得如下柱状图：（图）（图）由图可知最优路径为V1→V3→V4→V8。(图)（图）收敛性分析：两个多项式时间算法之和还为多项式时间算法，其复杂性比列举的低。当问题的规模趋向无穷大时，时间复杂度的数量级将表现为渐进性态。即当路径K趋于无穷时，该模型一定收敛。4.3问题三4.3.1问题分析在问题三中，要求进一步考虑各路段之间行驶时间的相关性。我们用概率论中的协方差来表示这种耦合关系，并建立推广的NPK模型。4.3.2模型建立与求解在问题二中我们已得出最短预留时间的数学表达式：t（min）=[Tkpath为方便模型的建立与求解。在此我们假设仅相邻两条路段之间具有相关性。根据协方差的性质Var（T1+T2）=Var(T1)+Var(T2)+2cov(T1,T2);可以得出t=ETKpath称t为可靠时间。[4]以下图为例：图图为从A到B的一条路径。可靠时间t=E[Tkpath]+Φ-1（ρ）Var[由于问题二中，我们没有给出任意两条路段其时间随机向量（t1，t2）的密度函数。无法具体求出协方差cov[t1+t2]。对此我们假设任意两条相邻路段组成的协方差矩阵为一个随机的正半定矩阵。在matlab中随机函数rand()基础上得出如下协差阵（附录三）：表对问题二找出的预留时间最小的前16条路分别求其可靠时间（表）:表运用excel制图直观比较可靠时间大小（表）：表故得出最优路径为第三条V1->V3->V5->V8.（图）表五．模型评价本论文针对在不确定性条件下求解最优路径的问题，建立了以求最小预留时间t为目标的NP模型，并对问题一给出了合理的解答。在此基础上运用双目标规划的思想,结合求k短路径的方法，在没有考虑非相邻路段间相关性基础上，针对更一般的问题建立了推广的NPK模型。复杂性低，并随k增大具有较强收敛性。但在求均值与方差最短的并路径时，没有设计出相应的算法，且本文只针对k较大时有效。六．参考文献[1]袁东,肖广冰.详解matlab快速入门与应用.北京：电子工业出版社,2011,73-80.[2]邵虎,林兴强,孟强,谭美琳.基于出行问题可靠性的交通分配流问题.管理科学学报.2009.12(5):1-4.[3]百度百科,多目标规划，/,2015-8-11.[4]桂云丽.变分不等式的算法研究.西安电子科技大学硕士论文.2010:1-8.附录一.%二维正太分布图functionY=fun1(x);Y=(1/(2*pi*1))*exp(-(x-33)^2/(2*1*1));Y=(1/(2*pi*15))*exp(-(x-30)^2/(2*15*15));subplot(1,2,1);[xy]=meshgrid(25:0.1:40);z=1/(2*pi*1).*exp(-(x-33).^2/(2*1*1));h=mesh(x,y,z);set(h,'edgecolor','none','facecolor','interp');subplot(1,2,2);[xy]=meshgrid(-50:0.1:100);z=1/(2*pi*15).*exp(-(x-30).^2/(2*15*15));h=mesh(x,y,z);set(h,'edgecolor','none','facecolor','interp');附录二%第K短路算法function[shortestPaths,totalCosts]=kShortestPath(netCostMatrix,source,destination,k_paths)ifsource>size(netCostMatrix,1)||destination>size(netCostMatrix,1)warning('ThesourceordestinationnodearenotpartofnetCostMatrix');shortestPaths=[];totalCosts=[];elsek=1;[pathcost]=dijkstra(netCostMatrix,source,destination);%Pisacellarraythatholdsallthepathsfoundsofar:ifisempty(path)shortestPaths=[];totalCosts=[];elsepath_number=1;P{path_number,1}=path;P{path_number,2}=cost;current_P=path_number;%XisacellarrayofasubsetofP(usedbyYen'salgorithmbelow):size_X=1;X{size_X}={path_number;path;cost};%Spath_numberx1S(path_number)=path(1);%deviationvertexisthefirstnodeinitially%K=1istheshortestpathreturnedbydijkstra():shortestPaths{k}=path;totalCosts(k)=cost;while(k<k_paths&&size_X~=0)%removePfromXfori=1:length(X)ifX{i}{1}==current_Psize_X=size_X-1;X(i)=[];%deletecellbreak;endendP_=P{current_P,1};%P_iscurrentP,justtomakeiseasierforthenotations%Findwin(P_,w)insetS,wwasthedevvertexusedtofoundP_w=S(current_P);fori=1:length(P_)ifw==P_(i)w_index_in_path=i;endendforindex_dev_vertex=w_index_in_path:length(P_)-1%index_dev_vertexisindexinP_ofdeviationvertextemp_netCostMatrix=netCostMatrix;%RemoveverticesinPbeforeindex_dev_vertexandthereincidentedgesfori=1:index_dev_vertex-1v=P_(i);temp_netCostMatrix(v,:)=inf;temp_netCostMatrix(:,v)=inf;end%removeincidentedgeofvifvisinshortestPaths(K)UP_withsimilarsub_pathtoP_....SP_sameSubPath=[];index=1;SP_sameSubPath{index}=P_;fori=1:length(shortestPaths)iflength(shortestPaths{i})>=index_dev_vertexifP_(1:index_dev_vertex)==shortestPaths{i}(1:index_dev_vertex)index=index+1;SP_sameSubPath{index}=shortestPaths{i};endendendv_=P_(index_dev_vertex);orj=1:length(SP_sameSubPath)next=SP_sameSubPath{j}(index_dev_vertex+1);temp_netCostMatrix(v_,next)=inf;end%getthecostofthesubpathbeforedeviationvertexvsub_P=P_(1:index_dev_vertex);cost_sub_P=0;fori=1:length(sub_P)-1cost_sub_P=cost_sub_P+netCostMatrix(sub_P(i),sub_P(i+1));end%calldijkstrabetweendeviationvertextodestinationnode[dev_pc]=dijkstra(temp_netCostMatrix,P_(index_dev_vertex),destination);if~isempty(dev_p)path_number=path_number+1;P{path_number,1}=[sub_P(1:end-1)dev_p];%concatenatesubpath-to-vertex-to-destinationP{path_number,2}=cost_sub_P+c;S(path_number)=P_(index_dev_vertex);size_X=size_X+1;X{size_X}={path_number;P{path_number,1};P{path_number,2}};else%warning('k=%d,isempty(p)==true!\n',k);endend%Stepnecessaryotherwiseifkisbiggerthannumberofpossiblepaths%thelastresultswillgetrepeated!ifsize_X>0shortestXCost=X{1}{3};%costofpathshortestX=X{1}{1};%refnumberofpathfori=2:size_XifX{i}{3}<shortestXCostshortestX=X{i}{1};shortestXCost=X{i}{3};endendcurrent_P=shortestX;k=k+1;shortestPaths{k}=P{current_P,1};totalCosts(k)=P{current_P,2};else%k=k+1;endendendendfunction[shortestPath,totalCost]=dijkstra(netCostMatrix,s,d)n=size(netCostMatrix,1);fori=1:n%initializethefarthestnodetobeitself;farthestPrevHop(i)=i;%usedtocomputetheRTS/CTSrange;farthestNextHop(i)=i;end%allthenodesareun-visited;visited(1:n)=false;distance(1:n)=inf;%itstorestheshortestdistancebetweeneachnodeandthesourcenode;parent(1:n)=0;distance(s)=0;fori=1:(n-1),temp=[];forh=1:n,if~visited(h)%inthetree;temp=[tempdistance(h)];elsetemp=[tempinf];endend;[t,u]=min(temp);%itstartsfromnodewiththeshortestdistancetothesource;visited(u)=true;%markitasvisited;forv=1:n,%foreachneighborsofnodeu;if((netCostMatrix(u,v)+distance(u))<distance(v))distance(v)=distance(u)+netCostMatrix(u,v);%updatetheshortestdistancewhenashortershortestPathisfound;parent(v)=u;%updateitsparent;end;end;end;shortestPath=[];ifparent(d)~=0%ifthereisashortestPath!t=d;shortestPath=[d];whilet~=sp=parent(t);shortestPath=[pshortestPath];ifnetCostMatrix(t,farthestPrevHop(t))<netCostMatrix(t,p)farthestPrevHop(t)=p;end;ifnetCostMatrix(p,farthestNextHop(p))<netCostMatrix(p,t)farthestNextHop(p)=t;end;t=p;end;end;totalCost=distance(d);%return;function[]=TestKShortestPath(case_number)switch2case1netCostMatrix=[inf1020infinf4080inf;10inf2030infinf70inf;2020inf1050infinfinf;inf3010inf30infinf60;70inf5030inf30inf40;40infinfinf30inf2060;8070infinfinf20inf40;infinfinf60406040inf];source=1;destination=8;k=10;case2netCostMatrix=[inf425inf3622525inf;4inf10025infinf225inf;25100inf494infinfinf;inf2549inf144infinf100;36inf4144inf225inf4;225infinfinf225inf64100;25225infinfinf64inf25;infinfinf100410025inf];source=1;destination=8;k=10;otherwiseerror('Theonlycaseoptionsavailableare1,2or3');end%------Showcaseselected:------:fprintf('Youselectedcase#%d,a%dx%dnetwork:\n',1,size(netCostMatrix,1),size(netCostMatrix,1));disp(netCostMatrix);fprintf('Thepathrequestisfromsourcenode%dtodestinationnode%d,withK=%d\n',source,destination,k);%------CallkShortestPath------:[shortestPaths,totalCosts]=kShortestPath(netCostMatrix,source,destination,k);%------Displayresults------:fprintf('\nResultofFunctioncall:kShortestPath(netCostMatrix,source,destination,k)=\n\n');ifisempty(shortestPaths)fprintf('Nopathavailablebetweenthesenodes\n\n');elsefori=1:length(shortestPaths)fprintf('Path#%d:\n',i);disp(shortestPaths{i})fprintf('Costofpath%dis%5.2f\n\n',i,totalCosts(i));endendend附录三%矩阵的产生a=rand(17,17);fori=1:size(a,2)forj=1:size(a,2)c(i,j)=sum((a(:,i)-mean(a(:,i))).*(a(:,j)-mean(a(:,j))))/(size(a,1)-1);endendd=fix(c*100);d82200-20-50-20-10220-22100-1-122-3-30-20-223-1020600-3001-200110000-10510202000202100-101400-10-101