Bernard Brogliato
email: bernard.brogliato@inria.fr

Alexandre Rocca
email: alexandre.rocca@inria.fr

analysis of the implicit euler time-discretization of passive linear descriptor complementarity systems

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.

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 [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF] 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] [START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF], and in a more abstract framework in state-dependent switching DAEs [START_REF] Rocca | Index-2 hybrid DAE: a case study with well-posedness and numerical analysis[END_REF][START_REF] Rocca | Index-2 hybrid DAE: a case study with well-posedness and numerical analysis[END_REF], 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 [START_REF] Stechlinski | Nonsmooth DAEs with Applications in Modeling Phase Changes[END_REF][START_REF] Stechlinski | Nonsmooth differential-algebraic equations in chemical engineering[END_REF][START_REF] Stechlinski | Dependence of solutions of nonsmooth differential-algebraic equations on parameters[END_REF].

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 [START_REF] Attouch | Convergence of inertial dynamics and proximal algorithms governed by maximally monotone operators[END_REF][START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF][START_REF] Han | Convergence of Time-Stepping Schemes for Passive and Extended Linear Complementarity Systems[END_REF][START_REF] Greenhalgh | On preserving dissipativity properties of linear complementarity dynamical systems with the 𝜃 -method[END_REF][START_REF] Acary | Nonsmooth Modeling and Simulation for Switched Circuits[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Monteiro-Marques | Differential Inclusions in Nonsmooth Mechanical Problems[END_REF][START_REF] Peypouquet | Evolution Equations for Maximal Monotone Operators: Asymptotic Analysis in Continuous and Discrete Time[END_REF], 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.

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 discretetime 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 𝑥 ∈ ℝ 𝑛 and any matrix 𝑀 ∈ ℝ 𝑚×𝑛 , ∥𝑥 ∥ is the Euclidean norm and ∥𝑀 ∥ is the Frobenius norm, which are compatible norms [START_REF] Bernstein | Matrix Mathematics. Theory, Facts and Formulas[END_REF]Proposition 9.3.5], i.e., ∥𝑀𝑥 ∥ ≤ ∥𝑀 ∥ ∥𝑥 ∥. Let 𝑀 ∈ ℝ 𝑛×𝑚 , then Im(𝑀) is its range, Ker(𝑀) is its null space. We use ⟨𝑥, 𝑦⟩ = 𝑥 ⊤ 𝑦, so ⟨𝑥, 𝑥⟩ = ∥𝑥 ∥ 2 . Positive definite matrix: 𝑀 ≻ 0 if 𝑥 ⊤ 𝑀𝑥 > 0 for all 𝑥 ≠ 0, positive semidefinite matrix: 𝑀 ≽ 0 if 𝑥 ⊤ 𝑀𝑥 ≥ 0 for all 𝑥 (such 𝑀 is not necessarily symmetric). The maximum singular value is denoted as 𝜎 max (𝑀), and the minimum and maximum eigenvalues as 𝜆 min (𝑀) and 𝜆 max (𝑀). 𝑀 •𝑖 denotes the 𝑖-th column of 𝑀. More mathematical results and definitions are given in Appendices a to g.

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) 𝑃 𝑥 (𝑡) = 𝐴𝑥 (𝑡) + 𝐵𝜆(𝑡) + 𝐸 (𝑡), a.e. 𝑡 ≥ 0 0 ≤ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 𝐶𝑥 (𝑡) + 𝐷𝜆(𝑡) + 𝐹 (𝑡) ≥ 0, for all 𝑡 ≥ 0 with 𝑥 (𝑡) ∈ ℝ 𝑛 , 𝜆(𝑡) ∈ ℝ 𝑚 , 𝐴, 𝐵, 𝐶, 𝐷 constant matrices of appropriate dimensions, 𝑃 ∈ ℝ 𝑛×𝑛 has rank 𝑝 < 𝑛. Dissipativity is a fundamental property in Systems and Control [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF][START_REF] Willems | Dissipative dynamical systems, part I: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems, part II: Linear systems with quadratic supply rates[END_REF]. The material in this article relies on a result in [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF][START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF][START_REF] Freund | An extension of the positive real lemma to descriptor systems[END_REF] on a special Weierstrass form for passive descriptor variable systems [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF]Section 3.1.7] which possess a minimal state-space realization [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Theorem 3.1]. The passivity of the quintuple (𝑃, 𝐴, 𝐵, 𝐶, 𝐷) means the passivity of the operator 𝜆 ↦ → 𝑤, see Appendix d. Assuming that 𝐸 (𝑡) = 0 and 𝐹 (𝑡) = 0, passivity and minimality, this special form writes as (see [START_REF] Freund | An extension of the positive real lemma to descriptor systems[END_REF] [32, Equation ( 16)] [28, Proposition A.3]):

(2.2)

             (𝑎) 𝑥 1 (𝑡) = 𝐴 1 𝑥 1 (𝑡) + 𝐵 1 𝜆(𝑡) (𝑏) 𝑥 2 (𝑡) = 𝑧 (𝑡) + 𝐵 2 𝜆(𝑡) (𝑐) 0 = 𝑥 2 (𝑡) + 𝐵 3 𝜆(𝑡) (𝑑) 0 ≤ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 𝐶 1 𝑥 1 (𝑡) + 𝐶 2 𝑥 2 (𝑡) + 𝐶 3 𝑧 (𝑡) + 𝐷𝜆(𝑡) ≥ 0,
where the quadruple

(𝐴 1 , 𝐵 1 , 𝐶 1 , D) is itself passive, D Δ = 𝐷 -𝐶 3 𝐵 2 -𝐶 2 𝐵 3 , 𝑥 1 (𝑡) ∈ ℝ 𝑛 1 , 𝑥 2 (𝑡) ∈ ℝ 𝑛 2 , 𝑧 (
𝑡) ∈ ℝ 𝑛 2 , 𝑛 = 𝑛 1 + 2𝑛 2 , 𝜆(𝑡) ∈ ℝ 𝑚 , and there exists 𝑋 3 = 𝑋 ⊤ 3 ≽ 0 such that 𝐵 ⊤ 3 𝑋 3 = -𝐶 3 . The matrix 𝑋 3 is unique when the state-space realization is minimal [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Theorem 5.4]. The LMI in (d.2) holds. Due to the passivity and the complementarity conditions which imply 𝜆(𝑡) ⊤ 𝑤 (𝑡) = 0 for all times, it follows that for all 𝑡 ≥ 0 one obtains 𝑉 (𝑥 1 (𝑡), 𝑥 2 (𝑡)) [START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF] and [28, Proof of Theorem 5.1], with 𝑋 1 ≻ 0 if the system is strongly SPR and minimal [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Theorem 5.4], see Appendix d. In the latter case 𝑥 1 (•) is bounded, 𝐴 1 is Hurwitz 1 , D ≻ 0 and 𝐵 3 is full row rank (⇒ 𝑛 2 ≤ 𝑚) [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Proposition A.4].

Δ = 1 2 (𝑥 ⊤ 1 𝑋 1 𝑥 1 + 𝑥 ⊤ 2 𝑋 3 𝑥 2 ) (𝑡) ≤ 𝑉 (𝑥 1 (0), 𝑥 2 (0)),
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 𝜆(𝑡) = 0) to recover the variable 𝑧 (𝑡). 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 [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF]) have been studied [START_REF] Vilches | Evolution inclusions governed by time-dependent maximal monotone operators with a full domain[END_REF][START_REF] Adly | Convex sweeping process in the framework of measure differential inclusions and evolution variational inequalities[END_REF][START_REF] Adly | On semicoercive sweeping process with velocity constraint[END_REF], 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) 𝑃 𝑥 (𝑡) -𝐴𝑥 (𝑡) -𝐸 (𝑡) ∈ 𝐵(N 𝐾 ★ + 𝐷) -1 (-𝐶𝑥 (𝑡) -𝐹 (𝑡))

with 𝐾 = 𝐾 ★ = ℝ 𝑚 + , 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)

𝐼 𝑛 1 0 0 0 𝐼 𝑛 2 0 0 0 0 𝑥 1 (𝑡) 𝑥 2 (𝑡) 𝑧 (𝑡) - 𝐴 1 0 0 0 0 𝐼 𝑛 2 0 𝐼 𝑛 2 0 𝑥 1 (𝑡) 𝑥 2 (𝑡) 𝑧 (𝑡) ∈ 𝐵 1 𝐵 2 𝐵 3 (N 𝐾 ★ + 𝐷) -1 (-𝑤 (𝑡)).
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 𝐷 = 0, both (2.4) and (2.3) can be rewritten, under some assumptions, as a singular first-order sweeping process (FOSwP [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF]). The optimal control of FOSwP as treated in [START_REF] Colombo | Optimal control of the sweeping process over polyhedral controlled sets[END_REF][START_REF] Colombo | Discrete approximations of a controlled sweeping process[END_REF][START_REF] Cao | Optimization of fully controlled sweeping processes[END_REF], 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.

the euler implicit scheme and the osnsp

In the sequel of the article, we let 𝑡 ∈ [0,𝑇 ], 0 < 𝑇 < +∞, ℎ = 𝑇 𝑁 , 𝑁 ∈ ℕ \ {0}, 𝑡 𝑘 = 𝑘ℎ, 𝑘 ∈ {0, 1, . . . , 𝑁 }, and 𝑓 𝑘 = 𝑓 (𝑡 𝑘 ). The studied implicit Euler time-discretization of (2.2) is as follows:

(3.1)              𝑥 1,𝑘+1 = 𝑥 1,𝑘 + ℎ𝐴 1 𝑥 1,𝑘 + ℎ𝐵 1 𝜆 𝑘+1 𝑥 2,𝑘+1 = 𝑥 2,𝑘 + ℎ𝑧 𝑘+1 + ℎ𝐵 2 𝜆 𝑘+1 0 = 𝑥 2,𝑘+1 + 𝐵 3 𝜆 𝑘+1 0 ≤ 𝜆 𝑘+1 ⊥ 𝑤 𝑘+1 = 𝐶 1 𝑥 1,𝑘+1 + 𝐶 2 𝑥 2,𝑘+1 + 𝐶 3 𝑧 𝑘+1 + 𝐷𝜆 𝑘+1 ≥ 0,
for all 𝑘 ∈ {0, 1, . . . , 𝑁 -1}, where 𝐴 1 𝑥 1,𝑘 could be replaced by 𝐴 1 𝑥 1,𝑘+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 ℎ𝑤 𝑘+1 is used instead of 𝑤 𝑘+1 ) as an LCP:

(3.2) 0 ≤ 𝜆 𝑘+1 ⊥ Dℎ 𝜆 𝑘+1 + ℎ𝛽 𝑘 ≥ 0 ⇔ Dℎ 𝜆 𝑘+1 + ℎ𝛽 𝑘 ∈ -N ℝ 𝑚 + (𝜆 𝑘+1 ) ⇔ 𝜆 𝑘+1 ∈ ( Dℎ + N ℝ 𝑚 + ) -1 (-ℎ𝛽 𝑘 ), with Dℎ Δ = ℎ D + ℎ 2 𝐶 1 𝐵 1 + 𝐵 ⊤ 3 𝑋 3 𝐵 3 , 𝛽 𝑘 Δ = 𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 -1 ℎ 𝐶 3 𝑥 2,𝑘 . It is noteworthy that 𝑧 𝑘+1 is eliminated and does not appear in 𝛽 𝑘 . The inclusion (3.2) is a VI as (b.2) with 𝜑 (•) = Ψ ℝ 𝑚 + (•).
Obviously the well-posedness of the OSNSP depends on the properties of Dℎ . From the passivity we know that D ≽ 0 and 𝐵 ⊤ 3 𝑋 3 𝐵 3 ≽ 0. Using Theorem a.1 and Corollary a.2 it is possible to characterize the set of admissible ℎ𝐶 1 𝐵 1 that do not destroy the positive definiteness of Dℎ in case 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 [START_REF] Greenhalgh | On preserving dissipativity properties of linear complementarity dynamical systems with the 𝜃 -method[END_REF] for the case of LCS discretized with a 𝜃 -method.

As alluded to in Remark 2.2, the cone ℝ 𝑚 + can be replaced in (2.1), (2.2) and (3.1) by any closed non empty convex cone 𝐾 ⊂ ℝ 𝑚 and its dual cone 𝐾 ★ . Then the complementarity conditions become

𝐾 ★ ∋ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 𝐶 1 𝑥 1 (𝑡) + 𝐶 2 𝑥 2 (𝑡) + 𝐶 3 𝑧 (𝑡) + 𝐷𝜆(𝑡) ∈ 𝐾, yielding the inclusion (3.3) Dℎ 𝜆 𝑘+1 + ℎ𝛽 𝑘 ∈ -N 𝐾 ★ (𝜆 𝑘+1 ) ⇔ 𝜆 𝑘+1 ∈ ( Dℎ + N 𝐾 ★ ) -1 (-ℎ𝛽 𝑘 ),
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 𝛽 𝑘 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 ℎ > 0 and 𝑘 ≥ 0 be given. Assume that Dℎ ≽ 0, not necessarily symmetric, and let 𝐾 ⊂ ℝ 𝑚 be a non empty closed convex cone (where the complementarity conditions in (2.2) are extended to

𝐾 ★ ∋ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 𝐶 1 𝑥 1 (𝑡) + 𝐶 2 𝑥 2 (𝑡) + 𝐶 3 𝑧 (𝑡) + 𝐷𝜆(𝑡) ∈ 𝐾). Let S = 𝐾 ★ ∩ Ker( Dℎ + D⊤ ℎ ) ∩ {𝑤 ∈ ℝ 𝑚 | Dℎ 𝑤 ∈ 𝐾 }. Then: 1. If S = {0}, for each ℎ𝛽 𝑘 the LCCP in (3.
3) has at least a solution.

2. If S ≠ {0}, and there exists 𝜉 0 ∈ ℝ 𝑚 such that ⟨ℎ𝛽 𝑘 -D⊤ ℎ 𝜉 0 , 𝑣⟩ + Ψ 𝐾 ★ (𝑣) > 0 for all 𝑣 ≠ 0, 𝑣 ∈ S, then the LCCP in (3.3) has at least a solution. 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 (𝑃, 𝐴, 𝐵, 𝐶, 𝐷) be SSPR. Then, the LCCP in (3.3) has a unique solution for any ℎ𝛽 𝑘 , for ℎ > 0 small enough.

If 𝜆

(•) = Ψ 𝐾 ★ (•), M = Dℎ , q = ℎ𝛽 𝑘 . Using Proposition b.1 c) d), it follows that (dom(𝜑)) ∞ = 𝐾 ★ . The
Proof. This follows from [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Theorem 5.4]

(⇒ D ≻ 0), the fact that 𝐵 ⊤ 3 𝑋 3 𝐵 3 ≽ 0, Corollary a.2 (⇒ Dℎ ≻ 0 for 0 < ℎ < ℎ * Δ = 1 ∥𝐶 1 𝐵 1 ∥ 2 ∥ ( D+ D⊤ 2 ) -1 ∥ 2
⇒ S = {0}), and items 1 and 3 in Proposition 3.2. □ It is noteworthy that the SSPRness is only sufficient to guarantee that D ≻ 0, which is itself a sufficient condition for the existence of solutions to the LCCP (3.3).

boundedness and convergence analysis

From now on, we assume that the complementarity conditions are given as in (2. [START_REF] Cottle | The Linear Complementarity Problem[END_REF]Theorem 3.1.7] and [57, Lemma 1]. □

(= ℎ D + ℎ 2 𝐶 1 𝐵 1 + 𝐵 ⊤ 3 𝑋 3 𝐵 3 ) such that ∥𝜆 min 𝑘+1 ∥ ≤ 𝜂ℎ∥𝛽 𝑘 ∥. Proof. Follows from
The constant 𝜂 generally depends on ℎ since it depends on Dℎ , and we may denote it as 𝜂 (ℎ) (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 Dℎ = 0 then 𝛽 𝑘,𝑖 > 0 ⇒ 𝜆 𝑘+1,𝑖 = 0, and 𝛽 𝑘,𝑖 = 0 ⇒ 𝜆 𝑘+1,𝑖 ∈ ℝ + , hence 𝜆 min 𝑘+1 = 0 so that 𝜂 = 0. If ℎ = 0 then D0 = 𝐵 ⊤ 3 𝑋 3 𝐵 3 , 𝛽 0 = -𝐶 3 𝑥 2,𝑘 , and 𝜂 = 𝜂 0 such that ∥𝜆 min 𝑘+1 ∥ ≤ 𝜂 0 ∥𝐶 3 𝑥 2,𝑘 ∥ is independent of ℎ. If the LCP( Dℎ , ℎ𝛽 𝑘 ) enjoys the stability of solutions property [START_REF] Cottle | The Linear Complementarity Problem[END_REF]Theorem 7.3.12], then a small perturbation of Dℎ implies a small variation of the solution, hence of 𝜆 min 𝑘+1 , and consequently of the constant 𝜂. Assume that 𝐶 3 = 0 (this is Assumption 2 below). Then solutions of the LCP( Dℎ , ℎ𝛽 𝑘 ) are also solutions of LCP( D + ℎ𝐶 1 𝐵 1 , 𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 ). If D + ℎ𝐶 1 𝐵 1 is a P-matrix (which occurs if D is a P-matrix and ℎ > 0 is small enough, see Corollary a.2), then it follows from [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]Theorem B.3] [34] that the unique solution of LCP( D

+ ℎ𝐶 1 𝐵 1 , 𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 ) (which is also the least-norm solution) satisfies ∥𝜆 min 𝑘+1 ∥ ≤ 𝜂 ∥𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 ∥ for some constant 𝜂 > 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).

convergence analysis

From (3.2) and under the conditions of Corollary 4.1, the mapping 𝛽 𝑘 ↦ → 𝜆 𝑘+1 (-ℎ𝛽 𝑘 ) = 𝜆 min 𝑘+1 with 𝜆 min 𝑘+1 the least-norm element of the convex polyhedral set ( Dℎ + N ℝ 𝑚 + ) -1 (-ℎ𝛽 𝑘 ), 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)                      𝑥 1,𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 + ℎ𝐵 1 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) 𝑥 2,𝑘+1 = 𝑥 2,𝑘 + ℎ𝑧 𝑘+1 + ℎ𝐵 2 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) 0 = 𝑥 2,𝑘+1 + 𝐵 3 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) 0 ≤ 𝜆 min 𝑘+1 ⊥ 𝑤 𝑘+1 = 𝐶 1 𝑥 1,𝑘+1 + 𝐶 2 𝑥 2,𝑘+1 + 𝐶 3 𝑧 𝑘+1 + 𝐷𝜆 min 𝑘+1 ≥ 0 ∥𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) ∥ ≤ 𝜂 ∥ℎ𝛽 𝑘 ∥, ℎ𝛽 𝑘 = ℎ𝐶 1 (𝑥 1,𝑘 + ℎ𝐴 1 𝑥 1,𝑘 ) -𝐶 3 𝑥 2,𝑘 , 𝜂 > 0,
for all 𝑘 ∈ {0, . . . , 𝑁 -1}. The above results explore the well-posedness of the OSNSP, for fixed ℎ > 0.

