# Model predictive eco–driving control for heavy–duty trucks using Branch and
Bound optimization ††thanks: This research was partially funded by the
European Union’s Horizon 2020 Research and Innovation Programme, Grant
Agreement no. 874972, Project LONGRUN.

Bart Wingelaar1, Gustavo R. Gonçalves da Silva1 and Mircea Lazar1 1All authors
were affiliated with the department of Electrical Engineering, Eindhoven
University of Technology, The Netherlands, at the time this research was
conducted<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Eco–driving (ED) can be used for fuel savings in existing vehicles, requiring
only a few hardware modifications. For this technology to be successful in a
dynamic environment, ED requires an online real–time implementable policy. In
this work, a dedicated Branch and Bound (BnB) model predictive control (MPC)
algorithm is proposed to solve the optimization part of an ED optimal control
problem. The developed MPC solution for ED is based on the following
ingredients. As a prediction model, the velocity dynamics as a function of
distance is modeled by a finite number of driving modes and gear positions.
Then we formulate an optimization problem that minimizes a cost function with
two terms: one penalizing the fuel consumption and one penalizing the trip
duration. We exploit contextual elements and use a warm–started solution to
make the BnB solver run in real–time. The results are evaluated in numerical
simulations on two routes in Israel and France and the long haul cycle of the
Vehicle Energy consumption Calculation Tool (VECTO). In comparison with a
human driver and a Pontryagin’s Minimum Principle (PMP) solution, 25.8% and
12.9% fuel savings, respectively, are achieved on average.

##  I INTRODUCTION

Heavy–duty trucks (HDTs), also known as Heavy–duty Vehicles (HDVs), are
responsible for 5% of the total emissions and a quarter of the road transport
emissions in the EU [1]. Intending to reduce these emissions for HDVs, the EU
has set goals for a 15% HDV emission reduction in 2025 and a 30% reduction in
2030 [2]. These emission reductions can be obtained as a consequence of
reduced fuel consumption, which in turn can be realized by applying
eco–driving (ED).

Eco–driving is a driving strategy and its goal is to minimize fuel consumption
by determining an optimal velocity profile, which can be defined as an optimal
control problem (OCP) [3]. Using ED is cost–effective [4] and does not require
substantial vehicle hardware modifications [5]. This makes ED a promising
fuel–saving method for existing and future HDTs with internal combustion
engines (ICEs).

Since HDTs have discrete gear transmissions, the velocity profile optimization
is generally a complex mixed–integer nonlinear programming problem. Several
approaches are available to solve the ED optimization problem, including
analytical solutions to simplified problems [6], numerical application of PMP

[7, 3, 8, 9, 10], dynamic programming (DP) [11] and some bi–level optimization
methods [12, 13]. The integrality constraint of the gear position can also be
relaxed, reducing the problem to a nonlinear programming problem as done in
[14, 15]. This ED optimization can also run in a model predictive control
(MPC) fashion (e.g., [16]). In this context, MPC can be particularly useful in
cases where spatial and time–varying elements are included in the optimization
problem, such as look–ahead information and vehicle–to–everything (V2X)
communication [17]. However, online MPC inherently requires real–time
optimization, which can be challenging [18].

In recent studies [3, 7, 8, 15], this computational burden was reduced by
optimization over driving modes. This novel method converts the
difficult–to–solve mixed–integer nonlinear programming problem with continuous
decision variables (e.g., engine and brake torques) and discrete decision
variables (e.g., gears), into an integer programming problem with only
discrete decision variables (driving modes). Thus, this formulation converts
the optimization problem into a combinatorial optimization problem. Studies
based on this approach, [3, 7, 8, 15], apply a numerical PMP solution to
achieve real-time implementation. Although a PMP algorithm is computationally
efficient, its use can lead to difficulties with constraints handling and
co–state initialization.

Therefore, to achieve a MPC solution for ED that is implementable in real-
time, it is more of interest to use solvers that are specifically aimed at
solving combinatorial problems, such as Branch and Bound (BnB). However, using
a general–purpose BnB reintroduces the high computational complexity issue,
compromising the real–time constraint. Hence, in this work we investigate and
develop a dedicated BnB MPC solver for the driving modes eco–driving problem
including slope and V2X preview information, such that real–time
implementability constraints are met. In relation with previous works [3, 15],
due to the inclusion of slope preview, we extend the velocity profile
generating capabilities with two extra driving modes: one to accelerate the
HDT and one for prolonged negative road gradients. Secondly, we introduce
three elements to make the solver real–time implementable: 1) application of
context–based (such as cost function, dynamics, constraints, and a look-ahead
horizon) node elimination—this element is not present in general–purpose BnB
solvers such as BONMIN [19]; 2) the BnB solver is warm–started with a PMP
solution [3], providing a locally optimal solution; 3) in contrast to general
BnB approaches, the lower bound is not formulated with a relaxation of the
integrality constraint, instead a custom heuristic based on the same context
information as 1) is used.

Furthermore, we include a realistic model for gearbox transmission efficiency
to improve fuel consumption assessment. Finally, we evaluate the fuel savings
in numerical simulations on three routes: one in Israel, one in France, and
the VECTO long haul cycle, and we show that the BnB solver outperforms a
human–driver model and a PMP–MPC-based solution.

The remainder of this paper is organized as follows: In section II we present
the longitudinal HDT dynamics. Then in Section III, the ED OCP is formulated
for the MPC algorithm, which can be solved with a BnB algorithm (Section IV).
In Section V, we present the implementation details and results are presented
in Section VI. Finally, conclusions are presented in Section VII.

##  II Dynamics modeling of heavy–duty trucks

In this section we introduce the main modeling components of the HDT and the
driving modes.

###  II-A Longitudinal dynamics

In this work, the longitudinal dynamics is given by modeling the HDT with a
positive velocity and considering the following powertrain components: the
rear wheels, final drive, discrete gearbox, transmission efficiencies, ICE,
(rotational) inertia and a continuous service brake. Moreover, the modeling
fundamentals are coherent with the modeling as used in [6, 15, 7, 20, 21].

Since in our problem we are also interested in including spacial elements such
as velocity limits, gradient information, and traffic signals and signs, let
us first introduce the transformation from time–domain to space–domain
$s\in\mathbb{S}$:

$\frac{\text{d}v}{\text{d}s}=\frac{\text{d}v}{\text{d}t}\frac{\text{d}t}{\text{d}s}=\frac{\sum
F}{m}\frac{1}{v},$ (1)

where $m$ is the HDT mass, state $v$ the velocity, and $\sum F$ the net force
acting on the HDT. The equivalent rotational inertia element is included in
$m$ since it is independent of the sign of $\frac{\text{d}v}{\text{d}s}$ [22].

Before defining the driving modes, we define three common elements of the
system dynamics: resistance forces, engine speed, and transmission
efficiencies.

####  II-A1 Resistance forces

The resistance forces $F_{res}$ are always present and are defined as follows:

$F_{res}=\tfrac{1}{2}\rho_{a}c_{d}A_{f}v^{2}+mgc_{rr}\,\text{cos}(\alpha)+mg\,\text{sin}(\alpha).$
(2)

In (2), $\rho_{a}$ is the air density, $c_{d}$ the aerodynamic drag
coefficient, $A_{f}$ the frontal area, $g$ the gravitational constant, and
$c_{rr}$ the roll resistance coefficient and $\alpha$ the slope.

