# Correlating thermodynamics, morphology, mechanics and thermal transport in
PMMA-PLA blends

Debashish Mukherji<EMAIL_ADDRESS>Quantum Matter Institute,
University of British Columbia, Vancouver BC V6T 1Z4, Canada  Tiago Espinosa
de Oliveira  Departamento de Farmacociências, Universidade Federal de Ciências
da Saúde de Porto Alegre, Porto Alegre 90050-170, Brazil  Céline Ruscher
Department of Mechanical Engineering, University of British Columbia,
Vancouver BC V6T 1Z4, Canada  Jörg Rottler  Quantum Matter Institute,
University of British Columbia, Vancouver BC V6T 1Z4, Canada  Department of
Physics and Astronomy, University of British Columbia, Vancouver BC V6T 1Z1,
Canada

###### Abstract

Thermodynamics controls structure, function, stability and morphology of
polymer blends. However, obtaining the precise information about their mixing
thermodynamics is a challenging task, especially when dealing with complex
macromolecules. This is partially because of a delicate balance between the
local concentration/composition fluctuations and the monomer level (multi-
body) interactions. In this context, the Kirkwood-Buff (KB) theory serves as a
useful tool that connects the local pairwise fluid structure to the mixing
thermodynamics. Using larger scale molecular dynamics simulations, within the
framework of KB theory, we investigate a set of technologically relevant
poly(methyl methacrylate)-poly(lactic acid) (PMMA-PLA) blends with the aim to
elucidate the underlying microscopic picture of their phase behavior.
Consistent with these experiments, we emphasize the importance of properly
accounting for the entropic contribution, to the mixing Gibbs free energy
change $\Delta{\mathcal{G}}_{\rm mix}$, that controls the phase morphology. We
further show how the relative microscopic interaction details and the
molecular level structures between different mixing species can control the
non-linear mechanics and ductility. As a direct consequence, we provide a
correlation that links thermodynamics, phase behavior, mechanics, and thus
also thermal transport in polymer blends. Therefore, this study provides a
guiding principle for the design of light weight functional materials with
extraordinary physical properties.

##  I Introduction

Polymers are an important class of soft materials that are extensively used
for the advanced applications [1, 2, 3], this reaches from common household
items to the advanced technologies for their use in mechanical [4, 5, 6, 7],
thermal [8, 9, 10] and optoelectronics devices [11, 12]. These systems consist
of, but are not limited to, linear polymers, copolymers, polymers with
different architectures and/or polymer blends. In this context, the latter is
of particular interest because a polymer blend provides an additional
flexibility to design functional materials with tunable properties [8, 10,
13]. Here, however, the microscopic phase behavior is important for the
usefulness of these materials for their desired applications [15, 14]. For
example, the micro-phase separation in these systems severely impacts the
materials performance by creating the weak spots that act as the nucleating
sites for the crack propagation and mechanical instabilities in these
materials [15]. Therefore, it is important to achieve a compatible chemical
composition that is not only simple from the processing point of view, but
also provides a reasonably stable homogeneous blend. Here, the standard
commodity polymers, such as poly(methyl methacrylate) (PMMA),
poly($N-$acryloyl piperidine) (PAP), poly(lactic acid) (PLA), poly(acrylic
acid) (PAA), and poly(acrylamide) (PAM), have been shown to provide
extraordinary phase, mechanical, and thermal behavior [8, 9, 10, 16, 17, 18,
13, 15, 14].

The miscibility of polymer blends at different mixing ratios is dictated by
the change in Gibbs free energy of mixing $\Delta{\mathcal{G}}_{\rm
mix}=\Delta H_{\rm mix}-T\Delta S_{\rm mix}$, where $\Delta H_{\rm mix}$ and
$T\Delta S_{\rm mix}$ are the enthalpic and the entropic contributions,
respectively. A positive value of $\Delta{\mathcal{G}}_{\rm mix}$ does not
favor miscibility, while miscibility is favored for $\Delta{\mathcal{G}}_{\rm
mix}<0$. Within the simple Flory-Huggins (FH) like representation [19, 20],
one can estimate the entropy change $\Delta S_{\rm mix}$ as,

$\Delta S_{\rm mix}=-k_{\rm B}\left[\frac{x_{\rm i}}{N_{\rm i}}\ln{x_{\rm
i}}+\frac{1-x_{\rm i}}{N_{\rm j}}\ln\left({1-x_{\rm i}}\right)\right].$ (1)

Here, $k_{\rm B}$ is the Boltzmann constant, ${x_{\rm i}}$ is mole fraction of
species i, and $N_{\rm i}$ is the degree of polymerization of species i. While
$\Delta S_{\rm mix}$ scales as $1/N$, it is always positive and thus favors
miscibility. On the contrary, when $\Delta H_{\rm mix}=\chi{x_{\rm
i}}\left({1-x_{\rm i}}\right)>>0$ (or $\chi>>0$), a system can phase separate,
with $\chi$ being the FH interaction parameter.

Following the discussion in the preceding paragraph, it becomes apparent that
a delicate balance between $\Delta H_{\rm mix}$ and $T\Delta S_{\rm mix}$
plays a key role in dictating the phase behavior of polymer blends. Therefore,
it is crucial to obtain these individual contributions to
$\Delta{\mathcal{G}}_{\rm mix}$, which is a nontrivial task in molecular
simulations. In this context, most computational studies usually deal with
enthalpy ($\chi-$based arguments), while ignoring $T\Delta S_{\rm mix}$
contributions. Here, we note in passing that most molecular simulations of
polymers commonly deal with the coarse-grained (CG) potentials that are by
construction free energies (not internal energies) and thus contain both the
entropic and the enthalpic contributions. Therefore, special attention should
be devoted when discussing results using different CG models.

Studying the thermodynamic properties of polymer mixtures is a nontrivial task
within molecular simulation protocol. Here, the fluctuation theory of
Kirkwood-Buff (KB) serves as an important method, where local fluctuations and
monomer-level interactions play a key role in dictating the mixing
thermodynamics. While the KB theory has been extensively used in studying
various solvent mixtures [21, 22, 23, 24, 25], polymer solutions [26, 27, 28],
and/or ionic systems [29], their application for the systems with long chain
blends are rather limited, if existing at all. This is in particular because
in polymers, chain connectivity introduces additional constraints that
significantly modifies the solvation structure, which is otherwise spherically
symmetric, and thus can significantly impact their solvation thermodynamics.

In this work, we investigate phase behavior and its links to the nonlinear
mechanics and the thermal transport in a set of PMMA-PLA blends at different
PLA mixing ratios $x_{\rm PLA}$. The specific system choice is motivated by
the recent experimental study [13], where it has been shown that the PMMA-PLA
blends exhibit extraordinary mechanical response and ultra-high toughness. We
individually identify the entropic and the enthalpic contributions to
highlight the thermodynamic origin of their phase behavior. We use the phase
behavior to establish a delicate link with the mechanics of polymer blends and
ultimately to the thermal transport.

The remainder of the paper is as follows: In Sec. II we sketch brief
description of the KB theory, model and method related details and the
materials used in our study. Results are presented in Sec. III and finally we
draw our conclusions in Sec. IV.

##  II Method, model and materials

###  II.1 Kirkwood-Buff theory of solution

KB theory is a well known method to investigate the thermodynamics properties
of complex mixtures. Here we only sketch the key ingredients relevant for our
work. When two different species i and j are blended in at different $x_{\rm
j}$, the fluctuation theory of Kirkwood and Buff connects fluctuations in the
grand-canonical ensemble to the pair distribution functions ${\rm g}_{\rm
ij}(r)$ via “so called” Kirkwood-Buff integral (KBI) [21, 22],

$G_{\rm ij}=4\pi\int_{0}^{\infty}\left[{\rm g}_{\rm ij}(r)-1\right]r^{2}dr.$
(2)

Here $G_{\rm ij}$ is known as the excess coordination and in polymer physics
$-G_{\rm ij}/2$ gives an estimate of the monomer excluded volume
${\mathcal{V}}_{\rm ij}$.

Following Eq. 2, the second derivative of $\Delta{\mathcal{G}}_{\rm mix}$ can
be written in terms of $G_{\rm ij}$ as [24],

$\frac{\partial^{2}\Delta{\mathcal{G}}_{\rm mix}}{\partial x_{\rm
j}^{2}}=\frac{k_{\rm B}T\left(\rho_{\rm i}+\rho_{\rm j}\right)}{x_{\rm
j}(1-x_{\rm j})\left(\rho_{\rm j}+\rho_{\rm i}+\rho_{\rm j}\rho_{\rm
i}\eta\right)},$ (3)

with $\rho_{\rm i}$ being the monomer number density of species i. The
preferential solvation parameter $\eta=G_{\rm ii}+G_{\rm jj}-2G_{\rm ij}$
gives further information about the relative (local) coordination between
different mixture components. For example, immiscibility is preferred when
$G_{\rm ii}+G_{\rm jj}\gg 2G_{\rm ij}$, weakly immiscible to miscible when
$G_{\rm ii}+G_{\rm jj}\geq 2G_{\rm ij}$ the system is miscible for $\eta\leq
0$. Note that $G$ values depend on the volume of the molecules.

Eq. 3 satisfies the boundary condition $\Delta{\mathcal{G}}_{\rm mix}(x_{\rm
j}=0.0)=0$ and $\Delta{\mathcal{G}}_{\rm mix}(x_{\rm j}=1.0)=0$. From
$\Delta{\mathcal{G}}_{\rm mix}$, one can estimate $T\Delta S_{\rm mix}=\Delta
H_{\rm mix}-\Delta{\mathcal{G}}_{\rm mix}$ with,

$\Delta H_{\rm mix}=H_{\rm mixture}-x_{\rm j}H_{\rm j}^{\rm
pure}-\left(1-x_{\rm j}\right)H_{\rm i}^{\rm pure}.$ (4)

Here, the enthalpy per monomer is calculated using $H=U+pV$, $U$ is the
internal energy including the kinetic energy, $p$ is the externally imposed
pressure and $V$ is the system volume.

###  II.2 All atom simulations

For the all-atom simulations, we have chosen two different polymers: namely
PMMA and PLA. The specific chemical structures are shown in Fig. 1.

Figure 1: Schematics showing the chemical structures of the monomer units
investigated in this work: namely poly(methyl methacrylate) (PMMA) and
poly(lactic acid) (PLA).

The explicit atom OPLS force field is used to model PLA[30], while we have
used a set of modified OPLS parameters for PMMA [31]. Earlier it has been
shown that the PMMA force field parameters reasonably capture the properties
of the solid, solvent-free PMMA, such as the specific volume, the macroscopic
elastic moduli, and the thermal transport coefficient [32]. In this work, we
have also validated that the PLA parameters give the correct estimates of
specific volume at temperature of $T=300$ K (i.e., $8.06\times 10^{-4}~{}{\rm
m^{3}/kg}$ for experiments and $8.11\times 10^{-4}~{}{\rm m^{3}/kg}$ for
simulations), elastic modulus, and thermal properties. The latter two will be
discussed at a later stage in this draft. The GROMACS molecular dynamics
package [33] is used for the simulations.

The simulations in this work are performed in four steps with different chain
lengths $N_{\ell}$ and number of chains $N_{\rm c}$. In the following we give
the details:

  * •

Set 1, the KB analysis of blends: For this step we have chosen $N_{\ell}=3$,
$N_{\rm c}=1000$ and six configurations within the range $0.0\leq x_{\rm
LA}\leq 1.0$. Each configuration is equilibrated for 500 ns and the last 200
ns data is used for the calculation of ${\rm g}_{\rm ij}(r)$ between different
solution components. Here, ${\rm g}_{\rm ij}(r)$s are calculated between the
center-of-masses of the central monomer of the individual trimers. The
simulations are performed at a melt temperature of $T=600$ K.

  * •