Our goal is now to analyse the boundedness of the iterates on the integration interval [0,𝑇 ], 𝑇 > 0. For all 𝑘 ≥ 1 we have 𝑥

                             𝑥 1,𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 + ℎ𝐵 1 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) 𝑥 2,𝑘+1 = -𝑥 2,𝑘 -𝐵 3 (𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) + 𝜆 min 𝑘 (𝑥 1,𝑘 -1 , 𝑥 2,𝑘 -1 )) 0 ≤ 𝜆 min 𝑘+1 ⊥ 𝑤 𝑘+1 = 𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 -(𝐶 2 + 1 ℎ 𝐶 3 )𝑥 2,𝑘 + 1 ℎ Dℎ 𝜆 min 𝑘+1 -𝐶 2 𝐵 3 𝜆 min 𝑘 ≥ 0 0 = 𝑥 2,𝑘+1 + 𝐵 3 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ), ℎ𝑧 𝑘+1 = -𝐵 3 𝜆 min 𝑘+1 (𝑥
𝑀 𝑖+1 = 𝑀 𝑖 𝑀 + 𝑀 𝑖 -1 𝑈 , 𝑀 0 = 𝐼 2 , 𝑀 -1 = 0, 0 ≤ 𝑖 ≤ 𝑘 -1, 𝑘 ≥ 1,
where

𝑀 = 𝑎 𝑏 𝑐 𝑑 , 𝑈 = 0 0 𝑒 𝑓 , 𝑎 = ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥ (1 + ℎ 2 𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥), 𝑏 = ℎ𝜂 ∥𝐵 1 ∥ ∥𝐶 3 ∥, 𝑐 = ℎ𝜂 ∥𝐵 3 ∥ ∥𝐶 1 ∥ ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥, 𝑑 = (1 + 𝜂 ∥𝐵 3 ∥ ∥𝐶 3 ∥), 𝑒 = 𝑐, 𝑓 = 𝜂 ∥𝐵 3 ∥ ∥𝐶 3 ∥. In particular (4.4) 𝛼 𝑘+1 ≤ (𝑀 𝑘 -1 𝑀 + 𝑀 𝑘 -2 𝑈 )𝛼 1 + 𝑀 𝑘 -1 𝑈 𝛼 0 for all 𝑘 ≥ 1.
Proof. The proof starts from (4.2). The following upper bounds can be derived:

∥𝜆 min 𝑘+1 ∥ ≤ 𝜂ℎ∥𝐶 1 ∥ ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥ ∥𝑥 1,𝑘 ∥ + 𝜂 ∥𝐶 3 ∥ ∥𝑥 2,𝑘 ∥ (4.5a) ∥𝑥 1,𝑘+1 ∥ = ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥ (1 + ℎ 2 𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥) ∥𝑥 1,𝑘 ∥ + ℎ𝜂 ∥𝐵 1 ∥ ∥𝐶 3 ∥ ∥𝑥 2,𝑘 ∥ = 𝑎∥𝑥 1,𝑘 ∥ + 𝑏 ∥𝑥 2,𝑘 ∥ (4.5b) ∥𝑥 2,𝑘+1 ∥ ≤ (1 + 𝜂 ∥𝐵 3 ∥ ∥𝐶 3 ∥) ∥𝑥 2,𝑘 ∥ + 𝜂 ∥𝐵 3 ∥ ∥𝐶 3 ∥ ∥𝑥 2,𝑘 -1 ∥ + ℎ𝜂 ∥𝐵 3 ∥ ∥𝐶 1 ∥ ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥ (∥𝑥 1,𝑘 ∥ + ∥𝑥 1,𝑘 -1 ∥) = 𝑑 ∥𝑥 2,𝑘 ∥ + 𝑓 ∥𝑥 2,𝑘 -1 ∥ + 𝑐 (∥𝑥 1,𝑘 ∥ + ∥𝑥
= ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥.
Moreover 𝑑 > 1 whenever 𝐵 3 ≠ 0 and 𝐶 3 ≠ 0. The term in position (2, 2) of 𝑀 𝑛 is equal to 𝑑 𝑛 + O (ℎ 1 ℎ 2 ). Before stating the next result, let us introduce the following assumptions: Assumption 1. The matrix 𝐴 1 is such that ℎ 1 ≤ 1 for ℎ > 0 small enough.

Let us notice that the induced norm ∥ • ∥ 2,2 = 𝜎 max (•) [16, Proposition 9.4.9] is compatible with the vector norm ∥ • ∥ 2 = ∥ • ∥ [16, pp.607-609], and can be used instead of the Frobenius norm to compute ℎ 1 . The next assumption is crucial in the body of this section. Assumption 2. One has 𝐵 3,•𝑖 ∈ Ker(𝑋 3 ) for all 1 ≤ 𝑖 ≤ 𝑚 (=⇒ 𝐶 3 = 0).

Let us remind that 𝑋 3 is uniquely defined under minimality and positive-realness [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF]Theorem 5.4]. If Assumption 2 holds then 𝑑 = 1 and 𝑏 = 𝑓 = 0. This assumption means that the variable 𝑧 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 𝛽 𝑘 no longer depends on 1 ℎ 𝐶 3 𝑥 2,𝑘 . The terms in (4.5) involve products with 𝜂 (ℎ), and it is important to guarantee that they are bounded, when Assumption 2 holds.

Assumption 3. The constant 𝜂 = 𝜂 (ℎ) is such that the products 𝜂 (ℎ)ℎ can be made arbitrarily small for small enough ℎ ≥ 0.

Notice that this assumption is verified in the examples of Section 6, where 𝜂 is a constant (but these examples do not verify Assumption 2, however). The matrices 𝑀 𝑘 then involve only products of powers of the terms 𝑎, 𝑏, 𝑐, 𝑑 = 1, 𝑒. These products are bounded for bounded 𝑘, and provided that both ℎ and ℎ 1 are small enough so that 𝑎 ≤ 1, 𝑏 ≤ 1, 𝑐 ≤ 1, 𝑑 = 1, 𝑒 ≤ 1, they remain bounded as 𝑘 → +∞. If the inequalities are strict (apart form 𝑑 = 1), the first rows in both 𝑀 𝑘 -1 𝑀 + 𝑀 𝑘 -2 𝑈 and 𝑀 𝑘 -1 𝑈 converge to zero, while the second rows converge to (0, 1) (the 1 being the consequence of a term 𝑑 𝑛 in matrices 𝑀 𝑛 ). 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 ℎ and ℎ 1 are small enough, so that 𝑎 ≤ 1, 𝑏 ≤ 1, 𝑐 ≤ 1, 𝑑 = 1, 𝑒 ≤ 1. Then the sequences {𝛼 1,𝑘 } 𝑘 and {𝛼 2,𝑘 } 𝑘 are uniformly bounded in 𝑘. If 𝐴 1 ≺ 0, and ℎ > 0 is small enough such that 𝑎 < 1, 𝑐 < 1, 𝑒 < 1, then 𝛼 1,𝑘 → 0 as 𝑘 → +∞.

Recall that 𝛼 1,𝑘 = ∥𝑥 1,𝑘 ∥ and 𝛼 2,𝑘 = ∥𝑥 2,𝑘 ∥. 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 ∥𝜆 min 𝑘+1 ∥ is uniformly bounded in 𝑘 and ℎ > 0, and it converges to zero as 𝑘 → +∞ if 𝛼 1,𝑘 → 0 as 𝑘 → +∞. Therefore the sequence

𝑥 1,𝑘+1 -𝑥 1,𝑘 ℎ
𝑘 is uniformly bounded also in 𝑘 and ℎ > 0.

Proof. The boundedness of ∥𝜆 min 𝑘+1 ∥ and its convergence follow from (4.5) (a), the boundedness of ∥𝑥 1,𝑘 ∥ and Assumption 2. From the first line in (4.2), one has

𝑥 1,𝑘+1 -𝑥 1,𝑘 ℎ = 𝐴 1 𝑥 1,𝑘 + 𝐵 1 𝜆 min 𝑘+1
, hence the result is proved.

□

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

(4.7)              𝑥 ℎ 1 (𝑡) = 𝑥 1,𝑘+1 + 𝑡 𝑘+1 -𝑡 ℎ (𝑥 1,𝑘 -𝑥 1,𝑘+1 ) 𝑥 ℎ 1 (𝑡) = 𝑥 1,𝑘+1 -𝑥 1,𝑘 ℎ 𝑥 ★ 1,ℎ (𝑡) = 𝑥 1,𝑘+1 for all 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 )
and similarly for the other variables. Corollary 4.6. There exists subsequences denoted as {𝑥 ℎ 1 (•)} ℎ≥0 which converge uniformly to a limit

𝑥 1 (•) in C 0 ([0,𝑇 ]; ℝ 𝑛 1 ), such that almost everywhere one has 𝑥 1 (𝑡) = 𝐴 1 𝑥 1 (𝑡) + 𝜆 1 (𝑡), for some function 𝜆 1 (•). The sequence {𝑥 ℎ★ 1 (•)} ℎ≥0 converges to 𝑥 1 (•) strongly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ). Moreover the sequences {𝑥 ★ 2,ℎ (•)} ℎ≥0 and {𝜆 min,★ ℎ (•)} ℎ≥0 converge weakly★ in L ∞ ( [0,𝑇 ]; ℝ 𝑙 ), 𝑙 = 𝑛 2 and 𝑙 = 𝑚, respectively, to limits 𝑥 2 (•) and 𝜆(•) which satisfy the continuous-time constraints. Proof. Both {𝑥 ℎ 1 (•)} ℎ and { 𝑥 ℎ 1 (•)} ℎ are uniformly bounded on [0,𝑇 ]. The sequence of C 0 ( [0,𝑇 ]; ℝ 𝑛 1 ) functions {𝑥 ℎ 1 (•)} ℎ is thus equicontinuous.
By the Arzela-Ascoli Theorem this sequence stays in a compact subset of C 0 ( [0,𝑇 ]; ℝ 𝑛 1 ) and converges uniformly towards a continuous limit 𝑥 1 (•) : [0,𝑇 ] → ℝ 𝑛 1 as ℎ → 0, ℎ > 0. Let us notice that 𝑥 1 (𝑡) → 0 as 𝑡 → +∞. Indeed from the above results which are stated for 𝑘 ∈ {0, . . . , 𝑁 } and ℎ > 0, 𝑘 → +∞ only if 𝑇 = 𝑁ℎ → +∞. Let us now define the step functions 𝑥 ℎ★

1 (𝑡) = 𝑥 1,𝑘+1 if 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ), 𝜆 min ℎ★ (𝑡) = 𝜆 min 𝑘+1 if 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ), and
(4.8) 𝑥 ℎ 1 (𝑡) = 𝐴 1 𝑥 ℎ★ 1 (𝑡) + 𝐵 1 𝜆 min ℎ★ (𝑡), almost everywhere.
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

In addition one has:

(4.9) ∥𝑥 ℎ 1 (𝑡) -𝑥 ℎ★ 1 (𝑡) ∥ 2 L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) = 𝑁 -1 ∑ 𝑘=0 ∫ 𝑡 𝑘+1 𝑡 𝑘 (𝑡 𝑘+1 -𝑡) 2 ∥ 𝑥 1,𝑘 ∥ 2 𝑑𝑡 ≤ 𝐶 2 𝑁 -1 ∑ 𝑘=0 ∫ 𝑡 𝑘+1 𝑡 𝑘 (𝑡 𝑘+1 -𝑡) 2 𝑑𝑡 ≤ 𝐶 2 𝑁 -1 ∑ 𝑘=0 1 3 (𝑡 𝑘+1 -𝑡 𝑘 ) 3 = 𝐶 2 𝑁ℎ 3 3 = 𝐶 2 𝑇ℎ 2 3
,

where 𝑥 1,𝑘 = 𝑥 1,𝑘+1 -𝑥 1,𝑘 ℎ
, and ∥ 𝑥 1,𝑘 ∥ ≤ 𝐶 for all 𝑘 ≥ 0 for some 𝐶 > 0 from Corollary 4.5. It follows that

𝑥 ℎ★ 1 (•) → 𝑥 1 (•) strongly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) as ℎ → 0. On the other hand since 𝑥 ℎ 1 (•) ∈ L ∞ ( [0,𝑇 ]; ℝ 𝑛 1 ) it follows from the Banach-Alaoglu-Bourbaki Theorem [18, Theorem III.15] that 𝑥 ℎ 1 (•) → 𝑥 1,lim (•) in the weak★ topology as ℎ → 0. Since the considered domain [0,𝑇 ] is bounded, one has L 1 ( [0,𝑇 ]; ℝ 𝑛 1 ) ⊆ L 2 ([0,𝑇 ]; ℝ 𝑛 1 ), thus 𝑥 ℎ 1 (•) converges weakly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ). Moreover since 𝑥 ℎ 1 (𝑡) = 𝑥 ℎ 1 (0) + ∫ 𝑡 0 𝑥 ℎ 1 (𝑡)𝑑𝑡, one has 𝑥 1,lim (•) = 𝑥 1 (•)
almost everywhere. Therefore we deduce using (4.8) that (𝐵 1 𝜆 min ℎ★ ) (•) converges weakly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) as ℎ → 0 towards some limit 𝜆 1 (•), and we have 𝑥

1 (𝑡) = 𝐴 1 𝑥 1 (𝑡) + 𝜆 1 (𝑡) almost everywhere on [0,𝑇 ].
Let us consider now the step-functions

𝑥 ★ 2,ℎ (•) and 𝜆 min,★ ℎ (•) defined as 𝑥 ★ 2,ℎ (𝑡) = 𝑥 2,𝑘+1
for all 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ), and similarly for the multiplier. It follows from Proposition 4.4 and from Corollary 4.5 that both sequences {𝑥 ★ 2,ℎ (•)} ℎ≥0 and {𝜆 min,★ ℎ (•)} ℎ≥0 are uniformly bounded. Hence from the Banach-Alaoglu-Bourbaki Theorem, there exist subsequences of both which converge weakly★ to limits

𝑥 2 ∈ L ∞ ([0,𝑇 ]; ℝ 𝑛 2 ) and 𝜆 ∈ L ∞ ( [0,𝑇 ]; ℝ 𝑚 , respectively. We infer that 0 = 𝑥 ★ 2,ℎ (𝑡) + 𝐵 3 𝜆 ★ ℎ (𝑡) → 𝑥 2 (𝑡) + 𝐵 3 𝜆(𝑡)
for almost all 𝑡 ≥ 0 in the weak★ sense. Next, using the second line in (4.1) it follows that the sequence {ℎ𝑧 𝑘 } ℎ≥0 is also uniformly bounded. Denoting 𝜎 𝑘 = ℎ𝑧 𝑘 , it follows that the sequence {𝜎 ★ ℎ } ℎ≥0 converges weakly★ to a limit denoted as

𝜎 (𝑡) ∈ L ∞ ( [0,𝑇 ]; ℝ 𝑛 2 ). Moreover one has 𝜎 (𝑡) = 0 almost everywhere since 𝑥 ★ 2,ℎ (𝑡) -𝑥 ★ 2,ℎ (𝑡 -𝑡 𝑘 ) -ℎ𝐵 2 𝜆 min,★ ℎ (𝑡) = 𝜎 ★ ℎ (𝑡) almost everywhere in [0,𝑇 ].
Finally, the complementarity conditions in (4.1) are equivalently rewritten as 0

≤ 𝜆 ★ ℎ (𝑡) ⊥ 𝐶 1 𝑥 ★ 1,ℎ (𝑡) + 𝐶 2 𝑥 ★ 2,ℎ (𝑡) + 𝐷𝜆 ★ ℎ (𝑡) ≥ 0 for all 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ), since 𝐶 3 = 0 by Assumption 2.
It is deduced that in the limit as ℎ → 0 one obtains 0 ≤ 𝜆(𝑡) ⊥ 𝐶 1 𝑥 1 (𝑡) + 𝐶 2 𝑥 2 (𝑡) + 𝐷𝜆(𝑡) ≥ 0. □

stability analysis via a lyapunov function