####  II-A2 Transmission efficiencies

The efficiency of the discrete transmission and final drive ratio,
$\eta_{t}$[-] and $\eta_{r}$[-] respectively, are described by

$\displaystyle\eta_{t}$ $\displaystyle=1-\frac{T_{loss}}{|T|}$ (3)
$\displaystyle\eta_{r}$ $\displaystyle=c_{0},$ (4)

where $T$ is any torque transferred by the transmission and $c_{0}$ a
constant. $T_{loss}$ is a function of gear ratio, torque and engine speed, and
is approximated by a polynomial which is fitted on the data from [23]:

$T_{loss}=c_{1}(i_{t})\omega_{e}+c_{2}(i_{t})T+c_{3}(i_{t}).$ (5)

In (5), $c_{1}$–$c_{3}$ are gear ratio $i_{t}$–dependent constants. The used
torque loss maps are illustrated for gears 1, 8, and 12 in Fig. 1. Note that
gear 12 is in a direct drive configuration.

Figure 1: Torque loss maps as function of any transferred torque and engine
speed for gears 1,8 and 12.

####  II-A3 Engine speed

Engine speed $\omega_{e}$ [RPM] is a function of the gear position $y$ [-] and
velocity. In neutral gear ($y\\!=\\!0$), the ICE is disengaged from the
drivetrain and idles at a low engine speed. $\omega_{e}$ is defined as:

$\omega_{e}=\begin{cases}\omega_{idle}&\text{if }y\\!=\\!0,\\\
\frac{30i_{t}(y)i_{r}}{\pi r_{w}}\,v&\text{otherwise},\\\ \end{cases}$ (6)

where $r_{w}$ is the radius of the rear wheels, $i_{r}$ the final drive ratio
and $\omega_{idle}$ the idling speed.

###  II-B Driving modes

Next, we assume that the velocity profile can be followed by switching between
a finite number of subsystems, called driving modes. The resulting switched
hybrid system is defined as:

$\frac{\text{d}v}{\text{d}s}=f(M(s),v(s)),$ (7)

with $f(\cdot)$ the system dynamics of driving mode $M\in\tilde{M}$ with
$|\tilde{M}|=6$. The system dynamics is now only dependent on the driving mode
as control variable.

####  II-B1 Cruising mode $M^{cr}$

In cruising mode, the truck drives with a constant velocity and $y>0$. The ICE
delivers an engine output torque (as measured by a dynamometer) $T_{e}$, to
overcome $F_{res}>0$ and transmission losses. The dynamics in this mode are
described by:

$\frac{\text{d}v}{\text{d}s}=\frac{1}{mv}\left(\frac{\eta_{r}\eta_{t}i_{r}i_{t}(y)}{r_{w}}T_{e}-F_{res}\right)=0.$
(8)

####  II-B2 Eco–roll mode $M^{ec}$

In eco–roll mode, the ICE is disengaged from the drivetrain with the gear
position in neutral ($y\\!=\\!0$) and it rotates at an idling speed
$\omega_{idle}$. The HDT dynamics are now fully dependent on $F_{res}$ and the
equivalent rotational inertia of the drivetrain $J_{dt}$. The dynamics in this
mode are described by:

$\frac{\text{d}v}{\text{d}s}=-\frac{r_{w}^{2}}{(mr_{w}^{2}+J_{dt})v}F_{res}.$
(9)

####  II-B3 Coasting mode $M^{co}$

In coasting mode, the resistance force, equivalent powertrain rotational
inertia, and internal engine friction $T_{in.fr}$, determine the HDT dynamics.
The dynamics are described by:

$\frac{\text{d}v}{\text{d}s}=-\frac{r_{w}^{2}}{(mr_{w}^{2}+J_{pt}(y))v}\\!\left(\frac{\eta_{r}\eta_{t}i_{r}i_{t}(y)}{r_{w}}T_{in.fr}+F_{res}\right).$
(10)

####  II-B4 Engine brake mode $M^{eb}$

Besides the internal engine friction, an ICE can also deliver a braking force
with an engine brake (a type of continuous service brake). In this mode, we
assume that maximum engine brake torque $T_{eb,m}$ is applied to decelerate
the vehicle. The dynamics are described by:

$\frac{\text{d}v}{\text{d}s}\\!=\\!\frac{-1}{(mr_{w}^{2}\\!+\\!\frac{J_{pt}(y)}{r_{w}^{2}})v}\\!\left(\\!\frac{\eta_{r}\eta_{t}i_{r}i_{t}(y)}{r_{w}}\\!(T_{in.fr}\\!+\\!T_{eb,m})\\!+\\!F_{res}\\!\right)\\!,$
(11)

where $T_{eb,m}$ is the maximum braking torque of the continuous service
brake. The modeling of this brake is based on an approximation of the data in
[24].

####  II-B5 Acceleration mode $M^{ac}$

In the acceleration mode, the engine torque is used to accelerate the truck
instead of just cruising. The dynamics are described by:

$\frac{\text{d}v}{\text{d}s}=\frac{r_{w}^{2}}{(mr_{w}^{2}\\!+\\!J_{pt}(y))v}\left(\frac{\eta_{r}\eta_{t}i_{r}i_{t}(y)}{r_{w}}T_{e,ef}-F_{res}\right),$
(12)

where $T_{e,ef}$ is a fixed engine output torque profile based on the optimal
brake specific fuel consumption (BSFC) line [22]. For a given engine output
power $P_{e}=\frac{30}{\pi}T_{e}\omega_{e}$, the corresponding point,
($T_{e,ef},\omega_{e}$), on the optimal BSFC line minimizes the following:

$T_{e,ef}(\omega_{e})=\min_{T_{e}}\frac{\dot{m}_{f}(T_{e},\frac{\pi
P_{e}}{30T_{e}})}{P_{e}}.$ (13)

This line is illustrated in Fig. 2.

####  II-B6 Downhill mode $M^{dh}$

In downhill mode, the HDT maintains a constant speed in negative gradient
scenarios. Unlike cruising mode, $F_{res}<0$ due to the negative gradient.
Moreover, this mode is only enabled when any gear position in coasting mode
cannot provide sufficient braking torque to maintain a constant speed or
decelerate the HDT. Therefore, the engine brake provides just enough
additional braking torque to maintain a constant speed. The engine brake
torque is computed such that:

$\frac{\text{d}v}{\text{d}s}\\!=\\!\frac{1}{mv}\\!\left(\\!\frac{\eta_{r}\eta_{t}i_{r}i_{t}(y)}{r_{w}}(T_{in.fr}\\!+\\!T_{eb})\\!+\\!F_{res}\\!\right)=0.$
(14)

###  II-C Fuel consumption

Each driving mode has a specific fuel consumption
$\frac{\text{d}m_{f}}{\text{d}s}$, which is presented in (15):

$\frac{\text{d}m_{f}}{\text{d}s}=\begin{cases}\frac{1}{v}\dot{m}_{f}(T_{e},\omega_{e})&M\in\\{M^{cr},M^{ac}\\}\\\
\frac{1}{v}\dot{m}_{f,idle}&M\in\\{M^{ec}\\}\\\
0&M\in\\{M^{co},M^{eb},M^{dh}\\},\\\ \end{cases}$ (15)

where $\dot{m}_{f,idle}$ is a small amount of fuel consumed by the idling ICE.
The fuel consumption $\dot{m}_{f}$ is a polynomial fit which maps $T_{e}$ and
$\omega_{e}$ to the diesel fuel map, given by:

