Scalable O(log2n) Dynamics Control for Soft Exoskeletons

1. Introduction

Exoskeletons are emerging as effective tools for facilitating motor rehabilitation in patients with impaired upper limb function [1,2]. In this regard, robotic exoskeletons designed specifically for the hand have the potential to enhance hand rehabilitation and assist individuals in activities of daily live (ADLs). The human hand poses a significant challenge for the design and control of such devices, due to its high bio-mechanical complexity [3,4]. The hand can perform precise manipulation tasks, thanks to the coordinated action of 27 bones, 27 joints, and 34 intrinsic and extrinsic muscles.

Hand exoskeletons commonly mimic the skeletal structure of the fingers, the phalanges, by means of joint-based rigid-bodies. Recently, the use of soft materials has been proposed, as they allow for more natural movements and better adaptation to the hand dexterity. The main efforts on soft exoskeletons have focused on deformable materials and real-time control of an increased number of degrees of freedom (DoF) [5,6,7,8]. However, robust closed-loop control strategies are difficult to apply since most of these methods require the calculation of the exoskeleton dynamics to compute the control signal in real-time. Among the control models for exoskeletons, position-based and force-based controls are widely used. Additionally, control strategies based on the evaluation of electromyographic (EMG) signals, combined with machine learning algorithms are increasingly proposed [9,10,11]. More recently, electroencephalographic (EEG) signals and cameras have also been proposed for motor control [12]. For instance, in [13], a rigid-soft combined mechanism is developed to support hand rehabilitation based on a EEG brain-controlled switch implemented in Matlab™, while an external VICON™motion capturing system is required to determine the fingers position in 3D space. The authors in [14] present a force myography (FMG) to allow grasping assistance. The controller is based on a switching algorithm to properly modulate a soft-extra-muscle (SEM) glove, while machine-learning models running offline were trained to estimate grasp forces. Approaches in [11,15,16] depend on machine-learning offline processing to support EEG/EMG control, force-based control or simple switch-driven modulation.

In terms of actuation, pneumatics actuators (PAs) [17], electroactive polymers (EAP) [18] and shape memory alloys (SMAs) [19] are the most predominant technologies used in soft robotic exoskeleton and gloves. As mentioned in [11,18], in most of the cases, finite elements methods (FEMs) are used to model flexible and deformable properties, requiring compute-intensive requirements. In previous efforts reported in [20,21], we presented the development of a rigid-body robotic exoskeleton to assist in the recovery of hand motor function for post-stroke subjects. Our study proposed a closed-loop model predictive control (MPC) method with the aim of compensating the effects of muscle fatigue during the rehabilitation therapies. To this purpose, EMG-based muscular effort must be detected within the control loop in order to modulate the exoskeleton’s velocity accordingly. Since the MPC method requires the computation of the exoskeleton dynamics model, applying this control algorithm implies significant challenges to soft exoskeleton mechanisms. In fact, few works in the literature [22] have applied robust control methods to modulate a soft robotic exoskeleton.

In this paper, we tackle the challenge of calculating the exoskeleton’s equations of motion (EoM) to support robust dynamics control for a hybrid soft-exoskeleton design based on compliant joints, as depicted in the inset of Figure 1a. We have adopted the parallel structure introduced by [23] to solve large-scale EoM dynamics as a foundation of the proposed soft-based highly-redundant joint exoskeleton robot. The proposed dynamics control method via inverse-dynamics EoM is based on the well-known Recursive Newton-Euler Algorithm (RNEA), presented by [24]. Figure 1b details the multi-body diagram to represent each digit of the exoskeleton. Spatial operators for velocity ($V_i$), acceleration (${\dot{V}}_i$) and force ($F_i$) led to six-dimensional physical quantities that combined the angular and linear motions of the system. The RNEA method computes two recurrences to propagate the aforementioned spatial operators from the exoskeleton’s base to each fingertip of the serial-chain. In this regard, solving the inverse-dynamics via the RNEA method allows for the calculation of applied forces to properly control the exoskeleton’s dynamics motion. Other works in [25,26,27] have used similar approaches to dynamics modeling and robust control for exoskeletons based on the RNEA method. Nonetheless, to the best of our knowledge, this is the first work to apply a parallel solution for the RNA as an alternative to control soft exoskeletons in real-time.

