# Analysis of the implicit Euler time-discretization of passive linear descriptor complementarity systems

Bernard Brogliato, Alexandre Rocca

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bernard Brogliato, Alexandre Rocca. Analysis of the implicit Euler time-discretization of passive linear descriptor complementarity systems. Journal of Nonsmooth Analysis and Optimization, 2022, Volume 3, paper 7269. 10.46298/jnsao-2022-7269. hal-03169123v3

HAL Id: hal-03169123

https://inria.hal.science/hal-03169123v3

Submitted on 12 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# ANALYSIS OF THE IMPLICIT EULER TIME-DISCRETIZATION OF PASSIVE LINEAR DESCRIPTOR COMPLEMENTARITY SYSTEMS

Bernard Brogliato*

Alexandre Rocca†

**Abstract** This article is largely concerned with the time-discretization of descriptor-variable systems coupled with complementarity constraints. They are named linear descriptor complementarity systems (LDCS). More specifically passive LDCS with minimal state space representation are studied. The Euler implicit discretization of LDCS is analysed: the one-step non-smooth problem (OSNSP), that is a generalized equation, is shown to be well-posed under some conditions. Then the convergence of the discretized solutions is studied. Several examples illustrate the applicability and the limitations of the developments.

## 1 INTRODUCTION

The analysis of non-smooth dynamical systems with set-valued right-hand sides satisfying maximal monotone properties has been the object of many contributions, see [26] for a survey. Linear and nonlinear complementarity dynamical systems, various kinds of projected dynamical systems, of differential variational inequalities, differential inclusions with maximal monotone right-hand sides, Moreau's sweeping processes, and some switching dynamics can be recast in this class of set-valued dynamical systems. In this article we focus on singular, or linear descriptor complementarity systems (LDCS). They may be viewed either as an extension of "classical" descriptor-variable systems which are widely used in Automatic Control and Circuit Theory, or of "classical" linear complementarity systems. They can also be seen as differential-algebraic equations (DAEs) coupled with complementarity constraints. Motivations for studying LDCS are in circuits with set-valued components [2, Chapters 3 and 5] [32], and in a more abstract framework in state-dependent switching DAEs [71, 70], after a kind of convexification (or filling-in the graph) procedure has been made. They are also closely related to DAEs whose equality constraints are nonsmooth, as they occur in chemistry or thermo-dynamical systems [78, 79, 77].

It has to be noted that the study of LDCS has not received a lot of attention yet, and is still at an embryonic stage of analysis. In particular the analysis of the implicit Euler method, which has received a lot of attention in several fields [10, 14, 57, 55, 2, 3, 67, 69], is an open issue for LDCS. Our goal is to analyze the well-posedness of the one-step nonsmooth problem (OSNSP) obtained after an implicit (backward) Euler discretization, and the convergence of the discrete-time solutions towards a continuous-time limit which is a solution to the continuous-time system. The second step is useful not only for proving the existence of solutions to the continuous-time system, but also, perhaps most importantly, to show that the proposed time-discretization makes sense for numerical simulations.

This work was partially supported by the INRIA Project Lab (IPL) *ModeliScale* https://team.inria.fr/modeliscale/, and by the Fonds Unique Interministériel (FUI) *ModeliScale*. It has not been presented at any conference or workshop. The authors also want to thank Dr Vincent Acary for insightful comments and discussions.

<footnote>*Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France. bernard.brogliato@inria.fr</footnote>

<footnote>†INRIA, Université de Rennes, Irisa 35000 Rennes, France. alexandre.rocca@inria.fr</footnote>
Passivity properties (passive systems, positive real transfer functions, Lur'e equations) are pivotal throughout the article. In Section 2, a class of passive LDCS (which extends descriptor variable systems with a positive real transfer function) is introduced. A specific Weierstrass form associated with these systems, when the state-space realization is minimal, is used. The well-posedness of the OSNSP associated with the implicit Euler discretization, is analysed in Section 3. Convergence of the discrete-time solutions is studied in Sections 4 and 5. Examples that illustrate the theoretical developments, are given in Section 6. Conclusions are drawn in Section 7. The appendix is dedicated to recall various mathematical tools.

**Notation and definitions:** For any vector $x \in \mathbb{R}^n$ and any matrix $M \in \mathbb{R}^{m \times n}$, $\|x\|$ is the Euclidean norm and $\|M\|$ is the Frobenius norm, which are compatible norms [16, Proposition 9.3.5], i.e., $\|Mx\| \le \|M\| \|x\|$. Let $M \in \mathbb{R}^{n \times m}$, then $\mathrm{Im}(M)$ is its range, $\mathrm{Ker}(M)$ is its null space. We use $\langle x, y \rangle = x^\top y$, so $\langle x, x \rangle = \|x\|^2$. Positive definite matrix: $M > 0$ if $x^\top M x > 0$ for all $x \neq 0$, positive semidefinite matrix: $M \ge 0$ if $x^\top M x \ge 0$ for all $x$ (such $M$ is not necessarily symmetric). The maximum singular value is denoted as $\sigma_{\max}(M)$, and the minimum and maximum eigenvalues as $\lambda_{\min}(M)$ and $\lambda_{\max}(M)$. $M_{\bullet i}$ denotes the $i$-th column of $M$. More mathematical results and definitions are given in Appendices A to G.

## 2 THE CLASS OF PASSIVE LDCS

Let us consider the following singular dynamical system, that we may name a linear descriptor complementarity system (LDCS):

$$ (2.1) \qquad \begin{cases} P\dot{x}(t) = Ax(t) + B\lambda(t) + E(t), & \text{a.e. } t \ge 0 \\ 0 \le \lambda(t) \perp w(t) = Cx(t) + D\lambda(t) + F(t) \ge 0, & \text{for all } t \ge 0 \end{cases} $$

with $x(t) \in \mathbb{R}^n$, $\lambda(t) \in \mathbb{R}^m$, $A, B, C, D$ constant matrices of appropriate dimensions, $P \in \mathbb{R}^{n \times n}$ has rank $p < n$. Dissipativity is a fundamental property in Systems and Control [25, 87, 88]. The material in this article relies on a result in [28, 32, 51] on a special Weierstrass form for passive descriptor variable systems [25, Section 3.1.7] which possess a minimal state-space realization [28, Theorem 3.1]. The passivity of the quintuple $(P, A, B, C, D)$ means the passivity of the operator $\lambda \mapsto w$, see Appendix D. Assuming that $E(t) = 0$ and $F(t) = 0$, passivity and minimality, this special form writes as (see [51] [32, Equation (16)] [28, Proposition A.3]):

$$ (2.2) \qquad \left\{ \begin{aligned} & (a) \quad \dot{x}_1(t) = A_1 x_1(t) + B_1 \lambda(t) \\ & (b) \quad \dot{x}_2(t) = z(t) + B_2 \lambda(t) \\ & (c) \quad 0 = x_2(t) + B_3 \lambda(t) \\ & (d) \quad 0 \le \lambda(t) \perp w(t) = C_1 x_1(t) + C_2 x_2(t) + C_3 z(t) + D \lambda(t) \ge 0, \end{aligned} \right. $$

where the quadruple $(A_1, B_1, C_1, \hat{D})$ is itself passive, $\hat{D} \stackrel{\Delta}{=} D - C_3B_2 - C_2B_3$, $x_1(t) \in \mathbb{R}^{n_1}$, $x_2(t) \in \mathbb{R}^{n_2}$, $z(t) \in \mathbb{R}^{n_2}$, $n = n_1 + 2n_2$, $\lambda(t) \in \mathbb{R}^m$, and there exists $X_3 = X_3^\top \ge 0$ such that $B_3^\top X_3 = -C_3$. The matrix $X_3$ is unique when the state-space realization is minimal [28, Theorem 5.4]. The LMI in (d.2) holds. Due to the passivity and the complementarity conditions which imply $\lambda(t)^\top w(t) = 0$ for all times, it follows that for all $t \ge 0$ one obtains $V(x_1(t), x_2(t)) \stackrel{\Delta}{=} \frac{1}{2}(x_1^\top X_1 x_1 + x_2^\top X_2 x_2)(t) \le V(x_1(0), x_2(0))$, [32] and [28, Proof of Theorem 5.1], with $X_1 > 0$ if the system is strongly SPR and minimal [28, Theorem 5.4], see Appendix D. In the latter case $x_1(\cdot)$ is bounded, $A_1$ is Hurwitz<fnref n="1" />, $\hat{D} > 0$ and $B_3$ is full row rank ($\Rightarrow n_2 \le m$) [28, Proposition A.4].

<footnote n="1">¹That is, all its eigenvalues have a negative real part.</footnote>
Remark 2.1. Notice that the differentiation index (shortly the index in the following) of the DAE part of the system (2.2) is equal to 2, since one needs to differentiate the third line twice (assuming $\lambda(t) = 0$) to recover the variable $\dot{z}(t)$. It is clear from (2.2) that the LDCS (2.1) may be seen as a DAE coupled with nonsmooth constraints that take the form of complementarity conditions. It is noteworthy that singular zero order sweeping processes (named in ZOSwP [26]) have been studied [85, 6, 9], but they form another class of singular systems than (2.1) or (2.2).

Remark 2.2. The system in (2.1) can be rewritten equivalently as:

$$ (2.3) \qquad P\dot{x}(t) - Ax(t) - E(t) \in B(\mathcal{N}_{K^*} + D)^{-1}(-Cx(t) - F(t)) $$

with $K = K^* = \mathbb{R}_+^m$, which is, as pointed out in [26, section 3.16], a class of singular systems that has not yet received much attention. The same applies to (2.2) which is rewritten as:

$$ (2.4) \qquad \begin{pmatrix} I_{n_1} & 0 & 0 \\ 0 & I_{n_2} & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \dot{z}(t) \end{pmatrix} - \begin{pmatrix} A_1 & 0 & 0 \\ 0 & 0 & I_{n_2} \\ 0 & I_{n_2} & 0 \end{pmatrix} \begin{pmatrix} x_1(t) \\ x_2(t) \\ z(t) \end{pmatrix} \in \begin{pmatrix} B_1 \\ B_2 \\ B_3 \end{pmatrix} (\mathcal{N}_{K^*} + D)^{-1}(-w(t)). $$

Exogenous terms can also be considered in (2.2) and (2.4), though the analyses in Sections 3 and 4 are made without such terms. Also, an alternative way of analysing (2.2) is presented in Section 5, incorporating exogenous signals. An example is treated in Section 6.1 where exogenous terms are taken into account in the special form (2.2), see (6.3) and (6.4). In the case $D = 0$, both (2.4) and (2.3) can be rewritten, under some assumptions, as a singular first-order sweeping process (FOSwP [26]). The optimal control of FOSwP as treated in [42, 41, 36], relies on the analysis of the time-discretization and the convergence of solutions. Such an important issue could be tackled for singular FOSwP, where the results in this article may serve as a starting point.

### 3 THE EULER IMPLICIT SCHEME AND THE OSNSP

In the sequel of the article, we let $t \in [0, T]$, $0 < T < +\infty$, $h = T/N$, $N \in \mathbb{N} \setminus \{0\}$, $t_k = kh$, $k \in \{0, 1, \dots, N\}$, and $f_k = f(t_k)$. The studied implicit Euler time-discretization of (2.2) is as follows:

$$ (3.1) \qquad \left\{ \begin{array}{l} x_{1,k+1} = x_{1,k} + hA_1x_{1,k} + hB_1\lambda_{k+1} \\ x_{2,k+1} = x_{2,k} + hz_{k+1} + hB_2\lambda_{k+1} \\ 0 = x_{2,k+1} + B_3\lambda_{k+1} \\ 0 \le \lambda_{k+1} \perp w_{k+1} = C_1x_{1,k+1} + C_2x_{2,k+1} + C_3z_{k+1} + D\lambda_{k+1} \ge 0, \end{array} \right. $$

for all $k \in \{0, 1, \dots, N-1\}$, where $A_1x_{1,k}$ could be replaced by $A_1x_{1,k+1}$ without modifying much the next developments. Nevertheless we choose to call the method in (3.1) an implicit method despite of the fact that the linear term is explicitly calculated, because our objective is to focus solely on the set-valued part of the system (i.e., the complementarity conditions). The set of equations in (3.1) makes a mixed linear complementarity problem (MLCP). After few calculations from (3.1) the following OSNSP is obtained (where the variable $hw_{k+1}$ is used instead of $w_{k+1}$) as an LCP:

$$ (3.2) \qquad \begin{array}{l} 0 \le \lambda_{k+1} \perp \bar{D}_h \lambda_{k+1} + h\beta_k \ge 0 \Leftrightarrow \bar{D}_h \lambda_{k+1} + h\beta_k \in -\mathcal{N}_{\mathbb{R}_+^m}(\lambda_{k+1}) \\ \phantom{0 \le \lambda_{k+1} \perp \bar{D}_h \lambda_{k+1} + h\beta_k \ge 0 \Leftrightarrow} \Leftrightarrow \lambda_{k+1} \in (\bar{D}_h + \mathcal{N}_{\mathbb{R}_+^m})^{-1}(-h\beta_k), \end{array} $$

with $\bar{D}_h \triangleq h\hat{D} + h^2C_1B_1 + B_3^\top X_3B_3$, $\beta_k \triangleq C_1(I_{n_1} + hA_1)x_{1,k} - \frac{1}{h}C_3x_{2,k}$. It is noteworthy that $z_{k+1}$ is eliminated and does not appear in $\beta_k$. The inclusion (3.2) is a VI as (B.2) with $\varphi(\cdot) = \Psi_{\mathbb{R}_+^m}(\cdot)$. Obviously the well-posedness of the OSNSP depends on the properties of $\bar{D}_h$. From the passivity we know that $\hat{D} \ge 0$ and $B_3^\top X_3 B_3 \ge 0$. Using Theorem A.1 and Corollary A.2 it is possible to characterize the set of admissible $hC_1B_1$ that do not destroy the positive definiteness of $\bar{D}_h$ in case $\hat{D} > 0$.
**Remark 3.1.** We do not investigate here whether or not the discrete time system (3.1) is dissipative. Such a study is made in [55] for the case of LCS discretized with a $\theta$-method.

As alluded to in Remark 2.2, the cone $\mathbb{R}_+^m$ can be replaced in (2.1), (2.2) and (3.1) by any closed non empty convex cone $K \subset \mathbb{R}^m$ and its dual cone $K^*$. Then the complementarity conditions become $K^* \ni \lambda(t) \perp w(t) = C_1x_1(t) + C_2x_2(t) + C_3z(t) + D\lambda(t) \in K$, yielding the inclusion

$$ (3.3) \qquad \begin{aligned} & \bar{D}_h \lambda_{k+1} + h\beta_k \in -\mathcal{N}_{K^*}(\lambda_{k+1}) \\ & \Leftrightarrow \lambda_{k+1} \in (\bar{D}_h + \mathcal{N}_{K^*})^{-1}(-h\beta_k), \end{aligned} $$

which is a linear cone complementarity problem (LCCP), instead of (3.2). One can also consider time-varying terms in (2.2), in both the dynamics and the complementarity conditions, and adapt $\beta_k$ accordingly. Thus in the next proposition a more general complementarity framework is considered where the positive orthant is replaced by a closed convex cone.

**Proposition 3.2.** Let $h > 0$ and $k \ge 0$ be given. Assume that $\bar{D}_h \ge 0$, not necessarily symmetric, and let $K \subset \mathbb{R}^m$ be a non empty closed convex cone (where the complementarity conditions in (2.2) are extended to $K^* \ni \lambda(t) \perp w(t) = C_1x_1(t) + C_2x_2(t) + C_3z(t) + D\lambda(t) \in K$). Let $S = K^* \cap \text{Ker}(\bar{D}_h + \bar{D}_h^\top) \cap \{w \in \mathbb{R}^m \mid \bar{D}_h w \in K\}$. Then:

1. If $S = \{0\}$, for each $h\beta_k$ the LCCP in (3.3) has at least a solution.

2. If $S \neq \{0\}$, and there exists $\xi_0 \in \mathbb{R}^m$ such that $\langle h\beta_k - \bar{D}_h^\top \xi_0, v \rangle + \Psi_{K^*}(v) > 0$ for all $v \neq 0$, $v \in S$, then the LCCP in (3.3) has at least a solution.

3. If $\lambda_{k+1}^1$ and $\lambda_{k+1}^2$ are two solutions of the LCCP, then $\lambda_{k+1}^1 - \lambda_{k+1}^2 \in \text{Ker}(\bar{D}_h + \bar{D}_h^\top)$.

*Proof.* Let us consider Proposition B.2. Comparing (3.3) and (B.2) one has $\varphi(\cdot) = \Psi_{K^*}(\cdot)$, $M = \bar{D}_h$, $q = h\beta_k$. Using Proposition B.1 c) d), it follows that $(\text{dom}(\varphi))_\infty = K^*$. The second set in $S$ is obvious. The third set necessitates the calculation of $(\text{dom}(\varphi_\infty))^*$, see (B.3). Again the result of the calculation follows from Proposition B.1 c) and d). Then the proof follows from Proposition B.2 a), b) and c). $\square$

Once the existence of $\lambda_{k+1}$ is established for any $k \ge 0$, then the existence of $x_{1,k+1}$, $x_{2,k+1}$ and $z_{k+1}$ can be inferred from (3.1) for any $h > 0$.

**Remark 3.3.** If $K = \mathbb{R}_+^m$ then other results like [43, Theorem 3.8.6] may be used to study the LCP (3.2), using for instance the copositivity of $\bar{D}_h$.

Let us state a consequence of Proposition 3.2. The definition of a strongly SPR (SSPR) system is recalled in Appendix D.

**Corollary 3.4.** Let the quintuple $(P, A, B, C, D)$ be SSPR. Then, the LCCP in (3.3) has a unique solution for any $h\beta_k$, for $h > 0$ small enough.

*Proof.* This follows from [28, Theorem 5.4] ($\Rightarrow \hat{D} > 0$), the fact that $B_3^\top X_3 B_3 \ge 0$, Corollary A.2 ($\Rightarrow \bar{D}_h > 0$ for $0 < h < h^* \triangleq \frac{1}{\|C_1B_1\|_2\|(\bar{D}_h+\bar{D}_h^\top)^{-1}\|_2}$) $\Rightarrow S = \{0\}$), and items 1 and 3 in Proposition 3.2. $\square$

It is noteworthy that the SSPRness is only sufficient to guarantee that $\hat{D} > 0$, which is itself a sufficient condition for the existence of solutions to the LCCP (3.3).

# 4 BOUNDEDNESS AND CONVERGENCE ANALYSIS

From now on, we assume that the complementarity conditions are given as in (2.2) and (3.1), i.e., $K = K^* = \mathbb{R}_+^m$, see Remark 4.2. The next result characterizes least-norm solutions of the LCP.
**Corollary 4.1.** Assume that the conditions for existence of a solution in Proposition 3.2 are satisfied. Then the set of solutions to the LCP in (3.2) is nonempty convex polyhedral and contains a unique least-norm element $\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k})$ such that $\|\lambda_{k+1}^{\min}\| \le \|\lambda_{k+1}\|$ for any other solution $\lambda_{k+1}$. Moreover there exists a constant $\eta > 0$ depending only on $\bar{D}_h(=h\hat{D} + h^2C_1B_1 + B_3^\top X_3B_3)$ such that $\|\lambda_{k+1}^{\min}\| \le \eta h\|\beta_k\|$.

*Proof.* Follows from [43, Theorem 3.1.7] and [57, Lemma 1]. □

The constant $\eta$ generally depends on $h$ since it depends on $\bar{D}_h$, and we may denote it as $\eta(h)$ (as the examples in Section 6 show, it may also be constant in some cases). This result certainly is important because it will subsequently allow us to get upper bounds on the iterates of the implicit scheme (3.1). Notice that if $\bar{D}^h = 0$ then $\beta_{k,i} > 0 \Rightarrow \lambda_{k+1,i} = 0$, and $\beta_{k,i} = 0 \Rightarrow \lambda_{k+1,i} \in \mathbb{R}_+$, hence $\lambda_{k+1}^{\min} = 0$ so that $\eta = 0$. If $h = 0$ then $\bar{D}^0 = B_3^\top X_3B_3$, $\beta_0 = -C_3x_{2,k}$, and $\eta = \eta_0$ such that $\|\lambda_{k+1}^{\min}\| \le \eta_0\|C_3x_{2,k}\|$ is independent of $h$. If the LCP($\bar{D}_h, h\beta_k$) enjoys the stability of solutions property [43, Theorem 7.3.12], then a small perturbation of $\bar{D}_h$ implies a small variation of the solution, hence of $\lambda_{k+1}^{\min}$, and consequently of the constant $\eta$. Assume that $C_3 = 0$ (this is Assumption 2 below). Then solutions of the LCP($\bar{D}_h, h\beta_k$) are also solutions of LCP($\hat{D} + hC_1B_1, C_1(I_{n_1} + hA_1)x_{1,k}$). If $\hat{D} + hC_1B_1$ is a P-matrix (which occurs if $\hat{D}$ is a P-matrix and $h > 0$ is small enough, see Corollary A.2), then it follows from [3, Theorem B.3] [34] that the unique solution of LCP($\hat{D} + hC_1B_1, C_1(I_{n_1} + hA_1)x_{1,k}$) (which is also the least-norm solution) satisfies $\|\lambda_{k+1}^{\min}\| \le \eta \|C_1(I_{n_1} + hA_1)x_{1,k}\|$ for some constant $\eta > 0$.

