HEV with a Three-Way Catalyst example

The scripts developed in this example can be opened by entering

open('example_hev_split')
open('hev_ext')
open('hev_int')

in the command window.

Consider the Hybrid-Electric vehicle example, where the objective wad to minimize the fuel consumption of an HEV driving through a given driving mission. However, instead of minimizing the fuel consumption only, we want to minimize a trade-off of the fuel consumption and tailpipe pollutants emissions. The vehicle is equipped with a three-way catalyst (TWC), which is used to reduce engine-out emissions.

Problem definition

As for the HEV example, the vehicle speed and acceleration in time is treated as an exogenous input, while the control variables are the gear number and the EM torque ratio.

The battery's state of charge (SOC) is set as a state variable, so that we must are able to assign an initial condition and a terminal constraint.

Engine-out emissions are evaluated based on the engine speed and torque with the pollutants maps

Then, tailpipe emissions must be evaluated by characterizing the pollutants conversion in the TWC. The TWC is characterized by the reduction efficiency maps

and

, which evaluate a conversion efficiency for each pollutant as a function of the mass flow rate of the exhaust gases flowing through the catalyst and its temperature.

For this reason, we must be able to evaluate and track the evolution of the TWC temperature in time, which means that we must assign it as a state variable. Its dynamics are

where

is the catalyst's heat capacity and

is the net heat flow, which includes convective heat transfer with the exhaust gases, chemical heat release due to the pollutants reduction reactions, radiative and convective heat transfer towards the surrounding environment.

These heat flows are in turn functions of the exhaust gases flow rate and temperature, the engine-out and emissions, and the TWC temperature itself.

To sum up, the TWC temperature dynamics are ultimately dependent on the TWC temperature itself and both the control variables, as well as the exogenous inputs.

Define the system and cost function

The problem presented in example involves a relatively complex model, with two state variables, two control variables and many equations involved. Therefore, it belongs to the class of problems that may benefit (in terms of optimization time) from splitting the model function.

First, we analyze what part of the model function we can evaluate without involving the state variables. In our case, this roughly coincides with evaluating the electrical machine's electrical power demand, the engine exhaust gases, fuel and pollutants mass flow rate and the exhaust gases temperature. These will therefore be our intermediate variables.

Start with the basic function signature. You can start by opening a template m-file:

open('sysfun_ext_template.m')

and creating a copy in your working folder.

function [intVar, unfeas] = sysname_ext(u, w)
    ...
end

For more information on the system and cost function signature for split models, see Splitting the model function.

Let's modify the function signature. Similarly to the HEV example, we want to include an exogenous input in the system dynamics and we want to pass the vehicle data as additional inputs. For the external funcion, we need the vehicle data stored in the structures veh, fd, gb, eng and em.

function [intVar, unfeas] = hev_ext(u, w, veh, fd, gb, eng, em)
    ...
    intVar{1} = emElPwr;
    intVar{2} = fuelFlwRate;
    intVar{3} = egTemp;
    intVar{4} = hcFlwRate;
    intVar{5} = coFlwRate;
    intVar{6} = noxFlwRate;
end

Then, we include all operations needed to evaluate the intermediate variables and all constraints that we can set at this stage. After evaluating our intermediate variables, which we will need in the interal function, we must return them in a cell array wich is the first output of the external function.

Next, we build the internal function. This will include all the operations which are also dependant on the state variables.

Start with the basic function signature. You can start by opening a template m-file:

open('sysfun_int_template.m')

and creating a copy in your working folder.

function [x_new, stageCost, unfeas] = sysname_int(x, u, w, intVar)
    ...
end

Once again, we modify the function signature. In this case, we no longer need the control variables and exogenous inputs, so we can suppress the second and third output with a tilde ~. Moreover, we want to include the battery and TWC data a additional inputs, so we specify them starting from the fifth input argument as the fourth argument is reserved to the intermiediate variables. Finally, we return the fuel consumption and tailpipe pollutants as additional outputs to make them available when analyzing the optimization results.