$\dot{m}_{f}=c_{4}+c_{5}\omega_{e}+c_{6}T_{e}+c_{7}\omega_{e}^{2}+c_{8}\omega_{e}T_{e}+c_{9}T_{e}^{2},$
(16)

where $c_{4}$–$c_{9}$ are constants. In Fig. 2, $\dot{m}_{f}$ is illustrated,
as well as the optimal BSFC line, maximum engine output torque constraint
$T_{e,max}$, internal engine friction $T_{in.fr}$, and engine speed
constraints.

Figure 2: 2D polynomial fit of $\dot{m}_{f}$($T_{e}$,$\omega_{e}$) [g/s] (data
taken from [15]); Optimal BSFC line, maximum $T_{e}$-constraint and
$T_{in.fr}$($\omega_{e}$)

##  III Formulation of the Eco–driving Optimal Control Problem

The ED MPC problem formulation can be described by the following OCP:

$\min_{M\in\tilde{M}}~{}J$ (17)
$J=\int_{s_{0}}^{s_{f}}\left(w_{f}\frac{\text{d}m_{f}(M,v)}{\text{d}s}+w_{t}\frac{1}{v}\right)\text{d}s+\beta(v(s_{f})-v_{f})^{2},$
(18)

where $J$ [-] is the cost to minimize, $s_{0}$ the current position, $s_{f}$
the end position, $w_{f}\\!\in\\!(0,1]$ the weight on fuel consumption,
$w_{t}\\!\in\\![0,1)$ the weight on trip duration, $\beta$ a terminal penalty
and $v_{f}$ an end velocity reference. For appropriate weight scaling, we
assume $w_{f}+w_{t}=1$. $w_{f}$ and $w_{f}$ can then be adjusted with one
parameter, $\phi$ [g/s]:

$\phi=\frac{w_{t}}{w_{f}}.$ (19)

Next, we discretize (18) in two steps. First, the distance domain $s$ is
discretized into $N$ equidistant samples with length $\delta
s=\frac{s_{f}-s_{0}}{N}$ on which the driving mode–gear combination is
constant. Secondly, the state dynamics $f(\cdot)$ are discretized with the
forward Euler method. This results in the following discrete OCP at time $k$:

$\displaystyle
J=\min_{M_{i|k}\in\tilde{M}}\sum_{i=0}^{N-1}L(\bar{v},M_{i|k})\delta
s+\beta(v_{N|k}-v_{f,N})^{2}$ (20a)
$\displaystyle\bar{v}=\frac{1}{2}(v_{i|k}+v_{i+1|k})$ (20b) $\displaystyle
L(\bar{v},M_{i|k})=w_{f}\frac{\text{d}m_{f}(\\!M_{i|k},\bar{v})}{\text{d}s}+w_{t}\frac{1}{\bar{v}}$
(20c) $\displaystyle\text{s.t.}\quad v_{0|k}=v_{0}$ (20d)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}v_{i+1|k}=v_{i|k}+f(M_{i|k},v_{i|k})\delta
s$ (20e) $\displaystyle~{}~{}~{}~{}~{}~{}~{}v_{min}(i)\leq v_{i|k}\leq
v_{max}(i)$ (20f)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\omega_{e,min}(T_{e})\leq\omega_{e}(y,v_{i|k})\leq\omega_{e,max}(T_{e})$
(20g) $\displaystyle~{}~{}~{}~{}~{}~{}~{}T_{e}\leq T_{e,max}(\omega_{e})$
(20h) $\displaystyle~{}~{}~{}~{}~{}~{}~{}T_{eb}\leq T_{eb,m}(\omega_{e}),$
(20i)

where $v_{0}$ is the current velocity of the HDT at time $k$ and (20e) is the
discretized dynamics. (20f)–(20i) represent, respectively, constraints on
minimum and maximum vehicle speed, engine speed, engine output torque and
engine brake torque.

In contrast to conventional ED OCPs, which optimize over discrete and
continuous (input) variables ($T$, $T_{eb}$ in (20)), the continuous inputs
are now fixed by the driving mode dynamics. Consequently, this new formulation
is converted into a combinatorial optimization problem, which is typically
solved using BnB methods.

##  IV Eco–driving MPC Branch and Bound algorithm

In this section we introduce the general idea of using MPC with Branch and
Bound and the main elements that compose the solution. In Section V we then
formalize the local elements that make it a dedicated solution.

###  IV-A The MPC control part

The MPC part of this algorithm consists of an outer loop where the
constraints, cost function and receding horizon are defined. An initial
warm–started solution via PMP is also obtained. The BnB algorithm is then used
to solve the constrained optimization problem. A schematic overview is
provided in Fig. 3.

Figure 3: ED MPC algorithm architecture flow chart diagram.

We assume that live route data is provided to the MPC control loop in the form
of $\alpha_{r}$, $v_{min,r}$ and $v_{max,r}$. Although temporal events
influence this data, for now it is assumed to be constant. Inside the control
loop, the incoming data is first pre-processed to extract the relevant data
for the MPC OCP ($\alpha$,$v_{min}$,$v_{max}$). Finally, a pseudo–code of the
ordered implementation steps of the ED MPC part is provided in Algorithm 1.
Notice that the algorithm already includes some context–based steps.

Inputs: $v_{0|k}$, $v_{min,r}$, $v_{max,r}$, $\alpha$, $\phi$.

1\.
$[\mathcal{X}_{PMP},{U\\!B}_{0}]=\text{WarmStart}(v_{0|k},v_{min},v_{max},\alpha,\phi)$
// run PMP algorithm and return upper bound and set of nodes
$\mathcal{X}_{PMP}$ corresponding to local optimal PMP mode sequence

2\. $[v_{min},v_{max}]=\text{NodeElim1}(v_{min,r},v_{max,r})$ // update for
feasible areas

3\. Determine $h$ // determine heuristic term in LB on a regular sized grid
for all stage–velocity combinations

4\.
$\mathcal{X}_{BnB}=\text{EDBnB}(\mathcal{X}_{PMP},{U\\!B}_{0},v_{0|k},v_{min},v_{max},\alpha,\phi,h)$
// run ED BnB and return set of visited nodes

5\. $v^{*}=\text{OptSeq}(\mathcal{X}_{BnB})$ // determine optimal Node path
from visited nodes

6\. $[M^{*},y^{*}]=\text{DetOptPath}(v^{*})$ // deduce optimal mode–gear
sequence

7\. Apply $M[0],y[0]$ // apply first input of optimal sequence

Algorithm 1 ED MPC algorithm

###  IV-B The BnB algorithm

Branch and Bound is an algorithm for solving discrete combinatorial
optimization problems, it is generally considered computationally intensive
[25], and it can converge to the global optimal solution [26]. BnB consists of
several fundamental elements: branching, bounding, a search strategy and
elimination. The branching term refers to the exploration of nodes in a tree
and the search strategy is the order in which nodes are branched (explored, or
visited). This branching operation partitions the problem into smaller sub-
problems (nodes) [27]. The search strategy can be categorized in informed
search (e.g., best–first search and heuristics) or uninformed search (e.g.,
depth–first search, breadth–first search) [28].

