Pseudospectral Model Predictive Control under Partially Learned Dynamics

TL;DR

This paper introduces a semi-parametric model predictive control method that combines physics-based models with Gaussian Processes for improved trajectory optimization under uncertain dynamics, demonstrated on robotics simulations.

Contribution

It presents a novel semi-parametric modeling approach integrated with pseudospectral control and online learning using Sparse Spectrum Gaussian Processes for robotic systems.

Findings

Effective handling of unmodeled damping and parametric errors

Successful online learning and obstacle avoidance in simulations

Improved trajectory optimization with semi-parametric models

Abstract

Trajectory optimization of a controlled dynamical system is an essential part of autonomy, however many trajectory optimization techniques are limited by the fidelity of the underlying parametric model. In the field of robotics, a lack of model knowledge can be overcome with machine learning techniques, utilizing measurements to build a dynamical model from the data. This paper aims to take the middle ground between these two approaches by introducing a semi-parametric representation of the underlying system dynamics. Our goal is to leverage the considerable information contained in a traditional physics based model and combine it with a data-driven, non-parametric regression technique known as a Gaussian Process. Integrating this semi-parametric model with model predictive pseudospectral control, we demonstrate this technique on both a cart pole and quadrotor simulation with unmodeled damping and parametric error. In order to manage parametric uncertainty, we introduce an algorithm that utilizes Sparse Spectrum Gaussian Processes (SSGP) for online learning after each rollout. We implement this online learning technique on a cart pole and quadrator, then demonstrate the use of online learning and obstacle avoidance for the dubin vehicle dynamics.

Tables

Table 1. Table 1 : Gaussian Process Example Problem Settings

Problem	Parameters	True	Nominal	Cost Function
Cart Pole	Cart Mass $m_c$ (kg)	0.5	0.5
	Link Mass $m_p$ (kg)	0.5	0.5	${\int }_{t_0}^{tf}\left(R{𝐮}^2\right)𝑑t$
	Link Length $l$ (m)	0.6	0.6
	Ground damping $b_1$ (N/s)	0.01	0
	Pivot damping $b_2$ (Nm/s)	0.01	0
Quadrotor	Quadrotor Mass $m$ (kg)	1	.9
	Inertias $J_x$ (kg m2)	8.1E-3	8.1E-3
	$J_y$ (kg m2)	8.1E-3	8.1E-3	${\int }_{t_0}^{tf}\left({𝐱}^{T}Q𝐱+{𝐮}^{T}R𝐮\right)𝑑t$
	$J_z$ (kg m2)	14.2E-3	14.2E-3
	Arm Length $L$ (m)	.24	.24

Table 2. Table 2 : Online Learning Example Problem Settings

Problem	Parameters	True	Nominal	Cost Function
Cart Pole	Cart Mass $m_c$ (kg)	1.0	.9
	Link Mass $m_p$ (kg)	0.5	0.45	${\int }_{t_0}^{tf}\left(R{𝐮}^2\right)𝑑t$
	Link Length $l$ (m)	0.5	0.45
	Ground damping $b_1$ (N/s)	0.3	0
	Pivot damping $b_2$ (Nm/s)	0.3	0
Quadrotor	Quadrotor Mass $m$ (kg)	1	.99
	Inertias $J_x$ (kg m2)	8.1E-3	8.1E-3
	$J_y$ (kg m2)	8.1E-3	8.1E-3	${\int }_{t_0}^{tf}\left({𝐱}^{T}Q𝐱+{𝐮}^{T}R𝐮\right)𝑑t$
	$J_z$ (kg m2)	14.2E-3	14.2E-3
	Arm Length $L$ (m)	.24	.2376

Equations

$$\displaystyle\underset{\mathbf{u}}{\text{min}}~{}J(\mathbf{x}(t),\mathbf{u}(t))=\underset{\mathbf{u}}{\text{min}}~{}\bigg{(}\underbrace{\psi(\mathbf{x}(t_{f}),t_{f})}_{\text{Cost to go}}+\underbrace{\int_{t_{0}}^{t_{f}}\mathcal{L}(\mathbf{x}(t),\mathbf{u}(t),t)dt}_{\text{Running Cost}}\bigg{)}$$

