基于matlab实现马尔科夫链的工况预测，包含多个标准工况的数据

clear all; close all; clc; format compact; %% unit conversion mph2mps = 1609/3600; mps2mph = 1/mph2mps; %% velocity data processing cd('CYC'); %% NEDC工况数据读取 % load nedc.mat; % cyc_name='nedc'; % cyc_time = 1:length(sch_cycle); % cyc_mps = sch_cycle(cyc_time,2); %% FUDDS工况数据读取 load CYC_FUDS.mat; cyc_name = 'FUDS'; cyc_mps = cyc_mph(:,2)*mph2mps; %% US06工况数据读取 % load CYC_US06.mat; % cyc_name = 'US06'; % cyc_mps = cyc_mph(:,2)*mph2mps; %% 工况数据拼接 % cyc_name='combine_cyc'; % load nedc.mat;cyc_time = 1:length(sch_cycle);cyc1_mps=sch_cycle(cyc_time,2); % load CYC_FUDS.mat;cyc2_mps=cyc_mph(:,2)*mph2mps; % load CYC_US06.mat; cyc3_mps=cyc_mph(:,2)*mph2mps; % load CYC_HWFET.mat;cyc4_mps=cyc_mph(:,2)*mph2mps; % load test_cycle_ag.mat;cyc5_mps=cyc_mph(:,2)*mph2mps; % load test_cycle_ag_rb.mat;cyc6_mps=cyc_mph(:,2)*mph2mps; % load test_cycle_agx2_idbetween.mat;cyc7_mps=cyc_mph(:,2)*mph2mps; % load test_cycle_const_rb.mat;cyc8_mps=cyc_mph(:,2)*mph2mps; % load test_cycle_constant_spd.mat;cyc9_mps=cyc_mph(:,2)*mph2mps; % %此处还可以继续添加工况信息！ % cyc_mps=[cyc1_mps;cyc2_mps;cyc3_mps;cyc4_mps;cyc5_mps;cyc6_mps;cyc7_mps;cyc8_mps;cyc9_mps]; %% 工况数据初步处理 cyc_mps(cyc_mps<0) = 0; mps=cyc_mps'; acc_row=[diff(mps) 0]; acc=acc_row'; cd .. %% data statistics %% velocity data statistics cyc_mps_rnd=roundn(cyc_mps,-1); %对速度值的小数位数进行限制，保证速度值可以相同。 v_sorted=unique(cyc_mps_rnd,'rows'); %对近似后的速度值进行重复删减，并将其按大小顺序排列。 v_no=length(v_sorted); v=zeros(v_no,2); v(:,1)=v_sorted'; cyc_mps_rnd(1,1)=1; %做标记，用于加速度统计。 for i=1:length(cyc_mps_rnd) %速度归类，相同的速度归为一组。外循环用于循环原速度值。 for j=1:length(v_sorted) %内循环表示已经分组的速度值。 if cyc_mps_rnd(i)==v(j,1) %如果循环到的原速度值与已经分组的速度值相等，则加入该分组。 v(j,2)=v(j,2)+1; cyc_mps_rnd(i,2)=j; %做标记，用于加速度统计。 end end end %% accelaration data statistics acc_rnd=roundn(acc,-1); %对加速度值的小数位数进行限制，保证加速度值可以相同。 v_acc_joint=[cyc_mps_rnd acc_rnd]; acc_sorted=unique(acc_rnd,'rows'); acc_no=length(acc_sorted); v_acc=zeros(v_no,acc_no); for m=1:length(cyc_mps_rnd) for r=1:v_no if v_acc_joint(m,2)==r for n=1:acc_no if v_acc_joint(m,3)==acc_sorted(n) v_acc(r,n)=v_acc(r,n)+1; end end end end end %% calculate the transitional matrix v_sum_update_row=sum(v_acc'); v_sum_update=v_sum_update_row'; for i=1:length(v_sum_update_row) if isnan(v_sum_update(i)) v_sum_update(i)=1; end end for i=1:v_no for j=1:acc_no v_acc(i,j)=v_acc(i,j)/v_sum_update(i); end end %% predefine the conditions to obtain the transitional value trans_val=zeros(v_no,1); trans_p=zeros(v_no,acc_no); for i=1:v_no %用于循环联合概率分布矩阵中每一速度行的的概率值 for j=1:acc_no if v_acc(i,j)>=0.8 trans_p(i,j)=1; %概率大于限定值直接取此概率对应的加速度作为下一时刻的预测值 trans_val(i)=acc_sorted(j); end end end for i=1:v_no if max(trans_p(i,:))==0 trans_val(i)=v_acc(i,:)*acc_sorted; %概率值均小于限定值则取加权值作为下一时刻的预测值 end end %% output the statistics and graphs x_max=max(v_sorted); x_min=min(v_sorted); y_max=max(acc_sorted); y_min=min(acc_sorted); % figure(1) % xlim([x_min x_max]); % ylim([y_min y_max]); % mesh(v_acc); % xlabel('速度（m/s）'); % ylabel('加速度（m/s^2）'); % zlabel('概率'); %% make sure that there are no zero elements in the arrow trans_val for i=2:(v_no-1) if trans_val(i)==0 trans_val(i)=(trans_val(i-1)+trans_val(i+1))/2; %为避免预测成直线，将转移数据中所有的零数据均处理为前一个值和后一个值的平均 end end if trans_val(1)==0 trans_val(1)=trans_val(2)/2; end if trans_val(v_no)==0 trans_val(v_no)=trans_val(v_no-1)/2; end % trans_val(14)=abs(trans_val(14)); % trans_val(15)=abs(trans_val(15)); % trans_val(16)=abs(trans_val(16)); % trans_val(17)=abs(trans_val(17)); % trans_val(18)=abs(trans_val(18)); %% prediction demo %% load data cd('CYC'); load nedc.mat; % load CYC_FUDS.mat; % load CYC_US06.mat; %% get the first data to carry on the next prediction step % cyc_sample=sch_cycle(cyc_time,2); cyc_sample= cyc_mph(:,2)*mph2mps; v_0=cyc_sample(1); %% start prediction figure(2) v_prediction=zeros(2000,1); plot(v_sorted,trans_val); xlabel('速度（m/s）'); ylabel('加速度（m/s^2）'); acc_0=interp1(v_sorted,trans_val,v_0,'linear'); v_prediction(1)=v_0; v_prediction(2)=v_0+acc_0; for i=2:1200 %define the length of the driving cycle acc_predict=interp1(v_sorted,trans_val,v_prediction(i),'linear'); if acc_predict~=0 v_prediction(i+1)=v_prediction(i)+acc_predict; else v_prediction(i+1)=v_prediction(i)+0.8; end end %% 强制给定一个高速以使其在高速工况之后继续预测 % v_prediction(81)=9; % for i=82:1200 %define the length of the driving cycle % % acc_predict=interp1(v_sorted,trans_val,v_prediction(i),'linear'); % v_prediction(i+1)=v_prediction(i)+acc_predict; % % % end markov_demo(); %% compare the result to the original data figure(3) plot(v_prediction/mph2mps); xlabel('时间（s）'); ylabel('速度（m/s）'); xlim([0 1200]); hold on; plot(cyc_mps,'r:'); legend('预测工况','标准工况');