**Remark 4.2.** As long as a minimal norm element exists and satisfies the property stated in the corollary, one may consider an LCCP as in (3.3) instead of an LCP (3.2).

## 4.1 CONVERGENCE ANALYSIS

From (3.2) and under the conditions of Corollary 4.1, the mapping $\beta_k \mapsto \lambda_{k+1}(-h\beta_k) = \lambda_{k+1}^{\min}$ with $\lambda_{k+1}^{\min}$ the least-norm element of the convex polyhedral set $(\bar{D}_h + N_{\mathbb{R}_+^m})^{-1}(-h\beta_k)$, is a single-valued mapping. The implicit scheme (3.1) may be modified as follows (which may be named the “minimal norm” implicit scheme):

$$ (4.1) \qquad \left\{ \begin{array}{l} x_{1,k+1} = (I_{n_1} + hA_1)x_{1,k} + hB_1\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) \\ x_{2,k+1} = x_{2,k} + hz_{k+1} + hB_2\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) \\ \qquad 0 = x_{2,k+1} + B_3\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) \\ \qquad 0 \le \lambda_{k+1}^{\min} \perp w_{k+1} = C_1x_{1,k+1} + C_2x_{2,k+1} + C_3z_{k+1} + D\lambda_{k+1}^{\min} \ge 0 \\ \|\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k})\| \le \eta\|h\beta_k\|, \quad h\beta_k = hC_1(x_{1,k} + hA_1x_{1,k}) - C_3x_{2,k}, \quad \eta > 0, \end{array} \right. $$

for all $k \in \{0, \dots, N-1\}$. The above results explore the well-posedness of the OSNSP, for fixed $h > 0$. Our goal is now to analyse the boundedness of the iterates on the integration interval $[0, T]$, $T > 0$. For all $k \ge 1$ we have $x_{2,k} = -B_3\lambda_k^{\min}$, and $hz_{k+1} = -B_3\lambda_{k+1}^{\min} - x_{2,k} - hB_2\lambda_{k+1}^{\min}$, thus it is inferred that $x_{2,k+1} = -B_3(\lambda_{k+1}^{\min} + \lambda_k^{\min}) - x_{2,k}$. Consequently (4.1) is rewritten equivalently as:

$$ (4.2) \qquad \left\{ \begin{array}{l} x_{1,k+1} = (I_{n_1} + hA_1)x_{1,k} + hB_1\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) \\ x_{2,k+1} = -x_{2,k} - B_3(\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) + \lambda_k^{\min}(x_{1,k-1}, x_{2,k-1})) \\ \qquad 0 \le \lambda_{k+1}^{\min} \perp w_{k+1} = C_1(I_{n_1} + hA_1)x_{1,k} - (C_2 + \frac{1}{h}C_3)x_{2,k} + \frac{1}{h}\bar{D}_h\lambda_{k+1}^{\min} - C_2B_3\lambda_k^{\min} \ge 0 \\ \qquad 0 = x_{2,k+1} + B_3\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}), \\ \qquad hz_{k+1} = -B_3\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) - x_{2,k} - hB_2\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k}) \\ \|\lambda_{k+1}^{\min}(x_{1,k}, x_{2,k})\| \le \eta\|h\beta_k\|, \quad h\beta_k = hC_1(I_{n_1} + hA_1)x_{1,k} - C_3x_{2,k}, \quad \eta > 0, \end{array} \right. $$

for all $k \in \{1, N-1\}$, given initial data $x_{1,0}$ and $x_{2,0}$. Starting from (4.2) the next result is obtained:
Lemma 4.3. Let $\alpha_{1,k} \triangleq \|x_{1,k}\|$, $\alpha_{2,k} \triangleq \|x_{2,k}\|$, $\alpha_k = (\alpha_{1,k} \ \alpha_{2,k})^{\top}$, then:

$$ (4.3) \qquad \left\{ \begin{array}{l} \alpha_{k+1} \le (M_i M + M_{i-1} U) \alpha_{k-i} + M_i U \alpha_{k-i-1} \\ M_{i+1} = M_i M + M_{i-1} U, \quad M_0 = I_2, M_{-1} = 0, \ 0 \le i \le k-1, \ k \ge 1, \end{array} \right. $$

where $M = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$, $U = \begin{pmatrix} 0 & 0 \\ e & f \end{pmatrix}$, $a = \|I_{n_1} + hA_1\| (1 + h^2\eta \|B_1\| \|C_1\|)$, $b = h\eta \|B_1\| \|C_3\|$, $c = h\eta \|B_3\| \|C_1\| \|I_{n_1} + hA_1\|$, $d = (1 + \eta \|B_3\| \|C_3\|)$, $e = c$, $f = \eta \|B_3\| \|C_3\|$. In particular

$$ (4.4) \qquad \alpha_{k+1} \le (M_{k-1}M + M_{k-2}U)\alpha_1 + M_{k-1}U\alpha_0 $$

for all $k \ge 1$.

*Proof.* The proof starts from (4.2). The following upper bounds can be derived:

$$ (4.5a) \qquad \|\lambda_{k+1}^{\min}\| \le \eta h \|C_1\| \|I_{n_1} + hA_1\| \|x_{1,k}\| + \eta \|C_3\| \|x_{2,k}\| $$

$$ (4.5b) \qquad \begin{aligned} \|x_{1,k+1}\| &= \|I_{n_1} + hA_1\| (1 + h^2\eta \|B_1\| \|C_1\|) \|x_{1,k}\| + h\eta \|B_1\| \|C_3\| \|x_{2,k}\| \\ &= a\|x_{1,k}\| + b\|x_{2,k}\| \end{aligned} $$

$$ (4.5c) \qquad \begin{aligned} \|x_{2,k+1}\| &\le (1 + \eta \|B_3\| \|C_3\|) \|x_{2,k}\| + \eta \|B_3\| \|C_3\| \|x_{2,k-1}\| \\ &\quad + h\eta \|B_3\| \|C_1\| \|I_{n_1} + hA_1\| (\|x_{1,k}\| + \|x_{1,k-1}\|) \\ &= d\|x_{2,k}\| + f\|x_{2,k-1}\| + c(\|x_{1,k}\| + \|x_{1,k-1}\|). \end{aligned} $$

One infers that:

$$ (4.6) \qquad \alpha_{k+1} \le M\alpha_k + U\alpha_{k-1} $$

The upper bound in (4.3) therefore holds for $i = 0$. Let us prove by induction that it holds for all $0 \le i \le k - 1$. Let it hold for $i > 0$, we have $\alpha_{k-i} \le M\alpha_{k-i-1} + U\alpha_{k-i-2}$. Thus $\alpha_{k+1} \le (M_i M + M_{i-1}U)M\alpha_{k-i-1} + (M_i M + M_{i-1}U)U\alpha_{k-i-2} + M_i U\alpha_{k-i-1} = (M_{i+1}M + M_i U)\alpha_{k-i-1} + M_{i+1}U\alpha_{k-i-2}$. The result is proved. $\square$

Our goal is now to examine (4.4). Let us come back to the matrices $M$ and $U$ in Lemma 4.3. It is clear from (4.4) that the boundedness of the sequence $\{\alpha_k\}_{k \ge 2}$ depends on the boundedness of $M_k$, $k \ge 1$. Apart from $d$ and $f$, all the other terms in $M$ and $U$ are proportional to $h$, $\eta h$ or to $h_1 \triangleq \|I_{n_1} + hA_1\|$. Moreover $d > 1$ whenever $B_3 \ne 0$ and $C_3 \ne 0$. The term in position (2, 2) of $M^n$ is equal to $d^n + O(h_1 h^2)$. Before stating the next result, let us introduce the following assumptions:

**Assumption 1.** The matrix $A_1$ is such that $h_1 \le 1$ for $h > 0$ small enough.

Let us notice that the induced norm $\| \cdot \|_{2,2} = \sigma_{\max}(\cdot)$ [16, Proposition 9.4.9] is compatible with the vector norm $\| \cdot \|_2 = \| \cdot \|_{[16, pp.607-609]}$, and can be used instead of the Frobenius norm to compute $h_1$. The next assumption is crucial in the body of this section.

**Assumption 2.** One has $B_{3,\bullet i} \in \text{Ker}(X_3)$ for all $1 \le i \le m (\implies C_3 = 0)$.

Let us remind that $X_3$ is uniquely defined under minimality and positive-realness [28, Theorem 5.4]. If Assumption 2 holds then $d = 1$ and $b = f = 0$. This assumption means that the variable $z$ does not enter the complementarity conditions, see (2.2). Hence in a sense the complementarity part of the LDCS and the algebraic part are decoupled, and $\beta_k$ no longer depends on $\frac{1}{h}C_3x_{2,k}$. The terms in (4.5) involve products with $\eta(h)$, and it is important to guarantee that they are bounded, when Assumption 2 holds.
**Assumption 3.** The constant $\eta = \eta(h)$ is such that the products $\eta(h)h$ can be made arbitrarily small for small enough $h \ge 0$.

Notice that this assumption is verified in the examples of Section 6, where $\eta$ is a constant (but these examples do not verify Assumption 2, however). The matrices $M_k$ then involve only products of powers of the terms $a, b, c, d = 1, e$. These products are bounded for bounded $k$, and provided that both $h$ and $h_1$ are small enough so that $a \le 1, b \le 1, c \le 1, d = 1, e \le 1$, they remain bounded as $k \to +\infty$. If the inequalities are strict (apart from $d = 1$), the first rows in both $M_{k-1}M + M_{k-2}U$ and $M_{k-1}U$ converge to zero, while the second rows converge to $(0, 1)$ (the 1 being the consequence of a term $d^n$ in matrices $M^n$). The next proposition is therefore proved:

**Proposition 4.4.** Assume that: a) the conditions of existence of solutions in Proposition 3.2 are verified, b) Assumptions 1 and 3 hold true, and c) both $h$ and $h_1$ are small enough, so that $a \le 1, b \le 1, c \le 1, d = 1, e \le 1$. Then the sequences $\{\alpha_{1,k}\}_k$ and $\{\alpha_{2,k}\}_k$ are uniformly bounded in $k$. If $A_1 < 0$, and $h > 0$ is small enough such that $a < 1, c < 1, e < 1$, then $\alpha_{1,k} \to 0$ as $k \to +\infty$.

Recall that $\alpha_{1,k} = \|x_{1,k}\|$ and $\alpha_{2,k} = \|x_{2,k}\|$. The last assertion of the proposition follows from the material in Appendix A. Then we have the next corollary.

**Corollary 4.5.** Let the conditions a), b), c) of Proposition 4.4 hold. Then $\|\lambda_{k+1}^{\min}\|$ is uniformly bounded in $k$ and $h > 0$, and it converges to zero as $k \to +\infty$ if $\alpha_{1,k} \to 0$ as $k \to +\infty$. Therefore the sequence $\{\frac{x_{1,k+1}-x_{1,k}}{h}\}_k$ is uniformly bounded also in $k$ and $h > 0$.

**Proof.** The boundedness of $\|\lambda_{k+1}^{\min}\|$ and its convergence follow from (4.5) (a), the boundedness of $\|x_{1,k}\|$ and Assumption 2. From the first line in (4.2), one has $\frac{x_{1,k+1}-x_{1,k}}{h} = A_1x_{1,k} + B_1\lambda_{k+1}^{\min}$, hence the result is proved. $\square$

Let us now study the convergence of the piecewise-linear and step-function approximations defined as:

$$ (4.7) \qquad \begin{cases} x_1^h(t) = x_{1,k+1} + \frac{t_{k+1} - t}{h}(x_{1,k} - x_{1,k+1}) \\ \dot{x}_1^h(t) = \frac{x_{1,k+1} - x_{1,k}}{h} \\ x_{1,h}^*(t) = x_{1,k+1} \end{cases} \qquad \text{for all } t \in [t_k, t_{k+1}) $$

and similarly for the other variables.

**Corollary 4.6.** There exists subsequences denoted as $\{x_1^h(\cdot)\}_{h \ge 0}$ which converge uniformly to a limit $x_1(\cdot)$ in $C^0([0, T]; \mathbb{R}^{n_1})$, such that almost everywhere one has $\dot{x}_1(t) = A_1x_1(t) + \lambda_1(t)$, for some function $\lambda_1(\cdot)$. The sequence $\{x_1^{h^*}(\cdot)\}_{h \ge 0}$ converges to $x_1(\cdot)$ strongly in $L_2([0, T]; \mathbb{R}^{n_1})$. Moreover the sequences $\{x_{2,h}^{h^*}(\cdot)\}_{h \ge 0}$ and $\{\lambda_h^{\min,h^*}(\cdot)\}_{h \ge 0}$ converge weakly$^*$ in $L_\infty([0, T]; \mathbb{R}^l)$, $l = n_2$ and $l = m$, respectively, to limits $x_2(\cdot)$ and $\lambda(\cdot)$ which satisfy the continuous-time constraints.

**Proof.** Both $\{x_1^h(\cdot)\}_h$ and $\{\dot{x}_1^h(\cdot)\}_h$ are uniformly bounded on $[0, T]$. The sequence of $C^0([0, T]; \mathbb{R}^{n_1})$ functions $\{x_1^h(\cdot)\}_h$ is thus equicontinuous. By the Arzela-Ascoli Theorem this sequence stays in a compact subset of $C^0([0, T]; \mathbb{R}^{n_1})$ and converges uniformly towards a continuous limit $x_1(\cdot) : [0, T] \to \mathbb{R}^{n_1}$ as $h \to 0, h > 0$. Let us notice that $x_1(t) \to 0$ as $t \to +\infty$. Indeed from the above results which are stated for $k \in \{0, \dots, N\}$ and $h > 0$, $k \to +\infty$ only if $T = Nh \to +\infty$. Let us now define the step functions $x_1^{h^*}(t) = x_{1,k+1}$ if $t \in [t_k, t_{k+1})$, $\lambda_h^{\min}(t) = \lambda_{k+1}^{\min}$ if $t \in [t_k, t_{k+1})$, and

$$ (4.8) \qquad \dot{x}_1^h(t) = A_1x_1^{h^*}(t) + B_1\lambda_h^{\min}(t), \quad \text{almost everywhere.} $$
In addition one has:

$$ (4.9) \quad \begin{aligned} \|x_1^h(t) - x_1^{h\star}(t)\|_{\mathcal{L}_2([0, T]; \mathbb{R}^{n_1})}^2 &= \sum_{k=0}^{N-1} \int_{t_k}^{t_{k+1}} (t_{k+1} - t)^2 \|\dot{x}_{1,k}\|^2 dt \le C^2 \sum_{k=0}^{N-1} \int_{t_k}^{t_{k+1}} (t_{k+1} - t)^2 dt \\ &\le C^2 \sum_{k=0}^{N-1} \frac{1}{3}(t_{k+1} - t_k)^3 = \frac{C^2 N h^3}{3} = \frac{C^2 T h^2}{3}, \end{aligned} $$

where $\dot{x}_{1,k} = \frac{x_{1,k+1} - x_{1,k}}{h}$, and $\|\dot{x}_{1,k}\| \le C$ for all $k \ge 0$ for some $C > 0$ from Corollary 4.5. It follows that $x_1^{h\star}(\cdot) \to x_1(\cdot)$ strongly in $\mathcal{L}_2([0, T]; \mathbb{R}^{n_1})$ as $h \to 0$. On the other hand since $\dot{x}_1^h(\cdot) \in \mathcal{L}_\infty([0, T]; \mathbb{R}^{n_1})$ it follows from the Banach-Alaoglu-Bourbaki Theorem [18, Theorem III.15] that $\dot{x}_1^h(\cdot) \to \dot{x}_{1,\lim}(\cdot)$ in the weak\star topology as $h \to 0$. Since the considered domain $[0, T]$ is bounded, one has $\mathcal{L}_1([0, T]; \mathbb{R}^{n_1}) \subseteq \mathcal{L}_2([0, T]; \mathbb{R}^{n_1})$, thus $\dot{x}_1^h(\cdot)$ converges weakly in $\mathcal{L}_2([0, T]; \mathbb{R}^{n_1})$. Moreover since $x_1^h(t) = x_1^h(0) + \int_0^t \dot{x}_1^h(t)dt$, one has $\dot{x}_{1,\lim}(\cdot) = \dot{x}_1(\cdot)$ almost everywhere. Therefore we deduce using (4.8) that $(B_1\lambda_h^{\min,\star})(\cdot)$ converges weakly in $\mathcal{L}_2([0, T]; \mathbb{R}^{n_1})$ as $h \to 0$ towards some limit $\lambda_1(\cdot)$, and we have $\dot{x}_1(t) = A_1x_1(t) + \lambda_1(t)$ almost everywhere on $[0, T]$.

Let us consider now the step-functions $x_{2,h}^{\star}(\cdot)$ and $\lambda_h^{\min,\star}(\cdot)$ defined as $x_{2,h}^{\star}(t) = x_{2,k+1}$ for all $t \in [t_k, t_{k+1})$, and similarly for the multiplier. It follows from Proposition 4.4 and from Corollary 4.5 that both sequences $\{x_{2,h}^{\star}(\cdot)\}_{h \ge 0}$ and $\{\lambda_h^{\min,\star}(\cdot)\}_{h \ge 0}$ are uniformly bounded. Hence from the Banach-Alaoglu-Bourbaki Theorem, there exist subsequences of both which converge weakly\star to limits $x_2 \in \mathcal{L}_\infty([0, T]; \mathbb{R}^{n_2})$ and $\lambda \in \mathcal{L}_\infty([0, T]; \mathbb{R}^m)$, respectively. We infer that $0 = x_{2,h}^{\star}(t) + B_3\lambda_h^{\star}(t) \to x_2(t) + B_3\lambda(t)$ for almost all $t \ge 0$ in the weak\star sense. Next, using the second line in (4.1) it follows that the sequence $\{hz_k\}_{h \ge 0}$ is also uniformly bounded. Denoting $\sigma_k = hz_k$, it follows that the sequence $\{\sigma_h^\star\}_{h \ge 0}$ converges weakly\star to a limit denoted as $\sigma(t) \in \mathcal{L}_\infty([0, T]; \mathbb{R}^{n_2})$. Moreover one has $\sigma(t) = 0$ almost everywhere since $x_{2,h}^{\star}(t) - x_{2,h}^{\star}(t - t_k) - hB_2\lambda_h^{\min,\star}(t) = \sigma_h^\star(t)$ almost everywhere in $[0, T]$. Finally, the complementarity conditions in (4.1) are equivalently rewritten as $0 \le \lambda_h^\star(t) \perp C_1x_{1,h}^\star(t) + C_2x_{2,h}^\star(t) + D\lambda_h^\star(t) \ge 0$ for all $t \in [t_k, t_{k+1})$, since $C_3 = 0$ by Assumption 2. It is deduced that in the limit as $h \to 0$ one obtains $0 \le \lambda(t) \perp C_1x_1(t) + C_2x_2(t) + D\lambda(t) \ge 0$. $\square$

## 4.2 STABILITY ANALYSIS VIA A LYAPUNOV FUNCTION

To complete the foregoing section, it is natural to analyze how the Lyapunov function $V(x_1, x_2)$ for the continuous-time system in (2.2) (see also Appendix D), could be used in the discrete-time context for (4.2). Indeed implicit Euler methods are known to be able to mimic the continuous-time passivity and stability properties, in some cases [55]. It is noteworthy that Assumption 2 implies that $X_3x_{2,k} = -X_3B_3\lambda_k = 0$ and $X_3x_2(t) = -X_3B_3\lambda(t) = 0$ for all $k$ and $t$. Therefore $V(x_1, x_2) = \frac{1}{2}x_1^\top X_1x_1$, so that this Lyapunov function is not a positive definite function of the whole state. For the sake of completeness let us nevertheless make the calculations in the general setting. We set $V_k \stackrel{\Delta}{=} \frac{1}{2}x_{1,k}^\top X_1x_{1,k} + x_{2,k}^\top X_3x_{2,k}$. It follows that:

$$ (4.10) \quad \begin{aligned} V_{k+1} - V_k &= \frac{1}{2}(x_{1,k+1} - x_{1,k})^\top X_1(x_{1,k+1} + x_{1,k}) + (x_{2,k+1} - x_{2,k})^\top X_3(x_{2,k+1} + x_{2,k}) \\ &= \frac{1}{2}(hA_1x_{1,k} + hB_1\lambda_{k+1})^\top X_1(2x_{1,k} + hA_1x_{1,k} + hB_1\lambda_{k+1}) \\ &\quad + (hz_{k+1} + hB_2\lambda_{k+1})^\top X_3(2x_{2,k} + hz_{k+1} + hB_2\lambda_{k+1}) \\ &= h(x_{1,k}^\top A_1^\top X_1x_{1,k} + \lambda_{k+1}^\top B_1^\top X_1x_{1,k}) + h^2(\frac{1}{2}x_{1,k}^\top A_1^\top X_1A_1x_{1,k} + \frac{1}{2}\lambda_{k+1}^\top B_1^\top X_1B_1\lambda_{k+1}) \\ &\quad + h^2x_{1,k}^\top A_1^\top X_1B_1\lambda_{k+1} + hz_{k+1}^\top C_3\lambda_{k+1} + h\lambda_{k+1}^\top C_3B_2\lambda_{k+1} + hz_{k+1}^\top X_3x_{2,k} \\ &\quad + h\lambda_{k+1}^\top B_2^\top X_3x_{2,k} + h\lambda_{k+1}^\top w_{k+1} - h\lambda_{k+1}^\top w_{k+1} \end{aligned} $$
Since $\lambda_{\max}(Q_3) < 0$ (which is the reason why the second inequality in (4.12) holds), one has $1 + \frac{h}{2} \frac{\lambda_{\max}(Q_3)}{\lambda_{\max}(X_1)} < 1$ for all $h > 0$. Therefore all $V_k$ are uniformly bounded for any bounded $V_0$, and $V_k \to 0$ as $k \to +\infty$ for any $h \in (0, h_{\max})$. Let us set $\lambda_{k+1} = \lambda_{k+1}^{\min}$. Since the dependence of $\lambda_{k+1}^{\min}$ in $x_{2,k}$ through $\beta_k$ is $C_3x_{2,k}$ (see (4.2)), it follows from Assumptions 2 and 3 that $\lambda_{k+1}^{\min} = \lambda_{k+1}^{\min}(x_{1,k})$, and $\lambda_{k+1}^{\min}(x_{1,k}) \to 0$ as $k \to +\infty$. Using the third line in (4.2), we infer that the iterates $x_{2,k}$ (which are uniformly bounded from Proposition 4.4) converge to zero asymptotically also. From the results in Section 4.1 and the above Lyapunov analysis, we have therefore proved the following:

**Lemma 4.7.** Let Assumptions 2 and 3 hold, and the quintuple $(P, A, B, C, D)$ be strongly SPR and minimal. Then there exists $h_{\max} > 0$ depending on $A_1, C_1, X_1, X_3, B_1, B_2$, such that for all $h \in (0, h_{\max})$, and for all bounded initial condition $x_{1,0}$, the sequences $\{x_{1,k}\}_k, \{x_{2,k}\}_k, \{\lambda_k^{\min}\}_k$ are uniformly bounded and $x_{1,k}, x_{2,k}, \lambda_{k+1}^{\min} \to 0$ as $k \to +\infty$.

**Remark 4.8.** Concerning Assumption 1, let us examine the constants $a, b, c, d, f$ in Lemma 4.3, and the inequalities in (4.5), as well as the above Lyapunov analysis. The term $h_1$ is present in both $a$ and $c$. The inequality $a \le 1$ implies $h_1 < 1$, but $c \le 1$ does not imply that $h_1 < 1$. The constant $a$ is present in (4.5) (b), however the Lyapunov analysis guarantees the boundedness of $x_{1,k}$ without this assumption (the analysis relies on the passivity of $(A_1, B_1, C_1, \hat{D})$ and $h \ge 0$ small enough). But, Lemma 4.7 relies on $Q_1 < 0$, which is a strong passivity condition and may be restrictive in applications. There exists some cases where SSPRness can be relaxed, as shown next.

Let $(A_1, B_1, C_1, \hat{D})$ be strictly state passive (⇒ $X_1A_1 + A_1^\top X_1 = -\mu X_1, \mu > 0, X_1 = X_1^\top > 0, \mu$ is called the dissipation constant) and $\hat{D} + \hat{D}^\top = 0 \Rightarrow X_1B_1 = C_1^\top$, then $Q_2$ in (4.11) and $Q_1$ can be rewritten as:

$$ (4.13) \qquad Q_2 = \begin{pmatrix} A_1^\top X_1 A_1 & 0 \\ 0 & B_2^\top X_3 B_2 - \frac{3}{2} B_1^\top X_1 B_1 \end{pmatrix}, \qquad Q_1 = \begin{pmatrix} -\mu X_1 & 0 \\ 0 & 0 \end{pmatrix} $$

Thus $Q_1 + hQ_2 < 0$ if and only if $\mu X_1 - hA_1^\top X_1 A_1 > 0$ and $B_2^\top X_3 B_2 - \frac{3}{2} B_1^\top X_1 B_1 < 0$. Let us recall that if $(A_1, B_1, C_1, \hat{D})$ is minimal, then strict state passivity and SPRness are equivalent, see Appendix D.

**Lemma 4.9.** Assume that: 1) the quadruple $(A_1, B_1, C_1, \hat{D})$ is strictly state passive with $\hat{D} + \hat{D}^\top = 0$ and dissipation constant $\mu > 0$, 2) $B_2^\top X_3 B_2 - \frac{3}{2} B_1^\top X_1 B_1 < 0$. Let Assumptions 2 and 3 hold. Then there exists $h_{\max} > 0$ such that for all $h \in [0, h_{\max})$, the sequences $\{x_{1,k}\}_k, \{x_{2,k}\}_k$ and $\{\lambda_k^{\min}\}_k$ are uniformly bounded.

**Proof.** Denote $\tilde{X}_1 = \mu X_1 - hA_1^\top X_1 A_1$, then $\tilde{X}_1 > 0$ for $h \ge 0$ small enough, using Corollary A.2. Using 2) and Assumption 2, and after few calculations, there exists $h_{\max} > 0$ such that:

$$ (4.14) \qquad V_{k+1} \le \left( 1 - \frac{h}{2} \frac{\lambda_{\min}(\tilde{X}_1)}{\lambda_{\max}(X_1)} \right) V_k. $$

for all $h \in [0, h_{\max})$. Thus $V_k \to 0$ as $k \to +\infty$, and $x_{1,k} \to 0$ as $k \to +\infty$. It can be inferred that $\|x_{1,k}\| \le \sqrt{\frac{\lambda_{\max}(X_1)}{\lambda_{\min}(X_1)}} \|x_{1,0}\| \alpha^{\frac{k}{2}}$ with $\alpha < 1$. From (4.5) (a) and Assumption 3 it follows that $\|\lambda_{k+1}^{\min}\| \to 0$ as $k \to +\infty$, and $\{\lambda_k^{\min}\}_k$ is uniformly bounded. From (4.1) it is inferred that $x_{2,k} \to 0$ as $k \to +\infty$ and $\{x_{2,k}\}_k$ is uniformly bounded. $\square$

Similar convergence results as in Corollary 4.6 can then be deduced, and are not developed here for the sake of briefness. One sees from the above and Corollary A.2 that the degree of strict SPRness (represented by the constant $\mu$) is related to the value of $h_{\max}$: the greater $\mu$, the greater $h_{\max}$. Finally it is noteworthy that the strict state passivity in item 1) of Lemma 4.9 can again be relaxed, provided

$$ \text{that the matrix } A_1^\top X_1 + X_1 A_1 + h A_1^\top X_1 A_1 < 0. \text{ Consider as an example } A_1 = \begin{pmatrix} -1 & 0 \\ 0 & 0 \end{pmatrix}, X_1 = I_2, \text{ then} \\ A_1^\top X_1 + X_1 A_1 + h A_1^\top X_1 A_1 = \begin{pmatrix} -2 + h & 0 \\ 0 & 0 \end{pmatrix} < 0 \text{ for all } 0 \le h < 2. $$
# 5 ANOTHER APPROACH FOR PASSIVE SYSTEMS

The LDCS in (2.2) are suitably transformed in [32] to enable their analysis. The transformation may be seen as a kind of relative degree reduction, involving the multiplier derivative. As will become clear later, this new dynamics allows us to somewhat relax the stringent Assumption 2 (see Assumption 4 c) below). More accurately, a close relationship is shown in [32, Lemma 3, Theorem 3] between (2.2) and the following system:

$$ (5.1) \qquad \left\{ \begin{array}{ll} (a) \begin{pmatrix} I_{n_1} & 0 \\ 0 & B_3^\top X_3 B_3 \end{pmatrix} \begin{pmatrix} \dot{x}_1(t) \\ \dot{\lambda}(t) \end{pmatrix} = \begin{pmatrix} A_1 & B_1 \\ -C_1 & -\hat{D} \end{pmatrix} \begin{pmatrix} x_1(t) \\ \lambda(t) \end{pmatrix} + \begin{pmatrix} 0 \\ I_m \end{pmatrix} w(t) & \\ (b) \quad 0 \le \lambda(t) \perp w(t) = C_1 x_1(t) + C_2 x_2(t) + C_3 z(t) + D \lambda(t) \ge 0 & \\ (c) \quad x_2(t) = -B_3 \lambda(t), \quad z(t) = -B_2 \lambda(t) - B_3 \dot{\lambda}(t). & \end{array} \right. $$

It is noteworthy that this dynamics is obtained from the original one by treating the multiplier $\lambda$ as a state variable. The “equivalence” between both systems is to be understood first as follows [32, Lemma 3]:

$$ \{(x_1, \lambda, w) \mid (x_1, x_2, \lambda, w) \text{ is a solution of (2.2) } (a, b, c)\} = \{(x_1, \lambda, w) \mid (x_1, \lambda, w) \text{ is a solution of (5.1) } (a) \text{ and } B_3\lambda \text{ is absolutely continuous}\}. $$

By solution of (2.2) it is meant that $(x_1, x_2)$ is absolutely continuous and the dynamics holds Lebesgue almost everywhere [32, Definition 3]. The initial conditions are chosen compatible with the constraints. It is noteworthy that the variables $z$ and $\lambda$ are not shown to be absolutely continuous, and examples in Sections 6.1 and 6.3 will illustrate this fact. One has $w(t) \in -N_{\mathbb{R}_+^m}(\lambda(t))$, therefore one obtains the differential inclusion (DI):

$$ (5.2) \qquad \left( \begin{array}{cc} I_{n_1} & 0 \\ 0 & B_3^\top X_3 B_3 \end{array} \right) \left( \begin{array}{c} \dot{x}_1(t) \\ \dot{\lambda}(t) \end{array} \right) - \left( \begin{array}{cc} A_1 & B_1 \\ -C_1 & -\hat{D} \end{array} \right) \left( \begin{array}{c} x_1(t) \\ \lambda(t) \end{array} \right) \in - \left( \begin{array}{c} 0 \\ I_m \end{array} \right) N_{\mathbb{R}_+^m}(\lambda(t)) \\ = - \left( \begin{array}{c} 0 \\ I_m \end{array} \right) N_{\mathbb{R}_+^m} \left( \begin{pmatrix} 0 & I_m \end{pmatrix} \begin{pmatrix} x_1(t) \\ \lambda(t) \end{pmatrix} \right) \\ = -N_\Phi \left( \begin{pmatrix} x_1(t) \\ \lambda(t) \end{pmatrix} \right) $$

with $\Phi \stackrel{\Delta}{=} \{(x_1, \lambda) \in \mathbb{R}^{n_1} \times \mathbb{R}^m \mid \lambda \ge 0\}$. If this DI has a solution $(x_1(\cdot), \lambda(\cdot))$, then $x_2(\cdot), z(\cdot)$, and $w(\cdot)$ can be obtained (provided indeed that $B_3\lambda(\cdot)$ is absolutely continuous). Multiplying the first line in (5.2) by $X_1$ and adding an exogenous excitation $G(t)$ to the left-hand side, one obtains

$$ (5.3) \qquad \underbrace{\begin{pmatrix} X_1 & 0 \\ 0 & B_3^\top X_3 B_3 \end{pmatrix}}_{\triangleq \mathcal{P}} \begin{pmatrix} \dot{x}_1(t) \\ \dot{\lambda}(t) \end{pmatrix} - \underbrace{\begin{pmatrix} X_1 A_1 & X_1 B_1 \\ -C_1 & -\hat{D} \end{pmatrix}}_{\triangleq \mathcal{A}} \begin{pmatrix} x_1(t) \\ \lambda(t) \end{pmatrix} + G(t) \in -N_\Phi \left( \begin{pmatrix} x_1(t) \\ \lambda(t) \end{pmatrix} \right) $$

which is similar to the form obtained in [32, Theorem 3]. Let us assume that an exogenous input $V(t) = (v_{(a)}(t), v_{(b)}(t), v_{(c)}(t), v_{(d)}(t))^{\top}$ is added to (2.2), then:

$$ (5.4) \qquad G(t) = \begin{pmatrix} v_{(a)}(t) \\ C_3(v_{(b)}(t) + \dot{v}_{(c)}(t)) + C_2v_{(c)}(t) - v_{(d)}(t) \end{pmatrix}. $$

Notice that $B_3^\top X_3 B_3 \ge 0$, therefore $\mathcal{P} = \mathcal{P}^\top \ge 0$. However $\mathcal{P}$ has a specific block-diagonal structure. Under strong SPRness and minimality, $\mathcal{A} < 0$, since $\mathcal{A} + \mathcal{A}^\top < 0$ from the KYP Lemma [25] (then $A_1$
is a Hurwitz matrix, $X_1 = X_1^\top > 0$ and $\hat{D} > 0$). The exogeneous term $G(t)$ stems from the exogeneous terms $E(t)$ and $F(t)$ in (2.1) and their time-derivatives (see (5.4) and (6.4) for an example).

In a second step the “equivalence” between both systems is to be understood as follows [32, Theorem 3]:

$$ \begin{aligned} & \{(x_1, \lambda) | (x_1, x_2, \lambda, w) \text{ is a solution of (2.2) (a,b,c,d)}\} \\ & = \{(x_1, \lambda) | (x_1, \lambda) \text{ is a solution of (5.3) and } B_3\lambda (= -x_2) \text{ is absolutely continuous}\}. \end{aligned} $$

As will become clear later, the results that are obtained using this approach, and the approach in the foregoing section, are not the same (especially concerning the convergence of the piecewise-linear discrete solutions).

## 5.1 THE TIME-DISCRETIZATION AND THE OSNSP

As alluded to in (5.3), we can enlarge the analysis by considering time-varying terms, so as to embed (5.3) into the DI: $P\xi(t) - A\xi(t) + G(t) \in -N_\Phi(\xi(t))$. The Euler implicit scheme associated with (5.3) is:

$$ (5.5) \qquad P\xi_{k+1} - \underbrace{(P + hA)}_{\Delta=\mathcal{P}_h} \xi_k + hG_k \in -N_\Phi(\xi_{k+1}), $$

It is noteworthy that this scheme is not equal to the implicit scheme in (3.1), as can be shown after few manipulations, see Appendix H. Though the difference between both approaches is not analysed further in this paper, it will be tackled numerically through an example in Section 6.3. In addition, a slightly more general convex set than in (5.5) is considered, for the sake of generality. It corresponds to changing the complementarity conditions in (5.1) (b) to cone complementarity conditions $K^* \ni \lambda(t) \perp w(t) \in K$, $K \subseteq \mathbb{R}^m$ a nonempty closed convex cone, similarly to what is done in Proposition 3.2. Thus $\Phi = \mathbb{R}^{n_1} \times K^*$.

**Proposition 5.1.** Let $h > 0$ and $k \ge 0$ be given. Let us consider the DI in (5.3) and its discretization in (5.5), with $\Phi = \mathbb{R}^{n_1} \times K^*$, with $K$ a nonempty closed convex cone. Let us denote $q_k \stackrel{\Delta}{=} -\mathcal{P}_h\xi_k + hG_k$. Let us consider the set:

$$ (5.6) \qquad S \stackrel{\Delta}{=} \{(x_1, \lambda) \in \mathbb{R}^{n_1} \times \mathbb{R}^m \mid x_1 \in \text{Ker}(X_1), \lambda \in \text{Ker}(B_3^\top X_3 B_3) \cap (B_3^\top X_3 B_3)^{-1}(K) \cap K^*\}. $$

Then:

1. If $S = \{0\}$, for each $q_k \in \mathbb{R}^{n_1 \times m}$ the problem $VI(P, q_k, \Psi_\Phi)$ has at least one solution.

2. Assume that $S \neq \{0\}$. If there exists $\xi^0 \in \text{dom}(\Phi)$ such that $\langle q_k - P^\top \xi^0, v \rangle > 0$ for all $v \in S, v \neq 0$, then the problem $VI(P, q_k, \Psi_\Phi)$ has at least one solution. In case $X_1 = X_1^\top$, one can take $\xi_0 = 0$.

3. If $\xi_{k+1}^1$ and $\xi_{k+1}^2$ are two solutions of $VI(P, q_k, \psi_\Phi)$, then $\xi_{k+1}^1 - \xi_{k+1}^2 \in \text{Ker}(P)$, hence $x_{1,k+1}$ is unique. Moreover $\langle q_k, \begin{pmatrix} 0 \\ \lambda_{k+1}^1 - \lambda_{k+1}^2 \end{pmatrix} \rangle = 0$.

4. If $X_1 > 0$ and $B_3^\top X_3 B_3 > 0$, then the problem $VI(P, q_k, \Psi_\Phi)$ always has a unique solution.

5. Assume that the quintuple $(P, A, B, C, D)$ is minimal and strongly SPR, so that the quadruple $(A_1, B_1, C_1, \hat{D})$ is minimal and strongly SPR. If $B_3^\top X_3 B_3 > 0$, then the problem $VI(P, q_k, \Psi_\Phi)$ always has a unique solution.

6. Let $K = K^* = \mathbb{R}_+^m$. Then $S = \{(x_1, \lambda) | x_1 \in \text{Ker}(X_1), \lambda \in \mathbb{R}_+^m \cap \text{Ker}(B_3^\top X_3 B_3)\}$.
Proof. From Proposition B.1 d) one has $(K^*)_{\infty} = K^*$ and $\Phi_{\infty} = \Phi$. Also $\varphi = \Psi_{\Phi}$ hence from Proposition B.1 c) and d), $\varphi_{\infty} = \Psi_{\mathbb{R}^{n_1 \times K^*}}$ and $(\mathrm{dom}(\varphi_{\infty}))^* = \Phi^* = \mathbb{R}^{n_1 \times K}$. Then one has $S = \Phi \cap \mathrm{Ker}(\mathcal{P}) \cap \{\xi \in \mathbb{R}^{n_1+m} | \mathcal{P}(\begin{pmatrix} x_1 \\ \lambda \end{pmatrix}) = \begin{pmatrix} 0 \\ \alpha \end{pmatrix}, \alpha \in K\} = \Phi \cap \mathrm{Ker}(\mathcal{P}) \cap \{\xi \in \mathbb{R}^{n_1+m} | x_1 \in \mathrm{Ker}(X_1), B_3^\top X_3 B_3 \lambda \in K, \lambda \in \mathrm{Ker}(B_3^\top X_3 B_3), \lambda \in K^*\}$. In item 2) one uses that $\Psi_{\Phi}(v) = 0$ when $v \in S$, and item b') in Proposition B.2. Let us now deal with item 4. In this case $\mathcal{P} = \mathcal{P}^\top > 0$ and the result follows from items 1 and 3. Let us define $\mathcal{R} = \mathcal{R}^\top > 0$, $\mathcal{R}^2 = \mathcal{P}$, and $\zeta_k = \mathcal{R}\zeta_k$. Then the VI in (5.5) is equivalently rewritten as:

$$ (5.7) \qquad \zeta_{k+1} - \mathcal{R}^{-1}\mathcal{P}_h\mathcal{R}^{-1}\zeta_k + \mathcal{R}^{-1}G_k \in -\mathcal{R}^{-1}\mathcal{N}_{\Phi}(\mathcal{R}^{-1}\zeta_{k+1}) = -\mathcal{N}_{\tilde{\Phi}}(\zeta_{k+1}), $$

