EMG Model Calibration and Trajectory Tracking by Remote Rehabilitation Exoskeleton

1. Introduction

From time to time, mankind has had to confront global infectious diseases with deadly impact, as has recently been the case with the Coronavirus, COVID-19. Family members have been compelled to take on in-home responsibilities due to quarantine, including patients following medical treatments. Pandemic-related robotic systematized solutions have been described for multiple scenarios [1]. Similarly, telesystems, communication networks, intelligent machines, and information technologies provide support for easy post-pandemic scenarios [2,3] by keeping us away from contagious scenarios [4,5]. The current technological market has provided a spectrum of modern communication technology such as distributed networked equipment, Internet-based telesystems, telepresence devices, either synchronous or asynchronous modalities [6]. Telesystems applications have evolved into Internet shopping, remote academic/job meetings, virtual scholarly classes, and healthcare support, such as tele-rehabilitation therapies.

An exoskeleton is a wearable robot fitted on a human’s body, designed to enhance biomechanical functions, with the human serving as the controller [12,13]. To be effective, an exoskeleton must provide a sense of lightness and agility, akin to the human bone structure. The advantages of exoskeleton technology extend to enhancing biomechanical capabilities, addressing diverse needs such as the substitution of body parts’ mobility, providing locomotion assistance [7,8], and aiding in rehabilitation [10,11,15,16]. Rehabilitation is particularly crucial in cases of post-corrective surgeries, where constant therapy is required to expedite recovery. Demand for exoskeletons spans various sectors, including athletics, sports, manual labor, and the military, where common surgeries involve conditions like rotator cuff injuries, acromioclavicular articulation injuries, elbow dislocations, hip dislocations, reconstruction of injured anterior cruciate ligaments, knee arthroscopy, and total knee replacement. Prolonged post-surgery rehabilitation plays a critical role in recovery, as prolonged limb immobilization can lead to a significant reduction in muscle strength and mobility. Additionally, various studies have explored control algorithms for lower-limb wearable robot systems [17], as well as for upper limbs [24], to provide power assistance [23]. Wearable robotic hands have been developed for purposes such as haptic feedback, assistive functions, and rehabilitation [34]. These devices also find applications in tasks like grasping [32], the development of multimodal rehabilitation exoskeletons [20,30], and walking assistance [27].

The work by [28] introduced an adaptive neural control system for rehabilitating lower limbs with unknown model dynamics, focusing on adapting to the complex dynamics of human-exoskeleton interactions. In this study, we propose a feedback kinematic tracking model for a networked exoskeleton, which is taught asynchronously by a remote therapist (Figure 1abc). The therapist’s rehabilitation exercises are visually assessed and transformed into analytical polynomial Cartesian functions through a least-square fitting approach. A noteworthy innovation presented in this work is the transmission of a concise stream of polynomial coefficients, as opposed to the entire set of Cartesian data. On the patient’s side, these Cartesian polynomials are employed to recreate the therapist’s experimental data with higher resolution. Additionally, by introducing a constant time period, denoted as T, we achieve the reconstruction of Cartesian components, x and y, as periodic functions of time using Fourier series. The integration of Cartesian and electromyography (EMG) data as feedback, similar to the approach in Kang et al. [9], represents another significant contribution of this work. Notably, our model is designed to infer the posture of Cartesian limbs from EMG signals, allowing for recalibration during new exercises. This feature offers a generalized approach, enabling the control law to adapt to the structure of different exoskeletons simply by readjusting their Jacobian matrices. Furthermore, a linear trajectory smoother model is introduced in cases where patients have biomechanical movement constraints. These smoothing functions facilitate a gradual convergence towards the desired reference trajectory shape. To evaluate the efficacy of the proposed exoskeleton, we conducted extensive simulations involving hundreds of tracking control rehabilitation exercises.

This study introduces a distinctive approach to remote therapy within the realm of tele-rehabilitation. Specifically, it proposes a feedback kinematic tracking model for a networked exoskeleton that undergoes asynchronous visual instruction from a remote therapist. The emphasis lies in transmitting concise streams of polynomial coefficients instead of the entire set of Cartesian data, subsequently transformed into analytical polynomial Cartesian functions. Additionally, the study underscores the integration of Cartesian and electromyography (EMG) data as feedback, facilitating the inference of limb posture from EMG signals and recalibration during new exercises. In the existing literature various approaches within closely related fields of study, including works by [36,37,38,39,40,41], partially intersect with the themes presented in this work. However, the authors remain uncertain as to whether fully analogous studies have been previously published.

The proposed solution specifically focuses on tele-rehabilitation, where a therapist’s exercises are visually assessed and transmitted to the patient’s exoskeleton. The key innovation of the proposed solution lies in the transmission of concise streams of polynomial coefficients, rather than transmitting the entire set of Cartesian data. This approach provides flexibility in tracking and control, making the rehabilitation process highly accurate and efficient. Overall, the proposed solution offers a generalized approach to tele-rehabilitation, making it more accessible and accommodating to a wide range of patient needs. It combines clinical precision with technological sophistication, providing an effective and efficient means of rehabilitation.

The proposed feedback kinematic tracking model for a networked exoskeleton taught asynchronously by a remote therapist is required to ensure accurate, personalized, and efficient rehabilitation for patients, regardless of their location or access to specialized facilities. In addition, it provides i) accuracy and precision: The feedback kinematic tracking model follows the desired trajectory and provides the necessary level of accuracy which is crucial for effective rehabilitation. ii) Adaptability and customization: The feedback model allows for adjustments and adaptations based on the patient’s specific biomechanical needs, and the use of electromyography (EMG) data as feedback allows for the inference of Cartesian posture from EMG signals, enabling recalibration during new exercises. iii) Remote guidance and supervision: Asynchronous teaching by a remote therapist regardless of their location, enables benefits for patients who may not have access to specialized rehabilitation facilities or therapists in their vicinity. iv) Flexibility and convenience: Asynchronous teaching allows patients to perform and follow the prescribed exercises and receive feedback from the therapist at a time that suits them best. v) Efficiency and cost-effectiveness: By utilizing a networked exoskeleton and asynchronous teaching, the rehabilitation process can save time and resources by avoiding frequent visits to rehabilitation centers, while still receiving expert guidance and support.

This paper is structured as follows: Section 2 delves into the therapist-side tele-rehabilitation, elucidating the teaching process facilitated by a vision system. In Section 3, we unveil the proposed EMG calibration and observation measurement models. Following this, Section 4 delves into the exoskeleton’s motion constraints and its kinematic control law. A comprehensive summary of the research’s results is provided in Section 5. Finally, Section 6 offers a set of conclusions drawn from this study.

2. Human therapist remote visual teaching

The proposed networked distributed system is structured as an asynchronous tele-system architecture, with the therapist’s side and the patients constituting the client side. The server plays a pivotal role by providing geometric data acquired visually from the therapist’s exercise sessions in front of a camera, transmitted over the network. An artificial vision system captures the data related to the training trajectories taught by the therapist (Figure 1d, left side).

Notably, the interaction session between the therapist and the patient occurs asynchronously, eliminating the need for simultaneous network connections. This approach stands in contrast to traditional tele-operation methods [19]. The exoskeleton has the flexibility to connect at any given date and time, allowing the download of therapy exercises in the form of analytic trajectories, represented by coefficients (Figure 1d, right side). The server’s computer is equipped with a monocular vision system that accurately measures the angles and positions of colored landmarks strategically placed at the therapist’s limb joints (Figure 3ab-left). Remarkably, the therapist does not require any specialized equipment during the exercise sessions.

In light of the evident benefits of rehabilitation exoskeletons ([22]), we selected commonly practiced exercises from a recognized exercise program ([18]) to serve as the basis for our study. These exercises, whether post-injury or post-surgery, were adopted as trajectories (Figure 1abc). For example, exercises such as the rotator cuff, shoulder and elbow flexion exercises were employed to teach various motions (Figure 1ac). Similarly, we incorporated exercises like knee extensions (Figure 1b) and the stand-up paddle exercise with an elastic band (Figure 1c) to broaden the scope of our research.

