An Online Learning Approach to Model Predictive Control

TL;DR

This paper introduces DMD-MPC, a novel online learning-based framework for model predictive control that unifies and extends existing MPC techniques, demonstrated through simulations and real-world experiments.

Contribution

It presents DMD-MPC, a new algorithm linking MPC with online learning, enabling principled design and extension of MPC methods.

Findings

DMD-MPC generalizes many existing MPC algorithms.

The approach is effective in simulated cartpole and real-world driving tasks.

Videos demonstrate practical applicability of the proposed methods.

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.

Tables

Table 1. TABLE I : Statistics for real-world experiments at target of 9 ​ m / s 9 m s 9~{}\mathrm{m/s} .

Samples	Step size ${\gamma }_t$	Lap time ($s$)	Avg. speed ($m/s$)	Max speed ($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 2. TABLE II : Statistics for real-world experiments at target of 11 ​ m / s 11 m s 11~{}\mathrm{m/s} .

Samples	Step size ${\gamma }_t$	Lap time ($s$)	Avg. speed ($m/s$)	Max speed ($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$

Table 3. TABLE III : Cost function settings for AutoRally experiments.

	$s_{\mathrm{tgt}}$ ($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$

Equations

$$x_{t+1}\sim f(x_t,u_t)$$

$${\hat{\mathbf{x}}}_t\sim \hat{\mathbf{f}}(x_t,{\hat{\mathbf{u}}}_t)$$

$$\underset{\mathbf{\theta }\in \mathbf{\Theta }}{\min }\hat{J}({\mathbf{\pi }}_{\mathbf{\theta }};x_t).$$

$$C({\hat{\mathbf{x}}}_t,{\hat{\mathbf{u}}}_t)\triangleq \sum _{h=0}^{H-1}c({\hat{x}}_{t+h},{\hat{u}}_{t+h})+c_{\mathrm{end}}({\hat{x}}_{t+H}),$$

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

$${\ell }_t(\cdot )=\hat{J}(\cdot ;x_t).$$

$$\mathrm{D}\text{-}\mathrm{Regret}=\sum _{t=1}^T{\ell }_t({\tilde{\mathbf{\theta }}}_t)-\sum _{t=1}^T{\ell }_t({\mathbf{\theta }}_t^{\star }),$$

$$\begin{array}{rl}{\mathbf{\theta }}_t& =\arg \underset{\mathbf{\theta }\in \mathbf{\Theta }}{\min }⟨{\gamma }_t{\mathbf{g}}_t,\mathbf{\theta }⟩+D_{\psi }(\mathbf{\theta }\Vert {\tilde{\mathbf{\theta }}}_t)\\ {\tilde{\mathbf{\theta }}}_{t+1}& =\Phi ({\mathbf{\theta }}_t),\end{array}$$

$$\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].$$

$${\ell }_t(\mathbf{\theta })$$

$$\nabla {\ell }_t(\mathbf{\theta })$$

$${\ell }_t(\mathbf{\theta })$$

$$\nabla {\ell }_t(\mathbf{\theta })$$

$$\nabla {\ell }_t(\mathbf{\theta })\approx -\sum _{i=1}^N{w}_i{\nabla }_{\mathbf{\theta }}\log {\mathbf{\pi }}_{\mathbf{\theta }}({\hat{\mathbf{u}}}_t^i),$$

$${\ell }_t(\mathbf{\theta })$$

$$\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)$$

$$\nabla {\ell }_t(\mathbf{\theta })$$

$${\ell }_t(\mathbf{\theta })$$

$$\nabla {\ell }_t(\mathbf{\theta })$$

$$\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})]$$

$${\mathbf{\pi }}_{\mathbf{\theta }}({\hat{\mathbf{u}}}_t)=\prod _{h=0}^{H-1}{\pi }_{{\theta }_h}({\hat{u}}_{t,h}),$$

$$\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]$$

$${\ell }_t(\mathbf{\theta })=\frac{1}{2}{\mathbf{\theta }}^{\mathsf{T}}{\mathbf{R}}_t\mathbf{\theta }+{\mathbf{r}}_t^{\mathsf{T}}\mathbf{\theta }+const.$$

$${\mathbf{\theta }}_t$$

$$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)$$

$${\theta }_{t,h}=\frac{1}{Z_{t,h}}{\tilde{\theta }}_{t,h}\odot \exp (-{\gamma }_t{g}_{t,h})(h=0,1,\dots ,H-1)$$

$${\mathbf{\theta }}_t$$

$$\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]$$

$$\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].$$

$$\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}$$

$$\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}$$

$$\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]},$$

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

Full text

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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.

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:

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

[TABLE]

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.

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]).

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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]:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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

[TABLE]

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).

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.

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.

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.

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

[TABLE]

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}$.

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

[TABLE]

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

[TABLE]

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:

[TABLE]

or, in matrix form:

[TABLE]

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

[TABLE]

We define the instantaneous and terminal costs as

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

D-B LEQR

The loss is defined as

[TABLE]

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}$:

[TABLE]

Proof.

We expand the expectation and complete the square:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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

[TABLE]

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.

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}$.

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.

F-B Figures for Real-World Experiments