Model-Based Policy Search for Automatic Tuning of Multivariate PID Controllers

TL;DR

This paper introduces a model-based policy search method that automatically tunes multivariate PID controllers using data, enabling efficient and effective control of complex systems like robotic arms balancing an inverted pendulum.

Contribution

It extends the PILCO framework to optimize multivariate PID controllers directly from data without prior knowledge, simplifying the tuning process for complex, coupled controllers.

Findings

Demonstrated fast, data-efficient policy learning on a robotic arm balancing an inverted pendulum.

Showed the method's effectiveness in tuning multivariate PID controllers automatically.

Validated the approach on a complex real-world control task.

Abstract

PID control architectures are widely used in industrial applications. Despite their low number of open parameters, tuning multiple, coupled PID controllers can become tedious in practice. In this paper, we extend PILCO, a model-based policy search framework, to automatically tune multivariate PID controllers purely based on data observed on an otherwise unknown system. The system's state is extended appropriately to frame the PID policy as a static state feedback policy. This renders PID tuning possible as the solution of a finite horizon optimal control problem without further a priori knowledge. The framework is applied to the task of balancing an inverted pendulum on a seven degree-of-freedom robotic arm, thereby demonstrating its capabilities of fast and data-efficient policy learning, even on complex real world problems.

Equations

$${\mathbf{x}}_{t+1}$$

$$J=\sum _{t=0}^T\mathbb{E}[c({\mathbf{x}}_t,{\mathbf{u}}_t;t)],$$

$$u_t$$

$$e_t$$

$$C(s)=\left[\begin{array}{ccc}{c}_{11}(s)& \cdots & c_{1D}(s)\\ ⋮& \ddots & ⋮\\ c_{F1}(s)& \cdots & c_{FD}(s)\end{array}\right],$$

$${\mathbf{z}}_t:=({\mathbf{x}}_t,{\mathbf{e}}_{t-1},\Delta T\sum _{\tau =0}^{t-1}{\mathbf{e}}_{\tau }),$$

$$\left[\begin{array}{c}{\mathbf{z}}_t\\ {\mathbf{x}}_{\text{des,t}}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_z\\ {\mu }_{\text{des},t}\end{array}\right],\left[\begin{array}{cc}{\Sigma }_z& 0\\ 0& {\Sigma }_{\text{des},t}\end{array}\right]\right).$$

$${\dot{\mathbf{e}}}_t$$

$${\int }_0^{t\cdot \Delta T}\mathbf{e}(\tau )d\tau$$

$${\tilde{\mathbf{z}}}_t:=({\mathbf{z}}_t,{\mathbf{x}}_{\text{des},t},{\mathbf{e}}_t,\Delta T\sum _{\tau =0}^t{\mathbf{e}}_{\tau },\frac{{\mathbf{e}}_t-{\mathbf{e}}_{t-1}}{\Delta T}).$$

$$\begin{array}{rl}{\mathbf{u}}_t& ={\mathbf{A}}_{\text{PID}}({\tilde{\mathbf{z}}}_t^{(3)},{\tilde{\mathbf{z}}}_t^{(4)},{\tilde{\mathbf{z}}}_t^{(5)})\\ & ={\mathbf{A}}_{\text{PID}}\left({\mathbf{e}}_t,\Delta T\sum _{\tau =0}^t{\mathbf{e}}_t,\frac{{\mathbf{e}}_t-{\mathbf{e}}_{t-1}}{\Delta T}\right),\end{array}$$

$${\mathbf{A}}_{\text{left}}$$

$${\mathbf{A}}_{\text{right}}$$

$${\mathbf{z}}_{t+1}=({\mathbf{x}}_{t+1},{\tilde{\mathbf{z}}}_t^{(3)},{\tilde{\mathbf{z}}}_t^{(4)})=({\mathbf{x}}_{t+1},{\mathbf{e}}_t,\Delta T\sum _{\tau =0}^t{\mathbf{e}}_{\tau }).$$