As mentioned, classical techniques used in soft-robotics that are mostly based on finite element theory are not suitable in our application, since we require online computation of the dynamics model within a closed-loop control algorithm capable to scale up. To this purpose, hardware-based parallel computation of the EoM will be processed within an embedded system with multi-core CPU/GPU capabilities, to support the real-time calculation of the exoskeleton’s large-scale dynamics model, considering up to $n=512$ degrees of freedom to emulate continuous flexible soft bending, and achieving a $O\left(lo{g}_{2}n\right)$ computational complexity. To support the evaluation of this approach, we propose a hardware-in-the-loop (HIL) architecture as represented in Figure 1a. This control scheme consists of two main components: (i) a software-side modeling the soft-exoskeleton under analysis, which is supported by Matlab®, directly running on a desktop computer; and (ii) a hardware-side supported by an embedded platform with highly parallel computing capabilities, in charge of calculating the inverse dynamics control to drive the exoskeleton.

By adopting this HIL-based architecture, it is possible to represent very complex models, such as the referenced soft-exoskeleton, without the necessity of actually building it and integrating all the required sensors and acquisition stages. Nonetheless, this approach also allows us to validate the proposed control scheme in real hardware platforms, such as the selected embedded systems. The software-side reproduces all the complex dynamics of the exoskeleton and generates the corresponding outputs (position, velocity and acceleration) that are sent to the hardware-side. With these signals coming from the exoskeleton’s digital twin, the hardware-side applies the corresponding control scheme and in turn generates the control signals to the virtual actuators (e.g., motors), also implemented in the software-side. As it has been explained before, the biggest advantage is the ability to test, validate and evaluate the performance and scalability (with respect to the number of available cores or processing units) of the control scheme as it would execute in a real deployment. Our main contributions can be summarized as follows:

• an alternative approach to calculate dynamics equations of motion for soft exoskeletons in real-time, adopting a parallel computation for highly-articulated systems with a $O\left(lo{g}_{2}n\right)$ computational complexity, and
• the integration and deployment of this approach on embedded systems, facilitating modular dynamics-based control of soft exoskeletons by following a hardware-in-the-loop (HIL) approach.

The rest of the paper is organized as follows: Section 2.1 briefly introduces several design approaches for the proposed soft exoskeleton, based on previous works reported in [20,21,28]. Also, the formulation for the equations of motion (EoM) based on the Newton-Euler formalism using spatial operators; Section 2.2 details the proposed parallel computing approach for the EoM, by applying the well-known RNEA method (Recursive Newton-Euler Algorithm) [24,25], allowing dynamics control of the soft exoskeleton with a computational performance $O\left(lo{g}_{2}n\right)$; Section 3 presents the computing results obtained from several multi-core CPU/GPU embedded systems running under a HIL approach; Section 4 deals with discussions and remarks about the proposed approach for controlling soft exoskeletons; and Section 5 concludes this work with some insights, limitations and future work.

2. Methods

Our first design approach for the exoskeleton was presented in previous works reported in [20,21]. The system was conceived as a branched physical model, in which each finger is composed by a chain of rigid bodies serially connected through joint links. Figure 2a details the former mechanical design. As mentioned, we are working on a novel design approach supported by a hybrid soft mechanism based on compliant joints, as shown in Figure 2c. In order to apply our proposed method to model such mechanism, we require to determine a sufficient amount of joints to emulate soft motion properties. Figure 2b details this approximation to a soft finger of the exoskeleton via 512 degrees of freedom.

2.1. Equations of Motion

Soft motion is described based on Newton-Euler dynamics extended to large-scale articulated systems. Equations of Motion (EoM) are presented using spatial operators in ${\mathbb{R}}^6$, allowing joint kinematics and dynamics calculations to mimic soft deformation. Next subsection reviews the fundamental concepts of the classic mechanics that are related to body dynamics modeling, establishing an appropriate mathematical representation of the physical quantities involved in the dynamics behavior of the robotic exoskeleton.

Assuming from Figure 1b that the joint frame i and the center of mass ($cm$) are two points located on the body i, the term ${\overrightarrow{s}}_{i,cm}\in {\mathbb{R}}^3$ is the position vector connecting frames i and $cm$, while ${\overrightarrow{p}}_{i,i+1}\in {\mathbb{R}}^3$ connects consequently joint frames that are serially coupled. Also, the term $r_{i,i-1}\in {\mathbb{R}}^{3\times 3}$ is the rotation matrix that relates frames $\left\{i\right\}$ with $\{i-1\}$. The translational ($v_i$) and angular velocities (${\omega }_i$) for each body, as well as and forces ($f_i$) and moments (${\tau }_i$) acting on the joint frame can be integrated into spatial operators that led to six-dimensional physical quantities, by combining the angular and linear variables of body motions. Using spatial notation, velocities V, accelerations $\dot{V}$, and forces F of a body i are expressed with respect to center of mas $\{i,cm\}$ using six-dimensional vectors in ${\mathbb{R}}^6$, as follows:

