# Moving Horizon Estimation of Human Kinematics and Muscle Forces

Amedeo Ceglia, François Bailly, Mickael Begon

► To cite this version:

Amedeo Ceglia, François Bailly, Mickael Begon. Moving Horizon Estimation of Human Kinematics and Muscle Forces. IEEE Robotics and Automation Letters, 2023, 8 (8), pp.5212 - 5219. 10.1109/LRA.2023.3291921. [hal-04151820](https://HAL Id: hal-04151820)

Submitted on 21 Aug 2023

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
# Moving horizon estimation of human kinematics and muscle forces

Amedeo Ceglia<fnref n="1" />, Francois Bailly<fnref n="2" />, Mickael Begon<fnref n="3" />

**Abstract**—Human-robot interaction based on real-time kinematics or electromyography (EMG) feedback improves rehabilitation using assist-as-needed strategies. Muscle forces are expected to provide even more comprehensive information than EMG to control these assistive rehabilitation devices. Measuring in vivo muscle force is challenging, leading to the development of numerical methods to estimate them. Due to their high computational cost, forward dynamics-based optimization algorithms were not viable for real-time estimation until recently. To achieve muscle forces estimation in real time, a moving horizon estimator (MHE) algorithm was used to track experimental biosignals. Two participants were equipped with EMG sensors and skin markers that were streamed in real time and used as targets for the MHE. The upper-limb musculoskeletal (MSK) model was composed of 10 degrees-of-freedom actuated by 31 muscles. The MHE relies on a series of overlapping trajectory optimization subproblems of which the following parameters have been adjusted: the fixed duration and the frame to export. We based this adjustment on the estimation delay, the muscle saturation, the joint kinematic mean power frequency, and errors to experimental data. Our algorithm provided consistent estimates of muscle forces and kinematics with visual feedback at 30 Hz with a 110 ms delay. This method is promising to guide rehabilitation and enrich assistive device control laws with personalized force estimations.

**Index Terms**—Modeling and Simulating Humans, Physical Human-Robot Interaction, Sensor-based Control.

## I. INTRODUCTION

REAL-TIME biofeedback providing biomechanical metrics is highly relevant in rehabilitation [1] or to enhance control strategies of assistive robots [2]. In particular, feedback on electromyography (EMG), kinetics, and joint angles has successfully improved post-stroke patient rehabilitation and has proven to be useful in human-robot interaction [3]. This feedback must be provided at least at 13 Hz [4] and with a maximum delay of 125 ms to be relevant [5]. Inverse kinematics algorithms, such as the extended Kalman filter (EKF) [6], estimate joint angles from marker positions and provide real-time kinematics feedback [7]. When it comes to kinetic data, muscle forces provide more information than

EMG as they depend on both muscle activations and force-length-velocity relationships [8]. However, measuring muscle forces in vivo can be challenging due to the invasive nature of some measurement techniques and instrumentation limitations. To address this issue, researchers have focused on developing numerical models to estimate muscle forces.

Inverse methods, especially static optimization, are the most used algorithms for real-time muscle force estimation, due to their computational efficiency [9]. However, they rely on joint torques calculated by inverse dynamics, which is affected by the noise amplification due to the double differentiation of kinematic data [10]. Moreover, static optimization is often based on minimizing least-square muscle activations without considering control time histories. It could likely lead to non-physiological activations, especially for patients with MSK disorder [11].

Forward methods, which require finding optimal trajectories for the state and control variables satisfying the MSK dynamics, could lead to more reliable estimation. Additionally, incorporating tracking terms on experimental data, such as EMG and skin markers, has been shown to improve estimation accuracy [12], particularly for patients with pathological conditions [13]. These methods address the noise introduced by double differentiation and provide more physiological activations. However, their computational cost and challenging convergence may compromise their use in real time applications. To address computational limitations, a combination of forward and inverse methods can be employed. This approach takes advantage of both methods to improve computation time and increase dynamic consistency. However, some methods rely on noisy data on joint moments, such as the CEINMS toolkit [14], while others do not consider the control time history, like the forward static optimization approach [15]. Receding horizon methods have demonstrated promising outcomes when leveraging forward dynamics [16], [17]. In [17], we proposed a real-time muscle force estimation algorithm based on: optimal control problem (OCP), forward dynamics, moving horizon estimator (MHE), and simulated data. An originally intractable OCP was split into a series of smaller subproblems solved at high frequency (15-30 Hz). We showed that forward methods are suitable for real-time applications, providing muscle forces and joint kinematics satisfying motion dynamics. Unfortunately, it was only tested on a simplified model (4 degrees-of-freedom (DoFs), 10 muscles) using simulated data. This study builds upon prior work and makes a three-fold contribution:

1) to stream experimental data to live feed our algorithm, and to provide real-time feedback of estimated data (joint kinematics and muscle forces) on two participants

Manuscript received: January 31, 2023; Revised May 24, 2023; Accepted June 23, 2023.

This paper was recommended for publication by Editor Aniket Bera upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the CREATE OPSIDIAN program and the discovery grant RGPIN/04978-2019

<footnote n="1">¹A. Ceglia is with the Institute of Biomedical Engineering, University of Montreal, Montreal, QC H3T 1J4 Canada amedeo.ceglia@umontreal.ca</footnote>

<footnote n="2">²F. Bailly is with INRIA, Univ Montpellier, Montpellier, France</footnote>

<footnote n="3">³M. Begon is with the School of Kinesiology and Human Kinetics, University of Montreal, Montreal, QC H3T 1J4 Canada.</footnote>
with distinct morphologies (e.g., male and female);

2) to show the applicability of real-time MHE using experimental biosignals and a more refined MSK model;

3) to develop a method to determine the best receding horizon parameters (problem duration, and at which time frame to keep a sample of the optimal trajectory).

## II. METHOD