Figure 1. Tele-training rehabilitation exercises. a)-b) Elbow, knee flexion-extension. c) Shoulder sagittal rotation. d) Telesystem architecture.

Figure 1. Tele-training rehabilitation exercises. a)-b) Elbow, knee flexion-extension. c) Shoulder sagittal rotation. d) Telesystem architecture.

The proposed tele-system was crafted on standard-capability computers running the Linux Gentoo operating system, version 4.4. The entire codebase was developed using the GNU C/C++ compiler, making extensive use of multiple open-source libraries, including openCV, Armadillo C++, the Open Dynamics Engine (ODE), and Gazebo. Communication between components was facilitated using the TCP/IP protocol via sockets. The vision sensor employed was a cost-effective USB camera (see Table 1). For the acquisition of EMG biosignals, we utilized the gMOBIlab+ device with two EMG channels. To bring our vision to life, we implemented a virtual exoskeleton prototype with original design mechanisms. This virtual exoskeleton operates under the control of human EMG signals, functioning as a master-slave system.

This work places some emphasis on vision-based teaching. To streamline the measurement process, we strategically position colored landmarks at the joints of the therapist’s limbs, as depicted in Figure 3. These landmarks serve as the key reference points for measuring angles and positions during exercise sessions. Our approach to measuring the metric trajectories of exercises involves a two-step process. First, we segment the joints’ landmarks based on their RGB color. Subsequently, we employ a nominal resolution geometric model to track the centroid positions of these metric landmarks. In this approach, we establish the reference coordinate system at the location of the visual sensor. In line with the extrinsic parameters outlined in Figure 2, we define an RGB image as a cubic matrix $\mathbf{I}\in {\mathbb{R}}^{m\times n\times 3}$, comprising the union of the three primary color channels: red ${\mathbf{I}}_R$, green ${\mathbf{I}}_G$, and blue ${\mathbf{I}}_B$. To extract the shoulder’s landmark, we apply a threshold ${\zeta }_S$ to the red channel ${\mathbf{I}}_R$ and incorporate a zeros matrix $\mathbf{0}\in {\mathbb{R}}^{m\times n}$. This process yields a binary segmented image denoted as ${\mathbf{I}}_S$ and is achieved through the following Boolean equation:

$${\mathbf{I}}_S=({\mathbf{I}}_R>{\zeta }_S)\cup ({\mathbf{I}}_G\cap \mathbf{0})\cup ({\mathbf{I}}_B\cap \mathbf{0}).$$

(1)

Similarly, the isolation of the elbow’s joint landmark is achieved through the utilization of a binary image denoted as ${\mathbf{I}}_E$. The effectiveness of this process hinges on the presence of a critical threshold denoted as ${\zeta }_E$, which assumes a pivotal role in the Boolean criterion:

$${\mathbf{I}}_E=({\mathbf{I}}_R\cap \mathbf{0})\cup ({\mathbf{I}}_G<{\zeta }_E)\cup ({\mathbf{I}}_B\cap \mathbf{0}).$$

(2)

In a similar manner, the identification of the wrist landmark is carried out through the utilization of a binary image represented as ${\mathbf{I}}_W$ The determination of this landmark relies on the precise specification of a threshold value denoted as ${\zeta }_W$ which is an integral part of the following criterion:

$${\mathbf{I}}_W=(({\mathbf{I}}_R\cap \mathbf{0})\cup ({\mathbf{I}}_G\cap \mathbf{0}))\cup ({\mathbf{I}}_B>{\zeta }_W).$$

(3)

For each segmented region in the images ${\mathbf{I}}_R$, ${\mathbf{I}}_G$, and ${\mathbf{I}}_B$, the position of the landmarks in spatial coordinates is determined by their centroids, denoted as column $p_c$ and row $p_r$, as expressed by:

$$p_c=\frac{1}{2}\sum _i{c}_i,p_r=\frac{1}{2}\sum _i{r}_i,$$

(4)

In this context, the spatial coordinates $c_i$ and $r_i$, where both $c_i$ and $r_i$ belong to a segmented region $\mathcal{R}$, represent the pixels within the region. We will now proceed to offer a comprehensive description of the geometric visual model, as depicted in Figure 2.

Figure 2. Monocular vision-based geometrical model on the therapist’s side.

Figure 2. Monocular vision-based geometrical model on the therapist’s side.

Geometrically, the equality for the vertical field of view ((5)) is valid.

$$tan\left(\frac{{\phi }_r}{2}\right)=\left(\frac{h/2}{z}\right),$$

(5)

Here, h represents the height of the object of interest, and z denotes the metric depth of the human therapist from the camera. Thus, building upon the previous expression, we can determine ${\phi }_r$, which represents the vertical camera angle quantized in terms of image rows.

$${\phi }_r=2{tan}^{-1}\left(\frac{h}{2z}\right).$$

(6)

With respect to (w.r.t.) the horizontal field of view,

$$tan\left(\frac{{\phi }_c}{2}\right)=\left(\frac{w/2}{z}\right),$$

(7)

then, dropping off the angle along image’s columns ${\phi }_c$

$${\phi }_c=2{tan}^{-1}\left(\frac{w}{2z}\right).$$

(8)

We estimate the vertical metric resolution factor, denoted as $f_h$, as well as the horizontal metric resolution factor, denoted as $f_w$1,

$$f_h=2tan\left(\frac{{\phi }_r}{2}·z\right),f_w=2tan\left(\frac{{\phi }_c}{2}·z\right).$$

(9)

Subsequently, we determine the vertical metric position ${\Delta }_r$ and the horizontal metric position ${\Delta }_c$ of the center pixel $(p_r,p_c)$ of a landmark.

$${\Delta }_r=\frac{f_h}{{\mathrm{p}}_{\mathrm{r}}},{\Delta }_c=\frac{f_w}{{\mathrm{p}}_{\mathrm{c}}}.$$

(10)

Next, in the therapist’s side, the superior limb kinematic motion model in 3D space is provided by the expression

$$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{c}{\ell }_{1}sin\left({\varphi }_0\right)cos\left({\varphi }_1\right)+{\ell }_{2}sin\left({\varphi }_0\right)cos({\varphi }_1+\varphi 2)\\ {\ell }_{1}sin\left({\varphi }_1\right)+{\ell }_{2}sin({\varphi }_1+\varphi 2)\\ {\ell }_{1}cos\left({\varphi }_0\right)cos\left({\varphi }_1\right)+{\ell }_{2}cos\left({\varphi }_0\right)cos({\varphi }_1+\varphi 2)\end{array}\right),$$

where depth is affected by the yaw angle ${\varphi }_0$ and double pitch joits ${\varphi }_1,{\varphi }_2$, and ${\ell }_1$ and ${\ell }_2$ to denote the actual metric lengths of the upper arm and forearm, respectively. When representing a 2D exercise scenario, the yaw angle ${\varphi }_0=\pi /2$. Therefore, the depth dimension becomes irrelevant since the arm’s length is not projected, and the trajectory is drawn on a plane in front of the camera. The links ${\ell }_1,{\ell }_2$ visual measurements,

$${\ell }_1=\sqrt[2]{{(x_S-x_E)}^2+{(y_S-y_E)}^2}$$

(11)

and

$${\ell }_2=\sqrt[2]{{(x_E-x_W)}^2+{(y_E-y_W)}^2}.$$

(12)

Similarly, we use the variables ${\theta }_1$ and ${\theta }_2$ to denote the joint angles of the limb.

$${\theta }_1={tan}^{-1}\left(\frac{y_E-y_S}{x_E-x_S}\right),{\theta }_2={tan}^{-1}\left(\frac{y_W-y_E}{x_W-x_E}\right)$$

(13)

We illustrate a shoulder exercise in Figure 3, and in Figure 9, we display the resulting Cartesian trajectory, with red points indicating the Cartesian positions of the blue wrist landmark.

Figure 3. Visual instruction of various exercises by the therapist.

