Region of Attraction Estimation Using Invariant Sets and
Rational Lyapunov Functions ⋆

arXiv:1604.01084v1 [math.OC] 4 Apr 2016

Giorgio Valmorbida a, James Anderson b,
a

Laboratoire des Signaux et Systèmes, CentraleSupélec, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 3 Rue Joliot-Curie,
Gif sur Yvette 91192, France.
b
Department of Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ, UK.

Abstract
This work addresses the problem of estimating the region of attraction (RA) of equilibrium points of nonlinear dynamical
systems. The estimates we provide are given by positively invariant sets which are not necessarily defined by level sets of a
Lyapunov function. Moreover, we present conditions for the existence of Lyapunov functions linked to the positively invariant
set formulation we propose. Connections to fundamental results on estimates of the RA are presented and support the search
of Lyapunov functions of a rational nature. We then restrict our attention to systems governed by polynomial vector fields
and provide an algorithm that is guaranteed to enlarge the estimate of the RA at each iteration.
Key words: Estimates of Region of Attraction, Polynomial systems, Invariant sets, Sum-of-squares.

1

Introduction

The problem of computing the region of attraction (RA)
of asymptotically stable equilibria, or inner estimates
to this set (ERA) [6], is central in several applications
and its relevance is immediately clear for many practical
nonlinear systems for which we can only guarantee local
properties of operating points.
With a converse Lyapunov theorem [22, Theorem 19],
Zubov answered the question “[...] Is it possible, with
the help of the Lyapunov function to find a region of
variation of the initial values x0 such that kφ(t, x0 )k →
0 as t → ∞ ? ” [22, p.3]. The theorem states that if S
is the RA of an equilibrium then the existence of a Lyapunov function (LF) satisfying some conditions on such
a set S is necessary and sufficient. However, computing the LF and the exact RA following Zubov’s theorem
requires the solution of a partial differential equation,
⋆ At the time of writing G. Valmorbida was affiliated to the
Department of Engineering Science and Somerville College,
University of Oxford, Oxford, U.K. J. Anderson is funded
via a Junior Research Fellowship from St. John’s College,
University of Oxford, Oxford, U.K. Work partly supported
by EPSRC grant EP/J010537/1.
Email addresses:
giorgio.valmorbida@l2s.centralesupelec.fr (Giorgio
Valmorbida), james.anderson@eng.ox.ac.uk (James
Anderson).

Preprint submitted to Automatica

which is difficult to obtain in all but simple cases. However, local solutions (in a compact set around the equilibrium point) to the conditions can be obtained more
easily and yield ERAs for the equilibrium point of interest. In this context, a method to approximate solutions
to the conditions of [22, Theorem 19] is obtained with a
series expansion of the LF [22, p.91] and is now referred
to as Zubov’s Method.
In [20], Zubov’s theorem was modified to consider Lyapunov functions mapping Rn to R≥0 (the original result is stated in terms of a map from Rn to the interval
[−1, 0]). One of the conditions in [20] imposes that the
LF V (x) satisfies V (x) → ∞ whenever x → ∂S (the
boundary of the RA) or whenever kxk → ∞. Such a
property is described by the observation that “the candidate must in effect ‘blow up’ near the boundary of the
domain of attraction”. These functions were called maximal Lyapunov functions (MLFs). One of the key observations was that rational functions could be used to
approximate MLFs and therefore be used to obtain estimates of the RA. As a matter of fact, the class of ratio(x)
nal functions of the form V (x) = VVN
where VN and
D (x)
VD are polynomials, were considered as LF candidates
in the algorithm proposed in [20] with the boundary of
the ERA characterised by the set {x ∈ Rn | VD (x) = 0}.
At this point, for the sake of clarity, it is important to
distinguish between two similar sounding yet very differ-

6 April 2016

the ERA (and can be approximated to arbitrary accuracy by rational functions) and sum-of-squares methods
for rational LF construction. Note that rational Lyapunov functions were considered in [5] to obtain MLSs.

ent objects: a maximal Lyapunov function (MLF) and a
maximal Lyapunov set (MLS). An MLF is a Lyapunov
function which satisfies a strict set of conditions (cf. Definition 1 in Section 3). In contrast, a MLS is defined as
the largest level set of a given LF contained in a specified set. Computing the MLS is of course of interest since
one might wish to compute the best ERA achievable for
a given Lyapunov function [3, 5]. Further to the choice
of the class of the LF, conservativeness is introduced by
imposing the level sets of the Lyapunov function to be
the ERA, as observed in [9, p.320]“Estimating the region
of attraction by Ωc = {x|V (x) ≤ c} is simple but usually
conservative. According to LaSalles’s theorem [...] we can
work with any compact set Ω ⊂ D provided we can show
that Ω is positively invariant.” The statement highlights
the fact that contractiveness of the function defining the
ERA is restrictive.

The paper is organised as follows: We present some definitions and the problem statement in Section 2 and describe the main theoretical results in Section 3. Narrowing our attention to systems described by polynomial
vector fields we describe a computational method for
constructing ERAs based on sum-of-squares programming in Section 4 which is illustrated by numerical examples in Section 5.
2

Preliminaries

Let R, R≥0 , R>0 and Rn denote the field of reals, nonnegative reals, positive reals and the n-dimensional Euclidean space respectively. The function f : Rn → R is
positive definite if f (x) > 0 for all non-zero x ∈ Rn ,
similarly if f (x) ≥ 0 for all x ∈ Rn then f is positive
semidefinite. The set of functions g : Rn → R which is
n-times continuously differentiable is denoted C n . co(X )
denotes the convex hull of the set X , X ◦ its interior, ∂X
its boundary, and X its closure. The minimum (maximum) of a scalar function S(x) in a compact set Y
is denoted min(S(x)) (max(S(x))). We also use max

In recent years, sufficient conditions for local stability
analysis, requiring invariance and contractiveness of a
set led to numerical methods for the estimation of the
RA with polynomial Lyapunov functions [15–17].These
methods rely on the solution of non-convex sum-ofsquares (SOS) constraints constructed with the Positivstellensatz [10, Theorem 2.14]. The solutions to these
problems require a coordinate-wise search since the
non-convex nature results from the fact that some polynomial variables appear multiplying the Lyapunov function which is itself a variable. For a detailed description
of sum-of-squares methods for RA estimation the reader
is referred to [4]. For the case of a given LF, the computation of the MLS was pursued in [5]. In [8] the theory of
moments is used to estimate the RA of uncertain polynomial systems. We also find in the literature numerical
methods exploiting topological properties of the boundary of the RA requiring the computation of trajectories
and equilibrium points. However the complexity of such
methods has restricted them to 2-dimensional examples [6]. Recently, in [21], set advection methods are
described for polynomial systems.