$$\text{subject to}\dot{\mathbf{x}}(t)=f_t(\mathbf{x},\mathbf{u},t)$$

$$0\ge h_s(\mathbf{x})$$

$$0\ge h_c(\mathbf{u})$$

$$f_t(\mathbf{x}(t),\mathbf{u}(t))$$

$$\displaystyle p~{}\bigg{(}\begin{matrix}\mathbf{~{}Y~{}}\\ \mathbf{~{}y}^{*}\end{matrix}\bigg{)}\sim\mathcal{N}\bigg{(}~{}0,\bigg{[}\begin{matrix}\mathbf{K}(\mathbf{Z},\mathbf{Z})+\sigma_{n}^{2}\mathbf{I}&\mathbf{K}(\mathbf{Z},\mathbf{z}^{*})\\ \mathbf{K}(\mathbf{z}^{*},\mathbf{Z})&\mathbf{K}(\mathbf{z}^{*},\mathbf{z}^{*})\end{matrix}\bigg{]}~{}\bigg{)}$$

$$\mathbf{K}({\mathbf{z}}_i,{\mathbf{z}}_j)={\sigma }_s^2\exp (-\frac{1}{2}({\mathbf{z}}_i-{\mathbf{z}}_j{)}^T\mathbf{W}({\mathbf{z}}_i-{\mathbf{z}}_j))+{\sigma }_n^2$$

$${\mathbf{y}}^{\ast }|{\mathbf{z}}^{\ast },\mathbf{Z},\mathbf{Y}$$

$$\mathbf{K}({\mathbf{z}}_i-{\mathbf{z}}_j)$$

$$=$$

$$\mathbf{K}({\mathbf{z}}_i-{\mathbf{z}}_j)$$

$$=$$

$${\varphi }_{\omega }(\mathbf{z})=\frac{{\sigma }_s}{\sqrt{r}}[\cos ({\omega }^T\mathbf{z})\sin ({\omega }^T\mathbf{z})],\omega \sim \mathcal{N}(0,{\mathbf{W}}^{-1})$$

$${\mathbf{y}}^{\ast }|{\mathbf{z}}^{\ast },\mathbf{Z},\mathbf{Y}$$

$$\varphi (\mathbf{z})$$

$$\mathbf{\Phi }$$

$$\mathbf{A}$$

$$\mathbf{w}$$

$${\mathbf{A}}_{i+1}$$

$${\mathbf{Y}}_{i+1}$$

$$\mathbf{w}=({{\mathbf{R}}^{\mathbf{T}}\mathbf{R}}^{-1})\mathbf{y}$$

$$J=\Psi (\mathbf{x}({\tau }_0),t_0,\mathbf{x}({\tau }_f),t_f)+\frac{t_f-t_0}{2}{\int }_{{\tau }_0}^{{\tau }_f}g(\mathbf{x}(\tau ),\mathbf{u}(\tau ),\tau ;t_0,t_f)d\tau$$

$$\frac{d\mathbf{x}}{dt}=\frac{t_f-t_0}{2}{\mathbf{f}}_t(\mathbf{x}(\tau ),\mathbf{u}(\tau ),\tau ;t_0,t_f)$$

$$\mathbf{C}(\mathbf{x}(\tau ),\mathbf{u}(\tau ),\tau ;t_0,t_f)\le 0$$

$$\tau =\frac{2t}{t_f-t_0}-\frac{t_f+t_0}{t_f-t_0}$$

$$P_{K}(\tau)=\frac{1}{2^{K}K!}\frac{d^{K}}{d\tau^{K}}\big{[}(\tau^{2}-1)^{K}\big{]}$$

$$\mathbf{x}(\tau )\approx \mathbf{X}(\tau )=\sum _{i=1}^K{\mathcal{M}}_i(\tau )\mathbf{X}({\tau }_i)$$

$${\mathcal{M}}_i(\tau )=\prod _{j=1,j\ne i}^K\frac{\tau -{\tau }_j}{{\tau }_i-{\tau }_j}=\frac{g(\tau )}{(\tau -{\tau }_i)\dot{g}(\tau )}$$