$$V_{i,cm}\equiv \left[\begin{array}{c}{w}_{i,cm}\\ v_{i,cm}\end{array}\right],{\dot{V}}_{i,cm}\equiv \left[\begin{array}{c}{\dot{w}}_{i,cm}\\ {\dot{v}}_{i,cm}\end{array}\right],F_{i,cm}\equiv \left[\begin{array}{c}{\tau }_{i,cm}\\ f_{i,cm}\end{array}\right]$$

(1)

Using the position vector ${\overrightarrow{s}}_{i,cm}\in {\mathbb{R}}^3$, shown in Figure 1b, the physical quantities can be also expressed with respect to the joint frame $\left\{i\right\}$, as:

$$\begin{array}{c}{V}_{i,cm}=S_{i,cm}^T{V}_i\hfill \\ F_{i,cm}=S_{i,cm}^T{F}_i\hfill \end{array}$$

(2)

Being $S_{i,cm}\in {\mathbb{R}}^{6\times 6}$ the spatial operator connecting both frames, defined by the skew-symmetric matrix ${\tilde{s}}_{i,cm}\equiv {\overrightarrow{s}}_{i,cm}$, as:

$$\begin{array}{c}{S}_{i,cm}=\left[\begin{array}{cc}U& {\tilde{s}}_{i,cm}\\ 0& U\end{array}\right]\hfill \end{array}$$

(3)

Where $U\in {\mathbb{R}}^{3\times 3}$ is the identity operator. From Equation (2), the spatial forces acting onto the joint frame ($F_i$) can be derived as $F_{i,cm}={\displaystyle \frac{d}{dt}}\left(I_{i,cm}{V}_{i,cm}\right)$, where $I_{i,cm}\in {\mathbb{R}}^{6\times 6}$ is the spatial inertia operator of the body i with respect to its center of mass ($cm$), and $V_{i,cm}\in {\mathbb{R}}^{6\times 1}$ corresponds to the spatial velocity defined in Equation (2). Differentiating as a function of time yields:

$$\begin{array}{c}{F}_{i,cm}={\displaystyle \frac{d}{dt}}\left(I_{i,cm}{V}_{i,cm}\right)\hfill \\ F_{i,cm}=I_{i,cm}{\dot{V}}_{i,cm}+{\dot{I}}_{i,cm}{V}_{i,cm}\hfill \end{array}$$

(4)

Using the spatial operation for translation defined in Equation (3), spatial forces with respect to the joint frame can be calculated as:

$$\begin{array}{c}{F}_i=S_{i,cm}{F}_{i,cm}\hfill \\ =S_{i,cm}[I_{i,cm}{\dot{V}}_{i,cm}+{\dot{I}}_{i,cm}{V}_{i,cm}]\hfill \end{array}$$

(5)

From Equation (2), spatial accelerations can be defined as ${\dot{V}}_{i,cm}=S_{i,cm}^T{\dot{V}}_i+{\dot{S}}_{i,cm}^T{V}_i$. Replacing into Equation (5) yields:

$$\begin{array}{c}{F}_i=S_{i,cm}[I_{i,cm}(S_{i,cm}^T{\dot{V}}_i+{\dot{S}}_{i,cm}^T{V}_i)+{\dot{I}}_{i,cm}{V}_{i,cm}]\hfill \\ =S_{i,cm}{I}_{i,cm}{S}_{i,cm}^T{\dot{V}}_i+S_{i,cm}{I}_{i,cm}{\dot{S}}_{i,cm}^T{V}_i+S_{i,cm}{\dot{I}}_{i,cm}{V}_{i,cm}\hfill \\ =I_i{\dot{V}}_i+I_i{\dot{S}}_{i,cm}{V}_i+{\dot{I}}_i{V}_i\hfill \\ =I_i{\dot{V}}_i+[I_i{\dot{S}}_{i,cm}+{\dot{I}}_i]V_i\hfill \end{array}$$

(6)

Since each finger of the exoskeleton is composed by a serial chain of connected joints, spatial forces can be propagated from the fingertip ($i=n$) to the base ($i=1$), by using specific operators for rotation and translation. Therefore, a propagated term can be included into Equation (6), such as:

