2.1. Kinematic analysis of robotic arm

This section mainly focuses on the kinematic analysis of the PUMA560 and xArm6 six-degree-of-freedom robotic arms. Figure 1 illustrates the link structure of the PUMA560 robotic arm.

Figure 1.

Link structure diagram of PUMA560 robotic arm.

The Denavit-Hartenberg (D-H) parameters of the PUMA560 robotic arm are constructed based on its kinematic parameters, as shown in Table I.

Table I.

DH parameters with symbolic variables.

The key parameters are as follows: $d_3 = 0.149$ , $d_4 = 0.433$ , $a_2 = 0.4318$ , $a_3 = 0.0203$ .

Figure 2 shows the xArm6 robotic arm. Like the PUMA560 robotic arm, this arm is also a 6-axis robotic arm, with three joints used for position adjustment and three joints used for posture adjustment.

Figure 2.

Link structure diagram of xArm6 robotic arm.

The Denavit-Hartenberg (D-H) parameters of the xArm6 robotic arm are constructed based on its kinematic parameters, as shown in Table II.

Table II.

DH parameters of the xArm robot.

Tx_offset is the offset angle from the mathematical position to the mechanical position.

Since both the PUMA560 and xArm6 robotic arms are six-degree-of-freedom robotic arms, their link transformation formulas are the same, as given by Equation (1).

(1) \begin{equation} ^{i-1}_iT=\begin{pmatrix}c\theta _i&-s\theta _ic\alpha _i&s\theta _is\alpha _i&a_ic\theta _i\\s\theta _i&c\theta _ic\alpha _i&-c\theta _is\alpha _i&a_is\theta _i\\0&s\alpha _i&c\alpha _i&d_i\\0&0&0&1\end{pmatrix} \end{equation}

where s $\theta _i$ is the abbreviation for sin $\theta _i$ , c $\theta _i$ for cos $\theta _i$ , s $\alpha _i$ for sin $\alpha _i$ , and c $\alpha _i$ for cos $\alpha _i$ . By substituting the D-H parameters into the link transformation equation, the transformation matrix from link 6 to link 0 can be obtained, as shown in Equation (2).

(2) \begin{equation} \left \{\begin{array}{c} { }_6^0 T={ }_1^0 T_2^1 T_3^2 T_4^3 T_5^4 T_6^5 T \\[5pt] { }_6^0 T=\left [\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} r_{11} & r_{12} & r_{13} & P_x \\[5pt] r_{21} & r_{22} & r_{23} & P_y \\[5pt] r_{31} & r_{32} & r_{33} & P_z \\[5pt] 0 & 0 & 0 & 1 \end{array}\right ]=\left [\begin{array}{c@{\quad}c} R_{3 \times 3} & P \\[5pt] 0 & 1 \end{array}\right ] \end{array}\right . \end{equation}

where $_6^0 T$ represents the end-effector pose matrix of the robotic arm, $R_{3 \times 3}$ is the rotation matrix of the end-effector, and P is the translation vector of the end-effector.

2.2. Fitness function for inverse kinematics solution

Although analytical solutions to the IK problem exist for certain robot configurations – such as manipulators with a spherical wrist structure like the PUMA – these solutions are highly dependent on the specific kinematic architecture and may not be applicable to robots with more complex or arbitrary structures. To improve the generality and flexibility of the solution method, this study reformulates the IK problem as a nonlinear optimization problem. The improved DBO algorithm is then employed to numerically solve the joint configurations that minimize the pose error of the robot end-effector.

For the minimization problem, a fitness function must first be constructed. Since the robotic arm’s pose parameters consist of a 3×3 rotation matrix and a 3×1 translation vector, the fitness function is defined as shown in Equation (3).

