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

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. 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 due to Bobet and Stein [2], Law and Shields [15] and a more sophisticated model was proposed by Ding et al. [6,7,8,9] where the force model is coupled to a fatigue model based on experimental results [20]. This led to a set of five differential equations with a sampled-data control which can be used to describe 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.

There is a limited literature that used this model to design an optimized train of pulses to control the force level [6,7,8,9,1]. Our aim is to produce 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 are presented in the framework 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 based on [11] to estimate the state variables.

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 investigated to describe the input-output properties. In Section 4, the force-fatigue model is analyzed in the framework of geometric optimal sampled-data control systems and preliminary results are presented with 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 with the algorithm implemented to compute in practice 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 the so-called tetania phenomenon in muscular responses. The first part of the model is the force response (output) to electrical stimulations pulses (input). The pulses are normalized Dirac impulses at times $ 0 = t_1<t_2<\ldots< t_n: $

and $ I_i = t_i-t_{i-1} $ is the interpulse and convexifying leads to apply the input:

with parameters $ \eta_i\in [0, 1], \; i\in \{1, \ldots, n\} $. Such pulses train feeds a first-order model to produce an output according to the dynamics

where $ R_i $ is a scaling factor associated to the phenomenon of tetania and corresponds to an accumulated effect of successive pulses and is modelled as

where the magnitude is characterized by $ R_0 $ which is the limit case to high-frequency pulse. The FES response to the input is denoted by $ E_s(t) $ and is given by

where $ H(t-t_i) = \left\{

$$\begin{array}{l} 0 \text{ if } t<t_i \\ 1 \text{ if } t\ge t_i \end{array}$$ \right. $ is the Heaviside function.

The force response to such a train is modelled by the two equations of the so-called 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 variable. Non-linear features of the model are described by the two mappings $ a, b $:

where $ K_m, \tau_1 $ are fatigue variables and $ \tau_2 $ an additional parameter.

The complete model is obtained by describing the evolutions of the variables 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 force-fatigue model.

For the purpose of our analysis the force-fatigue model is written as the single-input control system:

with $ x = (x_1, x_2, x_3, x_4, x_5) = (C_N, F, A, K_m, \tau_1) $ where $ u $ denotes the control $ u = E_s(t) $ corresponding to the sampled physical control data:

with constraints

corresponding to interpulses bounds and amplitude convexification. This leads to a non-linear model with sampled control data with prescribed convex control constraints $ (\eta, I)\in C;\; \eta = (\eta_1, \ldots, \eta_n), \; I = (I_1, \ldots, I_n) $.

3.   The force model

The force model can be briefly investigated to describe preliminary results.

3.1.   Parameterization

Observe that in the force model the static non-linearities described by the mappings $ a $ and $ b $ can be absorbed by time reparameterization to provide an explicit form for the force responses.

Proposition 1. The force model can be integrated by quadratures using a time reparameterization.

Proof. Integrating the linear system (3), the equation (4) can be written as

with $ \mathrm{d} s = b(t)\, \mathrm{d} t, \; c(s) = A(s)\, a(s)/b(s) $. It can be integrated using Lagrange formula.

Thus this gives an explicit force response in the time parameter $ s $

Clearly we have

Lemma 3.1. The above mapping is smooth with respect to $ I, \eta $ 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, \eta) $ the FES signal $ E_s(t) $ is taken as the input $ u(\cdot) $ of the control system. Using (2) one can write

with

Whence $ 0 = t_1<t_2<\ldots<t_n $ are given $ v(t) $ is a piecewise constant control depending upon the parameters $ \eta_1, \ldots, \eta_n $ and the dynamics of the force model can be analyzed in the frame of geometric control.

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 [13,14]

Consider a (smooth) control single-input control system.

