# Evolution of Fullerenes in Circumstellar Envelopes by Carbon Condensation:
Insights from Reactive Molecular Dynamics Simulations

Zhisen Meng  Laboratory for Relativistic Astrophysics, Department of Physics,
Guangxi University, 530004, Nanning, China  Zhao Wang<EMAIL_ADDRESS>Laboratory for Relativistic Astrophysics, Department of Physics, Guangxi
University, 530004, Nanning, China

###### Abstract

Fullerenes, including C60 and C70, have been detected in various astronomical
environments. Understanding how their structures evolve over time is essential
for gaining insights into their life cycle and making further observations. To
address this, we conducted reactive molecular dynamics simulations to
investigate the evolution of fullerenes in the circumstellar envelopes
surrounding carbon-rich asymptotic giant branch stars. Our simulations
employed a bottom-up chemistry scheme, wherein fullerenes grow by absorbing
and condensing small carbon-based molecules. The results revealed the
formation of different structures through heterogeneous reactions based on
hydrogen concentration, leading to the emergence of onion-like nanostructures
or single-layer fullerenes. To examine the impact of these structural changes
on the infrared emission characteristics of fullerenes, we performed quantum
chemical calculations. The results indicate that as fullerenes grow larger,
additional emission features are introduced in the infrared spectrum.
Moreover, two-layered fullerenes show noticeable blue-shift or weakening
effects on the bands associated with out-of-plane vibration modes.

##  I Introduction

Fullerenes, such as C60 and C70, were discovered by Kroto et al. (1985) in a
laser ablation experiment on graphite, and were then confirmed by
spectroscopic studies (Kratschmer et al., 1990a, b; Taylor et al., 1990).
Although they are considered to be stable in harsh interstellar conditions,
direct observational evidence of cosmic fullerenes was not obtained until 2010
(Cami et al., 2010). Fullerenes have since been reported to exist in various
astronomical environments, including asymptotic giant branch (AGB) stars
(Gielen et al., 2011), reflection nebulae (Berne et al., 2013; Sellgren et
al., 2010), planetary nebulae (Cami et al., 2010; Garcia-Hernandez et al.,
2010), young stellar objects (Roberts et al., 2012), photodissociation regions
(Castellanos et al., 2014), and diffuse interstellar medium (ISM) (Berné et
al., 2017). Fullerenes are vital to the life cycle of large carbonaceous
molecules/particles in interstellar space (Candian et al., 2018).
Understanding their formation and evolution is hence crucial not only for
their further observation but also for comprehending the evolution of carbon-
rich dust, which is intricately linked to the formation, life, and death of
stars (Otsuka et al., 2014).

There are two primary mechanisms for the formation of fullerenes: top-down and
bottom-up chemistry. Fullerenes have often been observed to coexist with
polycyclic aromatic hydrocarbons (PAHs) (Sellgren et al., 2010; Dunk et al.,
2013; García-Hernández et al., 2013), with abundance observed to decrease with
increasing concentration of PAHs (Berné and Tielens, 2012). This has led to
the exploration of the top-down formation mechanism, in which large PAH
molecules or other carbonaceous dust are transformed into fullerenes (Chuvilin
et al., 2010; Zhang et al., 2013; Pietrucci and Andreoni, 2014; Garcia-
Hernandez et al., 2010; Micelotta et al., 2012). Dehydrogenation of PAHs is
considered a critical step in this process, occurring before the folding of a
large PAH into a fullerene cage (Montillaud et al., 2013; Berné et al., 2015;
Patra et al., 2014).

On the other hand, the bottom-up approach, which is more commonly considered,
is assumed to take place in the dense and hot envelopes of evolved stars. It
typically involves reactions between carbon chain molecules followed by
hydrocarbon addition (Cherchneff et al., 2000; Pascoli and Polleux, 2000;
Bettens and Herbst, 1997; Jones et al., 2011; Parker et al., 2012). Based on
ab initio calculations (Bates and Scuseria, 1997), it has been demonstrated
that closed-shell molecules with adjacent pentagons exhibit thermodynamic
stability. Subsequently, Morokuma, Irle, and coworkers introduced the concept
of the “shrinking hot giant road” using density functional tight-binding
molecular dynamics simulations (Irle et al., 2003a, b; Zheng et al., 2004,
2005, 2007; Irle et al., 2006, 2007; Saha et al., 2009). This proposed
mechanism suggests that the formation of larger fullerene molecules
(comprising more than 100 atoms) initiates from carbon chains through self-
catalytic reactions within the thermal carbon vapor. As these larger
fullerenes are exposed to ultraviolet radiation, they undergo shrinkage
reactions, transforming into C60 and C70, while releasing smaller molecular
fragments like C2 (Irle et al., 2007; Berné and Tielens, 2012; Berné et al.,
2015).

Despite extensive efforts to elucidate the formation mechanisms of
interstellar fullerenes, their subsequent evolution in ISM remains poorly
understood. The scarcity of C60 detections in planetary nebulae suggests that
these molecules may undergo evolutionary reactions and structural
transformations in space, leading to significant changes in their spectral
features (Otsuka et al., 2014; Zhen et al., 2019). Laboratory experiments are
valuable for simulating interstellar reactions, but they often face challenges
in determining the chemical pathway due to the presence of highly reactive
intermediates with short lifetimes. Quantum chemical calculations (QCCs)
provide a solution by assisting in the identification of intermediates and the
elucidation of reaction pathways (Parker and Kaiser, 2017; Fioroni et al.,
2019). However, the computational cost associated with QCCs often limits their
applicability in studying large molecules such as fullerenes. In this context,
classical molecular dynamics (MD) simulations based on reactive force fields
offer a promising alternative (Meuwly, 2019). These simulations strike a
balance between computational accuracy and efficiency, enabling the
investigation of interstellar formation of carbon-based large molecules and
nanoparticles under complex conditions (Patra et al., 2014; Montillaud and
Joblin, 2014; Qi et al., 2018; Chen, 2018; Chen et al., 2020; Ostroumova et
al., 2019; Hanine et al., 2020; Orekhov et al., 2020; Anders and Urbassek,
2021; Tang et al., 2022).

In this study, we utilize reactive MD simulations to investigate the process
of fullerene evolution in the circumstellar envelopes (CSEs) of carbon-rich
AGB stars. The CSEs of AGB stars, including dust condensation zone, are well-
known for their abundance of small carbonaceous chains such as C2, making them
primary sites for the production of interstellar carbonaceous dust, and
earning them the designation of “cosmic dust factories” (Clayton et al., 1995;
Bakker et al., 1997; Millar et al., 2000; Santoro et al., 2020). The main
objective of this study is to provide valuable insights into the structural
evolution of fullerenes, driven by the adsorption of carbonaceous chains,
under various simulated conditions that mimic the CSEs of AGB stars.

##  II Methods

Surface adsorption reactions have been identified as a significant mechanism
for the growth of interstellar dust (Carpentier et al., 2012; Burke et al.,
2015). In the bottom-up chemistry approach, the growth of fullerenes is
primarily driven by the adsorption and enrichment of small carbonaceous
molecules. The dust condensation zone within the CSEs of carbon-rich AGB stars
offers a favorable astronomical environment for these adsorption reactions, as
it provides the necessary materials and thermal energy for the reactions to
take place (Decin et al., 2010; Carpentier et al., 2012).

Given that the primary objective of this study was to explore the enrichment
and evolution of C2 small molecules on the surface of stable C60 fullerene, we
conducted simulations to investigate the condensation process of C2 molecules
onto a C60 seed particle. The initial configuration consisted of a periodic
simulation cell with dimensions of $80\times 80\times 80$ Å3, with the C60
molecule placed at its center. A total of 200 C2 molecules
($N_{\text{C}}=400$) were gradually introduced into the cell at random
positions over a duration of $20$ ns. Additionally, a variable number of
neutral hydrogen atoms $N_{\text{H}}$ were included in the system, with
different H concentrations corresponding to $N_{\text{H}}/N_{\text{C}}$ ratios
of 0, 0.05, 0.15, or 0.50.

During the simulations, the majority of the added C2 and H particles adsorbed
onto the surface of the C60 seed particle through van der Waals interactions
(Wang, 2019, 2023). This adsorption process facilitated chemical reactions
between the C2 and H particles, leading to the formation of an initial
structure as depicted in Figure 1. The simulations were conducted at a
specified temperature $T_{\text{ad}}$ to mimic the thermal energy available in
the dust condensation zone within the CSEs of carbon-rich AGB stars. The
thermal energy in this region provides favorable conditions for the adsorption
reactions to take place and contributes to the growth of fullerenes through
the enrichment of carbonaceous molecules.

Figure 1: Schematics for the simulated evolution of C60 via adsorption
reactions.

To simulate the evolution of the formed initial structure, the temperature was
progressively increased from $T_{\text{ad}}$ to $T_{\text{trans}}$ over a
duration of $8$ ns. This temperature increase led to the transformation of the
initial structure into a final structure, whose characteristics depended on
the specific temperature values and the hydrogen concentration. These factors
are known to be influenced by the condition of the CSE and the distance from
the central AGB star (Nowotny et al., 2005; Freytag et al., 2017). Although
the representative temperature near AGB star regions is $2000$ to $3500$ K
(Ziurys, 2006), the dust formation regions in its CSEs are typically found at
distances ranging from $2.5\times 10^{14}$ to $5\times 10^{15}$ cm from the
star’s center, with temperatures estimated to be in the range of $100$ to
$1000$ K (Decin et al., 2010). Therefore, our choice of $T_{\text{ad}}$ fell
within this temperature range. However, we set $T_{\text{trans}}$ to be in the
range of $1500$ to $4000$ K in order to accelerate the reactions, as the
actual timescale of the reactions is much longer than what can be simulated
using MD. Hence, the specific value of $T_{\text{trans}}$ does not hold direct
thermodynamic significance but serves as a parameter to represent the reaction
rate within the constraints of the simulation time.

The interactions between atoms in our simulations were modeled using the
ReaxFF reactive force field (Ashraf and van Duin, 2017). ReaxFF is an
interatomic potential function that takes into account bond-order dependence
and consists of eight terms:
$E_{\text{bond}},E_{\text{over}},E_{\text{under}},E_{\text{lp}},E_{\text{val}},E_{\text{tors}},E_{\text{vdW}},$
and $E_{\text{Coulomb}}$, as shown in Equation 1.

$E_{\text{p}}=E_{\text{bond}}+E_{\text{over}}+E_{\text{under}}+E_{\text{lp}}+E_{\text{val}}+E_{\text{tors}}+E_{\text{vdW}}+E_{\text{Coulomb}}.$
(1)

These terms describe the potential energy of the system with respect to bond
length and order, overcoordination penalty, undercoordination stability, lone
pair, bond angle, single bond torsion, van der Waals interactions, and Coulomb
interactions, respectively. The detailed expression of each term,
parameterization, and benchmarking of the ReaxFF force field can be found in
Ashraf and van Duin (2017). We chose the ReaxFF force field for its proven
success in modeling the formation of fullerenes and related molecules or
nanostructures using MD simulations, as demonstrated in previous studies such
as Martin et al. (2017); Mainitz et al. (2016); Mao et al. (2017); Ostroumova
et al. (2019); Orekhov et al. (2020); Li et al. (2022).

In this study, the simulations were performed using the Large-scale
Atomic/Molecular Massively Parallel Simulator (LAMMPS) (Plimpton, 1995). The
temperature control was achieved using the canonical Nosé-Hoover thermostat
with a time step of $0.1$ fs. It is important to note that, for simplicity, we
considered C2 molecules as the initial reactants, although other carbon
sources, such as larger carbonaceous chains or PAHs, could also be present
(Millar et al., 2000; Li et al., 2008, 2019; Agúndez et al., 2020; Kovács et
al., 2020; Meng et al., 2021; McGuire, 2022; Meng et al., 2023). The ratio
between the number of C and H atoms, $N_{\text{H}}$/$N_{\text{C}}$, was
determined based on the abundance ratio of $\leq 0.03$ between C2 and C60 in
IRC+10216 (Clayton et al., 1995), while the number of hydrogen atoms was set
to balance the formation and photodissociation of C-H bonds (Patra et al.,
2014). In the Supplementary
Data111https://github.com/mengzss/Fullerene_Evolution.git, we have included an
example of the simulation code used in our study, as well as the corresponding
simulation outputs. This will facilitate readers who are interested in
replicating the simulations and conducting further investigations.

##  III Results and discussions

###  III.1 Initial Structures

