DIAL算法计算单位流率-python实现

dail算法

DIAL算法又称STOCH算法,它基于路段分配流量,并将OD量仅仅加载到连接OD对的有效路径上,成功避免了路径枚举,被公认为是相当高效的Logit型随机网络加载算法。也就是说,DIAL算法并非利用Logit随机加载模型计算路径流量和路段流量,但结果却与这些公式吻合。

dail算法步骤

对OD起点r,DIAL算法执行如下3个步骤:
步骤1:预处理
(1)对 ∀i∈N\forall i∈ N∀i∈N,计算从r到i的最短路径长度rir_iri;
(2)对 ∀(i,j)∈A\forall(i,j)∈ A∀(i,j)∈A,按下式计算路段似然值LijL_{ij}Lij
Lij={eθ(rj−ri−tij)，ri≤rj0，其他L_{ij}=\left\{ \begin{aligned} e^{\theta(r_j-r_i- t_{ij})}， & {r_i \leq r_j} \\ 0， & 其他 \end{aligned} \right. Lij={eθ(rj−ri−tij)，0，ri≤rj其他
步骤2:正向计算路段权重
从r开始,按rir_iri的升序依次处理每个节点i,即对 ∀j∈F(i)\forall j∈ F(i)∀j∈F(i),按下式计算路段(i,j)的权重wijw_{ij}wij
wij={Lij，i=rLij∑m∈B(i)wmi，其他w_{ij}=\left\{ \begin{aligned} L_{ij}， & i= r \\ \\ L_{ij}\sum_{m∈ B(i)}{w_{mi}} ， & 其他 \end{aligned} \right. wij=⎩⎪⎪⎪⎨⎪⎪⎪⎧Lij，Lijm∈B(i)∑wmi，i=r其他

步骤3:反向分配路段流量
从具有最大rjr_jrj值的节点开始,按rjr_jrj的降序依次处理每个节点j,即对 ∀i∈B(j)\forall i∈ B(j)∀i∈B(j),按下计算路段(i,j)的流量xijx_{ij}xij
xij={0，j=r(qrj+∑∀m∈F(j)xjm)wij∑∀m∈B(j)wmj，其他x_{ij}=\left\{ \begin{aligned} 0， & j= r \\ \\ (q_{rj}+ \sum_{\forall m∈ F(j)}x_{jm})\frac{w_{ij}}{\sum _{\forall m∈ B(j)}w_{mj}} ， & 其他 \end{aligned} \right. xij=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧0，(qrj+∀m∈F(j)∑xjm)∑∀m∈B(j)wmjwij，j=r其他

实例

求一下网络的单位流率（实心顶点为起讫点）

图形数据说明：顶点编号如图，从左到右从上到下，从1为起始号进行编号，在程序中顶点i编号则为i-1

准备

使用工具及环境：
编程语言：python3.6（Anaconda3发行版）
依赖库：numpy，pandas
编译环境：Ubantu -> Vim / Windows -> Pycharm /Spyder
程序运行要求：当使用相对位置读取文件时，解释器的工作文件夹应与数据文件位置保持一致，否则会出现无法读取文件的问题。（解决办法：import os）

流程图

代码

daildemo.py

#!/usr/bin/env python
import numpy as np
import pandas as pd
#import networkx as nxdef get_T_matrix(file):f = open(file)first_ele = Truefor data in f.readlines():data = data.strip('\n')nums = data.split("\t")if first_ele:nums = [int(x) for x in nums]T = np.array(nums,dtype = float)first_ele = Falseelse:nums = [int(x) for x in nums]T = np.c_[T,nums]f.close()nnode = T.shape[0]T0 = T.flatten()for index in range(T0.shape[0]):if T0[index] == 0:T0[index] = np.infT = T0.reshape((nnode,nnode))for i in range(nnode):T[i,i] = 0return Tdef get_information(T):print('请输入起点（1~',T.shape[0],")：",end = '')r = input()r = int(r)-1print("请输入讫点(1~",T.shape[0],")：",end = '')s = input()s = int(s)-1q = input("请输入总需求量（q>=0)：")q = int(q)thita = input("请输入参数（默认为1）：")thita = int(thita)info = (r,s,q,thita)return(info)def dial(T,r,s,q,thita):'''利用dial算法计算单位流率T 路网邻接矩阵r 标号起点s 标号讫点q 总需求量thita dial算法中的参数'''nnode = T.shape[0]print("step1  求最短距离","="*70,sep = '\n')R = all_shortest_path(T,r)# 根据到起点的最短距离进行标号connection = np.zeros((nnode,nnode)) #对于connection[i,j]合理路段为1，其他为0for i in range(nnode):for j in range(nnode):if(T[i,j] != 0 and T[i,j] != np.inf ):if(R[i] < R[j]):connection[i,j] = 1 #i为j的前项print("step2   计算合理路径使用概率","="*70,sep = '\n')P = np.zeros((nnode,nnode))for i in range(nnode):for j in range(nnode):if connection[i,j]:P[i,j] = np.exp(thita*(R[j]-R[i]-T[i,j]))P_f = pd.DataFrame(P)print("路段使用概率：",P_f.round(3),sep = '\n')print("step3 计算线权","="*70,sep = '\n')W = np.zeros((nnode,nnode)) # 初始化线权矩阵for i in range(nnode):for j in range(nnode):if connection[i,j]:if i == r:W[i,j] = 1*P[i,j]else:W[i,j] = P[i,j]*connection[:,i].dot(W[:,i].T)W_f = pd.DataFrame(W)print("线权：",W_f.round(3),sep = '\n')print("step4 计算单位流率","="*70,sep = '\n')Qw = np.zeros((nnode,nnode))for i in range(nnode-1,-1,-1):for j in range(nnode-1,-1,-1):if connection[i,j]:if (connection[:,j].dot(W[:,j].T)):if j==s:Qw[i,j] = q*W[i,j]/(connection[:,j].dot(W[:,j].T))else:Qw[i,j] = Qw[j].dot(connection[j].T)*W[i,j]/(connection[:,j].dot(W[:,j].T))print("单位流率:")Qw_f = pd.DataFrame(Qw)print(Qw_f.round(3))return Qw_fdef all_shortest_path(T,r):nnode = T.shape[0]path = np.zeros((nnode,nnode),dtype = int)Tmin = T.copy()for i in range(nnode):for j in range(nnode):if(T[i,j] != np.inf and T[i,j] != 0):path[i][j] = ifor a in range(nnode):for b in range(nnode):for c in range(nnode):if(Tmin[b,a]+Tmin[a,c]<Tmin[b,c]):Tmin[b,c] = Tmin[b,a]+Tmin[a,c]path[b][c] = path[a][c]print("result:")print("节点标号：",'\n',Tmin[r])return(Tmin[r])# 程序整体流程print("="*70)
print(" "*25,"dail计算单位流率")
print("使用环境：python3.6    依赖库：numpy pandas")
print("班级","姓名","学号",sep = '\t')
print("文件要求：数据文件以txt格式导入，当前工作文件夹与数据文件所在文件夹相同")
print("-"*70)
print('''选择例题对应的数据文件或自定义数据文件:（1.图1  2.图2  3.图3  4.自定义）:''',end = '\t')
flag = input()
file = {'1':'tu1.txt','2':'tu2.txt','3':'tu3.txt'}
if flag == '4':f = input("输入文件名：")
else:f = file[flag]T = get_T_matrix(f) #从文件或内置材料读取路网拓扑数据
info = get_information(T) #获取配流必要的起讫点及流量、参数信息
Qw = dial(T,info[0],info[1],info[2],info[3]) # 向核心算法传参
#数据后处理
flag1 = input("是否保存配流数据文件到工作文件夹（1.保存 2.离开）：")
if flag1 == '1':Qw.to_excel("单位流率.xls")

内置网络文件
2. tu1.txt

0    3   2   4   0   0   0   0   0   0   0
3   0   2   0   3   0   0   0   0   0   0
2   2   0   1   5   7   0   0   0   0   0
4   0   1   0   0   0   3   0   0   0   0
0   3   5   0   0   0   0   3   3   0   0
0   0   7   0   0   0   0   0   3   5   0
0   0   0   3   0   0   0   0   0   2   0
0   0   0   0   3   0   0   0   2   0   4
0   0   0   0   3   3   0   2   0   3   2
0   0   0   0   0   5   2   0   3   0   4
0   0   0   0   0   0   0   4   2   4   0

3.tu2.txt

0    2   2   0   3   0   0   0   0   0
2   0   0   3   0   0   0   0   0   0
2   0   0   4   1   5   0   0   0   0
0   3   4   0   0   2   0   0   0   0
3   0   1   0   0   2   5   0   7   0
0   0   5   2   2   0   0   5   0   0
0   0   0   0   5   0   0   4   6   3
0   0   0   0   0   5   4   0   0   4
0   0   0   0   7   0   6   0   0   2
0   0   0   0   0   0   3   4   2   0

4.tu3.txt

0    1   3   0   0   0   0   0   0   0   0   0   0   0   0
1   0   4   2   0   0   0   0   0   0   0   0   0   0   0
3   4   0   0   2   0   0   0   0   0   0   0   0   0   0
0   2   0   0   3   1   0   0   0   0   0   0   0   0   0
0   0   2   3   0   2   3   0   5   0   0   0   0   0   0
0   0   0   1   2   0   0   4   5   0   0   0   0   0   0
0   0   0   0   3   0   0   0   3   0   0   0   0   9   0
0   0   0   0   0   4   0   0   6   5   7   0   0   0   0
0   0   0   0   5   5   3   6   0   3   0   0   6   0   0
0   0   0   0   0   0   0   5   3   0   0   2   0   0   0
0   0   0   0   0   0   0   7   0   0   0   4   0   0   0
0   0   0   0   0   0   0   0   0   2   4   0   4   0   3
0   0   0   0   0   0   0   0   6   0   0   4   0   5   6
0   0   0   0   0   0   9   0   0   0   0   0   5   0   0
0   0   0   0   0   0   0   0   0   0   0   3   6   0   0

运行结果

参考资料

[1] 可直接计算转向流量的改进型DIAL交通分配算法，任刚,王炜，中国公路学报Vol.18 No.4 Oct. 2005
[2] 配流DAIL算法-刘志祥