# Model Predictive Path Integral Control Framework for Partially Observable
Navigation: A Quadrotor Case Study

Ihab S. Mohamed1 and Guillaume Allibert2 and Philippe Martinet1 This research
was supported by the ANR CLARA project (ANR-18-CE33-0004).1Ihab S. Mohamed and
Philippe Martinet are with the Université Côte d’Azur, Inria, France,
{ihab.mohamed<EMAIL_ADDRESS>Allibert is with the
Université Côte d’Azur, CNRS, I3S, France<EMAIL_ADDRESS>
###### Abstract

Recently, Model Predictive Path Integral (MPPI) control algorithm has been
extensively applied to autonomous navigation tasks, where the cost map is
mostly assumed to be known and the 2D navigation tasks are only performed. In
this paper, we propose a generic MPPI control framework that can be used for
2D or 3D autonomous navigation tasks in either fully or partially observable
environments, which are the most prevalent in robotics applications. This
framework exploits directly the 3D-voxel grid acquired from an on-board
sensing system for performing collision-free navigation. We test the
framework, in realistic RotorS-based simulation, on goal-oriented quadrotor
navigation tasks in a cluttered environment, for both fully and partially
observable scenarios. Preliminary results demonstrate that the proposed
framework works perfectly, under partial observability, in 2D and 3D cluttered
environments.

## Multimedia Material

The supplementary video attached to this work is available at:
https://urlz.fr/cs2L

##  I INTRODUCTION

Having a safe and reliable system for autonomous navigation of robotic systems
such as Unmanned Aerial Vehicles (UAVs) is a highly challenging and partially-
solved problem for robotics communities, especially for cluttered and GPS-
denied environments such as dense forests, crowded offices, corridors, and
warehouses. Such a problem is very important for solving many complex
applications, such as surveillance, search-and-rescue, and environmental
mapping. To do so, UAVs should be able to navigate with complete autonomy
while avoiding all kinds of obstacles in real-time. To this end, they must be
able to (i) perceive their environment, (ii) understand the situation they are
in, and (iii) react appropriately.

Obviously enough, this problem has been already addressed in the literature,
particularly those works related to dynamics and control, motion planning, and
trajectory generation in unstructured environments with obstacles [1, 2, 3, 4,
5]. Moreover, the applications of the path-integral control theory have
recently become more prevalent. One of the most noteworthy works is Williams’s
iterative path integral method, namely MPPI control framework [6]. In this
method, the control sequence is iteratively updated to obtain the optimal
solution on the basis of importance sampling of trajectories. In [7], authors
derived a different iterative method in which the control- and noise-affine
dynamics constraints, on the original MPPI framework, are eliminated. This
framework is mainly based on the information-theoretic interpretation of
optimal control using KL-divergence and free energy, while it was previously
based on the linearization of Hamilton-Jacob Bellman (HJB) equation and
application of Feynman-Kac lemma. Although different methods are adopted to
derive the MPPI framework, they are practically equivalent and theoretically
related111In the sense that the method in [7] can exactly recover that in [6]
if dynamics is considered to be affine in control. In other words, the
iterative method in [7] can be seen as the generalization of the latter
method.. An extension to Williams’s information-theoretic-based approach is
presented in [8], where a learned model is used to generate informed sampling
distributions.

The attractive features of MPPI controller, over alternative methods, can be
summarized as: (i) a derivative-free optimization method, i.e. no need for
derivative information to find the optimal solution; (ii) no need for
approximating the system dynamics and cost functions with linear and quadratic
forms, i.e., non-linear and non-convex functions can be naturally employed,
even that dynamics and cost models can be easily represented using neural
networks; (iii) planning and execution steps are combined into a single step,
providing an elegant control framework for autonomous vehicles. However, one
drawback of MPPI is that its convergence rate is empirically slow, which has
been addressed in [9].

In the context of autonomous navigation, it is observed that the MPPI
controller has been mainly applied to the tasks of aggressive driving and UAVs
navigation in 2D cluttered environments. To do so, MPPI requires a cost map,
as an environment representation, to drive the autonomous vehicle. Concerning
autonomous driving, the cost map is obtained either off-line [7] or from an
on-board monocular camera using deep learning approaches [10, 11]. Regarding
UAV navigation in cluttered environments, the obstacle (i.e., cost) map is
assumed to be available, and only static 2D floor-maps are used. Conversely,
in practice, the real environments are often partially observable, with
dynamic obstacles. Moreover, it is noteworthy that only 2D navigation tasks
are performed so far, which limits the applicability of the control framework.
For this reason, this paper focuses on MPPI for 2D and 3D navigation tasks in
a previously unseen and dynamic environment. In particular, the main
contributions of our work can be summarized as follows:

  1. 1.

We propose a generic MPPI framework for autonomous navigation in cluttered 2D
and 3D environments, which are inherently uncertain and partially observable.
To the best of our knowledge, this point has not been reported in the
literature, which opens up new directions for research.

  2. 2.

We demonstrate this framework on a set of simulated quadrotor navigation tasks
using RotorS simulator and Gazebo [12], assuming that: (i) there is a priori
knowledge about the environment (namely, fully observable case); (ii) there is
no a priori information (namely, partially observable case). In this case, the
robot is building and updating the map, which represents the environment,
online as it goes along. This allows the opportunity to navigate in dynamic
environments.

  3. 3.

To ensure a realistic simulation, our proposed framework is evaluated taking
into account the modelling errors, noisy sensors, and windy environments.

This paper is organized as follows. Section II describes the quadrotor
dynamics model which represents our case study for the framework validation,
whereas in Section III the real-time MPPI control strategy is explained. Our
proposed framework is then described in Section IV and evaluated in Section V.
Finally, concluding remarks are provided in Section VI.

##  II Quadrotor Dynamics Model

Considering a quadrotor vehicle model illustrated in Fig. 1, the dynamics
model can be defined by assigning a fixed inertial frame $\mathcal{W}$ and
body frame $\mathcal{B}$ attached to the vehicle. The origin of the body
frame, $\mathcal{B}$, is located at the center of mass of the quadrotor, where
$x_{B}$ and $y_{B}$ lie in the quadrotor plane defined by the centers of the
four rotors, and $z_{B}$ is perpendicular to this plane and points upward. The
inertial reference frame, $\mathcal{W}$, is defined by $x_{W},y_{W}$, and
$z_{W}$, with $z_{W}$ pointing upward. We assume that the first rotor (i.e.,
along the $+x_{B}$ axis) and the third rotor rotate counter clockwise (namely,
$+1$), whilst the second (i.e., along the $+y_{B}$ axis) and fourth rotors
rotate clockwise (namely, $-1$). Here, the vehicle is only subject to: (i) a
gravitational acceleration $g$ in $-z_{W}$ direction; (ii) the sum of the
forces generated by each rotor, $F=\sum_{i=1}^{4}F_{i}$, acts along the
$+z_{B}$ direction.

Figure 1: Schematic of the considered quadrotor model in conjunction with the
coordinate systems and forces acting on a vehicle frame.

Furthermore, the Euler angles, with $ZXY$ transformation sequence, are used to
model the rotation of the quadrotor in frame $\mathcal{W}$, where the roll
$\phi$, pitch $\theta$, and yaw $\psi$ angles refer to a rotation about
$x_{B}$, $y_{B}$, and $z_{B}$ axis, respectively. The rotation matrix from
$\mathcal{B}$ to $\mathcal{W}$ is accordingly expressed as

$^{\mathcal{W}}\mathcal{R}_{\mathcal{B}}=\left[\begin{array}[]{ccc}{c_{\psi}c_{\theta}-s_{\phi}s_{\psi}s_{\theta}}&{-c_{\phi}s_{\psi}}&{c_{\psi}s_{\theta}+c_{\theta}s_{\phi}s_{\psi}}\\\
{c_{\theta}s_{\psi}+c_{\psi}s_{\phi}s_{\theta}}&{c_{\phi}c_{\psi}}&{s_{\psi}s_{\theta}-c_{\psi}c_{\theta}s_{\phi}}\\\
{-c_{\phi}s_{\theta}}&{s_{\phi}}&{c_{\phi}c_{\theta}}\end{array}\right],$

where $s_{x}=\sin(x)$ and $c_{x}=\cos(x)$ $\forall$
$x\in\\{\phi,\theta,\psi\\}$. The transformation matrix from Euler angular
velocities, $\dot{\Phi}$, to body frame angular velocity, $\Omega$, is given
by

$T=\left[\begin{array}[]{ccc}{c_{\theta}}&{0}&{-c_{\phi}s_{\theta}}\\\
{0}&{1}&{s_{\phi}}\\\
{s_{\theta}}&{0}&{c_{\phi}c_{\theta}}\end{array}\right].$

Given all considerations above and according to the well-described models in
[1] and [13], the dynamics of the position $\xi$, linear velocity $v$,
orientation locally defined by Euler angles $\Phi$, and body angular rates
$\Omega$ can be written as

$\displaystyle\dot{\xi}$ $\displaystyle=v,$ $\displaystyle m\dot{v}$
$\displaystyle=-mge_{3}+^{\mathcal{W}}\mathcal{R}_{\mathcal{B}}Fe_{3},$ (1)
$\displaystyle\dot{\Phi}$ $\displaystyle=T^{-1}\Omega,$ $\displaystyle
J\dot{\Omega}$ $\displaystyle=\Gamma-\Omega\times J\Omega,$

where $m$ is the mass of the quadrotor and $J$ its inertia matrix expressed in
the body frame $\mathcal{B}$, $\xi=[x,y,z]^{T}$, $v=[v_{x},v_{y},v_{z}]^{T}$,
$\Phi=[\phi,\theta,\psi]^{T}$, $\Omega=[p,q,r]^{T}$,
$e_{3}=[0,0,1]^{T}\in\mathbb{R}^{3}$, $F\in\mathbb{R}^{+}$ is the accumulated
force (i.e., thrust) generated by all rotors and constitutes the first control
input to the system,
$\Gamma=[\tau_{x},\tau_{y},\tau_{z}]^{T}=[L(F_{2}-F_{4}),L(F_{3}-F_{1}),(M_{1}-M_{2}+M_{3}-M_{4})]^{T}\in\mathbb{R}^{3}$
is the total torque applied to the vehicle with its components expressed in
$\mathcal{B}$ which represents the second control input, and
$L\in\mathbb{R}^{+}$ is the distance from the center of the vehicle to the
axis of rotation of each rotor. In this paper, the state of the system is
defined as
$\mathbf{x}=[x,y,z,\phi,\theta,\psi,v_{x},v_{y},v_{z},p,q,r]^{T}\in\mathbb{R}^{12}$.
Each rotor produces a vertical force, $F_{i}$, and moment, $M_{i}$, according
to $F_{i}=k_{F}\omega_{i}^{2}$ and $M_{i}=k_{M}\omega_{i}^{2}$, where
$\omega_{i}$ is the angular velocity of the $i^{th}$ rotor, $k_{F}$ is the
rotor force constant, and $k_{M}$ is the rotor moment constant. Accordingly,
the mapping between the control inputs, namely $F$ and $\Gamma$, and the
system’s input, i.e., $\omega_{i}$, in order to control the quadrotor, can be
expressed as

$\left[\begin{array}[]{c}F\\\ \tau_{x}\\\ \tau_{y}\\\
\tau_{z}\end{array}\right]=\left[\begin{array}[]{cccc}{k_{F}}&{k_{F}}&{k_{F}}&{k_{F}}\\\
{0}&{k_{F}L}&{0}&{-k_{F}L}\\\ {-k_{F}L}&{0}&{k_{F}L}&{0}\\\
{k_{M}}&{-k_{M}}&{k_{M}}&{-k_{M}}\end{array}\right]\left[\begin{array}[]{c}{\omega_{1}^{2}}\\\
{\omega_{2}^{2}}\\\ {\omega_{3}^{2}}\\\ {\omega_{4}^{2}}\end{array}\right].$
(2)

##  III MPPI Control Strategy

The MPPI controller is a stochastic Model Predictive Control (MPC) method,
which can be applied to non-linear dynamics and non-convex cost objectives.
So, it is a sampling-based and derivative-free approach. The key idea of MPPI
is to sample thousands of trajectories, based on Monte-Carlo simulation, in
real-time from the system dynamics. Each trajectory is then evaluated
according to a predefined cost function. Consequently, the optimal control
sequence is updated over all trajectories. This can be easily done, in real-
time, by taking advantage of the parallel nature of sampling and using a
Graphics Processing Unit (GPU).

Let assume that the discrete-time stochastic dynamical system has a form of

$\mathbf{x}_{t+1}=f\left(\mathbf{x}_{t},\mathbf{u}_{t}+\delta\mathbf{u}_{t}\right),$
(3)

where $\mathbf{x}_{t}\in\mathbb{R}^{n}$ is the state vector of the system at
time $t$, $\mathbf{u}_{t}\in\mathbb{R}^{m}$ denotes a control input for the
system, and $\delta\mathbf{u}_{t}$ is a zero-mean Gaussian noise vector with a
variance of $\Sigma_{\mathbf{u}}$, i.e.,
$\delta\mathbf{u}_{t}\sim\mathcal{N}(\mathbf{0},\Sigma_{\mathbf{u}})$, which
represents the control input updates. Given a finite time-horizon
$t\in\\{0,1,2,\cdots,T-1\\}$, the objective of the stochastic optimal control
problem is to find a control sequence,
$\mathbf{u}=\left(\mathbf{u}_{0},\mathbf{u}_{1},\ldots\mathbf{u}_{T-1}\right)\in\mathbb{R}^{m\times
T}$, which minimizes the expectations over all generated trajectories taken
with respect to (3), i.e.,
$J=\min_{\mathbf{u}}\mathbb{E}\left[S\left(\tau\right)\right]$, where
$S\left(\tau\right)\in\mathbb{R}$ is the state-dependent cost-to-go of a
trajectory
$\tau=\left\\{\mathbf{x}_{0},\mathbf{u}_{0},\mathbf{x}_{1},\cdots,\mathbf{u}_{T-1},\mathbf{x}_{T}\right\\}$.
Accordingly, the optimal problem can be formulated as

$J=\min_{\mathbf{u}}\mathbb{E}\left[\phi\left(\mathbf{x}_{T}\right)+\sum_{t=0}^{T-1}\left(q\left(\mathbf{x}_{t}\right)+\frac{1}{2}\mathbf{u}_{t}^{T}R\mathbf{u}_{t}\right)\right],$
(4)

where $\phi\left(\mathbf{x}_{T}\right)$, $q\left(\mathbf{x}_{t}\right)$, and
$R\in\mathbb{R}^{m\times m}$ are a final terminal cost, a state-dependent
running cost, and a positive definite control weight matrix, respectively. To
solve this optimization problem, we consider the iterative update law derived
in [6], in which MPPI algorithm updates the control sequence, from $t$ onward,
as

$\mathbf{u}_{t}\leftarrow\mathbf{u}_{t}+\frac{\sum_{k=1}^{K}\exp\left(-(1/\lambda)\tilde{S}\left(\tau_{t,k}\right)\right)\delta\mathbf{u}_{t,k}}{\sum_{k=1}^{K}\exp\left(-(1/\lambda)\tilde{S}\left(\tau_{t,k}\right)\right)},$
(5)

where $K$ is the number of random samples (namely, rollouts),
$\lambda\in\mathbb{R}^{+}$ is a hyper-parameter so-called the inverse
temperature, and
$\tilde{S}\left(\tau_{t,k}\right)=\phi\left(\mathbf{x}_{T}\right)+\sum_{t=0}^{T-1}\tilde{q}\left(\mathbf{x}_{t},\mathbf{u}_{t},\delta\mathbf{u}_{t}\right)$
is the modified cost-to-go of the $k^{th}$ rollout from time $t$ onward.

Algorithm 1 Real-Time MPPI Control Scheme [6]

Given:
$K,T$: Number of rollouts (samples) & timesteps
$\left(\mathbf{u}_{0},\mathbf{u}_{1},\ldots,\mathbf{u}_{T-1}\right)\equiv\mathbf{u}$:
Initial control sequence
$f,\Delta t$: Dynamics & time-step size
$\phi,q,\lambda,\nu,\Sigma_{\mathbf{u}},R$: Cost functions/hyper-parameters
SGF: Savitzky-Galoy (SG) convolutional filter

1:while task not completed do

2: $\mathbf{x}_{0}\leftarrow$ StateEstimator(),
$\mathbf{x}_{0}\in\mathbb{R}^{n}$

3: $\delta\mathbf{u}\leftarrow$ RandomNoiseGenerator(),
$\delta\mathbf{u}\in\mathbb{R}^{K\times T}$

4: $\tilde{S}\left(\tau_{k}\right)\leftarrow$ TrajectoryCostInitializer(),
$\tilde{S}\left(\tau_{k}\right)\in\mathbb{R}^{K}$

5: for $k\leftarrow 0$ to $K-1$ do

6: $\mathbf{x}\leftarrow\mathbf{x}_{0}$

7: for $t\leftarrow 0$ to $T-1$ do

8:
$\mathbf{x}_{t+1}\leftarrow\mathbf{x}_{t}+f\left(\mathbf{x}_{t},\mathbf{u}_{t}+\delta\mathbf{u}_{t,k}\right)\Delta
t$

9:
$\tilde{S}\left(\tau_{t+1,k}\right)\leftarrow\tilde{S}\left(\tau_{t,k}\right)+\tilde{q}$

10: end for

11:
$\tilde{S}\left(\tau_{k}\right)\leftarrow\tilde{S}\left(\tau_{t+1,k}\right)+\phi\left(\mathbf{x}_{T}\right)$,
$\forall t=T-1$

12: end for

13: $\tilde{S}_{\min}\leftarrow\min_{k}[\tilde{S}\left(\tau_{k}\right)]$

14: for $t\leftarrow 0$ to $T-1$ do

15:
$\mathbf{u}_{t}\leftarrow\mathbf{u}_{t}+\frac{\sum_{k=0}^{K-1}\exp\bigl{(}\frac{-1}{\lambda}\left[\tilde{S}\left(\tau_{t,k}\right)-\tilde{S}_{\min}\right]\bigr{)}\delta\mathbf{u}_{t,k}}{\sum_{k=0}^{K-1}\exp\bigl{(}\frac{-1}{\lambda}\left[\tilde{S}\left(\tau_{t,k}\right)-\tilde{S}_{\min}\right]\bigr{)}}$

16: end for

17: $\mathbf{u}\leftarrow\text{SGF}(\mathbf{u})$

18: $\mathbf{u}_{0}\leftarrow$ SendToActuators($\mathbf{u}$)

19: for $t\leftarrow 1$ to $T-1$ do

20: $\mathbf{u}_{t-1}\leftarrow\mathbf{u}_{t}$

21: end for

22: $\mathbf{u}_{T-1}\leftarrow$
ControlSequenceInitializer($\mathbf{u}_{T-1}$)

23: Check for task completion

24:end while

In this work, the modified running cost
$\tilde{q}\left(\mathbf{x},\mathbf{u},\delta\mathbf{u}\right)$ is defined as

$\tilde{q}=q(\mathbf{x})+\frac{1}{2}\mathbf{u}^{T}R\mathbf{u}+\frac{1-\nu^{-1}}{2}\delta\mathbf{u}^{T}R\delta\mathbf{u}+\mathbf{u}^{T}R\delta\mathbf{u},$
(6)

where $\nu\in\mathbb{R}^{+}$ refers to the exploration noise which determines
how aggressively the state-space is explored. It is noteworthy that the low
values of $\nu$ result in the rejection of many sampled trajectories because
their cost is too high, while too large values result in that the controller
produces control inputs with significant chatter.

The real-time control cycle of MPPI algorithm is described in Algorithm 1 with
more detail. At each time-step $\Delta t$, the system current state is
estimated, and a $K\times T$ random control variations are generated on a GPU
using CUDA’s random number generation library (lines $2:3$). Then, based on
the parallel nature of sampling, all trajectory samples are executed
individually in parallel. For each trajectory, the dynamics are predicted
forward and its expected cost is computed (lines $5:12$), bearing in mind that
the cost of each trajectory is zero-initialized (line $4$). The control
sequence is then updated (lines $14:16$), taking into account the minimum
sampled cost $\tilde{S}_{\min}$ (line $13$). Due to the stochastic nature of
the sampling procedure which leads to significant chattering in the resulting
control, the control sequence is then smoothed using a Savitzky-Galoy (SG)
filter (line 17). Finally, the first control is executed (line 18), while the
remaining sequence of length $T-1$ is slid down to be used at next time-step
$\Delta t$ (lines $19:22$).

##  IV Generic MPPI Framework

In this section, we present a generic and elegant MPPI framework, as
illustrated in Fig. 2, in order to not only navigate autonomously in
previously unseen 2D or 3D environments while avoiding collisions with
obstacles but also to explore and map them. Moreover, we describe how the
environment is represented to be used by MPPI for partially and fully
observable navigation tasks. In Fig. 2, we present the block diagram of our
proposed control framework integrated into the Robot Operating System (ROS).
The individual components of our framework are described in detail in the
following sections.

Figure 2: The global architecture of our proposed framework.

###  IV-A Environment Representation

To clarify how MPPI can be used for 2D or 3D navigation in unseen
environments, we assume that the environment is discretized into a 2D or 3D
grid $\mathcal{M}_{g}$ (see Fig. 2), where each cell is labeled as free,
occupied, or unknown, i.e.,
$\mathcal{M}_{g}=\mathcal{M}_{\text{free}}\cup\mathcal{M}_{\text{occ}}\cup\mathcal{M}_{\text{unk}}$.
In practice, this labeling can be acquired from depth sensors using OctoMap
[14]. As the 3D grid is a grid of cubic volumes of equal size called voxels,
we can represent the environment by a set of layers $\ell_{N}$ along $z_{W}$
direction, where $\ell_{N}=(\frac{\ell_{z}}{r})$, $\ell_{z}$ refers to the
real environment’s height, and $r$ is the voxel size. Accordingly, each layer
represents a 2D occupancy grid, producing in-total $\ell_{N}$ 2D grids for a
given environment, as illustrated in Fig. 3. Since this work is mainly
concerned with the control framework, the perception problem is not presently
covered. Thus, the perception is here imitated by the so-called 2D or 3D mask
$\mathcal{F}_{fov}$, which represents the sensor’s field of view (FoV) as the
robot cannot normally perceive the entire environment.

$1^{\text{st}}$ Layer $\ell_{1}$$2^{\text{nd}}$ Layer
$\ell_{2}$$3^{\text{rd}}$ Layer $\ell_{3}$$4^{\text{th}}$ Layer
$\ell_{4}$$N^{\text{th}}$ Layer
$\ell_{N}$$\mathcal{M}_{\text{unk}}\equiv-1$$\mathcal{M}_{\text{occ}}\equiv
1$$\mathcal{M}_{\text{free}}\equiv 0$3D $\mathcal{F}_{fov}$, Dim: $(3\times
3\times 3)$,where robot is placed in $\ell_{3}$,and $\ell_{2}:\ell_{4}$ are
updated.Robotat $(0,0,\ell_{3})$ Figure 3: 3D environment representation of a
generic MPPI.

For the sake of simplicity, we assume that the 3D $\mathcal{F}_{fov}$ has
dimensions $f_{x}\times f_{y}\times f_{z}$ and an orientation $f_{\theta}$,
where $f_{x}$ and $f_{y}$ represent the number of cells in the 2D grid along
$x_{W}$ and $y_{W}$ axes. While $f_{z}$ represents the number of layers along
$z_{W}$, assuming that the robot’s vertical FoV $f_{N}$ is between
$[(\textit{int}(\frac{-f_{z}}{2})+\ell_{n}):(\textit{int}(\frac{f_{z}}{2})+\ell_{n})]$
where the robot is located in layer $\ell_{n}$ (see Fig. 3). Note that the 3D
$\mathcal{F}_{fov}$ is equivalent to $f_{N}$ 2D $\mathcal{F}_{fov}$.

###  IV-B Fully/Partially Observable 2D/3D Navigation Tasks

Let assume that $\mathcal{M}_{g}$ represents the global map representation of
a given environment. This map is initialized using a priori knowledge about
the environment. Let $\mathcal{M}_{mppi}$ be the local map of MPPI, which has
the same size as $\mathcal{M}_{g}$ and each cell is initialized with $-1$
referring to unknown cells. This local map is continuously updated using an
onboard sensing system, as depicted in Fig. 2.

####  IV-B1 Fully Observable Case

In a fully observable case, MPPI is directly provided with a global map
$\mathcal{M}_{g}$ to compute the trajectories’ cost for avoiding obstacles,
taking into account the current robot position obtained from a state
estimator. For a 2D navigation task, a 2D floor grid is sufficient; only the
first layer, $\ell_{1}$, of 3D voxel-grid will be accessed by MPPI. While, for
the 3D case, the 3D voxel-grid is required for performing collision-free
navigation. Clearly, given the robot’s position, particularly its
$z$-component, the corresponding layer, $\ell_{n}$, to $z$-component is used
for evaluating the cost-to-go of each sampled trajectory (see Fig. 3). At the
moment, the main limitation of the fully observable case is that our framework
is only able to handle static environments; this is an important motivation
for defining a partially observable navigation task that is prevalent in
robotics applications.

####  IV-B2 Partially Observable Case

In this case, the robot is located in a previously unseen and dynamic
environment and must navigate to: (i) a predefined goal, (ii) or explore and
map that area. As the environment is unknown and $\mathcal{M}_{g}$ is not
directly accessible by MPPI, the robot must perceive its environment (in our
case, through the predefined mask $\mathcal{F}_{fov}$) and react
appropriately. Here, the MPPI map $\mathcal{M}_{mppi}$, which is initialized
with $-1$, is fed directly into the control framework as a local map for
avoiding obstacles. This map is continuously updated as the robot moves
around. For the sake of clarity, let us consider the map in Fig. 3 as
$\mathcal{M}_{g}$, where the robot is located in $(0,0,\ell_{3})$ and
$\mathcal{F}_{fov}$ has dimensions of $3\times 3\times 3$. As a consequence,
based on the intersection between $\mathcal{M}_{g}$ and $\mathcal{F}_{fov}$,
the layers from $\ell_{2}$ to $\ell_{4}$ in $\mathcal{M}_{mppi}$ are updated,
while other layers have remained constant.

Algorithm 2 Real-Time Generic MPPI Framework

Given:
$\mathcal{M}_{g},\mathcal{M}_{mppi}$: Global and local map of MPPI
$\mathcal{F}_{fov}$: Sensor’s FoV and its parameters

1:while task not completed do

2: $\mathbf{x}_{0}\leftarrow$ StateEstimator(),
$\mathbf{x}_{0}\in\mathbb{R}^{n}$

3: $\mathcal{M}_{mppi}\leftarrow$
MapUpdate($\mathcal{M}_{g}\cap\mathcal{F}_{fov}$)

4: $\mathbf{u}_{0}\leftarrow$ MPPIController($\mathcal{M}_{mppi}$)

5: Check for task completion

6:end while

Generally speaking, the real-time implementation of our proposed generic MPPI
framework for 2D or 3D navigation in cluttered environments is better
described in Algorithm 2, which employs the MPPI control scheme described
above. At each time-step $\Delta t$, the current state is estimated (line 2).
The local map of MPPI $\mathcal{M}_{mppi}$ is then updated accordingly, given
the global map $\mathcal{M}_{g}$ (line 3). As discussed previously,
$\mathcal{M}_{g}$ in conjunction with $\mathcal{F}_{fov}$ is used for updating
the local map $\mathcal{M}_{mppi}$, as the perception modules have not been
considered in the current work. Thus, in practice, the robot must be equipped
with an on-board sensing system, e.g. depth camera or laser scanner, with a
maximum FoV $\mathcal{F}_{fov}$. This currently sensed data is used to obtain
a 2D or 3D occupancy map, which is continuously updated, i.e.,
$\mathcal{M}_{mppi}$. For instance, a quadrotor-based exploration algorithm is
proposed in [15] to build a real-time 3D map. Finally, this map enables the
controller to find the optimal control $\mathbf{u}_{0}$ to be executed,
resulting in collision-free navigation (line 4).

##  V Simulation Details and Results

In this section, we describe how we evaluate our approach in terms of
simulation scenarios and performance metrics. Moreover, we conduct realistic
simulations to evaluate and demonstrate the performance of the proposed
framework.

###  V-A Simulation Setup, Scenarios, and Metrics

####  V-A1 Simulation Setup

In order to evaluate the performance of our proposed MPPI framework in a
previously unseen and cluttered environment, simulation studies have been
performed using RotorS simulator and Gazebo [12]. The parameters of the real
simulated quadrotor (namely, Humminbird quadrotor) are tabulated in Table I.
To ensure a realistic RotorS-based simulation, all navigation tasks are
carried out by (i) using noisy sensors such as GPS and IMU, (ii) adding
external disturbances such as continuous wind with changing speed and
direction, as proposed in [16], and (iii) considering $\pm$10% of modelling
errors in the real values of mass $m$ and inertia $J$ of the prediction model
given in (1). The MPPI controller has a time horizon $t_{p}$ of
$3\text{\,}\mathrm{s}$, a control frequency of $50\text{\,}\mathrm{Hz}$, and
generates 2700 samples each time-step $\Delta t$. The rest of its hyper-
parameters are also listed in Table I, where the $2.5$ value in
$\Sigma_{\mathbf{u}}$ represents the noise in the thrust input $F$ and
$R=\lambda\Sigma_{\mathbf{u}}^{-1}$. For the SG filter, we set the length of
filter window and order of the polynomial function to 51 and 3, respectively.
The real-time execution of MPPI is performed on a GeForce GTX 1080 Ti, where
all algorithms are written in Python and have been implemented using ROS.

TABLE I: Parameters of quadrotor and MPPI Parameter | Value | Parameter | Value
---|---|---|---
$m$ [$\mathrm{kg}$] | $0.716$ |  $t_{p}$ [$\mathrm{s}$] | $3$
$L$ [$\mathrm{m}$] | $0.17$ | $T$ | $150$
$g$ [$\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$] | $9.81$ | $K$ | $2700$
$k_{F}$ [$\mathrm{N}\text{\,}{\mathrm{rpm}}^{-2}$] | $8.55\text{\times}{10}^{-6}$ | $\lambda$ | $0.02$
$k_{M}$ [$\mathrm{N}\text{\,}{\mathrm{rpm}}^{-2}$] | $1.6\text{\times}{10}^{-2}$ | $\nu$ | $1000$
$J$ [$\mathrm{k}\mathrm{g}\;\mathrm{m}\mathrm{{}^{2}}$] | $\operatorname{Diag}\left($7\text{\times}{10}^{-3}$,$7\text{\times}{10}^{-3}$,$12\text{\times}{10}^{-3}$\right)$
$\Sigma_{\mathbf{u}}$ | $\operatorname{Diag}\left(2.5,$5\text{\times}{10}^{-3}$,$5\text{\times}{10}^{-3}$,$5\text{\times}{10}^{-3}$\right)$
$R$ | $\operatorname{Diag}\left($8\text{\times}{10}^{-3}$,4,4,4\right)$

In our RotorS-based simulations, the full-state information,
$\mathbf{x}=[x,y,z,\phi,\theta,\psi,v_{x},$ $v_{y},v_{z},p,q,r]^{T}$, is
directly provided, using an odometry estimator based on an Extended Kalman
Filter (EKF), as an input to the controller. While, since the actual control
signal of the quadrotor is the angular velocity of each rotor $\omega_{i}$,
the controller outputs $F$ and $\Gamma$ are directly converted into
$\omega_{i}$ using (2), to be sent to the quadrotor as Actuators message. In
summary, the closed-loop of MPPI in ROS, at each $\Delta t$, can be summarized
as: (i) MPPI first receives the odometry message and the updated map
$\mathcal{M}_{mppi}$ obtained from the on-board sensor; (ii) the control
action is then computed and published. As mentioned before, the whole process
of our proposed framework is summarized in Fig. 2.

Since we are interested in goal-oriented quadrotor navigation tasks in
cluttered environments, the state-dependent cost function is defined as

$q(\mathbf{x})=(\mathbf{x}-\mathbf{x}^{\text{des}})^{\mathrm{T}}Q(\mathbf{x}-\mathbf{x}^{\text{des}})+${10}^{8}$C_{1}+${10}^{5}$C_{2},$

where:

$\displaystyle C_{1}$
$\displaystyle=\\{\mathbf{x}:\,(z<0)\,\textbf{or}\,(\mathcal{M}_{mppi}[x,y,z]=1)\\},$
$\displaystyle C_{2}$
$\displaystyle=\\{\mathbf{x}:\,(\|v\|>v_{\text{max}})\,\textbf{or}\,\left(\left((\lvert\cos{\phi}\lvert\right)\,\textbf{or}\,\left(\lvert\cos{\theta}\lvert\right)\right)<0.1)$
$\displaystyle\hskip 89.61519pt\textbf{or}\,(z>8.5)\\},$
$\displaystyle\mathbf{x}^{\text{des}}$
$\displaystyle=(x^{\text{des}},y^{\text{des}},z^{\text{des}},0,0,\psi^{\text{des}},\operatorname{zeros}(1,6)),$
$\displaystyle Q$
$\displaystyle=\\!\\!\left\\{\begin{array}[]{l}{\\!\\!\\!\\!\operatorname{Diag}(2.5,2.5,5,1,1,50,\operatorname{zeros}(1,6))},\,\forall
v_{\text{max}}\leq 1.5\frac{\text{m}}{\text{s}},\\\
{\\!\\!\\!\\!\operatorname{Diag}(5,5,15,30,30,50,\operatorname{zeros}(1,6)),\;\text{Otherwise}.}\end{array}\right.$

The first term $C_{1}$ indicates the collision with ground or obstacles, while
$C_{2}$ prevents the quadrotor from (i) going too fast, (ii) using too
aggressive roll and pitch angles, or (iii) colliding with the ceiling, where
$C_{1}$ and $C_{2}$ are boolean variables. $\mathbf{x}^{\text{des}}$ refers to
the desired position to be reached and its orientation, while $Q$ is a
weighing matrix.

####  V-A2 Simulation Scenarios

Two different scenarios are considered for evaluating the proposed framework:
2D scenario and 3D scenario. The 2D scenario refers to a $40\times 40\times
8.5$ $\mathrm{m}$ windy forest of cylindrical obstacles placed in a 2D grid
pattern, where each cylinder has a radius of $0.16\text{\,}\mathrm{m}$ with
equal spacing of $4\text{\,}\mathrm{m}$ apart. The 3D scenario refers to the
same environment described in the 2D scenario, while we have added two
horizontal layers of cylinders at $z=$ $3\text{\,}\mathrm{m}$ and
$6\text{\,}\mathrm{m}$ with the same spacing each, as shown in Fig. 4. The
voxel size $r$ is set to $0.2\text{\,}\mathrm{m}$, to ensure high accuracy and
to meet the reality when a 3D occupancy grid is involved. As a result, the 3D
voxel grid $\mathcal{M}_{g}$ has 43 layers, each layer represents a $40\times
40$ 2D grid.222The cell size of 2D grid is $1\text{\,}\mathrm{m}$ (i.e., 1:1
scale), since all obstacles, in our case, are placed in a 2D grid pattern with
integer numbers. The 3D scenario is used for performing fully and partially
observable 3D navigation tasks, while the former is used for the 2D tasks. For
both scenarios, two different cases are considered: (i) Fully Observable Case
(FOC), in which an a priori knowledge about the environment is used to
initialize the global map $\mathcal{M}_{g}$; in FOC, $\mathcal{M}_{mppi}$ is
exactly $\mathcal{M}_{g}$; and (ii) Partially Observable Case (POC), where it
is assumed that there is no a priori information about the environment; for
this reason, the local map $\mathcal{M}_{mppi}$ is discovered and built online
as the robot moves around.

Figure 4: The 3D cluttered environment used for 3D tasks.

####  V-A3 Performance Metrics

We define two metrics for evaluation. First, the performance of MPPI in POC is
compared to its performance in FOC, considering the latter case as a baseline.
To achieve a fair comparison, in all simulations, the robot navigates to the
same specified goals with a maximum velocity $v_{\text{max}}$ of
$1.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The predefined goals (in
$\mathrm{[}\mathrm{m}\mathrm{]}$) in order are:
$G_{1}=(23,38,1.5),G_{2}=(40,23,4),G_{3}=(22,0,8),G_{4}=(0,22,5)$, then the
UAV will land. While, at each $\Delta t$, the desired yaw angle
$\psi^{\text{des}}$ is updated, letting the front camera points towards the
next goal. For both cases (FOC/POC), we use a number of indicators, as a
second metric, to describe the general performance such as the number of
collisions $N_{\text{col}}$, task completion percentage $t_{\text{comp}}$,
average completion time $t_{\text{av}}$, average flying distance
$d_{\text{av}}$, average flying speed $v_{\text{av}}$, and average energy
consumption $E_{\text{av}}$ (for more details, we refer to [16]). In all
simulations, the task is considered to be terminated if the quadrotor reached
the given goals (i.e., $t_{\text{comp}}=100\%$) or crashed into an obstacle
(i.e., $N_{\text{col}}=1$).

###  V-B Simulation Results

The performance of our proposed framework is validated for both fully and
partially observable 2D/3D quadrotor navigation tasks, considering the
predefined goals. In all partially observable tasks, we set
$\mathcal{F}_{fov}$ to $5\times 5\times 3$ $\mathrm{m}$, while $f_{\theta}$
represents the angle between the current and next goal. To test whether the
quadrotor is able to navigate successfully through the cluttered environment,
we performed 5 trials for both FOC and POC. The reader is invited to watch the
whole simulation results at https://urlz.fr/cs2L.

####  V-B1 2D Navigation Results

Figure 5 shows an example of a final trajectory generated by MPPI in the case
of (i) using a 2D global map $\mathcal{M}_{g}$ for FOC, or (ii) using only a
2D local map $\mathcal{M}_{mppi}$ for POC. The 2D-floor $\mathcal{M}_{g}$
(with grey circles) and its updated $\mathcal{M}_{mppi}$ (with blue circles,
representing the obstacles within the robot’s $\mathcal{F}_{fov}$), for given
goals, are shown in Fig. 6, including the generated trajectories in both cases
and the robot’s FoV $\mathcal{F}_{fov}$. In both 2D FOC and POC, it can be
seen that MPPI is able to safely navigate through the windy obstacle field, in
spite of the presence of modelling errors and measurements noise.

Figure 5: Comparison of 3D trajectories in a 2D scenario. Figure 6: The
$40\text{\,}\mathrm{m}$ obstacle field map, i.e., 2D $\mathcal{M}_{g}$, and
its updated 2D $\mathcal{M}_{mppi}$ map using a $5\times 5$ 2D
$\mathcal{F}_{fov}$. Two lines represent the trajectories generated by MPPI in
both cases, i.e., FOC and POC. Figure 7: Comparison of 3D trajectories in a 3D
scenario.

####  V-B2 3D Navigation Results

Another example of successfully generated trajectories for fully and partially
observable 3D navigation tasks in a 3D scenario environment is depicted in
Fig. 7. In this scenario, since $\mathcal{M}_{mppi}$ consists of 43 2D grids,
it is very difficult to visualize clearly in this paper the updates over all
layers and insure that the robot performs collision-free navigation. So,
first, we show only the $11^{\text{th}}$ layer $\ell_{11}$, which has been
mainly updated while the robot was moving towards $G_{1}$ and while it was
landing after task completion, as an example of how the 3D
$\mathcal{M}_{mppi}$ is updated (see Fig. 8). Second, it is highly recommended
to use the indicators described previously, especially for 3D cases.

Figure 8: The update of the $11^{\text{th}}$ layer $\ell_{11}$ in the 3D
$\mathcal{M}_{mppi}$, showing the generated trajectory in POC.

####  V-B3 Overall Performance

Table II shows the performance statistics for the five trials in 2D and 3D
scenarios, considering the parameters of the controller tabulated in Table I.

TABLE II: General performance of MPPI in 2D and 3D scenarios for both FOC and POC Indicators | 2D Scenario | 3D Scenario
---|---|---
FOC | POC | FOC | POC
$N_{\text{col}}(t_{\text{comp}}[\%])$ | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
$t_{\text{av}}$ [$\mathrm{s}$] | $122.1\pm 1.2$ | $124.5\pm 3.0$ | $128.2\pm 3.0$ | $129.5\pm 1.4$
$d_{\text{av}}$ [$\mathrm{m}$] | $150.3\pm 0.8$ | $151.1\pm 1.9$ | $152.9\pm 1.8$ | $153.1\pm 0.87$
$v_{\text{av}}$ [$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$] | $1.23\pm 0.02$ | $1.22\pm 0.03$ | $1.19\pm 0.02$ | $1.18\pm 0.01$
$E_{\text{av}}$ [$\mathrm{W}\text{\,}\mathrm{h}$] | $6.63\pm 0.07$ | $6.76\pm 0.16$ | $6.96\pm 0.16$ | $7.01\pm 0.07$
$t_{\text{mppi}}$ [$\mathrm{ms}$] | $18.4\pm$0.51$$ | $18.8\pm$0.42$$ | $18.2\pm$0.21$$ | $18.8\pm$0.67$$

In all trials (i.e., 20 trials in total), we can observe that the quadrotor
navigates autonomously (i.e., $N_{\text{col}}=0$) while avoiding obstacles
with an average flying speed $v_{\text{av}}$ of about
$1.2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, which is closer to the
maximum specified velocity, i.e., $v_{\text{max}}=$
$1.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, regardless of the limited
FoV of the sensor in both 2D and 3D scenarios. Note that these 20 trials are
equivalent to more than $3\text{\,}\mathrm{km}$ of autonomous navigation.
Moreover, it can be seen that the average execution time of MPPI per
iteration, $t_{\text{mppi}}$, in the 3D scenario is approximately equal to
that in the 2D scenario (with very low standard deviation values), showing the
superiority of the proposed framework even with 3D-voxel grids and its
applicability to be used for 3D navigation tasks. In fact, this is not
surprising as the 2D- or 3D-voxel grids are stored and directly used by MPPI
on the GPU. For other indicators, we can observe that the performance of MPPI
in POC cases is slightly different than its performance in FOC cases, although
we have set the same parameters for all simulations. The reason behind this
difference is that the control variations $\delta\mathbf{u}$ are generated
randomly, each time-step $\Delta t$, on the GPU.

TABLE III: Effect of changing $T$ and $\nu$ on the behavior of the proposed controller Indicators | 2D Scenario | 3D Scenario
---|---|---
FOC | POC | FOC | POC
Tuning Case 1: $T=75$ (i.e., $t_{p}=$1.5\text{\,}\mathrm{s}$$) and $\nu=1000$
_N_ col (_t_ comp[%]) | 1 (57.4${\scriptstyle\pm}$30) | 1 (57.4${\scriptstyle\pm}$30) | 1 (40.3${\scriptstyle\pm}$6) | 1 (36.3${\scriptstyle\pm}$9)
Tuning Case 2: $T=100$ (i.e., $t_{p}=$2\text{\,}\mathrm{s}$$) and $\nu=1000$
_N_ col (_t_ comp[%]) | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
_d_ av [m] | 51.7${\scriptstyle\pm}$0.79 | 51.7${\scriptstyle\pm}$1.1 | 51.1${\scriptstyle\pm}$0.16 | 52.1${\scriptstyle\pm}$0.68
Tuning Case 3: $T=125$ (i.e., $t_{p}=$2.5\text{\,}\mathrm{s}$$) and $\nu=1000$
_N_ col (_t_ comp[%]) | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
_d_ av [m] | 50.4${\scriptstyle\pm}$0.86 | 50.7${\scriptstyle\pm}$0.25 | 50.2${\scriptstyle\pm}$0.49 | 50.5${\scriptstyle\pm}$0.36
Tuning Case 4: $T=150$ (i.e., $t_{p}=$3\text{\,}\mathrm{s}$$) and $\nu=300$
_N_ col (_t_ comp[%]) | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
_d_ av [m] | 52.4${\scriptstyle\pm}$0.78 | 52.5${\scriptstyle\pm}$0.36 | 52.8${\scriptstyle\pm}$0.74 | 53.0${\scriptstyle\pm}$0.45
Tuning Case 5: $T=150$ and $\nu=500$
_N_ col (_t_ comp[%]) | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
_d_ av [m] | 51.9${\scriptstyle\pm}$0.59 | 51.5${\scriptstyle\pm}$0.85 | 52.1${\scriptstyle\pm}$0.25 | 52.6${\scriptstyle\pm}$0.36
Tuning Case 6: $T=150$ and $\nu=800$
_N_ col (_t_ comp[%]) | 0 (100%) | 0 (100%) | 0 (100%) | 0 (100%)
_d_ av [m] | 50.9${\scriptstyle\pm}$0.25 | 51.0${\scriptstyle\pm}$0.46 | 51.5${\scriptstyle\pm}$0.93 | 51.9${\scriptstyle\pm}$0.33

The number of timesteps $T$ (which is related to the time horizon $t_{p}$ and
control frequency) and the exploration variance $\nu$ (which is related to the
number of sampled trajectories $K$) play an important role in determining the
behavior of MPPI. Therefore, we tested six different tuning cases, where
different values of $T$ and $\nu$ have been considered, as shown in Table III.
In the first three tuning cases, we study the influence of changing $T$ where
$\nu$ remains constant (i.e., $\nu=1000$ as defined before in Table I). While
the effect of changing $\nu$ is studied in the last three cases, where $\nu$
was varied between $300$ and $800$. For the sake of simplicity, each tuning
case was tested by letting the quadrotor navigates only to $G_{1}=(23,38,1.5)$
then lands. In each tuning case, we conducted $3$ trials. As a consequence,
the total trials for both FOC and POC cases are $72$. For Tuning Case 1, where
a short time horizon $t_{p}$ is chosen, it can be clearly noticed that the
MPPI controller is unable to complete the trials at a satisfactory rate, where
the success rate $t_{\text{comp}}$ varies from $36.3\%$ to a maximum of
$57.4\%$ in both 2D FOC and POC (1/3 successful trials). For other tuning
cases (namely, from case 2 to 6), we can observe that the controller performs
perfectly and is able to successfully complete all trials while avoiding
obstacles. We can also notice that as $T$ and $\nu$ increase, the performance
of the controller improves. Clearly, for Tuning Case 3 and 6 where high values
of $T$ and $\nu$ are chosen, we can see that the average flying distance
$d_{\text{av}}$ (for both 2D and 3D scenarios) is the shortest (compared to
other successful cases) with low standard deviation values, which means that
the quadrotor is taking a more direct route leading to the goal. However,
having too long time horizons increase the computational effort dramatically
because each trajectory takes more time to simulate. While, if the exploration
variance $\nu$ is too large, the controller produces control inputs with
significant chatter. We also notice during our simulations that higher values
of control weight matrix $R$ leads to empirically slow convergence of MPPI
towards the goal and fluctuated motion. For this reason, in all experiments,
we set $T$, $\nu$, and $R$ to the values given in Table I, where MPPI performs
very consistently for all given goals.

##  VI Conclusions and future work

Within this work, we proposed an extension to the classical MPPI framework
that enables the robot to navigate autonomously in 2D or 3D environments while
avoiding collisions with obstacles. The key point of our proposed framework is
to provide MPPI with a 2D or 3D grid representing the real world to perform
collision-free navigation. This framework has been successfully tested on
realistic simulations using quadrotor, considering both fully and partially
observable cases. The current simulations illustrate the efficiency and
robustness of the proposed controller for 2D and 3D navigation tasks. Although
the theoretical stability of MPPI has not been addressed and proofed in the
literature, its practical stability can be achieved by setting the MPPI
parameters carefully as we described previously. We will explore the possible
methods that may allow in the future to study the theoretical stability of
MPPI. Our future work will also include the implementation of the framework in
practical applications.

## References

  * [1] D. Mellinger, N. Michael, and V. Kumar, “Trajectory generation and control for precise aggressive maneuvers with quadrotors,” _The International Journal of Robotics Research_ , vol. 31, no. 5, pp. 664–674, 2012.
  * [2] D. González, J. Pérez, V. Milanés, and F. Nashashibi, “A review of motion planning techniques for automated vehicles,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 17, no. 4, pp. 1135–1145, Nov. 2015\.
  * [3] K. Mohta, M. Watterson, Y. Mulgaonkar, S. Liu, C. Qu, A. Makineni, K. Saulnier, K. Sun, A. Zhu, J. Delmerico, _et al._ , “Fast, autonomous flight in GPS-denied and cluttered environments,” _Journal of Field Robotics_ , vol. 35, no. 1, pp. 101–120, 2018.
  * [4] T. Baca, D. Hert, G. Loianno, M. Saska, and V. Kumar, “Model predictive trajectory tracking and collision avoidance for reliable outdoor deployment of unmanned aerial vehicles,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Madrid, Spain, Oct. 2018, pp. 6753–6760.
  * [5] M. Ryll, J. Ware, J. Carter, and N. Roy, “Efficient trajectory planning for high speed flight in unknown environments,” in _International Conference on Robotics and Automation_ , Montreal, Canada, May 2019, pp. 732–738.
  * [6] G. Williams, A. Aldrich, and E. A. Theodorou, “Model predictive path integral control: From theory to parallel computation,” _Journal of Guidance, Control, and Dynamics_ , vol. 40, no. 2, pp. 344–357, 2017.
  * [7] G. Williams, P. Drews, B. Goldfain, J. M. Rehg, and E. A. Theodorou, “Information-theoretic model predictive control: Theory and applications to autonomous driving,” _IEEE Transactions on Robotics_ , vol. 34, no. 6, pp. 1603–1622, 2018.
  * [8] R. Kusumoto, L. Palmieri, M. Spies, A. Csiszar, and K. O. Arras, “Informed information theoretic model predictive control,” in _International Conference on Robotics and Automation_ , Montreal, Canada, May 2019, pp. 2047–2053.
  * [9] M. Okada and T. Taniguchi, “Acceleration of gradient-based path integral method for efficient optimal and inverse optimal control,” in _IEEE International Conference on Robotics and Automation_ , Brisbane, Australia, Sept. 2018, pp. 3013–3020.
  * [10] P. Drews, G. Williams, B. Goldfain, E. A. Theodorou, and J. M. Rehg, “Aggressive deep driving: Combining convolutional neural networks and model predictive control,” in _Conference on Robot Learning_ , 2017, pp. 133–142.
  * [11] A. Buyval, A. Gabdullin, K. Sozykin, and A. Klimchik, “Model predictive path integral control for car driving with autogenerated cost map based on prior map and camera image,” in _2019 IEEE Intelligent Transportation Systems Conference_. IEEE, 2019, pp. 2109–2114.
  * [12] F. Furrer, M. Burri, M. Achtelik, and R. Siegwart, _Robot Operating System (ROS): The Complete Reference (Volume 1)_. Cham: Springer International Publishing, 2016, ch. RotorS—A Modular Gazebo MAV Simulator Framework, pp. 595–625.
  * [13] J.-M. Kai, G. Allibert, M.-D. Hua, and T. Hamel, “Nonlinear feedback control of quadrotors exploiting first-order drag effects,” in _IFAC World Congress_ , Toulouse, France, July 2017, pp. 8189–8195.
  * [14] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “OctoMap: An efficient probabilistic 3D mapping framework based on octrees,” _Autonomous robots_ , vol. 34, no. 3, pp. 189–206, 2013.
  * [15] T. Cieslewski, E. Kaufmann, and D. Scaramuzza, “Rapid exploration with multi-rotors: A frontier selection method for high speed flight,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Vancouver, Canada, Sept. 2017, pp. 2135–2142.
  * [16] T. Tezenas Du Montcel, A. Negre, J.-E. Gomez-Balderas, and N. Marchand, “BOARR: A benchmark for quadrotor obstacle avoidance based on ROS and RotorS,” in _French national conference on ROS (ROSCon Fr)_ , Paris, France, May 2019. [Online]. Available: https://hal.archives-ouvertes.fr/hal-02142571


responsive and sub-complete preferences, Theorem 6 showed that if the couples
additionally belong to one of several types, hrc is solvable in polynomial
time. This provides one type of boundary, and it remains open to identify
other types of couples that could give rise to additional polynomial-time
solvable cases.

Whilst Corollary 16 showed that hrc-dual market is NP-hard even if every
preference list is of length 3 and derived from a master list of single
doctors, couples and hospitals, Corollary 7 showed that hrc-dual market is
solvable in polynomial time for sub-responsive and sub-complete preferences.

Finally, whilst [6] showed that hrc is NP-hard even if every couple’s
preference list is of length at most 2, Theorem 6 showed that the problem is
solvable in polynomial time if every couple’s preference list is of length 1
(since every such couple is trivially of type-b in this case).

An important future direction is to continue to expand our knowledge of the
frontier between tractable and intractable variants of hrc. A key message from
this paper, however, is that hrc  _can_ be tractable in a range of different
scenarios, and perhaps much more so than might have been predicted previously,
given the inherent intractability of the problem in general and in so many
restricted settings.

##  6 Acknowledgements

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## References

  * [1] B. Aldershof and O.M. Carducci.  Stable matching with couples.  Discrete Applied Mathematics, 68:203–207, 1996.
  * [2] P. Berman, M. Karpinski, and A.D. Scott.  Approximation hardness of short symmetric instances of MAX-3SAT.  Electronic Colloquium on Computational Complexity Report, number 49, 2003\.
  * [3] P. Biró, T. Fleiner, and R.W. Irving.  Matching couples with Scarf’s algorithm.  Annals of Mathematics and Artificial Intelligence, 77:303–316, 2016\.
  * [4] P. Biró, R.W. Irving, and I. Schlotter.  Stable matching with couples: an empirical study.  ACM Journal of Experimental Algorithmics, 16, 2011.  Section 1, article 2, 27 pages.
  * [5] P. Biró and F. Klijn.  Matching with couples: a multidisciplinary survey.  International Game Theory Review, 15(2), 2013.  article number 1340008.
  * [6] P. Biró, D.F. Manlove, and I. McBride.  The Hospitals/Residents problem with Couples: Complexity and integer programming models.  In Proceedings of SEA 2014: the 13th International Symposium on Experimental Algorithms, volume 8504 of Lecture Notes in Computer Science, pages 10–21. Springer, 2014.
  * [7] Katarína Cechlárová and Tamás Fleiner.  On a generalization of the stable roommates problem.  ACM Transactions on Algorithms (TALG), 1(1):143–156, 2005.
  * [8] G. Csáji.  On the complexity of stable hypergraph matching, stable multicommodity flow and related problems.  Theoretical Computer Science, 931:1–16, 2022.
  * [9] M. Delorme, S. García, J. Gondzio, J. Kalcsics, D. Manlove, and W. Pettersson.  Stability in the hospitals/residents problem with couples and ties: Mathematical models and computational studies.  Omega, 103:102386, 2021.
  * [10] J. Drummond, A. Perrault, and F. Bacchus.  SAT is an effective and complete method for solving stable matching problems with couples.  In Proceedings of IJCAI ’15: the 24th International Joint Conference on Artificial Intelligence, pages 518–525. AAAI Press, 2015.
  * [11] T. Fleiner.  The stable roommates problem with choice functions.  Algorithmica, 58(1):82–101, 2010.
  * [12] D. Gale and L.S. Shapley.  College admissions and the stability of marriage.  American Mathematical Monthly, 69:9–15, 1962.
  * [13] D. Gale and M. Sotomayor.  Some remarks on the stable matching problem.  Discrete Applied Mathematics, 11:223–232, 1985.
  * [14] D. Gusfield and R.W. Irving.  The Stable Marriage Problem: Structure and Algorithms.  MIT Press, 1989.
  * [15] R.W. Irving.  An efficient algorithm for the “stable roommates” problem.  Journal of Algorithms, 6:577–595, 1985.
  * [16] R.W. Irving, D.F. Manlove, and G. O’Malley.  Stable marriage with ties and bounded length preference lists.  Journal of Discrete Algorithms, 7(2):213–219, 2009.
  * [17] R.W. Irving, D.F. Manlove, and S. Scott.  The stable marriage problem with master preference lists.  Discrete Applied Mathematics, 156(15):2959–2977, 2008.
  * [18] R.W. Irving and S. Scott.  The stable fixtures problem.  Discrete Applied Mathematics, 155:2118–2129, 2007.
  * [19] B. Klaus and F. Klijn.  Stable matchings and preferences of couples.  Journal of Economic Theory, 121(1):75–106, 2005.
  * [20] B. Klaus, F. Klijn, and T. Nakamura.  Corrigendum: Stable matchings and preferences of couples.  Journal of Economic Theory, 144(5):2227–2233, 2009.
  * [21] F. Kojima, P.A. Pathak, and A.E. Roth.  Matching with couples: Stability and incentives in large matching markets.  Quarterly Journal of Economics, 128(4):1585–1632, 2013.
  * [22] D.F. Manlove.  Algorithmics of Matching Under Preferences.  World Scientific, 2013.
  * [23] D.F. Manlove, I. McBride, and J. Trimble.  “Almost-stable” matchings in the Hospitals / Residents problem with Couples.  Constraints, 22:50–72, 2017.
  * [24] D. Marx and I. Schlotter.  Stable assignment with couples: parameterized complexity and local search.  Discrete Optimization, 8:25–40, 2011.
  * [25] E.J. McDermid and D.F. Manlove.  Keeping partners together: Algorithmic results for the hospitals / residents problem with couples.  Journal of Combinatorial Optimization, 19(3):279–303, 2010.
  * [26] C. Ng and D.S. Hirschberg.  Complexity of the stable marriage and stable roommate problems in three dimensions.  Technical Report UCI-ICS 88-28, Department of Information and Computer Science, University of California, Irvine, 1988.
  * [27] C. Ng and D.S. Hirschberg.  Three-dimensional stable matching problems.  SIAM Journal on Discrete Mathematics, 4:245–252, 1991.
  * [28] T. Nguyen and R. Vohra.  Near-feasible stable matchings with couples.  American Economic Review, 108(11):3154–3169, 2018.
  * [29] E. Ronn.  NP-complete stable matching problems.  Journal of Algorithms, 11:285–304, 1990.
  * [30] A.E. Roth.  The evolution of the labor market for medical interns and residents: a case study in game theory.  Journal of Political Economy, 92(6):991–1016, 1984.
  * [31] A.E. Roth.  On the allocation of residents to rural hospitals: a general property of two-sided matching markets.  Econometrica, 54:425–427, 1986.
  * [32] A.E. Roth and E. Peranson.  The redesign of the matching market for American physicians: Some engineering aspects of economic design.  American Economic Review, 89(4):748–780, 1999.
  * [33] A.E. Roth and M.A.O. Sotomayor.  Two-sided matching: a study in game-theoretic modeling and analysis, volume 18 of Econometric Society Monographs.  Cambridge University Press, 1990.
  * [34] H.E. Scarf.  The core of an N-person game.  Econometrica, 35(1):50–69, 1967.
  * [35] J.J.M. Tan.  A necessary and sufficient condition for the existence of a complete stable matching.  Journal of Algorithms, 12(1):154–178, 1991.
  * [36] https://www.nrmp.org (National Resident Matching Program website).
  * [37] https://www.nrmp.org/match-data-analytics/residency-data-reports (NRMP main residency match data and reports).
  * [38] https://www.carms.ca (Canadian Resident Matching Service website).
  * [39] https://tinyurl.com/4wsdudru (Description of the pairing procedure of the SiiMOP-Pairiement platform).
  * [40] https://jrmp.jp (Japan Resident Matching Program website).
  * [41] https://www.nobelprize.org/prizes/economic-sciences/2012/summary (The Sveriges Riksbank Prize in Economic Sciences in Memory of Alfred Nobel 2012).


# Systematic Improvement of Molecular Excited State Calculations by Inclusion
of Nuclear Quantum Motion

Timothy J. H. Hele  Department of Chemistry, University College London, 20,
Gordon Street, London WC1H 0AJ, United Kingdom  Bartomeu Monserrat  Department
of Materials Science and Metallurgy, University of Cambridge, 27 Charles
Babbage Road, Cambridge CB3 0FS, United Kingdom  Cavendish Laboratory,
University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United
Kingdom  Antonios M. Alvertis  Cavendish Laboratory, University of Cambridge,
J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom  e-mail:
<EMAIL_ADDRESS>
## Abstract

Many theoretical studies of molecular excited states aim to provide accurate
solutions to the electronic Schrödinger equation in order to produce energies
that can be compared to experiment. However, the role of nuclear quantum
effects is usually overlooked. Here we provide an intuitive picture for this
effect and find that nuclear quantum fluctuations almost always lower the
energies of excited states from the conventional vertical excitation energy.
We compute the vibration-induced corrections to the exciton energies of the
commonly used Thiel set of molecules by combining TD-DFT with Monte Carlo
sampling techniques based on finite difference methods, and find corrections
as large as $1.1$ eV. We compare our results for the exciton energies to
uncorrected TD-DFT and the accurate reference Thiel values, finding that
incorporating nuclear quantum motion substantially improves agreement with
experiment, at a low computational cost. Accounting for nuclear quantum motion
is therefore established as a critical factor towards the accurate calculation
of exciton energies in molecules.

TOC graphic

The optoelectronic properties of organic molecules are dominated by their low-
lying excited states known as excitons [1]. Exciton energies are critical to
several technologically-relevant processes in these systems, such as singlet
fission and thermally activated delayed fluorescence, which find applications
in photovoltaics and LEDs respectively [2, 3]. It is therefore desirable to
develop methods for the accurate prediction of exciton energies.

Typical methods for the calculation of excited state energies, such as time-
dependent density functional theory (TD-DFT) [4], complete active space self-
consistent ﬁeld (CAS-SCF) [5] and the Bethe-Salpeter equation (BSE) [6, 7, 8],
only account for electronic effects on exciton states, typically computing
these at a fixed geometry of the studied system. However, even at $0$ K,
atomic nuclei vibrate according to normal modes with a zero-point energy of
$\frac{1}{2}\hbar\omega$. In a recent study on solid state organic
semiconductors [9], two of us showed that nuclear quantum motion can
significantly change exciton energies computed at the ground state geometry of
a system, and that incorporating these effects leads to improved agreement
with experiment.

Here we systematically investigate the effect of nuclear quantum motion on the
exciton energies of organic molecules. We provide an intuitive explanation of
why this will lead to a red shift of exciton energies in the vast majority of
molecules, and based on a harmonic approximation we derive an analytic
expression quantifying this effect and relating it to the frequency difference
between the vibrational normal modes of the ground and excited states.
Although this simple model provides useful insights, quantitative calculations
often require a more sophisticated treatment due to the anharmonic character
of excited state potential energy surfaces.

We therefore use a Monte Carlo sampling technique, combining TD-DFT and finite
difference methods for the molecular vibrations [10, 11], allowing us to treat
excited state energy surfaces without any harmonic assumption. We apply this
method to the so-called Thiel set of organic molecules, for which highly
accurate exciton energies have been computed using wavefunction-based methods
[12]. The Thiel set of molecules consists of four broad categories of
structures: unsaturated aliphatic hydrocarbons, aromatic hydrocarbons and
heterocycles, carbonyl compounds, and nucleobases. A full list of the studied
structures is given in section S1 of the Supporting Information (SI). It has
become common the the literature to compare new methods for the calculation of
exciton energies to the Thiel set [13, 7, 14, 15, 8], a path that we also
take. Remarkably, we find that the inclusion of nuclear quantum effects into
our TD-DFT calculations leads to significantly better agreement with
experimental results than both the (uncorrected) TD-DFT and high-level Thiel
calculations. This highlights the fact that corrections to exciton energies
due to nuclear quantum motion can be as important, and occasionally more
important than those that are calculated via highly accurate beyond-TD-DFT
electronic structure calculations; an effect that is even more pronounced for
small molecules with very strong exciton-vibration interactions. The method we
use is computationally inexpensive compared to the highly accurate electronic
structure methods used to go beyond TD-DFT, making it a practical way of
achieving greater accuracy in exciton energy calculations. While here we
employ TD-DFT to solve the electronic structure problem, our methodology for
correcting exciton energies can be used in conjunction with any electronic
structure method, including high-level techniques such as CAS-SCF.

Figure 1: Schematic representation of the energy of a $\pi$ and $\pi^{*}$
orbital, as well as the exciton state that results from a transition between
them, along a generalized coordinate $u$ (panel a). Due to the excited state
potential energy surface being less steep, nuclear quantum motion leads to a
red-shift of the energy difference between the ground and excited states
compared to its static value (panel b).

Let us start with a qualitative discussion of the effects of nuclear quantum
motion on exciton energies. In organic molecules, the highest occupied
molecular orbital (HOMO) is a bonding $\pi$ orbital that lowers the energy of
the molecule once occupied and leads to nuclei that are closer to each other.
Therefore, the energy of such a $\pi$ orbital along a generalized nuclear
coordinate $u$ of the molecular system will look similar to a Morse potential,
as visualized in Figure 1a (black). The lowest unoccupied molecular orbital
(LUMO) is an anti-bonding $\pi^{*}$ orbital that raises the energy of the
molecule once occupied and forces the nuclei apart, its energy along $u$
showing an exponential decay (Figure 1a, red). For most conjugated organic
molecules, the lowest energy singlet exciton ($\text{S}_{1}$) is formed by
exciting an electron from the HOMO to the LUMO [16], with energy:

$\displaystyle
E(\text{S}_{1})=E(\text{S}_{0})+\epsilon_{L}-\epsilon_{H}-J_{HL}+2K_{HL}.$ (1)

Here $E(\text{S}_{0})$ is the energy of the ground state, $\epsilon_{H}$ and
$\epsilon_{L}$ are the energies of the HOMO and the LUMO respectively,
$J_{HL}$ the HOMO-LUMO Coulomb integral and $K_{HL}$ the HOMO-LUMO exchange
integral. By using equation 1, we schematically describe the energy of an
exciton along $u$ (Figure 1a, blue), having assumed that the dependence of the
integrals $J,K$ on $u$ does not alter its qualitative characteristics.

What can be observed from Figure 1a is that the excited state potential energy
curve (blue) is less steep compared to the ground state one (black).
Therefore, if we approximate the two within the harmonic approximation, we
obtain two parabolas of different curvature, as shown in Figure 1b. If the
structure was ‘frozen’ at its ground state configuration $u_{\text{GS}}$, then
the energy required to access the exciton state would be $E_{\text{static}}$,
which is the energy gap usually computed by electronic structure calculations.
However, due to the quantum motion of the nuclei, the system can explore a
distribution of configurations within the region denoted by green arrows, even
at $0$ K. If we excited the system from any given displaced configuration
within the available region, then the excitation energy becomes smaller than
$E_{\text{static}}$, since the lower parabola corresponding to the ground
state is steeper (has higher frequency) than the excited state one. Therefore,
it is expected that inclusion of the effects of nuclear quantum fluctuations
will lead to a red-shift of exciton energies. Naturally, thermal motion of the
nuclei will have a similar result on exciton energies, however the light mass
of the constituent elements of organic molecules leads to vibrational modes
with high frequencies, which are only weakly activated at room temperature.
Additionally, thermal nuclear motion mostly affects the energies of the
spatially extended excited states that crystalline systems can host [9], which
is not the scenario we are examining here. We therefore focus on nuclear
motion at $T=0$ K, however our presented method could be used to study
exciton-vibration interactions at any temperature.

Now let us quantify the change in exciton energies due to nuclear quantum
motion. At temperature $T$, we include the contribution of molecular
vibrations to the exciton energy $E_{\text{exc}}(T)$ by means of the standard
quantum mechanical expectation value:

$E_{\text{exc}}(T)=\frac{1}{\mathcal{Z}}\sum_{\mathbf{s}}\bra{\chi_{\mathbf{s}}(\boldsymbol{u})}E_{\text{exc}}(\boldsymbol{u})\ket{\chi_{\mathbf{s}}(\boldsymbol{u})}e^{-E_{\mathbf{s}}/k_{\mathrm{B}}T},$
(2)

where $\ket{\chi_{\mathbf{s}}(\boldsymbol{u})}$ is a vibrational eigenstate on
the ground state potential energy surface with energy $E_{\mathbf{s}}$,
$\mathcal{Z}=\sum_{\mathbf{s}}e^{-E_{\mathbf{s}}/k_{\mathrm{B}}T}$ is the
partition function, and $\boldsymbol{u}$ is the nuclear displacement. To
simplify the problem, we resort to the harmonic approximation for the ground
state potential, and by substituting $\ket{\chi_{\mathbf{s}}(\boldsymbol{u})}$
with the wavefunction of a quantum harmonic oscillator we obtain [11]:

$E_{\text{exc}}(T)=\int
d\boldsymbol{u}|\Phi(\boldsymbol{u};T)|^{2}E_{\text{exc}}(\boldsymbol{u}),$
(3)

where:

$|\Phi(\boldsymbol{u};T)|^{2}=\prod_{\nu}(2\pi\sigma^{2}_{\nu}(T))^{-1/2}\exp{\left(-\frac{u_{\nu}^{2}}{2\sigma^{2}_{\nu}(T)}\right)},$
(4)

is the harmonic density at temperature $T$, which in turn is a product of
Gaussian functions of width:

$\sigma^{2}_{\nu}(T)=\frac{1}{2\omega_{\nu}}\cdot\coth{\left(\frac{\omega_{\nu}}{2k_{B}T}\right)}.$
(5)

In the above, atomic units and mass-weighted coordinates have been used, and
$\nu$ is the index labeling the ground state vibrational modes of the studied
molecule. In order to make further progress, we assume that the excited state
surface is also parabolic, and taking the limit of $T=0$ K we find:

$E(T=0\hskip
2.84544ptK)=E_{\text{static}}-\frac{1}{4}\sum_{i=1}^{3N-6}\frac{\omega_{Gi}^{2}-\omega_{Ei}^{2}}{\omega_{Gi}}$
(6)

where $\omega_{Gi}$ and $\omega_{Ei}$ are the frequencies of the vibrational
mode $i$ of the ground and excited state respectively. For a molecule with $N$
atoms, there are $3N-6$ non-trivial vibrational modes in total, and equation 6
suggests that as expected $E(T=0\hskip 2.84544ptK)<E_{\text{static}}$, given
that the potential energy surface of the ground state is steeper than that of
the excited state and therefore $\omega_{Gi}^{2}>\omega_{Ei}^{2}$. The full
derivation of equation 6 is provided in section S4 of the SI.

While equation 6 confirms the intuitive expectation that nuclear quantum
motion leads to a red shift of exciton energies, the assumption that the
excited state surface along vibrational normal modes is harmonic is often
violated. Examples of strongly anharmonic excited state potentials are known
from the literature [17, 18], and here we exemplify this with the simulated
potential energy surfaces of formaldehyde, shown in Figure 2. Therefore,
equation 6 is of mostly qualitative value and one needs to go beyond this
model in order to achieve quantitative accuracy. However, we emphasize that
treating the ground state potential energy surface along vibrational
eigenstates as harmonic remains a good approximation in many situations.

Figure 2: Ground (black) and excited state (red) potential energy surfaces of
formaldehyde (left) along a vibrational mode with frequency $1190$
$\text{cm}^{-1}$.

To do so, we evaluate the integral in equation 3 using Monte Carlo techniques,
generating $100$ configurations at $T=0$ K distributed according to the
corresponding harmonic density. Then, the expectation value of equation 3 is
computed as the simple average of the computed values of the exciton energy
$E_{\text{exc}}(\boldsymbol{u})$ at the displaced configurations. To compute
the exciton energies, we employ TD-DFT using the popular B3LYP hybrid
functional [19] and the cc-pVDZ basis set, as implemented within the NWChem
code [20]. This computational approach relies on no adjustable parameters,
apart from those in the DFT functional, which are fixed throughout the entire
series of calculations. We hence obtain zero-point renormalized (ZPR) exciton
energies compared to those of a ‘static’ TD-DFT calculation for the Thiel set
of molecules. We generally find that approximately $50$ points are sufficient
to sample the expectation value of the exciton energies. The numerical results
of these calculations are given in SI section S2, not only at the B3LYP level,
but also for three additional DFT functionals as discussed below.

We now proceed to compare the energies of the singlet excitons obtained within
our TD-DFT+ZPR approach that includes nuclear quantum motion to the
experimental values collected in the original publication on the Thiel set
[12]. We also compare the values obtained from ‘static’ TD-DFT and from highly
accurate complete active space second-order perturbation theory (CASPT2) [21]
calculations in the original Thiel publication [12] to experiment. The
computed exciton energies always refer to single excitations as described
within TD-DFT, and we exclude any double-excitations (bi-excitons) from our
analysis, which are however accessible using wavefunction-based methods [12].
In Figure 3 we plot the computed versus experimental values for the three
approaches, for all the studied molecules for which experimental data exists.
The closer a point lies to the $y=x$ line, the better the agreement between
theory and experiment. We fit a linear model to the three sets of results, as
a means of visualizing the overall agreement with experiment, from where it
becomes evident that the TD-DFT B3LYP+ZPR (red) results provide a significant
improvement to the static TD-DFT B3LYP values (green). The parameters for all
the linear fits are given in SI section S3. Remarkably, we find that the
values we compute with the inclusion of nuclear quantum fluctuations are in
better agreement with experiment even compared to the benchmark Thiel values
(black) that were obtained using accurate wavefunction-based methods. While
naturally the Thiel values provide an improvement to static TD-DFT, the
correction to exciton energies induced by the quantum fluctuations of
molecular vibrations seems to be more significant than corrections due to
beyond-TD-DFT electronic effects. We now proceed to further quantity these
comparisons.

Figure 3: Comparison of the different methods used to obtain the singlet
exciton energies of the Thiel set to experiment. A ‘static’ TD-DFT approach at
the B3LYP level (green) performs the worst, and accurate wavefunction-based
methods reported in the literature (black) improve overall agreement (the
values in black refer to the CASPT2 Thiel values [12] for single excitations,
corresponding to the states computed within TD-DFT here). However, accounting
for the zero-point renormalization (ZPR) of exciton energies due to molecular
quantum fluctuations leads to the best agreement with experiment (red), even
when ignoring beyond-TD-DFT electronic effects.

We employ four different statistical measures in order to compare the accuracy
of the different methods presented in Figure 3, which we summarize in Table 1.
In particular, for our $n=24$ studied molecules with computed exciton energies
$b_{i}$ and experimental exciton energies $a_{i}$, we compute the average
values of the bias: $\text{bias}=\frac{1}{n}\sum_{i}(b_{i}-a_{i})$, the root
mean-squared error: $\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i}(b_{i}-a_{i})^{2}}$,
as well as the relative values to these quantities:
$(\text{rel.bias})=\frac{1}{n}\sum\frac{b_{i}-a_{i}}{a_{i}}$,
$(\text{rel.RMSE})=\sqrt{\frac{1}{N}\sum_{i}\frac{(b_{i}-a_{i})^{2}}{a_{i}^{2}}}$.
The reason that we also consider the relative quantities is that the largest
errors for the static TD-DFT and Thiel approaches are found in the region of
large exciton energies, hence these molecules dominate the bias and root mean-
squared error. For example, the static TD-DFT exciton energy of ethene is
renormalized from $8.8$ eV to $7.7$ eV once we account for nuclear quantum
effects. Using the relative bias and root mean-squared error, which are
calculated by dividing by $a_{i}$ and $a_{i}^{2}$ respectively, we achieve a
fairer comparison without under-representing the effect of molecules with
lower exciton energies. We see from Table 1 that the average value of all the
employed statistical measures is minimized when using TD-DFT at the B3LYP
level combined with the zero-point renormalization (B3LYP+ZPR) of exciton
energies due to molecular vibrations. We therefore quantitatively confirm the
observation of Figure 3: corrections to static exciton energies due to nuclear
quantum motion are generally larger than corrections due to beyond-TD-DFT
electronic effects that are included in the Thiel values. We also quantify the
spread of the four statistical quantities around their average value by
computing their standard deviation, which is given in SI Table S32. We find
that the spread of the B3LYP+ZPR values is in every case comparable to that of
the Thiel values, and significantly smaller than that of the static TD-DFT
(B3LYP) values.

