# A neural network-based construction of Lyapunov functions based on homogeneous approximations

Artemiy Zenkin<fnref n="1" />, Rosane Ushirobira², Denis Efimov², Alexey Bobtsov<fnref n="1" />

*Abstract*— This paper applies an artificial neural network approach to design a local Lyapunov function for a general class of nonlinear systems. To avoid the singularity at the origin in computing the Lyapunov function, it is complemented by using a homogeneous approximation at the origin.

## I. INTRODUCTION

Control theory and engineering are fundamentally concerned with the stability of nonlinear dynamic systems [12]. A significant body of literature exists that explores various techniques to address this issue, with a primary focus on Lyapunov theory. The existence of a Lyapunov function is recognized as a necessary and sufficient condition for a broad class of dynamical systems to achieve asymptotic stability. However, the design of such functions is not straightforward, as no universal technique is available. The methods for identifying Lyapunov function candidates can be categorized into two main approaches: a) *analytical methods* are typically limited to specific canonical forms of the system's model; b) *numeric approximation techniques* are effective for systems with a small number of state variables, they often incur high computational costs when applied to large-scale systems.

Artificial neural networks (ANNs) have emerged as a popular tool for addressing the challenges associated with Lyapunov functions. The universal approximation theorem plays a crucial role in this context, stating that a neural network with at least one hidden layer can approximate any continuous function arbitrarily well [3], [11]. This theorem implies that the functions represented by neural networks are dense within the set of continuous functions. Many notable contributions have been made to the construction of Lyapunov functions using ANNs. For example, [14] proposed a straightforward method for constructing a Lyapunov function modeled by a neural network, utilizing weights determined by a genetic algorithm. The paper [8] introduced a deep neural network architecture designed for training procedures that approximate Lyapunov functions for nonlinear ordinary differential equations. It was demonstrated that if a system allows for a compositional Lyapunov function, the number of neurons required to estimate this function with fixed accuracy increases only polynomially with the state dimension. A control Lyapunov function was constructed in [13] using a neural network, which ensured the stability of the closed-loop system post-training. The authors in [16] developed a technique, which enhances the training of ordinary differential equations (ODEs) by improving prediction performance, accelerating inference dynamics convergence, and increasing adversarial robustness. In [15], a deep neural network architecture was utilized to derive a numerical control Lyapunov function for nonlinear dynamics, providing an estimate of the region of attraction and ensuring stability.

Nowadays, the exploration of Lyapunov functions within the framework of nonlinear systems continues to evolve, particularly with the integration of ANNs. These advancements not only enhance the understanding of stability, but also improve the practical applications of control theory in complex systems. As research progresses, further innovations are anticipated in the methodologies for constructing and applying Lyapunov functions, paving the way for more robust and efficient control strategies.

In particular, the study of numerical methods for solving partial differential equations (PDEs) has gained significant attention in recent years, particularly with the advent of neural network-based approaches. A numerical method was first introduced for solving PDEs using neural network-based functions in [5] (this method employs an unconstrained minimization problem, which allows for the effective approximation of solutions to complex equations), then the Physics Informed Neural Network (PINN) approach was recently described by Raissi et al. (2019). This innovative framework utilizes a fully connected ANN that incorporates a loss function encoding the physical laws represented by nonlinear PDEs. The integration of physics into the training process enhances the model's ability to generalize and accurately predict solutions. The paper [2] proposed enhancements to the PINN framework by combining it with a coupled-automatic-numerical differentiation approach. This combination facilitates a robust and efficient training of neural networks, thereby improving the overall performance of the PINN method. In fact, [8] utilized PINNs to approximate a Lyapunov function in the context of small-gain theory for large-scale systems. However, it was noted that this approach is generally applicable only for local stability analysis, with no guarantees in the vicinity of the origin. Moreover, the questions remained for conditions when such an ANN-based representation of a Lyapunov function exists. While challenges remain, particularly regarding global stability and approximation, ongoing research continues to refine these techniques.

According to [9], under specific nonlinear changes of variables, the global asymptotic stability of the origin for a family of ordinary differential equations is equivalent to global exponential stability, when a quadratic Lyapunov function is employed. This finding highlights the importance of variable transformations in the approximation of Lyapunov functions. In a recent paper [4], following [8], [9], an ANN-based approximation of a Lyapunov function was proposed in a quadratic canonical form, showing the natural obstacles appearing close to the origin. Moreover, a PINN design that efficiently approximates a global Lyapunov function for homogeneous systems was also presented.

Continuing in this vein, the present note considers the class of locally homogeneous systems close to the origin, and an ANN-based design of Lyapunov functions is proposed uniting two PINNs candidates, one for the homogeneous approximation in a vicinity of the origin and another for the original system on the complement set. Such a combination allows us to avoid the problems with the

<footnote n="1">¹ITMO University, 49 av. Kronverkskiy, 197101 Saint Petersburg, Russia
² Inria, Univ. Lille, CNRS, UMR 9189, F-59000 Lille, France</footnote>
non-definiteness of ANN-approximated Lyapunov function close to zero. Auxiliary results on estimation of accuracy of the approximation of Lyapunov function from the result on the training dataset are presented.

### Notation

* Denote by $\mathbb{R}$ and $\mathbb{R}_+$ the sets of real and nonnegative real numbers, respectively.

* The Euclidean norm of a vector in $\mathbb{R}^n$ is denoted by $\| \cdot \|$.

* For $F \in \mathbb{R}$ define $[F]_+ := \max\{0, F\}$ and $[F]_- := \min\{0, F\}$. If $F$ takes matrix or vector values, then these operators are understood elementwise.

* A continuous function $\alpha : \mathbb{R}_+ \to \mathbb{R}_+$ belongs to the class $\mathcal{K}$ if $\alpha(0) = 0$ and it is strictly increasing; $\alpha$ belongs to the class $\mathcal{K}_\infty$ if $\alpha \in \mathcal{K}$ and it is increasing to infinity. A continuous function $\beta : \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}_+$ belongs to the class $\mathcal{KL}$ if $\beta(\cdot, t) \in \mathcal{K}$ for each fixed $t \in \mathbb{R}_+$ and $\beta(s, \cdot)$ is decreasing to zero for each fixed $s > 0$.

