<figure><img src="image_1.png" /></figure>

HAL
open science

# A guaranteed numerical method to prove the exponential stability of nonlinear discrete-time systems

Morgan Louedec, Luc Jaulin, Christophe Viel

<figure><img src="image_2.png" /></figure>

## To cite this version:

Morgan Louedec, Luc Jaulin, Christophe Viel. A guaranteed numerical method to prove the exponential stability of nonlinear discrete-time systems. 2024. hal-04465206

HAL Id: hal-04465206

https://hal.science/hal-04465206v1

Preprint submitted on 19 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A guaranteed numerical method to prove the exponential stability of nonlinear discrete-time systems

Morgan Louédec, Luc Jaulin, Christophe Viel

*Abstract*—This paper presents a guaranteed numerical method for proving the exponential stability of an n-dimensional nonlinear discrete-time system. This method also finds positive invariant ellipsoids. It is based on Lyapunov theory, set theory, interval analysis, and guaranteed ellipsoidal propagation. The method is also computationally tractable and can be used on high-dimensional systems where Lyapunov functions are difficult to find.

*Index Terms*—Discrete-time, Ellipsoids, Nonlinear, Numerical method, Positive Invariant, Stability Analysis

## I. INTRODUCTION

For nonlinear systems, the equilibrium point stability is usually studied by examining the eigenvalues of the Jacobian matrix and then searching for the region of attraction of the equilibrium point using Lyapunov functions, as presented in [2], [6] and applied in [10], [12], [13]. However, Lyapunov functions are difficult to find for complex high-dimensional system, such as in [11].

An alternative method can be found in set theory with the notion of positive invariance [1]. When the trajectory of the system enters a positive invariant set, it will remain in it. Set theory can be numerically implemented with guaranteed computation, especially on intervals, as pretended in [7], [8]. Ellipsoids are often positive invariant because their contour is a level set of a quadratic Lyapunov function.

In this paper, we propose a guaranteed numerical set method to prove the exponential stability of a n-dimensional nonlinear discrete-time system and to find a positive invariant ellipsoid.

This method first solves a discrete Lyapunov equation to find an ellipsoid that is likely to be positive invariant. From this ellipsoid, an ellipsoidal enclosure of the reachable set is computed using guaranteed propagation method. If this enclosure is strictly included in the initial ellipsoid, then the ellipsoid is positive invariant with respect to the system. In addition, the system is exponentially stable on the ellipsoid. If the inclusion is not verified, the process is repeated by shrinking the initial ellipsoid, which reduces the pessimism caused by the non-linearity and makes the inclusion more likely to be verified.

This method assumes that the evolution function of the discrete system is continuously differentiable and that the Jacobian matrix is not singular. It can be applied without knowing the analytical expression of the evolution function

of the discrete system as long as one can enclose the image of the inclusion function of the Jacobian matrix. This numerical method is also computationally tractable and can therefore be applied to high-dimensional problems for which Lyapunov functions are difficult to find.

Section II presents some notions about discrete systems and ellipsoidal propagation. The stability problem is described in Section III. Section IV presents some key theorems to prove stability using ellipsoidal propagation. Section V presents a method to find ellipsoids that are likely to be positive invariant. Section VI presents the guaranteed numerical method to prove stability. Finally Section VII presents an example to illustrate the method.

## II. NOTATION AND HYPOTHESIS

This paper deals with the stability of a discrete nonlinear system described by

$$ 
\begin{aligned}
x_{k+1} &= h(x_k), \\
0 &= h(0)
\end{aligned} 
\quad (1) 
$$

with the state vector $x_k \in \mathbb{R}^n$ and where $h$ is $C^1$. The origin $0$ of $\mathbb{R}^n$ is an equilibrium point of the system. Given a set $S \subseteq \mathbb{R}^n$, the image of this set by the mapping $h$ is written $h(S) \subseteq \mathbb{R}^n$. The Jacobian matrix of $h$ at the origin is written

$$ J = \frac{dh}{dx}(0), \quad (2) $$

and is assumed invertible.

In this paper, ellipsoids are used to study the system. Let $S_n^+$ be the set of all the real symmetric positive definite matrices. For each matrix $Q \in S_n^+$, there exists a unique square root matrix $P \in S_n^+$ such that $Q = P^2$. One can write $P = Q^{1/2}$. Moreover, the $Q$-norm is given by

$$ \forall x \in \mathbb{R}^n, \|x\|_Q = \sqrt{x^T Q x}. \quad (3) $$

Positive symmetric matrices can define ellipsoids as described in Definition 1 and as presented in [1, Chapter 3.2].

**Definition 1.** A non-degenerated centered ellipsoid is a subset of $\mathbb{R}^n$, described by a unique shape matrix $\Gamma \in S_n^+$ such that

$$ \mathcal{E}(\Gamma) := \{x \in \mathbb{R}^n \mid \|x\|_{\Gamma^{-2}} \le 1\}. \quad (4) $$

An ellipsoid can also be written $\mathcal{E}$ with a shape matrix $\Gamma$ such that $\mathcal{E} = \mathcal{E}(\Gamma)$ and $\Gamma = \mathcal{E}^{-1}(\mathcal{E})$. The border of the ellipsoid $\mathcal{E}$ is written $\partial\mathcal{E}$. For all, $\Gamma \in S_n^+$, all real scalar

