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

2015年暑期数学建模培训第一次模拟

承诺书

我们仔细阅读了数学建模联赛的竞赛规则。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与本队以外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其它公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们愿意承担由此引起的一切后果。

我们授权数学建模联赛赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号为（从A/B/C中选择一项填写）：

我们的参赛报名号为：

参赛组别（研究生或本科或专科）：本科

所属学校（请填写完整的全名）中国矿业大学南湖校区

参赛队员 (打印并签名) ：1. 赖增强

2. 兰卫旗

3. 李康杰

日期：2015年8月11日

获奖证书邮寄地址：中国矿业大学南湖校区桃4B5032邮政编码：221116

收件人姓名：赖增强 联系电话：183********

2015年暑期数学建模培训第一次模拟

编号专用页

竞赛评阅编号（由竞赛评委会评阅前进行编号）：

评阅记录

裁剪线裁剪线裁剪线

竞赛评阅编号（由竞赛评委会评阅前进行编号）：

参赛队伍的参赛号码：1（请各参赛队提前填写好）：

不确定条件下的最优路径问题

摘要

本文针对如何在复杂的交通环境下寻找一条可靠、快速、安全的最优路径的问题，考虑到交通堵塞、恶劣天气、路途成本等不确定因素对司机路径选择的影响，建立多个不确定条件下的最优路径模型。

对于问题一，我们在各个路段所用时间服从正态分布N(μ，δ2)的基础上,建立了在不确定条件下求最短路的NP模型，给每个路段设定一个预留到达的时间t，为了尽可能准确的到达目的地，选取95%的概率，满足P{T≤t}?95%,那么最优路径的定义就是预留时间最小的那个路径，将其转换为标准的正态分布，通过标准的正态分布得到了在不确定性条件下车辆从起点到终点预留时间的数学表达式：t=μ+δ。计算得对应的t（绕城）=34.645min，t（市区）=54.675min，那么最优路径为绕城快速路。

对于问题二，在第一问定义的基础上进一步引入Bool系数,在搜集得到的具体的交通网络中，建立了一个从起点到终点路径为的正态分布，通过求最小预留时间t(min)=E[]+ ,得出最优路径的算法。其中EE[]，=，但的根式不具有线性可加性。不能用经典的dijkstra算法求解。对此采用基于双目标规划的思路，利用第K短路径算法，分别对E[]，,运用matlab编程，找出各自前十条最短路径。之后在其并集中找出最优路径：V1→V3→V4→V8。由此建立了求最短路的NPK模型。最后从时间的渐进性态上分析模型的复杂性和收敛性。

对于问题三，我们只考虑各路段空间上的相关性，并用概率论中的协方差来表示这种耦合关系，建立了NPK模型。得出可靠时间的数学表达式t= E+；求解得最优路径：V1->V3->V5->V8。

关键词：NPK模型，K-短路，预留时间，正态分布，渐进性态。

一.问题重述

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

（一）针对一般的交通网络，假设已知每条路段行驶时间的均值和标准差，请建立相关的数学模型，定量地分析车辆行驶时间的不确定性，然后给出在不确定性条件下车辆从起点到终点的最优路径的定义和数学表达式。并将此模型应用到图1的例子中会选择哪条道路。

（二）根据第一问的定义，假设已知每条路段行驶时间的均值和标准差，设计算法搜索最优路径，并将该算法应用到具体的交通网络中，用计算结果验证算法的有效性。从理论上分析算法的收敛性、复杂性等性质。

（三）建立数学模型，描述下游路段发生交通拥堵使车辆减速或者排队，导致上游路段发生拥堵的交通路段之间驾驶时间的相关性，并将此相关性应用到第一问和第二问的最优路径搜索问题中，并设计算法解决考虑相关性的最优路径搜索问题，给出算例验证算法的有效性。从理论上分析算法的收敛性、复杂性等问题。

二.基本假设符号说明

3.1基本假设

。

3.1.2假设仅相邻两条路段之间具有相关性。

3.1.3假设任意两条相邻路段组成的协方差矩阵为一个随机的正半定矩阵。

3.2符号说明

3.2.1 两相邻节点之间的路段

3.2.2 从起点到终点的路径

3.2.3 Bool系数。当路径a通过此路段k时为1

3.2.4 对随机变量T求方差运算

3.2.5 （ρ） 标准正太分布累计概率函数的逆函数

3.2.6 随机向量t1,t2协方差

四.模型的建立与求解

4.1问题一

分析

本问题是对于题设的交通网络，已知每条路段行驶时间的均值μ和标准差δ条件下。定量分析车辆行驶时间的不确定性，以及给出在不确定性条件下最优路径的定义和数学表达式。

假设各个路段所用时间服从正态分布N(μ，δ2)模型，利用MATLAB模拟（附录一）两条路径的正态分布图[1]：）