Set 2, potential of mean force (PMF) using the umbrella sampling: For these
simulations, we have used monomeric fluid of MMA and LA at a 50-50 mixing
ratio. A configuration consists of 6000 monomers and at $T=600$ K. We have
deliberately chosen a monomeric fluid to avoid any possible conformational
fluctuation that might lead to bias in the calculation of PMF. Here, we note
in passing that PMF obtained from a monomeric system can be translated into
the monomers connected in a chain molecule using a treatment proposed by one
of us with collaborators [34].

  * •

Set 3, morphology of blends: In this case we have chosen $N_{\ell}=30$,
$N_{\rm c}=200$ and five configurations within the range $0.0\leq x_{\rm
LA}\leq 1.0$. Here, the individual configurations are equilibrated for
$1~{}\mu$s. The simulations are performed for the melt state at $T=600$ K.

  * •

Set 4, mechanics and thermal transport calculations: In this step, we take the
equilibrated configurations from the set 3 and then instantaneously quench
these systems to $T=300$ K where the stress-strain behavior and the thermal
transport coefficients $\kappa$ are calculated.

During $NpT$ simulations, the temperature is imposed using the velocity
rescale thermostat with a coupling constant of 1 ps [35]. The pressure is set
to 1 atm with a Berendsen barostat using a time constant 0.5 ps [36].
Electrostatics are treated using the particle-mesh Ewald method [37]. The
interaction cutoff for the nonbonded interactions is chosen as $r_{c}=1.0$ nm,
which is of the same range as the typical correlation length in these systems
[10]. The simulation time step is chosen as $\Delta t=1$ fs and the equations
of motion are integrated using the leap-frog algorithm.

###  II.3 Monomer-monomer potential of mean force

To investigate the relative interaction strengths between different solution
components, we have calculated PMF. Here, PMFs between two monomers is
calculated using the umbrella sampling method [38] over a series of three
independent simulations, i.e., between PMMA-PMMA, PLA-PLA and PMMA-PLA. For
each simulation, a configuration is selected after 50 ns where a monomer type
i is pulled apart from another monomer type j using the steered molecular
dynamics, with a spring constant $k=1000~{}{\rm kJ\ mol^{-1}\ nm^{-2}}$ and at
a pulling rate of $10^{-3}~{}{\rm nm\ ps^{-1}}$ over 1.5 ns. From these
simulations, we have selected twelve configurations within the range 0.35 nm
$<r<$ 1.55 nm and with $\triangle r=0.10$ nm. The force distributions is
obtained for each of the twelve distances over 5 ns long simulation each.
Subsequently, the weighted histogram analysis method (WHAM) was used to
combine the force distributions into a potential of mean force (PMF) [39, 40].

###  II.4 Thermal transport calculations

$\kappa$ is calculated using the approach-to-equilibrium method [41]. For this
purpose, we divide the simulation box into three compartments along the
$x-$direction. Here, the middle slab of width $L_{x}/2$ is sandwiched between
two side slabs each having the width of $L_{x}/4$. Here, we employ periodic
boundary conditions in all directions.

Simulations are performed in two stages: during the initial stage of canonical
simulations the middle slab is kept at $T_{\rm Hot}=350$ K, while the two side
slabs are maintained at $T_{\rm Cold}=250$ K. This simulation is performed for
5 ns each configuration. After this stage, the temperatures are allowed to
relax during a microcanonical run for 50 ps.

From the relaxation of $T_{\rm Hot}-T_{\rm Cold}$ during the microcanical
simulations, we extract a bi-exponential fit and obtain the longitudinal time
constants $\tau_{x}$. $\kappa$ can then be calculated using
$\kappa=C_{v}L_{x}/4\pi^{2}L_{y}L_{z}\tau_{x}$. Here, we consider the Dulong-
Petit classical estimate of the specific heat, i.e., $C_{v}=3N_{\rm
total}k_{\rm B}$, where $N_{\rm total}$ is the total number of atoms in the
simulation box.

##  III Results and discussions

###  III.1 Radial distribution function and Kirkwood-Buff integral

In Fig. 2 we show ${\rm g}_{\rm ij}(r)$ and ${G}_{\rm ij}(r)$ for different
$x_{\rm LA}$.

Figure 2: The pair distribution function ${\rm g}_{\rm ij}(r)$ (a-c) and the
cumulative Kirkwood-Buff integral $G_{\rm ij}(r)$ (d-f) between different
solution components of poly(methyl methacrylate) (PMMA) and poly(lactic acid)
(PLA) blends for different PLA mole fractions $0.0\leq x_{\rm LA}\leq 1.0$.
Data is shown from the simulations of trimers and ${\rm g}_{\rm ij}(r)$ and
$G_{\rm ij}(r)$ are calculated between the central monomers, namely MMA-MMA (a
and d), MMA-LA (b and e) and LA-LA (c and f). The temperature is set at
$T=600$ K and at 1 atm pressure.

It can be appreciated that the LA$-$LA structures are mostly affected by the
change in $x_{\rm LA}$, see Figs. 2(c) and 2(f), while the structures between
the other pairs remain mostly unaffected within the error bar. To obtain
$G_{\rm ij}$ we average individual $G_{\rm ij}(r)$ within the range 1 nm $<r<$
2 nm. The deviation from the plateau for $r>$ 2 nm is a typical behavior
observed for $G_{\rm ij}(r)$ in the mid-sized closed boundary simulation
setups [28]. This is because the aggregation of one molecular species, in one
region of the simulation domain, leads to the depletion of the same species
elsewhere. In this context, ideally ${G}_{\rm ij}(r)$ should be calculated in
the grand canonical ensemble [21, 22, 28]. However, since most simulations are
performed in the mid-sized closed boundary domains, ${G}_{\rm ij}(r)$ often
suffers from the severe system size effects. Here, we find that $N_{\rm
c}=1000$ is a good estimate to obtain the converged ${G}_{\rm ij}(r)$, where
the equilibrium box sizes for all $x_{\rm LA}$ values are over 8 nm. These box
dimensions are over five times the typical correlation lengths in these
systems, i.e., $r\simeq 1.5$ nm. We also note in passing that there are more
elegant methods to estimate ${G}_{\rm ij}$ within a semi-grand canonical
simulation setup [28] and from the particle fluctuations [23], which is,
however, beyond the scope of the present study.

Using the estimates of ${G}_{\rm ij}$ we obtain $\Delta{\mathcal{G}}_{\rm
mix}$ by double integrating Eq. 3 using the boundary condition, i.e.,
$\Delta{\mathcal{G}}_{\rm mix}=0$ at $x_{\rm LA}=0.0$ and 1.0. We have also
calculated $\Delta{H}_{\rm mix}$ per monomer from $NpT$ simulations using Eq.
4. From these two estimates we obtain $T\Delta S_{\rm mix}$.

Figure 3: Parts (a) and (b) show the entropic $T\Delta S_{\rm mix}$ and the
enthalpic $\Delta{H}_{\rm mix}$ contributions to the change in the Gibbs free
energy of mixing $\Delta{\mathcal{G}}_{\rm mix}$, respectively. Part (c) shows
the preferential solvation parameter $\eta$. The data sets are shown for
different lactic acid (LA) mole fractions $x_{\rm LA}$. The lines in parts (a)
and (c) are drawn to guide the eye, while the line in the part (b) is a fit to
the data following $\chi x_{\rm LA}\left(1-x_{\rm LA}\right)$ with
$\chi=0.021k_{\rm B}T$.

In Figs. 3(a-b) we show $T\Delta S_{\rm mix}$ and $\Delta{H}_{\rm mix}$ as a
function of $x_{\rm LA}$. It can be seen that the entropic contribution (that
favors mixing) dominates over the enthalpy (that favors demixing because it is
positive). This indicate that the mixing of PMMA-PLA systems is entropically
stabilized. This is also consistent with the recent experimental results [13].
From the data in Fig. 3(b), we also obtain $\chi-$parameter using the FH
function $\Delta{H}_{\rm mix}=\chi x_{\rm LA}\left(1-x_{\rm LA}\right)$. Here,
we find $\chi=0.021~{}k_{\rm B}T$ at the monomer level, which might look as a
small value. However, for large $N_{\ell}$ this can lead to a significantly
large $N_{\ell}\chi$. Additionally, $T\Delta S_{\rm mix}\propto 1/N_{\ell}$,
see Eq. 1. Therefore, within the long chain limit, one is expected to observe
a clear demixing in these systems, as predicted by the FH theory [19, 20].
However, for the moderate $N_{\ell}$ values, weak miscibility is expected. We
will come back to this point in the next subsection.

A fit to the data in Fig. 3(b) gives $\chi$. This is a good estimate so long
as the mixing of two species does not introduce any nonlinear variation in the
total number density $\rho$. However, in most chemical systems, as in the case
of PMMA-PLA blends, $\rho$ values usually show a weak non-linear variation
with $x_{\rm LA}$ (data not shown). Therefore, a more reasonable picture will
be to have a concentration dependent $\chi-$parameter. In this context, a
quantity that can give an estimate of $x_{\rm LA}$ dependent effective
relative interaction strengths is the preferential solvation parameter $\eta$.
In Fig. 3(c) we shown $\eta$ as a function of $x_{\rm LA}$, which also
indicates that $\eta>0$ for all $x_{\rm LA}$.

Figure 4: Potential of mean force $u^{\rm PMF}_{\rm ij}(r)$ between different
solution components of methyl methacrylate (MMA) and lactic acid (LA) at a LA
mole fractions $x_{\rm LA}=0.5$ and at a temperature $T=600$ K. Symbols are
only shown for every 10th data point. The data is shown from the simulations
of the monomeric fluids, where $u^{\rm PMF}_{\rm ij}(r)$ can be translated
from an isolated monomer system to a monomer in a chain using
${\tilde{u}}_{\rm ij}^{\rm PMF}(r)=-k_{\rm B}T[\ln\\{e^{-u_{ij}^{\rm
PMF}(r)/k_{\rm B}T}+1\\}-\ln\left(2\right)]$.

The data in Figs. 3(b-c) also indicate that the interactions between the like
species are preferred over the cross-interactions. To quantify these relative
interaction strengths, we have also calculated the monomer-monomer PMFs
between different solution species $u^{\rm PMF}_{\rm ij}(r)$, see Fig. 4. It
can be appreciated that the interactions between LA$-$LA and MMA$-$MMA is
preferred over the cross-interactions. It is important to highlight that we
have only calculated $u^{\rm PMF}_{\rm ij}(r)$ for $x_{\rm LA}=0.5$ because
pair structures are mostly unaffected by $x_{\rm LA}$, except for ${\rm
g}_{\rm LA-LA}(r)$, and thus $u^{\rm PMF}_{\rm ij}(r)$ is a good estimate at
$x_{\rm LA}=0.5$.

###  III.2 Morphology

Experiments have shown that the PMMA-PLA blends form “so called” co-continuous
phases [16, 13]. The typical length scales in these phases were found to be of
the order of $1-2~{}\mu$m [13]. It has been argued that$-$ while $\Delta
H_{\rm mix}$ favors demixing, the morphologies are entropically stable. This
trend is consistent with our simulation data presented in Fig. 3. Therefore,
it will also be important to study how the observed trends in Fig. 3 can be
used to understand the morphologies of PMMA-PLA mixtures in simulations. Here,
however, it is important to emphasize that by no means we aim to claim
anything close to the experimental length scales (both in terms of the system
dimensions and the chain lengths) due to the obvious limitations associated
with the molecular simulations. The goal of this work is rather to see if the
simulations can show some qualitative features observed in the experiments.