x∈Y

x∈Y

to denote the function taking the maximum of its arguments. For x ∈ Rm the ring of polynomials in m
variables is denoted by R[x]. For p ∈ R[x], deg(p) denotes the degree of p. A polynomial p(x) is said to be a
sum-of-squares if there exists a finite set of polynomials
P
g1 (x), . . . , gk (x) such that p(x) = ki=1 gi2 (x). The set of
SOS polynomials in x is denoted by Σ[x1 , . . . , xm ] which
can be abbreviated to Σ[x]. Equivalently p(x) is SOS if
there exists a positive semidefinite matrix Q such that
p(x) = Z T (x)QZ(x) where Z(x) is a vector of monomials [13]. Note that the search for Q can be formulated as
a semidefinite programme and thus solved using convex
optimization techniques [19].

In this paper we derive conditions based on Lyapunov
stability results that guarantee that trajectories initiated
from an positively invariant set converge to a level set of
the LF which is contractive and invariant therefore guaranteeing such a postively invariant set to be an ERA. In
addition to the positively invariant estimates, we present
conditions to obtain LF certificates of a specific form
which specializes to rational functions in case of polynomial data. We then propose a numerical method based
on the solution of SOS constraints for the case of polynomial systems and estimates in the form of semi-algebraic
sets (sets defined by polynomial constraints). The work
in this paper extends the work of [18] and connects the
concept of maximal Lyapunov functions [20] to polynomial optimization techniques based on sum-of-squares
programming. To the best of the authors knowledge this
is the first work to offer a theoretical link between maximal Lyapunov functions, which completely characterise

Consider the dynamical system
ẋ = f (x)

(1)

where f : D → Rn is a locally Lipschitz map from a domain D ⊂ Rn to Rn , with 0 ∈ D. Let us assume x = 0
is an equilibrium point, i.e. 0 ∈ {x ∈ Rn |f (x) = 0}.
Denote by φ(t, x(0)) the solution to (1) that is initiated from the point x(0) at time t = 0, the set L is
said to be invariant with respect to (1) provided x(0) =
φ(0, x(0)) ∈ L ⇒ x(t) = φ(t, x(0)) ∈ L, ∀t ∈ R. Furthermore, L is said to be positively invariant with respect to (1) if the previous implication holds for all
t ≥ 0. Given a function R : Rn → R we define the set
E(R, γ) := {x ∈ Rn | R(x) ≤ γ} for some γ > 0 and
E ◦ (R, γ) := {x ∈ Rn | R(x) < γ} for some γ > 0. Additionally, provided that a function V : Rn \ {0} → R>0

2

satisfies V̇ (x) = ∂V
∂x f (x) < 0 on E(V, γ) then the set
E(V, γ) is said to be contractive and invariant, furthermore the function V is said to be a Lyapunov function [9,
Chapter 4]. We assume throughout this work that any
function used to define a contractive set is in C 1 . The
region of attraction of an asymptotically stable equilibrium point x∗ of (1) is defined as the set
S :=

(

x∈R

n

φ(t, x) is defined ∀t ≥ 0,
limt→∞ φ(t, x) = x∗

)

,

hold and there exists a positive scalar ǭ such that 0 ∈
/
E(R, γ) \ E ◦ (R, γ − ǭ) then
(II) the function
V (x) =

Proof. Proof of (I). Part 1 (Contractiveness of a level set
of VN ): by assumption E(R, γ) is compact thus we can
compute α =
min VN (x), then from (4a), α > 0.
x∈∂E(R,γ)

Define U = E(VN , α) ∩ E(R, γ), from (4a), (4b) we have

The focus of this paper is to construct inner estimates
of S by computing positively invariant sets.

VN (x) > 0 ∀x ∈ U \ {0} ,
− h∇VN (x), f (x)i > 0 ∀x ∈ U \ {0} .
Following [9, Theorem 4.1] we have that the origin is
asymptotically stable and an inner approximation of
its region of attraction is given by U, that is trajectories φ(t, x0 ) with x0 ∈ U exist, are unique, and satisfy
φ(t, x0 ) → 0 as t → ∞.

Main results

In this section we present conditions to certify that a
compact set is a postively invariant set and provides an
estimate of the RA for the origin of (1). It is also shown
how, under a different condition, to obtain an alternative
Lyapunov certificate. We then extend these results to
the case where the system under study is affected by
parametric uncertainty.
3.1

Part 2 (Positively Invariance of E(R, γ)): since E(R, γ)
is compact and f (x) is locally Lipschitz in any compact
set, we have existence and uniqueness of solutions to
ẋ = f (x), for all x0 ∈ E(R, γ), provided every solution
lies in E(R, γ). Let us prove that for all T ∈ [0, ∞) we
have x(T ) ∈ E(R, γ). Assume there exists x0 ∈ E(R, γ)
for which the solution leaves the set, then, there must
exist a T ∗ such that x(T ∗ ) = φ(T ∗ , x0 ) that satisfies
R(x(T ∗ )) > γ. From the continuity of solutions and continuity of R(x) there must exist T̄ , 0 < T̄ < T ∗ such that
R(x(T̄ )) = γ and Ṙ(x(T̄ )) = h∇R(x(T̄ )), f (x(T̄ ))i ≥ 0,
which contradicts (3b). Hence E(R, γ) is a positively invariant set.

Region of Attraction Estimates

The following theorem is used to verify that a compact
set is positively invariant and defines an estimate of the
RA of the equilibrium point at the origin and to obtain
functions of which the denominator provides the RA estimate.
Theorem 1 Given R : Rn → R, R ∈ C 1 and γ > 0,
satisfying
E(R, γ) ⊂ D is compact and 0 ∈ E(R, γ),
− h∇R(x), f (x)i > 0 ∀x ∈ ∂E(R, γ),

(6)

is a Lyapunov function for (1) and gives E ◦ (R, γ) as an
estimate of S.

(2)

without loss of generality, we will assume throughout
this paper that the equilibrium point of interest is at the
origin, i.e. x∗ = 0.

3

VN (x)
γ − R(x)

Part 3 (Convergence of trajectories starting in E(R, γ) \
E ◦ (VN , α) to U): Finally let us now prove that every
trajectory satisfying x(0) ∈ E(R, γ) \ E ◦ (VN , α) enters
the positively invariant and contractive set U, that is,
that there exists a T such that x(T ) ∈ U. Let β =
max VN (x). Since (4b) holds in E(R, γ) \ E ◦ (VN , α)

