# An Online Learning Approach to
Model Predictive Control

Nolan Wagener,1$\\#$ Ching-An Cheng,1$\\#$ Jacob Sacks,2 and Byron Boots1
Georgia Institute of Technology
{nolan.wagener, cacheng<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Model predictive control (MPC) is a powerful technique for solving dynamic
control tasks. In this paper, we show that there exists a close connection
between MPC and online learning, an abstract theoretical framework for
analyzing online decision making in the optimization literature. This new
perspective provides a foundation for leveraging powerful online learning
algorithms to design MPC algorithms. Specifically, we propose a new algorithm
based on dynamic mirror descent (DMD), an online learning algorithm that is
designed for non-stationary setups. Our algorithm, Dynamic Mirror Descent
Model Predictive Control (DMD-MPC), represents a general family of MPC
algorithms that includes many existing techniques as special instances. DMD-
MPC also provides a fresh perspective on previous heuristics used in MPC and
suggests a principled way to design new MPC algorithms. In the experimental
section of this paper, we demonstrate the flexibility of DMD-MPC, presenting a
set of new MPC algorithms on a simple simulated cartpole and a simulated and
real-world aggressive driving task. Videos of the real-world experiments can
be found at https://youtu.be/vZST3v0_S9w and https://youtu.be/MhuqiHo2t98.

##  I Introduction

††1Institute for Robotics and Intelligent Machines
2School of Electrical and Computer Engineering
$\\#$Equal contribution

Model predictive control (MPC) [20] is an effective tool for control tasks
involving dynamic environments, such as helicopter aerobatics [1] and
aggressive driving [30]. One reason for its success is the pragmatic principle
it adopts in choosing controls: rather than wasting computational power to
optimize a complicated controller for the full-scale problem (which may be
difficult to accurately model), MPC instead optimizes a simple controller
(e.g., an open-loop control sequence) over a shorter planning horizon that is
just sufficient to make a sensible decision at the current moment. By
alternating between optimizing the simple controller and applying its
corresponding control on the real system, MPC results in a closed-loop policy
that can handle modeling errors and dynamic changes in the environment.

Various MPC algorithms have been proposed, using tools ranging from
constrained optimization techniques [7, 20, 27] to sampling-based techniques
[30]. In this paper, we show that, while these algorithms were originally
designed differently if we view them through the lens of _online learning_
[16], many of them actually follow the same general update rule. Online
learning is an abstract theoretical framework for analyzing online decision
making. Formally, it concerns iterative interactions between a learner and an
environment over $T$ rounds. At round $t$, the learner makes a decision
$\tilde{\bm{\theta}}_{t}$ from some decision set $\bm{\Theta}$. The
environment then chooses a loss function $\ell_{t}$ based on the learner’s
decision, and the learner suffers a cost $\ell_{t}(\tilde{\bm{\theta}}_{t})$.
In addition to seeing the decision’s cost, the learner may be given additional
information about the loss function (e.g., its gradient evaluated at
$\tilde{\bm{\theta}}_{t}$) to aid in choosing the next decision
$\tilde{\bm{\theta}}_{t+1}$. The learner’s goal is to minimize the accumulated
costs $\sum_{t=1}^{T}\ell_{t}(\tilde{\bm{\theta}}_{t})$, e.g., by minimizing
regret [16].

We find that the MPC process bears a strong similarity with online learning.
At time $t$ (i.e., round $t$), an MPC algorithm optimizes a controller (i.e.,
the decision) over some cost function (i.e., the per-round loss). To do so, it
observes the cost of the initial controller (i.e.,
$\ell_{t}(\tilde{\bm{\theta}}_{t})$), improves the controller, and executes a
control based on the improved controller in the environment to get to the next
state (which in turn defines the next per-round loss) with a new controller
$\tilde{\bm{\theta}}_{t+1}$.

In view of this connection, we propose a generic framework, _DMD-MPC_ (Dynamic
Mirror Descent Model Predictive Control), for synthesizing MPC algorithms.
DMD-MPC is based on a first-order online learning algorithm called dynamic
mirror descent (DMD) [14], a generalization of mirror descent [4] for dynamic
comparators. We show that several existing MPC algorithms [31, 32] are special
cases of DMD-MPC, given specific choices of step sizes, loss functions, and
regularization. Furthermore, we demonstrate how new MPC algorithms can be
derived systematically from DMD-MPC with only mild assumptions on the
regularity of the cost function. This allows us to even work with
discontinuous cost functions (like indicators) and discrete controls. Thus,
DMD-MPC offers a spectrum from which practitioners can easily customize new
algorithms for their applications.

In the experiments, we apply DMD-MPC to design a range of MPC algorithms and
study their empirical performance. Our results indicate the extra design
flexibility offered by DMD-MPC does make a difference in practice; by properly
selecting hyperparameters which are obscured in the previous approaches, we
are able to improve the performance of existing algorithms. Finally, we apply
DMD-MPC on a real-world AutoRally car platform [13] for autonomous driving
tasks and show it can achieve competent performance.

Notation: As our discussions will involve planning horizons, for clarity, we
use lightface to denote variables that are meant for a single time step, and
boldface to denote the variables congregated across the MPC planning horizon.
For example, we use $\hat{u}_{t}$ to denote the planned control at time $t$
and $\hat{\bm{u}}_{t}\triangleq(\hat{u}_{t},\ldots,\hat{u}_{t+H-1})$ to denote
an $H$-step planned control sequence starting from time $t$. We use a
subscript to extract elements from a congregated variable; e.g., we use
$\hat{u}_{t,h}$ to the denote the $h$${}^{\text{th}}$ element in
$\hat{\bm{u}}_{t}$ (the subscript index starts from zero). All the variables
in this paper are finite-dimensional.

##  II An Online Learning Perspective on MPC

###  II-A The MPC Problem Setup

Let $n,~{}m\in\mathbb{N}_{+}$ be finite. We consider the problem of
controlling a discrete-time stochastic dynamical system

$\displaystyle x_{t+1}\sim f(x_{t},u_{t})$ (1)

for some stochastic transition map
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$. At time $t$, the
system is in state $x_{t}\in\mathbb{R}^{n}$. Upon the execution of control
$u_{t}\in\mathbb{R}^{m}$, the system randomly transitions to the next state
$x_{t+1}$, and an instantaneous cost $c(x_{t},u_{t})$ is incurred. Our goal is
to design a state-feedback control law (i.e., a rule of choosing $u_{t}$ based
on $x_{t}$) such that the system exhibits good performance (e.g., accumulating
low costs over $T$ time steps).

In this paper, we adopt the MPC approach to choosing $u_{t}$: at state
$x_{t}$, we imagine controlling a stochastic dynamics model $\hat{f}$ (which
approximates our system $f$) for $H$ time steps into the future. Our planned
controls come from a control distribution $\bm{\pi}_{\bm{\theta}}$ that is
parameterized by some vector $\bm{\theta}\in\bm{\Theta}$, where $\bm{\Theta}$
is the feasible parameter set. In each simulation (i.e., rollout), we
sample111This can be sampled in either an open-loop or closed-loop fashion. a
control sequence $\hat{\bm{u}}_{t}$ from the control distribution
$\bm{\pi}_{\bm{\theta}}$ and recursively apply it to $\hat{f}$ to generate a
predicted state trajectory
$\hat{\bm{x}}_{t}\triangleq(\hat{x}_{t},\hat{x}_{t+1},\ldots,\hat{x}_{t+H})$:
let $\hat{x}_{t}=x_{t}$; for $\tau=t,\dots,t+H-1$, we set
$\hat{x}_{\tau+1}\sim\hat{f}(\hat{x}_{\tau},\hat{u}_{\tau})$. More compactly,
we can write the simulation process as

$\displaystyle\hat{\bm{x}}_{t}\sim\hat{\bm{f}}(x_{t},\hat{\bm{u}}_{t})$ (2)

in terms of some $\hat{\bm{f}}$ that is defined naturally according to the
above recursion. Through these simulations, we desire to select a parameter
$\bm{\theta}_{t}\in\bm{\Theta}$ that minimizes an MPC objective
$\hat{J}(\bm{\pi}_{\bm{\theta}};x_{t})$, which aims to predict the performance
of the system if we were to apply the control distribution
$\bm{\pi}_{\bm{\theta}}$ starting from $x_{t}$.222$\hat{J}$ can be seen as a
surrogate for the long-term performance of our controller. Typically, we set
the planning horizon $H$ to be much smaller than $T$ to reduce the
optimization difficulty and to mitigate modeling errors. In other words, we
wish to find the $\bm{\theta}_{t}$ that solves

$\min_{\bm{\theta}\in\bm{\Theta}}\hat{J}(\bm{\pi}_{\bm{\theta}};x_{t}).$ (3)

Once $\bm{\theta}_{t}$ is decided, we then sample333This setup can also
optimize deterministic policies, e.g., by defining $\bm{\pi}_{\bm{\theta}}$ to
be a Gaussian policy with the mean being the deterministic policy.
$\hat{\bm{u}}_{t}$ from $\bm{\pi}_{\bm{\theta}_{t}}$, extract the first
control $\hat{u}_{t}$, and apply it on the real dynamical system $f$ in (1)
(i.e., set $u_{t}=\hat{u}_{t}$) to go to the next state $x_{t+1}$. Because
$\bm{\theta}_{t}$ is determined based on $x_{t}$, MPC is effectively state-
feedback.

The motivation behind MPC is to use the MPC objective $\hat{J}$ to reason
about the controls required to achieve desirable long-term behaviors. Consider
the statistic

$\displaystyle
C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\triangleq\sum_{h=0}^{H-1}c(\hat{x}_{t+h},\hat{u}_{t+h})+c_{\mathrm{end}}(\hat{x}_{t+H}),$
(4)

where $c_{\mathrm{end}}$ is a terminal cost function. A popular MPC objective
is
$\hat{J}(\bm{\pi}_{\bm{\theta}};x_{t})=\mathbb{E}[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\mid
x_{t},\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}]$, which estimates the expected
$H$-step future costs. Later in Section III-A, we will discuss several MPC
objectives and their properties.

Although the idea of MPC sounds intuitively promising, the optimization can
only be approximated in practice (e.g., using an iterative algorithm like
gradient descent), because (3) is often a stochastic program (like the example
above) and the control command $u_{t}$ needs to be computed at a high
frequency. In consideration of this imperfection, it is common to
heuristically _bootstrap_ the previous approximate solution as the
initialization to the current problem. Specifically, let $\bm{\theta}_{t-1}$
be the approximate solution to the previous problem and
$\tilde{\bm{\theta}}_{t}$ denote the initial condition of $\bm{\theta}$ in
solving (3). The bootstrapping step can then written as

$\displaystyle\tilde{\bm{\theta}}_{t}=\Phi(\bm{\theta}_{t-1})$ (5)

by effectively defining a _shift operator_ $\Phi$ (see Appendix A for
details). Because the subproblems in (3) of two consecutive time steps share
all control variables except for the first and the last ones, shifting the
previous solution provides a warm start to (3) to amortize the computational
complexity.

###  II-B The Online Learning Perspective

$x_{t}$$x_{t-1}$$x_{t+1}$$u_{t-2}$$u_{t-1}$$u_{t}$$\tilde{\bm{\theta}}_{t-1}$$\bm{\theta}_{t-1}$$\tilde{\bm{\theta}}_{t}$$\bm{\theta}_{t}$$\tilde{\bm{\theta}}_{t+1}$$\bm{\theta}_{t+1}$$\ell_{t-1}$$\ell_{t}$$\ell_{t+1}$round
$t-1$round $t$round $t+1$$\Phi$$\Phi$ Figure 1: Diagram of the online learning
perspective, where blue and red denote the learner and the environment,
respectively.