给每个路段设定一个预留到达的时间t，为了尽可能准确的到达目的地，选取95%的概率，满足P{T≤t}≥95%。那么最优路径的定义就是预留时间最小的那个路径。

建立与求解

已知到达目的地所用时间和预留时间满足：P{T≤t}≥95%，将其转换为标准的正态分布：P{＜ }≥95%，得到Φ0（）≥95%，≥Φ0-1（ρ），（其中）即可解得每条路段的t≥μ+1.645δ，取t=μ+1.645δ。我们将上诉模型称之为不确定条件下的NP模型。

应用已知的数学表达式，将图1所给的数据：

μ 1=33， δ1=1；

μ 2=30， δ2=15；

带入公式计算出：

t（绕城）=34.645min，t（市区）=54.675min，最优路径为绕城快速路。

4.2问题二

根据第一问中的最优路径定义和相关数学表达式，对于一般交通网络，可以结合K短路径算法建立NPK模型，最后从时间的渐进性态上分析其复杂性和收敛性。

与求解

对于一般交通网络，为了方便设计算法找到最优路径，

我们可以将其转化为图论问题。将八个城市看做节点构成一个节点集：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}

八个城市之间交通网络数据图

图

在第一问定义的基础上，针对每条路段引入Bool系数,当该路段被选择时为1，否则为0。那么从起点到终点的路径可表示为Τkpath，可知其服从正态分布。通过求该路径的最小预留时间t(min)=E[]+ ,得出最优路径。[2]

对于t（min）= ]+，其中：

]=] =Var

由于的根式不具有线性可加性。不能用经典的dijkstra算法求解。对此我们基于双目标规划的思路[3]，分别将]，运用matlab编程（附录二），求出前十条最短路径。

表

表

根据上述两表数据，运用公式:

t（min）=E[]+]

求出并集中的前16条最优路径的预留时间（表

为了直观进行对比，将上表用excel制得如下柱状图：

（图）

由图可知最优路径为V1→V3→V4→V8。

（图）

收敛性分析：两个多项式时间算法之和还为多项式时间算法，其复杂性比列举的低。当问题的规模趋向无穷大时，时间复杂度的数量级将表现为渐进性态。即当路径K趋于无穷时，该模型一定收敛。

4.3问题三

分析

在问题三中，要求进一步考虑各路段之间行驶时间的相关性。我们用概率论中的协方差来表示这种耦合关系，并建立推广的NPK模型。

与求解

在问题二中我们已得出最短预留时间的数学表达式：

t（min）= ]+]

为方便模型的建立与求解。在此我们假设仅相邻两条路段之间具有相关性。根据协方差的性质Var（T1+T2）=Var(T1)+Var(T2)+2cov(T1,T2);可以得出

t= E+]

称t为可靠时间。[4]

以下图为例：

图为从A到B的一条路径。可靠时间t = ]+，其中]=E[T1]+E[T2]，。

由于问题二中，我们没有给出任意两条路段其时间随机向量（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.

附录一.

%二维正太分布图

function Y=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);

[x y]=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);

[x y]=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)

if source > size(netCostMatrix,1) || destination > size(netCostMatrix,1)

warning('The source or destination node are not part of netCostMatrix');

shortestPaths=[];

totalCosts=[];

else

k=1;

[path cost] = dijkstra(netCostMatrix, source, destination);

%P is a cell array that holds all the paths found so far:

if isempty(path)

shortestPaths=[];

totalCosts=[];

else

path_number = 1;

P{path_number,1} = path; P{path_number,2} = cost;

current_P = path_number;

%X is a cell array of a subset of P (used by Yen's algorithm below):

size_X=1;

X{size_X} = {path_number; path; cost};

%S path_number x 1

S(path_number) = path(1); %deviation vertex is the first node initially

% K = 1 is the shortest path returned by dijkstra():

shortestPaths{k} = path ;

totalCosts(k) = cost;

while (k < k_paths && size_X ~= 0 )

%remove P from X

for i=1:length(X)

if X{i}{1} == current_P

size_X = size_X - 1;

X(i) = [];%delete cell

break;

end

end

P_ = P{current_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_P);

for i = 1: length(P_)

if w == P_(i)

w_index_in_path = i;

end

end

for 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 edges

for i = 1: index_dev_vertex-1

v = P_(i);

temp_netCostMatrix(v,:)=inf;

temp_netCostMatrix(:,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_sameSubPath{index}=P_;

for i = 1: length(shortestPaths)

if length(shortestPaths{i}) >= index_dev_vertex

if P_(1:index_dev_vertex) == shortestPaths{i}(1:index_dev_vertex)

index = index+1;

SP_sameSubPath{index}=shortestPaths{i};

end

end

end

v_ = P_(index_dev_vertex);

or j = 1: length(SP_sameSubPath)

next = SP_sameSubPath{j}(index_dev_vertex+1);

temp_netCostMatrix(v_,next)=inf;

end

%get the cost of the sub path before deviation vertex v

sub_P = P_(1:index_dev_vertex);

cost_sub_P=0;

for i = 1: length(sub_P)-1

cost_sub_P = cost_sub_P + netCostMatrix(sub_P(i),sub_P(i+1));

end

%call dijkstra between deviation vertex to destination node

[dev_p c] = 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] ; %concatenate sub path- to -vertex -to- destination

P{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);

end

end

%Step necessary otherwise if k is bigger than number of possible paths

%the last results will get repeated !

if size_X > 0

shortestXCost= X{1}{3}; %cost of path

shortestX= X{1}{1}; %ref number of path

for i = 2 : size_X

if X{i}{3} < shortestXCost

shortestX= X{i}{1};

shortestXCost= X{i}{3};

end

end

current_P = shortestX;

k = k+1;

shortestPaths{k} = P{current_P,1};

totalCosts(k) = P{current_P,2};

else

%k = k+1;

end

end

end

end

function [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-visited;

visited(1:n) = false;

distance(1:n) = inf; % it stores the shortest distance between 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)];

