Optimal fuzzy inverse dynamics control of a parallelogram mechanism based on a new multi-objective PSO

: This work presents a multi-objective optimization method based on high exploration particle swarm optimization, called MOHEPSO, for optimization problems with multiple objectives. In order to convert the single-objective (HEPSO) algorithm to the multi-objective one, its fundamentals should be changed. The leaders’ selection in the proposed algorithm is based on the neighborhood radius concept for the global best position and the Sigma method for the personal best position. Also, a fuzzy elimination technique is used for pruning the archive. The numerical results of the MOHEPSO algorithm on mathematical test functions are compared with those of other multi-objective optimization algorithms for the performance evaluation of the algorithm. Finally, the proposed algorithm is implemented to find the optimum values of controller coefficients for a parallelogram five-bar linkage mechanism. The introduced control strategy is designed based on the inverse dynamics concepts, improved by fuzzy systems and optimized by regarding two objective functions. The simulation results are presented to demonstrate the efficiency and accuracy of this approach.


PUBLIC INTEREST STATEMENT
The Particle Swarm Optimization (PSO) algorithm, originally was introduced by Kennedy (an American social psychologist) and Eberhart (an American electrical engineer), and is one of the modern heuristic algorithms. It was developed through simulation of simplified social systems as a stylized representation of the movement of organisms in a bird flock or fish school, and is robust in solving nonlinear optimization problems. In comparison with other evolutionary methods, the PSO technique can generate a high-quality solution with shorter calculation time and a more stable convergence characteristic. On the other hand, the relative simplicity of PSO, and the fact that it is a population-based technique have made it a natural candidate to be extended for multiobjective optimization. Multi-objective optimization which is also called multi-criteria optimization or vector optimization has been defined as finding a vector of decision variables satisfying constraints to give acceptable values to all objective functions.

Introduction
In engineering, the concept of optimization is noticed when complex problems are going to be solved. Such complexity may be associated with the kind of problem which an operator wants to solve (i.e. to be nonlinear or not) or the kind of solution which a person wishes to address (i.e. whether it is an exact or an approximate solution). Furthermore, the optimization problems fall into two categories: the single-objective and the multi-objective. The latter has been gaining increasing attention in the stochastic optimization community. A great variety of multi-objective optimizers have been developed, and their performance has been tested on many problems with different characteristics. In fact, the goal is to find a series of non-dominated solutions that represent a trade-off among the objectives, because unlike single-objective optimization problems, there is no single solution available for these problems (Donoso & Fabregat, 2016;Hughes, 2008).
In the practical multi-objective problems, the design criteria called objective functions may conflict with each other so that improving one of them will deteriorate another. When solving multiobjective problems, due to the presence of two or more conflicting objectives, it is often impossible to obtain a solution vector that simultaneously optimizes all of the objective functions and satisfies all of the constrains. In such cases, the concept of Pareto-optimal solutions or equivalently Pareto front arises. A state of solutions is named Pareto front if and only if there is no alternative state that would make an objective function better off without making anyone worse off (non-dominated solutions) (Cheng, Chen, Wai, & Wang, 2014;Liu, Chen, Deb, & Goodman, 2016;Meza, Espitia, Montenegro, Giménez, & González-Crespo, 2017;Razmi, Jafarian, & Amin, 2016).
Among different multi-objective optimization algorithms, multi-objective evolutionary algorithms (MOEAs), which make use of the strategy of the population evolution, are some effective methods for solving multi-objective problems. In recent years, several approaches such as multi-objective particle swarm optimization algorithms (MOPSO) (Mahmoodabadi, Bagheri, Mostaghim, & Bisheban, 2011), vortex multi-objective particle swarm optimization (MOVPSO) (Meza et al., 2017), multi-objective artificial bee colony algorithms (MOABC) (Atashkari, NarimanZadeh, Ghavimi, Mahmoodabadi, & Aghaienezhad, 2011) and etc. have been proposed for the performance improvement of MOEAs. Particle swarm optimization (PSO), first introduced by Kennedy and Eberhart, is one of the modern heuristic algorithms and was inspired by natural flocking and swarming behavior of birds and fish (Bratton & Kennedy, 2007;Valle, Venayagamoorthy, Mohagheghi, Hernandez, & Harley, 2008). Further, it was developed through simulation of simplified social systems, and has shown robust performance, in solving nonlinear optimization problems. This method is able to quickly produce a proper solution which also has better convergence properties in comparison with the other evolutionary approaches (Ding, 2017).
Nonetheless, solving multi-objective problems by PSO is involved with two matters. The first one is to sustain an archive in order to make a trade-off between convergence and diversity. The nondominated solutions are sustained by an external elitist archive which also can update itself so as to preserve the diversity. There are several methods used to update the archive. For example, many multi-objective algorithms adopt the crowding distance (Al Moubayed, Petrovski, & McCall, 2014) to prune the archive; the clustering mechanism for maintaining an archive is also applied to the multiobjective algorithms (Abbasian, Nezamabadi-pour, & Amoozegar, 2015;Padhye, Branke, & Mostaghim, 2009) to keep the size of external archive constant. The other matter in MOPSOs is updating the global best (gbest) and the personal best (pbest) positions for each particle, since there is no absolute best solution, but rather a set of non-dominated solutions. In recent years, several methods such as the Ranking methods implemented to find global best positions (Wang & Yang, 2009); the decomposition approach employed to select global and personal best positions for each particle (Gong, Cai, Chen, & Ma, 2014), etc. are proposed to determine the gbest and pbest.
In order to have an effective multi-objective algorithm with the good population diversity, an external archive strategy is often used. In fact, the external elitist archive is implemented to store the non-dominated solutions obtained by the algorithm and these solutions are filtered by a certain quality measure, such as fuzzy elimination technique, and etc. In this work, a new multi-objective particle swarm optimization algorithm based on the high exploration particle swarm optimization (HEPSO) which is one of the evolutionary methods introduced by one of the authors (Mahmoodabadi, Salahshoor Mottaghi et al., 2014) is proposed. The HEPSO seems particularly suitable for multi-objective optimization, mainly because of convergence speed, global optimality, solution accuracy, and algorithm reliability, in comparison with well-known and recent evolutionary algorithms for singleobjective optimization. It is in fact the combination of PSO with two different operators, so as to increase the exploration capability of the PSO algorithm. The first operator was inspired by the multi-crossover mechanism of the genetic algorithm (Chang, 2007), and the other uses the bee colony mechanism to update the position of the particles (Akay & Karaboga, 2012).
The introduced MOHEPSO algorithm is tested on the mathematical benchmark functions and challenged by optimally design of fuzzy inverse dynamics controller for a parallelogram mechanism. The objective functions which often conflict with each other are appropriately defined as weighting normalized summation of angle errors and weighting normalized summation of control efforts. The constant parameters of the introduced controller are regarded as the design variables used for the optimization process. The obtained Pareto front is shown and three optimum design points are selected for computer simulation and verifying the effectiveness and robustness of the proposed strategy.
The inverse dynamics controller has a long history in the control engineering and is accepted in a lot of real applications due to the efficient performance. Peng, Lin, and Su (2009) applied a computed torque control-based composite nonlinear feedback controller for robot manipulators with bounded torques. Zeinali and Notash (2010) introduced fuzzy logic-based inverse dynamic modelling for robot manipulators. Wang (2011) implemented adaptive inverse dynamics for free-floating space manipulators. Chen, Mei, Ma, Lin, and Gao (2014) proposed robust adaptive inverse dynamics control for uncertain robot manipulator. Giusti, Malzahn, Tsagarakis, and Althoff (2017) presented a combined inverse-dynamics/passivity-based controller for robots with elastic joints.
The rest of this paper is organized as follows: Section 2 gives a brief review on the multi-objective optimization. High exploration particle swarm optimization is presented in Section 3. In Section 4, the proposed algorithm is introduced. Experimental results and comparison studies to verify the capability of proposed algorithm are shown in Section 5. The description of a parallelogram controllable mechanism is presented in Section 6. Its fuzzy inverse dynamics controller is given in Section 7. The optimal fuzzy inverse dynamics controller of the parallelogram robot is simulated in Section 8. Finally, Section 9 concludes the paper.

Multi-objective optimization
From a classical standpoint, optimization of a single function simply entails determining a set of stationary points, identifying a local maximum or minimum, and possibly finding the global optimum. In contrast, the process of determining a solution for a multi-objective optimization problem is slightly more complex and less definite than a single-objective problem.
Multi-objective optimization problems (MOPs) are problems involving more than one objective to be optimized. In these problems, there are several objectives or cost functions (a vector of objectives) which are to be optimized (minimized or maximized) simultaneously. Unlike single objective optimization algorithms, the performance of multi-objective optimization solutions cannot be improved without sacrificing the performance of at least another one. Therefore, there is not a single optimal solution as the best solution with respect to all the objective functions (Kukkonen & Coello, 2017). Instead, there is a set of optimal solutions, known as Pareto-optimal solutions or Pareto front (Khishtandar & Zandieh, 2017) for such problems. The Pareto-optimal set is defined based on Pareto dominance. Each new addition to this set is compared with all the objective function values of potential solution points in order to determine if the new point is dominated. If it is non-dominated, then it is kept in the set of potential solution points. In this section, the main concepts related to the multi-objective optimization problems are described.
A general multi-objective optimization problem is defined as follows: where N is the number of objective functions involving N ≥ 2, k is the number of inequality constraints, and p is the number of equality constraints. f(x) is a N-dimensional vector of objective functions.
Definition 1 (dominance): A decision vector x 1 dominates another vector x 2 (denotes as x 1 < x 2 ) if: Definition 2 (Pareto-optimal): A vector of decision variables x* ∊ s ⊂ R n (s is the feasible region) is Pareto-optimal if it is non-dominated with respect to s.

High exploration particle swarm optimization
In nature, birds seek food by considering their personal experience and the knowledge of the other birds in the flock. This idea was first introduced by Kennedy and Eberhart to propose the PSO method. The original version of PSO suffers from trapping in local minima and premature convergence. Several approaches, such as FIPSO (Wang, Wang, Yan, & Shen, 2017), APSO (Zhan, Zhang, Li, & Chung, 2009), ACOR-PSO (Huang, Huang, Chang, Yeh, & Tsai, 2013), HEPSO (Mahmoodabadi, Salahshoor Mottaghi et al., 2014) and etc., have been thus far proposed to modify the performance of PSO. The HEPSO algorithm has shown more successes in comparison with well-known and recent evolutionary algorithms for single-objective optimization.
In the PSO algorithm, each candidate solution to the problem is considered a particle and the set of solutions is called "the swarm". Moreover, every candidate solution is associated with a velocity, and the position of each particle in population is changed according to its own experience and the other particles experience (velocity) in the flock. HEPSO is a combination of PSO with two different operators so as to increase the exploration capability of the PSO algorithm. The first one is based on the multi-crossover mechanism of the genetic algorithm (Chang, 2007), and the second one performs similarly to the bee colony mechanism (Akay & Karaboga, 2012).These formulae are used to improve the converging process and to escape from local minima in the original version of the PSO algorithm.
In PSO, the position of each particle is changed according to its own experience and that of its neighbors. Let ⃗ x i (t) denote the position of particle i, at time step t. The position of particle i is then changed by adding a velocity ⃗ v i (t) to its current position, i.e.: and the velocity vector changes in the following way: where ⃗ r 1 and ⃗ r 2 are vectors which contain the random elements from the interval [0,1]; C1 is the cognitive learning factor and represents the attraction of a particle toward its own success, and C 2 is the social learning factor and represents the attraction of a particle toward the success of the entire swarm; and finally, w is the inertia weight which is used to balance local and global searches.
In HEPSO, the first added operator uses the global best position of the population (⃗ x gbest i ) as the premier parent, and the personal best position (⃗ x pbest i ) as the second parent. It will randomly generate the new velocity for the selected particle ⃗ x i (t) as below: where ⃗ x i (t) denotes the position of particle i at iteration (time step) t, ⃗ r ∈ [0, 1] is a random vector which contains values, c2 is the social learning factor, ⃗ x gbest is the global best position, and ⃗ x pbest i is the personal best position. Furthermore, the second operator is inspired by the foraging behavior of honey bees (Karaboga, 2005), and the food source obtaining operator of Artificial Bee Colony (ABC) algorithm is used for the selected particles in the HEPSO method. The position of the randomly selected particle x i (t) will change in the dimensions as below: This is an important issue, because, the archive has to be updated where d is a random integer in the range [1, dimension], ⃗ r ∈ [0, 1] is a random value, and j is also a random integer in the range [1, number of particles]. After calculation of Equation (7), the superior value between ⃗ x i (t + 1) and ⃗ x i (t) should be selected. The flowchart of this algorithm is shown in

