{'fraction_non_alphanumeric': 0.049992117915976984, 'fraction_numerical': 0.012946323007803262, 'mean_word_length': 4.170028524857376, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cone-Copositive Lyapunov Functions for Complementarity Systems: Converse Result and Polynomial Approximation

This article establishes the existence of Lyapunov functions for analyzing the stability of a class of state-constrained systems, and it describes algorithms for their numerical computation. The system model consists of a differential equation coupled with a set-valued relation which introduces discontinuities in the vector field at the boundaries of the constraint set. In particular, the set-valued relation is described by the subdifferential of the indicator function of a closed convex cone, which results in a cone-complementarity system. The question of analyzing stability of such systems is addressed by constructing cone-copositive Lyapunov functions. As a first analytical result, we show that exponentially stable complementarity systems always admit a continuously differentiable cone-copositive Lyapunov function. Putting some more structure on the system vector field, such as homogeneity, we can show that the aforementioned functions can be approximated by a rational function of cone-copositive homogeneous polynomials. This later class of functions is seen to be particularly amenable for numerical computation as we provide two types of algorithms for precisely that purpose. These algorithms consist of a hierarchy of either linear or semidefinite optimization problems for computing the desired cone-copositive Lyapunov function. Some examples are given to illustrate our approach.


I. INTRODUCTION
Lyapunov functions provide a useful tool for the stability analysis of dynamical systems. Several advances have been made on the theoretical side to establish existence of Lyapunov functions for various classes of dynamical systems, see e.g. [30], [33], [23] for examples of standard expositions. The fundamental question in most of these works boils down to checking the positivity of certain functions over the state space, which is a challenging problem numerically [36]. Modern developments in the field of real algebraic geometry [46], [48] provide certificates of positivity of (polynomial) functions with Positivstellensätze relying on sums-of-squares (SOS) decompositions. Since it has been observed in [43] that finding SOS decompositions is equivalent to semidefinite programming (SDP) or linear matrix inequalities (LMI), numerical tools based on SOS optimization have been developed extensively over the past two decades to compute Lyapunov functions, see e.g. [40], [45], [26], [15].
Stability analysis of hybrid, or nonsmooth dynamical systems, where the vector field is set-valued with possible discontinuities, is of particular relevance with respect to several applications. Naturally, Lyapunov functions for such systems provide a potent tool for studying stability related properties as well. When the system is modeled by switching vector fields over the whole state space, then the construction of Lyapunov functions using SOS is studied in [39], [3], [1]. However, we are concerned with a certain class of differential inclusions which is useful in modeling systems with state constraints, where the vector field exhibits discontinuous behaviour on the boundary of the constraints so that the state trajectory is forced to evolve within the prespecified set. In the literature, there are several frameworks for modeling this behaviour, such as sweeping processes, projected dynamical systems or complementarity systems [7], [8], [13]. The relevance of these systems is seen in many practical systems encountered in engineering, physics and biology. For example, in mechanics, the interaction between multiple rigid bodies or between a rigid body and the environment can be modelled using nonsmooth force laws for contact, impact and friction. Some variants of these systems are also studied in [49], [50] in the context of control-theoretic problems.
The stability analysis of complementarity systems using Lyapunov functions has received some attention in the literature. Since the state of such systems essentially evolves in a closed convex cone, often chosen to be the positive orthant, it is naturally desirable to consider Lyapunov functions which are positive definite over the positive orthant; the functions satisfying this latter property are called copositive functions. The need to search such functions for stability analysis of complementarity systems was presented as an open problem in [14]. The papers [25], [24], [13] investigate sufficient stability conditions for linear complementarity systems, or conewise linear systems [28] in terms of copositive Lyapunov functions. The paper [24] also provides examples of systems where a positive definite Lyapunov function does not exist, but the system is nonetheless asymptotically stable and it admits a copositive Lyapunov function.
While these existing works have shown the utility of enlarging the search space of Lyapunov functions from positive definite to copositive functions, and cone-copositive functions when considering systems with state trajectories constrained to a cone rather than the positive orthant, none of the existing works has addressed the converse question: Does every asymptotically stable complementarity system admit a cone-copositive Lyapunov function?
The first objective of this paper is to answer this question in the affirmative by constructing a Lyapunov function as a functional of the solution trajectories, thereby concluding that one does not need to go beyond cone-copositive functions to find Lyapunov functions for complementarity systems. By putting more structure on the system dynamics, and using the appropriate density results, we are able to prove the existence of a cone-copositive Lyapunov function which can be expressed as a ratio of homogeneous polynomials. Converse stability results for dynamical systems have been studied for a long time in control community, see the recent survey article [29]. Moreover, due to discontinuities in the vector field at the boundary of the constraint set (which can be seen as an example of constrained switching), establishing the existence of Lyapunov functions within cone-copositive functions becomes difficult.
The second objective of this paper is to propose computationally tractable algorithms for finding the Lyapunov functions. The interesting aspect of our problem lies in computing Lyapunov functions which satisfy certain inequalities over a given set. For example, in linear complementarity systems, one needs to check the positivity of a function over the positive orthant only, and if the function we seek is of the form x P x, then finding such a function boils down to finding a copositive matrix P that satisfies certain inequalities. However, checking whether a given matrix is copositive is an NP-hard problem [5]. The papers [10], [11], [37], [21] propose algorithms for detecting copositivity of a matrix or tensor. Moreover, we will show with the help of an example that, even in the case of linear complementarity systems, such functions cannot be computed by solving a linear set of equations, as is done for unconstrained linear systems. Another challenging aspect of these problems is that, when dealing with conic constraints which are unbounded sets, there are no readily available Positivstellensatz that guarantee SOS decompositions of a positive polynomial over the sets of our interest. The field of copositive programming has been active area of research over the past decade which addresses some of these challenges [4]. In computing the Lyapunov functions for complementarity systems which evolve on unbounded cones with positivity constraints, we are faced with similar challenges.
Motivated by such questions, we propose two approaches for computing homogeneous cone-copositive Lyapunov functions numerically. The first one is a discretization method which is based on finding an inner approximation of the cone of cone-copositive polynomials by using simplicial partitions and evaluating inequalities over a set of points taken on the simplex. It is shown that, as the partition gets finer, we can approximate any cone-copositive polynomial function. The second approach is an SOS method where we show that the positivity of polynomial over the given cone can be checked by expressing it as an SOS function. By increasing the degree of the approximating SOS polynomial, we again obtain a hierarchy of SDP problems to compute the desired Lyapunov function. Then, we derive the corresponding algorithms for those two techniques, which can be seen as an adaptation of tools available in the literature on polynomial optimization. The illustration of some academic examples is provided using standard Matlab toolboxes.

II. SYSTEM CLASS
We begin this section by introducing some basic notions from convex analysis which will be used for describing the class of dynamical systems studied in this paper.

A. Dynamical System with Constrained Trajectories
We are interested in studying a class of dynamical systems described by the variational inequalitieṡ where f : R n → R n is a given vector field, x(t) ∈ R n denotes the state, S is a given closed convex subset of R n containing the origin. The normal cone to S at x is defined by If x ∈ int(S), the interior of S, then N S (x) = 0 and by convention, we let N S (x) := ∅ for all x ∈ S. The differential inclusion (1) is a particular case of the following class of variational inequalitieṡ where ϕ : R n → R ∪ {+∞} is a proper, convex and lower semicontinuous function, and the subdifferential of ϕ at x ∈ R n is defined as ∂ϕ(x) := {λ ∈ R n | λ, z − x ϕ(z) − ϕ(x), ∀z ∈ dom(ϕ)}, with dom(ϕ) := {x ∈ R n | ϕ(x) < +∞}. Indeed, if one denotes the indicator function of a closed convex set S ⊂ R n by ψ S (·), that is, ψ S (x) = 0, if x ∈ S and ψ S (x) = +∞ if x ∈ S, then (1) is obtained from (3) by choosing ϕ = ψ S . Inclusion (1) captures the class of complementarity systems studied in this paper, but the framework of (3) is necessary for a broader class of complementarity systems such as the ones studied in [13], [50]. The formalism of system (1) with inclusion naturally allows us to describe dynamics constrained to evolve in set S. Using the depiction in Figure 1, it is seen that, during the evolution of a trajectory, if x(t) is in interior of S, then N S (x(t)) = 0 and the motion of the trajectory continues according to the differential equationẋ(t) = f (x(t)). While x(t) is at the boundary, we add a vector from the set −N S (x(t)), which restricts the motion of the state trajectory in tangential direction on the boundary of the constraint set S.

B. Complementarity Systems
In this article, we focus on the particular class of constrained systems where the admissible set S is a cone, denoted by K. Hence for each x ∈ K, we have λx ∈ K for each λ ∈ R 0 , and for all α, β ∈ R 0 and x, y ∈ K, we have αx + βy ∈ K. The dual cone K ⊂ R n of a cone K ⊂ R n is defined as We recall a basic result from convex analysis [22, Proposition 1.1.3]: where the notation K η ⊥ x ∈ K is the short-hand for three statements: i) x ∈ K, ii) η ∈ K , and iii) x η = 0. With these basic definitions, we introduce the following class of systems for which we develop the main results of this paper. Definition 1 (Complementarity System). Given a function f : R n → R n and a cone K ⊂ R n , a complementarity system is described by the following differential equation: Several works exist in the literature which deal with existence and numerical construction of the solution to system (C-Sys). A recent reference [12] contains results in this direction, along with pointers to earlier works. Motivated by these works, it is stipulated in the remainder of this paper that the data of (C-Sys) satisfy the following assumption.
Assumption 1. Function f : R n → R n is locally Lipschitz continuous, f (0) = 0 and K ⊂ R n is a closed convex cone.