Figure 3. Visual instruction of various exercises by the therapist.

In summary, the rehabilitation trajectory, denoted as $\mathcal{S}={\Theta }_0,{\Theta }_1,{\Theta }_2,\cdots ,{\Theta }_N$, is comprised of a set of joint angles ${\Theta }_i={({\theta }_1,{\theta }_2)}^{\top }$ and limb parameters ${\ell }_1$ and ${\ell }_2$.

In Figure 4, a series of repetitions stemming from exercises performed by the therapist is presented. These exercises encompass a single experimental repetition involving shoulder flexion (depicted in Figure 4a), along with three experimental repetitions focusing on elbow flexion (illustrated in Figure 4b). Additionally, three experimental repetitions of a rowing-like motion are depicted (Figure 4c), wherein both the shoulder and elbow joints operate in unison during the exercise. It is imperative to emphasize that the reference coordinates, denoted as $(0,0)$, are firmly anchored at the centroid of the shoulder.

Additionally, in Figure 5, we present angular motion measurements for all joints utilizing visual odometry, with a focus on the $x-y$ plane. These angular measurements correspond to the Cartesian motion data depicted in Figure 4. In Figure 5a (left), we observe an experiment exclusively involving shoulder flexion, where a linear rigid arm is used, and there are a total of five repetitions. It’s important to note that this exercise is the simplest, resulting in data points with a higher degree of overlap. On the right side of Figure 5a, while still involving shoulder flexion, we observe slight angular perturbations in the wrist, averaging approximately $9^{\circ }$. These perturbations were primarily caused by minor limb vibrations from the therapist and were further exacerbated by contrast and brightness noise in the measurements. Figure 5b displays angular measurements for elbow flexion exercises. To ensure clarity without excessive data overlap, we have chosen to present three distinct experiments, each with multiple repetitions. On the left side of Figure 5b, we observe angular perturbations in the wrist during elbow flexion exercises. Conversely, the right side of Figure 5b depicts three separate experiments involving elbow-flexion exercises, displaying multiple repetitions that result in periodic oscillations. This presentation highlights variations in speed during the experiments, reflecting biomechanical consistency and differences in execution.

Figure 5c, in its entirety, represents three distinct experiments of rowing-like motions, each of which consists of six repetitions. This comprehensive view enables us to examine the angular dynamics in these exercises in greater detail. On the left side of Figure 5c, we can observe the angular magnitudes of elbow flexions. Meanwhile, the right side of the same figure portrays shoulder flexions, which synchronize their motion with the elbow flexions. On the therapist’s side, the process involves employing a least square polynomial fitting approach to derive a Cartesian function $y\left(x\right)$. This method is used to process the vision-based 3D measured data, resulting in the determination of polynomial coefficients of the n-th degree, as described in the following equation:

$$\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right)={\left(\begin{array}{cccc}n& \sum _i{x}_i& \cdots & \sum _i{x}_i^n\\ \sum _i{x}_i& \sum _i{x}_i^2& \cdots & \sum _i{x}_i^{n+1}\\ ⋮& & & \\ \sum _i{x}_i^j& \sum _i{x}_i^{j+1}& \cdots & \sum _i{x}_i^{2n}\end{array}\right)}^{-1}·\left(\begin{array}{c}\sum y_i\\ \sum x_i{y}_i\\ ⋮\\ \sum _i{x}_i^n{y}_i\end{array}\right),$$

(14)

In this manner, we obtain a parametric function that accurately fits the exercise’s Cartesian data, following a polynomial form:

$$y\left(x\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+\cdots +a_n{x}^n.$$

(15)

This approach significantly reduces the trade-off during data streaming to only a few coefficients and the time period T of one exercise oscillation. These coefficients require a minimal amount of bytes for transmission over the network, thereby mitigating significant time delays. Figure 6 illustrates experimental data from teaching trajectories that have been effectively fitted by various polynomials, exemplifying the practicality and efficiency of this approach.

By representing the exercise trajectory as a continuous polynomial function, denoted as $y\left(x\right)$, available on the exoskeleton’s side, we open the door to the replication of experimental data at higher metrical resolutions. With the time constant of a period T at our disposal, the exoskeleton’s side can autonomously convert the high-resolution Cartesian data into new components as functions of time. Using Fourier series, we generate new Cartesian components, $x\left(t\right)$ and $y\left(t\right)$, as functions of time. Figure 7 effectively demonstrates this process, showcasing a Cartesian reference trajectory as $x\left(t\right)$ versus $y\left(t\right)$, realistically reproducing three experiments involving elbow flexion. This approach, as presented in this work, offers sophisticated analytical functions that serve as effective substitutes for the therapist’s exercises, which typically involve large amounts of experimental data. In doing so, it efficiently addresses the challenges of long transmission time delays and the handling of extensive data tables, ultimately yielding parametric functions that streamline the entire process.

Moreover, in scenarios where biomechanical constraints affecting a patient’s movements are detected during exoskeleton deployment, we employ smoothing functions to facilitate gradual transitions that accommodate these constraints. One such smoothing function, based on a linear approach, is detailed below:

$$y_{ref}={\kappa }_{\tau }+{\kappa }_{s}y\left(x\right),$$

(16)

In this equation, the original therapist’s trajectory is regarded as the ideal Cartesian reference motion, denoted as $y\left(x\right)$. Therefore, we introduce two key parameters: ${\kappa }_{\tau }$, representing the translation factor of a family of functions, and ${\kappa }_s$, representing the scaling function’s coefficient. When ${\kappa }_{\tau }=0$ and ${\kappa }_s=1$, the exercise function aligns precisely with the therapist’s reference trajectory. However, by adjusting the values of these smoothing coefficients, we achieve gradual convergence toward the therapist’s ideal trajectory. Figure 8 purposely illustrates various smoothing teaching trajectories that progressively approach the original therapist’s reference trajectory, represented by the continuous line, underscoring the flexibility and efficacy of this approach.

Indeed, it’s important to note that both ${\kappa }_{\tau }$ and ${\kappa }_s$ are coefficients specific to each exercise trajectory, differences of patient’s limbs size or patient’s motion conditions. The coefficients are manually adjusted individually for different exercises based on the patient’s comfort. This adaptability ensures that the smoothing coefficients are tailored to the characteristics of each exercise, optimizing their effectiveness.

3. Rehabilitation trajectory model

From the perspective of the patient-exoskeleton interface, the instructional trajectory $\mathcal{S}$ is inherently expressed in polar form. To facilitate processing, it is subsequently transformed into Cartesian points at the wrist location, represented as $x_i^{ref}={\ell }_{1}cos\left({\theta }_1\right)+{\ell }_{2}cos({\theta }_1+{\theta }_2)$ and $y_i^{ref}={\ell }_{1}sin\left({\theta }_1\right)+{\ell }_{2}sin({\theta }_1+{\theta }_2)$, commonly referred to as the ideal or reference trajectory ${\mathbf{p}}^{ref}$. In lieu of managing an extensive dataset comprising hundreds or thousands of points, this study derives an adaptive analytical reference model expressed as functions of time, akin to the approach proposed by [25]. Given the periodic nature of a rehabilitation exercise assumed to be a Cartesian function, the suggested methodology involves obtaining xx and yy separately with respect to time. The utilization of a general Fourier series is deemed apt for representing a suitable function of time, providing a concise and analytically manageable characterization of the patient’s movement trajectory.

$$x\left(t\right)={\alpha }_0^x+{\alpha }_k^{x}cos\left(k\omega t\right)+{\beta }_k^{x}sin\left(k\omega t\right)+\cdots$$

(17a)

and

$$y\left(t\right)={\alpha }_0+{\alpha }_k^{y}cos\left(k\omega t\right)+{\beta }_k^{y}sin\left(k\omega t\right)+\cdots$$

(17b)

In the general formulation (18)-(20) the Fourier series estimate the coefficients for (17). The parameter n represents the number of data to be fitted through the Fourier series. However, the $k_{th}$ number of Fourier coefficients are predetermined with a sufficiently large value to appropriately accommodate the trajectory data. Equations 18, 19, and 20 comprise the Fourier series tailored to fit the trajectory data provided by the therapist, expressing it as a trigonometric function.

A C++ program was implemented to autonomously compute the discrete Fourier coefficients. Regarding the independent terms, the average of the discrete data is presented as follows:

$${\displaystyle {\alpha }_0^x=\frac{1}{n}\sum _{k=1}^n{x}_k,\mathrm{and}{\alpha }_0^y=\frac{1}{n}\sum _{k=1}^n{y}_k.}$$

(18)

Similarly, even though initially presented in continuous-time, the discrete Fourier series was formulated to derive the remaining coefficients for the Cartesian component x vs t by

$${\alpha }_k^x=\frac{2}{n}{\int }_{t}x\left(t\right)cos\left(k\omega t\right)\mathrm{d}t,{\beta }_k^x=\frac{2}{n}{\int }_{t}x\left(t\right)sin\left(k\omega t\right)\mathrm{d}t,$$

(19)

likewise, for the Cartesian component y vs t,

$${\alpha }_k^y=\frac{2}{n}{\int }_{t}y\left(t\right)cos\left(k\omega t\right)\mathrm{d}t,{\beta }_k^y=\frac{2}{n}{\int }_{t}y\left(t\right)sin\left(k\omega t\right)\mathrm{d}t.$$

(20)

Therefore, we have now acquired the Cartesian components of an $xy$ trajectory as reference functions of time, denoted as $x\left(t\right)$ and $y\left(t\right)$. The numerical coefficients related to time obtained to align with the experimental reference data of shoulder flexion-extension, illustrated in Figure 9, are:

$$\begin{array}{c}\hfill x_t^{ref}=0.17-0.05cos\left(\omega t\right)-0.06cos\left(2\omega t\right)+0.01\left(sin\left(2\omega t\right)-sin\left(3\omega t\right)\right)\end{array}$$

(21)

and

$$\begin{array}{c}\hfill y_t^{ref}=-0.30-0.22cos\left(\omega t\right)+0.05sin\left(\omega t\right)-0.02sin\left(2\omega t\right)+0.01\left(cos\left(3\omega t\right)-sin\left(3\omega t\right)\right).\end{array}$$

(22)

The experimental coefficients obtained are presented for illustrative purposes, specifically related to Figure 9. It’s essential to note that their numeric values may vary for different exercises or when administered by a human therapist, accounting for differences in limb size.

Figure 9. Periodic positions of shoulder flexion-extension fitted by Fourier series.

Figure 9. Periodic positions of shoulder flexion-extension fitted by Fourier series.

3.1. EMG sensing model calibration

This study suggests deriving a sensing model for estimating the Cartesian position of the limb through EMG measurements. Given the challenge of achieving repeatability in EMG data from the same exercises and subjects, an EMG calibration model is established at the onset of training. Real-time motion prediction becomes feasible using raw EMG data, typically encompassing combined muscle motions [33]. Table 2 provides a comprehensive overview, detailing the type of therapy and the associated muscles involved in obtaining EMG measurements.

Based on the electrodes illustrated in Figure 10c, the experimental raw EMG data is presented in Figure 14. This data pertains to the shoulder-elbow-wrist region and corresponds to multiple periodic flexion-extension exercises.

The series of plots illustrated in Figure 11 portrays three distinct experiments involving the same exercise, specifically elbow flexions, each comprising approximately five repetitions. In this configuration, Channel Ch5 was situated at the bicep brachii muscle, while concurrently, Channel Ch6 was positioned at the tricep brachial muscle. Notably, the magnitudes of amplitude ($\mu V$) exhibit coherence and consistency across both channels.

The subsequent Figure 12 illustrates an experiment involving shoulder flexion. In this instance, Channel Ch5 was positioned at the Pectoralis Major muscle, and Channel Ch6 was situated at the Deltoid Anterior muscle. The plots represent approximately five oscillations (flexions) with amplitude magnitudes measured in $\mu V$. It is noteworthy that, in this case, both frequency and amplitudes consistently differ from those observed in other exercises.

Figure 13 presents experimental data related to rowing-like movements, showcasing three distinct experiments involving the same exercise with approximately six rowing flexions each. Channel Ch5 was positioned at the Trapezius middle muscle, while Channel Ch6 was situated at the Deltoid anterior muscle. The amplitude magnitudes ($\mu V$) and oscillation frequency exhibit distinct behavioral patterns, maintaining substantial consistency across both channels.

Figure 13. Three different experiments of rowing exercise using EMG channels Ch5 (Trapezius middle) and Ch6 (Deltoid anterior).

Figure 13. Three different experiments of rowing exercise using EMG channels Ch5 (Trapezius middle) and Ch6 (Deltoid anterior).

Figure 14. Native EMG signals $\tilde{E}$ during simultaneous shoulder-elbow flexion extension and wrist.

Figure 14. Native EMG signals $\tilde{E}$ during simultaneous shoulder-elbow flexion extension and wrist.

The activity level of muscles is inherently noisy due to the interaction of bi-joint muscles with movements in other joints. The EMG calibration measurement is acquired by initially filtering the noisy native EMG signal $\tilde{E}$ (Figure 14) into a linear envelope ${\overline{E}}_t$ within segments of $N_k$ samples. Subsequently, a 2^nd order low-pass filter with a cutoff frequency $f_c=20Hz$ is applied,

$${\overline{E}}_t=\left|\tilde{E}-\frac{1}{N_k}\sum _1^{N_k}{\tilde{E}}_i\right|,E_t^{cal}=\frac{{\overline{E}}_t}{\sqrt[2]{1+{\left(\frac{f_t}{f_c}\right)}^4}}.$$

(23)

Consequently, an instantaneous EMG calibration measurement value, denoted as $E_t^{cal}$, is deduced from the periodicity of the filtered EMG signal ${\overline{E}}_t$. To obtain an analytical representation of $E^{cal}$ over a period of a specified exercise, a fitness model employing a Fourier series in terms of time, denoted as $E_t^{cal}$, is derived.

$$\begin{array}{c}{E}_t^{cal}=56.60-29.65cos\left(\omega t\right)-15.53sin\left(\omega t\right)\hfill \\ \hfill -2.18cos\left(2\omega t\right)+7.73sin\left(2\omega t\right)-1.37cos\left(3\omega t\right)\\ \hfill +3.26sin\left(3\omega t\right)+1.76cos\left(4\omega t\right)-2.15sin\left(4\omega t\right)\\ \hfill +0.13cos\left(5\omega t\right)+2.58sin\left(5\omega t\right).\end{array}$$

(24)

The numerical coefficients in (24) are derived from the experimental shoulder-elbow flexion-extension exercise shown in Figure 14. Let’s introduce a polynomial function $\widehat{E}\left(t\right)$ as a time-dependent function capable of estimating EMG measurements over time. This function can then be compared against the calibrated model $E_t^{cal}$, which serves as a reference function over time. To approximate $E_t^{cal}$, we employ the Taylor series (Maclaurin), representing it as $\widehat{E}\left(t\right)$.

$$\widehat{E}\left(t\right)\approx \frac{E_0^{\left(0\right)}}{0!}{t}^0+\frac{E_0^{\left(1\right)}}{1!}t+\frac{E_0^{\left(2\right)}}{2}{t}^2+\cdots \frac{E_0^{\left(m\right)}}{m!}{t}^m.$$

(25)

A polynomial approximation is advantageous for an observer due to its simplicity in terms of differentiation and integration with respect to time. In Figure 15a, a gradual polynomial fitting up to the 80^th degree aligns with our calibration model (24). Additionally, Figure 15b visually represents the correlation between raw EMG data, overlaid with calibrated signals (filtered out), and their corresponding Cartesian positions over time during periodic exercises.

3.2. Observability based on real-time EMG

Given that the polynomial observation model $\widehat{E}\left(t\right)$ is expressed as a function of time t, its inverse model $t\left({\widehat{E}}_t\right)$ provides a time value consistent with representing instantaneous Cartesian positions $x\left(t\right)$ and $y\left(t\right)$. In the case of the experimental EMG data from Figure 16, an 80^th grade polynomial was obtained. As obtaining a non-linear analytical solution is not trivial, its inverse result for time is numerically calculated through the Newton-Raphson method for nonlinear equations. The instantaneous measurement ${\widehat{E}}_t$ substitutes the left side of the following expression,

$$\widehat{E}\left(t\right)=a_0+a_1\widehat{E}+a_2{\widehat{E}}^2+a_3{\widehat{E}}^3+\cdots +a_i{\widehat{E}}^i$$

(26)

where the time $t_k$ is the unknown value. Therefore, the Newton-Raphson method necessitates the first-order derivative with respect to time $d\widehat{E}/dt$,

$$\widehat{\dot{E}}\left(t\right)={\displaystyle \frac{\mathrm{d}\widehat{E}\left(t\right)}{\mathrm{d}t}}=a_1+a_2{\widehat{E}}_t+2a_3{\widehat{E}}_t^2+\cdots +a_i(i-1){\widehat{E}}_t^{i-1}$$

(27)

hence, a solution for $t_k$ is provided in the following manner,

$${\widehat{t}}_k={\widehat{t}}_{k-1}-{\displaystyle \frac{\widehat{E}\left(t\right)}{\widehat{\dot{E}}\left(t\right)}},$$

(28)

recursively, the best solution approximation is iteratively refined until the relative numeric error falls below ${ϵ}_t$.

$$\left|\frac{{\widehat{t}}_k-{\widehat{t}}_{k-1}}{{\widehat{t}}_k}\right|<{ϵ}_t$$

(29)

and now Cartesian models are satisfied as a functions of time $t_k$ as well

$$\widehat{x}\left(\widehat{t}\right)=b_0+b_1\widehat{t}+b_2{\widehat{t}}^2+\cdots +b_i{\widehat{t}}^i,$$

(30)

$$\widehat{y}\left(\widehat{t}\right)=c_0+c_1\widehat{t}+c_2{\widehat{t}}^2+\cdots +c_i{\widehat{t}}^i.$$

(31)

The solution for $t_k$ establishes synchronicity between ${\widehat{E}}_t$ and the expected Cartesian positions ($x_t^{ref}$, $y_t^{ref}$) in (21) and (22), respectively. Figure 16a-b display the native data and filtered signals for two EMG channels. Similarly, Figure 16c presents only the filtered signals.

Figure 16. Shoulder flexion-extension EMG measurement.

Figure 16. Shoulder flexion-extension EMG measurement.

4. Exoskeletons morphology and control law

In this study, a proposed exoskeleton model incorporates a well-designed mechanism to facilitate standard rehabilitation exercises (Figure 18). Alternative methodologies, such as adaptive geometrical design, have been documented in previous literature [35]. The design of speed reducers (e.g., the shoulder mechanism shown in Figure 17a) primarily incorporates pulse encoders for joint rotation sensing. These mechanisms are intended to enforce limb tracking along the reference trajectories.

The design assumes all encoders have a constant angular resolution R [pulses/rev]. The number of pulses ${\eta }_t$ measures the motor’s instantaneous angle ${\theta }_0$ [rad], and its angular speed $w_0$ [rad/s].

$${\theta }_0=\frac{2\pi }{R}{\eta }_t,w_0=\frac{2\pi }{R}\left(\frac{{\eta }_{t_2}-{\eta }_{t_1}}{t_2-t_1}\right).$$

(32)

All gears’ rotary speeds are described by (33) and illustrated in Figure 17b. Gears with a radius $r_i$ yield speeds $w_i$ as functions of the motor’s velocity $w_0$,

$$\left(\begin{array}{c}{w}_1\\ w_{23}\\ w_5\\ w_7\\ w_8\end{array}\right)=w_0\left(\begin{array}{c}\frac{r_0}{r_1}\\ \frac{r_0}{r_1}\\ \frac{r_0{r}_{23}}{r_1{r}_5}\\ \frac{r_0{r}_{23}}{r_1{r}_5}\\ \frac{r_0{r}_{23}{r}_7}{r_1{r}_5{r}_8}\end{array}\right).$$

(33)

By expressing (33) in vector notation, the time-invariant linear system (34) emerges,

$$\mathbf{w}=w_0\mathbf{k}.$$

(34)

This model encapsulates the functional representation of the joints, providing a framework to deduce Cartesian positions over time for feedback purposes. In Figure 18, the kinematic parameters of the upper and lower limbs are illustrated, featuring a structured organization of joint axes.

Figure 18. Exoskeleton mechanism structure.

Figure 18. Exoskeleton mechanism structure.

Therefore, it follows a dynamic controller at the level of the exoskeleton’s joint gearing, rather than at the level of the limbs’ global kinematic-EMG calibration. Employing an Euler-Lagrange approach, the general Lagrangian model ${\mathcal{L}}_i$ for the $i^{th}$ gear, transmitting rotary kinetic and potential energy, is:

$${\mathcal{L}}_i=r_i^2{m}_i{\dot{\varphi }}_i-m_i\mathrm{g}{r}_{i}sin\left({\varphi }_i\right),$$

(35)

and developing the following differential equation

$${\tau }_i=\frac{d}{\mathrm{d}t}\left({\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}{\dot{\varphi }}_i}}\right)-{\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}{\varphi }_i}}.$$

