LyZNet: A Lightweight Python Tool for Learning and Verifying Neural Lyapunov Functions and Regions of Attraction

The proposed tool solves Zubov’s equation using physics-informed neural networks (PINNs) (Lagaris et al., 1998; Raissi et al., 2019) and verifies the neural solutions using satisfiability modulo theories (SMT) solvers. The tool currently relies on dReal (Gao et al., 2013) for verification through interval analysis.

Computation of Lyapunov functions using Zubov’s equation has been previously explored, for example, with radial basis functions (Giesl, 2007) and SOS techniques (Jones and Peet, 2021) (though the authors did not explicitly mention Zubov’s equation). More recently, the authors of (Kang et al., 2023) employed a data-driven approach to solve Zubov’s equation using neural networks; however, the computed neural solutions were not verified, and Zubov’s equation was not directly encoded in the loss function. The authors of (Liu et al., 2023b) solved Zubov’s equation using neural networks and verified them with dReal, but did not contrast their results with SOS techniques.

Many authors have recently investigated the use of neural networks for computing Lyapunov functions (see, e.g., (Chang et al., 2019; Grüne, 2021b; Gaby et al., 2022; Abate et al., 2020; Gaby et al., 2022; Kang et al., 2023), and (Dawson et al., 2023) for a recent survey). In fact, such efforts date back to as early as the 1990s (Long and Bayoumi, 1993; Prokhorov, 1994). Unlike SDP-based synthesis of SOS Lyapunov functions, neural network Lyapunov functions obtained by training are not guaranteed to be Lyapunov functions. Subsequent verification is required, e.g., using satisfiability modulo theories (SMT) solvers (Chang et al., 2019; Ahmed et al., 2020). The use of SMT solvers for searching and refining Lyapunov functions has been explored previously (Kapinski et al., 2014). Counterexample-guided search of Lyapunov functions using SMT solvers is investigated in (Ahmed et al., 2020) and the associated tool (Abate et al., 2021), which supports both Z3 (De Moura and Bjørner, 2008) and dReal (Gao et al., 2013) as verifiers. Neural Lyapunov functions with SMT verification are explored in (Zhou et al., 2022) for systems with unknown dynamics. SMT verification is often time-consuming, especially when seeking a maximal Lyapunov function (Liu et al., 2023b) or dealing with high-dimensional systems. Recent work has also focused on learning neural Lyapunov functions and verifying them through optimization-based techniques, e.g., (Chen et al., 2021a, b; Dai et al., 2021, 2020). Such techniques usually employ (leaky) ReLU networks and use mixed integer linear/quadratic programming (MILP/MIQP) for verification.

In contrast with existing approaches, the proposed tool offers the ability to compute regions of attraction (ROA) close to the domain of attraction and provides support for compositional verification to handle high-dimensional systems. We demonstrate its usage and effectiveness on both low-dimensional systems with non-trivial dynamics and high-dimensional systems to showcase compositional verification. The tool achieves state-of-the-art results on ROA approximations for nonlinear systems and their verification for high-dimensional systems.

The tool consists of a set of Python modules built around the machine learning framework PyTorch (Paszke et al., 2019) and the numerical satisfiability modulo theories (SMT) solver dReal (Gao et al., 2013). It also relies on the NumPy library (Harris et al., 2020) for numerical computation, SciPy library (Virtanen et al., 2020) for scientific computing tasks, and the SymPy library (Meurer et al., 2017) for symbolic calculation.

We describe the main modules as follows:

• •

  dynamical_systems.py defines the DynamicalSystems class and associated methods such as linearization, computation of quadratic Lyapunov function by solving Lyapunov’s equation, and providing interfaces for numerical and symbolic function valuations related to the vector field.

• •

  local_verifier.py provides functions for verifying local stability using quadratic Lyapunov functions, with interval analysis from the SMT solver dReal. The key component is using the mean value theorem to bound the remainder of the linear approximation. This enables exact verification of local asymptotic stability, despite dReal’s numerical conservativeness. This sets it apart from existing methods that often ignore verification around a neighbourhood of the origin.

• •

  quadratic_verifier.py performs reachability analysis using a quadratic Lyapunov function. The target is the local ROA obtained above. The set guaranteed to reach this target is formulated as a sublevel set of the quadratic Lyapunov function.

• •

  generate_data.py generates data for training neural Lyapunov functions. It involves solving system (1), augmented with an additional state to compute value for the solution to Zubov’s equation, which we accomplish using SciPy and “embarrassing parallelization” with the joblib library (Joblib Development Team, 2023).

• •

  neural_learner.py learns a neural Lyapunov function by solving Zubov’s PDE using physics-informed neural networks (PINN) (Lagaris et al., 1998; Raissi et al., 2019). We allow customization of neural network structure as well as selection of different loss function modes.

• •

  neural_verifier.py verifies the learned neural Lyapunov function indeed satisfies the Lyapunov condition required for reach a target invariant set within the domain of attraction. The target can be set as the region of attraction already verified by local analysis or quadratic Lyapunov analysis above. Our current implementation uses dReal as the verifier.

• •

  network_verifier.py includes functions for compositional verification of local stability using quadratic and neural Lyapunov functions. The underlying theory involves vector Lyapunov functions, differential inequalities, and comparison techniques for nonlinear systems.

• •

  sos_verifier.py, plot.py, and utils.py provide support for comparisons with other Lyapunov functions, such as sums-of-squares (SOS) Lyapunov functions, visualization of verified level sets, and other routines used in learning and verification.

