# Data-driven stabilizer design and closed-loop analysis of general nonlinear
systems via Taylor’s expansion

Meichen Guo<EMAIL_ADDRESS>Claudio De Persis<EMAIL_ADDRESS>Pietro
Tesi<EMAIL_ADDRESS>Delft Center for Systems and Control, Delft
University of Technology, 2628 CD Delft, The Netherlands  Engineering and
Technology Institute Groningen, University of Groningen, Nijenborgh 4, 9747 AG
Groningen, The Netherlands  Department of Information Engineering, University
of Florence, 50139 Florence, Italy

###### Abstract

For data-driven control of nonlinear systems, the basis functions
characterizing the dynamics are usually essential. In existing works, the
basis functions are often carefully chosen based on pre-knowledge of the
dynamics so that the system can be expressed or well-approximated by the basis
functions and the experimental data. For a more general setting where explicit
information on the basis functions is not available, this paper presents a
data-driven approach for stabilizer design and closed-loop analysis via the
Lyapunov method. First, based on Taylor’s expansion and using input-state
data, a stabilizer and a Lyapunov function are designed to render the known
equilibrium locally asymptotically stable. Then, data-driven conditions are
derived to check whether a given sublevel set of the found Lyapunov function
is an invariant subset of the region of attraction. One of the main challenges
is how to handle Taylor’s remainder in the design of the local stabilizers and
the analysis of the closed-loop performance.

###### keywords:

data-driven control , nonlinear control , region of attraction estimation ,
non-polynomial systems

##  1 Introduction

Most control approaches of nonlinear systems are based on well-established
models of the system constructed by pre-knowledge or system identification.
When the models are not explicitly constructed, nonlinear systems can be
directly controlled using input-output data. Controlling a system via input-
output data without explicitly identifying the model is called the direct
data-driven control method, and it has been gaining more and more attentions
for both linear and nonlinear systems. An early survey of data-driven control
methods can be found in [1]. More recently, the authors of [2] developed an
online control approach, the work [3] utilized the dynamic linearization data
models for discrete-time non-affine nonlinear systems, the authors of [4] and
[5] considered feedback linearizable systems, and the works [6] and [7]
designed data-driven model predictive controllers. Inspired by Willems _et
al._ ’s fundamental lemma, [8] proposed data-driven control approaches for
linear and nonlinear discrete-time systems. Using a matrix Finsler’s Lemma,
[9] applied data-driven control to Lur’e systems. The authors of [10] used
state-dependent representation and proposed an online optimization method for
data-driven stabilization of nonlinear dynamics. For polynomial systems, [11]
designed global stabilizers using noisy data, and [12] synthesized data-driven
safety controllers. The recent work [13] investigated dissipativity of
nonlinear systems based on polynomial approximation.

_Related works._ Some recent works related to nonlinear data-driven control
and the region of attraction (RoA) estimation are discussed in what follows.

Deriving a data-based representation of the dynamics is one of the important
steps in data-driven control of unknown nonlinear systems. If the controlled
systems are of certain classes, such as polynomial systems having a known
degree, the monomials of the state can be chosen as basis functions to design
data-driven controllers such as presented in [14, 11]. By integrating noisy
data and side information, [15] showed that unknown polynomial dynamics can be
learned via semidefinite programming. When the nonlinearities satisfy
quadratic constraints, data-driven stabilizer was developed in [16]. With
certain knowledge and assumptions on the nonlinear basis function, systems
containing more general types of nonlinearities have also been studied in
recent works. For instance, under suitable conditions, some nonlinear systems
can be lifted into polynomial systems in an extended state for control, such
as the results shown in [17, Section IV] and [18, Section 3.2]. Using
knowledge of the basis functions, [19] designed data-driven controllers by
(approximate) cancellation of the nonlinearity. When the system nonlinearities
cannot be expressed as a combination of known functions, [19] presented data-
driven local stabilization results by choosing basis functions carefully such
that the neglected nonlinearities are small in a known set of the state. On
the other hand, if the knowledge on the basis functions is not available,
approximations of the nonlinear systems are often involved. The previous work
[8] tackled the nonlinear data-driven control problem by linearizing the
dynamics around the known equilibrium and obtaining a local stability result.
According to these existing results, it is clear that the efficiency and the
performance of data-driven controllers can be improved via pre-known knowledge
such as specific classes of the systems or the nonlinear basis functions.
Nonetheless, there is still a lack of comprehensive investigation of the more
general case where the nonlinear basis functions cannot be easily and
explicitly obtained.

The RoA estimation is another relevant topic in nonlinear control. For general
nonlinear systems, it is common that the designed controllers only guarantee
local stability. Hence, it is of importance to obtain the RoA of the closed-
loop systems for the purpose of theoretical analysis as well as engineering
applications. Unfortunately, for general nonlinear systems, it is extremely
difficult to derive the exact RoA even when the model is explicitly known. A
common solution is to estimate the RoA based on Lyapunov functions. Using
Taylor’s expansion and considering the worst-case remainders, [20] estimated
the RoA of uncertain non-polynomial systems via linear matrix inequality (LMI)
optimizations. RoA analysis for polynomial systems was presented in [21] using
polynomial Lyapunov functions and sum of squares (SOS) optimizations. The
authors of [22] studied uncertain nonlinear systems subject to perturbations
in certain forms and used the SOS technique to compute invariant subsets of
the RoA. It is noted that, in these works, the RoA estimation winds up in
solving bilinear optimization problems, and techniques such as bisection or
special bilinear inequality solver are required to find the solutions. For
nonlinear systems without explicit models, there are also efforts devoted to
learning the RoA by various approaches. The authors of [23] developed a
sampling-based approach for a class of piecewise continuous nonlinear systems
that verifies stability and estimates the RoA using Lyapunov functions. Based
on the converse Lyapunov theorem, [24] processed system trajectories to lift a
Lyapunov function whose level sets lead to an estimation of the RoA. Using the
properties of recurrent sets, [25] proposed an approach that learns an inner
approximation of the RoA via finite-length trajectories. It should be pointed
out that, all the aforementioned works focus on stability analysis of
autonomous nonlinear systems, _i.e._ , the control design is not considered.

_Contributions._ For general nonlinear systems, this paper presents a data-
driven approach to simultaneously obtaining a Lyapunov function and designing
a state feedback stabilizer that renders the known equilibrium locally
asymptotically stable. Using Taylor’s expansion, the unknown dynamics are
approximated by linear systems or polynomial systems. Then, linear stabilizers
and polynomial stabilizers are designed for the approximated models using
finite-length input-state data collected in an off-line experiment. To handle
the remainder resulting from the approximation, we conduct the experiment
close to the known equilibrium such that the remainder is small with a known
bound. An over-approximation of all the feasible dynamics is then found using
the collected data, and Petersen’s lemma [26] is used for the controller
design. For polynomial approximations, as the conditions characterizing the
stabilizers are positive conditions for polynomial matrices, the SOS technique
[27] is utilized to make the conditions easily solvable. The data-driven
stabilizer design can be seen as a generalization of the nonlinear control
result in the previous work [8, Section V.B]. Specifically, this paper
considers both the linear and the polynomial approximations of continuous-time
systems. Analogous results can also be derived for discrete-time systems. On
the other hand, the focus of this paper is on general nonlinear systems, while
the previous works [11] and [28] dealt with only polynomial systems. In
comparison with [19] where the influence of the remaining nonlinearity is
attenuated by careful choices of the basis functions, this work achieves the
objective by using Taylor’s polynomials as basis functions and conduct the
experiment close to the equilibrium.

For estimating the RoA with the designed data-driven controller, we first
derive an estimation of the remainder by assuming a known bound on the high-
order derivatives of the unknown functions. With the help of the
Positivstellensatz [29] and the SOS relaxations, we derive data-driven
conditions that verify whether a given sublevel set of the obtained Lyapunov
function is an invariant subset of the RoA. This is achieved by deriving a
sufficient condition for the negativity of the derivative of the Lyapunov
function based on the estimation of the remainder. The conditions are derived
via data and some prior knowledge on the dynamics, and can be easily solved by
software such as Matlab. The estimated RoA gives insights to the closed-loop
system under the designed data-driven controller, and is relevant for both
theoretical and practical purposes. Note that alternatively, the RoA can be
estimated based on other data-driven methods, such as the ones developed in
[23, 24, 25]. Simulations results on the inverted pendulum show the
applicability of the design and estimation approach. For future works, it is
of our interests to thoroughly investigate the influence of certain parameters
and further improve the design and analysis.

The rest of the paper is arranged as follows. Section 2 formulates the problem
and then presents relevant techniques and existing results for the subsequent
sections. Data-driven control designs with different orders of approximations
are presented in Section 3. The data-driven characterization of the RoA is
derived in Section 4. Numerical results and analysis on the inverted pendulum
are illustrated in Section 5. Finally, Section 6 concludes the paper.

_Notation._ Throughout the paper, $A\succ(\succeq)0$ denotes that matrix $A$
is positive (semi-)definite, and $A\prec(\preceq)0$ denotes that matrix $A$ is
negative (semi-)definite. For vectors $a,b\in\mathbb{R}^{n}$, $a\preceq b$
means that $a_{i}\leq b_{i}$ for all $i=1,\dots,n$. $\|\cdot\|$ denotes the
Euclidean norm.

Moreover, we list some important symbols used in subsequent sections in the
table below.

$R(x,u)$ | Taylor’s remainder of $f(x,u)$
---|---
$S=\begin{bmatrix}B&A\end{bmatrix}$ | True dynamics of the first-/high-order
| approximation
$\widehat{S}=\begin{bmatrix}\widehat{B}&\widehat{A}\end{bmatrix}$ | Dynamics consistent with the collect-
| ed data
$\mathcal{C}_{k}$ | Feasible set of the dynamics at each
| time $t_{k}$
$\textbf{A}_{k}$, $\textbf{B}_{k}$, $\textbf{C}_{k}$ | Matrices describing set $\mathcal{C}_{k}$
$\mathcal{I}$ | Intersection of the sets $\mathcal{C}_{k}$,
| $k=0,\dots,T-1$
$\overline{\mathcal{I}}$ | Over-approximation of the set $\mathcal{I}$
$\overline{\textbf{A}}$, $\overline{\textbf{B}}$, $\overline{\textbf{C}}$ | Matrices describing set $\overline{\mathcal{I}}$

##  2 Problem Formulation and Preliminaries

Consider a general nonlinear continuous-time system

$\displaystyle\dot{x}=f(x,u)$ (1)

where the state $x\in\mathbb{R}^{n}$ and the input $u\in\mathbb{R}^{m}$. The
function $f$ is of class $C^{r}$ for some integer $r\geq 0$. Assume that
$(x_{e},u_{e})$ is a known equilibrium of the system to be stabilized. For
simplicity and without loss of generality, in this paper we let
$(x_{e},u_{e})=(0,0)$, as any known equilibrium can be converted to the origin
by a change of coordinates.

For a general nonlinear discrete-time system, [8, Section V.B] approximates
the unknown dynamics as a linear system around a known equilibrium and
developed a data-driven local stabilizer. In this work, we show that the
result in therein can be extended to approximations with any order via
Taylor’s expansion. Moreover, for the closed-loop system under the designed
data-driven controller, a data-driven estimation of the RoA will also be
presented.

To gather information regarding the system, we perform $E$ experiments on the
system, where $E$ is an integer satisfying $1\leq E\leq T$ and $T$ an integer
equalling the total number of collected samples. On one extreme, one could
perform $1$ single experiment during which a total of $T$ samples are
collected. On the other extreme, one could perform $T$ independent
experiments, during each one of which a single sample is collected. The
advantage of short multiple experiments is that they allow the designer to
collect information about the system at different points in the state space
without incurring in problems due to the free evolution of the system.

Either way, a dataset
$\mathcal{DS}_{c}:=\\{\dot{x}(t_{k});x(t_{k});u(t_{k})\\}^{T-1}_{k=0}$ for the
continuous-time system can be obtained. Organize the data collected in the
experiment(s) as

$\displaystyle X_{0}$
$\displaystyle=\begin{bmatrix}x(t_{0})&\cdots&x(t_{T-1})\end{bmatrix},$
$\displaystyle U_{0}$
$\displaystyle=\begin{bmatrix}u(t_{0})&\cdots&u(t_{T-1})\end{bmatrix},$
$\displaystyle X_{1}$
$\displaystyle=\begin{bmatrix}\dot{x}(t_{0})&\cdots&\dot{x}(t_{T-1})\end{bmatrix}.$

The problem studied in this work is formulated as follows.

######  Problem 1 (_Data-driven stabilizer design and RoA estimation_).

For system (1), design a state feedback controller $u=k(x)$ using the input-
state data $X_{0}$, $X_{1}$, and $U_{0}$, such that the origin is locally
asymptotically stable for the closed-loop system, and an inner estimation of
the RoA of the closed-loop system is derived.

