Toufik Bakir
email: toufik.bakir@u-bourgogne.fr

Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr

Loïc Bourdin
email: loic.bourdin@unilim.fr

Jérémy Rouot
email: jeremy.rouot@grenoble-inp.org

Pontryagin-Type Conditions for Optimal Muscular Force Response to Functional Electrical Stimulations

Keywords: Functional electrical stimulation, Muscle mechanics, Optimal control problems, Sampled-data controls, Pontryagin-type necessary optimality conditions Mathematics Subject Classification (2000) 49K15, 93B07, 92B05

In biomechanics, recent mathematical models allow one to predict the muscular force response to functional electrical stimulations. The main concern of the present paper is to deal with the computation of optimized electrical pulses trains (for example in view of maximizing the final force response). Using the fact that functional electrical stimulations are modeled as Dirac pulses, our problem is rewritten as an optimal sampled-data control problem, where the control parameters are the pulses amplitudes and the pulses times. We establish the corresponding Pontryagin first-order necessary optimality conditions and we show how they can be used in view of numerical simulations.

Introduction

Predicted muscular force response to functional electrical stimulations (in short, FES) is utilized in biomechanics for muscular reeducation and in case of paralysis. The main concern of the present paper is to deal with the mathematical computation of optimized electrical pulses trains (for example in view of maximizing the final force response). A simplified model (see [START_REF] Law | Mathematical models of human paralyzed muscle after long-term training[END_REF][START_REF] Ding | Two-step, predictive, isometric force model tested on data from human and rat muscles[END_REF][START_REF] Ding | Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains[END_REF]) is a non-fatigue model derived from the Hill equation 1 and used in the context of biochemistry and pharmacology [START_REF] Gesztelyi | The Hill equation and the origin of quantitative pharmacology[END_REF]. More complete models (see [START_REF] Ding | A predictive model of fatigue in human skeletal muscles[END_REF][START_REF] Ding | Mathematical models for fatigue minimization during functional electrical stimulation[END_REF]), taking into account the muscle fatigue, were obtained recently in the framework of model identification and produce a dynamics described by a set of five differential equations. In the present article, we will use the so-called Ding et al. force-fatigue model [START_REF] Ding | Mathematical models for fatigue minimization during functional electrical stimulation[END_REF]. We refer to [START_REF] Wilson | Force response of locust skeletal muscle[END_REF] for a general discussion about all these models.

In the Ding et al. force-fatigue model, the physical FES input is modeled by Dirac pulses and is integrated using a linear dynamics, which leads to a sampled-data control system, where the control parameters are the pulses amplitudes and the pulses times. In this paper, we consider the problem of minimizing a general cost function of Mayer form, which depends on the final values of the force response and of the fatigue variables. The resulting optimization problem is clearly related to optimal sampled-data control problems, but it does not fit exactly with the framework considered in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF], in which a Pontryagin maximum principle with first-order necessary optimality conditions is derived. As a consequence, our main objective in the present paper is to adapt the techniques of [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] to our specific problem in order to derive the corresponding Pontryagin-type conditions.

As illustration, we propose some preliminary numerical simulations in the context of maximization of the final force response (with fixed pulses amplitudes) determining the corresponding optimal pulses times. Precisely, we first implement a direct method using Bocop software [START_REF] Bonnans | Bocop -A collection of examples[END_REF], which allows us, in a second time, to initialize and implement an indirect method (shooting method) based on the Pontryagin-type conditions derived in our main result and using HamPath software [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF].

Let us mention that the present work is related to an industrial research project, whose aim is to design a smart electrical muscle stimulator. In this context, preliminary suboptimal numerical strategies were designed using a Model Predictive Control (MPC) coupled with online estimation of the fatigue variables using a nonlinear observer; see [START_REF] Bakir | A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model[END_REF]. Additional information on the practical issues about this project may be found in [START_REF] Yochum | Contribution à la conception d'un électromyostimulateur intelligent[END_REF][START_REF] Bakir | Contribution à la modélisation, l'estimation et la commande de systèmes non linéaires dans les domaines de la cristallisation et de l'électrostimulation musculaire[END_REF].

The paper is organized as follows. Section 2 is dedicated to a recap about the Pontryagin maximum principle for optimal control problems. We distinguish two situations: the permanent control case versus the sampled-data control case. After recalling briefly the Ding et al. force-fatigue model, Section 3 contains the main mathematical contribution of the present article. Precisely, the physical FES system is set into the framework of sampled-data control system and a general Mayer optimal control problem is considered. A discussion is initiated in order to explain why it does not fit exactly with the framework investigated in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF]. Then, we adapt the previous works to our particular problem and we establish the corresponding Pontryagin first-order necessary optimality conditions (see Theorem 3.1). Finally, we show in Section 4 how they can be used in view of numerical simulations.

Recap on the Pontryagin Maximum Principle: Permanent Control Case versus Sampled-Data Control Case

In this section, we propose a brief recap on the Pontryagin maximum principle for optimal control problems and we compare two situations: the permanent control case versus the sampled-data control case. We consider, in this section, a very simple framework (smooth dynamics, Mayer cost, fixed initial condition, no final state constraint) and we give some recalls on the main techniques (related to the classical calculus of variations [START_REF] Gelfand | Calculus of Variations[END_REF]) leading to each version of the Pontryagin maximum principle. Let T > 0 and let d, m ∈ N * be two fixed positive integers. Let us consider a general nonlinear control system of the form ẋ