else

temp=[temp inf];

end

end;

[t, u] = min(temp); % it starts from node 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); % update the shortest distance when a shorter shortestPath is found;

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 ~= s

p = parent(t);

shortestPath = [p shortestPath];

if netCostMatrix(t, farthestPrevHop(t)) < netCostMatrix(t, p)

farthestPrevHop(t) = p;

end;

if netCostMatrix(p, farthestNextHop(p)) < netCostMatrix(p, t)

farthestNextHop(p) = t;

end;

t = p;

end;

end;

totalCost = distance(d);

%return;

function [] = TestKShortestPath(case_number)

switch 2

case 1

netCostMatrix = [

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;

70 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 inf];

source=1;

destination=8;

k = 10;

case 2

netCostMatrix = [

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=1;

destination=8;

k = 10;

otherwise

error('The only case options available are 1, 2 or 3');

end

%------Show case selected:------:

fprintf('Youselected case #%d, a %dx%d network:\n',1,size(netCostMatrix,1),size(netCostMatrix,1));

disp(netCostMatrix);

fprintf('The path request is from source node %d to destination node %d, with K = %d \n',source,destination, k);

%------Call kShortestPath------:

[shortestPaths, totalCosts] = kShortestPath(netCostMatrix, source, destination, k);

%------Display results------:

fprintf('\nResult of Function call: kShortestPath(netCostMatrix, source, destination, k) = \n\n');

if isempty(shortestPaths)

fprintf('No path available between these nodes\n\n');

else

for i = 1: length(shortestPaths)

fprintf('Path # %d:\n',i);

disp(shortestPaths{i})

fprintf('Cost of path %d is %5.2f\n\n',i,totalCosts(i));

end

end

end

附录三

%矩阵的产生

a =rand(17,17);

for i=1:size(a,2)

for j=1:size(a,2)

c(i,j)=sum((a(:,i)-mean(a(:,i))).*(a(:,j)-mean(a(:,j))))/(size(a,1)-1);

end

end

d=fix(c*100);

d

8 2 2 0 0 -2 0 -5 0 -2 0 -1 0 2 2 0 -2

2 10 0 -1 -1 2 2 -3 -3 0 -2 0 -2 2 3 -1 0

2 0 6 0 0 -3 0 0 1 -2 0 0 1 1 0 0 0

0 -1 0 5 1 0 2 0 2 0 0 0 2 0 2 1 0

0 -1 0 1 4 0 0 -1 0 -1 0 1 0 0 -1 0 0

-2 2 -3 0 0 10 0 0 2 2 -4 0 -3 2 1 0 1

0 2 0 2 0 0 4 0 1 0 0 -1 3 0 1 0 0

-5 -3 0 0 -1 0 0 10 0 2 2 -1 0 -4 0 -2 2

0 -3 1 2 0 2 1 0 7 1 -1 -2 0 1 0 2 0

-2 0 -2 0 -1 2 0 2 1 7 0 -2 -3 0 2 2 3

0 -2 0 0 0 -4 0 2 -1 0 6 2 2 -3 -1 -1 0

-1 0 0 0 1 0 -1 -1 -2 -2 2 9 1 1 -5 -2 -2

0 -2 1 2 0 -3 3 0 0 -3 2 1 11 -1 0 -1 -3

2 2 1 0 0 2 0 -4 1 0 -3 1 -1 5 0 0 0

2 3 0 2 -1 1 1 0 0 2 -1 -5 0 0 9 2 0

0 -1 0 1 0 0 0 -2 2 2 -1 -2 -1 0 2 5 0

-2 0 0 0 0 1 0 2 0 3 0 -2 -3 0 0 0 4