# Fast Optimization Scheme for the Muscular Response to FES Stimulation to Design a Smart Electrostimulator

Sandrine Gayrard

► To cite this version:

Sandrine Gayrard. Fast Optimization Scheme for the Muscular Response to FES Stimulation to Design a Smart Electrostimulator. 2023 European Control Conference (ECC), Jun 2023, Bucharest, France. pp.1-6, 10.23919/ECC57647.2023.10178228 . hal-03891925

HAL Id: hal-03891925

https://inria.hal.science/hal-03891925

Submitted on 13 Dec 2022

**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
# Fast Optimization Scheme for the Muscular Response to FES Stimulation to Design a Smart Electrostimulator

Sandrine Gayrard<fnref n="1" />*

**Abstract**— In this article, we present the design of a smart electrostimulator for muscle rehabilitation or reinforcement, using fast computations, in order to control the muscular force. The Ding and al. model allows to predict and to optimize the muscular force response to functional electrical stimulation. [1]. We analyze the estimation of the Ding and al. parameters using an approximation of the force response [6] which depends upon the 6 parameters of the Ding's model and we derive optimization scheme which bypass the time computational expensive integration of the dynamics of the Ding and al. equations.

## I. INTRODUCTION

### A. Industrial Project

The aim of the study project is the design of a smart electrostimulator, in the isometric case (without displacement of the muscle), for muscle rehabilitation or reinforcement based, in particular, on modern sensors and fast computations using embedded electronics. The idea is to design training programs to track a force of reference [5] or to maximize the muscular force. Each training program can be converted into an optimization problem. The mathematical equations presented below are the tools to the design of this smart electrostimulator.

### B. Mathematical equations

The Ding and al. model [1] allows to predict and to optimize the muscular force response to functional electrical stimulations and was validated experimentally [2]. The system is described as follow.

The input $u$ of a pulse train is defined, for $t \in [0, T]$, by: $u(t) = \sum_{i=0}^{n} \delta(t - t_i)$, with $0 = t_0 < t_1 < \dots < t_n < T$ the impulse times with $n \in \mathbb{N}$ fixed. This physical control allows us to obtain the signal $E(t)$, which permits to calculate the force response following a train of electrical pulses. The dynamics is given by (1) :

$$ \dot{E}(t) + \frac{E(t)}{\tau_c} = \frac{1}{\tau_c} \sum_{i=0}^{n} \eta_i R_i \delta(t - t_i), \quad (1) $$

with $E(0) = 0$, and depending on the response time $\tau_c$ and the function $R_i$ defined by $R_i = 1$ if $i = 0$ and $R_i = 1 + (R_0 - 1) \exp(-\frac{t_i - t_{i-1}}{\tau_c})$ elsewhere. It allows the visualization of the tetany phenomenon, which corresponds to the memory of the effect of successive pulses. $\eta_i$ is the amplitude of the electric pulse stimulation.

The FES (Functional Electrical Stimulation) signal allows to determine the evolution of electrical conduction as a function of the dynamics describing the evolution of the concentration $c_N$ of the calcium ions ($Ca^{2+}$), which plays an important role in muscle contraction.

$$ \dot{c}_N(t) + \frac{c_N(t)}{\tau_c} = E(t) \quad (2) $$

The others equations of the Ding and al. model are described below and provide the muscular force response.

$$ m_1(t) = \frac{c_N(t)}{K_m + c_N(t)}, \quad m_2(t) = \frac{1}{\tau_1 + \tau_2 m_1(t)}, \quad (3) $$

with $m_1$ the Michaelis-Menten (1913) function. The muscular force response satisfied the Hill-Huxley's models dynamics.

$$ \dot{F}(t) = -m_2(t)F(t) + m_1(t)A, \quad (4) $$

with $\tau_c, R_0, A_0, K_m, \tau_1, \tau_2$ the set of parameters described in [1].

### C. Motivation

