Efficient Verification of CPA Lyapunov Functions
Sigurdur Freyr Hafstein
Science Institute, University of Iceland, Dunhagi 3, 107 Reykjav
´
ık, Iceland
Keywords:
Lyapunov Function, CPA Verification, Efficient Algorithms.
Abstract:
Lyapunov functions can be used to characterize the stability and basins of attraction for dynamical systems,
whose dynamics are defined by ordinary differential equations. Since the analytic generation of Lyapunov
functions for nonlinear systems is a formidable task, one often resorts to numerical methods. In this paper
we study the efficient verification of the conditions for a Lyapunov function using affine interpolation over
a triangulation; the values of the Lyapunov function candidate at the vertices of the triangulation can be
generated using various different formulas from converse theorems in the Lyapunov stability theory. Further,
we give an implementation in C++ and demonstrate its efficiency and applicability.
1 INTRODUCTION
In applications of dynamical systems in science and
engineering the stability of equilibria and other in-
variant sets is often a necessary requirement. In par-
ticular, control theory is concerned with the design-
ing of controllers and observers, such that the re-
sulting closed-loop systems are stable, e.g. have ex-
ponentially stable equilibria. The Lyapunov stabil-
ity theory is a much used tool in this regard and has
been intensively studied since its introduction by Lya-
punov in 1892 (Lyapunov, 1992); see e.g. the text-
books (Zubov, 1964; Yoshizawa, 1966; Hahn, 1967;
Sastry, 1999; Vidyasagar, 2002; Khalil, 2002). For a
physical system the obvious candidate for a Lyapunov
function is the system’s (free) energy and a dissipa-
tive physical system approaches a local minimum of
the energy.
The analytical generation of a Lyapunov func-
tion for a system is in general a formidable problem.
Therefore numerous numerical methods have been
developed, e.g. parameterizing rational (Vannelli and
Vidyasagar, 1985; Valmorbida and Anderson, 2017)
or polynomial (Parrilo, 2000; Chesi, 2011; Ander-
son and Papachristodoulou, 2015; Ratschan and She,
2010; Kamyar and Peet, 2015) Lyapunov functions;
for an overview of numerical methods see (Giesl and
Hafstein, 2015b).
In (Julian, 1999; Julian et al., 1999; Marin
´
osson,
2002a) linear programming was used to parameter-
ize continuous and piecewise affine (CPA) Lyapunov
functions. In this approach, a subset of the state space
is first triangulated, i.e. subdivided into simplices, and
then a number of constraints are derived for a given
nonlinear system, such that a feasible solution to the
resulting linear programming problem allows for the
parametrization of a CPA Lyapunov function for the
system. We refer to this method as the CPA algorithm
and the constraints of the linear programming prob-
lem as the CPA constraints.
In (Hafstein, 2004; Giesl and Hafstein, 2014) it
was proved that the CPA algorithm is always able
to compute a Lyapunov function for systems with an
exponentially stable equilibrium. Another approach
uses the CPA constraints, but instead of solving the
linear programming problem, a faster method to com-
pute values for the variables of the problem that likely
constitute a solution is used. Thus, the CPA con-
straints of the slow but rigorous CPA method are
combined with a faster but less rigorous method to
compute Lyapunov functions to deliver a fast and
rigorous method. The fast method has been based
on an integration- or sum formula for a Lyapunov
function from a converse theorem, see e.g. (Hafstein
et al., 2014a; Hafstein et al., 2014b; Bj
¨
ornsson et al.,
2014; Li et al., 2015; Hafstein et al., 2015; Bj
¨
ornsson
et al., 2015; Doban and Lazar, 2016; Doban, 2016;
Bj
¨
ornsson and Hafstein, 2017; Hafstein and Valfells,
2017; Hafstein and Valfells, 2019; Gudmundsson and
Hafstein, 2015; Hafstein, 2019), or on collocation us-
ing radial basis functions, see (Giesl and Hafstein,
2015a).
In this paper we present a fast verification of the
conditions of the linear programming problem from
120
Hafstein, S.
Efficient Verification of CPA Lyapunov Functions.
DOI: 10.5220/0011231700003271
In Proceedings of the 19th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2022), pages 120-129
ISBN: 978-989-758-585-2; ISSN: 2184-2809
Copyright
c
 2022 by SCITEPRESS – Science and Technology Publications, Lda. All rights reserved