C. Solution of Complementarity Systems
Before proceeding with the problem formulation and the corresponding results, it is instructive to recall how a solution to (C-Sys) evolves with time, and the underlying optimization problem which may be solved to compute η. For a fixed s 0, if x(s) ∈ int(K), then N K (x(s)) = {0}, and we let η(s) = 0. As a result, for some ε > 0 and t ∈ [s, s + ε), we haveẋ(t) = f (x(t)) and x(t) ∈ int(K). However, if for somes 0, we have that x(s) ∈ bd(K), the boundary of K, then we essentially compute η(s) satisfying the relation 1 If the boundary constraint remains active over an interval [s,s + ε] for some ε > 0, that is, for each t ∈ [s,s + ε], x(t) ∈ bd(K), then η(t) satisfies the complementarity relation in (5). Using the notation in Appendix A, we say that η(t) ∈ LCCP(f (x(t)), I, T K (x(t))), which is equivalently described as the solution to the optimization problem stated in (36).
In what follows, it is also important to recall how we interpret the solution to (C-Sys) if x(0) = x 0 ∈ K. In such a case, we let and then propagate the solution with x + 0 , the projection of x 0 on K with respect to Euclidean norm. We can thus formally define the solution to (C-Sys) as follows: Definition 2. For a given initial condition x 0 ∈ R n and an interval [0, T ], a solution to (C-Sys) is an absolutely continuous function x : [0, T ] → R n , such that x(t) ∈ K for each t > 0, and x + 0 = proj K (x 0 ) := arg min z∈K |x 0 − z|.
Under Assumption 1, there exists a unique solution to (C-Sys) in the sense of Definition 2. We denote by x(t; x 0 ) the solution of (C-Sys), at time t 0 starting with initial condition x 0 at time 0. Assumption 1 also guarantees that the origin is an equilibrium and x(t; 0) = 0 is the unique trivial solution starting from x 0 = 0. Indeed, with K being a closed convex cone, we have 0 ∈ K. Under the condition f (0) = 0, we have η(t) = 0 andẋ(t) = 0, for all t 0.