where $ x\in \mathbb{R}^n, \, y = h(x)\in \mathbb{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 $ \mathbb{R}^n $:

and a vector field $ U $ acts on (smooth) mappings $ f $ with the Lie derivative:

Let $ D^1 = \mathrm{span}\{X, Y\} $ and define recursively : $ D^k = \mathrm{span} \{D^k\cup [D, D^{k-1}] \}, \; k>1. $ The length of $ [X_{i_1}, \ldots, [X_{i_{k-1}}, X_{i_k}]\ldots] $ is $ k $. Hence $ D^k $ represents Lie brackets of $ X, Y $ with lengths smaller than $ k+1 $ and denote $ D_{L.A.} = \mathrm{span} \cup_{k\ge 0} D^k $ the Lie algebra generated by $ \{X, Y\} $. The system is called weakly controllable if for each $ x\in \mathbb{R}^n, \; D_{L.A.}(x) = \mathbb{R}^n $.

The observation space is the set of mappings: $ \Theta = \{\mathcal{L}_{G}\, h;\; G\in D_{L.A.}\} $ and the system is called observable if for each $ x_1, x_2\in \mathbb{R}^n $, $ x_1\neq x_2 $ there exists $ \theta\in \Theta $ such that $ \theta(x_1)\neq \theta(x_2) $.

Taking $ x_0\in \mathbb{R}^n $, a frame at $ x_0 $ is a set of elements $ X_1, \ldots, X_n $ of $ D_{L.A.} $ such that: $ X_1, \ldots, X_n $ are linearly independent at $ x_0 $ and $ \sum_{i = 1}^n \text{length}(X_i) $ is minimal.

The system is called feedback linearizable in an open set $ \mathcal{V}\subset \mathbb{R}^n $ if $ \dot x = X(x)+u\, Y(x) $ is feedback equivalent to the linear system $ \dot x = Ax+ub $, that is there exists a diffeomorphism $ \varphi $ on $ \mathcal{V} $, $ \varphi(0) = 0 $ and a feedback $ u = \alpha(x)+\beta(x)\, v, \; \beta(x)\neq 0 $ such that $ g\cdot (X, Y) = (A, b) $ with $ g = (\varphi, \alpha, \beta) $ acting by change of coordinates and (affine) feedback.

Fix $ x(0) = x_0 $ and $ T>0 $ and consider the extremity mapping: $ E^{x_0, T}: u(\cdot)\in L^\infty([0, T]) \mapsto x(T, x_0, u) $ where $ x(\cdot) $ is the response of $ \dot x = X(x)+u\, Y(x) $ to $ u $ defined on $ [0, T] $.

The control $ u(\cdot) $ is called singular on $ [0, T] $ if the extremity mapping $ E^{x_0, T} $ is not of maximal rank $ n $ when evaluated at $ u(\cdot) $.

Geometric analysis of an observed 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 control 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

The control in the force-fatigue model (1) falls into the framework of sampled-optimal control problem and we refer to [3] for more details. We use the following terminology.

Consider a control system $ \dot x = f(x, u) $. When the state $ x(\cdot) $ and the control $ u(\cdot) $ evolve continuously in time, we speak of a continuous-time control problem and the control is said permanent. The sampled-data control case is when $ u(\cdot) $ is piecewise constant.

4.2.2.   Application to the force-fatigue model

For simplicity consider the reduced model in dimension $ 3 $ related to the force-fatigue model.

Let the minimization problem be defined by

subject to the dynamics in $ \tilde x = (C_N, F, A) $ given by

and with the initial conditions

The parameters $ T, \tau_2, \tau_{fat}, \alpha_A, A_{rest} $ and the variables $ K_m, \tau_1 $ are constant and fixed to some prescribed values.

The system (15) can be written into the form

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/\tau_c)/\tau_c $ appearing in (2). Control constraints $ |u|\le 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 can be summarized as a permanent optimal control problem as follows

The pseudo-Hamiltonian of the system is

where $ (p, p^0):[0, T]\mapsto \mathbb{R}^4 $ is the adjoint vector.

We denote by $ H_i = p\cdot \tilde F_i, \; i = 0, 1 $ the Hamiltonian lifts of the vector fields $ \tilde F_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 $ p\cdot \tilde F_1 = 0 $. Differentiating twice with respect to $ t $ leads to

Eliminating $ p $ we obtain $ D + u\, D' = 0 $ where

Computing we have

and $ D = 0 $ is equivalent to $ A (C_N (\tau_1+\tau_2)+K_m\tau_1)^2+F \tau_2 (C_N+K_m)^2 = 0. $ This implies $ F = A = 0 $ and the following lemma.