$$g(\tau )=(1+\tau )P_k(\tau )$$

$${\int }_a^{b}f(\tau )d\tau \approx \sum _{i=1}^K{w}_{i}f({\tau }_i)$$

$$w_i={\int }_{-1}^1{\mathcal{M}}_i(\tau )d\tau =\frac{2}{(1-{\tau }_i^2)[{\dot{P}}_k({\tau }_i){]}^2},(i=1,...,K)$$

$$\dot{\mathbf{x}}({\tau }_k)\approx \dot{\mathbf{X}}({\tau }_k)=\sum _{i=1}^K\dot{\mathcal{M}}({\tau }_k)\mathbf{X}({\tau }_i)=\sum _{i=1}^K{D}_{ki}\mathbf{X}({\tau }_k),(k=1,...,K)$$

$$D_{ki}=\begin{cases}\frac{(1+\tau_{k})\dot{P}_{K}(\tau_{k})+P_{K}(\tau_{k})}{(\tau_{k}-\tau_{i})\big{[}(1+\tau_{i})\dot{P}_{K}(\tau_{i})+P_{K}(\tau_{i})\big{]}}&\text{if }k\neq i\\ \\ \frac{(1+\tau_{i})\ddot{P}_{K}(\tau_{i})+2\dot{P}_{K}(\tau_{i})}{2\big{[}(1+\tau_{i})\dot{P}_{K}(\tau_{i})+P_{K}(\tau_{i})\big{]}}&\text{if }k=i\par\end{cases}$$

$${\mathbf{R}}_k=\sum _{i=1}^K{D}_{ki}\mathbf{X}({\tau }_k)-\frac{t_f-t_0}{2}\mathbf{f}(\mathbf{X}({\tau }_k),\mathbf{U}({\tau }_k),{\tau }_k;t_0,t_f)=0,(k=1,...,K)$$

$$J=\Psi ({\mathbf{X}}_0,t0,{\mathbf{X}}_f,t_f)+\frac{t_f-t_0}{2}\sum _{k=1}^K{w}_{k}g({\mathbf{X}}_k,{\mathbf{U}}_k,{\tau }_k;t_0,t_f)$$

Full text

Nomenclature

[TABLE]

1 Introduction

Guidance and navigation in aerospace has the capability to leverage high fidelity physics based models when generating a trajectory or designing a controller. These models can be expensive to evaluate or sometimes unavailable; thus an increasing trend in robotics and machine learning is to use a data-driven approach. In these frameworks, the measurements taken during system identification and control tasks are used to build a dynamical model for the system. The drawbacks to purely data-driven approaches are that they often require multiple, repeated trials on the real system which is not always practical. Thus, we hope to leverage the advantages of both techniques through the use of a semi-parametric representation of the dynamics: a combination of a physics based model and a Gaussian Process.

Gaussian processes are a data efficient, non-parametric regression technique that can be utilized to learn any nonlinear function, including dynamics [1]. The method is non-parametric because there is no explicit, parameter-based characterization of the function, just a kernel function between training data and new inputs. In addition to data-efficiency, a Gaussian Process is a probabilistic inference method, meaning that we can determine both the mean and variance of our prediction based on a test input. This has lead to a variety of applications for a Gaussian process. One example is learning the model behind human motion capture data, where each pose is up to 50 dimensions [2]. An extension to this method is to use a semi-parametric approach to learn the inverse dynamics for a 7-DoF BARRET WAM arm [3]. Here they characterize the error dynamics between their model and rollouts to the system with a Gaussian Process. In the context of control, Gaussian processes have been used to characterize the uncertainty of a dynamical model for a learning based robust control algorithm [4]. Finally, in the context of trajectory optimization, there has been work to combine a Gaussian process an algorithm known as Differential Dynamic Programming. As a result, the new algorithm, coined Probabilistic Differential Dynamic Programming, is capable of completely learning the dynamics of a system, while simultaneously performing trajectory optimization on it[5].