(36)

Subsequently, by considering and substituting joint kinematics from (34), and simplifying in terms of controlled motor motion ${\varphi }_0$, the following dynamic vector $\mathit{\tau }={({\tau }_1,{\tau }_{23},{\tau }_5,{\tau }_7,{\tau }_8)}^{\top }$ is deduced,

$$\mathit{\tau }=\left(\begin{array}{c}{r}_0{r}_1{m}_1{\ddot{\varphi }}_0+m_1\mathrm{g}{r}_{1}cos\left(\frac{r_0}{r_1}{\varphi }_0\right)\\ \frac{r_{23}^2{r}_0}{r_1}{m}_{23}{\ddot{\varphi }}_0+m_{23}\mathrm{g}{r}_{23}cos\left(\frac{r_0}{r_1}{\varphi }_0\right)\\ \frac{r_0{r}_{23}{r}_5}{r_1}{m}_5{\ddot{\varphi }}_0+m_5\mathrm{g}{r}_{5}cos\left(\frac{r_0{r}_{23}}{r_1{r}_5}{\varphi }_0\right)\\ \frac{r_0{r}_{23}{r}_7^2}{r_1{r}_5}{m}_7{\ddot{\varphi }}_0+m_7\mathrm{g}{r}_{7}cos\left(\frac{r_0{r}_{23}}{r_1{r}_5}{\varphi }_0\right)\\ \frac{r_0{r}_{23}{r}_7{r}_8}{r_1{r}_5}{m}_8{\ddot{\varphi }}_0+m_8\mathrm{g}{r}_{8}cos\left(\frac{r_0{r}_{23}{r}_7}{r_1{r}_5}{\varphi }_0\right)\end{array}\right).$$

(37)

Expressing $\mathit{\tau }$ in terms of the explicit vectors for inertia $\mathbf{m}$ and gravity force $\mathit{g}$,

$$\mathit{\tau }=\left(\begin{array}{c}{r}_0{r}_1{m}_1{\ddot{\varphi }}_0\\ \frac{r_{23}^2{r}_0}{r_1}{m}_{23}\\ \frac{r_0{r}_{23}{r}_5}{r_1}{m}_5\\ \frac{r_0{r}_{23}{r}_7^2}{r_1{r}_5}{m}_7\\ \frac{r_0{r}_{23}{r}_7{r}_8}{r_1{r}_5}{m}_8\end{array}\right){\ddot{\varphi }}_0+\left(\begin{array}{c}{m}_1\mathrm{g}{r}_{1}cos\left(\frac{r_0}{r_1}{\varphi }_0\right)\\ m_{23}\mathrm{g}{r}_{23}cos\left(\frac{r_0}{r_1}{\varphi }_0\right)\\ m_5\mathrm{g}{r}_{5}cos\left(\frac{r_0{r}_{23}}{r_1{r}_5}{\varphi }_0\right)\\ m_7\mathrm{g}{r}_{7}cos\left(\frac{r_0{r}_{23}}{r_1{r}_5}{\varphi }_0\right)\\ m_8\mathrm{g}{r}_{8}cos\left(\frac{r_0{r}_{23}{r}_7}{r_1{r}_5}{\varphi }_0\right)\end{array}\right)$$