Listing 1: Simple thread-parallelization in C++.
using bint = long long ; // b i nt = big inte g e r
void Pa r all e lFo r ( b i n t _beg , bi n t _end , fun c t ion < void ( b i nt ) > func , bint NrT h r ead s ){
for ( bint i = _beg ; i < _end ; i += N r T hre a d s ) {
vector < t hread > threads ( N r Th r e a ds );
for ( bint j = i; j < i + NrThr e a ds && j < _ end ; j + +) {
threa d s [ j % Nr T h rea d s ] = thread ( func , j );
}
for ( bint j = i ; j < i + NrThr e a ds && j < _ end ; j + +) {
threa d s [ j % Nr T h rea d s ]. join ();
}
}
}
void Pa r all e lFo r ( bint _end , functio n < v oid ( bint ) > func , bint NrT h r ea d s ) {
Par a lle l Fo r (0 , _end , func , N r T hre a d s );
}
(Marin
´
osson, 2002a) and its implementation in C++
using the Armadillo linear algebra library (Sanderson,
2010; Sanderson and Curtin, 2016). The implemen-
tation does not actually construct the linear program-
ming problem, but verifies the constraints on the fly,
and the procedure is very fast and memory efficient.
The paper is organized as follows. After present-
ing the necessary notation and prerequisites we dis-
cuss and develop CPA constraints that are particularly
well suited for an efficient verification in Section 2.
Then we discuss the implementation of the verifica-
tion in C++ in Section 3 and conclude the paper in
Section 4.
1.1 Notation and Prerequisites
We denote the set {1, 2,.. ., } by N and set N
0
:=
N ∪ {0}. For a vector x ∈ R
n
and p ≥ 1 we de-
fine the norms ∥x∥
p
= (
∑
n
i=1
|x
i
|
p
)
1/p
and ∥x∥
∞
=
max
i∈{1,...,n}
|x
i
|.
We utilize a bold-face font for (column) vectors,
e.g. x ∈ R
n×1
= R
n
, and we denote by e
1
,e
2
,. .. ,e
n
the standard orthonormal basis of R
n
. For a vector x
we write x
i
for its ith component (U
i
for the ith com-
ponent of U) and for a matrix A ∈ R
n×n
we write a
i j
for its (i, j)th element. Their transposes are denoted
x
T
and A
T
. The zero vector in R
n
is written 0 and 1 :=
(1,1, .. .,1)
T
∈ R
n
. For two vectors x, y ∈ R
n
the in-
equality x ≤ y is understood component-wise, i.e. x
i
≤
y
i
for i = 1,2, ... ,n. B = diag(a
1
,a
2
,. .. ,a
n
) defines a
matrix B ∈ R
n×n
with b
11
= a
1
,b
22
= a
2
,. .. ,b
nn
= a
n
and b
i j
= 0 for i ̸= j, i.e. B is a diagonal matrix.
Our C++ implementation makes heavy use of the
STL library and the Armadillo linear algebra library,
that is very well documented on its webpage http://
arma.sourceforge.net. Just a few comments: the STL
library delivers a routine to iterate through permuta-
tions called next_permutation; if one starts with the
identity permutation int v[i]=i one iterates through
all permutations. In Armadillo the most important
types, at least here, are vec for a column vector of
double and ivec for a column vector of integers; mat
and imat are the matrix versions. Since C++11 mul-
tithreading has been made very easy in C++, we use
the simple code in Listing 1 for multithreading.
2 CPA CONSTRAINTS
We consider systems, whose dynamics are given by
an ODE of the form
˙
x = f(x), f ∈ C
2
(R
n
;R
n
). (1)
Typically one assumes that the system in question has
an equilibrium w.l.o.g. at the origin, but that is not
necessary for our discussion. We assume there is a
Lyapunov function candidate V : D → R given, to-
gether with a triangulation of its domain D ⊂ R
n
. Our
intention is to verify where the continuous and affine
interpolation of the function V over the simplices of
the triangulation has a negative orbital derivative. Re-
call that a negative orbital derivative implies that the
function is decreasing along solution trajectories of
the system (1). For concretizing our ideas a few defi-
nitions are necessary.
Let C = {x
0
,x
1
,. .. ,x
n
} ⊂ R
n
be a set of affinely
independent vectors, i.e. the augmented vectors
(x
0
,1), (x
1
,1), .. .,(x
n
,1) ∈ R
n+1
are linearly inde-
pendent. The convex hull of the vectors in C, i.e. the
Efficient Verification of CPA Lyapunov Functions
121

set
coC :=
(
n
∑
k=0
λ
k
x
k
: x
k
∈ C, λ
k
∈ [0,1],
n
∑
k=0
λ
k
= 1
)
,
is called a proper n-simplex and the vectors in C are
said to be its vertices.
Let
{
S
ν
}
ν∈T
= T , T an index set, be a set of
proper n-simplices in R
n
, such that different simplices
S
ν
,S
µ
∈ T intersect in a common face or not at all
and such that the interior of the set D
T
=
S
ν∈T
S
ν
is a
simply connected set. The set T is said to be a shape-
regular triangulation in R
n
. We refer to the set
V
T
:= {x
i
: x
i
is a vertex of a simplex S
ν
∈ T }
as the vertex set of the triangulation T .
Let a shape-regular triangulation T =
{
S
ν
}
ν∈T
in
R
n
be given, together with the system (1) and the Lya-
punov functions candidate V : D → R with D = D
T
.
The following estimates are of essential importance
for the CPA algorithm, because they allow us to check
certain inequalities at the vertices V
T
of the simplices
in T to obtain estimates on the entire domain D
T
.
For a proper n-simplex S
ν
= co{x
0
,x
1
,. .. ,x
n
} in
T and the C
2
vector field f = ( f
1
, f
2
,. .. , f
n
)
T
define a
constant B
ν
r,s
such that
B
ν
r,s
≥ max
x∈S
ν
m=1,2,...,n




