**1. Introduction**

Functional rehabilitation aims to recover as much as possible of the patient's locomotion independence. It requires the assistance of a therapist to perform repetitive exercises for an injured member [1]. Task-oriented protocols, where the patient is assisted to perform a specific prescribed movement, such as kicking a ball or standing up and walking, show promising outcomes compared to the conventional training based on passively moving the impaired joints in the limits of their range of motion [2].

Rehabilitation sessions can last up to several weeks [3]. In addition, to guarantee a better quality of the followed protocol, one-to-one assistance is needed [4]. However, the limited number of available therapists influences the high-intensity and the repetition of the assistance. Given these issues, researchers have developed robotic devices to assist practitioners' tasks [5–9]. They also provide the opportunity to assess the patient's recovery progress and monitor protocol efficiency; for instance, by using ARMin [10], the ARM Guide [11], and the MIME [12].

**Citation:** Ennaiem, F.; Chaker, A.; Laribi, M.A.; Sandoval, J.; Bennour, S.; Mlika, A.; Romdhane, L.; Zeghloul, S. Task-Based Design Approach: Development of a Planar Cable-Driven Parallel Robot for Upper Limb Rehabilitation. *Appl. Sci.* **2021**, *11*, 5635. https://doi.org/10.3390/ app11125635

Academic Editor: Manuel Armada

Received: 24 May 2021 Accepted: 14 June 2021 Published: 18 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The robotic devices used for clinical practice allow mainly the rehabilitation of the knee, the shoulder and the elbow [13]. Cable-driven parallel robots (CDPRs) [14] can be used to extend the coverage of existing robotic platforms since they allow rehabilitation of other joints [13]. They have also a larger translational workspace, less dynamic inertia and higher flexibility compared to serial manipulators [13]. Thanks to their simple reconfiguration and light weight, they have no setup constraints and can be used at the patients' personal spaces to reduce the need to move to rehabilitation centers.

Various studies have been carried out dealing with the design [15], control [16–18] and the structural optimization of CDPRs. In [19], Lorenzo et al. studied the optimal position of the cable exit points allowing minimization of CDPR size for fully constrained and suspended configurations. In this work, the location of the cable exit points was supposed to be fixed, which was an unrealistic assumption. In [20], Hussein et al. optimized the CDPR geometry by minimizing the maximum cable tensions. This criterion could not guarantee minimum cable tensions since only the maximum value was optimized. Abbas et al. [21] studied the optimal design of a suspended CDPR using first the workspace area, then the global condition index, which describes the robot dexterity, as two separated objective functions and, thus, only one criterion in each optimization was used. Yangmin et al. [22] studied first, the optimal design of a CDPR taking into consideration the dexterity then the elastic stiffness as criteria, then a multiobjective optimization approach was used mixing the two characteristics. The selection of each criterion weight was chosen in a way that either the system was preferred to be more dexterous or stiffer. The authors chose to give the same importance to the two criteria. Such a choice can be improved by studying the influence of the variation of each coefficient on the objective function and on the criteria.

The gaps presented above are taken into consideration in this paper. The goal is to design a planar CDPR for upper limb rehabilitation. The prescribed exercise consists of tracing the number eight shape with the hand. This form was used among others since it involves shoulder and elbow motions [23]. The rehabilitation of the wrist joint was also considered since the hand orientation was not fixed. The cable tensions, the elastic stiffness and the dexterity were selected as the problem criteria. The choice of each criterion weight was justified, and the design parameters were set in a way to obtain a realistic design. This design problem can be formulated as a constrained optimization problem. Thus, this paper aims to design an optimal planar CDPR, based on the above-mentioned criteria, able to mobilize the patient's affected upper limb along the prescribed task. The design optimization problem was formulated after considering the disadvantages of some proposed approaches introduced above. Once the optimal structure was selected, an experimental prototype was developed allowing a check of solution feasibility and reliability.

The paper is organized as follows: Section 2 details the experimental protocol followed to record the prescribed task. The planar CDPR model is presented as well as the criteria and the constraints adopted to seek the optimal design. In Section 3, the optimization results, using a multiobjective formulation then an adapted mono-objective formulation, are discussed and the experimental validation using a developed prototype is shown. The last section concludes the paper.