The bounding term refers to exploring the lower and upper bounds of the nodes.
Upper bounds (of which only the smallest is of interest in a minimization
problem) indicate a maximum of the domain on which the optimal solution of $J$
can be found. Likewise, the lower bound is a minimum on this domain. In a
typical optimization problem, each node has a lower bound, and the best
solution thus far is represented by an upper bound. The inclusion of these
bounds makes BnB inherently an informed search strategy by itself. As the
parent node is branched into child nodes in an iterative way, the lower bounds
and upper bounds are evaluated (bounding) for each child node. All child nodes
with a lower bound lower than the upper bound have the potential to yield a
better solution (a better upper bound). If the lower bound of a node is higher
than the upper bound, then branching this node further has no potential to
yield a better upper bound. Thus the node is not further explored (elimination
or pruning).

In the ED context, the BnB elements have the following meaning. A node
corresponds to the HDT velocity at a certain distance stage (or depth level)
after having applied a driving mode–gear combination for a certain distance (a
multiple of $\delta s$). The stage costs of reaching a node represent the sum
of fuel consumption and time duration. Nodes are eliminated when constraints
(for instance maximum speed) are exceeded or when the final stage is reached
with a lower total cost.

Inputs: $\mathcal{X}_{PMP},{U\\!B}_{0},v_{0|k},v_{min},v_{max},\alpha,\phi,h$

1\. Set upper bound $U\\!B=U\\!B_{0}$

2\. Set $\mathcal{X}_{BnB}=\mathcal{X}_{PMP}$

3\. Set $\mathcal{B}=v_{0|k}$ and $k=0$

4\. Explore nodes until queue is empty or until time limit is exceed:

while $k<N$ AND time$<$time_limit do

4.1 Set $\mathcal{C}=\\{\\}$

4.2 for ( $node$ in $\mathcal{B}$ ) $\\{$

$[\mathcal{C},U\\!B]=\text{singleLevelBBSearch}(\mathcal{C},node,h,k,N,U\\!B)$
// Algorithm 3

$\\}$

4.3 Sort ($\mathcal{C}$) // from lowest LB to highest LB

4.4 $U\\!B=\text{GreedyProbing}(U\\!B,\mathcal{C}[0],k,N)$ // Algorithm 4

4.5 Remove nodes in $\mathcal{C}$ with $L\\!B\geq U\\!B$

4.6
$\mathcal{C}=\text{ContextNodeElim2}(\mathcal{C},v_{min,k},v_{max,k},\varepsilon)$
// Algorithm 5

4.7 $\mathcal{B}=\mathcal{C}$

4.8 Add $\mathcal{C}$ to $\mathcal{X}_{BnB}$

4.9 $k=k+1$

end while

Output: $\mathcal{X}_{BnB}$

Algorithm 2 ED Branch and Bound algorithm

A better BnB performance can be achieved by initializing an initial upper
bound that is close to the optimal solution. This initial upper bound can be
provided by a warm-started feasible solution, such as PMP. Having a feasible
solution also means that the BnB algorithm can be terminated at any time.

Likewise, better lower bounds can also yield better BnB performance and even
speed–up computations. These lower bounds are usually composed by a stage cost
at stage $j$ and some heuristic term for the cost to go from stage $j$ to the
final stage $N$. Its main implication is that it limits the generation of
child nodes. In Section V-C we formalize the computation of this bound.

In Algorithm 2, the general idea of the BnB optimization implementation is
presented.

##  V Implementation of the BnB MPC algorithm for eco–driving

In this section we describe the elements in MPC and BnB pseudo–codes that make
this a dedicated context–based solver. First, we describe node elimination
methods based on exploiting the nature of the problem (this is used in Step 2
of MPC and Step 4.6 of BnB). We then introduce the idea of warm–starting with
a fast PMP solver and further we develop an heuristic for the lower–bound
computation. Finally we present the two search strategies used with their
respective modifications and a sorting algorithm for the lower–bound.

###  V-A Warm start with PMP

A BnB algorithm can be warm–started with a feasible initial solution: a
feasible sequence of driving modes and the cost thereof as initial upper bound
$U\\!B_{0}=J_{PMP}$. In this work, the fast PMP algorithm as developed in [3]
provides that initial solution in case a feasible solution is found
($v_{min}\leq v_{PMP}\leq v_{max}$). If such a solution is found, then the BnB
algorithm can be stopped at any time and the final result is feasible.

