ENFORCING ROBUST CONTROL GUARANTEES WITHIN NEURAL NETWORK POLICIES
Year: 2021

Abstract:
When designing controllers for safety-critical systems, practitioners often face a challenging tradeoff between robustness and performance. While robust control methods provide rigorous guarantees on system stability under certain worst-case disturbances, they often yield simple controllers that perform poorly in the average (non-worst) case. In contrast, nonlinear control methods trained using deep learning have achieved state-of-the-art performance on many control tasks, but often lack robustness guarantees. In this paper, we propose a technique that combines the strengths of these two approaches: constructing a generic nonlinear control policy class, parameterized by neural networks, that nonetheless enforces the same provable robustness criteria as robust control. Specifically, our approach entails integrating custom convex-optimization-based projection layers into a neural network-based policy. We demonstrate the power of this approach on several domains, improving in average-case performance over existing robust control methods and in worst-case stability over (non-robust) deep RL methods.

heading: 1 INTRODUCTION
text: The field of robust control, dating back many decades, has been able to provide rigorous guarantees on when controllers will succeed or fail in controlling a system of interest. In particular, if the uncertainties in the underlying dynamics can be bounded in specific ways, these techniques can produce controllers that are provably robust even under worst-case conditions. However, as the resulting policies tend to be simple (i.e., often linear), this can limit their performance in typical (rather than worst-case) scenarios. In contrast, recent high-profile advances in deep reinforcement learning have yielded state-of-the-art performance on many control tasks, due to their ability to capture complex, nonlinear policies. However, due to a lack of robustness guarantees, these techniques have still found limited application in safety-critical domains where an incorrect action (either during training or at runtime) can substantially impact the controlled system.
In this paper, we propose a method that combines the guarantees of robust control with the flexibility of deep reinforcement learning (RL). Specifically, we consider the setting of nonlinear, time-varying systems with unknown dynamics, but where (as common in robust control) the uncertainty on these dynamics can be bounded in ways amenable to obtaining provable performance guarantees. Building upon specifications provided by traditional robust control methods in these settings, we construct a new class of nonlinear policies that are parameterized by neural networks, but that are nonetheless provably robust. In particular, we project the outputs of a nominal (deep neural network-based) controller onto a space of stabilizing actions characterized by the robust control specifications. The resulting nonlinear control policies are trainable using standard approaches in deep RL, yet are guaranteed to be stable under the same worst-case conditions as the original robust controller.
We describe our proposed deep nonlinear control policy class and derive efficient, differentiable projections for this class under various models of system uncertainty common in robust control. We demonstrate our approach on several different domains, including synthetic linear differential inclusion (LDI) settings, the cart-pole task, a quadrotor domain, and a microgrid domain. Although these domains are simple by modern RL standards, we show that purely RL-based methods often produce unstable policies in the presence of system disturbances, both during and after training. In contrast, we show that our method remains stable even when worst-case disturbances are present, while improving upon the performance of traditional robust control methods.

heading: 2 RELATED WORK
text: We employ techniques from robust control, (deep) RL, and differentiable optimization to learn provably robust nonlinear controllers. We discuss these areas of work in connection to our approach.
Robust control. Robust control is concerned with the design of feedback controllers for dynamical systems with modeling uncertainties and/or external disturbances (Zhou and Doyle, 1998; Başar and Bernhard, 2008), specifically controllers with guaranteed performance under worst-case conditions. Many classes of robust control problems in both the time and frequency domains can be formulated using linear matrix inequalities (LMIs) (Boyd et al., 1994; Kothare et al., 1996); for reasonably-sized problems, these LMIs can be solved using off-the-shelf numerical solvers based on interior-point or first-order (gradient-based) methods. However, providing stability guarantees often requires the use of simple (linear) controllers, which greatly limits average-case performance. Our work seeks to improve performance via nonlinear controllers that nonetheless retain the same stability guarantees.
Reinforcement learning (RL). In contrast, RL (and specifically, deep RL) is not restricted to simple controllers or problems with uncertainty bounds on the dynamics. Instead, deep RL seeks to learn an optimal control policy, represented by a neural network, by directly interacting with an unknown environment. These methods have shown impressive results in a variety of complex control tasks (e.g., Mnih et al. (2015); Akkaya et al. (2019)); see Buşoniu et al. (2018) for a survey. However, due to its lack of safety guarantees, deep RL has been predominantly applied to simulated environments or highly-controlled real-world problems, where system failures are either not costly or not possible.
Efforts to address the lack of safety and stability in RL fall into several main categories. The first tries to combine control-theoretic ideas, predominantly robust control, with the nonlinear control policy benefits of RL (e.g., Morimoto and Doya (2005); Abu-Khalaf et al. (2006); Feng et al. (2009); Liu et al. (2013); Wu and Luo (2013); Luo et al. (2014); Friedrich and Buss (2017); Pinto et al. (2017); Jin and Lavaei (2018); Chang et al. (2019); Han et al. (2019); Zhang et al. (2020)). For example, RL has been used to address stochastic stability in H∞ control synthesis settings by jointly learning Lyapunov functions and policies in these settings (Han et al., 2019). As another example, RL has been used to address H∞ control for continuous-time systems via min-max differential games, in which the controller and disturbance are the “minimizer” and “maximizer” (Morimoto and Doya, 2005). We view our approach as thematically aligned with this previous work, though our method is able to capture not only H∞ settings, but also a much broader class of robust control settings.
Another category of methods addressing this challenge is safe RL, which aims to learn control policies while maintaining some notion of safety during or after learning. Typically, these methods attempt to restrict the RL algorithm to a safe region of the state space by making strong assumptions about the smoothness of the underlying dynamics, e.g., that the dynamics can be modeled as a Gaussian process (GP) (Turchetta et al., 2016; Akametalu et al., 2014) or are Lipschitz continuous (Berkenkamp et al., 2017; Wachi et al., 2018). This framework is in theory more general than our approach, which requires using stringent uncertainty bounds (e.g. state-control norm bounds) from robust control. However, there are two key benefits to our approach. First, norm bounds or polytopic uncertainty can accommodate sharp discontinuities in the continuous-time dynamics. Second, convex projections (as used in our method) scale polynomially with the state-action size, whereas GPs in particular scale exponentially (and are therefore difficult to extend to high-dimensional problems).
A third category of methods uses Constrained Markov Decision Processes (C-MDPs). These methods seek to maximize a discounted reward while bounding some discounted cost function (Altman, 1999; Achiam et al., 2017; Taleghan and Dietterich, 2018; Yang et al., 2020). While these methods do not require knowledge of the cost functions a-priori, they only guarantee the cost constraints hold during test time. Additionally, using C-MDPs can yield other complications, such as optimal policies being stochastic and the constraints only holding for a subset of states.
Differentiable optimization layers. A great deal of recent work has studied differentiable optimization layers for neural networks: e.g., layers for quadratic programming (Amos and Kolter, 2017), SAT solving (Wang et al., 2019), submodular optimization (Djolonga and Krause, 2017; Tschiatschek et al., 2018), cone programs (Agrawal et al., 2019), and other classes of optimization problems (Gould et al., 2019). These layers can be used to construct neural networks with useful inductive bias for particular domains or to enforce that networks obey hard constraints dictated by the settings in which they are used. We create fast, custom differentiable optimization layers for the latter purpose, namely, to project neural network outputs into a set of certifiably stabilizing actions.

heading: 3 BACKGROUND ON LQR AND ROBUST CONTROL SPECIFICATIONS
text: In this paper, our aim is to control nonlinear (continuous-time) dynamical systems of the form
ẋ(t) ∈ A(t)x(t) +B(t)u(t) +G(t)w(t), (1)
where x(t) ∈ Rs denotes the state at time t; u(t) ∈ Ra is the control input; w(t) ∈ Rd captures both external (possibly stochastic) disturbances and any modeling discrepancies; ẋ(t) denotes the time derivative of the state x at time t; and A(t) ∈ Rs×s, B(t) ∈ Rs×a, G(t) ∈ Rs×d. This class of models is referred to as linear differential inclusions (LDIs); however, we note that despite the name, this class does indeed characterize nonlinear systems, as, e.g., w(t) can depend arbitrarily on x(t) and u(t) (though we omit this dependence in the notation for brevity). Within this class of models, it is often possible to construct robust control specifications certifying system stability. Given such specifications, our proposal is to learn nonlinear (deep neural network-based) policies that provably satisfy these specifications while optimizing some objective of interest. We start by giving background on the robust control specifications and objectives considered in this work.

