Empirical Study of Constraint-Handling Techniques in the Optimal Synthesis of Mechanisms for Rehabilitation

Currently, rehabilitation systems with closed kinematic chain mechanisms are low-cost alternatives for treatment and health care. In designing these systems, the dimensional synthesis is commonly stated as a constrained optimization problem to achieve repetitive rehabilitation movements, and metaheuristic algorithms for constrained problems are promising methods for searching solutions in the complex search space. The Constraint Handling Techniques (CHTs) in metaheuristic algorithms have different capacities to explore and exploit the search space. However, the study of the relationship in the CHT performance of the mechanism dimensional synthesis for rehabilitation systems has not been addressed, resulting in an important gap in the literature of such problems. In this paper, we present a comparative empirical study to investigate the influence of four CHTs (penalty function, feasibility rules, stochastic-ranking, and ϵ-constraint) on the performance of ten representative algorithms that have been reported in the literature for solving mechanism synthesis for rehabilitation (four-bar linkage, eight-bar linkage, and cam-linkage mechanisms). The study involves analysis of the overall performance, six performance metrics, and evaluation of the obtained mechanism. This identified that feasibility rules usually led to efficient optimization for most analyzed algorithms and presented more consistency of the obtained results in these kinds of problems.


Introduction
The increasing number of people with motor deficiencies has been a crucial factor in developing rehabilitation devices. Motor injuries frequently occur in the upper and lower extremities due to degenerative diseases, neuronal disorders, sports injuries, traffic, and work accidents. Presently, rehabilitation devices for upper and lower limbs have been developed [1,2]. Stationary rehabilitation devices are used more frequently in the rehabilitation of the patient. However, these devices are generally expensive because the control of several degrees of freedom is required to generate a rehabilitation routine. Rehabilitation systems based on closed-chain mechanisms are a low-cost alternative for health care [3]. These systems can provide a pre-established rehabilitation routine that controls a degree of freedom. Therefore, these devices can be considered for rehabilitation at home.
In the design of closed-chain mechanisms, a synthesis process of the mechanism is carried out [4]. In this process, the length of links is determined to achieve the desired movements for rehabilitation. The graphic, analytic, and optimal methods [5] are the used techniques for the synthesis of closed-chain mechanisms. Graphic synthesis of mechanisms is performed by a visual representation of the mechanism, where Euclidean geometry concepts are used to obtain the design solutions. This method is characterized by providing quick solutions with a low precision rate.
The analytical synthesis is based on getting algebraic expressions that describe the movement of the mechanism. This method offers solutions with a high precision rate. However, the number of precision points is limited by the number of independent variables of the closed-chain mechanism because it results in an overdetermined system of equations [6] requiring a numeric solution to solve them. The optimal synthesis of mechanisms is performed through the optimization problem solution. This synthesis method supports multiple-precision points, but the precision rate can decrease caused by the numerical optimization technique used. Currently, the last method is frequently used for the synthesis of mechanisms.
In the optimal synthesis of closed-chain mechanisms for rehabilitation, indirect and direct search methods (numerical optimization techniques) [7] have been used to find suitable solutions. In the indirect search methods, derivatives of the objective function and constraints are required. For instance, Ávila-Hernández and Cuenca-Jiménez [8] designed a thumb prosthesis using a nine-bar spatial mechanism. For the solution of the synthesis problem, the function "findminimum" of Mathematica ® was used. Zhiming Ji and Yazan Manna [9] designed a lower limb rehabilitation system using a four-bar linkage. In the syntheses process, the method "lsqnonlin" of MATLAB ® was implemented.
Wang et al. [10] designed a rehabilitation system using a four-bar linkage mechanism in the synthesis process. In that work, the search method "fmincon" of MATLAB ® was implemented. Tsuge et al. [11] designed a rehabilitation system using a ten-bar linkage mechanism. The synthesis problem was solved by Bertini™ software.
Tsuge et al. [12] designed a rehabilitation system using a Stephenson III linkage mechanism. The synthesis problem was solved using algorithms, such as Fletcher-Reeves, Newton, conjugate gradient, interior point, and Quasi-Newton.
On the other hand, direct search methods do not require such derivatives, and the objective function is considered in its original form to find a solution. The use of metaheuristic algorithms as the direct search method is gaining more attention in the solution of mechanism synthesis problems [13] because they do not depend on the problem features (nonlinear, discontinuous, etc.), they can endow different operators to find the most promising solution, and they are easy to implement. A summary of various studies found in the literature that use metaheuristic algorithms in the synthesis problem of mechanisms is presented in Table 1.
The particular interest in this study is in the works related to the mechanism design for rehabilitation presented at the end of such table. For instance, Bataller et al. [14] proposed a mechanism for finger rehabilitation, where the direct search algorithm Málaga University Mechanism Synthesis Algorithm (MUMSA) was used, and a penalty function was implemented as a constraint-handling technique in the synthesis process.
Shao et al. [15] designed a rehabilitation system using a cam-linkage mechanism for the lower limb rehabilitation to perform the gait. The synthesis problem was solved by a genetic algorithm, where a penalty function was considered as the constraint-handling technique.
Calva-Yáñez et al. [16] designed a gait rehabilitation system for children using a fourbar linkage mechanism. The indirect search algorithm Sequential Quadratic Programming (SQP) and the direct one given by Differential Evolution (DE) algorithm were used in the synthesis process. In this work, feasibility rules were considered as a constraint-handling technique.
In this work, mechanical constraints in the synthesis process are not considered. Leal-Naranjo et al. [18] presented the design of a wrist prosthesis using a spherical parallel manipulator. In the synthesis process, Nondominated Sorting Genetic Algorithm (NSGA-II), Multiobjective Particle Swarm Optimization (MOPSO), and Multiobjective Evolutionary Algorithm based on Decomposition (MOEA/D) algorithms were considered, and the feasibility rules were implemented as constraint-handling technique.
Muñoz-Reina et al. [3] designed a lower limb rehabilitation machine. The differential evolution algorithm was considered in the synthesis process, and the feasibility rules were used as the constraint-handling technique. Table 1. Summary of investigations tackling metaheuristic algorithms in the kinematic synthesis problem. "-" indicates that the use of CHT is not required.