Gauss pseudospectral optimal control is our trajectory optimization algorithm of choice because of three reasons: natural incorporation of state and control constraints, integration accuracy of the Gauss quadrature, and efficiency of dynamics function calls. We utilize GPOPS-II for an implementation of the Gauss Pseudospectral optimal control method [6], as well as a in-house python implementation that naturally allows for data driven dynamics. We use a model predictive approach to incorporate feedback into our algorithm and improve robustness against disturbances [7]. MPC also sets up a framework to incorporate online learning. The rest of the paper is organized as follows. First we introduce the mathematical background with the problem formulation, and presentation of Gaussian processes and Gauss Pseudospectral Optimal Control. Next we summarize the base algorithm and present the simulations on the Cart Pole and Quadrotor systems. Finally we introduce an online learning algorithm and compare perform the trajectory optimization tasks for even more uncertain systems, as well as for an obstacle avoidance problem.

2 Mathematical Background

First we present the problem and necessary background material on Gaussian Processes and Gauss Pseudospectral optimal control. The Gaussian process derivation can be found in detail by the Gaussian Process text by Rasmussen [8], and Gauss Pseudospectral has its roots in Hunting’s thesis [16] and has been further extended by Patterson and Rao with GPOPS-II [6].

2.1 Problem Formulation

We would like to minimize a user-defined cost function $J$ subject to dynamics $f_{t}(\mathbf{x},\mathbf{u},t)$ and inequality constraints $h_{s}(\mathbf{x})$ and $h_{c}(\mathbf{u})$.

[TABLE]

$\Psi(\mathbf{x}(t_{f}),t_{f})$ is the terminal cost and $\mathcal{L}(\mathbf{x}(t),\mathbf{u}(t),t)$ is the running cost.

$f_{t}$ represents the true dynamics, a function that we do not have complete knowledge of. Instead, we assume that $f_{t}$ can be represented as the sum of two functions that we do know: a parametric physics model $f_{p}$ and a non-parametric Gaussian Process $f_{err}$. In our current formulation we will only utilize the predictive mean of the Gaussian process as part of our dynamics, there are other frameworks that incorporate a predictive variance [5],[9] .

[TABLE]

Next we go into further detail about Gaussian Processes and how they fit into our framework.

2.2 Gaussian Processes

A Gaussian Process (GP) is defined as “a collection of random variables, any finite number of which have a Gaussian distribution” [8] . We can collect $N$ pairs of state-control pairs, $\mathbf{Z}=\{(\mathbf{x}_{1},\mathbf{u}_{1}),...,(\mathbf{x}_{N},\mathbf{u}_{N})\}$, and the corresponding outputs from our true dynamics $\mathbf{Y}=\{(\mathbf{y}_{1}),...,(\mathbf{y}_{N})\}$ ,when $\mathbf{y}_{i}=f_{t}(\mathbf{x}_{i}(t),\mathbf{u}_{i}(t))-f_{p}(\mathbf{x}_{i}(t),\mathbf{u}_{i}(t))~{}~{}i=1,..,N$. Classical Gaussian processes are scalar-valued function approximators, meaning that we will have an independent Gaussian process for each dimension of $\mathbf{y}_{i}$. Methods exist for correlated vector outputs of Gaussian Processes [10] , however for our purposes we have not found significant evidence for including correlated vector outputs. The output of a Gaussian process is completely defined as a normal distribution with some predictive mean and variance. The joint distribution of the observed output (from our test cases) and the output corresponding to a test case $\mathbf{z^{*}}=(\mathbf{x^{*}},\mathbf{u^{*}})$ can be written as the following equation:

[TABLE]

Note that the joint distribution is characterized with a user defined Kernel function $\mathbf{K}(\mathbf{z}_{i},\mathbf{z}_{j})$. We consider the kernel function corresponding to the squared exponential or Gaussian kernel:

[TABLE]

where $\sigma_{s},\sigma_{n}$, and $\mathbf{W}$ are hyper-parameters that are must be trained with the data. $\sigma_{s}$ is a measure of the variance of the data itself, $\sigma_{n}$ is a measure of noise of the observations, and $\mathbf{W}$ is a diagonal matrix containing the characteristic length scales of each input. We see that the covariance function between outputs is written in terms of the inputs, implying that we are making some assumptions about the smoothness of the underlying function we are trying to learn. The text by Rasmussen [8] shows that this particular covariance function corresponds to a Bayesian linear regression model with an infinite number of basis functions.