######  Remark 1 (_On the experimental data_).

The derivative data $\dot{x}(t_{k})$, $k=1,\dots,T-1$ can be approximated
using numerical differentiation. For instance, using the forward difference
approximation gives

$\dot{x}_{i}(t_{k})=\frac{x_{i}(t_{k+1})-x_{i}(t_{k})}{t_{k+1}-t_{k}}+e_{i}(t_{k}),~{}~{}i=1,\dots,N,$

where the approximation error $e_{i}(t_{k})$ is proportional to
$t_{k+1}-t_{k}$ and can be handled as noise.

In this work, the data is assumed to be noiseless for the sake of simplicity.
If the data is corrupted by noise, to decrease the effect of noise in the
derivative approximation, numerical approaches such as the total variation
regularization [30] can be used as shown in [31] and [32].
$\hfill\blacksquare$

The rest of this section presents techniques and existing results needed for
the subsequent data-driven control design and analysis.

###  2.1 Taylor’s expansion

To approximate the unknown nonlinear dynamics, we use Taylor’s expansion to
represent the unknown function by the sum of the $r$th Taylor polynomial and
the remainder for any integer $r\geq 0$ at a given equilibrium. As explained
previously, we assume that the functions in consideration have an equilibrium
at the origin.

Consider a function $\phi:\mathbb{R}^{\sigma}\rightarrow\mathbb{R}$ that is of
class $C^{r}$ on an open convex set $\mathbb{D}\subseteq\mathbb{R}^{n}$
containing the origin, and $\phi(0)=0$. Define
$\alpha=(\alpha_{1},\dots,\alpha_{\sigma})$ as a $\sigma$-tuple of nonnegative
integers and

$\displaystyle|\alpha|$ $\displaystyle=\alpha_{1}+\cdots+\alpha_{\sigma},$
$\displaystyle\alpha!$ $\displaystyle=\alpha_{1}!\cdots\alpha_{\sigma}!,$
$\displaystyle~{}~{}\partial^{\alpha}\phi$
$\displaystyle=\frac{\partial^{|\alpha|}\phi}{\partial
z_{1}^{\alpha_{1}}\cdots\partial z_{\sigma}^{\alpha_{\sigma}}},$
$\displaystyle z^{\alpha}$ $\displaystyle=z_{1}^{\alpha_{1}}\cdots
z_{\sigma}^{\alpha_{\sigma}}$

for any $z=[z_{1}~{}\cdots~{}z_{\sigma}]^{\top}\in\mathbb{R}^{\sigma}$. As
given in [33], for $z\in\mathbb{D}$, Taylor’s expansion of the function $\phi$
at the origin is

$\displaystyle\phi(z)$ $\displaystyle=\sum_{|\alpha|\leq
r}\frac{\partial^{\alpha}\phi(0)}{\alpha!}z^{\alpha}+R_{r}(z)$ (2)

for any integer $r\geq 0$, and the remainder can be expressed in integral form
as

$\displaystyle
R_{r}(z)=(r+1)\sum_{|\alpha|=r+1}\frac{z^{\alpha}}{\alpha!}\int^{1}_{0}(1-t)^{r}\partial^{\alpha}\phi(tz)dt.$
(3)

In [34], the author has shown that the estimate of the modulus of
$\partial\phi^{\alpha}(z)$, $|\alpha|=r$ at the origin can be used to obtain
an estimate for the remainder $R_{r}(z)$.

######  Lemma 1.

Consider all $z\in\mathbb{D}$ where $\mathbb{D}\subseteq\mathbb{R}^{\sigma}$
is an open convex set containing the origin. If for $|\alpha|=r$ with any
integer $r\geq 0$, $\phi$ is of class $C^{r}$ for all $z\in\mathbb{D}$, and
there exists $C\geq 0$ such that

$\displaystyle|\partial^{\alpha}\phi(z)-\partial^{\alpha}\phi(0)|\leq
C\|z\|~{}~{}\forall z\in\mathbb{D},$ (4)

then the remainder $R_{r}(z)$ in (3) satisfies

$\displaystyle\left|R_{r}(z)\right|\leq\frac{\sigma^{r/2}C\|z\|^{r+1}}{(r+1)!}~{}~{}\forall
z\in\mathbb{D}.$ (5)

The proof of Lemma 1 can be found in A.

###  2.2 Petersen’s lemma

In the section for data-driven controller design, Petersen’s lemma is
essential for deriving the sufficient condition characterizing the controller.
Due to the space limit, the proof of the lemma is omitted and one may refer to
works such as [26, 35, 28] for more details.

######  Lemma 2 (Petersen’s lemma [26]).

Consider matrices $\mathcal{G}=\mathcal{G}^{\top}\in\mathbb{R}^{n\times n}$,
$\mathcal{M}\in\mathbb{R}^{n\times m}$, $\mathcal{M}\neq 0$,
$\mathcal{N}\in\mathbb{R}^{p\times n}$, $\mathcal{N}\neq 0$, and a set $F$
defined as

$F=\\{\mathcal{F}\in\mathbb{R}^{m\times
p}:\mathcal{F}^{\top}\mathcal{F}\preceq\overline{\mathcal{F}}\\}$

where $\overline{\mathcal{F}}=\overline{\mathcal{F}}^{\top}\succeq 0$. Then,
for all $\mathcal{F}\in F$

$\mathcal{G}+\mathcal{M}\mathcal{F}\mathcal{N}+\mathcal{N}^{\top}\mathcal{F}^{\top}\mathcal{M}^{\top}\preceq
0$

holds if and only if there exists $\mu>0$ such that

$\mathcal{G}+\mu\mathcal{M}\mathcal{M}^{\top}+\mu^{-1}\mathcal{N}^{\top}\overline{\mathcal{F}}\mathcal{N}\preceq
0.$

###  2.3 Sum of squares relaxation

As solving positive conditions of multivariable polynomials is in general NP-
hard, the SOS relaxations are often used to obtain sufficient conditions that
are tractable. The SOS polynomial matrices are defined as follows.

######  Definition 1.

(_SOS polynomial matrix_ [36])
$M:\mathbb{R}^{n}\rightarrow\mathbb{R}^{\sigma\times\sigma}$ is an SOS
polynomial matrix if there exist
$M_{1},\dots,M_{k}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{\sigma\times\sigma}$
such that

$\displaystyle M(x)=\sum^{k}_{i=1}M_{i}(x)^{\top}M_{i}(x)~{}~{}\forall
x\in\mathbb{R}^{n}.$ (6)

Note that when $\sigma=1$, $M(x)$ becomes a scalar SOS polynomial.

It is straightforward to see that if a matrix $M(x)$ is an SOS polynomial
matrix, then it is positive semi-definite, _i.e._ , $M(x)\succeq 0~{}\forall
x\in\mathbb{R}^{n}$. Relaxing the positive polynomial conditions into SOS
polynomial conditions makes the conditions tractable and easily solvable by
common software.

###  2.4 Positivstellensatz

In the RoA analysis, we need to characterize polynomials that are positive on
a semialgebraic set and the Positivstellensatz plays an important role in the
characterization.

Let $p_{1},\dots,p_{k}$ be polynomials. The multiplicative monoid, denoted by
$\mathcal{S}_{M}(p_{1},\dots,p_{k})$, is the set generated by taking finite
products of the polynomials $p_{1},\dots,p_{k}$. The cone
$\mathcal{S}_{C}(p_{1},\dots,p_{k})$ generated by the polynomials is defined
as

$\displaystyle\mathcal{S}_{C}(p_{1},\dots,p_{k})$
$\displaystyle=\\{s_{0}+\sum_{i=1}^{j}s_{i}q_{i}:s_{0},\dots,s_{j}\text{ are
SOS polynomials},$ $\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad
q_{1},.\dots,q_{j}\in\mathcal{S}_{M}(p_{1},\dots,p_{k})\\}.$

The ideal $\mathcal{S}_{I}(p_{1},\dots,p_{k})$ generated by the polynomials is
defined as

$\mathcal{S}_{I}(p_{1},\dots,p_{k})=\left\\{\sum_{i=1}^{k}r_{i}p_{i}:r_{1},\dots,r_{k}\text{
are polynomials}\right\\}.$

Stengle’s Positivstellensatz [29] is presented as follows in [36].

######  Theorem 1 (_Positivstellensatz_).

Let $f_{1},\dots,f_{k}$, $g_{1},\dots,g_{l}$, and $h_{1},\dots,h_{m}$ be
polynomials. Define the set

$\displaystyle\mathcal{X}=\\{x\in\mathbb{R}^{n}:$ $\displaystyle f_{1}(x)\geq
0,\dots,f_{k}(x)\geq 0,$ $\displaystyle g_{1}(x)=0,\dots,g_{l}(x)=0,$
$\displaystyle\text{and }h_{1}(x)\neq 0,\dots,h_{m}(x)\neq 0\\}.$

Then, $\mathcal{X}=\emptyset$ if and only if

$\exists
f\\!\in\mathcal{S}_{C}(f_{1},\dots,f_{k}),g\\!\in\mathcal{S}_{I}(g_{1},\dots,g_{l}),h\\!\in\mathcal{S}_{M}(h_{1},\dots,h_{m})$

such that

$f(x)+g(x)+h(x)^{2}=0.$

For the subsequent RoA analysis, we will use the following result derived from
the Positivstellensatz.

######  Lemma 3.

Let $\varphi_{1}$ and $\varphi_{2}$ be polynomials in $x$. If there exist SOS
polynomials $s_{1}$ and $s_{2}$ in $x$ such that

$\displaystyle-(s_{1}\varphi_{1}(x)+s_{2}\varphi_{2}(x)+x^{\mbox{\tiny\sf
T}}x)\text{ is SOS }~{}\forall x\in\mathbb{R}^{n}$ (7)

then the set inclusion condition

$\displaystyle\\{x\in\mathbb{R}^{n}:\varphi_{1}(x)\geq 0,x\neq
0\\}\subseteq\\{x\in\mathbb{R}^{n}:\varphi_{2}(x)<0\\}$ (8)

holds.

The proof of Lemma 1 can be found in B.

##  3 Data-driven controller design

To approximate the system (1), we write Taylor’s expansion of the function
$f(x,u)$ as

$\displaystyle f(x,u)=f_{appr}(x,u)+R(x,u),$

where $f_{appr}(x,u)$ contains Taylor’s polynomials up to a certain degree and
$R(x,u)$ represents the truncated remainder containing high-order terms. This
remainder constitutes the main uncertainty of the data-based representation of
the closed-loop system and brings difficulties in the controller design. To
attenuate the influence of the remainder, we collect state data close to the
equilibrium such that the maximum amplitude of the remainder during the
experiment is instantaneously bounded by a known constant. While the first-
order approximation of the system leads to a solution to Problem 1, to further
minimizing the remainder in the neighborhood of the equilibrium, we can
truncate Taylor’s polynomials at a higher degree and obtain a high-order
approximation.

Using the collected data, we first find an over-approximation of the set
containing all dynamics that are consistent with the data. Then, stabilizers
are designed such that the origin is locally asymptotically stable for all the
dynamics in the approximated set.

###  3.1 First-order approximation

First, for completeness, we will address the local data-driven stabilizer
design via first-order approximation.

Consider the continuous-time nonlinear system (1). Denote each element of $f$
as $f_{i}$ and let $f_{i}\in C^{1}(\mathbb{R}^{n}\times\mathbb{R}^{m})$,
$i=1,\dots,n$. The first-order approximation of (1) through Taylor’s expansion
of $f(x,u)$ is

$\displaystyle\dot{x}=Ax+Bu+R(x,u)$ (9)

where $R(x,u)$ denotes the remainder and

$\displaystyle A=\left.\frac{\partial f(x,u)}{\partial
x}\right|_{(x,u)=(0,0)},~{}~{}B=\left.\frac{\partial f(x,u)}{\partial
u}\right|_{(x,u)=(0,0)}.$

One can treat the remainder as a disturbance that affects the data-driven
characterization of the unknown dynamics, and focus on controlling the
approximated linear dynamics to obtain a local linear stabilizer. Then, by
tackling the impact of $R(x,u)$, the RoA of the closed-loop system can be also
be characterized.

Similar to our previous works such as [8] and [11], some bound on the
experimental data of the remainder is needed to design a controller. For this
purpose, the following assumption is posed.

######  Assumption 1.

For $k=0,\dots,T$ and a known $\gamma$,

$\displaystyle R(x(t_{k}),u(t_{k}))^{\top}R(x(t_{k}),u(t_{k}))\leq\gamma^{2}.$
(10)

######  Remark 2 (_Instantaneous remainder bound_).