However, the PMP algorithm is not guaranteed to yield a feasible solution due
to constraint handling difficulties. If no feasible PMP solution is found,
$U\\!B_{0}=\infty$. To minimize constraint handling issues, the implementation
of the used PMP algorithm is described as follows. First, $v_{max}$ is
subdivided into sections where $v_{max}$ is constant. Then the final stage and
velocity of each of these sections form a velocity event. Then, the velocity
events are solved sequentially with the PMP algorithm. Finally, a terminal
penalty is added to the cumulative cost of the individual simulations to
formulate the initial upper bound
($J_{PMP}=J_{PMP}+\beta(v_{PMP,N}-v_{f,N})^{2}$ since the PMP solution solves
a two–point boundary value problem without the terminal penalty.

###  V-B Update maximum and minimum velocities

With the BnB search tree exploring the distance–velocity domain, it is of
interest to exclude infeasible node regions. Avoiding these regions leads to a
faster solver. To be more specific, real–world $v_{min}$ and $v_{max}$ are
discontinuous and create such regions. The maximum engine power and
deceleration power limit the HDT dynamics, which thereby can violate $v_{min}$
and $v_{max}$ constraints if such discontinuities occur. Although the
predictive element takes care of this problem and ensures finding a feasible
sequence of driving modes, regions of nodes whose child nodes violate these
constraints are still explored.

To overcome this problem, $v_{min}$ and $v_{max}$ were updated such that their
derivatives match the maximum acceleration and minimum deceleration of the
driving modes. This means that at the beginning of an acceleration event (for
example, after a red light), $v_{max}$ is updated forwards as using the
acceleration mode, so no node above this line is explored. For a deceleration
(so, when reaching a red light), then $v_{max}$ is updated backwards as
following engine brake mode. Similarly for $v_{min}$.

###  V-C Lower bound and admissible heuristic

A lower bound is a minimum on the domain on which the optimal solution of $J$
can be found. There are usually no formal closed–form definition of the lower
bounds (as they would imply knowledge of the optimal solution), so they are
usually composed by a sum of stage costs and some heuristic term. We formalize
this definition next. Let us first define the total cost $g$ to reach a node
at stage $j$ form stage 0:

$g=\sum_{k=0}^{j-1}L_{k}$ (21)

where $L_{k}$ is the stage cost given in (20c).

Next, the lower bound $L\\!B$ is formulated as:

$LB=g+h,$ (22)

where $h$ is a heuristic term that is an estimate on the cost from node $j$ to
the final stage $N$. Typically, this heuristic term is formulated by
relaxation of the integrality constraint [27]—which turns (20) into a
nonlinear programming (NLP) problem—, but it can also be simply 0.
Furthermore, $h$ must be admissible to retrieve the global optimal solution of
$J$.

In [15], a single case of such an NLP was solved, but this process is not
real–time implementable if applied to individual nodes. Moreover, if the
resulting NLP solution yields a local optimum, the heuristic is not
necessarily admissible and the globally optimal solution is no longer
guaranteed to be retrieved. To find a real–time implementable and admissible
heuristic, we again use context–specific elements by applying an energy and
time–based heuristic similar to [29, 30]. The added novelty is that for each
lower bound, a 2-parameter optimization problem is solved.

Consider node $i$ at stage $j$ and $\bar{v}_{h}$ as the average velocity from
stage $j$ to stage $N$. The minimum work done by the air drag can then be
described by:

$W_{air}(\bar{v}_{h})=(N-j)\delta
s\frac{1}{2}\rho_{a}c_{d}A_{f}\bar{v}_{h}^{2}.$ (23)

This minimum work does not require the velocity profile to be known. The work
of the other resistance forces $W_{d}$ from stage $j$ to stage $N$ is
described by:

$W_{d}=\sum_{k=j}^{N-1}mg(c_{rr}\text{cos}(\alpha_{k})+\text{sin}(\alpha_{k})),$
(24)

which can be determined when slope preview is given. Next, we consider the
kinetic energy (excluding rotating parts) of the HDT at the node $i$ with
velocity $v_{i}$ as:

$E_{kin,i}=\frac{1}{2}mv_{i}^{2}.$ (25)

By determining the energy difference $\Delta E_{kin}$ for going from $j$ to
$N$, the energy consumption can be estimated as:

$\Delta E_{kin}=E_{kin,N}-E_{kin,i}+W_{d}+W_{air}(\bar{v}_{h}).$ (26)

If $\Delta E_{kin}>0$, then energy is added to the HDT by means of burning
fuel. If not, energy is subtracted due to braking forces. With both the energy
consumption and time duration defined, we now formulate the total cost–to–go
$J_{iN}$ from node $i$ to the final stage:

$J_{iN}=\frac{w_{f}\text{max}(0,\Delta
E_{kin})}{Q\eta_{opt}}+\frac{w_{t}\delta
s(N\\!-\\!j)}{\bar{v}_{h}}\\!+\\!\beta(v_{N}\\!-\\!v_{f})^{2},$ (27)

where $Q$ [J/g] is the specific energy of diesel fuel and $\eta_{opt}$ an
estimated maximum combined efficiency of all conversion efficiencies from fuel
to work. Notice that the actual value of $J_{iN}$ is not known beforehand, as
it depends on the two free variables $\bar{v}_{h}$ and $v_{N}$. However, for
each node, we can solve a small optimization problem with these two variables
(all other variables are fixed) in order to provide a minimal admissible
value, which is defined as the heuristic $h$:

$h(j,v_{i})=\min_{\bar{v}_{h},v_{N}}J_{iN}(\bar{v}_{h},v_{N},j,v_{i}).$ (28)

This simple two–variable optimization is solved easily in a (partially)
analytic way and is provided to the BnB solver in a 2D lookup table (LUT) with
dimensions: node stage and node velocity on a regular sized grid. For each
point in the LUT, (28) is calculated and subsequently $h$ can be sampled from
the LUT during the BnB optimization.

Both $\bar{v}_{h}$, $v_{N}$ are constraint by $v_{max}$ and $v_{min}$.
Moreover, the search domain of $\bar{v}_{h}$ can be further reduced since
acceleration to $v_{max}$ or deceleration to $v_{min}$ is not instantaneous.
Therefore, we update the $v_{max}$ and $v_{min}$ profiles by constraining
their derivatives with maximum acceleration and deceleration. To do this
efficiently, we assume a maximum ideal acceleration which is defined as
applying maximum engine power $P_{e,max}$ to the HDT, defined as:

$P_{e,max}=\max\frac{\pi}{30}T_{e}\omega_{e},$ (29)

subjected to constraints (20g) and (20h). Furthermore, the transmission
efficiency of the gearbox and equivalent rotational inertia are neglected. The
dynamics are described by:

$\left(\frac{dv}{ds}\right)_{max}=\frac{\eta_{t}P_{e,max}-F_{res}v}{mv^{2}}.$
(30)

Similarly, the maximum braking power $P_{b,max}$ is defined as:

$P_{b,max}=\max\frac{\pi\omega_{e}}{30}(T_{in.fr}+T_{eb,m}).$ (31)

Similarly to (30), the minimum ideal deceleration is:

$\left(\frac{dv}{ds}\right)_{min}=\frac{-(\eta_{t}P_{b,max}+F_{res}v)}{mv^{2}}.$
(32)

###  V-D Search strategy

Several options exist for the search strategy. For instance, breadth–first
search explores all the nodes of a depth level (distance stage) first, before
advancing to a deeper level. Consequently, the upper bound is only updated at
the final depth level of the search tree, resulting in many visited nodes
which is memory intensive. On the other hand, depth–first search updates the
upper bound in the earliest search iterations. However, depth–first search can
explore subtrees that do not hold an optimal solution. In this regard,
informed search strategies could be a better alternative, although informed
search strategies could still end up in an exhaustive search with a large
number of visited nodes (e.g., when many nodes have a similar lower bound).

Inputs: $\mathcal{C},node,h,k,N,U\\!B$

1\. Branch $node$ in $ChildNodes$ // applying driving mode–gear combinations

2\. for ( $childNode$ in $ChildNodes$ ) $\\{$

$L\\!B=\text{DetermineLB}(childNode,h,k)$ //determine $g$ and $L\\!B$ using
(25) and (26).

2.1 if $L\\!B<U\\!B$ AND $childNode$ is feasible then

Add $childNode$ to $\mathcal{C}$

if $k==N-1$ then

$U\\!B=L\\!B$

end if

end if

$\\}$

Output: $\mathcal{C},U\\!B$

Algorithm 3 Single level bounded breadth–first search

Therefore it is of interest to utilize the strengths of both informed and
uninformed search strategies. To this extend, we formulate the search strategy
for the ED BnB algorithm as follows. We apply a custom breadth–first search
strategy (Algorithm 3) which, after each branched depth level, probes the best
child node (node with the lowest lower bound) with a greedy depth policy
(Algorithm 4). This combination has shown to perform well with context–based
node elimination (Section V-E).

The late upper bound update that inherently comes with pure breadth-first
search means that few nodes are eliminated based on their lower bounds. This
problem is magnified in case of a too low $h$ (see (22)). To overcome this, an
additional depth probing strategy in combination with a greedy policy was
used. A probing strategy keeps a certain branching policy constant and
explores the consequences thereof [31]. Only the best node (node with lowest
$L\\!B$) is depth probed at each stage after breadth-first branching. This
probing is then done with a greedy policy, where only the child node with the
lowest $L\\!B$ is branched further. This way, the upper bound can be
efficiently improved at the earliest stages yielding less node exploration.
Finally, it is important to mention that the heuristic term in the lower bound
plays an important role in the performance of these search strategies.

Inputs: $U\\!B,\mathcal{D},k,N$

2\. while $k<N$ do

2.1 branch $\mathcal{D}$ into $chilNodes$

2.2 $\mathcal{E}$ is $childNode$ with lowest $L\\!B$

if $\mathcal{E}==\O$ or $L\\!B$ of $\mathcal{E}>U\\!B$ then

2.3 return // no $U\\!B$ update

else

2.4 $\mathcal{D}=\mathcal{E}$

end if

2.5 $k=k+1$

end while

3\. $U\\!B$ = min($U\\!B$, $L\\!B$ of $\mathcal{E}$)

Output: $U\\!B$

Algorithm 4 Greedy probing algorithm

###  V-E ED context–based node elimination

Finally, the BnB problem formulation tends to have an exponential growth in
the amount of explored node regardless of the constraints and tight lower
bounds. Given the nature of the ED problem, this exponential growth can be
avoided. For example, consider two nodes at the same distance stage (or tree
depth level) with an almost identical velocity. The cost of reaching these
nodes from the root node can significantly deviate. However, the minimum cost
of reaching the final distance stage from those two nodes is approximately the
same. Therefore it makes no sense to branch both nodes into sub–trees as it
has negligible benefit. Instead, only the node with the lowest $g$ (see (21))
can be branched.

We applied this procedure as follows. For a given stage, the velocity range
$v_{max}-v_{min}$ is discretized into small domains of width $\varepsilon$
($\varepsilon=0.01$ was used). Only the node with the lowest $g$ in that
domain is kept while other nodes with a higher $g$ are eliminated. This
process is illustrated in Fig. 4 and presented in Algorithm 5.

Figure 4: Node elimination process. (left) nodes before context base node
elimination; (right) nodes after context based node elimination.

This node elimination is an important factor for making the ED MPC algorithm
real–time implementable. It effectively assumes that different nodes in the
search tree have the same physical meaning. Although this method is problem
specific, it can be applied to general problems where the nodes have a
contextual relation.

Inputs: $\mathcal{A},v_{min,k},v_{min,k},\varepsilon$

1\. Discretize $v$ into $n$ sections, yielding:
$v_{min,k},v_{min,k}\\!+\\!\varepsilon,v_{min,k},\dots,v_{max,k}\\!-\\!\varepsilon,v_{max,k}$

2\. Define nodes set $\mathcal{C}$

for ( $i=0;\ i<n;\ i=i+1$ ) $\\{$

3.1 $\mathcal{B}=\\{nodes\,|\,v_{min,k}+i\varepsilon<nodes\leq
v_{min,k}+(i+1)\varepsilon\\}\subseteq\mathcal{A}$

3.2 add $node$ with lowest $L\\!B$ of $\mathcal{B}$ to $\mathcal{C}$

$\\}$

Output: $\mathcal{C}$

Algorithm 5 Context based node elimination 2/2

##  VI Results and evaluation

In this section we assess the results obtained with the proposed ED MPC BnB
algorithm against a human driver behavior control baseline and an MPC PMP
solution for three different routes. The results are evaluated based on: total
fuel consumption, trip duration, and the real–time computation constraint. We
emphasize a comparison with the PMP solution since a comparison with a human
driver has already been presented in [3].

Three routes are assessed to avoid making route–specific conclusions. An
overview of these routes is presented in Table I. The Israel route is near
Neve Zohar, and the France route is from Paris to Rouen. Lastly, the third
route is the standardized VECTO long hauling cycle [32]. The $v_{max}$ profile
of the routes are shaped to include stop signs, curvatures, and traffic lights
events, while $v_{min}=10$ [km/h].

For standstill events we also set $v_{max}=10$ [km/h] for the following
reasons: 1) the dynamics discretization (node spacing) is too coarse for the
lowest velocities, yielding less feasible nodes; 2) the lowest velocities are
incompatible with the $\omega_{e,min}$ constraint; and 3) the dynamics are
infeasible at $v=0$. The results in this section are focused on the Israel
route due to space limitations.