This work has been supported by the French Defence Innovation Agency (AID) and the Brittany Region
$\alpha \in \mathbb{R}$ with $\alpha \neq 0$, and all invertible matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$ one has

$$ \alpha \cdot \mathcal{E}(\mathbf{\Gamma}) = \mathcal{E}(\alpha \cdot \mathbf{\Gamma}), \qquad (5) $$

$$ \mathbf{A} \cdot \mathcal{E}(\mathbf{\Gamma}) = \mathcal{E} \left( \left( \mathbf{A} \mathbf{\Gamma}^2 \mathbf{A}^T \right)^{\frac{1}{2}} \right). \qquad (6) $$

For two ellipsoids $\mathcal{E}_1$ and $\mathcal{E}_2$, the inclusion and the strict inclusion are respectively denoted by the symbols $\subseteq$ and $\subset$ such that

$$ (\mathcal{E}_1 \subset \mathcal{E}_2) \Leftrightarrow (\mathcal{E}_1 \subseteq \mathcal{E}_2 \text{ and } \mathcal{E}_1 \cap \partial \mathcal{E}_2 = \emptyset). \qquad (7) $$

As presented in [3], with $\mathbf{\Gamma}_1 = \mathcal{E}^{-1}(\mathcal{E}_1)$ and $\mathbf{\Gamma}_2 = \mathcal{E}^{-1}(\mathcal{E}_2)$, one has the following equivalence

$$ (\mathcal{E}_1 \subset \mathcal{E}_2) \Leftrightarrow \left( (\mathbf{\Gamma}_1^{-2} - \mathbf{\Gamma}_2^{-2}) \in \mathcal{S}_n^+ \right). \qquad (8) $$

Ellipsoids can be positive invariant with respect to the discrete system as described in Definition 2. It means that if the state of the system enters the ellipsoid $\mathcal{E}$, then it will remain in it.

**Definition 2.** (Positive Invariance) An ellipsoid $\mathcal{E}$ of $\mathbb{R}^n$ is positive invariant with respect to the discrete system (1) if and only if

$$ \mathbf{h}(\mathcal{E}) \subseteq \mathcal{E} \qquad (9) $$

Using the norm from (4), the inclusion (9) is equivalent to

$$ \forall \mathbf{x}_0 \in \mathcal{E}, \quad \|\mathbf{h}(\mathbf{x}_0)\|_{\mathbf{\Gamma}^{-2}} \le 1, \qquad (10) $$

with $\mathbf{\Gamma} = \mathcal{E}^{-1}(\mathcal{E})$. When the mapping $\mathbf{h}$ is linear, positive invariant ellipsoids are found by solving a Lyapunov equation as described in Theorem 1 and presented in [2, Section 3].

**Theorem 1.** Let $\mathbf{J} \in \mathbb{R}^n$ be a Schur matrix and let $\mathbf{Q} \in \mathcal{S}_n^+$. Let $\mathbf{P} \in \mathcal{S}_n^+$ be the unique solution to the discrete Lyapunov equation

$$ \mathbf{A}^T \mathbf{P} \mathbf{A} - \mathbf{P} = -\mathbf{Q}, \qquad (11) $$

and let $\mathbf{\Gamma} = \mathbf{P}^{-\frac{1}{2}}$. Then, the ellipsoid $\mathcal{E}(\mathbf{\Gamma})$ is positive invariant with respect to the linear mapping $\mathbf{x}_{k+1} = \mathbf{J} \cdot \mathbf{x}_k$.

To perform this study, some tools from Interval analysis are also used. For further information on interval analysis, one may refer to [7]. An interval $[x]$ is defined as a connected subset of $\mathbb{R}$. The function $\mathrm{ub}(\cdot)$ gives the upper bound of an interval. A box $[x]$ is defined as a subset of $\mathbb{R}^n$ and a Cartesian product of $n$-intervals. The set of boxes over $\mathbb{R}^n$ is written $\mathbb{R}^n$. The Unit box is written

$$ [1_n] = [-1, 1]^n. \qquad (12) $$

An $(m \times n)$-dimensional interval matrix $[\mathbf{X}]$ is a subset of $\mathbb{R}^{n \times m}$ that can be defined as the Cartesian product of $mn$ intervals.

Considering the function $\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m$, the interval function $[\mathbf{f}]: \mathbb{R}^n \to \mathbb{R}^m$ is an inclusion function for $\mathbf{f}$ if

$$ \forall [\mathbf{x}] \in \mathbb{R}^n, \mathbf{f}([\mathbf{x}]) \in [\mathbf{f}]([\mathbf{x}]), \qquad (13) $$

where $\mathbf{f}([\mathbf{x}]) = \{\mathbf{f}(\mathbf{x}) | \mathbf{x} \in [\mathbf{x}]\}$ is the image of the interval $[\mathbf{x}]$ by the function $\mathbf{f}$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Propagation of ellipsoids</figcaption></figure>

Then, the evolution of ellipsoids with the evolution function $\mathbf{h}$ can be studied using outer ellipsoidal enclosures as defined by Definition 3.