#### **2. Materials and Methods**

#### *2.1. Prescribed Exercise Analysis*

Rehabilitation aims to recover the functional abilities of the affected member by performing intensive and repetitive training [1]. The chosen exercise for this study was commonly performed for upper limb rehabilitation. It consisted of tracing with the hand an "8" curve. This exercise allows the rehabilitation of the three upper limb joints.

A healthy subject was asked to perform the prescribed drill in order to obtain a normal trajectory to use as a reference. This exercise involved five joints' movements, namely the three rotations of the shoulder and the flexion/extension motion of the elbow and the wrist. The participant's gestures were recorded using a Qualisys motion capture system

with five infra-red cameras and five reflective markers attached to his hand as illustrated in Figure 1. Figure 2 illustrates the steps followed for the data acquisition.

**Figure 1.** (**a**) Motion capture setup; (**b**) Marker locations.

#### **Figure 2.** Data acquisition steps.

The hand trajectory was defined by tracking the successive positions taken by the marker H3 during the prescribed exercise. Hand orientation was delimited by computing the rotation angle between a local frame (H3, **x**, **y**) attached to the participant's hand and a global frame (O, **X**, **Y**) attached to the table as illustrated in Figure 3. Since the patient's hand was attached to the robot end-effector, the recorded data, given in Figure 4, were used to define the robot's prescribed workspace.

**Figure 3.** Local and global frames.

**Figure 4.** (**a**) Hand trajectory (composed of 331 equidistant points); (**b**) hand orientation. "+" denotes the starting and the ending position.

#### *2.2. Optimal Synthesis Problem and Its Formulation*

A planar cable-driven parallel robot with three degrees of freedom (DOFs) is considered in this paper. Since at least one more cable than the DOFs was needed to fully constrain the robot, four cables were used [24]. Their mass and elasticity were neglected. They were modeled as straight lines. The design process consisted of finding the optimal position of each actuator and the end-effector size satisfying a set of criteria and constraints. The motors' positions were defined using the parameters ai and bi, which represent the coordinates of the center of the pulley P*i*, fixed on the *ith* actuator. The mobile platform was considered as a square of side c (see Figure 5). The design vector is given by Equation (1). The global and the local frames (**X**,**<sup>Y</sup>**) and (**<sup>x</sup>**, **<sup>y</sup>**), respectively, matched those represented in Figure 2, used for the motion capture analysis.

$$\mathbf{I} = [\mathfrak{a}\_1, \mathfrak{b}\_1, \mathfrak{a}\_2, \mathfrak{b}\_2, \mathfrak{a}\_3, \mathfrak{b}\_3, \mathfrak{a}\_4, \mathfrak{b}\_4, \mathfrak{c}],\tag{1}$$

n**i** denotes the unit vector along the *ith* cable. It is defined as follows, where *ψi* and *θi* are expressed as given in Equation (2). *Li* is the length of the *ith* cable and *Rp* is the pulleys' radius. *Rp* is considered to be fixed, since the cable radius and the coiling effects are neglected.

**Figure 5.** Geometric representation of the planar robot.

$$\mathfrak{n}\_{1} = \begin{bmatrix} -\sin\left(\psi\_{1} + \theta\_{1}\right) \\ -\cos\left(\psi\_{1} + \theta\_{1}\right) \end{bmatrix}, \mathfrak{n}\_{2} = \begin{bmatrix} \sin\left(\psi\_{2} - \theta\_{2}\right) \\ -\cos\left(\psi\_{2} - \theta\_{2}\right) \end{bmatrix}, \mathfrak{n}\_{3} = \begin{bmatrix} \sin\left(\psi\_{3} + \theta\_{3}\right) \\ \cos\left(\psi\_{3} + \theta\_{3}\right) \end{bmatrix}, \mathfrak{n}\_{4} = \begin{bmatrix} -\sin\left(\psi\_{4} - \theta\_{4}\right) \\ \cos\left(\psi\_{4} - \theta\_{4}\right) \end{bmatrix}, \tag{2}$$

