使用遗传算法标定VISSIM 4.3

GITHUB真是个好东西，最近在GITHUB的帮助下，我第一次接触了使用Python调用VISSIM 4.3 COM接口的方法，尽管在其官方文档中没有对此进行详细说明，不过根据测试，VISSIM API中可以通过AttValue()和SetAttValue()函数调用API里边的函数。

环境准备

首先准备环境

• VISSIM 4.30
• Pycharm 2019.03
• Python 3.5

需要的包

• geatpy 2.2.3 用来调用遗传算法
• win32com 用来调用com接口
• numpy 是该遗传算法包的依赖
• xlrd 用来读Excel表的

标定方式概述

目前针对传统微观交通跟驰模型的标定方法主要有两种： 一种是宏观标定，即根据道路中的固定检测器，包括微波、线圈等统计的流量、速度、行程时间等信息，另一种是微观（轨迹）标定，根据航拍（最近几年一般用无人机）得到车辆的轨迹，通过车辆的轨迹数据进行标定。很显然，轨迹标定方法更为精确，不仅可以标定出跟驰的特征，也可以包含换道的信息。但是无人机也有缺陷，一次只能拍一个地点，而且续航时间有限。具体而言，选取哪一种标定方法还是由手头的数据决定的。

数据处理

小时交通量处理

我手头有一个检测器1日的数据，原始数据是每5分钟统计各条车道的车速、以及通过车辆数。首先对该数据进行集计，统计成15分钟车速、车道通过的车辆数数据、大车比例，再通过集合各车道通过车辆数数据，加以扩样得到每15分钟统计一次的 断面小时交通量 这里不再是15分钟流量，是由于VISSIM中设置车辆输入都是以小时交通量为单位。

车道速度处理

随后对车辆车速分车道进行集合统计，得到分车道的速度频率分布累计曲线，这个曲线还是蛮重要的 ，决定了自由流状态下，车辆运行的期望车速分布。查看原始数据，我们可以发现原始数据不仅包含了自由流状态，还包含了拥堵状态，显然拥堵状态是不能作为期望车速输入的。

将拥堵状态和自由流状态分开

首先使用车道流量Q、速度V计算该车道的车辆密度k，利用的是交通流基本关系Q=KV

在Excel中直接选中所有车道速度，绘制K-V散点图，绘制之后可以用直线拟合趋势线，拟合后我们可以发现用直线拟合程度就比较好，各个车道K-V关系拟合R²都在0.85以上。于是我们可以拍一拍脑袋，认为该路段各车道交通流服从格林希尔德模型。拟合出来的直线显示截距，截距就是自由流车速vf（K=0时的速度），vf/2就是拥堵和自由流的临界。

统计速度在各区间内的频数

使用vf/2为临界，筛选出大于vf/2的速度数据，放到另一张表中，进行分段统计，分段长度可以采用10，如果你比较细心也可以用5。

这里有一个小技巧，分段需要用到EXCEL中的FREQUENCY函数，首先选中一列，然后输入函数公式后按ctrl+shift+enter可以填充该列。

最后计算各个区间内的比例，并进行累加，就可以得到速度频率分布曲线，我们在VISSIM中新建4个期望速度分布曲线，把我们得到的4个数据按图索骥画进去。

VISSIM 设置

设置不同车道的期望车速曲线

这里又有需要注意的地方，在VISSIM中，车辆输入的时候是不能输入车道级的期望车速分布，那么问题就来了，我的车辆在一条4车道高速上，怎么样保证他在不同的车道上有不同的期望车速分布呢？我们可以通过设置期望速度决策点来修改不同车道的期望车速。然而，由于路网中的车辆会进行选择性换道，进行选择性变道后，车辆通过上一个期望速度决策点的期望车速会保留到换道后，这就和我们期望的不一样。想要解决也很简单，一个期望速度决策点不够，就来几个。

设置LINK的类型

这个要特别拿出来说一下，天坑！！原来通过一番设置之后代码都跑起来了，可是发现程序每次变换随机参数跑出来的结果都一样（不管调成啥都没用），整整想了一个晚上挠掉N多头发才发现是我LINK的类型选成了城市道路！！！也就是说我调了一晚上高速路的模型参数，压根就没在路网中生效，还我头发啊啊啊啊！！！！

就是这个，要选成3 Freeway 才有用

设置LINK车道限行