(t) = f (t, x(t), u(t)), a.e. t ∈ [0, T ], (1)
where

f : [0, T ] × R d × R m → R d is of class C 1 , together with the fixed initial condition x(0) = x 0 ∈ R d and the control constraint u(t) ∈ U , where U ⊂ R m is a nonempty subset of R m .
We focus in this section on the general Mayer optimal control problem given by

min u∈U ϕ(x(T )), (2)
where ϕ : R d → R is a differentiable function and where U stands for the set of admissible controls (see details in each subsection below). Recall that the Hamiltonian function

H : [0, T ] × R d × R d × R m → R associated to Problem (2) is defined by H(t, x, p, u) := p, f (t, x, u) R d for all (t, x, p, u) ∈ [0, T ] × R d × R d × R m , where •, • R d stands for the standard scalar product in R d .

Permanent Control Case

If the set U of admissible controls in Problem ( 2) is the set of all essentially bounded measurable functions u : [0, T ] → U , then there is no restriction on the modification of the value of the control and thus it can occur at any time in [0, T ]. In such case, Problem ( 2) is said to be with permanent control. This situation corresponds to the very well-known framework deeply studied and developed in the literature (see, e.g., [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Trélat | Contrôle optimal: théorie et applications[END_REF][START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF] and references therein). From now, our aim is to briefly recall the derivation of the Pontryagin maximum principle for Problem (2) in the case of permanent control. Let x * be a reference optimal curve associated to the control u * . Take a L 1 -perturbation (or needlevariation) of u * defined by u ε (t

) := v ∈ U on [s, s + ε[, where s ∈ [0, T [ is a Lebesgue time of the function t → f (t, x * (t), u * (t))
, and u ε (t) := u * (t) elsewhere. The corresponding variation vector satisfies the linear equation

ẇ(t) = ∇ x f (t, x * (t), u * (t)) × w(t), a.e. t ∈ [s, T ],
where the notation × stands for the usual matrix-vector product (also matrixmatrix product in the sequel), and the initial condition

w(s) = f (s, x * (s), v) -f (s, x * (s), u * (s)).
From optimality, it holds that ϕ(xε(T ))-ϕ(x * (T )) ε ≥ 0, where x ε denotes the response to u ε . Taking the limit ε ↓ 0, one gets

∇ϕ(x * (T )), w(T ) R d ≥ 0. ( 3
)
Write the adjoint equation as

ṗ(t) = -∇ x f (t, x * (t), u * (t)) × p(t), a.e. t ∈ [0, T ], (4)
where stands for the standard transpose of a matrix, with the final condition

p(T ) = -∇ϕ(x * (T )). (5)
Using in (3) the equalities w(T ) = Φ(T, s) × w(s) and p(s) = Φ(T, s) × p(T ), where Φ(

can be identified to an absolutely continuous on [0, T ], which moreover satisfies

Ḣ(t) = ∇ t H(t, x * (t), p(t), u * (t)), a.e. t ∈ [0, T ].
In particular, if the control system (1) is autonomous, then H remains constant over [0, T ].

Sampled-Data Control Case

At the opposite of the permanent control case, if the set U of admissible controls authorizes the value of the control u : [0, T ] → U to be modified at most n times, where n ∈ N is fixed, Problem (2) is said to be with sampled-data control. In such case, for any u ∈ U, there exists a finite set of n times 0 < t 1 < . . . < t n < T (called sampling times) such that u(t) = u i ∈ U over [t i , t i+1 [ for all i = 0, . . . , n, where t 0 := 0 and t n+1 := T . We refer to [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] for the statement of a Pontryagin maximum principle handling sampled-data controls. From now, we assume that U is convex (see Remark 2.2 for details) and our aim is to briefly recall the derivation of the Pontryagin maximum principle for Problem (2) in the case of sampled-data control. Let x * be a reference optimal curve associated to the control u * and let us denote by t * i the corresponding sampling times. Consider the convex L ∞perturbation u ε := u * + ε(u -u * ), where u ∈ U has the same sampling times than u * . The corresponding variation vector satisfies the affine equation

ẇ(t) = ∇ x f (t, x * (t), u * (t)) × w(t) + ∇ u f (t, x * (t), u * (t)) × (u(t) -u * (t)), a.e. t ∈ [0, T ],
with the initial condition w(0) = 0 R d . Introducing the adjoint vector p as in ( 4)-( 5) and using the equality

w(T ) = T 0 Φ(T, s) × ∇ u f (s, x * (s), u * (s)) × (u(s) -u * (s)) ds, in Inequality (3) give T 0 p(s), ∇ u f (s, x * (s), u * (s)) × (u(s) -u * (s)) R d ds ≤ 0. Taking u(t) := v ∈ U over [t * i , t * i+1
[ and u(t) := u * (t) elsewhere, we exactly recover the nonpositive averaged Hamiltonian gradient condition derived in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] given by

t * i+1 t * i ∇ u H(s, x * (s), p(s), u * i ) ds , v -u * i R m ≤ 0, (8)
for all v ∈ U and all i = 0, . . . , n, where u * i corresponds to the value of u * over the interval [t i , t i+1 [. Remark 2.1 The Pontryagin conditions ( 6) and ( 8) are clearly different according to the nature of the admissible controls set U. We refer to the works [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] for detailed discussions about this feature.

Remark 2.2

The convexity assumption on U is used in the sampled-data control context in order to guarantee that the convex L ∞ -perturbation u ε is admissible. This assumption is made here for simplicity, but it can be removed by using the concept of U -dense directions (see, e.g., [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF]).

Remark 2.3

In Problem (2) with sampled-data control, the sampling times t i can be fixed or not. Note that the works [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF] focus only on fixed sampling times. A recent study can be found in [START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] about the differences between fixed and free sampling times. In the case of fixed sampling times, the maximized Hamiltonian function H defined in ( 7) is piecewisely smooth of class C 1 , but is not continuous on [0, T ] in general. It may admit a discontinuity at each sampling time t i . At the opposite, if the sampling times are free in Problem (2), then additional necessary optimality conditions are established and they exactly coincide with the continuity of H on the whole interval [0, T ]. We refer to [START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] for more details.

Main Result

Ding et al. Force-Fatigue Model

The FES input v (or pulses train) is of the form

v(t) := n i=0 η i δ(t -t i ), for all t ∈ [0, T ],
modeled as a finite sum of Dirac impulses at times 0 = t 0 < t 1 < . . . < t n < T , where n ∈ N and T > 0 are fixed, and where η i ∈ [0, 1] are the amplitudes of each pulse. To describe the phenomenon of tetania (that is, the memory effect of successive pulses), we introduce the scaling factor R i defined by

R i :=    1, for i = 0, 1 + ( R -1) exp - t i -t i-1 τ c , for i = 1, . . . , n,
where the constants τ c > 0 and R > 1 are given in Table 3 of the Appendix.

We define the FES signal E to be the solution in the distributional sense of the linear scalar dynamics

Ė(t) = - E(t) τ c + 1 τ c n i=0 R i η i δ(t -t i ), over [0, T ], (9)
with E(0) = 0. We get that

E(t) = 1 τ c n i=0 R i exp - t -t i τ c η i h(t -t i ), for all t ∈ [0, T ], (10)
where h stands for the left-continuous Heaviside function. The FES signal drives the evolution of Ca 2+ -concentration C N according to the linear scalar dynamics

ĊN (t) = - C N (t) τ c + E(t), a.e. t ∈ [0, T ]. ( 11
)
Integrating the (resonant) system [START_REF] Bonnans | Bocop -A collection of examples[END_REF] with C N (0) = 0 leads to

C N (t) = 1 τ c n i=0 R i exp - t -t i τ c η i (t -t i )h(t -t i ), for all t ∈ [0, T ].
Introducing the Hill functions β(t) :=

C N (t)
Km+C N (t) and γ(t) := 1 τ1+τ2β(t) for all t ∈ [0, T ], the non-fatigue model [START_REF] Law | Mathematical models of human paralyzed muscle after long-term training[END_REF][START_REF] Ding | Two-step, predictive, isometric force model tested on data from human and rat muscles[END_REF][START_REF] Ding | Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains[END_REF] describes the force response F by the dynamics

Ḟ (t) = -γ(t)F (t) + Aβ(t), for all t ∈ [0, T ], (12)
with F (0) = 0, where A, K m , τ 1 , τ 2 are given positive constants. We refer to Figure 1 for a summary diagram of the non-fatigue model ( 11)-( 12). The complete force-fatigue model [START_REF] Ding | Mathematical models for fatigue minimization during functional electrical stimulation[END_REF] is obtained by considering A, K m , τ 1 as fatigue variables following the linear dynamics

Linear Filter R i Linear Filter Static Nonlinearity Filter E(t) C N (t) β(t) F (t) i η i δ(t -t i )

Pulses train

Ȧ(t) = - A(t) -A rest τ f at + α A F (t), (13)
Km (t) = - K m (t) -K m,rest τ f at + α Km F (t), (14)
τ1 (t) = - τ 1 (t) -τ 1,rest τ f at + α τ1 F (t), (15)
for all t ∈ [0, T ], with the initial conditions A(0) = A rest , K m (0) = K m,rest , and τ 1 (0) = τ 1,rest . Table 3 of the Appendix contains the definitions and details on the variables and constants introduced in this section. We also refer to [START_REF] Wilson | Force response of locust skeletal muscle[END_REF] for a general discussion about these models.

A General Rewritting of the Model and Some Comments

In this paper, we consider the complete force-fatigue model ( 11)-( 12)-( 13)-( 14)-( 15) and we denote by x = (C N , F, A, K m , τ 1 ) . The model can be rewritten as the control system

ẋ(t) = g(x(t)) + b(t) n i=0 G(t i-1 , t i )η i h(t -t i ) e, a.e. t ∈ [0, T ], ( 16
)
with the initial condition x(0) = x 0 , where x 0 , e ∈ R 5 are the vectors defined by x 0 := (0, 0, A rest , K m,rest , τ 1,rest ) and e := (1, 0, 0, 0, 0) , where b(t) := 1 τc e -t τc for all t ∈ [0, T ] and G(t i-1 , t i ) := ( R-1)e ti-1/τc +e ti/τc , where t -1 := -∞, t 0 := 0 and t n+1 := T , and finally where the explicit expression of the function g : R 5 → R 5 is given by

g(x) :=         -x1 τc - x1+x4 x5(x1+x4)+τ2x1 x 2 + x 3 x1 x1+x4 -x3-Arest τ f at + α A x 2 - x4-Km,rest τ f at + α Km x 2 - x5-τ1,rest τ f at + α τ1 x 2         , for all x = (x 1 , x 2 , x 3 , x 4 , x 5 ) ∈ R 5 .
Remark 3.1 Due to the definition of the vector e, note that the control parameters η i and t i intervene explicitly in the control system [START_REF] Gelfand | Calculus of Variations[END_REF] only in the differential equation of the first component x 1 of x. Furthermore note that this first differential equation is linear and independent of the rest of the system. This remark is of particular interest when dealing with the sensitivity analysis of the control system (16) in Section 3.4.

In this section, we consider the general optimal control problem of Mayer form given by min ϕ(x(T ))

s.t.                              ẋ(t) = g(x(t)) + b(t) n i=0 G(t i-1 , t i )η i h(t -t i ) e, a.e. t ∈ [0, T ], x(0) = x 0 , (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ R 2n+1 , η i ∈ [0, 1], ∀i = 0, . . . , n, 0 < t 1 < t 2 < . . . < t n < T, t i -t i-1 ≥ I min , ∀i = 1, . . . , n, (17)
where n ∈ N and T > 0 are fixed, where ϕ : R 5 → R is a differentiable function and where I min ≥ 0 is the minimal interpulse authorized. The term ϕ(x(T )) plays the role of a general cost of Mayer form to minimize (for instance the opposite of the final force response, as illustrated in Section 4).

One should note that the above optimal control problem (17) does not fit exactly with the framework of optimal sampled-data control problems considered in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] (and recalled in Section 2.2) for two reasons. The first difference lies in the fact that the value of each amplitude η i intervenes in the control system over the interval [t i , T ] and not only on the interval [t i , t i+1 [. The second difference lies in the fact that the values of the sampling times t i intervene explicitly in the control system. As a consequence, the first-order necessary optimality conditions from [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF] cannot be applied and have to be adapted. This is the major aim of the next subsections. Precisely, after introducing the notions of admissible tuples of parameters and admissible perturbations in Section 3.3, we proceed in Section 3.4 to the sensitivity analysis of the control system [START_REF] Gelfand | Calculus of Variations[END_REF] with respect to the parameters η i and t i . Finally, the Pontryagin first-order necessary optimality conditions for Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF] are deduced in Section 3.5 (see Theorem 3.1).

Remark 3.2

In the model considered in this paper, the first time t 0 of pulse is fixed at t 0 = 0, while the corresponding amplitude η 0 ∈ [0, 1] is not fixed. As a consequence, in Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF], there are exactly 2n + 1 control parameters given by η 0 , η 1 , . . . , η n , t 1 , . . . , t n .

Admissible Tuples of Parameters and Perturbations

We introduce in this section the following notions of admissible tuples of parameters and admissible perturbations. Definition 3.1 (Admissible tuple of parameters) We say that the tuple δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ R 2n+1 of parameters is admissible for Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF] if

       η i ∈ [0, 1], ∀i = 0, . . . , n, 0 < t 1 < t 2 < . . . < t n < T, t i -t i-1 ≥ I min , ∀i = 1, . . . , n.
We denote by P the set of all admissible tuples of parameters for Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF].

Definition 3.2 (Admissible perturbations) Let us consider an admissible tuple δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P of parameters for Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF].

(i) Let i ∈ {0, . . . , n}. We say that ηi ∈ R is an admissible perturbation of η i if there exists ε > 0 such that η i + ε ηi ∈ [0, 1] for all 0 ≤ ε ≤ ε. (ii) Let i ∈ {1, . . . , n -1}. We say that ti ∈ R is an admissible perturbation of t i if there exists ε > 0 such that (t i + ε ti ) -t i-1 ≥ I min and t i+1 -

(t i + ε ti ) ≥ I min for all 0 ≤ ε ≤ ε. (iii) Let i = n. We say that tn ∈ R is an admissible perturbation of t n if there exists ε > 0 such that (t n + ε tn ) -t n-1 ≥ I min for all 0 ≤ ε ≤ ε. Remark 3.3 Let δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P and i ∈ {0, . . . , n}. If η i = 0 (resp. η i = 1), then the set of all admissible perturbations of η i is R + (resp. R -). If η i ∈]0, 1[
, then the set of all admissible perturbations of η i is the entire line R.

Remark 3.4 Let δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P and i ∈ {1, . . . , n -1}. If

t i -t i-1 = I min < t i+1 -t i (resp. t i -t i-1 > I min = t i+1 -t i )
, then the set of all admissible perturbations of t i is R + (resp. R -). If t i -t i-1 > I min and I min < t i+1 -t i , then the set of all admissible perturbations of t i is the entire line R. If t i -t i-1 = I min = t i+1 -t i , then the set of all admissible perturbations of t i is reduced to the singleton {0}.

Remark 3.5 Let δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P and i = n. If t n -t n-1 > I min (resp. t n -t n-1 = I min ), then the set of all admissible perturbations of t n is the entire line R (resp. R + ).

Sensitivity Analysis of the Control System and Variation Vectors

For all δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P, we denote by x δ the unique solution to the Cauchy problem

     ẋ(t) = g(x(t)) + b(t) n i=0 G(t i-1 , t i )η i h(t -t i ) e, a.e. t ∈ [0, T ], x(0) = x 0 .
Recall that the first component x δ 1 of x δ satisfies

x δ 1 (t) = b(t) n i=0 G(t i-1 , t i )η i (t -t i )h(t -t i ), (18)
for all t ∈ [0, T ] (see Section 3.1, in which x 1 δ = C N ). Proposition 3.1 (Perturbation of η i ) Let δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P and i ∈ {0, . . . , n}. Let ηi be an admissible perturbation of η i and let w : [t i , T ] → R 5 be the unique solution to the linear Cauchy problem given by

   ẇ(t) = ∇g(x δ (t)) × w(t) + ηi b(t)G(t i-1 , t i ) e, a.e. t ∈ [t i , T ], w(t i ) = 0 R 5 .
Then, it holds that

w(t) = lim ε↓0 x δ+ε δ (t) -x δ (t)
ε , for all t ∈ [t i , T ], where δ := (0, . . . , 0, ηi , 0, . . . , 0, 0

R n ) ∈ R 2n+1 .
Proof We reason component by component of x δ = (x δ 1 , . . . , x δ 5 ). Note that the first component w 1 of w = (w 1 , . . . , w 5 ) is the unique solution to the linear Cauchy problem given by

   ẇ1 (t) = -w1(t) τc + ηi b(t)G(t i-1 , t i ) , a.e. t ∈ [t i , T ], w 1 (t i ) = 0,
which is given by

w 1 (t) = ηi b(t)G(t i-1 , t i )(t -t i ) ,
for all t ∈ [t i , T ], which exactly corresponds to the derivative of x δ 1 with respect to η i in the direction ηi over the interval [t i , T ] (see Equality [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]). Then, denoting by g = (g 1 , . . . , g 5 ), the second component x δ 2 of x δ satisfies the differential equation

( ẋδ 2 )(t) = g 2 (x δ 1 (t), . . . , x δ 5 (t)), a.e. t ∈ [t i , T ],
which is, as noted in Remark 3.1, independent of the parameters η i and t i . As a consequence, from the previous result on the first component x δ 1 of x and from standard techniques of sensitivity analysis of differential equations (essentially by invoking the classical Gronwall lemma), we get that:

lim ε↓0 (x δ+ε δ 2 (t) -x δ 2 (t))/ε = w 2 (t) for all t ∈ [t i , T ]
, where w 2 is the unique solution to the linearized differential equation

ẇ2 (t) = ∇g 2 (x δ 1 (t), . . . , x δ 5 (t)) × w(t), a.e. t ∈ [t i , T ],
with the initial condition w 2 (t i ) = 0. The arguments are exactly the same for the components x δ 3 , x δ 4 and x δ 5 .

Proposition 3.2 (Perturbation of t i )

Let δ = (η 0 , η 1 , . . . , η n , t 1 , . . . , t n ) ∈ P and i ∈ {1, . . . , n}. Let ti be an admissible perturbation of t i and let w : [t i , T ] → R 5 be the unique solution to the linear Cauchy problem given by

           ẇ(t) = ∇g(x δ (t)) × w(t) + ti b(t)b(-t i ) η i + ( R -1)η i+1 h(t -t i+1 ) e, a.e. t ∈ [t i , T ], w(t i ) = -ti b(t i )G(t i-1 , t i )η i e.
Then, it holds that

w(t) = lim ε↓0 x δ+ε δ (t) -x δ (t) ε ,
for all t ∈]t i , T ], where δ := (0 R n+1 , 0, . . . , 0, ti , 0, . . . , 0) ∈ R 2n+1 .

Proof As in the proof of Proposition 3.1, we reason in this proof component by component of x δ = (x δ 1 , . . . , x δ 5 ). Note that the first component w 1 of w = (w 1 , . . . , w 5 ) is the unique solution to the linear Cauchy problem given by

           ẇ1 (t) = -w1(t) τc + ti b(t)b(-t i ) η i + ( R -1)η i+1 h(t -t i+1 ) , a.e. t ∈ [t i , T ], w 1 (t i ) = -ti b(t i )G(t i-1 , t i )η i ,
which is given by

w 1 (t) = -ti b(t)η i G(t i-1 , t i ) + b(t)b(-t i )η i (t -t i ) + b(t)b(-t i )( R -1)η i+1 (t -t i+1 )h(t -t i+1 ) ,
for all t ∈ [t i , T ], which exactly corresponds to the derivative of x δ 1 with respect to t i in the direction ti over the interval ]t i , T ] (see Equality [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]).

To see this, one may use the equalities

∂ 1 G(t i , t i+1 ) = b(-t i )( R -1) and ∂ 2 G(t i-1 , t i ) = b(-t i ).
We conclude this proof by invoking the same arguments as in the proof of Proposition 3.1 for the rest of the components of x δ = (x δ 1 , . . . , x δ 5 ).

First-Order Necessary Optimality Conditions

The main result of the present paper is stated as follows. It corresponds to the Pontryagin first-order necessary optimality conditions for Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF].

Theorem 3.1 Let x * stand for an optimal state of Problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF] and let δ * = (η * 0 , η * 1 , . . . , η * n , t * 1 , . . . , t * n ) ∈ P be the corresponding optimal tuple of parameters. Then, the adjoint vector p = (p 1 , . . . , p 5 ) (also called costate), defined as the unique solution to the backward linear Cauchy problem given by

ṗ(t) = -∇g(x * (t)) × p(t), a.e. t ∈ [0, T ],
p(T ) = -∇ϕ(x * (T )), satisfies:

(i) the inequality T t * i p 1 (s)b(s) ds ηi ≤ 0, ( 19
)
for all i = 0, . . . , n and all admissible perturbation ηi of η * i ; (ii) and the inequality

-p 1 (t * i )b(t * i )G(t * i-1 , t * i )η * i + b(-t * i )η * i T t * i p 1 (s)b(s) ds + b(-t * i )( R -1)η * i+1 T t * i+1 p 1 (s)b(s) ds ti ≤ 0, ( 20
)
for all i = 1, . . . , n and all admissible perturbation ti of t * i .

Proof In this proof, we denote by Φ(•, •) the state-transition matrix associated to the matrix function t → ∇g(x * (t)). Let us prove the first item of Theorem 3.1. Let i ∈ {0, . . . , n} and ηi be an admissible perturbation of η * i . From optimality of x * and following the same general idea detailed in Section 2, we derive that p(T ), w(T ) R 5 ≤ 0, where the variation vector w is defined in Proposition 3.1. Since p(t) = Φ(T, t) × p(T ) for all t ∈ [0, T ] and

w(T ) = T t * i Φ(T, s) × b(s)G(t * i-1 , t * i )e ds ηi ,
we deduce that Inequality [START_REF] Trélat | Contrôle optimal: théorie et applications[END_REF] holds true. To see this, one may use the equality p 1 (t) = p(t), e R 5 for all t ∈ [0, T ]. Let us prove the second item of Theorem 3.1. Let i ∈ {1, . . . , n} and ti be an admissible perturbation of t * i . From optimality of x * and following the same general idea detailed in Section 2, we derive that p(T ), w(T ) R 5 ≤ 0, where the variation vector w is defined in Proposition 3.2. Since p(t) = Φ(T, t) × p(T ) for all t ∈ [0, T ] and

w(T ) = -Φ(T, t * i ) × b(t * i )G(t * i-1 , t * i )η * i e + T t * i Φ(T, s) × b(s)b(-t * i )η * i e ds + T t * i+1 Φ(T, s) × b(s)b(-t * i )( R -1)η * i+1 e ds ti ,
we deduce that Inequality (20) holds true. To see this, one may use the equality p 1 (t) = p(t), e R 5 for all t ∈ [0, T ]. The proof is complete.

Remark 3.6 As mentioned in Section 3.2, the optimal control problem [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF] does not fit exactly with the framework of optimal sampled-data control problems considered in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF][START_REF] Bourdin | Pontryagin Maximum Principle for optimal sampled-data control problems[END_REF][START_REF] Bourdin | Continuity/constancy of the Hamiltonian function in a Pontryagin maximum principle for optimal sampled-data control problems with free sampling times[END_REF]. As a consequence, we had to adapt in Sections 3.3 and 3.4 the techniques in order to derive our main result Theorem 3.1. It appears that introducing and dealing with a very general framework, allowing to cover in particular the Ding et al. force-fatigue model considered in this paper, constitute an interesting challenge for future works.

Preliminary Numerical Results

Our objective here is to illustrate our main result with some preliminary numerical simulations. In this section, we will focus on the problem of maximizing the final force response F (T ), with fixed pulses amplitudes η i = 1, for all i = 0, . . . , n, and with no minimal interpulse (that is, with I min = 0). In that context, we will thus focus only on the necessary optimality condition [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF] provided in Theorem 3.1. We denote by σ = (t 1 , . . . , t n ) the n-tuple of pulses times and we define

NC i (σ) := -p 1 (t i )b(t i )G(t i-1 , t i ) + b(-t i ) T ti p 1 (s)b(s) ds + b(-t i )( R -1) T ti+1 p 1 (s)b(s) ds,
for all i = 1, . . . , n. Assuming that the considered optimization problem admits a solution σ * = (t * 1 , . . . , t * n ), and since the set of admissible perturbations of each optimal pulses time t * i is the entire real line R, Inequality (20) reduces to the equality

NC i (σ * ) = 0, (21)
for all i = 1, . . . , n. Using the notations introduced in Section 3.2, let us recall that the force variable is denoted by F = x 2 and let us precise that it will be expressed in milliNewton (mN). The unit of time will be the second (s).

We implement the combination of several numerical methods. First, we use in Section 4.1 a direct method in order to compute a numerical approximation σ = ( t1 , . . . , tn ) of σ * using Bocop software [START_REF] Bonnans | Bocop -A collection of examples[END_REF], which takes into account the constraint t i ≤ t i+1 . Indeed, note that the numerical inversion of the elements t i is not allowed in our context since they explicitly intervene in the dynamics ( 16) and they do not play symmetrical roles. Then, we compute the initial costate p(0) using specific and suitable numerical integrators from the Julia programming language [START_REF] Rackauckas | DifferentialEquations.jl -A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia[END_REF]. Finally, the couple (σ, p(0)) is used in Section 4.2 as initialization of an indirect method (shooting method) based on the necessary condition (21) and using HamPath software [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] (which does not take into account the constraint t i ≤ t i+1 ) in order to numerically compute the optimal pulses times t * i .

Numerical Computation of Optimal Pulses Times With Bocop to Initialize the Shooting Method

Bocop software [START_REF] Bonnans | Bocop -A collection of examples[END_REF] is mainly dedicated to solve optimal permanent control problems by implementing a so-called direct transcription method. Namely, a time discretization is used in order to rewrite the optimal permanent control problem as a finite-dimensional optimization problem (i.e. nonlinear programming), solved by the interior point method (Ipopt). It is also designed to solve finite-dimensional optimization problems, where there is a finite number of parameters, in which the optimization problem considered in this section is plugged into. For our problem, this method is robust with respect to the number n of parameters (n < 20) and to T (T < 1s).

Once we have determined a numerical approximation σ = ( t1 , . . . , tn ) of σ * with Bocop , we compute the corresponding initial costate p(0), which will be used in order to initialize our shooting method in the next section. The numerical integration of the costate p is numerically intricate when T > 0.2, n > 4 and we need low tolerances (< 10 -20 ) and suitable ODE solvers using specific algorithms of the Julia programming language [START_REF] Rackauckas | DifferentialEquations.jl -A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia[END_REF]. For this reason, we use the package DifferentialEquations.jl, namely, a 4th order A-stable stiffly stable Rosenbrock method with a stiff-aware 3rd order interpolant (Rodas42), coupled with the package ArbNumerics.jl for multiprecision numerical computing. The numerical results of three cases are detailed in Table 1. We also provide the time evolutions of the corresponding state x and costate p (and also of the associated FES signal Ẽ) in Figures 2 and3 for the last two cases.

Note that we are able to evaluate the values NC i (σ) for all i = 1, . . . , n and, in order to check the correctness of the initialization couple (σ, p(0)), we represent in Figures 4 and5 the values of (NC 1 (σ), . . . , NC n (σ)) R n in relation with the final force response x 2 (T ) for many perturbations σ = σ + εχ of σ, where ε > 0 is a small parameter and χ is a uniform random variable in [-1, 1] n . It shows (see Figure 4) that the sampling times σ computed with the direct method yield the one with maximum cost and for which the quantity NC 1 (σ), . . . , NC n (σ) is the smallest among the perturbed sampling times σ + εχ. Such σ is a good initialization for the shooting algorithm (the latter converges typically in less than 10 iterations). Indirect methods are generally more precise than direct methods, but they are also very sensitive with respect to the initialization. Given T and n, Bocop yields many different sampling times (t 1 , . . . , t n ) and one needs to compute 1b). Time evolutions of state x, costate p and FES signal Ẽ.

