A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model

. The objective of this article is to make the analysis of the muscular force response to optimize electrical pulses trains using Ding et al. force-fatigue model. A geometric analysis of the dynamics is provided and very preliminary results are presented in the frame of optimal control using a simpliﬁed input- output model. In parallel, to take into account the physical constraints of the problem, partial state observation and input restrictions, an optimized pulses train is computed with a model predictive control, where a nonlinear observer is used to estimate the state-variables.

A case study of optimal input-output system with sampled-data control: Ding  1. Introduction. Functional Electrical Stimulation (FES) consists of applying an electrical stimulation to the muscle, in order to produce functional movements. It can be used for the muscular reinforcement, reeducation of the muscle and in the case of paralysis to activate the paralyzed muscles to produce movements. Mathematically FES leads to a sampled-data control problem which can be analyzed in this framework. The simulations of muscular response to electrical stimulations are based on dynamics models. The origin comes from the Hill-Langmuir equation in the context of biochemistry and pharmacology, see [12]. More recent models in the framework of model identification in non linear control are to due Bobet and Stein [2] and Law and Shields [15] and in this context a more sophisticated model was proposed by Ding et al. [6,7,8,9] where the force model is coupled to a fatigue model. This led to a set of five differential equations with a sampled-data control which can be used to stimulate and to control the force response to a pulses train of electrical stimulations and we shall refer to this model as the Ding et al. force-fatigue model in the sequel.
Only few works used this model to design an optimized train of pulses to control the force level [6,7,8,9]. Our aim being to make a more complete study of the problem in the framework of non linear optimal control, using sampled-data controls.
A geometric analysis of the model is provided and preliminary results in the frame of optimal control with sampled-data control [3]. It is based on a simplified dynamics using the geometric properties of the force-fatigue model and control reduction to simplify the physical control constraints. The complete system is analyzed in details using a control predictive strategy (MPC), see [17,19] coupled with a non linear observer technique based on [11] to estimate the state variable. This study completes [1] and gives non trivial optimized pulses trains which can be implemented in practise up to specific devices to speed on line numerical computations.
The article is organized as follows. In Section 2, we make a brief presentation of Ding et al. force-fatigue model based on [20]. In Section 3, the dynamics of the force model is briefly investigated to describe the input-output properties. In Section 4, the force-fatigue model is analyzed in the framework of geometric optimal sampleddata control systems and preliminary results are presented using a simplified model using a model reduction and an input transformation. Section 5 is devoted to the observer description. In Section 6, MPC method is presented using a further discretization of the dynamics to conclude by the algorithm used to compute in practise the optimized pulses trains. Numerical results are presented in the final Section 7.

2.
Mathematical force-fatigue model. We refer to [20] for a complete description and discussions of the model. The Ding et al. model studied in this article is presented next in the framework of model dynamics construction based on input-output observation and using the so-called tetania phenomenon in muscles responses. The first part of the model is the force (output) to pulses electrical stimulations (input). The pulses are normalized Dirac impulses at times 0 = t 1 < t 2 < . . . < t n : and I i = t i − t i−1 is the interpulse and convexifying leads to apply the input: for some amplitudes η i ∈ [0, 1], i ∈ {1, . . . , n}. Such pulses train feeds a first order model to produce an output according to the dynamicṡ where R i is a scaling factor associated to the phenomenon of tetania and corresponds to a accumulated effect of successive pulses and is modeled as OBSERVED WITH SAMPLED INPUT OF THE DING'S MODEL   3 where the magnitude is characterized by R 0 which is the limit case to high-frequency pulse. We denote E s (t) = u(t) the response to the pulses train and which is the (physical) control variable. It is of the form where The force response to such a train is modelled by the two equations of the socalled force model: which corresponds to a first order (resonant) linear dynamics which can be integrated with C N (0) = 0 and the force response F (t) is described by the equation where A is a fatigue parameter. Non-linear features of the model are described by the two mappings a, b: where K m , τ 1 are fatigue parameters and τ 2 an additional constant. The complete model is obtained by describing the evolutions of the parameters associated to fatigue and corresponds to the linear dynamics The full set of equations (3)-(4)-(6)-(7)-(8) are the force and fatigue model. We refer to Table 1 for the definitions and details of the symbols of the forcefatigue model.
For the purpose of our analysis the force fatigue model is written as the singleinput control system: corresponding to the sampled physical control data: (I 1 , I 2 , . . . , I n , η 1 , . . . , η n ) with constraints I min ≤ I i ≤ I max , 0 ≤ η i ≤ 1 corresponding to interpulses bounds and amplitude convexification. This leads to a non linear model with sampled control data with prescribed convex control constraints (η, I) ∈ C; η = (η 1 , . . . , η n ), I = (I 1 , . . . , I n ).
3. The force model. The force model can be briefly investigated to describe preliminary results.
Clearly we have Lemma 3.1. The above mapping is smooth with respect to I, η and piecewise smooth with respect to s.