(3a)
(3b)

if there exists VN : Rn → R, VN (0) = 0, VN ∈ C 1 , such
that
VN (x) > 0 ∀x ∈ E(R, γ) \ {0},
(4a)
− h∇VN (x), f (x)i > 0 ∀x ∈ E(R, γ) \ {0} ,
(4b)
then

x∈∂E(R,γ)

let λ = −

max

x∈E(R,γ)\E ◦ (VN ,α)

h∇VN , f (x)i, which exists be-

cause the continuous function has a maximum over any
compact set, from (4b) we also get

(I) the solutions x(t) = φ(t, x0 ) to (1), with x0 ∈
E(R, γ) satisfy x(t′ ) ∈ S ∀t′ ∈ [0, ∞).

VN (x(t)) = VN (x(0))+

Z

t

V̇N (x(τ ))dτ ≤ VN (x(0))−λt.

0

Moreover, if (3a), (4a) and

Since VN (x(0)) ≤ β we have VN (x(t)) ≤ β −λt in the set
E(R, γ) \ E ◦ (VN , α). This implies that ∃T ≥ 0 satisfying
T ≤ β−α
λ such that VN (x(T )) = α and hence x(T ) ∈ U.

− h[∇VN (x)(γ − R(x)) + VN (x)∇R(x)] , f (x)i > 0
∀x ∈ E(R, γ) \ {0} (5)

3

then it is immediate that δ > 1. As ǫ ≤ ǭ we
have ∂E(R, γ − ǫ) ⊂ E(R, γ) \ E ◦ (R, γ − ǭ), thus
min
VN (x) ≥
min◦
VN (x) and hence

Proof of (II): From (4a) we have that V (x) > 0 ∀x ∈
E ◦ (R, γ) \ {0}. The time-derivative of (6) along the trajectories of (1) is given by
V̇ (x) =

h∇VN (x)(γ − R(x)) + VN (x)∇R, f (x)i
(γ − R(x))2

x∈E(R,γ)\E (R,γ−ǭ)

x∈∂E(R,γ−ǫ)

(7)

δ>

VN (x)
, ∀x ∈ E(R, γ).
min
VN (x)

x∈∂E(R,γ−ǫ)

which satisfies −V̇ (x) > 0 ∀x ∈ E ◦ (R, γ)\{0} if (5) holds
true. Thus V is a Lyapunov function for the equilibrium
point at the origin of (1). What is left to show is that the
level curves of V define the stated estimate of the ERA.

Since ǫ is positive, the set containment


E(R, γ−ǫδ) ⊂ x ∈ Rn | R(x) ≤ γ − ǫ


We have that V (x) and V̇ (x) are not defined in ∂E(R, γ),
thus it is not possible to compute scalars α and β Parts 1
and 3 in the above Proof of (I). Consider the arbitrarily
small scalar, ǫ that defines the set E(R, γ −ǫ) and assume
it satisfies ǫ ≤ ǭ. Then, as E(R, γ) is compact it follows
that so is E(R, γ − ǫ). We can then follow the same steps
of Proof of (I), Part 1 and compute a positive scalar
αǫ =
min
V (x). Define Uǫ = E(V, αǫ )∩E(R, γ−ǫ),




VN (x)
min
VN (x) 

x∈∂E(R,γ−ǫ)

holds true. Using (8) we have E(R, γ − ǫδ) ⊂ E(V, αǫ ),
thus E(R, γ−ǫδ) = E(R, γ−ǫδ)∩E(R, γ−ǫ) ⊂ E(V, αǫ )∩
E(R, γ − ǫ) = Uǫ which, using the fact that Uǫ ⊂ E(R, γ),
we obtain

x∈∂E(R,γ−ǫ)

since V (x) > 0 ∀x ∈ E ◦ (R, γ) \ {0} and −V̇ (x) > 0 ∀x ∈
E ◦ (R, γ) \ {0} we have

E(R, γ − ǫδ) ⊂ Uǫ ⊂ E(R, γ).

V (x) > 0 ∀x ∈ Uǫ \ {0} ,
− h∇V (x), f (x)i > 0 ∀x ∈ Uǫ \ {0} .

Thus E(R, γ −ǫδ) is an estimate of the RA of (1). Since δ
is bounded and ǫ can be chosen to be arbitrarily small,
we have that the set E ◦ (R, γ) is an estimate of the RA
of (1).


Then following [9, Theorem 4.1] the origin is asymptotically stable and an inner approximation of its region
of attraction is given by Uǫ , that is trajectories φ(t, x0 )
with x0 ∈ Uǫ exist, are unique, and satisfy φ(t, x0 ) → 0
as t → ∞.

A consequence of the the assumption imposed in Theorem 1 that the set E(R, γ) is compact, is that E(R, γ) will
be connected. This follows from the fact that (3b) makes
E(R, γ) a positively invariant set, and by (4a), VN (x) is
strictly positive on E(R, γ) \ {0}, and finally (4b) ensures its derivative is strictly negative on E(R, γ) \ {0}.

Using (6), we obtain


VN (x)
n
E(V, αǫ ) = x ∈ R | R(x) ≤ γ −
.
αǫ

Remark 1 Note that R(x) is not required to be positive definite, however it is required that E(R, γ) is compact and contains the origin. This requirement guarantees that, in the proof of Part 1, min∂E(R,γ) (V (x)) is
well-defined, such that E(V, α) ⊆ E(R, γ).
y

Since
αǫ =

min

V (x)

x∈∂E(R,γ−ǫ)

VN (x)
− R(x)
VN (x)
=
min
x∈∂E(R,γ−ǫ)
γ
−
(γ − ǫ)


1
min
VN (x)
=
ǫ x∈∂E(R,γ−ǫ)

=

min

Remark 2 The boundedness of the set E(R, γ) also
guarantees the uniqueness of solutions in the set if the
vector field is not globally Lipschitz (as for instance, the
polynomial vector fields).
y

x∈∂E(R,γ−ǫ) γ

we then have


E(V, αǫ ) = x ∈ Rn | R(x) ≤ γ − ǫ


Remark 3 In the case R(x) = VN (x), (4b) implies (3b)
and the set E(VN , γ) has to be compact as required
by (3a).
y




VN (x)
.
min
VN (x) 

x∈∂E(R,γ−ǫ)

The proposition below presents sufficient conditions to
satisfy the constraints of Theorem 1 formulated in terms
of inequalities and the definition of E(R, γ).

(8)

Define a positive scalar δ as follows:

Proposition 1 Given R ∈ C 1 , R : Rn → R, γ > 0,
satisfying