Study
Application of the Metaheuristic Algorithms CHT Synthesis Problem [19] Four-bar mechanism Genetic Algorithm - [20] Four-bar mechanism Genetic Algorithm -Stephenson's six-bar mechanism Watt's six-bar mechanism [5] Four-bar mechanism Genetic Algorithm Penalty Function [21] Hand robot mechanism Pareto Optimum Evolutionary Feasibility rules Multiobjective Algorithm (POEMA) [22] Four-bar mechanism Differential Evolution Penality Function [23] Six-bar mechanism Differential Evolution Penality Function [24] Four-bar mechanism Differential Evolution Penality Function [25] Four-bar mechanism Genetic algorithm-fuzzy logic Penality Function [26] Four-bar and six-bar mechanisms MUMSA Penality Function [27] Four-bar mechanism Genetic Algorithm, Penality Function Differential Evolution, Particle Swarm Optimization [28] Four-bar mechanism Ant-gradient Penality Function [29] Four-bar mechanism GA-DE Penality Function [30] Six-bar mechanism Cuckoo Search Penality Function [31] Four-bar mechanism NSGA-II Feasibility rules [32] Four-bar mechanism Imperialist competitive algorithm, Penality Function Genetic Algorithm, Differential Evolution, Particle Swarm Optimization [33] Four-bar mechanism Modified Krill Herd Penality Function [34] Four-bar mechanism TLBO Penality Function Genetic Algorithm, Particle Swarm Optimization [35] Four-bar mechanism Hybrid Lagrange Interpolation DE Penality Function (HLIDE) [36] Four-bar and six-bar mechanisms Hybridization Differential Evolution Penality Function with Generalized Reduced Gradient Mechanisms for rehabilitation [18] Spherical parallel manipulator NSGA-II, MOPSO, MOEA/D Feasibility rules in prosthetic wrist [14] Six-bar mechanism in MUMSA Feasibility rules finger rehabilitation [15] Cam-linkage mechanism Genetic Algorithm Penalty Function. in gait rehabilitation [16] Four-bar mechanism Differential Evolution Feasibility rules in gait rehabilitation [17] Four-bar mechanism in Particle Swarm Optimization and TLBO gait rehabilitation and orthotic devices [3] Eight-bar mechanism in Differential Evolution Feasibility rules. lower limb rehabilitation We observed in the literature that metaheuristic algorithms in their original versions were proposed to solve unconstrained optimization problems. Thus, it is important to find the most promising techniques for including real-word constraints [37][38][39]. Currently, for the solution of constrained optimization problems, CHTs have been proposed. According to [40], these can be classified into two generations, where the penalty functions, decoders, special operators, and separation of the objective function and constraints are included in the first generation. On the other hand, the feasibility rules, stochastic ranking, and -constraint involve the second generation. The first generation is characterized by a high calibration complexity and a high probability of convergence to local minimums. The second generation is characterized by overcoming the mentioned limitations in the first generation. From the literature review regarding the design of closed-chain mechanisms for rehabilitation and the general framework of kinematic synthesis of mechanisms, both included in Table 1, the metaheuristic algorithms incorporate only one constraint-handling technique to solve the constraint optimization problem. Figure 1 shows the use of CHTs reported in the literature for the synthesis of mechanisms using metaheuristic algorithms. As observed in that figure, 64% of the reviewed literature uses the penalty functions, followed by the feasibility rules with a 24% of use.
According to [37,39,41], the constraint-handling technique is an essential factor to solve constrained optimization problems because these techniques influence the exploration and exploitation capabilities of algorithms for the search of the most promising solution. However, other CHTs in the second generation, such as stochastic ranking and -constraint have not been evaluated in solving these types of problems. 12% 24% 64% Unconstrained Feasibility Rules Penalty function Figure 1. Use of the constraint-handling techniques in the optimal synthesis of mechanisms with metaheuristic algorithms according to Table 1.
Depending on the problem nature in the kinematic synthesis for rehabilitation (multimodality, types of functions, design variable space, and the feasible space ratio), the most frequently used algorithms and the constraint-handling techniques can be suitable during the search process of the metaheuristic algorithms. Nevertheless, the high interactions in the problem nature and the metaheuristic algorithms make it difficult to suggest a constraint-handling technique that is more likely to perform well against different algorithms.
Motivated by these facts, and to the best of the author's knowledge, there is a gap in the comparative study of the performance of different constraint-handling techniques in the dimensional synthesis of mechanisms. This research can aid researchers or practitioners interested in applying metaheuristic algorithms in the kinematic synthesis of this kind of problem with insights about the quality and consistency of the obtained results in the studied CHTs and can identify the most promising CHT alternatives for the synthesis problems.
Thus, researchers and practitioners can have some basic knowledge in the analyzed CHTs about the likelihood of increasing the optimization efficiency and the consistency of the obtained results when those CHTs are included in different metaheuristic algorithms. Hence, this work can suggest the use of a particular CHT in the optimizers to achieve optimum kinematic synthesis.
Therefore, in this work, the performance of the penalty function, feasibility rules, stochastic-ranking, and -constraint is studied to solve three real-world engineering problems reported in the literature related to the synthesis of mechanisms for rehabilitation by using metaheuristic algorithms. In this study, ten metaheuristic algorithms commonly used in the synthesis of mechanisms for rehabilitation are implemented in such CHTs. These are eight variants of a differential evolution algorithm, a Genetic Algorithm (GA), and a Particle Swarm Optimization (PSO) algorithm.
The parameters of the algorithms are tuned using the irace package to make a fair comparison and provide useful insights about the studied CHTs. On the other hand, the four-bar linkage mechanism, the cam-linkage mechanism, and the eight-bar linkage mechanism are considered in the mechanism synthesis for rehabilitation. This paper is organized as follows: In Section 2, the general problem of mechanism synthesis for rehabilitation is presented. Section 3 describes the metaheuristic algorithms and the constraint-handling techniques used in this study. Section 4 states three state-of-art synthesis problems for rehabilitation to be solved. In Section 5, the study of the constrainthandling techniques is presented and discussed. In addition, the best and worst solution are analyzed and compared. Finally, Section 6 presents our conclusions.

Optimization Problem Statement
In this work, optimization problems for the synthesis of rehabilitation mechanisms are studied. The synthesis problem usually requires satisfying more than one design objective and constraint. Therefore, the design problems can be considered as a Constrained Multiobjective Optimization Problem (CMOPs). According to literature related to the synthesis of mechanisms for rehabilitation presented in Table 1, it is observed that CMOPs are typically transformed into constrained mono-objective optimization problems by the weighted sum approach [42,43]. Hence, mono-objective optimization problems for the mechanism synthesis are only considered in this study. These are formally stated as follows: subject to: x min ≤ x ≤ x max (4) whereJ(x) = ∑ n i=1 w i J i (x) is the weighted objective function; J i is the i − th design objective of the system; x is the design variable vector; w i is the i − th weight attributed to each objective; n is the number of design objectives; g j (x) and h k (x) are the inequality and equality constraints, respectively; and x min , x max are the design variable bounds.

Constraint-Handling Techniques in Metaheuristic Algorithms
Metaheuristic algorithms are optimization strategies used to solve a wide range of hard optimization problems. The metaheuristics are inspired by biological systems, social behaviors, among others.
In recent years, metaheuristic algorithms have been used to provide satisfactory solutions to optimization problems related to the mechanism synthesis, see Table 1. The behavior of metaheuristic algorithms to search for better solutions can be improved through suitable constraint handling. In this work, the behavior of four Constraint-Handling Techniques in different metaheuristic algorithms is studied. The CHTs [44] involve the Penalty Functions (PF), Feasibility Rules (FR), Stochastic Ranking (SR), and -Constraint ( C) method.
The considered metaheuristic algorithms are eight Differential Evolution variants [45], the Genetic Algorithm [46], and the Particle Swarm Optimization [47]. Those techniques and algorithms are frequently used in mechanism synthesis for rehabilitation and have provided satisfactory solutions, as is shown in Table 1. However, according to that table, differential evolution is one of the most promising reported optimizers. Therefore, the main DE variants [48] are analyzed in this study.