3.2.
Smoothing process. For the sake of providing a smooth response for the observer it is sufficient to smooth the physical sampled control date as follows: use a bump function to smooth each Heaviside mapping H(t − t i ) at the sampling time t i .
3.3. Input-simplification. For the sake of the geometric analysis of the dynamics and to simplify the control constraints on (I, η) the FES signal E s (t) is taken as the input u(·) of the control system. Using (2) one can write Whence 0 = t 1 < t 2 < . . . < t n are given v(t) is a piecewise constant control depending upon the parameters η 1 , . . . , η n and the dynamics of the force model can be analyzed in the frame of geometric control. Fixing the interpulse led to a sample-data control model.

4.
Force-fatigue control model. First of all the control system (9) is analyzed in the framework of geometric control.
4.1. The concepts of geometric control system. Consider a (smooth) control single-input control system.
where x ∈ R n , y = h(x) ∈ R corresponds to a (smooth) single-observation mapping.
The following concepts rely on seminal results of geometric control, see [18,13]. We denote [U, V ] the Lie bracket of two (smooth) vector fields of R n : and a vector field U acts on (smooth) mappings f with the Lie derivative: Let D 1 = span{X, Y } and define recursively : Hence D k represents the Lie brackets of X, Y with lengths smaller than k + 1 and The observation space is the set of mappings: Taking x 0 ∈ R n , a frame at x 0 is a set of elements X 1 , . . . , X n of D L.A. such that: X 1 , . . . , X n are linearly independent at x 0 and acting by change of coordinates and (affine) feedback.
Fix x(0) = x 0 and T > 0 and consider the extremity mapping: The control u(·) is called singular on [0, T ] if the extremity mapping E x0,T is not of maximal rank n when evaluated at u(·).
Geometric analysis of an observation system of the form (14) amounts to compute D L.A. , the observation space, the singular controls and the feedback equivalence properties. Achievements of geometric optimal controls amounts to synthesize the optimal control in relation with the Lie brackets properties of a frame.

4.2.
Optimal control a force-fatigue model with sampled control data.

4.2.1.
Concepts. By the physical nature of the problem, the control in the force/fatigue model (1) is fixed at specific times and remains constant during some period of time. This falls into the framework of sampled-optimal control problem that we recall briefly and we refer to [3] for more details.
Consider a systemẋ = f (x, u). When the state x(·) and the control u(·) evolves continuously in time, we speak of a continuous-time optimal control problem and the control is said permanent. When the state x(·) evolves continuously in time whereas the control u(·) evolves in a discrete way, we speak of an optimal sampled-data control problem. Denoting by t f the final time, on [0, t f ] the control u takes its values in a discrete set. That means, for T > 0 fixed and k ∈ N, once the value of u(kT ) is chosen, then u(t) = u(kT ) for t ∈ [kT, (k + 1)T ]. T is called the sampling period and [kT, (k + 1)T ] is the sampling time interval where the control "is freezing".

