Efficient Implementation of Piecewise Quadratic Lyapunov Function
Computations for Switched Linear Systems
∗
Stefania Andersen
1 a
, Sigurdur Hafstein
1 b
, Juan Javier Palacios Roman
2 c
and Sebastiaan J.A.M. van den Eijnden
2 d
1
Faculty of Physical Sciences, University of Iceland, Dunhagi 5, 107 Reykjavik, Iceland
2
Eindhoven University of Technology, Department of Mechanical Engineering,
Groene Loper 3, 5612 AE Eindhoven, Netherlands
Keywords:
Common Lyapunov Function, Switched Linear Systems, Linear Programming, C++, MATLAB.
Abstract:
We describe a linear programming (LP) problem to parameterize continuous and piecewise quadratic (CPQ)
Lyapunov functions for switched linear systems. We discuss some algorithms and data-structures for its im-
plementation in C++ and compare the computational efficiency of our implementation to an analogous imple-
mentation in MATLAB.
1 INTRODUCTION
Switched systems play an important role in modelling
in science and engineering (Davrazos and Kous-
soulas, 2001; Liberzon, 2003; Shorten et al., 2007;
Sun and Ge, 2011). In control theory the stability
of an equilibrium point is commonly of central im-
portance and is most conveniently dealt with using
the Lyapunov stability theory of dynamical systems
(Hahn, 1967; Sastry, 1999; Khalil, 2002; Vidyasagar,
2002). It is well known that the existence of a Lya-
punov function for a switched system, which is a
common Lyapunov function for all the subsystems,
is equivalent to its stability and suitable classes of po-
tential Lyapunov functions is a thoroughly researched
subject (Dayawansa and Martin, 1999; Goebel et al.,
2006; Mason et al., 2006; Shorten et al., 2007; Ma-
son et al., 2022). In this paper we consider the al-
gorithm from (Palacios Roman et al., 2024) to pa-
rameterize continuous and piecewise quadratic (CPQ)
Lyapunov functions for switched linear systems along
with its efficient implementation in the programming
language C++. The class of CPQ Lyapunov func-
tions have been considered in (Johansson and Rantzer,
a
https://orcid.org/0000-0001-6747-775X
b
https://orcid.org/0000-0003-0073-2765
c
https://orcid.org/0009-0004-3145-1119
d
https://orcid.org/0000-0003-1273-109X
∗
This work was supported in part by the Icelandic Re-
search Fund under Grant 228725-051.
1998), and have recently been successfully used to
study the stability of, hybrid integrator-gain systems
(HIGS), see e.g. (van den Eijnden et al., 2020; Dee-
nen et al., 2021; van den Eijnden et al., 2022).
The main contribution of this work is the effi-
cient computation of CPQ Lyapunov functions for
switched linear systems based on a linear program,
where the classical constraints for CPQ Lyapunov
functions, which are often expressed in terms of lin-
ear matrix inequalities, are formulated as a linear pro-
gramming (LP) problem.
This paper is organized as follows: In Section 2
we describe how we write the LP problem in our im-
plementation, before we discuss in Section 3 trian-
gulations, CPQ functions, and how we parameterize
them. Section 4 deals with the linear constraints we
use to compute a suitable parameterization for a CPQ
Lyapunov function for a given switched linear system.
Then, in Section 5, we give some details on the imple-
mentation of the linear constraints and compare the
numerical efficiency of our implementation in C++
to a corresponding implementation in MATLAB. We
conclude the paper in Section 6.
2 WRITING OUR LP PROBLEM
An LP feasibility problem can be characterized using
a matrix A ∈ R
r×c
, a vector b ∈ R
r
, and an enum vec-
tor s, s
i
∈ {LE,GE,EQ}, having r elements. A feasible
Andersen, S., Hafstein, S., Roman, J. and van den Eijnden, S.
Efficient Implementation of Piecewise Quadratic Lyapunov Function Computations for Switched Linear Systems.
DOI: 10.5220/0012927500003822
Paper published under CC license (CC BY-NC-ND 4.0)
In Proceedings of the 21st International Conference on Informatics in Control, Automation and Robotics (ICINCO 2024) - Volume 1, pages 277-284
ISBN: 978-989-758-717-7; ISSN: 2184-2809
Proceedings Copyright © 2024 by SCITEPRESS – Science and Technology Publications, Lda.
277