### A. Musculoskeletal model

The MSK model derives from the upper-limb model developed by Wu & al. [18]. It consists of five segments (thorax, clavicle, scapula, arm, and lower arm) and 10 DoFs (two sternoclavicular, three acromioclavicular, and three glenohumeral rotations, elbow flexion, and pronation/supination). The model is actuated by 31 Hill-type muscles [19] around the shoulder and the elbow. Biceps brachii and triceps brachii muscle groups were added to the original model [18] based on Holzbaur's model [20]. Pure additional joint torques are included in the optimization to assist muscle actions on the sternoclavicular, acromioclavicular, and elbow (supination-pronation) joints [21].

### B. Muscle force estimation

The Hill model, which describes muscles by a contractile element and a passive element in parallel, was used for each muscle unit. To speed up computations, the tendons were considered rigid [22]. As the movements of interest were slow (around 70°/s), muscles were directly controlled by their activation (neglecting neural excitation dynamics with its electromechanical delay [23]) while trying to match online recorded EMG signals. For a muscle *m*, the force ($f^m$) computed from activation (a) and muscle kinematics (normalized muscle fiber length $\tilde{l}_m$ and velocity $\tilde{v}_m$) was expressed as:

$$f^m = f_o^m (af^{act}(\tilde{l}_m)f^v(\tilde{v}_m) + f^{pas}(\tilde{l}_m)) \cos \alpha, \quad (1)$$

where $f_o^m$ is the muscle isometric force, $f^{act}(\tilde{l}_m)$, $f^v(\tilde{v}_m)$ and $f^{pas}(\tilde{l}_m)$ are the active force-length relationship, force-velocity relationship and the passive force-length relationship, respectively; $\alpha$ is the muscle pennation angle [24].

### C. Problem formulation

Our framework consists in solving a series of OCPs of fixed duration and number of frames ($t_{mhe}$ and $n_{mhe}$, respectively), implemented as constrained nonlinear problems. Each sub-problem aims to find estimates of the control and of the state trajectories (respectively $\hat{u}$ and $\hat{x}$) that minimize the tracking error of experimental data (denoted by $*$), while satisfying the system dynamics and constraints according to the following formulation:

$$(\hat{\mathbf{x}}, \hat{\mathbf{u}}) = \arg \min_{\mathbf{u}} \mathcal{J}(\mathbf{x}, \mathbf{u}) \quad (2a)$$

$$\text{s.t.} \quad \forall t, \quad \dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{u}) \qquad \text{Dynamics constraints} \quad (2b)$$

$$\forall t, \mathbf{u} \in \mathcal{U}, \mathbf{x} \in \mathcal{X} \qquad \text{Path constraints} \quad (2c)$$

Equation (2a) is the running cost to be minimized (developed in (3a-3c)). This optimization problem is subject to multi-body dynamics driven by muscles and joint torques (Eq. 2b), as well as path constraints in the domains of $\mathbf{u}$ and $\mathbf{x}$ ($\mathcal{U}$

and $\mathcal{X}$, respectively) that include muscle activation ranging from 0 to 1 and joint kinematics within their physiological ranges of motion. A state constraint is enforced to dictate the solver's initial state $\mathcal{J}$ is the cost function to be minimized and expressed as:

$$\mathcal{J} = \sum_{n=1}^{n_{mhe}} (\omega_{m^*} \|\mathbf{m}_n(\mathbf{q}) - \mathbf{m}_n^*\|^2 + \omega_{a^*} \|\mathbf{a}_n - \mathbf{a}_n^*\|^2) \quad (3a)$$

$$+ \int_{t}^{t+t_{mhe}} (\omega_{\tau} \|\boldsymbol{\tau}\|^2 + \omega_{a} \|\boldsymbol{a}\|^2 + \omega_{q} \|\boldsymbol{q}\|^2 + \omega_{\dot{q}} \|\boldsymbol{\dot{q}}\|^2) dt \quad (3b)$$

$$+ \sum_{n=1}^{n_{mhe}} (\omega_{q^{\text{prev}}} \|\mathbf{q}_n - \mathbf{q}_n^{\text{prev}}\|^2 + \omega_{\dot{q}^{\text{prev}}} \|\dot{\mathbf{q}}_n - \dot{\mathbf{q}}_n^{\text{prev}}\|^2) \quad (3c)$$

here $\mathbf{m}$ is a vector of marker positions (obtained from $\mathbf{q}$ by forward kinematics), $\mathbf{q}$ and $\dot{\mathbf{q}}$ are joint angles and velocities respectively, with $\mathbf{x} = (\mathbf{q}, \dot{\mathbf{q}})$ and $\boldsymbol{\tau}$ is a vector of pure additional joint torques. Equation (3a) stands for the tracking of marker positions and EMG signals, weighted by $\omega_{m^*}$ and $\omega_{a^*}$, respectively. Experimental data are tracked forward in time allowing a more physiological estimation [25]. Equation (3b) represents controls minimization (additional joint torques and muscle activations), weighted by $\omega_{\tau}$ and $\omega_{a}$, respectively, and regularization terms on states $\mathbf{q}$ and $\dot{\mathbf{q}}$, weighted by $\omega_{q}$, and $\omega_{\dot{q}}$, respectively.

Equation (3c) is a correction-based term added for tracking estimated states from the previously solved problem, weighted by $\omega_{q^{\text{prev}}}$ and $\omega_{\dot{q}^{\text{prev}}}$, respectively. These terms act as a proportional derivative controller, to reduce steady-state error and noise [26]. The weights on these terms can be interpreted as the gain of the controller.

