Verifying Probabilistic Regions of Attraction with Neural Lyapunov Functions for Stochastic Systems

TL;DR

This paper introduces a neural network-based method to accurately estimate the largest probabilistic regions of attraction in stochastic systems, combining physics-informed learning with formal verification techniques.

Contribution

It develops a neural Lyapunov function approach using a stochastic Zubov's equation and SMT verification, improving over traditional methods in estimating probabilistic ROA.

Findings

More accurate probabilistic ROA estimates than SOS methods

Effective neural Lyapunov functions for stochastic stability

Validated on nonlinear stochastic systems

Abstract

Leveraging a stochastic extension of Zubov's equation, we develop a physics-informed neural network (PINN) approach for learning a neural Lyapunov function that captures the largest probabilistic region of attraction (ROA) for stochastic systems. We then provide sufficient conditions for the learned neural Lyapunov functions that can be readily verified by satisfiability modulo theories (SMT) solvers, enabling formal verification of both local stability analysis and probabilistic ROA estimates. By solving Zubov's equation for the maximal Lyapunov function, our method provides more accurate and larger probabilistic ROA estimates than traditional sum-of-squares (SOS) methods. Numerical experiments on nonlinear stochastic systems validate the effectiveness of our approach in training and verifying neural Lyapunov functions for probabilistic stability analysis and ROA estimates.

Equations

$$\{\begin{array}{l}d{X}_t=f(X_t)dt+\sigma (X_t)d{B}_t,t\ge t_0\\ X_{t_0}=x_0,\end{array}$$

$$B_t={\left(B_1(t),\dots ,B_m(t)\right)}^T$$

$$\mathbb{P}\bigl{\{}\bigl{|}X_{t}(x_{0})\bigr{|}<r\text{ for all }t\geq t_{0}\bigr{\}}\geq 1-\varepsilon,$$

$$\mathbb{P}\bigl{\{}\lim_{t\to\infty}X_{t}(x_{0})=0\bigr{\}}\geq 1-\varepsilon,$$

$$\mathbb{P}\bigl{\{}\lim_{t\to\infty}X_{t}(x_{0})=0\bigr{\}}=1.$$

$$L:=\sum _{i=1}^n{f}_i(x)\frac{\partial }{\partial x_i}+\frac{1}{2}\sum _{i,j=1}^n{\left[\sigma (x){\sigma }^T(x)\right]}_{ij}\frac{{\partial }^2}{\partial x_i\partial x_j}.$$

$$LV(x)=V_x(x)f(x)+\frac{1}{2}\mathrm{Tr}[{\sigma }^T(x)V_{xx}(x)\sigma (x)],$$

$$dV(X_t)=LV(X_t)dt+V_x(X_t)\sigma (X_t)d{B}_t.$$

$$V^c=\{x\in {\mathbb{R}}^n\mid V(x)\le c\}$$

$$p\text{-}\mathrm{ROA}:=\bigl{\{}x_{0}\in\mathbb{R}^{n}:\mathbb{P}\bigl{[}\lim_{t\to\infty}|X_{t}(x_{0})|=0\bigr{]}\geq p\bigr{\}},$$

$$\mathcal{C}=\left\{x_{0}\in\mathbb{R}^{n}:\mathbb{P}\bigl{[}\lim_{t\to\infty}|X_{t}(x_{0})|=0\bigr{]}>0\right\}.$$

$$W(x)=1-\mathbb{E}[\exp ({\int }_{t_0}^{\infty }-g(X_s(x))ds)]$$

$$\{\begin{array}{l}LW(x)=-g(x)(1-W(x)),\\ W(0)=0.\end{array}$$

$$\mathcal{C}=\{x\in {\mathbb{R}}^n:W(x)<1\}.$$

$$\mathcal{L}(\theta )=L_r(\theta )+L_b(\theta )+L_d(\theta ),$$

$$L_r(\theta )$$

$$\frac{1}{2}\mathrm{Tr}[{\sigma }^T(x_i)\frac{{\partial }^2{W}_{\mathrm{NN}}}{\partial x^2}(x_i)\sigma (x_i)]$$

$$\displaystyle\quad\quad+g(x_{i})(1-W_{\mathrm{NN}}(x_{i};\theta))\big{)}^{2}.$$

$$L_b(\theta )=W_{\mathrm{NN}}(0;\theta {)}^2.$$

$$L_d=\frac{1}{N_d}\sum _{i=1}^{N_d}(W_{\mathrm{NN}}\left(y_i;\theta \right)-\hat{W}\left(y_i\right){)}^2.$$

$$d{X}_t=\left(A{X}_t+m(X_t)\right)dt+\left(C{X}_t+n(X_t)\right)d{B}_t,$$

$$PA+A^{T}P+\sum _{i=1}^m{S}_i^{T}P{S}_i=-Q.$$

$$M(x)=D^{2}h(x)+2Q,$$

$$x\in V_P^c=\{x\in {\mathbb{R}}^n\mid V_P(x)\le c\}\Rightarrow \Vert M(x){\Vert }_F\le 2r,$$

$$L{V}_P(x)<0,\forall x\in V_P^c\setminus \{0\}.$$

$$L{V}_P(x)=$$

$$=$$

$$\mathrm{Tr}[(Cx+n(x){)}^{T}P(Cx+n(x))].$$

$$\sum _{i=1}^m{x}^T{S}_i^{T}P{S}_{i}x+2\mathrm{Tr}[n(x{)}^{T}PCx]+\mathrm{Tr}[n(x{)}^{T}Pn(x)].$$

$$L{V}_P(x)=$$

$$+\mathrm{Tr}[2n(x{)}^{T}PCx+n(x{)}^{T}Pn(x)].$$

$$L{V}_P(x)$$

$$+n(x{)}^{T}Pn(x)].$$