Assumption 1 is an instantaneous bound on the maximum amplitude of the
remainder during the experiment. The bound can be obtained by prior knowledge
of the model, such as the physics of the system. If such knowledge is
unavailable, one may resort to an over-estimation of $\gamma$.
$\hfill\blacksquare$

######  Remark 3 (_Noisy data._).

If the data is corrupted by additive measurement and/or actuator noise, _i.e._
, $x(t_{k})=x^{*}(t_{k})+d_{x}(t_{k})$ and
$u(t_{k})=u^{*}(t_{k})+d_{u}(t_{k})$ where $x^{*}(t_{k})$, $u^{*}(t_{k})$ are
the true data, and $d_{x}(t_{k})$, $d_{u}(t_{k})$ represent the noise, the
proposed design and analysis approach is still applicable to the noisy case,
provided that the remainder $R(x(t_{k}),u(t_{k}))$ in Assumption 1 is replaced
by $R(x^{*}(t_{k}),u(t_{k}))+d(t_{k})$, where $d(t_{k})$ is the total noise
due to $d_{x}(t_{k})$, $d_{u}(t_{k})$, and the derivative approximation error.
$\hfill\blacksquare$

####  3.1.1 Over-approximation of the feasible set

Under Assumption 1, an over-approximation of the set of dynamics that are
consistent with the experimental data can be derived as shown in [37].

Denote $S=[B~{}~{}A]$ as the true dynamics. Based on (9), at each time
$t_{k}$, $k=0,\dots,T-1$, the collected data satisfies

$\displaystyle\dot{x}(t_{k})=S\begin{bmatrix}u(t_{k})\\\
x(t_{k})\end{bmatrix}+R\big{(}x(t_{k}),u(t_{k})\big{)}.$

Under Assumption 1, one has that

$\displaystyle\quad
R\big{(}x(t_{k}),u(t_{k})\big{)}R\big{(}x(t_{k}),u(t_{k})\big{)}^{\top}$
$\displaystyle=\left(\dot{x}(t_{k})-S\begin{bmatrix}u(t_{k})\\\
x(t_{k})\end{bmatrix}\right)\left(\dot{x}(t_{k})-S\begin{bmatrix}u(t_{k})\\\
x(t_{k})\end{bmatrix}\right)^{\top}\preceq\gamma^{2}I.$

Hence, at each time $t_{k}$, $k=0,\dots,T-1$, the matrices
$\widehat{S}=[\widehat{B}~{}~{}\widehat{A}]$ consistent with the data belongs
to the set

$\displaystyle\mathcal{C}_{k}=\left\\{\widehat{S}:\textbf{C}_{k}+\widehat{S}\textbf{B}_{k}+\textbf{B}_{k}^{\top}\widehat{S}^{\top}+\widehat{S}\textbf{A}_{k}\widehat{S}^{\top}\preceq
0\right\\}$ (11)

where

$\displaystyle\textbf{A}_{k}$
$\displaystyle=l(t_{k})l(t_{k})^{\top},~{}~{}\textbf{B}_{k}=-l(t_{k})\dot{x}(t_{k})^{\top},$
$\displaystyle\textbf{C}_{k}$
$\displaystyle=\dot{x}(t_{k})\dot{x}(t_{k})^{\top}-\gamma^{2}I,~{}~{}l(t_{k})=\begin{bmatrix}u(t_{k})\\\
x(t_{k})\end{bmatrix}.$

Then, the feasible set of matrices $\widehat{S}$ that is consistent with all
data collected in the experiment(s) is the intersection of all the sets
$\mathcal{C}_{k}$, _i.e._ , $\mathcal{I}=\bigcap^{T-1}_{k=0}\mathcal{C}_{k}$.
Though the exact set $\mathcal{I}$ is difficult to obtain, an over-
approximation of $\mathcal{I}$ in the form of a matrix ellipsoid and of
minimum size can be computed. Denote the over-approximation set as

$\displaystyle\overline{\mathcal{I}}:=\left\\{\widehat{S}:\overline{\textbf{C}}+\widehat{S}\overline{\textbf{B}}+\overline{\textbf{B}}^{\top}\widehat{S}^{\top}+\widehat{S}\overline{\textbf{A}}\widehat{S}^{\top}\preceq
0\right\\}$ (12)

where $\overline{\textbf{A}}=\overline{\textbf{A}}^{\top}\succ 0$,
$\overline{\textbf{C}}$ is set as
$\overline{\textbf{C}}=\overline{\textbf{B}}^{\top}\overline{\textbf{A}}^{-1}\overline{\textbf{B}}^{\top}-\delta
I$ and $\delta>0$ is a constant fixed arbitrarily. Following [37, Section
5.1], the set $\overline{\mathcal{I}}$ can be found by solving the
optimization problem

$\displaystyle\underset{\overline{\textbf{A}},\overline{\textbf{B}},\overline{\textbf{C}}}{\text{minimize}}\quad$
$\displaystyle-\mathrm{log~{}det}(\overline{\textbf{A}})$ (13) subject to
$\displaystyle\overline{\textbf{A}}=\overline{\textbf{A}}^{\top}\succ 0$
$\displaystyle\tau_{k}\geq 0,~{}k=0,\dots,T-1$
$\displaystyle\begin{bmatrix}-\delta
I-\sum\limits^{T-1}_{k=0}\tau_{k}\textbf{C}_{k}&\overline{\textbf{B}}^{\top}-\sum\limits^{T-1}_{k=0}\tau_{k}\textbf{B}_{k}^{\top}&\overline{\textbf{B}}^{\top}\\\
\overline{\textbf{B}}-\sum\limits^{T-1}_{k=0}\tau_{k}\textbf{B}_{k}&\overline{\textbf{A}}-\sum\limits^{T-1}_{k=0}\tau_{k}\textbf{A}_{k}&0\\\
\overline{\textbf{B}}&0&-\overline{\textbf{A}}\end{bmatrix}\preceq 0.$

######  Remark 4 (_On the parameter $\delta$_).

To make sure that the over-approximation set $\overline{\mathcal{I}}$ is not
empty, we need that
$\overline{\textbf{B}}^{\top}\overline{\textbf{A}}^{-1}\overline{\textbf{B}}^{\top}-\overline{\textbf{C}}\succ
0$, which is guaranteed by setting
$\overline{\textbf{C}}=\overline{\textbf{B}}^{\top}\overline{\textbf{A}}^{-1}\overline{\textbf{B}}^{\top}-\delta
I$ for some positive constant $\delta$. As pointed out in [38, Section 3.7],
the description of the set $\overline{\mathcal{I}}$ is homogeneous, that is,
the matrices $\overline{\textbf{A}}$, $\overline{\textbf{B}}$, and
$\overline{\textbf{C}}$ can be scaled by any positive factor without affecting
$\overline{\mathcal{I}}$. Hence, the positive constant $\delta$ can be fixed
arbitrarily without changing the resulting $\overline{\mathcal{I}}$. In many
works, such as [38] and [37], the variables are normalized, _i.e._ , $\delta$
is set as $1$. In this work, we keep the parameter $\delta$ because although
it does not affect the over-approximation of the feasible set, it has effects
on the subsequent RoA estimation as shown in the simulation results. However,
it is still unclear to us how $\delta$ affects the RoA estimation
quantitatively or how to choose an optimal $\delta$ for the design and
analysis. These questions can only be answered after more thorough studies on
the parameters in the control design process and the RoA analysis, which is
out of the scope of this work. $\hfill\blacksquare$

######  Remark 5 (_Persistency of excitation_).

As pointed out in [37, Section 3.1], if the collected data is rich enough,
_i.e._ , $\begin{bmatrix}U_{0}\\\ X_{0}\end{bmatrix}$ has full row rank, then
the intersection set $\mathcal{I}$ is bounded, which allows the optimization
problem (13) to have a solution. Hence, $\begin{bmatrix}U_{0}\\\
X_{0}\end{bmatrix}$ having full row rank implies the feasibility of (13).

####  3.1.2 Data-driven stabilizer design

Stabilizing the linear approximation of the unknown system renders the origin
locally asymptotically stable as the remainder $R(x,u)$ contains higher-order
terms and converges to the origin faster than the linear part in a
neighborhood of the origin. Hence, the objective of the controller design is
to stabilize the origin for all dynamics belonging to the over-approximation
set $\overline{\mathcal{I}}$. This can be achieved in the same manner as done
in [28] via Petersen’s lemma. For the completeness of this work, we give the
following result on designing a data-driven local stabilizer. The proof shares
the same idea with Theorem 2 in [28], and thus is omitted in this paper.

######  Theorem 2.

Under Assumption 1, given a constant $w>0$, if there exist matrices $Y$ and
$P=P^{\top}$ such that

$\displaystyle\begin{bmatrix}wP-\overline{\textbf{C}}&\overline{\textbf{B}}^{\top}-\begin{bmatrix}Y\\\
P\end{bmatrix}^{\top}\\\ \overline{\textbf{B}}-\begin{bmatrix}Y\\\
P\end{bmatrix}&-\overline{\textbf{A}}\end{bmatrix}$ $\displaystyle\preceq 0$
(14) $\displaystyle P$ $\displaystyle\succ 0,$

then the origin is a locally asymptotically stable equilibrium for the closed-
loop system composed of (1) and the control law $u=YP^{-1}x$.

Consider the Lyapunov function $V(x)=x^{\top}P^{-1}x$. The controller designed
by Theorem 2 guarantees that the derivative of $V(x)$ along the trajectory of
the closed-loop system $\dot{x}=\widehat{A}x+\widehat{B}u$ satisfies that

$\dot{V}(x)\leq-wV(x)$

for any given constant $w>0$ and any
$\begin{bmatrix}\widehat{B}&\widehat{A}\end{bmatrix}\in\overline{\mathcal{I}}$.
Hence, by choosing the value of $w$, a certain decay rate of the closed-loop
solution is enforced.

###  3.2 High-order approximation

One may also approximate the nonlinear system as a polynomial system by
truncating Taylor’s expansion at a degree higher than $1$. This will result in
a smaller remainder in the neighborhood close to the origin. Similar to the
first-order approximation case, by making sure that the remainder converges to
the origin faster, a stabilizer for the approximated system can render the
origin locally asymptotically stable for the overall system.

To write the approximated system into a linear-like form, we consider the
nonlinear input-affine system

$\displaystyle\dot{x}=f(x)+g(x)u$ (15)

where $x\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$, and $f(0)=0$. Functions $f$
and $g$ are of class $C^{r_{f}}$ and $C^{r_{g}}$ respectively for some
integers $r_{f},r_{g}\geq 1$. Using Taylor’s expansion, we truncate the series
of polynomials at order $r_{f}$ for functions $f_{i}$, and at order $r_{g}$
for functions $g_{ij}$, $i=1,\dots,n$, $j=1,\dots,m$, respectively. Following
Section 2.1, the functions can be written as

$\displaystyle f_{i}(x)$ $\displaystyle=\sum_{|\alpha|\leq
r_{f}}\frac{\partial^{\alpha}f_{i}(0)}{\alpha!}x^{\alpha}+R_{f_{i}}(x),$
$\displaystyle g_{ij}(x)$ $\displaystyle=\sum_{|\alpha|\leq
r_{g}}\frac{\partial^{\alpha}g_{ij}(0)}{\alpha!}x^{\alpha}+R_{g_{ij}}(x)$

with remainders $R_{f_{i}}(x)$ and $R_{g_{ij}}(x)$. Then, one can write the
polynomial part into the linear-like form and obtain the system

$\displaystyle\dot{x}$ $\displaystyle=AZ(x)+BW(x)u+R(x,u)$ (16)

where $R(x,u)=[R_{1}(x,u)\dots R_{n}(x,u)]^{\top}$,

$\displaystyle R_{i}(x,u)=R_{f_{i}}(x)+\sum^{m}_{j=1}R_{g_{ij}}(x)u_{j},$ (17)

$Z(x)$ is a vector of monomials in $x$ having degree $1$ to $r_{f}$, $W(x)$ is
a matrix of monomials in $z$ having degree $0$ to $r_{g}$, constant matrices
$A$, $B$ are unknown. Denote the degree of $u$ as $r_{u}$, then the remainder
$R_{i}(x,u)$ has the degree $r_{R}:=\max\\{r_{f},r_{g}+r_{u}\\}$.

####  3.2.1 Over-approximation of the feasible set

Similar to the case of first-order approximation, when Assumption 1 holds, one
can obtain an over-approximation of the feasible set by solving the
optimization problem (13). Abusing the notations $\overline{\textbf{A}}$,
$\overline{\textbf{B}}$, and $\overline{\textbf{C}}$, we define the over-
approximation set as