III. PROBLEM FORMULATION
This article addresses some questions regarding the stability analysis of the trivial solution, the origin, for system (C-Sys), while assuming throughout that Assumption 1 holds. We first describe the appropriate notion of stability, and discuss some interesting properties that may arise due to the presence of constraints.

A. Stability Notions
We may now define as in [24], [25] the stability of the origin: it is stable if small perturbations of the initial condition at the origin lead to solutions remaining in the neighborhood of the origin for all forward times: Definition 3 (Stability). The origin is stable in the sense of Lyapunov if for every ε > 0 there exists δ > 0 such that The origin is locally asymptotically stable if it is stable in the sense of Lyapunov and there exists β > 0 such that The origin is globally asymptotically stable if the latter implication holds for arbitrary β > 0. The origin is globally exponentially stable if there exists c 0 > 0 and α > 0 such that x(t, x 0 ) c 0 e −αt x 0 , for every x 0 ∈ K.
Compared to the conventional definitions of stability for unconstrained dynamical systems, our domain of interest is reduced to the set K in system (C-Sys). Also, the vector field jumps instantaneously at the boundaries of the set K, which may have an impact on the stability of the system. The following examples motivate why it is not enough to analyze stability just by looking at the vector field f in (C-Sys).
Example 1 (Constraints make the system stable, even if the unconstrained system is unstable). Let f (x) = Ax with A = −1 −2 −1 −1 , and K = R 2 + . Matrix A is not Hurwitz stable since one of its eigenvalues is in the right-half complex plane. However, constrained system (C-Sys) is globally asymptotically stable, see our later Example 5 in Section VII for a proof based on a Lyapunov function.
Example 2 (Constraints make the system unstable, even if the unconstrained system is stable). Let f (x) = Ax with A = −1.5 −1 2 1 , and K = R 2 + . The matrix A is Hurwitz but the constrained system (C-Sys) is unstable because on the x 2 -axis, the vector field is pointing away from the origin.
Note that in the interior of K, system (C-Sys) follows the dynamicsẋ = f (x). The first example, however, shows that even if the constrained system is globally asymptotically stable, it is not possible to work with a Lyapunov function for the unconstrained system. In Example 1, the unconstrained system does not admit a positive definite function with negative definite time derivative over the entire state space. Consequently, one has to enlarge the search for Lyapunov functions to functions which are positive definite only on the admissible domain. The second example shows that even if one can find a Lyapunov function for the unconstrained system, it may not correspond to a Lyapunov function for the constrained system. Thus, the search of Lyapunov functions for the constrained system needs to be investigated differently from the unconstrained system.

B. Cone-Copositive Lyapunov Functions
Based on the above notions, one has to adapt the notion of Lyapunov functions when analyzing the stability of complementarity systems. It is thus of interest to introduce cone-copositive functions: Definition 4 (Cone-copositivity and copositivity). Let K ⊆ R n be a closed convex cone. A real-valued function h : R n → R is said to be cone-copositive with respect to K, if h(x) 0 for each x ∈ K. When K = R n + , we simply say that h is copositive.
Positive definite functions are obviously cone-copositive, regardless of the cone under consideration. However, in general, when the cone K is fixed, positive definite functions only form a subclass of the functions which are cone-copositive with respect to K. With this function class, the following definition of Lyapunov functions for (C-Sys) provides more flexibility: Definition 5 (Cone-copositive Lyapunov Function). System (C-Sys) has a continuously differentiable (global) conecopositive Lyapunov function V : R n → R with respect to K if 1) There exist class K ∞ functions 2 α, α such that 2) There exists a class K function α such that where η ∈ LCCP(f (x), I, T K (x)). Condition (7) is split into two parts because the complementarity variable resulting from a complementarity relation takes nonzero value only on the boundary of K.
Note that we require the inequalities in (7) to hold only for a particular selection of η. This aspect of our definition is in contrast with several existing works dealing with the existence of Lyapunov functions for differential inclusions [16], [51]. Our first major question relates to the existence of Lyapunov functions in the sense of Definition 5.
Problem 1: Does there exist a cone-copositive Lyapunov function for a stable complementarity system ? We address Problem 1 in Section IV and our first main result in Theorem 1 provides conditions which guarantees existence of a cone-copositive Lyapunov function for exponentially stable systems. Building on this result, and imposing further assumptions on the vector field f in (C-Sys), we are able to prove the existence of homogeneous Lyapunov functions, which are desired for computational reasons.