∂
2
f
m
∂x
r
∂x
s
(x)




. (2)
Further, for each (vertex) y of S
ν
define
C
ν
y,s
:= max
j=0,1,...,n
|e
s
•
(x
j
− y)|,
where
•
denotes the scalar-product, and set
E
y
ν,x
i
:=
1
2
n
∑
r,s=1
B
ν
r,s
|e
r
•
(x
i
− y)|(C
ν
y,s
+ |e
s
•
(x
i
− y)|)
(3)
for i = 0,1, .. .,n. The constants E
y
ν,x
i
are defined such
that for a fixed vector v ∈ R
n
we have that
v
•
f(x
i
) + E
y
ν,x
i
∥v∥
1
≤ 0 (4)
for i = 0, 1,. .., n implies v
•
f(x) ≤ 0 for all x ∈ coC.
This follows by the estimate (proved in
(Marin
´
osson, 2002b, Lemma 4.16))





f(x) −
n
∑
i=0
λ
i
f(x
i
)





∞
≤
n
∑
i=0
λ
i
E
y
ν,x
i
(5)
for all convex combinations x =
n
∑
i=0
λ
i
x
i
∈ S
ν
and
H
¨
older’s inequality:
v
•
f(x) =
n
∑
i=0
λ
i
v
•
f(x
i
) + v
•
"
f(x) −
n
∑
i=0
λ
i
f(x
i
)
#
≤
n
∑
i=0
λ
i
v
•
f(x
i
) + ∥v∥
1





f(x) −
n
∑
i=0
λ
i
f(x
i
)





∞
≤
n
∑
i=0
λ
i
v
•
f(x
i
) + ∥v∥
1
·
n
∑
i=0
λ
i
E
y
ν,x
i
=
n
∑
i=0
λ
i

v
•
f(x
i
) + ∥v∥
1
E
y
ν,x
i

≤ 0.
Note that in the condition (4) the B
ν
r,s
are just upper
bounds and that the vertex y of S
ν
for E
y
ν,x
i
is arbitrary
but fixed for i = 0,1,. .. ,n.
Given a triangulation T , a continuous and piece-
wise affine function, a so-called CPA function, can be
defined by fixing its values at V
T
.
Definition 2.1 (CPA function). For a shape-regular
triangulation T in R
n
we denote by CPA[T ] the set of
all continuous functions
V : D
T
→ R
that are affine on each simplex S
ν
∈ T . Hence, for
each S
ν
∈ T there exists a vector v
ν
∈ R
n
and a num-
ber a
ν
∈ R such that
V (x) = v
ν
•
x + a
ν
for all x ∈ S
ν
.
It is not difficult to see that V in the definition
above is completely determined by its values at the
vertices V
T
and with S
ν
= co{x
0
,x
1
,. .. ,x
n
},
w
ν
:=




V (x
1
) −V (x
0
)
V (x
2
) −V (x
0
)
.
.
.
V (x
n
) −V (x
0
)




∈ R
n
,
and
X
ν
:=





(x
1
− x
0
)
T
(x
2
− x
0
)
T
.
.
.
(x
n
− x
0
)
T





∈ R
n×n
,
the formula for v
ν
is
v
ν
= X
−1
ν
w
ν
. (6)
The system (1) together with the B
ν
r,s
is imple-
mented as in Listing 2, where we show it for the con-
crete example of the time reversed van der Pol system
˙
x = f(x) with f(x,y) =

−y
x + (x
2
− 1)y

. (7)
ICINCO 2022 - 19th International Conference on Informatics in Control, Automation and Robotics
122