4.2.2.
Application to the force-fatigue model. We present the geometric methods presented in [3] and numerical results on a reduced model in dimension 3 related to the force-fatigue model to validate the method.
Consider the minimization problem subject to the dynamic inx = (C N , F, A) defined by the equations and to the initial conditions The variables t f , K m , τ 1 , τ 2 , τ f at , α A , A rest are constant and fixed to some ad-hoc values.
The system (15) can be written into the forṁ Remark 1. The model (15) is a simplification of the force-fatigue model. In particular, without loss of generality, we don't take into account the factor exp(−t/τc) τc appearing in (2). Control constraints |u|≤ 1 are not the physical constraints, see Section 3.3. Also the fatigue dynamics is reduced to a single equation, motivated by the controllability properties of the system (6)- (7)- (8).
The permanent control case. The problem considered can be summarized as a permanent optimal control problem as follows The pseudo-Hamiltonian of the system is We denote by H i , i = 0, 1 the Hamiltonian lifts of the vector fieldsF i , i = 0, 1. Normal case: p 0 = −1/2. Applying the Pontryagin maximum principle (PMP), the optimal (permanent) control is given in the normal case by Abnormal case: p 0 = 0. Abnormal controls are characterized by the equation H 1 = 0, which differentiated twice by t leads to Also, we have Differentiating one more time led to: There are no admissible singular trajectories for the problem (15).
The sampled-data control case. The corresponding optimal sampled-data control problem is given by where T > 0 is a fixed sampling period such that t f = j T for some j ∈ N.
Normal case: p 0 = −1/2. Following [3], the optimal (sampled-data) control is for all controlling times kT, k ∈ {0, . . . , d}. WriteH then, the optimal sampled-data control is Numerical results. In Fig.1, we represent numerical results for several values of t f /T . The optimal permanent control is represented with thin continuous line. A Gaussian quadrature rule is used to computeH 1 using approximations of terms of the form p 1 (kT + T ), p 1 (kT + 1/2) . . . . We observe the convergence of the optimal sampled-data control to the permanent control as T tends to 0.
5. Observer. The model used for this study (five state variables) is based on force measurements collected from a set of subjects. Thus, the accuracy of the calculated parameters is directly related to these persons. In this study, we suppose that initial conditions are different following the subject under study. That is why, we need to estimate some of these five initial conditions. Indeed, from a first analysis we can deduce that the rest values C N (t 1 ) = 0 is well known and don't need estimation. The remaining parameters are A rest , K m,rest and τ 1,rest . Concerning K m,rest , a sensibility study is realized in order to determine its effect on the force variation.
5.1. Sensibility study of the force versus K m : The force evolution is compared for K m,rest and different values K m,rest (±30% of error) for I = 10ms, 50ms and 100ms (see Fig. 2 for I = 10ms). In the case of I = 10ms, the maximum force error is of −0.3%. Following Interpulse value, the maximum force error is obtained for I = 100ms (−1.3%) which means that a tolerance of ±30% gives force evolution with good accuracy.

5.2.
High-gain observer synthesis for the estimation of A and τ 1 . In this section, we design a modified version of the standard high-gain observer given in [11] taking into account the specific structure of the problem. The system defined by the force equation (4) and the fatigue model (6) and (8)   be rewritten as the single input-output system: , and m is a positive integer. Note that in (21), K m is not a state variable thanks to the weak sensibility of the solution with respect to this variable (see sensibility study of the force versus K m ). We introduce the change of variables φ: We have: where A sensibility study in (6) and (8) concerning respectively − A−Arest τ f at and − τ1−τ1,rest τ f at , shows that neglecting these two terms induces a maximum error of 7% for A and τ 1 , and 6% for the force value. Using the simplified model: with (4) gives: Considering the input E S (t) ≥ 0, we have In the second case of (27), we have Thus ( ∂φ ∂x (x(t)) = 0 for: • (C N , x 1 ) = 0 which correponds to the rest time (no control and then no need of state variable estimation).
In the simplified model (4),(24) and (25), there is no singularity of ∂φ ∂x (x(t)) −1 during stimulation period (C N , x 1 > 0), and ill-conditioning of the matrix is reduced. Based on the simplified model, the modified high-gain observer is defined as: Recall that m is a positive integer. It depends on the used pulses frequency for the stimulation of the muscle. S θ is a symmetric positive definite matrix given by the following Lyapunov equation: where θ is a tuning parameter, The terms of this matrix S θ = [S θ (l, k)] 1≤l,k≤3 have the following form:
• H1) ϕ n (u, z) is globally Lipschitz with respect to z and uniformly with respect to u.
• H2) ∃ U ∈ R of admissible controls, a compact K ∈ R n and two positive constants β min and β max such that: for all u ∈ U and all y(t) associated to u and initial conditions z(0) ∈ K, we have β min ≤ β(t) ≤ β max .
The observer has the following expression: where S θ is the solution of the Lyapunov equation: The solution of (34) can be rewritten as and with S 1 is the solution of (34) for θ = 1. Let ρ = β m and e =ẑ − z ⇒ė =ż −ż, thenė