these sampling times accurately. Although Figures 4 and5 give a rough idea of the quality of the sampling times, we check the accuracy of the solutions from Bocop by combining direct and indirect method. We describe here an indirect method (shooting method) in order to solve the optimization problem considered in this section.

Our shooting method is based on Equality (21) and computes numerically the optimal solution σ * by solving a boundary value problem associated to the state-costate system. Precisely, we are looking for a zero of the shooting function defined by S : R n+5 -→ R n+5 (σ, p 0 ) -→ NC 1 (σ), . . . , NC n (σ), p(T ) -p T , where p 0 ∈ R 5 is the initial condition considered in the adjoint system and p T = (0, 1, 0, 0, 0) (because we minimize -x 2 (T )). We use a Newtonlike algorithm in HamPath software [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] in order to compute a zero of S. In order to initialize this method, we need a guess on the optimal tuple σ * and on the initial value of the costate vector. These guesses are provided by the direct method using Bocop (see Section 4.1). Note that a good initialization is crucial in particular because we are not able to take into account the constraint t i ≤ t i+1 with the shooting method.

Figure 6 (resp. Figure 7) represents the time evolutions of the optimal state and costate (and also of the associated FES signal E) computed with 1c). Time evolutions of state x, costate p and FES signal Ẽ. the above shooting method, which was initialized by the values provided in Table 1a (resp. Table 1b). In both cases, the shooting method recovers the solutions found by the direct method (see Table 2). 1c and Figure 3). Plot of the norm (NC1(σ), . . . , NCn(σ)) R n with respect to the final force response x2(T ) for 600 perturbations σ of σ. The red dot indicates the solution σ computed by Bocop.

