[发明专利]用于可持续交通网络设计的多属性决策软件

权利要求书

1.本发明采用R语言开发了一种用于可持续交通网络设计的多属性决策软件，该软件的算法框架包括以下步骤：

步骤一：上层政策制定者使用Dirichlet分布生成一个随机的交通网络通行能力提升方案Δc，方案编号m＝0，下层出行者作出一系列行为反应；

步骤二：下层模型是交通产生、交通分布、交通方式划分和交通流分配的顺序模型，通过反馈迭代达到交通系统平衡，可以计算出平衡状态时的路段交通流量和通行时间；

步骤三：计算交通系统平衡时可持续交通的四个属性值：

·经济属性

·环境属性

·社会属性

·安全属性；

步骤四：返回步骤一，使用Dirichlet分成另外一个随机的交通网络通行能力提升方案Δc，方案编号m＝m+1，当方案数目达到预先定义的M(M≥200)时，即m＝M时，停止循环，转到步骤五；

步骤五：使用多属性决策方法找出最优的交通网络设计方案。

2.对于权利1中的步骤二，采用如下的计算过程：

步骤1：从Dirichlet分布Dir(α)得到通行能力提升模式Δc；

步骤2：通过均匀分布初始化交通分布矩阵设置n＝0，表示迭代次数；

步骤3：通过Frank-Wolfe算法基于用户均衡将交通分布矩阵分配给交通网络，以计算每个路段a上的交通流量和出行时间，之后，起点i和目的地j之间的最短出行时间，即可以通过Dijkstra算法求得；

步骤4：基于采用目的地选择模型来更新交通分布矩阵

步骤5：利用权重递减的MSA对交通分布矩阵和求平均

步骤6：使用相对根平方误差(RRSE)检查交通分布矩阵的收敛性

如果满足收敛条件，则转步骤8，否则转步骤7；

步骤7：令n＝n+1，然后通过Frank-Wolfe算法基于用户均衡将交通分布矩阵分配给交通网络，以计算每个路段a上的交通流量和出行时间，之后，起点i和目的地j之间的最短出行时间，即可以通过Dijkstra算法计算，反馈到步骤4；

步骤8：输出交通分布矩阵以及路段a上的交通流量v_a和起点i与目的地j之间的出行时间

3.对于权利1中的步骤五，采用如下的计算过程：

步骤1：计算每个方案下交通系统均衡的四个属性，即经济，环境，社会和安全；

步骤2：由于属性的量纲不同，它们必须在整合之前进行标准，对于数值越大越好的属性，标准化公式是

其中y_ij是模式i的属性j的值；是属性j的最坏情况；是属性j的最佳情况；z_ij是模式i的属性j的标准化值，对于数值越小越好的属性，标准化公式是

步骤3：采用特征向量法确定属性的权重w_j；

步骤4：网络设计方案i的综合得分是

步骤5：在计算出所有的M个得分之后，可以对网络设计方案进行排序，并且可以选择得分最高的网络设计，这就是最佳设计。

4.权利1、2、3中的算法以R语言编写，具体归纳如下：

#可持续交通网络设计的多属性决策问题：双层模型与算法

#步骤1：初始化，按格式输入数据和必要的包；

#1.1加载计算最短路径的包，准备调用diikstra最短路径算法，注意igraph包首次使用需要安装，然后才能调用；

#install.packages(″igraph″)#安装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，900，4.00，2，1，12，9.0，700，4.00，3，4，5，9.0，700，4.00，4，4，9，12.0，900，7.00，5，5，6，3.0，800，2.00，6，5，9，9.0，600，4.00，7，6，7，5.0，900，4.00，8，6，10，13.0，500，8.00，9，7，8，5.0，300，4.00，10，7，11，9.0，400，5.00，11，8，2，9.0，700，5.00，12，9，10，10.0，700，6.00，13，9，13，9.0，600，5，00，14，10，11，6.0，700，4.00，15，11，2，9.0，700，5.00，16，11，3，8.0，700，4.00，17，12，6，7.0，300，4.00，18，12，8，14.0，700，9.00，19，13，3，11.0，700，6.00)，ncol＝6，byrow＝T)