where the last equality stems from the chain rule, with $\tilde{\Phi} = \{z \in \mathbb{R}^{n_1+m} | \mathcal{R}^{-1}z \in \Phi\}$. One infers from (5.7) that $\zeta_{k+1} = \mathrm{proj}[\tilde{\Phi}; \mathcal{R}^{-1}\mathcal{P}_h\mathcal{R}^{-1}\zeta_k - \mathcal{R}^{-1}G_k] \Rightarrow \zeta_{k+1} = \mathcal{R}^{-1}\mathrm{proj}[\tilde{\Phi}; \mathcal{R}^{-1}\mathcal{P}_h\zeta_k - \mathcal{R}^{-1}G_k]$. The first assertion in item 5 about strong SPRness follows from [28, Theorem 5.4 (A')], which guarantees that $X_1 > 0$ (see Appendix D). The last assertion of item 5 is a consequence of item 4. Consider finally item 6. One has $(B_3^\top X_3 B_3)^{-1}(K) = \{\lambda \in \mathbb{R}^m | B_3^\top X_3 B_3 \lambda = \alpha, \alpha_i \ge 0, 1 \le i \le m\}$, see e.g., [53, Example 15]. Hence $\mathrm{Ker}(B_3^\top X_3 B_3) \cap (B_3^\top X_3 B_3)^{-1}(K) = \mathrm{Ker}(B_3^\top X_3 B_3)$. $\square$

In item 2, the sufficient condition $\langle q_k, v \rangle > 0$ for all $v \in S$ is equivalent to $q_k \in S^* \setminus S^\perp$. One notices that under the conditions of item 4 in Proposition 5.1, the DI in (5.3) is equivalent to a DI with maximal monotone set-valued term:

$$ (5.8) \qquad \dot{\zeta}(t) - \mathcal{R}^{-1}\mathcal{A}\mathcal{R}^{-1}\zeta(t) + \mathcal{R}^{-1}G(t) \in -\mathcal{N}_{\tilde{\Phi}}(\zeta(t)) $$

with $\zeta = \mathcal{R}\zeta$. This fits within the class of DIs analysed for instance in [14]. Under basic assumption on the regularity of $G(\cdot)$, (5.8) has Lipschitz continuous solutions with uniqueness.

**Proposition 5.2.** Assume that $X_1 > 0$, and $G(\cdot)$ is bounded. The problem VI($\mathcal{P}$, $q_k$, $\Psi_{\Phi}$) always has at least one solution for bounded $\xi_k$.

**Proof.** Consider item 2 in Proposition 5.1, with $\xi^0 = \begin{pmatrix} x^0 \\ \lambda^0 \end{pmatrix}$, $v = \begin{pmatrix} v_x \\ v_\lambda \end{pmatrix}$. Then:

$$ (5.9) \qquad \begin{aligned} \langle q_k - \mathcal{P}^\top \xi^0, v \rangle &= \langle -\mathcal{P}_h \begin{pmatrix} x_{1,k} \\ \lambda_k \end{pmatrix} + hG_k, v \rangle - \langle B_3^\top X_3 B_3 \lambda^0, v_\lambda \rangle \\ & \quad - \langle X_1 x^0, v_x \rangle \\ &= \langle - \begin{pmatrix} X_1 + hX_1A_1 & hX_1B_1 \\ -hC_1 & B_3^\top X_3 B_3 - h\hat{D} \end{pmatrix} \begin{pmatrix} x_{1,k} \\ \lambda_k \end{pmatrix} + hG_k, v \rangle \\ & \quad - \langle B_3^\top X_3 B_3 \lambda^0, v_\lambda \rangle - \langle X_1 x^0, v_x \rangle. \end{aligned} $$

Since $\xi^0 \in \Phi$ it follows that $x^0 \in \mathbb{R}^{n_1}$ while $\lambda^0 \in \mathbb{R}_+^m$, and similarly $v_x \in \mathbb{R}^{n_1}$, $v_\lambda \in \mathbb{R}_+^m$. Since $v \ne 0$ let us take $x^0 = -\alpha X_1^{-1} \frac{v_x}{\|v_x\|} - \beta X_1^{-1} \frac{v_\lambda}{\|v_\lambda\|} \|v_\lambda\|$ for some $\alpha \ge \beta \ge 0$, and $\lambda^0 = 0$. It follows that:

$$ (5.10) \qquad \begin{aligned} \langle q_k - \mathcal{P}^\top \xi^0, v \rangle &= - \langle \begin{pmatrix} X_1 + hX_1A_1 & hX_1B_1 \\ -hC_1 & B_3^\top X_3 B_3 - h\hat{D} \end{pmatrix} \begin{pmatrix} x_{1,k} \\ \lambda_k \end{pmatrix}, v \rangle + \langle hG_k, v \rangle \\ & \ge - \| \begin{pmatrix} X_1 + hX_1A_1 & hX_1B_1 \\ -hC_1 & B_3^\top X_3 B_3 - h\hat{D} \end{pmatrix} \begin{pmatrix} x_{1,k} \\ \lambda_k \end{pmatrix} \| \|v\| - \|hG_k\| \|v\| \\ & \quad + \alpha \|v_x\| + \beta \|v_\lambda\| \\ & \ge - \| \begin{pmatrix} X_1 + hX_1A_1 & hX_1B_1 \\ -hC_1 & B_3^\top X_3 B_3 - h\hat{D} \end{pmatrix} \begin{pmatrix} x_{1,k} \\ \lambda_k \end{pmatrix} \| (\|v_x\| + \|v_\lambda\|) \\ & \quad - \|hG_k\| (\|v_x\| + \|v_\lambda\|) + \alpha \|v_x\| + \beta \|v_\lambda\|. \end{aligned} $$
Since $G(\cdot)$ is bounded, for large enough $\alpha$ and $\beta$ it follows that $\langle q_k - \mathcal{P}^\top \xi^0, v \rangle > 0$ and this holds for any $v \neq 0$ inside $S$. $\square$

## 5.2 BOUNDEDNESS AND CONVERGENCE ANALYSIS

In the next analysis we assume that $B_3^\top X_3 B_3$ is singular, which is the most general case (see the example in Section 6.2, see also Section 5.3 for the invertible case), and that $K = \mathbb{R}_+^m$ and that $X_1 > 0$. From the structure of $\mathcal{P}$ and of the set-valued right-hand side, it follows from (5.5) that $x_{1,k+1}$ is calculable from $x_{1,k}$ and $\lambda_k$:

$$ (5.11) \qquad \begin{aligned} X_1 x_{1,k+1} &= (X_1 + hX_1A_1)x_{1,k} + hX_1B_1\lambda_k - hG_k^x \\ \Leftrightarrow x_{1,k+1} &= (I_{n_1} + hA_1)x_{1,k} + hB_1\lambda_k - hX_1^{-1}G_k^x, \quad \text{for all } k \ge 0. \end{aligned} $$

where $G_k = \begin{pmatrix} G_k^x \\ G_k^\lambda \end{pmatrix}$.

Lemma 5.3. *From (5.11), for all $k \ge 0$ one has:*

$$ (5.12) \qquad x_{1,k+1} = (I_{n_1} + hA_1)^{k+1}x_{1,0} + h \sum_{i=0}^{k} (I_{n_1} + hA_1)^{k-i} B_1 \lambda_i - h \sum_{i=0}^{k} (I_{n_1} + hA_1)^{k-i} X_1^{-1} G_i^x. $$

*Proof.* The relation is true for $k=0$. Let it be true for $k \ge 1$. Then using (5.11) it follows that:

$$ (5.13) \qquad \begin{aligned} x_{1,k+2} &= (I_{n_1} + hA_1)x_{1,k+1} + hB_1\lambda_{k+1} - hX_1^{-1}G_{k+1}^x \\ &= (I_{n_1} + hA_1)\{(I_{n_1} + hA_1)^{k+1}x_{1,0} + h \sum_{i=0}^{k} (I_{n_1} + hA_1)^{k-i} B_1 \lambda_i \\ &\quad - h \sum_{i=0}^{k} (I_{n_1} + hA_1)^{k-i} X_1^{-1} G_i^x\} + hB_1\lambda_{k+1} - hX_1^{-1}G_{k+1}^x \\ &= (I_{n_1} + hA_1)^{k+2}x_{1,0} + h \sum_{i=0}^{k+1} (I_{n_1} + hA_1)^{k+1-i} B_1 \lambda_i \\ &\quad - h \sum_{i=0}^{k+1} (I_{n_1} + hA_1)^{k+1-i} X_1^{-1} G_i^x. \end{aligned} $$

By induction the relation in (5.12) is true for all $k \ge 0$. $\square$

Using (5.5) we can also express $\lambda_{k+1}$ as the solution of a generalised equation (GE):

$$ (5.14) \qquad \begin{aligned} & B_3^\top X_3 B_3 \lambda_{k+1} + (h\hat{D} - B_3^\top X_3 B_3)\lambda_k + hC_1 x_{1,k} + hG_k^\lambda \in -\mathcal{N}_{\mathbb{R}_+^m}(\lambda_{k+1}) \\ & \phantom{B_3^\top X_3 B_3 \lambda_{k+1} + (h\hat{D} - B_3^\top X_3 B_3)\lambda_k + hC_1 x_{1,k} + hG_k^\lambda \in} \underbrace{\triangleq q_k^\lambda(\lambda_k, x_{1,k}, G_k^\lambda)}_{\Downarrow} \\ & \lambda_{k+1} \in (B_3^\top X_3 B_3 + \mathcal{N}_{\mathbb{R}_+^m})^{-1} \left( (B_3^\top X_3 B_3 - h\hat{D})\lambda_k - hC_1 x_{1,k} - hG_k^\lambda \right) \end{aligned} $$

for all $k \ge 0$. The second formulation shows the iterative feature of this GE for $\lambda_k$, which makes it quite different from the LCCP (3.3). It is equivalent to the LCP($B_3^\top X_3 B_3$, $q_k^\lambda(\lambda_k, x_{1,k}, G_k^\lambda)$). Under the conditions stated in Proposition 5.1 and Corollary 4.1, the set of solutions to this LCP is non empty, and it contains a least-norm element $\lambda_{k+1}^{\min}$, $k \ge 0$, satisfying $\|\lambda_{k+1}^{\min}(\lambda_k, x_{1,k}, G_k^\lambda)\| \le \eta \|q_k^\lambda(\lambda_k, x_{1,k}, G_k^\lambda)\|$, $k \ge 0$, where $\eta > 0$ depends only on $B_3^\top X_3 B_3$. Therefore $\|\lambda_{k+1}^{\min}(\lambda_k, x_{1,k}, G_k^\lambda)\| \le \eta \|h\hat{D} - B_3^\top X_3 B_3\| \|\lambda_k^{\min}\| + h\eta \|C_1 x_{1,k}\| + h\eta \|G_k^\lambda\|$. Let us denote $h_3 \triangleq \|h\hat{D} - B_3^\top X_3 B_3\|$.
Lemma 5.4. One has for all $k \ge 0$:

$$ (5.15) \qquad \begin{aligned} \|\lambda_{k+2}^{\min}(\lambda_{k+1}, x_{1,k+1}, G_{k+1}^{\lambda})\| & \le (\eta h_3)^{k+1} \|\lambda_1^{\min}(\lambda_0, x_{1,0}, G_0^{\lambda})\| \\ & \quad + h\eta \sum_{i=0}^{k} (\eta h_3)^{k-i} (\|C_1 x_{1,i+1}\| + \|G_{i+1}^{\lambda}\|) \end{aligned} $$

*Proof.* The inequality is satisfied for $k = 0$. Assume that it is true for some $k \ge 1$. Then

$$ (5.16) \qquad \begin{aligned} \|\lambda_{k+3}^{\min}(\lambda_{k+2}, x_{1,k+2}, G_{k+2}^{\lambda})\| & \le \eta \|h\hat{D} - B_3^\top X_3 B_3\| \|\lambda_{k+2}^{\min}\| + h\eta \|C_1 x_{1,k+2}\| + h\eta \|G_{k+2}^{\lambda}\| \\ & \le \eta \|h\hat{D} - B_3^\top X_3 B_3\| \left\{ \eta^{k+1} \|h\hat{D} - B_3^\top X_3 B_3\|^{k+1} \|\lambda_1^{\min}(\lambda_0, x_{1,0}, G_0^{\lambda})\| \right. \\ & \qquad \left. + h\eta \sum_{i=0}^{k} \eta^{k-i} \|h\hat{D} - B_3^\top X_3 B_3\|^{k-i} (\|C_1 x_{1,i+1}\| + \|G_{i+1}^{\lambda}\|) \right\} \\ & \quad + h\eta \|C_1 x_{1,k+2}\| + h\eta \|G_{k+2}^{\lambda}\| \\ & = \eta^{k+2} \|h\hat{D} - B_3^\top X_3 B_3\|^{k+2} \|\lambda_1^{\min}(\lambda_0, x_{1,0}, G_0^{\lambda})\| \\ & \quad + h\eta \sum_{i=0}^{k} \eta^{k+1-i} \|h\hat{D} - B_3^\top X_3 B_3\|^{k+1-i} (\|C_1 x_{1,i+1}\| + \|G_{i+1}^{\lambda}\|) \\ & \quad + h\eta \|C_1 x_{1,k+2}\| + h\eta \|G_{k+2}^{\lambda}\| \\ & = \eta^{k+2} \|h\hat{D} - B_3^\top X_3 B_3\|^{k+2} \|\lambda_1^{\min}(\lambda_0, x_{1,0}, G_0^{\lambda})\| \\ & \quad + h\eta \sum_{i=0}^{k+1} \eta^{k+1-i} \|h\hat{D} - B_3^\top X_3 B_3\|^{k+1-i} (\|C_1 x_{1,i+1}\| + \|G_{i+1}^{\lambda}\|). \end{aligned} $$

The result follows by induction. $\square$

The next step is to use Lemmas 5.3 and 5.4 to get a new upper bound. Let us denote $h_1 \triangleq \|I_{n_1} + hA_1\|$. Then using (5.12) and (5.15) one has:

$$ (5.17) \qquad \left\{ \begin{aligned} \|x_{1,k+1}\| & \le h_1^{k+1} \|x_{1,0}\| + h \sum_{i=0}^{k} h_1^{k-i} \|B_1\| \|\lambda_{i+1}^{\min}\| + h \sum_{i=0}^{k} h_1^{k-i} \|X_1^{-1}\| \|G_i^x\| \\ \|\lambda_{k+2}^{\min}\| & \le \eta^{k+1} h_3^{k+1} \|\lambda_1^{\min}\| + h\eta \sum_{i=0}^{k} \eta^{k-i} h_3^{k-i} (\|C_1 x_{1,i+1}\| + \|G_{i+1}^{\lambda}\|) \end{aligned} \right. $$

Therefore we obtain:

$$ (5.18) \qquad \begin{aligned} \|x_{1,k+1}\| & \le h_1^{k+1} \|x_{1,0}\| + h \sum_{i=0}^{k} h_1^{k-i} \|X_1^{-1}\| \|G_i^x\| \\ & \quad + h \sum_{i=0}^{k} h_1^{k-i} \|B_1\| \eta^i h_3^i \|\lambda_1^{\min}\| + h h_1^k \|B_1\| \|\lambda_1^{\min}\| \\ & \quad + h \sum_{i=1}^{k} h_1^{k-i} \|B_1\| \eta h \sum_{j=0}^{i-1} \eta^{i-1-j} h_3^{i-1-j} \|C_1\| \|G_{j+1}^{\lambda}\| \\ & \quad + h \sum_{i=1}^{k} h_1^{k-i} \|B_1\| h \eta \sum_{j=0}^{i-1} \eta^{i-1-j} h_3^{i-1-j} \|C_1\| \|x_{1,j+1}\|, \end{aligned} $$
which we rewrite as:

$$ (5.19) \qquad \|x_{1,k+1}\| \le h_1^{k+1} \|x_{1,0}\| + h f_k + h^2 \eta \|C_1\| \|B_1\| \sum_{i=1}^{k} h_1^{k-i} \sum_{j=0}^{i-1} (\eta h_3)^{i-1-j} \|x_{1,j+1}\| $$

with

$$ (5.20) \qquad f_k = \sum_{i=0}^{k} h_1^{k-i} \|X_1^{-1}\| \|G_i^x\| + \|B_1\| \|\lambda_1^{\min}\| \sum_{i=0}^{k} h_1^{k-i} (\eta h_3)^i + h_1^k \|B_1\| \|\lambda_1^{\min}\| \\ + h\eta \|B_1\| \|C_1\| \sum_{i=1}^{k} h_1^{k-i} \sum_{j=0}^{i-1} (\eta h_3)^{i-1-j} \|G_{j+1}^{\lambda}\|. $$

The third term in the right-hand side of (5.19) involves terms from $\|x_{1,1}\|$ to $\|x_{1,k}\|$. The inequality (5.19) can be rewritten as:

$$ (5.21) \qquad \|x_{1,k+1}\| \le h_1^{k+1} \|x_{1,0}\| + h f_k + h^2 \sum_{i=1}^{k} \beta_i \|x_{1,i}\|, $$

for some positive sequence of real numbers $\{\beta_i\}$, each $\beta_i$ being calculated from the expression in (5.19) $\eta \|C_1\| \|B_1\| \sum_{i=1}^{k} h_1^{k-i} \sum_{j=0}^{i-1} (\eta h_3)^{i-1-j} \|x_{1,j+1}\|$, and thus being made of powers of $h_1$ and of $h_3\eta$. Calculations show that (see Appendix F)

$$ (5.22) \qquad \beta_i = \eta \|C_1\| \|B_1\| \sum_{j=i}^{k} h_1^{k-j} (\eta h_3)^{j-i}, \quad i \ge 1. $$

Let us now apply a version of the discrete-time Gronwall inequality, due to Bachpatte [13] [47, Lemma 100], see Appendix E. It is deduced that for all $k \ge 0$:

$$ (5.23) \qquad \|x_{1,k+1}\| \le h_1^{k+1} \|x_{1,0}\| + h f_k + h^2 \sum_{i=0}^{k} (h_1^{i+1} \beta_i \|x_{1,0}\| + h \beta_i f_i) \prod_{j=i+1}^{k+1} (1+h^2 \beta_j). $$

Let us state the next assumption (recall that $h_3 = \|h\hat{D} - B_3^\top X_3 B_3\|$ and $h_1 = \|I_{n_1} + hA_1\|$).

**Assumption 4.** One has a) $h_1 \le 1$ for $h > 0$ small enough, b) $\|G(t)\|$ is uniformly bounded by $G_{\max} < +\infty$ on $[0, T]$ for any $T > 0$, and $G(\cdot)$ is Lipschitz continuous with Lipschitz constant $k_G \ge 0$, c) $\eta h_3 \le 1$ for $h > 0$ small enough.

Inequality a) is similar to **Assumption 1**. As illustrated through examples in Section 6, it may hold in practice for reasonable values of the time step. See also the developments in Appendix A for conditions on its validity, and **Remark 5.7** for its relaxation. **Assumption 4 b)** is a standard property of exogenous time-dependent terms. The inequality c) is the central assumption of **Assumption 4**, which quantifies the couplings between the DAE part and the complementarity part of the LDCS. Some examples of circuits in Section 6 will prove that it can be satisfied in practice. Notice from (5.14) that if $B_3^\top X_3 B_3 = 0$ (Assumption 2), then the minimal-norm element of the set of solutions to the LCP is always $\lambda_k^{\min} = 0$, hence $\eta = 0$ in that case. Using (5.20) the next upper bound can be found:

$$ (5.24) \qquad f_k \le (k+1)(G_{\max}\|X_1^{-1}\| + \|B_1\| \|\lambda_1^{\min}\|) + \|B_1\| \|\lambda_1^{\min}\| + h\eta \|B_1\| \|C_1\| G_{\max} \sum_{i=1}^{k} \sum_{j=0}^{i-1} j, $$

Since $h = \frac{T}{N}$ and $k \in \{0, N-1\}$:

$$ (5.25) \qquad h f_k \le T \frac{(k+1)}{N} (G_{\max} \|X_1^{-1}\| + \|B_1\| \|\lambda_1^{\min}\|) + \frac{T}{N} \|B_1\| \|\lambda_1^{\min}\| + \frac{T^2}{n^2} \eta \|B_1\| \|C_1\| G_{\max} \sum_{i=1}^{k} \sum_{j=0}^{i-1} j. $$
Now we have that $\sum_{i=1}^{k} \sum_{j=0}^{i-1} j = (k-1)+2(k-1)+3(k-3)+...+(k-3)3+(k-2)2+(k-1)$, where the sum consists of $k-1$ terms. Hence $\sum_{i=1}^{k} \sum_{j=0}^{i-1} j = O(k^2)$, and $\frac{T^2}{N^2}\eta\|B_1\| \|C_1\| G_{\max} \sum_{i=1}^{k} \sum_{j=0}^{i-1} j = T^2\eta\|B_1\| \|C_1\| G_{\max} \frac{O(k^2)}{N^2}$. Since $k \le N$ we infer that this term is uniformly bounded in $k$, for each fixed bounded $T$. The same holds for $T \frac{(k+1)}{N} (G_{\max}\|X_1^{-1}\| + \|B_1\| \|\lambda_1^{\min}\|)$. Thus $hf_k \le \alpha_1 \max(T, T^2)$ for some $\alpha_1$ which depends on $\eta$, $\|B_1\|$, $\|C_1\|$, $G_{\max}$, $\|\lambda_1^{\min}\|$.

The next step consists of computing an upper bound from (5.23) as:

$$ (5.26) \qquad \|x_{1,k+1}\| \le \left( \|x_{1,0}\| + \alpha_1 \max(T, T^2) \right) \left( 1 + \frac{T^2}{N^2} \sum_{i=0}^{k} \beta_i \Pi_{j=i+1}^{k+1} \left( 1 + \frac{T^2}{N^2} \beta_j \right) \right). $$