Full text

I Introduction

Recent advances in neural networks and machine learning have transformed computational research, including the computational aspects of control systems. For deterministic systems, the work [1] proposes a learning framework to simultaneously stabilize nonlinear systems with a linear neural controller and learn a neural Lyapunov function to certify a region of attraction (ROA). This approach was extended to cope with systems of unknown dynamics in [2]. In many applications, such as power systems [3], robotics [4], and biological systems [5], accurate estimates of the ROA are essential because they provide critical information for system stability. In [6], a maximal Lyapunov function was introduced that can exactly characterize the ROA. Recent works [7, 8, 9] have shown that neural network approximations to the solution of Zubov’s partial differential equation (PDE) can effectively capture maximal Lyapunov functions, thereby characterizing the entire ROA. Consider the potential approximation error of neural Lyapunov functions, the works [7, 8] also formally verify the satisfaction for Lyapunov condition by satisfiability modulo theories (SMT) solvers. While these techniques enable tasks such as constructing Lyapunov functions for stability analysis and control design in [1, 2, 10, 11, 12, 13, 14] within deterministic frameworks, many real-world systems are subject to random disturbances. Stochastic models can describe such systems more accurately than deterministic ones. Therefore, we are motivated to investigate probabilistic variants of ROA for stochastic systems.

In this paper, we extend the approach of using neural networks to approximate the solution of Zubov’s PDE and the SMT-based verification for neural Lyapunov functions from deterministic systems to stochastic systems. Our goal is twofold: (1) to develop neural Lyapunov functions that approximate probabilistic regions of attraction for stochastic systems, and (2) to formally verify that these neural Lyapunov functions satisfy the stochastic Lyapunov conditions, both locally and over the largest probabilistic region of attraction, i.e., the set of points which are attracted with positive probability to the equilibrium.

I-A Related work

In [15], V. I. Zubov introduced a method to characterize the region of attraction using a Lyapunov function that satisfies an appropriate first order partial differential equation, known as Zubov’s PDE. In the stochastic setting, the work [16] extends Zubov’s method to stochastic systems, defines the notion of probabilistic ROA, and characterizes probabilistic ROA using the stochastic Zubov PDE. Several follow-up works have further explored the characterization of the probabilistic ROA. For example, the work [17] builds on the theory in [16] and extends the framework to compute the almost sure ROA (i.e., attraction with probability one) for stochastic systems. Their Lyapunov function is computed by iterative numerical methods. Moreover, the work [18] employs numerical methods to compute a Lyapunov function that characterizes the probabilistic ROA for controlled stochastic systems with a prescribed probability.

In [19], although the work focuses on deriving quadratic Lyapunov functions for autonomous linear stochastic differential equations using sum-of-squares (SOS) methods, it is closely related to our study. This is because an essential step in accurately characterizing the probabilistic ROA is to compute a Lyapunov function for the linearized stochastic system. We extend their work by employing the SOS method for nonlinear stochastic systems, which serves as a comparative experiment later.

Several works [20, 21, 22] have employed neural networks to learn Lyapunov functions for stochastic systems. The work [20] proposes a learning framework for neural Lyapunov functions and a linear neural stochastic controller. In [21], neural network–based Lyapunov functions and controllers are learned using reinforcement learning methods, and almost sure asymptotic stability is formally verified for discrete-time nonlinear stochastic control systems. Meanwhile, [22] extends the methods of [21, 1] by employing a learner-verifier framework to jointly learn a control policy and a Lyapunov function. However, to the best of our knowledge, none of these studies analyze the largest probabilistic ROA for stochastic systems. This motivates us to further investigate this area. Using neural networks to solve stochastic Zubov’s equation allows us to capture the probabilistic ROA as accurately as possible. However, a drawback of neural network approximations is their lack of formal guarantees. To address this issue, SMT solvers such as [23] have been used to formally verify the trained neural Lyapunov functions and provide stability guarantees with verified ROAs in [7, 8]. In this work, we also employ an SMT solver to establish stochastic stability guarantees with formally verified probabilistic ROAs. To enable this, we derive sufficient conditions for the learned neural Lyapunov functions, which is the main contribution of this paper.

II Preliminaries

II-A Stochastic Differential Equation and Stochastic Stability

