Learning-based Model Predictive Control for Safe Exploration and Reinforcement Learning
Torsten Koller Felix Berkenkamp Matteo Turchetta Joschka Boedecker Andreas Krause










































































Learning-based Model Predictive Control for Safe Exploration and Reinforcement Learning
    Torsten Koller^*,
        Felix Berkenkamp^*,
        Matteo Turchetta,
        Joschka Boedecker
        and Andreas Krause
^* T. Koller and F. Berkenkamp contributed equally to this work.
T. Koller and J. Boedecker are with the Neurorobotics Group, University of Freiburg, Freiburg, Germany (email: {kollert, jboedeck}@informatik.uni-freiburg.de}).
F. Berkenkamp, M. Turchetta and A. Krause are with Learning and Adaptive Systems Group, ETH Zurich, Zurich, Switzerland (email: {befelix, matteotu, krausea}@inf.ethz.ch}).
This work was supported by SNSF grant 200020_159557, a fellowship within the FITweltweit program of the German Academic Exchange Service (DAAD), an Open Philantropy Project AI fellowship, the Max Planck ETH Center for Learning Systems, and the BrainLinks-BrainTools  Cluster  of  Excellence  (grant  number  EXC1086). T. Koller is funded through the State Graduate Funding Program of Baden-Wuerttemberg


     
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






Reinforcement learning has been successfully used to solve difficult tasks in complex unknown environments. However, these methods typically do not provide any safety guarantees during the learning process. This is particularly problematic, since reinforcement learning agents actively explore their environment. This prevents their use in safety-critical, real-world applications.

In this paper, we present a learning-based model predictive control scheme that provides high-probability safety guarantees throughout the learning process.

Based on a reliable statistical model, we construct provably accurate confidence intervals on predicted trajectories. Unlike previous approaches, we allow for input-dependent uncertainties.

Based on these reliable predictions, we guarantee that trajectories satisfy safety constraints. Moreover, we use a terminal set constraint to recursively guarantee the existence of safe control actions at every iteration.

We evaluate the resulting algorithm to safely explore the dynamics of an inverted pendulum and to solve a reinforcement learning task on a cart-pole system with safety constraints.



























§ INTRODUCTION


In model-based reinforcement learning (RL,<cit.>), we aim to learn the dynamics of an unknown system from data and, based on the model, derive a policy that optimizes the long-term behavior of the system. In order to be successful, these methods need to find an efficient trade-off between exploitation, i.e. finding and executing optimal policies based on the current knowledge of the system and exploration, where the system traverses regions of the state-space that are unknown or uncertain to collect observations that improve the model estimation. 
In real-world safety-critical systems, we further need to restrict exploration to actions that are safe to perform by satisfying state and control constraints at all times.
In contrast, current approaches often use exploration schemes that lead to unpredictable and possibly unsafe behavior of the system, e.g., based on the optimism in the face of uncertainty principle <cit.>. Consequently, such approaches have limited applicability to real-world systems. 

In this paper, we present a safe learning-based model predictive control (MPC) scheme for nonlinear systems with state-dependent uncertainty. The method guarantees the existence of feasible return trajectories to a safe region of the state space at every time step with high probability. Specifically, we estimate the system uncertainty based on observations through statistical modelling techniques. Based on the statistical model, we derive two reliable uncertainty propagation techniques that allow for state-dependent uncertainty. We show that the resulting MPC scheme provides high probability safety guarantees and combine it with exploration schemes from RL. We show that our algorithm can safely learn about the system dynamics (system identification) when it is combined with an exploration scheme that aims to reduce the uncertainty about the model through safe exploration. Moreover, we apply our algorithm to a reinforcement learning task where the algorithm only aims to collect task-relevant data in order to maximize expected performance on a RL task. In our experiments, we empirically show the importance of both, planning a safety trajectory in avoiding crashes and planning a performance trajectory in solving the task and how the design choices of our algorithms can improve the performance in a RL task. 

Preliminary results on one of the uncertainty propagation scheme for safe control, but without RL, were shown in <cit.>. This paper provides additional information, proposes an additional uncertainty propagation scheme, shows how the algorithm can be combined with RL, and provides a more thorough experimental evaluation.







  
*Outline
We provide a brief introduction to our work and outline our contributions in <ref>, relating them to previous work in <ref>. The problem statement is given in <ref> and the necessary background for our research is provided in <ref>. In <ref>, we derive our multi-step ahead prediction technique. In <ref>, we then use this technique to derive a MPC scheme and prove its safety. <ref> discusses how we can efficiently use this for exploration and to solve a safe model-based RL task. In <ref>, we briefly address practical application issues and in <ref>, we report our experimental results in safe exploration and safe RL on an inverted pendulum and a cart-pole system, respectively. We conclude our paper with a discussion in <ref>.




§ RELATED WORK


One area that considers safety guarantees in partially unknown systems is robust MPC. These methods iteratively optimize the performance along finite-length trajectories at each time step, based on a known model that incorporates uncertainties and disturbances acting on the system <cit.>. In a constrained robust MPC setting, we optimize these local trajectories under additional state and control constraints. Safety is typically defined in terms of recursive feasibility and robust constraint satisfaction.

In <cit.>, this definition is used to safely control urban traffic flow, while <cit.> guarantees safety by switching between a standard and a safety mode. However, these methods are conservative, since they do not update the model. 


In contrast, learning-based MPC approaches (LBMPC) adapt their models online based on observations of the system <cit.>. This allows the controller to improve over time, given limited prior knowledge of the system. Theoretical safety guarantees in LBMPC are established in <cit.>. Both approaches compute tubes around trajectories of a known nominal linear model, which accounts for all disturbances. While the former approach requires disturbances to be bounded in an a-priori specified  polytope, the latter relies on sampling to approximate model errors. Robust constraint satisfaction and closed-loop stability for LBMPC with nominal linear models and state-dependent uncertainty is shown in <cit.>. However, the method relies on extensive offline computations. In <cit.>, instead of updating a model of the system, samples are used to grow a safe set in a robust MPC setting over time.   

Another area that considers learning for control is model-based RL, where we learn a model from observations of our system and use it to find a policy for a given task based on an external cost signal <cit.>. One statistical model that is commonly used in model-based RL are Gaussian processes (GP, <cit.>). In RL, the GP uncertainty estimates are typically used to optimize the expected long-term cost of a feedback policy by propagating the state-dependent uncertainties in the forward predictions of the GP model over multiple time steps. However, trajectories have complex dependencies on states and unbounded stochastic uncertainties, which requires approximation techniques to get tractable solutions. The expected cost under the forward propagated approximate state distributions are then optimized w.r.t. parameterized or MPC-based policies <cit.>.

In contrast to MPC, where we optimize finite-length trajectories based on a pre-specified system model, in RL we typically aim to find an infinite horizon optimal policy under unknown dynamics of our system. Hence, enforcing hard constraints in RL is challenging. In <cit.>, a safe exploration heuristic for RL is introduced that uses the proximity to a finite set of safe states generated from pre-defined exploration policies as the safety criterion.
In <cit.>, a safety layer is introduced that acts on top of an RL agent's possibly unsafe actions to prevent immediate negative consequences, but does not deal with negative long-term consequences of an action. 
Safety based on constrained Markov decision processes <cit.> is considered in  <cit.>, where the RL agent additionally receives a constraint cost signal and the task is to keep the expected sum of constraint costs below a given threshold. The issue with approaches based on this notion of expected safety is that we can still have arbitrary large and frequent violations of our constraints as we do not account for higher moments of the accumulated constraint cost distribution. Control-theoretic safety properties such as Lyapunov stability or robust constraint satisfaction have received little attention <cit.>. In <cit.>, safety is guaranteed by optimizing parametric policies under stability constraints, while <cit.> guarantees safety in terms of constraint satisfaction through reachability analysis. In <cit.>, a safety filter is introduced that employs robust MPC to correct for possibly unsafe actions of an RL agent. While the authors show safety in terms of constraint satisfaction, the RL agent and the safety mechanism are fundamentally detached from each other. This could potentially lead to suboptimal and possibly oscillatory behavior, where the RL agent constantly tries to approach the constraints.
In <cit.>, methods from reinforcement learning are used to approximately solve the constrained linear quadratic regulator (LQR), but the system is assumed to be known and the approach is limited to linear systems. In <cit.>, a learning-based technique for the constrained LQR for unknown linear dynamics with sub-optimality bounds and statistical learning guarantees is proposed. 
In GP-based RL, safety is often considered through probabilistic chance constraints <cit.> based on approximate uncertainty propagation. While often empirically successful, these approaches do not theoretically guarantee safety of the underlying system, since constraints are typically only enforced over finite horizons and without incorporating the approximation errors of the multi-step ahead predictions.