$$\begin{array}{c}\sum _{i=n}^1{F}_i=\underset{propagated}{\underbrace{\left(P_{i,i+1}^T{R}_{i,i-1}\right)F_{i+1}}}+I_i{\dot{V}}_i+[I_i{\dot{S}}_{i,cm}+{\dot{I}}_i]V_i\hfill \end{array}$$

(7)

Equation (7) describes a backward propagation of spatial forces, where six-dimensional operators in ${\mathbb{R}}^{6\times 6}$ for rotation R and translation P are defined as follows:

$$\begin{array}{c}{R}_{i,i-1}=\left[\begin{array}{cc}{r}_{i,i-1}& 0\\ 0& r_{i,i-1}\end{array}\right],P_{i,i+1}=\left[\begin{array}{cc}U& {\tilde{p}}_{i,i+1}\\ 0& U\end{array}\right]\hfill \end{array}$$

(8)

Where $r_{i,i-1}\in {\mathbb{R}}^{3\times 3}$ is the classical basic rotation matrix and ${\tilde{p}}_{i,i+1}\in {\mathbb{R}}^{3\times 3}$ is the skew-symmetric matrix of the position vector ${\overrightarrow{p}}_{i,i+1}$. As depicted in Figure 1b, spatial forces are calculated by following a backward propagation, initiating at the fingertip frame ($i=n$) and finalizing at the base frame of the exoskeleton’s finger ($i=1$). To this purpose, the spatial terms $P_{i,i+1}^T{R}_{i,i-1}$ in Equation (4) allow the propagation of the spatial forces (F) through the exoskeleton, considering the external load ($F_{load}$) imposed by the hand or any manipulated object. Therefore, $F_{i+1}\equiv F_{load},\forall i=n$.

Applying a similar approach, both $R_{i,i-1}$ and $P_{i,i+1}$ are also used to propagate spatial velocities ($V_i$) and accelerations (${\dot{V}}_i$) by following a forward propagation from $i=1...n$, as:

$$\begin{array}{c}\sum _{i=1}^n{V}_i=\underset{propagated}{\underbrace{R_{i,i-1}^T{P}_{i,i+1}{V}_{i-1}}}+H{\dot{q}}_i\hfill \\ \\ \sum _{i=1}^n{\dot{V}}_i=\underset{propagated}{\underbrace{R_{i,i-1}^T{P}_{i,i+1}{\dot{V}}_{i-1}+{\dot{P}}_{i,i+1}{\dot{R}}_{i,i-1}^T{V}_{i-1}}}+H{\ddot{q}}_i+\dot{H}{\dot{q}}_i\hfill \end{array}$$

(9)

The mathematical expressions for the Equations of Motion (EoM): spatial velocities $V_i$, accelerations ${\dot{V}}_i$ and forces $F_i$ are presented in Equation (7) and (9) respectively, and consolidated in Algorithm 1. The inputs to the algorithm are the joint trajectories (therapy motion) determined by the terms $q_i,{\dot{q}}_i,{\ddot{q}}_i$, whereas the output is the required torque ${\tau }_i$ for each joint i. As mentioned, the method is divided into two recurrences: (i) a forward calculation and propagation of the spatial velocities and accelerations and (ii) a backward recurrence to calculate spatial forces and the torques acting on each joint of the fingers. Table 1 presents the list of parameters used for the EoM. Nonetheless, this $O\left(n\right)$ approach is insufficient to model large-scale soft systems composed by thousands of joints that enable body deformation, since the computing time increases linearly as a function of the number of degrees of freedom n.

Algorithm 1 $O\left(n\right)$ RNE serial computation

1. Input joint trajectory (therapy motion) for each joint i: $q_i,{\dot{q}}_i,{\ddot{q}}_i$   2. initialize parameters: fixed-base (no 6D motion of the wrist) and gravity acceleration acting on the carpal (along x-axis).   3. $V_0\leftarrow {\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right]}^T$,    4. ${\dot{V}}_0\leftarrow {\left[\begin{array}{cccccc}0& 0& 0& -9.81& 0& 0\end{array}\right]}^T$   5. Forward recurrence: spatial velocities and accelerations propagation from carpals to distal phalanges:   for$i=1\to n$do    (1) $V_i\leftarrow R_{i,i-1}^T{P}_{i,i+1}{V}_{i-1}+H{\dot{q}}_i$      (2) ${\dot{V}}_i\leftarrow R_{i,i-1}^T{P}_{i,i+1}{\dot{V}}_{i-1}+{\dot{P}}_{i,i+1}{\dot{R}}_{i+1,i}^T{V}_{i-1}+H{\ddot{q}}_i+\dot{H}{\dot{q}}_i$   end for   ——————————————————-    6. $if$ (payload at the distal phalange is included) do$F_{i+1}\leftarrow$ payload force $F_{load}$   $else$do$F_{i+1}\leftarrow {\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right]}^T$    7. Backward recurrence: spatial forces calculated from distal phalanges to the carpals.   for$i=n\to 1$do       (3) $F_i\leftarrow P_{i,i+1}^T{R}_{i,i-1}{F}_{i+1}+I_i{\dot{V}}_i+[I_i{\dot{S}}_{i,cm}+{\dot{I}}_i]V_i$    (4) ${\tau }_i\leftarrow H^T{F}_i$    end for   Return $\tau$