**Definition 3.** (Ellipsoidal outer enclosure) Let $\mathcal{E}$ be an ellipsoid of $\mathbb{R}^n$ and let $\mathbf{h}$ be a mapping over $\mathbb{R}^n$. The ellipsoid $\mathcal{E}_{\mathrm{out}}$ is an outer enclosure of the image set $\mathbf{h}(\mathcal{E})$ if it verifies $\mathbf{h}(\mathcal{E}) \subseteq \mathcal{E}_{\mathrm{out}}$.

A method to compute outer ellipsoidal enclosures is presented in Theorem 2 and illustrated in Figure 1. This Theorem is proved in [9].

**Theorem 2.** (Ellipsoidal Propagation) Let $\mathcal{E}$ be an ellipsoid of $\mathbb{R}^n$, let $\mathbf{h}$ be a $C^1$ mapping over $\mathbb{R}^n$ that verifies $\mathbf{0} = \mathbf{h}(\mathbf{0})$, let $\mathbf{\Gamma} = \mathcal{E}^{-1}(\mathcal{E})$, let $\mathbf{J}$ be the Jacobian matrix of $\mathbf{h}$ at the origin and let $[\mathbf{J}]$ be an inclusion function for $\frac{d\mathbf{h}}{dx}$ over $\mathbb{R}^n$. Consider the matrix $\mathbf{\Gamma}_{\mathrm{out}}$ given by

$$ \mathbf{\Gamma}_{\mathrm{out}} = (1 + \rho) \left( \mathbf{J} \mathbf{\Gamma}^2 \mathbf{J}^T \right)^{\frac{1}{2}} \qquad (14) $$

with the inflation gain

$$ \rho = \mathrm{ub} \left( \| (\mathbf{\Gamma}^{-1} \cdot \mathbf{J}^{-1} \cdot [\mathbf{J}]([\mathcal{E}]) \cdot \mathbf{\Gamma} - \mathbf{I}_n) \cdot [1_n] \| \right), \qquad (15) $$

where the box $[\mathcal{E}]$ that enclose $\mathcal{E}$ is given by

$$ [\mathcal{E}] = \begin{bmatrix} \|\mathbf{\Gamma}_1\|_2 & 0 & \cdots & 0 \\ 0 & \|\mathbf{\Gamma}_2\|_2 & & \vdots \\ \vdots & & \ddots & \vdots \\ 0 & \cdots & 0 & \|\mathbf{\Gamma}_n\|_2 \end{bmatrix} \cdot [1_n]. \qquad (16) $$

where $\mathbf{\Gamma}_i$ is the $i$-th column of $\mathbf{\Gamma}$. Then, the ellipsoid $\mathcal{E}(\mathbf{\Gamma}_{\mathrm{out}})$ is an outer outer enclosure of the set $\mathbf{h}(\mathcal{E})$.

Finally, let $\mathcal{P}_h$ be the operator that computes the outer enclosure of $\mathbf{h}(\mathcal{E})$ such that $\mathcal{E}_{\mathrm{out}} = \mathcal{P}_h(\mathcal{E})$.

### III. PROBLEM DESCRIPTION

This paper proposes a method to find an ellipsoid $\mathcal{E}$ that is positive invariant with respect to the system (1) and such that this system is exponentially stable in $\mathcal{E}$. In other words, one must prove that there exist $\alpha \in ]0, 1]$ and $\beta > 0$ and an ellipsoid $\mathcal{E}$ that verify

$$ \mathbf{x}_0 \in \mathcal{E} \Rightarrow \|\mathbf{x}_k\|_{\mathbf{\Gamma}^{-2}} \le \alpha \|\mathbf{x}_0\|_{\mathbf{\Gamma}^{-2}} e^{-\beta k}, \forall k \in \mathbb{N}, \qquad (17) $$

with $\mathbf{\Gamma} = \mathcal{E}^{-1}(\mathcal{E})$. Note that $\alpha < 1$ to respect the requirement for positive invariance from (10).
Thus the problem can be decomposed in two steps:

1) Having a sufficient condition for (17) that can be numerically verified, as presented in Section IV.

2) Finding an ellipsoid $\mathcal{E}$ that is likely positive invariant, as presented in Section V.

#### IV. PROOF OF THE EXPONENTIAL STABILITY

A sufficient condition for (17) is described in Theorem 3. The inclusion (18) can be numerically verified as presented in Section VI.

**Theorem 3.** Consider the discrete system (1) with the propagation operator $\mathcal{P}_h$. If there exist an ellipsoid $\mathcal{E}$ such that

$$ \mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}, \quad (18) $$

then $\mathcal{E}$ is positive invariant with respect to the system (1) and the system is exponentially stable on the ellipsoid $\mathcal{E}$.

The proof of Theorem 3 can be found in the Appendix. Theorem 3 is deduced from the Lemma 1 and the Theorems 4 and 5 described below. It can be summarized as follows.

If the propagation operator $\mathcal{P}_h$ contracts the ellipsoid $\mathcal{E}$ by a factor $\alpha$, then the mapping $\mathbf{h}$ contracts every shortened ellipsoid $\alpha \cdot \mathcal{E}$ with the same factor, as described in Theorem 4 and illustrated by Figure 2.

Then, multiple occurrences of the mapping $\mathbf{h}$ results in an exponential contraction, as described in Theorem 5 and illustrated by Figure 3.

Finally, from this exponential contraction, the exponential stability of the system is deduced, as described in Theorem 3.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of Theorem 4</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Illustration of Theorem 5</figcaption></figure>