Consider a probability space $\left(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq t_{0}},\mathbb{P}\right)$, where $\Omega$ is the sample space, $\mathcal{F}$ is a $\sigma$-algebra, $\{\mathcal{F}_{t}\}_{t\geq t_{0}}$ is a filtration satisfying the usual conditions, and $\mathbb{P}$ is a probability measure. On this space, we consider the stochastic differential equation (SDE):

[TABLE]

where $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a locally Lipschitz drift term, $\sigma:\mathbb{R}^{n}\to\mathbb{R}^{n\times m}$ is the diffusion coefficient, and

[TABLE]

is an $m$-dimensional Brownian motion. We denote the solution of (1) as $X_{t}(x_{0})$ or $X(t,x_{0})$. We refer the readers to [24] for the standard conditions ensuring existence and uniqueness of solutions to equation (1). In addition, we make the following assumption.

Assumption 1

The functions $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and $\sigma:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times m}$ are twice continuously differentiable functions, $f(0)=0$, and $\sigma(0)=0$.

The following definitions are standard and can be found in [25, Section 4.2].

Definition 2** (Stochastic Stability)**

The trivial solution of (1) is said to be stochastically stable (or stable in probability) if, for every $\varepsilon\in(0,1)$ and every $r>0$, there exists a $\delta=\delta(\varepsilon,r)>0$ such that

[TABLE]

whenever $\left|x_{0}\right|<\delta$. Otherwise, the trivial solution is said to be stochastically unstable.

Definition 3** (Stochastic Asymptotic Stability (SAS))**

The trivial solution of (1) is said to be stochastically asymptotically stable if it is stochastically stable and, for every $\varepsilon\in(0,1)$, there exists a $\delta_{0}=\delta_{0}(\varepsilon)>0$ such that

[TABLE]

whenever $\left|x_{0}\right|<\delta_{0}$.

Definition 4** (SAS in the Large)**

The trivial solution of (1) is said to be stochastically asymptotically stable in the large if it is stochastically stable and, for all $x_{0}\in\mathbb{R}^{n}$,

[TABLE]

Definition 5** (Differential Operator $L$)**

Define the differential operator $L$ associated with equation (1) by

[TABLE]

If $L$ acts on a function $V\in C^{2}\left(\mathcal{D};\mathbb{R}\right)$, then

[TABLE]

where $V_{x}=\left(\frac{\partial V}{\partial x_{1}},\cdots,\frac{\partial V}{\partial x_{n}}\right),$ V_{xx}=\left(\frac{\partial^{2}V}{\partial x_{i}\partial x_{j}}\right)_{n\times n}=\left(\begin{array}[]{ccc}\frac{\partial^{2}V}{\partial x_{1}\partial x_{1}}&\cdots&\frac{\partial^{2}V}{\partial x_{1}\partial x_{n}}\\ \vdots&&\vdots\\ \frac{\partial^{2}V}{\partial x_{n}\partial x_{1}}&\cdots&\frac{\partial^{2}V}{\partial x_{n}\partial x_{n}}\end{array}\right), and $\operatorname{Tr}[\cdot]$ indicates the trace of a matrix.

By Itô’s formula, if $X_{t}\in\mathcal{D}$, then

[TABLE]

Theorem 6

[25*, Theorem 2.2]**

If there exists a positive-definite function $V\in$ $C^{2}\left(\mathcal{D};\mathbb{R}_{+}\right)$ such that $LV(x)\leq 0$ for all $x\in\mathcal{D}$, where $\mathcal{D}$ is an open set containing the origin, then the trivial solution of (1) is stochastically stable.*

Theorem 7

[25*, Theorem 2.3]**

If there exists a positive-definite function $V\in$ $C^{2}\left(\mathcal{D};\mathbb{R}_{+}\right)$ such that $LV(x)<0$ for all $x\in\mathcal{D}\setminus\{0\}$, where $\mathcal{D}$ is an open set containing the origin, then the trivial solution of (1) is stochastically asymptotically stable.*

II-B Zubov’s Equation and Region of Attraction

Definition 8** (Sublevel-$c$ Set)**

Given $V:\mathbb{R}^{n}\rightarrow\mathbb{R}$, for any $c\in\mathbb{R}$, a set of the form

[TABLE]

is called a sublevel-$c$ set of $V$.

Definition 9** ($p$-Region of Attraction)**

For any probability level $p\in[0,1]$, we define the $p$-region of attraction for the origin:

[TABLE]

which is the largest probabilistic region of attraction at level $p$, consisting of all $x_{0}$ for which the solution $X_{t}(x_{0})$ converges to the origin with probability at least $p$.

Definition 10

Denote by $\mathcal{C}$ the set of points that are attracted to the equilibrium with positive probability, i.e.,

[TABLE]

Similar to how Zubov’s equation characterizes the domain of attraction for deterministic systems, Zubov’s equation can be used to characterizes probabilistic domains of attraction for stochastic systems [26, 18, 16]. We present the most relevant part as a preliminary to our work below.

Introduce a value function $W:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined as

[TABLE]

where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is any Lipschitz continuous function that is positive definite with respect to the origin. Under suitable assumptions [17, Theorem 2.2], it can be shown that $W$ is the unique solution (in viscosity sense) to the PDE

[TABLE]

