A step toward safe and reliable autopilots for flying

Solving Stabilize-Avoid Optimal Control via Epigraph Form and Deep Reinforcement Learning Oswin So, Chuchu Fan Massachusetts Institute of Technology {oswinso, chuchu}@mit.edu Abstract—Tasks for autonomous robotic systems commonly require stabilization to a desired region while maintaining safety specifications. However, solving this multi-objective problem is challenging when the dynamics are nonlinear and high- dimensional, as traditional methods do not scale well and are often limited to specific problem structures. To address this issue, we propose a novel approach to solve the stabilize-avoid problem via the solution of an infinite-horizon constrained optimal control problem (OCP). We transform the constrained OCP into epigraph form and obtain a two-stage optimization problem that optimizes over the policy in the inner problem and over an auxiliary variable in the outer problem. We then propose a new method for this formulation that combines an on-policy deep reinforcement learning algorithm with neural network regression. Our method yields better stability during training, avoids instabilities caused by saddle-point finding, and is not restricted to specific requirements on the problem structure compared to more traditional methods. We validate our approach on different benchmark tasks, ranging from low-dimensional toy examples to an F16 fighter jet with a 17-dimensional state space. Simulation results show that our approach consistently yields controllers that match or exceed the safety of existing methods while providing ten-fold increases in stability performance from larger regions of attraction. Project page can be found at https://mit-realm.github.io/efppo. I. INTRODUCTION Autonomous systems are becoming increasingly prevalent in our lives in recent years; however, designing controllers for complex dynamics systems is difficult. For example, an unmanned aerial vehicle may be required to to observe a target location while maintaining line of sight to a base station [53]. Another example is the problem of satellite docking that requires approaching a target satellite from a specific direction [23]. Robot control tasks often involve both stability and safety requirements, where the controller must both drive the system towards and remain stable within some goal region while avoiding unsafe regions. We denote this as the stabilize-avoid problem. However, synthesizing a policy that achieves both tasks in the presence of input constraints is challenging as these objectives can often be contradictory [22]. Reach-Avoid. Reachability analysis and the reach-avoid problem [50] are very closely related to the stabilize-avoid problem that we tackle in this paper. Given a dynamical system, the reach-avoid problem aims to solve for the set of initial conditions and the appropriate control policy to drive a system to a desired goal set while avoiding undesirable states. Hamilton-Jacobi analysis [27] provides a methodology for computing the solution to reach-avoid problems, and is Robotics: Science and Systems 2023 Fig. 1: Visualization of F16 ground collision avoidance within a low-altitude flight corridor using EFPPO. The system is required to stabilize to a target altitude near the ground (in green), avoid collision with the ground and stay within the flight corridor defined by the two walls and the ceiling. conventionally solved via numerical partial differential equation (PDE) techniques that use state space discretization. These methods are limited in practice to systems with up to 5 continuous state variables [30]. Recent works have applied Deep Reinforcement Learning (DeepRL) to solve reach-avoid problems for higher dimensional systems (e.g., 6-dimensional system in [21]). While reaching a goal set is related to stabilization to a goal set, the two objectives have important differences. The goal of the stabilizing controller is to induce stability of the system within a subset of the goal set. This need not be true for the reach controller, where the goal set may not contain any equilibrium point at all. In the worst case, unwanted oscillations could be introduced into the system. We discuss the relationship between the two formulations in Section II-A. Constrained Reinforcement Learning. Works that address the problem of task completion with safety constraints from the reinforcement learning community usually do so from the constrained Markov Decision Process (CMDP) [5] framework. Many of these works adopt techniques from constrained optimization to handle the additional safety constraints, of which the use of Lagrangian duality is popular due to its simplicity [1, 25, 26, 48]. However, the CMDP formulation considers the discounted sum over constraints instead of enforcing the constraint at each state. This allows the additional arXiv:2305.14154v1 [cs.RO] 23 May 2023 constraint terms to be treated in the same way as the objective function at the expense of allowing constraint violations. While we can modify the CMDP formulation to disallow all constraint violations, this can lead to an ill-conditioned problems. The relationship between our work and the CMDP setting is discussed in Section IV-D. Lyapunov Methods. Lyapunov theory provides an attractive option for synthesizing safe, stabilizing controllers by using control Lyapunov functions (CLFs) [7, 43] and control barrier functions (CBFs) [6, 52]. CLFs and CBFs provide conditions for synthesizing controllers that are certified to be stable and safe respectively. However, they are difficult to construct analytically for general nonlinear systems [16]. CLFs / CBFs can be synthesized via convex optimization (e.g., sum-of- squares programming [2, 12, 47]). However, such approaches are limited to systems with polynomial dynamics and rely on the use of solvers for semidefinite programs which can face numerical difficulties [32]. Alternatively, neural networks can be used to synthesize these certificate functions [11, 13, 34]. However, one problem that remains is that CLFs and CBFs cannot be easily combined to yield combined safety and stability guarantees when the set of feasible controls induced by the CLF and CBF do not intersect, forcing the controller to pick one and sacrifice either safety or stability. This can lead to the presence of unwanted local minima [36]. Although this can be resolved by learning a joint Control Lyapunov Barrier Function (CLBF) [14, 37], the training process requires the ability to sample from the control-invariant set, which is difficult in the case of complex nonlinear dynamics where the control-invariant set is not known. Our experiments show that this method is difficult to apply in practice. We provide discussions on this in Section VII. Model Predictive Control. Online optimization-based control methods such as model predictive control (MPC) have become increasingly popular for general-purpose control synthesis with the increase in computational power available for robotic systems. Moreover, they can be viewed as a finite-horizon approximation to an infinite horizon optimal control problem (OCP), which is closely linked to Lyapunov stability [18]. However, the constrained OCPs that need to be solved online are computationally expensive, making it difficult to achieve high frequency control updates in practice [35]. Moreover, accurate gradient information is typically necessary for solving nonlinear OCPs quickly, making it further difficult to use with dynamics that have expensive gradients or are non- differentiable. Finally, guaranteeing the recursive feasibility of MPC for general nonlinear systems is challenging [28], and in many cases, requires the solution (or approximation) of a control invariant set which can be difficult to find. Continuous-Time Constrained OCP. In addition to the discrete-time formulation used MPC, there exists a number of methods that investigate the problem of constrained OCP in continuous time. Works within this area mainly focus on investigating theoretical properties of the value function for the finite and infinite horizon problems [4, 10, 24, 29, 41, 42]. We note similarities of our problem formulation to the one discussed in [4], where the continuous-time constrained optimal control problem is solved via transformation into epigraph form. An associated Hamilton-Jacobi PDE and its properties are investigated and used as the basis of a numerical PDE solver for a two-dimensional finite-horizon problem. To the best of our knowledge, the transformation of the infinite-horizon constrained OCP problem into epigraph form for discrete-time problems has not been proposed before in literature. In this work, we solve the stabilize-avoid problem by formulating an infinite-horizon constrained OCP, inspired by the global asymptotic stability guarantees in the unconstrained case [33]. Our method for solving the constrained OCP departs from traditional Lagrangian duality based methods and uses an epigraph form which we denote as the Epigraph Form Constrained Optimal Control Problem (EF-COCP). We then solve the EF-COCP using DeepRL by deriving a corresponding policy gradient theorem and applying the proximal policy optimization (PPO) algorithm [39]. This allows us to tackle a wider range of systems compared to non-RL-based methods and handle general nonlinear non-differentiable black box dynamics with minimal computational cost online. Contributions. We summarize our contributions below. • We propose a new formulation of the safety-constrained OCP via an epigraphic reformulation (EF-COCP) which is easier to interpret and avoids the optimization instability of existing Lagrangian duality methods. • We derive a policy-gradient theorem for the inner problem of EF-COCP and propose the EFPPO algorithm for solving the stabilize-avoid problem using DeepRL. • The proposed EFPPO method is validated on a range of challenging systems, yielding promising empirical results on complex systems such as a 17-dimensional F16 fighter jet, visualized in Figure 1. II. THE STABILIZE-AVOID PROBLEM We consider arbitrary nonlinear discrete-time dynamical systems of the form xk+1 = f (xk, uk) (1) where x ∈ X ⊂ Rnx , u ∈ U ⊂ Rnu and f : X × U → X . In this paper, we consider the following control synthesis problem. Problem 1 (Stabilize-Avoid Problem). Given a nonlinear system with a goal set G ⊂ X and avoid set A ⊂ X , find a control policy u = π(x) that maximizes the size of the set R ⊂ X defined as the set of initial states x0 such that all trajectories started from x0 evolving under the dynamics xk+1 = f (xk, π(xk)), x0 ∈ R, (2) also satisfy the following two properties. Stabilize: lim sup k→∞ min y∈G ∥xk − y∥ = 0, Avoid: xk̸ ∈ A for all k ≥ 0. 0 1 2 p −1 0 1 v Stabilize 0 1 2 p Reach 0 1 V (Relative)Fig. 2: Trajectory and value functions comparing a solution to the stabilize and reach objectives of the double integrator in 1D. The reach controller can reach the goal region at p = 1 (red line) faster than stabilize (compare magenta triangles), but the controller never stabilizes to G and induces a periodic orbit. The stable controller reaches G slower but remains in the set. In short, the objective is to reach and (asymptotically) stabilize to a goal set G while avoiding the set of unsafe states A. A. Relationship with the Reach-Avoid Problem Note that we use lim sup when defining stabilize instead of minimizing over time as in the reach formulation min k min y∈G ∥xk − y∥ = 0. (3) A system that enters but then subsequently exits the goal set G will satisfy the reach-avoid problem but not the stabilize-avoid problem, as we illustrate in the following example. Example 1 (Stabilize vs Reach). Consider a double-integrator with states x = [p, v] ∈ R2, control u ∈ R1, and the following task specification: Constraints Goal |u| ≤ 1, A := ∅ G := { x | p = 1.0 } One solution to the reach problem (minimizing the first hitting time) yields a controller with periodic orbits, while minimizing the stabilize objective yields a globally stabilizing controller on G (Figure 2). We note that both formulations yield similar results when G consists only of equilibrium points. However, the specifications for G may include non-equilibrium points. Applying reach to this problem to obtain stability would require the set of equilibrium points for arbitrary nonlinear dynamics, which is a challenging task in itself to find and may not even exist. In this work, we consider dynamics for which such a set is not known a priori. Hence, applying reach-avoid methods to the stabilize-avoid problem may not give desireable results. III. STABILIZE-AVOID AS INFINITE-HORIZON CONSTRAINED OPTIMAL CONTROL PROBLEM We tackle Problem 1 by solving an infinite-horizon con- strained OCP. To motivate this problem formulation, let l : X → R≥0 denote a non-negative cost function that takes zero value on G and is positive outside G, and define the infinite-horizon undiscounted policy value function V l,π for an arbitrary policy π as V l,π (x0) := ∞X k=0 l(xk), xk+1 = f xk, π(xk)
. (4) Since l is non-negative, V l,π is also non-negative. Using dynamic programming principles, we obtain V l,π (xk) = l(xk) + V l,π f (xk, π(xk))
. (5) The above equations are very close to satisfying the conditions for a discrete-time Lyapunov function [18]. By imposing additional assumptions on l, we can show that V l,π is a Lyapunov function. We leave the proof in the Appendix A1 for conciseness. Note that, for a given π, the set over which stability holds may be very small or even empty. This motivates us to ask whether this set can be maximized. The answer is affirmative here: solving the undiscounted infinite-horizon OCP min π ∞X k=0 l(xk) (6a) s.t. xk+1 = f (xk, π(xk)), (6b) gives a globally asymptotically stabilizing controller under some mild assumptions on the cost function and the control- lability of the dynamics [33]. This provides an answer to Problem 1 when safety constraints are not considered. However, solving an unconstrained infinite-horizon OCP does not guarantee satisfaction of the safety constraints xk̸ ∈ A. Hence, we consider solving a constrained infinite-horizon OCP to obtain a policy that is safe by construction. Let the superlevel set of h : X → R describes the avoid set A := { x : h(x) > 0 } . (7) We then solve the following constrained infinite-horizon OCP min π ∞X k=0 l(xk) (8a) s.t. xk+1 = f (xk, π(xk)), (8b) h(xk) ≤ 0, k ≥ 0. (8c) In the constrained setting, the proof from [33] that the optimal policy is globally asymptotically stabilizing is not applicable here. However, the proof can be extended to handle this case under certain conditions, which we leave as future work. IV. SAFETY CONSTRAINED OPTIMAL CONTROL VIA EPIGRAPHIC REFORMULATION The previous section describes how solving the stabilize- avoid problem (Problem 1) can be cast as solving an infinite- horizon constrained OCP (6). In this section, we propose solving the constrained OCP by reformulating the problem into its epigraph form and then solving the resulting two-stage optimization problem. A. Epigraph Form For any optimization problem of the form min x J(x) (9a) s.t. h(x) ≤ 0, (9b) the epigraph form [8, pp 134] of the above (9) is the optimization problem min x, z z (10a) s.t. h(x) ≤ 0, (10b) J(x) ≤ z. (10c) where z ∈ R is an auxiliary optimization variable. It is a standard result in optimization [8, pp 134] that (9) and (10) are equivalent. Now, observe that the constraints (10b) and (10c) can be combined to yield the following: min x, z z (11a) s.t. max{h(x), J(x) − z} ≤ 0. (11b) We can further move the minimization of the x variable into the constraint (11b) (see Appendix A2 for proof) to yield min z z (12a) s.t. min x max{h(x), J(x) − z} ≤ 0. (12b) This form allows us to convert the original constrained problem (9) into an unconstrained inner problem over x (12b) and a constrained outer problem (12) over the scalar decision variable z. At the optimal point (x∗, z∗), optimality conditions imply that z∗ = J(x∗). Solving for z∗ can thus be thought of as solving for the cost J at the optimal solution. Consequently, if we can bound the value of J(x∗), then z∗ will lie within the same bound. This facilitates treating z as a “cost budget” (with units of J) for satisfying h. As z → ∞ (i.e., the “cost budget” for J(x) increases), h(x) will dominate the max, and x∗ will focus on minimizing h more. On the other hand, as z → −∞ (i.e., the “cost budget” for J(x) decreases), J(x) − z will dominate the max, and x∗ will focus on minimizing J. B. Epigraph Form Constrained OCP We now apply this to the constrained OCP (8). First, we express the safety constraint (8c) equivalently as max k≥0 h(xk) ≤ 0. (13) Using this, the epigraph form of the constrained OCP (8) reads min z z (14a) s.t. ˜V (x0, z) ≤ 0, (14b) where the auxiliary value function ˜V is the OCP analogue of the LHS of (12b) ˜V (x0, z) := min π ˜Jπ (x0, z) s.t. xk+1 = f (xk, π(xk)), (15) with ˜Jπ defined as ˜Jπ (x0, z) := max ( max k≥0 h(xk), ∞X k=0 l(xk) − z ) . (16) We denote this the epigraph form constrained OCP (EF-COCP). C. Dynamic Programming for EF-COCP Note that (16) has both a maximization and a sum and hence has a different structure compared to the single sum in the objective function of the typical OCP. Consequently, the Bellman equation cannot be used in this case. We derive the corresponding dynamic programming equations below as ˜V (xk, zk) = min uk max n h(xk), V xk+1, zk+1
o , (17) where zk+1 has the following “dynamics” zk+1 = zk − l(xk). (18) This can again be understood from the intuition of z as a “cost budget” in J for satisfying the constraints h. Moving from timestep k to k + 1 incurs the cost l(xk), which is subtracted from the current “budget” zk to yield the next budget zk+1. If the “budget” zk falls low enough, the cost term in (16) will dominate the max, and we will have “run out of budget” to focus on constraint satisfaction. D. Relationship with Lagrangian Duality The epigraph form (11a) shares some similarities with the La- grangian duality formulation commonly used in (undiscounted) constrained MDPs max λ≥0 min π ∞X k=0 l(xk) + λ ∞X k=0 [h(xk)]+ | {z } :=L(π,λ) , (19) where [·]+ = max(0, ·), is used to disallow constraints viola- tion. Both formulations (19), (14) are two-stage optimization problems, where the outer problem consists of an extra scalar variable (z, λ respectively) while the inner problem optimizes with respect to the policy. The inner problem for both formulations is shown in Figure 3. We note the following two differences. Optimization Stability. Due to the [·]+ in (19), gradients of the inner problem with respect to λ will always be non-negative. Consequently, as long as the constraints are not satisfied, λ will continue to increase. However, large values of λ are problematic when constraints are not satisfied, since the gradients ∇xL with respect to x scale linearly in λ. In [45], a solution to this problem is proposed by rescaling L by 1/(1 + λ). However, when λ is large, the gradients for the l terms will instead vanish. Moreover, since λ is a non-decreasing function of the number of optimization iterations, this problem will only become worse as optimization proceeds. In contrast, since z is additive within the max, the scale of gradients is not altered. Consequently, EF-COCP does not suffer from this issue of optimization instability. 0.0 0.5 1.0 1.5 2.0 2.5 h(x) J(x) = J + h = 1.0 = 4.0 = 2ln(2)(a) Lagrangian Dual2 1 0 1 2 h(x) J(x)-z max(h, J z) z = 0 z = 2 z = 2 z = J = 1 z = 1 (b) Epigraph Form Fig. 3: Comparison of the inner subproblem for the Lagrangian dual formulation used in CMDP (left) and the epigraphic formulation (right) for the problem minx J(x) s.t. h(x) ≤ 0. Note that the gradients of the full objective (purple) scale with λ (left) but are unaffected in scale by z (right). Intuition. The auxiliary variable z in the epigraph form is in units of cost and represents a cost budget, as shown earlier. On the other hand, the Lagrange multiplier λ is a ratio representing the cost per unit constraint, but this is harder to interpret when cost and constraints cannot be easily compared. Consequently, it is much easier to estimate upper bounds for z. We take advantage of this to bound the range of z used for solving the inner problem, which we discuss in the next section. V. SOLVING EF-COCP WITH DEEP REINFORCEMENT LEARNING The previous section introduces a new epigraph form of the constrained optimal control problem, but does not provide a method of solving this formulation. In this section, we tackle this problem via reinforcement learning and introduce a framework for learning controllers for complex, nonlinear, potentially non-smooth dynamics. Given that the inner optimization problem of the epigraph form (15) still retains many similarities with the original problem, we choose to solve for the policy and the value function using reinforcement learning. Specifically, we use Proximal Policy Optimization (PPO) [39] but with modified definitions of the value function, advantage functions, returns, and generalized advantage estimator (GAE) [38]. While we have treated the policy π as a deterministic function, we will use a stochastic π for the purpose of improved exploration while performing Deep RL, where π(u|x) now defines a distribution over controls. However, while E[a + b] = E[a] + E[b], the same does not hold over the max operator used to define ˜Jπ in (16). Consequently, we need to be careful when defining ˜V π for a stochastic policy such that an analogous dynamic programming equation to the deterministic case (17) can be applied. Consider the following nested expectation form of an OCP policy value function: V π (x0) = lim K→∞ E0:K " KX k=0 lk # , (20a) = lim K→∞ E0  l0 + · · · + EK−1 h lK−1 + EK [lK ] i , (20b) = E0[l0 + V π (x1)], (20c) where, for conciseness, we denote Ek := Euk |xk , Ek:t := Euk ,...ut|xk , lk := l(xk). (21) Note how the nested expectations of (20b) lends itself to the dynamic programming equations of (20c). In the case of EF-COCP, we define ˜V π analogously to obtain ˜V π (x0, z0) (22a) = lim K→∞E0  h0 ∨ · · · ∨ EK−1 h hK−1 ∨ EK [hK ∨ KX k=0 lk − z0] i , (22b) = E0 h0 ∨ ˜V π (x1, z1), (22c) where we have used a ∨ b := max(a, b) for conciseness and where xk, zk follow the dynamics (1) and (18). We also define the action-value function ˜Qπ : X × R × U → R≥0 as ˜Qπ (xk, zk, uk) = max  h(xk), ˜V π (xk+1, zk+1)  , (23) such that ˜V π (xk, zk) = Euk  ˜Qπ (xk, zk, uk). We can now de- rive a policy gradient theorem for the inner problem. Theorem 1 (Policy Gradient Theorem). The gradient of the policy value function ˜V πθ (22) for the inner subproblem satisfies ∇θ ˜V πθ (x0, z0) ∝ E(x,z,u)1:k ∼πθ h ˜Qπθ (xk, zk, uk)∇θ ln πθ (uk|xk, zk) ξ = 1 i , (24) where the binary random variable ξk is defined to be equal to 1 when h(xt) ≤ ˜V πθ (xt+1, zt+1) is true for all t = 0, . . . , k. Proof: The proof follows from the proof of the normal policy gradient theorem [46], differing only in the expression for ∇θ ˜Qπθ . In the normal setting [46], (20c) gives ∇θ Qπθ (xk, uk) = ∇θ  l(xk)+V πθ (xk+1)  = ∇θ V πθ (xk+1). (25) In the case of EF-OCP with (22c), ∇θ ˜Qπθ (xk, uk) = ∇θ max n h(xk), ˜V πθ (xk+1, zk+1) o , (26) = 1h(xk )≤ ˜V πθ (xk+1,zk+1)∇θ ˜V πθ (xk+1, zk+1). (27) ˜V (x, z) = min π max { max k≥0 h(xk), ∞∑ k=1 l(xk) − z } Agent (Inner Problem) xk+1 = f (xk, uk) Environment0 1 2 z(x) 0.0 0.2 0.4 0.6 V z∗(x) = minz { z ∣ ∣ ∣ ˜V (x, z) ≤ 0 } z∗(x) (Outer Problem) π∗(x) = π(x, z∗(x)) Stabilize-Avoid Policyuk l(xk), h(xk) π(x, z) ˜V (x, z) EFPPOFig. 4: Summary of the EFPPO algorithm. First, reinforcement learning is used to solve the inner problem (15) and learn ˜V (x, z) and π(x, z) over the entire state space. Then, the optimal z∗ which solves the outer problem (14) is regressed. Following the rest of the normal proof [46] then yields (24) From Theorem 1, we can construct a basic on-policy DeepRL algorithm to solve the inner problem (15) over all states x parametrized by a range of z, yielding a learned (stochastic) controller πθ (x, z). However, we can do better by performing variance reduction via subtracting the baseline ˜V πθ (x, z) to get the advantage ˜Aπθ (x, z, u) := ˜Qπθ (x, z, u) − ˜V (x, z), as Eπθ [∇θ ln πθ ] = 0. Following PPO, we also apply the GAE estimate [38], perform clipped importance sampling and add in an entropy bonus to arrive at an algorithm that is very similar to PPO [39] but with ˜Qπ and ˜V π defined as above. Stochastic policy considerations: While we have performed the above developments using a stochastic policy and derived a stochastic policy gradient theorem, our desired solution to the inner optimization problem is a deterministic controller. Hence, we only take the mode of the learned policy, and treat the stochasticity purely as a means of performing exploration, discouraging premature convergence to local minima and smoothing the optimization landscape [3, 19, 31]. We also fine-tune the learned value function at the end of the DeepRL training by freezing the obtained deterministic policy and performing policy evaluation. It can be appealing to consider using a deterministic policy and apply a deterministic policy gradient theorem [40]. However, the proof of this theorem requires the transition distribution p(xt+1|xt, ut) to be continuous [40]. This does not hold in our problem since the transition distribution is degenerate (and hence discontinuous) due to the use of deterministic dynamics, i.e., p(xt+1|xt, ut) = δxt+1 − f (xt, ut)
, (28) preventing the use of deterministic policy gradient in this case. After obtaining ˜V π and π for the deterministic policy, we now turn to the outer problem (14). Since the outer optimization problem is only 1 dimensional, we can solve for z∗ easily via classical scalar optimization methods such as the bisection method which run quickly. Instead of running bisection online, however, we choose to learn the optimal z∗ : X → R offline by using the result of bisection as the label for a regression problem. This gives us z∗(x), which in turn provides the optimal policy for the original constrained optimal control problem. Algorithm 1 EFPPO Inner Problem input: Estimate of maximum cost upper bound zmax repeat if reset environments then Sample random x0 ∈ X and z0 ∈ [0, zmax] end if Run policy π in environments for T timesteps Compute advantage estimates Aπ with GAE Update policy π and baseline ˜V π with PPO clipping and entropy until converged Fine-tune ˜V π using the mode of π via policy evaluation. Algorithm 2 EFPPO Outer Problem input: Estimate of maximum cost upper bound zmax Sample dataset of x randomly from the state space Bisect ˜V π (x, ·) over [0, zmax] to obtain labels z∗ for each x repeat Train network ˜z∗(x) to predict z∗ given x via regression until converged The proposed EFPPO algorithm is summarized in Algo- rithms 1 and 2 and illustrated in Figure 4. EFPPO solves the two-stage optimization problem of EF-COCP (14a) sequen- tially. The inner problem Equation (15) is solved via policy gradient using Theorem 1 for a range of z values and uses the improvements in PPO such as GAE estimates, clipped importance sampling, and an entropy bonus. Next, we extract the mode of the stochastic π and fine-tune ˜V π . Then, we fix the value function ˜V π and the deterministic policy π(x, z) and learn z∗(x) by randomly sampling states in the state space and minimizing the residual to the analytical solution of (14) found by applying the bisection method. Since this is a 1D optimization problem, the bisection method converges to almost machine precision within tens of iterations. The final policy is then obtained as π(x, z∗(x)). When solving the inner problem, we randomly sample states from the state space and random sample z within the range [0, zmax], where zmax is an upper- bound estimate of the total cost P∞ k=0 l(xk). To prevent using zmax = ∞ since this term may be unbounded, we introduce a small discount factor only in the total costs. More details on the effect of discounting are available in the Appendix C. VI. EXPERIMENTS To evaluate the performance of the proposed EFPPO al- gorithm, we compare EFPPO against related algorithms on simulated stabilize-avoid problems with increasing complexity. The last problem involves a non-differentiable, nonconvex, non-control-affine system and demonstrates the ability of our approach to both maintain safety and successfully stabilize the system within the goal region even in nontrivial high- dimensional environments. We compare EFPPO against the following baselines methods. • PPO [39], a popular on-policy unconstrained DeepRL method. Despite its unconstrained nature, it is common to apply “soft constraints” in the form of penalties incurred when an undesirable state is reached [9, 51]. Unlike constrained methods such as EFPPO, the scale of the penalties is a hyperparameter that must be tuned to trade- off between constraint violation and training stability. We denote by PPO(λ) the family of methods where λ denotes the penalty scale of the modified cost function ˜l, i.e., ˜l(x) = l(x) + λ[h(x)]+, (29) where PPO(0) solves the unconstrained problem. • CPPO [45], a representative algorithm among the family of constrained DeepRL methods which use Lagrangian duality, and is an improvement on PPO-Lagrangian from [1]. As noted in Section IV-D, the constrained MDP formulation allows constraint violations up to a cost thresh- old. To learn policies that strictly satisfy the constraints, we set the cost threshold to 0. We perform a manual hyperparameter search to select the PID parameters. • PPO-SIS [26], a constrained RL method that applies the safety constraint at each state instead of in expectation as in CMDP. Moreover, a safety certificate is learned jointly to improve the safety of the learned policy. • CLBF1 [14], which learns a Lyapunov function using a neural network via losses which penalize violations of the Lyapunov conditions. Safety can be guaranteed by enforcing that a sublevel set of the learned neural Lyapunov function lies outside the avoid set. Given a learned CLBF, we consider two methods for synthesizing the controller: (1) CLBF(QP) solves the CLBF-QP using the CLBF as a constraint; (2) CLBF(Opt) applies the (bang-bang) control that minimizes ˙V . As noted previously, while the reach-avoid problem is closely related to the stabilize-avoid problem considered in this work, the policies obtained from reach-avoid do not induce stability when the goal set is not an equilibrium point (see Section II-A). For a fair comparison, a feedforward neural network with tanh activations is used for both the policy (if used) and value functions. Additional networks are also defined using the same architecture above but with modified final activation function (e.g., softplus) when used in the original implementation. 1While the original work considers CLBFs that are robust to parametric uncertainties, the problems we consider here do not have parameters. Safety Rate ↑ Cost ↓ Stabilize Rate ↑ PPO(0) 0.55 1.251 1.0 PPO(1) 0.61 1.283 1.0 PPO(10) 1.0 1.299 1.0 CPPO [45] 1.0 1.314 1.0 PPO-SIS [26] 1.0 2.227 0.182 CLBF(QP) [14] 0.96 3.042 0.180 CLBF(Opt) [14] 0.99 3.180 0.035 EFPPO (Ours) 1.0 1.285 1.0 TABLE I: Comparison of controller performance on the double integrator example. Metrics are computed from states randomly sampled from the control-invariant set. To compare each method, we use the safety rate, cost and stabilize rate metrics computed by rolling out the learned policy on randomly sampled initial states. When the true control- invariant set that can guarantee safety is analytically known, the sampled states are sampled from the control-invariant set. However, for systems with complex dynamics where such a set cannot be found analytically, we sample from a crude box estimate. The safety rate and stabilize rate are defined as the fraction of trajectories that satisfy the safety constraints and can reach and stay within the goal set for the last 50 timesteps respectively. The top method for each metric is shown in bold. For the cost, we only highlight the best performing method among the methods that handle constraints (i.e., unconstrained PPO is not taken into account). More details regarding the specific dynamics and constraints for each task are provided in the Appendix D for brevity. A. 1D Double-Integrator We first consider a simple double-integrator dynamics in 1D, where the optimal policy and corresponding optimal control- invariant region for the problem can be computed. Given states [p, v] ∈ R2, controls u = a ∈ R the task here is to stabilize to the region pgoal := [0.65, 0.85] while satisfying box constraints on both states and controls |p| ≤ 1, |v| ≤ 1, |u| ≤ 1. (30) While all constrained methods are able to maintain safety for all states within the control-invariant region, only EFPPO is able to also stabilize all trajectories to the goal set G with a cost similar to the unconstrained solution from PPO(0) (see Table I). PPO-SIS has regions of the state space outside G which are equilibrium points (see Figure 5), while CPPO has a much larger cost compared to the unconstrained solution. The unconstrained PPO(λ) methods require the penalty weight λ to be large enough for a safe policy, but this trade-off comes at the cost of the policy’s performance and requires careful tuning of the penalty weight λ. In contrast, EFPPO is able to synthesizes a performant safe controller without requiring any cost function tuning. PPO(0) PPO(1) PPO(10) CPPO SIS CLBF(QP) EFPPOFig. 5: Trajectory rollouts (• → •) on the double-integrator system with box constraints on both position and velocity (shown in grey) for the avoid set. Unsafe trajectories are shown in red. For CLBF(QP), a contour plot of the learned CLBF is shown with red regions denoting higher values and the safe level set shown as the blue region. Due to control constraints |u| ≤ 1, the control-invariant region is smaller than the complement of the avoid set.PPO(0) PPO(1) PPO(10) CPPO SIS CLBF EFPPO Fig. 6: Trajectory rollouts (• → •) on the single-integrator system in 2D with sector constraints. For CLBF, a contour plot of the learned CLBF is shown. Safety Rate ↑ Cost ↓ Stabilize Rate ↑ PPO(0) 1.0 0.363 1.0 PPO(1) 1.0 0.364 1.0 PPO(10) 1.0 0.459 1.0 CPPO [45] 1.0 0.365 1.0 PPO-SIS [26] 1.0 2.626 0.0 CLBF(QP) [14] 1.0 12.613 0.0 CLBF(Opt) [14] 1.0 1.308 0.0 EFPPO (Ours) 1.0 0.363 1.0 TABLE II: Comparison of controller performance on the 2D single integrator with sector obstacles. Metrics are computed from states randomly sampled from the control-invariant set. B. 2D Single-Integrator with Sector Obstacle Next, we consider a single-integrator in 2D where the avoid set is defined as a sector of the circle and the goal set defined in the center. The state is defined as the positions [px, py ] ∈ R2 with controls denoting the velocities [vn, vt] ∈ [−1, 1]2 ⊂ R2, where vn and vt denote the normal and tangent components of the velocity vector to the center of the circle. Table II summarizes the results in this task. Due to the control parametrization and the shape of the avoid set, the optimal control at all states is the constant vector [1, 0]. Of the constrained methods, only EFPPO is able to learn this, with other methods learning suboptimal versions of the optimal policy. Note that the rollouts for CLBF [14] actually move away from the center when starting near the sector obstacle on the right (see Figure 6). This is because the CLBF learns a distorted metric near the obstacles (as seen from the value function) to ensure the level-set is contained within the control- invariant set. Moreover, we see that SIS also fails to reliably stabilize on this simple example. The safety index used in SIS [26] is taken from [54] and consists of only three parameters. We suspect that this is insufficient for tasks where the avoid set A is relatively complex. Finally, we see that PPO(10) learns a particularly poor performing controller due to the large penalty weights. In particular, the difference in scales between the large costs in the unsafe regions and the small costs near the goal set poses a challenge for learning the value function accurately, resulting in the learned policy putting more weight on constraint satisfaction. C. Hopper Stabilization The preceding benchmarks are low-dimensional toy- examples that provide intuition on how the different methods be- have for simple linear systems. Our next example demonstrates the ability of EFPPO to stabilize to goal sets while maintaining safety for more complex dynamics. We consider the hopper in the Brax simulator [15], a nonlinear non-differentiable system with a 12 dimensional state space and a 3 dimensional control space. Note that the original dynamics only consider a 11 dimensional state space as they discard the x-coordinate, which we keep. The goal set G in this task is defined as G := { x | px ∈ [2.8, 3.0] } , (31) while the constraints limit the height and rotation of the hopper’s torso pz ≥ 0.7, |θ| ≤ 0.2. (32) Unlike the normal setup for Hopper, stabilization to G requires keeping track of the x-position. Consequently, the optimal policy is no longer a limit-cycle and even requires the hopper to move backwards for some initial states. We report our results in Table III. On this problem, we see that both CPPO (a) PPO(0) (b) PPO(1) (c) PPO(10) (d) EFPPO (e) CPPO Fig. 7: The hopper task asks for stabilization of the hopper’s torso to the set px ∈ [2.8, 3.0] (dashed green lines, drawn when the hopper does not reach the goal set). Two rollouts are shown starting from px = 1.0 (blue) and px = 4.0 (red) respectively, with time progressing from opaque to transparent. Only EFPPO safely stabilizes the system to the goal set in these two rollouts. Safety Rate ↑ Cost ↓ Stabilize Rate ↑Safe All PPO(0) 0.000 1.360 0.893 PPO(1) 0.676 2.986 3.806 0.661 PPO(10) 0.037 15.193 15.503 0.087 CPPO [45] 0.724 8.353 9.458 0.087 PPO-SIS [26] 0.000 9.703 0.084 EFPPO (Ours) 0.833 1.568 3.695 0.843 TABLE III: Comparison of controller performance on the Hopper system on a set of random initial states that may lie outside the (unknown) optimal control-invariant set. and EFPPO achieve high safety rates, while PPO-SIS fails to stabilize. We see the same trend of PPO(λ) presenting a trade- off between stability and safety, except for PPO(10) which is more unsafe than PPO(1). We suspect this is due to larger costs destabilizing training. D. F16 Ground Collision Avoidance In a Low Altitude Flight Corridor Finally, we showcase the scalability of our method on a ground collision avoidance example involving the F16 fighter jet [20]. This system is non control-affine, non-smooth and involve lookup tables, making this a challenging system to solve stabilize-avoid on in addition to the high 17-dimensional state space and 4-dimensional control space. The task here is to stabilize the F16 to a target altitude defined by the set [50, 150]ft under box control constraints while staying within a flight corridor heading North (box constraints on the East-Up plane). For positions (pE , pN , pU ), this corresponds to −200 ft ≤ pE ≤ 200 ft, 0 ft ≤ pU ≤ 1000 ft. (33) The results are summarized in Table IV and Figure 8. Following the trend from before, we see that CPPO is able to achieve high safety rates at the expense of being unable to stabilize to the goal set. In contrast, EFPPO using the optimal z∗ achieves similar safety rates as CPPO but has a 70-fold increase in the stabilize rate. Moreover, note that EFPPO(0) and EFPPO(1.8) denote using a constant value of z for the policies from the inner EF-COCP problem (15). As expected, using a larger value of z results in improved constraint satisfaction at the cost of a lower stabilize rate. By using the optimal value of z∗, EFPPO(z∗) is able to achieve the high stabilize rate of EFPPO(0) while Safety Rate ↑ Cost ↓ Stabilize Rate ↑Safe All PPO(0) 0.000 0.843 0.917 PPO(1) 0.657 0.639 1.025 0.742 PPO(10) 0.000 4.032 0.092 CPPO [45] 0.827 4.874 4.822 0.012 EFPPO(0) (Ours) 0.792 0.735 1.102 0.837 EFPPO(1.8) (Ours) 0.856 2.806 2.945 0.022 EFPPO(z∗) (Ours) 0.823 0.724 1.139 0.843 TABLE IV: Comparison of controller performance on the F16 system on a set of random initial states that may lie outside the (unknown) optimal control-invariant set. maintaining the high safety rates of EFPPO(1.8). However, note that the safety rate of EFPPO(z∗) is actually slightly lower than that of EFPPO(1.8). This is because z∗ is approximated via a neural network which may have small approximation errors. We also see the same trends for PPO(λ) from Hopper carry over to this task. VII. DISCUSSION Many Lagrangian duality-based methods suffer from fo- cusing too much on satisfying the safety constraints and consequently suffer in terms of their stability. As noted in Section IV-D, the non-negative gradients of the Lagrange multipliers λ means that they will continue to increase as long as constraints are not being satisfied. However, a large λ causes the optimization problem to become badly conditioned, hindering the performance of these algorithms. Moreover, in systems with complex dynamics where the control-invariant set is not known, it may be impossible for constraints satisfaction to occur in all states. In this case, λ will grow unbounded until either the training algorithm becomes unstable or it hits a user-defined maximum limit. This is a failure that the PID mechanism in CPPO does not fix, as λ does not oscillate. We found all three coefficients of the PID played a similar role in adjusting the rate at which λ (monotonically) increases and eventually destabilizes training. Although the CMDP problem formulation can be adapted to solve stabilize-avoid problems, these methods are a poor fit. Additionally, we observe that the CLBF performs poorly even on both toy examples. We believe there are two reasons for this. (a) EF-PPO (b) CPPO Fig. 8: Ground collision avoidance and stabilization to a low altitude flight corridor for the F16 fighter jet system. 1) Learning of the CLBF assumes that the control-invariant set is known (or can be sampled from) a-priori. However, in the experiments considered in this work, we do not assume knowledge of these sets. While we do give CLBF a small region around the equilibrium subset of the goal region, information that was not available to other algorithms, this was not enough for CLBFs to perform well. 2) The learning problem formulated in [14] is under- parametrized. More specifically, given any CLBF V which satisfies the CLBF conditions, αV will induce the same set of feasible constraints and will perform identically to V . However, if V does not satisfy the CLBF conditions, the violation error can be reduced arbitrarily by taking α → 0. This yields a value function where regions close to the goal set have near-zero gradients which may violate the CLBF conditions and compromise the controller’s performance, which we have observed empirically in our experiments. In contrast, the proposed EFPPO algorithm can perform well even when large areas of the state space are not control- invariant as in our experiments. This is due to two reasons: (a) z does not affect the magnitude of the gradients directly, unlike λ which scales the gradients of the constraint (and also the cost, if the reweighing scheme from [45] is used). Consequently, while z → ∞ will make no difference to EFPPO (if the z feature variable is normalized correctly), taking λ → ∞ will cause the Lagrangian to diverge and cause training instabilities. (b) Since z can be interpreted as a “cost budget” and it is easier to estimate an upper-bound on costs, we can afford to solve the two-stage optimization problem sequentially. In contrast, it is much harder to bound the optimal λ since the units of “cost to constraint ratio” is more difficult to reason about. VIII. CONCLUSION We present a new method for synthesizing nonlinear feedback controllers for performing stabilization while maintaining safety under control constraints. By formulating the stabilize- avoid problem as an infinite-horizon epigraph-form constrained optimal control problem and applying deep reinforcement learning, our approach is able to sidestep numerical challenges that other methods face and achieve vastly larger regions of attraction while still maintaining safety for high dimensional complex systems. Limitations and future work: The EFPPO algorithm currently splits the task of learning V, π and z∗ into separate stages of optimization, and relies on random sampling of z in the first stage to cover the state-space. However, given the structure of the epigraph form constrained OCP, it should be possible to perform both optimizations simultaneously such that only a single stage of optimization is required. Moreover, the current method does not account for model errors and parametric uncertainties which may compromise the safety and stability of the learned controllers. An extension of the current method to consider this would allow for more robustness when deploying such a controller to the hardware systems in the real world. Additionally, PPO is an on-policy online reinforcement learning algorithm. Consequently, a simulator of the dynamics is necessary to solve the EF-COCP problem. Extending this work to off-policy offline reinforcement learning setting will allow our method to be applied to settings when a dynamics simulator is not available. Finally, it is difficult to provide useful statements about the convergence of practical DeepRL algorithms. Nevertheless, it is important to understand properties of the learned policy, especially as violations of safety constraints are serious and undesirable in any system. We leave this as future work. ACKNOWLEDGMENTS This work is partially supported by the MIT Lincoln Lab under the Safety in Aerobatic Flight Regimes (SAFR) program. However, this article solely reflects the opinions and conclusions of its authors and not the MIT Lincoln Lab. REFERENCES [1] Joshua Achiam, David Held, Aviv Tamar, and Pieter Abbeel. Constrained policy optimization. In International conference on machine learning, pages 22–31. PMLR, 2017. [2] Amir Ali Ahmadi and Anirudha Majumdar. Some applications of polynomial optimization in operations research and real-time decision making. Optimization Letters, 10:709–729, 2016. [3] Zafarali Ahmed, Nicolas Le Roux, Mohammad Norouzi, and Dale Schuurmans. Understanding the impact of en- tropy on policy optimization. In International conference on machine learning, pages 151–160. PMLR, 2019. [4] Albert Altarovici, Olivier Bokanowski, and Hasnaa Zidani. A general hamilton-jacobi framework for non-linear state-constrained control problems. ESAIM: Control, Optimisation and Calculus of Variations, 19(2):337–357, 2013. [5] Eitan Altman. Constrained Markov decision processes: stochastic modeling. Routledge, 1999. [6] Aaron D Ames, Xiangru Xu, Jessy W Grizzle, and Paulo Tabuada. Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control, 62(8):3861–3876, 2016. [7] Zvi Artstein. Stabilization with relaxed controls. Non- linear Analysis: Theory, Methods & Applications, 7(11): 1163–1173, 1983. [8] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004. [9] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba. Openai gym. arXiv preprint arXiv:1606.01540, 2016. [10] Italo Capuzzo-Dolcetta and P-L Lions. Hamilton-jacobi equations with state constraints. Transactions of the American mathematical society, 318(2):643–683, 1990. [11] Ya-Chien Chang, Nima Roohi, and Sicun Gao. Neural lyapunov control. Advances in neural information processing systems, 32, 2019. [12] Hongkai Dai and Frank Permenter. Convex synthesis and verification of control-lyapunov and barrier functions with input constraints. arXiv preprint arXiv:2210.00629, 2022. [13] Charles Dawson, Sicun Gao, and Chuchu Fan. Safe control with learned certificates: A survey of neural lyapunov, barrier, and contraction methods. arXiv preprint arXiv:2202.11762, 2022. [14] Charles Dawson, Zengyi Qin, Sicun Gao, and Chuchu Fan. Safe nonlinear control using robust neural lyapunov- barrier functions. In Conference on Robot Learning, pages 1724–1735. PMLR, 2022. [15] C. Daniel Freeman, Erik Frey, Anton Raichuk, Sertan Girgin, Igor Mordatch, and Olivier Bachem. Brax - a differentiable physics engine for large scale rigid body simulation, 2021. URL http://github.com/google/brax. [16] Peter Giesl and Sigurdur Hafstein. Review on compu- tational methods for lyapunov functions. Discrete & Continuous Dynamical Systems-B, 20(8):2291, 2015. [17] Gene Grimm, Michael J Messina, Sezai Emre Tuna, and Andrew R Teel. Model predictive control: for want of a local control lyapunov function, all is not lost. IEEE Transactions on Automatic Control, 50(5):546–558, 2005. [18] Lars Grüne, Jürgen Pannek, Lars Grüne, and Jürgen Pannek. Nonlinear model predictive control. Springer, 2017. [19] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. In International conference on machine learning, pages 1861– 1870. PMLR, 2018. [20] Peter Heidlauf, Alexander Collins, Michael Bolender, and Stanley Bak. Verification challenges in f-16 ground collision avoidance and other automated maneuvers. In ARCH@ ADHS, pages 208–217, 2018. [21] Kai-Chieh Hsu, Vicenç Rubies-Royo, Claire J Tomlin, and Jaime F Fisac. Safety and liveness guarantees through reach-avoid reinforcement learning. arXiv preprint arXiv:2112.12288, 2021. [22] Mrdjan Jankovic. Robust control barrier functions for constrained stabilization of nonlinear systems. Automatica, 96:359–367, 2018. [23] Christopher Jewison and R Scott Erwin. A spacecraft benchmark problem for hybrid control and estimation. In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 3300–3305. IEEE, 2016. [24] Paola Loreti. Some properties of constrained viscosity solutions of hamilton–jacobi–bellman equations. SIAM journal on control and optimization, 25(5):1244–1252, 1987. [25] Haitong Ma, Yang Guan, Shegnbo Eben Li, Xiangteng Zhang, Sifa Zheng, and Jianyu Chen. Feasible actor-critic: Constrained reinforcement learning for ensuring statewise safety. arXiv preprint arXiv:2105.10682, 2021. [26] Haitong Ma, Changliu Liu, Shengbo Eben Li, Sifa Zheng, and Jianyu Chen. Joint synthesis of safety certificate and safe control policy using constrained reinforcement learn- ing. In Learning for Dynamics and Control Conference, pages 97–109. PMLR, 2022. [27] Kostas Margellos and John Lygeros. Hamilton–jacobi formulation for reach–avoid differential games. IEEE Transactions on automatic control, 56(8):1849–1861, 2011. [28] David Mayne. An apologia for stabilising terminal conditions in model predictive control. International Journal of Control, 86(11):2090–2095, 2013. [29] Hiroyoshi Mitake. Asymptotic solutions of hamilton- jacobi equations with state constraints. Applied Mathe- matics and Optimization, 58:393–410, 2008. [30] Ian M Mitchell. The flexible, extensible and efficient toolbox of level set methods. Journal of Scientific Computing, 35:300–329, 2008. [31] Volodymyr Mnih, Adria Puigdomenech Badia, Mehdi Mirza, Alex Graves, Timothy Lillicrap, Tim Harley, David Silver, and Koray Kavukcuoglu. Asynchronous methods for deep reinforcement learning. In International con- ference on machine learning, pages 1928–1937. PMLR, 2016. [32] Frank Permenter and Pablo Parrilo. Partial facial reduc- tion: simplified, equivalent sdps via approximations of the psd cone. Mathematical Programming, 171:1–54, 2018. [33] Romain Postoyan, Lucian Bu¸soniu, Dragan Neši´c, and Jamal Daafouz. Stability analysis of discrete-time infinite- horizon optimal control with discounted cost. IEEE Transactions on Automatic Control, 62(6):2736–2749, 2016. [34] Zengyi Qin, Kaiqing Zhang, Yuxiao Chen, Jingkai Chen, and Chuchu Fan. Learning safe multi-agent control with decentralized neural barrier certificates. arXiv preprint arXiv:2101.05436, 2021. [35] Sasa V Rakovic and William S Levine. Handbook of model predictive control. Springer, 2018. [36] Matheus F Reis, A Pedro Aguiar, and Paulo Tabuada. Control barrier function-based quadratic programs intro- duce undesirable asymptotically stable equilibria. IEEE Control Systems Letters, 5(2):731–736, 2020. [37] Muhammad Zakiyullah Romdlony and Bayu Jayaward- hana. Stabilization with guaranteed safety using control lyapunov–barrier function. Automatica, 66:39–47, 2016. [38] John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel. High-dimensional continuous control using generalized advantage estimation. arXiv preprint arXiv:1506.02438, 2015. [39] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017. [40] David Silver, Guy Lever, Nicolas Heess, Thomas Degris, Daan Wierstra, and Martin Riedmiller. Deterministic policy gradient algorithms. In International conference on machine learning, pages 387–395. Pmlr, 2014. [41] Halil Mete Soner. Optimal control with state-space constraint i. SIAM Journal on Control and Optimization, 24(3):552–561, 1986. [42] Halil Mete Soner. Optimal control with state-space constraint. ii. SIAM journal on control and optimization, 24(6):1110–1122, 1986. [43] Eduardo D Sontag. A lyapunov-like characterization of asymptotic controllability. SIAM journal on control and optimization, 21(3):462–471, 1983. [44] Brian L. Stevens and Frank L. Lewis. Aircraft Control and Simulation. Aircraft Engineering and Aerospace Technology, 76(5), January 2004. ISSN 0002-2667. doi: 10.1108/aeat.2004.12776eae.001. URL https://doi.org/10. 1108/aeat.2004.12776eae.001. Publisher: Emerald Group Publishing Limited. [45] Adam Stooke, Joshua Achiam, and Pieter Abbeel. Respon- sive safety in reinforcement learning by pid lagrangian methods. In International Conference on Machine Learning, pages 9133–9143. PMLR, 2020. [46] Richard S Sutton and Andrew G Barto. Reinforcement learning: An introduction. MIT press, 2018. [47] Weehong Tan and Andrew Packard. Searching for control lyapunov functions using sums of squares programming. sibi, 1(1), 2004. [48] Chen Tessler, Daniel J Mankowitz, and Shie Mannor. Reward constrained policy optimization. arXiv preprint arXiv:1805.11074, 2018. [49] Emanuel Todorov, Tom Erez, and Yuval Tassa. Mujoco: A physics engine for model-based control. In 2012 IEEE/RSJ international conference on intelligent robots and systems, pages 5026–5033. IEEE, 2012. [50] Claire J Tomlin, John Lygeros, and S Shankar Sastry. A game theoretic approach to controller design for hybrid systems. Proceedings of the IEEE, 88(7):949–970, 2000. [51] Zhuang Wang, Hui Li, Zhaoxin Wu, and Haolin Wu. A pretrained proximal policy optimization algorithm with reward shaping for aircraft guidance to a moving destination in three-dimensional continuous space. In- ternational Journal of Advanced Robotic Systems, 18(1): 1729881421989546, 2021. [52] Peter Wieland and Frank Allgöwer. Constructive safety us- ing control barrier functions. IFAC Proceedings Volumes, 40(12):462–467, 2007. [53] Shuowen Zhang, Yong Zeng, and Rui Zhang. Cellular- enabled uav communication: A connectivity-constrained trajectory optimization perspective. IEEE Transactions on Communications, 67(3):2580–2604, 2018. [54] Weiye Zhao, Tairan He, and Changliu Liu. Model-free safe control for zero-violation reinforcement learning. In 5th Annual Conference on Robot Learning, 2021. APPENDIX A PROOFS A1 Proof that V l,π is a discrete-time Lyapunov function Before we begin the statement of the theorem and its proof, define K to be the class of functions γ : R≥0 → R≥0 that is continuous, zero at zero and strictly increasing. Let K∞ to be the class that are K and are also unbounded. We then have the following definition of a Lyapunov function for discrete-time systems, which we modify from [18, Ch.2] to use σ(·) instead of ∥·∥xref . Definition 1 (Discrete-Time Lyapunov Function). Suppose the system has discrete-time dynamics xk+1 = f (xk, uk), (A.1) for states xk ∈ X and controls uk ∈ U. Consider a reference set G and a subset of the state space R ⊆ X . Let σ : X → R≥0 be a state measure (as in [17]) that is continuous and positive-definite. A function V : R → R≥0 is a uniform Lyapunov function on R if the following conditions are satisfied. (i) There exist functions α1, α2 ∈ K∞ such that α1(σ(x)) ≤ V (x) ≤ α2(σ(x)) (A.2) holds for all x ∈ R. (ii) There exists a function αV ∈ K such that V (xk+1) ≤ V (xk) − αV (σ(x)) (A.3) holds for all xk ∈ R. Theorem 2 (Policy Value Function is Lyapunov). Let π : X → U be an arbitrary deterministic policy, and define V l,π : X → R≥0 ∪ {+∞} to be the policy value function V l,π (x0) := ∞X k=0 l(xk), xk+1 = f (xk, π(xk)) (A.4) for cost function l : X → R≥0 and discrete dynamics f : X × U → X . Let F denote the set where V l,π is finite, i.e., F :=  x V l,π (x) < ∞ , (A.5) and let σ : X → R≥0 be a state measure (as in [17]) that is continuous and positive-definite. Suppose that the following holds for the cost function l and the policy value function V l,π . (i) There exists α ∈ K∞ such that, for any x ∈ F, V l,π (x) ≤ α(σ(x)) (A.6) (ii) There exists ρ ∈ K∞ such that, for any x ∈ F, l(x) ≥ ρ(σ(x)) (A.7) Then, V l,π is a Lyapunov function on F. Proof. By the definition of V l,π , use of dynamic programming shows that V l,π (xk) = l(xk) + V l,π (xk+1). (A.8) Since V l,π ≥ 0 by definition (A.4), by using (A.7) and (A.8) we can conservatively lower bound V l,π in terms of σ on F as V l,π (x) ≥ l(x) ≥ ρ(σ(x)). (A.9) Combining the same two equations again without dropping V l,π (xk+1), we can also show that for x ∈ F, V l,π (xk+1) = V l,π (xk) − l(xk) ≤ V l,π (xk) − ρ(σ(x)). (A.10) Combining (A.6), (A.9) and (A.10) then gives us that for xk ∈ F, ρ(σ(x)) ≤ V l,π (x) ≤ α(σ(x)), (A.11a) V l,π (xk+1) ≤ V l,π (xk) − ρ(σ(x)). (A.11b) Since ρ ∈ K∞, α ∈ K∞, (A.11a) and (A.11b) thus show that V l,π is a Lyapunov function on F by Definition 1. From Theorem 3, we can then apply the standard proof of local asymptotic stability using Lyapunov functions [18] to show asymptotic stability. Corollary 1. Define the set Z to be the zero-set of V l,π , i.e., Z :=  x V l,π (x) = 0 . (A.12) Then, Z is also the zero-set of σ, i.e., Z = { x | σ(x) = 0 } . (A.13) Moreover, Z is locally asymptotically stable within F under the controller π on F. Proof. First, note that by the definition of K, (A.9) implies that for x ∈ Z, 0 ≤ ρ(x) ≤ V l,π (x) = 0. (A.14) Moreover, since ρ is strictly increasing, ρ(x) > 0 =⇒ V l,π (x) > 0. (A.15) Hence, the zero-set Z of V l,π is also the zero-set of σ. Applying Theorem 2.19 from [18] using the policy value function as the Lyapunov function as shown in Theorem 3 then gives us the result. Note. As noted in the main paper, while Theorem 3 and Corollary 1 show that we can use V l,π to show stability for any policy π within the region F under assumptions (A.7) and (A.6), we note that F may be a tiny set or even empty. Hence, the theorems above do not give us a direct method of constructing stable controllers. Nevertheless, the above theorems provide intuition on the relationship between the optimality of a policy (measured by the size of F) and its stability, which we use when solving the infinite-horizon constrained OCP in the main paper. A2 Equivalence of (11) and (12) Theorem 3. Let x ∈ Rn and z ∈ R, and let g : Rn × R → R be a continuous (potentially non-differentiable) function. Then, if a solution exists (i.e., an optimal x∗, z∗ exist, are finite), then the following optimization problems are equivalent. min x, z z s.t. g(x, z) ≤ 0, (A.16) min z z s.t. h min x g(x, z) i ≤ 0, (A.17) Proof. We begin by comparing the Lagrangian primal problem of eq. (A.16) and eq. (A.17). min z min x max λ≥0 z + λ g(x, z) = min z  z + min x max λ≥0 λ g(x, z)  (A.18) min z max λ≥0 z + λ h min x g(x, z) i = min z  z + max λ≥0 min x λ g(x, z)  (A.19) Comparing the two, the only difference is that the order of minx and maxλ are flipped. Hence, it is sufficient to show that, for any z where minx g(x, z) < 0, p∗ := min x max λ≥0 λ g(x, z) = max λ≥0 min x λ g(x, z) =: d∗. (A.20) Note that this is exactly equivalent to showing that strong duality holds for the following constraint satisfaction problem. min x 0 s.t. g(x, z) ≤ 0, (A.21) We now prove that strong duality holds for the above problem in a similar fashion to the proof that Slater’s condition is a sufficient condition for strong duality to hold in convex optimization problems [8]. Define the set A ⊆ Rn × R as A := { (u, t) | ∃x, g(x, z) ≤ u, 0 ≤ t } , (A.22) = n u inf x g(x, z) ≤ u o × [0, ∞). (A.23) Note that A is convex. Furthermore, since a feasible solution exists by assumption, we have that p∗ = min x max λ≥0 λ g(x, z) = min x (∞ g(x, z) > 0 0 g(x, z) ≤ 0 = 0. (A.24) We now define a second set B ⊆ Rn × R as B := { (0, s) | s < p∗ } , (A.25) = {0} × (−∞, 0). (A.26) Note that B is also convex, and that the sets A and B do not intersect. We can then invoke the separating hyperplane theorem to show that there exists a (˜λ, μ)̸ = 0 and a α that defines a hyperplane which separates the two sets, i.e., (u, t) ∈ A =⇒ ˜λu + μt ≥ α (A.27) (u, s) ∈ B =⇒ ˜λu + μs ≤ α (A.28) In (A.27), since both u and t are unbounded above, we must have ˜λ ≥ 0 and μ ≥ 0. Furthermore, in (A.28), since s < p∗, we have that μp∗ ≤ α. Combining both then gives us that for all x, 0 = p∗ = μp∗ ≤ α ≤ ˜λg(x, z). (A.29) Minimizing the RHS over x then maximizing over ˜λ then gives us that p∗ ≤ min x ˜λg(x, z) ≤ max λ min x λg(x, z) = d∗. (A.30) Finally, by weak duality, we have that p∗ ≥ d∗. (A.31) Combining the two then allows us to conclude that p∗ = d∗. APPENDIX B UNDERSTANDING THE ROLE OF z IN THE EFCOCP INNER PROBLEM In this section, we provide more intuition about the role of z on the learned policy π and the learned value function V π . We first restate the EFCOCP inner problem below. ˜Jπ (x0, z) := max ( max k≥0 h(xk), ∞X k=0 l(xk) − z ) . (B.1) As z → −∞, the cost (i.e., stability) related term dominates the max. Consequently, we should see that the optimal policy will prioritize stability. On the other hand, as z → ∞, the constraint (i.e., safety) related term dominates the max. In this case, the optimal policy will prioritize safety. Moreover, if the optimal policy is safe under the unconstrained minimizer, i.e., hmax := max k≥0 h(xk) ≤ 0, (B.2) then the second term will be larger than the first. Since the second term is non-negative, we have that for any z ∈ (−hmax, 0], ∞X k=0 l(xk) − z ≥ −z > hmax (B.3) Consequently, for such a choice of z, ˜Jπ (x0, z), = max ( max k≥0 h(xk), ∞X k=0 l(xk) − z ) , (B.4) = max ( hmax, ∞X k=0 l(xk) − z ) , (B.5) = ∞X k=0 l(xk), (B.6) and we recover the unconstrained optimizer. We now compare the policy rollouts for different values of z on different systems. The policy rollouts for the 1D double- integrator system are shown in Figure 9. Note that for z = 0, states whose unconstrained minimizer are safe follow the unconstrained optimal trajectory and converge to the goal region. As z increases, the policy focuses more on constraint satisfaction. Consequently, trajectories that were originally unsafe (red) become safe (blue). However, as z increases further, the policy focuses too much on minimizing constraint function (i.e., maxk≥0 h(xk)) and converges to the minimizer of h instead of the goal region (olive).z=0.0 z=0.5 z=1.0 z=1.25 z=1.375 z=1.5 z=2.0 z=3.0 stable, safe stable, unsafe unstable, safe Fig. 9: Comparison of the policy π(·, z) for different values of z on the 1D double integrator. (a) z = 0.0 (b) z = 0.75 (c) z = 2.8 Fig. 10: Comparison of the policy π(·, z) for z = 0, z = 0.75 and z = 2.8 respectively from left to right on the Hopper system. The direction of time follows the colors red, purple, blue, green, yellow. Each row represents a different initial condition. The policy for z = 0.0 is too aggressive and eventually topples, violating the safety constraints. In contrast, the policy for z = 2.8 prioritizes safety by keeping the torso vertical, but hops very slowly in doing so and also overshoots the goal region (compare with z = 2.8). Taking z to be a value between these two extremes (e.g., z = 0.75) stabilizes to the goal while maintaining safety. By training a policy π that is conditioned on z, we can maintain safety and obtain a stabilizing controller despite π being suboptimal by learning a proper value of z∗. We next show the policy rollouts on the Hopper system for different values of z in Figure 10. Again, we can see that larger values of z correlates to higher emphasis on safety. In the case of Hopper, note that the optimal unconstrained optimizer should be able to maintain safety. However, despite the learned policy being suboptimal (and hence unsafe), we are still able to obtain a safe final policy by using z > 0. APPENDIX C DISCOUNTING IN EFPPO As noted in the main text, we randomly sample z from [z, zmax] when solving the inner problem of EFCOCP, where zmax is an upper bound of the total cost under the optimal policy π∗, i.e., zmax ≥ ∞X k=0 l(xk), (C.1) for any trajectory {xk}∞ k=0. However, if the system under the optimal policy does not stabilize to the zero-set of l fast enough (or not at all) due to lack of controllability, then this may be infinite. While this is not a problem for the solution of the optimization problem if the system is not controllable from x0, it is problematic when we apply reinforcement learning to the problem and learn a neural network that approximates the policy value function V π . Such a term will dominate the loss function when training V π . To alleviate this, we apply a small discount factor γ ∈ (0, 1), taken to be 0.97 in all of our experiments. Consequently, we now consider the discounted EFCOCP inner problem, where the cost function ˜Jπ now takes the form ˜J(x0, z) := max n max k≥0 γkh(xk), ∞X k=0 γkl(xk) − z o . (C.2) The dynamic programming equations are modified correspondingly, which we derive below. ˜V (x0, z0) = min u0:∞ max n max k≥0 γkh(xk), ∞X k=0 γkl(xk) − z o , (C.3) = min u0:∞ max n h(x0), max k≥1 γkh(xk), ∞X k=1 γkl(xk) − z − l(x0)
o , (C.4) = min u0 max n h(x0), min u1:∞ max  max k≥1 γkh(xk), ∞X k=1 γkl(xk) − z − l(x0)
o , (C.5) = min u0:∞ max ( h(x0), γ min u1:∞ max max k≥0 γkh(xk+1), ∞X k=0 γkl(xk+1) − z − l(x0) γ !) , (C.6) = min u0:∞ max  h(x0), γV  x1, z − l(x0) γ  , (C.7) = min u0:∞ max {h(x0), γV (x1, z1)} , (C.8) where the “dynnamics” for z now read zk+1 = zk − l(xk) γ . (C.9) Following this, the policy value function V π and policy action-value function Qπ used for EFPPO are modified accordingly. With the discounted formulation, V π is now finite assuming h(xk) does not explode and l(xk) does not grow faster than γk, which is satisfied in most practical problems where the system has enough control authority. While we can find zmax analytically, in practice zmax is found empirically by running the inner loop of EFPPO for several iterations and then taking zmax to be a constant multiple (e.g., 1.5) of the largest value of P∞ k=0 γkl(xk) seen so far. Since the initial policy is generally worse than π∗ (i.e., has larger cost), this procedure yields a conservative over-estimate of zmax that we have found to be robust. In our experiments, this procedure only needs to be performed once for every new task to set zmax and does not require any tuning afterwards. APPENDIX D SIMULATION DETAILS Details for the simulation environments used are provided below. D1 1D Double-Integrator States Controls Index Symbol Description Index Symbol Description 0 p Position 0 a Acceleration 1 v Velocity TABLE V: States and Controls for the 1D Double-Integrator The 1D Double-Integrator is a system with 2 state and 1 control dimensions (see Table V). The dynamics are linear and take the form pk+1 vk+1  = 1 ∆t 0 1 pk vk  +  1 2 ∆t2 ∆t ak  (D.1) for timestep ∆t. We use ∆t = 0.025. The control constraints are box constraints within [−1, 1] |a| ≤ 1. (D.2) The state constraints (which define the avoid set A) are |p| ≤ 1, |v| ≤ 1 (D.3) To represent A, we define h(x) = max h1(x), h2(x), where h1(x) := |p| − 1, h2(x) := |v|3 − 1. (D.4) The goal set G is defined as the region G := { x | p ∈ [0.65, 0.85] } , (D.5) which we represent via the cost function l as l(x) := |p − 0.75| − 0.1+. (D.6) D2 2D Single-Integrator with Sector Obstacle States Controls Index Symbol Description Index Symbol Description 0 px Position along x 0 vn Velocity along the normal to the origin 1 py Position along y 1 vt Velocity along the tangent to the origin TABLE VI: States and Controls for the 2D Single-Integrator The 2D Single-Integrator is a system with 2 state and 2 control dimensions (see Table VI). The continuous-time dynamics are as d dt px py  = 1 q p2 x + p2 y −px −py py −px vn vt  , (D.7) where the denominator is clipped to prevent division by 0. The discrete-time dynamics are obtained by discretizing the above using Euler integration with timestep ∆t = 0.05. The control constraints are box constraints within [−1, 1]2 |vn| ≤ 1, |vt| ≤ 1. (D.8) The state constraints are represented as the set h(x) = max(h0(x), h1(x)) ≤ 0, where h0(x) := r − 1, h1(x) := 0.2(1 − √2) + √2px − r. (D.9) where r := q p2 x + p2 y denotes the distance to the origin. h0 defines a circle with radius 1, while h1 defines the sector obstacle. The goal set G is defined as a circle at the origin with radius R = 0.05, which we represent via the cost function l(x) := [r − 0.05]+. (D.10) D3 Hopper Stabilization States Controls Index Symbol Description Index Symbol Description 0 px x-coordinate of the torso 0 τt Torque applied to the thigh motor 1 pz z-coordinate of the torso 1 τl Torque applied to the leg motor 2 θ Angle of the torso 2 τf Torque applied to the foot motor 3 θt Joint Angle of the thigh 4 θl Joint Angle of the leg 5 θf Joint Angle of the foot 6 vx Velocity of x-coordinate of the torso 7 vz Velocity of z-coordinate of the torso 8 ω Angular velocity of the torso 9 ωt Joint Velocity of the thigh 10 ωl Joint Velocity of the leg 11 ωf Joint Velocity of the foot TABLE VII: States and Controls for the Hopper The Hopper is a system with 12 state and 3 control dimensions (see Table VII) implemented using the Brax [15] simulator. Note that the version of Hopper we use includes px. The original Hopper environment from Brax and other simulators such as Mujoco [49] exclude px, since the goal is to learn a limit-cycle that is independent of px. In Brax, the dynamics are defined using the rigid body equations, which are then discretized using the default integrator settings (Euler integration, ∆t = 0.008). The control constraints are box constraints within [−1, 1]3 taken from the default settings. The state constraints are represented as the set h(x) = max(h0(x), h1(x)) ≤ 0, where h0(x) := 0.7 − pz , h1(x) := |θ| − 0.2, (D.11) which represent maintaining a minimum height and preventing the torso from tipping over too much, and are taken from the default settings. The goal set G is defined as the set of states where px is within the set [2.8, 3.0]. This is represented via the cost function l(x) := [|px − 2.9| − 0.1]+. (D.12) D4 F16 Ground Collision Avoidance in a Low Altitude Flight Corridor States Controls Index Symbol Description Index Symbol Description 0 vT Air speed 0 N zd Setpoint for accleration 1 α Angle of attack 1 P sd Setpoint for stability roll rate 2 β Angle of sideslip 2 N ypRd Setpoint for side acceleration and yaw rate 3 ϕ Roll 3 δt Throttle 4 θ Pitch 5 ψ Yaw 6 P Roll rate 7 Q Pitch rate 8 R Yaw rate 9 pN Northward displacement 10 pE Eastward displacement 11 pU Altitude 12 pow Engine power lag 13 N z Upward acceleration 14 P s Stability roll rate 15 N ypR Side acceleration and yaw rate 15 V Valid mask TABLE VIII: States and Controls for the F16 The F16 is a system with 17 state and 4 control dimensions (see Table VIII) based on [20]. Note that the original system includes only 16 states. We have added the final state V to prevent the state from exiting the region where the F16 model is numerically accurate. The continuous-time dynamics are defined in [20] is a standard model used in aerospace engineering and described extensively in the textbook by Stevens and Lewis [44]. Notably, the dynamics makes use of look-up tables for aspects such as the engine model and aerodynamic coefficients. The discrete-time dynamics are defined by integrating the continuous-time dynamics using the RK4 integrator with a step size of ∆t = 0.05. Moreover, when the system exits the region where the model is valid, defined as the set n x α ∈ [α, α], β ∈ [β, β], θ ∈ [− π 2 , π 2 ] o , (D.13) we set V to 0 and stop integrating the dynamics to prevent the model from misbehaving. The bounds for α and β are taken as the limits of the aerodynamic data tables, while the bounds for the pitch θ are used to avoid the singularity due to the use of Euler angles. The control constraints are box constraints defined as N zd ∈ [−10, 15], P sd ∈ [−10, 10], N ypRd ∈ [−10, 10], δt ∈ [0, 1] (D.14) The state constraints are represented as the set h(x) = maxi hi(x) ≤ 0, where the constraint functions hi are defined as (Avoid ground and stay below ceiling) h0(x) := max(0 − pU , pU − 1000) / 200, (D.15a) (Keep α valid) h1(x) := max(α − α, α − α) / 0.2, (D.15b) (Keep β valid) h2(x) := max(β − β, β − β) / 0.2, (D.15c) (Keep θ valid) h3(x) := max(θ − θ, θ − θ) / 0.2, (D.15d) (Stay within the flight corridor) h4(x) := max(−200 − pE , pE − 200) / 50. (D.15e) The goal region G is defined as the set of states where the altitude pU is within the set [50, 150]. Note that this set is very close to the ground and thus the reason why we call this task “low-altitude flight corridor”. We implement this via the cost function l(x) := max 50 − pU , pU − 150
/ 250. (D.16)