Multi-objective high exploration particle swam optimization
The proposed MOHEPSO algorithm is designed, as briefly described with regard to such items as diversity and convergence here. Initially, a random population is generated and in each iteration, the learning factors (C 1 and C 2 ) and the inertia weight (w) would be allocated. Also, after calculation of the fitness values of all particles, the non-dominated solutions are determined and the archive is formed.
Now, it is possible to identify the ⃗ x gbest and ⃗ x pbest from the Pareto front based on the considered strategies mentioned in the following sections. The random values 1 and 2 ∊ [0,1] are allocated to every particle. If, for a particle, ρ 1 is not greater than the standard deviation of the fitness values; or if ρ 2 is not smaller than ; the operator in Equation (7) will generate a new particle. Note that, p B denotes the probability of the bee colony, t refers to the present iteration, and max iteration denotes the biggest possible iteration number. For the particles that are not selected for the bee colony operation, another random number [0, 1] should be assigned. If a particle has < P c , then the multi-crossover operator would generate a velocity for it Equation (6). Here, P c is the multi-crossover probability. Furthermore, the particles that are not chosen for these two operations will be enhanced by PSO, and this cycle should be repeated until the user-defined stopping criterion is satisfied.
In the following subsections, the strategies to restrain the archive size and select ⃗ x pbest and ⃗ x gbest are discussed.

Archive size and fuzzy elimination
In most of multi-objective optimization methods, the archive must contain a set of specified criteria solutions, while maintaining a good spread of solutions in the obtained Pareto front. However, if all of non-dominated solutions are maintained in the archive, then the size of the archive grows very quickly.
This is an important issue, because, the archive has to be updated in each iteration and with a growing archive size, the update may become computationally demanding. Herein, a fuzzy elimination technique (Mahmoodabadi, Taherkhorsandi et al., 2014) is used to prune the archive. It involves a fuzzy parameter ɛ Fuzzy , which sets the extent of fuzzy desired in a problem, In this approach all the particles in the archive have a neighborhood radius equal to ɛ Fuzzy , and if the Euclidean distances of particles (in the objective function space) from a certain particle are fewer than ɛ Fuzzy , then they will be simply eliminated from the objective function space. Figure 2 illustrates this technique as an example in a bi-objective space. The technique is formulated as:  Table 1 and Figure 3 respectively. The E variable , as the inference outcome of the consequent variables, can be evaluated by applying any of the simplified inference, the product-sum-gravity, and the min-max gravity approaches (Rao & Kumar, 2017).
The input of the membership function ought to be evaluated as follows: where fix t 7 is the rounding function of t 7 to the closest integer toward zero.
In fact, Equation (9) normalizes the input of the membership functions for the fuzzy variables. Equation (8), however, makes the archive retain more non-dominated solutions at first iterations due to low elimination radius, and consequently increases the algorithm convergence. An increase in the iteration number causes a rise in the elimination radius and omits more close solutions. Therefore, the algorithm spread becomes broader and the non-dominated solutions are more uniformly diversified.