max VN (x)

δ :=

x∈E(R,γ)

min

x∈E(R,γ)\E ◦ (R,γ−ǭ)

VN (x)

E(R, γ) is compact, 0 ∈ E(R, γ),

4

(9)

if there exist VN : Rn → R, VN ∈ C 1 , VN (0) = 0 and
m0 : Rn → R≥0 , m1 : Rn → R≥0 , p : Rn → R, such that

The main result from [20] is summarised by the following
theorem:

− h∇R(x), f (x)i > p(x)(γ − R(x)) ∀x ∈ D
VN (x) > m0 (x)(γ − R(x)) ∀x ∈ D′
− h∇VN (x), f (x)i > m1 (x)(γ − R(x)) ∀x ∈ D′

Theorem 2 Suppose we can find a set A ⊆ Rn which
contains x = 0 in its interior, a continuously differentiable function Vm : A → R>0 and a positive definite
function ψ(x) such that

(10a)
(10b)
(10c)

where D′ := D\{0}, then E(R, γ) is an ERA of the origin.
If (9), (10b), (10a) hold and there exist m2 : Rn → R≥0
such that

(1) Vm (0) = 0, V (x) > 0 if x ∈ A \ {0},
(2) V̇m (x) = −ψ(x) for all x ∈ A,
(3) Vm (x) → ∞ as x → ∂A and/or kxk → ∞,

−h∇VN (x), f (x)i+VN (x)p(x) > m2 (x)(γ−R(x)) (11)

then A = S.

holds, then (6) is a Lyapunov function for (1) and
E ◦ (R, γ) is an ERA of the origin.

Clearly an MLF satisfies the properties of Theorem 2,
however in [20] the authors constructively provide a
method for extending any LF for (1) into an MLF (the
assumption is that ψ(x) has been constructed). They
(x)
then show that rational LFs of the form V (x) = VVN
D (x)
with polynomial numerator and denominator can arbitrarily approximate an MLF. Furthermore the set
{x | VD (x) = 0} defines the boundary of S.

Proof. From the non-negativity of m0 (x), and m1 (x) we
have that (10b), (10c) imply that (4a), (4b) hold. Since
at ∂E(R, γ) we have γ − R(x) = 0, (10a) implies (3b),
and, according to Theorem 1 Claim (I), E(R, γ) is an
ERA of the origin. The time derivative of V (x), as in (6),
is given by (7). Since we have VN (x) > 0 in E(R, γ) if
(10b) holds, inequality (10a) provides a lower bound for
−V̇ as follows

In comparison to Theorem 2, Theorem 1 shows that the
zero level set of a rational Lyapunov function can be used
to compute an estimate of the RA. The main difference
between the two results is that our theorem provides an
estimate of the RA whilst the MLF provides the exact
RA (albeit at the cost of having to solve a partial differential equation, namely item 2 of Theorem 2). Thus the
LF (6) is not necessarily an MLF. In the remainder of
this section we develop further the results of Theorem 1.

h∇VN (x), f (x)i VN (x) h∇R(x), f (x)i
−
(γ − R)
(γ − R)2
h∇VN (x), f (x)i VN (x)p(x)
> −
+
,
(γ − R)
(γ − R)

−V̇ (x) = −

therefore, if there exists m2 : Rn → R>0 satisfying
(− h∇VN (x), f (x)i + VN (x)p(x))
> m2 (x),
(γ − R)

3.2

The following result parallels Theorem 1 and considers
positively invariant regions defined by the maximum of
a set of differentiable functions.

that is (11), then −V̇ (x) > 0 (and (5) holds true). Following Theorem 1, Claim (II), V (x) is a Lyapunov function for (1) in E ◦ (R, γ), thus providing an ERA for the
origin of (1).


Let d be a finite positive integer and functions
Ri (x), . . . , Rd (x) be given. The point-wise maximum of
these functions at x is defined as

Before developing the theory further in order to take
into account non-smooth set descriptions and ERAs for
uncertain systems we first compare the above conditions to the classical results of Zubov via Vannelli and
Vidyasagar’s maximal Lyapunov function [20] framework which characterises the region of attraction S.

RM (x) := max(R1 (x), . . . , Rd (x)).

(12)

Theorem 3 Given Ri : Rn → R, Ri ∈ C 1 i = 1, . . . , d
and a positive scalar γ, satisfying

Definition 1 A function Vm : Rn → R>0 ∪ {∞} that
for the system (1) satisfies
(1)
(2)
(3)
(4)

Piece-wise Lyapunov Functions

E(RM , γ) is compact and 0 ∈ E(R, γ),

Vm (0) = 0, Vm (x) > 0 if x ∈ S \ {0},
Vm (x) < ∞ iff x ∈ S,
Vm (x) → ∞ as x → ∂S and/or kxk → ∞,
V̇m (x) < 0 and well defined for all x ∈ S \ {0},

(13a)

∂RM (x)
(13b)
∂x
M (x)
denotes the subdifferential of RM (x) at x,
where ∂R∂x
if there exists a function VN : Rn → R≥0 , VN (0) = 0,
− hξ, f (x)i > 0 ∀x ∈ ∂E(RM , γ), ∀ξ ∈

is called a maximal Lyapunov function.

5

VN ∈ C 1 such that
VN (x) > 0 ∀x ∈ E(RM , γ) \ {0}
− h∇VN (x), f (x)i > 0 ∀x ∈ E(RM , γ) \ {0}
then

to satisfy the conditions of Theorem 3. These conditions
are formulated in terms of inequalities and the description of the set ∂E(RM , γ).

(14a)

Proposition 2 Given Ri : Rn → R, Ri ∈ C 1 i =
1, . . . , d and a positive scalar γ, if there exist VN : Rn →
R, VN ∈ C 1 , VN (0) = 0, and m0 : Rn → R≥0 , m1 :
Rn → R≥0 , p : Rn → R such that

(14b)

(I) all trajectories of (1) initiated from the set
E(RM , γ) converge to the origin.

E(RM , γ) is compact, 0 ∈ E(RM , γ)
∂RM (x)
− hξ, f (x)i > p(x)(γ − RM (x)) ∀ξ ∈
∂x
VN (x) > m0 (x)(γ − RM (x))
− h∇VN (x), f (x)i > m1 (x)(γ − RM (x))

Moreover, if (14a), and
−h[∇VN (x)(γ − RM (x)) + VN (x)∇RM (x)] , f (x)i > 0
∀x ∈ E(RM , γ) (15)
hold, and there exists a positive scalar ǭ such that 0 ∈
/
E(RM , γ) \ E ◦ (RM , γ − ǭ) then