这个还是要和实际相符，国内大部分多车道高速都禁止大货车在最内侧车道通行，如果是4车道高速，一般会只让大货车在最外侧两条车道运行。在VISSIM里上面那个界面点Lane Closure 然后把Lane3 和lane4 对HGV关闭（VISSIM最右边车道是1 所以应该关闭lane3 和lane 4）

设置评价

评价也是一个大坑，在4.30版本的VISSIM中，如果你不打开评价，仿真他也不会报错，而是会输出0.00这样令人摸不着头脑的数据，掉无数头发后，才能想到 哦是评价没打开。由于需要检测车辆数和速度，这里要点选上两个数据。同时要注意时间是从900s开始统计，集计间隔也是900s。

这里为了演示只搞了1个检测器，实际上检测器个数肯定不止一个，VISSIM里一个检测器只能看一条车道的数据

遗传算法准备

首先要介绍下遗传算法，有一段解释特别有趣形象：

从前，有一大群袋鼠，它们被莫名其妙的零散地遗弃于喜马拉雅山脉。于是只好在那里艰苦的生活。海拔低的地方弥漫着一种无色无味的毒气，海拔越高毒气越稀薄。可是可怜的袋鼠们对此全然不觉，还是习惯于活蹦乱跳。于是，不断有袋鼠死于海拔较低的地方，而越是在海拔高的袋鼠越是能活得更久，也越有机会生儿育女。就这样经过许多年，这些袋鼠们竟然都不自觉地聚拢到了一个个的山峰上，可是在所有的袋鼠中，只有聚拢到珠穆朗玛峰的袋鼠被带回了美丽的澳洲。

这里的山峰就是优化问题的最大化问题，当然问题也可以是找一个低谷（最小化问题）

我们使用遗传算法，目的是使我们的输入（模型中不同的参数）得到的形成速度、通过车辆数的仿真结果与真实结果差异最小，用数学公式表达，可以写成
$$
min z=sum(|v_s-v_t|/v_t+|Nveh_s-Nveh_t|/Nveh_t)
$$
这里我们是使用了geatpy上的一个模板，仅对编码方式从’BG’ 改为’RI’ （三个参数都是实数范围求解），其他改动不大

直接贴代码：

main.py

# -*- coding: utf-8 -*-
import numpy as np
import geatpy as ea # import geatpy
from MyProblem import MyProblem # 导入自定义问题接口

"""===============================实例化问题对象==========================="""
problem = MyProblem() # 生成问题对象
"""=================================种群设置==============================="""
Encoding = 'RI'       # 编码方式
NIND = 30             # 种群规模
Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 创建区域描述器
population = ea.Population(Encoding, Field, NIND) # 实例化种群对象（此时种群还没被初始化，仅仅是完成种群对象的实例化）
"""===============================算法参数设置============================="""
myAlgorithm = ea.soea_SEGA_templet(problem, population) # 实例化一个算法模板对象
myAlgorithm.MAXGEN = 25 # 最大进化代数
"""==========================调用算法模板进行种群进化======================="""
[population, obj_trace, var_trace] = myAlgorithm.run() # 执行算法模板
population.save() # 把最后一代种群的信息保存到文件中
# 输出结果
best_gen = np.argmin(problem.maxormins * obj_trace[:, 1]) # 记录最优种群个体是在哪一代
best_ObjV = obj_trace[best_gen, 1]
print('最优的目标函数值为：%s'%(best_ObjV))
print('最优的控制变量值为：')
for i in range(var_trace.shape[1]):
    print(var_trace[best_gen, i])
print('有效进化代数：%s'%(obj_trace.shape[0]))
print('最优的一代是第 %s 代'%(best_gen + 1))
print('评价次数：%s'%(myAlgorithm.evalsNum))
print('时间已过 %s 秒'%(myAlgorithm.passTime))

MyProblem.py ：用来描述问题

# -*- coding: utf-8 -*-
import numpy as np
import geatpy as ea
import xlrd
import win32com.client as com  # VISSIM COM
import time

"""

"""


