# Control and Estimation for the Design of a Smart Electrostimulator using Ding et al. Model

Toufik Bakir, Bernard Bonnard, Ilias Boualam, Mokrane Abdiche

► To cite this version:

Toufik Bakir, Bernard Bonnard, Ilias Boualam, Mokrane Abdiche. Control and Estimation for the Design of a Smart Electrostimulator using Ding et al. Model. FGS 2024 - French-German-Spanish Conference on Optimization, Jun 2024, Gijon, Spain. hal-04598660

HAL Id: hal-04598660

https://hal.science/hal-04598660

Submitted on 3 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Control and Estimation for the Design of a Smart Electrostimulator using Ding et al. Model

Toufik Bakir ¹, Bernard Bonnard ², Ilias Boualam ³, Mokrane Abdiche ⁴

1. toufik.bakir@u-bourgogne.fr ImVIA, université de Bourgogne, France

2. bernard.bonnard@u-bourgogne.fr Institut de Mathématique de Bourgogne, université de Bourgogne and Mc Tao team, Inria Sophia Antipolis, France

3. mohamedilias.boualam@segula.fr Entreprise Segula Matra Automotive, Trappes, France

4. mokrane.abdiche@segula.fr Entreprise Segula Matra Automotive, Trappes, France

## Abstract

Based on the A. V Hill's muscle model (Medicine Nobel prize 1922), mathematical models validated by experiments due to Ding et al. in the 2000's allow to describe the muscular force isometrical contraction due to electrostimulation, taking into account the fatigue. They serve as a model to control and to predict the effect of trains of electrical stimulations, with rest periods aiming to force rehabilitation or reinforcement. In this article we briefly present the main issues of the problem. Two typical training sessions are described related to increase the force or the endurance. Each program is translated into an optimization problem which is analyzed in the sample-data control frame. The parameters of the models split into parameters independent of each individual vs. parameters related mainly to the fatigue, which have to be online estimated. Geometric estimation theory leads to describe a software sensor to make explicit computations. NMPC algorithm vs MPC algorithm can be used to regulate the force.

## 1. Introduction

Recent mathematical models validate by experiments due mainly to Ding et al. [4–6] and based on the earliest work by A.G. Hill [9] allow to predict the isometric force response to external electrical stimulation, taking into account the fatigue phenomenon due to a long stimulation period. Such models contain two basic nonlinearities which constitute the intricate part of the dynamics. First of all, the ionic conduction and the nonlinear effect of successive pulses on the $Ca^{++}$ concentration. Second, the nonlinearity relating the muscular force response to such concentration, modeled by the Michaelis-Menten-Hill functions, which cause the force saturation called tetany. The control is formed by a sequence of trains of pulses which fit in the frame of sample-data control (digital controls) due to limitation on the interpulse. Our objective is to use the model to construct a smart electrostimulator for force rehabilitation or reinforcement based on two objectives : maximize the force response $F_{max}$ to a single train corresponding exactly to the tetany or an endurance session regulating the force to a reference force e.g. $\frac{F_{max}}{2}$ while minimizing the fatigue. Each training session is limited to 30 minutes since external stimulation causes severe fatigue and even during an endurance session rest periods have to be imposed. Besides those objectives they are computational limits related to on-board electronics and cost reduction. In particular integrating the nonlinear dynamics is time consuming and is bypassed by an approximation of the force response.

Hence the first part of this article is to briefly recall an off line formal approximation of the force dynamics to compute $F_{max}$. It is based on a piecewise linear approximations of the Michaelis-Menten-Hill functions and is fully described in [3]. The second part of this article is to describe an internal input-output model which is used to regulate the force-fatigue to a given level using Model Predictive Control [1, 13] and based on the Ding et al. (nonlinear)model. The final important issue of the project is to estimate the parameters based on preliminary experiments in the industrial realization of the electrostimulator. They are based mainly on the piezoelectric force sensor and the measurements being realized either at the beginning of each training session or during the "rest periods" where the muscle can be in reality stimulated with small intensity and frequency. Geometric control techniques developed in the 90s, see for instance [10], are used to study the observability of the system and to identify the "bad inputs" for which the systems is not observable. In the experiments they correspond roughly to the zero input where no force is produced. Recent geometric estimation techniques allow to identify the parameters of the model. The experiments they can be sorted into two types in the Ding et al. model: four parameters which are not depending on the individuals and four additional parameters related mainly to the fatigue phenomenon and which are depending of each individual and can be time varying. Such parameters can be estimated in the frame of geometric estimation developed in the 2000s and the general techniques presented in [7, 8, 11, 12] based on the construction of normal coordinates and Luenberger-type observers, where the effect of the inputs formed by trains of impulses on the $Ca^{++}$ concentration is described and leads to explicit estimation of the so-called observation space related to such inputs.

Numerical simulations are presented for the MPC algorithm based on the linear parametric model with online computation of the parameters and the train of pulses, taking into account the problem of severe fatigue caused by external stimulations.
# Control and Estimation for the Design of a Smart Electrostimulator using Ding et al. Model

Toufik Bakir ¹, Bernard Bonnard ², Ilias Boualam ³, Mokrane Abdiche ⁴

1*. toufik.bakir@u-bourgogne.fr ImVIA, université de Bourgogne, France

2. bernard.bonnard@u-bourgogne.fr Institut de Mathématique de Bourgogne, université de Bourgogne and Mc Tao team, Inria Sophia Antipolis, France

3. mohamedilias.boualam@segula.fr Entreprise Segula Matra, Trappes, France

4. mokrane.abdiche@segula.fr Entreprise Segula Matra, Trappes, France

## Abstract

Based on the A. V Hill's muscle model (Medicine Nobel prize 1922), mathematical models validated by experiments due to Ding et al. in the 2000's allow to describe the muscular force isometrical contraction due to electrostimulation, taking into account the fatigue. They serve as a model to control and to predict the effect of trains of electrical stimulations, with rest periods aiming to force rehabilitation or reinforcement. In this article we briefly present the main issues of the problem. Two typical training sessions are described related to increase the force or the endurance. Each program is translated into an optimization problem which is analyzed in the sample-data control frame. The parameters of the models split into parameters independent of each individual vs. parameters related mainly to the fatigue, which have to be online estimated. Geometric estimation theory leads to describe a software sensor to make explicit computations.

## 1. Introduction