(16)

3.3

Remark 4 To characterise the set E(RM , γ), notice
that RM (x) = Ri (x) ∀x ∈ {x ∈ Rn |Ri (x) − Rj (x) ≥
0, j = 1, . . . , d}. By defining

i=1

Mi (RM , γ).

Uncertain systems

Consider uncertain dynamical systems of the form
ẋ = f (x, θ),

Mi (RM , γ) := {x ∈ Rn |Ri (x) ≤ γ,
Ri (x) − Rj (x) ≥ 0, j = 1, . . . , d}.
Td

(17c)
(17d)

− h∇VN (x), f (x)i + VN (x)p(x) > m2 (x)(γ − RM (x))
(18)
hold, then (16) is a LF for (1) in the set E ◦ (RM , γ).

is a Lyapunov function for (1) and gives
E ◦ (RM , γ) as an estimate of S.

we can write E(RM , γ) =

(17b)

then E(RM , γ) is an ERA of the origin. If (17a), (17c),
(17b) hold and there exist m2 : Rn → R≥0 such that

(II) the function
VN (x)
V (x) =
γ − RM (x)

(17a)

θ ∈ Θ ⊂ Rnp

(19)

where f : D × Θ → Rn and Θ denotes the uncertainty
set. We assume that f satisfies conditions so as to provide
uniqueness and local existence of solutions 1 . We shall
also assume that x∗ = 0 is the equilibrium of interest,
and require that f (0, θ) = 0 ∀θ ∈ Θ.

y

The subdifferential for the function RM (x) is defined as
∂RM (x)
:= co{∇Rℓ (x), ℓ ∈ I(x)}, where I(x) = {i ∈
∂x
{1, . . . , d}|Ri (x) = RM (x)} denotes the set of “active”
functions at point x. Notice that RM (x) is not differentiable at points x where ∃i, j ∈ I(x), i 6= j, that is,
at points satisfying RM (x) = Ri (x) = Rj (x), i 6= j. At
M (x)
such points ∂R∂x
defines a set, hence (13b) describe a
set of inequalities. Whenever I(x) contains only one elM (x)
is a singleton given by
ement, say I(x) = {k}, ∂R∂x
∇Rk , which exists since Ri (x) ∈ C 1 , ∀i = 1, . . . , d.

We are interested in determining a robust estimate for
region of attraction i.e. an estimate of the RAs for all
dynamical systems of the form (19) with a fixed θ ∈ Θ

Sθ :=

(

x0 ∈ R

n

φ(t, x, θ) is defined ∀t ≥ 0
limt→∞ φ(t, x0 , θ) = x∗ , ∀θ ∈ Θ

)

,

where φ(t, x0 , θ) is a solution to (19) starting from x0
at time t with fixed θ ∈ Θ. We establish conditions for
positively invariant sets to be estimates of the region of
attraction for parametrically uncertain systems. Whilst
we will consider parameter dependent Lyapunov functions (PDLFs), our ERAs will be defined by positively
invariant sets which are not dependent on the system
parameters.

The proof of Theorem 3 follows closely the proof of Theorem 1 and is therefore omitted. The only difference is
related to the lackDof differentiabilityEof RM (x) which
(x(t))
, f (x(t)) . Therefore progives ṘM (x(t)) ∈ ∂RM∂x
vided (13b) holds we can use it to arrive at a contradiction as in the proof of Claim (I) of Theorem 1.

1
This requires f to be continuous in (x, θ, t) and locally
Lipschitz in x (uniformly in θ and t) on a bounded domain.
Exact conditions can be found in [9, Theorem 3.5]

The proposition below parallels Proposition 1 and is presented without proof. It introduces sufficient conditions

6

4.1

PDLFs have been shown to be an effective tool for certifying the stability of linear system with parametric uncertainties [2, 7, 11]. They have also been successfully
applied to obtain certificates for the local stability of
polynomial systems leading to parameter-dependent estimates of the RA. In those results, a robust estimate is
then obtained as the intersection of the estimates given
for each fixed parameter value [3,17], which is in contrast
to the result in this section where the estimate does not
depend on the parameters but the Lyapunov function
does, thus avoiding the computation of the intersection
of the parametrised estimates.

Recall that for a locally asymptotically stable equilibrium point of (1), converse Lyapunov theorems tells us
there exist a Lyapunov function V and a set D ⊂ Rn ,
0 ∈ D, satisfying V : D → R such that V (x) > 0
∀x ∈ D \ {0}, V (0) = 0 and h∇V (x), f (x)i < 0
∀x ∈ D \ {0}, see for example [9, Section 4.7]. Based on
this fact, a common (but conservative) approach to finding an ERA for systems of the form (1) is to compute a
Lyapunov function certifying the local asymptotic stability of the origin and obtain the largest (maximal) level
set of a Lyapunov function that is contained within the
set D which the LF is constructed on. We can describe a
general algorithm based on the Lyapunov function computations to obtain ERAs in the form E(V, γ) as follows:

The following result extends Theorem 1 to the case of
uncertain systems of the form (19).
Theorem 4 Consider the uncertain dynamical system
described by (19) where Θ is a compact set and x∗ = 0.
Given a function R : Rn → R, R ∈ C 1 and a positive
scalar γ, which defines a compact set E(R, γ), and satisfy

Algorithm 1
Input k = 0, a compact set D0 , {0} ∈ D0 .
Step 1 (Lyapunov function computation): Given Dk ,
compute a Lyapunov function Vk for (1).
Step 2 (Maximization of Lyapunov level set): Given Vk
and Dk , solve

− h∇R(x), f (x, θ)i > 0 ∀(x, θ) ∈ ∂E(R, γ) × Θ, (20)
if there exists a function VN : Rn ×Θ → R, VN (0, ·) = 0,
VN ∈ C 1 such that
VN (x, θ) > 0 ∀(x, θ) ∈ E(R, γ) \ {0} × Θ

maximize γ

Whilst the algorithm above may look simple enough,
observe that: i) in general, constructing Lyapunov functions for nonlinear dynamical systems is a non-trivial
task. ii) Existing formulations for the optimization problem in Step 2 are typically non-convex. iii) Determining the update for Dk+1 typically relies on some heuristic. As a general rule, by necessity {0} ∈ Dk+1 , the set
should be connected, and should contain points that are
not already in Dk . In Section 4.3 we will specialise the
above algorithm and explicitly describe how to update
the search domain, and in the case of polynomial systems construct all the required functions.

