Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Authors: Haohan Zou
Jie Feng
Hao Zhao
Yuanyuan Shi
Venue: ICML 2025 poster (2025)
Keywords: Control Theory
Lyapunov Function
Symbolic Transformer
Risk-seeking Policy Optimization
Analytical Function Discovery
Stability
Primary area: applications->everything_else

Abstract: Despite advances in learning-based methods, finding valid Lyapunov functions for nonlinear dynamical systems remains challenging. Current neural network approaches face two main issues:  challenges in scalable verification and limited interpretability. To address these, we propose an end-to-end framework using transformers to construct analytical Lyapunov functions (local), which simplifies formal verification, enhances interpretability, and provides valuable insights for control engineers. Our framework consists of a transformer-based trainer that generates candidate Lyapunov functions and a falsifier that verifies candidate expressions and refines the model via risk-seeking policy gradient. Unlike Alfarano et al. (2024), which utilizes pre-training and seeks global Lyapunov functions for low-dimensional systems, our model is trained from scratch via reinforcement learning (RL) and succeeds in finding local Lyapunov functions for *high-dimensional* and *non-polynomial* systems. Given the symbolic nature of the Lyapunov function candidates, we employ efficient optimization methods for falsification during training and formal verification tools for the final verification. We demonstrate the efficiency of our approach on a range of nonlinear dynamical systems with up to ten dimensions and show that it can discover Lyapunov functions not previously identified in the control literature. Full implementation is available on [Github](https://github.com/JieFeng-cse/Analytical-Lyapunov-Function-Discovery).

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Haohan Zou
Jie Feng"
Hao Zhao?
Yuanyuan Shi
Abstract
system trajectories. Lyapunov functions also play a central role in controller design, providing formal guarantees of closed-loop stability and robustness Khalil (2002). How- ever, designing a Lyapunov function for nonlinear systems has long been considered more of an 'art' than a science, even for stable dynamics, due to its inherent complexities. Motivated by this challenge, we have witnessed great in- terest in the development of computational algorithms for Lyapunov function construction. McGough et al. (2010) employed an evolutionary algorithm for the symbolic com- putation of Lyapunov functions, but the exponential growth search space impedes its scalability. Alternatively, sum-of- squares (SOS) methods reformulate the task as a semidefi- nite program (SDP) that certifies stability with polynomial candidates (Papachristodoulluu & Prajna, 2005a:b; Ahmadi & Majumdar, 2016; Dai & Permenter, 2023). However, handling local constraints non-polynomial dynamics re- quires auxiliary variables and extra (in jqquality constraints (Papachristodoulou & Prajna, 2005a;b), which leads to scal- ability issues of the SOS methods for real-world problems. Moreover, the theoretical result Ahmadi et al. (2011) on asymptotic Lyapunov stability shows that even some very simple globally asymptotically stable dynamics may not agree with polynomial Lyapunov function of any degree.
Despite advances in learning-basee methods, find- ing valid Lyapunov functions for nonlinear dy- namical systems remains challenging. Current neural network approaches face two main issues: challenges in scalable verification and limited interpretability. To address these, we propose an end-to-end framework using transformers to construct analytical Lyapunov functions (local), which simplifies formal verification, enhances in- terpretability, and provides valuable insights for control engineers. Our framework consists of transformer-basee based trainer that generates candi- date Lyapunov functions and a falsifier that veri- fies candidate expressions and refines the model via risk-seeking policy gradient. Unlike Alfarano et al. (2024), which utilizes pre-training and seeks global Lyapunov functions for low-dimensional systems, our model is trained from scratch via re- inforcement learning (RL) and succeeds in finding local Lyapunov functions for high-dimensional and non-pone polyoolllliall systems. Given the symbolic nature of the Lyapunov function candidates, we employ efficient optimization methods for falsi- fication during training and formal verification tools for the final verification. We demonstrate the efficiency of our approach on range of nonlinear dynamical systems with up to ten dimensions and show that it can discover Lyapunov functions not previously identified in the control literature. Full implementation is available on Github.
Recent advances in deep learning have enabled data-driven neural Lyapunov function with formal verification (Chang et al., 2019; Zhou et al., 2022; Wu et al., 2023; Dawson et al., 2023b; Edwards et al., 2024; Wang et al., 2024; Yang et al., 2024). However, these methods face two key chal- lenges: 1) lack of interpretability and high verification costs (Dawson et al., 2023b). Neural networks' black-box nature limits insights into the system's dynamical behavior. Additionally, over-paramettrization and nonlinear activa- tions complicate formal verification, leading to scalability issues. Tools like SMT Chang et al. (2019), MIP Wu et al. (2023), and a, 9-CROWN Yang et al. (2024) require small, specialized networks to ensure feasible verification times.
1. Introduction
A Lyapunov function is an energy-like function used to cer- tify stability of dynamical systems. A sufficient condition for stability is that the Lyapunov function decreases along
Compared with neural Lyapunov functions, analytical Lya- рunov functions offer two distinct advantages. First, their symbolic nature offers interpretability, and provides insights for designing stability-guaranteed guaranteed control policy (Sontag, 1989 Feng et al., 2023a;b; 2024b; Cui et al., 2023b). Sec- ond, analytical functions enable efficient verification due their simplicity and symbolic structure, which seamlessly
'Equal contribution 'Department of Elecctical and Computer Engineering, UC San Diego, La Jolla, United States "Swiss Federal Institute of Technology in Lausanne. Correspondence to: Haohan Zou <hazou edu>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
We introduce the first RL-based framework for directly discovering analytical Lyapunov functions for nonlin- ear dynamical systems, bypassing the need for super- vised learning with large-scale datasets.
350
Sample Analytical Formula:
300
sin(x)] - castx;) - cos(x,) - x, - Xx) 3
250
200
We propose novel and efficient policy optimization pipeline integrating Alobal-optimizztion- punov verification, reward design for candidate Lya- punov evaluation, and risk-seeking policy gradient to optimize the symbolic transformer, trainable on ma- chines with limited computation resources.
144.216
Lya-
106.06s
0.23s
Ours
Fassil 2.0
ANLC
We demonstrate the efficiency of our method on vari- ous systems, including non-polynomial dynamics like the pendulum, quadrotor, and power system frequency control. Notably, our approach scales to a 10-D sys- tem and discovers a valid local Lyapunov function for power system frequency control with lossy transmis- sion lines, that is previously unknown in the literature.
Figure 1. We present the statistics of runtime for single SMT ver- ification (precision e tolerance error e e on the 6-D system (Appendix F.4). Each method is validated once per epoch. The red line represents the median. The mean verification time is averaged over one complete training. and the value is displayed in bold. Neural Network in FOSSII. 2.0 Edwards al. (2024) has 10-D hidden layer, and Augmented Neural Lyapunov Control (ANI.C Grande et al. (2023)) has two 10D hidden layers. Our formulas have a complexity (number of tokens) of 20 or fewer. For ours, the best candidate in each epoch verified in this motivation example.
2. Related Work
2.1. Learning-based Lyapunov Function Construction
integrate with formal verification tools like SMT solvers. This reduces parameter complexity and eliminates the need to verify complex neural network elements like nonlinear activations. Figure highlights the efficiency of verifying analytical expressions compared to neural networks.
The field of learning-basse baaed Lyapunov function construction is advancing rapidly. Chang et (2019) formulated Lya- punov condition violations as the objective, jointly learning neural Lyapunov function and linear controller to guar- antec stability for given system, with stability verified via SMT solvers. Zhou et al. (2022) extended this to unknown dynamics with a neural controller. Dai et al. (2021) and Wu et al. (2023) focused on discrete-time systems, using MIP solvers for stability verification, requiring piecewise linear approximations. Yang et al. (2024) applied a, 3-CROWN for scalable neural network verification, extending state feedback to output feedback control. However, scalability remains challenge: SMT solvers handle up to 30 neurons, MIP solvers scale to 200 neurons Dawson et al. (2023a), and or 3-CRRWN Yang et al. (2024) is limited to a three-layer architecture (16 neurons per layer).
In this work, we aim to address the following question:
Can neural networks effectively discover valid analytical Lyapunov functions directly from complex system dynamics?
To tackle this challenge, we introduce an end-to-end frame- work designed find analytical Lyapunov functions for non- linear dynamical systems given in analytical form. Building on the transformer's i ability to model complex dependen- cies Vaswani et al. (2017) and the success of deep symbolic regression methods Holt et al. (2023), our framework de- ploys a symbolic transformer Kamienny et al. (2022) for Lyapunov function discovery. The transformer's encoder captures system dynamics represented token sequences derived from the ordinary differential equations (ODEs), while the decoder generates candidate Lyapunov functions by modeling symbolic token distributions. Given the lack of high-quality (local) Lyapunov function datasets, particularly for high-dimensional systems, we propose a reinforcement learning (RL)-based approach to search for Lyapunov func- tions on per-system basis, instead of pre-training like Alfarano et al. (2024). We verify Lyapunov conditions by localized sampling in the neighborhoods of minimizers of the candidate expressions, which are most likely to have violations. The identified counterexamples are then incor- porated into the training set for further optimization. Our main contributions are summarized as follows:
In contrast to neural Lyapunov functions, Feng et al. (2024c) and Alfarano et al. (2024) derived analytical Lyapunov func- tions. Feng et al. (2024c) combined a neural network with the symbolic regression package PySR (Cranmer, 2023), which approximates the network to produce analytical Lya- punov functions, but the lack of interaction between system dynamics and symbolic regression component limits its po- tential. Alfarano al. (2024) pre-trained a transformer on backward- and forward weward-generated global Lyapunov function datasets, relying beam search for candidate generation. However, their method cannot adaptively refine the candi- date Lyapunov functions if the beam search fails on specific dynamics, and it requires a dataset that is expensive to gen- crate (e.g., thousands of CPU hours for a 5-D dynamics dataset) to achieve adequate generalization during inference.
2

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Furthermore, its emphasis on global stability limits its ap- plicability to real-world, nonpolynomial control systems, which typically only admit local stability. Consequently, RL-based approach that directly searches Lyapunov func- tions for a given system is indispensable.
where f D + I is a Lipschitz-continuous vector field, and D C TR" is set containing the origin that defines the state space. Each srt) E (D is state vector.
Definition 3.2 (Asymptotic stability). A system of (1) is stable at the origin if V > O, there exists 6 = 8(e) > 0 such that |z(0)]| < d = ||x(()]] < e, Vt 2 0. The origin is asymptotically stable if is stable and 5 can be chosen such that | = lim at) 0 (Khalil, 2002).
2.2. Symbolic Regression with Generative Model
Symbolic regression is supervised learning task, searching for an analytical function f: Ik= + R that fits yi E IR from input Tj E R" (Petersen et al., 2020). With suitable exten- sions, it can scale to multi-input-multi-uuu
1000
Definition 3.3 (Lic derivative). The Lie derivative of continuously differentiable scalar function V D - R along the trajectory of (1) is given by
mappings.
RL-baased Symbolic Regression. RL-based symbolic re- gression algorithms Petersen et al. (2020); Costa et al. (2020); Landajuela et al. (2021) employed generative mod- els, typically RNNs, to generate distributions of symbolic tokens representing mathematical operations and variables, from which analytical expressions are sampled. Rewards, evaluating the quality of the sampled expressions, are measured by fitness metries like RMSE. Due to the non- differentiable step of converting token sequences into sym- bolic equations, policy gradients are used to optimize the output distributions. Mundhenk et al. (2021) extended this approach by integrating Genetic Programming (GP) to refine generated expressions, improving the overall performance.
L,V(x)= = Σ OV dr,
OV = Σ i=1 dz, f((x).
(2)
dz. dt
Proposition 3.4 (Lyapunov functions for asymptotic stabil-
ity). Let x =0 =0be be an equilibrium point for (1) and D C R" be a domain containing the x 0. Let V : D - R be continuously differentiable function such that
V(0) = and V(x) > 0 in DV(0),
(3a)
Vv in DV(0),
(3b)
then the origin is asymptotically stable.
Definition 3.5 (Lyapunov risk). Consider a candidate Lya- punov function V for system f from Definition 3.1. For a dataset X = {=r,... ,NN] where x; E DD the Lyapunov risk of V Chang et al. (2019) over D is defined by
Pre-trained Symbolic Regression methods. These meth- ods are inspired by the success of transformers in Natural Language Processing (NLP) tasks. These algorithms con- tain two steps: 1) pre-train an encoder-decoder network to model DD) on curated datasets by cross-entropy loss, and 2) sample from this distribution during inference via beam search or Monte Carlo Tree Search (MCTS). Methods like
C(V)
1 N
maxx(0 (z)(x.)) + max(0, (xx))) (4)
4. Proposed Framework
Biggio et al. (2021); Kamienny et al. (2022); Bendinelli et al. (2023) rely on beam search without gradient refinement, of- ten yielding suboptimal results for -ittt---iistiittionn data. In contrast, Holt et al. (2023) integrates RL-based policy gradient optimization with end-to-end RMSE loss for both pre-training and inference, allowing gradient refinement for unseen datascts during inference. Further, Shojace et al. (2023); Kamienny et al. (2023) enhance the decoding pro- cess (expression generation) by incorporating MCTS with feedback, such fitting accuracy and equation complexity.
This section introduces our RL-based generative approach, which aims to find an analytical Lyapunov function for given dynamical system that certifies asymptotic stability following conditions in Proposition 3.4. Successfully ident-- fying such a function guarantees system stability.
The framework, visualized in Figure 2. consists of three components: 1) a symbolic transformer, parameterized as (), for candidate analytical Lyapunov functions generation, where Ç and 0 denote the parameters of en- coder and decoder, 2) numerical verifier employing the SHGO Endres et al. (2018) 8) global optimization algorithm for Lyapunov conditions' checking (Proposition 3.4) and counterexamples' feedback, and 3) risk-seeking policy gradient algorithm optimizing the transformer's parameters based candidate Lyapunov functions' rewards. To tackle the challenges posed by the exponentially growing search space of complex, nigh-dimensional systems, our frame- work integrates Genetic Programming as expert guidance to improve expression quality and training efficiency. We denote X C as the training set for reward calculation.
3. Preliminary
Our framework searches for analytical Lyapunov functions for autonomous nonlinear dynamical systems at an equilib- rium point. Without loss of generality, we assume the origin to be the equilibrium point.
Definition 3.1 (Dynamical systems). An n-dimensional autonomous nonlinear dynamical system is formulated as
dx
dt = f(x),x(0)) = Zo.
(1)
3

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Vo(a)
LyVg(x)
Localized Ball
Counterexamples
Reward calculation
Maximizer
(a) x2)
Minimizer
Binary tree Candidate Lyapunov function
Expression
Risk-soeking policy
Violations
gradient
Lyapunov verification
Sampled tokens
2
Analytical system model
System embedding
Encoder
Decoder
Encoder
Decoder
Encoder
Decoder
di fi (x) i22
E1
S
Z
E1
S
θ
E1
Z
0
/2 (x)
e2
E2
E2
Start
Hierarchical information
Provious
Previous
outputs
outputs
Sibling
Figure 2. Framework overview: The transformer takes emheddings of the dynamical system model as input and generates candidate Lyapunov functions autoregressive manner. Hierarchical information is deployed to enhance the model input. For example, when generating the last token 22, its parent is +. and its sibling is z1. The output is the pre pre-order traversal of the expression's 's binary tree. Candidates are verified using coobal-ppimizzatin based verification process, with counterexamples added to the training set for reward calculation. The transformer is updated via the risk-seeking policy gradient. The program terminates once valid expression is found.
4.1. Candidate Lyapunov Function Generation from Symbolic Transformer
Dynamies Tokenization. By Definition 3.1, the symbolic transformer models' input f(x) is composed by n ordinary differential equations da- fi()) for i 1, n. Each dt
Expression Representation. Inspired by the deep sym- bolic regression frameworks, we use symbolic transformer model as the backbone. The transformer takes a dynami- cal system f(x) as input and generates candidate analytical Lyapunov functions Va such that: (x(x) > 0 and T V, < 0, æ E D).. Following the expression representation rules in Lample& & Charton (2020), framework repre- sents symbolic transformer models' input and output as sequences of symbolic tokens. Each mathematical expres- sion can be converted into a symbolic expression tree, a binary tree where internal nodes are symbolic operators and terminal nodes (leaves in the tree) are variables or con- stants. Symbolic operators can be either unary (i.e., one child), such as sin, cos, or binary (i.e., children), such +, Furthermore, each symbolic expression tree can be represented as sequence of node values, either symbolic tokens or numerical coefficients, by its pre-order traversal (i.e., first visiting the parent, then traversing the left child and right child). In this way, each expression obtains a pre-order traversal representation, which can uniquely re- construct the original expression Petersen al. (2020). We denote Vo, as the th node value in the pre-order traversal of Va 's expression tree, and C. the symbolic library, e.g. { x, log, sin, 2z J. where Vo, sampled from.
analytical expression f((x) can be represented as sequence of symbolic tokens and numerical coefficients by the pre- order traversal of its expression tree. Concatenated the se- quences of pre-order traversal for all (xx), with sos (start token) and EOS (end token) separators, tokenized dynamics, which is fed into the encoder of sym- bolic transformer and encoded as latent vector F E R". The numerical coefficients are tokenized in two schemes: an integer represented as sequence of digits in base 10 (e.g. 123 tokenized as [1, 2, 3)). and a real number is represented scientific notation rounded
obtain the
significant
digits (e.g. 3. 14 is tokenized [3, 1,4, 4,400 10%). A detailed
example is illustrated in Figure 3, Appendix A, where
present the symbolic representations of the simple pendulum dynamics in sequences of pre-order traversal.
Candidate Expression Generation. The decoder generates candidate expressions Vo an auto-regressive manner. That is, each token Vo, in the pre-order traversal of V. is sampled from the symbolic library L. according to conditional distri- bution p( IVe..Ve.. Φ, f (z)). where Vona-1 is the first (i - 1) selected symbolic tokens in the pre-order traversal of Vg. This conditional dependence decoder, which emits probability distribution v over the symbolic library C., conditioned on the previously selected
be achieved by the

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
tokens and input dynamics. Since analytical expression in its pre-order traversal is inherently hierarchical, we also deploy the hierarchical tree state representation (Petersen et al., 2020; Holt et al., 2023). This method enhances the decoder input by concatenating the representations of the parent and sibling nodes with previously selected outputs and the dynamics, Once the token sampling process for Va is complete, we evaluate the function value at origin Vo and subtract it from Vg, ensuring the Lyapunov condition (0) 0. Through this process, we can sample batch of Q candidates V/ {V3 ~ p(Valo. f (x))) o be verified according to the Lyapunov conditions.
small value relative to l/P||a, are sampled check Lya- punov conditions, i.e. Vo()) > 0 and - (o)) > 0 for E DV(0). This localized sampling scheme effectively identifies violations within D. Additional random sampling covering approximately 30% of the total data and condition checking across the state space are performed to comple- ment this localized sampling to capture additional potential violations and provide a global assessment. Identified coun- terexamples Nxe are added to the training set X' for reward calculations. Once a candidate Lyарunov function passes this verification process and does not encounter any viola- tion in X, it indicates a numerically valid solution is found, pending the final formal verification. Appendix B details the verification implementation.
1-1- which will
4.2. Verification and Falsification Feedback
Leveraging the analytical nature of candidate Lyapunov functions, efficient methods for symbolic expressions, such as root finding Feng et al. (2024c), can be applied Lya- punov condition verification and counterexample generation. In this work, we propose a global-optimiiztion-beoe based nu- merical verification algorithm, using Simplicial Homology Global Optimization (SHGO) (Endres et al., 2018), that effectively checks the Lyapunov conditions around mini- mizers and feedback counterexamples into the training set X' for reward calculation. SHGO is a constrained global optimization algorithm with theoretical guarantees for con- vergence. To make the paper self-contained we present the guarantees in Proposition 4. This algorithm identifies the global minimizer over the state space from a set of local minimums, each of which is obtained from a convex sub- domain in the feasible search space. Taking advantage of the theoretical results, we employ the SHGO algorithm for counterexample detection in our verification process.
4.3. Disk-Seeking Policy Gradient
The empirical Lyapunov risk [(V)) in Equation (4) quan- tifies the violation degree of Lyapunov conditions over a given dataset. Following Petersen al. (2020); Bastiani et al. (2024), we apply the continuous mapping g(x) = It= and define proposed Lyapunov risk reward as:
1 R(V,) = 9 (C(V))) 1 C(V))'
(5)
where [(V²) is measured over training set X'. The contin- uous mapping g(z) bounds the reward value to [o, 1]. For candidate expressions that do not incorporate all state vari- ables or are analytically incomplete, we assign their rewards to be 0 to ensure they are effectively penalized.
Given the violation measure for each sampled Vo is non- differentiable with respect to the transformer parameters o, we employ the risk-seeking policy gradient to update & end-to-end. The objective of standard policy gradi- ent Williams (1992) is defined to maximize Jsa(() = (e)) (RVV)), the expectation of the reward function R(.)) for candidates' quality evaluation based on the current parameter ø. This objective is desirable for prob- lems aiming to optimize the average performance of the policy network. In our task, final performance depends on finding at least one valid Lyapunov function that meets the conditions in Proposition 3.4, rather than optimizing for av- erage performance. Consequently, standard policy gradient methods are inadequate due to the misalignment.
Proposition 4.1 (Convergence Gurantees of SHGO Endres et al. (2018).. For a given continuous objective function f that is adequately sampled by sampling set of size N.. If the size of the minimizer pool M extracted from the directed simplex (a convex polyhedron) H. is MM] Then any further increase of the sampling size N₂ will not increase IM).
This result shows that if the initial points are adequately sampled, that is, the union of identified locally sub-convex domains initiated from starting points covers the feasible search space, then the minimizer pool M, which contains all local minimum extracted from the directed simplexes, will contain the global minimizer of the feasible search space. Notably, the required sampling size N, can be unbounded.
In our framework, we adopt risk-seeking policy gradient Petersen et al. (2020) that only focuses on maximizing best- case performance. Let R,, (Ø) as the 1- a quantile of the distribution of rewards of sampled candidates under the current policy d. The learning objective of risk-secking policy gradient, parameterized by a, is formulated as:
During verification, for a candidate Vas SHGO is first ap- plied to identify minimizers z and 2; of V. and its negated Lie derivative Vy in the state space D. These minimiz- ers highlight the regions where Lyapunov conditions are most likely to fail. Next, data points x from neighborhoods around these minimizers, (x]) and B, (x²)) where r is
Jna a) = Evo-m(Volo.ra [ R(V)) 2 Ro)
(6)
5

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
Algorithm 1 Trriniin Framework for Analytical Lyapunov Function Discovery via Reinforcement Learning
Lyapunov risk reward as the measure of fitness, and obtain batch of refined expressions. We select an 'elite set' of the refined expressions Vap = {V sp GP GP(V)) them as target expressions, and optimize the transformer model in a supervised learning manner, maximizing the probability that the generated token matches the reference tokens from Vgp with the following expert guidance loss:
Input: Dynamics f()), state space D, quantile a, symbolic library L, batch size Q, max complexity k, radius r.
regard
Output: Valid Lyapunov function V for system f(x).
1: Initialize the conditional generator with parameters o,
2: Randomly sample training datapoints X'
(ø1,... ".I,..,,N' , where æ; E DD
1 0 ((Vop) =o2 KG
Ni - log (p(V So
3: while no valid candidates found do
-
RVV.
18 +. (7)
4: Vo + [V3 ~ p(Vole, f()))
J-1
5: Vop + Genetic Programming(Vs).
where G' is the number of expressions in V sp- and k; is the complexity (number of symbolic tokens in the pre-order traversal) of Vi gp' Each expression V is weighted by its Lyapunov risk reward in C. Algorithm 1 summarizes the training process, with more details in Appendix D. The GP solutions explore the characteristies of Lyapunov functions that have not been captured by the transformer yet and effectively guide the transformer.
6: V+ U.P.
7: V., P", + Cerification(V,, r), [Verify candidates
V. return the valid Lyapunov function V (if any) and counterexamples No.. Details in App. B. 8: if V' is not empty then
9: Return V"
10:
end if
11:
R t (R(V)) EV)
v)
Remark 4.2 (Exploraton--explottttoo trade-off). our framework, exploration arises during candidate expression generation and through GP's evolutionary operations, while exploitation is driven by the risk-secking policy gradient and the expert-guided loss.
12:
(a( (1-a
R.
13:
@ + -- VaJria( a). (risk-seeking policy gradi- ent update. Equation (6) }
14:
@ - 2- V (V fexpert guidance loss based on the genetic programming refined Lyapunov functions
Vap to update policy. Equation (7).}
15:
Concatenate counterexamples N to dataset X.
5. Experiment
16: end while
We validate the proposed algorithm across a variety of non- linear dynamics by finding their local Lyapunov functions the equilibrium point to verify their stability, where the sys- tems are autonomous (or closed-loop systems with known feedback control laws). We use dReal Gao et al. (2013) SMT solver for final verification of found Lyapunov func- tions, with a numerical tolerance error E = e and pre- cision e e- over the state space, i.c. V(x) 8 and (x) < -ō over BiBB.() The excluded ball (0) is to avoid numerical issues, which is common practice for SMT-based formal verification (Chang et al., 2019). Global stability can be determined through further expert analysis; for instance, if the Lie derivative is a negation of SOS, it is sufficient to establish global stability.
This objective aims to optimize only the rewards of high- quality candidates from the top 1 - a quantile.
4.4. Automated Expert Guidance
While the risk-seeking policy gradient algorithm effectively optimizes the model, training efficiency can be enhanced by off-the-shelf tools that further explore the function space based on the transformers. Inspired by Mundhenk et al. (2021), we incorporate a Genetic Programming (GP) compo- nent (DEAP Fortin al. (2012)) into the training paradigm.
The GP algorithm starts with batch of initial populations (expression trees) and iteratively refines these populations through evolutionary operations: mutation, selection, and crossover, with a pre-defined metric to evaluate the fitness of populations to the task (e.g., MSE for symbolic regression task). Within our framework, we feed the latest batch of generated candidates (Vg ~ pVVlo. f()))) from the decoder into the GP module as the starting population', re- fine these expressions through evolutionary operations with
5.1. Experimental Setting
Test Dynamics. We categorize our collected test dynamical systems into two kinds: 1). Polynomial Systems, and 2). Non- non-polynomial Systems, where three polynomial systems are adopted from Alfarano et al. (2024) (Appendices F.2 & F33), and others come from real-world examples. Detailed information about systems is summarized in Appendices F and G. The dimension of these test systems ranges up to 10.
'Without good initial population from the transformer, GP algorithms itself face significant challenges in directly finding valid Lyapunov functions for high-dimensional systems due the exponentially growing search space. Appendix H.3.
Implementation Details of Framework. Detailed expla- nation for all components in our framework is presented in Appendices A (Transformer), B (Global-optimization based
6

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
Table 1. Performance and time consumption of our method dynamics. 'App.' refers Appendix.
Dynamics
Runtime
Ver."
Found Lyapunov Functions
Stab
Suce
2-D Polynomial Sys (Ape. F2)
68s
2 ms
V 9zî + 2x²
La.s.
100
2-D Van Der Pol (App.
126s
1 ms
V = x + x2
La.s.
100
2-D Simple Pendulum (App. G.I)
288s
1 ms
V = =2(1 - cos((1)) + +
La.s.
100
3-D Polynomial Sys (App. E3)
112s
ms
V 9x + x + +
La...
100
3-D Trig Dynamies (App. (.2)
157s
ms
V cos(c2)) x11) x + sin(xs)?
La.s.
100
4-D Lossy Power Sys (App. G.5)
3632s
621s
V = w/ + w/ + (w2 - sin( + sin
Lass.
100
6-D Polynomial Sys (App. F4)
1667s
7 ms
V G
la..s
100
6-D Quadrotor (App. G4)
3218s
ms
V = 21=1
g.a.s.
80
6-D Lossless Power Sys (App G.3)
18094s
2 ms
V =
Li=1 i=
w
l.a.s.
60
3
0.5
Σ
9-D Synthetic Sys (App. G.6)
27047s
6.6 s
V
Σ x²
+ sinn(2)) (27)
La.s.
60
+
2
cos((1)) +
10-D Polynomial Sys (App. FG)
64223s
2 ms
V
E 10 11=1 x
l.a.s.
60
a. 'Ver.' presents the time consumption for the final verification of the found Lyupuzov functions, All found Lyspunoe functions passed SMT solver's verification.
1. 'Stab' meuns siability this column, 'gass' represents globully asymputtiaally stable, and 'lass' represents locally asmpputiiaaly stable.
'Suce denotes the successful finding a valid Lyspunov function of random seeds.
Numerical Verification), C (Risk-secking Policy Gradient), and D (Genetic Programming). The symbolic library C, is defined as {+, x, sin, cos, x in all tests.
5.2. Performance Analysis
Table summarizes the runtime, success rate, and discov- ered Lyapunov functions for selection of tested nonlin- ear systems, ranging from 2-D to 10-D, demonstrating the robustness and scalability of our framework. As dimen- sionality increases, runtime grows with the exponentially expanding search space, reward calculations, SHGO opti- mization, and genetic programming.
Baseline Algorithms. We compare our proposed frame- work against four baseline algorithms for continuous non- linear dynamics. Neural methods: 1) Augmented Neural Lyapunov Control (ANLC) Grande et al. (2023), and 2) FOSSIL 2.0 Edwards al. (2024). Analytical methods: 3) the transformer-based global Lyapunov search of Alfarano et al. (2024), and 4) SOS methods via SOSTOOLS (Matlab) (Papachristodoollou et al., 2013). We use the formulation of Pupachristodoulou& & Prajna (2005a) on polynomial systems, and apply the recasting technique of Papachristodoulou & Prajna (2005b) convert non-pols polyyooollia dynamics into rational form they can also be handled by SOS.
Unlike existing methods that produce neural Lyapunov fune- tions, our framework yields interpretable analytical candi- dates. For example, it correctly identifies the energy func- tion as a valid Lyapunov function for the simple pendulum. Likewise, for the 3-bus power system (Appendix G.3), discovers the commonly used energy- storage function for incremental passive systems Weitenberg et al. (2018).
Both ANLC and FOSSIL 2.0 train a neural Lyapunov func- tion with the empirical Lyapunov risk loss and employ a counter-example guided inductive synthesis (CEGIS) loop for better generalization over the state space D. Alfarano et al. (2024) pre-trains a transformer on global Lyapunov function datasets with beam search for candidate generation. SOS methods formulate Lyapunov functions as the feasi- ble solutions of some semi-definite programming tasks and solve these tasks via convex optimization tools. Details of each baseline can be found in Appendix E.
Analytical Lyapunov functions can potentially bypass the need for formal verification. In the 3-D Trig system (Ap- pendix G.2), over the state space D = {(21 Z2, x3) E TR3 the positive definiteness of the
|r:| < 1.5, V E identified Lyapunov function is evident from its formulation.
Moreover, the Lie derivative LyV = -2x² - I'3 sin sin(2xg) is directly identifiable as non-positive in D, since x sin(æ) > for all (- By the invariance principle Khalil (2002), the discovered function certifies the asymptotic sta- bility of the origin in state space D. When direct identi- fication is non-trivial SMT solvers can efficiently verify
7

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
Table 2. Training time and success rate comparison between ours and the neural baselines. The Succ % is the successful rate of finding a valid Lyapunov function within 5 hours out of 5 random seeds. Runtime is the average training time for a successful trial.
2-D Dynamics
3-D Dynamies
6-D Dynamics
8-D Dynamics (App. F.5)
Frameworks
Succ %
Runtime
Succ %
Runtime
Suce %
Runtime
Suce %
Runtime
Ours
100
165s
100
124s
80
6290s
80
14358s
ANLC
53.3
1.409s
46.7
63.91s
0
FOSSIL 2.0
80
7.708s
66.7
221s
0
0
Lyapunov conditions given analytical formulations' simplic-
tlenecks for both baselines' counter-example feedback paradigm. Even networks with fewer than 15 neurons per layer may require hours to finish formal verification. In con- trast, our method remains robust up to 6-D and 8-D systems, as its simpler analytical formulations allow numerical verifi- cation efficiently identify violation regions. Moreover, for final verification, our candidates can pass the SMT solver in milliseconds, thanks to term cancellations and algebraic restructuring made possible by their analytical form.
ity.
5.3. Newly Discovered Lyapunov Function
Despite decades of effort in the control community to iden- tify Lyapunov functions, certain stable dynamics still lack a valid Lyapunov function to directly certify their stability. One example is the lossy frequeney dynamics in power sys- tems Chiang (1989); Cui & Zhang (2022), where "lossy" refers to the consideration of energy losses in transmission lines. For simplicity, we focus on 2-b (4-D) lossy system with the equilibrium point set at the origin, with detailed descriptions provided in Appendix G.5. Using the proposed method, we successfully discover two local Lyapunov func- tions valid in the considered region D = {(61, , 82, w1, W2) E R² [611 s 0.75 and |ws| S for i = 1,2]:
A Analytical Lyapunov Function Baselines. We evaluated the pre-trained model of Alfarano et al. (2024) on our bench- marks. Trained solely on globally stable systems with fewer than six states, it produced valid Lyapunov functions only for the low-dimensional examples in Appendices F.1, F.2, E.3, & G.I, which have global stability guarantees, and failed on every benchmark that is only locally stable
2
Table 3 compares the training and solving times of our method and the SOS approach on polynomial systems from Appendix F. For low-dimensional systems, SOS finds valid Lyapunov functions more efficiently, with shorter run- times. However, it fails to certify stability for the higher- dimensional systems in Appendices F.4, F.5, and F.6, which are only locally stable. While SOS methods can scale to 10-dimensional or higher systems for global stability ver- ification, local stability requires additional constraints on the Lie derivative, which significantly impact scalability. For example, in the 6D polynomial system from Appendix F.4, verifying local stability with a degree-2 polynomial candidate introduces hundreds of symbolic terms in the Lie derivative constraint, rendering the problem intractable for SOSTOOLS. More details are provided in Appendix E.4.
(1 82, w1, w2)
Σ
(in( - sin(61) + w2) w/ +
V2 62, w1, W2)
Σ
2
2 w + (sin( - sin(81) - wi)
Both functions are formally verified by the SMT solver within the defined state space. To the best of our knowledge, these are the first analytical Lyapunov functions used to certify the local stability of 2-bus lossy power system.
5.4. Comparisons with Baselines
Neural Lyapunov Function Baselines. Table 2 compares our success rate and training runtime with two neural Lya- punoy function baselines across various test dynamics. For low-dimensional systems, both baselines achieve notably shorter overall training runtime. However, to achieve ef- ficient verification, these methods rely on relatively small neural networks, which fail converge on more challenging tasks (e.g., the simple pendulum and 3-D Trig dynamics in- volving trigonometric terms). Consequently, both baselines exhibit lower overall success rates than our approach.
To use SOS on non-polynomial systems without relax- ing the dynamics, we follow the recasting scheme of Pa- pachristodoulou & Prajna (2005b), replacing each non- polynomial term with auxiliary variables linked by extra (in)equalities. We applied this procedure to the simple pen- dulum and the 3-D trig system (Appendices G.1 & G.2). SOS did recover Lyapunov certificates in both cases, but at a much higher computational cost than our framework for the 3-D trig example. Three main limitations are revealed during implementation: (i) recasting demands substantial
As dimensionality grows, the complexity of the Lie deriva- tive of neural Lyapunov function increases significantly following Definition 3.3, creating severe verification bot-
8

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
Table Training/solving time of ours and sum-of-squares (SOS) polynomial systems, averaged over successful trials. The stable
regions (local or global) of the considered dynamics are indicated in smaller font.
2-D Systems
3-D Systems
6-D System (App. F.4, local)
8-D System (App. B55, local.)
10-D System (App. F.6, local.)
Test Systems
(App. F.1, local; App. global.)
(App. E3, global))
Ours
97s
108s
1667s
14358s
64223s
Sun-of-squares
0.765s
1.503s
domain expertise and hand-crafted constraints; (ii) added variables and relations greatly increase the computational burden; and (iii) the formulation certifies only stability, not asymptotic stability. For the locally stable 3-D trig system, eliminating the trigonometric terms introduces four aux- iliary variables and two equality constraints, pushing the solve time beyond one hour, versus 157s for our method. These limitations undermine SOS's practicality for high- dimensional, non-polynomial dynamics.
refinement + expert guidance (ours). While the transformer only can achieve a 100% success rate, it requires triple the training time compared to 3) and 4). GP only fails to con- verge due to its limited understanding of system dynamics. In contrast, GP refinement and expert guidance learning efficiently accelerate transformer parameters' update and enable faster Lyapunov function discovery (Appendix H.3).
6. Conclusion
This work introduces an end-to-end framework for discov- cring analytical Lyapunov functions for nonlinear dynam- ical systems. A symbolic transformer, trained with risk- seeking policy gradient and augmented by genetic program- ming, proposes candidate expressions; the SHGO global op- timizer rapidly verifies them and generates counterexamples during training; and an SMT solver certifies the final Lya- punov candidates. The framework scales to 10-dimensional dynamics and has discovered previously unknown local Lya- punov functions for lossy power system dynamics. It can be extended to other certificate-functin discoveries, such control barrier functions for safety certificates.
Table Training/solving time of ours and sum-of-square (SOS) on two non-polynomial systems, averaged over successful trials.
Frameworks
App. G.1
App. G.2
Ours
288s
157s
Sum-of-squares
19.58s
6163s
5.5. Ablation Studies
Risk Sesk-Seking Quantile a. We compare the performance on the 3-D Trig dynamics without GP refinement under 0.1,1 0.5, 1. With a 0.1, the framework achieves steady convergence, lower variance, and the highest success rate of 66. 66.67%. outperforming
Several promising future research directions emerge from this work. One is efficiently incorporating physical con- stants, such as the gravitational constant and object mass, which could significantly improve generalization-ttouu directly introducing them as variables may add unnecessary complexity. Another is extending the framework to discover control Lyapunov functions, which essential for design- ing stabilizing feedback laws. Theoretical analysis of the framework, including proofs of completeness and conver- gence, would also be valuable. Finally, since constructing large-scale datasets for local Lyapunov functions remains challenging, using our approach to refine pre-trained mod- els (Alfarano et al., 2024) offers an exciting opportunity to further advance Lyapunov function discovery.
0.5 (33.33%) and the vanilla policy gradient a = 1 (0%). These findings confirm the importance of the risk-seeking strategy (Appendix H.I).
Verification Comparison. We compare three verification methods-1l root-finding (Feng et al.. 2024c), 2) the dreal SMT solver, and 3) random caaplling-ggiine our approach on the 6-D polynomial system. SHGO-based counterex- ample feedback is an adversarial reward, allowing fast re- finement and stronger final guarantees but risking training instability. In contrast, random sampling yields smoother rewards and more stable training, yet it identifies violations less efficiently. We thus blend both methods to balance final performance and training stability. Meanwhile, the SMT solver and root-finding produce mostly mild violations, of- fering limited optimization guidance (Appendix H.2).
A Acknowledgment
The authors thank Dr. Huan Zhang at the University of Illinois Chammaaa for helpful discussions. J. Feng is supported by the UC-National Laboratory In-Residence Graduate Fellowship L24GF7923. Y. Shi acknowledges sup- DOE grant DE-SC0025495.
Expert Guidance. We evaluate four settings on the 6-D polynomial system to assess the impact of GP refinement and expert guidance: 1) transformer only, 2) GP only, 3) transformer + GP refinement, and 4) transformer GP
port from NSF

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
and the Schmidt Sciences A12050 Early Career Fellowship.
Chang, Y.-C., Roohi, N., and Gao, S. Neural Yyapunov con- trol. Advances in neural information processing systems, 32, 2019.
Impact Statement
This paper advances learning-base based Lyapunov function dis- covery in control theory by proposing a systematic approach to finding Lyapunov functions for any dynamical system. Once valid Lyapunov function found, guarantees sys- tem stability in the state space and offers insights for control engineering. Our method can be applied in energy systems, robotics, transportation, and other control systems. While these improvements have the potential to enhance safety and efficiency, we do not foresee specific negative ethical or societal consequences arising directly from our approach.
Chiang, H.-D. Study of the existence of energy functions for power systems with losses. IEEE Transactions on Circuits and Systems, 36(L1):1433-1229 1989. doi: 10. 1109/31.4198.
Costa, A., Dangovski, R., Dugan, O., Kim, S., Goyal, P., Soljačić, M., and Jacobson, J. Fast neural mod- els for symbolic regression at scale. arXiv preprint arXiv:2007.10784 2020.
Cranmer, M. Interpretable Machine Learning for Sci- ence with PySR and SymbolicRegression.jj. May 2023. URL
References
org/ borg/abat absr//abs/333... 01582. arXiv:2305.01582 [astro-ph, physics:physics)
Ahmadi, A. A. and Majumdar, A. Some applications of polynomial optimization operations research and real- time decision making. Optimization Letters, pp. 709-729, 2016.
Cui, W. and Zhang, B. Lyapunov-regullaiized reinforce- ment learning for power system transient stability. IEEE Control Systems Letters, 6:974-979, 2022. doi: 10.1 109/
ACCSYS2022.308880
Ahmadi, A. A., Krstic, M., and Parrilo, P. A. A globally asymptotically stable polynomial vector field with no polynomial Σуарunov function. In 2011 50th IEEE Con- ference Decision and Control and European Control Conference, pp. 7579-7580, 2011. doi: 10.1 2011.6161499.
Cui, W., Jiang. Y., and Zhang, B. Reinforcement learning for optimal primary frequency control: A lyapunov approach. IEEE Transactions Power Systems, 1672))6676-166 2023a. doi: 2..11099/PWRS..222217776.6
Cui, W., Jiang, Y., Zhang, B., and Shi, Y. Structured neural- PI control with end-to-end stability and output tracking guarantees. Thirty-seventh Conference on Neural In- formation Processing Systems, 2023b.
Alfarano, A., Charton, F., and Hayat, A. Global іyapunov functions: a long-standing open problem in mathematics, with symbolic transformers. In The Thirth-eighth Annual Conference on Neural Information Processing Systems, nww./oruum forum? 2024. URL https::///........... : http::///pe..ee.eeeeeoeo net/ id=kOMrm4 4m.3.
Dai, H. and Permenter, F. Convex synthesis and verification of control-lyapunov and barrier functions with input con- straints. In 2023 American Control Conference (ACC), pp. 4116-4123. IEBE, 2023.
Bastiani, Z., Kirby, M., Hochhalter, J., and Zhe, S. Complexity-aware deep symbolic regression with ro- bust risk-seeking policy gradients, 2024. URL https: / //openreview. net wwtp:/ert/fofmrmmeemm i w
Dai, H., Landry, B., Yang, L., Pavone, M., and Tedrake, R. Lyapunov-stable neural- nettol- control. In Proceedings of Robotics: Science and Systems, 2021.
Bendinelli, T., Biggio, L., and Kamienny, P.-A. Controllable neural symbolic regression. In International Conference on Machine Learning, pp. 2063-2077. PMLR, 2023.
Dawson, C., Gao, S., and Fan, C. Safe control with learned certificates: A survey of neural Iyapunov, barrier, and contraction methods for robotics and control. IEEE Transactions on Robotics, 39(3):1749-1767, 2023a. doi:
Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., and Paras- candolo, G. Neural symbolic regression that scales. In International Conference on Machine Learning, pp. 936 945. PMLR, 2021.
10.
Dawson, C., Gao, S., and Fan, C. Safe control with learned certificates: A survey of neural lyapunov, barrier, and contraction methods for robotics and control. IEEE Transactions on Robotics, 39(3):1749-1767, 2023b. doi: 10.1109/TR..2022.33355
Bouabdallah, S., Murrieri, P., and Siegwart, R. Design and control of an indoor micro quadrotor. In IEEE Interna- tional Conference on Robotics and Automation, 2004. Proceedings. ICRA '04. 2004, volume 5, pp. 4393-4398 Vol.5, 2004. doi: 10.1 00.1100/RHHT.2224111022220
Edwards, A., Peruffo, A., and Abate, A. Fossil 2.0: Formal certificate synthesis for the verification and control of
10

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
dynamical models. In Proceedings of the 27th ACM In- ternational Conference Hybrid Systems: Computation and Control, pp. 1-10, 2024.
Holt, S., Qian, Z., and van der Schaar, M. Deep generative symbolic regression. In The Eleventh In- ternational Conference Learning Representations, 2023. URL / net. net/forun? id=o 7 KoREMA1bR.
Endres, S. C.,, Sandrock, C., and Focke, W. W. simplicial homology algorithm for lipschitz optimisation. Journal of Global Optimization, 72(2):181-217, 2018. doi: 10. 0777/1100888-11--64455-
Kamienny, P.-.., d' l'Ascoli, S., Lample, G., and Charton, F. End-to-end symbolic regression with transformers. In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K. (eds.), Advances in Neural Information Processing Systems, 2022. URL https: openreview .et// forum?id=GoOuIIDHGE.
Feng, J., Cui, W., Cortés, J., and Shi, Y. Bridging transient and steady-state performance in voltage control: A re- inforcement learning approach with safe gradient flow. IEEE Control Systems Letters, 7:2845-2850, 2023a. doi: 10.1 WWW..CSSYS.2223338842488
Kamienny, P.-.., Lample, G., Lamprier, S., and Virgolin, M. Deep generative symbolic regression with monte-carlo- tree-search. In International Conference Machine Learning, pp. 15655-15668, PMLR, 2023.
Feng, J., Shi, Y., Qu, G., Low, S. H., Anandkumar, A., and Wierman, A. Stability constrained reinforcement learning for decentralized real-time voltage control. IEEE Transactions on Control of Network Systems, pp. 1-12, 2023b. doi: (1/////NSS.023.3332..
Khalil, H. K. Nonlinear systems. Prentice Hall, 3nd edition, 2002.
Feng. J., Cui, W., Cortés, J., and Shi, Y. Online event- triggered switching for frequency control in power grids with variable inertia. IEEE Transactions on Power Sys- tems, pp. 1-13, 2024a. doi: 20.1099.PWRS.2024..
Koza, R. Genetic Programming: On the Programming Computers by Means of Natural Selection, volume 1, MIT press, 1992.
3523262..
Lample, G. and Charton, F. Deep learning for symbolic mathematics. In International Conference on Learning Representations, 2020. URL https: openrevie...e.eeeeeeee net/ Landajuela, M., Petersen, B. K., Kim, S., Santiago, C. P., Glatt, R., Mundhenk, N., Pettit, J. F., and Faissol, D. Discovering symbolic policies with deep reinforcement learning. In International Conference on Machine Learn- ing, pp. 5979-5989. PMLR, 2021.
Feng. J., Muralidharan, M., Henriquez-Auba, R., Hidalgo- Gonzalez, P., and Shi, Stability-constranned learning for frequency regulation in power grids with variable inertia. IEEE Control Systems Letters, 99994-99,, 2024b. doi: 10.1109/LCSSS 1..11999LSSSS.2224444000063
Feng. J.. Zou, H., and Shi. Y. Combining neural networks and symbolic regression for analytical łyapunov function discovery. In ICML 2024 Workshop: Foundations of Reinforcement Learning and Control Connections and Perspectives, 2024c. URL https: //ppn///eepe...e.....mmmm futt/forum?iie (inerummiiiXjjj. 77 78WY9T4.
Madry, A., Makelov, A., Schmidt, L.., Tsipras, D., and Vladu, A. Towards deep learning models resistant to adversarial attacks. CoRR, abs/1706.06083, 2017. 2017. URL www.orgg//b///...oe//e/oo 116.. 06083.
Fortin, F.-A., De Rainville, F.-M., Gardner, M.-A. G., Parizeau, M., and Gagné, C. Deap: Evolutionary al- gorithms made easy. The Journal of Machine Learning Research, 1311:221755 2012.
McGough, J. S., Christianson, A. W., and Hoover, R. Symbolic computation of Iyаpunov functions using evo- lutionary algorithms. In Proceedings of the 12th IASTED international conference, volume 15, pp. 508 515, 2010.
Gao, S., Kong. S., and Clarke, E. M. dreal: An smt solver for nonlinear theories over the reals. In Bonacina, M. P. (ed.), Automated Deduction CADE-24, 24, pp. 208-2111 214, Berlin, Heidelherg, 2013. Springer Berlin Heidelberg.
Mundhenk, T. N., Landajuela, M., Glatt, R., Santiago, C. P., faissol, D., and Petersen, B. K. Symbolic regression via deep reinforcement learning enhanced genetic pro- gramming seeding. In Advances in Neural Information Processing Systems, 2021.
Grande, D., Peruffo, A., Anderlini, E., and Salavasidis, G. Augmented Neural Lyapunov Control. IEEE Access, 11:67979-67986, 2023. doi: 10.1109/ACCESS.2000 2003. 3291349.
Papachristodoulou, A. and Prajna, S. A tutorial on sum of squares techniques for systems analysis. In Proceed- ings of the 2005, American Control Conference, 2005., pp. 2686-2700 vol. 4, 2005a. doi: 10.1 109/ACC.0055 1470374.
Grüne, L. Computing Iyapunov functions using deep neural networks. Journal of Computational Dynamics, 8, 2019. doi: ...3334/jcd.2021000

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Papachristodoulou, A. and Prajna, S. Analysis of Non- polynomial Systems Using the Sum of Squares Decom- position, pp. 23-43, Springer, Aug 2005b. ISBN 978-3- 540-23948-2. doi: 10. 10.1007/10997703.22
https////ppeeeeeeeee/eeeeee/eeee / 1d=3xFMM9JJRR..
2024. URL https:
forrrrr forum?
Zhou, R., Quartz, T., Sterek, H. D., and Liu, J. Neural yyapunov control of unknown nonlinear systems with stability guarantees. In Advances in Neural Information Processing Systems, 2022.
Papachristodoulou, A., Anderson, J., Valmorbida, G., Prajna, S., Seiler, P., and Parrilo, P. A. SOSTOOLS: Sum of squares optimization toolbox for MAT- LAB. 2013.
0ttss/111..11 from
Available http:/ httt:://ww..r cds. caaa.aalll caaltechh.dd//eeeee eddu//stoosssssss
Petersen, B. K., Larma, M. L., Mundhenk, T. N., Santi- ago, C. P., Kim, K., and Kim, T. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2020.
Shojace, P., Meidani, K., Barati Farimani, A., and Reddy, C. Transformer- baseed planning for symbolic regression. Advances in Neural Information Processing Systems, 36: 45907-45919, 2023.
Sontag, E. A universal'construction construction of artstein's theorem on nonlinear stabilization. Systems & control letters, (2):117-123, 1989.
Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, N., Kaiser, u., and Polosukhin, Attention is all you need. In Guyon, Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL
oc cc/paper_fiiee ppapee__fleesppppp////////Le 35502445577611111033111024252223 pdf.
Wang, Z., Andersson, S. B., and Tron, R. Lyapunov neural network with region of attraction search, 2024.
Weitenberg, E., De Persis, C., and Monshizadeh, N. Expo- nential convergence under distributed averaging integral frequency control. Automatica, 98:103-113, 2018
Williams, R. J. Simple statistical following algo- rithms for connectionist reinforcement learning. Machine learning, 8(3-4):229-25 256,2,,22.622.. 1992.
Wu, J.. Clark, A., Kantaros, Y., and Vorobeychik, Y. Neural Ιyapunov control for discrete-time systems. In Thirty- seventh Conference on Neural Information Processing Systems, 2023.
Yang. L.. Dai, H., Shi, Z., Hsich, C.-J,, Tedrake, R., and Zhang, H. Lyapunov-stable neural control for state and output feedback: A novel formulation. In Forty- first International Conference on Machine Learning,
12

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Contents of Supplementary Materials
Framework Details:
Appendix A: Symbolic Transformer Model
Appendix B: Global-optimizzaion---sedd Numerical Verification
Appendix C: Risk-seeking Policy Gradient
Appendix D: Genetic Programming
Baseline Descriptions:
Appendix E.1: Augmented Neural Lyapunov Control (ANLC)
Appendix E.2: FOSSIL 2.0
Appendix E.3: Global Lyapunov Function Discovery by Pre-trained Transformer
Appendix E.4:
Sum-of-Squares Methods
Tested Dynamical Systems:
Appendix F: Polynomial Nonlinear Dynamical System
Appendix G: Non-polynomial nonlinear Dynamical System
Ablation Studies:
Appendix H.1 Ablation Study 1 Risk-secking Policy Gradient
Appendix H.2 Ablation Study II Global-opptimization-aasss Numerical Verification Process
Appendix H.3 Ablation Study III Supervised Learning with Expert Guidance Loss
13

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
Simple Pendulum
Dynamies Formulation
Expression Trees
Sequences of Traversal
11
222
[xc]
sin (1)
mf² 22
Tree
txx, 9
sim,x,,,, m1 22
sin
22
Tree 2
Figure visualize dynamics tokenization process of simple pendulum system. Each analytical formula in the ODE representa- tion of input dynamics is first converted into expression tree and then represented by the pre-order traversal in symbolic tokens and constant coefficients.
Computation Resources: All experiments in this work were performed on a workstation with an AMD Ryzen Threadripper 2920X 12-Core Processor, an Nvidia RTX 2080Ti GPU 11 GB, and an Nvidia RTX 2080 GPU GB.
A. Symbolic Transformer Model
We outline the details of our symbolic transformer model, a conditional generator for analytical candidate expression generation, comprising two components: 1) an encoder, and 2) a decoder. Ø = {c denotes the transformer parameters, where Ç and 0 denote the parameters of encoder and decoder respectively.
A.1. Encoder Structure
Our framework employs the vanilla transformer encoder from Vaswani et al. (2017) to encode two types of input information: 1) the input system dynamics f(x), encoded into a latent vector F E IR, where E R, and 2) the hierarchical tree state representation (Petersen et al., 2020) of the selected tokens, encoded as latent vector W E R where k E R. The resulting representations are concatenated as Z. Both inputs are expressed as se- quences of symbolic tokens and numerical coefficients when fed into the encoder. For damped pendulum example 0.5 kg. 1 1 m, 9 9.81, and b 0.1, the dynamics , 9, 8, 1, , x1, x, - 2, 0, 10-1, 22, BOS]. In this work, we set the embedding dimension to 128, attention head to 2. and applied a 2-layer transformer encoder for system dynamics fir encoding and a 3-layer transformer encoder for hierarchical tree state representation encoding.
in Figure 3, suppose we have
be tokenized as:
[SOS, 22, EOS, SOS, x,
A.2. Decoder Structure
The decoder of the symbolic transformer model also uses the vanilla transformer decoder, with an additional linear layer to output the token probability ub over the symbolic library L. for token selection. Candidate I Lyapunov functions V,, are sampled sequences of symbolic tokens in pre-order traversal. Each symbolic token Va. is sampled autogressively from conditional distribution p( VV. VV Φ, (=)). Upon token sampling is complete for Vg. subtract (o from the candidate expression to enforce the Lyapunov condition (0) = 0. In each epach, a batch candidate I Lyapunov functions (V/ p(Vol, { is sampled candidates, which are verified by cotiaallotimizaiin-baaaaee numerical verification. In this work, we set the embedding dimension to 128, attention head to 2, and applied a 6-layer transformer decoder for the candidate expression generation. In each epoch, we sample Q = 500 expressions as candidates.
B. Global-optimizziion--ears
Numerical Verification
For given dynamics f(x), suppose V. is an invalid analytical candidate Lyapunov function. According to Lyapunov conditions defined in Proposition 3.4, for w1 2z E DD where w;, 2² are the global minimizers of V. and -T,, Vo in the state space D, the following two inequalities hold: (z(x)) < 0 and Vz 2 0. This implies that if V. is invalid, the neighborhoods of z and zū are highly likely to capture significant violations. Based on this observation, we propose a global-optimization based numerical verification. This verification identifies minimizers xi and x by Simplicial Homology
14

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Localized Ball
4.5
1
3.0
1.5
-1
2
0.0
3
1.0
1.0
0.5
0.5
-0.5
0.0
-1.0
0.0
0.5 +
-0.5
0.5 +
0.0
0.0
X1
0.5
-1.0
X1
0.5
1.0
1.0
1.0
V. (x1 + x2) xx2
LVV
2x, (x, + x,)
(x - xx 2 1)) +x2)) + +1)
Figure 4. This plot visualizes proposed verification process sampled candidate. V. (x1 x2)² 2 22 is sampled candidate during the training for Van Der Pol Oscillator. Using Simplicial Homology Global Optimization, we first identify the minimizer of the Lyapunov function and the maximizer of the Lie Derivative, the black dots in each graph. Next, data points sampled in the neighborhoods of the two points, the regions in red circles. For sampled data points that violate the Lyapunov conditions, we feed them into the training set X.
Global Optimization (SHGO), verifies Lyapunov conditions on localized samples in neighborhoods (=i) and (x2), and feeds counterexamples back into the training set X. We detail the sampling and manditin-chhcking procedures in Algorithm 2. Figure 4 illustrates this verification process sampled candidate, Vo (æ1 + z²2² I2, for the Van der Pol Oscillator. In the implementation, we initiate with 2048 starting points and iterate times in the SHGO algorithm for the minimizer detection. We tested the number of starting points with values [1024, 2048, 4096, 8196] on various high-dimensional continuous functions, and the setting with 2048 starting points achieves the best efficiency. In datapoint sampling for counter-exampe identification, for each candidate expression, 800 data points are sampled from each (x and B, (x)), and additional 800 data points are randomly sampled across the state space D.
C. Risk-seeking Policy Gradient
Objective. The standard policy gradient (a) = V..l [R(V)) aims to optimize the average performance of policy given the reward function R(-). However, for the task of Lyapunov function construction, the final performance is measured by identifying single few valid analytical Lyapunov functions that satisfy the Lyapunov conditions, Thus, (Ø) is not an appropriate objective, as there is mismatch between the objective being optimized and the final performance evaluation metric. To address this misalignment, we adopt risk-secking policy gradient (Petersen et al., 2020), optimizing the best-case performance via the objective (á, at), as defined in Equation (6). In implementation, we choose
= 1 in the training all tested dynamics.
Proposition C.1 (Petersen et al. (2020)). Let (n, denote the conditional expectation rewards above the (1 a)-
quantile R. ()) as in Equation (6). Then the gradient of Jrie( a) is given by:
V Vodrialth ()) = Beg-p(Vslo.((e)) [ (R) -R(Va)) V logp(Vo o. (.)()) R(Vo) 2 (1)
(8)
15

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Algorithm 2 Giobal-oppiimizztion----eed Numerical Verification
Input: A set of analytical expressions = =[V* = =, Q). radius r, and state space D.
Output: a set of numerically valid candidate y/* a set of encountered counterexample A'e.
V. eee < 0.00.
2: for i = to Q do
3: zi, z2 - --SHHOOVi D), V D), (Identify global minimizers within the state space D)
4: 12, t {zil 2; E B, (z)), {zyl 2j E (x2)), {zyl 24 E D)
5: Check Lyapunov conditions on 21 U22 U33,
6:
if R(V)) = and no counter example found in 2i U U then
7:
V + 1 [VV)
8: else
9:
Xec Nx U identified counterexamples in Xi U U22 U J. (Gather falsification)
10: end if
11:
end for
12:
Return V*, Nee.
The proposition suggests a Monte Carlo estimate of the gradient of Jrik(Ø, a) from a batch of N samples:
()))
Σ R. R(V( ( (o V log (V) f(x)), aN
(9)
where Ra (Ø) is the empirical (1-o) -uuannile of the batch of rewards, and 1g returns if condition is true and otherwise. Compared to standard REINFORCE algorithm (Williams, 1992), Equation (9) has two distinct features: (1) it has specific baseline, R. (ø), instead of an arbitrary baseline in standard policy gradients chosen by user; (2) the gradient computation only uses the top a fraction of samples.
Lemma C.2 (Bastiani et al. (2024)). When using the empirical Lyapunov risk in Equation (4) as the reward functionn, the risk-seeking policy gradient is not guaranteed to be unbiased.
Proof. For given dataset X' C 'D, the empirical Lyapunov risk is random variable whose value is determined by the random expression sampled from the symbolic transformer model. Let Z denote the negated empirical Lyapunov risk, with its probability density p(z(Ø) depending on transformer parameters db.
Z has range (-o0, 0J. Suppose negated empirical Lyapunov risk is used as the reward function for risk-secking policy gradient optimize the average performance of top - -qqunnteee samples. In Equation (9), the risk-seeking policy approximates the gradient of the expectation over the truncated random variable Z., Z Izzs, with respect to where Sa is the (1 -u-unnteile of distribution of Z.
84 inf(= CDF(=) a).
The probability density of Z, given by
1 -1 a
p(za)
To effectively penalize the invalid sampled expressions that miss some state variables or are symbolically invalid, their empirical Lyapunov risks are set to be oo, or a symbolically valid candidate that incorporates all state variables might be less preferable than symbolically invalid criaall-incomplee expression during the training. Consequently, 8, has
16

---

Analytical Lyapunov Function Discovery; An RL-based Generative Approach
non-zero probability density of being -00. If 8a -c, then the gradient of the expectation in Equation (8) is:
1 0 zp(2 10)d2
B|Za|
B|Za] 1
2p(2 p2(20))2,
a
V
V2BlZal = V zp(2|0)d2.
Since the integration lower- bound is -oo, the Leibniz rule does not apply, which means interchanging the gradient and the integration changes the result. Therefore, the policy gradient in Equation (8) is no longer guaranteed to be unbiased.
Reward Design. To optimize the symbolic transformer parameters & such that the decoder generates a valid candidate Lyapunov function V satisfying Lyapunov conditions, we employ empirical Lyapunov risk as the fitness metric to measure the violation degree of Lyapunov conditions within the state space following Chang et al. (2019). However, shown in Lemma C.2, directly using the unbounded empirical Lyapunov risk as the reward for risk-seeking policy gradient might introduce bias. To address this issue, we adopt bounded reward function using the continuous mapping g(x) (Petersen et al., 2020; Bastiani et al., 2024), defined as:
R(V=) - g(((vv) 1 L(V=)
where ((V)) measures the violation degree over the training set X. This design ensures the reward is bounded [o, 1]. avoiding bias in the risk-seeking policy gradient.
Reward Calculation. The reward function in Equation (5) quantifies the degree of violation based empirical data points. However, as the dimensionality of the input dynamics increases linearly, the training set X' must grow exponentially to maintain precision in the empirical measurement of violation degree, which is impractical for high-dimensional cases. To ensure the reward signal remains a reliable indicator of expression quality in the risk-seeking policy gradient without requiring an excessively large dataset, we incorporate Projected Gradient Descent (PGD) (Madry et al., 2017) into the reward caleulation process. Prior to evaluating the reward, PGD is employed on sampled expressions to efficiently identify a set of "minimizers" in the state space. These minimizers are computed from randomly sampled set of starting points through PGD for each sampled expression, and the process can be parallelized on a GPU, enabling efficient computation. Though PGD may converge to local minima, it remains effective in identifying a few counterexamples that violate Lyapunov conditions across invalid candidates. These counter-example are then added to the training set X' for reward calculation of all sampled expressions in the current epoch and are removed immediately after the calculation. This process leverages PGD to capture violation data points in advance, ensuring the quality of the reward signal without the need for an excessively large dataset.
D. Genetic Programming
In the field of symbolic regression, given the large, combinatorial search space, traditional approaches commonly utilize evo- lutionary algorithms, especially genetic programming (GP) (Koza, 1992). to retrieve analytical expressions that approximate the output values given input data æ. The GP-based symbolic regression operates by evolving the input population of mathematical expressions through evolutionary operations such as selection, crossover, and mutation. A pre-defined fitness metric serves as the objective function to guide the optimization the population over successive generations. However, for analytical Lyapunov function construction, GP algorithms lack the capability to directly generate Lyapunov functions from the given dynamics and require an initial population that represents potential Lyapunov functions.
As the search space grows exponentially with the expression complexity and the number of states in input dynamics, it is a challenging task even for the symbolic transformer model to search a valid Lyapunov function for complex, high high-dimensional systems. Inspired by Mundhenk et al. (2021), we incorporate a GP component into the training framework to complement the symbolic transformer model the symbolic transformer model outputs a well-behaved initial populations of expressions Vo. which serve as the starting points for the GP component, and GP component refines V. through evolutionary operations
17

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Algorithm 3 Expert Guidance Loss
Input: 'Elite set' of analytical expressions V input dynamics f*)) and transformer parameters d.
Oulput:
The weighted cross-entropy loss between the transformer model output probability distribution and given refined expressions V spp-
1: G - pvpl {Get the size of 'elite set'}
2: [ +00, 0,
3: for = 1 to G do
4: - + bi
R(V
ap Eラ= - log (p(V,n, 0Ds VV p.
Φ. [()))). {Calculate the expert guidance loss based on 'elite
set' Vap- Equation (7)}
5: end for
6: c(Pop)= t gC.
7: Return g()
to explore the characteristics of Lyapunov functions that might be overlooked by symbolic transformer. The fitness metric for the GP component is the same as the reward function used in the risk-seeking policy gradient. After each refinement, we select an 'elite set' of the top-performing refined expressions, Vgp. based on fitness values. These expressions are treated as ground-truth truth solutions for the transformer decoder, and transformer parameters á are optimized through the expert guidance loss introduced Subsection 4.4. In the implementation, the size of 'elite set' D is chosen to be 0.1Q, where Q the number of sampled candidate expressions Vo in each epoch.
In our framework, we employ three evolutionary operations: mutation, crossover, and selection, within our Genetic Programming component (DEAP (Fortin et al., 2012)). A mutation operator introduces random variations an expression, such replacing subtree of one expression with another randomly generated subtree. A crossover operator exchanges content between two expressions, e.g., by swapping subtree of one expression with a subtree of another expression, enabling the combination of their features. A selection operator determines which expressions persist into the next population. common method is tournament selection (Koza, 1992), where a set of candidate expressions is randomly sampled from the population, and the expression with the highest fitness value is selected. each iteration of GP evolution, each expression has probability of undergoing mutation and a probability of undergoing crossover; selection is performed until the new generation's is population has the same size the current generation's population. In empirical experiments, set the probability of undergoing mutation and crossover to be 0.5, and we adjust the size of the tournament and number of evolutions proportional to the dimension of the input system.
E. Baseline Descriptions
F. Augmented Neural Lyapunov Control (ANLC)
The Augmented Neural Lyapunov Control (ANLC) (Grande et al., 2023) combines Artificial Neural Networks (ANNs) with Satisfiability Modulo Theories (SMT) solvers to synthesize stabilizing control laws for the input dynamics f()) with formal guarantees. The neural network is trained over a dataset of state-space samples generate candidate control laws and Lyapunov functions, while the SMT solvers are tasked with certifying the Lyapunov conditions of the neural Lyapunov function over continuous domain and returning counterexample if the function is invalid. To ease the computationally inefficient verification process in the SMT module, ANLC proposed discrete falsifier, which discretized the state space for sample selection and evaluation, employed before the SMT call to avoid the frequent calling of the consuminn SMT falsifier. As the previous learning-based Lyapunov function construction approaches usually initialized the parameters of control policy with pre-computed gains from state-feedback controllers, c.g. Linear- Quadratic Linear-Quadratic Regulators, which requires
user time and control expertise to properly perform the initialization process, ANLC instead removes the need of control initialization by its proposed compositional control architecture containing both linear and nonlinear control laws that the proposed method allows the synthesis of nonlinear (as well as linear) control laws with the sole requirement being the knowledge of the system dynamics, For empirical experiments, we tested the ANLC algorithm for all system dynamics in Appendices F and G. We tested on the Van Der Pol Oscillator and 3D Trig dynamics to get the best hyperparameter setting. In bold, we show the chosen parameters, selected to have the best success discovery rate on Van Der Pol Oscillator and 3D Trig Dynamics.
18

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Ir [0.1,
activations 1(=2,=* _x² x=), (tanh, tanh, x?), (a²,, xx², tanh), (tanh, tanh, tanh)]
hidden neurons- 66, 12, 15, 20)
data = [500, 1000, 2000
iteration = [500, 1000, 2000]
E.2. FOSSIL 2.0
FOSSIL 2.0 (Edwards et al., 2024) is a software tool for robust formal synthesis of certificates (e.g., Lyapunov and barrier functions) for dynamical systems modelled ordinary differential and difference equations. FOSSIL 2.0 implements counterexampl--guieed inductive synthesis (CEGIS) for the construction of certificates alongside feedback control law. In the loop of CEGIS, the learner, based upon neural network templates, acts candidate to satisfy the conditions over finite set D of samples, while the verifier (formal verification tools) works in a symbolic environment that either confirms or falsifies whether the candidate from learner satisfies the conditions over the whole dense domain A. If the verifier falsifies the candidate, one or more counterexamples identified by the verifier are added to the sample set, and the network retrained. This loop repeats until the verification proves that no counterexamples exist or until timeout reached. Similar the ANLC, in the empirical experiment, we set the hyperparameters based on the Van Der Pol Oscillator and 3-D Trig dynamics and tested for all other dynamics in Appendices F and G.
Ir 0.1, 0.01,
activations = [(=²,=²), (tanh, tanh, ,²), (tanh, x²)
hidden neurons [6, 10, 12]
data [500, 1000, 1000,2000]
iteration [25, 50, 100]
E.3. Global Lyapunov Function Discovery by Pre-trained Transformer
Alfarano et al. (2024) pre-trained a transformer on backward backward-generated and forward forward-generated global Lyapunov function datasets. The backward-generated datasets involve sampling arbitrary positive definite functions and deriving corresponding stable dynamics through some specific symbolic designs, while the forward-generated polynomial datasets contain randomly generated dynamies with corresponding Lyapunov functions identified by SOS methods the system inherently globally stable. Candidate Lyapunov expressions are sampled using beam search in a bakeneb-tokken manner. However, their method cannot adaptively refine the candidate Lyapunov functions if the beam search fails specific dynamics, and it requires dataset that is expensive to generate (e.g., thousands of CPU hours for 5-D dynamics dataset) to achieve adequate generalization during inference. Furthermore, its emphasis on global stability limits its applicability to real-world, nonpolynomial control systems, which typically only admit local stability. Due to the lack of resources of multiple industrial- level GPUs, we contacted the authors of A Alfarano et al. (2024) to conduct the evaluation their rre-trained model on our test systems, which is shown Section 5.
E.4. Sum-of-Squares (SOS) Methods
SOS methods formulate Lyapunov functions discovery of given dynamics as a semi-definite programming task, where the
coefficients a pre-defined SOS candidate expressions are optimized to satisfy the Lyapunov conditions (hard constraints in the optimization problem) using convex optimization tools. SOS methods are generally applied to polynomial systems for stability analysis. With proper recasting techniques, SOS methods also be applied to non-polynomial systems.
Definition E.1 (Sum of Squares, Papachristodoulou & Prajna (2005a)). For E R" multivariate polynomial p(z) is sum of squares (SOS) if there exist some polynomials f((x), = 1.... . such that
M
p(x) = Σ f(xx2
19

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Polynomial systems. By Papachristodoulou & Prajna (2005a), for a given n-dimensional polynomial dynamics f(x)) and integer degree 2d, check the globally asymptotical stability of f())), SOS method aims find polynomial V (xx) of degree 2d, such that
1. V(x) - E i= Ej=1 Gja 2/ SOS, where Σj=1 Gj > 2, n. with y > O, and Gjj 2 0 A and j.
2. -
av - dr fz) is SOS.
For local stability analysis, consider ball of radius centered at origin (0), which can be represented by the semialgebraic set S (x g(z,r) r) O, where g(z,r) = 2 x We require that the stability condition holds in S. Retaining the same optimization objective and constraints on (ar) as before, a modified constraint on Lie derivative is imposed: - - fx) -ss-ss()))))) (x)) r) is a SOS for some SOS s(z). If such an s(x) exists, we can establish local stability.
In Section 5, we develop our code based on the findlyas yap function from SOSTOOLS (MATLAB) and issue- 16 of SOSTOOLSS official GitHub repo to examine the SOS method on polynomial systems in Appendix F. Table summarizes the experiment results of SOS approach on our polynomial test dynamics.
Table 5. Training) solving time of sum-of-squares (SOS) test polynomial systems.
Systems
App. F.1
App. F.2
App. F.3
App. F.3
App. F.4
App. F.5
App. F.6
Degree 2d
2 B (0)
2 Global 0.497s
4 Global 2.509s
2 (0)
2 B (0)
Region Runtime
Global 0.832s
(0)
0.697s
Pon-Polynomial System. To apply SOS method on non non-polynomial system á = f(z), E D1. define x1 2 as the original states and 22 as the newly introduced variables to recast non-polynomial terms in f()). Let x = (x1, 22). The system dynamics can then written in rational polynomial forms:
f =f((),
i2 = f2(x),
with constraints x2 F(x1), (z) 0, G2(x) 2 0, where F, G1, G2 are vectors of functions, restrict the states of recast dynamics (2)) equal the manifold of the original state
Define g(x) as the collective denominator of f1, f2, and local region of interest semialgebraic set:
{xx) 1=+ (Gr(z) 2 0}.
where Gz(z) is vector of polynomials designed to match the original state space. Let 22,0 F(0). Suppose there exist polynomial functions Vzx), (z) (2)), and SOS polynomials 0,(2), i 1,2,3, 4, of appropriate dimensions such that
V(O, 22,0) = 0,
(10)
V(x) Xf ([]]))))))))())(( - 6(x) SOS,
(11)
0V -g(z)( dz (x)f((z++
av
sos,
(12)
Ouz
where d(x) is some scalar polynomials with F(I,,F(I1)) > 0, Vr E Dí M00. then 2 = is stable (not necessarily asymptotically stable) (Papachristodoulluu & Prajna, 2005b). In Section 5, we develop our code from SOSTOOLS to examine the SOS methods on two non- pol-poollo systems: simple pendulum and 3-D trig dynamics (Appendices G. G.2).
If constraints (10), (11), and (12) hold,
20

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
F. Polynomial Nonlinear Dynamical System
F.I Van Der Pol Oscillator
Van Der Pol Oscillator a nonconservative, oscillating system with nonlinear damping (Zhou et al., 2022). The dynamics of the Van Der Pol Oscillator have two state variables, formulated as follows:
$ 22,
52 -1 1e(1 -x7)= 222
where 21 and represent the object's position in the Cartesian coordinate, parameter Rt indicates the strength of
the damping. Under the state space D = {(01, x22) E R | S 1 and setting = 1, our proposed method found valid local Lyapunov function (2) = = + z.. Other forms of Lyapunov functions for Van Der Pol Oscillator, for example, (2) 2 + (2 + 22), are also recovered during the experiments.
F.2. www.warralle--oyymoiiaasyyssss with higher degree
Here we have a polynomial system two variables with a higher degree, adopted from Alfarano et al. (2024), formulated as:
d1 = - 5
2x1 2zz 92 + 3x
d2
D2 4x
Under the state space D {(x1, 22) E IR2 function (x11,2) = 9r; - 22.
|a,| < 1}, our proposed method successfully found valid local Lyapunov
F.3. Tiree-arribbeepyyyymiii--yyssss with higher degree
Table 6 describes two polynomial systems of three variables with a higher degree, adopted from Alfarano al. (2024). Our framework successfully retrieves valid local Lyapunov functions on both examples under the state space D {(11, 22,
R3 1xa| S 1}.
Table 6. Three--imensional Polynomial Example with Higher Degree.
System
Lyapunov fuunction
in -3x + 3m. .3 9
d² = x -522 + 5x3
V(11, = 9x + + +
53 9x9
i 8x1 x - 10x1
82 803 + 3m2 8
((11,22,2a) = =z .2
43 -33
21

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
F.4. 6-D Polynomial Nonlinear System
This 6-D dynamics consists of three two-dimensional asymptotically stable linear subsystems that are coupled by three nonlinearities with small gains adopted from Grüne (2019). The dynamics are written as:
$ -1 + 0.502 0.1
i -0.501 22.
f 23 + 0.5zI 0.1s
44 -0.503 24:
i5 =5 + 0.5ar.
ie 0.520 I6 + 0.1xI.
Our proposed method is trained over the state space D = 22 24, 25, ma) E R |>,|| S 1, Vi (1,2, 6)}, and able to find valid Lyapunov function (xx) zî + z + + 2 + + + This dynamics is not globally asymptotically stable since if ,2, or I3 has significantly large value, the perturbations introduced by the small gains will shift the object by significant amount away from the equilibrium point. By empirical checking, our found Lyapunov function certifies the asymptotical stability of this system over the region D' = {(01, T2, T3, IA, Is, 26) E R. > I? < 500).
F.5. 8-D Polynomial Nonlinear System
This 8-D dynamics consists of four two-dimensional asymptotically stable linear subsystems that are coupled by four nonlinearities with small gains, modified from the above 6D polynomial dynamics. The dynamics are written as:
è = 1 + 0.5gg - 0.1x9
4 = 0.5zz - 22:
23 = Σ3 + 0.5xz 0.1a².
4 = 0.52g - T4:
25 = 25 + 0.5zz + 0.1x7,
èg = -0.325 - I6;
87 = 27 + 0.5.rg:
èg = 0.52y - Ig - 0.1
Our proposed method is trained over state space D = [[x[[]]] x2, 2z, Es, 25, To I7, xz) E R | S 1, e 8}}. and is able to find a valid Lyapunov function V(x)) + +f 2z. This Lyapunov function certifies the asymptotical stability of this system over the region D = (111, I4, H5, Xig 27, x8) E & Σ; < 450J. This dynamics not globally asymptotically stable since if x1, 24, or Ø7 has significantly large value, the perturbations introduced by the small gains will shift the object by significant amount away from the equilibrium point.
F.6. 10-D Polynomial Nonlinear System
Finally, we extend to the original 10-D polynomial dynamics proposed in Grüne (2019). This 10-D dynamics consists of five two-dimensional asymptotically stable linear subsystems that are coupled by four nonlinearities with small gains. The
22

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
dynamics written as:
2 = -1 + 0.522 - 0.1a
22 = 0.5z - 22:
23 = Σ3 + 0.5zz 0.1aI
4 = 0.52z J4,
i5 = =5 + 0.5ml + 0.1ag.
dg = 0.52% - J6.
Ż7 = 27 + 0.52g:
Fg 0.5.r Zg.
ig .g + 0.5æ00
310 = 0.5.g - 310 0.1
Our proposed method trained over the state space = {(x1, 22 T3. F4, 45, T5, 27, 27, Fg. T10) E
R 10
|
S
1. Vi E (1.2..., 10],, and able to find a valid Lyapunov function V(x) = of + 2ÿ + z + x + 2 + 2ÿ + z0 This Lyapunov function certifies the asymptotic stability of this system over the region D = {(x1, I3, F4, rs, Te: 27, 2g, Ty, I1o) E RIO | E xi S 400). This dynamics is not globally asymptoti- cally stable since if æ 24, v5, or 2g has a significantly large value, the perturbations introduced by the small gains will shift the object by a significant amount away from the equilibrium point.
2
+ 2 +
G. Non-polynomial Nonlinear Dynamical Systems
G.1. Simple Pendulum
The simple pendulum is well-known classical nonlinear system that contains two state variables. The dynamics are formulated as follows,
+ 22
4 9 sin(x1)
b
m 322
where 2 the angular position from the inverted position, x2 the angular velocity, and parameters g, m, m, 6 are the acceleration of gravity, the mass of the inverted object, the length of the string, and the coefficient of friction, respectively. In experiments, since don't incorporate the constant generation capability within the training framework,
set g 1, 2 cos(#1) + x over the state space: D = (x1 x2) E 2 |x|| < and |zz|| < 63. This found Lyapunov function has the same analytical structure as the energy function of the inverted pendulum.
kg, 1 = 1 m, and 0.1. Our proposed method finds the valid Lyapunov function V
G.2. Trigonometric System
3-D trig dynamics comes from exercise problems in textbook Khalil (2002) whose dynamics written follows,
i = 2/2 4 -h(x1) - x2 h(x3),
d/ D2 F33
where l((x) = sin(z)-cos(z) When the state space is D (11,,
E ER
| < 1.5 (1,2,3)}, the valid local
Lyapunov function found by our proposed method is V(I1, I2, xz) = 1 cos(21)² + + + sin(az)², which is consistent to the textbook solution of Lyapunov function for this particular dynamics.
G.3. N-bus Lossless Power System
We test our proposed framework on the N-bus power lossless system (Cui et al., 2023a: Feng et al., 2024a) to examine its ability to handle complex high-dimensional dynamics. Consider 9g, W; as the phase angle and the frequency of bus i,
23

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
respectively, the dynamics for each bus formulated follows,
0 Wis
N
mWW; Ps dww; - ui(w;) Σ
B sin(0; - -8)),
where m, is the generator inertia constant, de is the combined frequency response coefficient from synchronous generators and frequency sensitive load, and is the net power injection, for each bus i = 1.... NN. B RNNN is the susceptance matrix with B.j- = 0 for every pair ji such that bus and bus j are not connected, and (w())) is the controller at bus i that adjusts the power injection to stabilize the frequency.
Since the frequency dynamics of the system depends only on the phase angle differences, so we change the coordinates:
N Σ
6, = 0,
0
N 1=1
where 8, can be understood as the center-of-inertia coordinates of each bus. In our experiment, we test the proposed framework on the 3-bus power system. For simplicity, we set p;
m; 2, de 1, u (w:) wi, and Big
1 1 # j, Ba
In this case, the equilibrium point for our system is at the origin, i.e. of w/ 0,
3, The state space for our experiment is defined as: D {(11, 62, 03,w1, W2) R [6l S 0.75 and laal S
1.2 for i 1,2,33. 33. Through our method, we retrieved valid Lyapunov function V(ū,, 82, 83, 41, W2, wa) (Ei w
0.5
cos(o, 5,) -1): which is consistent to the known Lyapunov function presented in (Cui et al., 2023a).
The Lie derivative of the identified Lyapunov function can be simplified as L, 22( + w + w.). The analytical structure of this found Lyapunov function and invariance principle allows us to easily identify it as valid Lyapunov function by hand.
G.4. Indoor Micro Quadrotor
For the angular rotations subsystems of the quadrator from Bouabdallah et al. (2004), it has 6 states to describe the angular motion of the quadrotor. The states T2, I3, and Øs describe the roll, pitch, and yaw of the quadrator, and states 22 F4, and de represent their time derivatives. With perturbation terms Ω and control inputs U1, U2, and U3, the subsystem can be formulated as follows,
1 = I2
1 12
JR
82 x436(
x
UI,
1x
T
T
49 -44
I, I
JR
21 IzT6
U2.
ln
ly
Ig
45 J6
1 12- U3.
1
26 12
1 I
where Jay Iy, and represents the body inertia, I denotes the lever, and Jn is the rotor inertia. With the control policy
L
U
[ (x) x11) kj
Ty
Un 1
(0z xz)
k234,
U3 - Lo(25 x k306.
and restricting Is Iy. the angular rotations subsystems stabilized to the chosen equilibrium point X₃ (xf. 0, ag. 0, .g, 0}. In empirical experiments, we set X, (0,0,0,0,0,000 state space D ((x1, T2, F3, 3., Rs, T6)
24

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
2, 1, 5, k1 = 20, - our framework successfully found Lyapunov function (a 22, 23, Fa. 26) )) = -ا= xº. By examining the Lie deriva- tive L 2.533 - 10z - 0.8v and using the Invariance principle, we conclude that this subsystem globally asymptotically stable. Under other parameter settings, our framework can also retrieve valid analytical Lyapunov functions for this dynamics.
R$ [ S 3, Vi E (1, 1,
1, JR 1, and Ω sin(x)) cos(
G.5. N-bus Lossy Power System
Unlike N-bus lossless power systems in Appendix G.3 which has a well-known mergy-based storage function served as a valid Lyapunov function for asymptotical stability guarantee, the N-bus lossy power system (Cui & Zhang, 2022) does have a known analytical Lyapunov function to certify stability, though by passivity the system should be asymptotically stable at origin. Utilizing the proposed framework, we aim to discover valid analytical Lyapunov function for an N-bus lossy power system to formally certify its asymptotic stability.
The 0g and ó; are the angle and frequency deviation of bus i, the dynamics an N-bus lossy power system is represented by the swing equation, formulated as:
0,-
m,i: ps= diw, u:(wi)
Σ Bi sin(0; -0,)
Σ Gy cos 0z),
where m; is the generator inertia constant, de is the combined frequency response coefficient from synchronous generators and frequency-sensitive load, and p; is the net power injection, for each bus i = 1, ...... N. The susceptance and conductance of the line (i, are B.j- = By, and Gi5- = Gjt. respectively. The value is if the buses are not connected. In this work, we consider input uy to be a static feedback controller where only its local frequency measurement wi is available. Like the lossless power system, since the frequency dynamics of the system depends only on the phase angle differences, we change the coordinates:
-
O
N
where 6 can be understood as the center-of-inertic coordinates of each bus.
We test the proposed framework on b-s lossy power system. In experiment, we set p; 1,m; = 2, d; = 1,, uz (w;) Wi.
B., = 1, Giy = i + j. Ba = O, Gij = O. By this setting. the equilibrium point for this system is at the origin. i.e. 6;
w = 0. The state space for the experiment is defined as: D = {((1,, 82, w2) E R4 | < 0.75 and |لس| < 2 for
1, 2). The proposed method found two valid Lyapunov function (88 82, w122) = w + w + (w2 - sin(o,) + )22)
and Vō11, 62, w1, w2) = w + w3 + (-w1 sin(61) + sin(62))2 Both Lyapunov functions pass the formal verification by SMT solver in the state space D (0), where precision is set be -12 and - to avoid tolerable numerical error.
G.6. 9-D Synthetic Dynamics
Consider the synthetic dynamics adapted from Appendices F.4 & G.2 with linear interactions between two subsystems:
1 = -1 + 0.5x% - 0.10I
22 = -0.5z, - D2 + 0.128.
23 = -T3 + 0.5rz 0.1æ?.
4 = -0.523 - F4,
d5 = -x5 + 0.5x6,
de = -0.50g - ve + 0.1x2
Φ7 = Fs+
ds = sin (I7) cos(17) D's sin(rg) cos(xg) 0.1æ2
dbg Ig Zg.
25

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
1.05
a=0.1
1.0
0.1
1.00
0-5.5
0.95
Reward Values
Reward Values
0.90
0.6
0.85
0.75
0.70
0.65
0.0
Epoch
Epoch
Figure 5. Reward Trajectory Comparison for Ablation Study I: Test proposed framework (without GP component) 3-D Trig dynamics with three different choices of for 100 epochs. The figure on the left visualizes the highest reward the batch of samples obtained in each epoch, and the figure the right visualizes the 90% quantile of reward distributions of the hatch of samples obtained in each epoch.
To properly address the trigonometric terms in 2g. the Lyapunov function for this dynamics can't 't be a simple form like Ex . and should include some trigonometric terms. Setting the state space D {ir IR |x:|| < 1.5, V 1, 9). our method successfully identifies valid Lyapunov function V = Ein x² sin(x7)2 + x cos(z9) + 1, which passes formal verification following settings in Section 5.
H. Ablation Studies
H.1. Ablation Study I Risk-seeking Policy Gradient
To better understand the influence of a in policy gradient on the training performance of the symbolic transformer, we experiment with our framework without GP refinement on 3-D Trig dynamics under different choices of (the risk-seeking quantile) over 100 epochs. Specifically, we choose a = 0.1, 0.5, and compare the transformer model's performance by the reward trajectories and success discovery rate. Figure 5 plots the reward trajectories- both the highest reward and the 90%-quantile- uunetereleudddrreeeeee the three different choices of a (the risk-seeking quantile). For a 0.1, the setting we used on tested dynamics in experiments, the highest reward and the 90% quantile converge steadily to and 0.9, respectively, exhibiting lower variance and faster stabilization. For 0.5, the median-based -ased threshold increases steadily but shows greater instability. Finally, for around 0.6. Additionally, with a = 0.1, the framework achieves the highest success recovery rate (66.67%) compared to 33.33% and 0% for a = 0.5 and a = 1. These findings correspond to the analysis of misaligned objective between the standard policy gradient and Lyapunov function construction in Subsection 4.3, confirm the importance of risk-seeking policy optimization, and highlight the effectiveness of setting a 0.1 experiments.
vanilla policy gradient, the best reward reaches only 0.94, while the 90% quantile hovers
H.2. Ablation Study Global-optimizationn- based Numerical Verification Process
This work proposes a global-opiimization based numerical verification for candidate verification and countercxamples' feedback within the training paradigm. To identify the effectiveness the proposed verification algorithm, this ablation study compares the following verification methods on the proposed framework without the GP component on the 6-D polynomial dynamics (Appendix F.4):
Global optimization (SHGO): Proposed framework with global-optimizationn based numerical verification in Subsec- tion 4.2.
Root finding: Proposed framework with root-finding numerical verification introduced in Feng et al. (2024c).
SMT: Proposed framework with formal verification dreal (Gao al., 2013) SMT solver.
Ramdom sampling: Proposed framework with random sampling verification.
26

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Table Number of epoch, training time, success rate, and formal verification time comparison between four different settings in Ablation study II. Runtime is the average training time for successful trials. The Suec. % is the success rate of finding valid Lyapunov function out random seeds. Time the average formal verification time identified candidates on successful trials. The experiment lerminates either when valid Lyapunov function is found (passes formal verification), the training epoch exceeds 250, timeout (4 hours limit).
Stats.
Global optimization
Root finding
SMT
Random sampling
Epoch
168
160
Runtime
8399s
5684s
Suce. %
100
0
0
100
Ver. Time
1.67s
15.93s
1.1
1.0
0.9
Reward Values
Reward Values
0.8
MW
0.7
0.6
0.6
0.5
0.5
0.4
Root Finding
SHGO
0.4
Tanstomer
Random
0.3
25
50
75
100
125
150
175
200
20
60
80
100
120
140
160
Epoch
Epoch
Figure Best Reward Trajectory for Ablation Study H: Visual- ize the best reward trajectory of our framework with SHGO, root finding, and random sampling schemes in Ablation Study Il on 6-D polynomial dynamies. The SeT-veiifiattion based trajectory not included as it timeouts at the beginning training stage.
Figure 7. Best Reward Trajectory for Ablation Study III: Visualize the best reward trajectory the proposed framework under four different settings in Ablation Study III on 6-D polynomial dynamics.
Figure 6 compares the highest reward trajectory in each epoch for all four settings, and Table 7 summarizes the training and formal verification statistics. SMT is not included in Figure as the training stuck on the SMT solver for hours in the middle of the training stage. The proposed SHGO-hased counterexample feedhack makes the training reward an adversarial reward, as it can identify the most critical violations. As a 'challenger', the global optimization method finds the adversarial counterexamples, not just mild violations. Consequently, this leads to faster refinement the policy to mitigate the critical violations and results in stronger guarantees the final Lyapunov candidates. However, it risks training stability and may over-focus 'hard" violations. In this case, it identifies valid Lyapunov function V1 2t 2 + + x + = + a which is well-structured and easy verify by formal verification tools. Note that, unlike in the experiment section, the GP component is not used here, resulting in a less structured local Lyapunov function with few squared interaction terms (x: + )2 under one random seed. This increases verification time compared to Table 1, yet remains significantly faster than verifying the final result from the random sampling scheme.
Among the four settings, random sampling achieves the fastest convergence for this polynomial system, as it yields smoother rewards, covers the whole state space uniformly, and promotes training stability. However, because random sampling cannot actively search for violations and its search space grows exponentially with dimensionality for larger or more complex systems (e.g., power systems). Even with reward near 1, it may fail to find a valid Lyapunov function as the high reward only means it fails to locate counterexamples instead of nonexistence; or if it does, it can converge to a local minimum with poor generalization or verification properties. For instance, the discovered Lyapunov function V2(2)) = 2 - cos((2)) cos(xr) cos((4) cos(2s) cos(æ; + +x1) +4 has a complex shape that formal verification tools struggle handle, and it offers limited insight beyond the training region D. By contrast, V, (a) remains valid over the
becomes less effective
27

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
Table Number of epachs, training time, and success comparison between four settings of the proposed framework. Runtime is the average training time for successful trials. Suec. % is the successful rate of finding valid Lyapunov function out of 3 random seeds.
Training Stats.
Transformer only
Transformer + GP
Ours
GP only
Epoch
168
18
Runtime
8399s
2731s 100
2467s
Suce. %
100
100
much larger domain D' = ( T3: 24, 2i, va) E TR
1i=1 2 S 500).
Root-finding usually identifies counterexamples near function roots, producing mostly mild violations. Although an ideal root-finding algorithm might be efficient, off-the-shelf methods become inaccurate for dimensions above five. As a result, they fail to provide precise violation signals or strong optimization guidance to the transformer, causing the best reward stall 0.7. Similarly, dreal SMT often locates only mild violations rather than critical ones, limiting its impact on optimization. While these counterexamples help the transformer reduce violations, the refined transformer increasingly acts adversary the SMT solver, prolonging verification until timeout. In neural network formal verification settings, this effect is even exacerbated by the model's overparameterizztion. Consequently, we rely on global optimization instead of formal verification in the training loop.
H.3. Ablation Study III Generic Programming and Expert Guidance Loss
In the proposed framework, the Genetic Programming (GP) component refines sampled candidates V& from the symbolic transformer and optimizes the transformer parameters & using the expert guidance loss in Equation (7). To assess the impact of the GP component on the training paradigm, we evaluate the following four settings of the proposed framework on the 6-D polynomial dynamics (Appendix F.4):
Transformer only: A pure symbolic transformer model without GP component.
Transformer + GP: symbolic transformer model with GP component for candidate refinement but without expert-guided learning.
Transformer + GP + expert-guided learning (ours): A symbolic transformer model with GP component for both candidate refinement and expert-guided learning for transformer optimization.
GP only: GP component only. The initial population of the first epoch comes from randomly initialized symbolic transformer, and the initial population of other epochs comes from the latest refined expressions.
Figure 7 shows the highest reward in each epoch for the four settings, and Table compares their training statistics. The transformer-only approach achieves a 100% success rate across three random seeds, proving its capability to generate Lyapunov functions on high-dimensional systems. However, without the GP component, it requires significantly more training time and epochs identify a valid Lyapunov function because it must independently explore all necessary characteristics of valid solution. Conversely, the GP-only approach fails to find any valid Lyapunov functions when starting from a random population, it lacks the transformer's deeper understanding of system dynamics and relies on less capable evolutionary operations for exploration.
Figure explicitly compares the reward trajectories with and without expert-guided supervised learning, where ours (with expert-guided learning) achieves the fastest convergence. Comparing the reward trajectory. ours achieves faster and more stable training. Overall, these visualizations and comparisons emphasize benefits of GP refinement and the expert-guided learning on 'elite set' V 3p accelerate transformer parameters' update and enable faster Lyapunov function discovery.
28

---

Analytical Lyapunov Function Discovery: An RL-based Generative Approach
1.1
1.1
Values
Reward Values
0.9
0.8
0.8
0.7
0.6
0.6
Transfomer GP (refine)
Transfomeer GP trefina)
10
15 Epoch
20
25
30
10
15 Epoch
20
25
30
Figure 8. The figure the left visualizes the overall best reward (best candidate among V. U Vpp)) trajectory. The figure the left visualizes the best reward trajectory for transformer candidates (best candidate among Vo only, excludes GP-generated solutions). These plots compare the performance of transformer GP refinement and method Ablation Study III.
29


Reviews:
[{"review_id": "YxpXNXYVFB", "reviewer": "Reviewer_SwEU", "summary": "This paper presents a novel and promising approach to discovering analytical Lyapunov functions for nonlinear dynamical systems using reinforcement learning and transformer-based generative models. The work addresses a fundamental challenge in control theory with significant practical implications. The framework consists of three key components: (1) a symbolic transformer that generates candidate Lyapunov functions, (2) a global-optimization-based numerical verifier that checks Lyapunov conditions and provides counterexamples, and (3) a risk-seeking policy gradient algorithm that optimizes the transformer parameters. The approach is enhanced with genetic programming for expert guidance.\n\n## update after rebuttal: I have decided not to change my assessment and maintain my score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer SwEU for the valuable time and constructive feedback. We provide point-by-point responses below.\n\n**Q1: Alfarano et al. (2024) was excluded from direct comparison**\n\nTo complete the comparison, we contacted the authors of Alfarano et al. (2024), who conducted the evaluation of their pre-trained model on our test systems. Due to constraints related to global stability and the dimensionality of their training data, their model can only successfully find Lyapunov functions for systems in Appendixs F.1, F.2, F.3, \\& G.1. These systems have dimensions lower than 6 and the found Lyapunov functions offer guarantees of global stability. We will include these results in the experiment section of the final paper.\n\n**Q2: No comparison with traditional symbolic regression or optimization-based methods.**\n\n1. We conducted a comparison to the Genetic Programming (GP) algorithm, a representative symbolic regression method, on 6-D polynomial dynamics. See Appendix H.3 of the manuscript. Ours can achieve a 100\\% success rate within 20 epochs, while GP alone failed all trials. As shown in Figure 7, GP initially improved generation quality (measured by Lyapunov risk reward) rapidly from random initialization, its evolutionary exploration failed to identify any valid Lyapunov functions due to insufficient incorporation of system dynamics.\n\n2. We added the sum-of-square method as a baseline for optimization-based method. In our experiments, SOS techniques successfully identify Lyapunov functions for polynomial systems of dimension up to 3 (Appendix F.1, F.2, \\& F.3) but failed to retrieve any valid local Lyapunov functions for polynomial test systems in Appendix F.4, F.5, \\& F.6 of dimension $\\geq6$. Due to the character limit, please refer to Response to Reviewer s9Wv (Q1) for the detailed comparison setup and results.\n\n3. In addition, we apply the SOS method on non-polynomial systems by recasting the non-polynomial dynamics to polynomial dynamics. It successfully identifies valid Lyapunov function for the pendulum system (Appendix G.1) but suffers scalability issues for complex high-dimensional non-polynomial systems and requires substantial domain expertise in the recasting process and (in)-equality constraints design. Please refer to our response to reviewer uz6L (Q1) for detailed results and analysis.\n\nWe will highlight the comparisons to classic symbolic regression methods (GP) and add the comparison to optimization-based SOS method in the final version.\n\n**W1: It's difficult to understand which part provides the most gains.**\n\nThanks for the valuable feedback. The backbone of our framework is a symbolic transformer trained via a risk-seeking policy gradient algorithm to efficiently discover Lyapunov functions, which provides the most gains. To clarify the gains of other parts of the model design, we performed detailed ablation studies in the appendices. The global optimization module is deployed to evaluate the generated expressions and guide the training, as shown in Appendix H.2 (Figure 6). The GP module provides expert guidance to accelerate training, illustrated in Appendix H.3 (Figure 7 \\& 8). Together, the global optimization and GP modules assist the training of the symbolic transformer. \n\n**Q4: Exploration-exploitation trade-off in your approach:**\n\nThanks for this insightful question. In our framework, exploration occurs at the candidate expression generation step and GP's evolutionary operations, while risk-seeking policy gradient and expert guidance loss are responsible for the exploitation. \n\nAs symbolic tokens in candidates are sampled according to some learned conditional probability distribution (the output of the decoder), and GP algorithms employ mutation, crossover, and selection operations to introduce randomness for the search of higher quality candidates, exploration is encouraged in these two steps.\n\nFor exploitation, the risk-seeking policy gradient optimizes around high-quality expressions, reinforcing high-reward candidates to focus the search on promising regions. The expert guidance loss also facilitates exploitation by encouraging the learned distribution to match the best-known candidates. We tuned the risk-seeking hyperparameter $\\alpha$ to achieve the trade-off. In practice, we find $\\alpha=0.1$ achieves the best performance, details are available in Appendix H.1.", "claims_and_evidence": "Primary claim: The RL-based framework can discover valid Lyapunov functions for complex systems\n\nEvidence: \n1. The method achieves 100% success rate on most low-dimensional systems and maintains 60-80% success on higher-dimensional systems\n2. The framework successfully handles non-polynomial systems with trigonometric terms\n\nSecondary Claim 1: The approach discovers previously unknown Lyapunov functions\n\nEvidence: \n1. For the 4-D lossy power system, they discover two distinct valid Lyapunov functions\n2. The authors claim these are the first analytical Lyapunov functions for certifying local stability of a 2-bus lossy power system\n\n\nSecondary Claim 2: Risk-seeking policy gradient improves discovery performance\nEvidence: \n1. Ablation study (referenced in section 5.5) shows risk-seeking with α=0.1 achieves 66.67% success rate compared to 33.33% with α=0.5 and 0% with standard policy gradient (α=1)", "methods_and_evaluation_criteria": "The paper doesn't use standard benchmark datasets but rather a collection of dynamical systems of varying complexity. It compares against couple of baselines (ANLC and FOSSIL). Formal verification is done on discovered Lyapunov functions using various methods.", "theoretical_claims": "I did not check for the correctness of any proofs or theoretical claims.", "experimental_designs_or_analyses": "The experiments are quite thorough and the design seems valid. The evaluation metrics are clear and the ablation studies are well designed. \n\nLimitations of the experiments: \n1. Alfarano et al. (2024) was excluded from direct comparison due to its focus on global Lyapunov functions\n2. Limited to only two neural network-based baselines\n3. No comparison with traditional symbolic regression or optimization-based methods", "relation_to_broader_scientific_literature": "The topic is very much related to learning based search and can be applied to many areas of sciences and not just control theory. The approach could be adapted to discover invariants and correctness certificates for software verification. Instead of dynamical systems, the framework could analyze program state transitions and generate invariants that prove program correctness. Automated theorem proving could be another area of application.", "essential_references_not_discussed": "None, that i know of.", "other_strengths_and_weaknesses": "Weakness: \n1. The method is quite complicated and contains multiple parts. Its difficult to understand which part provides the most gains.", "other_comments_or_suggestions": "Could not find typos.", "questions_for_authors": "I am curious about the exploration-exploitation trade-off in your approach: Since the space of potential analytical expressions grows exponentially with expression complexity, how did you address the exploration challenge in your RL framework?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3HlFDSIbar", "reviewer": "Reviewer_uz6L", "summary": "The paper presents an RL framework for discovering analytical Lyapunov functions for nonlinear dynamical systems. The proposed approach trains a symbolic transformer from scratch. The transformer generates candidate Lyapunov functions in a symbolic form, which is then refined and verified via a combination of global optimization (using SHGO), risk-seeking policy gradient updates, and genetic programming (GP) as expert guidance. The authors compared their approach with two other tools over a series of benchmarks.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1 and Q2: Justify the claimed advantage over SOS techniques. Sensitivity issues of SOS. Completeness of their approach.**\n\nThanks for this valuable advice. We tested SOS using SOSTOOLS on polynomial (poly) and non-polynomial (non-poly) systems. Below is a summary of the setup and results.\n\nPolynomial systems: SOS techniques successfully identify Lyapunov functions for poly systems up to 3-D (Appendix F.1, F.2, \\& F.3) but failed to retrieve valid local Lyapunov functions for systems in Appendix F.4, F.5, \\& F.6 of dimension $\\geq6$. Please refer to our response to reviewer s9Wv (Q1) for detailed setup and results.\n\nNon-polynomial systems: Following [3], we apply the SOS method on two non-poly systems: the simple pendulum and 3-D trig dynamics (Appendix G.1 \\& G.2). Both are recast into poly form by introducing new variables and necessary (in)equality constraints. SOSTOOLS identifies valid Lyapunov functions for both cases (see Table). For the pendulum, it succeeds under 20s. However, for the locally stable 3-D trig system, the addition of one state introduces four extra constraints, and the resulting formulation leads to over an hour of solving time, in contrast to ours’ 157s.\n\n## Table: Performance of SOS Method on Non-Poly System\n| **System**| **Degree 2d** | **Runtime** | **Region**|\n|-|-|-|-|\n| App. G.1|2|19.58s|Global|\n| App. G.2|2|6163s|$\\\\{x\\in\\mathbb{R}^3\\|\\|x_i\\|\\leq1.5\\\\}$|\n\nTo apply SOS methods on non-poly system $\\dot{z}=f(z),z\\in\\mathcal{D}_1$, define $x_1=z$ as the original states and $x_2$ as the new variables. Let $x=(x_1,x_2)$. The system dynamics can then be written in rational poly forms: $\\dot{x}_1=f_1(x),\\dot{x}_2=f_2(x),$ with constraints: $x_2=F(x_1),G_1(x)=0,G_2(x)\\geq0,$\nwhere $F,G_1,G_2$ are vectors of functions.\n\nDefine $g(x)$ as the collective denominator of $f_1,f_2$, and local region of interest as a semialgebraic set: $$\\\\{(x)\\in\\mathbb{R}^{n+m}|G_D(x)\\geq0\\\\},$$ where $G_D(x)$ is a vector of poly designed to match the original state space.\n\nLet $x_{2,0}=F(0)$. Suppose there exist poly functions $V(x)$, $\\lambda_1(x),\\lambda_2(x)$, and SOS poly $\\sigma_i(x),i=1,2,3,4$ of appropriate dimensions, s.t.\n$$V(0,x_{2, 0})=0,\\\\;(1)$$\n$$ V(x)-\\lambda_1^T(x)G_1(x)-\\sigma_1^T(x)G_2(x)-\\sigma_3^T(x)G_D(x) -\\phi(x)\\in\\text{SOS},\\\\;(2)$$\n$$-g(x)(\\frac{\\partial V}{\\partial x_1}(x)f_1(x)+\\frac{\\partial V}{\\partial x_2}(x)f_2(x))-\\lambda_2^T(x)G_1(x)-\\sigma_2^T(x)G_2(x)-\\sigma_4^T(x)G_D(x)\\in\\text{SOS},\\\\;(3)$$\nwhere $\\phi(x)$ is some scalar poly with $\\phi(x_1,F(x_1))>0,\\forall x_1\\in\\mathcal{D}_1\\backslash \\\\{0\\\\}$. If (1), (2), and (3) hold, then $z=0$ is stable.\n\nCompared with ours, the SOS method on non-poly systems has a few limitations: a) it requires substantial domain expertise for recasting ($G_D,F,G_1,G_2$) and manual design of coefficient constraints (ensuring $V(0)=0,V(x)>0$); b) additional state variables and rapidly growing constraint complexity hinder scalability—constraints (1), (2), and (3) are more complex than those in polynomial systems, which already struggle at dimension $\\geq6$; c) the approach does not extend to certifying asymptotic stability.\n\nAs for completeness, symbolic transformers have shown strong empirical performance in expression search, supported by recent work (Alfarano et al., 2024; Holt et al., 2023). However, if our method fails to return a valid Lyapunov function, it does not imply the system is unstable as other constructive Lyapunov function methods. Formal completeness proofs remain an open question for future research. We will clearly state this limitation in the final paper.\n\nSensitivity issues are noted in Dawson et al. (2023b). Since our work lacks a thorough analysis of it, we agree with the reviewer and will remove the claim.\n\n**Q3: soft verification**\n\nWe introduce soft verification to improve training efficiency and candidate quality. Specifically, our framework employs SHGO to find the minimizer of the Lyapunov function $V$ and the maximizer of the Lie derivative $L_f V$. We sample around these points to find counterexamples, that avoids the frequent timeouts observed with SMT-based verification during training and enhances Lyapunov quality (Appendix H.2). Formal SMT verification is applied post-training. Similar strategies include adaptive sampling (Grande et al., 2023) and projected gradient descent (Yang et al., 2024) to reduce verification costs.\n\n**Q4: Discussion of [2]**\n\nThis paper uses NN poly expansions to reformulate barrier function search from nonconvex bilinear problems to CEGAR training with LMI verification. It handles hybrid systems up to 3-D and continuous systems up to 23-D, achieving scalability by using SOS only for verification. However, it is limited to polynomial dynamics and candidates. We will include [2] in the final paper.\n\n[3] P. Antonis, and S. Prajna. Analysis of non-polynomial systems using the sum of squares decomposition. Positive polynomials in control. Springer Berlin Heidelberg. 2005.", "claims_and_evidence": "The authors claimed to address the challenge: \"Can neural networks effectively discover valid analytical Lyapunov functions directly from complex system dynamics?\"\n\nIt is partially supported by the experiment results, which show an improvement in scalability compared to ANLC and FOSSIL. However, all the benchmarks are addressable for classical SDP-based Lyapunov function synthesis techniques, e.g., using [1] to convert ODEs with elementary functions to polynomial ODE and using SOSTools with Sedumi to address them. Classical SDP-based techniques can typically handle the systems up to about 10 dimensions, e.g., in the demos provided by SOSTools, sosdemo5 is an 8-dimensional example using SDP, though it is not for Lyapunov stability analysis. Thus, the given experiments cannot fully support the claim.\n\nSimilarly, the authors claimed in the introduction that \"despite the progress of SOS methods, several theoretical results (Ahmadi et al.,2011; Ahmadi, 2012; Ahmadi & El Khadir,2018) on asymptotic systems may not agree with a polynomial Lyapunov function of any degree. In addition, SOS methods suffer numerical sensitivity issues in practice.\" No experiment is given to demonstrate that the proposed technique outperforms SDP/SOS in sensitivity issues. In addition, the authors did not provide the proof of completeness of their approach. Thus, it is not fair to blame the incompleteness of the SDP/SOS techniques.\n\nFinally, except for the classical SDP/SOS techniques, there are existing learning-based techniques, e.g., [2], that can address harder systems (hybrid systems) with up to 23 dimensions to generate analytical (Lyapunov-like) barrier functions. Considering the fundamental similarity between the Lyapunov function and barrier function, it is necessary to compare with these works.\n\n[1] Liu, Jiang, et al. \"Abstraction of elementary hybrid systems by variable transformation.\" International Symposium on Formal Methods. Cham: Springer International Publishing, 2015.\n\n[2] Zhao, H., Liu, B., Dehbi, L., Xie, H., Yang, Z., & Qian, H. (2024). Polynomial Neural Barrier Certificate Synthesis of Hybrid Systems via Counterexample Guidance. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 43(11), 3756-3767.", "methods_and_evaluation_criteria": "Overall, the proposed technique follows the paradigm of data-driven or counter-example guided abstraction refinement (CEGAR) techniques. \n\nHowever, there is a concern regarding Section 4.2. Section 4.2 introduces an intermediate soft verification -- I refer to it as soft verification as it does not provide hard guarantees and an additional SMT-based verification step is needed. It is not clear why this step is necessary. In existing approaches, no soft verification is used during the training, for instance, Chang et al.(2019), and hard verification is conducted when the training converges. Any specific reason to introduce this step into the framework?", "theoretical_claims": "This paper does not involve theoretical claims.", "experimental_designs_or_analyses": "As mentioned in the Claims And Evidence, it is expected to see the direct comparison with SOS/SDP methods over the given benchmarks. It is also expected to see the experiment on higher dimensional systems (>10), which is beyond the capacity of SOS/SDP methods. It helps demonstrate the advantages of the proposed technique.", "relation_to_broader_scientific_literature": "None.", "essential_references_not_discussed": "There is a large amount of work [2] on the analytical neural barrier function learning that needs to be cited, considering the relation between barrier function and Lyapunov function. In [2], they handled a 13D continuous system and a 23D hybrid system, which is a significant result compared to the examples considered in this work.", "other_strengths_and_weaknesses": "The idea of using a symbolic transformer is quite interesting, and it is expected to see deeper analysis and discussion in this part, e.g., convergence, (statistical) completeness.", "other_comments_or_suggestions": "As mentioned in the Other Strengths And Weaknesses, it is expected to see deeper analysis and discussion in symbolic transformer, e.g., convergence, (statistical) completeness. While the whole framework follows the classical paradigm, it is suggested to mostly focus on the symbolic transformer in this work.", "questions_for_authors": "See previous sections.\n\nThanks for the additional effort. I have gone through the additional response. My concern is fairly addressed, with some reservations about the local robustness. I have updated my rate. I suggest the authors further compare the latest SOS-based techniques on local robustness in the later version.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kRqDfdAu8H", "reviewer": "Reviewer_ve2p", "summary": "1. Similarly to other lines of work, authors train a new model for a single nonlinear dynamical system. \n2. Key innovation lies on the symbolic formulation of the problem: instead of training a fixed model M (which would serve as the Lyapunov candidate) as per previous line of work, they work with symbolic transformer. \n3. They employ an RL loop to discover local Lyapunov candidates. A good component introduced by authors is the SHGO algorithm.\n4. The most real-world result is the fact that (a) the paper is not constrained by dimensions, (b) they discovered new local Lyapunov functions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer ve2p for the valuable time and constructive feedback. We provide point-by-point responses below.\n\n**Q1: Given the symbolic answer provided by the model, can this help to enlarge the stability area given you are not overfitting on a particular locality domain?**\n\nThanks for the insightful question. \n\nFirstly, our approach implicitly promotes generalization beyond the local sample region. Since the symbolic transformer directly operates on analytical expressions of dynamics rather than finite samples from the state space, the resulting Lyapunov functions naturally generalize beyond the sampled region. Empirically, as demonstrated by the 6-D polynomial system (Appendix F.4), our method successfully identified Lyapunov functions valid across expanded state-space domains. Additionally, the Lyapunov function obtained for the 6-D quadrator (Appendix G.4) is globally valid. \n\nFurthermore, to explicitly enlarge the stability region, a promising future direction is to characterize the region of attraction (ROA) of the learned symbolic Lyapunov function and incorporate it into the reward function of risk-seeking policy gradient. Related efforts to enlarge the ROA in learning-based Lyapunov function search include Chang et al. (2019), Wu et al. (2023), and Yang et al. (2024). The corresponding references can be found in the reference section of the original manuscript.\n\n**Q2: It's not clear how much general can the approach be: due to the nature of symbolic representation, this forces the candidate to be always the same function in the entire domain. Can this be problematic?**\n\nAs mentioned in Section 3.1 Line 156 of our manuscript, ``Without loss of generality, we assume the origin to be the equilibrium point.'' For a general nonlinear system with equilibrium not at the origin or more than one equilibrium, our approach applies in the following way. \n\nSuppose $\\dot{x}=f(x),x\\in\\mathbb{R}^n$ has two equilibrium points, $x_1^*,x_2^*\\in\\mathbb{R}^n,f(x_1^*)=f(x_2^*)=0$. For each equilibrium, define a change of variable. First for equilibrium point $x_1^*$, define $\\tilde{x}_1=x-x_1^*$. The dynamics of the transformed coordinate  $\\dot{\\tilde{x}}_1=f(\\tilde{x}_1+x_1^*):=f_1({\\tilde{x}}_1)$ has an equilibrium at the origin. We can apply the proposed method to find a local Lyapunov function for $f_1$ near the origin, which corresponds to a local Lyapunov for $f(x)$ near equilibrium point $x_1^*$. Similar procedures can be applied to $x_2^*$ to obtain a local Lyapunov for $f(x)$ near equilibrium point $x_2^*$.\n\nIn summary, our method can apply to general nonlinear systems. Users need to first compute the equilibrium points of the original nonlinear system. For each equilibrium, define the change of variables and transform system dynamics to have equilibrium point at origin. Then, users apply the proposed approach to find local Lyapunov function for each of the transformed system separately, which correspond to different local Lyapunov functions near each equilibrium. We will include the above discussion in the final version.  \n\nWe hope the above explanation clarifies this question. Should there be any additional question, we would be happy to address them in the discussion period.\n\n**Q3: Did you find some systems where the algorithm cannot find a Lyapunov candidate?**\n\nThanks for the valuable question. Our framework will fail to return a valid Lyapunov function if the dynamical system is unstable. Also, the framework may encounter difficulties in highly complex dynamics. For example, the success rate for high-dimensional systems is not consistently 100\\%, indicating occasional failures within a predefined number of training epochs under some random seeds. Importantly, such failures do not imply system instability. As shown in Table 1, one can successfully find valid Lyapunov functions from other random seeds (with different network initializations). In comparison, for unstable inverted pendulum without control, using the same experiment settings as in simple pendulum (Appendix G.1), our framework fails to recover an expression satisfying Lyapunov conditions over 200 epochs' training under different initializations.\n\n**Q4: To get a better sense of the proposed improvement, would it be possible to sample random dynamical systems (similar to Alfarano et al) and see how many you can recover compared to Alfarano, ANLC and FOSSIL 2.0?**\n\nThanks for the constructive suggestion. To complete the comparison, we contacted the authors of Alfarano et al. (2024), who conducted the evaluation of their pre-trained model on our test systems. Due to constraints related to global stability and the dimensionality of their training data, their model can only successfully find Lyapunov functions for systems in Appendix F.1, F.2, F.3, \\& G.1. These systems have dimensions lower than 6 and the found Lyapunov functions offer global stability guarantees. We will include these results in the experiment section of the final paper.", "claims_and_evidence": "The new framework can effectively discover valid analytical Lyapunov. Authors support the claim by considering different Lyapunov systems where a solution was not known in the past. \nTable 1 provides data on runtime, verification time, discovered Lyapunov functions, stability type (local or global asymptotic stability), and success rates for various systems.", "methods_and_evaluation_criteria": "Adding the symbolic component was really missing from previous line of work where they parametrized Lyapunov function as a non-linear neural network.", "theoretical_claims": "There is no particular theoretical claim made by authors.", "experimental_designs_or_analyses": "Authors compared the new framework with standard baselines (ANLC and FOSSIL 2.0).", "relation_to_broader_scientific_literature": "The current scientific literature relies on either parametrized model (Lyapunov candidate) and use backpropagation + falsification to find a local Lyapunov function or train a general model and ask the model to guess a global Lyapunov candidate.\nThe key contributions is the symbolic parametrization and the use of RL to successfully train the model.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "N/A", "questions_for_authors": "1. Given the symbolic answer provided by the model, can this help to enlarge the stability area given you are not overfitting on a particular locality domain? \n2. It's not clear how much general can the approach be: due to the nature of symbolic representation, this forces the candidate to be always the same function in the entire domain. Can this be problematic? \n3. Did you find some systems where the algorithm cannot find a Lyapunov candidate? \n4. To get a better sense of the proposed improvement, would it be possible to sample random dynamical systems (similar to Alfarano et al) and see how many you can recover compared to Alfarano, ANLC and FOSSIL 2.0?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "i11TBElwV4", "reviewer": "Reviewer_s9Wv", "summary": "This paper proposes to use transformers and reinforcement learning (RL) to construct local analytical Lyapunov functions for high-dimensional non-polynomial systems. The proposed framework consists of three components: 1) a symbolic transformer to generate candidate Lyapunov functions; 2) a numerical verifier for finding counterexamples; and 3) a risk-seeking policy gradient algorithm to optimize the symbolic transformer based on the Lyapunov risk. Furthermore, to improve the training efficiency, a genetic programming component is also included in the training framework. Experimental results on dynamical systems with up to 10 dimensions show that the proposed algorithm successfully found the analytical Lyapunov function in most cases, and even discovered new Lyapunov functions for the 2-bus lossy system.\n\n## update after rebuttal\n\nThe rebuttal addresses my concerns, and I think the paper can benefit from the additional results provided in the rebuttal. I have also read other reviewers' comments and the authors' responses. The comments on the scalability of the SOS methods raised by Reviewer uz6L were also my concerns on the additional experimental results, but I think the authors' reply is fair. Overall, I think this is a good paper that provides valid insights on learning-based Lyapunov function discovery. I have increased my score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer s9Wv for the valuable time and constructive feedback. We provide point-by-point responses below.\n\n**Q1: It would be great if more baselines for analytical Lyapunov functions were introduced, e.g., sum-of-squares methods (at least for the polynomial systems).**\n\nThanks for the valuable suggestion. We have added comparison to sum-of-squares methods as a baseline for both polynomial and non-polynomial systems detailed as below. \n\nPolynomial System: For a given $n$-dimensional polynomial (poly) dynamics $f(x)$ and a degree 2$d$, to check the globally asymptotical stability of $f(x)$, SOS method aims to find a poly $V(x)$ of degree $2d$, such that 1) $V(x)-\\sum_{i=1}^{n}\\sum_{j=1}^{d}\\epsilon_{ij}x_i^{2j}$ is SOS where $\\sum_{j=1}^{d}\\epsilon_{ij}>\\gamma,\\forall i=1,...,n$ with $\\gamma>0$, and $\\epsilon_{ij}\\geq0$ for all $i$ and $j$, and  2) $-\\frac{\\partial V}{\\partial x}f(x)$ is SOS [1]. \n\nFor local stability analysis, consider a ball of radius $r$ centered at origin $\\mathcal{B}_r(0)$, which can be represented by the semialgebraic set \n\n$S=\\\\{x:g(x,r)\\geq0$, where $g(x,r)=r-\\sum_{i=1}^{n}x_i^2\\\\}$. \n\nWe require that the stability condition holds in $S$. Retaining the same optimization objective and constraints on $V(x)$ as before, a modified constraint on Lie derivative is imposed: $-\\frac{\\partial V}{\\partial x}f(x)-s(x)g(x,r)$ is SOS for some SOS poly $s(x)$. If such an $s(x)$ exists, we can certify local stability.\n\nWe develop our code based on *findlyap* function from SOSTOOLS (MATLAB) and issue-16 of SOSTOOLS' official GitHub repo to examine the SOS method on poly systems in Appendix F. Table summarizes the results, and SOSTOOLS identifies valid Lyapunov functions for systems in Appendix F.1, F.2, and F.3 (up to 3-D). For example, $V(x) =0.86x_1^2-0.60x_1x_2-6.63e^{-7}x_1x_3+0.90x_2^2+4.49e^{-7}x_2x_3+0.79x_3^2$ certifies the global stability of 3-D systems in Appendix F.2.\n\n## Table: Performance of SOS Method on Poly Systems\n| **System**| **Degree 2d** | **Runtime** | **Region**|\n|---|---|--|---|\n| App. F.1| $2$|$0.697$s|$\\mathcal{B}_1(0)$|\n| App. F.2| $2$|$0.832$s|Global|\n| App. F.3 - I|$2$|$0.497$s|Global|\n| App. F.3 - II|$4$|$2.509$s|Global|\n\nHowever, for higher-dimensional systems considered in Appendix F.4, F.5,\\& F.6, SOSTOOLS fails to certify global stability - as expected, since these systems are inherently not globally stable. Furthermore, it is unable to verify local stability due to computational limitations arising from the complexity of the $s(x)g(x,r)$ term in the Lie derivative constraint, which grows quadratically with system dimension. For system in Appendix F.4, we have $g(x,1)=1-\\sum_{i=1}^6x_i^2$ containing 7 terms, and degree 2 SOS poly $s(x)$ has 36 terms at most. Thus, the constraint on lie derivative would have hundreds of terms, which is unsolvable by SOSTOOLS. \n\nPlease refer to Response to Reviewer uz6L (Q1) for SOS method on non-poly systems.\n\n**Q2: It would be more interesting to see if the ground truth Lyapunov function is not a simple form like $\\sum_{i=1}^nx_i^2$.**\n\nConsider the synthetic dynamics adapted from Appendix F.4 \\& G.2 with interactions between two subsystems:\n$$\n\\dot{x}_1=- x_1+0.5x_2-0.1x_5^2,\n$$\n$$\n\\dot{x}_2=-0.5x_1-x_2+0.1x_8,\n$$\n$$\n\\dot{x}_3=-x_3+0.5x_4-0.1x_1^2,\n$$\n$$\n\\dot{x}_4=-0.5x_3-x_4,\n$$\n$$\n\\dot{x}_5=-x_5+0.5x_6,\n$$\n$$\n\\dot{x}_6=-0.5x_5-x_6+0.1x_2^2,\n$$\n$$\n\\dot{x}_7=x_8,\n$$\n$$\n\\dot{x}_8=-\\sin(x_7)\\cos(x_7)-x_8-\\sin(x_9)\\cos(x_9)-0.1x_2,\n$$\n$$\n\\dot{x}_9=x_8-x_9.\n$$\n\nTo properly address the trigonometric terms in $\\dot{x_8}$, the Lyapunov function for this dynamics can't be a simple form like $\\sum_{i=1}^{n}x_i^2$ and should include some trigonometric terms. Setting the state space $\\mathcal{D}=\\\\{x\\in\\mathbb{R}^9| |x_i|\\leq1.5,\\forall i=1,\\cdots,9\\\\}$, our method successfully identifies a valid Lyapunov function $V=\\sum_{i=1}^6x_i^2+\\sin(x_7)^2+x_8^2-\\cos(x_9)+1$, which passes formal verification following settings in Section 5. This example and power system examples in Appendix G.3 \\& G.5 demonstrate our method's capability to recover the Lyapunov function with complex structures. \n\n**Q3: Messy citation format.**\n\nWe thank the reviewer for the helpful suggestion. We will solve the format issue in the final version.\n\n**Q4: In Figure 2, where does the $x_1$ in the bottom-right green box come from?**\n\nConsider the binary tree in the top-right corner, when generating the last token $x_2$, its parent is '+', and its sibling is $x_1$. As a result, the $x_1$ shows up in the bottom-right green box.\n\n**Q5: Can the authors discuss the completeness of the proposed framework? What if the algorithm is applied to an unstable system?**\n\nPlease refer to response to reviewer uz6L (Q1) for the completeness of proposed method and reviewer ve2p (Q3) for the analysis on unstable systems.\n\n[1] A. Papachristodoulou, and S. Prajna. A tutorial on sum of squares techniques for systems analysis. In Proceedings of the IEEE American Control Conference, 2005.", "claims_and_evidence": "The claims made in the submission are supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "The proposed methods and/or evaluation criteria (e.g., benchmark datasets) make sense for the problem or application at hand.", "theoretical_claims": "I checked the correctness of any proofs for theoretical claims. They are standard results from the literature, and there are no issues.", "experimental_designs_or_analyses": "I checked the soundness/validity of the experimental designs or analyses. Below are some potential improvements that could be made. \n\n1. The baselines of the experiments can be improved. The proposed framework specifically targets **analytical** Lyapunov functions, but the introduced baselines are for finding **neural** Lyapunov functions. It would be great if more baselines for analytical Lyapunov functions were introduced, e.g., sum-of-squares methods (at least for the polynomial systems). \n\n2. For the experiment to demonstrate the scalability of the proposed framework on the 10-D system, it would be more interesting to see if the ground truth Lyapunov function is not a simple form like $\\sum_{i=1}^{10}x_i^2$ since this is generally the first analytical form that a human would try. Potentially, the authors could change the variables a bit to design a more difficult ground truth Lyapunov function, which would introduce more impressive results.", "relation_to_broader_scientific_literature": "The key contributions of the paper are strongly related to the construction of Lyapunov functions using machine learning techniques. This is crucial for verifying the stability of dynamical systems and developing stable controllers for dynamical systems.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "**Other Strengths:**\n\n1. The writing of the paper is clear and easy to follow. \n\n2. All parts of the proposed framework are ablated in the experiments.\n\n3. The newly discovered Lyapunov function part is impressive.", "other_comments_or_suggestions": "I suggest the authors to use `\\citep` and `\\citet` carefully. The current citation format is messy.", "questions_for_authors": "1. In Figure 2, where does the $x_1$ in the bottom-right green box (the input to the encoder) come from?\n\n2. Can the authors discuss the completeness of the proposed framework? What if the algorithm is applied to an unstable system?\n\n3. Is it possible to add more experiments considering my comments in **Experimental Designs Or Analyses**?\n\nI will be willing to increase the score if the questions are addressed.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents a novel reinforcement learning (RL) framework that uses a symbolic transformer to discover analytical Lyapunov functions for nonlinear dynamical systems. The approach combines a risk-seeking policy gradient, a numerical verifier based on global optimization, and optional guidance from genetic programming. Reviewers found the paper potentially impactful, particularly due to its symbolic formulation, ability to handle non-polynomial systems, and success in discovering Lyapunov functions for challenging benchmarks such as the 2-bus lossy power network.

A key strength highlighted by reviewers is the use of a symbolic transformer, which contrasts favorably with prior approaches that rely on parameterized neural models. The method demonstrates strong empirical performance across a diverse suite of benchmarks, notably solving cases that were previously unsolved and overcoming limitations of dimensionality and polynomial constraints found in prior work.

Before the rebuttal, reviewers raised concerns about the lack of direct comparisons with established SDP/SOS-based techniques, which can scale up to 10-dimensional systems. The rebuttal partially addressed these concerns, but a more direct empirical comparison should be included in the camera-ready version. Additionally, several reviewers suggested that the paper would benefit from a theoretical or statistical analysis of the symbolic transformer’s behavior—e.g., convergence guarantees or completeness. The authors should clearly address this in the revision. Further motivation or ablation to justify the design choice of soft verification would also strengthen the paper.

In summary, this work introduces a promising new direction for learning-based discovery of analytical Lyapunov functions, effectively bridging symbolic synthesis and reinforcement learning. The rebuttal has helped alleviate major concerns raised by reviewers. The meta reviewer recommend acceptance, contingent on the authors incorporating the suggested additional comparisons and clarifications in the final version.
Decision: Accept (poster)