§ PROBLEM STATEMENT


We consider a deterministic, discrete-time  dynamical system

    x_t+1 = f(x_t,u_t) = h(x_t,u_t)_prior model + g(x_t,u_t)_unknown error,

where x_t ∈^p is the state and u_t ∈^q is the control input to the system at time step t∈. We assume that we have access to a known, Lipschitz continuous prior model  h(x_t,u_t) and that the a priori unknown model-error g(x_t,u_t) is Lipschitz continuous as well. The prior model h could be based on an idealized first principles physics model that may be inaccurately identified. By collecting transition data, we want to reduce these model uncertainties and learn the unknown model-error g from the observations of our system. In order to learn about g, we use a statistical model with mean and variance given by (μ_n(x_t,u_t),Σ_n(x_t,u_t)). This provides us with a point-wise approximation of g, given by μ_n, and corresponding input-dependent uncertainty estimates Σ_n(x_t,u_t) = ([σ_n,1^2(x_t,u_t),..,σ_n,p^2(x_t,u_t)]), where the subscript n ∈ denotes the number of observations we have collected from our system so far. As we collect more observations from the system during operation, we want to learn and to improve this estimate over time. 

We assume that the system is subject to polytopic state and control constraints

    = {x ∈^p | H_x x ≤ h_x,   h_x ∈^m_x}, 
     = {u ∈^q | H_u u ≤ h_u,   h_u ∈^m_u},

which are given by the intersection of m_x, m_u ∈ individual half-spaces, respectively. As a running example, we consider an autonomous driving scenario, where the state region could correspond to a highway lane or race track and the control constraints could represent the physical limits on acceleration and steering angle of the car.

In order to provide safety guarantees, we need reliable estimates of the model-error inside the operating region 𝒳×𝒰 of our system. Restricted to this region, we make the assumption that, with high probability jointly throughout the operating time of our system, at a finite number of inputs in each time step, values of our unknown but deterministic function  g are contained in the confidence intervals of our statistical model.



 For all n ∈ and a fixed sample size T ∈, pick arbitrary sets of inputs 𝒟_n = {z_1,..,z_T}⊂𝒳×𝒰. For any confidence level δ∈ (0,1], we have access to scaling factors β_n,T > 0 such that uniformly with probability greater than 1-δ, we have that ∀  n ∈,  1≤ j ≤ p, 1≤ k ≤ T: 

    |μ_n,j(z_k) - g_j(z_k)| ≤β_n,T·σ_n,j(z_k).



Intuitively, this means that we can trust our model to return well-calibrated estimates. For notational convenience, we drop the subscripts of β_n,T and refer to it as β in the remainder of the paper.
In general, providing these kind of guarantees is impossible for arbitrary functions g. We provide further analysis in <ref> on how we can guarantee this sort of reliability when using a Gaussian process (GP) statistical model to approximate a certain class of deterministic functions. 

Lastly, we assume access to a backup controller that guarantees that we never violate our safety constraints once we enter a safe subset of the state space. In the autonomous driving example, this could be a simple linear controller that stabilizes the car in a small region in the center of the lane.


We are given a controller π_safe(·) and a polytopic safe region

    _safe := { x∈^p | H_s x ≤ h_s, h_s ∈^m_s}⊆,

at the intersection of m_s ∈ individual half-spaces. We denote with
x_t+1 = f_π_safe = f(x_t,π_safe(x_t)) the closed-loop system under π_safe and assume for some arbitrary k ∈ that

    x_k ∈𝒳_safe⇒  f_π_safe(x_t) ∈𝒳,   ∀ t ≥ k.



This assumption allows us to gather initial data from the system even in the presence of significant model errors, since the system remains safe under the controller π_safe. Moreover, we can still guarantee constraint satisfaction asymptotically outside of _safe, if we can show that a finite sequence of control inputs eventually steers the system back to the safe set. In <cit.>, we stated a stricter assumption that required _safe to be robust control positive invariant (RCPI) under π_safe, i.e. the system is not allowed to leave the safe set once it enters it. Our relaxed definition only requires the safe set to be a subset of a RCPI set that is contained in the set of admissible states <cit.>.

Given a controller π, ideally we want to enforce the state- and control constraints at every time step,

    ∀ t ∈: f_π(x_t) ∈𝒳,  π(x_t) ∈𝒰.

Apart from π_safe, which trivially and conservatively fulfills this if x_0 ∈_safe, it is in general impossible to design a controller that enforces
(<ref>) without additional assumptions. This is due to the fact that we only have access to noisy observations of our system. Instead, we slightly relax this requirement to safety with high probability throughout its operation time.


Let π: ^p→^q be a controller for (<ref>) with the corresponding closed-loop system f_π.
Let x_0 ∈𝒳_safe and δ∈ (0,1]. A system is δ-safe under the controller π iff:

    [  ∀ t ∈:  f_π(x_t) ∈𝒳,  π(x_t) ∈𝒰] ≥ 1-δ.


As a comparison, in chance constrained MPC we typically get per time step high probability constraint satisfaction. Informally speaking, there we have definitions of the form ∀ t ∈: [No collision at time t] ≥ 1-δ. In contrast, <ref> requires high probability safety guarantees that are independent from the duration of system operation.

Aside from these mandatory safety requirements, we are provided with a mission objective and we want our controller to improve its performance in solving a given task as we learn about our system trough observations. In this paper, we are particularly interested in two intimately related tasks. Firstly, we are interested in safely and efficiently identifying the system by exploring the safe region of the state space.
Secondly, we want to find a controller that solves a RL task by minimizing the accumulated sum of costs our system receives in every state while still remaining safe at all times. This naturally entails the former task, namely,  a mechanism that can safely explore the state space while trying to solve the task. The next section elaborates on this.




§ BACKGROUND


In this section, we provide a brief introduction to model-based RL and introduce the necessary background on GPs and set-theoretic properties of ellipsoids that we use to model our system and perform multi-step ahead predictions.




 §.§ Model-based Reinforcement Learning

In reinforcement learning (RL), an agent interacts with an unknown environment according to a policy π: ^p →^q by applying an action u_t = π(x_t) ∈^q in a state x_t ∈^p. The environment defines a probability distribution over the next state X_t+1∼ p(x_t+1|x_t,u_t). After every time step the agent receives a scalar cost c(x_t,u_t) ∈. 
The performance of an agent is then measured through the expected cumulative cost

    V^π(x) = ∑_t=0^∞𝔼_x_t ∼ X_t[ γ^t c(x_t,π(x_t))| x_0 = x],   x ∈^p,

where the discount-factor γ∈ [0,1) guarantees convergence of the infinite summation. In our special case of a deterministic environment, we can remove the expectation over state trajectories. 
The goal of a reinforcement problem is to find 

    π^* = min_π  V^π(x),

the policy that minimizes the expected cost of interacting with its environment.

Approaches to solving these RL tasks can generally be subsumed under model-based and model-free RL algorithms. We are interested in the former, where we try to learn an approximate model of the environment from experience samples  (x_t,u_t,x_t+1) that allows us to directly approximate V^π(x) in equation (<ref>) for an agent π. This is typically done by replacing the infinite sum in (<ref>) with a finite horizon look-ahead and using the model to approximate the expectation over trajectories. Depending on the class of policies under consideration, various ways of optimizing (<ref>) are proposed in the literature, including MPC <cit.>.



 §.§ Gaussian Processes and RKHS

