本发明针对基于设计速度的城市离散交通网络模型，提供了一种求解算法并使用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。因此，双层规划问题表示为：
其中a为候选建设的路段集合；b为新建道路的资金约束；xa为路段a上的交通流量；
步骤2：使用迭代优化算法进行求解。其上层采用枚举法，下层采用frank-wolfe算法。算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案。
步骤3：设计求解程序如下：
#步骤1：初始化。按格式输入数据、函数和必要的包。
#1.1加载计算最短路径的包，准备调用dijkstra最短路径算法，注意igraph包首次需要安装，然后才能调用。
#install.packages(″igraph″)
library(igraph)
options(digits＝3)
#1.2创建图的距离矩阵，包含所有的候选路段。第一列为路段标号(road)，第二列为路段起点标号(roadorigin)，第三列为路段终点标号(roaddestination)，第四列为该路段自由流时间(freeflowtime)，第五列为道路通行能力(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″，″roadorigin″，″roaddestination″，″time″，″roadcapacity″，″roadlength″，″designedspeed″)
e
#1.3输入交通需求矩阵，第一列为起讫点对的标号(odpair)，第二列为起点标号(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(″odpair″，″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.ha(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(1ength(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(iin1：n)
{
con＝fw(e[plan[，i]，]，d)
nv[i]＝mean((con[con[，2]＞1，4]-e[plan[，i]，7][con[，2]＞11)＾2)#假设各路段的设计速度都为40km/h，注意对于流量小于1的路段不参与计算均方误差
}
nv#输出所有方案的均方误差
new[，which.min(nv)]#输出均方误差最小的建设方案
有益效果：本软件采用开源免费的r语言，并且方便交通工程师使用，特别是本发明建立了城市离散交通网络设计的双层规划模型，以实际行驶速度与设计速度的均方误差最小为上层政策目标，以用户均衡作为下层的政策反应。并且针对该模型提供了一种求解算法并使用r语言程序来实现。最后，采用交通网络分析中常用的nguyen-dupuis网络验证了该软件的有效性。
附图说明：
图1是迭代优化算法的流程图。
图2是nguyen-dupuis测试网络，其中虚线段为候选建设路段。
图3是不同方案下道路网络的速度均方误差。
具体实施方式：
下面结合说明书附图，对本发明作进一步说明：
步骤1：城市离散交通网络设计模型
交通网络设计问题是一个考虑使用者路径选择行为，并且在给定的各种约束条件下，选择改建或者新建的路段，从而使得交通网络某种性能最优化的问题。这个问题是交通规划领域研究的重点和难点问题，十分具有挑战，引起了许多学者的研究兴趣。交通网络设计问题属于典型的领导者-跟随者的双层规划问题，其上层问题为交通网络的某种性能最优化，下层问题通常为用户平衡问题。常见的网络性能包括总行驶时间、环境污染、投资费用，备用能力、可靠性等。根据政策变量的不同，交通网络设计问题又分为离散交通网络设计、连续交通网络设计问题和混合交通网络设计问题。其中，离散交通网络设计一般指在投入资金预算的情况下，采用定量方法研究在现有路网上新建某些路段以使系统达到最优的问题，属于交通规划的方案设计部分。
虽然交通安全广受关注，但以网络安全性能为目标的交通网络设计问题尚未见到。本发明采用城市交通网络中各路段实际行驶速度与相应设计速度的均方误差作为网络安全水平的替代指标，该均方误差越小网络安全水平越高，此时驾驶员面临最适宜的道路行驶环境，有利于减少交通事故的发生。因此，交通网络规划的政策目标就是路段行驶速度与设计速度的均方误差最小化，即
其中n为交通网络路段的数目，
其中α和β为阻滞系数，在美国公路局交通流分配程序中，取值分别为α＝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。因此，双层规划问题表示为：
其中a为候选建设的路段集合；b为新建道路的资金约束：xa为路段a上的交通流量；
步骤2：迭代优化算法
城市交通网络设计的双层规划模型是一个np-hard的问题，是一个具有十分挑战的问题。学者们提出了很多方法，主要有迭代优化算法(ioa)、基于(近似)梯度的算法、遗传算法、模拟退火法等。本发明采用迭代优化算法，其上层采用枚举法，下层采用frank-wolfe算法。算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案。详细的迭代优化算法可以归纳如下：
step1：生成一个建设方案。判断其是否满足上层约束，如果不满足再生成一个新的建设方案，直到得到一个可行的建设方案v，令m＝1。
step2：初始化：按照
step3：更新各路段的阻抗：
step4：寻找下一步的迭代方向：按照
step5：确定迭代步长：求满足下式的λ：
step6：确定新的迭代起点：
step7：收敛性检验：如果
step8：计算行驶速度和设计速度的均方误差：路段的行驶时间为
步骤3：r语言求解算法
#步骤1：初始化。按格式输入数据、函数和必要的包。
#1.1加载计算最短路径的包，准备调用dijkstra最短路径算法，注意igraph包首次需要安装，然后才能调用。
#install.packages(″igraph″)
library(igraph)
options(digits＝3)
#1.2创建图的距离矩阵，包含所有的候选路段。第一列为路段标号(road)，第二列为路段起点标号(roadorigin)，第三列为路段终点标号(roaddestination)，第四列为该路段自由流时间(freeflowtime)，第五列为道路通行能力(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″，″roadorigin″，″roaddestination″，″time″，″roadcapacity″，″roadlength″，″designedspeed″)
e
#1.3输入交通需求矩阵，第一列为起讫点对的标号(odpair)，第二列为起点标号(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(″odpair″，″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(iin1：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，则该算例的上层目标函数简化为：
表1nguyen-dupuis网络的路段参数
表2nguyen-dupuis网络的od对交通需求
对所有可行方案进行交通流平衡分配后，可以得到该方案下各路段的行驶速度，进而可得行驶速度与设计速度的均方误差。从表3和图3可以看出，不同方案下，均方误差相差较大，最大为修建路径23和24，速度均方误差为48.77，最小为修建路径20与24，速度均方误差为28.96。因此，为了设计最安全的交通网络，应该修建路径20与24。
表3不同方案下道路网络的速度均方误差