TABLE I: Details of simulated routes Route | start coordinate | end coordinate | length[m]
---|---|---|---
Israel route | 30°58’09.6”N | 31°17’25.8”N | 45900
| 35°18’26.2”E | 35°22’59.8”E |
France route | 48°51’25.7”N | 49°26’35.8”N | 132940
| 2°20’29.0”E | 1°05’58.8”E |
VECTO cycle | – | – | 100180

###  VI-A Baseline reference of human driver

To assess the relative fuel savings compared to existing methods, a comparison
is made with the driving behavior of a (simulated) human driver. This human
driving behavior forms a baseline reference and consists of a speed controller
to track the maximum velocity as velocity reference and a previewing time in
which the HDT can decelerate earlier in case the velocity reference decreases
in the previewed horizon. For the longitudinal speed controller we use a PI
controller with parameters $K_{p}=10000$ and $K_{i}=1$. The preview time is
the equivalent of the anticipating behavior of a real human driver and, based
on [33], is set to

$t_{prev}=2.8+0.25(v_{0}-\text{min}(v_{m,prev})),$

with $v_{m,prev}$ the maximum velocity profile of the previewed horizon. Next,
we model the gear shift heuristic as in [34] with a gear upshift if
$\omega_{e}>2000$ and a downshift if $\omega_{e}<1000$. In addition, a low
gear is kept in downhill situations to increase the braking torque [35]. The
baseline reference velocity profile for the Israel route is illustrated in
Fig. 10. For this scenario, the total fuel consumption was 2621.5 [g] with a
trip duration of 2886 [s].

###  VI-B Simulation results: influence of the tuning parameters

In this section we analyze the influence of the tuning parameters $\phi$,
$\beta$ and $N$ on the fuel consumption and time duration of the optimized
profiles. For that, the routes are simulated for combinations of
$\phi\in\\{4,6,8,10,20,30,40,50,60,\\\ 70,80,100\\}$,
$\beta\in\\{1,5,10,20,50,100,$ $200\\}$ and $N\in\\{100,\\\ 200,300,400\\}$.

The influence of the parameter $\phi$ on fuel consumption and trip duration
illustrated in Fig. 5. One can observe that for the optimal solution an
increase in $\phi$ correlates with an increased fuel consumption and shorter
trip duration; this behavior was consistent for the three routes. A higher
$\phi$ increases the velocity and thereby the fuel consumption. This
correlation trend is also stronger for the BnB solver. Additionally, a larger
spread in data points is visible for low $\phi$ values, suggesting that the
influence of $N$ and $\beta$ increases as $\phi$ decreases. However, no trend
in the influence of $\beta$ was observed.

Figure 5: Total fuel consumption and trip duration for various simulations of
the Israel route for the ED MPC Bnb and PMP solver.

Next we examined the influence of the prediction horizon $N$, as it
significantly affects computation times (see Section VI-C). From Fig. 6 we
notice that for a larger $\phi$, the prediction horizon has a relatively small
influence on fuel consumption. Therefore, using larger $\phi$-values enables
the BnB MPC solver to also run on less powerful hardware.

Figure 6: Influence of prediction horizon $N$ on the Israel route for
different $\phi$. Observe that a larger $N$ does not always yield a lower fuel
consumption.

On the other hand, a lower $\phi$ increases the influence of $N$ on fuel
consumption, although a larger horizon does not necessarily result in lower
fuel savings and $N=200$ seems to be optimal. The main cause for this is that
the average $v_{max}$–constraint (which is approximately followed) and the
relative potential energy of the HDT are prediction horizon–dependent, such
that the fuel consumption and trip duration do not increase proportionally
with the horizon increase.

###  VI-C Real–time constraint

The algorithms were simulated in Matlab, which has no explicit code
compilation step. Once compiled in C++, for instance, a 20x speedup can be
achieved [7]. In [8] the real–time constraint is quantified with a policy
update of 1 [Hz] for a driver advisory system, resulting in a time constraint
of 20 [s] on the optimization time in Matlab. The average optimization
computation times of the routes are presented in Table II, along with the
estimated computation time in C++. One can observe that the constraint is met
and that the computation time increases disproportionately with a horizon
increase.

