CN108647835B - 基于设计速度的城市离散交通网络设计r语言实现方法 - Google Patents
基于设计速度的城市离散交通网络设计r语言实现方法 Download PDFInfo
- Publication number
- CN108647835B CN108647835B CN201810527899.0A CN201810527899A CN108647835B CN 108647835 B CN108647835 B CN 108647835B CN 201810527899 A CN201810527899 A CN 201810527899A CN 108647835 B CN108647835 B CN 108647835B
- Authority
- CN
- China
- Prior art keywords
- road
- traffic
- speed
- network
- design
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000013461 design Methods 0.000 title claims abstract description 55
- 238000000034 method Methods 0.000 title claims abstract description 19
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 42
- 238000005457 optimization Methods 0.000 claims abstract description 10
- 230000006870 function Effects 0.000 claims description 38
- 238000010276 construction Methods 0.000 claims description 17
- 238000011156 evaluation Methods 0.000 abstract description 17
- 206010039203 Road traffic accident Diseases 0.000 abstract description 5
- 238000003012 network analysis Methods 0.000 abstract description 3
- 239000011159 matrix material Substances 0.000 description 33
- 238000004364 calculation method Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 238000011160 research Methods 0.000 description 3
- ICAKYBFVLALTJC-LQJZCPKCSA-N (2s)-2-[4-(2-hydroxyethyl)triazol-1-yl]-1-[4-[4-[4-[(2s)-2-[4-(2-hydroxyethyl)triazol-1-yl]-3-methylbutanoyl]piperazin-1-yl]-6-[2-[2-(2-prop-2-ynoxyethoxy)ethoxy]ethylamino]-1,3,5-triazin-2-yl]piperazin-1-yl]-3-methylbutan-1-one Chemical compound N1([C@@H](C(C)C)C(=O)N2CCN(CC2)C=2N=C(NCCOCCOCCOCC#C)N=C(N=2)N2CCN(CC2)C(=O)[C@H](C(C)C)N2N=NC(CCO)=C2)C=C(CCO)N=N1 ICAKYBFVLALTJC-LQJZCPKCSA-N 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000006399 behavior Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 101000713575 Homo sapiens Tubulin beta-3 chain Proteins 0.000 description 1
- 206010027146 Melanoderma Diseases 0.000 description 1
- 102100036790 Tubulin beta-3 chain Human genes 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003912 environmental pollution Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000001172 regenerating effect Effects 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 238000002922 simulated annealing Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
- G06Q10/047—Optimisation of routes or paths, e.g. travelling salesman problem
-
- G06Q50/40—
Abstract
本发明针对基于设计速度的城市离散交通网络模型，提供了一种求解算法并使用R语言程序来实现。考虑到设计速度为路段的最佳行驶速度，偏离设计速度容易发生交通事故，本发明采用行驶速度与设计速度的均方误差作为网络安全水平的替代指标，将主动安全评价提前到了道路网络规划阶段。主要步骤包括：(1)建立城市离散交通网络设计的双层规划模型，其中上层规划以网络安全为政策目标，下层规划以用户均衡为政策的行为反应；(2)采用迭代优化算法进行模型求解；(3)设计R语言程序来实现求解方法；(4)结合交通网络分析中常用的Nguyen‑Dupuis网络描述了具体实施方式和该方法的有效性。使用R语言来实现模型的求解具有开源免费、可操作性强等技术优势。
Description
技术领域：
本发明针对基于设计速度的城市离散交通网络设计问题，建立了数学模型，设计了求解算法，采用R语言实现了一种可操作的计算机软件，属于交通工程技术领域。
背景技术：
传统的交通安全研究是在道路系统建成后，通过对大量事故数据的统计分析识别危险路段和事故黑点，然后提出相应的改善措施。由于此时可以采取的手段非常有限，交通安全治理的效果往往并不显著，并不能从根本上提高道路系统的安全水平。目前国内外对道路安全评价方法研究较多，常用的典型定量安全评价方法主要有：设计一致性评价、事故预测模型、交通冲突等。其中设计一致性评价属于事前主动安全评价，而事故预测模型和交通冲突属于事后被动安全评价。设计一致性评价指标又通常分为三类：基于运行速度协调性的评价指标、基于道路几何线形的评价指标和基于驾驶员工作负荷的评价指标。当前国内外道路安全评价分析时，普遍采用运行速度协调性来检验道路几何线形设计一致性。基于速度协调性的交通安全评价就是以速度为媒介，认为速度及其波动性与交通事故密切相关，把速度当作交通安全的替代指标，通过观测路段的运行速度进行交通安全评价。一般来说，运行速度协调性有以下两个含义：①同一路段实际行驶速度与设计速度相一致；②相邻路段实际行驶速度相一致。
本发明使用R语言实现模型求解，总体来说具有以下技术特点：1)开源免费：目前国外的交通规划类软件普遍价格高昂，而R语言是开源免费的；2)具有可操作性：R语言简单易用，适于一般的交通工程师和规划师学习采用；3)主动性：主动安全是交通安全规划的本质特征，本程序将交通安全评价提前到道路网络规划阶段；4)系统性：本程序考虑整个交通系统的安全性，而非单一路段或事故黑点。
发明内容：
技术问题：现阶段，交通规划软件基本都是国外产品，价格高昂并且操作复杂不易使用，特别是在交通规划中缺乏对于设计速度的考虑。设计速度是道路的最佳行驶速度，偏离设计速度行驶就容易发生交通事故。本发明提供一种基于设计速度的城市离散交通网络设计方法，建立城市离散交通网络模型并且设计求解模型的R语言方法。目的是将交通安全评价提前到道路网络规划阶段，采用同一路段实际行驶速度与设计速度的一致性作为道路网络安全评价的替代指标，以从源头上减少交通事故的发生。
技术方案：本发明提供了一种基于设计速度的城市离散交通网络设计的R语言实现方法，主要包括以下步骤：
步骤1：建立城市离散交通网络设计模型。模型上层为行驶速度与设计速度的均方误差最小化，下层为用户平衡模型。上层决策变量为ya，表示是否修建某条候选路段a，为0-1变量，a∈A，所有的候选路段构成0-1决策向量y。上层决定新建道路方案后，下层形成平衡状态网络流xa，也就是说路段流量xa是决策向量y的函数，表示为xa(y)。另外，道路网的规划受到资本的约束。假设单位长度的路段修建成本为ua，则长度为la的路段修建成本为uala。因此，双层规划问题表示为：
其中n为网络上路段的数目；D(y)为上层规划的目标函数；Z(x)为下层规划的目标函数；va d为路段a的设计速度；A为候选建设的路段集合；B为新建道路的资金约束；xa为路段a上的交通流量；
为自由流行驶时间，即路段a为空净状态时车辆自由行驶所需要的时间；ca为路段a的通行能力，即单位时间内路段实际可通过的车辆数；ta(xa，ca)为路段a以交通流量为自变量的阻抗函数，也称为行驶时间函数；
为出发地为r目的地为s的OD间的第k条路径上的流量；
为路段-路径相关变量，即0-1变量，如果路段a属于从出发地为r目的地为s的OD间的第k条路径，则
否则
qrs为出发地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″，″Roadcapacity″，″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(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)]＝0
V
#OD对42的最短路径和流量
sp42＝as.vector(b42)#转化为路段标号(Road)
x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点
colnames(x42)＝c(″Road″，″V42″)
x42
V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#OD对43的最短路径和流量
sp43＝as.vector(b43)#转化为路段标号(Road)
x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点
colnames(x43)＝c(″Road″，″V43″)
x43
V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#当所有最短路径上的流量求和，得到初始流量
VS＝rowSums(V[，seq(ncol(V)-3，ncol(V))])
VS
#步骤2：更新各路段的阻抗
t0＝e[，4]
c＝e[，5]
a＝0.15
b＝4
tp＝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)]＝0
V
#OD对42的最短路径和流量
sp42＝as.vector(b42)#转化为路段标号(Road)
x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点
colnames(x42)＝c(″Road″，″V42″)
x42
V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#OD对43的最短路径和流量
sp43＝as.vector(b43)#转化为路段标号(Road)
x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点
colnames(x43)＝c(″Road″，″V43″)
x43
V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#当所有最短路径上的流量求和，得到迭代方向
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)break
VS＝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)
plan
nv＝numeric(n)
nv
for(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)]#输出均方误差最小的建设方案
有益效果：本软件采用开源免费的R语言，并且方便交通工程师使用，特别是本发明建立了城市离散交通网络设计的双层规划模型，以实际行驶速度与设计速度的均方误差最小为上层政策目标，以用户均衡作为下层的政策反应。并且针对该模型提供了一种求解算法并使用R语言程序来实现。最后，采用交通网络分析中常用的Nguyen-Dupuis网络验证了该软件的有效性。
附图说明：
图1是迭代优化算法的流程图。
图2是Nguyen-Dupuis测试网络，其中虚线段为候选建设路段。
图3是不同方案下道路网络的速度均方误差。
具体实施方式：
下面结合说明书附图，对本发明作进一步说明：
步骤1：城市离散交通网络设计模型
交通网络设计问题是一个考虑使用者路径选择行为，并且在给定的各种约束条件下，选择改建或者新建的路段，从而使得交通网络某种性能最优化的问题。这个问题是交通规划领域研究的重点和难点问题，十分具有挑战，引起了许多学者的研究兴趣。交通网络设计问题属于典型的领导者-跟随者的双层规划问题，其上层问题为交通网络的某种性能最优化，下层问题通常为用户平衡问题。常见的网络性能包括总行驶时间、环境污染、投资费用，备用能力、可靠性等。根据政策变量的不同，交通网络设计问题又分为离散交通网络设计、连续交通网络设计问题和混合交通网络设计问题。其中，离散交通网络设计一般指在投入资金预算的情况下，采用定量方法研究在现有路网上新建某些路段以使系统达到最优的问题，属于交通规划的方案设计部分。
虽然交通安全广受关注，但以网络安全性能为目标的交通网络设计问题尚未见到。本发明采用城市交通网络中各路段实际行驶速度与相应设计速度的均方误差作为网络安全水平的替代指标，该均方误差越小网络安全水平越高，此时驾驶员面临最适宜的道路行驶环境，有利于减少交通事故的发生。因此，交通网络规划的政策目标就是路段行驶速度与设计速度的均方误差最小化，即
其中n为交通网络路段的数目，
为路段a的设计速度，va为路段a的行驶速度，va＝la/ta，la为路段a的长度，ta为路段a的行驶时间。考虑到拥挤效应，路段行驶时间为交通流量的函数ta(xa，ca)，其中xa为路段交通流量。目前交通研究和实践中广泛应用由美国公路局(Bureau of Public Road，BPR)开发的行驶时间函数，被称为BPR函数，形式为：
其中α和β为阻滞系数，在美国公路局交通流分配程序中，取值分别为α＝0.15和β＝4，也可以由实际数据回归分析求得，ca为路段a的通行能力。因此，式(1)转化为：
不同于事后评价式的被动交通安全，主动交通安全规划需要预测不同网络设计下路段的运行速度va。本发明采用Wardrop用户均衡原理作为网络用户对不同政策的行为反应，建立了一个双层规划模型用于城市离散交通网络设计，其上层为行驶速度与设计速度的均方误差最小化，下层为用户平衡模型。上层决策变量为ya，表示是否修建某条候选路段a，为0-1变量，a∈A，所有的候选路段构成0-1决策向量y。上层决定新建道路方案后，下层形成平衡状态网络流xa，也就是说路段流量xa是决策向量y的函数，表示为xa(y)。另外，道路网的规划受到资本的约束。假设单位长度的路段修建成本为ua，则长度为la的路段修建成本为uala。因此，双层规划问题表示为：
其中n为网络上路段的数目；D(y)为上层规划的目标函数；Z(x)为下层规划的目标函数；va d为路段a的设计速度；A为候选建设的路段集合；B为新建道路的资金约束；xa为路段a上的交通流量；
为自由流行驶时间，即路段a为空净状态时车辆自由行驶所需要的时间；ca为路段a的通行能力，即单位时间内路段可通过的车辆数；ta(xa，ca)为路段a以交通流量为自变量的阻抗函数，也称为行驶时间函数；
为出发地为r目的地为s的OD间的第k条路径上的流量；
为路段-路径相关变量，即0-1变量，如果路段a属于从出发地为r目的地为s的OD间的第k条路径，则
否则
qrs为出发地r和目的地s之间的OD交通需求量。
步骤2：迭代优化算法
城市交通网络设计的双层规划模型是一个NP-hard的问题，是一个具有十分挑战的问题。学者们提出了很多方法，主要有迭代优化算法(IOA)、基于(近似)梯度的算法、遗传算法、模拟退火法等。本发明采用迭代优化算法，其上层采用枚举法，下层采用Frank-Wolfe算法。算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案。详细的迭代优化算法可以归纳如下：
Step 1：生成一个建设方案。判断其是否满足上层约束，如果不满足再生成一个新的建设方案，直到得到一个可行的建设方案y，令m＝1。
Step 2：初始化：按照
进行一次0-1交通流分配，得到各路段的流量
令n＝1。
Step 3：更新各路段的阻抗：
Step 4：寻找下一步的迭代方向：按照
再进行一次0-1分配，得到一组附加交通量流量
Step 5：确定迭代步长：求满足下式的λ：
Step 6：确定新的迭代起点：
Step 7：收敛性检验：如果
满足规定的收敛准则，如
其中ε是预先给定的误差限值，则
即为要求的平衡解，停止计算，否则令n＝n+1，返回Step1。
Step 8：计算行驶速度和设计速度的均方误差：路段的行驶时间为
则容易计算路段的行驶速度为va＝la/ta，进而计算该行驶速度与设计速度的均方误差Dm。转入Step 1，计算下一个可行方案。当所有可行方案计算完毕后，停止计算，从所有可行方案中寻找上层目标函数最优的建设方案。迭代优化算法的流程图如图1所示。
步骤3：R语言求解算法
#步骤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″，″Roadcapacity″，″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(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)]＝0
V
#OD对42的最短路径和流量
sp42＝as.vector(b42)#转化为路段标号(Road)
x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点
colnames(x42)＝c(″Road″，″V42″)
x42
V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#OD对43的最短路径和流量
sp43＝as.vector(b43)#转化为路段标号(Road)
x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点
colnames(x43)＝c(″Road″，″V43″)
x43
V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#当所有最短路径上的流量求和，得到初始流量
VS＝rowSums(V[，seq(ncol(V)-3，ncol(V))])
VS
#步骤2：更新各路段的阻抗
t0＝e[，4]
c＝e[，5]
a＝0.15
b＝4
tp＝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)]＝0
V
#OD对42的最短路径和流量
sp42＝as.vector(b42)#转化为路段标号(Road)
x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点
colnames(x42)＝c(″Road″，″V42″)
x42
V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#OD对43的最短路径和流量
sp43＝as.vector(b43)#转化为路段标号(Road)
x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点
colnames(x43)＝c(″Road″，″V43″)
x43
V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵
V[is.na(V)]＝0
V
#当所有最短路径上的流量求和，得到迭代方向
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)break
VS＝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)
plan
nv＝numeric(n)
nv
for(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)]#输出均方误差最小的建设方案
步骤4：算例分析
本算法采用交通网络分析中常用的Nguyen-Dupuis网络(图2)作为测试网络。道路参数列在表1中，其中路段编号20-24为候选的路段。为简化计算，本例假设政策制定者要求从5条候选道路中选择2条进行开工建设，原有各路段的设计速度都为30km/h，新建各路段的设计速度都为40km/h。交通需求矩阵如表2所示。假设候选路段的结合为A，则该算例的上层目标函数简化为：
表1 Nguyen-Dupuis网络的路段参数
表2 Nguyen-Dupuis网络的OD对交通需求
对所有可行方案进行交通流平衡分配后，可以得到该方案下各路段的行驶速度，进而可得行驶速度与设计速度的均方误差。从表3和图3可以看出，不同方案下，均方误差相差较大，最大为修建路径23和24，速度均方误差为48.77，最小为修建路径20与24，速度均方误差为28.96。因此，为了设计最安全的交通网络，应该修建路径20与24。
表3不同方案下道路网络的速度均方误差
Claims (1)
1.一种用于城市离散交通网络设计的方法，采用R语言实现，步骤如下：
步骤1：建立城市离散交通网络设计模型，上层为行驶速度与设计速度的均方误差最小化，下层为用户平衡模型，上层决策变量为ya，表示是否修建某条候选路段a，为0-1变量，所有的候选路段构成0-1决策向量y，上层决定新建道路方案后，下层形成平衡状态网络流xa，也就是说路段流量xa是决策向量y的函数，表示为xa(y)，另外，道路网的规划受到资本的约束，假设单位长度的路段修建成本为ua，则长度为la的路段修建成本为uala，因此，双层规划问题表示为：
其中va d为路段a的设计速度；D(y)为是上层规划的目标函数，表示路段行驶速度与设计速度的均方误差；A为候选建设的路段集合；n为所有路段的数目；B为新建道路的资金约束；xa为路段a上的交通流量；
为自由流行驶时间，即路段a为空净状态时车辆自由行驶所需要的时间；ca为路段a的通行能力，即单位时间内路段可通过的车辆数；ta(xa，ca)为路段a以交通流量为自变量的阻抗函数，也称为行驶时间函数；
为出发地为r目的地为s的OD间的第k条路径上的流量；
为路段-路径相关变量，即0-1变量，如果路段a属于从出发地为r目的地为s的OD间的第k条路径，则
否则
qrs为出发地r和目的地s之间的OD交通需求量；Z(x)为下层规划的目标函数，没有直观的物理解释，借助其来求解平衡分配问题；
步骤2：使用迭代优化算法进行求解，上层采用枚举法，下层采用Frank-Wolfe算法，算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810527899.0A CN108647835B (zh) | 2018-05-18 | 2018-05-18 | 基于设计速度的城市离散交通网络设计r语言实现方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810527899.0A CN108647835B (zh) | 2018-05-18 | 2018-05-18 | 基于设计速度的城市离散交通网络设计r语言实现方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108647835A CN108647835A (zh) | 2018-10-12 |
| CN108647835B true CN108647835B (zh) | 2021-10-08 |
Family
ID=63758264
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810527899.0A Expired - Fee Related CN108647835B (zh) | 2018-05-18 | 2018-05-18 | 基于设计速度的城市离散交通网络设计r语言实现方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108647835B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110457012A (zh) * | 2019-06-27 | 2019-11-15 | 东南大学 | 用于可持续交通网络设计的多属性决策软件 |
| CN116777120B (zh) * | 2023-08-16 | 2023-10-27 | 合肥综合性国家科学中心人工智能研究院(安徽省人工智能实验室) | 一种基于路网od对的城市道路移动源碳排放计算方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH11345396A (ja) * | 1998-06-02 | 1999-12-14 | Pub Works Res Inst Ministry Of Constr | 走行支援道路システムの合流制御システムにおける本線交通流予測方法 |
| CN101593320A (zh) * | 2009-03-27 | 2009-12-02 | 北京交通大学 | 一种基于运输需求特征的多种交通方式组合运能优化方法 |
| CN104537439A (zh) * | 2014-12-29 | 2015-04-22 | 中南大学 | 一种用于缓解轨道交通拥挤的最小成本路径和混合路径分配方法 |
| CN107766945A (zh) * | 2017-09-06 | 2018-03-06 | 北京交通发展研究院 | 城市路网承载力的计算方法 |
-
2018
- 2018-05-18 CN CN201810527899.0A patent/CN108647835B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH11345396A (ja) * | 1998-06-02 | 1999-12-14 | Pub Works Res Inst Ministry Of Constr | 走行支援道路システムの合流制御システムにおける本線交通流予測方法 |
| CN101593320A (zh) * | 2009-03-27 | 2009-12-02 | 北京交通大学 | 一种基于运输需求特征的多种交通方式组合运能优化方法 |
| CN104537439A (zh) * | 2014-12-29 | 2015-04-22 | 中南大学 | 一种用于缓解轨道交通拥挤的最小成本路径和混合路径分配方法 |
| CN107766945A (zh) * | 2017-09-06 | 2018-03-06 | 北京交通发展研究院 | 城市路网承载力的计算方法 |
Non-Patent Citations (1)
| Title |
|---|
| 《交通网络设计的优化模型及算法》;桂岚;《系统工程》;20061231;第26-32页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108647835A (zh) | 2018-10-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108470444B (zh) | 一种基于遗传算法优化的城市区域交通大数据分析系统与方法 | |
| CN108776722B (zh) | 一种基于设计速度的城市离散交通网络设计方法 | |
| CN108694278B (zh) | 一种基于道路负载均衡的城市离散交通网络设计方法 | |
| CN107944605B (zh) | 一种基于数据预测的动态交通路径规划方法 | |
| CN108804801B (zh) | 基于目标配流的城市离散交通网络设计r语言实现方法 | |
| CN109959388B (zh) | 一种基于栅格扩展模型的智能交通精细化路径规划方法 | |
| CN108681788B (zh) | 一种基于主动安全的城市离散交通网络设计方法 | |
| CN108647475B (zh) | 基于负载均衡的城市离散交通网络设计r语言实现方法 | |
| Zheng et al. | A stochastic simulation-based optimization method for equitable and efficient network-wide signal timing under uncertainties | |
| CN112629533B (zh) | 基于路网栅格化道路车流预测的精细化路径规划方法 | |
| CN108647835B (zh) | 基于设计速度的城市离散交通网络设计r语言实现方法 | |
| Ma et al. | Road maintenance optimization model based on dynamic programming in urban traffic network | |
| Zhang et al. | Integrated path controlling and subsidy scheme for mobility and environmental management in automated transportation networks | |
| Karimi et al. | Integration of dynamic route guidance and freeway ramp metering using model predictive control | |
| Zhang et al. | Bi-objective bi-level optimization for integrating lane-level closure and reversal in redesigning transportation networks | |
| Mehr et al. | A game theoretic macroscopic model of lane choices at traffic diverges with applications to mixed–autonomy networks | |
| CN110489871A (zh) | 考虑新能源汽车渗透的环境影响评估软件 | |
| Karimi et al. | Integrated model predictive control of dynamic route guidance information systems and ramp metering | |
| CN108776584B (zh) | 城市离散交通网络主动安全设计的r语言实现方法 | |
| Chen et al. | Network-based optimization modeling of manhole setting for pipeline transportation | |
| CN114842641B (zh) | 面向省域的多模式链式交通分配方法 | |
| CN114495529B (zh) | 一种基于分布式模型预测控制的信号配时优化系统 | |
| Azadi et al. | Estimation of Delay and Surrogate Conflicts in Automated Network Management Using Linear Regression and Multigene Genetic Programming | |
| Kolarich et al. | Stackelberg routing of autonomous cars in mixed-autonomy traffic networks | |
| Li et al. | An extended game–theoretic model for aggregate lane choice behavior of vehicles at traffic diverges with a bifurcating lane |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20211008 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |