Implementation of Simplicial Complexes for CPA Functions in C++11

using the Armadillo Linear Algebra Library

Sigurdur Freyr Hafstein

School of Science and Engineering, Reykjavik University, Menntavegur 1, 101 Reykjavik, Iceland

Keywords:

CPA Function, Lyapunov Function, Piecewise Linear, Nonlinear System, Triangulation, Simplicial Complex,

C++11, Armadillo Linear Algebra Library.

Abstract:

Continuous, piecewise afﬁne (CPA) functions can be algorithmically parameterized to deliver Lyapunov func-

tions for compact invariant sets. We discuss ﬂexible structures and algorithms to manipulate CPA functions

for these purposes and discuss their implementation in C++11 using the Armadillo linear algebra library. Es-

pecially, we discuss some of the new language features in C++11 that lead to simpler and more readable

code. The implementation was developed in the freeware Visual Studio Express 2012 for Windows Desktop

(VS2012). Apart from a detailed description and code examples for the construction and manipulation of

the simplicial complex that serves as a basis for CPA functions, this contribution includes some discussion

on practical implementation details when using VS2012, C++11, and the linking to and use of the excellent

Armadillo linear algebra library. Thus, some parts of this paper, especially Section 3, might be useful not only

for those interested in the implementation of the simplicial complex for computing CPA Lyapunov functions,

but also for those generally interested in using the free Armadillo library for computations in VS2012.

1 INTRODUCTION

Lyapunov functions are a fundamental concept in the

study of dynamical systems. Their central role in

studies of the stability behavior of dynamical systems

is well known. Their construction is, however, difﬁ-

cult in the general case, i.e. for nonlinear systems.

Several methods to numerically compute Lya-

punov functions for nonlinear systems have been sug-

gested. To name a few, in (Johansson and Rantzer,

1997) a construction method for piecewise quadratic

Lyapunov functions for piecewise afﬁne autonomous

systems is suggested. In (Eghbal, Pariz, and Karim-

pour, 2012) the computation of piecewise quadratic

Lyapunov functions for planar piecewise afﬁne sys-

tems is formulated as linear matrix inequalities. In

(Johansen, 2000) linear programming is used to pa-

rameterize Lyapunov functions for autonomous non-

linear systems. In (Rezaiee-Pajand and Moghad-

dasie, 2012) a different collocation method using two

classes of basis functions is suggested. In (Giesl,

2007) radial basis functions are used to solve numer-

ically a generalized Zubov equation. In (Peet and

Papachristodoulou, 2010) the numerical construction

of Lyapunov functions that are presentable as sum of

squares of polynomials is considered. The Lyapunov

functions are computed by means of convexoptimiza-

tion.

One method that has been studied in some de-

tail recently, uses linear programming to parameterize

CPA Lyapunov functions in compact neighbourhoods

of exponentially stable equilibria. This approach was

ﬁrst followed in (Julian, Guivant, and Desages, 1999)

and was enhanced in (Marinosson, 2002a and 2002b)

to compute true Lyapunov functions, rather than ap-

proximations requiring a posteriori analysis to deter-

mine their quality. In (Hafstein, 2004 and 2005) it was

proved that when an arbitrary small hypercube around

the equilibrium is excluded from the domain of the to

be computed CPA Lyapunov function, the computa-

tion would always succeed. The domain of the com-

puted CPA Lyapunov function is otherwise only lim-

ited to any compact subset of the equilibrium’s region

of attraction.

In (Giesl and Hafstein, 2012 and 2013) the neces-

sity of excluding an arbitrary small hypercube around

the equilibrium was removed, at the expense of need-

ing a more reﬁned simplicial complex than in pervi-

ous works. In this paper we will discuss the imple-

mentation of this novel simplicial complex that pos-

sesses a simplicial fan at the equilibrium.

The term simplicial fan seems natural, for math-

Freyr Hafstein S..

Implementation of Simplicial Complexes for CPA Functions in C++11 using the Armadillo Linear Algebra Library.

DOI: 10.5220/0004423400490057

In Proceedings of the 3rd International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH-2013),

pages 49-57

ISBN: 978-989-8565-69-3

 2013 SCITEPRESS (Science and Technology Publications, Lda.)

Figure 1: The simplicial complex T

std

N,K

in two dimensions

with K

= (−4,−4)

= (4,4)

= (−6,−6)

, and

= (6,6)

ematically it is a straightforward extension of the 3D

graphics primitive triangular fan to arbitrary dimen-

sions. For graphical examples of the simplicial com-

plexes discussed in this paper see Figure 1 and 2.

In Section 2 we deﬁne the simplicial complex

mathematically. In Section 3 we give a short descrip-

tion of how to include Armadillo in a VS2012 project

and discuss the basics of the Armadillo library and

then we deﬁne in Section 4 the data-structures

Grid,

zJs

, and

T std NK

used to describe the simplicial

complex. In Section 5 we implement the construc-

tion of the complex. We then discuss the efﬁcient im-

plementation of some non-trivial algorithms for the

simplicial complex in Section 6 before making some

conclusions at the end.

Figure 2: A schematic picture of the simplicial complex

std

N,K

in three dimensions. By adding the origin as a vertex

to all the simplices in the simplicial 2-complex subdividing

the boundary of the hypercube we get a fan-like simplicial

3-complex (tetrahedra) locally at the origin.

2 SIMPLICIAL COMPLEX T

std

N,K

To deﬁne the simplicial omplex T

std

N,K

we ﬁrst give

a few deﬁnitions. We denote by Z, N

, and R the

sets of the integers, the nonnegative integers, and the

real numbers respectively. We write vectors in bold-

face, e.g. x ∈ R

and y ∈ Z

, and their components

as x

,... ,x

and y

,... ,y

. All vectors are as-

sumed to be column vectors. An inequality for vectors

is understood to be component-wise, e.g. x < y means

that all the inequalities x

< y

,... ,x

< y

are fulﬁlled.

The convex combination of an (m + 1)−tuple

,... ,x

) of vectors x

,... ,x

∈

is deﬁned by co(x

,... ,x

) :=

{

∑

i=0

: 0 ≤ λ

≤ 1,

∑

i=0

= 1}. The set of

vectors x

,... ,x

∈ R

is called afﬁnely indepen-

dent if

∑

i=1

− x

) = 0 implies λ

= 0 for all

i = 1, . . . , m. This deﬁnition is independent of the

order of the vectors. If x

,... ,x

∈ R

are afﬁnely

independent the set co(x

,... ,x

) is called an

m-simplex.

A triangulation of a set C ⊂ R

is the subdivision

of C into n-simplices, such that the intersection of

any two different simplices in the subdivision is ei-

ther empty or a k-simplex, 0 ≤ k < n, and then its

vertices are the common vertices of the two different

n-simplices. Such a structure is often referred to as a

simplicial n-complex.

For the deﬁnition of T

std

N,K

we use the set S

of all

permutations of the numbers 1,2,..., n, the character-

istic functions χ

(i) equal to one if i ∈ J and equal

to zero if i /∈ J , the null vector 0 ∈ R

and the stan-

dard orthonormal basis e

,... ,e

of R

. Further,

we use the functions R

: R

→ R

, deﬁned for every

J ⊂ {1,2,..., n} by R

(x) :=

∑

i=1

(−1)

(i)

To construct the triangulation T

std

N,K

, we ﬁrst de-

ﬁne the triangulations T

std

and T

std

K,fan

as intermediate

steps.

1. For every z ∈ N

, every J ⊂ {1,2, ...,n}, and ev-

ery σ ∈ S

deﬁne the simplex

zJ σ

:= co(x

zJ σ

,... ,x

zJ σ

) (1)

where

zJ σ

:= R

∑

j=1

σ( j)

(2)

for i = 0,1,2,...,n.

2. Let N

∈ Z

, N

< 0 < N

, and deﬁne the

hypercube N := {x ∈ R

: N

≤ x ≤ N

}. The

simplicial complex T

std

is deﬁned by

std

:= {S

zJ σ

: S

zJ σ

⊂ N}. (3)

3. Let K

∈ Z

, N

≤ K

< 0 < K

≤ N

and consider the intersections of the n-simplices

z,J ,σ

in T

std

and the boundary of the hyper-

cube K := {x ∈ R

: K

≤ x ≤ K

}. We are

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

only interested in those intersections that are

(n − 1)-simplices, i.e. co(v

,... ,v

) with ex-

actly n-vertices. For every such intersection

add the origin as a vertex to it, i.e. consider

co(0,v

,... ,v

). The set of such constructed

n-simplices is denoted T

std

K,fan

. It is a triangulation

of the hypercube K.

4. Finally, we deﬁne our main simplicial complex

std

N,K

by letting it contain all simplices S

zJ σ

std

, that have an empty intersection with the in-

terior K

◦

of K, and all simplices in the simplicial

fan T

std

K,fan

. It is thus a triangulation of N having a

simplicial fan in K.

We have several remarks on this construction. First,

std

N,K

is indeed a simplicial complex, as can easily be

deducted from the proof of Lemma 3.6 in (Giesl and

Hafstein, 2013). Second, if K

= (−1, − 1,..., −1)

and K

= (1,1,...,1) the complexes T

std

N,K

and T

std

are identical. Third, when using the complex T

std

N,K

to compute CPA Lyapunov functions one most com-

monly uses a transformation F : R

→ R

to de-

form and scale down the simplices, i.e. every sim-

plex co(v

,... ,v

) ∈ T

std

N,K

is mapped to a sim-

plex co(F(v

),F(v

),.. . ,F(v

)). The transformation

F must be chosen such that the resulting set of sim-

plices is a simplicial complex.

3 VS2012 AND ARMADILLO

Before we come to our implementation of the sim-

plicial complex T

std

N,K

we explain how to get a project

using the Armadillo linear algebra library running in

VS2012 on a Windowscomputer. This is by no means

the only nor the most elegant way, but it is very simple

and it works.

First download and install Visual Studio Ex-

press 2012 for Windows Desktop. Then go to

http://arma.sourceforge.net and download and extract

Armadillo. Start VS2012 and choose “FILE→New

Project”. In the window that pops up choose

“Visual C++” and “Console Application” and in

the following check “Empty project”. We assume

for simplicity that the name given to the project

is “SIMP” and that the location is “c:\”. The

folder where our program will be running is then

“c:\SIMP\SIMP”. Where armadillo was extracted, in

the “include” folder, there is a ﬁle named “armadillo”

and a folder named “armadillo

bits”. Copy both to

“c:\SIMP\SIMP”. In the “examples” folder there is a

folder named “lib

win32”. Also copy its contents to

“c:\SIMP\SIMP”. Many functions in Armadillo use

the LAPACK and BLAS libraries and therefore we