$${\mathbf{z}}_0:=({\mathbf{x}}_0,{\mathbf{x}}_{\text{des},0}-{\mathbf{x}}_0,0).$$

$$\frac{dJ(\theta )}{d\theta }=\sum _{t=1}^H\frac{d}{d\theta }\underset{=:{\mathcal{E}}_t}{\underbrace{{\mathbb{E}}_{{\mathbf{z}}_t}[c({\mathbf{z}}_t)]}}=\sum _{t=1}^T\frac{d{\mathcal{E}}_t}{dp({\mathbf{z}}_t)}\frac{dp({\mathbf{z}}_t)}{d\theta }.$$

$$\displaystyle\frac{dp(\mathbf{z}_{t+1})}{d\theta}={\color[rgb]{0,0.4470,0.7410}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.4470,0.7410}\frac{\delta p(\mathbf{z}_{t+1})}{\delta p(\tilde{\mathbf{z}}_{t})}}\frac{dp(\tilde{\mathbf{z}}_{t})}{d\theta}+{\color[rgb]{0,0.4470,0.7410}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.4470,0.7410}\frac{\delta p(\mathbf{z}_{t+1})}{\delta p(\mathbf{x}_{t+1})}}\frac{dp(\mathbf{x}_{t+1})}{d\theta}\,.$$

$$\displaystyle\frac{dp(\mathbf{x}_{t+1})}{d\theta}={\color[rgb]{0.8500,0.3250,0.0980}\definecolor[named]{pgfstrokecolor}{rgb}{0.8500,0.3250,0.0980}\frac{\delta p(\mathbf{x}_{t+1})}{\delta p(\tilde{\mathbf{z}}_{t})}}\frac{dp(\tilde{\mathbf{z}}_{t})}{d\theta}+{\color[rgb]{0.8500,0.3250,0.0980}\definecolor[named]{pgfstrokecolor}{rgb}{0.8500,0.3250,0.0980}\frac{\delta p(\mathbf{x}_{t+1})}{\delta p(\mathbf{u}_{t})}}\frac{dp(\mathbf{u}_{t})}{d\theta}\,.$$

$$\frac{dp(\mathbf{u}_{t})}{d\theta}={\color[rgb]{0,0.4470,0.7410}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.4470,0.7410}\frac{\delta p(\mathbf{u}_{t})}{\delta p(\tilde{\mathbf{z}}_{t})}}\frac{dp(\tilde{\mathbf{z}}_{t})}{d\theta}+{\color[rgb]{0,0.4470,0.7410}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.4470,0.7410}\frac{\delta p(\mathbf{u}_{t})}{\delta\theta}}\,,$$

$$\displaystyle\frac{dp(\tilde{\mathbf{z}}_{t})}{d\theta}={\color[rgb]{0,0.4470,0.7410}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.4470,0.7410}\frac{dp(\tilde{\mathbf{z}}_{t})}{dp(\mathbf{z}_{t})}}\frac{dp(\mathbf{z}_{t})}{d\theta}\,$$

$$c({\mathbf{e}}_t,{\bar{\mathbf{u}}}_t)=1-\exp (-({\mathbf{e}}_t^{T}Q{\mathbf{e}}_t+{\bar{\mathbf{u}}}_t^{T}R{\bar{\mathbf{u}}}_t)/2).$$

$${\mathbf{x}}_{t+1}=\hat{f}({\mathbf{x}}_t,\dots ,{\mathbf{x}}_{t-n},{\mathbf{u}}_t,\dots ,{\mathbf{u}}_{t-m}).$$

$$\mathbf{Y}=\mathbf{AX}+\mathbf{b}=\mathcal{N}({\mu }_Y,{\Sigma }_Y)\in {\mathbb{R}}^P,$$