As discussed, the iterative update process of MPC resembles the setup of
online learning [16]. Here we provide the details to convert an MPC setup into
an online learning problem. Recall from the introduction that online learning
mainly consists of three components: the decision set, the learner’s strategy
for updating decisions, and the environment’s strategy for updating per-round
losses. We show the counterparts in MPC that correspond to each component
below. Note that in this section we will overload the notation
$\hat{J}(\bm{\theta};x_{t})$ to mean $\hat{J}(\bm{\pi}_{\bm{\theta}};x_{t})$.

We use the concept of per-round loss in online learning as a mechanism to
measure the decision uncertainty in MPC, and propose the following
identification (shown in Fig. 1) for the MPC setup described in the previous
section: we set the rounds in online learning to synchronize with the time
steps of our control system, set the decision set $\bm{\Theta}$ as the space
of feasible parameters of the control distribution $\bm{\pi}_{\bm{\theta}}$,
set the learner as the MPC algorithm which in round $t$ outputs the decision
$\tilde{\bm{\theta}}_{t}\in\bm{\Theta}$ and side information $u_{t-1}$, and
set the per-round loss as

$\displaystyle\ell_{t}(\cdot)=\hat{J}(\cdot\,;x_{t}).$ (6)

In other words, in round $t$ of this online learning setup, the learner plays
a decision $\tilde{\bm{\theta}}_{t}$ along with a side information $u_{t-1}$
(based on the optimized solution $\bm{\theta}_{t-1}$ and the shift operator in
(5)), the environment selects the per-round loss
$\ell_{t}(\cdot)=\hat{J}(\cdot;x_{t})$ (by applying $u_{t-1}$ to the real
dynamical system in (1) to transit the state to $x_{t}$), and finally the
learner receives $\ell_{t}$ and incurs cost
$\ell_{t}(\tilde{\bm{\theta}}_{t})$ (which measures the sub-optimality of the
future plan made by the MPC algorithm).

This online learning setup differs slightly from the standard setup in its
separation of the decision $\tilde{\bm{\theta}}_{t}$ and the side information
$u_{t-1}$; while our setup can be converted into a standard one that treats
$\bm{\theta}_{t-1}$ as the sole decision played in round $t$, we adopt this
explicit separation in order to emphasize that the variable part of the
incurred cost $\ell_{t}(\tilde{\bm{\theta}}_{t})$ pertains to only
$\tilde{\bm{\theta}}_{t}$. That is, the learner cannot go back and revert the
previous control $u_{t-1}$ already applied on the system, but only uses
$\ell_{t}$ to update the current and future controls
$\hat{u}_{t},\dots,\hat{u}_{t+H-1}$.

The performance of the learner in online learning (which by our identification
is the MPC algorithm) is measured in terms of the accumulated costs
$\sum_{t=1}^{T}\ell_{t}(\tilde{\bm{\theta}}_{t})$. For problems in non-
stationary setups, a normalized way to describe the accumulated costs in the
online learning literature is through the concept of _dynamic regret_ [14,
34], which is defined as

$\displaystyle\mathrm{D}\text{-}\mathrm{Regret}=\sum_{t=1}^{T}\ell_{t}(\tilde{\bm{\theta}}_{t})-\sum_{t=1}^{T}\ell_{t}(\bm{\theta}_{t}^{\star}),$
(7)

where
$\bm{\theta}_{t}^{\star}\in\arg\min_{\bm{\theta}\in\bm{\Theta}}\,\ell_{t}(\bm{\theta})$.
Dynamic regret quantifies how suboptimal the played decisions
$\tilde{\bm{\theta}}_{1},\dots,\tilde{\bm{\theta}}_{T}$ are on the
corresponding loss functions. In our proposed problem setup, the optimality
concept associated with dynamic regret conveys a _consistency_ criterion
desirable for MPC: we would like to make a decision $\bm{\theta}_{t-1}$ at
state $x_{t-1}$ such that, after applying control $u_{t-1}$ and entering the
new state $x_{t}$, its shifted plan $\tilde{\bm{\theta}}_{t}$ remains close to
optimal with respect to the new loss function $\ell_{t}$. If the dynamics
model $\hat{f}$ is accurate and the MPC algorithm is ideally solving (3), we
can expect that bootstrapping the previous solution $\bm{\theta}_{t-1}$
through (5) into $\tilde{\bm{\theta}}_{t}$ would result in a small
instantaneous gap
$\ell_{t}(\tilde{\bm{\theta}}_{t})-\ell_{t}(\bm{\theta}_{t}^{\star})$ which is
solely due to unpredictable future information (such as the stochasticity in
the dynamical system). In other words, an online learning algorithm with small
dynamic regret, if applied to our online learning setup, would produce a
consistently optimal MPC algorithm with regard to the solution concept
discussed above. However, we note that having small dynamic regret here does
not directly imply good absolute performance on the control system, because
the overall performance of the MPC algorithm is largely dependent on the form
of the MPC objective $\hat{J}$ (e.g., through choice of $H$ and accuracy of
$\hat{f}$). Small dynamic regret more precisely means whether the plan
produced by an MPC algorithm is consistent with the given MPC objective.

##  III A Family of MPC Algorithms Based on Dynamic Mirror Descent

The online learning perspective on MPC suggests that good MPC algorithms can
be designed from online learning algorithms that achieve small dynamic regret.
This is indeed the case. We will show that a range of existing MPC algorithms
are in essence applications of a classical online learning algorithm called
dynamic mirror descent (DMD) [14]. DMD is a generalization of mirror descent
[4] to problems involving dynamic comparators (in this case, the
$\\{\bm{\theta}_{t}^{\star}\\}$ in dynamic regret in (7)). In round $t$, DMD
applies the following update rule:

$\displaystyle\begin{split}\bm{\theta}_{t}&=\arg\min_{\bm{\theta}\in\bm{\Theta}}\,\langle\gamma_{t}\bm{g}_{t},\bm{\theta}\rangle+D_{\psi}(\bm{\theta}\|\tilde{\bm{\theta}}_{t})\\\
\tilde{\bm{\theta}}_{t+1}&=\Phi(\bm{\theta}_{t}),\end{split}$ (8)

where $\bm{g}_{t}=\nabla\ell_{t}(\tilde{\bm{\theta}}_{t})$ (which can be
replaced by unbiased sampling if $\nabla\ell_{t}(\tilde{\bm{\theta}}_{t})$ is
an expectation), $\Phi$ is called the _shift model_ ,444In [14], $\Phi$ is
called a _dynamical model_ , but it is not the same as the dynamics of our
control system. We therefore rename it to avoid confusion. $\gamma_{t}>0$ is
the step size, and for some $\bm{\theta},\bm{\theta}^{\prime}\in\bm{\Theta}$,
$D_{\psi}(\bm{\theta}\|\bm{\theta}^{\prime})\triangleq\psi(\bm{\theta})-\psi(\bm{\theta}^{\prime})-\langle\nabla\psi(\bm{\theta}^{\prime}),\bm{\theta}-\bm{\theta}^{\prime}\rangle$
is the Bregman divergence generated by a strictly convex function $\psi$ on
$\bm{\Theta}$.

The first step of DMD in (8) is reminiscent of the proximal update in the
usual mirror descent algorithm. It can be thought of as an optimization step
where the Bregman divergence acts as a regularization to keep $\bm{\theta}$
close to $\tilde{\bm{\theta}}_{t}$. Although
$D_{\psi}(\bm{\theta}\|\bm{\theta}^{\prime})$ is not necessarily a metric
(since it may not be symmetric), it is still useful to view it as a distance
between $\bm{\theta}$ and $\bm{\theta}^{\prime}$. Indeed, familiar examples of
the Bregman divergence include the squared Euclidean distance and KL
divergence555For probability distributions $p$ and $q$ over a random variable
$x$, the KL divergence is defined as
$\mathrm{KL}(p\,\|\,q)\triangleq\mathbb{E}_{x\sim
p}\mathopen{}\mathclose{{}\left[\log\frac{p(x)}{q(x)}}\right]$. [3].

The second step of DMD in (8) uses the shift model $\Phi$ to anticipate the
optimal decision for the next round. In the context of MPC, a natural choice
for the shift model is the shift operator in (5) defined previously in Section
II-A (hence the same notation), because the per-round losses in two
consecutive rounds here concern problems with shifted time indices. Hall and
Willett [14] show that the dynamic regret of DMD scales with how much the
optimal decision sequence $\\{\bm{\theta}_{t}^{\star}\\}$ deviates from $\Phi$
(i.e., $\sum_{t}\|\bm{\theta}_{t+1}^{\star}-\Phi(\bm{\theta}_{t}^{\star})\|)$,
which is proportional to the unpredictable elements of the problem.

for  _$t=1,2,\ldots,T$_ do

$\ell_{t}(\cdot)=\hat{J}(\cdot\,;x_{t})$

$\bm{\theta}_{t}=\displaystyle\arg\min_{\bm{\theta}\in\bm{\Theta}}\;\langle\gamma_{t}\nabla\ell_{t}(\tilde{\bm{\theta}}_{t}),\bm{\theta}\rangle+D_{\psi}(\bm{\theta}\|\tilde{\bm{\theta}}_{t})$

Sample $\hat{\bm{u}}_{t}\sim\bm{\pi}_{\bm{\theta}_{t}}$ and set
$u_{t}=\hat{u}_{t}$

Sample $x_{t+1}\sim f(x_{t},u_{t})$

$\tilde{\bm{\theta}}_{t+1}=\Phi(\bm{\theta}_{t})$

end for

Algorithm 1 Dynamic Mirror Descent MPC (DMD-MPC)

Applying DMD in (8) to the online learning problem described in Section II-B
leads to an MPC algorithm shown in Algorithm 1, which we call _DMD-MPC_. More
precisely, DMD-MPC represents a family of MPC algorithms in which a specific
instance is defined by a choice of:

  1. 1.

the MPC objective $\hat{J}$ in (6),

  2. 2.

the form of the control distribution $\bm{\pi}_{\bm{\theta}}$, and

  3. 3.

the Bregman divergence $D_{\psi}$ in (8).

Thus, we can use DMD-MPC as a generic strategy for synthesizing MPC
algorithms. In the following, we use this recipe to recreate several existing
MPC algorithms and demonstrate new MPC algorithms that naturally arise from
this framework.

###  III-A Loss Functions

We discuss several definitions of the per-round loss $\ell_{t}$, which all
result from the formulation in (6) but with different $\hat{J}$. These loss
functions are based on the statistic $C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$
defined in (4) which measures the $H$-step accumulated cost of a given
trajectory. For transparency of exposition, we will suppose henceforth that
the control distribution $\bm{\pi}_{\bm{\theta}}$ is open-loop666Note again
that even while using open-loop control distributions, the overall control law
of MPC is state-feedback.; similar derivations follow naturally for closed-
loop control distributions. For convenience of practitioners, we also provide
expressions of their gradients in terms of the likelihood-ratio
derivative777We assume the control distribution is sufficiently regular with
respect to its parameter so that the likelihood-ratio derivative rule holds.
[12]. For some function $L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$, all these
gradients shall have the form

$\nabla\ell_{t}(\bm{\theta})=\mathbb{E}_{\hat{\bm{u}}_{t}\sim\bm{\pi}_{\bm{\theta}}}\mathbb{E}_{\hat{\bm{x}}_{t}\sim\hat{\bm{f}}(x_{t},\hat{\bm{u}}_{t})}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})}\right].$
(9)

