# Optimal Control Framework for Personalized FES-cycling in Individuals with Spinal Cord Injury

Tiago Coelho-Magalhães, Christine Azevedo Coste, François Bailly

► To cite this version:

Tiago Coelho-Magalhães, Christine Azevedo Coste, François Bailly. Optimal Control Framework for Personalized FES-cycling in Individuals with Spinal Cord Injury. FESWS 2025 - 15th Vienna International Workshop on Functional Electrical Stimulation & 30 years IFESS Anniversary, Sep 2025, Vienna (AUSTRIA), Austria. **hal-05165086v2**

HAL Id: hal-05165086

https://inria.hal.science/hal-05165086v2

Submitted on 22 Jul 2025

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
# Optimal Control Framework for Personalized FES-cycling in Individuals with Spinal Cord Injury

Coelho-Magalhães T¹, Azevedo-Coste C¹, Bailly F¹

¹CAMIN, Inria, University of Montpellier, Montpellier, France

**Abstract:** Functional Electrical Stimulation (FES)-cycling is employed to facilitate pedalling in individuals with spinal cord injury (SCI), with the aim of promoting physical activity and enhancing overall health. However, the efficiency of FES-cycling remains limited due to the challenge of optimally timing muscle contractions to maximise power output while minimizing muscle fatigue. This study aims to explore the use of an optimal control framework to improve FES-cycling performance in simulation. We use a muscle model that accounts for the effects of electrical stimulation on muscle contraction dynamics and identify its parameters using experimental data collected from 1 individual with SCI. Isometric evoked contractions were measured using the crankset power meter of our tricycle instead of an isokinetic dynamometer. We then demonstrate the feasibility of using trajectory optimization to track a 4.4 rad/s (42 rpm) pedalling cadence and a 20 W crank power output by optimizing the electrical stimulation pulse duration applied on six muscles. Results show the power achieved is 17.53 ± 7.87 W and velocity 4.09 ± 0.36 rad/s. This study presents a first step towards improving the realism of FES-driven simulations, showcases the use of physiological model-based control of FES-cycling, and highlights the potential of predictive modelling and optimal control to personalize FES protocols.

**Keywords:** Biomechanics, Functional Electrical Stimulation, numerical optimal control, FES-cycling, musculoskeletal model

## Introduction

Functional Electrical Stimulation (FES) assisted cycling employs electrical stimulation pulses to evoke coordinated contractions of paretic or paralyzed muscles of the lower limbs to reproduce a pedaling movement, usually on an adapted tricycle [1]. It is a widely investigated rehabilitation strategy for people with neurological conditions, such as spinal cord injury (SCI), and has demonstrated to be efficient in improving muscle, bone and cardiovascular health, fitness, feelings of well-being, and motor functions [2]. Despite its benefits, effectively controlling evoked cycling movement through FES remains a major challenge due to the resulting rapid onset of muscle fatigue [3], and the non-selective activation of muscle groups that results in poor synergistic coordination of agonist and antagonist muscles around a joint [4]. To minimize these drawbacks, researchers aim to optimize the stimulation patterns (i.e., the pulse parameters and their activation timing within the pedalling cycle) and the control laws that dynamically adjust it to the activity goals. Yet, the solution is frequently sub-optimal and often obtained through trial and error or time-consuming procedures [1].

Numerical simulation can be used to explore stimulation patterns or advanced control strategies that generate optimized muscle forces and joint kinematics. In particular, trajectory optimization is an efficient method to explore novel control strategies, guided by the system's dynamics and task-related constraints. Recent advances in musculoskeletal modeling and numerical optimal control provide tools to simulate and optimize such movements [5]. In [6], a framework combining numerical optimal control and FES cycling to track different crank velocities and resistances is presented. The model in this study uses pure torques at the joints, avoiding muscle modeling, thus

lacking physiological accuracy related to muscle contraction dynamics. In a previous work [7], we adapted a physiological muscle model that predicts electrically-evoked muscular force to be used with gradient-based optimization, in particular with numerical optimal control/estimation problems. Experimental data obtained from individuals with SCI were used to identify the model's parameters and the feasibility of employing this adapted model in numerical optimal control problems was demonstrated.

In the following, this specific muscle formulation is used in a personalized FES-cycling model containing three different muscle groups in each leg. Their parameters are identified using experimental data from one individual with SCI. We then formulate a trajectory optimization problem to achieve a target pedaling cadence and crank power output, and simulate the FES-cycling activity to obtain optimal stimulation patterns.

## Material and Methods