Using the tool is straightforward. As shown in the code snippet below, we can define a dynamical system with SymPy variables and symbolic expressions for the vector field. We can also specify the domain for learning and training, as well as name the system for file management and result storage.

The code above defines the reversed Van der Pol equation with parameter $\mu=1.0$ on the domain $X=[-2.5,2.5]\times[-3.5,3.5]$.

The following lines of code can be used to call relevant functions for verifying local stability using quadratic Lyapunov function, as well as learning and verifying neural Lyapunov functions.

Most functions and parameters are self-explanatory. We remark that the argument loss_mode in neural_learner currently can be set to three different modes: "Zubov", "Data", and "Lyapunov". The latter two take a purely data-driven approach for solving Zubov’s equation or encode a Lyapunov inequality instead of Zubov’s equation, respectively.

The code above computes and verifies four sublevel sets:

satisfying $\mathcal{P}_{1}\subseteq\mathcal{P}_{2}$, $\mathcal{V}_{1}\subseteq\mathcal{P}_{2}$, and $\mathcal{V}_{1}\subseteq\mathcal{V}_{2}$, where $V_{P}$ is a quadratic Lyapunov function of the form $V_{P}(x)=x^{T}Px$ and $V_{N}(x)$ is a learned neural network Lyapunov function. The set $\mathcal{P}_{1}$ is a local region of attraction verified by linearization. The set $\mathcal{P}_{2}$ represents the set of initial conditions from which solutions of (1) are verified to reach $\mathcal{P}_{1}$ using a quadratic Lyapunov function. The sets $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ are verified using a neural Lyapunov function, where $\mathcal{V}_{1}$ represents a target invariant set contained in $\mathcal{P}_{2}$¹¹1This could be set to $\mathcal{P}_{1}$ if one wishes to skip the call of quadratic_reach_verifier for any reason, or any known region of attraction verified by other means., and $\mathcal{V}_{2}$ is an invariant set from which solutions of (1) are guaranteed to reach $\mathcal{V}_{1}$. All sets are also verified to be contained in the interior of $X$ to ensure set invariance. Put together, $\mathcal{V}_{2}$ is a verified ROA for (1).

Local stability analysis using quadratic Lyapunov functions is fairly standard. Suppose that $x=0$ is an exponentially stable equilibrium point for (1); i.e., $f(0)=0$ and the Jacobian matrix $A=Df(0)=\frac{df}{dx}\big{|}_{x=0}$ is a Hurwitz matrix, i.e., all eigenvalues of $A$ have negative real parts. Rewrite (1) as

where $g(x)=f(x)-Ax$ satisfies $\lim_{x\rightarrow 0}\frac{\left\|g(x)\right\|}{\left\|x\right\|}=0$. Given any symmetric positive definite real matrix $Q\in\mathbb{R}^{n\times n}$, there exists a unique solution $P$ to the Lyapunov equation

Let $V_{P}(x)=x^{T}Px$. Then

where we set $r=\lambda_{\min}(Q)-\varepsilon>0$ for some sufficiently small $\varepsilon>0$ and $\lambda_{\min}(Q)$ is the minimum eigenvalue of $Q$. The goal of local_stability_verifier is to determine the largest number $c1\_P$ such that

where $\mathcal{P}_{1}$ is as defind in (3).

While (6) can be easily formulated as a satisfiability condition in dReal (Gao et al., 2013), due to conservative use of interval analysis in dReal to account for numerical errors, verification of inequalities such as $2x^{T}Pg(x)\leq r\left\|x\right\|^{2}$ will return a counterexample close to the origin. To overcome this issue, we look at a higher-order approximation of $Pg(x)$. By the mean value theorem, we have

where $Dg$ is the Jacobian of $g$ given by $Dg=Df-A$, which implies

As a result, to verify (6), we just need to verify

If $\mathcal{P}_{1}$ is star-shaped with respect to the origin²²2A set $S$ is star-shaped with respect to the point $x_{0}$ if for all $y\in S$ the line segment joining $y$ and $x_{0}$ is a subset of $S$., then (7) is equivalent to

Since $Dg(0)=0$ and $Dg$ is continuous, one can always choose ${c1}\_P>0$ sufficiently small such that (8) can be verified. Furthermore, in rare situations, if (8) can be verified for $X=\mathbb{R}^{n}$, then global exponential stability of the origin is proved for (1).

The reversed Van der Pol equation is given by

where $\mu>0$ is a parameter that affects the stiffness of the equation.

Consider the two-machine power system (Vannelli and Vidyasagar, 1985) modelled by

where $\delta=\frac{\pi}{3}$.

Consider an inverted pendulum

where the linear gains are given by $k_{1}=[4.4142,2.3163]$.

Consider the 10-dimensional nonlinear system from (Grüne, 2021a):

We report the computational times required for the training and verification processes in Examples 5.1 and 5.2 on a 2020 MacBook Pro with a 2 GHz Quad-Core Intel Core i5, without any GPU. For training, we randomly chose 300,000 collocation points in the domain and trained for 20 epochs. The sizes of the neural networks were as reported in Tables 1 and 2. Verification in dReal can be easily parallelized by adjusting the config parameter number_of_jobs, which may correspond to the number of cores/threads on the computer. For the laptop used in these experiments, there are four cores, allowing for eight threads via hyperthreading. As discussed in Remark 3, computational time is largely dependent on computer hardware. The computational times provided in Table 4 serve as a reference point. We expect that training time can be significantly improved by using GPUs, and verification time can be further reduced by utilizing additional cores.