Constantin-Irinel Morărescu
email: constantin.morarescu@imag.fr

Bernard Brogliato
email: bernard.brogliato@inrialpes.fr

Joseph Fourier

Passivity-based switching control of flexible-joint complementarity mechanical systems $

Keywords:

In this study one considers the tracking control problem of a class of nonsmooth Lagrangian systems with flexible joints and subject to frictionless unilateral constraints. The task under consideration consists of a succession of free-motion and constrained-motion phases. Particular attention is paid to impacting and detachment phases. A passivity-based switching controller that allows one to extend the stability analysis described in our previous works to the case of systems with lumped flexibilities, is proposed. Numerical tests show the effectiveness of the controller.

Introduction

The control of mechanical systems with unilateral constraints is a challenging problem, see e.g. [START_REF] Galeani | Trajectory tracking for a particle in elliptical billiards[END_REF]; Leine and van de Wouw (2008a,b). In this work we extend the trajectory tracking control framework developed in [START_REF] Brogliato | On the control of finitedimensional mechanical systems with unilateral constraints[END_REF], [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF], Morărescu and Brogliato (2008) considering the multiconstraints flexible-joint case. The typical task representation and the stability analysis are more complex than in the rigid case because the ''free-motion'' controller originates from a passivitybased backstepping scheme proposed in [START_REF] Brogliato | Global tracking controllers for flexiblejoint manipulators: A comparative study[END_REF], and the study of the variations of the Lyapunov function V (•) arising from the backstepping method is more difficult than in the rigid case. Indeed V (•) cannot be forced to match the system's total energy at the impact times, and suitable upperbounds have to be computed to lead the stability analysis. The dynamics in (1) is a simplified dynamics obtained from more general Lagrangian systems using a generalized coordinate $ The material in this paper was partially presented at European Control Conference, 23-26 August 2009, Budapest, Hungary. This paper was recommended for publication in revised form by Associate Editor Zongli Lin under the direction of the Editor Andrew R. Teel. transformation as in [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF], that is supposed to hold globally in the configuration space. Notice that a nonlinear stiffness K Z(q, θ) may appear due to the transformation. Details on can be found in Morărescu, Brogliato, and Haad (2008). The class of systems under consideration is expressed as:

      
M(q)q + C (q, q)q + G(q) + K (qθ) = D λ J θ + K (θq) -K Z(q, θ) = U q 1 ≥ 0, (q 1 ) T λ = 0, λ ≥ 0 Collision rule (1) where q ∈ R n is the vector of rigid link angles, θ ∈ R n is the vector of motor shaft angles, M(q) = M T (q) ∈ R n×n is a positive definite inertia matrix, C (q, q) is the matrix of Coriolis and centripetal forces, G(q) contains conservative forces, λ ∈ R m is the vector of contact forces (or Lagrangian multipliers) associated with the constraints, J ∈ R n×n is the diagonal and constant matrix of the actuator inertia, K = K > 0, K ∈ R n×n represents the stiffness matrix, U ∈ R n is the vector of generalized torque inputs, and The time-derivative of a function f (•) is denoted by ḟ (•). For any function f (•) the limit to the right at the instant t will be denoted by f (t + ) and the limit to the left will be denoted by f (t -) when they exist. A simple jump of the function f (•) at the moment t = t is denoted σ f (t ) = f (t + ) -f (t -). For a real valued function f : R + → R one denotes by S(f ) the set of all real valued functions g : R + → R such that there exists a positive real constant 0 < c < ∞ satisfying g(t) ≤ cf (t), ∀t ≥ 0. One writes g ∈ S(1) ≡ L ∞ if f (t) = 1, ∀t ≥ 0. 0 n is the n-vector with entries 0, and 0 n×m is the n × m-zero matrix. I m is the m × m identity matrix A vector is considered positive if all its components are positive. A Linear Complementarity Problem (LCP) with unknown λ is a system: λ ≥ 0, Aλ + b ≥ 0, λ (Aλ + b) = 0, which is compactly rewritten as 0 ≤ λ ⊥ Aλ + b ≥ 0. Such an LCP has a unique solution for all b if and only if A is a P-matrix. Positive definite matrices (not necessarily symmetric) are P-matrices.

q 1 = Dq ∈ R m with D = [I m 0 m×(n-m) ]. A constraint i is said to be active if q 1 i = 0, and inactive if q 1 i > 0. General notations
Definition 1. Let ≥ 0 be a fixed real number. We say that a pimpact occurs at the instant t if

q 1 i i∈I (t) ≤ , i∈I q 1 i (t) = 0 where I ⊂ {1, . . . , m}, card(I) = p. If = 0 all p surfaces Σ i = ∂Φ i = {q ∈ R n
| q 1 i = 0}, i ∈ I are struck simultaneously. When > 0 the system collides ∂Φ in a neighborhood of the intersection i∈I Σ i .

A collision (or restitution) rule is a relation between the postimpact and the pre-impact velocities. Among the various models of collision rules, Moreau's rule is an extension of Newton's law which is energetically consistent [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF] and is numerically tractable [START_REF] Acary | Lecture notes in applied and computational mechanics[END_REF]. In the special coordinates of (1) this reads as q1 i (t + ) = -eq 1 i (t -) when q 1,i (t) = 0 and q1 i (t -) ≤ 0, where i ∈ {1, . . . , m} and e ∈ [0, 1]. Under mild conditions on the data, the solutions are such that positions q(•) and θ(•) are absolutely continuous functions of time, whereas the generalized velocity is right continuous of local bounded variation [START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF]. The continuity of θ (•) holds [START_REF] Brogliato | Nonsmooth mechanics[END_REF] and will be used in the stability analysis developed in Section 5.

Basic concepts

Since the system's dynamics does not change when the number of active constraints decreases one gets the following typical task representation:

R + = Ω ∅ 0 ∪ k≥1 I k ∪ m k i=1 Ω B k,i k B k,m k ⊂ B k+1,1 , B k,m k ⊂ B k,m k -1 ⊂ • • • ⊂ B k,1 (2)
where the superscript B k represents the set of active constraints phases 1 when the number of active constraints increases. We note that B k = ∅ corresponds to free-motion. When the number of active constraints decreases there is no impact, thus no other transition phases are needed. Throughout the paper, the sequence

(B k = {i ∈ {1, . . . , m} | F i (X) = 0})
I k ∪ m k i=1 Ω B k,i k
will be referred to as the cycle k of the system's evolution.

System properties

For kinematic chains with prismatic or revolute joints the following properties hold.

1 To simplify the presentation the phases of motion and the time intervals that correspond to the phases are denoted identically.

Property 1. The matrix d dt M(q) -2C (q, q) is skew-symmetric and Ṁ(q) d dt M(q) = C (q, q) + C (q, q). Furthermore the matrix C (q, q) is a smooth function of q and q with the well-known properties C (q, q) ∈ S( q ) and C (q, y)z = C (q, z)y, ∀ q, y, z ∈ R n .

Property 2. The conservative forces vector G(q) is such that ∂G(q) ∂q ∈ S(1) which implies by the mean value theorem G(q 1 ) -

G(q 2 ) ∈ S( q 1 -q 2 ), ∀q 1 , q 2 ∈ R n .
Property 3. The matrix C (q, q) is such that ∂C(q,q) ∂q ∈ S( q ) and ∂C(q,q) ∂ q ∈ S(1).

Design of the exogenous trajectories

Throughout the paper, the following trajectories will play a role in the closed-loop dynamics:

• q nc (•) denotes the desired trajectory that the system should track if there were no constraints. We suppose that q 1,nc (t) < 0 for some t, otherwise the problem reduces to the tracking control of a system with no constraints.

• q * d (•) denotes the signal entering the control input and playing the role of the desired trajectory during some parts of the motion.

• q d (•) represents the signal entering the Lyapunov function V (•). This signal is set on the boundary ∂Φ after the first impact of each cycle.

We always consider

I k = [τ k 0 , t k f ]
where τ k 0 is chosen by the designer as the start of the transition phase I k and t k f is the end of this phase. During the transition phases the system must be stabilized on the intersection of some surfaces Σ i . This will be done by mimicking the behavior of a ball falling on the ground under gravity. Therefore all the components except the ones that are normal to the constraints belonging to B k will be frozen. Moreover for robustness reasons one avoids a tangential approach and imposes some impacts defining an exogenous signal q * d that violates the constraints.

In order to simplify the presentation we introduce the following notations (where all superscripts (•) k will refer to the cycle k of the system motion):

• t k 0 is the first impact during the cycle k, • t k ∞ is the accumulation point of the sequence {t k } ≥0 of the impact instants during the cycle k (t k f ≥ t k ∞ ),
• τ k 1 will be explicitly defined later and represents the instant when the exogenous signal q * d reaches a given value chosen by the designer in order to impose a closed-loop dynamics with impacts during the transition phases,

• t k,i d is the desired detachment instant at the end of Ω J k,i k . It is noteworthy that t k 0 , t k ∞ are state-dependent whereas τ k 1 , τ k 0 and t k,i
d are exogenous and imposed by the designer.

Stability analysis criteria

Let us define Ω as the complement in R + of I = k≥0 I k and assume that the Lebesgue measure of Ω, denoted λ[Ω], equals infinity. Let x(•) be the state of the closed-loop system in (1) with some feedback controller U(q, q, θ, θ, t). Definition 2. The closed loop system is called practically weakly stable if there exists 0 < R < +∞ and t * < +∞ such that x(t) < R for all t > t * , t ∈ Ω.

Consider V (•) such that there exist strictly increasing functions α(•) and β(•) satisfying the conditions: α(0) = 0, β(0) = 0 and α( x ) ≤ V (x, t) ≤ β( x ). The following result was proved in Morărescu et al. (2008):

Proposition 1. Assume that the task admits the representation (2) and that

(a) t k f < t k,1 d , ∀k ∈ N, (b) outside the impact accumulation phases [t k 0 , t k ∞ ] one has V (x(t), t) ≤ -γ V (x(t), t) for some constant γ > 0, (c) the system is initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, (d) V (t k ∞ ) ≤ ρ * V (τ k 0 ) + ξ where ρ * , ξ ∈ R + . Then V (τ k 0 ) ≤ δ(γ , ξ ), ∀k ≥ 1 where δ(γ , ξ
) is a function that can be made arbitrarily small by increasing either the value of γ or the length of the time interval [t ∞ , t f ]. Thus, the system is practically weakly stable with R = α -1 (δ(γ , ξ )).

It is worth to point out that condition (d) means that the impacts may be considered as a kind of disturbance that can be suitably upper bounded. This is certainly the most crucial point in Proposition 1.