Listing 2: Implementation of system (7).
extern c o n s t unsi g n ed int n = 2;
struct S ystem {
virtu a l ve c f( c o n s t v ec & x ) = 0; // the vec tor - field f in x ' = f ( x )
// B (r ,s , xl , xu ) >= max_ { i =1 ,.. , n } sup_ { xl <= x <= xu } | D_ { rs } f_i (x )|
virtu a l double B( int r , int s , const vec &xl , const vec & xu ) = 0;
};
struct I VDP : p u b l i c S y s t e m {
vec f ( const vec & x ) over r i d e {
vec fx ( n );
fx (0) = - x (1);
fx (1) = x (0) + ( pow ( x (0) ,2) - 1.0) * x ( 1);
return fx ;
}
double B( int r , int s , const vec &xl , const vec & xu ) overri d e {
if( r == 0 && s == 0){
return 2 .0* max ( abs ( xl (1)) , abs ( xu (1)));
}
else if ( r ==0 && s = =1 || r ==1 && s ==0){
return 2 .0* max ( abs ( xl (0)) , abs ( xu ( 0 )));
}
else {
return 0 .0;
}
}
};
A concrete system inherits the virtual class System
and delivers implementations for the member func-
tions System.f and System.B. The latter gives an up-
per bound as in (2), however not for x ∈ S
ν
but for x
in the cube {x ∈ R
n
: xl ≤ x ≤ xu} . Note that r,s = 0
corresponds to the variable x(0) = x and r, s = 1 to the
variable x(1) = y and that the only non-zero second-
order derivatives of the components of f = ( f
1
, f
2
)
T
in
(7) are
∂
2
f
2
∂x
2
(x,y) = 2y and
∂
2
f
2
∂x∂y
(x,y) = 2x.
For an efficient implementation of the verification
it is advantageous to use particularly simple and reg-
ular triangulations, both in terms of speed but more
importantly in terms of size. We will use a stair case
triangulation, which is the so-called standard triangu-
lation T
std
, see e.g. (Albertsson et al., 2020) and for
extensions (Giesl and Hafstein, 2021b; Giesl and Haf-
stein, 2021a), but restricted to {x ∈ R
n
: 0 ≤ x ≤ U}
for a vector U ∈ N
n
. This triangulation is then scaled
along the axes and translated to the area of interest.
Definition 2.2 (Stair case Triangulation). The stair
case triangulation T
U
sc
for a vector U ∈ N
n
is a tri-
angulation T
U
sc
= {
e
S
ν
}
ν∈T
with indices ν = (z, σ),
for all z ∈ N
0
fulfilling 0 ≤ z ≤ U − 1 and all per-
mutations σ of {1,2, ... ,n}. The vertices of each
e
S
ν
= (
e
x
ν
0
,
e
x
ν
1
,. .. ,
e
x
ν
n
) are given by
e
x
ν
k
= z +
k
∑
j=1
e
σ( j)
= z + u
σ
k
(8)
where
u
σ
k
=
k
∑
j=1
e
σ( j)
.
Given a star case triangulation T
U
sc
and vectors
ℓ
ℓ
ℓ,u ∈ R
n
, ℓ
ℓ
ℓ < u, we define P
i
:= (u
i
− ℓ
i
)/U
i
for i =
1,2, .. .,n and the matrix P := diag(P
1
,P
2
,. .. ,P
n
) ∈
R
n×n
. The triangulation T
U
sc,ℓ
ℓ
ℓ,u
= {S
ν
}
ν∈T
is now de-
fined by mapping the simplices
e
S
ν
of T
U
sc
with the
mapping x 7→ Px + ℓ
ℓ
ℓ, i.e. S
ν
in T
U
sc,ℓ
ℓ
ℓ,u
is the image
of the simplex
e
S
ν
. The vertices of the simplex S
ν
are
given by
x
ν
k
:= P
e
x
ν
k
+ ℓ
ℓ
ℓ
for k = 0,1,.. ., n.
Note that since T
U
sc
is a triangulation of {x ∈
R
n
: 0 ≤ x ≤ U} we have that T
U
sc,ℓ
ℓ
ℓ,u
is a triangulation
of {x ∈ R
n
: ℓ
ℓ
ℓ ≤ x ≤ u}.
Using the triangulation T
U
sc,ℓ
ℓ
ℓ,u
the constraints (4)
for ν = (z,σ), S
ν
= co{x
ν
0
,x
ν
1
,. .. ,x
ν
n
} and with v
ν
:=
Efficient Verification of CPA Lyapunov Functions
123

v and E
ν
i
:= E
x
ν
0
ν,x
i
can be written in the form
0 ≥
n
∑
j=1
V (x
ν
j
) −V (x
ν
j−1
)
e
σ( j)
•
(x
ν
j
− x
ν
j−1
)
f
σ( j)
(x
ν
i
) (9)
+ E
ν
i
n
∑
j=1





V (x
ν
j
) −V (x
ν
j−1
)
e
σ( j)
•
(x
ν
j
− x
ν
j−1
)