One way to learn the unknown model-error g from data is by using a GP model.
A 𝒢𝒫(m,k) is a distribution over functions, which is fully specified through a mean function m: ^d→ and a covariance function k: ^d ×^d →, where d = p + q. Given a set of n noisy observations[For notational simplicity, we only outline the equations for p=1, i.e. a single-output system.] f̂_i = f(z_i) + w_i,   w_i ∼𝒩(0,λ^2),  z_i = (x_i,u_i),  i = 1,…,n,  λ∈, we choose a zero-mean prior on g as m ≡ 0 and regard the differences y_n = [f̂_1 - h(z_1),…,f̂_n - h(z_n)]^T between prior model h and observed system response at input locations Z = [z_1,..,z_n]^T. The posterior distribution at z is then given as a Gaussian 𝒩(μ_n(z),σ_n^2(z))
with mean and variance

    μ_n(z)    = k_n(z)^T [K_n+λ^2 I_n]^-1y_n, 
    σ^2_n(z)    = k(z,z) - k_n(z)^T [K_n + λ^2 I_n]^-1k_n(z),

where [K_n]_ij = k(z_i,z_j), [k_n(z)]_j = k(z,z_j), and I_n is the n-dimensional identity matrix.
In the case of multiple outputs p > 1, we model each output dimension with an independent GP, 𝒢𝒫(m_j,k_j), j = 1,..,p. We then redefine (<ref>) and (<ref>) as μ_n(·) = (μ_n,1(·),..,μ_n,n_x(·)) and σ_n(·) = (σ_n,1(·),..,σ_n,n_x(·)) corresponding to the predictive mean and variance functions of the individual models.



  §.§.§ RKHS

We consider the space of functions ℋ_k that is defined through the completion of all possible posterior mean functions (<ref>) of a 𝒢𝒫(0,k), i.e.

    ℋ_k = {∑_l=0^∞α_l k(z,z_l) |  α_l ∈, z,z_l ∈^d}.

This class of functions is well-behaved in the sense that they form a reproducing kernel Hilbert space (RKHS, <cit.>) equipped with an inner-product ⟨·,·⟩_k. The induced norm ||·||_k^2 is a measure of the complexity of a function in ℋ_k.
We consider the case, where our model-error g is indeed of the form g(z) = ∑_l=0^∞α_l k(z,z_l),  α_l ∈, z,z_l ∈^d, i.e. g ∈ℋ_k. Using results from the statistical learning theory literature, we can directly derive useful properties from this assumption.


<cit.>:
Let g have bounded RKHS norm ||g||_k ≤ B_g induced by a continuously differentiable kernel k. Then g is Lipschitz continuous.

That is, the regularity assumption on our model-error and the smoothness assumption on the covariance function k directly imply that the function g is Lipschitz.

More importantly, the assumption g ∈ℋ_k guarantees that using GPs to model the unknown part of the system (<ref>), provides us with reliable confidence intervals on the model-error g as required in <ref>.


<cit.>: Assume ||g||_k ≤ B_g and that measurements are corrupted by λ-sub-Gaussian noise. Let β_n = B_g + 4λ√(γ_n + 1 + ln(1/δ)), where γ_n is the information capacity associated with the kernel k. Then with probability at least 1-δ, with δ∈ (0,1), we have for all 1≤ j ≤ n_x,  z ∈𝒳×𝒰 that |μ_n-1,j(z) - g_j(z)| ≤β_n ·σ_n-1,j(z).


In combination with the prior model h(z), this allows us to construct reliable confidence intervals around the true dynamics of the system <ref>. The scaling β_n depends on the number of data points n that we gather from the system through the information capacity,
γ_n = max_𝒜⊂𝒵̃ , |A| = ñ I(g̃_A;g),  𝒵̃ = 𝒳×𝒰×ℐ,  ñ = n · p,  ℐ = {1,..,p}, i.e. the maximum mutual information I(g̃_A,g) between a finite set of samples A and the function g.
Exact evaluation of γ_n is NP-hard in general, but it can be greedily over-approximated with bounded error and has sublinear dependence on n for many commonly used kernels <cit.>. <ref> is stronger than required by <ref>, since the confidence level δ holds jointly for all z ∈𝒳×𝒰 and not only for a finite set of T samples. As a side note, we refer to a result in <cit.> showing that <ref> holds if our system transitions x_t+1 is represented by stochastic draws from the GP posterior <ref>, with z = (x_t, u_t).[Due to the difficulty of designing a safety controller and safe region that would fulfill <ref> in such a stochastic system (see <cit.> for a discussion), we nonetheless consider a deterministic system.]



 §.§ Ellipsoids

We use ellipsoids to bound the uncertainty of our system when making multi-step ahead predictions.
Due to appealing geometric properties, ellipsoids are widely used in the robust control community to compute reachable sets <cit.>. These sets intuitively provide an outer approximation on the next state of a system considering all possible realizations of uncertainties when applying a controller to the system at a given set-valued input.
We briefly review some of these properties and refer to <cit.> for an exhaustive introduction to ellipsoids and to the derivations for the following properties.

We use the basic definition of an ellipsoid,

    ℰ(p,Q) := { x ∈^n | (x-p)^T Q^-1(x-p) ≤ 1},

with  center p ∈^n and a symmetric positive definite (s.p.d) shape matrix Q ∈^n × n.
Ellipsoids are invariant under affine subspace transformations such that, for A ∈^n × r, r ≤ n with full column rank and b ∈^r, we have that

    A ·ℰ(p,Q) + b = ℰ(p+b,AQA^T).

The Minkowski sum ℰ(p_1,Q_1) ⊕ℰ(p_2,Q_2), i.e. the pointwise sum between two arbitrary ellipsoids, is in general not an ellipsoid anymore, but we have

    ℰ(p_1,Q_1) ⊕ℰ(p_2,Q_2) ⊂ℰ_c(p̃,Q̃),

where p̃ = p_1 + p_2,  Q̃ = (1+c^-1)Q_1 + (1+c)Q_2 for all c > 0. Moreover, the minimizer of the trace of the resulting shape matrix is analytically given as c = √(Tr(Q_1) / Tr(Q_2)). A particular problem that we encounter is finding the maximum distance r to the center of an ellipsoid ℰ(0,Q) under a special transformation, i.e.

    r(Q,S) = max_x ∈ℰ(p,Q) ||S(x-p)||_2 =  max_s^T Q^-1s ≤ 1 s^T S^T S s,

where S ∈^m × n with full column rank. This is a generalized eigenvalue problem of the pair (Q,S^TS) and the optimizer is given as the square-root of the largest generalized eigenvalue.











§ ROBUST MULTI-STEP AHEAD PREDICTIONS



In order to plan safe trajectories based on our statistical model, we need to reliably estimate the region of the state space that can be reached over multiple time steps under a sequence of control inputs. Based on <ref> and our prior model h(x_t,u_t), we directly obtain high-probability confidence intervals on f(x_t,u_t) uniformly for all t ∈ given a single control input u_t. We extend this to over-approximate the system after a sequence of inputs (u_t,u_t+1,..,u_T-1). The result is a sequence of set-valued confidence regions that contain the true trajectory of the system with high probability.



 §.§ One-step ahead predictions


We derive a function that takes as an input an ellipsoidal subset of the state space and outputs a second ellipsoid providing a confidence region around the next state of the system under a specific control input, given that the current state is contained in the input ellipsoid.


    

In order to approximate the system, first we linearize our prior model h(x_t,u_t) and use the affine transformation property (<ref>) to compute the ellipsoidal next state of the linearized prior model. Next, we approximate the unknown model-error g(x_t,u_t) using the confidence intervals of our statistical model. We propose a locally constant approximation of g in <ref>, while in <ref>, we locally linearize our statistical model.  We finally apply Lipschitz arguments to outer-bound the approximation errors. We sum up these individual approximations, which result in an ellipsoidal approximation of the next state of the system. This is illustrated in <ref>. We formally derive the necessary equations in the following paragraphs. The reader may choose to skip the technical details of these approximations, which result in  <ref>. To guide the interested reader through the derivations, we provide a list of important variables functions and constants in <ref>.





  §.§.§ Predictions via locally constant model approximation


We first regard the system f in (<ref>)
for a single input vector z = (x,u),  f(z) = h(z) + g(z).
We linearly approximate f around z̅ = (x̅,u̅) via

    f(z) ≈ h(z̅) + J_h(z̅) (z-z̅) + g(z̅) = f̃(z),

where J_h(z̅) = [A_h,B_h] is the Jacobian of h at z̅.

