Neural Lyapunov Control for Discrete-Time Systems
Authors: Junlin Wu
Andrew Clark
Yiannis Kantaros
Yevgeniy Vorobeychik
Venue: NeurIPS 2023 poster (2023)
Keywords: nonlinear systems
Lyapunov stability
neural Lyapunov control
TL;DR: We propose a novel approach to learn stabilizing controllers with provable guarantees for nonlinear discrete systems.

Abstract: While ensuring stability for linear systems is well understood, it remains a major challenge for nonlinear systems. A general approach in such cases is to compute a combination of a Lyapunov function and an associated control policy. However, finding Lyapunov functions for general nonlinear systems is a challenging task. To address this challenge, several methods have been proposed that represent Lyapunov functions using neural networks. However, such approaches either focus on continuous-time systems, or highly restricted classes of nonlinear dynamics. We propose the first approach for learning neural Lyapunov control in a broad class of discrete-time systems. Three key ingredients enable us to effectively learn provably stable control policies. The first is a novel mixed-integer linear programming approach for verifying the discrete-time Lyapunov stability conditions, leveraging the particular structure of these conditions. The second is a novel approach for computing verified sublevel sets. The third is a heuristic gradient-based method for quickly finding counterexamples to significantly speed up Lyapunov function learning. Our experiments on four standard benchmarks demonstrate that our approach significantly outperforms state-of-the-art baselines. For example, on the path tracking benchmark, we outperform recent neural Lyapunov control baselines by an order of magnitude in both running time and the size of the region of attraction, and on two of the four benchmarks (cartpole and PVTOL), ours is the first automated approach to return a provably stable controller. Our code is available at: https://github.com/jlwu002/nlc_discrete.

Neural Lyapunov Control for Discrete- Systems
Junlin Wu
Andrew Clark
Computer Science & Engineering Washington University in St. Louis
Electrical Systems Engineering
Washington University in St. Louis
St. Louis, MO 63130 junlin. ww@wustl. edu
St. Louis, MO 63130
andrewclark@wusta edu
Yiannis Kantaros
Yevgeniy Vorobeychik
Electrical & Systems Engineering
Computer Science & Engineering Washington University in St. Louis
Washington University in St. Louis
St. Louis, MO 63130 ioanni wwe.inikwwwstt.ccmmmm edu
St. Louis, MO 63130
yworobeychikkvustt.. edu
Abstract
While ensuring stability for linear systems is well understood, it remains major challenge for nonlinear systems. A general approach in such cases is to compute a combination of Lyapunov function and an associated control policy. However, finding Lyapunov functions for general nonlinear systems a challenging task. To address this challenge, several methods have been proposed that represent Lyapunov functions using neural networks. However, such approaches either focus on continuous- time systems, or highly restricted classes of nonlinear dynamics. We propose the first approach for learning neural Lyapunov control in broad class of discrete-time systems. Three key ingredients enable to effectively learn provably stable control policies. The first is novel mixed-integer linear programming approach for verifying the discrete-time Lyapunov stability conditions, leveraging the particular structure of these conditions. The second is novel approach for computing verified sublevel sets. The third a heuristic gradient-based method for quickly finding counterexamples to significantly speed up Lyapunov function learning. Our experiments on four standard benchmarks demonstrate that our approach significantly outperforms state-of- of-t-offthhe-rt baselines. For example, the path tracking benchmark, we outperform recent neural Lyapunov control basclines by order of magnitude in both running time and the size of the region of attraction, and on two of the four benchmarks (cartpole and PVTOL), ours is the first automated approach to return provably stable controller. Our code is available at: htt////iti..cmm///////////ccccl/
Introduction
Stability analysis for dynamical systems aims to show that the system state will return to equi- librium under small perturbations. Designing stable control in nonlinear systems commonly relies
constructing Lyapunov functions that certify the stability of equilibrium points and estimate
their region of attraction. However, finding Lyapunov functions for arbitrary nonlinear systems is challenging task that requires substantial expertise and manual effort KKaliil 2015
Lavaei and
Bridgeman 2023). To address this challenge, recent progress has been made in learning Lyapunov functions continuous-tim nonlinear dynamics [Abate al. 2020 Chang al. 2019
Zhou et al.
2022]. However, few approaches exist for doing so in discrete-time systems [Dai et al.
2021]. and
none for general Lipschitz-continuous dynamics. Since modern learning-based controllers often
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

take non-negligible time for computation, the granularity of control is effectively discrete-time, and developing approaches for stabilizing such controllers is major open challenge.
We propose a novel method to learn Lyapunov functions and stabilizing controllers, represented as neural networks (NNs) with ReL.U activation functions, for discrete-time nonlinear systems. The proposed framework broadly consists of a learner, which uses a gradient-based method for updating the parameters of the Lyapunov function and policy, and verifier, which produces counterexamples (if any exist) to the stability conditions that are added as training data for the learner. The use of a verifier in the learning loop is critical enable the proposed approach to return a provably stable policy. However, no prior approaches enable sound verification of neural Lyapunov stability for general discrete-tine time dynamical systems. The closest is Dai al. (2021]. who assume that dynamies are represented by a neural network, an assumption that rarely yy holds in real systems. On the other hand, approaches for continuous-timg systems [Chang et al. 2019 Zhou et al. [022] have limited efficacy in discrete-time environments, as our experiments demonstrate. 'o address this gap, we develop novel verification tool that checks if the candidate NN Lyapunov function satisfies the Lyapunov conditions by solving Mixed Interger Linear Programs (MILPs). This approach, which takes advantage of the structure of discrete-time Lyapunov stability conditions, can soundly verify a broad class of dynamical system models.
However, using sound verification tool in the learning loop makes it significant bottleneck, severely limiting scalability. We address this problem by also developing a highly effective gradient- based technique for identifying counterexamples, resorting to the full MILP-based verifier only as last resort. The full learning process thereby iterates between learning and verification steps, and returns only when the sound verifier is able to prove that the Lyapunov function and policy satisfy the stability conditions.
The final technical challenge stems from the difficulty of verifying stability near the origin [Chang iet al. 2019]. typically addressed heuristically by either adding fixed tolerance to a stability condi- Tion Dai et al. 202T]. J. or excluding small area around the origin from verification [Chang et al. 2019 Zhou et al. 2022]. We address it by adapting Lyapunov stability theory to ensure convergence to small region near the origin, thereby achieving the first (to our knowledge) sound approach for computing a stable neural controller that explicitly accounts for such approximations near the origin.
We evaluate the proposed approach in comparison to state-of-the-art baselines on four standard nonlinear control benchmarks. On the two simpler domains (inverted pendulum and path following), our approach outperforms the state-of-the-art continuous-time neural Lyapunov control approaches by at least several factors and up to an order of magnitude in both running time and the size of the region of attraction. On the two more complex cartine and PVTOL-ours the first automated approach that returns a provably slable controller. Moreover, our ablation experiments demonstrate that both the MILP-based verifier and heuristic counterexample generation technique we propose are critical to the success of our approach.
In summary, make the following contributions:
novel MILPP-baeed approach for verifying a broad class of discrete-time controlled nonlinear systems.
A novel approach for learning provably verified stabilizing controllers for a broad class of discrete- time nonlinear systems which combines our MILP-based verifier with heuristic gradient-based approximate counterexample generation technique.
A novel formalization of approximate stability in which the controller provably converges to small ball near the origin in finite time.
Extensive experiments
four standard benchmarks demonstrate that by leveraging the special structure of Lyapunov stability conditions for discrete- tecrette-timeeeeee system, our approach significantly outperforms prior art.
Related Work Much prior work on learning Lyapunov functions focused on continuous- time
systems [Abate et al.] [220 Ravanbakhsh and Sankaranarayanan, 2019 Chang et al., 2019| Zhou
et al. 2022]. Common approaches have assumed that dynamics are linear Donti et al. [021 Tedrake 2009] or polynomial (Ravanbakhsh and Sankaranarayanan, 2019). Recently, Chang et al. (001))
Rego and de Araújo| (2022). Zhou et al. [202]] proposed learning Lyapunov functions represented as neural networks while restricting policies to be linear. These were designed for continuous-time
dynamics, and are not effective in discrete-time settings, we show in the experiments. Chen et al.
2

---