$\displaystyle\widetilde{\mathcal{I}}:=\left\\{\widehat{S}:\overline{\textbf{C}}+\widehat{S}\overline{\textbf{B}}+\overline{\textbf{B}}^{\top}\widehat{S}^{\top}+\widehat{S}\overline{\textbf{A}}\widehat{S}^{\top}\preceq
0\right\\},$ (18)

where
$\overline{\textbf{C}}=\overline{\textbf{B}}^{\top}\overline{\textbf{A}}^{-1}\overline{\textbf{B}}^{\top}-\delta
I$ for some arbitrarily fixed constant $\delta>0$, $\overline{\textbf{A}}$ and
$\overline{\textbf{B}}$ are the solutions to the optimization problem (13),
with

$\displaystyle\textbf{A}_{k}$
$\displaystyle=\ell(t_{k})\ell(t_{k})^{\top},~{}~{}\textbf{B}_{k}=-\ell(t_{k})\dot{x}(t_{k})^{\top},$
$\displaystyle\textbf{C}_{k}$
$\displaystyle=\dot{x}(t_{k})\dot{x}(t_{k})^{\top}-\gamma^{2}I,~{}~{}\ell(t_{k})=\begin{bmatrix}W\big{(}x(t_{k})\big{)}u(t_{k})\\\
Z\big{(}x(t_{k})\big{)}\end{bmatrix}.$

Similar to the analysis in Remark 5, if the data is rich enough, _i.e._ ,
$\begin{bmatrix}\overline{U}_{0}\\\ Z_{0}\end{bmatrix}$ has full row rank,
then the set $\mathcal{I}$ is bounded, allowing (13) to have a solution.

Note that in the linear like forms (16), the vector $Z(x)$ contains all
monomials in $x$ having degree from $1$ to $r_{f}$. Hence, the size of $Z(x)$
can be substantially large due to high truncating degree or high system order.
As pointed out in our previous work [11], a large vector $Z(x)$ tends to cause
computational issues in solving the sufficient conditions characterizing the
stabilizers. Therefore, to avoid bringing a large $Z(x)$ directly into the
conditions, we use a smaller vector $\widehat{Z}(x)$ for the Lyapunov function
and controller design, same as in [11]. Specifically, $\widehat{Z}(x)$ is
chosen as a $p\times 1$ vector on $x$ such that $\widehat{Z}(x)=0$ if and only
if $x=0$ and $Z(x)=H(x)\widehat{Z}(x)$ where the matrix $H(x)$ is non-unique.
Then, the linear-like system is rewritten as

$\displaystyle\dot{x}$ $\displaystyle=AH(x)\widehat{Z}(x)+BW(x)u+R(x,u).$

Define the Lyapunov function as
$V(x)=\widehat{Z}(x)^{\top}P^{-1}\widehat{Z}(x)$ where $P\succ 0$ and the
controller as $u=K(x)\widehat{Z}(x)$. Note that the vector $\widehat{Z}(x)$
needs to be chosen such that $V(x)$ is radially unbounded, and a simple choice
of $\widehat{Z}(x)$ is to make its first $n$ components coincide with $x$.

$\displaystyle\Upsilon(x)=-\frac{\partial\widehat{Z}(x)}{\partial
x}\left(-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)^{\top}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}-\\!\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\\!\left(-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)\frac{\partial\widehat{Z}(x)}{\partial
x}^{\top}\\!-\mu(x)\delta\frac{\partial\widehat{Z}(x)}{\partial
x}\frac{\partial\widehat{Z}(x)}{\partial x}^{\top}$ (19)
$\displaystyle\frac{\partial\widehat{Z}(x)}{\partial
x}\left(-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)^{\top}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}+(\star)^{\top}+\epsilon(x)I+\mu(x)\delta\frac{\partial\widehat{Z}(x)}{\partial
x}\frac{\partial\widehat{Z}(x)}{\partial
x}^{\top}\\!+\mu(x)^{-1}\\!\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\\!\overline{\textbf{A}}^{-1}\\!\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}\preceq 0$ (20)

####  3.2.2 Data-driven stabilizer design

Based-on the over-approximated set $\widetilde{\mathcal{I}}$, we can
characterize sufficient conditions for a local stabilizer using Lyapunov’s
second method and Petersen’s lemma. To make the condition tractable, the SOS
technique is applied to relax the positivity conditions for the polynomial
matrices. The design of the data-driven stabilizer is given in the following
result.

######  Theorem 3.

Under Assumption 1, if there exist polynomial $\mu(x)>0$ $\forall
x\in\mathbb{R}^{n}$, a positive definite polynomial $\epsilon(x)$, _i.e._ ,
zero at the origin and positive elsewhere, and matrices $Y(x)$ and $P$, such
that

$\displaystyle\begin{bmatrix}\Upsilon(x)-\epsilon(x)I_{p}&\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\\\ \begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}&\mu(x)\overline{\textbf{A}}\end{bmatrix}\text{ is SOS}$
(21)

where $\Upsilon(x)$ is defined in (19) with the parameter $\delta$ used to
define $\overline{\textbf{C}}$ for the set $\widetilde{\mathcal{I}}$, then the
state feedback controller $u=Y(x)P^{-1}\widehat{Z}(x)$ makes the origin a
locally asymptotically stable equilibrium for the system $\dot{x}=f(x)+g(x)u$.

###### Proof.

First, following the description (18) of set $\widetilde{\mathcal{I}}$ and the
definition of $\overline{\textbf{C}}$, it holds that, for all
$\widehat{S}\in\widetilde{\mathcal{I}}$,

$\displaystyle\quad~{}~{}\overline{\textbf{C}}+\widehat{S}\overline{\textbf{B}}+\overline{\textbf{B}}^{\top}\widehat{S}^{\top}+\widehat{S}\overline{\textbf{A}}\widehat{S}^{\top}$
$\displaystyle=\left(\widehat{S}^{\top}+\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)^{\top}\overline{\textbf{A}}\left(\widehat{S}^{\top}+\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)-\overline{\textbf{B}}^{\top}\overline{\textbf{A}}^{-1}\overline{\textbf{B}}+\overline{\textbf{C}}$
$\displaystyle=\left(\widehat{S}^{\top}+\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)^{\top}\overline{\textbf{A}}\left(\widehat{S}^{\top}+\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)-\delta
I\preceq 0.$

Define
$\Delta=\overline{\textbf{A}}^{1/2}\left(\widehat{S}^{\top}+\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)$,
and it follows that $\Delta^{\top}\Delta\preceq\delta I$ and

$\widehat{S}^{\top}=-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}+\overline{\textbf{A}}^{-1/2}\Delta.$

The objective is to find a control gain $K(x)=Y(x)P^{-1}$ that stabilizes the
approximated systems for all $\widehat{S}\in\widetilde{\mathcal{I}}$. Hence,
the closed-loop system of the controlled approximated system is

$\displaystyle\dot{x}$
$\displaystyle=\widehat{A}Z(x)+\widehat{B}W(x)Y(x)P^{-1}\widehat{Z}(x)$
$\displaystyle=\widehat{S}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}P^{-1}\widehat{Z}(x).$

Taking the time derivative of the Lyapunov function
$V(x)=\widehat{Z}(x)^{\top}P^{-1}\widehat{Z}(x)$ along the trajectory of the
closed-loop approximated system gives

$\displaystyle\dot{V}(x)=\widehat{Z}(x)^{\top}P^{-1}\frac{\partial\widehat{Z}(x)}{\partial
x}\widehat{S}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}P^{-1}\widehat{Z}(x)+(\star)^{\top}.$

Recall that
$\widehat{S}^{\top}=-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}+\overline{\textbf{A}}^{-1/2}\Delta$.
Let $\Phi(x)$ be such that
$\dot{V}(x)=\widehat{Z}(x)^{\top}P^{-1}\Phi(x)P^{-1}\widehat{Z}(x)$. Then, it
holds true that

$\displaystyle\Phi(x)$ $\displaystyle=\frac{\partial\widehat{Z}(x)}{\partial
x}\left(-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)^{\top}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}$ $\displaystyle\quad+\frac{\partial\widehat{Z}(x)}{\partial
x}\left(\overline{\textbf{A}}^{-1/2}\Delta\right)^{\top}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}+(\star)^{\top}.$

If the SOS condition in Theorem 3 is satisfied, using the Schur complement,
one derives (20). By Petersen’s Lemma (Lemma 2), the inequality (20) is
equivalent to $\Phi(x)\preceq-\epsilon(x)I_{p}$.

For the original systems with remainders, one has that

$\displaystyle\dot{V}(x)$
$\displaystyle\leq-\epsilon(x)\widehat{Z}(x)^{\top}P^{-2}\widehat{Z}(x)$
$\displaystyle\quad+2\widehat{Z}(x)^{\top}P^{-1}\frac{\partial\widehat{Z}(x)}{\partial
x}R\big{(}x,Y(x)P^{-1}\widehat{Z}(x)\big{)}.$

Note that the remainder $R\big{(}x,Y(x)P^{-1}\widehat{Z}(x)\big{)}$ having the
form (17) is of degree $r_{R}=\max\\{r_{f},r_{g}+r_{u}\\}$. As the degrees of
$\epsilon(x)$ and $\widehat{Z}(x)$ are fixed by design, we can guarantee that
the term $2\widehat{Z}(x)^{\top}P^{-1}\frac{\partial\widehat{Z}(x)}{\partial
x}R\big{(}x,Y(x)P^{-1}\widehat{Z}(x)\big{)}$ is of higher order and thus
converges to $0$ faster than
$-\epsilon(x)\widehat{Z}(x)^{\top}P^{-2}\widehat{Z}(x)$ for all $x$ in a
neighborhood of the origin. Moreover, as $\epsilon(x)>0$ $\forall x\neq 0$,
the origin is locally asymptotically stable for the closed-loop system
$\dot{x}=f(x)+g(x)Y(x)P^{-1}\widehat{Z}(x)$. ∎

######  Remark 6 (_Dimension of $\widehat{Z}(x)$_).

The feasibility of the condition $\Phi(x)\preceq-\epsilon(x)I_{p}$ when $x\neq
0$ implies that rank$(\Phi(x))=p$ for $x\neq 0$. Meanwhile, since
$\widehat{Z}(x)=0$ if and only if $x=0$, one has that $p\geq n$. Hence,
rank$(\frac{\partial\widehat{Z}(x)}{\partial x})\leq n$, and it holds that
rank$(\Phi(x))\leq 2n$ $\forall x\in\mathbb{R}^{n}$. Therefore, the
feasibility of $\Phi(x)\preceq-\epsilon(x)I_{p}$ implies that the dimension
$p$ of vector $\widehat{Z}(x)$ is such that $n\leq p\leq 2n$.
$\hfill\blacksquare$

$\displaystyle\begin{bmatrix}-P+\epsilon(x)I_{n}&\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\left(-\overline{\textbf{A}}^{-1}\overline{\textbf{B}}\right)\\\
\star&-P\end{bmatrix}+\mu(x)\delta\begin{bmatrix}0\\\
I_{n}\end{bmatrix}\begin{bmatrix}0\\\
I_{n}\end{bmatrix}^{\top}+\mu(x)^{-1}\begin{bmatrix}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\overline{\textbf{A}}^{-1/2}\\\
0\end{bmatrix}\begin{bmatrix}\begin{bmatrix}W(x)Y(x)\\\
H(x)P\end{bmatrix}^{\top}\overline{\textbf{A}}^{-1/2}\\\
0\end{bmatrix}^{\top}\preceq 0$ (22)

######  Remark 7 (_Comparison of Theorem 3 with previous results_).

Theorem 3 provides a local data-driven controller design for general nonlinear
dynamics. It is a continuous-time generalization of the previous result [8,
Theorem 6] where the first-order approximation and a linear controller is
considered. Compared to results that relies on specific choices of the basis
functions, such as [19], Theorem 3 synthesizes a data-driven controller using
Taylor polynomials as basis functions. The approach for the control of
polynomial systems in Theorem 3 is an alternative to the results presented in
[11] and [28], where [11] handles the additional noisy term in a different
way, and [28] searches for a Lyapunov function without restricting to a
special form. $\hfill\blacksquare$

##  4 Region of attraction estimation

In the previous section, we have shown that data-driven stabilizers can be
designed for nonlinear systems using first-order or high-order approximations.
The resulting controllers make the origin locally asymptotically stable.
Besides this property, it is of paramount importance to estimate the RoA of
the closed-loop system. The definition of the RoA is given as follows.

######  Definition 2 (_Region of attraction_).

For the systems $\dot{x}=f(x)$, if for every initial condition
$x(t_{0})\in\mathcal{R}$, it holds that $\lim_{t\rightarrow\infty}x(t)=0$,
then $\mathcal{R}$ is a region of attraction of the system with respect to the
origin. If there exists a $C^{1}$ function
$V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and a positive constant $c$ such that