For these simulations, we have chosen $N_{\ell}=30$ for both PMMA and PLA.
Here, PMMA has a persistence length $\ell_{\rm p}\simeq 0.65$ nm and PLA has
$\ell_{\rm p}\simeq 0.50$ nm. This will give $N_{\ell}=12\ell_{\rm p}$ for
PMMA and $N_{\ell}=15\ell_{\rm p}$ for PLA. While $N_{\ell}=12-15\ell_{\rm p}$
might look as too short to study the polymer properties, it is known to
reasonably capture the conformational behavior of a single chain in solution
[31] and in melt [32].

In Fig. 5 we show simulation snapshots for three different $x_{\rm PLA}$.

Figure 5: Simulation snapshots of the poly(methyl methacrylate) (PMMA) and
poly(lactic acid) (PLA) blends for three different PLA mole fractions $x_{\rm
LA}$. The data is shown for a chain length $N_{\ell}=30$ and for the melt
state at a temperature of $T=600$ K. PMMA and PLA chains are rendered in blue
and yellow, respectively. The system dimensions $L$ for $x_{\rm LA}=0.25$,
0.50, and 0.75 are 9.76 nm, 9.45 nm and 9.12 nm, respectively.

Configurations are shown after $1~{}\mu$s equilibration for each system. We
would also like to highlight that the typical end-to-end relaxation time of a
chain in a PMMA-PLA melt is only about 20$-$25 ns for $N_{\ell}=30$ and thus
the individual chains have moved 40$-$50 times its own size $R_{\rm g}\simeq
1.5$ nm over $1~{}\mu$s. From the simulations, it is evident that the systems
do not phase separate with the well-defined PLA and PMMA domains, rather form
phases with PLA forming percolating network-like structures within the blends.
While it is quite evident that the study of phase separation within an all-
atom setup is nontrivial task especially for the long chain molecules, our
system dimensions are only about $6R_{\rm g}$ and thus the systems show weak
miscibility that is entropically stabilized. Another quantity that speaks in
the favor that the systems do not phase separate (as predicted by Fig. 3 and
Fig. 4) is that $G_{\rm ij}(r)$ (in Figs. 2(d-f)) show convergence. For the
phase separated systems, $G_{\rm ij}(r)$ should diverge for the large $r$
values.

We also note in passing that a completely immiscible polymer blend, especially
consisting of the brittle polymers, would have been extremely unstable under
any external deformation. Therefore, to further consolidate our argument of
PMMA-PLA phase behavior with percolating PLA chains, we will now investigate
the mechanical response of the configurations presented in Fig. 5.

###  III.3 Mechanical response and thermal transport

In Fig. 6 we show the stress-strain behavior for three different $x_{\rm PLA}$
in their glassy state at $T=300$ K. It can be appreciated that the individual
data sets clearly show the typical behavior: namely a linear response for
$\varepsilon<3\%$, a yield peak at around $\varepsilon\simeq 15\%$, a plastic
deformation for $\varepsilon>20\%$, and an ultimate failure with $x_{\rm
PLA}-$dependent strain-to-break $\varepsilon_{\rm f}$, which is defined as
strain where the $\sigma$ values depart from the plateau in Fig. 6. Here,
$\varepsilon_{\rm f}=0.34$ for $x_{\rm PLA}=0.25$, $\varepsilon_{\rm f}=0.55$
for $x_{\rm PLA}=0.50$ and $\varepsilon_{\rm f}=1.38$ for $x_{\rm PLA}=0.75$
[42].

Figure 6: Stress $\sigma$ as a function of strain $\varepsilon$ of the
poly(methyl methacrylate) (PMMA) and poly(lactic acid) (PLA) blends for three
different PLA mole fractions $x_{\rm PLA}$. The simulation box is subjected to
an uniaxial tensile deformation along the $z-$axis with a velocity of
$10^{-3}$ nm/ps. The background shaded data sets are from the raw data and for
the representation purpose we also draw the block averaged data (solid lines).
The data is shown at a temperature $T=300$ K.

Fig. 6 also shows that the systems become more ductile with increasing $x_{\rm
PLA}$, attaining the largest ductility for $x_{\rm LA}=0.75$ with
$\varepsilon_{\rm f}>100\%$. This is also evident from the simulation
snapshots in Fig. 7.

Figure 7: Simulation snapshots showing the undeformed and at the onset of
fracture configurations for three different PLA mole fractions $x_{\rm PLA}$.
The data is shown at a temperature $T=300$ K.

The observation in Fig. 6 is in qualitative agreement with the experimentally
observed trends [13].

What causes this increasing ductility with the increase in $x_{\rm PLA}$? To
discuss this point, we will now put forward different scenarios:

It is a common practice that a brittle polymeric material (such as PMMA) is
usually plasticized by blending in another polymer with a lower glass
transition temperature $T_{\rm g}$ [43, 44]. In that case, the enhancement of
ductility would then be due to the possible enhancement of mobility in the
vicinity of $T_{\rm g}$. Here, a pure PLA material has a lower $T_{\rm g}$
than a pure PMMA sample [45]. Furthermore, $T_{\rm g}$ changes monotonically
with $x_{\rm PLA}$ [18, 13], i.e., the larger the $x_{\rm PLA}$, the lower the
$T_{\rm g}$. Within this argument, if $T$ of a glassy state is chosen such
that $T_{\rm g}-T$ is small (as in the case of the PMMA-PLA blends, i.e., 35 K
$<\left(T_{\rm g}-T\right)<$ 75 K) and thus the elastic response gets highly
dependent on $x_{\rm PLA}$ (or $T_{\rm g}-T$). Within this picture, if $T_{\rm
g}-T$ is increased, a system might become less ductile and can even show
$\varepsilon_{\rm f}$ closer to that of as observed for the smaller $x_{\rm
PLA}$. To test this hypothesis, we have performed one set of simulations with
$T=200$ K and $x_{\rm PLA}=0.75$. In this case $\varepsilon_{\rm f}$ changes
only about 5% with respect to $T=300$ K. Here, it is important to highlight
that a system deep in the glassy state, as in the case for $T=200$ K, there
still exists a contrast in the elastic modulus between the two pure phases,
which is, however, only about 15-20% between a pure PMMA and a pure PLA
systems and thus is not so significant to make a large difference in the
ductility. These observations speak against any explanation of increased
ductility purely on $T_{\rm g}$.

Another possible scenario is based on the details of chemical structures of
individual monomer units, see Fig. 1. In this case, the side groups of PMMA
can strongly interact and thus increase inter-chain interlocking, thus makes
PMMA a brittle system. On the contrary, from a purely structural point of
view, a PLA chain can slide past another PLA or a PMMA chain with rather small
friction. Therefore, increasing $x_{\rm PLA}$ can withstand larger
$\varepsilon$ upon deformation, while the chains percolate to form the
elongated network-like structures. This is seen by comparing two snapshots in
Fig. 7(c). This suggests a delicate interplay between the monomer structures,
the local packing and the mechanical response.

Here it is also important to highlight that upon deformation small free spaces
appear that coalesce to fracture the samples [5]. This coalescence happens at
a larger $\varepsilon$ with increasing $x_{\rm PLA}$ values. Furthermore, PLA
chains can also join hands (or intermingle) because of their preferred
monomer-monomer interaction, see Fig. 4. Therefore, a PLA chain can move more
easily with respect to a PMMA chain and thus can rearrange itself more easily
when put under a tensile deformation. To investigate this aspect, we have also
performed two additional simulations with shorter PLA chains $N_{\rm PLA}=10$
and 20, while keeping the chain length of PMMA fixed, i.e., $N_{\rm PMMA}=30$.
Note that our chain lengths are short and therefore we are deliberately
avoiding to use “entanglement” in our nomenclature above (instead used
intermingle), which is a well defined concept in the long chain polymer melts.

In Fig. 8, we show the typical stress-strain behavior for three different
$N_{\rm PLA}$ and for $x_{\rm PLA}=0.75$.

Figure 8: Same as Fig. 6, however, for one poly(lactic acid) (PLA) mole
fraction $x_{\rm PLA}=0.75$ and for three different PLA chain lengths $N_{\rm
PLA}$.

It can be appreciated that $\varepsilon_{\rm f}$ shows a non-monotonic
dependence with $N_{\rm PLA}$, i.e., going from $N_{\rm PLA}=30$ to 20,
$\varepsilon_{\rm f}$ first decreases by about 10-15% and then for $N_{\rm
PLA}=10$, $\varepsilon_{\rm f}$ increases again. A simple picture would then
be$-$ increasing $N_{\rm PLA}$ mean that a chain will read a fully stretched
conformation at a larger $\varepsilon_{\rm f}$. This, however, contradicts the
increased ductility at $N_{\rm PLA}=10$. Here it noteworthy that in the case
of perfect miscibility this non-monotonic ductility would be quite obvious
because of the chain bridging/gluing the broken interfaces and chain
rearrangements, a similar scenario as observed in the thermoset-thermoplastic
blends [46]. On the contrary, however, PLA-PMMA blends are partially miscible
as expected from the solvation thermodynamics, see Fig. 3. In this context, we
speculate that the significant increase in $\varepsilon_{\rm f}$ for $N_{\rm
PLA}=10$ (which is only about 3 persistence lengths $\ell_{\rm p}$) might be
due to the fact that the shorter chains act like single particles and thus can
rearrange faster, thus hindering the coalescence of the free volumes within
the samples. These results highlight the importance of chain friction in the
brittle-to-ductile transition of PLA-PMMA blends. We also note in passing
that$-$ the numerical validation of the above scenarios will require a set of
careful investigations with a variety of monomer structures, monomer-monomer
interactions, polymer blends, chain lengths, and mixing ratios. In particular,
to investigate how the delicate balance between the molecular flexibility, the
monomer structure, the detailed interactions and the blend morphologies can
dictate the mechanical response of the commodity blends, which will, however,
be presented elsewhere [47].

Figure 9: The main panel shows the thermal transport coefficient $\kappa$ as a
function of poly(lactic acid) (PLA) mole fraction $x_{\rm PLA}$. For
comparison, we have also included experimental data for PMMA [48] and PLA [45]
and also simulation data for PMMA from an earlier work using the Müller-Plathe
method [32]. Inset shows the variation of $\kappa$ with the elastic modulus
$E$. Lines are drawn to guide the eye. The data is shown for a temperature
$T=300$ K.

Another quantity that can be extracted from Fig. 6 is the elastic modulus $E$
from the linear regime, i.e., for $\varepsilon<3\%$. In our simulations, going
from pure PLA to PMMA, we find 3.0 GPa $<E<$ 4.3 GPa. These values are also
within the similar range as observed in the experiments [45, 48]. A quantity
that is directly related to $E$ is the thermal transport coefficient $\kappa$
of the amorphous polymers, as predicted by the minimum thermal conductivity
model $\kappa\propto E$ [49]. Therefore, we have also calculated $\kappa$ and
the corresponding data is shown in the main panel of Fig. 9 and the dependence
of $\kappa$ on $E$ is shown in the inset of Fig. 9. As expected from the
blending of two materials, $\kappa$ with $x_{\rm PLA}$ shows a monotonic
variation between the $\kappa$ values of a pure PLA and a pure PMMA systems.
While $\kappa$ values are known for the pure phases on PMMA and PLA, which is
also consistent with the earlier experiments [45, 48] and simulations [32], we
report $\kappa$ data for PMMA-PLA blends here. This shows an internal
correlation between the linear elasticity, mechanics, and thermal transport of
commodity polymer blends. Here, we also note in passing that the earlier
studies have reported that the symmetric blends of the long chain
macrmolecules can show the monotonic variation in $\kappa$ between two pure
phases [48], while a weak non-monotonic variation in $\kappa$ with relative
concentration was observed for the asymmetric blends [10]. Here, PMMA-PLA
blends fall under the former catagory.