have to uncomment (remove “//” in front of)

#define

ARMAUSELAPACK

and

#define ARMAUSEBLAS

in “con-

ﬁg.hpp” in the folder “armadillo

bits” if we want to

use the full functionality of Armadillo.

To actually use the functionality from LAPACK

and BLAS we have to link to these libraries dy-

namically. To enable that choose “DEBUG→SIMPL

Properties”. In the window that pops up choose

“Conﬁguration Properties→Linker→Input” and add

“lapack

win32 MT.lib;blas win32 MT.lib;” (without

the quotation marks) to “Additional Dependencies”.

Do this both with “Conﬁguration:” on “Release” and

“Debug”.

VS2012 has the unexpected feature (error?)

that it does not search for .dll ﬁles in the di-

rectory where the program generated is running,

in our case “c:\SIMP\SIMP”. To change this go

to “Conﬁguration Properties→Debugging” and add

“PATH=%PATH%;$(ProjectDir)” (without the quota-

tion marks) to “Environment”. As before do this both

with “Conﬁguration:” on “Release” and “Debug”.

Now everything should be ready to use Armadillo.

Right-click on “Source ﬁles” in the “Solution Ex-

plorer” and choose “Add New Item”. For simplicity

we use the default, which is aﬁle named “Source.cpp”

in “c:\SIMP\SIMP”.

To test if everything is in place we can e.g. try to

compile and run the following program:

#include "armadillo"

#include<list>

// any other headers we might want to include

using namespace arma;

using namespace std;

int main(int argc, char **argv){

mat A=randu<mat>(5,5);

det(A);

}

For our implementation of the simplicial complex be-

low we need to include

list.

We also use

vector

and

algorithm

from the Standard Template Library

(STL), but they are already included in

armadillo.

To compile and run a console application it is ad-

vantageous to choose “DEBUG→Start Without De-

bugging” (or press Ctrl+F5), for otherwise the con-

sole closes immediately when the program has ﬁn-

ished running. This procedure above has been tested

to work with Armadillo 3.8.0.

Now a few comments on Armadillo: Very

good documentation on the library is available at

http://arma.sourceforge.net and in (Sanderson, 2010).

The vector and matrix types we will use in this paper

are

ivec

vec

, and

mat

, which are column vector of

int

, column vector of

double

, and matrix of

double

respectively. Armadillo starts indexing of vectors and

matrices at zero and not at one, just as in C and C++.

ImplementationofSimplicialComplexesforCPAFunctionsinC++11usingtheArmadilloLinearAlgebraLibrary

Note that Armadillo does not support implicit or ex-

plicit conversions between vector and matrix types

only because they might make sense mathematically.

If e.g.

is a function expecting a

vec

as an argument

we cannot call it with an

ivec vi.

We have to use

convto<vec>::from(vi)

to explicitly convert

vec.

The compiler expects the result of a matrix multi-

plication to be a matrix. If we know that it is a scalar

(1 × 1 matrix) the function

as scalar

can be used,

e.g.

double y=asscalar(x.t()*x);

for a vector

In debug modus

as scalar

will report an error if the

argument is not a 1× 1 matrix, in release modus it will

simply give incorrect results. Using the “

” operator

is a short and readable way to assign values to vec-

tors and matrices (

endr

stands for end row). It, how-

ever, does not work like

pushback

in the STL. Thus

x<<1<<2;

makes

= (1,2)

. But if this is followed

x<<3<<4;

then

= (3,4)

and not

= (1,2,3,4)

Lambda functions in C++11, functions that can

be written within other functions and have access

to their data, are a very nice addition to C++, but

there are some pitfalls when using Armadillo. It is

safer to specify the return value of a lambda function,

for e.g.

[](vec v)

{

return 1*v;

} otherwise returns

an object of type

const eOp<vec,eopscalartimes>,

but

[](vec v)->vec

{

return 1*v;

} returns a

vec.

A further nice addition in C++11 are auto types. Thus,

if the compiler can determine the type of an entity at

compile time, it will assign that type to the entity if it

is declared

auto.

For vectors

x,y

of type

vec

ivec

compar-

isons like

x > y

return a vector with 1 in the entities

where the inequality holds true and 0 otherwise. Thus

(1,2,3)

< (2,2, 4)

results in the vector (1,0,1)

. If

we want a simple true or false answer to whether the

inequality is true for all components we can e.g. use

min(x-y) > 0

With the compiler set to “debug” operations in Ar-

madillo are orders of magnitude slower than with the

compiler set to “release”.

4 THE DATA STRUCTURES

We use the data structures

Grid, zJs

, and

T std NK

to implement the simplicial complex T

std

N,K

Grid

used to enumerate the vectors in {z ∈ Z

: z ∈ N}

and other similar grids,

zJs

is a simple container for

z ∈ N

, J ⊂ {1,2,...,n}, and σ ∈ S

and is used to lo-

cate the simplex S

zJ σ

in the data structure

T std NK

which contains all necessary information on the com-

plex T

std

N,K

. For simplicity and to shorten the code we

use very short variable names and do not care about

data encapsulation at all. Further, we omit declaring

functions const and using const references and un-

signed types, even where it would be more natural and

efﬁcient. We make heavy use of the STL in C++ and

assume that the dimension, i.e. n in R

, is already de-

ﬁned by e.g.

int n=4;

The data structure

Grid

, initialized with

ivec Nm

and

ivec pN

contains all the vertices in G(Nm,Np) :=

{z ∈ Z

: Nm ≤ z ≤ Np}. It assigns a unique integer

to each of these vertices and can calculate the corre-

sponding vertex from this number and vise versa. It is

deﬁned as follows:

struct Grid {

ivec mN, pN;

int EndI;

int V2I(ivec);

ivec I2V(int);

vector<int> V2I(vector<ivec> v);

bool InGrid(ivec v);

Grid(ivec _mN,ivec _pN);

˜Grid() {};

};

The numbers assigned to the vectors are 0,1, . . .,N,

where N = EndI − 1. Thus, the constructor can be

coded

Grid::Grid(ivec _Nm,ivec _Np):Nm(_Nm),Np(_Np){

EndI=1;

for(int i=0;i<n;i++){EndI *= Np(i)-Nm(i)+1;}

}

A simple method to assign unique numbers to the ver-

tices is to translate them with the vector

-Nm

and then

enumerate them starting at the origin. The reverse

process is done by using repeated division with re-

minder. The following implementation of the pair

V2I

(vertex to index) and

I2V

(index to vertex) should il-

luminate the strategy:

int Grid::V2I(ivec v){

int i, Index, Mul;

v -= Nm;

for(i=0,Mul=1,Index=0;i<n;i++){

Index += v(i)*Mul;

Mul *= Np(i)-Nm(i)+1;

}

return Index;

}

ivec Grid::I2V(int Index){

ivec v(n);

for(int i=0;i<n;i++){

v(i) = Index%(Np(i)-Nm(i)+1);

Index /= Np(i)-Nm(i)+1;

}

return v += Nm;

}

Further, it is advantageous to be able to pass a

vector of vertices to

V2I

. This is implemented by

vector<int> Grid::V2I(vector<ivec> v),

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

where it can be seen how lambda functions can lead

to efﬁcient and readable code:

vector<int> Grid::V2I(vector<ivec> v){

vector<int> iv;

for_each(v.begin(),v.end(),

[&](ivec &val){iv.push_back(V2I(val));}

);

return iv;

}

The

[&]

allows the lambda function access to all vari-

ables of the enclosing function by reference, in this

case

and

this

. Note that the call to

V2I(val)

is an abbreviation for

this->V2I(val)

and thus the

this

pointer is implicitly used. If we want the lambda

function to use copies of the variables by default we

should replace

[&]

[=]

. We could also have listed

their access mode individually by

[iv&,this]

, be-

cause we need to modify

in the lambda function

but not

this

Only one more member function is needed for

Grid,

bool InGrid(ivec v),

which returns

true

v ∈ G(Nm,Np) and

false

otherwise. Here the Ar-

madillo functions

min

and

max,

which deliver the

minimum and maximum values of a vector respec-

tively, are useful:

bool Grid::InGrid(ivec v){

return min(v-Nm) >= 0 && max(v-Np) <= 0;

}

We now come to the structure

zJs

. It is a simple con-

tainer, on which we deﬁne an ordering relation “<”.

The ordering is used by

T std NK

to sort and then ﬁnd

simplices referred to by z ∈ N

, J ⊂ {1,2,...,n}, and

σ ∈ S

quickly. The variable

int Pos

zJs

is the

positioning used by

T std NK

struct zJs {

int J,Pos;

ivec z,sig;

zJs(ivec _z,int _J,ivec _sig,int _Pos=-1):

z(_z), J(_J), sig(_sig), Pos(_Pos) {};

};

The set J ⊂ {1,2,...,n} is stored as an integer

The idea is to use the representation of

as a binary

number to mark which elements of {1,2,... , n} are

in J and which are not. This is best shown by ex-

amples. The number 0 = (00...0000)

is the empty

set, 1 = (0...0001)

is the set {1}, 2 = (0. .. 0010)

is the set {2}, 3 = (0...0011)

is the set {2,1}, and

e.g. 12 = (0. . . 01010)

is the set {4,2}. In general,

j ∈ J if and only if the j-th bit in the binary rep-

resentation of

is 1. To check whether j ∈ J one

can use bit-shifts and the bitwise and-operation “&”,

i.e.

(J>>(j-1))&1

is one if j ∈ J and zero other-

wise. For

int J

this works for n ≤ 31, for

unsigned

long long J

this works for n ≤ 64. For n > 64 this

strategy has to be reﬁned.

The permutation σ ∈ S

is stored as an

ivec

sigma

in its one-line notation, i.e. it is deﬁned through

sigma[i]

= σ(i). Here the fact that Armadillo starts

indexing of vectors at zero is a little confusing, be-

cause

sigma

is a reordering of the indices. Thus

sigma[0], sigma[1],

...

, sigma[n-1]

is actually a

permutation of the numbers 0,1, . . .,n− 1 rather than

the numbers 1,2,...,n. We discuss the interplay be-

tween

and

sigma

in more detail in the next sec-

tion, when we give the implementation of

x zJs i

that computes the vertices x

zJ σ

according to the for-

mula (2).

The ordering relation on

zJs

is rather arbitrary,

it should just order objects of type

zJs

according to

z,J,

and

sig

adequate to the STL functions

sort

and

equalrange

somehow.

Pos

should not be considered

in the ordering. The following deﬁnition does the job

just ﬁne:

bool operator<(zJs lhs,zJs rhs){

if(lhs.J != rhs.J) return lhs.J < rhs.J;

int i;

for(i=0;i<n && lhs.z(i)==rhs.z(i);i++);

if(i!=n){return lhs.z(i)<rhs.z(i);}

for(i=0;i<n && lhs.sig(i)==rhs.sig(i);i++);

if(i!=n){return lhs.sig(i)<rhs.sig(i);}

return false; // they are equal

}

We come to the main structure

T std NK

that de-

scribes the simplicial complex T

std

N,K

. It is deﬁned as

follows:

struct T_std_NK {

ivec Nm,Np,Km,Kp;

Grid G;

int Nr0;

vector<ivec> Ver;

vector<vector<int>> Sim;

vector<zJs> NrInSim;

vector<int> Fan;

int InSimpNr(vec x); // -1 if not found

bool InSimp(vec x,int ind);

T_std_NK(ivec Nm,ivec Np,ivec Km,ivec Kp);

};

= N

and

= N

deﬁne the hypercube N and

= K

and

= K

deﬁne the hypercube K from

Section 2.

is a grid deﬁned by

and

and is used

to have a coherent enumerationof all vertices possibly

used by

T std NK

Ver

is a vector containing all the

vertices of all the simplices in the complex and

Nr0

is the position of the zero vector/vertex in this vector,

i.e.

Ver[Nr0]

is the zero-vector.

Sim

is a vector con-

taining all the simplices of the complex. A simplex is

basically (n+ 1) vertices. Each simplex is stored as a

vector of (n+1)-integers, the integers refereing to the

positions of the corresponding vertices in

Ver

The remaining members are not used for the con-

struction of the simplicial complex. They are, how-

ImplementationofSimplicialComplexesforCPAFunctionsinC++11usingtheArmadilloLinearAlgebraLibrary

ever, advantageous if one wants to use the simplicial

complex as a basis to deﬁne CPA functions, because

given a vector x in the triangulated hypercube N, they

enable the fast search of a simplex S such that x ∈ S.

NrInSim

contains all simplices of thekind S

zJ σ

in the

complex sorted according to the ordering on

zJs

Fan

contains the rest of the simplices, i.e. the simplices

in the simplicial fan at the origin. We discuss this in

more detail after the next section, in which we discuss

the construction of the simplicial complex

T std NK

5 CONSTRUCTION OF T

std NK

To construct the simplicial complex

T std NK

need a function to compute the vertices x

zJ σ

as in for-

mula (2), i = 0, 1, . . .,n, for the simplices S

zJ σ

. As

mentioned in the last section the interplay between

and

sigma

here play a little confusing role. Be-

cause the set

is supposed to contain the indices of

those coordinates of a vector v, whose coordinates

should be multiplied with minus one, and Armadillo

starts indexing of vectors at zero, we should multi-

ply the coordinate

v[j],

which corresponds to the

coordinate v

j+1

of v, by minus one, if and only if

(J>>((j+1)-1))&1,

i.e.

(J>>j)&1,

is equal to one.

Further,

sigma

is actually a permutation of the num-

bers 0,1,... , n − 1 as discussed above. The formula

(2) for x

zJ σ

, i = 0,1,..., n, can thus be implemented

as follows:

ivec x_zJs_i(ivec z,int J,ivec sigma,int i){

ivec x_s_i=zeros<ivec>(n), v(n);

for(int j=0;j<i;j++){

x_s_i(sigma(j))=1;

}

for(int j=0;j<n;j++){

v(j)=((J>>j)&1 ? -1:1)*(z(j)+x_s_i(j));

}

return v;

}

We now have everything we need to actually con-

struct

T std NK

. The code for the construction can

be partitioned into three parts. In the ﬁrst part some

variables and class members are initialized. This is

done in an initializer list and at the beginning of the

function. In the second part we actually construct

the simplicial complex. This involves a triple loop,

for we have to iterate over all relevant z ∈ N

, all

J ⊂ {1, 2, . ..,n}, and all σ ∈ S

, cf. formulas (1) and

(2). In the third part we tidy up, which includes sort-

ing some vectors to make them eligible for binary

search, removing duplicates, etc. The body of the im-

plementation for the constructor is as follows:

T_std_NK::T_std_NK(ivec _Nm,ivec _Np,

ivec _Km,ivec _Kp) : Nm(_Nm), Np(_Np),

Km(_Km),Kp(_Kp),G(_Nm,_Np) {

// FURTHER INITIALIZATION

int EndSet=1<<n;

ivec ZV=zeros<ivec>(n),pQ1N(n),

IdPerm(n),sigma(n),*pivec,z;

int N=max(max(Np),max(-Nm));

pQ1N.fill(N-1);

for(int i=0;i<n;i++) IdPerm(i)=i;

vector<ivec> sver(n+1);

Grid Q1(ZV,pQ1N);

Grid Ki(Km+1,Kp-1);

// ACTUAL CONSTRUCTION OF THE COMPLEX

for(int J=0;J<EndSet;J++){

for(int zNr=0;zNr<Q1.EndI;zNr++){

z=Q1.I2V(zNr);

sigma = IdPerm;

auto sb=sigma.begin(),se=sigma.end();

do{

// CODE BLOCK 1

// ...

}while(next_permutation(sb,se))

;

}

// TIDY UP

// CODE BLOCK 2

// ...

}

We ﬁrst concentrate on the initialization, the imple-

mentation of

CODE BLOCK 1

and

CODE BLOCK 2

given below. In the initializer list we assign values

to the pairs

Nm, Np

and

Km, Kp

. They correspond to

the vectors N

and K

respectively, that de-

ﬁne the hypercubes N and K as in Section 2. N \ K

◦

is triangulated using the simplices S

zJ σ

and K is tri-

angulated using a simplicial fan. The grid

G(Nm,Np)

includes all vectorsz ∈ Z

that might be vertices in the

triangulation.

EndSet

:= 2

is chosen such that every

subset J of {1, 2, . ..,n} has a unique representation as

a number

= 0,1,..., EndSet −1 as described above.

The grid

is deﬁned with just enough vectors z∈ N

to sufﬁce for the construction of all S

zJ σ

relevant for

std

, cf. (3). The grid

is deﬁned such that the rel-

evant intersections of simplices S

zJ σ

⊂ N with the

boundary of K := {x ∈ R

: K

≤ x ≤ K

} are char-

acterized by having exactly one vertex in

. That we

get any relevant intersection by this characterization is

quite clear. The fact that we get every relevant inter-

section no more than once can be deduced by consid-

ering the intersection of two different such simplices,

which would clearly not be an allowed intersection of

two different simplices in a simplicial complex.

IdPerm

is deﬁned to be the permutation

IdPerm[i]=i

for i = 0,1, ...,n − 1. The func-

tion

nextpermutation

from the STL considers

this to be the ﬁrst permutation. Successive calls to

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

nextpermutation

then iterates through all possible

permutations.

For the actual construction of the simplicial com-

plex we iterate over all z ∈

⊂ N

, all permutations

sigma

of the numbers 0,1,... , n− 1, and all subsets J

of {1, 2, . . .,n}. The z are represented through their

unique numbers in

, the permutations are repre-

sented through their one-line form, and the subsets

through numbers 0,1,..., 2

− 1. The code for the ac-

tual construction is as follows:

// CODE BLOC 1 - IMPLEMENTATION

for(int i=0;i<=n;i++){

sver[i] = x_zJs_i(z,J,sigma,i);

}

int NrInN=0,NrInKi=0;

for_each(sver.begin(),sver.end(),

[&] (ivec &v) {

if(G.InGrid(v)) NrInN++;

if(Ki.InGrid(v)){

pivec=&v; NrInKi++;

}

);

if(NrInN == n+1){

if(NrInKi == 0){

Sim.push_back(G.V2I(sver));

int SLE=Sim.end()-Sim.begin()-1;

NrInSim.push_back(zJs(z,J,sigma,SLE));

} else if(NrInKi == 1){

*pivec=sver[0];

sver[0]=ZV;

Sim.push_back(G.V2I(sver));

int SLE=Sim.end()-Sim.begin()-1;

Fan.push_back(SLE);

}

First we construct the simplex S

zJ σ

co(x

zJ σ

,... ,x

zJ σ

) by writing its vertices to

sver,

i.e.

sver[i]

:=x

zJ σ

for i = 0,1,2...,n.

Then we count how many of the vertices are in

N and K

◦

. Note that x

zJ σ

∈ N if and only if

G.InGrid(sver[i]) == true

and x

zJ σ

∈ K

◦

and only if

Ki.InGrid(sver[i]) == true.

zJ σ

∈ K

◦

we tactically store a pointer to its corre-

sponding

sver[i].

Then we verify if S

zJ σ

∈ T

std

i.e. if S

zJ σ

⊂ N, which holds true if and ony if

NrInN == n+1.

Now there are two relevant cases.

One is if no vertex is in K

◦

, i.e.

NrInKi == 0.

Then

the simplex is added as is to

Sim

, however, using the

unique numbers given to its vertices by

. We then

record its position in

NrInSim

to give fast access to

it later through the data structure

zJs

. If exactly one

vertex, say

sver[i]

= x

zJ σ

, is in K

◦

, i.e.

NrInKi

== 1,

then we modify this simplex and add it to

the simplicial fan. We ﬁrst copy

sver[0]

= x

zJ σ

sver[i]

and then the zero vector

sver[0]

Then we add it to

Sim

and record its position in

Fan.

Now that we have constructed the simplicial com-

plex we tidy up and prepare the simplicial complex

for efﬁcient application. This is implemented as fol-

lows:

// CODE BLOC 2 - IMPLEMENTATION

// record all vertices

list<int> lv;

for_each(Sim.begin(),Sim.end(),

[&](vector<int> &v){

for_each(v.begin(),v.end(),

[&](int iv){

lv.push_back(iv);

}

);

}

);

lv.sort();

lv.unique();

vector<int> vID(lv.size());

vID.assign(lv.begin(),lv.end());

Nr0=equal_range(vID.begin(),vID.end(),

G.V2I(ZV)).first - vID.begin();

// let the simplices in "Sim" refer to the

// vertices by their positions in vID rather

// than their ID-number from "Grid G"

for_each(Sim.begin(),Sim.end(),

[&](vector<int> &v){

for_each(v.begin(),v.end(),

[&](int &iv){

iv=equal_range(vID.begin(),vID.end(),

iv).first - vID.begin();

}

);

}

);

// record the vertices in "Ver" in the same

// order as in vID

for_each(vID.begin(),vID.end(),

[&](int vID){

Ver.push_back(G.I2V(vID));

}

);

// sort "NrInSim" and "Fan" for binary search

sort(NrInSim.begin(),NrInSim.end());

sort(Fan.begin(),Fan.end());

To record all the vertices in the complex we ﬁrst add

all vertex ID-numbers of all simplices to the list

sort it and use

unique

to remove duplicates. Then we

copy the contents of

to the vector

vID

to be able

to apply efﬁcient binary search, i.e. the STL function

equalrange.

Further, we record the position

Nr0

the zero vertex in

vID.

Then we go through all vertices of all the sim-

plices and replace the ID-number of every vertex by

its position in

vID

. Then we actually construct the ver-

tices as

ivec

and write them in

Ver

in the same order

as in

vID

. Thus

Ver[Sim[k][i]]

is the i-th vertex

of the k-th simplex in

Sim

. Finally, we sort

NrInSim

ImplementationofSimplicialComplexesforCPAFunctionsinC++11usingtheArmadilloLinearAlgebraLibrary

and

Fan,

again to be able to efﬁciently apply the STL

function

equalrange.

Thus given

z,J,

and

sigma

for a simplex S

zJ σ

we can efﬁciently locate it in

Sim

and given the position of a simplex in

Sim

we can efﬁ-

ciently check whether it is in the simplicial fan, i.e. in

Fan,

or not. We take advantage of the former prop-

erty in the function

int InSimplexNr(vec x),

dis-

cussed in the next section. The second property is not

important for the applications described in this paper,

but is useful for other applications.

6 ALGORITHMS FOR T

std NK

If the data structure

T std NK

is to be useful for serv-

ing as a basis for CPA functions, we have to be able to

efﬁciently solve the following problem: For an arbi-

trary x ∈ N ﬁnd a simplex S ∈ T

std

N,K

such that x ∈ S.

In this section we implement this efﬁciently given that

the simplicial fan contains a small fraction of the total

number of simplices in the complex. If this is not the

case a different strategy should be used, e.g. storing

an appropriate

zJs

for simplices in the simplicial fan

and project x ∈ K to the boundary of K and search

for this appropriate

zJs.

For demonstrating our ideas

the following is, however, more informative, because

it includes ideas necessary to solve this problem when

std

N,K

has been deformed as explained in Section 2.

Given a simplex S = co(v

,... ,v

) and a

vector x ∈ R

, x ∈ S if and only if x can be written

as a convex combination of the vertices of the sim-

plex. This means that there are nonnegative numbers

,λ

,... ,λ

, such that

x =

∑

i=0

, where

∑

i=0

= 1, (4)

which in turn is equivalent to

x− v

∑

i=1

− v

), where

∑

i=1

≤ 1. (5)

We construct the matrix X by writing the vectors

− v

,... ,v

− v

in its columns conse-

quently. Because the vertices of a simplex are afﬁnely

independent, the equation x − v

= XL always has

a solution L = (λ

,λ

,... ,λ

)

. If the solution ful-

ﬁlls λ

≥ 0 for all i = 1,2,... , n and

∑

i=1

≤ 1, then

x ∈ S, otherwise x /∈ S. Thus we can implement

bool T std NK::InSimp(vec x,int ind),

which

returns

true

vec x

is in the simplex

Sim[ind]

and

false

otherwise, as follows:

bool T_std_NK::InSimp(vec x,int ind){

vector<vec> v(n+1);

for(int i=0;i<=n;i++){

ivec t=Ver[Sim[ind][i]];

v[i]=conv_to<vec>::from(t);

}

mat X(n,n);

for(int i=1;i<=n;i++){

X.col(i-1)=v[i]-v[0];

}

vec L=solve(X,x-v[0]);

return min(L) >= 0 && sum(L) <= 1;

}

The code is self explaining. The connection to CPA

functions f : N → R, deﬁned by giving its values f(v)

at every vertex v of every simplex of T

std

N,K

is as fol-

lows: If x =

∑

i=0

∈ co(v

,... ,v

), then

f(x) = f(

∑

i=0

) =

∑

i=0

f(v

), (6)

as can be easily veriﬁed. Because λ

= 1 −

∑

i=1

the solution L thus gives us a formula for the function

value.

Going through all simplices S ∈ T

std

N,K

to check

whether a given x ∈ S is not very efﬁcient and if x ∈

N \ K

◦

we can do much better. In this case we know

that x ∈ S

zJ σ

for some z ∈ N

, J ⊂ {1,2,...,n}, and

σ ∈ S

, and if we calculate z, J , and σ directly from

x we can ﬁnd the simplex S

zJ σ

using binary search

in the vector

NrInSim.

To compute z and J we ﬁrst

construct a vector y = (y

,... ,y

)

and an inte-

ger J. We do this by going through the entities x

x. If x

≥ 0 we set y

= x

and the i-th bit of the bi-

nary representation of J equal to zero. If x

< 0 we

set y

= −x

and the i-th bit in the binary representa-

tion of J equal to one. After this procedure the inte-

ger J characterizes the set J as discussed in Section 4

and z = (z

,... , z

)

can be computed by z

= ⌊y

⌋

(largest integer ≤ y

) for i = 1,2, ...,n. Now clearly

y ∈ S

0σ

, i.e. y ∈ S

∗

with J

∗

0 the empty set,

and it is easily veriﬁed that

w := y− z =

∑

i=1

∑

j=1

σ( j)

∑

i=1

∑

j=i

σ(i)

. (7)

Because all the λ

are nonnegative we have the rela-

tion

1 ≥ w

σ(1)

≥ w

σ(2)

≥ ... ≥ w

σ(n)

≥ 0 (8)

for w = (w

,... , w

)

. Thus, if we sort the en-

tities of w in decreasing order and record the corre-

sponding indices we have computed σ. The function

sortindex(w,1)

in Armadillo does exactly this (the

optional argument

speciﬁes descending sort).

In the implementation of

int

std NK::InSimpNr(vec x),

that returns the

position in

Sim

of a simplex containing

vec x

possible and −1 otherwise, we ﬁrst check if

in the domain of the simplicial complex. Then we

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

check if

is (only) in the domain of the simplicial

fan. If it is we go through all simplices in the fan to

ﬁnd a simplex containing

. If

is in the domain of

the simplicial complex, but not in the fan, we use the

efﬁcient strategy of computing a simplex containing

it as described above. The code has, obviously, to

be adapted to the fact that Armadillo indexes vectors

from zero. The implementation is as follows:

int T_std_NK::InSimpNr(vec x){

if(!(min(Np-x)>0 && min(x-Nm)>0)){

return -1;

}

if(min(Kp-x)>0 && min(x-Km)>0){

for(int i=0;i<Fan.size();i++){

if(InSimp(x,Fan[i])){

return Fan[i];

}

// WE CAN COMPUTE THE POSITION OF THE SIMPLEX

int J=0;

ivec z(n),sig; // sig=sigma

for(int i=0;i<n;i++){

if(x(i)<0){

x(i)=-x(i);

J |= 1<<i;

}

z(i)=static_cast<int>(x(i));

}

sig=conv_to<ivec>::from(sort_index(x-z,1));

return equal_range(NrInSim.begin(),

NrInSim.end(),zJs(z,J,sig)).first->Pos;

}

We havea fewcomments on this implementation. The

command

J |= 1<<i;

sets the (i+ 1)-th bit of the bi-

nary representation of J equal to 1. Because the enti-

ties of

are all nonnegativewhen we want to compute

their ﬂoor, we can simply cast from

double

int

The Armadillo function

sortindex

returns a vector

of unsigned integers that describes the sorted order

of the given vector’s elements. The optional second

parameter can be set to 1 to let

sortindex

use de-

scending sort, otherwise it uses the default, which is

ascending sort.

7 SUMMARY

We described the implementation of a simplicial com-

plex with a simplicial fan at the origin. Such com-

plexes allow for the parameterizations of continuous,

piecewise afﬁne (CPA) functions, with an arbitrary

rich structure at a singularity. Such CPA functions

have been shown to be irreplaceable in the computa-

tion of true CPA Lyapunov functions for general non-

linear systems. We used C++11 and the Armadillo

linear algebra library for the implementation and we

discussed some of the advantages of doing so in the

paper. Thus, the paper might be of interest to sci-

entists and engineers interested in modern numerical

programming in C++11 under Windows, even if they

are not necessarily interested in our particular prob-

lem of implementing simplicial complexes for CPA

functions.

REFERENCES

Eghbal, N., Pariz, N., and Karimpour, A. (2012). Discon-

tinuous piecewise quadratic Lyapunov functions for

planar piecewise afﬁne systems. J. Math. Anal. Appl.,

399, pp. 586–593.

Giesl, P. (2007). Construction of Global Lyapunov Func-

tions Using Radial Basis Functions. Lecture Notes in

Mathematics, 1904, Springer.

Giesl, P. and Hafstein, S. (2012). Existence of piecewise lin-

ear Lyapunov functions in arbitary dimensions, Dis-

crete Contin. Dyn. Syst., 32, pp. 3539–3565.

Giesl, P. and Hafstein, S. (2013). Revised CPA method to

compute Lyapunov functions for nonlinear systems.

(submitted)

Hafstein, S. (2004). A constructive converse Lyapunov the-

orem on exponential stability. Discrete Contin. Dyn.

Syst., 10, pp. 657–678.

Hafstein, S. (2005). A constructive converse Lyapunov

theorem on asymptotic stability for nonlinear au-

tonomous ordinary differential equations. Dynamical

Systems, 20, pp. 281–299

Johansen, T. (2000). Computation of Lyapunov Functions

for Smooth Nonlinear Systems using Convex Opti-

mization. Automatica, 36, pp. 1617–1626.

Johansson, M. and Rantzer, A. (1997). On the computa-

tion of piecewise quadratic Lyapunov functions. In:

Proceedings of the 36th IEEE Conference on Decision

and Control.

Julian, P., Guivant, J., and Desages, A. (1999). A

parametrization of piecewise linear Lyapunov func-

tion via linear programming Int. Journal of Control,

72, pp. 702–715.

Marinosson, S. (2002a). Lyapunov function construction

for ordinary differential equations with linear pro-

gramming. Dynamical Systems, 17, pp. 137–150.

Marinosson, S. (2002b). Stability analysis of nonlinear sys-

tems with linear programming: A Lyapunov functions

based approach. Ph.D. Thesis: Gerhard-Mercator-

University, Duisburg, Germany.

Peet, M. and Papachristodoulou, A. (2010). A converse

sum-of-squares Lyapunov result: An existence proof

based on the Picard iteration. In: 49th IEEE Confer-

ence on Decision and Control, pp. 5949–5954.

Rezaiee-Pajand, M. and Moghaddasie, B. (2012). Estimat-

ing the Region of Attraction via collocation for au-

tonomous nonlinear systems. Struct. Eng. Mech., 41-

2, pp. 263–284.

Sanderson, C. (2010). Armadillo: An Open Source C++

Linear Algebra Library for Fast Prototyping and

Computationally Intensive Experiments. Technical

Report, NICTA.

ImplementationofSimplicialComplexesforCPAFunctionsinC++11usingtheArmadilloLinearAlgebraLibrary