Controller design

Let us introduce the notations:

q θ = q -q d θ -θ d , s = s 1 s 2 = q + γ 2 q
θ + γ 2 θ , qr = qdγ 2 q, q = qq * d and s 1 = q + γ 2 q, where γ 2 > 0 is a scalar gain. The tracking problem is solved using a generalization of the controller proposed in Brogliato et al. (1995, Equ. (28)) and the closed-loop stability analysis of the system is based on Proposition 1. The controller is defined by

U = J θr + K (θ d -q d ) -γ 1 s 2 -K Z(q, θ ) θ d = q d + K -1 U r (3)
where U r is given by:

U r =                  U ∅ c U nc = M(q)q r + C (q, q)q r + G(q) -γ 1 s 1 for t ∈ Ω ∅ 2k U B k c = U nc -P d + K f (P q -P d ) for t ∈ Ω B k k U B k c for t ∈ I k before the first impact U B k t = M(q)q r + C (q, q)q r + G(q) -γ 1 s 1 for t ∈ I k after the first impact (4)
where γ 1 > 0 is a scalar gain, K f > 0, P q = D T λ and P d = D T λ d is the desired contact force during the persistently constrained motion. It is clear that during Ω B k k not all the constraints are active and, therefore, some components of λ and λ d are zero.

Notice that on impacting phases no force feedback is applied.

Also U is a function of q, θ , q, θ only (no acceleration feedback).

The rationale behind the change of structure of U r after the first impact, is that it facilitates the calculation of some upper-bounds which are necessary to recast the closed-loop stability analysis into Proposition 1 (see Section 5 and the Appendix). In order to prove the stability of the closed-loop system (1), ( 3) and (4) we will use the following positive definite function:

V (t, s, q, θ ) = 1 2 s T 1 M(q)s 1 + 1 2 s T 2 Js 2 + γ 1 γ 2 qT q + γ 1 γ 2 θT θ + 1 2 (q -θ ) T K (q -θ ).
(5)

Tracking control strategy

The design of q * d (•) and q d (•) during the phases I k , the design of the desired contact force during constraint phases and the strategy for take-off at the end of constraint phases are similar to those proposed in the rigid case (see Morărescu & Brogliato, 2008). Therefore, we briefly present the results and some discussions without giving the proof details.

During the impacting transition phases the system must be stabilized on ∂Φ. Obviously, this does not mean that all the constraints have to be activated (i.e. q 1 i (t) = 0, ∀i = 1, . . . , m). Let us consider that only the first p constraints (eventually reordering the coordinates) define the border of Φ where the system must be stabilized. The signal q * d (•) will be then defined as follows:

• choosing ν > 0 and denoting t =

t-τ k 0 τ k 1 -τ k 0 , the components q i d * , i = 1, . . . , p of q *
d p are defined as:

q i d * (t) = a 3 (t ) 3 + a 2 (t ) 2 + a 0 , t ∈ [τ k 0 , min{τ k 1 ; t k 0 }] -νV 1/3 (τ k 0 ), t ∈ (min{τ k 1 ; t k 0 }, t k f ] (6)
where V (•) is defined in (5) and τ k 1 is chosen by the designer such that the limit condition

q i d * (τ k 1 ) = -νV 1/3 (τ k 0 ), qi d *
(τ k 1 ) = 0 holds, which allows the computation of the previous coefficients as:

a 3 = 2[ q i nc (τ k 0 ) + νV 1/2 (τ k 0 )] a 2 = -3[ q i nc (τ k 0 ) + νV 1/2 (τ k 0 )] a 0 = q i nc (τ k 0 ) (7)
• all the other components of q * d (•) are frozen on I k , i.e.:

q * d n-p (t) = q nc n-p (τ k 0 ), t ∈ (τ k 0 , t k f ] (8)
In order to limit the deformation of the desired trajectory q * d (•) w.r.t. q nc (•) during the I k phases, we impose in the sequel

q nc p (τ k 0 ) ≤ ν 1 (9)
where ν 1 > 0 is chosen by the designer. A smaller ν 1 leads to a smaller deformation of the desired trajectory and to a smaller deformation of the real trajectory. Nevertheless, due to the tracking error, ν 1 cannot be chosen to be zero. We also note that ( 9) is a practical way to choose τ k 0 .

During the transition phases

I k we define (q d ) n-p (t) = q * d n-p (t).
Assuming a finite accumulation period, the impact process can be considered in some way equivalent to a plastic impact. Therefore, (q d ) p (•) and (q d ) p (•) are set to zero on the right of t k 0 . It is worth recalling that the first impact time t k 0 of each cycle k, is unknown.

   M(q)q + F = (1 + K f )D p (λ -λ d ) J ṡ2 + γ 1 s 2 + K ( θ -q) = 0 0 ≤ q p ⊥ λ p ≥ 0 (10)
where F = F (q, q, q, q, θ ) = -M(q)q r +C(q, q)s 1 +γ 1 s 1 +K (qθ)

and D p = [I p . . . O p×(n-p) ] ∈ R p×n . On Ω B k
k the system has to be permanently constrained which is equivalent to q p (•) = 0 and qp (•) = 0. In order to assure these conditions it is sufficient to have λ p > 0.

We denote

M -1 (q) = [M -1 (q)] p,p [M -1 (q)] p,n-p [M -1 (q)] n-p,p [M -1 (q)] n-p,n-p
and C (q, q) = C (q, q) p,p C (q, q) p,n-p C (q, q) n-p,p C (q, q) n-p,n-p where the meaning of each component is obvious. Let us also denote by K p the matrix made of the first p rows and p columns of K .

Proposition 2. On Ω B k k the constraint motion of the closed-loop system (10), ( 3) and ( 4) is assured if the desired contact force is defined by

(λ d ) p ν p + K p θp 1 + K f - Mp,p (q) 1 + K f [M -1 (q)] p,p C p,n-p (q, q) + [M -1 (q)] p,n-p (C n-p,n-p (q, q) + γ 1 I n-p ) (s 1 ) n-p where Mp,p (q) = [M -1 (q)] p,p -1 = D p M -1 (q)D T p -1 and ν p ∈ R p , ν p > 0.
Let us now consider that the first h constraints (h < p) have to be deactivated. Thus, the detachment takes place at

t k d if qh (t k+ d ) > 0 which requires λ h (t k- d ) = 0. The last p -h constraints remain active which means λ p-h (t k- d ) > 0.
To simplify the notation we drop the time argument in many equations of this section. We decompose the LCP matrix as:

(1

+ K f )D p M -1 (q)D T p = A 1 (q) A 2 (q) A 2 (q) T A 3 (q) (11) with A 1 ∈ R h×h , A 2 ∈ R h×(p-h) and A 3 ∈ R (p-h)×(p-h) .
Proposition 3. The closed-loop system (10), ( 3) and (4) is perma-

nently constrained on [t k f , t k d ) and a smooth detachment is guaran- teed on [t k d , t k d + ) ( is a small positive real number chosen by the designer) if (λ d ) h (t k d ) (λ d ) p-h (t k d ) = A 1 -A 2 A -1 3 A T 2 -1 b h -A 2 A -1 3 b p-h -C 1 (t -t k d ) C 2 + A -1 3 b p-h -A T 2 (λ d ) h where b p b(q, q, U ∅ c ) -D p M -1 (q)F ≥ 0 and C 1 ∈ R h , C 2 ∈ R p-h such that C 1 ≥ 0, C 2 > 0. Furthermore on [t k d , t k d + ) q * d (t) = q d (t) = q * h (t) q nc n-h (t)
,

where q * h (•) is a twice-differentiable function such that q * h (t k d ) = 0, q * h (t k d + ) = q nc h (t k d + ), q * h (t k d ) = 0, q * h (t k d + ) = qnc h (t k d + ) (12) and q * h (t k+ d ) = a > max 0, -A 1 (q)(λ d ) h (t k- d ) .

Closed-loop stability analysis

For convenience V (t, s(t), q, θ(t)) is written as V (t).

Lemma 1. Consider the closed-loop system (1), ( 3) and (4) with (q * d ) p (•) defined on the interval [τ k 0 , t k 0 ] as in (6)-( 8). Let us also suppose that condition (b) of Proposition 1 is satisfied. The following inequalities hold:

q(t k- 0 ) ≤ V (τ k 0 ) γ 1 γ 2 , s 1 (t k- 0 ) ≤ 2V (τ k 0 ) λ min (M(q)) , θ(t k- 0 ) ≤ V (τ k 0 ) γ 1 γ 2 , s 2 (t k- 0 ) ≤ 2V (τ k 0 ) λ min (J) , q(t k- 0 ) ≤ 2 λ min (M(q)) + γ 2 γ 1 V 1/2 (τ k 0 ) θ(t k- 0 ) ≤ 2 λ min (J) + γ 2 γ 1 V 1/2 (τ k 0 ). Furthermore, if t k 0 ≤ τ k 1 one has (q d ) p (t k- 0 ) ≤ + V (τ k 0 ) γ 1 γ 2 , (q d ) p (t k- 0 ) ≤ k + k * V 1/6 (τ k 0 ) (q d ) p (t k- 0 ) ≤ 6 √ 2 q nc p (τ k 0 ) + √ pνV 1/2 (τ k 0 ) (q (3) d ) p (t k- 0 ) ≤ 6 √ 2 q nc p (τ k 0 ) + √ pνV 1/2 (τ k 0 ) (13)
where is the real constant fixed in Definition 1,

k = 6 √ pν 1 τ k 1 -τ k 0 and k * = 6 √ p τ k 1 -τ k 0 1 √ γ 1 γ 2 + ν (ν + ν 1 ) + ν.
Proof. The proof closely follows the ideas presented in Morărescu and Brogliato (2008) and is omitted. The details can be found in Morărescu et al. (2008).

Lemma 2. The controller U in (3) and (4)

guarantees that σ θ (•) , σ˙θ (•) ∈ S(1) ≡ L ∞ .
Proof. See Appendix A.

Theorem 1. Let e = 0 and q * d (•) defined as in (6)-( 8). The closedloop system (1), ( 3) and (4) initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, satisfies the requirements of Proposition 1 and is therefore practically weakly stable with the closed-loop state x(

•) = [q, θ (•), s(•)] and R = e -γ (t k f -t k
∞ ) (ρ * + ξ )/ ρ where ρ * , ρ and ξ are defined in the proof.

Proof. See Appendix B.

Illustrative example

