











\section{Introduction}
\label{sec:introduction}
\IEEEPARstart{C}{omplementarity} conditions emerge naturally in a variety of contexts, from economics to engineering, from operations research to optimization.
Electrical networks with semiconductor devices~\cite{acary2011,adly2017,goeleven2017}, 
mechanical systems with unilateral constraints~\cite{brogliato1999}, price equilibrium, 
traffic network~\cite{nagurney1993}, and portfolio selection problems~\cite{ferris1997}
can all be modeled using a formalism in which decision variables must satisfy a set of 
complementarity conditions.
Linear complementarity problems (LCPs) are nonsmooth problems in which a set of complementarity conditions is coupled with a set of linear relations among the problem variables.
Due to the broad spectrum of applications, complementarity problems (and LCPs in 
particular) have been extensively studied, see, 
e.g.,~\cite{murty1988,isac1992,leenaerts1998,billups2000,facchinei2003,cottle}
and references therein. Complementarity problems appear in nonlinear control theory
in the form of linear complementarity systems (LCSs), which are defined by coupling a set of
complementarity conditions to a linear input--output control system, see, e.g.,~\cite{schaft1998,heemels2000,camlibel2002,brogliato2003}.

Because of their nonlinear nature, LCSs can undergo bifurcations at
which equilibria, limit cycles, or other types of steady-state solutions appear, 
disappear, or change stability. From a control-theoretical perspective, knowing 
which bifurcations can occur in a controlled system is instrumental for two
complementary reasons. Suppose the control objective is to preserve the stability
of a desired set point. In that case, bifurcations must be avoided, and knowing
the types of bifurcation that the system can undergo guides the design of
stabilizing control
laws~\cite{chen2003bifurcation,krener2004control,abed1986local,abed1987local,chen2000bifurcation}.
On the other hand, if the goal is to design rich nonlinear temporal and distributed
behaviors, then bifurcations must be sought because, in their vicinity, the nonlinear 
behavior can be modulated with enhanced sensitivity. This second use of bifurcations in
control systems has been explored
in~\cite{castanos2017implementing,franci2014realization,franci2019sensitivity} for the
design of neuro-inspired control systems and in~\cite{gray2018multiagent,franci2021analysis,bizyaeva2020general,franci2019model}
for the design of bio-inspired collective decision making in multi-agent systems.

As opposed to smooth bifurcation theory, a nonsmooth bifurcation theory for LCSs is entirely undeveloped.
Here, we elaborate on our preliminary results~\cite{castanos2020}
and introduce novel theoretical tools to analyze and design the possible ways in which the set of
solutions to an LCP changes as a function of control and system parameters. We propose these
tools as candidates for the development of a nonsmooth bifurcation theory for LCSs. 

The first contribution of our paper is to show that steady-state bifurcations (i.e., those bifurcations solely involving equilibrium points) of 
an LCS are fully characterized by the (nonsmooth)
singularities of an auxiliary LCP. The definition of the auxiliary LCP depends both
on the linear part and on the complementarity conditions of the LCS.
The problem of analyzing, classifying, and detecting steady-state bifurcation in LCSs
thus reduces to the problem of analyzing, classifying, and detecting nonsmooth singularities in LCPs.

The second main contribution of our paper is the introduction of a notion of equivalence between
LCPs. Various notions of equivalence between vector fields are at the basis of smooth bifurcation theory in its different
flavors~\cite{arnold,golubitsky1985,guckenheimer}. In line with the smooth theory, we declare two LCPs to be equivalent if one can be transformed into the other by suitable 
coordinate changes. Roughly speaking, according to our definition, two LCPs are equivalent if the geometry
of their solution space is the same, in the sense that regular (singular) solutions to the first LCP can be put in one-to-one correspondence
with regular (singular) solutions to the second LCP via a homeomorphism. Crucially, this notion of equivalence allows us to characterize
which LCPs are \emph{LCP-stable}, i.e., those LCPs for which the geometry of the solution space is not destroyed
by perturbations. This class of LCPs is particularly relevant in a control-theoretical setting because an
associated LCS inherits this robustness.

Again in line with the smooth theory, the coordinate transformations defining our notion of LCP-equivalence
are required to be continuous, i.e., they are topological isomorphisms. However,
topological equivalence alone turns out to be too crude a concept in the LCP setting. 
This is akin to the situation in smooth singularity theory where, for having a meaningful notion
of equivalence, additional structure (smoothness) is imposed on
the coordinate transformations.
The intrinsic nonsmoothness of LCPs places these problems at a different
stage, so a novel framework is required. It has long been recognized that the so-called
complementary cones associated with an LCP play a central role in their analysis and
reflect its rich geometric structure. Thus, the additional restriction 
we impose on the coordinate transformations is to preserve the Boolean
algebra generated by the complementary cones. The proposed equivalence
relation proves to be a powerful framework for the analysis, classification, and
design of LCPs and LCSs.

The applicability of our theory is illustrated in two extended examples.
In the first (more applied) example, we consider the LCS associated with a general 
electrical circuit with two transistors in feedback interconnection. It is known that
this circuit is essential to achieve negative resistance and thus bistability and
switches~\cite{chua1983}, but a rigorous and constructive characterization
of its steady states has remained elusive. Our theory allows us to derive constructive and rigorous
conditions to ensure LCP-stability, and thus robustness to disturbances and parametric uncertainties
of the circuit, and ensure bistability of its input-output behavior.
In the second (more theoretical) example, we provide an exhaustive classification 
of two-dimensional LCPs. Inspired by the concept of organizing center in smooth
singularity theory, we derive the classification by first characterizing
two-dimensional LCP-unstable complementarity problems and then characterizing
the open complement of this unstable set. 

The manuscript is organized in the following way: 
Section~\ref{sec:notation} introduces the notation 
used in this work, whereas linear complementarity problems and useful
properties are recalled in Section~\ref{sec:lcp}. Section \ref{sec:motivation} motivates
the proposed framework of LCP-equivalence in the context of LCSs. 
Afterward, Sections \ref{sec:equivalence} and \ref{sec:stability} constitute the main 
body of the paper, where the notions of LCP-equivalence and LCP-stability are studied, and 
the main results are presented. Subsequently, the proposed framework is applied in 
Section~\ref{sec:application} to the design of bistable LCSs and in Section~\ref{sec:classification} 
to the classification of LCPs.

\section{Notation} \label{sec:notation}

Let $X$ be a non-empty set. The power set of $X$ is denoted as $\mathcal{P}(X)$. 
Whenever $X$ is finite, $\vert X \vert$ denotes its cardinality. 
For a set $Y \subset X$, $Y^{c}$ denotes its complement in $X$, $Y^{c} := X \setminus Y$.
The notations $\interior Y$, $\boundary Y$ and $\closure{Y}$ stand for the interior,
the boundary and the closure of $Y$, respectively.

Given an integer $n$, we use $\Zint{n}$ to denote the interval 
of all integers between $1$ and $n$. We endow all subsets of $\mathbb{N}$
with the natural order. Let $A \in \RE^{n \times m}$, $\alpha \subseteq \Zint{n}$, and
$\beta \subseteq \Zint{m}$. The submatrix 
$A_{\alpha, \beta} \in \RE^{\vert \alpha \vert \times \vert \beta \vert}$ 
is the matrix $[A_{ij}]_{i \in \alpha,j \in \beta}$. The notation 
$A_{\bscdot, \beta}$ stands for the submatrix $A_{\Zint{n}, \beta}$.
We omit the braces whenever $\beta$ is a singleton, that is, 
$A_{\bscdot, j} = A_{\bscdot, \{j\}}$. The (polyhedral) cone generated by 
$A$ is the set of all positive linear combinations 
of the columns of $A$, that is,
\opt{twoside}{
 \begin{multline*}
 	\pos A := \Big\{ x \in \RE^{n} \mid x = \sum_{j = 1}^{m} p_{j} A_{\bscdot, j}, \\
 	p_{j} \geq 0, \, j \in \Zint{m} \Big\} \;.
 \end{multline*}
}
\opt{onecolumn}{
 \begin{displaymath}
 	\pos A := \Big\{ x \in \RE^{n} \mid x = \sum_{j = 1}^{m} p_{j} A_{\bscdot, j}, \\
 	p_{j} \geq 0, \, j \in \Zint{m} \Big\} \;.
 \end{displaymath}
}
The columns of $A$ are the \emph{generators} of $\pos A$.
The cone $\pos A$ is said to be \emph{pointed} if there is no non-trivial 
subspace contained in it, whereas it is said to be \emph{strictly pointed} 
if $\ker(A) \cap \RE_{+}^{m} = \{ 0 \}$, where 
\begin{displaymath}
 \RE_{+}^{m} := \{ x\in \RE^{m} \mid x_{i} \geq 0, \, i \in \Zint{m} \}
\end{displaymath}
denotes the positive orthant of $\RE^{m}$ and $\ker(A)$
is the null space of $A$. 
Finally, we denote the projection of $x \in \RE^{n}$ onto the positive orthant $\RE_{+}^{n}$
by $[x]^{+}$, that is,
\begin{displaymath}
	[x]^{+} := \argmin_{w \in \RE_{+}^{n}} \frac{1}{2} \Vert x - w \Vert^{2} \;.
\end{displaymath}

\section{The linear complementarity problem}
\label{sec:lcp}
We start with a brief introduction to the LCP and its geometric properties. We also recall some 
equivalent formulations that will be useful in the sequel.

\begin{defn}
	\label{defn:lcp}
	Given a matrix $M \in \RE^{n \times n}$ and a vector $q \in \RE^{n}$, the 
   \emph{linear complementarity problem} $\lcp(M,q)$ consists in finding a vector 
	$z \in \RE^{n}$ such that

	\begin{equation} \label{eq:complementarity:condition}
		\RE_{+}^{n} \ni M z + q \perp z \in \RE_{+}^{n} \;.
	\end{equation}

	Condition~\eqref{eq:complementarity:condition} is called the 
	complementarity condition, it is the compact form of the following three 
	conditions: $1)$ $M z + q \in \RE_{+}^{n}$; $2) z \in \RE_{+}^{n}$; and 
$3)$ $(Mz + q)^{\top} z = 0$.
\end{defn}
 
Linear complementarity problems have a rich geometric structure that determines the 
properties of their solution set. In what follows, we introduce the necessary elements
for the geometric analysis of LCPs.
 
\begin{defn}
	\label{defn:complementary:cone}
	Given an $\lcp(M,q)$ and a subset $\alpha \subseteq \Zint{n}$. 
   The \emph{complementary matrix} $C_{M}(\alpha) \in \RE^{n \times n}$ is the matrix
   whose columns are

	\begin{equation}
		\label{eq:complementary:matrix}
		C_{M}(\alpha)_{\bscdot, j}
		=
		\begin{cases}
			-M_{\bscdot, j} & j \in \alpha
			\\
			I_{\bscdot, j} & j \notin \alpha
		\end{cases} \;.
	\end{equation}
	The cone generated by $C_{M}(\alpha)$, $\pos C_{M}(\alpha)$, is called a 
	\emph{complementary cone} of $M$.
\end{defn}
 
A \emph{facet} of $\pos C_M(\alpha)$ is an $(n-1)$-dimensional face of the form
$\pos C_M(\alpha)_{\bscdot, i^c}$ for some $i \in \Zint{n}$ and $i^c$ the complement
of $\{i\}$ in $\Zint{n}$.
In what follows,  $\mathcal{K}(M)$ denotes the union of all facets of the complementary cones
of $M \in \RE^{n \times n}$.

The complementary cone $\pos C_{M}(\alpha)$ is called \emph{non-degenerate} if
it is solid, that is, if $C_{M}(\alpha)$ is nonsingular. 
Otherwise, $\pos C_{M}(\alpha)$ is called \emph{degenerate}. 

\begin{defn}
 A matrix $M \in \RE^{n\times n}$ is said to be \emph{non-degenerate} if all 
 its complementary cones, $\pos C_M(\alpha)$, $\alpha \subseteq \Zint{n}$ 
 are non-degenerate.
\end{defn}

