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

权利要求书

1.一种基于道路负载均衡的城市离散交通网络设计R语言实现方法，步骤如下：

步骤1：建立城市离散交通网络设计模型，采用Wardrop用户均衡原理作为网络用户对不同政策的行为反应，建立了一个双层规划模型用于城市离散交通网络设计，其上层为道路负载均衡，下层为用户平衡模型，上层决策变量为y_a，表示是否修建某条候选路段a，为0-1变量，a∈A，所有的候选路段构成0-1决策向量y，上层决定新建道路方案后，下层形成平衡状态时的路段流量x_a，x_a是决策向量y的函数，表示为x_a(y)，另外，道路网的规划受到资本的约束，假设单位长度的路段修建成本为u_a，则长度为l_a的路段修建成本为u_al_a，因此，双层规划问题表示为：

其中A为候选建设的路段集合；B为新建道路的资金约束；x_a为路段a上的交通流量；为自由流行驶时间，即路段a为空净状态时车辆自由行驶所需要的时间；c_a为路段a的通行能力，即单位时间内路段可通过的车辆数；n为交通网络中路段的数目；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最短路径算法

#install.packages(″igraph″)

library(igraph)

options(digits＝3)

#1.2创建图的距离矩阵，包含所有的候选路段：第一列为路段标号Road，第二列为路段起点标号Road origin，第三列为路段终点标号Road destination，第四列为该路段自由流时间free flow time，第五列为道路通行能力capacity，第六列为道路长度length

#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″，″Road origin″，″Road destination″，″Time″，″Roadcapacity″，″Road length″)

e

#1.3输入交通需求矩阵，第一列为起讫点对的标号OD pair，第二列为起点标号origin，第三列为终点标号destination，第四列为交通需求demand

#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]＝var(con[，2]/e[plan[，i]，5])#假设各路段的流量与通行能力的比值，即路段负载，计算所有路段负载的方差；

}

nv#输出所有方案的方差

new[，which.min(nv)]#输出方差最小的建设方案。