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

HAL
open science

# Convex computation of regions of attraction from data using Sums-of-Squares programming

Oumayma Khattabi, Matteo Tacchi, Sorin Olaru

► To cite this version:

Oumayma Khattabi, Matteo Tacchi, Sorin Olaru. Convex computation of regions of attraction from data using Sums-of-Squares programming. IFAC Journal of Systems and Control, In press, pp.100361. 10.1016/j.ifacsc.2026.100361. hal-05379529v2

HAL Id: hal-05379529

https://hal.science/hal-05379529v2

Submitted on 11 Dec 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Convex computation of regions of attraction from
data using Sums-of-Squares programming

Oumayma Khattabi, Matteo Tacchi-Bénard, Sorin Olaru

December 11, 2025

## Abstract

The paper concentrates on the analysis of the finite horizon region of attraction for un-
known autonomous dynamical systems. The aim is to explore a data-driven approach based on
moment-sums-of-squares hierarchy, which enables novel region of attraction outer approxima-
tions despite the reduced information on the structure of the dynamics. The main contribution
of this work is bypassing the system model and, consequently, the recurring constraint on its
polynomial structure. Numerical experimentation showcases the influence of data on learned
approximating sets, offering a promising outlook on the potential of this method.

**Keywords** : Data-driven, SoS optimisation, Region of attraction approximation, Linear programming.

# 1 Introduction

Stability analysis of dynamical systems is one of the fundamental pillars of control theory. It is traditionally based on the system model, leveraging approaches such as Lyapunov's direct method [1], LaSalle's invariance principle [2], or the comparison principle in dissipativity theory [3]. While relying on a system model gives strong theoretical guarantees, defining a precise model is becoming harder with the growing complexity of modern systems. In recent years, data-driven methods increased in popularity due to the abundance of data, increased computational power, and the need to adapt to nonlinear real systems. As a result, many stability analysis approaches have been developed or adapted to prove stability directly from data.

Model-based and data-driven techniques have contributed significantly to the stability analysis of nonlinear systems. In both settings, several works stand out. For instance, [4] uses linear programming to learn piecewise affine Lyapunov functions for uncertain systems; while [5] uses second-order cone programming to approximate such functions from data. A related data-based approach is the continuous piecewise affine method proposed in [6].

Neural networks have also become practical tools in stability analysis. In [7], neural Lyapunov candidates are identified jointly with the system model from data to establish stability. They are further employed in [8] to generate control laws constrained by Lyapunov functions learned online, in [9] while mitigating the curse of dimensionality in Lyapunov function identification, and in [10] to compute constraint-admissible invariant sets.

Another approach is the use of polynomials: in [11], dissipativity is analysed from noisy data using piecewise Taylor approximations. On the model-based side, Sum of Squares (SoS) programming has been used to find control Lyapunov functions in [12], estimate the Region of Attraction (RoA) via polynomial Lyapunov functions in [13], and compute inner and outer Finite Horizon Region of Attraction (FH-RoA) approximations in [14, 15].

While SoS-based methods provide strong theoretical guarantees for stability analysis, they are typically limited to model-based polynomial systems. In contrast, data-driven approaches enable the
study of more general nonlinear dynamics. Some works use data to approach the unknown part of the model, and SoS methods to learn Lyapunov candidates [16] or controllers with guarantees [17]. This work, in contrast, examines the feasibility and effectiveness of a fully data-driven approach under minimal model information.

The remainder of the paper is organised as follows. In section 2, the theoretical framework of the data-driven concept is provided, as well as a trivial example in $\mathbb{R}$ to illustrate. Section 3 details the optimisation problem based on data and its complementary case of inner approximations. Numerical results are presented in section 4, along with an analysis of the computational difficulty and a proposed solution. An appendix is provided at the end of the paper to explain the moment-SoS hierarchy and prove the convergence of the optimisation problem.

**Notation:** Let $n \in \mathbb{N}$. $\|\cdot\|$ is the Euclidean norm. For $X \subset \mathbb{R}^n$, $\partial X$ is the boundary of $X$, $X^c$ its complement and vol($X$), its $n$-dimensional volume. We define $[n] \triangleq \{1, \dots, n\}$. $\mathbb{R}[x]^n$ is the set of $n$-dimensional vectors of real polynomials in $x$. $C^n(X)$ is the set of functions on $X$ with continuous derivatives up to order $n$. $W^{n,\infty}(X)$ is the set of functions with bounded derivatives up to order $n$ on $X$.

# 2 Data-driven region of attraction

## 2.1 Problem statement

We consider a dynamical system:

$$ \dot{x} = f(x) \quad (1) $$

defined by an unknown vector field $f: \mathbb{R}^n \to \mathbb{R}^n$. We seek the set $X_0(f)$ of initial states whose trajectories stay in $X \subset \mathbb{R}^n$ over $[0, T]$, and reach $X_T \subset X$ at time $T > 0$. We call $X_0(f)$ the FH-RoA, defined for a horizon $T$, an admissible set $X$, and a target set $X_T$. To underline its dependence on the original model (1), the set $X_0$ will be denoted as such:

$$ X_0(f) \triangleq \left\{ x_0 \in X \left| \begin{array}{l} \exists x(\cdot) \in C^1([0, T])^n \text{ s.t.} \\ x(0) = x_0, \quad x(T) \in X_T, \\ \forall t \in [0, T], \quad x(t) \in X \\ \text{and } \dot{x}(t) = f(x(t)) \end{array} \right. \right\}. \quad (2) $$

Although unknown, $f$ is accessible as follows:

**Assumption 1** The function $f$ is Lipschitz continuous, and the following information is available:

* an upper bound $M > 0$ on its Lipschitz constant,

* a finite and fixed sample of noise-free evaluation points: $D = \{(x_i, y_i = f(x_i))\}_{i \in [N]} \subset \mathbb{R}^{2n}$.

**Remark 1** Bounding arguments are common in data-driven analysis: partial derivative bounds estimated from data [11, 18], Lipschitz bounds [5], Hilbert norm bounds [19], etc.

**Remark 2** Assumption 1 implies that $\forall x \in X$, $f(x)$ lies within a state-dependent semi-algebraic uncertainty set:

$$ F_D(x) = \left\{ y \in \mathbb{R}^n \left| \begin{array}{l} \forall i \in [N], \\ \|y - y_i\|^2 \le M^2 \|x - x_i\|^2 \end{array} \right. \right\}. \quad (3) $$

**Remark 3** In this work, the data are assumed noise-free. Noisy datasets can still be incorporated by assuming $y_i = f(x_i) + z_i$, with $\|z_i\| \le \delta$ a bounded disturbance. This simply enlarges the uncertainty set by the disturbance bound $\delta$: $\|z_i\|^2 \le \delta^2 \Rightarrow \|y - y_i - z_i\|^2 \le M^2 \|x - x_i\|^2$.
**Remark 4** Even when $f$ is known, and $X$ and $X_T$ are simple sets, exact computation of $X_0(f)$ is rarely tractable. For instance, if $f \in \mathbb{R}[x]^n$, and $X$ and $X_T$ are semi-algebraic, the moment-SoS hierarchy enables the computation of certified inner ([15]) and outer ([14]) approximations of $X_0(f)$ via convex optimisation. The present contribution extends these works to the setting where $f$ can only be accessed through the dataset $D$, and the Lipschitz bound $M$.

To compute certified approximations of $X_0(f)$, we intend to bridge the gap between model-based moment-SoS frameworks in [15, 14], and the data-based approach in [5]. As a result, we make the following assumption on the structure of the admissible and target sets $X$ and $X_T$:

**Assumption 2** $X$ and $X_T$ are compact and there exist $n_X, n_T \in \mathbb{N}$, vectors $g_X \in \mathbb{R}[x]^{nx}$ and $g_T \in \mathbb{R}[x]^{n_T}$, s.t.:

$$X = \{x \in \mathbb{R}^n | g_X(x) \ge 0\} \qquad (4a)$$

$$X_T = \{x \in \mathbb{R}^n | g_T(x) \ge 0\} \qquad (4b)$$

where vector inequalities are considered component-wise.

## 2.2 Finite time region of attraction

We revisit the formulation proposed in [14] for outer FH-RoA approximation. The authors approximate the horizon $T$ FH-RoA of the target $X_T$ for a (model-based) polynomial control system with polynomial input constraints, seen as a differential inclusion. The resulting FH-RoA is defined identically to (2), with the additional requirement that there exists a control law generating the converging trajectory. In this contribution, we work with uncontrolled but unknown dynamics. The key observation is that [14] does not require the input $u$ to be a control input: the methodology applies if $u$ is replaced with an unknown disturbance $y$, and our framework corresponds to $\dot{x} = y$ with $y \in F_D(x)$ as described in (3). Hence, one can formulate the model uncertainty as a fixed implicit semi-algebraic set inclusion:

$$(x, \dot{x}) \in \Gamma_D \quad \text{with} \\ \Gamma_D \triangleq \{(x, y) \in \mathbb{R}^{2n} | g_X(x) \ge 0 \text{ and } y \in F_D(x)\}. \qquad (5)$$

This way, the data-based FH-RoA can be expressed as follows:

$$X_0 \triangleq \left\{ x_0 \in X \left| \begin{array}{l} \exists x(\cdot) \in W^{1,\infty}([0,T])^n \text{ s.t.} \\ x(0) = x_0, x(T) \in X_T \text{ &} \\ \forall t \in [0,T], (x(t), \dot{x}(t)) \in \Gamma_D \end{array} \right. \right\}. \qquad (6)$$

In words, $X_0$ is the set of all initial conditions of trajectories $x(\cdot)$, such that $x(T) \in X_T$ and at all times $t \in [0, T]$, $x(t) \in X$, and $\dot{x}(t) \in F_D(x(t))$.

**Remark 5** The description (6) of $X_0$ requires the existence of at least one dynamical system that maps $x(0) = x_0 \in X_0$ to $x(T) \in X_T$, subject to state constraints and compliant with the data. We refer to this case as the best-case FH-RoA, which provides an outer approximation of the FH-RoA of the real system. In contrast, the worst-case FH-RoA denotes the set $X_0^*$ corresponding to all dynamical systems compliant with the data that map $x(0) = x_0 \in X_0^*$ to $x(T) \in X_T$. This yields an inner approximation of the FH-RoA of the real system. We show that the worst-case inner approximation can be obtained from the best-case outer approximation of $X_0^*$, the complementary set of $X_0$.

**Remark 6** From a different perspective, the trajectories initiated in $X_0$ are viable ([20]) over the horizon $T$, with target set $X_T$, according to the uncertainty that can be inferred from the available data.
**Remark 7** It is well known that a control input and an unknown disturbance are mathematically equivalent; the only difference is that the control value can be chosen. The interplay between the existence of a control steering the system to the target and the requirement of reaching the target under all possible perturbations is standard in $H^\infty$ control. The originality of our approach lies in applying this observation to the notions of inner and outer approximations of the FH-RoA in a data-driven framework.

**Remark 8** Our approach differs from [14] in four aspects: (i) it does not require the underlying dynamics $f$ to be polynomial, they can be any Lipschitz continuous vector field; (ii) it requires no model at all, and is instead fully data-driven; (iii) it accounts for the more general coupled constraint $y \in F_D(x)$ instead of the decoupled $x \in X, u \in U$; (iv) it can provide either best-case outer approximations or worst-case inner approximations.

## 2.3 Illustration on a toy example

The concept is explained using the following 1D system defined for $x \in X = [-1, 1] \subset \mathbb{R}$ by:

$$ \dot{x} = \begin{cases} 2x(x^2 - 0.5^2) & \text{if } |x| \le 0.5 \\ x - 0.5 & \text{if } x \ge 0.5 \\ x + 0.5 & \text{if } x \le -0.5 \end{cases} \qquad (7) $$

For illustration purposes, the data is generated by the model $f$ given in (7) that is assumed unknown, and we show how the dynamics of (7) can be analysed with only data. Using the following dataset with three points: $D_1 = \{(-1; -0.5), (0; 0), (1; 0.5)\}$ and the Lipschitz constant $M = 1$ of the system, we plot the area $F_{D_1}(X)$ where all the possible function values exist (see green region in figure 1). The largest possible FH-RoA corresponds to the piecewise affine function constructed out of the Lipschitz inequality limits with values in quadrants II and IV. We call this function $f_{best}$, associated with the best-case RoA.