TABLE II: Average computation time of the optimization intervals in Matlab and C++. The C++ times are estimated and based on a 20x speedup [7]. Values must remain under 20[s] and 1[s], respectively for Matlab and C++. | Computation time [s]
---|---
| Israel route | France route | VECTO cycle
horizon $N$[-] | Matlab | C++ | Matlab | C++ | Matlab | C++
$N$=100 | 0.1729 | 0.0086 | 0.2590 | 0.0130 | 0.2399 | 0.0120
$N$=200 | 0.5202 | 0.0260 | 0.7364 | 0.0368 | 0.6598 | 0.0330
$N$=300 | 1.0387 | 0.0519 | 1.3930 | 0.0697 | 1.2316 | 0.0616
$N$=400 | 1.7609 | 0.0880 | 2.2655 | 0.1133 | 2.0182 | 0.1009
Figure 7: Israel route. Respectively from top to bottom:
velocity–trajectories; driving mode; gear position. Simulation settings:
$N=200$, $\phi=10$, $\beta=10$. Results of this particular simulation (total
fuel consumed[g], trip duration[s]): PMP ED MPC: 2901[s],9722[g]; BnB ED MPC:
2927[s],8601[g]; Performance: BnB vs. PMP: 0.90% time increase, 11.53% fuel
saved; BnB vs. baseline reference: 1.43 % time increase, 31.85% fuel saved.

###  VI-D Simulation results: final route simulation

Based on the observations made in the previous sections and given that in
practical applications, especially in delivery situations, only a certain in
increase in trip duration is admissible, we now set our choice of parameters
as follows. A maximum trip duration increase of 5% with respect to the
baseline reference trip duration is allowed, leading to a choice of $\phi=10$.
Based on the trend in Fig. 6 and a trade–off with computational time, we also
selected $N=200$. Finally, since for this range of parameters the final
parameter $\beta$ has no clear influence, we used $\beta=10$.

The total fuel consumption, trip duration, and relative fuel savings for the
three routes are presented in Table III. One can observe that the BnB solver
consistently outperformed the PMP solution and the human driver baseline
reference. In addition, it can also be observed in Table III that the highest
fuel savings occurred on the Israel route, while this is also the shortest
route (Table I). It is therefore likely that extending the Israel route will
yield similar fuel savings as the other two. The resulting velocity profile
with corresponding driving mode and gear position for the Israel route is
illustrated in Fig. 7.

One can observe that the BnB and PMP velocity profiles are similar and
relatively close to $v_{max}$. The main differences are mostly observed at
road gradient changes and $v_{max}$ transients. Consider for instance the
velocities around $11000$ [m]: around this point, the BnB velocity is lower
than the PMP velocity and the altitude profile indicates a local maximum. This
means that part of potential energy is temporarily converted into kinetic
energy and that BnB is able to capture this advantage more efficiently.

This performance difference can be explained by two main differences: while
PMP is only guaranteed to be locally optimal [36] (by the the initial co-state
guess) and focuses in matching the boundary constraints, the BnB algorithm is
approximately globally optimal (from a driving modes perspective) and
therefore finds a lower cost.

Finally, one can observe that the driving modes and gear switching occur can
rather frequently since it is not penalized in the cost function. For a more
realistic outcome, a constraint on the number of switches per distance could
be added.

TABLE III: Total fuel consumption and trip duration of the baseline reference, and relative fuel savings of ED MPC BnB and PMP. | baseline reference | fuel savings [%]
---|---|---
| Fuel con. | Trip dur. | BnB vs. | PMP vs. | BnB vs.
route [-] | [g] | [s] | b. ref. | b. ref. | PMP
Israel route | 12621.5 | 2886.0 | 31.85 | 22.97 | 11.53
France route | 42904.6 | 7545.6 | 23.18 | 9.95 | 14.87
VECTO cycle | 27898.7 | 4674.2 | 18.56 | 12.03 | 7.44

##  VII CONCLUSIONS

In this work, we proposed a dedicated MPC BnB solver for a driving modes ED
problem including slope and V2X preview information to meet real–time
implementability constraints. We introduced two new driving modes for the
inclusion of look–ahead slope information in the OCP. To make the solver
real–time implementable, context–based node elimination was applied, a
context–based heuristic in the lower bound was used and the solver was
warm–started with a PMP solution. Then, we evaluated the results of the BnB
solver, as well as a human driver model and a MPC PMP–based solution in a
numerical simulation of three routes. In comparison with a human driver and a
PMP solution, respectively, up to 25.8% and 12.9% fuel savings are achieved on
average. The BnB solver also consistently outperforms the human driver and the
MPC PMP-–based solution. These fuel savings were consistently achieved at the
lowest time/fuel weight ratios in the cost function, but not necessarily at
the largest prediction horizon.

In future work, it can be investigated whether the developed driving modes
method approached the optimal solution that is based on continuous control
variables. Moreover, three things are expected to improve performance: 1) the
heuristic can be improved using a feedforward artificial neural network; 2)
the nodes can be branched in parallel; 3) the warm–started initial solution
can be improved by also using the optimal driving mode sequence of the
previous optimization interval.