The strategy for the globel best position (gbest)
To update the position of particles in MOPs, one should employ the whole swarm's best particle (⃗ x gbest ) as a leader. Every particle can only has a single leader, which must be suitably selected so as  to raise the diversity and the convergence of the solutions. To this end, the selection of ⃗ x gbest is done through assigning a neighborhood radius (R neighborhood ) to all of the non-dominated solutions with regard to density measures . Accordingly, when the Euclidean distance between two non-dominated solutions is smaller than this radius, they are said to be the neighbor of together. Hence, after counting the neighbors of every single non-dominated solution at every iteration, the one with fewer neighbors is considered to be the leader.

The strategy for the personal best position (pbest)
By selecting a proper technique to find ⃗ x pbest i for particle i, the diversity within the swarm is maintained. In the proposed algorithm, each particle of population selects a personal best position from the archive during each iteration via the Sigma method (Mostaghim & Teich, 2003). These selections are distinct from each other; in the first step, the value σ i is calculated for particle i in the archive; and in the second step, σ j for particle j of the population is assigned. Then, the distance between σ i and σ j is evaluated for all the particles of both the population and the archive. Finally, the particle k of the archive is selected as the personal best position of the particle j of the population if the distance between their sigmas is the minimum. In other words, each particle with a closer sigma value (in the case of two dimensional objective space, closer means having less difference of sigma values) to the sigma value of an archive member must select that archive member as its best local guide (⃗ x pbest ).
For a two-objective space, the sigma parameter is defined as follows: where f 1 and f 2 are the first and second objective functions. Figure 4 shows the strategy for selecting the personal best position among the archive members for each particle, based on the Sigma method, in a bi-objective space.

Numerical experiments on multi-objective high exploration particle swarm optimization
This section contains the computational results obtained by the MOHEPSO algorithm using mathematical test functions compared to the other multi-objective methods.

Test functions
The test functions for multi-objective optimization should impose sufficient difficulty to challenge the algorithm in searching for the Pareto-optimal solutions. Here, SCH, FON, ZDT1, and ZDT6 benchmark problems, all of them have two objective functions, are utilized to examine the performance of the proposed method (Mahmoodabadi, Taherkhorsandi et al., 2014). Table 2 shows the formulae of the functions along with the number of dimensions, the admissible bound of their variables, the Pareto-optimal solutions, and the nature of the Pareto-optimal front for each problem.

Comparison metrics
In order to provide a quantitative assessment of the performance of multi-objective optimization on the test functions, two goals are often taken into consideration: a good convergence toward the Pareto-optimal set; and the maintenance of diversity in obtained solutions existing in the Paretooptimal set. Since both of them are conflicting in nature, comparing two sets of trade-off solutions require different performance measures. Herein, we consider five performance metrics for evaluating each of the above two goals in the set of Pareto-optimal solutions.
(1) The first is the diversity metric (∆). This measures the extent of spread achieve among the nondominated solutions (Mahmoodabadi, Taherkhorsandi et al., 2014). It is given as: Nonconvex, nonuniformly and spaced where the parameters d f and d l are the Euclidean distances between the extreme solutions and the boundary solutions of the non-dominated set; d i is the Euclidean distance between consecutive solutions in the obtained non-dominated set; d denotes the average of these distances; and finally, n is the number of members in the set of non-dominated solutions. A good distribution would make all distances d i equal to d , and would make: d f = d l = 0. Thus, for the most extensive spread set of non-dominated solutions, Δ = 0.
(2) The metric of generational distance (GD) (Azzouz, Bechikh, & Said, 2017) gives a good indication of the distance between the true Pareto-optimal front and non-dominated solution members as follows: where n is the number of members in the set of non-dominated solutions archive and d i is the minimum Euclidean distance between the member i in the non-dominated set and the Pareto-optimal front. If all members of non-dominated set are in the true Pareto-optimal front, then GD = 0.
(3) The metric of the maximum spread (MS), calculates the normalized Euclidean distance of boundary solutions in the objective space. It is formulated for a two-objective problem as: where n is the numbers in the discovered Pareto front; f i m is the mth objective of member i; and F max m and F min m are respectively the maximum and minimum of the mth objective function in the Paretooptimal front. The greater the maximum spread is, the more area of Pareto-optimal front is covered by the non-dominated solutions. If the above metric equals 1, a widely spread set of solution is obtained (Khishtandar & Zandieh, 2017).
(4) The spacing metric (SP) aims at assessing the spread distribution of vectors throughout the set of non-dominated solutions. It is calculated with a relative distance measure between consecutive solutions in the obtained non-dominated set (Mirjalili, Saremi, Mirjalili, & Coelho, 2016): where n is the number of members in the set of non-dominated solutions and d i is the Euclidean distance between the member i and its nearest neighbor in the non-dominated set. The desired value for this metric is zero, which means that the elements of the set of non-dominated solutions are equidistantly spaced.
(5) The metric of convergence (Υ) (Pal & Bandyopadhyay, 2016) measures the proximity of the generated set of Pareto-optimal solutions to the true Pareto front. First, a set of uniformly spaced solutions of the true Pareto-optimal front in the objective space is defined; and then, for each solution obtained by the algorithm, its minimum Euclidean distance (d i ) from the chosen solutions on the Pareto front is calculated. Also, n is the number of members in the set of nondominated solutions. The average of these distances is computed, and the lower the value of Υ, the closer the generated Pareto front to the true Pareto front.