We simulate the behavior of a planar two-link flexible-joint manipulator in presence of two constraints. The admissible domain is Φ = {(x, y) | y ≥ 0, 0.7 -x ≥ 0}. The unconstrained desired trajectory is given by the circle {(x, y) | (x -0.7) 2

+ y 2 = 0.5} and violates both constraints. In other words, the two-link planar manipulator must track a quarter-circle; stabilize on and then follow the line Σ 1 = {(x, y) | y = 0}; stabilize on the intersection of Σ 1 and Σ 2 = {(x, y) | x = 0.7}; detach from Σ 1 and follow Σ 2 until the unconstrained circle re-enters Φ and finally take-off from Σ 2 in order to repeat the previous steps. Therefore, we have:

R + = Ω ∅ 0 ∪I 1 ∪Ω J 1,1 1 ∪I 2 ∪Ω J 2,1 2 ∪Ω J 2,2 2 ∪Ω J 2,3 2 ∪I 3 ∪Ω J 3,1 2 ∪I 4 ∪• • • with J 1,1 = {1}, J 2,1 = {1, 2}, J 2,2 = {2}, J 2,3 = ∅, etc.
We note that during I 2k+1 the system is stabilized on Σ 1 (1-impacts) while during I 2k the system is stabilized on Σ 1 ∩ Σ 2 (2 k -impacts). The numerical values used for the dynamical model are l 1 = l 2 = 0.5 m, m 1 = m 2 = 1 kg, I 1 = I 2 = 0.5 kg m 2 , J 1 = J 2 = 0.1 kg m 2 and the impacts are imposed by ν = 10 in ( 6) and ( 7). The stiffness matrix is defined by K = diag (2000 N m, 2000 N m). Let us say that the quarter-circle is completely tracked in one round. We set the period of each round to 10 s and we simulate the dynamics during 6 rounds using the Moreau's time-stepping algorithm of the Siconos software platform [START_REF] Acary | Lecture notes in applied and computational mechanics[END_REF]. We set the controller gains γ 1 = 10, γ 2 = 1 and we choose ν 1 = 0.1 (like this we implicitly set τ k 0 see ( 9)) in order to better point out the deformation of q d (•) on the transition phases (Fig. 1 (left)). In

Conclusions

In this paper we have proposed a solution for the control of nonsmooth Lagrangian systems with flexible joints. All the ingredients entering the dynamics (desired trajectories, desired contact forces, exogenous instants playing a role in the definition of the control law) are explicitly defined. Numerical simulations are done with the software platform in order to illustrate the results.

On the other hand s (t k+ 0 )M(q)s(t k+ 0 ) -s (t k- 0 )M(q)s(t k- 0 ) = σ s 1 M(q)s 1 (t k 0 ) + σ s 2 Js 2 (t k 0 ). As proved in Morărescu and Brogliato (2008) there exists a real positive constant c 3 such that σ s 1 M(q)s 1 (t k 0 ) + γ 1 γ 2 σ q q(t k 0 ) ≤ c 3 , ∀k ≥ 0. 3) and ( 4) is practically weakly stable with R = α -1 e -γ (t k f -t k ∞ ) (1 + ξ ) . Let us consider ρ = min{λ min (M(q))/2; γ 1 γ 2 }. Defining α : R + → R + , α(ω) = ρω 2 we get α(0) = 0, α( [s(t), q(t)] ) ≤ V (t, s, q) and the proof is finished.

  and definitions. • is the Euclidean norm, b p ∈ R p and b n-p ∈ R n-p are the vectors formed with the first p and the last np components of b ∈ R n , respectively. λ min (•) and λ max (•) represent the smallest and the largest eigenvalues, respectively.

  during the corresponding motion phase Ω k , and I B k k denotes the transient between two Ω k

  The desired contact force P d = D λ d must be designed such that it is large enough to assure the constraint motion on the Ω B k kphases. Some contact force components have also to be decreased at the end of the Ω B k k -phases in order to allow the detachment. Therefore we need a lower bound of the desired force which assures both contact (without any undesired detachment which can generate other impacts) during the Ω B k k phases and a smooth detachment at the end of Ω B k k . Dropping the time argument, the dynamics of the system on Ω B k k can be written as

Fig. 1 .Fig. 2 .

 12 Fig. 1. Left: The trajectory of the system during 6 rounds; Right: The variation of the almost non-increasing Lyapunov function during the first round.

Fig. 2

 2 Fig. 2 we have shifted the desired trajectory backwards on I 3 to highlight that the Lyapunov function at the instant τ k 0 is smaller when k increases. The behavior of the system during one round is emphasized in Fig. 1 (right).

(

  2 assures the existence of c4 ∈ R + such that γ 1 γ 2 σ θ θ (t k 0 ) ≤ c 4 , ∀k ≥ 0. (B.5)In conclusion, inserting (B.2)-(B.5) in (B.1) one getsV (t k+ 0 ) -V (t k- 0 ) ≤ c 1 + c 2 + c 3 + c 4 , ∀k ≥ 0. (B.6) Thus condition (d) of Proposition 1 is verified for ρ * = 1, ξ = c 1 + c 2 + c 3 + c 4and the closed-loop system (1), (

  It is easy to see thatσ s 2 Js 2 (t k 0 ) = 2s 2 (t k- 0 )Jσ s 2 (t k 0 ) + σ s 2 (t k 0 )Jσ s 2 (t k 0 )and using Lemma 1, Lemma 2 and the relation σ s 2 (t k 0 ) = σ˙θ (t k 0 ) + γ 2 σ θ (t k 0 ) one deduces that there exists a real positive constant c 2 such that σ s 2 Js 2 (t k 0 ) ≤ c 2 , ∀k ≥ 0.

	(B.3)

Acknowledgement

The first author's work was supported by the European project SICONOS FP5 IST-2001-37172.

Appendix A. Proof of Lemma 2

Since θ(•), θ(•) are continuous on R + and θ d (•), θd (•) are continuous on R + \ {t k 0 | k ∈ Z} one deduces that σ θ (t) = 0 = σ˙θ (t), ∀t = t k 0 . Therefore Lemma 2 holds if there exist some real constants that upper bound σ θ (t k 0 ) , σ˙θ (t k 0 ) , ∀k ∈ Z. The definition of θ d (•) (see (3)) allows us to write

Using a (4) one obtains

Moreover, as we have mentioned at the end of Section 3, (q d ) p (•), (q d ) p (•) and implicitly (q d ) p (•) are set to zero on (t k 0 , t k f ]. Thus taking into account that q(t k+ 0 ) ≤ w q(t k- 0 ) (where w = λ max (M) λ min (M) ) one arrives at

are bounded by some constants. Thus all the quantities in (A.2) are bounded by some constants independent of the cycle index k. This means that σ U r (t k 0 ) is bounded by a constant independent of the cycle index, which implies the same for σ θ (t k 0 ) . In other words σ θ (t) ∈ S(1). Differentiating (4) one obtains Ur (t) = M(q)q (3) r (t) + Ṁ(q)q r (t) + C (q, q)q r (t)

where Ṁ, Ċ stand for dM dt and dC dt respectively. It is clear that Ċ(q, q)(t) = ∂C ∂q (q, q)q(t) + ∂C ∂ q (q, q)q(t)

and using Properties 1 and 3 one derives Ċ(q, q)(t) ∈ S( q(t) 2 + q(t) ).

Furthermore, Lemma 1 and the first equation in (1) assure that q(t) 2 , q(t) ∈ S(1). Thus Ċ(q, q)(•) , σĊ (q,q) (•) ∈ S(1) and one obtains σĊ (q,q) (t k 0 )q r (t k 0 ) ≤ σĊ (q,q) (t k 0 ) . qr (t k+ 0 )

+ Ċ(q, q)(t k- 0 ) . σ qr (t k 0 ) ∈ S(1).

(A.4)

Property 1 allows us to replace Ṁ(q) by C (q, q) + C (q, q) which leads to Ṁ(q)q r (t) + C (q, q)q r (t) ≤ 3 C (q, q) . qr (t)

⇒ Ṁ(q)q r (t) + C (q, q)q r (t) ∈ S( q . qr (t) ).

Since qr (t) ≤ qd (t) + γ 2 q(t) , using Lemma 1 one gets Ṁ(q)q r (t) + C (q, q)q r (t) ∈ S(1).

(A.5)

The definitions ( 6)-( 8) and the first equation in (1) assure that q

(3)

Property 2 states that ∂G ∂q ∈ S(1), which implies ∂G ∂q q(t) ∈ S( q(t) ) 1). (b) Using that Ṁ(q) -2C (q, q) is a skew-symmetric matrix (see Property 1), straightforward computations show that on R + \ k≥0 [t k 0 , t k f ] the time derivative of the Lyapunov function is given by

where we have used the fact that (q d ) p ≡ 0, (q d ) p ≡ 0, q p ≡ 0, qp ≡ 0, thus (s 1 ) p ≡ 0 on constraint phases and λ p ≡ 0, (λ d ) p ≡ 0 on free-motion phases. On the other hand

(d) There is only one impact during each transition phase since e = 0 and with the choice of U B t in (4). Therefore

. We compute now the jump of the Lyapunov function at the impact time t k

Replacing ψ(t k+ 0 ) by ψ(t k- 0 ) + σ ψ (t k 0 ), the second term of the right hand side of (B.1) becomes 1 2 2 ψ, (t k- 0 )Kσ ψ (t k 0 ) + σ , ψ (t k 0 )Kσ ψ (t k 0 )

which is upper bounded by

Therefore Lemmas 1 and 2 imply that there exists a real positive constant c 1 such that 1 2 ψ, (t k+ 0 )K ψ(t k+ 0 ) -ψ, (t k- 0 )K ψ(t k- 0 ) ≤ c 1 , ∀k ≥ 0.

(B.2)


Laetitia Laguzet

High order variational numerical schemes with application to Nash -MFG vaccination games

Keywords: Variational scheme, Mean Field Games, Nash equilibrium, Vaccination games, SIR model

This paper introduces high-order explicit Runge-Kutta numerical schemes in metric spaces. We show that our approach reduces to the corresponding Runge-Kutta schemes if the ambient space is Hilbert.

We apply these schemes to compute the Nash equilibrium in a Mean Field vaccination Game. Numerical simulations show improvement in the speed of convergence towards the Nash equilibrium; the numerical scheme has high order (two to four) in time.

Introduction 1.State of the art and motivation

The games with a continuum of agents have been widely studied during the last decade thanks to the Mean Field Games (MFG) theory, introduced by Lasry-Lions ( [START_REF] Lasry | Jeux à champ moyen. I: Le cas stationnaire[END_REF][START_REF] Lasry | Jeux à champ moyen. II: Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF]) and Huang-Caines-Malhamé ( [START_REF] Huang | Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle[END_REF][START_REF] Huang | An invariance principle in large population stochastic dynamic games[END_REF][START_REF] Huang | Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized -Nash equilibria[END_REF]). An application of this framework is the modeling of the vaccination decision in a population where an individual has the choice between only two pure strategies: to get vaccinated or not.