We obtain preliminary numerical results of optimal pulses times in the context of maximization of the final force response (with fixed pulses amplitudes and no minimal interpulse) based on a shooting method implemented with the HamPath software. Since the Newton-like algorithm is very sensitive with respect to its initialization, the shooting method is initialized using a direct method based on the Bocop software and some Julia's packages for differential equations and multiprecision computing.

Further studies have to be conducted in many directions. First, calculations should be obtained in the frame of geometric control described in [START_REF] Hermes | Lie algebras of vector fields and local approximation of attainable sets[END_REF] using Baker-Campbell-Hausdorff formulae and Lie brackets. Second, a research perspective is to provide a theoretical and numerical study of the shooting equations (with possibly optimized pulses amplitudes and in presence of a positive minimal interpulse, with a large final time T > 0 and/or a large number n of pulses times, etc.). In particular, comparisons should be provided with respect to suboptimal numerical strategies [START_REF] Bakir | A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model[END_REF] evoked in Introduction, which concern closed loop controls. Finally, another challenge should be to take into account state constraints on the fatigue variables [START_REF] Bourdin | Optimal sampled-data controls with running inequality state constraints -Pontryagin maximum principle and bouncing trajectory phenomenon[END_REF].

Appendix: Physical Descriptions of the Force-Fatigue Model