Differential Evolution Algorithm
The DE algorithm was proposed by Storn and Price [45], and it is a direct and stochastic search method. DE was inspired by the theory of natural selection and is shown in the Algorithm 1.
Algorithm 1 Differential evolution pseudocode 1: Generate an initial population X 0 with NP individuals. 2: Evaluate X 0 . 3: Initialize the best individual x best . 4: G ← 0 5: while G ≤ G max do 6: for all x i ∈ X G do 7: Generate a child individual u i based on (5)-(12). 8: Evaluate u i and U G ← u i . 9: end for 10: Select the new population X G+1 between X G and U G according to CHT.
[X G+1 , x best swarm ] ← f ncCHT(X G , U G , x best ) (This function is associated to the CHT.) 11: G ← G + 1 12: end while DE starts with a population of individuals called parents X 0 = {x 1 , x 2 , ..., x NP }, where NP is the number of individuals in the population. At each generation G, the individuals in the population X G perform a crossover and mutation process to generate a population of child individuals U G . The crossover and mutation process depend on the used DE variant. In this case, eight DE variants are chosen, and these are: DE/rand/1/bin (DER1B), DE/rand/1/exp (DER1E), DE/best/1/bin (DEB1B), DE/best/1/exp (DEB1E), DE/current to rand (DECR), DE/current to best (DECB), DE/current to rand/1/bin (DECR1B), and DE/current to rand/1/exp (DECR1E).
The crossover and mutation process for each variant is defined by Equations (5)- (12), where CR is the crossover factor, F and K are scale factors, r 1 , r 2 , and r 3 are three random individuals from the current population, D is the number of design variables and j rand ∈ {1, 2, . . . , D} is a random value that represents the crossover point. Binomial recombination: Rand/1/Bin: Best/1/Bin: Exponential recombination: Rand/1/Exp: Best/1/Exp: Arithmetic recombination: Current-to-Rand/1: Current-to-Best/1: Arithmetic-binomial recombination: Current-to-Rand/1/Bin: Otherwise (11) Current-to-Rand/1/Exp: The population of parents X G and children U G compete at each generation to remain to a new population of individuals X G+1 for the next generation. The selection process between the parents and children is based on the CHT, as defined in Section 3.2. Then, the best individual in the population X G max is considered the best solution for the optimization problem.

Particle Swarm Optimization
The PSO algorithm was proposed by Kennedy and Eberhart [47]. This algorithm is inspired by the collaborative behavior of species in search of food and is shown in Algorithm 2.
Algorithm 2 Particle swarm optimization pseudocode 1: Initialize the swarm position u i ∈ U 0 with NP particles. 2: Evaluate the swarm U 0 . 3: Initialize the best known position X 0 = U 0 . 4: Initialize the best position of the swarm x best swarm . 5: Initialize the velocity of each particle v i ∈ V 0 . 6: G ← 0 7: while G ≤ G max do 8: Update the inertial weight w based on (13). 9: for all u i ∈ U G do 10: Update the velocity v i based on (14). 11: Update the position u i based on (15). 12: Evaluate u i . 13: end for 14: Update the best known position X G+1 and the best particle in the swarm x best swarm with the use of CHT between X G and U G . [X G+1 , x best swarm ] ← f ncCHT(X G , U G , x best swarm ) (This function is associated to the CHT.) 15: The PSO algorithm initializes with a swarm U G=0 with NP particles. Each particle is located in a position u i with a velocity v i ∀ i = {1, 2, ..., NP}. The best known positions of the particles are saved in the swarm X G=0 . The best particle in the swarm X 0 is stored in x best swarm . In the time G (iteration), the i − th position and velocity of each particle u i are updated by Equations (13)- (15), where w is the inertial weight, v max and v min are the maximum and minimum velocity values, x best is the best neighborhood particle of x i in a ring topology, and C 1 , C 2 are weighting factors.
Finally, the best particle x best swarm and the best particles positions X G are updated considering the CHT as is defined in Section 3.2. Then, the best solution is found in x best swarm .

Genetic Algorithm
The GA was proposed by John Henry Holland [49]. This algorithm is inspired by genetic-molecular evolution. The essential operation of GA is shown in Algorithm 3.

Algorithm 3
Genetic algorithm pseudocode 1: Generate a initial population X 0 with NP chromosomes. 2: Evaluate X 0 3: G ← 0 4: while G ≤ G max do 5: for all x i ∈ X G do 6: Obtain x r 1 and x r 2 in X G by tournament. 7: Generate a child v i by (16) 8: Generate a mutant u i by (17) 9: Evaluate u i 10: end for 11: Replace the population X G to X G+1 considering X G and U G in the CHT. X G+1 ← f ncCHT(X G , U G , ∼) (This function is associated to the CHT.) 12: G ← G + 1 13: end while GA initializes with a group X 0 = {x 1 , x 2 , ..., x NP } of NP chromosomes. In the group, two parent chromosomes x r 1 and x r 2 are selected by tournament [50] and recombined to generate a child chromosome v i as shown in (16). The recombination process uses a piece-wise multipoint crossover [5], where the crossover points are given by a crossover rate CR ∈ [0, 1] . The child chromosome v i can mute using uniform mutation [51] as shown in (17), where MR ∈ [0, 1] is a mutation rate, and x j min , x j max are the minimum and maximum values of the design variables, respectively.
Between the mutant chromosomes U G and the parent one X G , the best chromosomes are selected based on the replacement mechanism using the CHT defined in Section 3.2.

Constraint-Handling Techniques
Metaheuristic algorithms have been used to solve complex optimization problems. In their original versions, these have been designed to solve unconstrained optimization problems. However, optimization problems in engineering usually have constraints. Currently, constraint-handling techniques, such as the PF, FR, SR, and C method, have been the most popular techniques to solve constrained optimization problems [44]. Each one of these techniques is described below.

Penalty Function
The PF [52] has been frequently used to solve constrained optimization problems in the mechanism synthesis. The PF transforms an optimization problem with constraints into an unconstrained problem. The most used transformation equation is: whereJ (x) is the new function to minimize,J is the weighted objective function, x is the vector of design variables, m is the number of constraints, γ j (x) is the infeasible constraint distance defined by (19), and j is the penalty factor that can be static, dynamic, or adaptive.
In engineering optimization problems, the term j is usually defined with a static value.
Algorithm 4 describes the operation of PF. This constraint-handling technique requires three sets of solutions (three input arguments). At each iteration G, those solutions are compared according to the new functionJ (x) (18). The best solutions are stored in the output arguments.