We refer to this PDE as the stochastic Zubov equation for (1). Furthermore, $\mathcal{C}$ is completely characterized by the strict sublevel-$1$ set of $W$, i.e.,

[TABLE]

Since $W$ satisfies $LW(x)=-g(x)(1-W(x))$ and $g$ is positive definite, $W(x)<1$ if and only if $LW(x)<0$ when $x\neq 0$. Furthermore, $W(0)=0$ and $LW(0)=0$. Therefore, the sets $\{x\in\mathbb{R}^{n}:W(x)<1\}$ and $\{x\in\mathbb{R}^{n}:LW(x)<0\}$ are equivalent, provided that $W$ is positive definite. Consequently, by Theorem 7, accurately solving the stochastic Zubov PDE (4) can lead to a valid Lyapunov function for (1) and an approximation of the probabilistic domain of attraction $\mathcal{C}$. We will also provide sufficient conditions for verifying $p$-regions of attraction in Section IV.

III Training Neural Lyapunov Functions via PINNs

Following [8], we approximate the unique viscosity solution $W$ to the stochastic Zubov equation (4) with a neural network $W_{\mathrm{NN}}(x;\theta)$ through a physics-informed neural network (PINN) framework.

Let $W_{\mathrm{NN}}(x;\theta)$ denote a neural network approximation for solving the stochastic Zubov PDE (4). The loss function consists of three terms:

[TABLE]

which are defined as follows:

Residual loss $L_{r}$: This term measures the error in satisfying the PDE over a set of collocation points $S=\{x_{i}\}_{i=1}^{N}\subseteq X$ (with $X$ being a compact training domain). For example, one may define

[TABLE] 2. 2.

Boundary loss $L_{b}$: This term enforces the boundary conditions. We require that $W(x)=0$ for $x=0$. The boundary loss is defined as

[TABLE] 3. 3.

Data loss $L_{d}$: By an optional set of data points $S_{d}=\left\{y_{i}\right\}_{i}^{N_{d}}$ on which approximate values $\hat{W}(y_{i})$ for $W(y_{i})$ can be obtained. According to the analysis in Section II-B, $W\left(y_{i}\right)$ can be approximated by simulating $X_{t}(y_{i})$ for a sufficiently long period of time and use (3) to compute $\hat{W}\left(y_{i}\right)$. Based on the values of $\left\{\hat{W}\left(y_{i}\right)\right\}_{i}^{N_{d}}$, the data loss $L_{d}$ can be defined as

[TABLE]

This term incorporates additional data from stochastic simulations to further improve accuracy.

IV Formal Verification of the Learned Neural Lyapunov Function

Due to potential approximation errors and the lack of formal guarantees in neural network approximations, as highlighted in [7, 8, 2], we employ SMT solvers for formal verification. In the following sections, we provide sufficient stochastic Lyapunov conditions for probabilistic attraction, which are the main contribution of this paper.

IV-A Verification of Local Stability

In this section, we assume that Assumption 1 holds. Rewrite equation (1) as

[TABLE]

where $A=Df(0)\in\mathbb{R}^{n\times n},$ and for $k=1,\ldots,m$, $S_{k}=D\sigma_{k}(0)\in\mathbb{R}^{n\times n},$ with $\sigma_{k}(x)$ denoting the $k$-th column of $\sigma(x)$. The matrix $C$ is given by $C=[S_{1},S_{2},\ldots,S_{m}].$ Here, the remaining terms are defined as $m(x)=f(x)-Ax\quad\text{and}\quad n(x)=\sigma(x)-Cx,$ which satisfy $\lim_{x\to 0}\frac{\|m(x)\|}{\|x\|}=0\quad\text{and}\quad\lim_{x\to 0}\frac{\|n(x)\|}{\|x\|}=0.$ Here, $\|\cdot\|$ denotes the matrix 2-norm. When it applies to vectors, $\|\cdot\|$ coincides with vector 2-norm $|\cdot|$ (Euclidean norm).

In [27, Theorem 3.4], local Lyapunov functions for stochastic systems were studied and a corresponding theorem was established. In this paper, we offer a simpler proof under suitable assumptions.

Proposition 11

Let $A$ be a Hurwitz matrix. For any symmetric positive definite matrix $Q\in\mathbb{R}^{n\times n}$, suppose there exists a symmetric positive definite matrix $P\in\mathbb{R}^{n\times n}$ satisfying the stochastic Lyapunov equation

[TABLE]

Define $V_{P}(x)=x^{T}Px.$ Denote $LV_{P}(x)$ by $h(x)$ and define

[TABLE]

with $D^{2}h(x)$ being the Hessian matrix of $h(x)$. For some sufficiently small $\varepsilon>0,$ set $r=\lambda_{\min}(Q)-\varepsilon>0.$ If there exists a constant $c>0$ such that

[TABLE]

then

[TABLE]

Proof:

We begin by computing the generator $LV_{P}(x)$ for the function $V_{P}(x)=x^{T}Px$. By definition, using Itô’s formula, we have

[TABLE]