$\displaystyle\Omega_{c}:=\\{x\in\mathbb{R}^{n}:V(x)\leq c\\}$

is bounded and

$\displaystyle V(0)=0,~{}~{}V(x)>0~{}~{}\forall x\in\mathbb{R}^{n}$
$\displaystyle\\{x\in\mathbb{R}^{n}:V(x)\leq c,x\neq
0\\}\subseteq\\{x\in\mathbb{R}^{n}:\dot{V}(x)<0\\},$

then $\Omega_{c}$ is an invariant subset, or called an estimation, of the RoA.

In this section, for the designed data-driven controllers in Section 3, we
derive data-driven conditions to determine whether a given sublevel set of the
Lyapunov function is an invariant subset of the RoA. The derived conditions
are data-driven because they are obtained using the over-approximated set
$\overline{\mathcal{I}}$. We note that once the controller is computed, it is
possible to use any other data-driven method to estimate the RoA, see for
example [23, 24, 25].

###  4.1 First-order approximation

By the controller design method in Section 3, the Lyapunov function $V(x)$,
and thus the set $\Omega_{c}$, are available for the analysis. To characterize
the set $\\{x\in\mathbb{R}^{n}:\dot{V}(x)<0\\}$, we need a bound on the
remainder for $x$ in a neighborhood of the origin. This is achievable by
posing the following assumption on the partial derivative of each $f_{i}$.

######  Assumption 2.

For all $z\in\mathbb{D}\subseteq\mathbb{R}^{n+m}$ where $\mathbb{D}$ is a
star-convex neighborhood of the origin, $f_{i}$ is continuously differentiable
and

$\displaystyle\left|\frac{\partial f_{i}}{\partial z_{j}}(z)-\frac{\partial
f_{i}}{\partial z_{j}}(0)\right|\leq L_{i}\|z\|~{}~{}\forall j=1,\dots,m+n,$
(23)

for $i=1,\dots,n$ with known $L_{i}>0$.

######  Remark 8 (_Existence of $L_{i}$_).

Assumption 2 is the weakest condition needed for deriving a bound on the
remainder using Lemma 1. A stronger condition, such as the Lipschitz
continuity of $\frac{\partial f_{i}}{\partial z_{j}}$, guarantees the
existence of $L_{i}$. It is also noted that $L_{i}$ can be estimated using a
data-based bisection procedure as shown in [13, Section III.C].
$\hfill\blacksquare$

Under Assumption 2, it follows from Lemma 1 that the first order approximation
remainder $R(x,u)$ of $f(x,u)$ satisfies that, for all $(x,u)\in\mathbb{D}$

$\displaystyle|R_{i}(x,u)|\leq\frac{\sqrt{m+n}L_{i}}{2}\|(x,u)\|^{2},~{}~{}i=1,\dots,n,$
(24)

where $R_{i}(x,u)$ is the $i$th element of vector $R(x,u)$.

######  Remark 9 (_On Assumptions 1 and 2_).

Under Assumption 2, using Lemma 1, the bound $\gamma^{2}$ in Assumption 1 can
be derived for the experimental data. During the experiment(s), suppose that
the smallest ball containing $(x(t),u(t))$ has radius $R_{e}$, _i.e._ ,
$\|(x(t_{k}),u(t_{k}))\|\leq R_{e}$ for all $k=0,\dots,T$. Then, for
$k=0,\dots,T$,

$\displaystyle\quad~{}R(x(t_{k}),u(t_{k}))^{\top}R(x(t_{k}),u(t_{k}))$
$\displaystyle=\sum^{n}_{i=1}R_{i}(x(t_{k}),u(t_{k}))^{2}$
$\displaystyle\leq\sum^{n}_{i=1}\frac{(m+n)L_{i}^{2}}{4}\|x(t_{k}),u(t_{k})\|^{4}$
$\displaystyle\leq\sum^{n}_{i=1}\frac{(m+n)L_{i}^{2}}{4}R_{e}^{4}.$

Hence, if some prior knowledge on the dynamics is known such that Assumption 2
holds, $\gamma^{2}$ can be chosen as
$\sum^{n}_{i=1}\frac{(m+n)L_{i}^{2}}{4}R_{e}^{4}$ to satisfy Assumption 1.
$\hfill\blacksquare$

After finding an estimate of the remainder $R(x,u)$, we analyze the RoA of the
closed-loop system.

Under Assumption 2, for the closed-loop system with the controller $u=Kx$
designed via Theorem 2, the derivative of the Lyapunov function can be
described in the following lemma.

######  Lemma 4.

Consider system (1) and the linear controller $u=YP^{-1}x$ where $Y$ and $P$
are designed to satisfy (14) with any given constant $w>0$. Under Assumption
2, the derivative of the Lyapunov function $V(x)=x^{\top}P^{-1}x$ satisfies
that, for all $x\in\mathbb{D}$,

$\displaystyle\dot{V}(x)\leq-wx^{\top}P^{-1}x+2\kappa(x)\rho(x)$ (25)

where

$\displaystyle\kappa(x):=\begin{bmatrix}x^{\top}Q_{1}\|(x,Kx)\|^{2}&\cdots&x^{\top}Q_{n}\|(x,Kx)\|^{2}\end{bmatrix},$
(26)

$Q_{i}(x)$ is the $i$th column of $P^{-1}$, and the vector $\rho(x)$ is
contained in the polytope

$\displaystyle\mathcal{H}:=\\{\varrho:-\bar{h}\preceq\varrho\preceq\bar{h}\\}$
(27)

with

$\bar{h}=\begin{bmatrix}\bar{h}_{1}~{}\cdots~{}\bar{h}_{n}\end{bmatrix}^{\top}=\begin{bmatrix}\frac{\sqrt{m+n}L_{1}}{2}&\cdots&\frac{\sqrt{m+n}L_{n}}{2}\end{bmatrix}^{\top}.$

The proof of Lemma 4 can be found in C.

Denote the number of distinct vertices of $\mathcal{H}$ as $\nu$ and each
vertex of $\mathcal{H}$ as $h_{k}$, $k=1,\dots,\nu$. Using the
Positivstellensatz result in Lemma 3, we present the following result.

######  Proposition 1.

Suppose that the controller $u=Kx$ renders the origin a locally asymptotically
stable equilibrium for (1) with the Lyapunov function $V(x)=x^{\top}P^{-1}x$.
Under Assumption 2, given a $c>0$ such that
$\Omega_{c}=\\{x\in\mathbb{R}^{n}:V(x)\leq c\\}\subseteq\mathbb{D}$, if there
exist SOS polynomials $s_{1k},s_{2k}$ in $x$, $k=1,\dots,\nu$ such that

$\displaystyle-\left[s_{1k}(c-V(x))+s_{2k}\left(-wx^{\top}P^{-1}x+2\kappa(x)h_{k}\right)+x^{\top}x\right]$
(28)

is SOS, where $\kappa(x)$ is defined as in (26) and $h_{k}$ are the distinct
vertices of the polytope $\mathcal{H}$ defined in (27), then $\Omega_{c}$ is
an invariant subset of the RoA of the system $\dot{x}=f(x,Kx)$ relative to the
equilibrium $x=0$.

###### Proof.

According to [39, Page 87], the polytope $\mathcal{H}$ can be expressed as

$\displaystyle\mathcal{H}=\left\\{\varrho=\sum^{\nu}_{k=1}\lambda_{k}(x)h_{k},~{}\sum^{\nu}_{k=1}\lambda_{k}(x)=1,~{}\lambda_{k}(x)\geq
0\right\\}$

for any fixed $x\in\Omega_{c}$. Then, the derivative of the Lyapunov function
satisfies that

$\displaystyle\dot{V}(x)$ $\displaystyle\leq-
wx^{\top}P^{-1}x+2\kappa(x)\rho(x)$
$\displaystyle=-wx^{\top}P^{-1}x+2\kappa(x)\sum^{\nu}_{k=1}\lambda_{k}(x)h_{k}$
$\displaystyle=\sum^{\nu}_{k=1}\lambda_{k}(x)(-wx^{\top}P^{-1}x)+\sum^{\nu}_{k=1}\lambda_{k}(x)2\kappa(x)h_{k}$
$\displaystyle=\sum^{\nu}_{k=1}\lambda_{k}(x)\left(-wx^{\top}P^{-1}x+2\kappa(x)h_{k}\right)$

As $\lambda_{k}(x)\geq 0$ and $\sum^{\nu}_{k=1}\lambda_{k}(x)=1$, if

$-wx^{\top}P^{-1}x+2\kappa(x)h_{k}<0$

holds for all $k=1,\dots,\nu$, then $\dot{V}(x)<0$.

By Lemma 3, for each $k=1,\dots,\nu$, if there exist SOS polynomials
$s_{1k},s_{2k}$ such that

$\displaystyle-\left[s_{1k}(c-V(x))+s_{2k}\left(-wx^{\top}P^{-1}x+2\kappa(x)h_{k}\right)+x^{\top}x\right]$

is SOS, then the set inclusion condition

$\displaystyle\\{x\in\mathbb{R}^{n}:$ $\displaystyle V(x)\leq c,~{}x\neq 0\\}$
$\displaystyle\subseteq\\{x\in\mathbb{R}^{n}:-wx^{\top}P^{-1}x+2\kappa(x)h_{k}<0\\}$

holds. This leads to the set inclusion condition

$\displaystyle\\{x\in\mathbb{R}^{n}:V(x)\leq c,~{}x\neq
0\\}\subseteq\\{x\in\mathbb{R}^{n}:\dot{V}(x)<0\\},$

and hence $\Omega_{c}$ is an inner estimate of the ROA. ∎

######  Remark 10 (_Numerical method for RoA estimation_).

In [19], the RoA is estimated by a numerical method, _i.e._ , a sufficient
condition of $V(x^{+})-V(x)<0$ is found using data, and the grids in a compact
region are tested to see whether the sufficient condition is satisfied so that
the sublevel sets of $V(x)$ can be found as the RoA estimation. In this work,
we derive SOS conditions for the RoA estimation that is an alternative to the
mesh method used in [19].

###  4.2 High-order approximation

For high-order approximation, to analyze the closed-loop system, we also need
a bound on the high-order remainder $R_{f}(x)+R_{g}(x)u$. To simplify the
analysis, this subsection considers system (15) with a single input, that is
$u\in\mathbb{R}$. In this case, the remainder takes the form

$R(x,u)=R_{f}(x)+R_{g}(x)u$

where

$\displaystyle R(x,u)$
$\displaystyle=[R_{1}(x,u)~{}\cdots~{}R_{n}(x,u)]^{\top},$ $\displaystyle
R_{f}(x)$ $\displaystyle=[R_{f_{1}}(x)~{}\cdots~{}R_{f_{n}}(x)]^{\top},\text{
and}$ $\displaystyle R_{g}(x)$
$\displaystyle=[R_{g_{1}}(x)~{}\cdots~{}R_{g_{n}}(x)]^{\top}.$

To obtain a bound on the remainder using Lemma 1, we pose an assumption on the
high-order partial derivatives of the functions $f(x)$ and $g(x)$.

######  Assumption 3.

For all $x\in\mathbb{D}\subseteq\mathbb{R}^{n}$ where $\mathbb{D}$ is a star-
convex neighborhood of the origin, $\partial^{\alpha}f_{i}$,
$\partial^{\beta}g_{i}$, $|\alpha|=r_{f}$, $|\beta|=r_{g}$, are absolutely
continuous for $x\in\mathbb{D}$ and

$\displaystyle|\partial^{\alpha}f_{i}(x)-\partial^{\alpha}f_{i}(0)|$
$\displaystyle\leq L_{i}\|x\|,$
$\displaystyle|\partial^{\beta}g_{i}(x)-\partial^{\beta}g_{i}(0)|$
$\displaystyle\leq M_{i}\|x\|~{}~{}i=1,\dots,n,$

with known $L_{i},M_{i}>0$.

Under Assumption 3, by Lemma 1, it holds that for all $x\in\mathbb{D}$

$\displaystyle|R_{f_{i}}(x)|$
$\displaystyle\leq\frac{\sqrt{n}L_{i}}{(r_{f}+1)!}\|x\|^{r_{f}+1},$ (29)
$\displaystyle|R_{g_{i}}(x)|$
$\displaystyle\leq\frac{\sqrt{n}M_{i}}{(r_{g}+1)!}\|x\|^{r_{g}+1}.$ (30)

Denote the degree of the control input $u$ as $r_{u}$ and note that there is
no constant terms in $u$. Hence, it holds that

$\displaystyle|u|\leq\sum^{r_{u}}_{j=1}\overline{K}_{j}\|x\|^{j}$

for some positive constants $\overline{K}_{j}$ obtained from the designed
$u=K(x)\widehat{Z}(x)$. The term $R_{g_{i}}(x)u$ can be bounded by