**Subject:** One male participant with spinal cord injury (SCI) for over 14 years (J.S.M., 66 kg, 42 years, male) and a body mass index of 26.6 (±2.9) was included in the study. Written consent was obtained from the participant for the study procedures. The research was approved by a local Ethical Committee (CAAE: 30989620.6.0000.5149, Ethical Approval number 4.190.128) in agreement with the Declaration of Helsinki. The participant exhibited unrestricted joint mobility and did not present any lower limb fractures within the 12 months preceding the start of the protocol. The individual had previously undergone neuromuscular electrical stimulation (NMES) muscle strengthening and FES-cycling sessions for about four years prior to the commencement of the study.
**Experimental Setup:** A commercially recumbent tadpole trike (ICE Adventure, ICE, Falmouth, Cornwall, UK) was fixed on an instrumented cycling stationary trainer (KICKR Smart Trainer, Wahoo, Atlanta, Georgia, US). The original crankset was replaced by a power-meter crank (2INPOWER, Rotor Bike Components, Madrid, Spain) to quantify torque and pedaling power, with the resultant data transmitted via ANT+ communication protocol to a personal computer at a sample rate of 50 Hz, as detailed in [8]. Orthotic calf support (HASE Bikes, Waltrop, Germany) was fixed to the pedals to keep the ankle joint at 90° and restrict the leg movement to the sagittal plane. A personal electrical stimulator system with eight channels (detailed in [1]) was employed in the experiments. The stimulation interval and the stimulation parameters were adjusted via a mobile application.

**Measurements for parameter identification:** electrical stimulation pulse were applied to the following muscles: *rectus femoris* and *vastus lateralis* muscles; *vastus medialis*; and *hamstrings* muscle group (see [1] for electrode positioning). A warming-up phase preceded the isometric measurements phase, which are detailed below:

a) **Warm-up phase:** the volunteer was seated on the tricycle and hydrogel electrodes were attached to the legs over the motor points previously identified. The experimenter propels the participant legs passively (no FES) for 3 min at a maximum speed of 30 rpm.

b) **Static measurement:** the maximum evoked torque was first measured with the pedal at 90° (right crank arm in the vertical top position). Balanced biphasic pulses were delivered over the quadriceps muscle group with pulse amplitude of 50 mA, frequency of 35 Hz, and pulse duration varying from 500 to 1200 µs in 50 µs steps, for a stimulation duration of 500 ms. We then measured the range of motion (ROM) of the volunteer's knee while seated on the tricycle using a goniometer, which varied from 70° to 120°, and performed six isometric evoked contractions within this range, with knee angle varying in 10° steps, for both knee extension (quadriceps) and flexion (hamstrings). At this point, we set the pulse duration to ~40% of the value that elicited the strongest evoked contraction in the previous step (resulting in 800 µs). The experimenter moved the crank to the designated position, locked it in place by holding the pedal opposite the leg being measured, and provided a minimum of 5-second settling period before the stimulation began. A subsequently resting time of 2 s was then given. For each knee angle, two stimulations were performed before the crank was moved to the next position. This novel procedure allowed for the evaluation of evoked muscle contractions at different knee angles as an alternative to conventional dynamometers. Finally, 20 s of passive pedalling at 20 rpm was performed by the researcher.

**Musculoskeletal Model:** We adapted the musculoskeletal model of the lower limbs from the OpenSim core models Leg6Dof9Musc to be used with Bioptim [5]. The rotational movement was kept limited to the sagittal plane. The model has seven degrees of freedom (DOF) (gear, pedal right, pedal left, knee right, knee left, hip right, hip left). The ankles were kept at 90° and the feet were fixed to the

pedals. In addition, muscles in the model are *biceps femoris* long head for knee flexion, *vastus medialis* for knee extension and *rectus femoris* for both knee extension and hip flexion (Fig. 1).

<figure><img src="image_2.png" /><figcaption>Figure 1: FES-cycling musculoskeletal model</figcaption></figure>

The dynamics of our adapted muscle model captures the temporal influence of consecutive electrical pulses on the formation of the calcium-troponin complex and the regulation of muscle contraction leading to force development. The activation dynamics is given by:

$$ \frac{dc_N}{dt} = \frac{1}{\tau_c} \sum_{i}^{n(t)} \alpha_i R_i \exp\left(-\frac{t-t_i}{\tau_c}\right) - \frac{c_N(t)}{\tau_c}, \quad (1) $$

where T represents a constant pulse period and $n(t) = |t/T|$ is the largest integer $n$, such that $n \le t/T$, and $i = \max(1, n(t) - p)$ is within a limited window train of pulses containing the history of the last $p \in \mathbb{N}^*$ pulses. The activate or deactivate states of each pulse is given by $\alpha_i = 0,1$ during the $i^{th}$ pulse $t_i = iT$. The parameter $R_i$ accounts for the nonlinear summation of the cN transient upon the time response $\tau_c$ (ms) – a constant controlling the rise and decay of $c_N$, which is given by:

$$ R_i = \begin{cases} 1 & \text{if } i = 1, \\ 1 + (R_0 - 1) \exp\left(-\frac{t-t_i}{\tau_c}\right) & \text{if } i > 1, \end{cases} \quad (2) $$

with $R_0$ (unitless) representing the magnitude of enhancement in $c_N$ [9]. The activation dynamics is then related to the muscle force dynamics defined by:

$$ \frac{dF}{dt} = A(p_d, \theta) \frac{c_N(t)}{K_m + c_N(t)} - \frac{F(t)}{\tau_1 + \tau_2 \frac{c_N(t)}{K_m + c_N(t)}}, \quad (3) $$

where $K_m$ is the sensitivity parameter of strongly bound cross bridges to $c_N$, $\tau_1$ (ms) is the time constant of force decline in the absence of strongly bound cross-bridges and $\tau_2$ (ms) is the time constant of force decline due to friction between actin and myosin [9]. The resulting force $F$ (N) is parameterized by five parameters ($\tau_1$, $\tau_2$, $\tau_c$, $K_m$, and $R_0$) and one function ($A$) that accounts for the influence of the stimulation pulse duration and the knee angle. It is defined by:

$$ A(p_d(t), \theta) = a(\theta)(1 - \exp(-(p_d(t) - p_{d_0})/p_{d_t})), \quad (4) $$
where $a(\theta)$ represents a force scaling factor that is dependent on the knee joint angle $\theta$. The parameter $p_d(t)$ (µs) is the pulse duration, $p_{d_0}$ (µs) is a constant threshold offset characterizing the muscle sensibility to the pulse intensity, and $p_{d_t}$ (µs) describes how the scaling parameter varies with the increasing of pulse duration [9].

**Optimal Control Problem Formulation:** a multiple shooting trajectory optimization problem was formulated to both identify the model parameters (estimation) and to control the model in simulation. The resulting nonlinear programs are then solved using the IPOPT, with CasADi used for algorithmic differentiation. For a system whose dynamics is governed by an ordinary differential equation of the form $\dot{\mathbf{x}} = f(\mathbf{x}(t), \mathbf{u}(t))$, with $\mathbf{x}$ and $\mathbf{u}$ the state and control vectors, a generic trajectory optimization problem can be written as:

$$ \min_{\mathbf{x}(t),\mathbf{u}(t),\mathbf{p}} J(\mathbf{x}(t),\mathbf{u}(t),\mathbf{p}) \quad (5a) $$

$$ \text{s.t.} \quad \forall t, \quad \dot{\mathbf{x}} = f(\mathbf{x}(t), \mathbf{u}(t), \mathbf{p}) \quad (5b) $$

$$ \forall t, \quad \mathbf{g}(\mathbf{x}(t), \mathbf{u}(t)) \le 0 \quad (5c) $$

$$ \mathbf{x}(0) = \mathbf{x}_0 \quad (5d) $$

with $J$ the cost function, $\mathbf{p}$ a vector of model parameters and $\mathbf{x}_0$ the initial state. The inequality constraints are gathered in function $\mathbf{g}$. The objective is to find the control input $\mathbf{u}(t)$, the states $\mathbf{x}(t)$, and the parameters $\mathbf{p}$ that minimize the cost function while satisfying the dynamic and boundary constraints. The control input $\mathbf{u}(t)$ contains the pulse duration and the resistive torque. The latter needed only for the parameter identification (see **Model Personalization** steps 1 and 2) process; the states $\mathbf{x}(t)$ includes $c_N(t)$ and force $F(t)$ for each muscle, the joint positions $\mathbf{q}(t)$ and velocities $\dot{\mathbf{q}}(t)$, and pulse state $\alpha_i(t)$.

**Model Personalization:** the model is described with one function $A(p_d, \theta)$ and five parameters: $\tau_1, \tau_2, \tau_c, K_m$, and $R_0$. In Eq. 4, we defined $A$ as a function of $\theta$ and as function of pulse duration $p_d(t)$. Parameters $p_{d_0}$ and $p_{d_t}$ are the average values presented in [9]. The parameters $R_0$ and $\tau_c$ are fixed at 5 and 11 ms, respectively, as suggested in [9] for the quadriceps femoris muscles of individuals with SCI. Thus, only three parameters and one function remain to be identified for each subject and leg: $\tau_1, \tau_2, K_m$ and $a(\theta)$. Finally, the following steps were used to identify the model's parameters and obtain predicted torques.