Table 3 provides the values of the constants and the physical descriptions of the variables for the force-fatigue model presented in Section 3.

Fig. 1 :

 1 Fig. 1: Non-fatigue model description from the FES input v to the output force response F .

Fig. 2 :

 2 Fig. 2: Case n = 5, T = 0.5 (see Table1b). Time evolutions of state x, costate p and FES signal Ẽ.

Fig. 3 :

 3 Fig. 3: Case n = 10, T = 0.8 (see Table1c). Time evolutions of state x, costate p and FES signal Ẽ.

  (a) Case n = 4, T = 0.2, x(T ) = (0.3308, 0.2752, 2.996, 0.1036, 0.0516), p(0) = (0.0210, 0.2614, 0.0914, -0.6225, 0.8796).

  i

Fig. 4 :

 4 Fig. 4: Case n = 5, T = 0.5 (see Table1band Figure2). Plot of the norm (NC1(σ), . . . , NCn(σ)) R n with respect to the final force response x2(T ) for 600 perturbations σ of σ. The red dot indicates the solution σ computed by Bocop.

Fig. 6 :

 6 Fig. 6: Case n = 4, T = 0.2. Time evolutions of state, costate and FES signal computed from HamPath. We have: σ * = (0.0484, 0.0895, 0.1257, 0.1501).

