# A Hybrid Strategy of Differential Evolution and Modified Particle Swarm Optimization for Numerical Solution of a Parallel Manipulator

Authors: Bingyan Mao; Zhijiang Xie; Yongbo Wang; Heikki Handroos; Huapeng Wu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents a hybrid strategy combined with a differential evolution (DE) algorithm and a modified particle swarm optimization (PSO), denominated as DEMPSO, to solve the nonlinear model of the forward kinematics. The proposed DEMPSO takes the best advantage of the convergence rate of MPSO and the global optimization of DE. A comparison study between the DEMPSO and the other optimization algorithms such as the DE algorithm, PSO algorithm, and MPSO algorithm is performed to obtain the numerical solution of the forward kinematics of a 3-RPS parallel manipulator. The forward kinematic model of the 3-RPS parallel manipulator has been developed and it is essentially a nonlinear algebraic equation which is dependent on the structure of the mechanism. A constraint equation based on the assembly relationship is utilized to express the position and orientation of the manipulator. Five configurations with different positions and orientations are used as an example to illustrate the effectiveness of the proposed DEMPSO for solving the kinematic problem of parallel manipulators. And the comparison study results of DEMPSO and the other optimization algorithms also show that DEMPSO can provide a better performance regarding the convergence rate and global searching properties.


---

## Body

## 1. Introduction

Differential evolution (DE) is a heuristic and straightforward strategy with the prominent features of global optimization and only a few control parameters [1, 2]. Particle swarm optimization (PSO) is a group theory based algorithm that was inspired by fish schools, bird flocks, and others. The disadvantage of PSO is that the individuals are easily influenced by the best particle and best position so it may only get the local optimum [3, 4]. The control parameters of PSO can be modified according to the specific optimization problems and applications [5–7]. This paper presents a new method integrated with differential evolution (DE) and modified particle swarm optimization (PSO). In particular, this strategy aims to combine the advantages of DE and the modified PSO together and then apply the hybrid algorithm to the numerical solution of a parallel manipulator.Parallel robots have been a hot research topic for many years due to their superior performance such as high response speed, high stiffness, high payload/weight ratio, low inertia, and good dynamic performance [8–10]. This paper focuses on a spatial parallel manipulator with 3 DOFs which was developed by Lee and Shah [11]. The 3-RPS parallel manipulator has three legs and each branch is a serial kinematic chain [12]. To obtain the position and orientation of the moving platform, it is necessary to define the forward kinematics of the parallel manipulator based on the length of the branches [13]. It should be noted that forward kinematics of parallel manipulators is more complicated than that of serial ones, and vice versa [14–16].Analysis of kinematics can be divided into two approaches: analytical methods and numerical methods [17, 18]. In numerical methods, the forward kinematic solution is found by solving a nonlinear global optimization problem to find the numerical solution [19]. Many numerical calculation methods for solving nonlinear equations have been utilized in kinematic problems. Newton’s iteration method, together with its improvements, is a common method that is very efficient as regards convergence speed. However, Newton’s iteration method is an arduous procedure which is sensitive to initial values and involves a large number of calculation steps [8]. Compared to PSO, the differential evolution (DE) algorithm has received much more attention due to its capability of global optimization [20, 21].The nonlinearity and multiple solution properties of the forward kinematics of a parallel manipulator make the analytical method more difficult and sophisticated than the numerical one. The numerical solution of the forward kinematics can be obtained using optimization methods such as DE and PSO or other hybrid optimization methods [22, 23]. In this paper, some modern intelligent optimization methods will be utilized and compared with each other. A time-saving hybrid strategy combined by differential evolution and modified particle swarm optimization is developed for the numerical solution of the forward kinematics of a 3-DOF parallel manipulator.


## 2. A Hybrid Strategy of Differential Evolution and Particle Swarm Optimization

