# Stability Analysis of Discrete-Time Linear Complementarity Systems

Arvind U. Raghunathan Mitsubishi Electric Research Laboratories
<EMAIL_ADDRESS>Jeffrey T. Linderoth Department of Industrial and
Systems Engineering, University of Wisconsin-Madison (linderoth@wisc.edu).

###### Abstract

A Discrete-Time Linear Complementarity System (DLCS) is a dynamical system in
discrete time whose state evolution is governed by linear dynamics in states
and algebraic variables that solve a Linear Complementarity Problem (LCP). The
DLCS is the hybrid dynamical system that is the discrete-time counterpart of
the well-known Linear Complementarity System (LCS). We derive sufficient
conditions for Lyapunov stability of a DLCS when using a quadratic Lyapunov
function that depends only on the state variables and a quadratic Lyapunov
function that depends both on the state and the algebraic variables. The
sufficient conditions require checking the feasibility of a copositive program
over nonconvex cones. Our results only assume that the LCP is solvable and do
not require the solutions to be unique. We devise a novel, exact cutting plane
algorithm for the verification of stability and the computation of the
Lyapunov functions. To the best of our knowledge, our algorithm is the first
exact approach for stability verification of DLCS. A number of numerical
examples are presented to illustrate the approach. Though our main object of
study in this paper is the DLCS, the proposed algorithm can be readily applied
to the stability verification of LCS. In this context, we show the equivalence
between the stability of a LCS and the DLCS resulting from a time-stepping
procedure applied to the LCS for all sufficiently small time steps.

##  1 Introduction

A Discrete-Time Linear Complementarity System (DLCS) is linear dynamical
system in a discrete time where the state evolution is governed by
_complementarity_ constraints. Mathematically, a DLCS is written as

$\displaystyle x_{k+1}=$ $\displaystyle\;Ax_{k}+C\lambda_{k+1}$ (1a)
$\displaystyle 0\leq\lambda_{k+1}\perp$
$\displaystyle\;Dx_{k}+F\lambda_{k+1}\geq 0,$ (1b)

where $x_{k}\in\mathbb{R}^{n_{x}}$ is the state of the system at the time-step
$k$ and _complementarity variables_ $\lambda_{k}\in\mathbb{R}^{n_{c}}$ are the
algebraic variables satisfying complementarity constraints in (1b). The
dimensions of the matrices in (1) are $A\in\mathbb{R}^{n_{x}\times n_{x}}$,
$C\in\mathbb{R}^{n_{x}\times n_{c}}$, $D\in\mathbb{R}^{n_{c}\times n_{x}}$,
and $F\in\mathbb{R}^{n_{c}\times n_{c}}$. The complementarity constraints (1b)
distinguish the DLCS from a standard linear dynamical system. A DLCS can be
naturally derived by time sampling a Linear Complementarity System (LCS) [18,
19, 8]. A LCS is a continuous-time linear dynamical systems governed by
_complementarity_ constraints and can be written mathematically as follows:

$\displaystyle\frac{dx}{dt}=$
$\displaystyle\;\widetilde{A}x(t)+\widetilde{C}\lambda(t)$ (2a) $\displaystyle
0\leq\lambda(t)\perp$
$\displaystyle\;\widetilde{D}x(t)+\widetilde{F}\lambda(t)\geq 0,$ (2b)

where the matrices $\widetilde{A},\widetilde{C},\widetilde{D},\widetilde{F}$
have conformal dimensions. An appropriate time-stepping scheme such as
_explicit_ or _implicit_ Euler [31] can be shown to yield the DLCS (1). For
example, Pang and Stewart [31, §3] propose simulating a LCS using time-
stepping formulations of the form

$\displaystyle x_{k+1}=$ $\displaystyle\;x_{k}+\Delta
t(\widetilde{A}x_{k+\theta}+\widetilde{C}\lambda_{k+1})$ (3a) $\displaystyle
0\leq\lambda_{k+1}\perp$
$\displaystyle\;\widetilde{D}x_{k+1}+\widetilde{F}\lambda_{k+1}\geq 0,$ (3b)

where $\Delta t$ is the time-step, $\theta\in[0,1]$ and
$x_{k+\theta}=(1-\theta)x_{k}+\theta x_{k+1}$. The dynamics of the time-
stepping formulations can be transformed to a DLCS by rearranging (3a) as

$\displaystyle x_{k+1}=$ $\displaystyle\;(I_{n_{x}}+\Delta
t\widetilde{A}_{\theta}^{-1}\widetilde{A})x_{k}+\Delta
t\widetilde{A}_{\theta}^{-1}\widetilde{C}\lambda_{k+1}$ (4a) $\displaystyle
0\leq\lambda_{k+1}\perp$
$\displaystyle\;\widetilde{D}x_{k+1}+\widetilde{F}\lambda_{k+1}\geq 0$ (4b)

where $\widetilde{A}_{\theta}=(I_{n_{x}}-\theta\Delta t\widetilde{A})^{-1}$.
Specifically, substituting for $x_{k+1}$ in (4b) using (4a), we obtain the
DLCS (1). The use of $\theta=0,\frac{1}{2},1$ results in the explicit Euler,
trapezoidal, or implicit Euler time-stepping schemes, respectively.
Convergence of the solutions of the time-stepping formulations (3) as $\Delta
t\rightarrow 0$ to a solution of the LCS (2) are provided in [31].

In dynamical systems, the discrete-time sampled variant plays an important
role in implementing control algorithms in practice. For example, discrete-
time control algorithms, such as model predictive control (MPC) [26], rely on
the time-stepping formulations of continuous-time linear dynamical systems. In
an analogous manner, the DLCS naturally arises in applications that can be
modeled with LCS. The modeling and study of dynamical systems with
complementarity conditions has been steadily increasing, with application in
electrical circuit simulation, robotics, nonsmooth mechanics, economics,
bioengineering [1, 5, 9, 10, 29, 30, 35, 36, 39, 40, 41]. Another important
recent application of interest is robotics, where complementarity is required
to model the friction forces that arise in contact-based manipulation tasks
[12, 25, 27, 32, 34, 33].

Stability is a fundamental issue in dynamical systems, and stability
verification by computing stability certificates is a key to guaranteeing
performance and safety in real-world systems. Lyapunov stability is a widely
used concept for analyzing the stability of dynamical systems [23], and its
extension to hybrid and switched systems has been considered by several
authors [4, 15, 20, 21, 22, 24, 38]. The papers extending stability to hybrid
systems assume that a Lyapunov-like function exists for each mode’s vector
field and for the entire state space. In many hybrid and switched systems,
however, each mode is active only over a subset of the state space, especially
for those systems whose switchings are triggered by state evolution, such as
in the LCS and DLCS. Hence, the above results are rather restrictive, even for
linear switched systems. To rectify this, Çamlibel and Schumacher [11]
proposed copositive Lyapunov functions for conewise linear systems in which
the feasible region of each mode is a polyhedral cone. Bundfuss and Dür [6]
present necessary and sufficient conditions for the existence of such Lyapunov
functions when the cone is polyhedral. Çamlibel, Pang and Shen [8] derived
sufficient conditions for the stability of LCS (2). Extending standard
approaches that employ piecewise-quadratic Lyapunov functions depending only
on state $x(t)$, the authors in [8] proposed an extended quadratic Lyapunov
function that depends on both $x(t)$ and $\lambda(t)$. However, the authors
did not provide an algorithm for computing such a Lyapunov function. Recently,
Aydinoglu, Preciado, and Posa [2] extended the sufficient conditions of [8]
under the assumption of an existence of a feedback controller that depends on
both $x(t)$ and $\lambda(t)$. The authors also proposed an algorithm based on
bilinear matrix inequalities for jointly determining the feedback controller
and the Lyapunov function. However, we are not aware of similar results for
the DLCS.

In this paper, we focus on the stability analysis of the DLCS and the
development of algorithms for computing a Lyapunov function certifying
stability. The paper makes three main contributions.

First, in §3 we derive sufficient conditions for the stability of a DLCS using
two quadratic Lyapunov functions. The first Lyapunov function is a quadratic
function that depends only on the states of the DLCS (1). Such a Lyapunov
function is referred to as a _Common Quadratic Lyapunov Function_ (CQLF) [24,
38]. The second Lyapunov function is a quadratic function that depends on both
the states and complementarity variables of the DLCS (1). This is the
discrete-time analogue of the _Extended Quadratic Lyapunov Function_ (EQLF)
introduced in [8]. The stability conditions are provided under the assumption
that $F$ is a Q-matrix and a R0-matrix [14]. The Q-matrix is the largest class
of matrices for which the complementarity constraints in (1b) have a solution
for all $x_{k}$. Our assumptions on the matrix $F$ are weaker than those used
in [2, 8].

Second, in §4 we present exact algorithms for computing the CQLF and EQLF for
a DLCS. The computation of Lyapunov functions requires finding a feasible
solution for a copositive program that is defined over a nonconvex cone—in
particular, the union of convex cones. It is well-known that checking whether
a given matrix is copositive is a co-NP-complete problem [28]. We pose the
problem of computing the Lyapunov functions as semi-infinite programming
problem and devise a cutting-plane algorithm. The cutting plane algorithm
alternates a linear programming master problem and a separation problem
formulated as a Mixed Integer Quadratically Constrained Program (MIQCP).
Numerical experiments demonstrating the utility of the algorithm are presented
in §5.

Third, in §6 we consider the case where the DLCS is obtained by a time-
stepping scheme applied to a LCS. An important question as yet unanswered in
the literature is whether stability of the LCS is equivalent to stability of
the DLCS when the time-step $\Delta t$ is sufficiently small. We show
equivalence for both the case of CQLF and EQLF under the assumption that the
matrix $\widetilde{F}$ in the LCS (2) is a P-matrix.

_Notation._ The set of reals, nonnegative reals, integers and nonnegative
integers is denoted by $\mathbb{R}$, $\mathbb{R}_{+}$, $\mathbb{Z}$ and
$\mathbb{Z}_{+}$ respectively. The 1-norm ball is denoted as
$\mathbb{B}^{n}_{1}=\\{w\in\mathbb{R}^{n}\,|\,\|w\|_{1}=1\\}$. For a vector
$v\in\mathbb{R}^{n}$ ${\left[v\right]}_{i}$ denotes the $i$-th component of
the vector. For $v\in\mathbb{R}^{n}$ the notation ${\left[v\right]}_{\alpha}$
for $\alpha\subset\\{1,\ldots,n\\}$ denotes the subvector obtained by removing
the elements ${\left[v\right]}_{i}$ for $i\notin\alpha$. The notation $(x,y)$
of two column vector $x,y$ is also used to represent a vertical stacking of
the vectors $\begin{pmatrix}x\\\ y\end{pmatrix}$. For a matrix
$M\in\mathbb{R}^{n\times m}$ ${\left[M\right]}_{ij}$ denotes the $(i,j)$-th
entry in the matrix. For a matrix $M$ and index sets
$\alpha\subseteq\\{1,\ldots,n\\}$, $\beta\subseteq\\{1,\ldots,m\\}$ the
notation ${\left[M\right]}_{\alpha\beta}$ denotes the submatrix of $M$ formed
by removing rows not in $\alpha$ and removing columns not in $\beta$. The
notation ${\left[M\right]}_{\alpha\bullet}$ and
${\left[M\right]}_{\bullet\alpha}$ refers to the submatrix of $M$ obtained by
retaining only the rows and columns in $\alpha$ respectively. The notation
$I_{n}$ represents the $n\times n$ identity matrix. The set $\mathbb{S}^{n}$
denotes the $n\times n$ real, symmetric matrices and
$\mathbb{S}^{n}_{+},\mathbb{S}^{n}_{++}$ denote the set of symmetric positive
semidefinite and positive definite matrices. A symmetric positive
(semi)definite matrix $N\in\mathbb{S}^{n}$ is denoted as
$N\succ(\succcurlyeq)0$.

##  2 Background

In this section, we present relevant results and notation that are used
throughout the paper.

###  2.1 Linear Complementarity Problem (LCP)

The Linear Complementarity Problem (LCP), denoted as LCP($q$,$M$), is to find
a solution $\lambda\in\mathbb{R}^{n_{c}}$ of

$0\leq\lambda\perp M\lambda+q\geq 0,$ (5)

where $q\in\mathbb{R}^{n_{c}},M\in\mathbb{R}^{n_{c}\times n_{c}}$. The
solution set of the LCP (5) is denoted as SOL($q$,$M$)
$=\\{\lambda\,|\,\lambda\text{ satisfies }\eqref{lcp}\\}$. If $q=Nx$ for some
matrix $N$, the graph of $\text{SOL}(Nx,M)$ is denoted as
$\text{Gr}\,\text{SOL}(Nx,M)$
$=\\{(x,\lambda)\,|\,\lambda\in\text{SOL}(Nx,M)\\}$. The monograph [14]
provide an in-depth treatment describing conditions for the existence and
uniqueness of solutions to (LCP). We provide a brief summary of definitions
and results that will be useful for subsequent developments.

Given a solution $\lambda\in$ SOL($q$,$M$), we define index sets
$\alpha(\lambda)$, $\beta(\lambda)$, $\gamma(\lambda)$ partitioning
$\\{1,\ldots,n_{c}\\}$ as

$\displaystyle\alpha(\lambda)$
$\displaystyle:=\;\\{i\,|\,{\left[\lambda\right]}_{i}>0={\left[M\lambda+q\right]}_{i}\\}$
(6a) $\displaystyle\beta(\lambda)$
$\displaystyle:=\;\\{i\,|\,{\left[\lambda\right]}_{i}=0={\left[M\lambda+q\right]}_{i}\\}$
(6b) $\displaystyle\gamma(\lambda)$
$\displaystyle:=\;\\{i\,|\,{\left[\lambda\right]}_{i}=0<{\left[M\lambda+q\right]}_{i}\\}.$
(6c)

If $\beta(\lambda)=\emptyset$, then $\lambda$ is called a _strict
complementarity solution_ of LCP($q$,$M$). Otherwise, it is called a
_nonstrict complementarity solution_. We will suppress $\lambda$ when the
dependence is clear from the context.

A matrix $M$ is a Q-matrix if LCP($q$,$M$) has a solution for all
$q\in\mathbb{R}^{n_{c}}$. A matrix $M$ is an R0-matrix if SOL($0$,$M$)
$=\\{0\\}$. A matrix $M$ is said to be a P-matrix if the principal minors are
all positive, i.e. $\det({\left[M\right]}_{JJ})>0$ for all
$J\subseteq\\{1,\ldots,n_{c}\\}$. It is known that a P-matrix is both an
R0-matrix and a Q-matrix [14, Theorem 3.9.22]. The following lemma stating
properties of solutions to LCP($q$,$M$) will be useful in the remainder of the
paper.