* The notation $DV(x)f(x)$ stands for the directional derivative of a continuously differentiable function $V : \mathbb{R}^n \to \mathbb{R}$ with respect to a vector field $f : \mathbb{R}^n \to \mathbb{R}^n$ evaluated at the point $x$, then $DV(x)$ is the gradient; and this notation represents the generalized directional Dini derivative of $V$ at $x$ in the direction $f(x)$ if $V$ is just locally Lipschitz continuous.

* For a set $A \subset \mathbb{R}^n$, denote its boundary and interior by $\partial A$ and $\text{int}(A)$, respectively.

* The finite set of the first $n$ positive integers is denoted by $\overline{1, n} := \{1, 2, ..., n\}$.

## II. PRELIMINARIES

Consider a continuous-time system:

$$ \dot{x}(t) = f(x(t)), \forall t \in \mathbb{R}_+, \quad (1) $$

where $x(t) \in \mathbb{R}^n$ represents the state function<fnref n="1" /> and $f : \mathbb{R}^n \to \mathbb{R}^n$ is a continuous vector field, locally Lipschitz continuous on $\mathbb{R}^n \setminus \{0\}$, with $f(0) = 0$. We suppose that the system (1) with initial conditions $x_0 \in \mathbb{R}^n$ has a solution $x(t, x_0)$ defined for all $t \in \mathbb{R}_+$ (the system is forward complete, see [12] for basic definitions of properties of dynamical systems).

### A. Stability notions

We say that $x = 0$ is a globally asymptotically stable (GAS) equilibrium of (1) if there exists function $\beta \in \mathcal{KL}$ such that

$$ \|x(t, x_0)\| \le \beta(\|x_0\|, t) $$

for all $t \in \mathbb{R}_+$ and all $x_0 \in \mathbb{R}^n$. Also, $x = 0$ is globally exponentially stable (GES), if there exist $c \ge 1$ and $\alpha > 0$ such that

$$ \|x(t, x_0)\| \le c e^{-\alpha t} \|x_0\|, \quad (2) $$

for all $t \in \mathbb{R}_+$ and all $x_0 \in \mathbb{R}^n$. If $\mathcal{O} \subset \mathbb{R}^n$ is an open set in the domain of attraction of the origin for (1), then the same properties can be defined locally for $x_0 \in \mathcal{O}$.

**Definition 1.** A Lyapunov function on $\mathcal{O}$ is a continuously differentiable ($C^1$) function $V : \mathcal{O} \to \mathbb{R}_+$ such that $V(0) = 0$; $V(x) > 0$ for all $x \in \mathcal{O} \setminus \{0\}$, and the derivative of $V$ in the direction of

the vector field $f$ satisfies $DV(x)f(x) = -g(x) \le 0$, for some function $g : \mathbb{R}^n \to \mathbb{R}_+$, with $g(x) > 0$ for all $x \in \mathcal{O} \setminus \{0\}$.

In the case $\mathcal{O} = \mathbb{R}^n$ and $V(x) \to +\infty$ when $x \to +\infty$, we call $V$ a global Lyapunov function.

The existence of a (global) Lyapunov function for (1) is a necessary and sufficient condition for $0 \in \mathbb{R}^n$ to be (globally) asymptotically stable equilibrium, see for instance [12].

A significant result in [9], Theorem 2, claims that the system (1), supposed to be GAS, can be transformed into a GES system:

**Theorem 1** ([9]). Let $n \ne 4, 5$, and consider any system (1) on $\mathbb{R}^n$, which is GAS. Then (1) can be transformed into a system that is GES. In particular, in (2), $c = 1$ and $\alpha = 1$ can be chosen.

The proposed transformation is realized through what the authors in [9] call a change of variables, that is, a homeomorphism $T : \mathbb{R}^n \to \mathbb{R}^n$ satisfying $T(0) = 0$, $T$ is $C^1$ on $\mathbb{R}^n$, and a diffeomorphism on $\mathbb{R}^n \setminus \{0\}$. These transformations are infinitely differentiable everywhere, except at the origin, where they are just continuously differentiable. Their respective inverses are continuous globally and infinitely differentiable away from the origin. Then, for $y = T(x)$, in the new coordinates the dynamics of (1) takes the form:

$$ \dot{y} = DT(x)f(x) = DT(T^{-1}(y)) f(T^{-1}(y)) =: \tilde{f}(y). $$

Moreover, this change of variables can be used for definition of a Lyapunov function, as suggested by the following proposition.

**Proposition 1** ([9]). For $n \ne 4, 5$, let $V : \mathbb{R}^n \to \mathbb{R}$ be a proper<fnref n="2" />, positive definite, $C^1$ global Lyapunov function of (1), which is smooth on $\mathbb{R}^n \setminus \{0\}$. Then there exists a change of variables $T : \mathbb{R}^n \to \mathbb{R}^n$ with $DT(0) = 0$ such that $\tilde{V}(y) = V(T^{-1}(y)) = \|y\|^2$ is the transformed Lyapunov function with $D\tilde{V}(y)\tilde{f}(y) = -\|y\|^2$ for all $y \in \mathbb{R}^n$.

Therefore, in the original coordinates, we have the expressions:

$$ V(x) = \|T(x)\|^2 \quad \text{and} \quad DV(x)f(x) = -\|T(x)\|^2 \quad (3) $$

for all $x \in \mathbb{R}^n$.

### B. Homogeneity

1) **Definitions of homogeneity:** For any $r_i > 0$, $i \in \overline{1, n}$, and $\lambda > 0$, define the vector of weights $\mathbf{r} = (r_1, \dots, r_n)$ and the dilation matrix $\Lambda_r(\lambda) = \text{diag}\{\lambda^{r_i}\}_{i=1}^n$. Also denote $r_{\min} = \min_{i \in \overline{1, n}} \{r_i\}$ and $r_{\max} = \max_{i \in \overline{1, n}} \{r_i\}$.