**Lemma 1.** Consider the discrete system (1) with the propagation operator $\mathcal{P}_h$. For every real $\alpha \in ]0, 1]$ and every ellipsoid $\mathcal{E}$, one has

$$ \mathcal{P}_h (\alpha \cdot \mathcal{E}) \subseteq \alpha \cdot \mathcal{P}_h (\mathcal{E}). \quad (19) $$

*Proof.* Consider the ellipsoid $\mathcal{E}$, the matrix $\mathbf{\Gamma} = \mathcal{E}^{-1}(\mathcal{E})$ and $\alpha \in ]0, 1]$. From (5), one has

$$ \alpha \cdot \mathcal{E} = \mathcal{E} (\alpha \cdot \mathbf{\Gamma}). \quad (20) $$

Then, from Theorem 2, one has

$$ \left\{ \begin{array}{l} \mathcal{P}_h (\mathcal{E} (\alpha \cdot \mathbf{\Gamma})) = \mathcal{E} (\mathbf{\Gamma}_{\text{out},\alpha}), \\ \mathcal{P}_h (\mathcal{E} (\mathbf{\Gamma})) = \mathcal{E} (\mathbf{\Gamma}_{\text{out}}), \end{array} \right. \quad (21) $$

with the matrices

$$ \left\{ \begin{array}{l} \mathbf{\Gamma}_{\text{out},\alpha} = (1 + \rho_{\alpha}) \left( \mathbf{J} (\alpha \mathbf{\Gamma})^2 \mathbf{J}^T \right)^{\frac{1}{2}}, \\ \mathbf{\Gamma}_{\text{out}} = (1 + \rho) \left( \mathbf{J} \mathbf{\Gamma}^2 \mathbf{J}^T \right)^{\frac{1}{2}}, \end{array} \right. \quad (22) $$

with the inflation gains

$$ \left\{ \begin{array}{l} \rho_{\alpha} = \text{ub} (\| (\mathbf{M}_{\alpha} - \mathbf{I}_n) \cdot [1_n] \|), \\ \rho = \text{ub} (\| (\mathbf{M} - \mathbf{I}_n) \cdot [1_n] \|), \end{array} \right. \quad (23) $$

with the interval matrices

$$ \begin{aligned} \mathcal{M}_{\alpha} &= (\alpha \mathbf{\Gamma})^{-1} \cdot \mathbf{J}^{-1} \cdot [\mathbf{J}] ([\mathcal{E}_{\alpha}]) [\mathbf{J}]_{\alpha} \cdot (\alpha \mathbf{\Gamma}), \\ &= \mathbf{\Gamma}^{-1} \cdot \mathbf{J}^{-1} \cdot [\mathbf{J}] ([\mathcal{E}_{\alpha}]) [\mathbf{J}]_{\alpha} \cdot \mathbf{\Gamma}, \end{aligned} \quad (24) $$

and

$$ \mathcal{M} = \mathbf{\Gamma}^{-1} \cdot \mathbf{J}^{-1} \cdot [\mathbf{J}] ([\mathcal{E}]) \cdot \mathbf{\Gamma}, \quad (25) $$

where $[\mathcal{E}_{\alpha}]$ and $[\mathcal{E}]$ are given by (16). From (16), one has $[\mathcal{E}_{\alpha}] = \alpha \cdot [\mathcal{E}]$. Since $\alpha \in ]0, 1]$, one deduces

$$ [\mathcal{E}_{\alpha}] \subseteq [\mathcal{E}]. \quad (26) $$

Thus one gets

$$ \mathcal{M}_{\alpha} \subseteq \mathcal{M}, \quad (27) $$

which, using (23), results in an inequality of the upper bound, $\rho_{\alpha} \le \rho$. Then, from (22) one also has

$$ \mathbf{\Gamma}_{\text{out},\alpha} = \frac{\alpha \cdot (1 + \rho_{\alpha})}{(1 + \rho)} \cdot \left( (1 + \rho) \left( \mathbf{J} \mathbf{\Gamma}^2 \mathbf{J}^T \right)^{\frac{1}{2}} \right), \quad (28) $$

$$ = \gamma \cdot \mathbf{\Gamma}_{\text{out}}, \quad (29) $$

with

$$ \gamma = \frac{\alpha \cdot (1 + \rho_{\alpha})}{(1 + \rho)}. \quad (30) $$

Therefore,

$$ \begin{aligned} \mathcal{P}_h (\alpha \cdot \mathcal{E}) &= \mathcal{P}_h (\mathcal{E} (\alpha \cdot \mathbf{\Gamma})) \\ &= \mathcal{E} (\mathbf{\Gamma}_{\text{out},\alpha}) \\ &= \mathcal{E} (\gamma \cdot \mathbf{\Gamma}_{\text{out}}) \\ &= \gamma \cdot \mathcal{E} (\mathbf{\Gamma}_{\text{out}}) \\ &= \gamma \cdot \mathcal{P}_h (\mathcal{E}) \end{aligned} \quad (31) $$
Moreover, since $\alpha \in ]0, 1]$ and $\rho_\alpha \le \rho$, one gets $\gamma \in ]0, \alpha]$. As a result,

$$ P_h(\alpha \cdot \mathcal{E}) \subseteq \alpha \cdot P_h(\mathcal{E}). \quad (32) $$