Next we can condition the joint prior Gaussian distribution on the observations to obtain the following:

[TABLE]

This is the one step prediction of our output $\mathbf{y}^{*}$ given our test input $\mathbf{z}^{*}$ and our training data $(\mathbf{Z},\mathbf{Y})$. In some applications this prediction can be extended to a multi-step prediction, where we propagate the probabilistic dynamics, however, in our framework, we do not need to do this since the Gauss pseudospectral method approximates the dynamics at each collocation point separately, meaning we only perform a one step prediction. Note that performing a prediction requires inversion of the kernel function $\mathbf{K}(\mathbf{Z},\mathbf{Z})+\sigma_{n}^{2}\mathbf{I}$. This inversion has computational complexity of $\mathcal{O}(N^{3})$ where $N$ is the number of data points. Clearly using a Gaussian process in an online learning fashion is not scalable with respect to the number of data points, as you gather more data you must perform a larger and more difficult inversion at each update. There are techniques such as the inclusion of forgetting factors and local Gaussian Processes to improve the computational performance [11] . Additionally we can perform online learning via an approximate Gaussian Process [12],[13] . This is known as a sparse approximation to a Gaussian kernel or a Sparse Spectrum Gaussian Process.

Sparse Spectrum Gaussian Processes

The key aspect of the Sparse Spectrum Gaussian Process (SSGP) is employing a random Fourier feature approximation of the kernel function [14]. Based on Bochner’s theorem, shift-invariant kernel functions can be represented as the Fourier transform of a unique measure [15].

[TABLE]

where $\mathbf{\phi}(\mathbf{z})=[\cos(\mathbf{\omega}^{T}\mathbf{z})\sin(\mathbf{\omega}^{T}\mathbf{z})]$. The term $p(\mathbf{\omega})$ is a probability distribution. The number of samples $r$ taken from the distribution $p(\mathbf{\omega})$ allows us to construct an unbaised approximate squared exponential function.

[TABLE]

where

[TABLE]

Now that we have this efficient feature mapping, we can represent our underlying dynamics function $f_{err}(\mathbf{z})$ as a weighted sum of the feature functions. We assume a prior distribution of $\mathbf{w}\sim\mathcal{N}(0,\Sigma_{p})$ on the weight matrix. The complete posterior distribution can be calculated as follows:

[TABLE]

where

[TABLE]

With both standard Gaussian Processes and Sparse Spectrum Gaussian Processes, the derivative with respect to the input can be computed analytically, thus enabling a significant cost savings when utilizing trajectory optimization techniques. The whole purpose behind introducing SSGP’s was to incorporate online learning, which we will introduce in the next section.

Online model learning

Incremental online learning provides a way to update the SSGP. The computational complexity of has become $\mathcal{O}(Nr^{2}+r^{3})$ where $r$ is the number of random features of the SSGP. To update the weights $\mathbf{w}$ as new samples come in, we store the upper triangular Cholesky factor $\mathbf{A}=\mathbf{R^{T}R}$ [13]. With each new sample, $\mathbf{R}$ is updated with computational complexity $\mathcal{O}(r^{2})$. We also utilize a forgetting factor $\lambda\in(0,1)$ to assign a lower weight to data that is no longer relevant (particularly in the case of time varying systems).

[TABLE]

The updates to $\mathbf{R}$ are rank-1, and the weights $\mathbf{w}$ can be computed by

[TABLE]

2.3 Gauss Pseudospectral Optimal Control

We present an overview of pseudospectral methods based off our previous work [7]. Pseudospectral optimal control is in a class of algorithms known as ”direct” methods. The general outline [16] is to do the following: first, we will convert the dynamic system into a problem with a finite set of variables and transform the differential constraints into algebraic constraints with spectral approximation for the derivatives. This will allow us to solve with a non-linear program (NLP) solver. Second, we will solve this finite dimensional problem with the NLP and optimize the parameters with which we approximated the original system. Third, we can assess the accuracy of our approximation utilizing the costate and the Karush-Kuhn-Tucker (KKT) conditions. The advantage of spectral differentiation for smooth problems is that it exhibits a fast convergence rate. The Gauss Pseudospectral Method is described in much further detail in the referenced thesis by Huntington [16].

