A 6-DOF industrial robot kinematic inverse solution solving method based on an improved COA optimization algorithm

Abstract

Six Degrees of Freedom (6-DOF) industrial robot inverse kinematics solutions are constrained by high-dimensional nonlinearity, multi-solution problems, and traditional algorithms’ local optima tendency, affecting operational accuracy and real-time performance. Both traditional analytical methods and standard swarm intelligence algorithms are confronted with numerous challenges. To address these issues, this paper proposes an inverse solution method based on an improved Coati Optimization Algorithm (COA). The Modified Denavit-Hartenberg (MD-H) kinematic model is introduced to solve the joint singularity problem, and an improved COA based on an elite retention strategy and dynamic weights is proposed. The diversity of the population is enhanced through a mutation mechanism to improve convergence efficiency, and dynamic weights are used to balance position and orientation errors. The experimental results show that when performing inverse kinematics calculations on 10 groups of spatial points, this algorithm outperforms the standard algorithms (Whale Optimization Algorithm (WOA), standard COA, Particle Swarm Optimization Algorithm (PSO), and Genetic Algorithm (GA)) in terms of average number of iterations (40.4 times), average position error (0.0021 mm), and average orientation error (8.5 × 10⁻⁵ rad), meeting the accuracy requirements of industrial robot inverse kinematics and providing a new method for robot inverse solution.

Introduction

Under the backdrop of the rapid development of intelligent manufacturing and industrial automation, 6-DOF industrial robots have become core equipment in fields such as automotive manufacturing, electronic assembly, and aerospace due to their flexible motion capabilities and wide applicability. The inverse kinematics solution of robots, as the foundation for precise trajectory planning and control, ¹ directly affects the operational performance of industrial robots in terms of its solving efficiency and accuracy.^2–4 The inverse kinematics problem is essentially a highly nonlinear and strongly coupled complex solution process,^5–7 with challenges such as multiple solutions, large computational load, and sensitivity to initial values.^8,9 In recent years, various swarm intelligence optimization algorithms have been applied to the inverse kinematics solution of 6-DOF manipulators, such as the Boomerang algorithm based on swarm optimization proposed by Duymazlar and Engin, ¹⁰ and the research by Galan-Uribe and Morales-Velazquez ¹¹ on kinematic optimization of 6-DOF serial robot arms using bio-inspired algorithms.

The existing inverse kinematics solution methods are mainly divided into three categories: analytical methods, numerical iterative methods, and intelligent optimization algorithms. The analytical method is an algebraic solution or geometric decomposition method based on the D-H parameter model.^7,12 Although this method can obtain an exact solution, it has strict requirements for the robot configuration. ⁶ For example, the Pieper criterion requires that the axes of three adjacent joints intersect ¹³ and it requires a large number of trigonometric function operations, resulting in poor real-time performance. Numerical methods such as the Newton iterative method are universal, but they are dependent on the initial value selection and are prone to getting stuck in local optima. The singularity problem of the Jacobian matrix during the iterative process can lead to calculation failure.^14,15 With the rise of intelligent optimization algorithms, PSO,^16,17 GA,^18,19 WOA, ²⁰ and others have been widely applied in inverse solution solving. In the research related to intelligent algorithms and inverse kinematics of robots, He et al. ²¹ proposed an radial basis function-proportional integral differentiation (RBF-PID) adaptive controller that combines radial basis function neural networks with PID control, which improves the solution accuracy but the convergence speed still needs to be optimized. Li et al. ²² proposed an extreme learning machine-sparrow search algorithm-sine cosine algorithm (ELM-SSA-SCA) hybrid algorithm that integrates extreme learning machines with swarm intelligence algorithms, achieving a breakthrough in solution efficiency, but the multi-solution screening mechanism still needs to be improved.