,
for i = 0, 1,. .., n, where
E
ν
i
=
1
2
n
∑
r,s=1
B
ν
r,s
A
σ
r,i
(A
σ
s,i
+ A
σ
s,n
) (10)
and
A
σ
k,i
:= |e
k
•
(x
ν
i
− x
ν
0
)|
for k = 1,2,.. ., n and i = 0, 1,.. ., n.
Note that A
σ
k,i
= P
k
if k ∈ {σ(1),σ(2),. .. ,σ(i)}
and A
σ
k,i
= 0 otherwise and that A
σ
s,n
= C
ν
x
ν
0
,s
for all
ν = (z, σ). Further,
e
σ( j)
•
(x
ν
j
− x
ν
j−1
) = e
σ( j)
•
Pe
σ( j)
= P
σ( j)
,
which simplifies the formulas even further, and the
A
σ
k,i
only depend on σ in ν = (z,σ). Hence, the
A
σ
k,i
, for all permutations σ of {1,2, ... ,n}, all k =
1,2, .. .,n and all i = 0,1,.. ., n, can be computed once
and for all before the verification of (9) for all the sim-
plices S
ν
∈ T
U
sc,ℓ
ℓ
ℓ,u
.
3 EFFICIENT CPA
CONSTRAINTS VERIFICATION
For the implementation of the simplicial complexes
we use two grids. The first one is of type
struct ZGrid and models the vertices of T
U
sc
and the
second is of type struct xGrid and models the ver-
tices of T
U
sc,ℓ
ℓ
ℓ,u
. Their declarations and the defini-
tion of xGrid are given in Listing 3. The imple-
mentation of ZGrid was discussed in (Hafstein, 2013)
(as class Grid). Its main purpose here is to en-
able linear indexing of the grid points, i.e. G.V2I(i)
for i=0,1,. .., G.NrPoints()-1 iterates through all the
grid points of G. The implementation of xGrid is es-
sentially trivial, just note that / and % are component-
wise division and multiplication of vectors in Ar-
madillo, respectively; similar to ./ and .* in Mat-
lab. The correspondence between the variables and
the mathematical symbols is xL = ℓ
ℓ
ℓ, xU = u, iU = U
and hv = (P
1
,P
2
,. .. ,P
n
).
The verification is then implemented in the func-
tion VerifyLya, see Listing 4. The arguments of the
function are System *psys, a pointer to a concrete sys-
tem that implements psys->f and psys->B, see Listing
2. The argument const vec &V is a vector containing
the values of the function V at the vertices of the trian-
gulation T
U
sc,ℓ
ℓ
ℓ,u
modelled by the argument xGrid xG.
More exactly, V(i) is the value of V at the vertex
xG.I2vec(i) (index to vector). If V fails the condi-
tion (9) for some permutation σ in ν = (z, σ), then the
index of the vertex x = Pz + ℓ
ℓ
ℓ is written to the vec-
tor vector<bint> &Failed, that is also an argument to
the function. We do not keep track of the permutation
σ, and thus the exact simplices, where the condition
fails. This is faster and requires less memory and usu-
ally the information about the cube P(z + [0.1]
n
) + ℓ
ℓ
ℓ
containing a simplex where the Lyapunov function
candidate V fails to have a negative orbital derivative
is detailed enough. The vector Failed is sorted before
the function returns.
A few comments on the implementation:
• The code is more optimized for size than for
speed. The reason is that the computation of the
values in V typically takes much more time than
the verification anyways. Storing the values of
psys->f at all vertices of T
U
sc,ℓ
ℓ
ℓ,u
would speed up
the verification, but the storing needs n times the
memory needed to store the values in V, where n is
the dimension of the system, and would introduce
a new limiting factor.
• It is important when verifying the condition (9)
to not use orbder + errbound > 0.0 to check
if the condition fails on a simplex, but to
instead use !(orbder + errbound <= 0.0). If
orbder+errbound is not NaN (not a number) then
both are equivalent. However, since the methods
used to compute the values in V, i.e. the values of
V at the vertices, sometimes produce NaN, there
is an important difference. A comparison with
NaN always produces false, even NaN == NaN is
false, but we want orbder + errbound > 0.0 to
deliver true if orbder+errbound is NaN and this is
achieved with !(orbder + errbound <= 0.0).
• We use the infamous goto to break out of nested
loops when we have identified a cube
P(z + [0.1]
n
) + ℓ
ℓ
ℓ ⊃ S
ν
containing a simplex S
ν
where the condition (9)
fails; the alternatives are simply more messy.
• Finally, it might seem at first glance that the
indices of the vertices where the condition (9)
fails are already sorted, and thus an extra sorting
at the end of the function is unnecessary. A closer
inspection, however, reveals that they are already
sorted with respect to the order of the vertices in
ZGrid ROG(xG.G.L, xG.G.U - ones<ivec>(n)),
over which the algorithm iterates to obtain the zs
ICINCO 2022 - 19th International Conference on Informatics in Control, Automation and Robotics
124

Listing 3: Code for the grids.
// for L , U in Z ˆ n all z in Z ˆ n fulf i lli n g L <= z <= U componen t - w i s e
struct Z G r i d {
ivec L , U ;
bint End I n dex ;
ZGrid ( const ivec & _L , const ivec & _U );
bint V2 I ( ivec ) const ; // ive c to index
ivec I2 V ( bint ) const ; // index to ivec
bint NrP o i nts ( void );
};
struct x G r i d {
ZGrid G ;
vec xL , xU , hv ;
vec I2vec ( bint I n d e x ) c o n s t { // i n d e x to vec
return xL + c onv_to < vec >:: from ( G . I2V ( I n d e x )) % hv ;
}
xGrid ( const vec & _xL , const vec & _xU , const ivec & iU )
: xL ( _xL ) , xU ( _xU ) , G ( zeros < ivec >( n ) , iU ){
hv = ( xU - xL ) / conv_to < vec >:: from ( iU );
}
};
in ν = (z,σ), but not with respect to the order of
the vertices in xG modeling T
U
sc,ℓ
ℓ
ℓ,u
.
To give an idea of the performance of the al-
gorithm we computed on a Threadripper 3990X a
Lyapunov function candidate using the RBF method
(Giesl, 2007). We used 12,480 collocation points
and the computation of the parameters for the func-
tion took 3.5 sec, i.e. writing and solving a system
of 12,480 linear equations in 12,480 unknowns. The
values of V are then computed at 1,002,001 vertices
in 4.9 sec. Finally, the negativity of the orbital deriva-
tive in 2,000,000 simplices was verified in 0.026 sec;
see Figure 1 for the results.
Figure 1: Lyapunov function candidate for system (7) and
the area of its domain where the function is not decreasing
along solution trajectories (red).
4 CONCLUSIONS
We presented an algorithm to verify the negativity
of the orbital derivative of CPA Lyapunov function
candidates interpolated over shape-regular triangula-
tions. By carefully choosing a regular triangulation
the algorithm is very fast and memory efficient. Sub-
level sets of Lyapunov functions are positively invari-
ant sets and this can be efficiently computed for CPA
Lyapunov functions (Giesl et al., 2020). We gave an
implementation in C++ of the algorithm and demon-
strated its applicability and effectiveness for an exam-
ple.
ACKNOWLEDGEMENT
The research done for this paper was partially sup-
ported by the Icelandic Research Fund in the grant
228725-051, Linear Programming as an Alternative
to LMIs for Stability Analysis of Switched Systems,
which is gratefully acknowledged.
Efficient Verification of CPA Lyapunov Functions
125