The controls are the frequency and amplitudes of the electrical stimulations. They fit in the sampled-data control category [6]. The goal is to maximize the force response or to track a force of reference [5].

In previous references [3], [4] the problem was analyzed using direct or indirect optimal control schemes which are computationally expensive due to the numerical integration of the dynamics. In a more recent work we derive an approximation of the force response in [6] which depends upon the 6 parameters of the Ding's model and which can be used to derive a fast optimization scheme in this category by bypassing the time computational expensive integration of the dynamics.

## II. APPROXIMATION OF THE MUSCULAR FORCE RESPONSE TO FES

We present a brief recap of the [6]. We derive an approximation of the force response, which gives an explicit form of the dynamics of the force as a function of time.

1) Definition 1: The muscular force, with $F(0) = 0$, is defined in $[t_k, t_{k+1}]$, $k = 0, \dots, n$, by

$$ F(t) = A_0 M(t) \int_0^t M^{-1}(s) m_1(s) s, \quad (5) $$

with $M(t) = \exp(-\int_0^t m_2(s) s)$.

<footnote>*This work was supported by Segula Matra Technologies organization, Trappes (France)</footnote>

<footnote n="1">¹PhD student, Burgundy Unviversity of Dijon (France), INRIA McTAO (France) sandrine.gayrard@inria.fr</footnote>
We define a finer partition of $(t_i)_{1 \le i \le n}$ noted $(t_{i+j/p})$, $i = 0, \dots, n$, $j = 0, \dots, p-1$, with $p \in \mathbb{N}^*$ the number of the new points inserted between $t_i$ and $t_{i+1}$. Using (5), we construct the approximation of the force $F$ on $[t_{k+j/p}, t_{k+(j+1)/p}]$, $k = 0, \dots, n$, $j = 0, \dots, p-1$, by :

$$F(t) \approx A_0 \tilde{M}(t) \int_0^t \tilde{M}^{-1}(s) \tilde{m}_1(s) s, \quad \forall t \in [t_{k+j/p}, t_{k+(j+1)/p}]$$

with $\tilde{M}(t) = \exp(-\int_0^t \tilde{m}_2(s) s)$, et $\tilde{m}_i$, $i = 1, 2$ an approximation of $m_i$ defined by :

$$\tilde{m}_i(t) = a_{ij,k} (t - t_{k+j/p}) + b_{ij,k}, \quad \text{for } t \in [t_{k+j/p}, t_{k+(j+1)/p}]$$

and $k = 0, \dots, n$, $j = 0, \dots, p-1$.

The constants $a_{ij,k}$, $b_{ij,k}$, $j = 0, \dots, p-1$ are calculated at specific times, depending on the choice on the approximation of $m_i$.

As $\dot{m}_2 = \tau_2 \dot{m}_1 / (\tau_1 + \tau_2 m_1)^2$, then $\dot{m}_2$ is equal to 0 when $c_N$ is maximal. On $[t_k, t_{k+1}]$, $m_2$ can be approximated by $\tilde{m}_2^k(t) = a_{ij,k} (t - t_{k+j/2}) + b_{ij,k}$, $k = 0, \dots, n$ and $j = 0, 1$, with $t_{k+1/2} = \underset{t \in [t_k, t_{k+1}]}{\operatorname{argmax}} \tilde{c}_N(t)$. By calculation, we end up with

$$a_{ij,k} = \frac{m_i^k(t_{k+1/2}) - m_i^k(t_k)}{t_{k+1/2} - t_k}, \quad b_{ij,k} = m_i^k(t_k).$$

2) **Proposition 1:** The approximation of the muscular force $F$ on $[t_{k+j/p}, t_{k+(j+1)/p}]$, $k = 0, \dots, n$, $j = 0, \dots, p-1$ is given by :