2.2. $O\left(lo{g}_{2}n\right)$ Parallel Computing

Embedded systems with multi-core CPU/GPU architectures has enabled an increase in computational power that allows complex algorithms to run on the edge in real-time. This novel edge computing paradigm can be applied to our system, allowing a modular and compact design for the exoskeleton’s control electronics. In this regard, we propose a reformulation of Algorithm 1 to compute the EoM in a parallel fashion.

As previously mentioned, the classical formulation of the EoM based on the spatial Newton Euler formalism requires an $O\left(n\right)$ computational complexity. It was [30] that improved the computational efficiency by computing the EoM recursively (forward dynamics) by referring the forces to generalized coordinates in the joint space. This was called the Recursive Newton-Euler Algorithm (RNEA), which was later modified by [31], by using a new approach called CRBA (Composed Rigid Body Algorithm). In [32] presented the first parallel approach to solve the CRBA for both inverse/foward dynamics problems with an $O\left(lo{g}_{2}n\right)$ complexity, but strictly dependent on the overlapping costs associated to calculation of the inertial operators. In [23], the overlapping costs of Fijany’s parallel approach was solved, adapting the RNEA algorithm by following the well-known Kogge and Stone recursive doubling technique [33].

Here, we adopted the solution proposed by [23], in order to solve the inverse dynamics model to compute the exoskeleton’s EoM into two computational recurrences: i) a forward propagation of spatial velocities (V) and accelerations ($\dot{V}$), ii) a backward propagation of spatial forces (F). Applying the Kogge and Stone method for each recurrence yields the structure presented in Figure 3, which reduces the equation set $\forall i$ of Equation (7) and (9) for a total of $lo{g}_{2}n+1$ computational steps (e) by powers of 2. The key element for this parallel approach is to identify the $i-1$ dependent term of each equation, which can be propagated to each computer node/process of the embedded system.

Note in Equation (9), the spatial velocity is composed by two terms: $V_i\leftarrow R_{i,i-1}^T{P}_{i,i+1}{V}_{i-1}+H{\dot{q}}_i$, in which $V_{i-1}$ can be propagated whereas $H{\dot{q}}_i$ can be computed in parallel for each computer node at the same computational step (e). The left plot of Figure 3 illustrates this process. While the serial $O\left(n\right)$ Algorithm 1 would require $e=8$ computational steps to solve $n=8$ equations for the spatial velocities $V_i,\forall i=1..n$, the parallel approach solves the same $n=8$ equations in only $e=3$ computational steps, since $O\left(lo{g}_2\left(8\right)\right)=3$. Likewise, spatial accelerations follow the same structure. In the case of the spatial forces, the propagation is backwards, as detailed by the right plot of Figure 3.

Therefore, scaling up to higher degrees of freedom (n), e.g. $n=512$ will required only $e=9$ computational steps, assuming the embedded systems incorporates 512 processing nodes/cores (p), i.e. $n=p$. If $n>p$, multi-threading can be applied for the parallel computation of the EoM. Algorithm 2 presents the $O\left(lo{g}_2\left(n\right)\right)$ parallel computation of the dynamics equations of motion for the soft exoskeleton.

In order to evaluate the performance of Algorithm 2, the Amdahl’s law is used [34]. Equation (10) describes the speedup metric, which is defined by the fraction of the algorithm that is inherently serial ($f_s$) and the portion that is parallelized ($f_p$):

Algorithm 2 $O\left(lo{g}_{2}n\right)$ RNE parallel computation