Weights ($\omega_{m^*}$, $\omega_{a^*}$, $\omega_{\tau}$, $\omega_{a}$, $\omega_{q}$, $\omega_{\dot{q}}$, $\omega_{q^{\text{prev}}}$, and $\omega_{\dot{q}^{\text{prev}}}$) were identified offline by finding the best trade-off between marker tracking and the quality of the EMG using recorded data from the first participant. Each time a sub-problem is solved, a new one is set up with the most recent measurements used as targets. In addition, as the sub-problems overlap, previously estimated points are used to provide an initial guess close to the optimal solution of the current problem. To build the initial guess, the oldest estimated frame is removed from the previous solution, and the most recent is repeated as many times as needed to match the current problem size (Fig. 1). This overlap requires selecting which sample to keep as the optimal solution, expressed as a percentage of $n_{mhe}$, referred to as frame$_x$ ($x$ represents the percentage of $n_{mhe}$). The final delay is dependent on both the solving delay and which frame$_x$ is exported (Fig. 2). To catch up with upcoming experimental data, the estimation horizon must move forward in time, by a number of samples denoted $s$, named as shift size, computed through the division between sampling and solving frequencies (Fig. 1). Extensive tests were conducted on the receding horizon parameters offline using the recorded data for the first participant to determine the optimal trade-off between performance and accuracy among:

* $t_{mhe}$: large values lead to higher accuracy but lower solving frequencies;

* frame$_x$: the closer the frame is to the end (100% of $n_{mhe}$), the less it is filtered by the algorithm.
These receding horizon parameters were fine tuned while varying $t_{mhe}$ from 4 ms to 120 ms and $frame_x$ amongst 0%, 25%, 50%, 75%, and 100% of $n_{mhe}$, for three movements (shoulder abduction, flexion, and one that mimics hand cycling (see video in supplementary material)). Six criteria (italicized below) were used to determine the best parameters. As we targeted a delay below 125 ms [5], the *overall estimation delay* was computed as the sum of the solving duration and the delay introduced by the chosen kept $frame_x$. In addition, as the movements were not executed at a maximal contraction level, muscle activations should not reach saturation. Therefore, we evaluated the total number of frames where muscles were saturated (activation > 95% of the maximal voluntary contraction (MVC)) for each set of receding horizon parameters. To avoid non-physiologic results on the kinematics, the *angular velocity's mean power frequency* was computed for all joints, representing their inner variation [27]. Then, the *root-mean-square errors (RMSE)* were calculated between estimated and experimental data (markers, EMG). Finally, the inverse dynamics based on experimental markers data was computed. The *root-mean-square difference (RMSD)* was calculated between the computed inverse dynamics and the estimated torques (sum of muscle-induced torques and additional torques).

Results are shown in Section III-A. The problem was implemented using Bioptim [28], a Python library dedicated to solving biomechanical OCP, leveraging the MHE pre-implemented class. Non-linear programs were solved using Acados [29], with its real-time iteration scheme [30] to improve the robustness and computational cost of the convergence of the subproblems. In particular, the OSQP [31] solver was used to solve the quadratic program. It allows for stricter boundary constraints in contrast to HPIPM [32] used in our previous study [17]. To avoid indefinite Hessian in the OSQP, a regularization term was added through the Levenberg-Marquardt coefficient [33]. Then, a damped Newton method was used by adjusting the step length to 0.9 to improve convergence rate [34]. Also, the Jacobian matrix was reused within Newton iterations to reduce computational cost [35]. Since the acquisition rate was set to 100 Hz, the number of available point frames for each sub-problem (e.g., 10 points over 0.1 s) was insufficient to solve the problem accurately. To address this, the experimental data were linearly interpolated, resulting in twice the number of points (e.g., 20 points for 0.1 s). The code is available on GitHub.

## D. Participants and instrumentation

Two healthy right-handed participants were involved after signing an informed consent: one male (23 years old, 72 kg, 1.85 m) and one female (29 years old, 70 kg, 1.68 m). This experiment was approved by the ethics committee of the University of Montreal (#14-122-D). A 12-camera motion capture system (Vicon, Oxford, UK) was used to record 16 skin marker trajectories at 100 Hz. Markers were placed on the trunk and right upper limb at predefined locations to reduce soft tissue artifacts [36] (Fig. 3). Muscle activities were recorded at 2 kHz using 10 surface EMG sensors (Trigno

<figure><img src="image_2.png" /><figcaption>Figure 1. Series of two consecutive subproblems at time $t-1$ and $t$ (grey and black boxes, respectively), shifted by a *shift size*. The red curve represents the ground-truth trajectory for one joint angle. The solid grey and black lines correspond to the previous and current solutions, respectively, which aim to estimate the true trajectory. The dashed black line represents the initial guess of the subproblem at $t$, built from the solution at $t-1$. The grey circle denotes the initial state constraint between the two consecutive subproblems.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Representation of a series of MHE subproblems (from $t-4$ to $t$) of five frames size with a one frame shift size. Experimental data are used as targets. The red pie charts represent the exported frame for each sub-problem corresponding to the estimation for the current sample.</figcaption></figure>

EMG Wireless System, Delsys, USA) put on the Pectoralis (major), Biceps brachii (long head), Triceps brachii (long heads), Trapezius (upper, medial, and lower heads), Deltoid (anterior, medial, posterior heads), and Latissimus dorsi, after shaving and cleaning the skin. Data were synchronized using Nexus 2.8 (Vicon, Oxford, UK). Ten MVC were collected following the recommendations of [37]. Offline EMG signals were band-pass filtered (10-425 Hz), rectified, and lowpass filtered (5 Hz) using a fourth-order Butterworth filter. Then, maximal activations were obtained by computing the mean of the 2000 highest values (corresponding to 1 s) for each muscle during all 10 MVC trials [37]. A static pose was used to scale the MSK model using OpenSim [38], which was then converted into a Biorbd model format [39], as Bioptim MSK dynamics relies on Biorbd [39]. Finally, the participants were asked to perform two arm elevations in the frontal and sagittal planes and movements mimicking hand cycling in the sagittal plane. For further validation, their were asked to perform motions while holding a 2 kg dumbbell which result in different muscle activity levels.