$$\begin{aligned} \tilde{F}(t)/A_0 = & \sum_{i=0}^{k-1} \sum_{j=0}^{p-1} \int_{t_{i+j/p}}^{t_{i+(j+1)/p}} \tilde{M}(t) \tilde{M}^{-1}(s) \tilde{m}_1(s) s \\ & + \sum_{j=0}^{j-1} \int_{t_{k+j/p}}^{t_{k+(j+1)/p}} \tilde{M}(t) \tilde{M}^{-1}(s) \tilde{m}_1(s) s \\ & + \int_{t_{k+j/p}}^{t} \tilde{M}(t) \tilde{M}^{-1}(s) \tilde{m}_1(s) s. \end{aligned} \quad (6)$$

3) **Proposition 2:** From a practical point of view, we just need to estimate the values of the force $F$ for $t_{k+j/p}$, which gives :

$$\begin{aligned} \tilde{F}(t_{k+j/p})/A_0 = & \sum_{i=0}^{k-1} \sum_{j=0}^{p-1} \int_{t_{i+j/p}}^{t_{i+(j+1)/p}} \tilde{M}(t_{k+j/p}) \tilde{M}^{-1}(s) \tilde{m}_1(s) s \\ & + \sum_{j=0}^{j-1} \int_{t_{k+j/p}}^{t_{k+(j+1)/p}} \tilde{M}(t_{k+j/p}) \tilde{M}^{-1}(s) \tilde{m}_1(s) s \end{aligned} \quad (7)$$

# III. SMART ELECTROSTIMULATOR

## A. Industrial Project

Let's remind that the goal is the design of smart electrostimulator for muscle rehabilitation or reinforcement, using modern sensors and fast computations using embedded electronics, in order to control the muscular force [5].

<figure><img src="image_2.png" /><figcaption>Fig. 1. Construction of the approximation of $\tilde{F}(t)$ and $\tilde{F}(t_{k+j/p})$ (red curve) with $k = 0, \dots, n$, $j = 0, \dots, p-1$ and $p \in \mathbb{N}^*$ the number of the new points added between $t_k$ and $t_{k+1}$. The black curve is the force calculated with an Euler integration scheme. The constants are $p = 2$ and $n = 6$ impulses.</figcaption></figure>

## B. Materials and methods

The measurements of the muscular force were carried out under isometric conditions, i.e. without displacement of the muscle. The measurements were performed at the biceps brachii (BB) level. The stimulation protocol and data acquisition are also presented. The subject is placed on an adjustable chair. His/her arm is placed on a L-shaped platform so that the arm is in a 90 degree position and the chair is adjusted so that the arm is completely on the platform and is "glued" to it to allow stabilization with a strap hooked around the biceps. The experimental setup is presented Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Experimental setup</figcaption></figure>

1) **Sensor:** The force sensor is an FSR sensor (ref : Interlink Eletronics FSR 402, 406). The resistance of the sensor decreases as the force applied to its surface increases. The pressure sensor is placed between the strap and the
biceps, in such way that when the biceps contracts, the biceps compresses the sensor, in a perpendicularly direction.

The range of force sensitivity varies according on the chosen model. The sensor is connected to an Arduino UNO to sample the output voltage and to process it.

Finally, a premilinary calibration to grade the relation between the sensor output voltage and the induced force, is performed.

2) Stimulation: Electrical stimulation is provided by a commercial electrostimulation device (ref : Beurer EM 49 Digital TENS/EMS unit, see Fig. 2) designed for EMS (Electrical Muscle Stimulation) to promote muscle regeneration. Its adjustable parameters are the duration of the stimulation train $t_{stim}$ (fixed at 15 ms), the duration of rest $t_{rest}$ (fixed at 10ms), the width of the pulse (fixed at 150µs). The stimulation frequency f is chosen according to the objective of the experiment and the amplitude is modifiable during the experiment, as the controls.

### C. Data acquisition

During the muscle contraction, determined by its stimulation amplitude, the biceps compresses the force sensor. Isometric force is sampled at a sampling rate of 1kHz and recorded through a serial port driven by a Matlab program.