Next, we use the Lagrangian remainder theorem <cit.> on the linearization of h and apply a continuity argument on our locally constant approximation of g. This results in an upper-bound on the approximation error,


    |f_j(z) - f̃_j(z) | ≤L_∇ h,j/2 ||z - z̅||_2^2 + L_g ||z-z̅||_2,


where f_j(z) is the jth component of f, 1≤ j ≤ p, L_∇ h,j is the Lipschitz constant of the gradient ∇ h_j, and L_g is the Lipschitz constant of g, which exists by <ref>.

The function f̃ depends on the unknown model error g. We approximate g with the statistical model, μ(z̅) ≈ g(z̅). From <ref> we have

    |g_j(z̅)-μ_j(z̅)| ≤βσ_j(z̅),   1≤ j ≤ p,

with high probability.
We combine (<ref>) and (<ref>) to obtain

    | f_j(z) - f̃_μ,j(z)| ≤βσ(z̅) + L_∇ h,j/2 ||z - z̅||_2^2 + L_g ||z-z̅||_2,

where 1 ≤ j ≤ p and

f̃_μ(z) = h(z̅) + J_h(z̅) (z-z̅) + μ_n(z̅).
We can interpret  (<ref>) as the edges of the confidence hyper-rectangle

    m̃(z)  = f̃_μ(z) ±[ βσ(z̅)+L_∇ h/2 ||z - z̅||_2^2 + L_g ||z-z̅||_2 ],

where L_∇ h = [L_∇ h,1,..,L_∇ h,p] and we use the shorthand notation a ± b := [a_1 ± b_1] × [a_p± b_p],  a,b ∈^p.

We are now ready to compute a confidence region based on an ellipsoidal state ℛ = ℰ(p,Q) ⊂^p and a fixed input u ∈^q, by over-approximating the output of the system f(ℛ,u) = { f(x,u)| x ∈ R} for all inputs contained in an ellipsoid ℛ. Here, we choose p as the linearization center of the state and choose u̅ = u, i.e. z̅ = (p,u). Since the function f̃_μ is affine, we can make use of (<ref>) to compute

    f̃_μ(ℛ,u) = ℰ(h(z̅) + μ(z̅),A_hQA_h^T),

which results in an ellipsoid.
This is visualized in <ref> by the upper ellipsoid in the center.
To upper-bound the confidence hyper-rectangle on the right hand side of (<ref>), we upper-bound the term z - z̅_2 by










    l(ℛ,u) := max_z(x) = (x,u), 
     x ∈ℛ ||z(x) - z̅||_2,

which leads to

    d̃(ℛ, u)  := βσ(z̅) + L_∇ h l^2(ℛ,u) / 2 + L_g l(R,u).


Due to our choice of z, z̅, we have that ||z(x) - z̅||_2 = ||x-p||_2 and  we can use (<ref>) to get
l(ℛ,u) = r(Q,I_p), which corresponds to the largest eigenvalue of Q^-1.
Using (<ref>), we can now over-approximate the right side of (<ref>) for inputs R by an ellipsoid

    0 ±d̃(ℛ,u) ⊂ℰ(0,Q_d̃(ℛ,u)),


where we obtain Q_d̃(ℛ,u)  by over-approximating the hyper-rectangle d̃(ℛ, u) with the ellipsoid ℰ(0,Q_d̃(ℛ,u)) through a ± b ⊂ℰ(a,√(p)·diag([b_1,..,b_p])),  ∀ a,b ∈^p. This is illustrated in <ref> by the lower ellipsoid in the center.
Combining the previous results, we can compute the final over-approximation using (<ref>),

    ℛ_+ = m̃(ℛ,u) = f̃_μ(R,u) ⊕ℰ(0,Q_d̃(ℛ,u)).

Since we carefully incorporated all approximation errors and extended the confidence intervals around our model predictions to set-valued inputs, we get the following generalization of the reliability <ref> for single inputs. 



Let δ∈ (0,1], choose β and T as in <ref> and D_n = {(x_1,n,u_1,n),..,(x_T,n,u_T,n)}⊂𝒳×𝒰 for each n ∈. Then, with probability greater than 1-δ, we have that:

    x_k,n∈ℛ_k,n⇒ f(x_k,n,u_k,n) ∈m̃(ℛ_k,n,u_k,n),

uniformly for all n ∈,   1 ≤ k ≤ T,  R_k,n = ℰ(p_k,n,Q_k,n) ⊂𝒳.


Define m(x,u) = h(x,u)+μ_n(x,u) ±βσ_n(x,u). From <ref> we have ∀ ℛ_k,n⊂𝒳,  u_k,n∈𝒰 with x_k,n∈ℛ_k,n that, uniformly with high probability, f(x_k,n,u_k,n) ∈ m(x_k,n,u_k,n). Due to the over-approximations, we have m(x_k,n,u_k,n) ⊂m̃(ℛ_k,n,u_k,n). 


<ref> allows us to compute confidence ellipsoid around the next state of the system, given that the current state of the system is known to be contained in an ellipsoidal confidence region.







  §.§.§ Predictions via model linearization


The previous derivations that lead to the reliable ellipsoidal one-step ahead predictions (<ref>) use Lipschitz arguments on the unknown model-error g to give an outer-bound bound on the linearization errors of our locally linear approximations. A different way to arrive at a reliable outer-bound, that can potentially reduce the conservatism of the proposed uncertainty propagation technique, is by using a first-order Taylor approximation of the statistical model μ_n.
Following previous derivations, we can approximate our unknown model-error with g(z) ≈μ_n(z̅) + A_μ_n (x-x̅) + B_μ_n (u - u̅) = P_μ_n^z̅(z), where z̅ = (x̅,u̅)^T and J_μ_n(z̅) = [A_μ_n,B_μ_n] is the Jacobian of μ_n at z̅.
Following (<ref>) - (<ref>), we can get an alternative outer bound by making the approximation

    f(z) ≈  h(z̅) + J_h(z̅) (z-z̅) + P_μ_n,j^z̅(z) = f̃_P_μ_n(z).

By adding and subtracting with g(z) and μ(z), the approximation error between f and f̃_P_μ can now be outer-bound via 

    | f_j(z) -  f̃_P_μ_n,j(z)|    ≤L_∇ h,j/2 ||z - z̅||_2^2 
       +  |P_μ_n,j^z̅(z) - μ_n,j(z)| + |μ_n,j(z) - g_j(z)|,

for z ∈𝒳×𝒰, 1 ≤ j ≤ p.
We use the reliability <ref> on our statistical model combined with Lipschitz arguments to obtain

    |P_μ_n,j^z̅(z) - μ_n,j(z)|    ≤L_∇μ_n,j/2 ||z - z̅||_2^2, 
    
         |μ_n,j(z) - g_j(z)|    ≤β(σ_n(z̅) + L_σ_n||z - z̅||_2),

where the second inequality holds uniformly
with high probability ∀ z ∈ D_n ⊂𝒳×𝒰, 1 ≤ j ≤ p, n ∈


We regard this over-approximation for ellipsoidal inputs ℛ=ℰ(p,Q). Closely following the derivations (<ref>),(<ref>), we can decompose our over-approximation into an affine part

    f̃_P_μ(ℛ,u) = ℰ(h(z̅) + μ(z̅),(A_h+A_μ)Q(A_h+A_μ)^T)

and combine the Lipschitz approximations into the hyper-rectangle

    0 ±β_n(σ(z̅) + L_σ l(R,u))+ L_∇ h+L_∇μ/2 l^2(R,u),

which results in the ellipsoidal over-approximation ℰ(0,Q_d̃_P(ℛ,u)). We arrive at the alternative one-step ahead prediction 

    m̃_P_μ(ℛ,u) = f̃_P_μ(ℛ,u) ⊕ℰ(0,Q_d̃_T(ℛ,u)).

Applying the arguments of <ref>, we can directly see that the function m̃_P_μ also provides reliable confidence ellipsoids for arbitrary ellipsoidal inputs R ⊂𝒳.

Apart from being possibly less conservative in many applications, we note that the one-step ahead predictions (<ref>) only require Lipschitz constants L_∇μ, L_∇ h, L_σ of the known functions ∇μ, ∇ h, σ. These could be identified using tools from the global optimization literature <cit.>.


 §.§ Multi-step ahead predictions