##  IV Conclusions

Using large scale molecular dynamics simulation of an all-atom model, we have
established an internal correlation between the solvation thermodynamics, the
blend morphology, the non-linear mechanics and the thermal transport in a
commodity polymer blend. For this purpose, we have chosen a technologically
relevant blends of the commodity polymers, namely poly(methyl methacrylate)
(PMMA) and poly(lactic acid) (PLA). Studying the thermodynamics in polymer
mixtures is rather nontrivial task. Therefore, we make use of the fluctuation
theory of Kirkwood and Buff [21] that connects the pair-wise structure to the
mixture thermodynamics. We highlight the importance of the delicate balance
between the entropic $T\Delta S_{\rm mix}$ and the enthalpic $\Delta H_{\rm
mix}$ contributions to the mixing Gibbs free energy change
$\Delta{\mathcal{G}}_{\rm mix}$. We further show how these contributions can
be individually extracted within the molecular dynamics simulations,
especially for the long chain macromolecular mixtures. These observations are
correlated with the morphologies of PLA-PMMA blends and finally to the
mechanical and thermal responses.

Even when we have only presented results for one system, our work achieves
three key highlights: (1) Establish a structure-property relationship within
one unified simulation framework, (2) a possible explanation of the unexpected
mechanical response of the PMMA-PLA blends, as observed in a recent experiment
[13], and (3) provide a set of possible scenarios that can serve as a guiding
tool to understand polymer blends. Moreover, further works need to perform on
a broad range of polymer blends with different microscopic interactions,
miscibility, monomeric structures and chain conformations. However, the hand-
in-hand correlation presented in this work between different physical
properties is rather rare and the protocol presented may be used to design the
multi-component materials for their advanced functional applications.

##  V Acknowledgement

D.M. further thanks the ARC Sockeye facility of the University of British
Columbia, the Compute Canada facility and the Quantum Matter Institute LISA
cluster, where major part of the simulations were performed. This research was
undertaken thanks, in part, to the Canada First Research Excellence Fund
(CFREF), Quantum Materials and Future Technologies Program. T.E.O. thanks the
Centro Nacional de Supercomputção (CESUP) for providing HPC resources, where
part of the simulations were performed. C.R. thanks MITACS for the financial
support.

## References

  * [1] M. Köger, Physics Reports 390, 453 (2004).
  * [2] M. Müller, Progress in Polymer Science 101, 101198 (2020).
  * [3] D. Mukherji, C. M. Marques, and K. Kremer, Annual Reviews of Condensed Matter Physics 11, 271 (2020).
  * [4] M. J. Stevens, Macromolecules 34, 2710 (2001).
  * [5] D. Mukherji and C. F. Abrams, Phys. Rev. E 78, 050801 (2008).
  * [6] B. J. Blaiszik, S. L. B. Kramer, S. C. Olugebefola, J. S. Moore, N. R. Sottos, S. R. White, Annual review of materials research 40, 179 (2010).
  * [7] M. Sharifi, C. W. Jang, C. F. Abrams and G. R. Palmese, J. Mater. Chem. A 2, 16071 (2014).
  * [8] G. Kim, D. Lee, A. Shanker, L. Shao, M. S. Kwon, D. Gidley, J. Kim, and K. P. Pipe, Nat. Mater. 14, 295 (2015).
  * [9] X. Xie, D. Li, T. Tsai, J. Liu, P. V. Braun, and D. G. Cahill, Macromolecules 49, 972 (2016).
  * [10] D. Bruns, T. E. de Oliveira, J. Rottler, and D. Mukherji, Macromolecules 52, 5510 (2019).
  * [11] N. Kim, B. Domercq, S. Yoo, A. Christensen, B. Kippelen, and S. Graham, Appl. Phys. Lett. 87, 241908 (2005).
  * [12] M. K. Smith, V. Singh, K. Kalaitzidou, and B. A. Cola, ACS Appl. Mater. Int. 8, 14788 (2016).
  * [13] X. Hou, S. Chen, J. J. Koh, J. Kong, Y.-W. Zhang, J. C. C. Yeo, H. Chen, and C. He ACS Macro Lett. 10, 406 (2021).
  * [14] J. C. C. Yeo, J. K. Muiruri, J. J. Koh, W. Thitsartarn, X. Zhang, J. Kong, T. T. Lin, Z. Li, C. He, Adv. Funct. Mater. 30, 2001565 (2020).
  * [15] P. Pötschke and D. R. Paul, J. Macromol. Sci., Polym. Rev. 43, 87 (2003).
  * [16] K.-P. Le, R. Lehman, J. Remmert, K. Vanness, P. M. L. Ward, and J. Idol, J. Bio. Sci. Pol. Ed. 17, 121 (2006).
  * [17] C. Samuel, J.-M. Raquez, and P. Dubois, AIP Conf. Pro. 1664, 030005 (2015).
  * [18] M. Gonzalez-Garzon, S. Shahbikian, and M. A. Huneault, J. Pol. Res. 25, 58 (2018).
  * [19] M. Doi, Soft Matter Physics Oxford University Press (2013).
  * [20] P.-G. de Gennes, Scaling Concepts in Polymer Physics Cornell University Press, London (1979).
  * [21] J. G. Kirkwood and F. P. Buff, J Chem. Phys. 19, 774 (1951).
  * [22] A. Ben Naim, Molecular Theory of Solutions Oxford University Press, London (2006).
  * [23] R. Cortes-Huerto, K. Kremer, and R. Potestio, J Chem. Phys. 145, 141103 (2006).
  * [24] A. A. Galata, S. D. Anogiannakis, and D. N. Theodorou, Fluid Phas. Equil. 470, 25 (2018).
  * [25] M. Sevilla and R. Cortes-Huerto, arXiv:2110.14536 (2021).
  * [26] M. Kang and P. E. Smith, Fluid Phas. Equil. 256, 14 (2007).
  * [27] V. Pierce, M. Kang, M. Aburi, S. Weerasinghe, and P. E. Smith, Cell Biochem. Biophys. 50, 1 (2008).
  * [28] D. Mukherji and K. Kremer, Macromolecules 46, 9158 (2013).
  * [29] P. Kumari, V. V. S. Pillai, D. Gobbo, P. Ballone, and A. Benedetto, Phys. Chem. Chem. Phys. 23, 944 (2021).
  * [30] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J Am. Chem. Soc. 118, 11225 (1996).
  * [31] D. Mukherji, C. M. Marques, T. Stühn, and K. Kremer, Nat. Commun. 8 1374 (2017).
  * [32] C. Ruscher, J. Rottler, C. E. Boott, M. J. MacLachlan, and D. Mukherji, Phys. Rev. Mat. 3, 125604 (2019).
  * [33] S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl, Bioinformatics 29, 845 (2013).
  * [34] C. De Silva, P. Leophairatana, T. Ohkuma, J. T. Koberstein, K. Kremer, and D. Mukherji, J Chem. Phys. 147, 064904 (2017).
  * [35] G. Bussi, D. Donadio, and M. Parrinello. Canonical sampling through velocity rescaling. J Chem. Phys. 126, 014101 (2007).
  * [36] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J Chem. Phys. 81, 3684 (1984).
  * [37] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. A smooth particle mesh ewald method. J Chem. Phys. 103, 8577 (1995).
  * [38] G. M. Torrie, J. P. Valleau J. Comput. Phy. 23, 187 (1977).
  * [39] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, J. M. Rosenberg J. Comput. Chem. 13, 1011 (1992).
  * [40] J. S. Hub, B. L. de Groot, D. van der Spoel J. Chem. Theory Comput. 6, 3713 (2010).
  * [41] E. Lampin, P. L. Palla, P. A. Francioso, and F. Cleri, Thermal conductivity from approach-to-equilibrium molecular dynamics. J. Appl. Phys. 114 033525 (2013).
  * [42] Since the simulations domains are rather small (only about 8$-$9 nm), we have also performed simulations when the systems are elongated along the $y-$ and the $z-$directions to investigate if the deformation along only one direction introduces any morphological bias on the mechanical response. The mechanical responses along all three different directions only varry by less than 2%.
  * [43] A. Bouzouita, Elaboration of polylactide-based materials for automotive application : study of structure-process-properties interactions PhD thesis (2016).
  * [44] R. Long, C.-Y. Hui, J. P. Gong, and E. Bouchbinder, Annual Reviews of Condensed Matter Physics 12, 71 (2021).
  * [45] J. Brandrup, E. H. Immergut, and E. A. Grulke, Polymer Handbook, 1079 2 Volumes Set, 4th ed. Wiley (2003).
  * [46] D. Mukherji and C. F. Abrams, Euro Phys. Lett. 88, 56001 (2009).
  * [47] D. Mukherji et al, in preparation (2021).
  * [48] X. Xie, D. Li, T. Tsai, J. Liu, P. V. Braun, and D. G. Cahill, Macromolecules 49, 972 (2016).
  * [49] D. G. Cahill, S. K. Watson, and R. O. Pohl, Phys. Rev. B 46, 6131 (1992).


Jasper Hoffmann<EMAIL_ADDRESS>
Diego Fernandez Clausen<EMAIL_ADDRESS>
Julien Brosseit<EMAIL_ADDRESS>
Neurorobotics Lab, University of Freiburg, Germany. and Julian Bernhard
<EMAIL_ADDRESS>
Klemens Esterle<EMAIL_ADDRESS>
Moritz Werling<EMAIL_ADDRESS>
Michael Karg<EMAIL_ADDRESS>
BMWGroup, Unterschleissheim, Germany. and and Joschka Bödecker
<EMAIL_ADDRESS>
Neurorobotics Lab, University of Freiburg, Germany.

# PlanNetX: Learning an Efficient Neural Network Planner from MPC for
Longitudinal Control

###### Abstract

Model predictive control (MPC) is a powerful, optimization-based approach for
controlling dynamical systems. However, the computational complexity of online
optimization can be problematic on embedded devices. Especially, when we need
to guarantee fixed control frequencies. Thus, previous work proposed to reduce
the computational burden using imitation learning (IL) approximating the MPC
policy by a neural network. In this work, we instead learn the whole planned
trajectory of the MPC. We introduce a combination of a novel neural network
architecture PlanNetX and a simple loss function based on the state trajectory
that leverages the parameterized optimal control structure of the MPC. We
validate our approach in the context of autonomous driving by learning a
longitudinal planner and benchmarking it extensively in the CommonRoad
simulator using synthetic scenarios and scenarios derived from real data. Our
experimental results show that we can learn the open-loop MPC trajectory with
high accuracy while improving the closed-loop performance of the learned
control policy over other baselines like behavior cloning.

###### keywords:

imitation learning, model predictive control, autonomous driving

##  1 Introduction

Motion planning for autonomous driving is challenging due to the high
uncertainty of the surrounding traffic requiring constant re-planning to
adjust the plan to newly obtained information. Model predictive control (MPC)
is an optimization-based approach that relies on rapidly resolving an optimal
control problem (OCP) numerically, planning a trajectory at each step.
However, one challenge lies in the worst-case solving time of the optimization
problem which is crucial for safety critical systems such as autonomous
driving, where a control update is required after a fixed time interval.