######  Lemma 1 ([14]).

Consider the LCP($q$,$M$) where the matrix $M$ is a P-matrix. Then

  1. (a)

LCP($q$,$M$) has an unique solution $\lambda(q)$ for every $q$.

  2. (b)

$\lambda(q)$ is a piecewise linear function of $q$, is globally Lipschitz
continuous, and is directionally differentiable.

  3. (c)

There exists a constant $c_{M}>0$ such that $\|\lambda(q)\|\leq c_{M}\|q\|$
for all $q\in\mathbb{R}^{n_{c}}$.

###  2.2 Lyapunov Stability

Consider the discrete-time nonlinear dynamical system

$x_{k+1}=f(x_{k});\;\;x_{0}=\hat{x},$ (7)

where $f:\mathbb{R}^{n_{x}}\rightarrow 2^{\mathbb{R}^{n_{x}}}$ is a mapping
into subsets of $\mathbb{R}^{n_{x}}$. Let $x^{e}$ be an equilibrium point of
the system (7); i.e. $f(x^{e})=0$, and let $x_{k}(\hat{x})$,
$k\in\mathbb{Z}_{+}$, be a trajectory of the system initialized at $\hat{x}$.
We do not assume the trajectory is unique.

######  Definition 1.

The equilibrium point $x^{e}$ is

  1. (i)

stable (in the sense of Lyapunov) if for each $\epsilon>0$, there is
$\delta_{\epsilon}>0$ such that

$\|\hat{x}-x^{e}\|<\delta_{\epsilon}\implies\|x_{k}(\hat{x})-x^{e}\|<\epsilon\quad\forall
k\in\mathbb{Z}_{+}$

for all trajectories $x_{k}(\hat{x})$;

  2. (ii)

asymptotically stable if it is stable and there exists $\delta>0$ such that

$\|\hat{x}-x^{e}\|<\delta\implies\lim_{k\rightarrow\infty}x_{k}(\hat{x})=x^{e}\quad\forall
k\in\mathbb{Z}_{+}$

for all trajectories $x_{k}(\hat{x})$;

  3. (iii)

exponentially stable if there exist scalars $\delta>0$, $c>0$ and $\mu>0$ such
that

$\|\hat{x}-x^{e}\|<\delta\implies\|x_{k}(\hat{x})-x^{e}\|\leq
c\mu^{k}\|\hat{x}-x^{e}\|\quad\forall k\in\mathbb{Z}_{+}$

for all trajectories $x_{k}(\hat{x})$.

If the function $f(x)$ is positively homogeneous in $x$, i.e. $f(\tau x)=\tau
f(x)$ then the stability results hold _globally_ for all initial conditions
$\hat{x}\in\mathbb{R}^{n_{x}}$. Note that the functions defining the dynamics
in both LCS and DLCS are positively homogeneous. We discuss how stability
analysis for the homogeneous case can be extended to some inhomogeneous
systems in § 3.3.

###  2.3 Semidefiniteness and Matrix Copositivity

Semidefiniteness with respect to a closed, not necessarily convex, cone ${\cal
K}$ is called copositivity with respect to ${\cal K}$. Specifically, a matrix
$M\in\mathbb{R}^{n\times n}$ is said to be _copositive_ with respect to a cone
${{\cal K}}\subseteq\mathbb{R}^{n}$ if $u^{T}Mu\geq 0$ for all $u\in{{\cal
K}}$. In this case, we say $M$ is a ${\cal K}$-copositive matrix, and denote
this fact by $M\succcurlyeq_{{\cal K}}0$. Similarly, if $u^{T}Mu>0$ for all
$u\in{\cal K}$ then $M$ is a strict ${\cal K}$-copositive matrix, which we
denote by $M\succ_{{\cal K}}0$. Note that the standard semidefinite and
copositivity condtions are obtained by choosing ${\cal K}=\mathbb{R}^{n}$ and
${\cal K}=\mathbb{R}^{n}_{+}$ respectively. For additional details, the reader
is referred to [7, 16].

##  3 Stability of DLCS

In this section we present stability conditions for the DLCS at $x^{e}=0$. We
assume that $(x^{e},\lambda^{e})=(0,0)$ is the unique equilibrium point of
DLCS, i.e. $(x^{e},\lambda^{e})$ is the only solution of

$\displaystyle 0=$ $\displaystyle\;Ax+C\lambda$ (8a) $\displaystyle
0\leq\lambda\perp$ $\displaystyle\;Dx+F\lambda\geq 0.$ (8b)

The uniqueness of the equilibrium is guaranteed if $F$ is a R0-matrix i.e.
LCP($0$,$F$) $=\\{0\\}$. We first show sufficient conditions for a CQLF in
§3.1 and extend these to sufficient conditions for an EQLF in §3.2. A
discussion on the applicability of stability results to inhomogeneous DLCS is
presented in §3.3 under the assumption that $F$ is a P-matrix.

###  3.1 Common Quadratic Lyapunov Function (CQLF)

We will first present sufficient conditions for stability of the equilibrium
point using a quadratic Lyapunov function that depends only on the states, or
a Common Quadratic Lyapunov function (CQLF) [24, 38].

######  Theorem 1.

Consider the DLCS in (1). Assume that $F$ is a Q-matrix and R0-matrix. Let
$\psi:\mathbb{R}^{n_{x}+n_{c}}\rightarrow\mathbb{R}$ denote a quadratic
function

$\displaystyle\psi(x,\lambda)$ $\displaystyle=\begin{pmatrix}x\\\
\lambda\end{pmatrix}^{T}M(P_{xx})\begin{pmatrix}x\\\
\lambda\end{pmatrix},\mbox{ where}$ (9) $\displaystyle M(P_{xx})$
$\displaystyle=\begin{pmatrix}A^{T}P_{xx}A-P_{xx}&A^{T}P_{xx}C\\\
C^{T}P_{xx}A&C^{T}P_{xx}C\end{pmatrix},$ (10)

and $P_{xx}\in\mathbb{S}^{n_{x}}$. If there exists $P_{xx}\succ 0$, then
$x^{e}=0$ is

  1. (i)

stable if $\psi(x,\lambda)\leq 0$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$,

  2. (ii)

asymptotically stable if $\psi(x,\lambda)\leq 0$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$, $\psi(x,\lambda)=0\implies x=0$,
and

  3. (iii)

exponentially stable if $\psi(x,\lambda)<0$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$.

###### Proof.

Let $x_{0}$ be an arbitrary initial condition for the DLCS and let
$(x_{k},\lambda_{k})$ for $k\geq 1$ be a trajectory from such an initial state
with $\lambda_{k+1}\in\text{Gr}\,\text{SOL}(Dx_{k},F)$. Consider the function
$V(x_{k})=x_{k}^{T}P_{xx}x_{k}$. By the positive definiteness of $P_{xx}$,
$V(x)\geq 0$ for all $x\in\mathbb{R}^{n_{x}}$ and $V(x)=0$ if and only if
$x=x^{e}$. Thus, $V(x)$ is a Lyapunov function. The Lyapunov function at
$x_{k+1}$ is

$V(x_{k+1})=\begin{pmatrix}x_{k}\\\
\lambda_{k+1}\end{pmatrix}^{T}\begin{pmatrix}A^{T}P_{xx}A&A^{T}P_{xx}C\\\
C^{T}P_{xx}A&C^{T}P_{xx}C\end{pmatrix}\begin{pmatrix}x_{k}\\\
\lambda_{k+1}\end{pmatrix}.$

It can be verified that $\psi(x_{k},\lambda_{k+1})=V(x_{k+1})-V(x_{k})$. If
$\psi(x,\lambda)\leq 0$ for all $(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$
then $V(x_{k+1})\leq V(x_{k})$ for all $x_{k}\in\mathbb{R}^{n_{x}}$. In other
words, along all trajectories of DLCS from $x_{0}$

$\displaystyle V(x_{k+1})\leq V(x_{k})\implies V(x_{k})\leq V(x_{0}).$ (11)

Since $P_{xx}\succ 0$ there exists $\omega>0$ such that
$V(x_{k})\geq\mu\|x_{k}\|^{2}$. Further, $V(x)\leq\|P_{xx}\|\|x\|^{2}$.
Combining the two observations in (11) yields
$\|x_{k}\|\leq(\|P_{xx}\|/\mu)\|x_{0}\|^{2}$ for all $k$ and all trajectories
from $x_{0}$. The claim on stability in (i) follows. If the condition in (ii)
holds then the decrease $\psi(x,\lambda)=0$ only if $x=0$, so the the sequence
$\\{V(x_{k})\\}$ is strictly decreasing unless $x=0$. Since $V(x)$ is lower
bounded by zero and $V(x)\geq\mu\|x\|^{2}$ we have that the claim in (ii)
holds. If $\psi(x,\lambda)<0$ then there exists a $\gamma>0$ such that
$\psi(x,\lambda)\leq-\gamma\|(x,\lambda)\|^{2}$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$. Then along all trajectories of
DLCS

$V(x_{k+1})-V(x_{k})\leq-\gamma\|(x_{k},\lambda_{k+1})\|^{2}\leq-\gamma\|x_{k}\|^{2}\leq-(\gamma/\|P_{xx}\|)V(x_{k})$

where the second inequality follows from
$\|x_{k}\|^{2}\leq\|(x_{k},\lambda_{k+1})\|^{2}$ and the third from
$V(x_{k})\leq\|P_{xx}\|\|x_{k}\|^{2}$. Hence,
$V(x_{k})\leq(1-\gamma/\|P_{xx}\|)^{k}V(x_{0})$ along all trajectories of DLCS
from $x_{0}$. Combining this with $V(x)\geq\omega\|x\|^{2}$ yields
$\|x_{k}\|^{2}\leq(\|P_{xx}\|/\omega)(1-\gamma/\|P_{xx}\|)^{k}\|x_{0}\|^{2}$
and the claim in (iii) holds. ∎

Theorem 1 shows that the sufficient conditions for exponential stability using
a CQLF are to

Find $\displaystyle\;\;P_{xx}\in\mathbb{S}^{n_{x}}$ (12a) such that
$\displaystyle\;P_{xx}\succ 0,$ (12b) $\displaystyle\;-M(P_{xx})\succ_{{\cal
K}}0,$ (12c)

where $M(P_{xx})$ is defined in (10), and the set ${\cal
K}\subseteq\mathbb{R}^{n_{x}+n_{c}}$ is defined as

${{\cal K}}=\text{Gr}\,\text{SOL}(Dx,F)$ (13)

###  3.2 Extended Quadratic Lyapunov Function (EQLF)

The use of a quadratic function in the states alone can limit the set of
systems that can be stabilized. Specifically, there exist linear switching
systems that can be stabilized, but for which the CQLF does not exist [24,
§2.1.5]. This has motivated the use of unique quadratic Lyapunov functions for
each mode of the switching system with continuity imposed across the switching
boundaries [4, 15, 20, 22, 24, 38]. As noted in the introduction, this
construction can also be restrictive since the Lyapunov decrease condition for
a mode is enforced over the entire state-space rather than over the subset of
the state-space corresponding to the mode. In a significant departure from
prior literature, Çamlibel, Pang and Shen [8] proposed a quadratic Lyapunov
function that depends on $(x(t),\lambda(t))$ for the LCS (2). This Lyapunov
function is piecewise quadratic since $\lambda(t)$ is a piecewise linear
function (see Lemma 1). We follow a similar path for the DLCS and derive
sufficient conditions for stability of $(x^{e},\lambda^{e})=0$.

######  Theorem 2.

Consider the DLCS in (1). Assume that $F$ is a Q-matrix and an R0-matrix. Let
$\widehat{\psi}:\mathbb{R}^{n_{x}+2n_{c}}\rightarrow\mathbb{R}$ denote a
quadratic function

$\displaystyle\widehat{\psi}(x,\lambda,\widehat{\lambda})$
$\displaystyle=\;\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix}^{T}\widehat{M}(P)\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix},$ (14) $\displaystyle\text{with
}\widehat{M}(P)$
$\displaystyle=\;\begin{pmatrix}A^{T}P_{xx}A-P_{xx}&A^{T}P_{xx}C-P_{x\lambda}&A^{T}P_{x\lambda}\\\
C^{T}P_{xx}A-P_{x\lambda}^{T}&C^{T}P_{xx}C-P_{\lambda\lambda}&C^{T}P_{x\lambda}\\\
P_{x\lambda}^{T}A&P_{x\lambda}^{T}C&P_{\lambda\lambda}\end{pmatrix},$ (15)

where $P=\begin{pmatrix}P_{xx}&P_{x\lambda}\\\
P_{x\lambda}^{T}&P_{\lambda\lambda}\end{pmatrix}$,
$P_{xx}\in\mathbb{S}^{n_{x}}$, $P_{\lambda\lambda}\in\mathbb{S}^{n_{c}}$ and
$P_{x\lambda}\in\mathbb{R}^{n_{x}\times n_{c}}$. If there exists
$P_{xx},P_{x\lambda},P_{\lambda\lambda}$ such that the function
$\widehat{V}(x,\lambda)$

$\widehat{V}(x,\lambda)=\begin{pmatrix}x\\\
\lambda\end{pmatrix}^{T}P\begin{pmatrix}x\\\ \lambda\end{pmatrix}$ (16)

satisfies $\widehat{V}(x,\lambda)\geq 0$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$ and $\widehat{V}(x,\lambda)=0$ if
and only if $(x,\lambda)=(x^{e},\lambda^{e})$, then $(x^{e},\lambda^{e})=0$ is

  1. (i)

stable if $\widehat{\psi}(x,\lambda,\widehat{\lambda})\leq 0$ for all
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$ and
$\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)$,

  2. (ii)

asymptotically stable if $\widehat{\psi}(x,\lambda,\widehat{\lambda})\leq 0$
for all $(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$,
$\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)$,
$\widehat{\psi}(x,\lambda,\widehat{\lambda})=0\implies x=0$

  3. (iii)

exponentially stable if $\widehat{\psi}(x,\lambda,\widehat{\lambda})<0$ for
all $(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$ and
$\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)$.

###### Proof.

Let $x_{0}$ be an arbitrary initial condition for the DLCS and let
$(x_{k},\lambda_{k})$ for $k\geq 1$ be a trajectory from such an initial state
with $\lambda_{k+1}\in\text{Gr}\,\text{SOL}(Dx_{k},F)$. Then