The differential evolution (DE) optimization algorithm is a simulation of the biological evolution process. It is capable of handling nondifferentiable, nonlinear, and multimodal objective functions. To start an optimization process, an initial population must be randomly generated within a predefined bound, and then a new population in the next generation is generated through mutation, crossover, and selection operations.Particle swarm optimization (PSO) is a computational algorithm that optimizes a problem by iteratively improving a candidate solution with regard to a given measure requirement. The movements of the particles mimic the movement of organisms in a bird flock or fish school and are guided by their own best known position in the search space as well as the entire swarm’s best known position.Modified particle swarm optimization (MPSO) has the user-defined constant parametersω, c1, and c2 which have a great impact on the optimization performance. The parameter ω is utilized to adjust the velocity, the parameter c1 is utilized to adjust pbest(t), which is the best position achieved so far by every individual, and the parameter c2 is utilized to adjust gbest(t), which is the best value obtained so far by any particle in the population. In order to improve the performance of PSO, time-varying acceleration coefficients and time-varying inertia weight can be utilized. The inertia weight can provide a balance between the local search and global search during the optimization process:(1)ω=ω1-ω2×Itermax-IterItermax+ω2,where ω1 and ω2 are the initial and the final inertia weight, respectively, Iter is the current iteration number, and Itermax is the maximum number of iterations.By increasing the number of iterations, the weights ofpbest(t) and gbest(t) can have different convergence rates:(2)c1=c1min-c1max×IterItermax+c1maxc2=c2max-c2min×IterItermax+c2min,where cmax and cmin are the maximum and minimum acceleration coefficient of c1 and c2, respectively.The hybrid strategy of DE and MPSO (DEMPSO) is the new strategy. A key merit of DE algorithm is the efficient global optimization. Furthermore, the diversity of the entire population can be always maintained during the whole evolution process, which can prevent the individuals from falling into a local optimum. PSO, on the other hand, has the advantage of fast convergence speed. The individual with the best history and the best individual among the entire iteration are saved to obtain the lowest fitness values.Combining the advantages of DE and PSO, a new hybrid DEMPSO strategy is proposed which aims to achieve both fast convergence speed and efficient global optimization. Since the population of PSO easily falls into a local optimum, the proposed DEMPSO method uses the DE algorithm to reduce the search space first, and then the obtained populations are taken over by the MPSO as an initial population to get a fast convergence rate to the final global optimum. The hybrid algorithm can obtain the global minimum value based on the fitness function, which is built for the numerical solution for forward kinematics of a parallel manipulator. The procedure of the DEMPSO algorithm is illustrated as follows.(  1) Population Initialization. The individual x, with the population number NP, is randomly generated to form an initial population in a D-dimensional space. All the individuals should be generated within the bounds of the solution space. The initial individuals are generated randomly in the range of the search space. And the associated velocities of all particles in the population are generated randomly in the D-dimension space. Therefore, the initial individuals and the initial velocity can be expressed as follows:(3)Xi0=xi,10,xi,20,…,xi,D0xi,j0=xmin+randi,j0,1×xmax-xminVi0=vi,10,vi,20,…,vi,Np0,where 1≤i≤Np, 1≤j≤D, [xmin,xmax] is the range of the search space, and rand(0,1) is a random number chosen between 0 and 1.(  2) Iteration Loop of DE. Let individual Xi(t) denote the mutation operation at time t. By randomly choosing three individuals from the previous population, the mutant individual Vi(t+1) can be generated as the following equation:(4)Vit+1=Xr1t+FXr2t-Xr3t,where F is a differential weight between 0 and 1. It is a zoom factor to control the amplification of the mutation operation. r1, r2, and r3 are random integers that have been selected from 1 to NP and r1, r2, r3, and i are not the same as each other.The crossover operation aims to construct a new populationui,j(t+1) which is chosen from the current individuals and mutant individuals in order to increase the diversity of the generated individuals:(5)Uit+1=ui,1t+1,ui,2t+1,…,ui,Dt+1ui,jt+1=vi,jt+1,ifrand0,1≤Crorj=jrandxi,jt,otherwise,where rand(0,1) is a random number chosen from 0 to 1, jrand is an integer chosen from 1 to D randomly, and Cr is a crossover parameter that is randomly chosen from 0 to 1.In the selection operation, the crossover vectorUi(t+1) is compared to the target vector Xi(t) by evaluating the fitness function value based on a greedy criterion, and the vector with a smaller fitness value is selected as the next generation vector:(6)Xit+1=Uit+1,iffXit≥fUit+1Xit,otherwise.Update the global best part with the minimum fitness value(gbest) and the personal-best part (pbest). Let the gbest value be F1 and perform a comparison with the stopping tolerance value E1. If F1 is less than E1, the iteration of DE has finished. All the population and positions will continue to the next loop of MPSO.(  3) Iteration Loop of MPSO. Set the time-varying parameters c1, c2, and ω as the MPSO defined. Renew the individual velocity as follows:(7)vi,jt+1=w·vi,jt+c1·rand0,1·pbesti,jt-xi,jt+c2·rand0,1·gbesti,jt-xi,jt,where the constant parameters ω, c1, and c2 are defined and rand(0,1) is a number randomly chosen between 0 and 1.The new individuals are generated as follows:(8)xi,jt+1=xi,jt+vi,jt+1.Let thegbest value of MPSO be F2 and perform a comparison with fitness value E2. If F2 is less than E2, the iteration of MPSO has finished.The initial population is generated by the DE algorithm, and the stopping criterion of DE is set as the fitness valueF1 less than a user-defined stopping tolerance value E1 which is dependent on the specific kinematics of a parallel manipulator. When the fitness value F1 is less than E1, the DE population will be taken over by the MPSO algorithm. The new velocity and new location of the population are updated in every generation until the fitness value F2 becomes less than the bound condition E2 satisfied.In order to depict clearly the population moving process, the position and difference vector distribution of the population for the DE-based algorithm and MPSO-based algorithm are shown in Figures1–4.Figure 1
Position of the population of DE. (a) At the beginning of the simulation; (b) on convergence.
(a) (b)Figure 2
Difference vector distribution of DE. (a) At the beginning of the simulation; (b) on convergence.
(a) (b)Figure 3
Position of the population of MPSO. (a) At the beginning of the simulation; (b) on convergence.
(a) (b)Figure 4
Difference vector distribution of MPSO. (a) At the beginning of the simulation; (b) on convergence.
(a) (b)The initial and final distributions of the DE-based population are shown in Figures1(a) and 1(b), respectively. The position of individuals of the DE-based algorithm has converged to 10−4. The initial and final distributions of the MPSO-based population are shown in Figures 3(a) and 3(b), respectively. The position of individuals of the MPSO-based algorithm has converged to 10−3.The individuals of the MPSO-based algorithm are not gathered in a small space like the individuals of the DE-based algorithm. Based on the fitness value, the best individual of the MPSO-based algorithm will be chosen at every generation, and the algorithm will converge when the best individual meets the fitness value set before.The initial and final difference vector distributions of the DE-based algorithm are shown in Figures2(a) and 2(b), respectively. The initial difference vector distribution of the MPSO-based algorithm shown in Figure 4(a) is similar to that of the DE-based algorithm. Comparison of the final difference vector of the MPSO-based algorithm, shown in Figure 4(b), with that of the DE-based algorithm, shown in Figure 2(b), indicates that the population of the MPSO-based algorithm stops the iteration without convergence of all the individuals, while the individuals of the DE-based algorithm converge to a small space as the global optimization.


