Direct and Indirect Methods to Optimize the Muscular Force Response to a Pulse Train of Electrical Stimulation

. Recent force-fatigue mathematical models in biomechanics [7] allow to pre-dict the muscular force response to functional electrical stimulation (FES) and leads to the optimal control problem of maximizing the force. The stimulations are Dirac pulses and the control parameters are the pulses amplitudes and times of application, the num-ber of pulses is physically limited and the model leads to a sampled data control problem. The aim of this article is to present and compare two methods. The ﬁrst method is a direct optimization scheme where a further reﬁned numerical discretization is applied on the dynamics. The second method is an indirect scheme: ﬁrst-order Pontryagin type necessary conditions are derived and used to compute the optimal sampling times.


Introduction
Optimized force response to FES is an important problem for muscular reeducation and in case of paralysis. An historical model is known as the Hill model [8] and more refined models are taking into account the muscular fatigue ( [3], [7], [9]). See [12] for a comparison of the models. In this article, we shall consider the Ding et al. force-fatigue model. The physical control amounts to apply on [0, T ] a finite sequence of Dirac pulses at times 0 = t 1 < t 2 < . . . < t n < T , with amplitudes η i ∈ [0, 1], i = 1, . . . , n and in the model they are integrated using a linear dynamics to produce the so-called E s input (see Fig.1), which drives the force response. This fits in the sampled-data control frame since between each interpulse t i − t i−1 , i = 1, . . . , n the control is constant. Our aim is to optimize a cost function to maximize the force response. The control parameters are the interpulse t i − t i−1 , i = 2, . . . , n and the Dirac amplitudes η i , i = 1, . . . , n. There is a minimal interpulse due to the digital control but the dynamics can be discretized over a refined interval and this led to a so-called direct optimization problem which can be handled with an optimization routine. Another method is an indirect optimal control scheme. The optimal problem fits in the frame of sampled-data control problem and some Pontryagin type necessary optimality conditions were obtained in [6] and they were refined to be applied to our specific problem [2]. This leads to numerical methods to compute the optimal control in which we can distinguish between control by interpulses and control by amplitudes. The objective of this article is to present and to give preliminary numerical simulations.

Linear Filter
Pulses train The article is organized as follows. In section 1, the Ding et al. model is presented and the sampled data optimal control problem is introduced. Using a numerical discretization method, the direct optimization scheme is presented. In section 2, we give the Pontryagin type necessary optimality condition for sampled-data optimal control problem, derived for our category of problem and coming from [2]. Both sections contain numeric simulations.
1. The mathematical model and the direct scheme  [7]. The control associated to electrical pulses is described by applied at times t 1 = 0 < t 2 < ... < t n < T with train T := train duration (fixed) and η i ∈ [0, 1] are the amplitudes of each pulse. The FES signal denoted by E s (t) is the solution of the linear dynamicsĖ with E s (0) = 0 and R i is the parameter describing the phenomenon of tetania defined by and τ c is a parameter. In the model, the input E s (t) is the control to drive the force response using the dynamics where C N is the concentration of Ca 2+ , F is the force, A, K m , τ 1 , τ 2 are parameters and are the Hill functions. The triplet (A, K m , τ 1 ) contain the fatigue variables following the linear dynamics: with (A rest , K m,rest , τ 1,rest ) being the equilibrium, while τ f at , α A , α Km and α τ1 are parameters.
The force-fatigue model can be written shortly as: with x = (C N , F, A, K m , τ 1 ), and We use the convention t 0 = −∞, t n+1 = T and g(x) is defined by Eq. (4) and Eq. (6). Note that once the sampling times are fixed, the control is constant on each [t i−1 , t i ] and we get a sampled control system. The optimal sampled-data control problem that we want to study is the Mayer problem of the form: min u(·) ϕ(x(T )) where ϕ : R 5 → R is the cost function and the minimum is taken over the set (η 1 , η 2 , · · · , , η n , t 2 , · · · , t n ) ∈ R 2n−1 with the constraints η i ∈ [0, 1], 0 = t 1 < t 2 < · · · < t n < T = t n+1 and the interpulse constraints t i+1 − t i ≥ I m , i = 1, . . . , n.

Direct method
Direct methods can solve efficiently optimal control problems in the permanent control case and for the sampled-data control case. In the permanent control case, the optimal control problem is transformed as a non linear finite dimensional optimization problem and this is done by a discretization in time of the state and control variables in the dynamics. For the sampled-data control problems, the dynamics is discretized in time.

Direct method using Maltab
We present the general frame (stimulation time and amplitude as control variables). To control the force and/or the fatigue levels, three cases are considered (t 1 = 0 and T = t n+1 being fixed): • Stimulation times t i as control variables (with fixed amplitudes η i ): • Stimulation amplitudes (η i = η(t i )) as control variables (with fixed stimulation times t i ): with η 1 = η(t 1 ) and t 1 = 0.
• Stimulations times and amplitudes as control variables: In the third case (stimulations times and amplitudes as control variables), and from (7) and (8), we define ζ n : which, using Euler discretizing scheme, becomes where dt i is the discretizing step: To maximize the force level, the problem is written as: min ϕ(T ) subject to equality constraints (13) and linear inequality constraints Au ≤ B.
In the case stimulation times as control variables, stimulation amplitudes are fixed for i = 1, 2, . . . , N . Stimulation times constraints are: In the case of stimulation amplitudes as control variables, t 1 , T and t i for i = 2, . . . , N are fixed. Stimulation amplitudes constraints are 0 ≤ η i ≤ 1. Simulations results following different optimization methods are presented in Table 1 Table 1). Time evolutions of state and FES signal E s obtained by Matlab.