### D. Training program : algorithm

Each training program can be translated into an optimization problem. The several steps are described below and resumed in Fig. 3

* A succession of pulse trains is send to the biceps. Each train (of frequency f) lasts $t_{stim}$ and is followed by a rest time of $t_{rest}$.

* During the pulse train sent to the biceps, the sensor measures the induced force.

* At the end of the pulse train, the resting time starts. During this resting time, the force measurements are recovered and processed, by a Matlab program running on a standard computer.

* The force model parameters are estimated by fitting the measured data to the force prediction model described by the muscle force approximation.

* From the force prediction model and the parameters estimated before, the optimal control is calculated according to an optimization problem that allows to track a force of reference or to maximize the muscular force.

* This new calculated optimal command is then applied to the muscle during the next impulses train, at the end of the resting time.

The two main parts are : 1) the muscle parameter's estimation; 2) the computation of the optimal control, using the approximation model force. They are both computed during the resting time $t_{rest}$, which reinforces the need for fast computations.

## IV. ESTIMATION OF THE DING AND AL. PARAMETERS

The muscle parameters of each person are estimated during the resting periods with low frequency stimulation using the piezoelectric force sensor.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Successive steps of the training program.</figcaption></figure>

### A. Data-fitting

1) Model: The explicit form of the force approximation gives the force in terms $\alpha$ a set 6 parameters of Ding and al. More precisely, (7) gives the approximation of the force in terms of a set $\alpha$ of 6 parameters and $t_{k+j/p}$ with $k = 0, ..., n$ ($n$ number of pulses during the train), $j = 0, ..., p - 1$ and $p \in \mathbb{N}^*$ the number of the new points added between $t_k$ and $t_{k+1}$. $\alpha = [\tau_c, R_0, A_0, K_m, \tau_1, \tau_2]$.

2) Data: Force data were recorded using the experimental set-up. The input, a $n^{th}$ pulse train with a stimulation frequency f gives the muscular response with as many bounces as pulses. For each bounce k, the average value is calculated using the data value $F_{data}$:

$$ \bar{F}_{[t_k, t_{k+1}]} = \frac{1}{t_{k+1} - t_k} \int_{t_k}^{t_{k+1}} F_{data}(t) dt \quad (8) $$

The vector of data force, which is going to be used to fit the model and estimate the 6 parameters of Ding and al. can be written for each impulse $k \in [1, n - 1]$ by $F_{data-fit} = [F_{data}(t_k), \bar{F}_{[t_k, t_{k+1}]}]$ (9) and is represented Fig. 4.

$$ F_{data-fit} = [F_{data}(t_1), \bar{F}_{[t_1, t_2]}, \dots, F_{data}(t_{n-1}), \bar{F}_{[t_{n-1}, t_n]}] \quad (9) $$

All the information of the model (the 6 parameters) is contained in the muscular response of the impulse train of size $n$. However, because of too high calculation time, only the muscular response to 5 or 6 impulses is taken into account, thus requiring a compromise between amount of information and computing time.

<figure><img src="image_5.png" /><figcaption>Fig. 4. At the right, $F_{fit}$ with the dashed curves corresponding to the times $t_{k+j/p}$. At the left, $F_{data}(t)$ and $F_{data-fit}$ the vector of values which contains the data that will be fitted to the model. $p = 2$ and $n = 6$.</figcaption></figure>
3) **Cost Function:** Each value of $F_{data-fit}$ fits the model $F(t_{k+j/p})$ which corresponds to the approximation force calculated for each impulse time $t_k$ and for $t_{k+j/p}$, the intermediate times between 2 pulses (see Fig. 4). For $k \in [1,n]$, the values $F_{data}(t_k)$ fit $\tilde{F}(t_k)$ and the values $\tilde{F}_{[t_k, t_{k+1}]}$ fit $\tilde{F}(t_{k+j/p})$ ($p=2$ and $j=1$).

