title: Safe Active Dynamics Learning and Control: A Sequential Exploration-Exploitation Framework
authors: Thomas Lew Apoorva Sharma James Harrison Andrew Bylard Marco Pavone
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
      
  
  
  
  
  
  theoremTheorem
  corollaryCorollary
  lemmaLemma
  assumptionAssumption
  
  
  definitionDefinition
   
  
  
  
  
   
  
  
   
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
   
   
   
  
  
  
  
  
  
  
  
  
   
   
   
   
   
   
   
   
  
   
  
   
   
   
   
  
  
  
   
   
   
   
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
   
   
   
  
  
  
  
  
   
   
   
   
   
   
  
  
   
   
  
  
  
  
  
  
  
  
   
  
  
  
   
   
   
   
   
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
   
   
  
  
  
  
  
  exampstyle
    1em plus .2em minus .1em
    1em plus .2em minus .1em
  
  
     
     
     
    . 
    .5em 
     
    
  
  
  
  mydefDefinition
  
  mylemmaLemma[section]
  
  exampstyle
  mycorollaryCorollary[section]
  
  exampstyle
  myproblemProblem
  
  exampstyle
  
  preremark3Theorem
  mytheorem
  
  
   
   
  exampstyle
  mypropositionProposition[section]
  
  exampstyle
  mypropertyProperty[section]
  
  exampstyle
  myassumptionAssumption
  
  
  
  
  
  
  
  
   
   
  
  
  
    
    
    
    
    
    
    
  
  
  
  
  
  
  
  
  C>p35mm
  
  L>p15mm
  S>p8mm
  
  H>p14mm
  G>p16mm
  F>p12mm
           
     
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    @red@titletitle
    title
      -2em
      titleheader1.5em@red@title
  
  
  
  
  
  This paper has been accepted for publication in the
  
  IEEE Transactions on Robotics (T-RO).
  2022 IEEE
  
  
  
  
  
  
  
  
  Safe Active Dynamics Learning and Control:  
  A Sequential Exploration-Exploitation Framework
      Thomas Lew, Apoorva Sharma, James Harrison, Andrew Bylard, and Marco Pavone
  
  
  
  
  The authors are with 
  
  Stanford University, Stanford, CA 94305-4035 USA. Contact: .
  
      December 30, 2023
  ========================================================================================================================================================================
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Lew et al.: Safe Active Dynamics Learning and Control: A Sequential Exploration-Exploitation Framework
  
  
  Safe Active Dynamics Learning and Control:  
  A Sequential Exploration-Exploitation Framework
      Thomas Lew, Apoorva Sharma, James Harrison, Andrew Bylard, and Marco Pavone
  
  
  
  
  The authors are with 
  
  Stanford University, Stanford, CA 94305-4035 USA. Contact: .
  
      December 30, 2023
  ========================================================================================================================================================================
  
  
  
  
  
  
  
  
  Safe deployment of autonomous robots in diverse scenarios requires agents that are capable of efficiently adapting to new environments while satisfying constraints. 
  In this work, we propose a practical and theoretically-justified approach to maintaining safety in the presence of dynamics uncertainty. 
  Our approach leverages Bayesian meta-learning with last-layer adaptation. The expressiveness of neural-network features trained offline, paired with efficient last-layer online adaptation, enables the derivation of tight confidence sets which contract around the true dynamics as the model adapts online. We exploit these confidence sets to plan trajectories that guarantee the safety of the system. Our approach handles problems with high dynamics uncertainty, where reaching the goal safely is potentially initially infeasible, by first exploring to gather data and reduce uncertainty, before autonomously exploiting the acquired information to safely perform the task. 
  Under reasonable assumptions, 
  we prove that our framework guarantees the high-probability satisfaction of all constraints at all times jointly, i.e. over the total task duration. 
  
  This theoretical analysis also motivates two regularizers of last-layer meta-learning models that improve online adaptation capabilities as well as performance by reducing the size of the confidence sets. 
  
  We extensively demonstrate our approach in simulation and on hardware.
  
  
  
  
  
  
  
  
  
  
  
  Robotics, meta-learning, chance-constrained planning, system identification, dynamics, reachability analysis.
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  Deploying autonomous robotic systems in safety-critical applications requires agents that are capable of adapting to partially unknown environments while satisfying constraints. 
  For example, an autonomous robot assisting astronauts in space must be able to handle a priori unknown payloads while respecting velocity constraints and avoiding collisions with obstacles. 
  This requires quantifying the uncertainty in the environment and factoring it into the decision process to take reliable actions. 
  As the initial levels of uncertainty may be very high as the agent encounters a new environment, this also requires choosing actions that actively reduce the level of uncertainty to eventually be capable to safely perform the given task. 
  For instance, a robot grasping a new payload should first identify the properties of its coupled dynamics before attempting to transport the payload to a destination. 
  
  
  Despite recent rapid progress on facets of this problem in the fields of data-driven control, machine learning, and reinforcement learning, the overall problem of safely controlling a robotic system in uncertain environments remains a challenge. Specifically, it is difficult to (1) obtain statistical models which are general yet offer a tight characterization of uncertainty, (2) design control laws that formally guarantee the safety of uncertain nonlinear dynamical systems, and (3) ensure that these guarantees hold even as the system adapts and reduces uncertainty using online data. 
  
  Designing a robotic system that is capable of actively and safely exploring its environment  without human intervention remains an open challenge. 
  
  
  
  
  
  
  
  
  
  
  
  
  Contributions: We propose a practical framework for the safe learning and control of uncertain nonlinear systems. 
  Our approach is capable of tackling chance-constrained problems that are initially infeasible due to high dynamics uncertainty by autonomously exploring and learning online 
  to reduce uncertainty before performing the task, as shown in Figure <ref>. 
  This framework is enabled by the following contributions:
  
  
    
      
    * Dynamics modeling: 
      We show that last-layer meta-learning models can learn tight and structured  representations of uncertain dynamics  from data and can efficiently reduce uncertainty online.  Critically, we derive adaption guarantees of such models in the form of confidence sets 
      guaranteed with high probability to contain the true dynamics throughout the online learning process. 
      These guarantees rely on two realistic and intuitive assumptions on the quality of the offline meta-training process which can be empirically validated before deploying the system. This theoretical analysis also suggests two regularizers to meta-train adaptive  models with improved performance and more expressive features. 
      
    * Autonomous exploration:  
      We combine these confidence sets with reachability-aware trajectory optimization to tackle problems with high levels of uncertainty that are initially infeasible. 
  
  Specifically, we use the current uncertainty level to autonomously switch from active learning, to reduce dynamics uncertainty (the exploration phase), to executing the desired task (the exploitation phase) once it is feasible. 
  
  
  
  
      
    * Safety: 
      In this work, we define safety in terms of the satisfaction of all constraints over the entire operation time. As such, 
      we prove that our framework guarantees the high-probability satisfaction of all constraints at all times jointly, i.e. over the total task duration. This contrasts with prior work which instead enforces pointwise chance constraints (see Section <ref>). 
      We also prove that the algorithm is continually recursively feasible with high probability,  
      which is crucial to enable autonomous operation. 
  
      Our guarantees are independent of the (unknown) total time required to solve the task, and safety still holds even if the perfect-information problem is impossible to solve, e.g., if the only path to the goal is blocked by an obstacle.
  
  
  We validate our approach on two challenging highly uncertain nonlinear systems 
  and show that we are able to autonomously and safely infer dynamics to accomplish goal-reaching tasks. 
  We also provide hardware experiments demonstrating the reliability and real-time capabilities of our approach, as well as the practical applicability of dynamics meta-learning models with last-layer adaptation.
  
  Organization: We first discuss related work in Section <ref>. Section <ref> describes our chance-constrained problem formulation. We present our contributions to the meta-learning model used to capture uncertain dynamics in Section <ref> and our control approach in Section <ref>.  
  Finally, we present results in Sections <ref> and <ref> and conclude in Section <ref>.
  
  
  
  
  Notation: 
  
  
  χ^2_d(p) denotes  the p-th quantile of the χ^2 distribution with d degrees of freedom. 
  For any ,∈^d and  a d× d positive definite matrix, _^2 =  ^⊤
  and λ̅() (λ()) denotes the maximum (minimum) eigenvalue of . 
  
  ⟨, ⟩ denotes the inner product of  with  (for ,∈^d,  ⟨, ⟩= ^⊤).  
  We denote a conjunction (logical 0.95AND) by ∧ and a disjunction (logical 0.95OR) by ∨.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § RELATED WORK
  
  
  
  
  Our approach is enabled by two key components: a dynamics meta-learning model with adaptation guarantees and a sequential exploration and exploitation control algorithm. 
  
  Choosing a model: In contrast to model-free approaches to learning-based control, model-based methods generally provide better sample efficiency while enabling guarantees on constraint satisfaction and stability <cit.>. 
  
  These model-based methods rely on the choice of a dynamics model—e.g., neural networks <cit.>, squared-exponential kernel Gaussian processes <cit.>, or linear models <cit.>—each with associated strengths and weaknesses.
  Recent work in the controls community 
  has leveraged behavioral systems theory to guarantee stability and probabilistic constraint satisfaction of a non-parametric model predictive control (MPC) scheme <cit.>. 
  Although such methods have been shown to perform well for nonlinear systems <cit.>, their guarantees currently do not extend beyond linear systems
  (see also <cit.> for a recent analysis of the nonlinear setting). 
  
  
  
  Alternatively, assuming 
  linearly parameterizable dynamics with known nonlinear basis functions allows the design of stable adaptive controllers <cit.> and planning algorithms with adaptation guarantees 
  <cit.>. 
  
  When such structure is not known a priori, purely data-driven  models such as neural networks can be used to derive nonlinear controllers that guarantee stability <cit.>. 
  However, these methods require collecting a dataset characterizing the system dynamics throughout the operational regime and would need retraining with a new dataset if the environment or the system changes. Generally, this 
  
  is prohibitively expensive in terms of data requirements. 
  
  
  
  To tackle problems where data is scarce, meta-learning <cit.> aims to train a model from data over different tasks, making it capable of rapid adaptation given limited data from a given new task.  
  
  In the context of dynamics learning, a task corresponds to a dynamical system and meta-learning consists of training a model capable of efficiently fitting the true system's dynamics with limited data. 
  
  Prior work has successfully applied meta-learning to learning-based control  <cit.>: performing gradient descent on meta-learned neural networks at run time enables fast adaptation with limited data. 
  
  However, the non-convexity of the online adaptation objective 
  makes it difficult to provide adaptation guarantees for such gradient descent approaches. Instead of performing gradient descent online, meta-learning has also been combined with adaptive controllers <cit.> yielding efficient closed-loop tracking performance. Further analysis is still necessary to provide formal performance guarantees for this approach in safety-critical applications.  
  
  
  
  
  To obtain performance guarantees, model adaptivity, and data efficiency,
  
  our approach 
  
  combines meta-trained neural network features with online last layer adaptation <cit.>:
  offline, data from related environments is used to train a neural network to produce features tailored to the structure of the system; 
  online, only adapting the last layer yields a sample-efficient adaptation process with strong guarantees on statistical performance. 
  
  By taking a Bayesian perspective on meta-learning, our model also yields a calibrated prior uncertainty representation that enables quantifying the initial uncertainty of the system. 
  
  Importantly, building on this initial uncertainty quantification, the linear structure of the last layer adaptation process enables constructing confidence sets for the parameters of the model that hold throughout the online learning process. 
  
  
  
  Within the realm of Bayesian modeling, a common modeling choice is Gaussian processes (GPs) <cit.>. GPs with squared-exponential kernels have been widely used to represent dynamics for learning-based control and exploration, as they can represent any bounded continuous function arbitrarily well given enough data. 
  While bounds providing similar guarantees as our confidence sets  
  can be derived for these models, such bounds are generally too conservative <cit.> 
  
  so that scaling constants are heuristically selected 
  in experiments <cit.>.  
  
  Alternatively, it is common to assume that the true system dynamics lie in a reproducing kernel Hilbert space (RKHS) determined by finite-dimensional basis functions 
  <cit.>. Importantly, this approach enables deriving tighter
  bounds over possible models <cit.> which we regularize for and use directly in this work. 
  
  In contrast to prior work, we explicitly learn these basis functions and regularize these features to improve generalization of the model.  Moreover, we quantify prior uncertainty in an offline meta-training procedure, creating a model which is both calibrated and expressive enough to represent possible systems. This also allows verifying that representation error is small offline before deploying the system, see Figure <ref>. 
  
  
  
  Safe dynamics learning and control  
  
  
  amounts to the design of a controller that guarantees satisfaction of all constraints. 
  Given a statistical model of the system subject to external random disturbances, it is common to use chance constraints to ensure that all specifications are satisfied with high probability (1 - δ), where δ > 0 is a small tolerable probability of failure. 
  Such chance constraints are often preferred to robust constraints (δ = 0), as statistical models predicting unbounded distributions over possible outcomes are typically used in such applications <cit.>, in which case requiring zero probability of failure is infeasible. 
  
  Typically, pointwise chance constraints  
  (_t ∈ _i)  ≥  (1-δ)
  are enforced for each i-th constraint _t ∈ _i at each time t ≥  0  <cit.>, where _t ∈ ^n denotes the state of the system at time t and _i ⊂ ^n denotes a subset of feasible states. 
  
  Unfortunately, this approach only guarantees that each i-th constraint is satisfied at each timestep t with high probability and does not guarantee safety over the entire operation time, i.e., that all constraints are satisfied at all times with high probability  <cit.>. 
  
  This is particularly an issue when the time to complete a given task is a priori unknown, which 
  
  motivates using a single joint chance constraint instead, taking the form (_t ∈ _i ∀ i ∀ t)  ≥  (1-δ). 
  
  Instead of allocating the total allowed probability of failure δ to each i-th constraint and each time t using Boole's inequality <cit.>, we take a frequentist viewpoint and use the confidence sets over the model parameters to plan robust trajectories that are safe with respect to all possible parameters in these confidence sets. 
  This approach is key to guaranteeing the satisfaction of all constraints trajectory-wise 
  with probability at least (1-δ) despite an unknown time to complete the task and  while switching between exploration to reduce dynamics uncertainty and exploitation to complete the given task.
   
  
  
  
  
  Reinforcement learning (RL) can also be effective for controlling uncertain systems <cit.>, and model-based methods in particular enable an agent to consider its uncertainty over dynamics when choosing actions 
  <cit.>. 
  However, standard model-based reinforcement learning methods do not provide sufficient guarantees for maintaining safety during operation as they commonly penalize constraint violation <cit.>, whereas safety-aware RL typically considers finite-dimensional state-action spaces, as using continuous state-action spaces remains an open area of research <cit.>. In contrast to traditional RL, we meta-learn a model and uncertainty representation offline, fine-tune it online via last-layer adaptation only, and use trajectory optimization and reachability analysis for planning and control. 
  
  This structure allows our framework to rapidly infer an accurate dynamics model while providing strong safety guarantees throughout the control task. 
  
  
  
  An outline of our active dynamics learning and control framework first appeared in <cit.>. However, it used pointwise chance constraints evaluated approximately and thus did not provide sufficient guarantees of safety. Subsequently, <cit.> developed a similar trajectory optimization approach to learn dynamics. 
  By replacing uncertain parameters with Brownian motion (which has independent increments), this approach makes an approximation that neglects time correlations of parameter uncertainty over the state trajectory. This limits the applicability of the theoretical analysis, as is discussed in more detail in <cit.>. 
  
  
  
  
  Learning-based control algorithms typically tackle the problems of stabilization around an operating point, 
  planning to a goal, 
  active system identification <cit.>, and 
  dual control <cit.>. 
  
  These approaches outperform methods that neither consider uncertainty nor leverage learning-based components. Still, they rely on the assumption that the problem is initially feasible given the current system's uncertainty characterization (e.g., learning-based MPC requires the first optimization problem to be feasible). 
  
  
  
  
  In contrast to prior work, we consider the problem of eventually performing a task which is potentially initially infeasible with the current information about the system. Actively learning the system's properties is thus necessary to reduce uncertainty before performing the task. Importantly, we explicitly guarantee the satisfaction of all safety constraints at all times jointly with high probability, 
  even if the problem is infeasible with perfect information. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § PROBLEM FORMULATION
  
  
  The goal of this work is to enable robots to safely navigate from an initial state (0)∈^n to a goal region ⊆^n despite highly uncertain dynamics while minimizing a chosen cost metric ℓ(·) (e.g., fuel consumption). 
  We write the state of the agent at time t ∈ ℕ as _t∈^ and the control input as _t∈^. 
  The system follows dynamics _t+1 = (_t, _t) + (_t,_t,) + _t, where  is known and corresponds to prior knowledge of the system, whereas  is unknown with unobserved parameters , 
  and _t are stochastic disturbances. 
  At the beginning of each episode j, the parameters are sampled ^j  ∼  p() and are fixed for the episode duration, so that the uncertain dynamics (·,·)+(·,·,^j) of the system are time-invariant. These parameters  correspond to unknown features that vary between episodes (e.g., the mass of a payload) and introduce correlations over the state trajectory that one must account for to guarantee safety. 
  We assume that the disturbances _t=(ϵ_1,t,,ϵ_n,t) are uncorrelated over t ∈ ℕ and i = 1,, 
  and that each ϵ_i,t is σ_i-subgaussian and bounded such that ϵ_i,t∈ℰ_i, where ℰ_i⊂ is a bounded subset.  
  
  Critically, this algorithm should guarantee
  trajectory-wise constraints satisfaction:  obstacles are avoided (_t∉ with ⊂^n the obstacle set) and system constraints are respected (_t∈, _t∈ with ⊂^n, ⊂^m the feasible state and control sets, which represent, for example, velocity bounds and torque-limited actuators, see Section <ref>) at all times, and the goal  is eventually reached. 
  Due to the system stochasticity and uncertain dynamics, enforcing all constraints with probability one 
  may be challenging without further assumptions (e.g., bounded model mismatch). 
  Instead, we enforce all constraints via a single joint chance constraint with probability level (1 - δ), where δ > 0 is a small tolerable probability of failure. 
  
  The full problem is expressed as
  
  
  
  Chance-Constrained Optimal Control Problem  ()
  
  
  
  
  
      min_,  ( ∑_t=0^T-1(_t, _t))
          s.t.    _0 = (0),
      
      _t+1 = (_t, _t) + (_t,_t,) + _t, 
        
      (
      ⋀_t=1^T_t ∈  ∧⋀_t=0^T-1_t ∈  ∧ _T ∈ ) ≥ 1 - δ
      ,
  
   
  where = (_0,,_T) and  
  = (_0,,_T-1) denote the state and control trajectories, 
  =  ∖  ⊆ ^n is the safe set, 
  ⊆  ⊆ ^n is the goal region, 
   T ∈ ℕ is the total problem duration (possibly infinite), 
   and (0) ∈ ^n is the initial state of the system.
  This problem formulation can be equivalently described as a constrained Markov decision process (CMDP) <cit.> with a continuous state and action space, nonlinear stochastic dynamics
  , and a non-convex cost function. 
  
  
  Satisfying trajectory-wise safety constraints with unknown dynamics is extremely difficult without further information. Thus, our approach relies on the following assumption <cit.>:
  
  
  
  [A1]
  The initial state (0) ∈ _0 ⊆, where _0 is a robust control invariant set for some feedback controller π:_0→: 
  
  if _t∈_0 and _t=π(_t), then 
  _t+1=(<ref>)∈_0 for all  ∼ p() and all _t with ϵ_i,t∈Ε_i.
  
  
  
  
  
  This assumption reflects prior knowledge that the system 
  is initially stable and satisfies all constraints 
  under a nominal controller π (e.g., regulated to a stable point using a simple linear feedback law). 
  
  
  Computing robust control invariant sets for general nonlinear uncertain systems (e.g., by leveraging additional prior knowledge of the structure and smoothness of the true system and of the support of the distribution p()) is an active area of research, see for instance <cit.>. 
  
  
  
  
  In this work, we make no assumptions on boundedness or smoothness properties of  (e.g., known Lipschitz constants). 
  Instead, we assume that we have access to a dataset of trajectories  generated from sampled parameters ^j. Such information may come from previous operation of a robot in similar environments or data generated from parameterized simulations. This motivates our use of meta-learning to encode this information and characterize the dynamics uncertainty.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § META-LEARNING & ADAPTATION GUARANTEES
  
  
  
   §.§ Bayesian meta-learning
  
  Our approach leverages a model for the unknown dynamics . To this end, we employ the Bayesian meta-learning architecture presented in <cit.> referred to as . Meta-learning (or “learning-to-learn”) <cit.> aims to train a model to be capable of rapid adaptation[The exact mechanism of “adaptation” to a task is at the core of the meta-learning algorithm. In MAML <cit.> it consists of a gradient step, in recurrent models it occurs via the hidden state dynamics <cit.>, and in <cit.> the update consists of Bayesian linear regression on the last layer.] within a distribution of tasks. 
  In this problem setting, a task corresponds to a dynamical system (·,·,) parameterized by the unknown parameters .
  models the unknown dynamics as =(ĝ_1,0.9,ĝ_n), where each ĝ_i: 0.94^n× ^m× ^d → is defined as
  
      ĝ_i(, , _̨i) = _̨i^⊤_i(,),    
      i = 1,,n,
  
  where 
  each 
  _i: ^n × ^m → ^ is a feed-forward neural network and _̨i ∈ ^ corresponds to the weights of the last layer. 
  The uncertainty in the space of dynamics functions is encoded through a normal distribution on each 
  last layer as _̨i ∼(_i, σ_i^2 _i), where _i∈^d denotes the mean parameters and _i∈^d× d the positive definite precision matrix. 
  As such, each component ĝ_i(,, _̨i) = _̨i^⊤_i(,) defines a Gaussian process (GP) with mean function 
  ̨̅_i^⊤_i(,) and kernel 
  σ_i^2 _i(,)^⊤_i_i(','). 
  As is common in the literature <cit.>, we opt for independent distributions for the parameters _̨i[I.e., we choose a Gaussian distribution for the parameters (_̨1,,_̨n), of mean (_1,,_n) and diagonal covariance with diagonal entries σ_i_i^-1. ]. This choice leads to substantial computational benefits during meta-training and to better performance at run-time over the alternative of  maintaining a joint distribution over all the parameters (_̨1,…,_̨n), see <cit.>. 
  
  
  
  
  
  The basis functions _i could be chosen as independently-parameterized neural networks. Instead, in this work, 
  we use fully-connected L-layer neural network basis functions 
  _i()=^L+1_itanh(^Ltanh∘⋯∘tanh(^1))[We omit biases in the notation for brevity; our implementation includes a bias term in each layer.]. Sharing the first L weight matrices across all n state dimensions allows for efficient scalability of the model and faster offline meta-training compared to training n separate neural networks _i. Choosing the bounded tanh activation functions yields smoother mean and variance predictions than common alternatives such as ReLU activations  <cit.>. 
  Finally, the last-layer matrices _i^L+1 provide an additional degree of freedom which enables each i-th dimension of the model to select relevant features
  [The matrices _i^L+1 
  also stabilize the offline meta-training process: by weighting intermediate features differently (e.g., when predicting quantities of different units, such as a position versus an angular velocity), this dimension-wise rescaling is not reflected in the precision matrices  _i which are then better conditioned. 
  By yielding a larger model class, this additional degree of freedom also allows further regularization to meta-train tighter uncertainty representations, see Section <ref>.].
  
  We found that using independent neural networks or sharing less than L weight matrices did not lead to better performance, although our analysis applies to such models as well. 
  
  
  
  
  
  
  
  
  
  
  
  
  The linear structure of this model allows for efficient online updates whose behavior is well understood.
  Given transitions from the system {(_0, _0, _1), …, (_t, _t, _t+1) } and prior parameters (_i,0, _i,0), 
  the normal distribution of the last-layer parameters is updated according to Bayesian linear regression: the posterior precision matrices and mean vectors are given by
  
      [              _i,t = Φ_i,t-1^⊤Φ_i,t-1 + _i,0,; _i,t = _i,t^-1 ( Φ^⊤_i,t-1_i,t + _i,0_i,0 ), ]   
          i = 1,,,
  
  where _i,t^⊤=  [x_i,1- h_i(_0,_0), , x_i,t- h_i(_t-1,_t-1)]  ∈ ^t and
  Φ_i,t-1^⊤=  [_i(_0, _0), , _i(_t-1,_t-1)]  ∈ ^× t. 
  Since each i-th model ĝ_i in (<ref>) is linear in the uncertain parameters _̨i, the posterior predictive state distribution 
  is also Gaussian. 
  Specifically, (_t+1|_0:t,_0:t) ∼(_t+1, _t+1), where the entries of _t+1 and the diagonal covariance matrix _t+1 are
  
      [ μ_i,t+1 = h_i(_t,_t) + _i,t^⊤_i(_t,_t) ,
         
      i = 1,,,;     Σ_i,t+1 = σ_i^2 ( 1 + _i(_t,_t)^⊤_i,t_i(_t,_t) ). ]
  
  
  
  
  
  
  Offline, this model is meta-trained on a dataset of trajectories corresponding to different system dynamics sampled from the distribution over possible systems. 
  By backpropagating through the posterior predictive distribution to learn the features _i and the prior parameters (_i,0, _i,0), 
  this model translates a dataset of uncertain trajectories into a learned feature space and a calibrated uncertainty characterization. 
  
  More precisely, given J different sampled dynamics and a meta-training horizon T_ML, we consider as meta-training objective the joint marginal log likelihood across the data
  
      ℒ(^0.6J_0.6T_ML; , , )
      =
      ∑_j=1^J ∑_t=0^T_ML-1log p(_t+1^j |_0:t^j,_0:t^j),
  
  where 
      ^J_T_ML = {(_t^j, _t^j, _t+1^j)_t=0^T_ML-1}_j=1^J is the meta-training dataset 
      and
      each j-th trajectory (_t^j, _t^j, _t+1^j)_t=0^T_ML-1 
      corresponds to the j-th sampled system (·,·,^j) parameterized by ^j  ∼  p().  
  Since each posterior distribution in (<ref>) is Gaussian as defined by (<ref>), the total meta-training loss 
  can be easily differentiated with respect to the model parameters using modern automatic differentiation packages <cit.>. During offline meta-training, we optimize this loss using stochastic gradient descent, evaluating (<ref>) on minibatches subsampled from ^J_T_ML.
  
  
  
  
  
  Online, using transition tuples { (_τ, _τ, _τ+1) }_τ=0^t-1 from the true system, we only adapt the last layer parameters _̨i using (<ref>). This restriction on the online behavior of the model allows deriving strong adaptation guarantees under reasonable assumptions on the quality of the offline training process. 
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Probabilistic adaptation guarantees
  
  Our first contribution is providing probabilistic adaptation guarantees for this model in the form of uniformly calibrated confidence sets which hold under the following assumptions, which are reasonable for common robotic applications. 
  
  
  
  
  [Capacity of meta-learned model)  (A2]
  For all , i  =  1,,n, there exists () ∈^ such that 
  
      ()^⊤_i(, )  
      =
      g_i(,,)   ∀∈, ∀∈.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  [Calibration of meta-learned prior) (A3]
  
  
  For ∼  p(), all i = 1,,, 
  and δ_i = δ/(2), 
  
  
      (
      () - _i,0__i,0^2 ≤σ_i^2 χ^2_(1-δ_i)
      )
      ≥
      (1 - δ_i).
   
  
  
  
  
  
  
  
  These two key assumptions on the quality of the offline meta-learning process state that 
  the meta-learning model is capable of fitting the true dynamics (A<ref>) and 
  that the prior uncertainty characterization is conservative (A<ref>).
  
   
   
   
   
   
   
   
   
   
   
   
  
  Since the features _i and the parameters (_i,0, _i,0) of the Gaussian prior 0.95_̨i  ∼ (_i,0, σ_i^2 _i,0) are chosen to maximize the log-likelihood of the dataset under ^j  ∼  p(), 
  
  the offline meta-training procedure presented in    <cit.> is well-suited to yield a model satisfying A<ref> and A<ref>. 
  
  Importantly, these assumptions can be empirically verified through predictive performance on a validation dataset, and techniques such as temperature scaling can be used to ensure calibration in a post-hoc manner <cit.>. 
  
  
  As discussed in details in Appendix  <ref>, A<ref> can be relaxed to hold approximately for some approximation error Δ_i>0: similar to <cit.>, the resulting analysis   would yield larger confidence sets but is beyond the scope of this work.  
  
  A<ref> states that offline meta-training yields a model with a prior variance σ_i^2 _i,0 that is large enough: with probability at least (1 - δ_i) over ^j ∼  p(), the true parameters () belong to the ellipsoidal (1 - δ_i)-credible sets predicted by the model. 
  We note that A<ref> is weaker than assuming that the model is calibrated (i.e., that A<ref> jointly holds for all values of δ_i), a topic extensively studied in the literature <cit.>. 
  
  
  
  
  
  In addition to enabling a theoretical analysis, we also use A<ref> and A<ref> to derive two regularizers improving the performance and adaptation capabilities of ; we refer the reader to Section <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
   
  These assumptions, together with the linear structure (<ref>) and adaptation rule (<ref>) of the meta-learned model, allow us to define adaptive confidence sets over the model parameters _̨i that are guaranteed to contain the true parameters  with high probability, even as the model adapts using subgaussian  process-noise-corrupted measurements.
  
  
  
  
  
  Consider the true system  (<ref>), modeled using the meta-learning model (<ref>), which is sequentially updated with online data from (<ref>) using (<ref>), leading to the updated parameters (_i,t, _i,t) for each dimension i=1,,. 
  
  Let
  
  
      β_i,t^δ = σ_i (
          √( 2 log(1/δ_i(_i,t)^1/2/(_i,0)^1/2) )+√(λ̅(_i,0)/λ(_i,t)χ^2_(1-δ_i)))
          
      and  _i,t^δ(_i,t,_i,t)  = {_̨i ∈ ^d  | _̨i -_i,t__i,t≤ β_i,t^δ}.
  
  
  Then, under Assumptions <ref> and <ref>, for δ_i=δ/(2),
  
      (
          ∈_i,t^δ(_i,t,_i,t)
            ∀ t ≥  0
          ) ≥ (1-2δ_i).
  
  
  
  
  
  
  
  
  
  
  A proof of this result, which follows from <cit.>, is available in the Appendix. 
  Using Boole's inequality,  the σ_i-subgaussianity of the disturbances _t, and Assumptions <ref> and <ref>, Theorem <ref> is derived by adapting results from the literature on linear contextual bandits <cit.>, which applies a stopping time argument to bound the probability of the bad event that the true parameters lie outside _i,t^δ at some time t. 
  
  
  
  Theorem <ref> provides confidence sets over model parameters which hold uniformly over all future times. 
  
  This is critical to ensure satisfying the trajectory-wise chance constraint 
  (<ref>) despite an unknown final time N. 
  
  To derive Theorem <ref>, we assume there exist fixed true parameters  corresponding to the true dynamics (A<ref>), and we take a frequentist viewpoint such that the confidence set _i,t^δ is a stochastic function of the transition tuples observed online. 
   By defining _i,t^δ through time-dependent values of β_i,t^δ, we ensure that the event that _i,t^δ excludes the true parameters  at any time t occurs with probability less than (δ/n). 
   
   
  
   
   
   
  
  
  Theorem <ref> is similar to results derived for kernel-based GPs <cit.> which rely on a good choice of kernel to satisfy assumptions similar to A<ref> and A<ref>; we discuss similarities in the next section. 
  Specifically, the scaling factor β_i,t^δ in Theorem <ref> also appears when analyzing kernel-based GPs, but the value of β of such models is often too large for practical use and is set to a lower value for experiments, see <cit.>. 
  In contrast, our 
  model operates within a finite-dimensional feature space and has values of β that are practically usable.
  
  
  
   §.§ Practical implications for better meta-training
  
  
  These theoretical results offer guidance on how to add regularization to the offline meta-learning process for the purposes of safe learning-based control. Specifically, we propose to add two regularizers to the offline meta-training of . 
  
  This regularization does not compromise on safety as long as the regularized model still meets Assumptions <ref> and <ref>.
  
  
  
  
  
  
  Orthogonal features are better: 
  Our theory relies on A<ref>, which holds as long as each component g_i(·,·,) of the true dynamics lies in the Hilbert space
  
  
      ℋ__i={ĝ_i,_̨i: ^n × ^m →  | ĝ_i,_̨i=_̨i^⊤_i,  _̨i∈^d}.
   
  Equivalently, A<ref> states that each g_i(·,·,) belongs to the RKHS with reproducing kernel 
  σ_i^2 _i(·,·)^⊤_i^-1_i(·,·). 
  Together, A<ref> and A<ref> state that with probability at least (1-δ_i), each component g_i(·,·,) belongs to a bounded subset of the  RKHS ℋ__i, defined as
  
  
  
  
  
  
  
  
  
  
      ℬ__i^ℋ = {ĝ_i,_̨i: ^n× ^m→  | 
      [            ĝ_i,_̨i=_̨i^⊤_i, 
         _̨i∈^d; _̨i  - _i,0__i,0^2 ≤ σ_i^2 χ^2_(1-δ_i) ]}.
  
  A<ref> and A<ref> are thus comparable to statements in related work on learning-based control with squared-exponential kernel GPs. For instance, <cit.> assume the true dynamics lie in a RKHS of known bounded norm, see also <cit.>.
  
  This observation suggests that we should meta-learn basis functions _i corresponding to a large space of models ℋ__i. 
  Since ℋ__i contains linear combinations of the elements of _i, A<ref> thus  motivates learning features that are orthogonal. To this end, we propose augmenting the meta-learning objective with the following regularizer:
  
      ℒ_⊥ reg() = α_⊥ reg∑_i=1^n ∑_ℓ=1^L+1I^ℓ-_i^ℓ⊤_i^ℓ_F^2,
  
  where α_⊥ reg > 0 controls the strength of this regularization, 
  the _i^ℓ are the weights of the neural network basis functions _i, and each I^ℓ is an identity matrix of appropriate dimensions. 
  
  As derived, theoretically analyzed, and empirically validated in <cit.>, this regularizer improves generalization.  
  
  Also, it incurs marginal computational overhead compared to the orthogonality regularizer for previously proposed in <cit.> to encourage interpretability of the model. 
  
  Beyond generalization, for our model with last layer adaptation, orthogonal features also lead to better statistical properties when adapting each _̨i through least-squares estimation with (<ref>), see <cit.>. 
  
  
  
  
  
  
  
  β-regularization: 
  Another aspect which can be regularized during offline meta-learning to improve downstream performance is the scaling factor β_i, which controls the size of the confidence sets for the model parameters _̨i. 
  With expressive neural network features _i, multiple settings of the neural network weights and of the prior parameters (_i,0, _i,0) can yield a model satisfying Assumptions <ref> and <ref>. As larger values of β_i yield more conservative confidence sets, it is preferable to select features and prior parameters that lead to lower values of β_i. To do so, we propose adding regularization to encourage meta-learning representations with lower values of β_i without compromising on safety.  
  
  
  
  
  We note from (<ref>) that 
  the value of each β_i depends on the ratio between the maximum and minimum eigenvalues of the prior and posterior precision matrices _i. 
  
  
  If λ̅(_i,0) ≤ 1  (which can be enforced by rescaling the features _i), then 
  
      λ̅(_i,0)/λ(_i,t)
      =
      λ̅(_i,t^-1)/λ(_i,0^-1) ≤λ̅(_i,t^-1)
      λ(_i,0^-1)
      ≤λ̅(_i,t^-1)
      λ̅(_i,0^-1),
   
  for any t∈ℕ. 
  Furthermore, λ̅()
      ≤√(Tr(^⊤)). 
  Combining with the above, we propose to regularize  an upper bound of the ratios λ̅(_i,0)/λ(_i,T_ML) during offline meta-training, namely: 
  
      ℒ_reg(_0) = α_reg∑_i=1^Tr(_i,T_ML^-T_i,T_ML^-1)
          Tr(_i,0^-T_i,0^-1),
  
  where α_reg>0 controls the strength of this regularization and is selected using a validation dataset and T_ML is the meta-training horizon. 
  As the model is directly parameterized by the inverse of the precision matrices _i <cit.>, this regularizer can easily be added to the standard training loss (<ref>).
  
  
  From (<ref>), we observe that β_i also depends on the ratio of determinants of the prior and posterior precision matrices (_i,t)/(_i,0). 
  Although a convex regularizer for this term can be derived, we found that including it  did not lead to performance improvements. 
  This ratio can be interpreted as capturing the amount of information that the model has gathered online, which is independent of the structure of the prior model. 
  Before learning, this ratio is 1, so the other term composed of the ratio of eigenvalues dominates β_i. 
  We observed that it is during these early stages that the meta-training model and its bounds β_i are most conservative, which could explain the importance of the regularizer in (<ref>), whereas regularizing the ratio of determinants appears to make little difference. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § SEQUENTIAL EXPLORATION AND EXPLOITATION FOR LEARNING SAFELY (SEELS)
  
  
  In order to ensure safety at all times with high probability and eventually reach , 
  we  provide a conservative deterministic reformulation of , 
  propose an algorithm tackling problems with high uncertainty, 
  analyze its safety and feasibility properties, and 
  discuss implementation.
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Relaxation of using reachability analysis
  
  
  In order to reason about how parameter uncertainty manifests in  terms of potential system behavior in the state space (and whether it might violate safety constraints), we translate confidence sets over parameters into confidence tubes over trajectories. 
  
  We stress that confidence tubes jointly depending on the parameters are required to satisfy the trajectory-wise chance constraint (<ref>), since enforcing a pointwise chance constraint at each timestep as in <cit.>  
  does not guarantee safety of the whole trajectory <cit.>. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We construct these tubes by leveraging the confidence sets _i,t^δ.
  
  Indeed, 
  a direct consequence of Theorem <ref> is that
  
      (
      ∃_̨i ∈_i,t^δ(_i,t,_i,t)
           ∀ t
          s.t._̨i^⊤_i  =  g_i(·,·,)
          ) 
           ≥ (1-2δ_i).
  
  This implies that a control trajectory respecting all constraints for all _̨i∈_i,t^δ and ϵ_i,t∈ℰ_i also satisfies (<ref>). 
  
  
  Thus, to conservatively reformulate (<ref>), we propose performing reachability analysis using the sets _i,t^δ and ℰ_i.
  
  Specifically, given a sequence of open-loop control inputs 
  = (_0,,_N-1), we define the sequence of reachable sets, for k=1,,N, as[This definition closely follows <cit.> for the specific case of a sequence of open-loop control inputs. One can account for a nominal feedback controller as a simple extension to reduce the tube size <cit.>. 
  In our simulated experiments, we omit feedback to better demonstrate the adaptation capabilities of the meta-training model, show the tightness of the confidence sets, and better verify the safety claims of the framework.]
  
      _k^t,δ()  = {_k =(·, _.7k-1, ,̨_.7k-1)
      	∘∘(_0, _0, ,̨_0)
      	
      | _0 = (t),  _̨i ∈ _i,t^δ,  ϵ_i,s ∈ ℰ_i,
       }
      ,
  
  where f_i(,,,̨_t)=h_i(, ) + _̨i^⊤_i(,) + ϵ_i,t. 
  
  
  
  
  
  
  
  
  
  By quantifying the uncertainty of the system dynamics and of the online learning process using confidence sets over the model parameters $̨ and by subsequently transforming these sets into reachable sets in the state space, we relax the original chance-constrained problem into a deterministic one:
  
  
  
  
  
  
  
  (N)
  
      min_, ∑_k=0^N-1(_k, _k) 
      s.t.[ ⋀_.8k=1^N _k^t,δ() ⊆ , 
       _N^t,δ() ⊆ ,;  ⋀_.8k=0^N-1_k ∈ ,  
       _0^t,δ = {(t)}, ]
  
  where=(_0,,_N)with _0=(t) are the centers of the reachable sets{_k^t,δ()}_k=1^Nin (<ref>), which depend on the control trajectory. 
  Since the cost functionℓtypically penalizes the deterministic open-loop control inputs, 
  we rewrite the expected cost in (<ref>) using a mean-equivalent reformulation as in <cit.>[As we use a sampling-based approach to compute the reachable sets <cit.>, the variance associated with these samples could be used to approximate the variance of the cost, or to minimize a given risk metric <cit.>.]. 
  This approach of considering a tube of reachable states and of evaluating the cost along a nominal trajectory is common in  the literature on robust MPC, see e.g., <cit.>. Note that the constraint-satisfaction guarantees of our approach do not depend on the chosen reformulation of (<ref>).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ : algorithm for safe learning and reaching 
  
  
  
  
  
  
  
  
  
  
  
  
  
  Due to high dynamics uncertainty, tight control constraints, and long planning horizons,  may be infeasible. 
  
  This motivates a safe learning-based exploration-exploitation framework to sequentially reduce uncertainty and eventually reach. 
  
  Our approach is based on a repeated two-phase approach: 
  
  when  is feasible, we enter the exploitation phase and plan a safe trajectory towith the current model uncertainty.
  
  In the exploration phase, we instead strictly perform safe exploration, planning an information-gathering trajectory that returns to the initial safe invariant set_0.  
  
  The corresponding robust optimal control problem is written as 
  
  
  (N)
  
      max_, ∑_k=0^N-1 (+)(_k, _k) 
      s.t.[  ⋀_.8k=1^N _k^t,δ ⊆ , 
       _N^t,δ ⊆ _0,; ⋀_.8k=0^N-1_k ∈ ,  
       _0^t,δ = {(t)}, ]
  
  where the sequence of reachable sets{_k^t,δ}_k=1^Nsatisfies (<ref>). 
  
  
  
  
  
   is similar to  but instead uses_0as the final set, thus ensuring the system is safe for the next phase with high probability. Also,  uses an information-gathering cost functionto encourage visiting states which reduce the remaining uncertainty in the dynamics 
      (, ) = 1/2∑_i=1^log(1+_i(,)^⊤_i,t^-1_i(,)).
  
  In the next section, we derivefrom the mutual information between the unknown dynamics and the observations. 
  
  
  
  
  
  
  
  
  SEELS: We summarize our approach in Algorithm <ref> () and show illustrative simulation results in Figure <ref>.  consists of iteratively learning a model of the dynamics by solving  before reachingwhen  admits a feasible solution. We perform a search over multiple horizons N_r and N_e, as the feasibility of each problem depends on the optimization  horizon (see Sections <ref> and <ref>)). 
  
  
  
  This split yields a tractable sequence of trajectory optimization problems, although it induces sub-optimality relative to the computationally intractable problem of simultaneous exploration and exploitation <cit.>. 
  
  Since safely learning nonlinear dynamics with minimal assumptions remains an open area of research (see Section <ref>), quantifying this sub-optimality gap (perhaps by computing the regret of the algorithm) is beyond the scope of this work.
  
  
  
  
  
   §.§ Information objective
  
  
  
  
  During the exploration phase, we perform trajectory optimization with an objective functionthat encourages visiting states and taking actions that reduce uncertainty over the unknown dynamics. To do so, a natural objective to maximize is the mutual information between the unknown dynamicsand the observations of the state. 
  This cost characterizes the information gain<cit.> about, gained from applying the control inputsto the true system from a stateand observing the next state corrupted by random disturbances. 
  
  
  
  
  
  
  To derive this objective, we leverage the Bayesian interpretation of  which is used during offline meta-training to obtain a model that has the correct structure (Assumption <ref>) and is calibrated (Assumption <ref>).  
  
  
  
  Specifically,  is trained assuming zero-mean Gaussian-distributed observation noiseϵ_i,tof varianceσ_i^2and a Gaussian prior over model parameters_̨i. 
  Using (<ref>) andttransitions from the system, this yields the posterior distribution_̨i∼(_i, σ_i^2_i,t^-1). 
  
  In this setting, the marginal distribution over observationsx̃^+_i =  x_i,t+1-h_i(_t,_t)=_̨i^⊤_i(_t,_t)+ϵ_i,tgiven an arbitrary state_tand control input_tis a normal distribution(_i _i, (1+_i^⊤^-1_i,t _i)σ_i^2), where_i = _i(_t,_t). 
  
  
  
  For this formulation, the mutual informationℐbetween the observation^+ =  _t+1 - (_t,_t)and the true unknown dynamicscan be derived in closed-form. 
  It is defined using the entropyℋ(·), which for a random variable^+∼𝒩(, )evaluates toℋ() = (1/2)log(det(2πe ) ). Hence, the information gain from observing the scalar random variablex^+_ican be expressed asℐ(x^+_i ; g_i) = ℋ(x^+_i) - ℋ(x^+_i | g_i)
  = 1/2( 
  log(var(x^+_i)) -log(var(x^+_i | g_i)
  )
  = 1/2( 
  log((1+_i^⊤_i,t^-1_i)σ_i^2)) -log(σ_i^2))
  )
  = 
  1/2( 
  log(1+_i^⊤_i,t^-1_i)).  
  
  This quantity expresses the information gain from observing each dimensioniof the state, which are modeled independently in our formulation. By summingℐ(x^+_i ; g_i)overi=1,,n, we arrive at the exploration objectivein (<ref>). 
  
  
  
  The objectiveis a function of the current information state of the model, specified by the updated precision matrices_1,t, , _,t. It explicitly encourages taking actions and visiting states in the feature space spanned bywhich have highest variance_i,t^-1. The resulting observations are then the most informative in terms of reducing uncertainty over.
  
  
  Note that the expected information gain along a trajectory is not simply the sum of the expected information gains per transition, as expressed in  when summing (<ref>) overk=0,,N-1. 
  However, correctly computing the expected information gain along the trajectory would require factoring in model updates along the trajectory; we find that considering the sum of single-transition information gain with the current precision matrices_i,tis sufficient in guiding exploration. 
  
  Also, while in this work we assume bounded (non-Gaussian) noise corrupting our measurements, 
  we find that using this information objective 
  works well in practice to encourage exploration. 
  
  The problem of optimal exploration is beyond the scope of this framework.
  
  
  
  Finally, we note that this cost functiondoes not suffer from computational complexity scaling with the amount of data, as is the case for similar objectives derived for squared-exponential kernel GPs <cit.>. This computational efficiency comes from choosing finite-dimensional basis functions to parameterize the model class for the true dynamics.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Probabilistic safety and feasibility guarantees
  
  
  
  Next, we prove that 
   enjoys probabilistically guaranteed safety and ensures that  remains feasible at all times with high probability, which is crucial to enable autonomous operation. Our next result also states that provided that  is eventually feasible, the agent will safely reachwith high probability. This result characterizes our proposed approach as a feasible solution to .
  
  
  
  
  With the confidence sets defined in (<ref>) and the reachable sets satisfying (<ref>), apply Algorithm <ref> 
  to sequentially solve  and  and to control the uncertain system (<ref>). 
  Then, under Assumptions <ref>-<ref>:
    * With probability at least (1-δ), there exists an optimization horizon N_e ensuring the feasibility of (N_e) in each 
  exploration phase in line 9 of Algorithm <ref>.
  
    * With probability at least (1-δ), during exploration, the state trajectory satisfies all constraints _t∈ at all times t and finishes in _0 at the end of each phase.
    * With probability at least (1-δ), assuming that  is eventually feasible,  
  the system is guaranteed to satisfy (<ref>), i.e., to be safe at all times and eventually reach . 
   
  First, we prove that  is recursively feasible with high probability, before proving that the system satisfies all constraints during exploration and that  (<ref>) is satisfied if  is eventually feasible.
  Probabilistic recursive feasibility: Let  be the number of exploration phases before  becomes feasible[ is a random variable that can be defined as a stopping time with respect to the natural filtration associated to the stochastic process _t. This result also holds if  is not feasible and the algorithm can never solve   (e.g., if  is surrounded by obstacles). Indeed, if the algorithm is stuck in an infinite number of exploration steps, the last inequality of this proof still holds for →∞, by Theorem <ref>.]. 
  
  
  Also, let ^j, and t_j =∑_l=1^j-1^l be, respectively, the planning horizon and the start time index of each _j, where j=1,,.  
  
  
  For conciseness, define _j for {_j is feasible}, corresponding to the event that the j-th exploration problem is feasible. 
  
  
  Then, 
  
  
      (
      ⋀_j=1^_j
      )
         ≥(
      ⋀_j=1^_j,  _t_∈_0,  )
      
         
      =
      (
      _ | ⋀_j=1^-1_j, _t_∈ _0
      )
      
         ·(
      ⋀_j=1^-1_j, _t_∈ _0
      )
      .
  
  
  By Assumption <ref>, given that _t_j∈_0, _j is feasible for any j-th exploration phase.  
  Indeed, choose ^j=1 for _j. 
  Then, _0^j=π(_t_j) is a feasible solution to _j. 
  Thus, the event {_j  | _t_j∈_0} holds with probability one.  
  
  In particular, this implies that the event  {_ | ⋀_j=1^-1_j,
  _t_∈ _0} holds with probability one.  
  
  
  
  Next, we 
  leverage our confidence sets:
  
      (
      ⋀_j=1^_j
      )
         ≥(
      ⋀_j=1^-1_j, _t_∈ _0
      )
      
         ≥(
      ⋀_j=1^-1_j, _t_∈ _0,  ∈_i,t_-1^δ ∀ i
      )
      
         
      =
      (
      _t_∈_0  | ⋀_j _j, ∈_i,t_-1^δ∀ i
      ) 
         ·(
      ⋀_j _j, ∈_i,t_-1^δ∀ i
      ).
  
  By construction of the reachable sets {_k^t_j,δ}_k=1^^j, by Assumption <ref> and by definition of _j (since _^j^t_j,δ⊆_0), we have that _t_j+1∈_0 given that _j is feasible and  ∈_i,t_j^δ ∀ i, for any j-th exploration problem. 
  
  Thus, the first term {_t_∈_0  | ⋀_j _j, ∈_i,t_-1^δ∀ i
  } holds with probability one, and 
  
      (
      ⋀_j=1^_j
      )
         ≥(
      ⋀_j=1^-1_j, ∈_i,t_-1^δ ∀ i
      ).
  
  Since _0 is feasible with probability one as _0∈_0, 
  and by reasoning by induction for all j=,,0, we obtain that 
  
      (
      ⋀_j=1^_j
      )
         ≥(
      ⋀_j=1^-1∈_i,t_j^δ ∀ i
      )
      ≥
      (1-δ),
  
  where the last inequality comes from Theorem <ref>. 
  
  
  
  Probabilistic safety: Before proving the second statement of Theorem <ref>, we prove the third statement by assuming that <∞, which holds if  is eventually feasible. Let also N_reach, and t_f be, respectively, the planning horizon and the start time index of .  
  
  For conciseness, define _k^t_j=_t_j+k, corresponding to the state at time (t_j+k) in the j-th phase. 
  
  Further, define the event that the trajectory during the j-th exploration phase (or exploitation phase) satisfies all constraints as
  
  
      {_info^j  ∈ _info^j}    = {⋀_k=1^^j(_k^t_j ∈ ) 
      ∧(_^j^t_j ∈ _0)
      }
      ,
      
      
      {_reach ∈ _reach}    = {⋀_k=1^N_reach(_k^t_f ∈ ) 
      ∧(_N_reach^t_f ∈ )
      }
      ,
  
  where j = 1,,. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  With this notation, we rewrite the probabilistic safety constraint of the original problem as
  
      (<ref>)   =
      0.9(
      
      j=1⋀{_info^j  ∈ _info^j}∧{_reach ∈ _reach})≜(
      {Success})
      
         ≥(
      {Success}  |  ∈ _i,t^δ ∀ t  ∀ i)
      ( ∈ _i,t^δ ∀ t  ∀ i)
      
      ,
  
  where t=t_1,,t_,t_f, and i=1,,n.[The inequality follows from (<ref>)=0.8(
  {Success}  |  ∈ _i,t^δ ∀ t∀ i)
  ( ∈ _i,t^δ ∀ t∀ i)0.8+  (
  {Success}  | ∉_i,t^δ ∀ t, ∀ i)
  (∉_i,t^δ ∀ t, ∀ i)≥ (<ref>).]
  
  
  
  Again, by Assumption <ref>, the meta-learning model can fit the true dynamics. Hence, if the true parameters are within the confidence sets _i,t^δ, then, the reachable sets _k^t,δ defined in (<ref>) necessarily contain the state trajectory of the true system. Hence,
  
      {_k^t,δ⊂}
      =
      {_k(^̨*)∈  | ∈_i,t^δ,  ∀ i}
      .
  
  By definition of  and , the reachable sets are subsets of the safe set (with probability one), and 
  
      (_k^t(^̨*) ∈ ,     k=1,,N   | ∈_i,t^δ,  i=1,,)
      
         =
      (_k^t,δ ⊂ ,  k=1,,N)
      = 1,
  
  which also holds for the final constraints _N^t∈_0 and _N^t_f∈. 
  
  Thus, 
  
   (
  {Success}  | ∈_i,t^δ ∀ t  ∀ i) = 1.
  
  Combining this result with (<ref>), we obtain that the system satisfies all constraints and eventually reaches  with probability at least that of the probability of the model parameters belonging to the confidence sets, i.e., 
  
  (<ref>)≥(∈_i,t^δ ∀ t  ∀ i).
  
  
  This last term holds with probability greater than  (1-δ). Indeed, using (a) a union bound (Boole's inequality) and (b) Theorem <ref>, we obtain
  
      (<ref>) ≥ ( ∈ _i,t^δ ∀ t ∀ i) 
      
      
      
      
      
      
      
       = 
      1 -
      (⋁_i=1^n⋁_t∉_i,t^δ) 
      
         ≥^(a)
      1  - ∑_i=1^n(⋁_t ∉ _i,t^δ) 
       =  
      1 -∑_i=1^n(
      1 -(⋀_t ∈ _i,t^δ) 
      )
      
         ≥^(b)
      1 -
      ∑_i=1^n(
      1 - 
      (1-2δ_i)
      )
      =
      1 -
      ∑_i=1^n(
      2δ_i
      )
      =
      (1-δ)
  
  since δ_i = δ/(2n). This concludes the third statement  of Theorem <ref>. 
  
  The second statement of Theorem <ref> follows by relaxing the assumption that <∞ and by redefining the event {Success}
  =
  ⋀_j=1^{_info^j  ∈ _info^j}. In this case, all inequalities follow similarly from Theorem <ref>.
  
  
  
  
  
  Key to this proof are Theorem <ref> and the definition of the reachable sets in (<ref>). 
  We stress that our safety probability(1 - δ)is independent of the (unknown) time to reach, which would not be the case if  pointwise chance constraints were used instead of (<ref>).
  
  
  
  To solve  by satisfying all constraints and reaching  with high probability, Theorem <ref> relies on  eventually becoming feasible. If the original problem is feasible with perfect knowledge of the dynamics, this assumption holds if the objective used for exploration leads to actions that continually reduce uncertainty in dynamics, see <cit.>. Related conditions on observability and persistence of excitation are available in the literature, see<cit.>. With full state information and long exploration horizons using our information cost with Algorithm <ref>, we observed in our experiments that  eventually becomes feasible if the original problem is truly feasible with perfect knowledge of the dynamics, see Sections <ref> and <ref>.
  
  In situations where the problem is not feasible in the first place (e.g., an obstacle blocks the only path to), our algorithm defaults to exploring within a feasible neighborhood around_0, in which case the system is guaranteed to satisfy all safety constraints at all times with probability greater than(1-δ). 
  
  
  By assuming bounded disturbances_tand exploiting confidence sets over the model parameters which hold jointly for all times with high probability (Theorem <ref>), we can guarantee the feasibility of  during exploration. This contrasts with related work in the MPC literature which provides probabilistic recursive feasibility over a finite horizon only  <cit.>. 
  
  This aspect is crucial to enable autonomous operation and reliability of the approach.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Implementation and practical considerations
  
  Implementation of the algorithm is complicated by challenges in reachability analysis and non-convex optimization. 
  
  Reachability Analysis: Computing the reachable sets in (<ref>) is difficult due to the non-convexity of the features. Methods reasoning about single-step set propagation (e.g., <cit.>) are generally too conservative <cit.> as they neglect time correlations induced by the parameters_̨i. Moreover, the online updates to the model parameters preclude exact methods which perform computations offline to compute reachable sets, e.g., <cit.>. 
  
  Neural network verification techniques <cit.> are tailored to offline computations and are not yet fast enough for online trajectory optimization. 
  
  
  
  Finally, methods using Lipschitz continuity can provide conservative reachable set approximations <cit.>, but such techniques would be so conservative for the systems we consider in this work that the agent would never deem reaching the goal to be feasible, even if one had access to the true Lipschitz constant of the system. The main limitation of <cit.> is that it adopts a single-step worst case analysis and does not account for the time correlations of the parameters <cit.>.  
  
  For these reasons, we leverage a recently-derived sampling-based uncertainty propagation scheme for reachability analysis () <cit.>. By sampling M parameter  and disturbance tuples(_̨i^j,ϵ_i,k^j, 0.851≤ i≤,   0≤ k≤ N-1)within their associated bounded sets _i,t^δ and Ε_i, 
  computing the resulting reachable states _k^j for these tuples according to x_i,k+1^j=h_i(_k^j,_k^j)+_̨i^j⊤_i(_k^j,_k^j)+ϵ_i,k^j, and approximating the reachable sets _k^t,δ in (<ref>)with the convex hull of the samples {_k^j}_j=1^M, provides a scalable approach to compute these tubes with minimal assumptions on the system's dynamics, enabling the use of arbitrary neural network features. 
  It is fast enough to be used with common robotic systems (computing a trajectory requires less than a second for a 13-dimensional spacecraft system with a Python implementation <cit.>) and can be further accelerated through parallelization on GPUs. 
  Although  lacks finite-sample guarantees of safety (see also <cit.> for a recent extension with finite-sample guarantees), 
  asymptotic guarantees can be derived using random set theory <cit.>, and finite-sample approximations are sufficient to ensure safety in practice, as shown in our results. 
  
  
  
  
  Optimization-based planning:
  Using , we reframe a generally intractable stochastic optimal control problem into  and , which are non-convex optimal control problems. Efficiently computing solutions to this class of problems is an active field of research. In this work, we use a direct method based on sequential convex programming (SCP) <cit.>. 
  By solving a sequence of convexified versions of the original problem, SCP-based methods can run in real time and provide theoretical guarantees of local optimality <cit.>. 
  In this work, we initialize each method with an infeasible straight-line trajectory, solve each convexified problem using OSQP<cit.>, and provide further details in Appendix <ref>. 
  
  
  Additionally, due to uncertainty, the feasibility of each problem depends on the optimization horizonN. 
  For this reason, Algorithm <ref> performs a search over a predefined range of planning horizons (lines2and8). 
  For exploitation, it selects the first feasible solution if one exists, although other criteria could be used, e.g., minimal control cost. 
  For exploration, we select the trajectory which leads to the largest expected information gain. 
  Indeed, due to tight control constraints and safety constraints, a longer horizon does not necessarily lead to higher information gain. 
  This heuristic works well in practice, and future work will adopt a continuous-time problem formulation with free final time, which is an active area of research  <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § SIMULATION STUDIES
  
  
   §.§ Robotic manipulation in a cluttered environment
  
  
  First, consider the three-link open-chain torque-controlled manipulator shown in Figure <ref><cit.>. The state of the system is = (q_1,q_2,q_3,q̇_1,q̇_2,q̇_3)and its control is=(τ_1,τ_2,τ_3), whereq_iis the angle of thei-th joint andτ_iis its torque. Its true dynamics follow
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
      M()+C(,)+N()=,
  
  where 0.95= (q_1,q_2,q_3). 
  
  The expressions forM(),C(,), andN()depend on the lengths 0.95l_1  =  1 m and 0.95l_2 = 
  l_3 = 0.5 m, massesm_i, and inertiasI_iof eachi-th link of the manipulator (see <cit.>). By introducing uncertainty in eachm_iandI_i, this experiment corresponds to a scenario where the robot must manipulate uncertain payloads and identify its inertial properties. 
  
  Although manipulators are typically controlled using trajectory generation and high-frequency low-level tracking (rather than open-loop control), 
  the instability of this system makes it a useful case study for our framework.  
  
  
  
  
  Starting at(0) = 0 ∈ _0
   = 0.9{  |  
  |q_1|≤π/8,
    |q_2,3|≤π/6, 
    |q̇_1,2,3|≤π/6}, 
  we consider the problem of 
  
  moving the arm to the goal region = 0.9{  |  
  |q_1 -π/2|≤0.9π/8,
    |q_2,3|≤π/10, 
    |q̇_1,2,3|≤π/8}. 
  
  We consider two spherical obstacles_i ⊂ ^3which should be avoided by the end-effector. 
  This results in two obstacle avoidance constraints∉_iwritten as_i = 0.95{  |  _e()-_i_2 ≤0.951/4}, 
   where 
   0.9_1 = (√(2)/2,-√(2)/2,1.4),  
  0.9_2 = (√(2)/2,-√(2)/2,0.6), and_e:^6 → ^3maps the stateto the end-effector position 
   <cit.>. 
  
  We also account for joint angle and velocity limits ∈ so that 
  0.95= ∖(∪_i_i), with= 0.9{  |  
  q_1∈[-π/2,π], |q_2,3|≤π/2,
    
  |q̇_i| ≤ π/8}. 
  
  Control bounds are set to𝒰 = 0.95{ | |τ_i| ≤  2  N·m}. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We simulate (<ref>) with an Euler discretization scheme with timestep50 ms. 
  
  To slow down the unstable dynamics of this system, we use a damping controller with approximate gravity compensation = N̂(,) - K_d + _seels, whereN̂is computed using nominal mass and inertia values,K_d = [0,4,1], and_seelsis the open-loop control trajectory computed with our framework, where each input lasts1 sand_seels ∈ 𝒰.
  
  We use a nominal modelof (<ref>) with 0.95C(q,q̇) = 0 and  0.95(m_i,I_xi = I_yi,I_zi) = (2,0.05,0.01). 
  
  We then train an  model withL=2layers of width128each and an output dimensiond=64for each_i. 
  
  We randomize parameters according to uniform distributions withm_i ∈ [1.75,2.25] kg,I_xi = I_yi∈[0.04,0.06] kg·m^2,I_zi ∈ [0.008,0.012] kg·m^2, with 0.9i = 1,2,3. We assumeσ_i = 10^-3, and|ϵ_i,t|≤(σ_i^2χ_1^2(0.95))^1/2, 0.9i = 1,,6. 
  
  We useβ_iand orthogonality regularization,M = 2500samples for reachability analysis 
  
  and directly use our theoretically computed bounds within our safe learning algorithm, i.e., we use (<ref>) to sample model parameters. 
  
  Further details about the trajectory optimization are provided in Appendix <ref>.
  
  
  
  
  
  
  
  
  
  Figure <ref> (left) shows the results for low mass and inertia 0.9(m_i,I_xi = I_yi,I_zi) = (1.85,45·10^-3,8.5·10^-3). Note that lower mass and inertia would cause the robot to leave_0in1s, making the problem infeasible. In this experiment, the end-effector initially rises (due to the lower-than-expected mass/inertia), does a second exploration maneuver, then reaches the goal, always satisfying constraints. 
  
  Figure <ref> (right) shows results for a heavier arm. 
  
  Again, all constraints are satisfied and the goal is reached safely.  
  
  For both scenarios, we compare to a method which only considers uncertainty from the additive disturbances_t, thereby deciding that directly reachingis safe. This naive approach fails to reach the goal and violates constraints, which demonstrates the need for an accurate model
  and sequential online learning to reliably solve this problem. 
  
  
  
  
  
  
  We plot control trajectories and angular velocities in 
  Figure <ref>. This shows that steering the system with lower mass values requires less control effort, 
  and 
   angular velocities(q̇_2,q̇_3)are larger during the second exploration phase (fromt = 1s untilt = 5s) in order to efficiently identify the unknown parameters of the manipulator. 
   
   Also, note the approximate symmetry between the two control trajectories due to the different mass and inertial properties.
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Safely transporting an uncertain payload
  
  
  
  Next, we verify our approach on a nonlinear six-dimensional planar free-flyer robot navigating in a cluttered environment. 
  
  We consider the problem of cargo transport, in which the robot is attached to an uncertain payload that results in changes to the inertial properties of the system. This system mimics a cargo-unloading scenario that is one plausible near-term application of autonomous robots onboard the International Space Station <cit.>. 
  
  
  
  
  
  
  The state of the system is given by= (,θ,,ω)∈^6, with,∈^2the planar position and velocity andθ,ω∈the heading and angular velocity, respectively. 
  For safety, we constrain|v_i| ≤  0.2m/s, and|ω| ≤  0.25rad/s. 
  
  The robot is controlled using two pairs of gas thrusters and a reaction wheel. The system's control inputs are denoted by = (F_x,F_y,M)∈^3, where=[-u̅_i,u̅_i]^×3represent the limited control authority, withu̅_1,2 = 0.15N andu̅_3 = 0.01Nm. 
  
  The payload causes a change in mass/inertia properties and causes the center of mass to be offset at_0∈^2.  
  
  
  
  
  The continuous-time nonlinear dynamics of the system (which we write as=𝐟_t(·)) 
  are given by  =  , =  ω, and  
  
      m = -ω̇[ -p_oy;  p_ox ]+ω^2 _o,
        
      Jω̇ = 
      M 
      - 
      p_oxF_y 
      +
      p_oyF_x
      .
  
  
  
  To capture payload uncertainty, we randomize the massm, inertiaJ, and center of mass_0of the coupled system as 
  
  
      m ∼Unif(25,60)  kg, 
        
      J ∼Unif(0.30, 0.70)  kg·m^2, 
      
      
      p_oi∼Unif(-7.5,7.5) cm,  i ∈ {x,y}
      .
  
  
  Using a zero-order hold on the controls and a forward Euler discretization scheme withΔt = 3   s, we discretize (<ref>) as
  
      _t+1 =  _t + Δ t ·𝐟_t(_t,_t, m, 𝐉,_o)
      + _t
      .
  
  We assume that theϵ_i,tareσ_i-subgaussian, each bounded as 
  0.95|ϵ_i,t| ≤ (σ_i^2χ_1^2(0.95))^1/2, where 
  0.95σ_1,2^2=10^-6,  
  0.95σ_3,6^2=10^-5, and 
  0.95σ_4,5^2=10^-7. 
  We use this discrete-time system in simulation experiments and to collect training data for offline meta-learning. 
  
  
  
  We use a nominal modelof the system using (<ref>) with(m̅,J̅,_0)=(35,0.4,0), which corresponds to a double-integrator model. 
  To represent the unknown model mismatch(·,·,), 
  we train an  model as described in <cit.> for6000iterations for all experiments, 
  
  withL=2layers of width64each and an output dimensiond=32for each_i. 
  
  
  
  
  
  
  
  For trajectory optimization, we use standard linear-quadratic final and stagewise costs on states and controls to minimize control cost and deviation to_0ordepending on the phase. 
  Specifically, we maximize the information cost  (<ref>) while minimizing control effort and   penalizing high velocities and the final distance to_g, the center of either_0, or: 
  
  
      min_,∑_k=0^N-1(
      -α_info(_k,_k)
         +
      _k^⊤_k + _k^⊤_k
      )
      
         + 
      (_N - _g)^⊤_N(_N - _g)
      .
  
  We setα_info = 0.025for exploration, whereasα_info = 0when reaching. 
  We set = 0.95diag([0,0,0,1,1,10]), = 0.95diag([10,10,10]), and_N = 0.9510^3diag([1,1,0.1,10,10,10])for both  and . 
  
  
  
  
  
  To validate our framework, 
  we run  on 250 randomized problems with different dynamics parameters, 
  four different obstacle configurations and initial/final conditions, as shown in Figure <ref>. 
  
  WithM=2500and both orthogonality andβ-regularization, we obtain a success rate (all constraints in (<ref>) are satisfied, the optimizer finds a feasible solution at each step of , and the system reaches) at93.2%forδ = 0.1(withδ = 0.1, the success rate should be at least90%according to Theorem <ref>), at90.5%forδ = 0.2, and at88.8%forδ = 0.5. This experiment shows that the trajectory-wise chance constraint (<ref>) is conservatively satisfied in practice and verifies the probabilistic safety and recursive feasibility results of Theorem  <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Sensitivity analysis: 
  we perform an ablation study to evaluate the sensitivity of our approach to different parameters: 
  toδ, 
  to the effect of theβ-regularizer, 
  to the number of samplesMfor reachability analysis with , 
  and 
  to the noise intensitiesσ_i.   
  
  For each parameter set, we perform 250 randomized experiments on the environments shown in Figure  <ref>. 
  
  As with the manipulator experiments, we compare our approach to a baseline.   
  
  Results are shown in Figure <ref>, in Tables <ref> and <ref> and are discussed below:
  
  
  
  
  
  
  
    * Mean-equivalent baseline: we compare to a method which only considers uncertainty from the additive disturbances _t and uses the mean parameters _i for planning, with meta-learned  features _i. This baseline thus assumes reaching  directly is safe. This naive approach violates the trajectory-wise chance constraint (<ref>), again demonstrating the need for uncertainty quantification and sequential active learning to reliably solve this problem.
      
      
    	
    * Random Fourier features (RFF) baseline: we replace  with a Bayesian model with d=128 random Fourier features <cit.> as described in <cit.>. 
    	  Offline, prior to running Algorithm <ref>, we select length-scale hyper-parameters and calibrate the priors (_i,0, σ_i^2 _i,0) over the weights to obtain an accurate uncertainty representation satisfying Assumption <ref>. 
    	    We observe that this approach does not achieve the desired success rates, which we attribute to substantial violations of Assumption <ref>. Making accurate multi-step open-loop  predictions appears challenging without expressive meta-learned neural network features.
    * δ: 
  	as a consequence of Theorem <ref>, the algorithm's conservatism can be tuned by the
  choice of δ, which is supported by our results. 
  In particular, by opting for a lower probability of safety 1-δ, the goal is reached faster on average. 
  
  Importantly,  conservatively satisfies the chance constraint (<ref>) for these different parameters. 
      
      
      
      
      
    * M: increasing M leads to increased success rate and safety probability, which is supported by the results shown in Figure <ref> (right). 
      
      Theoretically, under A<ref>-<ref>,  
  success with high probability (as stated in Theorem  <ref>) is guaranteed as long as the number of samples M is large enough (a consequence of <cit.>), which is supported by our experiments. 
  
  In practice, one should choose the largest value of M to satisfy a given computational budget and multiple techniques could be used to improve the quality of the approximation, e.g., parallelization on GPUs and importance-sampling approaches <cit.>, which will be considered in future work. We refer the reader to <cit.> and <cit.> for further details. 
  
  
  
      
    * β-regularization 
      reduces conservatism while still guaranteeing probabilistic safety in practice, i.e. the model still satisfies A<ref> and A<ref>. 
      
      Indeed, Figure <ref> shows that the goal is reached faster on average and the system still conservatively satisfies (<ref>). 
      
      This also implies that the model class with last-layer adaptation of  is expressive enough to enable meta-training of models fulfilling A<ref> and A<ref> with varying levels of conservatism, which could motivate future research on better meta-training of such models with different architectures and regularizers.
      
      
      
      
      
    * σ_i: we also consider two noise levels:
  
    * σ_i^2=10^-7 for  , and 
  σ_i^2=10^-6 for  .
  
    * σ_i^2=10^-6 for  ,  
  σ_i^2=10^-5 for  , and 
  σ_i^2=10^-7 for  .
  
  
  Results for these different noise levels for different δ are reported in Tables <ref> and  <ref>. 
  From Table <ref>, we see that the performance and overall probability of safety for small σ_i is not sensitive to the chosen value of δ. We conjecture that failures are mostly due to under-approximations from the approximate computation of the reachable sets with . For higher noise levels, it is evident that the conservatism of the algorithm can be tuned by choosing a different value for δ, since failures come from statistical errors from updating the model with noisy data (see Theorem <ref>). We also observe faster times to reach  when opting for lower probability of safety. In all scenarios,  solves the problem with probability at least (1 - δ), verifying the results of Theorem <ref>.
  
  
  § HARDWARE EXPERIMENTS
  
  
  Next, we verify our framework via hardware experiments. 
  
  
  
  We consider the planar spacecraft free-flyer platform controlled by eight thrusters as shown in Figure <ref>. 
  
  
  
  
  
  
  
  
  
  
  As the system weighs∼ 16kgand each thruster can only apply a maximum force of0.4N, an accurate model and dynamically-feasible trajectory are necessary to safely reach a goal region.  
  
  Other challenges include imperfect actuators,  time-varying mass of the system as the gas tanks deplete, and unmodeled friction and tilt of the table surface.  
  
  Our implementation follows Section <ref> (withδ = 0.1,  orthogonality andβ-regularization) with the following modifications:
  
    * Dynamics: we meta-train using (<ref>) discretized at 10Hz with an additional random force 𝐅_tilt to account for the tilt of the table. For offline meta-training, the parameters are randomized as 
  
      m ∼Unif(8,40)  kg, 
        
      J ∼Unif(0.08,0.30)  kg·m^2, 
      
      
      p_oi∼Unif(± 5) cm,
        
      F_tilt,i∼Unif(± 0.3) N, 
        
      	i ∈ {x,y}
      .
  
    * Regulation: we regulate the system with a PD controller  at the end of each phase as a new plan is computed.[We solve four  with different horizons N (, line 7) which takes a total of 5 seconds in average with our Python implementation, measured on a computer with a 3.70GHz Intel Core i7-8700K CPU. As the computed control inputs are directly used to steer the system, this delay may be problematic if the robot is in movement as a solution is computed.]  We do not use these intermediate trajectories to learn dynamics. 
      
    * Feedback: we use an LQR controller to stabilize the trajectory around reachable set centers _k and reduce uncertainty. This feedback controller is computed using the prior mean parameters _i,0 of , linearized around (0).   We account for this controller during planning as in <cit.> and apply the control _k = _seels,k + LQR(_k - _k) at 1/3Hz.
  
  
  
  We use a motion capture system and interface our Python implementation with ROS.  
  
  Using , we compute a desired wrench = [F_x,F_y,M]in the world frame, which we then convert to the local frame of the robot, yielding_body. 
  We then translate these forces and torque into corresponding PWM signals for the actuators:_PWM =
  saturate(A_body)∈^8.The constant matrixA ∈ ^8×3accounts for the maximum applicable force of each thruster of0.4N. With the control bounds for planning set conservatively to|F_x,y|≤0.4N and|M|≤0.03Nm, we can ensure that actuator saturation will neither occur during exploration nor exploitation.  
  
  
  This approach of using a lower-level control mapping to turn world-frame forces and torques to gas-thruster commands allows us to use the same dynamics model used in the simulation experiments (Section <ref>) for the hardware experiments (retraining it to also account for the slight tilt of the table).
  An alternative would be to learn a dynamics model in terms of gas-thruster inputs directly. While this modeling choice would allow factoring in dynamics uncertainty at the actuator level (e.g., actuator failure), this would also yield a more complex class of possible models and thus require a larger model to ensure the assumption A<ref> on model capacity is still satisfied.
  
  
  
  With these modifications, we apply  on multiple obstacle fields and analyze the representative experiment in Figure <ref>. 
  Initially, uncertainty is too large to safely reachand  is infeasible. Instead, the agent performs a short exploration trajectory with negligible control inputs to first infer the tilt of the table ( is necessarily feasible for a planning horizon short enough, see Theorem <ref>).  
  Indeed, applying larger control inputs would yield larger reachable sets (due to inertial uncertainty) and violate constraints. 
  Then,  still deems reaching to be unsafe and thus performs two additional exploration trajectories to infer the inertial properties of the system. Once  is feasible, the robot finally reaches the goal. All constraints are satisfied throughout the experiment.
  
  
  
  We present further results in the video available at 0.95https://youtu.be/5ux1Jg4FgEshttps://youtu.be/5ux1Jg4FgEs 
  and summarize them below:
  
    * Baseline: we compare with a mean-equivalent approach which collides with obstacles despite the LQR controller. This shows that considering uncertainty during planning and learning the dynamics is necessary to reliably solve this task.
      
    * Reliability: our approach works repeatably in multiple obstacle fields, including a payload-grasping experiment. 
      We also present two experiments where the agent infers its dynamics, reaches a first goal, and is then capable of reaching another destination directly without needing further exploration. All constraints are satisfied throughout the experiments (Theorem <ref> can easily be extended to also guarantee the success of such multi-setpoint problems).
      
    * Failure modes are also discussed.  
      We observed that failures to reach  are due to one of the following:
      (1) the model adapting incorrectly during the initial exploration phases due to noisy measurements and sim-to-real mismatch, resulting in large variations in the parameters _̨i, which is easy to detect and could be addressed by re-initializing parameters based on the variances σ_i_i^-1, using event-triggered learning <cit.>, or by increasing σ_i to learn less rapidly, 
      (2) infeasible trajectory optimization, or  
      (3) gas tanks that are too depleted to apply desired forces. 
      
      In these situations, we observed that no feasible trajectory to the goal is found and the algorithm defaults to exploring safely within _0. 
  
  
  § CONCLUSION
  
  
  To safely perform tasks under high initial uncertainty, 
  we presented a framework to sequentially explore and learn the properties of a dynamical system while guaranteeing safety at all times jointly as a single trajectory-wise chance constraint. 
  We demonstrated 
  our approach through simulations 
  
  
  and hardware experiments, where we showed that a free-flying robot emulating a spacecraft can reliably infer its dynamics and safely reach a goal. 
  
  Future work. 
  
  Problems with time-varying dynamics (e.g., the system is attached to a new payload during operation and gas tanks are depleting) could be addressed using continual meta-learning approaches <cit.>, which would require a new analysis to guarantee safety.
  
  Future work could also 
  
  
  derive algorithms to compute invariant sets satisfying A<ref> by leveraging the prior on the meta-learning model. 
  
  Extensions of <cit.> would strengthen the safety guarantees to the case of a finite number of samples for uncertainty propagation and 
  
  work on free-final-time trajectory optimization would replace the search over horizonsNin our algorithm. 
  
  Investigating regret bounds for such constrained problems would also help provide guarantees on the time required to perform a task given known geometric properties of the problem.   
  
  Finally, our framework could be combined with existing learning-based controllers to improve performance.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § ACKNOWLEDGMENTS
  The authors thank Edward Schmerling and the anonymous reviewers for their feedback which significantly helped improve the presentation of this work. This work was supported in part by NASA under the University Leadership Initiative (Grant #80NSSC20M0163) and the Early Stage Innovations program, and by DARPA under the Assured Autonomy program. This article solely reflects the opinions and conclusions of its authors. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  IEEEtran
  
   §.§ Proof of Theorem <ref>
  
  
  The proof of Theorem <ref> follows from the proof of <cit.>, by making substitutions accordingly for our meta-learning model. To do so, we use the following lemma, which corresponds to <cit.>. We restate this result below by changing notations to match our problem setting, where the index i=1,, corresponds to each state component.<cit.> (Self-Normalized Bound for Vector-Valued Martingales).
  Let {ℱ_t }_t=0^∞ be a filtration and i∈ℕ. 
  Let{ϵ_i,t}_t=1^∞be a  real-valued stochastic process such that ϵ_i,t is ℱ_t-measurable and conditionally σ_i-subgaussian, such that [exp(λϵ_i,t) | ℱ_t-1]≤exp(λ^2σ_i^2/2) for all λ∈. 
  Let {_t}_t=1^∞ be a ^d-valued stochastic process such that _t is ℱ_t-1-measurable. 
  
  Let _i,0 be a d× d positive definite matrix, and define
  _i,t as in (<ref>). 
  Further, for any t ≥ 0, define _t = ∑_s=1^t ϵ_i,s_s. 
  
  Then, for any δ_i > 0, with probability at least (1-δ_i)
      _t_^-1_i,t^2    ≤ 
          2 
          σ_i^2 
          log( 
              1/δ (_i,t)^1/2/ (_i,0)^1/2) 
                 ∀ t≥ 0.
  
  See <cit.>.
  In this work, the filtration{ℱ_t}_t=0^∞corresponds to the  natural filtration associated to the stochastic process _t satisfying (<ref>), which is defined by considering theσ-algebraℱ_t=σ(_1,,_t+1,ϵ_i,0,,ϵ_i,t), where_t=(_t,_t).  
  
  Lemma <ref> corresponds to <cit.> with slightly modified notations, 
  
  substituting(X,η,θ,V̅_t,V)with(,(ϵ_i,1, ,ϵ_i,t),_̨i,_i,t,_i,0)for any arbitrary i=1,,n. 
  
  
  
  We stress that (<ref>) holds jointly for all timest ≥  0such that((<ref>)) ≥  (1 - δ_i). This result is key to ensure that the trajectory-wise chance constraint (<ref>) is satisfied and to guarantee the safety and recursive feasibility of our framework. 
  
  Next, we restate and prove Theorem <ref>. 
  
  
  Consider the true system  (<ref>),
  
      _t+1 = (_t, _t) + (_t,_t,) + _t,
  
  where each ϵ_i,t is σ_i-subgaussian and bounded (ϵ_i,t∈ℰ_i).   
  Consider the meta-learning model  = (ĝ_1,0.9,ĝ_n) 
  defined for each 
  1i = 1,,n as ĝ_i(, ,_̨i) = _̨i^⊤_i(,), where _i: ^n × ^m → ^ and   
  _̨i ∈ ^. 
  Starting from (_i,0,_i,0), with _i,0∈^, and _i,0 is a d× d positive definite matrix, define the sequence {(_i,s,_i,s)}_s=0^t,  
  where (_i,t, _i,t) is computed with online data from (<ref>) using (<ref>):
  
      [              _i,t = Φ_i,t-1^⊤Φ_i,t-1 + _i,0,; _i,t = _i,t^-1 ( Φ^⊤_i,t-1_i,t + _i,0_i,0 ), ]   
          i = 1,,,
  
  where _i,t^⊤=  [x_i,1- h_i(_0,_0), , x_i,t- h_i(_t-1,_t-1)]  ∈ ^t, and
  Φ_i,t-1^⊤= [_i(_0, _0), , _i(_t-1,_t-1)]∈^× t. 
  
  
  
  
  
  
  
  
  Further, define δ_i=δ/(2), and let
  
      β_i,t^δ = σ_i (
          √( 2 log(1/δ_i(_i,t)^1/2/(_i,0)^1/2) )+√(λ̅(_i,0)/λ(_i,t)χ^2_(1-δ_i))),
          
      and    _i,t^δ(_i,t,_i,t)  = {_̨i  | _̨i -_i,t__i,t≤ β_i,t^δ}.
  
  Then, under Assumptions <ref> and <ref>, 
  
      (
          ∈_i,t^δ(_i,t,_i,t)
            ∀ t ≥  0
          ) ≥ (1-2δ_i).
  
  This proof is a straightforward extension of <cit.>, where we use Assumption <ref> to provide a probabilistic error bound for the model mismatch over the prior for , and Lemma <ref> to bound the estimation error due to _twith a union bound to obtain β_i.
  
  Define ^i = (ϵ_i,1, , ϵ_i,t)^⊤. 
  
  
  
  For conciseness, we drop the indices i and t and denote
  (,,,_0,_0,Φ,) 
  = 
  (,_i,t,_i,t,_i,0,_i,0,Φ_i,t-1,^i).  
  
  Under Assumption <ref>, we can write 
  _i,t=Φ +. Then we rewrite the mean estimate  of  at time t as
  
      = ( _0 + Φ^⊤Φ )^-1 (_0 _0 +  Φ^⊤ ( Φ +  ) )
         = ( _0 + Φ^⊤Φ )^-1Φ^⊤ + ( _0 + Φ^⊤Φ )^-1 ( _0 + Φ^⊤Φ ) 
           - ( _0 + Φ^⊤Φ )^-1_0 ( - _0) 
         = Φ^⊤ +  - _0 (  - _0 ),
  
  from which we obtain, for any ∈^, that 
  ^⊤ (  -  ) 
  = 
  ^⊤ (Φ^⊤)  -^⊤ (_0 ( - _0)).  
  
  
  
  
  
  With this result, by the Cauchy-Schwarz inequality,
  
      | ^⊤ (  -  ) | 
         ≤__t( Φ^⊤ϵ_ + _0 ( - _0)_) 
      
         ≤_( Φ^⊤_+√(_0/)-_0__0),
  
  where the second inequality is obtained as
  
      0.9_0 ( - _0)^2_≤/_0_0 ( - _0)^2__0
      
         0.9=
      _0/_0 ( - _0)^2__0 
      
      
      =
      _0/ - _0^2__0
      .
  
  
  
  
  
  
  
  
  
  
  
  
  By Lemma <ref>, for any δ_i, 
  with probability at least (1-δ_i), 
  
      Φ^⊤_^2    ≤
      2σ_i^2
      log( 
          1/δ_i ()^1/2/ (_0)^1/2)
             ∀ t≥ 0
           .
   
  By Assumption <ref>, for δ_i =δ/(2), with probability at least (1-δ_i), 
  
  
      - _0__0^2 ≤σ_i^2 χ^2_(1-δ_i)
           
          .
  
  
  
  
  From (<ref>), by a union bound,[((<ref>)∩(<ref>))
  =1-((<ref>)^C ∪(<ref>)^C), where A^C denotes the negation of A. 
  Then, by Boole's inequality, 
  1-((<ref>)^C ∪(<ref>)^C)
  ≥
  1-((<ref>)^C) - ((<ref>)^C)
  =
  -1 +((<ref>)) + ((<ref>)). 
  Finally, using the lower bounds on the probabilities that (<ref>) and (<ref>) occur, we obtain 
  ((<ref>)∩(<ref>))
  ≥ -1 +(1-δ_i) + (1-δ_i) = 1-2δ_i.
  ] we have that with probability at least (1-2δ_i), for all t≥ 0 and any ∈^, 
  
      | ^⊤ (  -  ) | 
         ≤_σ_i  
      ( 
         √( 2 log(  1/δ_i ()^1/2/ (_0)^1/2) )
                     +
          √(_0/χ^2_(1-δ_i))).
  
  Define β_i,t^δ as in (<ref>). Then letting =  ( - ) in the expression above, we obtain
  
      - ^2_   ≤ ( - )_β_i,t^δ.
  
  Since ( - )_ = _t - _, we divide both sides by - _ and obtain
  
      - _   ≤β_i,t^δ  ∀ t≥ 0,
  
  which holds with probability at least (1-2δ_i). Using the definition of the confidence sets _i,t^δ in  (<ref>) and the result above, this concludes our proof.
  
  
   §.§ Relaxing Assumption <ref> and extending Theorem <ref>
  Our probabilistic guarantees rely on the assumption that the meta-learning model is capable of fitting the true dynamics. 
  Specifically, in Assumption <ref>, we assume that for 
  all  and all i = 1,,n,  there exists () ∈^  such that 
   
  |()^⊤_i(, )  -
  g_i(,,)|≤Δ_i   ∀∈, 
  ∀∈,
  
  where Δ_i=0 for all i=1,…,n. 
  This assumption of the true model belonging to the model class is standard in the literature (see Sections <ref> and <ref>) and reasonable in practice. 
  Indeed, our simulated and hardware experiments verified the safety guarantees of the  under this assumption, which shows that assuming Δ_i=0 is not limiting when using expressive meta-learning models. 
  Nevertheless, our results can be generalized to the case where Δ_i>0. 
  Indeed, the proof of Theorem <ref> follows by adapting the proof of <cit.> to use the prior of our Bayesian model in Assumption <ref> and our particular model structure. 
  As such, below, we extend <cit.> to also account for bounded model mismatch. Due to space constraints, we leave the simple extensions of Theorem <ref> and <ref> (and the evaluation of the bounds Δ_i, e.g.,  depending on the boundedness of the features _i(·) and the chosen meta-learning architecture) for future work.
  [Confidence sets for linear models corrupted by subgaussian and bounded noise]
  
  
  On a probability space (Ω,𝒢,), 
  let (ℱ_t)_t∈ℕ be a filtration. 
  Let 
  (ϵ_t)_t∈ℕ and (γ_t)_t∈ℕ be real-valued stochastic processes adapted to ℱ_t such that
  
    *  [ϵ_t|ℱ_t-1]≤exp(λ^2R^2/2) for all λ∈ for all t∈ℕ.
  
    * |γ_t|≤Δ almost surely for some constant Δ<∞. 
  
  Let (x_t)_t∈ℕ be a ^d-valued stochastic process adapted to ℱ_t-1.  
  
  For fixed parameters θ∈^d, consider the linear model
  
  y_t = 
  θ^⊤ x_t 
  + γ_t + ϵ_t
  
  and assume that θ_2≤ S.  
  
  Given t observations tuples (x_s,y_s)_s=1^t, define
  X_t=[x_1,,x_t]^⊤∈^t× d, 
  Y_t=[y_1,,y_t]^⊤∈^t, 
   and 
  
      θ̅_t = (X_t^⊤ X_t+λ I)^-1X_t^⊤ Y_t,
            V̅_t = X_t^⊤ X_t + λ I.
  
  	
  	Given δ∈(0,1), 
  define the ellipsoidal confidence set
  
      𝒞̅_t^δ
      =
      {θ∈^d :
      θ-θ̅_t_V̅_t≤
          Rβ_t
          +
          λ^1/2S
          +
          D_t
      },
   
  where 
  	β_t^2 
  	=	
  	2 log(1/δ((V̅_t)/(λ I))^1/2) and 
  D_t = σ_max(L_t^⊤)  
   t^1/2Δ where L_t∈^d× d is such that 𝒳_t^⊤V̅_t 𝒳_t=L_t L_t^⊤ with 𝒳_t=(X_t^⊤ X_t+λ I)^-1X_t^⊤. 
  Then,
  
      (
          θ∈𝒞̅_t^δ ∀ t ≥  0
          ) ≥ 1-δ.
   
  Define 
  
  y̅_t = θ^⊤ x_t + ϵ_t so that  y_t = y̅_t + γ_t. 
  
  Denoting Y̅_t=[y̅_1,,y̅_t] 
   and 
   Γ_t=[γ_1,,γ_t], we rewrite θ̅_t as 
  
      θ̅_t    = 
          (X_t^⊤ X_t+λ I)^-1X_t^⊤ Y_t
          
         =
          (X_t^⊤ X_t+λ I)^-1X_t^⊤Y̅_t
          +
          (X_t^⊤ X_t+λ I)^-1X_t^⊤Γ_t
          
         =
          𝒳_tY̅_t
          +
          𝒳_tΓ_t
           ≜ θ̅_ϵ + Δθ̅_γ,
  
  where 
  𝒳_t=(X_t^⊤ X_t+λ I)^-1X_t^⊤. 
  We then obtain that 
  
      θ - θ̅_t_V̅_t = θ-(θ̅_ϵ+Δθ̅_γ)_V̅_t ≤ θ - θ̅_ϵ_V̅_t+Δθ̅_γ_V̅_t.
  
  The second term 
  Δθ̅_γ_V̅_t is bounded with probability one since |γ_t|≤Δ. 
  Indeed, denoting L_t∈^d× d for the Cholesky decomposition 𝒳_t^⊤V̅_t 𝒳_t=L_t L_t^⊤, we have that
  
      Δθ̅_γ_V̅_t^2
         =
      𝒳_t Γ_t _V̅_t^2
      =
      Γ_t^⊤𝒳_t^⊤V̅_t 𝒳_t Γ_t
      =
      Γ_t^⊤ L_tL_t^⊤Γ_t
      
      
      
         =
      (L_t^⊤Γ_t)^⊤ 
       (L_t^⊤Γ_t)
       =
        L_t^⊤Γ_t _2^2
      	
         ≤ 
       ( L_t^⊤_2 ·Γ_t _2)^2
      	
         ≤ 
       (σ_max(L_t^⊤)  
       √(t)Γ_t _∞)^2
      	
         ≤ 
       (σ_max(L_t^⊤)  
       t^1/2Δ)^2
       = 
       D_t^2.
  
  Combining this bound that holds almost surely with the result from
  Theorem <cit.> to bound the first term θ-θ̅_ϵ_V̅_t in (<ref>) at all times with probability at least 1-δ, this concludes our proof.
  
  
   §.§ Robust Obstacle Avoidance Constraints Implementation
  
  We describe our implementation of obstacle avoidance constraints in the special case where all obstacles are spheres in^3and the system is a manipulator. Other types of constraints (e.g., linear constraints), obstacles (e.g., polyhedral and non-convex obstacles), as well as point-mass systems, are particular cases of the following derivations, as we discuss next.
  
  
  
  Consider a function:^n→^3mapping the state of the system (its configuration) to a point that should not intersect with obstacles (e.g., for a manipulator,is the end-effector position as a function of the joint angles. For a point-mass system,maps to the positional variables only.). 
  
  
  
  Before treating the robust case under uncertainty, consider this constraint for a single state∈^n. 
  
  Given a spherical obstacle={∈^3  |   -_o_2< r_o,  _o∈^3,r_o> 0}, we define the obstacle avoidance constraint as()∉, which is equivalent to()-_o_2-r_o≥0. 
  Definingd()=()-_o_2-r_o, we observe thatd()≥0is a non-convex constraint. To tackle this type of constraint, sequential convex programming (SCP) consists of solving a sequence of convex programs where any non-convex constraint is linearized. Specifically, the linearization of this constraint at^j ∈ ^nisd(^j) + ((^j)-_o/(^j)-_o_2)^⊤(∂/∂(^j))(-^j)≥ 0,where the term multiplying(-^j)corresponds to∂d/∂(^j), 
  and∂/∂(^j) ∈ ^3×6denotes the Jacobian ofevaluated at^j. 
  This constraint is linear inand can be handled by off-the-shelf convex optimization algorithms. 
  Note that the case whereis an arbitrary non-convex set can be handled similarly by replacingd()with a signed distance function, see <cit.>. 
  
  
  Next, we treat the robust constraint()∉ ∀∈, with⊂^6a reachable set. 
  A conservative reformulation consists of defining an outer-bounding rectangular set={∈^6 |  |x_i-μ_i|≤δ_i, ∈^6,δ_i>0}such that⊆. 
  
  In this work, as we use a sampling-based algorithm to approximate each reachable set,is finite, i.e., = {^j}_j=1^M. Thus, we simply choose the nominal predicted state_k(see ) for the centerof each rectangular set, and we chooseδ_i=max_j|x_i^j-μ_i|for its half-width. 
  
  Then, we define an  outer-bounding ellipsoidal set= {∈^n |  (-)^⊤^-1(-)≤1}, where the centeris the same as the center of the rectangular setand=n·(δ_i^2, i = 1,,n)is a diagonal positive definite matrix.[Observe that  typically  only depends on positional and angular variables, so the entries of the last rows of ∂ d/∂(^j) are zero. Thus, one can reduce conservatism by formulating a constraint with =n_pos·(δ_i^2, i = 1,,n_pos) instead, where n_pos < n (e.g., for a three-joint manipulator, n_pos = 3 < 6 = n).] This construction guarantees that⊆⊆, so thatouter-bounds. Finally, we use this set to reformulate the robust contraintd()≥0  ∀∈asd() - √(^⊤)≥0, 
  where=∂d/∂(), 
  andis the nominal state of the trajectory, see . 
  For a justification of this reformulation, refer to <cit.>. 
  
  Note that this constraint is non-convex in, and also depends on the control trajectorydue to the dependence on(which depends on the sampled trajectory when performing reachability analysis). 
  
  We iteratively linearize this last constraint within each SCP iteration to obtain a quadratic program which is efficiently solved using OSQP <cit.>. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  [
      < g r a p h i c s >
  ]Thomas Lew is a Ph.D. candidate in Aeronautics and Astronautics at Stanford University. He received his BSc. degree in Microengineering from École Polytechnique Fédérale de Lausanne in 2017 and his MSc. degree in Robotics from ETH Zurich in 2019. His research focuses on the intersection between optimization, control theory, and machine learning techniques for robotics and aerospace applications.
  [
      < g r a p h i c s >
  ]Apoorva Sharma
  is a PhD candidate in the Autonomous Systems Lab at Stanford University. He received an M.S. degree in Aeronautics and Astronautics from Stanford University in 2018, and a B.S. degree in Engineering from Harvey Mudd College in 2016. His research interests center around the application of machine learning models in control and sequential decision making applications, with a focus on Bayesian methods and quantifying uncertainty to enable safe application of machine learning systems in robot autonomy and decision making. 
  [
      < g r a p h i c s >
  ]James Harrison
  is a Ph.D. candidate in the Autonomous Systems Lab at Stanford University. He received an M.S. degree from Stanford University in 2018 and a B.Eng. degree from McGill University in 2015, both in mechanical engineering. His research interests include few-shot, adaptive, and open-world learning, Bayesian deep learning, and applications in safe robot autonomy, decision-making, and control. 
  [
      < g r a p h i c s >
  ]Andrew Bylard is a Ph.D. candidate in Aeronautics and Astronautics at Stanford University. He received a B.Eng. with a dual concentration in Mechanical Engineering and Electrical Engineering from Walla Walla University in 2014 and a M.Sc. in Aeronautics and Astronautics from Stanford University in 2016. Andrew's research interests include real-time trajectory planning and optimization and unconventional space robotics, including using gecko-inspired adhesives for microgravity manipulation and repurposing rollable extendable booms for long-reach mobile manipulators in reduced gravity. He has been a lead developer at the Stanford Space Robotics Facility, where he designed robot test beds used to perform spacecraft contact dynamics experiments and develop autonomous space robotics capabilities under simulated frictionless, microgravity conditions.
  [
      < g r a p h i c s >
  ]Marco Pavone 
  is   an   Associate Professor of Aeronautics and Astronautics at Stanford University,  where he is the Director of the Autonomous Systems Laboratory. Before  joining  Stanford, he  was  a  Research  Technologist  within the  Robotics  Section  at  the  NASA  Jet Propulsion  Laboratory.   He  received  a Ph.D. degree in Aeronautics and Astronautics from the Massachusetts Institute of  Technology  in  2010.   His  main  research  interests  are  in the  development  of  methodologies  for  the  analysis,  design, and  control  of  autonomous  systems,  with  an  emphasis  on self-driving cars, autonomous aerospace vehicles, and future mobility systems.  He is a recipient of a number of awards, including  a  Presidential  Early  Career  Award  for  Scientists and Engineers, an ONR YIP Award, an NSF CAREER Award, and a NASA Early Career Faculty Award.  He was identified by the American Society for Engineering Education (ASEE) as one of America’s 20 most highly promising investigators under the age of 40.  He is currently serving as an Associate Editor for the IEEE Control Systems Magazine.