## E. Real-time software implementation

We developed a fast (C++ core) and open-source Python library to achieve real time muscle force estimation and biofeedback using live data. Two scripts were implemented: (1) sensor data streaming and processing, and (2) muscle force and kinematics estimation and visual biofeedback. To enhance performance, the two scripts were executed on two computers connected through TCP/IP connection, with the Python multiprocessing package used for parallelization. The software architecture is depicted in Fig. 4.
<figure><img src="image_4.png" /><figcaption>Figure 3. First participant equipped with 10 EMG sensors and 16 skin markers on the right upper limb.</figcaption></figure>

1) *Sensor data streaming and processing (1st computer; three tasks):* During Task #1a, skin markers and EMG signals were streamed in real time from Nexus using the Vicon DataStream SDK Python interface [40]. Concurrently, joint kinematics were computed from the marker 3D positions using EKF, and EMG envelopes were generated using a sliding window (2000 frames, 1 s), band-pass filtering (10-425 Hz), rectification, and moving average smoothing (200 frames) during Tasks #1b. In Task #1c, a TCP/IP server was started to transmit EMG envelopes, joint angles, and marker trajectories to the second computer on demand.

2) *Muscle force and kinematics estimation and biofeedback (2nd computer; two tasks in parallel):* Task #2a involved implementing the optimal estimator described in Section II-C and updating tracked data at each iteration via the TCP/IP server. EKF-based joint angles were used as an initial guess, and estimated data were made available for transmission to a visual application at the end of each iteration. In Task #2b, real-time biofeedback on optimal joint kinematics and muscle forces was provided through visualization (Fig.7).

## F. Data analysis

Once the receding horizon parameters were tuned using data from the first participant (Section II-C), the RMSE of the marker trajectories and the EMG, and the RMSD of the inverse dynamics were used to evaluate our results. To assess the robustness of our method, we used the same parameters without any changes to run the estimation on a second participant. As one of the objectives of this study was to demonstrate the real-time capabilities of our method we evaluated the time performance of our estimator. First, the solving frequency was computed for each iteration. Second, the total process latency between measured data and final biofeedback visualization was expressed as the sum of the following delays: $d_1$) Vicon SDK latency using the Python interface; $d_2$) data processing duration (Section II-E); $d_3$) latency of the TCP/IP connection; $d_4$) estimator solving duration (Section II-C); $d_5$) visualization refresh delay. Latency was computed for each solver iteration; we reported the mean value with standard deviation ($m \pm sd$) to evaluate the program's real-time capabilities.

# III. RESULTS

## A. Receding horizon parameters

The sum of non-compressible delays ($d_{1,2,3,5}$) was 50.22±10.92 ms (Tab. I). For a total delay under 125 ms [5] $d_4$

<table><caption>Table I
MEANS OF COMPUTATIONAL DELAY INTRODUCED BY EACH TASK.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Tasks</th><th colspan="3">Delay (ms)</th></tr><tr><th>name</th><th>mean</th><th>SD</th></tr></thead><tbody><tr><td rowspan="3">Computer 1</td><td>Nexus</td><td>d1</td><td>6.6</td><td>1</td></tr><tr><td>Processing</td><td>d2</td><td>9.5</td><td>1.07</td></tr><tr><td>TCP/IP</td><td>d3</td><td>6.12</td><td>1.25</td></tr><tr><td rowspan="3">Computer 2</td><td>Estimator</td><td>d4</td><td>33.9</td><td>1.1</td></tr><tr><td>Frame shown</td><td>d5</td><td>25.4</td><td>0.8</td></tr><tr><td>Visualization</td><td>d6</td><td>28</td><td>7.6</td></tr></tbody><tfoot><tr><td colspan="2">Total</td><td></td><td>109.6</td><td>12.2</td></tr></tfoot></table>

must be under 60 ms (125 - (50.22 + 10.92) = 63.86 ms, with one standard deviation conservatively included). Receding horizon parameters inducing delays longer than 60 ms (contained in the red shaded area, Fig. 5-a) were discarded. Fig. 5-b shows that every frame100% and frame75% with $t_{mhe}$ = 40 ms led to more than 4% of muscle activations saturation. As the saturation needs to be as low as possible, receding horizon parameters inducing saturation higher than 4% (contained in the shaded red area, Fig. 5-b) were discarded. Parameters discarded in Fig.5-a appear in grey in Fig.5-b. Similarly, parameters discarded in Fig.5-a and b appear in grey in Fig. 5-c. Fig. 5-c shows that the shorter the horizon length, the higher the mean power frequency of joint angular velocities. To remain physiologic the mean power frequency must remain as low as possible given the slow movements executed by the participants. It was observed that choosing an 80-ms window length resulted in kinematics as smooth as those obtained with the Kalman filter. No significant difference in errors with experimental data and difference in torques with inverse dynamics are observed (Fig.,5-d-f). Therefore, the set of parameters that introduced the smallest mean power frequency among the ones not previously discarded (frame75% and $t_{mhe}$ = 80 ms) was chosen, resulting in a delay under 60 ms, error lies within the midpoint of the maximum and minimum range, and good performance (30 Hz, which is twice the real-time threshold of 13 Hz [4]).

## B. Real-time performances

Data were streamed and processed (Nexus software included) on a Windows 7 computer (Intel Xeon® X5680 @ 3.33 GHz 6-core processor). The MHE and biofeedback were carried on eight cores of a second computer (AMD Ryzen® 9 5950x 16-core processor) running Ubuntu 20.04 LTS. The program estimated and displayed the 31 muscle forces, and the 10 DoFs 3D kinematics at 30 Hz with a total delay of 110 ms (Section II-E). The delays introduced by each step are summarized in Tab. I. Main delays come from the second computer for the estimation and visualization tasks (80%). Data streaming and processing on the first computer correspond to 20% of the total latency.