**Definition 2** ([6], [18]). 1) A function $h : \mathbb{R}^n \to \mathbb{R}$ is called r-homogeneous, if for some $\nu \in \mathbb{R}$ the relation $h(\Lambda_r(\lambda)x) = \lambda^\nu h(x)$ holds for all $x \in \mathbb{R}^n$ and any $\lambda > 0$.

2) A vector field $f : \mathbb{R}^n \to \mathbb{R}^n$ is called r-homogeneous, if for some $\nu \ge -r_{\min}$ the relation $f(\Lambda_r(\lambda)x) = \lambda^\nu \Lambda_r(\lambda)f(x)$ holds for all $x \in \mathbb{R}^n$ and any $\lambda > 0$.

In both cases, the constant $\nu$ is called the degree of homogeneity.

**Definition 3.** A dynamical system as in (1) is called r-homogeneous of degree $\nu$ if this property is satisfied for the vector field $f$ in the sense of Definition 2.

**Definition 4.** For $x \in \mathbb{R}^n$ and $\varpi \ge r_{\max}$, the r-homogeneous norm is defined as follows $\|x\|_r = (\sum_{i=1}^n |x_i|^{\varpi/r_i})^{1/\varpi}$.

<footnote n="1">¹In this note, we will sometimes omit the time dependency of the variables if it does not lead to a confusion.</footnote>

<footnote n="2">²That is, for each $\beta \ge 0$, the set $\{x \in \mathbb{R}^n | V(x) \le \beta\}$ is bounded.</footnote>
Notice that, for all $x \in \mathbb{R}^n$, its Euclidean norm $\|x\|$ is related with the homogeneous one, $\underline{\sigma}_r(\|x\|_{\mathbf{r}}) \le \|x\| \le \bar{\sigma}_r(\|x\|_{\mathbf{r}})$ for some $\underline{\sigma}_r, \bar{\sigma}_r \in \mathcal{H}_{\infty}$ [6].

In the following, due to this norm equivalence, the stability analysis with respect to the norm $\| \cdot \|$ can be replaced with the analysis for the norm $\| \cdot \|_{\mathbf{r}}$. The homogeneous norm has an important property: it is $\mathbf{r}$-homogeneous of degree 1, that is $\|\Lambda_{\mathbf{r}}(\lambda)x\|_{\mathbf{r}} = \lambda\|x\|_{\mathbf{r}}$ for all $x \in \mathbb{R}^n$ and $\lambda > 0$.

2) **Lyapunov functions:** One of the main features of homogeneous systems related to the Lyapunov function method is described in the following theorem:

**Theorem 2** ([17], [19]). For an $\mathbf{r}$-homogeneous system (1) of degree $\nu$, if the origin is GAS, then there exists an $\mathbf{r}$-homogeneous function $V: \mathbb{R}^n \to \mathbb{R}_+$ of degree $\mu > -\nu$, $C^1$ on $\mathbb{R}^n$, smooth on $\mathbb{R}^n \setminus \{0\}$, and such that for all $x \in \mathbb{R}^n$:

$$a \|x\|_{\mathbf{r}}^{\mu} \leq V(x) \leq b \|x\|_{\mathbf{r}}^{\mu}, \quad (4)$$

$$DV(x)f(x) \leq -cV^{1+\frac{\nu}{\mu}}(x) \leq -d\|x\|_{\mathbf{r}}^{\nu+\mu}, \quad (5)$$

for some $0 < a \le b$, $c > 0$ and $d > 0$.

We describe next an interesting robustness result for homogeneous Lyapunov functions established in [6]:

**Proposition 2.** For a locally Lipschitz continuous and $\mathbf{r}$-homogeneous function $V: \mathbb{R}^n \to \mathbb{R}_+$ of degree $\mu > -\nu$, assume that the estimates (4), (5) are satisfied for constants $0 < a \le b$ and $d > 0$. Let $\varepsilon: \mathbb{R}^n \to \mathbb{R}$ be a locally Lipschitz continuous and $\mathbf{r}$-homogeneous function of degree $\mu$ such that $-a < \varepsilon, d > \bar{\varepsilon}$, where $^3\varepsilon = \inf_{y \in S_r(1)} \varepsilon(y)$ and $\bar{\varepsilon} = \sup_{y \in S_r(1)} D\varepsilon(y)f(y)$. Then $V'(x) = V(x) + \varepsilon(x)$ is a locally Lipschitz continuous and $\mathbf{r}$-homogeneous Lyapunov function for (1).

3) **Homogeneous approximations:** The $\mathbf{r}$-homogeneity property used in Definition 2 and Theorem 2 is introduced for some $\mathbf{r} > 0$ and all $\lambda > 0$. Restricting the set of admissible values for $\lambda$ we can introduce local homogeneity (in the 0-limit or in the $\infty$-limit):

**Definition 5.** [1] The function $g: \mathbb{R}^n \to \mathbb{R}$ is called $(\mathbf{r}, \lambda_0, g_0)$-homogeneous ($r_i > 0$, $i = \frac{1}{n}$; $g_0: \mathbb{R}^n \to \mathbb{R}$; $\lambda_0 \in \{0, +\infty\}$) if uniformly in $x \in S_r(1)$, $\lim_{\lambda \to \lambda_0} \lambda^{-d_0} g(\Lambda_{\mathbf{r}}(\lambda)x) - g_0(x) = 0$, for some $d_0 \ge 0$.