The cost function to minimize is described below:

$$J(\alpha) = \sum_{i=1}^{N} \| \tilde{F}_i(t_{k+j/p}, \alpha) - F_{i,\text{data-fit}} \|_{2}^{2} \quad (10)$$

with $\tilde{F}$ the approximation force that fits the values of force $F_{data-fit}$. $\alpha = [\tau_c, R_0, A_0, K_m, \tau_1, \tau_2]$. N the size of the data and the fit vectors ($N=n.p$).

## B. Practical estimation

1) **Computation:** The use of a non-linear solver, helps to solve non-linear curve fitting problems of the form $\min J(\alpha)$ with $\alpha$ the elements and $\alpha_0$ the initial vector of the parameters. The solver starts at $\alpha_0$ and attempts to find a local minimum $\alpha$ of the function described by $J$.

The associated problem is solved using a quasi-newton method, with Matlab [7], on a standard computer. We initialize the parameters with the Ding and al. parameters [1].

2) **Experimental results:** The estimated parameters are calculated for one subject with 4 different frequencies and various amplitudes. They are given in Fig. 6 and in Table I. For each frequency, identical stimulation have been send to the muscle BB and the parameters have been estimated for each train: the represented values correspond to the average values calculated using all the data, for each frequency. The corresponding curves for 3 various frequencies are depicted Fig. 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. The black curve corresponds to the force data measured with the FSR sensor. The dashed curve is the average force data measured at $t_{k+j/p}$ with $k=0, ..., n$, $j=0, ..., p-1$ and $p \in \mathbb{N}^*$ the number of the new points added between $t_k$ and $t_{k+1}$. The red curve represents the approximation of the force that fits the measured data. Values of the constants are $n=5$ and $p=2$. The computation time is less than 2s.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Estimated experimental parameters for one individual biceps, non-fatigued case, compared to the Ding and al. parameters, for 4 different frequencies: $f = [1, 2, 5, 10]$ Hz.</figcaption></figure>

<table><caption>TABLE I
ESTIMATED EXPERIMENTAL PARAMETERS FOR ONE INDIVIDUAL BICEPS FOR DIFFERENT STIMULATION FREQUENCY.</caption><thead><tr><th>α</th><th>τ<sub>c</sub> (ms)</th><th>R<sub>0</sub></th><th>K<sub>m</sub></th><th>A<sub>0</sub> (N/ms)</th><th>τ<sub>1</sub> (ms)</th><th>τ<sub>2</sub> (ms)</th></tr></thead><tbody><tr><td>1Hz</td><td>19.83</td><td>1.14</td><td>4.38</td><td>1.31</td><td>50.91</td><td>124.39</td></tr><tr><td>2Hz</td><td>19.86</td><td>1.14</td><td>4.82</td><td>1.14</td><td>50.88</td><td>124.40</td></tr><tr><td>5Hz</td><td>19.87</td><td>1.13</td><td>4.87</td><td>0.77</td><td>50.83</td><td>124.39</td></tr><tr><td>10Hz</td><td>19.82</td><td>1.13</td><td>5.56</td><td>0.94</td><td>50.86</td><td>124.39</td></tr><tr><td>α<sub>Ding</sub></td><td>20</td><td>1.14</td><td>0.1030</td><td>3.009</td><td>50.95</td><td>124.40</td></tr></tbody></table>

3) **Observation and discussion:** The errors between measured data $F_{data-fit}$ and predicted force values increase at lower pulse frequencies because the approximating force model considers force at a given time $t_{k+j/p}$ but does not take into account the dynamic behavior of the force with bounces, whose peak-to-peak amplitude increases for lower frequency.

The parameters $\tau_c$, $R_0$, $\tau_1$ and $\tau_2$ show a standard deviation with the Ding and al. parameters, respectively less than 26%, 3%, 12% and 1.5%. The difference appears for the parameters $K_m$ and $A_0$ with standard deviations greater than 50% and increases with the frequency of stimulation.