Expanding the trace term111We have $C^{T}PC=[(S_{i}x)^{T}P(S_{j}x)]_{i,j=1}^{m},$ $\operatorname{Tr}[C^{T}PC]=\sum_{i=1}^{m}x^{T}S_{i}^{T}PS_{i}x.$ in (8) gives

[TABLE]

Thus, we can write

[TABLE]

By (6), it follows that

[TABLE]

Since $\sigma(0)=n(0)=0$, we have $LV_{P}(0)=0$. Moreover, since both $m(x)$ and $n(x)$ vanish faster than $x$ as $x$ approaches zero, the dominant term in the expansion of $LV_{P}(x)$ is quadratic with a coefficient of $-Q$. Define $h(x)=LV_{P}(x)$. Then $h(0)=0$, $\nabla h(0)=0$ and the Hessian at $x=0$ satisfies $D^{2}h(0)=-2Q.$

By Taylor’s theorem in several variables, we can express

[TABLE]

where the remainder $R(x)$ can be written as

[TABLE]

Define

[TABLE]

Then,

[TABLE]

Assume that, for all $x\in V_{P}^{c}$, (7) holds. Then, by estimating the remainder term, we have

[TABLE]

Thus, $LV_{P}(x)=h(x)<0$ for all $x\in V_{P}^{c}\setminus\{0\}$. ∎

Remark 12

For any matrix $A\in\mathbb{R}^{n\times n}$, we have

[TABLE]

which shows that the 2-norm and Frobenius norm are equivalent. In practice, we use the Frobenius norm due to its computational simplicity, particularly in the context of formal verification with SMT solvers.

Since $V_{P}^{c}$ is a convex set containing 0, the implication

[TABLE]

is equivalent to (7). For any given $r>0$, condition (7) can be verified by dReal [23]. Since $\|M(0)\|_{F}=0<2r$ and $P$ is positive definite, one can always choose $c>0$ sufficiently small so that (7) holds. In rare cases, if (7) holds for all $x\in\mathbb{R}^{n}$, then the solution of (1) is stochastically asymptotically stable in the large.

IV-B Verification of Probabilistic Regions of Attraction

In [28, Theorem 4.4], reach-avoid-stay and safe-stability for stochastic systems are established using a single Lyapunov function. Here, we reformulate the probabilistic safe stability by employing the concept of sublevel sets. This reformulation provides a precise definition of probabilistic attraction and yields a straightforward, verifiable condition for ensuring probabilistic safe stability.

Theorem 13

Let Assumption $1$ hold and $\mathcal{D}\subseteq\mathbb{R}^{n}$. Consider a function $V\in C^{2}\left(\mathbb{R}^{n};\mathbb{R}_{\geq 0}\right)$. Let $c_{2}>0$ be such that $V^{c_{2}}\subset\mathcal{D}$. The following holds:

1. i)

If $LV(x)\leq 0$ for all $x\in V^{c_{2}}$, then for any $x_{0}\in V^{c_{2}}$, we have

[TABLE] 2. ii)

If $LV(x)<0$ for all $x\in V^{c_{2}}\setminus\{0\}$ and $V$ is positive definite on $\mathcal{D}$, then, for any $x_{0}\in V^{c_{2}}$, we have

[TABLE] 3. iii)

If there exists $0<c_{1}<c_{2}$ and $\zeta>0$ such that $LV(x)\leq-\zeta$ for all $x\in V^{c_{2}}\setminus V^{c_{1}}$. Then, for any $x_{0}\in V^{c_{2}}\setminus V^{c_{1}}$,

[TABLE]

Proof:

i) For any initial condition $x_{0}\in V^{c_{2}}$, define $\tau$ be the first exit time of $X_{t}(x_{0})$ from $V^{c_{2}}$, that is

[TABLE]

If $X_{t}(x_{0})\in V^{c_{2}}$ for all $t\geq t_{0}$, then $\tau=\inf\emptyset=\infty$. In the following proof, we simplify the notation by writing $X_{s}(x_{0})$ as $X_{s}$. By Itô’s formula,

[TABLE]

For any fixed $t\geq t_{0}$, taking expectation on both sides and making use of $LV(X_{t}(x_{0}))\leq 0$ for all $X_{t}(x_{0})\in V^{c_{2}}\subseteq\mathcal{D}$, we can get

[TABLE]

Considering $\mathbb{E}V(X_{\tau\wedge t})=\mathbb{E}[(I_{\{\tau\leq t\}}+I_{\{\tau>t\}})V(X_{\tau\wedge t})]$ and $V(X_{t})\geq 0$ for all $X_{t}\in\mathcal{D}$, we can get

[TABLE]

By the definitions of $V^{c_{2}}$ and $\tau$ and continuity of $V(x)$ and $X_{t}$, if $\tau\leq t$, we know that $V(X_{\tau\wedge t}(x_{0}))=V(X_{\tau}(x_{0}))=c_{2}$. Then we can get

[TABLE]

By (9), (10) and (11), we can get

[TABLE]

