Robotic arm trajectory optimization based on multiverse algorithm

1.   Introduction

Recently, with the continuous development of robotic arm technology, the requirements for its working environment are becoming increasingly stringent. Therefore, the performance requirements for robotic arms are also getting higher and higher. Under different working conditions, the operating time, speed, acceleration, jerk, and other performance requirements of the manipulator are different. Optimizing the trajectory of the robotic arm has always been the goal of designers.

Trajectory optimization can be generally classified into single-objective optimization, and multi-objective optimization ^[1] according to the optimization objective. Single-objective optimization usually optimizes only one objective without considering other performance metrics. For example, in a trajectory optimization problem concerning time, increasing the speed to reduce the running time also increases the impact on the robot arm. Multi-objective optimization needs to consider two or more performance optimalities simultaneously, obtain the corresponding Pareto solution set, and then select the optimal solution based on the working conditions.

Rabab Benotsmane et al. ^[2] elaborated a new "whip-lashing" method that aims to realize an optimized trajectory for the five-degree-of-freedom RV-2AJ robot arm, and the method lowered energy consumption and increased productivity. Han et al. ^[3] proposed a particle swarm optimization (PSO) algorithm that can dynamically adjust the learning factor. This method effectively combines the piecewise polynomial interpolation function with PSO, avoids the complex process of particle swarm algorithm to construct the adaptation function, and improves the problem that the traditional PSO is more probably to fall into local extreme value in the early stage and convergence speed is slow in the later stage. Zhao et al. presented an efficient numerical method for trajectory optimization for motion planning involving complicated robot dynamics. Peng et al. ^[4] employed the non-dominated neighborhood immune genetic algorithm, obtaining each joint's optimal position, velocity, acceleration, and jerk planning curves. Kei Ota et al.^[5] proposed a reinforcement learning-based algorithm for trajectory optimization for constrained dynamical systems. Shi et al.^[6] proposed a multi-objective genetic algorithm (NSGA-II) and multi-objective particle swarm optimization algorithm (MOPSO) in order to settle trajectory planning with kinematic and dynamic constraints.

In this study, a 3-5-3 spline polynomial and multiverse optimization algorithm are used to plan the trajectory of the manipulator in the joint space. Nevertheless, the multiverse optimization algorithm performs poorly in solving large-scale optimization problems and lacks jumping out of local poles^[7]. The ability of value makes it impossible to find the optimal global solution. This paper proposes a multi-objective optimization based on MVO and a method to improve the wormhole probability curve, aiming to improve the convergence speed and enhance the global search ability through population mutation. The final result confirms the method's feasibility, improves the manipulator's efficiency and reduces energy consumption.

2.   Trajectory planning

2.1.   Construction of joint trajectory based on 3-5-3 spline polynomial

The parameters, including speed, acceleration, and jerk of the manipulator, have better real-time performance in joint space, do not need to consider particular problems, and are more convenient to control and describe. As a result, trajectory planning is generally carried out in joint space. In order to solve the relationship between each joint variable and time, it is necessary to discretize the task trajectory performed by the manipulator first and then perform inverse kinematics analysis. In the trajectory planning of the manipulator, the 3rd-degree polynomial interpolation can only ensure that the velocity and displacement are continuous, and the acceleration may jump, which is not necessarily continuous. Moreover, the discontinuous acceleration will make the joints of the robotic arm shake and impact, which will influence the working efficiency and service life of the robotic arm^[8]. Compared with the 3rd-degree polynomial, the 5th-degree polynomial interpolation can ensure continuous and smooth curves of joint displacement, velocity, and acceleration, even though the calculation amount is increased^[9]. Generally, the higher the order of the polynomial function, the more constraints it can satisfy. Therefore, a 5th-degree polynomial is used to fit each joint position sequence node, and a continuous and smooth trajectory of displacement, velocity, and acceleration can thus be obtained^[10,11].

Assuming that the robot's starting point, ending point, and two intermediate path points are known, each joint angle of each joint at the four interpolation points can be obtained by solving the inverse kinematics equation. The interpolation angle of joint $i$ is denoted by $\theta_{ji}$, where $i$ represents the number of joints $(i = 1, 2, 3, 4, 5, 6)$, and$j$ denotes the serial number of the interpolation point $(j = 1, 2, 3, 4)$. The general formula of the 3-5-3 spline polynomial for the $i$ joint is:

where $l_{j1}(t), l_{j2}(t)$ and $l_{j3}(t)$ represent the trajectory of the 3-5-3 spline polynomial, respectively. According to the constraints of the multi-degree polynomial, the polynomial is solved. The unknown coefficient $a$ can be obtained. The initial point of each segment of the joint, the displacement of the two path points and the termination point, the continuous displacement, velocity, and acceleration between the path points, as well as the velocity (usually 0) and acceleration of the initial point and the termination point, are known.

2.2.   Solving polynomial interpolated trajectories

The matrix $\boldsymbol{A}$ can be listed according to the constraints, including Eq (2.2). It can be seen from the matrix expression that the constraints and constraints are only related to the time $t$.The value of the coefficient $a$ can be obtained from Relational Eq (2.4). $t_{1}, t_{2}, t_{3}$ in Eq (2.2) respectively represent the time of the 3-segment polynomial interpolation of the $i$ joint, respectively.

2.3.   Establishment of objective function and kinematic constraints

Under the premise of considering various constraints in kinematics to improve the working efficiency of the manipulator, reduce energy consumption, reduce the impact, and achieve the multi-objective comprehensive optimization of time, energy consumption, and impact, the optimization objective function is defined as:

In the formula, $S_1$ is the total running time of the manipulator, the sum of the time intervals between the position sequence nodes, and the work efficiency index of the manipulator. $S_2$ is the average acceleration of the joint, which is the energy consumption index. $S_3$ is the average jerk in the joint motion and the pulsating impact index. $T$ is the total running time of the manipulator. $a_i$ is the acceleration of the joint, and $J_i$ is the jerk of the joint. The manipulator must be subject to constraints such as position, velocity, acceleration, and jerk during operation. $t_j$ is the time of the jth point, and $N$ is the total number of joints of the manipulator. Besides, the constraints are defined as

wherein, $g_{1}(t)$ is the joint displacement constraint; $g_{2}(t)$ represents the joint velocity constraint; $g_{3}(t)$ indicates joint acceleration constraints; $g_{4}(t)$ suggests joint acceleration constraints; $Q_{jmax}$ is the maximum displacement of the joint, $q_{i}(t)$ is the angle of the ith joint; $V_{jmax}$ is the maximum speed of the joint, $v_{i}(t)$ is the speed of the ith joint; $A_{jmax}$ is the maximum acceleration of the joint, $a_{i}(t)$ is the acceleration of the ith joint; and $J_{jmax}$ indicates the maximum jerk of the joint, $j_{i}(t)$ is the jerk of the ith joint.

2.4.   Robotic arm trajectory optimization based on multiverse algorithm

2.4.1.   Multiverse algorithm

The main idea of the multiverse optimization algorithm $(\mathrm{MVO})$ is to build a model based on the three main concepts of the multiverse theory, namely, white holes, black holes, and wormholes. A white hole is a celestial body that only emits but does not absorb. A black hole absorbs everything in the universe. A wormhole is like a space-time tunnel connecting a white hole and a black hole, sending individuals to any corner of the universe, even from one universe to another. The multiverse reaches a stable state through the three. The candidate solution is defined as the universe, and the fitness of the candidate solution is the universe's expansion rate. Additionally, the MVO algorithm can be divided into the following steps:

1) Set up a group consisting of $N$ universes to search in the $D$-dimensional target space, and initialize it.

2) To establish a mathematical model between white and black holes and exchange the objects of the universe, the roulette mechanism is adopted. In each iteration, the universe is sorted according to the expansion rate (fitness) of the universe. Besides, a roulette wheel is used to select one White hole.

where $x_{i}^{j}$ is the $j$ argument for the $i$ universe, $\mathrm{NI}\left(X_{i}\right)$ shows the normalized expansion rate of the $i$ universe, $r_{1}$ is a random number between $[0, 1]$, and $x_{k}^{j}$ is the $j$ parameter of the $k$ universe selected by the roulette mechanism.

3) The wormhole existence probability $WEP$ increases linearly, the travel distance rate $TDR$ decreases continuously in the iterative process. Thus, a more accurate local search can be performed within the obtained global optimum range.The adaptive formulas of two coefficients are Eqs (2.9) and (2.10):

where $WEP_{\mathrm{min}}$ is the minimum value of $WEP$, $WEP_{\mathrm{max}}$ is the largest value of $WEP$, $l$ is the current number of iterations, $L$ is the maximum number of iterations, and $p$ defines the development precision during iteration. It indicates the higher the value of $p$, the faster the local search.

4) Update the cosmic position and find the optimal individual.