solution to the LP problem is a vector x ∈ R
c
that sat-
isfies
[Ax]
i
[s
i
] b
i
, for all i = 1,2,...,r,
where [Ax]
i
is the i-th component of the vector Ax.
This means that
[Ax]
i
≤ b
i
, if s
i
is LE,
[Ax]
i
= b
i
, if s
i
is EQ,
and [Ax]
i
≥ b
i
, if s
i
is GE.
Since the matrix A is usually sparse, we implement
it using two vectors, rn and cn, of row and column
indices, respectively, and one vector of values val.
The vector b is implemented by the vector b and the
vector of symbols s by the enum TypeCon {LE,GE,EQ}
vector con.
The vectors rn, cn and val will all have the same
number of elements, namely the number of non-zero
elements of A, and to set the (i, j)-th element of A to
some non-zero value x, i.e. a
i, j
:=x, we write
rn[k]=i; cn[k]=j; val[k]=x;
A constraint of the LP problem is written in a row of
the matrix A. A column of A corresponds to a vari-
able. More exactly, the columns of A contain the co-
efficients with which the corresponing variable in the
vector Variables appears in the constrains of the LP
problem. The Variables vector is described below.
2.1 Storing the Variables
To implement a variable of the LP problem we use
an Armadillo (Sanderson, 2010) integer vector ivec.
For example, the variable ϕ
ν
k,ℓ
is ivec {'P',nu,k,l},
where nu, k and l are some numbers of type bint (big
integer, i.e. long long). All the variables of the prob-
lem are stored together in the sorted vector Variables.
In order to sort the Variables vector we need to
define an order relation. We first let the variable’s
length dictate the relation. If the variables have the
same length we compare the first non-equal elements
or last elements, which ever comes first. This is im-
plemented with the following binary function:
1 b ool v ar Cmp ( cons t ivec &x , c onst ...
iv ec & y) {
2 if ( x. si ze () < y. size ()) {
3 re tu rn tr ue ;
4 }
5 else if (x. s ize () > y. size () ) {
6 re tu rn fa ls e ;
7 }
8 in t i, len = x. size () ;
9 fo r (i =0; i < len -1 && x ( i) == ...
y(i ) ; i ++) {};
10 r eturn x(i ) < y ( i) ;
11 }
We also implement a function, VarID, that ob-
tains the index of a variable in Variables. That is,
if var=Variables(i) then VarID(var) returns i. If the
variable is not in Variables we return the impossible
value −1. The function is as follows:
1 b int V arID ( c on st ivec & v) {
2 auto fo und = ...
equal _r an ge ( V ar ia bl es . begi n () , ...
Variabl es . end () , v, varCmp ) ;
3 if ( f ound . first == fou nd . second ){
4 re tu rn -1;
5 }
6 else {
7 re tu rn fo un d . first - ...
Variabl es . begi n ();
8 }
9 }
2.2 Construction of the Matrix A
Since each row of the matrix A corresponds to a con-
straint, it is preferable to write this matrix in a row-
by-row manner. We implement the function, new_aij,
that writes a new nonzero element of A or modi-
fies a nonzero value, and a function, close_constr,
that “closes” the constraint, such that the next call to
new_aij will write an element in the next row. We
store the current index of rn, cn and val in the bint
variable Index, and the current row number in the int
variable ConstrNr. The corresponding functions are
as follows:
1 v oid n ew _a ij ( cons t ivec ...
& variabl e , do ub le valu e ){
2 bint VID = VarI D ( v ariable ) ;
3 a ssert ( VID != -1) ;
4 fo r ( bi nt i = Inde x - 1; i >= 0 ...
&& rn[ i ] == ConstrNr ; i - -) {
5 if ( cn [i] == VI D ) {
6 val[i ] += valu e ;
7 return ;
8 }
9 rn . push_ ba ck ( C on st rN r ) ;
10 cn . push_ ba ck ( VID );
11 val . p us h_back ( value ) ;
12 In dex ++;
13 }
14 }
1 v oid c los e_ co nst r ( T yp eC on Type , ...
do ub le _b ) {
2 co n . pus h_ ba ck ( Type );
3 b. push_ ba ck ( _b ) ;
4 C on st rN r ++;
5 }
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
278

