Traffic Distribution

这里整理一下四阶段法中交通分布预测阶段中的重要预测方法。具体包括：平均系数法、底特律法、Frater法、Funess法、单约束重力模型法、双约束重力模型法，并对比一下各种方法的求解收敛速度。

import numpy as np
import matplotlib.pyplot as plt
import math

首先输入所需数据，包括现状OD表、规划年发生吸引量以及小区间阻抗矩阵

#定义初始OD矩阵
init_OD=np.array([[200, 100, 100],[150, 250, 200],[100, 150, 150]],dtype=float)
init_O=np.array([400, 600, 400],dtype=float)
init_D=np.array([450, 500, 450],dtype=float)
#定义规划年数据
later_O=np.array([1000,1000,1250])
later_D=np.array([1250,900,1100])
#阻抗矩阵
cost_mat=np.array([[14, 32, 40],[32, 16, 22],[40, 22, 16]])
stop = 0.03

再通过几种方法求解远景年OD矩阵（各方法具体步骤及代码封装在最后）

并记录各种方法求解的收敛速度，进行比较。各个算法收敛效率比较如下图。

x_ave,y_ave=ave_para(init_OD, init_O, init_D, later_O, later_D, stop)
x_de,y_de=detroit(init_OD, init_O, init_D, later_O, later_D, stop)
x_fra,y_fra=frator(init_OD, init_O, init_D, later_O, later_D, stop)
x_fu,y_fu=funes(init_OD, init_O, init_D, later_O, later_D, stop)
x_gra_single,y_gra_single=gravity_model_single(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat)
x_gra_double,y_gra_double=gravity_model_double(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat)
plt.figure(figsize=(8, 6))
plt.title('收敛曲线(系数法)')  # 折线图标题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 显示汉字
plt.xlabel('迭代次数')  # x轴标题
plt.ylabel('指标')  # y轴标题
plt.xlim([0, 10])
plt.ylim([0, 2.5])
plt.plot(y_ave, label='平均系数法')  # 绘制折线图，添加数据点，设置点的大小
plt.plot(y_de,label='底特律法')
plt.plot(y_fra,label='福莱特法')
plt.plot(y_fu,label='福尼斯法')
plt.plot(y_gra_single,label='单约束重力')
plt.plot(y_gra_double,label='双约束重力')
plt.axhline(0.03, color='red',linestyle='--')
plt.legend()
plt.show()

平均系数法

该方法思路为：qijq_{ij}qij 的增长与 iii 区产生量的增长及jjj分区吸引量的增长同时相关，而且相关的程度也相同。即 fave(FpiFaj)=12(Fpi+Fai){f_{{\rm{ave}}}}\left( {{F_{pi}}{F_{aj}}} \right) = \frac{1}{2}\left( {{F_{pi}} + {F_{ai}}} \right)fave(FpiFaj)=21(Fpi+Fai)

def ave_para(init_OD, init_O, init_D, later_O, later_D, stop):f_O=later_O/init_Of_D=later_D/init_Drecord=[]a=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))while np.max(a)>stop:for i in range(init_OD.shape[0]):init_OD[i,:]=init_OD[i,:]*(f_O[i]+f_D)/2init_O=init_OD.sum(axis=1)init_D=init_OD.sum(axis=0)f_O=later_O/init_Of_D=later_D/init_Da=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))return init_OD,record

Detroit法

此方法思路为，qijq_{ij}qij的增长与i分区产生量增长率 FpiF_{pi}Fpi成正比，而且还与jjj分区吸引量增长占整个区域吸引量增长的相对比率成正比。全区域在现年和规划年的吸引总量分别为∑jAj/∑jAj0\sum_{j} A_{j} / \sum_{j} A_{j}^{0}∑jAj/∑jAj0。于是Detroit增长函数为 fD(Fpi,Faj)=Fpi⋅FajQ/Q0=PiPi0⋅Aj/Aj0∑jAj/∑jAj0f_{D}\left(F_{p i}, F_{a j}\right)=F_{p i} \cdot \frac{F_{a j}}{Q / Q^{0}}=\frac{P_{i}}{P_{i}^{0}} \cdot \frac{A_{j} / A_{j}^{0}}{\sum_{j} A_{j} / \sum_{j} A_{j}^{0}}fD(Fpi,Faj)=Fpi⋅Q/Q0Faj=Pi0Pi⋅∑jAj/∑jAj0Aj/Aj0