In our simulations, we observed the formation of different initial structures
with varying shapes under different conditions, as depicted in Figure 2. The
compactness of these structures was characterized using the radius of gyration
$RG$, which is a measure of the distribution of atoms from the center of mass
(Lobanov et al., 2008). A lower $RG$ indicates a more condensed structure for
a given number of atoms. In Panel (a) of Figure 2, the size of the circles
represents the $RG$ values of the initial structures. We observed that at low
concentrations of H, the initial structures exhibited a high degree of
condensation, as indicated by the small circles in the bottom of Panel (a). In
these cases, the adsorbed C2 molecules formed aligned long carbon chains at
low temperatures, as shown in Panel (b). At higher temperatures, we observed
that the chains interconnected to form cage-like structures, as shown in Panel
(c). This is in agreement with the observation of Anders and Urbassek (2021),
which suggests that warm organic cluster particles have a tendency to stick
together through chemical reaction when colliding.

Figure 2: (a) $RG$ of the initial structures formed with different
$N_{\text{H}}$/$N_{\text{C}}$ at various adsorption temperatures
$T_{\text{ad}}$. The bubbles’ size and colors represent the $RG$ value and
different classes of the initial structures, respectively. (b-e) Snapshots of
the formed initial structures.

On the contrary, at high H concentrations, it becomes challenging for the C2
molecules to form the same types of structures. This is due to many C2
molecules being saturated with H atoms, which occupy the chemically active
sites and hinder the formation of interconnected chains. This observation
aligns with the findings of Cherchneff (2011), who observed from the molecular
abundance of IRC+10216 that excess hydrogen led to the formation of
hydrocarbon molecules rather than carbon chains. As a result, the C2 form
small hydrocarbon molecules such as C2H2, C2H4, C4H6, and C5H5, which remain
as individual entities. These molecules are primarily held together by vdW
forces on the surface of the C60 at low temperatures. However, when the
temperature exceeds the binding energy between these molecules, they can no
longer remain bound, leading to their dissociation, as shown in Panel (e) of
Figure 2. The case depicted in Panel (d) represents an intermediate scenario
between the condensed structure and the dissociated state. It is worth noting
that previous studies have demonstrated that excessive hydrogenation can lead
to the dissociative fragmentation of molecules (Tang et al., 2022; Gatchell et
al., 2015).

###  III.2 Final Structures

We explore further possible evolution of the condensed initial structures. In
the CSEs of AGB stars, the dust particles experience various astronomical
events such as photoprocessing, stellar winds, and shocks, which could result
in a sharp increase in their temperature and induce heterogeneous reactions
(Jones et al., 2013). In our simulations, after the initial structures were
formed through adsorption, we subjected them to high temperatures to explore
their potential structural transformations. Our results revealed two distinct
evolution pathways leading to different types of final structures. In the case
of low H concentration, onion-like nanostructures consisting of two layers of
fullerenes were observed, as depicted in the top panels of Figure 3.
Conversely, for the case of high H concentration, single-layer fullerenes were
formed, as shown in the bottom panels of Figure 3.

Figure 3: Final structures formed with different H concentrations through
heterogeneous reactions at different temperatures. The gray, red, and blue
dots represent H and C atoms originating from C60 and C2, respectively.

It is worth noting that the specific temperature at which the transformation
occurs does not definitively determine the type of final structure, except in
the case of $2500$ K, which appears insufficient to trigger the formation of
fullerene within the simulated time scale. This finding is consistent with the
observation of Cami et al. (2018), which shows that PAHs become the main
product of evolution in the hydrogen-rich environment at temperatures between
$1000$ and $1700$ K. Moreover, it is important to acknowledge that fullerene
transformation could still potentially occur at that temperature, taking into
account that the actual reaction time is significantly longer than what is in
the MD study despite of a much lower concentration.

The C-H bonds played a key role in determining the final structures. The plot
in Figure 4 illustrates the evolution of the number of C-C and C-H bonds,
denoted as $N_{\text{\ce{C-C}}}$ and $N_{\text{\ce{C-H}}}$, revealing three
distinct stages of structural transformation. Initially, during the first
$0.8$ nanoseconds, the bond numbers remain relatively stable. In the second
stage, a notable increase in $N_{\text{\ce{C-C}}}$ occurs as the temperature
surpasses a specific threshold ($1500$ K), indicating a significant structural
transition. This transition involves the connection of adsorbed molecules
through unoccupied sites, with the highest bonding rate observed in the
absence of hydrogen. This process is similar to the ring-condensation process
outlined in Irle et al. (2003a), marking the phase in which carbon rings
commence formation. Subsequently, at around $1.5$ nanoseconds and at a
temperature of approximately $2500$ K, a substantial breakage of C-H bonds
commences. This leads to the generation of new unoccupied carbon sites
available for bonding, which typically promotes the formation of new carbon
rings (Saha and Schatz, 2012).

Figure 4: Evolution of C-C (a) and C-H bonds (b) during heating to
$T_{\text{trans}}=4000$ K for different H concentrations. The background
colors indicate three stages of structural transformation. (c) Configurations
of the simulated system at different time for the case of
$N_{\text{H}}/N_{\text{C}}=0.5$.

The selection between the two distinct types of final structures predominantly
occurred during the transition from the second to the third stage. During this
stage of evolution, the predominant process involves cage closure,
characterized by the collapse of a sizable ring that undergoes a closing
mechanism (Irle et al., 2003a), while the closure of a double-layered
fullerene was not observed. In scenarios with low hydrogen concentration, the
majority of adsorbed carbon atoms formed bonds with each other in the second
stage, resulting in the creation of a relatively stable layer that enveloped
the central C60. The temperature during this phase was insufficient to
facilitate the transformation of intra-layer $sp^{2}$ bonds into inter-layer
$sp^{3}$ bonds. Conversely, in high hydrogen concentration scenarios, the
breakage of C-H bonds in the third stage led to the emergence of numerous
unoccupied carbon sites. As the temperature increased in the third stage,
these sites interacted with the central C60, eventually causing it to dissolve
into a larger, single-layered fullerene cage. This process is depicted in
Panel (c) of Figure 4 taking the case of $T_{\text{trans}}=4000$ K for
example.

During the process of compact cage contraction, we observed events of small
molecule pop-out, involving species like C2 and H2, referred to as the
“shrinking hot giant” road in Irle et al. (2007). It is worth noting that in
cases of low hydrogen concentration, when the temperature was high enough, the
C60 could also react with atoms in the outer layer; however, the two-layered
structures were ultimately maintained. While direct evidence is lacking, it is
plausible that dissociated hydrogen atoms played a catalytic role in the
opening of the central C60 (Yang et al., 2023). We note that the formation of
multilayer fullerenes studied here differs from the previous high-density pure
carbon gas ultrafast cooling mechanism (Ostroumova et al., 2019). The
adsorption-driven process studied here is more likely to occur in the long-
term evolution of low-density interstellar environments.

The formation of H2 molecules in the ISM remains a challenging question in
astronomy (Wakelam et al., 2017). The exothermic nature of the combination of
hydrogen atoms to form H2 poses a dilemma, as H2 lacks an efficient mechanism
to dissipate the released heat in the low-density environment of the ISM.
Current understanding suggests that the primary route for interstellar H2
formation is through the Eley-Rideal reaction on dust or ice particles, which
act as heat reservoirs (Bauschlicher Jr., 1998; Mennella et al., 2012). These
particles play a crucial role in H2 formation by facilitating the dissipation
of heat (Gould and Salpeter, 1963; Lemaire et al., 2010; Thrower et al., 2012;
Mennella et al., 2012; Boschman et al., 2015; Foley et al., 2018; Barrales-
Martinez and Gutierrez-Oliva, 2019; Pantaleone et al., 2021). Conversely,
hydrogen can also facilitate the formation of complex organic molecules (Lu et
al., 2021; Yang et al., 2023). In the simulations conducted for this study,
the formation of H2 molecules was observed as by-products during the growth of
fullerenes. It was found that the formation of H2 is closely correlated with
the breaking of C-H bonds, as shown in Figure 5, suggesting a potential
catalytic role of evolving fullerenes in the H2 formation process. However, it
is important to note that the specific question regarding heat dissipation
during H2 formation was not directly addressed in these simulations, as a
global thermostat was applied. To fully understand the intricacies of heat
dissipation in the H2 formation process, a separate set of simulations
specifically designed to explore this aspect and more in-depth analyses would
be necessary.

Figure 5: (a) Ratio of the number of H2 molecules formed over the total number
of H atoms initially present in the simulation cell as a function of time, for
different H concentrations. (b) Time evolution of the C-H bonds in the system.

###  III.3 Infrared spectrum

To explore the impact of the evolution towards the two different types of
fullerenes on their infrared spectra, we employed density functional theory
(DFT) for geometry optimizations and frequency calculations (Banhatti et al.,
2021). The calculations were conducted at the B3LYP/6-31G(d) level of theory
for three selected samples: a C60 molecule, a C180 molecule, and a two-layered
C60@C180 composite, as depicted in the inset of Figure 6. We chose these
smaller molecules for simplicity, considering the significant computational
challenges associated with performing calculations on larger simulated final
structures. Given the large size of the system, we computed the spectra under
the assumption of harmonic vibrations. To account for anharmonicity effects,
we uniformly scaled the calculated spectra by a factor of 0.9613, as proposed
by Borowski (2012). This scaling factor enables us to approximate the
influence of anharmonicity on the spectra and improve their accuracy.
Additionally, for an accurate description of the interaction between different
layers of fullerenes (Saal et al., 2009), we included the D3 version of
Grimme’s empirical dispersion interaction (Grimme et al., 2010), along with
the Becke-Johnson damping function (Grimme et al., 2011).

Figure 6: (a) Calculated IR spectrum of C60, C180, and C60@C180. The dashed
line represents the sum of the blue and orange lines. (b) Visualizations of
the vibration modes of C180 at 14.73, 19.93, and 29.86$\,\upmu$m. The arrows
indicate the vibration directions of atoms.

Figure 6 (a) displays the computed IR spectra of the three samples. It is
noted that when comparing C180 with C60, additional emission features are
observed at specific wavelengths ($10.84$, $14.65$, and $29.86\,\upmu$m) for
C180. The complex vibrational behavior of C180 indicated by the additional
features aligned with previous observation of Bakowies and Thiel (1991); Cami
et al. (2010). Upon comparing the summed spectrum of C60 and C180 (purple,
dashed line) with that of C60@C180 (green line), notable observations are
made. Specifically, the $14.73$, $19.93$, and $29.86\,\upmu$m bands of the
summed spectrum exhibit hindrance or a blue-shifted pattern when the two
molecules are combined. This alteration in the IR features, characterized by a
decrease in intensity accompanied by an increase in frequency, can typically
be attributed to the increased rigidity of the molecule in the direction of
vibration when the C60 molecule is encapsulated within C180. This suggests
that these bands likely correspond to vibrations perpendicular to the surface.
A more detailed analysis of the vibrational modes associated with these bands,
as depicted in Figure 6 (b), provides further confirmation of this hypothesis.
Therefore, the blue-shifted behavior of these bands holds particular interest
for astronomical observations targeting onion-like fullerenes.

##  IV Conclusions

In conclusion, this study employed MD simulations to investigate the evolution
of fullerenes in environments resembling those found in the CSEs of AGB stars.
The simulations revealed two distinct types of transformations that occur when
fullerenes adsorb carbon chains, resulting in the formation of single- and
double-layered carbon nanostructures. The formation of single-layered
structures was observed under high H concentration, while the formation of
double-layered structures occurred under low H concentration. We note that the
simulations did not investigate the transformation of fullerenes into larger
sizes with additional layers due to system size limitations. Additionally, the
study found that the rate of H2 molecule formation is closely linked to the
breaking of C-H bonds in the evolving fullerene, implying a potential
catalytic role of fullerene in H2 formation.

Furthermore, DFT calculations were conducted on three fullerene molecules to
investigate the impact of structural evolution on IR emission features. The
results demonstrated that as fullerenes undergo structural evolution towards
larger sizes, additional emission features emerge at specific wavelengths,
such as $10.84$, $14.65$, and $29.86\,\upmu$m. Additionally, the presence of a
double layer in the fullerene structure led to noticeable blue-shift or
weakening effects on the bands at $14.73$, $19.93$, and $29.86\,\upmu$m. These
findings provide valuable insights into the structural changes and spectral
characteristics of fullerenes, suggesting the possibility of fullerene growth
in the ISM through the adsorption of small carbon molecules and subsequent
heterogeneous reactions.