[02 learns convex Lyapunov functions for discrete-timg hybrid systems. Their approach requires hybrid systems admit mixed-integed linear programming formulation, essentially restricting it to piecewise-affine systems, and does not learn stabilizing controllers for these. Our approach considers a much broader class of dynamical system models, and learns provably stable controllers and Lyapunov functions, allowing both to be represented as neural networks. Kolter and Manek [2019) learn stable nonlinear dynamics represented as neural networks, but do not provide stability guarantees with respect to the true underlying system. In addition, several approaches have been proposed that either provide only probabílistic stability guarantees [Berkenkamp et al.) 2017 Richards et al. 2018), or
do not guarantee stability (Choi et al.) 2020 Han et al.. 2020). Several recent approaches propose methods for verifying stability in dynamical discrete- systems. Chen et al. (2021b] compute an approximate region of attraction of dynamical systems with neural network dynamics, but assume that Lyapunov functions are quadratic. |Da et al. [2021] consider the problem of verifying Lyapunov conditions in discrete-time systems, well as learning provably stable policies. However, they
assume that dynamics are represented by neural networks with (leaky) ReLU activation functions.
Our verification approach, in contrast, is for arbitrary Lipschitz continuous dynamics.
2 Model
We consider a discrete-time nonlinear dynamical system
(I) where I, E X' is state in domain X' and 2, the control input at time t. We assume that f is Lipschitz continuous with Lipschitz constant L. This class of dynamical systems includes the vast majority of (non-hybrid) dynamical system models in prior literature. We assume that = 0 is an equilibrium point for the system, that is, f(0, uo) = 0 for some ut. Let (Tx) denote a control policy, with u, T(x)). For example, in autonomous vehicle path tracking, can measure path tracking error and u the steering angle. In the conventional setup, the goal is to learn a control policy such that the dynamical system in Equation ① converges to the equilibrium x 0, condition referred to as stability. To this end, we can leverage the Lyapunov stability framework (Tedrake] identify a Lyapunov function V (z) and controller π that satisfies the
T++1 [(I1, we),
2009]. Specifically, the goal Following conditions over a subset of the domain R S X:: (0) = O; 2)V(() > O, ve E R {0}; R. These conditions imply that the system (locally)
and 3)V (f(z, (x))) V(x) < O, Vr
stable in the Lyapunov sense [Bof et al. 2018 Tedrake (2009].
In practice, due to numerical challenges in verifying stability conditions near the origin, it is common to verify slight relaxations of the Lyapunov conditions. These typically fall into two categories: I) Lyapunov conditions are only checked for ||>|| 2 (Chang al. 2019 Zhou al. 2022] (our main baselines, and the only other methods for learning neural Lyapunov control for general nonlinear dynamics), and/or 2) small tolerance
is added in verification, allowing small violations of
Lyapunov conditions near the origin Dai et al. 2021].
Our goal now is to formally investigate the implications of numerical approximations of this kind. We next define a form of stability which entails finite-time convergence to B(0, c) = {x || < c). Definition 2.1 (e-stability). We call pair of (V, x) e-stable within region R when the following conditions are satisfied: (a) V(0)= = 0: (b) there exists Ç >> such that ((f(x, (x)) V(x) < - for all x E R B(0, <); and (c) V(x)= >> 0 for all x E RR (()
A key ingredient to achieving stability is to identify region of altraction (ROA), within which we are guaranteed to converge to the origin from any starting point. In the context of e-stability, our goal is to converge to a small ball near the origin, rather than the origin; let e-ROA be the set of initial inputs that has this property. In order to enable us to identify an e- E-OA,, we introduce a notion of an invariant sublevel set. Specifically, we refer to a set D(R, p) {a E R (x) S p} with the property that x E D(R. p) + fI, (()) E R as a R-invariant sublevel set. In other words, D(R, p) is a sublevel set which is additionally forward invariant with respect to R. We assume that (0, e) C D(R,
Next, we formally prove that c-stability combined with a sublevel set D(R. p) entails convergence in three senses: 1) that we reach an c- e-ball around the origin in finite time, 2) that we reach it infinitely often, and 3) we converge to an arbitrarily small ball around the origin. We refer to the first of these senses as reachability. What we show is that for sufficiently small, reachability implies convergence in all three senses. For this, a key additional assumption is that V and n are Lipschitz continuous, with Lipschitz constants L, and L,, respectively (c.g., in ReLU neural networks).
3

---

Theorem 2.2. Suppose V and πτ are e-stable on a compact R. D(R, p) is an R-invariant set, and Ee such that ~ollo I C,E. Then 2f E D(R. p) ((0,e) e):
(i) there exists a finite K such that J.K E B(( c)
(ii) Boz such that if e2e < p and ((0, c26) C R, then there exists a finite K such that Vk 2 K, Ip E D(R,, o26) and the sequence { Je20 reaches Bi0, e) infinitely often, and furthermore (iii) for any 1l > such that ((0, 11) C R, Ee and finite K such that [[=(0] uoll do S C1E +
| S vk 2 K.
Proof. We prove (i) by contradiction. Suppose that || Vk E (0. [V //1). Then, when & = [V( condition (b) of e-stability and R-invariance of DRR, p) implies that (Ix) < (zo) Çk- < 0, contradicting condition (c) of e-stability.
To prove (i), fix the finite K from (i), and let k 2 K be such that z E e). By Lipschitz continuity of V and the fact that V(0) = 0 and V(x) 2 0 (conditions (a) and (c) of c-stability), ((I++1) < Ly + where || || is the fe norm here and below. Moreover, since #x+1 = f(z,, T (Tx)), f0, uo) = 0 (stability of the origin), and by Lipschitz continuity of f and x,
| || S Ly (1=)) - uo)ll S Lg max { [| ll, [|((=)) woll
By Lipschitz continuity of rT, triangle inequality, and the condition that []) woll S cıc, we have [ - Holl S - [| c S (Lx + (2)c. Let 02 = max{ 1 max(1, L + t.). Then || S eae and (II+1I) I Czc. Thus, if e2E < p and (00,, czé) C RR and D(R. p) is R-invariant, then I8+1 D(R, which R-invariant. Additionally, either 2x+1 B(O.c). and by the argument above F&+2 E D(R, o2e), or FE+1 E D(R, C24) \ BB0,, e). and by e-stability Tk+2 E D(R, c2é). Thus, by induction, we have that for all 2 K, Fy E D(R, o2c). Finally, since for all & 2 K, Tg E D(R, c. if Ih & (B,, c). by (i) it must reach B(0, <) in finite time. Consequently, there is an infinite subsequence {aw) of {aa3=20 such that all TM E B(0, e). that is, the sequence {# reaches ((0, e) infinitely often.
We prove part (iii) by contradiction. Fix 7) > and define S = fo R | > n). Suppose that Ve > 0 there exists x E S such that V(z) S 62e where e2 is as in (ii). Then for any (infinite) sequence of {er} we have {x,} such that V(xx) S C2C,, where I₁ E S. Now, consider a set S= fx E R. : ||>||| 2 m)- Since S is compact and {2, } is an infinite sequence, there is an infinite subsequence (wt, } of {1x} such that limg St, = x² and X² E S. Since V continuous, we
have limg- (T = (x=). Now, we choose {{t] such that lim = 0. This means that limg->oo (( V(x*) S O, and since x" E S, this contradicts condition (c) of e-stability. Since
by (ii) there exists a finite K, V (xxx)= S C2e, Vk 2 K, we have vk 2 K, | 17.
The crucial implication of Theorem 2.2 is that as long as we choose > 0 sufficiently small, verifying that V and π are e-stable together with identifying an R-invariant set D(R. p) suffices for convergence arbitrarily close to the origin. One caveat is that we need to ensure that π (0) is sufficiently close to uo for any
choose. In most domains, this can be easily achieved: for example, if uo 0 (as is the case in many settings, including three of our four experimental domains), we can use a representation for ı with no bias terms, so that (0) = 0 = un by construction. In other cases, we can simply check this condition after learning.
Another caveat is that we need to define R to enable us to practically achieve these properties. To this end, we define R parametrically as R(7) = {a E X
|| S <>). Additionally, we introduce the < ||0||>0 S 7) Note that R(0, x)) = X, R(e, oo) = {x E X ||=]||>o 2 } and R(0, 7) = R(7). Thus, for sufficiently large compared to e, conditions such as B()0,) C D(R, p). B(0, r) C R. and B(0, c2t) C R. will be easy to satisfy. Following conventional naming, we denote R(e, 7) as e-valid region, i.e., the region that satisfies conditions (a) (c) of e-stability, and refer to a function V satisfying these conditions as an e-Lyapunov function.
following useful notation. Define R(e, 7) {x E X
We assume that the c-Lyapunov function as well as the policy can be represented in a parametric function class, such as by deep neural network. Formally, we denote the parametric c-Lyapunov function by Vz(z) and the policy by [2),,, where 0 and 3 are their respective parameters. Let D(1, p) = D(R( p). Our goal is to learn Vo and TH such we maximize the size of the R(+)- invariant sublevel set D(7;p) p) such that Vo and TT are provably e-stable on R(7). Armed with Theorem 2.2 we refer to this set simply as ROA below for simplicity and for consistency with prior
work, e.g., Chang et al. 2019 Zhou et al. 2022). noting all the caveats discussed above.
4

---

Define P(n) as the set (6, 5) for which the conditions (a)-(c) of e-stability are satisfied over the domain R(y). Our main challenge below is to find (0.3) E P(y) for a given y. That, in turn, entails solving the key subproblem of verifying these conditions for given Vo and mgg-
Next, in Section 3 we address the verification problem, and in Section f we describe approach for jointly learning Ve and Ty that can be verified to satisfy the e-stability conditions for given R(7).
Verifying Stability Conditions
Prior work on learning <-I-yapunvv functions for continuous-time nonlinear control problems has leveraged off-the-shelf SMT solvers, such as dReal [Gao et al.) 2013). However, these solvers scale poorly in our setting (see Supplement[[foor for details). In this section, we propose novel approach for verifying the e-Lyapunov conditions for arbitrary Lipschitz continuous dynamics using mixed-integer linear programming, through obtaining piecewise-linear bounds on the dynamics. We assume Vo and Ty are K- and N-layer neural networks, respectively, with ReLU activations.
We begin with the problem of verifying condition (c) of e-stability, which represent feasibility problem: to find if there is any point ă E R(E, 7) such that V I 0. We can formulate it as the following MILP:
=K 1 0
(2a)
Z1+1 = 9o, (21);
0 0SISK- 1
(2b)
S | S 7
20 = I,
(2c)
where refers to a layer in the neural network Vo. 2k Vo(x),, and the associated functions ge, (2))
are either W/2 + by for a linear layer (with 8: (W, by)) or max 2,, 0} for ReLU activation. Any feasible solution x² is then a counterexample, and if the problem is infeasible, the condition is
satisfied. ReLU activations g(=) can be linearized by introducing an integer variable E {0, 1), and replacing the z' g(2) terms with constraints S Ua and E' -1((1 a), where and U are specified so that T. < U (we deal with identifying T. and U below).
Next, we cast verification of condition (b) of e-stability, which involves the nonlinear control dynamics f, as the following feasibility problem:
Σ 25 2 -
(3a)
Vi+1 h1, (W:),
0 0sigN-1
(3b)
žt+1 -90, (ž)),
0 0SIsK-1
(3c)
5g = fax yN)
(3d)
yo constraints
2b
2cp
(3e)
where ha, () are functions computed at layers i of nT. zK-
(z(z), and žK Va(f T (z))).
At this point, all of the constraints can be linearized as before with the exception of Constraint (3d).
which involves the nonlinear dynamics f. To deal with this, we relax the verification problem by replacing f with linear lower and upper bounds. To obtain tight bounds, we divide R(y) into collection of subdomains {R For each R,, we obtain linear lower bound fix(2) and upper bound (x) on f, and relax the problematic Constraint (3d) into fIa S 2o S (x)) which
now pair of linear constraints. We can then solve Problem for each RR-
Computing Linear Bounds System Dynamics Recall that f R + R. is the Lipschitz-- continuous system dynamic x++1 = f(z,, u.). For simplicity we omit = Ja (z) below. Let A be the Loo Lipschitz constant f. Suppose that we are given region of the domain R2 represented
as a hyperrectangle, i.c., R₂ = x, [ Ti,u,, where [ Ti,al are the lower and upper bounds of x along coordinate i. Our goal is to compute a linear upper and lower bound on f(x) over this region. We bound fy (z) R + R along each j-th dimension separately. By X-Lipschiz-coontinuty we Ex -x--- xx,1) and (x) 2 f 22,,,.,,,,nI..
can obtain (z) S (11, Fy.............................. In.1)
(x; J. The full derivation is provided in the Supplement Section B
Alternatively, if fj is differentiable, convex Ii, and monotone over other dimensions, can restrict it to calculating one-dimensional bounds by finding coordinates - and x_1, such that (11, < (1,1,, <
for any given æ. Then fz(z) I f)(,u(())
fy (Z -.) + I. f,)-f,(a,,1e snd) (x+ IT.1) and fz( 2 fr.((1) = (2)))=111) +
20, 1.
5

---

( 2;), where we let 2 = ming, (( Ienn Se. (x1: - fg (Zx, x-1,) + f2 I.)) ddr,, to minimize the error area between fj (z) and f1,0(()). Note that the solution for z can be approximated when conducting experiments, and this approximation has no impact on the correctness of the bound. A similar result obtains for where fz(#) concave over zi. Furthermore, we can use these single-dimensional bounds to obtain tighter lower and upper bounds on [z) as follows: note that for any c²c² with Σ. = and Σ.e² = 1, E; (f) S fx( S E; & (z), which means we can optimize c' and c* to minimize the bounding error. In practice, we can typically partition R(1) so that the stronger monotonicity and convexity or concavity assumptions hold for each Rs-
Note that our linear bounds f1)(z) and fz)(2) introduce errors compared to the original nonlinear dynamics f. However, we can obtain tighter bounds by splitting R2 further into subregions, and computing tighter bounds in each of the resulting subregions, but at the cost of increased computation. To balance these considerations, we start with a relatively small collection of subdomains (R&), and only split a region R. if we obtain counterexample in Rg that is not an actual counterexample for the true dynamics f.
Computing Bounds on ReLU Linearization Constants In linearizing the ReLU activations, we supposed an existence of lower and upper bounds L and U. However, we cannot simply set them to some large negative and positive number, respectively, because (z(z) has no a priori fixed bounds (in particular, note that for any e-Lyapunov function V and constant a > 1, aV is also a e- e-Lyapunov function). Thus, arbitrarily setting L and U makes verification unsound. To address this issue, we use interval bound propogation (IBP) [Gowal al., 2018) to obtain M = max ( [Us]. [,]]: where U, is the upper bound, and L; is the lower bound returned by IBP for the i-th layer, with inputs for the first layer the upper and lower bounds of (R(())). Setting each = M and U = M then yields sound verification.
Computing Sublevel Sets The approaches above verify the conditions (a)-(c) of e-stability on R((). The final piece is to find the R(())-invariant sublevel set D(7, p), that is, to find p. Let B(x) 2 max (max [I ma(x)) (oo, ()))))))))))))))) 7). We find p by solving
min V(x).
my || los BB))
We can transform both Problem (4) and computation of B(7) into MILP for other problems above.
Theorem 3.1. Suppose that V and are e-stable on R(7)) ll((0) uol co I 6, and γ 2 Ls max (1, Lx + 1Je. Let V be the optimal value of the objective in Problem 4 and p V μ for any > O. Then the set D(T,P) p)
= (x 10:00 E R(7), S p} is an ((r)-invariam sublevel set.
V(z) S pps V by e-stability of V and ग. Suppose that . = [sx m(x)) & R(7). Since (TT) < V V must not be an optimal solution to contradiction. |o S e, the argument is similar to Theorem 2.2 (i)).
Proof. If ||=||<c > e, v(f(x,
4 Learning e-Lyapunov Function and Policy
We cast learning the e-I-yapunov function and policy as the following problem:
min 0,
Cla;: Tp).
(5)
FSS
where C()) is a loss function that promotes Lyapunov learning and the set S C R(7) is finite subset of points in the valid region. We assume that the loss function is differentiable, and consequently training follows sequence of gradient update steps 8') = (0, 3) Eres Ve,aC(z.; Vo, TT).
Clearly, the choice of the set S is pivotal to successfully learning Vo and F. with provable stability properties. Prior approaches for learning Lyapunov functions represented by neural networks use one of two ways to generate S. The first is to generate fixed set comprised of uniformly random samples from the domain [Chang and Gao 2021 Richards et al. 2018). However, this fails to learn verifiable I Lyapunov functions. starting with randomly generated samples in the domain, and then using the verification tool, such as dReal, in each step of the learning process to update (Chang et al. 2019 Zhou et al. 2022]. However, verification is often slow, and becomes a significant bottleneck in the learning process.
the second, S is not fixed, but changes in each learning iteration,

---

Our key idea is to combine the benefits of these two ideas with fast gradient- based heuristic approach for generating counterexamples. In particular, our proposed approach for training Lyapunov control functions and associated control policies (Algorithm[]) DITL.) involves five parts:
1. heuristic counterexample generation (lines 10-12),
2. choosing a collection S of inputs to update and & in each training (gradient update) iteration (lines 9-13),
3. the design of the loss function L.
4. initialization of policy Th (line 3), and
5. warm starting the training (line 4).
We describe each of these next.
Algorithm 1 DITL Lyapunov learning algorithm.
1: Input: Dynamical system model f' u) and target valid region R(x)
2: Output: Lyapunov function Ve and control policy n
3: T Initialize()
4: Ve PreTrainLLyapunoo(MW n&)
5: B = InitializeBuffer()
6: W <- /
7: while True do
8: for N iterations do
9:
Sample(r, B) //sample points from B
T Sample(g. //smpll //sample q points from R(7)
11:
T PGD(7)
T" FillerCouteEExamples(("))
13:
SUT"
14:
B = B B
15:
() + (0, B) - ErEs Vo,ah(an kp)
16:
end for
17: (success, W) = Veriffy(V, ng)
18:
if success then
19:
return Vo, n
20: else
21:
W + W WW
22:
end if
23:
end while
return FAILED // Timeout
Heuristic Counterexample Generation As discussed in Section verification in general involves two optimization problems: min Va(x) (o ne (x))) and minzeR(-)) Va().. We propose a projected gradient descent (PGD) approach for approximating solutions to either of these problems. PGD proceeds as follows, using the second minimization problem as an illustration; the process is identical in the first case, Beginning with a starting point 20 E R(y), we iteratively update Ik+1;
Tk+1 [1 + ansgg VVl
where [T projects the gradient update step onto the domain R(7) by clipping each dimension of Tk. and og is the PGD learning rate. Note that as is typically differentiable, directly through it when we apply PGD for the first verification problem above.
take gradients
We run this procedure for K iterations, and return the result (which may may not be actual counterexample). Let T" =PG((T) denote running PGD for a set of T' starting points for both of the optimization problems above, resulting in corresponding set T" of counterexamples.
Selecting Inputs for Gradient Updates A natural idea for selecting samples is simply add counterexamples identified in each iteration of the Lyapunov learning process. However, as S then grows without bound, this progressively slows down the learning process. addition, is important for S to contain relatively broad sample of the valid region to ensure that counterexamples we generate along the way do not cause undesired distortion of V, and n. subregions not well represented in S. Finally, need to retain the memory of previously generated counterexamples, as
7

---

these are often particularly challenging parts of the input domain for learning. We therefore construct S as follows.
We create an input buffer B, and initialize it with a random sample of inputs from R(y). Let W denote a set of counterexamples that we wish to remember through the entire training process, initialized to be empty. In our implementation, W consists of all counterexamples generated by the sound verification tools described in Section 3
At the beginning of each training update iteration, we randomly sample a fixed-size subset Š of the buffer B. Next, take a collection T of random samples from R to initialize PGD, and generate T = PGD(T). We then filter out all the counterexamples from T", retaining only those with either (=(x) < 0 or (z) - Vo(f To (z)x)) S where é is non-negative hyperparameter; this yields set T We then use S = S U T U W for the gradient update in the current iteration. Finally, we update the buffer B by adding to it 7" This process is repeated for N iterations.
After N iterations, check to see if have successfully learned a stable policy and a Lyapunov control function by using the tools from Section 3 to verify Ve and πg. If verification succeeds, training is complete. Otherwise, use the counterexamples W generated by the MILP solver to update W W U W, and repeat the learning process above.
Loss Function Design Next, we design of the loss function L(x; Vo, TT) for given input x. A key ingrediant in this loss function is term that incentivizes the learning process to satisfy condition (b) of c-stability: L1 (z; Va, R = ReL.U(V,
{=()))) - V*(x) + 7)), where the parameter v) 2 determines how aggressively we try satisfy this condition during learning.
There are several options for learning Vo satisfying conditions (a) and (c) of e-stability, The simplest is to set all bias terms in the neural network to zero, which immediately satisfies Vo = 0. An effective way to deal with condition (c) is to maximize the lower bound on (x)).. To this end, we propose to make use of the following loss term: Le Vo, Rg) = Ve L.B.
"), where V is the lower bound on the Lyapunov function over the target domain R(7). We use use a, 3-CRÓWN [Xu et al. 2020 to obtain this lower bound.
The downside of setting bias terms zero is that we lose many learnable parameters, reducing flexibility of the neural network. If we consider a general neural network, on the other hand, it is no longer the case that Vi(x)) 0 by construction. However, it is straigthforward to ensure this by defining the final Lyapunov function as Vz(2) = Vo(x) - (0). Now, satisfying condition (c) amounts to satisfying the condition that (z(z) 2 Vo(O), which we accomplish via the following pair of loss terms: L. (2; Vo, πT) = (e) - Vz(x) + j min([]=][a, v)), where μ and are hyperparameters of the term u min( ( which effectively penalizes Vo(x) for being too large, and L4 (z; Vo, лT = (0) /.
The general loss function is then a weighted sum of these loss terms,
C(x; Ve, RT) L1 (a;; Ve, mp) + c2; Vo, xp) + caCs(æ; TB) cxC1(z; Vo, TB). When we set bias terms to zero, we would set = = O; otherwise, we set Co 0.
Initialization We consider two approaches for initializing the policy Ng. The first is to linearize the dynamics/ f around the origin, and use a policy computing based on a linear quadratic regulator (LQR) to obtain simple linear policy ng. The second approach is use deep reinforcement, such as PPO (Liu et al., 2019), where rewards correspond to stability (e.g., reward is the negative value of the l2 distance of state to origin). To initialize Vo, fix T3 after its initialization, and follow learning procedure above using solely heuristic counterexample generation to pre-train Vo-
The next result now follows by construction.
Theorem 4.1. If Algorithm returns Ve, Tg. they are guaranteed to satisfy e-stability conditions.
5 Experiments
5.1
Experiment Setup
Benchmark Domains Our evaluation of the proposed DITL approach uses four benchmark control domains: inverted pendulum, path tracking, cartpole, and drone planar vertical takeoff and landing (PVTOL). Details about these domains are provided in the Supplement.
8

---

Baselines We compare the proposed approach to four baselines: linear quadratic regulator (LQR), sum-of-squares (SOS), neural Lyapunov control (NLC) [Chang et al.) 2019), and neural Lyapunov control for unknown systems (UNI.) [Zhou et al. 2022] The first two, L.QR and SOS, are the more traditional approaches for computing Lyapunov functions when the dynamies are either linear (LQQR) or polynomial (SOS) [Tedrake] 2009). LQQ solutions (when the system can be stabilized) are obtained through matrix multiplication, while SOS entails solving a semidefinite program (we solve it using YALMIP with MOSEK solver in MATLAB 2022b). The next two baselines, NLC and UNL, are recent approaches for learning Lyapunov functions in continuous-time systems (no approach for learning provably stable control using neural network representations exists for discrete time systems). Both NLC and UNL yield provably stable control for general non-linear dynamical systems, and are thus the most competitive baselines to date. In addition to these baselines, consider ablation in which PGDDbased counterexample generation for our approach is entirely replaced by the sound MlLL--aased method during learning (we refer to it as DITL-MILLP)
Efficacy Metrics We compare approaches terms of three efficacy metries. The first is (serial) runtime (that is, if no parallelism is used), which we measure wall clock time when only single task is running on a machine. For inverted pendulum and path tracking, all comparisons were performed machine with AMD Ryzen 9 5900X 12-Core Processor and Linux Ubuntu 20.04.5 LTS oS. All cartpole and PVTOL experiments were run on a machine with a Xeon Gold 6150 CPU (64-bit 18-core x86), Rocky Linux 8.6. UNL and RL training for Path Tracking are the only two cases that make use of GPUs, and was run on NVIDIA GeForce RTX 3090. The second metric was the size of the valid region, measured using l norm for NLC and UNL, and fxx norm (which dominates (2) for LQR, SOS, and our approach. The third metric is the region of attraction (ROA). Whenever verification fails, set ROA to O. Finally, compare all methods whose results stochastic (NLC, UNL, and ours) in terms of success rate.
Verification Details For LQR, SOS, NLC, and UNL, we used dReal as the verification tool, as done by Chang et al. 2019 and Zhou et al. 2022 For DITL verification we used CPLEX version 22.1.0.
5.2 Results
160
10
140
120
Angular Velectt
100
Angle
-3-
40
12
-0
-2
-10-
0
0
-20
-3-
-10
10
Angle (rad)
Distance Error
Our Approach
NLC
UNL
505
LOR
Volid Region
Figure I: ROA plot of inverted pendulum (left) and path tracking (right). We select the best result for each method.
Inverted Pendulum For the inverted pendulum domain, we initialize the control policy using the LQR solution (see the Supplement for details). We train Vo with non-zero bias terms. We set e == 0.1 (<C 0.007% of the valid region), and approximate ROA using a grid of 2000 cells along each coordinate using the level set certified with MILP ④ All runtime capped at 600 seconds, except the UNL baseline, which we cap at 16 minutes, as it tends to be considerably slower than other methods. Our results are presented in Table While LQR and SOS are fastest, our approach (DITL) is the next fastest, taking average ~ seconds, with NLC and UNL considerably slower. However, DITL yields an ROA a factor >4 larger than the nearest baseline (LQR), and 100% success rate. Finally, the DITL-MILP ablation is two orders of magnitude slower (runtime 300 seconds) and less effective (average ROA=42, 80% success rate) than DITL. We visualize maximum ROA produced by all methods in Figure [ (left).
Path Tracking In path tracking, we initialize our approach using both the RL and LQR solutions, drawing direct comparison between the two (see the Supplement for details). We set 0.1. The running time of RI. is ~ 155 seconds. The results are provided in Table 2 We can observe that both RL- and LQR- initialized variants of our approach outperform all prior art, with RL exhibiting

---

Table I: Inverted Pendulum
Valid Region
Runtime (s)
ROA
Max ROA
Success Rate
NLC (free)
x 2 S 6.0
28 +
11 + 4.6
22
100%
NLC (max torque 6.0)
|x] 2 S 6.0
519 = 184
13 + 27
66
20%
UNL (max torque 6.0)
xx 2 S 4.0
821 227
1
30%
LQR
[x ac S 5.8
<1
14
14
success
SOS
|x 1.7
<1
6
6
success
DITL
| I 12
8.1 + 4.7
61
+ 31
123
100%
Table 2: Path Tracking
Valid Region
Runtime (s)
ROA
Max ROA
Success Rate
NLC
2 S
1.0
109 + 81
0.5 + 0.2
0.76
100%
NLC
z 2 S 1.5
151 238
1.4 149
2.8
80%
UNL
2 2 I 0.8
925 + 110
0.1
0.56
10%
LQR
z S 0.7
1.02
1.02
success
SOS
>o S 0.8
<
1.8
1.8
success
DITL (LQR)
X0 S 3.0
9.8 4
8
12.5
100%
DITL (RL)
0 S 3.0
14 11
9 3.5
16
100%
factor of 5 advantage over the next best (SOS, in this case) in terms of ROA, and nearly factor of advantage in terms of maximum achieved ROA (NLC is the next best in this case). Moreover, our approach again has a 100% success rate. Our runtime is order of magnitude lower than NLC or UNL. Overall, the RL-initialized variant slightly outperforms LQR initialization. The DITL-MILP ablation again performs far worse than DITL: running time several orders of magnitude slower (at >> 550 seconds), with low efficacy (ROA is 1.1, success rate 10%). We visualize comparison of maximum ROA produced by all methods Figure (right).
Cartpole For cartpole, we used LQR for initialization, and set bias terms of Vo to zero. We set 0.1 (0.01% of the valid region area) and the running time limit hours for all approaches. None of the baselines successfully attained provably stable control policy and associated Lyapunov function for this problem. Failure was either because could find counterexamples within the target valid region, or because verification exceeded the time limit. DITL found valid region of ||~||> S 1.0, < 1.6 hours with a 100% success rate, and average ROA of 0.021 + 0.012.
PVTOL The PVTOL setup was similar to cartpole. We set = 0.1 (0.0001 % of the valid region area), and maximum running time to 24 hours. None of the baselines successfully identified provably stable control policy. In contrast, DITL found one within 13 6 hours on average, yielding 100% success rate. We identified valid region of ||a]|| oc < 1.0, and ROA of 0.011 + 0.008.
Conclusion
We presented novel algorithmic framework for learning Lyapunov control functions and policies for discrete-time nonlinear dynamical systems. Our approach combines miixed-nntege linear pro- gramming verification tool with training procedure that leverages gradient-based based approximate verification. This combination enables a significant improvement in effectiveness compared to prior art: experiments demonstrate that approach yields several factors larger regions of attraction in inverted pendulum and path tracking, and ours is the only approach that successfully finds stable policies in cartpole and PVTOL.
Acknowledgments
This research was partially supported by the NSF (grants CNS-1941670, ECCS-2020289, IIS- 1905558, I1S-2214141, and CNS-2231257), AFOSR (grant 119955022--100555 ARO (grant W911NF-I1-110041D, and NVIDIA.
10

---

References
Alessandro Abate, Daniele Ahmed, Mirco Giacobbe, and Andrea Peruffo. Formal synthesis of 1yapunov neural networks. IEEE Control Systems Letters, 5(3):773-778 2020.
Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe model-based reinforcement learning with stability guarantees. Advances neural information processing systems, 30, 2017.
Nicoletta Bof, Ruggero Carli, and Luca Schenato. Lyapunov theory for discrete time systems. arXiv preprint
2018.
Ya-Chien Chang and Sicun Gao. Stabilizing neural control using self-learned almost Iyapunov critics. In IEEE International Conference on Robotics and Automation, pages 1803-1809, 2021.
Ya-Chien Chang, Nima Roohi, and Sicun Gao. Neural Iyapunov control. In Neural Information Processing Systems, 2019.
Shaoru Chen, Mahyar Fazlyab, Manfred Morari, George J Pappas, and Victor M Preciado. Learning Yуарunov functions for hybrid systems. In Proceedings of the 24th International Conference on Hyybrd Systems: Computation and Control, pages 1-11, 2021a.
Shaoru Chen, Mahyar Fazlyab, Manfred Morari, George J Pappas, and Victor M Preciado. Learning region of attraction for nonlinear systems. In 2021 60th IEEE Conference on Decision and Control
(CDC), pages 6477-6484. IEEE, 2021b.
Jason Choi, Fernando Castaneda, Claire Tomlin, and Koushil Sreenath. Reinforcement learning for safety-critical control under model uncertainty, using control Iyapunov functions and control barrier functions. arXiv preprint wrXiv:2004.07584, 2020.
Hongkai Dai, Benoit Landry. Lujie Yang, Marco Pavone, and Russ Tedrake. Lyapunov-stable neural-network control. In Robotics: Science and Systems (RSS), 2021.
Priya L Donti, Melrose Roderick, Mahyar Fazlyab, and Zico Kolte. Enforcing robust control guar- antees within neural network policies. In International Conference Learming Representations, 2021.
Sicun Gao, Soonho Kong, and Edmund M Clarke. dreal: An smt solver for nonlinear theories over the reals. In International conference automated deduction, pages 208-214. Springer, 2013.
Sven Gowal, Krishnamurthy Dvijotham, Robert Stanforth, Rudy Bunel, Chongli Qin, Jonathan Uesato, Relja Arandjelovic, Timothy Mann, and Pushmeet Kohli. On the effectiveness of interval bound propagation for training verifiably robust models. arXiv preprint arXiv: 118801127111 2018.
Minghao Han, Lixian Zhang, Jun Wang, and Wei Pan. Actor-critic reinforcement learning for control with stability guarantee. IEEE Robotics and Automation Letters, 5(4):6217-6224, 2020.
Hassan K Khalil. Nonlinear control, volume 406. Pearson New York, 2015.
J Zico Kolter and Gaurav Manek. Learning stable deep dynamics models. Advances neural
information processing systems, 32, 2019.
Reza Lavaei and Leila J Bridgeman. Systematic, lyapunov-based, safe and stabilizing controller synthesis for constrained nonlinear systems. IEEE Transactions on Automatic Control, 2023.
Boyi Liu, Qi Cai, Zhuoran Yang, and Zhaoran Wang. Neural trust region/proximal policy optimization attains globally optimal policy. In neural information processing systems, 2019.
Antonin Raffin, RI baselines3 zoo. https://gi http://gttttu.... www.//DL-RR///1-baaeeeellllllll- weweeliiieeeeee -00 2020.
Hadi Ravanbakhsh and Sriram Sankaranarayanan. Learning control lyapunov functions from coun-
terexamples and demonstrations. Autonomous Robots, 43(2):275-307, 2019.
Rosana CB Rego and Fábio MU de Araújo. Learning- robust neuro-control: A method to
compute control Yyapunov functions. International Journal of Robust and Nonlinear Control, 32 (5):2644-266) 2022.
11

---

Spencer M Richards, Felix Berkenkamp, and Andreas Krause. The Iyapunov neural network: Adaptive stability certification for safe learning of dynamical systems. In Conference on Robot Learning, pages 4666476 PMLR, 2018.
John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017.
Sumeet Singh, Spencer M Richards, Vikas Sindhwani, Jean-Jacques E Slotine, and Marco Pavone. Learning stabilizable nonlinear dynamics with contraction- based regularization. The International Journal of Robotics Research, 40(10-11): 123-1 150, 2021.
Jarrod M. Snider. Automatic steering methods for autonomous automobile path tracking. Robotics Institute, Pittsburgh, PA, Tech. Rep. CMU-RITR-09-08, 2009.
Russ Tedrake. Underactuated roboties: Learning, planning, and control for efficient and agile machines course notes for mit 6.832. Working draft edition, 3:4, 2009.
Kaidi Xu, Zhouxing Shi, Huan Zhang, Yihan Wang, Kai- Wei Chang, Minlie Huang. Bhavya Kailkhura, Xue Lin, and Cho-Jui Hsich. Automatic perturbation analysis for scalable certi- fied robustness and beyond. Advances in Neural Information Processing Systems, 2020.
Ruikun Zhou, Thanin Quartz, Hans De Sterck, and Jun Liu. Neural Iyapunov control of unknown nonlinear systems with stability guarantees. In Neural Information Processing Systems, 2022.
12

---

A Further Details about Experiments
A.1 Dynamics
Since these are described as continuous-time domains, we set the time between state updates to be 0.05s in all cases for our discretized versions of these systems.
Inverted Pendulum We use the following standard dynamics model for inverted pendulum (Chang al. 2019 Richards et al. 2018 Zhou et al. 2022]:
i mggssn(( sin + 60
mf²
with gravity constant 9.81, ball weight m = 0.15, friction b = 0.1, and pendulum length I = 0.5.
We set the maximum torque allowed to 6N m, meaning u E [-6.0, 6.0). and is applied to all experiments except for the NLC (free) case.
Path Tracking
Our path tracking dynamics follows Snider 2009
Chang et 2019
V cos (fc) éra
éra usin (0.) ,0 tan(6) uk((s cos (0c)
L
(8)
where the control policy determines the steering angle era is the distance error and 0g is the angle error. We can simplify the dynamics by defining tan( that non-linearity of dynamics is only in terms of state variables. Since in practice |6| S 40", we set u E (lan(-40%) lan(40)) )]; this applied to all experiments. We set speed 2.0, while the track circle with radius 10.0 (and, thus, x 0.1) and L, 1.0.
For RL in path tracking, the reward function is - to (c?a + 02) and we use PPO Schulman et al. 2017]. training over 50K steps.
Cartpole For cartpole we use the dynamics from Tedrake
2009
The original dynamics is
1
#
S m, sin 0(102 cos cos00) 0)
me + my sin"
0 1(m + mp sin 2 0)
fu cos 8 m1160 cos 0 sin O (m. mp)g sin
We change the variable and rename 0 as 0, so that in setting, 0 represents the angle between the pole and the upward horizontal direction, the pole angle is positive if it to the right. The purpose of this transformation is that the equilibrium point is the origin.
Using the second order chain rule
d
d2
du
2
dy duu
d
du
da
du d.22
let y = π 0,2 = t, 0, we have the new dynamics as
[ + mp sin 0(102 - cos cos 0)
H me mp sin 2
(8)
0 1(m + mn sin" 0)
fa cos 12, cos 0 sin 0 + (m, + m,)g sin
(9)
Here m. 1.0 is the weight the cart, m, position of the cart, / = 1.0 is the length of the pole, 9.81 is the gravity constant and far is the controller (force applied to the cart). Furthermore, the maximum force allowed the cart is set be 30.., meaning fo E |-30, 30]; this is applied to all experiments.
is the weight of the pole, x is the horizontal
13

---

PVTOL For PVTOL, our dynamics follows Singh et al. [2021]. Specifically, the system has the following state representation: x Pxx Ps, Uz, U2 Q, Φ (Px: Pz) and (0z, Us) are the 2D position
and velocity, respectively, and (Ø) σ) are the roll and angular rate. Control u E 1 corresponds to the controlled motor thrusts. The dynamics f(z, u) are described by
U2 cos Φ - 12 sin Φ
U2 sin Us Φ
0
&(t)
0 9 sin の ux g 008 Ø
0
u
(1/m)
(1/m) (-1/J)
1/J
where g = 9.8 is the acceleration due to gravity, m = 4.0 is the mass, = 0.25 is the moment-arm of the thrusters, and 0.0475 is the moment of inertia about the roll axis. The maximum force allowed to both ii and is set to be 39. 39.2N, which is mg. meaning "1, U2 E [0, 39.2]; this is applied all experiments.
A.2 Computation of the LQR Solution
Given the system dx/di = Ar + Bu, the objective function of the LQR is to compute the optimal controller Kr that minimizes the quadratic cost J0 (a' Qa: + u Rot)dt, where R and Q are identity matrices. The details of the system linearization as well the LQR solution are described below.
Inverted Pendulum
We linearize the system as
A =
ml2
The final LQR solution is u -1.977252340 0.97624064M.
Path Tracking We linearize the system as
0
0
0.04
The final LQR solution is u' 2u 0.1 -Era 2.196425720.
Cartpole We linearize the system
0
0
07
0
0
0.98
0
1
B
0
0
0
1
0
0
0
10.78
0
The final LQR solution is u = x + 2.41094614 34.362039470 + 10.70094483/
PVTOL
We linearize the system
0
0
07
0
0
0
0
0
1
0
0
0
0
0
0
4
0
0
0
0
B
0
0
0
0
0
0
0
0
0
0
The final LQR solution is u1, 0.70710678r
0.70710678y 5.039548718 1.10781077
1.82439774g 1.20727555d, and -0.70710678
0.70710678g
+ 5.039548718
1.10781077à 1.82439774j 1.207275550.
14

---

A.3 PPO Training for Path Tracking
We use the implementation in Stable-Baselines [Rafin 2020 for PPO training. The hyperparameterss of are fine-tuned using their auto-tuning script with budget of 1000 trials with a maximum 50000 steps. The final hyperparameters are in Table
Table 3: Hyperparameter for PPO (Path Tracking)
Parameter
Value
policy __timesteps
MIpPolicy !!float 100000
batch size
32
n_steps
64
gamma
0.95
learning_ rate
0,000208815
ent_ _nt_coef
1.90E-06
elip_range
0.1
n_epochs
10
gae_lambda
0.99
max max_grad_nme norm
0.8
vf_coef activation fn log std_ssd_m_ initt ortho artho_init init net arch vf
0.550970466
Rn.ReLU
-0.388380542
False
[8, 8]
[8, 8]
sde sample_freeee
128
A.4 Benchmark Models
When relevant, benchmark models (as well as our approach) are run for 10 seeds (seed 0 to 9). Mean + standard deviation are reported in Table and Table For SOS benchmark, we use polynomials of degree < for all environments.
NLC and UNL For Cartpole and PVTOL, we train against diameter 1.0 under the le norm for NLC and UNL. The main reason for the failure of training is that the dReal verifier becomes incredibly show after some iterations. For example, for seed 0, NLC only finished 167 certifications within hours limit for Cartpole, and 394 certifications within 24 hours limit for PVTOL.
LQR and SOS For Cartpole environment, we first certify against target region 0.1 S |||| S 0.16 for LQR benchmark where dReal returns the counterexample; we later attempt to certify against 0.1 < ||>||/o < 0. 15, dReal did not return any result within the 2-hour limit. For SOS benchmark.
dReal verifier did not return any result within the 2 hours limit against target region 0.1 S || S 0.5. For PVTOL environment, LQR benchmark returns the counterexample against target region
0.1 < || < 0.5 after 10 hours, then attempt to certify against 0.1 < || < 0.4 and fail to return any result within the remaining hours time limit. For SOS benchmark. dReal verifier did not return any result within the 24 hours limit against target region 0.1 < ||0||00 0.5.
A.5 Bounding the Functions
There are three types of bounding in our paper: 1) bound the dynamics fæx u) with linear upper/lower bound functions for verification; 2) bound the dynamics with constants for the calculation of M for ReI.U; 3) bound the dynamics with constants for the calculation of ROA. Note that the bounding for 2) and 3) are similar, except that for 2) is calculated for each subgrid, and for 3) is calculated for the entire valid region.
15