<figure><img src="image_3.png" /><figcaption>Figure 1: Toy example with 3 data points</figcaption></figure>

Given the same dataset $D_1$ and $X_T = [-0.25, 0.25]$ for $T = 1s$, one can analytically compute $f_{best}$ as in figure 2a, as well as its corresponding FH-RoA. In this case, the FH-RoA of system (7) for
<figure><img src="image_4.png" /><figcaption>Figure 2: Identification of the FH-RoA with 3 data points</figcaption></figure>

$T$ and $X_T$ is $X_0(f) = [-0.34, 0.34]$, and the FH-RoA of $f_{best}$ for the same parameters is $X_0(f_{best}) = [-0.408, 0.408]$ (see figure 2b), making it an outer approximation of $X_0(f)$.

Adding 2 points $D_2 = D_1 \cup \{(-0.3; 0.096), (0.3; -0.096)\}$ imposes tighter constraints (see figure 3a), reducing the feasible solution space (see figure 3b). However, such an analytical approach is cumbersome even for scalar dynamics.

<figure><img src="image_5.png" /><figcaption>Figure 3: Identification of the FH-RoA with 5 data points</figcaption></figure>

**Remark 9** The description of $X_0$ in (6) implicitly allows non-Lipschitz behaviour in the dynamics $\dot{x}$ as a function of $x$, requiring the Lipschitz condition to hold only on the data points (formally, $f_{best}$ in Fig. 1 is allowed to jump within the green envelope). However, the toy example suggests that the $f_{best}$, whose FH-RoA corresponds to the best-case FH-RoA, always lies on the boundary of the uncertainty set $F_D(x)$, and is exactly M-Lipschitz continuous.
# 3 The optimisation framework

## 3.1 Linear program on functions

We propose a generalisation of [14] from explicit semi-algebraic differential inclusions (and polynomial control systems) to implicit semi-algebraic set inclusions (and systems with bounded uncertainty).

Let $\mathcal{L}_y$ denote the operator acting on $v \in C^1([0, T] \times X)$, defined for $(t, x, y) \in [0, T] \times \Gamma_D$ by:

$$ \mathcal{L}_y v(t, x, y) \triangleq \frac{\partial v}{\partial t}(t, x) + y^\top \frac{\partial v}{\partial x}(t, x) \quad (8) $$

The Linear Program (LP) reads as follows:

$$ \underset{v \in C^1([0, T] \times X) \atop w \in C^0(X)}{\text{minimise}} \quad \int_X w(x) \, dx \qquad (9a) $$

$$ \text{s.t.} \quad \mathcal{L}_y v(t, x, y) \leq 0 \qquad \forall (t, x, y) \in [0, T] \times \Gamma_D \qquad (9b) $$

$$ w(x) \geq 0 \qquad \forall x \in X \qquad (9c) $$

$$ w(x) \geq v(0, x) + 1 \qquad \forall x \in X \qquad (9d) $$

$$ v(T, x) \geq 0 \qquad \forall x \in X_T \qquad (9e) $$

**Remark 10** The difference with [14] is that (9b) is defined on $[0, T] \times \Gamma_D$ instead of $[0, T] \times X \times U$ and controlled dynamics $\dot{x} = f(t, x, u)$ are replaced with unknown speed of variation $\dot{x} = y \in F_D(x)$.

**Lemma 1** Let $(v, w)$ be feasible for problem (9) and define

$$ \hat{X}_0(w) \triangleq \{x \in X \mid w(x) \geq 1\}. $$

Recalling the definition (6) of $X_0$, it holds $X_0 \subset \hat{X}_0(w)$. Moreover, considering a minimising sequence $(v_\epsilon, w_\epsilon)$ s.t.

$$ 0 \le \int w_{\epsilon}(x) \, dx - w^{*} \le \epsilon, $$
where $w^*$ is the optimal value of problem (9), it holds
$$ \mathrm{vol}(\hat{X}_0(w_{\epsilon}) \setminus X_0) \xrightarrow{\epsilon \to 0} 0. \qquad (10) $$

**Proof.** Almost identical to the ones in [14, 15, 21, 22]: Let $x_0 \in X_0$ and consider an admissible trajectory $x(t)$ such that $x(0) = x_0$, $x(T) \in X_T$ and at all times $(x(t), \dot{x}(t)) \in \Gamma_D$. Such a trajectory exists by definition of $X_0 \ni x_0$. Now, consider the map $t \mapsto v(t, x(t))$; by design (8) of the operator $\mathcal{L}_y$, constraint (9b) ensures that

$$ 0 \ge \mathcal{L}_y v(t, x(t), \dot{x}(t)) = \frac{d}{dt} \left[ v(t, x(t)) \right] \quad (*) $$

$$ \text{so that} \quad w(x_0) \stackrel{(9d)}{\ge} 1 + v(0, x_0) \stackrel{(*)}{\ge} 1 + v(T, x(T)) \stackrel{(9e)}{\ge} 1, $$

i.e. $x_0 \in \hat{X}_0(w)$. Eventually, (10) follows from

$$ \mathrm{vol}(\hat{X}_0(w) \setminus X_0) = \mathrm{vol}(\hat{X}_0(w)) - \mathrm{vol}(X_0) \le \int w(x) dx - w^* \quad \square $$

In practice, these converging approximations are computed via the SoS hierarchy: certificates $(v, w)$ are parame-trised as $v(t, x) = c_v^\top \varphi(t, x)$ and $w(x) = c_w^\top \psi(x)$, where $\varphi(t, x)$ and $\psi(x)$ are polynomial vectors, and $c_v, c_w$ are real vectors; inequalities are parametrised by Linear Matrix Inequality (LMI), for instance constraint (9c) reads $w(x) = \psi(x)^\top Q_1 \psi(x) + g_X(x) \psi(x)^\top Q_g \psi(x)$, $Q_1, Q_g \ge 0$. See Appendix A for details on the SoS hierarchy.
## 3.2 Toy example revisited