Recent mathematical models validate by experiments due mainly to Ding et al. [4–6] and based on the earliest work by A.G. Hill [9] allow to predict the isometric force response to external electrical stimulation, taking into account the fatigue phenomenon due to a long stimulation period. Such models contain two basic nonlinearities which constitute the intricate part of the dynamics. First of all, the ionic conduction and the nonlinear effect of successive pulses on the $Ca^{++}$ concentration. Second, the nonlinearity relating the muscular force response to such concentration, modeled by the Michaelis-Menten-Hill functions, which cause the force saturation called tetany. The control is formed by a sequence of trains of pulses which fit in the frame of sample-data control (digital controls) due to limitation on the interpulse. Our objective is to use the model to construct a smart electrostimulator for force rehabilitation or reinforcement based on two objectives : maximize the force response $F_{max}$ to a single train corresponding exactly to the tetany or an endurance session regulating the force to a reference force e.g. $\frac{F_{max}}{2}$ while minimizing the fatigue. Each training session is limited to 30 minutes since external stimulation causes severe fatigue and even during an endurance session rest periods have to be imposed. Besides those objectives they are computational limits related to on-board electronics and cost reduction of the commercial device. In particular integrating the nonlinear dynamics is time consuming and is bypassed by an approximation of the force response.

Hence the first part of this article is to briefly recall an off line formal approximation of the force dynamics to compute $F_{max}$. It is based on a piecewise linear approximations of the Michaelis-Menten-Hill functions and is fully described in [3]. The second part of this article is to describe an internal input-output model which is used to regulate the force-fatigue to a given level using Model Predictive Control [1, 13] and based on the Ding et al. (nonlinear)model. The final important issue of the project is to estimate the parameters based on preliminary experiments in the industrial realization of the electrostimulator. They are based mainly on the piezoelectric force sensor and the measurements being realized either at the beginning of each training session or during the “rest periods” where the muscle can be in reality stimulated with small intensity and frequency. Geometric control techniques developed in the 90s, see for instance [10], are used to study the observability of the system and to identify the “bad inputs” for which the systems is not observable. In the experiments they correspond roughly to the zero input where no force is produced. Recent geometric estimation techniques allow to identify the parameters of the model. Based on the experiments they can be sorted into two types in the Ding et al. model: four parameters which are not depending on the individuals and four additional parameters related mainly to the fatigue phenomenon and which are depending of each individual and can be time varying. Such parameters can be estimated in the frame of geometric estimation developed in the 2000s and the general techniques presented in [7, 8, 11, 12] based on the construction of normal coordinates and Kalman or Luenberger-type observers, and adaptive control where the effect of the inputs formed by trains of impulses on the $Ca^{++}$ concentration is described and leads to explicit estimation of the so-called observation space related to such inputs.

Numerical simulations are presented for the MPC algorithm based on the linear parametric model with online computation of the parameters and the train of pulses in the industrial prototype, taking into account the problem of severe fatigue caused by external stimulations.
## 2. Force-fatigue model

### 2.1. Ding et al. force-fatigue model

The FES input over a pulse train [0, T] is modelled as a sum of Dirac pulses by

$$t \to \sum_{i=0,n} \eta_i \delta(t - t_i), \qquad (2.1)$$

<figure><img src="image_1.png" /><figcaption>Fig. 1 stimulation period, stimulation and rest sub-periods</figcaption></figure>

where $0 = t_0 < t_1... < t_n < t_{n+1} = T$ are the impulses times with $n \in \mathbb{N}$ being fixed (see figure 1 for constant interpulse and stimulation amplitude), $\eta_i$ being the amplitude of each pulse, which are convexified by taking $\eta_i \in [0, 1]$ and $\delta(\cdot - t_i)$ denoting the Dirac at time $t_i$. We denote by $I_i = t_i - t_{i-1}$ the interpulse and we have a digital constraint $I_i \ge I_m$ in the problem e.g. $I_m \ge 30ms$ for a train of 10 impulses of around $T = 500ms$. Such a control provides the FES signal taken as the physical input, using a linear filter (first-order linear dynamics).

$$\frac{dE}{dt}(t) + \frac{1}{\tau_c}E(t) = \sum R_i \eta_i \delta(t - t_i) \qquad (2.2)$$

so that it takes the form

$$E(t) = \frac{1}{\tau_c} \sum_{i=0,n} R_i e^{-\frac{t-t_i}{\tau_c}} \eta_i H(t - t_i), \qquad (2.3)$$

where $H$ is the Heaviside function. $E(t)$ depends upon the time response parameter $\tau_c$ and the scaling function $R_i$ depending on parameter $R(0)$ as following:

$$R_0 = 1, R_i = 1 + (R(0) - 1)e^{-(t_i - t_{i-1})/\tau_c}, i = 1, \dots, n, \qquad (2.4)$$

which codes the memory effect of successive muscle contractions.

The FES signal drives the evolution of the electrical conduction describing the evolution of Ca$^{++}$-concentration $c_N$ which is related to the force response $F$. The dynamics being described by

$$\frac{dc_N}{dt}(t) = E(t) - \frac{c_N(t)}{\tau_c}, \qquad (2.5)$$

$$\frac{dF}{dt}(t) = -m_2(t)F(t) + m_1(t)A(t) \qquad (2.6)$$

where

$$m_1(t) = \frac{c_N(t)}{K_m + c_N(t)}, m_2(t) = \frac{1}{\tau_1 + \tau_2 m_1(t)}. \qquad (2.7)$$

Hence six parameters are introduced in the model ($\tau_c$, $R(0)$, $\tau_1$, $\tau_2$, $K_m$, $A(t)$), where to simplify ($\tau_c$, $R(0)$, $\tau_1$, $\tau_2$, $K_m$) are fixed parameters and the time variable parameter $A(t)$ is the scaling force parameter which is used to model the fatigue dynamics according to

$$\frac{dA}{dt}(t) = -\frac{A(t) - A_0}{\tau_{fat}} + \alpha_A F(t). \qquad (2.8)$$
## 2. Force-fatigue model

### 2.1. Ding et al. force-fatigue model

The FES input over a pulse train [0, T] is modelled as a sum of Dirac pulses by

$$t \rightarrow \sum_{i=0,n} \eta_i \delta(t - t_i), \qquad (2.1)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1 stimulation period, stimulation and rest sub-periods</figcaption></figure>

where $0 = t_0 < t_1... < t_n < t_{n+1} = T$ are the impulses times with $n \in \mathbb{N}$ being fixed (see figure 1 for constant interpulse and stimulation amplitude), $\eta_i$ being the amplitude of each pulse, which are convexified by taking $\eta_i \in [0, 1]$ and $\delta(\cdot - t_i)$ denoting the Dirac at time $t_i$. We denote by $I_i = t_i - t_{i-1}$ the interpulse and we have a digital constraint $I_i \ge I_m$ in the problem e.g. $I_m \ge 30ms$ for a train of 10 impulses of around $T = 500ms$. Such a control provides the FES signal taken as the physical input, using a linear filter (first-order linear dynamics).