2.3 Minimal LP Example
As a minimal example for our format, consider the
constraints
3x − y ≤ 5 (1)
x + 2y = 7.
We need the following variables:
1 us ing b int = long lo ng ;
2 us ing namespa ce std;
3 us ing namespa ce ar ma ;
4 b int I ndex ;
5 int Con st rN r ;
6 vector <ivec > V ar ia bl es
7 vector <bint > rn , cn ;
8 vector < double > val , b;
9 e num T yp eC on { LE , EQ , GE };
10 vector < Typ eCo n > con ;
The constraints can be implemented as:
1 vector <ivec > V ar ia bl es ={ ...
iv ec {' x' }, ive c {' y' }};
2 rn [ 0]=0; cn [0]=0; val [ 0] =3;
3 rn [ 1]=0; cn [1]=1; val [1]= -1;
4 con [0]= LE; b [0 ]= 5;
5 rn [ 2]=1; cn [2]=0; val [ 2] =1;
6 rn [ 3]=1; cn [3]=1; val [ 3] =2;
7 con [1]= EQ; b [1 ]= 7;
or alternatively with new_aij and close_constr as:
1 vector <ivec > V ar ia bl es ={ ...
iv ec {' x' }, ive c {' y' }};
2 new_aij ( i vec {'x' } ,3) ;
3 new_aij ( i vec {'y' } , -1) ;
4 c lo se _cons tr ( LE ,5);
5 new_aij ( i vec {'x' } ,1) ;
6 new_aij ( i vec {'y' } ,2) ;
7 c lo se _cons tr ( EQ ,7);
Note that 'x' and 'y' are interpreted as their ASCII
dec equivalents, that is 120 and 121 respectively.
The purpose of the structure discussed above is
to act as a parser between the problem formulation
and the linear solver we use to solve the LP problems,
Gurobi (Gurobi Optimization, LLC, 2023). The same
example using the Gurobi interface for C would be
given by
1 int va rs = 2;
2 int constrs = 2;
3 si ze _t vb eg [] = {0, 2};
4 int vl en [] = {2 , 2};
5 int vi nd [] = {0 , 1, 0 , 1};
6 do ub le vv al [] = {3.0 , 1.0 , -1.0 , ...
2. 0};
7 c har s ense [] = ...
{ G RB _L ES S_ EQUAL , G RB _E QU AL };
8 do ub le r hs [] = {5. 0 , 7.0} ;
9 er ror = G RB Xlo ad mo del ( env , ...
& model , " e xa mp le " , vars , ...
constrs , 1, 0.0 , NULL , sen se , ...
rhs , vbeg , vlen , vind , vval , ...
NULL , NULL , NULL , NULL , NU LL ) ;
The syntax of this interface has two main draw-
backs; it is rather opaque and one can easily overwrite
a previously declared nonzero element of the matrix
A accidentally. If, for example, we obtained the 3x
in (1) via a sum (x + 2x) and wrote
1 int vl en []={3 , 2};
2 int vi nd []={0 , 0 , 1, 0, 1};
3 do ub le vv al []={1.0 , 2.0 , 1.0 , ...
-1.0 , 2.0};
we would overwrite the first instance 1.0·x with 2.0·x
rather than adding them together. A parser is easily
implemented:
1 // v ar_size = V ar ia bl es . size ()
2 // rnc nv al _si ze = rn . s ize ()
3 for ( in t v =0 ,j =0; v < v ar _s iz e ; v ++) {
4 vbeg [v]=j;
5 fo r ( int k =0 ; k < rn cnv al _s ize ;k++) {
6 if ( cn[ k ] == v){
7 vv al [ j ]= val [k];
8 vi nd [ j ]= rn[k ];
9 j++;
10 }
11 }
12 }
3 CPQ FUNCTIONS
In order to parameterize a CPQ function, we use a
triangulation of the domain of the function. We will
first discuss the specifics of a suitable triangulation T .
Then, a prelude on the parameterization of continuous
and piecewise affine (CPA) functions will be given,
from which the parameterization of CPQ functions
using the triangulation T will follow naturally. Note
that our parameterization of CPQ Lyapunov functions
largely follows (Johansson, 1999).
3.1 Triangulation
A triangulation T with vertices {v
1
,v
2
,..., v
p
} ⊂ R
n
is a subdivision of a subset of R
n
into simplices S. A
simplex S
ν
is defined as
S
ν
:=co{x
ν
0
,x
ν
1
,..., x
ν
n
}
=
(
x ∈ R
n
: x =
n
∑
i=0
λ
i
x
ν
i
,λ
i
≥ 0,
n
∑
i=0
λ
i
= 1
)
,
Efficient Implementation of Piecewise Quadratic Lyapunov Function Computations for Switched Linear Systems
279

S
0
S
1
S
3
S
2
S
6
S
7
S
5
S
4
v
6
v
9
v
8
v
7
v
4
v
1
v
2
v
3
v
5
Figure 1: Triangulation T
1
with simplices and vertices in-
dexed.
where x
ν
i
= v
g
ν
(i)
for i = 0,1,. .. ,n and g
ν
:
{0,1,... ,n} −→ {1, 2,. .., p} is an index-function. For
our purposes, we require the triangulation T to
be shape-regular, i.e., every two different simplices
S
ν
,S
µ
∈ T either intersect in a common lower-
dimensional face or do not intersect at all. Further, we
demand that x
ν
0
= 0
0
0 and that the vectors x
ν
1
,x
ν
2
,..., x
ν
n
are linearly independent for all S
ν
∈ T . Finally, the
set theoretic union of all S
ν
∈ T , denoted D
T
, must
be a neighbourhood of the origin of R
n
.
An efficient implementation of a triangulation
that satisfies these requirements is the triangular fan
(T
std
K,fan
)
F
discussed in (Hafstein, 2019), from here on
denoted simply by T
K
, where a formula for all x
ν
i
is
given. The vertices {v
1
,v
2
,..., v
p
} of the triangula-
tion T
K
are
{0
0
0} ∪

K
∥
z
∥
2
z : z ∈ Z
n
,
∥
z
∥
∞
= K