Numerical results and discussion
The MOHEPSO algorithm is tested with a population of size 150. The maximum iteration, E constant and R neighborhood are respectively set at 250, 200 and 0.1. Moreover, C 1 is linearly decreased from C 1i = 2.5 to C 1f = 0.5, while C 2 is linearly increased from C 2i = 0.5 to C 2f = 2.5 over iteration. Furthermore, the bee colony probability is set at P B = 0.02, and the multi-crossover probability is set for SCH and FUN at P c = 0.95 and for ZDT1 and ZDT6 at P c = 0.1.
The non-dominated solutions obtained, in one arbitrary run, by the new considered multi-objective algorithm for the test functions given in Table 2 are shown in Figure 5. As these figures show, good diversity and convergence have been achieved by the proposed algorithm.
In order to evaluate the performance of the proposed algorithm on the test functions via the metrics, three well-known versions of multi-objective optimization algorithms are used for comparison, as detailed in Table 3.
Each algorithm is implemented in 20 independent runs, for which the average values of the performance metrics corresponding to each test function are brought in Table 3.
In other words, this table includes the average values of the metrics obtained from the four algorithms (MOHEPSO, MOPSO Ingenious, MOEA-D, and NSGA-II) for the test functions given in Table 2. As it is specified from the table, the proposed algorithm delivers a better performance in all of the functions for the GD, MS and Υ metrics except ZDT1. The MOPSO Ingenious in turn has desirable Δ and SP values in all the test functions except ZDT6. The MOEA-D, however, performs poorly in MS and Υ metrics along with the SCH and FUN test functions. In addition, the MOHEPSO has desirable Δ values in ZDT6; and MOPSO Ingenious has a better performance in Υ metric for the ZDT1 test function (Table 4).