| bias (eV) | rel.bias | RMSE (eV) | rel.RMSE
---|---|---|---|---
B3LYP | $0.385$ | $0.073$ | $0.523$ | $0.098$
Thiel | $0.303$ | $0.060$ | $0.380$ | $0.072$
B3LYP+ZPR | $0.031$ | $0.005$ | $0.245$ | $0.051$
Table 1: Statistical measures of the accuracy of the different methods. The
average values of the (relative) bias and the (relative) root mean-squared
error (RMSE) as defined in the text are given.

Ideally, one would combine accurate electronic structure calculations with our
approach for accounting for nuclear quantum motion. However, for most of the
studied molecules we found that a minimum of approximately $50$ points is
required to converge the Monte Carlo sampling of the integral in equation 3.
While the cost of $50$ TD-DFT calculations of small molecules is reasonably
low, it quickly becomes very high as one moves to more accurate wavefunction-
based methods. It is therefore reasonable to wonder whether the vibrational
corrections to exciton energies could be computed at one (cheaper) level of
theory and applied to the static exciton energy obtained at another (more
accurate) level of theory. To check whether the vibrational correction to
exciton energies changes when using different electronic structure methods, we
repeated our Monte Carlo sampling within three additional levels of TD-DFT,
while using the vibrational modes obtained at the B3LYP level. In particular,
for the excited state calculations we used the local density approximation
(LDA), the generalized gradient approximation (GGA) as well as a functional
with pure Hartree-Fock exchange. The results are presented in SI section S3.
We find that the correction to static exciton energies can change
significantly depending on the employed functional, suggesting it is not
possible to accurately obtain the static and displaced exciton values at
different levels of theory. Additionally, we find for most molecules at the
LDA level that computing both the vibrational modes and excited states using
this functional leads to small changes to the corrections of exciton energies,
compared to using B3LYP vibrations along with LDA excited state calculations.
Therefore, the change in the exciton energy correction between the different
functionals seems to be due to changes in the excited state wavefunction and
not in the molecular vibrational properties.