colnames(e)＝c(″Road″，″Road origin″，″Road destination″，″Free Time″，″Roadcapacity″，″Road length″)

#e#用于检查程序的断点

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

tge＝3000#总的现状交通需求

d0＝matrix(c(1，1，2，0.2*tge，2，1，3，0.4*tge，3，4，2，0.3*tge，4，4，3，0.1*tge)，ncol＝4，byrow＝T)#初始分配方案

colnames(d0)＝c(″OD pair″，″Origin″，″Destination″，″Demand″)

#d0 #用于检查程序的断点

#自定义的Frank-Wolfe算法函数，注意输入的需求矩阵d形式如d0，交通网络e的形式如上面的e，相对误差0.001；

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])

st0＝numeric(4)#存放初始的各OD对最短行驶时间

colnames(V)＝″Road″

V

#OD对12的最短路径和流量

sp12＝as.vector(b12)#转化为路段标号(Road)

st0[1]＝sum(e[sp12，4])#各路段时间求和

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)

st0[2]＝sum(e[sp13，4])#各路段时间求和

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)

st0[3]＝sum(e[sp42，4])#各路段时间求和

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)

st0[4]＝sum(e[sp43，4])#各路段时间求和

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＝ad.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.ns(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的取值范围，步长不能太长，uniroot要求两端的函数值符号相反，有的函数不一定满足，采用optimize函数可以确保找到一元函数的最优值；

g＝function(lamda){step(lamda)^2}

lamda＝optimize(g，c(0，1))$minimum

lamda

#步骤5：确定新的迭代起点

VS3＝VS+lamda*(VS2-VS)

VS3

#步骤6：收敛性检验

if((sqrt(sum((VS3-VS)^2))/sum(VS))＜0.001)break

VS＝VS3#如果不满足收敛条件则用新点VS3替代原点VS，如此循环直到收敛

}

#步骤7：输出平衡状态的特征矩阵result和OD行驶时间矩阵u；

#步骤7.1：输出平衡状态各路径的流量、通行时间和速度；

result＝cbind(e[，1]，round(VS，0)，tp(VS)，e[，6]/(tp(VS)/60)，e[，5]，round(VS，0)/e[，5])

colnames(result)＝c(″Road″，″Volume″，″Time″，″Speed″，″Road Capacity″，″Levelof Service″)

#步骤7.2：输出各OD行驶时间矩阵u；

g＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝result[，3])#创建图，13为节点的个数，result为步骤7生成的矩阵

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对的行程时间，初始假设各OD行程时间为0

u＝matrix(c(1，1，2，0，2，1，3，0，3，4，2，0，4，4，3，0)，ncol＝4，byrow＝T)

#OD对12的行程时间

sp12＝as.vector(b12)#转化为路段标号(Road)

u[1，4]＝sum(result[sp12，3])#各路段时间求和

#OD对13的行程时间

sp13＝as.vector(b13)#转化为路段标号(Road)

u[2，4]＝sum(result[sp13，3])#各路段时间求和

#OD对42的行程时间

sp42＝as.vector(b42)#转化为路段标号(Road)

u[3，4]＝sum(result[sp42，3])#各路段时间求和

#OD对43的行程时间

sp43＝as.vector(b43)#转化为路段标号(Road)

u[4，4]＝sum(result[sp43，3])#各路段时间求和

u＝cbind(u，st0)#OD对间无障碍行驶时间st0

#以列表的形式输出result矩阵和OD行驶时间矩阵

list(result，u)

}

#fw(e，d0)#用于检查程序的断点

#步骤8：定义目的地选择的多项式Logit函数mlogit，输入为各OD行驶时间时间矩阵u和各地交通需求sg，输出为新的交通分布矩阵；

