Analytical Lyapunov Function Discovery: An RL-based Generative Approach

We thank reviewer s9Wv for the valuable time and constructive feedback. We provide point-by-point responses below.

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).

Thanks 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.

Polynomial 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].

For local stability analysis, consider a ball of radius $r$ centered at origin $\mathcal{B}_r(0)$, which can be represented by the semialgebraic set

S=\\{x:g(x,r)\geq0, where g(x,r)=r-\sum_{i=1}^{n}x_i^2\\}.

We 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.

We 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.

Table: Performance of SOS Method on Poly Systems

However, 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.

Please refer to Response to Reviewer uz6L (Q1) for SOS method on non-poly systems.

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$.

Consider the synthetic dynamics adapted from Appendix F.4 & G.2 with interactions between two subsystems:

$$\dot{x}_1=- x_1+0.5x_2-0.1x_5^2,$$ $$\dot{x}_2=-0.5x_1-x_2+0.1x_8,$$ $$\dot{x}_3=-x_3+0.5x_4-0.1x_1^2,$$ $$\dot{x}_4=-0.5x_3-x_4,$$ $$\dot{x}_5=-x_5+0.5x_6,$$ $$\dot{x}_6=-0.5x_5-x_6+0.1x_2^2,$$ $$\dot{x}_7=x_8,$$ $$\dot{x}_8=-\sin(x_7)\cos(x_7)-x_8-\sin(x_9)\cos(x_9)-0.1x_2,$$ $$\dot{x}_9=x_8-x_9.$$

To 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.

Q3: Messy citation format.

We thank the reviewer for the helpful suggestion. We will solve the format issue in the final version.

Q4: In Figure 2, where does the $x_1$ in the bottom-right green box come from?

Consider 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.

Q5: Can the authors discuss the completeness of the proposed framework? What if the algorithm is applied to an unstable system?

Please refer to response to reviewer uz6L (Q1) for the completeness of proposed method and reviewer ve2p (Q3) for the analysis on unstable systems.

[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.

We thank reviewer ve2p for the valuable time and constructive feedback. We provide point-by-point responses below.

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?

Thanks for the insightful question.

Firstly, 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.

Furthermore, 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.

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?

As 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.

Suppose $\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^*$.

In 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.

We hope the above explanation clarifies this question. Should there be any additional question, we would be happy to address them in the discussion period.

Q3: Did you find some systems where the algorithm cannot find a Lyapunov candidate?

Thanks 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.

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?

Thanks 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.

Q1 and Q2: Justify the claimed advantage over SOS techniques. Sensitivity issues of SOS. Completeness of their approach.

Thanks 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.

Polynomial 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.

Non-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.

Table: Performance of SOS Method on Non-Poly System

To 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,$ where $F,G_1,G_2$ are vectors of functions.

Define $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.

Let $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. $V(0,x_{2, 0})=0,\\;(1)$ $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)$ $-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)$ where $\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.

Compared 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.

As 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.

Sensitivity 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.

Q3: soft verification

We 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.

Q4: Discussion of [2]

This 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.

[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.

We thank reviewer SwEU for the valuable time and constructive feedback. We provide point-by-point responses below.

Q1: Alfarano et al. (2024) was excluded from direct comparison

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 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.

Q2: No comparison with traditional symbolic regression or optimization-based methods.

1. 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.

2. 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.

3. 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.

We will highlight the comparisons to classic symbolic regression methods (GP) and add the comparison to optimization-based SOS method in the final version.

W1: It's difficult to understand which part provides the most gains.

Thanks 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.

Q4: Exploration-exploitation trade-off in your approach:

Thanks 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.

As 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.

For 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.