Samir Adly
email: samir.adly@unilim.fr

Hadia Rammal
email: hadia.rammal@xlim.fr

A New Method for Solving Second-Order Cone Eigenvalue Complementarity Problems

Keywords: Lorentz cone, Second-order cone eigenvalue complementarity problem, Semismooth Newton method, Lattice Projection Method

In this paper, we study numerical methods for solving eigenvalue complementarity problems involving the product of second-order cones (or Lorentz cones). We reformulate such problem to find the roots of a semismooth function. An extension of the Lattice Projection Method (LPM) to solve the second-order cone eigenvalue complementarity problem is proposed. The LPM is compared to the semismooth Newton methods, associated to the Fischer-Burmeister and the natural residual functions. The performance profiles highlight the efficiency of the LPM. A globalization of these methods, based on the smoothing and regularization approaches, are discussed.

Introduction

The theory of Complementarity Problems plays an important role in nonlinear optimization. The complementarity appears, for example, in the optimality conditions for constrained nonlinear programming or variational inequalities when the set of constraints is a closed and convex cone. The particular case of linear complementarity problems was widely discussed in the literature from both the theoretical and numerical point of view. More recently, the subject of eigenvalue complementarity problems has become one of the most well-established discipline in nonlinear optimization. It consists, for a given square matrix, in finding a scalar (eigenvalue) and nonzero vector (eigenvector) satisfying a complementarity condition over a closed and convex cone. When this cone coincides with the non-negative orthant, the problem is called Pareto eigenvalue complementarity problem. In this case, the spectrum is finite and its cardinality grows exponentially with the size of the matrix (see [START_REF] Adly | A new method for solving eigenvalue complementarity problems[END_REF][START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF][START_REF] Seeger | On cardinality of Pareto spectra[END_REF] for more details).

A wide variety of applications in Sciences and Engineering requires the resolution of eigenvalue complementarity problems such as the dynamic analysis of structural mechanical systems, vibro-acoustic systems, electrical circuit simulation, fluid dynamic, contact problem in mechanics (see, for instance, [START_REF] Martins | Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and coulomb friction[END_REF][START_REF] Martins | Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction[END_REF][START_REF] Martins | Bifurcations and instabilities in frictional contact problems: theoretical relations, computational methods and numerical results[END_REF][START_REF] Martins | Computation of bifurcations and instabilities in some frictional contact problems[END_REF][START_REF] Martins | The directional instability problem in systems with frictional contacts[END_REF]). The development of robust and efficient algorithms for solving eigenvalue complementarity problems (EiCP, for short) has received more attention in recent years such as the spectral projected gradient algorithm [START_REF] Júdice | On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm[END_REF], the semismooth Newton method [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF], the scaling and projection algorithm [START_REF] Pindo Da Costa | Cone-constrained eigenvalue problems: theory and algorithms[END_REF] or the Path solver [START_REF] Dirkse | The PATH solver: a non-monotone stabilization scheme for mixed complementarity problems[END_REF].

More recently, a wide range of applications in engineering design, transportation science, game theory, and economic equilibrium, can be formulated as optimization problems involving second-order cone constrained [START_REF] Lobo | Applications of second-order cone programming[END_REF] (know also as the Lorentz cone). Equally important is the study of second-order cone complementarity problems (SOCCP, for short) since it contains a large class of problems such as, for instance, nonlinear complementarity problems and second-order cone programming problems. Fukushima et al. [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF], Chen et al. [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF][START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF] studied smooth and nonsmooth approaches. In [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF], it has been shown that the natural residual and the Fischer-Burmeister functions can be extended to the SOCCP by means of Jordan Algebra.

There are presently certain theoretical properties and a large variety of computational methods for solving SOCCP. We cite, for example, interior-points methods, smoothing methods, SQP-type methods ( [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF][START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF][START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF][START_REF] Bonnans | Perturbation analysis of second-order cone programming problems[END_REF][START_REF] Chen | Alternative proofs for some results of vector-valued functions associated with second-order cones[END_REF][START_REF] Chen | An unconstrained smooth minimization reformulation of the second-order cone complementarity problem[END_REF][START_REF] Hayashi | A combined smoothing and regularization method for monotone second-order cone complementarity problems[END_REF][START_REF] Kato | An SQP-type algorithm for nonlinear second-order cone programs[END_REF][START_REF] Sun | Strong semismoothness of the Fischer-Burmeister SDC and SOC complementarity functions[END_REF][START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF][START_REF] Yamashita | A primal-dual interior point method for nonlinear optimization over second order cones[END_REF]). Furthermore, the projection mapping onto second-order cone (SOC, for short) was studied by several authors. In [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF][START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF], it was proved that projection mapping onto SOC is a strongly semismooth function. Later, in 2005, Hayashi et al. [START_REF] Hayashi | A combined smoothing and regularization method for monotone second-order cone complementarity problems[END_REF] gave an explicit representation for the Jacobian of the projection onto SOC.

Motivated by these recent developments, we study throughout this paper, the second-order eigenvalue complementarity problems (SOCEiCP). We reformulate the SOCEiCP into a system of semismooth equations by using the second-order cone complementarity functions (SOCC-functions, for short), the most frequently used in the literature, namely, the min-function and the Fischer-Burmeister function. Our main goal in this paper is to introduce the extended Lattice Projection Method (LPM) for solving SOCEiCP. LPM has been first introduced in [START_REF] Adly | A new method for solving eigenvalue complementarity problems[END_REF] to solve the Pareto eigenvalue complementarity problems. The second-order cone eigenvalue complementarity problems is considered to be one of the most difficult problems to solve, and this is due to the structure of the Lorentz spectrum which may be continuous (not finite and not countable). Therefore, compute and detect all the eigenvalues of this problem are not an easy task. As far as the authors know, this kind of technique for solving SOCE-iCP has never been studied before. In addition, we study under which conditions the Jacobian matrix in the SNM algorithm (see Algorithm 1), at a solution, is nonsingular. LPM is then compared to the semismooth Newton methods: SNM min and SNM FB , by using the performance profiles [24,[START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] as a comparison tool. The numerical experiments highlight that the LPM solver is efficient for solving SOCEiCP. This paper is organized as follows. In Sect. 2, we recall some background material on Euclidean Jordan Algebras and review the spectral factorization associated with a Lorentz cone. Moreover, we state a number of preliminary results for the projection mapping onto a second-order cone. In Sect. 3, we reformulate SOCEiCP into a system of semismooth equations. We extend, in Sect. 4, LPM to solve the SOCEiCP and give its B-subdifferential. In Sect. 5, we study the nonsingularity conditions of the Jacobian matrix, defined in the SNM algorithm. Then, we compare the LPM, in Sect. 6, with SNM min using the min-function and SNM FB via the Fischer-Burmeister function. In Sect. 7, we discuss a globalization of the three methods studied in the last sections by constructing appropriate merit functions. Finally, some concluding remarks are presented in Sect. 8.

Preliminaries

In this section, we recall some background material and preliminary results on complementarity problems and Euclidean Jordan Algebras (see, e.g., [START_REF] Faraut | Analysis on Symmetric Cones[END_REF][START_REF] Koecher | The Minnesota Notes on Jordan Algebras and Their Applications[END_REF]) as well as some definitions and properties of the B-subdifferential and (strong) semismoothness, which will be used in next sections. We recall that the Pareto eigenvalue complementarity problem consists in finding a scalar λ ∈ R and a vector x ∈ R n \ {0} such that x ≥ 0, λx -Ax ≥ 0 and x, λx -Ax = 0, where A is real matrix of order n. The classical complementarity problem consists in finding z ∈ R n such that z ≥ 0, F(z) ≥ 0 and z, F(z) = 0, where F : R n → R n is a map. As mentioned in [START_REF] Júdice | On the asymmetric eigenvalue complementarity problem[END_REF], the EiCP can be formulated as the following variational inequality (VI):

Find x ∈ C = {x ∈ R n + : e, x = 1} such that F(x), y -x ≥ 0, ∀ y ∈ C,
where F : R n → R n is defined by

F(x) := x T Ax x T x I -A x,
I is the identity matrix and e is a column vector, whose components are all one. The second-order cone complementarity problem consists in finding

ζ ∈ R m satisfying f (ζ ) ∈ K, g(ζ ) ∈ K, f (ζ ), g(ζ ) = 0, (1)
where f, g : R m → R n are continuously differentiable mapping and K is the Cartesian product of second-order cones in R n , defined by

K := K n 1 × K n 2 × • • • × K n r ,
and

n = n 1 + n 2 + • • • + n r . K n i ⊂ R n i
is the n i -dimensional second-order cone, also called the Lorentz cone or ice cream cone, defined by

K n i := {x = (x 1i , x 2i ) ∈ R × R n i -1 : x 1i ≥ x 2i }, (2)
where • denotes the Euclidean norm defined by x := √ x T x for a vector x ∈ R n . The cone K is self-dual, i.e., K + = K, where K + is the dual cone of K, defined by

K + := {x ∈ R n : x, y ≥ 0, ∀ y ∈ K}.
Corresponding to the Cartesian structure of K, we write in [START_REF] Adly | A new method for solving eigenvalue complementarity problems[END_REF]

f = ( f 1 , . . . , f r ) and g = (g 1 , . . . , g r ) with f i , g i : R m → R n i .
The second-order cone eigenvalue complementarity problem consists in finding a scalar λ > 0 and a vector x ∈ R n \ {0} satisfying x ∈ K, λx -Ax ∈ K and x, λx -Ax = 0.

(

) 3
The scalar λ is called the Lorentz eigenvalue of the real matrix A ∈ M n (R). From (1), it is sufficient to take

m = n + 1, ζ = (x, λ), f (ζ ) = f (x, λ) = x and g(ζ ) = g(x, λ) = λx -Ax.
Then, (1) reduces to (3). For two n-dimensional vectors x = (x 1 , x 2 ) and y = (y 1 , y 2 ) ∈ R × R n-1 , we define their Jordan product as

x • y := (x T y, x 1 y 2 + y 1 x 2 ). ( 4
)
An important character of Jordan Algebra is its eigen-decomposition, also called the spectral factorization with respect to the second-order cone K n . For more details, see, e.g., [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF][START_REF] Faraut | Analysis on Symmetric Cones[END_REF].

For any vector x = (x 1 , x 2 ) ∈ R × R n-1 , its spectral factorization is defined as

x := λ 1 (x)u (1) (x) + λ 2 (x)u (2) (x), (5)
where λ 1 , λ 2 ∈ R and u (1) , u (2) ∈ R n are the spectral values and the associated spectral vectors of x, respectively, given by

λ i (x) = x 1 + (-1) i x 2 , i = 1, 2, (6)
and

u (i) (x) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 2 1, (-1) i x 2 x 2 , if x 2 = 0, i = 1, 2, 1 2 (1, (-1) i w), if x 2 = 0, (7)
with w ∈ R n-1 , such that w = 1.

For the proof of ( 5), ( 6) and [START_REF] Martins | Computation of bifurcations and instabilities in some frictional contact problems[END_REF], see, e.g., [START_REF] Alizadeh | Second-order cone programming[END_REF].

We recall now some definitions and properties of the B-subdifferential and (strong) semismoothness, which will be used later.

Let H : R n → R m be a locally Lipschitz function. Then, the set

∂ B H (z) := M ∈ R m×n : ∃ (z k ) ⊆ D H : z k → z and lim k→+∞ ∇ H (z k ) = M ,
where D H := {z ∈ R n : H is differentiable at z}, is nonempty and is called the Bsubdifferential of H . The Clarke subdifferential [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] of H is defined by the convex hull

∂ H (z) := conv ∂ B H (z).
A vector z is called a B D-regular for H iff all matrices in the B-subdifferential ∂ B H (z) are non-singular. By using the subdifferential and the directional derivative (denoted by H (z; d)), we define semismoothness and strong semismoothness, which are first introduced by Mifflin [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF] and extended later by Qi and Sun [START_REF] Qi | A nonsmooth version of Newton's method[END_REF].

A directionally differentiable and locally Lipschitz function

H : R n → R m is said to be semismooth at x, iff Md -H (x; d) = o(||d||), for all M ∈ ∂ H (z + d) and d → 0. Moreover, if o(||d||) can be replaced by O(||d|| 2 ),
then the function H is said to be strongly semismooth.

Next, we recall some properties of the projection of any z ∈ R n onto a second-order cone.

Projection onto (SOC)

For z ∈ R n , the projection of z onto the second-order cone K n is defined by

P K n (z) := arg min z ∈K n z -z . ( 8
)
It can be written explicitly as

P K n (z) = max{0, λ 1 }u 1 + max{0, λ 2 }u 2 , ( 9
)
where λ 1 , λ 2 and u 1 , u 2 are defined by ( 5) and ( 6) [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF]Proposition 3.3]. Hence, we get

P K n (z) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ z (λ 1 ≥ 0, λ 2 ≥ 0), 1 2 (z 1 + z 2 ) 1, z 2 z 2 (λ 1 < 0, λ 2 > 0), 0 (λ 1 ≤ 0, λ 2 ≤ 0), (10)
which is strongly semismooth at any z ∈ R n [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF].

Lemma 2.1 [START_REF] Kanzow | On the local convergence of semismooth Newton methods for linear and nonlinear second-order cone programs without strict complementarity[END_REF] The projection mapping

