Comparison of metaheuristic optimization algorithms for dimensional synthesis of a spherical parallel manipulator

Mechanism and Machine Theory 140 (2019) 586–600

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

Research paper

Comparison of metaheuristic optimization algorithms for dimensional synthesis of a spherical parallel manipulator José-Alfredo Leal-Naranjo a,∗, Jorge-Alberto Soria-Alcaraz b, Christopher-René Torres-San Miguel c, Juan-Carlos Paredes-Rojas d, Andrés Espinal b, Horacio Rostro-González a Departamento de Electrónica, DICIS, Universidad de Guanajuato, Carr. Salamanca-Valle de Santiago km 3.5+1.8, Comunidad de Palo Blanco, Salamanca, Guanajuato 36885, Mexico b DECEA, Universidad de Guanajuato, Fraccionamiento el establo 1, Guanajuato, Guanajuato 36250, Mexico c Sección de Estudios de Posgrado e Investigación, Instituto Politécnico Nacional, Unidad Profesional “Adolfo López Mateos”, Escuela Superior de Ingeniería Mecánica y Eléctrica Zacatenco, Edificio 5, 2° Piso Col. Lindavista, 07738 México City, Mexico d Laboratorio Nacional de Desarrollo y Aseguramiento de la Calidad en Biocombustibles del Centro Mexicano para la Producción más Limpia, Instituto Politécnico Nacional, Mexico City, Mexico a

a r t i c l e

i n f o

Article history: Received 27 February 2019 Revised 7 May 2019 Accepted 22 June 2019

Keywords: Multi-objective optimization Parallel manipulator NSGA-II OMOPSO MOEA/D

a b s t r a c t It is well known that the performance of parallel manipulators is highly dependent on their geometry, therefore the dimensional synthesis plays an important role in the design of a manipulator. In this work, the dimensional synthesis of a spherical parallel manipulator is carried out by means of a multi-objective optimization. This manipulator is planned to be a prosthetic wrist. The considered objective functions were the dexterity, the variability of the dexterity, the maximum required torque and the stiffness. The optimization was performed using three different metaheuristic algorithms: NSGA-II, OMOPSO, MOEA/D. After performing a set of simulations, the results show a clearly better performance of the MOEA/D algorithm. The optimized geometry of the manipulator allowed to design a virtual prototype of a prosthetic wrist. © 2019 Published by Elsevier Ltd.

1. Introduction The great mobility of the human wrist allows to the upper limb performing a wide range of activities of daily living (ADL). Some works have suggested that the implementation of the wrist movements in a prosthetic arm could improve the execution of ADL [1,2]. For these reasons, the proper design of a prosthetic wrist plays an important role. Many kinds of prosthetic wrist have been proposed, going from passive wrist with one or two degrees of freedom (DOF), to actuated wrist with one, two or three DOFs [3]. Three DOF allows to perform the three movements of the wrist, and among the existing solutions using such characteristic it can be found devices with serial and parallel configuration. The most common configuration is based on a 3-DOF serial mechanism, but generally this leads to bulky designs. On the other hand, the parallel mechanisms have the advantage of higher rigidity and the capacity of sharing the load whereas one of the main drawbacks of parallel manipulators is its reduced workspace [4] in comparison to serial manipulators.

∗

Corresponding author. E-mail address: [email protected] (J.-A. Leal-Naranjo).

https://doi.org/10.1016/j.mechmachtheory.2019.06.023 0094-114X/© 2019 Published by Elsevier Ltd.

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

587

Nomenclature _ || ⊥

actuated joint axis of a joint parallel to axis of the previous joint axis of a joint perpendicular to axis of the previous joint ∗ axis of a joint oblique to axis of the previous joint [Axi , Ayi , Azi ] coordinates of the point A in the limb i F forces applied to the end-effector fi objective function i GCI global condition index J Jacobian matrix k condition number K stiffness matrix La limb with a configuration R||R||R⊥R∗ R Lb limb with a configuration R⊥R⊥R Li_k length of the link i in the k limb Q orientation matrix of the end-effector VH instantaneous twist of the end-effector ws workspace σ standard deviation τ vector of actuators torques i joint velocity rate between the two links j-1 and j in the i limb. j−1 ω j j−1 $ j i j { j−1 $i ; l−1 $lk }

screw between the two links j-1 and j in the i limb reciprocal product between two screws

A problem related to the design of a parallel mechanism is the fact that its performance (workspace, dexterity among others) is dependent on its geometry and pose. Furthermore, the different performance measures are mutually dependent [5]. For this reason, a proper design of a parallel manipulator involves the dimensional synthesis considering all the performance indexes of interest simultaneously. Many authors have investigated the optimization of mechanism using different techniques. In [6], the weighted sum approach for multi-objective optimization is carried out to design a planar 3-DOF parallel manipulator. The objective function was established as a mixed performance of two different indexes and three parameters of the mechanism were defined. The population size was set as 30 and the maximum number of generations as 100. A parallel robot with a structure based on a Delta-like manipulator is optimized in [7]. The algorithm used is the evolutionary algorithm presented in [8]. Nine design parameters and four objective functions were selected. A population size of 90 and a maximum number of generations of 30 was set. In [9], a parallel manipulator with six degrees of freedom is synthetized using a genetic algorithm based on the Pareto front approach. Two objective functions and five design parameters were considered. The population size and the number of generations were set to 200 and 300, respectively. A multi-objective optimization based on the Control Elitist Non-dominated Sorting Genetic Algorithm (CENSGA) is carried out in [10] to synthetize a parallel manipulator. Five design parameters and three objective functions were taken into account. The population size was set as 50 and the maximum number of generations as 52. The particle swarm optimization (PSO) and the GA algorithm are implemented in [11] to optimize a compliant parallel manipulator. The weighted sum approach is used to merge three different performance indexes. Three design parameters were considered. A better performance was found using the PSO algorithm. In [5] the Non-dominated Sorting Genetic Algorithm (NSGA-II) is used to optimize the design of a 3-RRR spherical parallel manipulator (SPM). Two objective functions and five parameters were considered as well as 12 restrictions. The initial population and maximum number of generations were set as 40 and 200, respectively. A bio-inspired parallel manipulator is optimized by means of the PSO algorithm in [12]. Two objective functions and five design parameters were chosen. The number of particles was set to 80. A Delta robot is optimized using the genetic algorithm SPEA-II in [13]. Two objective functions and three design parameters were considered. The population size and the maximum number of generations were set to 200. In [14], a dynamic optimization of a Delta manipulator is carried out. The sequential quadratic programming (SQP) algorithm is implemented to solve a multi-objective problem. Two objective functions and four design parameters were considered. A parallel manipulator was designed using a multi objective optimization in [15]. The optimization was carried out using the GA Matlab® toolbox. Two objective functions and four design parameters were chosen for the optimization.