The numerical example introduced in (7) is reused with dataset $D_1$. In figure 4, we compare the outer approximation $\hat{X}_0(w)$ (in blue, obtained from LP (9)) with $X_0(f_{best,D_1})$ (in red). This provides a good approximation of the real FH-RoA, as $w(x)$ closely mimics its indicator function and reaches $w(x) = 1$ along the set's boundary.

<figure><img src="image_6.png" /><figcaption>Figure 4: LP based FH-RoA approximation with 3 data points.</figcaption></figure>

Applying the same optimisation procedure with dataset $D_2$ produces a narrower outer approximation $\hat{X}_0(w)$ of the FH-RoA $X_0(f)$ by approaching the new $X_0(f_{best,D_2})$ (indicator function in red), as illustrated in figure 5. This highlights how the available data affects the accuracy of the approximation, a point further investigated in section 4.

## 3.3 Worst-case inner approximations

As highlighted in Remark 5, the above framework concerns the best-case RoA outer approximation. The worst-case FH-RoA inner approximation is relevant for data-driven analysis, i.e., requiring that all dynamical systems consistent with the data steer $x_0$ to $X_T$ under state constraints. Below, we provide a procedure to transition from the former to the latter, in the spirit of [15]; though that reference considers closed-loop systems, whereas we retain the uncertainty set inclusion $\dot{x} \in F_D(x)$.

Consider the FH-RoA $\tilde{X}_0$ obtained by replacing $X_T$ with $X_T^c$ in (6), as well as, for any $\tau \in (0, T)$, the FH-RoA $\tilde{X}_\tau^\partial$ obtained by replacing $T$ with $\tau$ and $X_T$ with $\partial X$ in (6). Then, one can straightforwardly compute outer approximations of $\tilde{X}_0$ and $\tilde{X}_\tau^\partial$ doing the same replacements in (9). Moreover, defining the worst-case FH-RoA:

$$X_0^* \triangleq \left\{ x_0 \in X \left| \begin{array}{l} \forall x(\cdot) \in W^{1,\infty}(0,T)^n \text{ s.t.} \\ \forall t \in [0,T], \dot{x}(t) \in F_D(x(t)) \\ \text{and } x(0) = x_0, \text{ it holds} \\ \forall t \in [0,T], x(t) \in X \\ \text{and } x(T) \in X_T \end{array} \right. \right\}, \quad (11a)$$

one gets the following complement formula:
<figure><img src="image_7.png" /><figcaption>Figure 5: LP based FH-RoA approximation with 5 data points</figcaption></figure>

$$ (X_0^*)^c = \left\{ \begin{array}{l} x_0 \in X \quad \left| \begin{array}{l} \exists x(\cdot) \in W^{1,\infty}(0,T)^n \text{ s.t.} \\ \forall t \in [0,T], \quad \dot{x}(t) \in F_D(x(t)), \\ x(0) = x_0 \quad \text{and} \\ \left[ \begin{array}{l} \exists \tau \in (0,T); x(\tau) \in \partial X \\ \text{or} \quad x(T) \in X_T^c \end{array} \right] \end{array} \right\} \\ \\ = \bigcup_{0<\tau<T} \tilde{X}_\tau^\partial \cup \tilde{X}_0, \end{array} \right. \quad (11b) $$

so that outer approximations of $(X_0^*)^c$ can be deduced from a minimising sequence for the following problem:

$$ \begin{array}{ll} \text{minimise} & \displaystyle\int_{v \in C^1([0,T] \times X)} w(x) \, dx \qquad (12a) \\ & w \in C^0(X) \end{array} $$

$$ \text{s.t.} \quad \mathcal{L}_y v(t, x, y) \le 0 \qquad \forall (t, x, y) \in [0, T] \times \Gamma_D \quad (12b) $$

$$ w(x) \ge 0 \qquad \forall x \in X \qquad (12c) $$

$$ w(x) \ge v(0, x) + 1 \qquad \forall x \in X \qquad (12d) $$

$$ v(T, x) \ge 0 \qquad \forall x \in X \setminus X_T \qquad (12e) $$

$$ v(t, x) \ge 0 \qquad \forall (t, x) \in [0, T] \times \partial X \qquad (12f) $$

Since the complements of outer approximations of $(X_0^*)^c$ are inner approximations of $X_0^*$, and similarly to Lemma 1, any $(v, w)$ feasible for (12) is such that:

$$ X_0^* \supset \{x \in X \mid w(x) \le 1\} \triangleq \tilde{X}_0^*, \quad (13) $$

then the moment-SoS hierarchy can be implemented to compute such inner approximations $\tilde{X}_0^*$.
**Theorem 1** For a given dataset $D = \{(x_i, f(x_i))\} \subset \mathbb{R}^{2n}$, let $X_0^*(D)$ be the corresponding worst-case data-based FH-RoA as per (11a). Then, for all $N \in \mathbb{N}$ there exists a dataset $D_N = \{(x_i, f(x_i))\}_{i=1}^N \subset \mathbb{R}^{2n}$ such that:

$$ \forall N \ge 1, \quad X_0^*(D_N) \subset X_0(f) \qquad (14a) $$

$$ \text{and} \quad \text{vol}(X_0(f) \setminus X_0^*(D_N)) \underset{N \to \infty}{\longrightarrow} 0. \qquad (14b) $$

**Proof.** See appendix B. □

**Remark 11** The difference between this framework and the method proposed in [5] is twofold: here, the optimal solutions are proved to converge towards the actual FH-RoA, while [5] only guarantees convergence of certificates towards a Lyapunov function; plus, this framework uses polynomial certificates and SoS programming while the latter is based on piecewise affine Lyapunov functions and second order cone programming. While the uncertainty modeling is close, the methods and the guarantees they yield are completely different.

# 4 Numerical results

All the simulations were conducted on a laptop with 32GB of memory and an AMD Ryzen 7 7735U processor. The algorithm was coded using YALMIP [23], with MOSEK [24] being the selected solver.

## 4.1 Analysis of data position's influence