Dynamical equation of a parallelogram five-bar linkage mechanism
The mechanism in Figure 6 is called a five-bar linkage, where having all revolute joints. Clearly there are only four bars in the figure, but in the theory of mechanisms, it is a convention to count the ground as an additional linkage, which explains the terminology. In Figure 6, it is assumed that the lengths of links 1 and 3 are the same, and that the two lengths marked l 2 are also the same. In this way, the closed path in the figure is in fact a parallelogram, which greatly simplifies the computations. Notice, however, that the quantities l c1 , and l c3 need not to be equal. In other words, even though links 1 and 3 have the same length, they need not to have the same mass distribution (Spong & Vidyasagar, 2008).
It is clear from the figure that, even though there are four links being moved, there are in fact only two degrees of freedom, identified as q 1 andq 2 . At first, the coordinates of the centers of mass of the various links are written as a function of the generalized coordinates. This gives Next, with the aid of these expressions, the velocities of the various centers of mass are written as a function of q 1 and q 2 . For convenience, the third row of each of the following Jacobian matrices is dropped as it is always zero. The result is: x c1 y c1 = l c1 cos q 1 l c1 sin q 1 x c2 y c2 = l c2 cos q 2 l c2 sin q 2 x c3 y c3 = l c2 cos q 1 l c2 sin q 2 + l c3 cos q 1 l c3 sin q 1 x c4 y c4 = l 1 cos q 1 l 1 sin q 1 + l c4 cos(q 2 − ) l c4 sin(q 2 − ) l c3 cos q 1 l c3 sin q 1 = l 1 cos q 1 l 1 sin q 1 − l c4 cos q 2 l c4 sin q 2 (17) v c1 = −l c1 sin q 1 0 l c1 cos q 1 0 q v c2 = 0 −l c2 sin q 2 0 l c2 cos q 2 q v c3 = −l c3 sin q 1 −l 2 sin q 2 l c3 cos q 1 l 2 cos q 2 q v c3 = −l 1 sin q 1 −l 2 sin q 2 l 1 cos q 1 l 2 cos q 2 q The velocity Jacobians J v c1 , i = 1, … ,4 are defined as the four matrices appearing in the last equations. Next, it is clear that the angular velocities of the four links are simply given by Thus the inertia matrix is given by If Equation (17) is substituted into Equation (19) and the standard trigonometric identities are used then  IF then the inertia matrix is diagonal and constant, and as a consequence the dynamical equations will contain neither Coriolis nor centrifugal terms.
Turning now to the potential energy:

Hence
Notice that ∅ 1 depends only on q 1 but not on q 2 and similarly that ∅ 2 depends only on q 2 but not on q 1 . Hence, if Equation (21) is satisfied, then the rather complex looking manipulator in Figure 6 is described by the decoupled set of equations as follows.
where u 1 and u 2 are the control inputs. Equation (24) can be rewritten as follows.
The state-space formulation of the system can be described as Equation (26).

Fuzzy inverse dynamics controller for the mechanism
The inverse dynamic model of an arbitrary manipulator can be expressed as: where u = [u 1 , u 2 , … , u n ] T is the vector of input generalized forces, and n denotes the number of generalized coordinates of the manipulator. Where inertia matrix M(q) is an n × n symmetric positive definite matrix, C(q,q) is an n × n matrix of centripetal and Coriolis terms and g(q) is an n × 1 vector of gravitational terms; q = [q 1 , q 2 , … , q n ] T , q = [q 1 ,q 2 , … ,q n ] T and q = [q 1 ,q 2 , … ,q n ] T are the displacement, velocity and acceleration vectors of joints respectively. In Equation (27), if (20) d 11 (q) = m 1 l 2 c1 + m 3 l 2 c3 + m 4 l 2 1 + I 1 + I 3 d 12 (q) = d 21 (q) = m 3 l 2 l c3 − m 4 l 1 l c4 cos(q 2 − q 1 ) d 22 (q) = m 2 l 2 c2 + m 3 l 2 2 + m 4 l 2 c4 + I 2 + I 4 (21) m 3 l 2 l c3 = m 4 l 1 l c4 y ci = g sin q 1 (m 1 l c1 + m 3 l c3 + m 4 l 1 ) = g sin q 2 (m 2 l c2 + m 3 l 2 − m 4 l c4 ) � 1 = g cos q 1 (m 1 l c1 + m 3 l c3 + m 4 l 1 ) � 2 = g cos q 2 (m 2 l c2 + m 3 l 2 − m 4 l c4 ) then Finally, the tracking error e(t) = q d − q satisfies where K p and K d are diagonal matrices with diagonal elements consisting of position and velocity gains, respectively. The following choice for the gain matrices K p and K d causes an exponentially stability of the system.
The fuzzy inverse dynamics controller with input of the membership function e(t) = q d − q and output F for the robot manipulator could be calculated using Equations (30) and (31) for each link as the following equations.
where Then, in Equation (29): The membership functions are shown in Figures 7 and 8, and their rules are mentioned in Tables 5 and 6.