## Acknowledgments

The authors acknowledge financial support from the National Natural Science
Foundation of China (11964002). Prof. Yong Zhang is acknowledged for fruitful
discussion.

## Data Availability

The Supplementary Data include an example of the LAMMPS input script that was
used in our study, along with the corresponding simulation outputs. They are
available at https://github.com/mengzss/Fullerene_Evolution.git.

## References

  * (1)
  * Kroto et al. (1985) H. W. Kroto, J. R. Heath, S. C. O’Brien, R. F. Curl, and R. E. Smalley, Nature 318, 162 (1985).
  * Kratschmer et al. (1990a) W. Kratschmer, K. Fostiropoulos, and D. R. Huffman, Chem. Phys. Lett. 170, 167 (1990a).
  * Kratschmer et al. (1990b) W. Kratschmer, L. D. Lamb, K. Fostiropoulos, and D. R. Huffman, Nature 347, 354 (1990b).
  * Taylor et al. (1990) R. Taylor, J. P. Hare, A. K. Abdul-Sada, and H. W. Kroto, J. Chem. Soc., Chem. Commun. 24, 1423 (1990).
  * Cami et al. (2010) J. Cami, J. Bernard-Salas, E. Peeters, and S. E. Malek, Science 329, 1180 (2010).
  * Gielen et al. (2011) C. Gielen, J. Cami, J. Bouwman, E. Peeters, and M. Min, A&A 536, A54 (2011).
  * Berne et al. (2013) O. Berne, G. Mulas, and C. Joblin, A&A 550, L4 (2013).
  * Sellgren et al. (2010) K. Sellgren, M. W. Werner, J. G. Ingalls, J. D. T. Smith, T. M. Carleton, and C. Joblin, ApJL 722, L54 (2010).
  * Garcia-Hernandez et al. (2010) D. A. Garcia-Hernandez, A. Manchado, P. Garcia-Lario, L. Stanghellini, E. Villaver, R. A. Shaw, R. Szczerba, and J. V. Perea-Calderon, ApJL 724, L39 (2010).
  * Roberts et al. (2012) K. R. G. Roberts, K. T. Smith, and P. J. Sarre, MNRAS 421, 3277 (2012).
  * Castellanos et al. (2014) P. Castellanos, O. Berné, Y. Sheffer, M. G. Wolfire, and A. G. G. M. Tielens, ApJ 794, 83 (2014).
  * Berné et al. (2017) O. Berné, N. L. J. Cox, G. Mulas, and C. Joblin, A&A 605, L1 (2017).
  * Candian et al. (2018) A. Candian, A. G. G. M. Tielens, and J. Zhen, Phys. Today 71, 38 (2018).
  * Otsuka et al. (2014) M. Otsuka, F. Kemper, J. Cami, E. Peeters, and J. Bernard-Salas, MNRAS 437, 2577 (2014).
  * Dunk et al. (2013) P. W. Dunk, J.-J. Adjizian, N. K. Kaiser, J. P. Quinn, G. T. Blakney, C. P. Ewels, A. G. Marshall, and H. W. Kroto, PNAS 110, 18081 (2013).
  * García-Hernández et al. (2013) D. A. García-Hernández, F. Cataldo, and A. Manchado, MNRAS 434, 415 (2013).
  * Berné and Tielens (2012) O. Berné and A. G. G. M. Tielens, PNAS 109, 401 (2012).
  * Chuvilin et al. (2010) A. Chuvilin, U. Kaiser, E. Bichoutskaia, N. A. Besley, and A. N. Khlobystov, Nat. Chem. 2, 450 (2010).
  * Zhang et al. (2013) J. Zhang, F. L. Bowles, D. W. Bearden, W. K. Ray, T. Fuhrer, Y. Ye, C. Dixon, K. Harich, R. F. Helm, M. M. Olmstead, et al., Nat. Chem. 5, 880 (2013).
  * Pietrucci and Andreoni (2014) F. Pietrucci and W. Andreoni, J. Chem. Theory Comput. 10, 913 (2014).
  * Micelotta et al. (2012) E. R. Micelotta, A. P. Jones, J. Cami, E. Peeters, J. Bernard-Salas, and G. Fanchini, ApJ 761, 35 (2012).
  * Montillaud et al. (2013) J. Montillaud, C. Joblin, and D. Toublanc, A&A 552, A15 (2013).
  * Berné et al. (2015) O. Berné, J. Montillaud, and C. Joblin, A&A 577, A133 (2015).
  * Patra et al. (2014) N. Patra, P. Král, and H. Sadeghpour, ApJ 785, 6 (2014).
  * Cherchneff et al. (2000) I. Cherchneff, Y. H. Le Teuff, P. M. Williams, and A. G. G. M. Tielens, A&A 357, 572 (2000).
  * Pascoli and Polleux (2000) G. Pascoli and A. Polleux, A&A 359, 799 (2000).
  * Bettens and Herbst (1997) R. P. A. Bettens and E. Herbst, ApJ 478, 585 (1997).
  * Jones et al. (2011) B. M. Jones, F. Zhang, R. I. Kaiser, A. Jamal, A. M. Mebel, M. A. Cordiner, and S. B. Charnley, PNAS 108, 452 (2011).
  * Parker et al. (2012) D. S. N. Parker, F. Zhang, Y. S. Kim, R. I. Kaiser, A. Landera, V. V. Kislov, A. M. Mebel, and A. G. G. M. Tielens, PNAS 109, 53 (2012).
  * Bates and Scuseria (1997) K. R. Bates and G. E. Scuseria, J. Phys. Chem. A 101, 3038 (1997).
  * Irle et al. (2003a) S. Irle, G. Zheng, M. Elstner, and K. Morokuma, Nano Lett. 3, 465 (2003a).
  * Irle et al. (2003b) S. Irle, G. Zheng, M. Elstner, and K. Morokuma, Nano Lett. 3, 1657 (2003b).
  * Zheng et al. (2004) G. Zheng, S. Irle, M. Elstner, and K. Morokuma, J. Phys. Chem. A 108, 3182 (2004).
  * Zheng et al. (2005) G. Zheng, S. Irle, and K. Morokuma, J. Chem. Phys. 122, 014708 (2005).
  * Zheng et al. (2007) G. Zheng, Z. Wang, S. Irle, and K. Morokuma, J. Nanosci. Nanotechnol. 7, 1662 (2007).
  * Irle et al. (2006) S. Irle, G. Zheng, Z. Wang, and K. Morokuma, J. Phys. Chem. B 110, 14531 (2006).
  * Irle et al. (2007) S. Irle, G. Zheng, Z. Wang, and K. Morokuma, Nano 2, 21 (2007).
  * Saha et al. (2009) B. Saha, S. Shindo, S. Irle, and K. Morokuma, Acs Nano 3, 2241 (2009).
  * Zhen et al. (2019) J. Zhen, W. Zhang, Y. Yang, and Q. Zhu, MNRAS 490, 3498 (2019).
  * Parker and Kaiser (2017) D. S. N. Parker and R. I. Kaiser, Chem. Soc. Rev. 46, 452 (2017).
  * Fioroni et al. (2019) M. Fioroni, R. E. Savage, and N. J. DeYonker, Phys. Chem. Chem. Phys. 21, 8015 (2019).
  * Meuwly (2019) M. Meuwly, WIREs Comput. Mol. Sci. 9, e1386 (2019).
  * Montillaud and Joblin (2014) J. Montillaud and C. Joblin, A&A 567, A45 (2014).
  * Qi et al. (2018) H. Qi, S. Picaud, M. Devel, E. Liang, and Z. Wang, ApJ 867, 133 (2018).
  * Chen (2018) T. Chen, ApJ 866, 113 (2018).
  * Chen et al. (2020) T. Chen, Y. Luo, and A. Li, A&A 633, A103 (2020).
  * Ostroumova et al. (2019) G. Ostroumova, N. Orekhov, and V. Stegailov, Diam. Relat. Mater. 94, 14 (2019).
  * Hanine et al. (2020) M. Hanine, Z. Meng, S. Lu, P. Xie, S. Picaud, M. Devel, and Z. Wang, ApJ 900, 188 (2020).
  * Orekhov et al. (2020) N. Orekhov, G. Ostroumova, and V. Stegailov, Carbon 170, 606 (2020).
  * Anders and Urbassek (2021) C. Anders and H. M. Urbassek, A&A 647, L13 (2021).
  * Tang et al. (2022) Z. Tang, F. D. S. Simonsen, R. Jaganathan, J. Palotás, J. Oomens, L. Hornekær, and B. Hammer, A&A 663, A150 (2022).
  * Clayton et al. (1995) G. C. Clayton, D. M. Kelly, J. H. Lacy, I. R. Little-Marenin, P. A. Feldman, and P. F. Bernath, ApJ 109, 2096 (1995).
  * Bakker et al. (1997) E. Bakker, E. van Dishoeck, L. Waters, and T. Schoenmaker, A&A 323, 469 (1997).
  * Millar et al. (2000) T. Millar, E. Herbst, and R. Bettens, MNRAS 316, 195 (2000).
  * Santoro et al. (2020) G. Santoro, L. Martínez, K. Lauwaet, M. Accolla, G. Tajuelo-Castilla, P. Merino, J. M. Sobrado, R. J. Peláez, V. J. Herrero, I. Tanarro, et al., ApJ 895, 97 (2020).
  * Carpentier et al. (2012) Y. Carpentier, G. Féraud, E. Dartois, R. Brunetto, E. Charon, A.-T. Cao, L. d’Hendecourt, P. Bréchignac, J.-N. Rouzaud, and T. Pino, A&A 548, A40 (2012).
  * Burke et al. (2015) D. J. Burke, F. Puletti, P. M. Woods, S. Viti, B. Slater, and W. A. Brown, J. Phys. Chem. A 119, 6837 (2015).
  * Decin et al. (2010) L. Decin, M. Agúndez, M. J. Barlow, F. Daniel, J. Cernicharo, R. Lombaert, E. De Beck, P. Royer, B. Vandenbussche, R. Wesson, et al., Nature 467, 64 (2010).
  * Wang (2019) Z. Wang, J. Phys. Chem. C 123, 15166 (2019).
  * Wang (2023) Z. Wang, J. Phys. Chem. B 127, 6227 (2023).
  * Nowotny et al. (2005) W. Nowotny, B. Aringer, S. Höfner, R. Gautschy-Loidl, and W. Windsteig, A&A 437, 273 (2005).
  * Freytag et al. (2017) B. Freytag, S. Liljegren, and S. Höfner, A&A 600, A137 (2017).
  * Ziurys (2006) L. M. Ziurys, PNAS 103, 12274 (2006).
  * Ashraf and van Duin (2017) C. Ashraf and A. C. van Duin, J. Phys. Chem. A 121, 1051 (2017).
  * Martin et al. (2017) J. W. Martin, G. J. McIntosh, R. Arul, R. N. Oosterbeek, M. Kraft, and T. Söhnel, Carbon 125, 132 (2017).
  * Mainitz et al. (2016) M. Mainitz, C. Anders, and H. M. Urbassek, A&A 592, A35 (2016).
  * Mao et al. (2017) Q. Mao, A. C. T. van Duin, and K. Luo, Carbon 121, 380 (2017).
  * Li et al. (2022) S. Li, F. Bian, X. Meng, D. Zhai, H. Yang, and G. Qin, Carbon 189, 484 (2022).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Li et al. (2008) A. Li, J. H. Chen, M. P. Li, Q. J. Shi, and Y. J. Wang, MNRAS 390, L39 (2008).
  * Li et al. (2019) Q. Li, A. Li, and B. W. Jiang, MNRAS 490, 3875 (2019).
  * Agúndez et al. (2020) M. Agúndez, J. Martínez, P. de Andres, J. Cernicharo, and J. A. Martín-Gago, A&A 637, A59 (2020).
  * Kovács et al. (2020) P. Kovács, X. Zhu, J. Carrete, G. K. Madsen, and Z. Wang, ApJ 902, 100 (2020).
  * Meng et al. (2021) Z. Meng, X. Zhu, P. Kovács, E. Liang, and Z. Wang, ApJ 922, 101 (2021).
  * McGuire (2022) B. A. McGuire, ApJS 259, 30 (2022).
  * Meng et al. (2023) Z. Meng, Y. Zhang, E. Liang, and Z. Wang, MNRAS 525, L29 (2023).
  * Lobanov et al. (2008) M. Y. Lobanov, N. S. Bogatyreva, and O. V. Galzitskaya, Mol. Biol. 42, 623 (2008).
  * Cherchneff (2011) I. Cherchneff, EAS Publ. Ser. 46, 177 (2011).
  * Gatchell et al. (2015) M. Gatchell, M. H. Stockett, N. de Ruette, T. Chen, L. Giacomozzi, R. F. Nascimento, M. Wolf, E. K. Anderson, R. Delaunay, V. Vizcaino, et al., Phys. Rev. A 92, 050702 (2015).
  * Jones et al. (2013) A. P. Jones, L. Fanciullo, M. Köhler, L. Verstraete, V. Guillet, M. Bocchio, and N. Ysard, A&A 558, A62 (2013).
  * Cami et al. (2018) J. Cami, E. Peeters, J. Bernard-Salas, G. Doppmann, and J. De Buizer, Galaxies 6, 101 (2018).
  * Saha and Schatz (2012) B. Saha and G. C. Schatz, J. Phys. Chem. B 116, 4684 (2012).
  * Yang et al. (2023) S. Yang, P. Xie, E. Liang, and Z. Wang, RAA 23, 055019 (2023).
  * Wakelam et al. (2017) V. Wakelam, E. Bron, S. Cazaux, F. Dulieu, C. Gry, P. Guillard, E. Habart, L. Hornekær, S. Morisset, G. Nyman, et al., Mol. Astrophys. 9, 1 (2017).
  * Bauschlicher Jr. (1998) C. W. Bauschlicher Jr., ApJ 509, L125 (1998).
  * Mennella et al. (2012) V. Mennella, L. Hornekr, J. Thrower, and M. Accolla, ApJ 745, L2 (2012), ISSN 2041-8205, 2041-8213.
  * Gould and Salpeter (1963) R. J. Gould and E. E. Salpeter, ApJ 138, 393 (1963).
  * Lemaire et al. (2010) J. L. Lemaire, G. Vidali, S. Baouche, M. Chehrouri, and H. Mokrane, ApJL 725, L156 (2010).
  * Thrower et al. (2012) J. D. Thrower, B. Jorgensen, E. E. Friis, S. Baouche, V. Mennella, A. C. Luntz, M. Andersen, B. Hammer, and L. Hornekr, ApJ 752, 3 (2012).
  * Boschman et al. (2015) L. Boschman, S. Cazaux, M. Spaans, R. Hoekstra, and T. Schlathlter, A&A 579, A72 (2015), ISSN 0004-6361, 1432-0746.
  * Foley et al. (2018) N. Foley, S. Cazaux, D. Egorov, L. M. P. V. Boschman, R. Hoekstra, and T. Schlathlter, MNRAS 479, 649 (2018), ISSN 0035-8711, 1365-2966.
  * Barrales-Martinez and Gutierrez-Oliva (2019) C. Barrales-Martinez and S. Gutierrez-Oliva, MNRAS 490, 172 (2019), ISSN 0035-8711, 1365-2966.
  * Pantaleone et al. (2021) S. Pantaleone, J. Enrique-Romero, C. Ceccarelli, S. Ferrero, and P. Ugliengo, ApJ 917, 49 (2021).
  * Lu et al. (2021) S. Lu, Z. Meng, P. Xie, E. Liang, and Z. Wang, A&A 656, A84 (2021).
  * Banhatti et al. (2021) S. Banhatti, J. Palotás, P. Jusko, B. Redlich, J. Oomens, S. Schlemmer, and S. Brünken, A&A 648, A61 (2021).
  * Borowski (2012) P. Borowski, J. Phys. Chem. A 116, 3866 (2012).
  * Saal et al. (2009) A. Saal, T. Jarrosson, O. Ouamerali, and C. A. Daul, Chem. Phys. Lett. 480, 225 (2009).
  * Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132 (2010).
  * Grimme et al. (2011) S. Grimme, S. Ehrlich, and L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
  * Bakowies and Thiel (1991) D. Bakowies and W. Thiel, Chem. Phys. 151, 309 (1991).