588

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

A reconfigurable 6-DOF manipulator is optimized in [16]. Three performance criteria are combined using the weighted sum approach and the optimization is carried out using the sequential simplex method. In [17], the dimensional synthesis of a parallel manipulator is carried out through the implementation of the PSO algorithm. Four design parameters were considered. The weighted sum approach was used to combine two different performance indexes. The NSGA-II algorithm is implemented in [18] to optimize a closed-loop mechanism. Two objective functions and six design parameters were considered. A population size of 200 and the same number of generations were used. In [19], the dimensional synthesis of a 3-DOF prosthetic shoulder is carried out using the CENSGA algorithm. Three objective functions and a total of 13 design parameters were taken in to account. Here, 100 individuals and 200 generations were considered. In addition to the research works of optimization developed in the last few years, a couple of investigations have tried to assess the difference in performance using different algorithms. In [20], a comparison between the SQP, GA, DE and PSO algorithms is performed. To perform the comparison, the optimization of the workspace of a serial manipulator was carried out. The SQP performed better in terms of computational time but the results are affected due to the presence of local minima, and on the other hand the metaheuristic methods overcome the issues with local minima. In [4], the dimensional synthesis of a parallel manipulator is carried out using two different algorithms, GA and SQP. Only one objective function and four design parameters were considered. The algorithms were implemented using a Matlab® toolbox. As shown in the work of [20], the SQP performed better in terms of computation time but the results are influenced by the initial guessing and sometimes this leads to divergence. Considering the aforementioned, in this work, a comparison between three different optimization algorithms is performed with a case study of the dimensional synthesis of a spherical parallel manipulator that will be used as a prosthetic wrist. This research is organized as follows: Section 2 a brief description of the used algorithms is presented. Section 3 shows the architecture of the manipulator considered in this work and the formulation of the objective functions. In Section 4, the analysis and comparison of the algorithms is carried out, besides the synthetized manipulator and its application is shown. Finally, in Section 5 present our conclusions. 2. Algorithms description 2.1. Non-dominated Sorting Genetic Algorithm (NSGA-II) The Non-dominated Sorting Genetic Algorithm II, proposed by Deb et al. [21], is a Multiple Objective Optimization (MOO) algorithm and is an instance of an Evolutionary Algorithm. NSGA-II is an extension of the Genetic Algorithm for multiple objective function optimizations. NSGA-II is a fast and robust elitist sorted based algorithm. This well-known multi objective algorithm uses a selection operator that creates a pool of solutions by combining parent and offspring populations to select the best (non-dominated) solutions through its fast non-dominated sorting procedure to ranking solutions. One solution dominates another if it is no worse than others in all objectives and it is strictly better that other at least in one objective. The NSGA-II Algorithm also uses a crossover and mutation operator similar to a Genetic Algorithm in order produce solutions for further iterations also it has a diversity preserving mechanism which allows solutions to cover a wide range of the Pareto front. 2.2. Multiple Objective Particle Swarm Optimization (MOPSO) Multiple Objective Particle Swarm Optimization proposed by Coello and Salazar [22] is an algorithm based on the Particle swarm optimization proposed by Kennedy et al. [23], which uses the concept of Pareto dominance to determine the trajectory of a particle over the fitness landscape. This algorithm maintains previously non-dominated solution in a global repository to keep the best found solution and best local solution per particle. This algorithm have been successfully used in path planning optimization [24,25] MOPSO uses a population of solutions measured by a multiple objective fitness function. Each solution keeps a local repository as well as a global repository with the best non-dominated solution. These local and global solutions are used to guide the search trajectory of each particle according to a speed function Eq. (1) and a position function Eq. (2)

Vit+1 = wVit + c1 r1 (xLbest − Xit ) + c2 r2 (xGbest − Xit )

(1)

Xit+1 = Xit + Vit+1

(2)

where w is the inertia weight, c1 is the cognitive acceleration coefficient, c2 is the social acceleration coefficient, r1 and r2 are random values between 0 and 1. xLbest is the personal best of the particle and xGbest is the best known solution. Xit+1 is the current solution (position) of the ith particle at iteration t. 2.3. Multi-objective Evolutionary Algorithm based on Decomposition (MOEA/D) Multi-objective Evolutionary Algorithm based on Decomposition is an algorithm based on the notion of decomposing a multi-objective problem into several scalar single objective problems. The objective of each subproblem is then considered

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

589

Fig. 1. Main parameters of the manipulator. (a) Parameters of limb La . (b) Parameters of limb Lb. (c) Architecture of the manipulator.