−h[∇VN (x, θ)(γ − R(x)) + VN (x, θ)∇R(x)] , f (x)i > 0
∀(x, θ) ∈ ∂E(R, γ) × Θ (22)
hold, and there exists a positive scalar ǭ such that 0 ∈
/
E(R, γ) \ E ◦ (R, γ − ǭ) then the function
(23)

is a Lyapunov function for (19) for all θ ∈ Θ in the set
E ◦ (R, γ).

Step 1 asks for h∇Vk (x), f (x)i < 0 ∀x ∈ Dk . In general there is no guarantee that every x(0) ∈ Dk satisfies
φ(t, x(0)) ∈ Dk ∀t > 0. Extra conditions must hold for
Dk to be a positively invariant set. We will expand upon
this point in the next section.

The proof is similar to that of Theorem 1 and so is omitted.
4

Computational
the RA

Methods

for

(24)

Step 3 (Update of search domain): If stopping criteria
is satisfied then return E(Vk , γ ∗ ), with γ ∗ the solution
to (24), as the ERA else specify Dk+1 , set k ← (k + 1)
go to Step 1.


Moreover, if (21), and

VN (x, θ)
γ − R(x)

subject to E(Vk , γ) ⊂ Dk , γ > 0

(21)

− h∇VN (x, θ), f (x, θ)i > 0 ∀(x, θ) ∈ E(R, γ) \ {0} × Θ
then the solutions to (19), φ(t, x0 , θ) for any x0 ∈ E(R, γ)
and θ ∈ Θ lie in the set Sθ with respect to x∗ = 0.

V (x, θ) =

Maximal Lyapunov Sets

Estimating
Denote the ERA and the search domain from iteration k
of the above algorithm by E(Vk , γk ) and Dk respectively.
Note that Dk ⊂ Dk+1 does not necessarily guarantee
E(Vk , γk ) ⊂ E(Vk+1 , γk+1 ). Satisfying such constraints is
of desirable as it guarantees improvement of the ERA.
For polynomial systems, we describe next how such criteria can be enforced.

We now present a computational method for constructing positively invariant estimates of the RA. First a
method for estimating the RA via maximal Lyapunov
sets is reported and then algorithms that implement the
main results of the paper are described.

7

4.2

sum-of-squares polynomial m2 such that

Estimating the RA With Invariant Sets

−h∇VN (x), f (x)i+VN (x)p(x)−m2 (x)(γ −R(x)) ∈ Σ[x]
(26)
hold, then (6) is a LF for (1) in the set E ◦ (R, γ).

We now illustrate how Theorem 1, Claim (I) can be implemented in an algorithmic manner to obtain estimates
for the RA. We start by presenting a generic algorithm,
analogous to Algorithm 1 with the exception that the
obtained ERA is not given by Lyapunov level sets:

As proven in Proposition 1, the fact that (25a) and (26)
hold is a sufficient condition for (5) and, if satisfied, it
VN
certifies that the rational function V = γ−R
is a Lya◦
punov function on the set E (R, γ).

Algorithm 2
Input: k = 0, a function R0 satisfying the (3a) in Theorem 1.
Step 1 (Invariant set enlargement ): Maximize γ subject to (3b), (4) with R = Rk .
Step 2 (Update function R): If stopping criteria is satisfied then return ERA given by the set E(Rk , γ); else
compute R∗ satisfying E(R∗ , δ) ⊇ E(Rk , γ ∗ ) , where γ ∗ is
the optimal solution to Step 1 and δ > 0. Set k ← k + 1;
γ ← δ ; Rk ← R∗ ; go to Step 1.


Sum-of-squares constraints such as those above can be
formulated using freely available software such as SOSTOOLS [12] and solved using a semidefinite programme
solver. Note that in (25), VN and R appear affinely in the
SOS constraints, and the only product between VN and
the set of multipliers is in (26) where it multiplies polynomial p. This fact is central in developing coordinatewise search algorithms allowing VN to be a variable at
all steps of the algorithm iterations. Notice that is only
possible because in (25) there is no product between VN
and the multipliers m0 , m1 and m2 . When searching for
a rational LF certificate (by imposing (26)), the product
VN p is handled by fixing p from a solution to (25a).

Note that any function R0 that satisfies (3) can be used
to initiate Algorithm 2 provided the set E(R0 , γ) is invariant. A straightforward choice for R0 is any Lyapunov
function V , which provides a (possibly arbitrarily small)
level set E(V, γ) which is invariant and also contractive
(although such a property is not required for E(R0 , γ)).
We take such an approach in the examples in Section 5.

VN
Remark 5 Notice that, for a given rational V = γ−R
,
(assuming polynomial dependence of VN and R on x)
straightforward formulations to compute the MLS of the
function V , as the set E(V, C), lead to inequalities as

Observe that the update of function R parallels the update of the domain in Algorithm 1, with the difference
that it defines an ERA itself. In the next subsection we
describe a specific strategy for the update R∗ of function
R from step 2.
4.3

−

Polynomial systems

For the remainder of the paper it is assumed that the
vector field f in (1) and the functions VN , R and RM in
Theorems 1 and 3 are multivariate polynomials. For this
class of systems, [4] presents a comprehensive set of results on estimates of the domain of attraction with polynomial Lyapunov functions with Sum-of-Squares based
approaches [1].

h∇VN (x)(γ − R(x)) + VN (x)∇R(x), f (x)i
(γ − R)2


VN (x)
.
≥ m(x) C −
γ − R(x)

with m(x) ≥ 0. In order to avoid the rational inequalities
of the form above, one can restrict the attention to the
set where (γ − R(x)) > 0, thus formulating polynomial
inequalities as
− h∇VN (x)(γ − R(x)) + VN (x)∇R(x), f (x)i
≥ m(x)(C(γ − R(x))2 − VN (x)(γ − R(x))).

The following corollary to Proposition 1 provides sufficient conditions for Theorem 1 to hold which are verifiable using convex optimization.

(27)

Notice also that, in contrast to (27), (26) does not present
products VN (x)R(x) and R(x)2 . Such a property also
allows the degree of the SOS constraints in (26) to be
lower than a SOS constraint obtained with (27).
y

Corollary 1 Let VN and R be given multivariable polynomials and γ a given positive constant. Then, if there
exist sum-of-squares polynomials m0 , m1 and a polynomial p such that
(25a)
(25b)
(25c)

Remark 6 The analysis of polynomial systems with
polynomial LFs and level sets of the LF as ERAs are
a particular case of the conditions imposed and are obtained by imposing R(x) = VN (x), m0 (x) = 0, and
p(x) = m1 (x) = m2 (x).
y

then the inequalities (10) are satisfied and E(R, γ) is an
ERA of the origin. If (25a), (25b) hold and there exist a