\begin{defn}[\cite{danao1994,doverspike1982}]
	\label{defn:r0}
   We say that $M \in \RE^{n \times n}$ \emph{belongs to the class $R_{0}$} if,
   for each $\alpha \subseteq \Zint{n}$, we have $\ker(C_{M}(\alpha)) \cap 
   \RE_{+}^{n} = \{0\}$.
\end{defn}

The previous definition implies that all the complementary cones are
strictly pointed, that is, that $M$ does not have zero columns nor complementary cones that
are subspaces~\cite[Thm. 6.1.19]{cottle}. 

Now we define the piecewise linear map $f_M: \RE^{n} \to \RE^{n}$ as
\begin{equation}
	\label{eq:f:piecewise:linear}
	f_{M}(x) := C_{-M}(\alpha) x \; \text{for $x \in \pos C_{I}(\alpha)$} \;.
\end{equation}
Note that the cones $\pos C_I(\alpha)$ are simply the $2^n$ orthants in $\RE^n$ 
indexed by $\alpha \subseteq \Zint{n}$, and that 
\begin{equation} \label{eq:f:cones}
 f_M(\pos C_I(\alpha)) = \pos C_M(\alpha) \;.
\end{equation}

\begin{thm}[\cite{cottle}]
	\label{thm:lcp:z:x:equivalent}
	Let $z \in \RE_{+}^{n}$ be a solution to the $\lcp(M,q)$. 
	Then, $x = (M - I_{n})z + q$ is a solution to the equation

	\begin{equation}
		\label{eq:lcp:f}
		f_{M}(x) - q = 0
	\end{equation}

	Conversely, if $x \in \RE^{n}$ is a solution to \eqref{eq:lcp:f}, then 
	$z = [-x]^{+}$ is a solution to the $\lcp(M,q)$.
\end{thm}

It follows that the $\lcp(M,q)$ and the problem~\eqref{eq:lcp:f} are equivalent in the 
sense that we only need to know the solutions to one of them in order to know
the solutions to the other. 

The functional representation~\eqref{eq:lcp:f} lends itself to the application of degree theory.

\begin{defn}
	Let $M \in \RE^{n \times n}$ be an $R_{0}$-matrix and $x \in \RE^{n}$ be in 
	the interior of $\pos C_{I}(\alpha)$ for some $\alpha \subseteq \Zint{n}$.
	The \emph{index} of $f_{M}$ at $x$, denoted as $\ind_{M}(x)$, is given by

	\begin{displaymath}
		\ind_{M}(x) = \sgn (\det M_{\alpha, \alpha}) \;.
	\end{displaymath}

	By convention, if $\alpha = \emptyset$, then $\ind_{M}(x) = 1$.
\end{defn}

\begin{defn}
	Let $M \in \RE^{n \times n}$ be an $R_{0}$-matrix and let $q \in \RE^{n}$ be such
	that $f_{M}^{-1}(q)$ consists of finitely many points and $\ind_{M}(x)$ is 
	well-defined for all $x \in f_{M}^{-1}(q)$. The \emph{degree} of $f_{M}$ at $q$, 
	denoted as $\deg_{M}(q)$ is

	\begin{displaymath}
		\deg_{M}(q) = \sum_{x \in f_{M}^{-1}(q)} \ind_{M}(x) \;.
	\end{displaymath}

\end{defn}

Note that, for a given $\alpha \subseteq \Zint{n}$, the index of $f_{M}$ at $x$ is the same 
for all $x \in \interior \pos C_{I}(\alpha)$ so, there are as many indices as there are 
subsets in $\Zint{n}$. For matrices in the class $R_{0}$, the degree is the same for any 
$q \notin \mathcal{K}(M)$~\cite[Theorem 6.1.14]{cottle}. Thus, for $R_{0}$-matrices the degree 
is global and one only needs to test one point $q \notin \mathcal{K}(M)$. The reader is addressed
to~\cite{cottle,garcia1981,howe1981} and references therein for a detailed 
account of degree theory in the context of complementarity problems.

The solutions to the $\lcp(M,q)$ depend on the geometry of the complementary cones 
$\pos C_{M}(\alpha)$. More precisely, there exists at least one solution $x$ 
of~\eqref{eq:f:piecewise:linear} for every $\alpha \subseteq \Zint{n}$ such
that $q \in \pos C_{M}(\alpha)$. If $C_{M}(\alpha)$ is nonsingular the associated solution 
is unique, whereas there exists a continuum of solutions if $C_{M}(\alpha)$ is singular. Thus, 
for a given $q$, there can be no solutions, there can be one solution, multiple isolated solutions, 
or a continuum of solutions, depending on how many complementary cones $q$ belongs to 
and what the properties of these cones are.

We will later make use of an alternative representation of $f_{M}$.
\begin{lem}
	\label{lemma:pwl:f:projection}
	The map~\eqref{eq:f:piecewise:linear} can also be written as

	\begin{equation}
		\label{eq:pwl:f:projection}
		f_{M}(x) = [x]^{+} - M [-x]^{+} \;.
	\end{equation}
\end{lem}
\begin{proof}
	It follows from \eqref{eq:complementary:matrix} that, for $x \in \pos C_{I}(\alpha)$,
	we have $x_{j} \leq 0$ if $j \in \alpha$, whereas $x_{j} \geq 0$ if $j \notin \alpha$. 
   Then, equation~\eqref{eq:f:piecewise:linear} can be written as

	\begin{align*}
		f_{M}(x) &= \sum_{j = 1}^{n} C_{-M}(\alpha)_{\bscdot, j} x_{j}
		 = \sum_{j \notin \alpha} I_{\bscdot, j} x_{j} + \sum_{j \in \alpha} 
		 M_{\bscdot, j} x_{j} \\
		&= \sum_{j = 1}^{n} I_{\bscdot, j} [x_{j}]^{+} -
		 \sum_{j=1}^{n} M_{\bscdot, j} [-x_{j}]^{+} \;,
	\end{align*}
  
   from which~\eqref{eq:pwl:f:projection} follows directly.
\end{proof}

\section{Bifurcations in linear complementarity systems}
\label{sec:motivation}

We are interested in LCSs of the form
\begin{subequations}\label{eq:Sigma}
\begin{align}
 \dot{\xi}(t) &= A \xi(t) + B z(t) + E_{1} r \label{eq:lcs:1} \\
         w(t) &= C \xi(t) + D z(t) + E_{2} s \label{eq:lcs:2} \;, \\
 & \RE_{+}^{m} \ni w(t) \perp z(t) \in \RE_{+}^{m} \label{eq:lcs:3} 
\end{align}
\end{subequations}
where $\xi(t) \in \RE^{n}$ is the state of the system at time $t$,
$z(t)$, $w(t) \in \RE^{m}$ are the so-called complementary variables (they can 
be interpreted as external port variables), and $r \in \RE^{l}$, $s \in \RE^{l}$ are 
vectors of parameters. The matrices $A, B, C, D, E_{1}$ and $E_2$ are constant and of
appropriate dimensions. Allow us to provide an example. 

\subsection{LCS model of a circuit with bipolar transistors} \label{sec:circuit_model}

Motivated by the work of Chua \emph{et al.} on negative resistance devices~\cite{chua1983},
we consider the circuit in Fig.~\ref{fig:circuit:nr}. The circuit topology is fairly general
as, by opening or short-circuiting some of the resistors, it is possible to recover many of the
circuits considered in~\cite{chua1983}.

\begin{figure}[t]
	\centering
	\input{shockley_NegativeResistance.tex}
	\caption{Feedback interconnection of PNP and NPN 
		transistors showing differential
	negative resistance between terminals $P_{1}$ and $P_{0}$.}
	\label{fig:circuit:nr}
\end{figure}


The circuit contains two bipolar junction transistors (BJTs) that we describe using
a dynamical version of the Ebers-Moll model depicted in Fig.~\ref{fig:ebersMoll} 
(see, e.g., \cite{getreu1978}), where each capacitor models the charge 
storage effects of the P-N and N-P junctions. 
\begin{figure}[t]
	\centering
	\input{ebersMoll_model_pnp_dissipation.tex}
	\input{ebersMoll_model_npn_dissipation.tex}
	\caption{Ebers-Moll model of bipolar junction transistors. 
		PNP transistor (above), and NPN	transistor (below).}
	\label{fig:ebersMoll}
\end{figure}


Each diode in Fig.~\ref{fig:ebersMoll}
satisfies the complementarity condition
\begin{displaymath}
	0 \leq I_{d} \perp V^{*}-V_{d} \geq 0 \;,
\end{displaymath}
where $I_{d}$ is the current flowing through the diode, $V_{d}$ is the voltage across 
its terminals, and $V^{*}$ is the forward voltage of the diode (typically between $0.3$ V--$0.7$ V).
After the application of Kirchhoff's laws, we realize that the circuit
evolves according to~\eqref{eq:lcs:1}-\eqref{eq:lcs:3} with state $\xi = [x_{1}^{a}, x_{2}^{a}, 
x_{1}^{b}, x_{2}^{b}]^{\top} \in \RE^{4}$. Here, $x_{k}^{a}$ and $x_{k}^{b}$, 
$k \in \{1, 2\}$, are the voltages across the capacitors of the transistors 
$Q_{a}$ and $Q_{b}$, respectively (see Fig.~\ref{fig:ebersMoll}). The 
complementary variables are given by
$z = [I_{EC}^{a}, I_{CC}^{a}, I_{EC}^{b}, I_{CC}^{b}] \in \RE^{4}$
and $w = [V^{*} - V_{dR}^{a}, V^{*}-V_{dF}^{a}, V^{*} - V_{dR}^{b}, V^{*} - V_{dF}^{b}]^{\top} \in \RE^{4}$.
We have the parameters $s = V^{*}$ and $r = P_{1} - P_{0}$, the latter being the potential difference
across the terminals shown in Fig.~\ref{fig:circuit:nr}. The remaining parameters
are
\begin{equation} \label{eq:lcs:param:1}
\begin{aligned}
 A &= Q
  \begin{bmatrix}
  	A_{1,1} & A_{1,2}
  	\\
	A_{1, 2}^{\top} & A_{2,2}
 \end{bmatrix}, \quad B =  Q
 \begin{bmatrix}
	T & 0
	\\
	0 & T
 \end{bmatrix}, \\
 E_{1} &=  Q [-G_{2}, G_{1} + G_{2}, -G_{1}, G_{1} + G_{2}]^{\top}, \quad E_{2} = \mathbf{1}_{4}, \\
 Q &= \diag\left( \frac{1}{C_{1}^{a}}, \frac{1}{C_{2}^{a}}, \frac{1}{C_{1}^{b}}, 
   \frac{1}{C_{2}^{b}}\right), \quad C  = -I_{4}, \\
 D &= \diag \left( \frac{R_{R}}{\alpha_{R}}, \frac{R_{F}}{\alpha_{F}}, 
   \frac{R_{R}}{\alpha_{R}}, \frac{R_{F}}{\alpha_{F}} \right) ,
 \end{aligned}
\end{equation}
where $\mathbf{1}_{4} \in \RE^4$ is a vector of ones and
\begin{equation} \label{eq:lcs:param:2}
\begin{aligned}
	A_{1,1} & = 
	\begin{bmatrix}
		- (G_{0}^{a} + G_{2}^{a} + G_{2}) 
		& G_{2}^{a} + G_{2}
		\\
		G_{2}^{a} + G_{2} 
		& -(G_{1}^{a} + G_{2}^{a} + G_{1} + G_{2})
	\end{bmatrix}, \\
	A_{2, 2} & = 
	\begin{bmatrix}
		-(G_{0}^{b} + G_{2}^{b} + G_{1}) 
		& G_{2}^{b} + G_{1}
		\\
		G_{2}^{b} + G_{1} 
		& -(G_{1}^{b} + G_{2}^{b} + G_{1} + G_{2})
	\end{bmatrix}, \\
	A_{1,2} & = 
	\begin{bmatrix}
		0 & G_{2}
		\\
		G_{1} & - (G_{1} + G_{2})
	\end{bmatrix}, \quad 
	T = 
	\begin{bmatrix}
		-\frac{1}{\alpha_{R}} & 1
		\\
		1 & -\frac{1}{\alpha_{F}}
	\end{bmatrix},