$$\left[\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}{\mu }_X\\ {\mu }_Y\end{array}\right],\left[\begin{array}{cc}{\mathbf{\Sigma }}_X& {\mathbf{\Sigma }}_X\mathbf{C}\\ {\mathbf{C}}^T{\mathbf{\Sigma }}_X^T& {\mathbf{\Sigma }}_Y\end{array}\right])$$

$${\mu }_Y=\mathbf{A}{\mu }_X+\mathbf{b},{\Sigma }_Y=\mathbf{A}{\Sigma }_X{\mathbf{A}}^T,\mathbf{C}={\mathbf{A}}^T$$

$$\frac{\delta {\mu }_Y}{\delta {\mu }_X}$$

$$\frac{\delta {\mu }_Y}{\delta A}$$

$$\frac{\delta C_{kl}}{\delta A_{ij}}$$

$$\frac{\delta ({\Sigma }_Y{)}_{kl}}{\delta A_{ij}}$$

Full text

I Introduction

Proportional, Integral and Derivative (PID) control structures are still the main control tool being used in industrial applications, in particular in the process industry [1], but also in automotive applications [2] and in low-level control in robotics [3]. The large share of PID controlled applications is mainly due to the past record of success, the wide availability, and the simplicity in use of this technique.

In practice, control design is still often achieved by tedious manual tuning or by heuristic PID tuning rules [4]. More advanced tuning concepts are most frequently developed for Single-Input-Single-Output (SISO) systems [5, 6]. For Multi-Input-Multi-Output (MIMO) systems, popular tuning methods, such as biggest log-modulus and the dominant pole placement tuning method [7], strive to tune each control loop individually, followed by a collective de-tuning to stabilize the multi-loop system. These tuning methods, however, rely on linear process models and require stable processes. For general PID control structures where multiple controllers act on each input, controller design is usually conducted by decoupling the process, subsequently allowing the design of individual SISO PIDs. One example are online adjusted precompensators, which decouple the process transfer function matrix [8].

In this paper, we extend Probabilistic Inference for Learning COntrol (PILCO) [9], a framework for iteratively improving a controller based on its expected finite horizon cost as predicted from a learned system model. This model is learned in a fully Bayesian setting using Gaussian Process Regression (GPR) as a non-parametric function estimator [10]. Gaussian Processes (GPs) as probability distributions over a function space do not only provide a prediction but also an uncertainty measure of the process dynamics. The uncertainty component can be utilized to implement cautious control [11] or trigger further exploration to improve the model knowledge [12]. As it is based on the full nonlinear system model, this framework takes into account all process couplings for controller tuning. In contrast to some of the aforementioned PID tuning concepts, e.g. SISO or multi-loop tuning methods, this framework imposes no structural restrictions on the multivariate PID control design.

PILCO has been successfully applied to reinforcement learning tasks, such as inverted pendulum swing-up [9], control of low-cost manipulators [13], and real-world applications like throttle valve control [14]. These examples focus mostly on nonlinear state feedback policies. In industrial applications, however, it is desirable to obtain interpretable control designs, which is the case for PID control structures rather than for arbitrarily complex nonlinear control designs.

Contribution of the paper: We propose a general framework for multivariate PID controller tuning based on PILCO. In particular, the developed framework (i) is applicable to nonlinear MIMO systems with arbitrary couplings; (ii) allows for the tuning of arbitrary multivariate PID structures; (iii) involves a probabilistic treatment of recorded system data to iteratively improve control performance without a priori knowledge; and (iv) requires no process model. The auto-tuning method is demonstrated in pole balancing experiments on Apollo, a complex robot platform as shown in Fig. 1, coping with imperfect low-level tracking controllers and unobserved dynamics.

Outline of the paper: The paper continues with an introduction of the policy search problem and PILCO in Sec. II. The proposed multivariate PID tuning framework is then developed in Sec. III. Section IV presents the results of applying the framework for tuning coupled PID controllers on Apollo. The paper concludes with remarks and propositions for future work in Sec. V.