**Theorem 4.** Consider the ellipsoid $\mathcal{E}$ and the discrete system (1) with the propagation operator $P_h$. If there exists a real $\gamma \in ]0, 1[$ such that

$$ P_h(\mathcal{E}) \subset \gamma \cdot \mathcal{E}, \quad (33) $$

then, for all $\alpha \in ]0, 1]$ the ellipsoid $\mathcal{E}_\alpha = \alpha \cdot \mathcal{E}$ verifies

$$ h(\mathcal{E}_\alpha) \subset \gamma \cdot \mathcal{E}_\alpha. \quad (34) $$

*Proof.* Let $\gamma \in ]0, 1[, \alpha \in ]0, 1]$ and assume that (33) is true. Since $P_h(\mathcal{E}_\alpha)$ is an outer enclosure of $h(\mathcal{E}_\alpha)$, one has

$$ \begin{aligned} h(\mathcal{E}_\alpha) &\subseteq P_h(\mathcal{E}_\alpha), \\ &\subseteq P_h(\alpha \cdot \mathcal{E}). \end{aligned} \quad (35) $$

Moreover, since $\alpha \in ]0, 1]$, using Lemma 1, one gets

$$ P_h(\alpha \cdot \mathcal{E}) \subseteq \alpha \cdot P_h(\mathcal{E}). \quad (36) $$

Then, from (33), one deduces

$$ \begin{aligned} \alpha \cdot P_h(\mathcal{E}) &\subseteq \alpha \cdot \gamma \cdot \mathcal{E}, \\ &\subseteq \gamma \cdot \mathcal{E}_\alpha. \end{aligned} \quad (37) $$

So, by regrouping (35), (36) and (37), (34) is verified. □

**Theorem 5.** Consider the ellipsoid $\mathcal{E}$ and the discrete system (1) with the propagation operator $P_h$. If there exist a real $\gamma \in ]0, 1[$ such that (33) is true, then for all $\alpha \in ]0, 1]$ and all integer $k \ge 1$, the ellipsoid $\mathcal{E}_\alpha = \alpha \cdot \mathcal{E}$ verifies

$$ h^k(\mathcal{E}_\alpha) \subset \gamma^k \mathcal{E}_\alpha. \quad (38) $$

*Proof.* Let $\gamma \in ]0, 1[, \alpha \in ]0, 1]$ and assume that (33) is true. Let us prove (38) by recurrence on $k$. The case $k = 1$ is verified by Theorem 4. Then, let $k \ge 1$ and assume that

$$ h^k(\mathcal{E}_\alpha) \subset \gamma^k \mathcal{E}_\alpha. \quad (39) $$

Let us show that (38) holds for $k + 1$. Applying the mapping $h$ on (39) results in

$$ h^{k+1}(\mathcal{E}_\alpha) \subseteq h(\gamma^k \mathcal{E}_\alpha), \quad (40) $$

Then, since $P_h(\gamma^k \mathcal{E}_\alpha)$ is an outer enclosure of $h(\gamma^k \mathcal{E}_\alpha)$, one also has

$$ \begin{aligned} h(\gamma^k \mathcal{E}_\alpha) &\subseteq P_h(\gamma^k \mathcal{E}_\alpha) \\ &\subseteq P_h(\gamma^k \alpha \cdot \mathcal{E}) \end{aligned} \quad (41) $$