These two parameters correspond to the non-fatigued model. The experimental results show difficulties to predict results depending on the muscle fatigue. The biceps can be subjected to a significant effort without the subject having performed it, thus disturbing the measurements. Sources of error can be a lack of hydration, general fatigue or skin's condition. To overcome this problem, it is necessary to multiply the experiments, on the same subject and on other subjects.

The difference obtained with the parameter $A_0$ with the Ding's values comes from the fact that the latter depends on
the calibration of the force measurements. It is important to note that Ding et al. performed experiments mainly on the quadriceps. [8] shows that parameters varied with muscle length and fatigue.

A statistical analysis was done to understand the difference in data for each group of parameters for different frequencies. Using an R algorithm and after checking that the measurements fits either a normal distribution or a 'poisson' distribution, it was observed, using the ANOVA (analysis of variance) test that:

* For the parameter $\tau_1$, there are significant errors (p < 0.05) between the group of measurements for frequencies: $5Hz-1Hz$ with $p = 0.00014$ and $5Hz-2Hz$ with $p = 0.012$. In other words, the group of $\tau_1$ values at 5Hz is significantly different than at 1Hz and 2Hz.

* There are no significant errors between measurement groups of estimated parameters $\tau_c$, $\tau_2$, $A_0$, $R_0$ and $K_m$ ($p > 0.05$) for each frequency, which means that for each parameter, the groups of parameter measurements can be considered identical. This also shows that these 5 factors do not vary with stimulation frequency.

The most surprising result is the statistical analysis of $A_0$: this parameter describes muscle fatigue and decreases when the muscle is fatigued. The table of values shows a decrease in $A_0$ as the frequency increases, which is consistent with the fact the higher the frequency the more the muscle is subjected to electric impulses, the more it fatigues the muscle. However, the statistical analysis of these values describes that there is no significant differences between these values packages for the three different frequencies studied. It was also shown in the analysis of Ding et al. [9] that the parameters $\tau_1$ and $K_m$ change with muscle fatigue. A significant error was found between the parameters at 5 Hz versus 1 Hz and 2 Hz for $\tau_1$. But no significant errors ($p > 0.05$) were found for $K_m$. Another comment, only experiments with constant stimulation frequency were performed because, firstly, the stimulation equipment imposes this constraint and secondly it was shown that parameters vary significantly between constant and non-constant frequency stimulation, which is part of the limitation of the force model described by Ding in [8].

## V. OPTIMIZATION PROBLEM

The training programs are being translated into an optimization problem. Both problems fit in the sampled-data control category.

a) Endurance program: We consider a single train $[0, T]$ on which the corresponding problem is:

$$ \min_x \int_0^T |F(t) - F_{ref}|^2 dt \quad (11) $$

The amplitudes are optimization variables ($x = (\eta_1, \eta_2, ..., \eta_n)$), which belong to the finite dimensional input-space $I \in [0, 1]^{n+1}$ defined by amplitude constraints: $\eta_i \in [0, 1]$, $i = 0, ..., n$, with $n$ a fixed integer, representing the number of pulses in the train. $T$ is the time of the end of the pulse train. The reference force $F_{ref}$ is adjusted in relation with the user.

b) Punch program: We consider a single train $[0, T]$, $T$ free, on which the corresponding problem is:

$$ \max_x F(T) \quad (12) $$

Amplitudes are fixed at one. Impulse times are the optimization variables ($x = (t_1, t_2, ..., t_n)$). They belong to the finite dimensional input-space $I \in [0, T]^n$. They are subjected to two constraints:

* $0 = t_0 < t_1 < ... < t_n < t_{n+1} = T$, with $n$ a fixed integer, representing the number of pulses in the train.

* $t_{i+1} - t_i > \varepsilon = 10$ ms for $i = 0, ..., n$.

### A. Optimal sampled-data control problem