We use the previous results to compute a sequence of ellipsoids that contain a trajectory of the system with high-probability, by iteratively applying the one-step ahead predictions (<ref>). We note that the following line of arguments equivalently holds for the alternative one-step ahead prediction in (<ref>).

Given an initial ellipsoid ℛ_0 ⊂^p and control input u_t ∈^q, we iteratively compute confidence ellipsoids as

    ℛ_t+1 = m̃(ℛ_t,u_t).

We can directly apply <ref> to get the following result.


Let δ∈ (0,1] and choose β and T as in <ref>. For each n ∈ choose x_0,n∈ℛ_0,n⊂𝒳, {u_0,n,..,u_T-1,n}⊂𝒰. Then the following holds jointly for all n ∈, 0 ≤ t ≤ T-1 with probability at least 1-δ:

x_t,n∈ R_t,n,

where (x_t,n,u_t,n) ∈𝒳×𝒰, ℛ_1,n,..,ℛ_T,n⊂𝒳 is computed as in (<ref>) and {x_0,n,..,x_T,n} are trajectories of the system under input sequences {u_0,n,..,u_T-1,n} for each n ∈.
 
 
 Since <ref> holds for any set 𝒟_n, n ∈, we can choose 𝒟_n = {(x_0,n,u_0,n),..,(x_T-1,n,u_T-1,n)} and so that starting from tarting in x_0,n∈ℛ_0,n.  we get with high probability that x_i,n∈ℛ_i,n = m̃(ℛ_i-1,n,u_i-1,n) for all i=1,..T with n ∈.


<ref> guarantees that, with high probability, the system is always contained in the propagated ellipsoids (<ref>). However, this only holds if all confidence ellipsoids are inside the feasible state region 𝒳. Thus, if we provide safety guarantees for these sequences of ellipsoids, we obtain high-probability safety guarantees for the system (<ref>).


 §.§ Predictions under state-feedback control laws



When applying multi-step ahead predictions under a sequence of feed-forward inputs u_t ∈𝒰, the individual sets of the corresponding reachability sequence can quickly grow unreasonably large. This is because these open loop input sequences do not account for future control inputs that could correct deviations from the model predictions.
Hence, we extend (<ref>) to
affine state-feedback control laws of the form

    u_K,t(x_t) := K_t(x_t-p_t) + k_t,

where K_t ∈^q × p is a feedback matrix and k_t ∈^q is the open-loop input. The parameter p_t is determined through the center of the current ellipsoid ℛ_t = ℰ(p_t,Q_t). Given an appropriate choice of K_t, the control law actively contracts the ellipsoids towards their center. This technique is commonly used in tube-based MPC, to reduce the size of tubes around a nominal trajectory of the system that incorporate uncertainties and disturbances <cit.>.
Similar to the derivations (<ref>)-(<ref>), we can compute the function m̃ for affine feedback controllers (<ref>) and ellipsoids R_t =ℰ(p_t,Q_t).
The resulting ellipsoid is


    m̃(ℛ_t,u_K,t) =  ℰ(h(z̅_t)+μ(z̅_t),H_tQ_tH_t^T) ⊕ℰ(0,Q_d̃(ℛ_t,u_K,t)),


where z̅_̅t̅ = (p_t,k_t)^T and H_t = A_h + B_h K_t. The set ℰ(0,Q_d̃(ℛ_t,u_K,t)) is obtained similarly to (<ref>) as the ellipsoidal over-approximation of

    0 ± [β_n σ(z̅) + L_∇ hl^2(ℛ_t,S_t)/2  + L_gl(ℛ_t,S_t)],

with S_t = [I_n_x,K_t^T] and l(ℛ_t,S_t) = max_x ∈ℛ_t||S_t(z(x)-z̅_̅t̅)||_2  , z(x) = (x,u_K_t(x))^T. The theoretical results of <ref> and <ref> directly apply to the case of the uncertainty propagation technique (<ref>). For the remainder of this paper, we assume K_t is pre-specified, while k_t is assumed to be a decision variable. For the sake of generality, we drop the subscript K and the functional dependency on x in u_K,t(x) unless required and refer to (<ref>) when writing u_t.



 §.§ Safety constraints

The derived multi-step ahead prediction technique provides a sequence of ellipsoidal confidence regions around trajectories of the true system f through (<ref>). We can guarantee that the system is safe by verifying that the computed confidence ellipsoids are contained inside the polytopic constraints (<ref>) and (<ref>).
That is, given a sequence of feedback controllers u_K,t,   t = 0,..,T-1 we need to verify

    ℛ_t+1⊂𝒳,   u_t(ℛ_t) ⊂𝒰,   t=0,..,T-1,

where (R_0,..,R_T) is given through (<ref>) and u_t(R_t) := {u_K,t(x) | x ∈ R_t}. Since our constraints are polytopes, we have that

𝒳 = ⋂_i=1^m_x𝒳_i,

𝒳_i = {x ∈^p |  [H_x]_i,·  x - h_i^x≤ 0 },
where [H_x]_i,· is the ith row of H^x.
We can now formulate the state constraints through the condition ℛ_t = ℰ(p_t,Q_t) ⊂𝒳 as m_x individual constraints R_t ⊂𝒳_i,  i= 1,..,m_x, for which an analytical formulation exists <cit.>,

    [H_x]_i,·p_t + √([H_x]_i,·Q_t[H_x]_i,·^T)≤ h^x_i,

 ∀ i ∈{1,..,m_x}. Moreover, we can use the fact that u_t is affine  in x to obtain u_t(ℛ_t) = ℰ(k_t,K_tQ_t,K_t^T), using (<ref>). The corresponding control constraint u_t(ℛ_t) ⊂𝒰 is then equivalently given by

    [H_u]_i,·k_t + √([H_u]_i,·K_tQ_tK_t^T [H_u]_i,·^T)≤ h^u_i,

∀ i ∈{1,..,m_u}. This provides us with a closed-form expression of our safety constraints (<ref>) that deterministically guarantees the safety of our system over an arbitrary finite horizon T, given that the system is contained in the sequence of ellipsoids ℛ_t, t=0,..,T. Hence, these constraints are as reliable as our multi-step ahead prediction technique and, consequently, as reliable as our statistical model.






§ SAFE MODEL PREDICTIVE CONTROL


Based on the previous results, we formulate a MPC scheme that is guaranteed to satisfy the safety condition in <ref>:




    u_0,..,u_T-1minimize       J_t(ℛ_0,..,ℛ_T) 
       subject to       ℛ_t+1 = m̃(ℛ_t,u_t),   t=0,..,T-1   
             ℛ_t ⊂𝒳,   t = 1,..,T-1 
              u_t(ℛ_t) ⊂𝒰, t = 0,..,T-1 
             ℛ_T ⊂𝒳_safe,



where R_0 := {x_t} is the current state of the system and the intermediate state and control constraints are defined in (<ref>)and (<ref>), respectively. The terminal set constraint ℛ_T ⊂𝒳_safe has the same form as (<ref>) and can be formulated accordingly. For now, we assume an arbitrary objective function and discuss how to choose J_t to solve a RL task in <ref>.

Due to the terminal constraint ℛ_T ⊂𝒳_safe, a solution to  (<ref>) provides a sequence of feedback controllers u_0,..,u_T that steers the system back to the safe set 𝒳_safe.
We cannot directly show that a solution to MPC problem (<ref>) exists at every time step (this property is known as recursive feasibility) without imposing additional assumptions. This is mainly due to the fact that we cannot (not even with high probability) guarantee that, after transitioning to the next time step and shifting our horizon t ← t + 1, and hence, ℛ_T-1←ℛ_T, there exists a new state feedback controller u_T-1 such that m̃(R_T-1,u_T-1) = ℛ_T ⊂𝒳_safe, u_T-1(ℛ_T-1) ⊂𝒰. While it may be possible to enforce this in general, we would have to, among other difficulties, carefully deal with the nonlinearity and non-convexity of our MPC problem, as e.g. in <cit.>, and the fact that our terminal set is not necessarily RCPI, which is a pre-requisite in many robust MPC approaches <cit.>. However, employing a control scheme similar to standard robust MPC, we guarantee that such a sequence of feedback controllers exists at every time step as follows: Given a feasible solution Π_t = ( u_t^0,.., u_t^T-1 ) to (<ref>) at time t, we apply the first control u_t^0. In case we do not find a feasible solution to (<ref>) at the next time step, we shift the previous solution in a receding horizon fashion and append π_safe to the sequence to obtain Π_t+1 = ( u_t^1,.., u_t^T-1, π_safe ). We repeat this process until a new feasible solution exists that replaces the previous input sequence. This procedure is summarized in <ref>.