In short, we will denote
$\mathbb{E}_{\hat{\bm{u}}_{t}\sim\bm{\pi}_{\bm{\theta}}}\mathbb{E}_{\hat{\bm{x}}_{t}\sim\hat{\bm{f}}(x_{t},\hat{\bm{u}}_{t})}$
as $\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}$. These gradients in
practice are approximated by finite samples.

####  III-A1 Expected Cost

The most commonly used MPC objective is the $H$-step expected accumulated cost
function under model dynamics, because it directly estimates the expected
long-term behavior when the dynamics model $\hat{f}$ is accurate and $H$ is
large enough. Its per-round loss function is888In experiments, we subtract the
empirical average of the sampled costs from $C$ in (11) to reduce the
variance, at the cost of a small amount of bias.

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\right]$
(10) $\displaystyle\nabla\ell_{t}(\bm{\theta})$
$\displaystyle=\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})}\right].$
(11)

####  III-A2 Expected Utility

Instead of optimizing for average cost, we may care to optimize for some
preference related to the trajectory cost $C$, such as having the cost be
below some threshold. This idea can be formulated as a _utility_ that returns
a normalized score related to the preference for a given trajectory cost
$C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$. Specifically, suppose that $C$ is
lower bounded by zero999If this is not the case, let
$c_{\min}\triangleq\inf_{\hat{\bm{x}}_{t},\hat{\bm{u}}_{t}}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$,
which we assume is finite. We can then replace $C$ with
$\tilde{C}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\triangleq
C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})-c_{\min}$. and at some round $t$ define
the utility $U_{t}:\mathbb{R}_{+}\to[0,1]$ (i.e.,
$U_{t}:C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\mapsto
U_{t}(C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t}))$) to be a function with the
following properties: $U_{t}(0)=1$, $U_{t}$ is monotonically decreasing, and
$\lim_{z\to+\infty}U_{t}(z)=0$. These are sensible properties since we attain
maximum utility when we have zero cost, the utility never increases with the
cost, and the utility approaches zero as the cost increases without bound. We
then define the per-round loss as

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[U_{t}(C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t}))}\right]$
(12) $\displaystyle\nabla\ell_{t}(\bm{\theta})$
$\displaystyle=-\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[U_{t}(C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t}))\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})}\right]}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[U_{t}(C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t}))}\right]}.$
(13)

The gradient in (13) is particularly appealing when estimated with samples.
Suppose we sample $N$ control sequences
$\hat{\bm{u}}_{t}^{1},\dots,\hat{\bm{u}}_{t}^{N}$ from
$\bm{\pi}_{\bm{\theta}}$ and (for the sake of compactness) sample one state
trajectory from $\hat{\bm{f}}$ for each corresponding control sequence,
resulting in $\hat{\bm{x}}_{t}^{1},\dots,\hat{\bm{x}}_{t}^{N}$. Then the
estimate of (13) is a convex combination of gradients:

$\displaystyle\nabla\ell_{t}(\bm{\theta})\approx-\sum_{i=1}^{N}w_{i}\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t}^{i}),$

where $w_{i}=\frac{U_{t}(C_{i})}{\sum_{j=1}^{N}U_{t}(C_{j})}$ and
$C_{i}=C(\hat{\bm{x}}_{t}^{i},\hat{\bm{u}}_{t}^{i})$, for $i=1,\dots,N$. We
see that each weight $w_{i}$ is computed by considering the relative utility
of its corresponding trajectory. A cost $C_{i}$ with high relative utility
will push its corresponding weight $w_{i}$ closer to one, whereas a low
relative utility will cause $w_{i}$ to be close to zero, effectively rejecting
the corresponding sample.

We give two examples of utilities and their related losses.

##### Probability of Low Cost

For example, we may care about the system being below some cost threshold as
often as possible. To encode this preference, we can use the threshold utility
$U_{t}(C)\triangleq\mathbf{1}\mathopen{}\mathclose{{}\left\\{C\leq
C_{t,\max}}\right\\}$, where
$\mathbf{1}\mathopen{}\mathclose{{}\left\\{\cdot}\right\\}$ is the indicator
function and $C_{t,\max}$ is a threshold parameter. Under this choice, the
loss and its gradient become

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}}\right]$ (14)
$\displaystyle=-\log\mathbb{P}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left(C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right)$ $\displaystyle\nabla\ell_{t}(\bm{\theta})$
$\displaystyle=-\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})}\right]}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}}\right]}.$ (15)

As we can see, this loss function also gives the probability of achieving cost
below some threshold. As a result (Fig. 2(a)), costs below $C_{t,\max}$ are
treated the same in terms of the utility. This can potentially make
optimization easier since we are trying to make good trajectories as likely as
possible instead of finding the best trajectories as in (10).

However, if the threshold $C_{t,\max}$ is set too low and the gradient is
estimated with samples, the gradient estimate may have high variance due to
the large number of rejected samples. Because of this, in practice, the
threshold is set adaptively, e.g., as the largest cost of the top _elite
fraction_ of the sampled trajectories with smallest costs [6]. This allows the
controller to make the best sampled trajectories more likely and therefore
improve the controller.

##### Exponential Utility

We can also opt for a continuous surrogate of the indicator function, in this
case the exponential utility $U_{t}(C)\triangleq\exp(-\frac{1}{\lambda}C)$,
where $\lambda>0$ is a scaling parameter. Unlike the indicator function, the
exponential utility provides nonzero feedback for any given cost and allows us
to discriminate between costs (i.e., if $C_{1}>C_{2}$, then
$U_{t}(C_{1})<U_{t}(C_{2})$), as shown in Fig. 2(b). Furthermore, $\lambda$
acts as a continuous alternative to $C_{t,\max}$ and dictates how quickly or
slowly $U_{t}$ decays to zero, which in a soft way determines the cutoff point
for rejecting given costs.

Under this choice, the loss and its gradient become

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\lambda}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\right)}\right]$
(16) $\displaystyle\nabla\ell_{t}(\bm{\theta})$
$\displaystyle=-\frac{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\lambda}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\right)\nabla_{\bm{\theta}}\log\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})}\right]}{\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\lambda}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\right)}\right]}.$
(17)

The loss function in (16) is also known as the risk-seeking objective in
optimal control [28]; this classical interpretation is based on a Taylor
expansion of (16) showing

$\displaystyle\lambda\ell_{t}(\bm{\theta})\approx\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\right]-\frac{1}{\lambda}\mathbb{V}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})]$

when $\lambda$ is large, where
$\mathbb{V}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{f}}}[C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})]$
is the variance of $C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$. Here we derive (16)
from a different perspective that treats it as a continuous approximation of
(14). The use of exponential transformations to approximate indicators is a
common machine-learning trick (like the Chernoff bound [8]).

(a) Threshold utility

(b) Exponential utility

Figure 2: Visualization of different utilities.

###  III-B Algorithms

We instantiate DMD-MPC with different choices of loss function, control
distribution, and Bregman divergence as concrete examples to showcase the
flexibility of our framework. In particular, we are able to recover well-known
MPC algorithms as special cases of Algorithm 1.

Our discussions below are organized based on the class of Bregman divergences
used in (8), and the following algorithms are derived assuming that the
control distribution is a sequence of independent distributions. That is, we
suppose $\bm{\pi}_{\bm{\theta}}$ is a probability density/mass function that
factorizes as

$\displaystyle\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})=\prod_{h=0}^{H-1}\pi_{\theta_{h}}(\hat{u}_{t,h}),$
(18)

and $\bm{\theta}=(\theta_{0},\theta_{1},\ldots,\theta_{H-1})$ for some _basic
control distribution_ $\pi_{\theta}$ parameterized by $\theta\in\Theta$, where
$\Theta$ denotes the feasible set for the basic control distribution. For
control distributions in the form of (18), the shift operator $\Phi$ in (5)
would set $\tilde{\bm{\theta}}_{t}$ by identifying
$\tilde{\theta}_{t,h}=\theta_{t-1,h+1}$ for $h=0,\dots,H-2$, and initializing
the final parameter as either $\tilde{\theta}_{t,H-1}=\tilde{\theta}_{t,H-2}$
or $\tilde{\theta}_{t,H-1}=\bar{\theta}$ for some default parameter
$\bar{\theta}$.

####  III-B1 Quadratic Divergence

We start with perhaps the most common Bregman divergence: the quadratic
divergence. That is, we suppose the Bregman divergence in (8) has a quadratic
form 101010This is generated by defining
$\psi(\bm{\theta})\triangleq\frac{1}{2}\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\bm{A}\bm{\theta}$.
$D_{\psi}(\bm{\theta}\|\bm{\theta}^{\prime})\triangleq\frac{1}{2}(\bm{\theta}-\bm{\theta}^{\prime})^{\mkern-1.5mu\mathsf{T}}\bm{A}(\bm{\theta}-\bm{\theta}^{\prime})$
for some positive-definite matrix $\bm{A}$. Below we discuss different choices
of $\bm{A}$ and their corresponding update rules.

##### Projected Gradient Descent

This basic update rule is a special case when $\bm{A}$ is the identity matrix.
Equivalently, the update can be written as
$\bm{\theta}_{t}=\arg\min_{\bm{\theta}\in\bm{\Theta}}\,\|\bm{\theta}-(\tilde{\bm{\theta}}_{t}-\gamma_{t}\bm{g}_{t})\|^{2}$.

##### Natural Gradient Descent

We can recover the natural gradient descent algorithm [2] by defining
$\bm{A}=\mathcal{F}(\tilde{\bm{\theta}}_{t})$ where

$\mathcal{F}(\tilde{\bm{\theta}}_{t})=\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}}}\mathopen{}\mathclose{{}\left[\nabla_{\tilde{\bm{\theta}}_{t}}\log\bm{\pi}_{\tilde{\bm{\theta}}_{t}}(\hat{\bm{u}}_{t})\nabla_{\tilde{\bm{\theta}}_{t}}\log\bm{\pi}_{\tilde{\bm{\theta}}_{t}}(\hat{\bm{u}}_{t})^{\mkern-1.5mu\mathsf{T}}}\right]$

is the Fisher information matrix. This rule uses the natural Riemannian metric
of distributions to normalize the effects of different parameterizations of
the same distribution [25].

##### Quadratic Problems

While the above two update rules are quite general, we can further specialize
the Bregman divergence to achieve faster learning when the per-round loss
function can be shown to be quadratic. This happens, for instance, when the
MPC problem in (3) is an LQR or LEQR problem111111The dynamics model $\hat{f}$
is linear, the step cost $c$ is quadratic, the per-round loss $\ell_{t}$ is
(10), and the basic control distribution is a Dirac-delta distribution. [11].
That is, if

$\ell_{t}(\bm{\theta})=\frac{1}{2}\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\bm{R}_{t}\bm{\theta}+\bm{r}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{\theta}+\mathrm{const.}$

for some constant vector $\bm{r}_{t}$ and positive definite matrix
$\bm{R}_{t}$, we can set $\bm{A}=\bm{R}_{t}$ and $\gamma_{t}=1$, making
$\bm{\theta}_{t}$ given by the first step of (8) correspond to the optimal
solution to $\ell_{t}$ (i.e., the solution of LQR/LEQR). The particular values
of $\bm{R}_{t}$ and $\bm{r}_{t}$ for each of LQR and LEQR are derived in
Appendix D.

####  III-B2 KL Divergence and the Exponential Family