def detroit(init_OD, init_O, init_D, later_O, later_D, stop):f_O=later_O/init_Of_D=later_D/init_Drecord=[]a=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))while np.max(a)>stop:for i in range(init_OD.shape[0]):init_OD[i,:]=init_OD[i,:]*f_O[i]*f_D/(np.sum(later_D)/np.sum(init_D))init_O=init_OD.sum(axis=1)init_D=init_OD.sum(axis=0)f_O=later_O/init_Of_D=later_D/init_Da=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))return init_OD,record

Frator法

1954年Frator提出了分别从产生区和吸引区两个角度分析计算qijq_{ij}qij，然后平均的方法。
先从产生区考虑：
① Frator认为，qijq_{ij}qij与 iii区出行量中j分区的“相对吸引增长率” bijb_{ij}bij成正比。由于规划年从 iii区产生的出行量中被分区j吸引去的出行量为，因此，这个相对吸引增长率为：
bij=qij0Faj∑jqij0Faj{b_{ij}} = \frac{{q_{ij}^0{F_{aj}}}}{{\sum\limits_j {q_{ij}^0{F_{aj}}} }}bij=j∑qij0Fajqij0Faj
② Frator还认为：qijq_{ij}qij也应与i分区规划年的产生量成正比；
③ 综上两点，得qij1′=Pi⋅bij=Pi0⋅Fpi0⋅qij0⋅Faj0∑jqij0⋅Faj0q_{ij}^{1'} = {P_i} \cdot {b_{ij}} = P_{i}^{0} \cdot F_{pi}^{0} \cdot \frac{{q_{ij}^{0} \cdot F_{aj}^{0}}}{{\sum\limits_j {q_{ij}^{0} \cdot F_{aj}^{0}} }}qij1′=Pi⋅bij=Pi0⋅Fpi0⋅j∑qij0⋅Faj0qij0⋅Faj0
再从吸引区j分区的角度分析，同样，类似于以上三步，得qij1"=Aj0⋅Faj0⋅qij0⋅Fpi0∑jqij0⋅Fpi0q_{ij}^{1"} = A_{j}^{0} \cdot F_{aj}^{0} \cdot \frac{{q_{ij}^{0} \cdot F_{pi}^{0}}}{{\sum\limits_j {q_{ij}^{0} \cdot F_{pi}^{0}} }}qij1"=Aj0⋅Faj0⋅j∑qij0⋅Fpi0qij0⋅Fpi0
由于上两式的、是表示同一个量，故预测值应取其平均值：

def frator(init_OD, init_O, init_D, later_O, later_D, stop):f_O=later_O/init_Of_D=later_D/init_Drecord=[]a=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))while np.max(a)>stop:for i in range(init_OD.shape[0]):init_OD[i,:]=init_OD[i,:]*(f_O[i]*f_D)*(init_O[i]/(init_OD[i,:]*f_D).sum()+init_D/np.array([(init_OD[:,0]*f_O).sum(),(init_OD[:,1]*f_O).sum(),(init_OD[:,2]*f_O).sum()]))/2init_O=init_OD.sum(axis=1)init_D=init_OD.sum(axis=0)f_O=later_O/init_Of_D=later_D/init_Da=np.vstack((f_O,f_D))-1a=np.maximum(a, -a)record.append(np.max(a))return init_OD,record

Fueness约束条件法