This decision is a part of the behavioral epidemiology, which takes into account the interplay between human behavior and spread of infectious diseases (see [START_REF] Manfredi | Modeling the interplay between human behavior and the spread of infectious diseases[END_REF] for a detailed presentation and [START_REF] Ibuka | Free-riding behavior in vaccination decisions: An experimental study[END_REF][START_REF] Verelst | Behavioural change models for infectious disease transmission: a systematic review (2010-2015)[END_REF][START_REF] Nardin | Planning horizon affects prophylactic decision-making and epidemic dynamics[END_REF] for other examples). In particular, risk perception and spread of information (on vaccination risk or probability of infection) can influence at a large scale the population behavior and thus the spread of the disease. As detailed in [START_REF] Manfredi | Modeling the interplay between human behavior and the spread of infectious diseases[END_REF], a well-known example is Measles-Mumps-Rubella (MMR) vaccination coverage decrease between 1998 and 2003. This caused a decline in herd immunity and a measles resurgence after the publication of an hypothesis linking MMR vaccination and autism. Behavioral responses to the threat of a pandemic event can influence the spread of the disease in several ways: people can choose preventive protection (by vaccination as presented below) or limit the disease spread by isolation and quarantine (e.g., school closures). Since then, many epidemiological models have been developed to include individual response (see [START_REF] Verelst | Behavioural change models for infectious disease transmission: a systematic review (2010-2015)[END_REF] for a review); on the contrary see [START_REF] Abakuks | Optimal immunisation policies for epidemics[END_REF][START_REF] Behncke | Optimal control of deterministic epidemics[END_REF][START_REF] Morton | On the optimal control of a deterministic epidemic[END_REF][START_REF] Piunovskiy | An explicit optimal intervention policy for a deterministic epidemic model[END_REF] for the treatment of compulsory (as opposed to voluntary) vaccination.

We use the family of compartmental models in order to describe the spread of the disease in the population. This supposes that all individuals are the same and the population is divided into several exclusive classes; the transmission between two classes is modeled through a time-dependent rate. As each compartment represents the status of an individual, it means that the structure of model with immunity is not the same as that of a disease without immunity. SIR terminology usually represents a disease with immunity against re-infection, indicating that the passage of an individual is from the susceptible class S to the infectious class I and then to the recovered class R. If the re-infection is possible, then the terminology used can be SIS, which means that an individual in the infected class I goes back to the susceptible class S. If temporary immunity is given by the disease, the compartmental structure can be SIRS, meaning that after infection, individual stays in R class before going back to the susceptible class. Other well-known compartmental models can include an exposed class (E) that represents the period between being infected and becoming contagious, like in the SEIR or SEIS models. The evolution of the size of the compartments is given by ordinary differential equations (for more details see for instance [START_REF] Hethcote | The mathematics of infectious diseases[END_REF][START_REF] Brauer | Mathematical Epidemiology[END_REF][START_REF] Diekmann | Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation[END_REF]).

In order to model individual vaccination decision in the spread of the disease, vaccination coverage can depend on available information (see [START_REF] Buonomo | Global stability of an SIR epidemic model with information dependent vaccination[END_REF][START_REF] Onofrio | Vaccinating behaviour, information, and the dynamics of SIR vaccine preventable diseases[END_REF][START_REF] Onofrio | Bifurcation thresholds in an sir model with information-dependent vaccination[END_REF]) or be determined by a vaccination game. The evolution of the vaccination coverage could be determined by an imitation dynamics depending (or not) of the available information (see [START_REF] Onofrio | The impact of vaccine side effects on the natural history of immunization programmes: An imitation-game approach[END_REF][START_REF] Onofrio | The interplay of public intervention and private choices in determining the outcome of vaccination programmes[END_REF][START_REF] Bauch | Imitation dynamics predict vaccinating behaviour[END_REF][START_REF] Oraby | The influence of social norms on the dynamics of vaccinating behaviour for paediatric infectious diseases[END_REF][START_REF] Junyuan | Imitation dynamics of vaccine decision-making behaviours based on the game theory[END_REF][START_REF] Onofrio | The impact of vaccine side effects on the natural history of immunization programmes: An imitation-game approach[END_REF]), or be the solution of a Nash equilibrium determined at the beginning of the game, as described below. The Nash-MFG equilibrium corresponds in this case to a strategy where nobody has interest to change his own vaccination strategy. Therefore the existence and the determination of an equilibrium is central.

Before the development of the Nash-MFG theory, some work have been proposed to couple game and epidemiology. Three well-known methods that introduced specific models are presented in the next paragraph with their advantages and limitations. The compartmental model used in [START_REF] Bauch | Group interest versus self-interest in smallpox vaccination policy[END_REF] includes five classes (Susceptible, Exposed, Infected, Recovered, Vaccinated), the dynamics is presented in equation [START_REF] Abakuks | Optimal immunisation policies for epidemics[END_REF] where the population is divided into the following classes: S (susceptible), I (infected), R (recovered) and V (vaccinated). In this case preventive vaccination for the smallpox is available: the model is used to determine the probability of an individual to choose preventive vaccination. To do so, explicit individual payoffs for the two pure strategies (vaccination before or after the outbreak) are defined. Then the existence and uniqueness of the Nash equilibrium is obtained using properties of such payoffs.

dS dt = -βSI -f (S, t), dE dt = βSI -ρE, dI dt = ρE -γI, dR dt = γI, dV dt = f (t, S), (1)
Another method introduced in [START_REF] Bauch | Group interest versus self-interest in smallpox vaccination policy[END_REF] uses the same methodology as previously but on the stationary state of the system (2) (or equivalently system (2bis) that distinguishes immunity given by the disease and immunity given by the vaccine) because of the possible non extinction of the disease. This allows to determine the probability of vaccination when disease is already present in the population with a compartmental system where vital dynamics is considered. As in previous case, the probability p to vaccinate is not time dependent and properties of payoff are used to prove existence and uniqueness of the equilibrium. Finally, solutions for other incidence rate (term βI is replaced by a function λ(I)) are given in [START_REF] Reluga | A general approach for population games with application to vaccination[END_REF] by using Markov Decision Process. Analytical solutions for the vaccination probability at the Nash equilibrium can be found with these approaches but it supposes constant probability of vaccination.

dS dt = µ(1 -p) -βSI -µS dI dt = βSI -γI -µI dR dt = µp + γI -µR (2)
dS dt = µ(1 -p) -βSI -µS dI dt = βSI -γI -µI dR dt = γI -µR dV dt = µp -µV (2bis)
The third and last model being detailed here describes numerical approach when two SEIR models (see for example [START_REF] Abakuks | Optimal immunisation policies for epidemics[END_REF]) are used in order to separate two age classes (see [START_REF] Galvani | Long-standing influenza vaccination policy is in accord with individual self-interest but not with the utilitarian optimum[END_REF]) to obtain the vaccination probability of each age class. The evolution of the probability of vaccination can be modeled by introducing a rule of thumb (see [START_REF] Bauch | Imitation dynamics predict vaccinating behaviour[END_REF][START_REF] Bauch | Evolutionary game theory and social learning can determine how vaccine scares unfold[END_REF]) to study the epidemiological stationary state and the role of parameters like vaccination and infection costs. But reaching the stationary state of the epidemiological system can take a long time as explained in [START_REF] Hubert | Nash-MFG equilibrium in a SIR model with time dependent newborn vaccination[END_REF]. To avoid this problem, [START_REF] Doutor | Optimal vaccination strategies and rational behaviour in seasonal epidemics[END_REF] considers a SIRS model where the transmission rate β is a periodic function, causing the same property for the incidence. Economic approach (see [START_REF] Francis | Dynamic epidemiology and the market for vaccinations[END_REF]) based on the epidemiological model (3) (where notations are the same as above and V(t) is the instantaneous vaccination rate) allows to determine the vaccination and no vaccination region (depending on the state of the disease) by equalization of the two payoffs but no proofs of optimality (existence or uniqueness) are provided.

dS(t) dt = -βS(t)I(t) -V(t) dI(t) dt = βS(t)I(t) -γI(t) dR(t) dt = γI(t)dt dV (t) dt = V(t) (3)
Within the MFG framework, article [START_REF] Laguzet | Individual Vaccination as Nash Equilibrium in a SIR model with application to the 2009-2010 influenza A (H1N1) epidemic in France[END_REF] proposes analytic solutions in some particular cases. Furthermore, finding an equilibrium becomes more complicated if the model has specific behaviour with network interaction or epidemic spread for example ( [START_REF] Fu | Imitation dynamics of vaccination behaviour on social networks[END_REF][START_REF] Colizza | Modeling the worldwide spread of pandemic influenza: Baseline case and containment interventions[END_REF]). Thus finding the equilibrium can quickly become time consuming. One commonly used method is based on finding a fixed point of some function, which can be formulated as the attractor of an evolution system.

The existence of an equilibrium in the general case of MFG is of utmost interest (see for instance [START_REF] Blanchet | Optimal transport and Cournot-Nash equilibria[END_REF][START_REF] Blanchet | Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case[END_REF] for an entry point to this literature). Several approaches were developed, inspired by the general framework of gradient flows (see [START_REF] Jordan | The Variational Formulation of the Fokker-Planck Equation[END_REF]) where the procedure is iterative.

Mathematical framework

Let us recall that in a gradient flow framework (which is not necessarily related to a MFG model !) the equation to be solved is written symbolically (see [START_REF] Ambrosio | Gradient Flows In Metric Spaces and in the Space of Probability Measures[END_REF] for an introduction):

d dτ ξ + ∇F (ξ) = 0 (4)
where the function F is given but ∇F cannot be computed. The numerical counterpart of ( 4) is the famous JKO numerical scheme (introduced by Jordan, Kinderleher and Otto in [START_REF] Jordan | The Variational Formulation of the Fokker-Planck Equation[END_REF]) which, for ξ k an approximation of ξ τ k , can be written as:

ξ k+1 ∈ argmin η d(ξ k , η) 2 2∆ τ + F (η) . (5)
This scheme can be seen as the variational implicit Euler scheme and has been used for the heat equation as a gradient flow of the entropy in the Wasserstein space. Second order in time methods have been proposed very recently (see [START_REF] Legendre | Second-order in time schemes for gradient flows in Wasserstein and geodesic metric spaces[END_REF]).

