{'fraction_non_alphanumeric': 0.04205607476635514, 'fraction_numerical': 0.011828271028037383, 'mean_word_length': 4.193377148634985, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite Dimensional Approximation to Muscular Response in Force-Fatigue Dynamics using Functional Electrical Stimulation

Recent dynamical models, based on the seminal work of V. Hill, allow to predict the muscular response to functional electrostimulation (FES), in the isometric and non-isometric cases. The physical controls are modeled as Dirac pulses and lead to a sampled-data control system, sampling corresponding to times of the stimulation, where the output is the muscular force response. Such a dynamics is suitable to compute optimized controls aiming to produce a constant force or force strengthening, but is complex for real time applications. The objective of this article is to construct a finite dimensional approximation of this response to provide fast optimizing schemes, in particular for the design of a smart electrostimulator for muscularreinforcement or rehabilitation. It is an on-going industrial project based on force-fatigue models, validated by experiments.Moreover it opens the road to application of optimal control to track a reference trajectory in the joint angular variable to produce movement in the non-isometric models.


Introduction
Based on the seminal work of V. Hill [10], recent mathematical models (validated by experiments) allow to predict the force response to external stimulation. They are presented and discussed in details in [19] in the non-fatigue isometric case. They were extended in particular by Ding et al. [7,8,9] to take into account the muscular fatigue due to a long stimulation period and later in [14] to analyze the joint angular variable response, in the non-isometric case aiming to produce movements. Such models contain two basic nonlinearities to model the complexity of the dynamics. Fist of all, the ionic conduction and the nonlinear effect of successive pulses on the Ca 2+ -concentration. Second, the nonlinear dynamics relating the muscular force response to such concentration, modeled by the Michaelis-Menten-Hill functions [15].
For each train of pulses, due to digital constraints, only a finite number of pulses can be applied and from the optimal control point of view, the problem fits into the frame of optimal sampled-data control problems, studied in particular in [4] to derive Pontryagin necessary conditions. They can be analyzed to determine optimized train pulses and compared with direct optimizing schemes. A previous series of articles described the optimal control problems related to track a reference force or force strengthening, the control being either the interpulse I i = t i − t i−1 between two successive pulses or the amplitude of each pulse. In particular, model predictive control (MPC) method is presented in [1] aiming the use of online optimized closed loop control in the applications using force-fatigue model, as suggested in [6]. Direct methods vs indirect methods based on Pontryagin type necessary conditions are discussed and numerically implemented in [2] for the isometric case or in [3] for the non-isometric case.
The conclusion of aforementioned articles is that the nonlinear dynamics is computationally expensive in the numerical integration procedure and a challenging task is to reduce this time for real time computation in the applications. This article is motivated by the design of a smart electrostimulator, where the Ding et al. model is used to adjust automatically the frequency and the amplitude of the stimulations and to compute the sequence of stimulations and rest periods adapted to the task of the training program, e.g. endurance program or force strengthening program. The objective of this article being to bypass the computational difficulty by constructing a finite dimensional approximation of the force response, depending upon the parameters of each individual, which can be online estimated, aiming a real time computation of the optimized amplitudes and times, for each training program. Note that this approximation has been coded and the application scheme to the smart electrostimulator is briefly presented in the final section.
The article is organized as follows. In section 2, the mathematical model called the Ding et al. model [7,8,9] is presented and the main properties of the dynamics are described, reflecting the features of the muscular activity. Hence our analysis can be applied to different models discussed in [19]. The section 3 presents the optimization problems, in relation with muscular dynamics and oriented towards the design of a smart electrostimulator, where each training program must be translated into an optimization problem. Section 4 is the technical contribution of this article, that is the construction of the approximation for real time computation. In section 5, we present some numerical simulations aiming to validate the approximation and the optimizing scheme. In the final section 6, we outline the application to the design of the smart electrostimulator. It is based on a nonlinear output tracking [11,12,13] as a general theoretical frame and is applied to produce a constant force in our situation, but it can be extended to the non-isometric case to obtain any reference force. The conclusion indicates directions to complete our analysis, related to online parameters estimation of the problems [17,19] and MPC-methods [16,18] suitable for practical applications.

Mathematical model and main properties
We present the Ding et al. model force-fatigue model [7,8,9], extension of the original Hill model [10]. [7,8,9] The FES input u over a pulse train [0, T ] is given by

Ding et al. force-fatigue model
where 0 = t 0 < t 1 < · · · < t n < T are the impulsion times with n ∈ N being fixed and η i being the amplitudes of each pulses, which are convexified by taking η i ∈ [0, 1], δ(t − t i ) denoting the Dirac function at time t i .
Such physical control will provide the FES-signal denoted by E(t), which drives the force response using electrical conduction and its dynamics is given bẏ with E(0) = 0, depending upon the time response τ c and the scaling function R i defined by which codes the memory effect of successive muscle contractions and is associated to tetanus [19].
The first result is: where H is the Heaviside function and the FES signal depends upon two parameters (τ c ,R).
Definition 2 Consider a control system of the form: It is said permanent if u is a measurable bounded mapping valued in U . It is called a sampled-data control system if the set of controls is restricted to the set of piecewise constant mappings [u 0 , u 1 , . . . , u n ], u i ∈ U over a set of times t 0 = 0 < t 1 < · · · < t n < T , where n is a fixed integer.
Our problem can be formulated in the sampled-data control frame. One can write from (1), where u i (t) is the effect of the pulse η i δ(t − t i ) on the linear dynamics (2). One introduces the following.
is called the head and the restriction to [t i+1 , T ] is called the tail.
Clearly the FES-input is in the generalized frame of sampled data control system, provided we take into account the time-dependence and the phenomenon of tetanus. Observe also that each impulse have an effect on the whole train [0, T ].
The FES signal drives the evolution of the electrical conduction according to the linear dynamics describing the evolution of Ca 2+ -concentration c N :ċ and integrating the (resonant) system with c N (0) = 0 yields the following: which are the superposition of lobes of the form whose restriction to [t i , t i+1 ] forms the head of the corresponding lobe.
PROOF. Apply a time translation to the initial lobe with R 0 = 1.
Introducing the functions , m 2 (t) = 1 where m 1 is the Michaelis-Menten-Hill function [15], the force response satisfies the Hill dynamicṡ and where A, K m , τ 1 , τ 2 being additional parameters and we denote by Λ = (R, τ c , A, K m , τ 1 , τ 2 ) the whole set of parameters.
The model can be extended to take into account the fatigue. Using sensitivity analysis from [3], we shall restrict our study to the case of the force-fatigue Ding et al. model with the single equation: for all t ∈ [0, t f ], where t f is the total time and A(0) = A rest corresponds to the fixed value of A for the non fatigue model. This leads to introduce additional parameters τ f at , α A . Typical parameters values used in this article to validate numeric simulations are reported in Table 1 .

Mathematical rewriting
For the previous force-fatigue model and for the sake of the analysis, the model is rewritten as the control system: This leads to a control system of the forṁ The variable σ = (t 1 , . . . , t n , η 0 , η 1 , . . . , η n ) denotes the finite dimensional input-space with the constraints where I min is the smallest admissible interpulse.
Moreover the control is observed using the following observation function and h : x → (F, A) serves as a direct measure of the muscular force and the fatigue variable.
The following properties are straightforward but crucial in our analysis.

Proposition 6
For the non fatigue model, the force response can be integrated up to a time reparameterization as PROOF. Hill dynamics (7) is rewritten as and this linear dynamics can be integrated using Lagrange formula with F (0) = 0. This proves the assertion. 3 Optimization problems related to the design of the electrostimulator

Standard electrostimulators vs smart electrostimulators
The standard commercial electrostimulators apply a sequence of pulses trains and rest periods, where on each train [0, T ] the user only imposes the amplitude of the pulses trains and the constant frequency is related to training program, typically low frequency for endurance program and high frequency for force strengthening program. Our aim is to introduce optimization problems related to the design of a smart electrostimulator, which will be discussed in Section 6.

The punch program
In this case, our aim is to optimize the force at the end of the train over each train [0, T ]. This leads to: In this case, the amplitudes can be held at the constant maximal values η i = 1, i = 0, . . . , n and the optimization variables are the impulse times: 0 = t 0 < t 1 < · · · < t n < T. Since one considers a single train, the force model is sufficient.

The train endurance program
We consider a single train [0, T ] on which the model is the force model and the corresponding problem is Here, the amplitudes are appended to the impulse times to form the optimization variables and we use the convexified amplitudes constraints: η i ∈ [0, 1], i = 0, . . . , n.
The force reference has to be adjusted in relation with the user and can be set to F max /k, where k is a suitable positive number greater than 1 and F max is deduced from OCP1.

The endurance program
We consider an interval [0, t f ], where t f is the total training period formed by sequences of stimulation and rest periods. In this case, one must use a force-fatigue model and we take into account the constraint A ∈ [A rest , A rest /k ], where S = A rest /k corresponds to a fatigue threshold since, as reported in [7], if the user is exhausted, the force signal is totally noisy. Moreover in the case of exhaustion, a large rest period is required. This constraint can be penalized as follows OCP3: minσ where A S is related to S, while w 1 is a weight parameter.

Mathematical analysis of c N
A pulses train is defined by a finite sequence of impulse times σ = (t i ) 0≤i≤n such that t 0 < · · · < t n and we extend it to the left by t −1 = −∞ and to the right by t n+1 = T . The response c N can be decomposed as a sum of lobes defined as follows.
• defines a density probability function and more than 95% where Definition 10 The polynomial-exponential category for (piecewise) smooth functions [0, T ] → R is the category generated by sums, products of polynomials P (t) and exponential mappings to generate exponential-polynomials: n P n (t)e λnt . This category is stable with respect to derivation and integration.
Using proposition 4 one has: is in the polynomial-exponential category. Moreover, the coefficients are linear with respect to η i and polynomial-exponential with respect to t i .
We introduce the notion of p-persistent pulses related to the case where pth successive lobes have an influence on the (p + 1)th lobe.
Definition 12 Let t 1 < · · · < t p be p successive impulses times of a pulses train satisfying for any i ∈ 1, p , t i ≤ t i−1 + 5 and t i − t i−1 ≥ I min . When such integer p is maximal then the pulses train is said p-persistent.
(the factor R † j has a negligible effect).
We define now an approximation of c N denoted asc N that will be used to construct an approximation of the force F limiting the number of terms and the error between c N andc N is analyzed in the following proposition.
and define its (lower) approximation bỹ Then, we have: where · is the ceiling function and κ = min p, 5τc Imin is independent of k.
The number of indices i ∈ {0, . . . , k−2} for which t−t i ≤ 5τ c is at most κ := min p, 5τc Imin since (t i ) i is p-persistent and 5τc Imin stands for the maximum number of impulse times satisfying the constraint t i − t i−1 ≥ I min in an interval of length 5τ c .

Approximations of F
Integrating (7), the force with F (0) = 0 can be written as where M (t) = exp − t 0 m 2 (s) ds .
The following properties show that it is natural to approximate m 1 and m 2 by polynomial functions.

Remark 19
To constructF , the functions m 1 and m 2 were considered independently in the sense that the approximation does not rely on the relation (6) between m 1 and m 2 . A direct consequence is that an upper approximation of the force can be obtained from an upper approximation of m 1 . Outside the scope of this paper, this method may be applied for more general non-autonomous models.
Remark 21 A naive approach is to use classic integration schemes to define an explicit expression forF . Namely, using an explicit Euler scheme for the force equation (7) adapted to the partition (t i+j/p ) ij gives: where c ij = (1 − h i,j m 2 (t i+j/p )), h i,j = t i+(j+1)/p − t i+j/p and d ij = m 1 (t i+j/p ) for i = 0, . . . , n and j = 0, . . . , p − 1.
We deduce the following explicit expression forF (t kt+jt/p ), k t = 0, . . . , n, j t = 0, . . . p − 1: However, such method is not adapted for the design of our electrostimualtor (see Section 6). Indeed, it does not exploit the structure of the Hill functions m 1 and m 2 and yields worse results -in terms of time complexity and approximation error -compared to the approximation (18).

Error estimate.
We give a bound on the error between the approximationF and F in the case wherem 2|[t k ,t k+1 ] is the average of m 2 on [t k , t k+1 ].
(iii) There exists an unique We have: . As in the case (i), the first integral in the right hand side of (22) is bounded by (θ i − t i ) θi ti |ξ(u)| du. Likewise, the second integral in the right hand side of (22) is bounded by This concludes the proof.

Functional specification for the computation of a pulses train
The aim is to compute a local minimum σ * = (η * 0 , . . . , η * n , t * 1 , . . . , t * n , T ) ∈ R 2n+2 + of a cost function denoted as Θ. The free final time T adjusts automatically the optimal frequency of the pulses train. The functional specification of the electrostimulator imposes efficient computation of this minimum (real time computation) and this prevents us (at least when Θ involves the force) from using direct or indirect methods such as those presented in [2], mainly because these methods are based on a numerical scheme to approximate the variable F .

Finite dimensional optimization methods
We recall basic facts about finite dimensional optimization, see [5] for details, to emphasize that an optimal sampleddata control problem can be viewed as an instance of such optimization problem.
The optimization problems, associated to the optimal sampled-data control problems OCP1 and OCP2 presented in Section 3.2, can be written in the form: min where (σ) = (Ξ 1 (σ), . . . , Ξ 3n+5 (σ)) is the vector of constraints defined by: The cost Θ : σ → Θ(σ) related to the endurance or the force strengthening program is smooth with respect to σ. Consider the Lagrangian defined for all (σ, µ) ∈ R 2n+1 × R 3n+5 + by: The problem (23) is equivalent to the primal problem and the first order necessary optimality conditions for σ * to be a local minimizer, assuming the vectors Ξ i (σ * ), i ∈ {i, Ξ i (σ * ) = 0} to be linearly independent, state that there exists a Lagrange multiplier λ ∈ R 3n+5 such that We usually do not solve directly these optimality conditions to compute an optimal pair (σ * , λ * ), but a relaxation of these conditions can lead to efficient algorithm, namely the primal-dual interior point method [5].

Force optimization
We consider the problems the endurance and force strengthening optimization problems OCP1 and OCP2. For each problem, we give the approximationΘ of the cost functions Θ based on the approximationF of the variable F described in section 4. We solve the associated problem (23) -where Θ is replaced byΘ -using an interior point method on a standard computer 1 . Note thatΘ may consist of million of bytes, for that reason it is crucial to use an approximation of the gradient ofΘ with respect to t i , i = 1, . . . , n, computed via finite differences (vs formal computation). We initialize the pulses train to a regular partition of [0, 1] and the initial amplitudes being equal to 1.
Numerical result: The optimal solution σ * , the force response F and its approximationF are depicted in Fig.2.

Ca 2+ concentration optimization
Uniqueness of the optimal solution:  In Fig.4, we plot for n = 1 the objective function Θ(t 1 ) and for T fixed at some specific values.
The well-posedness of this optimization problem for any number n of impulsions times can be shown by inductive reasoning.

Cost: Θ(σ) = −c N (T ).
True cost. We have an explicit expression for c N , this problem can be easily solved numerically. We consider the finite dimensional optimization problem 23 where the cost is (Note that the amplitudes are fixed to 1).
Numerical result: Fig. 5 represents the time evolution of c N associated to (locally) optimal impulse times.

Cost:
In this case, we approximate Θ byΘ where the amplitudes and T are free andc N is given by Proposition 9.
Numerical result: The optimal solution σ * and its response c N are depicted in Fig. 6.
6 Isometric case: design of a smart muscular electrostimulator In this section we apply our study in the isometric case associated to the conception of a smart electrostimulator.
In this case, the task is to assign a constant reference force, but the general frame is suboptimal motion planning, see [11,12] for the theoretical foundations.
Advanced commercial muscular electrostimulator for training or reeducation purposes are based on the following. First of all, the user defines a program training. Basically endurance program with low frequency sequences of trains  (with constant interpulse) or force strengthening with high frequency trains. The program is a sequence of trains or rest periods. Before starting, the muscle is scanned to determine the parameters. A smart electrostimulator based on our study aims to design automatically such sequence, each program being translated into an optimization problem. Besides, in our framework, one can used VFT (Variable Frequency Trains) vs CFT (Constant Frequency Trains) in the standard case to complete the tuning of the amplitude.
One needs the following proposition associated to the endurance program presented in Fig.7 to illustrate the smart electrostimulator conception.

Conclusion
In this short article, we have mainly presented a finite dimensional approximation of the muscular force response to FES-input exploiting the mathematical structure of the model. The construction is based on the Ding et al. model but can be adapted to deal with the different models discussed in [19]. We have presented one application of our study related to motion planning in the isometric case in view to design a smart electrostimulator, which is an ongoing industrial project. Our approximation can be used to parameters estimation [19,17] and to design MPC-optimized sampled-data control schemes, applying standard algorithms [16,18] to this situation. Another application of our study is to track in the non-isometric case a path in the joint angle variable and this will be developed in a forthcoming article.  Fig. 7. Endurance program for a smart electrostimulator.