To address this problem, a lot of effort was put into computing the MPC
solutions offline and storing them in an appropriate representation to allow
fast online evaluation during deployment. One example of this is _explicit
MPC_. In the case of linear MPC, the MPC policy, defined as the first control
of the MPC solution, is a piece-wise linear function that can be stored as a
look-up table [Bemporad et al.(1999)Bemporad, Borrelli, and Morari]. However,
for non-linear dynamics or constraints representing the MPC policy exactly is
challenging and often done approximately [Johansen and Grancharova(2003)].
Another prominent example is imitation learning (IL). Here, we try to
approximate the MPC policy by using function approximators like neural
networks. Previous work focused on improving the safety aspects of the learned
policy or replacing the loss function of standard behavior cloning, which
tries to minimize the Euclidean distance to the first control of the MPC
solution, using objectives that are aware of the underlying OCP structure of
the MPC.

In this work, instead of just learning the first control we want to learn the
whole planned trajectory of the MPC solution. This is important in autonomous
driving where we use the planned trajectory as a tracking signal for low-level
controllers that might directly control the motors. Additionally, we find that
learning the whole planned trajectory allows us to define a simple loss
formulation that uses the underlying OCP structure during training being
potentially aware of the network’s approximation errors. Furthermore, learning
the trajectory could be used to warm start the MPC solver for significantly
faster convergence and better safety guarantees. The goal of this work is to
achieve a higher planning frequency in comparison to numerically solving the
OCP. Thus, we further investigate common neural network compression techniques
like pruning and quantization to reduce the inference time.

The contributions of this work are threefold: (1) We propose an imitation
learning framework called PlanNetX specifically for learning the planned
trajectory of an MPC planner, consisting of a simple loss function and a
specialized neural network architecture. (2) We validate the approach by
learning a longitudinal planner for autonomous driving and extensively testing
it on a handcrafted CommonRoad, [Althoff et al.(2017)Althoff, Koschi, and
Manzinger], benchmark set consisting of real and synthetic driving data. (3)
We investigate how common neural network compression techniques like pruning
and quantization can reduce the inference time even further.

###  1.1 Related Work

In the context of IL from an MPC policy, various aspects were considered: (1)
_Safety and Stability_ : [Mamedov et al.(2022)Mamedov, Reiter, Diehl, and
Swevers] use a safety filter formulation similar to [Wabersich and
Zeilinger(2021)] to guarantee the safety of the final control policy. Using
the concept of control-barrier functions, [Cosner et al.(2022)Cosner, Yue, and
Ames] show that under certain conditions one can transfer the safety
guarantees of a robust MPC expert to the learned policy. [Schwan et
al.(2023)Schwan, Jones, and Kuhn] provide an upper bound on the approximation
error of the learned policy to the original MPC policy by solving a mixed-
integer formulation. (2) _Improved Loss Function_ : Further research tries to
replace the surrogate loss of standard behavior cloning, which is often the
Euclidean distance, with an objective that takes into account the underlying
OCP problem: In [Carius et al.(2020)Carius, Farshidian, and Hutter, Reske et
al.(2021)Reske, Carius, Ma, Farshidian, and Hutter] a loss based on the
control Hamiltonian for continuous-time OCP is introduced, whereas in [Ghezzi
et al.(2023)Ghezzi, Hoffman, Frey, Boedecker, and Diehl] a Q-loss is
introduced which calculates the value of a fixed first control by solving the
rest of the discrete-time OCP.

Imitation learning from trajectory optimization was considered in [Mordatch
and Todorov(2014)] where an alternating direction method of multipliers (ADMM)
method jointly optimizes a task and policy reconstruction cost. Related to
this approach, is the work in [Xie et al.(2021)Xie, Li, Van Wyk, Dellaert,
Boots, and Ratliff] where an additional trajectory is optimized to follow the
expert trajectories, which is then used to optimize a policy in an optimal
tracking formulation. In [Levine and Koltun(2013a), Levine and Koltun(2013b)]
an iterative linear quadratic regulator (iLQR) [Li and Todorov(2004)] is used
for guiding a reinforcement learning policy to low costs regions. However, all
the mentioned approaches just focus on learning the control policy but not the
planned trajectory.

##  2 Preliminaries

In this section, we give a quick introduction to model predictive control and
imitation learning.

###  2.1 Model Predictive Control

In the following, we introduce the notion of a parameterized discrete-time
OCP. At each time step, an MPC planner observes a current state
$\bar{x}_{0}\in\mathbb{R}^{n_{x}}$ to generate an optimal control and state
trajectory $U\in\mathbb{R}^{n_{u}\times N}$ and
$X\in\mathbb{R}^{n_{x}\times{N+1}}$ by solving the OCP [2] X, UT(x_N) +
∑_k=0^N-1 L_k(x_k, u_k) x_k+1= f(x_k, u_k), k=0,…, N-1 h(x_k, u_k, p_k)≤0,
k=0,…, N-1 x_0= ¯x_0, with a planning horizon length $N$, stage costs $L_{k}$,
terminal costs $T$, and known, potentially simplified, system dynamics $f$.
Furthermore, in (2.1) we assume that we can parameterize the inequality
constraints $h$ at each planning step with additional parameters
$P\in\mathbb{R}^{n_{p}\times N}$. For the rest of the work, we denote
$X^{\star}\in\mathbb{R}^{n_{x}\times{N+1}}$ and
$U^{\star}\in\mathbb{R}^{n_{u}\times N}$ as the optimal solution found by the
MPC when solving the above OCP.

###  2.2 Imitation Learning

With imitation learning (IL) we try to approximate the behavior of an expert
by using a function approximator like a neural network. In the context of this
work, the MPC planner is the expert. We differentiate between learning the
_control policy_ or control law $\bar{x}_{0}\mapsto u_{0}^{\star}$ of the MPC,
which is already enough to control the system, and learning the _planning
policy_ of the MPC planner, where we want to approximate the mapping
$\bar{x}_{0}\mapsto\left(X^{\star},U^{\star}\right)$. Note, that the second
formulation contains the first one as a sub-problem. In more mathematical
terms, we want to approximate either the MPC control or plan by the following
parameterized policies

$\displaystyle\pi_{\theta}^{\textrm{policy}}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{p}\times
N}$ $\displaystyle\to\mathbb{R}^{n_{u}},$
$\displaystyle\pi_{\theta}^{\textrm{plan}}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{p}\times
N}$
$\displaystyle\to\mathbb{R}^{n_{x}\times(N+1)}\times\mathbb{R}^{n_{u}\times
N},$ (1) $\displaystyle\bar{x}_{0},P$
$\displaystyle\mapsto\hat{u}_{0}^{\theta}$ $\displaystyle\bar{x}_{0},P$
$\displaystyle\mapsto\hat{X}^{\theta}\times\hat{U}^{\theta}$

where $\theta$ denotes the parameter vector of the policy which could be for
example the weights of a neural network. For the planning policy
$\pi_{\theta}^{\textrm{plan}}$, we assume that
$\hat{x}_{0}^{\theta}=\bar{x}_{0}$. The goal of both learning formulations is
to minimize the following two objectives:

$\displaystyle\mathcal{L}^{\textrm{policy}}(\theta)$
$\displaystyle=\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\ell\left(\hat{u}_{0}^{\theta},\;u_{0}^{\star}\right)\right],$
$\displaystyle\mathcal{L}^{\textrm{plan}}(\theta)$
$\displaystyle=\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\ell\left(\hat{U}^{\theta},\;\hat{X}^{\theta},U^{\star},\;X^{\star}\right)\right],$
(2)

where $\mathcal{D}$ is a given state and parameter distribution over
$\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{p}\times(N+1)}$ whereas $\ell$ is an
arbitrary point wise loss functions.

#### Behavior Cloning:

For learning just the control policy, the simplest formulation is behavior
cloning (BC) which does not use any information about the underlying OCP.
Instead, here we minimize the distance of our predicted control
$\hat{u}_{0}^{\theta}$ to the optimal control $u^{\star}_{0}$ quadratically

$\displaystyle\mathcal{L}^{\textrm{BC}}(\theta)\coloneqq\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\|\hat{u}_{0}^{\theta}-u^{\star}_{0}\|_{2}^{2}\right].$
(3)

Recent work tries to replace the quadratic distance with a distance measure
that considers the underlying OCP structure like in [Ghezzi et
al.(2023)Ghezzi, Hoffman, Frey, Boedecker, and Diehl] or [Carius et
al.(2020)Carius, Farshidian, and Hutter, Reske et al.(2021)Reske, Carius, Ma,
Farshidian, and Hutter].

##  3 Learning a MPC Trajectory Planner

In the following, we want to discuss how we can use the knowledge of the
experts underlying OCP structure to learn the planned trajectory. We first
discuss different loss formulations and then continue with the PlanNetX
architecture tailored to learning the MPC planning policy
$\pi^{\textrm{plan}}$ of (1).

###  3.1 Loss Formulation

Given the initial state $\bar{x}_{0}$ and the dynamics model $f$, predicting
the controls $\hat{U}^{\theta}$ is already enough to also generate the
according state trajectory $\hat{X}^{\theta}$. One idea would be to minimize
the distance $\hat{U}^{\theta}$ between $U^{\star}$, which we define as the
_control trajectory loss_

$\mathcal{L}^{u}(\theta)\coloneqq\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\frac{1}{N}\sum_{k=0}^{N-1}\gamma^{k}\left\lVert\hat{u}_{k}^{\theta}-u_{k}^{\star}\right\rVert_{W}^{2}\right],$
(4)

where $\gamma$ is a discount factor, with $0<\gamma\leq 1$ allowing to put
more weights on early controls in the loss, and $W$ is a positive-definite
weight matrix, defining a norm via $\lVert x\rVert_{W}^{2}=x^{T}Wx$. However,
the loss function defined in (4) does not consider any approximation errors of
the neural network. Thus, if there is a mismatch early in the predicted
trajectory, there is no incentive for the network to correct it in the later
parts of the predicted trajectory.

Instead of minimizing the distance between $\hat{U}^{\theta}$ and $U^{\star}$
in control space, we instead minimize the distance in the state space, thus
between the predicted state trajectory $\hat{X}^{\theta}$ and $X^{\star}$. We
define the _state trajectory loss_ with

$\displaystyle\mathcal{L}^{x}(\theta)\coloneqq\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\frac{1}{N}\sum_{k=1}^{N}\gamma^{k}\left\lVert\hat{x}_{k}^{\theta}-x_{k}^{\star}\right\rVert_{W}^{2}\right].$
(5)

If there are specific constraints and costs regarding the control trajectory,
one could also potentially combine the two loss formulations (4) and (5).
Note, that using the above loss formulations we are implicitly able to learn
the underlying expert cost function and the respective constraints.

###  3.2 PlanNetX

$f$
---
$\bar{x}_{0}$
---
$\hat{u}_{0}$
---
$\pi_{\theta}$
---
$f$
---
$\pi_{\theta}$
---
$\hat{u}_{1}$
---
$\hat{x}_{1}$
---
$p_{0},t_{0}$
---
$p_{1},t_{1}$
---
$f$
---
$\pi_{\theta}$
---
$p_{N-1},t_{N-1}$
---
$\hat{x}_{N-1}$
---
$\hat{u}_{N-1}$
---
PlanNetX
---
PlanNetXEnc
---
$\bar{x}_{0}$
---
$f$
---
$\hat{u}_{0}$
---
$\pi_{\theta}$
---
$f$
---
$\pi_{\theta}$
---
$\hat{u}_{1}$
---
$\hat{x}_{1}$
---
$z_{0}$
---
$z_{1}$
---
$f$
---
$\pi_{\theta}$
---
$z_{N-1}$
---
$\small\textrm{Transformer Encoder}_{\theta}$
---
$p_{0},t_{0}$
---
$p_{1},t_{1}$
---
$p_{N-1},t_{N-1}$
---
$\hat{u}_{N-1}$
---
$\hat{x}_{N}$
---
$\hat{x}_{N-1}$
---
$\dots$
---
$\hat{x}_{N}$
---
Trainable Parts
---
Feature Normalization
---
$p_{k}$
---
Constraint Parameter
---
Figure 1: The two proposed architectures PlanNetX and PlanNetXEnc. The
PlanNetXEnc includes an additional Transformer encoder layer that contains
potentially the information of the parameters of all different time steps.