Lemma 4.1. There are no physically admissible singular trajectories for the problem (15).

The sampled-data control case. The corresponding optimal sampled-data control problem is given by

where $ T_s>0 $ is a fixed sampling period such that $ T = j\, T_s $ for some $ j\in \mathbb{N} $.

Normal case: $ p^0 = -1/2 $. Following [3], the optimal (sampled-data) control is

for all $ k\in \{0, \ldots, d\} $.

Write

then, the optimal sampled-data control is

Numerical results. In Fig. 1, we represent numerical results for several values of $ T/T_s $. The optimal permanent control is represented with thin continuous line. A Gaussian quadrature rule is used to compute $ \bar{H}_1 $. We observe the convergence of the optimal sampled-data control to the permanent control as $ T_s $ tends to $ 0 $.

5.   Observer

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}^{\prime} $ ($ \pm30\% $ of error) for $ I = 10ms $, $ 50ms $ and $ 100ms $ (see Fig. 2 for $ I = 10ms $). In the case of $ I = 10ms $, the maximum force error (see Fig. 3) is of $ -0.3\% $. Following interpulse value, the maximum force error is obtained for $ I = 100ms $ ($ -1.3\% $) which means that a tolerance of $ \pm30\% $ gives force evolution with a good accuracy.

5.2.   High-gain observer synthesis for the estimation of $ A $ and $ \tau_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) can be rewritten as the single input-output system:

with $ x = (F, A, \tau_1) \in \mathbb{R}^{3} $, $ y \in \mathbb{R} $, $ E_s \in \mathbb{R} $, $ 0<a<1 $ is given by (5) and $ m $ is a positive integer. Note that in (21), $ K_m $ is the solution of (7) thanks to the weak sensibility of the solution with respect to this variable (see Section 5.1 for the sensibility study of the force versus $ K_m $). We introduce the change of variables $ \phi $:

We have:

where

and

A sensibility study in (6) and (8) concerning respectively $ -(A-A_{rest})/\tau_{fat} $ and $ -(\tau_1-\tau_{1, rest})/\tau_{fat} $, shows that neglecting these two terms induces a maximum error of $ 7\% $ for $ A $ and $ \tau_1 $, and $ 6\% $ for the force value. Hence we use the simplified model:

Hence (4) gives:

Since $ E_s(t) \geq 0 $, we have:

The equation $ \det \left(\frac{\partial\phi}{\partial x}\right) = 0 $ yields

Since $ \alpha_{\tau_1} \sim 10^{-5} $, $ 0 \leq a \leq 1 $ and $ (x_3+\tau_2 a) \sim 10^2 $, we have from (28)

Thus $ (\frac{\partial\phi}{\partial x}(\hat{x}(t)) = 0, \; \hat x = (\hat F, \hat A, \hat \tau_1) $ for:

● $ (C_N, x_1) = 0 $ which corresponds to the period of rest,

● $ x_1 = 1/(2\alpha_{\tau_1}) $, $ x_1 \sim 10^{5} $ which is greater than the maximum force value.

In the simplified model (4), (24) and (25), $ \frac{\partial\phi}{\partial x}(\hat{x}(t)) $ is invertible during stimulation period $ (C_N, x_1>0) $, and the matrix is well-conditioned.

Based on the simplified model, the modified high-gain observer is defined as:

The experimental choice of $ m $ depends on the pulses frequency. $ S_{\theta} $ is a symmetric definite positive matrix given by the following Lyapunov equation:

where $ \theta $ is a high-gain tuning parameter introduced in [11], and

Hence the components of $ S_{\theta} = [S_{\theta}(l, k)]_{1 \leq l, k \leq 3} $ have the following form:

5.3.   High-gain observer convergence proof

Let

and

Then

Therefore

Take $ C = \left(

$$\begin{array}{cccc} \!1 \;\; 0 \;\; \cdots \;\; 0\! \end{array}$$ \right) $ and consider the system (33). We make the following assumptions.

5.3.1.   Assumptions

● H1) $ \varphi_n(u, z) $ is globally Lipschitz with respect to $ z $ and uniformly with respect to $ u $,