$\displaystyle|R_{g_{i}}(x)u|$ $\displaystyle\leq|R_{g_{i}}(x)||u|$
$\displaystyle\leq\frac{\sqrt{n}M_{i}}{(r_{g}+1)!}\|x\|^{r_{g}+1}\cdot\sum^{r_{u}}_{j=1}\overline{K}_{j}\|x\|^{j}$
$\displaystyle=\frac{\sqrt{n}M_{i}}{(r_{g}+1)!}\sum^{r_{u}}_{j=1}\overline{K}_{j}\|x\|^{r_{g}+1+j}.$

For any $r_{g}\geq 1$ and $r_{u}\geq 1$, one can always find a number $q\geq
1$ and positive constants $\widetilde{K}_{j}$ such that

$\displaystyle\sum^{r_{u}}_{j=1}\overline{K}_{j}\|x\|^{r_{g}+1+j}\leq\sum^{q}_{j=1}\widetilde{K}_{j}\|x\|^{2j}\leq\widetilde{K}_{M}\sum^{q}_{j=1}\|x\|^{2j}$
(31)

where $\widetilde{K}_{M}=\max_{j}{\widetilde{K}_{j}}$. (Note: for any odd
degree term $\|x\|^{a+b}$, we can use the triangular inequality to bound it as
$\|x\|^{a+b}\leq\epsilon\|x\|^{2a}+\epsilon^{-1}\|x\|^{2b}$
$\forall\epsilon>0$.) We then bound the remainder
$R_{i}(x,K(x)\widehat{Z}(x))$ as

$\displaystyle\quad~{}|R_{i}(x,K(x)\widehat{Z}(x))|$
$\displaystyle\leq|R_{f_{i}}(x)|+|R_{g_{i}}(x)u|$
$\displaystyle\leq\frac{\sqrt{n}L_{i}}{(r_{f}+1)!}\|x\|^{r_{f}+1}+\frac{\sqrt{n}M_{i}\widetilde{K}_{M}}{(r_{g}+1)!}\sum^{q}_{j=1}\|x\|^{2j}.$

To make sure the obtained bound is a polynomial, we need to set $r_{f}$ as an
odd number. Letting

$\widetilde{\varrho}_{i}=\max\left\\{\frac{\sqrt{n}L_{i}}{(r_{f}+1)!},\frac{\sqrt{n}M_{i}\widetilde{K}_{M}}{(r_{g}+1)!}\right\\},$

one can write

$R_{i}(x,K(x)\widehat{Z}(x))=\widetilde{\rho}_{i}(x)\left(\|x\|^{r_{f}+1}+\sum^{q}_{j=1}\|x\|^{2j}\right)$

where $|\widetilde{\rho}_{i}(x)|\leq\widetilde{\varrho}_{i}$ $\forall
x\in\mathbb{D}$.

Denote $Q_{i}(x)$ as the $i$th column of
$P^{-1}\frac{\partial\widehat{Z}(x)}{\partial x}$, and one has the following
term in the derivative of the Lyapunov function

$\displaystyle\quad~{}2\widehat{Z}(x)^{\top}P^{-1}\frac{\partial\widehat{Z}(x)}{\partial
x}R\big{(}x,K(x)\widehat{Z}(x)\big{)}$
$\displaystyle=2\begin{bmatrix}\widehat{Z}(x)^{\top}Q_{1}(x)&\cdots&\widehat{Z}(x)^{\top}Q_{n}(x)\end{bmatrix}R\big{(}x,K(x)\widehat{Z}(x)\big{)}$
$\displaystyle=2\sum^{n}_{i=1}\widehat{Z}(x)^{\top}Q_{i}(x)R_{i}\big{(}x,K(x)\widehat{Z}(x)\big{)}.$

Defining

$\displaystyle\widetilde{\kappa}(x)=\Big{[}$
$\displaystyle\widehat{Z}(x)^{\top}Q_{1}(x)(\|x\|^{r_{f}+1}+\sum^{q}_{j=1}\|x\|^{2j})~{}\cdots$
$\displaystyle\quad\widehat{Z}(x)^{\top}Q_{n}(x)(\|x\|^{r_{f}+1}+\sum^{q}_{j=1}\|x\|^{2j})\Big{]}$
(32)

and
$\widetilde{\rho}(x)=\begin{bmatrix}\widetilde{\rho}_{1}(x)&\cdots&\widetilde{\rho}_{n}(x)\end{bmatrix}^{\top}$
gives

$\displaystyle
2\widehat{Z}(x)^{\top}P^{-1}\frac{\partial\widehat{Z}(x)}{\partial
x}R\big{(}x,Y(x)P^{-1}\widehat{Z}(x)\big{)}=2\widetilde{\kappa}(x)\widetilde{\rho}(x),$

and thus

$\displaystyle\dot{V}(x)\leq-\epsilon(x)\widehat{Z}(x)^{\top}P^{-2}\widehat{Z}(x)+2\widetilde{\kappa}(x)\widetilde{\rho}(x)$
(33)

for all $x\in\mathbb{D}$, where the vector $\widetilde{\rho}(x)$ is contained
in the polytope

$\displaystyle\widetilde{\mathcal{H}}:=\left[-\widetilde{\varrho}_{1},\widetilde{\varrho}_{1}\right]\times\cdots\times\left[-\widetilde{\varrho}_{n},\widetilde{\varrho}_{n}\right].$
(34)

Following the same spirit of Proposition 1, we obtain the following result.

######  Proposition 2.

Suppose that the $u=Y(x)P^{-1}\widehat{Z}(x)$ renders the origin a locally
asymptotically stable equilibrium for (15) with the Lyapunov function
$V(x)=\widehat{Z}(x)^{\top}P^{-1}\widehat{Z}(x)$. Under Assumption 2, given
$c>0$ such that $\Omega_{c}=\\{x\in\mathbb{R}^{n}:V(x)\leq
c\\}\subseteq\mathbb{D}$, if there exist SOS polynomials $s_{1k},s_{2k}$ in
$x$, $k=1,\dots,\nu$ such that,

$\displaystyle-\left[s_{1k}(c-V(x))+s_{2k}(-\epsilon(x)\widehat{Z}(x)^{\top}P^{-2}\widehat{Z}(x)\right.$
$\displaystyle\left.+2\widetilde{\kappa}(x)\widetilde{h}_{k})+x^{\top}x\right]\text{
is SOS}$ (35)

where $\widetilde{\kappa}(x)$ is defined as in (4.2) and $\widetilde{h}_{k}$
are the distinct vertices of polytope $\widetilde{\mathcal{H}}$ defined in
(34), then $\Omega_{c}$ is an invariant subset of the RoA of the system
$\dot{x}=f(x)+g(x)Y(x)P^{-1}\widehat{Z}(x)$ relative to the equilibrium $x=0$.

We summarize the data-driven approach for controller design and RoA estimation
in Algorithm 1.

Algorithm 1 Data-driven controller design and RoA estimation

Collect data $X_{0}$, $X_{1}$, and $U_{0}$ from the offline open-loop
experiment(s)

Solve the optimization problem (13) for $\overline{\textbf{A}}$,
$\overline{\textbf{B}}$, and $\overline{\textbf{C}}$ to obtain
$\overline{\mathcal{I}}$

Solve the LMI (14) in Theorem 2 ( the SOS condition (21) in Theorem 3) for $Y$
($Y(x)$) and $P$ to obtain the controller $u=YP^{-1}x$
($u=Y(x)P^{-1}\widehat{Z}(x)$)

Bound $\dot{V}(x)$ as in (25) ((33)) for the closed-loop system using the
designed $u$ and $V(x)$

Determine if $\Omega_{c}=\\{x\in\mathbb{R}^{n}:V(x)\leq c\\}$ with a given
$c>0$ belongs to the RoA by checking the SOS condition (28) ((2)) for the
first-order (high-order) approximation

##  5 Numerical example

In this section, the proposed control design and closed-loop system analysis
is applied to the nonlinear benchmark, an inverted pendulum. Simulation
results are presented to illustrate the main result of this work.

###  5.1 First-order approximation

Consider the inverted pendulum having the dynamics

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2},$ $\displaystyle\dot{x}_{2}$
$\displaystyle=\frac{mgl}{J}\mathrm{sin}(x_{1})-\frac{r}{J}x_{2}+\frac{l}{J}\mathrm{cos}(x_{1})u.$
(36)

The system is in the form of a general nonlinear system with
$f(x,u)=[f_{1}(x,u)~{}~{}f_{2}(x,u)]^{\top}$ where

$f_{1}(x,u)=x_{2},~{}~{}f_{2}(x,u)=\frac{mgl}{J}\mathrm{sin}(x_{1})-\frac{r}{J}x_{2}+\frac{l}{J}\mathrm{cos}(x_{1})u.$

To collect data, an experiment is conducted with
$x(0)=[0.01~{}~{}-0.01]^{\top}$ and $u=0.1\mathrm{sin}(t)$ during the time
interval $[0,5]$. The data is sampled with fixed sampling period
$T_{s}=0.05$s. The system parameters are $m=0.1$, $g=9.8$, $r=1$, $l=1$, and
$J=1$. Collect the data and arrange it into data sets with length $T=10$. We
assume that the remainder is over-approximated by $100\%$; in other words, the
bound $\gamma$ is twice the largest instantaneous norm of the remainder during
the experiment. Then, for the experimental data, Assumption 1 holds with

$\gamma=3.3352\cdot 10^{-6}.$

Setting $\delta=0.01$, we first solve the optimization problem (13) to find
$\overline{\mathcal{I}}$, and then apply Theorem 2 with $w=1$. The solution
found by CVX is

$\displaystyle P$ $\displaystyle=10^{3}\cdot\begin{bmatrix}1.0152&-1.3289\\\
-1.3289&1.7727\end{bmatrix},$ $\displaystyle u$
$\displaystyle=-12.0432x_{1}-8.887x_{2}.$ (37)

The value of the constant is obtained as follows. We focus on $f_{2}$ whose
first order partial derivatives are

$\displaystyle\frac{\partial f_{2}}{\partial x_{1}}$
$\displaystyle=\frac{mgl}{J}\mathrm{cos}(x_{1})-\frac{l}{J}\mathrm{sin}(x_{1})u,$
$\displaystyle\frac{\partial f_{2}}{\partial x_{2}}$
$\displaystyle=-\frac{r}{J},~{}\frac{\partial f_{2}}{\partial
u}=\frac{l}{J}\mathrm{cos}(x_{1}).$

Assumption 2 holds for $f_{2}$ with

$L_{2}=\sqrt{2}~{}\max\left(\frac{mgl}{J},\frac{l}{J}\right).$

Assume that $L_{2}$ is over-estimated by $20\%$, which gives the estimated
bound on the remainder $R_{2}(x,u)$ as

$|R_{2}(x,u)|\leq
1.2\cdot\frac{\sqrt{3}L_{2}}{2}\|(x,u)\|^{2}=1.4697\|(x,u)\|^{2}.$

Applying Proposition 1, the largest $c$ found for the controller (5.1) is
$c^{*}=7.58\cdot 10^{-4}$.

To see the effect brought by the parameter $\delta$ in the over-approximation
of the feasible set, the same design and analysis approach is repeated for
different values of $\delta$. The simulation results associated with various
$\delta$ are recorded in Table 1. The resulting sets $\Omega_{c}$ are
illustrated in Figure 1. As observed from the simulation results, the control
gains $K$ do not show significant changes with different values of $\delta$,
while the sizes of the estimated RoA vary. The role of $\delta$ in the RoA
estimation remains unclear at the moment and interesting questions arise from
the simulation results. For instance, how $\delta$ affects the RoA estimation
and how to find an optimal $\delta$ for the RoA estimation. It is of interest
to answer these questions via a careful and thorough study in the future.

$\delta$ | $P$ | $K$ | $c^{*}$
---|---|---|---
$1$ | $10^{3}\cdot\\!\begin{bmatrix}3.13&-4.08\\\ -4.08&5.43\end{bmatrix}$ | $[-11.84~{}-8.75]$ | $0.90\cdot\\!10^{-6}$
$10^{-2}$ | $10^{3}\cdot\\!\begin{bmatrix}1.02&-1.33\\\ -1.33&1.77\end{bmatrix}$ | $[-12.04~{}-8.89]$ | $7.58\cdot\\!10^{-4}$
$10^{-4}$ | $10^{3}\cdot\\!\begin{bmatrix}0.427&-0.55\\\ -0.55&0.74\end{bmatrix}$ | $[-11.67~{}-8.63]$ | $7.54\cdot\\!10^{-4}$
Table 1: Simulation results on the continuous-time model with different values
of $\delta$.
Figure 1: The sets $\Omega_{c}$ found by Proposition 1 with different values
of $\delta$.