mlogit＝function(u，sg)

{

d＝numeric(4)

d[1]＝sg[1]*exp(-0.1*u[1，4])/(exp(-0.1*u[1，4])+exp(1-0.1*u[2，4]))

d[2]＝sg[1]*exp(1-0.1*u[2，4])/(exp(-0.1*u[1，4])+exp(1-0.1*u[2，4]))

d[3]＝sg[2]*exp(-0.1*u[3，4])/(exp(-0.1*u[3，4])+exp(1-0.1*u[4，4]))

d[4]＝sg[2]*exp(1-0.1*u[4，4])/(exp(-0.1*u[3，4])+exp(1-0.1*u[4，4]))

cbind(u[，1：3]，d)

}

#步骤9：定义一个给定交通需求d0和sg及交通网络e下综合的交通分布与交通分配交替迭代平衡函数，对于初始交通分布可以求得用户平衡状态时各OD的行驶时间矩阵，用户根据该矩阵重新选择目的地，对于新的交通分布又可以生成新的行驶时间矩阵，该过程一直循环进行，一直到交通分布矩阵不再变化为止；

cda＝function(e，d0，sg){

k＝3

repeat{

d1＝mlogit(fw(e，d0)[[2]]，sg)

k＝k+1

if(sqrt(sum((d1[4]-d0[，4])^2))/sum(d0[，4])＜0.01)break#满足一定的精度要求就停止

d0[，4]＝d0[，4]+(1/k)*(d1[，4]-d0[，4])#这里采用迭代加权法(Method ofSuccessive Average，MSA)，此处采用循环次数的倒数作为权重，随着循环次数的增加而减少；

#print(d1)#用于检查程序的断点

#print(k)#用于检查程序

if(k＝＝100)break#如果循环次数达到100次但还没有满足精度要求也跳出循环

}

#print(d1)

d2＝fw(e，d1)

#print(d2)

list(d1，d2)

}

#现状交通系统表现，用于政策的Before-After比较

ge1＝c(sum(d0[1；2，4])，sum(d0[3：4，4]))#用于检查程序的断点

before＝cda(e，d0，ge1)[[2]][[2]]#用于检查程序的断点

#步骤10：Dirichlet分配法主程序；

#install.packages(″DIRECT″)#安装生成Dirichlet的包，首次需要安装，再次不需要

library(DIRECT)#加载包

set.seed(100)#设定随机种子，这样每次生成的同样的随机数

hsd＝500 #表示Dirichilet样本的个数，对于每个进行比较，找出最优的一个

candicate＝7 #候选路段的数目

rD＝rDirichlet(nsd，rep(1，candicate))#生成Dirichlet随机分布

inv＝2000#投资预算

len＝inv*rD/(0.3*c(5，8，2，3，5，5，7))#备选路段通行能力的增加值

st＝matrix(numeric(nsd*4)，nc＝4)#用于存放每个Dirichlet方案的4个属性

#下面是用来测量安全性的参数

b1＝358.6

b2＝-407.7

b3＝175.3

#计算程序的开始运行时间！！前面都是在定义函数，并不占用计算时间；

timestart<-Sys.time()

for(i in(1：nsd))