Feasibility Rules
The FR were proposed by Deb [53]. Three rules define this constraint-handling technique, these are:

1.
Between two infeasible solutions, the solution with the fewest number of violated constraints is chosen.

2.
Between a feasible solution and an infeasible one, the feasible solution is chosen.

3.
Among two feasible solutions, the solution with the best objective function is chosen.
In this work, the number of violated constraints is considered the sum of infeasible constraint distance φ(x) and shown in (20).
Algorithm 5 describes the operation of FR. This constraint-handling technique requires three sets of solutions (three input arguments). At each iteration G, those solutions are compared according to the feasibility rules. Finally, the best solutions are stored in the output arguments. Input:ã ∈ X G andb ∈ U G Output: X G+1 and x best 2: X G+1 =ã and x best =c 3: for i = 1 to NP do 4: 7: x best ←b i 8: end if 9: end if 10: end for

-Constraint Method
The C method was proposed by Takahama et al. [54]. This method relaxes the constraints to select solutions close to the feasible region. The selection process is described by (21).
where is a relaxation factor defined by (22) and (23), cp is a parameter to control the speed of reducing relaxation of constraints, T c is the maximum number of iterations (generation or time) to relax such constraints, φ(x θ ) is θ-th best value of the objective function in the initial solution vector (population, swarm, or chromosomes).
The -constraint method finds feasible regions using the gradient of constraints at an infeasible region. Algorithm 6 describes the operation of the C method. This constrainthandling technique requires three sets of solutions (three input arguments). The gradientbased mutation process can improve an infeasible solution to obtain a new solution u new j . The new solution u new j and old solution are compared using (21), and the best solution is stored in the current solution b j . The set of current solutions are compared according to (21), and the best solutions are stored in the output arguments. Input:ã ∈ X G andb ∈ U G Output: X G+1 and x best 2: X G+1 =ã and x best =c 3: for j = 1 to NP do 4: if rand(0, 1) < Pg then 5: is the constraint vector, ∆C(b j ) −1 is the pseudo-inverse of the constraint's derivative obtained by Moore-Penrose pseudo-inverse using the singular value decomposition and Rg is the number of attempts to improve the solution. 8:  20: x best ← b j

21:
end if 22: end if 23 [55]. This constraint-handling technique makes the pseudo-random sort of solutions. The ranking of solutions is based on the bubble sort algorithm with a probabilistic factor P f ∈ [0, 1]. The values close to zero in P f mean a high probability that the sort is based on the constraint distance. On the contrary, this is based on the objective function.

Study Cases in the Mechanism Synthesis for Lower Limb Rehabilitation
In this work, three study cases in engineering related to the mechanism design for the rehabilitation of lower limbs are considered for testing four constraint-handling techniques in different metaheuristic algorithms. Case 1 is related to a four-bar linkage mechanism, Case 2 involves a cam-linkage mechanism, and Case 3 includes an eight-bar linkage mechanism. Each study case is described below.

Case 1: Four-Bar Linkage Mechanism
The four-bar linkage mechanism has been used in lower limb rehabilitation for adults and children with cerebral palsy [9,16]. The dimensional synthesis design problem consists of finding the link lengths that approximate the coordinates [x i P , y i P ] of the coupler link point P to the coordinates [x i P ,ȳ i P ] of the ankle trajectory in a standard gait cycle. The schematic representation of the four-bar mechanism is presented in Figure 2.
According to Figure 2, the vector of design variables comprises the lengths of the links r 1 , r 2 , r 3 , r 4 , the ground link angle θ 1 , the initial position of the mechanism [x 0 ,y 0 ] with respect to the coordinate system x − y, and the crank angles θ i 2 ∀i = {1, 2, ...n f }. The vector of design variables is defined as: The optimization problem for the dimensional synthesis of the mechanism is presented in (25)- (29), where the constant values w 1 = 1 and w 2 = 0.01 weight the design goals.

Min
subject to: The first design goal refers to the error between the desired precision points [x i P ,ȳ i P ] related to the ankle trajectory and the path [x i P , y i P ] of the mechanism's coupler link point. The second design goal refers to the angular displacement error between the desired angles θ i 2 and the crank angles θ i 2 . In this problem, the desired anglesθ i 2 are defined by (30).
The precision (desired) points [x i P ,ȳ i P ] are defined based on the standard ankle trajectory, these points are shown in Table A1 of Appendix A, and are obtained in [16]. On the other hand, according to Figure 2, the points generated by the coupler link [x i P , y i P ] are defined by (31), where θ 2 is the angle of the crank link (input link), and θ 3 is shown in (32). This is obtained by the Freudenstein equation. It is important to note that both solutions in θ 3 (related to the sign of the root) are computed, and the solution that provides a better objective function is chosen. Thus, with the selection of θ 3 , the resulting mechanism can be an open (+) and crossed (−) four-bar mechanism.
x Parameter x min x max

Case 2: Cam-Linkage Mechanism
The second optimization problem is related to the synthesis of a cam-linkage mechanism. The original design problem was presented by Yixin Shao et al. [15]. This mechanism is designed to provide rehabilitation routines for adults with motor problems in the lower limbs. The schematic diagram of the cam-linkage mechanism is shown in Figure 3. The design variable vector (33) is defined by the link lengths r 1 , r 2 , ..., r 8 , the angles β, γ, η, the initial position of the mechanism [x 0 , y 0 ] with respect to the x − y coordinate system, and the position of the slide displacement axis defined by the length e and the angle α.
subject to: The first term of the weighted objective function (34) refers to the error between the crank angles θ i 1 ∀i = {1, 2, ..., n f } and the desired anglesθ i 1 ∀i = {1, 2, ..., n f }, where n f is the number of precision points, andθ i 1 is defined by Equation (47).
On the other hand, the second term refers to the base radius R 0 of the cam (see Figure 3). Both design objectives are weighted by w 1 = 0.3 and w 2 = 0.7, and normalized by Θ re f = 0.4246 [rad] and R re f = 0.3782 [m] according to the best results presented in [15].
To guarantee the functionality of the mechanism, the design problem has inequality constraints. The first three constraints (35)-(37) guarantees a complete revolution of the crank angle θ 2 (Grashof criterion), where s min and s max are the minimum and maximum values of s(θ 1 ) in the precision points. The s(θ 1 ) is given in (48) and is visualized in Figure 3.
The next inequality constraints (38)-(41) provide high efficiency of the force transmission from the input link (crank link) to the output link, where the distances between point A to C d AC , and between point D to F d DF are obtained by (49) and (50), respectively.
The inequality constraints (42)-(44) guarantee the path feasibility in the output link, where the distances between point D to F d OD , and between point B to e d Be are obtained by (51) and (52), respectively.
The last inequality constraint (45) defines the minimum curvature radius of the cam profile. Moreover, the design variable vector x is bounded by x min and x max . The limit values used in this work are shown in Table 3. Table 3. The maximum and minimum parameters of the vector x for the cam-linkage mechanism. x x min x max x x min x max The precision points (x i P ,ȳ i P ) are presented in Table A2 of Appendix A. These coordinates define the standard ankle trajectory in the human gait. The precision points are obtained by the leg kinematics (53) The computation of the displacement angle θ i 1 and the cam radius R 0 in the optimization problem are obtained by inverse kinematics analysis of the mechanism. Therefore, R 0 is given in (54) and θ i 1 in (55).
The above equations are obtained according to the original design problem [15]. Except for the Equation (55), which is rewritten in this work from a two-argument arctangent. Equation (55) has two solutions that define the direction of crank rotation θ 1 . In this case, the clockwise turn in the crank link is chosen. Then, the solution in (55) that presents a positive increment between two different consecutive crank positions is selected.
On the other hand, the cam profile is obtained in (56) with respect to the coordinate system