Listing 4: Code for the verification.
bint fa c t ori a l ( int k ) { b i n t f =1; for ( i nt i = 1; i <= k ; i ++) f *= i ; return f; }
void Ve r i fyL y a ( System * psys , const vec &V , xGrid xG , vector < bint > & Faile d ,
int NrTh r e ads ) {
vector < int > iv ( n );
for ( int i = 0; i < n ; i ++) iv [ i ] = i ;
vector < v ector < int > > Sigma ;
Sigma . rese r v e ( fa c t ori a l ( n ));
do {
Sigma . pu s h _ba c k ( iv );
} w h i l e ( ne x t_ p er m ut a ti o n (& iv [0] , & iv [ n ]));
vector < mat > A ;
A. rese r v e ( fa c t ori a l ( n ));
for ( auto p s i g m a = Sigma . b e g i n (); psigma != S i g m a . end (); psigma ++) {
mat rA (n , n + 1);
vec sx = zeros < vec >( n );
for ( int r = 0; r < n ; r ++) {
rA(r , 0) = 0.0;
double v a l r A ;
for ( int i = 0; i < n ; i ++) {
if ((* ps i g m a )[ i ] == r) v a l r A = xG . hv ( r );
rA(r , i + 1) = valrA ;
}
}
A. pu s h _ba c k ( rA );
}
wal l _ cl o c k c l o c k ; clock . tic ( );
ZGrid ROG ( xG . G .L , xG . G . U - ones < ivec >( n ));
cout << " Ver i f yin g " << RO G . N r Poin t s () * f act o r ial ( n ) < < " simp l i ces " < < endl ;
vector < v ector < int > > fa i l ed P a rt ;
fai l e dP a r t . re s i z e ( NrT h r ead s );
functio n < v oid ( bint ) > parf o r = [ & ]( b i nt Th ) {
bint lowB = ( ROG . N rPo i n t s () * Th ) / Nr T h rea d s ;
bint up B = ( ROG . N rPoi n t s () * ( Th + 1)) / NrThre a d s ;
ivec z , zx ;
vec gradV ( n );
vec f xi ( n ) , fx0 ( n) , x0 ( n ) , xi ( n );
mat B ( n , n );
double V0 , Vold , Vnew ;
int r , s , i , j , k ;
bint I0 , Inew ;
double orbder , grbound , errbo u n d ;
for ( bint pnr = low B ; pnr < upB ; pnr ++) { // z
z = ROG . I2V ( pnr );
I0 = xG . G . V2I ( z );
x0 = xG . I2vec ( I0 );
fx0 = psys - > f ( x0 );
V0 = V ( I0 );
for ( r = 0; r < n ; r + +) {
for ( s = 0; s < r ; s + +) {
B(s , r ) = B (r , s ) = psys - >B (r , s , x0 , x0 + xG . hv );
}
B(r , r ) = psys - >B (r , r , x0 , x0 + xG . hv );
}
ICINCO 2022 - 19th International Conference on Informatics in Control, Automation and Robotics
126

Listing 4: Code for the verification (cont.).
for ( k = 0; k < Sigma . s i ze (); k + +) { // sigma
vector < int > psigma = Sigma [ k ];
Vold = V0 ;
zx = z ;
for ( int j = 0; j < n ; j ++) { // nabla V
zx( psigma [j ]) += 1;
Inew = xG . G . V2I ( zx );
Vnew = V ( Inew );
gradV ( psigma [ j ]) = ( Vnew - Vold ) / xG . hv ( psigm a [ j ] );
Vold = Vnew ;
}
// i =0
orbder = dot ( gradV , fx 0 );
if (!( or b d e r <= 0.0)) { // e r r boun d =0.0
goto CU B E FA I L ED ;
}
grbou n d = no r m ( gradV , 1) ;
zx = z ;
for ( int i = 0; i < n ; i ++) {
zx( psigma [i ]) += 1;
Inew = xG . G . V2I ( zx );
xi = xG . I2vec ( Inew );
fxi = psys - > f ( xi );
orbder = dot ( gradV , fx i );
errb o u nd = 0.0;
for ( r = 0; r < n ; r + +) {
for ( s = 0; s < n ; s + +) {
errb o u nd += B ( r , s ) * A [ k ]( r , i + 1) *
(A [k ]( s , i + 1) + A [ k ](s , n ));
}
}
errb o u nd *= 0 .5 * g r b o u nd ;
if (!( or b d e r + errbo u n d <= 0 . 0 ) ) {
goto CU B E FA I L ED ;
}
}
}
cont i n ue ;
CUB E F AI L E D :
fai l e dP a r t [ Th ]. p u sh_ b a ck ( I0 );
}
};
Par a lle l Fo r ( N r T h reads , parfor , NrThr e a ds );
Failed . clear ();
bint NrF a i led = 0;
for ( int i = 0; i < NrTh r e ads ; i ++) {
NrFa i l ed += fail e dPa r t [ i ]. size ();
}
Failed . r e s erve ( N r Fail e d );
for ( int i = 0; i < NrTh r e ads ; i ++) {
Failed . i n s e r t ( Failed . end () , f a il e d Par t [ i ]. b e g i n () , fai l e dP a r t [ i ]. end ());
}
sort ( F a i l ed . begin () , F a i l e d . e nd ());
}
Efficient Verification of CPA Lyapunov Functions
127