OPTIMAL OBSERVED WITH SAMPLED INPUT OF THE DING'S MODEL 13
Consider the Lyapunov function: V =ē T S 1ē . Theṅ Then, we havė We deduce that Using (41) and assymption H1: For m and θ sufficiently large: γ > β m θ n K 1 ⇒V ≤ −γ 1 V, γ 1 > 0. In the particular case of the force fatigue model, β(t) is piecewise smooth, the lack of regularity is numerically bypassed by the choice of the integer m. For example, for I = 10ms, m = 3 is sufficient to estimate the whole variables. However, for I = 25ms, m must be at least equal to 5 (see observer simulations in Section 7). 6. Model predictive control (MPC). MPC computes a sequence of future controls to optimize "future" plant behavior (see Fig. 4). This computation is based on the dynamic model of the system, and must respect a certain cost and associated constraints (optimal control). In this framework, Model Predictive Heuristic Control (MPHC) was introduced by Richalet et al. [17] using impulse response type dynamic model. Dynamic Matrix Control (DMC) followed in 1980 (Cutler et al. [5]) using step response type dynamic model. State space formulation of MPC was introduced by Li et al. [16].
In the state space framework, the system takes the form: Consider the cost function to be minimized and Without constraints, Bellman's principle of optimality (1950) gives the solution. In the case of constrained problem, instead of J ∞ , we minimize the following cost: To explain the method, we consider the discrete linear system: where x(t) ∈ R n , u(t) ∈ R m , y(t) ∈ R p are respectively the state, input and output vectors, and t = kT s . This system (64) results from a discrete system modeling or continuous-discrete transformation. Practical considerations of implementation make this form more appropriate in the framework of MPC. We suppose that the system is both controllable and observable, and we call x(k + i/k) with i ≥ 0, the future state vector x(k + i) starting from t = kT s . x(k + i/k) is assumed to be available through using linear state variables estimator. In the case of unconstrained problem, a quadratic cost is minimized in order to obtain optimal control sequence. For instance, choose the following cost: The input-output relation along the receding horizon is then: Write (53) asȳ = Ψx(k/k) + Γū (54) whereȳ k ∈ R pNr ,ū k ∈ R mNr , Ψ ∈ R pNr×n , Γ ∈ R pNr×mNr .
Most frequently used constraints concern the minimum and maximum input, state and output variables: • The control variables: • The change rate of the control variables: (57) • Soft output variables constraints (relaxed output constraints using large slack variable s v to avoid constraints conflicts when solving control problem): • Soft state variables constraints (for the same raison as output constraints): In the case of constrained problem, the cost function (51) becomes: subject to (54) and a set of inequality constraints among (55) -(59).
This minimization problem could be expressed as a quadratic cost function subject to linear equality and inequality constraints: Clearly, optimization problem (61) is a convex Quadratic Programming (QP) problem; X * = [x * ,ū * ] T is than a global minimizer at each iteration.ū * being calculated, only u * (k/k) is applied at time t = kT s and optimization algorithm is repeated at t = (k+1)T s . In the case of NMPC (Nonlinear Model Predictive Controller), the dynamic is defined using a discretization of the general form (43). Using the same above cost function, optimization problem could be defined as follow: with X = [x,ū] T . Equality constraint in the optimization problem (62) is nonlinear. Convexity condition in (62) is not guaranteed, then X * could be a local minimizer. NMPC (Nonlinear Model Predictive controller) is used to solve this problem. Numerical solution are computed using Active Set, Primal-Dual or SQP methods (Wang [19], Fletcher [10] and Boyd and Vandenberghe [4]).
In the case of force or force-fatigue model, we consider the more general case of a varying T s (T s → I(i)), then u(i) = [I(i) η(i)] T , where I(i) = t(i) − t(i − 1) is the interpulse between two successive pulses and η(i) is the pulse amplitude applied at time t i . In this case,ū(k) = [u(k/k) . . . u((k + N r − 1)/k)] is the discrete receding horizon control vector (of length N r ) to be calculated by the NMPC at time t corresponding to iteration k. Then: u(k) = I(k/k) I(k + 1/k) · · · I((t(N r ))/k) η(t k /k) η((t k + I(k/k))/k) · · · η((t k + Nr−1 j=0 I(k + j))/k) (63) t(N r ) being the final time of the optimization horizon which is unknown a priori. For instance, using single move strategy u(k/k) = u(k + 1/k) = . . . = u((k + N r − 1)/k), (see Fig. 5).
The force-fatigue model being nonlinear, an NMPC is used to solve the problem with Nonlinear Programming Algorithm (NPA). In this case, criterion to be minimized is: subject to: 0 ≤ η(i) ≤ 1 and 0.01ms ≤ I(i) ≤ 0.1.
The estimation of the state variables vector (by the high-gain observer) is used as an initial variables vector to perform the NMPC over the horizon N r : ALGORITHM 1. Give F inal t , k = 1 Figure 5. (left half plane) E s and force profile for applied amplitude and interpulse stimulation, (right half plane) Predicted E s and force using single move strategy to be optimized.
7.1. High-gain observer. To exhibit the interest of the use of the power m in the nonlinear observer, let's consider MPC based nonlinear observer using only amplitude as control variable. We consider also the worst case of +30% of error of K m .
The following simulation results follow the stimulation protocol used in practise (a set of periods with stimulation and rest time slots in each period). During the stimulation time slot, the control is calculated to bring the force to F ref and in the rest time slot, the stimulation amplitude is set to 0. In this section, only two stimulation periods are considered.
In the cases of interpulse I = 10ms (100Hz) and I = 25ms (40Hz), the power values are m = 3 and m = 5, respectively. Fig.6 and Fig.7 represent the A estimates for I = 10ms and I = 25ms, respectively. Fig.8 and Fig.9 are the τ 1 estimates for I = 10ms and I = 25ms.Â converges after 50ms when I = 10ms and 100ms when I = 25ms. Concerningτ 1 , it converges after 75ms when I = 10ms and 200ms when I = 25ms. Large I seems to delay the convergence of the observer.  Figure 7. Evolution of A andÂ for I = 25, 30% error of K m . Fig.10 represents the force response for amplitude control strategy (for a receding horizon N r = 10) based on the proposed observer for I = 25ms and a force reference of 250N . Force mean value converges to the force reference after 200ms. 7.2. MPC with interpulse and amplitude as control variables. In this section, five stimualtion periods are presented (due to the experiment protocol). Figures 11, 14 and 15 are the force response in the case of F ref = 425N and (N r = 3, 5, 10), the interpulse and the amplitude controls for N r = 10, respectively. As expected, N r = 10 gives the best regulation performances (response time and overshoot). Before t = 6000ms with N r = 10, the force is correctly maintained at F ref = 425N . For t ≥ 6000ms (starting from the forth period), the fatigue level A is very high so that maximum value of amplitude control (see Fig. 15) and computed interpulse control (see Fig. 14) cannot maintain F at F ref . Maximum interpulse frequency is not used at the beginning of the forth and the fifth periods. Higher value of N r could correct this problem (supposing that global minimizer of MPC algorithm is reached at each iteration). However, increasing N r will render MPC algorithm time consuming which causes a problem for real time implementation.    of stimulations is dedicated to the control of the force level using both stimulation interpulse and amplitude controls. These simulations show the effect of the receding horizon on the control efficiency. Reasonable value of N r is however suitable to guarantee a short computation time for real time implementation.