Abstract
In reinforcement learning (RL), an autonomous agent learns to perform complextasks by maximizing an exogenous reward signal while interacting with itsenvironment. In realworld applications, test conditions may differsubstantially from the training scenario and, therefore, focusing on purereward maximization during training may lead to poor results at test time. Inthese cases, it is important to tradeoff between performance and robustnesswhile learning a policy. While several results exist for robust, modelbasedRL, the modelfree case has not been widely investigated. In this paper, wecast the robust, modelfree RL problem as a multiobjective optimizationproblem. To quantify the robustness of a policy, we use delay margin and gainmargin, two robustness indicators that are common in control theory. We showhow these metrics can be estimated from data in the modelfree setting. We usemultiobjective Bayesian optimization (MOBO) to solve efficiently thisexpensivetoevaluate, multiobjective optimization problem. We show thebenefits of our robust formulation both in simtoreal and pure hardwareexperiments to balance a Furuta pendulum.
Quick Read (beta)
Robust Modelfree Reinforcement Learning with Multiobjective Bayesian Optimization
Abstract
In reinforcement learning (RL), an autonomous agent learns to perform complex tasks by maximizing an exogenous reward signal while interacting with its environment. In real world applications, test conditions may differ substantially from the training scenario and, therefore, focusing on pure reward maximization during training may lead to poor results at test time. In these cases, it is important to tradeoff between performance and robustness while learning a policy. While several results exist for robust, modelbased RL, the modelfree case has not been widely investigated. In this paper, we cast the robust, modelfree RL problem as a multiobjective optimization problem. To quantify the robustness of a policy, we use delay margin and gain margin, two robustness indicators that are common in control theory. We show how these metrics can be estimated from data in the modelfree setting. We use multiobjective Bayesian optimization (MOBO) to solve efficiently this expensivetoevaluate, multiobjective optimization problem. We show the benefits of our robust formulation both in simtoreal and pure hardware experiments to balance a Furuta pendulum.
equation \crefnamesectionSec.Sec. \crefnamefigureFig.Fig. \crefnamelemmaLemmaLemma
I Introduction
In reinforcement learning (RL) [1], the goal is to learn a controller to perform a desired task from the data produced by the interaction between the learning agent and its environment. In this framework, autonomous agents are trained to maximize their return. It is common to assume that such agents will be deployed in conditions that are similar, if not equal, to those they were trained in. In this case, a returnmaximizing agent performs well at test time. However, in real world applications, this assumption may be violated. For example, in robotics, we can use RL to learn to fly a drone indoor. However, later on we may use the same drone to carry a payload in a windy environment. The new environmental conditions and the possible deterioration of the drone components due to their usage may result in a poor, if not catastrophic, performance of the learned controller. Another scenario where training and testing conditions differ substantially is the simtoreal setting, i.e., when we deploy a controller trained in simulation on a realworld agent.
Considering robustness alongside performance when learning a controller can limit performance degradation due to different training and testing environments. In special cases, these goals may be aligned, and a highperforming controller can also be robust. This is the case for the Linear Quadratic Regulator (LQR), a linear statefeedback controller that is optimal for the case of linear dynamics, quadratic cost, and perfect state measurements. It is wellknown that the LQR exhibits strong robustness indicators, such as gain and phase margins [2]. While performance and robustness go hand in hand for the LQR, they are often conflicting in other cases. For example, a celebrated result in control theory shows that the Linear Quadratic Gaussian (LQG) regulator  the noisy counterpart of the LQR  can be arbitrarily close to instability, despite being optimal [3]. Thus, in general, we need to tradeoff between performance and robustness [4].
Contributions. While many works investigating the performance/robustness tradeoff exist in both the RL and control theory literature for the modelbased setting, few results are known for the modelfree scenario. However, there are several realworld scenarios where models are not available, inaccurate, or too expensive to use, but robustness is fundamental. Thus, in this paper, we introduce the first dataefficient, robust, modelfree RL method based on policy optimization with multiobjective Bayesian optimization (MOBO). In particular, these are our contributions:

•
We formulate the robust, modelfree RL as a multiobjective optimization problem.

•
We propose a modelfree, datadriven evaluation of delay and gain margins, two common robustness indicators from the modelbased setting (where they are computed analytically).

•
We solve this problem efficiently with expected hypervolume improvement (EHI).

•
We introduce the first method that can learn robust controllers directly on hardware in a modelfree fashion.