The vector function $f: \mathbb{R}^n \to \mathbb{R}^n$ is called $(\mathbf{r}, \lambda_0, f_0)$-homogeneous ($r_i > 0$, $i = \frac{1}{n}$; $f_0: \mathbb{R}^n \to \mathbb{R}^n$; $\lambda_0 \in \{0, +\infty\}$) if uniformly in $x \in S_r(1)$, $\lim_{\lambda \to \lambda_0} \lambda^{-d_0} \Lambda_{\mathbf{r}}^{-1}(\lambda)f(\Lambda_{\mathbf{r}}(\lambda)x) - f_0(x) = 0$, for some $d_0 \ge -r_{\min}$. The system (1) is called $(\mathbf{r}, \lambda_0, f_0)$-homogeneous if $f$ has this property.

The simplest example of local approximation at zero is the linearization.

**Theorem 3.** [10], [17], [18] Let the system (1) be $(\mathbf{r}, 0, f_0)$-homogeneous with the continuous functions $f: \mathbb{R}^n \to \mathbb{R}^n$ and $f_0: \mathbb{R}^n \to \mathbb{R}^n$. If the system $\dot{x} = f_0(x)$ is (locally) asymptotically stable, then the system (1) is also locally asymptotically stable.

**Theorem 4.** [1], [18] Let the system (1) be $(\mathbf{r}, +\infty, f_\infty)$-homogeneous with the continuous functions $f: \mathbb{R}^n \to \mathbb{R}^n$ and $f_\infty: \mathbb{R}^n \to \mathbb{R}^n$. If the system $\dot{x} = f_\infty(x)$ is globally asymptotically stable, then there exists a compact

invariant set $X_\infty \subset \mathbb{R}^n$ containing the origin such that the system (1) is globally asymptotically stable with respect to the set $X_\infty$.

The theorems 3 and 4 present results on stability of the system (1) derived from the corresponding properties of the approximating systems for $\lambda_0 = 0$ or $\lambda_0 = +\infty$. The converse Lyapunov theorem similar to Theorem 2 for the homogeneous in the bi-limit systems can also be found in [1].

### C. ANNs & PINNs

1) **Artificial neural networks:** An ANN is a graph-like machine learning model. Each node is an artificial neuron, and each edge connecting nodes is characterized by a weight. The number of incoming edges defines the degree of connection between neurons. Each neuron may be considered as a function that takes the output of nearby neurons as its input. Neurons are aggregated into L hidden layers, compressed between the input and output layers.

An ANN with more than one hidden layer ($L > 2$) between the input and output is called a deep neural network. For instance, in a feedforward neural network, the input is processed successively through the layers 1, 2, ..., $L$. We will use networks with the input vector $x = (x_1 \dots x_n)^\top \in \mathbb{R}^n$ and output $W(x; \theta) \in \mathbb{R}^m$, where $\theta \in \mathbb{R}^p$ is the vectors of parameters to be learned. There are $N_\ell$ neurons in the layer numbered $\ell$ and $y_k^\ell$ is the output of the neuron at position $k$ in the layer $\ell$, giving a vector $y^\ell = (y_1^\ell \dots y_{N_\ell}^\ell)^\top \in \mathbb{R}^{N_\ell}$. At the lowest level, $y^0 = x$ is the input, and the outputs of the neurons in the highest layer $y^L \in \mathbb{R}^m$ are used to compute the network's output $W(x; \theta)$:

$$y_k^\ell = \sigma^\ell \left( w_k^{\ell \top} y^{\ell-1} + b_k^\ell \right),$$

for $k \in \overline{1, N_\ell}$ and $\ell \in \overline{1, L}$, where $\sigma^\ell: \mathbb{R} \to \mathbb{R}$ is an activation function, $w_k^\ell = (w_{k,1}^\ell \dots w_{k,N_{\ell-1}}^\ell)^\top \in \mathbb{R}^{N_{\ell-1}}$ are the weights and $b_k^\ell \in \mathbb{R}$ are the biases (their elements are contained in $\theta$). Usually it is assumed that in the output layer ($\ell = L$) the activation function is a linear map, while other activation functions are sigmoidal<fnref n="4" />.

2) **Approximation of a continuous function by a neural network:** In this subsection, we recall a popular method of approximating a continuous function.

The works of [11] and [3] have established that, given enough hidden units, standard multilayer feedforward networks with any number of hidden layers and arbitrary squashing function are capable of approximating, with any desired degree of accuracy, any continuous function of $n$ variables with support in a compact set:

**Theorem 5** ([3]). Let $\phi: \mathbb{I}_n \to \mathbb{R}$ be a continuous function defined on the $n$-dimensional unit cube $\mathbb{I}_n$ in $\mathbb{R}^n$, and $\varepsilon > 0$. Then, there is $N_1 > 0$ such that

$$W(x; \theta) = \sum_{k=1}^{N_1} w_k^2 \sigma^1 \left( w_k^{1\top} x + b_k^1 \right) + b^2 \quad (6)$$

with $\sigma^1$ a sigmoidal function, $w_k^2, b_k^1, b^2 \in \mathbb{R}$ and $w_k^1 \in \mathbb{R}^n$ for all $k \in \overline{1, N_1}$ (components of $\theta$), satisfies

$$|W(x; \theta) - \phi(x)| < \varepsilon, \quad \forall x \in \mathbb{I}_n.$$

This result indicates that a single hidden layer is sufficient, or $L = 2$ is a common choice. An approximation of a function that incorporates information about its gradient is discussed in [7].

<footnote n="3">³For $a > 0$, we use the notation $S_r(a) := \{x \in \mathbb{R}^n \mid \|x\|_{\mathbf{r}} = a\}$.</footnote>

<footnote n="4">⁴A non-decreasing function $\sigma: \mathbb{R} \to \mathbb{R}$ is a sigmoidal function if $\lim_{t \to +\infty} \sigma(t) = 1$ and $\lim_{t \to -\infty} \sigma(t) = 0$.</footnote>
### 3) Application of ANNs for approximating Lyapunov functions:
We wish to select $\theta$ such that $W(x; \theta)$ approximates a Lyapunov function $V(x)$, for all $x \in K$, with $K$ being a fixed compact set including the origin for a GAS (at zero) system (1). To this end, according to Theorem 5, we can use one single hidden layer, so $L = 2$. The formula of $W(x; \theta)$ is given by (6).