in a linear or not linear aggregation function with all the objectives of the original multi-objective problem. Each subproblem is simultaneously optimized using mainly the information of its neighboring subproblems. These subproblems may require different amounts of computational resources so it is reasonable to apply a strategy to share efficiently the available computational resources to each part of the problem. Experimental results have demonstrated that MOEA/D with simple decomposition methods can outperforms or performs similarly to NSGA-II on multi-objective knapsack problems [26]. The main idea of MOEA/D is to decompose a multiple objective problem into N scalar optimization subproblems. This approach solves single objective problems by means of the evolution of a population of solutions. At a given iteration the population contains the best solution found so far since the start of the algorithm for each subproblem. Each subproblem is optimized using information of neighboring subproblems. This strategy allows MOEA/D to have lower computational complexity when calculating each iteration than NSGA-II because the use of a Tchebycheff decomposition approach. Importantly, normalization techniques can be fit into MOEA/D when dealing with disparately scaled objectives. 3. Formulation of the problem The parallel manipulator considered in this work is an asymmetrical parallel spherical manipulator. The SPM is made up by three limbs, both of them having the same configuration. The first limb (named La ), have the configuration R||R||R⊥R∗ R and the remaining one has the configuration R⊥R⊥R (Lb ), where R stands for revolute joint, _ stands for actuated joints, and ||, ⊥ and ∗ means that the axes of the joints are parallel, perpendicular and oblique, respectively. In the limb La , the first three revolute joints are parallel, and the axis of the joint is parallel to the unit vector ui , the fourth joint is perpendicular to the previous one, and the axes of the fourth and fifth revolute joints converges to the same center of rotation, as shown in Fig. 1a. The limb Lb consists of three revolute joints whose axes of rotation converges to the same point, as shown in Fig. 1b. The mechanism is generated assembling two limbs La and one limb Lb . The axis of the distal joint of each limb converges in the same point and forms the end-effector. The axes of the first revolute joint of both limbs La are orientated in such a way that they are perpendicular to the first revolute joint of the limb Lb . This configuration allows to obtain a mechanism with a reduced volume, Fig. 1c. The inverse kinematics problem involves the computation of the coordinates of the actuated joints knowing beforehand the orientation of the end-effector. For each distal joint of the limbs, unit vector wi shown in Fig. 1, is defined as

wi = Q wi 

(3)

where wi ’ is the initial direction of the vector when the manipulator is in its home configuration and Q is the orientation matrix of the end-effector. Referring to the limb La , for the fourth link this relation follows

vi  wi = cos(αi )

(4)

Since the unit vector vi lies on the plane YZ (Fig. 1a), the component of this unit vector in the X axis is 0, therefore the magnitude of vi is expressed as

vY 2 + vZ 2 = 1

(5)

590

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

Substituting vY from Eq. (5) into Eq. (4) and making algebraic simplifications it yields to

vZ =

−C2 +



C2 2 − 4C1C3 2C1

(6)

where

C1 = wY 2 + wZ 2

(7)

C2 = −2 cos(αi )wZ

(8)

C3 = cos2 αi − wY 2

(9)

Then, the first two links of the limb La can be analyzed as a 2R serial chain and the solution is given by

θi = 2 arctan(

−b ±

√ b2 − 4ac ) 2a

(10)

where

a = [(L3_i + ri )vZ − AZ ] + [(L3_i + ri )vY − AY ] + L1_i 2 − L2_i 2 + 2L1_i [(L3_i + ri )vZ − AZ ]

(11)

b = −4L1_i [(L3_i + ri )vY − AZ ]

(12)

2

2

c = [(L3_i + ri )vZ − AZ ] + [(L3_i + ri )vY − AY ] + L1_i 2 − L2_i 2 − 2L1_i [(L3_i + ri )vZ − AZ ] 2

2

(13)

where AY and AZ represents the coordinates of the first revolute joint. Regarding the limb Lb , for the second link, this relation follows

v3  w3 = cos(α3 )

(14)

Furthermore, the unit vector v3 can be expressed as



v3 =

cos θ3 0 − sin θ3

0 1 0

sin θ3 0 cos θ3





−1 0 0



=

− cos θ3 0 sin θ3



Substituting v3 into Eq. (14), considering the relations sin θ = fications we get



θ3 = 2 arctan

−B ±

(15) 2 tan θ2

1+tan2 θ2

and cos θ =

1−tan2 θ2 1+tan2 θ

, and making algebraic simpli-

2

 √ B2 − 4AC 2A

(16)

where

A = w3x − cos α3

(17)

B = 2 w 3z

(18)

C = −w3x − cos α3

(19)

With the previous analysis, it is possible to obtain in a closed form the coordinates of the actuated joints. 3.1. Mobility analysis In order to obtain the input-output velocity equation of the manipulator, screw theory approach is applied in this work [27]. Fig. 2 shows the screws associated to the kinematic joints in the limbs of the parallel manipulator. To generate fullrank matrices, in this work, pseudo revolute joints were added to the limbs; two pseudo revolute joints for the limb La and three for the limb Lb The instantaneous twist VH of the end-effector can be expressed as a linear combination of five instantaneous twist as o

ω1i o$1i + 1 ω2i 1 $2i + 2 ω3i 2 $3i + 3 ω4i 3 $4i + 4 ω5i 4 $5i = V H

for i = 1:3. Where limb.

i j−1 ω j

and

j−1 $ j represents i

(20)

the joint velocity rate and the screw between the two links j-1 and j in the i

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

591

Fig. 2. Screws of the kinematic joints. (a) Limb La . (b) Limb Lb .

Taking into account that the mechanism is a spherical manipulator, the three terms of the translational velocity of the end-effector vanish and its velocity state has the form

VH =



ωy

ωx

ωz

0

0

0

T

(21)

In order to solve the velocity analysis, it is needed to eliminate the intermediate joint velocity rates. For the limb La the calculation of an intermediate velocity it is needed. The intermediate velocity chosen is 1 ω2i The calculation of 1 ω2i is done applying the reciprocal product between both sides of the Eq. (20) and a reciprocal screw [28], in this case, the line vector that goes from point E to C. This line vector is expressed as:

$Bi = [(RC − RE )/|RC − RE |; RE × (RC − RE )/|RC − RE |]

(22)

And applying the reciprocal product it gives





V H ; $Bi =



o



ω1i o$1i ; $Bi + 1 ω2i 1 $2i ; $Bi

Besides, the screw 3 $4i intersects Eq. (20) with the screw 3 $4i leads to





V H ; 3 $4i =



o

the screws

(23) 2 $3 i

and

4 $5 . i

Then applying the reciprocal product to both sides of the



ω1i o$1i ; 3 $4i + 1 ω2i 1 $2i ; 3 $4i

(24)

Solving Eqs. (23) and (24) for 1 ω2i yields to 1

ω2i =

{V H ; 3 $4i } − {0 ω1i 0 $1i ; 3 $4i } {V H ; $Bi } − {0 ω1i 0 $1i ; $Bi } = {1 $2i ; 3 $4i } {1 $2i ; $Bi }

(25)

Eq. (25) can be rewritten as 0

ω1i {0 $1i ; $Ai } = {V H ; $Ai }

(26)

where

$Ai = 3 $4i {1 $2i ; $Bi } − $Bi {1 $2i ; 3 $4i }

(27)

For the remaining limb, let us consider the line vector that goes from C to B

$C3 = [(RB − RC )/|RB − RC |; RB × (RB − RC )/|RB − RC |]

(28)

592

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

Fig. 3. Wrist movements. (a) Pronation-supination. (b) Flexion-extension. (c) Ulnar-radial deviation. Table 1 Range of motion design requirements. Movement

Range of motion

Supination-pronation Flexion-extension Ulnar-radial deviation

−55°/15° −40°/40° −20°/20°

Knowing that for the limb Lb the velocities 3 ω4i and 4 ω5i are 0, and applying the reciprocal product to Eq. (20) it gives 0

ω13 {0 $13 ; $C3 } = {V H ; $C3 }

(29)

Furthermore, as the end-effector only presents rotations around a fixed center of rotation, the following relations can be established

{V H ; $4 } = vx = 0

(30)

{V H ; $5 } = vy = 0

(31)

{V H ; $6 } = vz = 0

(32)

where $4 = [1 0 0 ; 0 0 0] , $ 5 = [ 0 1 0 ; 0 0 0], and $6 = [0 Eq. (26) (with i = 1, 2), and (29)–(33), can be rewritten in matrix form as

Jq q˙ = JxV H where

0

1

;

0

0

0]

(33)

⎡0 1 A { $1 ; $1 } ⎢ 0 ⎢ 0 Jq = ⎢ ⎢ 0 ⎣

{

0 0 1 $2 ; $A2

}

0 0 0 0

0 0

⎡ T ⎤ $A1 ⎡ ⎢ A T ⎥ 0 ⎢ $2 ⎥ 0 ⎢  T ⎥⎢ ⎢ $C3 ⎥⎢0 ⎥⎢ Jx = ⎢ ⎢$ T ⎥⎢1 ⎢  4  ⎥⎣ ⎢ T⎥ 0 ⎣ $5 ⎦ 0  T

0 0 0 0 1 0

0 0 0 0 0 1

0 0 {0 $1i ; $C3 } 0 0 0

1 0 0 0 0 0

0 0 0 1 0 0

0 0 0 0 1 0

⎤

0 0⎥ 0⎥ ⎥ 0⎥ ⎦ 0 1

(34)

⎤

0 1 0 0 0 0

0 0⎥ 1⎥ ⎥ 0⎥ ⎦ 0 0

(35)

$6

q˙ =

 0

ω11

0

ω12

0

ω13

0

0

0

T

(36)

Eq. (33) represents the input-output velocity equation that allows to calculate the instantaneous velocity of the endeffector knowing the actuators velocity or vice versa. Furthermore, the Jacobian matrix can be calculated as

J = Jq −1 Jx

(37)

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

593

Fig. 4. Hyper-volume for a two objective problem considering the origin as reference point. Table 2 Two-way ANOVA F, pairwise t with Bonferroni correction and Tukey HSD tests for Algorithms using 10 0,0 0 0 function calls. Pr(>F) ANOVA t Tests NSGA-II – MOEAD NSGA-II – OMOPSO MOEAD – OMOPSO

<2.2e−16

Tukey HSD NSGA-II – MOEAD NSGA-II – OMOPSO MOEAD – OMOPSO

p adjusted 0.00 0.0 2.3e−15

2.2e−16 1.9e−6 2.2e−16

3.2. The objective functions One of the requirements for the mechanism is that it must be able to perform the movements of the wrist, Fig. 3. According to [29] and [30], the range of motion of the wrist needed to perform ADL is −40°/40° for wrist flexion-extension, 40° for ulnar and radial deviation combined, and for pronation-supination is −53/13. For this reason, a restriction imposed to this mechanism is that it can be able to perform the motions listed in Table 1. Even though, the manipulator would be able to perform a broader range of motion, for the sake of simplicity the movements and its ROMs considered in Table 1 will be considered as the workspace of the manipulator, so, the performance indicators considered in this work are only calculated in these points. In this work, a simplification is done considering the unit vector of the axis of rotation of the actuators as: u1 = [0,1,0], u2 = [1,0,0] and u3 = [0,0,1]. Considering the reference system shown in Fig. 3, the SPM must be able to perform a rotation of −20/20 around Y axis (ulnar-radial deviation), a rotation of −40/40 around X axis (flexion-extension) and a rotation of −55/15 around Z axis (supination-pronation). A widely used indicator of the kinematic performance of a manipulator is the condition number. The condition number provides a measure of the accuracy of the end-effector velocity, and its reciprocal represents a measure of the manipulator dexterity. The condition number is given by:





k = J J −1 

(38)

where ||J|| represents the norm of the matrix J, that in this work is considered as:

J  =



tr (JW J T )

(39)

and W is the identity matrix multiplied by the scalar 1/n, where n is the dimension of the matrix J. Since the condition number is a local property, the global condition index (GCI) was proposed in [31]. The GCI is based on the distribution of the condition number over a workspace ws and it is calculated as:



GCI = ws ws

1 dws k

dws

(40)

594

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

Fig. 5. Hyper-volume box plot for all algorithms with 50,0 0 0 and 10 0,0 0 0 function calls. Table 3 Optimized parameters of the parallel manipulator. Param.

Value [mm]

Param.

Value

Param.

Value

Param.

Value

L1_1 L2_1 L3_1

14.09 mm 39.96 mm 10.30 mm 83.84°

Ax1 Ay1 Az1

33.80 mm 0 38.38 mm

L1_2 L2_2 L3_2

19.11 mm 30.44 mm 11.38 mm 90°

Ax2 Ay2 Az2

0 34.55 mm −31.57 mm

α1

α2

And in a discrete form, the GCI is given by

nws

f1 = GCI =

1 j=1 k j

nws

(41)

where nws is the number of points considered in the workspace. The index 1/k is bounded to takes values in range [0,1] where 0 represents a singular condition and 1 is an isotropic condition. Therefore, the global condition index should be maximized. Therefore, the first objective function is the GCI. Part of the success of a prosthetic wrist is related to its actuators, since a compact size is wanted. It is desired a minimization of the maximum required torque by the actuators with the aim of using low size actuators. In order to evaluate the required torque a couple of assumptions are made. First, if the accelerations are low, then the inertial effects can be neglected. Furthermore, if the joints are considered frictionless, the principle of virtual work can be applied in a straightforward way. Using the virtual work approach, the external forces applied to the end-effector are given by [27]:

F = JT τ

(42)

thus, the motor torques are calculated as

τ = (JT )−1 F

(43)

The external forces applied to the end-effector considered in this work are as follows: When performing the pronationsupination movement, a torque of 0.1 Nm around the axis of pronation-supination was applied. This torque is greater than the value of 0.06 Nm required to perform ADL [32]. In the flexion-extension movement, a torque of 0.6 Nm was applied along the axis of flexion-extension. This value represents the sum of the torque of 0.35 Nm required for ADL plus the torque due to the weight of a prosthetic hand. In the ulnar-radial deviation movement a torque of 0.35 Nm was applied along the axis of ulnar-radial deviation that is the torque required for ADL and the torque due to the weight of the prosthetic hand along the axis of flexion-extension. The torques of the actuators was calculated using Eq. (43) for the different positions when the SPM performs the three required movements. The second objective function (f2 ) was defined as the value of the maximum required torque over these positions, and it needs to be minimized. The discrete calculation of the GCI is similar to the mean value of the CI and it does not provide information about regions of poor performance. For the third objective function, an indicator of the distribution of the CI is considered. In [33], an index named Global level and distribution ratio index is proposed. This index penalizes the non-uniformity of the index over the workspace. In [10], the maximum gradient of the dexterity is considered to evaluate the variation of the

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

595

Fig. 6. Evolution of the global condition index (f1 ) and the maximum required torque (f2 ) for different number of iterations.

dexterity (and the condition number) over the workspace. The gradient is evaluated numerically because there is no analytic expression to do it. In this work it is used the standard deviation of the CI, which is given by:

   2 nws  1  1  f3 = σ = − CI nws

j=1

kj

(44)

where CI is the mean of the CI. In this way, it is easily assessed the variability of the condition number. Since a lower value represents a more uniform behavior, f3 should be minimized. Another important feature commonly considered for the design of a manipulator is the stiffness. The stiffness of a parallel mechanism can be described by the stiffness matrix. Many approaches have been proposed to evaluate the stiffness of a manipulator, here it is used the method that relies in the Jacobian matrix. The stiffness matrix if given by:

K = J T KJ J

(45)

596

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

where KJ is the joint stiffness matrix of the parallel mechanism. If each actuator in the parallel mechanism is idealized as an elastic component with stiffness ki , then KJ = diag[k1 ,..,kn ]. In the case where all the actuators have the same stiffness, kj , Eq. (45) is simplified as:

K = k j JT J

(46)

Similarly as implemented in [34], to assess the stiffness of the manipulator, the global stiffness index (GSI) that is the inverse of the condition number of the stiffness matrix over ws and divided by the volume of the ws. The GSI is given by:





kk = K K −1 

(47)

nws f4 = GSI =

1 j=1 kk j

nws

(48)

Considering that in this case, the stiffness needs to be as high as possible, f4 should be maximized. 3.3. Optimization problem After the objective functions have been established, the optimization problem can be stated as

min(− f1 , f2 , f3 , − f4 )

(49)