C. Computations Using Numerical Approximations
Our next target in the paper is to address the computational aspects of the Lyapunov functions for complementarity systems (C-Sys). While working with homogeneous vector fields, we restrict our search to rational functions of homogeneous polynomials. It is observed that such functions are dense within the class of homogeneous differentiable functions, see [2, Lemma 2.1]. Moreover, one can adapt the algorithms from the literature on copositive programming to compute these rational functions.
Problem 2: If there exists a homogeneous rational cone-copositive Lyapunov function for a stable complementarity system, how can we construct it?
The answer to this question essentially boils down to finding certain polynomials which satisfy some nonnegativity condition. Such questions have again received a lot of attention in real algebraic geometry, and Positivstellensätze provide guidance in writing algorithms for the search of Lyapunov functions. We explore two possible routes.
The first method corresponds to a discretization in the set K by taking points over a simplex. We evaluate the inequalities over that set of discrete points, and we solve for the desired coefficients of the polynomials. This discretization method provides an inner approximation of the cone of copositive polynomials, and it is seen that as the size of discretization step goes to zero, this inner approximation converges to the actual cone.
The second method relies on SOS decomposition of our function. Let R[x] denote the vector space of real polynomials in the variables The existence of an SOS decomposition is an algebraic certificate for nonnegativity of a polynomial. It is obvious that every SOS polynomial is nonnegative on R n . But the converse is not always true, that is, a nonnegative polynomial is not necessarily SOS. Dealing with positivity of a polynomial is hard but with SOS, it becomes easier as the problem boils down to semidefinite programming (SDP) or linear matrix inequalities (LMIs), a particular class of convex optimization problems for which efficient algorithms are available, as explained already in Section I. For detailed accounts on SOS and positive polynomials and the algebraic concepts, we refer to [32], [31].

IV. CONVERSE COPOSITIVE LYAPUNOV RESULT
In this section, we will establish an existence result for Lyapunov function, that is, if the system is exponentially stable then there exists a Lyapunov function, with certain properties, for such system. There exist several results in the literature on converse Lyapunov theorems for systems where the vector fields are discontinuous, see [19], [35] for switched systems, and [13] for complementarity systems. The results in [19], [35] use linearity of the flows, and the results in [13] are restricted to the class of complementarity systems where the right-hand side is Lipschitz continuous (and hence not discontinuous). Here, we study the converse result where the flow maps are not necessarily linear, and the complementarity relations may induce discontinuities in the vector field. In essence, we generalize the converse results on differential inclusions presented in [16], [51]. An essential difference compared to these results is that our system does not satisfy the regularity assumptions imposed in those works, and instead of strong stability, we address weak stability. Moreover, the structure of the system only allows construction over the admissible domain, which is a closed convex cone in our case. Our main result in this direction appears below. The proof of this theorem is a rather lengthy and technical affair and is carried out in the remainder of this section. Theorem 1. Under Assumption 1, if the origin is globally exponentially stable for system (C-Sys), then there exists a continuously differentiable cone-copositive Lyapunov function.
To prove Theorem 1, we start with the following lemma.
Lemma 2. If Assumption 1 holds and the origin is globally exponentially stable for system (C-Sys), then there exists a globally Lipschitz function f : R n → R n such that the systeṁ has a globally exponentially stable equilibrium and a continuously differentiable cone-copositive Lyapunov function V . Moreover, V is a Lyapunov function for (C-Sys).
The proof of Lemma 2 appears in Appendix B. Based on Lemma 2, it can be assumed for the proof of Theorem 1, without loss of generality, that f in (C-Sys) is a globally Lipschitz continuous vector field with modulus L and this assumption is assumed to hold in the remainder of this section. Note that, in Theorem 1, we assume the origin to be globally exponentially stable and our proof (appearing next) indeed uses that property. It remains to be seen if the proof can be adapted to the case where the origin is only asymptotically stable.
For the proof of Theorem 1, we construct the Lyapunov function for (C-Sys) by introducing a function V : where x(τ ; z) denotes the solution to system (C-Sys) at time τ 0 with x(0 + ) = z ∈ K. Note that V is defined for each z ∈ R n and not just for z ∈ K. When z ∈ K, the term x(τ ; proj K (z)) can be interpreted as the solution obtained by projecting the initial condition on K, and then propagating it continuously according to the system vector field. Thus, for τ > 0, we have x(τ ; proj K (z)) = x(τ ; z) for each z ∈ R n .

A. Bounds on Solutions
The following lemma demonstrates the continuity of solutions with respect to the initial conditions and plays an important role in the remainder of the proof.
Lemma 3. Let L be the Lipschitz modulus of f . If x and x are two solutions to system (C-Sys) that satisfy x(0) = z ∈ K and x(0) = z ∈ K, then it holds that, for each τ > 0, and for some C > 0, Proof. It will be assumed without loss of generality that z ∈ K and z ∈ K since proj K (z)−proj K ( z) z − z . By definition of the solution to (C-Sys) and monotonicity of the normal cone operator, it follows that, for each and similarly, for each where we have suppressed the dependence of x and x on the initial condition for brevity. Letting y = x(t) ∈ K, and y = x(t) ∈ K, we get the following by adding the last two inequalities: Because of the Lipschitz continuity assumption, The bound in (10a) now follows by integrating both sides, or invoking the so-called comparison lemma [30,Lemma 3.4]. To get the bound in (10b), we make use of Proposition 19 in the Appendix A which ensures that there exists a constant C η > 0 such that |η| C η |f (x(t))|. We therefore get In particular, d dt x(t) 2 −2L(1 + C η ) x(t) 2 , and hence, the inequality in (10b) follows by taking C = L(1 + C η ).
To show that V satisfies item 1) of Definition 5, let us first use the bound in (10b) from Lemma 3, so that for some C > 0. Also, exponential stability of the origin implies that x(τ ; z) c 0 e −ατ proj K (z) and hence there exists To show that V is locally Lipschitz continuous, we need the following two properties [17]: 3 The Dini subderivative of V at x in the direction v is defined as for every v ∈ R n , every z ∈ R n , and some locally bounded function φ : R n → R, with φ(z) > 0 for z = 0. The continuity of V follows directly from Lemma 3 as the exponential bound on the solutions of the system makes V a composition of continuous functions. These properties can again be shown using Lemma 3. Fix v ∈ R n . Consider a sequence of initial conditions z k = z + ε k v. We get Using the mean-value theorem, for each s 0, there exists ξ(s) between x k (s; z k ) and x(s; z) such that It follows from Lemma 3 that x k (s; z k ) − x(s; z) e Ls ε k v . Substituting these bounds in (12), we get Due to the exponential stability assumption, ξ(s) c e −αs z , for some c > 0, and hence we choose so that the bound (11) is seen to hold. Thus, V is locally Lipschitz continuous.