On the other hand, a MFG equilibrium has no associated function 'F ' and it is likely that no such function exists. For each τ , ξ τ is a probability law. Rather, the search for an equilibrium takes the form of a fixed point:

ξ = J(ξ). (6)
In practice, inspired by the Fictitious Game (see [START_REF] Cardaliaguet | Learning in Mean Field Games: the Fictitious Play[END_REF]) and Best Reply (see [START_REF] Blanchet | Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case[END_REF][START_REF] Blanchet | From Nash to Cournot-Nash equilibria via the Monge-Kantorovich problem[END_REF]) procedures, a new algorithm has been introduced in [44, Eq 3.2] and [START_REF] Turinici | Metric gradient flows with state dependent functionals: the Nash-MFG equilibrium flows and their numerical schemes[END_REF] with the following form:

ξ k+1 ∈ argmin η d(η, ξ k ) 2 2∆τ + P(η, ξ k ) . (7)
In the limit τ → 0 the equation ( 7) describes a curve which under some assumptions on P will be called a solution of the evolution equation:

d dτ ξ τ + ∇ 1 P(ξ τ , ξ τ ) = 0. ( 8
)
On the other hand, the MFG framework provides a natural function P(•, •) even when ξ does not belong to a Hilbert space but only to a metric space. Thus we are led to consider evolution equations in metric spaces with semi-explicit numerical schemes (see [START_REF] Turinici | Metric gradient flows with state dependent functionals: the Nash-MFG equilibrium flows and their numerical schemes[END_REF][START_REF] Blanchet | Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case[END_REF]).

L. LAGUZET

Scope and structure of the paper

The purpose of this work is to introduce variational Runge-Kutta explicit methods of high order in a metric space using a generalization of linearity (presented in [START_REF] Cheeger | Differentiability of Lipschitz functions on metric measure spaces[END_REF]). These approaches are applied to the control of epidemic spread with voluntary vaccination as in [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF].

The paper is structured as follows: Section 2 is dedicated to the Runge-Kutta method recalling its general form (subsection 2.1) and then the introduction of variational Runge-Kutta methods (subsection 2.2). The equivalence of the two approaches in a Hilbert space is proved in subsection 2.3. Then Section 3 gives the epidemiological model used in the numerical application given in Section 4. An economical application is proposed in Section 5. Some considerations are discussed in Section 6.

2 Runge-Kutta methods

In a Hilbert space

Let H be a Hilbert space, f : [0, T ]×H → H a regular function and y : [0, T ] → H satisfying (in some sense to be specified):

d dτ y(τ ) = f (τ, y(τ )) (9)
The Runge-Kutta approach is an iterative method using temporal discretization in order to obtain a numerical solution of an ordinary differential equation of type [START_REF] Blanchet | Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case[END_REF]. The time horizon T is supposed to be finite and can be discretized in (N T + 1) time instants τ 0 = 0, τ 1 = ∆τ, τ 2 = 2∆τ, ..., τ N T = T where ∆τ is the time step. Let y k ∈ H be an approximation of y(τ k ).

Runge-Kutta method is based on the evaluation at intermediate points in time. Equation [START_REF] Brauer | Mathematical Epidemiology[END_REF] presents the general form of the Runge-Kutta method for (9).

y k+1 = y k + ∆τ s i=1 b i p i , p i = f τ n + c i τ, y k + ∆τ i-1 j=1 a i,j p j . (10)
Figure 1: Runge-Kutta method. Left: Butcher tableau. Right: Equation.

A method is defined by the values of the coefficients b i , c i and a i,j , often presented in the Butcher tableau (see table in figure 1). Consistency of such a method is ensured if and only if s i=1 b i = 1. We will present below some particular schemes that will later be formulated in a metric space.

Explicit Euler scheme

The Explicit Euler (denoting EE) scheme is a Runge-Kutta method with only one step. It is defined in figure 2 by the Butcher tableau or equivalently by equation [START_REF] Buonomo | Global stability of an SIR epidemic model with information dependent vaccination[END_REF]. In this case, this scheme is of order one in time. 0 0 1

y k+1 = y k + ∆τ f (τ k , y k ). ( 11
)
Figure 2: Explicit Euler method. Left: Butcher tableau. Right: Equation.

Heun scheme

Heun scheme uses an approximate value p 1 , in order to compute y k+1 , see the figure 3. This scheme is of order two in time.

0 0 0 1 1 0 1/2 1/2 p 1 = y k + ∆τ f (τ k , y k ),
y k+1 = y k + ∆τ 2 f (τ k , y k ) + f (τ k+1 , p 1 ) . (12)
Figure 3: Heun method. Left: Butcher tableau. Right: Equation.

Explicit Runge-Kutta 3 scheme

In the same way, we introduce the Runge-Kutta 3 (denoted RK3) method. It calculates two intermediate values (p 2 , p 3 ) to obtain y k+1 . See figure 4 for the value of coefficients in the Butcher table and equation ( 13) for the explicit equations. In this case, the method is of order three in time.

L. LAGUZET

0 0 0 0 1/2 1/2 0 0 1 -1 2 0 1/6 2/3 1/6 y k+1 = y k + ∆τ 6 p 1 + 4p 2 + p 3 , p 1 = f (τ k , y k ), p 2 = f (τ k+1/2 , y k + ∆τ 2 p 1 ), p 3 = f (τ k+1 , y k + ∆τ (-p 1 + 2p 2 )). ( 13
)
Figure 4: Runge-Kutta 3 method. Left: Butcher tableau. Right: Equation.

Explicit Runge-Kutta 4 scheme

The last method is the celebrated Runge-Kutta 4 (denoted RK4) where three intermediate values (p 2 , p 3 , p 4 ) are computed to determine y k+1 . The value of coefficients is presented in figure 5 and its application to equation ( 9) is given in equation ( 14). This method is of order four in time.

0 0 0 0 0 1/2 1/2 0 0 0 1/2 0 1/2 0 0 1 0 0 1 0 1/6 1/3 1/3 1/6 y k+1 = y k + ∆τ 6 p 1 + 2p 2 + 2p 3 + p 4 , p 1 = f (τ k , y k ), p 2 = f (τ k+1/2 , y k + ∆τ 2 p 1 ), p 3 = f (τ k+1/2 , y k + ∆τ 2 p 2 ), p 4 = f (τ k+1 , y k + ∆τ p 3 ). ( 14
)
Figure 5: Runge-Kutta 4 method. Left: Butcher tableau. Right: Equation.

Variational approach

Unfortunately, none of the schemes presented in Section 2.1 can be used for finding a MFG equilibrium. Firstly because the problem is not really presented that way, and secondly because the space of the unknowns is, in general, not a vector space but a metric space. Therefore we need first to generalize the Runge-Kutta schemes to metric spaces and then to apply them to the MFG setting.

Let (A, d) be a metric space. The definition 2.1 below gives the properties of a geodesic space (see [START_REF] Legendre | Second-order in time schemes for gradient flows in Wasserstein and geodesic metric spaces[END_REF] for additional details).

Definition 2.1. A curve γ : [0, 1] → A is called a (constant speed) geodesic provided that d(γ(s), γ(τ )) = |s -τ | • d(γ(0), γ(1)
) and the space is called geodesic if, for any couple of points (X, Y ) ∈ A 2 , there exists at least a geodesic γ connecting them, that is, γ(0) = X, γ(1) = Y .

From now on we assume that (A, d) is a geodesic space. Consider P : A × A → R a function with the following properties:

(H1) ∀n ≤ 4, ∀r > 0, ∀ξ, Y 1 , . . . , Y n ∈ A, the set of vectors {(P(z, Y 1 ), . . . , P(z, Y n )); z ∈ A, d(z, ξ) ≤ r} is compact as a subset of R n , (H2) for any point Y ∈ A and any X 0 , X 1 ∈ A a constant speed geodesic γ exits such that γ(0) = X 0 , γ(1) = X 1 and the function τ → P(γ(τ ), Y ) from [0, 1] to R is linear.
Examples:

• If H is a Hilbert space, f : H → H a function then (x, y) → x, f (y)
satisfies hypothesis (H1) and (H2), see proof of Proposition 2.2. • Let A be the set of probability measures on [0, T ]; F be the set of Lipschitz functions with constant less than 1 and the metric of weak convergence d is defined, for a 1 , a 2 ∈ A by (see [START_REF] Shiryaev | Probability-1[END_REF]): [START_REF] Shiryaev | Probability-1[END_REF]) and the application z