Therefore, $\mathbb{P}\{\tau\leq t\}\leq\frac{V\left(x_{0}\right)}{c_{2}}.$ By letting $t\rightarrow\infty$, we obtain $\mathbb{P}\{\tau<\infty\}\leq\frac{V\left(x_{0}\right)}{c_{2}}$. It follows that, for all $x_{0}\in V^{c_{2}}$,

[TABLE]

Since $\tau=\infty$ is equivalent to $X_{t}(x_{0})\in V^{c_{2}}$ for all $t\geq t_{0}$. We can rewrite (12) as

[TABLE]

ii) Fix an arbitrary $r\geq 0$ such that $\mathcal{B}_{r}(0)\subseteq V^{c_{2}}$. Note that $V$ satisfies the conditions in Theorem 6. By the definition of stochastic stability, it follows that there exists $\eta\in\left(0,r\right)$ such that for any $\varepsilon\in(0,1)$, $x(\tau^{*},x_{0})\in\overline{\mathcal{B}}_{\eta}(0)$ implies $\mathbb{P}\{X_{t}(x_{0})\in\mathcal{B}_{r}(0),\,\forall t\geq\tau^{*}\}\geq 1-\varepsilon$. Next, we define $\tau^{*}$ to be the first entry time of $X_{t}(x_{0})$ to $\mathcal{B}_{\eta}(0)$, that is

[TABLE]

Let $\xi=-\max_{x\in V^{c_{2}}\setminus\mathcal{B}_{\eta}(0)}LV(x)$. By Itô’s formula and $LV<-\xi$ on $V^{c_{2}}\setminus\mathcal{B}_{\eta}(0)$, we have

[TABLE]

Since on $\left\{\tau^{*}\wedge\tau\geq t\right\}$ we have $\tau^{*}\wedge\tau\wedge t=t$,

[TABLE]

Combining this with (14), we have

[TABLE]

By letting $t\rightarrow\infty$, we get $\mathbb{P}\left[\tau^{*}\wedge\tau\geq\infty\right]=0$, which implies $\mathbb{P}\left[\tau^{*}\wedge\tau<\infty\right]=1$. On $\{\tau=\infty\}$, we have $\mathbb{P}\left[\tau^{*}<\infty\right]=1$ and

[TABLE]

Since $\varepsilon$ and $r$ are arbitrary, we can obtain

[TABLE]

Since $\tau=\infty$ is equivalent to $X_{t}(x_{0})\in V^{c_{2}}$ for all $t\geq t_{0}$. We can rewrite (14) as

[TABLE]

Furthermore, combining this with (13), we can obtain

[TABLE]

iii) The proof follows the same strategy as in ii). For any initial condition $x_{0}\in V^{c_{2}}\setminus V^{c_{1}}$, let $\tau$ be as defined in the proof of part i) and $\sigma$ be the first entry time of $X_{t}(x_{0})$ to $V^{c_{1}}$, given by

[TABLE]

For $t\leq\tau\wedge\sigma$, $X_{t}(x_{0})\in V^{c_{2}}\setminus V^{c_{1}}$, where $LV(X_{t}(x_{0}))\leq-\zeta$. We have

[TABLE]

Since on $\left\{\sigma\wedge\tau\geq t\right\}$ we have $\sigma\wedge\tau\wedge t=t$,

[TABLE]

Combining this with (16), we have

[TABLE]

By letting $t\rightarrow\infty$, we get $\mathbb{P}\left[\sigma\wedge\tau\geq\infty\right]=0$, which implies $\mathbb{P}\left[\sigma\wedge\tau<\infty\right]=1$. On $\{\tau=\infty\}$, we have $\mathbb{P}\left[\sigma<\infty\right]=1$, that is,

[TABLE]

Since $\tau=\infty$ is equivalent to $X_{t}(x_{0})\in V^{c_{2}}$ for all $t\geq t_{0}$, we can rewrite (17) as

[TABLE]

With a slightly modified argument as in the proof of part ii), we can show that

[TABLE]

Thus,

[TABLE]

∎

In our computation, both quadratic Lyapunov functions and neural Lyapunov functions are used to certify probabilistic attraction. For this purpose, we need to formulate the following result. Given one Lyapunov function, we seek to use another to enlarge the probabilistic region of attraction, as stated in the following theorem.

Theorem 14

Let Assumption $1$ hold and $\mathcal{D}\subseteq\mathbb{R}^{n}$. Let $V\in C^{2}\left(\mathcal{D};\mathbb{R}_{\geq 0}\right)$ satisfy the conditions in Theorem 13 iii) for some $c_{2}>0$. Let $c_{1}\in(0,c_{2})$. Assume that there exists a function $W\in C^{2}\left(\mathcal{D};\mathbb{R}_{\geq 0}\right)$ and some $\zeta>0$ satisfying

[TABLE]

where $\beta_{2}>\beta_{1}>0$ and $W^{\beta_{2}}\subseteq\mathcal{D}$. Furthermore, assume that $W^{\beta_{1}}\subseteq V^{c_{1}}\subseteq V^{c_{2}}\subseteq W^{\beta_{2}}\subseteq\mathcal{D}$. Then, for any $x_{0}\in W^{\beta_{2}}$, we have

[TABLE]