,
where the scaling parameter K ∈ N determines the
fineness of the triangulation around zero. The num-
ber of simplices in the triangulation T
K
is given by
the formula 2
n
· K
n−1
· n!. The triangulation T
1
in R
2
is depicted in Figure 1.
3.2 Prelude on CPA Functions
To parameterize CPA functions on the whole space
R
n
, we associate a cone C
ν
to each simplex S
ν
∈ T ,
defined as
C
ν
:
=cone{x
ν
1
,x
ν
2
,..., x
ν
n
}
=
(
x ∈ R
n
: x =
n
∑
i=1
λ
i
x
ν
i
,λ
i
≥ 0
)
.
We refer to the vector λ
λ
λ = (λ
i
)
i=1,2,...,n
∈ R
n
+
as the
cone coordinates of x in C
ν
. Additionally, we define
the matrix X
ν
:= [x
ν
1
,..., x
ν
n
] ∈ R
n×n
. Note that for
this definition we must assume that the order of the
vertices of each S
ν
is fixed.
For a CPA function W : R
n
→ R defined using the
triangulation T we have for every x ∈ C
ν
that W (x) =
w
T
ν
x, where
w
T
ν
=
h
W (x
ν
1
),W (x
ν
2
),...,W(x
ν
n
)
i
X
−1
ν
.
Just note that for all x ∈ C
ν
we have x = X
ν
λ
λ
λ and thus
W (x) = w
T
ν
x (2)
=
h
W (x
ν
1
),W (x
ν
2
),...,W(x
ν
n
)
i
X
−1
ν
X
ν
λ
λ
λ
=
n
∑
i=1
λ
i
W (x
ν
i
).
To implement the computation of W using LP, it is
advantageous to define the vector of variables Φ ∈ R
p
,
where ϕ
j
= W (v
j
), j = 1,2, ..., p. Then W (x
ν
i
) =
ϕ
g
ν
(i)
for every vertex x
ν
i
= v
j
, i = 0, 1,..., n, of some
S
ν
∈ T . Next, we define the vector
Ψ
ν
:=

ϕ
g
ν
(1)
,ϕ
g
ν
(2)
,..., ϕ
g
ν
(n)

T
for every S
ν
∈ T . Then
W (x) = (Ψ
ν
)
T
λ
λ
λ
by (2), where λ
λ
λ ∈ R
n
+
are the cone coordinates of
x ∈ C
ν
.
That W is well-defined and continuous now fol-
lows easily from the fact that the vertices x
ν
1
,..., x
ν
n
of every S
ν
∈ T are linearly independent and that T
is shape-regular: Since T is shape-regular, C
ν
∩ C
µ
=
cone{y
1
,y
2
,..., y
j
}, where j < n and the y
i
are the
common nonzero vertices of S
ν
and S
µ
, i.e., y
i
=
x
ν
k
i
= x
µ
ℓ
i
for some k
i
,ℓ
i
∈ {1,..., n} and g
ν
(k
i
) =
g
µ
(ℓ
i
) for all i = 1, ..., j. Then, for all x ∈ C
ν
∩ C
µ
x =
j
∑
i=1
λ
i
y
i
=
j
∑
i=1
λ
i
x
ν
k
i
=
j
∑
i=1
λ
i
x
µ
ℓ
i
,
for some unique λ
λ
λ, because of the linear indepen-
dence of the vertices. Hence
(Ψ
ν
)
T
λ
λ
λ =
j
∑
i=1
λ
i
ϕ
g
ν
(k
i
)
=
j
∑
i=1
λ
i
ϕ
g
µ
(ℓ
i
)
= (Ψ
µ
)
T
λ
λ
λ
and W is well-defined on C
ν
∩ C
µ
and continuous,
since it is the restriction of the continuous functions
x 7→ w
T
ν
x and x 7→ w
T
µ
x to C
ν
∩ C
µ
. Hence, the vector
ϕ
ϕ
ϕ ∈ R
p
and the functions g
ν
connect the different for-
mulas W (x) = w
T
ν
x, S
ν
∈ T , such that the resulting
function
W (x) = w
T
ν
x if x ∈ C
ν
is well-defined and continuous.
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
280