The motivation for the PlanNetX architecture comes from taking a look at the
original OCP formulation of (2.1) that we solve for MPC inference. The idea is
to replace the free decision variables $X,U$ with predictions of a neural
network policy $\pi$ and use the state trajectory loss (5), by starting a
roll-out from the initial state $\bar{x}_{0}$ and iteratively applying our
learned policy $\pi_{\theta}$ and the dynamics $f$ resulting into the
following formulation:

$\displaystyle\min_{\theta}\quad\mathbb{E}_{\bar{x}_{0},P\sim\mathcal{D}}\left[\frac{1}{N}\sum_{k=1}^{N}\gamma^{k}\left\lVert\hat{x}_{k}^{\theta}-x_{k}^{\star}\right\rVert_{W}^{2}\right]$
(6a)
$\displaystyle\hat{x}_{k+1}=f(\hat{x}_{k},\hat{u}_{k}),\;\hat{u}_{k}=\pi^{\theta}(\hat{x}_{k},z_{k}),\;z_{k}=\psi_{k}^{\theta}(P,\;T),\hat{x}_{0}=\bar{x}_{0}.$
(6b)

The policy $\pi^{\theta}$ takes an information vector $z_{k}$ containing the
necessary information of the constraint parameter $P$. The information vector
is derived from an additional constraint network $\psi^{\theta}_{k}$ that
bundles all the information of the inequality constraint parameter $P$. We
additionally introduce a time vector
$T=(0,t_{d},2t_{d},\dotsc,(N-1)t_{d})^{T}\in\mathbb{R}^{N}$ to give the policy
the information, for which time step it predicts a control, where $t_{d}$
denotes the discretization time of the MPC planner. Note, that optimizing (6)
requires backpropagating through the dynamics function $f$.

We derive two neural network architectures, see figure 1. The PlanNetX
architecture does not use an additional parameter network $\psi_{\theta}$.
Instead, we directly give the parameter $p_{k}$ and time $t_{k}$ as an input,
assuming that to predict $\hat{u}_{k}$ we do not need the full parameter
matrix $P$, but $\hat{x}_{k}$, $p_{k}$ and $t_{k}$ are enough. The second
architecture PlanNetXEnc adds a Transformer encoder as described in [Vaswani
et al.(2017)Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser, and
Polosukhin] to compress the information of the whole parameter matrix $P$ and
the time vector $T$ into a latent vector $Z\in\mathbb{R}^{n_{z}\times N}$.
Here, each parameter and time tuple $p_{k}$ is one token, which is first
embedded with a linear layer like in [Dosovitskiy et al.(2021)Dosovitskiy,
Beyer, Kolesnikov, Weissenborn, Zhai, Unterthiner, Dehghani, Minderer,
Heigold, Gelly, Uszkoreit, and Houlsby] and then processed by the self-
attention layer. The time information is added by a learned positional
encoding, as often used in Transformer models [Devlin et al.(2019)Devlin,
Chang, Lee, and Toutanova]. We further normalize the input features by min-max
normalization where the minimum and maximum values are derived from the values
that would be obtained when solving the OCP numerically.

###  3.3 Compression of the Network

The goal of this work is to reduce the inference time of an MPC planner by
replacing it with a neural network planner. To reduce the inference time
further, we considered pruning and quantization: Pruning removes the weights
of the neural network. We used iterative $\ell^{1}$-_structured pruning_ to
prune the columns of the linear layers of the neural network, [Li et
al.(2017)Li, Kadav, Durdanovic, Samet, and Graf], which prunes the columns
with the smallest $\ell^{1}$-norm of one linear layer. Neural network
quantization is a technique that involves reducing the precision of the
numerical values in a neural network, typically from 32-bit floating-point
numbers to lower-bit fixed-point numbers typically 16-bit floating point or
8-bit integers. This is often done by linearly scaling the calculations into a
range where they are more precise. In this work, we focus on post-training
_static quantization_ which pre-calculates the scaling and re-scaling factors
on the validation dataset $\mathcal{D_{\mathrm{valid}}}$ before online
inference [Gholami et al.(2021)Gholami, Kim, Dong, Yao, Mahoney, and Keutzer].

##  4 Experiments

In the following, we present the empirical results of this work. After
introducing the longitudinal MPC Planner, we continue with our benchmarking
setup, and then present and discuss the results in terms of open-loop planning
performance, closed-loop performance, and worst-case inference time.

###  4.1 Longitudinal Planner

Velocity | Acceleration | Jerk | Speed Limits | Cost Weights | Slack Weights
---|---|---|---|---|---
$\begin{aligned} v_{\textrm{min}}&=0\\\ v_{\textrm{max}}&=50\\\ \end{aligned}$ | $\begin{aligned} a_{\textrm{min}}&=-6.5\\\ a_{\textrm{max}}&=1\\\ \end{aligned}$ | $\begin{aligned} j_{\textrm{min}}&=-8\\\ j_{\textrm{max}}&=16\\\ \end{aligned}$ | $\begin{aligned} v_{\textrm{max},1}&\in\\{50,80,120,130\\}\\\ v_{\textrm{max},2}&\in\\{50,80,120,130\\}\\\ s_{\textrm{change}}&\in[0,\infty)\\\ \end{aligned}$ | $\begin{aligned} w_{s}&=0.4\\\ w_{a}&=0.8\\\ w_{j}&=0.2\\\ w_{u}&=0.2\\\ \end{aligned}$ | $\begin{aligned} w_{\zeta,a_{N}}&=100\\\ w_{\zeta,\textrm{dist}}&=500\\\ \end{aligned}$
Table 1: Longitudinal Planner Specifications

The following MPC formulation is adopted from [Gutjahr et al.(2017)Gutjahr,
Gröll, and Werling] and [Pek and Althoff(2021)]. We define the longitudinal
state of the ego car with $x=[s,v,a,j]^{T}\in\mathbb{R}^{4}$, where $s$ is the
longitudinal position of the front of the car, $v$ is the velocity, $a$ is the
acceleration, and $j$ is the jerk. The input of the system is the snap given
by $u(t)=\ddot{a}(t)\in\mathbb{R}$ which is the derivative of the jerk
resulting in the linear time-invariant system:

$\frac{d^{4}}{dt^{4}}s(t)=u(t).$ (7)

Integrating the dynamics with a time discretization $t_{\textrm{d}}$ and
assuming piecewise constant inputs we get the discrete-time dynamics system

$x_{k+1}=f(x_{k},u_{k})=A_{d}x_{k}+B_{d}u_{k}.$ (8)

The cost function $L$ is a mixture of quadratic costs penalizing high
accelerations, jerk, and controls with weights $w_{a},w_{j},w_{u}>0$ for
comfort and encourages the progress made in the longitudinal direction with a
linear weight $w_{s}>0$, resulting in the following loss

$L_{k}(x,u)=\gamma^{k}\left(w_{a}a^{2}+w_{j}j^{2}+w_{u}u^{2}-w_{s}s\right).$
(9)

An additional discount factor $\gamma$, can be used to favor cost minimization
of short-term over long-term costs.

Figure 2:  Exemplary cut-in scenario of the CommonRoad benchmark created from
recordings of the HighD dataset. The green car is the ego car and the yellow
region is the goal position the ego car needs to reach. The surrounding
traffic is depicted in blue.

To guarantee the feasibility of the planned trajectory, we assume the
following box constraints for the velocity $v_{\textrm{min}}\leq v_{k}\leq
v_{\textrm{max}}$, acceleration $a_{\textrm{min}}\leq a_{k}\leq
a_{\textrm{max}}$ and jerk $j_{\textrm{min}}\leq j_{k}\leq j_{\textrm{max}}$
for $k=0,\dotsc,N$. We further added a terminal acceleration box constraint,
$0\leq a_{N}\leq 0$, such that the final acceleration is $0$. We describe the
rear position of the _lead vehicle_ that is in front of us by
$\bar{x}^{\textrm{lead}}_{0}=[s^{\textrm{lead}},v^{\textrm{lead}},a^{\textrm{lead}}]^{T}$
and use a simple forward prediction
$\hat{X}^{\textrm{lead}}\in\mathbb{R}^{n_{x,\textrm{lead}}}$ that assumes that
the lead car keeps its acceleration for $t_{\textrm{acc}}$ seconds and then
has an acceleration of $0$. With the forward prediction of the lead vehicle
$\hat{X}^{\textrm{lead}}$, a minimum distance $d_{\textrm{min}}$ and a
reaction time $t_{\textrm{brake}}$, the safety distance to the leading vehicle
is assured by the constraint

$h_{k}^{\textrm{dist}}(x_{k},\hat{x}_{k}^{\textrm{lead}})\coloneqq\max\left(\frac{v_{k}^{2}-(\hat{v}_{k}^{\textrm{lead}})^{2}}{2|a_{\textrm{max}}|}+v_{k}\;t_{\textrm{brake}},\;d_{\textrm{min}}\right)-\left(\hat{s}_{k}^{\textrm{lead}}-s_{k}\right).$
(10)

Additionally, to respect speed limits we use the following constraint

$h_{k}^{v}(x_{k},v^{\textrm{max}})\coloneqq
v_{k}-v^{\textrm{max}}(s_{k}),\quad\text{with}\;v^{\textrm{max}}(s_{k})\coloneqq\begin{cases}v_{\textrm{max,
1}}&\quad\text{if }s_{k}<s_{\textrm{change}}\\\ v_{\textrm{max,
2}}&\quad\text{else}\end{cases}$ (11)

where $v^{\textrm{max}}$ is a position-based function encoding the current
speed limit, with $v_{\textrm{max, 1}}$ being the current speed limit,
$v_{\textrm{max, 2}}$ the next speed limit, and $s_{\textrm{change}}$ the
position where the speed limit changes. Note, that in the context of (2.1) the
parameter $p_{k}$ is given by the forward prediction
$\hat{x}^{\textrm{lead}}_{k}$ and the constant speed limit parameters
$v_{\textrm{max},1}$, $v_{\textrm{max},2}$, $s_{\textrm{change}}$.
Additionally, we introduced a quadratic slack variable $\zeta_{\textrm{dist}}$
to loosen the distance constraint of (10) with a weight
$w_{\zeta,\textrm{dist}}>0$ and a slack variable $\zeta_{a_{N}}$ for the
terminal acceleration constraint with a weight $w_{\zeta,\textrm{a}_{N}}>0$.
The exact values of the different parameters and cost weights can be seen in
table 1.

###  4.2 CommonRoad Benchmark