In this part, the influence of the position of data points on the outer approximation (LP (9)) of the FH-RoA is explored. For this, we work with the same 1D system defined in (7) on $X = [-1, 1]$, and we keep $X_T = [-0.25, 0.25]$ for $T = 1$s. We consider the dataset $D(p) = D_1 \cup \{(-p; f(-p)), (p; f(p))\}$ of 5 data points, of which three are fixed (in $D_1$) and two vary in position according to $p \in [0.1, 0.5]$. In Figure 6, and in light of the uncertainty set representation of Figure 1, the more the two data points of positions $p$ and $-p$ restrain $F_D(X)$, the closer the outer approximation is to the FH-RoA of the real function.

<figure><img src="image_8.png" /><figcaption>Figure 6: Simulation using the variable dataset $D(p)$</figcaption></figure>
## 4.2 Numerical example in $\mathbb{R}^2$

The SoS reformulations of the LP (9) is solved for a dynamical system with a Lipschitz constant $M = 1$ defined on $x \in X = [-0.8, 0.8]^2 \subset \mathbb{R}^2$ as follows:

$$ \dot{x} = \begin{cases} 2x(\|x\|^2 - 0.5^2) & \text{if } \|x\| \le 0.5 \\ x\left(1 - \frac{1}{2\|x\|}\right) & \text{else.} \end{cases} \qquad (15) $$

We seek to find the FH-RoA of the system with $x \in X_T = \{x \in X \mid (0.25^2 - \|x\|^2) \ge 0\}$ at time $T = 1$s. For this, we fix the degree of the polynomials to 10, and the dataset is generated randomly in $X$ with $N = 50$. See the results in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7: Simulation results for $X_T$</figcaption></figure>

In 7c, the obtained approximation $\hat{X}_0(w)$ in blue approaches the actual FH-RoA of the system in red. The number of data points and the chosen degree of the polynomials influence the algorithm's performance, the approximation's accuracy, and can lead to a better approximation of the true function FH-RoA. The dataset shown in Figure 7a influences the shape of the obtained $\hat{X}_0(w)$ through its effect on the derived best-case function. By providing the LP with a different dataset, we can generate a new $\hat{X}_0(w)$ that addresses the shape disfigurement of the previous result without increasing the number of data points, as illustrated in figure 8.
<figure><img src="image_10.png" /><figcaption>(a) new generated dataset</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) polynomial $w(x)$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) resulting FH-RoA</figcaption></figure>

<caption>Figure 8: Alternative simulation results for $X_T$</caption>
To test the validity of the algorithm for non-symmetric FH-RoAs, the previous $X_T$ is changed into a non-symmetric set for the same time $T = 1s$:

$$X'_{T} = \{ x \in X \mid \begin{array}{l} (0.25^2 - x^2) \geq 0 \\ x_{(1)} + x_{(2)} \geq 0 \end{array} \}$$

The simulation results from such a change are shown in Figure 9. The shape of the approximating set $\hat{X}_0$ in Figure 9c changes to accommodate the new target set $X'_T$ so that it still approaches the real FH-RoA.

<figure><img src="image_13.png" /><figcaption>Figure 9: Simulation results for $X'_T$</figcaption></figure>

## 4.3 Computational challenges and proposed resolution

This proposed approach offers several advantages, including that $f$ does not need to be polynomial or known like it is supposed in [14, 15]. It is also proven to converge, as shown in Appendix B. It faces, however, some numerical problems when the value of $M$ is high or when $T \gg 1s$. In the first case, a high overestimate of $M$ is not informative enough for the LP to approach the FH-RoA (see the influence of $M$ on the convergence proof in Appendix B). When $T \gg 1s$, even model-based approaches are hampered by numerical conditioning issues in the SDP approximations: this stems
from the geometric properties of Putinar's Positivstellensatz as detailed in [25]; this is worse in the data-based setting where model approximation errors accumulate along the time horizon.

In future work, we aim to address this issue by introducing a partitioning strategy for the set X, inspired by the approach in [26, 5]. This method is expected to enhance the results by introducing local Lipschitz constants and data constraints for each cell in the set. It would also allow for a lower degree of polynomials with no significant drawbacks to accuracy.

Another potential improvement to the results can be achieved by changing the polynomial basis of the SoS formulation. An SoS decomposition expresses a polynomial as a quadratic form in a chosen polynomial basis with a semi-definite positive weight matrix, which can be checked using LMI. The polynomial basis influences the accuracy of the numerical approximation as shown in [27]. This idea has been used in [28] to drastically improve the scalability of SoS programming, which is known to be its main practical weakness.

# 5 Conclusion

This paper introduces a novel data-driven approach to approximate the FH-RoA via the moment-SoS hierarchy. This methodology effectively provides outer approximations of the real FH-RoA using a dataset and a Lipschitz constant of the unknown system dynamics. In contrast to model-based approaches, using data circumvents the limitation of the system's polynomiality. While the results are promising, further improvements are needed to address the scalability issue arising from numerical challenges. Future work includes introducing set partitioning to address the scalability issue, further developing the worst-case inner approximation part, investigating the influence of data on conservatism, and testing a polynomial basis for the SoS formulation other than monomials.

# References

[1] A. M. Lyapunov, “The general problem of the stability of motion,” *Int. J. Control*, vol. 55, no. 3, pp. 531–534, 1992.

[2] J. La Salle and S. Lefschetz, *Stability by Liapunov's direct method with applications by Joseph L Salle and Solomon Lefschetz*, vol. 4. Elsevier, 2012.

[3] Bitsoris, *Comparison methods in control*. Springer, 2025.

[4] P. Julian, J. Guivant, and A. Desages, “A parametrization of piecewise linear Lyapunov functions via linear programming,” *Int. J. Control*, vol. 72, pp. 702–715, Jan 1999.

[5] M. Tacchi, Y. Lian, and C. Jones, “Robustly learning regions of attraction from fixed data,” *IEEE Transactions on Automatic Control*, 2025.

[6] S. F. Hafstein, C. M. Kellett, and H. Li, “Computing continuous and piecewise affine Lyapunov functions for nonlinear systems,” *J. Comput. Dyn.*, vol. 2, no. 2, pp. 227–246, 2016.

