1.本发明采用r语言开发了一种用于可持续交通网络设计的多属性决策软件，该软件的算法框架包括以下步骤：
步骤一：上层政策制定者使用dirichlet分布生成一个随机的交通网络通行能力提升方案δc，方案编号m＝0，下层出行者作出一系列行为反应；
步骤二：下层模型是交通产生、交通分布、交通方式划分和交通流分配的顺序模型，通过反馈迭代达到交通系统平衡，可以计算出平衡状态时的路段交通流量和通行时间；
步骤三：计算交通系统平衡时可持续交通的四个属性值：
·经济属性
·环境属性
·社会属性
·安全属性；
步骤四：返回步骤一，使用dirichlet分成另外一个随机的交通网络通行能力提升方案δc，方案编号m＝m+1，当方案数目达到预先定义的m(m≥200)时，即m＝m时，停止循环，转到步骤五；
步骤五：使用多属性决策方法找出最优的交通网络设计方案。
2.对于权利1中的步骤二，采用如下的计算过程：
步骤1：从dirichlet分布dir(α)得到通行能力提升模式δc；
步骤2：通过均匀分布初始化交通分布矩阵
步骤3：通过frank-wolfe算法基于用户均衡将交通分布矩阵分配给交通网络，以计算每个路段a上的交通流量和出行时间，之后，起点i和目的地j之间的最短出行时间，即
步骤4：基于
步骤5：利用权重递减的msa对交通分布矩阵
步骤6：使用相对根平方误差(rrse)检查交通分布矩阵的收敛性
如果满足收敛条件，则转步骤8，否则转步骤7；
步骤7：令n＝n+1，然后通过frank-wolfe算法基于用户均衡将交通分布矩阵
步骤8：输出交通分布矩阵
3.对于权利1中的步骤五，采用如下的计算过程：
步骤1：计算每个方案下交通系统均衡的四个属性，即经济，环境，社会和安全；
步骤2：由于属性的量纲不同，它们必须在整合之前进行标准，对于数值越大越好的属性，标准化公式是
其中yij是模式i的属性j的值；
步骤3：采用特征向量法确定属性的权重wj；
步骤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)，第二列为路段起点标号(roadorigin)，第三列为路段终点标号(roaddestination)，第四列为该路段自由流时间(freeflowtime)，第五列为道路通行能力(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″，″roadorigin″，″roaddestination″，″freetime″，″roadcapacity″，″roadlength″)
#e#用于检查程序的断点
#1.3输入初始交通需求矩阵d0，第一列为起讫点对的标号(odpair)，第二列为起点标号(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(″odpair″，″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″，″roadcapacity″，″levelofservice″)
#步骤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])#这里采用迭代加权法(methodofsuccessiveaverage，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(iin(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)#输出运行时间。