•
We show how our approach outperforms nonrobust policy optimization in evaluations on a Furuta pendulum for both a simtoreal and a pure hardware setting.
Related work. Robustness has been widely investigated in control theory [5], and standard robust control techniques for linear systems include loop transfer recovery [6], ${H}_{\mathrm{\infty}}$ control, and $\mu $ synthesis [5, 7]. However, these methods typically assume the availability of a model, and none of these includes a learning component. Recently, robustness has drawn attention in datadriven settings, giving rise to the field of robust, modelbased RL. Robust Markov decision processes study the RL problem when the transition model is subject to known and bounded uncertainties. For example, [8] studies the dynamic programming recursion in this setting. Other methods that consider parametric uncertainties include [9, 10]. All the previous methods are modelbased.
Robustness and performance are typical objectives in control design, which often conflict each other, thus requiring design tradeoffs [4, 11]. In the model free literature, this tradeoff is often fixed a priori and the resulting problem is solved with standard optimization methods. In [12] a weighted cost that balances performance and robustness is optimized. In [13] robust controllers are learned via gradient ascent with random multiplicative noise on the control action. In [14, 15] external, adversarial disturbances are used instead. In these works, the upper bound on the magnitude of the disturbance implicitly balances robustness and performance. However, setting this tradeoff is often not intuitive and, in case the requirements are misspecified or updated, a new controller must be learned. Alternatively, robust control design methods based on multiobjective optimization explore the spectrum of such tradeoffs. The work in [16] gives a review of such methods, with a focus on genetic algorithms, which, due to their low data efficiency require the model to compute the robustness indices.
Modelfree RL algorithms are typically validated in simulations due to their high sample complexity. However, in robotics, it is crucial to test these methods on hardware. Bayesian optimization (BO) [17, 18] has been successfully applied to learn lowdimensional controllers for hardware systems. For example, [19] learns to control the $x$coordinate for a quadrotor hovering task with a linear controller, [20] learns a linear state feedback controller for a cartpole system in a simtoreal setting and [21, 22] tune the parameters of adhoc controllers for locomotion tasks. However, none of these methods considers robustness, making ours the first one to learn robust controllers from data directly on hardware.
MOBO is the branch of BO that solves multiobjective problems. MOBO algorithms include EHI [23], $\u03f5$PAL [24], and PESMO [25]. They have been applied to several tasks including trading off prediction speed and accuracy in machine learning models. However, they have rarely been applied to RL. To the best of our knowledge, this has been done only in [26, 27], where a tradeoff between frontal camera movement and forward speed is found for a snakelike robot, for homoschedastic and heteroschedastic noise respectively. Robustness is not explicitly treated in these works.
II Problem Statement
In this section, we introduce our formulation of robust modelfree RL as a multiobjective optimization problem. For ease of exposition, we limit ourselves to two objectives. However, this approach naturally extends to the any number of objectives, for example, multiple robustness indicators.
We assume we have a system with unknown dynamics, $h$, and unknown observation model, $g$,
$${\mathbf{x}}_{t+1}=h({\mathbf{x}}_{t},{\mathbf{u}}_{t},{\mathbf{w}}_{t}),{\mathbf{o}}_{t}=g({\mathbf{x}}_{t},{\mathbf{v}}_{t})$$  (1) 
where $\mathbf{x}$ is the state, $\mathbf{u}$ is the control input, $\mathbf{o}$ is the observation and $\mathbf{w}$ and $\mathbf{v}$ are the process and sensor noise. An RL agent aims at learning a controller ${\mathbf{u}}_{t}=\pi ({\mathbf{o}}_{t}\theta )$, i.e., a mapping parametrized by $\theta $ from an observation ${\mathbf{o}}_{t}$ to an action ${\mathbf{u}}_{t}$ that allows it to complete its task. Policy optimization algorithms are a class of modelfree RL methods that solve this problem by optimizing the performance of a given controller for the task at hand as a function of the parameters $\theta $. Concretely, given a performance metric ${f}_{1}:\mathrm{\Theta}\to \mathbb{R}$, standard, nonrobust policy optimization algorithms aim to find ${\theta}^{*}\in argmax{f}_{1}(\theta )$. In this work, we consider regulation tasks, i.e., bringing and keeping the system in a desired goal state $\overline{\mathbf{x}}$. This includes common problems like stabilization, setpoint tracking, or disturbance rejection. The performance indicator ${f}_{1}$ encodes these objectives.
To extend this framework to the robustnessaware case, we use a second function ${f}_{2}:\mathrm{\Theta}\to \mathbb{R}$ that measures the robustness of a controller. Since both the dynamics $h$ and the observation model $g$ are unknown, we must evaluate or approximate the value of ${f}_{2}$ from data. In creftype IIIB, we introduce the gain and the delay margin, two alternatives for ${f}_{2}$ that are commonly used in modelbased control and we discuss how to evaluate them in the modelfree setting.
We aim at finding the best controller in terms of performance and robustness, as measured by ${f}_{1}$ and ${f}_{2}$. However, since we compare controllers based on multiple, and possibly conflicting, criteria, we cannot define a single best controller. Given a controller $\theta $, we denote with ${\mathbf{f}}_{\theta}=[{f}_{1}(\theta ),{f}_{2}(\theta )]$ the array containing its performance and robustness values. To compare two controllers ${\theta}_{1}$ and ${\theta}_{2}$, we use the canonical partial order over ${\mathbb{R}}^{2}$: ${\mathbf{f}}_{{\theta}_{1}}\u2ab0{\mathbf{f}}_{{\theta}_{2}}$ iff ${f}_{i}({\theta}_{1})\ge {f}_{i}({\theta}_{2})$ for $i=1,2$. This induces a relation in the controller space $\mathrm{\Theta}$: ${\theta}_{1}\u2ab0{\theta}_{2}$ iff ${\mathbf{f}}_{{\theta}_{1}}\u2ab0{\mathbf{f}}_{{\theta}_{2}}$. If ${\theta}_{1}\u2ab0{\theta}_{2}$, we say that ${\theta}_{1}$ dominates ${\theta}_{2}$. The Pareto set ${\mathrm{\Theta}}^{*}\subseteq \mathrm{\Theta}$ is the set of nondominated points in the domain, i.e., ${\theta}^{*}\in {\mathrm{\Theta}}^{*}$ iff $\exists i=1,2$ such that ${f}_{i}({\theta}^{*})>{f}_{i}(\theta )$ for all $\theta \in \mathrm{\Theta}$. The Pareto front is the set of function values corresponding to the Pareto set. The Pareto set is optimal in the sense that, for each point ${\theta}^{*}$ in it, it is not possible to find another point in the domain that improves the value of one objective without degrading another [28]. The goal of this paper is to approximate ${\mathrm{\Theta}}^{*}$ from data.
creftype 1 represent our problem graphically: we suggest a controller, we evaluate its performance and robustness on the system and we select a new controller based on these observations to find an approximation of the Pareto front.
III Learning the Performancerobustness Tradeoff
For the robust, modelfree RL setting we consider, we propose to learn the Pareto front characterizing the performancerobustness tradeoff of a given system with MOBO. Here, we describe the necessary components to solve our problem in a data efficient way: MOBO and the robustness and performance indicators used in our experiments. Moreover, we discuss how to evaluate such indicators from data in a modelfree fashion.
IIIA Multiobjective Bayesian optimization
MOBO algorithms solve multiobjective optimization problems by sequentially querying the objective at different inputs and obtaining noisy evaluations of the corresponding values. They build a statistical model of the objectives to capture the belief over them given the data available. They measure how informative a point in the domain is about the problem solution with an acquisition function. At every iteration, they evaluate the objective at the most informative point, as measured by the acquisition function. Thus, the complex multiobjective optimization problem is decomposed into a sequence of simpler scalarvalued optimization problems. In the following, we describe the surrogate model and the acquisition function used in this work.
Intrinsic Model of Coregionalization A singleoutput Gaussian process (GP) [29] is a probability distribution over the space of functions of the form $f:\mathrm{\Theta}\to \mathbb{R}$, such that the joint distribution of the function values computed over any finite subset of the domain follows a multivariate Gaussian distribution. A GP is fully specified by a mean function $\mu :\mathrm{\Theta}\to \mathbb{R}$, which, w.l.o.g., is usually assumed to be zero, $\mu (\theta )=0$ for all $\theta \in \mathrm{\Theta}$, and a covariance function, or kernel, $k:\mathrm{\Theta}\times \mathrm{\Theta}\to \mathbb{R}$. The kernel encodes the strength of statistical correlation between two latent function values and, therefore, it expresses our prior belief about the function behavior.
Similarly, a $D$output GP is a probability distribution over the space of functions of the form $\mathbf{f}:\mathrm{\Theta}\to {\mathbb{R}}^{D}$. The difference with respect to singleoutput GPs is that, in this case, the kernel must capture the correlation across different output dimensions in addition to the correlation of function values at different inputs. The simplest way of doing this is by assuming that each output is independent. However, this model disregards the fundamental tradeoff between robustness and performance that we are considering. For a review on kernels for multioutput GPs, see [30]. In this work, we use the intrinsic model of coregionalization (ICM), which defines the covariance between the ${i}^{\text{th}}$ value of $\mathbf{f}(\theta )$ and the ${j}^{\text{th}}$ value of $\mathbf{f}({\theta}^{\prime})$ by separating the input and the output contribution as follows, ${b}_{ij}k(\theta ,{\theta}^{\prime})$. In this case, we say $\mathbf{f}\sim \mathcal{G}\mathcal{P}(\bm{\mu}(\cdot ),\mathbf{K}(\cdot ,\cdot )=\mathbf{B}k(\cdot ,\cdot ))$, where $\bm{\mu}:\mathrm{\Theta}\to {\mathbb{R}}^{D}$ is a $D$dimensional mean function, $k:\mathrm{\Theta}\times \mathrm{\Theta}\to \mathbb{R}$ is a scalarvalued kernel and $\mathbf{B}\in {\mathbb{R}}^{D\times D}$ is a matrix describing the correlation in the output space (more details on $\mathbf{B}$ in creftype IV). Given $N$ noisy observations of $\mathbf{f}$, $\mathcal{D}=\{({\theta}_{1},{\mathbf{y}}_{1}),\mathrm{\cdots}({\theta}_{N},{\mathbf{y}}_{N})\}$, with ${\mathbf{y}}_{i}=\mathbf{f}({\theta}_{i})+{\omega}_{i}$, where ${\omega}_{i}\sim \mathcal{N}(\mathrm{\U0001d7ce},\mathrm{\Sigma})$ is i.i.d. Gaussian noise, we can compute the posterior distribution of the function values conditioned on $\mathcal{D}$ at a target input ${\theta}^{\star}$ in closed form as $p(\mathbf{f}({\theta}^{\star})\mathcal{D},{\theta}^{\star})\sim \mathcal{N}({\mathbf{f}}^{\star}({\theta}^{\star}),{\mathbf{K}}^{\star}({\theta}^{\star},{\theta}^{\star})))$. We denote with $\bm{\theta}$ the inputs contained in $\mathcal{D}$ and with $\mathbf{K}(\bm{\theta},\bm{\theta})$ the $ND\times ND$ matrix with entries ${(\mathbf{K}({\theta}_{u},{\theta}_{v}))}_{i,j}$ for $u,v=1,\mathrm{\dots},N$ and $i,j=1,\mathrm{\dots},D$, then
${\mathbf{f}}^{\star}({\theta}^{\star})={\mathbf{K}}_{{\theta}^{\star}}^{\top}{(\mathbf{K}(\bm{\theta},\bm{\theta})+\mathbf{\Sigma})}^{1}\overline{\mathbf{y}},$  (2)  
${\mathbf{K}}^{\star}({\theta}^{\star},{\theta}^{\star})=\mathbf{K}({\theta}^{\star},{\theta}^{\star}){\mathbf{K}}_{{\theta}^{\star}}{(\mathbf{K}(\bm{\theta},\bm{\theta})+\mathbf{\Sigma})}^{1}{\mathbf{K}}_{{\theta}^{\star}}^{\top},$  (3) 
where $\mathbf{\Sigma}=\mathrm{\Sigma}\otimes {I}_{N}$, with $\otimes $ denoting the Kronecker product, ${\mathbf{K}}_{{\theta}^{\star}}\in {\mathbb{R}}^{D\times ND}$ has entries ${(\mathbf{K}({\theta}^{\star},{\theta}_{v}))}_{i,j}$ for $v=1,\mathrm{\dots},N$ and $i,j=1,\mathrm{\dots},D$ and $\overline{\mathbf{y}}$ is the $ND$dimensional vector containing the concatenation of the observations in $\mathcal{D}$.
Expected Hypervolume Improvement EHI is an acquisition function introduced in [23], which selects inputs to evaluate based on a notion of improvement with respect to the incumbent solution. In multiobjective optimization, incumbent solutions take the form of approximations of the Pareto set, ${\mathcal{X}}^{*}$, whose quality is measured by the hypervolume indicator induced by the corresponding front, ${\mathcal{Y}}^{*}$ with respect to a reference $r$. Formally, the hypervolume indicator of a set of points $A$ with respect to a reference $r$, $\text{HV}(A;r)$, is the Lebesgue measure of the hypervolume covered by the boxes that have an element in $A$ as upper corner and the reference as lower corner. It quantifies the size of the portion of the output space that is Paretodominated by the points in $A$. Given an estimate of the Pareto front, ${\mathcal{Y}}^{*}$, the hypervolume improvement of $\theta \in \mathrm{\Theta}$ is defined as the relative improvement in hypervolume obtained by adding the function value at $\theta $, $f(\theta )$, to ${\mathcal{Y}}^{*}$, $\text{HI}(f(\theta );{\mathcal{Y}}^{*},r)=\text{HV}({\mathcal{Y}}^{*}\cup \mathbf{f}(\theta );r)\text{HV}({\mathcal{Y}}^{*};r).$ However, we do not know $\mathbf{f}(\theta )$. Instead, we have a belief over its value expressed by the posterior distribution of the GP, which, in turn, induces a distribution over the hypervolume improvement corresponding to an input $\theta $. The EHI acquisition function quantifies the informativeness of an input $\theta $ toward the solution of the multiobjective optimization problem through the expectation of this distribution,
$$\alpha (\theta \mathcal{D},{\mathcal{Y}}^{*},r)={\int}_{f(\theta )\in {\mathbb{R}}^{n}}\text{HI}(\mathbf{f}(\theta );{\mathcal{Y}}^{*}),r)p(\mathbf{f}(\theta )\mathcal{D})d\mathbf{f}(\theta ).$$  (4) 
[23] shows how to compute the integral in creftype 4 in closed form.
IIIB Robustness
œ¬In general, robustness can have very different meanings. One may desire to ensure robustness to a certain class of disturbances, imperfections in the control system, or uncertainty in the process, for example. In control theory, the latter is often understood as robustness in the stricter sense. Specifically, robust stability assures that a controller stabilizes every member from a set of uncertain processes [5]. Such processes can, for example, be defined through a nominal process and variations thereof. Different variations lead to different robustness characterizations. Likewise, there are different notions of stability that are meaningful depending on the context. For example, for a deterministic system, asymptotic stability, i.e., ${\mathbf{x}}_{t}\to \overline{\mathbf{x}}$ as $t\to \mathrm{\infty}$, where $\overline{\mathbf{x}}$ is an equilibrium of the system, is often used; for systems that are continuously excited, e.g., through noise, and thus cannot approach $\overline{\mathbf{x}}$, one may seek the above limit to hold in expectation or practical stability in the sense of a bounded state, i.e., $\parallel {\mathbf{x}}_{t}\overline{\mathbf{x}}\parallel \le {\mathbf{x}}_{\mathrm{max}}$ for all $t\ge \overline{t}$. A controller is unstable when the respective condition does not hold (e.g., no asymptotic convergence, or ${\mathbf{x}}_{t}$ grows beyond any bounds).
While many sophisticated robustness metrics have been developed, stability margins such as gain and delay margins are some of the most common and intuitive ones [11, Sec. 9.3]. We consider these in this work and comment on alternatives in creftype V. Below, we formally introduce them and we explain how to evaluate them in a modelfree setting. Notice that, our datadriven definitions can be extended to any setting where a success/failure outcome can be defined and, therefore, are not limited to stability considerations.
Gain margin. In classical control, the upper (lower) gain margin is defined for singleinputsingleoutput (SISO) linear systems as the largest factor ${\kappa}_{\text{max}}\in (1,\mathrm{\infty})$ (the smallest factor ${\kappa}_{\text{min}}\in (0,1)$) that can multiply the openloop transfer function so that the closedloop system is stable [31, Sec. 9.5]. As the openloop transfer function encodes both the process and the controller dynamics, the factor may represent uncertainty in the process gain or the actuator efficiency, for example. In this work, we consider a factor $\kappa $ to be multiplied by the control action (i.e., ${\mathbf{u}}_{t}=\kappa \times \pi ({\mathbf{o}}_{t}\theta )$), which is equivalent to the definition for linear SISO systems, but can also be used for nonlinear ones. It quantifies how much we can lower/amplify the control action before making the system unstable. In a way, it quantifies how “far” we are from instability and, thus, how much we can tolerate differences between training and testing.
Delay margin. Similarly, we define the delay margin as the largest time delay on the measurement ${\mathbf{o}}_{t}$ such that the controlled system is still stable. Formally, it is the largest value of $d\in (0,\mathrm{\infty})$ such that the closedloop system with the delayed control action ${\mathbf{u}}_{t}=\pi ({\mathbf{o}}_{td}\theta )$ is stable. As delay in data transmission between sensor, controller, and actuator, and in the control computation are present in most control systems, the delay margin is a very relevant measure.
Estimate from data. While the indicators above can be readily computed for linear systems, they are difficult to compute analytically if the model is nonlinear, or impossible if no model is available, as considered herein. We describe an experiment to estimate the delay margin from data in a modelfree setting (those for the gain margins are analogous). For general nonlinear systems, stability with respect to an equilibrium is a local property. Thus, we assume we can reset the system to a state in the neighborhood of the equilibrium of interest, i.e., we have ${\mathbf{x}}_{0}\in \mathcal{B}(\overline{\mathbf{x}},\rho )$, where $\mathcal{B}(\overline{\mathbf{x}},\rho )$ is a ball centered at $\overline{\mathbf{x}}$ of radius $\rho $. We can establish whether the delay margin, denoted with ${d}^{*}$, is larger or smaller than a delay $d$ by resetting the system near $\overline{\mathbf{x}}$, deploying the delayed controller $\pi ({\mathbf{o}}_{td}\theta )$ and evaluating the stability of the resulting trajectory.
[t] \[email protected]@algorithmic [1] \STATEInputs: Reference $r$ \STATE${\mathcal{D}}_{0}\leftarrow \mathrm{\varnothing},{\mathcal{Y}}_{0}^{*}\leftarrow \{r\}$ \FOR$k=0,1,\mathrm{\dots}$ \STATE${\theta}_{k+1}\leftarrow {\text{argmax}}_{\theta \in \mathrm{\Theta}}\text{EHI}(\theta {\mathcal{D}}_{k},{\mathcal{Y}}_{k}^{*},r)$ \STATE${y}_{k+1}^{1}\leftarrow \text{\U0001d699\U0001d68e\U0001d69b\U0001d68f\U0001d698\U0001d69b\U0001d696\U0001d68a\U0001d697\U0001d68c\U0001d68e}\mathrm{\_}\text{\U0001d68e\U0001d6a1\U0001d699\U0001d68e\U0001d69b\U0001d692\U0001d696\U0001d68e\U0001d697\U0001d69d}({\theta}_{k+1})$ \STATE${y}_{k+1}^{2}\leftarrow \text{\U0001d69b\U0001d698\U0001d68b\U0001d69e\U0001d69c\U0001d69d\U0001d697\U0001d68e\U0001d69c\U0001d69c}\mathrm{\_}\text{\U0001d68e\U0001d6a1\U0001d699\U0001d68e\U0001d69b\U0001d692\U0001d696\U0001d68e\U0001d697\U0001d69d}({\theta}_{k+1})$ \STATE${\mathcal{D}}_{k+1}\leftarrow {\mathcal{D}}_{k}\cup \{({\theta}_{k+1},[{y}_{k+1}^{1},{y}_{k+1}^{2}])\}$ \STATE${\mathrm{\Theta}}_{k+1}^{*},{\mathcal{Y}}_{k+1}^{*}\leftarrow \text{\U0001d67f\U0001d68a\U0001d69b\U0001d68e\U0001d69d\U0001d698}\mathrm{\_}\text{\U0001d69c\U0001d68e\U0001d69d}\mathrm{\_}\text{\U0001d68a\U0001d697\U0001d68d}\mathrm{\_}\text{\U0001d68f\U0001d69b\U0001d698\U0001d697\U0001d69d}({\mathcal{D}}_{k+1})$ \ENDFOR\STATEOutputs: Pareto set ${\mathrm{\Theta}}^{*}$, Pareto front ${\mathcal{Y}}^{*}$
In practice, two problems arise with this approach: (i) we can evaluate a finite number of delays with a finite number of experiments; and (ii) while stability is an asymptotic condition on the state, we do not know the state and we run finite experiments. The first problem requires us to select carefully the delays we evaluate. We know that increasing values of delay take a stable system closer to instability. Thus, given $m$ delays $[{d}_{1},\mathrm{\cdots},{d}_{m}]$, we do a binary search to find the largest one for which the closedloop system is stable. This allows us to approximate the delay margin with $\mathrm{log}m$ experiments. Not knowing the state $\mathbf{x}$ can be solved by estimating it from the noisy sensor measurements ${\mathbf{o}}_{t}$ or by introducing a new definition of stability based on ${\mathbf{o}}_{t}$ rather than ${\mathbf{x}}_{t}$. Concerning the finite trajectories we note that, in practical cases, it is rare to have small compounding deviations from $\overline{\mathbf{x}}$ resulting in a divergent behavior emerging only in the long run. Often, a controller makes the system converge to or diverge from $\overline{\mathbf{x}}$ within a short amount of time. In our experiments, we say that a controller stabilizes the system if, after a burnin time that accounts for the transient behavior, it keeps the state within a box around $\overline{\mathbf{x}}$. Controllers with good margins are investigated further with longer experiments to eliminate potential outliers due to the finite trajectory issue.
Reliably estimating robustness indicators and stability of a system without a model is challenging. The estimation technique we presented is intuitive and easy to implement. While it does not provide formal guarantees on the estimation error, we show in creftype IV that it is accurate enough to greatly improve the robustness of our algorithm with respect to the nonrobust policy optimization baseline.
IIIC Algorithm
In this section, we describe our robust policy optimization algorithm, for the pseudocode see creftype IIIB. At iteration $k$, we select the controller that maximizes the EHI criterion. Then, we run two experiments to estimate its performance and robustness. For the performance, we introduce a state and action dependent reward and we define the return as the average reward obtained over an episode. The performance index is defined as the expectation of the return, which we approximate with a Monte Carlo estimate over multiple episodes. To estimate the robustness, we use the experiments from creftype IIIB. We update the data set with the experiments results. Finally, we update the estimate of the Pareto front that is used to compute the EHI as the set of dominating points of the data set. Other options to compute such estimate from the posterior of the GP exist. However, they are computationally more expensive and they resulted in a similar performance in our experiments. In the end, the algorithm returns an estimate of the Pareto set and front. The choice of a controller from the Pareto set depends on the performancerobustness tradeoff required by the test applications and, therefore, the choice is left to the practitioner.
IV Experimental Results
We compare the robust policy optimization algorithm in creftype IIIB to its nonrobust counterpart based on scalar BO as, e.g., in [19, 32]. We use the scalar equivalent of EHI for the nonrobust case, i.e., the expected improvement (EI) algorithm [17]. We present two set of experiments: training controllers in simulation and directly on hardware, respectively. In both cases, the learned controllers are tested on the hardware in a set of different conditions.
System. We learn a controller for a Furuta pendulum [33] (see Figure 1), a system that is closely related to the wellknown cart pole. It replaces the cart with a rotary arm that rotates on the horizontal plane. In our experiments, we use the Qube Servo 2 by Quanser [34], a small scale Furuta pendulum. It uses a brushless DC motor to exert a torque on the rotary arm, and it is equipped with two incremental optical encoders with 2048 counts per revolution to measure the angle of the rotary arm and the pendulum. For simtoreal, we use the dynamics model provided in the Qube Servo 2 manual [34], which is a nonlinear rigid body model. A more detailed model is presented in [33].
Controller. We consider a state feedback controller to stabilize the pendulum about the vertical equilibrium. The system has four states, $\mathbf{x}=[\alpha ,\beta ,\omega ,\varphi ]$: the angular position of the rotary arm and the pendulum, $\alpha ,\beta $, with $\beta =0$ being the vertical position, and the corresponding angular velocities, $\omega ,\varphi $. We control the voltage applied to the motor, ${V}_{\mathrm{m}}$. We use the encoder readings as estimates of the angular positions, $\widehat{\alpha},\widehat{\beta}$. We apply a lowpass filter to the difference of consecutive angular positions to estimate the angular velocities, $\widehat{\omega},\widehat{\varphi}$. We aim to find a controller of the form ${u}_{t}={V}_{\mathrm{m}}=[{\theta}_{1},{\theta}_{2},{\theta}_{3},{\theta}_{4}]{[\widehat{\alpha},\widehat{\beta},\widehat{\omega},\widehat{\varphi}]}^{\top}$.
Scaling and reward. We define a state and action dependent reward as the negative of a quadratic cost, $r(\mathbf{x},\mathbf{u})=({\mathbf{x}}^{\top}Q\mathbf{x}+{\mathbf{u}}^{\top}R\mathbf{u})$ with $Q=\text{diag}(1,10,0,0)$ and $R=8$. The performance associated to a controller is the expected average reward it induces on the system, that is, for a trajectory of duration $T$, $\mathbb{E}[1/T{\int}_{0}^{T}r({\mathbf{x}}_{t},\pi ({\mathbf{x}}_{t}\theta ))\mathit{d}t]$. To prevent one of the objectives from dominating the contribution to the hypervolume improvement in the EHI algorithm, we must normalize them. We control the range of the robustness indicators, see creftype IIIB, and, therefore it is easy to rescale them to the $[0,1]$ range. We observe empirically that the unnormalized return ranges in $[500,0]$. Thus, we clip every return value to this range and we rescale it to the $[0,1]$ interval. Since the pendulum incurs substantially different returns when a stabilizing or destabilizing controller is used, we cannot rescale the range linearly. Instead, we use a piecewise linear function. In particular, since we observe empirically that stabilizing controllers have a performance between 20 and 0, we rescale linearly the range $[500,20]$ to $[0,0.5]$ and the range $[20,0]$ to $[0.5,1]$. This differentiates coarsely the quality of unstable controllers, and it gives a more refined scale over stable ones.
Standard  Motor noise  Sensor noise  Add $2\text{g}$  

