VISSIM二次开发（Python）&大作业总结2

写在前面

前一小节已经概括了基本所有这次使用的二次开发的知识，在这一小节，我们主要介绍本文对.att文件的读取和绘图以及分析的工作。

在这一部分我们将展开对这一部分的介绍，由于这部分的内容专门性比较强，所以就与前面的二次开发部分分开撰写，供有需要的同学们参考。

.att文件示例

其中表头会介绍数据字节含义，这里需要点出第一列就是仿真运行次数。

需求1：依据仿真结果绘制contour图

#开发环境
matplotlib                         3.3.4
numpy                              1.19.5
pandas                             1.2.4

coutour图是指等高线图，直接上代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as pltS=8  #设置需要读取的仿真次数df1=pd.read_csv('file.att',sep=';',skiprows=21)#att中数据分隔符是“;”，然后需要跳过前面不需要转换成df的部分。
#读取数据需要跳过前面的不需要的行
#%%
# df2=df1.sort_values(by=['DATACOLLECTIONMEASUREMENT', 'TIMEINT'])dfn1=df1[df1['$DATACOLLECTIONMEASUREMENTEVALUATION:SIMRUN'].isin([S])]#这里其实在前面更改clomun后会更加方便。选择仅仅包含这第八次运行的结果
T=int(672*5/8)#这是其中非匝道段的线圈有效数据的行数，因为后面会有很多不需要的统计值，可能需要自己更改
#%%
dfn=dfn1[dfn1['DATACOLLECTIONMEASUREMENT'].isin([1,2,4,7,9])]#这里选取了其中需要分析的截面
pre_df=dfn[['TIMEINT','DATACOLLECTIONMEASUREMENT','SPEED(ALL)']].iloc[:T,:]#这里主要是选择需要分析的几列
pre_df1=pre_df[['SPEED(ALL)']]#最后发现还是只需要一列
pre_df1.index = range(T)#为方便，对index重新编号#下面就是把上面的一列数据转换
for i in range(84):if i ==0:A= pre_df1.iloc[0:5,0]C=pd.DataFrame(A)C.index = range(len(C))else:B=pd.DataFrame(pre_df1.iloc[i*5:i*5+5,0])B.index = range(len(B))C=pd.concat([C,B],axis=1,join='inner')#建立一个9*64的df作为等高线绘图的输入t=C.shape[1]#下面就是等高线图的绘制
C=C.iloc[::-1,:]
x=np.linspace(0,t-1,t)
y=np.linspace(37,41,5)
X,Y=np.meshgrid(x,y)   #上面三句话就是建立等高线图的网格
plt.figure()
plt.contourf(X,Y,np.array(C),7,cmap=plt.cm.RdYlGn,levels=np.linspace(0,90,50)) #绘制等高线图，中间需要以矩阵的形式，并且满足x列，y行
plt.colorbar()#显示色条
plt.xticks(range(0,84,12),['6:00','7:00','8:00','9:00','10:00','11:00','12:00'])#更改横纵坐标
plt.yticks([37,38,39,40,41],['NHNX41','NHNX40','NHNX39','NHNX38','NHNX37'])
# plt.savefig('C:\\Users\\13968\\Desktop\\0819-'+str(S)+'.jpg', dpi=1000)
plt.show()

其中控制绘图效果的主要就是levels=np.linspace(0,90,50)这个代码其中0，90是控制显示的纵坐标色彩范围（对应颜色），50是表述色彩过渡的趋势，可以参考下面链接中的方法进行选用。以及色板的使用都可以参考，当然官方文档写的也很清楚~

contourf的colorbar如何设置显示范围_solmn的博客-CSDN博客_colorbar设置想要的范围

【Python】绘制热力图seaborn.heatmap，cmap设置颜色的参数_小白兔de窝-CSDN博客_cmap

部分同学可能也有绘制真实场景的contour的需求，因为方法一致这里也不进行探讨了，主要就是构建这个绘图速度时空矩阵C，大家可以自己尝试，这部分中间的df传递有些混乱_{目前看来能达到效果就好}

最终绘制的效果如下~是不是还可以

需求2：依据结果进行分析

结果分析主要是参考这两篇文献来做的

Bottleneck Identification and Calibration for Corridor Management Planning

Calibration of a micro-traffic simulation model with respect to the spatial-temporal evolution of expressway on-ramp bottlenecks

评价指标	含义	要求
C1(Bottleneck Area Matching)	实际与仿真在瓶颈范围（时间和空间上）的匹配程度	>0.75
C2(Actual Speed Matching)	反映范围匹配以及实际速度的匹配	>0.70
GEH	由Geoffrey E. Havers提出，旨在比较两组流量数据的匹配程度	<5
速度相对误差		<15%

也简单罗列一下公式，其中C1和Devs还有一些其他的形式不过差别不大，具体怎么算可以参考前面推荐的文献~
Devs=∣Sr−SsjSr\\Devs= \frac { | S _ { r } - S _ { s } j } { S _ { r } } Devs=Sr∣Sr−Ssj

GH=(E−v)2(E+v)/2\sqrt { G H } = \sqrt { \frac { ( E - v ) ^ { 2 } } { ( E + v ) / 2 } } GH=(E+v)/2(E−v)2