$$\frac{dE}{dt}(t) + \frac{1}{\tau_c}E(t) = \sum R_i \eta_i \delta(t - t_i) \qquad (2.2)$$

so that it takes the form

$$E(t) = \frac{1}{\tau_c} \sum_{i=0,n} R_i e^{-\frac{t-t_i}{\tau_c}} \eta_i H(t - t_i), \qquad (2.3)$$

where $H$ is the Heaviside function. $E(t)$ depends upon the time response parameter $\tau_c$ and the scaling function $R_i$ depending on parameter $R(0)$ as following:

$$R_0 = 1, R_i = 1 + (R(0) - 1)e^{-(t_i - t_{i-1})/\tau_c}, i = 1, \dots, n, \qquad (2.4)$$

which codes the memory effect of successive muscle contractions.

The FES signal drives the evolution of the electrical conduction describing the evolution of Ca$^{++}$-concentration $c_N$ which is related to the force response $F$. The dynamics being described by

$$\frac{dc_N}{dt}(t) = E(t) - \frac{c_N(t)}{\tau_c}, \qquad (2.5)$$

$$\frac{dF}{dt}(t) = -m_2(t)F(t) + m_1(t)A(t) \qquad (2.6)$$

where

$$m_1(t) = \frac{c_N(t)}{K_m + c_N(t)}, m_2(t) = \frac{1}{\tau_1 + \tau_2 m_1(t)}. \qquad (2.7)$$

Hence six parameters are introduced in the model ($\tau_c$, $R(0)$, $\tau_1$, $\tau_2$, $K_m$, $A(t)$), where to simplify ($\tau_c$, $R(0)$, $\tau_1$, $\tau_2$, $K_m$) are fixed parameters and the time variable parameter $A(t)$ is the scaling force parameter which is used to model the fatigue dynamics according to

$$\frac{dA}{dt}(t) = -\frac{A(t) - A_0}{\tau_{fat}} + \alpha_A F(t). \qquad (2.8)$$
<table><caption>**Tab. 1** Ding et al. model parameters</caption><thead><tr><th>Symbol</th><th>Unit</th><th>Value</th><th>description</th></tr></thead><tbody><tr><td>c<sub>N</sub></td><td>—</td><td>—</td><td>Normalized amount of Ca<sup>2+</sup>-troponin complex</td></tr><tr><td>F</td><td>N</td><td>—</td><td>Force generated by muscle</td></tr><tr><td>t<sub>i</sub></td><td>ms</td><td>—</td><td>Time of the i<sup>th</sup> pulse</td></tr><tr><td>n</td><td>—</td><td>—</td><td>Total number of the pulses before time t</td></tr><tr><td>i</td><td>—</td><td>—</td><td>Stimulation pulse index</td></tr><tr><td>τ<sub>c</sub></td><td>ms</td><td>20</td><td>Time constant that commands the rise and the decay of C<sub>N</sub></td></tr><tr><td>R(0)</td><td>—</td><td>1.143</td><td>Term of the enhancement in C<sub>N</sub> from successive stimuli</td></tr><tr><td>A</td><td>N/ms</td><td>—</td><td>Scaling factor for the force and the shortening velocity of muscle</td></tr><tr><td>τ<sub>1</sub></td><td>ms</td><td>—</td><td>Force decline time constant when strongly bound cross-bridges absent</td></tr><tr><td>τ<sub>2</sub></td><td>ms</td><td>124.4</td><td>Force decline time constant due to friction between actin and myosin</td></tr><tr><td>K<sub>m</sub></td><td>—</td><td>—</td><td>Sensitivity of strongly bound cross-bridges to C<sub>N</sub></td></tr><tr><td>A<sub>rest</sub></td><td>N/ms</td><td>3.009</td><td>Value of the parameter A when muscle is not fatigued</td></tr><tr><td>K<sub>m,rest</sub></td><td>—</td><td>0.103</td><td>Value of the parameter K<sub>m</sub> when muscle is not fatigued</td></tr><tr><td>τ<sub>1,rest</sub></td><td>ms</td><td>50.95</td><td>The value of the parameter τ<sub>1</sub> when muscle is not fatigued</td></tr><tr><td>α<sub>A</sub></td><td>1/m<sub>s</sub><sup>2</sup></td><td>-4.0 10<sup>-7</sup></td><td>Coefficient for the force-model parameter A in the fatigue model</td></tr><tr><td>α<sub>K<sub>m</sub></sub></td><td>1/m<sub>s</sub>N</td><td>1.9 10<sup>-8</sup></td><td>Coefficient for the force-model parameter K<sub>m</sub> in the fatigue model</td></tr><tr><td>α<sub>τ<sub>1</sub></sub></td><td>1/N</td><td>2.1 10<sup>-5</sup></td><td>Coefficient for force-model parameter τ<sub>1</sub> in the fatigue model</td></tr><tr><td>τ<sub>fat</sub></td><td>s</td><td>127</td><td>Time constant controlling the recovery of (A, K<sub>m</sub>, τ<sub>1</sub>)</td></tr></tbody></table>

Values of the parameters are reported in the reference [4] (see table 1) in the frame of Ding et al. experiments, the system formed by (2.5) and (2.6) describing the non-fatigue model, while the additional equation (2.8) is describing the fatigue and depends on two parameters $α_A$ which defines the "slope" of the fatigue evolution while $τ_{fat}$ is the time constant controlling the recovery to the rest point $A_{rest} = A(0)$. The model provides a closed curve $t \to A(t)$ obtained from the fatigue dynamics associated to concatenation of two arcs: the first one associated to the application of the averaged force $F_{averaged} = \frac{1}{T} \int_{0}^{T} F(t)dt$ over a pulse train on $[0, T]$ and the recovery arc during the complete rest period where no force is applied.

The main properties of the dynamics of the non-fatigue model is resumed in two lemmas.

**Lemma 2.1** For a pulse train defined by $\sigma = (t_0 = 0, t_1, ..., t_n, t_{n+1} = T, \eta_0, \eta_1, ..., \eta_n)$ the concentration $c_N$ can be written as the superposition of $n + 1$ lobes

$$c_N(t) = \frac{1}{\tau_c} \sum_{i=0,n} R_i \eta_i(t - t_i) e^{-\frac{t-t_i}{\tau_c}} H(t - t_i) \quad (2.9)$$

which represents a piecewise polynomial-exponential mapping.

**Lemma 2.2** The force dynamics in the non-fatigue case can be written as

$$\frac{dF}{ds}(s) = c(s) - F(s),$$

using the time reparameterization $ds = m_2(t)dt$ and can be integrated by quadrature using Lagrange formula. This gives an explicit force response $s \to F(s)$ which is smooth with respect to the control parameters and $s$ at each time different of a pulse time $t_i$.