3.3 Representation of CPQ Functions
For representing a CPQ function V : R
n
→ R, i.e., V is
continuous and V (x) = x
T
P
ν
x for a symmetric matrix
P
ν
∈ R
n×n
on each S
ν
∈ T , we can proceed similarly
as in the CPA case. However, now we need a p × p
matrix Φ ∈ R
p×p
of variables, such that for a simplex
S
ν
∈ T we have
P
ν
= (X
−1
ν
)
T
Ψ
ν
X
−1
ν
,
where the (k,ℓ)-th entry of the matrix Ψ
ν
∈ R
n×n
is
equal to ϕ
g
ν
(k),g
ν
(ℓ)
, which is the (g
ν
(k), g
ν
(ℓ))-th en-
try of Φ.
Then, for an x =
∑
n
i=1
λ
i
x
ν
i
∈ C
ν
, i.e., x = X
ν
λ
λ
λ
with λ
λ
λ ∈ R
n
+
, we have
V (x) = x
T
P
ν
x
= λ
λ
λ
T
X
T
ν
(X
−1
ν
)
T
Ψ
ν
X
−1
ν
X
ν
λ
λ
λ
= λ
λ
λ
T
Ψ
ν
λ
λ
λ
=
n
∑
k,ℓ=1
λ
k
λ
ℓ
ϕ
g
ν
(k),g
ν
(ℓ)
.
(3)
That V is well-defined and continuous follows simi-
larly as for the CPA case: For all x ∈ C
ν
∩ C
µ
x =
j
∑
i=1
λ
i
y
i
=
j
∑
i=1
λ
i
x
ν
k
i
=
j
∑
i=1
λ
i
x
µ
ℓ
i
,
for some k
i
,ℓ
i
∈ {1, .. .,n} and a unique λ
λ
λ. Since
x
T
P
ν
x =
j
∑
r=1
λ
r
(x
ν
k
r
)
T
(X
−1
ν
)
T
Ψ
ν
X
−1
ν
j
∑
s=1
λ
s
x
ν
k
s
=
j
∑
r=1
λ
r
e
T
k
r
Ψ
ν
j
∑
s=1
λ
s
e
k
s
=
j
∑
r,s=1
λ
r
λ
s
ϕ
g
ν
(k
r
),g
ν
(k
s
)
,
where e
i
is the i-th unit vecotor of R
n
. Likewise,
x
T
P
µ
x =
∑
j
r,s=1
λ
r
λ
s
ϕ
g
µ
(ℓ
r
),g
µ
(ℓ
s
)
and it can be con-
cluded that
x
T
P
ν
x = x
T
P
µ
x
because g
ν
(k
i
) = g
µ
(ℓ
i
) for i = 1,2,. .. , j. Thus, V is
well-defined and continuous.
Remark 1. Note that g
ν
(0), i.e. the index of the zero
vertex v
k
= 0
0
0, is never needed in the formulas above.
4 COMPUTING A CPQ
LYAPUNOV FUNCTION
Consider a switched linear system with arbitrary
switching
˙
x(t) = A(t)x(t), A(t) ∈ A
:
= {A
1
,A
2
,..., A
N
},
where A : R
+
−→ A is an arbitrary right-continuous
piecewise constant mapping and N ∈ N is finite. As
discussed in the Introduction, the origin is asymptoti-
cally stable for the system, if and only if there exists a
Lyapunov function for the system, i.e., a locally Lip-
schitz function V : R
n
→ R that satisfies
V (x) > 0 ∀x ∈ R
n
\{0
0
0}, (4)
V (0
0
0) = 0,
⟨
∇V (x), A
i
x
⟩
< 0 ∀x ∈ R
n
\{0
0
0},
∀i ∈ {1,2, ...,N},
where
⟨
∇V (x), A
i
x
⟩
denotes the inner product of the
vectors ∇V (x) and A
i
x. Strictly speaking ∇V (x) is
the Clarke’s subdifferential (Clarke, 1990), which in
our CPQ case means that for x ∈ C
ν
∩ C
µ
we need

∇(x
T
P
ν
x),A
i
x

< 0 and

∇(x
T
P
µ
x),A
i
x

< 0.
Going into details would go beyond the scope of this
paper and we refer the interested reader to the litera-
ture (Clarke et al., 1998; Baier et al., 2012).
If V is CPQ and parameterized as in (3), then con-
ditions (4) are equivalent to
x
T
P
ν
x > 0 ∀x ∈ C
ν
\{0
0
0}
x
T
Q
ν
i
x < 0 ∀x ∈ C
ν
\{0
0
0}
∀i ∈ {1,2, ...,N}
for all S
ν
∈ T , where Q
ν
i
:
= A
T
i
P
ν
+ P
ν
A
i
. We will
refer to these conditions as P
ν
and Q
ν
i
being posi-
tive definite (p.d.) and negative definite (n.d.), respec-
tively, on the cone C
ν
.
4.1 Reexpression of Constraints
Now it is desired to reexpress the condition on Q
ν
i
into
one that can be verified with LP. Note that for x ∈ C
ν
we can write x = X
ν
λ
λ
λ, where λ
λ
λ ∈ R
n
+
are the cone
coordinates of x. Thus we can define a matrix
B
ν
i
= X
T
ν
Q
ν
i
X
ν
= X
T
ν
(A
T
i
(X
−1
ν
)
T
Ψ
ν
X
−1
ν
+ (X
−1
ν
)
T
Ψ
ν
X
−1
ν
A
i
)X
ν
=

Ψ
ν
X
−1
ν
A
i
X
ν

