（四）工况曲线构建
在运动学片段提取阶段，我们成功地提取到了每个文件下的运动学片段。下面，我们将把这些运动学片段放到一起，进行特征值计算，和工况曲线的构建。
理解工况曲线
明确工况曲线的基本含义
在构建工况曲线之前，首先得明白，什么是工况曲线？
一般来说，所谓的汽车行驶工况曲线，就是一段时长为1800秒以内的描述汽车行驶的速度-时间曲线。通过这个工况曲线，可以分析出汽车行驶过程中的一些运动学特征，进而用于评价汽车燃油经济性、污染排放量、以及交通控制方面的风险等等各种重要的研究。
因为汽车行驶过程中的一些参数（如速度、加速度、发动机转速等）都会因为汽车的型号（轻型/重型）、驾驶员驾驶习惯、驾驶路段、驾驶时间和地域等而产生区别。因此，在研究汽车行驶工况曲线构建时，一般都会固定在一个大的区域内，对同一种型号的汽车，不同驾驶员行驶过程做长时间的数据采集，目的在于真实地反映出当地汽车行驶的一般规律。最后再去对数据进行分析，构建工况曲线，那么这条工况曲线就能够反映出当地汽车的行驶的基本规律了。
明确工况曲线的基本内容
汽车行驶工况曲线都是怎么来的？为什么这样表示？
汽车行驶工况曲线，都是由一些运动学片段组合而构建成的。为了能够通过曲线表现出汽车行驶过程中的特性，一般所选取的运动学片段都是具有代表性的短行程片段。例如世界WLTC工况，该行驶工况中，就分别选取了低速、中速、高速、超高速四种行驶工况，一共1800秒运动学片段来构建轻型汽车的行驶工况。其中，各个工况片段出现的位置一般按照工况速度来排序，片段持续时长根据最终的评价结果来调整，使之满意度最大，或者直接按照一定的比例来安排时长。
根据查阅的文献资料，我们了解到，目前的汽车行驶工况构建方法大致为为“短行程法 + 主成分分析 + 聚类分析 + 随机选取”；此外，还有学者使用马尔科夫链的方式来构建工况曲线。
总的来说，最终的目的就是，能够构建出一条工况曲线，时长为1200~1300秒左右，他能够很好的呈现出所采集的所有汽车行驶数据的相应运动特征（如平均速度、平均加速度等），误差越小，则所构建的汽车行驶工况的代表性越好。
确定工况曲线构建模型
基于上面的分析，我们最终确定采取主成分分析+聚类分析+分工况选取和合成+两种方法综合评价的方法，利用提取到的1990个运动学片段进行工况曲线的构建。具体流程如下图所示：
然后，对速度片段数据计算运动学特征指标，那么将会得到 1990*15 的矩阵，其中1990是运动学片段的数量，15是15个运动学特征指标数量。对于上述的15个特征指标来说，直接聚类计算会比较复杂，因此首先进行主成分分析，用更少的指标代替原来的指标，这便是降维的思想。
主成分分析降维后，对其进行K-means聚类分析，即对群体进行聚类，最后聚成3大类（高速、中速、低速三种工况）。这样一来，运动学片段就被分成了三个小组了。
对于每个小组里的运动学片段，分别计算他们的总体分布特征指标（这里我们选取了12个指标），以及这个小组总体分布特征指标的均值。这些指标应该具有一定的规律，怎么样从这些片段中，挑选出最具有代表性的片段来合成工况？这里我们采取的是分别计算小组各个片段的总体分布特征指标与小组总体分布特征指标均值的相似度，并按照相似度降序排列片段，然后以1200-1300s的时长从前往后依次选取片段，来合成工况。于是对于每个小组来说，都能合成出一条1200-1300s的工况曲线。因为每个小组是聚类得到的，且聚类时的指标用的是速度相关的指标，所以各个小组的速度区间就代表了高中低三种工况。合成的曲线也就代表了三种工况下的行驶工况曲线。
为了更好的体现行驶工况，我们继续用三个小组中片段来合成一个组合工况。这里需要解决一个问题，每个小组选多少片段，持续多长时间？我们的做法是，首先按照聚类得到的小组片段，计算出三种工况的速度区间的片段总时长，然后按照这个时长去计算速度区间的时间占比，以此来划分1200-1300s时长。再根据划分的时长去对应的小组挑选片段（片段还是按照相似度最大降序排列）。最后组合得到一条组合工况曲线。
得到四条工况曲线后，需要对工况曲线作出评价。于是我们又重新定义了一些组合特征指标，除了题目中要求的指标外，还有一些附加的指标。对于四条工况曲线，以及原始数据，分别计算组合特征指标，计算构建的工况曲线与原始数据之间的误差。为了更清楚的对比所构建的工况曲线与原始数据之间的误差，我们还使用了速度-加速度联合概率分布图来对比。最后综合两种方式，确定采取哪一种工况曲线最合适。
模型代码
- 合并所有文件的运动学片段
function kinePart = kinePartCombine(filename)
% 载入数据
kinePart ={};
for i = 1:length(filename)
fprintf('正在导入%s运动学片段...\n',filename{i});
load([filename{i},'运动学片段'],'kinepart');
m = length(kinepart);
for j = 1:m
datatemp = cell2mat(kinepart{j,1}(2:end,2));
kinePart{end+1,1} = datatemp;
end
end
fprintf('所有文件运动学片段提取合并完成！\n');
end
此函数负责将三个运动学片段的速度矩阵合并到一起。
- 计算运动学特征指标
function [kineFeature,kineFeature_mat] = calKineFeature(kinePart)
% 定义存储特征值的结构体
kineFeature = struct(...
'run_time',{},... % 运行时间
'neglect_time',{},... % 怠速时间
'accelerate_time',{},... % 加速时间
'slowdown_time',{},... % 减速时间
'unispeed_time',{},... % 匀速时间
'av_speed',{},... % 平均速度
'av_runspeed',{},... % 平均行驶速度
'max_speed',{},... % 最大速度
'speed_std',{},... % 速度标准差
'av_accelerate',{},... % 平均加速度
'max_accelerate',{},... % 最大加速度
'accelerate_std',{},... % 加速度标准差
'av_slowdown',{},... % 平均减速度
'min_slowdown',{},... % 最小减速度
'slowdown_std',{}); % 减速度标准差
m = length(kinePart);
for i = 1:m
% 提取速度
data_temp = kinePart{i};
v_delta = (data_temp(2:end) - data_temp(1:end-1)); % 计算速度差
% 计算运动段特征值
% 时间特征值 单位：s
kineFeature(i).run_time = length( data_temp)-1; % 计算运行时间
kineFeature(i).neglect_time = sum(v_delta == 0); % 计算怠速时间
kineFeature(i).accelerate_time = sum( v_delta > 0.1*3600/1000); % 计算加速时间
kineFeature(i).slowdown_time = sum( v_delta < -0.1*3600/1000); % 计算减速时间
kineFeature(i).unispeed_time = sum( v_delta >= -0.1*3600/1000 &...
v_delta <= 0.1*3600/1000 & v_delta ~= 0); % 计算匀速时间
% 速度特征值 单位：km/h
kineFeature(i).av_speed = mean(data_temp); % 计算平均速度
kineFeature(i).av_runspeed = mean(data_temp(data_temp~=0)); % 计算平均行驶速度
kineFeature(i).max_speed = max(data_temp); % 计算最大速度
kineFeature(i).speed_std = std(data_temp); % 计算速度标准差
% 加速度特征值 单位：m/s^2
v_delta2 = v_delta*1000/3600; %单位换算
kineFeature(i).av_accelerate = mean(v_delta2(v_delta2 > 0.1)); % 计算平均加速度
kineFeature(i).max_accelerate = max(v_delta2(v_delta2 > 0.1)); % 计算最大加速度
kineFeature(i).accelerate_std = std(v_delta2(v_delta2 > 0.1)); % 计算加速度标准差
kineFeature(i).av_slowdown = mean(v_delta2(v_delta2 < -0.1)); % 计算平均减速度
kineFeature(i).min_slowdown = max(v_delta2(v_delta2 < -0.1)); % 计算最小减速度
kineFeature(i).slowdown_std = std(v_delta2(v_delta2 < -0.1)); % 计算减速度标准差
end
kineFeature_mat = cell2mat(reshape(struct2cell(kineFeature),15,length(kineFeature))'); % 运动学特征指标转矩阵存储
fprintf('运动学特征指标计算完成！\n');
end
此函数定义了15个运动学特征指标，主要包括各种工况的时长，以及速度、加速度等。
- 主成分分析降维
function [data2,mylatent] = myPCA(data,proportion)
% 数据标准化
data1 = zscore(data);
% 主成分分析
[coeff,~,latent] = pca(data1);
% 计算累计贡献率，确认维度
latent = latent/sum(latent);
sum_latent = cumsum(latent); % 累计贡献率
dimension = find(sum_latent>proportion);
dimension= dimension(1);
% 降维
data2 = data1 * coeff(:,1:dimension);
mylatent = [latent,sum_latent];% 贡献率 累计贡献率
% 画碎石图
figure;
set(gcf,'outerposition',get(0,'screensize'))
plot(mylatent(:,1),'bo-','linewidth',2)
xlabel('主成分序号'),ylabel('贡献率')
grid on
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24);
print(gcf,'-djpeg','-r300','主成分分析碎石图');
fprintf('主成分分析降维计算完成！\n');
end
此函数负责将1990个个体的15个特征指标降维，并绘制贡献率变化碎石图。
- K-Means聚类分析
function group =myK_means(data,k)
% 调用kmeans算法
% k = 3; % 聚类的类别
iteration =15000 ; % 聚类最大循环次数
distance = 'sqEuclidean'; % 距离函数
opts = statset('MaxIter',iteration);
index = kmeans(data,k,'distance',distance,'Options',opts);
group = cell(k,1);
for i = 1:k
group{i} = find(index==i);
end
save('聚类结果','group');
fprintf('K-Means聚类计算完成！\n');
end
此函数负责将降维后的个体进行聚类，聚成3大类。这里调用的是matlab自带的kmeans聚类函数。为了保证聚类结果稳定，将最大循环次数设置的比较大。
- 按聚类结果计算工况曲线
function [kinePartSorted,curve] = calCurve2(kinePart,group)
for i = 1:3
kinePart_g1 = kinePart(group{i,1},1);
% 计算总体分布特征指标
[~,partDistFeature_mat ]= calDistFeature(kinePart_g1);
% 求各指标均值
partDistFeature_mean = mean(partDistFeature_mat);
% 计算相似度，降序排列
m = length(partDistFeature_mat);
pearsonCorr = zeros(m,1);
for j = 1:m
pearsonCorr(j) = myPearson(partDistFeature_mat(j,:),partDistFeature_mean);
end
[~,mostCorr_index]= sort(pearsonCorr,'descend');
kinePart_g1_sorted = kinePart_g1(mostCorr_index);
% 计算片段时长
part_time = zeros(m,1);
for k = 1:m
part_time(k) = length(kinePart_g1_sorted{k,1});
end
% 选取 1200 ~ 1300s 片段
is_choose = choosePart(part_time,1200,1300);
kinePart_curve = kinePart_g1_sorted(is_choose==1,1);
% 得到最终的工况曲线
curve1 = cell2mat(kinePart_curve);
% 存储排序后的片段，以及对应的工况曲线
kinePartSorted{i} = kinePart_g1_sorted;
curve{i} = curve1;
end
% 画三类工况曲线图
fprintf('三类工况曲线计算完成！\n');
groupName = {'第一类工况曲线','第二类工况曲线','第三类工况曲线'};
for i = 1:length(groupName)
figure;
set(gcf,'outerposition',get(0,'screensize'))
plot(curve{i},'r-','linewidth',2)
xlabel('时间（s）'),ylabel('车速（km/h）')
title(groupName{i});
grid on
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24);
axis tight
print(gcf,'-djpeg','-r300',groupName{i});
fprintf('%s绘制完成！\n',groupName{i})
end
end
此函数负责按照聚类结果，进行相关的总体分布特征指标计算，并按照相似度挑选片段，合成工况。
- 总体分布特征指标
function [distFeature,distFeature_mat ]= calDistFeature(kinePart)
% 定义分存储分布特征的结构体
distFeature = struct(...
'neglect_time',{}, ... % 怠速时间比
'accelerate_time',{}, ... % 加速时间比
'slowdown_time',{},... % 减速时间比
'unispeed_time',{},... % 匀速时间比
'v0_10_time',{},... % 0-10km 速度段比例
'v10_20_time',{},... % 10-20km 速度段比例
'v20_30_time',{},... % 20-30km 速度段比例
'v30_40_time',{},... % 30-40km 速度段比例
'v40_50_time',{},... % 40-50km 速度段比例
'v50_60_time',{},... % 50-60km 速度段比例
'v60_70_time',{},... % 60-70km 速度段比例
'v70_inf_time',{}); % 70-infkm 速度段比例
m = length(kinePart);
for i = 1:m
% 提取速度
data_temp = kinePart{i};
v_delta = (data_temp(2:end) - data_temp(1:end-1)); % 计算速度差
% 计算运动段总体分布特征
distFeature(i).neglect_time = sum(v_delta == 0)/(length( data_temp)-1); % 计算怠速时间比
distFeature(i).accelerate_time = sum( v_delta > 0.1*3600/1000)/...
(length( data_temp)-1); % 计算加速时间比
distFeature(i).slowdown_time = sum( v_delta < -0.1*3600/1000)/...
(length( data_temp)-1); % 计算减速时间比
distFeature(i).unispeed_time = sum( v_delta >= -0.1*3600/1000 & ...
v_delta <= 0.1*3600/1000 & v_delta ~= 0)/(length( data_temp)-1); % 计算减速时间比
distFeature(i).v0_10_time = sum(data_temp>0 & data_temp<=10)/...
(sum( data_temp>0)); % 0-10km 速度段比例
distFeature(i).v10_20_time = sum(data_temp>10 & data_temp<=20)/...
(sum( data_temp>0)); % 10-20km 速度段比例
distFeature(i).v20_30_time = sum(data_temp>20 & data_temp<=30)/...
(sum( data_temp>0)); % 20-30km 速度段比例
distFeature(i).v30_40_time = sum(data_temp>30 & data_temp<=40)/...
(sum( data_temp>0)); % 30-40km 速度段比例
distFeature(i).v40_50_time = sum(data_temp>40 & data_temp<=50)/...
(sum( data_temp>0)); % 40-50km 速度段比例
distFeature(i).v50_60_time = sum(data_temp>50 & data_temp<=60)/...
(sum( data_temp>0)); % 50-60km 速度段比例
distFeature(i).v60_70_time = sum(data_temp>60 & data_temp<=70)/...
(sum( data_temp>0)); % 60-70km 速度段比例
distFeature(i).v70_inf_time = sum(data_temp>70)/(sum( data_temp)-1); % 70-infkm 速度段比例
end
distFeature_mat = cell2mat(reshape(struct2cell(distFeature),12,length(distFeature))'); % 总体分布特征指标转矩阵存储
end
此函数负责计算运动片段的总体分布特征指标，这里的指标主要包括各工况持续时间占比，以及速度区间占比
- 皮尔逊相关系数计算
function coeff = myPearson(X , Y)
% 本函数实现了皮尔逊相关系数的计算操作
%
% 输入：
% X：输入的数值序列
% Y：输入的数值序列
%
% 输出：
% coeff：两个输入数值序列X，Y的相关系数
%
if length(X) ~= length(Y)
error('两个数值数列的维数不相等');
return;
end
fenzi = sum(X .* Y) - (sum(X) * sum(Y)) / length(X);
fenmu = sqrt((sum(X .^2) - sum(X)^2 / length(X)) * (sum(Y .^2) - sum(Y)^2 / length(X)));
coeff = fenzi / fenmu;
end
此函数负责计算皮尔逊相关系数，源码来自网络
- 片段挑选函数
function is_choose = choosePart(part_time,time_min,time_max)
m = length(part_time);
is_choose = zeros(m,1);
for i = 1:m
is_choose(i) = 1;
if sum(part_time(is_choose==1)) < time_min
is_choose(i) = 1;
elseif sum(part_time(is_choose==1)) > time_min && sum(part_time(is_choose==1)) < time_max
break
else
is_choose(i) = 0;
end
end
此函数能够实现按照一定的时间长度区间，挑选片段。
- 计算合成工况曲线参数
function [sectionSpeed,sectionTime] = calSection(curve,kinePart)
% 分析三类工况曲线
curveDistFeature = calDistFeature(curve);
curvebardata = [...
curveDistFeature.neglect_time;...
curveDistFeature.accelerate_time;...
curveDistFeature.slowdown_time;...
curveDistFeature.unispeed_time,...
];
curvebarname = {'怠速模式','加速模式','减速模式','匀速模式'};
figure;set(gcf,'outerposition',get(0,'screensize'));
groupName = {'第一类工况曲线','第二类工况曲线','第三类工况曲线'};
bar(curvebardata),legend(groupName),grid on;
set(gca,'xticklabel',curvebarname,'fontsize',24,'GridLineStyle','--','GridColor','k','GridAlpha',1);
print(gcf,'-djpeg','-r300','三类工况曲线的运行模式对比');
% 计算速度区间
sectionSpeed = calSpeedSection(curveDistFeature);
% 计算速度区间持续长度
sectionTime = calSectionTime(sectionSpeed,kinePart);
end
此函数负责计算合成工况曲线的速度区间，以及对应的片段持续时长，并给出工况分析柱形图。
- 计算速度区间
function speedSection = calSpeedSection(curveDistFeature)
curvespeeddata = [...
curveDistFeature.v0_10_time;...
curveDistFeature.v10_20_time;...
curveDistFeature.v20_30_time;...
curveDistFeature.v30_40_time;...
curveDistFeature.v40_50_time;...
curveDistFeature.v50_60_time;...
curveDistFeature.v60_70_time;...
curveDistFeature.v70_inf_time;...
];
index = zeros(3,1);
for i = 1:3
indextemp = find(curvespeeddata(:,i)==0);
if isempty(indextemp)
index(i) = inf;
else
index(i) = indextemp(1)*10-10;
end
end
index_sorted = sort(index);
index_sorted = [0;index_sorted];
for i = 1:length(index)
ind = find(index_sorted==index(i));
speedSection(i,:) = [index_sorted(ind-1),index_sorted(ind)];
end
end
此函数实现根据总体分布特征指标，计算速度分布区间。
- 计算速度区间持续时长
function sectionTime = calSectionTime(sectionSpeed,kinePart)
kinePart_mat = cell2mat(kinePart);
speed_s = zeros(3,1);
for i = 1:length(kinePart_mat)
for j = 1:3
if kinePart_mat(i) > sectionSpeed(j,1) && kinePart_mat(i) <= sectionSpeed(j,2)
speed_s(j) = speed_s(j) +1;
end
end
end
time_s1 = (speed_s(1)/sum(speed_s))*[1200,1300];
time_s2 = (speed_s(2)/sum(speed_s))*[1200,1300];
time_s3 = (speed_s(3)/sum(speed_s))*[1200,1300];
sectionTime = [time_s1;time_s2;time_s3];
此函数负责按照速度区间计算合成工况曲线各各类工况片段的持续时长。
- 计算组合工况曲线
function curveComposed = calCurveComposed(kinePartSorted,sectionTime,sectionSpeed)
for i = 1:3
% 计算片段时长
kinePart_now = kinePartSorted{i};
m = length( kinePart_now);
part_time = zeros(m,1);
for j = 1:m
part_time(j) = length( kinePart_now{j});
end
% 选取片段
time_min = sectionTime(i,1);time_max= sectionTime(i,2);
is_choose = choosePart(part_time,time_min,time_max);
kinePart_curve = kinePart_now(is_choose==1,1);
% 得到片段工况曲线
curve_temp{i,1}= kinePart_curve;
end
[~,index]= sort(sectionSpeed(:,1));
curveComposed = [...
cell2mat( curve_temp{index(1),1});...
cell2mat( curve_temp{index(2),1});...
cell2mat( curve_temp{index(3),1});...
];
% 画合成工况曲线图
figure;
set(gcf,'outerposition',get(0,'screensize'))
plot(curveComposed,'r-','linewidth',2)
xlabel('时间（s）'),ylabel('车速（km/h）')
title('合成工况曲线');
grid on
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24);
axis tight
print(gcf,'-djpeg','-r300','合成工况曲线');
fprintf('%s绘制完成！\n','合成工况曲线')
end
此函数实现了按照速度区间和速度区间持续时长去选取片段并组合成合成工况曲线。
- 联合特征指标误差分析与计算
function [erro_abs,erro_rel,erro_mean] = calCombinFeatureErro(kinePart,curve,curveComposed)
% 计算所有运动学片段的联合特征指标
[~,combinFeature_mat] = calCombinFeature (kinePart);
combinFeature_all = mean(combinFeature_mat); % 取均值
% 计算所有工况曲线联合特征指标
curve_all = [curve,{curveComposed}];
[~,combinFeature_curve] = calCombinFeature(curve_all);
% 计算绝对误差
erro_abs = combinFeature_curve - combinFeature_all;
erro_rel = abs(erro_abs)./abs(combinFeature_all);
erro_mean = mean(erro_rel');
name = {'第一类工况曲线','第二类工况曲线','第三类工况曲线','合成工况曲线'};
[~,index]= min(erro_mean);
figure;set(gcf,'outerposition',get(0,'screensize'))
subplot(2,1,1)
plot(1:11,combinFeature_all,'ro-',1:11,combinFeature_curve(index,:),'bo-','linewidth',2);
xlabel('评价指标'),ylabel('指标取值'); legend({'实际工况',[name{index},'(误差最小工况)']})
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24);
grid on ,axis tight
subplot(2,1,2)
plot(1:11,erro_rel(index,:),'bo-','linewidth',2);
xlabel('评价指标'),ylabel('相对误差');
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24);
grid on ,axis tight
print(gcf,'-djpeg','-r300','联合特征评价指标误差分析');
fprintf('%s绘制完成！\n','联合特征评价指标误差分析')
end
此函数按照新的联合特征指标重新评价构建的工况曲线和原始数据，并计算其对应的绝对误差、相对误差、平均相对误差。绘制相应的工况误差曲线。
- 新的联合特征指标
function [combinFeature,combinFeature_mat] = calCombinFeature (kinepart)
% 定义存储特征值的结构体
combinFeature = struct(...
'av_speed',{},... % 平均速度
'av_runspeed',{},... % 平均行驶速度
... % 'max_speed',{},... % 最大速度
'speed_std',{},... % 速度标准差
'av_accelerate',{},... % 平均加速度
'accelerate_std',{},... % 加速度标准差
'av_slowdown',{},... % 平均减速度
'slowdown_std',{},... % 减速度标准差
'neglect_time',{},... % 怠速时间
'accelerate_time',{},... % 加速时间
'slowdown_time',{},... % 减速时间
'unispeed_time',{}); % 匀速时间
m = length(kinepart);
for i = 1:m
% 提取速度
data_temp = kinepart{i};
v_delta = (data_temp(2:end) - data_temp(1:end-1)); % 计算速度差
combinFeature(i).av_speed = mean(data_temp); % 计算平均速度
combinFeature(i).av_runspeed = mean(data_temp(data_temp~=0)); % 计算平均行驶速度
% combinFeature(i).max_speed = max(data_temp); % 计算最大速度
combinFeature(i).speed_std = std(data_temp); % 计算速度标准差
% 加速度特征值 单位：m/s^2
v_delta2 = v_delta*1000/3600; %单位换算
combinFeature(i).av_accelerate = mean(v_delta2(v_delta2 > 0.1)); % 计算平均加速度
combinFeature(i).accelerate_std = std(v_delta2(v_delta2 > 0.1)); % 计算加速度标准差
combinFeature(i).av_slowdown = mean(v_delta2(v_delta2 < -0.1)); % 计算平均减速度
combinFeature(i).slowdown_std = std(v_delta2(v_delta2 < -0.1)); % 计算减速度标准差
combinFeature(i).neglect_time = sum(v_delta == 0)/(length( data_temp)-1);% 计算怠速时间比
combinFeature(i).accelerate_time = sum( v_delta > 0.1*3600/1000)/...
(length( data_temp)-1); % 计算加速时间比
combinFeature(i).slowdown_time = sum( v_delta < -0.1*3600/1000)/...
(length( data_temp)-1); % 计算减速时间比
combinFeature(i).unispeed_time = sum( v_delta >= -0.1*3600/1000 & ...
v_delta <= 0.1*3600/1000 & v_delta ~= 0)/(length( data_temp)-1); % 计算匀速时间比
end
combinFeature_mat = cell2mat(reshape(struct2cell(combinFeature),11,length(combinFeature))'); % 总体分布特征指标转矩阵存储
end
此函数重新定义了一些特征指标，用来评估构建的工况曲线和原始数据。
- 速度-加速度联合概率分布误差分析
function [erro_abs,erro_rel,erro_mean] = calVAcombineProErro(kinePart,curve,curveComposed)
kinePart_all = cell2mat(kinePart);
kinePart_all_VAcombinePro = calVAcombinePro(kinePart_all);
curve1_VAcombinePro = calVAcombinePro(curve{1});
curve2_VAcombinePro = calVAcombinePro(curve{2});
curve3_VAcombinePro = calVAcombinePro(curve{3});
curveComposed_VAcombinePro = calVAcombinePro(curveComposed);
curve_VAcombinePro = {...
curve1_VAcombinePro,...
curve2_VAcombinePro,...
curve3_VAcombinePro,...
curveComposed_VAcombinePro
};
% 计算绝对误差
for i = 1:4
erro_abs{i} = curve_VAcombinePro{i} - kinePart_all_VAcombinePro;
erro_rel_temp = abs(curve_VAcombinePro{i} - kinePart_all_VAcombinePro)./...
abs(kinePart_all_VAcombinePro);
erro_rel_temp(isnan(erro_rel_temp)) = 0;
erro_rel{i} = erro_rel_temp;
erro_mean(i) = mean(mean(erro_rel{i}));
end
[~,index]= min(erro_mean);
figure;
set(gcf,'outerposition',get(0,'screensize'))
groupName = {'第一类工况','第二类工况','第三类工况','合成工况'};
for i = 1:4
subplot(2,2,i)
surf(curve_VAcombinePro{i})
xlabel('加速度(m/s^2)'),ylabel('速度(km/h)'),zlabel('概率');
if i == index
title([groupName{i},'（误差最小工况）']);
else
title(groupName{i});
end
grid on
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24,...
'xtick',1:2:8,'xticklabel',-3:2:3,'ytick',1:2:8,'yticklabel',0:20:80);
axis tight
end
print(gcf,'-djpeg','-r300','各类工况速度-加速度概率分布图');
figure;
set(gcf,'outerposition',get(0,'screensize'))
surf(kinePart_all_VAcombinePro)
xlabel('加速度(m/s^2)'),ylabel('速度(km/h)'),zlabel('概率');
title('原始工况');
grid on
set(gca,'GridLineStyle','--','GridColor','k','GridAlpha',1,'fontsize',24,...
'xtick',1:8,'xticklabel',-3:1:3,'ytick',1:8,'yticklabel',0:10:80);
axis tight
print(gcf,'-djpeg','-r300','原始工况速度-加速度概率分布图');
fprintf('%s绘制完成！\n','速度-加速度联合概率分布误差分析')
end
此函数运用速度-加速度联合概率分布计算构建的工况曲线和原始数据，绘制出对应的分布图，以此来评估构建的工况拟合效果。
- 速度-加速度联合概率分布计算
function VAcombinePro = calVAcombinePro(kinePart_all)
VAcombinePro = zeros(8,8);
Speed = kinePart_all(2:end);
Accelerate = (kinePart_all(2:end) - kinePart_all(1:end-1))*1000/3600; %单位换算
m = length(Speed);
for i = 1:m
row = findAccelerateSection(Accelerate(i));
col = findSpeedSection(Speed(i));
VAcombinePro(row,col) = VAcombinePro(row,col) + 1;
end
VAcombinePro = VAcombinePro/m;
end
</pre>
此函数实现了对一条工况曲线进行速度-加速度联合概率分布的计算。
- 加速度分布区间计算
<pre class="matlab-code">
function Accelerate_index = findAccelerateSection(Accelerate)
if Accelerate < -3
Accelerate_index = 1;
elseif Accelerate >=-3 && Accelerate < -2
Accelerate_index = 2;
elseif Accelerate >=-2 && Accelerate < -1
Accelerate_index = 3;
elseif Accelerate >=-1 && Accelerate < 0
Accelerate_index = 4;
elseif Accelerate >=0 && Accelerate < 1
Accelerate_index = 5;
elseif Accelerate >=1 && Accelerate < 2
Accelerate_index = 6;
elseif Accelerate >=2 && Accelerate < 3
Accelerate_index = 7;
elseif Accelerate >=3
Accelerate_index = 8;
end
end
此函数实现对于任意一个加速度，计算其对应分布区间。
- 速度分布区间计算
function Speed_index = findSpeedSection(speed)
if speed >=0 && speed < 10
Speed_index = 1;
elseif speed >=10 && speed < 20
Speed_index = 2;
elseif speed >=20 && speed < 30
Speed_index = 3;
elseif speed >=30 && speed < 40
Speed_index = 4;
elseif speed >=40 && speed < 50
Speed_index = 5;
elseif speed >=50 && speed < 60
Speed_index = 6;
elseif speed >=60 && speed < 70
Speed_index = 7;
elseif speed >=70
Speed_index = 8;
end
end
此函数实现对于任意一个速度，计算其对应分布区间。
- 主函数
% 准备存储空间
clc,clear,close all
% 合并所有文件运动学片段
kinePart = kinePartCombine({'文件1','文件2','文件3'});
% 计算运动学特征指标
[kineFeature,kineFeature_mat] = calKineFeature(kinePart);
% 主成分分析降维(按累计贡献率85%以上)
[kineFeature_dim,mylatent] = myPCA(kineFeature_mat,0.85);
% K-Means聚类，聚成三类
group =myK_means(kineFeature_dim,3);
% 按聚类结果，计算对应工况曲线
[kinePartSorted,curve] = calCurve2(kinePart,group);
% 计算合成曲线的片段参数
[sectionSpeed,sectionTime] = calSection(curve,kinePart);
% 合成工况曲线
curveComposed = calCurveComposed(kinePartSorted,sectionTime,sectionSpeed);
% 联合特征指标误差分析
[CombinFeature_erro_abs,CombinFeature_erro_rel,CombinFeature_erro_mean] = ...
calCombinFeatureErro(kinePart,curve,curveComposed);
% 速度-加速度联合概率分布误差分析
[VAcombinePro_erro_abs,VAcombinePro_erro_rel,VAcombinePro_erro_mean] = ...
calVAcombineProErro(kinePart,curve,curveComposed);
% 决定最终的工况曲线
curvename = {'第一类工况曲线','第二类工况曲线','第三类工况曲线','合成工况曲线'};
[~,index] = min(CombinFeature_erro_mean+VAcombinePro_erro_mean);
fprintf('经过两种误差分析与计算，最佳的工况曲线为%s\n',curvename{index});
主程序按照绘制的流程图一步步计算，并得出最终的最佳工况曲线。
运行结果
代码运行时，画面是这样的：
命令行窗口文字
正在导入文件1运动学片段...
正在导入文件2运动学片段...
正在导入文件3运动学片段...
所有文件运动学片段提取合并完成！
运动学特征指标计算完成！
主成分分析降维计算完成！
K-Means聚类计算完成！
三类工况曲线计算完成！
第一类工况曲线绘制完成！
第二类工况曲线绘制完成！
第三类工况曲线绘制完成！
合成工况曲线绘制完成！
联合特征评价指标误差分析绘制完成！
速度-加速度联合概率分布误差分析绘制完成！
经过两种误差分析与计算，最佳的工况曲线为第三类工况曲线
>>
运行结束后，将会得到这样几幅图：
- 联合特征指标相对误差分析表
| 第一类工况（低速工况） | 第二类工况（高速工况） | 第三类工况（中速工况） | 合成工况 | |
|---|---|---|---|---|
| 平均速度 | 0.682461 | 1.303494 | 0.093666 | 0.026168 |
| 平均行驶速度 | 0.594824 | 0.942352 | 0.148446 | 0.122579 |
| 速度标准差 | 0.486327 | 1.077564 | 0.33433 | 1.000774 |
| 平均加速度 | 0.047345 | 0.111123 | 0.041569 | 0.008425 |
| 加速度标准差 | 0.018052 | 0.135695 | 0.124836 | 0.051627 |
| 平均减速度 | 0.133704 | 0.062285 | 0.126807 | 0.08847 |
| 减速度标准差 | 0.151064 | 0.230756 | 0.198747 | 0.275743 |
| 怠速时间比 | 0.330738 | 0.471477 | 0.028332 | 0.166958 |
| 加速时间比 | 0.217645 | 0.281698 | 0.138994 | 0.055789 |
| 减速时间比 | 0.129899 | 0.215865 | 0.026914 | 0.117759 |
| 匀速时间比 | 0.158491 | 0.232193 | 0.202386 | 0.104164 |
| 平均误差 | 0.268232 | 0.460409 | 0.133184 | 0.183496 |
其中
- 低速工况速度区间：0-30 km/h ；
- 中速工况速度区间：30-50 km/h；
- 高速工况速度区间：50km/h 以上。
速度-加速度联合概率分布平均误差：
| 第一类工况（低速工况） | 第二类工况（高速工况） | 第三类工况（中速工况） | 合成工况 | |
|---|---|---|---|---|
| 平均相对误差 | 0.889031 | 6.839725 | 0.992904 | 1.145535 |
通过两种工况曲线评价方式，综合比较，得到结论为：
所采集的数据适合构建成行驶速度为30-50km/h的中速行驶工况。其工况曲线图如下：
参考文献
源码分享
整个题目三问的代码整理上传到了Github，点击查看
写在最后
对于建模的总结
其实，在整理代码过程中，我也发现了在第三问上，我们的做法存在一定的局限性。很明显，我们的模型误差在15%左右，相比于文献中的结果还是不够理想。不管怎么说，在短短的五天时间里，能够复现他人文献的内容，也算是一种能力。相信在后面的学习过程中，对于问题的宏观把控能力和关键细节信息的挖掘能力都能慢慢得到提高。
对于编程的总结
对于这次比赛的题目，其实要说模型的储备，也就只用到了主成分分析和聚类分析，其余的都得靠自己去想办法理解实现过程，再来编程，没有任何参考。这对我是一个非常大的考验。
通过这一个暑假的建模，我个人觉得，思维导图和流程图是分析问题的利器。思维导图能够细化知识点，让你的脑子能够聚焦问题的关键的同时，又能够对整体的框架有个宏观的把控。而流程图能够分析出，在设计函数时的输入和输出关系，以及各个模块之间的关联，所以要想编程思路清晰，必须对流程清晰，而流程的理解，画出来那就是流程图。
未经作者授权，禁止转载
THE END