Case 3: Eight-Bar Linkage Mechanism
Case 3 includes an eight-bar linkage mechanism for rehabilitation of the lower limbs. This mechanism has been designed for rehabilitation routines for the average anthropometry for the Mexican population [3]. The integration among the kinematic synthesis, the link shape design, and the mechanism dynamics are considered for its design. The schematic diagram of the mechanism is displayed in Figure 4a.
(b) Link shape. According to Figure 4a, the design variables vector x is defined by (57).
This vector is split into two kinds of parameters, the kinematic parameter vector p ki = [l 1 , l 2 , ..., l 9 , l 6 , l 8 , θ 1 , θ 9 ,θ 6 ,θ 8 The shape parameters define the octagonal link shape shown in Figure 4b and the link dynamic parameters, such as the masses, the mass centers, and the inertia of links.
The synthesis problem for the eight-bar linkage mechanism is presented in (58)-(84). The design objectives in (58) are weighted by parameters w 1 = 1 and w 2 = 1e − 6. The first term refers to the average error between the precision points (x i E ,ȳ i E ) and the points (x i E , y i E ) generated by the mechanism's coupler link. The desired precision points [x i E ,ȳ i E ] are defined as a semi-elliptical rehabilitation trajectory. Those points are shown in Table A3 of Appendix A. On the other hand, the points [x i E , y i E ] are set by the direct kinematics of the mechanism, as shown in (85). For more details of the analysis, see [3].
subject to: The second term in (58) represents the average of the applied torque τ in the crank link. The inverse dynamic modelX =Ǎ −1B is used to obtain the torque τ considering the position of the crank link known with a constant angular speed of . . , 8}, and the input torque τ. Moreover,Ǎ ∈ R 21×21 andB ∈ R 21 . For more details, see [3].
For the synthesis of the mechanism, a set of inequality constraints are established. Equations (59)-(64) define a crank-rocker mechanism configuration. The Equations (65)-(70) are related to the high efficiency of force transmission from the crank link (input link) to output link. In the coupler link, the patient's ankle is placed. Equation (71) avoids tear failure in the joint holes, where Fc min = 1.5 is the safety factor, σ adm = 6.2052 × 10 7 [Pa] is the allowable stress in the joint, σˆk s = || Fk s || a s is the joint stress, a s = π 4 φ s e s is the contact area of the joint, φ s = 0.01905 [m] is the constant diameter of the holes.
Equations (72)-(83) define the octagonal link shape of the mechanism, where ζ s = 1.5φ s is the distance between the edge of the hole and the link edge, and ψ s = 1.5φ s is the distance between holes. The final constraints are related to the limits of the design variable vector x represented by the vectors x min and x max . These values are defined in the Table 4.

Kinematic Parametersp ki
x min x max

Comparative Experimental Study
We analyzed four Constraint-Handling Techniques (CHTs) in lower limb rehabilitation mechanism synthesis. These CHTs are the Feasibility Rules (FR), the Stochastic-Ranking (SR), the -Constraint ( C) method, and the Penalty Function (PF). The constraint-handling techniques are included in ten Metaheuristic Algorithms (MAs) frequently used in the synthesis of mechanisms. These algorithms are: Particle Swarm Optimization (PSO), Genetic Algorithm (GA), and eight variants of Differential Evolution (DER1B, DER1E, DEB1B, DEB1E, DECR, DECB, DECR1B, and DECR1E).
Each metaheuristic algorithm solves three study cases of mechanism synthesis problems related to lower limb rehabilitation using four-bar linkage, cam-linkage, and eight-bar linkage mechanisms. A summary of the main attributes of the study cases is included in Table 5 because there is no single attribute to measure the problem's complexity. In such a table, the terms LI and N I mean the number of linear and nonlinear inequality constraints, respectively.
The ratioρ = |F|/|S| ∈ [0, 1] provides the relative size of the feasible region in the search space (a measure of how difficult is to generate feasible solutions) by using |S| = 1 × 10 6 random design variables (suggested in [58]), and the number of feasible solutions |F|. A high value inρ indicates that the MA would find the feasible region in an early number of generations. Based on those attributes, the less complex optimization problem is given by study case 1, followed by study case 2, and study case 3, in that order. The following comparative analysis is divided into three subsections. First, the conditions for making the parameter tuning of the CHTs in the algorithms are described in Section 5.1. In Section 5.2, the descriptive and inferential statistics of the overall performance for each algorithm with different CHTs are compared and discussed. Six performance metrics used in constrained optimization problems are included in Section 5.3 for evaluating the behavior of the CHTs in the metaheuristic algorithms. In Section 5.4, the best and worse solutions per each study case are visualized and evaluated to investigate the performance of the obtained mechanisms in the rehabilitation routine.

Parameter Tuning Conditions of the CHTs in the Algorithms
The parameter tuning for each metaheuristic algorithm with the different CHTs was performed for each study case. The irace package [59] was used for parameter tuning and to make the comparative study more fair and meaningful [60]. The parameters obtained by the irace package are presented in Tables A4-A6 of Appendix B. For solving each study case, the same number of objective function evaluations is also considered to make fair comparisons. The penalty factor k in the CHT by PF in (18) is taken into account as a constant value in all cases. Thus, these parameters are set in Table 6. 10,000 10,000 10,000

