A Human-like Inverse Kinematics Algorithm of an Upper Limb Rehabilitation Exoskeleton

Abstract

Powered exoskeleton rehabilitation is an effective way to help stroke patients recover their motor abilities. Bionic structures and human-like control strategies can be used to enhance both the safety and efficacy of exoskeletons. However, the motion characteristics of the shoulder complex are not sufficiently considered. In this paper, we designed a 7-degrees-of-freedom (DOF) upper limb rehabilitation exoskeleton, FREE (functional rehabilitation exoskeleton). The mechanical structures of the shoulder and forearm of FREE are in accordance with human anatomy, and can be used to perform a wide range of synergistic motion of multiple joints while keeping a safe distance from the patient’s head. A multiple-input-multiple-output (MIMO) shoulder girdle motion prediction model was developed to satisfy the synergy between humans and exoskeletons. Moreover, a constrained task priority and projected gradient-based inverse kinematics algorithm (CTPPG-IK) was proposed to achieve assistance with scapulohumeral rhythm. A motion capture system was used to collect different activities of daily life (ADL) motion data to validate the proposed algorithm. The experimental results show that the accuracy of the prediction model is higher than that of existing models, and the inverse kinematics algorithm can handle the end-effector task and joint space with a maximum angle error of $3.04\times {10}^{-3}$ rad.

1. Introduction

Stroke has been a widespread and serious global healthcare problem, with the highest rate of disability for a single disease. There were 13.7-million new cases of stroke worldwide in 2016, and this number has continued to rise [1]. The most common and widely recognized impairment following stroke is motor impairment, and rehabilitation is an effective way to help patients regain their motor abilities [2]. To date, several rehabilitation robots have been used to improve both motor control and strength in post-stroke upper extremity paralysis, which can help speed up neural remodeling significantly, and appear to be safe [3,4,5,6].

The upper limb is one of the many redundant subsystems of the human body and is commonly considered to have 7 DOF, consisting of three at the shoulder, one at the elbow, and three at the wrist. Positioning the wrist in space and orienting the palm is a task that requires only 6 DOF [7]. Mechanical engineers have drawn inspiration from human upper limb to design spherical–roll–spherical (SRS) robots such as KUKA IIWA and Franka Panda. However, the acromioclavicular joint in the shoulder actually provides additional degrees of freedom, a feature that is often overlooked in the design of exoskeletons. To ensure a safe and comfortable human–robot interface, precise alignment between the exoskeletal and anatomical joint axes throughout movement phases is critical. This prevents the generation of harmful tangential forces on soft tissues and joints, minimizing the risk of discomfort or injury. There are some exoskeletons that use passive self-alignment to achieve center alignment of the shoulders. ASSISTON-SE is equipped with a 3RRP parallel back-drivable mechanism, allowing for passive translational movements of the center of the GH joint and enabling independent active control of these DOFs [8]. NESM-$\gamma$ has an adjustable trunk link with passive ball joints at the two ends and connects the back of the user and the robot, allowing for unconstrained movements of the scapula [9]. Based on the concept of posture synergies, a semi-actuated upper limb exoskeleton Armule was designed, which achieves 5 DOF motion of the upper limb with two motors [10,11]. The above-mentioned exoskeletons are adjustable according to the user characteristics, provide repetitive active and passive assistance, accelerate neural remodeling, and record user motion data, but there are also limitations. For example, exoskeletons do not fully match to the upper limb kinetic chain, resulting in insufficient active degrees of freedom, safety, sensing and learning capabilities, and comfort.Therefore, the design objectives of the exoskeleton mainly include the following four points: (1) redundant active degrees of freedom, (2) bionic structure, (3) non-collision with the user, and (4) torque sensing and control [12,13].

Nevertheless, the passive alignment method may have problems, such as lack of adaptability, low comfort, and a complicated process for putting it on and taking it off, so only a few exoskeletons have active shoulder compensation mechanisms currently [14,15]. This limitation affects the patient’s ability to perform functional movements in rehabilitation training, potentially causing discomfort and compensatory movements. Furthermore, realizing the intuitive and ergonomic interaction between the user and exoskeleton can significantly enhance the rehabilitation outcomes. One key requirement is that the exoskeleton generates natural and anthropomorphic movements [16]. Human upper limb movements not only focus on the hand’s motion but also the coordination between joints. In the realm of robot kinematics, they correspond to task space and joint space movements, respectively [17].

Rehabilitation exoskeletons are highly coupled to human limbs; thus, the corresponding kinematics approach deserves more consideration. There are several approaches to the inverse kinematics (IK) of exoskeletons, including analytic geometry methods, learning-based methods, and differential kinematics. Analytic methods based on the task space are efficient and fast [18,19,20,21]. These methods require the parameterization of the upper limb in advance. Learning-based methods are usually computed in the joint space, enable exoskeletons to learn and reproduce the patient’s movements in an unstructured environment, and imitate human movement habits [22,23,24]. This method is applicable to nonstrict motion trajectories and require significant computational resources during the training phase. Differential kinematics can be combined with the gradient projection method to complete end-effector tasks and optimize the joint motions of redundant robots [25,26,27]. A kinematically redundant robot is a mechanism that has more DOF than are required to perform a given task. The upper limbs are redundant with respect to the 6 DOF of the hand, so upper limb exoskeletons should also be kinematically redundant for flexible manipulation. For these robots, the problem generally cannot be solved analytically and has an infinite number of solutions, so numerical methods are usually used to solve local optimization problems under joint angle constraints [28,29]. Most of the existing research studies assisted patients with tasks of the hand but did not take into account the movement patterns of the upper limb. Scapulohumeral rhythm is a kinematic interaction between the scapula and humerus, and results in non-coincidence between the center of the exoskeletal shoulder and the humeral origin. In the absence of kinematic restrictions, exoskeletons may compress or pull patients and cause pain. Therefore, the IK which unifies end-effector tasks and synergistic movement among links is worth researching. However, few studies have integrated scapulohumeral rhythm and kinematic approaches. The relative inclination angle and geodesics in Riemannian space were used to plan human-like paths [30].