## C. Muscle force and kinematics estimation

As described in Section III-A, receding horizon parameters have been adjusted (frame75% - $t_{mhe}$ = 80 ms). Although the algorithm parameters were tuned on data from the first participant, the algorithm converged for the second participant without any change. Results for both participants are
<figure><img src="image_5.png" /><figcaption>Figure 4. Overview of the tasks performed on two computers. The first one (in blue) is for data streaming and processing, and the other one/(in green) is for muscle forces and joint kinematics optimal estimation and visualization.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Average performances on all trials for the first participant (n=6) as a function of mean problem-solving duration for each set of receding horizon parameters. Exported frames (frame_x) are represented by pie charts from 0% (●) to 100% (○) while colors represent length of subproblems (t_mhe) from 40 ms to 120 ms. The figure shows: a) the delay introduced by the optimization; b) The per-trial average percentage of frames where muscle activations are saturated; c) The mean power frequency on estimated joint angular velocities; d) and e) The RMSE between estimated and experimental data, marker positions and EMG respectively; f) The RMSD between estimated torques and those computed from inverse dynamics; Discarded parameters during the parameters selection process (inside red areas of a and b) appear in grey in the subsequent plots (respectively b and c). A red circle surrounds the selected best parameters (frame_75% - 80 ms).</figcaption></figure>

presented separately to demonstrate the ability of the method to generalize without any extra tuning. After applying forward kinematics to estimate the marker positions from joint angles, the mean error across all trials for all markers was 24.8 ± 6.7 mm and 20.3 ± 7.5 mm for the first and second participant, respectively, when compared to the experimental positions. The RMSE between estimated activations and recorded EMG was 7.39 ± 1.95% and 3.4 ± 5.4% for the first and second participant, respectively. Estimated activations, corresponding EMG (if applicable), and muscle forces are shown in Fig. 6 for the first participant, for an abduction movement, with and without a 2-kg dumbbell. As expected, higher activations and forces were estimated when holding the dumbbell, even for non-tracked muscles. This shows the ability of the solver to consider experimental EMG but also to generate a physiological prediction for untracked muscles. Furthermore, the activation level of the infraspinatus and supraspinatus were of the same order as in [41] (Fig. 6, data were extracted from graphics using Engauge Digitalizer [42]).

The subscapularis shows the same activation pattern but with an activation almost three times higher. A loss of activation is observed (Fig. 6) for the supraspinatus around 90° of abduction which differs from [41]. The same observations were made on the seconds participant.

## IV. DISCUSSION

We presented an end-to-end method to estimate and visualize muscle forces and 3D kinematics applied to a 10-DoF upper limb model actuated by 31 muscles in real time. To our knowledge, this is the first time that a trajectory optimization algorithm coupled with MHE has been used to achieve such estimation in experimental conditions. We highlighted the robustness of our method by successfully applying it to two participants with different morphology's, using parameters based solely on the first participant. The results remained consistent for both participants, and notably, improved even further for the second participant. We developed a framework to stream and process experimental data in real time (100 Hz)
<figure><img src="image_7.png" /><figcaption>Figure 6. Muscle activations normalized by MVC (with corresponding recorded EMG, dashed lines) and forces normalized by body weight (BW) are shown for an abduction motion for the first participant. EMG from infraspinatus, subscapularis, and supraspinatus were not recorded in this study, therefore data from Wickham & al. [41] (red shaded areas and dashed lines) are used to represent the condition with 2 kg dumbbell. Data for abduction greater than 90° were removed from [41] to match with our abduction range of motion. Both conditions are represented without weight (blue) and with 2 kg dumbbell (orange/red). Only eight representative muscles are shown in the figure. The estimated arm elevation is duplicated under each column to facilitate the reading (with EKF estimate for comparison, dash-dotted lines).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Experimental setup with real-time visual biofeedback. The screen was put in front of the participant. An enlargement was made on the screen (left) and on the display of one muscle force (bottom).</figcaption></figure>

and used numerical optimization to fuse EMG and kinematics data for predicting kinematics and muscle forces. Biofeedback was provided at 30 Hz with an average delay of 110 ms. This method might directly be applied to the control of a patient force-informed assistive device, with lower delay, since the visualizer (25% of the delay) might not be required in this context.

Since our previous work, our method was adapted to a more refined and complex model [17]. Our previous model was insufficient in capturing the complexities of a patient's rehabilitation task, such as hand cycling, where an accurate representation of shoulder kinematics necessitates the inclusion of clavicle and scapula movements. Additionally, the limited number of muscles in the model hindered our understanding of muscle actions in the task. Real-time performance was maintained while transitioning from simulation-only to full experimental conditions. This was achieved through software development and formulation improvements of the MHE. First, we developed an open-source Python package (biosiglive [43]) designed for efficient (C++ based) real-time data streaming.

Second, a real-time iteration scheme was used inside Acados, and algorithm adjustments were chosen to allow robust and fast solving.

We have highlighted that receding horizon parameters influence the estimation results, leading to, for example, high-frequency kinematics trajectories or activation saturation. Fig. 5 shows that using frame<sub>100%</sub> led to inconsistent saturated muscle activations (> 4%). This non-physiological result suggests that, when keeping frame<sub>100%</sub> in the estimated trajectory (subproblem(t-4) in Fig. 2), the corresponding movement was not sufficiently filtered, as it went only once through the MHE. Fig. 5 also reveals that the saturation of muscle activations stays around 1% when the exported frame<sub>x</sub> (x < 100%) was filtered more than once. The frame<sub>75%</sub> (at least two times filtered) could therefore be chosen as a good compromise leading to physiological estimates and low optimization delay (Fig. 5-a). In addition, the mean power frequency of joint angular velocities increased as the horizon length decreased. As high mean power frequency could lead to non-physiological joint acceleration, the parameters which introduced the smallest mean power frequency were kept. Interestingly, the combined analysis of optimization delay, muscle activation saturation and mean joint velocity's power frequency led to a unique choice for receding horizon parameters: frame<sub>75%</sub> with t<sub>mhe</sub> = 80 ms. Concerning experimental data tracking, receding horizon parameters did not seem to have a significant impact overall.