In view of Assumption 4 and (5.22), one has $\beta_i \le \eta \|C_1\| \|B_1\| \sum_{j=i}^{k} 1 = \eta \|C_1\| \|B_1\|(k-i) \le \eta \|C_1\| \|B_1\|N, i \ge 1, k \le N$.

Consequently:

$$ (5.27) \qquad \begin{aligned} \|x_{1,k+1}\| &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \frac{T^2}{N} \sum_{i=0}^{k} \eta \|C_1\| \|B_1\| \Pi_{j=i+1}^{k+1} (1 + \frac{T^2}{N} \eta \|C_1\| \|B_1\|)) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} \sum_{i=0}^{k} \Pi_{j=i+1}^{k+1} (1 + \frac{T^2}{N} \eta \|C_1\| \|B_1\|)) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} \sum_{i=0}^{k} (1 + \frac{T^2}{N} \eta \|C_1\| \|B_1\|)^{k-i}) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} \sum_{i=0}^{k} (1 + \frac{\alpha_2}{N})^{k-i}) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} \sum_{i=0}^{N} (1 + \frac{\alpha_2}{N})^N) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} (N+1) (1 + \frac{\alpha_2}{N})^N) \\ &\le (\|x_{1,0}\| + \alpha_1 \max(T, T^2))(1 + T^2 \eta \|C_1\| \|B_1\| \frac{N+1}{N} (1 + \frac{\alpha_2}{N})^N), \end{aligned} $$

where $\alpha_2 = T^2\eta\|C_1\| \|B_1\|$. One has $(1 + \frac{\alpha_2}{N})^N \to e^{\alpha_2}$ as $N \to +\infty$ [73, p.196]. Therefore the sequence $\{x_{1,k}\}_k$ is uniformly bounded over $k \in \{0, N-1\}$. It directly follows from Lemma 5.4 that the sequence $\{\lambda_k^{\min}\}_k$ is uniformly bounded as well. It follows from (4.1) that the sequence $\{\frac{x_{1,k+1}-x_{1,k}}{h}\}_k$ is uniformly bounded also. Therefore we are ready to prove the following result. The piece-wise linear approximations are defined in (4.7), and the step approximations $x_1^{h*}(t) = x_{1,k+1}$ if $t \in [t_k, t_{k+1})$, $\lambda_{h*}^{\min}(t) = \lambda_k^{\min}$ if $t \in [t_k, t_{k+1})$.

**Proposition 5.5.** Let Assumption 4 hold true, as well as the conditions for the well-posedness of the OSNSP. Then the sequence $\{x_1^h(\cdot)\}_{h \ge 0}$ converges uniformly towards a continuous limit $x_1(\cdot)$. Moreover $x_1^{h*}(\cdot)$ converges strongly in $L_2([0, T]; \mathbb{R}^{n_1})$ towards $x_1(\cdot)$, and $\dot{x}_1^h(\cdot)$ converges weakly in $L_2([0, T]; \mathbb{R}^{n_1})$ towards $\dot{x}_1(\cdot)$. Therefore $\{(B_1\lambda_{h*}^{\min})(\cdot)\}_{h \ge 0}$ converges weakly in $L_2([0, T]; \mathbb{R}^{n_1})$ towards a limit $\lambda_1(\cdot)$, and $\dot{x}_1(t) = A_1x_1(t) + G^x(t) + \lambda_1(t)$ almost everywhere on $[0, T]$.

**Proof.** Both sequences $\{x_1^h(\cdot)\}$ and $\{\dot{x}_1^h(\cdot)\}$ are uniformly bounded on $[0, T]$, the sequence $\{x_1^h(\cdot)\}$ of $C^0([0, T]; \mathbb{R}^{n_1})$ functions is thus equicontinuous. From the Ascoli-Arzela Theorem, $\{x_1^h(\cdot)\}$ stays in a compact subset of $C^0([0, T]; \mathbb{R}^{n_1})$ and indeed it converges uniformly towards a continuous limit $x_1(\cdot)$. Let us now examine the sequences of piecewise-linear functions $\{G^h(\cdot)\}_{h \ge 0}$ and $\{\dot{G}^h(\cdot)\}_{h \ge 0}$,
designed similarly as (4.7). From Assumption 4 b) one infers that $G^h(t)$ is uniformly bounded on $[0, T]$, $\frac{\|G_{k+1}-G_k\|}{h} \le k_E$ implies that $\dot{G}^h(t)$ is also uniformly bounded in $h$ on $[0, T]$, thus $\{G^h(\cdot)\}_{h \ge 0}$ is a sequence of Lipschitz continuous, uniformly bounded, equicontinuous functions. From the Ascoli-Arzela Theorem, $\{G^h(\cdot)\}$ stays in a compact subset of $C^0([0, T]; \mathbb{R}^{n_1})$ and it converges uniformly towards the continuous limit $G(\cdot)$ as $h \to 0$. The remaining part of the proof is similar to the proof of Corollary 4.6 (adding $G^{x,h}(\cdot)$ with the required properties) and it is omitted. □

**Remark 5.6.** Assumption 4 c) is clearly less restrictive than Assumption 2. It still imposes some constraint on the matrix $C_3B_3$. Examples are treated next which show that it is satisfied in some practical cases.

Let us now prove the uniform boundedness of the sequences $\{\lambda_k^{\min}\}_k$ and $\{x_{2,k}\}_k$. To this aim we shall use first the second line in (5.17) and (5.27). Let us denote the right-hand side of (5.27) as $x_{1,\max}(T, \eta, x_{1,0})$. Then from (5.17) and Assumption 4 we obtain:

$$ (5.28) \qquad \begin{aligned} \|\lambda_{k+1}^{\min}\| &\le \|\lambda_0^{\min}\| + h\eta \sum_{i=0}^{k} (\|C_1\| x_{1,\max}(T, \eta, x_{1,0}) + \|G_{i+1}^{\lambda}\|) \\ &\le \|\lambda_0^{\min}\| + h\eta \|C_1\| x_{1,\max}(T, \eta, x_{1,0}) (k+1) + h\eta \sum_{i=0}^{k} \|G_{i+1}^{\lambda}\| \\ &\le \|\lambda_0^{\min}\| + \frac{T(k+1)}{N}\eta \|C_1\| x_{1,\max}(T, \eta, x_{1,0}) + \eta \frac{T(k+1)}{N} G_{\max} \\ &\le \|\lambda_0^{\min}\| + \frac{T(N+1)}{N}\eta \|C_1\| x_{1,\max}(T, \eta, x_{1,0}) + \eta \frac{T(N+1)}{N} G_{\max}. \end{aligned} $$

Therefore $\|\lambda_{k+1}^{\min}\|$ is uniformly bounded for any bounded $T$ and initial data. The boundedness of $x_{2,k}$ follows from the equality constraint in (4.1). A result similar to Corollary 4.6 holds, but this time verifying the complementarity conditions is less easy due to the presence of $z$ in the complementarity variable $w$.

**Remark 5.7 (relaxation of assumptions).** Assumption 4 a) is convenient to simplify the upper bounds calculations, however it is somewhat restrictive. One path to relax it, is to study the terms $h_1^k$ another way, since for any $a_1 > 0$ one has $(1 + \frac{T}{N}a_1)^N \to e^{a_1T}$ as $N \to +\infty$, while this sequence is increasing. Consider first (5.20). Taking into account Assumption 4 b) and c), one sees that the first two terms in $hf_k$ are upper bounded by constants times

$$ (5.29) \qquad h \sum_{i=0}^{k} h_1^{k-i} = h \sum_{i=0}^{k} h_1^i \le h k (1 + h a_1)^k, $$

where $h_1 \le \|I_{n_1}\| + h\|A_1\| = 1 + h a_1$ for some $a_1 \ge 0$. We have $hk(1 + h a_1)^k \le \frac{T}{N}k(1 + \frac{Ta_1}{N})^N \to Te^{a_1T}$ as $N \to +\infty$. The third term of $hf_k$ is upper bounded by a constant times $hh_1^k \le \frac{T}{N}(1 + h a_1)^N \to 0$ as $N \to +\infty$. The last term of $hf_k$ is upper bounded by a constant multiplied by

$$ (5.30) \qquad h^2 \sum_{i=1}^{k} h_1^{k-i} k \le \frac{T^2}{N^2} k^2 (1 + h a_1)^k \le \frac{T^2}{N^2} N^2 (1 + h a_1)^N \to T^2 e^{Ta_1} $$

as $N \to +\infty$. These first calculations allow us to upper bound the term $hf_k$ in (5.23). The first term $h_1^{k+1}\|x_{1,0}\| = (1 + h a_1)^{k+1}\|x_{1,0}\| \le (1 + h a_1)^N \|x_{1,0}\| \to \|x_{1,0}\|e^{a_1T}$ as $N \to +\infty$.

The third term in (5.23) is more complex: $h^2 \sum_{i=0}^{k} (h_1^{i+1} \beta_i \|x_{1,0}\| + h \beta_i f_i) \prod_{j=i+1}^{k+1} (1 + h^2 \beta_j)$. Since from Assumption 4 c) $\eta h_3 \le 1$, it follows that

$$ (5.31) \qquad \beta_i \le \eta \|C_1\| \|B_1\| \sum_{j=1}^{k} h_1^{k-j} \le \eta \|C_1\| \|B_1\| \sum_{j=1}^{k} (1 + h a_1)^{k-j} $$
for some $a_1 > 0$. Thus

$$ (5.32) \quad h^2 \beta_j \le \eta \|C_1\| \|B_1\| h^2 k (1 + ha_1)^{k-1} \le \eta \|C_1\| \|B_1\| \frac{T^2}{N^2} N (1 + \frac{T}{N} a_1)^N = \eta \|C_1\| \|B_1\| \frac{T^2}{N} (1 + \frac{T}{N} a_1)^N. $$

Hence

$$ (5.33) \quad \Pi_{j=i+1}^{k+1} (1 + h^2 \beta_j) \le (1 + \eta \|C_1\| \|B_1\| \frac{T^2}{N} (1 + \frac{T}{N} a_1)^N)^N \to e^{\eta \|C_1\| \|B_1\| T^2 e^{a_1 T}}. $$

We have

$$ (5.34) \quad h^2 \sum_{i=0}^{k} \beta_i h f_i = \sum_{i=0}^{k} (h^2 \beta_i) (h f_i) \le N \eta \|C_1\| \|B_1\| \frac{T^2}{N} (\alpha_1 T + \alpha_2 \frac{T}{N} + \alpha_3 T^2) (1 + \frac{T}{N} a_1)^{2N} \\ \to \eta \|C_1\| \|B_1\| T^2 (\alpha_1 T + \alpha_3 T^2) e^{2a_1 T}, $$

hence is bounded for any bounded $T$. We have also

$$ (5.35) \quad h^2 \sum_{i=0}^{k} h_1^{i+1} \beta_i \|x_{1,0}\| = \|x_{1,0}\| \sum_{i=0}^{k} h_1^{i+1} h^2 \beta_i \le \eta \|C_1\| \|B_1\| \|x_{1,0}\| \sum_{i=0}^{k} (1 + ha_1)^k h^2 k (1 + ha_1)^{k-1} \\ \le \eta \|C_1\| \|B_1\| \|x_{1,0}\| \sum_{i=0}^{N} (1 + ha_1)^{2N} \frac{T^2}{N^2} N \\ \le \eta \|C_1\| \|B_1\| \|x_{1,0}\| T^2 e^{2a_1 T}. $$

Consequently from (5.23) it is inferred that if Assumption 4 b) and c) holds, then $\|x_{1,k+1}\|$ is bounded by a term that depends only on $T$, initial conditions and constants.

Let us examine now Assumption 4 c) and the terms which which involve $h_3^k$. We have $h_3 = \|h\hat{D} - B_3^\top X_3 B_3\|$. Let us assume that $\hat{D}$ is invertible. Then (5.36)

$$ \begin{aligned} \|h\hat{D} - B_3^\top X_3 B_3\| = \|\hat{D}(h - \hat{D}^{-1} B_3^\top X_3 B_3)\| &\le \|\hat{D}\| \|(h I_m - \hat{D}^{-1} B_3^\top X_3 B_3)\| \le \|\hat{D}\| (h + \|\hat{D}^{-1} B_3^\top X_3 B_3\|) \\ &\le \|\hat{D}\| (h + \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|) \\ &= \|\hat{D}\| \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\| \left(1 + \frac{1}{\|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|} h\right). \end{aligned} $$

Therefore it follows that $h_3^k \le (\|\hat{D}\| \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|)^k \left(1 + \frac{T}{\|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|} \frac{1}{N}\right)^k$. Consequently $(\eta h_3)^k \le (\eta \|\hat{D}\| \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|)^k \left(1 + \frac{T}{\|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|} \frac{1}{N}\right)^k$ which converges as $N \to +\infty$ to $(\eta \|\hat{D}\| \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|)^N e^{\frac{T}{\|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\|}}$. A necessary condition for the limit to be bounded is $\eta \|\hat{D}\| \|\hat{D}^{-1}\| \|B_3^\top X_3 B_3\| \le 1$: it is unclear whether or not this inequality is less conservative that the one in Assumption 4 c).

We conclude that while the relaxation of Assumption 4 a) is reasonable, the relaxation of Assumption 4 c) is not obvious.

## 5.3 CASE $X_1>0$ AND $B_3^\top X_3 B_3>0$

This is the case studied in item 4 of Proposition 5.1 and may be verified in practice, see the example in Section 6.1. Then the DI in (5.3) is equivalently rewritten as in (5.8). Consequently one can apply the implicit Euler scheme and all the existence, convergence and order results presented in [14] (see that
article for details on the discretization method and theoretical results), to (5.8). It is inferred that in this case (5.8) possesses unique solutions in $C^0([0, T]; \mathbb{R}^{n_1+m})$ with essentially bounded derivatives [14, Proposition 2.6]. It is deduced that $x_2(\cdot)$ is also unique in $C^0([0, T]; \mathbb{R}^{n_1+m})$ with essentially bounded derivative, while the algebraic variable $z(\cdot)$ (which depends on $\dot{\lambda}(\cdot)$, see (5.1)) is essentially bounded (hence a priori it may be discontinuous).

# 6 EXAMPLES AND NUMERICAL SIMULATIONS

Two examples of circuits with ideal diodes, and one academic example, are developed in this section. Let us comment that the analysis of circuits applies to hydraulics, due to the strong analogy between hydraulic and electrical circuits [68, 56] [52, Chapter 1], where check valves play the role of diodes.

## 6.1 FIRST EXAMPLE

<figure><img src="image_3.png" /><figcaption>Figure 1: An RLCD passive circuit.</figcaption></figure>

Let us consider the circuit in Figure 1. Let us assume first that $R_1 = R_3 = 0 \, \Omega$. Let $\xi_1 = i_2$, $\xi_2 = \int_0^t i_3(t)dt$, $\xi_3 = i_1$. The dynamics is given by:

$$ (6.1) \qquad \left\{ \begin{array}{l} \dot{\xi}_1(t) = -\frac{R_2}{L}\xi_1(t) + \frac{1}{L}u_D + \frac{1}{L}u_1(t) \\ \dot{\xi}_2(t) = \xi_1(t) - \xi_3(t) \\ 0 = \frac{1}{C}\xi_2(t) + u_D(t) + u_2(t) \\ 0 \le u_D \perp \xi_3(t) \ge 0. \end{array} \right. $$

It is noteworthy that this circuit dynamics agrees with [2, Theorem 3.6] [48] about the index (indeed it cannot be of index 1 due to the loop with a capacitive element). In our notations we have $x = (\xi_1, \xi_2)^T$, $z = \xi_3$, $C = (0 \ 0 \ 1)$, $B = (\frac{1}{L} \ 0 \ 1)^T$. The transfer function $H(s) = \frac{z(s)}{u_D(s)} = \frac{Cs^2 + \frac{CR}{L}s + \frac{1}{L}}{s + \frac{R}{L}} = Cs + \frac{1}{Ls + R}$. It has a relative degree $r = -1$. It can be checked that $H(s)$ is positive real (PR) as long as $R_2 > 0$ [25, Theorem 2.45], but it is not strictly positive real (SPR) (the first condition in item 3.b in [25, Theorem 2.61] does not hold, despite of the fact that the second condition in the same item holds), and hence it is not strongly SPR. Let us verify the minimality of (6.1):

$$ 1. \ rank((A B) - sP B) = rank \begin{pmatrix} \frac{-R}{L} - s & 0 & 0 & \frac{1}{L} \\ 1 & -s & -1 & 0 \\ 0 & \frac{1}{C} & 0 & 1 \end{pmatrix} = 3 \ for \ all \ s \in \mathbb{C}. $$

$$ 2. \ rank(P B) = rank \begin{pmatrix} 1 & 0 & 0 & \frac{1}{L} \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} = 3. $$
$$3. \operatorname{rank}(A^\top - sP^\top C^\top) = \operatorname{rank}\begin{pmatrix} \frac{-R}{L} - s & 1 & 0 & 0 \\ 0 & -s & \frac{1}{C} & 0 \\ 0 & -1 & 0 & 1 \end{pmatrix} = 3 \text{ for all } s \in \mathbb{C}.$$

$$4. \operatorname{rank}(P^\top C^\top) = \operatorname{rank}\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} = 3.$$

$$5. \text{ (A B) Ker}(P) = \begin{pmatrix} \frac{-R}{L} & 0 & 0 \\ 1 & 0 & -1 \\ 0 & \frac{1}{C} & 0 \end{pmatrix} \text{ Ker}(P) = \{\xi \in \mathbb{R}^3 \mid \xi_1 = \xi_3 = 0, \xi_2 \in \mathbb{R}\} \subset \text{Im}(P) = \{\xi \in \mathbb{R}^3 \mid \xi_3 = 0, \xi_2 \in \mathbb{R}, \xi_1 \in \mathbb{R}\},$$

Therefore we infer that (6.1) admits a special Weiertrass form as (2.2). Let us choose $L = 4.7 \ 10^{-4}$

$$H, R = 1\Omega, C = 10^{-1} \text{ F. Let us define } S = \begin{pmatrix} -1.414 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 10^{-4} \end{pmatrix}, T = \begin{pmatrix} -0.707 & 0 & 0 \\ 0 & 0 & 1 \\ -0.707 & -1 & 0 \end{pmatrix}. \text{ Then } SAT = \begin{pmatrix} -2.12 & 10^3 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}, SB = \begin{pmatrix} -3 & 10^{-3} \\ 0 \\ 10^{-4} \end{pmatrix}, CT = (-0.707 & -1 & 0), S \begin{pmatrix} \frac{1}{L} & 0 \\ 0 & 0 \\ 0 & 1 \end{pmatrix} = \begin{pmatrix} 0.301 & 10^4 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 10^{-4} \end{pmatrix}. \text{ Let us make}$$

the state variable change $\zeta = T^{-1}\xi$, then we obtain the transformed dynamics:

$$\text{(6.2)} \qquad \left\{ \begin{aligned} \dot{\zeta}_1(t) &= -2.12 \ 10^3 \zeta_1(t) - 3 \ 10^{-3} \lambda(t) + 0.301 \ 10^4 u_1(t) \\ \dot{\zeta}_3(t) &= \zeta_2(t) \\ 0 &= \zeta_3(t) + 10^{-4} \lambda(t) + 10^{-4} u_2(t) \\ 0 \le \lambda(t) \perp w(t) &= -0.707 \zeta_1(t) - \zeta_2(t) \ge 0. \end{aligned} \right.$$

Let us now define $x_1 = \zeta_1, x_2 = \zeta_3, x_3 = \zeta_2$, then we obtain

$$\text{(6.3)} \qquad \left\{ \begin{aligned} \dot{x}_1(t) &= -2.12 \ 10^3 x_1(t) - 3 \ 10^{-3} \lambda(t) + 0.301 \ 10^4 u_1(t) \\ \dot{x}_2(t) &= x_3(t) \\ 0 &= x_2(t) + 10^{-4} \lambda(t) + 10^{-4} u_2(t) \\ 0 \le \lambda(t) \perp w(t) &= -0.707 x_1(t) - x_3(t) \ge 0, \end{aligned} \right.$$

which is the canonical form in (2.2) with $z = x_3, B_1 = -3.10^{-3}, B_2 = 0, B_3 = 10^{-4}, C_1 = -0.707, C_2 = 0, C_3 = -1, D = 0$. We have $X_3 = 10^4, A_1 = -2.12 \ 10^3, X_1 = \frac{C_1}{B_1} > 0, \hat{D} = D - C_3B_2 - C_2B_3 = 0, \bar{D}_h = h^2C_1B_1 + B_3^\top X_3B_3 = 2.121 \ 10^{-3}h^2 + 10^{-4}, h_3 = B_3^\top X_3B_3 = 10^{-4}, h_1 = |1 + hA_1| = |1 - 2.12 \ 10^3 h|$ which verifies Assumption 4 for $h > 0$ small enough. Notice that the set $S$ in Propositions 5.1 and 3.2 verifies $S = \{0\}$, hence the OSNSP is well-posed with unique solution. Provided that Assumption 4 is satisfied, the material in Section 5 applies. In a similar way to (5.1), the system (6.3) is transformed into:

$$\text{(6.4)} \qquad \left\{ \begin{aligned} \dot{x}_1(t) &= -2.12 \ 10^3 x_1(t) - 3 \ 10^{-3} \lambda(t) + 0.301 \ 10^4 u_1(t) \\ 10^{-4} \dot{\lambda}(t) &= 0.707 x_1(t) - 10^{-4} \dot{u}_2(t) + w(t) \\ x_2(t) &= -10^{-4} \lambda(t) - 10^{-4} u_2(t) \\ z(t) &= -10^{-4} \dot{\lambda}(t) - 10^{-4} \dot{u}_2(t). \end{aligned} \right.$$

From (5.14), the following LCP is obtained:

$$\text{(6.5)} \qquad \left\{ \begin{array}{l} 0 \le 10^{-4} \lambda_{k+1} + q_k \perp \lambda_{k+1} \ge 0 \\ q_k = -10^{-4} \lambda_k - 0.707 h x_{1,k} - 10^{-4} h u_{2,k}, \end{array} \right.$$
and $x_{1,k+1}$, $x_{2,k+1}$, $z_{k+1}$ as in (H.1). Let us check Assumptions 1 and 4 c). Since $\hat{D} = 0$ it is easy to show that $\eta = 10^4$ while $h_3 = 10^{-4}$, using the LCP in (5.14) that can be solved explicitly. Thus $\eta h_3 \le 1$ for all $h \ge 0$. However, $h_1 = \|1 - 2.12 10^3 h\| \le 1$ if and only if $h < 4.7 10^{-4}$. It is also noteworthy that this system verifies the conditions of Section 5.3, and item 4 of Proposition 5.1, since $X_1 > 0$ and $B_3^T X_3 B_3 = 10^{-4} > 0$. Hence the state space form (5.8), that is equivalent to (5.3) since $\mathcal{P} > 0$, can be calculated. Also the generalized equation (5.14) always has a unique solution $\lambda_{k+1}$. As we shall see through numerical simulations, the variable $z(\cdot) = x_3(\cdot)$ may be discontinuous. Actually this is in agreement with the theoretical conclusions from Section 5.3, which do not state that this variable is continuous. It also shows that a complete convergence analysis, including $z$, should rely on other tools than those used in this paper.

In Figure 2, the results of the numerical simulation of (6.3) with an implicit Euler scheme are depicted<fnref n="2" />. Numerical simulations are performed with the INRIA software package SICONOS [1]. At each time step, the Lemke solver [65, 64] is applied to solve the LCP (3.1) resulting from the discretization scheme. The initial conditions are set to $x(0) = (0.01, -10^{-4}, -7.07 \cdot 10^{-3})^T$, and the time step is $h = 5 \cdot 10^{-6}$s. The input functions $u_1(t)$ and $u_2(t)$ are taken as sinusoidal functions: $u_1(t) = \sin(t)$ and $u_2(t) = \cos(10t)$. Let us observe that the variable $z(t) = x_3(t)$ appears to be discontinuous at $t \approx 5.9$s. The multiplier $\lambda(t)$ remains continuous, see Figure 3. Finally, an order of convergence equal to one is deduced from the data in Figure 4. The error formula used is:

$$ (6.6) \qquad \text{error} = \|\mathbf{x}(t) - \mathbf{x}^h(t)\|_{\infty}. $$

Remark 6.1 (modelling). One can obtain a relative degree one LCS by eliminating the variable $u_D(t)$ in (6.1). However, for the sake of generality, and since in practice equality constraints are not always eliminated, we consider here a poorly formulated version of the circuit's dynamics, where the variables are badly matched to the constraints. Such bad matching could be avoided using structural analysis.

<figure><img src="image_4.png" /><figcaption>Figure 2: Implicit Euler numerical simulation of (6.3) with $h = 5 \cdot 10^{-6}$s. The numerical simulations of the variables $x_1(t)$, $x_2(t)$ and $x_3(t) = z(t)$ are depicted.</figcaption></figure>

<footnote n="2">²Source code of the examples can be found at https://github.com/siconos/siconos-tutorials/tree/master/.sandbox/mlcspaper</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Implicit Euler numerical simulation of (6.3) with $h = 5 \cdot 10^{-6}$s. The numerical simulations of the variables $\lambda(t)$ and $w(t)$ are depicted.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Numerical convergence of the implicit Euler scheme applied to (6.3). The reference curve is taken for $h = 5 \cdot 10^{-6}$s. The error formula used is given in (6.6).</figcaption></figure>
## 6.2 SECOND EXAMPLE

<figure><img src="image_7.png" /><figcaption>Figure 5: An RLCD passive circuit.</figcaption></figure>

This circuit example depicted in Figure 5 is inspired from [32]. The so-called modified nodal analysis is used to derive the dynamics [2, chapter 3]. It is assumed that the voltage source delivers a constant voltage $u(t) = u_c$, while the other voltage source $u_1(t)$ is time-varying Lipschitz continuous and bounded. Notice that the current/voltage characteristic of this device can be written as $i_5(t) \in -\partial\psi_{\{u_c\}}(u(t))$, or equivalently as $\mathbb{R} \ni i_5(t) \perp u(t) - u_c \in \{0\}$. The dynamics is given by:

$$ (6.7) \qquad \left\{ \begin{array}{l} C\dot{e}(t) = -R^{-1}e(t) + i_4(t) - i_1(t) - i_5(t) \\ L\frac{di_1}{dt}(t) = e(t) + u_1(t) \\ 0 = e(t) + u(t) \\ 0 \le u_D(t) \perp i_4(t) \ge 0 \end{array} \right. \quad \mathbb{R} \ni i_5(t) \perp u(t) - u_c \in \{0\}. $$

The state vector is $\xi = (\xi_1, \xi_2, \xi_3)^T = (e, i_1, i_5)^T$, and $u_D(t) = e(t)$. The complementarity conditions can be written as $K = \mathbb{R}^+ \times \mathbb{R} \ni (u_D(t), i_5(t)) \perp (i_4(t), u(t) - u_c) \in \mathbb{R}^+ \times \{0\} = K^*$. Setting $\lambda_1 = i_4$ and

$$ \lambda_2 = u(t) - u_c, \text{ we obtain } P = \begin{pmatrix} C & 0 & 0 \\ 0 & L & 0 \\ 0 & 0 & 0 \end{pmatrix}, A = \begin{pmatrix} -R^{-1} & -1 & -1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix}, B = \begin{pmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{pmatrix}, C = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}, D = 0, $$

and the LDCS form in (2.1) is equal to:

$$ (6.8) \qquad \left\{ \begin{array}{l} \dot{\xi}_1(t) = -\frac{1}{RC}\xi_1(t) - \frac{1}{C}\xi_2(t) - \frac{1}{C}\xi_3(t) + \frac{1}{C}\lambda_1(t) \\ \dot{\xi}_2(t) = \frac{1}{L}\xi_1(t) + \frac{u_1(t)}{L} \\ 0 = \xi_1(t) + \lambda_2(t) + u_c \\ 0 \le \xi_1(t) \perp \lambda_1(t) \ge 0 \\ \mathbb{R} \ni \xi_3(t) \perp \lambda_2(t) \in \{0\}. \end{array} \right. $$

$$ \text{Let } S = \begin{pmatrix} 0 & \frac{1}{L} & -\frac{1}{L} \\ -1 & 0 & -\frac{1}{R} \\ 0 & 0 & -C \end{pmatrix}, T = \begin{pmatrix} 0 & 0 & -\frac{1}{C} \\ 1 & 0 & 0 \\ -1 & 1 & 0 \end{pmatrix}, \text{ then } \begin{array}{l} \text{SAT} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \text{ SPT} = \begin{pmatrix} \frac{1}{L} & 0 & 0 \\ 0 & 0 & \frac{1}{C} \\ 0 & 0 & 0 \end{pmatrix}, \text{ CT} = \begin{pmatrix} 0 & 0 & -\frac{1}{C} \\ -1 & 1 & 0 \end{pmatrix}, \end{array} $$
$$SB = \begin{pmatrix} 0 & -\frac{1}{L} \\ 1 & -\frac{1}{R} \\ 0 & -C \end{pmatrix}. \text{ Using the state variable change } \zeta = T^{-1}\xi \text{ with } T^{-1} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 1 & 1 \\ -C & 0 & 0 \end{pmatrix}, \text{ one obtains:}$$

$$ (6.9) \qquad \left\{ \begin{aligned} \dot{\zeta}_1(t) &= -(\lambda_2(t) + u_c(t)) + \frac{u_1(t)}{L} \\ \frac{1}{C}\dot{\zeta}_3(t) &= \zeta_2(t) + \lambda_1(t) - \frac{1}{R}(\lambda_2(t) + u_c(t)) \\ 0 &= -C(\lambda_2(t) + u_c(t)) + \zeta_3(t) \\ K^* &\ni \lambda(t) \perp w(t) = \begin{pmatrix} 0 \\ -1 \end{pmatrix} \zeta_1(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} \zeta_2(t) + \begin{pmatrix} -\frac{1}{C} \\ 0 \end{pmatrix} \zeta_3(t) \in K. \end{aligned} \right. $$

Let $x_1 = \zeta_1, x_2 = \frac{1}{C}\zeta_3, z = \zeta_2, B_1 = (0 \ -1), B_2 = (1 \ -\frac{1}{R}), B_3 = (0 \ -1), C_1 = \begin{pmatrix} 0 \\ -1 \end{pmatrix}, C_3 = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, C_2 = \begin{pmatrix} -1 \\ 0 \end{pmatrix}$, $D=0$. This dynamics is rewritten in the canonical form as:

$$ (6.10) \qquad \left\{ \begin{aligned} \dot{x}_1(t) &= -\lambda_2(t) - u_c(t) + \frac{u_1(t)}{L} \\ \dot{x}_2(t) &= z(t) + \lambda_1(t) - \frac{1}{R}(\lambda_2(t) + u_c(t)) \\ 0 &= -\lambda_2(t) - u_c(t) + x_2(t) \\ K^* &\ni \lambda(t) \perp w(t) = \begin{pmatrix} 0 \\ -1 \end{pmatrix} x_1(t) + \begin{pmatrix} -1 \\ 0 \end{pmatrix} x_2(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} z(t) \in K. \end{aligned} \right. $$

In the simulations depicted in Figures 6a and 6b, it has been chosen $u_c(t) = \sin(10t) - 1 \le 0$, which complies with the system's intrinsic constraints. Then (when $u_1(t) = 0$) one recovers the standard form in (2.2), with a different class of sets K, however. One has $A_1 = 0, B_3^T X_3 = -C_3$ with $X_3 = 1$, and $X_1 B_1 = C_1^T$ (this follows necessarily from (D.2) and $Q_1 \le 0$ and using [25, Lemma A.69]), with $X_1 = 1$ (hence the quadruple $(A_1, B_1, C_1, \hat{D})$ is passive but it is not strictly passive and its associated transfer function is not strongly SPR). Also $\hat{D} = -C_3 B_2 - C_2 B_3 = \begin{pmatrix} 0 & 1 \\ -1 & \frac{1}{R} \end{pmatrix} \ge 0$. Also $B_3^T X_3 B_3 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} \ge 0$.

It is computed that $\bar{D}_h = \begin{pmatrix} 0 & h \\ -h & \frac{h}{R} + h^2 + 1 \end{pmatrix}$. The various sets for the set S in Proposition 3.2 are calculated as: $K^* = \mathbb{R}^+ \times \{0\}, \text{Ker}(\bar{D}_h + \bar{D}_h^T) = \begin{pmatrix} 0 & 0 \\ 0 & 2(\frac{h}{R} + h^2 + 1) \end{pmatrix} = \mathbb{R} \times \{0\}, \{w | \bar{D}_hw \in K = \mathbb{R}^+ \times \mathbb{R}\} = \mathbb{R} \times \mathbb{R}^+$. So $S = \mathbb{R}^+ \times \{0\} \ne 0$ and item 2 of Proposition 3.2 has to be applied. It is calculated that $\langle h\beta_k - \bar{D}_h\xi_0, v \rangle = \langle \begin{pmatrix} -h\xi_0^1 \\ -hx_{1,k} - x_{2,k} + h\xi_0^1 - (\frac{h}{2} + h^2 + 1)\xi_0^2 \end{pmatrix}, v \rangle$. Also $v \in K^*$ implies that $v_1 \ge 0$ and $v_2 = 0$, so it follows that $\langle h\beta_k - \bar{D}_h\xi_0, v \rangle + \psi_{K^*}(v) > 0$ is guaranteed, if $v_1 > 0$, by taking $\xi_0^1 < 0$. It is inferred that item 2 in Proposition 3.2 is verified and the OSNSP has at least one solution.

Any two solutions $\lambda_{k+1}^1$ and $\lambda_{k+1}^2$ satisfy $\lambda_{2,k+1}^1 = \lambda_{2,k+1}^2 = \lambda_{k+1}^2$ for all $h \ge 0$, in view of the form of $\bar{D}_h$. Then both $x_{1,k+1}$ and $x_{2,k+1}$ are unique due to the form of $B_1$ and $B_3$. For the convergence analysis, let us note that this system does not satisfy Assumption 2, therefore let us use the material in Section 5. Using the Frobenius norm one finds that $h_3 = \sqrt{2h^2 + (\frac{h}{R} - 1)^2}$, using the Hölder induced norm $\| \cdot \|_{2,2} = \sigma_{\max}(\cdot)$ one finds a smaller $h_3 = \sqrt{h^2 + (\frac{h}{R} - 1)^2}$, and $0 \le h_3 \le 1$ for all $h \in [0, \frac{2R}{1+R^2}]$. The material in Section 5 is based on the LCP in (5.14), for which the existence of a minimal norm solution satisfying the inequalities stated below (5.14) is guaranteed. However in this example $K = \mathbb{R}_+ \times \mathbb{R} \ne \mathbb{R}_+^2$. The following problem
has to be studied instead (it is noteworthy that the rest of the analysis remains unchanged as long as the minimal norm solution can be characterised):

$$ (6.11) \qquad K \ni B_3^\top X_3 B_3 \lambda_{k+1} + q_k^\lambda \perp \lambda_{k+1} \in K^*, $$

(where $q_k^\lambda$ is defined in (5.14)), i.e.,

$$ (6.12) \qquad \mathbb{R}_+ \times \{0\} \ni \begin{pmatrix} \lambda_{1,k+1} \\ \lambda_{2,k+1} \end{pmatrix} \perp \begin{pmatrix} 0 \\ \lambda_{2,k+1} \end{pmatrix} + q_k^\lambda \in \mathbb{R}_+ \times \mathbb{R}. $$

Therefore $\lambda_{2,k+1} = 0$, and $\lambda_{1,k+1} = 0$ is always a solution (notice that the external excitation is constrained to be non negative to comply with the system's constraints). This result agrees with the numerical simulations in Figure 6b.

<figure><img src="image_8.png" /><figcaption>Figure 6: Implicit Euler numerical simulation of (6.9) with $h = 5 \cdot 10^{-6}$s. The solutions $x_1(t)$ are not displayed as they are independent of $z(t)$, $\lambda_1(t)$. The variable $\lambda_2(t)$ trivially equals 0.</figcaption></figure>

## 6.3 THIRD EXAMPLE

Let us consider now an academic toy-example directly in the form of (2.2), whose dynamics are given as:

$$ (6.13) \qquad \left\{ \begin{array}{l} \dot{x}_1(t) = -2x_1(t) + 2\lambda_1(t) - \lambda_2(t) \\ \dot{x}_2(t) = z(t) \\ 0 = x_2(t) + \lambda_1(t) - 2\lambda_2(t) \\ 0 \le \begin{pmatrix} \lambda_1(t) \\ \lambda_2(t) \end{pmatrix} \perp \begin{pmatrix} 2 \\ -1 \end{pmatrix} x_1(t) + \begin{pmatrix} 1 \\ 3 \end{pmatrix} x_2(t) + \begin{pmatrix} -1 \\ 2 \end{pmatrix} z(t) + \begin{pmatrix} 2 & -1 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} \lambda_1(t) \\ \lambda_2(t) \end{pmatrix} \ge 0 \end{array} \right. $$

It can be verified that the transfer matrix $H(s) = D + C(sP - A)^{-1}B$ is positive real [25, Definition 2.34], and the state-space realization is minimal (see Appendix G). In fact $H(s)$ is even strongly (or extended) positive real [25, Definition 2.78]. Thus [28, Theorem 5.4 (A') (B')] applies. It is noteworthy
that since $D = \begin{pmatrix} 2 & -1 \\ 0 & 1 \end{pmatrix} > 0$, it follows that the multipliers $\lambda$ are Lipschitz continuous functions of the state $(x_1, x_2, z)$. Thus the LDCS is a DAE with Lipschitz continuous nonlinear and nonsmooth right-hand side. From (6.13) it follows that: $A_1 = -2 < 0$, $C_1 = X_1 B_1^T = (2, -1)^T$ with $X_1 = 1 > 0$, $X_3 = I_2$, $B_3 = (1 - 2) = -C_3^T$, and $\hat{D} = \begin{pmatrix} 1 & 1 \\ -3 & 7 \end{pmatrix} > 0$ with $(\hat{D} + \hat{D}^T) > 0$. It follows from [28, Theorem 5.4 (A'), Theorem 4.6] that the quadruple $(A_1, B_1, C_1, \hat{D})$ is strictly passive.

It can be checked that the set $S = \{0\}$ in Proposition 3.2 for any $h > 0$, thus there exists a unique solution $\lambda_{k+1}$ to the implicit Euler discretization of (6.13). However $C_3 \neq 0$ and Assumption 2 does not hold, preventing the convergence analysis in Section 4 and the stability analysis in Section 4.1 to apply. Consequently, the “equivalent” system (5.1) from Section 5 is analysed:

$$ (6.14) \qquad \left\{ \begin{array}{l} \left( \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & -2 \\ 0 & -2 & 4 \end{pmatrix} \begin{pmatrix} \dot{x}_1(t) \\ \dot{\lambda}_1(t) \\ \dot{\lambda}_2(t) \end{pmatrix} \right) = \begin{pmatrix} -2 & 2 & -1 \\ -2 & -1 & -1 \\ 1 & 3 & -7 \end{pmatrix} \begin{pmatrix} x_1(t) \\ \lambda_1(t) \\ \lambda_2(t) \end{pmatrix} + \begin{pmatrix} 0 \\ w_1(t) \\ w_2(t) \end{pmatrix} \\ w(t) \in -\mathcal{N}_{\mathbb{R}_+^2}(\lambda(t)) \end{array} \right. $$

Let us check if Assumption 4 (c) is verified. Using the Frobenius norm, one can obtain that $0 \le h_3 = \sqrt{60h^2 + 25 - 66h} \le 1$ if $0 < h \le 0.7$. In the context of (6.14), the LCP (5.14) becomes:

$$ (6.15) \qquad \begin{array}{l} 0 \le \begin{pmatrix} \lambda_{1,k+1} \\ \lambda_{2,k+1} \end{pmatrix} \perp \begin{pmatrix} 1 & -2 \\ -2 & 4 \end{pmatrix} \begin{pmatrix} \lambda_{1,k+1} \\ \lambda_{2,k+1} \end{pmatrix} + \begin{pmatrix} (h-1) & (h+3) \\ (2-3h) & (7h-4) \end{pmatrix} \begin{pmatrix} \lambda_{1,k} \\ \lambda_{2,k} \end{pmatrix} + \begin{pmatrix} 2hx_{1,k} \\ -hx_{1,k} \end{pmatrix} \ge 0 \\ \hfill \Leftrightarrow \\ 0 \le \begin{pmatrix} \lambda_{1,k+1} \\ \lambda_{2,k+1} \end{pmatrix} \perp \begin{pmatrix} \lambda_{1,k+1} - 2\lambda_{2,k+1} \\ 4\lambda_{2,k+1} - 2\lambda_{1,k+1} \end{pmatrix} + \begin{pmatrix} q_{1,k} \\ q_{2,k} \end{pmatrix} \ge 0. \end{array} $$

From Corollary 4.1, there exists a solution to (6.15). We can compute $\lambda_{k+1}$ in each mode of the LCP where $\lambda_{k+1} \neq 0$:

1. If $\lambda_{1,k+1} = 0$ and $\lambda_{2,k+1} \ge 0$ then $\|\lambda_{2,k+1}\| \le \frac{1}{2}\|q_{1,k}\|$.

2. If $\lambda_{2,k+1} = 0$ and $\lambda_{1,k+1} \ge 0$ then $\|\lambda_{1,k+1}\| \le \frac{1}{2}\|q_{2,k}\|$.

3. If $\lambda_{1,k+1} > 0$ and $\lambda_{2,k+1} > 0$ there is infinitely many solutions. Indeed, we obtain:

$$ (6.16) \qquad \left\{ \begin{array}{l} \lambda_{1,k+1} = 2\lambda_{2,k+1} - q_{1,k} \\ 0 = 2q_{1,k} + q_{2,k} \end{array} \right. $$

From [16, Proposition 6.1.7], there exists a minimal solution $\lambda_{\min,k+1}$. Then, (6.16)'s particular solution $\lambda_{2,k+1} = -\lambda_{1,k+1}$ leads to $\|\lambda_{\min,k+1}\| \le (1/3)\|q_{k+1}\|$. It follows that $\eta = \frac{1}{2} < 1$ for all $h > 0$, and $\eta h_3 \le 1$ for all $0 < h \le 1.4$, and finally, the results from Section 5.2 hold if small enough time steps are chosen.

Problem (6.14) is simulated using the discretization (6.15). The LCP (6.15) is solved using the numerical solver Lemke provided in the SICONOS Toolbox<fnref n="3" /> [1]. In Figures 7a and 7b the simulated variables $(x_1, x_2, z)$ and $(\lambda_1, \lambda_2)$ are displayed for time step $h = 5 \cdot 10^{-6}$s. The initial conditions are $(x_{1,0} = 5, \lambda_{1,0} = 5, \lambda_{2,0} = 0)$ and the discrete variables $x_{2,k}$ and $z_k$ are given by the discretization of (5.1) (c):

$$ (6.17) \qquad \left\{ \begin{array}{l} x_{2,k} = -B_3\lambda_k \\ z_k = -B_3 \frac{\lambda_{k+1} - \lambda_k}{h}, \end{array} \right. $$

<footnote n="3">³https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos/index.html</footnote>
with $h$ the time step. Let us notice in Figure 7a that the variable $z$ possesses a discontinuity. In Figure 8a we plot the error provided in (6.6), with respect to a reference solution taken for $h = 5 \cdot 10^{-6}$s. We confirm the order one convergence rate for $h > 0$ sufficiently small.

Finally, let us also notice that the solutions of the schemes in (3.1) and in (5.5), are not exactly the same (recall that the continuous systems (2.2) and (5.1) are “equivalent” in the sense explained after (5.3), see also Appendix H). This can be observed in Figure 8b where the differences between the two schemes are shown for $(x_1, x_2)$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Implicit Euler numerical simulation of (6.13) with $h = 5 \cdot 10^{-6}$s</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Numerical convergence studies, with error formulas (6.6), of schemes (3.1) and (5.5). The reference is computed with the scheme (5.5) and $h = 5 \cdot 10^{-6}$s.</figcaption></figure>
# 7 CONCLUSION

This article deals with the time-discretization of a class of singular linear complementarity systems (descriptor variable linear complementarity systems) which satisfy a passivity constraint. The implicit Euler method is analysed. First the one-step nonsmooth problem is studied, and conditions which guarantee its well-posedness are given. Then the approximate solutions convergence is analysed, and it is shown that under some reasonable assumptions, convergence holds. Our results are however partial in the sense that we are not able to prove the convergence of all the variables. One academic example and two examples of electrical circuits with ideal diodes (equivalently hydraulic circuits with check valves) illustrate the developments, with numerical simulations.

# APPENDIX A PRESERVATION OF POSITIVE DEFINITENESS WITH PERTURBATION

We give here an excerpt of [37, Theorem 2.11], and a corollary of it. Let us recall that for a given $M \in \mathbb{R}^{n \times n}$, $\|M\|_{2,2}$ is the induced matricial norm such that $\|M\|_{2,2} = \sigma_{\max}(M)$ (the largest singular value).

Theorem A.1 ([37]). Let $M \in \mathbb{R}^{n \times n}$ be a positive definite matrix. Then every matrix

$$A \in \{A \in \mathbb{R}^{n \times n} \mid \left\| \left( \frac{M + M^{\top}}{2} \right)^{-1} \right\|_{2,2} \|M - A\|_{2,2} < 1\}$$

is positive definite.

Corollary A.2 ([21]). Let $D = P + N$, where $D$, $P$ and $N$ are $n \times n$ real matrices, and $P > 0$, not necessarily symmetric. If

$$(A.1) \qquad \|N\|_{2,2} < \frac{1}{\left\| \left( \frac{P+P^{\top}}{2} \right)^{-1} \right\|_{2,2}}$$

then $D > 0$.

Consider the term $h_1$ in Assumptions 1 and 4. If $h < \frac{1}{\sigma_{\max}(A_1)}$ then $I_{n_1} + hA_1 > 0$ from Corollary A.2. If in addition $A_1 \ll (<) 0$ then $0 < I_{n_1} + hA_1 \ll (<) I_{n_1}$. Consider now $\tilde{A}_1 = (I_{n_1} + hA_1)^{\top}(I_{n_1} + hA_1) = I_{n_1} + h(A_1 + A_1^{\top}) + h^2 A_1^{\top} A_1 = I_{n_1} + h\tilde{A}_1$, with $\tilde{A}_1 = A_1 + A_1^{\top} + hA_1^{\top} A_1 = \tilde{A}_1^{\top}$. If $A_1 < 0$ ($\Rightarrow A_1 + A_1^{\top} < 0$), then $\tilde{A}_1 < 0$ for $h > 0$ small enough, still using Corollary A.2. Thus $\tilde{A}_1 < I_{n_1}$ for $h > 0$ small enough. Morover for $h > 0$ small enough it follows that $\tilde{A}_1 > 0$ by Corollary A.2. It is deduced that there exists $h_{\max} > 0$ such that $0 < \tilde{A}_1 < I_{n_1}$ for all $h \in (0, h_{\max})$. From [16, Lemma 8.4.1], $0 < \lambda_{\min}(\tilde{A}_1)$ and $\lambda_{\max}(\tilde{A}_1) < 1$. It is inferred that $\sigma_{\max}(I_{n_1} + hA_1) = \|I_{n_1} + hA_1\|_{2,2} = h_1 < 1$.

# APPENDIX B WELL-POSEDNESS OF VARIATIONAL INEQUALITIES

The next results use the notions of recession functions and cones, which we briefly introduce now (see [24, 53] for illustrating examples), [72, 53, 89]. Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper convex and lower semi-continuous function, we denote by $\text{dom}(f) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | f(x) < +\infty\}$ the domain of the function $f(\cdot)$. The Fenchel transform $f^*(\cdot)$ of $f(\cdot)$ is the proper, convex and lower semicontinuous function defined by

$$(B.1) \qquad (\text{for all } z \in \mathbb{R}^n) \mid f^*(z) = \sup_{x \in \text{dom}(f)} \{\langle x, z \rangle - f(x)\}.$$

The subdifferential $\partial f(x)$ of $f(\cdot)$ at $x \in \mathbb{R}^n$ is defined by

$$\partial f(x) = \{\omega \in \mathbb{R}^n | f(v) - f(x) \geq \langle \omega, v - x \rangle, \forall v \in \mathbb{R}^n\}.$$
We denote by $\text{Dom}(\partial f) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | \partial f(x) \neq \emptyset\}$ the domain of the subdifferential operator $\partial f : \mathbb{R}^n \to \mathbb{R}^n$. Recall that (see, e.g., Theorem 2, Chapter 10, Section 3 in [11]): $\text{Dom}(\partial f) \subset \text{dom}(f)$.

Let $x_0$ be any element in the domain $\text{dom}(f)$ of $f(\cdot)$, the recession function $f_\infty(\cdot)$ of $f(\cdot)$ is defined by

$$ (\text{for all } x \in \mathbb{R}^n) : f_\infty(x) = \lim_{\lambda \to +\infty} \frac{1}{\lambda} f(x_0 + \lambda x). $$

The function $f_\infty : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ is a proper convex and lower semicontinuous function which describes the asymptotic behavior of $f(\cdot)$.

Let $K \subset \mathbb{R}^n$ be a nonempty closed convex set. Let $x_0$ be any element in $K$. The recession cone of $K$ is defined by [72] [89, Definition 1.11]:

$$ K_{\infty} = \bigcap_{\lambda > 0} \frac{1}{\lambda} (K - x_0) = \{u \in \mathbb{R}^n | x + \lambda u \in K, \forall \lambda \ge 0, \forall x \in K\}. $$

The set $K_\infty$ is a nonempty closed convex cone that is described in terms of the directions which recede from $K$. The indicator function of a set $K \subset \mathbb{R}^n$ is $\Psi_K(x) = 0$ if $x \in K$, $\Psi_K(x) = +\infty$ if $x \notin K$. If $K$ is closed non empty convex, we have $\partial \Psi_K(x) = N_K(x)$, the so-called normal cone to $K$ at $x$, defined as $N_K(x) = \{v \in \mathbb{R}^n | v^\top(s-x) \le 0 \text{ for all } s \in K\}$. When $K$ is finitely represented, i.e., $K = \{x \in \mathbb{R}^n | k_i(x) \ge 0, 1 \le i \le m\}$, and if the functions $k_i(\cdot)$ satisfy some constraint qualification (like, independency, or extensions like the MFCQ), then $N_K(x)$ is generated by the outwards normals at the active constraints $k_i(x) = 0$, i.e., $N_K(x) = \{v \in \mathbb{R}^n | v = -\lambda_i \nabla k_i(x), k_i(x) = 0, \lambda_i \ge 0\}$.

Let us here recall some important properties of the recession function and recession cone (see e.g., [17, Proposition 1.4.8]):

**Proposition B.1.** *The following statements hold:*

a) Let $f_1 : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ and $f_2 : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be two proper, convex and lower semicontinuous functions. Suppose that $f_1 + f_2$ is proper. Then for all $x \in \mathbb{R}^n$: $(f_1 + f_2)_\infty(x) = (f_1)_\infty(x) + (f_2)_\infty(x)$.

b) Let $f : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper, convex and lower semicontinuous function and let $K$ be a nonempty closed convex set, such that $f + \Psi_K$ is proper (equivalently $\text{dom}(f) \cap K$ is non empty). Then for all $x \in \mathbb{R}^n$: $(f + \Psi_K)_\infty(x) = f_\infty(x) + (\Psi_K)_\infty(x)$.

c) Let $K \subset \mathbb{R}^n$ be a nonempty, closed and convex set. Then for all $x \in \mathbb{R}^n$: $(\Psi_K)_\infty(x) = \Psi_{K_\infty}(x)$. Moreover for all $x \in K$ and $e \in K_\infty$: $x + e \in K$.

d) If $K \subset \mathbb{R}^n$ is a nonempty closed and convex cone, then $K_\infty = K$.

e) Let $K = P(A, b) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | Ax \ge b\}$ for $A \in \mathbb{R}^{m \times n}$ and $b \in \mathbb{R}^m$. If $K \neq \emptyset$ then $K_\infty = P(A, 0) = \{x \in \mathbb{R}^n | Ax \ge 0\}$.

f) $K \subset \mathbb{R}^n$ is a non-empty closed convex bounded set if and only if $K_\infty = \{0_n\}$.

g) Let $f : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper, convex and lower semicontinuous function. Then $\text{epi}(f_\infty) = (\text{epi}(f))_\infty$.

Sets as in item e) are called *H-polyhedra*, and there is an equivalence between sets $P(A, 0)$ and finitely generated convex cones [89, Theorem 1.3]. Let us now concatenate [5, Theorem 3, Corollaries 3 and 4]. They concern variational inequalities (VI s) of the form: Find $u \in \mathbb{R}^n$ such that

$$ (B.2) \qquad \langle \mathbf{M}u + \mathbf{q}, v - u \rangle + \varphi(v) - \varphi(u) \ge 0, \text{ for all } v \in \mathbb{R}^n $$
where $\mathbf{M} \in \mathbb{R}^{n \times n}$ is a real matrix, $\mathbf{q} \in \mathbb{R}^n$ a vector and $\varphi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ a proper convex and lower semicontinuous function. The VI in (B.2) is equivalent to the inclusion $\mathbf{M}\mathbf{u} + \mathbf{q} \in -\partial\varphi(\mathbf{u}) \Leftrightarrow \mathbf{u} \in (\mathbf{M} + \partial\varphi)^{-1}(-\mathbf{q})$.

The problem in (B.2) is denoted as $VI(\mathbf{M}, \mathbf{q}, \varphi)$ in the next proposition. We also set:

$$ (B.3) \qquad \mathcal{K}(\mathbf{M}, \varphi) = \{x \in \mathbb{R}^n | \mathbf{M}x \in (\mathrm{dom}(\varphi_\infty))^*\}. $$

Note that $(\mathrm{dom}(\varphi_\infty))^*$ is the dual cone of the domain of the recession function $\varphi_\infty$ while $(\mathrm{dom}(\varphi))^*_{\infty}$ (that we may denote also as $\mathrm{dom}(\varphi)^*_{\infty}$) is the recession cone of $\mathrm{dom}(\varphi)$.

**Proposition B.2** ([5]). Let $\varphi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper, convex and lower semicontinuous function with closed domain, $\mathbf{M} \in \mathbb{R}^{n \times n}$, and suppose that $\mathbf{M} \ge 0$ (not necessarily symmetric).

a) If $(\mathrm{dom}(\varphi))^*_{\infty} \cap \mathrm{ker}\{\mathbf{M} + \mathbf{M}^\top\} \cap \mathcal{K}(\mathbf{M}, \varphi) = \{0\}$ then for each $\mathbf{q} \in \mathbb{R}^n$, problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at least one solution.