In order to obtain a low-volume design, the parameters were constrained to take a value from an established range. Taking into account that the unit vector of the joint of the first actuator is given by u1 = [0,1,0], therefore this joint is located in the XZ plane. The distance of the actuator location and the center of rotation of the SPM in the X direction (Ax1 ) is constrained to take values between 25 and 35 mm to avoid having a wrist with a big diameter. The distance of the actuator location and the center of rotation of the SPM in the Z direction (Az1 ) is constrained to take values between −45 and −20 mm to avoid having a long wrist. In a similar fashion, the location of the second actuator was constrained. The links lengths were constrained in order to avoid collisions among them. Therefore, the optimization was subjected to the following design parameters constraints: 10 ≤ L1_1 ≤ 20, 10 ≤ L2_1 ≤ 40, 10 ≤ L3_1 ≤ 20, 50° ≤ α 1 ≤ 85°, 25 ≤ Ax1 ≤ 35, Ay1 = 0, −45 ≤ Az1 ≤ −20, 10 ≤ L1_2 ≤ 20, 10 ≤ L2_2 ≤ 40, 10 ≤ L3_2 ≤ 20, α 2 = 90°, Ax2 = 0, 25 ≤ Ay2 ≤ 35, −45 ≤ Az2 ≤ −20. Where Ax, Ay and Az are the coordinates that represent the location of the actuators. All dimensions are expressed in mm 4. Results The proposed optimization problem was solved in the Jmetal Framework [35]. This framework allowed us to implement the three algorithms in order to compare and assess their performances. For the NSGA-II implementation it was used a Simulated Binary Crossover (SBX) as main intensification operator. This operator produces two offspring from a given pair of parents. The difference between offspring and parents depends on an index ηc . A big value in "ηc " gives a higher probability of creating an offspring similar to a parent and a small value produces an offspring distant to its parent. Preliminary results show that 0.95 was a good value for this parameter. We use Polynomial Mutation as mutation or diversification operator with the value of 0.15 and Binary Tournament as selection method i.e., a higher fitness is assigned to individuals located on a sparsely populated Pareto front. In every iteration, the N existing individuals (parents) generate N new individuals (offspring). Both, parents and offspring compete with each other for inclusion in the next iteration. Finally, our population size was 100 across 40 independent experiments and the number of iterations applied were 10,0 0 0, 50,0 0 0 and 10 0,0 0 0. For the MOPSO implementation in Jmetal, our parameters were: w = 0.3, c1 = 1.1, c2 = 1.5, the population size was 100. The experiment was executed 40 times and the number of iterations applied were 10,0 0 0, 50,0 0 0 and 10 0,0 0 0. Like previous algorithms, it was used the MOEA/D implementation in Jmetal framework. For this configuration, we used q decomposition in 4 subproblems (each one for a given objective) and a population of 100 individuals. Finally, we executed 40 independent experiments with 10,0 0 0, 50,0 0 0 and 10 0,0 0 0 iterations. The hyper-volume was used as a metric of performance to evaluate the algorithms the hyper-volume indicator (also known as Lebesgue Measure or S-Metric) [36] is a metric applied to a pareto hyper-front to measure the generated area under a set of non-dominated solutions and a reference point (in this work the origin of the hyper plane). This metric is popular since it produces a single real value from an n-dimensional pareto front. This metric calculates the normalized volume of the target space  covered by the Pareto set (P∗ ) and limited by a reference point r ∗ = (r1∗ , r1∗ , ..., rn∗ ) ∈ . Fig. 4 shows an example in two dimensions of a Hyper-Volume. In order to compare the performance of the algorithms, we used several statistical tests over a set of 40 independent experiments obtained from the execution of each algorithm. The objective of this comparison is to identify if there exists a difference in the performance between two or more algorithms with statistical significance. It is important to address that these statistical tests were applied to the set of Hyper-volume values generated for every Pareto hyper-front obtained by means of the application of a multi-objective algorithm to our problem. The advantages to do this is that we have a single normalized value that represents the quality of a whole pareto hyper front.

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

597

Fig. 7. Evolution of the objective function.

As a first step, we apply a Shapiro and Wilk [37] test to determine normality of the samples. Such test showed that data came indeed from normal distributions. This is an important step in order to identify the type of further statistical test to use. Once we determined that our data can be considered as normal we apply an Analysis of variance (ANOVA [38]) to determine if, firstly, implementing different algorithms impacts on the design of the parallel manipulator, and secondly, to identify which of these methodologies offers the best result. Table 2 shows the results obtained by one-way ANOVA test, observing as independent variables the number of algorithms in this work. ANOVA’s null hypothesis (H0 ) considers that all samples belong to one unique normal distribution. As P-values (Pr(>F) in Table 1) are lower than the significance value of 0.05, there is not enough evidence to accept H0 i.e., results come from different distributions, this means that the performance of the algorithms are different with statistical significance. Pairwise t and Tukey HSD [39] test were applied next. As in ANOVA test, null hypothesis in both tests assumes that two-given samples share a single distribution (Two-given algorithms have the same performance under our configuration.). Table 1 shows t-test P-values with a Bonferroni correction. Based on these results, it can be seen that every given pair of algorithms produce a P-value less than 0.05, this evidence rejects the null hypothesis. We can conclude with statistical significance that the algorithms have different performance when dealing with the problem of the dimensional synthesis of the parallel manipulator. Once the previous results were found, higher performance for MOEAD algorithm is noticeable in the three leftmost configurations showed in Fig. 5 hyper-volume performance box plots. The hyper-volume is defined as the area under the

598

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