From which we deduce the following, see [3] for the complete details and numerical simulations.

# 3. Construction of the approximation of the force response to a single train and the Punch Program in non-fatigue case

## 3.1. Approximation

According to (2.9) each lobe $l_k$ is given by

$$l_k = R_k \eta_k \frac{t - t_k}{\tau_c} e^{-(t - t_k)/\tau_c} H(t - t_k),$$
<table><caption>**Tab. 1** Ding et al. model parameters</caption><thead><tr><th>Symbol</th><th>Unit</th><th>Value</th><th>description</th></tr></thead><tbody><tr><td>C<sub>N</sub></td><td>—</td><td>—</td><td>Normalized amount of Ca<sup>2+</sup>-troponin complex</td></tr><tr><td>F</td><td>N</td><td>—</td><td>Force generated by muscle</td></tr><tr><td>t<sub>i</sub></td><td>ms</td><td>—</td><td>Time of the i<sup>th</sup> pulse</td></tr><tr><td>n</td><td>—</td><td>—</td><td>Total number of the pulses before time t</td></tr><tr><td>i</td><td>—</td><td>—</td><td>Stimulation pulse index</td></tr><tr><td>τ<sub>c</sub></td><td>ms</td><td>20</td><td>Time constant that commands the rise and the decay of C<sub>N</sub></td></tr><tr><td>R(0)</td><td>—</td><td>1.143</td><td>Term of the enhancement in C<sub>N</sub> from successive stimuli</td></tr><tr><td>A</td><td>N/ms</td><td>—</td><td>Scaling factor for the force and the shortening velocity of muscle</td></tr><tr><td>τ<sub>1</sub></td><td>ms</td><td>—</td><td>Force decline time constant when strongly bound cross-bridges absent</td></tr><tr><td>τ<sub>2</sub></td><td>ms</td><td>124.4</td><td>Force decline time constant due to friction between actin and myosin</td></tr><tr><td>K<sub>m</sub></td><td>—</td><td>—</td><td>Sensitivity of strongly bound cross-bridges to C<sub>N</sub></td></tr><tr><td>A<sub>rest</sub></td><td>N/ms</td><td>3.009</td><td>Value of the parameter A when muscle is not fatigued</td></tr><tr><td>K<sub>m,rest</sub></td><td>—</td><td>0.103</td><td>Value of the parameter K<sub>m</sub> when muscle is not fatigued</td></tr><tr><td>τ<sub>1,rest</sub></td><td>ms</td><td>50.95</td><td>The value of the parameter τ<sub>1</sub> when muscle is not fatigued</td></tr><tr><td>α<sub>A</sub></td><td>1/m<sub>s</sub><sup>2</sup></td><td>-4.0 10<sup>-7</sup></td><td>Coefficient for the force-model parameter A in the fatigue model</td></tr><tr><td>α<sub>K<sub>m</sub></sub></td><td>1/m<sub>s</sub>N</td><td>1.9 10<sup>-8</sup></td><td>Coefficient for the force-model parameter K<sub>m</sub> in the fatigue model</td></tr><tr><td>α<sub>τ<sub>1</sub></sub></td><td>1/N</td><td>2.1 10<sup>-5</sup></td><td>Coefficient for force-model parameter τ<sub>1</sub> in the fatigue model</td></tr><tr><td>τ<sub>fat</sub></td><td>s</td><td>127</td><td>Time constant controlling the recovery of (A, K<sub>m</sub>, τ<sub>1</sub>)</td></tr></tbody></table>

Values of the parameters are reported in the reference [4] (see table 1) in the frame of Ding et al. experiments, the system formed by (2.5) and (2.6) describing the non-fatigue model, while the additional equation (2.8) is describing the fatigue and depends on two parameters $α_A$ which defines the "slope" of the fatigue evolution while $τ_{fat}$ is the time constant controlling the recovery to the rest point $A_{rest} = A(0)$. The model provides a closed curve $t \to A(t)$ obtained from the fatigue dynamics associated to concatenation of two arcs: the first one associated to the application of the averaged force $F_{averaged} = \frac{1}{T} \int_{0}^{T} F(t)dt$ over a pulse train on $[0, T]$ and the recovery arc during the complete rest period where no force is applied.

The main properties of the dynamics of the non-fatigue model is resumed in two lemmas.

**Lemma 2.1** For a pulse train defined by $\sigma = (t_0 = 0, t_1, ..., t_n, t_{n+1} = T, \eta_0, \eta_1, ..., \eta_n)$ the concentration $C_N$ can be written as the superposition of $n + 1$ lobes

$$c_N(t) = \frac{1}{\tau_c} \sum_{i=0,n} R_i \eta_i (t - t_i) e^{-\frac{t-t_i}{\tau_c}} H(t - t_i) \quad (2.9)$$

which represents a piecewise polynomial-exponential mapping.

**Lemma 2.2** The force dynamics in the non-fatigue case can be written as

$$\frac{dF}{ds}(s) = c(s) - F(s),$$

using the time reparameterization $ds = m_2(t)dt$ and can be integrated by quadrature using Lagrange formula. This gives an explicit force response $s \to F(s)$ which is smooth with respect to the control parameters and $s$ at each time different of a pulse time $t_i$.

From which we deduce the following, see [3] for the complete details and numerical simulations.

# 3. Construction of the approximation of the force response to a single train and the Punch Program in non-fatigue case

## 3.1. Approximation

According to (2.9) each lobe $l_k$ is given by

$$l_k = R_k \eta_k \frac{t - t_k}{\tau_c} e^{-(t - t_k)/\tau_c} H(t - t_k),$$
the lobe reaches its maximum at $t = t_k + \tau_c$ which is equal to $R_k\eta_k/e$ and is concave on $[t_k, t_k + 2\tau_c]$ and can be approximated by its restriction to $[t_k, t_k + 5\tau_c]$. The restriction of $m_1$ to one lobe is maximal when the concentration $c_N$ is maximal and we denote $t_k^*$ the corresponding time. Let $\sigma$ be the sequence defined in (2.9) and assume that the minimal interpulse is such that $I_m \ge \tau_c$.

We divide the subdivision $t_0 < t_1 < \dots < t_n < T$ introducing the intermediate times $t_k^*$ where $m_1$ and $m_2$ are respectively approximated by a piecewise linear mapping and a piecewise constant mapping on each subinterval.

This leads to an explicit formula for the force response $F$ on $[0, T]$. Note that this basic partition can be refined to improve the approximation, see [3] for the complete description. The approximation contains the parameters of Ding et al. model.

## 3.2. Punch program

Using the previous force approximation denoted $F_{approximation}$ one compute a local minimum $\sigma^*$ over the set of pulse trains $\sigma$. The details of the optimization algorithm and the numerical simulations are presented in [3].