To complete the foregoing section, it is natural to analyze how the Lyapunov function 𝑉 (𝑥 1 , 𝑥 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 [START_REF] Greenhalgh | On preserving dissipativity properties of linear complementarity dynamical systems with the 𝜃 -method[END_REF]. It is noteworthy that Assumption 2 implies that 𝑋 3 𝑥 2,𝑘 = -𝑋 3 𝐵 3 𝜆 𝑘 = 0 and 𝑋 3 𝑥 2 (𝑡) = -𝑋 3 𝐵 3 𝜆(𝑡) = 0 for all 𝑘 and 𝑡. Therefore 𝑉 (𝑥 1 , 𝑥 2 ) = 1 2 𝑥 ⊤ 1 𝑋 1 𝑥 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

𝑉 𝑘 Δ = 1 2 𝑥 ⊤ 1,𝑘 𝑋 1 𝑥 1,𝑘 +𝑥 ⊤ 2,𝑘 𝑋 3 𝑥 2,𝑘 . It follows that: (4.10) 𝑉 𝑘+1 -𝑉 𝑘 = 1 2 (𝑥 1,𝑘+1 -𝑥 1,𝑘 ) ⊤ 𝑋 1 (𝑥 1,𝑘+1 + 𝑥 1,𝑘 ) + (𝑥 2,𝑘+1 -𝑥 2,𝑘 ) ⊤ 𝑋 3 (𝑥 2,𝑘+1 + 𝑥 2,𝑘 ) = 1 2 (ℎ𝐴 1 𝑥 1,𝑘 + ℎ𝐵 1 𝜆 𝑘+1 ) ⊤ 𝑋 1 (2𝑥 1,𝑘 + ℎ𝐴 1 𝑥 1,𝑘 + ℎ𝐵 1 𝜆 𝑘+1 ) + (ℎ𝑧 𝑘+1 + ℎ𝐵 2 𝜆 𝑘+1 ) ⊤ 𝑋 3 (2𝑥 2,𝑘 + ℎ𝑧 𝑘+1 + ℎ𝐵 2 𝜆 𝑘+1 ) = ℎ(𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝑥 1,𝑘 + 𝜆 ⊤ 𝑘+1 𝐵 ⊤ 1 𝑋 1 𝑥 1,𝑘 ) + ℎ 2 ( 1 2 𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝐴 1 𝑥 1,𝑘 + 1 2 𝜆 ⊤ 𝑘+1 𝐵 ⊤ 1 𝑋 1 𝐵 1 𝜆 𝑘+1 ) + ℎ 2 𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝐵 1 𝜆 𝑘+1 + ℎ𝑧 ⊤ 𝑘+1 𝐶 3 𝜆 𝑘+1 + ℎ𝜆 ⊤ 𝑘+1 𝐶 3 𝐵 2 𝜆 𝑘+1 + ℎ𝑧 ⊤ 𝑘+1 𝑋 3 𝑥 2,𝑘 + ℎ𝜆 ⊤ 𝑘+1 𝐵 ⊤ 2 𝑋 3 𝑥 2,𝑘 + ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 -ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

= ℎ 2 𝑥 ⊤ 1,𝑘 (𝐴 ⊤ 1 𝑋 1 + 𝑋 1 𝐴 1 + ℎ𝐴 ⊤ 1 𝑋 1 𝐴 1 )𝑥 1,𝑘 -ℎ𝜆 ⊤ 𝑘+1 ( D + ℎ(𝐶 1 𝐵 1 -𝐵 ⊤ 2 𝑋 3 𝐵 2 - 1 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 ))𝜆 𝑘+1 + ℎ𝑥 ⊤ 1,𝑘 (𝐵 ⊤ 1 𝑋 1 -𝐶 1 -ℎ𝐴 ⊤ 1 𝐶 ⊤ 1 + ℎ𝐴 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 + ℎ 2 𝑥 ⊤ 2,𝑘 𝑋 3 𝐵 2 𝜆 𝑘+1 + ℎ 2 𝑧 ⊤ 𝑘+1 𝑋 3 𝑥 2,𝑘 + ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 = ℎ 2 (𝑥 ⊤ 1,𝑘 , 𝜆 ⊤ 𝑘+1 ) Δ =𝑄 1 𝐴 ⊤ 1 𝑋 1 + 𝑋 1 𝐴 1 (𝐵 ⊤ 1 𝑋 1 -𝐶 1 ) ⊤ 𝐵 ⊤ 1 𝑋 1 -𝐶 1 -( D + D⊤ ) 𝑥 1,𝑘 𝜆 𝑘+1 + ℎ 2 2 𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝐴 1 𝑥 1,𝑘 -ℎ 2 𝜆 ⊤ 𝑘+1 (𝐶 1 𝐵 1 -𝐵 ⊤ 2 𝑋 3 𝐵 2 - 1 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 + ℎ 2 𝑥 ⊤ 1,𝑘 (-𝐴 ⊤ 1 𝐶 ⊤ 1 + 𝐴 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 + ℎ 2 𝑥 ⊤ 2,𝑘 𝑋 3 𝐵 2 𝜆 𝑘+1 + ℎ 2 𝑧 ⊤ 𝑘+1 𝑋 3 𝑥 2,𝑘 + ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 = ℎ 2 (𝑥 ⊤ 1,𝑘 , 𝜆 ⊤ 𝑘+1 )𝑄 1 𝑥 1,𝑘 𝜆 𝑘+1 + ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 + ℎ 2 1 2 𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝐴 1 𝑥 1,𝑘 -𝜆 ⊤ 𝑘+1 (𝐶 1 𝐵 1 -𝐵 ⊤ 2 𝑋 3 𝐵 2 - 1 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 +𝑥 ⊤ 1,𝑘 (-𝐴 ⊤ 1 𝐶 ⊤ 1 + 𝐴 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 + ℎ 2 𝑥 ⊤ 2,𝑘 𝑋 3 𝐵 2 𝜆 𝑘+1 + ℎ 2 𝑧 ⊤ 𝑘+1 𝑋 3 𝑥 2,𝑘 = ℎ 2 (𝑥 ⊤ 1,𝑘 , 𝜆 ⊤ 𝑘+1 )𝑄 1 𝑥 1,𝑘 𝜆 𝑘+1 + ℎ𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 + ℎ 2 1 2 𝑥 ⊤ 1,𝑘 𝐴 ⊤ 1 𝑋 1 𝐴 1 𝑥 1,𝑘 -𝜆 ⊤ 𝑘+1 (𝐶 1 𝐵 1 -𝐵 ⊤ 2 𝑋 3 𝐵 2 - 1 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 +𝑥 ⊤ 1,𝑘 (-𝐴 ⊤ 1 𝐶 ⊤ 1 + 𝐴 ⊤ 1 𝑋 1 𝐵 1 )𝜆 𝑘+1 - 1 2 𝑥 ⊤ 2,𝑘 𝑋 3 𝑥 2,𝑘 + 𝜆 ⊤ 𝑘+1 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 𝑘 =0 if Assumption 2 holds.
, where the last two terms are obtained using the second equality in (3.1). Remind that 𝜆 ⊤ 𝑘+1 𝑤 𝑘+1 = 0. Assume that 𝑄 1 ≺ 0, then let Assumption 2 hold. Let us rewrite (4.10) as 𝑉 𝑘+1 -

𝑉 𝑘 = ℎ 2 (𝑥 ⊤ 1,𝑘 , 𝜆 ⊤ 𝑘+1 ) [𝑄 1 + ℎ𝑄 2 ] 𝑥 1,𝑘 𝜆 𝑘+1
, where the indefinite matrix

(4.11) 𝑄 2 = 𝐴 ⊤ 1 𝑋 1 𝐴 1 -1 2 𝐴 ⊤ 1 𝐶 ⊤ 1 + 1 2 𝐴 ⊤ 1 𝑋 1 𝐵 1 (-1 2 𝐴 ⊤ 1 𝐶 ⊤ 1 + 1 2 𝐴 ⊤ 1 𝑋 1 𝐵 1 ) ⊤ -𝐶 1 𝐵 1 + 𝐵 ⊤ 2 𝑋 3 𝐵 2 -1 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 ,
and we remind that Assumption 2 implies that

𝑉 𝑘 = 1 2 𝑥 ⊤ 1,𝑘 𝑋 1 𝑥 1,𝑘 .
The matrix 𝑄 2 is treated as a perturbation to 𝑄 1 . Using Corollary a.2, there exists ℎ max > 0 such that for all ℎ ∈ [0, ℎ max ) one has

𝑄 3 Δ = 𝑄 1 + ℎ𝑄 2 ≺ 0.

Using classical bounds on quadratic forms and vector/matrices products (𝑥

⊤ 𝑦 ≤ ∥𝑥 ∥ ∥𝑦 ∥ ≤ 1 2 ∥𝑥 ∥ 2 + 1 2 ∥𝑦 ∥ 2 [16, Corollary 9.1.7], if 𝑄 = 𝑄 ⊤ then 𝜆 min (𝑄) ∥𝑥 ∥ 2 ≤ 𝑥 ⊤ 𝑄𝑥 ≤ 𝜆 max (𝑄) ∥𝑥 ∥ 2 [16, Corollary 8.4.2]
, ∥𝑄𝑥 ∥ ≤ ∥𝑄 ∥ ∥𝑥 ∥ for compatible norms [16, section 9.3]), and noting that we can replace 𝑄 3 by its negative definite symmetric part in the quadratic form, it follows that: (4.12)

𝑉 𝑘+1 ≤ 𝑉 𝑘 + ℎ 2 𝜆 max (𝑄 3 ) ∥ (𝑥 1,𝑘 , 𝜆 𝑘+1 ) ∥ 2 ≤ 𝑉 𝑘 + ℎ 2 𝜆 max (𝑄 3 ) ∥𝑥 1,𝑘 ∥ 2 ≤ 𝑉 𝑘 + ℎ 2 𝜆 max (𝑄 3 ) 𝜆 max (𝑋 1 ) 𝑉 𝑘 = 1 + ℎ 2 𝜆 max (𝑄 3 ) 𝜆 max (𝑋 1 )
𝑉 𝑘 , for all 𝑘 ≥ 0.

Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

Since 𝜆 max (𝑄 3 ) < 0 (which is the reason why the second inequality in (4.12) holds), one has 1

+ ℎ 2 𝜆 max (𝑄 3 )
𝜆 max (𝑋 1 ) < 1 for all ℎ > 0. Therefore all 𝑉 𝑘 are uniformly bounded for any bounded 𝑉 0 , and 𝑉 𝑘 → 0 as 𝑘 → +∞ for any ℎ ∈ (0, ℎ max ). Let us set 𝜆 𝑘+1 = 𝜆 min 𝑘+1 . Since the dependence of 𝜆 min 𝑘+1 in 𝑥 2,𝑘 through 𝛽 𝑘 is 𝐶 3 𝑥 2,𝑘 (see (4.2)), it follows from Assumptions 2 and 3 that 𝜆 min 𝑘+1 = 𝜆 min 𝑘+1 (𝑥 1,𝑘 ), and 𝜆 min 𝑘+1 (𝑥 1,𝑘 ) → 0 as 𝑘 → +∞. Using the third line in (4.2), we infer that the iterates 𝑥 2,𝑘 (which are 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 (𝑃, 𝐴, 𝐵, 𝐶, 𝐷) be strongly SPR and minimal. Then there exists ℎ max > 0 depending on 𝐴 1 , 𝐶 1 , 𝑋 1 , 𝑋 3 , 𝐵 1 , 𝐵 2 , such that for all ℎ ∈ (0, ℎ max ), and for all bounded initial condition 𝑥 1,0 , the sequences {𝑥 1,𝑘 } 𝑘 , {𝑥 2,𝑘 } 𝑘 , {𝜆 min 𝑘 } 𝑘 are uniformly bounded and 𝑥 1,𝑘 , 𝑥 2,𝑘 , 𝜆 min 𝑘+1 → 0 as 𝑘 → +∞. Remark 4.8. Concerning Assumption 1, let us examine the constants 𝑎, 𝑏, 𝑐, 𝑑, 𝑓 in Lemma 4.3, and the inequalities in (4.5), as well as the above Lyapunov analysis. The term ℎ 1 is present in both 𝑎 and 𝑐. The inequality 𝑎 ≤ 1 implies ℎ 1 < 1, but 𝑐 ≤ 1 does not imply that ℎ 1 < 1. The constant 𝑎 is present in (4.5) (b), however the Lyapunov analysis guarantees the boundedness of 𝑥 1,𝑘 without this assumption (the analysis relies on the passivity of (𝐴 1 , 𝐵 1 , 𝐶 1 , D) and ℎ ≥ 0 small enough). But, Lemma 4.7 relies on 𝑄 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. (4.11) and 𝑄 1 can be rewritten as:

Let (𝐴 1 , 𝐵 1 , 𝐶 1 , D) be strictly state passive (⇒ 𝑋 1 𝐴 1 + 𝐴 ⊤ 1 𝑋 1 = -𝜇𝑋 1 , 𝜇 > 0, 𝑋 1 = 𝑋 ⊤ 1 ≻ 0, 𝜇 is called the dissipation constant) and D + D⊤ = 0 ⇒ 𝑋 1 𝐵 1 = 𝐶 ⊤ 1 , then 𝑄 2 in
(4.13) 𝑄 2 = 𝐴 ⊤ 1 𝑋 1 𝐴 1 0 0 𝐵 ⊤ 2 𝑋 3 𝐵 2 -3 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 , 𝑄 1 = -𝜇𝑋 1 0 0 0 Thus 𝑄 1 + ℎ𝑄 2 ≺ 0 if and only if 𝜇𝑋 1 -ℎ𝐴 ⊤ 1 𝑋 1 𝐴 1 ≻ 0 and 𝐵 ⊤ 2 𝑋 3 𝐵 2 -3 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 ≺ 0. Let us recall that if (𝐴 1 , 𝐵 1 , 𝐶
constant 𝜇 > 0, 2) 𝐵 ⊤ 2 𝑋 3 𝐵 2 -3 2 𝐵 ⊤ 1 𝑋 1 𝐵 1 ≺ 0.
Let Assumptions 2 and 3 hold. Then there exists ℎ max > 0 such that for all ℎ ∈ [0, ℎ max ), the sequences {𝑥 1,𝑘 } 𝑘 , {𝑥 2,𝑘 } 𝑘 and {𝜆 min 𝑘 } 𝑘 are uniformly bounded.

Proof. Denote X1 = 𝜇𝑋 1 -ℎ𝐴 ⊤ 1 𝑋 1 𝐴 1 , then X1 ≻ 0 for ℎ ≥ 0 small enough, using Corollary a.2. Using 2) and Assumption 2, and after few calculations, there exists ℎ max > 0 such that:

(4.14) 𝑉 𝑘+1 ≤ 1 - ℎ 2 𝜆 min ( X1 ) 𝜆 max (𝑋 1 ) 𝑉 𝑘 .
for all ℎ ∈ [0, ℎ max ). Thus 𝑉 𝑘 → 0 as 𝑘 → +∞, and 𝑥 1,𝑘 → 0 as 𝑘 → +∞. It can be inferred that

∥𝑥 1,𝑘 ∥ ≤ √ 𝜆 max (𝑋 1 )
𝜆 min (𝑋 1 ) ∥𝑥 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 𝜇) is related to the value of ℎ max : the greater 𝜇, the greater ℎ max . Finally it is noteworthy that the strict state passivity in item 1) of Lemma 4.9 can again be relaxed, provided that the matrix

𝐴 ⊤ 1 𝑋 1 + 𝑋 1 𝐴 1 + ℎ𝐴 ⊤ 1 𝑋 1 𝐴 1 ≺ 0. Consider as an example 𝐴 1 = -1 0 0 0 , 𝑋 1 = 𝐼 2 , then 𝐴 ⊤ 1 𝑋 1 + 𝑋 1 𝐴 1 + ℎ𝐴 ⊤ 1 𝑋 1 𝐴 1 = -2 + ℎ 0 0 0 ≺ 0 for all 0 ≤ ℎ < 2.
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

another approach for passive systems

The LDCS in (2.2) are suitably transformed in [START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF] 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)

             (𝑎) 𝐼 𝑛 1 0 0 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝑥 1 (𝑡) 𝜆(𝑡) = 𝐴 1 𝐵 1 -𝐶 1 -D 𝑥 1 (𝑡) 𝜆(𝑡) + 0 𝐼 𝑚 𝑤 (𝑡) (𝑏) 0 ≤ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 𝐶 1 𝑥 1 (𝑡) + 𝐶 2 𝑥 2 (𝑡) + 𝐶 3 𝑧 (𝑡) + 𝐷𝜆(𝑡) ≥ 0 (𝑐) 𝑥 2 (𝑡) = -𝐵 3 𝜆(𝑡), 𝑧 (𝑡) = -𝐵 2 𝜆(𝑡) -𝐵 3 𝜆(𝑡).
It is noteworthy that this dynamics is obtained from the original one by treating the multiplier 𝜆 as a state variable. The "equivalence" between both systems is to be understood first as follows [32, Lemma 3]:

{(𝑥 1 , 𝜆, 𝑤) | (𝑥 1 , 𝑥 2 , 𝜆, 𝑤) is a solution of (2.2) (𝑎, 𝑏, 𝑐)} = {(𝑥 1 , 𝜆, 𝑤) | (𝑥 1 , 𝜆, 𝑤
) is a solution of (5.1) (𝑎) and 𝐵 3 𝜆 is absolutely continuous} .

By solution of (2.2) it is meant that (𝑥 1 , 𝑥 2 ) is absolutely continuous and the dynamics holds Lebesgue almost everywhere [START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF]Definition 3]. The initial conditions are chosen compatible with the constraints. It is noteworthy that the variables 𝑧 and 𝜆 are not shown to be absolutely continuous, and examples in Sections 6.1 and 6.3 will illustrate this fact. One has 𝑤 (𝑡) ∈ -N ℝ 𝑚 + (𝜆(𝑡)), therefore one obtains the differential inclusion (DI):

(5.2)

𝐼 𝑛 1 0 0 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝑥 1 (𝑡) 𝜆(𝑡) - 𝐴 1 𝐵 1 -𝐶 1 -D 𝑥 1 (𝑡) 𝜆(𝑡) ∈ - 0 𝐼 𝑚 N ℝ 𝑚 + (𝜆(𝑡)) = - 0 𝐼 𝑚 N ℝ 𝑚 + (0 𝐼 𝑚 ) 𝑥 1 (𝑡) 𝜆(𝑡) = -N Φ 𝑥 1 (𝑡) 𝜆(𝑡) with Φ Δ = {(𝑥 1 , 𝜆) ∈ ℝ 𝑛 1 × ℝ 𝑚 | 𝜆 ≥ 0}. If this DI has a solution (𝑥 1 (•), 𝜆(•)), then 𝑥 2 (•), 𝑧 (•)
, and 𝑤 (•) can be obtained (provided indeed that 𝐵 3 𝜆(•) is absolutely continuous). Multiplying the first line in (5.2) by 𝑋 1 and adding an exogenous excitation 𝐺 (𝑡) to the left-hand side, one obtains (5.3)

𝑋 1 0 0 𝐵 ⊤ 3 𝑋 3 𝐵 3 Δ =P 𝑥 1 (𝑡) 𝜆(𝑡) - 𝑋 1 𝐴 1 𝑋 1 𝐵 1 -𝐶 1 -D Δ =A 𝑥 1 (𝑡) 𝜆(𝑡) + 𝐺 (𝑡) ∈ -N Φ 𝑥 1 (𝑡) 𝜆(𝑡)
which is similar to the form obtained in [32, Theorem 3]. Let us assume that an exogenous input

𝑉 (𝑡) = 𝑣 (𝑎) (𝑡), 𝑣 (𝑏 ) (𝑡), 𝑣 (𝑐 ) (𝑡), 𝑣 (𝑑 ) (𝑡) ⊤ is added to (2.
2), then:

(5.4)

𝐺 (𝑡) = 𝑣 (𝑎) (𝑡) 𝐶 3 𝑣 (𝑏 ) (𝑡) + 𝑣 (𝑐 ) (𝑡) + 𝐶 2 𝑣 (𝑐 ) (𝑡) -𝑣 (𝑑 ) (𝑡)
.

Notice that 𝐵 ⊤ 3 𝑋 3 𝐵 3 ≽ 0, therefore P = P ⊤ ≽ 0. However P has a specific block-diagonal structure. Under strong SPRness and minimality, A ≺ 0, since A + A ⊤ ≺ 0 from the KYP Lemma [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF] (then 𝐴 1 is a Hurwitz matrix, 𝑋 1 = 𝑋 ⊤ 1 ≻ 0 and D ≻ 0). The exogeneous term 𝐺 (𝑡) stems from the exogeneous terms 𝐸 (𝑡) and 𝐹 (𝑡) 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]:

{(𝑥 1 , 𝜆) | (𝑥 1 , 𝑥 2 , 𝜆, 𝑤) is a solution of (2.2) (a,b,c,d)} = {(𝑥 1 , 𝜆) | (𝑥 1 , 𝜆
) is a solution of (5.3) and 𝐵 3 𝜆 (= -𝑥 2 ) is absolutely continuous}.

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).

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 𝜉 (𝑡) -A𝜉 (𝑡) + 𝐺 (𝑡) ∈ -N Φ (𝜉 (𝑡)). The Euler implicit scheme associated with (5.3) is:

(5.5)

P𝜉 𝑘+1 -(P + ℎA) Δ =P ℎ 𝜉 𝑘 + ℎ𝐺 𝑘 ∈ -N Φ (𝜉 𝑘+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.

(5.6) S Δ = {(𝑥 1 , 𝜆) ∈ ℝ 𝑛 1 × ℝ 𝑚 | 𝑥 1 ∈ Ker(𝑋 1 ), 𝜆 ∈ Ker(𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∩ (𝐵 ⊤ 3 𝑋 3 𝐵 3 ) -1 (𝐾) ∩ 𝐾 ★ }.
Then:

1. If S = {0}, for each 𝑞 𝑘 ∈ ℝ 𝑛 1 ×𝑚 the problem VI(P, 𝑞 𝑘 , Ψ Φ ) has at least one solution.

2. Assume that S ≠ {0}. If there exists 𝜉 0 ∈ dom(Φ) such that ⟨𝑞 𝑘 -P ⊤ 𝜉 0 , 𝑣⟩ > 0 for all 𝑣 ∈ S, 𝑣 ≠ 0, then the problem VI(P, 𝑞 𝑘 , Ψ Φ ) has at least one solution. In case

𝑋 1 = 𝑋 ⊤ 1 , one can take 𝜉 0 = 0. 3. If 𝜉 1 𝑘+1 and 𝜉 2 𝑘+1 are two solutions of VI(P, 𝑞 𝑘 ,𝜓 Φ ), then 𝜉 1 𝑘+1 -𝜉 2 𝑘+1 ∈ Ker(P), hence 𝑥 1,𝑘+1 is unique. Moreover ⟨𝑞 𝑘 , 0 𝜆 1 𝑘+1 -𝜆 2 𝑘+1 ⟩ = 0.
4. If 𝑋 1 ≻ 0 and 𝐵 ⊤ 3 𝑋 3 𝐵 3 ≻ 0, then the problem VI(P, 𝑞 𝑘 , Ψ Φ ) always has a unique solution. 5. Assume that the quintuple (𝑃, 𝐴, 𝐵, 𝐶, 𝐷) is minimal and strongly SPR, so that the quadruple (𝐴 1 , 𝐵 1 , 𝐶 1 , D) is minimal and strongly SPR. If 𝐵 ⊤ 3 𝑋 3 𝐵 3 ≻ 0, then the problem VI(P, 𝑞 𝑘 , Ψ Φ ) always has a unique solution.

Let 𝐾

= 𝐾 ★ = ℝ 𝑚 + . Then S = {(𝑥 1 , 𝜆) | 𝑥 1 ∈ Ker(𝑋 1 ), 𝜆 ∈ ℝ 𝑚 + ∩ Ker(𝐵 ⊤ 3 𝑋 3 𝐵 3 )}. Proof. From Proposition b.1 d) one has (𝐾 ★ ) ∞ = 𝐾 ★ and Φ ∞ = Φ. Also 𝜑 = Ψ Φ hence from Proposi- tion b.1 c) and d), 𝜑 ∞ = Ψ ℝ 𝑛 1 ×𝐾 ★ and (dom(𝜑 ∞ )) ★ = Φ ★ = ℝ 𝑛 1 ×𝐾. Then one has S Δ = Φ∩Ker(P) ∩ {𝜉 ∈ ℝ 𝑛 1 +𝑚 | P 𝑥 1 𝜆 = 0 𝛼 , 𝛼 ∈ 𝐾 } = Φ ∩ Ker(P) ∩ {𝜉 ∈ ℝ 𝑛 1 +𝑚 | 𝑥 1 ∈ Ker(𝑋 1 ), 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 ∈ 𝐾 },
where we used the block diagonal structure of P.