II Probabilistic Inference for Control

In this section, we introduce the policy search problem and notation. Subsequently, the main ideas of PILCO [9, 15] are briefly discussed.

II-A Problem Statement

We consider discrete time dynamic systems of the form

[TABLE]

with continuously valued state $\mathbf{x}_{t}\!\in\!\mathbb{R}^{D}$ and input $\mathbf{u}_{t}\!\in\!\mathbb{R}^{F}$. The system dynamics $f$ is not known a priori. We assume a fully measurable state, which is corrupted by zero-mean independent and identically distributed (i.i.d.) Gaussian noise, i.e. $\epsilon_{t}\!\sim\!\mathcal{N}({\bf 0},\Sigma_{\epsilon})$.

One specific reinforcement learning formulation aims at minimizing the expected cost-to-go given by

[TABLE]

where an immediate, possibly time dependent cost $c(\mathbf{x}_{t},\mathbf{u}_{t};t)$ penalizes undesired system behavior. Policy search methods optimize the expected cost-to-go $J$ by selecting the best out of a range of policies $\mathbf{u}_{t}\!=\!\pi(\mathbf{x}_{t};\theta)$ parametrized by $\theta$. Particularly in model-based policy search frameworks, a model $\hat{f}$ of the system dynamics $f$ is utilized to predict the system behavior and to optimize the policy.

II-B PILCO

PILCO as a specific model-based policy search framework emphasizes data-efficiency and consistent handling of uncertainty when constructing the system dynamics model $\hat{f}$. To incorporate all available data from policy rollouts (experiments) on the actual system, a Gaussian Process (GP) [16] is utilized as a non-parametric, probabilistic model.

The PILCO framework is outlined in Alg. 1. In the inner loop, a simulated rollout is conducted based on the dynamics model $\hat{f}$ and the current policy $\pi(\mathbf{x}_{t};\theta)$. The system’s state $\mathbf{x}_{t}$ is propagated over a finite prediction horizon $H$ starting at the system’s initial state $\mathbf{x}_{0}\!\sim\!\mathcal{N}(\mu_{0},\Sigma_{0})$ as visualized for one time step in Fig. 2 in black for the standard PILCO rollout. The posterior distribution $p(\mathbf{x}_{t+1}|\mathbf{x}_{t},\mathbf{u}_{t})$ is approximated in each time step by a Gaussian distribution using moment matching [17] yielding an approximately Gaussian marginal distribution of the long-term predictions $p(\mathbf{x}_{0}),\ldots,p(\mathbf{x}_{H})$. The expected long-term cost (2) of this rollout as well as its gradients with respect to the policy parameters $\theta$ can be computed analytically. Policy optimization can then be conducted based on this prediction method using standard gradient-based optimization techniques.

Starting with an initial random policy111A random policy might be a parametrized policy with randomly chosen parameters (cf. [18]) or a random input signal exciting the system to generate initial dynamics data., the algorithm optimizes the cost-to-go by repeatedly executing the policy on the system, thereby gathering new data and building the dynamics model, subsequently improving the policy iteratively, until the task has been learned.

III Policy Search for PID Control

We first introduce the general description of the considered PID control structure in Sec. III-A; the necessary modifications to the policy search framework are derived in Sec. III-B to Sec. III-D, followed by an analytic derivation of the required cost function gradients in Sec. III-E.

III-A PID Control Policies

The control output of a scalar PID controller is given by

[TABLE]

where $e(\tau)$ denotes the continuous error signal. The current desired state $x_{\text{des},t}$ can be either a constant set-point or a time-variable goal trajectory. This controller is agnostic to the system dynamics and depends only on the system’s error. Each controller is parametrized by its proportional, integral and derivative gain ($\theta_{\text{PID}}\!=\!(K_{\text{p}},K_{\text{i}},K_{\text{d}})$).