● H2) Let $ U $ be the control domain, a compact $ K\subset \mathbb{R}^n $ and two positive constants $ a_{\min} $ and $ a_{\max} $ such that: for all $ u $ valued in $ U $ and initial conditions $ z(0) $ $ \in $ $ K $, we have $ a_{\min}\leq a(t)\leq a_{\max} $,

● H3) $ a(\cdot) $ is $ C^1 $ (see Remark 2 below),

● H4) $ m\, \dot{a}(t)/ a(t)<1 $.

The observer has the following expression:

where $ S_{\theta} $ is the solution of the Lyapunov equation:

The solution of (36) can be rewritten as

and

with $ S_1 $ is the solution of (36) for $ \theta = 1 $. Let $ \rho = a^{m} $ and $ e = \hat{z}-z $ $ \Rightarrow $ $ \dot{e} = \dot{\hat{z}}-\dot{z} $, then

Write $ \bar{e} = \rho D_{\theta}e $, then

Consider the Lyapunov function: $ V = \bar{e}^TS_1\bar{e} $. Then

$ 2\bar{e}^T S_1 A \bar{e} = \bar{e}^T(S_1 A+A^T S_1)\bar{e} $ and $ S_1 A+A^T S_1 = -S_1+C^T C $.

Then, we have

We deduce that

Using (43) and assumption H1:

For $ m $ and $ \theta $ sufficiently large: $ \gamma>2\, a^m/\theta^n\, K_1 \Rightarrow \dot{V}\leq -\gamma_1\, V, \; \; \; \gamma_1>0 $.

Remark 2. In the force-fatigue model, $ a(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 predicted controls to optimize plant behavior (see Fig. 4).

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].

Consider a system written as

and a general cost function to be minimized

with

and

Without constraints, Bellman's principle of optimality (1950) gives the solution. In the case of constrained problem, instead of $ J_\infty $, we minimize the following cost:

with

and

This optimization over a finite horizon follows the algorithm:

1. Solve $ \min \; J(u(.), x_0, T) $ and find $ u^*(., x_0, T) $.