Therefore S = {(𝑥 1 , 𝜆) ∈ ℝ 𝑛 1 × ℝ 𝑚 | 𝑥 1 ∈ Ker(𝑋 1 ), 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 ∈ 𝐾, 𝜆 ∈ Ker(𝐵 ⊤ 3 𝑋 3 𝐵 3 ), 𝜆 ∈ 𝐾 ★ }.
In item 2) one uses that Ψ Φ (𝑣) = 0 when 𝑣 ∈ S, and item b') in Proposition b.2. Let us now deal with item 4. In this case P = P ⊤ ≻ 0 and the result follows from items 1 and 3. Let us define R = R ⊤ ≻ 0, R 2 = P, and 𝜁 𝑘 = R𝜉 𝑘 . Then the VI in (5.5) is equivalently rewritten as:

(5.7)

𝜁 𝑘+1 -R -1 P ℎ R -1 𝜁 𝑘 + R -1 𝐺 𝑘 ∈ -R -1 N Φ (R -1 𝜁 𝑘+1 ) = -N Φ (𝜁 𝑘+1 ),
where the last equality stems from the chain rule, with

Φ = {𝑧 ∈ ℝ 𝑛 1 +𝑚 | R -1 𝑧 ∈ Φ}. One infers from (5.7) that 𝜁 𝑘+1 = proj[ Φ; R -1 P ℎ R -1 𝜁 𝑘 -R -1 𝐺 𝑘 ] ⇒ 𝜉 𝑘+1 = R -1 proj[ Φ; R -1 P ℎ 𝜉 𝑘 -R -1 𝐺 𝑘 ].
The first assertion in item 5 about strong SPRness follows from [28, Theorem 5.4 (A')], which guarantees that 𝑋 1 ≻ 0 (see Appendix d). The last assertion of item 5 is a consequence of item 4. Consider finally item 6. One has

(𝐵 ⊤ 3 𝑋 3 𝐵 3 ) -1 (𝐾) = {𝜆 ∈ ℝ 𝑚 | 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 = 𝛼, 𝛼 𝑖 ≥ 0, 1 ≤ 𝑖 ≤ 𝑚}, see e.g., [53, Example 15]. Hence Ker(𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∩ (𝐵 ⊤ 3 𝑋 3 𝐵 3 ) -1 (𝐾) = Ker(𝐵 ⊤ 3 𝑋 3 𝐵 3 ). □
In item 2, the sufficient condition ⟨𝑞 𝑘 , 𝑣⟩ > 0 for all 𝑣 ∈ S is equivalent to 𝑞 𝑘 ∈ S ★ \ S ⊥ . 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) 𝜁 (𝑡) -R -1 AR -1 𝜁 (𝑡) + R -1 𝐺 (𝑡) ∈ -N Φ (𝜁 (𝑡))
with 𝜁 = R𝜉. This fits within the class of DIs analysed for instance in [START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF]. Under basic assumption on the regularity of 𝐺 (•), ( . Then:

(5.9)

⟨𝑞 𝑘 -P ⊤ 𝜉 0 , 𝑣⟩ = ⟨-P ℎ 𝑥 1,𝑘 𝜆 𝑘 + ℎ𝐺 𝑘 , 𝑣⟩ -⟨𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 0 , 𝑣 𝜆 ⟩ -⟨𝑋 1 𝑥 0 , 𝑣 𝑥 ⟩ = ⟨- 𝑋 1 + ℎ𝑋 1 𝐴 1 ℎ𝑋 1 𝐵 1 -ℎ𝐶 1 𝐵 ⊤ 3 𝑋 3 𝐵 3 -ℎ D 𝑥 1,𝑘 𝜆 𝑘 + ℎ𝐺 𝑘 , 𝑣⟩ -⟨𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 0 , 𝑣 𝜆 ⟩ -⟨𝑋 1 𝑥 0 , 𝑣 𝑥 ⟩. Since 𝜉 0 ∈ Φ it follows that 𝑥 0 ∈ ℝ 𝑛 1 while 𝜆 0 ∈ ℝ 𝑚 + , and similarly 𝑣 𝑥 ∈ ℝ 𝑛 1 , 𝑣 𝜆 ∈ ℝ 𝑚 + . Since 𝑣 ≠ 0 let us take 𝑥 0 = -𝛼𝑋 -1 1 𝑣 𝑥 ∥𝑣 𝑥 ∥ -𝛽𝑋 -1 1 𝑣 𝑥 ∥𝑣 𝑥 ∥ 2 ∥𝑣 𝜆 ∥
for some 𝛼 ≥, 𝛽 ≥ 0, and 𝜆 0 = 0. It follows that:

(5.10)

⟨𝑞 𝑘 -P ⊤ 𝜉 0 , 𝑣⟩ = -⟨ 𝑋 1 + ℎ𝑋 1 𝐴 1 ℎ𝑋 1 𝐵 1 -ℎ𝐶 1 𝐵 ⊤ 3 𝑋 3 𝐵 3 -ℎ D 𝑥 1,𝑘 𝜆 𝑘 , 𝑣⟩ + ⟨ℎ𝐺 𝑘 , 𝑣⟩ ≥ -∥ 𝑋 1 + ℎ𝑋 1 𝐴 1 ℎ𝑋 1 𝐵 1 -ℎ𝐶 1 𝐵 ⊤ 3 𝑋 3 𝐵 3 -ℎ D 𝑥 1,𝑘 𝜆 𝑘 ∥ ∥𝑣 ∥ -∥ℎ𝐺 𝑘 ∥ ∥𝑣 ∥ + 𝛼 ∥𝑣 𝑥 ∥ + 𝛽 ∥𝑣 𝜆 ∥ ≥ -∥ 𝑋 1 + ℎ𝑋 1 𝐴 1 ℎ𝑋 1 𝐵 1 -ℎ𝐶 1 𝐵 ⊤ 3 𝑋 3 𝐵 3 -ℎ D 𝑥 1,𝑘 𝜆 𝑘 ∥ (∥𝑣 𝑥 ∥ + ∥𝑣 𝜆 ∥) -∥ℎ𝐺 𝑘 ∥ (∥𝑣 𝑥 ∥ + ∥𝑣 𝜆 ∥) + 𝛼 ∥𝑣 𝑥 ∥ + 𝛽 ∥𝑣 𝜆 ∥.
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

Since 𝐺 (•) is bounded, for large enough 𝛼 and 𝛽 it follows that ⟨𝑞 𝑘 -P ⊤ 𝜉 0 , 𝑣⟩ > 0 and this holds for any 𝑣 ≠ 0 inside S. □

boundedness and convergence analysis

In the next analysis we assume that 𝐵 ⊤ 3 𝑋 3 𝐵 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 𝐾 = ℝ 𝑚 + and that 𝑋 1 ≻ 0. From the structure of P and of the set-valued right-hand side, it follows from (5.5) that 𝑥 1,𝑘+1 is calculable from 𝑥 1,𝑘 and 𝜆 𝑘 :

(5.11)

𝑋 1 𝑥 1,𝑘+1 = (𝑋 1 + ℎ𝑋 1 𝐴 1 )𝑥 1,𝑘 + ℎ𝑋 1 𝐵 1 𝜆 𝑘 -ℎ𝐺 𝑥 𝑘 ⇔ 𝑥 1,𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 + ℎ𝐵 1 𝜆 𝑘 -ℎ𝑋 -1 1 𝐺 𝑥 𝑘 , for all 𝑘 ≥ 0.
where

𝐺 𝑘 = 𝐺 𝑥 𝑘 𝐺 𝜆 𝑘 .
Lemma 5.3. From (5.11), for all 𝑘 ≥ 0 one has:

(5.12)

𝑥 1,𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘+1 𝑥 1,0 + ℎ 𝑘 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘 -𝑖 𝐵 1 𝜆 𝑖 -ℎ 𝑘 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘 -𝑖 𝑋 -1 1 𝐺 𝑥 𝑖 .
Proof. The relation is true for 𝑘 = 0. Let it be true for 𝑘 ≥ 1. Then using (5.11) it follows that:

(5.13)

𝑥 1,𝑘+2 = (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘+1 + ℎ𝐵 1 𝜆 𝑘+1 -ℎ𝑋 -1 1 𝐺 𝑥 𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 ){(𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘+1 𝑥 1,0 + ℎ 𝑘 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘 -𝑖 𝐵 1 𝜆 𝑖 -ℎ 𝑘 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘 -𝑖 𝑋 -1 1 𝐺 𝑥 𝑖 } + ℎ𝐵 1 𝜆 𝑘+1 -ℎ𝑋 -1 1 𝐺 𝑥 𝑘+1 = (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘+2 𝑥 1,0 + ℎ 𝑘+1 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘+1-𝑖 𝐵 1 𝜆 𝑖 -ℎ 𝑘+1 ∑ 𝑖=0 (𝐼 𝑛 1 + ℎ𝐴 1 ) 𝑘+1-𝑖 𝑋 -1 1 𝐺 𝑥 𝑖 .
By induction the relation in (5.12) is true for all 𝑘 ≥ 0. □ Using (5.5) we can also express 𝜆 𝑘+1 as the solution of a generalised equation (GE):

(5.14)

𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 𝑘+1 + (ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 )𝜆 𝑘 + ℎ𝐶 1 𝑥 1,𝑘 + ℎ𝐺 𝜆 𝑘 Δ =𝑞 𝜆 𝑘 (𝜆 𝑘 ,𝑥 1,𝑘 ,𝐺 𝜆 𝑘 ) ∈ -N ℝ 𝑚 + (𝜆 𝑘+1 ) ⇕ 𝜆 𝑘+1 ∈ (𝐵 ⊤ 3 𝑋 3 𝐵 3 + N ℝ 𝑚 + ) -1 (𝐵 ⊤ 3 𝑋 3 𝐵 3 -ℎ D)𝜆 𝑘 -ℎ𝐶 1 𝑥 1𝑘 -
) ∥ ≤ 𝜂 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ ∥𝜆 min 𝑘 ∥ + ℎ𝜂 ∥𝐶 1 𝑥 1,𝑘 ∥ + ℎ𝜂 ∥𝐺 𝜆 𝑘 ∥. Let us denote ℎ 3 Δ = ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥.
Lemma 5.4. One has for all 𝑘 ≥ 0:

(5.15)

∥𝜆 min 𝑘+2 (𝜆 𝑘+1 , 𝑥 1,𝑘+1 , 𝐺 𝜆 𝑘+1 ) ∥ ≤ (𝜂ℎ 3 ) 𝑘+1 ∥𝜆 min 1 (𝜆 0 , 𝑥 1,0 , 𝐺 𝜆 0 ) ∥ + ℎ𝜂 𝑘 ∑ 𝑖=0 (𝜂ℎ 3 ) 𝑘 -𝑖 (∥𝐶 1 𝑥 1,𝑖+1 ∥ + ∥𝐺 𝜆 𝑖+1 ∥)
Proof. The inequality is satisfied for 𝑘 = 0. Assume that it is true for some 𝑘 ≥ 1. Then (5.16)

∥𝜆 min 𝑘+3 (𝜆 𝑘+2 , 𝑥 1,𝑘+2 , 𝐺 𝜆 𝑘+2 ) ∥ ≤ 𝜂 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ ∥𝜆 min 𝑘+2 ∥ + ℎ𝜂 ∥𝐶 1 𝑥 1,𝑘+2 ∥ + ℎ𝜂 ∥𝐺 𝜆 𝑘+2 ∥ ≤ 𝜂 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝜂 𝑘+1 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘+1 ∥𝜆 min 1 (𝜆 0 , 𝑥 1,0 , 𝐺 𝜆 0 ) ∥ + ℎ𝜂 𝑘 ∑ 𝑖=0 𝜂 𝑘 -𝑖 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘 -𝑖 (∥𝐶 1 𝑥 1,𝑖+1 ∥ + ∥𝐺 𝜆 𝑖+1 ∥) + ℎ𝜂 ∥𝐶 1 𝑥 1,𝑘+2 ∥ + ℎ𝜂 ∥𝐺 𝜆 𝑘+2 ∥ = 𝜂 𝑘+2 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘+2 ∥𝜆 min 1 (𝜆 0 , 𝑥 1,0 , 𝐺 𝜆 0 ) ∥ + ℎ𝜂 𝑘 ∑ 𝑖=0 𝜂 𝑘+1-𝑖 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘+1-𝑖 (∥𝐶 1 𝑥 1,𝑖+1 ∥ + ∥𝐺 𝜆 𝑖+1 ∥) + ℎ𝜂 ∥𝐶 1 𝑥 1,𝑘+2 ∥ + ℎ𝜂 ∥𝐺 𝜆 𝑘+2 ∥ = 𝜂 𝑘+2 ∥ℎ D -⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘+2 ∥𝜆 min 1 (𝜆 0 , 𝑥 1,0 , 𝐺 𝜆 0 ) ∥ + ℎ𝜂 𝑘+1 ∑ 𝑖=0 𝜂 𝑘+1-𝑖 ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ 𝑘+1-𝑖 (∥𝐶 1 𝑥 1,𝑖+1 ∥ + ∥𝐺 𝜆 𝑖+1 ∥).

The result follows by induction. □

The next step is to use Lemmas 5. = ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥. Then using (5.12) and (5.15) one has:

(5.17)

                 ∥𝑥 1,𝑘+1 ∥ ≤ ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ + ℎ 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 ∥𝐵 1 ∥ ∥𝜆 min 𝑖+1 ∥ + ℎ 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 ∥𝑋 -1 1 ∥ ∥𝐺 𝑥 𝑖 ∥ ∥𝜆 min 𝑘+2 ∥ ≤ 𝜂 𝑘+1 ℎ 𝑘+1 3 ∥𝜆 min 1 ∥ + ℎ𝜂 𝑘 ∑ 𝑖=0 𝜂 𝑘 -𝑖 ℎ 𝑘 -𝑖 3 (∥𝐶 1 𝑥 1,𝑖+1 ∥ + ∥𝐺 𝜆 𝑖+1 ∥)
Therefore we obtain:

(5.18)

∥𝑥 ,1,𝑘+1 ∥ ≤ ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ + ℎ 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 ∥𝑋 -1 1 ∥ ∥𝐺 𝑥 𝑖 ∥ + ℎ 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 ∥𝐵 1 ∥𝜂 𝑖 ℎ 𝑖 3 ∥𝜆 min 1 ∥ + ℎℎ 𝑘 1 ∥𝐵 1 ∥ ∥𝜆 min 1 ∥ + ℎ 𝑘 ∑ 𝑖=1 ℎ 𝑘 -𝑖 1 ∥𝐵 1 ∥ ℎ𝜂 𝑖 -1 ∑ 𝑗=0 𝜂 𝑖 -1-𝑗 ℎ 𝑖 -1-𝑗 3 ∥𝐶 1 ∥ ∥𝐺 𝜆 𝑗+1 ∥ + ℎ 𝑘 ∑ 𝑖=1 ℎ 𝑘 -𝑖 1 ∥𝐵 1 ∥ℎ𝜂 𝑖 -1 ∑ 𝑗=0 𝜂 𝑖 -1-𝑗 ℎ 𝑖 -1-𝑗 3 ∥𝐶 1 ∥ ∥𝑥 1,𝑗+1 ∥,
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . which we rewrite as:

(5.19) ∥𝑥 1,𝑘+1 ∥ ≤ ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ + ℎ 𝑓 𝑘 + ℎ 2 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 ∑ 𝑖=1 ℎ 𝑘 -𝑖 1 𝑖 -1 ∑ 𝑗=0 (𝜂ℎ 3 ) 𝑖 -1-𝑗 ∥𝑥 1,𝑗+1 ∥ with (5.20) 𝑓 𝑘 = 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 ∥𝑋 -1 1 ∥ ∥𝐺 𝑥 𝑖 ∥ + ∥𝐵 1 ∥ ∥𝜆 min 1 ∥ 𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 (𝜂ℎ 3 ) 𝑖 + ℎ 𝑘 1 ∥𝐵 1 ∥ ∥𝜆 min 1 ∥ + ℎ𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥ 𝑘 ∑ 𝑖=1 ℎ 𝑘 -𝑖 1 𝑖 -1 ∑ 𝑗=0 (𝜂ℎ 3 ) 𝑖 -1-𝑗 ∥𝐺 𝜆 𝑗+1 ∥.
The third term in the right-hand side of (5.19) involves terms from ∥𝑥 1,1 ∥ to ∥𝑥 1,𝑘 ∥. The inequality (5. [START_REF] Brézis | Image d'une somme d'opérateurs monotones et applications[END_REF]) can be rewritten as:

(5.21)

∥𝑥 1,𝑘+1 ∥ ≤ ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ + ℎ 𝑓 𝑘 + ℎ 2 𝑘 ∑ 𝑖=1 𝛽 𝑖 ∥𝑥 1,𝑖 ∥,
for some positive sequence of real numbers {𝛽 𝑖 }, each 𝛽 𝑖 being calculated from the expression in (5. [START_REF] Brézis | Image d'une somme d'opérateurs monotones et applications[END_REF])

𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 𝑖=1 ℎ 𝑘 -𝑖 1 𝑖 -1 𝑗=0 (𝜂ℎ 3 ) 𝑖 -1-𝑗 ∥𝑥 1,𝑗+1
𝛽 𝑖 = 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 ∑ 𝑗=𝑖 ℎ 𝑘 -𝑗 1 (𝜂ℎ 3 ) 𝑗 -𝑖 , 𝑖 ≥ 1.
Let us now apply a version of the discrete-time Gronwall inequality, due to Bachpatte [START_REF] Bachpatte | On the discrete generalisations of Gronwall's inequality[END_REF] [47, Lemma 100], see Appendix e. It is deduced that for all 𝑘 ≥ 0:

(5.23)

∥𝑥 1,𝑘+1 ∥ ≤ ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ + ℎ 𝑓 𝑘 + ℎ 2 𝑘 ∑ 𝑖=0 (ℎ 𝑖+1 1 𝛽 𝑖 ∥𝑥 1,0 ∥ + ℎ𝛽 𝑖 𝑓 𝑖 ) Π 𝑘+1 𝑗=𝑖+1 (1 + ℎ 2 𝛽 𝑗 ).
Let us state the next assumption (recall that ℎ 3 = ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ and ℎ 1 = ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥). Assumption 4. One has a) ℎ 1 ≤ 1 for ℎ > 0 small enough, b) ∥𝐺 (𝑡) ∥ is uniformly bounded by 𝐺 max < +∞ on [0,𝑇 ] for any 𝑇 > 0, and 𝐺 (•) is Lipschitz continuous with Lipschitz constant 𝑘 𝐺 ≥ 0, c) 𝜂ℎ 3 ≤ 1 for ℎ > 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 𝐵 ⊤ 3 𝑋 3 𝐵 3 = 0 (Assumption 2), then the minimal-norm element of the set of solutions to the LCP is always 𝜆 min 𝑘 = 0, hence 𝜂 = 0 in that case. Using (5.20) the next upper bound can be found:

(5.24)

𝑓 𝑘 ≤ (𝑘 + 1) (𝐺 max ∥𝑋 -1 1 ∥ + ∥𝐵 1 ∥ ∥𝜆 min 1 ∥) + ∥𝐵 1 ∥ ∥𝜆 min 1 ∥ + ℎ𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥ 𝐺 max 𝑘 ∑ 𝑖=1 𝑖 -1 ∑ 𝑗=0 𝑗,
Since ℎ = 𝑇 𝑁 and 𝑘 ∈ {0, 𝑁 -1}:

(5.25) ℎ𝑓 𝑘 ≤ 𝑇

(𝑘 + 1) 𝑁 (𝐺 max ∥𝑋 -1 1 ∥ + ∥𝐵 1 ∥ ∥𝜆 min 1 ∥) + 𝑇 𝑁 ∥𝐵 1 ∥ ∥𝜆 min 1 ∥ + 𝑇 2 𝑛 2 𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥ 𝐺 max 𝑘 ∑ 𝑖=1 𝑖 -1 ∑ 𝑗=0 𝑗 .
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

Now we have that

𝑘 𝑖=1 𝑖 -1 𝑗=0 𝑗 = (𝑘 -1) + 2(𝑘 -1) + 3(𝑘 -3) + . . . + (𝑘 -3)3 + (𝑘 -2)2 + (𝑘 -1)
, where the sum consists of 𝑘 -1 terms. Hence 𝑘 𝑖=1 𝑖 -1

𝑗=0 𝑗 = O (𝑘 2 ), and 𝑇 2 𝑁 2 𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥ 𝐺 max 𝑘 𝑖=1 𝑖 -1 𝑗=0 𝑗 = 𝑇 2 𝜂 ∥𝐵 1 ∥ ∥𝐶 1 ∥ 𝐺 max O (𝑘 2 ) 𝑁 2 .
Since 𝑘 ≤ 𝑁 we infer that this term is uniformly bounded in 𝑘, for each fixed bounded 𝑇 . The same holds for 𝑇 (𝑘+1) 𝑁 (𝐺 max ∥𝑋 -1 1 ∥ + ∥𝐵 1 ∥ ∥𝜆 min 1 ∥). Thus ℎ𝑓 𝑘 ≤ 𝛼 1 max(𝑇 ,𝑇 2 ) for some 𝛼 1 which depends on 𝜂, ∥𝐵 1 ∥, ∥𝐶 1 ∥, 𝐺 max , ∥𝜆 min 1 ∥. The next step consists of computing an upper bound from (5.23) as:

(5.26) ∥𝑥 1,𝑘+1 ∥ ≤ ∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝑇 2 𝑁 2 𝑘 ∑ 𝑖=0 𝛽 𝑖 Π 𝑘+1 𝑗=𝑖+1 (1 + 𝑇 2 𝑁 2 𝛽 𝑗 ) .
In view of Assumption 4 and (5.22), one has

𝛽 𝑖 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 𝑗=𝑖 1 = 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ (𝑘 -𝑖) ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥𝑁 , 𝑖 ≥ 1, 𝑘 ≤ 𝑁 .
Consequently:

(5.27)

∥𝑥 1,𝑘+1 ∥ ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝑇 2 𝑁 𝑘 ∑ 𝑖=0 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ Π 𝑘+1 𝑗=𝑖+1 (1 + 𝑇 2 𝑁 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥)) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 𝑘 ∑ 𝑖=0 Π 𝑘+1 𝑗=𝑖+1 (1 + 𝑇 2 𝑁 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥)) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 𝑘 ∑ 𝑖=0 (1 + 𝑇 2 𝑁 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥) 𝑘 -𝑖 ) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 𝑘 ∑ 𝑖=0 (1 + 𝛼 2 𝑁 ) 𝑘 -𝑖 ) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 𝑁 ∑ 𝑖=0 (1 + 𝛼 2 𝑁 ) 𝑁 ) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 (𝑁 + 1) (1 + 𝛼 2 𝑁 ) 𝑁 ) ≤ (∥𝑥 1,0 ∥ + 𝛼 1 max(𝑇 ,𝑇 2 )) (1 + 𝑇 2 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑁 + 1 𝑁 (1 + 𝛼 2 𝑁 ) 𝑁 ),
where [73, p.196]. Therefore the sequence {𝑥 1,𝑘 } 𝑘 is uniformly bounded over 𝑘 ∈ {0, 𝑁 -1}. It directly follows from Lemma 5.4 that the sequence {𝜆 min 𝑘 } 𝑘 is uniformly bounded as well. It follows from (4.1) that the sequence