We show that for control distributions in the exponential family [23], the
Bregman divergence in (8) can be set to the KL divergence, which is a natural
way to measure distances between distributions. Toward this end, we review the
basics of the exponential family. We say a distribution $p_{\eta}$ with
natural parameter $\eta$ of random variable $u$ belongs to the _exponential
family_ if its probability density/mass function satisfies
$p_{\eta}(u)=\rho(u)\exp\mathopen{}\mathclose{{}\left(\langle\eta,\phi(u)\rangle-A(\eta)}\right)$,
where $\phi(u)$ is the sufficient statistics, $\rho(u)$ is the carrier
measure, and
$A(\eta)=\log\int\rho(u)\exp(\langle\eta,\phi(u)\rangle)\,\mathrm{d}u$ is the
log-partition function. The distribution $p_{\eta}$ can also be described by
its expectation parameter
$\mu\triangleq\mathbb{E}_{p_{\eta}}\mathopen{}\mathclose{{}\left[\phi(u)}\right]$,
and there is a duality between the two parameterizations: $\mu=\nabla
A(\eta)\text{ and }\eta=\nabla A^{*}(\mu)$, where
$A^{*}(\mu)=\sup_{\eta\in\mathcal{H}}\,\langle\eta,\mu\rangle-A(\eta)$ is the
Legendre transformation of $A$ and $\mathcal{H}=\\{\eta:A(\eta)<+\infty\\}$.
That is, $\nabla A=(\nabla A^{*})^{-1}$. The duality results in the property
below.

######  Fact 1.

[23]
$\mathrm{KL}(p_{\eta}\,\|\,p_{\eta^{\prime}})=D_{A}(\eta^{\prime}\|\eta)=D_{A^{*}}(\mu\|\mu^{\prime})$.

We can use 1 to define the Bregman divergence in (8) to optimize a control
distribution $\bm{\pi}_{\bm{\theta}}$ in the exponential family:

  * •

if $\bm{\theta}$ is an expectation parameter, we can set
$D_{\psi}(\bm{\theta}\|\tilde{\bm{\theta}}_{t})\triangleq\mathrm{KL}(\bm{\pi}_{\bm{\theta}}\,\|\,\bm{\pi}_{\tilde{\bm{\theta}}_{t}})$,
or

  * •

if $\bm{\theta}$ is a natural parameter, we can set
$D_{\psi}(\bm{\theta}\|\tilde{\bm{\theta}}_{t})\triangleq\mathrm{KL}(\bm{\pi}_{\tilde{\bm{\theta}}_{t}}\,\|\,\bm{\pi}_{\bm{\theta}})$.

We demonstrate some examples using this idea below.

##### Expectation Parameters and Categorical Distributions

We first discuss the case where $\bm{\theta}$ is an expectation parameter and
the first step in (8) is

$\displaystyle\bm{\theta}_{t}$
$\displaystyle=\arg\min_{\bm{\theta}\in\bm{\Theta}}\,\langle\gamma_{t}\bm{g}_{t},\bm{\theta}\rangle+\mathrm{KL}(\bm{\pi}_{\bm{\theta}}\,\|\,\bm{\pi}_{\tilde{\bm{\theta}}_{t}}).$
(19)

To illustrate, we consider an MPC problem with a _discrete_ control space
$\\{1,2,\ldots,m\\}$ and use the categorical distribution as the basic control
distribution in (18), i.e., we set
$\pi_{\theta_{h}}=\mathrm{Cat}(\theta_{h})$, where $\theta_{h}\in\Delta^{m}$
is the probability of choosing each control among $\\{1,2,\ldots,m\\}$ at the
$h$${}^{\text{th}}$ predicted time step and $\Delta^{m}$ denotes the
probability simplex in $\mathbb{R}^{m}$. This parameterization choice makes
$\bm{\theta}$ an expectation parameter of $\bm{\pi}_{\bm{\theta}}$ that
corresponds to sufficient statistics given by indicator functions. With the
structure of (9), the update direction is

$g_{t,h}=\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})e_{\hat{u}_{t,h}}\oslash\tilde{\theta}_{t,h}}\right]\quad(h=0,1,\dots,H-1)$

where $\tilde{\theta}_{t,h}$ and $g_{t,h}$ are the $h$${}^{\text{th}}$
elements of $\tilde{\bm{\theta}}_{t}$ and $\bm{g}_{t}$, respectively,
$e_{\hat{u}_{t,h}}\in\mathbb{R}^{m}$ has $0$ for each element except at index
$\hat{u}_{t,h}$ where it is $1$, and $\oslash$ denotes elementwise division.
Update (19) then becomes the exponentiated gradient algorithm [16]:

$\theta_{t,h}=\frac{1}{Z_{t,h}}\tilde{\theta}_{t,h}\odot\exp(-\gamma_{t}g_{t,h})\quad(h=0,1,\ldots,H-1)$
(20)

where $\theta_{t,h}$ is the $h$${}^{\text{th}}$ element of $\bm{\theta}_{t}$,
$Z_{t,h}$ is the normalizer for $\theta_{t,h}$, and $\odot$ denotes
elementwise multiplication. That is, instead of applying an additive gradient
step to the parameters, the update in (19) exponentiates the gradient and
performs elementwise multiplication. This does a better job of accounting for
the geometry of the problem, and makes projection a simple operation of
normalizing a distribution.

##### Natural Parameters and Gaussian Distributions

Alternatively, we can set $\bm{\theta}$ as a natural parameter and use

$\displaystyle\bm{\theta}_{t}$
$\displaystyle=\arg\min_{\bm{\theta}\in\bm{\Theta}}\,\langle\gamma_{t}\bm{g}_{t},\bm{\theta}\rangle+\mathrm{KL}(\bm{\pi}_{\tilde{\bm{\theta}}_{t}}\,\|\,\bm{\pi}_{\bm{\theta}})$
(21)

as the first step in (8). In particular, we show that, with (21), the
structure of the likelihood-ratio derivative in (9) can be leveraged to design
an efficient update. The main idea follows from the observation that when the
gradient is computed through (9) and $\tilde{\bm{\theta}}_{t}$ is the natural
parameter, we can write

$\displaystyle\bm{g}_{t}=\nabla\ell_{t}(\tilde{\bm{\theta}}_{t})=\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})(\phi(\hat{\bm{u}}_{t})-\tilde{\bm{\mu}}_{t})}\right]$
(22)

where $\tilde{\bm{\mu}}_{t}$ is the expectation parameter of
$\tilde{\bm{\theta}}_{t}$ and $\phi$ is the sufficient statistics of the
control distribution. We combine the factorization in (22) with a property of
the proximal update below (proven in Appendix C) to derive our algorithm.

######  Proposition 1.

Let $g_{t}$ be an update direction. Let $\mathcal{M}$ be the image of
$\mathcal{H}$ under $\nabla A$. If $\mu_{t}-\gamma_{t}g_{t}\in\mathcal{M}$ and
$\eta_{t+1}=\arg\min_{\eta\in\mathcal{H}}\,\langle\gamma_{t}g_{t},\eta\rangle+D_{A}(\eta\|\eta_{t})$,
then $\mu_{t+1}=\mu_{t}-\gamma_{t}g_{t}$.121212A similar proposition can be
found for (19).

We find that, under the assumption131313If $\mu_{t}-\gamma_{t}g_{t}$ is not in
$\mathcal{M}$, the update in (21) needs to perform a projection, the form of
which is algorithm dependent. in 1, the update rule in (21) becomes

$\displaystyle\bm{\mu}_{t+1}=(1-\gamma_{t})\tilde{\bm{\mu}}_{t}+\gamma_{t}\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\phi(\hat{\bm{u}}_{t})}\right].$
(23)

In other words, when $\gamma_{t}\in[0,1]$, the update to the expectation
parameter $\bm{\mu}_{t}$ in (8) is simply a convex combination of the
sufficient statistics and the previous expectation parameter
$\tilde{\bm{\mu}}_{t}$.

We provide a concrete example of an MPC algorithm that follows from (23). Let
us consider a continuous control space and use the Gaussian distribution as
the basic control distribution in (18), i.e., we set
$\pi_{\theta_{h}}(\hat{u}_{t,h})=\mathcal{N}(\hat{u}_{t,h};m_{h},\Sigma_{h})$
for some mean vector $m_{h}$ and covariance matrix $\Sigma_{h}$. For
$\pi_{\theta_{h}}$, we can choose sufficient statistics
$\phi(\hat{u}_{t,h})=(\hat{u}_{t,h},\hat{u}_{t,h}\hat{u}_{t,h}^{\mkern-1.5mu\mathsf{T}})$,
which results in the expectation parameter $\mu_{h}=(m_{h},S_{h})$ and the
natural parameter
$\eta_{h}=(\Sigma_{h}^{-1}m_{h},-\frac{1}{2}\Sigma_{h}^{-1})$, where
$S_{h}\triangleq\Sigma_{h}+m_{h}m_{h}^{\mkern-1.5mu\mathsf{T}}$ is the second
moment of $\pi_{\theta_{h}}$. Let us set $\theta_{h}$ as the natural
parameter. Then (21) is equivalent to the update rule for $h=0,\dots,H-1$:

$\displaystyle\begin{split}m_{t,h}&=(1-\gamma_{t})\tilde{m}_{t,h}+\gamma_{t}\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\hat{u}_{t,h}}\right]\\\
S_{t,h}&=(1-\gamma_{t})\tilde{S}_{t,h}+\gamma_{t}\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[L_{t}(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\hat{u}_{t,h}\hat{u}_{t,h}^{\mkern-1.5mu\mathsf{T}}}\right].\end{split}$
(24)

Several existing algorithms are special cases of (24).

  * •

Cross-entropy method (CEM) [6]:
If $\ell_{t}$ is set to (14) and $\gamma_{t}=1$, then (24) becomes

$\displaystyle\begin{split}m_{t,h}&=\frac{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}\hat{u}_{t,h}}\right]}{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}}\right]}\\\
S_{t,h}&=\frac{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}\hat{u}_{t,h}\hat{u}_{t,h}^{\mkern-1.5mu\mathsf{T}}}\right]}{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[\mathbf{1}\mathopen{}\mathclose{{}\left\\{C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})\leq
C_{t,\max}}\right\\}}\right]},\end{split}$ (25)

which matches the update rule of the cross-entropy method for Gaussian
distributions [6].141414Though CEM is typically presented as updating the mean
and _covariance_ of a Gaussian distribution, the update rule is derived by
matching the first and second moments between the Gaussian distribution and a
uniform distribution over trajectories whose costs are at most $C_{t,\max}$,
which is identical to (25).

  * •

Model-predictive path integral (MPPI) [31]:

If we choose $\ell_{t}$ as the exponential utility, as in (16), and do not
update the covariance, the update rule becomes

$\displaystyle
m_{t,h}=(1-\gamma_{t})\tilde{m}_{t,h}+\gamma_{t}\frac{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[e^{-\frac{1}{\lambda}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}\hat{u}_{t,h}}\right]}{\mathbb{E}_{\bm{\pi}_{\tilde{\bm{\theta}}_{t}},\hat{\bm{f}}}\mathopen{}\mathclose{{}\left[e^{-\frac{1}{\lambda}C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})}}\right]},$
(26)

which reduces to the MPPI update rule [31] for $\gamma_{t}=1$. This connection
is also noted in [24].

###  III-C Extensions

In the previous sections, we discussed multiple instantiations of DMD-MPC,
showing the flexibility of our framework. But they are by no means exhaustive.
In Appendix B, we discuss variations of DMD-MPC, e.g., imposing constraints
and different ways to approximate the expectation in (9).

##  IV Related Work