$\displaystyle\widehat{V}(x_{k+1},\lambda_{k+2})=$
$\displaystyle\;\begin{pmatrix}x_{k+1}\\\
\lambda_{k+2}\end{pmatrix}^{T}\begin{pmatrix}P_{xx}&P_{x\lambda}\\\
P_{x\lambda}^{T}&P_{\lambda\lambda}\end{pmatrix}\begin{pmatrix}x_{k+1}\\\
\lambda_{k+2}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\;\begin{pmatrix}Ax_{k}+C\lambda_{k+1}\\\
\lambda_{k+2}\end{pmatrix}^{T}\begin{pmatrix}P_{xx}&P_{x\lambda}\\\
P_{x\lambda}^{T}&P_{\lambda\lambda}\end{pmatrix}\begin{pmatrix}Ax_{k}+C\lambda_{k+1}\\\
\lambda_{k+2}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\;\begin{pmatrix}x_{k}\\\ \lambda_{k+1}\\\
\lambda_{k+2}\end{pmatrix}^{T}\begin{pmatrix}A^{T}P_{xx}A&A^{T}P_{xx}C&A^{T}P_{x\lambda}\\\
C^{T}P_{xx}A&C^{T}P_{xx}C&C^{T}P_{x\lambda}\\\
P_{x\lambda}^{T}A&P_{x\lambda}^{T}C&P_{\lambda\lambda}\end{pmatrix}\begin{pmatrix}x_{k}\\\
\lambda_{k+1}\\\ \lambda_{k+2}\end{pmatrix}.$

It can be verified that
$\widehat{V}(x_{k+1},\lambda_{k+2})-\widehat{V}(x_{k},\lambda_{k+1})=\widehat{\psi}(x_{k},\lambda_{k+1},\lambda_{k+2})$.
The assumptions on $\widehat{V}$ ensure that
$\widehat{V}(x_{k},\lambda_{k+1})\geq 0\mbox{ and
}\widehat{V}(x_{k+1},\lambda_{k+2})\geq 0$. If the condition in (i) holds then

$\displaystyle\widehat{V}(x_{k+1},\lambda_{k+2})\leq\widehat{V}(x_{k},\lambda_{k+1})\implies\widehat{V}(x_{k},\lambda_{k+1})\leq\widehat{V}(x_{0},\lambda_{1})$
(17)