Direct method using Bocop
We use the Bocop software [4] to implement a direct approach, where the NLP is solved using a primal-dual interior point algorithm and the derivatives required for the optimization are computed by automatic differentiation with CppAD, this is crucial when dealing with the singularity of the Heaviside function appearing in the FES signal E s .
We give in Fig.4 the results with Bocop for n = 8, I m = 10ms and T = 400ms. Matlab toolbox gives 315.6N as a maximum force versus Bocop result 339.3N.

Pontryagin maximum principle with sampled-data control and necessary optimality conditions, indirect method
In this section we propose a first brief recap on the Pontryagin maximum principle for optimal control problems and we compare two situations: the permanent control case versus the sampleddata control case. We consider in this section a very simple framework (smooth dynamics and Mayer cost, fixed initial condition and no final state constraint) and we give some recalls on the main techniques (related to the classical calculus of variations) leading to each version of the Pontryagin maximum principle. Let T > 0 and let d, m ∈ N * be fixed positive integers. Let us consider a general nonlinear control system of the forṁ 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) ∈ Ω ⊂ R m . We focus in this section on the Mayer optimal control problem given by   Table 2).

Permanent control case
If the set U of admissible controls in Problem (15) is the set of all bounded measurable functions u : [0, T ] → Ω, 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 (15) 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., [10] and references therein). From now our aim is to briefly recall the derivation of the Pontryagin maximum principle for Problem (15) in the case of permanent control. Let x * be a reference optimal curve associated to the control u * . Take a L 1 -perturbation (or needle-variation) defined by u ε (t) = v ∈ Ω on [s, s+ε), where s 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 equatioṅ with 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 Write the co-state equation asṗ with the final condition p(T ) = −∇ϕ(x * (T )). Using in (16) the equalities w(T ) = Φ(T, s)×w(s) and p(s) = Φ(T, s) × p(T ), where Φ(·, ·) stands for the state-transition matrix associated to the matrix function t → ∂f ∂x (t, x * (t), u * (t)), one gets the inequality p(s), f (s, x * (s), v) − f (s, x * (s), u * (s)) ≤ 0 which corresponds exactly to the standard Hamiltonian maximization condition of the Pontryagin maximum principle given by We conclude this section by recalling that the maximized Hamiltonian function H : [0, T ] → R defined by H(t) = H(t, x * (t), p(t), u * (t)), can be shown to be absolutely continuous on [0, T ] witḣ In particular, if the control system (14) is autonomous, then H(·) is 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 ] → Ω to be modified at most n − 1 times, where n ∈ N * is fixed, the problem (15) 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 2 < . . . < t n < T (called sampling times) such that u(t) = u i ∈ Ω on each interval [t i , t i+1 ). We refer to [5,6] for the statement of Pontryagin maximum principle handling sampled-data controls. From now we assume that Ω is convex and our aim is to briefly recall the derivation of the Pontryagin maximum principle for Problem (15) in the case of sampleddata controls. 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) = A(t) × w(t) + B(t) × (u(t) − u * (t)), a.e. t ∈ [0, T ], with the initial condition w(0) = 0 R d , where A(t) = ∂f ∂x (t, x * (t), u * (t)) and B(t) = ∂f ∂u (t, x * (t), u * (t)). Introducing the co-state vector p as in the previous section and using the equality ) and u(t) = u * (t) elsewhere, we exactly recover the nonpositive averaged Hamiltonian gradient conditions derived in [5,6]   the initial co-state vector p(0) because p(T ) can be fixed to (0, 1, 0, 0, 0) ). Note that this method is indirect since we do not optimize any cost: the computations of the optimal sampling times rely on the necessary conditions (20). It differs from the direct method presented in Section 1.2.2. The inequalities (20) are handled easily introducing new state variables, e.g., to compute the integral I := T ti p 1 (s)b(s) ds, we can consider the variable y such thatẏ(t) = H(t ≥ t i ) p 1 (t)b(t) and we have y(T ) = I.
The finite-dimensional problem is solved using an interior point method or a stable trust-region method (without any cost to minimize) coupled with an auto-differentiation method to compute the derivatives. We represent in Fig.5 the solution obtained in the case where n = 8 and T = 400ms. The final force F (T ) =337.8 N is close to the one computed by the direct method (F (T ) =339.3 N) and the optimal sampling times for both methods are given in Table 2.
This is a different indirect method than the one presented in [2], where the indirect method was based on a shooting algorithm and was not able to handle constraints on the sampling times.

Conclusion
In this article we present and compare two methods (direct and indirect schemes) to optimize the FES-input in the Ding et al. force-fatigue model to maximize the force response at the end of the train. We gave preliminary numerical simulations proving the ability in both cases to implement the two methods. The same techniques can be used to consider Mayer problems where the cost takes into account the fatigue variables. This approach fits in the frame of openloop optimal controls problems. For the application since only the force is observed one must use an adaptive algorithm Table 2. Numerical values for the optimal sampling times obtained by the direct approach (Bocop) and the indirect approach and for the necessary conditions Θ i , i = 2, . . . , n in the case: n = 8, T = 400ms and I m = 10ms. like MPC method [1] where we optimize over a larger time with several pulses trains and rest periods and where the fatigues variables are online estimated using an observer. The challenge is to couple an optimization technique with on-line estimation. Our study in this frame allows to compute the maximal response to a train. In the direct case, BOCOP leads to better results than Matlab toolbox. In the indirect case, a nice contribution is to handle the inequalities (20), which is important to investigate, in the future, first and second order optimality conditions.