Statistical Analysis of the Overall Performance
In this work at each study case, thirty executions of the four CHTs are carried out per each of the ten metaheuristics. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). Each sample contains the best thirty objective functions of an algorithm. Descriptive statistics analyzes those samples. The descriptive statistical results are shown in Tables A7-A9 of Appendix C, and these consider the following measures: the mean, the standard deviation (std), the median, and the maximum (max) and minimum (min) values of the samples. The graphical visualization of the distribution of the samples is shown in the box plot of Figure 5.
Based on the results presented in such tables, it is possible to define the CHT that delivers the best performance in search of the most promising solution for each study case. Thus, the best behavior of a metaheuristic algorithm is related to the ability to search for the best solutions and provide more reliable results in different runs. The mean, the standard deviation (std), the median, and the maximum (max) and minimum (min) values of the samples are chosen as the statistical measure to indicate the algorithm capacity to find suitable solutions through runs in the particular samples. The best CHT between two CHTs, called A and B, included in the algorithms, is obtained by the following comparison procedure: • The CHT A can be compared to the CHT B if their CHTs are different, and the algorithms that implement such CHTs are the same. On the contrary, they cannot be compared. • The statistical measure from the thirty executions of the CHT A is better than the CHT B for the same algorithm if the former presents less value than the latter. When this happens, the CHT A obtains the point (a win) in favor. Thus, the maximum number of points for comparing each CHT is 150 (three CHT comparisons per five measures per ten algorithms). • The CHT with the highest number of points is the best constraint-handling technique.  In the comparison method, the lower values of the statistical measures (the mean, the standard deviation (std), the median, and the maximum (max) and minimum (min) values of the samples) are preferable in each algorithm. All measures are considered in the comparison procedure. The results obtained by the comparison process are presented in Table 7. The best performances of metaheuristic algorithms include the FR for the solution of the three study cases. For study case 1, it wins 115, for study case 2, it wins 104, and for study case 3, it wins 105. The next CHT alternative is the PF because it obtains the second-best performance in two study cases (study cases 1 and 2). The inferential statistical analysis was performed to make general conclusions about the performance of the CHTs incorporated in the metaheuristic algorithms for the solution of the lower limb rehabilitation mechanism synthesis. Nonparametric statistical tests [61] were used because the independence, normality, and homoscedasticity assumptions in the samples were not fulfilled due to the stochastic nature of the algorithms.
The inferential statistical analysis consists of performing a multi-comparative post-hoc analysis with Holm post-hoc error correction method [61] in all CHTs incorporated in the algorithms. The first step is to guarantee that at least one comparison of the constrainthandling technique in the metaheuristic algorithms is different. The 95%-confidence Friedman test is applied for that purpose, and the results are presented in Table A10 of Appendix D. The obtained p-value near to zero proves the existence of such differences with confidence close to 100%.
Once the first step was done, the Friedman test for multiple comparisons with the Holm post-hoc error correction method was performed to make the general conclusions of comparisons and confirm the results' confidence. The forty samples (four CHTs included in ten algorithms), conformed to each one by the best thirty objective functions of the executions, can be compared if both CHTs are incorporated in the same metaheuristic algorithm.
In Table A11 of Appendix D, the results obtained by the multi-comparative analysis are shown. The two-sided alternative hypothesis was selected, which means that, in the case that the p-value of the test is less than the 5% statistical significance, the comparison between the CHT A and the one in B for the same algorithm presents significant differences (the alternative hypothesis is accepted). On the contrary, if the results are not convincing, and significant differences are not observed in the comparative result then the null hypothesis is accepted. Once the alternative hypothesis is accepted, the Friedman z-value sign defines the winner between both CHTs in the algorithms. The negative sign of z indicates that the CHT in the algorithm A was better than the one in B, and vice versa.
Three comparisons related to the four CHTs are made per each metaheuristic algorithm. A total of thirty comparisons are carried out per each CHT at each study case (considering the ten algorithms), and the highest number of wins defines the most suitable constraint-handling technique in the synthesis of lower limb rehabilitation mechanisms for each study case. A total of ninety comparisons was carried out, considering the three CHT comparisons with the ten metaheuristic algorithm for the three study cases.
The number of wins per CHT in algorithms from the comparisons is presented in Table 8. The results indicate that the feasible rules improve the algorithm performance in the comparisons at 66.66% for study case 1, 40% for study case 2, and 53.33% for study case 3. The overall improvement for the ninety comparisons (in all study cases) with the FR is around 53.3%. It is also observed in Table A11 of Appendix D that only eight comparisons with other CHTs (DECB EC, DECR1B SR, DECR1E EC, GA SR, GA PF, PSO SR, PSO C, and PSO PF) outperformed the values obtained in FR (around 8.8% of the total comparisons).
Nevertheless, the obtained solutions in such comparisons did not surpass the best solutions provided by the FR as is observed visually in Figure 5 and numerically in Tables A7-A9 of Appendix C. Then, the improvement in the eight algorithms using the CHTs in such comparisons did not influence a better overall behavior. In addition, in the remaining ninety comparisons (around 29.44%) with respect to the FR, there were no conclusive results about improving an algorithm. The latter means that there are no other CHTs in the algorithms that present a better behavior than the FR.
On the other hand, the penalty function outperformed at 43.33% in study case 1, 53.33% in study case 2, and 33.33% in study case 3 of the comparisons with respect to the other CHTs. The -constraint method improved at 36.66%, 3.33%, and 46.66% in the corresponding study cases. The stochastic ranking outperformed at 0%, 33.33%, and 13.33% of the comparisons. Nevertheless, in the last three CHTs, several comparisons exist where the CHTs in the algorithm to which it is compared outperformed such results as shown in Table A11 of Appendix D .
Thus, it is confirmed, with 95% confidence, that the feasibility rules improved the search capabilities of the algorithm with respect to the other CHTs. This implies that the optimal synthesis of the lower limb rehabilitation mechanism can be benefited by using the FR to solve such a constrained optimization problem because this can improve the quality and the consistency of results in most of the algorithms reviewed in this work.