[7] J. Z. Kolter and G. Manek, “Learning stable deep dynamics models,” *NeurIPS*, vol. 32, 2019.

[8] Y. Min, S. M. Richards, and N. Azizan, “Data-driven control with inherent Lyapunov stability,” in *IEEE 62nd CDC*, pp. 6032–6037, 2023.

[9] L. Grüne, “Overcoming the curse of dimensionality for approximating Lyapunov functions with deep neural networks under a small-gain condition,” *IFAC-PapersOnLine*, vol. 54, no. 9, pp. 317–322, 2021.
[10] D. Kim and H. J. Kim, "Estimation of constraint admissible invariant set with neural Lyapunov function," in *IEEE 63rd CDC*, pp. 5032-5039, 2024.

[11] T. Martin and F. Allgöwer, "Data-driven system analysis of nonlinear systems using polynomial approximation," *IEEE Trans. Autom. Control*, vol. 69, pp. 4261-4274, Jul 2024.

[12] W. Tan and A. Packard, "Searching for control Lyapunov functions using SoS programming," *sibi*, vol. 1, no. 1, 2004.

[13] B. Biswas, D. Ignatyev, A. Zolotas, and A. Tsourdos, "Region of attraction estimation using union theorem in sum-of-squares optimization," *arXiv:2305.11655*, 2023.

[14] D. Henrion and M. Korda, "Convex computation of the region of attraction of polynomial control systems," *IEEE Trans. Autom. Control*, 2014.

[15] M. Korda, D. Henrion, and C. Jones, "Inner approximations of the region of attraction for polynomial dynamical systems," in *IFAC NOLCOS*, 2013.

[16] B. K. Colbert and M. M. Peet, "Using trajectory measurements to estimate the region of attraction of nonlinear systems," in *IEEE 57th CDC*, pp. 2341-2347, 2018.

[17] D. Han and H. Huang, "Sum-of-squares program and safe learning on maximizing the region of attraction of partially unknown systems," in *IEEE ASCC*, pp. 875-881, 2022.

[18] A. Makdesi, A. Girard, and L. Fribourg, "Efficient data-driven abstraction of monotone systems with disturbances," *IFAC-PapersOnLine*, vol. 54, no. 5, pp. 49-54, 2021.

[19] E. T. Maddalena, P. Scharnhorst, and C. N. Jones, "Deterministic error bounds for kernel-based learning techniques under bounded noise," *Automatica*, vol. 134, 2021.

[20] J.-P. Aubin, A. M. Bayen, and P. Saint-Pierre, *Viability theory: new directions*. Springer Science, 2011.

[21] M. Korda, D. Henrion, and C. Jones, "Convex computation of the maximum controlled invariant set for polynomial control systems," *SIAM J. Control Optim.*, 2014.

[22] A. Oustry, M. Tacchi, and D. Henrion, "Inner approximations of the maximal positively invariant set for polynomial dynamical systems," *IEEE L-CSS*, 2019.

[23] J. Löfberg, "Yalmip: A toolbox for modeling and optimization in matlab," in *CACSD Conf.*, 2004.

[24] MOSEK ApS, *The MOSEK optimization toolbox for MATLAB manual. Version 10.1.*, 2024.

[25] C. Schlosser, M. Tacchi-Bénard, and A. Lazarev, "Convergence rates for the moment-sos hierarchy," *Numerical Algebra, Control and Optimization*, vol. 16, pp. 105-156, 2026.

[26] V. Cibulka, M. Korda, and T. Hanis, "Spatio-temporal decomposition of sum-of-squares programs for the region of attraction and reachability," *IEEE L-CSS*, 2021.

[27] D. Henrion, J. B. Lasserre, and C. Savorgnan, "Approximate volume and integration for basic semialgebraic sets," *SIAM review*, vol. 51, no. 4, pp. 722-743, 2009.

[28] A. A. Ahmadi and G. Hall, "Sum of Squares basis pursuit with linear and second order cone programming," in *Algebraic and Geometric Methods in Discrete Mathematics*, American Mathematical Society, 2017.
[29] M. Putinar, “Positive polynomials on compact semi-algebraic sets,” Indiana University Mathematics Journal, 1993.

[30] J. B. Lasserre, Moments, positive polynomials and their applications. Imperial College Press, 2010.

# A The sums-of-squares hierarchy

Introducing the cone

$$ \Sigma[z] = \{p_1(z)^2 + \dots + p_K(z)^2 \mid K \in \mathbb{N}, p_1, \dots, p_K \in \mathbb{R}[z]\} $$

of polynomial SoS in the variable $z \in \{x, (t, x, y)\}$, the univariate polynomial $q(t) = (T - t)t$, the size $N$ vector

$$ \gamma(x, y) = (M^2 \|x - x_i\|^2 - \|y - y_i\|^2)_{i=1}^N $$

and, for any vector of polynomials $g(z) \in \mathbb{R}[z]^p$ in the variable $z$, the quadratic module defined as

$$ \mathcal{Q}(g(z)) = \left\{ \sigma(z) + s(z)^{\top} g(z) \mid \begin{array}{l} \sigma(z) \in \Sigma[z], \\ s(z) \in \Sigma[z]^p \end{array} \right\}, \qquad (16) $$

the LP (9) can be recast into the following:

$$ \underset{\substack{v \in \mathbb{R}[t, x] \\ w \in \mathbb{R}[x]}}{\text{minimise}} \int w(x) \, dx \qquad (17a) $$

$$ \text{s.t.} - \mathcal{L}_y v(t, x, y) \in \mathcal{Q} \left( q(t), g_X(x), \gamma(x, y) \right) \qquad (17b) $$

$$ w(x) \in \mathcal{Q}(g_X(x)) \qquad (17c) $$

$$ w(x) - v(0, x) - 1 \in \mathcal{Q}(g_X(x)) \qquad (17d) $$

$$ v(T, x) \in \mathcal{Q}(g_T(x)) \qquad (17e) $$

**Proposition 1** *Problem (17) is equivalent to problem (9) in the sense that any $(v, w)$ feasible for (17) is feasible for (9) and minimising sequences $(v_\epsilon, w_\epsilon)_\epsilon$ for (9) can be approximated by minimising sequences $(\hat{v}_\epsilon, \hat{w}_\epsilon)_\epsilon$ for (17).*