𝛼 2 = 𝑇 2 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥. One has (1 + 𝛼 2 𝑁 ) 𝑁 → 𝑒 𝛼 2 as 𝑁 → +∞
𝑥 1,𝑘+1 -𝑥 1,𝑘 ℎ
𝑘 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

𝑥 ℎ★ 1 (𝑡) = 𝑥 1,𝑘+1 if 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ), 𝜆 min ℎ★ (𝑡) = 𝜆 min 𝑘 if 𝑡 ∈ [𝑡 𝑘 , 𝑡 𝑘+1 ).
Proposition 5.5. Let Assumption 4 hold true, as well as the conditions for the well-posedness of the OSNSP. Then the sequence {𝑥 ℎ 1 (•)} ℎ≥0 converges uniformly towards a continuous limit

𝑥 1 (•). Moreover 𝑥 ℎ★ 1 (•) converges strongly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) towards 𝑥 1 (•), and 𝑥 ℎ 1 (•) converges weakly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) towards 𝑥 1 (•). Therefore {(𝐵 1 𝜆 min ℎ★ ) (•)} ℎ≥0 converges weakly in L 2 ( [0,𝑇 ]; ℝ 𝑛 1 ) towards a limit 𝜆 1 (•), and 𝑥 1 (𝑡) = 𝐴 1 𝑥 1 (𝑡) + 𝐺 𝑥 (𝑡) + 𝜆 1 (𝑡) almost everywhere on [0,𝑇 ]. Proof. Both sequences {𝑥 ℎ 1 (•)} and { 𝑥 ℎ 1 (•)} are uniformly bounded on [0,𝑇 ], the sequence {𝑥 ℎ 1 (•)} of C 0 ([0,𝑇 ]; ℝ 𝑛 1 )
functions is thus equicontinuous. From the Ascoli-Arzela Theorem, {𝑥 ℎ 1 (•)} stays in a compact subset of C 0 ( [0,𝑇 ]; ℝ 𝑛 1 ) and indeed it converges uniformly towards a continuous limit 𝑥 1 (•). Let us now examine the sequences of piecewise-linear functions {𝐺 ℎ (•)} ℎ≥0 and { 𝐺 ℎ (•)} ℎ≥0 , Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . designed similarly as (4.7). From Assumption 4 b) one infers that 𝐺 ℎ (𝑡) is uniformly bounded on [0,𝑇 ], ∥𝐺 𝑘+1 -𝐺 𝑘 ∥ ℎ ≤ 𝑘 𝐸 implies that 𝐺 ℎ (𝑡) is also uniformly bounded in ℎ on [0,𝑇 ], thus {𝐺 ℎ (•)} ℎ≥0 is a sequence of Lipschitz continuous, uniformly bounded, equicontinuous functions. From the Ascoli-Arzela Theorem, {𝐺 ℎ (•)} stays in a compact subset of C 0 ( [0,𝑇 ]; ℝ 𝑛 1 ) and it converges uniformly towards the continuous limit 𝐺 (•) as ℎ → 0. The remaining part of the proof is similar to the proof of Corollary 4.6 (adding 𝐺 𝑥,ℎ (•) 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 𝐶 3 𝐵 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 {𝜆 min 𝑘 } 𝑘 and {𝑥 2,𝑘 } 𝑘 . 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 𝑥 1,max (𝑇 , 𝜂, 𝑥 1,0 ). Then from (5.17) and Assumption 4 we obtain:

(5.28)

∥𝜆 min 𝑘+1 ∥ ≤ ∥𝜆 min 0 ∥ + ℎ𝜂 𝑘 ∑ 𝑖=0 (∥𝐶 1 ∥𝑥 1,max (𝑇 , 𝜂, 𝑥 1,0 ) + ∥𝐺 𝜆 𝑖+1 ∥) ≤ ∥𝜆 min 0 ∥ + ℎ𝜂 ∥𝐶 1 ∥𝑥 1,max (𝑇 , 𝜂, 𝑥 1,0 ) (𝑘 + 1) + ℎ𝜂 𝑘 ∑ 𝑖=0 ∥𝐺 𝜆 𝑖+1 ∥ ≤ ∥𝜆 min 0 ∥ + 𝑇 (𝑘 + 1) 𝑁 𝜂 ∥𝐶 1 ∥𝑥 1,max (𝑇 , 𝜂, 𝑥 1,0 ) + 𝜂 𝑇 (𝑘 + 1) 𝑁 𝐺 max ≤ ∥𝜆 min 0 ∥ + 𝑇 (𝑁 + 1) 𝑁 𝜂 ∥𝐶 1 ∥𝑥 1,max (𝑇 , 𝜂, 𝑥 1,0 ) + 𝜂 𝑇 (𝑁 + 1) 𝑁 𝐺 max .
Therefore ∥𝜆 min 𝑘+1 ∥ is uniformly bounded for any bounded 𝑇 and initial data. The boundedness of 𝑥 2,𝑘 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 𝑧 in the complementarity variable 𝑤. 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 ℎ 𝑘 1 another way, since for any 𝑎 1 > 0 one has (1 + 𝑇 𝑁 𝑎 1 ) 𝑁 → 𝑒 𝑎 1 𝑇 as 𝑁 → +∞, 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 ℎ𝑓 𝑘 are upper bounded by constants times (5.29) ℎ

𝑘 ∑ 𝑖=0 ℎ 𝑘 -𝑖 1 = ℎ 𝑘 ∑ 𝑖=0 ℎ 𝑖 1 ≤ ℎ𝑘 (1 + ℎ𝑎 1 ) 𝑘 ,
where ℎ 1 ≤ ∥𝐼 𝑛 1 ∥ + ℎ∥𝐴 1 ∥ = 1 + ℎ𝑎 1 for some 𝑎 1 ≥ 0. We have ℎ𝑘 (1 + ℎ𝑎 1 ) 𝑘 ≤ 𝑇 𝑁 𝑘 (1 + 𝑇 𝑎 1 𝑁 ) 𝑁 → 𝑇 𝑒 𝑎 1 𝑇 as 𝑁 → +∞. The third term of ℎ𝑓 𝑘 is upper bounded by a constant times ℎℎ 𝑘 1 ≤ 𝑇 𝑁 (1 + ℎ𝑎 1 ) 𝑁 → 0 as 𝑁 → +∞. The last term of ℎ𝑓 𝑘 is upper bounded by a constant multiplied by (5.30)

ℎ 2 𝑘 ∑ 𝑖=1 ℎ 𝑘 -𝑖 1 𝑘 ≤ 𝑇 2 𝑁 2 𝑘 2 (1 + ℎ𝑎 1 ) 𝑘 ≤ 𝑇 2 𝑁 2 𝑁 2 (1 + ℎ𝑎 1 ) 𝑁 → 𝑇 2 𝑒 𝑇 𝑎 1
as 𝑁 → +∞. These first calculations allow us to upper bound the term ℎ𝑓 𝑘 in (5.23). The first term

ℎ 𝑘+1 1 ∥𝑥 1,0 ∥ = (1 + ℎ𝑎 1 ) 𝑘+1 ∥𝑥 1,0 ∥ ≤ (1 + ℎ𝑎 1 ) 𝑁 ∥𝑥 1,0 ∥ → ∥𝑥 1,0 ∥𝑒 𝑎 1 𝑇 as 𝑁 → +∞.
The third term in (5.23) is more complex:

ℎ 2 𝑘 𝑖=0 (ℎ 𝑖+1 1 𝛽 𝑖 ∥𝑥 1,0 ∥ + ℎ𝛽 𝑖 𝑓 𝑖 ) Π 𝑘+1 𝑗=𝑖+1 (1 + ℎ 2 𝛽 𝑗 ). Since from Assumption 4 c) 𝜂ℎ 3 ≤ 1, it follows that (5.31) 𝛽 𝑖 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 ∑ 𝑗=1 ℎ 𝑘 -𝑗 1 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑘 ∑ 𝑗=1 (1 + ℎ𝑎 1 ) 𝑘 -𝑗
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

for some 𝑎 1 > 0. Thus

(5.32) ℎ 2 𝛽 𝑗 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ℎ 2 𝑘 (1 +ℎ𝑎 1 ) 𝑘 -1 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 2 𝑁 (1 + 𝑇 𝑁 𝑎 1 ) 𝑁 = 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 (1 + 𝑇 𝑁 𝑎 1 ) 𝑁 .

Hence

(5.33)

Π 𝑘+1 𝑗=𝑖+1 (1 + ℎ 2 𝛽 𝑗 ) ≤ (1 + 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 (1 + 𝑇 𝑁 𝑎 1 ) 𝑁 ) 𝑁 → 𝑒 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥𝑇 2 𝑒 𝑎 1 𝑇 .
We have

(5.34) ℎ 2 𝑘 ∑ 𝑖=0 𝛽 𝑖 ℎ𝑓 𝑖 = 𝑘 ∑ 𝑖=0 (ℎ 2 𝛽 𝑖 ) (ℎ𝑓 𝑖 ) ≤ 𝑁 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ 𝑇 2 𝑁 (𝛼 1 𝑇 + 𝛼 2 𝑇 𝑁 + 𝛼 3 𝑇 2 ) (1 + 𝑇 𝑁 𝑎 1 ) 2𝑁 → 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥𝑇 2 (𝛼 1 𝑇 + 𝛼 3 𝑇 2 )𝑒 2𝑎 1 𝑇 ,
hence is bounded for any bounded 𝑇 . We have also

(5.35) ℎ 2 𝑘 ∑ 𝑖=0 ℎ 𝑖+1 1 𝛽 𝑖 ∥𝑥 1,0 ∥ = ∥𝑥 1,0 ∥ 𝑘 ∑ 𝑖=0 ℎ 𝑖+1 1 ℎ 2 𝛽 𝑖 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ ∥𝑥 1,0 ∥ 𝑘 ∑ 𝑖=0 (1 + ℎ𝑎 1 ) 𝑘 ℎ 2 𝑘 (1 + ℎ𝑎 1 ) 𝑘 -1 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ ∥𝑥 1,0 ∥ 𝑁 ∑ 𝑖=0 (1 + ℎ𝑎 1 ) 2𝑁 𝑇 2 𝑁 2 𝑁 ≤ 𝜂 ∥𝐶 1 ∥ ∥𝐵 1 ∥ ∥𝑥 1,0 ∥ 𝑇 2 𝑒 2𝑎 1 𝑇 .
Consequently from (5.23) it is inferred that if Assumption 4 b) and c) holds, then ∥𝑥 1,𝑘+1 ∥ is bounded by a term that depends only on 𝑇 , initial conditions and constants. Let us examine now Assumption 4 c) and the terms which which involve ℎ 𝑘 3 . We have