# Eco-Coasting Strategies Using Road Grade Preview: Evaluation and Online
Implementation Based on Mixed Integer Model Predictive Control

Yongjun Yan, Nan Li, Jinlong Hong, Bingzhao Gao, Hong Chen, Jing Sun, and
Ziyou Song Y. Yan is with the State Key Laboratory of Automotive Simulation
and Control, Jilin University, Changchun 130022, China. This work was done
when he is a visiting Ph.D. student in the Department of Naval Architecture
and Marine Engineering, University of Michigan, Ann Arbor, MI 48109 USA
(e-mail: yanyj18@mails.jlu.edu.cn).N. Li is with the Department of Aerospace
Engineering, University of Michigan, Ann Arbor, MI 48109 USA (e-mail:
nanli@umich.edu)).J. Hong, and B. Gao are with the Clean Energy Automotive
Engineering Center, Tongji University, Shanghai, 201804, China and also with
the State Key Laboratory of Automotive Simulation and Control, Jilin
University, Changchun 130022, China (e-mail<EMAIL_ADDRESS>gaobz@jlu.edu.cn).H. Chen is with the College of Electronic and Information
Engineering, Tongji University, Shanghai 201804, China (e-mail:
chenhong2019@tongji.edu.cn).J. Sun is with the Department of Naval
Architecture and Marine Engineering, University of Michigan, Ann Arbor, MI
48109 USA (e-mail: jingsun@umich.edu).Z. Song is with the Department of
Mechanical Engineering, National University of Singapore, Singapore 117575,
Singapore and also with the Department of Naval Architecture and Marine
Engineering, University of Michigan, Ann Arbor, MI 48109 USA (e-mail:
Ziyou@nus.edu.sg).

###### Abstract

Coasting has been widely used in the eco-driving guidelines to reduce fuel
consumption by profiting from kinetic energy. However, the comprehensive
comparison between different coasting strategies and online performance of the
eco-coasting strategy using road grade preview are still unclear because of
the oversimplification and the integer variable in the optimal control
problems. Herein, two different coasting strategies (fuel cut-off and engine
start/stop) are proposed to reveal the potential benefit of eco-coasting using
the road grade preview. Engine drag torque and energy cost used for engine
restart are considered in the modeling to give a fair evaluation of the
offline and online performance. The offline performance of these two coasting
methods is evaluated through dynamic programming (DP) under various driving
scenarios with different slope profiles. Offline simulation shows that the
engine start/stop method outperforms the fuel cut-off method in terms of fuel
consumption and travel time by getting rid of the engine drag torque. Then,
online performance of these two coasting methods is evaluated using Mixed
Integer Model Predictive Control (MIMPC). A novel operational constraint on
the minimum off steps is added in the MIMPC formulation to avoid frequent
switch of the integer variables which represent the fuel cut-off and the
engine start/stop mechanism. Simulation results show that, for both fuel cut-
off and engine start/stop coasting methods, the MPC controller reduces fuel
consumption to a level comparable to DP without sacrificing the travel time.

###### Index Terms:

Mixed integer model predictive control; Road grade; Eco-coasting strategy.

## Nomenclature

$A$ |  Frontal projected area of vehicle, [$m^{2}$].
---|---
$C$ |  Vector of constraints, [$-$].
$C_{d}$ |  Coefficient of drag, [$-$].
$d_{\text{min}}$ |  Minimum off steps, [$-$].
$d$ |  Engine start/stop signal, [$-$].
$F_{t}$ |  Traction force, [$N$].
$F_{r}$ |  Resistance force, [$N$].
$f_{r}$ |  Coefficient of rolling resistance, [$-$].
$g$ |  Gravitational constant, [$-$].
$I_{\text{engine}}$ |  Engine inertia, [$kg\,m^{2}$].
$I_{\text{final}}$ |  Gear ratio of final reduction drive, [$-$].
$I_{g}$ |  Gear ratio of gearbox, [$-$].
$M_{f}$ |  Fuel consumption model in the distance
|  domain, [$-$].
$m_{\text{eff}}$ |  Effective mass of the vehicle, [$kg$].
$\dot{m}_{{f}}$ |  Instantaneous fuel consumption rate, [$g/s$].
$N_{h}$ |  Prediction horizon, [$-$].
$r$ |  Wheel radius, [$m$].
$s$ |  Space coordinate, [$m$].
$T_{b}$ |  Brake torque, [$Nm$].
$T_{e}$ |  Engine torque, [$Nm$].
$T_{\text{final,1}}$ |  Final driveshaft torque of baseline
|  method, [$Nm$].
$T_{\text{final,2}}$ |  Final driveshaft torque of fuel cut-off
|  method, [$Nm$].
$T_{\text{final,3}}$ |  Final driveshaft torque of engine
|  start/stop method, [$Nm$].
$v$ |  Vehicle longitudinal velocity, [$m/s$].
$\bar{v}$ |  Partial velocity variation, [$m/s$].
$z$ |  fuel cut-off signal, [$-$].
$z_{\text{current}}$ |  Optimal binary variables among the current
|  horizon, [$-$].
$z_{\text{history}}$ |  History binary variables, [$-$].
$\alpha_{1\sim 4}$ |  Coefficients of fuel consumption model, [$-$].
$\beta$ |  Weight of weighted sum method, [$-$].
$\Delta E_{\text{engine }}$ |  Energy required for cranking the engine, [$J$].
$\Delta{s}$ |  Discretization distance, [$-$].
$\Delta{v}$ |  Velocity reduction, [$m/s$].
$\delta E_{v}$ |  Partial inertia energy, [$J$].
$\delta$ |  Extended 0/1 variables, [$-$].
$\varepsilon$ |  Tolerance constant, [$-$].
$\eta$ |  Gearbox efficiency, [$-$].
$\phi$ |  Slope data, [∘].
$\rho$ |  Air density, [$[kg/m^{3}]$].
$\omega_{e}$ |  Rotational velocity of the engine, [$rpm$].
$\omega_{w}$ |  Rotational velocity of wheel, [$rpm$].

## Acronyms

$BSFC$ |  Brake Specific Fuel Consumption.
---|---
$DP$ |  Dynamic Programming.
$FCO$ |  Fuel Cut-Off.
$MPC$ |  Model Predictive Control.
$MIMPC$ |  Mixed Integer Model Predictive Control.
$MIOCP$ |  Mixed Integer Optimal Control Problems.
$SUV$ |  Sports Utility Vehicle.
$TCU$ |  Transmission Control Unit.

##  I Introduction

For on-road vehicles, the key to further improve fuel economy with eco-driving
is comprehensive incorporation and full integration of the powertrain system
with the driving environment represented by traffic conditions and digital
maps. Eco-coasting, as a category of eco-driving, refers to the strategies to
roll the vehicle with kinetic energy without traction force [1]. For a
conventional car with a gasoline engine and automatic transmission, there are
three methods to coast during the deceleration phases: setting the gear in the
neutral position and turning the engine to idle [2, 3]; shutting fuel
injection off when no torque is requested [4, 5, 6]; manipulating the lock-up
clutch to disconnect the engine from powertrain when the engine is off [7, 8,
9, 10, 11].

The majority of studies on eco-coasting have focused on the fuel economy
benefit evaluation of different coasting strategies with backward vehicle
simulation approaches [2, 12, 13], where the required engine load is
calculated from the pre-defined vehicle speed profile. However, it does not
fully explore the potential of eco-coasting as the generated coasting action
only relied on current vehicle state and road information. Previous
investigation has proven that the look-ahead control can achieve better
performance when the slope information on the road ahead is available [14]. To
the author’s best knowledge, no research has been reported on the comparison
of the fuel saving potential among different coasting methods under the
cruising condition with road grade preview, let alone the online
implementation of these coasting strategies.

Model predictive control (MPC) is a promising method for the online
implementation of eco-coasting strategies. However, the discrete nature of the
engine start/stop and fuel cut-off coasting methods render the mixed-integer
optimal control problems (MIOCPs). Similar to the gear shift optimization [15,
16] or engine start/stop of hybrid vehicles [15, 17] in the online energy
management field. There are two open problems for the online implementation of
the eco-coasting strategies with MPC controller: (1) frequent switch of the
integer control variables; and (2) computational burden. To mitigate the
frequent switch of the integer signal, an additional penalty of the switching
frequency [15] and inequality constraints implying that the limitation of
total switch numbers of the integer variables is imposed in the optimization
problems [18]. These formulations yield a reduction in the switching frequency
for a finite horizon open-loop optimal control problem. However, neither the
additional penalty nor the inequality constraints is effective for MPC when
the first element of the control sequence will change at different sampling
instant if the solution is numerically sensitive to different driving
conditions. In [19, 20], mixed-integer linear inequalities are applied in the
MPC formulation, modeling the minimum amount of time for which a on/off
control variable must be kept on/off.

