Localized active learning of Gaussian process state space models

Autonomous systems often need to operate in complex environments, of which a model is difficult or even impossible to derive from first principles. Learning-based techniques have become a promising paradigm to address these issues (Pillonetto et al., 2014). In particular, Gaussian processes (GPs) have been increasingly employed for system identification and control (Umlauft et al., 2018; Capone and Hirche, 2019; Berkenkamp and Schoellig, 2015; Deisenroth et al., 2015). GPs possess very good generalization properties (Rasmussen and Williams, 2006), which can be leveraged to obtain data-efficient learning-based approaches (Deisenroth et al., 2015; Kamthe and Deisenroth, 2018). By employing a Bayesian framework, GPs provide an automatic trade-off between model smoothness and data fitness. Moreover, GPs provide an explicit estimate of the model uncertainty that is used to derive probabilistic bounds in control settings (Capone and Hirche, 2019; Beckers et al., 2019; Umlauft and Hirche, 2020).

A performance-determining factor of data-driven techniques is the quality of the available data. In settings where data is insufficient to achieve accurate predictions, new data needs to be gathered via exploration (Umlauft and Hirche, 2020). In classical reinforcement learning settings, exploration is often enforced by randomly selecting a control action with a predetermined probability that tends to zero over time (Dayan and Sejnowski, 1996). However, this is generally inefficient, as regions of low uncertainty are potentially revisited in multiple iterations. These issues have been addressed by techniques that choose the most informative exploration trajectories (Alpcan and Shames, 2015; Ay et al., 2008; Burgard et al., 2005; Schreiter et al., 2015). The goal of these methods is to obtain a globally accurate model. While this is reasonable for systems with a bounded state-action space, it is unsuited for systems with unbounded ones, particularly if a non-parametric model is used. This is because a potentially infinite number of points is required to achieve a globally accurate model. Furthermore, in practice a model often only needs to be accurate locally, e.g., for stabilization tasks.

In this paper, we propose a model predictive control-based exploration approach that steers the system towards the most informative points within a bounded subset of the state-action space. By modeling the system with a Gaussian process, we are able to quantify the information inherent in each data point. Our approach chooses actions by approximating the mutual information of the system trajectory with respect to a discretization of the region of interest. This is achieved by first selecting the single most informative data point within the region of interest, then steering the system towards that point using model predictive control. Through this approximation, the solution approach is rendered computationally tractable.

The remainder of this paper is structured as follows. Section 2 describes the general problem. Section 3 discusses how GPs are employed for modeling and exploration. Section 4 presents the LocAL algorithm, and is followed by a numerical simulation example, in Section 5.

We consider the problem of exploring the state and control space of a discrete-time nonlinear system with Markovian dynamics of the form

where $t\in\mathbb{N}_{0}$, ${{\bm{x}}_{t}\in\mathcal{X}\subseteq\mathbb{R}^{{d_{x}}}}$ and ${{\bm{u}}_{t}\in\mathcal{U}\subseteq\mathbb{R}^{{d_{u}}}}$ are the system’s state vector and control vector at the $t$-th time step, respectively. The system is disturbed by multivariate Gaussian process noise ${{\bm{w}}_{t}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{w})}$ with $\bm{\Sigma}_{w}=\text{diag}(\sigma_{\text{w},1}^{2},\ldots,\sigma_{\text{w},{d_{x}}}^{2})$, $\sigma_{\text{w,i}}^{2}\in\mathbb{R}_{+,0}$. The concatenation ${{\tilde{\bm{x}}}_{t}:=\begin{pmatrix}{\bm{x}}_{t}^{\text{T}}&{\bm{u}}_{t}^{\text{T}}\end{pmatrix}^{\text{T}}\in{\tilde{\mathcal{X}}}}$, where ${{\tilde{\mathcal{X}}}:=\mathcal{X}\times\mathcal{U}}$ is employed for simplicity of exposition. The nonlinear function ${{\bm{f}}:{\tilde{\mathcal{X}}}\rightarrow\mathcal{X}}$ represents the known component of the system dynamics, e.g., a model obtained using first principles, while ${{\bm{g}}:{\tilde{\mathcal{X}}}\rightarrow\mathcal{X}}$ corresponds to the unknown component of the system dynamics.

We aim to obtain an approximation of the function ${\bm{g}}(\cdot)$, denoted $\widehat{\bm{g}}(\cdot)$, which provides an accurate estimate of ${{\bm{g}}(\cdot)}$ on a predefined bounded subset of the augmented state space ${\tilde{\mathcal{X}}}_{B}\subset{\tilde{\mathcal{X}}}$. This is often required in practice, e.g., for local stabilization tasks.

In order to faithfully capture the stochastic behavior of (1), we model the system as a Gaussian process (GP), where we employ measurements of the augmented state vector ${\tilde{\bm{x}}}_{t}$ as training inputs, and the differences ${\bm{x}}_{t+1}-{\bm{f}}({\tilde{\bm{x}}}_{t})={\bm{g}}({\tilde{\bm{x}}}_{t})+{\bm{w}}_{t}$ as training targets.

A GP is a collection of dependent random variables variables, for which any finite subset is jointly normally distributed (Rasmussen and Williams, 2006). It is specified by a mean function $m:{\tilde{\mathcal{X}}}\rightarrow\mathbb{R}$ and a positive definite covariance function $k:{\tilde{\mathcal{X}}}\times{\tilde{\mathcal{X}}}\rightarrow\mathbb{R}$, also known as kernel. In this paper, we set $m\equiv 0$ without loss of generality, as all prior knowledge is already encoded in $f(\cdot)$. The kernel $k(\cdot,\cdot)$ is a similarity measure for evaluations of ${\bm{g}}(\cdot)$, and encodes function properties such as smoothness and periodicity.