for all $k$. Further, since $\widehat{V}(x_{k},\lambda_{k+1})>0$ for all
$(x_{k},\lambda_{k+1})\in\text{Gr}\text{SOL}(Dx_{k},F)$ we have that
$\widehat{V}(x_{k},\lambda_{k+1})\geq\omega\|(x_{k},\lambda_{k+1})\|^{2}$ for
some $\omega>0$. Further,
$V(x_{0},\lambda_{1})\leq\|P\|\|(x_{0},\lambda_{1})\|^{2}$. Applying these two
in (17) yields
$\|(x_{k},\lambda_{k+1})\|^{2}\leq(\|P\|/\omega)\|(x_{0},\lambda_{1})\|^{2}$
and the claim on stability in (i) follows. If condition in (ii) holds then
$\\{V(x_{k},\lambda_{k+1})$ is strictly decreasing along the trajectories
unless $x_{k}=0$ and the claim in (ii) follows. If condition in (iii) holds
then there exists $\gamma>0$ such that
$\widehat{\psi}(x_{k},\lambda_{k+1},\lambda_{k+2})\leq-\gamma\|(x_{k},\lambda_{k+1},\lambda_{k+2})\|^{2}$.
Then

$\displaystyle\widehat{V}(x_{k+1},\lambda_{k+2})-\widehat{V}(x_{k},\lambda_{k+1})$
$\displaystyle\leq-\gamma\|(x_{k},\lambda_{k+1},\lambda_{k+2})\|^{2}\leq-\gamma\|(x_{k},\lambda_{k+1})\|^{2}$
$\displaystyle\leq-(\gamma/\|P\|)\widehat{V}(x_{k},\lambda_{k+1})$

where the second inequality follows from
$\|(x_{k},\lambda_{k+1})\|^{2}\leq\|(x_{k},\lambda_{k+1},\lambda_{k+2})\|^{2}$
and the third inequality from
$\widehat{V}(x_{k},\lambda_{k+1})\leq\|P\|\|(x_{k},\lambda_{k+1})\|^{2}$.
Thus, $\widehat{V}(x_{k},\lambda_{k+1})$ $\leq(1-\gamma/\|P\|)^{k}$
$\widehat{V}(x_{0},\lambda_{1})$ along all trajectories of DLCS from the
initial point $x_{0}$. Combining with
$\widehat{V}(x_{k},\lambda_{k+1})\geq\omega\|(x_{k},\lambda_{k+1}\|^{2}$
yields

$\|(x_{k},\lambda_{k+1})\|^{2}\leq(\|P\|/\omega)(1-\gamma/\|P\|)^{k}\|(x_{0},\lambda_{1})\|^{2}$

and the claim in (iii) follows. ∎

Note that setting $P_{x\lambda}=0$ and $P_{\lambda\lambda}=0$ in (15), we
obtain the sufficient conditions for stability with a CQLF. Theorem 2 shows
that sufficient conditions for exponential stability using an EQLF are to

Find $\displaystyle\;\;P\in\mathbb{S}^{n_{x}+n_{c}}$ (18a) such that
$\displaystyle\;P\succ_{{\cal K}}0,$ (18b)
$\displaystyle\;-\widehat{M}(P)\succ_{\widehat{{\cal K}}}0,$ (18c)

where $\widehat{M}(P)$ is defined in (15), and the set $\widehat{{\cal
K}}\subseteq\mathbb{R}^{n_{x}+2n_{c}}$ is defined as

$\displaystyle\widehat{{\cal K}}=$
$\displaystyle\;\\{(x,\lambda,\widehat{\lambda})\,|\,(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F),\,\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)\\}.$
(19)

###  3.3 Extension to Inhomogeneous DLCS

In this section, we consider a generalization of the DLCS, the autonomous
inhomogeneous DLCS which is given as

$\displaystyle x_{k+1}=$ $\displaystyle\;Ax_{k}+C\lambda_{k+1}+g$ (20a)
$\displaystyle 0\leq\lambda_{k+1}\perp$
$\displaystyle\;Dx_{k}+F\lambda_{k+1}+h\geq 0,$ (20b)

where $g\in\mathbb{R}^{n_{x}}$ and $h\in\mathbb{R}^{n_{c}}$.

In this section, we show how the stability results in §3.1 and §3.2 are
applicable to an inhomogeneous DLCS (20). For the results in this section, we
must make the stronger assumption that $F$ is a P-matrix.

A point $(x^{e},\lambda^{e})$ is said to be an equilibrium of (20) if

$\displaystyle x^{e}=$ $\displaystyle\;Ax^{e}+C\lambda^{e}+g$ (21a)
$\displaystyle 0\leq\lambda^{e}\perp$
$\displaystyle\;Dx^{e}+F\lambda^{e}+h\geq 0.$ (21b)

Let $(\alpha^{e},\beta^{e},\gamma^{e})$ denote the index partition according
to (6) at the equilibrium $x^{e}$. Suppose that $F$ is a P-matrix, and let
$\lambda(x)$ be the unique solution of the LCP($Dx+h$,$F$). From the global
Lipschitz continuity of $\lambda(x)$ (Lemma 1(b)), we have that there exists a
neighborhood ${\cal N}(x^{e})$ around $x^{e}$ such that for all $x\in{\cal
N}(x^{e})$, ${\left[\lambda(x)\right]}_{i}>0$ for all $i\in\alpha^{e}$ and
${\left[Dx+F\lambda(x)\right]}_{i}>0$ for all $i\in\gamma^{e}$. Hence,
$\alpha(\lambda(x))\supseteq\alpha^{e}$ and
$\gamma(\lambda(x))\supset\gamma^{e}$ holds for all $x\in{\cal N}(x^{e})$.

Thus, for all $x\in{\cal N}(x^{e})$ the $\text{LCP}(Dx,F)$ can be restated as

$\displaystyle 0$
$\displaystyle=\;{\left[Dx+F\lambda(x)+h\right]}_{\alpha^{e}}$ (22a)
$\displaystyle 0\leq{\left[\lambda(x)\right]}_{\beta^{e}}$
$\displaystyle\perp\;{\left[Dx+F\lambda(x)+h\right]}_{\beta^{e}}\geq 0$ (22b)
$\displaystyle{\left[\lambda(x)\right]}_{\gamma^{e}}$ $\displaystyle=\;0.$
(22c)

Introducing $\delta x_{k+1}=x_{k+1}-x^{e}$,
$\delta\lambda_{k+1}=\lambda_{k+1}-\lambda^{e}$, and using the observation in
(22), the DLCS in (20) can be written for all $x\in{\cal N}(x^{e})$ as

$\displaystyle\delta x_{k+1}=$ $\displaystyle\;A\delta
x_{k}+C\delta\lambda_{k+1}$ (23a) $\displaystyle 0=$
$\displaystyle\;{\left[D\delta
x_{k}+F\delta\lambda_{k+1}\right]}_{\alpha^{e}}$ (23b) $\displaystyle
0\leq{\left[\delta\lambda_{k+1}\right]}_{\beta^{e}}\perp$
$\displaystyle\;{\left[D\delta
x_{k}+F\delta\lambda_{k+1}\right]}_{\beta^{e}}\geq 0$ (23c) $\displaystyle 0=$
$\displaystyle\;{\left[\delta\lambda_{k+1}\right]}_{\gamma^{e}}.$ (23d)

The system (23) can be further simplified to the homogeneous DLCS

$\displaystyle\delta x_{k+1}=$ $\displaystyle\;\widehat{A}\delta
x_{k}+\widehat{C}{\left[\delta\lambda_{k+1}\right]}_{\beta^{e}}$ (24a)
$\displaystyle 0\leq{\left[\delta\lambda_{k+1}\right]}_{\beta^{e}}\perp$
$\displaystyle\;\widehat{D}\delta
x_{k}+\widehat{F}{\left[\delta\lambda_{k+1}\right]}_{\beta^{e}}\geq 0,$ (24b)

where
$\widehat{A}=(A-{\left[C\right]}_{\bullet\alpha^{e}}({\left[F\right]}_{\alpha^{e}\alpha^{e}})^{-1}{\left[D\right]}_{\alpha^{e}\bullet}$,
$\widehat{C}={\left[C\right]}_{\bullet\beta^{e}}-{\left[C\right]}_{\bullet\alpha^{e}}({\left[F\right]}_{\alpha^{e}\alpha^{e}})^{-1}{\left[F\right]}_{\alpha^{e}\beta^{e}}$,
$\widehat{D}={\left[D\right]}_{\beta^{e}\bullet}-{\left[F\right]}_{\beta^{e}\alpha^{e}}({\left[F\right]}_{\alpha^{e}\alpha^{e}})^{-1}{\left[D\right]}_{\alpha^{e}\bullet}$,
and
$\widehat{F}={\left[F\right]}_{\beta^{e}\beta^{e}}-{\left[F\right]}_{\beta^{e}\alpha^{e}}({\left[F\right]}_{\alpha^{e}\alpha^{e}})^{-1}{\left[F\right]}_{\alpha^{e}\beta^{e}}$.
Thus, the stability analysis of the previous sections can be applied to the
homogenous DLCS in (24) to verify local stability of $x^{e}$.

##  4 A Cutting-Plane Algorithm for Computing the Lyapunov Function

The sufficient conditions (12) and (18) for exponential stability using a CQLF
or an EQLF can be written in a general fashion as

Find $\displaystyle\;\mathbf{P}\in\mathbb{S}^{n}$ (25a) such that
$\displaystyle\;\mathbf{P}\succ_{{\cal K}_{1}}0,$ (25b)
$\displaystyle\;-\mathbf{M}(\mathbf{P})\succ_{{\cal K}_{2}}0$ (25c)

where ${\cal K}_{1}$ and ${\cal K}_{2}$ are closed (not necessarily convex)
cones and the matrix $\mathbf{M}\in\mathbb{S}^{m}$ depends on the matrix
$\mathbf{P}$. The feasibility conditions require that $\mathbf{P}$ is strict
${\cal K}_{1}$-copositive and $-\mathbf{M}$ is strict ${\cal
K}_{2}$-copositive.

A typical approach to handling the strict ${\cal K}_{1}$-copositivity and
strict ${\cal K}_{2}$-copositivity is to relax the conditions into linear
matrix inequalities using the S-Lemma [2, 3]. The transformation to such a
semidefinite formulation for the CQLF and EQLF is provided in Appendix A. Our
intent is to provide an exact algorithm which we describe next.

###  4.1 Cutting Plane Algorithm

We state by recalling a result from [7] that bounds the set of test vectors
for copositivity with respect to a cone.

######  Lemma 2.

[7, Lemma 1] Let $\|\cdot\|_{1}$ be the 1-norm on $\mathbb{R}^{n}$ and ${\cal
C}\subset\mathbb{R}^{n}$ be a closed cone. Then

$N\succcurlyeq_{{\cal C}}0\iff z^{T}Nz\geq 0\text{ for all }z\in{\cal
C},\|z\|_{1}=1.$

An analogous results holds for $N\succ_{\cal C}0$.

Thus, the feasibility problem (25) can be equivalently cast as the semi-
infinite optimization problem

$\displaystyle\max\limits_{\mu,\mathbf{P}\in\mathbb{S}^{n}}$
$\displaystyle\;\mu$ (26a) s.t.
$\displaystyle\;u^{T}\mathbf{P}u\,\,\,\,\geq\mu\,\forall\,u\in{\cal U}$ (26b)
$\displaystyle\;-v^{T}\mathbf{M}(\mathbf{P})v\geq\mu\,\forall\,v\in{\cal V}$
(26c)

where ${\cal U}={\cal K}_{1}\cap\mathbb{B}^{n}_{1}$ and ${\cal V}={\cal
K}_{2}\cap\mathbb{B}^{m}_{1}$. If the optimal objective value of (26) is
positive then the origin is exponentially stable for (1). The presence of
infinite number of constraints poses a difficulty for computation. We propose
to solve (26) using a cutting plane algorithm which alternates between a
_master problem_ and two _separation problems_. At an iteration $l$ of the
cutting plane algorithm the master problem attempts to find a matrix
$\mathbf{P}^{l}$ for which the inequalities (26b)-(26c) hold over finite sets
${\cal U}^{l},{\cal V}^{l}$. The separation problem verifies if the computed
$\mathbf{P}^{l}$ satisfies the inequalities in (26b)-(26c) for all vectors in
${\cal U}$ and ${\cal V}$. If satisfied, then $\mathbf{P}^{l}$ is a matrix for
which (25) holds. If not, then the separation problem identifies points
$u^{l},v^{l}$ for which the inequalities (25b)-(25c) fail to hold. These
points are appended to the sets ${\cal U}^{l},{\cal V}^{l}$, and the algorithm
continues.

The _master problem_ associated with (26) is

$\displaystyle\max\limits_{\mu,\mathbf{P}\in\mathbb{S}^{n}}$
$\displaystyle\;\mu$ (27a) s.t.
$\displaystyle\;u^{T}\mathbf{P}u\geq\mu\,\forall\,u\in{\cal U}^{l}$ (27b)
$\displaystyle\;-v^{T}\mathbf{M}(\mathbf{P})v\geq\mu\,\forall\,v\in{\cal
V}^{l}$ (27c) $\displaystyle\;-1\leq{\left[\mathbf{P}\right]}_{ij}\leq
1\,\forall\,i,j=1,\ldots,n.$ (27d)

The constraint (27d) is a normalization constraint and does not affect the
feasible region due to the homogeneity of (25). Specifically, if
$\hat{\mathbf{P}}$ satisfies (25) then
$\hat{\mathbf{P}}/(\max_{i,j}{\left[\hat{\mathbf{P}}\right]}_{ij})$ also
satisfies (25). From a computational standpoint (27d) serves to bound the
feasible region and ensures that an optimal solution to the linear program
(27) always exists as long as ${\cal U}^{l}\neq\emptyset$. Note that from (27)
we have that

$\displaystyle\mu=\max_{\mathbf{P}:|{\left[\mathbf{P}\right]}_{ij}|\leq
1}\min\left\\{\min_{u\in{\cal U}^{l}}u^{T}\mathbf{P}u,\min_{v\in{\cal
V}^{l}}-v^{T}\mathbf{M}(\mathbf{P})v\right\\}.$ (28)

Thus, if ${\cal U}^{l}\neq\emptyset$, an optimal solution to (27) must have
optimal value at most $n$, since

$\displaystyle\mu\leq$
$\displaystyle\;\max_{\mathbf{P}:|{\left[\mathbf{P}\right]}_{ij}|\leq
1}\max_{u\in{\cal
U}^{l}}u^{T}\mathbf{P}u\leq\max_{\mathbf{P}:|{\left[\mathbf{P}\right]}_{ij}|\leq
1,u\in\mathbb{B}^{n}_{1}}u^{T}\mathbf{P}u\leq\max_{u\in\mathbb{B}^{n}_{1}}\sum\limits_{ij}u_{i}u_{j}\leq
n.$

Since $\mathbf{P}=0,\mu=0$ is a feasible solution to (27), the master linear
program must have an optimal solution, which we denote as
$(\mu^{l},\mathbf{P}^{l})$.

The _separation problem_ associated with (27) for the cone ${\cal K}_{1}$ is

$\displaystyle\min\limits_{\nu_{1}\in\mathbb{R},u\in\mathbb{R}^{n}}$
$\displaystyle\;\nu_{1}$ (29a) s.t.
$\displaystyle\;u^{T}\mathbf{P}^{l}u\leq\nu_{1}$ (29b)
$\displaystyle\;u\in{\cal K}_{1}\cap\mathbb{B}^{n}_{1}.$ (29c)

The _separation problem_ associated with (27) for the cone ${\cal K}_{2}$ is

$\displaystyle\min\limits_{\nu_{2}\in\mathbb{R},v\in\mathbb{R}^{m}}$
$\displaystyle\;\nu_{2}$ (30a) s.t.
$\displaystyle\;-v^{T}\mathbf{M}(\mathbf{P}^{l})v\leq\nu_{2}$ (30b)
$\displaystyle\;v\in{\cal K}_{2}\cap\mathbb{B}^{m}_{1}.$ (30c)

In the separation problems in (29) and (30) we have included the ${\cal
K}_{1}\cap\mathbb{B}^{n}_{1}$ and ${\cal K}_{2}\cap\mathbb{B}^{m}_{1}$
constraint without specifying a formulation. We provide a mixed integer
programming formulation for the CQLF conditions (12) in § 4.2 and for the EQLF
conditions (18) in § 4.3.

Algorithm 1 summarizes the steps in the cutting plane algorithm. The steps of
solving the master and separation problems at each iteration of the algorithm
are described in Lines 5-18. The algorithm aims to find $\mathbf{P}$,
$u\in\mathbb{B}^{n}_{1}$, $v\in\mathbb{B}^{m}_{1}$ satisfying

$u^{T}\mathbf{P}u\geq\epsilon\ \forall u\in{\cal K}_{1},\
-v^{T}\mathbf{M}(\mathbf{P})v\geq\epsilon\ \forall v\in{\cal K}_{2},\
\eqref{sicpPgen-mp:c3}.$ (31)

By appropriate scaling, this is equivalent to

$\mathbf{P}\succ_{{\cal K}_{1}}0,\,-\mathbf{M}(\mathbf{P})\succ_{{\cal
K}_{2}}0,\,\eqref{sicpPgen-mp:c3}.$ (32)

Algorithm 1 A cutting plane algorithm to find $\mathbf{P}$ satisfying (25)

1: Input: $\epsilon>0$, MaxIter.

2: Output: $\mathbf{P}$ satisfying (25) if one exists.

3: Set $l=0$. Initialize the sets ${\cal U}^{0}\subset{\cal U}$ and ${\cal
V}^{0}\subset{\cal V}$.

4: while $l<\text{MaxIter}$ do

5: Solve (27) to obtain optimal solution $(\mu^{l},\mathbf{P}^{l})$.

6: if $\mu^{l}<\epsilon$ then

7: No solution exists satisfying user-specified tolerances. Terminate.

8: end if

9: Solve (29) and (30) using $\mathbf{P}^{l}$ to obtain optimal solutions
$(\nu_{1}^{l},u^{l})$ and $(\nu_{2}^{l},v^{l})$ respectively.

10: if $\min(\nu_{1}^{l},\nu_{2}^{l})\geq\epsilon$ then

11: Return $\mathbf{P}^{l}$ as a feasible solution. Terminate.

12: end if

13: if $\nu_{1}^{l}<\epsilon$ then

14: ${\cal U}^{l+1}\leftarrow{\cal U}^{l}\cup\\{u^{l}\\}$

15: end if

16: if $\nu_{2}^{l}<\epsilon$ then

17: ${\cal V}^{l+1}\leftarrow{\cal V}^{l}\cup\\{v^{l}\\}$

18: end if

19: Set $l\leftarrow l+1$.

20: end while

It is easy to see that if on solving (27) (in Line 5) $\mu^{l}<\epsilon$ then
(26) is infeasible.

######  Lemma 3.

Suppose $\mu^{l}<\epsilon$ at some iteration $l$ of Algorithm 1. Then there
exists no $\mathbf{P}$ in (27d) satisfying (31).

###### Proof.

If $\mu^{l}<\epsilon$ then by (28)

$\displaystyle\;\max_{\mathbf{P}:|{\left[\mathbf{P}\right]}_{ij}|\leq
1}\min\left\\{\min_{u\in{\cal U}^{l}}u^{T}\mathbf{P}u,\min_{v\in{\cal
V}^{l}}-v^{T}\mathbf{M}(\mathbf{P})v\right\\}<\epsilon$
$\displaystyle\implies$
$\displaystyle\;\max_{\mathbf{P}:|{\left[\mathbf{P}\right]}_{ij}|\leq
1}\min\left\\{\min_{u\in{\cal U}}u^{T}\mathbf{P}u,\min_{v\in{\cal
V}}-v^{T}\mathbf{M}(\mathbf{P})\right\\}<\epsilon$

where the second inequality follows by noting that the minimization is over
larger sets ${\cal U}^{l}\subseteq{\cal U}$ and ${\cal V}^{l}\subseteq{\cal
V}$. Hence, there exists no $\mathbf{P}$ satisfying (31), which proves the
claim. ∎

Thus, if the Algorithm 1 terminates, it either terminates with a
$\mathbf{P}\in\mathbb{S}^{n}$ such that $|{\left[\mathbf{P}\right]}_{ij}|\leq
1$, $\mathbf{P}\succcurlyeq_{{\cal K}_{1}}\epsilon$ and
$-\mathbf{M}(\mathbf{P})\succcurlyeq_{{\cal K}_{2}}\epsilon$, or a proof that
no such matrix $\mathbf{P}$ exists.

###  4.2 Common Quadratic Lyapunov Function

To state the sufficiency conditions (12) required for a CQLF in the general
form (25), we set $n=n_{x}$, $\mathbf{P}=P_{xx}$, ${\cal
K}_{1}=\mathbb{R}^{n_{x}}$, $m=(n_{x}+n_{c})$, $\mathbf{M}=M(P_{xx})$ in (10),
and ${\cal K}_{2}={\cal
K}=\text{Gr}\,\text{SOL}(Dx,F)\subseteq\mathbb{R}^{n_{x}+n_{c}}$. In the
separation problem for a CQLF, we are given a matrix $P_{xx}^{l}$, and we wish
to verify if the conditions (12) hold for this matrix. Thus, we seek a vector
$w=(x,\lambda)\in\mathbb{B}^{n_{x}+n_{c}}_{1}$ such that either
$w^{T}P_{xx}^{l}w\leq 0$ or $-w^{T}M(P_{xx}^{l})w\leq 0$ for some
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$. If no such vector exits, then we
have verified that $P_{xx}^{l}$ satisfies conditions (12), else the associated
vector may be added to the appropriate set ${\cal U}^{l}$ or ${\cal V}^{l}$ in
the master problem (27).

The separation problem for the constraint $P_{xx}\succ 0$ can be posed as the
following mixed integer quadratic program, which minimizes the value of the
quadratic form $x^{T}P_{xx}^{l}x$ over a 1-norm ball:

$\displaystyle\min\limits_{x,y,x^{+},x^{-},\nu}$ $\displaystyle\;\nu$ (33a)
$\displaystyle\text{s.t}\qquad x^{T}P_{xx}^{l}x$ $\displaystyle\leq\nu$ (33b)
$\displaystyle x$ $\displaystyle=x^{+}-x^{-}$ (33c)
$\displaystyle\mathbf{1}^{T}(x^{+}+x^{-})$ $\displaystyle=1$ (33d)
$\displaystyle 0\leq x^{+}$ $\displaystyle\leq y$ (33e) $\displaystyle 0\leq
x^{-}$ $\displaystyle\leq\mathbf{1}-y$ (33f) $\displaystyle y$
$\displaystyle\in\\{0,1\\}^{n_{x}},$ (33g)

where $\mathbf{1}$ is a vector of all ones.

The separation problem for the constraint $-M(P_{xx}^{l})\succcurlyeq_{{\cal
K}}0$ can be posed as the following mixed integer quadratic program in
variables $w=(x,\lambda)\in\mathbb{R}^{n_{x}+n_{c}}$, binary variables
$y\in\\{0,1\\}^{n_{x}}$, $z\in\\{0,1\\}^{n_{c}}$ and auxiliary variables
$x^{+},x^{-}\in\mathbb{R}^{n_{x}}$:

$\displaystyle\min\limits_{x,\lambda,y,x^{+},x^{-},\nu}$ $\displaystyle\;\nu$
(34a) $\displaystyle\text{s.t}\qquad-(x,\lambda)^{T}M(P_{xx}^{l})(x,\lambda)$
$\displaystyle\leq\nu$ (34b) $\displaystyle x$ $\displaystyle=x^{+}-x^{-}$
(34c) $\displaystyle\mathbf{1}^{T}(x^{+}+x^{-}+\lambda)$ $\displaystyle=1$
(34d) $\displaystyle 0\leq x^{+}$ $\displaystyle\leq y$ (34e) $\displaystyle
0\leq x^{-}$ $\displaystyle\leq\mathbf{1}-y$ (34f) $\displaystyle y$
$\displaystyle\in\\{0,1\\}^{n_{x}}$ (34g) $\displaystyle 0\leq\lambda$
$\displaystyle\leq z$ (34h) $\displaystyle
0\leq\sum_{j=1}^{n_{x}}{\left[D\right]}_{ij}{\left[x\right]}_{j}+\sum_{j=1}^{n_{c}}{\left[F\right]}_{ij}{\left[\lambda\right]}_{j}$
$\displaystyle\leq\theta_{i}(1-{\left[z\right]}_{i})\quad\forall
i=1,\ldots,n_{c}$ (34i) $\displaystyle z$ $\displaystyle\in\\{0,1\\}^{n_{c}},$
(34j)

where ${\bf 1}$ is a vector of all ones of conformal dimension, and the
constant

$\theta_{i}=\max\left(\max_{j}\left|{\left[D\right]}_{ij}\right|,\max_{k}\left({\left[F\right]}_{ik}\right)\right)$

provides an upper bound on the value of the left hand side of constraint (34i)
for any $(x,\lambda)\in\mathbb{B}^{n_{x}+n_{c}}_{1}$. The objective (34a) and
constraint (34b) seek to minimize the value of the (possibly nonconvex)
quadratic function $(x,\lambda)^{T}M(P_{xx}^{l})(x,\lambda)$. Constraints
(34c)—(34g) ensure that $(x,\lambda)\in\mathbb{B}^{n_{x}+n_{c}}_{1}$, and
constraints (34h)—(34j) enforce the conditions that
$(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$.

###  4.3 Extended Quadratic Lyapunov Function

To state the sufficiency conditions for the EQLF in Theorem 2 in terms of the
general setting of (25), we can set $n=(n_{x}+n_{c})$, $\mathbf{P}=P$, ${\cal
K}_{1}={\cal K}$, $m=(n_{x}+2n_{c})$, $\mathbf{M}=\widehat{M}(P)$, and ${\cal
K}_{2}=\widehat{{\cal K}}$ where $\widehat{{\cal
K}}=\\{(x,\lambda,\widehat{\lambda})\,|\,(x,\lambda)\in\text{Gr}\text{SOL}(Dx,F),\,\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)\\}$.
In the separation problem for an EQLF, we are given a matrix $P^{l}$, and we
seek a vector $v=(x,\lambda,\hat{\lambda})\in\mathbb{B}^{n_{x}+2n_{c}}_{1}$
such that either

$\displaystyle(i)\quad(x,\lambda)^{T}P^{l}(x,\lambda)<0$ $\displaystyle\mbox{
for some }(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F);\mbox{ or }$
$\displaystyle(ii)-(x,\lambda,\hat{\lambda})^{T}\widehat{M}(P^{l})(x,\lambda,\hat{\lambda})<0$
$\displaystyle\text{ for some }(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$
$\displaystyle\text{ and }\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F).$

The separation problem for the condition (i) ($P\succ_{{\cal K}}0$) is very
similar to (34) and thus not repeated here. The separation problem for
condition (ii) ($\widehat{M}(P)\succ_{\widehat{{\cal K}}}0$) can be posed as
the following quadratic mixed integer program:

$\displaystyle\min\limits_{x,\lambda,\hat{\lambda},y,z,w,x^{+},x^{-},\nu}$
$\displaystyle\;\nu$ (35a)
$\displaystyle\text{s.t}\qquad(x,\lambda,\hat{\lambda})^{T}\widehat{M}(P^{l})(x,\lambda,\hat{\lambda})$
$\displaystyle\leq\nu$ (35b) $\displaystyle x$ $\displaystyle=x^{+}-x^{-}$
(35c) $\displaystyle\mathbf{1}^{T}(x^{+}+x^{-}+\lambda+\hat{\lambda})$
$\displaystyle=1$ (35d) $\displaystyle 0\leq x^{+}$ $\displaystyle\leq y$
(35e) $\displaystyle 0\leq x^{-}$ $\displaystyle\leq\mathbf{1}-y$ (35f)
$\displaystyle y$ $\displaystyle\in\\{0,1\\}^{n_{x}}$ (35g) $\displaystyle
0\leq\lambda$ $\displaystyle\leq z$ (35h) $\displaystyle
0\leq\sum_{j=1}^{n_{x}}{\left[D\right]}_{ij}{\left[x\right]}_{j}+\sum_{j=1}^{n_{c}}{\left[F\right]}_{ij}{\left[\lambda\right]}_{j}$
$\displaystyle\leq\theta_{i}(1-{\left[z\right]}_{i})\quad\forall
i=1,\ldots,n_{c}$ (35i) $\displaystyle z$ $\displaystyle\in\\{0,1\\}^{n_{c}}$
(35j) $\displaystyle 0\leq\hat{\lambda}$ $\displaystyle\leq w$ (35k)
$\displaystyle
0\leq\sum_{j=1}^{n_{x}}{\left[DA\right]}_{ij}{\left[x\right]}_{j}+\sum_{j=1}^{n_{c}}{\left[DC\right]}_{ij}{\left[\lambda\right]}_{j}+\sum_{j=1}^{n_{c}}{\left[F\right]}_{ij}{\left[\hat{\lambda}\right]}_{j}$
$\displaystyle\leq\Theta_{i}(1-{\left[w\right]}_{i})\quad\forall
i=1,\ldots,n_{c}$ (35l) $\displaystyle w$ $\displaystyle\in\\{0,1\\}^{n_{c}},$
(35m)

where
$\Theta_{i}=\max\left(\max_{j}\left|{\left[DA\right]}_{ij}\right|,\max_{j}\left({\left[F\right]}_{ij}\right),\max_{j}\left({\left[DC\right]}_{ij}\right)\right)$
provides an upper bound on the value of the left-hand-side of (35l) for any
feasible solution to the problem. The objective (35a) and constraint (35b)
minimize the nonconvex quadratic function
$(x,\lambda,\hat{\lambda})^{T}\widehat{M}(P^{l})(x,\lambda,\hat{\lambda})$.
Constraints (35c)—(35g) ensure that
$(x,\lambda,\hat{\lambda}\in\mathbb{B}^{n_{x}+2n_{c}}_{1}$. The constraints
(35h)—(35j) enforce that $(x,\lambda)\in\text{Gr}\,\text{SOL}(Dx,F)$, and the
constraints (35k)— (35m) force that
$\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F).$

##  5 Numerical Experiments

We implemented Algorithm 1 for the computation of CQLF and EQLF in Python and
solved the master problems and separation problemss using Gurobi 9.0.2 [17].
In our experiments, we chose separation tolerance $\epsilon=10^{-6}$ in
Algorithm 1. In our implementation, we do not always solve the separation
problems for the CQLF and EQLF to optimality. We terminate the separation
early, if the solver finds a feasible solution whose objective value is less
than or equal to zero. In this case, the incumbent solution is a vector that
can be used to separate the current master problem solution and can be added
to the set of vectors (${\cal U}^{l}$ or ${\cal V}^{l}$) in the master problem
(27).

###  5.1 DLCS derived from LCS with $\widetilde{F}$ a P-matrix

We consider the 3 LCS examples from [8] (Examples 3.1-3.3) and another example
from Heemels, Schumacher and Weiland [19, §2].

  * •

cam31 is Example 3.1 [8] with data matrices:
$\widetilde{A}=\begin{pmatrix}1.0\end{pmatrix}$,
$\widetilde{C}=\begin{pmatrix}2&-2\end{pmatrix}$,
$\widetilde{D}=\begin{pmatrix}1\\\ -1\end{pmatrix}$,
$\widetilde{F}=\begin{pmatrix}1&3\\\ 0&1\end{pmatrix}$.

  * •

cam32 is Example 3.2 [8] with data matrices:
$\widetilde{A}=\begin{pmatrix}-1.0\end{pmatrix}$,
$\widetilde{C}=\begin{pmatrix}0&1\end{pmatrix}$,
$\widetilde{D}=\begin{pmatrix}1\\\ 1\end{pmatrix}$,
$\widetilde{F}=\begin{pmatrix}1&3\\\ 0&1\end{pmatrix}$.

  * •

cam33 is Example 3.3 [8] with data matrices:
$\widetilde{A}=\begin{pmatrix}-5&-4&0\\\ -1&-2&0\\\ 0&0&1\end{pmatrix}$,
$\widetilde{C}=\begin{pmatrix}-3&0&0\\\ -21&0&0\\\ 0&2&-2\end{pmatrix}$,
$\widetilde{D}=\begin{pmatrix}1&0&0\\\ 0&0&1\\\ 0&0&-1\end{pmatrix}$,
$\widetilde{F}=\begin{pmatrix}1&0&0\\\ 0&1&3\\\ 0&0&1\end{pmatrix}$.

  * •

hem2 is Example in §2 [19] with data matrices:
$\widetilde{A}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\ -2&1&0&0\\\
1&-1&0&0\end{pmatrix}$, $\widetilde{C}=\begin{pmatrix}0\\\ 0\\\ 1\\\
0\end{pmatrix}$, $\widetilde{D}=\begin{pmatrix}1&0&0&0\end{pmatrix}$,
$\widetilde{F}=\begin{pmatrix}1\end{pmatrix}$.

Note that in all the examples above the $\widetilde{F}$ matrix is a P-matrix.
In the following we consider DLCS that are obtained from explicit ($\theta=0$
in (3)) and implicit Euler ($\theta=1$ in (3)) time-stepping schemes. We
considered two possible time-steps of $\Delta t=0.1$ and $0.05$. For these
choices of time-stepping schemes and $\Delta t$, it can be verified that $F$
is also a P-matrix.

Table 1 summarizes the results for the CQLF computation while Table 2 provides
the results for the EQLF. Table 1 shows that the algorithm is successful in
identifying a CQLF such that $\mu\geq\epsilon$ for cam31 and cam32 for both
time-steps. The algorithm shows that a CQLF does not exist such that
$\mu\geq\epsilon$ for cam33 and hem2 when using the explicit Euler time-
stepping scheme and either time-step choices. When using the implicit Euler
scheme cam31, cam32 and cam33 permit identifying a CQLF for both time-step
choices. As for hem2, a CQLF does not exist for either time-step choices.
Table 2 shows that an EQLF with $\mu\geq 10^{-6}$ can be computed for all the
examples except for hem2 for either time-steps when using the explicit Euler.
When using the implicit Euler an EQLF is computed for all systems.

In all cases that are stabilized using CQLF (EQLF), we also simulated the
system from 100 random initial and verified that the CQLF (EQLF) is indeed
decreasing over the trajectories as shown in Theorems 1 and 2.

From the tables it can be seen that the best $\mu$ for different time-steps
(whenever feasible) is in proportion to the time-steps. The next section shows
that this is indeed true for the CQLF (see Lemma 5).

|  | $\Delta t=0.1$ | $\Delta t=0.05$
---|---|---|---
| Name | Best $\mu$ | # iters. | Time(s) | Best $\mu$ | # iters. | Time(s)
Explicit | cam31 | $4.348\cdot 10^{-2}$ | 1 | 0.006 | $2.326\cdot 10^{-2}$ | 1 | 0.006
cam32 | $1.0\cdot 10^{-1}$ | 1 | 0.003 | $5.0\cdot 10^{-2}$ | 1 | 0.003
cam33 | $<10^{-6}$ | 3 | 0.171 | $<10^{-6}$ | 10 | 0.477
hem2 | $<10^{-6}$ | 7 | 0.157 | $<10^{-6}$ | 6 | 0.118
Implicit | cam31 | $4.762\cdot 10^{-2}$ | 1 | 0.003 | $2.439\cdot 10^{-2}$ | 1 | 0.010
cam32 | $9.091\cdot 10^{-2}$ | 1 | 0.004 | $4.762\cdot 10^{-2}$ | 1 | 0.003
cam33 | $3.941\cdot 10^{-3}$ | 9 | 0.624 | $4.373\cdot 10^{-4}$ | 15 | 1.106
hem2 | $<10^{-6}$ | 22 | 0.720 | $<10^{-6}$ | 16 | 0.320
Table 1: Summary of results from applying Algorithm 1 for the CQLF computation with $\epsilon=10^{-6}$. Best $\mu$ \- refers to the best $\mu$ obtained at termination; # iters. - number of times the while-loop. in Algorithm 1 is executed; Tot. CPU time (s) - total time for algorithm. |  | $\Delta t=0.1$ | $\Delta t=0.05$
---|---|---|---
| Name | Best $\mu$ | # iters. | Time(s) | Best $\mu$ | # iters. | Time(s)
Explicit | cam31 | $4.188\cdot 10^{-2}$ | 2 | 0.048 | $4.133\cdot 10^{-2}$ | 1 | 0.025
cam32 | $1.766\cdot 10^{-1}$ | 1 | 0.014 | $4.870\cdot 10^{-2}$ | 1 | 0.012
cam33 | $4.859\cdot 10^{-3}$ | 9 | 1.630 | $1.112\cdot 10^{-2}$ | 13 | 2.500
hem2 | $<10^{-6}$ | 8 | 0.365 | $<10^{-6}$ | 6 | 0.321
Implicit | cam31 | $8.453\cdot 10^{-2}$ | 1 | 0.033 | $4.334\cdot 10^{-2}$ | 1 | 0.029
cam32 | $1.607\cdot 10^{-1}$ | 1 | 0.014 | $8.453\cdot 10^{-2}$ | 1 | 0.017
cam33 | $1.2\cdot 10^{-2}$ | 4 | 0.990 | $6.723\cdot 10^{-3}$ | 8 | 1.501
hem2 | $9.629\cdot 10^{-5}$ | 35 | 7.658 | $2.831\cdot 10^{-5}$ | 33 | 23.152
Table 2: Summary of results from applying Algorithm 1 for the EQLF computation
with $\epsilon=10^{-6}$. Best $\mu$ \- refers to the best $\mu$ obtained at
termination; # iters. - number of times the while-loop. in Algorithm 1 is
executed; Tot. CPU time (s) - total time for algorithm.

###  5.2 DLCS with $F$ a Q-matrix and R0-matrix

Stability results in the literature assume that $F$ is a $P$-matrix [2, 8],
which ensures that the state trajectory of the system is unique. The
sufficient conditions for stability of DLCS given in Section 3 assume the
weaker condition that $F$ is a $Q$-matrix and an $R_{0}$ matrix. Here, we
demonstrate that the cutting-plane method can find stability conditions for
these more general systems by considering the DLCS

$A=\begin{pmatrix}0.5&0.25\\\ -0.25&0.5\end{pmatrix},\,C=\begin{pmatrix}3\\\
5\end{pmatrix},\,D=\begin{pmatrix}1.0&0.0\\\
0.0&0.0\end{pmatrix},\,F=\begin{pmatrix}1&-1\\\ 1&0\end{pmatrix}.$ (36)

Every principal minor of the matrix $F$ is nonnegative, and hence $F$ is a
P0-matrix. Further, $F$ is a R0-matrix since $\text{SOL}(0,F)=\\{0\\}$. Thus,
$F$ is also a Q-matrix. Running our cutting-plane algorithm proves that the
system cannot be stabilized using a CQLF. However, the system in (36) can be
stabilized using the EQLF with $\mu=4.057\cdot 10^{-3}$ and

$P=\begin{pmatrix}5.33238468\cdot 10^{-1}&-1.79873472\cdot
10^{-1}&-1.0&-4.13688777\cdot 10^{-2}\\\ -1.79873472\cdot
10^{-1}&1.70818637\cdot 10^{-1}&2.78264033\cdot 10^{-1}&-8.65180851\cdot
10^{-4}\\\ -1.0&2.78264033\cdot 10^{-1}&1.0&-1.0\\\ -4.13688777\cdot
10^{-2}&-8.65180851\cdot 10^{-4}&-1.0&7.53616763\cdot 10^{-2}\end{pmatrix}.$

##  6 LCS and DLCS Stability

Çamlibel, Pang and Shen [8, Theorem 3.1] present sufficient conditions for
stability using an EQLF for the LCS (2). In this section, we consider the DLCS
(1) obtained from the discretization of the LCS (2) for some $\theta\in[0,1]$
in (4). Recall that by defining the matrix
$\widetilde{A}_{\theta}=(I_{n_{x}}-\theta\Delta t\widetilde{A})^{-1}$, the
time-stepping formulation of the LCS (4) can be identified with the DLCS as

$\displaystyle A=I_{n_{x}}+\Delta
t\widetilde{A}_{\theta}^{-1}\widetilde{A},\,$ $\displaystyle C=\Delta
t\widetilde{A}_{\theta}^{-1}\widetilde{C},\,$ (37) $\displaystyle
D=\widetilde{D}(I_{n_{x}}+\Delta
t\widetilde{A}_{\theta}^{-1}\widetilde{A}),\,$ $\displaystyle
F=\widetilde{F}+\Delta
t\widetilde{D}\widetilde{A}_{\theta}^{-1}\widetilde{C}.$

The stability results in [8, Theorem 3.1] are obtained under the assumption
that $\widetilde{F}$ is a P-matrix. Accordingly, we assume in the rest of the
section that $\widetilde{F}$ is a P-matrix. The main result presented in this
section is to show that the LCS is exponentially stable using a CQLF (EQLF) if
and only if the DLCS is exponentially stable using a CQLF (EQLF) for all
$\Delta t$ sufficiently small and any $\theta\in[0,1]$.

###  6.1 CQLF Equivalence

We will begin by stating the conditions for exponential stability of LCS using
CQLF. Note that [8, Theorem 3.1] states the asymptotic stability of LCS using
EQLF. Just as in the DLCS (see remarks following Theorem 2) we can obtain the
stability results for the CQLF after setting certain matrices in the Lyapunov
function to zero.

######  Definition 2.

[8, Theorem 3.1] The LCS (2) is asymptotically stable using the CQLF if there
exists $\widetilde{P}_{xx}\in\mathbb{S}^{n_{x}}$ such that

$\displaystyle\widetilde{P}_{xx}\succ 0,\mbox{ and }$ (38a)
$\displaystyle-\underbrace{\begin{pmatrix}\widetilde{A}^{T}\widetilde{P}_{xx}+\widetilde{P}_{xx}\widetilde{A}&\widetilde{P}_{xx}\widetilde{C}\\\
\widetilde{C}^{T}\widetilde{P}_{xx}&0\end{pmatrix}}_{=:Q(\widetilde{P}_{xx})}\succ_{{\cal
C}}0,$ (38b) where the set ${{\cal C}}\subseteq\mathbb{R}^{n_{x}+n_{c}}$ is
defined as ${\cal C}=\text{Gr}\,\text{SOL}(\widetilde{D}x,\widetilde{F}).$
(38c)

In the following we assume that $\Delta t$ is chosen such that

$\Delta t\cdot\theta\cdot\|\widetilde{A}\|<1,\,\Delta
t\cdot\|\widetilde{A}_{\theta}^{-1}\widetilde{A}\|<1,\text{ and }F\text{ is a
Q-matrix and R${}_{0}$-matrix}.$ (39)

Note that the first inequality in (39) ensures that
$\widetilde{A}_{\theta}^{-1}$ exists while the second inequality ensures that
$A$ is invertible. Both of these conditions can be ensured for all $\Delta t$
sufficiently small. Finally since we assume that $\widetilde{F}$ is a P-matrix
then there must exist $\Delta t$ sufficiently small so that
$F=\widetilde{F}+\Delta
t\widetilde{D}\widetilde{A}^{-1}_{\theta}\widetilde{C}$ is also a P-matrix, so
it is also both a Q-matrix and an R0-matrix. For $\Delta t$ satisfying (39) we
have by (37) that

$A^{-1}=(I_{n_{x}}+\Delta
t\cdot\widetilde{A}_{\theta}^{-1}\widetilde{A})^{-1}=I_{n_{x}}-\Delta
t\cdot\widetilde{A}+\Delta t^{2}\cdot\widetilde{A}_{r}$ (40)

where $\widetilde{A}_{r}$ is a bounded, residual term.

We will begin the analysis by relating the cones ${\cal K}$ and ${\cal C}$.

######  Lemma 4.

Suppose $\widetilde{F}$ is a P-matrix, let $\Delta t$ be such that (39) holds,
and let the matrix $N$ be defined as

$\displaystyle{N}=\begin{pmatrix}A&C\\\ 0&I_{n_{c}}\end{pmatrix}.$ (41)

Then the following hold:

  1. (a)

If $(x,\lambda)\in{\cal K}$ then $(Ax+C\lambda,\lambda)\in{\cal C}$.

  2. (b)

If $(x,\lambda)\in{{\cal C}}$ then $(\widetilde{x},\lambda)\in{\cal K}$ where
$\widetilde{x}=A^{-1}x-A^{-1}C\lambda$.

  3. (c)

There exist constants $c_{1},c_{2}>0$ such that
$c_{1}\|(\widetilde{x},\lambda)\|\leq\|(x,\lambda)\|\leq
c_{2}\|(\widetilde{x},\lambda)\|$.

###### Proof.

Let $(x,\lambda)\in\mathbb{R}^{n_{x}+n_{c}}$ and define
$\widetilde{x}=Ax+C\lambda$. From the definition of (41), we have that

$\displaystyle\begin{pmatrix}\widetilde{x}\\\
\lambda\end{pmatrix}=N\begin{pmatrix}x\\\ \lambda\end{pmatrix}.$ (42)

The matrix $N$ is invertible with $N^{-1}=\begin{pmatrix}A^{-1}&-A^{-1}C\\\
0&I_{n_{c}}\end{pmatrix}$. The claims in (a) and (b) follow immediately. From
the upper triangular representation of $N$, we have that the eigenvalues of
$N$ are the eigenvalues of $A$ and $I_{n_{c}}$. Hence,
$\|N\|\leq\max(\|A\|,1)$ and $\|N^{-1}\|\leq\max(\|A^{-1}\|,1)$. Combining
this with (42) yields the claim in (c). ∎

We now present an alternative representation of the matrix $M(P_{xx})$ (10) in
terms of $Q(\widetilde{P}_{xx})$ in (38b).

######  Lemma 5.

Suppose $\Delta t$ is sufficiently small so that (39) holds. Let
$P_{xx}\in\mathbb{S}^{n_{x}}$ be given. Then

$M(P_{xx})=N^{T}\left(\Delta t\cdot Q(P_{xx})+\Delta t^{2}\cdot
Q_{r}(P_{xx})\right)N,$ (43)

where $Q_{r}(P_{xx})$ is a bounded matrix and $N=\begin{pmatrix}A&C\\\
0&I_{n_{c}}\end{pmatrix}$.

###### Proof.

Consider the definition of the matrix $M$ in (10)

$\displaystyle M(P_{xx})=$
$\displaystyle\;N^{T}N^{-T}\begin{pmatrix}A^{T}P_{xx}A-P_{xx}&A^{T}P_{xx}C\\\
C^{T}P_{xx}A&C^{T}P_{xx}C\end{pmatrix}N^{-1}N$ (44) $\displaystyle=$
$\displaystyle\;N^{T}\begin{pmatrix}P_{xx}-A^{-T}P_{xx}A^{-1}&A^{-T}P_{xx}A^{-1}C\\\
C^{T}A^{-T}P_{xx}A^{-1}&-C^{T}A^{-T}P_{xx}A^{-1}C\end{pmatrix}N,$

where the final equality is obtained by substituting for $N^{-1}$, multiplying
through, and simplifying. Substituting (37) into the blocks of the matrix in
(44) and simplifying, we obtain for all $\Delta t$ satisfying (39) that

$\displaystyle P_{xx}-A^{-T}P_{xx}A^{-1}=$ $\displaystyle\;\Delta
t(\widetilde{A}^{T}P_{xx}+P_{xx}\widetilde{A})+\Delta t^{2}\cdot
Q_{r,11}(P_{xx})$ (45a) $\displaystyle A^{-T}P_{xx}A^{-1}C=$
$\displaystyle\;\Delta tP_{xx}\widetilde{C}+\Delta t^{2}\cdot
Q_{r,12}(P_{xx})$ (45b) $\displaystyle-C^{T}A^{-T}P_{xx}A^{-1}C=$
$\displaystyle\;\Delta t^{2}\cdot Q_{r,22}(P_{xx})$ (45c)

where $Q_{r,11}(P_{xx})$, $Q_{r,12}(P_{xx})$, and $Q_{r,22}(P_{xx})$ are
bounded residual matrices. Define the matrix
$Q_{r}(P_{xx})=\begin{pmatrix}Q_{r,11}(P_{xx})&Q_{r,12}(P_{xx})\\\
Q_{r,12}(P_{xx})^{T}&Q_{r,22}(P_{xx})\end{pmatrix}$. Plugging (45) into the
expressions in (44) obtains (43), which proves the claim. ∎

Combining the Lemmas (4) and (5) yields the main result that there exists a
CQLF for the LCS if and only if there exists a CQLF for DLCS.

######  Theorem 3.

Suppose $\widetilde{F}$ is a P-matrix and $\Delta t$ is such that (39) holds.
Then the LCS is exponentially stable using a CQLF if and only if the DLCS
obtained by the time-stepping formulation (4) is exponentially stable using a
CQLF.

###### Proof.

Let $P_{xx}\in\mathbb{S}^{n_{x}}$ be given. Pre-multiply and post-multiply the
equation in (43) by $(x,\lambda)\in{\cal K}$ to obtain

$\begin{pmatrix}x\\\ \lambda\end{pmatrix}^{T}M(P_{xx})\begin{pmatrix}x\\\
\lambda\end{pmatrix}=\begin{pmatrix}\widetilde{x}\\\
\lambda\end{pmatrix}^{T}\left(\Delta t\cdot Q(P_{xx})+\Delta t^{2}\cdot
Q_{r}(P_{xx})\right)\begin{pmatrix}\widetilde{x}\\\ \lambda\end{pmatrix},$
(46)

where $\widetilde{x}=Ax+C\lambda$. Further, by Lemma 4, we have that
$(\widetilde{x},\lambda)\in{\cal C}$.

Consider the only if part of the claim. Suppose the LCS is exponentially
stable using a CQLF. Since $-Q(P_{xx})\succ_{\cal C}0$ we have that
$-(\widetilde{x},\lambda)^{T}Q(P_{xx})(\widetilde{x},\lambda)\geq\gamma\|(\widetilde{x},\lambda)\|^{2}$
for some $\gamma>0$. Choose $\Delta t$ smaller than dictated by (39) so that
$\Delta t\cdot\|Q_{r}(P_{xx})\|\leq\gamma/2$. Then we have for such $\Delta t$
from (46) that $-(x,\lambda)^{T}M(P_{xx})(x,\lambda)\geq(\gamma\Delta
t)/2\|(\widetilde{x},\lambda)\|^{2}$. Combining this with Lemma 4(c) proves
the claim. The if part of the claim can be shown using similar arguments. This
completes the proof. ∎

###  6.2 EQLF Equivalence

We now consider the the EQLF and derive similar results between the LCS and
DLCS. The condition for exponential stability derived in [8] is stated next.

######  Definition 3.

[8, Theorem 3.1] The LCS (2) is exponentially stable if there exists
$\widetilde{P}_{xx}\in\mathbb{S}^{n_{x}}$,
$\widetilde{P}_{x\lambda}\in\mathbb{R}^{n_{x}\times n_{c}}$,
$\widetilde{P}_{\lambda\lambda}\in\mathbb{S}^{n_{c}}$ such that

$\displaystyle\underbrace{\begin{pmatrix}\widetilde{P}_{xx}&\widetilde{P}_{x\lambda}\\\
\widetilde{P}_{x\lambda}^{T}&\widetilde{P}_{\lambda\lambda}\end{pmatrix}}_{=:\widetilde{P}}\succ_{{\cal
K}}0$ (47a)
$\displaystyle-\underbrace{\begin{pmatrix}\widetilde{A}^{T}\widetilde{P}_{xx}+\widetilde{P}_{xx}\widetilde{A}&\widetilde{P}_{xx}\widetilde{C}+\widetilde{A}^{T}\widetilde{P}_{x\lambda}&\widetilde{P}_{x\lambda}\\\
\widetilde{C}^{T}\widetilde{P}_{xx}+\widetilde{P}_{x\lambda}^{T}\widetilde{A}&\widetilde{P}_{x\lambda}^{T}\widetilde{C}+\widetilde{C}^{T}\widetilde{P}_{x\lambda}&\widetilde{P}_{\lambda\lambda}\\\
\widetilde{P}_{x\lambda}^{T}&\widetilde{P}_{\lambda\lambda}&0\end{pmatrix}}_{=:\widehat{Q}(\tilde{P})}\succ_{{\cal
C}}0$ (47b) where the set $\widehat{{\cal
C}}\subseteq\mathbb{R}^{n_{x}+2n_{c}}$ is defined as
$\displaystyle\widehat{{\cal
C}}=\left\\{(\widetilde{x},\lambda,\lambda^{\prime})\in\mathbb{R}^{n_{x}+2n_{c}}\,\left|\,\begin{aligned}
&(\widetilde{x},\lambda)\in\text{Gr}\,\text{SOL}(\widetilde{D}\widetilde{x},\widetilde{F})\\\
&{\left[\widetilde{D}(\widetilde{A}\widetilde{x}+\widetilde{C}\lambda)+\widetilde{F}\lambda^{\prime}\right]}_{\alpha(\lambda)}=0\\\
&0\leq{\left[\lambda^{\prime}\right]}_{\beta(\lambda)}\perp{\left[\widetilde{D}(\widetilde{A}\widetilde{x}+\widetilde{C}\lambda)+\widetilde{F}\lambda^{\prime}\right]}_{\beta(\lambda)}\geq
0\\\
&{\left[\lambda^{\prime}\right]}_{\gamma(\lambda)}=0\end{aligned}\right.\right\\}$
(47c)

with $\alpha(\lambda),\beta(\lambda),\gamma(\lambda)$ defined according to (6)
for $\text{SOL}(\widetilde{D}\widetilde{x},\widetilde{F})$.

Note that $\lambda^{\prime}$ in (47c) is the directional derivative of
$\lambda(\widetilde{x})$ along the direction
$(\widetilde{A}\widetilde{x}+\widetilde{C}\lambda)$.

Suppose $(x,\lambda,\widehat{\lambda})\in\widehat{{\cal K}}$ with
$\widehat{\lambda}\in\text{SOL}(D(Ax+C\lambda),F)$. By Lemma 4 we have that
$\lambda=\text{SOL}(\widetilde{D}\widetilde{x},\widetilde{F})$. In other
words, $\widehat{\lambda}$ is the one-time step evolution of the
complementarity variables from the initial state $\widetilde{x}$ using the
time-stepping DLCS (3). The next result relates $\widehat{\lambda}$ with
$\lambda$ and $\lambda^{\prime}$. The proof relies on results from Shen and
Pang [37] and Chen and Wang [13] which we state next.

######  Lemma 6.

([37, Lemma 14]) Suppose $\widetilde{F}$ is a P-matrix. Let
$\widetilde{D}\widetilde{x}\neq 0$ for all $\widetilde{x}\neq 0$. Let
$\lambda=\text{SOL}(\widetilde{D}\widetilde{x},\widetilde{F})$ and
$\lambda(\Delta t)$ denote the complementarity variable obtained by evolving
the LCS (2) over $\Delta t$ from initial state $\widetilde{x}$. Then

$\lambda(\Delta t)=\lambda+\Delta
t\cdot\lambda^{\prime}+h_{1}(\widetilde{x})\cdot o(\Delta t)$ (48)

where the function $h(\widetilde{x})$ satisfies
$\|h_{1}\|=O(\|\widetilde{x}\|)$.

###### Proof.

The claim follows from Lemma 14 in [37]. ∎

######  Lemma 7.

([13, Theorem 3.1]) Suppose $\widetilde{F}$ is a P-matrix. Let $\lambda(\Delta
t)$ denote the complementarity variable resulting from evolving the LCS (2)
from the initial state $\widetilde{x}$ over a time duration of $\Delta t$ and
$\lambda=\text{SOL}(\widetilde{D}\widetilde{x},\widetilde{F})$. Let
$\widehat{\lambda}$ be the one time-step evolution of (3) from the same
initial state $\widetilde{x}$. Then there exists a $\overline{\Delta t}>0$
such that for all $\Delta t<\overline{\Delta t}$

$\|\widehat{\lambda}-\lambda(\Delta t)\|\leq\Delta t^{2}\cdot
O(\|\widetilde{x}\|).$ (49)

###### Proof.

Under the assumptions of the lemma it is readily verified that Theorem 3.1 in
[13] holds. Hence, $\|\widehat{\lambda}-\lambda(\Delta t)\|_{\infty}\leq\Delta
t^{2}L\|\widetilde{A}\widetilde{x}+\widetilde{C}\lambda\|_{\infty}$ for some
constant $L>0$ independent of $\Delta t$, $x,\lambda$. Using Lemma 1 we know
that $\|\lambda\|$ can be bounded using $\|\widetilde{x}\|$. Combining this
with the above and using the equivalence of norms yields (49). ∎

Combining Lemma 6 and 7 gives a relation between $\widehat{\lambda}$ and the
pair $\lambda$, $\lambda^{\prime}$.

######  Lemma 8.

Suppose the assumptions of Lemmas 6 and 7 hold. Then

$\widehat{\lambda}=\lambda+\Delta t\cdot\lambda^{\prime}+o(\Delta t)\cdot
h_{2}(\widetilde{x})$ (50)

where $\|h_{2}(\widetilde{x})\|\leq O(\|\widetilde{x}\|)$.

###### Proof.

From Lemma 7 we have that $\widehat{\lambda}=\lambda(\Delta t)+\Delta
t^{2}\cdot h_{2}(\widetilde{x})$ where
$\|h_{2}(\widetilde{x})\|=O(\|\widetilde{x}\|)$. Substituting for
$\lambda(\Delta t)$ from (48) yields (50). ∎

We next must relate the set $\widehat{{\cal K}}$ in (19) with $\widehat{{\cal
C}}$.

######  Lemma 9.

Suppose $\widetilde{F}$ is a P-matrix and $\widetilde{D}\widetilde{x}\neq 0$
for all $\widetilde{x}\neq 0$. Let the matrix $\widehat{N}$ and vector
$\widehat{h}$ be defined as

$\displaystyle\widehat{N}=\begin{pmatrix}A&C&0\\\ 0&I_{n_{c}}&0\\\
0&-\frac{1}{\Delta t}I_{n_{c}}&\frac{1}{\Delta t}I_{n_{c}}\end{pmatrix}\text{
and }\widehat{h}(\widetilde{x})=\begin{pmatrix}0\\\ 0\\\
h_{2}(\widetilde{x})\end{pmatrix}.$ (51a) Then there exists a
$\overline{\Delta t}>0$ such that for all $\Delta t\in(0,\overline{\Delta t})$
the following hold:

  1. (a)

If $(\widetilde{x},\lambda,\lambda^{\prime})\in\widehat{{\cal C}}$ then
$\widehat{N}^{-1}\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}+o(\Delta
t)\widehat{h}(\widetilde{x})\in\widehat{{\cal K}}$.

  2. (b)

If $(x,\lambda,\widehat{\lambda})\in\widehat{{\cal K}}$ then
$\widehat{N}\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix}-\frac{o(\Delta t)}{\Delta
t}\widehat{h}(Ax+C\lambda)\in\widehat{{\cal C}}$.

  3. (c)