We have shown that, once the appropriate parameters were found, the results were consistent with other studies [18], [41] for the EMG and muscle forces patterns (Fig. 6), in addition to being dynamically consistent. Our method can estimate the movement kinematics with an average for both participants around 22 mm error on the markers, which is below the uncertainty due to soft tissue artifacts [36]. Despite the model's muscle parameters being uncalibrated and the absence of passive joint torque, the estimator still managed to consistently
estimate muscle contributions. First, as expected, the estimation provided higher muscle activations with a 2-kg dumbbell than without for both participants. Then, the RMSE for both participants around 5% when compared to EMG signals was low, particularly considering that a small change of the sensor position can induce up to 35% variation in signal amplitude [44]. The muscle activation and force level aligned with data from Wickham & al. [41] for untracked muscles such as the infraspinatus or supraspinatus. This consistency illustrates the estimator's capability to harness the system's dynamics and provide relevant estimates of concealed or unmeasured quantities. Additionally, the levels and patterns of muscle force were consistent with those found in [18]; for instance, for the deltoid (medial) (40% at 60°) or the infraspinatus (25% at 60°). Finally, the discrepancy between the superior trapezius's activation level (50% for a 2-kg dumbbell) and its force contribution, which is nearly 10 times lower (6%), underscores the relevance of opting for muscle force feedback over EMG. Our results must be put into perspective with inverse methods that allowed resolution at 120 Hz with a 44-DoF model actuated by 300 muscles [9]. The resolution is four times faster but relies on inverse dynamics that could lead to nonphysiological results. Previous study [45] that used hybrid methods on a 23-DoF model actuated by 34 muscles was able to estimate the joint contact force and display feedback in a total delay of 115 ms which is similar to ours. However, the estimation part, from inverse kinematics to muscle forces estimation, took around 57 ms (around 34 ms in our method, almost two times faster). We have shown that using forward dynamics optimization based can challenge hybrid methods while giving more reliable dynamics estimation. The visualization delay is the second-most resource-demanding feature in our study (28 ms - 25%) and could be improved for better performances. That said, our total latency (≈110 ms) is low enough for efficient visual biofeedback (≤125 ms).

In the continuity of this work, we will adapt our framework to a patient-specific assistive control of a hand-bike device dedicated to rehabilitation. The optimal control program implementation can be customized to add other measurements such as reaction forces. Moreover, as estimated values are available for each subproblem (30 Hz), they could be used to control any system requiring these values as input to interact with a patient (functional electrical stimulation system, for instance, [46]). In this study, we decided to present a method that provides muscle forces while acknowledging that additional parameters such as joint reaction force or moment arm can be estimated based on the system's requirements during patient interaction.

This study has limitations that should be considered. First, to improve convergence, the experimental data were interpolated from 100 to 200 Hz. The size of the MHE windows being small and the motion slow, interpolation effects should be low. However, the data should be acquired at a higher frequency in future work to avoid interpolation. The muscle parameters were not calibrated based on the participants, due to the yet unresolved challenge of calibrating muscles actuating more than one joint [14]. So, the amplitude of the estimated activation and force might lack accuracy. Nevertheless, we

evaluated the muscle force and activation time histories to show the relevance of the estimated muscle contributions. Only the subscapularis shows an activation almost three times higher than measured activity [41], which could be explained by the lack of scapulothoracic kinematics constraints or the non-calibrated muscle parameters (e.g., isometric force or optimal length). Also, we observed a loss of muscle contribution for the supraspinatus around 90° of abduction, which could be attributed to non physiological moment arm. Our model lacks features like joint capsules and the scapulothoracic pseudo-joint, which naturally enforce internal constraints. The addition of pure joint torques prevents the drawback to induce muscle contributions overestimation. However, additional glenohumeral torques were set to zero. As obtaining non-invasive in vivo measurements of muscle forces presents challenges, we were unable to directly validate our results. However, indirect validation was provided by comparing muscle activations and forces to literature data, even for non-tracked muscles (Section III-C). To further validate our method, we aim to track only 60-80% of recorded EMG and compare the estimated muscle activation with the experimental non-tracked EMG. Additionally, we tested our algorithm on only three motions and two participants. While the results are encouraging, further validation of our method's generality is necessary by applying it, in future work, to different contexts with more participants.

## V. CONCLUSIONS

We showed that trajectory optimization coupled with a MHE formulation can estimate muscle forces in real-time from experimental EMG-marker data. Real-time capabilities are good enough to be used in a clinical context. Biofeedback of estimated value was provided at 30 Hz with a 110 ms delay, on a complex model in experimental conditions. The estimated joint kinematics and muscle forces were consistent with the experimental motion and the literature. We are enthusiast that our algorithm provide consistent estimation for two different participants despite fine-tuning the parameters exclusively based on the first participant. Such a method could be used in clinical or industrial applications to give relevant real-time input data to assistive devices, leading to biomechanically informed interactions.

## REFERENCES

[1] J. R. Franz, M. Maletis, and R. Kram, "Real-time feedback enhances forward propulsion during walking in old adults," *Clinical biomechanics*, vol. 29, no. 1, pp. 68-74, 2014.