When $r_{2} < WEP$,

When $r_{2}\geq WEP$,

where $x_{j}$ represents the $j$ parameter of the best universe currently formed, $\mathrm{ub}_{j}$ and $\mathrm{lb}_{j}$ represents the upper and lower bounds of the $j$ variable, and $r_{2}, r_{3}$ and $r_{4}$ suggest random numbers between $[0, 1]$.

2.4.2.   Multiverse algorithm improvements

Wormhole existence probability $WEP$ increases linearly in the iterative process. In order to speed up the convergence speed of the algorithm, $WEP$ is currently improved on the basis of traditional MVO, and its expression is:

where $\alpha_0$ is the initialization learning rate, the learning rate at the $t$ iteration is $\alpha_t$, and $\beta$ is the attenuation rate.

MVO has great advantages in local search, but it is not difficult to fall into local optimum, resulting in insufficient accuracy, whereas the genetic algorithm has strong global search ability. In order to overcome the above shortcomings, the population mutation strategy is currently used to enhance the diversity of the population and avoid falling into local optimum. Improvement strategies are described below:

where $\operatorname{Mut}\left(X_{i}\right)$ is the rate of variability of the $i$ universe, $r_4$ is a random number between $[0, 1]$, and $k$ refers to a random integer between $[1, d]$.

In this paper, an adaptive adjustment strategy for parameter p is proposed. The specific steps of the strategy are as follows:

1) Calculate the average value of the normalized expansion rate $f_{m}$.

2) Take out the universes whose expansion rate is greater than the average expansion rate, and calculate their average expansion rate $f_{m1}$. Assign a universe with an expansion rate greater than $f_{m1}$ to a maximum $p$.

3) Take out the universes whose expansion rate is less than the average expansion rate, and calculate their average expansion rate $f_{m2}$. Particles with an expansion rate less than $f_{m2}$ are assigned the smallest $p$.

4) According to the adaptive adjustment strategy, the particles whose fitness value is between $f_{m1}$ and $f_{m2}$ are given values that vary linearly between the maximum and minimum $p$ of $f_{m1}$ and $f_{m2}$.

When there are multiple objectives, due to the conflict and incomparability between objectives, a solution may be the best on one objective and the worst on other objectives. These solutions, which, while improving any objective function, necessarily weaken at least one other objective function, are called non-dominated solutions or Pareto solutions.

The comparison in the single objective only compares the value of the single objective function. In the multi-objective, if all the objective function values of the A solution are better than the B solution, then A is then said to dominate B. If A solution cannot satisfy all target values that are better than B, and no target values of B are better than A, then A and B belong to the same level.

1) First, generate the population, obtain the first dominant solution, and use the roulette wheel to select the optimal universe in the first dominant solution.

2) Update the universe according to the original MVO algorithm, the expansion rate belonging to the first dominant solution takes the smaller value, and the expansion rate which does not belong to the feasible solution takes 1.

3) The newly obtained universe is merged with the parent population, and the first dominant solution is updated. Additionally, redundant and inferior solutions are deleted.

2.4.3.   Test function and validation

After testing the ZDT1 function and ZDT2 function, the effect is shown in Figures 1 and 2.

Ten classical benchmark functions, as shown in Table 1, are used to evaluate the performance of MVO, PSO and IMVO are the algorithms for comparison. In the ten functions, $f_1\sim f_4$ are unimodal test functions that can test the optimization accuracy of the algorithms, $f_5 \sim f_8$ are multimodal test functions that can test the global optimization ability and convergence speed of the algorithms, $f_9 \sim f_{10}$ is ill-conditioned test functions that can test the exploration and exploitation capabilities of the algorithms.

In Figure 4, it can be seen that the optimization speed and accuracy of IMVO added with the parameter adaptive adjustment and wormhole probability curve are significantly improved compared with MVO, which shows that improving the population diversity and the ability to jump out of the local solutions is incredibly beneficial. $f_1\sim f_4$ prove that the optimization accuracy of IMVO is not weaker than that of PSO, and the convergence speed is faster than that of MVO. $f_5 \sim f_8$ prove that the global optimization ability and convergence speed of IMVO is higher than that of MVO and PSO. $f_9 \sim f_{10}$ prove that the exploration capabilities of IMVO can be maintained at a high level.

All of the above show that the IMVO algorithm is effective in dealing with unimodal, multimodal, and ill-conditioned test functions, and it has better optimization accuracy and speed, which is very helpful in solving the problems to be optimized in engineering practice.