heading: 3.1 ROBUST CONTROL SPECIFICATIONS
text: In the continuous-time, infinite-horizon settings we consider here, the goal of robust control is often to construct a time-invariant control policy u(t) = π(x(t)), alongside some certification that guarantees that the controlled system will be stable (i.e., that trajectories of the system will converge to an equilibrium state, usually x = 0 by convention; see Haddad and Chellaboina (2011) for a more formal definition). For many classes of systems,1 this certification is typically in the form of a positive definite Lyapunov function V : Rs → R, with V (0) = 0 and V (x) > 0 for all x 6= 0, such that the function is decreasing along trajectories – for instance,
V̇ (x(t)) ≤ −αV (x(t)) (2)
for some design parameter α > 0. (This particular condition implies exponential stability with a rate of convergence α.2) For certain classes of bounded dynamical systems, time-invariant linear control policies u(t) = Kx(t), and quadratic Lyapunov functions V (x) = xTPx, it is possible to construct such guarantees using semidefinite programming. For instance, consider the class of norm-bounded LDIs (NLDIs)
ẋ = Ax(t) +Bu(t) +Gw(t), ‖w(t)‖2 ≤ ‖Cx(t) +Du(t)‖2, (3)
where A ∈ Rs×s, B ∈ Rs×a, G ∈ Rs×d, C ∈ Rk×s, and D ∈ Rk×a are time-invariant and known, and the disturbance w(t) is arbitrary (and unknown) but obeys the norm bounds above.3 For these systems, it is possible to specify a set of stabilizing policies via a set of linear matrix inequalities (LMIs, Boyd et al. (1994)):[
AS + SAT + µGGT +BY + Y TBT + αS SCT + Y TDT
CS +DY −µI
] 0, S 0, µ > 0, (4)
where S ∈ Rs×s and Y ∈ Ra×s. For matrices S and Y satisfying (4), K = Y S−1 and P = S−1 are then a stabilizing linear controller gain and Lyapunov matrix, respectively. While the LMI above is specific to NLDI systems, this general paradigm of constructing stability specifications using LMIs applies to many settings commonly considered in robust control (e.g., settings with normbounded disturbances or polytopic uncertainty, or H∞ control settings). More details about these types of formulations are given in, e.g., Boyd et al. (1994); in addition, we provide the relevant LMI constraints for the settings we consider in this work in Appendix A.
1In this work, we consider sub-classes of system (1) that may indeed be stochastic (e.g., due to a stochastic external disturbance w(t)), but that can be bounded so as to be amenable to deterministic stability analysis. However, other settings may require stochastic stability analysis; please see Astrom (1971).
2See, e.g., Haddad and Chellaboina (2011) for a more rigorous definition of (local and global) exponential stability. Condition (2) comes from Lyapunov’s Theorem, which characterizes various notions of stability using Lyapunov functions.
3A slightly more complex formulation involves an additional term in the norm bound, i.e.,Cx(t)+Du(t)+ Hw(t), which creates a quadratic inequality in w. The mechanics of obtaining robustness specifications in this setting are largely the same as presented here, though with some additional terms in the equations. As such, as is often done, we assume that H = 0 for simplicity.

heading: 3.2 LQR CONTROL OBJECTIVES
text: In addition to designing for stability, it is often desirable to optimize some objective characterizing controller performance. While our method can optimize performance with respect to any arbitrary cost or reward function, to make comparisons with existing methods easier, for this paper we consider the well-known infinite-horizon “linear-quadratic regulator” (LQR) cost, defined as∫ ∞
0
( x(t)TQx(t) + u(t)TRu(t) ) dt, (5)
for some Q ∈ Ss×s 0 and R ∈ Sa×a 0. If the control policy is assumed to be time-invariant and linear as described above (i.e., u(t) = Kx(t)), minimizing the LQR cost subject to stability constraints can be cast as an SDP (see, e.g., Yao et al. (2001)) and solved using off-the-shelf numerical solvers – a fact that we exploit in our work. For example, to obtain an optimal linear time-invariant controller for the NLDI systems described above, we can solve
minimize S,Y
tr(QS) + tr(R1/2Y S−1Y TR1/2) s. t. Equation (4) holds. (6)

heading: 4 ENFORCING ROBUST CONTROL GUARANTEES WITHIN NEURAL NETWORKS
text: We now present the main contribution of our paper: A class of nonlinear control policies, potentially parameterized by deep neural networks, that is guaranteed to obey the same stability conditions enforced by the robustness specifications described above. The key insight of our approach is as follows: While it is difficult to derive specifications that globally characterize the stability of a generic nonlinear controller, if we are given known robustness specifications, we can create a sufficient condition for stability by simply enforcing that our policy satisfies these specifications at all t. For instance, given a known Lyapunov function, we can enforce exponential stability by ensuring that our policy sufficiently decreases this function (e.g., satisfies Equation (2)) at any given x(t).
In the following sections, we present our nonlinear policy class, as well as our general framework for learning provably robust policies using this policy class. We then derive the instantiation of this framework for various settings of interest. In particular, this involves constructing (custom) differentiable projections that can be used to adjust the output of a nominal neural network to satisfy desired robustness criteria. For simplicity of notation, we will often suppress the t-dependence of x, u, and w, but we note that these are continuous-time quantities as before.