Using the equation of the approximation force (7), we can generate the explicit expression of the approximation of the force, in terms of the optimization variables $x$. $\tilde{J}_1$ and $\tilde{J}_2$ constitute the approximated cost functions of the optimization routine.

$$ \min_{\eta_i} \tilde{J}_1 = \min_{\eta_i} \sum_{k=1}^{n} \| \tilde{F}_k(t_{k+j/p}, \alpha, \eta_i) - F_{ref} \|_2^2 \quad (13) $$

$$ \min_{t_i} \tilde{J}_2 = \min_{t_i} (-\tilde{F}_k(T, \alpha, t_i)) \quad (14) $$

with $\alpha = [\tau_c, R_0, A_0, K_m, \tau_1, \tau_2]$.

For each case, the associated problem which gives the approximation $\tilde{J}$ of the cost function $J$ based on the approximation $\tilde{F}$ of the variable $F$ described before, is solved using an interior point method, with the function "fmincon" in Matlab, on a standard computer. For the first case, we initialize the amplitudes of the pulse train to random values on $[0, 1]$ and for the second case, we initialize the times of the pulse train to random values, respecting the two constraints described before.

1) Numerical results: 1.1) The optimal solution $\eta_i^*$ and the force response $\tilde{F}$ corresponding are depicted in Fig. 7.

<figure><img src="image_8.png" /><figcaption>Fig. 7. The optimal solution $(\eta_0^*, ..., \eta_n^*)$ minimizes the cost $\tilde{J}_1$ from (13). $\tilde{F}^+$ and $\tilde{F}^-$ are respectively the upper and the lower approximation of $F$, described in the article [6]. Values of the constants are $\tau_c = 20$ms, $K_m = 0.1030$, $R_0 = 1.14$, $A_0 = 3.14$ N/m, $\tau_1 = 54$ms, $\tau_2 = 127$ms, $n = 8$ and $F_{ref} = 0.1$ kN.</figcaption></figure>

1.2) The optimal solution $t_i^*$ and the force response $\tilde{F}$ corresponding are depicted in Fig. 8.
<figure><img src="image_9.png" /><figcaption>Fig. 8. The optimal solution $(t_0^*, ..., t_n^*)$ maximize the cost $\tilde{J}_2$ from (14). $\tilde{F}^+$ and $\tilde{F}^-$ are respectively the upper and the lower approximation of $F$, described in the article [6]. Values of the constants are $\tau_c = 20$ms, $K_m = 0.1030$, $R_0 = 1.14$, $A_0 = 3.14$ N/m, $\tau_1 = 54$ms, $\tau_2 = 127$ms, $n = 8$.</figcaption></figure>

## B. Smart electrostimulator

Only the endurance program has been implemented into the smart electrostimulator because of the hardware constraint which imposes a constant frequency during the stimulation train. Besides, maximizing muscle strength during the experimental tests generates pain. The training program, translated into an optimization problem, has been tested using the algorithm presented in Fig. 3. Optimal controls are the amplitudes $\eta$ of the stimulation train. In the case of the smart electrostimulator and for simplification purposes, we impose optimal solution $\eta_i^*$ constant for all $i \in [0, n]$ with n the number of stimulation impulses.

1) **Experimental results:** The optimal solution $\eta_i^*$, the predicted force response $\tilde{F}$ and the measured force $F_{data}$ are depicted in Fig. 9 and 10.

<figure><img src="image_10.png" /><figcaption>Fig. 9. The optimal solution $(\eta_0^*, ..., \eta_n^*) = 0.6615$ minimizes the cost $J$ from (13). The black curve is the approximation force $\tilde{F}$, the blue curve is $F_{measured}$ the data force measured with FSR sensor. The dashed curve is $F_{ref}$. Values of the constants are $\tau_c = 19.74$ms, $K_m = 6.86$, $R_0 = 1.10$, $A_0 = 0.68$, $\tau_1 = 50.84$ms, $\tau_2 = 124.39$ms, $n = 10$, $F_{ref} = 0.7$ N, the impulse frequency is constant and $f = 10$ Hz. Computation time is less than 3s.</figcaption></figure>

