Tomáš Ligurský
email: tomas.ligursky@upol.cz

Yves Renard
email: yves.renard@insa-lyon.fr

Ku

G Ku

N Ku

Z Ku

Νzl Ku

Τcl Ku

Τzl Ku

τzl KG,ci KN,u KN,G KN,N KN,Z KN,νzl KN,τcl KN,τzl KN,ci KZ Kg

G Kg

N Kg

Z Kg

G Kz

Z Kz,N Kz

Kz

Ci Kz

Kνzl

G Kνzl

N Kνzl

Z Kνzl

Νzl Kνzl

Τcl Kνzl

Τzl Kνzl

Ci Kνzl

Kτcl

G Kτcl

N Kτcl

Z Kτcl

Νzl Kτcl

Τcl Kτcl

Τzl Kτcl

Ci Kτcl

Kτzl

G Kτzl

N Kτzl

Z Kτzl

Νzl Kτzl

Τcl Kτzl

Τzl Kτzl

Ci Kτzl

Kci

G Kci

N Kci

Z Kci

Νzl Kci

A Continuation Problem for Computing Solutions of Discretised Evolution Problems. Part II -Application to Quasi-Static Contact Problems with Friction

Keywords: contact, Coulomb friction, incremental problem, continuation method, firstorder system

A continuation problem for finding successive solutions of discretised plane quasistatic contact problems with friction is proposed and studied. Various formulations of the first-order system are derived and results of existence and uniqueness of its solutions are obtained. Conditions guaranteeing continuation of a solution curve along a direction from the first-order system are given. Numerical continuation of solution curves is performed on examples with non-smooth folds.

Introduction

In the companion paper [16], we have proposed an abstract continuation problem and we have laid foundations for continuation of its solution curves provided that the problem is piecewise C 1 (P C 1 ). The present paper concerns a special case of this problem, namely, a continuation problem for discretised quasi-static plane contact problems with Coulomb friction in large deformations. As far as we know, the only papers dealing with non-smooth continuation of fully discretised contact problems with friction are [11,8,7], but those are all restricted to elastostatics and the small-deformation framework.

After formulating a continuous quasi-static contact problem in large deformations in Section 2, we discretise it in Section 3. The spatial discretisation is done with a finiteelement method applied to a mixed variational formulation of the problem while the time derivative is approximated by the backward difference. Following [2], we write the corresponding continuation problem so that it fits exactly the general framework.

Section 4 deals with the first aim of this paper, which is to employ a specific structure of the present problem to do more precise analysis of its first-order system. Let us recall that the first-order system consists in computing directional derivatives of a curve solving the continuation problem. In particular, we formulate the first-order system in such a way that it is close to a rate problem of a quasi-static contact problem in a sense. In the rate problem, one computes first-order right rates of change, that is, right-hand side time derivatives, of a solution of a quasi-static problem for a given rate of input data when time derivatives are not approximated. This problem can serve for direct study of tangent behaviour of solutions of the quasi-static problem. Making use of the similarity between these two problems, we adapt the analysis of the rate problem from [17] to our first-order system. Moreover, we investigate the result of continuation of solution curves in directions solving the first-order system from the abstract frame [16].

The second objective of this paper is to test the general method of numerical continuation from [16]. A similar algorithm has already been tested on static contact problems in [11,8], but only on finite-element models with very small number of degrees of freedom. Section 5 here presents results for more realistic models.

Notation and Preliminaries

The following notation is used throughout the paper: For a vector x ∈ R N , a matrix A ∈ R M ×N and index sets I ⊂ {1, . . . , M } and J ⊂ {1, . . . , N }, x i stands for the ith component of x, x J is the sub-vector of x composed from the components x i , i ∈ J, and A I,J is the sub-matrix of A with rows and columns specified by I and J, respectively. Furthermore, x • y = x ⊤ y is the scalar product of vectors x and y. The gradients of a real-valued function f and a vector-valued function f at a point x are denoted by ∇f (x) and ∇f (x), respectively, the partial gradients of f and f with respect to y at (x, y) are denoted by ∇ y f (x, y) and ∇ y f (x, y), respectively, and ∇ 2 f (x) is the Hessian of f at x. For reader's convenience, we recall the essentials from theory of P C 1 -functions [20]: Definition 1. (i) A function H : R M → R N is P C 1 if it is continuous and for every ȳ ∈ R M , there exist an open neighbourhood O of ȳ and a finite family of C 1 -functions H (i) : O → R N , i ∈ I, such that ∀y ∈ O : H(y) ∈ {H (i) (y); i ∈ I}.

The functions H (i) are termed selection functions of H at ȳ.

(ii) A selection function H (i) of a P C 1 -function H is called active at ȳ if

H (i) (ȳ) = H(ȳ).
The set of indices of active selection functions of H at ȳ is denoted by I H (ȳ).

Obviously, a composition of two P C 1 -functions is again a P C 1 -function and a vectorvalued function is a P C 1 -function if and only if all its component functions are P C 1 . The following chain rule holds:

Proposition 1. Let G : R M → R N and H : R N → R P be P C 1 -functions. Then the composite function H • G is directionally differentiable at any point and

(H • G) ′ (y; z) = H ′ (G(y); G ′ (y; z)).
If the inner function is a function of a real variable, one can also prove:

Proposition 2. If G : [s 0 , s 0 +δ) → R N is differentiable at s 0 from the right and H : R N → R P is a P C 1 -function, then H • G has a right-hand side derivative at s 0 and

(H • G) ′ (s 0 ) = H ′ (G(s 0 ); G ′ (s 0 ))
(for any function f of a real variable, f ′ denotes its right-hand side derivative for brevity of notation here and in what follows).

Plane Contact Problems with Coulomb Friction

We start our exposition with a brief introduction into quasi-static plane contact problems with Coulomb friction in large deformations. The reader is referred to [4] for a more detailed introduction to large deformations of elastic materials.

Let us consider deformation of an elastic body in a time interval [0, T ]. The reference configuration of the body is represented by the closure of a bounded domain Ω ⊂ R 2 with a Lipschitz-continuous boundary ∂Ω. The boundary ∂Ω is split into three disjoint relatively open subsets Γ D , Γ N and Γ c such that ∂Ω = Γ D ∪ Γ N ∪ Γ c . By u, we denote the displacement of the body and the deformation ϕ of the body is defined by ϕ = id + u, where id is the identity map. At each time t ∈ [0, T ], the deformed body ϕ(t, Ω) is subject to body forces with the density f ϕ while surface forces with the density h ϕ act on ϕ(t, Γ N ). The displacement U D is imposed on Γ D and points from Γ c may come into contact with a fixed curved rigid obstacle represented by a closed set O ⊂ R 2 . It is supposed that there exist g : R 2 → R of the class C 1 and a neighbourhood O of the boundary ∂O such that

g(x) > 0, ∀x ∈ Int O ∩ O, g(x) = 0, ∀x ∈ ∂O, g(x) < 0, ∀x ∈ Ext O ∩ O, ∇g = 0 in O and ϕ(t, Γ c ) ⊂ O for any t ∈ [0, T ]
. This enables us to extend the unit inward normal and tangent vectors to the obstacle to the vicinity O of the obstacle as follows (see Figure 1):

ν(x) = ∇g(x) ∇g(x) , τ (x) = (-ν 2 (x), ν 1 (x)). ϕ(t, Ω) ϕ(t, ΓD) ϕ(t, ΓN ) ϕ(t, ΓN ) ϕ(t, Γc) O n ν τ Figure 1: Geometry of the problem.
Finally, an initial displacement u 0 is given at the time t = 0. The classical formulation of this problem is the following:

Find u : [0, T ] × Ω → R 2 such that u(t, .) ∈ U ad for any t ∈ [0, T ] and div σ(t, x) + f (t, x) = 0, (t, x) ∈ (0, T ) × Ω, σ(t, x) = σ(x, I + ∇u(t, x)), (t, x) ∈ (0, T ) × Ω, u(t, x) = U D (t, x), (t, x) ∈ (0, T ) × Γ D , σ(t, x)n(x) = h(t, x), (t, x) ∈ (0, T ) × Γ N , g(x + u(t, x)) ≤ 0, T ν (t, x) ≤ 0, g(x + u(t, x))T ν (t, x) = 0, (t, x) ∈ (0, T ) × Γ c , |T τ (t, x)| ≤ -F (x)T ν (t, x), (t, x) ∈ (0, T ) × Γ c , T τ (t, x) = F (x)T ν (t, x) uτ (t, x) | uτ (t, x)| if uτ (t, x) = 0, (t, x) ∈ (0, T ) × Γ c , u(0, x) = U 0 (x), x ∈ Ω.                                        (1)
Here, I denotes the identity matrix, ∇u is the partial gradient of u with respect to x and U ad is the set of kinematically admissible displacements:

U ad = {v : Ω → R 2 "smooth enough"; id + v is injective on Ω, det(I + ∇v) > 0 in Ω}.
The first Piola-Kirchhoff stress tensor σ is defined from the Cauchy stress tensor σ ϕ by

σ(t, x) = (det(I + ∇u(t, x)))σ ϕ (t, x + u(t, x))(I + ∇u(t, x)) -⊤ (2)
and σ is a given response function for the first Piola-Kirchhoff stress of the elastic material of the body. By n, we denote the unit outward normal vector along ∂Ω and

f (t, x) = det(I + ∇u(t, x))f ϕ (t, x + u(t, x)), h(t, x) = det(I + ∇u(t, x)) (I + ∇u(t, x)) -⊤ n(x) h ϕ (t, x + u(t, x))
are the densities of volume and surface forces related to the reference configuration. Further,

T ν (t, x) = T (t, x) • ν(x + u(t, x)), T τ (t, x) = T (t, x) • τ (x + u(t, x))
are components of the first Piola-Kirchhoff stress vector T at (t, x) in the directions ν and τ , T (t, x) = σ(t, x)n(x). The friction coefficient F is a non-negative function of x and

uτ (t, x) = ∂u ∂t (t, x) • τ (x + u(t, x))
is the tangent velocity of a point x at a time t.

To interpret the boundary conditions on Γ c , let us consider (t, x) ∈ (0, T ) × Γ c fixed and assume that ϕ is sufficiently smooth so that ϕ(t, Ω) is a bounded domain with a Lipschitz-continuous boundary ∂ϕ(t, Ω) and ∂ϕ(t, Ω) = ϕ(t, ∂Ω). If x ϕ := x + u(t, x) denotes the current position of the point x at the time t, then the unit outward normal vector n ϕ along ∂ϕ(t, Ω) can be expressed according to [4] as

n ϕ (x ϕ ) = (I + ∇u(t, x)) -⊤ n(x) (I + ∇u(t, x)) -⊤ n(x) .
From here and (2), the Cauchy stress vector T ϕ , T ϕ (t, x ϕ ) = σ ϕ (t, x ϕ )n ϕ (x ϕ ), can be written as

T ϕ (t, x ϕ ) = 1 det(I + ∇u(t, x)) (I + ∇u(t, x)) -⊤ n(x) T (t, x)
and it is readily seen that the contact conditions (1) 6-8 at (t, x) are equivalent to

g(x ϕ ) ≤ 0, T ϕ ν (t, x ϕ ) ≤ 0, g(x ϕ )T ϕ ν (t, x ϕ ) = 0, |T ϕ τ (t, x ϕ )| ≤ -F ϕ (x ϕ )T ϕ ν (t, x ϕ ), T ϕ τ (t, x ϕ ) = F ϕ (x ϕ )T ϕ ν (t, x ϕ ) v ϕ τ (t, x ϕ ) |v ϕ τ (t, x ϕ )| if v ϕ τ (t, x ϕ ) = 0.
Here, v ϕ τ (t, x ϕ ) = uτ (t, x) is the tangent velocity of the point with the current position x ϕ , F ϕ (x ϕ ) = F (x) and

T ϕ ν (t, x ϕ ) = T ϕ (t, x ϕ ) • ν(x ϕ ), T ϕ τ (t, x ϕ ) = T ϕ (t, x ϕ ) • τ (x ϕ ).
In this way, we obtain the unilateral condition and the Coulomb friction law in the part ϕ(t, Γ c ) of the boundary ∂ϕ(t, Ω) of the deformed configuration.

Discretisation and the Continuation Problem

Spatial discretisation of our problem is done by applying a Lagrange finite-element method to a mixed variational formulation of (1) with Lagrange multipliers enforcing the Dirichlet and contact boundary conditions. In particular, the contact conditions are discretised similarly as in the example of the "almost conformal" discretisation of the displacement in [12, Section 4].

Let {φ i } 1≤i≤2n Ω be a basis of shape functions for displacements and {x i } 1≤i≤n Ω be the set of the corresponding finite-element nodes such that

φ i (x j ) = (δ i,2j-1 , δ i,2j ),
where δ denotes the Kronecker delta. Let {ξ i } 1≤i≤2n D be a basis of shape functions on Γ D . By {y i } 1≤i≤nc , we denote the nodes from {x i } 1≤i≤n Ω lying on Γ c \ Γ D and by {η i } 1≤i≤nc , the basis of shape functions on Γ c formed by the restrictions of appropriate component functions of

φ i onto Γ c such that η i (y j ) = δ ij .
The local and global numeration of the nodes on Γ c \ Γ D are linked by the mapping Θ : {1, . . . , n c } → {1, . . . , n Ω } such that

y i = x Θ(i) .
In what follows, we shall assume that the finite-element partitions corresponding to the bases {φ i } 1≤i≤2n Ω and {ξ i } 1≤i≤2n D are compatible with the decomposition of ∂Ω into Γ D , Γ N and Γ c . We introduce the following vectors and matrices for w = (w i ) 1≤i≤2n Ω and t ∈ [0, T ]:

A(w) = (A i (w)) 1≤i≤2n Ω , A i (w) = Ω σ x, I + 2n Ω j=1 w j ∇φ j (x) : ∇φ i (x) dx, B D = (B D,ij ) 1≤i≤2n D 1≤j≤2n Ω , B D,ij = Γ D ξ i (x) • φ j (x) ds, B ν (w) = (B ν,ij (w)) 1≤i≤nc 1≤j≤2n Ω , B ν,ij (w) =      ν 1 (y i + (w 2Θ(i)-1 , w 2Θ(i) )) if j = 2Θ(i) -1, ν 2 (y i + (w 2Θ(i)-1 , w 2Θ(i) )) if j = 2Θ(i), 0 otherwise, B τ (w) = (B τ,ij (w)) 1≤i≤nc 1≤j≤2n Ω , B τ,ij (w) =      τ 1 (y i + (w 2Θ(i)-1 , w 2Θ(i) )) if j = 2Θ(i) -1, τ 2 (y i + (w 2Θ(i)-1 , w 2Θ(i) )) if j = 2Θ(i), 0 otherwise, l(t, w) = (l i (t, w)) 1≤i≤2n Ω , l i (t, w) = Ω det I + 2n Ω j=1 w j ∇φ j (x) f ϕ t, x + 2n Ω j=1 w j φ j (x) • φ i (x) dx + Γ N det I + 2n Ω j=1 w j ∇φ j (x) I + 2n Ω j=1 w j ∇φ j (x) -⊤ n(x) • h ϕ t, x + 2n Ω j=1 w j φ j (x) • φ i (x) ds, u D (t) = (u D,i (t)) 1≤i≤2n D , u D,i (t) = Γ D U D (t, x) • ξ i (x) ds, g(w) = (g i (w)) 1≤i≤nc , g i (w) = g(y i + (w 2Θ(i)-1 , w 2Θ(i) )), F F F = (F i ) 1≤i≤nc , F i = F (y i ).
Furthermore, we take u 0 = (u 0 i ) 1≤i≤2n Ω such that 2n Ω i=1 u 0 i φ i approximates U 0 and we define the sets:

Λ ν = R nc -≡ {µ ν = (µ ν,i ) 1≤i≤nc ; µ ν,i ≤ 0}, Λ τ (F F F µ ν ) = {µ τ = (µ τ,i ) 1≤i≤nc ; |µ τ,i | ≤ -F i µ ν,i }, µ ν ∈ Λ ν .
Then, the spatial semi-discretisation of (1) can be formulated as follows:

Find u : [0, T ] → R 2n Ω , λ D : [0, T ] → R 2n D , λ ν , λ τ : [0, T ] → R nc such that u(0) = u 0 and for any t ∈ (0, T ), it holds that λ ν (t) ∈ Λ ν , λ τ (t) ∈ Λ τ (F F F λ ν (t)) and A(u(t)) -B ⊤ D λ D (t) -B ν (u(t)) ⊤ λ ν (t) -B τ (u(t)) ⊤ λ τ (t) = l(t, u(t)), B D u(t) = u D (t), (µ ν -λ ν (t)) ⊤ g(u(t)) ≥ 0, ∀µ ν ∈ Λ ν , (µ τ -λ τ (t)) ⊤ B τ (u(t)) u(t) ≥ 0, ∀µ τ ∈ Λ τ (F F F λ ν (t)),                      (3)
where u is the vector of nodal displacements, λ D is the Lagrange multiplier corresponding to the Dirichlet condition and λ ν , λ τ are the normal and tangential Lagrange multipliers on the contact zone, respectively. Approximating u(t) by a backward difference and employing the idea from the abstract frame [16, Section 2], one obtains the continuation problem:

Find γ ∈ R, u ∈ R 2n Ω , λ D ∈ R 2n D , λ ν ∈ Λ ν , λ τ ∈ Λ τ (F F F λ ν ) such that A(u) -B ⊤ D λ D -B ν (u) ⊤ λ ν -B τ (u) ⊤ λ τ = L(γ, u), B D u = U D (γ), (µ ν -λ ν ) ⊤ g(u) ≥ 0, ∀µ ν ∈ Λ ν , (µ τ -λ τ ) ⊤ B τ (u)u -B τ (u)w + (1 -γ)V τ ≥ 0, ∀µ τ ∈ Λ τ (F F F λ ν ), with L(γ, u) = γl(t k+1 , u) + (1 -γ)l(t k , u), U D (γ) = γu D (t k+1 ) + (1 -γ)u D (t k ), w = u k , V τ = (t k+1 -t k ) B τ (u k )u k -B τ (u k )u k-1 t k -t k-1 .
As in the abstract frame, the benefit of this continuation problem is that any solution from the kth time level of the corresponding fully discretised problem yields always a solution of the continuation problem, which can be taken as a starting point for continuation. Moreover, any solution of the continuation problem with γ = 1 yields a solution of the discretised problem on the (k + 1)st time level.

Due to variational characterisation of a projection onto a closed convex set, this problem is equivalent to

Find (γ, u, λ D , λ ν , λ τ ) ∈ R 1+2(n Ω +n D +nc) such that H(γ, u, λ D , λ ν , λ τ ) = 0, (P)
where H : R 1+2(n Ω +n D +nc) → R 2(n Ω +n D +nc) is defined by

H(γ, u, λ D , λ ν , λ τ ) =       A(u) -L(γ, u) -B ⊤ D λ D -B ν (u) ⊤ λ ν -B τ (u) ⊤ λ τ B D u -U D (γ) λ ν,i -(λ ν,i -g i (u)) -, i = 1, . . . , n c λ τ,i -P [F i (λ ν,i -g i (u)) -,-F i (λ ν,i -g i (u)) -] λ τ,i -(B τ (u)u -B τ (u)w) i -(1 -γ)V τ,i , i = 1, . . . , n c       (see also [2, Section 3]).
Here, (a) -is the non-positive part of a, or equivalently, the projection of a onto the interval (-∞, 0], and 2) and a composition of two P C 1 -functions is again P C 1 , H is a P C 1 -function under the following assumption: Assumption 1. Let A and l(t, .), t ∈ [0, T ], be C 1 -functions and g demarcating the rigid obstacle be of the class C 2 .

P [a,b] stands for the projection onto an interval [a, b]. Since the functions (a, b) → a -(a -b) -and (a, b, c) → a -P [kc -,-kc -] (a -b), k ≥ 0, are both P C 1 (see Figure
In this case, the general theory from [16] can be applied.

Remark 1. In the present discretisation, we treat the Dirichlet condition with the Lagrange multiplier so that it can be parametrised in a simple way. However, if the Dirichlet condition does not depend on the parameter (or more precisely, on the time), it can be prescribed directly in the discrete problem without any principal changes.

First-Order System

Let (γ, ū, λD , λν , λτ ) be a fixed solution of (P). We shall now give a particular expression for the first-order system

H ′ ((γ, ū, λD , λν , λτ ); (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ )) = 0
(compare [16, problem (P ′ )]) provided that Assumption 1 is satisfied. The dependence upon γ, ū, λD , λν and λτ will not be explicitly indicated in what follows usually.

After tedious applications of the chain rule, the first-order equilibrium equation corresponding to the first block of H reads

Ku ′ -γ ′ P -B ⊤ D λ ′ D -B ⊤ ν λ ′ ν -B ⊤ τ λ ′ τ = 0, (P ′ 1)
where we denote B ν = B ν (ū), B τ = B τ (ū), P = ∂L ∂γ (γ, ū), and

K = ∇A(ū) -∇ u L(γ, ū) -K ν (ū, λν ) -K τ (ū, λτ )
is the 2n Ω × 2n Ω tangent stiffness matrix with

K ν,ij (ū, λν ) =                ∂ν 1 ∂x 1 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λν,k if i = j = 2Θ(k) -1, ∂ν 1 ∂x 2 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λν,k if i = 2Θ(k) -1, j = 2Θ(k), ∂ν 2 ∂x 1 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λν,k if i = 2Θ(k), j = 2Θ(k) -1, ∂ν 2 ∂x 2 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λν,k if i = j = 2Θ(k), 0 otherwise, K τ,ij (ū, λτ ) =                ∂τ 1 ∂x 1 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λτ,k if i = j = 2Θ(k) -1, ∂τ 1 ∂x 2 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λτ,k if i = 2Θ(k) -1, j = 2Θ(k), ∂τ 2 ∂x 1 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λτ,k if i = 2Θ(k), j = 2Θ(k) -1, ∂τ 2 ∂x 2 y k + (ū 2Θ(k)-1 , ū2Θ(k) ) λτ,k if i = j = 2Θ(k), 0 otherwise, ∂ν k ∂x m = 1 ∇g 2 l=1 ∂ 2 g ∂x l ∂x m τ k τ l , ∂τ k ∂x m = - 1 ∇g 2 l=1 ∂ 2 g ∂x l ∂x m ν k τ l .
Remark 2. The partial derivatives of the normal and tangent are connected with the obstacle curvature κ: one has

∂ν 1 ∂x 1 + ∂ν 2 ∂x 2 = -κ
(curvature formulae for implicit planar curves can be found, for instance, in [6, Section 3]).

Furthermore, the first-order Dirichlet condition is

B D u ′ -γ ′ U ′ D = 0 (P ′ 2) with U ′ D = U ′ D (γ).
To derive the first-order contact conditions, one can show by other applications of the chain rule that

g ′ i (ū; u ′ ) = (B ν u ′ ) i , (B τ (.) . -B τ (.)w) i ′ (ū; u ′ ) = (C τ u ′ ) i ,
where

C τ = C ū τ (ū) + B τ -C w τ (ū) ∈ R nc×2n Ω , C w τ,ij (ū) =                ∂τ 1 ∂x 1 y i + (ū 2Θ(i)-1 , ū2Θ(i) ) w 2Θ(i)-1 + ∂τ 2 ∂x 1 y i + (ū 2Θ(i)-1 , ū2Θ(i) ) w 2Θ(i) if 2Θ(i) -1 = j, ∂τ 1 ∂x 2 y i + (ū 2Θ(i)-1 , ū2Θ(i) ) w 2Θ(i)-1 + ∂τ 2 ∂x 2 y i + (ū 2Θ(i)-1 , ū2Θ(i) ) w 2Θ(i) if 2Θ(i) = j, 0 otherwise. a -b = 0 a b 0 F (a, b) = a F (a, b) = b (a) F (a, b) = a -(a -b) -. a b c 0 G k (a, b, c) = a G k (a, b, c) = b G k (a, b, c) = a -kc G k (a, b, c) = a + kc a -b = -kc- a -b = kc- (b) G k (a, b, c) = a -P [kc-,-kc-] (a -b) for k > 0.
Figure 2: Subdivisions for the P C 1 -functions F and G k defined by ( 4) and ( 5).

Moreover, we define P C 1 -functions F : R 2 → R and G k : R 3 → R for any k ≥ 0 by

F (a, b) := a -(a -b) -= a if a -b ≥ 0, b if a -b ≤ 0, (4)
G k (a, b, c) := a -P [kc -,-kc -] (a -b) =          a if c ≥ 0, b if c ≤ 0, |a -b| ≤ -kc, a -kc if c ≤ 0, a -b ≤ kc, a + kc if c ≤ 0, a -b ≥ -kc (5) (see Figure 2). It is readily seen that for any a, b, c, a ′ , b ′ , c ′ ∈ R, F ′ ((a, b); (a ′ , b ′ )) =      a ′ if a -b > 0, b ′ if a -b < 0, a ′ -(a ′ -b ′ ) -if a -b = 0, G ′ 0 ((a, b, c); (a ′ , b ′ , c ′ )) = a ′ and if k > 0, G ′ k ((a, b, c); (a ′ , b ′ , c ′ )) =                                    a ′ if c > 0, b ′ if c < 0, |a -b| < -kc, a ′ -kc ′ if c < 0, a -b < kc, a ′ + kc ′ if c < 0, a -b > -kc, a ′ -k(c ′ ) - if c = 0, a -b < 0, a ′ + k(c ′ ) - if c = 0, a -b > 0, a ′ -P [kc ′ ,+∞) (a ′ -b ′ ) if c < 0, a -b = kc, a ′ -P (-∞,-kc ′ ] (a ′ -b ′ ) if c < 0, a -b = -kc, a ′ -P [k(c ′ ) -,-k(c ′ ) -] (a ′ -b ′ ) if c = a -b = 0.
Hence, Proposition 1 gives the first-order contact conditions in the form:

λ ′ ν,i = 0 if λν,i -g i (ū) > 0, (B ν u ′ ) i = 0 if λν,i -g i (ū) < 0, λ ′ ν,i -(λ ′ ν,i -(B ν u ′ ) i ) -= 0 if λν,i -g i (ū) = 0,      (P ′ 3) λ ′ τ,i = 0 if F i = 0 or F i > 0, λν,i -g i (ū) > 0, (C τ u ′ ) i -γ ′ V τ,i = 0 if F i > 0, λν,i -g i (ū) < 0, | λτ,i -vτ,i -(1 -γ)V τ,i | < -F i ( λν,i -g i (ū)), λ ′ τ,i -F i (λ ′ ν,i -(B ν u ′ ) i ) = 0 if F i > 0, λν,i -g i (ū) < 0, λτ,i -vτ,i -(1 -γ)V τ,i < F i ( λν,i -g i (ū)), λ ′ τ,i + F i (λ ′ ν,i -(B ν u ′ ) i ) = 0 if F i > 0, λν,i -g i (ū) < 0, λτ,i -vτ,i -(1 -γ)V τ,i > -F i ( λν,i -g i (ū)), λ ′ τ,i -F i (λ ′ ν,i -(B ν u ′ ) i ) -= 0 if F i > 0, λν,i -g i (ū) = 0, λτ,i -vτ,i -(1 -γ)V τ,i < 0, λ ′ τ,i + F i (λ ′ ν,i -(B ν u ′ ) i ) -= 0 if F i > 0, λν,i -g i (ū) = 0, λτ,i -vτ,i -(1 -γ)V τ,i > 0, λ ′ τ,i -P [F i (λ ′ ν,i -(Bν u ′ ) i ),+∞) (λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i ) = 0 if F i > 0, λν,i -g i (ū) < 0, λτ,i -vτ,i -(1 -γ)V τ,i = F i ( λν,i -g i (ū)), λ ′ τ,i -P (-∞,-F i (λ ′ ν,i -(Bν u ′ ) i )] (λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i ) = 0 if F i > 0, λν,i -g i (ū) < 0, λτ,i -vτ,i -(1 -γ)V τ,i = -F i ( λν,i -g i (ū)), λ ′ τ,i -P [F i (λ ′ ν,i -(Bν u ′ ) i ) -,-F i (λ ′ ν,i -(Bν u ′ ) i ) -] (λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i ) = 0 if F i > 0, λν,i -g i (ū) = λτ,i -vτ,i -(1 -γ)V τ,i = 0,                                                                              (P ′ 4)
where vτ,i = (B τ (ū)ū -B τ (ū)w) i and i passes {1, . . . , n c }.

To sum up, the complete first-order system for the problem (P) consists in finding ) satisfying (P ′ 1), (P ′ 2), (P ′ 3) and (P ′ 4).

(γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) ∈ R 1+2(n Ω +n D +nc
Unfortunately, Proposition 2 from [16] guaranteeing existence and uniqueness of solutions of a general P C 1 first-order system does not seem to be very useful for analysing this particular system. We can refer the reader to [1] and [15] for attempts to employ some results from theory of P C 1 -functions to discrete contact problems. For this reason, we shall present some other forms of the first-order system.

To start with, we shall reformulate the contact conditions with the aid of the following lemma, which can be proved by a brief discussion. By Sgn, we denote the multifunction defined by

Sgn(b) = b |b| if b = 0, [-1, 1] if b = 0.
Lemma 1. For any a, b, c ∈ R and any k ≥ 0, it holds that

(i) a -(a -b) -= 0 ⇐⇒ a ≤ 0, b ≤ 0, ab = 0; (ii) a -P [kc,+∞) (a -b) = 0 ⇐⇒ a -kc ≥ 0, b ≥ 0, (a -kc)b = 0; (iii) a -P (-∞,-kc] (a -b) = 0 ⇐⇒ a + kc ≤ 0, b ≤ 0, (a + kc)b = 0; (iv) a -P [kc -,-kc -] (a -b) = 0 ⇐⇒ a ∈ kc -Sgn(b).
Next, we introduce the index sets (compare (P ′ 3) and (P ′ 4)):

I = {1, . . . , n c }, I f = {i ∈ I; λν,i -g i (ū) > 0}, I c = {i ∈ I; λν,i -g i (ū) < 0}, I z = {i ∈ I; λν,i -g i (ū) = 0}, I 0 = {i ∈ I; F i = 0}, I s = {i ∈ I; F i > 0, | λτ,i -vτ,i -(1 -γ)V τ,i | < -F i ( λν,i -g i (ū)) -}, I l = {i ∈ I; F i > 0, | λτ,i -vτ,i -(1 -γ)V τ,i | > -F i ( λν,i -g i (ū)) -}, I i = {i ∈ I; F i > 0, | λτ,i -vτ,i -(1 -γ)V τ,i | = -F i ( λν,i -g i (ū)) -}.                    (6)
Observe that the contact conditions corresponding to the last two blocks of components of H can be rewritten for (γ, ū, λD , λν , λτ ) according to Lemma 1 as

λν,i ≤ 0, g i (ū) ≤ 0, λν,i g i (ū) = 0, i ∈ I, λτ,i ∈ F i λν,i Sgn(v τ,i + (1 -γ)V τ,i ), i ∈ I,
which can be viewed as contact conditions at each contact node. Consequently, one obtains

I f = {i ∈ I; λν,i = 0, g i (ū) < 0} (nodes not in contact -free), I c = {i ∈ I; λν,i < 0, g i (ū) = 0} (nodes in strong contact), I z = {i ∈ I; λν,i = g i (ū) = 0}
(nodes in grazing contact -with zero contact forces),

I s = {i ∈ I; F i > 0, | λτ,i | < -F i λν,i , vτ,i + (1 -γ)V τ,i = 0} (nodes in strong stick), I l = {i ∈ I; F i > 0, | λτ,i | = -F i λν,i , vτ,i + (1 -γ)V τ,i = 0} (nodes in non-zero slip), I i = {i ∈ I; F i > 0, | λτ,i | = -F i λν,i , vτ,i + (1 -γ)V τ,i = 0} (nodes in impending slip) and I s ⊂ I c .
To conclude, (P ′ 3) and (P ′ 4) are equivalent to

λ ′ ν,i = 0, i ∈ I f , (B ν u ′ ) i = 0, i ∈ I c , λ ′ ν,i ≤ 0, (B ν u ′ ) i ≤ 0, λ ′ ν,i (B ν u ′ ) i = 0, i ∈ I z ,      (P ′ 3') λ ′ τ,i = 0, i ∈ I f ∪ I 0 , (C τ u ′ ) i = γ ′ V τ,i , i ∈ I s , λ ′ τ,i = sgn(v τ,i + (1 -γ)V τ,i )F i λ ′ ν,i , i ∈ (I c ∪ I z ) ∩ I l , sgn( λτ,i )((C τ u ′ ) i -γ ′ V τ,i ) ≤ 0, sgn( λτ,i )λ ′ τ,i + F i λ ′ ν,i ≤ 0, sgn( λτ,i )((C τ u ′ ) i -γ ′ V τ,i )(sgn( λτ,i )λ ′ τ,i + F i λ ′ ν,i ) = 0, i ∈ I c ∩ I i , λ ′ τ,i ∈ F i λ ′ ν,i Sgn((C τ u ′ ) i -γ ′ V τ,i ), i ∈ I z ∩ I i .                      (P ′ 4')
The system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4') is close to a rate problem of a quasi-static contact problem discretised in space only (this problem can be found in [17], for example). The main difference between the two problems is that in our continuation problem, we suppress solutions in which slipping contact nodes become stuck when loading of the body changes to unloading suddenly. Indeed, there is no equivalent of the index set I l in the rate problem and all nodes in contact are treated either as in stick or as in possible slip in the near future. Of course, this causes the rate problem to be more complicated (compare (P ′ 4') for i ∈ I l and i ∈ I i ).

Let us also mention that due to nodal discretisation of the contact conditions considered in this paper, the continuation problem could be formulated without the term (1γ)V τ approximating a multiple of an initial velocity in the frictional condition and (0, u k , λ k D , λ k ν , λ k τ ) would remain its solution (which is not always the case for other types of discretisation). However, when starting continuation from this point, one would run into an analogous situation as in the rate problem. More precisely, one would have vτ = 0 for (γ, ū, λD , λν , λτ ) = (0,

u k , λ k D , λ k ν , λ k τ ), the indices i with (B τ (u k )u k -B τ (u k )u k-1 ) i = 0
, which correspond to slipping nodes at the time t k , would fall into I i instead of I l and I l would be empty.

By the way, the situation for continuation from (0,

u k , λ k D , λ k ν , λ k τ )
in the problem without (1γ)V τ is very similar to what happens in the semi-smooth (or piecewise-smooth) Newton method applied directly to a projection formulation of the fully discretised problem for the time t k+1 with the initial approximation (u k , λ k D , λ k ν , λ k τ ). Regarding the Newton method as an active-set strategy, one can employ adaptations of the index sets defined above to determine an element from the generalised gradient used in each Newton's iteration and one can see that any slipping node at the time t k makes the initial approximation (u k , λ k D , λ k ν , λ k τ ) a non-smooth point of the equation being solved! Example 1. Let us consider the simple example from [9, Section 4] modified to the quasistatic case (which is nothing else than a special case of the example with Coulomb's law

Γ D Γ N h Γ c O x 1 x 2 0 Ω Figure 3: Geometry from Example 1.
of friction in [13]): contact of a single P 1 triangle with a flat foundation in the framework of linearised elasticity (Figure 3). The volume forces and the initial state are supposed to be null, the triangle is fixed along Γ D and the Dirichlet condition is prescribed directly in this example.

The problem can be formulated as follows (compare (3) and (P)):

Find u 1 , u 2 , λ ν , λ τ : [0, T ] → R such that (u 1 (0), u 2 (0), λ ν (0), λ τ (0)) = (0, 0, 0, 0)
and for any t ∈ (0, T ), it holds that

a -b -b a u 1 (t) u 2 (t) - l 1 (t) l 2 (t) - λ ν (t) λ τ (t) = 0 0 , λ ν (t) -(λ ν (t) -u 1 (t)) -= 0, λ τ (t) -P [F (λν (t)-u 1 (t)) -,-F (λν (t)-u 1 (t)) -] (λ τ (t) -u2 (t)) = 0,
where a := (λ + 3µ)/2 and b := (λ + µ)/2 are given by Lamé constants λ and µ. The loading is considered to be (l 1 (t), l 2 (t)) = (0, t).

According to the analysis in [13], the unique solution of this problem is

u 1 (t) = 0, u 2 (t) = l 2 (t) -F l 1 (t) a + F b = t a + F b , λ ν (t) = - al 1 (t) + bl 2 (t) a + F b = - bt a + F b , λ τ (t) = F λ ν (t) = - F bt a + F b
and it corresponds to slip of the only contact node to the right. The rate problem at each solution point for t > 0 is:

Find ( u1 , u2 , λν , λτ ) ∈ R 4 such that a -b -b a u1 u2 - 0 1 -λν λτ = 0 0 , u1 = 0, u2 ≥ 0, -λτ + F λν ≤ 0, u2 (-λτ + F λν ) = 0,                (7)
that is, one has to deal with a complementarity between slip and stick at each time (in fact, the contact node would get stuck if the loading changed to unloading in the opposite direction).

On the other hand, solutions of the fully discretised (also called incremental) problems, which can be derived by analysing all possible cases similarly as in [9] for the static case, are:

u k 1 = 0, u k 2 = l 2 (t k ) -F l 1 (t k ) a + F b = t k a + F b , λ k ν = - al 1 (t k ) + bl 2 (t k ) a + F b = - bt k a + F b , λ k τ = F λ k ν = - F bt k a + F b
and the continuation problem proposed in this paper reads as follows:

Find (γ, u 1 , u 2 , λ ν , λ τ ) ∈ R 5 such that a -b -b a u 1 u 2 - L 1 (γ) L 2 (γ) - λ ν λ τ = 0 0 , λ ν -(λ ν -u 1 ) -= 0, λ τ -P [F (λν -u 1 ) -,-F (λν -u 1 ) -] (λ τ -(u 2 -u k 2 ) -(1 -γ)V τ ) = 0.                (8)
Here, (L 1 (γ), L 2 (γ)) = (0, t k + γ(t k+1t k )) and

V τ = t k+1 -t k t k -t k-1 (u k 2 -u k-1 2 )
is positive for any k ≥ 1. Solutions of this problem can be computed analogously as the solutions of the incremental problems and they form the set:

(γ, u 1 , u 2 , λ ν , λ τ ) ∈ R 5 ; γ ∈ [0, 1], u 1 = 0, u 2 = L 2 (γ) -F L 1 (γ) a + F b = t k + γ(t k+1 -t k ) a + F b , λ ν = - aL 1 (γ) + bL 2 (γ) a + F b = - b(t k + γ(t k+1 -t k )) a + F b , λ τ = F λ ν = - F b(t k + γ(t k+1 -t k )) a + F b .
Then, the corresponding first-order system for any γ ∈ [0, 1] becomes:

Find (γ ′ , u ′ 1 , u ′ 2 , λ ′ ν , λ ′ τ ) ∈ R 5 such that a -b -b a u ′ 1 u ′ 2 - 0 γ ′ (t k+1 -t k ) - λ ′ ν λ ′ τ = 0 0 , u ′ 1 = 0, λ ′ τ = F λ ′ ν .
In comparison with the rate problem (7), the complementarity conditions are replaced by an equality here, which is apparently easier to handle.

Finally, solutions of the continuation problem without (1 -γ)V τ would remain the same as with it, however, the first-order system would contain the following frictional condition at γ = 0:

u ′ 2 ≥ 0, -λ ′ τ + F λ ′ ν ≤ 0, u ′ 2 (-λ ′ τ + F λ ′ ν ) = 0.
Therefore, it would be more complicated to find a suitable initial tangent vector for numerical continuation.

Due to the similarity of the first-order system to the rate problem mentioned above, we shall follow steps in [17] and we shall reformulate it as a mixed complementarity inclusion and an explicit complementarity problem. To this end, we shall need the following assumption.

Assumption 2. Let the matrix

B :=   B D B ν B τ   have full row rank.
One can easily verify that this assumption is equivalent to the following inf-sup condition:

∃β > 0 : inf 0 =(λ D ,λν ,λτ )∈R 2(n D +nc) sup 0 =u∈R 2n Ω u ⊤ (B ⊤ D λ D + B ⊤ ν λ ν + B ⊤ τ λ τ ) u (λ D , λ ν , λ τ ) ≥ β.
This is a reasonable assumption on the used spatial discretisation because it prevents the discretised problem from being over-constrained. If it is satisfied,

dim Im(B ⊤ ) = 2(n D + n c ), dim ker(B) = 2n Ω -2(n D + n c )
and one can define a non-singular matrix

T B := N B + ,
where the columns of the matrix N ∈ R 2n Ω ×(2n Ω -2(n D +nc)) form a basis of ker(B) and