heading: 4.1 A PROVABLY ROBUST NONLINEAR POLICY CLASS
text: Given a dynamical system of the form (1) and a quadratic Lyapunov function V (x) = xTPx, let
C(x) := {u ∈ Ra | V̇ (x) ≤ −αV (x) ∀ẋ ∈ A(t)x+B(t)u+G(t)w} (7) denote a set of actions that, for a fixed state x ∈ Rs, are guaranteed to satisfy the exponential stability condition (2) (even under worst-case realizations of the disturbance w). We note that this “safe” set is non-empty if P satisfies the relevant LMI constraints (e.g., system (4) for NLDIs) characterizing robust linear time-invariant controllers, as there is then some K corresponding to P such that Kx ∈ C(x) for all states x. Using this set of actions, we then construct a robust nonlinear policy class that projects the output of some neural network onto this set. More formally, consider an arbitrary nonlinear (neural networkbased) policy class π̂θ : Rs → Ra parameterized by θ, and let P(·) denote the projection operator for some set (·). We then define our robust policy class as πθ : Rs → Ra, where
πθ(x) = PC(x)(π̂θ(x)). (8) We note that this policy class is differentiable if the projections can be implemented in a differentiable manner (e.g., using convex optimization layers (Agrawal et al., 2019), though we construct efficient custom solvers for our purposes). Importantly, as all policies in this class satisfy the stability condition (2) for all states x and at all times t, these policies are certifiably robust under the same conditions as the original (linear) controller for which the Lyapunov function V (x) was constructed.
Given this policy class and some performance objective ` (e.g., LQR cost), our goal is to then find parameters θ such that the corresponding policy optimizes this objective – i.e., to solve
minimize θ ∫ ∞ 0 ` (x, πθ(x) ) dt s. t. ẋ ∈ A(t)x+B(t)πθ(x) +G(t)w. (9)
Algorithm 1 Learning provably robust controllers with deep RL 1: input performance objective ` // e.g., LQR cost 2: input stability requirement // e.g., V̇ (x) ≤ −αV (x) 3: input policy optimizer A // e.g., a planning or RL algorithm 4: compute P , K satisfying LMI constraints // e.g., by optimizing (6) 5: construct specifications C(x) using P // as defined in Equation (7) 6: construct robust policy class πθ using C // as defined in Equation (8) 7: train πθ via A to optimize Equation (9) 8: return πθ
Since πθ is differentiable, we can solve this problem via a variety of approaches, e.g., a modelbased planning algorithm if the true dynamics are known, or virtually any (deep) RL algorithm if the dynamics are unknown.4
This general procedure for constructing stabilizing controllers is summarized in Algorithm 1. While seemingly simple, this formulation presents a powerful paradigm: by simply transforming the output of a neural network, we can employ an expressive policy class to optimize an objective of interest while ensuring the resultant policy will stabilize the system during both training and testing.
We instantiate our framework by constructing “safe” sets C(x) and their associated (differentiable) projections PC(x) for three settings of interest: NLDIs, polytopic linear differential inclusions (PLDIs), andH∞ control settings. As an example, we describe this procedure below for NLDIs, and refer readers to Appendix B for corresponding formulations for the additional settings we consider.

heading: 4.2 EXAMPLE: NLDIS
text: In order to apply our framework to the NLDI setting (3), we first compute a quadratic Lyapunov function V (x) = xTPx by solving the optimization problem (6) for the given system via semidefinite programming. We then use the resultant Lyapunov function to compute the system-specific “safe” set C(x), and then create a fast, custom differentiable solver to project onto this set.

heading: 4.2.1 COMPUTING SETS OF STABILIZING ACTIONS
text: Given P , we compute CNLDI(x) as the set of actions u ∈ Ra that, for each state x ∈ Rs, satisfy the stability condition (2) at that state under even a worst-case realization of the dynamics (i.e., in this case, even under a worst-case disturbance w). The form of the resultant set is given below.
Theorem 1. Consider the NLDI system (3), some stability parameter α > 0, and a Lyapunov function V (x) = xTPx with P satisfying Equation (4). Assuming P exists, define
CNLDI(x) := { u ∈ Ra | ‖Cx+Du‖2 ≤
−xTPB ‖GTPx‖2 u− x T (2PA+ αP )x 2‖GTPx‖2 } for all states x ∈ Rs. For all x, CNLDI(x) is a non-empty set of actions that satisfy the exponential stability condition (2). Further, CNLDI(x) is a convex set in u. Proof. We seek to find a set of actions such that the condition (2) is satisfied along all possible trajectories of (3). A set of actions satisfying this condition at a given x is given by
CNLDI(x) := { u ∈ Ra | sup
w:‖w‖2≤‖Cx+Du‖2 V̇ (x) ≤ −αV (x)
} .
Let S := {w : ‖w‖2 ≤ ‖Cx+Du‖2}. We can then rewrite the left side of the above inequality as
sup w∈S V̇ (x) = sup w∈S ẋTPx+ xTPẋ = 2xTP (Ax+Bu) + sup w∈S 2xTPGw
= 2xTP (Ax+Bu) + 2‖GTPx‖2‖Cx+Du‖2,
by the definition of the NLDI dynamics and the closed-form minimization of a linear term over an L2 ball. Rearranging yields an inequality of the desired form. We note that by definition of the
4While this problem is infinite-horizon and continuous in time, in practice, one would optimize it in discrete time over a large finite time horizon.
specifications (4), there is some K corresponding to P such that the policy u = Kx satisfies the exponential stability condition (2); thus, Kx ∈ CNLDI, and CNLDI is non-empty. Further, as the above inequality represents a second-order cone constraint in u, this set is convex in u.
We further consider the special case where D = 0, i.e., the norm bound on w does not depend on the control action. This form of NLDI arises in many common settings (e.g., where w characterizes linearization error in a nonlinear system but the dynamics depend only linearly on the action), and is one for which we can compute the relevant projection in closed form (as described shortly).
Corollary 1.1. Consider the NLDI system (3) with D = 0, some stability parameter α > 0, and Lyapunov function V (x) = xTPx with P satisfying Equation (4). Assuming P exists, define
CNLDI-0(x) := { u ∈ Ra | 2xTPBu ≤ −xT (2PA+ αP )x− 2‖GTPx‖2‖Cx‖2 } for all states x ∈ Rs. For all x, CNLDI-0(x) is a non-empty set of actions that satisfy the exponential stability condition (2). Further, CNLDI-0(x) is a convex set in u.
Proof. The result follows by setting D = 0 in Theorem 1 and rearranging terms. As the above inequality represents a linear constraint in u, this set is convex in u.

heading: 4.2.2 DERIVING EFFICIENT, DIFFERENTIABLE PROJECTIONS
text: For the general NLDI setting (3), we note that the relevant projection PCNLDI(x) (see Theorem 1) represents a projection onto a second-order cone constraint. As this projection does not necessarily have a closed form, we must implement it using a differentiable optimization solver (e.g., Agrawal et al. (2019)). For computational efficiency purposes, we implement a custom solver that employs an accelerated projected dual gradient method for the forward pass, and employs implicit differentiation through the fixed point equations of this solution method to compute relevant gradients for the backward pass. Derivations and additional details are provided in Appendix C.
In the case where D = 0 (see Corollary 1.1), we note that the projection operation PCNLDI-0(x) does have a closed form, and can in fact be implemented via a single ReLU operation. Specifically, defining ηT := 2xTPB and ζ := −xT (2PA+ αP )x− 2‖GTPx‖2‖Cx‖2, we see that
PCNLDI-0(x) (π̂(x)) = { π̂(x) if ηT π̂(x) ≤ ζ π̂(x)− η
T π̂(x)−ζ ηT η
η otherwise = π̂(x)− ReLU
( ηT π̂(x)− ζ
ηT η
) η. (10)

heading: 5 EXPERIMENTS
text: Having instantiated our general framework, we demonstrate the power of our approach on a variety of simulated control domains.5 In particular, we evaluate performance on the following metrics:
• Average-case performance: How well does the method optimize the performance objective (i.e., LQR cost) under average (non-worst case) dynamics?
• Worst-case stability: Does the method remain stable even when subjected to adversarial (worst-case) dynamics?
In all cases, we show that our method is able to improve performance over traditional robust controllers under average conditions, while still guaranteeing stability under worst-case conditions.

heading: 5.1 DESCRIPTION OF DYNAMICS SETTINGS
text: We evaluate our approach on five NLDI settings: two synthetic NLDI domains, the cart-pole task, a quadrotor domain, and a microgrid domain. (Additional experiments for PLDI and H∞ control settings are described in Appendix I.) For each setting, we choose a time discretization based on the speed at which the system evolves, and run each episode for 200 steps over this discretization. In all cases except the microgrid setting, we use a randomly generated LQR objective where the matrices Q1/2 and R1/2 are drawn i.i.d. from a standard normal distribution.
5Code for all experiments is available at https://github.com/locuslab/ robust-nn-control
101
105
Lo ss N LD I (D = 0
)
Non-robust Methods Robust Methods
101
105 Lo ss N LD I (D 0 )
100
103
Lo ss
Ca rt
po le
100
103
Lo ss
Q ua
dr ot
or
0 250 500 750 1000 Training epochs
10 1
102
Lo ss
M ic
ro gr
id
0 250 500 750 1000 Training epochs
Setting: MBP PPO RARL Robust MBP * Robust PPO *
Original Adversarial Figure 1: Test performance over training epochs for all learning methods employed in our experiments. For each training epoch (10 updates for the MBP model and 18 for PPO), we report average quadratic loss over 50 episodes, and use “X” to indicate cases where the relevant method became unstable. (Lower loss is better.) Our robust methods (denoted by ∗), unlike the non-robust methods and RARL, remain stable under adversarial dynamics throughout training.
Synthetic NLDI settings. We generate NLDIs of the form (3) with s = 5, a = 3, and d = k = 2 by generating matrices A,B,G,C and D i.i.d. from normal distributions, and producing the disturbance w(t) using a randomly-initialized neural network (with its output scaled to satisfy the norm-bound on the disturbance). We investigate settings both where D 6= 0 and where D = 0. In both cases, episodes are run for 2 seconds at a discretization of 0.01 seconds.
Cart-pole. In the cart-pole task, our goal is to balance an inverted pendulum resting on top of a cart by exerting horizontal forces on the cart. For our experiments, we linearize this system as an NLDI with D 6= 0 (see Appendix D), and add a small additional randomized disturbance satisfying the NLDI bounds. Episodes are run for 10 seconds at a discretization of 0.05 seconds.
Planar quadrotor. In this setting, our goal is to stabilize a quadcopter in the two-dimensional plane by controlling the amount of force provided by the quadcopter’s right and left thrusters. We linearize this system as an NLDI withD = 0 (see Appendix E), and add a small disturbance as in the cart-pole setting. Episodes are run for 4 seconds at a discretization of 0.02 seconds.
Microgrid. In this final setting, we aim to stabilize a microgrid by controlling a storage device and a solar inverter. We augment the system given in Lam et al. (2016) with LQR matrices and NLDI bounds (see Appendix F). Episodes are run for 2 seconds at a discretization of 0.01 seconds.

heading: 5.2 EXPERIMENTAL SETUP
text: We demonstrate our approach by constructing a robust policy class (8) for each of these settings, and optimizing this policy class via different approaches. Specifically, we construct a nominal nonlinear control policy class as π̂θ(x) = Kx+ π̃θ(x), where K is obtained via robust LQR optimization (6), and where π̃θ(x) is a feedforward neural network. To construct the projections PC , we employ the value of P obtained when solving for K. For the purposes of demonstration, we then optimize our robust policy class πθ(x) = PC(π̂θ(x)) using two different methods:
• Robust MBP (ours): A model-based planner that assumes the true dynamics are known.
• Robust PPO (ours): An RL approach based on PPO (Schulman et al., 2017) that does not assume known dynamics (beyond the bounds used to construct the robust policy class).
Robust MBP is optimized using gradient descent for 1,000 updates, where each update samples 20 roll-outs. Robust PPO is trained for 50,000 updates, where each update samples 8 roll-outs; we choose the model that performs best on a hold-out set of initial conditions during training. We note that while we use PPO for our demonstration, our approach is agnostic to the particular method of training, and can be deployed with many different (deep) RL paradigms.
We compare our robust neural network-based method against the following baselines:
• Robust LQR: Robust (linear) LQR controller obtained via Equation (6).
• Robust MPC: A robust model-predictive control algorithm (Kothare et al., 1996) based on state-dependent LMIs. (As the relevant LMIs are not always guaranteed to solve, our implementation temporarily reverts to the Robust LQR policy when that occurs.)
• RARL: The robust adversarial reinforcement learning algorithm (Pinto et al., 2017), which trains an RL agent in the presence of an adversary. (We note that unlike the other robust methods considered here, this method is not provably robust.)
• LQR: A standard non-robust (linear) LQR controller.
• MBP and PPO: The non-robust neural network policy class π̂θ(x) optimized via a modelbased planner and the PPO algorithm, respectively.
In order to evaluate performance, we train all methods on the dynamical settings described in Section 5.1, and evaluate them on two different variations of the dynamics:
• Original dynamics: The dynamical settings described above (“average case”).
• Adversarial dynamics: Modified dynamics with an adversarial test-time disturbance w(t) generated to maximize loss (“worst case”). We generate this disturbance separately for each method described above (see Appendix G for more details).
Initialization states are randomly generated for all experiments. For the synthetic NLDI and microgrid settings, these are generated from a standard normal distribution. For both cart-pole and quadrotor, because our NLDI bounds model linearization error, we must generate initial points within a region where this linearization holds. In particular, the linearization bounds only hold for a specified L∞ ball, BNLDI, around the equilibrium. We use a simple heuristic to construct this ball and jointly find a smaller L∞ ball, Binit, such that there exists a level set L of the Robust LQR Lyapunov function with Binit ⊆ L ⊆ BNLDI (details in Appendix H). Since Robust LQR (and by extension our methods) are guaranteed to decrease the relevant Lyapunov function, this guarantees that these methods will never leave BNLDI when initialized starting from any point inside Binit – i.e., that our NLDI bounds will always hold throughout the trajectories produced by these methods.

heading: 5.3 RESULTS
text: Table 1 shows the performance of the above methods. We report the integral of the quadratic loss over the prescribed time horizon on a test set of states, or indicate cases where the relevant method became unstable (i.e., the loss became orders of magnitude larger than for other approaches). (Sample trajectories for these methods are also provided in Appendix H.)
These results illustrate the basic advantage of our approach. In particular, both our Robust MBP and Robust PPO methods show improved “average-case” performance over the other provably robust methods (namely, Robust LQR and Robust MPC). As expected, however, the non-robust
LQR, MBP, and PPO methods often perform better within the original nominal dynamics, as they are optimizing for expected performance but do not need to consider robustness under worst-case scenarios. However, when we apply allowable adversarial perturbations (that still respect our disturbance bounds), the non-robust LQR, MBP, and PPO approaches diverge or perform very poorly. Similarly, the RARL agent performs well under the original dynamics, but diverges under adversarial perturbations in the generic NLDI settings. In contrast, both of our provably robust approaches (as well as Robust LQR) remain stable under even “worst-case” adversarial dynamics. (We note that the baseline Robust MPC method goes unstable in one instance, though this is due to numerical instability issues, rather than issues with theoretical guarantees.)
Figure 1 additionally shows the performance of all neural network-based methods on the test set over training epochs. While the robust and non-robust MBP and PPO approaches both converge quickly to their final performance levels, both non-robust versions become unstable under the adversarial dynamics very early in the process. The RARL method also frequently destabilizes during training. Our Robust MBP and PPO policies, on the other hand, remain stable throughout the entire optimization process, i.e., do not destabilize during either training or testing. Overall, these results show that our method is able to learn policies that are more expressive than traditional robust methods, while guaranteeing these policies will be stable under the same conditions as Robust LQR.

heading: 6 CONCLUSION
text: In this paper, we have presented a class of nonlinear control policies that combines the expressiveness of neural networks with the provable stability guarantees of traditional robust control. This policy class entails projecting the output of a neural network onto a set of stabilizing actions, parameterized via robustness specifications from the robust control literature, and can be optimized using a model-based planning algorithm if the dynamics are known or virtually any RL algorithm if the dynamics are unknown. We instantiate our general framework for dynamical systems characterized by several classes of linear differential inclusions that capture many common robust control settings. In particular, this entails deriving efficient, differentiable projections for each setting, via implicit differentiation techniques. We show over a variety of simulated domains that our method improves upon traditional robust LQR techniques while, unlike non-robust LQR and neural network methods, remaining stable even under worst-case allowable perturbations of the underlying dynamics.
We believe that our approach highlights the possible connections between traditional control methods and (deep) RL methods. Specifically, by enforcing more structure in the classes of deep networks we consider, it is possible to produce networks that provably satisfy many of the constraints that have typically been thought of as outside the realm of RL. We hope that this work paves the way for future approaches that can combine more structured uncertainty or robustness guarantees with RL, in order to improve performance in settings traditionally dominated by classical robust control.

heading: ACKNOWLEDGMENTS
text: This work was supported by the Department of Energy Computational Science Graduate Fellowship (DE-FG02-97ER25308), the Center for Climate and Energy Decision Making through a cooperative agreement between the National Science Foundation and Carnegie Mellon University (SES00949710), the Computational Sustainability Network, and the Bosch Center for AI. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1745016. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
We thank Vaishnavh Nagarajan, Filipe de Avila Belbute Peres, Anit Sahu, Asher Trockman, Eric Wong, and anonymous reviewers for their feedback on this work.

heading: A DETAILS ON ROBUST CONTROL SPECIFICATIONS
text: As described in Section 3.1, for many dynamical systems of the form (1), it is possible to specify a set of linear, time-invariant policies guaranteeing infinite-horizon exponential stability via a set of LMIs. Here, we derive the LMI (4) provided in the main text for the NLDI system (3), and additionally describe relevant LMI systems for systems characterized by polytopic linear differential inclusions (PLDIs) and for H∞ control settings.
A.1 EXPONENTIAL STABILITY IN NLDIS Consider the general NLDI system (3). We seek to design a time-invariant control policy u(t) = Kx(t) and a quadratic Lyapunov function V (x) = xTPx with P 0 for this system that satisfy the exponential stability criterion V̇ (x) ≤ −αV (x), ∀t. We derive an LMI characterizing such a controller and Lyapunov function, closely following and expanding upon the derivation provided in Boyd et al. (1994).
Specifically, consider the NLDI system (3), reproduced below:
ẋ = Ax+Bu+Gw, ‖w‖2 ≤ ‖Cx+Du‖2. (A.1)
The time derivative of this Lyapunov function along the trajectories of the closed-loop system is
V̇ (x) = ẋTPx+ xTPẋ
= (Ax+Bu+Gw)TPx+ xTP (Ax+Bu+Gw)
= ((A+BK)x+Gw)TPx+ xTP ((A+BK)x+Gw)
= [ x w ]T [ (A+BK)TP + P (A+BK) PG GTP 0 ] [ x w ] .
(A.2)
The exponential stability condition V̇ (x) ≤ −αV (x) is thus implied by inequality[ x w ]T M1 [ x w ] := [ x w ]T [ (A+BK)TP + P (A+BK) + αP PG GTP 0 ] [ x w ] ≤ 0. (A.3)
Additionally, the norm bound on w can be equivalently expressed as[ x w ]T M2 [ x w ] := [ x w ]T [ (C +DK)T (C +DK) 0 0 −I ] [ x w ] ≥ 0. (A.4)
Using the S-procedure, it follows that for some λ ≥ 0, the following matrix inequality is a sufficient condition for exponential stability:
M1 + λM2 0. (A.5) Using Schur Complements, this matrix inequality is equivalent to
(A+BK)TP + P (A+BK) + αP + λ(C +DK)T (C +DK) + 1
λ PGGTP 0. (A.6)
Left- and right-multiplying both sides by P−1, and making the change of variables S = P−1, Y = KS, and µ = 1/λ, we obtain
SAT +AS + Y TBT +BY + αS + 1
µ
( SCT + Y TDT ) (CS +DY ) + µGGT 0. (A.7)
Using Schur Complements again on this inequality, we obtain our final system of linear matrix inequalities as[
AS + SAT + µGGT +BY + Y TBT + αS SCT + Y TDT
CS +DY −µI
] 0, S 0, µ > 0, (A.8)
where then K = Y S−1 and P = S−1. Note that the first matrix inequality is homogeneous; we can therefore assume µ = 1 (and therefore, λ = 1), without loss of generality.
A.2 EXPONENTIAL STABILITY IN PLDIS Consider the setting of polytopic linear differential inclusions (PLDIs), where the dynamics are of the form
ẋ(t) = A(t)x(t) +B(t)u(t), (A(t), B(t)) ∈ Conv{(A1, B1), . . . , (AL, BL)}. (A.9)
Here, A(t) ∈ Rs×s andB(t) ∈ Rs×a can vary arbitrarily over time, as long as they lie in the convex hull (denoted Conv) of the set of points above, where Ai ∈ Rs×s, Bi ∈ Rs×a for i = 1, . . . , L. We seek to design a time-invariant control policy u(t) = Kx(t) and quadratic Lyapunov function V (x) = xTPx with P 0 for this system that satisfy the exponential stability criterion V̇ (x) ≤ −αV (x), ∀t. Such a controller and Lyapunov function exist if there exist S ∈ Rs×s 0 and Y ∈ Ra×s such that
AiS +BiY + SA T i + Y TBTi + αS 0, ∀i = 1, . . . , L, (A.10)
where then K = Y S−1 and P = S−1. The derivation of this LMI follows similarly to that for exponential stability in NLDIs, and is well-described in Boyd et al. (1994).
A.3 H∞ CONTROL Consider the following H∞ control setting with linear time-invariant dynamics
ẋ(t) = Ax(t) +Bu(t) +Gw(t), w ∈ L2, (A.11)
where A, B, and G are time-invariant as for the NLDI case, and where we define L2 as the set of time-dependent signals with finite L2 norm.6
In cases such as these with larger or more unstructured disturbances, it may not be possible to guarantee asymptotic convergence to an equilibrium. In these cases, our goal is to construct a robust controller with bounds on the extent to which disturbances affect some performance output (e.g., LQR cost), as characterized by the L2 gain of the disturbance-to-output map. Specifically, we consider the stability requirement that this L2 gain be bounded by some parameter γ > 0 when disturbances are present, and that the system be exponentially stable in the disturbance-free case. This requirement can be characterized via the condition that for all t and some σ ≥ 0,
E(x, ẋ, u) := V̇ (x) + αV (x) + σ ( xTQx+ uTRu− γ2‖w‖22 ) ≤ 0. (A.12)
We note that when E(x(t), ẋ(t), u(t)) ≤ 0 for all t, both of our stability criteria are met. To see this, note that integrating both sides of (A.12) from 0 to∞ and ignoring the non-negative terms on the left hand side after integration yields∫ ∞
0 (x(t)TQx(t)+u(t)TRu(t))dt ≤ γ2 ∫ ∞ 0 ‖w(t)‖22dt+ (1/σ)V (x(0)). (A.13)
This is precisely the desired bound on the L2 gain of the disturbance-to-output map (see Khalil and Grizzle (2002)). We also note that in the disturbance-free case, substituting w = 0 into (A.12) yields
V̇ (x) ≤ −αV (x)− σ ( xTQx+ uTRu ) ≤ −αV (x), (A.14)
where the last inequality follows from the non-negativity of the LQR cost; this is precisely our condition for exponential stability.
We now seek to design a time-invariant control policy u(t) = Kx(t) and quadratic Lyapunov function V (x) = xTPx with P 0 that satisfies the above condition. In particular, we can write
E (x(t), (A+BK)x(t) +Gw(t),Kx(t))= [ x(t) w(t) ]T M1 [ x(t) w(t) ] , (A.15)
where
M1 :=
[ (A+BK)TP + P (A+BK) + αP + σ(Q+KTRK) PG
GTP −γ2σI
] . (A.16)
6The L2 norm of a time-dependent signal w(t) : [0,∞)→ Rd is defined as √∫∞
0 ‖w(t)‖22dt.
Therefore, we seek to find a P ∈ Rs×s 0 and K ∈ Rs×a that satisfy M1 0, for some design parameters α > 0 and σ > 0. Using Schur complements, the matrix inequalityM1 0 is equivalent to
(A+BK)TP + P (A+BK) + αP + σ(Q+KTRK) + PGGTP/(γ2σ) 0. (A.17) As in Appendix A.1, we left- and right-multiply both sides by P−1, and make the change of variables S = P−1, Y = KS, and µ = 1/σ to obtain
SAT +AS+Y TBT +BY +αS+ 1
µ
( (SQ1/2)(Q1/2S) + (Y TR1/2)(R1/2Y ) ) +µGGT /γ2 0.
Using Schur Complements again, we obtain the LMISAT +AS + Y TBT +BY + αS + µGGT /γ2 [SQ1/2 Y TR1/2][Q1/2S R1/2Y ] −µI  0, S 0, µ > 0, (A.18) where then K = Y S−1, P = S−1, and σ = 1/µ.

heading: B DERIVATION OF SETS OF STABILIZING POLICIES AND ASSOCIATED PROJECTIONS
text: We describe the construction of the set of actions C(x), defined in Equation (7), for PLDI systems (A.9) and H∞ control settings (A.11). (The relevant formulations for the NLDI system (3) are described in the main text.)
B.1 EXPONENTIAL STABILITY IN PLDIS For the general PLDI system (A.9), relevant sets of exponentially stabilizing actions CPLDI are given by the following theorem.
Theorem B.1. Consider the PLDI system (A.9), some stability parameter α > 0, and a Lyapunov function V (x) = xTPx with P satisfying (A.10). Assuming P exists, define
CPLDI(x) := u ∈ R a |  2xTPB1 2xTPB2
... 2xTPBL
u ≤ −  xT (αP + 2PA1)x xT (αP + 2PA2)x
... xT (αP + 2PAL)x
 
for all states x ∈ Rs. For all x, CPLDI(x) is a non-empty set of actions that satisfy the exponential stability condition (2). Further, CPLDI(x) is a convex set in u.
Proof. We seek to find a set of actions such that the condition (2) is satisfied along all possible trajectories of (A.9), i.e., for any allowable instantiation of (A(t), B(t)). A set of actions satisfying this condition at a given x is given by
CPLDI(x) := {u ∈ Ra | V̇ (x) ≤ −αV (x) ∀(A(t), B(t)) ∈ Conv{(A1, B1), . . . , (AL, BL)}.
Expanding the left side of the inequality above, we see that for some coefficients γi ∈ R ≥ 0, i = 1, . . . , L satisfying ∑L i=1 γi(t) = 1,
V̇ (x) = ẋTPx+ xTPẋ = 2xTP (A(t)x+B(t)u)
= 2xTP ( L∑ i=1 γi(t)Aix+ γi(t)Biu ) = L∑ i=1 γi ( 2xTP (Aix+Biu) ) by definition of the PLDI dynamics and of the convex hull. Thus, if we can ensure
2xTP (Aix+Biu) ≤ −αV (x) = −αxTPx, ∀i = 1, . . . , L, then we can ensure that exponential stability holds. Rearranging this condition and writing it in matrix form yields an inequality of the desired form. We note that by definition of the specifications (A.10), there is some K corresponding to P such that the policy u = Kx satisfies all of the above inequalities; thus, Kx ∈ CPLDI(x), and CPLDI(x) is non-empty. Further, as the above inequality represents a linear constraint in u, this set is convex in u.
We note that the relevant projection PCPLDI(x) represents a projection onto an intersection of halfspaces, and can thus be implemented via differentiable quadratic programming (Amos and Kolter, 2017).
B.2 H∞ CONTROL For the H∞ control system (A.11), relevant sets of actions satisfying the condition (A.12) are given by the following theorem.
Theorem B.2. Consider the system (A.11), some stability parameter α > 0, and a Lyapunov function V (x) = xTPx with P satisfying Equation (A.18). Assuming P exists, define
CH∞(x) := { u ∈ Ra | uTRu+ (2BTPx)Tu+ xT ( PA+ATP+αP+Q+γ−2PGGTP ) x ≤ 0 } for all states x ∈ Rs. For all x, CH∞(x) is a non-empty set of actions that guarantee condition (A.12), i.e., that the L2 gain of the disturbance-to-output map is bounded by γ and that the system is exponentially stable in the disturbance-free case. Further, CH∞(x) is convex in u.
Proof. We seek to find a set of actions such that the condition E(x, ẋ, u) ≤ 0 is satisfied along all possible trajectories of (A.11), where E is defined as in (A.12). A set of actions satisfying this condition at a given x is given by
CH∞(x) := {u ∈ Ra | sup w∈L2 E(x, ẋ, u) ≤ 0, ẋ = Ax+Bu+Gw}.
To begin, we note that
E(x,Ax+Bu+Gw, u) = xTP (Ax+Bu+Gw) + (Ax+Bu+Gw)TPx+ αxTPx + σ ( xTQx+ uTRu− γ2‖w‖22 ) We then maximize E over w:
w? = arg max w E(x,Ax+Bu+Gw, u) = GTPx/(σγ2). (B.1)
Therefore,
CH∞(x) = {u | E(x,Ax+Bu+Gw?, u, w?) ≤ 0}. (B.2)
Expanding and rearranging terms, this becomes CH∞(x)={u | uT (σR)u+ (2BTPx)Tu+ xT ( PA+ATP+αP+σQ+PGGTP/(σγ2) ) x ≤ 0}.
(B.3)
We note that by definition of the specifications (A.18), there is some K corresponding to P such that the policy u = Kx satisifies the conditions above (see (A.17)); thus, Kx ∈ CH∞ , and CH∞ is non-empty. We note further that CH∞ is an ellipsoid in the control action space, and is thus convex in u.
We rewrite the set CH∞(x) such that the projection PCH∞ (x) can be viewed as a second-order cone projection, in order to leverage our fast custom solver (Appendix C). In particular, defining P̃ = σR, q̃ = BTPx, and r̃ = xT ( PA+ATP+αP+σQ+PGGTP/(σγ2) ) x, we can rewrite the ellipsoid above as
CH∞(x) = {u | u>P̃ u+ 2q̃>u+ r̃ ≤ 0}. (B.4)
We note that as P̃ 0 and r̃ − q̃>P̃−1q̃ < 0, this ellipsoid is non-empty (see, e.g., section B.1 in Boyd and Vandenberghe (2004)). We can then rewrite the ellipsoid as
CH∞(x) = {u | ‖Ãu+ b̃‖2 ≤ 1} (B.5) where Ã = √
P̃ q̃>P̃−1q̃−r̃ and b̃ =
√ P
q>P−1q−rP −1q. The constraint ‖Ãu + b̃‖2 ≤ 1 is then a
second-order cone constraint in u.

heading: C A FAST, DIFFERENTIABLE SOLVER FOR SECOND-ORDER CONE PROJECTION
text: In order to construct the robust policy class described in Section 4 for the general NLDI system (3) and the H∞ setting (A.11), we must project a nominal (neural network-based) policy onto the second-order cone constraints described in Theorem 1 and Appendix B.2, respectively. As this projection operation does not necessarily have a closed form, we implement it via a custom differentiable optimization solver.
More generally, consider a set of the form
C = {x ∈ Rn | ‖Ax+ b‖2 ≤ cTx+ d} (C.1)
for some A ∈ Rm×n, b ∈ Rm, c ∈ Rn, and d ∈ R. Given some input y ∈ Rn, we seek to compute the second-order cone projection PC(y) by solving the problem
minimize x∈Rn
1 2 ‖x− y‖22
subject to ‖Ax+ b‖2 ≤ cTx+ d. (C.2)
Let F denote the `2 norm cone, i.e., F := {(w, t) | ‖w‖2 ≤ t}. Introducing the auxiliary variable z ∈ Rm+1, we can then rewrite the above optimization problem equivalently as
minimize x∈Rn, z∈Rm+1
1 2 ‖x− y‖22 + 1F (z)
subject to z = [ Ax+ b cTx+ d ] =: Gx+ h,
(C.3)
where for brevity we define G = [ A cT ] and h = [ b d ] , and where 1F denotes the indicator function for membership in the set F . We describe our fast solution technique for computing this projection, as well as our method for obtaining gradients through the solution.
C.1 COMPUTING THE PROJECTION We construct a fast solver for problem (C.3) using an accelerated projected dual gradient method. Specifically, define µ = Rm+1 as the dual variable on the equality constraint in Equation (C.3). The Lagrangian for this problem can then be written as
L (x, z, µ) = 1
2 ‖x− y‖22 + 1F (z) + µT (z −Gx− h), (C.4)
and the dual problem is given by maxµ minx,z L (x, z, µ). To form the dual problem, we minimize the Lagrangian with respect to x and z as
inf x,z L (x, z, µ) = inf x
1
2
{ ‖x− y‖22 − µTGx } + inf
z {µT z + 1F (z)} − µTh. (C.5)
We note that the first term on the right side is minimized at x?(µ) = y +GTµ. Thus, we see that
inf x
1 2 {‖x− y‖22 − µTGx} = − 1 2 µTGGTµ− µTGy. (C.6)
For the second term, denote µ = (µ̃, s) and z = (z̃, t). We can then rewrite this term as
inf z {µT z + 1F (z)}} = inf t≥0 inf z̃ {t · s+ µ̃T z̃ | ‖z̃‖2 ≤ t}. (C.7)
For a fixed t ≥ 0, the above objective is minimized at z̃ = −tµ̃/‖µ̃‖2. (The problem is infeasible for t < 0.) Substituting this minimizer into (C.7) and minimizing the result over t ≥ 0 yields
inf z {µT z + 1F (z)} = inf t≥0 t(s− ‖µ̃‖2) = −1F (µ) (C.8)
where the last identity follows from definition of the second-order cone F . Hence the negative dual problem becomes
minimize µ
1 2 µTGGTµ+ µT (Gy + h) + 1F (µ). (C.9)
We now solve this problem via Nesterov’s accelerated projected dual gradient method (Nesterov, 2013). For notational brevity, define f(µ) := 12µ
TGGTµ + µT (Gy + h). Then, starting from arbitrary µ(−1), µ(0) ∈ Rm+1 we perform the iterative updates
ν(k) = µ(k) + β(k)(µ(k) − µ(k−1)) µ(k+1) = PF ( ν(k) − 1
Lf ∇f(ν(k))
) ,
(C.10)
where Lf = λmax(GGT ) is the Lipschitz constant of f , and PF is the projection operator onto F (which has a closed form solution; see Bauschke (1996)). Letting mf = λmin(GGT ) denote the strong convexity constant of f , the momentum parameter is then scheduled as (Nesterov, 2013)
βk =  k − 1 k + 2 if mf = 0√ Lf − √ mf√
Lf + √ mf
if mf > 0. (C.11)
After computing the optimal dual variable µ?, i.e., the fixed point of (C.10), the optimal primal variable can be recovered via the equation x? = y +GTµ? (as can be observed from the first-order conditions of the Lagrangian (C.4)).
C.2 OBTAINING GRADIENTS In order to incorporate the above projection into our neural network, we need to compute the gradients of all problem variables (i.e., G, h, and y) through the solution x?. In particular, we note that x? has a direct dependence on both G and y, and an indirect dependence on all of G, h, and y through µ?.
To compute the relevant gradients through µ?, we apply the implicit function theorem to the fixed point of the update equations (C.10). Specifically, as these updates imply that µ? = ν?, their fixed point can be written as
µ? = PF ( µ? − 1
Lf ∇f(µ?)
) . (C.12)
DefineM := ∂PF (·)∂(·) ∣∣ (·)=µ?− 1Lf ∇f(µ ?) , and note that∇f(µ?) = GGTµ?+Gy+h. The differential of the above fixed-point equation is then given by
dµ? = M × (
dµ? − 1 Lf
( dGGTµ? +GdGTµ? +GGTdµ? + dGy +Gdy + dh )) . (C.13)
Rearranging terms to separate the differentials of problem outputs from problem variables, we see that(
I −M + 1 Lf
MGGT )
dµ? = − 1 Lf
M ( dGGTµ? +GdGTµ? + dGy +Gdy + dh ) , (C.14)
where I is the identity matrix of appropriate size.
As described in e.g. Amos and Kolter (2017), we can then use these equations to form the Jacobian of µ? with respect to any of the problem variables by setting the differential of the relevant problem variable to I and of all other problem variables to 0; solving the resulting equation for dµ? then yields the value of the desired Jacobian. However, as these Jacobians can be large depending on problem size, we rarely want to form them explicitly. Instead, given some backward pass vector ∂` ∂µ? ∈ R 1×(m+1) with respect to the optimal dual variable, we want to directly compute the gradient
of the loss with respect to the problem variables: e.g., for y, we want to directly form the result of the product ∂`∂µ? ∂µ? ∂y ∈ R 1×n. We do this via a similar method as presented in Amos and Kolter (2017), and refer the reader there for a more in-depth explanation of the method described below.
Define J := I − M + 1Lf MGG T to represent the coefficient of dµ? on the left side of Equation (C.14). Given ∂`∂µ? , we then compute the intermediate term
dµ := −J−T ( ∂`
∂µ?
)T . (C.15)
We can then form the relevant gradient terms directly as( ∂`
∂µ? ∂µ? ∂G
)T = 1 Lf M ( dµ(G Tµ?)T + µ?(GTdµ) T + dµy T )
( ∂`
∂µ? ∂µ? ∂h
)T = 1
Lf Mdµ
( ∂`
∂µ? ∂µ? ∂y
)T = 1
Lf GTMdµ.
(C.16)
In these computations, we note that as our solver returns x?, the backward pass vector we are given is actually ∂`∂x? ∈ R 1×n; thus, we compute ∂`∂µ? = ∂` ∂x? ∂x? ∂µ? = ∂` ∂x?G T for use in Equation (C.15).
Accounting additionally for the direct dependence of some of the problem variables on x? (recalling that x? = y +GTu?), the desired gradients are then given by( ∂`
∂G
)T = ( ∂`
∂x? ∂x? ∂G + ∂` ∂x? ∂x? ∂u? ∂u? ∂G
)T = µ? ∂`
∂x? +
1 Lf M ( dµ(G Tµ?)T + µ?(GTdµ) T + dµy T )
( ∂`
∂h
)T = 
* 0
∂` ∂x? ∂x? ∂h + ∂` ∂x? ∂x? ∂u? ∂u? ∂h T = 1 Lf Mdµ
( ∂`
∂y
)T = ( ∂`
∂x? ∂x? ∂y + ∂` ∂x? ∂x? ∂u? ∂u? ∂y
)T = ( ∂`
∂x?
)T + 1
Lf GTMdµ.
(C.17)