\end{aligned}
\end{equation}
with $G_{k}^{j} = 1/R_{k}^{j}$, $k \in \{1,2\}$, $j \in \{a, b\}$. Finally, 
$\alpha_{F}$, $\alpha_{R} \in (0,1)$ are the forward and reverse gains of the bipolar transistors.

\subsection{Steady-state LCS bifurcations are determined by LCP singularities}

The main issues addressed in the literature about LCS are well-posedness and
asymptotic behavior of the resulting trajectories (see, 
e.g.,~\cite{brogliato2003,brogliato2020,camlibel2002,heemels2000,schaft1998}).
If $D$ is a $P$-matrix (i.e., all its principal minors are positive) then, for every $\xi$ and $s$, 
there exists a unique $z$ satisfying \eqref{eq:lcs:2},\eqref{eq:lcs:3} and given by
\begin{equation} \label{eq:z}
	z = \left[ -f_D^{-1}(C\xi + E_2 s) \right]^+ \;.
\end{equation}
Substituting~\eqref{eq:z} in~\eqref{eq:lcs:1},~\eqref{eq:Sigma} can be written explicitly as
\begin{equation} \label{eq:explicit}
	\dot{\xi}(t) = F(\xi(t),r,s)
\end{equation}
with
\begin{equation}\label{eq:F def}
	F(\xi,r,s) = A \xi + B \left[ -f_D^{-1}(C\xi + E_2 s) \right]^+ + E_{1} r \;.
\end{equation}
Furthermore, it can be readily shown that $F$ is globally Lipschitz-continuous in $\xi$,
so that~\eqref{eq:lcs:1}-\eqref{eq:lcs:3} is well posed, i.e., its solutions
are global and unique. A less-explored problem is understanding how the equilibria of~\eqref{eq:Sigma} change, appear, and 
disappear as a function of the input parameters $r$ and $s$, that is, what is the nonsmooth steady-state bifurcation 
structure of~\eqref{eq:Sigma}.

The following definition formalizes the notion of non-smooth (steady-state) bifurcation. 
It uses a notion of regularity that is closely related to that proposed in 
\cite{robinson1980} and it is inspired by the notion of smooth stead-state bifurcation 
used in~\cite{golubitsky}.

\begin{defn}
 Given a Lipschitz-continuous map depending on $m$ parameters, $g:\RE^n \times \RE^m\to \RE^n$, 
 a point $(y_0,p_0)$ is a \emph{regular solution} of $g$ if $g(y_0,p_0)=0$ and there exists a
 neighborhood $W=U \times V$ of $(y_0,p_0)$ such that the map $p \mapsto \{y \in U \mid g(y,p)=0 \}$ 
 is single-valued and Lipschitz-continuous over $V$\footnote{To simplify the exposition, a
 multi-valued map whose image is a singleton will be treated as a single-valued map.}. Solutions 
 of $g(y,p)=0$ that are not regular are called \emph{bifurcations} of $g$.
\end{defn}

Let $\partial_y g(y,p)$ denote the Clarke generalized Jacobian of $g$ with respect
to $y$ at $(y,p)$, that is,
\begin{displaymath}
	\partial_y g(y,p) = \co \Big\{ \lim_{i \to \infty} D_y g(y_i,p) \mid y_i \to y,
	y_i \not\in \Omega_g \Big\} \;,
\end{displaymath}
where $\Omega_g$ is the set of measure zero where the 
Jacobian $D_y g$ of the Lipschitz-continuous map $g$ does not exist. We say that $\partial_y g(y,p)$ 
is of \emph{maximal rank}, or \emph{nonsingular}, if every matrix in $\partial_y g(y,p)$ is nonsingular. 
Conversely, we say that $\partial_y g(y,p)$ is singular if it contains a singular matrix. The following
proposition is a direct consequence of the nonsmooth Implicit Function Theorem~\cite[p. 256]{clarke1990}.
\begin{prop}\label{prop: bifu to singu}
	If $(y_0,p_0)$ is a bifurcation of $g$ then $\partial_y g(y_0,p_0)$ is singular.
\end{prop}

The following theorem shows that equilibria of~\eqref{eq:Sigma} can be put in one-to-one 
correspondence with the solution of an associated LCP.

\begin{thm} \label{thm:equilibria}
 Consider~\eqref{eq:Sigma} with $A, B, C$ nonsingular and $D$ a $P$-matrix. Define
 \begin{align}
             M &= D - C A^{-1} B \label{eq:lcp:M} \\
  \bar{q}(r,s) &= E_{2}s - C A^{-1} E_{1} r \label{eq:lcp:q} \\
	   X(\xi,s) &= f_D^{-1}(C\xi + E_2 s) \label{eq:xi to x}
 \end{align}
 and denote by
 \begin{equation} \label{eq:x to xi}
  X^{-1}(x,s)=C^{-1}(f_D(x)-E_2 s)
 \end{equation}
 the inverse of $X(\cdot,s)$. Then, $F(\xi,r,s)=0$ implies that $f_M(X(\xi,s)) = \bar{q}(r,s)$. 
 Conversely, $f_M(x)=\bar q(r,s)$ implies that $F(X^{-1}(x,s),r,s) = 0$.
\end{thm}

Note that, since $D$ is a $P$-matrix, $X$ and $X^{-1}$ are well-defined, see e.g., \cite[Theorem 3.3.7]{cottle}.

\begin{proof}
	It follows from Lemma \ref{lemma:pwl:f:projection} that $f_{D}(x) = [x]^{+} - D [-x]^{+}$ and 
   $f_{M}(x) = [x]^{+} - (D - CA^{-1}B)[-x]^{+}$, leading to
	\begin{equation} \label{eq:fM:fD}
		f_{M}(x) = f_{D}(x) + C A^{-1}B[-x]^{+} \;.
	\end{equation}
	Suppose that $F(\xi, r, s) = 0$ and let $x=X(\xi,s)$. By~\eqref{eq:xi to x}, we have
   \begin{displaymath}
	 f_{D}(x) = C \xi + E_{2} s
   \end{displaymath}
   which, when substituted in \eqref{eq:fM:fD}, yields
	\begin{equation} \label{eq:fM:1}
	 f_{M}(x) = C \xi + E_{2} s + C A^{-1} B [-x]^{+} \;.
	\end{equation}
   Since $A$ is nonsingular, we can solve for $\xi$ in~\eqref{eq:F def} to obtain
	\begin{equation} \label{eq:xi:0}
	 \xi = -A^{-1}(B[-x]^{+} + E_{1}r) \;.
	\end{equation}	
	Substitution of \eqref{eq:xi:0} into \eqref{eq:fM:1} leads us to $f_{M}(x) = \bar q(r,s)$.
	
	Conversely, suppose that $f_{M}(x) = \bar{q}(r,s)$ and set $\xi = X^{-1}(x,s)$.  
   It follows from \eqref{eq:fM:fD} and~\eqref{eq:x to xi} that

	\begin{equation} \label{eq:kernel:C}
		\xi + A^{-1}\left( B[-x]^{+} + E_{1} r \right) \in \kernel C \;,
	\end{equation}

	where $\kernel C$ denotes the null space of $C$. Since $C$ is nonsingular by hypothesis, then \eqref{eq:xi:0} holds and therefore $F(\xi, r, s) = 0$.
\end{proof}

\begin{rem}
	In the more general case where $B \in \RE^{n \times m}$, $C \in \RE^{m \times n}$, and 
	$A \in \RE^{n \times n}$ is nonsingular, all equilibria of~\eqref{eq:Sigma} satisfy 
	$f_{M}(x) = \bar{q}(r,s)$ for $x = X(\xi, s)$. However,
   not all solutions of $f_{M}(x) = \bar{q}(r,s)$
   can be mapped to equilibrium points of~\eqref{eq:Sigma}. 
\end{rem}

\begin{cor}
 Under the same assumptions and definitions used in Theorem~\ref{thm:equilibria}, 
 the point $(\xi_0,r_0,s_0)$ is a bifurcation of $F$ if, and only if,
 $(x_0,q_0)=(X(\xi_0,s_0),\bar{q}(r_0,s_0))$ is a bifurcation of 
 $(x,q) \mapsto f_M(x)-q$.
\end{cor}

\begin{proof}
By Theorem~\ref{thm:equilibria},
 \begin{displaymath}
  \{x \mid  f_M(x)=\bar q(r,s) \} = \{ X(\xi,s) \mid F(\xi,r,s)=0 \}
 \end{displaymath}
and, similarly,
 \begin{displaymath}
 \{ \xi \mid F(\xi,r,s)=0 \} =	\{ X^{-1}(x,s) \mid  f_M(x)=\bar q(r,s) \} \;.
\end{displaymath}
Then, the map $(r,s)\mapsto\{x \mid  f_M(x)=\bar q(r,s)\}$ can be expressed as the
composition of the map $(r,s)\mapsto\{ \xi \mid F(\xi,r,s)=0 \}$ and the Lipschitz map
$X(\cdot,s)$, and, vice versa,  the map $(r,s)\mapsto\{ \xi \mid F(\xi,r,s)=0 \}$ can be expressed
as the composition of the map $(r,s)\mapsto\{x \mid  f_M(x)=\bar q(r,s)\}$ and the Lipschitz map
$X^{-1}(\cdot,s)$. Recalling that the composition of Lipschitz-continuous maps is Lipschitz-continuous,
the map $(r,s)\mapsto\{ \xi \mid F(\xi,r,s)=0\}$ is single-valued and Lipschitz-continuous 
in a neighborhood of $(\xi_0,r_0,s_0)$ if, and only if, the map 
$(r,s)\mapsto\{x \mid  f_M(x)=\bar q(r,s)\}$ is single-valued and Lipschitz-continuous in a neighborhood 
$(X(\xi_0,s_0),r_0,s_0)$.
\end{proof}

It follows that the steady-state bifurcations of an LCS are fully characterized by the bifurcations of an associated LCP. Furthermore, by Proposition~\ref{prop: bifu to singu}, we can make two key observations. First, LCS bifurcations are non-hyperbolic equilibria of~\eqref{eq:Sigma}, in the sense that Jacobian $\partial_{\xi} F(\xi,r,s)$ at a bifurcation point contains matrices with some zero eigenvalues. As such, system~\eqref{eq:Sigma} can lose stability at bifurcations. Second, bifurcations of an LCP can be found constructively by looking for its {\em singularities}, i.e., points $(x_0,q_0)$ such that $f_M(x_0)-q_0=0$ and $\partial_{x} f_M(x_0)$ is singular. Observing that 
\begin{displaymath}
	\partial f_M(x)=\co\left\{ C_{-M}(\alpha)\mid x\in\pos C_I(\alpha) \right\}  \;,
\end{displaymath}
it follows that $\partial f_M(x)$ is a singleton whenever $x$ belongs to the interior 
of an orthant and a set of matrices whenever $x$ belongs to the common boundary of 
two or more orthants. 
Hence, for a point $(x_{0}, q_{0})$ to be a nonsmooth singularity of the $\lcp(M,q)$ it 
is necessary that, for some $\alpha \subseteq \Zint{n}$, either 
$x_{0} \in \boundary \pos C_{I}(\alpha)$ 
or $x_{0} \in \interior \pos C_{I}(\alpha)$ with $C_{-M}(\alpha)$ singular. 
In the former case $\partial f_{M}$ is a set for which the maximal rank condition may 
not hold, whereas in the latter case $\partial f_{M}$ is a singular matrix and therefore
it is not of maximal rank.
Note that, in both cases, $q_{0} \in \boundary C_{M}(\alpha)$.
Therefore, nonsmooth bifurcations in LCPs are essentially determined by the 
configuration of the complementarity cones and how $q$ moves across them. 
The notion of equivalence for LCPs introduced in the next section helps in the task of analyzing
and classifying complementarity-cone configurations, as illustrated through both applied and theoretical extended examples.

\section{Equivalence of LCPs} \label{sec:equivalence}

Our notion of equivalence between an $\lcp(M,q)$ and an $\lcp(N,r)$ involves both topological
and algebraic properties. The algebraic properties deal with the relations among the 
complementary cones of $M$ and $N$. The relevant algebraic structure is that of a 
Boolean algebra, a subject that we now briefly recall (see~\cite{givant,sikorski} 
for more details).