B + = B ⊤ (BB ⊤ ) -1 is the right inverse of B. Then T ⊤ B B ⊤ D B ⊤ ν B ⊤ τ   λ ′ D λ ′ ν λ ′ τ   = BN BB + ⊤   λ ′ D λ ′ ν λ ′ τ   =     0 2n Ω -(2n D +2nc) λ ′ D λ ′ ν λ ′ τ     . (9)
Remark 3. For the nodal discretisation of contact conditions considered in this paper, the rows of Bν Bτ are always linearly independent and they are orthogonal to the rows of B D . Therefore, Assumption 2 is satisfied provided that B D has full row rank, which is equivalent to the inf-sup condition on the Dirichlet boundary solely:

∃β > 0 : inf 0 =λ D ∈R 2n D sup 0 =u∈R 2n Ω u ⊤ B ⊤ D λ D u λ D ≥ β.
If it is the case,

T B = N B + D Bν Bτ +
and the columns of N and Bν Bτ + form a basis of ker(B D ).

Similarly, we shall need the following additional assumption.

Assumption 3. Let the matrix

C :=   B D B ν C τ   have full row rank.
Let us note that if the rigid obstacle is plane, that is, ∇ 2 g = 0, then C τ = B τ and this assumption is always satisfied provided that Assumption 2 is fulfilled. Under Assumption 3, one can introduce a non-singular transformation matrix

T C := N C +
with the columns of N constituting a basis of ker(C) and C + = C ⊤ (CC ⊤ ) -1 (invoking the introduction of C, one can easily verify that N may be taken the same as above). Consequently, any displacement direction u ′ ∈ R 2n Ω in the first-order system can be represented by a vector ũ′ of the same length in the following way:

u ′ = T C ũ′ , ũ′ =     u ′ u u ′ D u ′ ν u ′ τ     , (10)
where u ′ D ∈ R 2n D is a part imposed to the first-order Dirichlet condition, u ′ ν , u ′ τ ∈ R nc correspond to the normal and tangential components on the contact zone and u ′ u ∈ R 2n Ω -2(n D +nc) is free of the first-order Dirichlet and contact conditions (unconstrained). In fact,

  B D u ′ B ν u ′ C τ u ′   =   B D B ν C τ   N C + ũ′ =   u ′ D u ′ ν u ′ τ   . ( 11
)
Making use of ( 9) and ( 10) and introducing further partitioning of ũ′ , λ ′ ν and λ ′ τ , we can write (P ′ 1) as

                   u ′ u u ′ G u ′ N u ′ Z u ′ νzl u ′ τ cl u ′ τ zl u ′ ci             =             0 λ ′ G λ ′ N λ ′ Z λ ′ νzl λ ′ τ cl λ ′ τ zl λ ′ ci             + γ ′              P u P G P N P Z P νzl P τ cl P τ zl P ci              , (12)
where K = T ⊤ B KT C , P = T ⊤ B P , G denotes jointly the indices D, νs, νc0, νcl and τ s, N the indices νf, τ f, τ c0 and τ z0, ci the indices νci and τ ci, Z the indices νz0, νzi and τ zi, and s stands for the indices from I s , c0 for the indices from I c ∩ I 0 , cl for the indices from I c ∩ I l and so forth.

The values of the components of u ′ G and u ′ νci are given by (P ′ 2), (P ′ 3') and (P ′ 4') combined with (11):

u ′ D = γ ′ U ′ D , u ′ νs = 0, u ′ νc0 = 0, u ′ νcl = 0, u ′ τ s = γ ′ V τ s , u ′ νci = 0
and λ ′ G can be expressed directly from (12). In particular,

λ ′ νcl = Kνcl,u u ′ u + Kνcl,N u ′ N + Kνcl,Z u ′ Z + Kνcl,νzl u ′ νzl + Kνcl,τcl u ′ τ cl + Kνcl,τzl u ′ τ zl + Kνcl,τci u ′ τ ci -γ ′ ( P νcl -Kνcl,D U ′ D -Kνcl,τs V τ s ).
Further, one has from (P ′ 3') and (P ′ 4') that

λ ′ N = 0, λ ′ τ cl = S cl F cl λ ′ νcl , λ ′ τ zl = S zl F zl λ ′ νzl ,
where

S cl = Diag(sgn(v τ,i + (1 -γ)V τ,i )), i ∈ I c ∩ I l , F cl = Diag(F i ), i ∈ I c ∩ I l , S zl = Diag(sgn(v τ,i + (1 -γ)V τ,i )), i ∈ I z ∩ I l , F zl = Diag(F i ), i ∈ I z ∩ I l .
Hence, ( 12) can be transformed into

             Ku,
                         u ′ u u ′ N u ′ Z u ′ νzl u ′ τ cl u ′ τ zl 0 u ′ τ ci             =           0 0 λ ′ Z λ ′ νzl 0 S zl F zl λ ′ νzl λ ′ ci           + γ ′              P u P N P Z P νzl P τ cl P τ zl P ci              (13) with              Ku,.
KN,.

KZ,.

Kνzl,.

Kτcl,.

Kτzl,.

Kci,.

             =              Ku,.
KN,.

KZ,.

Kνzl,.

Kτcl,. -S cl F cl Kνcl,.

Kτzl,.

Kci,.

             , P =            P u -Ku,D U ′ D -Ku,τs V τ s P N -KN,D U ′ D -KN,τs V τ s P Z -KZ,D U ′ D -KZ,τs V τ s P νzl -Kνzl,D U ′ D -Kνzl,τs V τ s P τ cl -Kτcl,D U ′ D -Kτcl,τs V τ s -S cl F cl ( P νcl -Kνcl,D U ′ D -Kνcl,τs V τ s ) P τ zl -Kτzl,D U ′ D -Kτzl,τs V τ s P ci -Kci,D U ′ D -Kci,τs V τ s           
.

Denoting the indices u, N, Z, νzl, τ cl and τ zl jointly by E, we want to eliminate u ′ E now. This leads us to the following assumption.

Assumption 4. Let the matrix KE,E be invertible.

As we shall show in the proof of Corollary 1, this holds, for example, for a model of linearised elasticity with a properly constrained body in contact with a plane rigid obstacle provided that the coefficient of friction is sufficiently small. The assumption enables us to express u ′ E from the first six rows of ( 13):

u ′ E = E E,Z λ ′ Z + ĒE,νzl λ ′ νzl + γ ′ E E,E P E -E E,E KE,τci u ′ τ ci , (14)
where

E E,E = K-1 E,E , ĒE,νzl = E E,νzl + E E,τ zl S zl F zl .
Then, substituting for it into the last row of ( 13), we obtain

λ ′ ci = -γ ′ P ci + Kci,E E E,Z λ ′ Z + Kci,E ĒE,νzl λ ′ νzl + Kci,τci u ′ τ ci with P ci = P ci -Kci,E E E,E P E , Kci,τci = Kci,τci -Kci,E E E,E KE,τci .
From here, we have

S ci λ ′ τ ci -F ci λ ′ νci = -γ ′ P τ ci + Kτci,τci S 2 ci (u ′ τ ci -γ ′ V τ ci ) + Kτci,E E E,Z λ ′ Z + Kτci,E ĒE,νzl λ ′ νzl , (15) where
S ci = -Diag(sgn( λτ,i )), i ∈ I c ∩ I i , F ci = Diag(F i ), i ∈ I c ∩ I i , Kτci,τci = S ci Kτci,τci -F ci Kνci,τci , Kτci,E = S ci Kτci,E -F ci Kνci,E , P τ ci = S ci ( P τ ci -Kτci,τci V τ ci ) -F ci ( P νci -Kνci,τci V τ ci ).
Taking ( 15), ( 14) with the indices from νz0, νzl, νzi and τ zi, and the remaining conditions from (P ′ 3') and (P ′ 4'), we arrive at the mixed complementarity inclusion:

Find (x, y) ∈ R M × R M such that Dx = y + γ ′ b, ζ τ ci ≥ 0, ψ τ ci ≥ 0, ψ ⊤ τ ci ζ τ ci = 0, ζ νz0 ≥ 0, ψ νz0 ≥ 0, ψ ⊤ νz0 ζ νz0 = 0, ζ νzl ≥ 0, ψ νzl ≥ 0, ψ ⊤ νzl ζ νzl = 0, ζ νzi ≥ 0, ψ νzi ≥ 0, ψ ⊤ νzi ζ νzi = 0, ψ τ zi ∈ -S zi (ζ τ zi )F zi ψ νzi ,                          (16)
where

M = n ci + n z0 + n zl + 2n zi with n ci = #(I c ∩ I i ), n z0 = #(I z ∩ I 0 ), n zl = #(I z ∩ I l ) and n zi = #(I z ∩ I i ), x =       ζ τ ci ψ νz0 ψ νzl ψ νzi ψ τ zi       =       S ci (u ′ τ ci -γ ′ V τ ci ) -λ ′ νz0 -λ ′ νzl -λ ′ νzi λ ′ τ zi       , y =       ψ τ ci ζ νz0 ζ νzl ζ νzi ζ τ zi       =       S ci λ ′ τ ci -F ci λ ′ νci -u ′ νz0 -u ′ νzl -u ′ νzi u ′ τ zi -γ ′ V τ zi       , D =         Kτci,τci S ci -Kτci,E E E,νz0 -Kτci,E ĒE,νzl -Kτci,E E E,νzi Kτci,E E E,τ zi E νz0,E KE,τci S ci E νz0,νz0 Ēνz0,νzl E νz0,νzi -E νz0,τ zi E νzl,E KE,τci S ci E νzl,νz0 Ēνzl,νzl E νzl,νzi -E νzl,τ zi E νzi,E KE,τci S ci E νzi,νz0 Ēνzi,νzl E νzi,νzi -E νzi,τ zi -E τ zi,E KE,τci S ci -E τ zi,νz0 -Ēτzi,νzl -E τ zi,νzi E τ zi,τ zi         , b =         P τ ci E νz0,E ( P E -KE,τci V τ ci ) E νzl,E ( P E -KE,τci V τ ci ) E νzi,E ( P E -KE,τci V τ ci ) -E τ zi,E ( P E -KE,τci V τ ci ) + V τ zi         , F zi = Diag(F i ), i ∈ I z ∩ I i ,
and S zi : R n zi ⇉ R n zi ×n zi is a matrix-set-valued map defined by