Moreover, since $\gamma^k \alpha \in ]0, 1[$, using Lemma 1, one deduces

$$ P_h(\gamma^k \alpha \cdot \mathcal{E}) \subseteq \gamma^k \alpha P_h(\mathcal{E}). \quad (42) $$

In addition, from (33), one gets

$$ \begin{aligned} \gamma^k \alpha P_h(\mathcal{E}) &\subseteq \gamma^k \alpha \cdot (\gamma \cdot \mathcal{E}) \\ &\subseteq \gamma^{k+1} \mathcal{E}_\alpha \end{aligned} \quad (43) $$

As a result, by regrouping (40), (41), (42) and (43) one obtains

$$ h^{k+1}(\mathcal{E}_\alpha) \subseteq \gamma^{k+1} \mathcal{E}_\alpha. \quad (44) $$

So (38) holds for $k + 1$. Therefore, by recurrence, (38) is true for all $k \ge 1$. □

## V. AXIS-ALIGNED LYAPUNOV EQUATION

To verify the requirement of Theorem 3, one must find an ellipsoid $\mathcal{E}$ that verifies

$$ P_h(\mathcal{E}) \subset \mathcal{E}. \quad (45) $$

In practice, Theorem 6 can find an ellipsoid that will likely verify (45). With this theorem, one can solve a Lyapunov equation to find an ellipsoid $\mathcal{E}(\Gamma)$ such that

$$ J \cdot \mathcal{E}(\Gamma) \subseteq \mathcal{E}(\Gamma). \quad (46) $$

Knowing that $P_h(\mathcal{E}(\Gamma)) = (1 + \rho) \cdot J \cdot \mathcal{E}(\Gamma)$, if $\rho$ is small enough, then $\mathcal{E}(\Gamma)$ will likely verify $P_h(\mathcal{E}(\Gamma)) \subseteq \mathcal{E}(\Gamma)$. In practice, this inclusion is more likely to be verified when the shape matrices have the same eigenvectors, as the ellipdoids then have the same axis. Note that if $J$ is not Schur, the equation has no solution but the system is unstable anyway.

**Theorem 6.** (Axis-aligned discrete Lyapunov equation) Let $J \in \mathbb{R}^{n \times n}$ be a Schur matrix and let $P \in S_n^+$ be the unique solution to the discrete Lyapunov equation

$$ J^T P J - P = -Q, \quad (47) $$

with $Q = J^T J$. Let $\Gamma = P^{-\frac{1}{2}}$ and $\Gamma_J = (J \Gamma^2 J^T)^{\frac{1}{2}}$. Then $\Gamma$ and $\Gamma_J$ have the same eigenvectors and $\mathcal{E}(\Gamma_J) \subseteq \mathcal{E}(\Gamma)$.

*Proof.* From Theorem 1, $\mathcal{E}(\Gamma)$ is positive invariant with respect to the linear mapping $x_{k+1} = J \cdot x_k$. Thus $(J \cdot \mathcal{E}(\Gamma) \subseteq \mathcal{E}(\Gamma))$. Then, one has

$$ \Gamma_J^{-2} = J^{-T} \Gamma^{-2} J^{-1}. \quad (48) $$

Thus, from (47), one has

$$ \begin{aligned} J^T P J - P &= -J^T J, \\ P - J^{-T} P J^{-1} &= -I_n, \\ \Gamma^{-2} - J^{-T} \Gamma^{-2} J^{-1} &= -I_n, \\ \Gamma^{-2} - \Gamma_J^{-2} &= -I_n \end{aligned} \quad (49) $$

Therefore, $\Gamma^{-2}$ and $\Gamma_J^{-2}$ have the same eigenvector. So $\Gamma$ and $\Gamma_J$ also have the same eigenvector. □

## VI. IMPLEMENTATION

To solve the problem introduced in Section III, we propose the Algorithm 1. This algorithm is described by the following steps.

Firstly, the algorithm solves the discrete Lyapunov equation (47) to compute matrix $P \in S_n^+$. The solution can be approximated. To solve this equation, the recommended method depends on the dimension $n$ of the problem. With a small $n$, one can use matrix factorisation as in [5]. In high dimensions, one can use bilinear transformation as described in [4].

Secondly, the algorithm tests the inclusion in (18) with the ellipsoid $\mathcal{E}(10^{-\alpha} P^{-\frac{1}{2}})$ with scale factor $\alpha > 0$. If the inclusion is not verified, the ellipsoid is shrunk with higher values of $\alpha$, until the inclusion is verified.

In the algorithm, the most complex operations are matrix inversion and matrix square root which have polynomial complexity.
**Algorithm 1** Stability test for nonlinear discrete-time system

Inputs $J$, $[\frac{dh}{dx}]$

Outputs res, $\Gamma$

1: $P$= solve_discrete_Lyapunov($J$)

2: $\alpha = 1$

3: res = False

4: while res = False do

5: $\Gamma = 10^{-\alpha} \cdot P^{-\frac{1}{2}}$

6: $[\mathcal{E}]$= diag($\|\Gamma_1\|_2, \|\Gamma_2\|_2, \dots, \|\Gamma_n\|_2) \cdot [1_n]$

7: $[J] = [\frac{dh}{dx}] ([\mathcal{E}])$

8: $\rho = \text{ub} (\| (\Gamma^{-1}J^{-1} [J] \Gamma - I_n) \cdot [1_n] \|)$

9: $\Gamma_{out} = (1 + \rho) (J\Gamma^2J^T)^{\frac{1}{2}}$

10: $M = \Gamma_{out}^{-2} - \Gamma^{-2}$

11: res = is_definite Wand ($M$)

12: $\alpha = \alpha + 1$

13: end while

VII. APPLICATION

This Section presents an application of Algorithm 1 on the
nonlinear discrete-time system defined by

$$
\begin{align}
x_{k+1} &= h(x_k), \nonumber \\
&= A \cdot x_k + B \cdot \tanh(x), \tag{50}
\end{align}
$$

with $A \in \mathbb{R}^{n \times n}$ and $B \in \mathbb{R}^{n \times n}$. The evolution function $h$ is $C^1$ and the origin of $\mathbb{R}^n$ is an equilibrium point of the system. The Jacobian matrix of $h$ is

$$
\frac{dh}{dx}(\mathbf{x}) = \mathbf{A} + \mathbf{B} (\mathbf{I}_n - \tanh(\mathbf{x})^T \tanh(\mathbf{x})), \quad (51)
$$

and the Jacobian matrix at the origin is

$$
\mathbf{J} = \frac{dh}{dx}(0) = \mathbf{A} + \mathbf{B}. \qquad (52)
$$

Consider the inclusion function $[\frac{dh}{dx}]$ such that $\forall [\mathbf{x}] \in \mathbb{R}^n$,

$$
\left[ \frac{dh}{dx} \right] ([\mathbf{x}]) = \mathbf{A} + \mathbf{B} \left( \mathbf{I}_n - [\tanh] ([\mathbf{x}])^T [\tanh] ([\mathbf{x}]) \right) \quad (53)
$$

where $[\tanh]$ is an inclusion function of the hyperbolic tan-
gent.

The Algorithm 1 is applied with $J$, $[\frac{dh}{dx}]$ and ($\alpha_{max} = 10$) to compute a matrix $\Gamma \in S_n^+$ so that the system (50) is exponentially stable in $\mathcal{E} = \mathcal{E}(\Gamma)$. The algorithm is implemented in Python language with the libraries `Numpy` for matrix operations, `Codac` for Interval analysis and `Scipy` for solving Lyapunov equations. The code and further documentation is available at https://morgan-louedec.fr/index.php/discrete-system/

A. Results and Discussion

Figure 4 illustrates the case with $n = 2$ and

$$
\mathbf{A} = \begin{bmatrix} -0.8 & 2 \\ -1 & 1.6 \end{bmatrix}, \mathbf{B} = \begin{bmatrix} 0.1 & 0.01 \\ 0.01 & 0.1 \end{bmatrix}, \quad (54)
$$

such that **B** is small compared to **A**. The algorithm has
computed the matrix

$$
\Gamma \approx \begin{bmatrix} 0.17 & 0.08 \\ 0.08 & 0.10 \end{bmatrix}, \quad (55)
$$

that defines the ellipsoid $\mathcal{E} = \mathcal{E}(\Gamma)$ verifying $\mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}$.
So the system is exponentially stable in $\mathcal{E}$. Note that, the term $-\mathbf{J}^T\mathbf{J}$ in the Lyapunov Equation (47) gives both $\mathcal{P}_h(\mathcal{E})$ and $\mathcal{E}$ the same axis but with different lengths.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Ellipsoid $\mathcal{E}$ with $n = 2$ with the inclusion $\mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}$ and the box $[\mathcal{E}]$.</figcaption></figure>