We now state the main result of the paper that guarantees the safety of our system under the proposed algorithm.


Let π be the controller defined through <ref>, x_0 ∈𝒳_safe and let the planning horizon T be chosen as in <ref>.
Then the system (<ref>) is δ-safe under the controller π.


From <ref>, the ellipsoidal outer approximations (and by design of the MPC problem, also the constraints (<ref>)) hold uniformly with high probability for all closed-loop systems f_Π, where Π is a feasible solution to (<ref>), over the corresponding time horizon T given in <ref>.

Hence we can show δ-safety by induction.

Base case: If (<ref>) is infeasible, we are δ-safe using the backup controller π_safe of <ref>, since x_0 ∈𝒳_safe. Otherwise the controller returned from (<ref>) is δ-safe as a consequence of <ref> and the terminal set constraint that leads to x_t+T∈𝒳_safe.

Induction step: let the previous controller be δ-safe. At time step t+1, if (<ref>) is infeasible then Π_t leads to a state x_t+T∈𝒳_safe, from which the backup-controller is δ-safe by <ref>. If (<ref>) is feasible, then the return path is δ-safe by <ref>.






In <ref>, we implicitly make the assumption that we have access to the true state of the system during planning, while updating the model with a noisy observation of the state. This is a standard assumption when acting in a Markov decision process (MDPs). The theoretical guarantees of <ref> can be extended to noisy observations of the state by either assuming that the true state is contained in a compact set B_0 and setting R_0 ⊃ B_0, or by jointly bounding the probability of the state to be contained in an ellipsoid around the current observation for all time steps.





§ MPC-BASED SAFE REINFORCEMENT LEARNING


A sequence of ellipsoids computed with the proposed uncertainty propagation technique contains a trajectory of the system with high probability. While we can, by using this technique in combination with <ref>,  guarantee the safety of our system, we are ultimately interested in safely identifying our system through exploration and to solve a task in a safe reinforcement learning setting. To this end, we have to find answers to the following questions: Firstly, how can we best approximate and minimize the long-term cost (<ref>) of our model predictive controller given our current knowledge of our system? And secondly, how can we gather task-relevant observations from our system that help us in solving the given RL task?



 §.§ Safety and Performance


We assume that we are given an analytic, twice differentiable cost-function c(x_t,u_t), e.g. the squared distance to a desired set-point, for which we try to optimize the accumulated long term cost defined in equation (<ref>). 
We could directly apply <ref> to solve this RL task, e.g., by choosing the objective function J_t = ∑_t=0^T-1γ^t c(p_t,k_t) for the MPC problem (<ref>), where p_t is the center of the ellipsoid R_t =  ℰ(p_t,Q_t) in every step of the uncertainty propagation.

The issue with this approach is two-fold. Firstly, the proposed uncertainty propagation technique is designed to robustly outer-approximate the reachability of our system in order to guarantee safety. This stands in contrast to many of the techniques proposed in the literature which provide a stochastic estimate of the state of the system in terms of probability distributions using our statistical model <cit.>. These estimation techniques give rise to more powerful stochastic planning approaches.

Secondly, we argue that a safety maneuver acts on a different (typically smaller) time scale than is necessary to estimate the long-term cost (<ref>) of applying an action in the current time step.
Imagine a race car driving on a track. Estimating the impact of a steering maneuver before a turn on the velocity of the car after the turn typically requires a much longer look-ahead horizon than planning an emergency braking maneuver that decelerates the car towards a modest velocity.  



 §.§ Performance trajectory planning

To avoid this issue, we propose to plan a probabilistic performance trajectory X_0,..,X_H under a sequence of inputs u^perf_0,..,u^perf_H-1 using a stochastic performance-model m_perf. The goal of the performance trajectory is to minimize an approximation of the cumulative cost (<ref>),

    V^u^perf_0,..,u^perf_H-1(x_0)    = 𝔼[∑_t=0^H-1γ^t c(X_t,u^perf_t(X_t)) | X_0 = x_0], 
       ≈∑_t=0^∞γ^t c(x_t,u^perf_t(x_t)),

where x_t+1 = f(x_t,u_t) is a trajectory of the real, deterministic system (<ref>) and X_t+1∼ m_perf(X_t,u^perf_t(X_t)), t=0,..,H-1 is a sequence of random variables[While in our case we assume the system to be deterministic and our performance model to be stochastic, we can treat a deterministic performance model as a special case thereof.] that represents our incomplete, stochastic belief about the future states of our system.


 §.§ Maintaining safety through simultaneous planning

In order to maintain the safety of our system, we need to combine the performance trajectory with our provably safe MPC <ref>. A simple solution would be given by the following two-stage process at every time-step:

    
  * Find a sequence of performance controls û^perf_0,..,û^perf = min_u^perf_0,..,u^perf_H-1   V^u^perf_0,..,u^perf_H-1(x_0)
    
  * Solve the MPC problem (<ref>) with objective function J_t = ∑_i = 0^min(H,T)dist(u_i,û^perf_i),

where dist(·,·) is an appropriate distance function (e.g. the squared l_2 distance). By applying this scheme, e.g. similar to the one proposed in <cit.>, we could easily maintain the safety of our system. However, proximity of two actions u_i,û^perf_i in the action space does not necessarily imply similar performance, as we show in <ref>.


Consider a simple system f(x_t,u_t) = x_t + u_t, u_t ∈{-1,0,1}, x ∈ in a discrete state-action space without model error, i.e. g ≡ 0. Let the cost-function be given by c(x) = 0,  ∀ x ∈∖{-1,1},  c(-1) = -2,   c(1) = - 1. We regard the case where our state constraints and safe set are given by 𝒳 = 𝒳_safe =. Starting in x_0 = 0 and using planning horizons T=H≥1, the optimal action in terms of performance would lead to the unsafe but low-cost state x_1 = -1, i.e. u_0^perf = -1. The closest safe action we can take is to apply u_0 = 0, resulting in the system being stuck in  x_t = 0  ∀ t ∈, although u_0 = 1 would directly lead to the best possible safe state.  

While the given example is highly simplified, it illustrates the fundamental issues that arise when a safety mechanism without an understanding of the task is combined with an RL agent that is ignoring the constraints. 
To approach this problem, we propose to simultaneously plan the performance and safety trajectory. We maintain the safety of our system by coupling both trajectories, enforcing that u^perf_0 = u^safe_0. That is, we apply the best possible action that is able to return the system to a safe state by means of a feasible safety trajectory.
This extended optimization problem is given by

    u_0,..,u_T-1
     u^perf_0,..,u^perf_H-1minimize       V^u^perf_0,..,u^perf_H-1(x_0) 
       subject to       (<ref>)-(<ref>),   t = 0,..,T-1 
              X_t+1  = m_perf(X_t,u^perf_t), t = 0,..,H-1 
              u_t = u^perf_t,   t = 0,..,r-1,

with r ≥ 1, possibly allowing even more than the first planned input to be identical.



 §.§ Exploration versus Exploitation

A solution to the MPC problem (<ref>) provides us with a control input that minimizes an estimate (<ref>) of the long-term consequences of applying this input to the real system. In order to improve this estimate, we need a better approximation of our system based on observations that are used to update our statistical model. However, by just following a control policy that maximizes this estimated objective, we may never obtain crucial task-relevant information about our system in previously unseen regions of the state space. Hence, finding a good exploration strategy is crucial for many real-world RL applications.

Since our statistical model knows what it does not know through the predictive variance Σ_n, we get a good estimate on which regions of the state space need to be explored. Consequently, a wide range of trajectory-based RL approaches and their corresponding exploration techniques attempt to find an efficient trade-off between exploration and exploitation by combining an estimate of the expected cost (<ref>) and the confidence intervals in the objective function, e.g. <cit.>. We can incorporate any of these exploration techniques in the objective function of our extended MPC problem (<ref>).