There exist constants $c_{3},c_{4}>0$ such that $c_{3}\Delta
t\|(\widetilde{x},\lambda,\lambda^{\prime})\|\leq\|(x,\lambda,\widehat{\lambda})\|\leq
c_{4}\|(\widetilde{x},\lambda,\lambda^{\prime})\|$.

###### Proof.

From Lemma 4 we have that
$(A^{-1}\widetilde{x}-A^{-1}C\lambda,\lambda)\in{\cal K}$. Combining this with
Lemma 8 and noting that $\widehat{N}^{-1}=\begin{pmatrix}A^{-1}&-A^{-1}C&0\\\
0&I_{n_{c}}&0\\\ 0&I_{n_{c}}&\Delta t\cdot I_{n_{c}}\end{pmatrix}$ yields
$\displaystyle\begin{pmatrix}x\\\ \lambda\\\ \widehat{\lambda}\end{pmatrix}=$
$\displaystyle\;\widehat{N}^{-1}\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}+o(\Delta t)\widehat{h}(\widetilde{x})$ (52a)
$\displaystyle\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}=$ $\displaystyle\;\widehat{N}\begin{pmatrix}x\\\
\lambda\\\ \widehat{\lambda}\end{pmatrix}-\frac{o(\Delta t)}{\Delta
t}\widehat{h}(\widetilde{x}).$ (52b) where the second equality is obtained by
left multiplying (52a) with $\widehat{N}$. The claims in (a) and (b) follow
from (52a) and (52b) respectively. To prove (c), note that the first and
second rows and columns of $\widehat{N}$ can be swapped to produce a lower
triangular matrix $\begin{pmatrix}I_{n_{c}}&0&0\\\ C&A&0\\\ -\frac{1}{\Delta
t}I_{n_{c}}&0&\frac{1}{\Delta t}I_{n_{c}}\end{pmatrix}$. Based on the lower
triangular representation, the eigenvalues of $\widehat{N}$ can be obtained as
the eigenvalues of $A$, $I_{n_{c}}$, and $(1/\Delta t)I_{n_{c}}$. Hence,
$\|\widehat{N}\|\leq 1/(\Delta t)$ and
$\|\widehat{N}^{-1}\|\leq\max(\|A^{-1}\|,1)$. By choosing $\Delta t$ to be
smaller if necessary we have from (52a) and (52b) that
$\|(x,\lambda,\widehat{\lambda})\|\leq
2\|\widehat{N}^{-1}\|\|(\widetilde{x},\lambda,\lambda^{\prime})\|$ and
$\|(\widetilde{x},\lambda,\lambda^{\prime})\|\leq
2\|\widehat{N}\|\cdot\|(x,\lambda,\widehat{\lambda})\|$. Substituting for
$\|\widehat{N}\|$ and $\|\widehat{N}^{-1}\|$ yields the claim in (c). This
completes the proof.