b) Suppose that $(\mathrm{dom}(\varphi))^*_{\infty} \cap \mathrm{ker}\{\mathbf{M} + \mathbf{M}^\top\} \cap \mathcal{K}(\mathbf{M}, \varphi) \neq \{0\}$. If there exists $x_0 \in \mathrm{dom}(\varphi)$ such that

$$ (B.4) \qquad \langle \mathbf{q} - \mathbf{M}^{\top} x_0, v \rangle + \varphi_{\infty}(v) > 0, \quad \forall v \in (\mathrm{dom}(\varphi))^*_{\infty} \cap \mathrm{ker}\{\mathbf{M} + \mathbf{M}^{\top}\} \cap \mathcal{K}(\mathbf{M}, \varphi), \quad v \neq 0, $$
then problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at least one solution.

b') If $\mathbf{M} = \mathbf{M}^\top$ then one can take $x_0 = 0$ in b).

c) If $u_1$ and $u_2$ denote two solutions of problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ then $u_1 - u_2 \in \mathrm{ker}\{\mathbf{M} + \mathbf{M}^\top\}$.

d) If $\mathbf{M} = \mathbf{M}^\top$ and $u_1$ and $u_2$ denote two solutions of problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$, then $\langle \mathbf{q}, u_1 - u_2 \rangle = \varphi(u_2) - \varphi(u_1)$.

e) If $\mathbf{M} = \mathbf{M}^\top$ and $\varphi(x + z) = \varphi(x)$ for all $x \in \mathrm{dom}(\varphi)$ and $z \in \mathrm{ker}\{\mathbf{M}\}$ and $\langle \mathbf{q}, e \rangle \neq 0$ for all $e \in \mathrm{ker}\{\mathbf{M}\}$, $e \neq 0$, then problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at most one solution.

f) If $\mathbf{M} = \mathbf{M}^\top$, then $u$ is a solution of $VI(\mathbf{M}, \mathbf{q}, \varphi)$ if and only if it is a solution of the optimization problem $\min_{x \in \mathbb{R}^n} \frac{1}{2} x^\top \mathbf{M} x + \langle \mathbf{q}, x \rangle + \varphi(x)$.

Notice that the function $\varphi(\cdot)$ will never be strictly convex in our case (it is an indicator function) so that the strict convexity argument of [5, Theorem 5] which applies when $\mathbf{M}$ is a $P_0$-matrix never holds. The study of VIs as in (B.2) can be traced back to [76].

## APPENDIX C SOME CONVEX ANALYSIS AND COMPLEMENTARITY THEORY TOOLS

If $K \subset \mathbb{R}^n$ is a set, then $K^* = \{z \in \mathbb{R}^n | \langle z, x \rangle \ge 0 \text{ for all } x \in K\}$ is its dual cone. Its closure is denoted $\bar{K}$. Let $K$ be a nonempty closed convex cone, then:

$$ (c.1) \qquad K^* \ni x \perp y \in K \iff x \in -\mathcal{N}_K(y) \iff y \in -\mathcal{N}_{K^*}(x). $$

Let $M = M^\top > 0$, $x$ and $y$ two vectors, then

$$ (c.2) \qquad M(x - y) \in -\mathcal{N}_K(x) \Leftrightarrow x = \mathrm{proj}_M[K; y] \Leftrightarrow x = \min_{z \in K} \frac{1}{2}(z - y)^\top M (z - y). $$

We note that this is a particular case of (B.2), so that Proposition B.2 can be considered as the characterization of a generalized projection operator $VI(\mathbf{M}, \mathbf{q}, \varphi)$. The first equivalence is [15, Proposition 6.46].
# APPENDIX D PASSIVE DESCRIPTOR VARIABLE SYSTEMS

Dissipative systems have been formalized rigorously by J.C. Willems in his two seminal articles [87, 88]. Passive systems are a particular case of dissipative systems. A quadruple $(A, B, C, D)$ is said passive if the next LMI is satisfied:

$$ (D.1) \qquad \begin{pmatrix} -XA - A^\top X & -XB + C^\top \\ -B^\top X + C & D + D^\top \end{pmatrix} \ge 0, \quad X = X^\top \ge 0. $$

It is said strictly passive if the inequalities are satisfied with $> 0$. This is not to be confused with the strict state passivity [25, Definition 4.54], since strict passivity implies strict state passivity, but not the inverse in general. In fact strict passivity is directly related to strongly SPR transfer matrices [25, p.62 and section 3.12.2], while SPR transfer matrices with minimal realisations are strictly state passive [25, Theorem 4.73] [66], i.e., $-XA - A^\top X + \mu X \ge 0$ for some $\mu > 0$, and one may have $D + D^\top \ge 0$ but not $> 0$.

Positive real transfer matrices [25, Definition 2.34] with a state space realization represented by quintuples $(P, A, B, C, D)$ which are minimal (see Appendix G for conditions of minimality of a state-space realization for descriptor systems) can be rewritten as in (2.2) [28, 51], and they have the positive semidefinite storage function $V(x_1, x_2) = \frac{1}{2}(x_1^\top X_1 x_1 + x_2^\top X_2 x_2)$, with supply rate $\lambda^\top w$. Indeed, by differentiating $V(x_1(t), x_2(t))$ along the trajectories of (2.2), adding $\pm \lambda^\top w$ and completing the squares, one can compute that:

$$ (D.2) \qquad \dot{V}(t) = \frac{1}{2}(x_1^\top, \lambda^\top) \underbrace{\begin{pmatrix} X_1 A_1 + A_1^\top X_1 & X_1 B_1 - C_1^\top \\ B_1^\top X_1 - C_1 & -(\hat{D} + \hat{D}^\top) \end{pmatrix}}_{\triangleq Q_1} \begin{pmatrix} x_1 \\ \lambda \end{pmatrix} + \lambda^\top(t) w(t), $$

where $X_1 = X_1^\top \ge 0$, $Q_1 \le 0$, and if the system is strongly SPR then $X_1 > 0$, $Q_1 < 0$ [28, Theorem 5.4]. The equality in (D.2) is an infinitesimal dissipation equality, showing that the system is passive [25, Definition 2.1, Theorem 4.53]. Such systems possess particular structural properties [30, 33] [25, section 3.8]. In the SSPR case it follows from (D.2) and the complementarity conditions that $x_1(\cdot)$ and $x_2^\top X_3 x_2(\cdot)$ are bounded for any bounded initial condition.

Let us recall some definitions, where $A^*$ denotes the conjugate transpose of $A$ (see [25, Chapter 2] for details):

**Definition D.1.** A transfer matrix $H(s) \in \mathbb{C}^{m \times m}$, $s \in \mathbb{C}$, is positive real (PR) if: 1) $H(s)$ has no poles in $\mathbf{Re}(H(s)) > 0$, 2) $H(s)$ is real for positive real $s$, 3) $H(s) + H^*(s) \ge 0$ for all $\mathbf{Re}(s) > 0$. Let $H(s)$ be regular, then it is strictly positive real (SPR) if $H(s - \epsilon)$ is PR for some $\epsilon > 0$. It is strongly SPR (SSPR) if 1) $H(s)$ is analytic in $\mathbf{Re}(s) \ge 0$, 2) $\mathbf{Re}(H(j\omega)) \ge \delta > 0$ for all $\omega \in [-\infty, +\infty]$ and some real $\delta$.

# APPENDIX E GRONWALL-BACHPATTE INEQUALITY [13, 47]

Let us consider a positive real numbers sequence $\{x_k\}$ such that $x_k \le \alpha_k + \sum_{i=0}^{k-1} \beta_i x_i$, where $\{\alpha_k\}$ and $\{\beta_k\}$ are two real numbers sequences with $\{\beta_k\}$ positive. Then one has