Figure 5 illustrates a case with $n = 20$ where $\mathbf{A}$ and $\mathbf{B}$ were randomly chosen to make $\mathbf{J}$ a Schur matrix. This figure represents the projection of the ellipsoids $\mathcal{E}$ and $\mathcal{P}_h(\mathcal{E})$ on the two first dimensions of the system. Since these ellipsoids have a dimension 20, the inclusion $\mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}$ is not fully visible but the Algorithm has verified it. So the system is exponentially stable in $\mathcal{E}$. Note that Whereas $\mathcal{P}_h(\mathcal{E})$ and $\mathcal{E}$ share the same axis, their projection may have a different axis.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Projection onto the first two dimensions of the ellipsoid $\mathcal{E}$ with $n = 20$ and the box $[\mathcal{E}]$. the inclusion $\mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}$ is partially represented.</figcaption></figure>

Moreover, in practice, the size of $\mathcal{E}$ is affected by
• The dimension of the problem. The ellipsoid $\mathcal{E}$ is smaller in high dimension as illustrated by Figures 4 and 5. This is explained by the fact that $\rho$ is computed by adding the pessimism of every dimension.

• The non-linearity of the system. Whereas a linear system gives no pessimism $\rho = 0$, the more non-linear the system, the bigger the gain $\rho = 0$, the smaller the ellipsoid $\mathcal{E}$.

Finally, although the analytical expression of $h$ is given in this example, the algorithm only uses the Jacobian matrix $J$ and the inclusion function $[\frac{dh}{dx}]$. Thus the method can be used to study synchronised hybrid systems. These types of systems contain both continuous-time and discrete-time variables. They can be discretized into nonlinear discrete-time systems. Whereas the resulting evolution functions $h$ may have no analytical expressions it is possible to evaluate their Jacobian matrices. Thus one can use the Algorithm 1 on these discretized systems.

## VIII. CONCLUSION

In this paper, a guaranteed numerical method was proposed to prove the exponential stability of nonlinear discrete-time systems and to find a positive invariant ellipsoid. This method can be used when the analytical expression of the evolution function of the discrete system is unknown, as long as one can provide a function that encloses the values of the Jacobian matrix. This numerical method is also computationally tractable. Therefore, it can be applied to high-dimensional problems, where Lyapunov functions are difficult to find.

For future work, this method will be applied in the study of nonlinear synchronous hybrid systems. It will also be adapted to study the stability of systems with small perturbations.

## REFERENCES

[1] Franco Blanchini and Stefano Miani. Set-Theoretic Methods in Control. Systems & Control: Foundations & Applications. Springer International Publishing, Cham, 2015.

[2] Nicoletta Bof, Ruggero Carli, and Luca Schenato. Lyapunov Theory for Discrete Time Systems. Technical report, arXiv, 2018.

[3] Julien Calbert, Lucas N. Egidio, and Raphaël M. Jungers. An Efficient Method to Verify the Inclusion of Ellipsoids. IFAC-PapersOnLine, 56(2):1958–1963, January 2023.

[4] Zoran Gajic and Muhammad Tahir Javed Qureshi. Lyapunov matrix equation in system stability and control. Courier Corporation, 2008.

[5] James D. Hamilton. Time Series Analysis. Princeton University Press, princeton university press edition, 2020.

[6] Morris W. Hirsch, Stephen Smale, and Robert L. Devaney. Differential Equations, Dynamical Systems, and an Introduction to Chaos. Elsevier, 2013.

[7] Luc Jaulin, Michel Kieffer, Olivier Didrit, and Eric Walter. Applied Interval Analysis. Springer, London, 2001.