## 3. Forward Kinematics of a Parallel Manipulator

A parallel manipulator with 3 DOFs was proposed by Fang and Huang [12], which has been widely used in airplane simulators, walking machines, and others. The 3-RPS parallel robot is composed of a base platform, three legs, and a moving platform, as shown in Figure 5. Each leg is a serial chain consisting of a revolute (r) joint, a prismatic (p) joint, and a spherical (s) joint. The manipulator has three degrees of freedom: two rotations are about the x-axis and y-axis and one translation is along the z-axis. The p joints are driven by linear actuators so that the moving platform can achieve the required position and orientation.Figure 5
A 3-RPS parallel manipulator.The global coordinateO-xyz is built at the center of the base platform. The orientation of the x-axis is from point O to point A1, and the orientation of the y-axis is parallel to the line A2A3→. The r joints are evenly distributed around the base platform as an equilateral triangle. The moving coordinate O1-x1y1z1 is built at the center of the moving platform. The orientation of the x1-axis is from point O1 to point B1, and the orientation of the y1-axis is parallel to the line B2B3→. The s joints are distributed as r joints. For simplicity and without loss of generality, all coordination systems abide by the right-hand rule. The geometric parameters of the manipulator are the radius ra of the base platform and the radius rb of the moving platform.LetR and P denote the rotation matrix and position vector which move from the global coordinate system to the moving coordinates, respectively. Then, the typical kinematic chain can be denoted as a mathematical formula:(9)AiBi→=RO1Bi→+P-OAi→i=1,2,3,where AiBi→ is a vector from point Ai to point Bi. The vectors O1Bi→ and OAi→ belong to the moving coordinates and the global coordinates, respectively. R is Z-Y-X Euler rotation transformation matrix and P is position transformation vector.(10)R=cγcβsγcβ-sβcγsβsα-sγcαsγsβsα+cγcαcβsαcγsβcα+sγsαsγsβcα-cγsαcβcα,P=xeyeze,where sα=sin⁡α and cα=cos⁡α.The position vectors of pointsAi and Bi with respect to the global coordinate system and the moving coordinates can be expressed as follows:(11)OAi→=racos⁡ϕi,sin⁡ϕi,0TO1Bi→=rbcos⁡ϕi,sin⁡ϕi,0Tϕi=2π3i-1i=1,2,3.A constraint equation, based on the assembly relationship where the revolute joint axis is perpendicular to the prismatic joint axis, can be written as follows:(12)RO1Bi→+PT·ei=0ei=racos⁡φi,sin⁡φi,0Tφi=ϕi+π2,where ei is the axis of the revolute joint and its coordinate values are shown in Table 1.Table 1
Coordinate values of each rotation axis.
e 1 e 2 e 3 x 0 - 3 / 2 3 / 2 y 1 - 1 / 2 - 1 / 2The length of each linkAiBi→, that is, the inverse kinematics of the manipulator, can be calculated as a closed chain as follows:(13)di=RO1Bi→+P-OAi→TRO1Bi→+P-OAi→i=1,2,3.Given a set of lengths of the prismatic joints, the forward kinematics is to get the position and orientation of the moving platform. The numerical forward kinematic solutions of the 3-RPS parallel manipulator are a nonlinear algebraic equation. There are three rotations and three translations in the transformation matrix, but only two rotations and one translation are left in the forward kinematics since there are three extra constraints equations. Based on the specific structure of the parallel manipulator, the constraint equations can be obtained as follows:(14)xe=fα,β,γ=rb2cos⁡γcos⁡β-sin⁡γsin⁡βsin⁡α-cos⁡γcos⁡αye=fβ,γ=-rb·sin⁡γcos⁡βγ=fα,β=arctan⁡sin⁡αsin⁡βcos⁡α+cos⁡β.Substituting (14) into (13), the inverse kinematics of (13) can be simplified as(15)di=fiα,β,zei=1,2,3.For the parallel manipulator, the analytical inverse kinematic solution of the above equation is straightforward. The nominal leg lengthli can be obtained easily based on the assumed real pose αr, βr, and zer. However, given li, the analytical solution for the simulated pose αs, βs, and zes becomes very complicated and it may have multiple solutions. The given li is obtained based on the real pose αr, βr, and zer. So, the forward kinematics based on the given li is to let the simulated pose αs, βs, and zes infinitely approach the real pose αr, βr, and zer. The simulated pose converges to the real pose with the convergence of the algorithm to the global minimum. The numerical solution of this problem can be solved by minimizing the difference between the given length li and the predicted length di calculated from (15) to find a set of optimized αs, βs, and zes. The fitness function based on the least squares method can be written as(16)fss=min⁡∑i=13li-fiαs,βs,zes2,where li is a known leg length which can be acquired from the linear actuator and fi(αs,βs,zes) is used to calculate the predicted leg length during the optimization process.