§ PRACTICAL CONSIDERATIONS


<ref> theoretically guarantees that the system remains safe, while actively optimizing for performance via the MPC problem (<ref>). This problem can be solved by commonly used, nonlinear programming (NLP) solvers, such as the Interior Point OPTimizer (Ipopt, <cit.>).
We consider possible design choices that could improve the performance in a practical application.



 §.§ Optimizing over affine feedback policies


In practice, the affine feedback control structure introduced in <ref> improves performance significantly. However, optimizing over K_t in (<ref>) seems to be challenging, both in terms of numerical stability and computational complexity. Hence, we pre-specify the feedback terms in all of our experiments and only optimize over the feed-forward terms.



 §.§ Lipschitz constants and eigenvalue computations


In our multi-step ahead predictions (<ref>), we need to solve a generalized eigenvalue problem for every step in the planning horizon. We use the inverse power iteration, an iterative method that asymptotically converges to the largest generalized eigenvalue of a pair of matrices <cit.>. We run the algorithm for a fixed number of p^2 iterations to solve these intermediate eigenvalue problems. In practice, this seems to result in sufficiently accurate estimations.  


Due to the generalized eigenvalue problem, the uncertainty propagation (<ref>) is not analytic.
However, we can still obtain exact derivative information by means of algorithmic differentiation, that is provided in many state-of-the-art optimization software libraries <cit.>.







§ EXPERIMENTS