---

Inverted Pendulum Split the region: we split the region over the domain of 0, please refer to our code for the details of the splitting.
(1) Bound for verification: since the only non-linear part of the function is sin(O), we use the
M upper/lower bounds in a, B.CROWN. (2) Bound for M: |ö] S 4 + |l + ml where max
is calculated using IBP bound. (3) Bound for ROA: (x; u) loo y max(7; mgg+ fo + dl, where |d] max is the maximum force allowed.
ml
Path Tracking Split the region: we split the region over the domain of de. please refer our code for the details of the splitting.
(I) Bound for verification: for éNa since the only nonlinear function is sin(O), we use the bound in
o, B-CROWN; for 0g we bound the term vx((s cos(O.) which concave in interval x- /2, π 1/2). and
1-érar(()]
convex otherwise. (2) Bound for M: |éral S v, 10el I e/ulme + R- where |24 |mmu is calculated
V
L
v. using IBP bound. (3) Bound for ROA: [x u)ll = max(v, elwlms + R- ddt, where |u max is the maximum force allowed.
Cartpole Split the region: we split the region over the domain of 0 and á with interval 0.1 when either 0 or 0 is greater than 0.1. In case when both 0 and 0 are smaller than 0.1, we split 0, into intervals of 0.05. Please refer to our code for the details of the splitting.
(1) Bound for verification: we split the dynamics into six functions and bound them separately for both one variable and multiple variables: f1 m.tmy g f2 +m sin" fss -mm, sin cos 0
cos @
t sin
/4 (m--mp)e 29 sin
m.tm, sin So 10, cos 0 sin
mylsin 062
7mm++m, sin 85 (2) Bound for M: [ S |1/max
mp sin 85.
( max +
myg j0 S +(/u/ + I 0l + (me + mp)g) where |u/max is calculated using IBP my. 1
2
bound. (3) Bound for ROA: [I, o) + max(7, |u/max + +ppl?
2 [[ + mpdr?
2
(me + dt, where |u| mmx is the maximum force allowed.
PVTOL Split the region: we split the region over the domain of 0 with interval 0.25 and split the region over the domain of àt, U, f with interval 0.5. Furthermore, in the MILP we add a constraint to ensure the loo norm of state value is no less than I. Please refer to our code for the details of the splitting.
(1) Bound for verification: split the dynamics into three functions and bound them separately for multi-variables fi = x cos fa = x sin y. fa = ry. (2) Bound for M: we bound each term of the dynamics with the min/max values using the monotonic property. (3) Bound for ROA:
[x u)ll = + max(y, x2 + 2/w/am 27, allowed.
2 wess dtt, where |m/max is the maximum force
m
s
A.6 Speeding up Verification
When there is a limit on parallel resources, so that verification must be done serially for the subgrids k. we can reduce the number of grids to verify in each step follows. The first time verifier called, we perform it over all subgrids; generally, only a subset will return a counterexample. Subsequently, we only verify these subgrids until none return counterexample, and only then attempt full verification.
A.7 ROA Calculation
For Inverted Pendulum and Path Tracking, we approximate ROA using grid of 2000 cells along each coordinate. For Cartpole, we use a grid of 150 cells along each coordinate, and for PVTOL, we grid of 50 cells along each coordinate.
For our baselines where the systems are continuous, the set levels are {a E R V(x)) < p) where p minzeoR V(r). Note that for discrete-time systems, having p minzeaR V(x) is not sufficient.
For example, assume R is the valid region, any point outside of R. satisfy V(f(a, u)) > V(x), and
Wx E ÖR, (zz p. Assume we have R ÖR and f(z, u) y. y & R (this is possible since
16

---

the system is discrete). Since is in the valid region, we have V④y) < V((x) p. However, since no longer within the valid region, assume (y,, u)) p (yy), the system will not converge to equilibrium point starting from point y. This means (xr R V ((x) < p} is not subset of ROA.
B Computing Linear Bounds on Lipschitz--Continuuss Dynamics
Theorem B.1. Given function g(x) R" + where the Lipschitz constant is X, we have g(x) _ g( 22,,,.. + X Σ, (I, z1,1) and (2.1,,,, X (2 21,) g(z),
where Dj E Ti,l E [n].
Proof. For any E [m]. given z; (a,,,, for arbitrary fixed values along dimension [n] \ {j) (we denote them as x i), we have (g( - g( it x--i] & Alas
vi.tl.
we show that
S ( + X(x; xz.1)
if g( x-1) 2 ( 2):,, then ( x.,)) - 1 g( S A(z; -2,,), meaning g( $ [z. x-5) + M(a; 2(,1).
if g(1,, < g( x-1), since M(x, - x2,1) 2 0, we have g(x;, S g( + A(z; Tt,)).
we show that g( -1i) M(1a x,,)) S g(11,--1)
if g( x-1) 2 g x-i,, then since - -M(x; - -1i,) S O, we have g
Maz {i,) S g( T-)).
- glg, < g( xT then g( x- g( x-) S A(xa Tt,)), meaning g(T1,1 A(x; x4,1) S g(
By applying the above result along all n dimensions recursively, we have the upper bound as g(x) S (211, En.) + Σ, X(z; 24,1), and the lower bound as g( .2. (n.i)
Σ, A(x: di,i) S g(x).
C Ablation Study
In this section, conduct ablation study the MILP solver. Table 4 are the results for DITL- dReal, in which we use our framework, but replace all the MILP components with dReal. As we can see, using MILP is an essential ingredient in the success of the proposed approach:
Table 4: Ablation study on MILP solver.
Environment
Runtime
ROA
Max ROA
Success Rate
Inverted Pendulum (DITL-dReal)
6.0 1.7 (8)
57 24
75
100%
Inverted Pendulum (DITL, main paper)
8.1 + 4.7(s)
61 - 31
123
100%
Path Tracking (RI., DITL-dReal)
600 + 0 (s)
+
0
0%
Path Tracking (RT., DITL,, main paper)
14 + 11 (s)
+ 3.5
16
100% %
Path Tracking (LQQR, DITL-ddRall)
420.9 + 182 (s)
+
11
60%
Path Tracking (LQQR, DITL, main paper)
9.8 + (s)
8 + 3
12.5
100 %
Cartpole (DITL-dRcal)
>2 (hours)
N/A
N/A
0%
Cartpole (DITL, main paper)
0.9 + 0.3 (hours)
0.021 + 0.012
0.045
100 %
PVTOL (DITL-dReal)
>24 (hours)
N/A
N/A 0.028
0%
PVTOL (DITL, main paper)
13 - (hours)
0.011 + 0.008
100 %
17


Reviews:
[{"review_id": "z6oH3z5mih", "reviewer": "Reviewer_NZiv", "summary": "This paper presents an algorithmic framework for learning neural Lyapunov control with stablizing control policy for discrete-time dynamic systems. It proposes to verify the stability condition with mixed-integer linear programming and speed up Lyapunov function learning with  gradient based approximation. On four benckmarks, it outperforms SOTA baselines in terms of running time, size of regions of attraction and successful rate.", "questions": "- The paper provides 2 options for setting the bias term in the training loss function. However, it is not crystal clear the impact of this term over the learning.  An ablation study over the bias setting could help. \n\n- To better motivate the use of mixed-integer linear programming for stability verification in discrete-time systems, it would be good to show a case by comparing DITL-MILP with DITL-dReal. \n\n- Does the run time include the initialization time (for example by LQR)? ", "limitations": "- The method is built upon the assumption that both the verification function and control policy can be simply represented with MLPs with ReLU activations. \n", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is well structured and nicely presented\n- The approach has extended the learning of Lyapunov functions to discrete-time dynamic systems, and has been shown effective on 4 benchmarks with a significant improvement. \n", "weaknesses": "See questions and limitations. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's thoughtful comments and suggestions.  Our detailed responses are below.\n\n1. ***Comment:*** *The paper provides 2 options for setting the bias term in the training loss function. However, it is not crystal clear the impact of this term over the learning. An ablation study over the bias setting could help.*\n\n    ***Response:*** We appreciate this suggestion.  In our method the bias term is treated as a hyperparameter for learning, thus the one presented in the main paper is the optimal one. Below is the ablation result for using the alternative method (i.e., without bias if the paper uses bias, or with bias when the results in the paper do not):\n\n\n| Environment |  Runtime  |ROA |Max ROA |Success Rate |\n| -------- | -------- | -------- | -------- |-------- |\n| Inverted Pendulum (without bias)         | 11.1 $\\pm$ 0.2  (s)   |45 $\\pm$ 17     |72 |100\\%     |\n| Path Tracking (RL, without bias)         | 17.8 $\\pm$ 0.1  (s)   | 8 $\\pm$ 0     |8 |100\\%     |\n| Path Tracking (LQR, without bias)       | 12.1 $\\pm$ 2.9  (s)   |9 $\\pm$ 2     |11 |100\\%     |\n| Cartpole (with bias)       | >2 (hours)    | N/A  |N/A | 0\\%     |\n| PVTOL (with bias)        | 17 $\\pm$ 10 (hours)  |   $\\sim$ 0 $\\pm$ 0 |0.0001 |20\\%     |\n\n$\\ $\n\n2. ***Comment***: *To better motivate the use of mixed-integer linear programming for stability verification in discrete-time systems, it would be good to show a case by comparing DITL-MILP with DITL-dReal.*\n\n    ***Response:*** Thank you for the suggestion!  Below are the results for DITL-dReal, in which we use our framework, but replace all the MILP components with dReal.  As we can see, using MILP is an essential ingredient in the success of the proposed approach:\n\n| Environment |   Runtime  |ROA|Max ROA |Success Rate |\n| -------- | -------- | -------- | -------- | -------- |\n| Inverted Pendulum (DITL-dReal)         | **6.0 $\\pm$ 1.7** (s)  | 57 $\\pm$ 24  |75|**100\\%**     |\n| **Inverted Pendulum (DITL, main paper)**        |   8.1 $\\pm$ 4.7(s)  |  **61 $\\pm$ 31** |**123**|  **100\\%**  |\n| Path Tracking (RL, DITL-dReal)        | 600 $\\pm$ 0 (s)   | 0 $\\pm$ 0  |0 | 0\\%     |\n| **Path Tracking (RL, DITL,  main paper)**         | **14 $\\pm$ 11** (s)  | **9 ± 3.5** |**16**|   **100\\%** |\n| Path Tracking (LQR, DITL-dReal)         | 420.9 $\\pm$ 182 (s) | 4 $\\pm$ 4   |11 |60\\%     |\n| **Path Tracking (LQR, DITL,  main paper)**         | **9.8 $\\pm$ 4** (s)  | **8 ± 3**  |**12.5**|  **100\\%**  |\n| Cartpole (DITL-dReal)       | >2  (hours)   | N/A  |N/A | 0\\%     |\n| **Cartpole (DITL, main paper)**       | **0.9 $\\pm$ 0.3**  (hours) | **0.021 $\\pm$ 0.012** | **0.045** | **100\\%**     |\n| PVTOL (DITL-dReal)       | >24  (hours)   |  N/A   |N/A|0\\%     |\n| **PVTOL (DITL, main paper)**       | **13 $\\pm$ 6** (hours)  | **0.011 $\\pm$ 0.008** | **0.028** | **100\\%**     |\n\n$\\ $\n\n3. ***Comment:*** *Does the run time include the initialization time (for example by LQR)?*\n\n    ***Response:*** Run time does not include initialization time as it is negligible in most cases. The run time for LQR for all instances is less than 0.1s. One exception is the run time for RL (path tracking), which is 154.96 seconds as  documented in the results section (Line 330).\n\n$\\ $\n\n4. ***Comment:*** *The method is built upon the assumption that both the verification function and control policy can be simply represented with MLPs with ReLU activations.*\n\n    ***Response:*** This is a good point. Nevertheless, we do not view this assumption as very limiting for our purposes, since NN with ReLU activation layers can approximate arbitrary continuous functions. Furthermore, policies with such structure are common in deep reinforcement learning. In any case, we will comment on this in the revision."}, {"review_id": "l8FPXR7WYY", "reviewer": "Reviewer_qMzc", "summary": "This paper studies learning stabilizing controllers with neural Lyapunov functions in discrete-time nonlinear systems. Previous works on neural Lyapunov control focused on continuous-time systems, and this is the first work that studies discrete-time systems. The outline of the proposed approach is similar with prior works on continuous-time systems, while the authors propose several key novel techniques for each step: They only use ReLU activations in the neural Lyapunov functions so that the constraints can be verified by mixed-integer programming. The authors also proposed a heuristic way to generate counterexamples efficiently via projected gradient descent. The simulation results show that the proposed method can find larger ROA and be implemented efficiently in several examples.", "questions": "Please see my comments in the previous section.", "limitations": "I didn't see the discussion about limitations or future directions in my read. But I don't think there is any potential negative societal impact.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Learning neural Lyapunov control in discrete-time systems is important for many applications. This work proposes many novel techniques to handle the difficulties that have been found in prior works. The simulation results are promising.", "weaknesses": "I have some minor concerns about the presentation. After reading this paper, it is not clear to me which challenges are specific to the discrete-time system.\n\nIn the first part of verifying Lyapunov constraints, the authors take the motivation from SMT solvers. However, this is also a challenge for neural Lyapunov control in continuous-time systems. I hope the authors can add a discussion about what efforts have been made to improve efficiencies in the continuous-time systems and why we cannot directly apply them here. And for the proposed ReLU network plus mixed-integer programming, does this method also apply to continuous-time systems? Or is it limited to discrete-time systems because of some special properties?\n\nThere is a similar question for proposed heuristic for generating counterexamples. Is this heuristic limited to discrete-time system? If not, can we replace the counterexample generation step in existing algorithms for continuous-time system and see if the performance improves?\n\nAnother minor concern I have with the PGD heuristic is that it might fail to find a counterexample because gradient-based search can get stuck in local minimums. Thus, in the final step when the algorithm claims success, shall we use a less efficient way to verify the Lyapunov conditions are actually satisfied?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's thoughtful and detailed comments and questions, to which we respond in detail below.\n\n1. ***Comment:*** *It is not clear to me which challenges are specific to the discrete-time system.*\n\n    ***Response:*** There is a consequential distinction in the definition of Lyapunov functions between continuous and discrete-time settings.  In continuous-time settings, the key condition is $\\forall x \\in D; \\nabla_{f_u} V(x) \\leq 0$, where $\\nabla_{f_u} V(x)$ is Lie derivative of function $V$ w.r.t. vector field $f_u$ at point $x$.  In discrete-time settings, it is: $\\forall x \\in D; V(f(x)) - V(x) \\leq 0$. If $V$ is a neural network with ReLU activations, it is not continuously differentiable, which makes it non-obvious how we can construct a MILP to verify the condition on the Lie derivative uniformly over the state space that is needed in continuous-time systems. However, this issue does not come up for the Lyapunov condition in discrete-time systems. An alternative could be to replace ReLU units with continuously differentiable neurons, but these must now be approximated in an MILP formulation using piecewise linear functions, which is a major challenge to scalability and efficacy, as it introduces an approximation gap that makes the Lyapunov condition even more difficult to achieve as we get closer to the origin. Another alternative is to use dReal for verification, but as we demonstrate above (see the data in response to Reviewer NZiv, Comment 2), MILP is in fact a critical ingredient in enabling high efficacy of the proposed approach. \n  \n$\\ $\n\n2. ***Comment:*** *In the first part of verifying Lyapunov constraints, the authors take the motivation from SMT solvers. However, this is also a challenge for neural Lyapunov control in continuous-time systems. I hope the authors can add a discussion about what efforts have been made to improve efficiencies in the continuous-time systems and why we cannot directly apply them here.*\n\n    ***Response:*** Thank you for this suggestion! We note that our baselines already include state of the art approaches for continuous-time neural Lyapunov control; these still make use of dReal, and no one has proposed effective heuristic techniques for improving verification during training for continuous-time neural Lyapunov control.  Indeed, the conceptual advances we make, such as heuristic verification techniques during training, could in principle be considered also in continuous-time control settings. However, it does not seem evident how to develop MILP-based approaches in that setting (see our response above), and dReal remains a significant bottleneck.\n    \n$\\ $\n\n3. ***Comment:*** And for the proposed ReLU network plus mixed-integer programming, does this method also apply to continuous-time systems? Or is it limited to discrete-time systems because of some special properties?\n\n    ***Response:*** As we discuss in our response to Comment 1 above, it is not clear how to adopt our MILP approach to verify Lyapunov conditions in continuous-time systems.\n\n$\\ $\n\n4. ***Comment:*** There is a similar question for proposed heuristic for generating counterexamples. Is this heuristic limited to discrete-time system? If not, can we replace the counterexample generation step in existing algorithms for continuous-time system and see if the performance improves?\n\n    ***Response:*** We believe that the proposed heuristic counterexample generation approach would indeed extend to continuous-time systems. However, we wish to emphasize that the key bottleneck for continuous-time systems remains the scalability of dReal. For instance, even in verifying LQR solutions (which are easy to compute), dReal was unable to complete verification within the specified time frame for both the cartpole (2 hours) and PVTOL (24 hours) control systems. This implies that even if we were to apply the suggested heuristic to continuous-time systems, dReal would still encounter challenges in cases like cartpole and PVTOL (our final step for stopping the algorithm is to pass the verification). It is ultimately the combination of the use of heuristic counterexample generation with MILP-based verification that enables us to make significant progress.\n\n$\\ $\n\n5. ***Comment:*** *Another minor concern I have with the PGD heuristic is that it might fail to find a counterexample because gradient-based search can get stuck in local minimums. Thus, in the final step when the algorithm claims success, shall we use a less efficient way to verify the Lyapunov conditions are actually satisfied?*\n\n    ***Response:*** We would like to clarify that we achieve provable stability precisely by resorting to full MILP-based verification whenver the PGD heuristic can no longer identify counterexamples. Moreover, the MILP-based verifier can generate a counterexample, which then generates additional training iterations (e.g., often PGD now again generates counterexamples once this one is added to the training dataset and can be used to bootstrap the heuristic). Consequently, there are typically several times during training that we have to call the MILP-based verifier. We only claim success when training returns with MILP proving that stability conditions hold. The soundness of our approach is formally asserted in Theorem 4.1, and holds by construction.\n"}, {"review_id": "EKGZwR1mfu", "reviewer": "Reviewer_yt77", "summary": "This paper extends the neural Lyapunov control method to the discrete-time dynamical systems and propose two tricks to quickly find the counterexamples to accelerate the training process. The proposed method is tested in several standard tasks with comparison to the SOTA methods.", "questions": "[1] Does the proposed method scale to the high-dimensional systems such as controlling the molecular dynamics?\n\n[2] The neural Lyapunov method for unknown system (UNL) focuses on the realistic model-free settings, while the proposed method relies on the fully known dynamics, is the comparison in terms of size of ROA fair?\n\n[3] The problematic Lemma 3.1 cannot guarantee the soundness of the proposed method.   ", "limitations": "The authors only test the proposed method in toy models; however, they do not show the current method can be applicable in the real-world scenarios.", "rating": 2, "confidence": 5, "soundness": 1, "presentation": 2, "contribution": 1, "strengths": "[1] The proposed mixed-integer LP and PGD method can efficiently find the counterexamples in the training process.\n\n[2] The proposed method can find the largest ROA compared to the existing methods.\n\n\n", "weaknesses": "[1] The extension of neural Lyapunov control from continuous-time dynamical systems to discrete-time dynamical systems has little contribution to the existing neural control community because the Lyapunov theories for continuous/discrete system are essentially the same.\n\n[2] The current method has huge computational cost for the case where the high-dimensional tasks are taken into account.\n\n[3] Lemma 3.1 is NOT mathematically correct.\n\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful comments. We respond to these below.\n\n1. ***Comment:*** *The extension of neural Lyapunov control from continuous-time dynamical systems to discrete-time dynamical systems has little contribution to the existing neural control community because the Lyapunov theories for continuous/discrete system are essentially the same.*\n\n    ***Response:*** As we highlight in our response to Reviewer qMzc (in particular, our response to Comment 1), there is a very consequential difference in the Lyapunov conditions between continuous-time and discrete-time systems. It is this difference that enables us to exploit the structure of the condition in discrete-time systems to enable a substantial increase in efficacy that we exhibit in the experiments.  We also wish to emphasize that stability for continuous-time systems does not imply stability in discrete-time systems, nor is the converse true. For example, in our inverted pendulum experiment, all the the continuous-time controller and Lyapunov function generated by the NLC method fails the discrete-time Lyapunov condition. Moreover, as neural network control takes non-negligible computational time in a real robotic system, the controller is effectively discrete time, and stability guarantees for continuous-time systems would in general be unsound in such settings. \n    \n$\\ $\n\n2. ***Comment:*** *The current method has huge computational cost for the case where the high-dimensional tasks are taken into account. Does the proposed method scale to the high-dimensional systems such as controlling the molecular dynamics?*\n\n    ***Response:*** Today, no method exists that can work on systems of such scale with general non-linear dynamics and provable stability guarantees. As our experiments demonstrate, our approach is an order-of-magnitude improvement over the state of the art in higher-dimensional settings, but clearly there is much research that remains to further improve scalability. The key bottleneck in this setting remains sound verification of stability.\n\n$\\ $\n\n3. ***Comment:*** *The neural Lyapunov method for unknown system (UNL) focuses on the realistic model-free settings, while the proposed method relies on the fully known dynamics, is the comparison in terms of size of ROA fair?*\n\n    ***Response:*** We appreciate this point. We wish to clarify that the comparison is not intended to make claims about efficacy of UNL, which is indeed targeted at a far more challenging problem; UNL simply serves as another state of the art baseline. Not including it would have exposed us to a criticism that UNL is a more recent approach and could be a more competitive baseline than others (which, as we show, is not the case).\n\n$\\ $\n\n4. ***Comment:*** *The authors only test the proposed method in toy models; however, they do not show the current method can be applicable in the real-world scenarios.*\n\n    ***Response:*** We observe that both of our baseline papers, UNL and NLC, utilize such control dynamics. Indeed, the benchmarks we use are common within the control community. Our aim is for our work to serve as a stepping stone towards practical application in real-world scenarios."}, {"review_id": "vSaQp4K9AO", "reviewer": "Reviewer_ruNW", "summary": "The stability of nonlinear systems has long been a challenge in the field of control systems, with current methodologies employing Lyapunov stability theory to derive control policies. However, finding suitable Lyapunov functions for these systems is notably complex. To address this, researchers have started using neural networks to approximate Lyapunov functions, though primarily for continuous-time systems. The paper introduces a first-of-its-kind approach for learning neural Lyapunov control in discrete-time systems. This approach has three critical components:A unique mixed-integer linear programming method to verify stability conditions.A new method for computing sub-level sets to identify the region of attraction.A heuristic gradient-based method to rapidly find counterexamples, which improves the learning process of Lyapunov functions.Experimental results show a significant improvement over current techniques, outperforming neural Lyapunov control baselines by an order of magnitude in both running time and the size of the region of attraction. For the 'cartpole' and 'PVTOL' benchmarks, this is the first automated method to yield a provably stable controller.", "questions": "How does your method compare in terms of advantages to similar works like \"Reinforcement Learning Control of Constrained Dynamic Systems with Uniformly Ultimate Boundedness Stability Guarantee\" by Han M, Tian Y, Zhang L, et al. (Automatica, 2021, 129: 109689)?", "limitations": "This method leverages a mixed-integer linear programming approach, which is known to be NP- hard.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This method outperform recent neural Lyapunov control baselines by an order of magnitude in both running time and the size of the region of attraction.", "weaknesses": "The loss function is then a weighted sum of several terms, which is hard to balance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the comments. Our detailed responses are below.\n\n1. ***Comment:*** *The loss function is then a weighted sum of several terms, which is hard to balance.* \n\n    ***Response:*** It is common for loss function to be the weighted sum of several terms, as those weights are hyperparameters for learning we need to calibrate. For example, in both of our baselines, NLC and UNL, their loss function are the weighted sum of several terms. Similarly, standard RL approaches such as PPO also involve many hyperparameters.\n\n    $\\ $\n\n2. ***Comment:*** *How does your method compare in terms of advantages to similar works like \"Reinforcement Learning Control of Constrained Dynamic Systems with Uniformly Ultimate Boundedness Stability Guarantee\" by Han M, Tian Y, Zhang L, et al. (Automatica, 2021, 129: 109689)?*\n\n    ***Response:*** We would like to clarify that the referenced paper is addressing a fundamentally different problem from ours. Although they also refer to a notion of stability, it is quite distinct, and is more useful for providing safety than stability guarantees. For example, in our paper we want the cartpole to converge to and remain at the origin, while in their experiment they want the carpole position to remain within a target safe interval ($x<4$). \n    \n    $\\ $\n\n    Due to the formal distinction, the goal of Han et al. is to achieve what is commonly referred to as safety, rather than stability in the Lyapunov sense. We would also like to point out that  Han et al. does not provide any strict verification guarantees like ours (for example, their Theorem 1 only guarantees safety approximately), which makes their evaluation essentially empirical, while our approach returns a controller that is provably stable. Moreover, while Han et al. also utilize Lyapunov functions, the associated conditions are in expectation, rather than uniform over all states, as required in the notion of stability that we use.\n    \n    $\\ $\n\n    To be more precise, in Han et al. the key concept is uniformly ultimate boundedness (UUB), which is defined as follows (Defn. 1 on page 3; also see Thowsen, 1983): A system is said to be uniformly ultimately bounded (UUB) with ultimate bound $\\eta$, if there exist positive constants $b$, $\\eta$, where $\\forall \\epsilon < b$ there exists $T(\\epsilon, \\eta)$, such that $\\|x_{0}\\| < \\epsilon$ $\\Rightarrow$ $\\|x_t\\| < \\eta$, $\\forall t > T(\\epsilon, \\eta)$. At the high level, this means that trajectories that start from inside the safe set (defined by $\\|x_{0}\\| < \\epsilon$) will remain inside a subset of the safe set (defined by $\\|x_t\\| < \\eta$), for suitably chosen $b$ and $\\eta$.\n    \n    $\\ $\n\n    In contrast, the stability notion we are concerned with is defined as follows: a system is stable within a region of attraction $\\mathcal{R}$ (which includes the origin) if $\\forall x_0 \\in \\mathcal{R}$, $\\lim_{t\\rightarrow\\infty} x_t = 0$, where $x_{t+1} = f(x_t)$ (we omit the control input here to simplify discussion). At the high level, we require convergence to the origin from any starting point in the region of attraction, which we also aim to make as large as possible.\n\n    $\\ $\n\n3. ***Comment:*** *This method leverages a mixed-integer linear programming approach, which is known to be NP-hard.*\n\n    ***Response:*** Tight verification is typically NP-hard, so that the main avenues for advances in scalability are either in the direction of relaxation (which loses tightness) or taking advantage of problem structure. MILP is an approach that takes advantage of any linearity structre in the problem, and MILP solvers can indeed scale well, although they, too, have their limits. While clearly there is much research that remains in further scaling the learning of provably stable control, we note that our approach is an order-of-magnitude advance in efficacy/scalability over prior art for the more challenging problems."}, {"review_id": "fkZA9y01YT", "reviewer": "Reviewer_KEba", "summary": "The paper proposes a Lyapunov control method for discrete-time systems, in contrast to previous approaches targeted at continuous-time systems. The paper proposes a mixed-integer linear programming approach for verifying stability conditions in discrete-time systems, a technique for computing sub-level sets which define the region of attraction, and a heuristic gradient-based approach for finding counterexamples that can accelerate the learning of Lyapunov functions.", "questions": "See weaknesses.", "limitations": "See weaknesses.", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "I appreciate the counterexample generation method proposed in this work. It leverages the key ideas in adversarial training in the computer vision domain to use gradients to accelerate the discovery of counterexamples. ", "weaknesses": "- The Heuristic Counterexample Generation cannot guarantee all counterexamples are found. It relies on sampling in the state space, but sampling, the method itself, cannot exhaust the state space. Therefore, the verification is more like an empirical evaluation. I am not saying being “empirical” is not good, but it just simply does not match the verification objective which should be very strict.\n\n- The authors did not explain if it is truly challenging to train in discrete time system compared to training in continuos time system. In fact, the continuous time system and discrete time system can share the same framework in training. See https://arxiv.org/pdf/2101.05436.pdf how they computed dot{h} in the paragraph below (7).\n\n- I would like to discuss with the authors if verification is truly necessary in Learning for Control. As we know, real-world dynamics are much more complex than the dynamics in the neural Lyapunov studies. The verification methods presented in these papers are very difficult to be applied to the real world and make real impacts. Often, a method that pursues excellent empirical performance instead of 100% rigorous verification will be closer to real-world scenarios in control (other domains such as planning or chip design might be different).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and engaging comments. Below we provide our responses.\n\n1. ***Comment:*** *The Heuristic Counterexample Generation cannot guarantee all counterexamples are found. It relies on sampling in the state space, but sampling, the method itself, cannot exhaust the state space. Therefore, the verification is more like an empirical evaluation.*\n\n    ***Response:*** This appears to be a misunderstanding. In fact, our approach uses both heuristic and sound counterexample generation (using MILP). It is this combination that enables both improved efficacy, as well as soundness of the overall approach in that in the end the Lyapunov conditions are provably satisfied.  In particular, the controller that our algorithm returns is always verified to be stable using the sound MILP verifier that we describe, and it is the stopping criterion of our algorithm.  Consequently, as we assert in Theorem 4.1, our algorithm always returns provable stable control (in the sense of satisfying Lyapunov stability conditions) by construction.\n\n     $\\ $\n\n2. ***Comment:*** *The authors did not explain if it is truly challenging to train in discrete time system compared to training in continuos time system.*\n\n    ***Response:*** As we emphasize in our Responses to Reviewer qMzc above, we can exploit structure in discrete-time settings that cannot be so easily exploited in continuous-time domains. As a result, we actually attain order-of-magnitude improvement compared to prior art, which is mostly in continuous time. Since neural network controllers take non-negligible computation to output control, systems using these are effectively discrete time, and the structure we exploit is therefore highly relevant.\n\n     $\\ $\n\n3. ***Comment:*** *In fact, the continuous time system and discrete time system can share the same framework in training. See https://arxiv.org/pdf/2101.05436.pdf how they computed dot{h} in the paragraph below (7).*\n    \n    ***Response:*** Note that in the referenced paper the authors are not doing strict verification, and their evaluation is entirely empirical. Thus, Qin et al. (the referenced paper) approximate the continuous time derivative $\\dot{h}$ numerically and still achieve good empirical results.  In contrast, we achieve fully verified stability (see our response to Comment 1 above). Since the definition of continuous and discrete-time Lyapunov stability condition are different, we cannot  use the numerical method for derivatives in the paper, as doing so will not allow the controller-Lyapunov function pair to pass our MILP verifier. More precisely, in continuous-time settings, the key condition is $\\forall x \\in D; \\nabla_{f_u} V(x) \\leq 0$, where $\\nabla_{f_u} V(x)$ is Lie derivative of function $V$ w.r.t. vector field $f_u$ at point $x$.  In discrete-time settings, the condition is: $\\forall x \\in D; V(f(x)) - V(x) \\leq 0$. Satisfying the former condition does not, in general, satisfy the latter (because discretized continuous-time systems imply that control is effectively piecewise constant, rather than continuously adapting to state feedback). \n\n    $\\ $\n\n    In addition the nature of our training framework is itself quite different from Qin et al. (and more akin to the NLC framework) in that we leverage a form of counterexample-guided abstraction and refinement (CEGAR) paradigm, iteratively generating counterexamples and adding these to training data. Our key advances within this framework involve the use of gradient-based heuristic counterexample generation and MILP-based verification, which in combination achieve significant improvements over prior art.\n\n     $\\ $\n\n    A final distinction is that Qin et al. are primarily concerned with safety (avoiding unsafe sets), while we focus on Lyapunov stability (ensuring that the system converges to the origin). \n\n     $\\ $\n\n4. ***Comment:*** *I would like to discuss with the authors if verification is truly necessary in Learning for Control. As we know, real-world dynamics are much more complex than the dynamics in the neural Lyapunov studies. The verification methods presented in these papers are very difficult to be applied to the real world and make real impacts. Often, a method that pursues excellent empirical performance instead of 100\\% rigorous verification will be closer to real-world scenarios in control (other domains such as planning or chip design might be different).*\n\n    ***Response:*** This is wonderful discussion to have! Indeed, the current state of the art of approaches that yield fully verified properties such as stability cannot be applied to the full complexity of real control problems. Nevertheless, just because we cannot do this now does not mean that we will never have this ability. We believe that our work serves as an important stepping stone in the line of research that continues to improve scalability. It can often takes several years of concerted effort to achieve a high-impact goal, and we feel that the aspiration to scalable synthesis of fully-verified control in real systems is a worthy goal to aspire to. In the meantime, we must indeed engage in a combination of empirical evaluation and testing with verification of smaller subcomponents, but the goal is to be able to verify larger and larger system components over time, until we can verify the full system."}]

Author rebuttal:

Meta-review: This paper proposes an approach to learning Lyapunov stable neural network-based controllers for discrete-time systems. The reviewers generally appreciated the novelty of the approach and the strength of the simulation results (though with some desire that they be run at larger scale). While several raised questions about the magnitude of difference between this paper's contribution in the discrete setting vs. prior papers' contributions in the continuous setting, the authors successfully argued that the distinction in the definition of Lyapunov functions between continuous and discrete-time settings yields consequential differences in the actual methodology for learning stable neural network-based controllers, and that their experimental results validate the improved performance of their method against continuous baselines. One reviewer objected to the inclusion of the condition $\|\|x\|\|_\infty \geq \epsilon$ within Lemma 3.1 on theoretical grounds, whereas the authors argued for it on numerical grounds - while this disagreement seems to stem from a difference in perspective between communities, the paper would benefit from further clarification of this as well as discussion on how $\epsilon$ is chosen in practice.
Decision: Accept (poster)