(38)

and expressing this relationship in vector notation, the following dynamic equation is formulated,

$$\mathit{\tau }=\mathbf{m}{\ddot{\varphi }}_0+\mathit{g}\left({\varphi }_0\right).$$

(39)

It follows that a controller is deduced from the dynamic equation (39). Thus, rewriting the torque vector with a differential equation notation,

$$\mathit{\tau }=\mathbf{m}\frac{\mathrm{d}\dot{\varphi }}{\mathbf{d}t}+\mathit{g}\left({\varphi }_0\right)$$

(40)

and by separating different differentials,

$$\left(\mathit{\tau }-\mathit{g}\left({\varphi }_0\right)\right)\mathrm{d}t=\mathbf{m}\mathrm{d}{\dot{\varphi }}_0$$

(41)

solving by definite integration,

$$\left(\mathit{\tau }-\mathit{g}\left({\varphi }_0\right)\right){\int }_{t_2}^{t_1}\mathrm{d}t=\mathbf{m}{\int }_{{\varphi }_1}^{{\varphi }_2}\mathrm{d}{\dot{\varphi }}_0,$$

(42)

then by algebraically arranging and changing subindex notation, a model-based controller is yielded

$${\dot{\varphi }}_{0_{t+1}}={\dot{\varphi }}_{0_t}+{\mathbf{m}}^{+}\left({\mathit{\tau }}^{ref}-\mathit{g}\left({\varphi }_0\right)\right)\left(t_2-t_1\right),$$

(43)

where ${\mathbf{m}}^{+}$ is the inertia vector pseudoinverse2. Thus, the speeds of the exoskeleton’s joint mechanisms are recursively controlled according to a specified reference torque ${\mathit{\tau }}^{ref}$.

4.1. Limbs kinematic model

Based on the mechanical model presented in Figure 18, a set of Cartesian equations was deduced for kinematic analysis, forming the foundation for position control. The non-linear and non-square systems obtained in (44) and (45) evidently lack inverse analytical solutions. Therefore, an algebraic numeric recursive approach was developed. Expressions (44a)-(44c) describes the upper limb kinematics of the 5 rotary joints. Notice that ${\varphi }_0$ and ${\varphi }_3$ are rolling axis. Thus, for x,

$$x=l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_3{\mathrm{s}}_{123},$$

(44a)

for coordiante y

$$y=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(44b)

and for coordinate z,

$$z=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0.$$

(44c)

Expressions (45a)-(45c) are the kinematic equations of the lower limb with 4 rotary joints,

$$x=l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{123}$$

(45a)

and

$$y=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123})$$

(45b)

and

$$z=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123}).$$

(45c)

The instantaneous Jacobian matrix represents the rate of change of coordinates $xyz$ w.r.t. the joint-space vector $\mathsf{\Phi }={({\varphi }_0,{\varphi }_1,{\varphi }_2,\cdots {\varphi }_n)}^{\top }$ of length n. Such that,

$$\mathbf{J}\left(\mathsf{\Phi }\right)=\left(\begin{array}{ccccc}\frac{\partial x}{\partial {\varphi }_0}& \frac{\partial x}{\partial {\varphi }_1}& \frac{\partial x}{\partial {\varphi }_2}& \cdots & \frac{\partial x}{\partial {\varphi }_n}\\ \frac{\partial y}{\partial {\varphi }_0}& \frac{\partial y}{\partial {\varphi }_1}& \frac{\partial y}{\partial {\varphi }_2}& \cdots & \frac{\partial y}{\partial {\varphi }_n}\\ \frac{\partial z}{\partial {\varphi }_0}& \frac{\partial z}{\partial {\varphi }_1}& \frac{\partial z}{\partial {\varphi }_2}& \cdots & \frac{\partial z}{\partial {\varphi }_n}\end{array}\right).$$

(46)

For the upper limb the Jacobian matrix elements are

$$\frac{\partial x}{\partial {\varphi }_0}=0$$

(47)

$$\frac{\partial x}{\partial {\varphi }_1}=l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}$$

(48)

$$\frac{\partial x}{\partial {\varphi }_2}=l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}$$

(49)

$$\frac{\partial x}{\partial {\varphi }_3}=-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}$$

(50)

$$\frac{\partial x}{\partial {\varphi }_4}=l_3{\mathrm{s}}_3{\mathrm{s}}_{124},$$

(51)

and for y,

$$\frac{\partial y}{\partial {\varphi }_0}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(52)

$$\frac{\partial y}{\partial {\varphi }_1}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(53)

$$\frac{\partial y}{\partial {\varphi }_2}={\mathrm{s}}_0(l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(54)

$$\frac{\partial y}{\partial {\varphi }_3}=l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(55)

$$\frac{\partial y}{\partial {\varphi }_4}=l_3{\mathrm{s}}_{124}{\mathrm{s}}_0+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0,$$

(56)

and for z,

$$\frac{\partial z}{\partial {\varphi }_0}=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})+l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(57)

$$\frac{\partial z}{\partial {\varphi }_1}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(58)

$$\frac{\partial z}{\partial {\varphi }_2}={\mathrm{s}}_0(l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{124})-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(59)

$$\frac{\partial z}{\partial {\varphi }_3}=-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(60)

$$\frac{\partial z}{\partial {\varphi }_4}=l_3{\mathrm{s}}_{124}{\mathrm{s}}_0-l_3{\mathrm{s}}_3{\mathrm{s}}_{124}{\mathrm{s}}_0$$

(61)

Likewise, for the lower limb, the Jacobian matrix terms are

$$\frac{\partial x}{\partial {\varphi }_0}=0;\frac{\partial x}{\partial {\varphi }_1}=l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123};$$

(62)

$$\frac{\partial x}{\partial {\varphi }_2}=l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123};\frac{\partial x}{\partial {\varphi }_3}=-l_3{\mathrm{s}}_{123}.$$

(63)

and for y

$$\frac{\partial y}{\partial {\varphi }_0}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123})$$

(64)

$$\frac{\partial y}{\partial {\varphi }_1}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{123})$$

(65)

$$\frac{\partial y}{\partial {\varphi }_2}={\mathrm{s}}_0(l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{123});\frac{\partial y}{\partial {\varphi }_3}=l_3{\mathrm{s}}_0{\mathrm{s}}_{123}.$$

(66)

and for z

$$\frac{\partial z}{\partial {\varphi }_0}=-{\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}-l_3{\mathrm{s}}_{123})$$

(67)

$$\frac{\partial z}{\partial {\varphi }_1}={\mathrm{s}}_0(l_1{\mathrm{s}}_1+l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{123})$$

(68)

$$\frac{\partial z}{\partial {\varphi }_2}={\mathrm{s}}_0(l_2{\mathrm{s}}_{12}+l_3{\mathrm{s}}_{123});\frac{\partial z}{\partial {\varphi }_3}=l_3{\mathrm{s}}_0{\mathrm{s}}_{123}$$

(69)

4.2. Recursive time-variant control law

Other exoskeleton control approaches have addressed nonlinearity using sliding surfaces [29] or adaptive methods, implementing impedance with feedback using biosignals [30], and neural approximation [31].

The present work takes an approach involving a model-based algebraic linear recursive system of equations. This work deduces a time-variant redundant kinematics controller since ${\mathbf{J}}_t$ is considered. The model incorporates the reference trajectory provided by the therapist, denoted as ${\mathbf{p}}^{ref}$. The real-time feedback utilizes the EMG observation ${\widehat{\mathbf{p}}}_t={({\widehat{x}}_t,{\widehat{y}}_t)}^{\top }$. The innovation model, representing the local reference, is the inverse model ${\mathsf{\Phi }}_{t+1}$, and the updating model involves the encoders’ observation ${\widehat{\mathsf{\Phi }}}_t$.