d(a 1 , a 2 ) = sup f ∈F | T 0 f (t)a 1 (dt) - T 0 f (t)a 2 (dt)|. For Y ∈ A, t ∈ [0, T ] we define g Y (t) = sin(t) T 0 cos(s)Y (ds) 2 and for z, Y ∈ A, Q(z, Y ) = T 0 g Y (t)z(dt). Let (z k ) k∈N be a sequence of elements in A converging weakly to z ∈ A then, for f ∈ F, T 0 f (t)z k (dt) converges to T 0 f (t)z(dt). As (A, d) is compact (see
∈ A → (Q(z, Y 1 ), . . . , Q(z, Y n )) ∈ R n continuous, Q satisfies (H1).
Consider a 1 and a 2 ∈ A and for s ∈ [0, 1],

γ s = sa 1 + (1 -s)a 2 . Then d(γ s , γ τ ) = sup f ∈F | T 0 f (u)γ s (du) - T 0 f (u)γ τ (du)| = sup f ∈F | T 0 f (u)(τ - s)a 1 (du)- T 0 f (u)(τ -s)a 2 (du)| = |τ -s| sup f ∈F | T 0 f (u)a 1 (du)- T 0 f (u)a 2 (du)| = |τ -s| • d(γ(0), γ(1)). It shows that γ is a (constant speed) geodesic. Then, the function s → Q(γ s , Y ) of the form s → s T 0 g Y (u)a 1 (du) + (1 - s) T 0 g Y (u)a 1 (du) is linear.
• Note that previous example also works with the total variation distance.

The evolution equation that we want to solve has the form of equation ( 8) (but note that this is only a formal expression as d dτ ξ τ does not have a well-defined meaning in a general metric space, see [START_REF] Turinici | Metric gradient flows with state dependent functionals: the Nash-MFG equilibrium flows and their numerical schemes[END_REF]).

This section develops a variational Runge-Kutta method in a metric space. We propose a method that does not assume any vector calculus but uses the hypothesis (H2) which defines the notion of linear application in a metric L. LAGUZET space (see also [13, Section 8, Definition 8.1 page 480] for more details). As ∆τ tends to zero equation ( 7) describes a curve. The purpose of the scheme is to describe this evolution faster and more precisely that is to say to obtain high order schemes.

Variational Explicit Euler scheme

The first scheme of that form was introduced in [44, Eq 3.2]. Formally, we define the variational Explicit Euler scheme (VEE ) by:

ξ k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + P(η, ξ k ) . (15)

Variational Heun scheme

The variational Heun scheme (VH ) introduces, as in vector space, an intermediate value ξ k+1 (equation ( 16a)) and uses it to compute ξ k+1 (equation (16b)).

ξ k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + P(η, ξ k ) , (16a)
ξ k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 2 P(η, ξ k ) + 1 2 P(η, ξ k+1 ) . (16b)
Note that two minimizations are required in order to obtain ξ k+1 .

Variational Runge-Kutta 3 scheme

The variational Runge-Kutta 3 scheme (VRK3 ) is defined by equation [START_REF] Onofrio | The interplay of public intervention and private choices in determining the outcome of vaccination programmes[END_REF]. Two intermediate values are computed in this case ( ξ 1 k+1 with equation (17a) and ξ 2 k+1 with equation (17b)) in order to obtain ξ k+1 using equation (17c).

ξ 1 k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 2 P(η, ξ k ) , (17a)
ξ 2 k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ -P(η, ξ 1 k+1 ) + 2P(η, ξ 2 k+1 ) , (17b)
ξ k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 6 P(η, ξ k ) + 2 3 P(η, ξ 1 k+1 ) (17c) + 1 3 P(η, ξ 2 k+1 ) .

Variational Runge-Kutta 4 scheme

The variational Runge-Kutta 4 scheme is presented in equation [START_REF] Onofrio | Vaccinating behaviour, information, and the dynamics of SIR vaccine preventable diseases[END_REF]. Three intermediate values ( ξ 1 k+1 with equation (18a), ξ 2 k+1 with equation (18b) and ξ 3 k+1 defined by equation (18c)) are computed in order to obtain ξ k+1 by using equation (18d).

ξ 1 k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 2 P(η, ξ k ) , (18a)
ξ 2 k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 2 P(η, ξ 1 k+1 ) , (18b)
ξ 3 k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 2 P(η, ξ 2 k+1 ) , (18c)
ξ k+1 ∈ argmin η∈A d(η, ξ k ) 2 2∆τ + 1 6 P(η, ξ k ) + 1 3 P(η, ξ 1 k+1 ) (18d) + 1 3 P(η, ξ 2 k+1 ) + 1 6 P(η, ξ 3 k+1 ) .
Note that four minimizations are needed in order to obtain ξ k+1 .

See section 4.2 for a numerical application and the study of the order in time of these schemes.

Property of variational scheme

In this subsection we prove that using hypotheses (H1) and (H2), the variational schemes introduced above are well-defined (Prop. 2.1) and correspond, when the ambient space is Hilbert, to the Runge-Kutta method presented in section 2.1 (Prop. 2.2).

Prop 2.1. Under the hypothesis (H1) and (H2), the equations (15), ( 16), ( 17) and (18) defining respectively the VEE, VH, VRK3 and VRK4 schemes admit a solution, i.e., the schemes are well defined.

Proof. As the argumentation is very similar for all schemes, we present the general case where ξ k+1 ∈ argmin η∈A {F (η)} with

F : A → R, for fixed ξ ∈ A, F (η) = d(η,ξ) 2
2∆τ +l(η) where l satisfies hypothesis (H2) and the following hypothesis (H1bis):

(H1bis) ∀r > 0, ∀ξ ∈ A, the set {l(z); z ∈ A, d(z, ξ) ≤ r} is compact as a subset of R.

We will prove that the application F attains its minimum in A.

Let (z k ) k∈N ∈ A be a sequence such that lim k→∞ F (z k ) = m with m = inf η∈A F (η). We show that d(z k , ξ) is bounded then we find a minimizer.

Suppose d(z k , ξ) = λ > 1 and let γ be the geodesic connecting ξ and z k (γ(0) = ξ and γ(1) = z k ) satisfying (H2). Consider an element z = γ(t) on the geodesic γ such that d(z, ξ) = 1. By definition of a geodesic,

d(z, ξ) = d(γ(t), γ(0)) = |t -0| • d(γ(0), γ(1)) = td(z k , ξ) = λ
, implies necessarily t = 1/λ and using hypothesis (H2), we obtain:

l(z) = l(γ(t)) = tl(γ(1)) + (1 -t)l(γ(0)), = (1 - 1 λ )l(ξ) + 1 λ l(z k ). ( 19
)
For

k large enough, z k is such that F (z k ) < 2m that is λ 2 2∆τ + l(z k ) < 2m.
Rearranging equation [START_REF] Onofrio | Bifurcation thresholds in an sir model with information-dependent vaccination[END_REF] gives l(z k ) = λl(z) -(λ -1)l(ξ), note that l(ξ) is constant and as d(z, ξ) = 1, hypothesis (H1bis) insures that l(z) is bounded, so l(z k ) is also bounded (by M ∈ R). Consequently, there exists a constant C, depending on m, τ and M such that λ < C.

As d(z k , ξ) is bounded, there exists a converging subsequence, (

x k ) k∈N ∈ A such that d(x k , ξ) → d 1 , and l(x k ) converges to m - d 2 1
2∆τ . Starting from some index k, x k ∈ B(ξ, d 1 ) (the ball with center ξ and radius d 1 ) and if it is not the case, we replace x k by the element at the intersection of the ball and the geodesic γ linking x k and ξ. This operation does not change the limit of l(z k ): noting temporarily xk the obtained sequence, as l satisfies hypothesis (H2), we have l(

x k ) = (1-d 1 d(x k ,ξ) )l(ξ)+ d 1 d(x k ,ξ) l(x k ) with d(x k , ξ) → d 1 .
By hypothesis (H1bis), the set l(x k ) is closed, so there exists an element Z ∈ B(ξ, d 1 ) such that l(Z) = m -

d 2 1 2∆τ with, of course, d(Z, ξ) ≤ d 1 . In this case, F (Z) = d(Z,ξ) 2 2∆τ + m - d 2 1 2∆τ ≤ m, so Z is a minimizer. Let n ≤ 4, Y 1 , . . . , Y n ∈ A, λ 1 , .
. . , λ n ∈ R and P satisfying (H1) and (H2). We define the linear combination l of P by: l(η) = n i=1 λ i P(η, Y i ). It is immediate that l satisfies (H2); we prove below that if P satisfies (H1) then l satisfies (H1bis). As all elements of the sum are in a compact set, the sum is bounded. Let (z k ) k∈N ∈ A be a sequence such that lim k→∞ l(z k ) = α. Using (H1), there exists a subsequence of z k (also noted z k ) such that the element (P(z k , Y 1 ), . . . , Proof. We prove hypothesis (H1) in the case n = 1, the proof for the general case is left to the reader. We show that, for ξ ∈ H, D = { z, -f (ξ) ; z ∈ H, d(z, ξ) ≤ r} is bounded and closed. Denoting ||•|| the norm of H, by the triangular inequality when x ∈ B(ξ, r), ||x|| ≤ r +||ξ||. For

P(z k , Y n )) ∈ R n converges to (α 1 , . . . α n ) so l(z k ) converges to l(z) with l(z) = n i=1 λ i α i = α. Prop 2.
x ∈ D, || x, f (ξ) || ≤ ||x||•||f (ξ)|| ≤ (r +||ξ||)•||f (ξ)|| < ∞, so D is bounded. Let z be an element of B(ξ, r) and x ∈ H such that z = ξ + x with ||x|| ≤ r, by linearity z, -f (ξ) = ξ, -f (ξ) + x, -f (ξ) . Furthermore, x, -f (ξ) ≤ ||x|| • ||f (ξ)|| ≤ r||f (ξ)|| so { z, -f (ξ) , z ∈ B(ξ, r)} is included in the closed interval [ ξ, -f (ξ) - r||f (ξ)||, ξ, -f (ξ) + r||f (ξ)||]. If we take x = -λ f (ξ) ||f (ξ)|| with λ ∈ [-r, r] then x, -f (ξ) = λ||f (ξ)|| which implies z, -f (ξ) = ξ, -f (ξ) +λ||f (ξ)||. There- fore { z, f (ξ) , z ∈ B(ξ, r)} is exactly [ ξ, f (ξ) -r||f (ξ)||, ξ, f (ξ) + r||f (ξ)||]
and is closed as a closed interval of R.

To prove hypothesis (H2), consider x 0 , x 1 ∈ H and for t ∈ [0, 1], γ(t) = tx 1 + (1 -t)x 0 the segment between x 0 and x 1 . For t, s ∈ [0, 1], with straightforward computation, we have d(γ(t), γ(s)) = |t -s| • d(γ(0), γ(1)), showing that γ is a (constant speed) geodesic. In this case, we can also prove the uniqueness of the geodesic. In fact, assume that there exists another geodesic γ 2 linking x 0 and x 1 . Let z = γ 2 (t) be an element only on γ 2 ; using the triangle inequality and that γ is the segment linking x 0 and x 1 we have d(x 0 , z) + d(z, x 1 ) > d(x 0 , x 1 ). Furthermore as γ 2 is a geodesic, left member of previous equation is td(x 0 , x 1 )+(1-t)d(x 0 , x 1 ) = d(x 0 , x 1 ). We obtain a contradiction thus there does not exist another geodesic. The function t → P(γ(t), Y ) has the form t → tx 1 + (1 -t)x 0 , -f (Y ) , so is linear.

We now prove the equivalence of schemes. The scheme VEE is defined by the following equation:

ξ k+1 ∈ argmin η∈H ||η -ξ k || 2 2∆τ + η, f (ξ k ) .
The application

F : H → R, F (η) = ||η-ξ k || 2 2∆τ + η, f (ξ k ) is differentiable and its derivative F : H → H is F (η) = η-ξ k ∆τ + f (ξ k ).
The minimizer of F is a critical point, the minimum noted ξ k+1 , satisfies F (ξ k+1 ) = 0. After trivial rearrangement, we obtain ξ k+1 = ξ k -∆τ f (ξ k ), that is the Explicit Euler scheme presented in [START_REF] Buonomo | Global stability of an SIR epidemic model with information dependent vaccination[END_REF] for equation [START_REF] Blanchet | Remarks on existence and uniqueness of Cournot-Nash equilibria in the non-potential case[END_REF]. Arguments are similar for other schemes.

Epidemiological model

From now on, all variables can have a double dependence: on the real time t which will be subscript and on the iteration time τ (previously k in equation L. LAGUZET [START_REF] Blanchet | Optimal transport and Cournot-Nash equilibria[END_REF]) which will be a superscript l.

This section reproduces a short presentation of the epidemiological model and of the cost computing as introduced in [44, Sec 2]. Please refer to the original article for more details and proof.

Spread of the disease

The model simulates the dynamics of an epidemic in a population. The final time period T is supposed finite, and the time horizon can be discretized in (N + 1) time instants, noted t 0 = 0, t 1 = ∆T, ..., t N = T .

The model is compartmental, and the population is divided into several classes: susceptible individuals (S n ) is the proportion of individuals susceptible to catch the disease at time instant t n ; infected individuals regroups individuals who are infected, more precisely, I ω n is the proportion of individuals who have been infected at time instant t n-ω (with ω ∈ {0, 1, ..., Ω}) and I n is the sum of all I ω n ; recovered individuals (R n ) is the proportion of individuals who have exited from one infected class; vaccinated individuals (V θ n where θ counts the time between vaccination time instant and current time instant with θ ∈ {0, 1, ..., Θ}) is the proportion of vaccinated individuals at time instant t n and not infected since; failed vaccinated individuals (F n ) is the proportion of individuals who were vaccinated at some time t ≤ t n but whose vaccination failed and have not been infected since. As individuals can only take vaccination once, V Ω represents individuals who have lost immunity. Lastly U n represents the proportion of individuals vaccinated between t n and t n+1 , a proportion f of them will never develop any immunity and go to the failed vaccination class.

The specificities of disease are represented by the vector γ = (γ 0 , ..., γ Ω ) ∈ (R + ) Ω+1 (which reflects the recovery of an infected individual); the function β(t) (which characterizes the contact between an infected individual and a susceptible one at time t); and the vector α ∈ [0, 1] which quantifies the protection provided by the vaccine in terms of the probability of infection.

The model is defined by the following system of equations:

S n+1 = (1 -β n ∆T I n ) (S n -U n ) ( 20a
)
I 0 n+1 = β n ∆T F n + S n + N -1 θ=0 α θ V θ n I n ( 20b
)
I ω+1 n+1 = (1 -γ ω ∆T )I ω n ω = 0, . . . , Ω -1 (20c) V 0 n+1 = (1 -f ) • (1 -β n ∆T I n ) U n (20d) V θ+1 n+1 = (1 -β n ∆T α θ I n ) V θ n , θ = 0, . . . , Θ -2 (20e) V Θ n+1 = (1 -β n ∆T α Θ-1 I n ) V Θ-1 n + (1 -β n ∆T I n ) V Θ n (20f) F n+1 = f • (1 -β n ∆T I n ) U n + F n (1 -β n ∆T I n ) (20g)
with initial conditions:

S 0 = S 0 -, I ω 0 = I ω 0 -, V θ 0 = 0, ∀θ ≥ 0. ( 20h
)
The continuous version of the previous discrete model is given by the following equations:

S(t) = -βS(t)I(t) -U (t) (21a)
I(t) = β S(t) + F (t) + +∞ 0 A(θ)V (t, θ)dθ I(t) -γI(t) (21b) ∂ t V (t, θ) + ∂ θ V (t, θ) = -βA(θ)V (t, θ)I(t) (21c)
F (t) = f U (t) -βF (t)I(t) (21d)
with initial and boundary conditions:

S(0) = S 0 -, I(0) = I 0 -, ∀t ≥ 0 : S(t) ≥ 0, F (0 -) = 0, ∀θ ≥ 0 : V (0 -, θ) = 0, V (t, 0) = (1 -f )U (t). ( 21e
)
In this continuous time model, we use the following notation (refer to [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF]Appendix] for details):