Statistical Analysis of the CHT Behavior through Metrics
This work is also interested in knowing the search capabilities of metaheuristic algorithms with the constraint-handling techniques related to different performance metrics. Then, six performance metrics [62] found in the specialized literature were used. Each of these metrics is described below: As commented previously, thirty executions of each constraint-handling technique included in the metaheuristic algorithms were carried out for each study case. The values of the metrics FP, P, AFES, and SP for the thirty executions for each algorithm at each study case are presented in Tables A12-A15 of Appendix E. The summaries (averages values) of those metrics are displayed in Table 9. We observed that, for the metric FP, the C provided the best results in the three study cases, followed by the FR and PF, in that order. For the metric P, the FR obtained the best results in cases 1 and 3. In the second study case, the result in FR was close to the best one provided by the PF. In the AFES metric, the FR, C, and PF won in only one study case. For the SP value, the FR provided the best result in two study cases.
In the case of the metrics EVALS and PR, inferential statistics using the non-parametric Friedman test for multiple comparisons are included to confirm the performance. The distribution of the samples related to such metric values per each algorithm at each study case is shown in the box plots of Figure 6 and Figure 7, respectively, and the numerical results of the EVALS and PR metrics are presented in Tables A16, A17, A18, A21, A22 and A23 of Appendix E, respectively.
The multi comparative Friedman test is presented in Appendixes E.6 and E.8. Summaries of the winning numbers of these comparisons are displayed in Tables 10 and 11.    Box plot of the sample distribution obtained by PR in ten metaheuristics at each study case of mechanism synthesis for lower limb rehabilitation.
The EVALS metric for study case 1 indicates that feasible solutions were found in the initialization of the 200 individuals (particles) of the population (swarm). As a result, all executions required at least 200 evaluations to find a feasible solution (the feasible region). For this reason, the samples were not significant to determine a winner in the CHT in such a metric. In study case 2, the FR confirmed the superior performance in the EVALS metric. In study case 3, the -constraint method improved the convergence to the feasible region in fewer generations. Finally, in the PR metric, it was confirmed in two study cases that the -constraint method can improve the solutions through the feasible region when the generations occur, followed by the SR that won in only one study case.
According to the results obtained by the previous analysis in the performance metrics, the following findings of the CHTs in the optimal synthesis of the lower limb rehabilitation mechanism synthesis were observed: • The FR had one of the best performance and the best one of the search for feasible solutions and satisfactory regions for the three study cases, respectively, because this constraint-handling technique has a high probability value in the FP and P metrics. The relationship between the AFES and P gave by SP is the best among the CHTs. Therefore, among the studied CHTs, the Feasibility Rules technique is the most reliable one for the synthesis of mechanisms for rehabilitation. • The PF can be considered the second option for solving mechanism synthesis for rehabilitation problems because it showed a suitable performance to find satisfactory solutions (P and SP metrics) compared to Stochastic-Ranking and -constraint. This establishes acceptable reliability to find feasible solutions (FP metric). • C method presented the best performance in the search for feasible regions as shown in FP. Inside the feasible region, this CHT produces improved solutions based on the PR metric. Nevertheless, the probability of this CHT to search for satisfactory solutions was low (P and SP metrics). This behavior is attributed to the gradient-based mutation where fast convergence to feasible local regions is promoted.
• SR is the less reliable CHT because some algorithm executions can not find feasible solutions (FP metric). This presents a low performance in the search for suitable solutions (P metric).
According to the results in this section, the Feasibility Rules presented the best performance for handling constraints in the mechanism synthesis for rehabilitation problems. The Penalty Functions gave the second best performance. Those CHTs improved the quality and the consistency in most of the algorithms reviewed in this work.

Evaluation of the Obtained Mechanisms
In this section, the best solutions obtained by the feasibility rules (according to the overall performance presented in Table 8), and the worst constraint-handling technique are presented.
In study case 1, the best constraint-handling technique was the feasibility rules, and the worst was stochastic-ranking. The Computer-Aided Design (CAD) and the generated trajectory of the best mechanisms obtained by both constraint-handling techniques are presented in Figure 8, where the objective function values of the mechanisms obtained by the feasibility rules (with the algorithm DECR) is J = 0.002098, and by stochastic-ranking (with the algorithm DECR1E) is J = 0.003726. Those values are given in Table A7 of Appendix C.
According to the design parameters obtained in both CHTs (see Table A26 of Appendix F), the feasibility rules present a Pearson correlation coefficient of 0.9939 with respect to the stochastic-ranking. This indicates a high degree of similarity (see Figure 8). On the other hand, the performance of the design goals is presented in Figure 8c,d. Figure 8c shows the trajectory of the mechanisms, the current, and the desired precision points. The average error between the desired precision points and the current ones isJ 1 Figure 8d shows the angles obtained by the crank and the desired angles, where the average angular error isJ 2 = 1.6270 [rad] by feasibility rules, andJ 2 = 3.2663 [rad] by stochastic-ranking. We observed that the design goals in both mechanisms related toJ 1 and J 2 , present a trade-off, i.e., in one goal, it is better but in the other not. Then, solutions obtained by the feasible rules and the stochastic-ranking are suitable for the rehabilitation mechanism. Nevertheless, the feasible rules provide more reliable results through different algorithm executions and with a better overall performance (related to the objective function J).   In study case 2, the feasibility rule was the best constraint-handling technique, and the worst was the epsilon-constraint. The CAD and the generated trajectory of the best mechanisms obtained by both constraint-handling techniques are presented in Figure 9. The objective function values given by the mechanisms obtained by the feasibility rules (DEB1B) is J = 0.6557 and by epsilon-constraint (DEB1E) is J = 0.7588. Those values are given in Table A8 of Appendix C.
According to the design parameters obtained in both CHTs (see Table A27 of Appendix F), the feasibility rules present a Pearson correlation coefficient of 0.9991 with respect to the epsilon-constraint. This correlation is larger than the obtained mechanisms in the previous study case 1, which implies a fairly strong relationship. In order to highlight the difference, the performance of the design goals is presented in Figure 9c,d for the cam shapes and the angles obtained by the crank, respectively.
The normalized average error between the desired crank angles and the current ones isJ 1 = 1.2217 by feasibility rules andJ 1 = 1.7122 by epsilon-constraint. The normalized base radius of the cam isJ 2 = 0.4137 by FR, andJ 2 = 0.3502 by epsilon-constraint. It is observed again that the design goals present a trade-off, such that both mechanisms are suitable for the rehabilitation routine, and it depends on the application that one selects. However, the feasible rules can search for solutions with a similar performance through different executions and with better quality in the overall performance in the weighted objective function J.  Finally, for study case 3, the best constraint-handling technique was the feasibility rules, and the worst one was stochastic-ranking. The CAD and the generated trajectory are presented in Figure 10, where the objective function values of the mechanisms obtained by the feasibility rules (DECR) is J = 0.0001903, and by stochastic-ranking (DECR) is J = 0.0001884. Those values are given in Table A9 of Appendix C.
According to the design parameters obtained in both CHTs (see Table A28 of Appendix F), the feasibility rules present a Pearson correlation coefficient of 0.7231 with respect to the stochastic-ranking. As the problem complexity in study case 3 is higher than in the other cases, the obtained solution presents differences, as is observed in Figure 10a,b. The behavior of the design goals related to the generated trajectory in the mechanisms and the torque applied in the crank link are displayed in Figure 10c and 10d, respectively.
In this case, the average error between the desired precision points and the current ones isJ 1 = 8.1216 × 10 −5 [m] by feasibility rules andJ 1 = 5.9633 × 10 −5 [m] by stochasticranking. The average torque isJ 2 = 5.4531 [Nm] by feasibility rules, andJ 2 = 6.4395 [Nm] by stochastic-ranking. Neither mechanism can be considered the best because they present a trade-off in the design goals. Thus, design solutions obtained by the feasible rules and the stochastic-ranking are suitable for the rehabilitation mechanism. Nevertheless, the feasible rule provides more reliable results through different executions and the second-best overall performance (related to the objective function J).  . CAD representation and trajectory generated of the mechanism (related to the design goals) given by the best and the worst constraint-handling technique for solving study case 3.