function [x_new, stageCost, unfeas, fuelFlwRate, tpEmissions] = hev_int(x, ~, ~, intVar, batt, twc)
    fuelFlwRate = intVar{2};
    ...
    % Additional output
    tpEmissions.hc = hcTailpipeFlwRate;
    tpEmissions.co = coTailpipeFlwRate;
    tpEmissions.nox = noxTailpipeFlwRate;
end

To see the full code for this example, type

open('hev_int')
open('hev_ext')
open('example_hev_split')

in the command window.

Set up the problem structure

Create another script where you will define and solve the optimization problem.

Now, create all the required information that define an optimization problem. See also How to set up a basic problem and Splitting the model function.

%% Set up the problem
clear
% State variables grids
SVnames = ["SOC", "TWC Temperature"];
x1_grid = 0.4:0.001:0.7;
x2_grid = 200:10:600;
x_grid = {x1_grid, x2_grid};
% Initial state 
x1_init = 0.6;
x2_init = 300;
x_init = {x1_init, x2_init};
% Final state constraints
x1_final = [0.599 0.601];
x2_final = [];
x_final = {x1_final, x2_final};
% Control variables grids
CVnames = ["Gear Number", "Torque split"];
u1_grid = [1 2 3 4 5];
u2_grid = -1:0.1:1;
u_grid = {u1_grid, u2_grid};
 
% Load a drive cycle
load UDDS % contains velocity and time vectors
dt = time_s(2) - time_s(1);
% Create exogenous input
w{1} = speed_kmh./3.6; 
w{2} = [diff(w{1})/dt; 0]; 
 
% Number of stages (time intervals)
Nint = length(time_s);
 
% Generate and store vehicle data
[veh, fd, gb, eng, em, batt, twc] = data();

Create the problem structure object by passing all relevant information to the DynaProg function.

% Create DynaProg object
prob = DynaProg(x_grid, x_init, x_final, u_grid, Nint, ...
 @(u, w) hev_ext(u, w, veh, fd, gb, eng, em), ...
 @(x, u, w, intVar) hev_int(x, u, w, intVar, batt, twc), ...
 'ExogenousInput', w);

Note that, since we are using a split model, we specified two function handles: the first points to the external function, while the second points to the internal function.

Also, since we are using additional inputs (structures veh, fd, gb, eng, em, batt and twc), we specify the system and cost functions using function handles which parametrizes them:

@(u, w) hev_ext(u, w, veh, fd, gb, eng, em)
@(x, u, w, intVar) hev_int(x, u, w, intVar, batt, twc)

Run the optimization and visualize results

Run the problem by using run:

prob = run(prob);

Make sure you return the problem structure as an output. Once the optimization algorithm is done running, the updated problem structure prob will also contain the simulation results. You can access them by querying prob.StateProfile, prob.ControlProfile and prob.CostProfile.

Since we defined additional outputs in the system and cost function, prob also contains their time profiles in prob.AddOutputsProfile.

Retrieve the additional outputs

fc_grams = prob.AddOutputsProfile{1};
hcTailpipeFlwRate = [prob.AddOutputsProfile{2}.hc];
coTailpipeFlwRate = [prob.AddOutputsProfile{2}.co];
noxTailpipeFlwRate = [prob.AddOutputsProfile{2}.nox];

Pot results and replace the cumulative cost with the fuel consumption and tailpipe pollutants

figure
t = plot(prob);
nexttile(t, 5)
plot(prob.Time(1:end-1), fc_grams.*1e-2)
hold on
plot(prob.Time(1:end-1), hcTailpipeFlwRate)
plot(prob.Time(1:end-1), coTailpipeFlwRate)
plot(prob.Time(1:end-1), noxTailpipeFlwRate)
legend('Fuel consumption, g/s \cdot 10^{-2}', 'Tailpipe HC, g/s', 'Tailpipe CO, g/s', ...
    'Tailpipe NOx, g/s', 'FontSize', 10)