$$\psi\_i = \tan^{-1}(\frac{|\mathcal{P}\_{ix} - \mathcal{B}\_{i\_x}|}{|\mathcal{P}\_{iy} - \mathcal{B}\_{i\_y}|}),\tag{3}$$

$$\theta\_i = \tan^{-1} \left( \frac{R\_p}{L\_i} \right),\tag{4}$$

The end-effector dynamic model was obtained using Newton-Euler formulation with the assumption of neglecting the cables mass and the dynamics of the pulleys. Its expression is given by Equation (5).

$$
\begin{bmatrix}
\sum \mathbf{F} \\
\sum \mathcal{M}
\end{bmatrix} = \mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} = \mathbf{J}^{\mathrm{T}}\mathbf{T} + \mathbf{f}\_{\mathbf{g}} + \mathbf{F}\_{\mathrm{ext}/\mathrm{EE}}\tag{5}
$$

where **M** = *mp*I3×<sup>3</sup> 03∗<sup>3</sup> 03×<sup>3</sup> **RΦ** I**p <sup>R</sup>Φ**<sup>T</sup> and **C** = 03×<sup>1</sup> ω × **<sup>R</sup>Φ**I**p <sup>R</sup>ΦT**<sup>ω</sup> are the mass and the Coriolis matrices, respectively. *mp* and ω are the mass and the angular velocity of the end-effector, **RΦ** is the rotation matrix, I**p** denotes the inertial matrix of the mobile platform written in its center of mass, χ = [*x y* Φ]<sup>T</sup> is the pose vector of the endeffector, **JT** is the transpose of the Jacobian matrix given in Equation (6), and **T**, **fg**, and **Fext**/**EE** are the cable tensions vector, the gravity force and the external forces applied on the end-effector, respectively.

$$\mathbf{J}^{\mathrm{T}} = \begin{bmatrix} \ \mathfrak{n}\_{\mathrm{i}} & \mathbf{Z} \cdot (\mathbf{R} \boldsymbol{\Phi} \cdot \mathbf{B}\_{\mathrm{i}}) \times \mathfrak{n}\_{\mathrm{i}} \end{bmatrix}^{\mathrm{T}}, \ \mathrm{i} = 1.4,\tag{6}$$

**Bi** is the vector containing the anchor points coordinates (B*i*) in the local frame.

Thus, the expression of the cable tensions vector **T** can be deduced using Equation (5) as follows:

$$\mathbf{T} = \mathbf{T}\_{\mathsf{P}} + \mathbf{T}\_{\mathsf{h}} = \left(\mathbf{J}^{\mathsf{T}}\right)^{+} \left(\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} - \mathbf{f}\_{\mathsf{g}} - \mathbf{F}\_{\mathsf{ext}/\mathsf{O}\_{\mathsf{L}}}\right) + \boldsymbol{\lambda} \text{Null}\left(\mathbf{J}^{\mathsf{T}}\right) \tag{7}$$

where **Tp** and **Th** are the particular and the homogenous solutions, **J<sup>T</sup>**<sup>+</sup> is the Moore-Penrose pseudoinverse and *λ* is an arbitrary scalar. The cable tensions must be bounded between a minimum positive and a maximum value in order to avoid slack and overtensioned cables. This condition forms the first problem constraint, which is formulated as follows:

$$0 < T\_{\rm min} \le T\_i(j) \le T\_{\rm max}, \ i = 1..4, \ j = 1..n \tag{8}$$

where *n* is the number of points composing each trajectory, *Ti*(*j*) is the tension of the *ith* cable at the *jth* position.

The second constraint concerns the collisions between the cables and the end-effector. to prevent this issue, the angle *αi* between the mobile platform and the *ith* cable is computed for each position of the prescribed trajectories. *αi* must remain higher than a limit angle *αlim*. The formulation of the collision constraint is given as follows:

$$a\_i = \cos^{-1} \left( n\_{\mathbf{i}} \cdot \frac{\mathbf{m}\_{\mathbf{i}}}{||\mathbf{m}\_{\mathbf{i}}||} \right) > a\_{\lim} \qquad i = 1..4,\tag{9}$$

where m**1** = −m**4** = *B***1***B***<sup>4</sup>**, m**2** = −m**3** = *B***2***B***<sup>3</sup>**, and n**i** is the unit vector along the *ith* cable as illustrated in Figure 5.

The last constraint concerns the location of the points P*i*. They must be located on the edges of the square forming the robot fixed frame. This constraint facilitates the robot reconfiguration when any other trajectory, included in the robot workspace, is selected.

The 4 × 3 Jacobian matrix **J** of the considered robot includes components of different physical units since the latter has mixed DOFs. To have a meaningful value of the condition number, the Jacobian matrix must be normalized allowing assessment of the closeness of a pose to a singularity. In order to settle the dimensional inhomogeneity, several methods have been suggested based on dividing the rotational elements of the Jacobian matrix **J** by a conventional length *L*. Lee et al. defined *L* in [25] as a nominal length represented by the distance between the origins of the global and the local frames. Angeles introduced in [26] the notion of natural length, which is the value of *L* that minimizes the condition number. This length is approximated to the radius of the end-effector [27]. The latter method was adopted in this paper. The homogenous Jacobian matrix, **J***<sup>h</sup>*, of size 4 × 3, is calculated as given by Equation (10), where c is the side length of the mobile platform.

$$\mathbf{J\_h} = \mathbf{J} \cdot \text{diag}\left(1, \begin{array}{c} 2 \\ \hline \text{c} \sqrt{2} \end{array}\right),\tag{10}$$

Three criteria are considered for the optimization problem formulation, namely minimizing the tension in each cable and maximizing both the elastic stiffness and the dexterity of the end-effector. The first, the second, and the third criteria are characterized using the parameters C1, C2, and C3, respectively, as follows:

$$\mathcal{C}\_1 = \frac{1}{4} \sum\_{i=1}^4 \frac{\sum\_{j=1}^n T\_{i(j)}}{n \max\_{j=1 \ldots n} T\_{i(j)}},\tag{11}$$

$$\mathcal{C}\_2 = 1 - \frac{1}{n} \sum\_{j=1}^{n} \frac{\lambda\_{\text{min}}\left(j\right)}{\lambda\_{\text{max}}\left(j\right)}\tag{12}$$

$$\mathcal{C}\_3 = 1 - \frac{1}{n} \sum\_{j=1}^n \frac{\sigma\_{\text{min}}\left(j\right)}{\sigma\_{\text{max}}\left(j\right)},\tag{13}$$

where *λ* and *σ* are the eigenvalue of the cable stiffness matrix and the singular value of the normalized Jacobian matrix **Jh**. The maximum value of the ratios *λmin λmax* and *σmin σmax* is 1, thus, maximizing the elastic stiffness and the dexterity leads to minimizing the parameters C2 and C3.

#### **3. Results and Discussions**

#### *3.1. Multiobjective Formulation*

The multiobjective formulation is given in Equation (14). A penalty formulation was adopted to handle the problem constraints.

$$\begin{cases} \min(\mathcal{C}\_1(\mathbf{I}) + \wp\_1 + \wp\_2 + \wp\_3), \\ \min(\mathcal{C}\_2(\mathbf{I}) + \wp\_1 + \wp\_2 + \wp\_3), \\ \min(\mathcal{C}\_3(\mathbf{I}) + \wp\_1 + \wp\_2 + \wp\_3), \end{cases} \tag{14}$$

$$\wp\_1 = \begin{cases} 0 & \text{if } \quad 0 < T\_{\text{min}} \le T\_i(j) \le T\_{\text{max}}\\ \psi & \text{otherwise} \end{cases},\tag{15}$$

$$\wp\_2 = \begin{cases} 0 & \text{if } \quad \alpha\_i > \alpha\_{lim} \\ \psi & \text{otherwise} \end{cases} \tag{16}$$

$$
\wp\_3 = \begin{cases} 0 & \text{if } \quad \mathcal{P}\_i \text{ é fixed frame} \\ \psi & \text{otherwise} \end{cases} \tag{17}
$$

where ℘*i*, *i*=1,2,3 are the penalty functions and *ψ* is a large scalar.

Multiple solutions coexist forming the Pareto front displayed in Figure 6. Without adding supplementary information about the desired result, all the nondominated points form a potential optimal solution.

**Figure 6.** The Pareto front.

In order to facilitate its interpretation, two dimensions representations of the Pareto front are given. Figure 7 shows the cable tension vs the elastic stiffness criteria, the cable tension vs dexterity criteria, the dexterity vs the elastic stiffness criteria and the histogram of each criterion. The three criteria are normalized, and their values vary from 0 to 1.

**Figure 7.** Bi-objective visualizations of the Pareto front.

The nondominated solutions for each two criteria present the preferable solutions if only the two corresponding objectives are considered. For the overall problem, the most desirable results are those which lay within the nondominated solutions of all the pairs illustrated in Figure 7. For the Pareto front presented above, the selection of the optimal solution was complicated since each solution minimized no more than two criteria and maximized those remaining. Additional subjective preference information was, therefore, needed to facilitate the decision-making.

#### *3.2. Mono-Objective Formulation*

Another method can be used for the resolution of a multiobjective optimization problem. This technique, called the scalarization method [28], consists of attributing a weight to each criterion that corresponds to the priority given to this objective in the optimization process. This method allows a single solution. Based on this technique, the modified optimization aims to find the optimal design vector **I**∗ which minimizes the new objective function F, defined as a weighted sum of the three subfunctions introduced above. The optimization problem is then formulated as follows:

$$\min(\mathcal{F}(\mathbf{I}))\_{\prime}$$

Subject to

$$0 < T\_{\min} \le T\_i(j) \le T\_{\max}, \ i = 1..4, \ j = 1..n \tag{18}$$

$$\alpha\_i > \alpha\_{lim} \qquad i = 1..4$$

$$\mathcal{P}\_i \text{ } \epsilon \text{ fixed frame, } i = 1..4$$

The problem constraints are handled using a penalty approach. F is then expressed as follows:

$$\mathcal{F}(\mathbf{I}) = \beta\_1 \mathcal{C}\_1(\mathbf{I}) + \beta\_2 \mathcal{C}\_2(\mathbf{I}) + \beta\_3 \mathcal{C}\_3(\mathbf{I}) + \wp\_1 + \wp\_2 + \wp\_3. \tag{19}$$

where ℘*i*, *i*=1,2,3 are the penalty functions given by the Equations (15)–(17) and *β*1, *β*2 and *β*3 define the weight of each criterion. These coefficients must verify Equation (20).

$$\sum\_{i=1}^{3} \beta\_i = 1,\tag{20}$$

Three strategies have been used in the literature to compute the values of the coefficients *βi* [29], namely, the equal weights method [30], the rank order centroid weights method [31] and the rank-sum weights method [31].

The equal weights method was used in our previous work [32]. In this paper, a deeper study of the appropriate values of these weighting coefficients is conducted in order to identify the impact of each coefficient on the quality of the final solution and the optimization process. A mapping of the global objective function and those corresponding to the three criteria is computed for all the possible coefficient values inside the interval (0,1) with respect to Equation (20), using the Particle Swarm Optimization algorithm (PSO) [33]. The computed variations are illustrated in Figure 8.

**Figure 8.** Variation of (**a**) the cable tension criterion, (**b**) the elastic stiffness criterion, (**c**) the dexterity criterion and (**d**) the objective function along the prescribed movement vs *β*1 and *β*2 (*β*3 = 1 − (*β*1 + *β*2)).

Referring to Figure 8d, the objective function decreases with *β*1. This means that the cable tension criterion C1 has the greatest impact on the objective function F compared to the other coefficients *β*2 and *β*3. Thus, having a structure with a minimum cable tension distribution is easier than generating good dexterity and elastic stiffness.

As illustrated in Figure 8c, the dexterity criterion C3 is nearly equal to 0.3 (C3 ∈ [0.3, 0.33] for *β*3 = 0 and C3 = 0.28 for *β*3 = <sup>1</sup>). Thus, according to Equation (13), the robot dexterity remains between 0.67 and 0.72. In other words, it takes good values even when *β*1 and *β*2 are large (*β*3 is close to zero). This is contrary to the elastic stiffness criterion C2,whose minimum value is equal to 0.58 for *β*2 = 1, which corresponds to an elastic stiffness equal to 0.42 (according to Equation (12)).

Based on this mapping study, since the dexterity has acceptable values for any chosen *β*3, the lowest coefficient is given to the criterion C3. The remaining weight is divided equally between the two other criteria. Thus, *β*3 = 0.1, *β*1 = *β*2 = 0.45.

Several methods exist in the literature allowing the resolution of optimization problems [34,35]. Particle Swarm Optimization (PSO) was used in this paper to select the optimal design vector. The different parameters used for the problem resolution are listed in Table 1. The lower and the upper boundaries of each design parameter are listed in Table 2.

**Table 1.** Optimization problem parameters.


**Table 2.** Boundaries of the design parameters.


The obtained solution, as well as the corresponding values of the objective function and the three criteria, are given in Table 3.

**Table 3.** Optimization results.


The maximum cable tension, the elastic stiffness and the dexterity of the optimal robot structure computed along the task workspace are illustrated in Figure 9. Figure 10 displays the free collision static workspace of the optimal structure for different values of Φ.

**Figure 9.** *Cont*.

**Figure 9.** Optimization results (**a**) the cable tension distribution, (**b**) the elastic stiffness distribution, (**c**) the dexterity distribution along the task workspace and (**d**) the arrows signification represented in (**<sup>a</sup>**–**<sup>c</sup>**).

**Figure 10.** The free collision static workspace (in blue) for (**a**) Φ = Φ*min* = <sup>−</sup>13◦, (**b**) Φ = 0◦ and (**c**) Φ = Φ*max* = 24◦.

## **4. Experimental Validation**

Following the optimization study detailed above, a CDPR prototype, presented in Figure 11a, was developed. The actuators' positions and the end-effector size were adjusted according to the optimal design vector resulting in Table 2. It was composed of a rigid frame with four actuated pulleys controlling the end-effector poses by means of four cables. The patient sits near to the robot and grabs the mobile platform with his hand. Figure 11b shows the prototype representation using QTM (Qualisys Track Manager) software recording the passive markers' locations along the end-effector movement.

**Figure 11.** (**a**) Robot prototype and (**b**) its representation using QTM (Qualisys Track Manager) software. Markers are represented in green.

The prototype was actuated using the Dynamixel MX-106T and equipped with a feedforward, a PID controller and a gear reducer. The pulleys and the end-effector were produced using 3D printing. The cables were made of Dacron.

The extended position control mode was used for the actuator commands. The desired angular positions of the motors were computed using the inverse kinematic model. Figure 12 presents the block diagram of the robot control.

**Figure 12.** Control block diagram.

The angular positions of the pulleys were controlled so that the robot end-effector could replicate the prescribed exercise. A motion capture system was also used in this phase to evaluate the performance of the robot by tracking the mobile platform poses along its movement so that a comparison, illustrated in Figure 13, between the target and the performed trajectories could be made. This external and global approach allowed

consideration of all error sources inside and outside the control loop, such as anchor point errors and the cable's behavior.

**Figure 13.** Variations between the desired and the real (**a**) trajectories and (**b**) orientations.

In order to assess the robot prototype performance, the mean absolute, the root-meansquare and the standard deviation position and orientation errors were computed as given in Table 4.


**Table 4.** Measured errors between the desired and the performed end-effector poses.

Where **Xt**(*i*) and **Xm**(*i*) are the true and the measured values corresponding to the *ith* position, (*i*) = **Xt**(*i*) − **Xm**(*i*) is the measured error at the *ith* position, is the mean of and *n* is the number of points forming the trajectory.