heading: D WRITING THE CART-POLE PROBLEM AS AN NLDI
text: In the cart-pole task, our goal is to balance an inverted pendulum resting on top of a cart by exerting horizontal forces on the cart. Specifically, the state of this system is defined as x = [px, ṗx, ϕ, ϕ̇]
T , where px is the cart position and ϕ is the angular displacement of the pendulum from its vertical position; we seek to stabilize the system at x = ~0 by exerting horizontal forces u ∈ R on the cart. For a pendulum of length ` and mass mp, and for a cart of mass mc, the dynamics of the system are (as described in Tedrake (2009)):
ẋ =  ṗx u+mp sinϕ(`ϕ̇ 2−g cosϕ) mc+mp sin2 ϕ ϕ̇
(mc+mp)g sinϕ−u cosϕ−mp`ϕ̇2 cosϕ sinϕ
l(mc+mp sin2 ϕ)
 , (D.1)
where g = 9.81 m/s2 is the acceleration due to gravity. We rewrite this system as an NLDI by defining ẋ = f(x, u) and then linearizing the system about its equilibrium point as
ẋ = Jf (0, 0) [ x u ] + Inw, ‖w‖ ≤ ‖Cx+Du‖, (D.2)
where Jf is the Jacobian of the dynamics, w = f(x, u)− Jf (0, 0) [x u] T is the linearization error, and In is the n × n identity matrix. We bound this linearization error by numerically obtaining the matrices C and D, assuming that x and u are within a neighborhood of the origin. We describe this process in more detail below. As a note, while we employ an NLDI here to characterize the linearization error, it is also possible to characterize this error via polytopic uncertainty (see Appendix J); we choose to use an NLDI here as it yields a much smaller problem description than a PLDI in this case.
D.1 DERIVING Jf (0, 0) For ẋ = f(x, u), we see that
Jf (x, u) = 0 1 0 0 00 0 ∂p̈x/∂ϕ ∂p̈x/∂ϕ̇ ∂p̈x/∂u0 0 0 1 0 0 0 ∂ϕ̈/∂ϕ ∂ϕ̈/∂ϕ̇ ∂ϕ̈/∂u,  , (D.3) where
∂p̈x ∂ϕ = mp cosϕ
( ϕ̇2l − g cosϕ ) + gmp sin 2 ϕ
mc +mp sin 2 ϕ
− 2mp sinϕ cosϕ
( mp sinϕ ( ϕ̇2l − g cosϕ ) + u )(
mc +mp sin 2 ϕ )2 ,
∂p̈x ∂ϕ̇ = 2ϕ̇lmp sinϕ mc +mp sin 2 ϕ ,
∂p̈x ∂u = 1 mc +mp sin 2 ϕ ,
∂ϕ̈ ∂ϕ =
g(mc +mp) cosϕ+ ϕ̇ 2lmp sin 2 ϕ −ϕ̇2lmp cos2 ϕ+ u sinϕ l ( mc +mp sin 2 ϕ ) − 2mp sinϕ cosϕ(g(mc +mp) sinϕ −ϕ̇2lmp sinϕ cosϕ− u cosϕ l ( mc +mp sin 2 ϕ )2 ,
∂ϕ̈ ∂ϕ̇ = −2ϕ̇mp sinϕ cosϕ mc +mp sin 2 ϕ ,
∂ϕ̈ ∂u = − cosϕ l(mc +mp sin 2 ϕ) .
We thus see that
Jf (0, 0) = 0 1 0 0 00 0 −mpg/mc 0 1/mc0 0 0 1 0 0 0 g(mc+mp)/lmc 0 −1/mc  . (D.4) D.2 OBTAINING C AND D We then seek to construct matrices C and D that bound the linearization error w between the true
dynamics ẋ and our first-order linear approximation Jf (0, 0) [ x u ] . To do so, we bound the error of
this approximation entry-wise: that is, for each entry i = 1, . . . , s, we want to find Fi such that for all x in some region x ≤ x ≤ x̄, and all u in some region u ≤ u ≤ ū,
w2i = ( ∇fi(0) [ x u ] − ẋi )2 ≤ [ x u ]T Fi [ x u ] . (D.5)
Then, given the matrix M = [ F T/2 1 F T/2 2 F T/2 3 F T/2 4 F T/2 5 F T/2 6 ]T (D.6)
we can then obtain C = M1:s and D = Ms:s+m (where the subscripts indicate column-wise indexing).
We solve separately for each Fi to minimize the difference between the right and left sides of Equation (D.5) (while enforcing that the right side is larger than the left side) over a discrete grid of points within x ≤ x ≤ x̄ and u ≤ u ≤ ū. By assuming that Fi is symmetric, we are able to cast this as a linear program in the upper triangular entries of Fi.
To obtain the matrices C and D used for the cart-pole experiments in the main paper, we let x̄ = [1.5 2 0.2 1.5]
T , ū = 10, x = −x̄, and u = −ū. As each entry-wise difference in Equation (D.5) contained exactly three variables (i.e., a total of three entries from x and u), we solved each entry-wise linear program over a mesh grid of 50 points per variable.

heading: E WRITING QUADROTOR AS AN NLDI
text: In the planar quadrotor setting, our goal is to stabilize a quadcopter in the two-dimensional plane by controlling the amount of force provided by the quadcopter’s right and left thrusters. Specifically, the state of this system is defined as x = [px pz ϕ ṗx ṗz ϕ̇]
T , where (px, pz) is the position of the quadcopter in the vertical plane and ϕ is its roll (i.e., angle from the horizontal position); we seek to stabilize the system at x = ~0 by controlling the amount of force u = [ur, ul]
T from right and left thrusters. We assume that our action u is additional to a baseline force of [mg/2 mg/2]T provided by the thrusters by default to prevent the quadcopter from falling. For a quadrotor with mass m, moment-arm ` for the thrusters, and moment of inertia J about the roll axis, the dynamics of this system are then given by (as modified from Singh et al. (2020)):
ẋ =  ṗx cosϕ− ṗz sinϕ ṗx sinϕ+ ṗz cosϕ ϕ̇ ṗzϕ̇− g sinϕ
−ṗxϕ̇− g cosϕ+ g 0
+  0 0 0 0 0 0 0 0
1/m 1/m `/J −`/J u, (E.1) where g = 9.81 m/s2. We linearize this system via a similar method as for the cart-pole setting, i.e., as in Equation (D.2). We describe this process in more detail below. We note that since the dependence of the dynamics on u is linear, we have that D = 0 for our resultant NLDI. As for cart-pole, while we employ an NLDI here to characterize the linearization error, it is also possible to characterize this error via polytopic uncertainty (see Appendix J); we choose to use an NLDI here as it yields a much smaller problem description than a PLDI in this case.
E.1 DERIVING Jf (0, 0) For ẋ = f(x, u), we see that
Jf (x, u) =  0 0 −ṗx sinϕ− ṗz cosϕ cosϕ − sinϕ 0 0 0 0 ṗx cosϕ− ṗz sinϕ sinϕ cosϕ 0 0 0 0 0 0 0 1 0 0 0 −g cosϕ 0 ϕ̇ ṗz 0 0 0 g sinϕ −ϕ̇ 0 −ṗx 0 0 0 0 0 0 0 0  , (E.2) and thus
Jf (0, 0) =  0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 −g 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  . (E.3)
E.2 OBTAINING C AND D We obtain the matricesC andD via a similar method as described in Appendix D, though in practice we only consider the linearization error with respect to x (i.e., since the dynamics are linear with respect to u, we have D = 0). We let x̄ = [1 1 0.15 0.6 0.6 1.3] and x = −x̄. As for cart-pole, each entry wise difference in the equivalent of Equation (D.5) contained exactly three variables (i.e., a total of three entries from x and u), and each entry-wise linear program was solved over a mesh grid of 50 points per variable.