C. Infinitesimal Decrease in V
As the next step, we show that the function V decreases along the system vector field. In what follows, we will denote the right-hand side of (C-Sys) by F (z), so that The function V in (9) is differentiable almost everywhere because it is locally Lipschitz continuous. We next show that the Dini subderivative of V , along F (z) is negative definite.
Lemma 4. For the function V : R n → R in (9), and z ∈ K, Proof of Lemma 4. To prove (14), we need a bound on V (z) − V (z + tF (z)) for t 0 sufficiently small. We will get the desired bounds by rewriting the difference as and getting a bound on each of the two difference terms on the right-hand side. The first term V (z + tF (z)) − V (x(t; z)) can be analyzed from the following lemma: Lemma 5. For t > 0 sufficiently small, it holds that for each z ∈ K.
The proof of Lemma 5 will follow momentarily. Using the estimate (16), and the inequalities (12) and (13), we get for a fixed z ∈ K, and some C φ > 0. For the second term on the right-hand side of (15), it follows from the definition of V in (9), with x(0) = z, that Substituting the bounds from (16) and (17) in (15), we get and hence the Dini subderivative of V is negative definite for almost every z ∈ K.
Proof of Lemma 5. By definition, the solution x of system (C-Sys), with x(0) = z ∈ K, satisfies for almost all t 0. For h > 0 small enough, introduce the function z : [0, h] → R n given by where η z is such that η z = 0 for z ∈ int(K) and η z ∈ LCCP(f (z), I, T K (z)) for z ∈ bd(K). It is readily checked For t > 0 small enough, we have (x(t) − z(t) + t z(t)) ∈ K. Since K is a cone, we can take y Taking y = z(t) in (18), and adding it to the last inequality, we get To bound the terms on the right-hand side, we observe that Using Proposition 19 in the Appendix A, there is some constant C η such that Consequently, we get where we used Young's inequality for the product term Lt F (z) · z(t) − x(t) , and chose C 1 = (L + 0.5L 2 ) and C 2,z = max{0.5 F (z) 2 , C η F (z) }. Solving the differential inequality, and using the fact that, z(0) = x(0), we get Solving the integral on the right, we get for some C 3,z > 0, whence the estimate in (16) follows.

D. Regularization of V
The next step is to regularize V so that we obtain a continuously differentiable Lyapunov function.

Lemma 6.
Under Assumption 1, if the origin is globally exponentially stable for system (C-Sys), then there exists a continuously differentiable cone-copositive Lyapunov function.
Proof of Lemma 6. Using the function V in (9) as a template, we introduce where ψ σ , σ ∈ (0, 1), is the so-called mollifier that satisfies: ψ σ ∈ C ∞ (R n , R + ), supp(ψ σ ) ⊂ B(0, σ), and R n ψ σ (y)dy = 1. It follows from standard texts in functional analysis, see for example [6,Proposition 4.21], that V σ is continuously differentiable and for every ε > 0 and a compact set U c , there exists σ > 0, such that for every σ ∈ (0, σ), we get |V (x) − V σ (x)| < ε for each x ∈ U c . Next, we show that ∇V σ (z), F (z) approximates DV (z, F (z)), for z ∈ K. Indeed, for a given y ∈ R n , and z ∈ K, letz y = proj K (z − y). It then follows that 4 where the bound on the first integral is due to Lemma 4, and the bound on the second integral is obtained from the Lipschitz continuity of f and that of η given in Proposition 19 in Appendix A. Thus, on each compact set excluding the origin, we can find a function V σ such that ∇V σ (z), F (z) is negative definite.
Let us now consider {U i } i∈N to be a locally finite open cover of R n \ {0} with U i bounded and 0 ∈ cl(U i ), for each i ∈ N. Let {χ i } i∈N be a subordinated C 1 partition of unity. For each i ∈ N, and ε i > 0, we can choose the function V i such that |V (x) − V i (x)| < ε i , and ∇V i (x), F (x) is negative, for each x ∈ cl(U i ). Let V : R n → R + be such that V (0) = 0 and V (x) := i∈N χ i (x)V i (x) for x = 0, then following the analysis in [16, Pages 106-108], it is seen that V is a cone-copositive Lyapunov function which is C 1 on R n \ {0}, and continuous at {0}. Finally, to achieve differentiability at the origin, we can introduce a positive definite function β : R + → R + with β (s) > 0 for each s > 0 such that W (x) = β(V (x)) is a continuously differentiable cone-copositive Lyapunov function with respect to K.
Remark 7. The construction given in the proof of Lemma 6 actually gives a C ∞ (R n , R) Lyapunov function. This regularization technique is inspired by [16], and has also been used for smoothening of locally Lipschitz Lyapunov functions for hybrid systems [23,Chapter 7] and switched systems [20].

V. HOMOGENEOUS AND POLYNOMIAL LYAPUNOV APPROXIMATIONS
For numerical purposes, it is useful to show the existence of homogeneous Lyapunov functions, and if possible, polynomial Lyapunov functions. In this section, we address the question whether there exist Lyapunov functions with these additional properties.

A. Homogeneous Lyapunov Functions
First, we show that the previous developments can be generalized to construct a homogeneous Lyapunov function when the vector field f in the system description (C-Sys) is homogeneous.
for each x ∈ K and λ 0.
The next two statements are generalizations of results given in [47].