###  5.2 High-order approximation

Consider again the inverted pendulum (5.1). Setting degrees of the Taylor
polynomials as $r_{f}=5$ and $r_{g}=2$, one can write the dynamics as

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2},$ $\displaystyle\dot{x}_{2}$
$\displaystyle=\frac{mgl}{J}x_{1}-\frac{mgl}{6J}x_{1}^{3}+\frac{mgl}{120J}x_{1}^{5}-\frac{r}{J}x_{2}$
$\displaystyle\quad+\left(\frac{l}{J}-\frac{l}{2J}x_{1}^{2}\right)u+R_{f_{2}}(x)+R_{g_{2}}(x)u$

where the remainders are

$\displaystyle R_{f_{2}}(x)$
$\displaystyle=\left.\frac{mgl}{J}\cdot\sin(\sigma
x_{1})\right|_{\sigma\in(0,1)}\cdot\frac{x_{1}^{6}}{6!},$ $\displaystyle
R_{g_{2}}(x)$ $\displaystyle=\left.\frac{l}{J}\cdot\cos(\sigma
x_{1})\right|_{\sigma\in(0,1)}\cdot\frac{x_{1}^{3}}{3!}.$

We again assume that the remainder data is over-approximated by $100\%$, and
hence Assumption 1 holds with

$\gamma=2.1602\cdot 10^{-4}.$

Set $\delta=1$ and the degree of the controller $r_{u}$ as $3$. The data-
driven controller designed by Theorem 3 is

$\displaystyle u$ $\displaystyle=x_{1}(1.5x_{1}^{2}-0.098x_{2}^{2}-11.4)$
$\displaystyle\quad+x_{2}(0.35x_{1}^{2}-0.036x_{2}^{2}-2.0).$

For all $x_{1},x_{2}\in\mathbb{R}$, one can bound the control input by

$|u|\leq 2.8\sqrt{2}\|x\|^{3}+11\sqrt{2}\|x\|.$

On the other hand, for all $x_{1},x_{2}\in\mathbb{R}$, it holds that

$|R_{f_{2}}(x)|\leq\frac{mgl}{720J}\|x\|^{6}\text{ and
}|R_{g_{2}}(x)|\leq\frac{l}{6J}|x_{1}^{3}|\leq\frac{l}{6J}\|x\|^{3}.$

Then, the remainder is bounded as

$\displaystyle|R_{f_{2}}(x)+R_{g_{2}}(x)u|$
$\displaystyle\leq|R_{f_{2}}(x)|+|R_{g_{2}}(x)||u|$
$\displaystyle=\\!\left(\frac{mgl}{720J}+\frac{7\sqrt{2}l}{15J}\right)\\!\|x\|^{6}\\!+\frac{11\sqrt{2}l}{6J}\|x\|^{4}$
$\displaystyle\leq 2.5927(\|x\|^{6}+\|x\|^{4}).$

For the RoA estimation, we suppose that the remainder bound is over-estimated
by $20\%$. More specifically, for any $x_{1},x_{2}\in\mathbb{R}$, there exists
a polynomial $\widetilde{\rho}_{2}(x)$ such that

$R_{f_{2}}(x)+R_{g_{2}}(x)u=\widetilde{\rho}_{2}(x)(\|x\|^{6}+\|x\|^{4})$

where $|\widetilde{\rho}_{2}(x)|\leq 1.2\cdot 2.5927=3.1112$. Let $Q_{2}$ be
the second column of $P^{-1}$ and we obtain that
$\widetilde{\kappa}(x)=x^{\top}Q_{2}(\|x\|^{6}+\|x\|^{4})$.

Using this bound and Proposition 2, we find an estimation of the RoA that is
$\\{x\in\mathbb{R}^{2}:x^{\top}P^{-1}x\leq 1.58\\}$ which is illustrated as
the darkest area in Figure 2. In the same figure, the light grey area is the
estimated RoA by checking point-by-point of a mesh of initial conditions using
explicit dynamics, and the medium dark area is the largest sublevel set of the
Lyapunov function contained in the RoA.

Figure 2: Estimations of the RoA using the Lyapunov function by the high-order
approximation ($\delta=1$). The lightest grey area is the estimated RoA via
the numerical method using explicit dynamics; the medium grey area is the
largest sublevel set of the Lyapunov function contained in the numerically
estimated RoA; the dark grey area is the estimated RoA obtained by Proposition
2.
Figure 3: Estimations of the RoA using the Lyapunov function by the first-
order approximation ($\delta=10^{-2}$). The lightest grey area is the
estimated RoA using the numerical method with explicit dynamics; the medium
grey area is the largest sublevel set of the Lyapunov function contained in
the numerically estimated RoA; the dark grey area is the estimated RoA
obtained by Proposition 1.
Figure 4: The sets $\Omega_{c}$ found by Proposition 2 with different values of $\delta$. $\delta$ | $P$ | $u$ | $c^{*}$
---|---|---|---
$0.5$ | $10^{-6}\cdot\begin{bmatrix}0.45&-0.63\\\ -0.63&4.57\end{bmatrix}$ | $x_{1}(0.64x_{1}^{2}-0.039x_{2}^{2}-5.15)+x_{2}(0.22x_{1}^{2}-0.027x_{2}^{2}-1.02)$ | $0.999$
$1$ | $10^{-4}\cdot\begin{bmatrix}0.47&-1.28\\\ -1.28&8.71\end{bmatrix}$ | $x_{1}(1.5x_{1}^{2}-0.098x_{2}^{2}-11.4)+x_{2}(0.35x_{1}^{2}-0.036x_{2}^{2}-2.0)$ | $1.58$
$2$ | $10^{-4}\cdot\begin{bmatrix}0.92&-4.18\\\ -4.18&25.9\end{bmatrix}$ | $x_{1}(5.09x_{1}^{2}-0.41x_{2}^{2}-36.96)+x_{2}(0.97x_{1}^{2}-0.090x_{2}^{2}-6.33)$ | $0.242$
Table 2: Simulation results on the high-order approximation with different
values of $\delta$.

In comparison, for the first-order approximation, the estimated RoA via the
numerical method, the largest sublevel set of the Lyapunov function contained
in the RoA, and the RoA estimated by Proposition 1 are illustrated in Figure
3.

Comparing the largest RoA estimation with known dynamics and the numerical
method, we see that the nonlinear controller leads to a larger RoA than the
linear controller. However, if the RoA is estimated using only data, the
nonlinear controller gives a much smaller area than the linear one. This could
be caused by the conservativeness introduced in the process of bounding the
remainder $R(x,u)$ for the high-order approximation.

For the high-order approximation, the data-driven controller design and RoA
estimation are also analyzed for different values of $\delta$, and the
simulation results are present in Table 2 and Figure 4. In this case, we could
not test values of $\delta$ in a larger range, as the SOS condition in Theorem
3 is solvable by the SOSTOOLS only with $\delta$ in a small range. As shown by
the simulation results, the estimations of the RoA vary with different values
of $\delta$, and similar to the previous cases, how $\delta$ affects the
estimation still remains as an open problem.

It should be pointed out that, besides $\delta$, there are other parameters
that may affect the solving of the SOS conditions in Propositions 1 and 2, for
example, the degree of the SOS polynomials, which could also affect the RoA
estimation. At the moment, it remains to be understood what is the optimal
choices of the parameters that will lead to the largest estimation of the RoA.
Hence, we expect that the RoA estimation for the high-order approximation case
can be enlarged by carefully choosing the parameters. In this work, we aim at
showing that one can use data to estimate the RoA of a nonlinear systems under
a data-driven controller. Later works may look into how to enlarge the RoA
estimation in the data-based setting.

##  6 Conclusion

For general nonlinear dynamics without explicit information on the
nonlinearities, this paper proposes data-driven stabilizer designs and RoA
analysis by approximating the unknown functions using Taylor’s expansion.
Using finite-length input-state data, linear and nonlinear stabilizers are
designed for continuous-time nonlinear systems that render the known
equilibrium locally asymptotically stable. Then, by estimating a bound on the
Taylor remainder, data-driven conditions are given to find an invariant subset
of the RoA. Simulation results on the inverted pendulum show the designed
data-driven controllers and the RoA estimations for both first-order and high-
order approximations. The estimation of the RoA can be conservative especially
for the high-order approximation, and may be further enlarged by optimizing
the choice of some parameters in the design and estimation steps. Topics such
as enlarging the RoA estimation and case studies on more complicated nonlinear
benchmarks are all interesting directions to be considered in future works.

## References

  * [1] Z. Hou, Z. Wang, From model-based control to data-driven control: Survey, classification and perspective, Information Sciences 235 (2013) 3–35.
  * [2] M. Tanaskovic, L. Fagiano, C. Novara, M. Morari, Data-driven control of nonlinear systems: An on-line direct approach, Automatica 75 (2017) 1–10.
  * [3] Z. Hou, S. Xiong, On model-free adaptive control and its stability analysis, IEEE Transactions on Automatic Control 64 (11) (2019) 4555–4569.
  * [4] P. Tabuada, W. L. Ma, J. Grizzle, A. D. Ames, Data-driven control for feedback linearizable single-input systems, in: Proceedings of the 56th IEEE Conference on Decision and Control, Melbourne, VIC, Australia, 2017, pp. 6265–6270.
  * [5] L. Fraile, M. Marchi, P. Tabuada, Data-driven stabilization of SISO feedback linearizable systems, arXiv preprint arXiv:2003.14240 (2021).
  * [6] J. Berberich, J. Köhler, M. A. Müller, F. Allgöwer, Linear tracking MPC for nonlinear systems - Part II: The data-driven case, IEEE Transactions on Automatic Control 67 (9) (2022) 4406 – 4421.
  * [7] W. Liu, J. Sun, G. Wang, F. Bullo, J. Chen, Data-driven self-triggered control via trajectory prediction, arXiv preprint arXiv:2207.08596 (2022).
  * [8] C. De Persis, P. Tesi, Formulas for data-driven control: Stabilization, optimality, and robustness, IEEE Transactions on Automatic Control 65 (3) (2020) 909–924.
  * [9] H. J. van Waarde, M. K. Camlibel, A matrix Finsler’s lemma with applications to data-driven control, in: Proceedings of the 60th IEEE Conference on Decision and Control, Austin, TX, USA, 2021, pp. 5777–5782.
  * [10] T. Dai, M. Sznaier, Nonlinear data-driven control via state-dependent representations, in: Proceedings of the 60th IEEE Conference on Decision and Control, Austin, TX, USA, 2021, pp. 5765–5770.
  * [11] M. Guo, C. De Persis, P. Tesi, Data-driven stabilization of nonlinear polynomial systems with noisy data, IEEE Transactions on Automatic Control 67 (8) (2022) 4210 – 4217.
  * [12] A. Nejati, B. Zhong, M. Caccamo, M. Zamani, Data-driven controller synthesis of unknown nonlinear polynomial systems via control barrier certificates, Proceedings of The 4th Annual Learning for Dynamics and Control Conference 168 (2022) 763–776.
  * [13] T. Martin, F. Allgöwer, Data-driven system analysis of nonlinear systems using polynomial approximation, arXiv preprint arXiv:2108.11298 (2022).
  * [14] M. Guo, C. De Persis, P. Tesi, Learning control for polynomial systems using sum of squares relaxations, in: Proceedings of the 59th IEEE Conference on Decision and Control, Jeju Island, Republic of Korea, 2020, pp. 2436–2441.
  * [15] A. A. Ahmadi, B. E. Khadir, Learning dynamical systems with side information, Proceedings of the 2nd Conference on Learning for Dynamics and Control 120 (2020) 718–727.
  * [16] A. Luppi, C. De Persis, P. Tesi, On data-driven stabilization of systems with nonlinearities satisfying quadratic constraints, Systems & Control Letters 163 (2022) 105206(1)–105206(11).
  * [17] R. Strasser, J. Berberich, F. Allgöwer, Data-driven control of nonlinear systems: Beyond polynomial dynamics, in: Proceedings of the 60th IEEE Conference on Decision and Control, Austin, TX, USA, 2021, pp. 4344–4351.
  * [18] E. Kaiser, J. N. Kutz, S. L. Brunton, Data-driven discovery of Koopman eigenfunctions for control, Machine Learning: Science and Technology 2 (3) (2021) 035023(1)–035023(31).
  * [19] C. De Persis, M. Rotulo, P. Tesi, Learning controllers from data via approximate nonlinearity cancellation, IEEE Transactions on Automatic Control, doi: 10.1109/TAC.2023.3234889 (2023).
  * [20] G. Chesi, Estimating the domain of attraction for non-polynomial systems via LMI optimizations, Automatica 45 (6) (2009) 1536–1541.
  * [21] W. Tan, A. Packard, Stability region analysis using polynomial and composite polynomial Lyapunov functions and sum-of-squares programming, IEEE Transactions on Automatic Control 53 (2) (2008) 565–571.
  * [22] U. Topcu, A. Packard, Local stability analysis for uncertain nonlinear systems, IEEE Transactions on Automatic Control 54 (5) (2009) 1042–1047.
  * [23] R. Bobiti, M. Lazar, Automated sampling–based stability verification and DOA estimation for nonlinear systems, IEEE Transactions on Automatic Control 63 (11) (2018) 3659–3674.
  * [24] B. K. Colbert, M. M. Peet, Using trajectory measurements to estimate the region of attraction of nonlinear systems, in: Proceedings of the 57th IEEE Conference on Decision and Control, Miami, FL, USA, 2018, pp. 2341–2347.
  * [25] Y. Shen, M. Bichuch, E. Mallada, Model-free learning of regions of attraction via recurrent sets, in: Proceedings of the 61st IEEE Conference on Decision and Control, Cancún, Mexico, 2022, pp. 4714–4719.
  * [26] I. R. Petersen, A stabilization algorithm for a class of uncertain linear systems, Systems & Control Letters 8 (4) (1987) 351–357.
  * [27] A. Papachristodoulou, S. Prajna, A tutorial on sum of squares techniques for systems analysis, in: Proceedings of 2005 American Control Conference, Portland, OR, USA, 2005, pp. 2686–2700.
  * [28] A. Bisoffi, C. De Persis, P. Tesi, Data-driven control via Petersen’s lemma, Automatica 145 (2022) 110537(1)–110537(14).
  * [29] G. Stengle, A nullstellensatz and a positivstellensatz in semialgebraic geometry, Mathematische Annalen 207 (1974) 87–97.
  * [30] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268.
  * [31] R. Chartrand, Numerical differentiation of noisy, nonsmooth data, International Scholarly Research Notices 2011 (2011).
  * [32] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113 (15) (2016) 3932–3937.
  * [33] G. B. Folland, Higher-order derivatives and Taylor’s formula in several variables, https://sites.math.washington.edu/ folland/Math425/taylor2.pdf.
  * [34] G. B. Folland, Remainder estimates in Taylor’s theorem, The American Mathematical Monthly 97 (3) (1990) 233–235.
  * [35] P. Shcherbakov, M. Topunov, Extensions of Petersen’s lemma on matrix uncertainty, IFAC Proceedings Volumes 41 (2) (2008) 11385–11390.
  * [36] G. Chesi, LMI techniques for optimization over polynomials in control: A survey, IEEE Transactions on Automatic Control 55 (11) (2010) 2500–2510.
  * [37] A. Bisoffi, C. De Persis, P. Tesi, Trade-offs in learning controllers from noisy data, Systems & Control Letters 154 (2021) 104985(1) – 104985(8).
  * [38] S. Boyd, L. El Ghaoui, E. Feron, V. Balakrishnan, Linear matrix inequalities in system and control theory, Vol. 15 of SIAM Studies in Applied Mathematics, SIAM, 1994.
  * [39] F. Blanchini, S. Miani, Set-Theoretic Methods in Control, Systems & Control: Foundations & Applications, Birkhäuser Boston, 2008.