**Proof.** Since by design any $\sigma(z) \in \Sigma[z]$ is non-negative, it is clear that any $q(z) \in \mathcal{Q}(g(z))$ is non-negative on $\{z \mid g(z) \ge 0\}$. Hence, for $x \in \{\text{b, c, d, e}\}$, (17x) implies (9x) and feasibility for (17) implies feasibility for (9).

Next, let $(v_\epsilon, w_\epsilon)$ be a minimising sequence for (9) that is strictly feasible (meaning that all the inequality constraints are strict). Since by Assumption 2 X is compact, the Weierstrass theorem yields a polynomial approximation $(\hat{v}_\epsilon, \hat{w}_\epsilon)$ of $(v_\epsilon, w_\epsilon)$ that is also a strictly feasible minimising sequence for (9). Eventually, we state Putinar's Positivstellensatz [29]: under mild assumptions (that are implied by Assumption 2 in our case), if $q(z) > 0$ on $\{z \mid g(z) \ge 0\}$, then $q(z) \in \mathcal{Q}(g(z))$. Then, since the $(\hat{v}_\epsilon, \hat{w}_\epsilon)$ are strictly feasible for (9), they are also feasible for (17). Since (17) is a strengthening of (9), any minimising sequence of (9) that is also feasible for (17) is a minimising sequence of (17). $\square$

The quadratic module $\mathcal{Q}(g(z))$ can be truncated into a bounded degree quadratic module $\mathcal{Q}_d(g(z))$ where $\sigma(z)$ and each term $s_j(z) \cdot g_j(z)$, $j \in [p]$ of $s(z)^\top g(z)$ have degree at most $d \in \mathbb{N}$. Such a truncated quadratic module has the particularity of being LMI-representable (see e.g. Prop. 2.1 in [30]). Then, the SoS hierarchy consists of the sequence of problems obtained by replacing $\mathcal{Q}$ with
$Q_d, d \in \mathbb{N}$ in (17); since the $Q_d$ are nested (in the sense that $Q_d(g(z)) \subset Q_{d+1}(g(z)) \subset Q(g(z)))$, this produces a monotonic sequence of convex, finite-dimensional strengthenings of (17).

Problem (17) has a linear cost and any of its finite-dimensional truncations has a non-empty, compact feasible set, so they all have an optimal solution. Moreover, since $Q(g(z)) = \cup_{d \in \mathbb{N}} Q_d(g(z))$, these optimal solutions form a minimising sequence for (17) and hence for (9), yielding a converging sequence of outer approximations of $X_0$ as described in (10) when the degree bound $d$ goes to infinity.

The main drawback of approximating (9) with truncated instances of (17) is the computational burden: at degree $d$, the involved LMI have size growing like $(2n+d+1)d$, i.e. combinatorially in the state dimension as well as in the degree of the involved polynomials. Furthermore, the LMI representation of $\Sigma[z]$ heavily depends on the choice of the polynomial basis for $\mathbb{R}[z]$, which can have significant effects on the numerical behaviour of the resulting Semi Definite Programming (SDP) problems.

## B Convergence of data-based regions of attraction

The proof of Theorem 1 is as follows:

**Proof.** (14a) holds by definition of the worst-case RoA $X_0^*(D)$: any of its elements only initialises trajectories that reach the target under dynamics compatible with the dataset $D$. By design, $f$ is compatible with $D$, so that its trajectories initialized in $X_0^*(D)$ indeed reach the target, i.e. $X_0^*(D) \subset X_0(f)$. Then, we prove (14b). Both $X_0(f)$ and $X_0^*(D)$ are respectively approximated by inside by:

$$ \tilde{X}_0(w_0) = \{x \in X \mid w_0(x) \le 1\}, $$

and:

$$ \tilde{X}_0(w_D) = \{x \in X \mid w_D(x) \le 1\}, $$

with $w_D$ feasible for the data-based LP in (12), and $w_0$ feasible for the model-based LP as described in [14]:

$$
\begin{array}{ll}
\text{minimise} & \displaystyle\int_X w_0(x) \, dx \\
\text{s.t.} & \mathcal{L}_f v_0(t, x) \le 0 \quad \forall (t, x) \in [0, T] \times X \\
& w_0(x) \ge 0 \quad \forall x \in X \\
& w_0(x) \ge v_0(0, x) + 1 \quad \forall x \in X \\
& v_0(T, x) \ge 0 \quad \forall x \in X \setminus X_T \\
& v_0(t, x) \ge 0 \quad \forall (t, x) \in [0, T] \times \partial X
\end{array}
\qquad (18)
$$

where $\mathcal{L}_f$ denotes the operator:

$$ v(t, x) \mapsto \partial_t v(t, x) + f(x)^{\top} \partial_x v(t, x). $$

A feasibility-preserving direction for (18) is defined by $v_\delta(t, x) = T - t$ and $w_\delta(x) = T \ge 0$, with:

$$
\begin{array}{ll}
w_\delta(x) - v_\delta(0, x) = 0 & \forall x \in X \\
v_\delta(T, x) = 0 & \forall x \in X \setminus X_T \\
v_\delta(t, x) \ge T - t = 0 & \forall (t, x) \in [0, T] \times \partial X \\
\mathcal{L}_f v_\delta(t, x) = -1 & \forall (t, x) \in [0, T] \times X
\end{array}
$$

Moreover, the following holds:

$$ \int_X w_\delta(x) \, dx = T \operatorname{vol}(X) \triangleq C. $$
From [15], we know that, for $\varepsilon > 0$, there exists an $\varepsilon/2$-optimal feasible solution $(v_\varepsilon, w_\varepsilon)$ for (18) in the sense that:

$$ \tilde{X}_0(w_\varepsilon) \subset X_0(f), \quad \text{and} $$

$$ \mathrm{vol}(X_0(f) \setminus \tilde{X}_0(w_\varepsilon)) \le \int_X w_\varepsilon(x)dx - \mathrm{vol}(X \setminus X_0(f)) \le \frac{\varepsilon}{2}. $$