∎

The matrices $\widehat{Q}$ and $\widehat{M}$ are related next.

######  Lemma 10.

Suppose $\Delta t$ is chosen to be sufficiently small such that (39) holds.
Let $P_{xx}\in\mathbb{S}^{n_{x}}$, $P_{x\lambda}\in\mathbb{R}^{n_{x}\times
n_{c}}$, and $P_{\lambda\lambda}\in\mathbb{S}^{n_{c}}$ be given. Then

$\widehat{M}(P)=\widehat{N}^{T}\left(\Delta t\cdot\widehat{Q}(P)+\Delta
t^{2}\cdot\widehat{Q}_{r}(P)\right)\widehat{N}$ (53)

where $\widehat{Q}_{r}(P)$ is a bounded matrix and $\widehat{N}$ is as defined
in (51a).

###### Proof.

Consider the definition of the matrix $\widehat{M}(P)$ in (15)

$\displaystyle\widehat{M}(P)=$
$\displaystyle\;\widehat{N}^{T}\widehat{N}^{-T}\begin{pmatrix}A^{T}P_{xx}A-P_{xx}&A^{T}P_{xx}C-P_{x\lambda}&A^{T}P_{x\lambda}\\\
C^{T}P_{xx}A-P_{x\lambda}^{T}&C^{T}P_{xx}C-P_{\lambda\lambda}&C^{T}P_{x\lambda}\\\
P_{x\lambda}^{T}A&P_{x\lambda}^{T}C&P_{\lambda\lambda}\end{pmatrix}\widehat{N}^{-1}\widehat{N}$
(54) $\displaystyle=$
$\displaystyle\;\widehat{N}^{T}\begin{pmatrix}P_{xx}-A^{-T}P_{xx}A^{-1}&U_{x\lambda}&\Delta
tP_{x\lambda}\\\ U_{x\lambda}^{T}&U_{\lambda\lambda}&\Delta
tP_{\lambda\lambda}\\\ \Delta tP_{x\lambda}^{T}&\Delta
tP_{\lambda\lambda}&0\end{pmatrix}\widehat{N},$ where
$\displaystyle\;U_{x\lambda}=A^{-T}P_{xx}A^{-1}C-A^{-T}P_{x\lambda}+P_{x\lambda}\text{
and }$
$\displaystyle\;U_{\lambda\lambda}=-C^{T}A^{-T}P_{xx}A^{-1}C+P_{x\lambda}^{T}A^{-1}C+C^{T}A^{-T}P_{x\lambda}.$

