{"id":"4d616b67-afaa-4cf9-af03-6765fd8a12c1","arxiv_id":"2502.02014","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"high","formal_verification":"none","parameter_count":7,"one_line_summary":"A reinforcement-learning-driven symbolic transformer generates and verifies analytical local Lyapunov functions for nonlinear systems up to ten dimensions, including a claimed new certificate for a lossy power system.","lead":"The paper trains a transformer with reinforcement learning to generate symbolic, human-readable Lyapunov functions that certify stability of nonlinear systems, then verifies them with a formal solver. It reports success on systems from two to ten dimensions, including a claimed first Lyapunov certificate for a lossy power system.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Formal verification domain for the 2-bus lossy system omits the center-of-inertia constraint delta1+delta2=0; on the stated box V1=0 at delta1=delta2≠0, so the 'previously unknown Lyapunov function' claim is not supported as written.","rationale":"The central claim is that the transformer-RL pipeline discovers a valid, previously unknown analytical Lyapunov certificate for a 2-bus lossy power system, formally verified by dReal. For that claim to hold, the certificate must be positive definite and have negative Lie derivative on the exact domain passed to the verifier. The manuscript defines the domain as a box in R^4, while its own coordinate change in G.5 makes delta1+delta2=0. On the box, the displayed V1 is not positive definite: with omega=0 and delta1=delta2 any nonzero value in [-0.75,0.75], V1=0, and these points are outside the excluded epsilon-ball. A dReal query over D cannot certify such a function. The reader's weakest assumption is exactly this hidden-domain issue, and I agree. This is a concrete correctness/reproducibility gap in the headline result, not a stylistic disagreement with the control community. It can be repaired: on the manifold, V1 becomes omega1^2+omega2^2+(omega2−2 sin(delta))^2 and may well be valid; the method itself is plausible and the other experiments (pendulum, 3-D trig, polynomial systems) provide independent support for the pipeline. The conditional verdict remains appropriate: the paper should state the verification domain, encode the sum-to-zero constraint, and re-run dReal, or weaken the 'previously unknown' claim.","tokens_in":29667,"tokens_out":8139,"duration_ms":80678,"concrete_test":"Use dReal with the paper's stated delta=1e−12 and epsilon=1e−3 on D exactly as written for V1: check whether there exists x in D\\B_epsilon(0) with V1(x)≤1e−12 or LfV1(x)≥−1e−12; the point (delta1,delta2,omega1,omega2)=(0.5,0.5,0,0) is an immediate counterexample. Then repeat with the manifold constraint delta1+delta2=0 encoded explicitly (equivalently delta2=−delta1) and with the 2-bus swing equations rewritten in the reduced coordinate delta=(theta1−theta2)/2, verifying both V1 and V2 on the 3D domain. If the constrained queries pass and the unconstrained queries fail, the concern is resolved but the paper must state the manifold domain.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 5.3 states D = {|delta_i|≤0.75, |omega_i|≤2} ⊂ R^4 and reports both V1 and V2 pass dReal. But G.5 defines delta_i = theta_i − (1/N)Σtheta_i, so for N=2, delta1+delta2=0 identically. The reported D is a box, not the 3D manifold, and V1 = omega1^2+omega2^2+(omega2−sin(delta1)+sin(delta2))^2 vanishes at delta1=delta2=c≠0, omega1=omega2=0, a point outside B_epsilon(0). Thus V1 is not positive definite on D, and any dReal query over the stated box should find this counterexample instead of certifying V1. If the verifier implicitly added delta1+delta2=0, the paper never states it; if the domain is the manifold, the free direction in the box is unphysical and the certificate is not for the stated region. The same pattern appears in G.3: the lossless 6-D function is negative at delta1=delta2=delta3=c≠0 unless sum-to-zero is enforced. This hidden-domain ambiguity is load-bearing because the headline 'previously unknown Lyapunov function for lossy power systems' depends entirely on the certificate being valid on the verified domain.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes an end-to-end framework that uses a symbolic transformer trained with risk-seeking policy gradient and genetic-programming expert guidance to generate candidate analytical Lyapunov functions, with SHGO-based global optimization for counterexample feedback during training and dReal for final verification. The authors report success on polynomial and non-polynomial systems up to 10 dimensions, and claim to have discovered previously unknown local Lyapunov functions for a 2-bus lossy power system. The manuscript also compares against neural Lyapunov baselines, SOS methods, and the pre-trained transformer of Alfarano et al. (2024), and provides ablations of the risk-seeking quantile, verification strategy, and GP guidance.","tokens_in":29994,"tokens_out":7016,"duration_ms":71935,"significance":"If the results are correct, the framework is a useful contribution: it produces concise, interpretable, analytical certificates, trains from scratch without a large pre-training dataset, scales to 10-dimensional dynamics, and ships reproducible code. The ablations are informative, and the use of an external SMT solver for final verification is a genuine safeguard against reward hacking. However, the correctness of the headline claims depends on the reported Lyapunov functions being certified on the stated state-space domains, and several of the displayed certificates appear not to satisfy the stated Lyapunov conditions on those domains. These verification/domain mismatches are load-bearing for the central contribution, so the manuscript needs substantial revision before acceptance.","major_comments":[{"comment":"The claimed 'previously unknown' Lyapunov functions for the 2-bus lossy power system are not validated as positive definite on the stated domain. The paper defines D = {(δ1,δ2,ω1,ω2) ∈ R^4 : |δi| ≤ 0.75, |ωi| ≤ 2}, but Appendix G.5 defines δi = θi − (1/N)Σθj, which for N=2 forces δ1+δ2=0. The function V1 = ω1^2 + ω2^2 + (ω2 − sin δ1 + sin δ2)^2 vanishes at δ1=δ2=c≠0, ω1=ω2=0, a point in the box D that is not in a small neighborhood of the origin. Hence V1 is not positive definite on D. If the dReal query included the implicit constraint δ1+δ2=0, that constraint is never stated; if it did not, the dReal verification cannot have succeeded as reported. The same hidden-domain issue affects the 6-D lossless power-system certificate in Appendix G.3. This must be corrected and the actual verified domain stated explicitly before the central claim can be assessed.","section":"§5.3, Appendix G.5"},{"comment":"The Van der Pol certificate V = x1^2 + x2^2 does not satisfy the strict Lyapunov inequality required by Proposition 3.4 on the stated domain D = {|x1|≤1, |x2|≤1}. Its Lie derivative is LfV = −2(1 − x1^2)x2^2, which vanishes on the entire line x2=0 and on the boundary x1=±1, not merely at the origin. Since the reported dReal verification uses strict inequalities V(x)>δ and LfV(x)<−δ over D\\B_ε(0), points such as (0.5,0) are counterexamples to the claimed verification. The function may still certify asymptotic stability via the invariance principle, but then the verification claim and the reported SMT certification need to be revised accordingly.","section":"Table 1, Appendix F.1"},{"comment":"The same center-of-inertia domain issue appears in the 6-D lossless power system. Appendix G.3 introduces δi = θi − (1/N)Σθj, so the physical state space is the manifold Σδi=0, but the experiment region D is written as a box in R^6 without this constraint. On the full box, the reported function V = Σωi^2 − 0.5(Σ_{i≠j} cos(δi−δj) − 1) takes negative values, e.g. at δ1=δ2=δ3=c≠0 and ω1=ω2=ω3=0. Thus the certificate is not positive definite on the stated D. Either the verification was performed on a restricted domain that must be documented, or the table overstates the validity of the certificate.","section":"Table 1, Appendix G.3"}],"minor_comments":[{"comment":"The mapping g is defined as g(x)=1/(1+x) in Equation (5), but Appendix C states 'we adopt a bounded reward function using the continuous mapping g(x)=1/x'. These are inconsistent; the intended formula is presumably 1/(1+L).","section":"Equation (5) vs. Appendix C"},{"comment":"The displayed definitions of V1 and V2 contain a typo: the intended leading term appears to be Σ_{i=1}^2 ω_i^2, but the manuscript writes 'Σ_{i=1}^2 ω1^2' for both functions.","section":"§5.3"},{"comment":"The paper should state explicitly whether the SHGO verification during training and the final dReal query used the box D or the manifold δ1+δ2=0; the current text makes this impossible to determine from the reported experiments.","section":"§5.3, Appendix G.5"},{"comment":"The table reports '−' for SOS on the 6-D, 8-D, and 10-D polynomial systems, but the text does not specify whether this is a timeout, a solver failure, or an infeasibility result. A brief note would make the comparison reproducible.","section":"§5.4, Table 3"},{"comment":"The phrase 'Ramdom sampling' should be corrected to 'Random sampling' in the table header.","section":"Appendix H.2"}],"recommendation":"major_revision","confidential_remarks":"The domain/verification mismatch in the power-system examples and the Van der Pol certificate is the core issue; it directly affects the abstract's claims of discovering previously unknown Lyapunov functions. I do not see this as a circularity problem, because the final certificates are checked by dReal, but the SMT verification cannot be trusted until the verified domain is precisely specified and the reported certificates actually satisfy the inequalities on it. The literature-support for the 'previously unknown' claim is also thin (Chiang 1989; Cui & Zhang 2022), so the authors should clarify what exactly they checked against the existing literature."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The framework is a real step forward: a symbolic transformer trained from scratch with risk-seeking RL plus SHGO falsification and GP expert guidance, producing analytical Lyapunov certificates that are cheap to verify. That integration is new compared to Alfarano et al. and Feng et al., and the ablations are careful. The authors clearly know the surrounding literature, and shipping code is a plus.\n\nThe serious problem is domain mismatch in the formal verification. For the 2-bus lossy system, Appendix G.5 defines delta_i as center-of-inertia coordinates, so delta_1 + delta_2 = 0 identically. Yet Section 5.3 sets D as a box in R^4 with no such constraint. On that box, V1 vanishes at (c,c,0,0) with c != 0, so it is not positive definite. The only way the reported dReal verification could pass is if the constraint was silently added. That is load-bearing for the 'previously unknown Lyapunov function' headline. The same pattern appears in Appendix G.3. This is fixable by stating the invariant manifold and running SMT on it, but as written the claim is unsupported.\n\nThe same mismatch shows up in the simpler systems. For Van der Pol, V = x1^2 + x2^2 has Lie derivative -2(1-x1^2)x2^2, which is zero on the boundary |x1|=1 and on the line x2=0. For the pendulum, LfV = -0.2 x2^2 is zero whenever x2=0. Strict inequality LfV < -delta on D\\B_epsilon(0) cannot hold in those places, so the SMT verification statement in Section 5 is not accurate for the displayed functions and domains. These could be re-scoped with open domains or re-verified with weak inequalities plus LaSalle, but the paper currently overclaims.\n\nOne place I disagree with the reader: the risk-seeking objective sign in Equation (6) is internally consistent. They define Jrisk as the negative conditional reward and then run gradient descent; Proposition C.1 is the gradient of that negative objective. No bug there.\n\nThere is also a sketchy claim that the 6D polynomial certificate is valid on a ball of radius sqrt(500); the cubic coupling terms dominate the quadratic decay long before that. That deserves a check.\n\nThe paper deserves a serious referee. The method is novel, the flaws are concrete and addressable, and a major revision can turn this into a solid contribution. I would engage with a revised version rather than recommend rejection.","headline":"The method is genuinely new and well-engineered, but the headline certificates don't satisfy the stated verification as written; the lossy power system claim hinges on an unstated manifold constraint.","tokens_in":30540,"tokens_out":9460,"would_cite":true,"duration_ms":87630,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["93D30","37N35","68T07"],"pacs":[],"model":"deepseek-v4-flash","headline":"A reinforcement-learning-trained transformer writes analytical Lyapunov functions for nonlinear systems up to ten dimensions.","keywords":["Lyapunov functions","reinforcement learning","symbolic transformer","formal verification","nonlinear dynamical systems","power system stability","genetic programming","symbolic regression"],"falsifier":"Run dReal with the paper's stated settings ($\\epsilon=10^{-3}$, tolerance $\\delta=10^{-12}$) over the full box $D=\\{(\\delta_1,\\delta_2,\\omega_1,\\omega_2): |\\delta_i|\\le 0.75,\\ |\\omega_i|\\le 2\\}$ and query whether $V_1(\\delta_1,\\delta_2,\\omega_1,\\omega_2)\\le \\delta$ or $L_fV_1(\\cdot)\\ge -\\delta$ at some point with $\\|x\\|>\\epsilon$. At $(\\delta_1,\\delta_2,\\omega_1,\\omega_2)=(0.5,0.5,0,0)$, $V_1=0$, so if the solver returns this as a counterexample the 'verified on $D$' statement is false; if no counterexample is found, the verifier must be restricting to the sum-to-zero manifold, and that restriction should be stated as part of the certificate.","tokens_in":29418,"feed_emoji":"⚡","tokens_out":9710,"duration_ms":84408,"temperature":0.7,"pith_summary":"The paper tries to show that a generative transformer, trained from scratch on a single system via reinforcement learning, can produce short symbolic Lyapunov functions for nonlinear dynamics, including non-polynomial and high-dimensional cases where neural-network certificates are hard to verify and sum-of-squares methods struggle. Because the candidates are analytical expressions rather than neural networks, the framework can check them with global optimization during training and certify them with an SMT solver in milliseconds. If true, this turns Lyapunov-function design from a manual art into an automatable search, and it supplies interpretable stability certificates that engineers can inspect. The headline application is a 2-bus lossy power system, for which the paper claims the first analytical local Lyapunov functions in the control literature.","feed_headline":"RL pipeline finds Lyapunov functions for 10-D systems","feed_subtitle":"Symbolic certificates verify in milliseconds and reveal a previously unknown stability proof for lossy power grids.","key_machinery":"The carrying object is a symbolic transformer that represents both dynamics and candidate Lyapunov functions as pre-order traversals of expression trees; a decoder samples candidates token by token, conditioned on the dynamics and on parent–sibling tree context. Because candidates are short analytic expressions, SHGO can be run over the state space to find minimizers of $\\tilde{V}$ and of $-L_f\\tilde{V}$, the natural violation points, and the resulting counterexamples are fed back into the reward set. The reward is the bounded empirical Lyapunov risk $R(\\tilde{V})=1/(1+L(\\tilde{V}))$, and the risk-seeking policy gradient with $\\alpha=0.1$ updates the transformer only on the best quantile of samples; a genetic-programming module further refines candidates into an elite set that drives a supervised guidance loss.","core_discovery":"The central claim is that an end-to-end pipeline—a symbolic transformer that tokenizes the ODEs into prefix expressions and decodes candidate certificates, a risk-seeking policy gradient that rewards only the best candidates, SHGO-based counterexample search, and genetic-programming refinement—finds valid local Lyapunov functions for systems up to ten dimensions, and that the found expressions pass formal SMT verification. On the 2-bus lossy power-system benchmark with state space $D=\\{(\\delta_1,\\delta_2,\\omega_1,\\omega_2)\\in\\mathbb{R}^4: |\\delta_i|\\le 0.75,\\ |\\omega_i|\\le 2\\}$, the paper reports $V_1=\\omega_1^2+\\omega_2^2+\\bigl(\\sin(\\delta_2)-\\sin(\\delta_1)+\\omega_2\\bigr)^2$ and $V_2=\\omega_1^2+\\omega_2^2+\\bigl(\\sin(\\delta_2)-\\sin(\\delta_1)-\\omega_1\\bigr)^2$ as the first analytic local Lyapunov functions for that system. The paper also recovers textbook energy functions for the pendulum and a known storage function for a 3-bus lossless power system, while a sum-of-squares baseline fails on the higher-dimensional local-stability cases.","pith_inferences":["One domain subtlety the paper leaves implicit: the power-system states are center-of-inertia coordinates, so the angle deviations sum to zero, while the search domain is written as a box with no such constraint; on the full box the reported $V_1$ is not positive definite, so the headline claim should be read as holding on the manifold and the verifier configuration should state that constraint exp","The framework's dependence on SHGO to find true minimizers suggests that as dimension grows, verification quality may degrade unless the sampling budget grows correspondingly; the paper's own ablation shows random sampling converges faster on 6-D but yields worse certificates, so a hybrid or adaptive budget is a natural next test.","Since dynamics constants are tokenized with rounded numerals and the symbolic library excludes general real constants, the discovered expressions are likely tailored to the exact parameter values in the benchmarks; a parameterized variant could be tested by training on a family of power-system operating points and asking whether the same template remains valid."],"forward_implications":["If the claim holds, stability certification for nonlinear control systems becomes a search problem that returns a human-readable formula rather than a black-box network, so the certificate itself can be audited and reused.","Final SMT verification of a found symbolic candidate takes milliseconds in the reported experiments, which would make formal guarantees cheap enough to run inside a training loop.","The reported success on a 10-D polynomial system and on non-polynomial dynamics suggests the approach covers systems where SOS methods need expensive recasting and where neural-Lyapunov baselines fail to converge.","The lossy power-system discovery, if valid, fills a specific gap: stability had been known by passivity, but no analytical Lyapunov function was available to certify it.","The same pipeline is proposed in the paper as extendable to other certificate functions, such as control barrier functions for safety."],"supporting_citations":[{"why":"Defines the empirical Lyapunov risk used as reward and the SMT-based verification loop that motivates the analytical-candidate design.","marker":"Chang et al. (2019)"},{"why":"Provides the risk-seeking policy-gradient objective that optimizes only the top reward quantile.","marker":"Petersen et al. (2020)"},{"why":"Supplies SHGO, the global optimizer used to locate minimizers of the candidate and its negated Lie derivative for counterexample feedback.","marker":"Endres et al. (2018)"},{"why":"Supplies dReal, the SMT solver used for final formal verification of candidate Lyapunov functions.","marker":"Gao et al. (2013)"},{"why":"The pre-trained transformer baseline for global Lyapunov functions that this paper contrasts with from-scratch local search.","marker":"Alfarano et al. (2024)"},{"why":"Provides the Lyapunov stability definitions, the asymptotic-stability proposition, and the invariance principle used throughout.","marker":"Khalil (2002)"},{"why":"Provides the genetic-programming implementation used for candidate refinement and expert guidance.","marker":"Fortin et al. (2012)"},{"why":"SOS baseline formulation for polynomial stability certificates used in comparisons.","marker":"Papachristodoulou & Prajna (2005a)"},{"why":"SOS recasting of non-polynomial dynamics used as the comparison baseline on trig and pendulum systems.","marker":"Papachristodoulou & Prajna (2005b)"}],"fun_headline_variants":["RL discovers analytical Lyapunov functions for 10D systems","Symbolic RL writes verifiable stability proofs for power grids","Transformer plus RL yields interpretable stability certificates","AI finds new analytic Lyapunov functions for nonlinear systems up to 10D","Closed-form Lyapunov functions found via RL, verified by SMT"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The claim depends on the verifier checking the Lyapunov conditions on the system's true state-space domain: the power-system examples use center-of-inertia angle coordinates whose sum is fixed at zero, while the paper states the domain as a box with no such constraint, and if the verifier silently uses the manifold the reported functions are not positive definite on the box.","fun_headline_variants_meta":{"raw":{"variants":["RL discovers analytical Lyapunov functions for 10D systems","Symbolic RL writes verifiable stability proofs for power grids","Transformer plus RL yields interpretable stability certificates","AI finds new analytic Lyapunov functions for nonlinear systems up to 10D","Closed-form Lyapunov functions found via RL, verified by SMT"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000237,"raw_usage":{"total_tokens":1565,"prompt_tokens":1059,"completion_tokens":506,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":675,"completion_tokens_details":{"reasoning_tokens":419}},"tokens_in":675,"tokens_out":506,"duration_ms":5142,"temperature":1.0,"reasoning_tokens":419,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-09T13:40:22.255952+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run dReal with the paper's stated settings ($\\epsilon=10^{-3}$, tolerance $\\delta=10^{-12}$) over the full box $D=\\{(\\delta_1,\\delta_2,\\omega_1,\\omega_2): |\\delta_i|\\le 0.75,\\ |\\omega_i|\\le 2\\}$ and query whether $V_1(\\delta_1,\\delta_2,\\omega_1,\\omega_2)\\le \\delta$ or $L_fV_1(\\cdot)\\ge -\\delta$ at some point with $\\|x\\|>\\epsilon$. At $(\\delta_1,\\delta_2,\\omega_1,\\omega_2)=(0.5,0.5,0,0)$, $V_1=0$, so if the solver returns this as a counterexample the 'verified on $D$' statement is false; if no counterexample is found, the verifier must be restricting to the sum-to-zero manifold, and that restriction should be stated as part of the certificate.","supporting_citations":[],"review_version":1}