heading: F DETAILS ON THE MICROGRID SETTING
text: For our experiments, we build upon the microgrid setting given in Lam et al. (2016). In this system, the state x ∈ R3 captures voltage deviations, frequency deviations, and the amount of power generated by a diesel generator connected to the grid; the action u ∈ R2 describes the current associated with a storage device and a solar PV inverter; and the disturbance w ∈ R describes the difference between the amount of power demanded and the amount of power produced by solar panels on the grid. The authors also define a performance index y ∈ R2 which captures voltage and frequency deviations (i.e., two of the entries of the state x).
To construct an NLDI of the form (3) for this system, we directly use the A, B, and G matrices given in Lam et al. (2016). We generate C i.i.d. from a normal distribution and let D = 0, to represent the fact that the disturbance w and the entries of the state x are correlated, but that w is likely not correlated with the actions u. Finally, we let Q and R be diagonal matrices with 1 in the entries corresponding to quantities represented in the performance index y, and with 0.1 in the rest of the diagonal entries, to emphasize that the variables in y are the most important in describing the performance of the system.

heading: G GENERATING AN ADVERSARIAL DISTURBANCE
text: In the NLDI settings explored in our experiments, we seek to construct an “adversarial” disturbance w(t) that obeys the relevant norm bounds ‖w(t)‖2 ≤ ‖Cx(t)+Du(t)‖2 while maximizing the loss. To do this, we use a model predictive control method where the actions taken are w(t). Specifically, for each policy π, we model w(t) as a neural network specific to that policy. Every 10 steps of a roll-out, we optimize w(t) through gradient descent to maximize the loss over a horizon of 40 steps, subject to the constraint ‖w(t)‖2 ≤ ‖Cx(t) +Du(t)‖2.