T
+ Ψ
ν
X
−1
ν
A
i
X
ν
= (Ψ
ν
ˆ
A
ν
i
)
T
+ Ψ
ν
ˆ
A
ν
i
,
where
ˆ
A
ν
i
= X
−1
ν
A
i
X
ν
. Then we define a matrix C
ν
i
such that c
ν,i
k,k
= b
ν,i
k,k
and c
ν,i
k,ℓ
= max{0,b
ν,i
k,ℓ
} for all
k, ℓ ∈ {1,2,. .. ,n}, where b
ν,i
k,ℓ
and c
ν,i
k,ℓ
are the (k, ℓ)-
th entries of B
ν
i
and C
ν
i
, respectively.
Consider the next proposition.
Proposition 1. If the sum
∑
n
ℓ=1
c
ν,i
k,ℓ
< 0 for all k ∈
{1,2,... ,n}, then Q
ν
i
is n.d. on the cone C
ν
.
Efficient Implementation of Piecewise Quadratic Lyapunov Function Computations for Switched Linear Systems
281

Sketch of proof. We will provide a sketch of the proof
here; a detailed proof will be presented in an upcom-
ing paper (Palacios Roman et al., 2024).
Because the off-diagonal elements of C
ν
i
are posi-
tive, the conditions
∑
n
ℓ=1
c
ν,i
k,ℓ
< 0 force the matrix C
ν
i
to be a diagonally dominant and n.d. matrix. Because
of how C
ν
i
is defined from B
ν
i
, λ
λ
λ
T
C
ν
i
λ
λ
λ ≥ λ
λ
λ
T
B
ν
i
λ
λ
λ for all
λ
λ
λ ∈ R
n
+
. Thus, B
ν
i
is n.d. on the cone R
n
+
, from which
it follows that Q
ν
i
is n.d. on C
ν
.
Finally, we can implement the max function in
c
ν,i
k,ℓ
= max{0,b
ν,i
k,ℓ
} with the constraints c
ν,i
k,ℓ
≥ 0 and
c
ν,i
k,ℓ
≥ b
ν,i
k,ℓ
. Then, the conditions on C
ν
i
become
c
ν,i
k,ℓ
≥
n
∑
r=1

ψ
ℓ,r
ˆa
ν,i
r,k
+ ψ
k,r
ˆa
ν,i
r,ℓ

c
ν,i
k,ℓ
≥ 0
for all k, ℓ ∈ {1,2, ...n} and k ̸= ℓ and
2
n
∑
r=1
ψ
k,r
ˆa
ν,i
r,k
+
n
∑
ℓ=1
ℓ̸=k
c
k,ℓ
< 0.
for all k ∈ {1, 2,. ..,n}.
4.2 Constraints of the LP Problem
Equipped with the above results, the final constraints
for the LP problem can now be summarized. On
every simplex S
ν
:= co{0
0
0,x
ν
1
,x
ν
2
,..., x
ν
n
}, we have
X
ν
= [x
ν
1
,x
ν
2
,..., x
ν
n
] and index-function g
ν
such that
x
ν
i
= v
g
ν
(i)
. The constraints are then given by:
1. V (v
i
) > 0 for all v
i
∈ T .
2. For each system define
ˆ
A
ν
i
:
= X
−1
ν
A
i
X
ν
and C
ν
i
such that
c
ν,i
k,k
= 2
n
∑
r=1
ϕ
g
ν
(r),g
ν
(k)
ˆa
ν,i
r,k
, (5)
for all k ∈ {1, .. .,n} and
c
ν,i
k,ℓ
≥
n
∑
r=1

ϕ
g
ν
(ℓ),g
ν
(r)
ˆa
ν,i
r,k
+ ϕ
g
ν
(k),g
ν
(r)
ˆa
ν,i
r,ℓ

,(6)
c
ν,i
k,ℓ
≥ 0
for all k,ℓ ∈ {1,. ..,n} and k ̸= ℓ. Then require for
each S
ν
∈ T that
n
∑
ℓ=1
c
ν,i
k,ℓ
< 0 ∀k ∈ {1, 2,. ..,n}.
The problem variables are the elements of the matri-
ces Φ and C
ν
i
.
S
0
v
5
= x
0
1
v
8
= x
0
2
v
4
= x
0
0
Figure 2: Simplex S
0
with its vertices indexed. For
example, with simNo=0 and vertNo=1 the functions
VertexNr(simNo,VertNo) and Vertex(simNo,VertNo)
return the int 5 and Armadillo vector [1,0]
T
respectively.
Remark 2. Note that V(v
i
) > 0 for all nonzero ver-
tices v
i
of T does not ensure that V (x) > 0 for all
x ∈ R
n
\ {0
0
0}. This can however be checked a posteri-
ori to solving the LP problem by checking whether
there exists a solution to the linear matrix inequal-
ity (LMI)
P
ν
− (X
−1
ν
)
T
U
ν
X
−1
ν
≻ 0,
where U
ν
∈ R
n×n
+
is the problem variable and A ≻ B
means that A − B is a symmetric and positive definite
matrix. This a posteriori check was proposed in (Jo-
hansson, 1999, Proposition 4.4).
5 IMPLEMENTATION AND
COMPARISON
In our implementation the object corresponding to
the simplicial complex T
K
provides the functions
VertexNr(simNo,VertNo), to obtain the vertex num-
ber analogous to the indexing function g
ν
(i), and
Vertex(simNo,VertNo) to obtain the vertex vector
analogous to x
ν
i
(see Figure 2). With these functions
we can define our variables and construct our con-
straints discussed in Section 4.2.
Note that since Φ is a sparse symmetric matrix, we
only need to define variables for the upper triangular
half. Further, we only need ϕ
k,ℓ
if v
k
and v
ℓ
are non-
zero vertices of the same simplex and then we define
var = ivec {'P', min(k,l), max(k,l)}. For exam-
ple for the simplex in Figure 2 we only need ϕ
5,8
, ϕ
5,5
and ϕ
8,8
. To avoid unnecessary checks and multiple
definitions of variables, we first define the variables
above the diagonal in this way, and subsequently the
diagonal ϕ
k,k
, k = 1,.. ., p. Note that where the index
for the zero vertex can be skipped, as it is not used in
the constraints.
Furthermore, each C
ν
i
is a symmetric n × n matrix
so we only need the upper triangular half of the ele-
ments. In fact, we only need the elements above the
diagonal since the diagonal is equal to B
ν
i
, so rather
than using c
ν,i
k,k
in the constraint we can directly use
the right hand side of (5).
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
282