# 4. Endurance session using the force-fatigue model and the MPC algorithm

## 4.1. Notation 1

The force fatigue model described by (2.5), (2.6), (2.8) is shortly written as

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)), \qquad (4.1) $$

where $x = (c_n, F, A)^T$ is the state variable and $u(t)$ denotes the general input corresponding to the FES signal. Restriction on $u$ are imposed by the physical device: bounds implied by the constraints $\eta_i \in [0, 1]$, sampling times and interpulse constraints. They will be considered as soft constraints in the MPC algorithm and relaxed in the control computations using a quadratic optimization method.

We assume that in the dynamics two variables are observed $y = h(x) = (F, A)$ which defines the construction of the input-output dynamics.

In the MPC algorithm we consider the following discrete linear input-output system

$$ \left\{ \begin{array}{l} VX(k + 1) = MA_kVX(k) + MB_kU(k) \\ y(k) = CX(k) \end{array} \right. \qquad (4.2) $$

where:

$$ MA_k = \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix}_k, MB_k = \begin{pmatrix} b_1 \\ b_2 \end{pmatrix}_k, VX(k) = \begin{pmatrix} F_{m_k} \\ A_{m_k} \end{pmatrix}. \qquad (4.3) $$

The parametric model (4.2) results from the identification routine minimizing the criterion:

$$ J = \min_{N_i=N_{i_1}, \dots, N_{i_{max}}} \frac{1}{N_i} \sum_{j=k-N_i}^{k} \left( \left( \begin{pmatrix} F_{mean_j} \\ A_{mean_j} \end{pmatrix} - \begin{pmatrix} F_{m_j} \\ A_{m_j} \end{pmatrix} \right)^2 \right). \qquad (4.4) $$

<figure><img src="image_2.png" /><figcaption>Fig. 2 Force, mean force and identification horizons to identify the parametric model</figcaption></figure>
the lobe reaches its maximum at $t = t_k + \tau_c$ which is equal to $R_k \eta_k / e$ and is concave on $[t_k, t_k + 2\tau_c]$ and can be approximated by its restriction to $[t_k, t_k + 5\tau_c]$. The restriction of $m_1$ to one lobe is maximal when the concentration $c_N$ is maximal and we denote $t_k^*$ the corresponding time. Let $\sigma$ be the sequence defined in (2.9) and assume that the minimal interpulse is such that $I_m \ge \tau_c$.

We divide the subdivision $t_0 < t_1 < ... < t_n < T$ introducing the intermediate times $t_k^*$ where $m_1$ and $m_2$ are respectively approximated by a piecewise linear mapping and a piecewise constant mapping on each subinterval.

This leads to an explicit formula for the force response $F$ on $[0, T]$. Note that this basic partition can be refined to improve the approximation, see [3] for the complete description. The approximation contains the parameters of Ding et al. model.

## 3.2. Punch program

Using the previous force approximation denoted $F_{approximation}$ one compute a local minimum $\sigma^*$ over the set of pulse trains $\sigma$. The details of the optimization algorithm and the numerical simulations are presented in [3].

# 4. Endurance session using the force-fatigue model and the MPC algorithm

## 4.1. Notation 1

The force fatigue model described by (2.5), (2.6), (2.8) is shortly written as

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)), \qquad (4.1) $$

where $x = (c_n, F, A)^T$ is the state variable and $u(t)$ denotes the general input corresponding to the FES signal. Restriction on $u$ are imposed by the physical device: bounds implied by the constraints $\eta_i \in [0, 1]$, sampling times and interpulse constraints. They will be considered as soft constraints in the MPC algorithm and relaxed in the control computations using a quadratic optimization method.

We assume that in the dynamics two variables are observed $y = h(x) = (F, A)$ which defines the construction of the input-output dynamics.

In the MPC algorithm we consider the following discrete linear input-output system

$$ \begin{cases} VX(k+1) = MA_k VX(k) + MB_k U(k) \\ y(k) = CX(k) \end{cases} \qquad (4.2) $$

where:

$$ MA_k = \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix}_k, MB_k = \begin{pmatrix} b_1 \\ b_2 \end{pmatrix}_k, VX(k) = \begin{pmatrix} F_{m_k} \\ A_{m_k} \end{pmatrix}. \qquad (4.3) $$

The parametric model (4.2) results from the identification routine minimizing the criterion:

$$ J = \min_{N_i=N_{i_1},...,N_{i_{max}}} \frac{1}{N_i} \sum_{j=k-N_i}^{k} \left( \left( \begin{pmatrix} F_{mean_j} \\ A_{mean_j} \end{pmatrix} - \begin{pmatrix} F_{m_j} \\ A_{m_j} \end{pmatrix} \right)^2 \right). \qquad (4.4) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2 Force, mean force and identification horizons to identify the parametric model</figcaption></figure>
$N_i$ being the backward identification horizon. Figure 2 represents the force, the mean force and the backward identification horizon to be found in order to get the best parametric model. The same figure could be constructed for A. $F_{mean_k}$ and $A_{mean_k}$ are calculated as following:

$$ 
\begin{aligned}
F_{mean_k} &= \frac{1}{t_{k+1}-t_k} \int_{t_k}^{t_{k+1}} F(\xi) d\xi \\
A_{mean_k} &= \frac{1}{t_{k+1}-t_k} \int_{t_k}^{t_{k+1}} A(\xi) d\xi
\end{aligned}
\qquad (4.5) $$

The criterion (4.4) traduces the fact that MA and MB are updated at each iteration ($MA_k, MB_k$). The couple is used in MPC strategy to calculate the control value $\eta_k$. The frequency of the stimulation being fixed.

## 4.2. MPC Algorithm

We present a version of the algorithm to illustrate the procedure which is classical, and the quadratic cost can be modified. The control constraints have been relaxed and they have to be introduced later to define the true feedback control.

We fix a sequence $k = 1, ..., K$ where $N_p$ is the prediction horizon and an output reference trajectory $y_{ref}$ associated to regulation of the force response $F$ to a fixed level $F_{max}/\rho$ with $\rho > 1$ and a fatigue reference $A_{ref}(.)$. Denoting by $e(k) = (y(k) - y_{ref}(k))$, we minimize a cost of the form

$$ J(y, u) = \sum_{k=1,K} \lambda_1 \| e(k) \|_2^2 + \lambda_2 \| \Delta u_k \|_2^2, \qquad (4.6) $$

where $\Delta u(k)$ is the control increment and $\lambda_i$ are weighting parameters.

The feedback control $u(k)$ is computed on the horizon $K$ solving the LQ-problem defined by the linear dynamics (4.2) with the quadratic cost (4.6). We implement $u(1)$ and we restart the computations.