In the case where the state is a scalar, i.e., ${d_{x}}=1$, given $n$ training input samples ${\tilde{\bm{X}}}=\left\{{\tilde{\bm{x}}}_{1},\ldots,{\tilde{\bm{x}}}_{n}\right\}\subset{\tilde{\mathcal{X}}}$ and training outputs $\bm{y}_{{\tilde{\bm{X}}}}=\begin{pmatrix}g({\tilde{\bm{x}}}_{1})+{w}_{1}&\ldots&g({\tilde{\bm{x}}}_{n})+{w}_{n}\end{pmatrix}^{\text{T}}$, the posterior mean and variance of the GP corresponds to a one-step transition model. Starting at a point ${\tilde{\bm{x}}}_{t}$, the difference between the subsequent state and the known component is normally distributed, i.e.,

with mean and variance given by

respectively, where ${\bm{k}(\cdot)=\begin{pmatrix}k({\tilde{\bm{x}}}_{1},\cdot)&\ldots&k({\tilde{\bm{x}}}_{n},\cdot)\end{pmatrix}^{\text{T}}}$, and the entries of the covariance matrix $\bm{K}$ are computed as $K_{ij}=k({\tilde{\bm{x}}}_{i},{\tilde{\bm{x}}}_{j})$, $i,j=1,\ldots,n$.

In the case where the state is multidimensional, we model dimension of the state transition function using a separate GP. This corresponds to the assumption that the state transition function entries are conditionally independent. For simplicity of exposition, unless stated otherwise, we henceforth assume ${d_{x}}=1$. However, the methods presented in this paper extend straightforwardly to the multivariate case.

The system dynamics (1) considerably limit the decision space at every time step $t$. Furthermore, after a data point is collected, both the GP model and mutual information change. Hence, we employ a model predictive control (MPC)-based approach to steer the system towards areas of high information. Ideally, at every MPC-step $t$, we would like to minimize (3) with respect to a series of $N_{H}$ inputs $\bm{U}\coloneqq\{{\bm{u}}_{t},\ldots,{\bm{u}}_{t+N_{H}-1}\}$. However, this is generally infeasible, limiting its applicability in an MPC setting. Hence, we consider an approximate solution approach that sequentially computes the most informative data point by minimizing (4) separately from the MPC optimization. This is achieved as follows. At every time step $t$, an unconstrained minimizer $\bm{\xi}^{*}$ of (4) is computed. Afterwards, the MPC computes the approximate optimal inputs $\bm{U}^{*}$ by minimizing a constrained optimization problem that penalizes the weighted distance to the reference point $\bm{\xi}$ using a positive semi-definite weight matrix $\bm{Q}$. The ensuing state is then measured, the GP model is updated, and the procedure is repeated. These steps yield the Localized Active Learning (LocAL) algorithm, which is presented in Algorithm 1.

The square weight matrix $\bm{Q}\in\mathbb{R}^{{d_{x}}+{d_{u}}}\times\mathbb{R}^{{d_{x}}+{d_{u}}}$ should be chosen such that the MPC steers the system as close to $\bm{\xi}$ as possible. This represents a system-dependent task. Alternatively, $\bm{Q}$ can be chosen such that the MPC cost function corresponds to a quadratic approximation of the mutual information, e.g., such that $I(\bm{y}_{{\tilde{\bm{x}}}_{t}},\bm{y}_{\bm{\xi}})\approx\sum_{t=1}^{N_{H}}\left(\bm{\xi}-{\tilde{\bm{x}}}_{t}\right)^{\text{T}}\bm{Q}\left(\bm{\xi}-{\tilde{\bm{x}}}_{t}\right)$ holds for ${\tilde{\bm{x}}}_{t}\approx\bm{\xi}$.

In this section, we apply the proposed algorithm to four different dynamical systems. We begin with a toy example, with which we can easily illustrate the explored portions of the state space. Afterwards, we apply the proposed approach to a pendulum, a cart-pole, and a synthetic model that generalizes the mountain car problem. The exploration is repeated $50$ times for each system using different starting points $\bm{x}_{0}$ sampled from a normal distribution. To quantify the performance of each approach, we compute the root mean square model error (RMSE) on $500$ points sampled from a uniform distribution on the region of interest ${{\tilde{\mathcal{X}}}_{\text{ref}}}$.

We employ a squared-exponential kernel in all examples, and train the hyperparemters online using gradient-based log likelihood maximization (Rasmussen and Williams, 2006). We employ an MPC horizon of $N_{H}=10$, and choose weight matrix for the MPC optimization step as

where $\sigma_{g_{d}}$ denotes the standard deviation of the GP kernel corresponding to the $d$-th dimension, and $l^{-2}_{1,g_{d}},\ldots,l^{-2}_{{d_{x}},g_{d}}$ denote the corresponding lengthscales. In order to ease the computational burden, we apply the first $7$ inputs computed by the LocAL algorithm before computing a new solution.

We additionally explore each system using a greedy entropy-based cost function, as suggested in Koller et al. (2018) and Schreiter et al. (2015), and compare the results. In all three cases, the LocAL algorithm yields a better model in the regions of interest than the entropy-based algorithm.

A technique for efficiently exploring bounded subsets of the state-action space of a system has been presented. The proposed technique aims to minimize the mutual information of the system trajectories with respect to a discretization of the region of interest. It employs Gaussian processes both to model the unknown system dynamics and to quantify the informativeness of potentially collected data points. In numerical simulations of four different dynamical systems, we have demonstrated that the proposed approach yields a better model after a limited amount of time steps than a greedy entropy-based approach.

This work was supported by the European Research Council (ERC) Consolidator Grant ”Safe data-driven control for human-centric systems (CO-MAN)” under grant agreement number 864686 .