2. Apply $ u^*(., x_0, T) $ for $ \tau \in [0, T_s[ $, $ 0\leq T_s \leq T $, ($ T_s $: sampling period).

3. Repeat using $ x(T_s) $ instead of $ x_0 $.

6.1.   Discrete linear system (basic case)

To explain the method, we consider the discrete linear system:

where $ x(t) $ $ \in $ $ \mathbb{R}^n $, $ u(t) $ $ \in $ $ \mathbb{R}^m $, $ y(t) $ $ \in $ $ \mathbb{R}^p $ are respectively the state, input and output vectors, and $ t = kT_s $. This system (66) results from a discrete system modelling 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 denote by $ x(k+i/k), \, i\geq 0 $, the predicted state vector $ x(k+i) $ with initial condition $ x(kT_s) $.

Consider the unconstrained problem and the quadratic cost:

with $ Q = Q^T $ $ \succ $ $ 0 $, $ R = R^T $ $ \succeq $ $ 0 $ and $ S = S^T $ $ \succeq $ $ 0 $. The couple ($ y_{ref} $, $ u_{ref} $) is solution of (52) and $ \Delta u(k) = u(k)-u(k-1) $. The solution of this problem is obtained using the LQ controller.

In the case of constrained problem, LQR could not be used to solve control problem. System (52) allows to calculate $ x(k+i/k) $, $ i = 1, ..., N_r $ at time $ kT_s $, ($ N_r $ being the receding horizon length):

The input-output relation along the receding horizon is then:

Write (55) as

where $ \bar{y}_k $ $ \in $ $ \mathbb{R}^{p N_r} $, $ \bar{u}_k $ $ \in $ $ \mathbb{R}^{m N_r} $, $ \Psi $ $ \in $ $ \mathbb{R}^{p N_r \times n} $, $ \Gamma $ $ \in $ $ \mathbb{R}^{p N_r \times m N_r} $.

Constraints can take the form:

● Control variable:

● Change rate of the control variables:

with

● 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 reason as output constraints):

The cost function (53) becomes:

subject to (56) and a set of constraints among (57)-(61).

This minimization problem takes the form:

and $ \bar{y} = C\bar{x} $ with $ X = [\bar{x}, \bar{u}]^T $. Optimization problem (63) is a convex Quadratic Programming (QP) problem. Denote $ X^* = [\bar{x}^*, \bar{u}^*]^T $ the global minimizer at each iteration. Once $ \bar{u}^* $ is calculated, only $ u^*(k/k) $ is applied at time $ t = kT_s $ and we iterate the algorithm.

In the case of NMPC (non-linear Model Predictive Controller), the dynamics is defined using a discretization of the equation (45). Using (62), optimization problem becomes:

with $ X = [\bar{x}, \bar{u}]^T $. Equality constraint in the optimization problem (64) is non-linear. Now convexity condition in (64) is not guaranteed, then $ X^* $ could be a local minimizer. NMPC (non-linear Model Predictive controller) is used to solve this problem.

Numerical solutions are computed using Active Set, Primal-Dual or SQP methods (Wang [19], Fletcher [10] and Boyd and Vandenberghe [4]).

In the case of force-fatigue model, one must consider the more general case of a varying $ T_s $ ($ T_s $ $ \rightarrow $ $ I(i) $), then $ u(i) = [I(i) \; \eta(i)]^T $, where $ I(i) = t(i)-t(i-1) $ is the interpulse between two successive pulses and $ \eta(i) $ is the pulse amplitude applied at time $ t_i $. In this case, $ \bar{u}(k) = [u(k/k) \ldots 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:

$ t(N_r) $ being the final time of the optimization horizon which is unknown a priori. For instance, using single move strategy, we have $ u(k/k) = u(k+1/k) = \ldots = u((k+N_r-1)/k) $, (see Fig. 5).

The force-fatigue model being non-linear, an NMPC is used to solve the problem with Non-linear Programming Algorithm (NPA). In this case, criterion to be minimized is:

subject to: $ 0 \leq \eta(i) \leq 1 $ and $ 0.01 $ms$ \leq I(i) \leq 0.1 $ms.

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 $.

7.   Numerical results

7.1.   High-gain observer

To exhibit the interest of the use of the power $ m $ in the non-linear observer, consider a MPC based on a non-linear observer with only amplitude as control variable. We consider also the worst case of $ +30\% $ of error of $ K_m $.

The following simulation results follow the protocol used in practice (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.

For $ I = 10ms $ ($ 100Hz $) and $ I = 25ms $ ($ 40Hz $), we take $ m = 3 $ and $ m = 5 $ respectively. Fig. 6 and Fig. 7 represent the $ \hat A $ for $ I = 10ms $ and $ I = 25ms $, respectively. Fig. 8 and Fig. 9 are the $ \hat \tau_1 $ for $ I = 10ms $ and $ I = 25ms $. $ \hat{A} $ converges after $ 50ms $ when $ I = 10ms $ and $ 100ms $ when $ I = 25ms $. Concerning $ \hat{\tau_1} $, it converges after $ 75ms $ when $ I = 10ms $ and $ 200ms $ when $ I = 25ms $. Large $ I $ seems to delay the convergence of the observer.

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 stimulations periods are presented (due to the experimental protocol). Figures 11, 12, and 13, 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 \geq 6000ms $ (starting from the forth period), the fatigue level $ A $ is very high so that the maximum value of the amplitude control (see Fig. 13) and the interpulse control (see Fig. 12) 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 and cause a problem for real time implementation.

8.   Conclusion

This work deals with the control and estimation of the state variables of the Ding et al. force-fatigue model where the control is the interpulse and/or amplitude of the electrical stimulation. Preliminary geometric analysis of the force control controlling force level is presented with the aim of a future PMP control strategy. In the case of a fixed interpulse, the proposed high-gain observer using the force measurements exhibits the relation between the interpulse and the parameter $ m $ to perform accurate variables estimation. Model Predictive Control (MPC) strategy is presented in the case of both stimulation interpulse and amplitude as control variables. Simulation results are presented concerning the high-gain observer and the MPC force control. 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.