The nonlinear system (4.1) is used as a simulation of the data which will be replaced by the experimental data during the endurance session.

# 5. Estimation of the parameters in the design of the electrostimulator

## 5.1. Notations and definitions

The force fatigue model is written shortly

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)), \qquad (5.1) $$

where $x = (x_1, x_2, x_3, x_4, x_5)^T = (c_N, F, A, \alpha_A, \tau_{fat})^T$ and $u$ represents the FES input which can be smoothed as $u=u_{smooth}$.

The full system is defined by extending the dynamics with $\frac{d\alpha_A}{dt} = \frac{d\tau_{fat}}{dt} = 0$. We denote by $h = (h_1, h_2) = (F, A)$ the observation mapping.

Fixing a smooth control $u(t)$, the system with $x(0) = 0$ defines a control trajectory pair $(x(\cdot), u(\cdot))$ and we denote in short the Lie derivative $L_{X+uY}h(x(t)) = \frac{d}{dt}h(x(t))$. We denote by $O(x)$ the observation space formed by the iterated functions $\{L_{X+uY}^{k}h_i; i = 1, 2, k = 0, +\infty\}$. The system is called $u$-(weakly) observable if $x \to dO(x)$ is of full rank=dimension of the state space. Given a smooth input $u$ the system is called locally observable if there exists a sequence $0, ..., k_1, 0, ..., k_2$ so that the mapping $x \to \Phi(x, u) = [h_1(x), ..., L_{X+uY}^{k_1}h_1(x), h_2(x), ..., L_{X+uY}^{k_2}h_2(x)]$ is a diffeomorphism with respect to $x$ for all $(x, u)$ in an nonempty set $\chi \times U$ where $U$ contains the $k-1$ derivatives of $u$, with $k = \max(k_1, k_2)$. We say that $\chi \times U$ is an observable set. The construction of the observer is described in full details in [12] and is presented shortly in the next section.

## 5.2. Construction of the observer

Assume that the control trajectory pair $(x, u) \in \Omega_x \times \Omega_u \subset \chi \times U$. Perform the nonlinear change of coordinates $z = \Phi(x, u)$ and construct the observer

$$ \frac{dz}{dt}(t) = \mathrm{Pro}(Az + \rho(z, u) + S^{-1}K_0(y - Cz), z). \qquad (5.2) $$

The triple $(A, C, \rho(z, u))$ is obtained writing the system is the coordinates $z$. The matrix $K_0$ is chosen so that $(A - K_0C)$ is Hurwitz and $\mathrm{Proj}(y, z)$ is the projection operator associated to

$$ p(z) = \frac{\|z - z_0\|^2 - r_{\Omega}^2}{\alpha^2 + 2\alpha r_{\Omega}}. $$
$N_i$ being the backward identification horizon. Figure 2 represents the force, the mean force and the backward identification horizon to be found in order to get the best parametric model. The same figure could be constructed for A. $F_{mean_k}$ and $A_{mean_k}$ are calculated as following:

$$ 
\begin{aligned}
F_{mean_k} &= \frac{1}{t_{k+1}-t_k} \int_{t_k}^{t_{k+1}} F(\xi) d\xi \\
A_{mean_k} &= \frac{1}{t_{k+1}-t_k} \int_{t_k}^{t_{k+1}} A(\xi) d\xi
\end{aligned}
\qquad (4.5) $$

The criterion (4.4) traduces the fact that MA and MB are updated at each iteration ($MA_k, MB_k$). The couple is used in MPC strategy to calculate the control value $\eta_k$. The frequency of the stimulation being fixed.

## 4.2. MPC Algorithm

We present a version of the algorithm to illustrate the procedure which is classical, and the quadratic cost can be modified. The control constraints have been relaxed and they have to be introduced later to define the true feedback control.

We fix a sequence $k = 1, ..., K$ where $N_p$ is the prediction horizon and an output reference trajectory $y_{ref}$ associated to regulation of the force response F to a fixed level $F_{max}/\rho$ with $\rho > 1$ and a fatigue reference $A_{ref}(.)$. Denoting by $e(k) = (y(k) - y_{ref}(k))$, we minimize a cost of the form

$$ J(y, u) = \sum_{k=1,K} \lambda_1 \| e(k) \|_2^2 + \lambda_2 \| \Delta u_k \|_2^2, \qquad (4.6) $$

where $\Delta u(k)$ is the control increment and $\lambda_i$ are weighting parameters.

The feedback control $u(k)$ is computed on the horizon $K$ solving the LQ-problem defined by the linear dynamics (4.2) with the quadratic cost (4.6). We implement $u(1)$ and we restart the computations.

The nonlinear system (4.1) is used as a simulation of the data which will be replaced by the experimental data during the endurance session.

# 5. Estimation of the parameters in the design of the electrostimulator

## 5.1. Notations and definitions

The force fatigue model is written shortly

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)), \qquad (5.1) $$

where $x = (x_1, x_2, x_3, x_4, x_5)^T = (c_N, F, A, \alpha_A, \tau_{fat})^T$ and $u$ represents the FES input which can be smoothed as $u = u_{smooth}$.

The full system is defined by extending the dynamics with $\frac{d\alpha_A}{dt} = \frac{d\tau_{fat}}{dt} = 0$. We denote by $h = (h_1, h_2) = (F, A)$ the observation mapping.