Figure 4: Size dependence of the correction due to nuclear quantum motion to
the singlet exciton energies $\Delta\text{E}(\text{S}_{1})$ of the acene
(panel a) and aliphatic hydrocarbon (panel b) families of molecules.

We found in Ref. [9] for solid state systems that the correction to exciton
energies due to nuclear quantum motion becomes more important for smaller
molecules. For the single-molecule systems studied here we find that the same
trend holds, as shown in Figure 4 for the acene and aliphatic hydrocarbons
families of molecules. This is because smaller molecules have a greater
electronic density in the vicinity of localized carbon motions that are
activated due to quantum fluctuations, hence leading to stronger exciton-
vibration interactions. Another way of reaching the same conclusion is using
Hückel theory. Let us consider the example of linear polyenes, such as ethene,
butadiene, hexatriene and octatetraene that are studied here and are included
among other aliphatic hydrocarbons in Figure 4b, showing a decrease in the
exciton energy correction with increasing system size. For a linear polyene
consisting of $N$ carbon atoms, Hückel theory predicts the energy of the
$i^{\text{th}}$ orbital is: $E_{i}=\alpha+2\beta\cos{(2\pi n/(N+1))}$, with
$\alpha,\beta$ the Hückel parameters, and $i=1,2,...,N$. From this, we find
for the HOMO-LUMO gap:

$E_{\text{gap}}\propto\sin(\frac{\pi}{2(N+1)}).$ (7)

In the limit of $N\rightarrow\infty$, this gives $E_{\text{gap}}\propto 1/N$.
A similar result can be obtained for cyclic polyenes. Assuming that the same
holds for any general conjugated molecule, this suggests that in larger
systems the HOMO is less bonding and the LUMO less anti-bonding than in
comparably smaller systems. Therefore, from the intuitive picture of Figure 1,
the curvature of the ground and excited state surfaces of large molecules is
similar, and the difference in their normal mode frequencies is small. Hence,
according to the qualitative expression 6, the correction to the exciton
energy is also smaller than that for a smaller molecule. This is encouraging,
in the sense that corrections due to nuclear quantum motion are mostly
important for small systems, for which they are cheaper to compute.

In conclusion, we have presented an intuitive picture and a computational
approach for estimating the effects of nuclear quantum motion on molecular
exciton energies, which we show can indeed be substantial, particularly for
smaller molecules. The corrections that molecular vibrations induce to exciton
energies computed at a static level are found to be greater than corrections
due to beyond-TD-DFT electronic effects as these were identified in previous
studies. One could further improve upon these results by combining our
computational methodology for nuclear quantum effects with a high-level theory
of electronic structure, such as CAS-SCF. Overall, we underline the need for
the community to explore nuclear quantum motion as a promising avenue for
achieving quantitative accuracy when computing molecular excited state
energies. While here we do not provide a mode-resolved picture, exploring
which vibrational modes provide the largest correction to exciton energies
could be of great interest and will be the subject of a future work.

