[發明專利]用于可持續交通網絡設計的多屬性決策軟件

在審

• 申請號: 201910627204.0
• 申請日: 2019-06-27
• 公開（公告）號: CN110457012A
• 公開（公告）日: 2019-11-15
• 發明（設計）人: 林宏志;趙宇軒
• 申請（專利權）人: 東南大學
• 主分類號: G06F8/20
• 分類號: G06F8/20
• 代理公司: 暫無信息
• 代理人: 暫無信息

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)#輸出運行時間。