**Step 1:** Identification of $\tau_1, \tau_2$, and $K_m$. For each muscle, we gathered six torque trajectories containing the first evoked contraction (see **Static measurements**) in an experimental torque vector $\tau^{ex}$. An optimal estimation problem was formulated to track these experimental torque data by optimizing $\tau_1, \tau_2, K_m$ and $A$ (Eq. 6), while satisfying the system's dynamics. Hence, a problem of the form of (Eq. 5) was solved with the following equation as a transcription of (Eq. 5a):

$$ \min_{\mathbf{x}(t),\mathbf{u}(t),\mathbf{p}} J = \frac{1}{N} \sum_{k=0}^{N-1} (\tau_k^{pr}(A, \tau_1, \tau_2, K_m) - \tau_k^{ex})^2, \quad (6) $$

where $\tau^{pr}$ is the torque trajectory predicted by the model and $N$ is the total number of shooting nodes (35 per second). To enforce isometric contractions during the optimization, the crank was locked by applying an opposite resistive torque through a null velocity constraint. We defined a constant $p_d(t) = 400 + p_{d_0}$ (µs) when stimulation was on, or $p_{d_0}$ otherwise. The control bounds were set according to the stimulation state (i.e., on or off) of the experimental data. As the parameters $\tau_1, \tau_2, K_m$ are not dependent on $\theta$, the resulting identified values for these parameters were kept fixed in the next steps.

**Step 2:** Determination of the $a(\theta)$ polynomial. In this step, we defined $a$ to be an optimization parameter and used six contractions from the dataset at $\theta = [70^\circ, 80^\circ, 90^\circ, 100^\circ, 110^\circ, 120^\circ]$ to identify $a(\theta)$. Similarly to the previous step, for each $\theta$ value, 6 separate problems were formulated to track the experimental torque data by optimizing $a$, while satisfying the system's dynamics. As a result, a problem of the form of Eq. 5 was solved with the following equation (Eq. 6) as a transcription of Eq. 5a:

$$ \min_{\mathbf{x}(t),\mathbf{u}(t),\mathbf{p}} J_{\theta} = \frac{1}{N} \sum_{k=0}^{N-1} (\tau_k^{pr}(\hat{\mathbf{a}}_{\theta}) - \tau_k^{ex})^2, \quad (6) $$

with $\hat{\mathbf{a}}_{\theta}$ the optimized value of $a$ at angle $\theta$. Finally, the values of $\hat{\mathbf{a}}_{\theta}$ were fitted with a 5th-degree polynomial, yielding an analytical expression for $a(\theta)$.

**FES-cycling optimal stimulation pattern:** to synthesize the optimal stimulation pattern, pulse duration was optimized as a piecewise constant control variable within 86-1200 µs at a constant frequency of 35 Hz. In Eq. 4, the

<figure><img src="image_3.png" /><figcaption>Figure 2: Solutions for the FES-cycling simulation: a) resulting stimulation pattern and the corresponding range of activation for each muscle; b) the crank angular velocity obtained while tracking ~4.4 rad/s. Due to the counter clockwise pedalling motion, the velocity appears to be negative; c) the power output produced at the crank axis of rotation while tracking 20 W.</figcaption></figure>
values of $p_{d_0}$ and $p_{d_t}$ were considered as the average values presented by [9]. The parameters for each muscle obtained in the previous steps were used. An optimal control problem was formulated to track a desired pedaling cadence (42 rpm or ~4.4 rad/s) by optimizing $p_d(t)$, while satisfying the system's dynamics (Eq. 5). A constant Coulomb friction torque ($-\mu_0 \cdot \text{sign}(\dot{q}_0)$, with $\mu_0 = 0.2$ Nm) and a viscous friction torque ($-B \cdot \dot{q}_0$, with $B = 1.0378$ Nm·s/rad) were applied at the crank axis to simulate the resistive environment opposing the bike's forward motion, where $\dot{q}_0$ is the crank velocity. The value of $\mu_0$ was set manually, while $B$ was computed based on the desired cadence and power. Hence, a problem of the form of Eq. 5 was solved with the following equation (Eq. 7) as a transcription of Eq. 5a:

$$ \min_{\mathbf{x}(t),\mathbf{u}(t)} \sum_{k=0}^{N-1} [\omega_1(\dot{q}_{0k}^{\text{sim}} - \dot{q}_{0k}^{\text{target}})^2 + \omega_2 \mathbf{p}_{d_k}^2 + \omega_3(c(\mathbf{u}(t)))_k^2], \quad (7) $$