• U : [0, ∞] → [0, 1] represents the vaccination, U (t)
is the fraction of the population that has been vaccinated by time t;

• constants β, γ, f and function A, have the same meaning as in the discrete model;

• V (t, θ) is the fraction of individuals who have been vaccinated at time t -θ and have not been infected since vaccination.

L. LAGUZET

Individual vaccination: as minimization of individual cost

The system [START_REF] Doutor | Optimal vaccination strategies and rational behaviour in seasonal epidemics[END_REF] presents the spread of the disease at the population level. For an individual, we consider the Markov chain M n which describes, in terms of transition probabilities (see [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF]Eq 2.10]), the state of the individual (M n ∈ {S, R, F, I j , V l } with j ∈ {0, . . . , Ω} and l ∈ {0, . . . , Θ}).

In this subsection, we look for the global decision of vaccination U n which corresponds to the sum of individual decisions, as shown by equation [START_REF] Galvani | Long-standing influenza vaccination policy is in accord with individual self-interest but not with the utilitarian optimum[END_REF]. In order to reflect the decision of an individual, we introduce the probability law µ defined on {t 0 , ..., t N -1 } ∪ ∞. In practice, before the epidemic starts, the individual chooses his probability of vaccination at each time step t k , assuming his non-infected status.

The collection of conditional rates λ = (λ n ) N n=1 is given by the law µ:

∀n ≤ N -1 : λ n =    µ n µ n + • • • + µ ∞ , if µ n + • • • + µ ∞ > 0 0, otherwise (22)
The individual in the susceptible class at time t 0 (M 0 = S) has the following cost: J indi (µ; U ) = µ, C U where C U (the vector representing the cost of all pure strategies "vaccination sure at time t") is defined by: 20) is attained when:

C U (t n ) = r I ϕ I n + (1 -ϕ I n )(r V + (1 -f )r I ϕ V,I n ) + r I f (ϕ I ∞ -ϕ I n ) for n < N, r I ϕ I ∞ for n = N. ( 23
) Here ϕ V,I n = 1-Θ k=n 1 -β k ∆T α k-n-1 I k , with the convention α -1 = 1 and ϕ I n = 1 -n k=0 1 -β k ∆T I k (for n = 0, ..., N -1), ϕ I ∞ = 1 -N k=0 1 -β k ∆T I k .

The equilibrium between individual dynamics and global dynamics (

U n = λ n S n . ( 24
)
The purpose of an individual is to minimize

J indi (µ; U ) under the constraint µ ∈ Σ N +1 with Σ N +1 = x ∈ R N +1 N k=0 x k = 1 and x k ≥ 0, 0 ≤ k ≤ N .
The numerical schemes give a sequence of probability laws µ l ; here the index "l" means that µ l is a numerical approximation of the solution ξ of (8) at the time τ l = l • ∆τ . For instance, if we consider the discrete model [START_REF] Doutor | Optimal vaccination strategies and rational behaviour in seasonal epidemics[END_REF], each µ l = (µ l 0 , µ l 1 , . . . , µ l N -1 , µ l ∞ ) ∈ Σ N +1 is a probability law on {t 0 ; . . . , t N -1 } ∪ {∞}, while if we consider the continuous time model [START_REF] Francis | Dynamic epidemiology and the market for vaccinations[END_REF], each µ l is a probability law on [0, T ]∪{∞}. Note that here a step of the numerical schemes increments the "l" index in order to converge to a Nash-MFG equilibrium. For the epidemic parameters, we consider a total simulation time at one year (T = 1); an initial proportion of susceptibles S 0 = 0.94 and infected I 0 = 2.0 -6 ; three time instants by day (N = 365 * 3); a recovery rate γ ω = 365/3.2 (Ω = 20); the reproduction number R 0 = 1.35 thus β = γR 0 ; t β 2 = 1/2 such that β(t) = β for t ≤ t β 2 and then β(t) = β min for t > t β 2 where β min = γ/S 0 ; the relative cost of the epidemic is r I = 1.

Numerical values

Finally, we introduce t α 1 , t α 2 and set

α θ = 1 -1 [t α 1 ,t α 2 ]
. For the case Short persistence, large efficacy, t α 1 = 5/365 and t α 2 = 1/12 (Ω = 93). For the case Long persistence, 100% efficacy, t α 1 = 0 and t α 2 = 1/2 (Ω = 549). We suppose a failure rate f = 0, and the relative cost of the vaccination r V = 0.005.

Numerical illustration 4.1 Framework

F. Salvarani and G. Turinici introduced an epidemiological model with possibility of vaccine which has imperfect efficiency and limited persistence recalled in the previous section 3. To model the spread of the disease in the population, they use a compartmental model and a probability distribution to reflect the individual decisions.

To find stable individual decision they define the problem as a Nash equilibrium. They find the probability distribution µ making an individual indifferent to change his vaccination decision if all individuals have the same µ (see [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF]Theorem 2.1] for the proof of equilibrium existence).

Let C µ be the cost of pure strategies "vaccination happens at time t" when all individuals choose as vaccination strategy µ. In that case η, C µ represents the cost of an individual with strategy η when others use strategy µ.

The definition of the mapping E(µ) as introduced in [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF] is the maximum gain obtained by an individual if he changes unilaterally his strategy and everybody else remains with the strategy µ; with mathematical notation:

E(µ) = µ, C µ -min η∈Σ N +1 η, C µ .
A minimum of the mapping µ → E(µ) is a Nash equilibrium. To find it they introduce an iterative method depending on a step ∆τ (algorithm 1). In L. LAGUZET this algorithm, equation ( 25) is used. The following intuitive interpretation is also provided: an individual in a population with strategy µ l will, if necessary, adjust his strategy to minimize his cost function η → η, C µ l while at the same time keeping small the distance between the previous strategy µ l and the new µ l+1 . But algorithm 1 can also be seen as numerical resolution of an evolution e.g., describing a curve in the metric space of the admissible strategies. That is why we apply variational methods in the metric space Σ N +1 with the standard Euclidean distance to obtain faster convergence to the Nash equilibrium. Here P(x, y) = x, C y and for the same reasons as in the proof of Prop. 2.2, P satisfies hypotheses (H1) and (H2). Note that with high order schemes VH, VRK3 and VRK4, we lose the intuitive idea but we increase the Nash equilibrium computing.

Algorithm 1 Finding a Nash equilibrium (as introduced in [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF])

1: Choose a step ∆τ > 0 and a starting distribution µ 0 ∈ Σ N +1 .

2: Set iteration count l = 1.

3: Compute µ l+1 with following formula:

µ l+1 ∈ argmin η∈Σ N +1 d(η, µ l ) 2 2∆τ + η, C µ l . (25)
4: If E(µ l+1 ) is smaller than a given tolerance then stop and exit, otherwise set l → l + 1 and go back to step 3.

Results

In this section we test the variational schemes proposed in Section 2.2. The reader is invited to refer to [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF]Sec 4.6] or Section 3 for the numerical values of the two studied cases Short persistence, large efficacy (with vaccine persistence at one month with a delay of action at five days) and Long persistence, 100% efficacy (with persistence of the vaccine at six months with no delay of action).