heading: H ADDITIONAL EXPERIMENTAL DETAILS
text: Initial states. To pick initial states in our experiments, for the synthetic settings, we sample each attribute of the state i.i.d. from a standard Gaussian distribution. For cart-pole and planar quadrotor, we sample uniformly from bounds chosen such that the non-robust LQR algorithm (under the original dynamics) did not go unstable. For cart-pole, these bounds were chosen to be px ∈ [−1, 1], ϕ ∈ [−0.1, 0.1], ṗx = ϕ̇ = 0. For planar quadrotor, these bounds were px, pz ∈ [−1, 1], ϕ ∈ [−0.05, 0.05], ṗx = ṗz = ϕ̇ = 0.
Constructing NLDI bounds. Given these initial states, for the cart-pole and quadrotor settings, we needed to construct our NLDI disturbance bounds such that they would hold over the entire trajectory of the robust policy; if not, the robustness specification (A.8) would not hold, and our agent might in fact increase the Lyapunov function. To ensure this approximately, we used a simple heuristic: we ran the (non-robust) LQR agent for a full episode with 50 different starting conditions, and constructed an L∞ ball around all states reached in any of these trajectories. We then used these L∞ balls on the states to construct the matrices C and D for our disturbance bounds, using the procedure described in Appendices D and E.
Computing infrastructure and runtime. All experiments were run on an XPS 13 laptop with an Intel i7 processor. The planar quadrotor and synthetic NLDI experiment with D = 0 took about 1 day to run (since the projections were simple half-space projections), while all the other synthetic domains and cart-pole took about 3 days to run. The majority of the run-time was in computing the adversarial disturbances for test-time evaluations.
Hyperparameter selection. For our experiments, we did not perform large parameter searches. The learning rate we chose for our model-based planner, (both robust and non-robust) remained constant for the different domains; we tried learning rates of 1× 10−3, 1× 10−4, 1× 10−5 and found 1× 10−3 worked best for the non-robust version and 1× 10−4 worked best for the robust version. For our PPO hyperparameters, we simply used those used in the original PPO paper.
One parameter we had to tune for each environment was the time step. In particular, we had to pick a time step high enough that we could run episodes for a reasonable total length of time (within which the non-robust agents would go unstable), but low enough to reasonably approximate a continuoustime setting (since, for our robustness guarantees, we assume the agent’s actions evolve in continuous time). Our search space was small, however, consisting of 0.05, 0.02, 0.01, and 0.005 seconds.
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(a) LQR
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(b) Robust LQR
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(c) MBP
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(d) Robust MBP
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(e) PPO
0 2 4 6 8 10 t
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0 x
x
(f) Robust PPO
Figure H.1: Trajectories of 6 different methods on the cart-pole domain under adversarial dynamics.
Trajectory plots. Figure H.1 shows sample trajectories of different methods in the cart-pole domain under adversarial dynamics. The non-robust LQR and model-based planning approaches both diverge and the non-robust PPO doesn’t diverge, but doesn’t clearly converge after 10 seconds. The robust methods, on the other hand, all clearly converge after 10 seconds.
Runtime comparison. Tables H.1 and H.2 show the evaluation and training time of our methods and the baselines over 50 episodes run in parallel. In the NLDI cases where D = 0, i.e., Generic NLDI (D = 0) and Quadrotor, our projection adds only a very small computational cost. In the other cases, the additional computational cost is more significant, but our method is still far less expensive than the Robust MPC method.
Environment LQR MBP PPO RobustLQR Robust MPC RARL Robust MBP∗ Robust PPO∗
Generic NLDI (D = 0) 0.63 0.61 0.84 0.57 718.06 0.71 0.73 0.94 Generic NLDI (D 6= 0) 0.64 0.62 0.83 0.58 824.86 0.81 15.13 25.38
Cart-pole 0.55 0.67 0.84 0.53 646.90 0.84 10.12 13.37 Quadrotor 0.95 0.98 1.19 0.88 3348.68 1.14 1.15 1.30 Microgrid 0.58 0.61 0.79 0.57 601.90 0.74 8.14 10.25 Generic PLDI 0.57 0.54 0.76 0.51 819.24 0.73 69.35 64.03 Generic H∞ 0.84 0.80 1.03 0.76 N/A 1.00 47.81 63.67
Table H.1: Time (in seconds) taken to run each method on the test set of every environment for 50 episodes run in parallel.
Environment MBP PPO RARL RobustMBP∗ Robust PPO∗
Generic NLDI (D = 0) 26.36 101.77 102.37 30.78 114.60 Generic NLDI (D 6= 0) 26.46 100.79 82.53 221.35 1158.28
Cart-pole 25.49 87.04 98.90 146.34 689.93 Quadrotor 41.24 131.48 112.95 46.13 159.06 Microgrid 23.03 112.52 87.71 113.61 436.64
Table H.2: Time (in minutes) taken to train each method in every environment.