1. Input joint trajectory (therapy motion) for each joint i: $q_i,{\dot{q}}_i,{\ddot{q}}_i$   $\forall i$: do parallel   2. $B_{local}\leftarrow H_i{\dot{q}}_i$,   $M_{local}\leftarrow R_{i,i-1}^T{P}_{i,i+1}$   3. Spatial velocity: $V_i\leftarrow M_{local}{V}_{i-1}+B_{local}$   4. $B{A}_{local}\leftarrow H_i{\ddot{q}}_i+{\dot{H}}_i{\dot{q}}_i$,   $M{A}_{local}\leftarrow R_{i,i-1}^T{P}_{i,i+1}$,   $C{A}_{local}\leftarrow {\dot{R}}_{i,i-1}^T\dot{P_{i,i+1}}$   5. Spatial acceleration: ${\dot{V}}_i\leftarrow M{A}_{local}{\dot{V}}_{i-1}+C{A}_{local}{V}_i+B{A}_{local}$   6. Forward recurrence: spatial velocities and accelerations propagation:  for$e=1\to lo{g}_2(n-1)$do      if $i+2^{e-1}<n$          if $i-1<2^{e-1}$             send $V_i,{\dot{V}}_i$ to node $i+2^{e-1}$          else             send $M_{local},B_{local},M{A}_{local},B{A}_{local}$ to node $i+2^{e-1}$          endif      endif      if $i-1>2^{e-1}<n$          if $i-1<2^{e-1}$             receive $V_{i-1},{\dot{V}}_{i-1}$ from node $i-2^{e-1}$          else             receive $M_{local,i-1},B_{local,i-1},M{A}_{local},B{A}_{local}$ from node $i-2^{e-1}$          endif      endif      if $i>1$          if $e=lo{g}_2(i-1)+1$             $V_i=M_{local}{V}_{i-1}+B_{local}$             ${\dot{V}}_i=M{A}_{local}{\dot{V}}_{i-1}+C{A}_{local}{V}_i+B_{local}$          else             $M_{local}=M_{local}{M}_{local,i-1}$,   $B_{local}=M_{local}{B}_{local,i-1}+B_{local}$             $M{A}_{local}=M{A}_{local}M{A}_{local,i-1}$,   $C{A}_{local}=C{A}_{local}C{A}_{local,i-1}$,   $B{A}_{local}=B{A}_{local}B{A}_{local,i-1}+B{A}_{local}$          endif      endif   end for  7. $B{F}_{local}\leftarrow I_i{\dot{V}}_i+[I_i{\dot{S}}_{i,cm}+{\dot{I}}_i]V_i$,   $M{F}_{local}\leftarrow {P_{i,i+1}}^T{R}_{i,i-1}$   8. Spatial force: $F_i\leftarrow M{F}_{local}{F}_{i+1}+B{F}_{local}$   8. Backward recurrence: spatial forces propagation:  for$e=lo{g}_2\left(n\right)\to 1$do      if $i-2^{e-1}>0$          if $i-1>n-2^{e-1}$             send $F_i$ to node $i-2^{e-1}$          else             send $M{F}_{local},B{F}_{local}$ to node $i-2^{e-1}$          endif      endif    if $i-1<n-2^{e-1}$          if $i-1>n-2^{e-1}$             receive $F_{i+1}$ from node $i+2^{e-1}$          else             receive $M{F}_{local,i+1},B{F}_{local,i+1}$ from node $i+2^{e-1}$          endif      endif      if $i<n$          if $e=lo{g}_2(n-i)+1$             $F_i=M{F}_{local}{F}_{i+1}+B{F}_{local}$          else             $M{F}_{local}=M{F}_{local}M{F}_{local,i+1}$             $B{F}_{local}=M{F}_{local}B{F}_{local,i+1}+B{F}_{local}$          endif      endif   end for  Return $F_i$

$$\begin{array}{c}Speedup={\displaystyle \frac{1}{f_s+\frac{f_p}{n}}}={\displaystyle \frac{1}{(1-f_p)+\frac{f_p}{n}}}\hfill \end{array}$$

(10)

if $n=p$ the load distribution of processor (p) can be assigned to each individual joint of the exoskeleton (n: degrees of freedom). Otherwise, if $n>p$ communication packing and multi-threading processing are required.

3. Results

Algorithm 2 was coded in ANSI C programming language and implemented using the standard message-passing-interface (MPI) supported by the mpich library [35]. In this section, we analyze the performance of the $O\left(lo{g}_{2}n\right)$ RNE algorithm running in several computing architectures:

• Intel®Core™i7 Processor with 6 CPU cores at 2.6GHz and 512 GPU NVIDIA Quadro P1000 cores | 32GB RAM | Ubuntu 22.04.3 LTS.
• Hardware-in-the-loop (HIL) with the Embedded Nvidia Jetson™Nano with 1024 cores at 625MHz | 8GB RAM | Ubuntu 22.04.3 LTS.