Fig. 7 :

 7 Fig. 7: Case n = 5, T = 0.5. Time evolutions of state, costate and FES signal computed from HamPath. We have: σ * = (0.2258, 0.3028, 0.3616, 0.4093, 0.4470).

- 4 .0 10 - 1 1 s•mN 1 .9 10 - 2 1 mN 2 .1 10 - 2

 41011110212102 Coefficient for the force-model parameter A in the fatigue model α Km Coefficient for the force-model parameter Km in the fatigue model ατ 1 Coefficient for force-model parameter τ 1 in the fatigue model τ f at s 127 Time constant controlling the recovery of (A, Km, τ 1 )

Table 1 :

 1 Numerical results from Bocop in three cases. × 10 -4 , 3.1 × 10 -3 , 0.0959, -0.8293, 1.246). Indirect Method Using HamPath and Initialized by Bocop

	(a)	Case n = 4, T = 0.2, x(T ) = (0.3308, 0.2752, 2.996, 0.1036, 0.0516),
	p(0) = (0.0210, 0.2614, 0.0914, -0.6224, 0.8795).
		i	1			2		3		4
		ti =	0.0485		0.0895		0.1248		0.1499
		NC i (σ) =	-0.010		-0.007		-0.008	-0.003
	(b)	Case n = 5, T = 0.5, x(T ) = (0.2652, 0.2890, 2.984, 0.104, 0.0522),
	p(0) = (4.6 i	1		2		3	4		5
	ti =	0.241		0.314		0.369	0.414	0.448
	NC i (σ) =	-0.06		-0.03		-0.01	-0.02	-0.007
	(c)	Case n = 10, T = 0.8, x(T ) = (0.297, 0.354, 2.95, 0.106, 0.0538),
	p(0) = (-1.51 × 10 -4 , -2.03 × 10 -3 , 0.118, -0.98, 1.91).
	i	1	2	3	4	5	6	7	8	9	10
	ti =	0.273 0.363 0.434 0.498 0.554 0.602 0.645 0.684 0.719 0.748
	NC i (σ) = 0.026 0.016 0.036 0.028 0.018 0.026 0.033 0.018 0.034 0.018
	4.2

Table 2 :

 2 Numerical results obtained by HamPath in the cases illustrated by Figures6 and 7.

Table 3 :

 3 Values of the constant parameters in the Ding et al. model Value of the parameter A when muscle is not fatigued Km,rest -0.103 Value of the parameter Km when muscle is not fatigued τ 1,rest s 0.05095 The value of the parameter τ 1 when muscle is not fatigued α A

	Symbol Unit Value	Description
	C N	--	Normalized amount of Ca 2+ -troponin complex
	F	mN -	Force generated by muscle
	t i	s	-	Time of the i th pulse
	n	-	-	Total number of the pulses before time t
	i	-	-	Stimulation pulse index
	τc R	s -1.143 0.02	Time constant that commands the rise and the decay of C N Term of the enhancement in C N from successive stimuli
	A	mN s	-	Scaling factor for the force and the shortening velocity of muscle
	τ 1	s	-	Force decline time constant when strongly bound cross-bridges
				absent
	τ 2	s	0.1244	Force decline time constant due to friction between actin and
				myosin
	Km	-	-	Sensitivity of strongly bound cross-bridges to C N
	Arest	mN s	3.009
		1
		s 2

Acknowledgements This research paper benefited from the support of the FMJH Program PGMO and from the support of EDF, Thales, Orange. T. Bakir, B. Bonnard and J. Rouot are partially supported by the Labex AMIES.