As to the computational complexity, generally, there are two approaches to
solve the MINLP: the first one is to solve the MINLP directly with branch &
bound method or the related solvers that require excessive computational
resource [21]; the second one is the heuristic method which reformulates the
original MINLP into NLP by relaxing the integer variables as continuous ones,
then round the feasible solution to the nearest integer. This relax-and-round
method is computation efficient but can not guarantee the optimality of the
final solution [22, 23]. In [24, 25, 26, 27], engine start/stop signal is
calculated in the upper layer with Dynamic Programming (DP). With a heuristic
engine start/stop strategies, the original mixed-integer hybrid powertrain
optimization problem can be converted into a simple convex optimization
problem. Therefore, the problems referring to the frequent switch and
computational burden are mitigated, but the robustness issues arising from the
heuristic feature is challenging.

In this paper, instead of assuming the drag torque is zero in the fuel cut-off
method [28] or ignoring the extra cost to restart the engine [29] for energy
management research, two coasting models are proposed to capture the fuel cut-
off and engine start/stop mechanisms. Comprehensive performance concerning the
travel time and fuel consumption of these two coasting methods are evaluated
by solving a multi-objective optimization problem with DP using different
slope profiles which are known as a priori. Finally, the online implementation
of the aforementioned coasting methods is investigated using an MPC framework.

The contributions of this paper are as follows:
1) We quantify the potential benefits of these two coasting methods based on
results derived by DP under different slope profiles from digital maps.
2) We assess the online implementation for these two coasting methods with the
MIMPC scheme compared with the PI controller.
3) A novel MIMPC formulation with inequality constraints applied to the
integer variables is used to mitigate the frequent switch of the coasting
signal.
4) Through simulations, we demonstrate that the computation burden can be
reduced dramatically using offline calculated coasting signal when the slope
profile and reference speed are fixed in advance.

The remainder of this paper is organized as follows: Section II discusses
different vehicle coasting models; Section III describes the offline
evaluation of the proposed coasting strategies; Section IV presents the online
implementation of two proposed coasting methods with MIMPC controller;
Finally, conclusions and future work are concluded in Section V.

###  I-A Vehicle Model

Consider the longitudinal vehicle model described as

$\displaystyle\dot{s}$ $\displaystyle=v,$ (1) $\displaystyle\dot{v}$
$\displaystyle=\frac{1}{m_{\text{eff}}}\left[F_{t}-F_{r}(s,v)\right],$

where $s$ is the distance traveled and $v$ is the velocity, $\dot{s}$ and
$\dot{v}$ are the derivative with respect to time $t$, $m_{\text{eff}}$ is the
total mass of the vehicle, $F_{t}$ is the traction force and $F_{r}$ is the
total resistance force.

Table 1: Parameters of Vehicle Symbol | Description | Value
---|---|---
$m_{\text{eff}}$ | Effective mass of the vehicle |  1870 $kg$
$C_{d}$ | Coefficient of drag | 0.373
$A$ | Frontal projected area of vehicle |  2.58 $m^{2}$
$\rho$ | Air density |  1.205 $kg/m^{3}$
$f_{r}$ | Rolling resistance coefficient | 0.011
$g$ | Gravitational constant |  9.8 $m^{3}kg^{-1}s^{-2}$
$\eta$ | Gearbox efficiency | 0.94
$r_{w}$ | Wheel radius |  0.364 $m$
$I_{g}$ | Gear ratio of gearbox | 0.672
$I_{\text{final}}$ | Final driver ratio | 4.103
$T_{\text{drag}}$ | Engine drag torque |  30 $Nm$
$I_{\text{engine}}$ | Engine inertia |  0.15 $kg\,m2$
$\alpha_{1}$ | Coefficient | 0.2159
$\alpha_{2}$ | Coefficient | 0.005676
$\alpha_{3}$ | Coefficient | 0.0004349
$\alpha_{4}$ | Coefficient | 8.899E-07

The total resistance force, $F_{r}$, consists of the aerodynamic resistance,
the grade resistance, and the rolling resistance, as calculated by

$F_{r}(s,v)=m_{\text{eff}}g\left(\sin(\phi(s))+f_{r}\cos(\phi(s))\right)+\frac{C_{d}\rho
A}{2}v^{2},$ (2)

where $g$ is the gravitational constant, $\phi(s)$ is the slope data in the
distance domain, $f_{r}$ is the coefficient of rolling resistance, $C_{d}$ is
the drag coefficient, $\rho$ is the air density, and $A$ is the frontal
projected area of the vehicle. Assuming no slip on the wheel, the wheel
rotational speed $\omega_{w}$ and traction force $F_{t}$ can be calculated as:

$\displaystyle\omega_{w}$ $\displaystyle=\frac{v}{r_{w}},$ (3) $\displaystyle
F_{t}$ $\displaystyle=\frac{1}{r_{w}}T_{\text{final}},$

where $T_{\text{final }}\in\left\\{T_{\text{final },1},T_{\text{final
},2},T_{\text{final, }3}\right\\}$ is the final driveshaft torque for
baseline, fuel cut-off, and engine start/stop methods respectively, and
$r_{w}$ is the wheel radius.

The vehicle dynamic equations in (2) are defined in the time coordination.
However, the slope information represented in the digital map is location-
specific, and is more naturally correlated to spatial information. To that
end, we adopt the spatial-based model, as shown below

$\frac{\text{d}v}{\text{d}t}=\frac{\text{d}v}{\text{d}s}\frac{\text{d}s}{dt}=\frac{\text{d}v}{\text{d}s}v.$
(4)

The model in the spatial domain can be given as

$\begin{split}\frac{\text{d}v}{\text{d}s}=\frac{1}{m_{\text{eff}}v}\left[F_{t}-F_{r}(s,v)\right].\end{split}$
(5)

By applying Euler forward method and sampling in the $s$-domain, the vehicle
dynamic model is discretized as

$\begin{split}v(k+1)&=g\left(v(k),T_{e}(k),T_{b}(k)\right)\\\
&=v(k)+\frac{1}{m_{\text{eff}}v(k)}\left[F_{t}-F_{r}(kds,v(k))\right]\Delta
s,\end{split}$ (6)

where $\Delta s$ is the discretization distance in each step.

###  I-B Powertrain Model

A quasi-static engine model is considered to quantify fuel consumption. A
polynomial fuel consumption model is fitted using the measurements on an
engine dynamometer provided by the manufacturer. The idle fuel mass flow rate,
which is absent from the engine dynamometer measurements, is measured with
fuel meter during on-road tests. The following model is constructed
considering the computational complexity and fuel cut-off mechanism:

$\dot{m}_{f}=\alpha_{1}+\alpha_{2}\omega_{e}T_{e}+\alpha_{3}\omega_{e}^{2}T_{e}+\alpha_{4}\omega_{e}T_{e}^{2}$
(7)

where $\dot{m}_{f}$ $[g/s]$ is the instantaneous fuel consumption rate, and
$\alpha_{1\sim 4}$ are the constant coefficients as shown in Table 1.

Note that:

  * $\blacktriangleright$

A constant term is included in the polynomial equation to represent the idle
fuel consumption;

  * $\blacktriangleright$

A 3rd-order polynomial is used to avoid over-fitting and to reduce computation
burden;

  * $\blacktriangleright$

Engine torque impacts every term besides the constant one.

The reformulated fuel consumption in the distance domain is represented as

$M_{f}(k)=\frac{\dot{m}_{f}\left(T_{e}(k)\right)}{v(k)}\Delta s$ (8)

where $M_{f}(k)$ is the fuel consumption over the discretization step.

In this study, three different coasting methods are considered and compared.

####  I-B1 Baseline Model

When the vehicle cruises on the highway, the transmission is set at the
highest gear which corresponds to the high efficient area in the Brake
Specific Fuel Consumption (BSFC) map. The lock-up clutch keeps engaged.
Therefore, the drivetrain model is given by

$\displaystyle T_{\text{final},1}(k)$ $\displaystyle=\eta
I_{g}I_{\text{final}}T_{e}(k)-T_{b}(k),$ (9) $\displaystyle\omega_{w}(k)$
$\displaystyle=\frac{\omega_{e}(k)}{I_{g}I_{\text{final}}},$

where $T_{\text{final}}$, $T_{e}$, $T_{b}$ are the torque of final driveshaft,
engine, and brake, respectively, $\omega_{e}$ is the rotational velocity of
engine, $\eta$ is the transmission efficiency, $I_{g}$, and $I_{\text{final}}$
are the gear ratios of gearbox and final driveshaft, respectively.

####  I-B2 Fuel Cut-Off

Figure 1: Field test data in Wuhan.

The fuel cut-off (FCO) strategy, as an efficient way to reduce fuel
consumption, has been applied in production engines. As an example, Fig. 1
shows vehicle driving data collected on a test platform, FCO is enabled when
the request engine torque is zero. For FCO, the clutch is engaged and holds
the transmission in the highest gear while coasting. The negative engine
torque in Fig. 1, measured based on the experimental data, indicates the drag.

Following the analysis of field test data, the coasting dynamic with FCO
mechanism in the distance domain after the discretization is modeled as

$\displaystyle T_{\text{final},2}(k)=\eta
I_{g}I_{\text{final}}\left(T_{e}(k)-(1-z)T_{\text{drag }}\right)-T_{b}(k)$
(10)

where $z\in\\{0,1\\}$ is FCO signal, with 0 and 1 representing the activation
and inactivation of fuel cut-off respectively. $T_{\text{drag }}$ is a
constant representing the engine drag torque when fuel is cut off.

The fuel consumption model is modified as follow:

