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

arXiv:1703.02899v1 [cs.LG] 8 Mar 2017

Andreas Doerr1,2 , Duy Nguyen-Tuong1 , Alonso Marco2 , Stefan Schaal2,3 , Sebastian Trimpe2
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.

I. I NTRODUCTION
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 im*This research was supported in part by Robert Bosch GmbH, the
Max Planck Society, National Science Foundation grants IIS-1205249, IIS1017134, EECS-0926052, the Office of Naval Research, and the Okawa
Foundation.
1 Bosch Center for Artificial Intelligence, Renningen, Germany.
{andreas.doerr3, duy.nguyen-tuong}@de.bosch.com
2 Autonomous Motion Department at the Max Planck Institute for Intelligent Systems, Tübingen, Germany.
{amarco, strimpe}@tue.mpg.de, sschaal@is.mpg.de
3 Computational Learning and Motor Control lab at the University of
Southern California, Los Angeles, CA, USA.

Fig. 1. Humanoid upper-body robot, Apollo, balancing an inverted pendulum. Using the proposed framework, coupled PID and PD controllers are
trained to stabilize the pole in the central, upright position without requiring
prior knowledge of the pendulum system dynamics.

proving 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 multiloop 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;

Accepted final version. To appear in 2017 IEEE International Conference on Robotics and Automation.
c 2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of
any copyrighted component of this work in other works.

π(xt ; θ)
x0 , z0 (15)

(11)
x t , zt

(10)

ut

GP fˆ

z̃t

xt+1

(14)

zt+1

iterate for t ∈ [0, H]

Fig. 2. Computations for one time step in the long-term prediction. Black:
original PILCO state propagation steps. Red: Augmented state propagation
to accommodate PID policy optimization.