2) **Observation and discussion:** The error between measured data and predicted force values increases when the pulse frequency is lower. Indeed, the dynamic behavior of the force with bounces whose peak-to-peak amplitude increases for a lower frequency, makes the tracking more difficult. The most problematic fact comes from the computation of the

<figure><img src="image_11.png" /><figcaption>Fig. 10. The optimal solution $(\eta_0^*, ..., \eta_n^*) = 0.8817$ minimizes the cost $J$ from (13). The black curve is the approximation force $\tilde{F}$, the blue curve is $F_{measured}$ the data force measured with FSR sensor. The dashed curve is $F_{ref}$. Values of the constants are $\tau_c = 19.86$ms, $K_m = 4.94$, $R_0 = 1.14$, $A_0 = 1.20$, $\tau_1 = 50.88$ms, $\tau_2 = 124.39$ms, $n = 8$, $F_{ref} = 1.26$ N, the impulse frequency is constant and $f = 2$ Hz. Computation time is less than 3s.</figcaption></figure>

error between measured force data and $F_{ref}$ reference force. The error has to be computed either from the lower peak of force bounces or the higher peak of force bounces. If we introduce a total error margin of 5% around the reference value, the measured force are out of this range from 28% at the bottom and 71% with a low frequency. This error decreases when the stimulation frequency is higher. But a higher frequency means the perception of pain increasing and more importantly, the fatigue increases rapidly which means that the muscle can't maintain the reference force.

## VI. CONCLUSIONS

Our study shows how to estimate the set of 6 parameters of the Ding's model using an approximation of the force response, that allows to bypass the time computational expensive integration of the dynamics. This study is about the design of a smart electrostimulator with fast computations and modern sensors, allowing us to fit the model to the sensors values and to compute the optimal trajectory.

## ACKNOWLEDGMENT

The work presented in this article was performed in the frame of a PhD "Design of a smart muscle stimulator", as a CIFRE agreement between the Burgundy university (Dijon, France) and Segula Matra Technologies (Trappes, France) under the direction of Pr. Toufik Bakir, Pr. Bernard Bonnard and Mr. Mokrane Abdiche.

## REFERENCES

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

[2] E. Wilson, Force response of locust skeletal muscle, Southampton University, Ph.D. thesis, 2021.

[3] T. Bakir, B. Bonnard, L. Bourdin, J. Rouot, Pontryagin-type conditions for optimal muscular force response to functional electrical stimulations, J. Optim. Theory Appl., 2020, pp. 184:581-602.

[4] T. Bakir, B. Bonnard, L. Bourdin, J. Rouot, Direct and indirect methods to optimize the muscular force response to a pulse train of electrical stimulation, ESAIM: Proceedings and Survey, EDP Sciences, vol. 71, 2021, pp. 1-10.
[5] B. Doll, N. Kirsch, N. Sharma, Optimization of a stimulation train based on a predictive model of muscle force and fatigue, IFAC-PapersOnLine, vol. 48(20), 2015, pp. 338-342.

[6] T. Bakir, B. Bonnard, S. Gayrard, J. Rouot, Finite Dimensional Approximation to Muscular Response in Force-Fatigue Dynamics using Functional Electrical Stimulation, Automatica, 2022.

[7] The MathWorks. Curve Fitting Toolbox User's Guide. The MathWorks, Inc, Natick, 2001.

[8] J. Ding, A.S. Wexler, S.A Binder-Macleod, Two-step, predictive, isometric force model tested on data from human and rat muscles, J. Appl. Physiol., vol. 85, 1998, pp. 2176-2189.

[9] J. Ding, A.S. Wexler and S.A. Binder-Macleod, A predictive model of fatigue in human skeletal muscles, J. Appl. Physiol., 89, 2000, pp. 1322-1332.