[发明专利]基于设计速度的城市离散交通网络设计R语言实现方法

摘要

本发明针对基于设计速度的城市离散交通网络模型，提供了一种求解算法并使用R语言程序来实现。考虑到设计速度为路段的最佳行驶速度，偏离设计速度容易发生交通事故，本发明采用行驶速度与设计速度的均方误差作为网络安全水平的替代指标，将主动安全评价提前到了道路网络规划阶段。主要步骤包括：(1)建立城市离散交通网络设计的双层规划模型，其中上层规划以网络安全为政策目标，下层规划以用户均衡为政策的行为反应；(2)采用迭代优化算法进行模型求解；(3)设计R语言程序来实现求解方法；(4)结合交通网络分析中常用的Nguyen‑Dupuis网络描述了具体实施方式和该方法的有效性。使用R语言来实现模型的求解具有开源免费、可操作性强等技术优势。

主权项

1.本发明针对基于设计速度的城市离散交通网络设计模型，设计了一种R语言求解方法，涉及步骤如下：步骤1：建立城市离散交通网络设计模型，上层为行驶速度与设计速度的均方误差最小化，下层为用户平衡模型，上层决策变量为ya，表示是否修建某条候选路段a，为0‑1变量，a∈A，所有的候选路段构成0‑1决策向量y，上层决定新建道路方案后，下层形成平衡状态网络流xa，也就是说路段流量xa是决策向量y的函数，表示为xa(y)，另外，道路网的规划受到资本的约束，假设单位长度的路段修建成本为ua，则长度为la的路段修建成本为uala，因此，双层规划问题表示为：其中A为候选建设的路段集合；B为新建道路的资金约束；x_a为路段a上的交通流量；为自由流行驶时间，即路段a为空净状态时车辆自由行驶所需要的时间；c_a为路段a的通行能力，即单位时间内路段可通过的车辆数；t_a(x_a，c_a)为路段a以交通流量为自变量的阻抗函数，也称为行驶时间函数；为出发地为r目的地为s的OD间的第k条路径上的流量；为路段‑路径相关变量，即0‑1变量，如果路段a属于从出发地为r目的地为s的OD间的第k条路径，则否则q_rs为出发地r和目的地s之间的OD交通需求量；步骤2：使用迭代优化算法进行求解，上层采用枚举法，下层采用Frank‑Wolfe算法，算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案；步骤3：设计求解的具体程序如下：#步骤1：初始化，按格式输入数据、函数和必要的包#1.1加载计算最短路径的包，准备调用dijkstra最短路径算法，注意igraph包首次需要安装，然后才能调用：#install.packages(″igraph″)library(igraph)options(digits＝3)#1.2创建图的距离矩阵，包含所有的候选路段：第一列为路段标号(Road)，第二列为路段起点标号(Road origin)，第三列为路段终点标号(Road destination)，第四列为该路段自由流时间(free flow time)，第五列为道路通行能力(capacity)，第六列为道路长度(length)；此处以交通配流中常用的Nguyen‑Dupuis网络为例，详细的参数设置可参考程序文档：#也可以在Excel中复制，然后执行#e＝read.delim(″clipboard″，header＝F)e＝matrix(c(1，1，5，7.0，800，4.00，2，1，12，9.0，400，6.00，3，4，5，9.0，200，5.00，4，4，9，12.0，800，8.00，5，5，6，3.0，350，2.00，6，5，9，9.0，400，5.00，7，6，7，5.0，800，3.00，8，6，10，13.0，250，8.00，9，7，8，5.0，250，3.00，10，7，11，9.0，300，6.00，11，8，2，9.0，550，5.00，12，9，10，10.0，550，6.00，13，9，13，9.0，600，5.00，14，10，11，6.0，700，4.00，15，11，2，9.0，500，6.00，16，11，3，8.0，300，5.00，17，12，6，7.0，200，4.00，18，12，8，14.0，400，6.00，19，13，3，11.0，600，7.00，20，1，6，10.0，600，8.00，21，5，10，10.0，600，8.00，22，6，11，10.0，600，8.00，23，7，2，10.0，600，8.00，24，10，3，10.0，600，8.00)，ncol＝6，byrow＝T)e＝cbind(e，cbind(c(rep(30，19)，rep(40，5))))#添加各个路段的设计速度colnames(e)＝c(″Road″，″Road origin″，″Road destination″，″Time″，″Road capacity″，″Road length″，″Designed speed″)e#1.3输入交通需求矩阵，第一列为起讫点对的标号(OD pair)，第二列为起点标号(origin)，第三列为终点标号(destination)，第四列为交通需求(demand)#也可以在Excel中复制，然后执行#d＝read.delim(″clipboard″)d＝matrix(c(1，1，2，400，2，1，3，800，3，4，2，600，4，4，3，200)，ncol＝4，byrow＝T)colnames(d)＝c(″OD pair″，″Origin″，″Destination″，″Demand″)d#自定的Frank‑Wolfe算法函数fw＝function(e，d){#1.4根据路径自由流时间计算各个OD对的最短路径和路径流量g＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝e[，4])#创建图，13为节点的个数b12＝get.shortest.paths(g，from＝″1″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点2的最短路径b13＝get.shortest.paths(g，from＝″1″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点3的最短路径b42＝get.shortest.paths(g，from＝″4″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点2的最短路径b43＝get.shortest.paths(g，from＝″4″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点3的最短路径#创建一个临时矩阵，用于保存各个OD对的最短路径和流量V＝cbind(e[，1])colnames(V)＝″Road″V#OD对12的最短路径和流量sp12＝as.vector(b12)#转化为路段标号(Road)x12＝cbind(e[sp12，1]，rep(d[1，4]，length(sp12)))#路段标号和流量，算法中的迭代起点colnames(x 12)＝c(″Road″，″V12″)x12V＝merge(V，x12，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#OD对13的最短路径和流量sp13＝as.vector(b13)#转化为路段标号(Road)x13＝cbind(e[sp13，1]，rep(d[2，4]，length(sp13)))#路段标号和流量，算法中的迭代起点colnames(x13)＝c(″Road″，″V13″)x13V＝merge(V，x13，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#OD对42的最短路径和流量sp42＝as.vector(b42)#转化为路段标号(Road)x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点colnames(x42)＝c(″Road″，″V42″)x42V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#OD对43的最短路径和流量sp43＝as.vector(b43)#转化为路段标号(Road)x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点colnames(x43)＝c(″Road″，″V43″)x43V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#当所有最短路径上的流量求和，得到初始流量VS＝rowSums(V[，seq(ncol(V)‑3，ncol(V))])VS#步骤2：更新各路段的阻抗t0＝e[，4]c＝e[，5]a＝0.15b＝4tp＝function(v){  t0*(1+a*(v/c)^b)}repeat{  #步骤3：寻找下一个迭代方向  g2＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝tp(VS))#构造图，13为节点的个数，更新路段阻抗  b12＝get.shortest.paths(g2，from＝″1″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点2的最短路径  b13＝get.shortest.paths(g2，from＝″1″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点3的最短路径  b42＝get.shortest.paths(g2，from＝″4″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点2的最短路径  b43＝get.shortest.paths(g2，from＝″4″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点3的最短路径  #创建一个临时矩阵，用于保存各个OD对的最短路径和流量  V＝cbind(e[，1])  colnames(V)＝″Road″  V  #OD对12的最短路径和流量  sp12＝as.vector(b12)#转化为路段标号(Road)  x12＝cbind(e[sp12，1]，rep(d[1，4]，length(sp12)))#路段标号和流量，算法中的迭代起点  colnames(x12)＝c(″Road″，″V12″)  x12  V＝merge(V，x12，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵  V[is.na(V)]＝0  V  #OD对13的最短路径和流量  sp13＝as.vector(b13)#转化为路段标号(Road)  x13＝cbind(e[sp13，1]，rep(d[2，4]，length(sp13)))#路段标号和流量，算法中的迭代起点  colnames(x13)＝c(″Road″，″V13″)  x13  V＝merge(V，x13，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#OD对42的最短路径和流量sp42＝as.vector(b42)#转化为路段标号(Road)x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点colnames(x42)＝c(″Road″，″V42″)x42V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#OD对43的最短路径和流量sp43＝as.vector(b43)#转化为路段标号(Road)x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点colnames(x43)＝c(″Road″，″V43″)x43V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵V[is.na(V)]＝0V#当所有最短路径上的流量求和，得到迭代方向VS2＝rowSums(V[，seq(ncol(V)‑3，ncol(V))])VS2#步骤4：计算迭代步长step＝function(lamda){  x2＝VS2  x1＝VS  q＝x1+lamda*(x2‑x1)  sum((x2‑x1)*tp(q))}lamda＝uniroot(step，c(0，1))$root#注意lamda的取值范围，步长不能太长lamda#步骤5：确定新的迭代起点VS3＝VS+lamda*(VS2‑VS)VS3#步骤6：收敛性检验if((sqrt(sum((VS3‑VS)^2))/sum(VS))＜1e‑5)breakVS＝VS3#如果不满足收敛条件则用新点VS3替代原点VS，如此循环直到收敛}#步骤7：输出平衡状态各路径的流量、通行时间和速度result＝cbind(e[，1]，round(VS，0)，tp(VS)，e[，6]/(tp(VS)/60))colnames(result)＝c(″Road″，″Volume″，″Time″，″Speed″)result}#步骤8：对可行方案，调用自定义的fw()函数，即Frank‑Wolfe算法：输入为可行方案的图矩阵e和交通需求矩阵d；输出用户平衡状态时各路段的交通流量、通行时间和速度，并据此计算该方案下的网络行驶速度与设计速度的均方误差，对所有方案计算速度均方误差后，输出速度均方误差最小的方案，即为最优方案de＝c(20，21，22，23，24)n＝choose(length(de)，2)new＝combn(de，2)old＝matrix(rep(c(1：19)，each＝n)，byrow＝T，nrow＝19)#所有的方案plan＝rbind(old，new)plannv＝numeric(n)nvfor(i in 1：n){  con＝fw(e[plan[，i]，]，d)  nv[i]＝mean((con[con[，2]＞1，4]‑e[plan[，i]，7][con[，2]＞1])^2)#假设各路段的设计速度都为40km/h，注意对于流量小于1的路段不参与计算均方误差}nv#输出所有方案的均方误差new[，which.min(nv)]#输出均方误差最小的建设方案。