## References

  * [1] E. E. Agency, “Carbon dioxide emissions from Europe’s heavy-duty vehicles,” https://www.eea.europa.eu/themes/transport/heavy-duty-vehicles/carbon-dioxide-emissions-europe, 2018, accessed: 2021-3-30.
  * [2] E. Commission, “Reducing co2 emissions from heavy-duty vehicles,” https://ec.europa.eu/clima/policies/transport/vehicles/heavy\\_en, 2020, accessed: 2021-04-02.
  * [3] B. Wingelaar, G. R. Gonçalves da Silva, M. Lazar, Y. Chen, and J. T. B. A. Kessels, “Design and assessment of an eco–driving pmp algorithm for optimal deceleration and gear shifting in trucks,” _IEEE Conference on Control Technology and Applications (CCTA)_ , pp. 8–13, 2021.
  * [4] Y. Xu, H. Li, H. Liu, M. O. Rodgers, and R. L. Guensler, “Eco-driving for transit: An effective strategy to conserve fuel and emissions,” _Applied Energy_ , vol. 194, pp. 784–797, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0306261916314040
  * [5] D. Shen, D. Karbowski, and A. Rousseau, “A minimum principle-based algorithm for energy-efficient eco-driving of electric vehicles in various traffic and road conditions,” _IEEE Transactions on Intelligent Vehicles_ , vol. 5, no. 4, pp. 725–737, 2020.
  * [6] A. Sciarretta, G. De Nunzio, and L. Ojeda, “Optimal ecodriving control: Energy-efficient driving of road vehicles as an optimal control problem,” _IEEE Control Systems Magazine_ , vol. 35, pp. 71–90, 2015.
  * [7] N. Rozkvas, “Real-time velocity and driving mode optimization for hybrid trucks using previewed road and traffic data,” Master’s thesis, Eindhoven University of Technology, 2019.
  * [8] E. C. H. Thomassen, “Anticipative velocity optimization for heavy duty trucks using road geometry information,” Master’s thesis, Eindhoven University of Technology, 2017.
  * [9] J. Zhu, C. Ngo, and A. Sciarretta, “Real-time optimal eco-driving for hybrid-electric vehicles,” _IFAC-PapersOnLine_ , vol. 52, no. 5, pp. 562–567, 2019, 9th IFAC Symposium on Advances in Automotive Control AAC 2019\.
  * [10] J. Hu, Y. Shao, Z. Sun, M. Wang, J. Bared, and P. Huang, “Integrated optimal eco-driving on rolling terrain for hybrid electric vehicle with vehicle-infrastructure communication,” _Transportation Research Part C-emerging Technologies_ , vol. 68, pp. 228–244, 2016.
  * [11] E. Hellström, M. Ivarsson, J. Åslund, and L. Nielsen, “Look-ahead control for heavy trucks to minimize trip time and fuel consumption,” _Control Engineering Practice_ , vol. 17, pp. 245–254, 2009.
  * [12] A. Hamednia, N. K. Sharma, N. Murgovski, and J. Fredriksson, “Computationally efficient algorithm for eco-driving over long look-ahead horizons,” _IEEE Transactions on Intelligent Transportation Systems_ , pp. 1–15, 2021.
  * [13] L. Guo, B. Gao, Y. Gao, and H. Chen, “Optimal energy management for hevs in eco-driving applications using bi-level mpc,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 18, no. 8, pp. 2153–2162, 2017.
  * [14] N. Robuschi, C. Zeile, S. Sager, and F. Braghin, “Multiphase mixed-integer nonlinear optimal control of hybrid electric vehicles,” _Automatica_ , vol. 123, p. 109325, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0005109820305252
  * [15] Y. Chen and M. Lazar, “Real-time driving mode advice for eco-driving using MPC,” in _IFAC World Congress_ , Germany, 2020.
  * [16] H. Chen, L. Guo, H. Ding, Y. Li, and B. Gao, “Real-time predictive cruise control for eco-driving taking into account traffic constraints,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 20, no. 8, pp. 2858–2868, 2019.
  * [17] M. Held, “Fuel-efficient look-ahead control for heavy-duty vehicles with varying velocity demands,” Ph.D. dissertation, KTH Royal Institute of Technology, 2020.
  * [18] J. Zhang and T. Shen, “Real-Time Fuel Economy Optimization With Nonlinear MPC for PHEVs,” _IEEE Transactions on Control Systems Technology_ , vol. 24, no. 6, pp. 2167–2175, 2016.
  * [19] P. Bonami and J. Lee, “Bonmin user’s manual,” 2008.
  * [20] J. Holmberg, “Heavy duty vehicle simulation tool for the calculation of fuel consumption and development of a new vecto gearshift algorithm,” Master’s thesis, KTH Royal Institute of Technology, 2018.
  * [21] B. Saerens, “Optimal control based eco-driving,” Ph.D. dissertation, KU Leuven, 2012.
  * [22] D. Forsberg and M. Hall, “Reduced fuel consumption of heavy-duty vehicles using pulse and glide,” Master’s thesis, Linköping University, 2019.
  * [23] F. Rodriguez, “Fuel consumption simulation of HDVS in the EU: comparisons and limitations,” in _International Council on Clean Transportation (ICCT)_ , 2018.
  * [24] DAF, “Retarder systems,” https://www.daf.com/-/media/files/document-library/infosheets/general/brakes/daf-retarder-systems-71053-hqp.pdf, 2021, accessed: 2021-10-15.
  * [25] P. Somol, P. Pudil, F. Ferri, and J. Kittler, “Fast branch and bound algorithm in feature selection,” in _Proceedings of World Multiconference on Systemics, Cybernetics and Informatics_ , vol. 7, 2000, pp. 1646–1651.
  * [26] L. F. Mendonça, J. M. Sousa, and J. S. da Costa, “Optimization problems in multivariable fuzzy predictive control,” _Int. J. Approx. Reason._ , vol. 36, pp. 199–221, 2004.
  * [27] D. R. Morrison, S. H. Jacobson, J. J. Sauppe, and E. C. Sewell, “Branch-and-bound algorithms: A survey of recent advances in searching, branching, and pruning,” _Discrete Optimization_ , 2016.
  * [28] Baeldung, “The informed vs. uninformed search algorithms,” https://www.baeldung.com/cs/informed-vs-uninformed-search, 2021, accessed: 2021-11-2.
  * [29] O. Chevrant-Breton, T. Guan, and C. W. Frey, “Search space reduction in dynamic programming using monotonic heuristics in the context of model predictive optimization,” in _17th International IEEE Conference on Intelligent Transportation Systems (ITSC)_ , 2014.
  * [30] Z. Ajanovic, M. Stolz, and M. Horn, “Novel model-based heuristics for energy optimal motion planning of an autonomous vehicle using A*,” in _Proc. 15th IFAC Symposium on Control in Transportation Systems (CTS 2018)_ , 2018.
  * [31] M. W. P. Savelsbergh, “Preprocessing and probing techniques for mixed integer programming problems,” _ORSA Journal on Computing_ , vol. 6, pp. 445–454, 1994.
  * [32] E. Comission, “Vehicle energy consumption calculation tool - vecto,” https://ec.europa.eu/clima/eu-action/transport-emissions/road-transport-reducing-co2-emissions-vehicles/vehicle-energy_nl, 2021, accessed: 2021-11-2.
  * [33] D.-w. Xing, X.-s. Li, X.-l. Zheng, Y.-y. Ren, and Y. Ishiwatari, “Study on driver’s preview time based on field tests,” in _2017 4th International Conference on Transportation Information and Safety (ICTIS)_. IEEE, 2017, pp. 575–580.
  * [34] W. Enang, C. Bannister, C. Brace, and C. Vagg, “Modelling and heuristic control of a parallel hybrid electric vehicle,” _Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering_ , vol. 229, no. 11, pp. 1494–1513, 2015.
  * [35] B. Song, J. Lv, Y. Liu, and F. Kong, “The simulation and analysis on engine and hydraulic retarder continual braking performance of the tracked vehicle on long downhill,” in _2009 9th International Conference on Electronic Measurement Instruments_ , 2009, pp. 3–928–3–931.
  * [36] M. Kelly, “An introduction to trajectory optimization: How to do your own direct collocation,” _SIAM Rev._ , vol. 59, no. 4, p. 849–904, 2017.

## APPENDIX

In this section we present the results of the France route and VECTO cycle.
The results of the Israel route are presented in Section VI.

Figure 8: France route. Respectively from top to bottom: velocity-
trajectories; driving mode; gear position. Simulation settings: $N=200$,
$\phi=10$, $\beta=10$. Results of this particular simulation (total fuel
consumed[g], trip duration[s]): PMP ED MPC: 7621[s],38909[g]; BnB ED MPC:
7706[s],35976[g]; Performance: BnB vs. PMP: 1.12% time increase, 7.54% fuel
saved; BnB vs. baseline reference: 2.13% time increase, 16.15% fuel saved.
Figure 9: VECTO long haul cycle. Respectively from top to bottom: velocity-
trajectories; driving mode; gear position. Simulation settings: $N=200$,
$\phi=10$, $\beta=10$. Results of this particular simulation (total fuel
consumed[g], trip duration[s]): PMP ED MPC:4686[s],24547[g]; BnB ED MPC:
4719[s],23671[g]; Performance: BnB vs. PMP: 0.70% time increase, 3.57% fuel
saved; BnB vs. baseline reference: 0.96% time increase, 15.15% fuel saved.
Figure 10: Simulated human driver velocity profiles. Respectively from top to
bottom: Israel route, France route, VECTO long haul cycle. Figure 11: Total
fuel consumption and trip duration for various simulations of the France route
(top) and VECTO cycle (bottom). This figure complements Fig. 5.  Figure 12:
Influence of horizon $N$ on the France route (left) and VECTO cycle (right).
This figure complements Fig. 6.