A general PID control structure $C(s)$ for MIMO processes (1) can be described in transfer function notation by a $D\times F$ transfer function matrix

[TABLE]

where $s$ denotes the complex Laplace variable and $c_{ij}(s)$ are of PID type. The multivariate error is given by $\mathbf{e}_{t}\!=\!\mathbf{x}_{\text{des},t}\!-\!\mathbf{x_{t}}\!\in\!\mathbb{R}^{D}$ such that the multivariate input becomes $\mathbf{u}(s)\!=\!C(s)\mathbf{e}(s)$. PID literature typically distinguishes between tuning methods for multi-loop PID control and multivariable PID control systems. The former have a diagonal transfer function matrix $C(s)$ whereas the latter allows PID controllers on all elements of $C(s)$, i.e. all combinations of errors and inputs can be controlled, thus allowing additional cross couplings.

With this framework, we address the general class of multivariable PID control systems with no restrictions on the elements of $C(s)$. Examples for possible PID structures are shown in Fig. 3.

III-B System State Augmentation

We present a sequence of state augmentations such that any multivariable PID controller (5) can be represented as a parametrized static state feedback law. A visualization of the state augmentation integrated into the one-step-ahead prediction is shown in red in Fig. 2 in comparison with the standard PILCO setting (in black). Given a Gaussian distributed initial state $\mathbf{x}_{0}$, the resulting predicted states will remain Gaussian for the presented augmentations.

To obtain the required error states for each controller (3), we define a new system state $\mathbf{z}_{t}$ that keeps track of the error at the previous time step and the accumulated error,

[TABLE]

where $\Delta T$ is the system’s sampling time. For simplicity, we denote vectors as tuples $(\mathbf{v}_{1},\ldots,\mathbf{v}_{n})$, where $\mathbf{v}_{i}$ may be vectors themselves. The following augmentations are made to obtain the necessary policy inputs:

III-B1 Desired Goal State

The desired set-point or target trajectory state222Drawing the desired state from a Gaussian distribution yields better generalization to unseen targets as discussed in [13]. $\mathbf{x}_{\text{des},t}\!\sim\!\mathcal{N}(\mu_{\text{des},t},\Sigma_{\text{des},t})$ is independent of $\mathbf{z}_{t}$ yielding:

[TABLE]

III-B2 Error States

The current error is a linear function of $\mathbf{z}_{t}$ and $\mathbf{x}_{\text{des},t}$. The current error derivative and integrated error are approximated by

[TABLE]

Both approximations are linear transformations of the augmented state. The resulting augmented state distribution remains Gaussian as it is a linear transformation of a Gaussian random variable (see Appendix).

The computation of the derivative error is prone to measurement noise. Yet, this framework can readily be extended to incorporate a low-pass filtered error derivative, which we omit for notational simplicity. In this case, additional historic error states would be added to the state $\mathbf{z}_{t}$ to provide the input for a low-pass Finite Impulse Response (FIR) filter.

The fully augmented state is given by

[TABLE]

III-C PID as Static State Feedback

Based on the augmented state $\tilde{\mathbf{z}}_{t}$, the PID control policy for multivariate controllers can be expressed as a static state feedback policy:

[TABLE]

where $\tilde{\mathbf{z}}_{t}^{(i)}$ indicates the i-th term of (10). The specific structure of the multivariate PID control law is defined by the parameters in $\mathbf{A}_{\text{PID}}$. For example, PID structures as shown in Fig. 3 would be represented by

[TABLE]

III-D State Propagation

Given the Gaussian distributed state and control input as derived in Sec. III-B and Sec. III-C, the next system state is computed using the GP dynamics model $\hat{f}$. PILCO approximates the predictive distribution $p(\mathbf{x}_{t+1})$ by a Gaussian distribution using exact moment matching. From the dynamics model output $\mathbf{x}_{t+1}$ and the current error stored in $\tilde{\mathbf{z}}_{t}$, the next state is obtained as