S zi (ζ τ zi ) = Diag(Sgn((ζ τ zi ) 1 ), . . . , Sgn((ζ τ zi ) n zi ).
It is readily seen from the construction that there is a one-to-one correspondence between solutions of the first-order system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4') and of ( 16) for γ ′ ∈ R given.

To transform (16) into the announced explicit complementarity problem, we start with an easy observation.

Lemma 2. For any a, c ∈ R and any b, k ≥ 0, the following two conditions are equivalent:

(i) a ∈ -kb Sgn(c); (ii) -|a| + kb ≥ 0 & ac + kb|c| = 0.
Next, we define closed convex cones of admissible values of x and y from ( 16) (all inequalities have to be understood component-wise): Clearly, these objects enjoy the properties summarised below (see also Figure 4). Recall that a dual cone C * of a cone C ⊂ R M is introduced as

C x := {x = (ζ τ ci , ψ νz0 , ψ νzl , ψ νzi , ψ τ zi ) ∈ R M ; ζ τ ci ≥ 0, ψ νz0 ≥ 0, ψ νzl ≥ 0, ψ νzi ≥ 0, -|ψ τ zi | + F zi ψ νzi ≥ 0}, C y := {y = (ψ τ ci , ζ νz0 , ζ νzl , ζ νzi , ζ τ zi ) ∈ R M ; ψ τ ci ≥ 0, ζ νz0 ≥ 0, ζ νzl ≥ 0, ζ νzi ≥ 0, ζ τ zi ∈ R n zi } and a mapping h : R M → R M by h(ψ τ ci , ζ νz0 , ζ νzl , ζ νzi , ζ τ zi ) =       ψ τ ci ζ νz0 ζ νzl ζ νzi + F zi |ζ τ zi | ζ τ zi       . F i 1 (ψ νzi )i (ψ τ zi )i 0 Cx (a) 1 F i (ζ νzi )i (ζ τ zi )i 0 h(Cy) Cy (b)
C * = {y ∈ R M ; y ⊤ x ≥ 0, ∀x ∈ C }. Lemma 3. It holds that (i) h(C y ) = {y = (ψ τ ci , ζ νz0 , ζ νzl , ζ νzi , ζ τ zi ) ∈ R M ; ψ τ ci ≥ 0, ζ νz0 ≥ 0, ζ νzl ≥ 0, ζ νzi ≥ 0, ζ νzi -F zi |ζ τ zi | ≥ 0} = (C x ) * ; (ii) h is a bijection. Finally, we define a map f γ ′ : R M → R M by f γ ′ (x) = h(Dx -γ ′ b) for any γ ′ ∈ R fixed.
Altogether, one can easily verify that there is a one-to-one correspondence between solutions of ( 16) and the following explicit complementarity problem:

Find x ∈ C x such that f γ ′ (x) ∈ (C x ) * and x ⊤ f γ ′ (x) = 0. ( 17
)

Results of Existence and Uniqueness

We are prepared to adapt the existence and uniqueness results from [17].

Proposition 3. Let the matrix D be strictly co-positive in C x , that is,

x ⊤ Dx > 0, ∀x ∈ C x , x = 0. (18)
Then (17) possesses a solution for any γ ′ ∈ R and any b ∈ R M .

Proof. The assertion follows directly from [10, Corollary 4.13]. Since C x is a closed convex cone and f γ ′ is a bounded continuous finite-dimensional operator, it suffices to verify that f γ ′ is coercive in the following sense:

lim x →+∞ x∈Cx x ⊤ f γ ′ (x) x = +∞.
But this is ensured by ( 18) as it gives D 0 > 0 such that

x ⊤ Dx x 2 ≥ D 0 , ∀x ∈ C x , x = 0. Hence, ∀x ∈ C x , x = 0 : x ⊤ f γ ′ (x) x = 1 x x ⊤ (Dx -γ ′ b) + ψ ⊤ νzi F zi |(Dx -γ ′ b) τ zi | ≥ D 0 x -γ ′ b ,
where we have used the non-negativity of the elements of ψ νzi and F zi , and (Dxγ ′ b) τ zi denotes the last sub-vector of Dxγ ′ b.

Before stating a consequence of this result, we recall the concept of positive definiteness of a matrix on a subspace: Let M be a square matrix of order n and V be an m-dimensional subspace of R n . We say that M is positive definite on V if u ⊤ M u > 0 for any non-zero u ∈ V , or equivalently, if the matrix T ⊤ M T is positive definite for any matrix T ∈ R n×m whose columns form a basis of V .

Corollary 1 (Existence). Let (P) correspond to a model of linearised elasticity where the body is properly constrained by the Dirichlet condition, that is, the rows of B D are linearly independent and the stiffness matrix A is positive definite on ker(B D ). If g satisfies Assumption 1 and ∇ 2 g / ∇g and F i , i ∈ ((I c ∪ I z ) ∩ I l ) ∪ (I c ∩ I i ), are small, then the first-order system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4') has at least one solution for any γ ′ ∈ R, P ∈ R 2n Ω and U ′ D ∈ R 2n D . Proof. Suppose for a while that the assumptions are satisfied even with ∇ 2 g = 0 and

F i = 0 for i ∈ ((I c ∪ I z ) ∩ I l ) ∪ (I c ∩ I i ). Then K = A, C τ = B τ
and both Assumptions 2 and 3 hold according to Remark 3. Furthermore, K is a principal sub-matrix of N ⊤ D KN D , where the columns of N D are formed by the columns of N and Bν Bτ + constituting a basis of ker(B D ). Thus, K is symmetric positive definite and so are KE,E and its inverse E E,E . Consequently, Assumption 4 is fulfilled, the first-order system can be transformed into (17) by the procedure described above and

x ⊤ Dx = ζ ⊤ τ ci S ci Kτci,τci S ci ζ τ ci + ψ ⊤ νz0 ψ ⊤ νzl ψ ⊤ νzi -ψ ⊤ τ zi     E νz0,νz0 E νz0,νzl E νz0,νzi E νz0,τ zi E νzl,νz0 E νzl,νzl E νzl,νzi E νzl,τ zi E νzi,νz0 E νzi,νzl E νzi,νzi E νzi,τ zi E τ zi,νz0 E τ zi,νzl E τ zi,νzi E τ zi,τ zi         ψ νz0 ψ νzl ψ νzi -ψ τ zi     .
The matrix S ci Kτci,τci S ci = S -1 ci Kτci,τci S ci in the first term above is similar to Kτci,τci = Kτci,τci -Kτci,E E E,E KE,τci , which is a principal sub-matrix of the Schur complement of KE,E in K. Thus, it is positive definite and so is the principal sub-matrix of E in the second term above. This implies the positive definiteness of D, a fortiori fulfilment of ( 18) and application of the previous proposition is justified under the stronger assumptions.

Continuous dependence of the elements of C τ , K and D on the elements of ∇ 2 g/ ∇g and F i , i ∈ ((I c ∪ I z ) ∩ I l ) ∪ (I c ∩ I i ), completes the claim.

The uniqueness result uses the notion of a P -matrix, which is a matrix whose all principal minors are positive. Proposition 4. A sufficient condition for uniqueness of a solution of (16) for γ ′ ∈ R and b ∈ R M given is that the matrix

D       I 0 0 0 0 0 I 0 0 0 0 0 I 0 0 0 0 0 I 0 0 0 0 -Szi F zi I       is a P -matrix for any Szi = Diag(s 1 , . . . , sn zi ) with si ∈ [-1, 1].
Proof. Let (x 1 , y 1 ) and (x 2 , y 2 ) be two solutions of ( 16) for the same γ ′ b and

ψ i τ zi = - Si zi F zi ψ i νzi , Si zi ∈ S zi (ζ i τ zi ), i = 1, 2, (19)
(compare (16) 7 ). One can write

ψ 1 τ zi -ψ 2 τ zi = - S1 zi F zi (ψ 1 νzi -ψ 2 νzi ) -( S1 zi - S2 zi )F zi ψ 2 νzi ( 20
)
and taking the difference for the two solutions in (16) 2 , one obtains

      ψ 1 τ ci -ψ 2 τ ci ζ 1 νz0 -ζ 2 νz0 ζ 1 νzl -ζ 2 νzl ζ 1 νzi -ζ 2 νzi ζ 1 τ zi -ζ 2 τ zi       = D       I 0 0 0 0 0 I 0 0 0 0 0 I 0 0 0 0 0 I 0 0 0 0 - S1 zi F zi I             ζ 1 τ ci -ζ 2 τ ci ψ 1 νz0 -ψ 2 νz0 ψ 1 νzl -ψ 2 νzl ψ 1 νzi -ψ 2 νzi -( S1 zi - S2 zi )F zi ψ 2 νzi      
.

By the complementarity conditions (16) 3-6 , one has and by (19), monotonicity of the multifunction Sgn and non-negativity of F zi and ψ 2 νzi ,

(ψ 1 τ ci -ψ 2 τ ci ) ⊤ (ζ 1 τ ci -ζ 2 τ ci ) ≤ 0, (ζ 1 νz0 -ζ 2 νz0 ) ⊤ (ψ 1 νz0 -ψ 2 νz0 ) ≤ 0, (ζ 1 νzl -ζ 2 νzl ) ⊤ (ψ 1 νzl -ψ 2 νzl ) ≤ 0, (ζ 1 νzi -ζ 2 νzi ) ⊤ (ψ 1 νzi -ψ 2 νzi ) ≤ 0
-(ζ 1 τ zi -ζ 2 τ zi ) ⊤ ( S1 zi - S2 zi )F zi ψ 2 νzi ≤ 0.
Since a P -matrix does not reverse the sign of any non-zero vector ([5, 3.3.4 Theorem]), the imposed assumption and (20) imply that x 1 = x 2 and y 1 = y 2 . Example 3 (Large ∇ 2 g / ∇g ). Let us modify the simple continuation problem (8) to the case where the rigid foundation is demarcated by

Γ D Γ N h Γ c O x 1 x 2 0 Ω ν τ
g : x = (x 1 , x 2 ) → r 2 -(x 1 -r) 2 -x 2 2
for r > 0 fixed (Figure 5). Then

∇g(x) = -2 x 1 -r x 2 , ∇g(x) = 2 (x 1 -r) 2 + x 2 2 , ν(x) = -1 (x 1 -r) 2 + x 2 2 x 1 -r x 2 , τ (x) = 1 (x 1 -r) 2 + x 2 2 x 2 -(x 1 -r)
and the continuation problem becomes

a -b -b a u - L 1 (γ) L 2 (γ) -λ ν ν(u) -λ τ τ (u) = 0 0 , λ ν -(λ ν -g(u)) -= 0, λ τ -P [F (λν -g(u)) -,-F (λν -g(u)) -] λ τ -(τ (u) ⊤ u -τ (u) ⊤ w) -(1 -γ)V τ = 0.
Here, we consider (L 1 (γ), L 2 (γ)) = (1, -γ), w = 0, V τ = 0 as in (i) of the previous example and F > 0 arbitrarily fixed.

It is readily seen that the set given by ( 21) is a part of the solution set in this case too. If one takes (γ, ū1 , ū2 , λν , λτ ) = (F , 0, 0, -1, F ) (strong contact with impending slip), one has

ν(ū) = 1 0 , τ (ū) = 0 1 , ∇ 2 g(ū) = -2 0 0 -2 , ∇ 2 g(ū) ∇g(ū) = - 1 r 1 0 0 1 , K ν (ū, λν ) = 1 r 0 0 0 1 , K τ (ū, λτ ) = 1 r 0 F 0 0
and the corresponding first-order system reads

a -b -F r -b a -1 r u ′ 1 u ′ 2 + 0 γ ′ - λ ′ ν λ ′ τ = 0 0 , u ′ 1 = 0, u ′ 2 ≤ 0, λ ′ τ + F λ ′ ν ≤ 0, u ′ 2 (λ ′ τ + F λ ′ ν ) = 0.
then there are δ > 0 and a C

1 -curve c : [s, s + δ) → R × R N such that (j) c(s) = (γ, x); (jj) ∀s ∈ [s, s + δ) : H(c(s)) = 0; (jjj) c ′ (s) = (γ ′ , x ′ ).
Moreover, if (γ ′ , x ′ ) = (0, 0), the image of any other curve c : [s, s + δ) → R × R N with δ > 0 and such that

( j) c(s) = (γ, x); ( jj) ∀s ∈ [s, s + δ) : H(c(s)) = 0; ( jjj) c′ (s) ∈ r>0 r(γ ′ , x ′ ) (23)
coincides with the image of c in a neighbourhood of (γ, x). Now, let (γ, ū, λD , λν , λτ )(= (γ, x)) solve (P), that is, H(γ, ū, λD , λν , λτ ) = 0. Taking into account the continuity of H and using the index sets defined by ( 6), one can construct the set {H (j) } of selection functions of H at (γ, ū, λD , λν , λτ ) such that, after eventual reordering of components, it consists of functions (compare the functions F and G k introduced by ( 4) and ( 5)):

H (j) (γ, u, λ D , λ ν , λ τ ) =           A(u) -L(γ, u) -B ⊤ D λ D -B ν (u) ⊤ λ ν -B τ (u) ⊤ λ τ B D u -U D (γ) λ ν,i , i ∈ I f ∪ I f z g i (u), i ∈ I c ∪ I c z λ τ,i , i ∈ I f ∪ I f z ∪ I 0 (B τ (u)u -B τ (u)w) i + (1 -γ)V τ,i , i ∈ I s ∪ I s ci ∪ I cs zi λ τ,i + s i F i (λ ν,i -g i (u)), i ∈ (I c ∩ I l ) ∪ I c zl ∪ I l ci ∪ I cl+ zi ∪ I cl- zi           .
Here, in a similar spirit as in [3], the index sets I c z , I l ci , I cl+ zi and I cl- zi pass all possible combinations such that

I c z ⊂ I z , I l ci ⊂ (I c ∩ I i ), I cl+ zi , I cl- zi ⊂ (I c z ∩ I i ), I cl+ zi ∩ I cl- zi = ∅ (24)
and

I f z = I z \ I c z , I s ci = (I c ∩ I i ) \ I s ci , I cs zi = (I c z ∩ I i ) \ (I cl+ zi ∪ I cl- zi ), I c zl = I c z ∩ I l , (25)
s i =          -sgn(v τ,i + (1 -γ)V τ,i ) if i ∈ (I c ∩ I l ) ∪ I c zl , sgn( λτ,i ) if i ∈ I l ci , -1 if i ∈ I cl+ zi , 1 if i ∈ I cl- zi . (26)
It is worth noticing that I f z and I c z constitute a decomposition of I z , I s ci and I l ci constitute a decomposition of I c ∩ I i and I cs zi , I cl+ zi and I cl- zi constitute a decomposition of I c z ∩ I i .

Gradients of these selection functions at (γ, ū, λD , λν , λτ ) are given by

∇H (j) (γ, ū, λD , λν , λτ )(γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) =           Ku ′ -γ ′ P -B ⊤ D λ ′ D -B ⊤ ν λ ′ ν -B ⊤ τ λ ′ τ B D u ′ -γ ′ U ′ D λ ′ ν,i , i ∈ I f ∪ I f z (B ν u ′ ) i , i ∈ I c ∪ I c z λ ′ τ,i , i ∈ I f ∪ I f z ∪ I 0 (C τ u ′ ) i -γ ′ V τ,i , i ∈ I s ∪ I s ci ∪ I cs zi λ ′ τ,i + s i F i (λ ′ ν,i -(B ν u ′ ) i ), i ∈ (I c ∩ I l ) ∪ I c zl ∪ I l ci ∪ I cl+ zi ∪ I cl- zi           (27)
(compare (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4')).

Next, assume that (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) solves the first-order system at (γ, ū, λD , λν , λτ ), or equivalently,

H ′ ((γ, ū, λD , λν , λτ ); (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ )) = 0. (28)
If we want to find a selection function H (j) with

∇H (j) (γ, ū, λD , λν , λτ )(γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) = H ′ ((γ, ū, λD , λν , λτ ); (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ )) = 0, (29)
one can see from (P ′ 3) and (P ′ 4) that it suffices to choose the determining index sets so that the following conditions are fulfilled:

λ ′ ν,i -(B ν u ′ ) i ≥ 0, i ∈ I f z , λ ′ ν,i -(B ν u ′ ) i ≤ 0, i ∈ I c z , s i (λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i ) ≤ -F i (λ ′ ν,i -(B ν u ′ ) i ), i ∈ I s ci , |λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i | ≤ -F i (λ ′ ν,i -(B ν u ′ ) i ), i ∈ I cs zi , s i (λ ′ τ,i -(C τ u ′ ) i + γ ′ V τ,i ) ≥ -F i (λ ′ ν,i -(B ν u ′ ) i ), i ∈ I l ci ∪ I cl+ zi ∪ I cl- zi ,                (30)
with s i introduced by ( 26) and s i = sgn( λτ,i ) for i ∈ I s ci . We are prepared to state the continuation result based on Theorem 1.

Theorem 2. Let Assumption 1 be satisfied, (γ, ū, λD , λν , λτ ) be a solution of (P) with non-empty I z ∪ (I c ∩ I i ), and (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) solve the first-order system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4'). If one can find index sets I c z , I l ci , I cl+ zi , I cl- zi , I f z , I s ci , I cs zi and I c zl in accordance with (24) and (25) and such that (i) all the inequalities in (30) are satisfied in strict sense;

(ii) the gradient of the corresponding selection function H (j) at (γ, ū, λD , λν , λτ ) given by (27) has the maximal rank, then there exist δ > 0 and a C 1 -curve c : [s, s + δ) → R 1+2(n Ω +n D +nc) such that Moreover, the image of any other curve c : [s, s + δ) → R 1+2(n Ω +n D +nc) with δ > 0 and such that ( j) c(s) = (γ, ū, λD , λν , λτ );

( jj) H(c(s)) = 0, ∀s ∈ [s, s + δ);

( jjj) c′ (s) = 0 satisfies (30)

       (31)
coincides with the image of c around (γ, ū, λD , λν , λτ ).

Proof. Let (γ, ū, λD , λν , λτ ) and (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) meet the imposed assumptions and H (j) be the selection function corresponding to the particular index sets from the assertion. We shall verify the assumptions of Theorem 1 for i 0 = j.

First, the construction of the selection functions above guarantees that all selection functions corresponding to various combinations of the index sets restricted by ( 24) and ( 25) are active at (γ, ū, λD , λν , λτ ). In particular, H (j) (γ, ū, λD , λν , λτ ) = H(γ, ū, λD , λν , λτ ) = 0.

Second, since (30) holds for (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ), ( 29) is valid for our selection function H (j) . Finally, if a selection function H ( j) corresponds to Ĩc z , Ĩl ci , Ĩcl+ zi , Ĩclzi , Ĩf z , Ĩs ci , Ĩcs zi and Ĩc zl that fulfil the same conditions as in ( 24) and ( 25) and at least one of the sets Ĩc z , Ĩl ci , Ĩcl+ zi , Ĩclzi differs from I c z , I l ci , I cl+ zi or I cl- zi , respectively, then we claim that ∇H ( j) (γ, ū, λD , λν , λτ )(γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) = 0. Indeed, making use of (29) and the assumption (i), one gets either

i ∈ Ĩc z ∩ I f z =⇒ (B ν u ′ ) i = -λ ′ ν,i + (B ν u ′ ) i < 0 or i ∈ Ĩf z ∩ I c z =⇒ λ ′ ν,i = λ ′ ν,i -(B ν u ′ ) i < 0 if Ĩc z = I c z .
Similarly, assuming that Ĩc z = I c z but another couple of the index sets does not coincide, one can show that

i ∈ Ĩl ci ∩ I s ci =⇒ si λ ′ τ,i + F i (λ ′ ν,i -(B ν u ′ ) i ) < 0, (i ∈ Ĩs ci ∩ I l ci ) ∨ (i ∈ Ĩcs zi ∩ (I cl+ zi ∪ I cl- zi )) =⇒ s i ((C τ u ′ ) i -γ ′ V τ,i ) < 0, (i ∈ Ĩcl+ zi ∩ I cl- zi ) ∨ (i ∈ Ĩcl- zi ∩ I cl+ zi ) ∨ (i ∈ ( Ĩcl+ zi ∪ Ĩcl- zi ) ∩ I cs zi ) =⇒ λ ′ τ,i + si F i (λ ′ ν,i -(B ν u ′ ) i ) = 0,
(compare (P ′ 3') and (P ′ 4')). This means that the direction (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) satisfying the assumption (i) in Theorem 2 leads to solutions with strict contact modes of all nodes, without any grazing contact or impending slip.

Similarly, non-strict inequalities in (30) correspond to (33) with non-strict inequalities above. Thus, the uniqueness part of Theorem 2 means that there is only one branch of solutions with the nodal contact modes corresponding to (30) according to the interpretation above with non-strict modes. From this mechanical point of view, the claim of Theorem 2 resembles the one of [17, Proposition 8].

Proposition 5 (Satisfaction of the assumption (ii) in Theorem 2). Let (P) correspond to a model of linearised elasticity with a body properly constrained by the Dirichlet condition. If g satisfies Assumption 1 and ∇ 2 g / ∇g and F i , i ∈ (I c ∪ I z ) ∩ (I l ∪ I i ) are small, then the gradients of all selection functions H (j) given by ( 27 Suppose first that ∇ 2 g = 0 and F i = 0 for any i ∈ (I c ∪ I z ) ∩ (I l ∪ I i ). Then K is positive definite on ker(B D ), C τ = B τ , SF = 0 and we claim that the partial gradient ∇ (u,λ D ,λν ,λτ ) H (j) (γ, ū, λD , λν , λτ ) is non-singular for any H (j) . Indeed, one has ∇ (u,λ D ,λν ,λτ ) H (j) (γ, ū, λD , λν , λτ ) =

Figure 4 :

 4 Figure 4: Components of the cones for the indices i ∈ I z ∩ I i : (a) C x ; (b) C y and h(C y ) = (C x ) * .

Figure 5 :

 5 Figure 5: Geometry of Example 3.

  (j) c(s) = (γ, ū, λD , λν , λτ ); (jj) H(c(s)) = 0, ∀s ∈ [s, s + δ); (jjj) c ′ (s) = (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ).

  ) have the maximal rank for any P ∈ R 2n Ω and any U ′ D ∈ R 2n D . Proof. We shall argue similarly as in [3, proof of Proposition 3.2]. By (27), ∇H (j) (γ, ū, λD , λν , λτ ) = {1, . . . , 2n Ω }, F denotes jointly the indices from I f and I f z , C the indices from I c and I c z , N the indices from I f , I f z and I 0 , S the indices from I s , I s ci and I cs zi , and L the indices from (I c ∩ I l ), I c zl , I l ci , I cl+ zi and I cl- zi . Further, S = Diag(s 11 , . . . , s ncnc ), s ii = s i if i ∈ (I c ∩ I l ) ∪ I c zl ∪ I l ci ∪ I cl+ zi ∪ I cl- zi , 0 otherwise, F = Diag(F 1 , . . . , F nc ).

  Kτcl,N Kτcl,Z Kτcl,νzl Kτcl,τcl Kτcl,τzl Kτcl,νci Kτcl,τci Kτzl,u Kτzl,N Kτzl,Z Kτzl,νzl Kτzl,τcl Kτzl,τzl Kτzl,νci Kτzl,τci Kci,u Kci,N Kci,Z Kci,νzl Kci,τcl Kci,τzl Kci,νci Kci,τci

	u	Ku,N	Ku,Z	Ku,νzl	Ku,τcl	Ku,τzl	Ku,νci	Ku,τci
	KN,u KN,N KN,Z KN,νzl KN,τcl KN,τzl KN,νci KN,τci
	KZ,u	KZ,N	KZ,Z	KZ,νzl	KZ,τcl	KZ,τzl	KZ,νci	KZ,τci
	Kνzl,u Kνzl,N Kνzl,Z Kνzl,νzl Kνzl,τcl Kνzl,τzl Kνzl,νci Kνzl,τci
	Kτcl,u

It consists of a continuous branch of solutions for any γ ′ > 0 fixed.To sum up, we see that the first-order system can lack solutions for some γ ′ if F i are large for i ∈ I c ∩ I i , i ∈ I c ∩ I l , or i ∈ I z ∩ I l . On the other hand, it can possess multiple solutions for some γ ′ if F i are large for i ∈ I c ∩ I i , i ∈ I c ∩ I l , i ∈ I z ∩ I l or i ∈ I z ∩ I i .

Invoking that a symmetric positive definite matrix is a P -matrix and arguing as for Corollary 1, one gets the following:

Corollary 2 (Uniqueness). Let (P) correspond to a model of linearised elasticity where the body is properly constrained by the Dirichlet condition. If g satisfies Assumption 1 and ∇ 2 g / ∇g and F i , i ∈ (I c ∪ I z ) ∩ (I l ∪ I i ), are small, then the first-order system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4') has a unique solution for any γ ′ ∈ R, P ∈ R 2n Ω and any U ′ D ∈ R 2n D . Example 2 (Large coefficient of friction). (i) Let us consider the simple continuation problem (8) from Example 1 with F = a/b, (L 1 (γ), L 2 (γ)) = (1, -γ) and u k 2 = V τ = 0. Similarly as in [9], one can show that this problem has solutions corresponding to stick of the only contact node for γ ∈ [0, F ], which form the set

solutions corresponding to non-positive slip for γ = F forming the set

and solutions corresponding to separation for γ ≥ F , which form the set

Taking (γ, ū1 , ū2 , λν , λτ ) = (F , 0, 0, -1, F ) (strong contact with impending slip), we get the first-order system:

that is, there is no solution for γ ′ > 0 but infinity of solutions for γ ′ = 0. The first-order system for (γ, ū1 , ū2 , λν , λτ ) = (F , 0, β, -1bβ, F (1 + bβ)), β ∈ (-1/b, 0) (strong contact with positive slip) consists of ( 22) and

It has no solution for γ ′ = 0 but infinity of solutions for γ ′ = 0.

Finally, the first-order system for (γ, ū1 , ū2 , λν , λτ ) = (F , 0, -1/b, 0, 0) (grazing contact with positive slip) is composed of ( 22) and

There is no solution for γ ′ < 0 but infinity of solutions for γ ′ = 0.

(ii) Now, let us consider the problem (8) with F = a/b, (L 1 (γ), L 2 (γ)) = (bγ, -aγ) and u k 2 = V τ = 0 and choose (γ, ū1 , ū2 , λν , λτ ) = (0, 0, 0, 0, 0) (grazing contact with impending slip). In this case, the first-order system comprises ( 22) with (P 1 , P 2 ) = (b, -a) and

. Trying all possible cases, one can find the following solution set for γ ′ ≥ 0:

As in the previous example, the case (u

) and r is small enough so that a -bF -(1 + F 2 )/r < 0, one obtains solutions forming the set

Hence, this first-order system has no solution for γ ′ > 0 and two different solutions for γ ′ < 0. This shows that large ∇ 2 g / ∇g can cause both lack and multiplicity of solutions for some γ ′ .

Let us point out that the first-order system can lack solutions for some γ ′ for larger coefficients of friction even in the case of linearised elasticity with the body properly constrained by the Dirichlet condition and in contact with a flat obstacle. This is in stark contrast to the corresponding incremental problem, where existence of a solution is guaranteed for any coefficient of friction [14, Section 6]. Nevertheless, it is in accordance with the analysis of the rate and quasi-static problem. It was shown in [17,Subsection 5.2] that solvability of the rate problem requires smallness of the coefficient of friction for nodes in strong contact with possible slip in near future and an example of non-existence of solutions of a quasi-static problem for high values of the coefficient of friction can be found in [13].

To add, the uniqueness result for the rate problem presented in [17] involves an assumption on smallness of the coefficient of friction for nodes in contact with possible slip in near future and for nodes in grazing contact. In [13], an example with multiple solutions of a quasi-static problem with high values of the coefficient of friction was shown.

Continuation in a Direction Solving the First-Order System

To show the possibility of continuation of a solution curve in a direction solving the firstorder system, we shall apply Theorem 2 from the companion paper [16]. First of all, we restate the theorem, taking into account Remark 2 from the same paper.

functions of H at (γ, x) and there exists i 0 ∈ I such that

where s i are defined by ( 26) and

Therefore, the assumptions of Theorem 1 hold and existence of the announced curve c follows.

Now, let c be a curve from the second part of the assertion. Comparing ( 23) and ( 31), we have to show (23)( jjj):

From ( 31)( jj), we see that (H • c) ′ (s) = 0. In addition, we deduce from Proposition 2 that (H • c) ′ (s) = H ′ (c(s); c′ (s)). Taking into account (31)( j) and ( jjj), we thus have similarly as before that ∇H (j) (γ, ū, λD , λν , λτ )c ′ (s) = H ′ (c(s); c′ (s)) = 0. This combined with (29) and the assumption (ii) implies that

) satisfies all the inequalities in ( 30) strictly, so it is clear that these conditions are violated for any r(γ

, which is assumed to be non-empty. This together with (31)( jjj) lead to (32) and the uniqueness part of Theorem 1 applies.

In the trivial case when I z ∪(I c ∩I i ) = ∅, H is smooth in some vicinity of (γ, ū, λD , λν , λτ ) and the classical implicit function theorem can be employed.

Theorem 3. Let Assumption 1 hold, (γ, ū, λD , λν , λτ ) be a solution of (P) with

be a non-zero vector solving the first-order system (P ′ 1), (P ′ 2), (P ′ 3') and (P ′ 4'). If ∇H(γ, ū, λD , λν , λτ ) has the maximal rank, then there exist δ > 0 and a

Moreover, the image of any other solution curve passing through (γ, ū, λD , λν , λτ ) is in some neighbourhood of (γ, ū, λD , λν , λτ ) contained in the image of c.

Remark 4. Assume that the vector (γ ′ , u ′ , λ ′ D , λ ′ ν , λ ′ τ ) satisfies (28). Then one can easily verify that (30) holds strictly if and only if

which is obviously a non-singular matrix iff the matrix

) = 0, one gets from the last three blocks of lines that

This together with the first block of lines of the system gives 0

and positive definiteness of K on ker(B D ) implies that u ′ = 0. Consequently,

and taking into account Remark 3, we have necessarily λ ′ D = 0, λ ′ νC = 0 and λ ′ τ S = 0 due to full row rank of

. This shows that J is non-singular and the conclusion is valid in the special case.

Continuous dependence on the elements of ∇ 2 g/ ∇g and F i completes the assertion.

Numerical Experiments

This section presents results computed with the predictor-corrector continuation method described in [16, Section 4] and performed with the finite-element library GetFEM++ [19].

V-Shaped Body

The first experiment is done for a V-shaped body from [17, Subsection 8.2], which is depicted in Figure 7(a). It is approximately 0.82 m wide, 0.8 m high and we consider a plane-strain model of linearised elasticity with the Lamé constants λ = 100 GN/m 2 and µ = 82 GN/m 2 . The body is subject to volume forces of the density f (γ) = (-5 GN/m 3 , 1-8(1γ) GN/m 3 ) and it is fixed along both parts of Γ D and free of surface tractions on Γ N , that is, U D = 0 and h = 0. The continuation problem is taken with w = 0 and V τ = 0. Three unstructured meshes with 2220, 8880 and 35520 linear triangles (M2220, M8880 and M35520) are used for the finite-element discretisation. The numerical continuation is initialised with ε = ε ′ = 5 • 10 -12 , c min = 0.99999, j thr = 4, j max = 5, h dec = 0.5, h inc = 1.3, h min = 5 • 10 -7 , h max = 0.1 and after a change of the tangent performed by Algorithm 2 from [16], the predictor-corrector is restarted with h = 5 • 10 -4 .


Jean Grégoire

Silvère Bonnabel

Arnaud De

La Fortelle

Optimal cooperative motion planning for vehicles at intersections

We consider the problem of cooperative intersection management. It arises in automated transportation systems for people or goods but also in multi-robots environment. Therefore many solutions have been proposed to avoid collisions. The main problem is to determine collision-free but also deadlockfree and optimal algorithms. Even with a simple definition of optimality, finding a global optimum is a problem of high complexity, especially for open systems involving a large and varying number of vehicles. This paper advocates the use of a mathematical framework based on a decomposition of the problem into a continuous optimization part and a scheduling problem. The paper emphasizes connections between the usual notion of vehicle priority and an abstract formulation of the scheduling problem in the coordination space. A constructive locally optimal algorithm is proposed. More generally, this work opens up for new computationally efficient cooperative motion planning algorithms.

I. INTRODUCTION

In this paper, we consider the problem of designing an intelligent transportation system meant to improve the safety at intersections, and to decrease the average time spent by vehicles going through the intersection. Due to the achievements and promises in autonomous cars design, this subject has attracted much interest, and a variety of systems have been proposed to address this issue [START_REF] Dresner | Multiagent traffic management: a reservationbased intersection control mechanism[END_REF], [START_REF] Hafner | Automated vehicleto-vehicle collision avoidance at intersections[END_REF], [START_REF] Mehani | Trajectory planning in a crossroads for a fleet of driverless vehicles[END_REF], [START_REF] Fraichard | Planning movements for several coordinated vehicles. In Intelligent Robots and Systems '89. The Autonomous Mobile Robots and Its Applications[END_REF], [START_REF] Colombo | Efficient algorithms for collision avoidance at intersections[END_REF]. The clear improvements in the intersection flow compared to the use of traffic lights [START_REF] Dresner | Multiagent traffic management: a reservationbased intersection control mechanism[END_REF], [START_REF] Kowshik | Provable systemwide safety in intelligent intersections[END_REF], and in terms of safety [START_REF] Dresner | Mitigating catastrophic failure at intersections of autonomous vehicles[END_REF], have become strong drivers for research. Moreover, conflict resolution is also an intensive field of research in railway and air transportation systems for instance [START_REF] Ismail | Railway traffic control and train scheduling based oninter-train conflict management[END_REF], [START_REF] Tomlin | Conflict resolution for air traffic management: a study in multiagent hybrid systems[END_REF].

As the task has a potentially very high complexity (essentially one tries to compactify the vehicles trajectories in space-time), the problem is often decomposed into two parts. The first part consists of determining fixed paths along which vehicles cross the intersection. The second part consists of adapting the velocity of each vehicle along the path in order to avoid collisions -generally resulting from human mistakes -and to optimize as well traffic by increasing flows, i.e. by decreasing the average time spent in the intersection. This approach was initiated in [START_REF] Kant | Toward efficient trajectory planning: The path-velocity decomposition[END_REF] for multi-robot applications, and has become standard for vehicle management at intersections, e.g., [START_REF] Fraichard | Planning movements for several coordinated vehicles. In Intelligent Robots and Systems '89. The Autonomous Mobile Robots and Its Applications[END_REF], [START_REF] Hafner | Automated vehicleto-vehicle collision avoidance at intersections[END_REF], [START_REF] Leroy | Multiple path coordination for mobile robots: A geometric algorithm[END_REF], [START_REF] Akella | Coordinating the motions of multiple robots with specified trajectories[END_REF].

The problem of planning velocities of robots along fixed paths has attracted much attention over the past decades, and was essentially motivated by applications for robot manipulators and automated guided vehicles in factories. A usual approach in motion planning consists of transferring a Mines ParisTech, Centre de Robotique, Mathématiques et Systèmes, 60 Bd St Michel 75272 Paris Cedex 06, France † Inria Paris -Rocquencourt, IMARA team, Domaine de Voluceau -Rocquencourt, B.P. 105 -78153 Le Chesnay, France planning problem in the physical space to a planning problem in an abstract space: the configuration space. For the problem of planning velocities along fixed paths, the configuration space is called coordination space, first introduced in [START_REF] O'donnell | Deadlock-free and collision-free coordination of two robot manipulators[END_REF] for robot manipulators, and possesses a specific cylindrical structure, as first noticed by [START_REF] Lavalle | Optimal motion planning for multiple robots having independent goals[END_REF]. In this paper, we propose to use the coordination space approach as a rigorous mathematical framework particularly well-suited to the problem of motion planning at intersections. We suppose the paths of the intersection fixed and we only focus on the coordination problem.

The process is then inherently cooperative: some vehicles need to brake to let some others go ahead until the conflict configuration is exited. The purpose of the present paper is to cast the problem in a well-defined mathematical framework, and to provide a provably collision-free solution that minimizes the average exit time over the vehicles. To do so we revisit the notion of priority introduced in [START_REF] Buckley | Fast motion planning for multiple moving robots[END_REF] and [START_REF] Fraichard | Planning movements for several coordinated vehicles. In Intelligent Robots and Systems '89. The Autonomous Mobile Robots and Its Applications[END_REF]. It appears to be especially meaningful for the problem considered as we recover the usual and intuitive notion of vehicle priority at intersections. We show it is a powerful tool to classify the locally optimal trajectories that have been studied as special representatives of homotopy classes in [START_REF] Ghrist | Nonpositive curvature and pareto optimal coordination of robots[END_REF]. We introduce a priority graph and prove several properties of this graph with respect to our problem. We use it to build an algorithm for the design of collision-free and deadlock-free optimal trajectories. The paper has two main contributions. From an applicative viewpoint, it shows the coordination space is a powerful mathematical tool that is especially suited to the problem of coordinating vehicles at intersections. In particular, it opens a new research avenue for the design of efficient algorithms. From a theoretical viewpoint, our contribution consists of drawing links between a discrete scheduling problem (priority graph) and a continuous problem formulated in the abstract coordination space. It illustrates those links in an intuitive way and relates those notions to the existence of collisions or deadlocks and to the optimality problem. New results are derived (see in particular Propositions 7 and 9) and a constructive optimal algorithm is built upon the proposed methodology.

Section II formulates the problem and introduces the fixedpath approach. Section III defines a priority relation between vehicles underlying the scheduling problem and presents some new results. In section IV we present algorithms for building an optimal motion planning based on schedule decoupling.

II. PROBLEM FORMULATION

A. Modeling assumptions

The main assumption we make is to constrain the vehicles to follow predefined paths to go through the intersection. This assumption fits well intersections in a road network which are a highly constrained environment. Thus, every vehicle i follows a particular path γ i and we denote s i ∈ R its curvilinear coordinate on the path. The configuration of the system of vehicles is s = (s i ) i∈{1...n} and we denote φ(t) the evolution of s through time. The curvilinear coordinates are normalized and we let χ = [0, 1] n be the coordination space, where n is the number of vehicles going through the intersection (possibly changing through time) and {e i } 1≤i≤n the canonical basis of χ. The boundedness condition on χ is rather technical but ensures the whole intersection lies in a bounded region of R n (somehow interactions are limited to a bounded area).

To perform motion planning optimization, we need to make hypothesis on the technical constraints of the vehicles. These include kinematic constraints (maximum velocity, maximum curve radius, etc.) and dynamic constraints (limited acceleration, adherence, jerk, etc.). In this paper, as a first step, we assume that the velocity of every vehicle on its path is limited to v max = 1 without accounting for dynamic constraints. Moreover we assume the velocities are always positive (the cars do not go in reverse in the intersection). These assumptions are standard [START_REF] Fraichard | Planning movements for several coordinated vehicles. In Intelligent Robots and Systems '89. The Autonomous Mobile Robots and Its Applications[END_REF].

As every vehicle occupies a non-empty geometric region, some states must be excluded to avoid collisions between vehicles. We define the obstacle region χ obs as the open set of all collision states. A collision occurs when two vehicles occupy a same region of space, so that the obstacle region can be described as the union of n(n -1)/2 open cylinders χ ij obs corresponding to as many collision pairs: χ obs = ∪ i>j χ ij obs [START_REF] Lavalle | Planning Algorithms[END_REF]. Indeed obviously χ ij obs = χ ji obs and χ ii obs = ∅. Figure 2 displays the obstacle region and a collision-free path for a two-path intersection.

Moreover, we assume that each cylinder χ ij obs has an open bounded convex cross-section (in the plane generated by e i and e j ). This assumption implies cross-sections are simply connected and excludes some cases of real intersections such as a pair of opposite turn-left, but the results we present in this paper could be easily extended with a more relaxed assumption. Cross-sections are open so that the complementary set is closed and hence complete. It also ensures all cross-sections are included in the interior of [0, 1] 2 : no collision can occur for a vehicle at coordinates 0 or 1. Such assumptions are very realistic, except for the case of vehicles following each other on a same path, but our framework can be extended to this case, as we plan to prove in future work.

B. Formulation in the coordination space

Let χ free = χ\χ obs denote the obstacle-free space and 1 = (1 • • • 1). The initial condition x init belongs to χ free , and the goal region is χ goal = {1} ⊂ χ free . In this paper we will consider continuous paths φ : [0, T ] → χ which are piecewise right differentiable. A collision-free path will be such that φ(t) ∈ χ free for all t.

Definition 1 (Feasible paths). A feasible path for the considered problem is a right-differentiable collision-free path such that

• φ(0) = x init and φ(T ) ∈ χ goal • φ (t) ≥ 0 and φ (t) ∞ ≤ 1 ∀t ∈ [0, T ].
In order to define an optimality criterion for the motion planning problem, we need to define a cost function c(φ) enabling to compare feasible paths. As explained in the introduction, the purpose is to minimize the average exit time over the vehicles. As a consequence, we define the cost function c as:

c(φ) = 1 n n i=1 T i = 1 n n i=1 φ -1 i (1)
with φ -1 i (s i ) denoting the first date at which φ i reaches s i and T i = φ -1 i (1) being the exit time for the vehicle i. The optimality problem consists of finding a feasible path φ * that minimizes the cost, i.e. for any feasible path φ, we have c(φ * ) ≤ c(φ). Figure 1 depicts three feasible trajectories with different costs in a two-dimensional coordination space for a given initial condition. The left drawing depicts two feasible trajectories in a twodimensional coordination space. With the trajectory φ, the vehicle 2 leaves the intersection earlier than with φ, while exit time is not changed for the vehicle 1. Hence, c( φ) < c(φ). φ * is an optimal trajectory since every vehicle is always at maximum speed: φ * = arg min c(φ).

Note that the cost can be easily bounded from above and from below. Suppose for simplicity's sake that x init = 0 and χ goal = {1}. If χ obs = ∅ the optimal path is given by φ ≡ 1. The corresponding cost c(φ) = 1 n n i=1 1 = 1 is thus a lower bound on the optimal time in the general case. On the other hand, a suboptimal trajectory is easily found moving the vehicles one after the other, i.e. moving along the edges of the hypercube. This trajectory is collision-free, because by assumption a vehicle cannot collide another one if its coordinate is 0 or 1.

It yields c(φ) = 1 n n i=1 i = n+1
2 . Note that vehicle i is free after T i when its coordinate is 1: there is no further interaction with any other vehicle. We hereafter take φ i (t) = 1 for t ≥ T i in order to stay in χ = [0, 1] n but any further trajectory φ i (t) ≥ 1 would be suitable had we taken χ = R n .

C. Geometry of the coordination space

In this section, we propose several geometric concepts in the coordination diagram.

1) Completing the obstacle region: For each collision pair (i, j) we define the south of χ obs along the plane e i × e j as:

S ij (χ obs ) = χ ij obs -R + e j (1) = s -λe j : ∀s ∈ χ ij obs , ∀λ ∈ R +
We say the south of χ ij obs along the plane e i × e j is the west of χ ij obs along the plane e j × e i and denote W ij (χ obs ) = S ji (χ obs ). These sets are clearly open. They play the role of gates: a feasible trajectory must go through one or the other and they are exclusive. As a result, the intersection of the south and the west is forbidden for feasible trajectories, even if there is no collision. As the velocities must be positive, being in the intersection at some time would necessarily lead to a future collision unless the cars stop moving (see Figure 3). This situation is referred to as a deadlock. To avoid such situations, an extended obstacle region is defined, so that no basic one-to-one dead-lock situation can occur [START_REF] O'donnell | Deadlock-free and collision-free coordination of two robot manipulators[END_REF].

Definition 2 (SW-Completion). The South-West completion of χ obs along the plane e i × e j is given by

SW ij (χ obs ) = S ij (χ obs ) ∩ W ij (χ obs )
where S ij (χ obs ) is defined in [START_REF] Akella | Coordinating the motions of multiple robots with specified trajectories[END_REF]. The completed obstacle region χ obs is defined by

χ obs = i =j SW ij (χ obs )
We call χ free = χ\χ obs the obstacle-free completed region.

The South-West closure is applied on each cylinder χ ij obs as depicted on Figure 3. An important fact -not explicitly mentioned in the literature yet -is that this method cannot be readily extended to avoid dead-locks in multi-dimensional planning problems while ensuring optimality. For this reason, the next sections rather present tools based on gates' positions to avoid dead-locks caused by multiple vehicles' interaction.

2) Gates in the coordination space: Let define the region H i j as follows:

H i j = S ij (χ obs ) \ SW ij (χ obs )
The top drawings of Figure 2 display H i j and H j i . Over the convexity hypothesis of the cylinders of χ obs , we can assert that any feasible trajectory will necessarily cross H i j or H j i exclusively. This is why these two boxes can be viewed as gates in the coordination space. In Figure 2, the path φ i j intersects the gate H i j : the vehicle j slows down to let the vehicle i go ahead. Each gate favours a vehicle over another i.e. a priority is assigned, and passing through a gate possibly prevents from passing through many other gates because of the positive speed assumption. H i>j SW ij (χ obs ) Fig. 2. The top-left drawing represents in the plane (s i , s j ) the obstacle region χ obs and a feasible trajectory φ. The top-right drawing visualizes the two gates H i j and H j i . For each gate, a feasible trajectory passing through the gate is given. The bottom drawings depict χ i j obs and χ j i obs .

We propose to define the obstacle region χ i j obs as a set that is incompatible with the crossing of gate H i j , as illustrated on Figure 2:

χ i j obs = χ ij obs -R + e i + R + e j

III. PRIORITY RELATIONS AND ORDERING

The geometry of the coordination space leads us to define a natural binary relation corresponding to priority relations between vehicles: a very familiar and intuitive notion in real life. Indeed, we say the vehicle i has priority over the vehicle j in the intersection if the associated path goes through the gate H i j in the coordination space.