Conclusions
In this work, the behavior of four constraint-handling techniques for the synthesis of rehabilitation mechanisms is studied in ten metaheuristic algorithms.
According to the results, we present the following conclusions: • Through the statistical analysis, we observed that the feasibility rules present the best overall performance for the solution of the three mechanism synthesis cases under study. This is because the feasibility rules show a high probability of convergence towards feasible solutions and satisfactory regions based on the performance metrics studied in this work. • In the solution of mechanism synthesis problems for rehabilitation, there exists a high probability that the search capability of an algorithm is improved by incorporating the FR because, in most of the algorithms reviewed in this work, the FR enhanced the quality and consistency of results. • The penalty functions can be used as the second option for solving mechanism synthesis for rehabilitation problems. These establish acceptable reliability to find feasible solutions based on the performance metrics. • The -constraint method and stochastic-ranking presented the worst overall performance in the solution of the mechanism synthesis for lower limb rehabilitation. This feature is due to the mechanism synthesis problem, which presents a complex search space and wide, with highly nonlinear functions (multimodality), several design variables, and feasible regions that are difficult to find. The -constraint presents a fast convergence to feasible local regions because the gradient-based mutation increases the speed in the search of feasible solutions. • We confirmed that all study cases for the assessment of the obtained mechanisms, by using the feasibility rules and the worst CHT, presented different trade-offs in the design goals of the mechanism synthesis problem for lower limb rehabilitation. Nevertheless, FR can endow the algorithm with a better search capability to find reliable results through different executions with a high probability of finding the best overall performance.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors acknowledge support from the Secretaría de Investigación y Posgrado (SIP) and the Comisión de Operación y Fomento de Actividades Académicas (COFAA) from Instituto Politécnico Nacional (IPN). The first author acknowledges support from the Mexican Consejo Nacional de Ciencia y Tecnología (CONACyT) through a scholarship to pursue graduate studies at CIDETEC-IPN.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Precision Points for Synthesis of Mechanisms
In this section, the precision points for the synthesis of mechanisms for lower limb rehabilitation are presented for each study case.

Appendix A.1. Precision Points for the Case 1: Four-Bar Linkage Mechanism
For the synthesis of the four-bar mechanism, the standard path of the ankle of a 12-year-old child is considered. In Table A1, those data are shown and are obtained by [16].

. Precision Points for Case 2: Cam-Linkage Mechanism
For the synthesis of the cam-linkage mechanism, the data in Table A2 are considered. These data represent the ankle path in the standard gait for an adult with a thigh-length L thigh = 0.3984 [m] and a shank-length of L shank = 0.3956 [m].  Table A3 for the synthesis of this mechanism.

Appendix B. Parameter Tuning
This section shows the tuning of the parameters for the algorithms per each study case. For this process, the irace program was used, where 5000 experiments are established, and the default configuration is considered.

Appendix B.1. Parameter Tuning of the Case 1: Four-Bar Linkage Mechanism
The parameters obtained by the irace program for the four-bar mechanism are presented in Table A4.  The parameters obtained by the irace program for study case 2 related to the camlinkage mechanism are stated in Table A5.  The parameters obtained by the irace program for the eight-bar mechanism are displayed in Table A6.

Appendix C. Descriptive Statistics of the Overall Performance
This section shows the statistical data collected from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). Each sample contains the best thirty objective functions of an algorithm. The following measures are considered: the mean, the standard deviation (std), the median, and the minimum and maximum values of the sample.

Appendix C.1. Descriptive Statistics for the Case 1: Four-Bar Linkage Mechanism
The measures of the descriptive statistics for the four-bar mechanism are presented in Table A7.  The measures of the descriptive statistics for the cam-linkage mechanism are shown in Table A8.  The measures of the descriptive statistics for the eight-bar mechanism are displayed in Table A9.

Appendix D. Inferential Statistics of the Overall Performance
This section shows the inferential statistics of the data collected from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). Each sample contains the best thirty objective functions of an algorithm. For carrying out these statistics, the R program and the library scmamp are used.

Appendix D.1. Ranks Achieved by the Friedman Test
The ranks obtained by the Friedman test for the three study cases are shown in Table A10.

. Multiple Comparison Friedman Test
The collected data from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms are used to make the multiple comparison Friedman test with the post-hoc correction method. Holm correction is used in this analysis. This test is presented in Table A11. Table A11. Multiple comparison Friedman test with the post-hoc correction method for the overall performance of the four CHTs included into the ten metaheuristic algorithms. "-" indicates that there is no winning algorithm.

Appendix E. Performance Metrics of the CHT Behavior
This section shows the statistical data collected from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). In each sample, six performance metrics, FP, P, AFES, SP, EVALS, and PR, are computed. The results are shown below.

Appendix E.1. FP Metric
The average values for FP metric are shown in Table A12.

Appendix E.2. P Metric
The average values for P metric are shown in Table A13.     The descriptive statistics for the four-bar linkage mechanism are presented in Table A16. The measures of the descriptive statistics for the cam-linkage mechanism are presented in Table A17. The measures of the descriptive statistics for the eight-bar mechanism are presented in Table A18. Table A18. EVALS descriptive statistics for the study case 3: Eight-bar linkage mechanism. IS means the number of infeasible solutions through the executions. "-" indicates that it is not possible to compute such a statistic. This section shows the inferential statistics of the EVALS data collected from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). For carrying out these statistics, the R program and the library scmamp are used.

CHT
Appendix E.6.1. Ranks Achieved by the Friedman Test The ranks obtained by the Friedman test for the three study cases are shown in Table A19. The collected data from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms are used to make the multiple comparison Friedman test with the post-hoc correction method. Holm correction is used in this analysis. This test is presented in Table A20. Table A20. Multiple comparison Friedman test with the post-hoc correction method for the overall performance of the four CHTs included into the ten metaheuristic algorithms for the EVALS metric. "-" indicates that there is no winning algorithm. The measures of the descriptive statistics for the four-bar linkage mechanism are presented in Table A21. The measures of the descriptive statistics for the cam-linkage mechanism are presented in Table A22. The measures of the descriptive statistics for the eight-bar mechanism are presented in Table A23. This section shows the inferential statistics of the PR data collected from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms. Each study case groups a total of forty samples (four CHTs per ten metaheuristics). For carrying out these statistics, the R program and the library scmamp are used.
Appendix E.8.1. Ranks Achieved by the Friedman Test The ranks obtained by the Friedman test for the three study cases are shown in Table A24. The collected data from the thirty executions of the four CHTs per each one of the ten metaheuristic algorithms are used to make the multiple comparison Friedman test with the post-hoc correction method. Holm correction is used in this analysis. This test is presented in Table A25. Table A25. Multiple comparison Friedman test with the post-hoc correction method for the overall performance of the four CHTs included into the ten metaheuristic algorithms for the PR metric. "-" indicates that there is no winning algorithm.

Appendix F. Optimum Design Variables Per Each Study Case
This section presents the design variable vector x of the optimum solutions (the best solutions) of the best and the worst CHTs obtained from the thirty executions. The optimum solutions of the best and the worst CHTs are presented in Table A26. The optimum solutions of the best and the worst CHTs are presented in Table A27. The optimum solutions of the best and the worst CHTs are presented in Table A28.