To test the closed-loop performance, we benchmarked our learned models with
the composable benchmark simulator CommonRoad [Althoff et al.(2017)Althoff,
Koschi, and Manzinger] on real and synthetic scenarios. The real scenarios
were created from the HighD [Krajewski et al.(2018)Krajewski, Bock, Kloeker,
and Eckstein] dataset, a collection of vehicle trajectories recorded on German
highways, and then transformed into the CommonRoad format. As the real
scenarios do not cover the whole space of potential driving maneuvers we
additionally created synthetic scenarios: (1) Breaking, involving uniformly
sampled lead cars with decreasing acceleration, (2) Synthetic speed limit
changes, which include uniformly sampled lead cars with either decreasing or
increasing speed limit changes and (3) Synthetic cut-ins, where cut-ins in
front of the ego car are generated with velocities uniformly sampled. The
trajectories of the lead car were generated via the intelligent driver model
(IDM) [Treiber et al.(2000)Treiber, Hennecke, and Helbing]. The average
scenario duration is around 6.5 seconds and we tested on $1540$ scenarios. In
figure 3 we can further see that the predicted open-loop trajectory takes into
account constraints such as the terminal acceleration constraint $a_{N}=0$.

Figure 3: The closed-loop control plots show a trajectory generated in the
CommonRoad simulator. The open-loop plan shows a predicted trajectory with
sampled inputs from the test set $\mathcal{D_{\mathrm{test}}}$. For the open-
loop behavior of BC, we only report the first control, marked with an orange
cross. The position constraints are either the forward prediction (open-loop)
or the driven trajectory of the leading vehicle (closed-loop). In the closed-
loop control example, the dashed line jumps because of a cut-in, thus a
different car becomes the leading vehicle.

###  4.3 Experimental Setup

#### Dataset

For solving the MPC formulation we used the NLP solver acados [Verschueren et
al.(2021)Verschueren, Frison, Kouzoupis, Frey, van Duijkeren, Zanelli,
Novoselnik, Albin, Quirynen, and Diehl] with HPIPM [Frison and Diehl(2020)].
For creating training data, we uniformly sampled the initial state
$\bar{x}_{0}$ from the box constraints of the longitudinal planner and
filtered out scenarios where a crash was inevitable. With a probability of
$1/3$, we sample a speed limit change and with a probability of $1/3$ we
sample that the car makes a cut-in, leading to a different forward prediction.
We created a training dataset $\mathcal{D}$ with $10^{5}$ samples, a
validation dataset $\mathcal{D_{\mathrm{valid}}}$ with $33333$ samples and a
test set $\mathcal{D_{\mathrm{test}}}$ with $33333$ samples.

#### Metrics

We use the following metrics averaged over three training seeds: 1) _Plan MSE_
: Is the mean squared error between the predicted trajectory and the ground
truth $\hat{X}^{\theta}$ and $X^{\textrm{MPC}}$ on the test set
$\mathcal{D_{\mathrm{test}}}$. 2) _Policy MSE_ : Is the mean squared error
between $\hat{u}_{0}^{\theta}$ and $u_{0}^{\textrm{MPC}}$ on the test set
$\mathcal{D_{\mathrm{test}}}$. 3) _Planning Time_ : We randomly sample $1000$
inputs and estimate their inference times. As there is some additive CPU noise
e.g. from background processes, we repeat the inference benchmark $20$ times
for each input, taking the minimum. Across all samples, we take the $95\%$
quantile to estimate the worst-case runtime. 4) _Policy Time_ : Similar to the
planning time estimation, we estimate the inference time, but only for the
first control. 5) _Collisions_ : Number of collisions encountered during
simulation. We checked for collisions with the CommonRoad feasibility checker
6) _Avg. $\lVert s\rVert_{2}$_, _Avg. $\lVert v\rVert_{2}$_ and _Avg. $\lVert
a\rVert_{2}$_: Are the average $\ell^{2}$-distances between the driven
trajectories of the closed-loop MPC policy and the learned policy, in terms of
the position $s$, velocity $v$, and acceleration $a$.

#### Hyperparameters

For each neural network architecture and loss combination, we did a
hyperparameter search using the Optuna library [Akiba et al.(2019)Akiba, Sano,
Yanase, Ohta, and Koyama]. The search space was defined by $[1\times
10^{-5},1\times 10^{-2}]$ for the learning rate, $\\{1,2,3\\}$ for the number
of hidden layers in $\pi_{\theta}$, $\\{32,64,128,256,512\\}$ for the width of
the hidden layers and $\\{8,16,32,64\\}$ for the batch size. For the policy
part of the PlanNetX architectures and the baselines, we found that a width of
$512$ and a depth of $3$ worked the best. Thus, this was also the model that
we used for the pruning and quantization experiments. For the PlanNetXEnc
Experiments we used three self-attention Transformer layers with embedding
dimension $32$ and $4$ attention heads. For the hyperparameter search, we used
$100$ training epochs, whereas for the full training, we used $300$ epochs. In
all experiments, we used the Adam optimizer [Kingma and Ba(2015)]. For the
weight matrix $W$ in (4) and (5), we always used the identity matrix, as this
gave us the best performance. For the discount factor $\gamma$ we always used
a value of $0.98$. For measuring the inference time, we used an AMD Ryzen 9
5950X @ 3.40 GHz CPU.

###  4.4 Results

#### Open-loop performance

Method | Policy MSE | Plan MSE
---|---|---
BC | 0.66 | -
PlanNetX | 107.68 | 0.01
PlanNetXEnc | 96.72 | 0.02
PlanNetX with $\mathcal{L}^{\textrm{u}}$ | 9.37 | 0.84
Table 2: Open-loop performance

For the open-loop planning performance, we measured the _Policy MSE_ and _Plan
MSE_ on the test dataset $\mathcal{D_{\mathrm{test}}}$. From table 2, we can
see that the proposed methods PlanNetX and PlanNetXEnc can approximate the MPC
plan policy very well especially when compared to the PlanNetX network with
the $\mathcal{L}^{u}$ loss. Interestingly, the _Policy MSE_ of the PlanNetX
methods that use the state-trajectory loss is significantly higher than of BC,
indicating that sometimes the PlanNetX deviates from the MPC plan to find its
control sequence that leads to better tracking in state space. We noticed that
this was especially true when the MPC planner showed bang-bang behavior for
the control sequence.

#### Closed-loop performance

When applying the approach PlanNetX in simulation for closed-loop control, we
find that the driven trajectories are closer to the trajectories of the MPC
policy than when just using the BC baseline, see figure 3 and table 2. The
PlanNetXEnc achieves a similar performance then PlanNetX, whereas the control
trajectory loss from (4) leads to significantly worse results, indicating that
learning the full control trajectory $U^{\star}$ does not necessarily help for
learning the first control $u^{\star}_{0}$. For the inference time of the
control policy, we achieve a significant speed-up over the MPC solver. Even
though, for the planning time, the speed up is not as significant, we were
able to reduce the planning time using structured $\ell^{1}$-pruning, without
losing significantly in closed-loop performance. For quantization, we found a
reduced inference time at the expense of the accuracy of our predictions.
Since inference is done over multiple time steps, errors accumulate, leading
to a large decrease in accuracy, making this architecture unsuitable for
static quantization after training. Further experiments could be done with
quantized aware training, which introduces quantization during training to get
better accuracy. To further speed up the inference time, we could also predict
only a couple of steps of the planned trajectory.

Method |  Avg. $\lVert\Delta s\rVert_{2}$ [m] |  Avg. $\lVert\Delta v\rVert_{2}$ [m/s] |  Avg. $\lVert\Delta a\rVert_{2}$ [m/$\textrm{s}^{2}$] | # Collisions | Policy Time [ms] | Planning Time [ms]
---|---|---|---|---|---|---
MPC | - | - | - | 0 | 6.44 | 6.44
MPC warm start | - | - | - | - | 5.08 | 5.08
BC | 0.65 | 0.17 | 0.05 | 0 | 0.06 | -
PlanNetX | 0.24 | 0.05 | 0.04 | 0 | 0.10 | 2.62
PlanNetXEnc | 0.29 | 0.05 | 0.03 | 0 | 0.31 | 2.79
PlanNetX with $\mathcal{L}^{u}$ | 2.00 | 0.56 | 0.13 | 0 | 0.10 | 2.61
PlanNetX quant. | 1.19 | 0.39 | 0.19 | 0 | 0.07 | 1.94
PlanNetX $\ell^{1}$-str | 0.27 | 0.11 | 0.06 | 0 | 0.06 | 1.83
Table 3: The closed-loop performance for the proposed PlanNetX and PlanNetXEnc
approaches and some further ablations, where we tested the control trajectory
loss $\mathcal{L}^{u}$, a pruned and quantized network as well as warm
starting the MPC.

#### When to use the encoder?

Figure 4: Trajectory predictions for varying forward pred.
$\hat{X}^{\textrm{lead}}$.

We further tested the capabilities of the encoder architecture to adapt to
varying forward predictions of the lead car $\hat{X}^{\textrm{lead}}$. In such
a scenario, the PlanNetX misses information about the whole parameter matrix
$P$. For this, we create a dataset, where the forward predictions are
randomized by incorporating a cut-in in the forward prediction, such that the
forward prediction contains a jump at a random time point. The results can be
seen in figure 4. Unsurprisingly, only the encoder architecture is able to fit
the MPC plan in that scenario. We defer a more detailed analysis for future
research, as in our example, there was no need for varying forward
predictions.

##  5 Conclusion

We found that our proposed PlanNetX framework can imitate the trajectories of
the longitudinal MPC planner to high accuracy while achieving a better closed-
loop control in simulation than BC and having a significantly lower inference
time than the MPC planner. Further, directions could be to extend the method
to MPC with non-linear dynamics where stability issues due to vanishing or
exploding gradients might arise. Here, reducing the roll-out horizon, or using
an MPC-based value function [Ghezzi et al.(2023)Ghezzi, Hoffman, Frey,
Boedecker, and Diehl] might further stabilize the training. Another
interesting direction would be to investigate the generalization capabilities
of the PlanNetXEnc architecture for more MPC formulations with a potentially
more complex parameterization of the MPC.

The authors thank Rudolf Reiter for providing helpful feedback throughout this
work. This project has been supported by the European Union’s Horizon 2020
research and innovation programme, Marie Skłodowska-Curie grant agreement
953348, ELO-X, and by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) with grant number 428605208.