3.   Algorithm flow

Details as follows:

1) Initialization parameters $\boldsymbol{X}$.

2) Set Normalized Expansion Ratio $\mathrm{NI}(X)$ and initialize $\alpha_t$, WEP and $p$.

3) Calculate and normalize various fitness of population individuals to find the optimal fitness and normalize it.

4) Black and white hole transfer and update position, $\alpha_t$ and $p$.

5) Individual variation in populations.

6) If the termination condition is not met, go to Step 3. otherwise, end.

4.   Simulation

4.1.   Testing environment

Test environment: the hard disk running environment is AMD Ryzen 7 5800X 8-Core Processor, the memory is 16GB, the software running environment is Windows11 system, and the running software is MATLAB 2021a.

4.2.   Model building

The D-H parameters of the PUMA560 robot arm are shown in Table 2.

As shown in Table 2, $\theta_{i}$ is the joint angle around the $Z$ axis; $d_{i}$ is the joint distance; $a_{i}$ is the connecting rod length, and $\alpha_{i}$ is the torsion around the $X$ axis.

Using the D-H parameters in Table 2, create the model of the PUMA560 manipulator in the four postures in Table 3 in Matlab simulation software. The models are shown in Figure 5.

The joint constraints are shown in Table 4. In order to objectively compare the optimization performance of each algorithm, all algorithms select the same initial population, and the average fitness value, optimal fitness value, standard deviation of fitness value, and average running time of each algorithm running independently for ten times are counted, and the algorithms are comprehensively evaluated. Specific parameter settings are shown in Table 5.

4.3.   Simulation results and analysis

Trajectory optimization simulation is carried out by MOPSO and MOMVO, and Pareto optimal solution is set with time, Energy, and impact as optimization objectives are obtained. As shown in Figure 6, Energy is positively correlated with impact, and both of them are negatively correlated with time, which jointly restricts the optimal performance of time. Because it is necessary to increase the speed to pursue the shortest time, the energy consumption and impact will increase. In practice, it is necessary to select an appropriate solution in the Pareto optimal solution set according to the specific operating conditions of the manipulator.

Since a single objective cannot compare the pros and cons of the solutions in the Pareto optimal solution set, it is essential to define a weighted objective function in order to compare the iterative results. The normalized weighted objective function is defined as

$N_1$, $N_2$, and $N_3$ take different values, respectively, and establish different objective functions for comparison. The comparison results of each algorithm are shown in Tables 6 and 7.

As seen from Figure 9, the improved MVO algorithm quickly obtained iterative convergence in the early stage, and the convergence was stable later.

In addition, as shown in Tables 6 and 7, in order to minimize the error, in a statistical sense, IMVO is better than MVO, PSO, and SAPSO in terms of mean and optimal values of four sets of models constructed with different values of $N_1$, $N_2$ and $N_3$ respectively. The standard deviation of the four data sets is also the smallest IMVO, meaning that the results are more concentrated and the algorithm is more stable. Only the running time of MVO is optimal in the four sets of data because IMVO is an improvement on the MVO, the long running time is normal, and the running time of IMVO is shorter than PSO and SAPSO. It is only about 1s more than the MVO, which is allowed in engineering. Figure 7 shows that the parameters optimized by IMVO can meet the constraints in Table 4.

The simulation results show that when the constraint conditions are satisfied, the convergence speed and accuracy of the improved MVO are obviously higher than those of MVO, PSO and SAPSO.

5.   Conclusions and future works

The kinematics analysis of the joint space of the 6-DOF manipulator is performed, and the trajectory model of the manipulator in the joint space is constructed by 3-5-3 polynomial interpolation.

Taking the running time, energy consumption, and impact of each joint of the manipulator as the optimization goal, the optimization is conducted under the specified constraints, and the Pareto optimal solution set with good diversity and convergence is obtained. Customizing the normalized weighted objective function and selecting the optimal solution meeting the conditions from the Pareto optimal solution set is essential.

Compared with MVO, PSO, and SAPSO, the improved MVO algorithm optimizes the iteration result, improves the optimization structure, and improves the algorithm efficiency.

This paper only discusses 3-5-3 polynomial interpolation. As future work, IMVO can optimize other interpolation functions such as cubic polynomial and quintic polynomial. IMVO can also be applied to optimization in other fields, such as parameter structure optimization of machine learning, fuzzy control optimization, inverse kinematics solution of the manipulator, etc.

Conflict of interest

The authors declare there is no conflict of interest.