P K n is differentiable at a point z = (z 1 , z 2 ) ∈ R × R n-1 if and only if z 1 = z 2 holds. In fact, the projection map- ping is continuously differentiable at every z such that z 1 = z 2 , i.e., z / ∈ bd(K n ), where bd(K n ) is the boundary of K n defined by bd(K n ) := {z = (z 1 , z 2 ) ∈ R × R n-1 : z 1 = z 2 }. Lemma 2.2 For any z = (z 1 , z 2 ) ∈ R × R n-1 , an element V of the B-subdifferential ∂ B P K n (z) has the following representation: V = I (λ 1 > 0, λ 2 > 0), i.e., z ∈ int(K n ), V = λ 2 λ 2 -λ 1 I + Z (λ 1 < 0, λ 2 > 0), i.e., z ∈ R n \ (K n ∪ -K n ), V = 0 n (λ 1 < 0, λ 2 < 0), i.e., z ∈ -int(K n ), V ∈ {I, I + Z } (λ 1 = 0, λ 2 > 0), i.e., z ∈ bd(K n ) \ {0}, V ∈ {0 n , Z } (λ 1 < 0, λ 2 = 0), i.e., z ∈ -bd(K n ) \ {0}, V = 0 n or V = I or V ∈ S (λ 1 = 0, λ 2 = 0), i.e., z ∈ {0}, (11) where
(r 1 , r 2 ) := (z 1 , z 2 ) z 2 , Z = 1 2 -r 1 r T 2 r 2 -r 1 r 2 r T 2 and S := 1 2 (1 + β)I + 1 2 -β ω T ω -βωω T : -1 ≤ β ≤ 1, ω = 1 .
Proof We can find a similar representation for the elements of

∂ B P K n (z)

Second-Order Eigenvalue Complementarity Problem

Throughout this section, we deal with the second-order cone eigenvalue complementarity problem defined in (3), which consists in finding a scalar λ > 0 and a vector x ∈ R n \{0} satisfying

x ∈ K, λx -Ax ∈ K, x, λx -Ax = 0.
The set

σ (A, K) := {λ > 0 : (λ, x) solves (3) for some x ∈ K \ {0}} (12)
is said to be the Lorentz spectrum of A. The component λ is called a Lorentz eigenvalue of A, while the component x is called a Lorentz eigenvector of A. The eigenvalue complementarity problem in the the case of the second-order cone has recently been studied. For this reason, the number of publications is rather limited (see, for instance, [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF]).

In [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF], Seeger and Torki studied the main properties and the structure of the Lorentz spectral mapping σ (., K). They recovered σ (A, K), defined in [START_REF] Lobo | Applications of second-order cone programming[END_REF], by joining two pieces and treated each one:

σ (A, K) = σ int (A, K) ∪ σ bd (A, K),
where σ int (A, K) (respectively, σ bd (A, K)) is the set of all eigenvalues λ of A associated with eigenvectors in the interior of K (respectively, in the boundary of K).