## 4. Case Simulation and Results

For the manipulator studied in this work, the task of the simulation is to obtain the end-effector pose when the length of the legs is already known. In the simulation, the geometric parameters of the manipulator arera = 274 mm and rb = 158 mm. The workspace of the moving platform of the manipulator is given in Table 2.Table 2
Workspace of the moving platform.
Direction ofx-axis Direction ofy-axis Direction ofz-axis Range of translation (mm) 0 0 [ 400,700 ] Range of rotation (deg.) [ - 20,20 ] [ - 20,20 ] 0For a specific pose selected in the workspace, for instance,αr = 5°, βr = 12°, and zer = 517 mm, through an inverse calculated as l1 = 499.0178 mm, l2 = 557.7314 mm, and l3 = 534.3032 mm. On the contrary, if these three-leg lengths have been obtained from linear actuators, then the numerical optimization task is to search for an optimal combination of αs, βs, and zes to minimize the fitness function. In our first simulation, the DE-, PSO-, and MPSO-based optimization algorithms will be employed to solve this numerical problem. The aim is to investigate the performance of each algorithm and find a possible solution with the characteristics of fast convergence rate and global optimization. The simulations were implemented using Matlab R2012b on a computer with an Intel® Core™ i7-4510U CPU @ 2.00 GHz and 7.71 GB RAM. During simulation, the classical DE algorithm DE/rand/bin with the control parameters of F = 0.85 and Cr = 1 is chosen; the PSO control parameters are chosen as c1=2, c2=2, and ω=0.9. For the control parameters of MPSO, c1 was set to gradually decrease from 2.5 to 0.5 and c2 was set to gradually increase from 0.5 to 2.5. The inertia weight factor ω was set in two ways: (1) the inverse way (here, we called it MPSO1, where the control parameter ω was set to gradually decrease from 0.9 to 0.6) and (2) the direct way (here, we called it MPSO2, where the control parameter ω was set to gradually increase from 0.6 to 0.9).Using the above parameters and given the same population number of 30 and the same stopping criterion of1e-08, the simulation results of the DE-, PSO-, and MPSO-based algorithms with different iterations and computation times are listed in Table 3. From the simulation results, it can be seen that the computation time of the DE-based algorithm is less than of other algorithms, but it takes more iterations to converge; on the other hand, the MPSO1 algorithm has the fastest convergence rate with only 560 iterations, but the computation time is a little greater than in the DE algorithm.Table 3
Simulation results of forward kinematics with DE, PSO, and MPSO.
Pop. number Fitness value Iteration Time (sec) DE 30 1 e - 08 3120 0.859791 PSO 30 1 e - 08 1137 1.608181 MPSO1 30 1 e - 08 560 1.119108 MPSO2 30 1 e - 08 2433 3.561454Figure6 illustrates the logarithm-based fitness function values of DE, PSO, and MPSO with respect to the generations. It can be seen that the convergence rate of the DE-based algorithm declines gradually. The fitness value of PSO and MPSO1 drops very fast at the beginning, but the decline becomes gentler after a number of iterations. The fitness value of MPSO2 has not got any changes before 2000 generations but it suddenly converged after that period; this is a very interesting characteristic which would be utilized. In our proposed hybrid method, that is, DEMPSO, the DE algorithm part is employed to bypass the steady-state part of MPSO2, and the MPSO2 part is used to obtain a fast convergence rate.Figure 6
Logarithm of fitness function values with DE, PSO, and MPSO.For the proposed DEMPSO algorithm, besides the control parameter selection, one important issue is to decide the break fitness valueE1 for the DE algorithm since this value will influence the total computation time after MPSO2 takes over the optimization process. Table 4 presents the DEMPSO simulation results of the forward kinematics with different break fitness values of DE. For different break fitness values E1 of DE and the same stopping fitness value E2 of MPSO2, the function ran 30 times to get average results of the total iterations of DE and MPSO2, the computation time of DE at break point E1, and the total computation time of DEMPSO at terminal point E2. From the table, it can be seen that the optimization will spend the least amount of time when the fitness value E1 is equal to 70.Table 4
Simulation results of forward kinematics with DEMPSO for different fitness values of DE.
E 1 E 2 Iteration of DE Iteration of MPSO2 Break time of DE (sec) Total time of DEMPSO (sec) 80 1 e - 08 225 1030 0.017000 0.182367 75 1 e - 08 224 77 0.032000 0.099333 70 1 e - 08 235 81 0.031000 0.096767 60 1 e - 08 246 84 0.047000 0.098933 50 1 e - 08 281 62 0.050000 0.104567 10 1 e - 08 400 68 0.060000 0.104100 1 1 e - 08 688 42 0.060000 0.127433 0.5 1 e - 08 796 42 0.080000 0.142100By selecting the fitness valueE1 as a break point of DEMPSO and DEPSO, the comparison results of forward kinematics with DEMPSO, DEPSO, DE, PSO, and MPSO are listed in Table 5. The logarithm-based fitness function values with respect to iterations are plotted in Figure 7. It can be clearly seen that DEMPSO has great advantage not only as regards convergence speed but also as regards the number of iterations for solving the forward kinematic problem of the parallel manipulator. From Figure 7, it is also shown that the proposed DEMPSO has successfully bypassed the steady state of MPSO2 and inherited its steepest descent properties.Table 5
Comparison results of forward kinematics with DEMPSO, DEPSO, DE, and PSO.
Pop. number Fitness error Iteration Time (sec) DEMPSO 30 1 e - 08 267 0.076401 DEPSO 30 1 e - 08 2628 2.244032 DE 30 1 e - 08 3120 0.859791 PSO 30 1 e - 08 1137 1.608181 MPSO1 30 1 e - 08 560 1.119108 MPSO2 30 1 e - 08 2433 3.561454Figure 7
Logarithm-based fitness values with DE, PSO, MPSO, DEPSO, and DEMPSO.To validate the effectiveness of the proposed DEMPSO algorithm, in our second simulation, we randomly select five poses in the workspace and calculate the relevant leg lengths through inverse kinematics. DEMPSO was employed to search for the optimum pose of the 3-RPS moving platform for the given leg lengths. Table6 shows the final results for the five different pose situations of the 3-RPS parallel manipulator. From the table, it can be seen that the searched pose value has approached the real pose value with a very small error value (≈1e-5) when the termination fitness value E2 is set as 1e-8. The computation time is almost the same for different pose situations.Table 6
Simulation results for different situations with DEMPSO.
Number Length of leg (mm) Real pose DEMPSO E r r o r ∗ 10 - 4 Fitness error Time (sec) 1 L 1: 499.018 α r: 5° α s: 4.9999° 0.4142 7.6010 e - 08 0.081441 L 2: 557.731 β r: 12° β s: 12.0000° 0.8340 L 3: 534.303 Z e  r: 517 mm Z e s: 517.0000 0.2991  2 L 1: 489.714 α r: −3° α s: −3.0000° 0.0303 3.4300 e - 08 0.135116 L 2: 506.550 β r: 6° β s: 6.0000° 0.1344 L 3: 520.611 Z e  r: 492 mm Z e s: 491.9999 1.0628  3 L 1: 701.251 α r: 13° α s: 13.0000° - 0.1310 6.2684 e - 08 0.087767 L 2: 656.309 β r: −18° β s: −18.0000° 0.09136 L 3: 600.050 Z e  r: 641 mm Z e s: 640.9999 - 1.4288  4 L 1: 613.344 α r: −12° α s: −12.0000° - 0.2467 5.7662 e - 08 0.099704 L 2: 558.392 β r: −7° β s: −7.0000° 0.4117 L 3: 613.051 Z e  r: 583 mm Z e s: 583.0001 1.0263  5 L 1: 655.862 α r: 16° α s: 16.0000° 0.2104 8.8969 e - 08 0.116486 L 2: 767.441 β r: 19° β s: 19.0001° 0.8855 L 3: 694.392 Z e r: 695 mm Z e s: 694.9999 - 0.8913


## 5. Conclusions

In this paper, a hybrid strategy combined with DE and MPSO, termed DEMPSO, is developed to get the numerical solution of parallel manipulator forward kinematics. The proposed hybrid method benefits from the efficient global optimization of DE and the fast convergence rate of MPSO; meanwhile, it avoids the possible local optimization of MPSO. Using this method, the search bounds can be narrowed by the DE-based algorithm subtly; afterwards, the MPSO, with its fast rate of convergence, can obtain the global optimization in this narrowed search space.To validate the proposed hybrid optimization method, a 3-RPS parallel manipulator was used as an example to numerically solve the nonlinear forward kinematics. In order to get the position and orientation of the moving platform of the manipulator, the algorithms are utilized to solve the constructed fitness function. The comparison study of the proposed DEMPSO algorithm and the DE, PSO, and MPSO algorithms showed that DEMPSO can obtain the best performance for the numerical solution of parallel manipulator forward kinematics. Regarding the general case of nondifferentiable, nonlinear, and multimodal objective functions, further benchmark simulation should be carried out to verify the universality of the proposed hybrid DEMPSO algorithm.



---

*Source: 9815469-2018-02-22.xml*