Following [8], to train the network and approach a Lyapunov function by $W(x; \theta)$, a loss function $\mathcal{L}: \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$ must be defined. Then the training phase consists of finding parameters $\theta$ that search for the minimum of

$$J(\theta) = \frac{1}{N} \sum_{i=1}^{N} \mathcal{L} \left( W(x^{(i)}; \theta), DW(x^{(i)}; \theta), x^{(i)} \right)$$

where $x^{(i)} \in K$ are test points (elements of a learning/test dataset). Unlike conventional ANN learning tasks, the loss function $\mathcal{L}$ additionally depends on the test points' values of the derivative $DW$ of $W$ with respect to $x$. Such a constraint to establish if $W$ is a Lyapunov function [8]:

$$DW(x; \theta)f(x) \leq -\alpha(\|x\|), \forall x \in K,$$

for a suitably chosen $\alpha \in \mathcal{K}$. The appearance of $DW$ naturally opens a way of interpreting this training problem within a partial differential equation (PDE) framework. In such a case, the positivity of the Lyapunov function can be introduced through the boundary conditions, $W(0, \theta) = 0$, and $W(x; \theta) > 0, \forall x \in K \setminus \{0\}$. However, it can be challenging to implement the boundary conditions in the above form numerically: either because the first expression involves a single point, or because of the strict positive condition in the second one. We may then replace it with the following to palliate this:

$$\alpha_1 (\|x\|) \leq W(x; \theta) \leq \alpha_2 (\|x\|), \forall x \in K,$$

for some functions $\alpha_1$ and $\alpha_2 \in \mathcal{K}$.

Based on the above observations, in [8] the following loss function is considered:

$$\mathcal{L}(w, p, x) = [pf(x) + \alpha(\|x\|)]^2 + [w - \alpha_1 (\|x\|)]^2 + [w - \alpha_2 (\|x\|)]^2.$$

The drawback of this functional $\mathcal{L}$ is that its minimization and the obtained solution $W$ clearly depend on the choice of the tuning functions $\alpha, \alpha_1$, and $\alpha_2$, and it could be the case that, for a chosen set of these functions, there is no Lyapunov function for (1).

### 4) PINN: Physics-Informed Neural Networks (PINNs) are a kind of ANN models that use a priori information from an ODE or PDE system, which governs the dynamics under study and approximates a solution of this ODE/PDE for given initial conditions and domain of the time and/or space. This requires introduction of terms into the loss function of the ANN that are somehow proportional to the considered physics/dynamics information. A typical PINN uses a fully connected ANN architecture.

For some domain $\Omega \subset \mathbb{R}^n$ and $T > 0$, we consider solutions $\mathcal{V}(x, t)$ of a dynamical system described by differential equations of a general form:

$$\begin{aligned} \mathcal{D}_x(\mathcal{V}(x, t)) + \mathcal{D}_t(\mathcal{V}(x, t)) &= 0, \forall x \in \Omega, t \in (0, T]; \\ \mathcal{V}(x, 0) &= \mathcal{V}_0(x), \forall x \in \Omega; \\ \mathcal{B}(\mathcal{V}(x, t)) &= g(x, t), \forall x \in \partial\Omega, t \in (0, T], \end{aligned}$$

where $\mathcal{D}_x$ and $\mathcal{D}_t$ are differential operators in $x$ and $t$, $\mathcal{B}$ is a suitable boundary functional, and the functions $\mathcal{V}_0$ and $g$ determine the

initial and the boundary conditions. The PINN generates $W(x, t; \theta)$ predicting $\mathcal{V}(x, t)$ for $x \in \Omega$ and $t \in [0, T]$, where $\theta$ are properly learned parameters of ANN as before. The accuracy of the PINN output is determined by $\theta$, which is optimized with respect to the PINN loss function during the training:

$$\begin{aligned} \mathcal{L}_1(W, x, t) = & \nu_1 \|\mathcal{D}_t(W(x, t; \theta)) + \mathcal{D}_x(W(x, t; \theta))\|^2 \\ & + \nu_2 \|W(x, 0; \theta) - \mathcal{V}_0(x)\|^2 \\ & + \nu_3 \|\mathcal{B}(W(x, t; \theta)) - g(x, t)\|^2, \end{aligned}$$

where $\nu_i > 0, i \in \overline{1, 3}$ are tuning parameters. As we can conclude, the functional $\mathcal{L}_1$ depends on the derivative of $W$ with respect to $x$ and $t$, which are the inputs of ANN producing $W(x, t; \theta)$, then obviously the problem of approximation of a Lyapunov function for the system (1) by an ANN solved in [8] is an example of exploitation of PINN. The functional $\mathcal{L}$ from the previous subsection contains the structure of PDE describing the properties of a Lyapunov function.

### D. Canonical form approximation for Lyapunov functions [4]

Consider the vector field in (1). We wish to employ an ANN or PINN-based structure with an input vector $x = (x_1 \dots x_n)^\top \in \mathbb{R}^n$ and a scalar output $W(x; \theta) \in \mathbb{R}$ to approximate an unknown Lyapunov function $V(x)$ of (1) at the origin for $x \in K$, with $K$ being a compact set in $\mathbb{R}^n$ on which $V$ must be computed, $0 \in \text{int}(K)$. The vector $\theta \in \mathbb{R}^p$ represents the network's free parameters that must be learned to get the desired result. To this end, we will refine the results of [8] by utilizing the outcomes of Proposition 1 and the corresponding PDE (3). Therefore, in our case, the ANN has to approximate the change of variables $T$ defining $\mathcal{W}$ and $DW$ due to (3). For brevity, let us use one single hidden layer ANN consisting of $N_1$ neurons:

$$T_{\theta}(x) = w^2 \bar{\sigma}_1 (w^1 x + b^1) + b^2, \forall x \in K, \quad (7)$$