Fueness1956年提出另一种增长率法。他认为，两个分区之间出行分布量qijq_{ij}qij的预测值与此两个分区之间出行分布的现状值 qij0q_{ij}^0qij0成正比，还与产生分区的规划年产生量预测值、吸引分区的规划年吸引量预测值有关，这种关系可用两个系数uiu_{i}ui、vjv_{j}vj表示（分别称之为产生系数、吸引系数），即qij=qij0⋅uivjq_{ij} = q_{ij}^0 \cdot {u_i}{v_j}qij=qij0⋅uivj
但这两个系数不是简单地等于产生量或吸引量的增长率Pi / P0i、Aj / A0j，而是必须满足两个约束条件：
∑iqij=ui∑jqij0⋅vj=Pi\sum\limits_i {q_{ij}} = {u_i}\sum\limits_j {q_{ij}^0 \cdot {v_j}} = {P_i}i∑qij=uij∑qij0⋅vj=Pi
∑jqij=vj∑iqij0⋅ui=Aj\sum\limits_j {q_{ij}} = {v_j}\sum\limits_i {q_{ij}^0 \cdot {u_i}} = {A_j}j∑qij=vji∑qij0⋅ui=Aj
算法:
Step1：初始化。令所有的ui0=1.0u_{i}^{0}=1.0ui0=1.0，k=0k=0k=0。
Step2：用方程组求解vjv_{j}vj。此时方程组简化为:∑jqij0⋅vj=Pi\sum\limits_j {q_{ij}^0 \cdot {v_j}} = {P_i}j∑qij0⋅vj=Pi,这其实是一个线性方程组，用线性代数的知识不难求解
Step3：用vjkv_{j}^{k}vjk代入方程组（5-16b）求解uik+1u_{i}^{k+1}uik+1，这仍是一个解线性方程组的问题。
Step4：再用新求出的uik+1u_{i}^{k+1}uik+1，代入方程组求解vjk+1v_{j}^{k+1}vjk+1。
Step5：检验收敛性。对所有的iii、jjj，考察ujku_{j}^{k}ujk与ujk+1u_{j}^{k+1}ujk+1、vjkv_{j}^{k}vjk与vjk+1v_{j}^{k+1}vjk+1的相对偏差<3%？若是，ujku_{j}^{k}ujk与vjkv_{j}^{k}vjk为之所求，停止；否则返回第3步。

def funes(init_OD, init_O, init_D, later_O, later_D, stop):record=[]f_O=later_O/init_Ofor i in range(init_OD.shape[0]):init_OD[i,:]=init_OD[i,:]*f_O[i]init_D=init_OD.sum(axis=0)a=later_D/init_D-1a=np.maximum(a, -a)record.append(np.max(a))while np.max(a)>stop:f_D=later_D/init_Dfor j in range(init_OD.shape[1]):init_OD[:,j]=init_OD[:,j]*f_D[j]init_O=init_OD.sum(axis=1)a=later_O/init_O-1a=np.maximum(a, -a)record.append(np.max(a)) if np.max(a)<stop:breakf_O=later_O/init_Ofor i in range(init_OD.shape[0]):init_OD[i,:]=init_OD[i,:]*f_O[i]init_D=init_OD.sum(axis=0)a=later_D/init_D-1a=np.maximum(a, -a)record.append(np.max(a))    return init_OD,record

单约束引力模型

以下以 f(Rij)=e−bRijf\left( {{R_{ij}}} \right) = {e^{ - b{R_{ij}}}}f(Rij)=e−bRij 为例说明单约束引力模型的参数的标定步骤，用被称为“试算法”的算法。意思是首先试探性地给参数bbb取一个初值，用现状PA表和阻抗矩阵进行检验，若不合乎精度要求，分析其原因是因为 bbb值太大还是太小了，据此调整 bbb值，进一步再作检验，直到合乎精度要求为止。

Step1：给 bbb一个初值，如 b=1b=1b=1

Step2: 从模型qij∧=Pi⋅Aj⋅e−Rij∑jAj⋅e−Rij\mathop{{q_{ij}}}\limits^ \wedge = \frac{{{P_i} \cdot {A_j} \cdot {e^{ - {R_{ij}}}}}}{{\sum\limits_j {{A_j}} \cdot {e^{ - {R_{ij}}}}}}qij∧=j∑Aj⋅e−RijPi⋅Aj⋅e−Rij 算得现状的出行量“理论值”qij∧{\mathop q\limits^ \wedge _{ij}}ijq∧（现状PA表中的qijq_{ij}qij被称为“实际值”），得现状理论分布表。