$\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  
Scalar BO  0.150  80%  0.92  0.151  80%  0.97  0.151  80%  1.05  0.185  100%  1.03 
MOBODM  0.044  32%  4.61  0.038  20%  4.07  0.063  20%  4.32  0.126  84%  3.21 
MOBOGM  0.003  0%  $\mathbf{\infty}$  0.013  0%  $\mathbf{\infty}$  0.057  0%  $\mathbf{\infty}$  0.004  0%  $\mathbf{\infty}$ 
Add $5\text{g}$  Add $9\text{g}$  Add $10\text{g}$ and $8\text{cm}$  Add $5\text{g}$, motor + sensor noise  

$\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  $\mathbb{E}\left[R\right]$  Fail  Fail time (s)  
Scalar BO  0.101  0%  $\mathbf{\infty}$  0.669  100%  4.07  0.711  100%  3.59  0.259  0%  $\mathbf{\infty}$ 
MOBOGM  0.031  0%  $\mathbf{\infty}$  0.026  0%  $\mathbf{\infty}$  0.0259  0%  $\mathbf{\infty}$  0.366  0%  $\mathbf{\infty}$ 
Surrogate models. For the nonrobust algorithm, we use a standard GP model with a zero mean prior and a Matérn kernel with $\nu =5/2$ with automatic relevance determination (ARD). We set the hyperprior over the lengthscales to $\mathrm{Lognormal}(1,3)$ and over the standard deviation to $\mathrm{Lognormal}(0.35,1)$. We use a Gaussian likelihood with no hyperprior. Similarly, for the robust algorithm, we use a zero prior mean. The correlation in the input space in the ICM model is captured by an ARD Matérn kernel with $\nu =5/2$, with the same hyperprior as the nonrobust case. For the correlation in the output space, we set a Gaussian hyperprior over each entry of the matrix $\mathbf{B}$, $\mathcal{N}(0,1)$. We use a Gaussian likelihood with a diagonal covariance matrix. In both cases, we udpate the hyperparameters using a maximum a posteriori estimate after every new data point is acquired.
Training. In the simtoreal setting, we train 5 different controllers for each of these methods: scalar BO (nonrobust), MOBO with performance and delay margin (DM), and MOBO with performance and lower gain margin (GM). The training consists of 200 BO iterations evaluated in simulation. In the hardware training setting, we train one controller for scalar BO and one for MOBOGM using 70 BO iterations evaluated on hardware. In both settings, MOBO requires fewer iterations than the given budget to find satisfactory solutions. Thus, using a stopping criterion in creftype IIIB would reduce the total number of iterations. We estimate performance by averaging the return over 10 independent runs. To estimate robustness, we require that the controller stabilizes the system for a given delay or gain for 5 independent runs. A trial is deemed stable if ${\alpha}_{t}\in [{8}^{\circ},{8}^{\circ}]$ and ${\beta}_{t}\in [{4}^{\circ},{4}^{\circ}]$ for all $t\in [4,5]$. Every training run lasts for 5 seconds. creftypeplural 3\crefpairconjunction2 show the fronts obtained by the MOBODM and MOBOGM simtoreal training, respectively. The gray circles correspond to controllers that appeared stabilizing at first, but that were ruled out with longer simulations, cf. creftype IIIB. The green squares indicate controllers tested on hardware. To emphasize the generality of our method, they were selected to be approximately at the elbow of the front without further tuning.
Simtoreal test. We test each controller learned in simulation on the hardware 5 times in 4 scenarios: (i) standard simtoreal, (ii) simtoreal adding Gaussian noise $\mathcal{N}(0,0.5)$ to the motor voltage, (iii) simtoreal adding noise to the encoder readings following a multinomial distribution over the integers in $[4,4]$ with $p(0)=0.6$ and 0.05 everywhere else, and (iv) simtoreal with the pendulum mass increased by $2\mathrm{g}$. A run is a failure if $\beta >{20}^{\circ}$. In creftype I, we compare the controllers in terms of average return, failure rate, and failing time, averaged over the runs that resulted in a failure. The robust methods consistently outperform the nonrobust policy optimization across all test scenarios. It appears that the lower gain margin is a more suitable robustness indicator in this setting. This may be due to the fact that, in our experience, the gain margin is less noisy to estimate.
Hardware test. We test each controller learned in hardware 5 times in 4 scenarios: (i) extra mass of $5\mathrm{g}$, (ii) extra mass of $9\mathrm{g}$, (iii) extra mass of $10\mathrm{g}$ and extra pendulum length of $8\text{cm}$, and (iv) extra mass of $5\mathrm{g}$ with the actuation and sensor noise used in the simtoreal experiments. creftype II summarizes the test results. Similarly to the simtoreal setting, the robust algorithm consistenlty outperforms its nonrobust counterpart.
V Concluding Remarks
We present a dataefficient algorithm for robust policy optimization based on multiobjective Bayesian optimization. We suggest a datadriven evaluation of two common robustness indicators, which is suitable to modelfree settings. Our hardware experiments on a Furuta pendulum show that (i) our method facilitates simulation to real transfer, and (ii) consistently increases robustness of the learned controllers as compared to BO with a single performance objective. Our results indicate a promising avenue toward robust learning control by leveraging robustness measures from control theory and multiobjective Bayesian optimization and point to several directions for extensions. While we show that gain and delay margings are effective in practice on a mildly nonlinear system, they may not fully characterize robust stability in general [31, 4]. Thus, investigating other relevant robustness indicators that can efficiently be estimated from data in a modelfree setting is a topic for future research. Also, using multiple robustness indicator simultaneously is relevant, which our method could do at the expense of a more complex scaling to balance robustness and performance.
References
 [1] R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction. MIT press, 2018.
 [2] B. D. Anderson and J. B. Moore, Optimal control: linear quadratic methods. Courier Corporation, 2007.
 [3] J. Doyle, “Guaranteed margins for LQG regulators,” IEEE Transactions on Automatic Control, vol. 23, no. 4, pp. 756–757, 1978.
 [4] B. Boulet and Y. Duan, “The fundamental tradeoff between performance and robustness,” IEEE Control Systems Magazine, vol. 27, no. 3, pp. 30–44, June 2007.
 [5] K. Zhou and J. C. Doyle, Essentials of robust control. Prentice Hall, 1998.
 [6] A. Saberi, S. Peddapullaiah, and B. M. Chen, Loop transfer recovery: analysis and design. SpringerVerlag, 1993.
 [7] M. Green and D. J. Limebeer, Linear robust control. Courier Corporation, 2012.
 [8] A. Nilim and L. El Ghaoui, “Robust control of markov decision processes with uncertain transition matrices,” Operations Research, vol. 53, no. 5, pp. 780–798, 2005.
 [9] R. Sharma and M. Gopal, “A robust markov game controller for nonlinear systems,” Applied Soft Computing, vol. 7, no. 3, pp. 818–827, 2007.
 [10] E. Delage and S. Mannor, “Percentile optimization for markov decision processes with parameter uncertainty,” Operations research, vol. 58, no. 1, pp. 203–213, 2010.
 [11] K. J. Astrom and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers. Princeton, NJ, USA: Princeton University Press, 2008.
 [12] M. NeumannBrosig, A. Marco, D. Schwarzmann, and S. Trimpe, “Dataefficient autotuning with bayesian optimization: An industrial control study,” IEEE Transactions on Control Systems Technology, 2019.
 [13] H. K. Venkataraman and P. J. Seiler, “Recovering robustness in modelfree reinforcement learning,” in 2019 American Control Conference (ACC). IEEE, 2019, pp. 4210–4216.
 [14] L. Pinto, J. Davidson, R. Sukthankar, and A. Gupta, “Robust adversarial reinforcement learning,” in International Conference on Machine Learning, 2017, pp. 2817–2826.
 [15] J. Morimoto and K. Doya, “Robust reinforcement learning,” in Advances in Neural Information Processing Systems, 2001, pp. 1061–1067.
 [16] A. Gambier and M. Jipp, “Multiobjective optimal control: An introduction,” in Control Conference (ASCC), 2011 8th Asian. IEEE, 2011, pp. 1084–1089.
 [17] J. Mockus, V. Tiesis, and A. Zilinskas, “The application of bayesian methods for seeking the extremum,” Towards global optimization, vol. 2, no. 117129, p. 2, 1978.
 [18] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas, “Taking the human out of the loop: A review of Bayesian optimization,” Proceedings of the IEEE, vol. 104, no. 1, pp. 148–175, 2016.
 [19] F. Berkenkamp, A. P. Schoellig, and A. Krause, “Safe controller optimization for quadrotors with gaussian processes,” in IEEE International Conference on Robotics and Automation (ICRA), May 2016, pp. 491–496.
 [20] A. Marco, F. Berkenkamp, P. Hennig, A. P. Schoellig, A. Krause, S. Schaal, and S. Trimpe, “Virtual vs. real: Trading off simulations and physical experiments in reinforcement learning with bayesian optimization,” in 2017 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2017, pp. 1557–1563.
 [21] R. Calandra, A. Seyfarth, J. Peters, and M. P. Deisenroth, “Bayesian optimization for learning gaits under uncertainty,” Annals of Mathematics and Artificial Intelligence, vol. 76, no. 1, pp. 5–23, Feb 2016.
 [22] R. Antonova, A. Rai, and C. G. Atkeson, “Deep kernels for optimizing locomotion controllers,” in Proceedings of the 1st Annual Conference on Robot Learning, ser. Proceedings of Machine Learning Research, vol. 78, 2017, pp. 47–56.
 [23] M. Emmerich and J.w. Klinkenberg, “The computation of the expected improvement in dominated hypervolume of pareto front approximations,” Rapport technique, Leiden University, vol. 34, 2008.
 [24] M. Zuluaga, A. Krause, and M. Püschel, “ePAL: An Active Learning Approach to the MultiObjective Optimization Problem,” Journal of Machine Learning Research, vol. 17, no. 104, pp. 1–32, 2016.
 [25] D. HernándezLobato, J. HernandezLobato, A. Shah, and R. Adams, “Predictive entropy search for multiobjective bayesian optimization,” in International Conference on Machine Learning, 2016, pp. 1492–1501.
 [26] M. Tesch, J. G. Schneider, and H. Choset, “Expensive multiobjective optimization for robotics,” 2013 IEEE International Conference on Robotics and Automation, pp. 973–980, 2013.
 [27] R. Ariizumi, M. Tesch, H. Choset, and F. Matsuno, “Expensive multiobjective optimization for robotics with consideration of heteroscedastic noise,” in 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. IEEE, sep 2014, pp. 2230–2235.
 [28] Y. Collette and P. Siarry, Multiobjective optimization: Principles and case studies. Springer Science & Business Media, 2013.
 [29] C. E. Rasmussen, “Gaussian processes in machine learning,” in Advanced lectures on machine learning. Springer, 2004, pp. 63–71.
 [30] M. A. Alvarez, L. Rosasco, N. D. Lawrence et al., “Kernels for vectorvalued functions: A review,” Foundations and Trends in Machine Learning, vol. 4, no. 3, pp. 195–266, 2012.
 [31] K. Zhou, J. C. Doyle, K. Glover et al., Robust and optimal control. Prentice hall New Jersey, 1996, vol. 40.
 [32] 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). IEEE, May 2016, pp. 270–277.
 [33] B. S. Cazzolato and Z. Prime, “On the dynamics of the furuta pendulum,” Journal of Control Science and Engineering, vol. 2011, p. 3, 2011.
 [34] Quanser, Qube Servo 2  Student workbook, Quanser Consulting Inc.