Let us define the state linear equation (70) as the first-order derivative forward kinematics.

$$\dot{\mathbf{p}}=\mathbf{J}·\dot{\mathsf{\Phi }}$$

(70)

and stating this expression in the form of differential equation, its derivative-order decreases

$$\frac{d\mathbf{p}}{dt}=\mathbf{J}·\frac{\mathsf{\Phi }}{dt},$$

(71)

completing all differentials with their respective integrals at both sides of the equality

$${\int }_{{\mathbf{p}}_1}^{{\mathbf{p}}_2}d\mathbf{p}=\mathbf{J}·{\int }_{{\mathsf{\Phi }}_1}^{{\mathsf{\Phi }}_2}d\mathsf{\Phi },$$

(72)

therefore, the expression (73) is the general recursive solution,

$${\mathbf{p}}_{t2}-{\mathbf{p}}_{t1}=\mathbf{J}·({\mathsf{\Phi }}_{t2}-{\mathsf{\Phi }}_{t1}).$$

(73)

Therefore, from such an equality the expressions (74) lead to forward/backward general solutions.

The recursive forward:

$${\mathbf{p}}_{t+1}={\mathbf{p}}_t+{\mathbf{J}}_t·({\mathsf{\Phi }}_{t+1}-{\mathsf{\Phi }}_t)$$

(74a)

and the recursive backward:

$${\mathsf{\Phi }}_{t+1}={\mathsf{\Phi }}_t+{\mathbf{J}}_t^{+}·({\mathbf{p}}_f-{\mathbf{p}}_t).$$

(74b)

Therefore, the inverse general model for the recursive control law is set by an ideal global trajectory ${\mathbf{p}}_t^{ref}$,

$${\mathsf{\Phi }}_{t+1}={\mathsf{\Phi }}_t+{\mathbf{J}}_t^{+}\left({\mathbf{p}}_t^{ref}-{\widehat{\mathbf{p}}}_t\right).$$

(75)

Some terms in the model are expanded. Since the Jacobian is non-squared, a pseudoinverse matrix solution is applied,

$${\mathsf{\Phi }}_{t+1}={\mathsf{\Phi }}_t+{\left({\mathbf{J}}_t{\mathbf{J}}_t^{\top }\right)}^{-1}{\mathbf{J}}_t^{\top }\left[\left(\begin{array}{c}{x}_t^{ref}\\ y_t^{ref}\end{array}\right)-\left(\begin{array}{c}\widehat{x}\left(\widehat{t}\right)\\ \widehat{y}\left(\widehat{t}\right)\end{array}\right)\right]$$

(76)

and the forward solution has the next local reference ${\mathsf{\Phi }}_{t+1}$,

$${\mathbf{p}}_{t+1}={\mathbf{p}}_t+{\mathbf{J}}_t\left({\mathsf{\Phi }}_{t+1}-{\widehat{\mathsf{\Phi }}}_t\right).$$

(77)

Where the exoskeleton’s joints direct measurements are obtained from encoder devices,

$${\mathbf{p}}_{t+1}={\mathbf{p}}_t+{\mathbf{J}}_t\left({\mathsf{\Phi }}_{t+1}-\frac{2\pi }{R}\left(\begin{array}{c}{\eta }_0\\ {\eta }_1\\ ⋮\\ {\eta }_k\end{array}\right)\right).$$

(78)

When the error condition $\parallel {\mathbf{p}}_f-{\mathbf{p}}_t\parallel >ϵ$ is met, the control loop achieves goal tracking. Figure 19a illustrates a complex trajectory and its tracking control performance concerning the joint angles (Figure 19b). Additionally, Figure 19c depicts the evolution of position metric errors throughout the entire reference trajectory.

In the context of numerical simulations for numerous controlled trajectories, the efficacy of the control approach was validated through a comparative analysis against a reference ideal function. The trajectory in Figure 20-above illustrates an $xy$ path with an arbitrary reference defined as $y\left(x\right)=k_5{x}^5+k_3{x}^3-k_0$. Along the same row of the figure, the plots elucidate the error between the reference function and the controlled function. Ideally, a precise match is anticipated between the two functions. The observed errors in the x and y components, in general, are relatively small. As for the z component, its error is nearly zero since this trajectory is intentionally planar, confined solely to the $xy$ plane. Consequently, any variation along the z axis should be virtually null. This illustration underscores the expectation that the controlled functions accurately align with a segment of a 5^th-grade polynomial function.

Moreover, Figure 20-below displays the resulting errors in joints’ angles induced by the trajectory in the same figure (expressed in physical units as radians).

A comprehensive depiction akin to Figure 20, showcasing the entire joint performance, is presented in Figure 21. A distinct range of coordinates from a different reference function is established as the reference trajectory, with the framework’s origin situated at the shoulder centroid. Similarly, Figure 21-above exhibits a Cartesian plot along with the errors along the three axes (reference versus controlled function). Figure 21-below illustrates the joints’ angle performance resulting from the controlled tracking trajectory.

5. Results and discussion

Summarizing the core focus of this work, we direct our attention to the following three figures: Figure 22, Figure 23, and Figure 24. These figures encapsulate the fundamental purpose and achievement of this study. They provide a visual representation of the work’s primary objective, which is the estimation of joint angles from EMG measurements. For the sake of clarity and simplicity, the plots in these figures display data for two key joint angles, highlighting the remarkable success of the approach. Figure 22 presents data from three distinct experiments involving elbow flexion. Within this figure, the shoulder angle is denoted as ${\varphi }_1$ (depicted in red), and the elbow angle is represented as ${\varphi }_2$ (in blue). Two specific channels, Ch5 and Ch6, correspond to the electromyography measurements of the biceps brachii and triceps brachii, respectively. The selection of these two channels reflects the authors’ recognition of their suitability for assessing electromyography activity during limb motion.

Following the calibration process, the results in Figure 22 demonstrate that the activity associated with elbow flexion (depicted in blue) reveals a strong correlation between the elevation angle and the measurements from Ch5, emphasizing the crucial role of the biceps brachii in this motion. Conversely, the decrease in the elbow’s angle is found to be correlated with the measurements from Ch6, highlighting the impact of the triceps brachii. Interestingly, the motion of the shoulder joint (${\varphi }_1$) is not detected within this specific configuration, as indicated by the red-colored plot. This observation sheds light on the distinct contributions of the two channels to the control of different joint angles.

Figure 22. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in elbow-flexion performance.

Figure 22. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in elbow-flexion performance.

Figure 23 offers a close look at an experiment comprising five repetitions of shoulder flexion. Within this figure, the shoulder angle continues to be represented as ${\varphi }_1$, depicted in red, while the elbow angle is identified as ${\varphi }_2$ and shown in blue. The electromyography measurements are captured by channels Ch5 and Ch6, corresponding to the activity of the pectoralis major and anterior pectoralis muscles, respectively. Following the calibration process, the results in Figure 23 reveal that during shoulder flexion activity (as depicted in the red-colored plot), the elevation angle is strongly influenced, with substantial contributions from both Ch6 and Ch5 measurements. This highlights the role of both the anterior pectoralis and the pectoralis major muscles in controlling the shoulder’s elevation angle.

However, as the exercise primarily involves shoulder motion and maintains the elbow at a fixed position, the motion of joint ${\varphi }_2$ is not detected in this particular configuration, as indicated by the blue-colored plot. This distinction further underscores the remarkable accuracy of the approach in distinguishing the contributions of different muscle groups to specific joint movements. Overall, Figure 23 demonstrates the successful application of electromyography measurements in estimating joint angles, exemplifying the adaptability and precision of this approach for diverse exercises involving multiple joints.

Figure 23. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in shoulder-flexion performance.

Figure 23. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in shoulder-flexion performance.