ℎ 3 = ∥ℎ D - 𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥. Let us assume that D is invertible. Then (5.36) ∥ℎ D -𝐵 ⊤ 3 𝑋 3 𝐵 3 ∥ = ∥ D (ℎ -D -1 𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ ≤ ∥ D ∥ ∥ (ℎ𝐼 𝑚 -D -1 𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ ≤ ∥ D ∥ (ℎ + ∥ D -1 𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥) ≤ ∥ D ∥ (ℎ + ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥) = ∥ D ∥ ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ 1 + 1 ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ ℎ . Therefore it follows that ℎ 𝑘 3 ≤ (∥ D ∥ ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥) 𝑘 1 + 𝑇 ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ 1 𝑁 𝑘 . Consequently (𝜂ℎ 3 ) 𝑘 ≤ (𝜂 ∥ D ∥ ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥) 𝑘 1 + 𝑇 ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ 1 𝑁 𝑘 which converges as 𝑁 → +∞ to (𝜂 ∥ D ∥ ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥) 𝑁 𝑒 𝑇 ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ . A necessary condition for the limit to be bounded is 𝜂 ∥ D ∥ ∥ D -1 ∥ ∥𝐵 ⊤ 3 𝑋 3 𝐵 3 ) ∥ ≤ 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.

case 𝑋

1 ≻0 and 𝐵 ⊤ 3 𝑋 3 𝐵 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 [START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF] (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,𝑇 ]; ℝ 𝑛 1 +𝑚 ) with essentially bounded derivatives [START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF]Proposition 2.6]. It is deduced that 𝑥 2 (•) is also unique in C 0 ( [0,𝑇 ]; ℝ 𝑛 1 +𝑚 ) with essentially bounded derivative, while the algebraic variable 𝑧 (•) (which depends on 𝜆(•), see (5.1)) is essentially bounded (hence a priori it may be discontinuous).

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 [START_REF] Pan | Digital switched hydraulics[END_REF][START_REF] Guglielmino | Power hydraulics-Switched mode control of hydraulic actuation[END_REF] [52, Chapter 1], where check valves play the role of diodes.

first example

i 3 i 2 i 1 u D u 2 u 1 C R 3 R 1 R 2 L Figure 1: An 𝑅𝐿𝐶𝐷 passive circuit.
Let us consider the circuit in Figure 1. Let us assume first that 𝑅

1 = 𝑅 3 = 0 Ω. Let 𝜉 1 = 𝑖 2 , 𝜉 2 = ∫ 𝑡 0 𝑖 3 (𝑡)𝑑𝑡, 𝜉 3 = 𝑖 1 .
The dynamics is given by: (6.1)

                   𝜉 1 (𝑡) = - 𝑅 2 𝐿 𝜉 1 (𝑡) + 1 𝐿 𝑢 𝐷 + 1 𝐿 𝑢 1 (𝑡) 𝜉 2 (𝑡) = 𝜉 1 (𝑡) -𝜉 3 (𝑡) 0 = 1 C 𝜉 2 (𝑡) + 𝑢 𝐷 (𝑡) + 𝑢 2 (𝑡) 0 ≤ 𝑢 𝐷 ⊥ 𝜉 3 (𝑡) ≥ 0.
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 𝑥 = (𝜉 1 , 𝜉 2 ) ⊤ ,

𝑧 = 𝜉 3 , 𝐶 = (0 0 1), 𝐵 = ( 1 𝐿 0 1) ⊤ . The transfer function 𝐻 (𝑠) = 𝑧 (𝑠 ) 𝑢 𝐷 (𝑠 ) = C𝑠 2 + C𝑅 𝐿 𝑠+ 1 𝐿 𝑠+ 𝑅 𝐿 = C𝑠 + 1 𝐿𝑠+𝑅 .
It has a relative degree 𝑟 = -1. It can be checked that 𝐻 (𝑠) is positive real (PR) as long as 𝑅 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((𝐴 𝐵) -𝑠𝑃 𝐵) = rank -𝑅 𝐿 -𝑠 0 0 1 𝐿 1 -𝑠 -1 0 0 1 C 0 1 = 3 for all 𝑠 ∈ ℂ. 2. rank(𝑃 𝐵) = rank 1 0 0 1 𝐿 0 1 0 0 0 0 0 1 = 3.

Brogliato, Rocca

Analysis of the implicit Euler time-discretization of passive . . .

rank(𝐴

⊤ -𝑠𝑃 ⊤ 𝐶 ⊤ ) = rank -𝑅 𝐿 -𝑠 1 0 0 0 -𝑠 1 C 0 0 -1 0 1 = 3 for all 𝑠 ∈ ℂ. 4. rank(𝑃 ⊤ 𝐶 ⊤ ) = rank 1 0 0 0 0 1 0 0 0 0 0 1 = 3. 5. (𝐴 𝐵) Ker(𝑃) = -𝑅 𝐿 0 0 1 0 -1 0 1 C 0 Ker(𝑃) = {𝜉 ∈ ℝ 3 | 𝜉 1 = 𝜉 3 = 0, 𝜉 2 ∈ ℝ} ⊂ Im(𝑃) = {𝜉 ∈ ℝ 3 | 𝜉 3 = 0, 𝜉 2 ∈ ℝ, 𝜉 1 ∈ ℝ},
Therefore we infer that (6.1) admits a special Weiertrass form as (2.2). Let us choose 𝐿 = 4. . Let us make the state variable change 𝜁 = 𝑇 -1 𝜉, then we obtain the transformed dynamics: (

               𝜁 1 (𝑡) = -2.12 10 3 𝜁 1 (𝑡) -3 10 -3 𝜆(𝑡) + 0.301 10 4 𝑢 1 (𝑡) 𝜁 3 (𝑡) = 𝜁 2 (𝑡) 0 = 𝜁 3 (𝑡) + 10 -4 𝜆(𝑡) + 10 -4 𝑢 2 (𝑡) 0 ≤ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = -0.707𝜁 1 (𝑡) -𝜁 2 (𝑡) ≥ 0. 6.2)
Let us now define 𝑥 1 = 𝜁 1 , 𝑥 2 = 𝜁 3 , 𝑥 3 = 𝜁 2 , then we obtain (6.3)

            
𝑥 1 (𝑡) = -2.12 10 3 𝑥 1 (𝑡) -3 10 -3 𝜆(𝑡) + 0.301 10 4 𝑢 1 (𝑡)

𝑥 2 (𝑡) = 𝑥 3 (𝑡) 0 = 𝑥 2 (𝑡) + 10 -4 𝜆(𝑡) + 10 -4 𝑢 2 (𝑡) 0 ≤ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = -0.707𝑥 1 (𝑡) -𝑥 3 (𝑡) ≥ 0, which is the canonical form in (2.2) with 𝑧 = 𝑥 3 , 𝐵 1 = -3.10 -3 , 𝐵 2 = 0, 𝐵 3 = 10 -4 , 𝐶 1 = -0.707, 𝐶 2 = 0, 𝐶 3 = -1, 𝐷 = 0. We have 𝑋 3 = 10 4 , 𝐴 1 = -2.12 10 3 , 𝑋 1 = 𝐶 1 𝐵 1 > 0, D = 𝐷 -𝐶 3 𝐵 2 -𝐶 2 𝐵 3 = 0, Dℎ = ℎ 2 𝐶 1 𝐵 1 + 𝐵 ⊤ 3 𝑋 3 𝐵 3 = 2.121 10 -3 ℎ 2 + 10 -4 , ℎ 3 = 𝐵 ⊤ 3 𝑋 3 𝐵 3 = 10 -4 , ℎ 1 = |1 + ℎ𝐴 1 | = |1 -2.
12 10 3 ℎ| which verifies Assumption 4 for ℎ > 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:

(6.4)                𝑥 1 (𝑡) = -2.12 10 3 𝑥 1 (𝑡) -3 10 -3 𝜆(𝑡) + 0.301 10 4 𝑢 1 (𝑡) 10 -4 𝜆(𝑡) = 0.707𝑥 1 (𝑡) -10 -4 𝑢 2 (𝑡) + 𝑤 (𝑡) 𝑥 2 (𝑡) = -10 -4 𝜆(𝑡) -10 -4 𝑢 2 (𝑡)
𝑧 (𝑡) = -10 -4 𝜆(𝑡) -10 -4 𝑢 2 (𝑡) .

From (5.14), the following LCP is obtained: ). Since D = 0 it is easy to show that 𝜂 = 10 4 while ℎ 3 = 10 -4 , using the LCP in (5.14) that can be solved explicitly. Thus 𝜂ℎ 3 ≤ 1 for all ℎ ≥ 0. However, ℎ 1 = ∥1 -2.12 10 3 ℎ∥ ≤ 1 if and only if ℎ < 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 𝑋 1 > 0 and 𝐵 ⊤ 3 𝑋 3 𝐵 3 = 10 -4 > 0. Hence the state space form (5.8), that is equivalent to (5.3) since P ≻ 0, can be calculated. Also the generalized equation (5.14) always has a unique solution 𝜆 𝑘+1 . As we shall see through numerical simulations, the variable 𝑧 (•) = 𝑥 3 (•) 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 𝑧, should rely on other tools than those used in this paper.

(6.5) 0 ≤ 10 -4 𝜆 𝑘+1 + 𝑞 𝑘 ⊥ 𝜆 𝑘+1 ≥ 0 𝑞 𝑘 = -10 -4 𝜆 𝑘 -0.707ℎ𝑥
In Figure 2, the results of the numerical simulation of (6.3) with an implicit Euler scheme are depicted 2 . Numerical simulations are performed with the INRIA software package Siconos [START_REF] Acary | An introduction to Siconos[END_REF]. At each time step, the Lemke solver [START_REF] Lemke | Equilibrium points of bimatrix games[END_REF][START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF] 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 • 10 -3 ) T , and the time step is ℎ = 5 • 10 -6 s. The input functions 𝑢 1 (𝑡) and 𝑢 2 (𝑡) are taken as sinusoidal functions: 𝑢 1 (𝑡) = sin(𝑡) and 𝑢 2 (𝑡) = cos(10𝑡). Let us observe that the variable 𝑧 (𝑡) = 𝑥 3 (𝑡) appears to be discontinuous at 𝑡 ≈ 5.9s. The multiplier 𝜆(𝑡) 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) error = ∥x(𝑡)x ℎ (𝑡) ∥ ∞ .

Remark 6.1 (modelling). One can obtain a relative degree one LCS by eliminating the variable 𝑢 𝐷 (𝑡) 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. Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . This circuit example depicted in Figure 5 is inspired from [START_REF] Camlibel | Differential-algebraic inclusions with maximal monotone operators[END_REF]. 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 𝑢 (𝑡) = 𝑢 𝑐 , while the other voltage source 𝑢 1 (𝑡) is time-varying Lipschitz continuous and bounded. Notice that the current/voltage characteristic of this device can be written as 𝑖 5 (𝑡) ∈ -𝜕𝜓 {𝑢 𝑐 } (𝑢 (𝑡)), or equivalently as ℝ ∋ 𝑖 5 (𝑡) ⊥ 𝑢 (𝑡) -𝑢 𝑐 ∈ {0}. The dynamics is given by: (6.7)

second example

L C R i 1 i 5 i 2 i 3 i 4 e(t) u D u(t) u 1 (t)
                       C 𝑒 (𝑡) = -𝑅 -1 𝑒 (𝑡) + 𝑖 4 (𝑡) -𝑖 1 (𝑡) -𝑖 5 (𝑡) 𝐿 𝑑𝑖 1 𝑑𝑡 (𝑡) = 𝑒 (𝑡) + 𝑢 1 (𝑡) 0 = 𝑒 (𝑡) + 𝑢 (𝑡) 0 ≤ 𝑢 𝐷 (𝑡) ⊥ 𝑖 4 (𝑡) ≥ 0 ℝ ∋ 𝑖 5 (𝑡) ⊥ 𝑢 (𝑡) -𝑢 𝑐 ∈ {0}.
The state vector is 𝜉 = (𝜉 1 , 𝜉 2 , 𝜉 3 ) ⊤ = (𝑒, 𝑖 1 , 𝑖 5 ) ⊤ , and 𝑢 𝐷 (𝑡) = 𝑒 (𝑡). The complementarity conditions can be written as 𝐾 = ℝ + × ℝ ∋ (𝑢 𝐷 (𝑡), 𝑖 5 (𝑡)) ⊥ (𝑖 4 (𝑡), 𝑢 (𝑡) -𝑢 𝑐 ) ∈ ℝ + × {0} = 𝐾 ★ . Setting 𝜆 1 = 𝑖 4 and

𝜆 2 = 𝑢 (𝑡) -𝑢 𝑐 , we obtain 𝑃 = C 0 0 0 𝐿 0 0 0 0 , 𝐴 = -𝑅 -1 -1 -1 1 0 0 1 0 0 , 𝐵 = 1 0 0 0 0 1 , 𝐶 = 1 0 0 0 0 1 , 𝐷 = 0,
and the LDCS form in (2.1) is equal to: (6.8)

                         𝜉 1 (𝑡) = - 1 𝑅C 𝜉 1 (𝑡) - 1 C 𝜉 2 (𝑡) - 1 C 𝜉 3 (𝑡) + 1 C 𝜆 1 (𝑡) 𝜉 2 (𝑡) = 1 𝐿 𝜉 1 (𝑡) + 𝑢 1 (𝑡) 𝐿 0 = 𝜉 1 (𝑡) + 𝜆 2 (𝑡) + 𝑢 𝑐 0 ≤ 𝜉 1 (𝑡) ⊥ 𝜆 1 (𝑡) ≥ 0 ℝ ∋ 𝜉 3 (𝑡) ⊥ 𝜆 2 (𝑡) ∈ {0}. Let 𝑆 = 0 1 𝐿 -1 𝐿 -1 0 -1 𝑅 0 0 -C ,𝑇 = 0 0 -1 C 1 0 0 -1 1 0 , then 𝑆𝐴𝑇 = 0 0 0 0 1 0 0 0 1 ,𝑆𝑃𝑇 = 1 𝐿 0 0 0 0 1 C 0 0 0 ,𝐶𝑇 = 0 0 -1 C -1 1 0 ,
Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

𝑆𝐵 = 0 -1 𝐿 1 -1 𝑅 0 -C
. Using the state variable change 𝜁 = 𝑇 -1 𝜉 𝑇 -1 = 0 1 0 0 1 1 -C 0 0 , one obtains: (6.9)

                         𝜁 1 (𝑡) = -(𝜆 2 (𝑡) + 𝑢 𝑐 (𝑡)) + 𝑢 1 (𝑡) 𝐿 1 C 𝜁 3 (𝑡) = 𝜁 2 (𝑡) + 𝜆 1 (𝑡) - 1 𝑅 (𝜆 2 (𝑡) + 𝑢 𝑐 (𝑡)) 0 = -C(𝜆 2 (𝑡) + 𝑢 𝑐 (𝑡)) + 𝜁 3 (𝑡) 𝐾 ★ ∋ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 0 -1 𝜁 1 (𝑡) + 0 1 𝜁 2 (𝑡) + -1 C 0 𝜁 3 (𝑡) ∈ 𝐾 . Let 𝑥 1 = 𝜁 1 , 𝑥 2 = 1 C 𝜁 3 , 𝑧 = 𝜁 2 , 𝐵 1 = (0 -1), 𝐵 2 = (1 -1 𝑅 ), 𝐵 3 = (0 -1), 𝐶 1 = 0 -1 , 𝐶 3 = 0 1 , 𝐶 2 = - 1 0
, 𝐷 = 0. This dynamics is rewritten in the canonical form as: (6.10)

                       𝑥 1 (𝑡) = -𝜆 2 (𝑡) -𝑢 𝑐 (𝑡) + 𝑢 1 (𝑡) 𝐿 𝑥 2 (𝑡) = 𝑧 (𝑡) + 𝜆 1 (𝑡) - 1 𝑅 (𝜆 2 (𝑡) + 𝑢 𝑐 (𝑡)) 0 = -𝜆 2 (𝑡) -𝑢 𝑐 (𝑡) + 𝑥 2 (𝑡) 𝐾 ★ ∋ 𝜆(𝑡) ⊥ 𝑤 (𝑡) = 0 -1 𝑥 1 (𝑡) + -1 0 𝑥 2 (𝑡) + 0 1 𝑧 (𝑡) ∈ 𝐾 .
In the simulations depicted in Figures 6a and6b, it has been chosen 𝑢 𝑐 (𝑡) = sin(10𝑡) -1 ≤ 0, which complies with the system's intrinsic constraints. Then (when 𝑢 1 (𝑡) = 0) one recovers the standard form in (2.2), with a different class of sets 𝐾, however. One has 𝐴 1 = 0, 𝐵 ⊤ 3 𝑋 3 = -𝐶 3 with 𝑋 3 = 1, and 𝑋 1 𝐵 1 = 𝐶 ⊤ 1 (this follows necessarily from (d.2) and 𝑄 1 ≼ 0 and using [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF]Lemma A.69]), with 𝑋 1 = 1 (hence the quadruple (𝐴 1 , 𝐵 1 , 𝐶 1 , D) is passive but it is not strictly passive and its associated transfer function is not strongly SPR). Also D

= -𝐶 3 𝐵 2 -𝐶 2 𝐵 3 = 0 1 -1 1 𝑅 ≽ 0. Also 𝐵 ⊤ 3 𝑋 3 𝐵 3 = 0 0 0 1 ≽ 0. It is computed that Dℎ = 0 ℎ -ℎ ℎ 𝑅 + ℎ 2 + 1
. The various sets for the set S in Proposition 3.2 are calculated as:

𝐾 * = ℝ + ×{0}, Ker( Dℎ + DT ℎ ) = 0 0 0 2( ℎ 𝑅 + ℎ 2 + 1) = ℝ×{0}, {𝑤 | Dℎ 𝑤 ∈ 𝐾 = ℝ + ×ℝ} = ℝ × ℝ + . So 𝑆 = ℝ + × {0}
≠ 0 and item 2 of Proposition 3.2 has to be applied. It is calculated that

⟨ℎ𝛽 𝑘 -Dℎ 𝜉 0 , 𝑣⟩ = ⟨ -ℎ𝜉 1 0 -ℎ𝑥 1,𝑘 -𝑥 2,𝑘 + ℎ𝜉 1 0 -( ℎ 2 + ℎ 2 + 1)𝜉 2 0
. Also 𝑣 ∈ 𝐾 ★ implies that 𝑣 1 ≥ 0 and 𝑣 2 = 0, so it follows that ⟨ℎ𝛽 𝑘 -Dℎ 𝜉 0 , 𝑣⟩ + 𝜓 𝐾 ★ (𝑣) > 0 is guaranteed, if 𝑣 1 > 0, by taking 𝜉 1 0 < 0. It is inferred that item 2 in Proposition 3.2 is verified and the OSNSP has at least one solution.

Any two solutions 𝜆

ℝ + × {0} ∋ 𝜆 1,𝑘+1 𝜆 2,𝑘+1 ⊥ 0 𝜆 2,𝑘+1 + 𝑞 𝜆 𝑘 ∈ ℝ + × ℝ.
Therefore 𝜆 2,𝑘+1 = 0, and 𝜆 1,𝑘+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.

third example

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

                   𝑥 1 (𝑡) = -2𝑥 1 (𝑡) + 2𝜆 1 (𝑡) -𝜆 2 (𝑡) 𝑥 2 (𝑡) = 𝑧 (𝑡) 0 = 𝑥 2 (𝑡) + 𝜆 1 (𝑡) -2𝜆 2 (𝑡) 0 ≤ 𝜆 1 (𝑡) 𝜆 2 (𝑡) ⊥ 2 -1 𝑥 1 (𝑡) + 1 3 𝑥 2 (𝑡) + -1 2