The second equality is obtained from the first by multiplying through and
simplifying. Substituting (37) into the blocks of the matrix in (54) and
simplifying, we obtain for all $\Delta t$ satisfying (39) that

$\displaystyle P_{xx}-A^{-T}P_{xx}A^{-1}=$ $\displaystyle\;\Delta
t(\widetilde{A}^{T}P_{xx}+P_{xx}\widetilde{A})+\Delta
t^{2}\cdot\widehat{Q}_{r,11},$ (55a) $\displaystyle U_{x\lambda}=$
$\displaystyle\;\Delta tP_{xx}\widetilde{C}+\Delta
t^{2}\cdot\widehat{Q}_{r,12},\text{ and}$ (55b) $\displaystyle
U_{\lambda\lambda}=$ $\displaystyle\;\Delta
t\left(P_{x\lambda}^{T}D+D^{T}P_{x\lambda}\right)+\Delta
t^{2}\cdot\widehat{Q}_{r,22},$ (55c)

where $\overline{M}_{r,11}$, $\overline{M}_{r,12}$, and $\overline{M}_{r,22}$
are bounded residual matrices. Define the matrix
$\overline{M}_{r}=\begin{pmatrix}\overline{M}_{r,11}&\overline{M}_{r,12}&0\\\
\overline{M}_{r,12}^{T}&\overline{M}_{r,22}&0\\\ 0&0&0\end{pmatrix}$. Plugging
(55) into the expressions in (54) yields (53), which proves the claim. ∎

We are now ready to show the equivalence between the EQLF conditions for the
LCS (47) and the DLCS (18).

######  Theorem 4.

Suppose $\widetilde{F}$ is a P-matrix. Let $\widetilde{D}\widetilde{x}\neq 0$
for all $\widetilde{x}\neq 0$. Then for all $\Delta t$ sufficiently small, the
LCS is exponentially stable using a EQLF if and only if the DLCS obtained by
the time-stepping formulation (3) is exponentially stable using a EQLF.

###### Proof.

Pre-multiply and post-multiply the equation in (53) by
$(x,\lambda,\widehat{\lambda})\in\widehat{{\cal K}}$ to obtain

$\displaystyle\;\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix}^{T}\widehat{M}\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\;\left(\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}+\frac{o(\Delta t)}{\Delta
t}\widehat{h}(\widetilde{x})\right)^{T}\left(\Delta t\cdot\widehat{Q}+\Delta
t^{2}\cdot\widehat{Q}_{r}\right)\left(\begin{pmatrix}\widetilde{x}\\\
\lambda\\\ \lambda^{\prime}\end{pmatrix}+\frac{o(\Delta t)}{\Delta
t}\widehat{h}(\widetilde{x})\right)$ $\displaystyle=$ $\displaystyle\;\Delta
t\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}^{T}\widehat{Q}\begin{pmatrix}\widetilde{x}\\\
\lambda\\\ \lambda^{\prime}\end{pmatrix}+o(\Delta
t)\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}^{T}\widehat{Q}\widehat{h}(\widetilde{x})+\frac{o(\Delta
t)^{2}}{\Delta
t}\widehat{h}(\widetilde{x})^{T}\widehat{Q}\widehat{h}(\widetilde{x})$
$\displaystyle\;+\Delta t^{2}\left(\begin{pmatrix}\widetilde{x}\\\ \lambda\\\
\lambda^{\prime}\end{pmatrix}+\frac{o(\Delta t)}{\Delta
t}\widehat{h}(\widetilde{x})\right)^{T}\widehat{Q}_{r}\left(\begin{pmatrix}\widetilde{x}\\\
\lambda\\\ \lambda^{\prime}\end{pmatrix}+\frac{o(\Delta t)}{\Delta
t}\widehat{h}(\widetilde{x})\right)$ (56)

where $\widetilde{x}=Ax+C\lambda$. The first equality follows from Lemma 8 for
all $\Delta t$ sufficiently small. Further by Lemma 9, we have that
$(\widetilde{x},\lambda,\lambda^{\prime})\in\widehat{{\cal C}}$.

Consider the only if part of the claim. Suppose LCS is exponentially stable
using an EQLF. Since $-\widehat{Q}\succ_{\widehat{\cal C}}0$ there exists
$\gamma>0$ such that
$-(\widetilde{x},\lambda,\lambda^{\prime})^{T}\widehat{Q}(\widetilde{x},\lambda,\lambda^{\prime})\geq\gamma\|(\widetilde{x},\lambda,\lambda^{\prime})\|^{2}$.
Then for all $\Delta t$ sufficiently small the terms other than the leading
term in (56) can be made smaller than or equal to $(\gamma\Delta
t)/2\|(\widetilde{x},\lambda,\lambda^{\prime})\|^{2}$. Then for all for such
$\Delta t$ we have from (56) and Lemma 9(c) that

$-(x,\lambda,\widehat{\lambda})^{T}\widehat{M}(x,\lambda,\widehat{\lambda})\geq(\gamma\Delta
t)/2\|(\widetilde{x},\lambda,\lambda^{\prime})\|^{2}\geq(\gamma\Delta
t)/(2c_{4})\|(x,\lambda,\widehat{\lambda})\|^{2}.$

The proof of the if part of the claim can be obtained in a similar manner.
This completes the proof. ∎

##  7 Conclusions & Future Work

In this work, we derived sufficient conditions for the Lyapunov stability of a
Discrete-Time Linear Complementarity System (DLCS), using both common and
extended quadratic Lyapunov functions. We also showed the equivalence between
the stability conditions of a Linear Complementarity System and its discrete-
time analog for all sufficiently small time-steps. We proposed a cutting plane
algorithm to find a Lyapunov function verifying exponential stability by
separating points from nonconvex copositive cones, and the algorithm was
demonstrated on small DLCS instances. The current investigation opens up a
number of avenues for future exploration.

  * •