Definition 3 (Priority relation). A feasible path φ induces a binary relation

on the set {1...n} as follows. For i = j, i j if φ intersects H i j .

The following proposition proves a result that is not totally obvious from an intuitive viewpoint. In the right drawing, the red vehicle goes ahead the blue vehicle and the green vehicle goes ahead the red vehicle, while the blue vehicle goes ahead the green vehicle. It yields in a non transitive binary relation: the set {1...n} is not partially ordered.

A. Priority graph

In this subsection, we propose to define the priority relations in terms of a graph. As for any collision pair i, j we have either i j or j i, the priority relation can be defined by an oriented graph G whose vertices are {1...n}. Definition 5. Given a feasible trajectory φ, we call the priority graph the oriented graph G whose vertices are {1...n} and such that i G -→ j if and only if i j.

Note the graph is complete in the sense that for all i = j there is one, and only one, arc linking i and j. There is no arc linking i to itself. Generally speaking, there are potentially 2 n(n-1) 2 possible priority graphs. Actually, one can easily understand that some priorities are not acceptable. The next section deals with this issue and characterizes what we call feasible priorities.

B. Feasible priorities

Definition 6 (Feasible Priority). A priority graph G is feasible if there exists a feasible trajectory φ whose priority graph is G.

In the intersection depicted on the left drawing of Figure 5, we can assert that some priorities are not feasible. Indeed, if the the vehicle 2 goes ahead the vehicle 1 and the vehicle 2 goes ahead the vehicle 3, the vehicle 1 must go ahead the vehicle 3, because the three paths have a common intersection point. Hence, for this particular situation, the priority graph cannot be cyclic and the priority relation is an order relation. The proposition below provides a necessary condition for priorities feasibility based on the relative position of the fixed-priority collision cylinders χ i j obs . The proof is omitted due to space limitations.