Step3: 计算现状实际PA分布表的平均交通阻抗R‾=1Q⋅∑i∑jqijRij\overline R = \frac{1}{Q} \cdot \sum\limits_i {\sum\limits_j {{q_{ij}}} } {R_{ij}}R=Q1⋅i∑j∑qijRij 再计算理论分布表的平均交通阻抗
R^^=∑i∑j∧qij⋅RijQ\hat{\hat{R}}=\frac{\sum_{i} \sum_{j}{ }^{\wedge} q_{i j} \cdot R_{i j}}{Q}R^^=Q∑i∑j∧qij⋅Rij求两者之间相对误差：

Step4：当δ<0\delta < 0δ<0 ，即R^‾<Rˉ\overline{\hat{R}}<\bar{R}R^<Rˉ，这说明理论分布量小于实际分布量，这是因为参数b太大的缘故，因此应该减少bbb值，令b=b/2b=b/2b=b/2；反之增加b值，令b=2bb=2bb=2b。返回第2步。

def gravity_model_single(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat):b=0.5record=[]delt_R=1while abs(delt_R)>0.03:if delt_R<0:b=b/2else:b=b*2cost_mat_e=np.zeros([3,3],dtype=float)for i in range(init_OD.shape[0]):for j in range(init_OD.shape[1]):cost_mat_e[i,j]=math.exp((-b)*cost_mat[i,j])temp_OD=np.zeros([3,3],dtype=float)for i in range(init_OD.shape[0]):for j in range(init_OD.shape[1]):temp_OD[i,j]=init_O[i]*init_D[j]*cost_mat_e[i,j]/(cost_mat_e[i,:]*init_D).sum()R_temp=(temp_OD*cost_mat).sum()/temp_OD.sum()R_init=(init_OD*cost_mat).sum()/init_OD.sum()delt_R=(R_temp-R_init)/R_initrecord.append(abs(delt_R))#求出了 b，通过远景年P、A求下远景年OD分布吧for i in range(init_OD.shape[0]):for j in range(init_OD.shape[1]):init_OD[i,j]=later_O[i]*later_D[j]*cost_mat_e[i,j]/(cost_mat_e[i,:]*later_D).sum()return init_OD,record

双约束引力模型