Fig. 8. Prototype of the prosthetic wrist. (a) Optimized spherical parallel manipulator (In red limb La1 , in green limb La2 and in blue limb Lb ). (b) Prosthetic wrist fitted with a hand (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .

Pareto hyper front, this means that the algorithm that produces the bigger hyper-volume (and the Maximum quality nondominated solutions) is MOEAD. Fig. 6 shows two comparisons of the evolution of the objective functions f1 and f2 using the three algorithms for different number of iterations. It can be seen that for each case, the solutions converge faster to the Pareto front using the MOEA/D, however the solutions are spread over a smaller range than with the other algorithms. When it is used 10,0 0 0 iterations, the NSGA-II is farther from the solution than the other algorithms. As the number of iterations increases, the three algorithms exhibit solutions closer to the Pareto front. Fig. 7 shows the impact of the maximum number of iterations when using the MOEA/D algorithm. The objective function f1 was plotted vs the other objective functions in three different charts. From the graphs, it can be concluded that the algorithm converges pretty fast to the solution since there is little to none changes as the number of iterations increases. After performing the experiments and comparing the solutions, a suitable solution was chosen. The considered solution consists of the optimized parameters shown in Table 3. This solution yields to the following values of the objective functions: f1 = 0.7934, f2 = 0.3477, f3 = 0.1088, f4 = 0.4752. Fig. 8 shows a virtual prototype of the prosthetic wrist designed using the software SolidWorks®. It can be seen that due to the reduction of the torque achieved with the optimization, the mechanism can be fitted inside a structure with reduced volume, Fig. 8b. The maximum torque required by the actuators is about 0.35 Nm. 5. Conclusions This paper has contrasted several Multi-objective algorithms carrying out the optimization of a spherical parallel manipulator. Our results indicate that the best performance was obtained by MOEA/D algorithm (Fig. 5), which is based on decomposition. The state of art has demonstrated that MOEA/D has lower computational complexity than NSGA-II and that the usage of simply decomposition methods have outperformed or performed similarly to NSGA-II test instances. Importantly, MOEA/D has reported being able to generate a very uniform distribution of representative Pareto optimal solutions on continuous multi-objective problems. The better performance of MOEA/D in dimensional synthesis of mechanisms is an important outcome since to the best knowledge of the authors, this is the first time that the MOEA/D algorithm has been used for this purpose, whereas NSGA-II and PSO have been widely used. The design of a prosthetic wrist with reduced volume has been designed. The mechanism of the wrist and its actuators can be allocated inside the forearm. This device was designed considering the range of motions and torques that a wrist needs to perform ADLs. In further work, the use of small actuators could allow to obtain a functional device with low weight.

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

599

Funding statement This research has been supported by the CONACYT project FC2016-1961 “Neurociencia Computacional: de la teoría al desarrollo de sistemas neuromórficos”. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmachtheory. 2019.06.023. References [1] F. Montagnani, M. Controzzi, C. Cipriani, Is it finger or wrist dexterity that is missing in current hand prostheses? IEEE Trans. Neural Syst. Rehabil. Eng. 23 (4) (2015) 600–609, doi:10.1109/TNSRE.2015.2398112. [2] T. Bertels, T. Schmalz, E. Ludwigs, Objectifying the functional advantages of prosthetic wrist flexion, JPO J. Prosthetics Orthot. 21 (2) (2009) 74–78, doi:10.1097/JPO.0b013e3181a10f46. [3] N.M. Bajaj, A.J. Spiers, A.M. Dollar, State of the art in prosthetic wrists: commercial and research devices, in: Proceedings of the IEEE International Conference on Rehabilitation Robotics 2015–Septe, 2015, pp. 331–338, doi:10.1109/ICORR.2015.7281221. [4] G. Liu, Y. Chen, Z. Xie, X. Geng, GA\SQP optimization for the dimensional synthesis of a delta mechanism based haptic device design, Robot. Comput. Integr. Manuf. 51 (2018) 73–84, doi:10.1016/j.rcim.2017.11.019. [5] G. Wu, Multiobjective optimum design of a 3-RRR spherical parallel manipulator with kinematic and dynamic dexterities, J. Model. Identif. Control. 33 (2012) 111–122, doi:10.4173/2012.3.3. [6] Y. Li, Q. Xu, GA-Based multi-objective optimal design of a planar 3-DOF cable-driven parallel manipulator, in: IEEE International Conference on Robotics and Biomimetics, IEEE, 2006, 2006, pp. 1360–1365, doi:10.1109/ROBIO.2006.340127. [7] E. Courteille, D. Deblaise, P. Maurine, Design optimization of a delta-like parallel robot through global stiffness performance evaluation, in: 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2009, 2009, pp. 5159–5166, doi:10.1109/IROS.2009.5353906. [8] C.M. Fonseca, P.J. Fleming, Multiobjective optimization and multiple constraint handling with evolutionary algorithms - Part I: a unified formulation, IEEE Trans. Syst. Man, Cybern. Part A Syst. Hum. 28 (1) (1998) 26–37, doi:10.1109/3468.650319. [9] Z. Gao, D. Zhang, Y. Ge, Design optimization of a spatial six degree-of-freedom parallel manipulator based on artificial intelligence approaches, Robot. Comput. Integr. Manuf. 26 (2) (2010) 180–189, doi:10.1016/j.rcim.20 09.07.0 02. [10] F.A. Lara-Molina, J.M. Rosario, D. Dumur, Multi-objective optimization of Stewart-Gough manipulator using global indices, 2011 IEEE/ASME Int. Conf. Adv. Intell. Mechatron. IEEE (2011) 79–85, doi:10.1109/AIM.2011.6026996. [11] Y. Yun, Y. Li, Optimal design of a 3-PUPU parallel robot with compliant hinges for micromanipulation in a cubic workspace, Robot. Comput. Integr. Manuf. 27 (6) (2011) 977–985, doi:10.1016/j.rcim.2011.05.001. [12] D. Zhang, Z. Gao, Forward kinematics, performance analysis, and multi-objective optimization of a bio-inspired parallel manipulator, Robot. Comput. Integr. Manuf. 28 (4) (2012) 484–492, doi:10.1016/j.rcim.2012.01.003. [13] R. Kelaiaia, O. Company, A. Zaatri, Multiobjective optimization of a linear delta parallel robot, Mech. Mach. Theory. 50 (2012) 159–178, doi:10.1016/j. mechmachtheory.2011.11.004. [14] Y. Zhao, Dynamic optimum design of a three translational degrees of freedom parallel robot while considering anisotropic property, Robot. Comput. Integr. Manuf. 29 (4) (2013) 100–112, doi:10.1016/j.rcim.2012.11.004. [15] Z. Chi, D. Zhang, L. Xia, Z. Gao, Multi-objective optimization of stiffness and workspace for a parallel kinematic machine, Int. J. Mech. Mater. Des. 9 (3) (2013) 281–293, doi:10.1007/s10999- 013- 9219- 9. [16] G. Coppola, D. Zhang, K. Liu, A 6-DOF reconfigurable hybrid parallel manipulator, robot, Comput. Integr. Manuf. 30 (2) (2014) 99–106, doi:10.1016/j. rcim.2013.09.011. [17] X.J. Wan, Q. Li, K. Wang, Dimensional synthesis of a robotized cell of support fixture, Robot. Comput. Integr. Manuf. 48 (2017) 80–92, doi:10.1016/j. rcim.2017.03.001. [18] A. Hassan, M. Abomoharam, Modeling and design optimization of a robot gripper mechanism, Robot. Comput. Integr. Manuf. 46 (2017) 94–103, doi:10. 1016/j.rcim.2016.12.012. [19] J.-A. Leal-Naranjo, M. Ceccarelli, C.-R. Torres-San Miguel, L.-A. Aguilar-Perez, G. Urriolagoitia-Sosa, G. Urriolagoitia-Calderón, Multi-objective optimization of a parallel manipulator for the design of a prosthetic arm using genetic algorithms, Lat. Am. J. Solids Struct. 15 (3) (2018) e26. [20] P.R. Bergamaschi, S. de, F.P. Saramago, L. dos, S. Coelho, Comparative study of SQP and metaheuristics for robotic manipulator design, Appl. Numer. Math. 58 (9) (2008) 1396–1412, doi:10.1016/j.apnum.2007.08.003. [21] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput. 6 (2) (2002) 182–197, doi:10.1109/4235.996017. [22] C.A. Coello Coello, M.S. Lechuga, MOPSO: a proposal for multiple objective particle swarm optimization, in: Proceedings of the 2002 Congress on Evolutionary Computation, CEC, 2002, 2002, doi:10.1109/CEC.2002.1004388. [23] J. Kennedy, Y. Shi, Swarm intelligence. The Morgan Kaufmann series in evolutionary computation, Elsevier Science & Technology, 2001. [24] Y. Zhang, D.W. Gong, J.H. Zhang, Robot path planning in uncertain environment using multi-objective particle swarm optimization, Neurocomputing 103 (2013) 172–185, doi:10.1016/j.neucom.2012.09.019. [25] N. Geng, D.W. Gong, Y. Zhang, PSO-based robot path planning for multisurvivor rescue in limited survival time, Math. Probl. Eng. 2014 (2014) 1–10, doi:10.1155/2014/187370. [26] Q. Zhang, H. Li, MOEA/D: a multiobjective evolutionary algorithm based on decomposition, IEEE Trans. Evol. Comput. (2007), doi:10.1109/TEVC.2007. 892759. [27] L.W. Tsai, Robot analysis: the mechanics of serial and parallel manipulators, John Wiley & Sons, 1999. [28] S.A. Joshi, L.-W. Tsai, Jacobian analysis of Limited-DOF parallel manipulators, J. Mech. Des. 124 (2002) 254, doi:10.1115/1.1469549. [29] J. Ryu, W.P. Cooney, L.J. Askew, K.N. An, E.Y.S. Chao, Functional ranges of motion of the wrist joint, J. Hand Surg. Am. 16 (3) (1991) 409–419, doi:10. 1016/0363-5023(91)90 0 06-W. [30] D.H. Gates, L.S. Walters, J. Cowley, J.M. Wilken, L. Resnik, Range of motion requirements for upper-limb activities of daily living, Am. J. Occup. Ther. 70 (1) (2016) 70, doi:10.5014/ajot.2016.015487. [31] C. Gosselin, J. Angeles, A global performance index for the kinematic optimization of robotic manipulators, J. Mech. Des. 113 (3) (1991) 220, doi:10. 1115/1.2912772. [32] J.C. Perry, J. Rosen, S. Burns, Upper-limb powered exoskeleton design, IEEE/ASME Trans. Mechatron. 12 (4) (2007) 408–417, doi:10.1109/TMECH.2007. 901934. [33] L.J. Puglisi, R.J. Saltaren, H.A. Moreno, P.F. Cárdenas, C. Garcia, R. Aracil, Dimensional synthesis of a spherical parallel manipulator based on the evaluation of global performance indexes, Rob. Auton. Syst. 60 (8) (2012) 1037–1045, doi:10.1016/j.robot.2012.05.013. [34] X.J. Liu, Z.L. Jin, F. Gao, Optimum design of 3-DOF spherical parallel manipulators with respect to the conditioning and stiffness indices, Mech. Mach. Theory. 35 (9) (20 0 0) 1257–1267, doi:10.1016/S0 094-114X(99)0 0 072-5.

600

J.-A. Leal-Naranjo, J.-A. Soria-Alcaraz and C.-R. Torres-San Miguel et al. / Mechanism and Machine Theory 140 (2019) 586–600

[35] J.J. Durillo, A.J. Nebro, JMetal: a Java framework for multi-objective optimization, Adv. Eng. Softw. 42 (10) (2011) 760–771, doi:10.1016/j.advengsoft. 2011.05.014. [36] K. Lwin, R. Qu, G. Kendall, A learning-guided multi-objective evolutionary algorithm for constrained portfolio optimization, Appl. Soft Comput. J. 24 (2014) 757–772, doi:10.1016/j.asoc.2014.08.026. [37] S. Shaphiro, M.B. Wilk, An analysis of variance test for normality, Biometrika 52 (3/4) (1965) 591–611. [38] F.J. Anscombe, The validity of comparative experiments, J. R. Stat. Soc. 111 (3) (1948) 181, doi:10.1002/jae.l. [39] J.W. Tukey, Comparing individual means in the analysis of variance, Biometrics 5 (2) (1949) 99, doi:10.2307/3001913.