[2] J. M. Prendergast, S. Balvert, T. Driessen, A. Seth, and L. Peternel, "Biomechanics aware collaborative robot system for delivery of safe physical therapy in shoulder rehabilitation," *IEEE Robotics and Automation Letters*, vol. 6, no. 4, pp. 7177-7184, 2021.

[3] T. Teramae, T. Noda, and J. Morimoto, "Emg-based model predictive control for physical human-robot interaction: Application for assist-as-needed control," *IEEE Robotics and Automation Letters*, vol. 3, no. 1, pp. 210-217, 2017.

[4] O. A. Kannape and O. Blanke, "Self in motion: sensorimotor and cognitive mechanisms in gait agency," *Journal of neurophysiology*, vol. 110, no. 8, pp. 1837-1847, 2013.

[5] T. Waltemate, I. Senna, F. Hülsmann, M. Rohde, S. Kopp, M. Ernst, and M. Botsch, "The impact of latency on perceptual judgments and motor performance in closed-loop interaction in virtual reality," in *Proceedings of the 22nd ACM conference on virtual reality software and technology*, 2016, pp. 27-35.
[6] K. Halvorsen, T. Söderström, V. Stokes, and H. Lanshammar, "Using an extended kalman filter for rigid body pose estimation," *Journal of Biomechanical Engineering*, 2005.

[7] B. J. Borbély and P. Szolgy, "Real-time inverse kinematics for the upper limb: a model-based algorithm using segment orientations," *Biomedical engineering online*, vol. 16, no. 1, pp. 1-29, 2017.

[8] A. D. Vigotsky, I. Halperin, G. J. Lehman, G. S. Trajano, and T. M. Vieira, "Interpreting signal amplitudes in surface electromyography studies in sport and rehabilitation sciences," *Frontiers in physiology*, p. 985, 2018.

[9] A. J. Van den Bogert, T. Geijtenbeek, O. Even-Zohar, F. Steenbrink, and E. C. Hardin, "A real-time system for biomechanical analysis of human movement and muscle function," *Med. Biol. Eng. Comput.*, vol. 51, no. 10, pp. 1069-1077, 2013.

[10] R. Dumas, R. Brånemark, and L. Frossard, "Gait analysis of transfemoral amputees: errors in inverse dynamics are substantial and depend on prosthetic design," *IEEE Tran. on Neural Sys. and Rehab. Eng.*, vol. 25, no. 6, pp. 679-685, 2016.

[11] A. Kian, C. Pizzolato, M. Halaki, K. Ginn, D. Lloyd, D. Reed, and D. Ackland, "Static optimization underestimates antagonist muscle activity at the glenohumeral joint: A musculoskeletal modeling study," *J. Biomech.*, vol. 97, p. 109348, 2019.

[12] C. Bélaise, F. Dal Maso, B. Michaud, K. Mombaur, and M. Begon, "An emg-marker tracking optimisation method for estimating muscle forces," *Multibody Sys. Dynamics*, vol. 42, no. 2, pp. 119-143, 2018.

[13] F. Moissenet, C. Bélaise, E. Piche, B. Michaud, and M. Begon, "An optimization method tracking emg, ground reactions forces, and marker trajectories for musculo-tendon forces estimation in equinus gait," *Front in neurorobotics*, vol. 13, p. 48, 2019.

[14] C. Pizzolato, D. G. Lloyd, M. Sartori, E. Ceseracciu, T. F. Besier, B. J. Fregly, and M. Reggiani, "Ceinms: A toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks," *J. Biomech.*, vol. 48, no. 14, pp. 3929-3936, 2015.

[15] M. S. Shourijeh, N. Mehrabi, and J. McPhee, "Forward static optimization in dynamic simulation of human musculoskeletal systems: a proof-of-concept study," *Journal of Computational and Nonlinear Dynamics*, vol. 12, no. 5, 2017.

[16] N. Mehrabi, R. Sharif Razavian, B. Ghannadi, and J. McPhee, "Predictive simulation of reaching moving targets using nonlinear model predictive control," *Frontiers in computational neuroscience*, vol. 10, p. 143, 2017.

[17] F. Bailly, A. Ceglia, B. Michaud, D. M. Rouleau, and M. Begon, "Real-time and dynamically consistent estimation of muscle forces using a moving horizon emg-marker tracking algorithm—application to upper limb biomechanics," *Front. in Bioeng. and Biotech.*, vol. 9, p. 112, 2021.

[18] W. Wu, P. V. Lee, A. L. Bryant, M. Galea, and D. C. Ackland, "Subject-specific musculoskeletal modeling in the evaluation of shoulder muscle and joint function," *J. Biomech.*, vol. 49, no. 15, pp. 3626-3634, 2016.

[19] J. M. Winters, *Hill-Based Muscle Models: A Systems Engineering Perspective*. New York, NY: Springer New York, 1990, pp. 69-93.

[20] K. R. Holzbaur, W. M. Murray, and S. L. Delp, "A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control," *Annals of biomedical engineering*, vol. 33, no. 6, pp. 829-840, 2005.

[21] S. Li, P. Dario, and Z. Song, "Prediction of passive torque on human shoulder joint based on bpann," *Applied Bionics and Biomechanics*, vol. 2020, 2020.

[22] M. Millard, T. Uchida, A. Seth, and S. L. Delp, "Flexing computational muscle: modeling and simulation of musculotendon dynamics," *Journal of biomechanical engineering*, vol. 135, no. 2, 2013.

[23] R. Neptune and S. Kautz, "Muscle activation and deactivation dynamics: the governing properties in fast cyclical human movement performance?" *Exercise and sport sciences reviews*, vol. 29, no. 2, pp. 76-81, 2001.

[24] F. De Groote, A. L. Kinney, A. V. Rao, and B. J. Fregly, "Evaluation of direct collocation optimal control problem formulations for solving the muscle redundancy problem," *Annals of biomedical engineering*, vol. 44, no. 10, pp. 2922-2936, 2016.