## Data availability

The data underlying this publication can be found under [URL added in proof].

## Acknowledgments

The authors acknowledge useful discussions with Stuart Althorpe (Cambridge)
and the support of the Winton Programme for the Physics of Sustainability.
T.J.H.H. acknowledges a Royal Society University Research Fellowship Ref:
URF\R1\201502. B.M. acknowledges support from the Gianna Angelopoulos
Programme for Science, Technology, and Innovation. A.M.A. acknowledges the
support of the Engineering and Physical Sciences Research Council for funding
under grant EP/L015552/1.

## References

  * [1] J. Frenkel.  On the transformation of light into heat in solids. I.  Physical Review, 37:17–44, 1931.
  * [2] Akshay Rao and Richard H Friend.  Harnessing singlet exciton fission to break the Shockley–Queisser limit.  Nature Reviews Materials, 2:17063, October 2017.
  * [3] Sebastian Reineke.  Organic light-emitting diodes: Phosphorescence meets its match.  Nature Photonics, 8(4):269–270, 2014.
  * [4] M. Petersilka, U. J. Gossmann, and E. K.U. Gross.  Excitation energies from time-dependent density-functional theory.  Physical Review Letters, 76(8):1212–1215, 1996.
  * [5] Per Åke Malmqvist and Björn O. Roos.  The CASSCF state interaction method.  Chemical Physics Letters, 155(2):189–194, 1989.
  * [6] Michael Rohlfing and Steven G. Louie.  Electron-hole excitations and optical spectra from first principles.  Physical Review B - Condensed Matter and Materials Physics, 62(8):4927–4944, 2000.
  * [7] Denis Jacquemin, Ivan Duchemin, and Xavier Blase.  Benchmarking the Bethe-Salpeter Formalism on a Standard Organic Molecular Set.  Journal of Chemical Theory and Computation, 11(7):3290–3304, 2015\.
  * [8] Denis Jacquemin, Ivan Duchemin, Aymeric Blondel, and Xavier Blase.  Benchmark of Bethe-Salpeter for Triplet Excited-States.  Journal of Chemical Theory and Computation, 13(2):767–783, 2017\.
  * [9] Antonios M. Alvertis, Raj Pandya, Loreta A. Muscarella, Nipun Sawhney, Malgorzata Nguyen, Bruno Ehrler, Akshay Rao, Richard H. Friend, Alex W. Chin, and Bartomeu Monserrat.  Impact of exciton delocalization on exciton-vibration interactions in organic semiconductors.  Physical Review B - Condensed Matter and Materials Physics, 102:081122(R), 2020.
  * [10] K. Kunc and Richard M. Martin.  Ab initio force constants of GaAs: A new approach to calculation of phonons and dielectric properties.  Physical Review Letters, 48(6):406–409, 1982.
  * [11] Bartomeu Monserrat.  Electron – phonon coupling from finite differences.  Journal of Physics Condensed Matter, 30, 2018.
  * [12] Marko Schreiber, Mario R. Silva-Junior, Stephan P.A. Sauer, and Walter Thiel.  Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3.  Journal of Chemical Physics, 128(13), 2008.
  * [13] Subhayan Roychoudhury, Stefano Sanvito, and David D. O’Regan.  Neutral excitation density-functional theory: an efficient and variational first-principles method for simulating neutral excitations in molecules.  Scientific Reports, 10(1):1–12, 2020.
  * [14] Fabien Bruneval, Samia M. Hamed, and Jeffrey B. Neaton.  A systematic benchmark of the ab initio Bethe-Salpeter equation approach for low-lying optical excitations of small organic molecules.  Journal of Chemical Physics, 142(24), 2015.
  * [15] Tonatiuh Rangel, Samia M. Hamed, Fabien Bruneval, and Jeffrey B. Neaton.  An assessment of low-lying excitation energies and triplet instabilities of organic molecules with an ab initio Bethe-Salpeter equation approach and the Tamm-Dancoff approximation.  Journal of Chemical Physics, 146(19), 2017.
  * [16] P. W. (Peter William) Atkins.  Molecular quantum mechanics.  Oxford University Press, New York, 4th ed. / peter atkins, ronald friedman. edition, 2005.
  * [17] Franco Egidi, David B. Williams-Young, Alberto Baiardi, Julien Bloino, Giovanni Scalmani, Michael J. Frisch, Xiaosong Li, and Vincenzo Barone.  Effective Inclusion of Mechanical and Electrical Anharmonicity in Excited Electronic States: VPT2-TDDFT Route.  Journal of Chemical Theory and Computation, 13(6):2789–2803, 2017\.
  * [18] Antonios M Alvertis, Steven Lukman, Timothy J H Hele, Eric G Fuemmeler, Jiaqi Feng, Jishan Wu, Neil C Greenham, Alex W Chin, and Andrew J Musser.  Switching between Coherent and Incoherent Singlet Fission via Solvent-Induced Symmetry Breaking.  Journal of the American Chemical Society, 141:17558–17570, 2019\.
  * [19] Axel D. Becke.  A new mixing of Hartree-Fock and local density-functional theories.  The Journal of Chemical Physics, 98(2):1372–1377, 1993.
  * [20] E. Aprà _et al._ NWChem: Past, present, and future.  The Journal of chemical physics, 152(18):184102, 2020.
  * [21] Kerstin Andersson, Per Åke Malmqvist, Björn O. Roos, Andrzej J. Sadlej, and Krzysztof Wolinski.  Second-order perturbation theory with a CASSCF reference function.  Journal of Physical Chemistry, 94(14):5483–5488, 1990.