That is, for all $x_{0}\in W^{\beta_{2}}$, $X_{t}(x_{0})$ will remain in $W^{\beta_{2}}$ and converge to the equilibrium point with probability at least $(1-\dfrac{c_{1}}{c_{2}})(1-\dfrac{W(x_{0})}{\beta_{2}})$.

Proof:

For any initial condition $x_{0}\in W^{\beta_{2}}\setminus W^{\beta_{1}}$, denote by $E_{1}$ the following set:

[TABLE]

It follows from (18), Theorem 13 iii) (applied to $W$), and $W^{\beta_{1}}\subseteq V^{c_{1}}$ that, for any $x_{0}\in W^{\beta_{2}}\setminus W^{\beta_{1}}$,

[TABLE]

Define $\tau$ be the first time $X_{t}(x_{0})$ enters the set $V^{c_{1}}$. The state $X_{\tau}(x_{0})$ serves as the initial condition for the evolution of $X_{t}(X_{\tau}(x_{0}))$.

Denote by $E_{2}$ the following set:

[TABLE]

By the assumptions on $V$ and Theorem 13 ii), we obtain

[TABLE]

It follows that

[TABLE]

∎

Taking the computed quadratic Lyapunov function as $V$ and the neural Lyapunov function as $W$, by Theorem 14, we obtain verified probabilistic ROAs for the neural Lyapunov function. This approach effectively resolves the issue that a neural Lyapunov function may not be well-trained and fail to satisfy the Lyapunov condition around the origin. We derive a general attraction probability function using a piecewise formulation that characterizes a conservative lower bound on the attraction probability for an initial point in $\mathcal{C}$, as stated in the following proposition.

Proposition 15

Let $V$ and $W$ be the functions and $\beta_{1}$, $\beta_{2}$, $c_{1}$, and $c_{2}$ be the constants defined in Theorem 14. Then, the attraction probability for an initial state $x_{0}\in W^{\beta_{2}}$, $\mathbb{P}(x_{0})$, defined as

[TABLE]

is at least

[TABLE]

Furthermore, $p(x_{0})$ is continuous with respect to $x_{0}$ and serves as a conservative lower bound on the attraction probability.

Proof:

The fact that $\mathbb{P}(x_{0})\geq p(x_{0})$ follows immediately from Theorems 13 and 14. Indeed, for any $x_{0}$, we have $\mathbb{P}(x_{0})\geq 1-V(x_{0})/c_{2}$ by Theorem 13 iii) and $\mathbb{P}(x_{0})\geq(1-W(x_{0})/\beta_{2})(1-c_{1}/c_{2})$ by Theorem 14. Note that these bounds hold even if $V(x_{0})>c_{2}$ or $W(x_{0})>\beta_{2}$ (although then the bounds will be trivial).

We now show that $p(x_{0})$ is continuous in $x_{0}$. We only need to verify this at every $x_{0}$ such that $W(x_{0})=\beta_{1}$. Since $W^{\beta_{1}}\subset V^{c_{1}}$, we have $V(x_{0})\leq c_{1}$. Hence, $1-V(x_{0})/c_{2}\geq 1-c_{1}/c_{2}>(1-W(x_{0})/\beta_{2})(1-c_{1}/c_{2})$. By continuity of both $V(x)$ and $W(x)$ on $\mathcal{D}$, we have $1-V(x_{0})>(1-W(x_{0})/\beta_{2})(1-c_{1}/c_{2})$ holds in a small neighborhood of $x_{0}$, on which we have $p(x_{0})=1-V(x_{0})/c_{2}$. Hence, $p(x_{0})$ is continuous at $x_{0}$ because $V$ is. ∎

Remark 16

To apply Proposition 15, we need to determine the constants $\beta_{1}$, $\beta_{2}$, $c_{1}$, and $c_{2}$ as described in Theorem 14. We can choose them as follows to optimize the lower bound $p(x_{0})$ given by Proposition 15. Suppose we already have a Lyapunov function $V$ such that $LV(x)<0$ for all $x\in V^{c_{2}}\setminus\{0\}$, where $c_{2}$ is chosen as the largest such value. A premise for applying Theorem 14 to further enlarge the probabilistic region of attraction is the existence of another Lyapunov function $W$ such that $LW(x)<0$ holds on $W^{\beta_{2}}\setminus W^{\beta_{1}}$, where $W^{\beta_{1}}\subseteq V^{c_{1}}\subseteq V^{c_{2}}\subseteq W^{\beta_{2}}$. We need $V^{c_{2}}\subseteq W^{\beta_{2}}$ so that the region of attraction can be enlarged using the Lyapunov function $W$, and we need $W^{\beta_{1}}\subseteq V^{c_{1}}$ to draw conclusions about probabilistic attraction using the Lyapunov function $V$. Hence, $\beta_{1}$ is chosen as the smallest and $\beta_{2}$ as the largest such that $LW(x)<0$ holds on $W^{\beta_{2}}\setminus W^{\beta_{1}}$. The constant $c_{1}$ is the smallest value ensuring that $\{x:W(x)\leq\beta_{1}\}\subseteq\{x:V(x)\leq c_{1}\}$. Clearly, minimizing $\beta_{1}$ and $c_{1}$, and maximizing $\beta_{2}$, leads to less conservative results in terms of $p$-ROA estimates.

