  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.