Recent work on MPC has studied sampling-based approaches, which are flexible
in that they do not require differentiability of a cost function. One such
algorithm which can be used with general cost functions and dynamics is MPPI,
which was proposed by Williams et al. [31] as a generalization of the control
affine case [30]. The algorithm is derived by considering an optimal control
distribution defined by the control problem. This optimal distribution is
intractable to sample from, so the algorithm instead tries to bring a
tractable distribution (in this case, Gaussian with fixed covariance) as close
as possible in the sense of KL divergence. This ends up being the same as
finding the mean of the optimal control distribution. The mean is then
approximated as a weighted sum of sampled control trajectories, where the
weight is determined by the exponentiated costs. Although this algorithm works
well in practice (including a robust variant [33] achieving state-of-the-art
performance in aggressive driving [10]), it is not clear that matching the
mean of the distribution should guarantee good performance, such as in the
case of a multimodal optimal distribution. By contrast, our update rule in
(26) results from optimizing an exponential utility.

A closely related approach is the cross-entropy method (CEM) [6], which also
assumes a Gaussian sampling distribution but minimizes the KL divergence
between the Gaussian distribution and a uniform distribution over low cost
samples. CEM has found applicability in reinforcement learning [19, 21, 26],
motion planning [17, 18], and MPC [9, 32].

These sampling-based control algorithms can be considered special cases of
general derivative-free optimization algorithms, such as covariance matrix
adaptation evolutionary strategies (CMA-ES) [15] and natural evolutionary
strategies (NES) [29]. CMA-ES samples points from a multivariate Gaussian,
evaluates their fitness, and adapts the mean and covariance of the sampling
distribution accordingly. On the other hand, NES optimizes the parameters of
the sampling distribution to maximize some expected fitness through steepest
ascent, where the direction is provided by the natural gradient. Akimoto et
al. [2] showed that CMA-ES can also be interpreted as taking a natural
gradient step on the parameters of the sampling distribution. As we showed in
Section III-B, natural gradient descent is a special case of DMD-MPC
framework. A similar observation that connects between MPPI and mirror descent
was made by Okada and Taniguchi [24], but their derivation is limited to the
KL divergence and Gaussian case.

##  V Experiments

We use experiments to the validate the flexibility of DMD-MPC. We show that
this framework can handle both continuous (Gaussian distribution) and discrete
(categorial distribution) variations of control problems, and that MPC
algorithms like MPPI and CEM can be generalized using different step sizes and
control distributions to improve performance. Extra details and results are
included in Appendices E and F.

###  V-A Cartpole

We first consider the classic cartpole problem where we seek to swing a pole
upright and keep it balanced only using actuation on the attached cart. We
consider both the continuous and discrete control variants. For the continuous
case, we choose the Gaussian distribution as the control distribution and keep
the covariance fixed. For the discrete case, we choose the categorical
distribution and use update (20). In either case, we have access to a biased
stochastic model (uses a different pole length compared to the real cart).

(a) Continuous controls

(b) Discrete controls

Figure 3: Varying step size and number of samples (same legends for (a) and
(b)). EC = expected cost (10). PLC = probability of low cost (14) with elite
fraction $=10^{-3}$. EU = exponential utility (16) with $\lambda=1$.

(a) Continuous controls

(b) Discrete controls

Figure 4: Varying loss parameter and step size (1000 samples).

We consider the interaction between the choice of loss, step size, and number
of samples used to estimate (9),151515For our experiments, we vary the number
of samples from $\bm{\pi}_{\bm{\theta}}$ and fix the number of samples from
$\hat{\bm{f}}$ to ten. Furthermore, we use common random numbers when sampling
from $\hat{\bm{f}}$ to reduce estimation variance. shown in Figs. 3 and 4. For
this environment, we can achieve low cost when optimizing the expected cost in
(10) with a proper step size ($10^{-2}$ for both continuous and discrete
problems) while being fairly robust to the number of samples. When using
either of the utilities, the number of samples is more crucial in the
continuous domain, with more samples allowing for larger step sizes. In the
discrete domain (Fig. 3(b)), performance is largely unaffected by the number
of samples when the step size is below $10$, excluding the threshold utility
with 1000 samples. In Fig. 4(a), for a large range of utility parameters, we
see that using step sizes above $1$ (the step size set in MPPI and CEM) give
significant performance gains. In Fig. 4(b), there’s a more complicated
interaction between the utility parameter and step size, with huge changes in
cost when altering the utility parameter and keeping the step size fixed.

###  V-B AutoRally

####  V-B1 Platform Description

We use the autonomous AutoRally platform [13] to run a high-speed driving task
on a dirt track, with the goal of the task to achieve as low a lap time as
possible. The robot (Fig. 6) is a 1:5 scale RC chassis capable of driving over
$20~{}\mathrm{m/s}$ ($45~{}\mathrm{mph}$) and has a desktop-class Intel Core
i7 CPU and Nvidia GTX 1050 Ti GPU. Our code for the control algorithm is based
on modifications of code available on the AutoRally
repository.161616https://github.com/AutoRally/autorally For real-world
experiments, we estimate the car’s pose using a particle filter from [10]
which relies on a monocular camera, IMU, and GPS. In both simulated and real-
world experiments, the dynamics model is a neural network which has been
fitted to data collected from human demonstrations. We note that the dynamics
model is deterministic, so we don’t need to estimate any expectations with
respect to the dynamics.

Figure 5: Simulated AutoRally performance with different step sizes and number
of samples. Though many samples coupled with large steps yield the smallest
lap times, the performance gains are small past 1920 samples. With fewer
samples, a lower step size helps recover some lost performance.

####  V-B2 Simulated Experiments

We first use the Gazebo simulator (Fig. 9 in Section E-B) from the AutoRally
repo to perform a sweep of algorithm parameters, particularly the step size
and number of samples, to evaluate how changing these parameters can affect
the performance of DMD-MPC. For all of the experiments, the control
distribution is a Gaussian with fixed covariance, and we use update (26)
(i.e., the loss is the exponential utility (16)) with $\lambda=6.67$. The
resulting lap times are shown in Fig. 5.171717The large error bar for 64
samples and step size of 0.8 is due to one particular lap where the car
stalled at a turn for about 60 seconds. We see that although using more
samples does result in smaller lap times, there are diminishing returns past
1920 samples per gradient. Indeed, with a proper step size, even as few as 192
samples can yield lap times within a couple seconds of 3840 samples and a step
size of 1. We also observe that the curves converge as the step size decreases
further, implying that only a certain number of samples are needed for a given
step size. This is a particularly important advantage of DMD-MPC over methods
like MPPI: by changing the step size, DMD-MPC can perform much more
effectively with fewer samples, making it a good choice for embedded systems
which can’t produce many samples due to computational constraints.

Figure 6: Rally car driving during an experiment. Figure 7: Real-world AutoRally task. TABLE I: Statistics for real-world experiments at target of $9~{}\mathrm{m/s}$. Samples | Step size $\gamma_{t}$ | Lap time ($\mathrm{s}$) | Avg. speed ($\mathrm{m/s}$) | Max speed ($\mathrm{m/s}$)
---|---|---|---|---
$1920$ | $1$ | $31.76\pm 0.55$ | $5.70\pm 0.16$ | $9.21\pm 0.30$
| $0.8$ | $31.81\pm 0.21$ | $5.75\pm 0.03$ | $9.03\pm 0.19$
| $0.6$ | $32.83\pm 0.31$ | $5.60\pm 0.05$ | $8.62\pm 0.12$
$64$ | $1$ | $33.74\pm 0.78$ | $5.45\pm 0.16$ | $9.50\pm 0.22$
| $0.8$ | $33.84\pm 0.80$ | $5.46\pm 0.11$ | $9.12\pm 0.26$
| $0.6$ | $33.61\pm 0.74$ | $5.50\pm 0.13$ | $9.14\pm 0.42$
TABLE II: Statistics for real-world experiments at target of $11~{}\mathrm{m/s}$. Samples | Step size $\gamma_{t}$ | Lap time ($\mathrm{s}$) | Avg. speed ($\mathrm{m/s}$) | Max speed ($\mathrm{m/s}$)
---|---|---|---|---
$64$ | $1$ | $31.05\pm 0.67$ | $5.80\pm 0.26$ | $10.17\pm 0.30$
| $0.6$ | $30.30\pm 0.56$ | $5.98\pm 0.15$ | $10.30\pm 0.05$

####  V-B3 Real-World Experiments

In the real-world setting (Fig. 7), the control distribution is a Gaussian
with fixed covariance, and we use update (26) with $\lambda=8$. We ran two
sets of experiments, each with a different target speed: one at
$9~{}\mathrm{m/s}$ and the other at $11~{}\mathrm{m/s}$.181818The conference
version of this paper does not have the second set of experiments (target of
$11~{}\mathrm{m/s}$). Those experiments were conducted after the camera-ready
deadline of the conference.

For the first set of experiments, we used the following configurations: each
of 1920 and 64 samples, and each of step sizes 1 (corresponding to MPPI), 0.8,
and 0.6.191919Due to weaker batteries used with 64 samples, results should not
be compared across number of samples. Overall (Table I), there’s a mild
degradation in performance when decreasing the step size at 1920 samples, due
to the car taking a longer path on the track (Fig. 12(a) vs. Fig. 12(c) in
Section F-B). With 64 samples, the results seem unaffected by the step size.
This could be because, despite the noisiness of the DMD-MPC update, the
setpoint controller in the car’s steering servo acts as a filter, smoothing
out the control signal and allowing the car to drive on a consistent path
(Fig. 13 in Section F-B). Videos of this experiment can be found at
https://youtu.be/vZST3v0_S9w.

For the second set of experiments, we fixed the number of samples at 64 and
used step sizes of 1 (corresponding to MPPI) and 0.6. The statistics slightly
improve with a decreased step size, but qualitatively there is a larger
difference between the step sizes. With a step size of 1, the car often
wobbles while driving, turns around at one point, and crashes in one of the
trials (Fig. 14(a)). On the other hand, with a step size of 0.6, the car
drives much more smoothly and achieves the aggressive driving task with no
issues (Fig. 14(b)). Despite the smoothing effect of the low-level controllers
in the car, the more stringent costs associated with the larger target speed
cause the noisiness of the DMD-MPC update to manifest in the car’s performance
when using a step size of 1. A smaller step size mitigates this noisiness.
Videos of this experiment can be found at https://youtu.be/MhuqiHo2t98.

##  VI Conclusion

We presented a connection between model predictive control and online
learning. From this connection, we proposed an algorithm based on dynamic
mirror descent that can work for a wide variety of settings and cost
functions. We also discussed the choice of loss function within this online
learning framework and the sort of preference each loss function imposes. From
this general algorithm and assortment of loss functions, we show several well
known algorithms are special cases and presented a general update for members
of the exponential family.

We empirically validated our algorithm on continuous and discrete simulated
problems and on a real-world aggressive driving task. In the process, we also
studied the parameter choices within the framework, finding, for example, that
in our framework a smaller number of rollout samples can be compensated for by
varying other parameters like the step size.

We hope that the online learning and stochastic optimization viewpoints of MPC
presented in this paper opens up new possibilities for using tools from these
domains, such as alternative efficient sampling techniques [5] and accelerated
optimization methods [22, 24], to derive new MPC algorithms that perform well
in practice.

## Acknowledgements

This material is based upon work supported by NSF NRI award 1637758, NSF
CAREER award 1750483, an NSF Graduate Research Fellowship under award No.
2015207631, and a National Defense Science & Engineering Graduate Fellowship.
We thank Aravind Battaje, Nathan Hatch, and Hemanth Sarabu for assisting in
AutoRally experiments.

