ON COMPUTATIONAL METHODS FOR LYAPUNOV FUNCTIONS

Lyapunov functions are an essential tool in the stability analysis of dynamical systems, both in theory and applications. They provide sufficient conditions for the stability of equilibria or more general invariant sets, as well as for their basin of attraction. The necessity, i.e. the existence of Lyapunov functions, has been studied in converse theorems, however, they do not provide a general method to compute them. Because of their importance in stability analysis, numerous computational construction methods have been developed within the Engineering, Informatics, and Mathematics community. They cover different types of systems such as ordinary differential equations, switched systems, non-smooth systems, discrete-time systems etc., and employ different methods such as series expansion, linear programming, linear matrix inequalities, collocation methods, algebraic methods, set-theoretic methods, and many others. This review brings these different methods together. First, the different types of systems, where Lyapunov functions are used, are briefly discussed. In the main part, the computational methods are presented, ordered by the type of method used to construct a Lyapunov function.


1.
Introduction. In 1892 Lyapunov published his famous doctoral dissertation [158], where he introduced a sufficient condition for the stability of a nonlinear system, namely, the existence of a positive definite function decreasing along the solution trajectories. In honor of his work, such functions are habitually referred to as Lyapunov functions. Although his work was motivated by problems in astronomy such as the stability of the motion of the planets, the concept of a Lyapunov function has turned out to be extremely fruitful in many other different contexts. In the mathematical theory of dynamical systems the Lyapunov theory is commonly considered the most useful general theory for stability and it is a mainstay in the engineering discipline of control theory. Since dynamical systems are a main tool for modeling in the applied sciences, Lyapunov functions appear in as various The general problem of constructing a Lyapunov function is a very hard problem. There have been numerous attempts and methods in the literature of how to compute Lyapunov functions for various kinds of systems. Some of them use a physical insight into the system to have a good intuition about a candidate for a Lyapunov function, others use more systematic methods, including numerical algorithms. These methods have come from different communities in engineering, mathematics, and computer science. To complicate things, the terminology is often not consistent between the various communities.
In this review, we give a short description of the different types of Lyapunov and Lyapunov-like functions for various systems, that have been discussed in the literature, see Section 2. In the main part, Section 3, we introduce the various methods that have been proposed to compute such functions and try to explain the main ideas behind them. We have put much effort into covering the subject as broadly as possible, but because of the vast number of publications on the subject, we are bound to miss some important contributions. We apologize in advance.
It is our hope that this review may help bringing together the different communities that are working on the numerical computation of Lyapunov functions and may be of some use to researchers interested in numerical methods for dynamical systems.
2. Lyapunov functions and tools in stability theory. In this section we will present the concept of a Lyapunov function in different contexts before we discuss computational methods for its construction in Section 3. After introducing the concept of a (strict) Lyapunov function in the original setting of an equilibrium for an autonomous Ordinary Differential Equation (ODE), we first consider more general settings and systems, before we then describe more general types of Lyapunov functions. The first group of systems starts with general invariant sets for general dynamical systems in Section 2.1, discusses non-autonomous ODEs (Section 2.2), switched systems (Section 2.3), infinite-dimensional systems, e.g. generated by PDEs (Section 2.4), control systems and ISS (input-to-state stability) in Section 2.5 and finally random dynamical systems (Section 2.6). The second group, including generalized notions of Lyapunov functions, begins with non-smooth Lyapunov functions in Section 2.7 and discusses complete Lyapunov functions (Section 2.8) and vector, matrix, and multiple Lyapunov functions (Section 2.9). Furthermore, contraction metrics and Finsler-Lyapunov functions as a related method to study stability are introduced in Section 2.10 and Section 2 ends with general remarks on equations and inequalities as main conditions on Lyapunov functions (Section 2.11). There are certainly more generalizations and related concepts; we have focussed on the above mentioned ones as they have led to numerous numerical construction methods, which will be discussed in Section 3.
Lyapunov 1892 [158] introduced what he called his second method, and what later became known as Lyapunov functions, as a method to show (asymptotic) stability of an equilibrium. Moreover, they can be used to determine its basin of attraction, i.e. to determine which trajectories have a certain long-time behavior. Lyapunov, due to his main applications to astronomy, was more interested in stability, i.e. solutions staying close to the equilibrium. Hence, he used Lyapunov functions, which are non-increasing along solutions, a notion of passivity. Due to computational advantages, we are rather interested in asymptotic stability, where solutions stay close to the equilibrium and additionally converge towards the equilibrium. This corresponds to strict Lyapunov functions, which are a generalization of the energy in a dissipative physical system, where the energy is decreasing along trajectories and thus the solutions of the system converge to a (local) minimum of the energy.
Let us consider the following simple case to explain the general ideas. We consider the autonomous Ordinary Differential Equation (ODE) where f ∈ C 1 (R n , R n ), i.e. f : R n → R n is continuously differentiable. The system (1) with an initial value x(0) = ξ has a unique solution, denoted by x(t; ξ). The Picard-Lindelöf Theorem guarantees only local existence in time; we, however, assume that solutions exist for all t ≥ 0. We assume further that (1) has an equilibrium, i.e. a point x 0 ∈ R n such that f (x 0 ) = 0. To present the implications of a Lyapunov function, let us also introduce some basic stability concepts. Denote by · the Euclidian norm on R n . An equilibrium x 0 is called stable if for all > 0 there exists a δ > 0 such that for all ξ with ξ − x 0 < δ we have x(t; ξ) − x 0 < for all t ≥ 0. It is called asymptotically stable if it is stable and, in addition, there exists a δ > 0 such that for all ξ with ξ − x 0 < δ we have lim t→∞ x(t; ξ) − x 0 = 0. It is called exponentially stable if, additionally, the rate of convergence to zero is exponential. There are many further notions of stability such as global, or uniform stability, which we will not discuss here in detail.
The basin of attraction (or domain/region of attraction) of an asymptotically stable equilibrium x 0 is defined by A(x 0 ) = {ξ ∈ R n : lim t→∞ x(t; ξ) − x 0 = 0}.

PETER GIESL AND SIGURDUR HAFSTEIN
A strict Lyapunov function for the equilibrium x 0 of (1) is a real-valued, continuously differentiable function V : U → R, defined on a neighborhood U of x 0 which satisfies • Minimum: V has a minimum at x 0 , i.e. V (x) ≥ 0 for all x ∈ U and V (x) = 0 ⇔ x = x 0 . • Decrease: V is strictly decreasing along solution trajectories of (1) in U except for the equilibrium. A sufficient condition isV (x) < 0 for all x ∈ U \ {x 0 }, whereV denotes the orbital derivative, see below.
If we relax the last requirement toV (x) ≤ 0 for all x ∈ U, then this defines a non-strict Lyapunov function (or just Lyapunov function), in contrast to a strict Lyapunov function, i.e. a function that is non-increasing along solution trajectories of (1). The decrease property was expressed using the orbital derivativeV , the derivative along solution trajectories of (1), which is defined by d dt V (x(t; ξ)) t=0 and can be calculated by the chain rule aṡ where a · b = a T b denotes the scalar product of the vectors a, b ∈ R n (vectors are assumed to be column vectors). The key observation is that one does not need to know the explicit solution to calculate the orbital derivative; the necessary information is given by f . It is, however, not necessary for a Lyapunov function to be differentiable; this will be further discussed in Section 2.7: the property that V strictly decreases/does not increase along solution trajectories can be expressed without using derivatives. The information provided by a Lyapunov function, as detailed below, is independent of the smoothness.
A (strict) Lyapunov function provides the following information: • Stability: The existence of a (strict) Lyapunov function implies the (asymptotic) stability of the equilibrium. • Basin of attraction: Compact sublevel sets of a strict Lyapunov function, which are completely contained in U, are subsets of the basin of attraction of the equilibrium.
To show that the equilibrium is globally asymptotically stable, i.e. the basin of attraction is the whole space, one needs the additional assumption on the strict Lyapunov function that sublevel sets of all levels attained by V are compact. A sufficient condition is lim x →∞ V (x) = ∞, i.e. V is radially unbounded, see Barbašin and Krasovskiȋ 1954 [22].
When characterizing stability and conditions on Lyapunov functions it is often advantageous to resort to so-called comparison functions of the classes K, K ∞ , and L of continuous functions R + 0 → R + 0 , where R + 0 := [0, ∞). Their definitions are as follows: γ ∈ K, if γ(0) = 0 and γ is strictly monotonically increasing. If γ ∈ K and lim r→∞ γ(r) = ∞, then γ ∈ K ∞ . We write γ ∈ L, if γ is monotonically decreasing and lim r→∞ γ(r) = 0. Additionally, it is convenient to define a continuous function ψ : R + 0 × R + 0 → R + 0 to be of class KL, if r → ψ(r, s) is of class K for every fixed s ∈ R + 0 and s → ψ(r, s) is of class L for every fixed r ∈ R + 0 .
With the help of such functions, stability characterization by Lyapunov functions can be written in a convenient and useful form. For example, the global asymptotic stability of the equilibrium x 0 = 0 forẋ = f (x), f locally Lipschitz, can be characterized as: There exists a ψ ∈ KL such that x(t; ξ) ≤ ψ( ξ , t) for all ξ ∈ R n and t ∈ R + 0 , if and only if there exists a V ∈ C ∞ (R n , R), α, β ∈ K ∞ , and γ ∈ K, such that These comparison functions were first introduced by Massera 1956 [167]. For a compendium for comparison function results together with historical remarks see Kellett 2014 [131]. Lyapunov [220] as well as the review in this volume Kellett 2015 [132]. The first main converse theorem for asymptotic stability was given by Massera 1949 [166] and it has been extended by many authors in several directions. However, most existence theorems offer no method to explicitly construct Lyapunov functions, apart from special cases such as linear autonomous ODE, and even for those cases, the actual computation for large systems requires sophisticated numerical methods. Hence, numerical construction methods have been developed, which will be discussed in detail in Section 3.
2.1. Generalizations: dynamical systems and invariant sets. The definition of a Lyapunov function has been generalized in different directions: Firstly, instead of solutions of ODEs one can consider more general dynamical systems. Secondly, instead of an equilibrium one can consider more general compact invariant sets, e.g. periodic orbits. We illustrate these generalizations in more detail.
A dynamical system (X, T, S t ) is a triple of a complete metric space X, called phase space or state space, the time set T, and the map (flow-operator) S t : X → X, defined for all t ∈ T, which is continuous with respect to both t ∈ T and x ∈ X and satisfies the semigroup properties S 0 = id (identity mapping) and S t+s = S t • S s for all t, s ∈ T.
Depending on the time set, we distinguish between continuous-time systems, where T = R + 0 , and discrete-time systems, where T = N 0 . Discrete-time systems are equivalent to iterations of a function g : X → X. The map g is given by the time-1 map of the dynamical system g = S 1 , while the dynamical system is defined through iterations of the map An example for a dynamical system is the ODE system (1), where S t (ξ) = x(t; ξ) defines a continuous-time dynamical system on the phase space X = R n . A discretetime dynamical system (R n , N 0 , S t ) may be defined by the iteration x k+1 = g(x k ) of a continuous function g : R n → R n , where S k is then understood to be the mapping S k x = g •k (x).
A set I ⊂ X is called positively (negatively) invariant if S t I ⊂ I (S t I ⊃ I) for all t ∈ T and invariant if I is both positively and negatively invariant, i.e. S t I = I for all t ∈ T. Examples of compact invariant sets are equilibria or periodic orbits, but they can be more complicated such as, e.g., the well-known Lorenz attractor. A strict Lyapunov function for a compact invariant set I of a general dynamical system is a real-valued function V : U → R, defined on a neighborhood U of I, which satisfies • Minimum: V has a minimum at I, i.e. V (x) ≥ 0 for all x ∈ U and V (x) = 0 ⇔ x ∈ I. • Decrease: V is strictly decreasing along trajectories of the dynamical system in U except the ones in I. A sufficient condition for this isV (x) < 0 for all x ∈ U \ I. A strict Lyapunov function provides the following information: • Stability: The existence of a strict Lyapunov function implies the asymptotic stability of I. • Basin of attraction: Compact sublevel sets of a strict Lyapunov function, which are completely contained in U, are subsets of the basin of attraction of I. Note that, as in the case of an equilibrium, the decrease property can be formulated differently, if V is not differentiable. For a discrete-time dynamical system x k+1 = g(x k ) the decrease property can be expressed by the discrete orbital "derivative" ∆V (x), which is not actually a derivative, defined by

2.2.
Non-autonomous systems. Non-autonomous ODEs of the forṁ where f depends explicitly on the time t, do not define a dynamical system. Rather, they can be described by a skew-product flow, see e.g. Rasmussen 2007 [197] or Kloeden and Rasmussen 2011 [137].
In the non-autonomous setting, equilibria or, more generally, invariant sets are replaced by non-autonomous invariant sets. There are different notions of stability and attractivity of non-autonomous sets depending on the times considered: future, past, or a finite time interval. The stability for future times is called forward stability and is similar to the classical notion. Stability in the past is characterized by so-called pullback stability. When considering a finite time interval, stability can be studied by using Lagrangian coherent structures, see Haller 2000 [113], stable and unstable manifolds, see Berger, Doan, and Siegmund 2009 [27], or finite-time Lyapunov exponents, see Berger 2011 [26].
Examples for such systems are non-autonomous differential and difference equations over infinite-and finite-time as well as random dynamical systems, see also Section 2.6. For converse theorems on Lyapunov functions for non-autonomous systems see Grüne, Kloeden, Siegmund, and Wirth 2007 [106], for finite-time Lyapunov functions see Giesl and Rasmussen 2012 [95].
To clarify the concept of a Lyapunov function for a non-autonomous system, we give a typical Lyapunov theorem for the non-autonomous ODE (2). Assume f ∈ C 1 (R n × R, R n ) and denote by x(t, t 0 ; ξ) the unique solution to (2) with initial REVIEW ON COMPUTATIONAL METHODS FOR LYAPUNOV FUNCTIONS 7 condition x(t 0 ) = ξ. We assume that x(t, t 0 ; 0) = 0 is a solution and we study its stability by a Lyapunov function V (x, t).
Then we have: There exists a ψ ∈ KL and a constant B 1 > 0 such that if and only if there exist V ∈ C 1 (R n × R, R), α, β, γ ∈ K, and a constant B 2 > 0, such that for all s ∈ R and all ξ < B 2 .
In words, the zero solution to (2) where J is some index set. To this end a family of switching signals σ taking values in J is often introduced. Let us consider a simple example of a switched system to explain the concept. Assume σ : R + 0 → J is piecewise constant and let t 0 = 0 < t 1 < t 2 < . . . < ∞ be the times where σ changes its value. The solution x σ (t; ξ) to the switched systemẋ = f σ (x), with x σ (0) = ξ as an initial-value, is then obtained by gluing together solution trajectory pieces of initial-value problems. On the interval [t 0 , t 1 ] we define x σ (t; ξ) to be the solution to the initial-value probleṁ x = f p (x), where p := σ(t 0 ) and x(t 0 ) = ξ, on the interval [t 1 , t 2 ] to be the solution to the initial-value probleṁ x = f p (x), where p := σ(t 1 ) and x(t 0 ) = x σ (t 1 ; ξ), and so on.
One can consider different kinds of switching, e.g. arbitrary switching, where the only limitation on σ : R + 0 → J is that (t i ) i∈N has no accumulation point, statedependent switching, hysteresis switching, etc. For a very readable discussion on the different kinds of switched systems we refer to Liberzon 2003 [151, Chapter 1]. For a concrete example of a switched system see Section 3.3 below.
Closely related to switched systems are differential inclusions. Given a setvalued mapping F from R n to the power-set of R n , often written F : R n ⇒ R n , a solution to the differential inclusionẋ is any absolutely continuous function x(t; ξ), fulfilling x(0; ξ) = ξ andẋ(t; ξ) ∈ F (x(t; ξ)) for almost all t ≥ 0. Often F is given as the convex hull of some functions f p : R n → R n , p ∈ J, i.e. F (x) = co{f p (x) : p ∈ J}. In this case any solution x(t; ξ) to the inclusionẋ ∈ co{f p (x) : p ∈ J} can be approximated by a solution x σ (t; ξ) to the switched systemẋ = f σ (x), with an appropriate switching σ. This result is often referred to as the Filippov-Wažewski Theorem, see also Ingalls, Sontag, and Wang 2002 [121].
Switching between systems with unstable equilibria can result in a system with a stable equilibrium and switching between systems with stable equilibria can result in a system with an unstable equilibrium. An example of the latter is, e.g., given in Agrachev and Liberzon 2001 [2], where even though every convex combination x with an asymptotically stable equilibrium at the origin, the origin is not a stable equilibrium of the inclusionẋ ∈ co{A 1 x, A 2 x}. Switched discrete-time systems and difference inclusions are defined analogously. For a detailed study of switched systems and inclusions, including their stability, we refer to the books Aubin and Cellina 1984 [16], Filippov 1988 [76], Liberzon 2003 [151], and Sun and Ge 2011 [218].