Let us begin with a few equations to describe the problem. The Pseudospectral methods are defined with the Transformed Continuous Bolza problem. The problem formulation given by Equations (LABEL:ProblemCost) and (3) must be put into the following form:

[TABLE]

This cost will be minimized subject to the following constraints.

[TABLE]

[TABLE]

[TABLE]

The function $g$ represents the integrated cost similar to the purpose of $\mathcal{L}$ in DDP. Equation (22) is the dynamic constraint to the problem. Equation (23) is the boundary condition at the final state, similar to our constraints $h_{c}$ and $h_{p}$ in the original problem. Equation (24) shows the time transformation required for this derivation of pseudospectral optimal control. Typically, we require a fixed time interval such as $[-1,1]$, thus such a transformation can be used that is still valid with free initial and final times $t_{0}$ and $t_{f}$.

Collocation Point Search

Collocation points are where we set the approximation to the nonlinear dynamic constraints equal to the true function across the interval. Below is the equation for the Legendre-Gauss points. The points are defined as the roots of Equation (25). Equation (25) is the $K^{th}$ degree Legendre polynomial.

[TABLE]

This choice of points mitigates the Runge phenomenon and improves accuracy towards the boundaries.

State, Control, and Costate Approximation

Next we will approximate the state, control, and the costate. Note that we do not have to use the same collocation points and basis functions for each state, control, and costate, but we do use the same process.

[TABLE]

[TABLE]

[TABLE]

Equation (26) represents the state approximation in terms of a set of Lagrange interpolating polynomials. Equation (27) represents the actual polynomials, where the $P_{K}$ is the Legendre Polynomial.

Integral Approximation via Gauss Quadrature

Within the dynamic constraints and cost function we also have integrals that must be approximated. One way to do this is via a Gauss quadrature. The generic form for this is in Equation (29).

[TABLE]

The points $\tau_{k}$ are the quadrature points on the interval $[-1,1]$ and the weights $w_{i}$ are the quadrature weights. Equation (30) is for finding these weights.

[TABLE]

Algrebriac Representation of Dynamics

The next step is to convert the differential equation constraint on the optimization problem into an algebraic constraint. This is done in Equation(31)

[TABLE]

The differentiation matrix which is given in Equation (32).

[TABLE]

We can now implement this algebraic representation in terms of a residual function by equating the derivative of our approximation to the vector field.

[TABLE]

Gauss Pseudospectral Discretization of the Continuous Bolza Problem

After these transformations, we collect terms for our new problem. We want to minimize the following cost function Equation (34) subject to the algebraic collocation constraints in Equation (35). We also have the quadrature constraint in Equation (36), and path constraint in Equation (37).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The procedure to solve is as follows.

Choose the initial and final times, then choose your basis functions and use the Legendre polynomials to get the collocation points. 2. 2.

Translate the dynamic system into the the discretized and approximated version of the Continuous Bolza Problem. Also translate the constraints into the equivalent forms. 3. 3.

Solve the finite-dimensional problem with a Non-linear Program Solver. 4. 4.

Check the accuracy of the problem by the costate approximation, and then run again with increased collocation points, ideally moving towards the optimum.

One program utilized to implement this method is called GPOPS II [6]. Note we directly substitute $f_{p}+f_{err}=f_{t}$ into the optimizer. We also implement this technique in python, the main difference in our programming being that we allow the representation of our error dynamics $f_{e}rr$ to be an arbitrary function with arbitrary analytical derivatives that are based on data. This way machine learning techniques can be more efficiently implemented with this powerful framework. Also of note is that we can naively incorporate the variance of our Gaussian Process prediction in the instance of obstacle avoidance. Since Pseudospectral Optimal Control allows for continuous nonlinear constraints that are satisfied explicitly in the problem formulation, we can write these as obstacles and perform path planning and control simultaneously with a probabilistic bound.

3 Results