## References

  * Abbeel et al. [2010] Pieter Abbeel, Adam Coates, and Andrew Y. Ng.  Autonomous helicopter aerobatics through apprenticeship learning.  _The International Journal of Robotics Research_ , 29(13):1608–1639, 2010.
  * Akimoto et al. [2012] Youhei Akimoto, Yuichi Nagata, Isao Ono, and Shigenobu Kobayashi.  Theoretical foundation for CMA-ES from information geometry perspective.  _Algorithmica_ , 64(4):698–716, 2012.
  * Banerjee et al. [2005] Arindam Banerjee, Srujana Merugu, Inderjit S. Dhillon, and Joydeep Ghosh.  Clustering with Bregman divergences.  _Journal of machine learning research_ , 6(Oct):1705–1749, 2005.
  * Beck and Teboulle [2003] Amir Beck and Marc Teboulle.  Mirror descent and nonlinear projected subgradient methods for convex optimization.  _Operations Research Letters_ , 31(3):167–175, 2003.
  * Bellman and Casti [1971] Richard Bellman and John Casti.  Differential quadrature and long-term integration.  _Journal of Mathematical Analysis and Applications_ , 34(2):235–238, 1971.
  * Botev et al. [2013] Zdravko I. Botev, Dirk P. Kroese, Reuven Y. Rubinstein, and Pierre L’Ecuyer.  The cross-entropy method for optimization.  In _Handbook of statistics_ , volume 31, pages 35–59. Elsevier, 2013\.
  * Camacho and Alba [2013] Eduardo F. Camacho and Carlos Bordons Alba.  _Model predictive control_.  Springer Science & Business Media, 2013.
  * Chernoff [1952] Herman Chernoff.  A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations.  _The Annals of Mathematical Statistics_ , 23(4):493–507, 1952.
  * Chua et al. [2018] Kurtland Chua, Roberto Calandra, Rowan McAllister, and Sergey Levine.  Deep reinforcement learning in a handful of trials using probabilistic dynamics models.  In _Advances in Neural Information Processing Systems_ , pages 4754–4765, 2018.
  * Drews et al. [2019] Paul Drews, Grady Williams, Brian Goldfain, Evangelos A. Theodorou, and James M. Rehg.  Vision-Based High-Speed Driving With a Deep Dynamic Observer.  _IEEE Robotics and Automation Letters_ , 4(2):1564–1571, 2019.
  * Duncan [2013] Tyrone E. Duncan.  Linear-exponential-quadratic Gaussian control.  _IEEE Transactions on Automatic Control_ , 58(11):2910–2911, 2013.
  * Glynn [1990] Peter W. Glynn.  Likelihood ratio gradient estimation for stochastic systems.  _Communications of the ACM_ , 33(10):75–84, 1990\.
  * Goldfain et al. [2019] Brian Goldfain, Paul Drews, Changxi You, Matthew Barulic, Orlin Velev, Panagiotis Tsiotras, and James M. Rehg.  AutoRally: An Open Platform for Aggressive Autonomous Driving.  _IEEE Control Systems Magazine_ , 39(1):26–55, 2019.
  * Hall and Willett [2013] Eric Hall and Rebecca Willett.  Dynamical models and tracking regret in online convex programming.  In _International Conference on Machine Learning_ , pages 579–587, 2013.
  * Hansen et al. [2003] Nikolaus Hansen, Sibylle D. Müller, and Petros Koumoutsakos.  Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES).  _Evolutionary Computation_ , 11(1):1–18, 2003\.
  * Hazan [2016] Elad Hazan.  Introduction to online convex optimization.  _Foundations and Trends® in Optimization_ , 2(3-4):157–325, 2016.
  * Helvik and Wittner [2002] Bjarne E. Helvik and Otto Wittner.  Using the Cross-Entropy Method to Guide/Govern Mobile Agent’s Path Finding in Networks.  _International Workshop on Mobile Agents for Telecommunication Applications_ , 2164, 2002.
  * Kobilarov [2012] Marin Kobilarov.  Cross-entropy randomized motion planning.  In _Robotics: Science and Systems_ , volume 7, pages 153–160, 2012\.
  * Mannor et al. [2003] Shie Mannor, Reuven Y. Rubinstein, and Yohai Gat.  The cross entropy method for fast policy search.  In _Proceedings of the 20th International Conference on Machine Learning (ICML-03)_ , pages 512–519, 2003.
  * Mayne [2014] David Q. Mayne.  Model predictive control: Recent developments and future promise.  _Automatica_ , 50(12):2967–2986, 2014.
  * Menache et al. [2005] Ishai Menache, Shie Mannor, and Nahum Shimkin.  Basis Function Adaptation in Temporal Difference Reinforcement Learning.  _Annals of Operations Research_ , 134:215–238, 2005.
  * Miyashita et al. [2018] Megumi Miyashita, Shiro Yano, and Toshiyuki Kondo.  Mirror descent search and its acceleration.  _Robotics and Autonomous Systems_ , 106:107–116, 2018.
  * Nielsen and Garcia [2009] Frank Nielsen and Vincent Garcia.  Statistical exponential families: A digest with flash cards.  _arXiv preprint arXiv:0911.4863_ , 2009.
  * Okada and Taniguchi [2018] Masashi Okada and Tadahiro Taniguchi.  Acceleration of Gradient-based Path Integral Method for Efficient Optimal and Inverse Optimal Control.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 3013–3020. IEEE, 2018.
  * Rattray et al. [1998] Magnus Rattray, David Saad, and Shun-ichi Amari.  Natural gradient descent for on-line learning.  _Physical review letters_ , 81(24):5461, 1998\.
  * Szita and Lörincz [2006] István Szita and András Lörincz.  Learning Tetris using the noisy cross-entropy method.  _Neural computation_ , 18(12):2936–2941, 2006\.
  * Tassa et al. [2014] Yuval Tassa, Nicolas Mansard, and Emo Todorov.  Control-limited differential dynamic programming.  In _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1168–1175. IEEE, 2014.
  * van den Broek et al. [2010] LJ van den Broek, WAJJ Wiegerinck, and HJ Kappen.  Risk sensitive path integral control.  In _Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence (UAI 2010)_ , pages 1–8. AUAI Press, 2010.
  * Wierstra et al. [2014] Daan Wierstra, Tom Schaul, Tobias Glasmachers, Yi Sun, Jan Peters, and Jürgen Schmidhuber.  Natural evolution strategies.  _The Journal of Machine Learning Research_ , 15(1):949–980, 2014.
  * Williams et al. [2016] Grady Williams, Paul Drews, Brian Goldfain, James M. Rehg, and Evangelos A. Theodorou.  Aggressive driving with model predictive path integral control.  In _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1433–1440. IEEE, 2016.
  * Williams et al. [2017] Grady Williams, Nolan Wagener, Brian Goldfain, Paul Drews, James M. Rehg, Byron Boots, and Evangelos A. Theodorou.  Information Theoretic MPC for Model-Based Reinforcement Learning.  In _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1714–1721. IEEE, 2017.
  * Williams et al. [2018a] Grady Williams, Paul Drews, Brian Goldfain, James M. Rehg, and Evangelos A. Theodorou.  Information-Theoretic Model Predictive Control: Theory and Applications to Autonomous Driving.  _IEEE Transactions on Robotics_ , 34(6):1603–1622, 2018a.
  * Williams et al. [2018b] Grady Williams, Brian Goldfain, Paul Drews, Kamil Saigol, James M. Rehg, and Evangelos A. Theodorou.  Robust sampling based model predictive control with sparse objective information.  In _Robotics Science and Systems_ , 2018b.
  * Zhang et al. [2018] Lijun Zhang, Shiyin Lu, and Zhi-Hua Zhou.  Adaptive Online Learning in Dynamic Environments.  In _Advances in Neural Information Processing Systems_ , pages 1323–1333, 2018.

##  Appendix A Shift Operator

We discuss some details in defining the shift operator. Let
$\bm{\theta}_{t-1}$ be the approximate solution to the previous problem and
$\tilde{\bm{\theta}}_{t}$ denote the initial condition of $\bm{\theta}$ in
solving (3), and consider sampling
$\hat{\bm{u}}_{t}\sim\bm{\pi}_{\tilde{\bm{\theta}}_{t}}$ and
$\hat{\bm{u}}_{t-1}\sim\bm{\pi}_{\bm{\theta}_{t-1}}$. We set

$\tilde{\bm{\theta}}_{t}=\Phi(\bm{\theta}_{t-1})$

by defining a _shift operator_ $\Phi$ that outputs a new parameter in
$\bm{\Theta}$. This $\Phi$ can be chosen to satisfy desired properties, one
example being that when conditioned on $\hat{u}_{t-1}$ and $x_{t}$, the
marginal distributions of $\hat{u}_{t},\dots,\hat{u}_{t+H-2}$ are the same for
both $\hat{\bm{u}}_{t}$ of $\bm{\pi}_{\tilde{\bm{\theta}}_{t}}$ and
$\hat{\bm{u}}_{t-1}$ of $\bm{\pi}_{\bm{\theta}_{t-1}}$. A simple example of
this property is shown in Fig. 8. Note that $\hat{\bm{u}}_{t}$ also involves a
new control $\hat{u}_{t+H-1}$ that is not in $\hat{\bm{u}}_{t-1}$, so the
choice of $\Phi$ is not unique but algorithm dependent; for example, we can
set $\hat{u}_{t+H-1}$ of $\bm{\pi}_{\tilde{\bm{\theta}}_{t}}$ to follow the
same distribution as $\hat{u}_{t+H-2}$ (cf. Section III-B). Because the
subproblems in (3) of two consecutive time steps share all control variables
except for the first and the last ones, the “shifted” parameter
$\Phi(\bm{\theta}_{t-1})$ to the current problem should be almost as good as
the optimized parameter $\bm{\theta}_{t-1}$ is to the previous problem. In
other words, setting $\tilde{\bm{\theta}}_{t}=\Phi(\bm{\theta}_{t-1})$
provides a warm start to (3) and amortizes the computational complexity of
solving for $\bm{\theta}_{t}$.

Figure 8: A simple example of the shift operator $\Phi$. Here, the control
distribution $\bm{\pi}_{\bm{\theta}}$ consists of a sequence of $H=5$
independent Gaussian distributions. The shift operator moves the parameters of
the Gaussians one time step forward and replaces the parameters at $h=4$ with
some default parameters.

##  Appendix B Variations of DMD-MPC

The control distributions in DMD-MPC can be fairly general (in addition to the
categorical and Gaussian distributions that we discussed) and control
constraints on the problem (e.g., control limits) can be directly incorporated
through proper choices of control distributions, such as the beta
distribution, or through mapping the unconstrained control through some
squashing function (e.g., $\tanh$ or clamp). Though our framework cannot
directly handle state constraints as in constrained optimization approaches, a
constraint can be relaxed to an indicator function which activates if the
constraint is violated. The indicator function can then be added to the cost
function in (4) with some weight that encodes how strictly the constraint
should be enforced.

Moreover, different integration techniques, such as Gaussian quadrature [5],
can be adopted to replace the likelihood-ratio derivative in (9) for computing
the required gradient direction. We also note that the independence assumption
on the control distribution in (18) is not necessary in our framework; time-
correlated control distributions and feedback policies are straightforward to
consider in DMD-MPC.

##  Appendix C Proofs

###### Proof of 1.

We prove the first statement; the second one follows directly from the duality
relationship. The statement follows from the derivations below; we can write

