clc;
clear all;
format long;
% 数据加载
%nedc = cell2mat(struct2cell(load('D:\matlab2022\nedc.mat')));
%nedc = cell2mat(struct2cell(load('D:\画图指标数据\CLTC-WLTC-NEDC-CHTC等七种主流工况数据MATLAB格式\cycles_chtc_ht.mat')));
cltc= readtable('D:/matlab2022/降噪CLTC-HT.xlsx');
nedc= table2array(cltc(:, 1));
time=table2array(cltc(:,2));
x_tot=trapz(time,nedc./3.6)/1000;
data = readtable('C:/Users/10058061/Desktop/发动机map2.xlsx');
engine_RPM_Axis = table2array(data(:, 1));
Ttq = table2array(data(:, 2));
engine_Torque_max = [2423.5 2430.6 2437.8 2434.2 2421.8 2429.8 2375.2 2317.4 2279.8 2251.2];
engine_RPM_Axis_max = [800 900 1000 1100 1200 1300 1400 1500 1600 1700];
mf_kn = table2array(data(:, 3));
engine_be_Data = table2array(data(:, 6));
% 确保mf_kn是二维矩阵
if ~ismatrix(mf_kn) || size(mf_kn, 2) == 1
mf_kn = reshape(mf_kn, length(engine_RPM_Axis), []);
end
% 整车参数
igt = [12.26, 9.556, 7.356, 5.734, 4.521, 3.524, 2.712, 2.114, 1.627, 1.268, 1, 0.779];
i0 = 3.083;
Rw = 0.522;
Mass = 8200;
Cx = 0.6;
S = 9.09;
f = 0.015;
g = 9.8;
det = 1.05;
Jw = 14.87;
itatgb = 0.96;
itatfd = 0.95;
itat = 0.912;
n_engine = [800, 1700];
shift = 20;
% 计算加速度和车轮转矩
K = length(nedc);
acc = zeros(K, 1);
Twheel = zeros(K, 1);
for kk = 1:K - 1
acc(kk + 1) = (nedc(kk + 1) - nedc(kk))/2;
end
for kk = 1:K - 1
if nedc(kk) == 0 && acc(kk) == 0
Twheel(kk) = 0;
else
Twheel(kk) = (Mass * g * f + 0.5 * 1.225 * Cx * S * (nedc(kk)/3.6)^2 + det * Mass * acc(kk)) * Rw;
end
end
% 确定车辆行驶状态变化点
S_r = [];
nnn = 1;
for kk = K:-1:2
if nedc(kk) == 0 && acc(kk) < 0 && nedc(kk - 1) > 0
S_r(nnn) = kk;
nnn = nnn + 1;
end
end
S_r = fliplr(S_r);
S_1 = [];
mmm = 1;
for kk = 2:K
if nedc(kk) > 0 && nedc(kk - 1) == 0 && acc(kk) > 0
S_1(mmm) = kk;
mmm = mmm + 1;
end
end
S_rr = [];
www = 1;
for hh = 1:length(S_1)
for tt = S_r(hh):-1:S_1(hh)
if Twheel(tt) > 0
S_rr(www) = tt;
break;
end
end
www = www + 1;
end
% 初步筛选适合的档位
igt_fea = zeros(length(igt), length(S_rr));
for h = 1:length(S_rr)
N_ice = 30 * nedc(S_rr(h)) * igt * i0 / (pi * Rw);
for pp = 1:12
if N_ice(pp) <= 1700 && N_ice(pp) >= 800
igt_fea(pp, h) = pp;
end
end
end
% 根据转矩限制进一步筛选档位
S_rr_igt = zeros(length(igt), length(S_rr));
for hh = 1:length(S_rr)
for pp = 1:12
if igt_fea(pp, hh) ~= 0
nee = 30 * nedc(S_rr(hh)) * igt(igt_fea(pp, hh)) * i0 / (pi * Rw);
te = Twheel(S_rr(hh)) / (i0 * igt(igt_fea(pp, hh)) * itat);
temax = interp1(engine_RPM_Axis_max, engine_Torque_max, nee, 'cubic');
if te <= temax
S_rr_igt(pp, hh) = igt_fea(pp, hh);
end
end
end
end
% 计算每个档位下的燃油消耗（改进版）
S_rr_fuel = ones(length(igt), length(S_rr)) * 100000;
for hh = 1:length(S_rr)
for pp = 1:12
if S_rr_igt(pp, hh) == 0
continue;
end
gear = S_rr_igt(pp, hh);
nee = 30 * nedc(S_rr(hh)) * igt(gear) * i0 / (pi * Rw);
te = Twheel(S_rr(hh)) / (i0 * igt(gear) * itat);
ne_clamped = max(min(nee, max(engine_RPM_Axis)), min(engine_RPM_Axis));
te_clamped = max(min(te, max(Ttq)), min(Ttq));
try
sfuel = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'cubic');
catch
sfuel = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'linear');
end
if isnan(sfuel)
[~, ne_idx] = min(abs(engine_RPM_Axis - ne_clamped));
[~, te_idx] = min(abs(Ttq - te_clamped));
ne_idx = max(1, min(ne_idx, size(mf_kn, 1)));
te_idx = max(1, min(te_idx, size(mf_kn, 2)));
sfuel = mf_kn(ne_idx, te_idx);
end
S_rr_fuel(pp, hh) = sfuel;
end
end
% 选择每个阶段的最佳档位
S_rr_f = zeros(length(S_rr), 1);
for hh = 1:length(S_rr)
[~, y2] = min(S_rr_fuel(:, hh));
S_rr_f(hh) = y2;
end
% 动态规划算法确定最佳换挡策略
ne_final = zeros(length(nedc), 1);
te_final = zeros(length(nedc), 1);
gear_final = zeros(length(nedc), 1);
fuel_final = zeros(length(nedc), 1);
N = length(igt);
f_n = zeros(length(S_1), 1);
for hh = 1:length(S_1)
f_n(hh) = S_rr(hh) - S_1(hh) + 1;
end
for hh = 1:length(S_1)
Fuel = zeros(N, f_n(hh));
for ii = 1:N
if ii == S_rr_f(hh)
Fuel(ii, f_n(hh)) = 0.05;
else
Fuel(ii, f_n(hh)) = 10900;
end
end
Gear_con = [-1, 0, 1];
Gear = zeros(N, f_n(hh));
for kk = (f_n(hh) - 1):-1:1
N_ice = 30 * nedc(kk - 1 + S_1(hh)) * igt * i0 / (pi * Rw);
tr_cost_1 = zeros(12, 1);
for ii = 1:12
if N_ice(ii) > 1699.99999 || N_ice(ii) < 799.9999
tr_cost_1(ii) = 10007;
else
tr_cost_1(ii) = 0.0001;
end
end
if tr_cost_1(1) == 10007
for ii = 1:12
if tr_cost_1(ii) >= 1699.99
x1 = ii;
Fuel(ii, kk) = tr_cost_1(ii);
Gear(ii, kk) = 17;
end
end
x1 = x1 + 1;
x2 = 12;
elseif tr_cost_1(1) == 0.0001 && tr_cost_1(12) == 0.0001
x1 = 1;
x2 = 12;
elseif tr_cost_1(12) == 10007
for ii = 12:-1:1
if tr_cost_1(ii) >= 1699.999
x2 = ii;
Gear(ii, kk) = 17;
Fuel(ii, kk) = tr_cost_1(ii);
end
end
x1 = 1;
x2 = x2 - 1;
else
x1 = 0;
x2 = 0;
end
if x1 > 0 && x2 > 0
Ne = 30 * nedc(kk - 1 + S_1(hh)) * igt(x2) * i0 / (pi * Rw);
Te = Twheel(kk - 1 + S_1(hh)) / (i0 * igt(x2) * itat);
Temax = interp1(engine_RPM_Axis_max, engine_Torque_max, Ne, 'cubic');
nn = 0;
while Te > Temax
x2 = x2 - 1;
Ne = 30 * nedc(kk - 1 + S_1(hh)) * igt(x2) * i0 / (pi * Rw);
Te = Twheel(kk - 1 + S_1(hh)) / (i0 * igt(x2) * itat);
Temax = interp1(engine_RPM_Axis_max, engine_Torque_max, Ne, 'cubic');
nn = nn + 1;
end
for ii = (x2 + nn):-1:(x2 + 1)
Fuel(ii, kk) = 20007;
Gear(ii, kk) = 777;
end
end
if Twheel(kk - 1 + S_1(hh)) >= 0.01
for ii = x1:x2
fuel = zeros(3, 1);
for jj = 1:3
if ii + Gear_con(jj) < 1 || ii + Gear_con(jj) > 12
fuel(jj) = Fuel(ii, kk + 1) + 10500;
else
v_current = nedc(kk - 1 + S_1(hh));
gear_candidate = ii + Gear_con(jj);
% 修正：使用连续的发动机转速范围计算车速范围
ne_range = n_engine(1):n_engine(2);
v_gear = 0.377 * ne_range * Rw / (igt(gear_candidate) * i0);
v_min = min(v_gear);
v_max = max(v_gear);
penalty = 0;
if v_current < v_min || v_current > v_max
penalty = 1000; % 车速超出档位合理范围的惩罚
end
if 30 * v_current * igt(gear_candidate) * i0 / (pi * Rw) < n_engine(1) || ...
30 * v_current * igt(gear_candidate) * i0 / (pi * Rw) > n_engine(2)
for pp = 1:12
if ii + Gear_con(jj) >= pp && ii + Gear_con(jj) <= pp + 1
fuel(jj) = Fuel(pp, kk + 1) + 10500 + penalty;
else
break;
end
end
else
ne = 30 * v_current * igt(gear_candidate) * i0 / (pi * Rw);
te = Twheel(kk - 1 + S_1(hh)) / (i0 * igt(gear_candidate) * itat);
ne_clamped = max(min(ne, max(engine_RPM_Axis)), min(engine_RPM_Axis));
te_clamped = max(min(te, max(Ttq)), min(Ttq));
try
tr_cost = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'cubic');
catch
tr_cost = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'linear');
end
if isnan(tr_cost)
[~, ne_idx] = min(abs(engine_RPM_Axis - ne_clamped));
[~, te_idx] = min(abs(engine_RPM_Axis - ne_clamped));
ne_idx = max(1, min(ne_idx, size(engine_be_Data, 1)));
te_idx = max(1, min(te_idx, size(engine_be_Data, 2)));
tr_cost = engine_be_Data(ne_idx, te_idx);
end
for pp = 1:12
if gear_candidate >= pp && gear_candidate < pp + 1
fuel(jj) = Fuel(pp, kk + 1) + tr_cost + penalty;
break;
end
end
end
end
end
[x22, x3] = min(fuel);
Fuel(ii, kk) = x22;
Gear(ii, kk) = Gear_con(x3);
end
else
for ii = 1:12
Fuel(ii, kk) = Fuel(ii, kk + 1);
Gear(ii, kk) = 0;
end
end
end
re_gear = zeros(f_n(hh), 1);
re_G = zeros(f_n(hh), 1);
re_Fuel = zeros(f_n(hh), 1);
re_ne = zeros(f_n(hh), 1);
re_te = zeros(f_n(hh), 1);
re_gear(1) = 1;
for ii = 1:12
if re_gear(1) == ii
re_G(1) = Gear(ii, 1);
re_G(1) = max(min(re_G(1), 1), -1);
break;
end
end
current_gear = min(max(re_gear(1) + re_G(1), 1), 12);
re_ne(1) = 30 * nedc(S_1(hh)) * igt(current_gear) * i0 / (pi * Rw);
re_te(1) = Twheel(S_1(hh)) / (i0 * igt(current_gear)*itat);
re_Fuel(1) = 0;
for kk = 2:f_n(hh)
current_gear = min(max(re_gear(kk - 1) + re_G(kk - 1), 1), 12);
re_ne(kk) = 30 * nedc(kk - 1 + S_1(hh)) * igt(current_gear) * i0 / (pi * Rw);
if Twheel(kk - 1 + S_1(hh)) > 0.001
re_te(kk) = Twheel(kk - 1 + S_1(hh)) / (i0 * igt(current_gear)*itat);
else
re_te(kk) = re_te(kk - 1);
end
ne_clamped = max(min(re_ne(kk), max(engine_RPM_Axis)), min(engine_RPM_Axis));
te_clamped = max(min(re_te(kk), max(Ttq)), min(Ttq));
try
re_fuel = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'cubic');
catch
re_fuel = griddata(engine_RPM_Axis, Ttq, engine_be_Data, ne_clamped, te_clamped, 'linear');
end
if isnan(re_fuel)
[~, ne_idx] = min(abs(engine_RPM_Axis - ne_clamped));
[~, te_idx] = min(abs(Ttq - te_clamped));
ne_idx = max(1, min(ne_idx, size(engine_be_Data, 1)));
te_idx = max(1, min(te_idx, size(engine_be_Data, 2)));
re_fuel = mf_kn(ne_idx, te_idx);
end
re_Fuel(kk) = re_Fuel(kk - 1) + re_fuel;
for pp = 1:12
if igt(current_gear) == igt(pp)
x4 = pp;
break;
end
end
x4 = min(max(x4, 1), 12);
re_G(kk) = Gear(x4, kk);
re_G(kk) = max(min(re_G(kk), 1), -1);
re_gear(kk) = current_gear;
end
for pp = S_1(hh):S_rr(hh)
idx = pp - S_1(hh) + 1;
if idx <= length(re_ne) && idx <= length(re_te) && idx <= length(re_gear) && idx <= length(re_Fuel)
ne_final(pp) = re_ne(idx);
te_final(pp) = re_te(idx);
gear_final(pp) = re_gear(idx);
fuel_final(pp) = re_Fuel(idx);
end
end
end
regions_to_smooth = [129,140;358,383;425,436;437,465;587,593;594,611;733,737;738,744;745,753;787,811;812,842;843,877;953,957;958,964;965,993;1176,1182;1390,1404];
% 保存原始数据的副本，用于回退
original_gear = gear_final;
original_ne = ne_final;
original_te = te_final;
for r = 1:size(regions_to_smooth, 1)
start_idx = regions_to_smooth(r, 1);
end_idx = regions_to_smooth(r, 2);
% 确保索引在有效范围内
if start_idx < 1 || end_idx > length(gear_final)
continue;
end
% 获取区域内的两个档位及出现次数
region_gears = gear_final(start_idx:end_idx);
unique_gears = unique(region_gears);
% 确认只有两个档位
if length(unique_gears) == 2
gear1 = unique_gears(1);
gear2 = unique_gears(2);
% 统计两个档位的出现次数
count1 = sum(region_gears == gear1);
count2 = sum(region_gears == gear2);
% 选择出现次数多的档位（若相同则选高档）
if count1 > count2
final_gear = gear1;
elseif count2 > count1
final_gear = gear2;
else
final_gear = max(gear1, gear2);
end
% 标记变量，记录区域是否被完全应用
region_applied = true;
% 尝试应用最终档位到整个区域
for i = start_idx:end_idx
new_ne = 30 * nedc(i) * igt(final_gear) * i0 / (pi * Rw);
% 放宽转速限制，允许短暂低于最低转速
if new_ne >= max(n_engine(1) - 100, 700) && new_ne <= n_engine(2)
gear_final(i) = final_gear;
ne_final(i) = new_ne;
te_final(i) = Twheel(i) / (i0 * igt(final_gear) * itat);
else
% 只要有一个点不满足条件，标记区域为未完全应用
region_applied = false;
end
end
% 如果区域未被完全应用，则整体回退到平滑档位，但使用原始转速和扭矩
if ~region_applied
for i = start_idx:end_idx
gear_final(i) = final_gear; % 应用平滑后的档位
ne_final(i) = original_ne(i); % 使用原始转速
te_final(i) = original_te(i); % 使用原始扭矩
end
end
end
end
figure(1);
time = (1:length(nedc));
plot(time,nedc);
title('NEDC工况车速');
xlabel('时间 (s)');
ylabel('车速 (km/h)');
grid on;
figure(2);
plot(gear_final);
title('改进的动态规划优化换挡策略');
xlabel('时间 (s)');
ylabel('挡位');
ylim([0 15]);
grid on;
figure(3);
plot(ne_final/1.5, te_final, 'b.', engine_RPM_Axis_max, engine_Torque_max, 'r-');
title('发动机工作点分布');
xlabel('发动机转速 (rpm)');
ylabel('发动机转矩 (N·m)');
grid on;
legend('优化工作点', '外特性曲线');
figure(4);
plot(ne_final/1.5);
title('转速变化');
xlabel('时间(s)');
ylabel('转速');
grid on;
figure(5);
plot(te_final);
title('扭矩变化');
xlabel('时间(s)');
ylabel('扭矩');
grid on;
fuel_final=diff(fuel_final);
% 计算总燃油消耗
% 计算总燃油消耗（排除fuel_final为0的点）
valid_idx = fuel_final > 0; % 找出fuel_final中不为0的索引
valid_time = time(valid_idx); % 获取对应的时间点
valid_fuel = fuel_final(valid_idx); % 获取对应的燃油消耗值
Fuel = trapz(valid_time, valid_fuel/3600); % 使用有效数据计算总燃油量
fuel_economy = Fuel/x_tot*100; % 计算百公里油耗
fprintf('总燃油消耗: %.2f L/100km\n', fuel_economy);将上述基于动态规划的换档策略改进成dp-mpc的换挡策略，给出改进后的完整代码
扫一扫
706
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