Proposition 7 (Characterization of feasible priority graphs). A priority graph G is feasible if and only if for any cycle

C in G, i C -→j χ i j obs = ∅
As χ ij obs ⊂ χ i j obs it readily implies: Proposition 8 (Necessary condition for priorities feasibility). Suppose G is a feasible graph. Then, for any cycle C in G,

i C -→j χ ij obs = ∅
In the intersection of Figure 5, the state (s i = 1 2 ) i=1,2,3 corresponding to a common intersection point is in each of the three collision cylinders. Hence, the necessary condition applies and it can be asserted that no cycle between the 3vehicle priority graph can exist. The underlying reason in the coordination space is that if there is a cycle in the priorities, then the trajectory must stay in the box [0, 1 2 ] 3 at any time to respect these priorities, i.e. cannot reach χ goal .

The characterization of feasible priorities shows that the feasibility of a priority graph relies on the relative position of the collision cylinders in the coordination space.

IV. THE MOTION PLANNING PROBLEM

The last section shows that the set of feasible trajectories can be divided into as many disjoint subsets as there are feasible priority graphs. In this section, we shall present an algorithm that builds an optimal trajectory for a fixed priority graph. We will then discuss the computational issues for determining the globally optimal trajectory.

A. Optimal trajectories with fixed priority graph

The main idea is that when priorities are fixed, if at any time every vehicle runs as much distance as possible in a time step, then we can assert that the trajectory is optimal (i.e. minimizes the average time spent in the intersection) for the fixed priorities. It means the optimal solution can be defined reactively.

Let H G denote the union of gates corresponding to the priorities fixed by G. After setting φ(0) = x init , the algorithm proceeds as follows. If G is not a feasible priority graph, the algorithm stops, else, until χ goal is reached, at each date t, the trajectory is defined via φ (t) which is calculated as a function of φ(t) by the following algorithm, that essentially consists of projecting the vector 1 consecutively: Given a feasible priority graph G, the trajectory φ * (t) generated by the algorithm described above is an optimal trajectory for the given priorities.

• φ (t) ← 1 • Stay in χ: For i such that φ i (t) = 1, do φ i (t) ← 0 • Avoid
It can be noted that the vehicles may brake for two reasons: avoiding the obstacle region or avoiding a gate that is not compatible with the priorities defined by G. Figure 6 displays a trajectory generated by the algorithm for a 3-paths intersection. For simplicity's sake, for each collision pair, the one-to-one obstacle region is approximated to a square.

Proof. Suppose there exists a trajectory ψ respecting the priorities G and the initial condition ψ(0) = φ * (0) such that ψ i0 overtakes φ * i0 for the first time at the date t 0 > 0. Consider the trajectory defined by φ i (t) = max(ψ i , φ * i )(t). The trajectory φ admits G as priority graph, it is collisionfree (because ψ and φ * respect the same priorities, see Figure 2) and it respects the maximum velocity constraint. Moreover, φ i (t 0 ) = ψ i (t 0 ) and for small enough > 0, φ i0 (t 0 + ) ≥ φ * i0 (t 0 ). By construction of φ * , for small enough , φ i (t 0 + ) ≤ φ * i (t 0 ), so that the last two conditions require φ i0 (t 0 + ) = φ * i0 (t 0 + ) = ψ i0 (t 0 + ) to respect both inequalities. As a consequence, for all i ∈ {1...n}, the trajectory φ * i is never strictly overtaken by a vehicle in another feasible trajectory ψ i : φ * is an optimal trajectory. The optimal solution constructed by the algorithm is known as the "left-greedy" solution [START_REF] Ghrist | Nonpositive curvature and pareto optimal coordination of robots[END_REF]. In this latter paper the authors address a more general problem of optimization in coordination space, and regroup the trajectories in equivalence classes based on homotopy (two trajectories in the

3 1 2 s 1 s 2 κ ij s 2 s 3 s 1 κ ij s 3
Fig. 6. The bottom-left drawing represents a 3-path intersection where priorities are fixed to 1 2, 2 3 and 1 3. The three other drawings visualize in each plane (s i , s j ) the obstacle region (in blue), the forbidden gates (red), and the optimal trajectory generated by the algorithm. The color of the trajectory changes when a collision pair leaves or reaches a gate or the boundary of the obstacle region. same class can be deduced from each other by a continuous deformation). In our formalism, each admissible graph defines an homotopy class. In our particular application, we recover their optimality result with simple means (as opposed to their proof that is based on concepts of advanced differential geometry, namely Gromov's link condition for CAT(0) spaces).

B. Optimal trajectories

As the algorithm of section IV-A generates an optimal trajectory for a given priority graph, the motion planning optimization boils down to a combinatorial problem. By computing and comparing the costs of the optimal trajectories among the feasible priority graphs, we can conclude on the optimal planning solution by choosing the least costly trajectory. This method is suitable when there are few vehicles crossing the intersection, but its high algorithmic complexity (due to the comninatorics of priorities) can become a prohibitive numerical cost when the number of vehicle grows. This motivates the following heuristic algorithm.

C. A heuristic algorithm with unconstrained graph

The last section proves that once the priority graph is fixed, the optimum is reached following algorithms belonging to the "bug" family, emanating from the work of [START_REF] Lumelsky | Path-planning strategies for a point mobile automaton moving amidst unknown obstacles of arbitrary shape[END_REF]. Indeed, the optimal paths go at maximum speed until they reach the frontier of a forbidden region. Then, they follow the boundary as long as necessary. In the absence of fixed priority graph, we propose to have the same reactive approach, moving the vehicles as fast as possible. A graph is thus iteratively built, and each time a gate is reached, the retained priority must be compatible with the graph already built.

After setting φ(0) = x init and initializing G to an empty graph, the algorithm proceeds as follows until χ goal is reached:

• φ (t) ← 1

• Stay in χ: For i such that φ i (t) = 1, do φ i (t) ← 0 • Update G: While there exists a collision pair i, j whose relative priorities are not defined in G and such that φ(t) belongs to the boundary of H i j , do i G -→ j if it yields a feasible priority graph, else do j G -→ i

• Avoid Forbidden Gates: While there exists i G -→ j such that φ(t) belongs to the boundary of the gate H j i , do: φ (t) ← orthogonal projection of φ (t) on the boundary of H j i

• Avoid Obstacle Region: While there exists i G -→ j such that φ(t) belongs to the boundary of χ obs ij , If 1 points towards χ ij obs , do: φ (t) ← orthogonal projection of φ (t) on any tangential hyperplane to χ obs ij at φ(t) The trajectory is defined recursively, the algorithmic complexity remains relatively low, and yields a generated priority graph G a posteriori. Moreover, we have guarantees the built trajectory is optimal for this priority graph.

V. CONCLUSION

The mathematical framework we propose to address the cooperative motion planning problem at intersections is based on path-velocity decomposition. Hence, the problem boils down to searching an optimal path in a highdimensional coordination space. We define a priority relation between the vehicles that provides a tool to classify the feasible trajectories, and propose an algorithm to construct an optimal trajectory for given priorities.

As the search for a global optimum is inherently related to the determination of all optima given a graph, it is of high algorithmic complexity. As a result, a good heuristics for fixing priorities is at the core of a simple quasi-optimal solution, as suggested by the intuition.

In the future, we plan to find such heuristics, and to take into account more realistic assumptions: varying number of vehicles, more general coordination space, and above all dynamic constraints. The next step will be the integration of the selected algorithms in a real intersection crossed by numerous vehicles such as cooperative cybercars developed at INRIA.

  Fig. 1.The left drawing depicts two feasible trajectories in a twodimensional coordination space. With the trajectory φ, the vehicle 2 leaves the intersection earlier than with φ, while exit time is not changed for the vehicle 1. Hence, c( φ) < c(φ). φ * is an optimal trajectory since every vehicle is always at maximum speed: φ * = arg min c(φ).

Fig. 3 .

 3 Fig.3. Two representations of a dead-lock situation involving two circleshaped vehicles. The left drawing depicts the real intersection with its two orthogonal paths and the two circle-shaped vehicles. Both of them have moved forward too much, so that now they cannot exit the intersection without turning around. The right drawing represents what happens in the two-dimensional coordination space and highlights the South-West region that must be added to χ obs to avoid the one-to-one dead-lock.

Proposition 4 .

 4 The binary relation induced by a feasible trajectory φ does not necessarily defines an order.

Figure 4

 4 Figure4depicts a 3-paths intersection with 3 intersection points and as many vehicles going through the intersection. In the right drawing, the red vehicle goes ahead the blue vehicle and the green vehicle goes ahead the red vehicle, while the blue vehicle goes ahead the green vehicle. It yields in a non transitive binary relation: the set {1...n} is not partially ordered.

Fig. 4 .

 4 Fig.[START_REF] Dresner | Multiagent traffic management: a reservationbased intersection control mechanism[END_REF]. Two representations of priority relations. In each drawing relation is represented in two ways: as a complete oriented graph, where orientation yields the priority; and as trajectories over time, foreground being first, background later. The left drawing represents a relation that is an order (even a total order). The right drawing shows a relation that is not an order.

Fig. 5 .

 5 Fig. 5. The left drawing represents a 3-path intersection with a common intersection point and 3 vehicles going through the intersection. The right drawing depicts what happens in the coordination space when cyclic priorities are set. Any feasible trajectory respecting cyclic priorities should stay in the red box [0, 1 2 ] 3 , but reaching the goal requires to cross this box.

  Forbidden Gates: While there exists i G -→ j such that φ(t) belongs to the boundary of the gate H j i , do: φ (t) ← orthogonal projection of φ (t) on the boundary of H j i• Avoid Obstacle Region: While there exists i

			G -→ j such
	that φ(t) belongs to the boundary of χ obs	ij , if 1 points
	towards χ obs	ij , do:
	φ (t) ← orthogonal projection of φ (t) on any tangential
	hyperplane to χ obs	ij at φ(t)
	Proposition 9 (Optimal motion planning for fixed priorities).