Proposition 8.
Under Assumption 1, if the origin is locally exponentially stable for system (C-Sys) with f homogeneous, then it is also globally exponentially stable.
Proof. We first show that if x : [0, ∞) → K is a solution that satisfies (C-Sys) starting with initial condition x 0 , then for each λ 0 and t 0, the function y(t) = λx(λ d−1 t) is also a solution to system (C-Sys) starting with initial condition λx 0 . It follows by inspection that y(t) ∈ K, for each t 0. Noting that for each z ∈ K, and λ > 0, there exists z ∈ K such that z = λz, we get and henceẏ(t) − f (y(t)) ∈ −N K (y(t)) for almost every t 0.
Since the origin is locally exponentially stable, there is an open set relative to K, say R 0 , such that for each x(0) ∈ R 0 , the corresponding solution x converges to the origin. For an initial condition y(0) ∈ R 0 , there is a constant λ > 0 such that y(0) = λx(0), with x(0) ∈ R 0 . Since the solutions are unique, the above reasoning shows that the solution starting from y(0) stays within a bounded set and converges to the origin.
The next result allows us to construct a homogeneous Lyapunov function under local exponential stability. The proof is inspired from [47] and is provided in Appendix C just for the sake of completeness.
Proposition 9. Consider dynamical system (C-Sys) with f homogeneous and the origin locally exponential stable. Let W ∈ C ∞ (R n , R) be a cone-copositive Lyapunov function for (C-Sys). Let a ∈ C ∞ (R, R) be such that and ∇a(s) 0, for each s ∈ R. Let k be a positive integer. Then the function is a cone-copositive Lyapunov function of class C k−1 on R n \{0}, and it satisfies for all x ∈ R n \{0} and s > 0.

B. Polynomial Approximation
For the class of numerical algorithms that we will propose in the next section, it is important to show that the cone-copositive Lyapunov functions of (C-Sys) can actually be approximated by polynomial functions. Among the existing results in this direction, it is seen that the existence of polynomial Lyapunov functions has been shown under certain restrictions only. In [42], the authors use generalizations of the Weierstrass approximation theorem for nonlinear systems with smooth vector fields to show existence of polynomial Lyapunov functions on compact sets for exponentially stable systems. In case of switched systems, the existence of polynomial Lyapunov functions has been proven in [35] when the solution maps (parameterized by time) are linear functions of the initial condition. Such methods cannot be generalized here because our vector fields are not even continuous, and even with f linear in (C-Sys), the resulting solution maps for the complementarity systems are nonlinear and hence nonconvex. As an example of this last observation, we consider the following example: Example 3 (Constraints make the solution space nonlinear). Let f (x) = Ax with A = 0 1 −1 0 and K = R 2 + and let x 1 (0) = (a, 0) and x 2 (0) = (0, b) . Let x i : R → R 2 be the solution starting with initial condition x i (0), i = 1, 2, and z denote the solution starting with initial condition x 1 (0) + x 2 (0). It can be checked that z(t) does not equal x 1 (t) + x 2 (t), for any t > 0 because we have x 1 (t) = x 1 (0) for t 0, x 2 (t) = e At x 2 (0) = (b sin(t), b cos(t)) which gives us x 1 (t) + x 2 (t) = (a + b sin(t), b cos(t)), but we have z(t) = e At (x 1 (0) + x 2 (0)) = (a cos(t) + b sin(t), −a sin(t) + b cos(t)) which is not equal to x 1 (t) + x 2 (t) for a, b = 0. These discussions and the example suggest that it may not be possible to find a homogeneous polynomial approximation to the Lyapunov function proposed in Theorem 1. Due to lack of any known results on density of homogeneous polynomials in the class of differentiable functions, we enlarge our search to rational functions whose numerator and denominator are homogeneous polynomials. For such functions, we have the following density result [2, Lemma 2.1]: Proposition 10. Let W ∈ C 1 (R n ; R + ) be a homogeneous function of degree d and > 0 be a given scalar. There exist an even integer r and a homogeneous polynomial p of degree r + d, such that where S n−1 denotes the unit sphere in R n and W (x) = W (x) − p(x) x r . With such a rational function in hand which approximates the homogeneous function from Proposition 9 (in terms of value and gradient) to desired accuracy, one can establish the existence of a rational homogeneous cone-copositive Lyapunov function.

VI. NUMERICAL CONSTRUCTION USING CONVEX OPTIMIZATION
In the previous section, we motivated the need for computing cone-copositive homogeneous Lyapunov functions for the class of constrained dynamical systems (C-Sys). Proposition 10 suggests that for a certain class of complementarity systems, we can reduce our search of Lyapunov functions to the space of rational polynomial functions, where the denominator has a certain structure. By fixing the denominator, we reformulate our problem as finding the numerator in the form of polynomial which satisfies certain inequalities. We carry out the steps by specifying the inequalities that need to be satisfied, and the algorithms using convex optimization methods that can be implemented for computing such functions.
Just as a quick motivation for what follows, we remark that contrary to unconstrained linear systems, the following example shows that copositive Lyapunov functions cannot be simply obtained by solving a linear equation, and hence there is a need to develop tools for computing them. by solving A G + GA = − H. This example shows that for a given matrix A, we can have a copositive matrix H without the existence of G copositive verifying A G + GA = −H , but with existence of some G such that −A G − GA is copositive.
We now establish the inequalities which will be used in our algorithms to find copositive homogeneous Lyapunov functions. We restrict our attention to full-dimensional polyhedral cones, that is, K = {x ∈ R n |F x 0} with nonempty interior. By using Proposition 10, let where r is a nonnegative integer, and h(·) is a homogeneous polynomial of degree at least 2r + 1, copositive on K. Here, we used the notation that x = (x 1 , x 2 , . . . , x n ) ∈ K.
As we know, finding such Lyapunov function is equivalent to finding V that satisfies the inequalities: with η = LCCP(f (x), I, T K (x)). The numerator is denoted by x, f (x) + η which is a homogeneous polynomial if h and f are homogeneous polynomials. So we have Thus, finding a copositive V for system (C-Sys) with the structure imposed in Proposition 10 boils down to finding h and s such that Since η is nonzero only on the boundaries of K, we replace the second inequality in (23) by inequalities with respect to each face of polyhedron K. Let S i := {x ∈ K | (F x) i = 0}, i ∈ {1, . . . , n K } denote the faces of K. Let for all x ∈ S i where η i = LCCP(f (x), I, T K (x)). In the interior of K, we have η = 0 so let . Consequently, the inequalities used for finding V can be rewritten as follows: While computing cone-copositive Lyapunov functions V using inequalities (24), we notice that we are faced with two problems, which prevent us from using conventional SOS techniques. The first problem is that there is no readily available Positivstellensatz for unbounded domains like cones. The second problem is that our Lyapunov functions are not necessarily SOS because a SOS polynomial is in particular positive definite but our systems require searching for a Lyapunov function beyond positive definite functions.
To overcome these problems, we study two techniques for finding polynomials that satisfy (24). In what follows, we assume that K = R n + , that is, K is the positive orthant with n faces so that S i = {x ∈ R n + | x i = 0}, i = 1, . . . , n. The more general case of polyhedral cones can be covered with state transformations or decompositions but such details are being avoided for the ease of exposition.