For the following simulation results, we implement GP’s with GPOPS-II on two primary systems: a cart pole and quadrotor model. On the cart pole we introduce unmodeled viscous friction and learn the error dynamics through an untuned feedback controller. For the quadrotor model, we have a $10\%$ mass difference between the true model and the parametric model, and we hope to absorb the resulting error into the Gaussian process.

3.1 Cart Pole

The first system is the classic cart pole. The equations of motion are given in Equation (42), where the task is to swing up the cart pole from a downwards stable equilibrium to an unstable upwards equilibrium. Here $x$ is the position of the cart, $\theta$ is the angle of the pendulum with respect to the cart, and $f$ is the control force pushing the cart from side to side. The following equation is the true dynamics.

[TABLE]

We trained our Gaussian process where the input was the entire state vector and the control, and the output was $\ddot{x},\ddot{\theta}$. We generated 200 points of data using a sinusoidal desired trajectory from a crude feedback controller. However, even with this relatively small amount of data we were able to swing up, even though there was error in the desired position.

We see that even minor unmodeled dynamics can have a drastic effect on the performance a nonlinear system in the episodic case. However, using this regression technique we were able to compensate. In the receding horizon case we see that performance is much improved, however the GP case is marginally more accurate. Neither achieve exact convergence, illustrating a need for a more advanced algorithm.

3.2 Quadrotor

Now we explore the quadrotor, whose equations of motion are provided in Equation (55). The dynamical system has 12 states corresponding to position velocity, orientation, and angular rates. We use a simplified model of the quadrotor that assumes small deviations in the attitude. The position in 3D space is defined by $x$, $y$, and $z$, while the symbols $\phi$, $\theta$, and $\psi$ represent the roll, pitch, and yaw angles, respectively. There are now four controls in this system as opposed to singular control in the other systems, these four controls correspond to the thrust of each rotor. The initial position was $[-1,1,.5]$ and the target state was $[.5,-1,1.5]$. The cost function now has an additional term which penalizes excessive deflection of the pose, preventing the quadrotor system from rolling, pitching, or yawing large amounts to retain the validity of our model. In this case we did not introduce any new terms, we perturbed the mass of the quadrotor. The quad is defined WRT a north east down frame so a positive z indicates a downward movement.

[TABLE]

In this case, we see that there is a primary deviation in the state from the target in our z position. The Gaussian process is trained from sinusoidal data in the x and y direction (i.e flying in circles) but had remained level in the z direction. We used 200 data points, but clearly the GP does not have enough data to make accurate predictions of the behavior in the z direction. The marginal difference between the nominal and semi-parametric dynamics further reaffirms the fact that our Gaussian Process is not contributing significantly due to a lack of data.

In the MPC case, we again see that the z position of both systems are off the target, and the GP does not have enough data to compensate. As a result, another approach is required, where we can incorporate data as we perform rollouts of the trajectory optimizer. This is the reason we introduced Sparse Spectrum Gaussian Processes and Online Model Learning.

3.3 Online Learning with Sparse Spectrum Gaussian Processes

In this framework, we aim to update our Gaussian Process as we receive new data in order to improve performance. Note that since a Full Gaussian Process is data-driven, it retains all of the training data whenever a prediction is made. Thus adding more data online is not computationally efficient, as we must invert an $N$ dimensional kernel matrix each time we update. Sparse Spectrum GP’s circumvent this problem by approximating this kernel matrix with randomly sampled Fourier features. Each data point that is an update to the SSGP constitutes a rank-1 matrix which makes the update computationally efficient. Note that we must optimize the hyperparameters of the system with an initial set of data. This means that we must have a representative set of the data that we will see, otherwise our predictions will be biased, even with the online learning. Next we present the algorithm for Model Predictive Pseudospectral Control with Online Learning:

This simple addition to Model Predictive Pseudospectral Control will incorporate the new data as we roll it out. Currently, the issue with the implementation is that the SSGP must be defined and updated in a vectorized fashion in order for the computational to be practical. Since the number of calls to the SSGP is equal to the number of meshpoints, as the grid becomes more fine, we must perform inference many times. If the inference is done sequentially, then it is unfeasible to run this algorithm.

Cart Pole, and Quadrotor Simulations