The existing intelligent algorithms still face three core problems when dealing with the inverse kinematics of 6-DOF robots: Firstly, traditional swarm intelligence algorithms such as PSO have the phenomenon of “premature convergence,” and are prone to getting stuck in local optima in high-dimensional solution spaces. Secondly, the multi-solution problem may cause the solution results to exceed the joint movement range, requiring additional constraint processing mechanisms. Thirdly, the real-time requirements of industrial scenarios demand a balance between accuracy and computational efficiency in the algorithm. For instance, Wu et al. ²³ proposed the geese migration optimization (GMO) algorithm, which enhances the global search ability by simulating the migration behavior of geese. However, its convergence accuracy in complex poses still needs improvement. Li et al. ²⁴ proposed the bidirectional opposite differential evolution–cuckoo search (BODE-CS) algorithm, which combines differential evolution and cuckoo search and optimizes the solution process through dynamic parameter adjustment. However, the parameter self-adaptive mechanism still relies on empirical settings. Recent research has further advanced the development of this field. Dereli ²⁵ proposed an improved grey wolf optimization algorithm by modifying the control parameters of the traditional GWO algorithm, which was successfully applied to the inverse kinematics solution of a 7-DOF redundant robot manipulator, significantly enhancing the algorithm’s convergence accuracy. Dereli and Köker, ²⁶ inspired by the motion of a golf ball, introduced the Random Descent Velocity (RDV) inertia weight PSO technique. This method simulates the parabolic descent structure of a golf ball as it approaches the target, converting the particle velocity into a more optimal movement pattern. In the inverse kinematics calculation of a 7-joint robot arm, it achieved an optimization effect thousands of times better than the traditional PSO algorithm.

The COA is an advanced heuristic algorithm inspired by biological intelligence, proposed by Dehghani et al. ²⁷ The algorithm design is inspired by the unique behavior patterns of raccoon groups in nature. When hunting for prey, some members of the raccoon group are responsible for disturbing the prey to make it flee, while others wait at key positions. This collaborative behavior enables them to effectively search for prey in complex environments. This hunting behavior corresponds to the global search process in COA, allowing the algorithm to widely explore the solution space and find potential optimal solutions. When encountering predators, raccoons quickly search for escape routes within their local area, which corresponds to the local exploitation mechanism in the algorithm and enables fine search near better solutions, further improving the quality of the solution. It has shown good performance in function optimization and engineering problems. However, the traditional COA has deficiencies such as insufficient maintenance of population diversity and weak local exploitation ability when dealing with high-dimensional and strongly nonlinear problems like the inverse kinematics of robots.^28,29

To address these issues, this paper proposes a 6-DOF industrial robot kinematic inverse solution method based on the improved COA. By introducing a mutation mechanism based on elitism strategy and a dynamic logarithmic weight strategy, the global search ability and convergence accuracy of the algorithm are enhanced. Specifically: (1) The MD-H model is introduced to establish a transformation matrix with redundant parameters, completing the modeling of the transformation relationship of the industrial robot link coordinate system and laying the foundation for the solution of the kinematic inverse solution. (2) The COA and the kinematic inverse solution problem are combined, mapping the joint angles of the industrial robot to individual position coordinates and establishing an appropriate fitness function. (3) Elitism strategy and dynamic weighting factor are introduced to optimize the mutation, weight, and search mechanisms, improving the performance of COA. (4) The improved COA is applied to the kinematic inverse solution experiment of industrial robots, and the results show that the algorithm can effectively obtain the inverse solution.

The remaining parts of this paper are as follows. Section “Kinematic modeling of 6-DOF industrial robots” elaborates on the kinematic modeling principle of 6-DOF industrial robots; Section “A kinematic inverse solution method for the MD-H model based on the improved COA” elaborates on the specific principle of the kinematic inverse solution of industrial robots based on the improved COA; Section “Experiment” verifies the performance of the proposed method through experiments; finally, Section “Conclusion” provides a summary of the content of this paper.

Kinematic modeling of 6-DOF industrial robots

In the field of kinematic analysis of industrial robots, the D-H model, as a classic modeling method, is widely used. This model describes the pose transformation relationship between the base and the end effector of an industrial robot by defining four geometric parameters between adjacent links. However, when the axes of adjacent joints are parallel or nearly parallel, the D-H model may encounter kinematic singularity problems due to the deficiency of its parameter system: the coordinate system establishment process lacks continuity, resulting in sudden changes in some parameters. To address this issue, Hayati ³⁰ proposed an improved solution, which introduced an additional rotation angle $\beta$ around the y-axis to establish an extended kinematic model with redundant parameters. Besides Hayati, other scholars have successively proposed various alternative modeling schemes, including the complete and parametrically continuous (CPC) parameter model ³¹ and the product of exponentials (POE) formula. ³² Notably, the current mainstream controllers in the industrial field still generally develop algorithms based on the standard D-H model or its derivative models. Against this technical background, this paper adopts the MD-H modeling method to construct a complete kinematic model of a 6-DOF industrial robot.