{

#i＝2#用于检查程序的

#下面是对每一个分配模型构建一个新的网络图e

Ren＝cbind(c(3，4，5，9，13，16，19)，len[i，])#一个更新方案

colnames(Ren)＝c(″Road″，″Enhancement″)#和路段编号结合起来

Ren2＝merge(e，Ren，by＝″Road″，all＝TRUE)#合并

Ren2[，7][is.na(Ren2[，7])]＝0#把NA用0替换下来，因为NA参与计算的结果还是NA

Ren2[，5]＝Ren2[，5]+Ren2[，7]#与原来的通行能力相加，变成改进后的通行能力

e1＝Ren2[，1：6]#一个新的网络e1

d3＝cda(e1，d0，ge1)#对固定交通需求(1200，800)和交通分布d0调用前面定义的cda函数

#下面求4个属性的值，并存放在st矩阵中

#1.经济属性：出行总时间(分钟)

st[i，1]＝d3[[2]][[1]][，2]％*％d3[[2]][[1]][，3]

#2.环境属性：CO的排放量

st[i，2]＝(0.2038*d3[[2]][[1]][，3]*exp(0.7962*e[，6]/d3[[2]][[1]][，3]))％*％d3[[2]][[1]][，2]

#3.社会属性：政策前后的变化

st[i，3]＝max(d3[[2]][[2]][，4]/before[，4])

#4.安全属性：事故个数

st[i，4]＝sum((b1*d3[[2]][[1]][，6]^2+b2*d3[[2]][[1]][，6]+b3)*d3[[2]][[1]][，2]*e1[，6]*365*10^(-8))

print(paste(″第″，i，″个方案的4个属性分别是″，round(st[i，]，2)))

}

#步骤11：将各属性标准化，这4个属性都是越小越好

zst＝st

zst[，1]＝(max(st[，1])-st[，1])/(max(st[，1])-min(st[，1]))

zst[，2]＝(max(st[，2])-st[，2])/(max(st[，2])-min(st[，2]))

zst[，3]＝(max(st[，3])-st[，3])/(max(st[，3])-min(st[，3]))

zst[，4]＝(max(st[，4])-st[，4])/(max(st[，4])-min(st[，4]))

#步骤12：计算各属性的权重

A＝c(1，1/3，1/2，1/4，3，1，2，1，2，1/2，1，1/2，4，1，2，1)

A＝matrix(A，nc＝4，byrow＝T)

A

a＝apply(A，1，prod)

ai＝a^(1/4)

w＝ai/sum(ai)

w#权重向量

lamda＝sum((A％*％w)/w)/4

CI＝(lamda-4)/3

CR＝CI/0.89

CR#一致性检验

#步骤13：输出各个方案的最终得分

zst％*％w

si＝zst％*％w

which.max(zst％*％w)

si[which.max(zst％*％w)，]

#步骤14：输出最终需要的结果

Ren＝cbind(c(3，4，5，9，13，16，19)，len[which.max(zst％*％w)，])#最优更新方案

colnames(Ren)＝c(″Road″，″Enhancement″)#和路段编号结合起来

Ren2＝merge(e，Ren，by＝″Road″，all＝TRUE)#合并

Ren2[，7][is.na(Ren2[，7])]＝0#把NA用0替换下来，因为NA参与计算的结果还是NA

Ren2[，5]＝Ren2[，5]+Ren2[，7]#与原来的通行能力相加，变成改进后的通行能力

e1＝Ren2[，1：6]#最优的网络e

d3＝cda(e1，d0，ge1)#对交通需求ge1和交通分布d01调用前面定义的cda函数

#print(st[which.max(zst％*％w)，]，digits＝7)#输出此时网络的4个属性，保留7位小数；

formatC(st[which.max(zst％*％w)，]，format＝′f′，digits＝3)

#在formatC()函数中可以用format＝参数指定C格式类型，如″d″(整数)，″f″′(定点实数)，″e″(科学记数法)，″E″，″g″(选择位数较少的输出格式)，″G″，″fg″(定点实数但用digits指定有效位数)，″s″(字符串)，可以用width指定输出宽度，用digits指定有效位数(格式为e，E，g，G，fg时)或小数点后位数(格式为f)write.csv(d3[[1]]，file＝″交通分布矩阵.csv″)#以csv格式保持到当前工作目录

write.csv(d3[[2]][[1]]，file＝″路段平衡结果.csv″)

write.csv(d3[[2]][[2]]，file＝″OD出行时间.csv″)

write.csv(Ren2，file＝″交通网络设计方案.csv″)

getwd()#查看输出文件的保存地址

###计算程序的运行时间

timeend<-Sys.time()

runningtime<-timeend-timestart

print(runningtime)#输出运行时间。