A. Discretization Method
The basic idea behind the discretization methods is to select a certain number of points in the cone R n + and evaluate the inequalities (24) with a certain polynomial function parameterized by finitely many unknowns. This allows us to construct an inner approximation of copositive polynomials with respect to cone R n + . In the literature, there exist algorithms for checking copositivity of a matrix using discretization methods [10], [11], [21], and using a moment relaxation hierarchy [37]. Here, we restrict ourselves to discretization schemes and generalize the existing algorithms for arbitrary nonlinear polynomials (not necessarily quadratic functions).
To describe this discretization algorithm, let us first consider the convex cone of copositive polynomials g is homogeneous and where R d [x] denotes the ring of polynomials of degree d, over the field of reals, in x ∈ R n . We will establish an inner approximation of C based on simplicial partitions inside cone R n + . To do so, we first need to introduce tensors, which generalize the notion of a matrix, and will be used for compact representation of polynomials of our interest.

Definition 7.
A tensor A of order d over R n is a multilinear form and a i1,i2,...,id corresponds to a real number from a table with n d entries, indexed by i 1 , i 2 , . . . , i d ∈ {1, . . . , n}.
A classic matrix A ∈ R n×n describes a tensor of order 2 over R n , also called a quadratic form, where the coefficients of the quadratic form belong to a table with n 2 entries a i,j with i, j = {1, . . . , n}. A general homogeneous polynomial g ∈ R d [x], with d 2, can be written as Using the tensor representation, g can also be compactly written in the form where G is a symmetric tensor. The following lemma shows an equivalent expression for copositivity which we will consider all along this section.
Lemma 11. Consider a homogenous polynomial g ∈ R d [x] of degree d and let · denote any norm in R n . We have : Take x ∈ R n + with x = 1. If x = 0 then x = 0 and g(0) = 0 because of the homogeneity of g. If x > 0 then x := x x fulfills x = 1, therefore g(x) = g( x x) = x d g( x) 0, for all x ∈ R n + which means g ∈ C. If we choose the 1-norm, then the set ∆ S := x ∈ R n + | x 1 = 1 is the standard simplex. Because of Lemma 11, copositivity of a homogenous polynomial g is then expressed as Our goal is to discretize the simplex ∆ S and obtain a hierarchy of linear inequalities with respect to the discretization points which allow us to approximate the set C.
For a given partition P m = ∆ 1 , . . . , ∆ m of ∆ S and a homogeneous polynomial g defined as in (26), let us consider the set Q k , which contains all the vertices of ∆ k , and moreover, let the set I pm,d Pm be defined as The following proposition shows that I p,d Pl l∈N is a sequence of inner approximation which approximates the cone of copositive polynomials under the condition that the diameter of the simplicial partition goes to zero. Proposition 12. Let {P l } l∈N be a sequence of simplicial partitions of ∆ S such that δ(P l ) → 0. Then, we have Proposition 12 ensures that if we construct a hierarchy of linear programs by making the partition finer, we can find a rational polynomial Lyapunov function for homogenous systems if the origin is exponentially stable. A pseudocode for implementing this method is given in the form of Algorithm 1 in Appendix D.
To prove Proposition 12, we need the following two lemmas. The first one gives us sufficient conditions for copositivity and the second one a necessary condition for strict copositivity.
This gives For x ∈ ∆, we have λ i 0, and by the assumption (28), we get g(x) 0 for all x ∈ ∆.
Proof of Proposition 12. Take g ∈ int C, which means that g is strictly copositive. Lemma 14 implies that there exists l 0 ∈ N, such that g ∈ I p,d Pl 0 . Then g ∈ l∈N I d Pl , and int C ⊆ l∈N I d Pl .
Next, for proving l∈N I d Pl ⊆ C, we have to show that I d Pl ⊆ C for some l ∈ N. So take g ∈ I d Pl for some l ∈ N. To prove g ∈ C, it is sufficient to prove nonnegativity of g(x) for x ∈ ∆ S . Let us choose an arbitrary x ∈ ∆ S , then x ∈ ∆ k for some ∆ k ∈ P l . By direct use of Lemma 13, we get g(x) = G[x, x, . . . , x] 0 for all x ∈ ∆ S .
Lastly, since C = int C, we get C = l∈N I d Pl .