heading: I EXPERIMENTS FOR PLDIS AND H∞ CONTROL SETTINGS
text: In addition to the NLDI settings explored in the main text, we test the performance of our method on PLDI and H∞ control settings. As for the experiments in the main text, we choose a time discretization based on the speed at which the system evolves, and run each episode for 200 steps over this discretization. In both cases, we use a randomly generated LQR objective where the matrices Q1/2 and R1/2 are drawn i.i.d. from a standard normal distribution.
Synthetic PLDI setting. We generate PLDI instances (A.9) with s = 5, a = 3, and L = 3. Specifically, we generate convex hull matrices (A1, B1), . . . , (A3, B3) i.i.d. from normal distributions, and generate (A(t), B(t)) by using a randomly-initialized neural network with softmax output to weight the convex hull matrices. Episodes were run for 2 seconds at a discretization of 0.01 seconds.
Synthetic H∞ setting. We generate H∞ control instances (A.11) with s = 5, a = 3, and d = 2 by generating matricesA,B andG i.i.d. from normal distributions. The disturbancew(t) was produced using a randomly-initialized neural network, with its output scaled to satisfy the L2 bound on the disturbance. Specifically, we scaled the output of the neural network to satisfy an attenuating normbound on the disturbance; at time t, the norm-bound was given by 20× f(2× t/T ), where T is the time horizon and f is the standard normal PDF function. Episodes were run for T = 2 seconds at a discretization of 0.01 seconds.
Results are given in Figure I.1 and Table I.1.

