数学建模竞赛优秀选2004-2019年

““题 基于改进K-means聚类和隐马尔可链的汽车行驶工况 要数据预处理（问题一I行驶轨迹图。对行驶轨迹图进行分析可得：本文测于S信号在车辆过隧道或进入建筑时存xyz三轴的加速度，并通过计算相邻数据的速度间隔和时间间隔之比，得到该车辆的加速度。在完成数据初步处理后，对问题一要求中的不良数据类型进行分情况处理：（1）对于时间不连续的数据，做删除和修正处理，其32秒的数据进行插值补充；（2）对于加、度异常数据，设定加、度的阈值对其进行删除处理；（）对于长时间180秒）到三个数据文件的剩余记录数，分别为17497131141145319约85的原始数据，有助于构建准确的汽车行驶工况。汽车的速度满足v1km/h，并且加速度满足ai0.1ms2，本文采取该标准来判定怠速约束条件，用分别标记出怠速状态的起始点位和终止点位；（2）对时间间隔不1秒的数据进行删除操作，保证了运动学片段的连续性；（3）根据标记点位，对时间连续的运动学片段进行提取。得到三个数据文件的运动学片段数量674个，413个，374个。汽车行驶工况的构建（问题三）首先，通过查阅文献对福州的发展程度、交通状况有了初步了解，为问题三的模型构建奠定了理论基础；其次，为了构建一条时间间隔为12001300行工线本出两种汽车行驶工况模型（1基于改进的kmns类和隐马尔链的汽车行驶工况模型，（2）标01线性划的汽车行驶工况模型；然后，确定汽车运动特征评估体系，综合各方考量，本文提出了16个特征指标，即均速度均行驶、怠速比、片续时间速时间、时间比、大 均 度标准差、加速度标准差、度标准差、行驶距离；最后，对特征指标计算其与实际数据的相对误差进行模型检验，通过实验发现：本文两种模型，其相对误差均小于8.032%0-1线性规划的汽车行驶工况模型相较于改进的k-eans聚类和隐马尔链的汽车行驶工况模型性能较弱，与真实值间的相对误的k-means聚和隐马尔链的汽车行驶工况模型作为本最后，关于构建汽车行驶工况模型，本文提出了几点想法：（1）可以从曲线拟合的（2可以根据汽车行驶工况曲线和道路路径分类（由S。：改 k-means；隐马尔链；转移矩阵；多目标线性规划；汽车工况一、问题重 研究背 问题重 二、模型假 三、模型符号说 四、问题一的模型建立与求 问题分 数据预处 GPS速度异常数 时间不连续数 缺失数据补 删除不良数 涉及的计算...................................................................................................预处理后的统计分 五、问题二的模型建立与求 问题分 数据处 工况状态划分标 时间不连续处 提取运动学片 六、问题三的模型建立与求 问题分 运动特征提 数据处 标记工况状 计算特征指标 确定候选片段的特征矩 汽车行驶工况模型的建立与求 基于改进k-means聚类与隐马尔链的行驶工况模型——模型 模型1的求 基于多目标0-1线性规划的行驶工况模型——模型 模型2的求 模型求解与检 模型求 模型检 七、模型评 模型的优 模型的不 模型的改 八、参考文 九、附 一、问题重研究背平均速度、平均行驶速度、平均加速度、平均度、怠速时间比、加速时间比、时WLTC工况等行驶工况构建。2000NEDC行驶工况对汽车产品能耗/排放的认证。然而经过多年实践发现：欧洲的NEDC行驶工况实际油耗与我国认证结果偏差越来越大。因而我国开始转变采用世界轻型车测试循环问题重1：数据预处理由于建筑覆盖或过隧道等，GPS信号丢失，造成所提供数据中的时间不于7秒，紧急刹车最大度在7.5-8180180秒处理。问题2：运动学片段1问题3：汽车行驶工请根据上述经处理后的数据，构建一条能体现参与汽车行驶特征的汽车行驶工况曲线（1200-1300秒），该曲线的汽车运动特征能代表所数据源（经处理后的数按照所构建的汽车行驶工况及汽车运动特征评估体系，分别计算出汽车行驶工况与该城市所数据源（经处理后的数据）的各指标（运动特征）值，并说明二、模型假隔超过180s即认定为长时间停车、堵车。V/aVx/Ti/av/mad/Ta/Td/V/a/mL2v/1v2/t1/%t3/%t4/%t2/V/amax/amin/a/1/1/2平均度最大度问题分根据问题一要求，由于汽车行驶设备直接记录的原始数据会包含一些不良所涉及到的异常数据处理较复杂，需要考虑路段交通状况等因素。题目已知的条件是某城市轻型汽车实际道路行驶的数据，3个数据文件分别为同一辆车在不同时段内的数据。首先，对数据进行初步分析，发现存在不良数据（数据缺失、数据异常、出现零值）等问题，通过制定具体、有效的汽车工况数据的处理加速度、平均度以及发动机转速等指标参数，对不良数据按汽车工况类型进行数据处精度。数据处理的流程图如下图2所示：2数据预处理流数据预处3汽车行驶轨迹（4汽车行驶轨迹（5汽车行驶轨迹（况以及福州市的信息，经度为119.30，纬度为26.08。以下4.2.1节、4.2.2节、4.2.3节GPS速度异常Excel1，2，3GPS1，2GPS车速120km/h3298GPS120km/h，需要对这部分的超速数汽车在高速公行驶的最高车速限制为120km/h，所以文件3中出现的GPS车速大120km/h处理第（1）种情况的异常，即：建筑覆盖或过隧道导致GPS信号丢失，数据存在时间不连续，对于这类不良数据，其相邻两条数据的时间间隔为t1GPS信3s2s的数据，的数据异常值，筛选出300条为0的异常记录条，如下表1截取了文件3的13时经纬2017/11/02002017/11/04002017/11/04002017/11/04002017/11/04002017/11/04002017/11/0400剔除该300条数据。在进行模型优化的时候再通过取平均值加入具体研究中。缺失数据处理第（1）种情况的异常，即：建筑覆盖或过隧道导致GPS信号丢失，数据存(1)单个数据由于汽车经过建筑或隧道的情况下，出现数据通道的不顺畅而导致数据缺失的情的补充，具体计算如式（4-2-1）所示：v(i)v(i1)v(i L(i2L(iL(i) L(i)L2(i1)L2(i vi1vi1i个时刻和第i1GPSL1i1L2i1)L1(i1),L2(i1)分别表示为缺失的时刻位置点前后两点的 数据，L1(i),L2(i)为取均值后的当前缺失时刻的数据。由于汽车行驶设备时间间隔为1秒，通过Excel筛选，得到相邻时间间隔为2（缺失单个数据）的记录数为111条，下表2为（以文件1为例）截取的缺失单一数据的部分数据。2缺失单一数据时GPS速经纬2017/12/182017/12/182017/12/182017/12/182017/12/18(2)多组数据删除不良4.2.2中实现了对第一类不良数据的预处理。该节主要从删除不良数据的角 度在7.5~81是加工况下的不良数据处理，2、3、4是怠速工况下的不良数据处理加、度异车最大度在7.5-8ms2。由于不同车型在紧急刹车的最大度存在差异，可以取最大度为7.75m/s2（度为7.5和8的均值）。这里假设加、度是根据速度差与对x,y,z轴加速度的数据列采取删除操作后，进行加、度的不良数据处理。

titi 7100km/ 0a10m 对加速度设定阈值范围，则加速度的最大值过10ms2 aivivi ti表示从第i1个时刻到第iai表示第i表3加、度删除记录数据文件编加速度异常记录度异常记录123在约束条件下，得到不满足该约束的加度对应的记录行，对加度异常的数据进行删除处理，下表4为加、度异常删除的部分数据，以文件1为例。4度异时GPS车时间度2017/12/201-2017/12/201-2017/12/201-2017/12/201-2017/12/211-2017/12/211-2017/12/2101-2017/12/221-2017/12/2331-2017/12/231-2017/12/241-180s即认定为长时间停车、堵车。因此不良数据类型（3）、（4）、（5）统一归为怠速类型（3）数据处熄火但设备仍在运行。满足ti180s的情况为长期停车。对于1、停车不熄火的行驶指标状态进行分析，汽车停车不熄火状态下，速度v0，不熄火则汽车发动机仍在运行，加速度ai0，以此为约束条件，进行删除处理。后汽车发动机不运行，汽车不会发生位置移动，有：ai0，si0。si表示汽车在i积计算得到的。10km/h的数据记录进行删除操作，具体约束条件maxvi10km/ t 怠速状态的判断方法：据题目所给的关于怠速的可知，怠速状态下，汽车是停止运动的，即满足vi00；同时，根据一般下的怠速时间超过180为异常情况，可以得到下式（4-2-6）的约束条件：0.1m|ai

t 4.2.4涉及的计算得出最终处理后的数据记录数。 （1）换算距离DarccossinL2isinL2i1cosL2(icosL2i1cosL1iL2iL1i表示iL2i表示ir

PiMi Pi表示第iMi表示第ini表示第i个时刻的汽车发动机转速。空燃比iM1iM2

其中，M1i表示第i个时刻的汽车 燃烧所需的空气量，M2i表示第i个时刻的汽车行驶燃烧的汽油量，i表示第i个时刻的汽车空燃比。 iMiM

其中，i表示第i个时刻的汽车扭矩百分比，Mmax表示第i个时刻的汽车在该转速下4.2.5数据记录数见下表5所示：5经处理后的记录数据文件编原删除数据量数据的剩余录123问题分于汽车在行驶过程中，受到路面交通情况的影响，会出现多种怠速、加速、匀速和状态。故定义了运动学片段为汽车从怠速状态开始至下一个怠速状态开始之间的车速区间。数据处工况状态划分标怠速工况：汽车的速度满足v1km/hai0.1ms2加速工况：汽车的速度满足v0kmh，并且加速度满足a0.1ms2工况：汽车的速度满足v0km/h，并且加速度满足a0.1ms2的工作状态；匀速工况：汽车的速度满足v1km/h，并且加速度满足a0.1ms2间隔是否为1秒，作为判断数据是否连续的指标。利用进行时间间隔的处理，仅保留时间间隔t1s的数据集，对于时间间隔t1s的数据，进行相应的删除。首先根5.1.1工况状态划分标准，确定了怠速情况下的汽车速度满足v1kmh，并且加速度满足ai0.1ms2。其次，对一个完整的运动学片段过程进行分析，可知运动学片段的第一个初始时刻和运动学片段的结束时刻，其速度都满足v0kmh，且a0。以此为约束条件，用分别标记出运动学片段的起始点位和终止点位。61的某个运动72的某个运动83的某个运动数据文件的运动学片段数量，分别为674个，413个，374个。问题分1200-1300s1200s作为要绘制的行运动特征6用于描述运动学片16个特征指19最大度平均度234n5678s1或数据处标记工况，工 记为3，匀速状态标记为4，统计得到第i个运动学片段的怠速状态个数为li，加速状态个数为ki，匀速状态个数为hi，状态个数为gi。现截取文件1预处理后的任意一个运动学片段，下表7是该部分运动学片段的标记情况和各个文件工况状态的个数。7工况状态标记和个数统计（部分具体时va编2017/11/10014-32017/11/1-3-3232

9汽车工况状态频下面是这16个特征指标的计算：片段持续 n表示第i个运动学片段的持续时间；ti表示第i个运动学片段的结束时刻与初平均速

v1

nnvjj1平均行驶vv2m

怠速时间1tli1其中，t1表示第i个运动学片段的怠速时间比，li为第i加速时间tki2

时间tgi3

其中，t3表示第i个运动学片段的时间比，gi为第i个运动学片段的状态个数匀速时间thi4

其中，t4表示第i个运动学片段的匀速时间比，hi为第i最高速j其中，vmax表示第i个运动学片段中的速度最大值最大加速

j其中，amax表示第i个运动学片段中的加速度最大值最大aminminaj其中，amin表示第i个运动学片段中的度最大值平均加速

a1

其中，a1表示第i个运动学片段中的平均加速度，ki为第i平

a2其中，a2表示第i个运动学片段中的平均度，gi为第i个运动学片段的状态速度标准

其中，1表示第i个运动学片段中的速度标准差（包括怠速状态）加速度标ajaj其中，2表示第i度标ajaj其中，3表示第i个运动学片段中

s1Di或s2v1 其中，s1是通过换算距离的方式；s2是通过平均速度乘以时间差计算的距离。由于换算距离计算出来的s1是直线距离，s2计算出来的是矢量位移距离，是符合汽车行驶轨迹的，所以在计算行驶距离时，应当使用s2的求解。81的运动学片段特征矩片t21234 片n123458692的运动学片段特征矩片1-234-5-6----续上片n1243456…103的运动学片段特征矩片1-2-3-4-5-6----续上片n1-2-3-4-5-6---- 11文1、2、3的总体特征指文1-2-3-文n 4.878810- - 3.315010- - 4.637510- -性，传达的信息存在一定的，这就给分析和解释特征指标上带来麻烦。因此在利用这16个特征指标进行因子分析，基于改进k-means聚类与隐马尔链的行驶工况模型——模型（一）改进k-means 改进k-means聚类的过程是把不同的运动学片段的特J 分类到不同的聚类群中。改进的k均值聚类是一种非监督而且快速聚类的方法，能够使特针对汽车的行驶工况数据特点确定k个初始聚类中心。2计算第iJij中的第kzj的相对距离，确定出各个运动学片段的归属类别。相对距离的计算如下式（6-4-1）所示：2D

k 2 2

Jk（二）基于隐马夫链的运动学片段合由于汽车在行驶过程中，车速是随城市的交通情况而随机的变化。而马尔链过程工况状态 ，记为YY1,Y2,,Y|1,2,,，把运动 链的状态，记作KK1,K2,,1,2,,k 到状态满足条件概率PYk|Y1p 'k| 条件概率PYk|Y1p定义为状态1到状态的转移概率，在马 链过中，转移概率的物理意义是当前时间为1时状态p的模型，转移到下个时间为时状态k的模型的概率。转移概率矩阵P'为 p1k pP1

2k ⁝ pk

pkk

N

1 聚类类 当前状

12工况状态的转移概加速工工加速工工匀速工怠速工加速工00工0聚类 匀速工0怠速工0加速工00工00聚类 匀速工00怠速工1000基于多目0-1线性规划的行驶工况模型——模型1是从概率的角度构建的行驶工况模型。现在从数据2——特征指标值。基于这法，建立了基于多目标0-1线性规划的行驶工况模型，以下是具0-1线性规划是有两个及两个以上的目标函数，同时满足目标函数和约束条件01，而决策变量取值10-1线性规划的行驶工况征指标尽可能的靠近，16个特征指标的误差都尽可能趋近0。首先，定义汽车实际行驶的特征指标集为 W1,W2,W3,,W15,W16；每一个运动

Ji,1,

,Ji,3,,

, i

, W Wmin W15i,15

W W W W

W

i,ji, 段的匀速时间之和，小于汽车实际行驶中的匀速时间比，运动学片段的时间之和，小于汽车实际行驶中的时间比，决策变量i,j01。2对于模型2多目标0-1线性规划模型的求解，可以采用模糊数学模型，考虑将多对每一个目标Wii,1Ji,1i1,2,ni给出了每个目标的模糊伸缩指标ii0缩指标i的选取要根据各个子目标的重要性进行选择，具体的取值原则为：越是重要的目构建模糊对子目标Wii,1Ji,1构造一个模糊目标Ui，定义其隶属函数为 U

h J i,

ii, i, 0

Ji,ji,jWii,jJi,ji 1 i

1Wii,jJi,jJi,ji,j,Wii,jJi,jiJi,ji,jWii,jJi, i

Wii,jJi,jJi,ji,

，则称UR| W 0or1i, i, i, i,模糊最优 为：MMU，则称满足MmaxMUmaxU i, i, i, i, i,maxWii,jJi,jni UW , 1,2,,

i, i, i, i,

ni

i,

Ji, Wj则上式（6-4-4）被转化为了一个普通的单目标0-1线性规划问题，通过 函数bintprog求解0-1线性规划问题。模型求16741313聚类与聚点距个案聚距个案聚距1 2223114215223111211 12121213141,2,3特征指123-n-15距离聚点最近个案聚距1114222333这一指标，分配1200秒，聚类1、聚类2、聚类3的样本时间长度比例分别为12.9310%，时间区间长度1200秒，下图10是得到的汽车行驶工况曲线。图 汽车行驶工况曲10（a），随机截取一段相同长度的原始信号，两个信号（橙色为）对比可以模型检（1）17特征的相对误特征指实际行驶数代表工相对误差----n4.878810-------（3）0-118为两种模型的相对误差效果。n--2的161的差距较大。故利用的kmns利用多目标线性规划求解汽车行驶工况模型使用较少，是本文一大创新点。后期可通过先进行01线性规划模型求解可能效果会有所提高。七、模型评模型的优多目标0-1线性规划求解汽车行驶工况的模型，以及改进的k-means聚类和隐马尔链模型0-1线性规划模型，目前使用该模型进行汽车行驶工况研究的实践甚少，是本文一大创新点。16个特征指标，实现了从隐信息（工况分析）120km/h的记录数进行删除，并且将时间间隔大于1秒的记录也进行全部删除，保留了85%的原始数据，对后续汽车行模型的不学片段划分。模型的改可以从分析、SOM神经网络模型、三次样条插值，高阶拟线性方程组（从a以及车辆工况状态等数可以进一步地扩展，根据汽车行驶工况曲线和道路路径划分（GPS根据经纬八、参考文LinJ,NiemeierDA.ExploratoryysiscomparingastochasticdrivingcycletoCalifornia'sregulatorycycle[J].AtmosphericEnvironment,2002,36(38):5759-5770.Karande,S.,Olson,M.,andSaha,B.DevelopmentofRepresentativeVehicleDriveCyclesforHybridApplications[J].SAETechnicalPaper2014-01-1900,2014,:10.4271/2014-01-KnezM,MuneerT,JerebB,etal.TheestimationofadrivingcycleforCeljeandacomparisontootherEuropeancities[J].SustainableCitiesandSociety,2014,11:56-60.Ho,Sze-Hwee,Wong,Yiik-Diew,Chang,VictorWei-Chung.DeveloSingaporeDrivingCycleforpassengercarstoestimatefuelconsumptionandvehicularemissions[J].Atmospheric,,.基于组合聚类法的行驶工况构建与精度分析[J].汽车工程,20123(2):164-169.郑殿宇,,陈汉,等.哈尔滨城区乘用车行驶工况的构建[J].公路交通科技,2017,34(4):101-107.[6],,.基于运动学片段的城市道路行驶工况的研究[J].汽车工程,20119(3):256-261.苗强,,白书战,等.基于聚类和马尔可夫链的车典型行驶工况构建[J].中国公路学报,2016,29(11):161-169.+余曼,,,等.基于FCM聚类算法的电动车城市循环工况构建[J].公路交通科技,2018,35(10):140-149.,,陈无畏,等.基于分析的城市道路行驶工况构建的研究[J].汽车工程,201118(1):70-73.杜爱民,,陈礼璠,等.市车行驶工况的和研究[J].同濟大學學(自然科學版2006,34(7943-张璇,刘鹏.一种汽车工况构建中的数据预处理方法[J].汽车实用技术,2017(.城市混合道路行驶工况的构建研究[D].合肥:合肥工业大学,2011.城市道路汽车行驶工况的构建与研究[D].合肥:合肥工业大学,200934(4):,,陈无畏.聚类和马尔方法结合的城市汽车行驶工况构建[J].中国机械工程,201021(23):2893-2897. .合肥市汽车行驶工况的研究[J].汽车技术201035(10):34-,刘海, .基于三次样条插值的车辆行驶数据分析.汽车技术.2013(8):54-7.弘,刘海, .基于三次样条插值的车辆行驶数据分析.汽车技术.2013(8):54-7.,,葛渭高.二阶拟线性微分方程组边值问题的三个对称正解.系统科学与数学.2004;24(4):513-519.,,陈无畏.基于马尔的城市道路行驶工况构建方法[J].农业机械学报,200924(11):26-30.程序1：表格文savename=strcat(num2str(i),'_handle.mat');%data.Properties.VariableNames{1}=data.Properties.VariableNames{2}=data.Properties.VariableNames{3}=data.Properties.VariableNames{4}=data.Properties.VariableNames{6}=data.Properties.VariableNames{7}=data.Properties.VariableNames{8}=data.Properties.VariableNames{9}=data.Properties.VariableNames{11}=data.Properties.VariableNames{12}=data.Properties.VariableNames{13}=data.Properties.VariableNames{14}=data=removevars(data,data=removevars(data,data=removevars(data,%%format %日期型时间变量数据相差较小，需要修改数据类datetime=datenum(data.date(:,1));%data.t(1,1)=1;%假设第一段时间间隔为data.a(1,1)=0;%假设第一段加速度为unit_t=1/60/60/24;%计算单位时间的数值forifdata.t(j,1)=round((datetime(j,1)-datetime(j-%计算上下数据前后时间间data.a(j,1)=(data.v(j,1)-data.v(j-1,1))/data.t(j,1);%计算加速度程序三：异常数据的标43.forj=1:n-flag=ifdel.v(errNumV,1)=errNumV=flag=ifdata.t(j,1)~=del.time(errNumTime,1)=errNumTime=flag=ifdata.a(j,1)>20&flagdel.jiansu(errNumJiasu,1)=errNumJiasu=flag=%%度异ifdata.a(j,1)<a_warning&flagdel.jiansu(errNumJiansu,1)=errNumJiansu=flag=%%长期停车、堵车等长期怠速时间处ifdataflag(j,1)==1&flagifdata.v(j,1)==0(data.mottoV<700|datadel.tingche(errNumTingche,1)=flag=ifsum(data.flag(j:j+179,1)==1)==180&flagdel.taisu(errNumTaisu,1)=errNumTaisu=&data.flag(j+179)==0)&flagfordel.taisu(errNumTaisu,1)=errNumTaisu=errNumTaisu+1;程序四：异常数据的删除ifdata(dell,:)=[];%删除数%data(del.jiasu,:)=[];%删除数formatdatetime=datenum(data.date(:,1));%data.t(1,1)=1;%假设第一段时间间隔为data.a(1,1)=0;%假设第一段加速度为unit_t=1/60/60/24;%计算单位时间的数值fordata.t(j,1)=round((datetime(j,1)-datetime(j-1,1))/unit_t);%计算上下数据前后时间间data.a(j,1)=(data.v(j,1)-data.v(j-1,1))/data.t(j,1);%计算加速度 distance(data.latitude(j,1),data.longitude(j,1),data.latitude(j-1,1),data.longitude(j- 程序五：车辆工况状态标forifabs(data.a(j,1))<0.1|abs(data.a(j,1))ifdata.flag(j,1)=4;%|a|<=0.1data.flag(j,1)=1;%|a|<=0.1&v<1elseifdataflag(j,1)elseifdata.a(j,1)<-dataflag(j,1)=3;%a<-0.1&v>1 17.程序六：Gps<!doctype<htmllang=" <!--原始地址 /ui/1.0/ui/misc/PathSimplifier/examples/simple-demohtml--<base /ui/1.0/ui/misc/PathSimplifier/examples/"<metacharset="utf-<metaname="viewport"content="initial-scale=1.0,user-scalable=no,width=device-<title>简单路径#containerwidth:height: margin: <div<scripttype="text/javascript" /maps?v=1.4.15&key=您申请的key<script <script//创建地图（基于高德地图varmap=newAMap.Map('container',zoom: {if(!PathSimplifier.supportCanva