We modify the parameters further for the SSGP test to show the increased robustness of this algorithm. We increase the unmodeled damping for the cart pole system and add a one percent perturbation in the length of the quadrotor arm.

The following simulations were run with an in-house implementation of Pseudospectral Optimal Control (implemented in python) that implements the SSGP framework as a natural extension. The jacobian and hessian calculations were incompatible with ADIGATOR (the GPOPS-II autodifferentiation software). The following table has the nominal and perturbed parameters for both the cart pole and the quadrotor case. In this, the quadrotor has an additional arm length perturbation of one percent, and the cart pole has a much larger amount of unmodeled viscous friction. Again we see that without any compensation of for the dynamics difference, we the cart pole cannot swing up at all. However, using the SSGP online learning algorithm, the cart pole is able to swing up close to the equilibrium point, then actively stabilize itself.

For the quadrotor case the difference becomes even more drastic, as we see exact covergence with the SSGP algorithm fairly quickly. The vanilla case had an error in the height dimension which is clearly not present here.

Obstacle Avoidance using SSGP

In our formulation of a semi-parametric representation, we use a probabilistic inference technique known as a Gaussian Process in order to represent unknown dynamics. A benefit of using this probabilistic machine learning tool is that we can get an estimate of the uncertainty of our prediction. In this section we will show the equations to combine this uncertainty estimate along with the Model Predictive Pseudospectral Control (MPPC) to perform obstacle avoidance using hard state constraints. We can formulate the constraint with the following equation:

[TABLE]

This inequality constraint can be naturally incorporated into the pseudospectral framework, making obstacle avoidance an ideal application. The uncertainty radius is determined by the predictive covariance of the SSGP. We propagate the covariance forward using the following equation, then use this value as a one step prediction on the uncertainty of the system.

[TABLE]

$\Sigma_{x_{1}}$ defines an ellipse around the obstacle body, using the angle between the body and the object, we can find the uncertainty radius $U_{r}$. We applied this on a dubin vehicle that uses bicycle dynamics, however for a deterministic system the performance of the with the variance was identical to performance without using the variance. We apply the technique on the following dynamics system excluding the $U_{r}$ term:

[TABLE]

The true values for $b_{1}$ and $b_{2}$ are 0.3, however the nominal vales are a Gaussian distribution $\mathcal{N}(0.3,0.3)$. Thus we can apply our algorithm over 50 trajectories and look at the distribution over states and over the path. We see that this constraint is adhered to at all points in the trajectory. The problem is to move from the point $(5,5)$ to the point $(2,-2)$, while avoiding an obstacle of radius $2$.

4 Conclusion

In conclusion, we see that even a small amount of viscous friction in the Cart Pole is enough to prevent an episodic trial of PSOC from reaching its target. The Gaussian process model is able to swing up the cart pole, however there is still error in the final position. In the MPC trials, both systems are able to swing up the cart pole and we see a marginal improvement using the Gaussian process. The reason why the GP model does not completely solve the issue is because it is data dependent, and if the Gaussian process input is far away from the training data, then the machine reverts back to the prior, in this case zero. These results are echoed in the quadrotor model. The Gaussian process is trained with sinusoidal tracking data from a MPC controller, however this data does not provide enough information about the z direction dynamics of the true system. It is unlikely that we can perfectly train a learned model a priori before we attempt a task. As a result, it is clear that we must go a step further and incorporate the new knowledge we gain from each rollout. A naive way of incorporating data would be to use the Full Gaussian Process at the expense of heavy computational load. Instead we opt for a Sparse Spectrum Gaussian process and present an algorithm to perform trajectory optimization with online learning. In these examples we see an improvement in the task completion for both the Cart Pole and Quadrotor systems, and we go a step further to perform obstacle avoidance in the Dubin Vehicle. The key points that must be stressed here is that performance of the algorithm hinges on the optimization of the Sparse Spectrum Gaussian Process. If there is not enough data, or if the hyperparameters of the data are not correct, then it is possible to have worse performance than in the vanilla case. Overall this represents a powerful tool in the area of trajectory optimization, and as further work is completed, this algorithm can be implemented on real systems that can simultaneously perform path planning, system identification, and control.

Acknowledgments

This work was was supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program.