class MyProblem(ea.Problem):  # 继承Problem父类
    flag = False

    def __init__(self):
        name = 'MyProblem'  # 初始化name（函数名称，可以随意设置）
        M = 1  # 初始化M（目标维数）
        maxormins = [1]  # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
        Dim = 3  # 初始化Dim（决策变量维数）
        varTypes = [0] * Dim  # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
        lb = [1.2, 0.9, 2]  # 决策变量下界 CC0推荐：1.2-1.7 cc1推荐：0.9-3.0 CC2推荐：2-7
        ub = [1.7, 3.0, 7]  # 决策变量上界
        lbin = [1, 1, 1]  # 决策变量下边界
        ubin = [1, 1, 1]  # 决策变量上边界
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)

    def aimFunc(self, pop):  # 目标函数
        x = pop.Phen  # 得到决策变量矩阵 x=30
        x0 = x[:, [0]]
        x1 = x[:, [1]]
        x2 = x[:, [2]]
        # x1 = x[:, 1]
        # x2 = x[:, 2]
        # print(x)
        if self.flag == False:
            self.Data = xlrd.open_workbook('resultToCali2.xlsx')
            self.Table = self.Data.sheet_by_name(u'Sheet1')
            self.flag = True
            self.V1 = self.Table.col_values(0)  # spd for segment 1
            self.V2 = self.Table.col_values(1)  # spd for segment 2
            self.V3 = self.Table.col_values(2)  # spd for segment 3
            self.V4 = self.Table.col_values(3)  # spd for segment 4
            self.V5 = self.Table.col_values(6)  # spd for segment 5
            self.V6 = self.Table.col_values(7)  # spd for segment 6
            self.V7 = self.Table.col_values(8)  # spd for segment 7
            self.V8 = self.Table.col_values(9)  # spd for segment 8
            self.relativeFlow = self.Table.col_values(4)  # relative flow
            self.volume = self.Table.col_values(5)  # 流量
        totalResult = []
        for k in range(30):
            # # =================================================================================
            # # VISSIM Configurations
            # # Load VISSIM Network
            self.Vissim = com.Dispatch("Vissim.Vissim");
            dir = "E:\\PycharmProjects\\GA-calibration\\test.inp"
            self.Vissim.LoadNet(dir)
            # Define Simulation Configurations
            graphics = self.Vissim.Graphics
            graphics.SetAttValue("VISUALIZATION", False)  ## 设为 不可见 提高效率
            self.Sim = self.Vissim.Simulation
            self.Net = self.Vissim.Net

            # G = self.Vissim.Graphics
            dbpss = self.Net.DrivingBehaviorParSets  # Driving behavior module
            dbps = dbpss(3)
            # # Set Simulation Parameters
            TotalPeriod = 55802  # Define total simulation period
            WarmPeriod = 900  # Define warm period 10 minutes
            Random_Seed = k  # Define random seed
            step_time = 1  # Define Step Time
            self.Sim.Period = TotalPeriod
            self.Sim.RandomSeed = 42
            # self.Sim.RunIndex= 1
            self.Sim.Resolution = step_time

            # Each scenario run 5 times
            # =================================================================================
            # Data Collection Variables
            t1 = []
            t2 = []
            t3 = []
            t4 = []
            t5 = []
            t6 = []
            t7 = []
            t8 = []
            nVeh = []
            dbps.SetAttValue('CC0', x0[k][0])  ## 天坑！ 不能写x0[k]！！！！！！
            dbps.SetAttValue('CC1', x1[k][0])
            dbps.SetAttValue('CC2', x2[k][0])
            print("第", k, "组参数： CC0=", x0[k][0], ",CC1=", x1[k][0], ",CC2=", x2[k][0])
            eval = self.Vissim.Evaluation
            eval.SetAttValue("TRAVELTIME", True)
            eval.SetAttValue("DATACOLLECTION", True) #打开检测器评价，不然获得的结果是0.00
            TT1 = self.Net.TravelTimes(1)
            TT2 = self.Net.TravelTimes(2)
            dataCollections = self.Vissim.Net.DataCollections
            dt1 = dataCollections(1)   #设置在K0+530断面
            dt2 = dataCollections(2)
            dt3 = dataCollections(3)
            dt4 = dataCollections(4)
            dt5 = dataCollections(5)   #设置在K4+600断面
            dt6 = dataCollections(6)
            dt7 = dataCollections(7)
            dt8 = dataCollections(8)
            composition = self.Net.TrafficCompositions(1)
            vehicleInput = self.Net.VehicleInputs(1)
            # self.Sim.RunContinuous()
            for j in range(1, TotalPeriod):
                if (j % 900 == 0) and (j >= 1801):
                    composition.SetAttValue1("RELATIVEFLOW", 100, 1 - self.relativeFlow[int(j / 900 - 3)])
                    composition.SetAttValue1("RELATIVEFLOW", 200, self.relativeFlow[int(j / 900 - 3)])
                    vehicleInput.SetAttValue("VOLUME", self.volume[int(j / 900 - 3)])
                    t1.append(dt1.GetResult("speed", "mean", 0))  # 车道1
                    t2.append(dt2.GetResult("speed", "mean", 0))  # 车道2
                    t3.append(dt3.GetResult("speed", "mean", 0))  # 车道3
                    t4.append(dt4.GetResult("speed", "mean", 0))  # 车道4
                    t5.append(dt5.GetResult("speed", "mean", 0))  # 车道1
                    t6.append(dt6.GetResult("speed", "mean", 0))  # 车道2
                    t7.append(dt7.GetResult("speed", "mean", 0))  # 车道3
                    t8.append(dt8.GetResult("speed", "mean", 0))  # 车道4  检测器5-8是后来新增的
                    nVeh.append(4 * (dt1.GetResult("NVEHICLES", "sum", 0)
                                     + dt2.GetResult("NVEHICLES", "sum", 0)
                                     + dt3.GetResult("NVEHICLES", "sum", 0)  
                                     + dt4.GetResult("NVEHICLES", "sum", 0))) ## 车辆数的检测
                self.Sim.RunSingleStep()
            spdTotal = np.array(0.125 * sum(abs(np.array(t1) - np.array(self.V1)) / np.array(self.V1)
                                            + abs(np.array(t2) - np.array(self.V2)) / np.array(self.V2)
                                            + abs(np.array(t3) - np.array(self.V3)) / np.array(self.V3)
                                            + abs(np.array(t4) - np.array(self.V4)) / np.array(self.V4)
                                            + abs(np.array(t5) - np.array(self.V5)) / np.array(self.V5)
                                            + abs(np.array(t6) - np.array(self.V6)) / np.array(self.V6)
                                            + abs(np.array(t7) - np.array(self.V7)) / np.array(self.V7)
                                            + abs(np.array(t8) - np.array(self.V8)) / np.array(self.V8)))  
            # NVehTotal = np.array(sum(abs(np.array(nVeh) - np.array(self.volume)) / np.array(self.volume)))  ##
            totalResult.append(spdTotal)  # +NVehTotal)
            # print("spd误差：", spdTotal)
            # print("NVeh误差：", NVehTotal)
            print("误差为：", totalResult[k])
            self.Sim.Stop()
        pop.ObjV = np.vstack(totalResult)  # 计算目标函数值，赋值给pop种群对象的ObjV属性