2.4.
Infinite-dimensional dynamical systems. Infinite-dimensional dynamical systems are dynamical systems as described in Section 2.1, where the metric space X is an infinite-dimensional Banach space, i.e. a complete normed vector space, often a space of functions. They can arise from evolution Partial Differential Equations (PDEs), such as the reaction-diffusion equation where u(t, x) is a function depending on the time variable t and the space variable x ∈ R n . Some references, where Lyapunov functions are used to study the stability of PDEs are Xu and Sallet 2002 [226], Coron, d'Andrea-Novel, and Bastin 2007 [61], and Prieur and Mazenc 2012 [194].
Another source for infinite-dimensional dynamical systems are dynamical systems with time delay. As a simple example we can take the time-delayed differential equation with delay τ > 0ẋ = f (x(t − τ ), t), where the initial condition is given by specifying x on the interval [−τ, 0]. The stability of such systems can be studied using Lyapunov-like functions, e.g. based on Razumikhin's theorem or the Lyapunov-Krasovskiȋ Theorem. For a reference on such methods cf. e.g. Krasovskiȋ 1959 [139], Malisoff and Mazenc 2009 [159], or Kharitonov 2013 [134].
2.5. Control-and ISS Lyapunov functions. In a simple setting, a control system is a system of the kindẋ where typically u ∈ U ⊂ R m is interpreted as a control. If u = u(t) is a function of time one speaks of an open-loop control and if u = k(x) is a function of the state one speaks of closed-loop or feedback control. If a feedback k(x) has been fixed, one is left with an autonomous systemẋ = f (x, k(x)). If, for this system, the equilibrium at the origin has some desired stability property, we say that the control system has been stabilized by feedback.
The system is called locally, asymptotically null-controllable, if for every ξ in a neighborhood of the origin there is an open-loop control u such that the solution to the system with initial value ξ approaches the origin asymptotically, i.e. if it is possible to steer the system state asymptotically to the origin. A control Lyapunov function for such a system, introduced by Sontag 1983 [208], is a positive definite function V such that inf for a K function γ, see Section 2 for the definition of K functions. Asymptotic null-controllability cannot be characterized by smooth control Lyapunov functions and one must resort to more general definitions of differentiability like the Dini-or the proximal subdifferential, cf. e.g. Sontag and Sussman 1995 [212] or Clarke 2011 [58], see also Section 2.7. The question whether asymptotic null-controllability is equivalent to stabilizability by feedback has given rise to a lot of research. For a more detailed discussion on control systems we refer to Sontag 1998 [211] and the references therein. For asymptotically null-controllable systems the equilibrium at the origin is sometimes referred to as weakly asymptotically stable, in contrast to strongly asymptotically stable equilibria, where every choice of u leads to states being attracted asymptotically to the equilibrium. Strong asymptotic stability can be characterized by smooth Lyapunov functions as shown in Clarke, Ledyaev, and Stern 1998 [59], i.e. there exists a smooth, positive definite function V such that A different concept is the so-called input-to-state stability or ISS stability, introduced in Sontag 1989 [209], and later characterized by ISS Lyapunov functions in Sontag and Wang 1995 [213] and Sontag 1996 [210]. The essential idea is the characterization of a certain kind of stability of the equilibrium at the origin with an ISS Lyapunov function V , fulfilling where γ and α are K functions. The origin is thus an asymptotically stable equilibrium of the systemẋ = f (x, 0) and a practically stable equilibrium (see end of Section 2.11 for a definition) of the systemẋ = f (x, u) for u ≤ u max , where u max > 0 is a (not too large) constant. Moreover, the smaller u max is, usually interpreted as a bound on the perturbation u, the closer solutions of the system will be to zero in the long run.
2.6. Random dynamical systems. The mathematical discipline of random dynamical systems is very general and includes many kinds of dynamical systems that are subject to some kind of randomness. There are many sensible ways of defining stability for random dynamical systems, e.g., convergence in probability, recurrence, almost sure exponential stability, p-th moment stability, etc. For a thorough treatment of random dynamical systems see Arnold 2002 [14].
Arnold and Schmalfuss 2001 [15] developed the theory of Lyapunov functions for random dynamical systems, including a converse theorem. In a series of papers Liu 2006 even developed a theory of complete Lyapunov functions for random dynamical systems. We will discuss complete Lyapunov functions in Section 2.8 below.
There have been numerous publications on Lyapunov stability for more specific kinds of systems subject to randomness, like stochastic differential equations, stochastic difference inclusions, Markov chains, etc. To name a few, a Lyapunov-like theory for the stability of stochastic differential equations has been developed in Arnold 1974 [13], Mao 2008 [161], and Khasminskii 2012 [135]. The existence of Lyapunov functions for discrete-time, stochastic difference inclusions was proved by Subbaraman and Teel 2013 [215].
To give a concrete idea of how a Lyapunov function for a random dynamical system might look like we consider the stochastic differential equation where f : R n × R → R n and G : R n × R → R n×n are functions, and B : R → R n is an n-dimensional Brownian motion. A Lyapunov function for this system, asserting (uniform) stochastic asymptotic stability of the zero solution (cf. e.g. [161, Definition For other similar theorems cf. e.g. Mao 2008 [161, Chapter 4].
After presenting various different systems, where Lyapunov functions play an important role, we discuss different types of Lyapunov functions in the following sections.

2.7.
Non-smooth Lyapunov functions. For continuous-time dynamical systems, the decrease property can be expressed as a condition on the orbital deriva-tiveV (ξ) = ∇V (ξ) · f (ξ) if both V and f are differentiable. Several computational methods, however, construct Lyapunov functions, which are not differentiable but merely locally Lipschitz. In this case, the orbital derivative can, e.g. be defined by the Dini-derivative, see Dini 1878 [67], aṡ For differential inclusions, appropriate concepts exist in the context of non-smooth calculus, e.g. the Clarke subdifferential, the Dini subdifferential, and the proximal subdifferential, cf. e.g. Clarke 2011 [58]. The main requirement is only that the Lyapunov function decreases in some sense along solution trajectories.

Complete Lyapunov functions.
In contrast to classical Lyapunov functions [158], which are defined on the basin of attraction of just one attractor, a complete Lyapunov function characterizes the flow on the whole phase space and distinguishes between the chain-recurrent set and the gradient-like flow. The decomposition of the flow of a dynamical system into a chain-recurrent part and a part, where the flow is gradient-like, is characterized by a so-called complete Lyapunov function (or sometimes global Lyapunov function) for the system, see Conley 1978 [60]. Originally, the phase space was assumed to be compact, or the (compact) global attractor was considered as the phase space. A complete Lyapunov function is related to the Morse decomposition, producing a partially ordered family of chain-recurrent isolated invariant sets. Later, complete Lyapunov functions were also considered on separable metric spaces for discrete-time Hurley 1998 [120] and continuous-time systems Patrão 2011 [187]. For a thorough study see Akin 2010 [6]. A dynamical system with such a complete Lyapunov function is also called a gradient semigroup, see Aragão-Costa, Caraballo, Carvalho, and Langa 2011 [11]. For generalizations to non-autonomous systems see Rasmussen 2007 [198] and Aragão-Costa, Caraballo, Carvalho, and Langa 2013 [12].
To present the essential ideas, we will consider the case of a general, autonomous ODE on a compact set C ⊂ R ṅ A complete Lyapunov function is a real-valued function V : if and only if ξ ∈ I i , where I i denotes the isolated invariant sets. The construction of a complete Lyapunov function on a compact set C can be achieved as follows: Consider an attractor A i and its dual repeller is then a Lyapunov function for A i , which is 1 on A * i , 0 on A i , and strictly decreasing along solution trajectories on C \ (A i ∪ A * i ). From any series i α i = 1 of numbers α i > 0 we can now construct a complete Lyapunov function V by summing over all, at most countably many, attractors A i , 2.9. Vector, matrix, and multiple Lyapunov functions. Vector Lyapunov functions are a generalization of scalar-valued Lyapunov functions, see Bellman 1962 [24]. Essentially, the components of the vector-valued Lyapunov function act as scalar Lyapunov functions on different parts of the phase space. Vector-valued Lyapunov functions have been applied to and connected with large and interconnected systems, e.g., inŠiljak 1979 [207] as well as Lakshmikantham, Matrosov and Sivasundaram 1991 [141].
Matrix Lyapunov functions formulate the conditions (minimum and decrease) on scalar-valued functions derived from the matrix-valued functions. An application of matrix-valued Lyapunov functions to interconnected systems is given in Martynyuk 1990 [165]. A hybrid dynamic system is a system that combines both continuous-time (flow) and discrete-time (jumps) behaviour. They are closely related to switched systems, but more general, since the solution trajectories are not necessarily continuous. In the context of hybrid dynamic systems, multiple Lyapunov functions were introduced by Branicky 1998 [44]. A typical example would be to specify a Lyapunov function for a set of continuous-time systems and a decrease condition when switching between systems. For reviews of the stability analysis of hybrid systems, that include numerous further references, we refer the reader to Davrazos and Koussoulas 2001 [63] and Shorten, Wirth, Mason, Wulff, and King 2007 [206]. Matrosov 1963 [168] introduced an auxiliary function, the so-called Matrosov function, that, combined with a non-strict Lyapunov function, can be used to prove asymptotic stability, whereas the non-strict Lyapunov function itself merely shows stability. These results have been extended to multiple Matrosov functions, and such Matrosov functions can be used to compute a strict Lyapunov function for the system. For a review on Matrosov functions we refer to the book Malisoff and Mazenc 2009 [159] and the introduction of Subbaraman and Teel 2013 [216], where the authors additionally prove a Matrosov theorem for a stochastic system.

Contraction metric.
A classical Lyapunov function measures the (decreasing) generalized distance between a point and an equilibrium or local attractor. A different way to study stability and attraction, which does not require knowledge about the attractor, is to compare adjacent trajectories with respect to a certain metric. This type of stability is called incremental stability; the metric with respect to which the distance between adjacent trajectories decreases is called contraction metric. In contrast to Lyapunov functions, no information about the attractor is required and the contraction metric is robust under perturbations of the system, even on the attractor.
We include contraction metrics in this review, as some of the computational methods, which will be discussed in Section 3, have also been used to construct contraction metrics. We present the ideas first for the simple case of the phase space X = R n and the autonomous ODĖ Denote by S n the set of symmetric n×n matrices with real entries. For P, Q ∈ S n we write P Q if P − Q is strictly positive definite and P ≺ Q if P − Q is strictly negative definite. Especially, P 0 means that P is strictly positive definite.
A contraction metric M satisfies: where Df (x) ∈ R n×n is the Jacobian matrix of f at x and denotes the orbital derivative of M along solution trajectories. • Positively invariant, compact set: the contraction property (6) is required to hold for all x in a positively invariant, compact set K. A contraction metric provides the following information: • Existence and uniqueness: There exists a unique equilibrium in K.
• Stability: The equilibrium is exponentially stable.
• Basin of attraction: K is a subset of the basin of attraction of the equilibrium. As in the case of Lyapunov functions, conditions replacing (6) have been derived to study different notions of stability. For global stability one needs to assume that the Riemannian metric is uniformly positive definite, i.e. M (x) I for all x with an > 0, and by also assuming a stronger contraction property than (6), see Aghannan and Rouchon 2003 [1]. This concept has been generalized in several directions. Firstly, the Euclidean space R n can be replaced by a general Riemannian manifold M with a Riemannian metric, a family of scalar products on each tangent space T x M. Secondly, the concept of a Riemannian contraction metric can be generalized to a Finsler-Lyapunov function, which is a scalar-valued function on the tangential bundle, see Forni and Sepulchre 2014 [78]. While a Riemannian metric is a scalar product on each tangent space T x M, a Finsler metric is given by a Finsler function on the tangent bundle T M, which defines a Minkowskii norm on each tangent space. A Finsler-Lyapunov function V : T M → R + 0 , defined for every (x, δx) ∈ T M, is such a Finsler function with the following contraction property Depending on the function α, one can deduce certain types of stability. For example, if α is a K function, then the system is incrementally asymptotically stable, i.e. the distance of adjacent solutions converges to 0. Thirdly, instead of the above contraction assumptions in all tangential directions or negative definiteness of a matrix, contraction conditions in certain tangential directions have been studied. These conditions are related to other attractors such as periodic orbits where no contraction along the periodic orbit takes place, so we only have contraction in certain directions; further applications are systems with symmetry. These conditions take, for example, the following form: • Horizontal Finsler-Lyapunov function: the tangent space is decomposed into a horizontal and vertical subspace at each point x: where both V x and H x are spanned by finitely many vectors in the tangent space. For δx ∈ T x M, the decomposition is denoted by only depends on the horizontal component, the condition that V is a Finsler function is relaxed to being a horizontal Finsler function on H x only, and condition (7) still applies. This implies contraction only with respect to the horizontal component, for details see Forni and Sepulchre 2014 [78]. For another example of contraction conditions in certain directions, but otherwise similar to (6), we consider the phase space R n . If the contraction property is only required for directions v perpendicular to the flow, then the long-time behavior is governed by a periodic orbit.
• Riemannian metric: M is a smooth matrix-valued function R n → S n such that M (x) 0 for all x ∈ R n . • Contraction property perpendicular to the flow: • Positively invariant, compact set: the contraction property (8) is required to hold for all x in a positively invariant, compact set K which contains no equilibrium. Such a contraction metric as above provides the following information: • Existence and uniqueness: There exists a unique periodic orbit in K.
• Stability: The periodic orbit is exponentially stable.
• Basin of attraction: K is a subset of the basin of attraction of the periodic orbit. In the literature, contraction metrics appear under different names. Already in Lewis 1949 [147], the distance between solutions with respect to a Riemannian or Finsler metric was studied. The case of a periodic orbit with a contraction condition similar to (8) and the Euclidean metric was developed in Borg 1960 [40]. It was extended by Hartman and Olech 1962 [116], see also Hartman 1964 [115] and Stenstroem 1962 [214], by employing a general Riemannian metric. Further contributions have been made in Krasovskiȋ 1963 [139], Kravchuk, Leonov, and Ponomarenko 1992 [140], see also the books Leonov, Burkin, and Shepelyavyi 1996 [146] and Boichenko, Leonov, and Reitmann 2005 [39]. For global attraction one needs to assume that the Riemannian metric is uniformly positive definite, see Aghannan and Rouchon 2003 [1]. This can be expanded to study the geodesic distance between solutions, and thus to construct a Lyapunov function.
Converse theorems have been proved, showing that a contraction metric exists for an exponentially stable periodic orbit Giesl 2004 [82] and for an exponentially stable equilibrium Giesl 2015 [89]. Generalizations to non-smooth systems Giesl 2005 [83] and finite-time systems Giesl and Rasmussen 2012 [95] have been made.
Incremental stability was studied in Lohmiller and Slotine 1998 [156] and Angeli 2002 [10]. Finsler-Lyapunov functions were introduced by Forni and Sepulchre 2014 [78]. Rüffer, van de Wouw, and Mueller 2013 [202] compare incremental stability with a related notion, namely convergent systems, and establish that they are equivalent on compact sets; moreover, non-autonomous systems are considered.
Other methods to study the stability of periodic orbits include angular one-forms, inspired by Birkhoff 1966 [30], see the expository article Byrnes 2010 [49].
2.11. Equations and Inequalities. The conditions on a Lyapunov function are inequalities outside the invariant set and equations on the invariant set. For contraction metrics, inequalities have to hold both on and outside of the invariant set. There is a fundamental difference between these, in particular when finding solutions algorithmically.
Inequalities are robust under perturbations. In particular, a numerical approximation of a function will satisfy the same inequality when the error is small enough. Also perturbations on the system can be allowed, i.e. a Lyapunov function for one system is also a Lyapunov function for a perturbed system. Hence, numerical methods to construct Lyapunov functions are in principle able to construct a true Lyapunov function (i.e. satisfying an inequality) outside the invariant set, not only in the limit, but even in a finite number of steps, as a small error still produces a Lyapunov function.
This argument does not apply, however, on and near an invariant set I. Therefore many construction methods consider a relaxed version of a Lyapunov function by only assuming the decrease condition outside of I := {x ∈ X : dist(x, I) < }, i.e.
• Decrease: V is strictly decreasing along trajectories of the dynamical system in U except in I ; here, U is a neighborhood of I . A sufficient condition for this isV (x) < 0 for all x ∈ U \ I . Such a Lyapunov function allows only for the conclusion that trajectories will reach I after a finite time and come back infinitely often as time goes to infinity.
This gives rise to an alternative stability concept called practical stability, ultimate boundedness, or invariance kernels with target regions in the different communities. For systems subject to randomness a weak notion of stochastic stability called recurrence is an analogous concept, cf. e.g. Meyn and Tweedie 2009 [170] or Subbaraman and Teel 2013 [216].
3. Numerical methods to compute Lyapunov functions. In this part, the main one of the review, we will present algorithms to compute Lyapunov functions, ordered by the method employed.
A general method to compute Lyapunov functions exists for linear, autonomous ODEs of the formẋ = Ax, where the origin is an exponentially stable equilibrium. Here, the quadratic form V (x) = x T P x is a Lyapunov function, where P is the solution of a matrix equation, the so-called Lyapunov equation (see Section 3.3); even in this case, numerical methods are needed, in particular for high-dimensional systems. For nonlinear systems, the linearization at the equilibrium is a linear system, and the quadratic Lyapunov function of the linearized system is a Lyapunov function for the nonlinear system, however, in general only in a small neighborhood. Hence, it is usually not well-suited to determine large subsets of the basin of attraction.
For nonlinear systems, there exists no general method to compute a Lyapunov function for a given system. In the 1950's, Zubov introduced a particular Lyapunov function as solution of a first-order PDE, the so-called Zubov's equation, see Section 3.1. Solutions of this PDE have been computed using series expansion and other methods.
The numerical computation of Lyapunov functions in the engineering community was originally largely motivated by the desire to prove the Lyapunov stability of systems modeled as difference or differential inclusions of linear systems. Some of these approaches also work for more general systems, but then usually in rather special cases. In these methods, Linear Programming (LP) or Linear Matrix Inequalities (LMI) are generally used to compute a norm on the state space, which is a (strict) Lyapunov function, Section 3.2. Such methods were first introduced around 1980, but have been considerably improved over time and are still a source of active research. There are so numerous publications concerned with these methods that our review is bound to miss some of them.
Around 2000, several new methods were proposed that were directly targeted at nonlinear difference or differential equations. In these methods, a Lyapunov function is usually computed on a compact domain, but in contrast to linearization methods, the size of the domain does not have to be small. In these methods collocation, LP, LMI, algebraic methods, graph theoretic methods, etc., are used to compute various kinds of Lyapunov functions.
Many of the proposed methods to compute Lyapunov functions use various forms of relaxation. In mathematical optimization and related fields the term relaxation is used for the concept of substituting a problem with a different problem that is easier to solve. A solution of the relaxed problem then gives information on the solutions to the more difficult problem. Somewhat counterintuitively, the relaxation problem may be lossy, in the sense that the relaxed problem is less general than the original problem, i.e. the conditions in the relaxed problem are stricter. An example is the SOS method, where the condition of a multivariable polynomial being non-negative is replaced by the condition that the polynomial is representable as the sum of squared polynomials. Clearly a polynomial that is representable as the sum of squared polynomials is non-negative, but there are non-negative polynomials that are not representable as the sum of squared polynomials. The drawback of this relaxation is, however, often outweighed by the fact that verifying the non-negativity in the general case is NP-hard, whereas the decomposability of a polynomial into a sum of squared polynomials can be formulated as an LMI problem. For a more detailed discussion on this we refer to Chapter 4 in Parrilo's PhD thesis 2000 [186].
In the next four sections we will list numerous methods developed to compute Lyapunov functions for various kinds of systems. We will also include some methods to compute contraction metrics, or the basin of attraction with other methods than Lyapunov functions. In some of these methods, the computed function does not fulfill all the properties of a Lyapunov function in the classical sense, but if the essential idea is to compute a real-valued function that is decreasing along solution trajectories of a dynamical system, we will still use the name Lyapunov function. We will discuss the relevant weakening of the classical properties and their implications as we deem appropriate when discussing the methods. For some of these methods one has a good understanding of when the method works. For some of the methods it has been proved that they always work given a certain set of parameters and a certain property of the system in question, e.g. smoothness, special algebraic structure, or a certain kind of stability, thus establishing a "constructive" converse result.
The presentation is ordered by the method used to compute the Lyapunov function. We start with methods related to Zubov's equation and collocation methods in Section 3.1, discuss methods using LP in Section 3.2, present LMI methods in Section 3.3, and end with methods in Section 3.4 that do not fit into the categories of the former sections.
We do not go into much detail on the numerical complexity of the methods, mainly because the complexity, both in time and space, is most often not discussed in much detail in the relevant publications. Note also, that even for linear programming (LP) the complexity is not a closed problem. On a general note, solving linear equations is fast (of complexity O(n 3 ) on one CPU core) in comparison to solving LP problems, and solving LP problems is faster than solving LMI problems. The complexity of LP and LMI is, however, polynomial in the size of the problem. Even for collocation, where linear equations are solved, the so-called curse of dimensionality applies. That is, if the number of parameters needed to determine a Lyapunov function grows exponentially with the dimension of the system, and this is to be expected for nonlinear systems, the number of variables and equations grows exponentially and so does the time needed to compute a Lyapunov function. For the computational complexity of the existence of a Lyapunov function for polynomial systems and other computational problems in dynamical systems see Ahmadi, Majumdar, and Tedrake 2013 [5]. For some examples of how the number of parameters needed to parameterize a Lyapunov function of a given template, e.g. piecewise linear, piecewise quadratic, or a polynomial of a fixed degree, cf. e.g. Giesl and Hafstein [90] and Ahmadi and Jungers 2014 [3] where ϕ(x) is an auxiliary function, satisfying ϕ(x) ≥ 0 and ϕ(x) = 0 if and only if x is the equilibrium. Moreover, the function ϕ has to take into account the approximate order of attraction of solutions to the equilibrium. Then V (x) ∈ [0, 1) is defined on the basin of attraction and V (x) approaches 1 if x approaches the boundary of the basin of attraction.
Different forms of Zubov's equation such aṡ can be obtained by reparameterizing the time, anḋ by the change of variables v(x) = − ln(1 − V (x)) from (9), where v(x) approaches ∞ if x approaches the boundary of the basin of attraction. Note that the inverse Zubov's work became known in the west in a translation of his monograph in 1964 [229], see also a description of the method in Hahn 1967 [112]. For an overview of Zubov's work in the theory of stability and control see Aleksandrov, Martynyuk, and Zhabko 2010 [7]. A generalization of Zubov's method to the case of periodic orbits in autonomous systems and systems of the formẋ = f (x, t), where f is almost periodic with respect to t, are given in Aulbach 1983 [17].
Although in general, finding a solution V for the Zubov equation, a PDE, is not an easier problem than solving the original ODE, Zubov's equation has served as a starting point for several computational methods. Zubov 1964 [229] noted that if f is analytic and ϕ can be chosen as a quadratic form, then V can be expanded into a series where V m (x) is a homogeneous form of m-th degree, and the V i can be determined recursively. When truncating the series after a finite number of terms, one obtains an approximation to V ; this approximation is a Lyapunov function that can determine a subset of the basin of attraction. Moreover, if the basin of attraction is bounded, the approximation also provides an upper bound on the boundary of the basin of attraction. Margolis and Vogts 1963 [162] used a computer to compute both the coefficients of the series recursively up to homogeneous terms of order 14 and the approximate boundary of the basin of attraction. This was the first method to compute Lyapunov functions with a computer, see also Fallside, Patel, Etherton, Margolis, and Vogt 1965 [73]. This method was then adapted to nonlinear autonomous difference equations in O' Shea 1964 [183]. For the Lyapunov equation (which is the equation for V 2 ), the method of Bellmann 1960 [25] using the Kronecker matrix product was employed, which can be extended to higher order homogeneous terms, see Ferguson 1970 [75]. Kinnen and Chen 1968 [136] unified certain approaches to construct Lyapunov functions. Instead of solving a PDE for the orbital derivative, i.e. essentially integrating along solution trajectories, the idea is to consider an exact differential equation, so that the integration is independent of the path. In particular, the integration can be performed along the coordinate axes.
In more detail, a function g is constructed from f such that g(x) · f (x) = 0, for example, . . − f n for i = 1, 2, . . . , n, where x ∈ R n . Now an auxiliary function h has to be found, satisfying the following conditions: (A) the n × n matrix given by the elements m ij = ∂(gi+hi) ∂xj is symmetric, i.e.
These conditions on the existence of a function h are necessary and sufficient for the existence of a Lyapunov function V .
Kinnen and Chen relate their method to several existing methods by choosing h to satisfy the conditions (A) to (C) in different orders, or parameterizing a family. Chen and Kinnen 1970 [54] express h as a sum of two components, one determined by the differential equation, and the other one being an undetermined scalar-valued function ψ. An algorithm is presented to determine ψ by polynomial series expansion.
Gurel and Lapidus 1969 [107] give a detailed survey on methods to construct Lyapunov functions until 1969, including a time-line and extensive references, ordered by date. They distinguish between different types of methods, of which we will discuss three here: the Chetaev-type methods try to build a Lyapunov function by constructing first integrals (so non-strict Lyapunov functions). For a recent book, summarizing methods to construct strict Lyapunov functions from non-strict ones, called strictification methods, for control systems see Malisoff and Mazenc 2009 [159].
The Krasovskiȋ-type methods construct a Lyapunov function as a quadratic form using the right-hand side f , e.g. V (x) = f (x) T Af (x) with suitable positive definite matrix A; the requirement for the decrease property is then that note the similarity to the condition on contraction metrics. Generalizations include an ansatz of the form V (x) = x T A(x)x.
Zubov-type methods construct v as the solution of the PDE (10). The variable gradient method makes an ansatz for ∇v(x) = A(x)x and then determines A(x) to makev(x) negative (semi-)definite. v is obtained by integration, where, as in Kinnen and Chen 1968 [136], the condition is imposed that the integration of v is independent of the path.
In the 1980's, a series of publications presented construction methods for Lyapunov functions based on Zubov's equation. Vannelli and Vidyasagar 1985 [223] consider Zubov's equation (10) and construct a solution with the rational ansatz where R i and Q i are homogeneous functions of degree i. A rational function is chosen to model the fact that v(x) tends to infinity as x approaches the boundary of the basin of attraction. A recursive algorithm is presented to compute the function and numerical examples are given.
Hassan and Storey 1981 [117] use Zubov's equation to calculate the inverse function x(v), for which the equation The equation is solved numerically, e.g. by Runge-Kutta methods, using initial points along a radial line (in two dimensions) from the origin to approximate the points x(v), where v is close to ∞ in the above form. This gives an estimate of the boundary of the basin of attraction.
Genesio, Targaglia, and Vicino 1985 [81] give a survey of the methods to estimate basins of attraction, in particular from the view of the engineering community. They organize the methods into two main classes: Zubov and La Salle, where the latter is a reference to a generalized version of Lyapunov functions: A (non-strict) Lyapunov function for an equilibrium still allows the same conclusions on the basin of attraction if no trajectory (apart from the equilibrium) lies entirely in the region {x :V (x) = 0}. Most of the discussed methods apply to systems with a specific structure. The original contribution of the paper is to start with a small neighborhood of (each) equilibrium, and then use backward integration to extend the respective basin of attraction.
Grujić  [51] and Grüne 2002 [105]; it is worth noting that it includes as a special case the unperturbed, autonomous ODE.
In more detail, the authors consider the systeṁ where u : R + 0 → U is an arbitrary, time-varying, measurable function and U ⊂ R m is a compact set of perturbations, cf. Section 2.5. Zubov's equation in this case generalizes to the following PDE where V approaches 1 if x approaches the uniform robust domain of attraction. The general idea to solve this equation numerically is to use a first order numerical scheme from Falcone 1997 [72]. A bounded domain Ω ⊂ R n is equipped with a simplicial grid with edges x i , and the approximation V solves where h > 0 is the time step. However, as the equilibrium is a singularity, the convergence fails. Hence, the authors regularize the equation by replacing ϕ by the strictly positive function ϕ (x, u) = max{ϕ(x, u), } for a small > 0. This method was generalized to other systems, including control systems Camilli, Grüne, and Wirth 2008 [52]. Zubov-based construction method for the design of controls were also considered in Dubljević and Kazantzis 2002 [68].
As Lyapunov functions can be described as solution of a PDE, one can use collocation methods to find an approximate solution to the PDE, which turns out to be a Lyapunov function. A collocation method in general is a numerical method to solve ODEs or PDEs by finding an approximate solution in a fixed finitedimensional function space, which satisfies the differential equation on finitely many so-called collocation points. Giesl 2007 [84] considers Zubov's equation in the form (10) and solves it numerically using mesh-free collocation, in particular using Radial Basis Functions (RBF). This is a general collocation method to solve linear PDEs, see e.g. Buhmann 2003 [47] and Wendland 2005 [225]. In particular, no triangulation of the state space is needed, i.e. collocation points can be scattered, and the approximation is a smooth function.
Given a set of N scattered collocation points X = {x 1 , . . . , x N } ⊂ Ω, where Ω ⊂ R n is a domain, the ansatzṽ is given bỹ where Lv =v is the operator of the orbital derivative, the superscript y denotes the application of the operator with respect to y, δ x k v(x) = v(x k ) denotes Dirac's delta distribution and Φ is an RBF, or more generally a kernel of a Reproducing Kernel Hilbert Space (RKHS) H. Note that among all functions in the RKHS, which satisfy the PDE at the collocation points,ṽ is the function with the minimal RKHS norm. The coefficient vector α ∈ R N in (11) is determined as the solution of a system of linear equations. Error estimates of the form have been shown Giesl and Wendland 2007 [96], where h is the fill distance, measuring how dense the collocation points X lie in Ω, and k is related to the smoothness of both f and Φ. Estimates on the level sets ofṽ are also available, showing that each compact subset of the basin of attraction can be covered by a sublevel set ofṽ. Thus, apart from an arbitrarily small neighborhood of the equilibrium, the computed functions itself is a Lyapunov function. The neighborhood can either be covered by standard estimates, or by a modification of the approach Giesl 2008 [86]. Refinement of the grid is discussed in Mohammed and Giesl 2015 [176]. The method has also been applied to time-periodic ODEs Giesl and Wendland 2009 [97], discrete-time dynamical systems Giesl 2007 [85], asymptotically autonomous systems Giesl and Wendland 2012 [98], and finite-time Lyapunov functions Giesl 2012 [88]. A construction method for a contraction metric for periodic orbits in autonomous ODEs using RBF was given in Giesl 2009 [87]. A different approach using RKHS to construct a control is used in Bouvrie and Hamzi 2011 [41]; the method assumes that the nonlinear system behaves linearly when lifted to the RKHS. Rezaiee-Pajand and Moghaddasie 2012 [200] proposed a collocation method using two classes of basis functions in Cartesian and polar coordinates.

LP based
Methods. An LP (Linear Programming) problem is a problem that has the following structure minimize c T x subject to x ≥ 0 and Ax ≤ 0.
Here c, x ∈ R s and A ∈ R r×s , where often r is much larger than s, and the inequalities are to be understood componentwise, e.g. x ≥ y means x i ≥ y i for all i. LP problems form a subclass of convex optimization problems and can be solved in polynomial time in the number of the variables and the number of the constraints, i.e. the number of the inequalities. There is a very rich literature on convex optimization, e.g. the book Boyd and Vandenberghe 2004 [43] is a good reference for convex optimization in control. Numerous software packages to solve various kinds of convex optimization problems are available, e.g. CPLEX, Gurobi, or GLPK to name a few. Brayton and Tong 1979 [45] published an influential paper where they proposed using LP to prove the stability of equilibria of dynamical systems. Their idea was the following: Given a set of n × n-matrices A = {M 0 , M 1 , . . . , M m−1 } and a matrix norm · W such that M i W ≤ 1 for all M i ∈ A, the discrete-time system , must be stable. In a modern presentation, they were studying the stability of the arbitrarily switched system x k+1 = M σ x k . Their constructive approach was to start with an initial bounded polyhedral (alias polytopic) neighborhood W 0 of the origin. They then proceeded to compute iteratively by LP the neighborhoods They proved that if W := ∞ k=0 W k is bounded, then it is convex and can be taken to be symmetric, i.e. x ∈ W ⇒ −x ∈ W , and thus its Minkowskii functional defines a vector norm on R n .
The vector norm · W then serves as a (non-strict) Lyapunov function for the dynamical system. They also used this method to study continuous-time dynamical systems by considering their discretization. In Brayton and Tong 1980 [46] they extended and improved their original results in several directions. Especially, they considered asymptotic stability and a set of matrices A compact in the R n×n topology, that is not necessarily finite. In this case the induced matrix norm fulfills M W < 1 for all M ∈ A.
By results in Filippov 1988 [76], the asymptotic stability of the origin for the continuous-time switched system under arbitrary switchinġ where A is any compact set of matrices, is equivalent to the asymptotic stability of the origin for the differential inclusioṅ The same holds for the discrete-time switched system and the difference inclusion As shown by Pyatnitskiȋ 1986, 1989 [178, 177], the asymptotic stability of both (14) and (13) is equivalent to the existence of a continuous, piecewise quadratic Lyapunov function of the form where the vectors l i are used to parameterize V . Since the square-root of a Lyapunov function is also a Lyapunov function, V (x) = max i=1,2,...,m |l i · x| is a Lyapunov function for the same system and by defining the m × n-matrix P by writing the l i s in its rows, V can be written as a weighted infinity norm Lyapunov functions of the form (15) and (16) are often referred to as piecewise quadratic and piecewise linear Lyapunov functions, respectively. Some authors, like Brayton and Tong whose approach was discussed above, use LP to compute convex and compact sets, of which the Minkowskii functional serves as a Lyapunov function. These Lyapunov functions are somewhat more general than piecewise linear Lyapunov function, but if the convex set C is a symmetric polyhedron, then its Minkowskii functional can be written as a weighted infinity norm and thus is a piecewise linear Lyapunov function, cf. Blanchini and Carabelli 1994 [37].
There is a vast number of publications concerned with the computation of piecewise quadratic and piecewise linear Lyapunov functions for various kinds of systems. Most usually, piecewise linear Lyapunov functions are computed using LP and piecewise quadratic Lyapunov functions are computed by solving LMI, see also Section 3.3. Some generalizations and extensions of the results by Brayton and Tong are: Michel, Miller, and Nam 1982 [172] and Michel, Nam, and Vittal 1984 [173] considered the applicability of the method of Brayton and Tong to interconnected continuous-time systems. Wang and Michel 1996 [224] advanced the theory of Brayton and Tong and reduced its computational complexity. Erickson and Michel 1985 [70, 71] utilize the method of Brayton and Tong in the stability analysis of fixed-point digital filters.
Bernussou and Peres 1989 [28] considered the computation of a quadratic Lyapunov function for differential inclusions. They solved an LMI by iteratively solving LP problems and using the cutting plane technique. Blanchini 1991Blanchini , 1994 showed that by using similar ideas as Brayton and Tong he could verify ultimate boundedness (practical stability) solving an LP problem to compute a so-called set-induced Lyapunov function, which is the Minkowskii functional of a convex, compact set. For symmetric polyhedral sets, set-induced Lyapunov functions are representable as weighted infinity norm Lyapunov functions as discussed above. In Blanchini 1995 [36] he extended his ideas to the robust control of uncertain linear systems. Polanski 1997Polanski , 2000, Ohta 2001 [181], Ohta and Tsuji 2003 [182], and Yfoulis and Shorten 2004 [227] considered the parameterization of piecewise linear Lyapunov functions for differential inclusions by LP. They were especially interested in how to compute an appropriate conic partition of the state space for their computations.
Blanchini and Miani published a book 2008 [38], where they treat polyhedral Lyapunov functions. Lazar 2010 [142] and Lazar and Doban 2011 [143] studied piecewise linear Lyapunov function, also for quadratic dynamics. Lazar and Jokić 2010 [145] discussed the control of switched systems with state dependent switching of linear systems defined on cones using LP and Lazar, Doban, and Athanasopoulos 2013 [144] treat discrete-time systems of the same kind.
Around 2000, approaches using LP to parameterize Lyapunov functions for systemsẋ = f (x), where f is a general nonlinear function that is not supposed to be of a specific algebraic structure, appeared in the literature. As one of the first of such approaches, Johansen 2000 [122] used LP to parameterize partitions of unity to weight multiple quadratic Lyapunov functions for such systems. By using more general convex optimization, his approach can be extended to parameter-dependent systems.
Julian, Guivant, and Desages 1999 [126] and Julian in his PhD thesis 1999 [125], presented an LP problem to parameterize Lyapunov functions that are continuous and affine on each simplex of a triangulation of the state space. Henceforth we refer to such functions as CPA (continuous piecewise affine) functions and the computation of CPA Lyapunov function as the CPA method. The values of a CPA function V : R n → R are determined on a simplex S ν by its values V ν,i = V (x ν,i ) at the simplex vertices x ν,i , i = 0, 1, . . . , n. Further, the gradient ∇V ν of V inside of S ν , a constant vector, is a linear function of the V ν,i . In its simplest form, an LP problem to compute a CPA Lyapunov function approximately for a systemẋ = f (x) on a triangulation (S ν ) ν of a compact subset of R n can be formulated as for all vertices x ν,i of all simplices S ν of the triangulation.
Hafstein (aka Marinósson) 2002 [163] and in his PhD thesis [164] included error terms in the LP problem, such that the computed CPA function is guaranteed to fulfill the conditions for a Lyapunov function exactly. The idea is to first compute constants E ν,i , that depend on the geometry of the simplices and upper bounds on the second derivatives of the components of f . The LP problem then becomes for all vertices x ν,i of all simplices S ν of the triangulation. As the 1-norm ∇V ν 1 = n i=1 |(∇V ν ) i | can be implemented as linear constraints, this results in an LP problem. An essential result to prove that the resulting function V is a Lyapunov function for the system is [164,Lemma 4.16], which in a more general presentation [19] states that for every x = n i=0 λ i x ν,i ∈ S ν we have Hafstein 2004Hafstein , 2005 proved that for nonlinear systems with an exponentially stable or asymptotically stable equilibrium this method will always result in the construction of a Lyapunov function if the triangulation is fine enough and an arbitrary small neighborhood of the origin is excised from the domain of the Lyapunov function. If such a neighborhood is excised, then the resulting Lyapunov function secures practical stability of the equilibrium, but not necessarily asymptotic stability. Hafstein 2012, 2014 [91, 94] proved that by using a fan-like triangulation at the equilibrium, a revised CPA method can compute a CPA Lyapunov function for any nonlinear system with exponentially stable equilibrium, fulfilling all conditions of a Lyapunov function, also in a neighborhood of the equilibrium. Especially, a system with f ∈ C 2 and an exponentially stable equilibrium at the origin possesses a local conewise linear Lyapunov function.

PETER GIESL AND SIGURDUR HAFSTEIN
The CPA method has been extended to more general systems: in Hafstein 2007 [110] to arbitrarily switched, non-autonomous systems, in Baier, Grüne, and Hafstein 2010 [19] to differential inclusions, in Giesl and Hafstein 2014 [93] to discretetime systems, and in Li, Baier, Hafstein, Grüne, and Wirth 2014, 2015 [148,149] to ISS Lyapunov functions. Baier and Hafstein 2014 [20] compute control Lyapunov functions, but the optimization problem that has to be solved is in this case a mixed-integer problem, a problem that is much more computationally demanding to solve. Further, the CPA method was adapted to compute CPA contraction metrics for nonlinear systems with periodic orbits in Giesl and Hafstein 2013 [92]; the resulting optimization problem is an LMI. LMIs are, just as LP problems, a subclass of convex optimization problems that can be solved effectively. Methods to compute Lyapunov functions using LMIs are the subject of next section.

LMI (Linear Matrix Inequality) based methods.
Recall that for linear systems, both continuous-time and discrete-time, the asymptotic stability of the origin is equivalent to its exponential stability. The origin is an asymptotically stable equilibrium of the linear systemẋ = Ax, if and only if there is a symmetric, positive definite matrix P such that P A + A T P is a negative definite matrix. In that case V (x) = x T P x is a quadratic Lyapunov function for the system. Recall, that for a symmetric P ∈ R n×n we write P 0 if P is strictly positive definite and P ≺ 0 if P is strictly negative definite. For symmetric P, Q ∈ R n×n we write P Q if P − Q 0 and P ≺ Q if P − Q ≺ 0.
The statement above can thus be written: The origin is asymptotically stable, if and only if there exists a P 0 such that P A + A T P ≺ 0. The inequality P A + A T P ≺ 0 is called the continuous-time Lyapunov inequality and we can reformulate the problem once again: the origin is asymptotically stable, if and only if the continuous-time Lyapunov inequality possesses a positive definite solution P .
In general, a problem of the form: find an x T = (x 1 , x 2 , . . . , x m ) ∈ R m such that is called a (strict) LMI problem. The problem of finding a solution P 0 to the Lyapunov inequality can be written as an LMI problem as follows: LetF i ∈ R n×n , i = 1, 2, . . . , m, be a basis of the vector space of symmetric n × n-matrices. Define  Bartels and Stewart 1972 [23], but are still a subject of active research, cf. e.g. the PhD thesis by Mikkelsen 2009 [175], but then for very large systems.
Slight generalizations of the problem, however, lead to LMI problems. As a simple example we can consider the differential inclusion of linear systemsẋ ∈ co{A i x : i ∈ J}, where J is some index set. To prove the asymptotic stability of the origin for this system one could search for a quadratic Lyapunov function V (x) = x T P x, where P 0 and P A T i + A i P ≺ 0 for all i ∈ J. This approach, however, is rather restrictive because there are many differential inclusions with an asymptotically stable equilibrium at the origin, that do not possess a quadratic Lyapunov function. Such an example is, e.g., given in Dayawansa and Martin 1999 [64].
For discrete-time, linear systems very similar considerations apply with the continuous-time Lyapunov inequality replaced by the discrete-time Lyapunov inequality We refer the interested reader to the book "Linear Matrix Inequalities in System and Control Theory" by Boyd, Ghaoui, Feron, and Balakrishnan 1997 [42]. Closely related to LMIs in stability theory is the so-called S-procedure. The S-procedure gives sufficient conditions on when one quadratic (in)equality is a consequence of another quadratic (in)equality. We refer the interested reader to the review paper Pólik and Terlaky 2007 [193] and let a simple example suffice to clarify its use. Consider the systeṁ One can think of the constant γ as a quantitative measure of the nonlinearity of f in some neighborhood of the origin. It is not difficult to see that V (x) = x T P x is a quadratic Lyapunov function for the system if P 0 and there exists a constant α > 0 such that for all g(x) ≤ γ x . The S-procedure now tells us that z T z − γ 2 x T x ≤ 0 implies where I n is the n × n identity matrix. The last inequality and the condition P 0 can now be written as an LMI problem similar to above, and its solution delivers a Lyapunov function. In this simple case, the S-procedure gives not only a sufficient, but also a necessary condition. In general, however, this is not true. The development of methods using LMI to compute Lyapunov functions is in some sense similar to the development of the methods using LP. Differential and difference inclusions of linear systems and other systems of specific algebraic structure, for example piecewise affine, have been considered for a long time. One can, for example, attempt to construct a piecewise quadratic Lyapunov function of the form (15) as discussed in the last section, see, for example Goebel, Teel, Hu, and Lin 2006 [100] or Ambrusino and Garone 2006 [8]. Further good references for such methods are the book Boyd, Ghaoui, Feron, and Balakrishnan 1994 [42], the review Shorten, Wirth, Mason, Wulff, and King 2007 [206] and the review Sun 2010 [217]. Another use of LMIs is to compute a good or even optimal basin of attraction given a quadratic Lyapunov function for a polynomial system as shown in Tesi, Villoresi, and Genesio 1994 [222].
In the mid 1990's studies of Branicky 1998 [44], Johansson and Rantzer 1998 [124], and Johansson in his PhD Thesis 1999 [123] on the stability of systems with state dependent switching received much attention. Among other things, they combined piecewise quadratic Lyapunov functions and the S-procedure with the continuous-time Lyapunov inequality.
To explain some of the more advanced ideas in LMI programming for the construction of Lyapunov functions we consider the following example from Johansson and Rantzer 1998 [124]: Set The origin is easily shown to be asymptotically stable for the systemsẋ = A i x for every i = 1, 2, 3, 4, however, the origin is an unstable equilibrium of the inclusion systemẋ ∈ co{A 1 , A 2 }x and there is thus no matrix P fulfilling P 0, P A T 1 + A 1 P ≺ 0, and P A T 2 + A 2 P ≺ 0. Consider the matrices i.e. the conical sectors E i are characterized by inequalities E i x ≥ 0. It is easy to verify that for every 2 × 2 matrix P the function For V to be a Lyapunov function for the inclusionẋ ∈ co{A i x : x ∈ E i } we need V to be positive definite and its orbital derivative to be negative definite. Thus, if we can find a matrix P such that x < 0 we are done. Let U 1 and W 1 be 2 × 2 matrices with non-negative entries. Then Thus, if we can find symmetric matrices P , U i and W i , i = 1, 2, 3, 4, such that the U i and W i have non-negative entries and such that the LMIs fulfilled, then V is a Lyapunov function. Note that it might be easier to fulfill is a Lyapunov function for the systemẋ ∈ co{A i x : x ∈ E i }.
One of the challenges of computing a piecewise quadratic Lyapunov function is to secure the continuity of the computed Lyapunov function, a non-trivial task for more complicated systems. This has led to the consideration of discontinuous piecewise quadratic Lyapunov functions. A recent publication on discontinuous Lyapunov functions is, for example, Eghbal, Pariz, and Karimpour 2013 [69].
From the year 2000 much effort has been put into LMI methods for general nonlinear systems. This was initiated by an influential PhD thesis by Parrilo 2000 [186], where he introduced the sum-of-squares relaxation for non-negative polynomials. The essential idea is as follows: Consider a systemẋ = f (x), where f : R n → R n is a polynomial in x. If V (x) is a polynomial function of degree m, we can write where a α ∈ R, α = (α 1 , α 2 , . . . , α n ) ∈ N n 0 is a multi-index, |α| := n i=1 α i is its length, and x α := x α1 1 x α2 2 · · · x αn n . Similarlẏ where M is the degree ofV (x) and the coefficients b α are linear functions of the coefficients a α . Now, V (x) is a (global) Lyapunov function for the system, if and only if V (x) and −V (x) are positive definite functions. We can also proceed by computing the coefficients a α , such that V (x) and −V (x) are positive definite. In general this is an NP-hard problem. Parrilo's idea was to relax this condition to computing the coefficients a α such that V (x) and −V (x) are the sum of squared (SOS) polynomials. A polynomial P (x) is said to be an SOS polynomial, if there are polynomials p i such that Obviously P (x) ≥ 0 for all x. Further, if we define a vector Z, whose entries are all monomials of the form x α , |α| ≤ N , then every SOS polynomial P (x) of degree 2N can be written as P (x) = Z T QZ, where Q is a positive definite matrix. As a clarifying example consider P (x 1 , x 2 ) of degree 4. Then P (x 1 , x 2 ) is an SOS 28 PETER GIESL AND SIGURDUR HAFSTEIN polynomial, if there is a Q 0 such that Just note that since Q 0, it can be factorized as Q = O T DO, where D = diag(d 1 , d 2 , . . . , d 6 ) is a diagonal matrix with non-negative entries and O is an orthogonal matrix. The entries p i = p i (x 1 , x 2 ) of the vector OZ = (p 1 , p 2 , . . . , p 6 ) T are polynomials in (x 1 , x 2 ) and Parrilo proposed an LMI program to search efficiently for Lyapunov functions that are SOS polynomials for polynomial systems. In the simplified setting above one could only find global Lyapunov functions, but by use of the Positivstellensatz it is possible to search for SOS Lyapunov functions on a compact domain. Peet and Papachristodoulou 2007 [189] proved that a polynomial system with an exponentially stable equilibrium possesses an SOS polynomial Lyapunov function for any compact subset of the basin of attraction. Note that even for polynomial systems there might not exist a global polynomial Lyapunov function, cf. e.g. Ahmadi, Krstic, and Parillo 2011 [4]. Peet 2007 [188] proved the existence of polynomial Lyapunov functions for systems possessing an exponentially stable equilibrium and with a smooth enough right-hand side. There is a vast literature on the SOS method to compute Lyapunov functions in various settings and for different kinds of systems. We refer the reader to Chesi 2010 [55] and Anderson and Papachristodoulou 2015 [9] for an overview. There is also a free toolbox for Matlab under the GNU General Public License called SOSTOOLS, meanwhile available in the version 3.00, to compute SOS Lyapunov functions. For further information see its user's guide Papachristodoulou, Anderson, Valmorbida, Pranja, Seiler, and Parrilo 2013 [185].
A construction method using SOS for a contraction metric to show global stability of an equilibrium of a polynomial system was presented in Aylward, Parrilo, and Slotine 2008 [18]. In Manchester and Slotine 2014 [160], SOS was used to construct a contraction metric for a periodic orbit of an autonomous system. There have been many other different approaches that use LMIs to compute Lyapunov functions. We end this section by listing a few of them: Pettersson and Lennartson 1996 [190] used LMI combined with the S-procedure to compute Lyapunov functions for hybrid systems. Hu  Optimal quadratic Lyapunov functions. Attempts have been made to maximize the estimate of the basin of attraction secured by a quadratic Lyapunov function for a nonlinear system, i.e. to compute the best quadratic Lyapunov function in this sense. The first such publication is Davison and Kurak 1971 [62], where they minimize the product of the eigenvalues of P 0 in V (x) = x T P x under the constraints max x T P x≤1V (x) < 0. In Loparo and Blankenship 1978 [157] the results were improved for systems with analytic right-hand sides. Michel, Sarabudla, and Miller 1982 [174] presented more effective algorithms for such computations. Panikhom and Sujitjorn 2012 [184] used adaptive tabu-search for the same purpose.
Neural networks and genetic algorithms. There have been several proposals of using neural networks to compute Lyapunov functions, for example by Prokhorov 1994 [195], Serpen 2005 [203], and Noroozi, Karimaghaee, Safaei, and Javadi 2008 [179]. Genetic algorithms have been applied, for example, by Grosman and Lewin 2009 [102]. Candidate Lyapunov functions receive a score based on the size of the subset of the basin of attraction. The probability of a Lyapunov function surviving into the next generation depends on this score. In the PhD thesis Hargrave 2008 [114], genetic algorithms are used to simultaneously tune the control law and the control Lyapunov function.
Simulation methods. Kapinski, Deshmukh, Sankaranarayanan, and Arechiga 2014 [130] suggested assuming a candidate Lyapunov function in a certain parameterized template form, a SOS polynomial of a fixed degree. They use linear programming to parameterize the candidate Lyapunov function, combined with a subsequent search of points in the domain of interest where the conditions of a Lyapunov function are violated. Such points are then used to add constraints to the the linear programming problem. This is done iteratively until their stochastic global optimizer is not able to find any further points where the conditions are violated.
Another recent simulation method was given in Menck, Heitzig, Marwan, and Kurths 2013 [169]. They do not use Lyapunov functions directly, but estimate the size of the basin of attraction of equilibria for nonlinear systems of very high dimension (> 300) using Monte-Carlo simulations.
Polynomial methods. There have been numerous proposals on the construction of polynomial Lyapunov functions outside the SOS framework. To name a few, Burchardt and Ratschan 2007 [48] used quantifier elimination and Ratschan and She 2010 [199] used a branch-and-relax algorithm to compute polynomial Lyapunov functions securing practical stability for polynomial systems. She, Li, Xue, Zheng, and Xia 2013 [204] proposed a real root classification based approach to compute polynomial Lyapunov functions for polynomial systems. Groebner bases, the analogue of Gaussian elimination to systems of polynomial equations, were used in Forsman 1991 [79]. For more information on this subject the reader is referred to Kamyar and Peet 2015 [129], which includes a discussion of the advantages and disadvantages of methods as different as quantifier elimination, reformulation linear techniques, blossoming, Groebner basis methods, algorithms defined by Polya's theorem, Bernstein's theorem, and Handelsman's theorem.
Set-oriented methods. The earliest application of set-oriented methods to study dynamical systems is given in Hsu 1987 [118] and Flashner 1988 [77]; here, the basin of attraction is computed by starting with a small neighborhood of the equilibrium and numerically solving the time-reversed equationẋ = −f (x).
In Grüne 2002 [105], set-oriented methods and a subdivision algorithm are used to construct Lyapunov functions. In this algorithm, a discrete-time system (possibly derived from a continuous-time one) is considered. Grüne considered an equilibrium (or more generally attractor), a small neighborhood S of the equilibrium, the socalled target set, and a compact set Ω, containing the basin of attraction. Ω is divided into cells and each cell has a status, e.g. in, out, or partially in according to whether or not the cell lies in the basin of attraction. The algorithm checks in each step whether a cell is mapped into a cell that is already identified as being in the basin of attraction; also, the subdivision has to be refined in each step. These checks can be made rigorous, using rigorous subdivision techniques, see Junge 1999 [127]. For an overview of more set-oriented numerical methods for dynamical systems see Dellnitz and Junge 2002 [66] and for a numerical software package see GAIO [65].
Koltai 2011 [138] used a set-oriented method to compute the basin of attraction, however, he does not use trajectory simulations. Instead, after dividing the area under consideration into boxes, he defines a finite state continuous-time Markov process (Markov jump process) which describes the flow between them. The absorption probabilities of the Markov jump process give an indication for the basin of attraction and its generator is computed by integration over box faces rather than by computing solution trajectories.
Kalies, Mischaikow, and VanderVorst 2005 [128] presented a set-theoretic method to compute complete Lyapunov functions. This method has been implemented in Ban and Kalies 2006 [21] and Björnsson, Giesl, and Hafstein 2014 [31], where it was combined with the CPA method, verifying the inequalities, and compared with a construction based on the RBF method. A thoroughly worked out efficient algorithm for the computation of such Lyapunov functions is presented in Goullet, Harker, Mischaikow, Kalies, and Kasti 2015 [101]. In this context, the international 'Computational Homology Project' (see http://chomp.rutgers.edu) of Mischaikow, Kokubu, Mrozek, Pilarczyk, Gedeon, and Lessard is worth mentioning.
Massera and Yoshizawa construction. The converse theorems in the Lyapunov theory, i.e. the theorems that assert the existence of a Lyapunov function, are usually proved by constructing the Lyapunov function explicitly from the solution trajectories. For a discrete-time, autonomous systems x k+1 = g(x k ), where g is continuous and the origin is a locally exponentially stable equilibrium, one can prove, for example, that for large enough N ∈ N V (ξ) = N k=1 x(k; ξ) 2 2 is a Lyapunov function for the system on a neighborhood of the origin; here, we denote the solution with initial condition x 0 = ξ by x(k; ξ) = g •k (ξ). We refer to such constructions as Massera constructions, cf. Massera 1949 [166]. For a continuous-time system an analogous construction is Another possibility is to use a Yoshizawa construction, see Yoshizawa 1966 [228], e.g. that for an appropriate α ∈ K ∞ the function V (ξ) := sup α( x(t; ξ) )e t is a Lyapunov function. Both can be used to compute the values of a Lyapunov function at a finite number of points. In a recent publication by Geiselhart, Gielen, Lazar, and Wirth 2014 [80] this idea is taken a step further for discrete-time systems x k+1 = g(x k ) by first constructing a so-called finite time Lyapunov function W fulfilling W (x(M ; ξ)) ≤ ρ(W (ξ)), where M ∈ N and ρ(r) < r for all r ∈ R + , and then a Lyapunov function V (ξ) = M −1 i=0 W (x(i; ξ)). In this setting g does not have to be continuous and they were able to prove that for a conewise linear system the exponential stability of the equilibrium at the origin is equivalent to the existence of a conewise linear Lyapunov function for the system, a previously open problem. Note that here the term finite time Lyapunov function denotes a different concept than Lyapunov function for non-autonomous systems on finite time intervals as in e.g. [88,95] 4. Concluding Remarks. Lyapunov functions, introduced well over a hundred years ago, are to this day an essential tool in the stability analysis of dynamical systems, both in theory and applications in science and engineering. Their construction is an ongoing problem in different communities, ranging from real world applications in engineering to theoretical mathematics. The methods to tackle the construction problem come from diverse areas and the dynamical systems considered vary considerably.
In this review we have brought together these different methods and described the state of the art of the vast variety of methods to compute Lyapunov functions for various kinds of systems. A main goal was in particular to bring the different scientific communities, usually publishing in separate journals, together in one review.
It will be interesting to see how these methods will evolve and how they will be generalized over the next decades, and which new methods will emerge. The simple, but powerful idea behind Lyapunov functions makes them still an up-to-date method to study dynamical systems in theory and practice, and their construction remains a challenging task, which is and will continue to be a very active area of research.