Overall, the main contribution of this paper is as follows:

• A multiple-input multiple-output (MIMO) shoulder girdle motion prediction model is proposed. It can be used to estimate elevation/depression and protraction/retraction angles of the shoulder girdle based on the humerus orientation.
• An IK algorithm that integrates the gradient projection method and task priority is proposed and validated. It enables FREE to achieve human-like shoulder movement.

2. System Design

The FREE exoskeleton was developed to assist patients with functional rehabilitation training after stroke. Seven motors shown in Figure 1 were configured to achieve shoulder joint, elbow joint and forearm motion after compromises. The exoskeleton is capable of shoulder girdle elevation/depression and forward extension/retraction; shoulder joint abduction/adduction and flexion/extension; glenohumeral (GH) joint self-rotation; elbow flexion/extension; and forearm pronation/supination. To increase the applicability of the exoskeleton, the upper arm and forearm were designed as a split structure for length adjustment. To improve the stiffness of the exoskeleton and reduce its mass, the body was mainly made of aluminum material and carbon fiber tubes.

J1–J5 were used to construct the exoskeleton shoulder joint. The GH joint is a visual characterization of shoulder motion. Almost all voluntary movements of the shoulder joint involve movement at the GH joint; this requires additional degrees of freedom to ensure that the GH joint is in central alignment with the exoskeleton shoulder during rehabilitation training [31]. To achieve synchronization between the exoskeleton and shoulder without overlap, a parallelogram was used for translation of the axis of rotation of the scapulohumeral out of the patient’s body.

The 3D model of FREE is shown in Figure 2a. In order to achieve a wide range of ROM, FREE adopts a serial chain with three revolute joints that encircle the upper arm for the ball-and-socket joint to match GH. The angle between the two adjacent axes is set at $4\pi /9$ rad to ensure sufficient manipulability of shoulder. In order to achieve a greater level of horizontal extension, the upper arm linkage exists a $\pi /6$ rad bend. This configuration places the singularities at the edge of the ROM while leaving a safety distance of more than 40 mm between the human and exoskeleton. A new linkage structure was developed to achieve pronation/supination of the forearm shown in Figure 2b. This structure consists of only five connecting links, which can translate the forearm rotation axis into the isometric plane beyond exoskeleton links. Compared to ARMin [32], this structure weighs only 210 g with the same range of motion (approximately $8\pi /9$ rad), dramatically reducing weight while having lower structural complexity.

The controller was designed based on the industrial PC CX6015-0100 (Beckhoff Automation Inc., Verl, Germany) under the software environment of TwinCAT3, and the hardware architecture is shown in Figure 3. Both the control algorithm and the communication cycle were strictly kept to 1 ms.

3. Kinematics Algorithm

Human upper limbs are able to perform various movements dexterously thanks to its kinematic redundancy. In order to match the kinematics of the patient’s upper extremity, the designed exoskeleton contains seven joints with redundant degrees of freedom. After planning the path of a set of rehabilitation motions, the joint angles are calculated by IK. Although the analytical geometry method is fast and always finds a usable solution, it is applicable only to robots of a specific configuration and is incompatible with the constraints of the joint space. For FREE, the ball-and-socket joint was placed at the midpoint of the kinematic chain; it is difficult to decouple the position from orientation. In addition, it is difficult to obtain the analytical solution of a redundant exoskeleton in the absence of an arm angle constraint, in advance. In this section, a model is constructed to predict the motion of the shoulder girdle for the constraint of scapulohumeral rhythm. A task priority and projected gradient-based IK solution method is proposed and prioritized for multiple tasks to improve the human–robot coordination during rehabilitation training.

3.1. Human Shoulder Kinematics

In order to coordinate movements between exoskeletons and patients, it is of the utmost importance to clarify the kinematic properties of the human upper limb, especially the shoulder complex. The shoulder complex consists of four sub-joints. Though each joint has independent degrees of freedom, the shoulder complex can be considered a compound mechanism, where all four joints must properly interact for normal shoulder motion to occur. A case in point for this interaction is the scapulohumeral rhythm, which is the relationship between the elevation of the scapulohumeral joint and the upward rotation of the scapula. Reference [30] used two polynomials to describe the relationship between the polar angle of the humerus and the elevation/depression and protraction/retraction of the scapulohumeral joint; however, this description is not comprehensive. In addition, movements of the shoulder girdle are associated with ${\phi }_{\mathrm{az}}$ in Figure 4. As the humerus of the upper arm moves, the scapulohumeral follows fixed patterns. Reference [33] shows that using several independent regression models may not be able to accurately describe the movement of the scapula and humerus bones. They maintain a non-linear relationship in multiple dimensions of motion.

The coupled movements of the shoulder joint result in a nonlinear relationship between the movement angles of the humerus and scapulohumeral. One common approach to dealing with multi-input multi-output regression problems is to build multiple multi-input single-output models and predict each variable separately. We trained and merged two BP neural networks but the prediction performance for ${\phi }_{\mathrm{pr}}$ was suboptimal. In order to improve prediction accuracy, the relationship can be abstracted into a MIMO model that takes into account the potential correlations between different outputs; however, it has not currently been applied in the human joint angles prediction domain. Xu presented a MLS-SVR multi-output algorithm, in a multi-output setting with a strong generalization capability [34]. Considering the limited amount of data and the requirement for a smooth fit, MLS-SVR was chosen as a predictor of the rotation angles of the shoulder girdle. A motion capture system can be used to record the movements of the upper limbs and calculate the four predefined sets of angles. Let $y=[{\phi }_{\mathrm{ed}},{\phi }_{\mathrm{pr}}]\in R^{n\times m}$, $x=[{\phi }_{\mathrm{ed}},{\phi }_{\mathrm{pr}}]\in R^{n\times p}$, where n is the sample size, and m and p represents the number of dimensions for the output and input variables, respectively. The multi-output regression is regarded as finding the mapping between an input matrix $[{\phi }_{\mathrm{po}},{\phi }_{\mathrm{az}}]$ and an output matrix $[{\phi }_{\mathrm{ed}},{\phi }_{\mathrm{pr}}]$.

In order to achieve accurate prediction of shoulder joint angles, the optimal parameters of prediction model need to be solved. The goal of the MLS-SVR algorithm is to find a function $f\left(x\right)=\phi {\left(x\right)}^{\mathrm{T}}W+b$ that has the maximum deviation from the actual observed output and is as flat as possible, where $W\in R^{h\times l}$, $b\in R^l$, l is the dimension of the output variable, and $\phi \left(x\right)$ is a a mapping function to achieve the linear regression of data in high-dimensional space. The element $w_i\in R^h$ of W can be written as $w_i=w_0+v_i$. $w_0$ represents the mean vector, while small values of the vector $v_i\in R^h$ indicate that the different outputs are similar to one another. The flatness of the function $f\left(x\right)$ can be obtained by minimizing the value of the parameters $w_0^{\mathrm{T}}{w}_0$ and $trace$$\left(V^{\mathrm{T}}V\right)$, which is equivalent to finding a small value for W, where $V=\left(v_1,v_2,\dots ,v_l\right)\in R^{h\times l}$. Therefore, the optimization problem could be attained by minimizing the following objective function while satisfying the constraint:

$$\begin{array}{cc}\hfill min\mathcal{J}\left(w_0,V,\Xi \right)=& \frac{1}{2}{w}_0^{\mathrm{T}}{w}_0+\frac{\gamma ″}{2l}trace\left(V^{T}V\right)+\frac{\gamma \prime }{2}trace\left({\Xi }^{\mathrm{T}}\Xi \right)\hfill \end{array}$$

(1)

$$\mathrm{s}.\mathrm{t}.y=Z^{\mathrm{T}}W+repmat\left(b^{\mathrm{T}},l,1\right)+\Xi$$

(2)

where $\Xi =\left({\xi }_1,{\xi }_2,\dots ,{\xi }_l\right)\in R^{n\times l}$ consists of slack variables, $b=\left(b_1,b_2,\dots ,b_l\right)\in R^l$, $\gamma \prime$ and $\gamma ″$ are hyper-parameters, $Z=\left(\phi \left(x_1\right),\phi \left(x_2\right),\dots ,\phi \left(x_n\right)\right)\in R^{h\times n}$. The optimization in Equations (1) and (2) can be made into a Lagrange function as follows:

$$\begin{array}{c}\hfill \mathcal{L}\left(w_0,V,b,\Xi ,A\right)=\mathcal{J}\left(w_0,V,\Xi \right)-trace\left(A^T(Z^{\mathrm{T}}W+repmat\left(b^{\mathrm{T}},n,1\right)+\Xi -y)\right)\end{array}$$

(3)