Figure 24 offers insight into three distinct experiments involving rowing-like flexion, wherein both shoulder and elbow joints exhibit synchronized motion activity. As in previous figures, the shoulder angle is denoted as ${\varphi }_1$, represented by the red-colored plot, while the elbow angle is identified as ${\varphi }_2$ and displayed in blue. The electromyography measurements for these experiments are derived from channels Ch5 and Ch6, corresponding to the activity of the middle trapezius and anterior deltoid muscles, respectively. After calibration, the results within Figure 24 reveal that elbow flexion activity (illustrated by the blue-colored plot) exhibits different oscillatory angular phases. The data suggests that the elevation of ${\varphi }_1$ appears to be significantly impacted by Ch5 measurements, while ${\varphi }_2$ is influenced by Ch6 data. The remarkable consistency observed across these three experiments serves to validate the proposed approach in this research, showcasing the method’s ability to accurately estimate joint angles during complex exercises that involve synchronized motion of multiple joints. This validation is a testament to the reliability and effectiveness of the approach, underscoring its potential for real-world applications in the field of tele-rehabilitation.

Figure 24. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in rowing-like flexions performance.

Figure 24. Calibrated EMG measurements, estimation joints’ angle ${\varphi }_1$ and ${\varphi }_2$ in rowing-like flexions performance.

Table 3. Physics-based computer simulation parameters^a.

Limbs:
m=1 kg	$I_x=0.1666$ kg ${\mathrm{m}}^2$	$I_y=0.1666$ kg ${\mathrm{m}}^2$	$I_z=0.1666$ kg ${\mathrm{m}}^2$
${ł}_{1,2}=0.3$ m	limbs $r_{1,2}=0.15$ m	${\tau }_{ref}=19.62$ Nm	${\ddot{\varphi }}_{ref}=36.14$$s^{-2}$
Gears:
$r_1=1.5r_0$ m	$r_{2,3}=2.6r_0$ m	$r_{4,5}=2r_0$ m	$r_{6,7}=0.66r_0$ m
$r_8=1.16r_0$ m	$n_{6,7}=2$ (threads)	$n_8=40$ (teeth)	$m_8=1.166m_0$ kg
$m_1=1.5m_0$ kg	$m_{2,3}=3m_0$ kg	$m_{4,5}=2m_0$ kg	$m_{6,7}=1.5m_0$ kg
$I_1=3.37r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{2,3}=1.72r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{4,5}=8r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{6,7}=0.2875r_0^2{m}_0$ kg ${\mathrm{m}}^2$
$I_8=1.58r_0^2{m}_0$ kg ${\mathrm{m}}^2$
aGazebo enabled with Object Dynamic Object (ODE).

Table 3. Physics-based computer simulation parameters^a.

Limbs:
m=1 kg	$I_x=0.1666$ kg ${\mathrm{m}}^2$	$I_y=0.1666$ kg ${\mathrm{m}}^2$	$I_z=0.1666$ kg ${\mathrm{m}}^2$
${ł}_{1,2}=0.3$ m	limbs $r_{1,2}=0.15$ m	${\tau }_{ref}=19.62$ Nm	${\ddot{\varphi }}_{ref}=36.14$$s^{-2}$
Gears:
$r_1=1.5r_0$ m	$r_{2,3}=2.6r_0$ m	$r_{4,5}=2r_0$ m	$r_{6,7}=0.66r_0$ m
$r_8=1.16r_0$ m	$n_{6,7}=2$ (threads)	$n_8=40$ (teeth)	$m_8=1.166m_0$ kg
$m_1=1.5m_0$ kg	$m_{2,3}=3m_0$ kg	$m_{4,5}=2m_0$ kg	$m_{6,7}=1.5m_0$ kg
$I_1=3.37r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{2,3}=1.72r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{4,5}=8r_0^2{m}_0$ kg ${\mathrm{m}}^2$	$I_{6,7}=0.2875r_0^2{m}_0$ kg ${\mathrm{m}}^2$
$I_8=1.58r_0^2{m}_0$ kg ${\mathrm{m}}^2$
aGazebo enabled with Object Dynamic Object (ODE).

Although this study does not address more intricate rehabilitation types, the showcased exercises are representative of typical rehabilitation routines, serving as a testament to the feasibility of the proposed method. Currently, the demonstration of the model’s applicability is achieved using only two EMG channels. Given the limitation of relying on a singular visual sensor, the integration of various technological sensing devices becomes imperative. This broader array of sensing devices aims to encompass diverse perspectives surrounding the therapist’s tasks, thereby enhancing the overall accuracy of measurements. In the context of trajectory adjustment for the linear trajectory smoother model, a manual trial-and-error approach is employed to fine-tune the trajectory parameters, specifically $k_t$ and $k_s$.

6. Conclusion

In this comprehensive research endeavor, a remarkable paradigm shift was introduced in the realm of tele-rehabilitation. The proposed approach has a feature: allowed for asynchronous interactions, removing the constraints of rigid scheduling, as exercises could be conducted at different dates and times, making rehabilitation more accessible and flexible for patients. The patient-exoskeleton interface was engineered to receive concise data structures, and the innovative use of coefficient sets emerged as a masterstroke. These sets proved invaluable in ensuring that data streaming remained swift and uninterrupted, unfazed by the passage of time. One of the standout features of this work was the treatment of all exercises as periodic functions. This approach enabled the seamless application of Fourier series fitting, providing an elegant and efficient means of data representation. The ability to translate complex exercises into Fourier series representations presented a new dimension of flexibility in tracking and control.

A cornerstone of this research was the method for reproducing rehabilitation exercises using Cartesian analytical polynomials. This method unlocked the door to producing high-resolution functions of time by synthesizing Cartesian functions of time with EMG data, with time itself serving as the universal synchronization variable. This synthesis represented the union of clinical precision and technological sophistication, making the rehabilitation process more accurate and efficient. A novel model for inferring Cartesian posture from EMG signals was an important issue of this work, showcasing remarkable reliability. The need for recalibration with each new exercise served as a testament to the adaptability and responsiveness of this approach, ensuring that the exoskeleton remained in sync with the patient’s unique characteristics and movements.

The proposed exoskeleton underwent rigorous simulation testing, encompassing hundreds of tracking control rehabilitation exercises. Its adaptability was underscored by a control law that could be seamlessly adjusted to suit various exoskeleton structures, simply by readjusting the Jacobian matrix. This adaptability introduced an exciting degree of versatility and expandability in the application of this technology, not limited by specific hardware constraints. In essence, this simulation prototype, anchored by a single reference trajectory, demonstrated remarkable potential in the calibration of numerous EMG experiments and the execution of multiple tracking control exercises. The system’s adaptability, efficiency, and the profound impact it could have on the field of tele-rehabilitation were evident throughout this research. The approach outlined in this study promises to redefine and enhance the way rehabilitation is conducted, making it not just more effective but also more accessible and accommodating to a wide range of patient needs. It represents a significant step forward in the integration of technology and healthcare, and its implications are bound to be felt for years to come.

As a future work, enhancing the visual teaching process can be accomplished by substituting the RGB camera with either an RGB-D or a ranging laser device. This alternative not only mitigates noise but is also impervious to constraints and illumination variations, offering the added advantages of infrared capabilities and 3D cloud data provision. In addition, to ensure the linear trajectory smoother model’s optimal performance, automatic optimization methods and tuning of rehabilitation trajectories, facilitated by parameters $k_t$ and $k_s$, will enable a systematic convergence toward the desired reference trajectory shape. In our forthcoming work, there is an envisioned redesign of the exoskeleton mechanism. This redesign incorporates elastic soft materials, mimicking artificial muscles activated by electrical/optical signals, preserving EMG feedback. In addition, addressing the expansion of the EMG device array is pivotal, allowing for the comprehensive mapping of additional body muscles while utilizing more compact technological devices. We posit that exploring alternative controllers is imperative for achieving enhanced performance, tailored to different rehabilitation tasks, fostering a more nuanced understanding of optimal controller-task relationships.