[TABLE]

Iterating (6) to (14), the long-term prediction can be computed over the prediction horizon $H$ as shown in Fig. 2. For the initial state, we define

[TABLE]

III-E Cost Function Derivatives

Given the presented augmentation and propagation steps, the expected cost gradient can be computed analytically such that the policy can be efficiently optimized using gradient-based methods. We summarize the high-level policy gradient derivation steps to point out the modifications to standard PILCO that are necessary to allow PID policy optimization. The expected cost333We only address cost on the state $\mathbf{z}_{t}$ for simplicity. For practical implementations, cost on $\mathbf{u}$ can be included by adding past inputs into the system state as shown in Sec. IV. derivative is obtained as

[TABLE]

To simplify the notation, we write $dp(\mathbf{z}_{t})$ to denote the sufficient statistics derivatives $d\mu_{t}$ and $d\Sigma_{t}$ of a Gaussian random variable $p(\mathbf{z}_{t})\!\sim\!\mathcal{N}(\mu_{t},\Sigma_{t})$ (analogous to the treatment in [15]). The gradient of the immediate loss with respect to the state distribution, $d\mathcal{E}_{t}/dp(\mathbf{z}_{t})$, is readily available for most standard cost functions like quadratic or saturated exponential terms and Gaussian input distributions (cf. [13]). The gradient for each predicted state in the long-term rollout is obtained by applying the chain rule to (14) resulting in

[TABLE]

The derivatives highlighted in blue are computed for the linear transformation in (14) according to the general rules for linear transformations on Gaussian random variables as summarized in the appendix. Based on the dynamics model prediction as detailed in Sec. III-D, the gradient of the dynamics model output $\mathbf{x}_{t+1}$ is given by

[TABLE]

The derivatives shown in red can be computed analytically for the specific dynamics model [17]. Applying the chain rule for the policy output $p(\mathbf{u}_{t})$ obtained by (11) yields

[TABLE]

The derivatives marked in blue are introduced by the linear control law (11) and can be computed as summarized in the appendix. The gradient of the augmented state is given by

[TABLE]

Again, the part marked in blue is computed for the linear transformation (10). Starting from the initial state where $dp(\mathbf{z}_{0})/d\theta\!=\!0$, we obtain the gradients for all states with respect to the policy parameters $dp(\mathbf{z}_{t})/d\theta$ by iteratively applying (17) to (20) for all time steps $t$.

IV Experimental Evaluation

To demonstrate the capabilities of the presented framework to automatically tune coupled PID controllers without prior system knowledge and in a data-efficient fashion, we consider the problem of balancing an inverted pendulum on the Apollo robot as shown in Fig. 1. The inverted pendulum is a well-known benchmark in the control and reinforcement learning communities [19, 20]. Demonstrations of the iterative learning process and the resulting optimized policy can be found in the supplementary video material.

IV-A Experimental Setup

We employ an imperfect inverse dynamics model of Apollo’s seven Degree-of-Freedom (DoF) robotic arm to compute the joint torques necessary to track the desired end effector acceleration $u_{t}$ [21]. Technical details concerning the hardware platform can be found in [22], where reinforcement learning on this platform has been addressed using Bayesian Optimization techniques. The state of the system $\mathbf{x}_{t}$ comprises the end effector position $x$, velocity $\dot{x}$, pendulum angle $\phi$, and angular velocity $\dot{\phi}$.

During the rollouts, the commanded acceleration is limited to $u_{\text{max}}\!=\!3$ m/s2 for safety reasons. Test policies are executed for 20 s or interrupted once safety limits ($x_{\text{max}}\!=\!0.3$ m, $\theta_{\text{max}}\!=\!30^{\circ}$) are violated. The control signal is computed at 100 Hz and low-pass filtered by a second order Butterworth filter having a cut-off frequency of 20 Hz.

The policy is optimized on a prediction horizon of $T\,=\,10\,s$ based on a saturated loss function [9] given by