V Numerical Examples

In this section, we present numerical experiments on stochastic systems. To demonstrate that the proposed PINN approach reliably approximates $W(x)$, we extend the LyZNet tool [13] to learn a neural Lyapunov function and use dReal [23] to formally verify the stochastic Lyapunov condition. This enables the characterization of the largest probabilistic region of attraction $\mathcal{C}$. The results confirm that the learned functions capture the probabilistic attraction accurately and quantify the attraction probabilities via Proposition 15.

Example 17** (Reversed Van der Pol)**

Consider the reversed Van der Pol oscillator

[TABLE]

which has a stable equilibrium point at the origin. The linearization at the origin is given by (5) with A=\left[\begin{array}[]{ll}0&-1\\ 1&-1\end{array}\right], S_{1}=\left[\begin{array}[]{ll}\alpha&0\\ 0&0\end{array}\right] and S_{2}=\left[\begin{array}[]{ll}0&0\\ 0&\beta\end{array}\right]. Taking the diffusion parameters as $\alpha=\beta=0.5$ and solving the stochastic Lyapunov equation (6) with $Q=I$ gives P=\left[\begin{array}[]{cc}2.2439&-0.7805\\ -0.7805&1.4634\end{array}\right]. With $r=0.9999$, we can verify (7) with $c=0.3317$. Next, we preform an additional step to check that $LV_{p}\leq-\epsilon$ holds on the ellipsoid $\Omega_{c}=\left\{x\in\mathbb{R}^{n}:x^{T}Px\leq c\right\}$ with $c=2.3151$, for some small $\epsilon>0$. Hence, $\Omega_{c}$ provides a verified local probabilistic ROA for the origin.

In order to better capture the largest probabilistic ROA $\mathcal{C}$, we employ a neural network with three hidden layers (10 neurons per layer) trained using Zubov’s PDE (4). For data generation, we estimate the expectation in (3) by running 100 simulations for each initial point. In this work, we set $g(x)=0.1\|x\|^{2}$, same as the setting for the deterministic system in [8]. The learned neural Lyapunov functions are formally verified using dReal. Additionally, we compare the verified probabilistic ROA $\mathcal{C}$ from our neural Lyapunov function with the Lyapunov function obtained via the SOS method.

The results are shown in Figure 1. In the right panel of Figure 1, the green curve represents the largest verified probabilistic ROA in terms of a sublevel set of the neural Lyapunov function on which the stochastic Lyapunov condition is satisfied. The pink curve shows the verified sublevel set computed via an SOS Lyapunov function, and the red curve corresponds to that obtained from a quadratic Lyapunov function. Note that the neural Lyapunov function yields a larger verified $\mathcal{C}$ than the other two approaches. Additionally, the black curve marks the lower bound of the region where the neural Lyapunov function satisfies the stochastic Lyapunov condition, while the blue curve indicates a lower bound for the quadratic Lyapunov function that encloses the region where the neural Lyapunov function fails to verify the condition.

Based on the verified level sets shown in Figure 1 and by applying Proposition 15, we generate a heatmap in the left image of Figure 2 to visualize the attraction probability for the initial point in $\mathcal{C}$. As depicted in the right panel of Figure 2, the noise term in the stochastic differential equation can sometimes enlarge the region of attraction. In our experiments, we sampled points from within the probabilistic ROA while excluding those on the deterministic limit cycle (indicated by the blue curve, which marks the boundary of the deterministic ROA). We then used this sampling point as the initial value and simulate 10,000 trajectories from it to estimate the attraction probability from this point. The results show that over $10\%$ of the trajectories converged. This indicates that random disturbances can, under certain conditions, have a stabilizing effect on the stochastic system.

VI Conclusions and Future Works

We employ a data-driven PINN to approximate the solution of stochastic Zubov PDE for stochastic system and use dreal to formally verify the stochastic Lyapunov condition for the learned neural Lyapunov function both locally and over the largest probabilistic region of attraction $\mathcal{C}$. We also derive sufficient stochastic Lyapunov conditions for probabilistic attraction and propose an attraction probability function based on both the local and neural Lyapunov functions. In our experiments, we generate a heatmap to visualize the attraction probability. Our results demonstrate that the neural approach yields a larger verified probabilistic ROA than the Lyapunov function computed by the SOS method.

Future research will focus on improving scalability for high-dimensional systems. A significant challenge remains in addressing high-dimensional dynamic systems, which are inherently difficult to solve even with neural networks. This difficulty is particularly pronounced during training data generation and verification, as increased dimensionality demands substantially more computational resources and time. Another promising direction is to try other verification tools, e.g., [29], which would help enlarge the probabilistic region of attraction $\mathcal{C}$. Additionally, Proposition 15 only provides a conservative estimate of the $p$-ROA. It would be of interest to investigate how to optimize the attraction probabilities.

VII ACKNOWLEDGMENTS

The authors are grateful to Dr. Yiming Meng for helpful discussions.