$\displaystyle\eta_{t+1}$
$\displaystyle=\arg\min_{\eta\in\mathcal{H}}\,\langle\gamma_{t}g_{t},\eta\rangle+D_{A}(\eta\|\eta_{t})$
$\displaystyle=\arg\min_{\eta\in\mathcal{H}}\,\langle\gamma_{t}g_{t},\eta\rangle+A(\eta)-\langle\nabla
A(\eta_{t}),\eta\rangle$
$\displaystyle=\arg\min_{\eta\in\mathcal{H}}\,\langle\gamma_{t}g_{t}-\mu_{t},\eta\rangle+A(\eta)$
$\displaystyle=\arg\max_{\eta\in\mathcal{H}}\,\langle\mu_{t}-\gamma_{t}g_{t},\eta\rangle-A(\eta)$
$\displaystyle=\nabla A^{*}(\mu_{t}-\gamma_{t}g_{t})$

where the last equality is due to the assumption that
$\mu_{t}-\gamma_{t}g_{t}\in\mathcal{M}$. Then applying $\nabla A$ on both
sides and using the relationship that $\nabla A=(\nabla A^{*})^{-1}$, we have
$\mu_{t+1}=\nabla A(\eta_{t+1})=\mu_{t}-\gamma_{t}g_{t}$. ∎

##  Appendix D Derivation of LQR and LEQR Losses

The dynamics in Equation (1) are given by

$x_{t+1}=Ax_{t}+Bu_{t}+w_{t}$

for some matrices $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times
m}$ and $w_{t}\sim\mathcal{N}(0,W)$, where $W\in\mathbb{S}_{++}^{n}$. For a
control sequence $\hat{\bm{u}}_{t}$, noise sequence $\hat{\bm{w}}_{t}$, and
initial state $x_{t}$, the resulting state sequence $\hat{\bm{x}}_{t}$ is
found through convolution:

$\begin{bmatrix}\hat{x}_{t}\\\ \hat{x}_{t+1}\\\ \hat{x}_{t+2}\\\ \vdots\\\
\hat{x}_{t+H}\end{bmatrix}=\begin{bmatrix}I\\\ A\\\ A^{2}\\\ \vdots\\\
A^{H}\end{bmatrix}x_{t}+\begin{bmatrix}0&0&\cdots&0\\\ B&0&\cdots&0\\\
AB&B&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
A^{H-1}B&A^{H-2}B&\cdots&B\end{bmatrix}\begin{bmatrix}\hat{u}_{t}\\\
\hat{u}_{t+1}\\\ \vdots\\\
\hat{u}_{t+H-1}\end{bmatrix}+\begin{bmatrix}0&0&\cdots&0\\\ I&0&\cdots&0\\\
A&I&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
A^{H-1}&A^{H-2}&\cdots&I\end{bmatrix}\begin{bmatrix}\hat{w}_{t}\\\
\hat{w}_{t+1}\\\ \vdots\\\ \hat{w}_{t+H-1}\end{bmatrix},$

or, in matrix form:

$\hat{\bm{x}}_{t}=\bm{F}x_{t}+\bm{G}\hat{\bm{u}}_{t}+\bm{L}\hat{\bm{w}}_{t},$

where $\bm{F}$, $\bm{G}$, and $\bm{L}$ are defined naturally from the
convolution equation above. Note that
$\hat{\bm{w}}_{t}\sim\mathcal{N}(0,\bm{W})$, where
$\bm{W}=\mathrm{diag}(W,W,\ldots,W,W)$. Thus, we also have that

$\hat{\bm{x}}_{t}\sim\mathcal{N}(\bm{F}x_{t}+\bm{G}\hat{\bm{u}}_{t},\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}).$

We define the instantaneous and terminal costs as

$\displaystyle c(x,u)$
$\displaystyle=\frac{1}{2}x^{\mkern-1.5mu\mathsf{T}}Qx+\frac{1}{2}u^{\mkern-1.5mu\mathsf{T}}Ru$
$\displaystyle c_{\mathrm{end}}(x)$
$\displaystyle=\frac{1}{2}x^{\mkern-1.5mu\mathsf{T}}Q_{\mathrm{end}}x,$

where $Q,Q_{\mathrm{end}}\in\mathbb{S}_{+}^{n}$ and $R\in\mathbb{S}_{++}^{m}$.
Thus, the statistic $C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})$ is

$C(\hat{\bm{x}}_{t},\hat{\bm{u}}_{t})=\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\hat{\bm{x}}_{t}+\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}\hat{\bm{u}}_{t},$

where $\bm{Q}=\mathrm{diag}(Q,Q,\ldots,Q,Q_{\mathrm{end}})$ and
$\bm{R}=\mathrm{diag}(R,R,\ldots,R,R)$.

Our control distribution is a Dirac delta distribution located at the given
parameter:
$\bm{\pi}_{\bm{\theta}}(\hat{\bm{u}}_{t})=\delta(\hat{\bm{u}}_{t}-\bm{\theta})$.

###  D-A LQR

The loss is defined as
$\ell_{t}(\bm{\theta})=\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{x}}_{t}}\mathopen{}\mathclose{{}\left[\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\hat{\bm{x}}_{t}+\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}\hat{\bm{u}}_{t}}\right]$.
Expanding this out gives:

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{x}}_{t}}\mathopen{}\mathclose{{}\left[\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\hat{\bm{x}}_{t}+\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}\hat{\bm{u}}_{t}}\right]$
$\displaystyle=\frac{1}{2}\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\mathopen{}\mathclose{{}\left(\bm{G}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{G}+\bm{R}}\right)\bm{\theta}+x_{t}^{\mkern-1.5mu\mathsf{T}}\bm{F}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{G}\bm{\theta}+\frac{1}{2}x_{t}^{\mkern-1.5mu\mathsf{T}}\bm{F}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{F}x_{t}+\frac{1}{2}\mathbb{E}\mathopen{}\mathclose{{}\left[\hat{\bm{w}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{L}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{L}\hat{\bm{w}}_{t}}\right]$
$\displaystyle=\frac{1}{2}\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\mathopen{}\mathclose{{}\left(\bm{G}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{G}+\bm{R}}\right)\bm{\theta}+x_{t}^{\mkern-1.5mu\mathsf{T}}\bm{F}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{G}\bm{\theta}+\frac{1}{2}x_{t}^{\mkern-1.5mu\mathsf{T}}\bm{F}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{F}x_{t}+\frac{1}{2}\mathrm{tr}\mathopen{}\mathclose{{}\left(\bm{Q}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}}\right).$

We see this is a quadratic problem in $\bm{\theta}$ by defining

$\displaystyle\bm{R}_{t}$
$\displaystyle=\bm{G}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{G}+\bm{R}$
$\displaystyle\bm{r}_{t}$
$\displaystyle=\bm{G}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{F}x_{t}.$

###  D-B LEQR

The loss is defined as

$\ell_{t}(\bm{\theta})=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{x}}_{t}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\lambda}\mathopen{}\mathclose{{}\left(\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\hat{\bm{x}}_{t}+\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}\hat{\bm{u}}_{t}}\right)}\right)}\right]$

for some parameter $\lambda>0$. For compactness, we define
$\bm{Q}^{\prime}=\frac{1}{\lambda}\bm{Q}$ and
$\bm{R}^{\prime}=\frac{1}{\lambda}\bm{R}$ so that the exponent contains
$-\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}^{\prime}\hat{\bm{x}}-\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}^{\prime}\hat{\bm{u}}_{t}$.
In expanding the loss, we use the following fact:

######  Fact 2.

For $x\sim\mathcal{N}(\mu,\Sigma)$, where $\Sigma\in\mathbb{S}_{++}^{n}$, and
constants $A\in\mathbb{S}_{+}^{n}$ and $b\in\mathbb{R}^{n}$:

$\mathbb{E}_{x}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}x^{\mkern-1.5mu\mathsf{T}}Ax-b^{\mkern-1.5mu\mathsf{T}}x}\right)}\right]=\frac{1}{\sqrt{|A\Sigma+I|}}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\mathopen{}\mathclose{{}\left(\mu^{\mkern-1.5mu\mathsf{T}}\Sigma^{-1}\mu-(\Sigma^{-1}\mu-b)^{\mkern-1.5mu\mathsf{T}}(A+\Sigma^{-1})^{-1}(\Sigma^{-1}\mu-b)}\right)}\right).$

###### Proof.

We expand the expectation and complete the square:

$\displaystyle\mathbb{E}_{x}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}x^{\mkern-1.5mu\mathsf{T}}Ax-b^{\mkern-1.5mu\mathsf{T}}x}\right)}\right]$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{n}|\Sigma|}}\int\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}(x-\mu)^{\mkern-1.5mu\mathsf{T}}\Sigma^{-1}(x-\mu)}\right)\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}x^{\mkern-1.5mu\mathsf{T}}Ax-b^{\mkern-1.5mu\mathsf{T}}x}\right)\,\mathrm{d}x$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{n}|\Sigma|}}\int\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\mathopen{}\mathclose{{}\left[x^{\mkern-1.5mu\mathsf{T}}(A+\Sigma^{-1})x+2(b-\Sigma^{-1}\mu)^{\mkern-1.5mu\mathsf{T}}x+\mu^{\mkern-1.5mu\mathsf{T}}\Sigma^{-1}\mu}\right]}\right)\,\mathrm{d}x$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{n}|\Sigma|}}\exp(c)\int\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}(x-\tilde{\mu})^{\mkern-1.5mu\mathsf{T}}\tilde{\Sigma}^{-1}(x-\tilde{\mu})}\right)\,\mathrm{d}x$
$\displaystyle=\frac{\sqrt{(2\pi)^{n}|\tilde{\Sigma}|}}{\sqrt{(2\pi)^{n}|\Sigma|}}\exp(c)$
$\displaystyle=\frac{1}{\sqrt{|A+\Sigma^{-1}||\Sigma|}}\exp(c)$
$\displaystyle=\frac{1}{\sqrt{|A\Sigma+I|}}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\mathopen{}\mathclose{{}\left(\mu^{\mkern-1.5mu\mathsf{T}}\Sigma^{-1}\mu-(\Sigma^{-1}\mu-b)^{\mkern-1.5mu\mathsf{T}}(A+\Sigma^{-1})^{-1}(\Sigma^{-1}\mu-b)}\right)}\right),$

where $\tilde{\mu}=(A+\Sigma^{-1})^{-1}(\Sigma^{-1}\mu-b)$,
$\tilde{\Sigma}=(A+\Sigma^{-1})^{-1}$, and
$c=-\frac{1}{2}\mathopen{}\mathclose{{}\left(\mu^{\mkern-1.5mu\mathsf{T}}\Sigma^{-1}\mu-(\Sigma^{-1}\mu-b)^{\mkern-1.5mu\mathsf{T}}(A+\Sigma^{-1})^{-1}(\Sigma^{-1}\mu-b)}\right)$.
∎

We now expand the loss:

$\displaystyle\ell_{t}(\bm{\theta})$
$\displaystyle=-\log\mathbb{E}_{\bm{\pi}_{\bm{\theta}},\hat{\bm{x}}_{t}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}^{\prime}\hat{\bm{x}}_{t}-\frac{1}{2}\hat{\bm{u}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{R}^{\prime}\hat{\bm{u}}_{t}}\right)}\right]$
$\displaystyle=-\log\mathbb{E}_{\hat{\bm{x}}_{t}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\hat{\bm{x}}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{Q}^{\prime}\hat{\bm{x}}_{t}-\frac{1}{2}\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\bm{R}^{\prime}\bm{\theta}}\right)}\right]$
$\displaystyle=-\log\Bigg{\\{}\frac{1}{\sqrt{|\bm{Q}^{\prime}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+I|}}\exp\bigg{(}-\frac{1}{2}\big{[}(\bm{F}x_{t}+\bm{G}\bm{\theta})^{\mkern-1.5mu\mathsf{T}}(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}(\bm{F}x_{t}+\bm{G}\bm{\theta})$
$\displaystyle\qquad~{}\qquad~{}\qquad~{}\qquad~{}\qquad~{}\qquad\quad\qquad~{}\;-(\bm{F}x_{t}+\bm{G}\bm{\theta})^{\mkern-1.5mu\mathsf{T}}(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}\bm{Q}^{\prime}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}(\bm{F}x_{t}+\bm{G}\bm{\theta})$
$\displaystyle\qquad~{}\qquad~{}\qquad~{}\qquad~{}\qquad~{}\qquad\quad\qquad~{}\;+\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\bm{R}^{\prime}\bm{\theta}\big{]}\bigg{)}\Bigg{\\}}$
$\displaystyle=\frac{1}{2}\mathopen{}\mathclose{{}\left[(\bm{F}x_{t}+\bm{G}\bm{\theta})^{\mkern-1.5mu\mathsf{T}}[(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}+(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}\bm{Q}^{\prime}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}](\bm{F}x_{t}+\bm{G}\bm{\theta})+\bm{\theta}^{\mkern-1.5mu\mathsf{T}}\bm{R}^{\prime}\bm{\theta}}\right]+\frac{1}{2}\log|\bm{Q}^{\prime}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+I|.$

We see this is a quadratic problem in $\bm{\theta}$ by defining

$\displaystyle\bm{R}_{t}$
$\displaystyle=\bm{G}^{\mkern-1.5mu\mathsf{T}}\mathopen{}\mathclose{{}\left[(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}+\mathopen{}\mathclose{{}\left(\frac{1}{\lambda}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}}\right)^{-1}}\right]\bm{G}+\frac{1}{\lambda}\bm{R}$
$\displaystyle\bm{r}_{t}$
$\displaystyle=\bm{G}^{\mkern-1.5mu\mathsf{T}}\mathopen{}\mathclose{{}\left[(\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}})^{-1}+\mathopen{}\mathclose{{}\left(\frac{1}{\lambda}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}\bm{Q}\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}+\bm{L}\bm{W}\bm{L}^{\mkern-1.5mu\mathsf{T}}}\right)^{-1}}\right]\bm{F}x_{t}.$

##  Appendix E Experimental Setup

###  E-A Cartpole

The state is $x_{t}=(p_{t},\varphi_{t},v_{t},\dot{\varphi}_{t})$, where
$p_{t}$ is the cart position, $\varphi_{t}$ is the pole’s angle, $v_{t}$ and
$\dot{\varphi}_{t}$ are the corresponding velocities, and the control $u_{t}$
is the force applied to the cart. We define the instantaneous cost and
terminal cost of the MPC problem as

$\displaystyle c(x_{t},u_{t})$
$\displaystyle=10p_{t}^{2}+500(\varphi_{t}-\pi)^{2}+v_{t}^{2}+15\dot{\varphi}_{t}^{2}+1000\cdot\mathbf{1}\mathopen{}\mathclose{{}\left\\{|\varphi_{t}-\pi|\geq\Delta}\right\\}$
$\displaystyle c_{\mathrm{end}}(x_{t})$ $\displaystyle=c(x_{t},0)$

where $\Delta$ is some threshold. For our experiments, we set
$\Delta=12^{\circ}=0.21$ radians.

In our experiments, the pole is massless except for some weight at the end of
the pole. The mass of the cart and pole weight are $0.711~{}\mathrm{kg}$ and
$0.209~{}\mathrm{kg}$, respectively. The true length of the pole is
$0.326~{}\mathrm{m}$, whereas the length used in the model is
$0.346~{}\mathrm{m}$. Each time step is modeled using an Euler discretization
of $0.02$ seconds. Each episode of the problem lasts 500 time steps (i.e, 10
seconds) and has episode cost equal to the sum of encountered instantaneous
costs. Both the true system and the model apply Gaussian additive noise to the
commanded control with zero mean and a standard deviation of $5$ newtons. For
the continuous system, the commanded control is clamped to $\pm 25$ newtons.
For the discrete system, the controller can either command $10$ newtons to the
left, $10$ newtons to the right, or $0$ newtons.

Both the discrete and continuous controller use a planning horizon of 50 time
steps (i.e., 1 second). For the continuous controller, we keep the standard
deviation of the Gaussian distribution fixed at $2$ newtons for each time step
in the planning horizon. When applying a control $u_{t}$ on the real cartpole,
we choose the mode of $\pi_{\theta_{t}}$ rather than sample from the
distribution.

All reported results were gathered using ten episodes per parameter setting.

###  E-B AutoRally

The state of the vehicle is
$x_{t}=(p_{x,t},p_{y,t},\varphi_{t},r_{t},v_{x,t},v_{y,t},\dot{\varphi}_{t})$,
where $(p_{x,t},p_{y,t})$ is the position of the car in the global frame,
$\varphi_{t}$ and $r_{t}$ are the yaw and roll angles, $v_{x,t}$ and $v_{y,t}$
are the longitudinal and lateral velocities in the car frame, and
$\dot{\varphi}_{t}$ is the yaw rate. The control $u_{t}$ we apply is the
throttle and steering angle. For some weights $w_{1},\ldots,w_{4}$, the cost
function is

$\displaystyle c(x_{t},u_{t})$
$\displaystyle=w_{1}|s_{t}-s_{\mathrm{tgt}}|^{k}+w_{2}M(p_{x,t},p_{y,t})+w_{3}S_{c}(x_{t})$
$\displaystyle c_{\mathrm{end}}(\bm{x}_{t})$
$\displaystyle=w_{4}C(\bm{x}_{t}).$

Here, $s_{t}$ and $s_{\mathrm{tgt}}$ are the current and target speed of the
car, respectively. Note the speed is calculated as
$s_{t}=\sqrt{v_{x,t}^{2}+v_{y,t}^{2}}$. $M(p_{x,t},p_{y,t})$ is the positional
cost of the car (low cost in center of track, high cost at edge of track),
$S_{c}(x_{t})$ is an indicator variable which activates if the slip
angle202020The slip angle is defined as $-\arctan\frac{v_{y,t}}{|v_{x,t}|}$,
which gives the angle between the direction the car is pointing and the
direction in which it is actually traveling. exceeds a certain threshold, and
$C(\bm{x}_{t})$ is an indicator function which activates if the car leaves the
track at all in the trajectory. Note that the terminal cost depends on the
trajectory instead of the terminal state. Each time step represents $0.02$
seconds for every experiment except the real-world experiment with a target of
$11~{}\mathrm{m/s}$ where each time step represents $0.025$ seconds. The
length of the planning trajectory is 100 time steps (i.e., either 2 seconds or
2.5 seconds depending on the length of the time step). The values for the cost
function parameters are given in Table III.

The control space for each of the throttle and steering angle is normalized to
the range $[-1,1]$. For our experiments, we clamp the throttle to $[-1,0.65]$.
In simulated experiments, the standard deviations of the throttle and steering
angle distributions were $0.3$ and $0.275$, respectively. In the real world
experiments, they were both set to $0.3$. When applying a control $u_{t}$ on
the car, we chose the mean of $\pi_{\theta_{t}}$ rather than sampling from the
distribution.

Figure 9: Simulated AutoRally task.

In simulation, the environment (Fig. 9) is an elliptical track approximately 3
meters wide and 30 meters across at its furthest point. The real-world dirt
track is about 5 meters wide and and has a track length of 170 meters. All
reported results for simulated experiments were gathered using 30 consecutive
laps in the counter-clockwise direction for each parameter setting. For real-
world experiments, results were gathered using ten laps for each parameter
setting when the target speed is $9~{}\mathrm{m/s}$ and five laps for
$11~{}\mathrm{m/s}$.

TABLE III: Cost function settings for AutoRally experiments. |  $s_{\mathrm{tgt}}$ ($\mathrm{m/s}$) | $k$ | $w_{1}$ | $w_{2}$ | $w_{3}$ | $w_{4}$ |  Slip angle threshold ($\mathrm{rad}$)
---|---|---|---|---|---|---|---
Gazebo simulator | $11$ | $1$ | $30$ | $250$ | $10$ | $10000$ | $0.275$
Real world |  $9$ or $11$ | $2$ | $4.25$ | $200$ | $100$ | $10000$ | $0.9$

##  Appendix F Extra Experimental Results

###  F-A Simulated Experiments

Adding onto the results from Section V-B2, we qualitatively evaluate two
particular extremes: few vs. many samples (64 vs. 3840) and small vs. large
step size (0.5 vs. 1) by looking at the path and speed of the car during the
episode (Fig. 10). At small step sizes (Figs. 10(a) and 10(c)), the path and
speed profiles are rather similar, while with few samples and a large step
size (Fig. 10(b)), the car drives much more slowly and erratically, sometimes
even stopping. In the ideal scenario with many samples and a large step size,
the car can achieve consistently high speed while driving smoothly (Fig.
10(d)).

We also experimented with instead optimizing the expected cost (10) and found
performance was dramatically worse (Fig. 11), even when using 3840 samples per
gradient. At best, the car would drive in the center of the track at speeds
below $4~{}\mathrm{m/s}$ (Fig. 11(c)), and at worst, the car would either
slowly drive along the track walls (Fig. 11(a)) or the controller would
eventually produce $\mathrm{NaN}$ controls that would prematurely end the
experiment (Fig. 11(d)). This poor performance is likely due to most samples
in the estimate of (11) having very high cost (e.g., due to leaving the track)
and contributing significantly to the gradient estimate. On the other hand,
when estimating (17), as in the experiments in Section V-B2, these high cost
trajectories are assigned very low weights so that only low cost trajectories
contribute to the gradient estimate.

(a) 64 samples, $\gamma_{t}=0.5$

(b) 64 samples, $\gamma_{t}=1$

(c) 3840 samples, $\gamma_{t}=0.5$

(d) 3840 samples, $\gamma_{t}=1$

Figure 10: Car speeds when optimizing the exponential utility (16). The speeds
and trajectories are very similar at step size $0.5$, irrespective of the
number of samples. At step size $1$, though, 64 samples result in capricious
maneuvers and low speeds, whereas 3840 samples result in smooth driving at
high speeds.

(a) $\gamma_{t}=0.025$

(b) $\gamma_{t}=0.05$

(c) $\gamma_{t}=0.075$

(d) $\gamma_{t}=0.1$

Figure 11: Car speeds when optimizing the expected cost (10). All tested step
sizes result in low speeds. At too low or too high of a step size, the car
will drive along the wall or crash into it.

###  F-B Figures for Real-World Experiments

(a) $\gamma_{t}=1$

(b) $\gamma_{t}=0.8$

(c) $\gamma_{t}=0.6$

Figure 12: Car speeds with 1920 samples per gradient estimate and target of
$9~{}\mathrm{m/s}$.

(a) $\gamma_{t}=1$

(b) $\gamma_{t}=0.8$

(c) $\gamma_{t}=0.6$

Figure 13: Car speeds with 64 samples per gradient estimate and target of
$9~{}\mathrm{m/s}$.

(a) $\gamma_{t}=1$

(b) $\gamma_{t}=0.6$

Figure 14: Car speeds with 64 samples per gradient estimate and target of
$11~{}\mathrm{m/s}$. In Fig. 14(a), note the crash and U-turn at the top of
the plot as well as the wider spread of the paths throughout the whole track.
By contrast, in Fig. 14(b), the resulting paths are more consistent, and there
are no failure points.