[TABLE]

For balancing the pendulum in the upright position the desired trajectory is given by $\mathbf{x}_{\text{des},t}=\mathbf{0}$. Weights are set to $Q=\text{diag}(1/0.2^{2},1/0.02^{2})$ for end effector and pendulum position error, and to $R=1/0.4^{2}$ for the control input. The selected cost function saturates quickly for $x\!>\!w$ if $x$ is weighted by $1/w^{2}$. End effector position and input are therefore only penalized lightly, which permits higher gains while the pendulum angle is stabilized as it is penalized much stronger.

IV-B PID Control Structure Setup

We employ a PID controller on the position error $e_{\text{x}}=x_{\text{des},t}-x$ and a PD controller acting on the pendulum angular error $e_{\theta}=\theta_{\text{des},t}-\theta$. The resulting control structure is shown in the right plot of Fig. 3. The PID/PD control structure with integral control on the end effector position serves to correct for any static bias in the angle measurement (e.g., from imperfect calibration) as is explained in [23, p. 67]. This structure has successfully been used for other balancing problems [22, 23].

The specific structure is chosen based on prior knowledge about the problem at hand as detailed in [22]. The integrator’s contribution is required to counteract any pendulum angle measurement bias introduced by imperfect sensor calibration. The policy parametrization is therefore given by $\theta\!=\!(K_{\text{p,x}},K_{\text{i,x}},K_{\text{d,x}},K_{\text{p},\theta},K_{\text{d},\theta})$. Assuming no prior knowledge, we initialize the policy to zero. Both controllers are coupled by the system dynamics and can therefore not be tuned independently, which makes the inverted pendulum a well suited benchmark for multivariate PID controller tuning.

IV-C Modifications for Dynamics Model Learning

When first training GP dynamics models $\mathbf{x}_{t+1}=\hat{f}(\mathbf{x}_{t},\mathbf{u}_{t})$ in the standard way (Sec. II), this did not lead to acceptable models for long-term predictions and thus successful controller learning. The problems are caused by imperfections in the inverse dynamics model, joint friction and stiction, as well as sensor and actuator delay. These factors add unobserved states and therefore additional dynamics to the system, corrupting the measured data. This is visible in Fig. 4, where the desired acceleration and the numerically computed actual acceleration are visualized for a policy rollout on the robot. Several adaptations to the standard GP dynamics model learning framework were required to obtain a good prediction model, which we explain next.

IV-C1 Gaussian Process Setup

In contrast to the model in (1), we train the GP models to predict the difference between the current and the next state $\Delta\mathbf{x}_{t}=\mathbf{x}_{t+1}-\mathbf{x}_{t}$. We compute a sparse GP using Snelson’s approximation [24] having a covariance parametrized by 400 inducing pseudo-input points. For both GP models, hyperparameters $\theta_{\text{GP},i}=(l_{1},l_{2},\sigma_{\text{f}},\sigma_{\text{n}})$ including lengthscales for each dimension, as well as signal and noise variance have to be chosen. We chose the maximum likelihood estimate (MLE) on the basis of the previously gathered system data to compute the hyperparameters. These hyperparameters are kept constant during the iterative policy learning.

IV-C2 NARX Dynamics Model

Instead of modeling the system’s full, four dimensional state and its dynamics, two independent GPs are trained to model the dynamics of the measured state parts; the end effector and pendulum position. The missing information about the system’s velocities and potential latent states is recovered by employing a Nonlinear AutoRegressive eXogenous model (NARX) [25] of the form

[TABLE]