Figure 4 details the computing time response of the proposed parallel approach running on the Intel®Core™i7 with NVIDIA Quadro GPU. The plot Figure 4a compares both serial $O\left(n\right)$ and parallel $O\left(lo{g}_{2}n\right)$ algorithms. As mentioned, both solutions for the inverse dynamics EoM were coded in ANSI C, using the high-performance Meschach framework; a comprehensive matrix-vector linear algebra for the operation of 6D spatial terms. The $O\left(lo{g}_{2}n\right)$ RNE approach significantly outperformed the serial algorithm $\forall n>16$ DoF, obtaining a real-time response of $t=0.127$ in the calculation of the inverse dynamics control for the soft exoskeleton with $n=512$ joints. In this test, we used an interpolated 3D spline trajectory with $pt=1,000$ knot-points. Contrary, the serial $O\left(n\right)$ RNE required $t=4.67s$ for the calculations. Considering the size of the input trajectory for the exoskeleton ($pt:$ knot-points), the net computational time is proportional to $pt*O\left(lo{g}_{2}n\right)$. In this sense, the inset within plot Figure 4b shows the resulting computing time of the parallel RNE approach with variations in the input trajectory knot-points.

The overall performance of the proposed solution is presented in Figure 5. The speedup exponentially increases after $n=p=64$, while achieving a $0.96$ degree of parallelism after approximately $n=p=256$, which demonstrates the required scalability for controlling highly-articulated and large-scale soft exoskeletons in real-time.

Another important element of the proposed $O\left(lo{g}_{2}n\right)$ RNE algorithm is related to the propagation structure of the EoM, in which the inherent sequential-time component of the algorithm can be drastically reduced if an embedded system-on-chip (SoC) is used. With an SoC-based embedded system, CPU/GPU cores are integrated with the memory in the same silicon, allowing higher data bandwidth to propagate data. This issue can be observed in Figure 6. Using the MPICH-library jumpshot feature, real execution with data propagation can be analyzed, by determining both serial and parallel computing time involved the calculation of the dynamics EoM. In this test, we used $n=p=8$ to clearly visualize the computational steps (e) followed by the $O\left(lo{g}_{2}n\right)$ propagation of spatial velocities ($V_i$), accelerations (${\dot{V}}_i$) and forces ($F_i$), as previously defined in Figure 3.

As previously mentioned, classical techniques used in soft-robotics that rely on finite element theory are not suitable in our application, since we require real-time computation running online while the exoskeleton. In this regard, we compared our $O\left(lo{g}_{2}n\right)$ RNE model against the Ansys®simulator, performing finite element analysis based on the Newton-Raphson iteration method applied to the proposed semi-soft structure driven by the compliant joint mechanism shown in Figure 2c. Also, as detailed in the inset of Figure 7a, the exoskeleton model is composed by one finger with 3 phalanges (the distal, middle, and proximal) performing a closing/opening trajectory with 90 knot-points with a step time of $0.0075s$. Under this topology Ansys used $42,220$ nodes for the semi-soft structure of the exoskeleton, requiring a simulation time of $1,255s$ and execution time of $4s$. Ansys determines the number of nodes and degrees of freedom (DoF) automatically. For the same trajectory with $1,000$ knot-points and a step time of $0.001s$, our $O\left(lo{g}_{2}n\right)$ RNE model required $0.12768s$ for 512 DoF, as shown in Figure 4b. Increasing up to $n=42,220$ DoF yielded a computing time of $6.8990s$, using a MPI configuration for $n>p$ (more degrees of freedom than processing cores). Even with this computing-intense configuration, our model dramatically outperformed Ansys results in terms of simulation time.

Also, Figure 7b presents computing time results for the soft model implemented in SoRoSim©Matlab™introduced by [29]. SoRoSim uses Geometric Variable Strain (GVS) approach for the dynamic simulation of soft systems. We tested several configurations from $n=10$ up to $n=512$ joints for the links assembly, as detailed in the insets of Figure 7b. For $n=512$, the method performed with a computational time of $0.918s$, whereas our method required $0.12768s$. Besides, the computational performance of SoRoSim exhibits a linear complexity $O\left(n\right)$, i.e., higher computing times will be obtained for large-scale articulated soft systems, compromising the scalability of the solution. In contrast, our model will scale with a complexity of $O\left(lo{g}_{2}n\right)$, as demonstrated in Figure 4.

