本发明针对基于道路负载均衡的城市离散交通网络模型，提供了一种求解算法并使用r语言程序来实现，属于交通工程技术领域。
背景技术：
：
交通网络设计问题是一个考虑使用者路径选择行为，并且在给定约束条件下，选择改建或者新建的路段，从而使得某种网络性能最优化的问题。该研究问题是交通规划领域研究的重点和难点，极具挑战性，引起了许多学者的研究兴趣。交通网络设计问题属于典型的领导者-跟随者的双层规划问题，其上层问题为交通网络的某种性能最优化，下层问题通常为用户平衡问题。常见的网络性能包括系统阻抗、环境污染、备用能力、公平性、可达性、不确定性等。根据决策变量的不同，交通网络设计问题又分为离散交通网络设计、连续交通网络设计问题和混合交通网络设计问题。其中，离散交通网络设计是指在投入资金有限的情况下，采用定量方法研究在现有路网上选择新建某些路段的问题，属于交通规划的方案设计部分。为了均衡利用道路资源，为了给道路使用者提供相似的驾驶环境，为了提高交通网络的可持续性和安全性，本发明建立了以道路负载均衡为目标的城市离散交通网络设计方法，并给出了求解算法。
本发明使用r语言实现模型求解，总体来说具有以下技术特点：1)开源免费：目前国外的交通规划类软件普遍价格高昂，而r语言是开源免费的；2)具有可操作性：r语言简单易用，适于一般的交通工程师和规划师学习采用；3)科学性：本软件的路段交通流量预测采用了广泛接受的wardrop用户均衡原理；4)系统性：本程序考虑整个交通系统的安全性，而非单一路段或事故黑点。
技术实现要素：
：
技术问题：现阶段，交通规划软件基本都是国外产品，价格高昂并且操作复杂不易使用，特别是实际交通网络中常常存在道路负载不均衡的问题，即部分路段交通流量很大，而部分路段交通流量很少，该问题导致交通网络的可持续性和安全性降低。首先，由于长期超负荷运行，拥挤道路的路基和路面都易受到损坏，降低该路段的通行能力，降低交通网络的可持续性。其次，由于各路段的驾驶环境不同，容易诱发交通事故，降低交通网络的安全性。因此，在交通网络规划阶段，应该建立以道路负载均衡为目标的交通网络设计方法，建立城市离散交通网络模型并且设计求解模型的r语言方法。
技术方案：本发明提供了一种基于道路负载均衡的城市离散交通网络设计的r语言实现方法，主要包括以下步骤：
步骤1：建立城市离散交通网络设计模型。本发明采用wardrop用户均衡原理作为网络用户对不同政策的行为反应，建立了一个双层规划模型用于城市离散交通网络设计，其上层为道路负载均衡，下层为用户平衡模型。上层决策变量为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)colnames(e)＝c(″road″，″roadorigin″，″roaddestination″，″time″，″roadcapacity″，″roadlength″)
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]＝var(con[，2]/e[plan[，i]，5])#假设各路段的流量与通行能力的比值，即路段负载，计算所有路段负载的方差；
}
nv#输出所有方案的方差
new[，which.min(nv)]#输出方差最小的建设方案
有益效果：本软件采用开源免费的r语言，并且方便交通工程师使用，特别是提供一种基于道路负载均衡的城市离散交通网络设计方法，目的是为了均衡利用道路资源，为了给道路使用者提供相似的驾驶环境，为了提高交通网络的可持续性和安全性。本发明建立了以道路负载均衡为目标的城市离散交通网络设计方法，给出了求解算法并设计了r语言实现方法。最后，采用交通网络分析中常用的nguyen-dupuis网络验证了该软件的有效性。
附图说明：
图1是迭代优化算法的流程图。
图2是nguyen-dupuis测试网络，其中虚线段为候选建设路段。
图3是不同方案下道路负载均衡的水平。
具体实施方式：
下面结合说明书附图，对本发明作进一步说明：
步骤1：城市离散交通网络设计模型
通常使用道路实际流量与该路段通行能力的比值作为道路的负载水平。假设路段a的流量为xa，通行能力为ca，则该路段的负载水平表示为xa/ca。本发明提出交通网络设计应该以各路段负载水平保持均衡作为政策目标，即以各路段负载水平的方差最小为政策目标，表示为：
其中n为交通网络中路段的数目。考虑到拥挤效应，路段行驶时间为交通流量的函数ta(xa)，其中xa为路段交通流量。
本发明采用wardrop用户均衡原理作为网络用户对不同政策的行为反应，建立了一个双层规划模型用于城市离散交通网络设计，其上层为道路负载均衡，下层为用户平衡模型。上层决策变量为ya，表示是否修建某条候选路段a，为0-1变量，a∈a，所有的候选路段构成0-1决策向量y。上层决定新建道路方案后，下层形成平衡状态网络流xa，也就是说路段流量xa是决策向量y的函数，表示为xa(y)。另外，道路网的规划受到资本的约束。假设单位长度的路段修建成本为ua，则长度为la的路段修建成本为uala。因此，双层规划问题表示为：
其中a为候选建设的路段集合；b为新建道路的资金约束；xa为路段a上的交通流量；
对于行驶时间函数ta(xa，ca)的研究，既有通过实测数据进行回归分析的，也有进行理论研究的。其中被广泛应用的是由美国公路局(bureauofpublicroad，bpr)开发的函数，被称为bpr函数，形式为：
其中α和β为阻滞系数，在美国公路局交通流分配程序中，取值分别为α＝0.15和β＝4，也可以由实际数据回归分析求得，ca为路段a的通行能力。
步骤2：迭代优化算法
城市交通网络设计的双层规划模型是一个np-hard的问题，是一个具有十分挑战的问题。学者们提出了很多方法，主要有迭代优化算法(ioa)、基于(近似)梯度的算法、遗传算法、模拟退火法等。本发明采用迭代优化算法，其上层采用枚举法，下层采用frank-wolfe算法。算法的基本思路是对上层满足约束的可行方案计算下层平衡网络流量和路段速度，再根据路段速度计算上层的目标函数，比较所有可行的方案，最后确定最优的目标函数方案。详细的迭代优化算法可以归纳如下：
step1：生成一个建设方案。判断其是否满足上层约束，如果不满足再生成一个新的建设方案，直到得到一个可行的建设方案y，令m＝1。
step2：初始化：按照
step3：更新各路段的阻抗：
step4：寻找下一步的迭代方向：按照
step5：确定迭代步长：求满足下式的λ：
step6：确定新的迭代起点：
step7：收敛性检验：如果
step8：计算道路负载均衡水平：根据公式(1)计算该方案下道路负载均衡水平dm。转入step1，计算下一个可行方案。当所有可行方案计算完毕后，停止计算，从所有可行方案中寻找上层目标函数最优的建设方案。迭代优化算法的流程图如图1所示
步骤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)colnames(e)＝c(″road″，″roadorigin″，″roaddestination″，″time″，″roadcapacity″，″roadlength″)
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(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]＝var(con[，2]/e[plan[，i]，5])#假设各路段的流量与通行能力的比值，即路段负载，计算所有路段负载的方差；
}
nv#输出所有方案的方差
new[，which.min(nv)]#输出方差最小的建设方案
步骤4：算例分析
本发明采用交通网络分析中常用的nguyen-dupuis网络(图2)作为测试网络。道路参数列在表1中，其中路段编号20-24为候选的路段。为简化计算，本例假设政策制定者要求从5条候选道路中选择2条进行开工建设。交通需求矩阵如表2所示。假设候选路段的结合为a，则该算例的上层目标函数简化为：
表1nguyen-dupuis网络的路段参数
表2nguyen-dupuis网络的od对交通需求
对所有可行方案进行交通流平衡分配后，可以得到该方案下各路段的实际交通流量，进而可得各路段的负载水平及整个交通网络的负载均衡水平。从表3和图3可以看出，不同方案下，网络负载均衡水平相差较大，最大为修建路径20和23，道路负载的方差为0.265，最小为修建路径20与24，道路负载的方差为0.095。因此，为了设计道路负载最为均衡的交通网络，应该修建路径20与24。
表3不同方案下道路网络的负载均衡水平