##  Appendix A Proof of Lemma 1

We prove this lemma for two cases, _i.e._ , $r=0$ and $r\geq 1$.

When $r=0$, Taylor’s expansion of $\phi(z)$ at $z=0$ is

$\displaystyle\phi(z)=\phi(0)+R_{0}(z).$

Then, (4) gives $|R_{0}(z)|\leq C\|z\|$.

For any integer $r\geq 1$, one can write the function $\phi(z)$ as

$\displaystyle\phi(z)$ $\displaystyle=\sum_{|\alpha|\leq
r}\frac{\partial^{\alpha}\phi(0)}{\alpha!}z^{\alpha}+R_{r}(z)~{}\text{ or}$
$\displaystyle\phi(z)$ $\displaystyle=\sum_{|\alpha|\leq
r-1}\frac{\partial^{\alpha}\phi(0)}{\alpha!}z^{\alpha}+R_{r-1}(z).$

As a consequence, one has that

$\displaystyle\quad~{}R_{r}(z)$
$\displaystyle=R_{r-1}(z)-\sum_{|\alpha|=r}\frac{\partial^{\alpha}\phi(0)}{\alpha!}z^{\alpha}$
$\displaystyle=r\sum_{|\alpha|=r}\frac{z^{\alpha}}{\alpha!}\int^{1}_{0}(1-t)^{r-1}\partial^{\alpha}\phi(tz)dt-\\!\sum_{|\alpha|=r}\frac{z^{\alpha}}{\alpha!}\partial^{\alpha}\phi(0)$
$\displaystyle=\sum_{|\alpha|=r}\frac{z^{\alpha}}{\alpha!}~{}r\int^{1}_{0}(1-t)^{r-1}\big{(}\partial^{\alpha}\phi(tz)-\partial^{\alpha}\phi(0)\big{)}dt.$

By the condition (4), one has

$\displaystyle|\partial^{\alpha}\phi(tz)-\partial^{\alpha}\phi(0)|\leq
C\|tz\|,~{}~{}t\in(0,1).$

Then, it holds that

$\displaystyle\left|R_{r}(z)\right|$
$\displaystyle\leq\sum_{|\alpha|=r}\frac{|z^{\alpha}|}{\alpha!}~{}r\int^{1}_{0}(1-t)^{r-1}Ct\|z\|dt$
$\displaystyle=rC\|z\|\sum_{|\alpha|=r}\frac{|z^{\alpha}|}{\alpha!}\int^{1}_{0}(1-t)^{r-1}tdt.$

Using integration by parts, one can show that

$\displaystyle\quad\int^{1}_{0}(1-t)^{r-1}tdt$
$\displaystyle=\frac{1}{2}\int^{1}_{0}(1-t)^{r-1}d(t^{2})$
$\displaystyle=\frac{1}{2}\left[\left.(1-t)^{r-1}t^{2}\right|^{1}_{0}-\int^{1}_{0}t^{2}d\big{(}(1-t)^{r-1}\big{)}\right]$
$\displaystyle=\frac{r-1}{2}\int^{1}_{0}(1-t)^{r-2}t^{2}dt$
$\displaystyle~{}~{}\vdots$ $\displaystyle=\frac{(r-1)!}{2\cdots
r}\int^{1}_{0}(1-t)^{r-r}t^{r}dt$ $\displaystyle=\frac{(r-1)!}{(r+1)!}.$

By the multinomial theorem, _i.e._ ,

$\displaystyle\sum_{|\alpha|=r}\frac{r!}{\alpha!}z^{\alpha}=(z_{1}+\cdots+z_{\sigma})^{r},$

and the fact that $|z_{1}+\cdots z_{\sigma}|\leq\sqrt{\sigma}\|z\|$, it holds
that

$\displaystyle\left|R_{r}(z)\right|$ $\displaystyle\leq
rC\|z\|\sum_{|\alpha|=r}\frac{|z^{\alpha}|}{\alpha!}\cdot\frac{(r-1)!}{(r+1)!}$
$\displaystyle=\frac{C\|z\|}{(r+1)!}\sum_{|\alpha|=r}\frac{r!}{\alpha!}|z^{\alpha}|$
$\displaystyle=\frac{C\|z\|}{(r+1)!}|z_{1}+\cdots z_{\sigma}|^{r}$
$\displaystyle\leq\frac{\sigma^{r/2}C\|z\|^{r+1}}{(r+1)!}.$

The proof is complete.∎

##  Appendix B Proof of Lemma 3

The set inclusion condition (8) can be equivalently written as

$\displaystyle\\{x\in\mathbb{R}^{n}:\varphi_{1}(x)\geq 0,\varphi_{2}(x)\geq
0,x\neq 0\\}=\emptyset.$

By Theorem 1, we know that this is true if and only if there exist
$\varphi(x)\in\mathcal{S}_{C}(\varphi_{1},\varphi_{2})$ and
$\zeta(x)\in\mathcal{S}_{M}(x)$, such that

$\displaystyle\varphi(x)+\zeta(x)^{2}=0.$ (38)

Let

$\varphi=s_{0}+s_{1}\varphi_{1}+s_{2}\varphi_{2}$

where $s_{j}$, $j=0,1,2$ are SOS polynomials. By the definition of the cone
$\mathcal{S}_{C}$, one has that
$\varphi\in\mathcal{S}_{C}(\varphi_{1},\varphi_{2})$. Choosing
$\zeta(x)^{2}=x^{\mbox{\tiny\sf T}}x$, we write the condition (38) as

$\displaystyle s_{0}+s_{1}\varphi_{1}+s_{2}\varphi_{2}+x^{\mbox{\tiny\sf
T}}x=0$ (39)

As $s_{0}=-(s_{1}\varphi_{1}+s_{2}\varphi_{2}+x^{\mbox{\tiny\sf T}}x)$ from
(39), if there exist SOS polynomials $s_{1}$ and $s_{2}$ such that the SOS
condition (7) holds, then there exist SOS polynomials $s_{j}$, $j=0,1,2$ such
that (39) is true, and hence the set inclusion condition (8) holds. ∎

##  Appendix C Proof of Lemma 4

For the closed-loop system with the controller $u=Kx$ designed via Theorem 2,
the derivative of the Lyapunov function $V(x)=x^{\top}P^{-1}x$ satisfies

$\displaystyle\dot{V}(x)$
$\displaystyle=x^{\top}P^{-1}(A+BK)x+x^{\top}(A+BK)^{\top}P^{-1}x$
$\displaystyle\quad+2x^{\top}P^{-1}R(x,Kx)$ $\displaystyle\leq-
wx^{\top}P^{-1}x+2x^{\top}P^{-1}R(x,Kx).$

Under Assumption 2, for all $x\in\mathbb{D}$ and $i=1,\dots,n$, the bounds of
the remainder can be found as

$\displaystyle|R_{i}(x,Kx)|\leq\frac{\sqrt{m+n}L_{i}}{2}\|(x,Kx)\|^{2}.$

Hence, for all $x\in\mathbb{D}$, there exists a continuous $\rho_{i}(x)$ for
each $i=1,\dots,n$ such that

$\displaystyle R_{i}(x,Kx)$ $\displaystyle=\rho_{i}(x)\|(x,Kx)\|^{2},$
$\displaystyle\rho_{i}(x)$
$\displaystyle\in\left[-\frac{\sqrt{m+n}L_{i}}{2},\frac{\sqrt{m+n}L_{i}}{2}\right].$

Define $\rho(x)=[\rho_{1}(x)~{}\dots~{}\rho_{n}(x)]^{\top}$. By the definition
of polytopes [39, Definition 3.21], the vector $\rho(x)$ belongs to the
polytope

$\displaystyle\mathcal{H}=\\{\varrho:-\bar{h}\preceq\varrho\preceq\bar{h}\\}$

where

$\bar{h}=[\bar{h}_{1}~{}\cdots~{}\bar{h}_{n}]^{\top}=\begin{bmatrix}\frac{\sqrt{m+n}L_{1}}{2}&\cdots&\frac{\sqrt{m+n}L_{n}}{2}\end{bmatrix}^{\top}.$

Denote $Q_{i}$ as the $i$th column of $P^{-1}$. It holds that

$\displaystyle\quad 2x^{\top}P^{-1}R(x,Kx)$
$\displaystyle=2\begin{bmatrix}x^{\top}Q_{1}&\cdots&x^{\top}Q_{n}\end{bmatrix}\begin{bmatrix}\rho_{1}(x)\|(x,Kx)\|^{2}\\\
\vdots\\\ \rho_{n}(x)\|(x,Kx)\|^{2}\end{bmatrix}$
$\displaystyle=2\sum^{n}_{i=1}x^{\top}Q_{i}\rho_{i}(x)\|(x,Kx)\|^{2}$
$\displaystyle=2\sum^{n}_{i=1}x^{\top}Q_{i}\|(x,Kx)\|^{2}\cdot\rho_{i}(x)$
$\displaystyle=2\begin{bmatrix}x^{\top}Q_{1}\|(x,Kx)\|^{2}&\cdots&x^{\top}Q_{n}\|(x,Kx)\|^{2}\end{bmatrix}\rho(x).$

Denote

$\displaystyle\kappa(x)=\begin{bmatrix}x^{\top}Q_{1}\|(x,Kx)\|^{2}&\cdots&x^{\top}Q_{n}\|(x,Kx)\|^{2}\end{bmatrix}.$

Then, the derivative of the Lyapunov function satisfies for all
$x\in\mathbb{D}$

$\displaystyle\dot{V}(x)\leq-wx^{\top}P^{-1}x+2\kappa(x)\rho(x)$

where $\rho(x)\in\mathcal{H}$. ∎