Optimal fuzzy inverse dynamics controller of the parallelogram mechanism
In this section, the proposed MOHEPSO algorithm introduced in Section 4 with a population of size, maximum iteration, E constant and R neighborhood 250, 300, 200 and 0.2 is used to determine the appropriate parameters of the fuzzy inverse dynamics controller for the five-bar linkage mechanism. In this paper, weighting normalized summation of angle errors and weighting normalized summation of control effort of the both links are regarded as the objective functions and calculated using Equation (34). q = a q , a q = −K p q − K dq + r and r =q d where W 1 and W 2 are the weighting coefficients. q 1 and q 2 are displacement of joints. u 1 and u 2 are the input generalized forces.
The vector [VB, PO, ZB, HS, HM, HB, M 1 , ΔM 1 , M 2 , ΔM 2 ] is the chosen parameters of fuzzy inverse dynamics controller. VB, PO, ZB, HS, HM and HB are input of the membership function. M 1 , ΔM 1 , M 2 and ΔM 2 are positive coefficients of the fuzzy inverse dynamics controller. The weighting normalized summation of angles errors and weighting normalized summation of control effort are the functions of this vector's components. This means that by selecting various values for the selective parameters, we can make changes in the weighting normalized summation of angles errors and weighting normalized summation of control effort. Indeed, this is an optimization problem with two objective functions and ten decision variables (VB, PO, ZB, HS, HM, HB, M 1 , ΔM 1 , M 2 , ΔM 2 ). The design variables corresponding to the optimal design points A, B, and C are illustrated in Table 7. The parameters of the multi-objective algorithms are chosen in Section 4 in a similar way.
In Figure 9, points A and C stand for the best weighting normalized summation of angles errors and weighting normalized summation of control effort, respectively. This figure illustrates that all the optimum design points in the Pareto front are non-dominated and can be chosen by the designer based on the design criteria. In Figure 9, point B can be regarded as a trade-off optimum choice when the minimum values of both the weighting normalized summation of angles errors and weighting normalized summation of control effort are considered. The real tracking trajectories of the optimal design points A, B, and C are shown in Figure 10. The tracking errors of the optimal design points A, B, and C are shown in Figure 11. Moreover, Figure 12 illustrates the control effort of the parallelogram mechanism for the optimal design points A, B and C.

Conclusions
This paper has proposed a new and simple generic evolutionary multi-objective optimization algorithm based on HEPSO which is called MOHEPSO. Different properties of the proposed MOHEPSO method have made it robust and competitive among the other existing methods in the literature. In order to supply a proper infrastructure for comparison, a fuzzy elimination technique is utilized to prune the archive instead of previous techniques in the literature. Also, the Sigma method is used to find the personal best positions of particles. Moreover, a neighborhood relationship among all the non-dominated solutions based on the Euclidean distance is defined to find and update the global best positions. Several experiments have been conducted to compare the performance of the proposed method with other recent and well-known multi-objective evolutionary algorithms. As the experiments demonstrate, the obtained optimal Pareto fronts in the proposed algorithm (MOHEPSO) outperform three effective multi-objective optimization algorithms (i.e. MOEA/D, MOPSO Ingenious and NSGA-II) for different complex multi-objective test functions in terms of both solution diversity and convergence. Finally, the proposed algorithm is used for the optimum design of a fuzzy inverse dynamics controller for a parallelogram five-bar linkage mechanism. Two conflicting objective functions, the weighting normalized summation of angle errors and weighting normalized summation of control effort, are regarded to design the optimal controller. Three points of Pareto front, obtained from MOHEPSO, are selected to work out the ten parameters of the optimal fuzzy inverse dynamic control; hence, the designer has the ample opportunity to select the finest point based on the design criteria. The result elucidated that the first selected point had minimum normalized summation of angles errors and maximum normalized summation of control effort. However, the third point was the direct opposite of the first one. Indeed, it had maximum normalized summation of angles errors and minimum normalized summation of control effort. Thus, the designer can regard the second well-chosen point as the trade-off optimum choice.