In this section, we evaluate the proposed safe MPC <ref> to safely explore the dynamics of an inverted pendulum system and the corresponding episode-based RL <ref> to solve a task in a cart-pole system with safety constraints. We provide the code to run all experiments detailed in this section on Github [<https://github.com/befelix/safe-exploration>.]


For our experiments we consider an episodic setting, where we interact with the system for n_steps time steps, e.g. the time it ideally takes for an autonomous race car to finish a lap, and then use the observations gathered from the system after each rollout to update our statistical model. This procedure is repeated n_ep times. We can directly replace (<ref>) with our new MPC problem (<ref>) in <ref> and apply it to our system during each rollout. This procedure is depicted in <ref>.

The difference between the episodic setting and the MPC <ref> lies mainly in the choice of the objective. However, a subtle but theoretically important difference is that we only update our model every n_steps time steps. Hence, we technically require that the scaling factor β for our statistical model holds for sample sizes T · n_steps in <ref>. Then, the safety guarantees of <ref> are maintained throughout our learning process in <ref> by noting that feasibility of the problem (<ref>) still guarantees existence of a safe return strategy and we can still use the same fall back strategy as before in case of infeasibility. We note that, as long as we account for a possible delay in updating our statistical model, our approach is not limited to episodic RL tasks only.

We use a GP with a mixture of linear and Matérn kernels for both systems. For the inverted pendulum system, we initially train our model with a dataset (𝒵_0,ỹ_0) sampled inside the safe set using the backup controller π_Safe. That is, we gather n_0 = 25 initial samples 𝒵_0 = { z_1^0,..,z^0_n_0} with z_i^0 = (x_i, π_safe(x_i)), x_i ∈𝒳_safe,   i = 1,..,n and observed next states f̂_0 ={f̂^1_0,..,f̂^n_0_0}. In the cart-pole RL task no prior observations from the system are available. The theoretical choice of the scaling parameter β_n,T for the confidence intervals in (<ref>) can be conservative and we choose a fixed value of β_n,T = 2 instead, following <cit.>. For improved numerical stability and computational efficiency, we limit the number of training points used to update the GP to 150 and use the maximum variance selection procedure (see e.g. <cit.>) to sub-select the most informative samples in cases where more samples are available, however there exist more sophisticated, provably near-optimal selection procedures <cit.>.



For experiments that employ a performance trajectory m_perf as introduced in <ref>, we choose a common uncertainty propagation technique, that approximates the intractable posterior of the GP statistical model with uncertain, Gaussian inputs  by another Gaussian. That is, given a Gaussian distributed probabilistic belief of our system X_t ∼𝒩(m_t,S_t), we approximate the intractable integral 

    p(·|m_t+1,S_t+1) = 𝔼_x_t ∼ X_t[ p_𝒢𝒫(·|x_t,u_t)]

with the Gaussian X_t+1∼𝒩(m_t+1,S_t+1), where p_𝒢𝒫(·|x_t,u_t) is the density function of the GP posterior for a deterministic input (<ref>), (<ref>). Under a sequence of control inputs[Many approaches in the literature, including the one we use, also allow for linear state-feedback controllers as inputs.], this results in a sequence of Gaussian random variables as a Bayesian approximation of a trajectory of the true system. We choose a technique where m_t+1 and S_t+1 are given by a Taylor approximation of the first and second moment of the posterior, respectively <cit.>, but we note that a number approaches exist in this category that could be readily used. 

We try to actively encourage exploration by employing the saturating cost function c_sc(x) = 1 - exp(-1/2(x-x_g)^TW(x-x_g)) as it behaves similarly to a squared distance function when close to the target x_g and prefers uncertain states when being far away from the target state <cit.>.
Since our stochastic performance trajectory is given by a sequence of Gaussians, we can compute the corresponding surrogate objective,

    𝔼[∑_t=0^H-1γ^t c_sc(X_t,x_g) | X_0 = x_0],

in closed form. We can now replace (<ref>) with this objective function in our MPC formulation. 



 §.§ Inverted pendulum exploration


As an initial experiment, we aim to safely identify a partially unknown inverted pendulum system through exploration, without trying to solve a specific task. To do so, we attempt to iteratively collect the most informative samples of the system, while preserving its safety. To evaluate the exploration performance, we use the mutual information I(g_𝒵_n,g) between the collected samples 𝒵_n = {z_1,..,z_n}∪𝒵_0 and the GP prior on the unknown model-error g, which can be computed in closed-form <cit.>.

The continuous-time dynamics of the pendulum are given by ml^2θ̈ = gmlsin(θ) - ηθ̇ + u, where m = 0.15 kg and l=0.5m are the mass and length of the pendulum, respectively, η = 0.1 Nmsrad is a friction parameter, and g = 9.81 ms^2 is the gravitational constant. The state of the system x = (θ,θ̇) consists of the angle θ and angular velocity θ̇ of the pendulum. The system is underactuated with control constraints 𝒰 = { u ∈ |-1 ≤ u ≤ 1}. Due to these limits, the pendulum becomes unstable and falls down beyond a certain angle. The origin of the system corresponds to the pendulum standing upright. The prior model h is given by the linearized and discretized system around the origin, albeit with friction neglected the mass the pendulum being lower than for the true system as in  <cit.>. The safety controller π_safe is a discrete-time, infinite horizon linear quadratic regulator (LQR,<cit.>) of the true system f linearized and discretized around the origin with cost matrices
Q = diag([1,2]), R = 20. The corresponding safety region 𝒳_Safe is given by a conservative polytopic inner-approximation of the true region of attraction of π_safe. We do not impose state constraints, i.e. 𝒳 = ^2. However the terminal set constraint (<ref>) of the MPC problem (<ref>) acts as a stability constraint and prevents the pendulum from falling.


  §.§.§ Static Exploration


For a first experiment, we assume that the system is static, so that we can reset the system to an arbitrary state x_n ∈^2 in every iteration. In the static case and without terminal set constraints, a provably close-to-optimal exploration strategy is to, at each iteration n, select state-action pair z_n+1 with the largest predictive standard deviation <cit.>


    z_n+1 = z ∈𝒳×𝒰max∑_1 ≤ j ≤ p σ_n,j(z),


where σ_n,j^2(·) is the predictive variance <ref> of the jth 𝒢𝒫(0,k_j) at the nth iteration. Inspired by this, at each iteration we collect samples by solving the MPC problem (<ref>) with cost function J_n = -∑_j=1^pσ_n,j(x_0,u_0), where we additionally optimize over the initial state x_0 ∈𝒳. Hence, we visit high-uncertainty states, but only allow for state-action pairs z_n that are part of a feasible return trajectory to the safe set 𝒳_Safe.

Since optimizing over the initial state is highly non-convex, we solve the problem iteratively with 25 random initializations to obtain a good approximation of the global minimizer. After every iteration, we update the sample set 𝒵_n+1 = 𝒵_n ∪{z_n}, collect an observation (z_n,f̂_n) and update the GP models. We apply this procedure for varying horizon lengths.

The resulting sample sets are visualized for varying horizon lengths T ∈{1,..,5} with 300 iterations in <ref>, while <ref> shows how the mutual information of the sample sets 𝒵_i,  i = 0,..,n for the different values of T. For short time horizons (T=1), the algorithm can only slowly explore, since it can only move one step outside of the safe set. This is also reflected in the mutual information gained, which levels off quickly. For a horizon length of T=4, the algorithm is able to explore a larger part of the state-space, which means that more information is gained. For larger horizons, the predictive uncertainty of the final state is too large to explore effectively, which slows down exploration initially, when we do not have much information about our system. The results suggest that our approach could further benefit from adaptively choosing the horizon during operation, e.g. by employing a variable horizon MPC approach <cit.>, or by increasing the horizon when the mutual information saturates for the current horizon.






  §.§.§ Dynamic Exploration


As a second experiment, we collect informative samples during operation; without resetting the system at every iteration. Starting at x_0 ∈𝒳_safe, we apply <ref> over 200 iterations.
We consider two settings. In the first, we solve the MPC problem (<ref>) with -J_n given by (<ref>), similar to the previous experiments. In the second setting, we additionally plan a performance trajectory as proposed in <ref>. We define the cost-function -J_t = ∑_t=0^Htrace(S_t^1/2) - ∑_t=1^T (m_t-p_t)^TQ_perf(m_t-p_t), which maximizes the sum of predictive confidence intervals along the trajectory m_1,..,m_H, while penalizing deviation from the safety trajectory. We choose r = 1 in the problem (<ref>), i.e. the first action of the safety trajectory and performance trajectory are the same. As in the static setting, we update our GP models after every iteration.

We evaluate both settings for varying T ∈{1,..,5} and fixed H=5 in terms of their mutual information in <ref>. We observe a similar behavior as in the static exploration experiments and get the best exploration performance for T=4, with a slight degradation of performance for T=5 after 200 iterations. By comparing the exploration performance between iteration 50 and 200, we can see that influence of longer return trajectories on the exploration performance only comes into play after a certain number of iterations. This can be seen by comparing the similar performance of T=3 and T=4 after 50 iterations with the significantly improved performance for T=4 after 200 iterations. The setting T=3 during the same period only sees modest performance improvements. We can see that, except for T=1, the performance trajectory decomposition setting consistently outperforms the standard setting. Planning a performance trajectory (green) provides the algorithm with an additional degree of freedom, which leads to drastically improved exploration performance.





 §.§ Cart-pole safe reinforcement learning


We tackle a cart-pole balancing task, where and underactuated cart learns to balance a pendulum in upright position while moving along a rail from starting position to a goal state. We initialize the pendulum in an upright position, with an initial position of the cart at x^cart_0 = -2  and we want to drive the system to the goal-position at x^cart_g = 2.6. We limit the length of the rail by x^cart∈ [-10,3.0] and simulate a floor, the pendulum is not allowed to hit, i.e. θ∈ [-90,90] degree, where θ is the angle of the pendulum.
 The equations of motion of the system are given by

    (M+m)ẍ^cart - mlθ̈cosθ+ mlθ̇^2sinθ = 0, 
    
         mlθ̈-mẋ^cartcosθ  -mgsinθ = u - ηẋ^cart,

where the velocity of the cart and the angular velocity are given by ẋ^cart and θ̇, respectively. The mass of the cart is M = 0.5 kg and the corresponding mass of the pendulum is given by m = 0.5 kg. The gravitational constant is again given by g = 9.81  ms^2 and the friction coefficient of the rail is η = 0.1. The origin of the system describes the cart at rest with the pendulum in an upright position The cart is strictly limited in its capability to overshoot the goal position due to the constraint x^cart < 3. The angular constraints θ∈ [-90,90] require the system to remain stable while driving to the goal-position. The actuators of the cart are limited to u ∈ [-5,5]. Again, we use a LQR controller based on the linearized true system with cost-matrices Q= ([4,8,12,2]), R = 40  as the safety controller π_safe.

We apply the episode-based safe RL <ref> over 8 episodes, where in each episode we interact with the system over n_steps = 50 time steps. We then measure the performance of a rollout x_0,..,x_n_steps via the accumulated squared distance to the goal C_ep =  ∑_t=0^n_steps 0.1 (x^cart_t-x_g)^2. All reported results are averaged over 6 repetitions of the corresponding experiment.

We compare the performance of the SafeMPC algorithm with and without the probabilistic planning module introduced in <ref>. In the settings where we use a performance trajectory, we set the number of planning steps to H=15 and choose the expected squared exponential cost function (<ref>) as the performance objective.
In the case of no performance trajectory, i.e. H=0, we choose the objective J_t = ∑_t=0^T-1 c_rl(p^cart_t-x_g)^2, where p_0,..,p_T-1 are the centers of the ellipsoids in the safety trajectory at each time step. 




  §.§.§ Influence of the performance trajectory


We report the performance of the last episode for T ∈{1,2,3,4} and H ∈{0,15} in <ref>.
In <ref>, we compare the summed cost C_ep for each episode with H = 15 and T ∈{1,2,3,4}.     

As for the exploration experiments, we can see that planning a performance trajectory greatly improves the performance of our system. In this RL task, the agent even fails to find any reasonable strategy when planning without using the additional performance trajectory to estimate the long-term utility of its actions. In the settings with T ∈{1,2}, the performance seems to improve quite steadily over the course of episodes. Whereas, for T ∈{3,4}, the performance seems to fluctuate more.  This could be due to the system switching to the safety trajectory more often, as longer safety trajectories allow more aggressive maneuvers from which the system then has to recover again. This seems to be harming overall performance. In future work, we want to investigate how to alleviate the negative effects of this switching control design. We note that, depending on the design of the RL problem, much longer horizons of the safety trajectory may be beneficial.

























  §.§.§ Baseline comparison

We compare our approach to a baseline, where we set T=0, i.e. we remove the safety trajectory and only consider the performance trajectory with the objective (<ref>). We use probabilistic chance constraints to approximately enforce the state and control boundaries, as detailed in <cit.>. We consider different lengths of the performance trajectory H ∈{5,10,15} and report the percentage of rollouts that were successful and the accumulated cost for the different settings in <ref>. We can see that the system crashes in a large number of rollouts in all three settings, where for H=5 the system crashes almost constantly. There also seem to be two types of shortcomings that lead to safety failures. For short horizons (H=5) the problem seems to lie in the fact that we do not account for the long-term effect on safety when taking an action, whereas for long horizons (H=15), finding a feasible solution to the MPC problem seems to be an issue, leading to safety failures due to the missing fallback strategy.

We note that the computational complexity of both approaches is dominated by the time to compute the GP predictive uncertainty (<ref>) along the planned trajectories in each iteration of the solver. Hence, when carefully optimized for runtime, it should be indeed possible to execute our proposed algorithm in real-time <cit.>.  








§ CONCLUSION


We developed a safe learning-based MPC framework that can solve reinforcement learning tasks under state and control constraints. We proposed two reliable uncertainty propagation techniques that underlie the algorithm. By combining the safety features of the learning-based MPC framework with techniques from model-based RL, we can guarantee the safety of the system while learning a given task. We experimentally showed that our proposed RL algorithm is capable of learning a task in a simulated cart-pole system without falling over or overshooting the end of the rail. 























 
















































































































































































plain        
 











Torsten Koller
Torsten Koller received the B.S. in mathematics from the Technical University Munich (TUM), Munich, Germany and the M.S. in computer science from the University of Freiburg, Freiburg, Germany, in 2014 and 2018, respectively. He is currently working towards the Ph.D. degree in computer science at the Neurorobotics research group with the University of Freiburg. He was a visiting scholar at ETH Zurich, where he conducted his master thesis research. 

His current research interests lie at the intersection of nonlinear robust (model predictive) control and machine learning.



Felix Berkenkamp
Biography text here.


Matteo Turchetta
Biography text here.


Joschka Boedecker
Biography text here.


Andreas Krause
Biography text here.

