The Lie bracket of two vector fields is defined by $[X, Y](x) = \frac{\partial X}{\partial x}(x)Y(x) - \frac{\partial Y}{\partial x}(x)X(x)$. Fixing a smooth control $u(t)$, the system with $x(0) = 0$ defines a control trajectory pair $(x(.), u(.)$ and we denote in short the Lie derivative $L_{X+uY}h(x(t)) = \frac{d}{dt}h(x(t))$. We denote by $O(x)$ the observation space formed by the iterated functions $\{L_{X+uY}^{k}h_i; i = 1, 2, k = 0, +\infty\}$. The system is called $u$- (weakly) observable if $x \to dO(x)$ is of full rank=dimension of the state space. Given a smooth input $u$ the system is called locally observable if there exists a sequence $0, ..., k_1, 0, ..., k_2$ so that the mapping $x \to \Phi(x, u) = [h_1(x), ..., L_{X+uY}^{k_1}h_1(x), h_2(x), ..., L_{X+uY}^{k_2}h_2(x)]$ is a diffeomorphism with respect to $x$ for all $(x, u)$ in an nonempty set $\chi \times U$ where $U$ contains the $k-1$ derivatives of $u$, with $k = \max(k_1, k_2)$. We say that $\chi \times U$ is an observable set. The construction of the observer is described in full details in [12] and is presented shortly in the next section.

## 5.2. Construction of the observer

Assume that the control trajectory pair $(x, u) \in \Omega_x \times \Omega_u \subset \chi \times U$. Perform the nonlinear change of coordinates $z = \Phi(x, u)$ and construct the observer

$$ \frac{dz}{dt}(t) = \mathrm{Pro}(Az + \rho(z, u) + S^{-1}K_0(y - Cz, z)). \qquad (5.2) $$

The triple $(A, C, \varrho(z, u))$ is obtained writing the system is the coordinates $z$. The matrix $K_0$ is chosen so that $(A - K_0C)$ is Hurwitz and $\mathrm{Proj}(y, z)$ is the projection operator associated to

$$ p(z) = \frac{\| z - z_0 \|^2 - r_{\Omega}^2}{\alpha^2 + 2\alpha r_{\Omega}}. $$
The point $z_0$ is the center of the domain $\Omega_z = B(z_0, r_\Omega)$ contained in $\Phi(\chi \times U)$ and $\alpha$ is an arbitrarily small positive constant. The construction involves block diagonal matrices including $S$ described in [12]. Note that it is introduced in relation with uniform linearization which is related in our construction to the uniform construction of the observable canonical form.

### 5.3. Geometric application

The experiments show that among the set of parameters the parameters $(\tau_c, R(0))$ are fixed and not depending upon the individual. Hence in particular the $Ca^{++}$ concentration $c_N$ can be taken as the control variable and can be chosen smooth according to a smooth FES-signal or taking the averaged value $c_{Naveraged}(t) = \frac{1}{t} \int_0^t c_N(s)ds$ over any subinterval of the training period. The bad input behavior is related to $c_N = 0$, in computing the inverse mapping of the map $z = \Phi(x, u)$. Hence note that the observer can be turned off imposing that: $\alpha \le c_N \le \beta$. At low level of stimulations corresponding to rest period one can rescaled $\alpha \to \infty$ $\alpha$ and expand the $F$-dynamics described by the Michaelis-Menten-Hill functions in Taylor Series at $c_N =0$, at a given order. This will reduce the computational complexity of the Lie derivatives which involve the derivative of the Michaelis-Menten-Hill functions with respect to $c_N$ and the time derivative of the concentration. A test input function of the form $c_N(t) = a + b\sin(\omega t)$ where a, b chosen so that the concentration stays in a arbitrarily band domain.

Additional parameters $A_{rest}, K_m$ are depending upon the individuals and can be estimated using the observer (5.2) during a single train $[0, T]$ using the force sensor only.

## 6. Simulation results

### 6.1. System identification

To identify the parametric model which will be used to calculate the MPC based control strategy, we use the Ding et al. model instead of real force and fatigue values (coming from experiments). The parametric model (linear model) will traduce locally the behaviour of the muscle, and needs to be updated for each new interpulse using, in our case, a variable identification moving horizon.

<figure><img src="image_3.png" /><figcaption>**Fig. 3** Evolution of the force (Ding et al. model) over a 1 second stimulation period with a 50ms interpulse interval, mean force values (for each interpulse) and identified mean force values over time</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 4** Evolution of the fatigue (Ding model), mean fatigue values and identified mean fatigue values over time</figcaption></figure>

Figures 3 and 4 show the evolution of the force and the fatigue based on the Ding et al. model over a 1 second stimulation period, with a 50ms interpulse interval, respectively. These figures clearly display the lobes generated by this stimulation. The mean values of this force for each interpulse, as well as those obtained by the least squares method, are also shown. The identified mean force value fits well the mean force value over identification horizon.
The point $z_0$ is the center of the domain $\Omega_z = B(z_0, r_\Omega)$ contained in $\Phi(\chi \times U)$ and $\alpha$ is an arbitrarily small positive constant. The construction involves block diagonal matrices including $S$ described in [12]. Note that it is introduced in relation with uniform linearization which is related in our construction to the uniform construction of the observable canonical form.

### 5.3. Geometric application

The experiments show that among the set of parameters the parameters $(\tau_c, R(0))$ are fixed and not depending upon the individual. Hence in particular the $Ca^{++}$ concentration $c_N$ can be taken as the control variable and can be chosen smooth according to a smooth FES-signal or taking the averaged value $c_{Naveraged}(t) = \frac{1}{t} \int_0^t c_N(s)ds$ over any subinterval of the training period. The bad input behavior is related to $c_N = 0$, in computing the inverse mapping of the map $z = \Phi(x, u)$. Hence note that the observer can be turned off imposing that: $\alpha \le c_N \le \beta$. At low level of stimulations corresponding to rest period one can rescaled $\alpha \to \infty$ $\alpha$ and expand the F-dynamics described by the Michaelis-Menten-Hill functions in Taylor Series at $c_N = 0$, at a given order. This will reduce the computational complexity of the Lie derivatives which involve the derivative of the Michaelis-Menten-Hill functions with respect to $c_N$ and the time derivative of the concentration. A test input function of the form $c_N(t) = a + b\sin(\omega t)$ where $a, b$ chosen so that the concentration stays in a arbitrarily band domain.

Additional parameters $A_{rest}, K_m$ are depending upon the individuals and can be estimated using the observer (5.2) during a single train $[0, T]$ using the force sensor only while the estimation of the fatigue parameters $(\alpha_A, \tau_{fat})$ needs to use the force and fatigue sensors.

## 6. Simulation results

### 6.1. System identification

To identify the parametric model which will be used to calculate the MPC based control strategy, we use the Ding et al. model instead of real force and fatigue values (coming from experiments). The parametric model (linear model) will produce locally the behaviour of the muscle, and needs to be updated for each new interpulse using, in our case, a variable identification moving horizon.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Evolution of the force (Ding et al. model) over a 1 second stimulation period with a 50ms interpulse interval, mean force values (for each interpulse) and identified mean force values over time</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Evolution of the fatigue (Ding model), mean fatigue values and identified mean fatigue values over time</figcaption></figure>

Figures 3 and 4 show the evolution of the force and the fatigue based on the Ding et al. model over a 1 second stimulation period, with a 50ms interpulse interval, respectively. These figures clearly display the lobes generated by this stimulation. The mean values of this force for each interpulse, as well as those obtained by the least squares method, are also shown. The identified mean force value fits well the mean force value over identification horizon.
## 6.2. Model predictive control

<figure><img src="image_5.png" /><figcaption>Fig. 5 Evolution of the force and fatigue over pulses at a reference of F=350N</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6 Evolution of the force and fatigue over pulses at a reference of F=230N</figcaption></figure>

Figures 5 and 6 illustrates the evolution of force and fatigue over pulses, computed from the Ding et al. model, over a stimulation period of 5 seconds with an interpulse interval of 30ms. The reference values for fatigue remain constant at 3.009N/ms, while the reference values for force are respectively 350N (5) and 230N (6). As expected, the MPC strategy allows to fit the force references while minimizing the difference between the fatigue and the fatigue rest value.

## 7. Conclusion

In this brief article we present the main steps in the design of a smart electrostimulator in relation with the construction of an industrial prototype: model, training sessions and estimation of the parameters using the physical sensors. Numerical simulations are presented for the MPC algorithm implemented to regulate the force and fatigue using a parametric model. The parameters are identified using the data of the Ding et al. model and will be replaced in fine by the experimental data. The Ding et al. model can be used to implement a NMPC algorithm where the parameters are estimated using an observer. But the method is computationally expensive and MPC
## 6.2. Model predictive control

<figure><img src="image_6.png" /><figcaption>Fig. 5 Evolution of the force and fatigue over pulses at a reference of F=350N</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 Evolution of the force and fatigue over pulses at a reference of F=230N</figcaption></figure>

Figures 5 and 6 illustrates the evolution of force and fatigue over pulses, computed from the Ding et al. model, over a stimulation period of 5 seconds with an interpulse interval of 30ms. The reference values for fatigue remain constant at 3.009N/ms, while the reference values for force are respectively 350N (5) and 230N (6). As expected, the MPC strategy allows to fit the force references while minimizing the difference between the fatigue and the fatigue rest value.

## 7. Conclusion

In this brief article we present the main steps in the design of a smart electrostimulator in relation with the construction of an industrial prototype: model, training sessions and estimation of the parameters using the physical sensors. Numerical simulations are presented for the MPC algorithm implemented to regulate the force and fatigue using a parametric model. The parameters are identified using the data of the Ding et al. model and will be replaced in fine by the experimental data. The Ding et al. model can be used to implement a NMPC algorithm where the parameters are estimated using an observer. But the method is computationally expensive for the design
algorithm based on linear parametric model can be chosen to tackle computational time while giving good results in terms of force and fatigue control.

References

[1] Toufik Bakir, Bernard Bonnard, Jérémy Rouot. A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model. Netw. Hetero. Media, 14(1) 79–100, 2019.

[2] Toufik Bakir, Bernard Bonnard, Loic Bourdin, Jérémy Rouot. Pontryagin-type conditions for optimal muscular force res-ponse to functional electrical stimulations. J. Optim. Theory Appl., 184 581–602, 2020.

[3] Toufik Bakir, Bernard Bonnard, Sandrine Gayrard, Jérémy Rouot. Finite Dimensional Approximation to Muscular Response in Force-Fatigue Dynamics using Functional Electrical Stimulation Automatica, 2022, (10.1016/j.automatica.2022.110464).

[4] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains. J. Appl. Physiol., 88 917–925, 2000.

[5] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. A predictive fatigue model. I. Predicting the effect of stimulation frequency and pattern on fatigue. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 10(1), 48–58, 2002.

[6] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. A predictive fatigue model. II. Predicting the effect of resting times on fatigue. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 10(1) 59–67, 2002.

[7] Joseph E. Gaudio et al.. Parameter estimation in adaptive control and time-varying systems under a range of excitation conditions, ArXiv:1911.03810v3, November 17th, 2021.

[8] Jean-Paul Gauthier, Ivan Kupka. Deterministic Observation Theory and Applications, Cambridge University Press, 2001.

[9] Rudolf Gesztelyi et al.. The Hill equation and the origin of quantitative pharmacology. Archive for history of exact sciences 66(4) 427–438, 2012.

[10] Alberto Isidori. Nonlinear Control Systems, Springer Verlag London, 1995.

[11] Jean-Baptiste Pomet and Laurent Praly. Adaptative nonlinear regulation: Estimation from the Lyapunov equation, IEEE Trans. on Automatic Control, 37:729-740, October 1992.

[12] Domitilla Del Vecchio and Richard M. Murray. Observability and Local Observer Construction for Unknown Parameters in Linearly and Nonlinearly Parameterized Systems, American Control Conference, 2003.

[13] Yang Wang, Stephen P. Boyd. Fast model predictive control using online optimization, IEEE Transactions on Control Systems Technology, 18 no. 2, 267–278, 2010.
of the prototype. MPC algorithm based on linear parametric model has been chosen in the industrial project to
tackle computational time while giving good results in terms of force and fatigue control.

References

[1] Toufik Bakir, Bernard Bonnard, Jérémy Rouot. A case study of optimal input-output system with sampled-data control: Ding et al. force and fatigue muscular control model. Netw. Hetero. Media, 14(1) 79–100, 2019.

[2] Toufik Bakir, Bernard Bonnard, Loic Bourdin, Jérémy Rouot. Pontryagin-type conditions for optimal muscular force res-ponse to functional electrical stimulations. J. Optim. Theory Appl., 184 581–602, 2020.

[3] Toufik Bakir, Bernard Bonnard, Sandrine Gayrard, Jérémy Rouot. Finite Dimensional Approximation to Muscular Response in Force-Fatigue Dynamics using Functional Electrical Stimulation Automatica, 2022, (10.1016/j.automatica.2022.110464).

[4] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains. J. Appl. Physiol., 88 917–925, 2000.

[5] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. A predictive fatigue model. I. Predicting the effect of stimulation frequency and pattern on fatigue. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 10(1), 48–58, 2002.

[6] Jun J. Ding, Anthony S. Wexler, Stuart A. Binder-Macleod. A predictive fatigue model. II. Predicting the effect of resting times on fatigue. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 10(1) 59–67, 2002.

[7] Joseph E. Gaudio et al.. Parameter estimation in adaptive control and time-varying systems under a range of excitation conditions, ArXiv:1911.03810v3, November 17th, 2021.

[8] Jean-Paul Gauthier, Ivan Kupka. Deterministic Observation Theory and Applications, Cambridge University Press, 2001.

[9] Rudolf Gesztelyi et al.. The Hill equation and the origin of quantitative pharmacology. Archive for history of exact sciences 66(4) 427–438, 2012.

[10] Alberto Isidori. Nonlinear Control Systems, Springer Verlag London, 1995.

[11] Jean-Baptiste Pomet and Laurent Praly. Adaptative nonlinear regulation: Estimation from the Lyapunov equation, IEEE Trans. on Automatic Control, 37:729-740, October 1992.

[12] Domitilla Del Vecchio and Richard M. Murray. Observability and Local Observer Construction for Unknown Parameters in Linearly and Nonlinearly Parameterized Systems, American Control Conference, 2003.

[13] Yang Wang, Stephen P. Boyd. Fast model predictive control using online optimization, IEEE Transactions on control systems technology, 18 no. 2, 267–278, 2010.