The industrial robot selected for this study is the IRB2600 industrial robot from ABB, which is a 6-DOF serial robot. The schematic diagram of the robot’s link model is shown in Figure 1. Due to the parallel relationship between the second and third axes of this robot, this special structure poses limitations in traditional kinematic modeling. Therefore, an additional angle parameter $\beta$ is introduced in the D-H model to avoid singularity problems during kinematic analysis. In the MD-H model, to describe the spatial relationship between two adjacent links $i-1$ and $i$ (as shown in Figure 2), five key parameters are required, namely the joint angle ${\theta }_i$ , the link length $a_{i-1}$ , the link twist angle ${\alpha }_{i-1}$ , the link offset $d_i$ , and the aforementioned additional parameter $\beta$ . The homogeneous transformation matrix ${}_0{}^{i}T$ from the base coordinate system 0 to the joint coordinate system i can be obtained by left-multiplying the transformation matrices ${}_{i-1}{}^{i}T$ between each pair of adjacent links, that is:

$${}_0{}^{i}T={}_0{}^{1}T{}_1{}^{2}T\cdots {}_{i-1}{}^{i}T$$

(1)

Based on the MD-H model, the kinematic parameters of the IRB2600 industrial robot are shown in Table 1. The values in parentheses represent the parameters at the zero point state (initial values). Among them, ${\alpha }_{i-1}$ , $a_{i-1}$ , and $d_i$ are all fixed parameters, while ${\theta }_i$ is a variable parameter related to the end effector pose.

Figure 1.

ABB robot link model.

Figure 2.

Adjacent connecting rod model.

Table 1.

Kinematic parameters of IRB2600 industrial robot.

Link	${\alpha }_{i-1}$ (rad)	$a_{i-1}$ (mm)	${\theta }_i$ (rad)	$d_i$ (mm)	${\beta }_i$ (rad)
1	0	0	${\theta }_1(0)$	445
2	$-\pi /2$	150	${\theta }_2(-\pi /2)$	0	${\beta }_2(0)$
3	0	700	${\theta }_3(0)$	0
4	$-\pi /2$	115	${\theta }_4(0)$	795
5	$\pi /2$	0	${\theta }_5(0)$	0
6	$-\pi /2$	0	${\theta }_6(0)$	85

As long as each parameter in the table is known, the correspondence between the joint angles and the end-effector pose can be confirmed. Therefore, in the kinematic analysis of the robot, as long as the parameters are given, a unique transformation matrix ${}_0{}^{6}T$ can be obtained. All the above variables are modeled and simulated through the initial values using the Matlab Robotics Toolbox, and the robot pose is shown in Figure 3.

Figure 3.

Simulation diagram of IRB2600 model.

A kinematic inverse solution method for the MD-H model based on the improved COA

The principle of COA

COA initialization

COA is a population-based meta-heuristic method, where each raccoon individual is regarded as a member of the population, and its position in the search space represents the values of decision variables. Therefore, in COA, the position $X_i$ of a certain raccoon represents a candidate solution to the problem. At the start of COA, the positions of raccoon individuals need to be randomly initialized in the search space based on the following formula.

$$\begin{array}{c}\multicolumn{1}{c}{X_i:{\theta }_{i,j}={\theta }_j^{\min }+r({\theta }_j^{\max }-{\theta }_j^{\min })}\end{array}$$

(2)

Here, $i=1,2,\cdots ,N$ , where N represents the total number of raccoons; $=1,2,\cdots ,m$ , with m being the number of decision variables; ${\theta }_{i,j}$ is the value of the $\mathrm{jth}$ decision variable; $r$ is a random real number within the interval [0, 1]; ${\theta }_j^{\min }$ and ${\theta }_j^{\max }$ are the lower and upper limits of the $\mathrm{jth}$ decision variable, respectively.

The population in COA can be represented by the population matrix $X$ :