where $\omega_1 = 5e2$, $\omega_2 = 1e3$, $\omega_3 = 1e3$, and $c$ is a function that penalizes co-activation of homologous muscles across both legs.

## Results

Figure 2.a presents the stimulation pattern, highlighting the intervals of the pedaling cycle in which each muscle is activated over a 3 s simulation, featuring coordinated agonist and antagonist activations. The resulting pedaling cadence was 4.09 ± 0.36 rad/s (Fig. 2.b), and the corresponding power output was 17.53 ± 7.87 (Fig. 2.c).

## Discussion

The proposed methodology demonstrates its ability in generating coordinated muscle activations to achieve specific target power output and pedaling cadence. As shown in Figure 2a, the stimulation pattern exhibits coherent phase-specific recruitment of agonist and antagonist muscles, closely aligned with the manually-tuned stimulation pattern obtained in [1]. Although the achieved cadence and power fall slightly below the target values defined, the results are encouraging given that it advances in the definition of the stimulation pattern relying entirely on muscle models that predict muscular force in response to electrical stimulation [8]. Our approach contrasts with prior FES-cycling simulation works [6], which are lacking muscle modeling.

Although significant nonlinearities and physiological constraints are introduced into the system, making the control problem more challenging, it offers the possibility of addressing realistic and personalized protocols needed for the generation of optimized muscle forces and joint kinematics. The ability of our model to produce near-target outputs under these constraints highlights the potential of personalized predictive simulations to support the design of FES protocols for individuals with SCI.

Despite the lack of experimental validation with the final stimulation pattern, this methodology can be further explored to improve FES-assisted activities performances, paving the way for more efficient rehabilitation strategies for people with SCI and other neurological conditions.

## Conclusions

This study presents a first step towards improving the realism of FES-driven simulations and showcases the possibility of using model-based control of FES-cycling, with possible application in other FES-assisted tasks.

## Acknowledgement

We thank the volunteer for his dedication and availability. This work was supported by Défi Clé Robotique centrée sur l'Humain funded by Région Occitanie, France.

## References

1. Coelho-Magalhães T., Fachin-Martins E., Silva A., Azevedo-Coste C., Resende-Martins H.. Development of a High-Power Capacity Open Source Electrical Stimulation System to Enhance Research into FES-Assisted Devices: Validation of FES Cycling. *Sensors*. 2022; no. 2:531.

2. van der Scheer, J.W., Goosey-Tolfrey, V.L., Valentino, S. E., Davis, G. M., & Ho, C. H.. Functional electrical stimulation cycling exercise after spinal cord injury: a systematic review of health and fitness-related outcomes. *J NeuroEngineering Rehabil.* 2021; 18(1),99

3. Maffiuletti, N.A. Physiological and methodological considerations for the use of neuromuscular electrical stimulation. *Eur J Appl Physiol* 110, 223–234 (2010).

4. Szeci J, Straube A, Fornusek C. A biomechanical cause of low power production during FES cycling of subjects with SCI. *J Neuroeng Rehabil.* 2014 Aug 16;11:123.

5. B. Michaud, F. Bailly, E. Charbonneau, A. Ceglia, L. Sanchez, and M. Begon, "Bioptim, a Python framework for musculoskeletal optimal control in biomechanics," *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, vol. 53, no.1, pp. 321–332, Jun. 2022.

6. de Sousa, A. C. C., Font-Llagunes, J. M.. Predictive Framework for Electrical Stimulation Cycling in Spinal Cord Injury. *IFAC-PapersOnLine*. 2024. Vol. 58, Issue 24, Pages 332-337.

7. T. Coelho-Magalhães, C. Azevedo-Coste and F. Bailly, "Numerical-Optimal-Control-Compliant Muscle Model for Electrically Evoked Contractions," in *IEEE Transactions on Medical Robotics and Bionics*, doi: 10.1109/TMRB.2025.3590453.

8. Schmoll Martin, Le Guillou Ronan, Fattal Charles, Coste Christine Azevedo. OIDA: An optimal interval detection algorithm for automatized determination of stimulation patterns for FES-Cycling in individuals with SCI. *J NeuroEngineering Rehabil.* 2022;19:39.

9. J. Ding, Li-Wei Chou, T. M. Kesar, S. C. K. Lee, T. E. Johnston, and A. S. Wexler, and S. A. Binder-Macleod, "Mathematical model that predicts the force-intensity and force-frequency relationships after spinal cord injuries," *Muscle Nerve*, vol. 36, no. 2, pp. 214–222, Aug. 2007.

## Author's Address

T. Coelho-Magalhães; tiago.coelho-magalhaes@inria.fr