[8] Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud. Introduction to Interval Analysis. Society for Industrial and Applied Mathematics, January 2009.

[9] Andreas Rauh and Luc Jaulin. A computationally inexpensive algorithm for determining outer and inner enclosures of nonlinear mappings of ellipsoidal domains. International Journal of Applied Mathematics and Computer Science, 31(3):399–415, 2021.

[10] Sergey Staroletov. Automatic proving of stability of the cyber-physical systems in the sense of lyapunov with keymaera. In 2021 28th Conference of Open Innovations Association (FRUCT), pages 431–438. IEEE, 2021.

[11] Christophe Viel, Sylvain Bertrand, Kieffer Michel, and Piet-Lahanier Helene. New state estimators and communication protocol for distributed event-triggered consensus of linear multi-agent systems with bounded perturbations. IET Control Theory & Applications, 11(11):1736–1748, 2017.

[12] Shen Yan, Zhou Gu, and Ju H Park. Lyapunov-Function-Based Event-Triggered Control of Nonlinear Discrete-Time Cyber-Physical Systems. IEEE Transactions on circuits and systems, 69(6), 2022.

[13] Haijuan Zhao. Security control of cyber-physical switched systems under Round-Robin protocol: Input-to-state stability in probability. Information Sciences, 2020.

# APPENDIX

## Proof of Theorem 3

Consider an ellipsoid $\mathcal{E}$ such that

$$ \mathcal{P}_h(\mathcal{E}) \subset \mathcal{E}, \qquad (56) $$

and let $\Gamma = \mathcal{E}^{-1}(\mathcal{E})$. So, there exist $\gamma \in ]0, 1[$ such that

$$ \mathcal{P}_h(\mathcal{E}) \subset \gamma \cdot \mathcal{E}. \qquad (57) $$

Thus, from Theorem 5, one gets

$$ \mathbf{h}^k(\mathcal{E}_\alpha) \subset \gamma^k \cdot \mathcal{E}_\alpha. \qquad (58) $$

Then, let $\mathbf{x}_0 \in \mathcal{E}$ such that $\mathbf{x}_0 \neq 0$ and let $\alpha = \|\mathbf{x}_0\|_{\Gamma^{-2}} \le 1$. Consider the ellipsoid $\mathcal{E}_\alpha = \alpha \cdot \mathcal{E}$ and $\Gamma_\alpha = \mathcal{E}^{-1}(\mathcal{E}_\alpha) = \alpha \Gamma$. From (3), one has

$$ \begin{aligned} \|\mathbf{x}_0\|_{\Gamma_\alpha^{-2}} &= \sqrt{\mathbf{x}_0^T \Gamma_\alpha^{-2} \mathbf{x}_0}, \\ &= \frac{1}{\alpha} \sqrt{\mathbf{x}_0^T \Gamma_\alpha^{-2} \mathbf{x}_0}, \\ &= \frac{1}{\alpha} \|\mathbf{x}_0\|_{\Gamma^{-2}}, \\ &= \frac{1}{\alpha} \cdot \alpha = 1 \end{aligned} \qquad (59) $$

Therefore $\mathbf{x}_0 \in \mathcal{E}_\alpha$. Let $k \in \mathbb{N}$. From (1), one gets

$$ \mathbf{x}_k = \mathbf{h}^k(\mathbf{x}_0) \in \mathbf{h}^k(\mathcal{E}_\alpha). \qquad (60) $$

Then, From (58), one deduces

$$ \begin{aligned} \mathbf{x}_k &\in \gamma^k \cdot \mathcal{E}_\alpha, \\ &\in \gamma^k \cdot \mathcal{E}(\alpha \Gamma), \\ &\in \mathcal{E}(\gamma^k \alpha \Gamma). \end{aligned} \qquad (61) $$

As a consequence, with $\Gamma_k = \gamma^k \alpha \Gamma$, one has

$$ \begin{aligned} \|\mathbf{x}_k\|_{\Gamma_k^{-2}} &\le 1, \\ &\Leftrightarrow \sqrt{\mathbf{x}_k^T \Gamma_k^{-2} \mathbf{x}_k} \le 1, \\ &\Leftrightarrow \sqrt{\mathbf{x}_k^T \frac{1}{(\gamma^k \alpha)^2} \Gamma_k^{-2} \mathbf{x}_k} \le 1, \\ &\Leftrightarrow \frac{1}{\gamma^k \alpha} \sqrt{\mathbf{x}_k^T \Gamma_k^{-2} \mathbf{x}_k} \le 1, \\ &\Leftrightarrow \frac{1}{\gamma^k \alpha} \|\mathbf{x}_k\|_{\Gamma^{-2}} \le 1, \\ &\Leftrightarrow \|\mathbf{x}_k\|_{\Gamma^{-2}} \le \gamma^k \alpha, \\ &\Leftrightarrow \|\mathbf{x}_k\|_{\Gamma^{-2}} \le e^{-\beta k} \cdot \|\mathbf{x}_0\|_{\Gamma^{-2}}, \end{aligned} \qquad (62) $$

with $\beta = -\ln(\gamma) > 0$. This inequality is also verified if $\mathbf{x}_0 = 0$ because $\mathbf{x}_k = 0$ for all $k$. So (17) is true. Therefore, the system is exponentially stable in the ellipsoid $\mathcal{E}$.