B. Sum-of-Squares (SOS) Method
A commonly employed tool for checking the positivity of a polynomial is to write it in the form of a sum of squares of other polynomials. While testing positivity is a computationally hard problem, the question of finding an SOS decomposition of a polynomial is actually a semidefinite program [43]. The crux of such ideas can be found in [41] and its application to copositivity is sketched in [40].
The basic idea is to get rid of the constraint x ∈ R n + . We let x i = y 2 i , i ∈ {1, · · · , n} be the change of variable where y 2 is the short-hand for (y 2 1 , . . . , y 2 n ). Clearly we have h(x) 0, ∀x ∈ R n + ⇐⇒ h(y 2 ) 0, ∀y ∈ R n . Then, the inequalities (24) are rewritten as follows P h (y) := h(y 2 ) 0, ∀y ∈ R n (29a) where h, s o , s i are homogeneous polynomials. Next, we define the polynomials where d is an integer. It is obvious that inequalities (29)  Proof. We will carry out the proof only for h and it will be similar for the other polynomials. Let has positive coefficients .
We notice that C d h (y) has only positive coefficients then the polynomial P h (y) = h(y 2 ) is SOS and since P h (y) is multiplied by y 2d , it follows that P The copositivity of h is equivalent to the positivity of P h . And since h is homogeneous, this will be also equivalent to the positivity of P h on the unit ball which means P h (y) 0, ∀y ∈ R n , n i=1 y 2 i = 1. By substituting y 2 i by z i , we obtain P h (z) 0 , ∀z 0, n i=1 z i = 1. Let us now recall Pólya's Theorem, see [32] and [44] for the proof.  Applying this Theorem to the homogeneous polynomial P h (z), we obtain that for sufficiently large d ∈ N, all the coefficients of the polynomial P (d) is SOS in view of the fact that C d n ⊆ K d n . To sum up this section, the foregoing result allows us to write an algorithm to compute the polynomials P

VII. EXAMPLES AND SIMULATIONS
In this section, we compute copositive polynomial Lyapunov functions for complementarity systems by implementing our two methods (discretization and SOS). In our examples, the YALMIP toolbox in Matlab is used to input the LP and SOS optimization problems and solve them with the conic solver MOSEK.
Example 5 (Quadratic Lyapunov function by the discretization method). Consider system (C-Sys) with f (x) = Ax and A = −1 −2 −1 −1 and K = R 2 + . We apply the discretization method on the standard simplex ∆ S := x ∈ R 2 + | x 1 = 1 of Algorithm 1. Starting with a degree 2 polynomial, we solve for its coefficients at the vertices of ∆ S . This procedure in Algorithm 1 is applied by partitioning all the simplices at each step by a factor of half and solving certain inequalities at the vertices of resulting simplices. For this example, we found V (x) = x 2 1 + x 1 x 2 + x 2 2 in four iterations, that is, Algorithm 1 terminates with δ = 1/16. Example 6 (Cubic Lyapunov function by the discretization method). Consider system (C-Sys) with K = R 2 + and Applying the discretization method of Algorithm 1, after 4 iterations with δ = 1 16 , we obtain Example 7 (Quadratic Lyapunov function by SOS method). Consider system (C-Sys) with f (x) = Ax and A = −1 10 0 −2 and K = R 2 + . Following Algorithm 2, we express positivity condition on desired polynomials by requiring them to be SOS, and use the YALMIP command solvesos which calls a semidefinite solver to yield the desired coefficients. We obtain V (x) = 0.1x 2 1 + 0.1916x 1 x 2 + 1. V (x) = 2.3234x 2 1 + 3.6729x 1 x 2 + 1.7352x 2 2 + 1.1273x 1 x 3 + 2.6769x 2 x 3 + 1.2820x 2 3 .
This polynomial is not positive definite since one of the eigenvalues of the corresponding matrix is negative. The unit level set of this polynomial Lyapunov function is shown in Figure 2.
These examples just provide an illustration of two classes of algorithms primarily used for checking positivity or copositivity of polynomials, and how they can be used for computing Lyapunov functions with constrained dynamics. The survey article [9] provides an overview of these methods, along with some other techniques, which appear in general in the literature on checking copositivity. Further questions such as using other algorithms or comparing computational complexity of different methods require further investigation.

VIII. CONCLUSIONS
This article addressed the stability analysis for a class of complementarity systems using the method of Lyapunov functions. Questions pertaining to the existence of cone-copositive Lyapunov functions were answered in the affirmative for exponentially stable systems. Some refinements of this result, under certain conditions on the vector field in the system dynamics, allow us to restrict our search for cone-copositive Lyapunov functions within the class of homogeneous and rational polynomials. These statements indeed bring tractability to the numerical methods that have been proposed in this paper for computing Lyapunov functions. In particular, two hierarchies of convex optimization problems are obtained using the methods based on discretization and SOS approximation, respectively.
Several immediate questions of interest emerge from our work which require further investigation. The first one among those is to extend our results to broader classes of complementarity systems. Systems of the form (C-Sys) are one particular class of relative degree one systems, but in applications, one sees more complex complementarity systems of the form studied in [50]. In such a wider class of systems, one sees different kinds of constraints on the state trajectories. Moreover, the constraints may vary with time in which case one has to consider the possibility of time-varying Lyapunov functions. It would be interesting to consider converse questions for this broader class of systems.
Some extensions at the level of designing algorithms are also of potential interest. At this moment our algorithms are specifically adapted to the constraint sets of the form positive orthant or an invertible linear transformation of such sets. Adapting this technique to more generic sets remains to be seen. Also, in our current treatment, we have considered discretization algorithms in R 2 and R 3 , where it is relatively straightforward to write algorithms for partition of simplices. It remains to be seen how the algorithms for simplical partition in higher dimensions perform in computing Lyapunov functions.