HEV example

The scripts developed in this example can be opened by enetering

open('example_hev')
open('hev')

in the command window.

Consider a p2 parallel HEV powertrain.

The vehicle must drive following a speed trace in time while minimizing the fuel consumption. This approach is particularly convenient for evaluating fuel consumption and emissions of a vehicle over a regulatory drive cycle.

The objective of this example is to design a control strategy which defines the gear number to be engaged by the gearbox and the torque that the electrical machine (EM) must provide or absorb. Furthermore, we must make sure that the battery's state of charge at the end of the drive cycle is equal to its initial value.

Problem definition

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, defined as:

where

is the torque provided (if positive) or absorbed (if negative) by the electrical machine and T_req

is the torque requested at the powertrain level to make the vehicle following the given speed trace. As such, T_req

is a function of the vehicle's speed and the gear number.

Since we must set a constraint on the terminal value of the battery's state of charge (SOC), we must be able to assign an initial condition to it and track its evolution in time (in other words, throughout the decision problem's stages). Thus, the battery SOC is set as a state variable, and we must define its dynamics. The SOC dynamics is derived using a simple battery internal resistance model

I_(batt) = (V_OC - sqrt(V_OC^2 - 4 R_eq P_batt)) / (2 R_eq)

Where

and

are the battery current, open-circuit voltage, equivalent resistance and capacity. The battery power P_batt

is evaluated as

$P_EM = eta_EM omega_EM T_EM if T_EM omega_EM >= 0 ; 1\eta_EM omega_EM T_EM if T_EM omega_EM < 0.$

with

and

being the efficiencies of the electrical machine and the inverter (or any other power electronics).

Finally, the engine's fuel consumption (which we will use as our cost function) is evaluated based on the engine speed and torque with its fuel consumption map

In order for the optimization problem to be meaningful, we must set a series of constraints that reflect the operational constraints of an actual powertrain.

Constraints

Engine

The engine's torque cannot exceed its limit torque (which is in turn dependant on its speed):

T_eng = (1-tau) T_req <= T_eng,lim(omega_eng)

Electrical Machine

The electrical machine's torque must stay within its limit torque in generation and motor mode:

When braking (i.e. T_req < 0

), the electrical machine should never operate in motor mode.

Battery

The terminal SOC must be equal to its initial value:

Also, the battery is characterized by a maximum discharge power and a maximum charge current that must not be exceeded:

Define the system and cost function

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

open('sysfun_template.m')

and creating a copy in your working folder.

function [x_new, stageCost, unfeas] = sysname(x, u, ~)
    ...
end

For more information on the system and cost function signature, see Set up a basic problem.

Let's start by modifying the function signature. For this problem, we want to include an exogenous input in the system dynamics: therefore, we replace the tilde ~ in the third input with a proper variable name, such as w.

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

The system dynamics and stage cost of this problem are particularly complex and they require a lot of data to define the powertrain components, such as the engine limit torque characteristic $T_eng,lim(\omega_eng)$ and the fuel consumption map

. Theoretically, we could create variables to define these characteristics, as well as all the remaining vehcile data, in the system and cost function itself, or we could store them in a .mat file and load them when they are needed. However, since the DynaProg optimization algorithm evaluates the system and cost function a very large number of times, this is strongly not recommended.

Instead, it is best to define all of the vehicle data in a separate function, independent of DynaProg, and the pass it to the system as cost function as additional inputs. In this example, the vehicle data is stored in six structures (veh, fd, gb, eng, em and batt).

function [x_new, stageCost, unfeas] = hev(x, u, w, veh, fd, gb, eng, em, batt)
    ...
end

Finally, it might be interesting to also include in the optimization results the time profiles of physical quantities other than the state variables, control variables and cost. To do this, we change the function signature to return additional outputs starting from the fourth postional output.

function [x_new, stageCost, unfeas, engTrq, emTrq] = hev(x, u, w, veh, fd, gb, eng, em, batt)
    ...
end

Now that the function signature has been defined, we can write the equations defining the state dynamics, stage cost and unfeasibilities. To see the full code for this example, type

open('hev')
open('example_hev')

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.

%% Set up the problem
clear
% State variable grid
SVnames = 'SOC';
x_grid = {0.4:0.001:0.7};
% Initial state 
x_init = {0.6};
% Final state constraints
x_final = {[0.599 0.601]};
 
% Control variable grid
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] = hev_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, ...
 @(x, u, w) hev(x, u, w, veh, fd, gb, eng, em, batt), 'ExogenousInput', w);

Note that, since we are using two exogenous inputs, we passed a cell array containing them using the name-value pair ('ExogenousInput', w).

Also, since we are using additional inputs (structures veh, fd, gb, eng, em and batt), we have specified the system and cost function using a function handle which parametrizes them:

@(x, u, w) hev(x, u, w, veh, fd, gb, eng, em, batt)

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 two additional outputs in the system and cost function, prob also contains their time profiles in prob.AddOutputsProfile.

Retrieve the additional outputs

engTrq = prob.AddOutputsProfile{1};
emTrq = prob.AddOutputsProfile{2};

Pot results and replace the cumulative cost with the engine and e-machine torque

figure
t = plot(prob);
nexttile(t, 5)
plot(time_s, engTrq, 'LineWidth', 1.5)
hold on
plot(time_s, emTrq, 'LineWidth', 1.5)
legend('Engine torque, Nm', 'EM torque, Nm', 'FontSize', 10)