C1=∑i=1N{(∑t=1T[BSs(i,t)ABSr(i,t)])⋅(xi+1−xi)}∑i=1N(∑t=1T[BSs(i,t)VBSr(i,t)])⋅(xi+1−xi)}\\C_1=\frac{\sum _ { i = 1 } ^ { N } \{ ( \sum _ { t = 1 } ^ { T } [ B S _ { s } ( i , t ) A B S _ { r } ( i , t ) ] ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}{\sum _ { i = 1 } ^ { N } ( \sum _ { t = 1 } ^ { T } [ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] ) \cdot ( x _ { i } + 1 - x _ { i } ) \}} C1=∑i=1N(∑t=1T[BSs(i,t)VBSr(i,t)])⋅(xi+1−xi)}∑i=1N{(∑t=1T[BSs(i,t)ABSr(i,t)])⋅(xi+1−xi)}

C2=2∑i=1N∑i=1T{[BSs(i,t)VBSr(i,t)]⋅∣Ss(i,t)−Sr(i,t)∣)⋅(xi+1−xi)}∑i=1N∑t=1T{[BSs(i,t)VBSr(i,t)]⋅(Ss(i,t)+Sr(i,t)))⋅(xi+1−xi)}\\C_2=\frac{2 \sum _ { i = 1 } ^ { N } \sum _ { i = 1 } ^ { T } \{[ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] \cdot| S _ { s } ( i , t ) - S _ { r } ( i , t ) | ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}{\sum _ { i = 1 } ^ { N } \sum _ { t = 1 } ^ { T } \{ [ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] \cdot ( S _ { s } ( i , t ) + S _ { r } ( i , t ) ) ) \cdot ( x _ { i } + 1 - x _ { i } ) \}} C2=∑i=1N∑t=1T{[BSs(i,t)VBSr(i,t)]⋅(Ss(i,t)+Sr(i,t)))⋅(xi+1−xi)}2∑i=1N∑i=1T{[BSs(i,t)VBSr(i,t)]⋅∣Ss(i,t)−Sr(i,t)∣)⋅(xi+1−xi)}

这部分代码为好兄弟xjt今天看论文了嘛的博客_CSDN博客提供

#文件导入部分前面有讲解，在此不再赘述，主要介绍主要其功能性模块
#计算BS矩阵（45代表阈值，可自行调节，索引自行依据矩阵调节）BS_wn[f'speed_{data}'] = BS_wn[f'speed_{data}'].apply(lambda x: 0 if x >= 45 else 1)   #这里直接利用pandas读取后的Dataframe类型来进行操作，其他类型也可以操作，只要索引对应上
#修正BS矩阵 怕出错这里可不带入，对结果影响不大
BS_s_1 = BS_wn.values   #转换成矩阵形式，遍历一个线圈的历史数据进行判断
for spa_ii in range(1,BS_s_1.shape[1]):for tt_3 in range(BS_s_1.shape[0]):if BS_s_1[tt_3][spa_ii] == 0:  #主要关注瓶颈区的设置是否合理，因此不合理的瓶颈区设置的异常值就是0，所以找到所有0（非瓶颈区位置）for tt_4 in range(tt_3-3,tt_3):data_listwn = []if tt_4 <0:passelse:for tt_5 in range(tt_4,tt_4+5):if tt_5 == tt_3 or tt_5 >=len(BS_s_1.shape[0])-1: #空间位置索引最大值passelse:data_listwn.append(BS_s_1[tt_5][spa_ii])if data_listwn == [1, 1, 1, 1]:BS_s_1[tt_3][spa_ii] = 1
#指标各参数含义可详见公式
#计算C1
def Ca_C1(Bs_s, Bs_r):C1_up = np.sum(np.logical_and(Bs_s,Bs_r))C1_down = np.sum(Bs_s+Bs_r)C1 = 2*C1_up/C1_downreturn C1
# 计算C2指标：
def Ca_C2(Bs_s1, Bs_r1, S_s, S_r):Bs_or = np.logical_or(Bs_s1,Bs_r1)S_abs = np.abs(S_r-S_s)S_add =S_r+S_sC2_up = np.sum(np.multiply(Bs_or, S_abs))C2_down = np.sum(np.multiply(Bs_or,S_add))C2 =1 - 2*C2_up/C2_downreturn C2
# 计算GEH
def Ca_GEH(V_ss,V_r):Geh_up = np.multiply((V_ss-V_r), (V_ss-V_r))Geh_down = (V_ss+V_r)/2Geh = np.sqrt(Geh_up/Geh_down)return np.percentile(Geh, 0.85) #返回85分位值
#计算速度相对偏差
def Ca_DevS(S_s1, S_r1):S_abs = np.abs(S_r1-S_s1)idxnonzeros = np.where(S_r1 != 0)Devs = S_abs[idxnonzeros]/S_r1[idxnonzeros]Devs_re = Devs[np.where(Devs <= np.percentile(Devs, 0.85))]  #剔除15%的较大的值后进行计算，也可以不剔除Devs_avr = np.average(Devs_re)return Devs_avr

通过计算这些指标再结合这个专栏的第一篇就能进行校正参数的标定~

当然，这些工作都是建立在经过前置的细致的检查之后的，比如软件检查-建模检查-动画检查，经过定性分析之后再进行这里介绍的方法进行定量分析。从而得到最佳的参数组合，完成模型的标定工作，对交通现象和微观行为有更深入的认知。

本系列第一篇：VISSIM二次开发（Python）&大作业总结1_tu_qing的博客-CSDN博客欢迎大家关注