In this way we can define each variable, push into
the Variables vector and then sort it.
All that is left to complete the setup of the LP
problem is defining the constraints. This is were
the benefit of our structure defined in Section 2 is
most apparent since there is no need to keep track of
whether a variable appears twice in a sum, as in (6), or
not. We can simply go through each constraint using
new_aij to add the coefficient of each summand and
close_constr before moving to the next constraint.
As an alternative, we could use a parser such as
YALMIP (L
¨
ofberg, 2004) in MATLAB (The Math-
Works Inc., 2023) to solve for Φ. This would admit-
tedly entail a less convoluted set-up process since we
could create matrix variables and don’t need to imple-
ment the max function with an intermediate variable
like C
ν
i
. However, since MATLAB is an interpreted
language it doesn’t benefit from the efficiency asso-
ciated with compiled languages, in our case the effi-
ciency of the GNU Compiler g++.
For comparison consider the arbitrarily switched
linear system discussed in (Andersen et al., 2023;
Polanski, 1997; Pyatnitskii, 1971; Brockett, 1966)
with subsystems given by
A
1
=

0 1
−0.01 −2

and A
2
=

0 1
−11.7 −2

.
We can compute a CPQ Lyapunov function with a
scaling factor K ≥ 4 and run the computations twenty
times for each K. As can be seen in Figure 3, the av-
erage time for each scaling factor using the C++ im-
plementation is roughly 1/200th of the time needed
when using MATLAB. We used the SDPT3 (Toh
et al., 1999) solver in the MATLAB implementation
and Gurobi version 11 for the C++ implementation.
Additional MATLAB tooling to speed up the compu-
tations like the MATLAB Coder can not be utilized to
minimize computation time without major changes to
YALMIP.
6 CONCLUSIONS
We presented the implementation of a method to com-
pute continuous and piecewise quadratic (CPQ) Lya-
punov functions for switched linear systems in C++;
such a Lyapunov function is a common Lyapunov
function for all the linear subsystems. We compared
our implementation to an analogous implementation
in MATLAB and showed that the C++ implementa-
tion is more than 100 times faster for an example sys-
tem from the literature. We expect this difference in
computational speed to be representative for the gen-
eral case, as the setup of the linear programming (LP)
5 10 15 20
10
-2
10
0
matlab
C++
Figure 3: Comparison of the computational time for param-
eterizing Lyapunov functions with our method in MATLAB
and C++; run on i9900K (8 cores) under Linux Mint. Aver-
age computation time of twenty tests as a function of trian-
gulation scaling factor K in T
K
. Note that the C++ code is
ca. 200 times faster than the MATLAB code.
problem used to parameterize the CPQ Lyapunov
function is quite involved and an interpreted computer
language like MATLAB is at a great disadvantage in
comparison to a compiled language like C++.
REFERENCES
Andersen, S., Giesl, P., and Hafstein, S. (2023). Common
Lyapunov functions for switched linear systems: Lin-
ear programming-based approach. IEEE Control Sys-
tems Letters, 7:901–906.
Baier, R., Gr
¨
une, L., and Hafstein, S. (2012). Linear pro-
gramming based Lyapunov function computation for
differential inclusions. Discrete Contin. Dyn. Syst. Ser.
B, 17(1):33–56.
Brockett, R. (1966). The status of stability theory for de-
terministic systems. IEEE Trans. Automat. Control,
11(3):596–606.
Clarke, F. (1990). Optimization and Nonsmooth Analysis.
Classics in Applied Mathematics. SIAM.
Clarke, F., Ledyaev, Y., and Stern, R. (1998). Asymptotic
stability and smooth Lyapunov functions. J. Differen-
tial Equations, 149:69–114.
Davrazos, G. and Koussoulas, N. (2001). A review of sta-
bility results for switched and hybrid systems. In Pro-
ceedings of 9th Mediterranean Conference on Control
and Automation, Dubrovnik, Croatia.
Dayawansa, W. and Martin, C. (1999). A converse Lya-
punov theorem for a class of dynamical systems which
undergo switching. IEEE Trans. Automat. Control,
(44):751–760.
Deenen, D., Sharif, B., van den Eijnden, S., Nijmeijer, H.,
Heemels, M., and Heertjes, M. (2021). Projection-
based integrators for improved motion control: For-
malization, well-posedness and stability of hybrid
integrator-gain systems. Automatica, 133:109830.
Goebel, R., Hu, T., and Teel, A. (2006). Current Trends
in Nonlinear Systems and Control. Systems and Con-
trol: Foundations & Applications, chapter Dual Ma-
trix Inequalities in Stability and Performance Analy-
sis of Linear Differential/Difference Inclusions, pages
103–122. Birkhauser.
Efficient Implementation of Piecewise Quadratic Lyapunov Function Computations for Switched Linear Systems
283