对于VISSIM 来说，一般Weidemann99跟驰模型，可调整的参数一共有9个CC0-CC9，不过一般调CC0 CC1 CC2三个参数即可，通过定性分析（pai nao dai），我们认为其他参数不太敏感。

对于6.0+版本，还可以调整换道参数，按理说换道参数有一些也是需要调整的，不过在4.3中没有开放换道参数的接口，这里简化起见，我们采用默认参数，仅对“同时观测的车辆数”进行调整，从2调整为4

上述代码常见错误会有： 输入的数组应该是n行1列的numpy array 矩阵，看着就令人头大，处理方法是不能用原来python的数组，应该用np.array系列的数组，经过一段时间的试验后发现还是numpy里的vstack函数靠谱

最终结果：

种群信息导出完毕。
最优的目标函数值为：17.6446518217351
最优的控制变量值为：
1.6210104341403664
1.0379648617231827
2.258785180700528
有效进化代数：25
最优的一代是第 17 代
评价次数：1250
时间已过 71968.43456172943 秒

Process finished with exit code 0

CC0=1.6210104341403664 CC1=1.0379648617231827 CC2=2.258785180700528

TO DO:

• 后续可能会出个6.0+版本，4.3和6的调用还是有很多不一样的地方，很多函数名都变了
• 最优化的函数一般用的是均方根相对误差
• 应该随机5次比较合理，但是耗时太长，可以通过缩小种群规模提高效率
• 总体来说效率还是低了些 不过并非不能接受吧

参考资料

• Geatpy官方文档 （https://github.com/geatpy-dev/geatpy）
• VISSIM 4.30 API 文档
• 杨神写的代码，使用遗传算法标定VISSIM9，学（膜）到了思想 （https://github.com/TransUMD/VISSIM-Calibration-via-LHD-GA）