where the vector $\theta$ is composed by the elements of $w^1 \in \mathbb{R}^{N_1 \times n}$, $b^1 \in \mathbb{R}^{N_1}$, $w^2 \in \mathbb{R}^{n \times N_1}$ and $b^2 \in \mathbb{R}^n$, the function $\bar{\sigma}_1: \mathbb{R} \to \mathbb{R}^{N_1}$ contains copies of a sigmoidal function $\sigma_1: \mathbb{R} \to \mathbb{R}$, i.e., $\bar{\sigma}_1(s) = [\sigma_1(s_1) \dots \sigma_{N_1}(s_{N_1})]^\top$ for $s \in \mathbb{R}^{N_1}$. The ANN $T_{\theta}(x)$ should approximate the transformation map $T$. Since $T(0) = 0$, the same requirement can be imposed for $T_{\theta}$ by assigning

$$b^2 = -w^2 \bar{\sigma}_1 (b^1),$$

which also reduces the size of the vector $\theta$ for tuning.

The Lyapunov function candidate can be represented by

$$W(x; \theta) = T_{\theta}(x)^{\top} T_{\theta}(x), \forall x \in K \quad (8)$$

provided that the conditions for positive definiteness are verified, which in such a case follow by

$$(w^1)^{\top} w^1 \geq I_n, \quad (w^2)^{\top} w^2 \geq I_{N_1},$$

where $I_n$ denotes the identity matrix in $\mathbb{R}^{n \times n}$. The remaining condition of (3) on the decay of $W(x; \theta)$ along the trajectories of the system (1) can be presented as follows:

$$\begin{aligned} DW(x; \theta) f(x) &= 2T_{\theta}^{\top}(x) DT_{\theta}(x) f(x) \\ &\leq -W(x; \theta) = -T_{\theta}(x)^{\top} T_{\theta}(x), \forall x \in K. \end{aligned}$$

To define a loss function, we combine the conditions formulated above and write, with respect to the map $T_{\theta}(x)$:

$$\mathcal{L}_2(w, p, x) = [2w^\top pf(x) + w^\top w]_+ - [w^\top w]_-, \quad (9)$$
where $w$ represents $T_\theta$ and $p$ corresponds to the derivative $DT_\theta$ with respect to $x$. Note that both terms in $\mathcal{L}_2$ are not squares of the argument in comparison with the formulation of [8]. Therefore, it results in a simpler gradient with respect to $\theta$. Moreover, according to the results of [9], such a presentation of a Lyapunov function always exists for $n \neq 4, 5$ (see Proposition 1).

Note that in such a case the Lyapunov function candidate (8) is a smooth function on $\mathbb{R}^n$, while in Proposition 1 the smoothness at the origin is not ensured for the studied class of systems in (1). Therefore, it is a mild limitation of this ANN approximation, which has however important consequences. Indeed, following the result of Theorem 5, for any $\varepsilon > 0$ and any compact set $K$, there is $N_1 > 0$ such that for a vector of parameters $\theta$ properly tuned by a learning algorithm, $\sup_{x \in K} \|V(x) - W(x; \theta)\| \le \varepsilon$, for some Lyapunov function $V$ for the system (1). However, the function $V$ vanishes at the origin by Definition 1, while the error function $V(x) - W(x; \theta)$, being upper bounded in the norm by $\varepsilon$ in the set $K$, does not necessarily vanishes around $x = 0$. The following worse-case estimates illustrate this idea, which connects the accuracy on a grid (measured after training) with the total accuracy in $K$:

**Lemma 1.** If the transformation map $T$ (with Lipschitz constant $\rho > 0$ in $K$: $\|T(x) - T(x')\| \le \rho\|x - x'\|$ for all $x, x' \in K$) is approximated by $T_\theta(x)$ in (7) on a grid $\{x_i\}_{i=1,N} \subset K$, with accuracy $\varepsilon > 0$, $\max_{i=1,N} \|T_\theta(x_i) - T(x_i)\| \le \varepsilon$, with maximal step $\Delta = \max_{i=1,N} \min_{j=1,N, j \neq i} \|x_i - x_j\|$, and uniformly sampled in $K$: $\sup_{\xi \in K} \min_{i=1,N} \|\xi - x_i\| \le \frac{\Delta}{2}$, then

$$ \|T_\theta(\xi) - T(\xi)\| \le (\rho + a\|w^2\|\|w^1\|) \frac{\Delta}{2} + \varepsilon \\ \text{for any } \xi \in K, \text{ where } \sup_{s \in \mathbb{R}} |\sigma'_1(s)| \le a. $$

The proof of this Lemma is skipped for the sake of brevity. A consequence of this result is that, if increasing the number of neurons in the hidden layer $N_1$ (decreasing $\varepsilon$) and increasing the density of the grid (i.e., decreasing $\Delta$), it is possible during the training to keep small the amplitudes of the ANN weights $\|w^2\|$ and $\|w^1\|$. Then any desired accuracy of approximation of $T(x)$ can be achieved by adjusting $N_1$ and $\Delta$. The same result can be demonstrated for approximation of $DT(x)$ for $x \in K$ in terms of a uniform upper bound on the error.

Therefore, $W(x; \theta)$ and $DW(x; \theta)f(x)$ could be sign-indefinite close to the origin, and there is no theoretical guarantee of having a solution to this drawback. Additional researches in this directions have a good perspective.

### E. Approximation for homogeneous systems

Thus, for the system (1) in its generality, the behavior at the origin of any Lyapunov function approximation can be a significant obstacle to ANN/PINN application. Fortunately, this obstruction can be avoided for homogeneous systems, where a homogeneous Lyapunov function can be selected (Theorem 2), whose design and analysis can be performed on a sphere only, leading to a global solution. Hence, the idea is to use a compact set that excludes the origin, estimate the Lyapunov function and then use dilations to produce a global expression.

For homogeneous systems, we consider an approximation $V_\theta$ of an unknown global Lyapunov function $V$, and where ANN unit $W(x; \theta)$ defined in (6) is needed on the set $S_r(1)$ only:

$$ V_{\theta}(x) = \|x\|_{\mathbf{r}}^{\mu} W(\Lambda_{\mathbf{r}}^{-1}(\|x\|_{\mathbf{r}}) x; \theta), \quad \forall x \in S_{\mathbf{r}}(1), \quad (10) $$

where $\mu > 0$ is a design parameter. Such a Lyapunov function $V_\theta$ is homogeneous by construction. Then the conditions to check can be obtained from the estimates (4), (5) and restricted to $S_r(1)$:

$$ DW(y; \theta)f(y) \le -1, \quad W(y; \theta) \ge 1, \quad \forall y \in S_r(1), $$

which can be transformed into the following loss function: $\mathcal{L}_r(w, p, y) = [pf(y) + 1]_+ - [w - 1]_-$.

It is easy to see that, in this case, there is no problem with the vanishing of $W$ somewhere on the set of approximation $S_r(1)$. Moreover, for homogeneous systems and the chosen form of Lyapunov function approximation (10), there exist theoretical guarantees coming from Theorem 5, contrarily the general case of the system (1) discussed in the previous section (10) may be not smooth at the origin):

**Theorem 6** ([4]). There are a number of hidden units $N_1 > 0$ and a vector of parameters $\theta$ such that (10) is an r-homogeneous global Lyapunov function for any degree $\mu > \max\{0, -\nu\}$ for an r-homogeneous system (1) of degree $\nu$ with GAS origin.

Following the result of Lemma 1, increasing the number of neurons in the hidden layer $N_1$ and density of the grid can guarantee the required accuracy of approximation.

Comparing this results and its applicability conditions with the corresponding ones from [6], it is easy to conclude that ANN-based approach from Theorem 6 is much simpler. In addition, it allows many extensions to be straightforwardly obtained, as coupling with a stabilizing control design presented in [4], or also observer synthesis.

### III. MAIN RESULT

To summarize the previously given material, for approximation of Lyapunov functions of general nonlinear systems, the existing theory cannot guarantee that close to the origin the obtained ANN-based function is positive definite and has a negative definite derivative. This drawback can be avoided for homogeneous systems, where derivation of a global Lyapunov function is possible according to Theorem 6.

Unfortunately not all systems are homogeneous, but many of them may have local homogeneous approximations, then in this work we are going to study such a class of systems in order to approximate a Lyapunov function on a given compact set $K \subset \mathbb{R}^n$. Our idea is to use a Lyapunov function candidate derived for the homogeneous approximation close to the origin, while a general NN-based structure (8) will be applied extending the domain of validity of the former one (out a vicinity of zero).

Hence, we start with assumption of existence of a homogeneous approximation.

**Assumption 1.** The system (1) is $(\mathbf{r}, 0, f_0)$-homogeneous with a continuous function $f_0: \mathbb{R}^n \to \mathbb{R}^n$, and the system

$$ \dot{x} = f_0(x) \quad (11) $$

is asymptotically stable.

Note that by construction $f_0$ is r-homogeneous of some degree $\nu_0 \ge -r_{\min}$. According to Theorem 3, the system (1) is locally asymptotically stable in this case, while according to Theorem 6, there exist $N_1 > 0$ and properly tuned weights $\theta_1$ such that (10) (we will denote it as $V_{\theta_1}(x)$) is a r-homogeneous Lyapunov function of some degree $\mu > \max\{0, -\nu_0\}$ for (11) (i.e., the properties of Theorem 2 are satisfied for $V_{\theta_1}(x)$ from (10) with the system (11)).
Denote by $\kappa > 0$ a level of $V_{\theta_1}$, where it is also a Lyapunov function for the original system (1):

$$DV_{\theta_1}(x)f(x) < 0, \forall x \in \mathcal{S}_{\kappa} = \{x \in \mathbb{R}^n : V_{\theta_1}(x) = \kappa\},$$

denote the respective ball by $\mathcal{B}_{\kappa} = \{x \in \mathbb{R}^n : V_{\theta_1}(x) \le \kappa\}$. Such a level $\kappa$ always exist by the properties of homogeneous approximations and continuity of $V_{\theta}$, $f$ and $f_0$. We assume that $\mathcal{B}_{\kappa} \subset \text{int} K$, where $K \subset \mathbb{R}^n$ is the set of interest, and it is required to find a Lyapunov function for (1) on the set $K$. Therefore, for $K \setminus \mathcal{B}_{\kappa}$ we are going to apply the Lyapunov function structure $W(x; \theta)$ from (8) (where the number of neurons in the hidden layer for the ANN approximating $T_{\theta}$ is denoted now by $N'_1 > 0$), where since the origin is excluded, the minimization of the loss (9) for all $x \in K \setminus \mathcal{B}_{\kappa}$ is feasible. Note also that we need to guarantee at least continuity of the combined Lyapunov function. To this end, we may introduce an additional boundary condition in order to unite $V_{\theta_1}(x)$ and $W(x; \theta)$ on the set $\mathcal{S}_{\kappa}$, or we can just extend the domain of approximation for $W(x; \theta)$ to $K \setminus \mathcal{B}_{\kappa/2}$ developing (9) to

$$\mathcal{L}_{\kappa}(w, p, x) = \left[ 2w^{\top}pf(x) + w^{\top}w \right]_{+} + \begin{cases} \left(w^{\top}w - \frac{\kappa}{4}\right)^2 & x \in \mathcal{S}_{\kappa/2} \\ -[w^{\top}w]_{-} & x \in \mathcal{B}_{\kappa} \setminus \mathcal{B}_{\kappa/2}, \\ -[w^{\top}w - \kappa]_{-} & x \in K \setminus \mathcal{B}_{\kappa} \end{cases}$$

and define a continuous Lyapunov function for (1) on $K$ as follows:

$$V_{\theta_1, \theta}(x) = \begin{cases} V_{\theta_1}(x) & x \in \mathcal{B}_{\kappa/2} \\ \max\{V_{\theta_1}(x), \chi W(x; \theta)\} & x \in K \setminus \mathcal{B}_{\kappa/2} \end{cases}, \quad (12)$$

where the constant $\chi > 0$ is chosen in a way to guarantee that $\chi W(x; \theta) \ge V_{\theta_1}(x)$ for all $x \in K \setminus \mathcal{B}_{\kappa}$. Such a coefficient $\chi$ exists due to introduced modification in the loss functional $\mathcal{L}_{\kappa}$, and we have the following properties of the proposed ANN-based Lyapunov function:

**Theorem 7.** Under Assumption 1, there exists $N'_1 > 0$ and $N'_1 > 0$ sufficiently large, $\chi > 0$ such that (12) is a locally Lipschitz continuous (piecewise continuously differentiable) Lyapunov function for (1) on the domain of attraction of the origin $K$.

For the sake of brevity, we skip the proof of the Theorem. Remark that we do not require that the derivatives of $V_{\theta_1}$ and $W$ coincide, then the final combined Lyapunov function is just locally Lipschitz continuous, and not continuously differentiable in $K$.

**Remark 1.** Replacing Assumption 1 with requirement of existence of a homogeneous approximation at infinity, and using the same ideas with Theorem 4, we can also construct a Lyapunov function for global practical asymptotic stability (global asymptotic stability of the set), then formulation of this result we left for a future work. In such a setting, application of the homogeneity framework allows us to design an ANN-based Lyapunov function on a non-compact domain.

## IV. CONCLUSION

Based on fundamental existing results, in this work, we developed an approach to ANN-based design of Lyapunov functions on a given compact set for asymptotically stable autonomous dynamical systems at the origin. To avoid losses of positive definiteness of an approximated by ANN Lyapunov function, and shortages of negative definiteness of its derivative close to the origin, we proposed to combine two Lyapunov functions: one derived in a

vicinity of the origin for a homogeneous approximation of the system (e.g., for the linearization) and another one based on the quadratic canonical form. The former does not possess the indicated drawback, while the latter can be rigorously constructed out of a neighborhood of zero. Future research will include the development of the proposed approach to simultaneous ANN-based control and observer design.

## REFERENCES

[1] V. Andrieu, L. Praly, and A. Astolfi. Homogeneous Approximation, Recursive Observer Design, and Output Feedback. *SIAM Journal of Control and Optimization*, 47(4):1814–1850, 2008.

[2] P.-H. Chiu, J. C. Wong, C. Ooi, M. H. Dao, and Y.-S. Ong. CAN-PINN: A fast physics-informed neural network based on coupled-automatic-numerical differentiation method. *Computer Methods in Applied Mechanics and Engineering*, 395:114909, 2022.

[3] G. Cybenko. Approximation by superpositions of a sigmoidal function. *Mathematics of Control, Signals and Systems*, 2(4):303–314, Dec 1989.

[4] D. R. de Lima, R. Ushirobira, and D. Efimov. A physics-informed neural network method to approximate homogeneous Lyapunov functions. In *Proc. European Control Conference (ECC)*, Stockholm, 2024.

[5] M. W. M. G. Dissanayake and N. Phan-Thien. Neural-network-based approximations for solving partial differential equations. *Communications in Numerical Methods in Engineering*, 10(3):195–201, 1994.

[6] D. Efimov, R. Ushirobira, J. A. Moreno, and W. Perruquetti. Homogeneous Lyapunov functions: from converse design to numerical implementation. *SIAM Journal on Control and Optimization*, 56(5):3454–3477, 2018.

[7] S. Ferrari and R.F. Stengel. Smooth function approximation using neural networks. *IEEE Transactions on Neural Networks*, 16(1):24–38, 2005.

[8] L. Grüne. Computing Lyapunov functions using deep neural networks. *Journal of Computational Dynamics*, 8(2):131–152, 2021.

[9] L. Grüne, E. D. Sontag, and F. R. Wirth. Asymptotic stability equals exponential stability, and ISS equals finite energy gain-if you twist your eyes. *Systems & Control Letters*, 38(2):127–134, 1999.

[10] H. Hermes. Nilpotent and high-order approximations of vector field systems. *SIAM review*, 33(2):238–264, 1991.

[11] K. Hornik, M. Stinchcombe, and H. White. Multilayer feedforward networks are universal approximators. *Neural Networks*, 2(5):359–366, 1989.

[12] H. K. Khalil. *Nonlinear Systems*. Pearson Education. Prentice Hall, 2002.

[13] Y. Long and M. M. Bayoumi. Feedback stabilization: control Lyapunov functions modelled by neural networks. In *Proceedings of 32nd IEEE Conference on Decision and Control*, pages 2812–2814 vol.3, 1993.

[14] V. Petridis and S. Petridis. Construction of neural network based Lyapunov functions. In *The 2006 IEEE International Joint Conference on Neural Network Proceedings*, pages 5059–5065, 2006.

[15] R. C. B. Rego and F. M. U. de Araújo. Learning-based robust neuro-control: A method to compute control Lyapunov functions. *International Journal of Robust and Nonlinear Control*, 32(5):2644–2661, 2022.

[16] I. D. Jimenez Rodriguez, A. Ames, and Y. Yue. LyaNet: A Lyapunov framework for training neural ODEs. In *Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning*, volume 162 of *Proceedings of Machine Learning Research*, pages 18687–18703. PMLR, 17–23 Jul 2022.

[17] L. Rosier. Homogeneous Lyapunov function for homogeneous continuous vector field. *Systems & Control Letters*, 19(6):467–473, 1992.

[18] V. I. Zubov. On systems of ordinary differential equations with generalized homogenous right-hand sides. *Izvestia vuzov. Mathematica.*, 1:80–88, 1958. in Russian.

[19] V. I. Zubov. *Methods of A. M. Lyapunov and Their Application*. P. Noordhoff Ltd, Groningen, 1964.