Gurobi Optimization, LLC (2023). Gurobi Optimizer Ref-
erence Manual.
Hafstein, S. (2019). Simulation and Modeling Method-
ologies, Technologies and Applications, volume 873
of Advances in Intelligent Systems and Computing,
chapter Fast Algorithms for Computing Continuous
Piecewise Affine Lyapunov Functions, pages 274–
299. Springer.
Hahn, W. (1967). Stability of Motion. Springer, Berlin.
Johansson, M. (1999). Piecewise Linear Control Systems.
PhD thesis: Lund University, Sweden.
Johansson, M. and Rantzer, A. (1998). Computation of
piecewise quadratic Lyapunov functions for hybrid
systems. IEEE Trans. Automat. Control, 43(4):555–
559.
Khalil, H. (2002). Nonlinear Systems. Pearson, 3. edition.
Liberzon, D. (2003). Switching in systems and control.
Systems & Control: Foundations & Applications.
Birkh
¨
auser.
L
¨
ofberg, J. (2004). YALMIP: A toolbox for modeling and
optimization in MATLAB. In In Proceedings of the
CACSD Conference, Taipei, Taiwan.
Mason, P., Boscain, U., and Chitour, Y. (2006). Common
polynomial Lyapunov functions for linear switched
systems. SIAM J. Control Optim., 45(1):226–245.
Mason, P., Chitour, T., and Sigalotti, M. (2022). On univer-
sal classes of Lyapunov functions for linear switched
systems. arXiv:2208.09179.
Palacios Roman, J., Hafstein, S., Giesl, P., van den Eijnden,
S., Andersen, S., and Heemels, M. (2024). Construct-
ing continuous piecewise quadratic Lyapunov func-
tions with linear programming. [Manuscript in prepa-
ration].
Polanski, A. (1997). Lyapunov functions construction by
linear programming. IEEE Trans. Automat. Control,
42:1113–1116.
Pyatnitskii, Y. (1971). Absolute stability criterion for the
2nd order nonlinear controlled systems with one non-
linear nonstationary element. Autom. Remote Control,
32(1):1–11.
Sanderson, C. (2010.). Armadillo: An open source C++
linear algebra library for fast prototyping and com-
putationally intensive experiments. Technical report,
NICTA.
Sastry, S. (1999). Nonlinear Systems: Analysis, Stability,
and Control. Springer.
Shorten, R., Wirth, F., Mason, O., Wulff, K., and King, C.
(2007). Stability criteria for switched and hybrid sys-
tems. SIAM Review, 49(4):545–592.
Sun, Z. and Ge, S. (2011). Stability Theory of Switched Dy-
namical Systems. Communications and Control Engi-
neering. Springer.
The MathWorks Inc. (2023). MATLAB version: 23.2.0
(r2023b).
Toh, K. C., Todd, M. J., and T
¨
ut
¨
unc
¨
u, R. H. (1999). SDPT3
- a matlab software package for semidefinite program-
ming, version 1.3. Optimization Methods and Soft-
ware, 11(1-4):545–581.
van den Eijnden, S., Heemels, M., Nijmeijer, H., and Heert-
jes, M. (2022). Stability and performance analysis of
hybrid integrator–gain systems: A linear matrix in-
equality approach. Nonlinear Analysis: Hybrid Sys-
tems, 45:101192.
van den Eijnden, S., Heertjes, M., Heemels, M., and Ni-
jmeijer, H. (2020). Hybrid integrator-gain systems:
A remedy for overshoot limitations in linear control?
IEEE Control Systems Letters, PP:1–1.
Vidyasagar, M. (2002). Nonlinear System Analysis. Clas-
sics in Applied Mathematics. SIAM, 2. edition.
ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics
284