同时引进行约束系数和列约束系数的引力模型叫双约束引力模型。双约束引力模型的形式是：
qij=Ki⋅Kj′⋅Pi⋅Aj⋅f(Rij){q_{ij}} = {K_i} \cdot K{'_j} \cdot {P_i} \cdot {A_j} \cdot f({R_{ij}})qij=Ki⋅Kj′⋅Pi⋅Aj⋅f(Rij)
Ki=(∑jKj′⋅Aj⋅f(Rij))−1{K_i} = {\left( {\sum\limits_j {K{'_j} \cdot {A_j} \cdot f({R_{ij}})} } \right)^{ - 1}}Ki=(j∑Kj′⋅Aj⋅f(Rij))−1
Kj′=(∑iKi⋅Pi⋅f(Rij))−1K{'_j} = {\left( {\sum\limits_i {{K_i} \cdot {P_i} \cdot f({R_{ij}})} } \right)^{ - 1}}Kj′=(i∑Ki⋅Pi⋅f(Rij))−1
式中，Ki{K_i}Ki、Kj′K{'_j}Kj′分别为行约束系数、列约束系数。
下面是以 f(Rij)=R−rijf(R_{ij})=R_{-r_ij}f(Rij)=R−rij为例的参数标定算法。

算法步骤：
Step1：给参数γγγ取初值，可参照已建立该模型的类似城市的参数作为估计初值，此处令：γ=1γ=1γ=1。
Step2：用迭代法求约束系数Ki{K_i}Ki、Kj′K{'_j}Kj′：
2-1、首先令各个列约束系数Kj′=1K{'_j}=1Kj′=1；
2-2、将各列约束系数Kj′（j=1,…,n）K{'_j}（j=1, …, n）Kj′（j=1,…,n）代入（2）式求各个行约束系数Ki{K_i}Ki；
2-3、再将求得的各个行约束系数Ki(i=1,…,n){K_i}(i=1, …, n)Ki(i=1,…,n)，代入（3）式求各个列约束系数Kj′K{'_j}Kj′；
2-4、比较前后两批列约束系数，考察：它们的相对误差小于3/%？若是，转至第3步；否则返回2-2步。

Step3：将求得的约束系数Ki{K_i}Ki、Kj′K{'_j}Kj′代入（1）式，用现状Pi{P_i}Pi、AjA{_j}Aj值求现状的理论分布表[q^ij{\hat q_{ij}}q^ij]。
Step4：计算现状实际PA分布表的平均交通阻抗 R‾=1Q⋅∑i∑jqijRij\overline R = \frac{1}{Q} \cdot \sum\limits_i {\sum\limits_j {{q_{ij}}} } {R_{ij}}R=Q1⋅i∑j∑qijRij；再计算理论分布表的平均交通阻抗：R^=∑i∑jqij⋅RijQ\hat{R}=\frac{\sum_{i} \sum_{j} q_{i j} \cdot R_{i j}}{Q}R^=Q∑i∑jqij⋅Rij 。求两者之间相对误差δδδ（见上式）。当∣δ∣<3%∣δ∣<3\%∣δ∣<3%, 接受关于γγγ值得假设，否则执行下一步。
步5：当 δ<0\delta < 0δ<0，即 $ \hat R < \overline R ，这说明理论分布量小于实际分布量，这是因为参数，这说明理论分布量小于实际分布量，这是因为参数，这说明理论分布量小于实际分布量，这是因为参数γ太大的缘故，因此应该减少太大的缘故，因此应该减少太大的缘故，因此应该减少γ值，令值，令值，令γ=γ/2；反之增加；反之增加；反之增加γ值，令值，令值，令γ=2γ$。返回第2步。
算法结束。

def gravity_model_double(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat):b=0.5record=[]delt_R=1while abs(delt_R)>0.03:if delt_R<0:b=b/2else:b=b*2 #初始Ki，Kj,b=1K_i=np.zeros((3),dtype=float)K_j=np.ones((3),dtype=float)temp_K_i=np.zeros((3),dtype=float)temp_K_j=np.zeros((3),dtype=float)delt_K=np.ones([2,3],dtype=float)#阻抗函数cost_mat_f=np.zeros([3,3],dtype=float)temp_K=np.zeros([3,3],dtype=float)for i in range(cost_mat.shape[0]):for j in range(cost_mat.shape[1]):cost_mat_f[i,j]=math.pow(cost_mat[i,j],-b)#求Ki1for i in range(K_i.shape[0]):for j in range(K_j.shape[0]):temp_K[i,j]=K_j[j]*init_D[j]*cost_mat_f[i,j]K_i[i]=1/temp_K[i,:].sum()#求kj1for i in range(K_i.shape[0]):for j in range(K_j.shape[0]):temp_K[i,j]=K_i[j]*init_O[j]*cost_mat_f[i,j]K_j[i]=1/temp_K[i,:].sum()while np.max(delt_K)>0.03:#存放系数Ki，Kjfor i in range(K_i.shape[0]):temp_K_i[i]=K_i[i]temp_K_j[i]=K_j[i]#求Kifor i in range(K_i.shape[0]):for j in range(K_j.shape[0]):temp_K[i,j]=K_j[j]*init_D[j]*cost_mat_f[i,j]K_i[i]=1/temp_K[i,:].sum()#求kjfor i in range(K_i.shape[0]):for j in range(K_j.shape[0]):temp_K[i,j]=K_i[j]*init_O[j]*cost_mat_f[i,j]K_j[i]=1/temp_K[i,:].sum()#求相对误差a=np.vstack(((K_i-temp_K_i)/temp_K_i,(K_j-temp_K_j)/temp_K_j))delt_K=np.maximum(a, -a)#得到Ki、Kjtemp_OD=np.zeros([3,3],dtype=float)for i in range(init_OD.shape[0]):for j in range(init_OD.shape[1]):temp_OD[i,j]=init_O[i]*init_D[j]*cost_mat_f[i,j]*K_i[i]*K_j[j]R_temp=(temp_OD*cost_mat).sum()/temp_OD.sum()R_init=(init_OD*cost_mat).sum()/init_OD.sum()delt_R=(R_temp-R_init)/R_initrecord.append(abs(delt_R))#求出了Ki、Kj、b，通过远景年P、A求下远景年OD分布吧for i in range(init_OD.shape[0]):for j in range(init_OD.shape[1]):init_OD[i,j]=later_O[i]*later_D[j]*cost_mat_f[i,j]*K_i[i]*K_j[j]return temp_OD,record

结语

这里方法也不多，方法对应的算法也不唯一，大家有什么问题，欢迎评论区或私信交流嗷

四阶段法-交通分布预测方法-算法复现相关推荐

四阶段法-出行分布计算底特律法
底特律法计算未来出行分布 % trans_planning 计算未来出行量分布 function main() clc;clear; confirm = 0; while confirm == 0in ...
基于四元素法的捷联惯导姿态更新算法
摘要本文主要介绍了机载捷联惯导系统常用的姿态更新算法--四元素法,并重点介绍了利用四元素法进行姿态更新的一般过程. 关键词:四元素法,连贯导,姿态 1 引言 ...
综合评价模型的缺点_浅谈交通影响评价中不同交通预测方法的特性
目前,据小编经验觉得,我国交通影响评价的工作中仍有很多问题未得到很好的解决,主要表现在:交通需求预测模型有待进一步精细.交通的影响程度评价内容与指标不够明确.报告中提出的交通改善措施很难落实到位等方面 ...
运筹学两阶段法编程c语言,运筹学上机实验 - 单纯形方法的两阶段法
理论部分不解释了, 就是粘个实验课的代码,按照书上的算法写的,仅仅是把课本上的样例过了,有bug可能难免,欢迎指出. Sample 1. $$ \left\{ \begin{aligned} min ...
运筹学上机实验 - 单纯形方法的两阶段法
理论部分不解释了, 就是粘个实验课的代码,按照书上的算法写的,仅仅是把课本上的样例过了,有bug可能难免,欢迎指出. Sample 1. $$ \left\{ \begin{aligned} min ...
【TransCad】采用重力模型法实现交通分布预测
目录前言回顾准备工作小区质心层的导入制作小区连杆建立质心选择集和路网文件重力模型法进行交通分布预测建立阻抗矩阵修改索引题外话重力模型运用写在最后前言明天就要上机练习重力模型 ...
四象限法推导lm曲线_试用几何方法推导IS曲线(四象限法)
[简答题]请根据自己的实际情况,撰写一份适合自己的"征婚启事".内容积极健康.文体,字数等不限 [填空题]由一棵二叉树的前序序列和序列可唯一确定这棵二叉树. (1.0分) [判断 ...
用四级法测量计算机土壤电阻率,四极法测量土壤电阻率的方法详解
1.四极法测量土壤电阻率的步骤 (1)在测量区沿直线等距离插入四根金属电极,电极彼此相距S米,S米为测量深度,插入地表深度为1/20S为宜. (2)将仪表稳定平放地面,调整仪表指针在标准黑线上. (3 ...
基于马尔可夫过程的一种新型混合PSO粒子群算法（SCI二区高被引文献）介绍及算法复现（使用chatgpt）
以下是一篇算法领域的SCI二区文献(原文见附件),介绍了一种使用Markov概率转移矩阵对种群拓扑结构进行加权的粒子群算法,相比于标准PSO算法该算法提高了全局覆盖率,更容易跳出局部最优,但是在局部最 ...
基于Aidlux平台的智慧交通AI安全算法实战
一. 通常来说,智慧交通行业主要是对交通场景内的行人,机动车,非机动车进行识别分析. 行人识别分析包括对行人的姿态,方向,外观,以及基于行人的交通事件识别分析,(行人闯红灯等). 机动车识别分析包括对 ...

四阶段法-交通分布预测方法-算法复现