\begin{defn}
	\label{defn:field:sets}
	Let $X$ be a set. A \emph{field of sets} is a pair $(X, \Fi)$, where 
   $\Fi \subseteq \mathcal{P}(X)$ is any non-empty family of sets that 
   is closed under the set operations of complement and finite union.
\end{defn}

It follows easily from the De Morgan rules that a field of sets $(X, \Fi)$ 
is also closed under finite intersections of sets. Note also that $\emptyset$ and $X$ 
are always members of $\Fi$. In the cases when the set $X$ is clear from the
context, we will denote the field of sets simply by the collection $\Fi$.
Fields of sets are concrete examples of Boolean algebras, see, e.g.,~\cite{sikorski}.
Therefore, the usual algebraic concepts apply to them.

\begin{defn}
	\label{defn:boolean:homeomorphism}
	Consider the fields $(X, \Fi)$ and $(\hat{X}, \hat{\Fi})$. 
   A mapping $h: \Fi \to \hat{\Fi}$ is said to be \emph{Boolean} 
   if, for all $P_{1}, P_{2} \in \Fi$,

	\begin{displaymath}
		h(P_{1} \bigcap P_{2}) = h(P_{1}) \bigcap h(P_{2})
	\end{displaymath}
	and
	\begin{displaymath}
		h(X \setminus P_{1}) = \hat{X} \setminus h(P_{1}) \;.
	\end{displaymath}
   A Boolean map that is also a bijection is said to be a \emph{Boolean isomorphism}. A
   Boolean isomorphism that maps a field of sets to itself is said to be a \emph{Boolean
   automorphism}.
\end{defn}

\begin{defn}
	\label{defn:boolean:map:induced}
   A map $\varphi_*: \Fi \to \hat{\Fi}$ is said to be
   \emph{induced} by a map $\varphi: \hat{X} \to X$ if, for every set 
   $P \in \Fi$,

	\begin{displaymath}
		\varphi_*(P) = \varphi^{-1}(P) \;.
	\end{displaymath}
\end{defn}

Let $\Ge \subseteq \mathcal{P}(X)$. The field of sets \emph{generated} by
$\Ge$, denoted $\Span \Ge$, is the intersection of all the fields of sets containing $\Ge$.
That is, it is the smallest field of sets containing $\Ge$.
In what follows, we consider the collection of complementary cones $\Ge_{M} = 
\{\pos C_{M}(\alpha)\}_{\alpha}$ and we denote the field of sets generated by 
$\Ge_{M}$ as $\Fi_{M}$.

We are now ready to introduce our main definition.

\begin{defn}
	\label{defn:lcp:equivalence}
	Two matrices $M, N \in \RE^{n \times n}$ are \emph{LCP-equivalent}, denoted
	$M \sim N$, if the fields of sets $\Fi_{M}$ and $\Fi_{N}$ are isomorphic
	and the Boolean isomorphism $\varphi_*: \Fi_{M} \to \Fi_{N}$ is induced
	by a homeomorphism $\varphi: \RE^{n} \to \RE^{n}$.
\end{defn}

In other words, two matrices are LCP-equivalent if the Boolean structures of their
associated complementary cones are isomorphic. That is, if there is a bijection
$\varphi_*: \Fi_{M} \to \Fi_{N}$ preserving intersections, unions, 
and complements of complementary cones. Since $\varphi_*$ is required to be induced by a 
homeomorphism, it suffices to verify that $\varphi$ induces a bijection from
$\Ge_{M}$ to $\Ge_{N}$~\cite[Cor. 12]{castanos2020}.

As an illustration, let us consider the matrices
\begin{equation}
	\label{eq:example:2:M:N}
	M = 
	\begin{bmatrix}
		1 & 1
		\\
		1 & 1
	\end{bmatrix} \quad \text{and} \quad N =
	\begin{bmatrix}
		1 & -1
		\\
		1 & 0
	\end{bmatrix} \;.
\end{equation}
The associated complementary cones are depicted in Fig.~\ref{fig:example:2:cones:M:N}.
Note that, although $M$ is singular and $N$ is not, they are in fact LCP-equivalent.
To see this, consider the homeomorphism $\varphi: \RE^{2} \to \RE^{2}$,
\begin{displaymath}
	\varphi(q') = \begin{bmatrix}
		-1 & 0
		\\
		-1 & 1
	\end{bmatrix} q' \;.
\end{displaymath}
Simple but lengthy computations show that 
\begin{displaymath}
 \varphi_* : \pos C_{M}(\alpha) \mapsto \pos C_{N}(\beta(\alpha))
\end{displaymath}
with $\beta: \mathcal{P}(\Zint{n}) \to \mathcal{P}(\Zint{n})$ given by
\begin{displaymath}
	\emptyset \mapsto \{1\} \;,  \quad 
	   \{1\}  \mapsto \emptyset \;, \quad
	  \{2 \}  \mapsto \{1, 2\} \;, \quad
	\{1, 2\}  \mapsto \{2\} \;.
\end{displaymath}
That is, $\varphi$ induces a bijection $\Ge_{M} \to \Ge_{N}$ and hence
an isomorphism $\Fi_M \to \Fi_N$.

\begin{figure}
	\centering
	\input{example_2_M_N.tex}
	\caption{Complementary cones associated with the matrices $M$ (left) and $N$ (right) 
	from \eqref{eq:example:2:M:N}. Red rays denote degenerate complementary
cones.}
	\label{fig:example:2:cones:M:N}
\end{figure}

\begin{thm}
	\label{thm:equivalence:suf}
	Consider two matrices $M, N \in \RE^{n \times n}$ and suppose there exists
 	homeomorphisms $\varphi, \psi: \RE^{n} \to \RE^{n}$ such that
	\begin{equation}
	\label{eq:equiv:f_M:f_N}
		f_{M} = \varphi \circ f_{N} \circ \psi \;,
	\end{equation}
	where $\psi$ induces a Boolean automorphism on $\Fi_{I}$. Then, $M \sim N$.
\end{thm}

\begin{proof}
   Since $\psi_*$ is an automorphism on $\Fi_{I}$, for any 
   $\alpha \subseteq \Zint{n}$, there exists a 
   unique $\beta \subseteq \Zint{n}$ such that
	$\pos C_{I}(\alpha)= \psi^{-1}(\pos C_{I}(\beta))$.
   It then follows from~\eqref{eq:f:piecewise:linear} and~\eqref{eq:equiv:f_M:f_N} that
   $\varphi^{-1}(C_{-M}(\alpha) x) = C_{-N}(\beta)\psi(x)$ for $x \in \pos C_{I}(\alpha)$.
   By~\eqref{eq:f:cones} we finally have
   \begin{displaymath}
    \varphi^{-1}(\pos C_{M}(\alpha)) = \pos C_{N}(\beta) \;,
   \end{displaymath}
   which shows that $\varphi$ induces a bijection from $\Ge_M$ to $\Ge_N$.
\end{proof}

There is a slightly weaker converse result.

\begin{thm}
	\label{thm:equivalence:nec}
	Consider two non-degenerate matrices $M, N \in \RE^{n \times n}$ and suppose 
   that $M \sim N$ with homomeorphism $\varphi$. Then, there exists another 
   homeomorphism $\psi: \RE^{n} \to \RE^{n}$ that induces an automorphism
   $\psi_*: \Fi_I \to \Fi_I$ and such that~\eqref{eq:equiv:f_M:f_N} holds.
\end{thm}

\begin{proof}
   The bijectivity of $\varphi_* : \Fi_M \to \Fi_N$ implies that, for each 
   $\alpha \subseteq \Zint{n}$, there exists a unique $\beta \subseteq \Zint{n}$ 
   such that
   \begin{equation} \label{eq:bool:iso}
    \varphi_*(\pos C_{M}(\alpha)) = \pos C_{N}(\beta) \;.
   \end{equation}
   We will explicitly construct a piecewise linear map $\psi:\RE^n \to \RE^n$ 
   and verify~\eqref{eq:equiv:f_M:f_N}. By~\eqref{eq:f:cones} and~\eqref{eq:bool:iso},
   for each $x \in \pos C_{I}(\alpha)$ there	exists an $x' \in \pos C_{I}(\beta)$ such that
	\begin{equation}
		\label{eq:psi:implicit}
		\varphi^{-1}(C_{-M}(\alpha) x) = C_{-N}(\beta) x' \;.
	\end{equation}
   Moreover, $x'$ is unique by the invertibility of $C_{-N}$ (which follows by the
   non-degeneracy of $N$). We define $\psi: x \mapsto x'$ with $x'$ so constructed, so that
   \begin{displaymath}
		\varphi^{-1}(C_{-M}(\alpha) x) = C_{-N}(\beta)\psi(x)
   \end{displaymath}
   for $x \in \pos C_{I}(\alpha)$. This is simply an explicit rewriting of~\eqref{eq:equiv:f_M:f_N}.
	Since $\varphi^{-1}$ is bijective, we can reverse the argument in order to find
   $\psi^{-1}$. The continuity of $\psi$ and $\psi^{-1}$ are simple consequences
   of the continuity of $\varphi$, $\varphi^{-1}$, $f_M$ and $f_N$.

   Finally, note that $\psi$ induces a permutation on $\Ge_I$, and hence an automorphism on $\Fi_I$.
\end{proof}

Condition~\eqref{eq:equiv:f_M:f_N} is the commutative diagram
\begin{displaymath}
\begin{tikzcd}
 \RE^n \arrow{r}{\psi} \arrow[swap]{d}{f_M} & \RE^n \arrow{d}{f_N} \\
 \RE^n                                      & \arrow{l}{\varphi} \RE^n