where $A=\left({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_m\right)\in R^{l\times m}$ is a matrix consisting of Lagrange multipliers. After eliminating W and $\Xi$ according to the KKT condition, the following linear equation system can be determined:

$$\left[\begin{array}{cc}{0}_{l\times l}& N^{\mathrm{T}}\\ N& M\end{array}\right]\left[\begin{array}{c}b\hfill \\ \alpha \hfill \end{array}\right]=\left[\begin{array}{c}{0}_l\hfill \\ y\hfill \end{array}\right]$$

(4)

where $N=blockdiag\left(1_n,1_n,\dots ,1_n\right)\in R^{ln\times l}$, $\Omega =repmat(\tilde{K},l,l)\in R^{ln\times ln}$, $M=\Omega +{\left(\gamma \prime \right)}^{-1}{I}_{ln}+\left(l/\gamma ″\right)Q\in R^{ln\times ln}$ is a positive definite matrix, $\tilde{K}=Z^{\mathbf{T}}Z\in R^{n\times n}$ is a matrix that contains the kernel function, and $Q=blockdiag(\tilde{K},\tilde{K},\dots ,\tilde{K})\in R^{ln\times ln}$, $\alpha ={\left({\alpha }_1^{\mathrm{T}},{\alpha }_2^{\mathrm{T}},\dots ,{\alpha }_m^{\mathrm{T}}\right)}^{\mathrm{T}}\in R^{ml}$. However, solving the linear equation system (4) directly is difficult because it is not positively definite. Therefore, the equation was reformulated into the following equation:

$$\left[\begin{array}{cc}G& 0_{l\times ln}\\ 0_{ln\times l}& M\end{array}\right]\left[\begin{array}{c}b\\ M^{-1}Nb+\alpha \end{array}\right]=\left[\begin{array}{c}{N}^{\mathrm{T}}{M}^{-1}y\\ y\end{array}\right]$$

(5)

where $G=N^{\top }{M}^{-1}N\in R^{l\times l}$ is a positive definite matrix. Then, the solution $\tilde{b}$ and $\tilde{\alpha }$ can be calculated by $\tilde{b}=G^{-1}{N}^{\mathbf{T}}{H}^{-1}y$ and $\tilde{\alpha }=M^{-1}(y-Nb)$, and the MLS-SVR decision function can be obtained as follows:

$$\begin{array}{cc}\hfill \tilde{f}\left(x\right)& =\phi {\left(x\right)}^{\mathrm{T}}\tilde{W}+{\tilde{b}}^{\mathrm{T}}\hfill \\ \hfill & =\phi {\left(x\right)}^{\mathrm{T}}repmat(\widetilde{w_0},1,l)+\phi {\left(x\right)}^{\mathrm{T}}\tilde{V}+{\tilde{b}}^{\mathrm{T}}\hfill \\ \hfill & =\phi {\left(x\right)}^{\mathrm{T}}repmat\left(\sum _{k=1}^l{Z}^{\mathrm{T}}{\tilde{\alpha }}_k,1,l\right)+\frac{l}{\gamma ″}\phi {\left(x\right)}^{\mathrm{T}}Z\tilde{A}+{\tilde{b}}^{\mathrm{T}}\hfill \\ \hfill & =repmat\left(\sum _{j=1}^l\sum _{i=1}^n{\tilde{\alpha }}_{ij}K\left(x,x_i\right),1,l\right)+\frac{l}{\gamma ″}\sum _{i=1}^n{\tilde{\alpha }}_{i}K\left(x,{\mathbf{x}}_i\right)+{\tilde{b}}^{\mathrm{T}}\hfill \end{array}$$

(6)

In this study, the RBF (radial basis function) kernel was used because it has only one hyper-parameter and is computationally efficient: $K\left(x,x_i\right)=exp{\left(-\sigma ∥x-x_i∥\right)}^2,\sigma >0$. This approach can be used to predict ${\phi }_{\mathrm{pr}}$ and ${\phi }_{\mathrm{ed}}$. They can be set as the target angles of ${\theta }_1$ and ${\theta }_2$, and ${\theta }_1$ and ${\theta }_2$ are actual joint values of exoskeletons. In addition, in order to reduce jitter in the prediction results and to prevent drastic changes in the angular velocities of the joints, a zero-phase low-pass filter was added after the prediction model. The cutoff frequency was set to 15 Hz.

3.2. Robot Kinematics

The standard Denavit–Hartenberg (SDH) convention is a widely used approach for robot modeling. However, it suffers from a significant drawback when dealing with robots that have tree-like or closed-chain structures, where some links may possess more than one axis of rotation. In such cases, employing the SDH method to establish link coordinate systems can lead to ambiguities. FREE features a closed-chain structure in its shoulder girdle. Consequently, the local coordinate system of each link of FREE was established according to the modified Denavit–Hartenberg (MDH) convention. Figure 5 shows the parallelogram structure of exoskeleton shoulder. Because of the presence of the parallelogram linkage at the shoulder girdle, the exoskeleton is modeled differently from the joint kinematics of the open-chain robot. Assuming a local coordinate system attached to the fixed end of the parallelogram, the linkage relative to the fixed end only changes its position but not its orientation during the motion of the mechanism.

The body twist can be defined as the spatial velocity in the body frame, expressed as

$${\mathcal{V}}_b=\left[\begin{array}{c}{\omega }_b\\ v_b\end{array}\right]=\left[\begin{array}{c}{}^o\omega \\ {}^o\omega \times {}^o\rho \end{array}\right]\dot{\theta }$$

(7)

where ${\mathcal{V}}_b$ is the spatial velocity of body B with respect to a fixed frame. ${\omega }_b$ and $v_b$ are the angular and linear velocity vectors of the instant point on the rigid body, respectively. ${}^o\omega$ and ${}^o\rho$ are the direction and location vectors of the rotational axis with respect to the fixed frame, respectively. $\dot{\theta }$ is the angular velocity of the link with respect to the rotational axis. When the parallelogram rotates, the link adjacent to the fixed end rotates at an angular velocity $\dot{\theta }$, and the link opposite to the fixed end rotates at an angular velocity $-\dot{\theta }$. Therefore, the kinematic rotation of the parallelogram mechanism can be rewritten as

$$\begin{array}{cc}\hfill {\mathcal{V}}_b& =\left[\begin{array}{c}{}^o{\omega }_{ad}\\ {}^o{\rho }_{ad}\times {}^o{\omega }_{ad}\end{array}\right]\dot{\theta }+\left[\begin{array}{c}{}^o{\omega }_{op}\\ {}^o{\rho }_{op}\times {}^o{\omega }_{op}\end{array}\right](-\dot{\theta })\hfill \\ \hfill & =\left[\begin{array}{c}0\\ ({}^o{\rho }_{op}-{}^o{\rho }_{ad})\times {}^o\omega \end{array}\right]\dot{\theta }\hfill \\ \hfill & ={}^o{\widehat{s}}_p\dot{\theta }\hfill \end{array}$$

(8)

where ${}^o{\omega }_{ad}$, ${}^o{\rho }_{ad}$, ${}^o{\omega }_{op}$, and ${}^o{\rho }_{op}$ are the direction and location vectors of the joint axis of each link adjacent to and opposite to the fixed end. ${}^o{\widehat{s}}_p$ is the spatial vector on the output side of the parallelogram. ${}^o\omega$ is the component of spatial velocity. Determining the end effector motion from a given joint motion is the subject of the forward kinematics problem. The MDH parameters of FREE are shown in

Table 1. MDH parameters table.

i	${\mathit{a}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$	${\mathit{\alpha }}_{\mathit{i}}$	${\mathit{\theta }}_{\mathit{i}}$
1	0	0	$\pi$	${\theta }_1$
$2\prime$	0	0	$-\pi /2$	${\theta }_2$
2	rotates with $-{\theta }_2$
3	0	0	$-5\pi /6$	${\theta }_3$
4	0	0	$4\pi /9$	${\theta }_4$
5	0	−0.563	$-4\pi /9$	${\theta }_5$
6	0	0.508	$\pi /6$	${\theta }_6$
7	0	0.333	$\pi /2$	${\theta }_7$

.

By using IK, the angle of each robot joint can be solved, given the position and direction of the endpoint. However, for redundant robots, the Jacobi matrix is not a square, and thus the inverse matrix cannot be directly solved. One solution is to use pseudo-inverse substitution. At the velocity level, they are expressed in a vector called the Jacobian IK (J-IK) form as

$$\dot{x}=J\dot{\theta }$$

(9)

$$\dot{\theta }=J^{†}\dot{x}$$

(10)

where $\dot{x}$ is the m-dimensional linear and rotational velocities vector of the end effector, and $\dot{\theta }$ is the n-dimensional joint angles vector. J is the $m\times n$ robot Jacobian matrix, and the $n\times n$ matrix $J^{†}$ is the pseudo-inverse of J, also called the Moore–Penrose inverse of J. It can be computed as $J^{†}=J^T{\left(J{J}^T\right)}^{-1}$. The pseudo-inverse gives the best possible solution to Equation (7) in the sense of least squares; however, the pseudo-inverse has stability problems in the region adjacent to the singularity. If the configuration is close to the singularity, then the pseudo-inverse method will lead to drastic changes in the joint angles, exceeding even the limit of joint velocity, even if the desired speed is extremely small.

The damped least squares method (DLS) avoids many of the singularity problems of the pseudo-inverse method and gives a numerically stable solution, which was first used for IK by Nakamura [35]. Thus, $J^{†}$ was rewritten as $J_D^{†}$:

$$J_D^{†}=J^{\mathrm{T}}{\left(J{J}^T+{\lambda }^{2}I\right)}^{-1}$$

(11)

The pseudo-inverse solution of Equation (8) minimizes the two-norm number of the joint velocity vector; however, the posture of the exoskeleton in this case may not match the current posture of the patient’s upper limb. For rehabilitation exoskeletons, it is clearly important to maintain kinematic consistency with the patient’s upper extremities. The redundant exoskeleton design can achieve this requirement. The projected gradient method was first proposed to solve IK with the pseudo-inverse [36]. The discrete form is used in place of the differential form for computational convenience as

$$\Delta \theta =J_D^{†}\Delta x+k(I-J_D^{†}J){\xi }_1$$

(12)

where ${\xi }_1$ is an arbitrary vector, and I is a $n\times n$ identity matrix. $(I-J^{†}J)$ is an operator which projects ${\xi }_1$ onto the null space of the Jacobian J such that ${\xi }_1$ can be interpreted as a desired velocity behavior that is only effective in the null space and does not interfere with task achievement. While FREE has 7 degrees of freedom, it has insufficient redundancy to perform the end-effector task after constraining both joint angles of the shoulder girdle mechanism. This is equivalent to controlling the remaining five joints to complete 6 DOF end-effector tasks. Under constraints of certain subtasks, if all the variables orientation and position are used to calculate the joint angles, it may cause a mismatch between the task and redundant degrees of freedom, and thus the solution will fail. Orientation-related tasks and position-related tasks are decoupled at the human wrist. In ADL motions, position-related tasks affect the distance error of the motions, while orientation-related tasks are related to human habits and comfort. For example, a user always grasps a bottle in the habitual posture, but this can be performed by grasping either the top or the bottom of the bottle to complete the task. Exoskeletons should perform shoulder assist and forearm posture tasks while maximizing the position of the end effector. Orientation errors between the end effector of the exoskeleton and patient’s wrist may lead to discomfort and injury. Therefore, exoskeletons should perform suitable orientation-related tasks while maximizing shoulder assistance and minimizing the position error of the end effector. For this scenario, we use priority-based task control with task orientation as the main priority and the task position as the secondary priority, while introducing constraints on the human shoulder. Therefore, the end-effector task can be divided into orientation-related and position-related tasks [37], with orientation-related tasks having a higher priority in this study:

$$\Delta \theta =J_{\omega }^{†}\Delta x_{\omega }+k(I-J_{\omega }^{†}{J}_{\omega })\left(\nabla H_{\mathrm{shr}}\left(\theta \right)+{\xi }_v\right)$$

(13)

$$H_{\mathrm{shr}}\left(\theta \right)=\frac{1}{2}\sum _{i=1}^2{\left({\theta }_i-{\widehat{\theta }}_i\right)}^2$$

(14)

where the first term on the right side corresponds to the high-priority orientation task, and the second term corresponds to the sub-priority task containing the end position and shoulder constraints. $J_{\omega }^{†}$ and $J_v^{†}$ are the first and last three rows of $J_D^{†}$, which are associated with the angular and linear velocities, respectively. $x_{\omega }$ and $x_v$ are the first and last three rows of x, which are associated with the orientation and position, respectively. ${\xi }_v=J_v^{†}\Delta x_v$ represents the spatial vector of the linear velocities subtask. $H_{\mathrm{shr}}\left(\theta \right)$ represents the joint space loss function, which is used to regulate the synchronized motion of the exoskeleton with the user’s shoulder; shr is short for shoulder. $\nabla H_{\mathrm{shr}}\left(\theta \right)=k(\theta -\widehat{\theta })$ represents the gradient of $H_{\mathrm{shr}}\left(\theta \right)$; when $H_{\mathrm{shr}}\left(\theta \right)\ne 0$, this term adjusts the angle of the shoulder movement in the direction of the gradient. ${\widehat{\theta }}_i$ represents the desired shoulder girdle movement angle, i.e., the previously mentioned ${\phi }_{\mathrm{ed}}$ and ${\phi }_{\mathrm{pr}}$. The weight factor is fixed to $k=1$ in this study. The pseudocode is given below, where $\mathrm{ori}$ is short for orientation and $tol$ is a small positive tolerance of task error. The complete pseudocode is given in Algorithm 1:

Algorithm 1 Constrained task-priority and projected-gradient IK algorithm (CTPPG-IK).

1: while  $e_{\mathrm{ori}}>{tol}_{\mathrm{ori}}\&\&e_{\mathrm{joint}}>{tol}_{\mathrm{joint}1}\&\&e_{\mathrm{joint}}>{tol}_{\mathrm{joint}2}$do 2:     $\mathrm{compute}\Delta x_i,J_D^{†}\left({\theta }_i\right),{\phi }_{\mathrm{ed}},{\phi }_{\mathrm{pr}},{\Delta x}_{\omega },{\Delta x}_v,J_{\omega }^{†},J_v^{†}$ 3:     $\widehat{\theta }={[{\phi }_{\mathrm{ed}},{\phi }_{\mathrm{pr}}]}^{\mathrm{T}}$ 4:     ${\xi }_v=J_v^{†}\Delta x_v$ 5:     $\nabla H_{\mathrm{shr}}\left(\theta \right)=[k{({\theta }_i-\widehat{\theta }),0,0,0,0,0]}^{\mathrm{T}}$ 6:     $\Delta {\theta }_i=J_{\omega }^{†}\Delta x_{\omega }+k(I-J_{\omega }^{†}{J}_{\omega })\left(\nabla H_{\mathrm{shr}}\left(\theta \right)+{\xi }_v\right)$ 7:     ${\theta }_{i+1}={\theta }_i+\Delta {\theta }_i$ 8:     $\mathrm{compute}{x}_{i+1}$ 9:     $e_{\mathrm{ori}}=\left|\right|x_{\omega (i+1)}-x_{\omega \left(i\right)}\left|\right|$ 10:     ${[e_{\mathrm{joint}1},e_{\mathrm{joint}2}]}^{\mathrm{T}}=[I_{2\times 2},0_{2\times 5}]{\theta }_{i+1}-\widehat{\theta }$ 11:     $i=i+1$ 12: end while

3.3. Convert Generated Motion to Exoskeleton Motion

In an ideal scenario, it is expected for the motion of the user’s upper limb and the exoskeleton to be perfectly synchronized. However, due to the presence of redundant DOF, using the markers set on the user’s wrist to obtain the Cartesian trajectory of the end effector and solving for the joint angles vector $\theta$ of the exoskeleton cannot guarantee motion synchrony. The goal of this study is to accurately predict the shoulder joint angle and enable human-like auxiliary movements for the exoskeleton; motions in the upper limb configuration space can be mapped to the exoskeleton configuration space by the analytical geometry method. The shoulder complex, upper arm and forearm, which make up the human upper limb, can be simplified as three vectors in series. The simplified systems are shown in Figure 6.

The origins of the local coordinate of link 1 and link 2 are coincident and denoted as $O_{12}$. Similarly, the origin of the local coordinate of link 3, link 4, and link 5 are coincident and denoted as $O_{345}$. When $z_1$ is aligned with the center of the ball on the shoulder girdle, the user’s upper arm and forearm lengths are parallel and congruent to those of the FREE. The vector set of FREE is equivalent to those of the upper limb in kinematics, and $O_{\mathrm{c}}$, $O_{\mathrm{s}}$, $O_{\mathrm{e}}$, and $O_{\mathrm{w}}$ coincide with $O_{12}$, $O_{345}$, $O_6$, and $O_7$, respectively. The coordinates of $O_{\mathrm{c}}$ can be calculated by the least squares method using the data of the GH center. The orientation and position transform can be obtained by a homogeneous rotation matrix of MDH as follows:

$$T_i^{i-1}=\left[\begin{array}{cc}{}_0^{n}R& {}_0^{n}P\\ 0& 1\end{array}\right]=\left[\begin{array}{cccc}c{\theta }_i& -s{\theta }_i& 0& a_i\\ c{\alpha }_{i}s{\theta }_i& c{\alpha }_{i}c{\theta }_i& -s{\alpha }_i& -d_{i}s{\alpha }_i\\ s{\alpha }_{i}s{\theta }_i& s{\alpha }_{i}c{\theta }_i& c{\alpha }_i& d_{i}c{\alpha }_i\\ 0& 0& 0& 1\end{array}\right]$$

(15)

The coordinates of points $O_{\mathrm{c}}$, $O_{\mathrm{s}}$, $O_{\mathrm{e}}$, and $O_{\mathrm{w}}$ in the world coordinate system can be directly obtained by the motion capture system. The vector $v_{\mathrm{cs}}$, $v_{\mathrm{se}}$, and $v_{\mathrm{ew}}$ can be calculated by taking the difference between adjacent reference points. According to the principle of motion synchrony, it is easy to observe that the shoulder girdle angles of FREE ${\theta }_{1b}$ and ${\theta }_{2b}$ correspond to those of the user ${\phi }_{\mathrm{ed}}$ and ${\phi }_{\mathrm{pr}}$, respectively. ${\theta }_{ib}$ represents the baseline angle of each joint of the exoskeleton. ${\theta }_{1b}$ and ${\theta }_{2b}$ can be calculated by bringing them into forward kinematics calculations as shown in the following formula:

$$v_{\mathrm{cs}}=T_2^0(1:3,4)$$

(16)

Similarly, ${\theta }_{3b}$, ${\theta }_{4b}$ and ${\theta }_{5b}$ can be calculated by $v_{\mathrm{se}}$ as follows:

$$v_{\mathrm{se}}=T_5^2(1:3,4)$$

(17)

${\theta }_{6b}$ is the angle between $v_{\mathrm{se}}$ and $v_{\mathrm{ew}}$, and ${\theta }_{7b}$ is the angle between $v_{\mathrm{se}}$ and $v_{\mathrm{ew}}$ as follows:

$${\theta }_{6b}=\pi /2-\frac{〈x_s,x_e〉}{∥x_s∥∥x_e∥},{\theta }_{7b}=\frac{〈z_e,z_w〉}{∥z_e∥∥z_w∥}$$

(18)

This calculation method has a clear geometric meaning and ensures both the consistency of the end trajectories and the synchronization of the kinematic chains. The calculated ${\theta }_{ib}$ is the baseline joint angles vector, which can be used to calculate the error of the shoulder angles prediction model. By calculating the forward kinematics through ${\theta }_{ib}$, the obtained trajectory can be considered the baseline trajectory. This trajectory can be used to calculate the end-effect error.

4. Experiments and Results

To validate the proposed human-like kinematics algorithm, our experimental objectives are twofold: (1) assess the precision of the shoulder girdle motion prediction model, and (2) examine the ability of FREE to achieve rehabilitation actions while simultaneously satisfying both shoulder and end-effector constraints by utilizing the joint angles calculated by the prediction model.

4.1. Experimental Protocol

Users of different body sizes have different shoulder rhythms, which lead to variability in the prediction model parameters. In order to make the prediction model as accurate as possible in predicting the shoulder girdle angle for a specific user, one right-handed healthy subject (male, age: 24 years, height: 177 cm, weight: 60.5 kg) volunteered to participate in the study. The testing was performed with the informed consent of the volunteer and in a protocol approved by the Declaration of Helsinki. This experiment is a benign behavioral intervention for adults. The experiment was completed in a quiet indoor environment and was completely free of charge. Volunteers were fully informed of the procedure and content of the experiment, and they signed an informed agreement. To ensure the rigor of the experiment and the accuracy of the results, the definition of a joint coordinate system for the shoulder, elbow, and wrist, proposed by the Standardization and Terminology Committee (STC) of the International Society of Biomechanics (ISB), was used [38]. An optical motion capture system with eight cameras (Mars2H, Nokov, sampling rate 60 Hz) was used to record the data of motions and reflectable marker points as the reference sensing method. The volunteers performed four ADL motions in Figure 7: reaching to the mouth (drinking/eating), cleaning the table, touching the head, and arm cycling.

The first three motions have a small ROM and a fixed pattern, and the last motion includes a larger ROM amplitude. Data containing information from a wide ROM helped identify relatively accurate humeral–clavicular kinematic patterns at the baseline. In order to avoid occlusion of the markers and ensure no data loss, movements in the experiment were performed at a slower pace compared to those in daily life. The volunteer was asked to complete each movement ten times in sequence within a two-minute time frame to ensure the continuity of data. During the whole process, the volunteers tried to keep the torso as immobile as possible. The experiment was conducted twice, with one dataset used for training the prediction model and the other dataset used to generate reference movements.

4.2. Data Processing and Analysis Criteria

To compute the values of $\alpha$ and b to predict cleaning, drinking, and arm circulation, the data of these three motions were used to train the predictive models as shown in Equation (3). After training was completed, to evaluate the performance of the shoulder girdle angles prediction model, R-square and RMSE indices were used. In order to validate the performance of the algorithms, the analytic geometry method was used to solve the predefined movement angles of the subject’s upper limbs at first. Next, forward kinematics calculations were performed, and the desired task trajectory was obtained based on the kinematic agreement between the subject and the exoskeleton. Lastly, J-IK and CTPPG-IK were used to calculate the joint space angle and compare the performance of both. The computation of kinematics and data processing in the simulation environment were completed offline. Then, we saved each set of joint angles in .csv format in the industrial PC and used the position control mode on the exoskeleton for experimental validation. The maximum joint space error of ${\theta }_1$ and ${\theta }_2$ are $E_{{\theta }_1}$ and $E_{{\theta }_2}$. The maximum orientation error of the end effector is $E_{\omega }$.

4.3. Results and Discussion

Figure 8 shows the rotation angles of the volunteer’s shoulder girdle during different motions relative to the angles predicted by MLS-SVR. The sky blue and dark blue lines represent the measurement data of the motion capture system, the orange and red lines represent the output angles of the prediction model, respectively. Figure 9 shows the estimation accuracy of ${\theta }_1$ and ${\theta }_2$. R-square values of ${\theta }_1$ exceeded 0.8 in all three movements, indicating that the prediction model was more accurate for elevation/depression than the shoulder girdle protraction/retraction.

The R-square of ${\theta }_2$ is 0.952 in cleaning, and 0.592 and 0.621 in drinking and arm cycling, respectively. Figure 9b also shows fluctuations in the prediction of ${\theta }_2$. Through the analysis of the raw data, it could be found that ${\theta }_1$ has uniform periodic variation in the drinking and arm cycling, while ${\theta }_2$ has a significant decrease in cycling consistency. Shoulder joint elevation is small, and the motion is dominated by anterior protraction and retraction during cleaning motion, thus causing the R-square deviation of ${\theta }_2$ in each group. The ROM of arm cycling is larger than the first two sets of motions, which is a reason for its higher RMSE.

CLEVERarm developed two polynomial models relating the movements of the GH center and arm elevations, employed through elevation/depression and protraction/retraction experiments [30]. However, there is a strong correlation between the elevation/depression angle and the type of movement. In addition, that study differed from this paper in that the center of motion of the shoulder girdle was set on the axis of symmetry of the coronal plane. The polynomial models were retained, and the parameters were re-fitted using the data for the touching head action, and used as a reference as follows:

$${\phi }_{ed}=4.33\times {10}^{-3}{\epsilon }^3-6.86\times {10}^{-2}{\epsilon }^2+0.062\epsilon +0.05$$

(19)

$${\phi }_{pr}=-1.89\times {10}^{-5}{\epsilon }^4+4.53\times {10}^{-5}{\epsilon }^3-3.72\times {10}^{-2}{\epsilon }^2+0.014\epsilon -0.194$$

(20)

where $\epsilon ={\phi }_{po}$. Figure 10 shows the comparison of the MLS-SVR model and the polynomials model during the touching head movement. The figure on the left shows the prediction performance of the two models of ${\theta }_1$. The RMSEs of ${\theta }_1$ errors of the MLS-SVR and polynomial model were 0.0201 rad and 0.0242 rad, and the maximum errors were 0.0275 rad and 0.0543 rad, respectively. Both models were able to predict ${\theta }_1$ relatively accurately, with the MLS-SVR model having slightly higher precision. In the right figure, the reference angle ${\theta }_2$ was between −0.15 rad and −0.03 rad. The RMSEs of ${\theta }_2$ errors of the MLS-SVR and polynomial model were 0.0125 rad and 0.0237 rad, and the maximum errors were 0.0218 rad and 0.0565 rad, respectively. The two indicators of the polynomial model are approximately twice those of the MLS-SVR model, and the maximum error is approximately 57.6% of the range of motion. Therefore, the MLS-SVR model has a significantly superior predictive performance. When considering both dimensions of motion, the MLS-SVR model shows better prediction precision.

Figure 11 shows the joint space angle errors using the J-IK algorithm and the CTPPG-IK algorithm for a segment of the volunteer head-touching trajectory, where ${\theta }_{error}={\theta }_{id}-{\theta }_i$. These errors represent the differences between these two algorithms and the baseline joint angles ${\theta }_{id}$ mentioned in Section 3, respectively. We set the joint space convergence threshold to 3.5 ×${10}^{-3}$ rad, so the exoskeleton is guaranteed to have a small error in tracking the shoulder girdle rotation angles predicted by the humerus. In contrast, the joint angles error calculated by J-IK is several orders of magnitude higher than that of CTPPG-IK, and the absolute values of the maximum errors for ${\theta }_1$ and ${\theta }_2$ are 0.4707rad and 0.1190rad, respectively. CTPPG-IK is a closed-loop method for adjusting the shoulder girdle orientation for the next cycle based on the predicted humeral spherical coordinate angles for ${\theta }_1$ and ${\theta }_2$ in the previous sampling cycle. In contrast, the joint rotation angles calculated by J-IK are the least two parametric solutions of the solution space, and each joint needs an angle as small as possible to achieve the end-effector task.

The performance metrics comparative experiments are shown in

Table 2. Performance metrics results.

Motion	Algorithm	${\mathit{E}}_{{\mathit{\theta }}_1}$	${\mathit{E}}_{{\mathit{\theta }}_2}$	${\mathit{E}}_{\mathit{\omega }}$
Cleaning	J-IK	1.52 × ${10}^{-1}$	1.47 × ${10}^{-1}$	2.77 × ${10}^{-6}$
	CTPPG-IK	2.84 × ${10}^{-4}$	6.65 × ${10}^{-4}$	7.58 × ${10}^{-6}$
Drinking	J-IK	2.47 × ${10}^{-1}$	1.27 × ${10}^{-1}$	3.19 × ${10}^{-6}$
	CTPPG-IK	1.61 × ${10}^{-3}$	4.15 × ${10}^{-4}$	1.46 × ${10}^{-5}$
Arm	J-IK	2.98 × ${10}^{-1}$	1.79 × ${10}^{-1}$	3.66 × ${10}^{-6}$
Cycling	CTPPG-IK	2.07 × ${10}^{-3}$	3.04 × ${10}^{-3}$	1.89 × ${10}^{-5}$

. The maximum of the joint space error of J-IK is much larger than that of CTPPG-IK. The maximum values of $E_{{\theta }_1}$ and $E_{{\theta }_2}$ by CTPPG-IK are 2.07 ×${10}^{-3}$ and 3.04 ×${10}^{-3}$ rad during arm cycling, respectively, which proves that CTPPG-IK has good performance in tracking the shoulder rhythm of the subject. The maximum values of $E_{{\theta }_1}$ and $E_{{\theta }_2}$ by J-IK are 2.98 ×${10}^{-1}$ rad and 1.79 ×${10}^{-1}$ rad, representing the motion mismatch of the shoulder joints of the subject with the exoskeleton’s shoulder girdle. The joint space error of J-IK increases as the shoulder girdle moves away from the initial position due to the cumulative error caused by the least squares solution.

5. Discussion

Nowadays, researchers in the literature have significantly explored intelligent techniques to solve the non-closed form equations for IK problems of more than 4-DOFs manipulators, such as genetic algorithms [39], Bayesian optimization [40], particle swarm optimization [41], and artificial neural networks [42]. In general, evolutionary algorithms are good choices for applications where a real-time solution is needed. Artificial neural networks are a good choice for applications where a more accurate solution is needed. However, these methods could be computationally expensive and may not always converge to the global optimum solution. Additionally, neural networks need a significant amount of training data and computing resources. Compared to these methods, the algorithm proposed in this paper takes into account the constraints of scapulohumeral rhythm. By projecting the predicted shoulder angles into null space, this algorithm typically requires only 2–5 iterations to determine joint angles that satisfy the constraints. The algorithm has a well-defined mathematical formulation and can meet the requirements for real-time computation.

While we employed the DLS method to mitigate excessive joint velocities near singular points and did not observe significant issues in the current experiments, constraints on velocity aspects should still be incorporated to enhance safety during practical use. Furthermore, 7 DOF is not redundant enough for the dexterous manipulation of the upper limb. We contemplate the addition of more joints to better align the exoskeleton with upper limb movements.

6. Conclusions

We designed a 7-degrees-of-freedom upper limb rehabilitation exoskeleton, FREE, which can be used to conduct a variety of rehabilitation training and improve the synergy of movement. A MIMO prediction model was used to predict shoulder girdle motion by scapulohumeral orientation. Only the elevation angle of the upper arm is the input in the existing model. In this study, both the impact of upper arm movement on the shoulder girdle elevation/depression and the protraction/retraction were simultaneously considered to establish a MIMO model to predict shoulder girdle angles. A human-like CTPPG-IK algorithm based on differential kinematics was proposed to handle multiple constraints in exoskeleton assistance. The CTPPG-IK algorithm combines the gradient-projection method with task priority-based control, which ensures both the degree of completion of the rehabilitation motions and the natural shoulder motion. The maximum errors of ${\theta }_1$ and ${\theta }_2$ are 2.07 ×${10}^{-3}$ and 3.04 ×${10}^{-3}$ rad, respectively. And most importantly, based on CTPPG-IK, FREE can be used to customize shoulder motion to match the characteristics of wearers.

In the future, an online parameter identification method will be developed to obtain both the shoulder girdle rotation center coordinate and scapulohumeral rhythm. In addition, the configuration parameters of FREE will be optimized according to the Asian body shape to match more wearers. Currently, experiments are not being conducted with users wearing FREE. We will conduct actual rehabilitation experiments after refining the compliant control.