$$ x_k \le \alpha_k + \sum_{i=0}^{k-1} \beta_i \alpha_i \prod_{j=i+1}^{k} (1 + \beta_j). $$

# APPENDIX F CALCULATION OF $\beta_i$

The starting point is the term $\eta \|C_1\| \|B_1\| \sum_{i=1}^{k} h_1^{k-i} \sum_{j=0}^{i-1} (\eta h_3)^{i-1-j} \|x_{1,j+1}\|$ in (5.19). Let us consider the factors of the terms $\|x_{1,j+1}\|$ separately:
• Factor of $\|x_{1,1}\|$ (all terms with $j=0$): $h_1^{k-1} + h_1^{k-2}(\eta h_3) + h_1^{k-3}(\eta h_3)^2 + \dots + h_1(\eta h_3)^{k-2} + (\eta h_3)^{k-1} = \sum_{i=1}^k h_1^{k-i}(\eta h_3)^{i-1}$.

• Factor of $\|x_{1,2}\|$ (all terms with $j=1$): $h_1^{k-2} + h_1^{k-3}(\eta h_3) + h_1^{k-4}(\eta h_3)^2 + \dots + h_1(\eta h_3)^{k-3} + (\eta h_3)^{k-2} = \sum_{i=2}^k h_1^{k-i}(\eta h_3)^{i-2}$.

• ...

• Factor of $\|x_{1,k-1}\|$ (all terms with $j=k-2$): $h_1 + (\eta h_3) = \sum_{i=k-1}^k h_1^{k-i}(\eta h_3)^{i-k+1}$.

• Factor of $\|x_{1,k}\|$ (all terms with $j=k$): $1 = \sum_{i=k}^k h_1^{k-i}(\eta h_3)^{i-k}$.

Therefore one infers (5.22).

## APPENDIX G MINIMALITY OF STATE SPACE REALIZATIONS

Let us report the results stated in [28, Theorem 3.1] [44, Theorems 2-2.1, 2-3.1, 2-6.3].

**Theorem G.1.** Let $H(s) = C(sP - A)^{-1}B + D$, $s \in \mathbb{C}$, be a rational transfer function where $P$ and $A$ are $n \times n$ matrices. Then $P\dot{x}(t) = Ax(t) + Bu(t)$, $y(t) = Cx(t) + Du(t)$ is a minimal realization of $H(s)$ if, and only if, the following conditions are satisfied:

• $\text{rank}(A - sP B) = n$ for all $s \in \mathbb{C}$,

• $\text{rank}(P B) = n$,

• $\text{rank}(A^\top - sP^\top C^\top) = n$,

• $\text{rank}(P^\top C^\top) = n$,

• A $\text{Ker}(P) \subseteq \text{Im}(P)$.

## APPENDIX H DIFFERENCES BETWEEN (3.1) AND IN (5.5)

Let us consider the discrete scheme (5.5) with $E(t) = 0$ and $X_1$ full rank, this is rewritten as:

$$ (H.1a) \qquad x_{1,k+1} = x_{1,k} + hA_1x_{1,k} + hB_1\lambda_k $$

$$ (H.1b) \qquad B_3^\top X_3 B_3 \lambda_{k+1} = B_3^\top X_3 B_3 \lambda_k - hC_1x_{1,k} - h\hat{D}\lambda_k + w_{k+1} $$

$$ (H.1c) \qquad 0 \le \lambda_{k+1} \perp w_{k+1} = C_1x_{1,k+1} + C_2x_{2,k+1} + C_3z_{k+1} + D\lambda_{k+1} \ge 0 $$

$$ (H.1d) \qquad x_{2,k} = -B_3\lambda_k $$

$$ (H.1e) \qquad z_k = -B_2\lambda_k - B_3 \frac{\lambda_{k+1} - \lambda_k}{h}. $$

This is the numerical scheme with which the examples have been simulated. Combining (H.1) (d) and (e), gives:

$$ (H.2) \qquad x_{2,k+1} = x_{2,k} + hz_k + hB_2\lambda_k. $$

Notice that this is the choice made in (6.17). The equations in (H.1) (a) and (H.2) are not equal to their counterparts in (3.1), because $\lambda_k$ and $z_k$ are used in (H.1) (a) and (H.2) instead of $\lambda_{k+1}$ and $z_{k+1}$ in (3.1). This explains that the numerical data obtained from both schemes, differ. This also explains why the convergence results obtained for both schemes, are not the same neither.
## REFERENCES

[1] V. Acary, O. Bonnefon, M. Bremond, O. Huber, F. Pérignon, and S. Sinclair, An introduction to *SICONOS* (2019). RT 0340, INRIA, https://hal.inria.fr/inria-00162911v3/document.

[2] V. Acary, O. Bonnefon, and B. Brogliato, *Nonsmooth Modeling and Simulation for Switched Circuits*, volume 69 of Lecture Notes in Electrical Engineering, Springer Science+Business Media, 2011. Erratum at https://hal.inria.fr/hal-01311078.

[3] V. Acary and B. Brogliato, *Numerical Methods for Nonsmooth Dynamical Systems. Applications in Mechanics and Electronics*, volume 35 of Lecture Notes in Applied and Computational Mechanics, Springer Verlag, Berlin Heidelberg, 2008.

[4] V. Acary, H. de Jong, and B. Brogliato, Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems, *Physica D: Nonlinear Phenomena* 269 (2014), 103-119.

[5] K. Addi, B. Brogliato, and D. Goeleven, A qualitative mathematical analysis of a class of linear variational inequabilities via semi-complementarity problems: applications in electronics, *Mathematical Programming A* 126 (2011), 31-67.

[6] S. Adly, T. Haddad, and L. Thibault, Convex sweeping process in the framework of measure differential inclusions and evolution variational inequalities, *Mathematical Programming A* 148 (2014), 5-47.

[7] S. Adly, A. Hantoute, and B. K. Le, Nonsmooth Lur'e dynamical systems in Hilbert spaces, *Set-Valued and Variational Analysis* 24 (2016), 13-35.

[8] S. Adly, A. Hantoute, and B. K. Le, Maximal monotonicity and cyclic monotonicity arising in nonsmooth Lur'e dynamical systems, *Journal of Mathematical Analysis and Applications* 448 (2017), 691-706.

[9] S. Adly and B. K. Le, On semicoercive sweeping process with velocity constraint, *Optimization Letters* 12 (2018), 831-843.

[10] H. Attouch and J. Peypouquet, Convergence of inertial dynamics and proximal algorithms governed by maximally monotone operators, *Mathematical Programming* 174 (2019), 391-432.

[11] J. P. Aubin, *Applied Functional Analysis*, John Wiley & Sons, New York, 1979.

[12] J. P. Aubin and A. Cellina, *Differential Inclusions. Set-Valued Maps and Viability Theory*, Springer-Verlag, Berlin Heidelberg, 1984.

[13] B. G. Bachpatte, On the discrete generalisations of Gronwall's inequality, *J. Indian Math. Soc.* 37 (1987), 147-156.

[14] J. Bastien and M. Schatzman, Numerical precision for differential inclusions with uniqueness, *ESAIM Mathematical Modelling and Numerical Analysis* 36 (2002), 427-460.

[15] H. H. Bauschke and P. L. Combettes, *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*, CMS Books in Mathematics. Canadian Mathematical Society. Société Mathématique du Canada, Springer Science+Business Media, 2011.

[16] D. S. Bernstein, *Matrix Mathematics. Theory, Facts and Formulas*, Princeton University Press, 2nd edition, 2009.
[17] D. P. Bertsekas, *Convex Optimization Theory*, Athena Scientific, Belmont, 2009.

[18] H. Brezis, *Analyse Fonctionnelle Théorie et Applications*, Mathématiques Appliquées pour la Maitrise, Masson, Paris, France, 1983.

[19] H. Brézis and A. Haraux, Image d'une somme d'opérateurs monotones et applications, *Israel Journal of Mathematics* 23 (1976), 165-186.

[20] B. Brogliato, Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings, *Systems and Control Letters* 51 (2004), 343-353.

[21] B. Brogliato, Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction, *Multibody System Dynamics* 32 (2014), 175-216.

[22] B. Brogliato and D. Goeleven, Well-posedness, stability and invariance results for a class of multivalued Lur'e dynamical systems, *Nonlinear Analysis. Theory, Methods and Applications* 74 (2011), 195-212.

[23] B. Brogliato and D. Goeleven, Existence, uniqueness of solutions and stability of nonsmooth multivalued Lur'e dynamical systems, *Journal of Convex Analysis* 20 (2013), 881-900.

[24] B. Brogliato and D. Goeleven, Singular mass matrix and redundant constraints in unilaterally constrained Lagrangian and Hamiltonian systems, *Multibody System Dynamics* 35 (2015), 39-61.

[25] B. Brogliato, R. Lozano, B. Maschke, and O. Egeland, *Dissipative Systems Analysis and Control. Theory and Applications*, Communications and Control Engineering, Springer Nature Switzerland AG, Cham, CH, third edition, 2020.

[26] B. Brogliato and A. Tanwani, Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability, *SIAM Review* 62 (2020), 3-129.

[27] B. Brogliato and L. Thibault, Existence and uniqueness of solutions for non-autonomous complementarity dynamical systems, *Journal of Convex Analysis* 17 (2010), 961-990.

[28] M. K. Camlibel and R. Frasca, Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems, *Systems and Control Letters* 58 (2009), 795-803.

[29] M. K. Camlibel, W. P. M. H. Heemels, and J. Schumacher, Consistency of a time-stepping method for a class of piecewise linear networks, *IEEE Transactions on Circuits and Systems I* 49 (2002), 349-357.

[30] M. K. Camlibel, W. P. M. H. Heemels, and J. M. Schumacher, On linear passive complementarity systems, *European Journal of Control* 8 (2002), 220-237.

[31] M. K. Camlibel, L. Iannelli, and A. Tanwani, Evolution inclusions with time-dependent maximal monotone operators (2019). arXiv:1903.10803v1.

[32] M. K. Camlibel, L. Iannelli, A. Tanwani, and S. Trenn, Differential-algebraic inclusions with maximal monotone operators, in *Proceedings of the 55th IEEE Conference on Decision and Control*, Las Vegas, USA, 2016, 610-615.

[33] M. K. Camlibel, L. Iannelli, and F. Vasca, Passivity and complementarity, *Mathematical Programming Series A* 145 (2014), 531-563.
[34] M. K. Camlibel, J. S. Pang, and J. Shen, Lyapunov stability of complementarity and extended systems, *SIAM Journal on Optimization* 17 (2006), 1056–1101.

[35] M. K. Camlibel and J. M. Schumacher, Linear passive and maximal monotone mappings, *Mathematical Programming Ser. B* 157 (2016), 397–420.

[36] T. H. Cao, G. Colombo, B. S. Mordukovich, and D. Nguyen, Optimization of fully controlled sweeping processes, *Journal of Differential Equations* 295 (2021), 138–186.

[37] X. Chen and S. Xiang, Perturbation bounds of P-matrix linear complementarity problems, *SIAM Journal on Optimization* 18 (2007), 1250–1265.

[38] R. Choupanzadeh and A. Zadehgol, Stability, causality, and passivity analysis of canonical equivalent circuits of improper rational transfer functions with real poles and residues, *IEEE Access* 8 (2020), 125149–125162.

[39] D. Cobb, On the solution of linear differential equations with singular coefficients, *Journal of Differential Equations* 46 (1982), 310–323.

[40] D. Cobb, Descriptor variable systems and optimal state regulation, *IEEE Transactions on Automatic Control* 28 (1983), 601–611.

[41] G. Colombo, R. Henrion, N. D. Hoang, and B. S. Mordukovich, Discrete approximations of a controlled sweeping process, *Set-Valued and Variational Analysis* 23 (2015), 69–86.

[42] G. Colombo, R. Henrion, N. D. Hoang, and B. S. Mordukovich, Optimal control of the sweeping process over polyhedral controlled sets, *Journal of Differential Equations* 260 (2016), 3397–3347.

[43] R. W. Cottle, J. S. Pang, and R. E. Stone, *The Linear Complementarity Problem*, Computer Science and Scientific Computing, Academic Press, 1992.

[44] L. Dai, *Singular Control Systems*, volume 118 of Lecture Notes in Control and Information Sciences, Springer Verlag, Berlin, Heidelberg, 1989.

[45] N. del Buono, C. Elia, and L. Lopez, On the equivalence between the sigmoidal approach and Utkin's approach for piecewise-linear models of gene regulatory networks, *SIAM Journal on Applied Dynamical Systems* 13 (2014), 1270–1292.

[46] A. Dontchev, T. Z. Dontchev, and I. Slavov, A Tikhonov theorem for singularly perturbed differential inclusions, *Nonlinear Analysis, Theory, Methods and Applications* 26 (1996), 1547–1554.

[47] S. S. Dragomir, *Some Gronwall Type Inequalities and Applications*, Nova Science Publishers, 2003.

[48] D. Estevez-Schwarz and C. Tischendorf, Structural analysis for electric circuits and consequences for MNA, *International Journal of Circuit Theory and Applications* 28 (2000), 131–162.

[49] F. Facchinei and J. S. Pang, *Finite-Dimensional Variational Inequalities and Complementarity Problems, Volume I*, Series in Operations Research, Springer-Verlag, New York, Berlin, Heidelberg, 2003.

[50] G. Fichera, Boundary value problems in elasticity with unilateral constraints, in *Encyclopedia of Physics*, S. Flügge (ed.), volume VIa/2 of Mechanics of Solids II, Springer, Berlin, 1972, 347–389.

[51] R. W. Freund and F. Jarre, An extension of the positive real lemma to descriptor systems, *Optimization Methods and Software* 19 (2004), 69–87.
[52] H. U. Fuchs, *The Dynamics of Heat: A Unified Approach to Thermodynamics and Heat Transfer*, Graduate Texts in Physics, Springer Science+Business Media, 2010.

[53] D. Goeleven, *Complementarity and Variational Inequalities in Electronics*, Mathematical Analysis and its Applications, Academic Press, 2017.

[54] G. Grammel, Singularly perturbed differential inclusions: an averaging approach, *Set-Valued Analysis* 4 (1996), 361–374.

[55] S. Greenhalgh, V. Acary, and B. Brogliato, On preserving dissipativity properties of linear complementarity dynamical systems with the $\theta$-method, *Numerische Mathematik* 125 (2013), 601–637.

[56] E. Guglielmino, C. Semini, H. Kogler, R. Scheidl, and D. G. Caldwell, Power hydraulics-Switched mode control of hydraulic actuation, in *Proc. IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, Taipei, Taiwan, 2010, 3031–3036.

[57] L. Han, A. Tiwari, M. K. Camlibel, and J. S. Pang, Convergence of Time-Stepping Schemes for Passive and Extended Linear Complementarity Systems, *SIAM Journal on Numerical Analysis* 47 (2009), 3768–3796.

[58] J. B. Hiriart-Urruty and C. Lemaréchal, *Fundamentals of Convex Analysis*, Grundlehren Text Editions, Springer-Verlag, Berlin Heidelberg, 2001. Corrected second printing 2004.

[59] S. Iemoto, Some results on approximate solutions of variational inequality problems for inverse strongly monotone operators, *Fixed Point Theory and Applications* (2015), 2015–86. DOI:10.1186/s13663-015-0337-y.

[60] H. Iiduka and W. Takahashi, Strong convergence theorems for nonexpansive mappings and inverse-strongly monotone mappings, *Nonlinear Analysis* 61 (2005), 341–350.

[61] A. Jourani and E. Vilches, Positively $\alpha$-far sets and existence results for generalized perturbed sweeping processes, *Journal of Convex Analysis* 23 (2016), 775–821.

[62] A. Jourani and E. Vilches, Galerkin-like method and generalized perturbed sweeping processes with nonregular sets, *SIAM Journal on Control and Optimization* 55 (2017), 2412–2436.

[63] B. K. Le, Well-posedness and nonsmooth Lyapunov pairs for state-dependent maximal monotone differential inclusions, *Optimization. A Journal of Mathematical Programming and Operations Research* (2020). https://doi.org/10.1080/02331934.2019.1686504.

[64] C. E. Lemke, Bimatrix equilibrium points and mathematical programming, *Management Science* 11 (1965), 681–689.

[65] C. E. Lemke and J. T. Howson, Jr, Equilibrium points of bimatrix games, *Journal of the Society for Industrial and Applied Mathematics* 12 (1964), 413–423.

[66] D. S. Madeira and J. Adamy, On the equivalence between strict positive realness and strict passivity of linear systems, *IEEE Transactions on Automatic Control* 61 (2016), 3091–3095.

[67] M. D. P. Monteiro-Marques, *Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction*, volume 9 of Prog. Nonlinear Differential Equations Appl., Birkhauser, Basel, CH, 1993.

[68] M. Pan and A. Plummer, Digital switched hydraulics, *Front. Mech. Eng.* 13 (2018), 225–231.
[69] J. Peypouquet and S. Sorin, Evolution Equations for Maximal Monotone Operators: Asymptotic Analysis in Continuous and Discrete Time, *Journal of Convex Analysis* 17 (2010), 1113–1163.

[70] A. Rocca, V. Acary, and B. Brogliato, Index-2 hybrid DAE: a case study with well-posedness and numerical analysis (2019). INRIA research Report, https://hal.inria.fr/hal-02381489/document.

[71] A. Rocca, V. Acary, and B. Brogliato, Index-2 hybrid DAE: a case study with well-posedness and numerical analysis, in *Proc. 21st IFAC World Congress*, Berlin, Germany, 2020.

[72] R. T. Rockafellar, *Convex Analysis*, Princeton Landmarks in Mathematics, Princeton University Press, 1970.

[73] W. Rudin, *Principles of Mathematical Analysis*, International Series in Pure and Applied Mathematics, McGraw Hill Inc., 3rd edition, 1976.

[74] J. Shen and J. S. Pang, Linear complementarity systems: Zeno states, *SIAM Journal on Control and Optimization* 44 (2005), 1040–1066.

[75] G. V. Smirnov, *Introduction to the Theory of Differential Inclusions*, volume 41 of Graduate Studies in Mathematics, American Mathematical Society, Providence, Rhode Island, 2001.

[76] G. Stampacchia, *Formes bilinéaires coercitives sur les ensembles convexes*, *C.R. Acad. Sci. Paris* 258 (1964), 4413–4416.

[77] P. Stechlinski and P. I. Barton, Dependence of solutions of nonsmooth differential-algebraic equations on parameters, *Journal of Differential Equations* 262 (2017), 2254–2285.

[78] P. Stechlinski, M. Patrascu, and B. I. Barton, Nonsmooth DAEs with Applications in Modeling Phase Changes, in *Applications of Differential-Algebraic Equations: Examples and Benchmarks*, C. S., I. A., M. V., and R. T. (eds.), Differential-Algebraic Equations Forum, Springer, 2018, 243–275.

[79] P. Stechlinski, M. Patrascu, and P. I. Barton, Nonsmooth differential-algebraic equations in chemical engineering, *Computers and Chemical Engineering* 114 (2018), 52–68.

[80] W. Sun, P. P. Khargonekar, and D. Shim, Solution to the positive real control problem for linear time-invariant systems, *IEEE Transactions on Automatic Control* 39 (1994), 2034–2046.

[81] W. Takahashi, Strong convergence theorems for maximal and inverse-strongly monotone mappings in Hilbert spaces and applications, *Journal of Optimisation Theory and Applications* 157 (2013), 781–802.

[82] A. Tanwani, B. Brogliato, and C. Prieur, Well-posedness and output regulation for implicit time-varying evolution variational inequalities, *SIAM Journal on Control and Optimization* 56 (2018), 751–781.

[83] V. I. Utkin, Sliding-mode control in discrete-time and difference systems, in *Variable Structure and Lyapunov Control*, A. S. I. Zinober (ed.), volume 193 of Lecture Notes in Control and Information Sciences, Springer Verlag, Berlin, Heidelberg, 1994, 87–107.

[84] V. Veliov, A generalization of the Tikhonov theorem for singularly perturbed differential inclusions, *Journal of Dynamical and Control Systems* 3 (1997), 291–319.

[85] E. Vilches and B. T. Nguyen, Evolution inclusions governed by time-dependent maximal monotone operators with a full domain, *Set-Valued and Variational Analysis* 28 (2020), 569–581.
[86] F. Watbled, On singular perturbations for differential inclusions on the infinite interval, *Journal of Mathematical Analysis and Applications* 310 (2005), 362-37.

[87] J. C. Willems, Dissipative dynamical systems, part I: General theory, *Archive of Rational Mechanics Analysis* 45 (1972), 321-351.

[88] J. C. Willems, Dissipative dynamical systems, part II: Linear systems with quadratic supply rates, *Archive of Rational Mechanics Analysis* 45 (1972), 352-393.

[89] G. M. Ziegler, *Lectures on Polytopes*, Graduate Texts in Mathematics, Springer Verlag, New York, USA, 1995. Revised first edition.