In order to appreciate the convergence scheme for the case Short persistence, large efficacy the evolution of the mapping E(•) is presented (see figure 6 for VEE with ∆τ = 0.1), where the evolution for the four schemes with lower ∆τ value is added to stabilize result of mapping E(•).

. In [START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF] VEE with ∆τ = 0.1 is used.

For the Short persistence, large efficacy case, figure 7 presents the numerical estimation of the scheme error and order: only the VEE is of order one, the others are of order two in ∆τ . Furthermore, figure 9 introduces the obtained errors for each schemes depending of the number of minimizations. Recall that a high order scheme needs a very regular functional to provide high order convergence, but in our case the regularity of the function is completely unknown. However, this case remains very interesting because it shows that even if the regularity of the function is not enough to obtain order three or four, the VRK3 and VRK4 still improve the regularity of the function and have a better order than the VEE scheme.

For the other example Long persistence, 100% efficacy, the convergence is faster: figure 11 compares the evolution of mapping E(•) for the two cases. In the second case, the graph in figure 8 shows that the scheme has high order as it was expected previously. VH is indeed of order 2 and VRK3 and VRK4 are respectively order three and four. Furthermore, figure 10 introduces the error for each scheme depending on the number of minimisations: this shows the accuracy of the schemes. Other steps used are ∆τ = 0.04, 0.06, 0.08, 0.10 and 0.12. Denote µ ∆τ n the probability of vaccination at time n/N given by the last iteration (10/∆τ ) of the numerical scheme (V EE, V H, V RK3 or V RK4) with parameter ∆τ . We plot:

(left) the error N n=1 µ ∆τ n -µ ∆τ ref n 2 + µ ∆τ ∞ -µ ∆τ ref ∞ 2
. Right: Order of convergence. This figure shows for each of the four detailed schemes the error depending on number of minimisations. Reference solution is given by a VRK4 scheme with ∆τ ref = 0.001. Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12. 5 Second application: an economical approach

Epidemiological and economic model

This numerical application is based on [START_REF] Francis | Optimal tax/subsidy combinations for the flu season[END_REF]. In order to model the spread of the disease, three classes (as explained in Section 1.1) are introduced: Susceptible (S), Infected (I), Recovered (R) with the following time evolution:

Ṡ(t) = -β S(t)I(t) N -r(t) İ(t) = β S(t)I(t) N -γI(t) Ṙ(t) = r(t) + γI(t) (26)
where S(t), I(t), R(t) represent respectively the number of susceptible, infected and recovered individuals at time t, N the total number of individuals and r(t) the rate of vaccination at time t. We also introduce the utility u H of a healthy individual and the utility u S when sick. The individual cost introduced is represented by the expected loss if individual is falling ill EL

= u H -u S γ 1 -S(∞) S(0)
and the cost of vaccination θ.

The computation that follows, taken from [START_REF] Francis | Optimal tax/subsidy combinations for the flu season[END_REF], is valid when the dynamics (26) traverses at most two regions of vaccination / non vaccination. Under this assumption, the expected payoff of the two pure strategies (vaccination or not) is equal on the frontier of the vaccination / non vaccination region, that is:

θ = u H -u S γ 1 -S(∞)
S(0) . Set w = θγ/(u H -u S ). Rewrite previous equation under the following form: w = 1 -S(∞)/S(0), that gives S(∞) = S(0)(1 -w). Using that, with ρ = N γ/β, d dt S(t)e R(t)ρ = 0 for r(t) = 0 ∀t ≥ 0, S(0)e R(0)ρ = S(∞)e R(∞)ρ (see [START_REF] Abakuks | Optimal immunisation policies for epidemics[END_REF]), and after replacing of S(∞) by S(0)(1 -w), we obtain: R(0) = ρ ln(1 -w) + R(∞). As the total number of the population is constant equal to N , we obtain equation ( 27). This last equation can be seen as a "switching curve" in the plane (S, I):

Γ = (S, I) I = -N γ β ln 1 - θγ u H -u S - θγ u H -u S S , (27)
where the individual is indifferent between being vaccinating and doing nothing.

With this economical approach, there is no insurance that the curve Γ generates a Nash-MFG equilibrium. In fact, [START_REF] Francis | Optimal tax/subsidy combinations for the flu season[END_REF] does not use Nash-MFG framework but assumes equality between two payoffs, false in general cases (see for example [START_REF] Galvani | Long-standing influenza vaccination policy is in accord with individual self-interest but not with the utilitarian optimum[END_REF]).

Nevertheless we seek to replicate these results in order to benchmark the performance of our numerical schemes. For each initial point, the algorithm reports if vaccination is optimal or not, which can be compared to the region given by the switching curve Γ.

We use the same approach as in Subsection 3.2 to determine the global vaccination response r(t). The proportion of individuals vaccinated up to time t is the cumulative function V (t) defined by:

V (t) = t 0 r(u)du, (28)
where r(u) is the instantaneous vaccination rate at time u with r(u) ∈ [0, r max ] and r max < ∞.

Numerical values

For the simulation, we use the following numerical values: an initial proportion of susceptibles S 0 = 0.75 and infected I 0 = 0.1; a cost of the vaccine θ = 0.5; an expected utility loss associated if illness γ u H -u S = 1; a recovery rate for an infected individual γ = 365/10 and a transmission rate of the disease β = 73. As in the Section 4, total simulation time is one year (T = 1) with L = 365 × 3 total time steps.

L. LAGUZET

Results

We obtain coherent results with the analytic solution by using the same approach as in Subsection 4.1. The equilibrium individual strategy is a mixed strategy with ξ 0 = 66% and ξ ∞ = 33% and its cost is 0.51. This results are coherent with the analytic result that provides a mixed strategy with ξ 0 = 34%, ξ ∞ = 66% and a cost at 0.5. Figure 12 presents the numerical estimation errors and order for the four schemes and figure 13 the errors for each schemes depending of the number of minimizations.

10 -1.4 10 -1.3 10 -1.2 10 -1.1 10 -1 10 -0.9 10 -0.8 10 -0. Other steps used are ∆τ = 0.04, 0.06, 0.08, 0.10 and 0.12. Denote µ ∆τ,l n the probability of vaccination at time n/N given by the l -th iteration of the numerical scheme (V EE, V H, V RK3 or V RK4) with parameter ∆τ . We plot: Figure 13: Numerical error obtained (see figure 12 for the definition of the error) for the four schemes depending on the number of minimizations for the model [START_REF] Huang | An invariance principle in large population stochastic dynamic games[END_REF]. Reference solution is given by a VRK4 scheme with ∆τ ref = 0.005. Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12.

Perspectives

This work introduces three high order schemes to find a Nash equilibrium and illustrates their use in an epidemiological application. We show numerically that the schemes VH, VRK3 and VRK4 exhibit better order than VEE and we can obtain, depending on the regularity of the function, up to order four convergence. The approach based on the Runge-Kutta method can be applied for other numerical schemes, for instance Mid Point, Leap-Frog or Adams-Bashford methods. The presented algorithms are not optimized, as minimizations can take a long time and need a high numerical precision, so the time execution can be significant. A perspective of this work is to provide an extension to bi-dimensional or tri-dimensional problems.

c 1 a

 1 1,1 a 1,2 . . . a 1,s c 2 a 2,1 a 2,2 . . . a 1,s . . . . . . . . . . . . . . . b 1 b 2 . . . b s

2 .

 2 Let H be a Hilbert space, f : H → H be a function and P(x, y) = x, -f (y) . Then the schemes VEE (equation (15)), VH (equation (16)), VRK3 (equation (17)), and VRK4 (equation (18)) correspond to Explicit Euler (equation (11)), Heun (equation (12)), Runge-Kutta 3 (equation (13)) and Runge-Kutta 4 (equation (14)) schemes respectively for the evolution equation (9).

  This section describes the numerical values used in the two cases presented in this paper: Short persistence, large efficacy and Long persistence, 100% efficacy. These values are sensibly similar to the ones used in [44, Subsec 4.2 and 4.3].

Figure 6 :

 6 Figure6: The evolution of the mapping E(•) for the four schemes introduced in section 2.2 in the case Short persistence, large efficacy for the model (20a) -(20h). In[START_REF] Salvarani | Optimal individual strategies for influenza vaccines with imperfect efficacy and limited persistence[END_REF] VEE with ∆τ = 0.1 is used.

Figure 7 :

 7 Figure7: Numerical estimated error and order obtained for the four schemes for the case Short persistence, large efficacy of the model (20a) -(20h). Reference solution (µ ∆τ ref ) is given by a VRK4 scheme with ∆τ ref = 0.001. Other steps used are ∆τ = 0.04, 0.06, 0.08, 0.10 and 0.12. Denote µ ∆τ n the probability of vaccination at time n/N given by the last iteration (10/∆τ ) of the numerical scheme (V EE, V H, V RK3 or V RK4) with parameter ∆τ . We plot:

Figure 11 :Figure 8 :

 118 Figure 11: Mapping E(•) for the case Short persistence, large efficacy and Long persistence, 100% efficacy (of the model (20a) -(20h)), generated with VEE scheme and ∆τ = 0.1.

Figure 9 :

 9 Figure 9: Numerical estimated error obtained (see figure 7 for the definition of the error) for the Short persistence large efficacy case of the model (20a) -(20h).This figure shows for each of the four detailed schemes the error depending on number of minimisations. Reference solution is given by a VRK4 scheme with ∆τ ref = 0.001. Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12.

Figure 10 :

 10 Figure 10: Numerical estimated error obtained (see figure8for the definition of the error) for the Long persistence, 100% efficacy case of the model (20a) -(20h). This figure shows for each of the four detailed schemes the error depending on number of minimisations. Reference solution is given by a VRK4 scheme with ∆τ ref = 0.005. Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12.

7

 7

Figure 12 :

 12 Figure12: Numerical error and order for the schemes and the model[START_REF] Huang | An invariance principle in large population stochastic dynamic games[END_REF]. Reference solution (µ ∆τ ref ) is given by a VRK4 scheme with ∆τ ref = 0.005. Other steps used are ∆τ = 0.04, 0.06, 0.08, 0.10 and 0.12. Denote µ ∆τ,l

( 2 + 2 .

 22 left) the error max l≤10/∆τ N n=1 µ ∆τ,l n -µ ∆τ ref ,l n µ ∆τ,l ∞ -µ ∆τ ref ,l∞Right: Order of convergence.

Acknowledgements

This work was supported by the French Research Agency, program Investments for the future, project ANR-10-BINF-07 (MIHMES) and by the European Union, through the European fund for the regional development of Pays-de-la-Loire (FEDER).