Let $\eta = \varepsilon/2C > 0$, and define:

$$ v_\eta(t, x) = v_\varepsilon(t, x) + \eta \cdot v_\delta(t, x), $$

$$ w_\eta(x) = w_\varepsilon(x) + \eta \cdot w_\delta(x) \ge T \cdot \eta > 0. $$

We check that, on the appropriate domains, it holds:

$$ \begin{aligned} w_\eta(x) - v_\eta(0, x) &= \underbrace{w_\varepsilon(x) - v_\varepsilon(0, x)}_{\ge 1} \\ &\quad + \underbrace{(w_\delta(x) - v_\delta(0, x)) \cdot \eta}_{=0} \ge 1 \\ v_\eta(T, x) &= \underbrace{v_\varepsilon(T, x)}_{\ge 0} + \underbrace{\eta \cdot v_\delta(T, x)}_{=0} \ge 0 \\ v_\eta(t, x) &= \underbrace{v_\varepsilon(t, x)}_{\ge 0} + \underbrace{\eta \cdot v_\delta(t, x)}_{\ge 0} \ge 0 \end{aligned} $$

$$ \mathcal{L}_f v_\eta(t, x) = \underbrace{\mathcal{L}_f v_\varepsilon(t, x)}_{\le 0} + \underbrace{\mathcal{L}_f v_\delta(t, x)) \cdot \eta}_{=-\eta} \le -\eta, $$

so that $(v_\eta, w_\eta)$ is feasible for (18). We are going to derive a condition on the dataset $D$ for $(v_\eta, w_\eta)$ to be feasible for the data-based problem (12). Notice that by construction, constraints (12c)–(12f) are satisfied by $(v_\eta, w_\eta)$. Hence we look at constraint (12b). Let $t \in [0, T]$, $(x, y) \in \Gamma_D$. We just need to prove:

$$ y^\top \partial_x v_\eta(t, x) \le f(x)^\top \partial_x v_\eta(t, x) + \eta, $$

since we already know that the RHS is not greater than $-\partial_t v_\eta(t, x)$. Hence, we look for a condition such that:

$$ (y - f(x))^\top \partial_x v_\eta(t, x) \le \eta. $$

Using the Cauchy-Schwarz inequality, we get:

$$ (y - f(x))^\top \partial_x v_\eta(t, x) \le \|y - f(x)\| \cdot \|\partial_x v_\eta(t, x)\|. $$

Moreover, for any $i \in [N]$, it holds:

$$ \|y - f(x)\| \le \underbrace{\|y - f(x_i)\|}_{\le M \|x - x_i\|} + \underbrace{\|f(x_i) - f(x)\|}_{\le M \|x - x_i\|} \le 2M \|x - x_i\| $$

where the first braced inequality comes from the definition of $\Gamma_D \ni (x, y)$ and the second comes from the fact that $M$ upper bounds the Lipschitz constant of $f$. Hence, defining:

$$ \rho = \frac{\eta}{2M \|\partial_x v_\eta\|_\infty}, $$
we have just proved that if for all $x \in X$, there exists an $i \in [N]$ with $\|x - x_i\| \le \rho$, then:

$$ 
\begin{aligned}
\mathcal{L}_y v_\eta(t, x) &= \underbrace{\mathcal{L}_f v_\eta(t, x)}_{\le -\eta} + (y - f(x))^{\top} \partial_x v_\eta(t, x) \\
&\le \|y - f(x)\| \cdot \|\partial_x v_\eta(t, x)\| - \eta \\
&\le 2M \|x - x_i\| \cdot \|\partial_x v_\eta\|_{\infty} - \eta \\
&\le 2M\rho \|\partial_x v_\eta\|_{\infty} - \eta = 0.
\end{aligned}
$$

Since $X$ is compact, it can be covered by a finite number of balls of radius $\rho$. Hence, sampling $f$ at the centres of these balls gives an appropriate dataset $D = \{(x_i, f(x_i))\}_{i=1}^N$ such that for all $x \in X$ there is an index $i \in [N]$ with $\|x - x_i\| \le \rho$. For this specific dataset, we have proved that $(v_\eta, w_\eta)$ is feasible for both problems (18) and (12).

Hence, it holds $\tilde{X}_0(w_\eta) \subset X_0^*(D) \subset X_0(f)$, and

$$ 
\begin{aligned}
\mathrm{vol}(X_0(f) \setminus \tilde{X}_0(w_\eta)) &\le \int_X w_\eta(x) \, dx - \mathrm{vol}(X \setminus X_0(f)) \\
&= \underbrace{\int_X w_\varepsilon(x) \, dx - \mathrm{vol}(X \setminus X_0(f))}_{\le \varepsilon/2} \\
&\quad + \eta \underbrace{\int_X w_\delta(x) \, dx}_{=C} \le \varepsilon/2 + C \cdot \eta = \varepsilon
\end{aligned}
$$

so that $X_0(f)$ is approximated arbitrarily closely (from inside) by a feasible function of the data-based problem (12), with no knowledge of the full dynamics. □

**Remark 12** What we actually proved is that, up to sufficiently informative datasets $D_N$ and high relaxation orders in the hierarchy (see Appendix A), one can approximate the true FH-RoA from inside arbitrarily closely, by solving SDP problems. Although less meaningful, the result also stands for best-case outer approximation (6) with:

$$ \mathrm{vol}(X_0(D_N) \setminus X_0(f)) \xrightarrow[N \to \infty]{} 0, $$

and the proof is identical to the above.

**Remark 13** The above proof highlights how crucial the validity and quality of the Lipschitz upper bound $M$ is. If $M$ underestimates the Lipschitz constant, then $f(x)$ can take values outside of $F_D(x)$, hence hampering the upper bound on $\|y - f(x)\|$ in the proof. If $M$ is very large, then the radius $\rho$ in the proof becomes very small, hence asking for a much larger dataset to obtain the appropriate covering of $X$. In short, the larger $M$, the more conservative the FH-RoA approximation. Methods are proposed to estimate and validate $M$ from data in [11, 18]. In this paper, we only assume that $M$ is known in advance.