The proof of the following claims are straightforward
and, therefore, omitted.

− h∇R(x), f (x)i − p(x)(γ − R(x)) ∈ Σ[x]
VN (x) − m0 (x)(γ − R(x)) ∈ Σ[x]
− h∇VN (x), f (x)i − m1 (x)(γ − R(x)) ∈ Σ[x]

8

choice for R0 is the Lyapunov function VN produced by
the algorithm (see Example 2 in Section 5).

Corollary 2 If (25a) holds with p ∈ R[x], that also satisfies p ∈ Σ[x] (i.e. a sum-of-squares polynomial), then
(26) is holds true with m2 = m1 + pm0 .

After using Algorithm 3 to construct an ERA, it is desirable to compute the corresponding rational LF of the
form (6). In order to do so, using the multiplier p, functions R and γ from the final iteration of Algorithm 3
solve the feasibility problem

b R and a scalar
Proposition 1 Given polynomials R,
γ > 0, if there exists a sum-of-squares polynomial m3 (x)
such that
b
(γ − R(x)) − m3 (x)(γ − R(x))
∈ Σ[x], m3 ∈ Σ[x]. (28)

find

b γ) ⊆ E(R, γ).
then E(R,

VN ,m0 ,m1 ,m2

γ

subject to

(25), (28), mi ∈ Σ[x], i ∈ {0, 1}.

VN ,R,mi ,p

5

From (28), any solution γ ∗ to the above problem guaranb γ∗)
tees that the set E(R, γ ∗ ) is contained in the set E(R,
which is an ERA since conditions in Theorem 1 are satisfied when (25) is satisfied.

Example 1: Consider system (1) with

Algorithm 3
Input: k = 0, a function R0 satisfying the condition (3a) of (I) in Theorem 1.
Step 1: With R = Rk , solve through a line search on γ:
maximize

γ

f (x) =

Step 2: If stopping criterion is satisfied, return E(Rk , γ)
else, using m0 , m1 , p and the optimal γ ∗ from Step 1, set
b
R(x)
= R(x) and solve with a line search on γ:
VN ,R,m3

"

−0.42x1 − 1.05x2 − 2.3x21 − 0.5x1 x2 − x31
1.98x1 + x1 x2

#

which satisfies f (0) = 0. This system was studied
in [20, Example 4] in the context of maximal LFs and
in [15] with composite Lyapunov functions. It describes
a Lotka-Volterra system with its stable equilibrium
point translated to the origin.

subject to (25b)-(25a).

maximize γ subject to (25), (28), γ ≥ γ ∗ ,

Numerical Examples

We now illustrate our results with three numerical examples, we use SOSTOOLS and the semidefinite programme solver SeDuMi [14].

The algorithm below exploits this fact to specialize Algorithm 2 and obtain ERAs for polynomial systems by
solving a sequence of sum-of-squares programmes.

VN ,m0 ,m1 ,p

subject to (25), (26), mi ∈ Σ[x],

i ∈ {0, 1, 2}. One important assumption in Theorem 1
is the compactness of set E(R, γ). In order to enforce
this property when computing R(x) ∈ R[x] we impose
a constraint of the form R(x) ≥ c(x) + κkxk2k where
c(x) ∈ R[x], deg(c) ≤ 2k − 1, κ ∈ R≥0 , k ∈ N.

The constraints (25) and (28) can be used to formulate
the following bilinear sum-of-squares programme
maximize

VN

With the initial function R obtained from the ERA
from [20], we apply Algorithm 3 (thus allowing the Lyapunov function VN to be a variable at each step). With
the obtained positively invariant set E(R, γ) defining the
ERA and a multiplier p solving (25a), we compute VN
satisfying (25b), (25c), (26) and mi ∈ Σ[x], i ∈ {0, 1, 2}
thus yielding a rational LF as (6). The boundary of the
ERA and a sequence of nested level sets of the rational
LF V (x) are depicted in Figure 1. Figure 2 depicts the
sets {x|V̇N = 0} and the set {x|Ṙ = 0}, illustrating that
constraints (25c) and (25a) hold, that is, the intersection
of ∂E(R, γ) with {x|Ṙ ≥ 0} is empty. For comparison
purposes, Figure 3 depicts the ERA, the maximal level
set obtained with the numerator VN , which is strictly
contained in the ERA and the set {x|V̇N = 0}. This feature illustrates the conservativeness of the estimate obtained by computing the MLS of a given polynomial LF.

(29)

Set k ← k + 1; Rk ← R∗ with R∗ the optimizer of (29);
go to Step 1.

Notice that the Lyapunov function VN is a decision variable in every step of the above algorithms. This is not the
case if one imposes R = VN as mentioned in Remark 6.
Algorithm 3 guarantees a sequence of non-decreasing
ERAs and a function R(x) satisfying (25) is required
for its initialisation. This function can be taken as any
Lyapunov function satisfying the local stability of the
origin, for example a quadratic function of the form
V (x) = xT P x for some positive definite P if the linearized system matrix A is Hurwitz. The estimate obtained from running Algorithm 3 depends on the initial function R0 , hence running Algorithm 3 with different initialisations R0 may lead to better estimates, one

Example 2: The following three-dimensional system
from [20, Example 5] presents a limit cycle and an stable

9

We apply Algorithm 3 starting with R given by
a quadratic LF for the linearised system and set
deg(VN ) = 4. We then use the obtained LF, VN , as the
initial invariant set function R and apply Algorithm 3
again. The obtained ERA of degree two and degree
four are depicted in Figure 4. We were unable to find
multipliers p and m2 that satisfy (26) for the computed
VN and R, hence we could not construct a rational LF
certificate of the form (6).

2

1

x2

0

−1

−2
−2

−1

0

1

2

x1
Fig. 1. Dark blue line dark depicts the boundary of the ERA,
i.e. the set ∂E (R, γ). Level sets of the function V (x) are
also depicted. Trajectories depicted in green converge to the
origin.

1

x3

0

2

−1
1

x2

−1
1

0

0

x1

−1
−2
−2

−1

0

1

1

−1

x2

Fig. 4. The boundary of the ERA, ∂E (R, γ), of degree
deg(R) = 4 corresponds to the white surface while the
black surface corresponds to the largest ERA obtained with
deg(R) = 2. Trajectories are depicted in green (converging)
and red (diverging).

2

x1
Fig. 2. The dashed light blue lines depict the set of points
satisfying Ṙ(x) = h∇R(x), f (x)i = 0, the dashed black lines
depict the set V̇ (x) = h∇V (x), f (x)i = 0. The boundary of
the RA lies in the set where Ṙ(x) is negative.

Example 3: In the following example we compute a
piece-wise positively invariant set. Consider now system (1) with

2

1

x2

0

f (x) =

"

−x1 (1 − x1 x2 )
−x2

#

.

(30)

0

Despite the fact that only the origin is an equilibrium
point, its RA is not the whole of Rn nor is its boundary defined by a limit-cycle. The boundary of the RA is
given by {x ∈ Rn |x1 x2 = 2} (obtained from analytical
solution to Zubov’s equation [22, p.73]).

−1

−2
−2

−1

0

x1

1

2

Fig. 3. The boundary of the ERA, ∂E (R, γ), is depicted in
dark blue. The dashed black lines depict the set V̇N (x) = 0.
The solid black line depicts the MLS obtained with VN as
the LF.

We fix the shape of the positively invariant sets by fix1
ing RM as in (12), d = 2 with R1 = 200
(x21 − 2x1 x2 +
2
x2 ), R2 = 2x1 x2 , and compute VN , deg(VN ) = 6 to
VN
define a rational LF V = γ−R
. We formulate SOS conM
straints analogous to (25) for conditions in Theorem 3.
When solving the constraints we keep RM constant, i.e.
only the multipliers mi , p and the LF function VN are
updated while increasing γ. As a final step, we solve constraints with fixed γ, p, RM to obtain a rational LF of
which the level sets are depicted in Figure 5.

equilibrium at the origin:



f (x) = 


−x2
−x3
−0.915x1 + (1 − 0.915x21 )x2 − x3




.

10

[6] H.-D. Chiang, M. W. Hirsch, and F. F. Wu. Stability
regions of nonlinear autonomous dynamical systems. IEEE
Transactions on Automatic Control, 33(1):16–27, 1988.

5

x2

[7] P. Gahinet, P. Apkarian, and M. Chilali. Affine parameterdependent Lyapunov functions and real parametric
uncertainty. IEEE Transactions on Automatic Control,
41(3):436–442, Mar 1996.

0

[8] D. Henrion and J.-B. Lasserre. Inner approximations for
polynomial matrix inequalities and robust stability regions.
IEEE Transactions on Automatic Control, 57(6):1456–1467,
2012.
−5
−5

0

[9] H. K. Khalil. Nonlinear systems, volume 3. Prentice hall
Upper Saddle River, NJ, 2002.

5

x1

[10] J. B. Lasserre. Moments, Positive Polynomials and Their
Applications. Imperial College Press, London, 2009.

Fig. 5. The boundary of the ERA, ∂E (R, γ), is depicted in
dark blue. Trajectories obtained with (30) are depicted in
green (converging) and red (diverging).

6

[11] R. C. L. Oliveira and P. L. D. Peres.
Parameterdependent LMIs in robust analysis: characterization of
homogeneous polynomially parameter-dependent solutions
via LMI relaxations. IEEE Transactions on Automatic
Control, 52(7):1334–1340, 2007.

Conclusion

[12] A. Papachristodoulou, J. Anderson, G. Valmorbida,
S. Prajna, P. Seiler, and P.A. Parrilo. SOSTOOLS: Sum of
squares optimization toolbox for MATLAB V3.00, 2013.

In this paper we presented conditions for a positively
invariant set to be an ERA of the origin and for Lyapunov
certificates given by quotient of two functions where the
denominator characterizes the ERA. The main feature
of our results is that the positively invariant set defining
the ERA is not necessarily a level set of a Lyapunov
function. Provided a stronger condition is satisfied we
obtain a Lyapunov function interpretation of the ERA,
which for polynomial systems is a rational function. We
subsequently proposed an algorithm for the estimation
of the RA that guarantees the increment of the estimate
at each iteration.

[13] P. A. Parrilo.
Structured Semidefinite Programs and
Semialgebraic Geometry Methods in Robustness and
Optimization. PhD thesis, Caltech, Pasadena, CA, 2000.
[14] J. F. Sturm. Using SeDuMi 1.02, a MATLAB toolbox for
optimization over symmetric cones. Optimization Methods
and Software, 11–12:625–653, 1999.
[15] W. Tan and A. Packard. Stability region analysis using
polynomial and composite polynomial Lyapunov functions
and sum-of-squares programming. IEEE Transactions on
Automatic Control, 53(2):565–571, 2008.
[16] U. Topcu and A. Packard.
Local stability analysis
for uncertain nonlinear systems. IEEE Transactions on
Automatic Control, 54(5):1042–1047, 2009.

We then applied the algorithm to the class of polynomial vector fields and semi-algebraic sets for which the
steps are performed via the solution to Sum-of-Squares
programmes. The extension of the results to the class
of systems with parametric uncertainties was also presented. For this case, the positively invariant ERA does
not depend on the uncertain parameters while the associated Lyapunov function can.

[17] U. Topcu, A. Packard, P. Seiler, and G. J. Balas. Robust
region-of-attraction estimation.
IEEE Transactions on
Automatic Control, 55(1):137–142, 2010.
[18] G. Valmorbida and J. Anderson. Region of attraction analysis
via invariant sets. In American Control Conference (ACC),
2014, pages 3591–3596, 2014.
[19] L. Vandenberghe and S. Boyd. Semidefinite programming.
SIAM review, 38(1):49–95, 1996.

References

[20] A. Vannelli and M. Vidyasagar. Maximal Lyapunov functions
and domains of attraction for autonomous nonlinear systems.
Automatica, 21(1):69 – 80, 1985.

[1] G. Blekherman, P. A. Parrilo, and R. R. Thomas.
Semidefinite optimization and convex algebraic geometry.
SIAM, 2013.

[21] T. C. Wang, S. Lall, and T. Y. Chiou. Polynomial method
for PLL controller optimization. Sensors, 11(7):6575–6592,
2011.

[2] P.-A. Bliman. A convex approach to robust stability for linear
systems with uncertain scalar parameters. SIAM Journal on
Control and Optimization, 42(6):2016–2042, 2004.

[22] V. I. Zubov. Methods of AM Lyapunov and their Application.
Noordhoff Groningen, 1964.

[3] G. Chesi. Estimating the domain of attraction for uncertain
polynomial systems. Automatica, 40(11):1981–1986, 2004.
[4] G. Chesi. Domain of Attraction: Analysis and Control via
SOS Programming. Number 415 in Lecture Notes in Control
and Information Sciences. Springer, 2011.
[5] G. Chesi. Rational Lyapunov functions for estimating and
controlling the robust domain of attraction. Automatica,
49(4):1051 – 1057, 2013.

11