\end{tikzcd} \;.
\end{displaymath}
It is standard in the literature of singularity theory~\cite{arnold}
and ensures that we can continuously map solutions of the problem $f_M(x) = q$
into solutions of the problem $f_N(x') = \varphi^{-1}(q)$. The requirement on $\psi$ being a
Boolean automorphism implies that $\psi$ maps orthants into orthants, intersections of 
orthants into intersections of orthants, and so forth; and this 
ensures that the complementarity condition is not destroyed by the homeomorphisms.

Linear homeomorphisms preserving the complementarity conditions include permutations of
coordinates and dilations of the coordinate axes~\cite{howe1981}.
\begin{cor}
	\label{cor:equivalent:permutation}
	Let $P \in \RE^{n \times n}$ be a permutation matrix. Then, for any matrix $M \in \RE^{n\times n}$,
   $M \sim N = P^{\top} M P$.
\end{cor}
\begin{proof}
   Let $\varphi(q') = Pq'$ and $\psi(x) = P^\top x$, and note that $\psi$ induces a Boolean
   automorphism on $\Fi_I$. Also, note that $P [x]^{+} = [P x]^{+}$. The proof then follows
   from~\eqref{eq:pwl:f:projection},
   \opt{twoside}{
   	\begin{multline*}
        f_M(x) = [x]^+ - M[-x]^+ = \\ P\left( [P^\top x]^+ - N[-P^\top x]^+ \right)
 		  = \varphi \circ f_{N}  \circ \psi(x) \;.
 	   \end{multline*}
   }
   \opt{onecolumn}{
   \begin{displaymath}
       f_M(x) = [x]^+ - M[-x]^+ = P\left( [P^\top x]^+ - N[-P^\top x]^+ \right)
		  = \varphi \circ f_{N}  \circ \psi(x) \;.
   \end{displaymath}
   }
\end{proof}

\begin{cor}
	Let $D \in \RE^{n \times n}$ be a diagonal matrix with strictly positive entries.
	Then, for any matrix $M \in \RE^{n\times n}$, $M \sim D^{-1} M D$ and 
	$M \sim M D$.
	\label{cor:equivalent:diagonal}
\end{cor}

\begin{proof}
	For the first equivalence, let $\varphi(q') = Dq'$ and $\psi(x) = D^{-1} x$, and note 
   that $\psi$ is a Boolean automorphism on $\Fi_I$. We also have $D [x]^{+} = [D x]^{+}$.
   As in the previous corollary, $f_M(x) = \varphi \circ f_{N}  \circ \psi(x)$. 
	For the second equivalence, let $\varphi$ be the identity map, $\varphi(q) = q$, and 
	$\psi(x) = C_{-D}(\alpha) x$ for $x \in \pos C_{I}(\alpha)$.
	Both $\varphi$ and $\psi$ are continuous and invertible.
	Hence, $f_{M D}(x) =  
	 C_{-M D}(\alpha) x = C_{-M}(\alpha) C_{-D}(\alpha) x = \varphi \circ f_{M} \circ
   \psi(x)$, for some $\alpha \in \Zint{n}$ and 
   the conclusion follows.
\end{proof}

\begin{defn}
	Let $M \in \RE^{n \times n}$ and $\beta \subseteq \Zint{n}$
	be such that $M_{\beta, \beta}$ is nonsingular. The \emph{principal
	pivotal transform} (PPT) of $M$ relative to $\beta$ is the matrix 
	$N \in \RE^{n\times n}$ defined by
	\begin{align*}
		N_{\beta, \beta} & = \left(M_{\beta, \beta}\right)^{-1}
		\\
		N_{\beta, \beta^{c}} & = -\left(M_{\beta, \beta}\right)^{-1} 
		M_{\beta, \beta^{c}}
		\\
		N_{\beta^{c}, \beta} & = M_{\beta^{c}, \beta} \left(M_{\beta, \beta}
		\right)^{-1}
		\\
		N_{\beta^{c}, \beta^{c}} & = M_{\beta^{c}, \beta^{c}} - M_{\beta^{c}, \beta}
		\left(M_{\beta, \beta}\right)^{-1} M_{\beta, \beta^{c}}
	\end{align*}
	By convention, the PPT of $M$ relative to $\emptyset$ is $M$.
	\label{defn:ppt:M}
\end{defn}

\begin{cor}
	Let $N$ be the PPT of $M$ with respect to $\beta$. Then,	$M \sim N$.
	\label{cor:ppt:equivalence}
\end{cor}

\begin{proof}
	If $\beta = \emptyset$ the result holds trivially, so we suppose that
	$\beta \neq \emptyset$. Let $P \in \RE^{n \times n}$ be the permutation 
	matrix such that $P^{\top} \xi = [ \xi_{\beta}^{\top}, \xi_{\beta^{c}}^{\top}]^{\top}$
	and define $\psi(x) = C_{I}(\beta) x$. 
	It follows from \eqref{eq:pwl:f:projection} that

	\begin{align*}
		& P^{\top} ( f_{M} \circ \psi )(x) 
		\\
		& = [P^{\top} C_{I}(\beta) x]^{+} - 
		P^{\top} M P [-P^{\top} C_{I}(\beta) x]^{+}
		\\
		& = 
		\begin{bmatrix}
			- x_{\beta} \\ x_{\beta^{c}}
		\end{bmatrix}^{+} - 
		\begin{bmatrix}
			M_{\beta ,\beta} & M_{\beta, \beta^{c}}
			\\
			M_{\beta^{c}, \beta} & M_{\beta^{c}, \beta^{c}}
		\end{bmatrix}
		\begin{bmatrix}
			x_{\beta} \\ -x_{\beta^{c}}
		\end{bmatrix}^{+}
		\\
		&  = 
		\begin{bmatrix}
				-M_{\beta, \beta} & 0
				\\
				-M_{\beta^{c}, \beta} & I_{\vert \beta^{c} \vert}
			\end{bmatrix}
		\begin{bmatrix}
			x_{\beta} \\ x_{\beta^{c}}
		\end{bmatrix}^{+}
		 -
		\begin{bmatrix}
			-I_{\vert \beta \vert} & M_{\beta, \beta^{c}}
			\\
			0 & M_{\beta^{c}, \beta^{c}}
		\end{bmatrix}
		\begin{bmatrix}
			-x_{\beta} \\ - x_{\beta^{c}}
		\end{bmatrix}^{+}
	\end{align*}
	\begin{align*}
	& = 
	\begin{bmatrix}
		-M_{\beta, \beta} & 0
		\\
		-M_{\beta^{c}, \beta} & I_{\vert \beta^{c} \vert}
	\end{bmatrix} \left( [P^{\top} x]^{+} - P^{\top} N P [- P^{\top} x]^{+} \right)
	\\
	& = \begin{bmatrix}
		-M_{\beta, \beta} & 0
		\\
		-M_{\beta^{c}, \beta} & I_{\vert \beta^{c} \vert}
	\end{bmatrix} P^{\top} f_{N}(x)
\end{align*}
Hence, $f_{N} = \varphi\circ f_{M} \circ \psi$ with
\begin{equation} \label{eq:varphi:pivot}
	\varphi(q') = P
	\begin{bmatrix}
		-M_{\beta, \beta}^{-1} & 0
		\\
		- M_{\beta^{c}, \beta} M_{\beta, \beta}^{-1} & I_{\vert \beta^{c}\vert}
	\end{bmatrix} P^{\top} q' \;.
\end{equation}
Therefore, $N \sim M$.
\end{proof}

\section{Stability of LCPs} \label{sec:stability}

We are now concerned with how stable the structure of the complementary cones of an
$\lcp(M,q)$ is with respect to changes in $M$.

\subsection{Definition and characterization of stability}

Henceforth, we consider the usual topology on the vector space $\RE^{n \times n}$.

\begin{defn}
	A matrix $M \in \RE^{n \times n}$ is \emph{LCP-stable} if there exists a neighborhood
   $U$ of $M$ such that $M \sim N$ for all $N \in U$.
	\label{defn:lcp:stable}
\end{defn}

In order to characterize LCP stability, we take a deeper look into the 
geometric properties of the facets  of the complementary cones.
Recalling that $\mathcal{K}(M)$ denotes the union of all facets of the complementary cones
of $M \in \RE^{n \times n}$. Then, $\mathcal{K}(M) = \bigcup_{\alpha} \boundary
C_{M}(\alpha)$ \cite[p. 511]{cottle} and it breaks down $\RE^{n}$ into 
a finite number of connected regions. 

\begin{defn}
	\label{defn:partition}
	Let $M \in \RE^{n \times n}$.  
	The \emph{partition induced by $M$}, denoted as $\Pa_{M}$, 	
	is the finite family of connected components of $\RE^{n} \setminus \mathcal{K}(M)$.
\end{defn}

A set $Q_{M}^{i} \in \Pa_{M}$ is called a \emph{cell} of $\Pa_{M}$. 
The number of cells, $r$, in the partition induced by $M$ depends on the number of solid 
complementary cones and how they intersect.
The following properties are easy to verify,
\begin{enumerate}[i)]
    \item $Q_{M}^{i} \bigcap Q_{M}^{j} = \emptyset$ for 
       $i \neq j$, $i,j \in \Zint{r}$,
 
    \item $\RE^{n} = \closure \bigcup_{i} Q_{M}^{i}$, \label{it:cover}

    \item $\mathcal{K}(M) = \bigcup_{i} \boundary Q_{M}^{i}$. \label{it:boundary}
\end{enumerate}

The cells of the partition are (non-necessarily convex) open cones 
delimited by the facets of the solid complementary cones of $M$. 
For $P$-matrices, the cells of the partition precisely agree with the
interior of the complementary cones of $M$, see e.g.,~\cite{samelson1958}.

\begin{lem}
	\label{lemma:complementary:partition:M}
	Let $M \in \RE^{n\times n}$ be an $R_{0}$ matrix and let $r$ be the number 
	of sets in $\Pa_{M}$.
   For every $\alpha \subseteq \Zint{n}$ there exists $J \subset \Zint{r}$ 
	such that either
	\begin{displaymath}
		\pos C_{M}(\alpha) = \closure \bigcup_{j \in J} Q_{M}^{j} 
	\end{displaymath}
   (when $\pos C_{M}(\alpha)$ has non-empty interior) or
   \begin{displaymath}
      \pos C_{M}(\alpha) \subset \bigcup_{j \in J} \boundary Q_{M}^{j}
   \end{displaymath}
   (when the interior of $\pos C_{M}(\alpha)$ is empty).
\end{lem}

\begin{proof}
 	 Assume first that $\pos C_{M}(\alpha)$ is a solid cone. Then,
	 there exists $j_{1} \in \Zint{r}$ such that $Q_{M}^{j_{1}} 
	 \subset \pos C_{M}(\alpha)$.
    Let $\mu(E)$ be the Lebesgue measure of $E \subset \RE^{n}$. If
    \begin{displaymath}
	  \mu( \pos C_{M}(\alpha) \setminus Q_{M}^{j_{1}}) = 0 \;,
    \end{displaymath}
    then $\pos C_{M}(\alpha) = \closure Q_{M}^{j_{1}}$, as desired.
    If
    \begin{displaymath}
  	  \mu( \pos C_{M}(\alpha) \setminus Q_{M}^{j_{1}}) > 0 \;,
    \end{displaymath}
    then there exists $j_{2} \in \Zint{r}$, $j_{2} \neq j_{1}$, such that
    $Q_{M}^{j_{2}} \subset \left(\pos C_{M}(\alpha) \setminus Q_{M}^{j_{1}}\right)$.
   As before, if 
   \begin{displaymath}
	   \mu \left(\pos C_{M}(\alpha) \setminus \left( Q_{M}^{j_{1}} 
	   \cup Q_{M}^{j_{2}} \right) \right) = 0 \;,
   \end{displaymath} 
   then $\pos C_{M}(\alpha) = \closure \left( Q_{M}^{j_{1}} \cup 
Q_{M}^{j_{2}} \right)$, as desired. Since $\pos C_{M}(\alpha)$ is strictly pointed
   and $\Pa_{M}$ is finite and the union is dense in $\RE^{n}$ 
   (property \ref{it:cover}) above), the procedure just introduced stops after a 
   finite number of iterations with $J \subset \Zint{r}$ and the result follows. 
   Suppose now that $\pos C_{M}(\alpha)$ is degenerate. Then,
   $\pos C_{M}(\alpha) \subset \bigcup_{j \in J} \boundary Q_{M}^{j} \subset \mathcal{K}(M)$
   for some $J \subset \Zint{r}$.
\end{proof}

\begin{cor}
	\label{cor:partition:cones} $M \sim N$ if, and only if, there exists a bijection 
   $\hat{\varphi}_*:\Pa_M \to \Pa_N$ induced by a homeomorphism $\hat{\varphi}:\RE^{n} \to \RE^{n}$.
\end{cor}

\begin{proof}
	Suppose that $M \sim N$. Then, there is a homeomorphism $\varphi: \RE^{n} \to \RE^{n}$ 
	such that 
   \begin{displaymath}
    \varphi^{-1}(\boundary \pos C_{M}(\alpha)) = \boundary \pos C_{N}(\beta(\alpha))
   \end{displaymath}
	and
   \begin{displaymath}
    \varphi^{-1}(\interior C_{M}(\alpha)) = \interior C_{N}(\beta(\alpha)) \;,
   \end{displaymath}
	for all $\alpha \subset \Zint{n}$ and some 
   \begin{displaymath}
    \beta: \mathcal{P}(\Zint{n}) \to \mathcal{P}(\Zint{n}) \;.
   \end{displaymath}
   Thus, $\varphi^{-1}$ maps connected regions of $\RE^{n} \setminus \mathcal{K}(M)$ into connected regions
   of $\RE^{n} \setminus \mathcal{K}(N)$, thereby inducing a bijection
   between $\Pa_{M}$ and $\Pa_{N}$.

   On the other hand, if $\hat{\varphi}$ induces a bijection $\hat{\varphi}_*$, it follows from 
   Lemma~\ref{lemma:complementary:partition:M} and the continuity of $\hat{\varphi}$, that
   $\hat{\varphi}$ also induces an isomorphism between $\Fi_{M}$ and $\Fi_{N}$.
\end{proof}

Corollary~\ref{cor:partition:cones} implies that, if $M \sim N$, then $\Pa_{M}$
and $\Pa_{N}$ must have the same number of cells.

Consider a matrix $M \in \RE^{n\times n}$ and a vector $-M_{\bscdot, k}$ with $k \in \Zint{n}$,
and let $\mathcal{T}_k(M)$ be the collection of facets for 
which $-M_{\bscdot, k}$ is not a generator,
\opt{twoside}{
 \begin{multline*}
  \mathcal{T}_{k}(M) := \Big\{ S \subset \RE^{n} \mid S = \pos C_{M}(\alpha)_{\bscdot, i^{c}}, \\
 	\alpha \subseteq \Zint{n}, \; k \notin \alpha, \; i \in \Zint{n}\Big\} \;.
 \end{multline*}
}
\opt{onecolumn}{
 \begin{displaymath}
  \mathcal{T}_{k}(M) := \Big\{ S \subset \RE^{n} \mid S = \pos C_{M}(\alpha)_{\bscdot, i^{c}}, \;
 	\alpha \subseteq \Zint{n}, \; k \notin \alpha, \; i \in \Zint{n}\Big\} \;.
 \end{displaymath}
}
\begin{defn}
 Given $M \in \RE^{n \times n}$, we say that $M$ is \emph{weakly degenerate} if
 it has a degenerate complementary cone or if 
 \begin{displaymath}
  -M_{\bscdot, k} \in \bigcup_{S \in \mathcal{T}_{k}(M)} S 
 \end{displaymath}
 for some $k \in \Zint{n}$.
\end{defn}

As an illustration, consider the matrix
\begin{equation}
	\label{eq:example:partition:M}
	M = 
	\begin{bmatrix}
		\frac{1}{2} & \frac{5}{3} & 0
		\\
		1 & 1 & 0
		\\
		-\frac{3}{10} & -1 & 1
	\end{bmatrix} \;.
\end{equation}
It is lengthy but straightforward to verify that there are no degenerate complementary cones,
so $M$ is non-degenerate. However, note that
\opt{twoside}{
 \begin{multline*}
 	\mathcal{T}_{2}(M) = \Big\{ 
 		\pos [I_{\bscdot, 2}, I_{\bscdot, 3}], 
 		\pos [I_{\bscdot, 2}, -M_{\bscdot, 3}],  
 		\pos [I_{\bscdot, 1}, I_{\bscdot, 3}], 
 		\\
 		\pos [-M_{\bscdot, 1}, I_{\bscdot, 3}], 
 		\pos [I_{\bscdot, 1}, -M_{\bscdot, 3}],
 		\\
 		\pos [-M_{\bscdot, 1}, -M_{\bscdot, 3}], 
 		\pos [I_{\bscdot, 1}, I_{\bscdot, 2}], 
  		\pos [-M_{\bscdot, 1}, I_{\bscdot, 2}]
 		\Big\}
 \end{multline*}
}
\opt{onecolumn}{
\begin{multline*}
	\mathcal{T}_{2}(M) = \Big\{ 
		\pos [I_{\bscdot, 2}, I_{\bscdot, 3}], 
		\pos [I_{\bscdot, 2}, -M_{\bscdot, 3}],  
		\pos [I_{\bscdot, 1}, I_{\bscdot, 3}], 
		\pos [-M_{\bscdot, 1}, I_{\bscdot, 3}], 
		\pos [I_{\bscdot, 1}, -M_{\bscdot, 3}],
		\\
		\pos [-M_{\bscdot, 1}, -M_{\bscdot, 3}], 
		\pos [I_{\bscdot, 1}, I_{\bscdot, 2}], 
		\pos [-M_{\bscdot, 1}, I_{\bscdot, 2}]
		\Big\}
\end{multline*}
}
and that $-M_{\bscdot, 2} \in \pos [-M_{\bscdot, 1}, I_{\bscdot, 2}]$, 
so $M$ is a weakly degenerate matrix.
Fig.~\ref{fig:partition} depicts the partitions that $M$ and a perturbation $\tilde{M}$
induce on the unit ball.

\begin{figure}[htpb]
	\centering
	\input{3dPlot.tex} 
	\input{3dPlot_2.tex}
	\caption{Projection onto the unit ball of complementary cones of the matrix 
		$M$ defined in~\eqref{eq:example:partition:M}(left) and a perturbation 
		$\tilde{M}$ (right).
		The generator $-M_{\bscdot, 2} \in \pos [-M_{\bscdot, 1}, I_{\bscdot, 2}]$ 
	(in red). Thus, $M$ is weakly degenerate whereas $\tilde{M}$ is not.}

	\label{fig:partition}
\end{figure}
\begin{lem}
	The set of weakly degenerate matrices is nowhere dense in $\RE^{n \times n}$.
	\label{lem:nowhereDense}
\end{lem}

\begin{proof}
	By definition, a weakly degenerate matrix $M$ belongs to the union of the set of 
	singular matrices and the set
   \opt{twoside}{
 	 \begin{multline} \label{eq:weak:H}
 	 	H : = \left\{ M \in \RE^{n \times n} \mid \text{there exist $k \in \Zint{n}$ and} \right. \\ 
        \left. S \in \mathcal{T}_{k}(M) \text{ such that $-M_{\cdot, k} \in S$} \right\} \;.
 	 \end{multline}
   }
   \opt{onecolumn}{
    \begin{equation} \label{eq:weak:H}
	  	H : = \left\{ M \in \RE^{n \times n} \mid \text{there exist $k \in \Zint{n}$ and
      $S \in \mathcal{T}_{k}(M)$ such that $-M_{\cdot, k} \in S$} \right\} \;.
	 \end{equation}
   }
	It is well know that the set of singular matrices is nowhere dense in 
	$\RE^{n \times n}$. Indeed, the same is true for the set $H$. 
   Consider $\hat{S} \in \RE^{n \times (n-1)}$ such that $\pos \hat{S} \in \mathcal{T}_{k}(M)$.
	It follows that $\det([\hat{S}, -M_{\cdot, k}]) = 0$ and
	$H$ is nowhere dense. To conclude, recall that the finite union
	of nowhere dense sets is again nowhere dense. 
\end{proof}

\begin{thm}
	A matrix $M \in \RE^{n \times n}$ is LCP-stable if, and only if,
	$M$ is not weakly degenerate.
	\label{thm:lcp:stable}
\end{thm}

\begin{proof}
	Let us assume first that some complementary cones of $M$ are degenerate. 
	Since the set of matrices with non-zero principal minors is dense in the 
   set of square matrices, there exists a matrix $N$, arbitrarily close to $M$,
   without degenerate complementary cones. A bijection
   $\varphi_* : \mathcal{G}_{M} \to \mathcal{G}_{N}$ necessarily maps a 
   complementary cone with empty interior onto a complementarity cone with
   non-empty interior. Clearly, $\varphi: \RE^{n} \to \RE^{n}$ cannot be 
   invertible and hence cannot be a homeomorphism. Therefore, $M$ is not LCP-stable.

   Assume now that $M \in H$. Then, there is $J \subset \Zint{n}$ such that,
   for each $k \in J$, there exists a matrix $\hat{S} \in \RE^{n \times n-1}$
   satisfying
   \begin{displaymath}
    \pos [\hat{S}, -M_{\cdot, k}]  = \pos \hat{S} = S \in \mathcal{T}_{k}(M) \;.
   \end{displaymath}
   From Lemma \ref{lem:nowhereDense}, there exists $N \in \RE^{n\times n}$ arbitrarily 
   close to $M$ such that 
   $N$ is not weakly degenerate. 
   In other words, $\det[\hat{S}, -N_{\cdot, k}] \neq 0$, so that 
   $\pos [\hat{S}, -N_{\cdot, k}]$ 
   has non-empty interior. Once again, there is no homeomorphism 
   mapping $\pos [\hat{S}, -M_{\cdot, k}]$ into $\pos [\hat{S}, -N_{\cdot, k}]$. 
   It follows that $M$ is not LCP-stable.

   For sufficiency, let us assume that $M$ is not weakly degenerate.
   There exists a neighborhood $U$ of $M$ such that, for
   every $N \in U$, the partitions $\Pa_M$ and $\Pa_N$ have the same number of cells. Indeed, each cell in 
	$\Pa_M$ is a solid cone by construction and, because there are no
	degenerate cones nor $k \in \Zint{n}$ such that $-M_{\cdot, k} \in S$ for
	some $S \in \mathcal{T}_{k}(M)$, the cells remain solid for sufficiently
   small perturbations. Thus, neither new cells appear nor existing cells disappear in $U$.

   Now, there exist refinements $\hat{\Pa}_M$ and $\hat{\Pa}_N$ of the partitions induced
   by $M$ and $N \in U$, respectively, such that each cell is a convex cone with exactly $n$ generators and
   the number of cells is the same for both partitions. There is thus a one-to-one map
   $\eta: \hat{\Pa}_{M} \to \hat{\Pa}_{N}$ such that
   neighbors are preserved, that is, any two cells in $\hat{\Pa}_{M}$ sharing a 
   $p$-dimensional face are mapped under $\eta$ into cells in $\hat{\Pa}_{N}$ sharing also
   a $p$-dimensional face.
   Let $\hat{Q}_{k}(M)$ and $\hat{Q}_{k}(N) = \eta(\hat{Q}_{k}(M))$ be, 
   respectively, matrices that generate the cells of the new partitions 
   $\hat{\Pa}_M$ and  $\hat{\Pa}_N$, where $k \in \Zint{\hat{r}}$ and $\hat{r} = \vert 
   \hat{\Pa}_{M} \vert = \vert \hat{\Pa}_{N} \vert$.
   Define $\hat{\varphi}: \RE^{n} \to \RE^{n}$ as the piecewise linear function

	\begin{equation}
		\label{eq:homeomorphism:partition}
		\hat{\varphi}(q') = \hat{Q}_{k}(M) \hat{Q}_{k}(N)^{-1} q' \; 
		\text{for $q' \in \pos \hat{Q}_{k}(N)$} \;.
	\end{equation}
	Notice that $\hat{\varphi}$ is continuous by construction, since neighbor cells 
	share $n-1$ generators. That is,
	\begin{displaymath}
		\lim_{\substack{q_i' \to q \\ q_i' \in \pos \hat{Q}_{k}(N)}} 
		\hat{\varphi} (q_i') = \lim_{\substack{q_i' \to q \\ q_i' \in 
		\pos \hat{Q}_{s}(N)}} \hat{\varphi} (q_i')
	\end{displaymath}
	for any $q' \in \pos \hat{Q}_{k}(N) \bigcap \pos \hat{Q}_{s}(N)$. 
   Finally, by similar arguments, $\hat{\varphi}^{-1}$ 
   exists and is also continuous. Furthermore,
   $\hat{\varphi}_*: \Pa_M \to \Pa_N$ is a bijection. 
   Therefore, $M$ and $N$ are LCP-equivalent and,
   since $N$ is arbitrary, we conclude that $M$ is LCP-stable.
\end{proof}

\subsection{Explicit conditions for stability}

Weak degeneracy is connected to the \emph{positive} linear dependence property
of particular subsets of generators of $M$, which 
is difficult to verify for large matrices \cite{regis2016}. On the other hand, efficient
algorithms for verifying linear dependence are readily available. 
This motivates the following.

\begin{cor}
	If all minors of $M \in \RE^{n \times n}$ are different from zero, 
   then $M$ is LCP-stable.
	\label{cor:lcp:stable:principal:minor}
\end{cor}

\begin{proof}
	Note that, if every minor of $M$ is non-zero, 
	then $M$ cannot be weakly degenerate. Indeed, using Laplace's expansion
   we can see that, for any non-empty sets $\alpha, \beta \subseteq \Zint{n}$ such that
	$\vert \alpha \vert = \vert \beta \vert$,
	\begin{equation}
		\label{eq:minor:formula}
		\vert \det M_{\alpha, \beta} \vert = \vert \det [I_{\cdot, \alpha^{c}},
		M_{\cdot, \beta}] \vert \;.
	\end{equation}
	Thus, our assumption on the minors of $M$ implies that any subset of $n$ vectors of 
	\begin{displaymath}
		\{I_{\cdot, 1}, \dots, I_{\cdot, n}, -M_{\cdot, 1}, \dots, -M_{\cdot, n}\}
	\end{displaymath}

	is linearly independent. Consequently, $M$ does not have degenerate cones and
	$-M_{\cdot, k}$ is not in the linear span of $S$ for all 
	$S \in \mathcal{T}_{k}(M)$. That is, $M$ is not weakly degenerate.
\end{proof}

\section{Application to negative-resistance circuits} \label{sec:application}

We will apply the results of the previous sections to the analysis of the circuit
described in Section~\ref{sec:circuit_model}. The application is motivated by~\cite{chua1983}, 
where the authors present a heuristic procedure for constructing negative differential-resistance 
devices using only BJTs and resistors. The design strategy outlined in~\cite{chua1983}
ensures the satisfaction of an analytic condition that is necessary for having multiple solutions.
Still, the question of whether a particular circuit actually exhibits negative resistance is 
settled by simulations. In this section, we employ the LCP formalism to derive necessary and 
sufficient conditions for having a negative resistance.
Our objective is to find values of parameters that 
guarantee bistability of the network, resulting in the hysteretic behavior intrinsic of differential
negative-resistance devices.

\subsection{Circuit equilibria}

\begin{figure*}
	\begin{equation}
		\label{eq:lcp:circuit:M}
	\hat{M} = 
	\begin{bmatrix}
		\frac{G_{0}^{a}+ (1-\alpha_{F}) (G_{2}^{a} + G_{2})}{\alpha_{F}} 
		& 
		\frac{\alpha_{F} (G_{1}^{a} + G_{1}) - 
		(1 - \alpha_{F})(G_{2}^{a} + G_{2})}{\alpha_{F}}
		&
		-G_{1}
		&
		\frac{\alpha_{F} G_{1} - (1 - \alpha_{F})G_{2}}{a_{F}}
		\\
		\frac{\alpha_{R} G_{0}^{a} - (1 - \alpha_{R})(G_{2}^{a} + G_{2})}{\alpha_{R}}
		&
		\frac{G_{1}^{a} + G_{1} + (1 - \alpha_{R})(G_{2}^{a} + G_{2})}{\alpha_{R}}
		&
		-\frac{G_{1}}{\alpha_{R}}
		& 
		\frac{G_{1} + (1 - \alpha_{R}) G_{2}}{\alpha_{R}}
		\\
		-G_{2}
		&
		\frac{\alpha_{F} G_{2}- (1 - \alpha_{F})G_{1}}{\alpha_{F}}
		&
		\frac{G_{0}^{b} + (1 - \alpha_{F}) (G_{2}^{b} + G_{1})}{\alpha_{F}}
		&
		\frac{\alpha_{F} (G_{1}^{b} + G_{2}) 
		- (1 - \alpha_{F})(G_{2}^{b} + G_{1})}{\alpha_{F}}
		\\
		-\frac{G_{2}}{\alpha_{R}}
		&
		\frac{G_{2} + (1 - \alpha_{R}) G_{1}}{\alpha_{R}}
		&
		\frac{G_{0}^{b} - (1 - \alpha_{R})(G_{2}^{b} + G_{1})}{\alpha_{R}}
		&
		\frac{G_{1}^{b} + G_{2} + (1 - \alpha_{R})(G_{2}^{b} + G_{1})}{\alpha_{R}}
	\end{bmatrix} \;. \tag{$\star$}
\end{equation}
\end{figure*}

Recall that the equilibria of~\eqref{eq:Sigma} are characterized by the solutions to the $\lcp(M, q)$,
where $M$ and $q$ are given in~\eqref{eq:lcp:M} and~\eqref{eq:lcp:q}. The matrix inverse
in~\eqref{eq:lcp:M} complicates  the computation of equilibria and the ensuing bifurcation
analysis. Fortunately, the results of the previous sections ease the computational burden and open the
way for an analytical treatment.

Let $A, B, C \in \RE^{4\times 4}$ be nonsingular and defined by~\eqref{eq:lcs:param:1} 
and~\eqref{eq:lcs:param:2}, and let $N := -C A^{-1} B$. By applying a pivotal transformation to $N$ 
with $\beta = \{1, 2, 3, 4\}$ we obtain
\begin{displaymath}
 \hat{M} := -B^{-1} A C^{-1} \;.
\end{displaymath}
By Corollary~\ref{cor:ppt:equivalence}, $N \sim \hat{M}$. 
Let $M$ and $q$ be defined by~\eqref{eq:lcp:M},~\eqref{eq:lcp:q},~\eqref{eq:lcs:param:1} and~\eqref{eq:lcs:param:2}.
If $N$ is stable, $M = N + D \sim \hat{M}$ for $\|D\|$ small enough. Moreover, the solutions of the $\lcp(M,q)$ 
can be continuously mapped into the solutions of the $\lcp(\hat{M},\hat{q})$ with 
\begin{displaymath}
 \hat{q} := B^{-1} \left( A C^{-1} E_{2} s - E_{1} r \right) \;,
\end{displaymath}
where the expression follows from $\hat{q} = \varphi^{-1}(q)$ with $\varphi$ as in~\eqref{eq:varphi:pivot}.

The interest of $M \sim \hat{M}$ is that we no longer require to compute $M$, which involves
the inverse of $A$, defined by~\eqref{eq:lcs:param:1} and~\eqref{eq:lcs:param:2}. Instead, we require to compute
$\hat{M}$, which is easy, since $B$ is block diagonal and $C = -I_{4}$. Indeed, $\hat{M}$ is given explicitly
in Equation~\eqref{eq:lcp:circuit:M}.

\subsection{Circuit bifurcations}

We regard the potential difference $r = P_{1} - P_{0}$ as a bifurcation or control parameter. We wish to find a set of
resistor values for which, by changing $r$, the circuit transitions from global stability to bistability.
The following corollary gives sufficient conditions for attaining this objective.

\begin{cor}[{\cite[Corollary 4.3]{howe1981}}]
Let $\hat{M} \in \RE^{n \times n}$ be an $R_{0}$-matrix such that $\deg_{\hat{M}}(q) = 1$ for some
$q \in \RE^{n}$ and $f_{\hat{M}}$ has only one negative index. 
Then, there is $\alpha_{0} \subseteq \Zint{n}$ such that,
for any $\hat{q} \in \interior \pos C_{\hat{M}}(\alpha_{0})$, the
$\lcp(\hat{M},\hat{q})$ has exactly three solutions, whereas if $\hat{q} \in \RE^{n} 
\setminus \pos C_{\hat{M}}(\alpha_{0})$, the solution is unique.
\label{cor:hysteresis}
\end{cor}

With an explicit expression for $\hat{M}$, it is straightforward to check the conditions of Corollary~\ref{cor:hysteresis}.
 
\begin{prop} \label{prop:hysteresis}
 Let $\hat{M} \in \RE^{4\times 4}$ be given by~\eqref{eq:lcp:circuit:M} with all parameters finite and positive,
 and with $\alpha_F, \alpha_R \in (0,1)$. Define
 \begin{multline} \label{eq:gamma}
  \gamma := G_{0}^{a} G_{0}^{b} + (1 - \alpha_{F}) \left( G_{0}^{a} G_{2}^{b} + G_{0}^{b} G_{2}^{a} + (1 - \alpha_{F})G_{2}^{a} G_{2}^{b} \right) \\
	+ (1 - \alpha_{F}) \left( G_{0}^{a} + (1-\alpha_{F}) G_{2}^{a} \right) G_{1} \\
	- \left( (2\alpha_{F} - 1) G_{1} - (1-\alpha_{F}) \left( G_{0}^{b} + 
	(1-\alpha_{F}) G_{2}^{b} \right) \right) G_{2} \;.
 \end{multline}
 If:
 \begin{itemize}
  \item $\gamma > 0$, then the $\lcp(\hat{M},\hat{q})$ has a unique solution, regardless of $\hat{q}$.
  \item $\gamma < 0$, then the $\lcp(\hat{M},\hat{q})$ may have a unique solution 
	  or three solutions, depending
   on the specific value of $\hat{q}$.
  \item $\gamma = 0$, then $\hat{M}$ is LCP-unstable.
 \end{itemize}  
\end{prop}

\begin{proof}
 Let $\hat{\alpha}_{0} = \{1,3\}$. After lengthy but simple computations, we find that the minor 
 $\det \hat{M}_{\hat{\alpha}_{0}, \hat{\alpha}_{0}} > 0$ if, and only if, $\gamma > 0$;
 whereas the remaining principal minors are always positive. Thus, for $\gamma > 0$, $\hat{M}$ is
 a matrix of class $P$ and the solution of the complementarity problem is unique, regardless of $\hat{q}$.
 If, on the other hand, $\gamma < 0$, then the conditions of Corollary~\ref{cor:hysteresis} hold
 and the result follows. Finally, if $\gamma = 0$, then $\hat{M}$ is degenerate and hence LCP-unstable.
\end{proof}

For real transistors we have $\alpha_{F} \in [0.8, 0.99]$ and $\alpha_{R} \approx 0.5$~\cite{getreu1978}. 
It thus follows that we can make $\gamma < 0$ by setting $G_{1}$ and $G_{2}$ large enough, whenever 
the remaining conductances stay finite. Regarding the condition
$\hat{q} \in \interior \pos C_{\hat{M}}(\hat{\alpha}_{0})$ in Corollary~\ref{cor:hysteresis},
we will verify it numerically once the resistor parameters have been chosen. 
For concreteness, we consider all capacitors with the same value of $100\,\mu$F, 
$\alpha_{F} = 0.99$, $\alpha_{R} = 0.5$, and $s=0.7$ V. The values of the resistances $R_{k}^{j}$, 
$k \in \{0,1,2\}$, $j \in \{a, b\}$ are selected arbitrarily, whereas $G_{1}$ is chosen in such a way that
the factor associated with $G_{2}$ in~\eqref{eq:gamma} is strictly negative.
The chosen values are shown in Table~\ref{tab:param:value}.
 
\begin{table}[t]
	\centering
	\begin{tabular}{ccccc}
		$R_{0}^{a} = 100\,\Omega$ 
		& 
		$R_{1}^{a} = 2.2\, \mathrm{k}\Omega$ 
		& 
		$R_{2}^{a} = 100\, \Omega$
		\\
		\hline
		$R_{0}^{b} = 100\, \Omega$ 
		& 
		$R_{1}^{b} = 100\, \Omega$ 
		& 
		$R_{2}^{b} = 10\, \mathrm{k}\Omega$
		\\
		\hline
		& $R_{1} = 10\, \Omega$ & $R_{2} = 0\,\Omega - 1\, \mathrm{k} \Omega$
		\\
	\end{tabular}
	\caption{Values of the parameters used for the circuit of Fig.~\ref{fig:circuit:nr}.}
	\label{tab:param:value}
\end{table}
With the setting described in the previous paragraph, we can use $G_2$ to set the sign
of $\gamma$. Fig.~\ref{fig:lcp:region:feasible} shows the region of the $(R_{2},r)$-plane for 
which $\gamma < 0$ and $C_{\hat{M}}^{-1}(\hat{\alpha}_{0}) \hat{q} \in \interior \RE_{+}^{4}$.
For instance, the parameters $R_{2} = 20\,\Omega$ and $r = 1.2$ V lie in the interior
of the feasible region. It is easy to verify that, for the chosen parameters, 
$\gamma < 0$ so, by Proposition~\ref{prop:hysteresis}, the $\lcp(\hat{M}, \hat{q})$
has three solutions. 

For the parameters of Table~\ref{tab:param:value} we have
\begin{equation}
	N = - C A^{-1} B= 
	\begin{bmatrix}
		60.12 & -4.83 & 87.66 & -45.17
		\\
		1.92 & 48.50 & 137.58 & -92.83
		\\
		39.80 & 3.16 & 104.60 & -49.62
		\\
		41.86 & -44.98 & -42.36 & 48.15
	\end{bmatrix}.
	\label{eq:circuit:N}
\end{equation}
It is a simple exercise to verify that all the minors of $N$ are different from zero. 
Thus, it follows from Corollary~\ref{cor:lcp:stable:principal:minor} that $N$ is LCP-stable. 
The solutions of the $\lcp(\hat{M},\hat{q})$ can be continuously mapped to the
solutions of the original $\lcp(M,q)$ whenever $\|D\|$ is small enough.
For instance, the parameters $R_{R} = R_{F} = 1\,\Omega$ yield an
$R_{0}$-matrix $M$ such that $\det M_{\alpha_{0}, \alpha_{0}} < 0$ and $q \in \interior
C_{M}(\alpha_{0})$ with $\alpha_{0} = \{2, 4\}$; while the rest of the principal minors are positive.
Therefore, by Corollary \ref{cor:hysteresis} there exists three equilibria for such parameters.

\begin{figure}
	\centering
	\input{feasible_region.tex}
	\caption{Region of parameters (black-filled) in the $(R_{2},r)$-plane
		for which $\ind_{\hat{M}}(x) = -1$, $x \in \interior \pos C_{I}(\hat{\alpha}_{0})$
      and $\hat{q} \in \interior \pos C_{\hat{M}}(\hat{\alpha}_{0})$. We have set
      $\hat{M}$ as in~\eqref{eq:lcp:circuit:M}, $\hat{\alpha}_{0} = \{1,3\}$ and used the  
      parameters shown in Table~\ref{tab:param:value}. The red line indicates the boundary of the
      feasible region for the original problem with $D$ as in~\eqref{eq:lcs:param:1} and $R_{R} = R_{F} = 1\,\Omega$.}
	\label{fig:lcp:region:feasible}
\end{figure}

\subsection{Asymptotic properties of the circuit}

It only rest to show that, with the chosen parameters, the network shows bistability.
To that end, we use the  recent approach of dominant systems developped in~\cite{forni2019,miranda2018}. 
Specifically, with our selection of parameters, the system~\eqref{eq:lcs:1}-\eqref{eq:lcs:2} is strictly $1$-passive from the 
``input'' $z$ to the ``output'' $w$. That is, it satisfies the $\lambda$-parameterized linear matrix inequality
\begin{displaymath}
	\begin{bmatrix}
		A^{\top} P + P A + 2 \lambda P + \varepsilon I & P B_{1} - C^{\top}
		\\
		B_{1}^{\top} P - C & -D - D^{\top}
	\end{bmatrix} \leq 0
\end{displaymath}
for values of $\lambda \in [105, 135]$ and
\begin{displaymath}
P = 
\begin{bmatrix}
	-1.87 & -5.96 & -2.42 & 3.65
	\\
	-5.96 & -5.66 & -5.90 & -0.18
	\\
	-2.42 & -5.90 & -0.86 & 4.42
	\\
	3.65 & -0.18 & 4.42 & 4.67
\end{bmatrix} \times 10^{-3} \;,
\end{displaymath}
which has inertia $(1, 0, 3)$.
Now, since the complementarity condition \eqref{eq:lcs:3} is incrementally passive ($0$-passive in the language of \cite{forni2019}), we conclude that the closed-loop 
\eqref{eq:lcs:1}-\eqref{eq:lcs:3}
is $1$-passive with rate $\lambda > 0$ and therefore $1$-dominant. 
Hence, the asymptotic behavior of the network 
is topologically equivalent to the asymptotic behavior of a $1$-dimensional system.
Further, standard computations show that the trajectories are bounded and there is
one unstable equilibrium, leading us to the desired conclusion.
Fig.~\ref{fig:sim:1} shows the time evolution of the states of the 
network for the choice of parameters described above when an additional current source 
is connected between the terminals $P_{1}-P_{0}$, confirming the desired bistable behavior. 

\begin{figure}
	\centering
	\input{circuit_simulation_1.tex}
	\caption{Time trajectories of states $\xi_{1}, \xi_{2}$ (top) and 
	$\xi_{3}, \xi_{4}$ (middle) 
	of the circuit of Fig.~\ref{fig:circuit:nr}. The extra current input $I$ (bottom)
	confirms the bistable nature of the circuit.}
	\label{fig:sim:1}
\end{figure}

\section{Classification of two-dimensional LCPs} \label{sec:classification}

In this section, we apply the results on LCP equivalence and stability in order to classify
all possible cone configurations of 2-dimensional LCPs. Our incentive for studying such a
low-dimensional problem in full detail comes from the fact that, in some cases, the singularities of an
LCP can be assesed by studying a lower-dimensional one. 

\begin{exmp}
 Let
 \begin{displaymath}
  M =
   \begin{bmatrix}
    M_{11} & M_{12} \\ 0 & M_{22}
   \end{bmatrix}
  \quad \text{and} \quad 
  q = 
   \begin{bmatrix}
    q_1 \\ q_2
   \end{bmatrix} \;,
 \end{displaymath}
 and consider a triangular $\lcp(M,q)$. Suppose that $M_{22}$ is a $P$-matrix. 
 Then, $x_2$ is uniquely given by 
 \begin{displaymath}
  x_2 = f_{M_{22}}^{-1}(q_2) \;.
 \end{displaymath}
 Moreover, since $M_{22}$ is a $P$-matrix, $f_{M_{22}}$ is regular.
 The problem can be then rewritten as
 \begin{displaymath}
  f_{M_{11}}(x_1) = \bar{q}_1 \;, \quad \bar{q}_1 = M_{12} [-f_{M_{22}}^{-1}(q_2)]^+ + q_1 \;,
 \end{displaymath}
 from where we see that the nonsmooth singularities of the $\lcp(M,q)$ are determined by the
 nonsmooth singularities of the lower-dimensional $\lcp(M_{11},\bar{q}_1)$.
\end{exmp}

Using the characterization of stable LCPs provided by Theorem~\ref{thm:lcp:stable},
we start by enumerating all the unstable matrices. Such enumeration is performed by 
writing down all the conditions leading to weakly degenerate matrices.
For $n = 2$, such classification is relatively easy since, in $\RE^{2 \times 2}$,
weakly degenerate matrices either satisfy $-M_{\cdot, k} 
\in \pos I_{\cdot, k}$, for $k \in \Zint{2}$, or they have degenerate complementary
cones. 
Moreover, it follows from Corollary~\ref{cor:equivalent:diagonal}
that any matrix $N \in \RE^{2 \times 2}$ is LCP-equivalent to the matrix
\begin{equation}
 \label{eq:polar:matrix}
 M (\theta_{1}, \theta_{2}) = 
 \begin{bmatrix}
	r_{1} \cos (\theta_{1} + \pi) & r_{2} \cos (\theta_{2} + \frac{3 \pi}{2})\\
   r_{1} \sin (\theta_{1} + \pi) & r_{2} \sin ( \theta_{2} + \frac{3 \pi}{2})
 \end{bmatrix} \;,
\end{equation}
where $r_{1}, r_{2} \in \{0, 1\}$ and $(\theta_{1}, \theta_{2}) \in [0, 2 \pi) \times 
[0, 2 \pi)$ are angles measured in counterclockwise direction with respect to the rays 
$\pos I_{\cdot, 1}$ and $\pos I_{\cdot, 2}$, respectively.
In what follows, we use the representation~\eqref{eq:polar:matrix} for enumerating all 
$2 \times 2$ unstable matrices. First we consider the class of strongly degenerate matrices, 
for which at least one column is zero. It is characterized by the conditions
\begin{equation} \label{eq:class:0}
\begin{aligned}
	r_{1} r_{2} & = 0 \\
	(\theta_{1}, \theta_{2}) & \in [0, 2\pi) \times [0, 2 \pi) 
\end{aligned} \;.
\end{equation}
The second class of unstable matrices encompasses the strongly degenerate matrices for 
which at least one complementary cone is a subspace of $\RE^{n}$. This class is 
characterized by the conditions
\begin{equation}
	\begin{aligned}
		r_{1} = r_{2} & = 1 \;,
		\\
		(\theta_{1}, \theta_{2}) & \in \Theta^+
	\bigcup \left( \left\{ \frac{3 \pi}{2} \right\} \times [0, 2 \pi]  \right)
	\\
	& \quad \bigcup \left( [0, 2 \pi) \times \left\{ \frac{\pi}{2} \right\} \right) \;,
	\end{aligned} 
	\label{eq:class:1}
\end{equation}
where 
\begin{displaymath}
	\Theta^+ = \left\{ (\theta_1, \theta_{2}) \mid \theta_{1} \in [0, 2\pi), 
	\theta_{2} = \left( \theta_{1} + \frac{\pi}{2} \right) 
\bmod{2 \pi} \right\} \;.
\end{displaymath}
Certainly, matrices satisfying either~\eqref{eq:class:0} or~\eqref{eq:class:1}, 
lie outside of the $R_{0}$-class. 

Finally, we consider the class of unstable $R_{0}$-matrices embracing the remaining degenerate 
and weakly degenerate matrices. This class is characterized by the conditions
\begin{equation}
	\begin{aligned}
		r_{1} = r_{2} & = 1 \;,
		\\
		(\theta_{1}, \theta_{2}) & \in \Theta^{-} 
		\bigcup
		\left( \left\{ \frac{\pi}{2} \right\} \times [0, 2 \pi) \right)
		\\
	& \quad \bigcup \left( [0,  2 \pi)] \times \left\{ \frac{3 \pi}{2} \right\} \right) \;,
	\end{aligned}
	\label{eq:class:2}
\end{equation}
where 
\begin{displaymath}
	\Theta^{-} := \left\{ (\theta_{1}, \theta_{2}) \mid \theta_{1} \in [0, 2\pi), 
	\theta_{2} = \left( \theta_{1} - \frac{\pi}{2} \right) \bmod{2 \pi} \right\} \;.
\end{displaymath} 

The list of unstable normal forms complements the stable normal forms enumerated
in~\cite{castanos2020}. However, the list~\eqref{eq:class:0}-\eqref{eq:class:2} is more 
fundamental, in the sense that representatives of all stable classes can be obtained
as perturbations of conditions \eqref{eq:class:0}-\eqref{eq:class:2}.

\begin{figure}[ht]
	\centering
	\input{bifDiagram_Complete.tex}
	\caption{Classification of $2 \times 2$ matrices~\eqref{eq:polar:matrix}. 
	Line segments and points with the same color, style and 
	shape are LCP-equivalent.
	In the central square, regions of stable matrices ($C_{1}$--$C_{5}$)
	are delimited by the unstable matrices 
	satisfying \eqref{eq:class:0}-\eqref{eq:class:2}. 
	}
	\label{fig:bifDiag}
\end{figure}
\begin{figure*}[tb]
	\centering
	\input{classificationForest.tex}
	\caption{Complementary cone configurations of representative members of the classes
	 depicted in the central square in Fig.~\ref{fig:bifDiag}. Highly-degenerate configurations (circles and 
    squares) break down into less-degenerate configurations (solid and dashed lines), which in turn break down
    into the stable configurations ($C_1-C_5$). Red rays denote degenerate cones.}
	\label{fig:unstable:classes}
\end{figure*}

Fig.~\ref{fig:bifDiag} illustrates the partition of the space of $2 \times 2$ matrices at 
different levels. Line segments and points with the same color, style and shape are LCP-equivalent.
Unstable matrices satisfying~\eqref{eq:class:0},~\eqref{eq:class:1}, and~\eqref{eq:class:2}
are depicted using dotted, dashed, and solid lines, respectively.
The double annular mark at the bottom left corner represents the zero matrix 
($r_{1} = r_{2} = 0$). Perturbations of this matrix will lead to any of the other three cases: a point on
the left vertical line characterizing the cases where $r_{1} = 0$, a point on the bottom horizontal line 
characterizing the cases where $r_{2} = 0$, or a point on the central square for which 
$r_{1} \neq 0$ and $r_{2} \neq 0$.
Connected white regions on the central square correspond to stable matrices. 
Notice that stable matrices belong to the $R_{0}$-class, which is 
open (see, e.g.,~\cite{doverspike1982}). It turns out that there are five different 
classes of stable matrices\footnote{In~\cite{castanos2020}, the classes 
$C_2$ and $C_4$ were mistakenly taken to be the same one.} in $\RE^{2 \times 2}$. The complementary cones of 
representative members of each class are shown in 
Fig.~\ref{fig:unstable:classes}.

\section{Conclusions} \label{sec:conclusions}

A notion of equivalence between linear complementarity problems of the same dimension was 
presented. At its fundamental level, two LCPs are equivalent if there are topological and Boolean-algebraic
isomorphisms that will take one into the other. The proposed notion allows us to analyze and classify
solutions of linear complementarity problems and hence to study local bifurcations in linear complementarity systems. 

A complete classification of LCPs in the case $n = 2$ was presented, where five stable classes were identified.
The results on LCP-equivalence were also applied to the design of bifurcations
of equilibria in a class of circuits in order to establish a pre-specified behavior.

Future work considers the extension of the notion of equivalence to allow complementarity problems 
of different dimensions, as well as extensions to other types of complementarity problems
and general bifurcations on linear complementarity systems.





\bibliographystyle{IEEEtran}