$M_{f}(k)=\left\\{\begin{array}[]{cl}\frac{\dot{m}_{f}\left(T_{e}(k)\right)}{v(k)}\Delta
s,&z(k)=1\quad(\text{fuel injected}),\\\ 0,&z(k)=0\quad(\text{fuel
off}).\end{array}\right.$ (11)

(11) reflects that there is no idle fuel consumption when fuel is cut off.

####  I-B3 Engine Start/Stop

For a conventional vehicle, another fuel efficient strategy during coasting is
to turn off the engine and disengage the clutch. Then, the engine will be
restarted leveraging the slipping of the clutch powered by vehicle inertia.
The clutch start method is superior to belt-driven starter systems or electric
motors mounted on the engine’s crankshaft in the high speed cruising condition
[13]. For the sake of simplicity, “Engine Start/Stop” refers to “Engine
Start/Stop coordinated with the engagement or disengagement of the lock-up
clutch”.The powertrain model is formulated as

$\displaystyle T_{\text{final},3}(k)=\eta
I_{g}I_{\text{final}}T_{e}(k)-T_{b}(k),$ (12)

the engine drag torque is eliminated thanks to the disengaged clutch.

Once the engine is turned off and the clutch is disengaged, the rotational
speed of the engine will reduce to zero and the rotational speed of the
transmission will be proportional to the vehicle speed due to the mechanical
connection. To achieve a smooth re-engagement process, the transmission is
assumed to stay in the highest gear to minimize the speed difference of the
clutch. High-speed cruise conditions normally call for the highest gear for
efficient engine operation area. Therefore, we assume that the gear is
constant during the process. This assumption is validated by the experiment
date in Fig. 1.

This restart process can be regarded as using the vehicle inertial energy to
“charge” the engine. The extra cost used to restart the engine is modeled
based on the energy balance

$\begin{gathered}\Delta E_{\text{engine }}=\frac{1}{2}I_{\text{engine
}}\left(\omega_{e}^{2}(k)-0\right)\\\
\omega_{e}(k)=\frac{I_{g}I_{\text{final}}v(k)}{2\pi r_{w}},\end{gathered}$
(13)

where $\Delta E_{\text{engine }}$ is the energy required for cranking the
engine from rest to the synchronization speed of transmission.
$I_{\text{engine }}$ is the engine inertia. $\omega_{e}(k)$ is the engine
speed after synchronization, which is proportional to the vehicle speed due to
the engaged lock-up clutch. In this study, it is assumed that the engine
start/stop and clutch engage/disengage process is finished within one sampling
interval. When the engine is turned on at the $k$th step. According to the
energy balance, the partial inertia energy of the vehicle will be transformed
into engine inertial energy. This transformation is formulated as

$\displaystyle\frac{1}{2}m_{\text{eff}}^{2}v(k+1)=$
$\displaystyle\frac{1}{2}m_{\text{eff}}^{2}v(k)+\left[F_{t}-F_{r}\right]\Delta
s+\delta E_{v}(k),$ (14) $\Delta E_{\text{engine }}(k)=\delta
E_{v}(k)=\frac{1}{2}m_{\text{eff}}\left(\bar{v}^{2}(k+1)-v^{2}(k)\right),$
(15)

where $\delta E_{v}$ is the partial inertia energy of the vehicle and
$\bar{v}(k+1)$ is the partial velocity variation caused by the engine start
process within this interval. Thus, the velocity reduction caused by the
engine start is formulated as

$\Delta v(k)=\left\\{\begin{array}[]{cc}\bar{v}(k+1)-v(k)&\text{ if
}d(k)=1\text{ and }d(k-1)=0\\\ 0&\text{ else }.\end{array}\right.$ (16)

where $d\in\\{0,1\\}$ is the engine start/stop signal, with 0 and 1
representing that the engine is turned off and turned on respectively.

The coasting dynamic with engine start mechanism in the distance domain after
the discretization can be expressed as

$\begin{split}v(k+1)&=g\left(v(k),T_{e}(k),T_{b}(k)\right)\\\
&=v(k)+\frac{1}{m_{\text{eff}}v(k)}\left[F_{t}-F_{r}\right]\Delta s-\Delta
v(k).\end{split}$ (17)

##  II Offline evaluation of eco-coasting methods

In this section, optimization problems regarding the three coasting methods,
including the baseline one, are summarized. Then, the performance of these
coasting strategies is evaluated using DP under different slope profiles from
the digital map. As we are mostly concerned with fuel economy and time spent
in traveling, fuel consumption and travel time are used in the cost function
and the weighted sum method is used to evaluate the trade off between these
competing objectives.

###  II-A Fuel-efficient Optimal Control of Baseline Method

Given a priori slope profile with a length of $N$, an optimization problem is
formulated to find an optimal control law, that minimizes the cost function
including fuel consumption and travel time over the entire driving cycle,

$\min_{T_{e}(.),T_{b}{(.)}}\sum_{i=0}^{N}(\beta({M}_{f}(k))+(1-\beta)\frac{1}{v(k)})\Delta
s$ (18)

subject to the system dynamic constraints (6), (9) and

$\displaystyle 0\leq T_{e}(k)\leq T_{e,\max},\quad k=0:N-1$ (19)
$\displaystyle 0\leq T_{b}(k)\leq T_{b,\max},\quad k=0:N-1$ $\displaystyle
v_{\min}\leq v(k)\leq v_{\max},\quad k=0:N$ $\displaystyle
v(0)=v_{0},v(N)=v_{0},$

where $T_{e,\max}$ are the engine torque constraints, $T_{b,\max}$ are the
brake torque constraints. $v_{\min}$ and $v_{\max}$ are the lower and upper
limits of cruise speed, implying that the vehicle maintains the highest gear
within this speed range. $\beta$ is the weight ranging from 0 to 1. The
initial and terminal speed is fixed as $v_{0}$ to make a fair comparison.

###  II-B Fuel-efficient Optimal Control of FCO Method

To avoid a high switching frequency of the FCO signal, which is associated
with driving comfort, an additional penalty is added to the cost function

$\displaystyle\min_{T_{e}(\cdot),T_{b}(\cdot),z(\cdot)}\sum_{i=0}^{N}(\beta
M_{f}(k)$ $\displaystyle+(1-\beta)\frac{1}{v(k)}$ (20)
$\displaystyle+\alpha\cdot(z(k+1)-z(k))^{2})\Delta s$

where the added term $\cdot(z(k+1)-z(k))^{2}$ reflects the FCO switch
frequency during the whole trip. With a weighted coefficient $\alpha$, this
term limits the switching numbers of the FCO signal. The cost of (20) subject
to the system dynamics (6), (10) and

$\displaystyle 0\leq T_{e}(k)\leq T_{e,\max}\cdot z(k),\quad k=0:N-1$ (21)
$\displaystyle 0\leq T_{b}(k)\leq T_{b,\max}\cdot(1-z(k)),\quad k=0:N-1$
$\displaystyle v_{\min}\leq v(k)\leq v_{\max},\quad k=0:N$ $\displaystyle
v(0)=v_{0},v(N)=v_{0}$ $\displaystyle z(k)\in\\{0,1\\},\quad k=0:N-1.$

The constraints above imply that $T_{e}(k)=0$ when FCO is activated, and
$T_{b}(k)=0$ when FCO is inactivated. This condition automatically assures
that $T_{e}$ and $T_{b}$ will not be non-zero simultaneously.

###  II-C Fuel-efficient Optimal Control of Engine Start/Stop Method

Since the cost to restart the engine is formulated in the dynamic equation
(12) and (17), the penalty term $(d(k+1)-d(k))^{2}$ intended to limit frequent
switch is redundant and the cost function for engine start/stop control is
comprised of fuel consumption and travel time as same as the baseline,

$\displaystyle\min_{T_{e}(\cdot),T_{b}(\cdot),d(\cdot)}\sum_{i=0}^{N}(\beta
M_{f}(k)$ $\displaystyle+(1-\beta)\frac{1}{v(k)})\Delta s$ (22)

subject to the system dynamics (12), (17). The same state variable and control
variable constraints as (LABEL:Constraints_of_FCO) for the fuel cut-off method
are imposed to the optimization problem formulation of the engine start/stop
method, except that the fuel cut-off signal $z(k)$ in
(LABEL:Constraints_of_FCO) is replaced by the engine start/stop signal $d(k)$.

###  II-D Dynamic Programming

To compare the fuel saving benefit among different eco-coasting strategies,
Dynamic Programming [14], as a global optimization approach, is used in this
research to derive the timing and duration of the coasting signal along with
the engine torque and brake torque.

The optimal control law is solved using backward iteration by minimizing the
cost-to-go function. In this study, the discretization distance is 5 $m$, the
other parameters corresponding to constraints are listed in Table 2. Note that
the maximum of the brake torque is set to 500 $Nm$, rather than the actual
upper limits as specified by the brake system characteristics, to avoid
unexpected sudden brake during this cruising condition.

Table 2: Parameters of optimization problem Parameter | Description | Value
---|---|---
$v_{\min}$ | Minimum of velocity | 50 $km/h$
$v_{\max}$ | Maximum of velocity | 90 $km/h$
$T_{e,\max}$ | Maximum of engine torque | 120 $Nm$
$T_{b,\max}$ | Maximum of brake torque | 500 $Nm$
$v_{0}$ | Initial and terminal velocity | 75 $m/s$

###  II-E Evaluation and Discussion

Figure 2: The performance trade off with FCO coasting method by varying
$\alpha$ (left). FCO signal with different $\alpha$ under Wuhan slope profile
(right). Figure 3: The Pareto Front with baseline, FCO, and engine start/stop
coasting methods under different slope profile: (a) and (d) The standard
uphill-downhill slope profile and the corresponding Pareto Front, (b) and (e)
The Chongqing slope profile and the corresponding Pareto Front, (c) and (f)
The Wuhan slope profile and the corresponding Pareto Front.  Figure 4: DP
simulation results under standard uphill-downhill slope profile with
$\beta=0.1$, $\beta=0.6$, and $\beta=0.9$. Figure 5: DP simulation results
under Chongqing slope profile with $\beta=0.6$, and $\beta=0.9$.

In this section, three slope profiles are used to assess the performance
metrics of different coasting strategies with DP algorithm. The first slope
profile is a typical uphill-downhill segment for 1000 m as shown in Fig. 3.
The other two profiles come from measurements on Chongqing and Wuhan highway
separately. There is no single solution that would simultaneously optimize all
criteria which are naturally in conflict. Varying the weighting factor $\beta$
allows us to put a different relative emphasis on each attribute to
investigate performance trade-offs (i.e., the Pareto front). The Pareto fronts
provide insights into the trade-off between travel time and fuel consumption
for different coasting methods. In such case, none of solutions from the
Pareto front can be improved, in the sense that improving one attribute will
lead to deterioration of the other [30].

Based on the simulation results in Fig. 2, for the FCO method, more FCO
mechanism is activated when the $\alpha$ is assigned as a larger value
consuming less fuel but extending the time.

In terms of travel time and fuel consumption as shown in Fig. 3, the baseline
strategy and FCO strategy have the same performance. The Engine start/stop
strategy is outperformed than the other two method due to the fact that the
performance of baseline strategy and FCO strategy is dominated by many
solutions from the engine start/stop strategy. The reason is that, as shown in
Fig. 4 with $\beta=0.1$, the vehicle will turn off the engine to coast in
advance considering the terminal speed constraints. The disengagement of the
lock-up clutch during the engine off period can eliminate the engine drag
torque. Therefore, the decoupled engine and transmission yield a moderate
deceleration and a longer coasting distance consequently. The longer coasting
distance, compared with the FCO approach, can decrease fuel consumption
without compromising travel time dramatically. This pattern is also observed
when the weight turns to 0.6 and 0.9, as shown in Fig. 4.

This optimization result leverages engine-off condition under downhill,
leading to better performance. Based on the simulation results under the
Chongqing slope profile with $\beta=0.6$ (Fig. 5), the engine-off is activated
in almost every downhill. In comparison, to reduce idle fuel consumption, the
fuel cut-off condition is only activated at the end of the driving cycle.
Furthermore, the terminal speed constraint is satisfied as well thanks to the
engine drag torque under the downhill condition.

Under Chongqing and Wuhan slope profiles, the FCO and engine start/stop
coasting methods can achieve better performance with increased weight $\beta$,
indicting that fuel consumption will be emphasized more in the cost function
and this will induce more fuel cut-off and engine-off conditions, as
illustrated in Fig. 5.

The cost for restarting the engine is included in the dynamic equation (17)
and will reduce the switching frequency of the engine start/stop signal. As
shown in Fig. 5, this formulation will not only provide a fair evaluation of
the overall performance (i.e., fuel consumption and travel time), but also
eliminate the driving discomfort associated with frequent start/stop.

##  III Online Implementation with MPC

Based on the DP simulation results of Section II, the engine start/stop
coasting strategy has an obvious advantage compared with the other two
methods. Compared with the baseline, FCO provides moderate improvement with
significant emphasis put on fuel consumption. However, this offline calculated
solution is not feasible for online implementation as it requires future
information of the driver’s inputs. To develop an online solution, MPC is
developed for these two proposed coasting methods and performance is
evaluated. To better assess the performance of these two coasting approaches
with different control methods (e.g., PI controller, MPC, DP), a reference
speed is extracted from the in-field experiment in which the vehicle is
controlled by a human driver [31]. In this section, the previous travel time
term in the cost function is replaced by the tracking error. The cost function
of DP for FCO and engine start/stop coasting strategies are reformulated as

$\displaystyle\min_{T_{e}(\cdot),T_{b}(\cdot),z(\cdot)}\sum_{i=0}^{N}($
$\displaystyle\beta\cdot\left(M_{f}(k)\right)+$ (23)
$\displaystyle(1-\beta)\cdot\left(v(k)-v_{\text{ref}}\right)^{2})\Delta s,$

and

$\displaystyle\min_{T_{e}(\cdot),T_{b}(\cdot),z(\cdot)}\sum_{i=0}^{N}(\beta\cdot
M_{f}(k)$ $\displaystyle+(1-\beta)\cdot\left(v(k)-v_{\text{ref}}\right)^{2}$
(24) $\displaystyle+\alpha\cdot(z(k+1)-z(k))^{2})\Delta s.$

In this section, a rule-based approach is first designed to provide a
reference. Then, a speed-tracking MPC formulation is proposed for these two
coasting strategies. A novel constraint is applied to satisfy the minimum off
steps constraints of the MIMPC and mitigate the switching frequency. After
that, the performance is discussed by comparing with the PI controller and DP
simulation results. Moreover, the effect of the length of the prediction
horizon is demonstrated. Finally, to mitigate the computation burden, a
heuristic MPC is proposed by using the engine start/stop sequence obtained
from the offline DP results.

###  III-A PI Controller

FCO is available in our in-field experiment platform in this study [31]. It is
activated when no torque is requested and the speed is higher than a
threshold. Inspired by this trigger condition, a rule-based FCO or Engine-
turn-off mechanism is proposed. The FCO or Engine turn-off signal is activated
if the current speed is higher than $75km/h$ and the required engine torque is
decreased to zero simultaneously. A PI controller is applied to track the
reference speed, as shown in Fig. 6.

Figure 6: The block diagram of PI controller. $T_{e}^{req}$ is the required
engine torque, $T_{b}^{req}$ is the required brake torque, $z(k-1)$ is the
binary variable at the last step.

The FCO or engine turn-off signal will hold for 4 discretization steps if it
is activated to avoid frequent switches, similar to the engine on command
filter in [32]. The operational constraints on minimum off discretization
steps chosen in this study is tuned offline considering the trade-off between
tracking accuracy and fuel consumption.

###  III-B MPC

MPC, also known as the receding horizon control, is a control technique that
embeds optimization within feedback to deal with systems subject to
constraints on inputs and states [33]. MPC determines the control action by
solving a finite horizon open-loop optimal control problem online at each
discrete instant. However, in this study, there exists discrete dynamics in
the engine start/stop or FCO coasting methods, resulting in a mixed-integer
optimal control problem (MIOCP).

Frequent switch of the engine start/stop or fuel cut-off signal will affect
driving comfort and cause carbon deposits in the engine. To avoid that, the
operational constraints on the minimum amount of steps for which an engine or
fuel injection must be kept off can be expressed by the following mixed
integer linear inequalities:

$\delta(\tau)\leq 1-(\delta(j-1)-\delta(j))$ (25)

with $\tau=j+1,\ldots,\min\left(j+d_{\min}-1,N\right)$, $d_{\min}=4$
represents the operational constraints on minimum off steps,
$j=2,\ldots,N+d_{\min}$. Let the extended state be represented by the vector
$\delta:=\left[z_{\text{history }},z_{\text{current }}\right]^{T}$, where
$z_{\text{history }}$ is the vector of the previous binary variables with
dimension $d_{\min}\times 1$ as

$\displaystyle z_{\text{history
}}:=\left[z_{(k-d_{\min}-1)},\ldots,z_{k-1}\right]^{T}.$ (26)

$z_{\text{current}}$ is the row vector of the current binary variables with
dimension $N_{h}\times 1$ as

$\displaystyle
z_{\text{current}}:=\left[z_{1|k},\ldots,z_{N_{h|k}}\right]^{T}.$ (27)

The implementation of the mixed integer linear inequalities in the MPC
formulation is shown in Algorithm 1. More details about the constraints can be
found in [19] and [34].

Algorithm 1 MIMPC with minimum off constraints

1:Initialize $v_{0|k}$, $C$ and $z_{\text{history }}$;

2:for $j=1,\ldots,N_{h}$ do

3: for $i=1,\ldots,(d_{\dim}+N_{h})$ do

4: $C:=[C,\tau(k:\min(N_{h},i+d_{\dim}-1))\leq\tau(i)-\tau(i-1)]^{T}$

5: end for

6: Impose constraints vector $C$ to original MIOCP

7: Calculate $T_{e}\in\mathcal{R}^{N_{h}}$, $T_{b}\in\mathcal{R}^{N_{h}}$,
$z\in\mathcal{R}^{N_{h}}$

8: Update
$z_{\text{history}}:=\left[z_{j-d_{\min}},\ldots,z_{k-1},z_{1|k}\right]^{T}$

9:end for

####  III-B1 FCO

Let $x(i|k)$ denotes the value of variable $x$ at distance $k+i$ while the
prediction is made at the current location $k$. The MIMPC for fuel cut-off is
defined over a finite-distance horizon ($N_{h}$) as:

$\displaystyle\min_{T_{e}(\cdot|k),T_{b}(\cdot|k),z(\cdot|k)}\sum_{i=0}^{N_{h}}($
$\displaystyle\beta\cdot\left(M_{f}(i|k)\right)+$ (28)
$\displaystyle(1-\beta)\cdot\left(v(i|k)-v_{\text{ref}}\right)^{2})\Delta s$

subject to:

$\displaystyle
v(i+1|k)=v(i|k)+\frac{1}{m_{\text{eff}}v(i|k)}\left[F_{t}-F_{r}\right]\Delta
s,i=0:N_{h}-1$ (29) $\displaystyle 50km/h\leq v(i|k)\leq 90km/h,i=0:N_{h}$
$\displaystyle 0Nm\leq T_{e}(i|k)\leq 120Nm\cdot z(i|k),i=0:N_{h}-1$
$\displaystyle 0Nm\leq T_{b}(i|k)\leq 500Nm\cdot(1-z(i|k)),i=0:N_{h}-1$
$\displaystyle z(i|k)\in\\{0,1\\},i=0:N_{h}-1$ $\displaystyle\delta(\tau)\leq
1-(\delta(j-1)-\delta(j)),j=2,\ldots,N_{h}+d_{\min}$ $\displaystyle
v(0)=75km/h.$

The optimal coasting action is achieved by solving the above minimization
problem repeatedly in a receding-horizon manner.

####  III-B2 Engine Start/Stop

The MIMPC for engine start/stop method is defined over a finite-distance
horizon ($N_{h}$) as:

$\displaystyle\min_{T_{e}(\cdot|k),T_{b}(\cdot|k),z(\cdot|k)}\sum_{i=0}^{N_{h}}($
$\displaystyle\beta\cdot\left(M_{f}(i|k)\right)+$ (30)
$\displaystyle(1-\beta)\cdot\left(v(i|k)-v_{\text{ref}}\right)^{2})\Delta s$

subject to the system dynamics:

$\displaystyle v(i+1|k)=v(i|k)$ (31)
$\displaystyle+\frac{1}{m_{\text{eff}}v(i|k)}\left[F_{t}-F_{r}\right]\Delta
s-\Delta v(i|k),i=0:N_{h}-1.$

Moreover, the same state variable and control variable constraints as
(LABEL:constraints_of_FCO) for the FCO method with the MPC framework are
imposed to the optimization problem formulation of the engine start/stop
method, except that the fuel cut-off signal $z(k)$ in
(LABEL:constraints_of_FCO) is replaced by the engine start/stop signal $d(k)$.

We implement the proposed MIMPC above in YALMIP [35] and solve the
corresponding non-convex non-linear optimal programming using the BARON
commercial solver [36]. The maximum number of allowable iteration is set as
200 since the difference of lower bound and upper bound can converge to zero
in most steps for this Branch $\&$ Bound based solver. If the solver can not
find a feasible solution within the max iteration time, the last feasible
solution will be applied to the plant to guarantee the recursive feasibility.

###  III-C Performance Evaluation and Discussion

Figure 7: Performance comparison with FCO method: DP
($\alpha=0.004,\beta=0.7$) vs MPC ($\beta=0.16$)vs PI controller
($P=0.5,I=0.0001$). Figure 8: Simulation results comparison with FCO method
under the Wuhan slope profile.

For the FCO coasting method under the Wuhan slope profile, compared with the
PI controller, DP achieves a 6.66% decrease in fuel consumption while the
travel time is decreased by $0.33\%$ (Fig. 7). The tracking error of the PI
controller is higher than DP considering the operational constraints on the
minimum off steps. There are two reasons for the fuel consumption reduction:
first, the timing and duration of the fuel cut-off results of the rule-based
PI controller show a significant difference compared with the global optimal
solution from the DP results. Second, more brake torque is involved in the PI
controller as shown in Fig. 8.

The proposed MPC achieves a substantial performance improvement with respect
to the fuel consumption with the FCO coasting method. There is a 6.15%
reduction in fuel consumption compared with the PI controller without
compromising the travel time, almost equal to the DP result. As shown in Fig.
8, the MPC repeats almost every fuel cut-off action of the global optimal
solution but will also take extra fuel cut-off action due to its finite
horizon.

Figure 9: Performance comparison with engine start/stop approach: DP
($\beta=0.7$) vs MPC ($\beta=0.16$) vs PI controller ($P=2,I=0.001$). Figure
10: Performance comparison with engine start/stop approach: DP
($\beta=0.9985$) vs MPC ($\beta=0.3$) vs PI controller ($P=0.18,I=0.01$).

For the engine start/stop coasting method, the fuel consumption of MPC is less
than the results of the PI controller, close to the DP results. The distinct
engine off decisions and duration among these three methods lead to the
different performance as shown in Fig. 9. As shown in Fig. 11, the engine is
turned off by the DP at 4500 m as the reference speed starts to decrease at
the downhill condition. The brake torque is involved to track the reference
speed precisely. MPC also predicts this engine-off action exactly, but this
engine start/stop signal switches multiple times when the prediction horizon
is chosen as 200 m. Regarding the results of the PI controller at the same
spot, as shown in Fig. 11, the engine is turned off for a while when its speed
is higher than 75 km/h. After that, the engine is turned on referring to its
rule-based trigger condition and fuel injection continues to maintain engine
idle speed.

Compared with the MPC performance between FCO (Fig. 7) and engine start/stop
(Fig. 9) coasting methods, the FCO has a better performance in terms of the
fuel consumption reduction compared with the results of PI controller.
However, if the tracking error of the engine start/stop method is relaxed to
the same level as the FCO method (Fig. 10), there will be a significant
improvement with respect to the fuel consumption for MPC. This observation
implies that the potential benefit from different coasting mechanisms will be
eliminated without proper speed tracking. As shown in Fig. 11, DP with
$\beta=0.9985$ calculates a longer engine-off period at the cost of a larger
tracking error compared with the DP results with $\beta=0.7$ in Fig. 12 at the
same location. More engine-off action is calculated to minimize the fuel
consumption while sacrificing the tracking accuracy a little bit.

Figure 11: Simulation results comparison with engine start/stop method zoomed
from 4000 $m$ to 5500 $m$ of the Wuhan slope profile.

###  III-D Sensitivity Analysis on the Length of the Prediction Horizon

Figure 12: Effect of increasing the prediction horizon with engine start/stop
method. (a) The fuel consumption; (b) The travel time; (c) The tracking error;
(d) The computation time for obtaining our numerical solution at each instant
Figure 13: Simulation comparison with engine start/stop method over different
length of prediction horizon Figure 14: Simulation comparison with engine
start/stop method: PI vs MPC vs Heuristic MPC vs DP. (a) The fuel consumption;
(b) The travel time; (c) The tracking error; (d) The computation time for
obtaining our numerical solution at each instant  Figure 15: Simulation
comparison with engine start/stop method zoomed from 12600 $m$ to 13600 $m$ of
the Wuhan slope profile: MPC vs Heuristic MPC vs DP  Figure 16: Robustness
analysis with engine start/stop method: MPC vs Heuristic MPC Figure 17:
Robustness analysis with engine start/stop method: MPC vs Heuristic MPC. (a)
The fuel consumption; (b) The travel time; (c) The tracking error; (d) The
computation time for obtaining our numerical solution at each instant

To demonstrate the effect of the prediction horizon on MPC performance, the
engine start/stop method is used as an example. The computations are performed
on the MATLAB R2019a platform running on an AMD Ryzen 5 2600X 3.60-GHz with
16.0-GB RAM. Fig. 12 shows that by increasing the prediction horizon from 50 m
to 100 m, the tracking error will decrease from $31.89\%$ to $8.21\%$ and the
fuel consumption will decrease by $1.21\%$. The reason is that increasing the
prediction horizon makes more reference speed information available for the
controller. MPC controller with a longer prediction horizon can accelerate
earlier to track the reference speed precisely and coast earlier on the
downhill to leverage the kinetic energy (Fig. 13). However, extending the
prediction horizon increases the computation footprint significantly as shown
in Fig. 12 (d). It also makes the MPC susceptible to speed prediction error.
The fuel consumption is increased by extending the horizon since the tracking
error is emphasized in the cost function. The controller will calculate more
engine torque for a better tracking performance with a longer prediction
horizon (Fig. 13).

###  III-E Heuristic MPC

Solving such an MINLP problem as shown in (30) and (31) is computationally
demanding. To eliminate the computational burden for real-time implementation,
the offline DP result of the engine start/stop sequence is utilized in the MPC
considering the fact that the engine start/stop decision is settled to some
extent for given road grade and reference speed [37]. Then the optimal problem
is tailored as

$\displaystyle\min_{T_{e}(\cdot|k),T_{b}(\cdot|k)}\sum_{i=0}^{N_{h}}($
$\displaystyle\beta\cdot\left(M_{f}(i|k)\right)+$ (32)
$\displaystyle(1-\beta)\cdot\left(v(i|k)-v_{\text{ref}}\right)^{2})\Delta s$

subject to:

$\displaystyle
v(i+1|k)=v(i|k)+\frac{1}{m_{\text{eff}}v(i|k)}\left[F_{t}-F_{r}\right]\Delta
s,i=0:N-1$ (33) $\displaystyle 50km/h\leq v(i|k)\leq 90km/h,i=0:N$
$\displaystyle 0Nm\leq T_{e}(i|k)\leq 120Nm,i=0:N-1$ $\displaystyle 0Nm\leq
T_{b}(i|k)\leq 500Nm,i=0:N-1$ $\displaystyle v(0)=75km/h.$

The integer variable $d$ and the operational constraints $\delta(\tau)\leq
1-(\delta(j-1)-\delta(j))$, which avoids the frequent switch, are dropped from
the original problem definition in (30) and (31). The simulation result shows
that the computation time will reduce dramatically. As shown in Fig. 14(d),
the computation time at each instance of the heuristic MPC is reduced to 7.0
s. Note that there is still a lot of work to do to further reduce the
computation times, for instance, by replacing the commercial NLP solver with
more efficient or tailored NLP solvers [38], exploiting real-time iteration
strategy [39], converting the MATLAB code to C code. These are left for future
research.

The fuel consumption of this heuristic MPC is better than the MPC result due
to the use of the offline engine start/stop results calculated with DP, and
only slightly worse than the DP results. As shown in Fig. 15, the heuristic
MPC will not even accelerate when the engine is turned on till the reference
speed increased dramatically over the prediction horizon.

However, the offline calculated engine start/stop sequence depends to a larger
extent on the road grade and driver’s demand (i.e., reference speed). When the
reference speed is changed from 11240 m to 14600 m as shown in Fig. 16, the
MPC can track the reference speed accurately. But the heuristic MPC, which
uses the pre-calculated engine start/stop sequence derived from the original
reference speed, shows a trajectory deviation around 12900 m due to the
inappropriate engine off signal from the results of DP. The tracking error of
heuristic MPC is increased by 46.72 % since the engine start/stop sequence is
not updated with the modified speed reference, while the fuel consumption and
travel time are almost the same (Fig. 17).

##  IV Conclusion

In this paper, two eco-coasting strategies are proposed using road grade
preview information. The evaluation and online implementation are investigated
with DP and MPC respectively. The dynamic equation considering the engine
restart cost will not only give a more precise evaluation but also avoid the
frequent switch of the engine start/stop signal. Finally, the online
performance is evaluated by reformulating the optimal control problems of eco-
coasting strategies as MIMPC problems. Simulation results show that the MIMPC
improves fuel consumption to a level comparable to DP without sacrificing the
travel time for both fuel cut-off and engine start/stop methods. The heuristic
MPC can reduce the computation footprint by utilizing the pre-calculated
engine start/stop sequence derived from DP results when the road grade and
driver’s demand (e.g., reference speed) are known in advance. As for the
implementation of the eco-coasting strategy in the future, the V2X information
can help to identify the reference speed, the pre-calculated coasting signal
can be derived with cloud computing before departure.

Since the mixed integer nonlinear programming in each discrete distance
instance requires high computational cost, an algorithm for solving the MIMPC
in real-time is a subject left for future research.

## Acknowledgment

Y. Yan acknowledges the support from the Chinese Scholarship Council that made
this work possible. This work was done when he was a visiting Ph.D. student in
the Department of Naval Architecture and Marine Engineering, University of
Michigan, Ann Arbor, MI 48109 USA. The work of B. Z. Gao and H. Chen was
supported by the China Automobile Industry Innovation and Development Joint
Fund under Grant No.U1864206, and in part by Jilin provincial Science $\&$
Technology Department No.20200301011RQ.

## References

  * [1] R. Varnhagen, “Electronic horizon: A map as a sensor and predictive control,” SAE Technical Paper, Tech. Rep., 2017.
  * [2] C. F. Salgueiredo, O. Orfila, G. Saint Pierre, P. Doublet, S. Glaser, S. Doncieux, and V. Billat, “Experimental testing and simulations of speed variations impact on fuel consumption of conventional gasoline passenger cars,” _Transportation Research Part D: Transport and Environment_ , vol. 57, pp. 336–349, 2017.
  * [3] S. E. Li, K. Deng, Y. Zheng, and H. Peng, “Effect of pulse-and-glide strategy on traffic flow for a platoon of mixed automated and manually driven vehicles,” _Computer-Aided Civil and Infrastructure Engineering_ , vol. 30, no. 11, pp. 892–905, 2015.
  * [4] A. Balluchi, M. Di Benedetto, C. Pinello, C. Rossi, and A. Sangiovanni-Vincentelli, “Cut-off in engine control: a hybrid system approach,” in _Proceedings of the 36th IEEE Conference on Decision and Control_ , vol. 5. IEEE, 1997, pp. 4720–4725.
  * [5] Y. Katsumata, S. Segawa, K. Adachi, A. Higashimata, and Y. Ochi, “Development of a slip speed control system for a lock-up clutch (part ii),” SAE Technical Paper, Tech. Rep., 2008.
  * [6] H. Koch-Groeber and J. Wang, “Criteria for coasting on highways for passenger cars,” SAE Technical Paper, Tech. Rep., 2014.
  * [7] J. Lee, “Vehicle inertia impact on fuel consumption of conventional and hybrid electric vehicles using acceleration and coast driving strategy,” Ph.D. dissertation, Virginia Tech, 2009.
  * [8] J. Bishop, A. Nedungadi, G. Ostrowski, B. Surampudi, P. Armiroli, and E. Taspinar, “An engine start/stop system for improved fuel economy,” SAE technical paper, Tech. Rep., 2007.
  * [9] P. Griefnow, J. Andert, and D. Jolovic, “Next-generation low-voltage power nets impacts of advanced stop/start and sailing functionalities,” _SAE International Journal of Fuels and Lubricants_ , vol. 10, no. 2, pp. 556–573, 2017\.
  * [10] P. Griefnow, J. Andert, F. Xia, S. Klein, P. Stoffel, M. Engels, and D. Jolovic, “Real-time modeling of a 48v p0 mild hybrid vehicle with electric compressor for model predictive control,” _SAE Tech. Paper_ , 2019\.
  * [11] C. Sohn, J. Andert, and D. Jolovic, “An analysis of the tradeoff between fuel consumption and ride comfort for the pulse and glide driving strategy,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 7, pp. 7223–7233, 2020.
  * [12] S. Choi, K. Ko, and I. Jeung, “Optimal fuel-cut driving method for better fuel economy,” _International Journal of Automotive Technology_ , vol. 14, no. 2, pp. 183–187, 2013.
  * [13] N. Mueller, S. Strauss, S. Tumback, G.-C. Goh, and A. Christ, “Next generation engine start/stop systems: “free-wheeling”,” _SAE International Journal of Engines_ , vol. 4, no. 1, pp. 874–887, 2011.
  * [14] E. Hellström, J. Åslund, and L. Nielsen, “Design of an efficient algorithm for fuel-optimal look-ahead control,” _Control Engineering Practice_ , vol. 18, no. 11, pp. 1318–1327, 2010.
  * [15] D. V. Ngo, “Gear shift strategies for automotive transmissions,” _Ph.D. Thesis_ , 2012.
  * [16] G. Li and D. Goerges, “Hybrid modeling and predictive control of the power split and gear shift in hybrid electric vehicles,” in _2017 IEEE Vehicle Power and Propulsion Conference (VPPC)_. IEEE, 2017, pp. 1–6.
  * [17] F. Yan, J. Wang, and K. Huang, “Hybrid electric vehicle model predictive control torque-split strategy incorporating engine transient characteristics,” _IEEE transactions on vehicular technology_ , vol. 61, no. 6, pp. 2458–2467, 2012.
  * [18] S. Esmaeili, S. Jadid, A. Anvari-Moghaddam, and J. M. Guerrero, “Optimal operational scheduling of smart microgrids considering hourly reconfiguration,” in _2018 IEEE 4th Southern Power Electronics Conference (SPEC)_. IEEE, 2018, pp. 1–6.
  * [19] A. Parisio, E. Rikos, and L. Glielmo, “Stochastic model predictive control for economic/environmental operation management of microgrids: An experimental case study,” _Journal of Process Control_ , vol. 43, pp. 24–37, 2016.
  * [20] B. Mayer, M. Killian, and M. Kozek, “Management of hybrid energy supply systems in buildings using mixed-integer model predictive control,” _Energy conversion and management_ , vol. 98, pp. 470–483, 2015.
  * [21] P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, and A. Mahajan, “Mixed-integer nonlinear optimization,” _Acta Numerica_ , vol. 22, pp. 1–131, 2013.
  * [22] R. Takapoui, N. Moehle, S. Boyd, and A. Bemporad, “A simple effective heuristic for embedded mixed-integer quadratic programming,” _International journal of control_ , vol. 93, no. 1, pp. 2–12, 2020.
  * [23] H. Yu, F. Zhang, J. Xi, and D. Cao, “Mixed-integer optimal design and energy management of hybrid electric vehicles with automated manual transmissions,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 11, pp. 12 705–12 715, 2020.
  * [24] J. Deng, F. Meier, and L. del Re, “Flexible predictive hybrid powertrain management with v2x information,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 3500–3505.
  * [25] P. Elbert, T. Nüesch, A. Ritter, N. Murgovski, and L. Guzzella, “Engine on/off control for the energy management of a serial hybrid electric bus via convex optimization,” _IEEE Transactions on Vehicular Technology_ , vol. 63, no. 8, pp. 3549–3559, 2014.
  * [26] S. Hadj-Said, G. Colin, A. Ketfi-Cherif, and Y. Chamaillard, “Convex optimization for energy management of parallel hybrid electric vehicles,” _IFAC-PapersOnLine_ , vol. 49, no. 11, pp. 271–276, 2016.
  * [27] H. Zhang, Q. Fan, S. Liu, S. E. Li, J. Huang, and Z. Wang, “Hierarchical energy management strategy for plug-in hybrid electric powertrain integrated with dual-mode combustion engine,” _Applied Energy_ , vol. 304, p. 117869, 2021.
  * [28] J. Han, A. Vahidi, and A. Sciarretta, “Fundamentals of energy efficient driving for combustion engine and electric vehicles: An optimal control perspective,” _Automatica_ , vol. 103, pp. 558–572, 2019.
  * [29] Y. Jia and D. Görges, “Energy-optimal adaptive cruise control based on hybrid model predictive control with mixed-integer quadratic programming,” in _2020 European Control Conference (ECC)_. IEEE, 2020, pp. 686–692.
  * [30] U. Sadek, A. Sarjaš, A. Chowdhury, and R. Svečko, “Fpga-based optimal robust minimal-order controller structure of a dc–dc converter with pareto front solution,” _Control Engineering Practice_ , vol. 55, pp. 149–161, 2016.
  * [31] H. Chu, L. Guo, B. Gao, H. Chen, N. Bian, and J. Zhou, “Predictive cruise control using high-definition map and real vehicle implementation,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 12, pp. 11 377–11 389, 2018.
  * [32] C. Hou, M. Ouyang, L. Xu, and H. Wang, “Approximate pontryagin’s minimum principle applied to the energy management of plug-in hybrid electric vehicles,” _Applied Energy_ , vol. 115, pp. 174–189, 2014.
  * [33] Z. Li, J. Sun, and S. Oh, “Path following for marine surface vessels with rudder and roll constraints: An mpc approach,” in _2009 American control conference_. IEEE, 2009, pp. 3611–3616.
  * [34] M. Carrión and J. M. Arroyo, “A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem,” _IEEE Transactions on power systems_ , vol. 21, no. 3, pp. 1371–1378, 2006.
  * [35] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in _2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508)_. IEEE, 2004, pp. 284–289.
  * [36] M. Tawarmalani and N. V. Sahinidis, “A polyhedral branch-and-cut approach to global optimization,” _Mathematical programming_ , vol. 103, no. 2, pp. 225–249, 2005.
  * [37] R. Biasini, S. Onori, and G. Rizzoni, “A near-optimal rule-based energy management strategy for medium duty hybrid truck,” _International Journal of Powertrains_ , vol. 2, no. 2-3, pp. 232–261, 2013.
  * [38] D. Liao-McPherson and I. Kolmanovsky, “Fbstab: A proximally stabilized semismooth algorithm for convex quadratic programming,” _Automatica_ , vol. 113, p. 108801, 2020.
  * [39] D. Henriksson and J. Åkesson, “Flexible implementation of model predictive control using sub-optimal solutions,” _Department of Automatic Control, Lund Institute of Technology (LTH)_ , 2004.