Information on latent states is implicitly encoded in the measured historic states and inputs. Different lengths of history might be required for individual parts of the system’s state depending on the dynamics’ time scales. We optimize the number of historic states individually for end effector position, pendulum position and control input. The NSGA II optimizer [26] is employed to compute the pareto front of the model’s prediction error (using the same dataset as previously used for hyperparameter tuning) and the size of the NARX history. For our problem, we ended up with a new state of dimensionality 14 given by $\mathbf{x}_{\text{NARX}}:=(x_{t},\ldots,x_{t-3},\phi_{t},\ldots,\phi_{t-2},u_{t},\ldots,u_{t-6})$. The hidden dynamics between commanded input and executed acceleration (cf. Fig. 4) requires a longer history in the input state to capture all relevant effects.

IV-C3 Data Preprocessing

The recorded data is downsampled to 25 Hz and non-causally low-pass filtered with a 2nd order Butterworth filter and cut-off frequency 12.5 Hz. The downsampled control input is obtained by averaging the input signal on each sampling interval.

IV-D PID Learning Results on Apollo

The iterative learning experiment is visualized in Fig. 5. To gather initial data about the system, four random rollouts are conducted, applying a white noise input $u_{\text{rnd},t}\!\sim\!\mathcal{N}(0,1$ m/s${}^{2})$ to the system. The first iteration is based on the dynamics model learned from these random rollouts. The left figure visualizes the optimized predicted loss for each iteration in comparison to the actual loss obtained from evaluating (21) on one sample rollout of the actual system. The drop of the predicted loss in iteration 5 shows that by that time, sufficient data about the system dynamics has been gathered to find a stabilizing policy parametrization. In the right figure, the predicted loss optimization is visualized for each iteration individually, as a function of the number of linesearches conducted by the BFGS optimizer.

In Fig. 6, we visualize the dynamics learning progress, showing predicted system behavior and actual rollout for the dynamics models and policies as obtained at the first, intermediate and final stage of the iterative learning process. The iterative improvement in model prediction accuracy and the improvement in the PID policy is clearly visible. The final dynamics model is accurately predicting the stabilization of the system by the optimized policy.

In this example, the total interaction time with the physical system is only 106 seconds, demonstrating fast and data-efficient learning. This model-based method outperforms a model-free, Bayesian optimization method (cf. [22]), with respect to the number of rollouts on the actual system. The policy optimization itself is carried out offline, and the predicted system behavior can be utilized to set appropriate safety boundaries to test new controllers without damaging the system.

To demonstrate the robustness of the learned policy, the system is manually deflected. Disturbances in pendulum angle (cf. Fig. 7) and end effector position are dispelled fast and without overshoot. By commanding a non-zero desired trajectory, the learned PID controller can be utilized for tracking tasks as demonstrated in the supplementary video for a sinusoidal end-effector trajectory.

V Conclusion

The proposed framework for multivariate PID tuning is flexible with respect to possible PID control structures and process dynamics. In particular, it is able to cope with general nonlinear MIMO processes and multivariate PID structures. The presented framework can readily be extended to cascaded PID structures and tracking controllers by considering multiple different [18] and time-varying goal states.

Appropriately dealing with hidden and low-level dynamics, problems found in almost all real-world applications, were major hurdles in the experimental application. In particular, we found that learning of dynamics models geared towards long-term predictions is key to successful finite horizon policy optimization, but largely unaddressed by current approaches. The presented dynamics model structure and learning framework helps to alleviate this problem. Principled ways for improving the learning of long-term prediction models shall be addressed in future work.

APPENDIX

A linear transformation of a Gaussian random variable $\mathbf{X}\sim\mathcal{N}(\mu_{X},\Sigma_{X})\in\mathbb{R}^{D}$ is given by

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{P\times D}$ and $\mathbf{b}\in\mathbb{R}^{P\times 1}$. The joint probability distribution of $\mathbf{X}$ and $\mathbf{Y}$ is given by

[TABLE]

where

[TABLE]

The non-zero partial derivatives of $Y$’s sufficient statistics are given by

[TABLE]

[TABLE]

where $\otimes$ is the Kronecker product and $\delta_{ij}$ is the Kronecker delta (cf. [27] for useful matrix derivatives).