REFERENCES
Albertsson, S., Giesl, P., Gudmundsson, S., and Hafstein,
S. (2020). Simplicial complex with approximate ro-
tational symmetry: A general class of simplicial com-
plexes. J. Comput. Appl. Math., 363:413–425.
Anderson, J. and Papachristodoulou, A. (2015). Advances
in computational Lyapunov analysis using sum-of-
squares programming. Discrete Contin. Dyn. Syst. Ser.
B, 20(8):2361–2381.
Bj
¨
ornsson, J., Giesl, P., Hafstein, S., Kellett, C., and Li,
H. (2014). Computation of continuous and piecewise
affine Lyapunov functions by numerical approxima-
tions of the Massera construction. In Proceedings
of the CDC, 53rd IEEE Conference on Decision and
Control, pages 5506–5511, Los Angeles (CA), USA.
Bj
¨
ornsson, J., Giesl, P., Hafstein, S., Kellett, C., and Li, H.
(2015). Computation of Lyapunov functions for sys-
tems with multiple attractors. Discrete Contin. Dyn.
Syst. Ser. A, 35(9):4019–4039.
Bj
¨
ornsson, J. and Hafstein, S. (2017). Efficient Lya-
punov function computation for systems with multiple
exponentially stable equilibria. Procedia Computer
Science, 108:655–664. Proceedings of the Interna-
tional Conference on Computational Science (ICCS),
Zurich, Switzerland, 2017.
Chesi, G. (2011). Domain of Attraction: Analysis and Con-
trol via SOS Programming. Lecture Notes in Control
and Information Sciences, vol. 415, Springer.
Doban, A. (2016). Stability domains computation and sta-
bilization of nonlinear systems: implications for bio-
logical systems. PhD thesis: Eindhoven University of
Technology.
Doban, A. and Lazar, M. (2016). Computation of Lyapunov
functions for nonlinear differential equations via a
Yoshizawa-type construction. IFAC-PapersOnLine,
49(18):29 – 34.
Giesl, P. (2007). Construction of Global Lyapunov Func-
tions Using Radial Basis Functions. Lecture Notes in
Math. 1904, Springer.
Giesl, P. and Hafstein, S. (2014). Revised CPA method to
compute Lyapunov functions for nonlinear systems. J.
Math. Anal. Appl., 410:292–306.
Giesl, P. and Hafstein, S. (2015a). Computation and veri-
fication of Lyapunov functions. SIAM J. Appl. Dyn.
Syst., 14(4):1663–1698.
Giesl, P. and Hafstein, S. (2015b). Review of computa-
tional methods for Lyapunov functions. Discrete Con-
tin. Dyn. Syst. Ser. B, 20(8):2291–2331.
Giesl, P. and Hafstein, S. (2021a). System specific trian-
gulations for the construction of CPA Lyapunov func-
tions. Discrete Contin. Dyn. Syst. Ser. B, 26(12):6027–
6046.
Giesl, P. and Hafstein, S. (2021b). Uniformly regular tri-
angulations for parameterizing Lyapunov functions.
In Proceedings of the 18th International Conference
on Informatics in Control, Automation and Robotics
(ICINCO), pages 549–557.
Giesl, P., Osborne, C., and Hafstein, S. (2020). Au-
tomatic determination of connected sublevel sets of
CPA Lyapunov functions. SIAM J. Appl. Dyn. Syst.,
19(2):1029–1056.
Gudmundsson, S. and Hafstein, S. (2015). Lyapunov func-
tion verification: MATLAB implementation. In Pro-
ceedings of the 1st Conference on Modelling, Identifi-
cation and Control of Nonlinear Systems (MICNON),
number 0235, pages 806–811.
Hafstein, S. (2004). A constructive converse Lyapunov the-
orem on exponential stability. Discrete Contin. Dyn.
Syst. Ser. A, 10(3):657–678.
Hafstein, S. (2013). Implementation of simplicial com-
plexes for CPA functions in C++11 using the Ar-
madillo linear algebra library. In Proceedings of the
3rd International Conference on Simulation and Mod-
eling Methodologies, Technologies and Applications
(SIMULTECH), pages 49–57, Reykjavik, Iceland.
Hafstein, S. (2019). Computational Science - ICCS 2019:
19th International Conference, Faro, Portugal, June
12-14, 2019, Proceedings, Part V, chapter Numerical
Analysis Project in ODEs for Undergraduate Students,
pages 412–434. Springer.
Hafstein, S., Kellett, C., and Li, H. (2014a). Computation
of Lyapunov functions for discrete-time systems using
the Yoshizawa construction. In Proceedings of 53rd
IEEE Conference on Decision and Control (CDC).
Hafstein, S., Kellett, C., and Li, H. (2014b). Continu-
ous and piecewise affine Lyapunov functions using the
Yoshizawa construction. In Proceedings of the 2014
American Control Conference (ACC), pages 548–553
(no. 0170), Portland (OR), USA.
Hafstein, S., Kellett, C., and Li, H. (2015). Computing con-
tinuous and piecewise affine Lyapunov functions for
nonlinear systems. J. Comput. Dyn, 2(2):227 – 246.
Hafstein, S. and Valfells, A. (2017). Study of dynamical
systems by fast numerical computation of Lyapunov
functions. In Proceedings of the 14th International
Conference on Dynamical Systems: Theory and Ap-
plications (DSTA), volume Mathematical and Numer-
ical Aspects of Dynamical System Analysis, pages
220–240.
Hafstein, S. and Valfells, A. (2019). Efficient computation
of Lyapunov functions for nonlinear systems by in-
tegrating numerical solutions. Nonlinear Dynamics,
97(3):1895–1910.
Hahn, W. (1967). Stability of Motion. Springer, Berlin.
Julian, P. (1999). A High Level Canonical Piecewise Lin-
ear Representation: Theory and Applications. PhD
thesis: Universidad Nacional del Sur, Bahia Blanca,
Argentina.
Julian, P., Guivant, J., and Desages, A. (1999). A
parametrization of piecewise linear Lyapunov func-
tions via linear programming. Int. J. Control, 72(7-
8):702–715.
Kamyar, R. and Peet, M. (2015). Polynomial optimization
with applications to stability analysis and control – an
alternative to sum of squares. Discrete Contin. Dyn.
Syst. Ser. B, 20(8):2383–2417.
Khalil, H. (2002). Nonlinear Systems. Pearson, 3. edition.
Li, H., Hafstein, S., and Kellett, C. (2015). Computation of
continuous and piecewise affine Lyapunov functions
ICINCO 2022 - 19th International Conference on Informatics in Control, Automation and Robotics
128