## References

  * [Akiba et al.(2019)Akiba, Sano, Yanase, Ohta, and Koyama] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A Next-generation Hyperparameter Optimization Framework.  In Ankur Teredesai, Vipin Kumar, Ying Li, Rómer Rosales, Evimaria Terzi, and George Karypis, editors, _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD 2019, Anchorage, AK, USA, August 4-8, 2019_, pages 2623–2631. ACM, 2019.  10.1145/3292500.3330701.
  * [Althoff et al.(2017)Althoff, Koschi, and Manzinger] Matthias Althoff, Markus Koschi, and Stefanie Manzinger.  CommonRoad: Composable benchmarks for motion planning on roads.  In _2017 IEEE Intelligent Vehicles Symposium (IV)_ , pages 719–726, June 2017.  10.1109/IVS.2017.7995802.
  * [Bemporad et al.(1999)Bemporad, Borrelli, and Morari] A. Bemporad, F. Borrelli, and M. Morari.  The explicit solution of constrained LP-Based receding horizon control.  In _Proceedings of the IEEE Conference on Decision and Control (CDC)_ , Sydney, Australia, 1999.
  * [Carius et al.(2020)Carius, Farshidian, and Hutter] Jan Carius, Farbod Farshidian, and Marco Hutter.  MPC-Net: A First Principles Guided Policy Search.  _IEEE Robotics and Automation Letters_ , 5(2):2897–2904, April 2020.  ISSN 2377-3766, 2377-3774.  10.1109/LRA.2020.2974653.
  * [Cosner et al.(2022)Cosner, Yue, and Ames] Ryan K. Cosner, Yisong Yue, and Aaron D. Ames.  End-to-End Imitation Learning with Safety Guarantees using Control Barrier Functions.  In _2022 IEEE 61st Conference on Decision and Control (CDC)_ , pages 5316–5322, Cancun, Mexico, December 2022. IEEE.  ISBN 978-1-66546-761-2.  10.1109/CDC51059.2022.9993193.
  * [Devlin et al.(2019)Devlin, Chang, Lee, and Toutanova] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  In Jill Burstein, Christy Doran, and Thamar Solorio, editors, _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT 2019, Minneapolis, MN, USA, June 2-7, 2019, Volume 1 (Long and Short Papers)_ , pages 4171–4186. Association for Computational Linguistics, 2019.  10.18653/V1/N19-1423.
  * [Dosovitskiy et al.(2021)Dosovitskiy, Beyer, Kolesnikov, Weissenborn, Zhai, Unterthiner, Dehghani, Minderer, Heigold, Gelly, Uszkoreit, and Houlsby] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale.  In _9th International Conference for Learning Representations_. International Conference for Learning Representations, 2021.
  * [Frison and Diehl(2020)] G. Frison and M. Diehl.  HPIPM: A high-performance quadratic programming framework for model predictive control.  In _Proceedings of the IFAC World Congress_ , Berlin, Germany, July 2020.
  * [Ghezzi et al.(2023)Ghezzi, Hoffman, Frey, Boedecker, and Diehl] Andrea Ghezzi, Jasper Hoffman, Jonathan Frey, Joschka Boedecker, and Moritz Diehl.  Imitation Learning from Nonlinear MPC via the Exact Q-Loss and its Gauss-Newton Approximation.  In _Proceedings of the IEEE Conference on Decision and Control (CDC)_. arXiv, 2023.  10.48550/arXiv.2304.01782.
  * [Gholami et al.(2021)Gholami, Kim, Dong, Yao, Mahoney, and Keutzer] Amir Gholami, Sehoon Kim, Zhen Dong, Zhewei Yao, Michael W. Mahoney, and Kurt Keutzer.  A Survey of Quantization Methods for Efficient Neural Network Inference.  In _Low-Power Computer Vision_ , pages 291–326. Chapman and Hall/CRC, June 2021.
  * [Gutjahr et al.(2017)Gutjahr, Gröll, and Werling] Benjamin Gutjahr, Lutz Gröll, and Moritz Werling.  Lateral Vehicle Trajectory Optimization Using Constrained Linear Time-Varying MPC.  _IEEE Transactions on Intelligent Transportation Systems_ , 18(6):1586–1595, June 2017.  ISSN 1558-0016.  10.1109/TITS.2016.2614705.
  * [Johansen and Grancharova(2003)] T.A. Johansen and A. Grancharova.  Approximate explicit constrained linear model predictive control via orthogonal search tree.  _IEEE Trans. Automatic Control_ , 48:810–815, 2003.
  * [Kingma and Ba(2015)] Diederik P. Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization.  In _3rd International Conference for Learning Representations_ , volume 3rd. International Conference for Learning Representations, 2015.  10.48550/arXiv.1412.6980.
  * [Krajewski et al.(2018)Krajewski, Bock, Kloeker, and Eckstein] Robert Krajewski, Julian Bock, Laurent Kloeker, and Lutz Eckstein.  The highD Dataset: A Drone Dataset of Naturalistic Vehicle Trajectories on German Highways for Validation of Highly Automated Driving Systems.  In _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_ , pages 2118–2125, 2018.  10.1109/ITSC.2018.8569552.
  * [Levine and Koltun(2013a)] Sergey Levine and Vladlen Koltun.  Guided Policy Search.  In _Proceedings of the 30th International Conference on Machine Learning_ , pages 1–9. PMLR, May 2013a.
  * [Levine and Koltun(2013b)] Sergey Levine and Vladlen Koltun.  Variational Policy Search via Trajectory Optimization.  In _Advances in Neural Information Processing Systems_ , volume 26. Curran Associates, Inc., 2013b.
  * [Li et al.(2017)Li, Kadav, Durdanovic, Samet, and Graf] Hao Li, Asim Kadav, Igor Durdanovic, Hanan Samet, and Hans Peter Graf.  Pruning Filters for Efficient ConvNets.  In _5th International Conference on Learning Representations_. International Conference on Learning Representations, 2017.
  * [Li and Todorov(2004)] Weiwei Li and Emanuel Todorov.  Iterative linear quadratic regulator design for nonlinear biological movement systems.  In _First International Conference on Informatics in Control, Automation and Robotics_ , volume 2, pages 222–229. SciTePress, 2004.
  * [Mamedov et al.(2022)Mamedov, Reiter, Diehl, and Swevers] Shamil Mamedov, Rudolf Reiter, Moritz Diehl, and Jan Swevers.  Safe Imitation Learning of Nonlinear Model Predictive Control for Flexible Robots.  _CoRR_ , abs/2212.02941, 2022.  10.48550/ARXIV.2212.02941.
  * [Mordatch and Todorov(2014)] Igor Mordatch and Emo Todorov.  Combining the benefits of function approximation and trajectory optimization.  In _Robotics: Science and Systems X_. Robotics: Science and Systems Foundation, July 2014.  ISBN 978-0-9923747-0-9.  10.15607/RSS.2014.X.052.
  * [Pek and Althoff(2021)] Christian Pek and Matthias Althoff.  Fail-Safe Motion Planning for Online Verification of Autonomous Vehicles Using Convex Optimization.  _IEEE Transactions on Robotics_ , 37(3):798–814, June 2021.  ISSN 1941-0468.  10.1109/TRO.2020.3036624.
  * [Reske et al.(2021)Reske, Carius, Ma, Farshidian, and Hutter] Alexander Reske, Jan Carius, Yuntao Ma, Farbod Farshidian, and Marco Hutter.  Imitation Learning from MPC for Quadrupedal Multi-Gait Control.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 5014–5020, May 2021.  10.1109/ICRA48506.2021.9561444.
  * [Schwan et al.(2023)Schwan, Jones, and Kuhn] Roland Schwan, Colin N. Jones, and Daniel Kuhn.  Stability Verification of Neural Network Controllers Using Mixed-Integer Programming.  _IEEE Transactions on Automatic Control_ , pages 1–16, 2023.  ISSN 1558-2523.  10.1109/TAC.2023.3283213.
  * [Treiber et al.(2000)Treiber, Hennecke, and Helbing] Martin Treiber, Ansgar Hennecke, and Dirk Helbing.  Congested Traffic States in Empirical Observations and Microscopic Simulations.  _Physical Review E_ , 62(2):1805–1824, August 2000.  ISSN 1063-651X, 1095-3787.  10.1103/PhysRevE.62.1805.
  * [Vaswani et al.(2017)Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser, and Polosukhin] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is All you Need.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * [Verschueren et al.(2021)Verschueren, Frison, Kouzoupis, Frey, van Duijkeren, Zanelli, Novoselnik, Albin, Quirynen, and Diehl] Robin Verschueren, Gianluca Frison, Dimitris Kouzoupis, Jonathan Frey, Niels van Duijkeren, Andrea Zanelli, Branimir Novoselnik, Thivaharan Albin, Rien Quirynen, and Moritz Diehl.  Acados – a modular open-source framework for fast embedded optimal control.  _Mathematical Programming Computation_ , October 2021.  ISSN 1867-2957.  10.1007/s12532-021-00208-8.
  * [Wabersich and Zeilinger(2021)] Kim Peter Wabersich and Melanie N. Zeilinger.  A predictive safety filter for learning-based control of constrained nonlinear dynamical systems.  _Automatica_ , 129:109597, July 2021.  ISSN 0005-1098.  10.1016/j.automatica.2021.109597.
  * [Xie et al.(2021)Xie, Li, Van Wyk, Dellaert, Boots, and Ratliff] Mandy Xie, Anqi Li, Karl Van Wyk, Frank Dellaert, Byron Boots, and Nathan Ratliff.  Imitation Learning via Simultaneous Optimization of Policies and Auxiliary Trajectories.  _CoRR_ , abs/2105.03019, June 2021.


A Sound Up-to-, Bisimilarity for PCTL

Massimo Bartoletti1, Maurizio Murgia2,3, Roberto Zunino2

Università degli Studi di Cagliari, Cagliari, Italy
Università degli Studi di Trento, Trento, Italy
Università degli Studi di Bologna, Bologna, Italy

We tackle the problem of establishing the soundness of approximate
bisimilarity with respect to PCTL and its relaxed semantics.

To this purpose, we consider a notion of bisimilarity similar to the
one introduced by Desharnais, Laviolette, and Tracol, which is
parametric with respect to an approximation error $\distD$, and to the
depth $\stepsN$ of the observation along traces.

Essentially, our soundness theorem establishes that, when a state
$\stateQ$ satisfies a given formula up-to error $\distD$ and steps
$\stepsN$, and $\stateQ$ is bisimilar to $\stateQi$ up-to error
$\distDi$ and enough steps, we prove that $\stateQi$ also satisfies
the formula up-to a suitable error $\distDii$ and steps $\stepsN$.

The new error $\distDii$ is computed from $\distD,\distDi$
and the formula, and only depends linearly on $\stepsN$.

We provide a detailed overview of our soundness proof, which exploits
the min-cut/max-flow theorem on graphs.

[1]
Bian, G., Abate, A.: On the relationship between bisimulation and trace
equivalence in an approximate probabilistic context. In: Esparza, J.,
Murawski, A.S. (eds.) Foundations of Software Science and Computation
Structures (FOSSACS). LNCS, vol. 10203, pp. 321–337 (2017).

[2]
Castiglioni, V., Gebler, D., Tini, S.: Logical characterization of bisimulation
metrics. In: Workshop on Quantitative Aspects of Programming Languages and
Systems (QAPL). EPTCS, vol. 227, pp. 44–62 (2016).

[3]
Dantzig, G.B., Fulkerson, D.R.: On the Max Flow Min Cut Theorem of Networks.
RAND Corporation, Santa Monica, CA (1955)

[4]
Desharnais, J., Edalat, A., Panangaden, P.: Bisimulation for labelled markov
processes. Inf. Comput.  179(2),  163–193 (2002).

[5]
Desharnais, J., Gupta, V., Jagadeesan, R., Panangaden, P.: Metrics for labeled
markov systems. In: CONCUR. LNCS, vol. 1664, pp. 258–273. Springer (1999)

[6]
Desharnais, J., Gupta, V., Jagadeesan, R., Panangaden, P.: Weak bisimulation is
sound and complete for pCTL\({}^{\mbox{*}}\). Inf. Comput.
208(2),  203–219 (2010). 10.1016/j.ic.2009.11.002

[7]
Desharnais, J., Laviolette, F., Tracol, M.: Approximate analysis of
probabilistic processes: Logic, simulation and games. In: Quantitative
Evaluaiton of Systems (QEST). pp. 264–273. IEEE Computer Society (2008).

[8]
D'Innocenzo, A., Abate, A., Katoen, J.: Robust PCTL model checking. In: Dang,
T., Mitchell, I.M. (eds.) Hybrid Systems: Computation and Control (HSCC).
pp. 275–286. ACM (2012). 10.1145/2185632.2185673

[9]
Hansson, H., Jonsson, B.: A logic for reasoning about time and reliability.
Formal Aspects Comput.  6(5),  512–535 (1994).

[10]
Larsen, K.G., Skou, A.: Bisimulation through probabilistic testing. Inf.
Comput.  94(1),  1–28 (1991). 10.1016/0890-5401(91)90030-6