(3) \begin{equation} \left \{\begin{array}{c} \operatorname {rot}_{\text{diff}}=\left \|\textrm{T}_{\text{rot}}-\textrm{At}_{\text{rot}}\right \| \\[1ex] \operatorname {pos}_{\text{diff}}=\left \|\textrm{T}_{\text{pos}}-\textrm{At}_{\text{pos}}\right \| \\[1ex] \textrm{F}=\textrm{rot}_{\text{diff}}+\textrm{pos}_{\text{diff}} \end{array}\right . \end{equation}

where ${T}_{rot}$ represents the actual rotational state of the robotic arm, ${At}_{rot}$ is the target rotational state, ${rot}_{diff}$ is the difference between the actual and target rotational states; ${T}_{pos}$ represents the actual position vector of the end-effector in world coordinates, ${At}_{pos}$ is the target position vector of the end-effector in world coordinates, ${pos}_{diff}$ is the difference between the actual and target positions; F represents the total difference between the actual pose matrix and the target pose matrix.

As derived from Equation (3), when the actual pose of the robotic arm matches or closely approximates the target pose, the fitness function F should equal 0 or approach 0.

3.1. Rolling behavior of dung beetles

In nature, dung beetles roll animal dung into balls and use polarized light for navigation, rolling the dung balls in a straight line toward an ideal location.

To simulate the behavior of dung beetles rolling dung balls, rolling beetles move in a specific direction across the entire search space. The position update equation for rolling beetles is given by Equation (4).

(4) \begin{equation} \left \{\begin{array}{c} x_i^{t+1}=x_i^t+\alpha \cdot k \cdot x_i^{t-1}+b \cdot \Delta x \\ \Delta x=\left |x_i^t-x_{\text{worst }}^t\right | \end{array}\right . \end{equation}

where $t$ is the current iteration, $x_i^t$ is the position of the $i$ -th beetle at iteration $t$ , $k$ is a constant deflection coefficient with a value range of $[0, 0.2]$ , $b$ is a constant between $[0, 1]$ , $\alpha$ is a natural coefficient, either $1$ or $-1$ , and $x_{\text{worst}}^t$ is the worst global position at iteration $t$ .

When rolling beetles encounter obstacles, they perform a “dancing” behavior on top of the dung ball to determine a new direction before continuing to roll. This behavior is simulated using the tangent function, as shown in Equation (5).

(5) \begin{equation} x_i^{t+1}=x_i^t+\tan (\theta )\left |x_i^t-x_i^{t-1}\right | \end{equation}

When $\theta = 0$ , $\pi /2$ or $\pi$ , $\tan (\theta )$ becomes undefined, meaning the beetle’s position remains unchanged.

3.2. Breeding behavior of dung beetles

Female dung beetles roll dung balls to safe locations and hide them to create a suitable environment for offspring. This behavior inspired the boundary constraint strategy used to simulate the egg-laying area of dung beetles, expressed in Equation (6).

(6) \begin{equation} \left \{\begin{array}{l} \textrm{Lb}^*=\max \left \{\textrm{x}^* \cdot (1-\textrm{R}), \textrm{Lb}\right \} \\ \textrm{Ub}^*=\min \left \{\textrm{x}^* \cdot (1+\textrm{R}), \textrm{Ub}\right \} \end{array}\right . \end{equation}

where $R=1-t/T$ , and T is the maximum number of iterations, Lb and Ub are the lower and upper bounds of the search space, $x^*$ is the global optimal position in the current generation, $Lb^*$ and $Ub^*$ are the lower and upper bounds of the egg-laying region.

Once the egg-laying region is determined, female beetles lay a single breeding ball within this area during each iteration. The position of the breeding ball changes dynamically with the egg-laying region, as given by Equation (7).

(7) \begin{equation} B_i^{t+1}=\textrm{x}^{\wedge }+\textrm{b}_1 \cdot \left (B_i^t-\textrm{Lb}^*\right )+\textrm{b}_2 \cdot \left (B_i^t-\textrm{Ub}^*\right ) \end{equation}

where $B_i^t$ is the position of the i-th breeding ball at iteration t, $b_1$ and $b_2$ are two independent random vectors of size 1xD, where D is the problem dimension.

As shown in Figure 3, the blue points represent the positions of the breeding balls, the black points represent the local optimum $x^*$ , and the red points represent the upper and lower boundaries of the region.

Figure 3.

Boundary selection strategy diagram.

3.3. Foraging behavior of juvenile beetles

Once the breeding balls successfully hatch, juvenile beetles emerge and forage within the local area. As the iterations progress, the foraging area of juvenile beetles dynamically updates, as shown in Equation (8).

(8) \begin{equation} \left \{\begin{array}{l} \textrm{Lb}^{\wedge }=\max \left \{\textrm{x}^{\wedge } \cdot (1-\textrm{R}), \textrm{Lb}\right \} \\ \textrm{Ub}^{\wedge }=\min \left \{\textrm{x}^{\wedge } \cdot (1+\textrm{R}), \textrm{Ub}\right \} \end{array}\right . \end{equation}

where $R=1-t/T$ , and T is the maximum number of iterations, $x^\wedge$ is the current local optimal position, $Lb^\wedge$ and $Ub^\wedge$ represent the lower and upper bounds of the foraging area.

The position update for juvenile beetles is given by Equation (9).

(9) \begin{equation} x_i^{t+1}=x_i^t+C_1 \cdot \left (x_i^t-\textrm{Lb}^{\wedge }\right )+C_2\left (x_i^t-\textrm{Ub}^{\wedge }\right ) \end{equation}

where $x_i^t$ is the position of the i-th juvenile beetle at iteration t, $C_1$ is a random number following a normal distribution $C_1 \sim N(0,1)$ , $C_2$ is a random vector with values in the range (0,1).

3.4. Stealing behavior of dung beetles

In the dung beetle world, the thieves are referred to as stealing beetles, which steal dung balls rolled by other beetles. From the equations provided earlier, $x^*$ and $x^\wedge$ represent the optimal positions for stealing dung balls. The position update for the stealing beetles is expressed as Equation (10).

(10) \begin{equation} x_i^{t+1}=\textrm{x}^{\wedge }+\textrm{S} \cdot \textrm{g} \cdot \left (\left |x_i^t-\textrm{x}^*\right |+\left |x_i^t-\textrm{x}^{\wedge }\right |\right ) \end{equation}

where $x_i^t$ represents the position of the i-th stealing beetle at iteration t, $g$ is a random vector of size 1xD drawn from a normal distribution $g \sim N(0,1)$ , $S$ is a constant.

Figure 4 illustrates the proportion of different types of dung beetles in the population. The blue section represents rolling beetles, accounting for 20% of the total population. The red section corresponds to breeding beetles, also comprising 20%. The yellow section indicates juvenile beetles, which make up 25%, and the green section represents stealing beetles, constituting 35% of the overall population.

Figure 4.

Various dung beetle types in the population.

4.1. Sobol sequence initialization

Traditional swarm intelligence optimization algorithms use generation of random numbers distributed in the solution space when generating the initial population. In contrast, random numbers generated based on modern CPUs are pseudo-random, with duplication and aggregation, resulting in reduced population diversity. Sobol sequences are low-discrepancy sequences that are more evenly distributed in space than pseudo-random sequences. Good overall coverage is ensured by precisely controlling the position in space of each newly generated point.

Figure 5.

Sobol sequence.

Figure 6.

Pseudorandom sequence.

As shown in Figures 5 and 6, a comparison of using Sobol sequences and pseudo-random sequences to generate population initiation locations in 2D space. Where Figure 4 shows the location of dung beetle populations generated using Sobol sequences, it can be seen that they are able to be evenly distributed and somewhat differentiated in any area. Figure 5 shows the distribution of populations generated using pseudo-random sequences, and it can be seen that multiple population clusters occur in the figure, as well as some areas that are not covered.

Sobol sequences are a class of low-discrepancy sequences based on binary reflections that are commonly used to generate uniformly distributed sampling points in high-dimensional spaces. For a Sobol sequence of dimension d, the generating process primitive polynomials and the number of directions. The formula for the mth point of the Sobol sequence is shown in Equation (11).

(11) \begin{equation} x_m=\sum _{i=1}^n\frac {v_{i,j}}{2^i} \end{equation}

where $v{(i,j)}$ denotes the number of directions in the $j$ dimension and $i$ denotes the order of the binary bits.

4.2. Boundary optimization strategy

From the convergence factor $R=1-t/T$ in the reproductive dung beetle boundary formula, we know that as the number of iterations increases, the area available for reproduction decreases. In the early iterations, the value of R is larger, leading to a larger optimal breeding area, enhancing the global search ability of the algorithm. As R decreases over time, the search area shrinks, which strengthens the local search ability in later iterations.

To balance global search in the early stages and local search in the later stages, we redefine the convergence factor as a piecewise function. Initially, the convergence remains unchanged, preserving the global search capability. Later, a logistic function is introduced to adjust the convergence speed of R, enhancing the algorithm’s local search capacity. The logistic function is described in Equation (12).

(12) \begin{equation} \begin{cases}\text{R}=1-\frac {\text{t}}{\text{T}} , \text{t}\leq \text{T}/2\\\text{R}=1-\displaystyle\frac {1}{1+\text{e}^{-\text{k}\cdot \left (\frac {\text{T}}{2}-\text{t}\right )}} , t\gt T/2\end{cases} \end{equation}

where $T$ is the number of iterations, $t$ is the current iteration, and $k$ is a constant, set to –0.01 in this study.

Figure 7.

Variation curve of convergence factor.

As shown in Figure 7, the logistic function maintains an appropriate convergence speed during the early iterations, resulting in a larger optimal search area and improved global search capability. In the later stages, the logistic function accelerates the search process, shrinking the optimal foraging area and enhancing local search, thereby improving the convergence speed.

4.3. Hybrid directional perturbation strategies

During the search process, the Rolling Dung Beetle and Stealing Dung Beetle will explore forward along a straight line path based on the paths generated by Equation (4) and Equation (10). Such a search method greatly diminishes their ability to explore, and they tend to miss opportunities to find better solutions. A hybrid directional perturbation strategy is proposed for this problem, which perturbs the forward direction of the rolling dung beetle and the stealing dung beetle to increase the chances of finding an optimal solution.

Conventional normally distributed perturbations can produce moderate perturbations that improve search. However, when faced with complex optimization problems, especially with multiple local optimal solutions, the normal distribution may not be able to jump out of the local optimum. T-distributed perturbations have thicker tails compared to normal-distributed perturbations, increasing the probability of generating large perturbations during the search process and improving the algorithm’s ability to jump out of local optima. T-distribution is a symmetric distribution with a probability density function as shown in Equation (13). where M is the degree of freedom, which determines the shape of the distribution, the t-distribution gradually converges to a normal distribution as the degree of freedom increases.

(13) \begin{equation} f(n)=\frac {\Gamma \left (\frac {m+1}{2}\right )}{\sqrt {m\pi }\Gamma \left (\frac {m}{2}\right )}\Bigg (1+\frac {n^2}{m}\Bigg )^{-\frac {m+1}{2}} \end{equation}

During the iteration process, the t-distribution perturbation is introduced into the update of the search position with a perturbation factor formulated as equation (14).

(14) \begin{equation} w=\mu _{w}+\sigma _{w}\cdot t_{m} \end{equation}

where $\mu _w$ is the offset. $\sigma _w$ is the scale parameter, which controls the magnitude of the perturbation. $t_m$ is a random variable drawn from a t-distribution with m degrees of freedom.

The variation curve of the disturbance factor at this time is shown in Figure 8.

Figure 8.

w perturbation curve.

Dynamic weight perturbations can change the perturbation factor as the number of iterations increases. Having a perturbation factor with a strong effect at the beginning of the iteration enables exploration over a large range. With the increase in the number of sum iterations, the perturbation factor gradually converges to 1, the perturbation effect is weakened, and a more detailed search can be carried out within a certain range, increasing the convergence accuracy. The dynamic weight perturbation factor formula is shown in equation (15). The dynamic perturbation factor w curve is shown in Figure 9.

(15) \begin{equation} w=e^{-\frac {T-t}{T}} \end{equation}

where $T$ is the maximum number of iterations and $t$ is the current number of iterations.

Figure 9.

Disturbance factor w.

T-distribution has thicker tails compared to the normal distribution, which means a greater probability of producing values far from the mean. The ability to introduce large random fluctuations allows the ball-rolling dung beetle to have a greater chance of jumping out of the current local optimum region and exploring unknown regions. The dynamic weight perturbations decay exponentially, and the magnitude of the perturbations is on the small side, which enables the Stealing Dung Beetle to explore the region in greater detail and helps the algorithm to approximate the global optimum more quickly. The two perturbation strategies were used in combination, with two consecutive t-distribution perturbations followed by a single dynamic weight perturbation on the parameters k and b for the ball-rolling dung beetle and on the parameter S for the stealing dung beetle, which in turn controlled the search routes of the two dung beetle types. As shown in Equation (16).

(16) \begin{equation} \begin{cases} k=k*w\\b=b*w\\S=S*w \end{cases} \end{equation}

Two consecutive perturbations of the t-distribution, while increasing the dung beetle’s global search ability, may also cause the search to be too far off course. A dynamic weight perturbation at this point can correct the search direction to avoid the problem of too much deviation. At the same time, the addition of dynamic weight perturbation makes the dung beetle search sufficiently in a small area and accelerates the convergence of the algorithm, and the hybrid strategy works better than any of the perturbation strategies alone.

4.4. Levy flight strategy

In order to solve the problem that the original DBO algorithm suffers from the tendency to fall into local optimums, the Levy flight strategy is introduced. The Levy flight strategy is a stochastic wandering method with heavy-tailed properties. Unlike Brownian motion, Levy flight allows for episodic large step-size jumps, enabling the dung beetle population to search over a larger area, improving the algorithm’s global search capability and enhancing the ability to jump out of local optima.

The step size s of the Levy flight is generated by Equation (17).

(17) \begin{equation} s=\frac {u}{v^{\frac {1}{\beta }}} \end{equation}

where $v$ and $u$ obey a normal distribution: $v\sim (0,1)$ , $u\sim N(0,\sigma _u^2)$ , $\sigma _u$ is the scale parameter, determined by Equation (18).

(18) \begin{equation} \sigma _u=\left (\displaystyle\frac {\Gamma (1+\beta )\sin \displaystyle\frac {\pi \beta }{2}}{\Gamma \left (\displaystyle\frac {1+\beta }{2}\right )\beta 2^{\displaystyle\frac {\beta -1}{2}}}\right )^{\displaystyle\frac {1}{\beta }} \end{equation}

$\Gamma$ is the gamma function. $\beta$ is the Levy distribution index, and $\beta$ = 1.5.

Figure 10.

Levy flight path.

During the iteration of the algorithm, the Levy flight strategy is used to generate new candidate solutions that jump out of the local optimum by perturbing the current optimal solution. In contrast to other perturbation methods, the Levy flight strategy has a perturbation rule for short-range perturbations and occasional long-range perturbations. The short-range perturbation ensures that the algorithm is able to perform a detailed search of the local space and avoids frequent perturbations that affect the convergence performance of the algorithm. And occasional long-range perturbations help the algorithm to go beyond the local optimum and explore a larger range. Levy flight jump trajectory is shown in Figure 10.

4.5. Complexity of ECDBO algorithm

In the initial phase of the algorithm, a Sobol sequence initiation strategy is introduced, which generates a Sobol sequence initial solution for some of the populations with a complexity of O(N * D), where N is the population size and D is the problem dimension. The bounded optimization strategy is used to change the complexity to O(N * D). Dynamic weight perturbations and t-distribution perturbations are used to update the parameters, and since the dynamic weights w are obtained by simple scalar computation and the generation of t-distributed random numbers is constant time in expectation, the complexity of this part is O(1), which has little impact on the overall time complexity. The Levy flight strategy is used to perturb the globally optimal solution, which involves operations on D-dimensional vectors with O(D) complexity. Taken together, the main time consumption of the ECDBO algorithm in each iteration remains the updating and fitness assessment of the population with a complexity of O(N * D), and the fitness ranking with a complexity of O(N log N). Thus, the overall time complexity of the algorithm is O(T * (N * D + N log N)), where T is the number of iterations. The overall time complexity of the DBO algorithm is also O(T * (N * D + N log N)). Compared to the DBO algorithm, the convergence speed and accuracy of the algorithm are substantially improved without increasing the time cost.

4.6. ECDBO algorithm workflow

The flow of Enhanced Chaotic Dung Beetle Optimization (ECDBO) algorithm can be summarized in the following pseudocode.

Algorithm 1

ECDBO Algorithm

5.1. CEC2017 function experimentation and analysis

The CEC2017 test set consists of 30 single-objective functions, including single-peak functions (F1-F3), simple multi-peak functions (F4-F10), hybrid functions (F11-F20), and combined functions (F21-F30). Figures 11 and 12 show the ranking of the average results of the eight algorithms after 30 independent runs in different dimensions.

Figure 11.

Average ranking under 30 dimensions.

Figure 12.

Average ranking under 100 dimensions.

5.1.1. Numerical analysis of CEC2017 results

The experimental results of ECDBO algorithm in 30 and 100 dimensions are recorded in Tables III and IV. It shows the mean, standard deviation, and ranking of each algorithm. The data is analyzed as follows:

• • ECDBO algorithm is ahead or even significantly ahead of the other seven algorithms in both the 30- and 100-dimensional tests for the single-peak functions F1-F3, which indicates that the ECDBO algorithm excels in its ability to find the global optimum. Overall, the ECDBO algorithm performs well in solving the single-peaked function problem.

• • The results of the simple multi-peak function F4-F10 test show that the ECDBO algorithm outperforms the other seven algorithms. In the 30-dimensional test, the ECDBO algorithm outperforms the other algorithms on the functions F4, F7, and F8, and ranks first. It is slightly weaker than the DBO algorithm on functions F5 and F9 and ranks second. It is only weaker than the iCSPM algorithm in functions F6 and F10, and ranks second. In the 100-dimensional test, the ECDBO algorithm outperforms the other seven algorithms in functions F4, F5, F7, F8 and F10, and ranks first, and is only slightly weaker than the DBO algorithm in function F6, and is weaker than the DBO, PSO and WOA algorithms in function F9, and ranks fourth. These results verify the superiority of ECDBO algorithm in solving multi-objective problems, showing very good performance.

• • ECDBO has very good convergence performance on hybrid functions F11-F20. In the 30-dimensional test, ECDBO tops the list of 7 hybrid functions, which is attributed to the improved search strategy and the ability to jump out of local optima. In the 100-dimensional test, ECDBO outperforms the other algorithms in 9 functions, and is almost equal to the iCPSM algorithm only in function F20.

Table III.

CEC2017 test results, dim = 30.

Table IV.

CEC2017 test results, dim = 100.

• • ECDBO also performs well in the combined functions of F21-F30. In the 30-dimension test, it ranks first in 8 test functions, and is only slightly weaker than the iCSPM algorithm in F23 and F29. In the 100-dimension test, it ranks first in 10 test functions ahead of the other algorithms. These results show that the ECDBO algorithm is highly competitive in solving difficult and complex problems. Its excellent global exploration and local exploitation capabilities, as well as its remarkable ability to jump out of local optima, provide strong support for solving complex problems.

5.1.2. CEC2017 convergence curve analysis

Algorithm convergence analysis is an important metric for evaluating the algorithm’s global and local search capabilities. Figures 13 and 14 show the convergence curves of ECDBO with the other seven comparison algorithms on the 30 functions of the CEC2017 test set.

In the 30-dimensional test, the ECDBO algorithm is significantly ahead of the other seven algorithms on functions F1, F2, F4, F12, F13, F14, F19, F21, F22, F24, F26, and F30, which can be attributed to the strong global exploration and local development capability of the ECDBO algorithm. Although the other algorithms may converge faster than ECDBO at the beginning of the iteration, they are soon overtaken by ECDBO, and the gap is getting wider and wider due to the lack of development. For functions F3, F7, F8, F11, F15, F18, F24, F25, F27, F28, and F29, the ECDBO algorithm and other better-performing algorithms (e.g., iCSPM, DBO, and PSO algorithms) alternately take the lead at the beginning of the iteration, and at the later stage of the iteration, the ECDBO algorithm eventually takes the lead due to its excellent development capability.

In the 100-dimensional test, the ECDBO algorithm performs well on functions F1, F4, F5, F10, F12, F13, F14, F15, F16, F18, F19, F21, F22, F23, F25, F28, and F30, and the final convergence result is significantly ahead of the other algorithms. On functions F2, F3, F5, F7, F11, F17, F20, F24, F26, F27, and F29, the ECDBO algorithm and other algorithms show alternating leads, but soon the ECDBO algorithm exceeds the other algorithms in terms of convergence speed and accuracy, and maintains a better convergence speed and accuracy in the later iterations. On functions F2, F17, F20, F24, and F26, the ECDBO algorithm and the iCSPM algorithm eventually converge with almost the same iteration speed and accuracy, ahead of the other six algorithms. On function F20, the iCSPM algorithm performs better, while on function F9, the DBO algorithm performs better.

The following results can be obtained from the above analysis:

• • ECDBO algorithm achieves the best results on a larger number of functions, both in terms of mean value and convergence speed. This shows that the ECDBO algorithm has a better ability to cope with a wide variety of optimization problems, reflecting its strong competitiveness.

• • Unlike the other compared algorithms, the ECDBO algorithm finds optimal solutions in multiple functions, reflecting the fact that ECDBO possesses greater adaptability to different problems.

• • When confronted with complex problems, ECDBO has more powerful exploration and exploitation capabilities, and is more likely to allow the algorithm to eventually converge to a globally optimal solution.

Figure 13.

30-dimensional CEC2017 convergence curve.

Figure 14.

100-dimensional CEC2017 convergence curve.

5.2. Parameter sensitivity analysis

In order to verify the parameter sensitivity of the proposed improvement strategy in the ECDBO algorithm, the parameter $k$ of the boundary optimization strategy and the parameter $\beta$ of the Lévy flight strategy are selected for analysis in this paper. To ensure the comprehensiveness and representativeness of the experimental conclusions, the better-performing functions in the single-peak, multi-peak, hybrid and combined functions (F1, F4, F13 and F22, respectively), as well as the corresponding worse-performing functions (F3, F9, F20 and F29, respectively), were selected from the set of test functions of the CEC 2017, and the stability and luberance of the above eight functions were deeply analyzed by solving the functions with different combinations of parameters under different conditions. By solving the above eight functions under different parameter combinations, the performance of the improved algorithm in terms of stability and robustness is analyzed in depth.

To investigate the sensitivity of the proposed algorithm to the parameter $k$ , a series of experiments was conducted using different values: $k$ = −0.01, −0.1, −0.3, −0.5, −0.9. The average fitness curves across eight representative benchmark functions (F1, F3, F4, F9, F13, F20, F22, and F29) are shown in Figure 15.

For functions F1, F3, F4, F13, and F29, the algorithm exhibits consistent convergence behavior for all values of k. The algorithm maintains a stable search dynamics regardless of the change in $k$ values. In F9 and F20, the setting $k = -0.01$ consistently outperforms others in terms of final fitness. In F22, the curves begin to diverge after 400 iterations, with $k = -0.01$ again achieving the best convergence, and $k = -0.9$ performing the worst. These results indicate that smaller-magnitude negative values of $k$ (closer to 0) are better suited to maintaining exploitation capability in complex search spaces, whereas larger negative values may hinder convergence precision. While small performance differences appear in certain cases, especially for complex multimodal functions, the algorithm remains effective under all tested settings. This highlights the robustness and adaptability of the algorithm with respect to parameter $k$ .

To further assess the influence of key parameters on algorithm performance, we conducted a sensitivity analysis on the parameter $\beta$ , with selected values $\beta = 1.0, 1.3, 1.5, 1.7, 1.9$ . The results of the average convergence curves for different settings in the benchmarking functions F1, F3, F4, F9, F13, F20, F22 and F29 are shown in Figure 16. The following insights can be drawn: In F1 and F13, $\beta = 1.5$ achieves superior convergence precision, while lower values converge more slowly and stagnate at suboptimal fitness levels. In F4, all curves show similar trends initially, but the later convergence performance is notably better with larger $\beta$ values ( $\beta$ = 1.5 to 1.9), especially beyond 300 iterations. In F3 and F29, the curves for all $\beta$ values rapidly converge and almost entirely overlap, suggesting the optimization difficulty is low and $\beta$ has negligible impact in these scenarios. Although $\beta$ = 1.7, 1.9 outperforms the default configuration ( $\beta$ = 1.5) for functions F9, F20, and F22, the differences between all functions are not significant. Therefore, $\beta$ = 1.5 provides a robust balance between the different problem types and was chosen as the default configuration. The convergence behavior is significantly more sensitive to the parameter $\beta$ compared to $k$ .

Figure 15.

Experimental results at different $k$ values.

5.3. Wilcoxon rank sum test

Wilcoxon rank sum test was used to test whether the ECDBO algorithm was significantly different from the other seven pairwise algorithms. Using the test results of the eight algorithmic CEC2017 functions as input samples, and tested at a significance level of 0.05. Tables IV and V demonstrate the results of the Wilcoxon rank sum test.

Table V.

30-dimensional CEC2017 Wilcoxon rank sum test.

Figure 16.

Experimental results at different $\beta$ values.

It is assumed that there is no significant difference between the two groups of data to be tested, and the hypothesis is rejected when P < 0.05 proves that there is a significant difference between the two groups of data. As can be seen from Table V and VI, the ECDBO algorithm realizes a significant performance difference in comparison with the results of the other seven algorithms, which indicates that the ECDBO algorithm outperforms the other algorithms in terms of the selected evaluation metrics, and the results are statistically significant. The next step will be to continue to explore the advantages of the algorithm in practical applications.

Table VI.

100-dimensional CEC2017 Wilcoxon rank sum test.

5.4. Inverse kinematics solution experiment

5.4.1. PUMA560 robotic arm experiment

To minimize the possibility of random experimental outcomes, fitness value F was used as the evaluation metric, with ECDBO, IDBO, DBO, and PSO algorithms as the subjects of study. The target pose matrix $A_T$ was used to set up the simulation experiment. At this time, the position and posture of the PUMA560 robotic arm are as shown in Figure 17. The experiment is set with a population size of 60, a maximum iteration count of 1000, and other parameters kept at default values. By comparing the mean value and standard deviation of the fitness values, superiority of the ECDBO algorithm in solving the IK problem was validated. In this experiment, the population size was consistently set to 60, with 1000 iterations, and each algorithm was independently run 30 times. The test results are shown in Table VII, and the fitness value curves of the 30 independent runs are illustrated in Figure 18.

\begin{equation*}A_T=\begin{bmatrix}0.5322&-0.1137&-0.8390&0.0833\\-0.4697&-0.8641&-0.1809&0.1890\\-0.7044&0.4903&-0.5133&-0.4925\\0&0&0&1\end{bmatrix}\end{equation*}

Table VII.

Statistical comparison of fitness values.

Figure 17.

Target position of robotic arm.

As can be seen from Table VII, the adaptation value of ECDBO is significantly smaller than the other three algorithms and is ranked first. This shows that the ECDBO algorithm has higher solution accuracy and is more suitable for IK solution problems.

Figure 18.

Comparison of F-value.

Figure 18 shows the fitness value F curves for 30 independent runs of the four algorithms. From the figure, it is evident that the PSO algorithm’s F values fluctuate within the range of 0 to 0.6, with large amplitude variations, indicating poor stability. The F-value of the DBO algorithm fluctuates within the range of 0 to 0.2, and the amplitude variation is smaller than that of the PSO, but it is still relatively large, indicating that the stability is average. The IDBO algorithm shows fluctuations in the F value within the range of 0 to 0.1, with smaller amplitude variations than PSO and DBO, demonstrating moderate stability. In contrast, the ECDBO algorithm maintains its F values within the range of 0 to 0.022, with minimal fluctuations, reflecting strong stability.

Figure 19.

Actual position of robotic arm.

To verify the feasibility and effectiveness of the ECDBO algorithm for IK solution of robotic arms, PUMA560 robotic arm is used as a research object, $A_T$ is the target position matrix to establish the simulation experiment. IK solution using ECDBO, IDBO, DBO and PSO algorithms, deriving six joint angles and actual pose matrix. PUMA560 robotic arm was modeled in MATLAB, and resulting positional errors were verified through the Robotics Toolbox.

The joint angles obtained by each algorithm are shown in Table VIII. Input the six joint angles obtained from different algorithms into the robot toolbox, and the actual pose of the robotic arm is shown in Figure 19.

Comparing Figures 17 and 19, it can be seen that the results sought by the ECDBO algorithm are the closest to the initially set target, both in terms of the end-effector position and the robotic arm attitude. ECDBO outperforms the other three compared algorithms in handling IK solution problems. This indicates that the improved ECDBO algorithm provides a solution much closer to the target values.

Table VIII.

Joint angle optimization results of three algorithms.

5.4.2. xArm6 robotic arm experiment

To comprehensively evaluate the performance of the proposed IK solving algorithm in continuous trajectory tasks, this paper designs a path tracking experiment based on the xArm6 robotic arm.

Firstly, a point cloud map of the workspace of the xArm6 robotic arm is generated according to the xArm6 robotic arm DH parameters, as shown in Figure 20. An initial path is generated in this point cloud map with a starting point of [30,0,50] and an end point of [−10, −20,50], and the generated initial path is a smooth and continuous 3D curve covering a representative range of attitude changes in the working space of the robotic arm. To facilitate calculation and comparison, 50 points are randomly and uniformly sampled from this path as the target end position sequence for IK solution, and four algorithms are used to perform continuous IK solution and record the experimental results.

The initial population number of each algorithm is set to 60, when the fitness value reaches 0.01 and below, the condition is met, the iteration stops and the results are recorded, or stops when the maximum number of iterations reaches 1,000, at which time the preset condition is not met, the fitness value and the number of times the condition is not met are recorded and used as a comparative analysis.

Figure 20.

Target position of robotic arm.

Table IX demonstrates the running time and number of failures of the four algorithms. Observing the data in the table, it can be learnt that the ECDBO algorithm outperforms the IDBO and DBO algorithms in terms of solution time and is the only algorithm among the four algorithms that passes all 50 IK solutions. Although the PSO algorithm outperforms the ECDBO algorithm in terms of solution time, it also has the relatively highest number of failures, with 11 failures out of 50 solutions.

Table IX.

Execution time and failure counts of different algorithms.

Figure 21 shows the variation of the fitness values for 50 IK solvers, in which ECDBO maintains very low and stable errors at all path points, indicating that its IK solver has the highest accuracy and is the most stable. IDBO has a smaller overall error, but large fluctuations occur at individual path points, showing a certain degree of instability. DBO has a good performance at most points, but there are equally some points where the DBO performs well in most points, but there are also some points with sudden increase in error, which indicates that its convergence ability is insufficient in some complex attitudes. PSO shows the greatest error variability and the poorest stability, suggesting a tendency to get trapped in local optima on complex inverse-kinematics tasks.

Figure 21.

Comparison of inverse kinematics error.

The Figure 22 shows the overlap between the reconstructed path (red dashed line) and the original path (blue solid line) in 3D space for each algorithm. As can be seen from the figure, the path generated by ECDBO is almost completely overlapped with the original path, with very high accuracy, and the overall trajectory is smooth and continuous, which indicates that the algorithm has strong convergence ability and global search ability. The paths generated by IDBO and DBO are also basically close to the target path, but there is a small deviation near some points, which shows a certain degree of instability. In contrast, the paths generated by PSO have obvious jitter and deviation phenomena, and the errors are large at several locations and even fail to accurately approximate the target point, which shows that it is easy to fall into the local optimum in the high-dimensional complex space.

Figure 22.

Comparison of path.

Figure 23.

Comparison of path.

This experiment clearly demonstrates the superior performance of the proposed ECDBO algorithm on the complex path inverse solution problem, and further proves its feasibility and practicability in practical scenarios such as industrial robot path control and attitude tracking.

In order to further verify the usability of the IK solution results in practical control, the 50 joint angle solution sequences computed by the ECDBO algorithm are sequentially imported into the solid xArm6 robotic arm. As shown in the Figure 23, in the experiment, the robotic arm moves point by point according to the inverse solution results in a set order, and finally realizes the overall trajectory tracking from the starting point of the path to the end point smoothly, and completes the path execution task stably under the condition of no external interference, with good continuity of trajectory, no obvious stutter or offset, and the execution process is smooth, which verifies its better engineering adaptability and control effect.