First, the cutting plane algorithm can be slow on large instances of the DLCS.
For large instances, the algorithm often requires performing many rounds of
cuts before terminating. Improving the behavior of the cutting plane algorithm
is an important future direction of this work.

  * •

Another natural area of improvement for this approach is that a Lyapunov
function is only found upon termination of the algorithm. In other words, we
do not yet have a heuristic approach for constructing a feasible matrix
$P_{xx}$ using information obtained from the cutting plane method. Such an
approach can be advantageous since for stability, we are only interested in a
feasible solution satisfying the copositivity conditions.

  * •

The incorporation of feedback control and the joint computation of stabilizing
controller and Lyapunov function is another avenue for investigation. The
inclusion of feedback controls complicates the separation problem since the
matrices in the LCP constraints are now dependent on the feedback matrices.
This necessitates the development of a different algorithm for the computation
of Lyapunov function.

##  Appendix A S-Lemma Formulations for CQLF and EQLF

###  A.1 Common Quadratic Lyapunov Function

The conditions for exponential stability in Theorem 1$(iii)$ can be written as
the feasibility problem (12). Moreover, the inequality in (12c) can be relaxed
into a linear matrix inequality using the S-Lemma [3]. Consider the matrix
$H\in\mathbb{R}^{(n_{x}+n_{c})}$ defined as

$H=\begin{pmatrix}D&F\\\ 0&I_{n_{c}}\end{pmatrix}.$ (57a) Consider the Linear
Matrix Inequality (LMI) $\displaystyle M+H^{T}WH\prec 0\text{ where
}W\in\mathbb{S}^{(n_{x}+n_{c})},W\geq 0.$ (57b) Let
$\textnormal{conv}(\text{Gr}\,\text{SOL}(Dx,F))$ denote the convex hull of
$\text{Gr}\,\text{SOL}(Dx,F)$, i.e.
$\textnormal{conv}(\text{Gr}\,\text{SOL}(Dx,F))=\\{(x,\lambda)\,|\,\lambda\geq
0,\,Dx+F\lambda\geq 0\\}.$

We show that satisfaction of (57b) implies that $\psi(x,\lambda)<0$ for all
$(x,\lambda)\in\textnormal{conv}(\text{Gr}\,\text{SOL}(Dx,F))$. Thus, we
obtain a relaxation of the copositivity conditions in (12) to a LMI (57b).
Multiplying the matrix inequality in (57b) from the left and right by
$(x,\lambda)\in\textnormal{conv}(\text{Gr}\,\text{SOL}(Dx,F))$ we obtain

$\displaystyle\begin{pmatrix}x\\\
\lambda\end{pmatrix}^{T}\left(M+H^{T}WH\right)\begin{pmatrix}x\\\
\lambda\end{pmatrix}<0$ $\displaystyle\implies$
$\displaystyle\psi(x,\lambda)+\begin{pmatrix}Dx+F\lambda\\\
\lambda\end{pmatrix}^{T}W\begin{pmatrix}Dx+F\lambda\\\ \lambda\end{pmatrix}<0$
$\displaystyle\implies$ $\displaystyle\psi(x,\lambda)<0$ (57c)

where the second implication follows from
$(x,\lambda)\in\textnormal{conv}(\text{Gr}\,\text{SOL}(Dx,F))$ and
$W_{xx},W_{\lambda\lambda}\geq 0$. Hence, satisfaction of (57b) $\implies$
(12c). The computation of a $P_{xx}$ satisfying (12) can be cast as the LMI
system

Find $\displaystyle\;P_{xx}\in\mathbb{S}^{n_{x}},W\in\mathbb{S}^{n_{x}+n_{c}}$
(57d) such that $\displaystyle\;P_{xx}\succ 0,W\geq 0,\eqref{feasPxxlmi}.$
(57e)

###  A.2 Extended Quadratic Lyapunov Function

Similar to the use of S-Lemma for the CQLF the copositivity conditions in (18)
can be turned into a feasibility problem involving LMIs. Consider the matrices
$J,J^{\prime}\in\mathbb{R}^{(n_{x}+n_{c})\times(n_{x}+2n_{c})}$, defined as
$J=\begin{pmatrix}D&F&0\\\
0&I_{n_{c}}&0\end{pmatrix},\,J^{\prime}=\begin{pmatrix}DA&DC&F\\\
0&0&I_{n_{c}}\end{pmatrix}.$ (58a) It can be verified that
$\left.\begin{aligned} W_{1}\in\mathbb{S}^{n_{x}+n_{c}},W_{1}\geq 0\\\
P-H^{T}{W}_{1}H\succ 0\end{aligned}\right\\}\implies\eqref{feasP:setK1}$ (58b)
$\left.\begin{aligned}
W_{2},{W}_{3}\in\mathbb{S}^{(n_{x}+n_{c})},W_{2},W_{3}\geq 0\\\
\widehat{M}+J^{T}W_{2}J+(J^{\prime})^{T}W_{3}(J^{\prime})\prec
0\end{aligned}\right\\}\implies\eqref{feasP:setK2}.$ (58c)

Note that (58b) can be derived in a manner similar to that in (57c). We derive
(58c) in the following. Multiply the matrices in the left hand side of (58c)
by $(x,\lambda,\widehat{\lambda})\in\widehat{{\cal K}}$ and rearrange to
obtain

$\displaystyle\begin{pmatrix}x\\\ \lambda\\\
\widehat{\lambda}\end{pmatrix}^{T}\left(\widehat{M}+J^{T}W_{2}J+(J^{\prime})^{T}W_{3}(J^{\prime})\right)\begin{pmatrix}x\\\
\lambda\\\ \widehat{\lambda}\end{pmatrix}<0$ $\displaystyle\implies$
$\displaystyle\widehat{\psi}(x,\lambda,\widehat{\lambda})+\begin{pmatrix}Dx+F\lambda\\\
\lambda\end{pmatrix}^{T}W_{2}\begin{pmatrix}Dx+F\lambda\\\
\lambda\end{pmatrix}+$
$\displaystyle\begin{pmatrix}DAx+DC\lambda+F\widehat{\lambda}\\\
\widehat{\lambda}\end{pmatrix}^{T}W_{3}\begin{pmatrix}DAx+DC\lambda+F\widehat{\lambda}\\\
\widehat{\lambda}\end{pmatrix}<0$ $\displaystyle\implies$
$\displaystyle\widehat{\psi}(x,\lambda,\widehat{\lambda})<0$

where the second implication follows from $\lambda\in\text{SOL}(Dx,F))$,
$\widehat{\lambda}\in\text{SOL}(DAx+DC\lambda,F)$ and $W_{2},W_{3}\geq 0$.
This prove the implication in (58c). Just as in the case of the CQLF we can
show that the LMI in (58c) implies that
$\widehat{\psi}(x,\lambda,\widehat{\lambda})<0$ for all
$(x,\lambda,\widehat{\lambda})\in\textnormal{conv}{\widehat{{\cal K}}}$.

## References

  * [1] V. Acary and B. Brogliato.  Numerical Methods for Nonsmooth Dynamical Systems.  Lecture Notes in Applied and Computational Mechanics. Springer-Verlag Berlin Heidelberg, first edition, 2008.
  * [2] Alp Aydinoglu, Victor M Preciado, and Michael Posa.  Contact-Aware Controller Design for Complementarity Systems.  In International Conference on Robotics and Automation (ICRA), 2020\.
  * [3] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan.  Linear Matrix Inequalities in System and Control Theory.  SIAM, 1994.
  * [4] M. Branicky.  Multiple lyapunov functions and other analysis tools for switched and hybrid systems.  IEEE Trans. Automat. Control, 43:475–482, 1998.
  * [5] B. Brogliato.  Some perspectives on analysis and control of complementarity systems.  IEEE Trans. Automat. Control, 48:918–935, 2003.
  * [6] Stefan Bundfuss and Mirjam Dür.  An Adaptive Linear Approximation Algorithm for Copositive programs.  SIAM J Optimization, 20(1):30–53, 2009.
  * [7] Stefan Bundfuss and Mirjam Dür.  Copositive Lyapunov functions for switched systems over cones.  Systems and Control Letters, 58(5):342–345, 2009.
  * [8] M. K. Çamlibel, J-S. Pang, and J. Shen.  Lyapunov stability of complementarity and extended systems.  SIAM Journal on Optimization, 17(4):1056–1101, 2006.
  * [9] M. K. Çamlibel.  Complementarity Methods in the Analysis of Piecewise Linear Dynamical Systems.  PhD thesis, Center for Economic Research, Tilburg University, 2001.
  * [10] M. K. Çamlibel, W. P. M. H. Heemels, A. J. van der Schaft, and J. M. Schumacher.  Switched networks and complementarity.  IEEE Trans. Circuits Systems I: Fund. Theory Appl., 50:1036–1046, 2003.
  * [11] M. K. Çamlibel and J. M. Schumacher.  Open Problems in Mathematical Systems and Control Theory, chapter Copositive Lyapunov functions, page 189–193.  Princeton University Press, Princeton, NJ, 2004.
  * [12] Kenneth Chao and Pilwon Hur.  A step towards generating human-like walking gait via trajectory optimization through contact for a bipedal robot with one-sided springs on toes.  In IEEE International Conference on Intelligent Robots and Systems, pages 4848–4853, 2017.
  * [13] Xiaojun Chen and Zhengyu Wang.  Computational error bounds for a differential linear variational inequality.  IMA Journal of Numerical Analysis, 32(3):957–982, 2012.
  * [14] R.W. Cottle, J-S. Pang, and R.E. Stone.  The Linear Complementarity Problem.  SIAM, 1992.
  * [15] R. A. Decarlo, M. S. Branicky, S. Pettersson, and B. Lennartson.  Perspectives and results on the stability and stabilization of hybrid systems.  Proceedings of the IEEE, 88:1069–1082, 2000.
  * [16] Gabriele Eichfelder and Janez Povh.  On the set-semidefinite representation of nonconvex quadratic programs over arbitrary feasible sets.  Optimization Letters, 7(6):1373–1386, 2013.
  * [17] LLC Gurobi Optimization.  Gurobi optimizer reference manual, 2020.
  * [18] W. Heemels.  Linear complementarity systems: A Study in Hybrid Dynamics.  Phd thesis, Technische Universiteit Eindhoven, 2000.
  * [19] W. P.M.H. Heemels, J. M. Schumacher, and S. Weiland.  Linear complementarity systems.  SIAM Journal on Applied Mathematics, 60(4):1234–1269, 2000.
  * [20] J. P. Hespanha.  Uniform stability of switched linear systems: Extension of lasalle’s invariance principle.  IEEE Trans. Automat. Control, 49:470–482, 2004.
  * [21] J. P. Hespanha, D. Liberzon, D. Angeli, and E. D. Sontag.  Nonlinear norm-observability notions and stability of switched systems.  IEEE Trans. Automat. Control, 50:154–168, 2005.
  * [22] M. Johansson and A. Rantzer.  Computation of piecewise quadratic lyapunov functions for hybrid systems.  IEEE Trans. Automat. Control, 43:555–559, 1998.
  * [23] H.K. Khalil.  Nonlinear Systems.  Pearson Education. Prentice Hall, 2002.
  * [24] D. Liberzon.  Switching in Systems and Control.  Birkhauser Basel, first edition, 2003.
  * [25] Zachary Manchester and Scott Kuindersma.  Variational Contact-Implicit Trajectory Optimization Variational Contact-Implicit Trajectory Optimization.  International Symposium on Robotics Research (ISRR), (December):1–16, 2017.
  * [26] D.Q. Mayne, J.B. Rawlings, C.V. Rao, and P.O.M. Scokaert.  Constrained model predictive control: Stability and optimality.  Automatica, 36(6):789–814, jun 2000.
  * [27] I. Mordatch, E. Todorov, and Z. Popović.  Discovery of complex behaviors through contact-invariant optimization.  ACM Transactions on Graphics (TOG), 31(4 Article: 43), 2012.
  * [28] K.G. Murty and S.N. Kabadi.  Some np-complete problems in quadratic and nonlinear programming.  Mathematical Programming, 39:117–129, 1987.
  * [29] Anna Nagurney and Ding Zhang.  Projected Dynamical Systems and Variational Inequalities with Applications.  International Series in Operations Research & Management Science. Springer, first edition, 1996.
  * [30] Jong-Shi Pang and Jinglai Shen.  Strongly regular differential variational systems.  IEEE Transactions on Automatic Control, 52, 2007.
  * [31] Jong Shi Pang and David E. Stewart.  Differential variational inequalities.  Mathematical Programming, 113(2):345–424, 2008.
  * [32] A. Patel, S. L. Shield, S. Kazi, A. M. Johnson, and LT Biegler.  Contact-implicit trajectory optimization using orthogonal collocation.  IEEE Robotics and Automation Letters, 4(2):2242–2249, 2019.
  * [33] M. Posa, C. Cantu, and R. Tedrake.  A direct method for trajectory optimization of rigid bodies through contact.  The International Journal of Robotics Research (IJRR), 33(1):69–81, 2014.
  * [34] Michael Posa.  Optimization for Control and Planning of Multi-contact Dynamic Motion.  PhD thesis, Massachusetts Institute of Technology, 2017.
  * [35] A U Raghunathan, J. R. Pérez-Correa, E. Agosin, and L. T. Biegler.  Parameter estimation in metabolic flux balance models for batch fermentation—formulation & solution using differential variational inequalities (dvis).  Annals of Operations Research, 148(1):251–270, 2006.
  * [36] J. M. Schumacher.  Complementarity systems in optimization.  Math. Program. Ser. B, 101:263–296, 2004.
  * [37] Jinglai Shen and Jong Shi Pang.  Linear complementarity systems: Zeno states.  SIAM Journal on Control and Optimization, 44(3):1040–1066, 2006\.
  * [38] Robert Shorten, Fabian Wirth, Oliver Mason, Kai Wulff, and Christopher King.  Stability criteria for switched and hybrid systems.  SIAM Review, 49(4):545–592, 2007.
  * [39] A. J. van der Schaft and J. M. Schumacher.  Complementarity modeling of hybrid systems.  IEEE Trans. Automat. Control, 43:483–490, 1998.
  * [40] A. J. van der Schaft and J. M. Schumacher.  An Introduction to Hybrid Dynamical Systems.  Springer-Verlag, London, 2000.
  * [41] F. Vasca, L. Iannelli, M. K. Çamlibel, and R Frasca.  A new perspective for modeling powerelectronics converters: Complementarity framework.  IEEE Transactions on Power Electronics, 24(1-2):456–468, 2009.