for discrete-time systems. J. Difference Equ. Appl.,
21(6):486–511.
Lyapunov, A. M. (1992). The general problem of the sta-
bility of motion. Internat. J. Control, 55(3):521–790.
Translated by A. T. Fuller from
´
Edouard Davaux’s
French translation (1907) of the 1892 Russian orig-
inal, With an editorial (historical introduction) by
Fuller, a biography of Lyapunov by V. I. Smirnov, and
the bibliography of Lyapunov’s works collected by J.
F. Barrett, Lyapunov centenary issue.
Marin
´
osson, S. (2002a). Lyapunov function construction
for ordinary differential equations with linear pro-
gramming. Dynamical Systems: An International
Journal, 17:137–150.
Marin
´
osson, S. (2002b). Stability Analysis of Nonlin-
ear Systems with Linear Programming: A Lyapunov
Functions Based Approach. PhD thesis: Gerhard-
Mercator-University Duisburg, Duisburg, Germany.
Parrilo, P. (2000). Structured Semidefinite Programs and
Semialgebraic Geometry Methods in Robustness and
Optimization. PhD thesis: California Institute of
Technology Pasadena, California.
Ratschan, S. and She, Z. (2010). Providing a basin of attrac-
tion to a target region of polynomial systems by com-
putation of Lyapunov-like functions. SIAM J. Control
Optim., 48(7):4377–4394.
Sanderson, C. (2010). Armadillo: An open source C++
linear algebra library for fast prototyping and com-
putationally intensive experiments. Technical report,
NICTA.
Sanderson, C. and Curtin, R. (2016). Armadillo: a template-
based C++ library for linear algebra. J. Open Source
Softw, 1(2):26.
Sastry, S. (1999). Nonlinear Systems: Analysis, Stability,
and Control. Springer.
Valmorbida, G. and Anderson, J. (2017). Region of attrac-
tion estimation using invariant sets and rational Lya-
punov functions. Automatica, 75:37–45.
Vannelli, A. and Vidyasagar, M. (1985). Maximal Lya-
punov functions and domains of attraction for au-
tonomous nonlinear systems. Automatica, 21(1):69–
80.
Vidyasagar, M. (2002). Nonlinear System Analysis. Clas-
sics in Applied Mathematics. SIAM, 2. edition.
Yoshizawa, T. (1966). Stability theory by Liapunov’s sec-
ond method. Publications of the Mathematical Society
of Japan, No. 9. The Mathematical Society of Japan,
Tokyo.
Zubov, V. I. (1964). Methods of A. M. Lyapunov and their
application. Translation prepared under the auspices
of the United States Atomic Energy Commission;
edited by Leo F. Boron. P. Noordhoff Ltd, Groningen.
Efficient Verification of CPA Lyapunov Functions
129