[25] C. Bélaise, B. Michaud, F. Dal Maso, K. Mombaur, and M. Begon, "Which data should be tracked in forward-dynamic optimisation to best predict muscle forces in a pathological co-contraction case?" *J. Biomech.*, vol. 68, pp. 99-106, 2018.

[26] E. M. Jafarov, M. N. A. Parlakçi, and Y. Istefanopulos, "A new variable structure pid-controller design for robot manipulators," *IEEE Tran. on Control Sys. Tech.*, vol. 13, no. 1, pp. 122-130, 2004.

[27] A. Venne, F. Bailly, E. Charbonneau, J. Dowling-Medley, and M. Begon, "Optimal estimation of complex aerial movements using dynamic optimisation," *Sports biomechanics*, pp. 1-16, 2022.

[28] B. Michaud, F. Bailly, E. Charbonneau, A. Ceglia, L. Sanchez, and M. Begon, "Bioptim, a python framework for musculoskeletal optimal control in biomechanics," *IEEE Tran. on Sys., Man, and Cybernetics: Sys.*, 2022.

[29] R. Verschueren, G. Frison, D. Kouzoupis, J. Frey, N. v. Duijkeren, A. Zanelli, B. Novoselnik, T. Albin, R. Quirynen, and M. Diehl, "acados—a modular open-source framework for fast embedded optimal control," *Math. Prog. Computation*, vol. 14, no. 1, pp. 147-183, 2022.

[30] M. Diehl, H. G. Bock, and J. P. Schlöder, "A real-time iteration scheme for nonlinear optimization in optimal feedback control," *SIAM J Control Optim.*, vol. 43, no. 5, pp. 1714-1736, 2005.

[31] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, "OSQP: An operator splitting solver for quadratic programs," *Math. Program. Comput.*, vol. 12, no. 4, pp. 637-672, 2020.

[32] G. Frison and M. Diehl, "Hpipm: a high-performance quadratic programming framework for model predictive control," *IFAC-PapersOnLine*, vol. 53, no. 2, pp. 6563-6569, 2020.

[33] R. Verschueren, M. Zanon, R. Quirynen, and M. Diehl, "A sparsity preserving convexification procedure for indefinite quadratic programs arising in direct optimal control," *SIAM Journal on Optimization*, vol. 27, no. 3, pp. 2085-2109, 2017.

[34] R. Hildebrand, "Optimal step length for the newton method: case of self-concordant functions," *Mathematical Methods of Operations Research*, vol. 94, no. 2, pp. 253-279, 2021.

[35] S. L. Campbell and Y. Zhong, "Jacobian reuse in explicit integrators for higher index daes," *Applied Num. Math.*, vol. 25, no. 4, pp. 391-412, 1997.

[36] Y. Blache, R. Dumas, A. Lundberg, and M. Begon, "Main component of soft tissue artifact of the upper-limbs with respect to different functional, daily life and sports movements," *J. Biomech.*, vol. 62, pp. 39-46, 2017.

[37] F. Dal Maso, P. Marion, and M. Begon, "Optimal combinations of isometric normalization tests for the production of maximum voluntary activation of the shoulder muscles," *Archives of physical medicine and rehabilitation*, vol. 97, no. 9, pp. 1542-1551, 2016.

[38] S. L. Delp, F. C. Anderson, A. S. Arnold, P. Loan, A. Habib, C. T. John, E. Guendelman, and D. G. Thelen, "Opensim: open-source software to create and analyze dynamic simulations of movement," *IEEE Tran. on Biomed. Eng.*, vol. 54, no. 11, pp. 1940-1950, 2007.

[39] B. Michaud and M. Begon, "biorbd: a c++, python and matlab library to analyze and simulate the human body biomechanics," *Journal of Open Source Software*, vol. 6, no. 57, p. 2562, 2021.

[40] "Datastream SDK, Stream 3rd Party Programs and Plug-Ins, Vicon." [Online]. Available: https://www.vicon.com/software/datastream-sdk/

[41] J. Wickham, T. Pizzari, K. Stansfeld, A. Burnside, and L. Watson, "Quantifying ‘normal’shoulder muscle activity during abduction," *J Electromyogr Kinesiol*, vol. 20, no. 2, pp. 212-222, 2010.

[42] M. Mitchell, B. Muftakhidinov, T. Winchen, A. Wilms, B. van Schaik, badshah400, Mo-Gul, T. G. Badger, Z. Jedrzejewski-Szmek, kensington, and kylesower, "markummitchell/engauge-digitizer: Nonrelease (v12.2.1)," July 2020.

[43] A. Ceglia, F. Verdugo, and M. Begon, "Biosiglive: an open-source python package for real-time biosignal processing," *Journal of Open Source Software*, vol. 8, no. 83, p. 5091, 2023.

[44] I. Campanini, A. Merlo, P. Degola, R. Merletti, G. Vezzosi, and D. Farina, "Effect of electrode location on emg signal envelope in leg muscles during gait," *Journal of Electromyography and Kinesiology*, vol. 17, no. 4, pp. 515-526, 2007.

[45] C. Pizzolato, M. Reggiani, D. J. Saxby, E. Ceseracciu, L. Modenese, and D. G. Lloyd, "Biofeedback for gait retraining based on real-time estimation of tibiofemoral joint contact forces," *IEEE Tran. on Neural Systems and Rehabilitation Eng.*, vol. 25, no. 9, pp. 1612-1621, 2017.

[46] Z. Li, D. Guiraud, D. Andreu, A. Gelis, C. Fattal, and M. Hayashibe, "Real-time closed-loop functional electrical stimulation control of muscle activation with evoked electromyography feedback for spinal cord injured patients," *Intern. J. of Neural Sys.*, vol. 28, no. 06, p. 1750063, 2018.