heading: J NOTES ON LINEARIZATION VIA PLDIS AND NLDIS
text: While we linearize the cart-pole and quadrotor dynamics via NLDIs in our experiments, we note that these dynamics can also be characterized via PLDIs. More generally, in this section, we show how we can use the framework of PLDIs to model linearization errors arising in the analysis of nonlinear systems.
Consider the nonlinear dynamical system
ẋ = f (x, u) with f(0, 0) = 0. (J.1)
for x ∈ Rs and u ∈ Ra. Define ξ = (x, u). We would like to represent the above system as a PLDI in the region R := {ξ | ξ ≤ ξ ≤ ξ̄} including the origin. The mean value theorem states that for each component of f , we can write
fi(ξ) = fi(0) +∇fi(z)T ξ, (J.2)
for some z = tξ, where t ∈ [0, 1]. Now, let p = s+ a. Defining the Jacobian of f as
Jf (z) = ∇f1(z) T
... ∇fp(z)T  , (J.3) and recalling that f(0) = 0, we can rewrite (J.2) as
f(ξ) = Jf (z)ξ. (J.4)
Now, suppose we can find component-wise bounds on the matrix Jf (z) overR, i.e,
M ≤ Jf (z) ≤ M̄ for all z ∈ R. (J.5)
We can then write Jf (z) = ∑
1≤i,j≤p
mij(t)Eij with mij(t) ∈ [mij , m̄ij ], (J.6)
where Eij = eieTj and ei is the i-th unit vector in Rp.
We now seek to bound the Jacobian using polytopic bounds. To do this, note that we can write
Jf (z) =
2p 2∑
κ=1 γκAκ γκ ≥ 0, ∑ κ γκ = 1, (J.7)
where Aκ’s are the vertices of the polytope in (J.6), i.e.,
Aκ ∈ V =  ∑ 1≤i,j≤p mijEij | mij ∈ {mij , m̄ij}  . (J.8) Together, Equations (J.2), (J.4), (J.7), and (J.8) characterize the original nonlinear dynamics as a PLDI.
We note that this PLDI description is potentially very large; in particular, the size of V is exponential in the square of the number of non-constant entries in the Jacobian Jf (z), which could be as large as 2p 2 = 2(s+a) 2 . This problem size may therefore become intractable for larger control settings.
We note, however, that we can in fact express this PLDI more concisely as an NLDI. More precisely, we would like to find matrices A,B,C parameterizing the form of NLDI below, which is equivalent to that presented in Equation (3) (see Chapter 4 of Boyd et al. (1994)):
Df(z) ∈ {A+B∆C | ‖∆‖2 ≤ 1} for all z ∈ R. (J.9)
It can shown that the solution to the SDP
minimize tr(V +W )
subject to W 0[ V (Aκ −A)T
Aκ −A W
] 0, ∀Aκ ∈ V
(J.10)
yields the matrices A, B, and C with V = CTC and W = BBT , which can be used to construct NLDI (J.9). While the NLDI here is more concise than the PLDI, the trade-off is that the NLDI norm bounds obtained via this method may be rather loose. As such, for our settings, we obtain NLDI bounds numerically (see Appendices D and E), as these are tighter than NLDI specifications obtained via the above method (though they are potentially slightly inexact). An alternative approach would be to examine how to tighten the conversion from PLDIs to NLDIs, which has been explored in other work (e.g. Kuiava et al. (2013)).

Review summaries:
In this paper, a neural control method is proposed with stability guarantees. The control is assumed to be from a neural network that takes in the state. Stability is guaranteed by projecting the control to the set that satisfies the Lyapunov stability condition for the LQR problem. In particular, minimizing the cost of LQR cost subject to stability constraints can be cast as an SDP for norm-bouned linear differential inclusions. Through making use of the convex optimization layers proposed in Agrawal et al. (2019), the SDP can be added as a layer after the neural policy and efficient projections can be derived such that implicit function theorem can be utilized to differentiate through the fixed point (the optimal conditions of the SDP), such that end to end learning is possible. The proposed approach is compared with the unconstrained method on various tasks.  Both model-based and model-free RL algorithms are used as the neural policy for comparison. The stability-guaranteed approach is able to remain stable even under bounded adversarial dynamics. In comparison, the non-robust methods fail to maintain stability.