$$\begin{array}{c}\multicolumn{1}{c}{X={\left[\begin{array}{c}{X}_1\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{X_i}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{X_N}\end{array}\right]}_{N\times m}={\left[\begin{array}{ccccc}{\theta }_{1,1}\hfill & \cdots \hfill & {\theta }_{1,j}\hfill & \cdots \hfill & {\theta }_{1,m}\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill & .\hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\theta }_{i,1}}& \cdots \hfill & {\theta }_{i,j}\hfill & \cdots \hfill & {\theta }_{i,m}\hfill \\ \multicolumn{1}{c}{⋮}& .\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\theta }_{N,1}}& \cdots \hfill & {\theta }_{N,j}\hfill & \cdots \hfill & {\theta }_{N,m}\hfill \end{array}\right]}_{N\times m}}\end{array}$$

(3)

To assess the feasibility of each candidate solution, a fitness function needs to be constructed. The definition of the fitness function $F$ is as follows:

$$\begin{array}{c}\multicolumn{1}{c}{F={\left[\begin{array}{c}{F}_1\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F_i}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F_N}\end{array}\right]}_{N\times 1}={\left[\begin{array}{c}F\left(X_1\right)\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F\left(X_i\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F\left(X_N\right)}\end{array}\right]}_{N\times 1}}\end{array}$$

(4)

In the COA, $F$ is the fitness function vector, where $F_i$ is the fitness function value calculated based on the $\mathrm{ith}$ individual. In most meta-heuristic algorithms, the quality of potential solutions is evaluated by the fitness function value. Therefore, the individual in the population that can produce the best fitness function value is recognized as the optimal individual. As the algorithm iterates, potential solutions are constantly updated, and the optimal individual in the population is also updated accordingly.

The mathematical model of COA

In COA, the update process of an individual position $X_i$ is based on the simulation of two natural behaviors of raccoons. Therefore, the population update process in COA can be divided into two stages.

• 1. The hunting strategy when encountering prey (global search behavior).

When simulating the predation behavior, COA sets the position of the best individual in the population as the position of the prey. Meanwhile, it is assumed that half of the individuals approach the prey while the other half wait in place. For the individuals that actively approach the prey (with the index $i=1,2,\cdots ,\left[\frac{N}{2}\right]$ , where [·] represents the floor function), the position update formula is:

$$\begin{array}{c}\multicolumn{1}{c}{X_i^{P_1}:{\theta }_{i,j}^{P_1}={\theta }_{i,j}+r(\mathrm{pre}{y}_j-I·{\theta }_{i,j})}\end{array}$$

(5)

Here, $\mathrm{pre}{y}_j$ represents the value of the $\mathrm{jth}$ decision variable of the prey’s position, and $I$ is an integer randomly selected from the set $\{1,2\}$ . This formula simulates the process of the raccoon approaching its prey. When the prey escapes, it will be randomly placed at a new position ${\mathrm{prey}}_j^C$ in the search space:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{prey}}_j^C={\theta }_j^{\min }+r({\theta }_j^{\max }-{\theta }_j^{\min })}\end{array}$$

(6)

Based on this random position ${\mathrm{prey}}_j^C$ , the previously stationary raccoon individuals (with index $i=\left[\frac{N}{2}\right]+1,\cdots ,N$ ) will update their positions in the search space according to the following formula:

$$\begin{array}{c}\multicolumn{1}{c}{X_i^{P1}:{\theta }_{i,j}^{P1}=\{\begin{array}{cc}{\theta }_{i,j}+r·({\mathrm{prey}}_j^C-I·{\theta }_{i,j}),\hfill & F_{\mathrm{pre}{y}^C}<F_i\hfill \\ \multicolumn{1}{c}{{\theta }_{i,j}+r·({\theta }_{i,j}-{\mathrm{prey}}_j^C),}& \mathrm{others}.\hfill \end{array}}\end{array}$$

(7)

For each individual, if the new fitness function value $F_i^{P1}$ calculated for the new position $X_i^{P1}$ is less than the original fitness function value $F_i$ (indicating an improvement in the fitness function value), the update process is accepted and individual i moves to the new position $X_i^{P1}$ . Otherwise, individual i remains at its previous position $X_i$ . In this way, the algorithm ensures that it searches in the direction of a better solution.

$$\begin{array}{c}{X}_i=\{\begin{array}{c}{X}_i^{P1},F_i^{P1}<F_i\hfill \\ X_i,\mathrm{others}.\hfill \end{array}\hfill \end{array}$$

(8)

• 1. Escape behavior upon encountering a predator (local search behavior).

The position update strategy in the second stage mimics the escape behavior of raccoons when encountering predators. In this stage, a position $X_i^{P_2}$ is randomly generated near the current position of each raccoon, as defined by equations (9) and (10).

$$\begin{array}{c}\{\begin{array}{c}{\theta }_j^{\mathrm{ma}{x}_{\mathrm{near}}}=\frac{{\theta }_j^{\max }}{t}\hfill \\ {\theta }_j^{\mathrm{mi}{n}_{\mathrm{near}}}=\frac{{\theta }_j^{\min }}{t}\hfill \end{array}\hfill \end{array}$$

(9)

In the formula, ${\theta }_j^{\mathrm{ma}{x}_{\mathrm{near}}}$ , ${\theta }_j^{\mathrm{mi}{n}_{\mathrm{near}}}$ represent the local upper and lower limits of the jth decision variable at this position, $t=1,2,\cdots ,T$ , where $T$ is the upper limit of the number of iterations, and $t$ indicates the current iteration number. The individual position update formula is:

$$\begin{array}{c}{X}_i^{P_2}:{\theta }_{i,j}^{P_2}={\theta }_{i,j}+(1-2r)·\hfill \\ ({\theta }_j^{\mathrm{mi}{n}_{\mathrm{near}}}+r·({\theta }_j^{\mathrm{ma}{x}_{\mathrm{near}}}-{\theta }_j^{\mathrm{mi}{n}_{\mathrm{near}}}))\hfill \end{array}$$

(10)

As in stage 1, if the newly calculated position $X_i^{P2}$ can improve the fitness function value, the position is acceptable and the raccoon moves to the new position $X_i^{P_2}$ ; otherwise, the raccoon remains at the original position $X_i$ unchanged.

$$\begin{array}{c}{X}_i=\{\begin{array}{cc}{X}_i^{P2},\hfill & F_i^{P2}<F_i\hfill \\ \multicolumn{1}{c}{X_i,}& \mathrm{others}.\hfill \end{array}\hfill \end{array}$$

(11)

Improved COA based on elite retention strategy and dynamic weighting factor

The standard COA is prone to getting stuck in local optimal solutions, having low convergence efficiency, and struggling to effectively balance different types of errors when dealing with complex nonlinear optimization problems. These issues limit its application effectiveness in complex optimization scenarios. To address these problems, this study proposes an improved COA based on an elite retention strategy and dynamic weighting factors. By optimizing the mutation mechanism, weight strategy, and search strategy, the performance of the algorithm has been significantly enhanced.

• 1. Mutation mechanism based on elite retention strategy

The random variation of the standard COA only exists in the random generation process of the prey as described in equations (7) and (10). Its variation mode is single and the result is easy to be ignored, leading to insufficient population diversity and making the algorithm prone to getting stuck in local optimal solutions during the search process. Therefore, this paper’s improved algorithm adopts the mutation strategy of the elite retention selection method. Let ${\eta }_e$ represent the number of elite retained individuals, and $\mu$ be the elite retention rate, that is, the elite individuals are retained in proportion to $\mu$ of the population size $N$ . Then, the number of elite retained individuals is:

$$\begin{array}{c}{\eta }_e=\left[\mu N\right]\hfill \end{array}$$

(12)

Non-elite individuals undergo random mutations with a noise attenuation mechanism, thereby ensuring population diversity while reducing the damage to excellent solutions.

$$x_j^{new}=x_j+\left(1-\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{$T$}\right.\right)r\left({\theta }_j^{max}-{\theta }_j^{min}\right)$$

(13)

For the new position $X_i^{\mathrm{new}}$ calculated for the mutant individual, if the new fitness function value $F_i^{\mathrm{new}}$ is less than the original fitness function value $F_i$ , then the update process is acceptable and individual i moves to the new position $X_i^{\mathrm{new}}$ . Otherwise, individual i remains at its previous position $X_i$ unchanged.

$$\begin{array}{c}{X}_i=\{\begin{array}{cc}{X}_i^{\mathrm{new}},\hfill & F_i^{\mathrm{new}}<F_i\hfill \\ \multicolumn{1}{c}{X_i,}& \mathrm{others}.\hfill \end{array}\hfill \end{array}$$

(14)

Mutation can dynamically adjust the search range through $\left(1-\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{$T$}\right.\right)$ , enhancing global exploration in the early iterations to avoid local optima, and focusing on local search exploitation in the later stages to improve convergence.

• 1. Dynamic logarithmic weighting strategy

To enable the algorithm to handle different types of errors more evenly, the improved algorithm introduces dynamic logarithmic weights to adjust each type of error value. The error formula is as follows:

$$\begin{array}{c}{w}_k={\omega }_k(\log (\frac{\mathrm{er}{r}_k}{\sum \mathrm{er}{r}_i-\mathrm{er}{r}_k}+1)+1)\hfill \end{array}$$

(15)

In the formula, ${\omega }_k$ represents the fixed weight corresponding to the $\mathrm{kth}$ type of error, $w_k$ represents the dynamic weight that changes with the ratio of errors, and represents the value of the ith type of error. Therefore, the fitness function considering dynamic errors is as follows:

$$\begin{array}{c}{F}_i=\sum _{k=1}^n{w}_k\mathrm{er}{r}_k\hfill \end{array}$$

(16)

Equation (16) obtains $F_i$ by weighted summation of the dynamic adjustment weights $w_k$ and the error values of each type $\mathrm{er}{r}_k$ , which can reflect the comprehensive effect of the algorithm in handling various types of errors. However, during the iterative optimization process of the algorithm, with the introduction of dynamic weights, the dynamic weights $w_k$ in each round may change, which leads to a lack of direct comparability between the $w_k$ of different rounds. As a result, the fitness function value cannot directly reflect the convergence effect of the algorithm. To address this issue, this paper additionally introduces an evaluation fitness function for assessing the convergence effect of the algorithm:

$$\begin{array}{c}F{s}_i=\sum _{k=1}^n{\omega }_k\mathrm{er}{r}_k\hfill \end{array}$$

(17)

Improving the process of solving the kinematic inverse problem using COA

The MD-H model describes the spatial relationship of robot links through joint angles ${\theta }_i$ , link lengths $a_{i-1}$ , link twist angles ${\alpha }_{i-1}$ , link offsets $d_i$ , and additional parameter $\beta$ . When solving the inverse kinematics problem, the general known condition is the target pose of the end effector $T_d={[x_d,y_d,z_d,{\varphi }_d,{\theta }_d,{\psi }_d]}^T$ , where $x_d$ , $y_d$ , $z_d$ represent the Cartesian coordinates of the target position, and ${\varphi }_d$ , ${\theta }_d$ , ${\psi }_d$ represent the Euler angles of the target orientation. Based on this target pose, the corresponding joint angles ${\theta }_i$ need to be obtained. Therefore, if the six joint angles of the industrial robot are combined and mapped to the decision variables of a raccoon individual, the number of decision variables $m=6$ .

$$\begin{array}{c}{\theta }_i=[{\theta }_1^{\left(i\right)},{\theta }_2^{\left(i\right)},{\theta }_3^{\left(i\right)},{\theta }_4^{\left(i\right)},{\theta }_5^{\left(i\right)},{\theta }_6^{\left(i\right)}]\hfill \end{array}$$

(18)

Here, ${\theta }_j^{(i)}\in [{\theta }_i^{\min },{\theta }_j^{\max }]$ represents the $\mathrm{jth}$ joint angle of the $\mathrm{ith}$ individual. For the kinematic inverse solution problem of the MD-H model, the fitness function $F$ can be constructed based on the error between the actual pose $T_c={[x_c,y_c,z_c,{\varphi }_c,{\theta }_c,{\psi }_c]}^T$ of the industrial robot’s end effector and the target pose $T_d$ . In this paper, the combination of the position error $\mathrm{er}{r}_p$ and the orientation error $\mathrm{er}{r}_o$ between the two is selected as the fitness function, and the corresponding formula is as follows:

$$\begin{array}{c}\{\begin{array}{c}\mathrm{er}{r}_p=\sqrt{{(x_c-x_d)}^2+{(y_c-y_d)}^2+{(z_c-z_d)}^2}\hfill \\ \mathrm{er}{r}_o=\sqrt{{({\varphi }_c-{\varphi }_d)}^2+{({\theta }_c-{\theta }_d)}^2+{({\psi }_c-{\psi }_d)}^2}\hfill \end{array}\hfill \end{array}$$

(19)

$$\begin{array}{c}F=w_p\mathrm{er}{r}_p+{\mathrm{w}}_o\mathrm{er}{r}_o\hfill \end{array}$$

(20)

Among them, the actual pose $T_c$ is calculated based on the current joint angles (obtained through the forward kinematics of the MD-H model), corresponding to the Cartesian coordinates of the end effector. $w_p$ is the dynamic position weighting coefficient, and ${\mathrm{\omega }}_o$ is the dynamic orientation weighting coefficient.

COA process

After all raccoon individuals’ positions in the search space are updated according to the first and second stage rules described in Section “Kinematic modeling of 6-DOF industrial robots,” the mutation is carried out using elitism strategy, and one iteration is completed. Thereafter, the process of updating the population based on equations (5)–(14) will be repeated continuously until the preset maximum number of iterations is reached, and the final optimization result is obtained. To further enhance the algorithm’s efficiency and avoid unnecessary computational waste, an early stopping condition is introduced in COA in this paper. Specifically, after each iteration, it is determined whether the current optimal evaluation fitness value $F{s}_{\min }$ is less than the threshold $\epsilon$ . If it is less than the threshold, the iteration is stopped early. The algorithm flowchart is shown in Figure 4.

Figure 4.

Improved COA process.

Experiment

For the computational tasks and code execution in this study, the experiments were carried out on a computer featuring a 12th Gen Intel^® Core™ i7-12700H processor (2.30 GHz), 40.0 GB of RAM, running Windows 11 Home Chinese Edition (Version 24H2), with Python 3.10 as the programming environment. The ABB industrial robot model used in this article is IRB2600-20/1.65. According to the official documentation, its repeatability accuracy is 0.04 mm. The working range of all joints is shown in Table 2.

Table 2.

Working ranges of each joint of IRB2600-20/1.65.

Joint number	Working range	Joint number	Working range
1	$+{180}^{°}$ to $-{180}^{°}$	4	$+{400}^{°}$ to $-{400}^{°}$
2	$+{155}^{°}$ to $-{95}^{°}$	5	$+{120}^{°}$ to $-{120}^{°}$
3	$+{75}^{°}$ to $-{180}^{°}$	6	$+{400}^{°}$ to $-{400}^{°}$

The improved COA parameters set in this paper are shown in Table 3.

Table 3.

COA parameter settings.

Parameter	Value	Parameter	Value
N	200	${\mathrm{\omega }}_o$	30
T	200	$\mu$	$0.2$
${\mathrm{\omega }}_p$	1	$\epsilon$	0.005

The random motion of industrial robots was controlled, and the Cartesian coordinates and Euler angles corresponding to any 10 points were read and recorded. The improved COA was used to calculate the kinematic inverse solution corresponding to each point, and the position error and orientation error corresponding to each set of inverse solutions were calculated. The results show that the improved COA can accurately obtain the kinematic inverse solution of industrial robots. The inverse solution error obtained by the algorithm is small, with position errors all less than 0.003 mm and orientation errors all less than 0.001 rad. These errors are within the allowable error range of the IRB2600 industrial robot. The inverse solution error obtained by this algorithm is small, with position errors all less than 0.005 mm and orientation errors all less than 0.005 rad. These errors are within the allowable error range of the IRB2600 industrial robot and fall within the normal error range of its operation.

For comparison, the algorithm proposed in this paper is compared with WOA, standard COA, PSO, and GA through experiments. The population size, maximum number of iterations, and evaluation function $F_s$ of the five algorithms are set to be the same. The inverse kinematics solutions for 10 points were solved. To reduce the influence of random factors, each point was calculated 10 times. The average results of each algorithm are shown in Table 4. It can be seen from the table that the improved COA is significantly superior to the comparison algorithms in terms of the number of iterations, fitness, and error indicators.

Table 4.

Results of optimization algorithm.

Algorithm	$\mathrm{average}\mathrm{iterations}$	$\mathrm{Fs}$	$\mathrm{er}{r}_p$ (mm)	$\mathrm{er}{r}_o$ (rad)
Improved COA	40.4	0.0046	0.0021	$8.5\times {10}^{-5}$
WOA	200	69.2211	15.7728	2.0620
COA	200	13.1753	8.6206	0.4291
PSO	136.2	17.2554	0.0011	0.5751
GA	200	25.7447	0.0015	0.8581

To visually present the performance differences among the five optimization algorithms, the average convergence curves of the algorithms were further plotted. Considering that the improved COA has relatively fewer iterations, after the algorithm’s iteration stops, the minimum $F_s$ it achieved was filled up to $T=200$ times. The results are shown in Figure 5. As can be seen from the figure, the fitness of the improved COA converges rapidly after the iteration begins, which is consistent with the results of the number of iterations in Table 4. In contrast, although the other four optimization algorithms also converge very quickly, they subsequently trapped in local optima until the maximum number of iterations are reached. In summary, the fitness value of this algorithm is smaller than that of the WOA, COA, PSO, and GA optimization algorithms, and it converges faster with a more stable curve. The final converged value of the improved COA is less than the threshold and much lower than the final stable values of other algorithms, demonstrating higher optimization accuracy. Overall, the improved COA outperforms the standard COA and other classical intelligent optimization algorithms in terms of convergence rate and optimization accuracy, verifying the effective enhancement of the algorithm’s performance through the improved strategy.

Figure 5.

Comparison of iteration effects of five optimization algorithms.

To further verify the robustness of the improved COA algorithm, this paper plots the detailed convergence curves of five algorithms at 10 different spatial points (Figure 6), where each point position has undergone 10 independent calculations, and the results are averaged. This figure visually demonstrates the convergence behavior patterns of each algorithm when facing different inverse kinematics solution scenarios. From the analysis of Figure 6, it can be seen that the improved COA shows the characteristics of rapid initial convergence and stable convergence in the later stage at all test points. It can quickly reduce the fitness value within the first 20 iterations, thanks to the synergy of the elite retention strategy and the dynamic weight mechanism, enabling the algorithm to identify and search toward the optimal solution, and eventually converge below the threshold to meet the truncation condition. In contrast, the other algorithms, although converging rapidly in the initial stage, are prone to getting stuck in local optima, showing convergence stagnation at some test points and unstable final convergence accuracy.

Figure 6.

The iterative curves of five algorithms at 10 points: (a) Improved COA, (b) WOA, (c) COA, (d) PSO, and (e) GA.

Conclusion

This paper focuses on the inverse kinematics problem of 6-DOF industrial robots, taking the ABB IRB2600 industrial robot as the research object. An MD-H kinematic model is constructed and an inverse kinematics solution method based on the improved COA is proposed. By introducing the MD-H model, a transformation matrix with redundant parameters is established to complete the modeling of the transformation relationship of the industrial robot’s link coordinate system. The COA is improved by introducing an elite retention strategy and dynamic weight factors, optimizing the mutation, weight, and search mechanism to enhance the algorithm performance. Experimental results show that the improved COA can effectively obtain the inverse solution, with a position error less than 0.003 mm and an orientation error less than 0.001 rad. In terms of performance indicators such as the number of iterations and the optimal evaluation fitness, it significantly outperforms comparison algorithms such as WOA, standard COA, PSO, and GA, providing a new intelligent algorithm solution for the inverse kinematics problem of 6-DOF industrial robots.

However, it must be admitted that the current algorithm version has limitations in terms of computational efficiency. The actual time consumption is higher than the expected requirements of the application. This is mainly due to the complexity of the iterative update rules, the high dimensionality of the intermediate variables, and the computational burden brought by the comprehensive evaluation mechanism. In actual experiments, when individuals in the population are provided with appropriate initial values, the algorithm can achieve satisfactory time efficiency under specific conditions. However, this approach still presents certain limitations in terms of broader applicability, and further optimization of the algorithm’s computational framework will be necessary to enhance its universal performance across diverse operational scenarios. Therefore, in subsequent improvements, we can focus on the multi-algorithm fusion strategy and integrate an optimization algorithm that can quickly obtain initial values. Besides the fusion of algorithms, in the future, we will further optimize the algorithm structure by simplifying redundant iterations and streamlining the population update mechanism to reduce unnecessary fitness evaluations.