Next, we present the performance of the $O\left(log2n\right)$ RNE algorithm running in a SoC-based embedded system with 1024 cores powered by the Nvidia Jetson Orin™Nano, with the ultimate goal of reducing the inherent serial time computation of Algorithm 2, and the advantage of compacting the proposed inverse-dynamics controller to support closed-loop robust control of the soft exoskeleton in real-time by following the proposed Hardware-in-the-loop (HIL) configuration. Figure 8 presents the experimental results obtained with this embedded system.

Both mpich and meschach libraries were compiled within the Nvidia Jetson Orin™Nano, allowing the computation of Algorithm 2 with similar computing times of those obtained from the Intel®Core™i7 Processor with 512 GPU NVIDIA Quadro cores, as presented in Figure 4b. For $n=p=512$, the parallel RNE algorithm required $0.1394s$ for $pt=1,000$ trajectory knot-points. As observed in Figure 8, the variance in computational time remains under $0.12s$, resulting feasible for in-loop real-time response.

Finally, Table 2 presents the computing time results for the aforementioned tests: $O\left(n\right)$ refers to Algorithm 1 running on the Intel®Core™i7 Processor with 512 GPU NVIDIA Quadro cores, SoroSim refers to the tests presented in Figure 7b running under Matlab™, $O\left(lo{g}_{2}n\right)$ also refers to Algorithm 2 running on the Intel®Core™i7 Processor, whereas HIL-$O\left(lo{g}_{2}n\right)$ refers to Algorithm 2 running on the embedded Nvidia Jetson Orin™Nano with 1024 cores. Furthermore, for the tests conducted in Ansys®, was not possible to configure the soft exoskeleton for different DoF values, since Ansys®determined $42,220$ DoF for the assembly shown in Figure 7a, requiring a simulation time of $1,255s$ ($20.91$ minutes).

4. Discussion

Algorithm 2 was implemented and tested on the embedded Nvidia Jetson Orin™Nano with 1024 cores, allowing the real-time calculation of the proposed inverse-dynamics controller based on the RNE method. As shown in Figure 8, the obtained computational complexity is of the form $pt*O\left(lo{g}_{2}n\right)$, being $pt$ the overall knot-points used to interpolate the reference trajectory for the exoskeleton.

In this regard, it is worth to highlight that those changes in the trajectory knot-points are associated with the precision and intensity of the therapy. In previous work reported in [28], we used a VICON™ motion capture system to extract the kinematics of motion associated to 6 hand gestures for rehabilitation. We implemented interpolation methods based on 3D splines to determine both cartesian and joint trajectories to drive the exoskeleton assistance. Here, the joint trajectory for positions ($q_i$), velocities (${\dot{q}}_i$) and accelerations (${\ddot{q}}_i$) are used as inputs to Algorithm 2, in order to estimate the computing time cost of the parallel inverse dynamics controller by analyzing the impact of increasing motion precision according to the computational complexity of the form $pt*O\left(lo{g}_{2}n\right)$.

Although similar computing times were obtained with the embedded platform in comparison to the desktop solution (c.f Table 2), the proposed HIL-driven embedded approach allow us to deploy a modular control system, facilitating the use of the soft exoskeleton in real rehabilitation scenarios. Also, we demonstrate that even with higher values of $pt$, the inherent logarithmic response $lo{g}_{2}n$ is maintained, supporting highly-articulated soft exoskeletons.

5. Conclusions

Classical finite element methods used in soft-robotics are not suitable for the real-time computation of model-dependent non-linear controllers, since most control algorithms require the dynamics calculation of the EoM within the closed-loop. As an alternative, our approach relies on the well-known Recursive Newton-Euler Algorithm (RNEA) calculated in a parallel fashion. As a result, a hardware-in-the-loop (HIL) structured allowed us to compute the proposed inverse-dynamics controller, by following a $O\left(lo{g}_{2}n\right)$ computational performance. Real-time response was achieved by testing our model on the Nvidia Jetson Orin™Nano embedded system, enabling to close the control-loop in $139ms$ ($n=512)$ for the hybrid soft-exoskeleton mechanism depicted in Figure 2c.

The EoM presented in Algorithm 1, allows us to incorporate external forces and perturbations for the robotic exoskeleton, that can be generated by the patient during real rehabilitation scenarios. Our non-linear control method based on the inverse-dynamics RNE is able to counteract these perturbations, as shown in step 3 from Algorithm 1, being the term $F_{i+1}$ the external payload to overcome with the control law ${\tau }_i$. Here, we have demonstrated that this control output can be generated in real-time thanks to the parallel computation of Algorithm 2.