𝑧 (𝑡) + 2 -1 0 1 𝜆 1 (𝑡) 𝜆 2 (𝑡) ≥ 0
It can be verified that the transfer matrix 𝐻 (𝑠) = 𝐷 + 𝐶 (𝑠𝑃 -𝐴) It can be checked that the set S = {0} in Proposition 3.2 for any ℎ > 0, thus there exists a unique solution 𝜆 𝑘+1 to the implicit Euler discretization of (6.13). However 𝐶 3 ≠ 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:

= -2 < 0, 𝐶 1 = 𝑋 1 𝐵 T 1 = (2, -1) T with 𝑋 1 = 1 > 0, 𝑋 3 = 𝐼 2 , 𝐵 3 = (1 -2) = -𝐶 ⊤ 3 ,
(6.14)              1 0 0 0 1 -2 0 -2 4 𝑥 1 (𝑡) 𝜆 1 (𝑡) 𝜆 2 (𝑡) = -2 2 -1 -2 -1 -1 1 3 -7 𝑥 1 (𝑡) 𝜆 1 (𝑡) 𝜆 2 (𝑡) + 0 𝑤 1 (𝑡) 𝑤 2 (𝑡) 𝑤 (𝑡) ∈ -N ℝ 2 + (𝜆(𝑡))
Let us check if Assumption 4 (c) is verified. Using the Frobenius norm, one can obtain that 0 ≤ ℎ 3 = √ 60ℎ 2 + 25 -66ℎ ≤ 1 if 0 < ℎ ≤ 0.7. In the context of (6.14), the LCP (5.14) becomes:

(6.15) 0 ≤ 𝜆 1,𝑘+1 𝜆 2,𝑘+1 ⊥ 1 -2 -2 4 𝜆 1,𝑘+1 𝜆 2,𝑘+1 + (ℎ -1) (ℎ + 3) (2 -3ℎ) (7ℎ -4) 𝜆 1,𝑘 𝜆 2,𝑘 + 2ℎ𝑥 1,𝑘 -ℎ𝑥 1,𝑘 ≥ 0 ⇐⇒ 0 ≤ 𝜆 1,𝑘+1 𝜆 2,𝑘+1 ⊥ 𝜆 1,𝑘+1 -2𝜆 2,𝑘+1 4𝜆 2,𝑘+1 -2𝜆 1,𝑘+1 + 𝑞 1,𝑘 𝑞 2,𝑘 ≥ 0.
From Corollary 4.1, there exists a solution to (6.15). We can compute 𝜆 𝑘+1 in each mode of the LCP where 𝜆 𝑘+1 ≠ 0 :

1. If 𝜆 1,𝑘+1 = 0 and 𝜆 2,𝑘+1 ≥ 0 then ∥𝜆 2,𝑘+1 ∥ ≤ 1 2 ∥𝑞 1,𝑘 ∥.

2. If 𝜆 2,𝑘+1 = 0 and 𝜆 1,𝑘+1 ≥ 0 then ∥𝜆 1,𝑘+1 ∥ ≤ 1 2 ∥𝑞 2,𝑘 ∥ .

3. If 𝜆 1,𝑘+1 > 0 and 𝜆 2,𝑘+1 > 0 there is infinitely many solutions. Indeed, we obtain : (6.16)

𝜆 1,𝑘+1 = 2𝜆 2,𝑘+1 -𝑞 1,𝑘 0 = 2𝑞 1,𝑘 + 𝑞 2,𝑘
From [16, Proposition 6.1.7], there exists a minimal solution 𝜆 min,𝑘+1 . Then, (6.16)'s particular solution 𝜆 2,𝑘+1 = -𝜆 1,𝑘+1 leads to ∥𝜆 min,𝑘+1 ∥ ≤ (1/3) ∥𝑞 𝑘+1 ∥. It follows that 𝜂 = 1 2 < 1 for all ℎ > 0, and 𝜂ℎ 3 ≤ 1 for all 0 < ℎ ≤ 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 3 [START_REF] Acary | An introduction to Siconos[END_REF]. In Figures 7a and7b the simulated variables (𝑥 1 , 𝑥 2 , 𝑧) and (𝜆 1 , 𝜆 2 ) are displayed for time step ℎ = 5•10 -6 s. The initial conditions are (𝑥 1,0 = 5, 𝜆 1,0 = 5, 𝜆 2,0 = 0) and the discrete variables 𝑥 2,𝑘 and 𝑧 𝑘 are given by the discretization of (5.1) (c): (6.17) 7a that the variable 𝑧 possesses a discontinuity. In Figure 8a we plot the error provided in (6.6), with respect to a reference solution taken for ℎ = 5.10 -6 We confirm the order one convergence rate for ℎ > 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 (𝑥 1 , 𝑥 2 ).

       𝑥 2,𝑘 = -𝐵 3 𝜆 𝑘 𝑧 𝑘 = -𝐵 3 𝜆 𝑘+1 -𝜆 𝑘 ℎ , 3

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 [

∥ 𝑃 +𝑃 ⊤ 2 -1 ∥ 2,2
then 𝐷 ≻ 0.

Consider the term ℎ 1 in Assumptions 1 and 4. If ℎ < 1 𝜎 max (𝐴 1 ) then 𝐼 𝑛 1 + ℎ𝐴 1 ≻ 0 from Corollary a.2. If in addition 𝐴 1 ≼ (≺) 0 then 0 ≺ 𝐼 𝑛 1 + ℎ𝐴 1 ≼ (≺) 𝐼 𝑛 1 . Consider now Ā1 = (𝐼 𝑛 1 + ℎ𝐴 1 ) ⊤ (𝐼 𝑛 1 + ℎ𝐴 1 ) = 𝐼 𝑛 1 +ℎ(𝐴 1 + 𝐴 ⊤ 1 ) +ℎ 2 𝐴 ⊤ 1 𝐴 1 = 𝐼 𝑛 1 +ℎ Ã1 , with Ã1 = 𝐴 1 + 𝐴 ⊤ 1 +ℎ𝐴 ⊤ 1 𝐴 1 = Ã⊤ 1 . If 𝐴 1 ≺ 0 (⇒ 𝐴 1 + 𝐴 ⊤ 1 ≺ 0), then Ã1 ≺ 0 for ℎ > 0 small enough, still using Corollary a.2. Thus Ā1 ≺ 𝐼 𝑛 1 for ℎ > 0 small enough. Morover for ℎ > 0 small enough it follows that Ā1 ≻ 0 by Corollary a.2. It is deduced that there exists ℎ max > 0 such that 0 ≺ Ā1 ≺ 𝐼 𝑛 1 for all ℎ ∈ (0, ℎ max ). From [START_REF] Bernstein | Matrix Mathematics. Theory, Facts and Formulas[END_REF]Lemma 8.4.1], 0 < 𝜆 min ( Ā1 ) and 𝜆 max ( Ā1 ) < 1. It is inferred that 𝜎 max (𝐼 𝑛 1 + ℎ𝐴 1 ) = ∥𝐼 𝑛 1 + ℎ𝐴 1 ∥ 2,2 = ℎ 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 [START_REF] Brogliato | Singular mass matrix and redundant constraints in unilaterally constrained Lagrangian and Hamiltonian systems[END_REF][START_REF] Goeleven | Complementarity and Variational Inequalities in Electronics, Mathematical Analysis and its Applications[END_REF] for illustrating examples), [START_REF] Rockafellar | Convex Analysis[END_REF][START_REF] Goeleven | Complementarity and Variational Inequalities in Electronics, Mathematical Analysis and its Applications[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF]. Let 𝑓 : ℝ 𝑛 → ℝ ∪ {+∞} be a proper convex and lower semi-continuous function, we denote by dom(𝑓 ) where M ∈ ℝ 𝑛×𝑛 is a real matrix, q ∈ ℝ 𝑛 a vector and 𝜑 : ℝ 𝑛 → ℝ ∪ {+∞} a proper convex and lower semicontinuous function. The VI in (b.2) is equivalent to the inclusion M𝑢 + q ∈ -𝜕𝜑 (𝑢) ⇔ 𝑢 ∈ (M + 𝜕𝜑) -1 (-q).

The problem in (b.2) is denoted as 𝑉 𝐼 (M, q, 𝜑) in the next proposition. We also set:

(b.3) K (M, 𝜑) = {𝑥 ∈ ℝ 𝑛 |M𝑥 ∈ (dom(𝜑 ∞ )) ★ }.

Note that (dom(𝜑 ∞ )) ★ is the dual cone of the domain of the recession function 𝜑 ∞ while (dom(𝜑)) ∞ (that we may denote also as dom(𝜑) ∞ ) is the recession cone of dom(𝜑). Proposition b.2 ( [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequaolities via semi-complementarity problems: applications in electronics[END_REF]). Let 𝜑 : ℝ 𝑛 → ℝ ∪ {+∞} be a proper, convex and lower semicontinuous function with closed domain, M ∈ ℝ 𝑛×𝑛 , and suppose that M ≽ 0 (not necessarily symmetric). a) If (dom(𝜑)) ∞ ∩ ker{M + M ⊤ } ∩ K (M, 𝜑) = {0} then for each q ∈ ℝ 𝑛 , problem 𝑉 𝐼 (M, q, 𝜑) has at least one solution. e) If M = M ⊤ and 𝜑 (𝑥 + 𝑧) = 𝜑 (𝑥) for all 𝑥 ∈ dom(𝜑) and 𝑧 ∈ ker{M} and ⟨q, 𝑒⟩ ≠ 0 for all 𝑒 ∈ ker{M}, 𝑒 ≠ 0, then problem 𝑉 𝐼 (M, q, 𝜑) has at most one solution.

f) If M = M ⊤ , then 𝑢 is a solution of 𝑉 𝐼 (M, q, 𝜑) if and only if it is a solution of the optimization problem min 𝑥 ∈ℝ 𝑛 1 2 𝑥 ⊤ M𝑥 + ⟨q, 𝑥⟩ + 𝜑 (𝑥).

Notice that the function 𝜑 (•) will never be strictly convex in our case (it is an indicator function) so that the strict convexity argument of [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequaolities via semi-complementarity problems: applications in electronics[END_REF]Theorem 5] which applies when M is a 𝑃 0 -matrix never holds. The study of VIs as in (b.2) can be traced back to [START_REF] Stampacchia | Formes bilinéaires coercitives sur les ensembles convexes[END_REF].

appendix c some convex analysis and complementarity theory tools

If 𝐾 ⊂ ℝ 𝑛 is a set, then 𝐾 ★ = {𝑧 ∈ ℝ 𝑛 |⟨𝑧, 𝑥⟩ ≥ 0 for all 𝑥 ∈ 𝐾 } is its dual cone. Its closure is denoted K. Let 𝐾 be a nonempty closed convex cone, then: (𝑧 -𝑦) ⊤ 𝑀 (𝑧 -𝑦).

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 𝑉 𝐼 (M, q, 𝜑). The first equivalence is [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 6.46].

• Factor of ∥𝑥 1,1 ∥ (all terms with 𝑗 = 0): ℎ 𝑘 -1 1 +ℎ 𝑘 -2 1 (𝜂ℎ 3 ) +ℎ 𝑘 -3 1 (𝜂ℎ 3 ) 2 +. . .+ℎ 1 (𝜂ℎ 3 ) 𝑘 -2 + (𝜂ℎ 3 ) 𝑘 -1 = 𝑘 𝑖=1 ℎ 𝑘 -𝑖 1 (𝜂ℎ 3 ) 𝑖 -1 .

• Factor of ∥𝑥 1,2 ∥ (all terms with 𝑗 = 1): ℎ 𝑘 -2 1 +ℎ 𝑘 -3 1 (𝜂ℎ 3 ) +ℎ 𝑘 -4 1 (𝜂ℎ 3 ) 2 +. . .+ℎ 1 (𝜂ℎ 3 ) 𝑘 -3 + (𝜂ℎ 3 ) 𝑘 -2 = 𝑘 𝑖=2 ℎ 𝑘 -𝑖 1 (𝜂ℎ 3 ) 𝑖 -2 .

• . . .

• Factor of ∥𝑥 1,𝑘 -1 ∥ (all terms with 𝑗 = 𝑘 -2): ℎ 1 + (𝜂ℎ 3 ) = 𝑘 𝑖=𝑘 -1 ℎ 𝑘 -𝑖 1 (𝜂ℎ 3 ) 𝑖 -𝑘+1 .

• Factor of ∥𝑥 1,𝑘 ∥ (all terms with 𝑗 = 𝑘): 1 = 𝑘 𝑖=𝑘 ℎ 𝑘 -𝑖 1 (𝜂ℎ 3 ) 𝑖 -𝑘 .

Therefore one infers (5.22).

appendix g minimality of state space realizations • rank(𝐴 -𝑠𝑃 𝐵) = 𝑛 for all 𝑠 ∈ ℂ,

• rank(𝑃 𝐵) = 𝑛,

• rank(𝐴 ⊤ -𝑠𝑃 ⊤ 𝐶 ⊤ ) = 𝑛,

• rank(𝑃 ⊤ 𝐶 ⊤ ) = 𝑛,

• 𝐴 Ker(𝑃) ⊆ Im(𝑃). 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 𝜆 𝑘 and 𝑧 𝑘 are used in (h.1) (a) and (h.2) instead of 𝜆 𝑘+1 and 𝑧 𝑘+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.

  second set in S is obvious. The third set necessitates the calculation of (dom(𝜑 ∞ )) ★ , 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). □ Once the existence of 𝜆 𝑘+1 is established for any 𝑘 ≥ 0, then the existence of 𝑥 1,𝑘+1 , 𝑥 2,𝑘+1 and 𝑧 𝑘+1 can be inferred from (3.1) for any ℎ > 0. Remark 3.3. If 𝐾 = ℝ 𝑚 + then other results like [43, Theorem 3.8.6] may be used to study the LCP (3.2), using for instance the copositivity of Dℎ .

Proposition 5 . 1 .

 51 1) (b) to cone complementarity conditions 𝐾 ★ ∋ 𝜆(𝑡) ⊥ 𝑤 (𝑡) ∈ 𝐾, 𝐾 ⊆ ℝ 𝑚 a nonempty closed convex cone, similarly to what is done in Proposition 3.2. Thus Φ = ℝ 𝑛 1 ×𝐾 ★ . Let ℎ > 0 and 𝑘 ≥ 0 be given. Let us consider the DI in (5.3) and its discretization in (5.5), with Φ = ℝ 𝑛 1 × 𝐾 ★ , with 𝐾 a nonempty closed convex cone. Let us denote 𝑞 𝑘 Δ = -P ℎ 𝜉 𝑘 + ℎ𝐺 𝑘 . Let us consider the set:

  3 and 5.4 to get a new upper bound. Let us denote ℎ 1 Δ

7 10 - 4 H 4 ,

 44 , 𝑅 = 1Ω, 𝐶 = 10 -1 F. Let us define 𝑆 = 𝐶𝑇 = (-0.707 -1 0), 𝑆

Figure 2 :

 2 Figure 2: Implicit Euler numerical simulation of (6.3) with ℎ = 5 • 10 -6 s. The numerical simulations of the variables 𝑥 1 (𝑡), 𝑥 2 (𝑡) and 𝑥 3 (𝑡) = 𝑧 (𝑡) are depicted.

Figure 3 :

 3 Figure 3: Implicit Euler numerical simulation of (6.3) with ℎ = 5 • 10 -6 s. The numerical simulations of the variables 𝜆(𝑡) and 𝑤 (𝑡) are depicted.

Figure 4 :

 4 Figure 4: Numerical convergence of the implicit Euler scheme applied to (6.3). The reference curve is taken for ℎ = 5.10 -6 s. The error formula used is given in (6.6).

Figure 5 :

 5 Figure 5: 𝑅𝐿𝐶𝐷 passive circuit.

( a )

 a Variables 𝑥 2 (𝑡) and 𝑥 3 (𝑡) = 𝑧 (𝑡). (b) Variables 𝜆 1 (𝑡) and 𝑤 1 (𝑡).

Figure 6 :

 6 Figure 6: Implicit Euler numerical simulation of (6.9) with ℎ = 5 • 10 -6 s. The solutions 𝑥 1 (𝑡) are not displayed as they are independent of 𝑧 (𝑡), 𝜆 1 (𝑡). The variable 𝜆 2 (𝑡) trivially equals 0.

  ( D + DT ) ≻ 0. It follows from [28, Theorem 5.4 (A'), Theorem 4.6] that the quadruple (𝐴 1 , 𝐵 1 , 𝐶 1 , D) is strictly passive.

  (a) Variables 𝑥 1 (𝑡), 𝑥 2 (𝑡) and 𝑧 (𝑡). (b) Variables 𝜆 1 (𝑡) and 𝜆 2 (𝑡).

Figure 7 :

 7 Figure 7: Implicit Euler numerical simulation of (6.13) with ℎ = 5 • 10 -6 s

Figure 8 :

 8 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 ℎ = 5 • 10 -6 s.

Δ=

  {𝑥 ∈ ℝ 𝑛 | 𝑓 (𝑥) < +∞} the domain of the function 𝑓 (•). The Fenchel transform 𝑓 ★ (•) of 𝑓 (•) is the proper, convex and lower semicontinuous function defined by (b.1)(for all 𝑧 ∈ ℝ 𝑛 )| 𝑓 ★ (𝑧) = sup 𝑥 ∈ dom(𝑓 ) {⟨𝑥, 𝑧⟩ -𝑓 (𝑥)}.The subdifferential 𝜕𝑓 (𝑥) of 𝑓 (•) at 𝑥 ∈ ℝ 𝑛 is defined by𝜕𝑓 (𝑥) = {𝜔 ∈ ℝ 𝑛 | 𝑓 (𝑣) -𝑓 (𝑥) ≥ ⟨𝜔, 𝑣 -𝑥⟩, ∀𝑣 ∈ ℝ 𝑛 }.Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

b)

  Suppose that (dom(𝜑)) ∞ ∩ ker{M + M ⊤ } ∩ K (M, 𝜑) ≠ {0}. If there exists 𝑥 0 ∈ dom(𝜑) such that (b.4) ⟨q -M ⊤ 𝑥 0 , 𝑣⟩ + 𝜑 ∞ (𝑣) > 0, ∀𝑣 ∈ (dom(𝜑)) ∞ ∩ ker{M + M ⊤ } ∩ K (M, 𝜑), 𝑣 ≠ 0, then problem 𝑉 𝐼 (M, q, 𝜑) has at least one solution. b') If M = M ⊤ then one can take 𝑥 0 = 0 in b). c) If 𝑢 1 and 𝑢 2 denote two solutions of problem 𝑉 𝐼 (M, q, 𝜑) then 𝑢 1 -𝑢 2 ∈ ker{M + M ⊤ }.d) If M = M ⊤ and 𝑢 1 and 𝑢 2 denote two solutions of problem 𝑉 𝐼 (M, q, 𝜑), then ⟨q, 𝑢 1 -𝑢 2 ⟩ = 𝜑 (𝑢 2 ) -𝜑 (𝑢 1 ).

(c. 1 )

 1 𝐾 ★ ∋ 𝑥 ⊥ 𝑦 ∈ 𝐾 ⇐⇒ 𝑥 ∈ -N 𝐾 (𝑦) ⇐⇒ 𝑦 ∈ -N 𝐾 ★ (𝑥).Let 𝑀 = 𝑀 ⊤ ≻ 0, 𝑥 and 𝑦 two vectors, then(c.2) 𝑀 (𝑥 -𝑦) ∈ -N 𝐾 (𝑥) ⇔ 𝑥 = proj 𝑀 [𝐾; 𝑦] ⇔ 𝑥 = min 𝑧 ∈𝐾 1 2

Let us report the results stated in [ 28 ,

 28 Theorem 3.1] [44, Theorems 2-2.1, 2-3.1, 2-6.3]. Theorem g.1. Let 𝐻 (𝑠) = 𝐶 (𝑠𝑃 -𝐴) -1 𝐵 + 𝐷, 𝑠 ∈ ℂ, be a rational transfer function where 𝑃 and 𝐴 are 𝑛 × 𝑛 matrices. Then 𝑃 𝑥 (𝑡) = 𝐴𝑥 (𝑡) + 𝐵𝑢 (𝑡), 𝑦 (𝑡) = 𝐶𝑥 (𝑡) + 𝐷𝑢 (𝑡) is a minimal realization of 𝐻 (𝑠) if, and only if, the following conditions are satisfied:

appendix h differences between ( 3 . 1 )

 31 and in(5.5) Let us consider the discrete scheme(5.5) with 𝐸 (𝑡) = 0 and 𝑋 1 full rank, this is rewritten as:𝑥 1,𝑘+1 = 𝑥 1,𝑘 + ℎ𝐴 1 𝑥 1,𝑘 + ℎ𝐵 1 𝜆 𝑘 (h.1a) 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 𝑘+1 = 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 𝑘 -ℎ𝐶 1 𝑥 1,𝑘ℎ D𝜆 𝑘 + 𝑤 𝑘+1 (h.1b) 0 ≤ 𝜆 𝑘+1 ⊥ 𝑤 𝑘+1 = 𝐶 1 𝑥 1,𝑘+1 + 𝐶 2 𝑥 2,𝑘+1 + 𝐶 3 𝑧 𝑘+1 + 𝐷𝜆 𝑘+1 ≥ 0 (h.1c) 𝑥 2,𝑘 = -𝐵 3 𝜆 𝑘 (h.1d) 𝑧 𝑘 = -𝐵 2 𝜆 𝑘 -𝐵 3 𝜆 𝑘+1 -𝜆 𝑘 ℎ . (h.1e)This is the numerical scheme with which the examples have been simulated. Combining (h.1) (d) and (e), gives: (h.2) 𝑥 2,𝑘+1 = 𝑥 2,𝑘 + ℎ𝑧 𝑘 + ℎ𝐵 2 𝜆 𝑘 .

  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 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) such that ∥𝜆 min 𝑘+1 ∥ ≤ ∥𝜆 𝑘+1 ∥ for any other solution 𝜆 𝑘+1 . Moreover there exists a constant 𝜂 > 0 depending only on Dℎ