(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. P ROBABILISTIC I NFERENCE FOR C ONTROL
In this section, we introduce the policy search problem
and notation. Subsequently, the main ideas of PILCO [9],
[15] are briefly discussed.
A. Problem Statement
We consider discrete time dynamic systems of the form
xt+1 = f (xt , ut ) + t ,

(1)

with continuously valued state xt ∈ RD and input ut ∈
RF . The system dynamics f is not known a priori. We
assume a fully measurable state, which is corrupted by zeromean independent and identically distributed (i.i.d.) Gaussian
noise, i.e. t ∼ N (0, Σ ).
One specific reinforcement learning formulation aims at
minimizing the expected cost-to-go given by
J=

T
X

E[c(xt , ut ; t)],

(2)

t=0

where an immediate, possibly time dependent cost
c(xt , ut ; 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 ut = π(xt ; θ)
parametrized by θ. Particularly in model-based policy search
frameworks, a model fˆ of the system dynamics f is utilized
to predict the system behavior and to optimize the policy.

(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 fˆ and the current policy π(xt ; θ). The system’s
state xt is propagated over a finite prediction horizon H
starting at the system’s initial state x0 ∼ N (µ0 , Σ0 ) as
visualized for one time step in Fig. 2 in black for the standard
PILCO rollout. The posterior distribution p(xt+1 |xt , ut ) 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(x0 ), . . . , p(xH ). The expected long-term cost (2) of this
rollout as well as its gradients with respect to the policy parameters θ 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 policy1 , 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.
Algorithm 1 PILCO
1: Experiment: Execute random policy;
2:
record {xt , ut }t=1...T
3: Train initial GP dynamics model xt+1 = fˆ(xt , ut )
4: repeat
5:
repeat
6:
Simulation: Predict J(θ) given fˆ, π(xt ; θ)
7:
Analytically compute gradient dJ(θ)/dθ
8:
Gradient-based policy update (e.g. CG, BFGS)
9:
until convergence: θ∗ = arg min J(θ)
10:
Experiment: Execute π(xt ; θ∗ ); record {xt , ut }t=1...T
11:
12:
13:

Update dynamics model fˆ using all recorded data
until task learned
return π ∗ (xt ; θ∗ )

III. P OLICY S EARCH FOR PID C ONTROL
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. IIIB to Sec. III-D, followed by an analytic derivation of the
required cost function gradients in Sec. III-E.
A. PID Control Policies
The control output of a scalar PID controller is given by
Z t·∆T
e(τ )dτ + Kd ėt
(3)
ut = Kp et + Ki

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
fˆ. To incorporate all available data from policy rollouts

0

et = xdes,t − xt ,

(4)

1 A 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.

e1
e2

PID1
PID2

u1
u2

e1

PID1

Sys 1

+
e2

u

Sys 2

PID2

Fig. 3. Possible PID structures, exemplified with two controllers. Left:
Individual PID controllers acting on different system inputs. Right: Combination of PID controllers acting on the same system input.

where e(τ ) denotes the continuous error signal. The current
desired state xdes,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 (θPID = (Kp , Ki , Kd )).
A general PID control structure C(s) for MIMO processes
(1) can be described in transfer function notation by a D ×F
transfer function matrix


c11 (s) · · · c1D (s)

..  ,
..
(5)
C(s) =  ...
.
. 
cF 1 (s) · · ·

cF D (s)

where s denotes the complex Laplace variable and cij (s)
are of PID type. The multivariate error is given by et =
xdes,t − xt ∈ RD such that the multivariate input becomes
u(s) = C(s)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.
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 x0 , 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 zt that keeps track of the error
at the previous time step and the accumulated error,
zt := (xt , et−1 , ∆T

t−1
X

eτ ) ,

(6)

τ =0

where ∆T is the system’s sampling time. For simplicity,
we denote vectors as tuples (v1 , . . . , vn ), where vi may be
vectors themselves. The following augmentations are made
to obtain the necessary policy inputs:

1) Desired Goal State: The desired set-point or target
trajectory state2 xdes,t ∼ N (µdes,t , Σdes,t ) is independent of
zt yielding:



 

zt
µz
Σ
0
∼N
, z
.
(7)
xdes,t
µdes,t
0 Σdes,t
2) Error States: The current error is a linear function of
zt and xdes,t . The current error derivative and integrated error
are approximated by
et − et−1
,
∆T
t−1
X
e(τ )dτ ≈ ∆T
eτ + ∆T et .
ėt ≈

Z

t·∆T

0

(8)
(9)

τ =0

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 zt to provide the input
for a low-pass Finite Impulse Response (FIR) filter.
The fully augmented state is given by
z̃t := (zt , xdes,t , et , ∆T

t
X

eτ ,

τ =0

et − et−1
).
∆T

(10)

C. PID as Static State Feedback
Based on the augmented state z̃t , the PID control policy
for multivariate controllers can be expressed as a static state
feedback policy:
(3)

(4)

(5)

ut = APID (z̃t , z̃t , z̃t )
!
t
X
et − et−1
= APID et , ∆T
et ,
,
∆T
τ =0

(11)

(i)

where z̃t indicates the i-th term of (10). The specific
structure of the multivariate PID control law is defined by the
parameters in APID . For example, PID structures as shown
in Fig. 3 would be represented by


Kp,1
0
Ki,1
0
Kd,1
0
Aleft =
, (12)
0
Kp,2
0
Ki,2
0
Kd,2


Aright = Kp,1 Kp,2 Ki,1 Ki,2 Kd,1 Kd,2 . (13)
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 fˆ. PILCO
approximates the predictive distribution p(xt+1 ) by a Gaussian distribution using exact moment matching. From the
2 Drawing the desired state from a Gaussian distribution yields better
generalization to unseen targets as discussed in [13].

(3)

(4)

zt+1 = (xt+1 , z̃t , z̃t ) = (xt+1 , et , ∆T

t
X

eτ ) . (14)

ẍ (m/s2 )

dynamics model output xt+1 and the current error stored in
z̃t , the next state is obtained as

1
0
−1

0

z0 := (x0 , xdes,0 − x0 , 0) .

(16)

=:Et

To simplify the notation, we write dp(zt ) to denote the sufficient statistics derivatives dµt and dΣt of a Gaussian random
variable p(zt ) ∼ N (µt , Σt ) (analogous to the treatment in
[15]). The gradient of the immediate loss with respect to
the state distribution, dEt /dp(zt ), 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
dp(zt+1 )
δp(zt+1 ) dp(z̃t )
δp(zt+1 ) dp(xt+1 )
=
+
. (17)
dθ
δp(z̃t )
dθ
δp(xt+1 )
dθ
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 xt+1 is given by
dp(xt+1 )
δp(xt+1 ) dp(z̃t ) δp(xt+1 ) dp(ut )
=
+
.
dθ
δp(z̃t )
dθ
δp(ut )
dθ

(18)

The derivatives shown in red can be computed analytically
for the specific dynamics model [17]. Applying the chain
rule for the policy output p(ut ) obtained by (11) yields
dp(ut )
δp(ut ) dp(z̃t ) δp(ut )
=
+
,
dθ
δp(z̃t ) dθ
δθ

(19)

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
dp(z̃t ) dp(zt )
dp(z̃t )
=
dθ
dp(zt ) dθ

Fig. 4. Comparison of commanded (blue) and executed (red) end effector
acceleration for a policy rollout on Apollo. Additionally to the tracking
errors introduced by the joint level controllers, the system shows further
deterioration in acceleration tracking due to errors in the inverse dynamics
model and due to noise and delays.

(15)

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 gradientbased 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 cost3 derivative is obtained as
H
T
X
dEt dp(zt )
dJ(θ) X d
=
Ezt [c(zt )] =
.
{z
}
|
dθ
dθ
dp(zt ) dθ
t=1
t=1

1.5

Time t (s)

τ =0

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

1

0.5

(20)

3 We only address cost on the state z for simplicity. For practical
t
implementations, cost on u can be included by adding past inputs into
the system state as shown in Sec. IV.

Again, the part marked in blue is computed for the linear
transformation (10). Starting from the initial state where
dp(z0 )/dθ = 0, we obtain the gradients for all states with
respect to the policy parameters dp(zt )/dθ by iteratively
applying (17) to (20) for all time steps t.
IV. E XPERIMENTAL E VALUATION
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.
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 ut [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 xt comprises the end effector position x, velocity ẋ, pendulum angle
φ, and angular velocity φ̇.
During the rollouts, the commanded acceleration is limited
to umax = 3 m/s2 for safety reasons. Test policies are executed
for 20 s or interrupted once safety limits (xmax = 0.3 m,
θmax = 30◦ ) 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
c(et , ūt ) = 1 − exp(−(eTt Qet + ūTt Rūt )/2) .

(21)

For balancing the pendulum in the upright position the
desired trajectory is given by xdes,t = 0. Weights are set to
Q = diag(1/0.22 , 1/0.022 ) for end effector and pendulum
position error, and to R = 1/0.42 for the control input.
The selected cost function saturates quickly for x > w if
x is weighted by 1/w2 . 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.

300

300

We employ a PID controller on the position error ex =
xdes,t −x and a PD controller acting on the pendulum angular
error eθ = θdes,t − θ. 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 θ =
(Kp,x , Ki,x , Kd,x , Kp,θ , Kd,θ ). 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.

250

250

200

200

150

150

100

100

50

50

0

0

C. Modifications for Dynamics Model Learning
When first training GP dynamics models xt+1 = fˆ(xt , ut )
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.
1) 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 ∆xt = xt+1 − xt .
We compute a sparse GP using Snelson’s approximation [24]
having a covariance parametrized by 400 inducing pseudoinput points. For both GP models, hyperparameters θGP,i =
(l1 , l2 , σf , σ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.
2) 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
xt+1 = fˆ(xt , . . . , xt−n , ut , . . . , ut−m ) .

(22)

Cost-to-go J (-)

B. PID Control Structure Setup

1 2 3 4 5 6 7
Iteration (-)

iteration
iteration
iteration
iteration

1
2
3
4

iteration 6
iteration 5 and 7
0

5
10
20
15
Gradient-based policy updates (-)

Fig. 5. Expected cost-to-go optimization results. Left: Iterative improvement in predicted loss (blue) compared to the cost observed in a single
robot experiment (red). Right: Optimization results for each iteration.

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 xNARX :=
(xt , . . . , xt−3 , φt , . . . , φt−2 , ut , . . . , ut−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.
3) 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.
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 urnd,t ∼ N (0, 1 m/s2 )
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

0

−0.25

Position x (m), angle θ (rad)

0.25

−0.5

0.5

0.5
Position x (m), angle θ (rad)

Position x (m), angle θ (rad)

0.5

0.25
0

0.5

1
Time t (s)

1.5

2

−0.5

0

−0.25

−0.25

0

0.25

1

0

2

3

−0.5

0

1

3
2
Time t (s)

Time t (s)

(a) Iteration 1

(b) Iteration 3

4

5

(c) Iteration 7

Fig. 6. Predicted system behavior (dashed lines, error-bars indicate 95 % confidence intervals) and experimental results (solid lines) visualized for the first
(a), an intermediate (b) and the final iteration (c) of the policy learning process. The end effector position (blue) and the angular position (red) are shown.

Position x (m)
angle θ (rad)

0.05
0
−0.05
−0.1

Disturbance in pendulum angle
0

2

4
6
Time t (s)

8

10

Fig. 7. Disturbance rejection of the optimized PID policy. End effector
position (blue) and pendulum angle (red) display the closed loop response
of the optimized PID policy to manually introduced disturbances.

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

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
X ∼ N (µX , ΣX ) ∈ RD is given by
Y = AX + b = N (µY , ΣY ) ∈ RP ,

(23)

where A ∈ RP ×D and b ∈ RP ×1 . The joint probability
distribution of X and Y is given by
 
  

ΣX
ΣX C
X
µ
∼ N( X ,
)
(24)
Y
µY
CT ΣTX
ΣY
where
µY = AµX + b,

ΣY = AΣX AT ,

C = AT

(25)

The non-zero partial derivatives of Y ’s sufficient statistics
are given by
δµY
= A ∈ RP ×D ,
δµX
δµY
= µTX ⊗ I ∈ RP ×DP ,
δA

2
2
δΣY
= A ⊗ A ∈ RP ×D
δΣX
(26)
δµY
= I ∈ RP ×P
(27)
δb

V. C ONCLUSION
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.

2 δ(ΣY )
2
δCkl
δC
= δil δkj ,
∈ RP D×D ,
∈ RP ×P D
δAij
δA
δA
(28)
δ(ΣY )kl
= δlj (AT ΣX )ki + δkj (ΣX A)il
(29)
δAij

where ⊗ is the Kronecker product and δij is the Kronecker
delta (cf. [27] for useful matrix derivatives).

R EFERENCES
[1] P. Cominos and N. Munro, “PID controllers: Recent tuning methods
and design to specification,” IET Proceedings on Control Theory and
Applications, vol. 149, no. 1, pp. 46–53, 2002.
[2] F. Jiang and Z. Gao, “An application of nonlinear PID control to a
class of truck ABS problems,” in IEEE Conference on Decision and
Control, vol. 1, 2001, pp. 516–521.
[3] B. Siciliano, L. Sciavicco, L. Villani, and G. Oriolo, Robotics: modelling, planning and control. Springer Science & Business Media,
2010.
[4] A. O’Dwyer, Handbook of PI and PID controller tuning rules. World
Scientific, 2009, vol. 57.
[5] K. J. Åström and T. Hägglund, Advanced PID control. ISA-The
Instrumentation, Systems and Automation Society, 2006.
[6] J. Berner, T. Hägglund, and K. J. Åström, “Asymmetric relay autotuning - practical features for industrial use,” Control Engineering
Practice, vol. 54, pp. 231–245, 9 2016.
[7] M. A. Johnson and M. H. Moradi, PID control. Springer, 2005.
[8] T. Yamamoto and S. Shah, “Design and experimental evaluation of a
multivariable self-tuning PID controller,” IET Proceedings on Control
Theory and Applications, vol. 151, no. 5, pp. 645–652, 2004.
[9] M. Deisenroth and C. E. Rasmussen, “PILCO: A model-based and
data-efficient approach to policy search,” in Proceedings of the International Conference on Machine Learning (ICML), 2011, pp. 465–
472.
[10] J. Kocijan, A. Girard, B. Banko, and R. Murray-Smith, “Dynamic
systems identification with Gaussian processes,” Mathematical and
Computer Modelling of Dynamical Systems, vol. 11, no. 4, pp. 411–
424, 2005.
[11] R. Murray-Smith, D. Sbarbaro, C. E. Rasmussen, and A. Girard,
“Adaptive, cautious, predictive control with Gaussian process priors,”
13th IFAC Symposium on System Identification, pp. 1195–1200, 2003.
[12] R. Murray-Smith and D. Sbarbaro, “Nonlinear adaptive control using
non-parametric Gaussian process prior models,” 15th Triennial World
Congress of the International Federation of Automatic Control (IFAC),
2002.
[13] M. P. Deisenroth, C. E. Rasmussen, and D. Fox, Learning to control
a low-cost manipulator using data-efficient reinforcement learning.
MIT Press Journals, 2011, vol. 7, pp. 57–64.
[14] B. Bischoff, D. Nguyen-Tuong, T. Koller, H. Markert, and A. Knoll,
“Learning throttle valve control using policy search,” in Machine

[15]

[16]
[17]

[18]
[19]
[20]
[21]

[22]

[23]
[24]
[25]
[26]
[27]

Learning and Knowledge Discovery in Databases. Springer, 2013,
pp. 49–64.
M. P. Deisenroth, D. Fox, and C. E. Rasmussen, “Gaussian processes
for data-efficient learning in robotics and control,” IEEE Transactions
on Pattern Analysis and Machine Intelligence, vol. 37, no. 2, pp. 408–
423, 2015.
C. K. Williams and C. E. Rasmussen, Gaussian processes for machine
learning. MIT Press, 2005.
J. Q. Candela, A. Girard, J. Larsen, and C. E. Rasmussen, “Propagation
of uncertainty in Bayesian kernel models-application to multiplestep ahead forecasting,” in Proceedings of the IEEE International
Conference on Acoustics, Speech, and Signal Processing (ICASSP),
vol. 2, 2003, pp. II–701.
M. P. Deisenroth and D. Fox, “Multiple-target reinforcement learning
with a single policy,” in ICML Workshop on Planning and Acting with
Uncertain Models. Citeseer, 2011.
R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction.
MIT press Cambridge, 1998, vol. 1, no. 1.
C. W. Anderson, “Learning to control an inverted pendulum using
neural networks,” IEEE Control Systems Magazine, vol. 9, no. 3, pp.
31–37, 1989.
L. Righetti, M. Kalakrishnan, P. Pastor, J. Binney, J. Kelly, R. C.
Voorhies, G. S. Sukhatme, and S. Schaal, “An autonomous manipulation system based on force control and optimization,” Autonomous
Robots, vol. 36, no. 1-2, pp. 11–30, 2014.
A. Marco, P. Hennig, J. Bohg, S. Schaal, and S. Trimpe, “Automatic
LQR tuning based on Gaussian process global optimization,” in
Proceedings of the IEEE International Conference on Robotics and
Automation (ICRA), 2016, pp. 270–277.
S. Trimpe and R. DAndrea, “The balancing cube: A dynamic sculpture
as test bed for distributed estimation and control,” IEEE Control
Systems, vol. 32, no. 6, pp. 48–75, Dec 2012.
E. Snelson and Z. Ghahramani, “Sparse Gaussian processes using
pseudo-inputs,” in Advances in Neural Information Processing Systems, 2005, pp. 1257–1264.
S. A. Billings, Nonlinear system identification: NARMAX methods in
the time, frequency, and spatio-temporal domains. John Wiley &
Sons, 2013.
K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: NSGA-II,” IEEE Transactions on
Evolutionary Computation, vol. 6, no. 2, pp. 182–197, 2002.
K. B. Petersen, M. S. Pedersen, et al., “The matrix cookbook,”
Technical University of Denmark, vol. 7, p. 15, 2008.