Moreover, they characterized the elements in the Lorentz spectrum and provided a valuable information on their number. They announced that σ (A, K) is not necessarily finite, but it can always be written as union of finitely many (at most 5n -4) mutually disjoint connected sets [35, Furthermore, they showed the classes of matrices over which the Lorentz cone can produce only a finite number of eigenvalues. Proposition 3.1 [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF] In each of the following cases, the real matrix A ∈ R n×n has a finite number of Lorentz eigenvalues:

(a) A is symmetric, (b) the leading (n -1) × (n -1) principal submatrix of A has no (real) eigenvalues with geometric multiplicity ≥ 2.
In what follows, we give some examples of discrete and continuous solutions.

Example 3.1 Consider the following matrix

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣
3 0 0 0 0 0 5 0 0 0 0 0 5 0 0 0 0 0 7 0 0 0 0 0 7

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ .
A has a discrete Lorentz spectrum given by

σ (A, K) = {3, 4, 5}. Example 3.2 Let B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣
3 2 3 0 2 3 0 0 5 0 0 0 0 0 5 0 0 0 0 0 7 0 0 0 0 0 7

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ .
We have

σ (B, K) = {3} ∪ [4.0009, 4.3333] ∪ [5.0016, 5.3333],
which is a continuous Lorentz spectrum of B.

As a first step toward a reformulation of the SOCEiCP as a system of equations, we write

x ∈ K, y ∈ K, x, y = 0, λx -Ax -y = 0, 1 n , x -1 = 0, (13)
where 1 n is a vector of ones. The normalization equation (the last equation in [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF]) is used to ensure that x is a nonzero vector. Without any loss of generality, we can suppose that λ is strictly positive. If the matrix A has a negative eigenvalue λ < 0, we set à = A + μI n , where μ > 0 is large enough. It is easy then to observe that the eigenvalues of à are strictly positive.

The SOCEiCP can be reformulated as the following nonsmooth system of equations

Φ(z) = Φ(x, y, λ) := ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ϕ 1 (x 1 , y 1 ) . . . ϕ r (x r , y r ) λx -Ax -y 1 n , x -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = 0, ( 14
)
where Φ : R n × R n × R → R 2n+1 and ϕ i : R n i × R n i → R n i is a second-order cone complementarity function (SOCC-function) associated with the cone K n i (n i ≥ 1) on Euclidean Jordan Algebras, i.e.,

ϕ i (x i , y i ) = 0 ⇐⇒ x i ∈ K n i , y i ∈ K n i , x T i y i = 0. ( 15
)
In this paper, we focus on the natural residual function, also called, the minfunction, denoted by ϕ soc min , and the Fischer-Burmeister function, denoted by ϕ soc FB , defined, respectively, by

ϕ soc min (x, y) := x -P K n (x -y), (16) ϕ soc FB (x, y) := x + y -(x 2 + y 2 ) 1/2 . ( 17
)
Note that the ϕ soc FB function defined above was proposed by Gowda et al. in [START_REF] Gowda | P-properties for linear transformations on Euclidean Jordan algebras[END_REF]. Xiu and al. established also the definition of ϕ soc FB and proved that it is Lipschitz continuous in [START_REF] Kong | Fischer-Burmeister complementarity function on Euclidean Jordan algebras[END_REF]. Furthermore, Sun and Sun in [START_REF] Sun | Strong semismoothness of the Fischer-Burmeister SDC and SOC complementarity functions[END_REF] showed that it is strongly semismooth everywhere, and so is the ϕ soc min function (see [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF] for more details). In what follows, we give the B-subdifferential of Φ in two cases. The first case is defined when the SOCC-function is given by the ϕ soc min function defined in [START_REF] Bonnans | Perturbation analysis of second-order cone programming problems[END_REF], and the second one is related to the ϕ soc FB function, defined in [START_REF] Chen | Alternative proofs for some results of vector-valued functions associated with second-order cones[END_REF].

Lemma 3.1 The function Φ min : R n × R n × R → R 2n+1
, defined in [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF] with ϕ being the ϕ soc min function, is semismooth. Moreover, each element H min of the B-subdifferential of Φ min at z = (x, y, λ) is given by

H min = ⎡ ⎣ I -V V 0 λI -A -I x 1 T n 0 0 ⎤ ⎦ , ( 18
)
where V ∈ ∂ B P K n (x -y).
Taking into account the proposition 3.1 in [START_REF] Pan | A damped Gauss-Newton method for the second-order cone complementarity problem[END_REF], we have the following. Given a vector (x, y) ∈ R n × R n , then each element V FB in ∂ B ϕ soc FB (x, y) is given by

V FB = [I -V x I -V y ], (19)
with V x and V y having the following representation:

(i) If x 2 + y 2 ∈ int(K n ), then V x = L -1 w L x and V y = L -1
w L y , where w := (x 2 + y 2 ) 1/2 , and for

x ∈ K n , L x = x 1 x T 2 x 2 x 1 I . If x ∈ int(K n ), where int(K n ) := {x = (x 1 , x 2 ) ∈ R × R n-1 : x 1 > x 2 }, then, L -1 x = 1 det(x) ⎡ ⎣ x 1 -x T 2 -x 2 det(x) x 1 I + x 2 x T 2 x 1 ⎤ ⎦ ,
where det(x) := x 2 1x 2 2 denotes the determinant of x.

(ii) If x 2 + y 2 ∈ bd(K n ) and (x, y) = (0, 0), then

V x ∈ 1 2 √ 2l 1 1 l 2 T l 2 4I -3 l 2 l 2 T L x + 1 2 1 -l 2 u T V y ∈ 1 2 √ 2l 1 1 l 2 T l2 4I -3 l 2 l 2 T L y + 1 2 1 -l 2 v T , where u = (u 1 , u 2 ), v = (v 1 , v 2 ) ∈ R × R n-1 satisfying |u 1 | ≤ u 2 ≤ 1 and |v 1 | ≤ v 2 ≤ 1, where l = (l 1 , l 2 ) := x 2 + y 2 and l 2 = l 2 / l 2 . (iii) If (x, y) = (0, 0), then V x ∈ {L x }, V y ∈ {L ỹ } for some x, ỹ with x 2 + ỹ 2 = 1, or V x ∈ 1 2 1 l2 ξ T + 1 2 1 -l 2 u T + 2 0 0 (I -l 2 l 2 T )s 2 (I -l 2 l 2 T )s 1 V y ∈ 1 2 1 l2 η T + 1 2 1 -l 2 v T + 2 0 0 (I -l 2 l 2 T )l 2 (I -l 2 l 2 T )l 1 , where u = (u 1 , u 2 ), v = (v 1 , v 2 ), ξ = (ξ 1 , ξ 2 ), η = (η 1 , η 2 ) ∈ R × R n-1 satisfying |u 1 | ≤ u 2 ≤ 1, |v 1 | ≤ v 2 ≤ 1, |ξ 1 | ≤ ξ 2 ≤ 1 and |η 1 | ≤ η 2 ≤ 1, where l 2 ∈ R n-1 such that l 2 = 1, and s = (s 1 , s 2 ), l = (l 1 , l 2 ) ∈ R × R n-1 satisfying s 2 + l 2 ≤ 1/2.
For more details, see [START_REF] Pan | A damped Gauss-Newton method for the second-order cone complementarity problem[END_REF]. Therefore, we get the following lemma.

Lemma 3.2 The function Φ

FB : R n × R n × R → R 2n+1 ,
H FB = ⎡ ⎣ I -V x I -V y 0 λI -A -I x 1 T n 0 0 ⎤ ⎦ , ( 20
)
where V x and V y are given by the above equations.

The Lattice Projection Method (LPM)

In this section, we generalize the LPM given in [START_REF] Adly | A new method for solving eigenvalue complementarity problems[END_REF] to solve SOCEiCP, which leads us to solve a nonlinear and nonsmooth system of (2n + 1) equations involving the same number of variables. The originality of this formulation is that it is not based on the complementarity approach.

Lemma 4.1 The SOCEiCP, with λ > 0, consists in finding the roots of the following nonlinear and nonsmooth function h : R n × R → R n , (x, λ) → h(x, λ), defined by h(x, λ)

:= P K (Ax) -λx. ( 21
)
Proof Suppose that λ > 0 and x ∈ K. By the definition of the normal cone, we have

K x ⊥ λx -Ax ∈ K ⇐⇒ Ax -λx ∈ N K (x).
Therefore SOCEiCP is equivalent to

x = P K ( 1 λ Ax).
Hence,

P K (Ax) = λx.
In short, we are led to solve the following system of (2n + 1) equations

P K ( ỹ) -λx = 0, Ax -ỹ = 0, 1 n , x -1 = 0. Remark 4.1 Lemma 4.1 showed that K x ⊥ (λx -Ax) ∈ K ⇐⇒ P K (Ax) = λx,
which means that the SOCEiCP is equivalent to solve the nonlinear eigenvalue problem

(P K • A)(x) = λx.
The SOCEiCP can be written equivalently to the following nonsmooth system of equations:

Φ LPM (z) := ⎡ ⎣ P K ( ỹ) -λx Ax - ỹ < 1 n , x > -1 ⎤ ⎦ = 0, ( 22
)
where

Φ LPM : R n × R n × R → R 2n+1 , z = (x, ỹ, λ) → Φ LPM (z).
Remark 4.2 We note that the solutions of the nonlinear equation ( 22) do not satisfy necessarily λ > 0. Maintaining the non-negativity of λ in the algorithm is not an easy task. If the matrix A has a nonpositive eigenvalue λ < 0, then we set à = A+μI n with μ > 0 and large enough. It is easy to see in this case that the spectrum σ ( Ã) ⊂ R * + .

The Algorithm 1 presented in the following part is used to solve [START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF]. This method will be called LPM.

The following lemma gives a brief description of the representation of the Bsubdifferential of Φ LPM . [START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF], is semismooth. Moreover, each element H of the B-subdifferential of Φ LPM at z = (x, ỹ, λ) is given by

Lemma 4.2 The function Φ

LPM : R n × R n × R → R 2n+1 , defined in
H = ⎡ ⎣ -λI Ṽ -x A -I 0 1 T n 0 0 ⎤ ⎦ , ( 23
)
where Ṽ ∈ ∂ B P K ( ỹ).

Proof Chen et al. [START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF] showed that P K (•) is strongly semismooth, so that Φ LPM is semismooth. The computation of the B-subdifferential of Φ LPM offers no difficulty and we get obviously the representation formula [START_REF] Yamashita | A primal-dual interior point method for nonlinear optimization over second order cones[END_REF].

The algorithm we are about to describe is the Semismooth Newton Method SNM for solving SOCEiCP.

Algorithm 1 (Semismooth Newton Method SNM)

Initialization. Choose an initial point z 0 and set k = 0. Iteration. 1 One has a current point z k . Choose M k ∈ ∂ B Φ(z k ) and compute h k by solving the linear system

M k h k = -Φ(z k ). ( 24
)
Then, set

z k+1 = z k + h k and k = k + 1,
where Φ is either Φ LPM or Φ min or Φ FB . To ensure that the linear system (24) admits a unique solution, the matrix M k must be nonsingular. Furthermore, the following theorem gives some conditions to ensure the local convergence of the above algorithm [START_REF] Qi | A nonsmooth version of Newton's method[END_REF].

Theorem 4.1 Let z be a zero of the function Φ. Suppose the following (i) Φ is semismooth (resp. strongly semismooth) at z ;

(ii) all matrices M k ∈ ∂ B Φ(z) are nonsingular. Then, there exists a neighborhood V of z such that the SNM, initialized at any z 0 ∈ V , generates a sequence (z k ) k∈N that converges superlinearly (resp. quadratically) to z.

Nonsingularity Conditions

We study, in this section, under which conditions the Jacobain matrix, defined in (24) of the Algorithm 1, is nonsingular at a solution. More precisely, we study the nonsingularity conditions of the elements of the B-subdifferential of Φ(z * ) and Φ LPM (z * ) defined, respectively, in ( 14) and [START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF]. Firstly, we give some notations, which will be used later.The index set {1, 2, ..., n} will be abbreviated by the capital letter I .

If M = (m i j ) ∈ M n (R) is a matrix and J , K ⊆ I , then M J K denotes the submatrix in M |J |,|K| (R) with elements m i j , i ∈ J , j ∈ K.
We introduce the following function

F of class C 1 F : R n × R n × R → R 2n+1 , z = (x, y, λ) → F(z)
given by

F(z) = ⎡ ⎣ y λx -Ax -y 1 n , x -1 ⎤ ⎦ . ( 25
)
Proposition 5.1 Any H ∈ ∂ B Φ(z), with ϕ being one of the SOCC-functions ϕ soc min or ϕ soc FB , can be written in the form

H = D a + D b F (z), ( 26
)
where D a , D b ∈ M 2n+1 (R) are positive semidefinite diagonal matrices, defined by

D a = ⎡ ⎣ I -V 0 0 0 0 0 0 0 0 ⎤ ⎦ and D b = ⎡ ⎣ V 0 0 0 I 0 0 0 1 ⎤ ⎦ , ( 27
)
in the case where ϕ(x, y) = ϕ soc min (x, y).

D a = ⎡ ⎣ I -L -1 w L x 0 0 0 0 0 0 0 0 ⎤ ⎦ and D b = ⎡ ⎣ I -L -1 w L y 0 0 0 I 0 0 0 1 ⎤ ⎦ , ( 28
)
in the case where ϕ(x, y) = ϕ soc FB (x, y) and x 2 + y 2 ∈ int(K n ). Proof If H ∈ ∂ B Φ min (z), then by Lemma 3.1, we have

H = H min = ⎡ ⎣ I -V V 0 λI -A -I x 1 T n 0 0 ⎤ ⎦ .
If H ∈ ∂ B Φ FB (z), then using Lemma 3.2, we get

H = H FB = ⎡ ⎣ I -L -1 w L x I -L -1 w L y 0 λI -A -I x 1 T n 0 0 ⎤ ⎦ .
On the other hand, we have

F (z) = ⎡ ⎣ 0 I 0 λI -A -I x 1 T n 0 0 ⎤ ⎦ . ( 29
)
Therefore, from these equalities and using [START_REF] Kanzow | A new class of semismooth Newton type methods for nonlinear complementarity problems[END_REF], we get the desired result. Now let z * = (x * , y * , λ * ) be a solution of ( 14), and from the vector function F defined in [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF], set

G : R n × R → R n+1 , w * = (x * , λ * ) → G(w * ) defined by G(w * ) = G(x * , λ * ) = λ * x * -Ax * 1 n , x * -1 , (30)
and

J = a + b G (w * ), (31)
where

a = I -V 0 0 0 , b = V 0 0 1 in the case ofϕ soc min . ( 32
)
a = I -L -1 w L x 0 0 0 , b = I -L -1 w L (λ * x * -Ax * ) 0 0 1 in the case of ϕ soc FB , (33)
and G (w * ) is given by

G (w * ) = λ * I -A x * 1 T n 0 . ( 34
)
Given a fixed solution z * = (x * , y * , λ * ) of ( 14) and let w * = (x * , λ * ).

The point w * is said to be nondegenerate iff

w * i + G i (w * ) ∈ int(K n i
), for all i ∈ {1, . . . , r + 1}. Then, the index sets

α := {i ∈ {1, . . . , r + 1} : w * i ∈ int(K n i ), G i (w * ) = 0}, β := {i ∈ {1, . . . , r + 1} : w * i ∈ bd + (K n i ), G i (w * ) ∈ bd + (K n i )}, γ := {i ∈ {1, . . . , r + 1} : w * i = 0, G i (w * ) ∈ int(K n i )},
where bd + (K n i ) = bd(K n i ) \ {0}, form a partition of {1, . . . , r + 1}. Moreover, we set

˜ a := I -Ṽ 0 0 0 , ˜ b := Ṽ 0 0 1 , where Ṽ ∈ ∂ B P K ( ỹ).
We note that λ * ˜ a = ˜ λ * a , and we set

J = ˜ λ * a + ˜ b G (w * ). ( 35
)
Proposition 5.2 Let z * = (x * , y * , λ * ) be a fixed solution of( 14) with w * = (x * , λ * ) being nondegenerate. If the submatrix G (w * ) αα is nonsingular and if its Schurcomplement

G (w * ) ββ -G (w * ) βα G (w * ) -1 αα G (w * ) αβ
is a P-matrix, then, the Jacobians J and J , defined, respectively, in [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF] and [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF], are nonsingular matrices.

Proof Solving SOCEiCP is equivalent to solving the following system

find w ∈ R n+1 \{0}, such that w ∈ K, G(w) ∈ K, w, G(w) = 0, ( 36
)
where G is given by [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. The desired result follows immediately from Theorem 2.8 in [START_REF] Kanzow | A new class of semismooth Newton type methods for nonlinear complementarity problems[END_REF].

Remark
(λ * I -A) ββ -(λ * I -A) βα (λ * I -A) -1 αα (λ * I -A) αβ is a P-matrix.
Therefore, we get the following generalized theorems to the second-order cone. Proof Let H = H min ∈ ∂ B Φ min (z * ) and let X = ( p, q, r ) ∈ R n × R n × R be an arbitrary vector such that

H X = 0 R 2n+1 . ( 37
)
The previous equation ( 37) can be written explicitly as

⎡ ⎣ I -V V 0 λ * I -A -I x * 1 T n 0 0 ⎤ ⎦ ⎡ ⎣ p q r ⎤ ⎦ = ⎡ ⎣ 0 R n 0 R n 0 R ⎤ ⎦ ,
where V ∈ ∂ B P K n (xy). Therefore,

(I -V ) p + V q = 0 R n , ( 38
) (λ * I -A) p -q + r x * = 0 R n , ( 39
) 1 n , p = 0. ( 40
)
Using [START_REF] Kanzow | A new class of semismooth Newton type methods for nonlinear complementarity problems[END_REF], we have q = (λ * I -A) p + r x * . (41)

Consequently, [START_REF] Pan | A damped Gauss-Newton method for the second-order cone complementarity problem[END_REF] becomes

(I -V ) p + V (λ * I -A) p + r V x * = 0 R n . ( 42
)
Therefore,

(I -V ) p + V (λ * I -A) p + r V x * = 0 R n , 1 n , p = 0,
which can be rewritten in the following matricial form

I -V + V (λ * I -A) V x * 1 T n 0 p r = 0 R n 0 R . ( 43
)
By setting Y = ( p, r ) ∈ R n × R, clearly, ( 43) is equivalent to

J Y = 0 R n+1 . ( 44
)
Hence,

H X = 0 R 2n+1 ⇐⇒ J Y = 0 R n+1 , q = (λ * I -A) p + r x * ,
which completes the proof. For the case of the elements H FB ∈ ∂ B Φ FB (z * ), it is sufficient to apply the same steps given above.

Theorem 5.2 Assume that z * = (x * , ỹ * , λ * ) is a fixed root of Φ LPM defined in [START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF] with w * = (x * , λ * ) nondegenerate. Then the Jacobian J defined in [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF] is nonsingular if and only if all elements in ∂ B Φ LPM (z * ) are nonsingular.

Proof Let H ∈ ∂ B Φ LPM (z * ) and let X = ( p, q, r ) ∈ R n × R n × R such that

H X = 0 R 2n+1 . ( 45
)
The previous equation can be written explicitly as

⎡ ⎣ -λ * I Ṽ -x * A -I 0 1 T n 0 0 ⎤ ⎦ ⎡ ⎣ p q r ⎤ ⎦ = ⎡ ⎣ 0 R n 0 R n 0 R ⎤ ⎦ ,
where Ṽ ∈ ∂ B P K ( ỹ). Hence,

-λ * p + Ṽ q -r x * = 0 R n , (46) Ap -q = 0 R n , (47) 1 n , p = 0. ( 48
)
Using (47), we have q = Ap.

(49)

Consequently, (46) becomes

(λ * I -Ṽ A) p + r x * = 0 R n . (50) Therefore, (λ * I -Ṽ A) p + r x * = 0 R n , 1 n , p = 0. ( 51
)
Note that

I = I -Ṽ + Ṽ .
Then, (51) becomes

λ * (I -Ṽ ) p + Ṽ (λ * I -A) p + r x * = 0 R n , 1 n , p = 0. ( 52
)
We note also that we have

Ṽ x * = x * .
Consequently, (52) can be rewritten in the following matricial form

λ * (I -Ṽ ) + Ṽ (λ * I -A) Ṽ x * 1 T n 0 p r = 0 R n 0 R . ( 53
) By setting Y = ( p, r ) ∈ R n × R, it is clear that (53) is equivalent to J Y = 0 R n+1 . (54)
Hence,

H X = 0 R 2n+1 ⇐⇒ J Y = 0 R n+1 , q = Ap.
The desired conclusion of Theorem 5.2 follows immediately from the last equivalence.

Numerical Experiments

In this section, we apply the semismooth Newton method SNM, defined in Algorithm 1, in order to solve the SOCEiCP and give some numerical results to verify the efficiency of the LPM. The latter will be compared with two other solvers. The first one is the semismooth Newton method, related to ϕ soc min defined in (16) (denoted by SNM min ), and the second one is related to ϕ soc FB defined in (17) (denoted by SNM FB ).

We choose the performance profiles developed in Dolan and Moré [24,[START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] as a tool to compare the solvers. Computing time, number of functions evaluations and number of failures are used as performance measures to compare the three solvers.

The numerical experiments are carried out in a Powerbook Mac OS 10.6.8 with a processor 2.33 GHz Intel Core 2 Duo and 2Go memory. All the program codes are written and executed in Matlab 7.7.

Let S be the set of the three solvers that will be compared. Before establishing the figures, we give some notations :

• r p,s denotes the performance ratio defined by r p,s := t p,s min{t p,s : s ∈ S} ,

• p is an element of the set of matrices P = 160 taken to compare the solvers, with 10,000 initial points,

• s denotes each solver,

• n s is the number of solvers,

• t is a real factor, and

• t p,s denotes the average number of functions evaluations, the average computing time and the number of failures, respectively, in Figs. 1, 2, and 3.

Remark 6.1 In our case, a failure is declared if the Jacobian matrix is ill-conditioned or in the SNM defined in Algorithm 1, the number of iterations exceeded 100. 3 Performance profiles where t p,s is the average number of failures Figure 1 shows the performance profiles given by the three solvers, i.e., LPM, SNM FB and SNM min . We compare the number of functions evaluations required to find a solution by each solver. Clearly, LPM has the most wins (since it has the highest probability). In the interval [0, 0.5], LPM can solve 90 % of the problems, while SNM min and SNM FB do not attain the 30 %. When t ≥ 2, the performance of the latters becomes interesting.

In Fig. 2, the computing time is the comparing tool. Clearly, the performance profile of the LPM is interesting since it has the highest probability to solve problems. On the other hand, SNM FB and SNM min require a more running time to find a solution and have a lower number of wins.

Figure 3 shows the efficiency of the LPM, where the average number of failures to find each solution is the comparing tool. Clearly, in the interval [0, 0.5], LPM is able to solve about 90 % of problems and detects the Lorentz eigenvalues while the others do not reach the 50 %. We see also that when t ≥ 2, the performance of SNM FB becomes interesting. This allows us to conclude that the number of failures (in the sense of Remark 6.1) given by SNM min is high, which prevents to find solutions.

Globally Convergent Methods

In this section, we will discuss briefly a globalization of the three methods studied in the last sections by constructing appropriate merit functions associated to Φ min , Φ FB and Φ LPM defined, respectively, in ( 14) and [START_REF] Tseng | Smoothing methods for second-order cone programs/complementarity problems[END_REF].

For the function Φ FB , we introduce the following classical merit function Ψ FB : R 2n+1 → R 2n+1 , defined by

Ψ FB (z) := 1 2 Φ FB (z) T Φ FB (z). ( 55
)
We note that the function Ψ FB is of class C 1 and that SOCEiCP is also equivalent to the following unconstrained global optimization problem min z∈R 2n+1

Ψ FB (z).

Since Ψ FB is continuously differentiable, it is easy to force global convergence of the algorithm by using the gradient of the merit function (see [START_REF] Facchinei | A semismooth Newton type method for variational inequalities: the case of box constraints[END_REF]).

Algorithm 2 (Semismooth Newton Method SNM)

Step 0. (Data) Choose parameters δ, ρ ∈ ]0, 1[, σ ∈ ]0, 0.5[, κ > 2 and initial point z 0 := (x 0 , y 0 , λ 0 ) ∈ R n × R n × R. Give the error bound 0 < << 1 and set k = 0.

Step 1. (Termination criteria) If ∇Ψ F B (z k ) ≤ , stop, where the function Ψ FB (•) is defined in (55).

Step 2. (Search direction calculation) At a current point z k , we select M k ∈ ∂Φ FB (z k ) and compute h k by solving the linear system M k h k = -Φ FB (z k ).

(56)

If system (56) is unsolvable or if h k does not satisfy the following condition

∇Ψ F B (z k ) T h k ≤ -ρ h k κ , ( 57
)
then set

h k = -∇Ψ F B (z k ).
the study of the equilibrium state of mechanical systems with unilateral contact [START_REF] Martins | The directional instability problem in systems with frictional contacts[END_REF].

Another interesting question would be the generalization of the LPM method for solving quadratic eigenvalue complementarity problems. This is out of the scope of this paper and will probably be the subject of a future project of research.

  Theorem 4.1 and Corollary 4.4] (see examples 3.1 and 3.2).

Fig. 1 Fig. 2

 12 Fig.1Performance profiles where t p,s is the average number of functions evaluations

  5.1 If r < n, i.e., α ⊂ {1, . . . , n + 1}, then G (w * ) αα is nonsingular if and only if (λ * I -A) αα is nonsingular since in this case, we have G (w * ) αα = (λ * I -A) αα . Moreover, its Schur-complement is a P-matrix if and only if

  ) nondegenerate. Then, all elements H in the B-subdifferential ∂ B Φ(z * ), with ϕ being one of the SOCC-functions ϕ soc min or ϕ soc FB , are nonsingular if and only if J is nonsingular.

Theorem 5.1 Assume that z * = (x * , y * , λ * ) is a fixed solution of

[START_REF] Chen | Complementarity functions and numerical experiments for second-order cone complementarity problems[END_REF]

with w * = (x * , λ *

As stopping criteria, we use Φ(z k ) < 10 -8 .

Step 3. (Line search) Let m k be the smallest non-negative integer m such that

and let α k := δ m k and z k+1 := z k + α k h k .

Step 4. (Update) Set k := k + 1 and go to Step 1.

We have the following convergence result (see [START_REF] Facchinei | A semismooth Newton type method for variational inequalities: the case of box constraints[END_REF] for more details).

Theorem 7.1 [START_REF] Facchinei | A semismooth Newton type method for variational inequalities: the case of box constraints[END_REF] Let F be a C We can construct also merit functions Ψ min and Ψ LPM associated to the natural residual function Φ min and Φ LPM defined, respectively, in ( 14) and ( 22):

The drawback of these two functions are their non-differentiability. As an immediate consequence, numerical methods based on the gradient of the function, such as steepest descent method and Newton's method, cannot be applicable directly to Ψ min and Ψ LPM .

To overcome this difficulty, Fukushima et al. [START_REF] Fukushima | Smoothing functions for second-order cone complementarity problems[END_REF] and Hayashi et al. [START_REF] Hayashi | A combined smoothing and regularization method for monotone second-order cone complementarity problems[END_REF] proposed smoothing and regularization methods by replacing the original non differentiable function by a sequence of differentiable approximations. The boundedness of the levelsets (or the coercivity) of the objective function plays an important role to ensure the global convergence of a descent method. Further investigations must be done in this direction for SOCEiCP.

Conclusions

In this paper, we have studied numerical methods for solving the second-order cone eigenvalue complementarity problem. We have reformulated such a problem into a system of semismooth equations. Moreover, we have extended the LPM, introduced in [START_REF] Adly | A new method for solving eigenvalue complementarity problems[END_REF], to solve the SOCEiCP. The performance profiles highlight the efficiency of the LPM method for solving SOCEiCP. It will be interesting to compare the globalization of the three methods in the future (using remark 7.1). Another natural question is to generalize the LPM to the mixed eigenvalue complementarity problems (where only some components of the vector x ∈ R n are cone-constrained and the identity I in (3) is replaced by some matrix of order n). This kind of problems appears in


Samir Adly
email: adly@unilim.fr

Emil Ernst
email: ernst@imar.ro

Michel Théra
email: michel.thera@unilim.fr

A Characterization of Convex and Semicoercive Functionals

Keywords: convex analysis, barrier cone, support functional, recession analysis, semicoercive functional

In this paper we prove that every proper convex and lower semicontinuous functional Φ defined on a real reflexive Banach space X is semicoercive if and only if every small uniform perturbation of Φ attains its minimum value on X.

Introduction

Throughout the paper X denotes a real reflexive Banach space, X * its topological dual and •, • the associated duality pairing on X × X * . We write "→ and "⇀ to denote respectively the strong and the weak convergence on X. We denote by B(x, r) (respectively by B(x, r)) the open ball (respectively closed ball) with center x and radius r > 0. Following the standard notations used in convex analysis, Γ 0 (X) stands for the set of all convex lower semicontinuous proper (not identically equal to +∞) extended-real-valued functionals Φ : X → R ∪ { +∞}.

By argmin Φ we mean the (possibly empty) set of all x ∈ X where Φ attains its minimum value, i.e. Φ(x) ≤ Φ(y), ∀ y ∈ X.

Recall that a functional Φ : X → R ∪ { +∞} is said to be semicoercive if there exists a closed subspace V of X such that Φ(x) = Φ(x + v), ∀x ∈ X, ∀ v ∈ V, and the quotient functional Φ : X/V → R ∪ { +∞}, defined by: Φ(x) = Φ(x), ∀x ∈ x, ∀ x ∈ X/V, is coercive, in the sense that Φ -1 ((-∞, M ]) is bounded for every M ∈ R.

Let us remark that every coercive functional is semicoercive. Let us also mention two of the most frequently encountered semicoercive functionals:

(i) The distance functional to a closed subspace F of an arbitrary Banach space X,

J 1 : X → R, J 1 (x) = dist(x, F ) = inf y∈F x -y .
(ii) If Ω ⊆ R n is a bounded subset with a smooth boundary, and H 1 (Ω) is the corresponding Sobolev space,

J 2 : H 1 (Ω) → R, J 2 (u) = Ω |∇u(x)| 2 dx.
By a uniform perturbation of Φ we mean a functional Ψ which satisfies for some ε > 0

Φ(x) -ε ≤ Ψ(x) ≤ Φ(x) + ε , ∀x ∈ X.
Obviously, every convex uniform perturbation of a semicoercive functional remains semicoercive. It is a well-known fact that every Γ 0 (X)-functional, which is semicoercive attains its minimum on X. Consequently, in the class of Γ 0 (X)-semicoercive functionals, the existence of a minimum on X is preserved under every uniform perturbation.

Our aim in this paper is to study if the converse of the preceding observation is true. More precisely, we prove (Theorem 4.2) that if Φ ∈ Γ 0 (X) and if every small uniform perturbation of Φ within the class Γ 0 (X) reaches its minimum value on X, then Φ is necessarily semicoercive.

The paper is organized as follows. Section 2 contains some classical convex analysis results which are needed throughout the paper, together with the proof of the technical Lemma 2.1. Section 3 is devoted to the proof of Theorem 3.1, which characterizes unbounded, linearly bounded convex closed sets (that is unbounded sets with a recession cone reduced to the singleton {0}). This theorem plays a central role in the proof of the main result which is given in Section 4.

The final section contains the conclusions and some open problems.

Background material and preliminary results

Let us recall some basic definitions from convex analysis. For Φ ∈ Γ 0 (X), we denote by Dom Φ, the effective domain of Φ which is defined by Dom Φ = {x ∈ X : Φ(x) < +∞} and by epi Φ the epigraph of Φ. Let K be a nonempty convex and closed subset of X.

Following Rockafellar [START_REF] Rockafellar | Convex Analysis[END_REF], the recession cone K ∞ of K is defined by:

K ∞ = t>0 t K -x 0 , x 0 ∈ K.
Note that K ∞ is independent of x 0 ∈ K and that K ∞ is a closed convex cone of X and describes the global behavior of the convex set K.

For each Φ ∈ Γ 0 (X), taking C = epi Φ we define the recession functional of Φ as the function Φ ∞ such that epi Φ ∞ = (epi Φ) ∞ . Equivalently, it amounts to saying that

Φ ∞ (x) := lim t→+∞ Φ(x 0 + tx) -Φ(x 0 ) t ,
where x 0 is any point in Dom Φ.

In the sequel, Ker Φ ∞ will denote the closed convex cone defined by

Ker Φ ∞ = {x ∈ X : Φ ∞ (x) = 0}.
The following properties of Φ ∞ will be used in the sequel. The proofs of ( 1)-( 3) are straightforward and were established in finite dimension in [START_REF] Rockafellar | Convex Analysis[END_REF]; the reader is for instance refered to [START_REF] Baiocchi | General existence results for unilateral problems in continuum mechanics[END_REF] for [START_REF] Auslender | Noncoercive optimization problems in reflexive Banach spaces[END_REF].

Φ ∞ ∈ Γ 0 (X) (1) Φ(x + y) ≤ Φ(x) + Φ ∞ (y), ∀x, y ∈ X (2) Φ ∞ (λx) = λΦ ∞ (x) ∀x ∈ X, ∀λ ≥ 0 (3) Φ ∞ (x) ≤ lim inf n→∞ Φ(x 0 + s n x n ) s n , (4)
where x 0 is any point in Dom Φ, (x n ) n∈N * is any sequence in X converging to x and (s n ) n∈N * is any real sequence converging to +∞.

Throughout the paper, x 0 ∈ X and M ∈ R are defined such that :

x 0 ∈ argmin Φ and M := Φ(x 0 ).

Denoting by co {(y, N ), epi Φ} the convex hull of (y, N ) and epi Φ, for y ∈ X and N < M , we set

F(x) = {s ∈ R : (x, s) ∈ co{(y, N ), epi Φ}},
the closure being taken with respect to the natural strong topology on X × R. Obviously, F(x) is a closed, convex, (possibly empty) subset of R; since

∅ = co{(y, N ), epi Φ} ⊆ X × [N, +∞) ⊂ X × R, it follows that F(x) ⊂ [N, ∞), ∀ x ∈ X.
Therefore, the relation

Φ y,N (x) = min F(x) if F(x) = ∅ +∞ if F(x) = ∅
defines a proper extended-real-valued functional, whose epigraph is co{(y, N ), epi Φ}. Some of the most important properties of Φ y,N are stated in the following Lemma.

Lemma 2.1. For every functional Φ y,N , the following holds:

(i) Φ y,N belongs to Γ 0 (X); (ii) Φ y,N is a minorant of Φ: Φ y,N (x) ≤ Φ(x), ∀ x ∈ X;
If, in addition, Φ attains its minimum value at y, then

Φ(x) + N -M ≤ Φ y,N (x), ∀ x ∈ X;
(iii) for every x ∈ X such that Φ y,N (x) ≤ M we have

Φ y,N (λy + (1 -λ)x) = λΦ y,N (y) + (1 -λ)Φ y,N (x), ∀ λ ∈ (0, 1), (5)
(iv) argmin Φ y,N = y + Ker Φ ∞ .

Proof of Lemma 2.1. (i) Since epi Φ y,N = co{(y, N ), epi Φ} is a non-empty, closed and convex set, Φ y,N ∈ Γ 0 (X).

(ii) As epi Φ ⊂ co{(y, N ), epi Φ},

we have Φ y,N (x) ≤ Φ(x) ∀ x ∈ X.
If, in addition, Φ(y) = M , we have (y,

N ) ∈ epi(Φ + N -M ) and epi Φ ⊂ epi(Φ + N -M ), so co{(y, N ), epi Φ} ⊂ epi(Φ + N -M ), that is Φ(x) + N -M ≤ Φ y,N (x), ∀ x ∈ X.
(iii) We shall distinguish the cases Φ y,N (x) < M , and Φ y,N (x) = M . Let us first suppose that Φ y,N (x) < M . Set z = λy + (1λ)x; for every δ > 0 set

α := min {δ(1 -λ), (M -Φ y,N (x))(1 -λ)} . (6)
Since (z, Φ y,N (z)) ∈ co{(y, N ), epi Φ}, there exist (w, θ) ∈ epi Φ and 0 ≤ µ ≤ 1 (depending on δ), such that

dist X×R µ(y, N ) + (1 -µ)(w, θ), (z, Φ y,N (z)) 2 = (µy + (1 -µ)w) -z 2 + |(µN + (1 -µ)θ) -Φ y,N (z)| 2 ≤ α 2 . Therefore µN + (1 -µ)θ ≤ Φ y,N (z) + α.
Since M ≤ θ and α ≤ (M -Φ y,N (x))(1λ), the previous relation yields

µN + (1 -µ)M ≤ Φ y,N (z) + (M -Φ y,N (x))(1 -λ). (7)
Using the convexity of Φ y,N , we deduce that

Φ y,N (z) + (M -Φ y,N (x))(1 -λ) (8) ≤ λΦ y,N (y) + (1 -λ)Φ y,N (x) + (1 -λ)(M -Φ y,N (x)) = λN + (1 -λ)M.
By [START_REF] Goeleven | On Noncoercive Variational Problems and Related Results[END_REF] and [START_REF] Laurent | Approximation et Optimisation[END_REF]

we get µN + (1 -µ)M ≤ λN + (1 -λ)M.
Consequently, λ ≤ µ, and therefore 0 ≤ λ 1 = (µλ)/(1λ) ≤ 1. This yields

λ 1 (y, N ) + (1 -λ 1 )(w, θ) ∈ co {(y, N ), epi Φ}.
Thus, since the relation

dist (x, 1 1 -λ Φ y,N (z) - λ 1 -λ Φ y,N (y)), λ 1 (y, N ) + (1 -λ 1 )(w, θ) = 1 1 -λ dist (z, Φ y,N (z)), µ(y, N ) + (1 -µ)(w, θ) ≤ α 1 1 -λ ≤ δ,
holds for every δ > 0, we derive that

x, 1 1 -λ Φ y,N (z) - λ 1 -λ Φ y,N (y) ∈ co{(y, N ), epi Φ}. Hence, Φ y,N (x) ≤ 1 1 -λ Φ y,N (λy + (1 -λ)x) - λ 1 -λ Φ y,N (y).
This relation combined with the convexity of Φ y,N implies (5).

Let us now consider the case Φ y,N (x) = M . Set

x n = 1 n y + n-1 n x; since Φ y,N is convex, Φ y,N (x n ) ≤ 1 n Φ y,N (y) + n -1 n Φ y,N (x) = 1 n N + n -1 n M < M.
Hence relation (5) holds for every x n . Therefore, for every λ in (0, 1) and every n in N * we have Φ y,N (λy

+ (1 -λ)x n ) = λΦ y,N (y) + (1 -λ)Φ y,N (x n ). (9)
On [0, 1], the mappings µ → Φ y,N (µy +(1-µ)x) and µ → Φ y,N (λy +(1-λ)(µy +(1-µ)x)) are convex and lower semicontinuous , thus continuous. Consequently, we obtain relation [START_REF] Baiocchi | General existence results for unilateral problems in continuum mechanics[END_REF] by taking in [START_REF] Rockafellar | Convex Analysis[END_REF] the limit as n goes to infinity.

(iv) Take x in X such that Φ y,N (x) = N . We may assume without loss of generality that x = y. Since (x, N ) belongs to co{(y, N ), epi Φ}, for every n in N * there are (x n , θ n ) in epi Φ and

µ n in [0, 1) such that dist µ n (y, N ) + (1 -µ n )(x n , θ n ), (x, N ) ≤ 1 n . (10)
By [START_REF] Valentine | Convex Sets[END_REF] we deduce, in particular, that

µ n N + (1 -µ n )θ n ≤ N + 1 n , and since θ n ≥ M , we obtain 1 -µ n ≤ 1 n(M -N ) . (11)
Relation [START_REF] Valentine | Convex Sets[END_REF] implies that lim

n→∞ µ n N + (1 -µ n )θ n -N = 0; as from (11) it follows that lim n→∞ µ n = 1, we deduce that lim n→∞ (1 -µ n )θ n = 0. Since Φ(x n ) ≤ θ n , the previous relation yields lim sup n→∞ (1 -µ n )Φ(x n ) ≤ 0. ( 12
)
On the other hand, relation [START_REF] Valentine | Convex Sets[END_REF] implies that

lim n→∞ µ n y + (1 -µ n )x n -x = 0; since lim n→∞ µ n = 1, we obtain that lim n→∞ (1 -µ n )x n -(x -y) = 0, which implies lim n→∞ (1 -µ n )(x n -z 0 ) -(x -y) = 0,
where z 0 is an arbitrary element of Dom Φ. We may therefore apply (4) for x -

y instead of x, 1/(1 -µ n ), and (1 -µ n )(x n -z 0 ) instead of x n and s n := 1/(1 -µ n ) (note that µ n = 1)
, and as

x n = z 0 + 1 1 -µ n [(1 -µ n )(x n -z 0 )],
we obtain

0 ≤ Φ ∞ (x -y) ≤ lim inf n→∞ (1 -µ n )Φ(x n ). (13)
By ( 12) and (13) it follows that xy belongs to Ker Φ ∞ , that is x ∈ y + Ker Φ ∞ .

Since obviously Φ y,N reaches its minimum value at every point of y + Ker Φ ∞ , the proof of Lemma 2.1 is established.

Unbounded linearly bounded closed convex sets

Recall that a closed convex set K is linearly bounded if K ∞ = {0}.

Theorem 3.1 below states an important property of convex, closed, unbounded and linearly bounded sets. This result shows that, contrary to the finite dimensional setting (see [START_REF] Rockafellar | Convex Analysis[END_REF], chapter 8), in a general linear space the recession cone does not characterize completely the behavior "at infinity of a convex set.

Theorem 3.1. Let K be a closed, convex, unbounded and linearly bounded subset of a reflexive Banach space X. Then, there exists h ∈ X * such that

inf w∈K h, w < h, u ≤ 1, ∀u ∈ K. ( 14
)
For every closed convex subset K of X we define the barrier cone B(K) of K as the domain of the support functional σ K of K defined by σ K (f ) := sup x∈K f, x . In other words,

B(K) = {f ∈ X * : σ K (f ) < +∞} = Dom σ K .
The following lemma is an elementary consequence of a known result.

Lemma 3.2. Suppose K is a closed convex and linearly bounded set of X. Then B(K) is dense in X * .

Proof of Lemma 3.2. It is a well-known fact that the recession cone K ∞ of K is the polar of B(K) (see for instance [START_REF] Rockafellar | Convex Analysis[END_REF] 14.2.1 in finite dimension and [START_REF] Aubin | Optima and Equilibria. An introduction to Nonlinear Analysis[END_REF], Proposition 3.10 in infinite dimension). Therefore, by the bipolar theorem we obtain X * = B(K) oo = B(K) (since X is reflexive) and the proof is complete.

Given R > 0 we define

B R (K) := {f ∈ B(K) : ∃x ∈ Ksuch that f, x > R}.
Lemma 3.3. Suppose K is a convex closed, unbounded and linearly bounded subset of X.

If B(K) is a linear space, then for every R > 0, B R (K) is dense in X * .
Proof of Lemma 3.3. Let us suppose that, for some R > 0, B R (K) is not dense in X * . Then, take f in X * and ε > 0 such that

(f + B(0, ε)) ∩ B R (K) = ∅.
Using the previous Lemma, we observe that

f + B(0, ε) ⊂ X * = B(K) = B R (K) ∪ B(K) \ B R (K),
and therefore, as

B(K) \ B R (K) = B(K) ∩ x∈K {f ∈ X * : f, x ≤ R} is convex and closed, f + B(0, ε) ⊆ B(K) \ B R (K). (15)
Accordingly, the linear space B(K) has a nonempty interior, from where it follows that B(K) = X * ,which by the principle of uniform boundedness implies that K is bounded, a contradiction.

We have proved that for every element f ∈ X * and for every R > 0, there exists a sequence (

f n ) n∈N * in B R (K) which converges to f and satisfies σ K (f n ) > R.
This result has an immediate consequence.

Lemma 3.4. Let f be in B(K). In the assumptions of Lemma 3.3, for every constant R, satisfying R > σ K (f ), and every γ > 0, there is a sequence

(g n ) n∈N * converging to f , such that σ K (g n ) = R, and g n , x < R, ∀x ∈ (K ∩ B(0, γ)), ∀n ∈ N * .
Proof of Lemma 3.4. According to Lemma 3.3, there is a sequence (f n ) n∈N * converging to

f such that R ≤ σ K (f n ) < +∞.
As, for every n, the function λ

→ σ K (λf + (1 -λ)f n ) is continuous on [0, 1], there is λ n ∈ [0, 1] such that σ K (g n ) = R,
where

g n = λ n f + (1 -λ n )f n .
Denoting by {λ n } n∈N a converging subsequence of {λ n } n∈N , the sequence {g n } n∈N obviously converges to f . Moreover, at least starting from a certain rank, this subsequence also satisfies the second condition of the lemma. Indeed, otherwise we would find a subsequence (g m

) m∈N * of (g n ) n∈N * such that x m ∈ K, x m ≤ γ and g m , x m = R. Consequently, |R -f, x m | = | g m -f, x m | ≤ γ f -g m ,
and, letting m → +∞, we would obtain

R > σ K (f ) ≥ lim m→∞ f, x m = R,
a contradiction and the result follows.

We have now all the ingredients which are necessary to prove the main result of this section.

Proof of Theorem 3.1. Let us first remark that, if B(K) is not a linear space, then Theorem 3.1 is established. Indeed, in this case, there is f ∈ B(K) such that -f / ∈ B(K), and h = f /σ K (f ) verifies (14), as inf w∈K h, w = -∞.

Let us now consider the case where B(K) is a linear space. Without loss of generality we may assume that 0 ∈ K. In order to define the element h of X * , we define by induction a sequence (h

n , γ n , ε n , x n ) n∈N * ⊂ X * × R × R × K as follows.
For n = 0, let us put

h 0 = 0, γ 0 = 1, ε 0 = 1, x 0 ∈ K arbitrary. For n = 1, take an element h 1 of X * satisfying sup x∈K h 1 , x = 3 4 and h 1 , x < 3 4
, ∀x ∈ K ∩ B(0, 1), (Lemma 3.4 applied for f = 0, R = 3/4 and γ = 1, ensures the existence of such an element). Consequently, there is some x 1 ∈ K such that h 1 , x 1 = 2/3; take γ 1 = max(2, x 1 ). Finally, set

ε 1 = min 3 4 -sup x∈K∩ B(0,1) h 1 , x , 1 6γ 1 . (16)
Let us now suppose that the sequence was defined for each i, 1 ≤ i ≤ n -1 in such a way that the following relations hold for every 1 ≤ i ≤ n -1:

h i -h i-1 < ε i-1 , (17) sup x∈K h i , x = i + 2 i + 3 and h i , x < i + 2 i + 3 , ∀x ∈ K ∩ B(0, γ i-1 ), h i , x i = (i + 1)/(i + 2), γ i = max(γ i-1 + 1, x i ), ε i = min i + 2 γ i-1 (i + 3) - sup x∈K∩ B(0,γ i-1 ) h i , x γ i-1 , 1 (i + 1)(i + 2)γ i , ε i-1 -h i -h i-1 .
Lemma 3.4, applied for f = h n-1 , R = (n + 2)/(n + 3) and γ = γ n-1 , yields the existence of an element h n of X * such that

h n -h n-1 < ε n-1 ,
and

sup x∈K h n , x = n + 2 n + 3 and h n , x < n + 2 n + 3 , ∀x ∈ K ∩ B(0, γ n-1 ).
Consequently, there exists

x n in K such that h n , x n = n + 1 n + 2 .
Set γ n = max(γ n-1 + 1, x n ), and put ε n for the following (strictly positive) expression:

min n + 2 γ n-1 (n + 3) - sup x∈K∩ B(0,γ n-1 ) h n , x γ n-1 , 1 (n + 1)(n + 2)γ n , ε n-1 -h n -h n-1 .
The sequence (h n , γ n , ε n , x n ) n∈N defined inductively satisfies the relations (17). The last relation in (17) implies that

h i + B(0, ε i ) ⊂ h i-1 + B(0, ε i-1 ), ∀i ≥ 2;
since γ n ≥ n + 1, we deduce from the previous relation that, for every m > n ,

h m -h n < ε n < n + 2 γ n-1 (n + 3) < 1 n . (18)
Relation (18) means that the sequence (h n ) n∈N * is a Cauchy sequence in X * and therefore converges to some h ∞ ∈ X * . By relation (18), we derive that

h n -h ∞ ≤ ε n ≤ 1 γ n-1 n + 2 n + 3 - sup x∈K∩ B(0,γ n-1 )
h n , x .

It follows that, for every x ∈ K ∩ B ¯(0, γ n-1 ),

h ∞ , x = h n , x + h ∞ -h n , x ≤ (19) ≤ h n , x + x • h ∞ -h n ≤ h n , x + γ n-1 1 γ n-1 n + 2 n + 3 - sup x∈K∩ B(0,R n-1 ) h n , x = n + 2 n + 3 + h n , x - sup x∈K∩ B(0,γ n-1 ) h n , x ≤ n + 2 n + 3 . Let x be in K; there exists n ∈ N * such that x ≤ γ n-1 . Hence relation (19) implies h ∞ , x ≤ n + 2 n + 3 < 1, ∀x ∈ K, (20)
and therefore sup

x∈K h ∞ , x ≤ 1. ( 21
)
Again from relation (18) we deduce that

h n -h ∞ ≤ ε n ≤ 1 (n + 1)(n + 2)γ n , that is h ∞ , x n = h n , x n + h ∞ -h n , x n (22)
≥ n + 1 n + 2 -h n -h ∞ x n ≥ n + 1 n + 2 - 1 (n + 1)(n + 2) = n n + 1 . Consequently σ K (h ∞ ) ≥ lim n→∞ h ∞ , x n = 1. (23)
Combining ( 20), ( 21) and ( 23), it follows that, for any

x ∈ K, h ∞ , x < 1 = σ K (h ∞ ).
The mapping h defined by

h =    h ∞ if inf w∈K h ∞ , w = -∞ h ∞ • min{-1, inf x∈K h ∞ , x } -1 if inf w∈K h ∞ , w > -∞ (24)
satisfies the conclusion of Theorem 3.1.

The main result

We state the following result; its proof will be presented afterwards.

Theorem 4.1. Let Φ be a Γ 0 (X)-functional which achieves its minimum value on X.

Assume that either (a) Ker Φ ∞ is not a linear subspace; or (b) Φ is non-semicoercive and Ker Φ ∞ is a linear space.

Then, for every ε > 0, there exists Φ ε ∈ Γ 0 (X), such that

Φ(x) -ε ≤ Φ ε (x) ≤ Φ(x), ∀x ∈ X,
and argmin Φ ε = ∅.

The result which follows is an immediate consequence of Theorem 4.1 and can be considered as the main result of the paper.

Theorem 4.2 (Main Result).

Let Φ be a Γ 0 (X)-functional. Suppose that Φ and every small uniform perturbation of Φ (in the class Γ 0 (X)) achieve its minimum value on X. Then, Φ is necessarily semicoercive.

The proof given below is based on preceding results established in Sections 2 and 3.

Proof of Theorem 4.1. Case (a): In this case, the functional Φ cannot be semicoercive. Let x 0 ∈ argmin Φ and M ∈ R be such that Φ(x 0 ) = M . Since Φ is bounded from below, Φ ∞ (x) ≥ 0 for every x in X. It follows from ( 2) and (3) that for every v ∈ Ker Φ ∞ , we have

Φ(x + tv) ≤ Φ(x), ∀ t ≥ 0 and x ∈ X. (25)
Since Ker Φ ∞ is not a linear subspace of X, there exists v in X such that Φ ∞ (v) = 0 and Φ ∞ (-v) > 0. Let σ : [0,1) → R be a convex and increasing function satisfying σ(0) = 0 and

lim t→1 -σ(t) = +∞. ( 26
)
For example take σ(t) = t 1t . We then define Φ ε by

Φ ε (x) = inf 0≤t<1 [Φ(x -σ(t)v) -εt] . ( 27
) 0 x M X R Φ Φ M- ε ε Φ-ε
Let us prove that Φ ε meets the requirements of Theorem 4.1.

Lemma 4.3. The functional Φ ε defined in (27) is a Γ 0 (X)-functional.

Proof of Lemma 4.3. (i) The functional Φ ε is an extended real-valued functional. In fact, since M is the minimum of Φ, we have

Φ(x -σ(t)v) -εt > M -ε, ∀x ∈ X, ∀t ∈ [0, 1).
By taking the infimum, in the previous relation, over t ∈ [0, 1), we obtain

Φ ε (x) ≥ M -ε > -∞, ∀x ∈ X. (28)
(ii) The functional Φ ε is convex. Let x 1 , x 2 be in X, t 1 , t 2 ∈ [0, 1) and 0 ≤ λ ≤ 1. Using the convexity of σ and relation (25) for

x = λ(x 1 -σ(t 1 )v) + (1 -λ)(x 2 -σ(t 2 )v), and t = λσ(t 1 ) + (1 -λ)σ(t 2 ) -σ(λt 1 + (1 -λ)t 2 )
(remark that t ≥ 0), we obtain,

Φ x 2 + λ(x 1 -x 2 ) -σ(t 2 + λ(t 1 -t 2 ))v = Φ(x + tv) (29) ≤ Φ(x) = Φ λ(x 1 -σ(t 1 )v) + (1 -λ)(x 2 -σ(t 2 )v) .
Combining the convexity of Φ and (29), we derive

Φ x 2 + λ(x 1 -x 2 ) -σ(t 2 + λ(t 1 -t 2 ))v ≤ (30) λΦ x 1 -σ(t 1 )v + (1 -λ)Φ x 2 -σ(t 2 )v .
Using ( 27) and (30), we derive for every t 1 , t 2 ∈ [0, 1):

Φ ε (λx 1 + (1 -λ)x 2 ) ≤ λΦ(x 1 -σ(t 1 )v) + (1 -λ)Φ(x 2 -σ(t 2 )v) -ε(λt 1 + (1 -λ)t 2 ).
Taking the infimum over t 1 and t 2 , in the previous relation yields

Φ ε (λx 1 + (1 -λ)x 2 ) ≤ λΦ ε (x 1 ) + (1 -λ)Φ ε (x 2 ).
Hence Φ ε is convex.

(iii) The functional Φ ε is lower semicontinuous. Let (x n ) n∈N * , be a sequence in X such that x n → x 0 as n → +∞. Consider also a sequence (t n ) n∈N * in [0, 1) such that lim

n→∞ t n = t * . If t * = 1, relation (26) implies that lim n→∞ σ(t n ) = +∞. Hence, v + (x 0 -x n ) σ(t n ) → v as n → +∞.
We may therefore apply (4) for -v,v + (x 0x n ) σ(t n ) and σ(t n ), and obtain

0 < Φ ∞ (-v) ≤ lim inf n→∞ Φ(x n -σ(t n )v) σ(t n ) . ( 31
)
As the sequence (σ(t n )) n∈N * tends to infinity, from (31) it follows that

lim inf n→∞ [Φ(x n -σ(t n )v) -εt n ] = +∞. Hence, lim inf n→∞ [Φ(x n -σ(t n )v) -εt n ] ≥ Φ ε (x 0 ). ( 32
)
If t * < 1, as Φ is lower semicontinuous , we have

lim inf n→∞ [Φ(x n -σ(t n )v) -εt n ] ≥ Φ(x 0 -σ(t * )v) -εt * ≥ Φ ε (x 0 ). ( 33
)
Relations ( 32) and (33) imply that for every sequence

(t n ) n∈N * , t n ∈ [0, 1) lim inf n→∞ [Φ(x n -σ(t n )v) -εt n ] ≥ Φ ε (x 0 ).
According to the definition of Φ ε , for every n in N * there is

t n in [0, 1) such that 0 ≤ Φ(x n -σ(t n )v) -εt n -Φ ε (x n ) ≤ 1 n . Consequently, lim inf n→∞ Φ ε (x n ) = lim inf n→∞ [Φ(x n -σ(t n )v) -εt n ] ≥ Φ ε (x 0 ),
that is, Φ ε is lower semicontinuous. Hence, this completes the proof of Lemma 4.3.

Lemma 4.4. The functional Φ ε satisfies:

Φ -ε ≤ Φ ε ≤ Φ.
Proof of Lemma 4.4. Let x be in X. As σ(0) = 0, using (27) we observe that

Φ ε (x) ≤ Φ(x -σ(0)v) -ε 0 = Φ(x). ( 34
)
Since σ is positive, relation (25) implies

Φ(x) -ε ≤ Φ(x -σ(t)v) -ε ≤ Φ(x -σ(t)v) -εt, ∀ t ∈ [0, 1).
By taking the infimum over t we derive

Φ(x) -ε ≤ Φ ε (x). ( 35
)
The conclusion of the Lemma 4.4 follows immediately by summing up (34) and ( 35).

Lemma 4.5. argmin Φ ε = ∅.

Proof of Lemma 4.5. Let x 0 be in argmin Φ. For every 0 < δ < 1, we have

Φ ε (x 0 + σ(1 -δ)v) ≤ Φ(x 0 + σ(1 -δ)v -σ(1 -δ)v) -ε(1 -δ) = M -ε + δε. Hence, inf x∈X Φ ε (x) ≤ M -ε.
Relation (28) implies now that the infimum of Φ ε is Mε. Suppose that this infimum is reached, i.e. there is x in X such that

Φ ε (x) = inf 0≤t<1 [Φ(x -σ(t)v) -εt] = M -ε; accordingly, there is a sequence (t n ) n∈N * in [0, 1) such that lim n→∞ [Φ(x -σ(t n )v) -εt n ] = M -ε. Since Φ(x -σ(t n )v) ≥ M ,
Φ(x -σ(t n )v) = M. Consequently, 0 < Φ ∞ (-v) = lim n→∞ Φ(x -σ(t n )v) σ(t n ) = 0,
a contradiction. Hence the functional Φ ε does not reach its infimum value.

Lemmata 4.3-4.5 show that, if Ker Φ ∞ is not a linear subspace of X, then the functional defined by (27) fulfills the conditions of Theorem 4.1. This thereby completes the proof of Case (a).

Case (b):

The construction of functionals of type Φ y,N allows us to deal with the case when Ker Φ ∞ is a linear (and closed) subspace of X. Let us consider first a particular case, namely when Ker Φ ∞ = {0}. The desired functional Φ ε will be in this case of the form

Φ ε (x) = max{Φ 1 (x), Φ 2 (x)}, (36)
where Φ 1 is defined by

Φ 1 (x) = Φ x 0 ,M -ε . 0 Φ ε Φ 1 Φ 2 x Φ M R X M-ε By Lemma 2.1, Φ 1 ∈ Γ 0 (X) and Φ(x) -ε ≤ Φ 1 (x) ≤ Φ(x), ∀ x ∈ X. (37)
In order to define Φ 2 , let us consider the closed and convex set

P = {x ∈ X : Φ 1 (x 0 + x) ≤ M }. (38)
Lemma 4.6. The closed convex set P is unbounded, nevertheless it is linearly bounded.

Proof of Lemma 4.6. Suppose P is bounded, i.e. there is a positive constant k such that x ≤ k for all x ∈ P. Pick x in Dom Φ and set

y = x 0 + ε Φ(x) -M + ε (x -x 0 ); since (y, M ) = ε Φ(x) -M + ε (x, Φ(x)) + 1 - ε Φ(x) -M + ε (x 0 , M -ε), we have (y, M ) ∈ co{(x 0 , M -ε), epi Φ}. Hence Φ 1 (y) ≤ M , that is (y -x 0 ) ∈ P . Thus ε Φ(x) -M + ε x -x 0 = y -x 0 ≤ k. Consequently, k ε Φ(x) ≥ x -x 0 + k M -ε ε ,
which means that Φ is coercive, a contradiction. Therefore P is unbounded.

Fix now u ∈ P ∞ ; then for every positive constant s we have su ∈ P . Relation (37) combined with (38) yields

Φ(x 0 + su) ≤ Φ 1 (x 0 + su) + ε ≤ M + ε.
Hence due to (4) Φ ∞ (u) = 0, that is u ∈ Ker Φ ∞ . Therefore u = 0 and, consequently,

P ∞ = {0}.
According to Theorem 3.1, take an element f ∈ X * such that

inf w∈P f, w < f, u ≤ 1, ∀ u ∈ P. (39) Define
Φ 2 (x) := ε 2 f, x -x 0 + M - 3ε 4 . (40)
It easily follows from the definition of Φ 1 and Φ 2 that the functional Φ ε defined in (36) belongs to Γ 0 (X), and satisfies Φ(x)ε ≤ Φ ε (x). In order to prove that Φ ε (x) ≤ Φ(x), we use the following result.

Lemma 4.7. Let x in X. If Φ 2 (x) > M then Φ 1 (x) > Φ 2 (x).
Proof of Lemma 4.7. Take

λ 0 = M -Φ 2 (x 0 ) Φ 2 (x) -Φ 2 (x 0 )
. We have

λ 0 Φ 2 (x) + (1 -λ 0 )Φ 2 (x 0 ) = Φ 2 (λ 0 x + (1 -λ 0 )x 0 ) = M. (41)
Since,

Φ 2 (x) = ε 2 f, x -x 0 + M - 3ε 4 ≤ ε 2 + M - 3ε 4 = M - ε 4 < M, ∀ x ∈ (x 0 + P ),
we deduce from (41) that

(λ 0 x + (1 -λ 0 )x 0 ) / ∈ (x 0 + P ).
Lemma 4.8 will be established if we show that Φ 2 does not reach its infimum value on B.

To this end, we prove that we can rewrite Φ 2 on B as

Φ 2 (x) = M -ε + ε 2 • 1 2 -f, T (x) , ( 44
)
where T is defined by

T (x) = ε(x -x 0 ) Φ 2 (x) -M + ε = 4(x -x 0 ) 1 + 2 f, x -x 0 .
Lemma 4.10. The operator T : B → P is one-to-one and onto.

0 + P x 0 + T(x) Φ ε x . M R X M-ε B x .
Proof of Lemma 4.10. Let us first prove that T is well-defined. As T (x 0 ) = 0 ∈ P, it suffices to show that T (x) ∈ P if x lies in B and x = x 0 .

By Lemma 2.1 (iv), argmin Φ 1 = {x 0 } and thus Φ

1 (x) > Φ 1 (x 0 ) = M -ε. Let us set s 0 = ε Φ 1 (x) -M + ε . As Φ 1 (x) ≤ Φ 2 (x)
, according to Lemma 4.7 we deduce that Φ 2 (x) ≤ M . Therefore,

s 0 ≥ ε Φ 2 (x) -M + ε = 4 1 + 2 f, x -x 0 . Suppose Φ 1 (x 0 + s 0 (x -x 0 )) > M . As the function [0, s 0 ] ∋ s → Φ 1 (x 0 + s(x -x 0 ))
is convex and lower semicontinuous , hence continuous, there is s 1 in [0, s 0 ) such that Φ 1 (x 0 + s 1 (xx 0 )) = M . Applying Lemma 2.1 (iii) with Φ x 0 ,M -ε for Φ, x 0 + s 1 (xx 0 ) for x and s 1 /(1 + s 1 ) for λ, and using the convexity of Φ 1 we obtain

M = Φ 1 (x 0 + s 1 (x -x 0 )) ≤ Φ 1 (x 0 ) + s 1 (Φ 1 (x) -Φ 1 (x 0 )), that is s 1 ≥ ε Φ 1 (x) -M + ε = s 0 , a contradiction. Consequently, Φ 1 (x 0 + s 0 (x -x 0 )) ≤ M. (45)
For simplicity, let us define

α 1 = s 0 - 4 1 + 2 f, x -x 0 s 0 and α 2 = 1 - s 0 - 4 1 + 2 f, x -x 0 s 0 .
The following relation

x 0 + T (x) = α 1 x 0 + α 2 (x 0 + s 0 (x -x 0 ))
together with the convexity of Φ 1 and (45) imply that

Φ 1 (x 0 + T (x)) ≤ α 1 (M -ε) + α 2 s 0 M ≤ M.
Hence, T (x) belongs to P . Now fix w in P and set y = x 0 + w,

λ 0 = 3 -2 f, y -x 0 4 -2 f, y -x 0 = 3 -2 f, w 4 -2 f, w ,
(remark that relation (14) implies 0 < λ 0 < 1), and x = λ 0 x 0 + (1λ 0 )y. After straightforward calculations we deduce that

Φ 2 (x) = M -ε + ε 2 • 1 2 -f, y -x 0 = M -ελ 0 . (46)
Since w ∈ P , we have Φ 1 (y) ≤ M . Lemma 2.1 (iii) implies that

Φ 1 (x) = λ 0 Φ 1 (x 0 ) + (1 -λ 0 )Φ 1 (y) ≤ M -ελ 0 = Φ 2 (x).
Thus x belongs to B, and since

T (x) = 4(x -x 0 ) 1 + 2 f, x -x 0 = 4(1 -λ 0 ) 1 + 2 f, x -x 0 w = w,
it follows that T is onto. Since the operator T is obviously one-to-one, this completes the proof of Lemma 4.10.

Let us return to the proof of Lemma 4.

f, w > -∞, then Φ 2 (x) = M -ε + ε 2 • 1 2 -f, T (x) (47) > M -ε + ε 2 • 1 2 -inf w∈P f, w = inf v∈B Φ 2 (v), ∀ x ∈ B, and that, if inf w∈P f, w = -∞, then Φ 2 (x) = M -ε + ε 2 • 1 2 -f, T (x) > M -ε = inf v∈B Φ 2 (v), ∀ x ∈ B.
In both cases, we have proved that Φ 2 does not reach its infimum value on B, which allows us to complete the proof of Theorem 4.1.

Let us return to the general case, where Ker Φ ∞ is an arbitrary closed subspace of X. As a consequence of (25), we have Φ(x + v) = Φ(x), for all x ∈ X, and all v ∈ Ker Φ ∞ .

We may therefore factorize X by Ker Φ ∞ ; the quotient functional Φ is a non-semicoercive Γ 0 (X/ Ker Φ ∞ )-functional which attains its minimum value, and satisfies Ker Φ ∞ = {0}. We may thus define Φ ε as before, and set Φ ε (x) = Φ ε (x), for all x ∈ x, and all x ∈ X.

This functional obviously satisfies the requirements of Theorem 4.1.

Following the lines of Theorem 4.1 we can also derive the following result:

Theorem 4.11. Suppose that Φ is a Γ 0 (X)-functional which achieves its minimum value on X. Moreover, assume that either (a) Ker Φ ∞ is not a subspace; or (b) Φ is non-semicoercive and Ker Φ ∞ is a linear space.

Then, for every ε > 0 and R > 0, there exists Φ ε,R ∈ Γ 0 (X) such that

• argmin Φ ε,R = ∅; • for each x ∈ X, Φ(x) -ε ≤ Φ ε,R (x) ≤ Φ(x); • if u ∈ argmin Φ ε,R , then u ≥ R.
Proof of Theorem 4.11. Case (a): If Ker Φ ∞ is not a linear subspace of X, then there is v in Ker Φ ∞ such that

-v / ∈ Ker Φ ∞ . Let δ : = dist(-v, Ker Φ ∞ ) and set y : = x 0 + R + x 0 δ v and define Φ ε,R (x) : = Φ y,M -ε (x); by virtue of Lemma 2.1, Φ ε,R lies in Γ 0 (X) and Φ(x) -ε ≤ Φ ε,R (x) ≤ Φ(x), ∀ x ∈ X. Moreover, if z ∈ argmin Φ ε,R , then z ∈ (y + Ker Φ ∞ ). Thus there is w ∈ Ker Φ ∞ such that z = x 0 + R + x 0 δ v + w Since δ = dist(-v, Ker Φ ∞ ), we have R + x 0 δ v + w ≥ R + x 0 δ -v - δ R + x 0 w ≥ R + x 0 ; accordingly, z ≥ R + x 0 δ v + w -x 0 ≥ R.
Hence Φ ε,R fulfills the conditions of Theorem 4.11.

Case (b): If Ker Φ ∞ is a linear subspace of X, consider Φ ε as constructed in the proof of Theorem 4.1, and let (x n ) n∈N * be a minimizing sequence for Φ ε . Since the functional Φ ε does not attain its infimum value and the space X is reflexive, the sequence (x n ) n∈N * is unbounded.

According to relation (47), it follows that

M -ε ≤ inf x∈X Φ ε ≤ M - 3ε 4 .
Thus take n 0 such that x n 0 ≥ R and Φ ε (x n 0 ) < M .

Consider now Φ xn 0 ,Φ ε (xn 0 ) . Lemma 2.1 implies that Φ xn 0 ,Φ ε (xn 0 ) is a Γ 0 (X/ Ker Φ ∞ )functional which reaches its minimum value only on x n 0 and satisfies

Φ xn 0 ,Φ ε (xn 0 ) (x) ≤ Φ(x). (49)
Since (x n 0 , Φ ε (x n 0 )) ∈ epi Φ ε , and Φ ε ≤ Φ, we obtain epi Φ xn 0 ,Φ ε (xn 0 ) = co{(x n 0 , Φ ε (x n 0 )), epi Φ} ⊂ epi Φ ε .

Therefore Φ -ε ≤ Φ ε ≤ Φ xn 0 ,Φ ε (xn 0 ) . (50)
The desired functional is now defined by setting Φ ε,R (x) = Φ xn 0 ,Φ ε (xn 0 ) (x), ∀ x ∈ x, ∀ x ∈ X/ Ker Φ ∞ .

By (49) and (50), it follows that Φ ε,R ∈ Γ 0 (X) and satisfies:

Φ(x) -ε ≤ Φ ε,R (x) ≤ Φ(x), ∀ x ∈ X.
The functional Φ ε,R reaches its minimum value only on the set x n 0 . Since x n 0 ≥ R, we have x ≥ R for every x in x n 0 , which means that Φ ε,R fulfills all the conditions of Theorem 4.11.

Concluding remarks

In this paper we have studied the stability under uniform perturbations of the existence of a solution for the simplest variational problem, namely the minimization of a proper, convex and lower semicontinuous functional. In summary, we established that the problem of finding a minimum point of a convex functional is stable under uniform perturbations only within the class of semicoercive Γ 0 (X)-functionals.

The same question may be raised in some other variational contexts, such as the theory of variational inequalities.

This case has already been considered and several results establishing sufficient existence conditions for noncoercive problems have been obtained recently, using the so-called recession analysis (see for instance the work by Adly et al [START_REF] Adly | Recession mappings and noncoercive variational inequality[END_REF] and Attouch et al [START_REF] Attouch | Recession operators and solvability of variational problems[END_REF]). However, even in the case of a positive operator, the question of the existence of solutions of the perturbed initial variational inequality remains partially open.

Another interesting direction of research is the nonconvex case. We remark that this case does not subsume the convex one, since, even if the class of functionals is broader, so is the uniform neighborhood composed of functionals which must attain their minima. Accordingly, no simple relation can be established between the two problems.

Sufficient stability conditions have been obtained in the non-convex setting (see [START_REF] Auslender | Noncoercive optimization problems in reflexive Banach spaces[END_REF] and [START_REF] Baiocchi | General existence results for unilateral problems in continuum mechanics[END_REF]), implying, inter alia, that semicoercive functionals are no longer the only functionals with a stable minimum. The problem of characterizing all lower semicontinuous functionals having a stable minimum is thus still open.

  and εt n ≤ ε, the previous relation implies that

	lim n→∞	t n = 1, and lim n→∞

by the Région Limousin under a PDZR European Grant.

Therefore, Φ 1 (λ 0 x + (1λ 0 )x 0 ) > M.

Since Φ 2 (x 0 ) = M -(3/4)ε < M (see (40)), we obtain 0 < λ 0 ≤ 1. As Φ 1 is convex, we have

Substracting (41) from (42) yields

and establishes Lemma 4.7.

By Lemma 4.7, Φ 2 (x) ≤ Φ(x) whenever Φ 2 (x) > M ; in other words,

This allows us to conclude that Φ ε (x) ≤ Φ(x).

We conclude the proof of the particular case Ker Φ ∞ = {0} by proving the following statement.

Lemma 4.8. The functional Φ ε does not attain its infimum value.

Proof of Lemma 4.8. Define B := {x ∈ X : Φ 2 (x) ≥ Φ 1 (x)} and let us show that, if Φ ε attains its minimum value over X, then Φ 2 also reaches its minimum value over B.

Lemma 4.9. For every x in X, there is b(x) in B such that Φ ε (x) ≥ Φ 2 (b(x)).

Proof of Lemma 4.9. If x is such that Φ ε (x) ≥ M , then b(x) := x 0 satisfies Lemma 4.9; and if x belongs to B, we set b(x) := x. It remains to define b(x) when M > Φ ε (x) and Φ 1 (x) > Φ 2 (x). In this case, set

and the conclusion of Lemma 4.9 follows by taking b(x) = λ 0 x 0 + (1λ 0 )x.


Matteo Ciocca
email: matteo.ciocca@inria.fr

Pierre-Brice Wieber

Thierry Fraichard
email: thierry.fraichard@inria.fr

Strong Recursive Feasibility in Model Predictive Control of Biped Walking

Realizing a stable walking motion requires satisfying a set of constraints. Model Predictive Control (MPC) is one of few suitable methods to handle such constraints. The capacity to satisfy constraints, which is usually called feasibility, is classically guaranteed recursively. In our applications, an important aspect is that the MPC scheme has to adapt continuously to the dynamic environment of the robot (e.g. collision avoidance, physical interaction). We aim therefore at guaranteeing recursive feasibility for all possible scenarios, which is called strong recursive feasibility. Recursive feasibility is classically obtained by introducing a terminal constraint at the end of the prediction horizon. Between two standard approaches for legged robot, in our applications we favor a capturable terminal constraint. When the robot is not really planning to stop and considers actually making a new step, recursive feasibility is not guaranteed anymore. We demonstrate numerically that recursive feasibility is actually guaranteed, even when a new step is added in the prediction horizon.

I. INTRODUCTION

Walking depends on contact forces between the feet and the ground. The unilateral nature [START_REF] Wieber | Viability and predictive control for safe locomotion[END_REF] of this interaction (the feet can only push on the ground) limits the motion that a legged robot can realize and plays a crucial role in its stability. In the case of walking on a flat ground, this corresponds to having the Center of Pressure (CoP) stay within the support polygon [START_REF] Herdt | Walking without thinking about it[END_REF]. Realizing a stable walking motion requires satisfying a set of constraints [START_REF] Wieber | Modeling and control of legged robots[END_REF]. Model Predictive Control (MPC) is one of few suitable methods to handle such constraints [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. It has been used therefore extensively for the control of legged robots [START_REF] Wieber | Model predictive control[END_REF], [START_REF] Gouaillier | Omni-directional closedloop walk for NAO[END_REF]. For example, the MPC scheme in [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF] generates a walking motion online with automatic footstep placement. It was expanded to ensure safe navigation in a crowd [START_REF] Bohorquez | Safe navigation strategies for a biped robot walking in a crowd[END_REF] and physical collaboration with humans [START_REF] Agravante | Walking pattern generators designed for physical collaboration[END_REF].

MPC solves an optimal control problem over a prediction horizon. The solution is a control sequence that satisfies a set of constraints. The first element of the sequence is applied to the system. The whole process is then repeated. The capacity to satisfy constraints, which is usually called feasibility, is classically guaranteed recursively [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. An important aspect of our applications is that the MPC scheme has to adapt continuously to the dynamic environment of the robot: collision avoidance [START_REF] Bohorquez | Safe navigation strategies for a biped robot walking in a crowd[END_REF], physical interaction with humans [START_REF] Agravante | Walking pattern generators designed for physical collaboration[END_REF], or visual feedback [START_REF] Dune | Vision based control for Humanoid robots[END_REF]. We aim therefore at guaranteeing recursive feasibility for all possible scenarios, which is called strong recursive feasibility [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF].

Recursive feasibility is classically obtained by introducing a terminal constraint at the end of the prediction horizon [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF] to make sure the system remains feasible indefinitely after the end of the horizon. Two standard approaches for legged robots are (i) to consider that the robot keeps repeating indefinitely the same cyclic motion [START_REF] Takenaka | Real time motion generation and control for biped robot: Walking gait pattern generation[END_REF], [START_REF] Tajima | Fast running experiments involving a humanoid robot[END_REF], [START_REF] Nagasaka | Integrated motion control for walking, jumping and running on a small bipedal entertainment robot[END_REF] or (ii) that it stops after a given number of steps [START_REF] Takanishi | Dynamic biped walking stabilized with optimal trunk and waist motion[END_REF], [START_REF] Lim | Online walking pattern generation for biped humanoid robot with trunk[END_REF], [START_REF] Bohorquez | Safe navigation strategies for a biped robot walking in a crowd[END_REF], what corresponds to capturability [START_REF] Pratt | Capture point: A step toward humanoid push recovery[END_REF], [START_REF] Koolen | Capturability-based analysis and control of legged locomotion, part 1: Theory and application to three simple gait models[END_REF]. Passive safety requires that the robot is able to stop before any collision occurs [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF]. For that reason we favor a capturable terminal constraint as in [START_REF] Bohorquez | Safe navigation strategies for a biped robot walking in a crowd[END_REF].

We can see in Fig. 1 how such a capturability terminal constraint makes sure that the system remains feasible indefinitely. This way, when the prediction horizon advances as in Fig. 2, we are sure that the MPC scheme remains feasible: strong recursive feasibility is guaranteed by construction. The problem is when the robot is not really planning to stop, and considers actually making a new step. With such a sudden change, recursive feasibility is not guaranteed anymore, as shown in Fig. 3. It seems that this issue has been overlooked in the literature of MPC for legged robots. Note that this issue is unrelated to the length of the prediction horizon. The objective of this paper is to investigate this issue and provide a numerical evidence that recursive feasibility is actually guaranteed, even when a new step is added in the prediction horizon.

Outline of the Paper: Section II introduces the definition of strong recursive feasibility in MPC. Section III describes how to provide a numerical evidence of strong recursive feasibility. The dynamics of walking is described in Section IV. Section V describes where we apply our numerical approach. Results of our numerical approach are presented in Section VI.

II. STRONG RECURSIVE FEASIBILITY IN MODEL PREDICTIVE CONTROL Consider a linear time-invariant discrete-time system

x i+1 = Ax i + Bu i , (1)
where x i , u i and x i+1 are respectively the state, control and successor state, and (A, B) are matrices of proper dimensions. For a state x i , an MPC sheme aims to find a control sequence U i = {u (i|i) , u (i+1|i) , . . . , u (i+N |i) } [1] that satisfies a set of constraints

E i x i + F i U i ≤ d i , (2)
where (E i , F i ) are time-varying matrices and d i is a timevarying vector of proper dimensions. Let the set of solutions be

W i {(x i , U i ) | E i x i + F i U i ≤ d i } . (3)
An MPC scheme classically chooses a sequence U i that minimizes an objective function [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], and its first element u (i|i) = κ is applied to the system. Two useful projections of W i are therefore:

(W κ ) i {(x i , κ) | ∃ U i s.t. (x i , U i ) ∈ W i ∧ u (i|i) = κ} , (4) and
X i {x i | ∃ U i s.t. (x i , U i ) ∈ W i } . (5)
With these sets, we can define strong recursive feasibility in the following way: Definition 1 (Strong Recursive Feasibility [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF]): The MPC scheme is strongly recursive feasible if and only if ∀i

∀ (x i , κ) ∈ (W κ ) i , x i+1 = A i , x i + B i κ ∈ X i . (6)
The sets W, W κ and X and Definition 1 are represented in Figs. 4 and5. By definition, these sets are convex polytopes [START_REF] Boyd | Convex optimization[END_REF]. In our case, they also happen to be closed. Thanks to these properties, it is sufficient to check property (6) only on the vertices of W κ , as can be seen on Fig. 5.

III. NUMERICAL EVIDENCE OF STRONG RECURSIVE FEASIBILITY

The number of vertices of these polytopes is finite, but enumerating all of them from definition (3) and ( 4) is actually an NP-hard problem [START_REF] Borwein | Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity[END_REF]. Various algorithms and approaches exist, e.g. [START_REF] Löfberg | Oops! I cannot do it again: Testing for recursive feasibility in MPC[END_REF], [START_REF] Fukuda | Double description method revisited[END_REF], [START_REF] Jones | Equality set projection: A new algorithm for the projection of polytopes in halfspace representation[END_REF]. We are going to use here a simple randomized approach. The method we propose is based on the observation that the simplex method for Linear Programming (LP) always terminate on a so-called basic solution, which is actually a vertex [START_REF] Bertsimas | Introduction to linear optimization[END_REF]. We propose therefore to solve LPs of the form

min xi,Ui γ T x i U i s.t. E i F i x i U i ≤ d i , (7)
for randomly chosen directions γ, which will provide a random selection of vertices of W i . Each vertex found in this way is projected in (W κ ) i and property ( 6) is checked simply verifying that the following LP

min Ui+1 ψ T U i+1 s.t. E i+1 (Ax i + Bκ) + F i+1 U i+1 ≤ d i+1 (8)
has a solution for any direction ψ, chosen arbitrarily.

IV. DYNAMICS OF WALKING

The horizontal motion of the Center of Mass (CoM) c ∈ R 2 is linearly related to the Center of Pressure (CoP) p ∈ R 2 when walking on a flat ground with constant height [START_REF] Wieber | Modeling and control of legged robots[END_REF]:

c = ω 2 (c -p) (9)
where ω 2 = g/h, h is the height of the CoM above the ground and g is the norm of the gravity vector. This assumes zero rate of change of centroidal angular momentum. Since contact forces with the ground are unilateral, the CoP is always constrained within the support polygon P

[27]: p -s j ∈ P , (10)
where s j ∈ R 2 is the j th footstep on the ground. The CoM position is constrained to a closed convex region C due to the maximal leg length of the robot [START_REF] Brasseur | A robust linear MPC approach to online generation of 3d biped walking motion[END_REF]:

c -s j ∈ C . (11)
Biped robots should not cross their legs while walking. The position of the (j + 1) th footstep with respect to the position of the j th footstep is therefore restricted to a region S j where the legs do not cross:

s j+1 -s j ∈ S j . ( 12
)
The robot is said to be 0-step capturable [START_REF] Koolen | Capturability-based analysis and control of legged locomotion, part 1: Theory and application to three simple gait models[END_REF] when it can stop without having to make any further step. We enforce this situation at the end of the prediction horizon by introducing the following terminal constraint:

ξ N = p N , ξN = 0 , (13)
where ξ ∈ R 2 is the Capture Point defined as:

ξ = c + ω -1 ċ . ( 14
)
Note that the constraints ( 10)-( 11) and ( 12) are time varying due to the introduction of the step position s j .

V. PROBLEM FORMULATION

In both horizontal coordinates x and y, the motion of the CoM c of a legged robot is commonly modeled as a triple integrator [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF]:

ĉi+1 =   1 T T 2 /2 0 1 T 0 0 1   ĉi +   T 3 /6 T 2 /2 T   ... c i , (15)
where T is the sampling time, ĉi = (c i , ċi , ci ) and ĉi+1 = (c i+1 , ċi+1 , ci+1 ) are two consecutive states, and ... c i is the CoM jerk. Thus, the CoM motion state is

x i = ĉx i ĉy i ∈ R 6 . ( 16
)
By evolving the state (16) N -times following the dynamics (15), we obtain:

Ĉ = Āc x i + Bc (U ... c ) i (17)
which relates the sequence Ĉ ∈ R 6N of N CoM position in x and y with a sequence U ... c ∈ R 2N of N CoM jerk in x and y. The same state evolution can be applied to the relation ( 9):

P = Āp x i + Bp (U ... c ) i ( 18
)
where P ∈ R 2N describes the evolution of the CoP position and ( Āc , Bc ) and ( Āp , Bp ) are described in the Appendix VIII. Similar relationship is applied to the selection of the footsteps as

S i = V c s c + V f (U s ) i ( 19
)
where s c ∈ R 2 is the current footstep on the ground, U s ∈ R 2m is a sequence of future footstep positions and (V c , V f ) are cyclic time-varying matrices (see Appendix VIII) that determine which foot is on the ground at what time. Equations ( 17)-( 18)-( 19) are used to formulate the time-varying constraints (10) ∧ (11) ∧ (12) ∧ (13) [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF] in the linear form as [START_REF] Herdt | Walking without thinking about it[END_REF] where

U i = (U ... c ) i (U s ) i T . ( 21
)
The set of constraints is used to investigate a numerical evidence that the recursive feasibility is actually guaranteed during the critical time transition when a new step appears at the end of the predicted horizon. We indicate the discrete time instants before and after the appearance of a new predicted step at the end of the finite horizon respectively t and t + 1. This critical time transition is shown for time t in Fig. 1 and for t + 1 in Fig. 3.

VI. EVALUATION OF RANDOMLY SELECTED VERTICES

The parameters of our legged robot were selected accordingly with the kinematics of the robot HRP-2 [START_REF] Kaneko | Design of prototype humanoid robotics platform for HRP[END_REF]. The constraints set ( 10)-( 11)-( 12) were defined with respect to the current footstep position, which in our case can be chosen arbitrarily: we chose s c = s t = (0, 0). The vertices in (W κ ) t were found with the linear programming problem [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF] with the set of constraints (E t , F t , d t ). Instead of projecting each solution W t onto (W κ ) t , it was decided to set to zero the entries of γ that do not multiply x t or κ and to choose randomly the others. Each new vertex was determined to be unique or a duplicate of the ones found before. The finite number of vertices at time t were used to check property [START_REF] Gouaillier | Omni-directional closedloop walk for NAO[END_REF] with the LP (8) with the set of constraints (E t+1 , F t+1 , d t+1 ).

A. Results

The choice of robot's model and constraint parameters is summarized in Table I. Three millions vertices were found with [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF]. Within the set of 3 million vertices, many of them were duplicates. Without repetitions, the total number of vertices was 180. In Fig. 6, after ≈1.8 million randomly found vertices, we found only duplicates of the 180 vertices. Thus, the search of the finite number of vertices in (W κ ) t was concluded. Successively, the LP (8) had a solution for all 180 vertices with the constraints at time t + 1. MATLAB R2016a was used to run both linear programming problems ( 7)-( 8) with the linprog function (simplex method).

VII. CONCLUSIONS

A terminal constraint at the end of the prediction horizon guarantees recursive feasibility [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. We favor a capturable terminal constraint to ensure passive safety of the robot in a dynamic environment [START_REF] Bohorquez | Safe navigation strategies for a biped robot walking in a crowd[END_REF]. In the case we do not want the robot to stop, we have to change the capturable terminal constraint region periodically in order to consider new steps. And with such sudden change of the terminal constraint, recursive feasibility is not guaranteed anymore.

In this paper we proposed a numerical evidence that shows recursive feasibility is actually guaranteed even when a new step is added in the prediction horizon. In our case the set of constraints for walking is linear and defines a close convex polytope, Section IV. Thanks to these properties we verified property (6) only for the vertices of the polytope when a new step is added in the prediction horizon. Despite the sudden change of the terminal region set, by construction the MPC for legged robot is actually strongly recursive feasible.

Aspects that can lead to further research are: (i) compare our method with one of the algorithms in the literature, e.g. the double description method [START_REF] Fukuda | Double description method revisited[END_REF], (ii) investigate whether the MPC scheme loses the strong recursive feasible property due to an inappropriate choice of one (or more) of these three parameters: the length of the finite horizon N , the number of predicted steps m or the duration of each predicted step t s . For example, the MPC problem could fail to find a solution when the change of the terminal constraint region is close to the present time (small N ) or when it is suddenly on the next predicted step (m = 1).

VIII. APPENDIX

From Section V, the matrices ( Āc , Bc ) are:

Āc =      A A 2
. . .

B = B 0 B 0 , B 0 =   T 3 /6 T 2 /2 T   (24)
and

C = C 0 0 0 C 0 , C 0 = 1 0 - 1 ω 2 . ( 26
)
And

V 0 =                  1 . . . 1 0 . . . 0 0 . . . 0                  ∈ R N ×1 , V 1 =                   0
                  ∈ R N ×m (27)
are cyclic time-varying matrices that compose (V c , V f ) as follows:

V c = V 0 0 0 V 0 , V f = V 1 0 0 V 1 ( 28
)
in which the ones determine which foot is on the ground at what time [START_REF] Herdt | Walking without thinking about it[END_REF].

Fig. 1 :

 1 Fig.1:A capturability terminal constraint (that defines a terminal constraint set X f ) makes sure that the motion of the legged robot remains feasible indefinitely.

Fig. 2 :Fig. 3 :

 23 Fig.2: When the prediction horizon advances, thanks to the fixed terminal constraint set X f , the MPC scheme remains feasible.

Fig. 4 :

 4 Fig. 4: The polytope W, defined in (3), and its projection W κ defined in (4).

Fig. 5 :

 5 Fig. 5: Strong Recursive Feasibility (Def. 1).

Fig. 6 :

 6 Fig.6: In a set of 3 million randomly selected vertices (x t , κ) ∈ (W κ ) t , only 180 were unique.

A

  N -1 B A N -2 B . . . B

TABLE I :

 I Robot's Parameters

	Parameter	Symbol	Value	Unit
	Height of the CoM	h	0.8	[m]
	Feet dimesions	(f l , fw) (0.24, 0.14)	[m]
	Leg length	L	0.3	[m]
	Feet separation	δ	0.2	[m]
	Step duration	ts	0.8	[s]
	Finite Horizon Dimension	N	16	-
	Number of predicted steps	m	2	-
	Sampling time	T	0.1	[s]

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK,

Grenoble, France pierre-brice.wieber@inria.fr

the double subscript notation (i + k|i) indicates a prediction k steps ahead from time i

ACKNOWLEDGMENT

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01), PSPC Romeo 2 project and EU H2020 Comanoid Research and Innovation Action (RIA). The authors would like to thank Stéphane Caron for the fruitful discussions.


The linear buncher of SPIRAL beam test of a prototype