2) and (3.1), i.e., 𝐾 = 𝐾 ★ = ℝ 𝑚 + , see Remark 4.2. The next result characterizes least-norm solutions of the LCP. Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . Corollary 4.1.

  2,𝑘 = -𝐵 3 𝜆 min 𝑘 , and ℎ𝑧 𝑘+1 = -𝐵 3 𝜆 min 𝑘+1 -𝑥 2,𝑘 -ℎ𝐵 2 𝜆 min

	𝑘+1 , thus it is inferred that
	𝑥 2,𝑘+1 = -𝐵 3 (𝜆 min 𝑘+1 + 𝜆 min 𝑘 ) -𝑥 2,𝑘 . Consequently (4.1) is rewritten equivalently as:
	(4.2)

  1,𝑘 , 𝑥 2,𝑘 ) -𝑥 2,𝑘 -ℎ𝐵 2 𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) ∥𝜆 min 𝑘+1 (𝑥 1,𝑘 , 𝑥 2,𝑘 ) ∥ ≤ 𝜂 ∥ℎ𝛽 𝑘 ∥, ℎ𝛽 𝑘 = ℎ𝐶 1 (𝐼 𝑛 1 + ℎ𝐴 1 )𝑥 1,𝑘 -𝐶 3 𝑥 2,𝑘 , 𝜂 > 0, for all 𝑘 ∈ {1, 𝑁 -1}, given initial data 𝑥 1,0 and 𝑥 2,0 . Starting from (4.2) the next result is obtained: 𝛼 𝑘+1 ≤ (𝑀 𝑖 𝑀 + 𝑀 𝑖 -1 𝑈 )𝛼 𝑘 -𝑖 + 𝑀 𝑖 𝑈 𝛼 𝑘 -𝑖 -1

	Lemma 4.3. Let 𝛼 1,𝑘	Δ = ∥𝑥 1,𝑘 ∥, 𝛼 2,𝑘	Δ = ∥𝑥 2,𝑘 ∥, 𝛼 𝑘 = (𝛼 1,𝑘 𝛼 2,𝑘 ) ⊤ , then:
	(4.3)
	Brogliato, Rocca		Analysis of the implicit Euler time-discretization of passive . . .

that the boundedness of the sequence {𝛼 𝑘 } 𝑘 ≥2 depends on the boundedness of 𝑀 𝑘 , 𝑘 ≥ 1. Apart from 𝑑 and 𝑓 , all the other terms in 𝑀 and 𝑈 are proportional to ℎ, 𝜂ℎ or to ℎ 1 Δ

  1,𝑘 -1 ∥). 𝑘 -1. Let it hold for 𝑖 > 0, we have 𝛼 𝑘 -𝑖 ≤ 𝑀𝛼 𝑘 -𝑖 -1 + 𝑈 𝛼 𝑘 -𝑖 -2 . Thus 𝛼 𝑘+1 ≤ (𝑀 𝑖 𝑀 + 𝑀 𝑖 -1 𝑈 )𝑀𝛼 𝑘 -𝑖 -1 + (𝑀 𝑖 𝑀 + 𝑀 𝑖 -1 𝑈 )𝑈 𝛼 𝑘 -𝑖 -2 + 𝑀 𝑖 𝑈 𝛼 𝑘 -𝑖 -1 = (𝑀 𝑖+1 𝑀 + 𝑀 𝑖 𝑈 )𝛼 𝑘 -𝑖 -1 + 𝑀 𝑖+1 𝑈 𝛼 𝑘 -𝑖 -2 .

	(4.5c)
	One infers that:
	(4.6)	𝛼 𝑘+1 ≤ 𝑀𝛼 𝑘 + 𝑈 𝛼 𝑘 -1
	The upper bound in (4.3) therefore holds for 𝑖 = 0. Let us prove by induction that it holds for all
	0 ≤ 𝑖 ≤ The
	result is proved.	□
	Our goal is now to examine (4.4). Let us come back to the matrices 𝑀 and 𝑈 in Lemma 4.3. It is clear
	from (4.4)

  1 , D) is minimal, then strict state passivity and SPRness are equivalent, see Appendix d. Lemma 4.9. Assume that: 1) the quadruple (𝐴 1 , 𝐵 1 , 𝐶 1 , D) is strictly state passive with D + D⊤ = 0 and dissipation

  1,0 ∥𝛼 𝑘 2 with 𝛼 < 1. Fom (4.5) (a) and Assumption 3 it follows that ∥𝜆 min 𝑘+1 ∥ → 0 as 𝑘 → +∞, and {𝜆 min 𝑘 } 𝑘 is uniformly bounded. From (4.1) it is inferred that 𝑥 2,𝑘 → 0 as 𝑘 → +∞ and {𝑥 2,𝑘 } 𝑘 is uniformly bounded. □

  ℎ𝐺 𝜆 𝑘 for all 𝑘 ≥ 0. The second formulation shows the iterative feature of this GE for 𝜆 𝑘 , which makes it quite different from the LCCP (3.3). It is equivalent to the LCP(𝐵 ⊤ 3 𝑋 3 𝐵 3 , 𝑞 𝜆 𝑘 (𝜆 𝑘 , 𝑥 1,𝑘 , 𝐺 𝜆 𝑘 )). 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 𝜆 min 𝑘+1 , 𝑘 ≥ 0, satisfiying ∥𝜆 min 𝑘+1 (𝜆 𝑘 , 𝑥 1,𝑘 , 𝐺 𝜆 𝑘 ) ∥ ≤ 𝜂 ∥𝑞 𝜆 𝑘 (𝜆 𝑘 , 𝑥 1,𝑘 , 𝐺 𝜆 𝑘 ) ∥, 𝑘 ≥ 0, where 𝜂 > 0 depends only on 𝐵 ⊤ 3 𝑋 3 𝐵 3 . Therefore ∥𝜆 min 𝑘+1 (𝜆 𝑘 , 𝑥 1,𝑘 , 𝐺 𝜆 𝑘

  ∥, and thus being made of powers of ℎ 1 and of ℎ 3 𝜂.

	Calculations show that (see Appendix f)
	(5.22)

  1,𝑘 -10 -4 ℎ𝑢 2,𝑘 , Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . and 𝑥 1,𝑘+1 , 𝑥 2,𝑘+1 , 𝑧 𝑘+1 as in (h.1). Let us check Assumptions 1 and 4 c

  1 𝑘+1 and 𝜆 2 𝑘+1 satisfy 𝜆 1 2,𝑘+1 = 𝜆 2 2,𝑘+1 for all ℎ ≥ 0, in view of the form of Dℎ . Then both 𝑥 1,𝑘+1 and 𝑥 2,𝑘+1 are unique due to the form of 𝐵 1 and 𝐵 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 ℎ 3 = 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 𝐾 = ℝ + × ℝ ≠ ℝ 2 + . The following problem Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . . 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) 𝐾 𝐵 ⊤ 3 𝑋 3 𝐵 3 𝜆 𝑘+1 + 𝑞 𝜆 𝑘 ⊥ 𝜆 𝑘+1 ∈ 𝐾 ★ ,

	√ 2ℎ 2 + ℎ 𝑅 -1 𝑅 -1 ℎ 2 + ℎ √ 2 , and 0 ≤ ℎ 3 ≤ 1 for all ℎ ∈ [0, 2𝑅 2 , using the Hölder induced norm ∥ • ∥ 2,2 = 𝜎 max (•) 𝑘 is defined in (5.14)), i.e., one finds a smaller ℎ 3 = 1+𝑅 2 ]. (where 𝑞 𝜆 (6.12)

  ≻ 0, it follows that the multipliers 𝜆 are Lipschitz continuous functions of the state (𝑥 1 , 𝑥 2 , 𝑧). Thus the LDCS is a with Lipschitz continuous nonlinear and nonsmooth right-hand side. From (6.13) it follows that: 𝐴 1

	that since 𝐷 =	2 -1 0 1

-1

𝐵 is positive real

[START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF] Definition 2.34]

, and the state-space realization is minimal (see Appendix g). In fact 𝐻 (𝑠) is even strongly (or extended) positive real

[START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF] Definition 2.78]

. Thus

[START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF] Theorem 5.4

(A') (B')] applies. It is noteworthy

Brogliato, Rocca

Analysis of the implicit Euler time-discretization of passive . . .

  37, Theorem 2.11], and a corollary of it. Let us recall that for a given 𝑀 ∈ ℝ 𝑛×𝑛 , ∥𝑀 ∥ 2,2 is the induced matricial norm such that ∥𝑀 ∥ 2,2 = 𝜎 max (𝑀) (the largest singular value). Theorem a.1 ([START_REF] Chen | Perturbation bounds of P-matrix linear complementarity problems[END_REF]). Let 𝑀 ∈ ℝ 𝑛×𝑛 be a positive definite matrix. Then every matrix𝐴 ∈ {𝐴 ∈ ℝ 𝑛×𝑛 |

	𝑀 + 𝑀 ⊤	-1
	2
		2,2

∥𝑀 -𝐴∥ 2,2 < 1} is positive definite. Corollary a.2 ([21]). Let 𝐷 = 𝑃 + 𝑁 , where 𝐷, 𝑃 and 𝑁 are 𝑛 × 𝑛 real matrices, and 𝑃 ≻ 0, not necessarily symmetric. If (a.1) ∥𝑁 ∥ 2,2 < 1

Brogliato, RoccaAnalysis of the implicit Euler time-discretization of passive . . .

Source code of the examples can be found at https://github.com/siconos/siconos-tutorials/tree/master/.sandbox/mlcspaper Brogliato, Rocca Analysis of the implicit Euler time-discretization of passive . . .

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.

We denote by Dom(𝜕𝑓 ) Δ = {𝑥 ∈ ℝ 𝑛 |𝜕𝑓 (𝑥) ≠ ∅} the domain of the subdifferential operator 𝜕𝑓 : ℝ 𝑛 → ℝ 𝑛 . Recall that (see, e.g., Theorem 2, Chapter 10, Section 3 in [START_REF] Aubin | Applied Functional Analysis[END_REF]): Dom(𝜕𝑓 ) ⊂ dom(𝑓 ).

Let 𝑥 0 be any element in the domain dom(𝑓 ) of 𝑓 (•), the recession function 𝑓 ∞ (•) of 𝑓 (•) is defined by (for all 𝑥 ∈ ℝ 𝑛 ) : 𝑓 ∞ (𝑥) = lim 𝜆→+∞ 1 𝜆 𝑓 (𝑥 0 + 𝜆𝑥).

The function 𝑓 ∞ : ℝ 𝑛 → ℝ ∪ {+∞} is a proper convex and lower semicontinuous function which describes the asymptotic behavior of 𝑓 (•).

Let 𝐾 ⊂ ℝ 𝑛 be a nonempty closed convex set. Let 𝑥 0 be any element in 𝐾. The recession cone of 𝐾 is defined by [START_REF] Rockafellar | Convex Analysis[END_REF] [89, Definition 1.11]:

The set 𝐾 ∞ is a nonempty closed convex cone that is described in terms of the directions which recede from 𝐾. The indicator function of a set

If 𝐾 is closed non empty convex, we have 𝜕Ψ 𝐾 (𝑥) = N 𝐾 (𝑥), the so-called normal cone to 𝐾 at 𝑥, defined as N 𝐾 (𝑥) = {𝑣 ∈ ℝ 𝑛 | 𝑣 ⊤ (𝑠 -𝑥) ≤ 0 for all 𝑠 ∈ 𝐾 }. When 𝐾 is finitely represented, i.e., 𝐾 = {𝑥 ∈ ℝ 𝑛 | 𝑘 𝑖 (𝑥) ≥ 0, 1 ≤ 𝑖 ≤ 𝑚}, and if the functions 𝑘 𝑖 (•) satisfy some constraint qualification (like, independency, or extensions like the MFCQ), then N 𝐾 (𝑥) is generated by the outwards normals at the active constraints 𝑘 𝑖 (𝑥) = 0, i.e., N 𝐾 (𝑥)

Let us here recall some important properties of the recession function and recession cone (see e.g., [START_REF] Bertsekas | Convex Optimization Theory[END_REF]Proposition 1.4.8]): Proposition b.1. The following statements hold: a) Let 𝑓 1 : ℝ 𝑛 → ℝ ∪ {+∞} and 𝑓 2 : ℝ 𝑛 → ℝ ∪ {+∞} be two proper, convex and lower semicontinuous functions. Suppose that 𝑓 1 + 𝑓 2 is proper. Then for all 𝑥 ∈ ℝ 𝑛 :

b) Let 𝑓 : ℝ 𝑛 → ℝ ∪ {+∞} be a proper, convex and lower semicontinuous function and let 𝐾 be a nonempty closed convex set, such that 𝑓 + Ψ 𝐾 is proper (equivalently dom(𝑓 ) ∩ 𝐾 is non empty). Then for all 𝑥 ∈ ℝ 𝑛 :

c) Let 𝐾 ⊂ ℝ 𝑛 be a nonempty, closed and convex set. Then for all 𝑥 ∈ ℝ 𝑛 :

appendix d passive descriptor variable systems

Dissipative systems have been formalized rigorously by J.C. Willems in his two seminal articles [START_REF] Willems | Dissipative dynamical systems, part I: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems, part II: Linear systems with quadratic supply rates[END_REF].

Passive systems are a particular case of dissipative systems. A quadruple (𝐴, 𝐵, 𝐶, 𝐷) is said passive if the next LMI is satisfied:

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 directy related to strongly SPR transfer matricess [25, p.62 and section 3.12.2], while SPR transfer matrices with minimal realisations are strictly state passive [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF]Theorem 4.73] [66], i.e., -𝑋𝐴 -𝐴 ⊤ 𝑋 + 𝜇𝑋 ≽ 0 for some 𝜇 > 0, and one may have 𝐷 + 𝐷 ⊤ ≽ 0 but not ≻ 0.

Positive real transfer matrices [25, Definition 2.34] with a state space realization represented by quintuples (𝑃, 𝐴, 𝐵, 𝐶, 𝐷) which are minimal (see Appendix g for conditions of minimality of a statespace realization for descriptor systems) can be rewritten as in (2.2) [START_REF] Camlibel | Extension of Kalman-Yakubovich-Popov Lemma to descriptor systems[END_REF][START_REF] Freund | An extension of the positive real lemma to descriptor systems[END_REF], and they have the positive semidefinite storage function

, with supply rate 𝜆 ⊤ 𝑤. Indeed, by differentiating 𝑉 (𝑥 1 (𝑡), 𝑥 2 (𝑡)) along the trajectories of (2.2), adding ±𝜆 ⊤ 𝑤 and completing the squares, one can compute that:

where Let us recall some definitions, where 𝐴 ★ denotes the conjugate transpose of 𝐴 (see [25, Chapter 2] for details): Definition d.1. A transfer matrix 𝐻 (𝑠) ∈ ℂ 𝑚×𝑚 , 𝑠 ∈ ℂ, is positive real (PR) if: 1) 𝐻 (𝑠) has no poles in Re(𝐻 (𝑠)) > 0, 2) 𝐻 (𝑠) is real for positive real 𝑠, 3) 𝐻 (𝑠) + 𝐻 ★ (𝑠) ≽ 0 for all Re(𝑠) > 0. Let 𝐻 (𝑠) be regular, then it is strictly positive real (SPR) if 𝐻 (𝑠 -𝜖) is PR for some 𝜖 > 0. It is strongly SPR (SSPR) if 1) 𝐻 (𝑠) is analytic in Re(𝑠) ≥ 0, 2) Re(𝐻 ( 𝑗𝜔)) ≥ 𝛿 > 0 for all 𝜔 ∈ [-∞, +∞] and some real 𝛿. appendix e gronwall-bachpatte inequality [START_REF] Bachpatte | On the discrete generalisations of Gronwall's inequality[END_REF][START_REF] Dragomir | Some Gronwall Type Inequalities and Applications[END_REF] Let us consider a positive real numbers sequence {𝑥 𝑘 } such that 𝑥 𝑘 ≤ 𝛼 𝑘 + 𝑘 -