Multiobjective optimization of collaborative robotic task sequence assignment problems under collision-free constraints

Abstract

This paper proposes a multiobjective optimization approach to address the challenge of collaborative manufacturing with multiple robot arms. Given the necessity for multiple robot arms to work together, the potential for collisions between robotic motions is a significant concern, and the automated task sequence assignment for robots becomes increasingly complex. Previous research has either simplified the collision-free conditions in a limited working area, or employed a master-slave approach to obtain only a local solution. Consequently, we propose a unified global optimization approach for simultaneously addressing various collaborative manufacturing issues, including robotic task sequence assignment (RTSA), multiple inverse kinematics (IK) selection, joint-space collision-free operations and multiple manufacturing objectives. As the optimal collaborative RTSA problem is a combinatorial optimization problem with non-deterministic polynomial-time hard (NP-hard) complexity, this paper presents a hybrid nondominated sorting genetic algorithm III (NSGA-III) method that integrates a Hamming-distance-based method and a greedy strategy within NSGA-III to improve population diversity and solution quality. To validate the efficacy of the proposed approach, simulation experiments were conducted on cooperative manufacturing scenarios, with two objectives: task completion time and task load balancing. The experimental results demonstrate that the proposed approach is effective in obtaining collision-free Pareto solutions. Furthermore, the proposed hybrid NSGA-III method obtains superior solutions compared to the original method in the studied problem, as measured by two performance indices.

Introduction

Collaborative robot arms are typically employed to boost automated production efficiency in the manufacturing industry. However, when multiple arms are operating simultaneously, the potential for collisions between their trajectories is significant and critical. ¹ The traditional design method relies on the experience of engineers, requiring significant investment of time and resources to construct the robotic task sequences and adjust the collision-free motion of robots through trial and error. Therefore, intelligent automation of collaborative robotic task sequence assignment (CRTSA) offers significant benefits. ² However, the process of automated searching for optimal task sequences and generating collision-free robotic trajectories is complex and challenging.

To obtain a solution to the CRTSA, it is necessary to discuss the problem of a large number of combinations between task sequences. If the problem is simply to split the multiple points into a number of independent point-to-point problems, this scheme will not be able to obtain the truly optimal solution for the overall operational efficiency. ³ Conversely, when the task sequences are considered in the manufacturing planning, each arm must begin at the starting point and visit all the intermediate points before returning to the original starting point to resume the manufacturing process in the next round. This is equivalent to solving the multi-travel salesman problem (MTSP). MTSP is a generalization of the single travel-salesman problem (TSP), ⁴ which considers multiple Hamiltonian circuits covering all the points to be visited. The problem is non-deterministic polynomial-time hard (NP-hard), and there is no polynomial-time algorithm for obtaining a truly optimal solution, making it highly challenging. In addition, at each task point on the task sequence, it is necessary to select the most suitable inverse kinematics (IK) solutions for the robot arms, which further increases the complexity of the problem.^2,5 Furthermore, multi-arm cooperation in the real world typically involves more than one objective. For instance, improving production efficiency often entails reducing the completion time and balancing the task loading of each arm. This avoids overusing some arms and ensures an optimal overall yield. ⁶ Therefore, in this paper, the CRTSA problem is studied and formulated into a multiobjective optimization problem (MOP).

Due to the high complexity of the robotics problem, several previous studies have developed the problem towards a single objective paradigm and used an evolutionary algorithm (EA) as a solution method. For instance, the modified differential evolution method in Pierezan et al. ⁷ and Weihmann et al. ⁸ and the beetle antennae olfactory (BAO) algorithm in Tan et al. ⁹ However, when the problem under investigation requires the consideration of multiple conflicting objectives, the resulting solution will not be a single optimal solution. Rather, it will be a set of trade-off solutions, a Pareto set, among the objectives. In such cases, a multiobjective optimization evolutionary algorithm (MOEA) has emerged as a popular method for solving MOP problems. ¹⁰

In general, MOEA can be categorized into three types: Pareto-based MOEA, indicator-based MOEA and decomposition-based MOEA. Classical Pareto-based MOEAs propose a variety of methods for Pareto-ranking ¹¹ in the objective space. These include the domination count in the multiobjective optimization genetic algorithm (MOGA), ¹² the strength-Pareto values and external archive set for Pareto-optimal solutions in the Strength Pareto Evolutionary Algorithm (SPEA), ¹³ a further improvement of SPEA with density evaluation in SPEA-II, ¹⁴ and the non-domination sort combined with a distance-crowding sort in the nondominated sorting genetic algorithm II (NSGA-II) method. ¹⁵ However, Pareto-based MOEAs are prone to generate a considerable number of non-dominated but similar solutions in the population at a later stage of evolution. This may result in poor convergence or only a partial coverage of the Pareto front. This phenomenon is particularly evident when solving problems with more than three objectives. Indicator-based MOEAs employ performance metrics for quality measurements to define genetic selection. ¹⁶ The hypervolume (HV) is the most representative metric and has been utilized in studies such as the indicator-based evolutionary algorithm (IBEA) ¹⁷ and the S metric selection evolutionary multiobjective algorithm (SMS-EMOA). ¹⁸ However, the scalability of the hypervolume is constrained by its high computational cost. Consequently, approximate methods are typically employed, such as the hypervolume estimation algorithm for multiobjective optimization (HypE). ¹⁹ Decomposition-based MOEAs address the convergence problem by decomposing the original multiobjective optimization problem into multiple single-objective optimization problems. This approach has been employed in the multiobjective evolutionary algorithm based on decomposition (MOEA/D), ²⁰ which was the first to propose this decomposition technique. Another notable example is the non-dominated sorting genetic algorithm III (NSGA-III), which generates uniform reference points as the convergence direction based on the objective hyperplane. In recent years, decomposition-based MOEAs have gained considerable popularity. ¹¹

There exist a multitude of studies in the literature which employ the MOEA as a solution for robotics-related optimization problems. In these studies, NSGA-II and NSGA-III are frequently employed methodologies, as evidenced by their applications to the optimization of energy consumption, weld quality and cable twist fluctuations in robotic welding, ²¹ the optimization of the motion planning of a 7-DOF assembly robot, ²² the optimization of proportional–integral–derivative (PID) controller parameters to minimize manipulated positional error and torque variations in a robotic manipulator of two-degree-of-freedom, ²³ and the optimization of the payload capacity and radial range in an in-pipe inspection robot. ²⁴ As previously demonstrated in a study, ²⁵ NSGA-III yields superior solution quality than NSGA-II, even in the case of a small number of objectives (in this instance, three objectives). Consequently, NSGA-III is employed as the MOP solution in this study. However, as the problem under investigation belongs to the class of discrete multiobjective optimization problems, the NSGA-III algorithm is susceptible to premature convergence. Consequently, we have developed a hybrid approach that combines the NSGA-III with a Hamming distance-based method, with the aim of enhancing the overall solution capability.

Furthermore, the aforementioned single-objective or multiobjective robotics problem only considers the application of a single robot. Given that the design of each individual robot already precludes the possibility of mutual interference between the robot rods, the collision between robot rods is not a concern. Nevertheless, when investigating the potential for collaboration between multiple industrial robots, it is essential to consider the spatial collision problem between the rods of different robots, given the overlap of their operating workspaces. Previous studies have examined the co-operation of multiple robot arms, ²⁶ with a particular focus on the analysis of gripping applications ²⁷ and the characteristics of the stress and force load distribution. ²⁸ However, there has been less discussion on the problem of searching collision-free trajectories in collaborative multi-joint robot arms.

In the literature, the main studies focused on the collision avoidance problem can be categorized as online and offline. ²⁹ Online collision avoidance problem has the advantage of real-time collision handling. ³⁰ Most of the solutions are implemented in the actual online processing, and do not need to set up the collision avoidance environment configuration and commands in advance. Its strategy is mainly aimed at the dynamic detection of collision-free point-to-point operation during the local operation, and to complement or modify the original path planning if collisions may occur in the subsequent operations. ³¹ However, since it focuses on the real-time collision in a local region, global collision-free planning is usually not considered. Conversely, the offline collision-free planning problem, in which all the settings and commands for the robot arms are planned before manufacturing begins, can deal with the global planning of more complex manufacturing tasks. This paper addresses the problem of offline collision avoidance optimization.

Due to the high complexity of 3D joint space collision avoidance, ² there is only a few academic research on global collision-free multi-robot cooperation. Chang et al. ³² used geometric modelling to calculate the collision avoidance path planning for two robot arms, which fitted the robot arms as polyhedrons according to the mechanical structure, and generated a collision map according to the paths. However, their study fixed the paths of one of the robot arms, and then adjusted the paths of the other robot arm, which limited the flexibility of adopting the dual-arms. Additionally, it was challenging to consider the paths for global optimization. Ćurković et al. ³³ used a cooperative evolutionary computational approach to investigate the collision-free paths of the two arms. However, they simplified the arm to linkage segments as the geometrical inference targets, and only focused on collision-free planning between single actions; Larsen et al. ³⁴ compared two obstacle avoidance planning schemes: a heuristic method combined with a sampling-based algorithm and an evolutionary algorithm. The former demonstrated good performance, but its applicability was limited to a single case. Redesign was required for different cases. The latter had greater generality, but a master-slave approach was used to plan the paths of the two robots. The main robot generates paths randomly. Subsequently, the other robot is paired with it, and after a fixed number of rounds of testing, the paths are re-planned if the path pairing fails. However, the master-slave planning is still a local solution. Additionally, the number of trials is usually difficult to set for more complex cases. Ni et al. ³ used a particle swarm optimization method to study the coordination of two-armed robots. Their approach avoids collisions with obstacles by planning the motion of the joints. However, it does not consider the self-collision-free requirements between the operations of two-armed robots. That is, the cooperation area of the two manipulators must not overlap on the time axis. Riboli et al. ³⁵ studied a motion planning approach for two-arm kinematics considering self-collision-free motions. However, their work is limited to solving the problem of collision-free motion planning given start and target points. It focuses mainly on 2D collision-free constraint solutions limited to their selected planar working area and spline-based motion trajectories. In this approach, 3D collision constraints are projected onto the 2D plane to simplify the derivation of geometric and timing relationships. In order to address the shortcomings of the aforementioned studies, our proposed hybrid NSGA-III solution method has integrated geometric inference and a greedy strategy into genetic flow. This framework resolves the critical issues that have been identified.

In summary, the primary contributions of our work are as follows: (1) Unlike the master-slave approach, which is limited to a local solution, the proposed multiobjective optimization approach offers more flexibility in planning the global CRTSA problem. This allows us to consider the efficiency and equilibrium workload of the multiple robot arms. (2) By integrating geometric inference with evolutionary computation, we have developed a multi-robot-arm spatial collision-free solution. This approach extends previous simplified collision-detection methods that either used the projected 2D geometric reasoning or were limited to single actions. (3) The proposed Hamming-distance-based NSGA-III has been validated to provide superior solution quality over the original NSGA-III in the CRTSA problem, as measured by performance measurement indices (PMIs).

The rest of the paper is organized as follows. Section ‘Problem definition and formulation’ provides an in-depth analysis of the forward and inverse kinematics of serial-link robots, defines the CRTSA problem in detail and transforms it into a multiobjective combinational optimization problem. Section ‘Solution using Hamming distance-based NSGA-III’ describes the proposed Hamming Distance-based NSGA-III as the optimization solver. Section ‘Experiments and discussions’ sets up and runs the simulated experiments and discusses the comparative results based on two performance metrics. The final section presents the conclusions of this work.

Problem definition and formulation

A multi-joint robot arm can theoretically have multiple IK solutions. This allows for evaluation of the robot’s kinematics model to obtain a variety of trajectory planning schemes that consider different IKs in a range of usage scenarios. ³⁶

Robotic forward kinematics

The forward kinematics (FK) of a serial-link robot arm is typically represented by a Denavit–Hartenberg (DH) formulation, ³⁷ which is a systematic approach to describe the FK by choosing the appropriate coordinate system for each joint and deriving the necessary parameters including $\alpha ,a,d$ and $\theta$ . To illustrate, consider a robot with six rotational joints ${\theta }_1\dots {\theta }_6$ as an example, ${{}^{i-1}A}_i$ is a homogeneous matrix which is used to represent the spatial transform matrix of joint i with respect to joint i−1 as follows:

$${{}^{i-1}A}_i=\left[\begin{array}{cccc}\cos \left({\theta }_i\right)\hfill & \multicolumn{1}{c}{-\sin \left({\theta }_i\right)}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{a_{i-1}}\\ \multicolumn{1}{c}{\sin \left({\theta }_i\right)\cos \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{\cos \left({\theta }_i\right)\cos \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{-\sin \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{-d_i\sin \left({\alpha }_{i-1}\right)}\\ \multicolumn{1}{c}{\sin \left({\theta }_i\right)\sin \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{\cos \left({\theta }_i\right)\sin \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{\cos \left({\alpha }_{i-1}\right)}& \multicolumn{1}{c}{d_i\cos \left({\alpha }_{i-1}\right)}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{1}\end{array}\right].$$

(1)

The matrix multiplication method can be used to obtain a coordinate transformation between the end of the manipulator and the base (or the endpoints of the joints) can be obtained. This results in a homogeneous matrix of the following form:

$${{}^{0}A}_6={{}^{0}A}_1{{}^{1}A}_2{{}^{2}A}_3{{}^{3}A}_4{{}^{4}A}_5{{}^{5}A}_6.$$

(2)

The derivation equations can be applied to all serial-link robots with only rotational or prismatic joints. The procedure for deriving the details, including the selection of coordinate axes and the determination of parameters, can be found in Craig. ³⁷ This paper examines two industrial robots, the PUMA 560 and the FANUC ER-4iA. The DH parameters are presented in Table 1.

Table 1.

PUMA 560 and FANUC ER-4iA DH forms.

(a) PUMA560						(b) FANUC ER-4iA
i	${\alpha }_{i-1}$ (deg)	$a_{\mathrm{i}-1}$ (mm)	$d_i$ (mm)	${\theta }_i$ (deg)	Range of ${\theta }_i$ (deg)	i	${\mathbf{\alpha }}_{\mathcal{i}-1}$ (deg)	$a_{\mathcal{i}-1}$ (mm)	${\mathbf{d}}_{\mathcal{i}}$ (mm)	${\mathbf{\theta }}_{\mathcal{i}}$ (deg)	Range of ${\mathbf{\theta }}_{\mathcal{i}}$ (deg)
1	0	0	0	${\mathrm{\theta }}_1$	[−160, 160]	1	0	0	0	${\theta }_1$	[−170, 170]
2	−90	0	0	${\mathrm{\theta }}_2$	[−225, 45]	2	−90	0	0	${\theta }_2$	[−110, 120]
3	0	431.8	149.09	${\mathrm{\theta }}_3$	[−45, 225]	3	0	260	0	${\theta }_3$	[−205, 69]
4	−90	20.32	433.07	${\mathrm{\theta }}_4$	[−140, 140]	4	−90	20	290	${\theta }_4$	[−190, 190]
5	90	0	0	${\mathrm{\theta }}_5$	[−100, 100]	5	90	0	0	${\theta }_5$	[−120, 120]
6	−90	0	0	${\mathrm{\theta }}_6$	[−266, 266]	6	−90	0	0	${\theta }_6$	[−360, 360]

Robotic inverse kinematics (IK)

However, since industrial robots operate in joint space, however, the practical manufacturing problem is defined in world coordinates (or Cartesian coordinate system). Therefore, the IK solution needs to be solved in reverse for the manipulation requirements in world coordinates. The IK formulation of PUMA560 has been previously reported in the literature, for example, in Craig. ³⁷ However, since no universal solution exists for the IK of a six-joint robot arm, ³⁸ and the IK formulation of the FANUC ER-4iA has not been fully discussed in academia in the past, we must derive it first.

The last three axes of an industrial robot structure are typically designed as wrist joints, aligning with the Pieper criteria. ³⁷ For specific end-effector poses, the last three joints can be used to determine the direction, in conjunction with the first three joints to determine the origin position of the wrist joints. This allows for the complete IK formulation to be algebraically deduced, resulting in a formula with a multiple-solution form. Suppose the pose of the end effector is as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Pose}=\left[\begin{array}{cccc}\multicolumn{1}{c}{}& R_{3\times 3}\hfill & \multicolumn{1}{c}{}& \multicolumn{1}{c}{T_{1\times 3}}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{1}\end{array}\right]}\\ \multicolumn{1}{c}{=\left[\begin{array}{cccc}{R}_{\mathrm{nx}}\hfill & \multicolumn{1}{c}{R_{\mathrm{ox}}}& \multicolumn{1}{c}{R_{\mathrm{ax}}}& \multicolumn{1}{c}{T_x}\\ \multicolumn{1}{c}{R_{\mathrm{ny}}}& \multicolumn{1}{c}{R_{\mathrm{oy}}}& \multicolumn{1}{c}{R_{\mathrm{ay}}}& \multicolumn{1}{c}{T_y}\\ \multicolumn{1}{c}{R_{\mathrm{nz}}}& \multicolumn{1}{c}{R_{\mathrm{oz}}}& \multicolumn{1}{c}{R_{\mathrm{az}}}& \multicolumn{1}{c}{T_z}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{1}\end{array}\right]={{}^{0}A}_1{{}^{1}A}_2{{}^{2}A}_3{{}^{3}A}_4{{}^{4}A}_5{{}^{5}A}_6,}\end{array}$$

(3)

where $R_{3\times 3}$ is the rotation matrix and $T_{1\times 3}$ is the translation vector. When the first three axes ${\theta }_1,{\theta }_2,{\theta }_3$ determine the origin position of the wrist, it can be formulated as follows:

$$\left[\begin{array}{c}{T}_x\hfill \\ \multicolumn{1}{c}{T_y}\\ \multicolumn{1}{c}{T_z}\\ \multicolumn{1}{c}{1}\end{array}\right]{=}^0{A}_1{·}^1{A}_2{·}^2{A}_3·\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{d_4}\\ \multicolumn{1}{c}{1}\end{array}\right].$$

(4)

Substituting the ${{}^{0}A}_1$ , ${{}^{1}A}_2$ and ${{}^{2}A}_3$ (calculated from the DH form) into the above equation, we can solve ${\theta }_1$ , ${\theta }_2$ and ${\theta }_3$ . The detailed solving procedure is shown in the appendix.

Next, since the last three axes ( ${\theta }_4,{\theta }_5,{\theta }_6$ ) determine the direction of the end of the manipulator, we multiply the left and right side of the equation (3) by the inverse matrix ${{{}^{0}A}_3}^{-1}$ . This results in the following equation:

$${{{}^{0}A}_3}^{-1}·\left[\begin{array}{cccc}{R}_{\mathrm{nx}}\hfill & \multicolumn{1}{c}{R_{\mathrm{ox}}}& \multicolumn{1}{c}{R_{\mathrm{ax}}}& \multicolumn{1}{c}{T_x}\\ \multicolumn{1}{c}{R_{\mathrm{ny}}}& \multicolumn{1}{c}{R_{\mathrm{oy}}}& \multicolumn{1}{c}{R_{\mathrm{ay}}}& \multicolumn{1}{c}{T_y}\\ \multicolumn{1}{c}{R_{\mathrm{nz}}}& \multicolumn{1}{c}{R_{\mathrm{oz}}}& \multicolumn{1}{c}{R_{\mathrm{az}}}& \multicolumn{1}{c}{T_z}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{1}\end{array}\right]{=}^3{A_4}^4{A_5}^5{A}_6,$$

(5)

Similarly, ( ${\theta }_4,{\theta }_5,{\theta }_6)$ can be found by substituting ${(}^3{A}_4{,}^4{A}_5{,}^5{A}_6)$ and ${{{}^{0}A}_3}^{-1}$ into the above equation. The solution procedure in this part is identical to that of PUMA 560, which can be found in Craig. ³⁷

It is worth noting that, based on the structure of the FANUC ER-4iA, up to eight sets of inverse kinematic solutions can be derived. However, some of these solutions may not be feasible solutions due to the limitation of the range of joint motion.

Task sequence assignment with multiple IK solutions

In a manufacturing task, the specified poses at each point may correspond to multiple IKs, which increases the flexibility of the manipulation. However, this also increases the complexity of the optimal path selection. Furthermore, if there are a large number of points to be manipulated, the overall manufacturing efficiency can be considerably affected by different visiting sequences and different selections of IKs. As an illustration, Figure 1(a) represents a sequence plan from left to right, with a relatively longer travel distance than that of Figure 1(b).

Figure 1.

Task assignment with different IKs and task sequences: (a) order of visits: A-B-C-D-E-F-G-H-J-A and (b) order of visits: A-B-C-E-G-D-F-H-I-J-A.

As previously stated, optimizing the sequence order is a NP-hard combinatorial optimization problem consisting of multiple IKs and distinct path sequences. The exhaustive search method is not a viable option for medium and large-scale manufacturing cases due to the size of the solution space. Consequently, heuristic algorithms are the commonly used option, and a comparison of solutions. This paper will discuss the comparison of solutions used in this context in subsequent chapters.

Collision-free constraints

When an industrial robot operates individually, the operation range of the joints has been appropriately limited in its original design to prevent self-interference. However, when different arms operate in cooperation, arm-to-arm collisions can occur.

In order to detect a collision, a larger bounding box, as illustrated in Figure 2, can be used to encompass a robotic link as the basis for collision detection. This approach offers two advantages. Firstly, it simplifies the free-form surfaces of the robot arm exterior into simple cuboids, which makes it easier to detect spatial intersection relationships by geometric inference. Secondly, although measurement errors may occur when simple cuboids are used instead of free-form shapes, this error can be used as an additional redundancy for collision avoidance, enhancing the safety of the operation in practical applications. ³⁹ Consequently, in the proposed scheme, the collision detection between any two arms will be converted to a spatial intersection problem between bounding boxes. The intersection of two boxes is a typical geometric problem. For further information, please refer to Schneider and Eberly. ⁴⁰ In order to complete this process, it is necessary to extend the collision calculation to each robotic component. The algorithm is shown in Algorithm 1.

Figure 2.

Bounding boxes for robotic collision detection: (a) a bounding box example and (b) a collision example in two bounding boxes.

Algorithm 1.

Collision_Detection(A, B): Integer.

Input:  A, B: Two multi-joint robots. Output:  Collision count. Begin  For each i in [1, (the link number of A)]   For each j in [1, (the link number of B)]    If collisions are detected between the bound boxes of the i-th link of A and the j-th link of B Then Increase the collision count by 1.    End   Next j  Next i  Return the collision count. End

Furthermore, since any two robot arms cannot interfere with each other throughout the timeline of cooperation, it is difficult to use a mathematical model to accurately represent the trajectories of each component of the robot arm scanned in 3D space during the motion. Consequently, this paper uses the timeline sampling method to divide the planned trajectory into discrete points and to ascertain that no collisions occur in these sampling points. The collision detection method described above requires the input of customized parameters based on the specific practical setup. Further details on the implementation method and algorithm will be provided in subsequent experiments.

The optimization problem

As previously stated, the CRTSA problem is to identify a collision-free task sequence assignment for robots that meets the two key objectives: task completion time and task load balancing. To model the optimization problem, some symbols are defined in Table 2. In general, the tasks are divided into sub-tasks and assigned for n robots. For the task assignment $R$ , we define the two-objective optimization problem:

$$\mathrm{Minimize}({\varphi }_1,{\varphi }_2)$$

$$\mathrm{subject}\mathrm{to}\mathrm{\varphi }(\mathrm{R},\mathrm{\delta })=0,$$

(6)

where ${\varphi }_1$  = Max{TR₁…TR_n}, ${\varphi }_2$  = Var(TR₁…TR_n), TR_i =  $\sum _{j=1}^{\left|r_i\right|-1}Ma{x}_{k=1}^6\frac{\left|r_{i.j+1}.{\theta }_j-r_{i.j}.{\theta }_k\right|}{{\stackrel{.}{\theta }}_k}$ , Var(TR₁…TR_n) = $\frac{1}{n}\sum _{k=1}^n(T{R}_k-\mu {)}^2$ , $\mu =\frac{1}{n}\sum _{k=1}^{n}T{R}_k$ , ${\stackrel{·}{\theta }}_k$ is the average angle velocity of ${\theta }_k$ , $\mathrm{\phi }(·)$ is the function evaluating the total collision count during performing the tasks. In particular, since $\mathrm{\phi }(·)$ is difficult to model mathematically, we have designed algorithms to calculate the complete amount of collisions based on Collision_Detection in Algorithm 1. This will be further discussed in the subsequent sections.

Table 2.

Symbols used in the CRTSA modelling.

$m$ :	the max task number in a sub-task;
$n$ :	the number of robots;
$l$ :	the number of IKs of each task point;
$\mathrm{R}$ :	task sequence assignment matrix $\mathrm{R}=\{r_{i,j}|1\le i\le n,1\le j\le m\}$ ; $r_{i,j}$ : the j-th task (point) assigned for i-th robot; $r_{i,j}.\mathrm{SelectedIK}$ : the selected IK index to multiple IK solutions of $r_{i,j}$ ; $r_{i,j}.{\mathrm{\theta }}_1\cdots r_{i,j}.{\mathrm{\theta }}_6$ : the joint values in the selected IK; $\left|r_i\right|$ : the assigned sub-task number for the i-th robot, $\left|r_i\right|\le m$ ; $R_i$ : the sub-task set assigned to the i-th robot, ${\mathrm{R}}_i=\{r_{i,1}\cdots r_{i,\left|r_i\right|}\}$ ;
$\mathrm{\delta }$ :	the timing resolution to detect the collision.

In addition to the collision constraints that must be satisfied first, the operating time of individual arms is evaluated by dividing the maximum joint change between consecutive poses by the average angular velocity. If the completion times of each arm are the same, a standard deviation of 0 represents the optimal equilibrium solution. Conversely, the larger the standard deviation, the greater the inconsistency in the completion times of each arm. However, the difficulty in solving this problem is that, in addition to optimally exploring the Pareto set for the dual-objective problem, it is necessary to identify the optimal task sequence assignment and to explore suitable IKs to generate collision-free trajectories. Given the highly constrained nature of the problem, a multiobjective evolutionary computation methodology is developed as the solution in the paper.

Solution using Hamming distance-based NSGA-III

This paper uses the NSGA-III scheme ¹⁸ to optimize this CRTSA problem. The overall architecture of NSGA-III is similar to that of NSGA-II. The main difference is that the last non-dominated set is now ranked using the reference point mechanism in the hyperplane of the objective vectors, rather than the original crowded distance sorting. This has been shown to be an effective approach in benchmark problems. ²⁷ However, this method is still primary focused on the domination in the objective space to maintain population diversity, rather than on ensuring genetic similarity within the population. When there is a high degree of population similarity, it is easy to fall into a local solution because of early maturation. We therefore propose a hybrid solution method to address this issue.

Figure 3 illustrates the overall flow. The initial population is randomly generated and the fittest individuals are selected for the mating pool. Crossover and mutation are then applied to obtain offspring according to the determined probability. The new generated offspring undergoes the fitness evaluations. After decoding, the overall task points are assigned to each arm in the order indicated by the gene code. If any arm is unable to reach the assigned task point, the repair strategy is activated. To ensure that the last assigned plan is available, the inapplicable task points are assigned to the other arm. A greedy search is then developed to select the most suitable IK for the task points assigned to robots. Finally, the objectives in equation (6) are calculated.

Figure 3.

Flowchart of proposed Hamming-distance-based NSGA-III.

Before entering the next generation, the offspring are merged with the parent to perform the NSGA-III non-domination sort and reference points ranking operation. Furthermore, in order to maintain the diversity of the population chromosomes, we have integrated the Hamming-distance-based method to preserve the diverse populations into the next generation, with the aim of improving the convergence quality.

Population initialization and coding schema

The proposed solution process involves calculating the IK solution for each task point firstly, with the multiple poses corresponding to each point then stored in the database. This allows for a quick query to be made when entering into the evolution process, avoiding the time-consuming and repeated solving of IK.

Chromosome coding is order-based coding, ⁴¹ and its genotype contains two elements: task point identifier and a task separation. For a case with n robot arms, each chromosome encoding will contain all task points and n−1 task separations. In the initialization phase, the order of the chromosomes is randomly shuffled. When decoding, the task points are assigned to each robot arm according to the task separation position. Figure 4 illustrates a two-arm (n = 2) cooperation example. The task order for arms 1 and 2 is (3, 6) and (2, 4, 1, 5), respectively.

Figure 4.

Chromosome coding schema.

Genetic selection, crossover and mutation

The algorithm employs customized genetic operators. First, the tournament selection is used to introduce some selection pressure (Figure 5). Once the superior parent individuals are selected for the mating pool, the partially mapped crossover (PMX) ⁴¹ and inverse mutation operators are applied according to the probability, respectively. The inverse mutation is implemented by randomly selecting two gene positions and reversing the order of the genes between two positions.

Figure 5.

Repair example: (a) an initial assignment where the circled points are outside the working area of their assigned robots and (b) when the assigned point cannot be reached, it will be redistributed.

Repair strategy

As the chromosome genes are generated in a random arrangement, the task points will be randomly assigned to each robot arm. Consequently, there is a possibility that a certain robot arm may be assigned to an unreachable task point. There are two potential solutions to this problem. The first one is that over more evolutionary generations, feasible chromosome may be generated through natural selection. The second is that, through the repairing method, unreachable tasks points may be randomly reassigned to arms that can reach them. We have chosen the latter option because it avoids the slow convergence of EA. The proposed repairing method employs genetic recombination by removing the genes that violate the kinematics of the robot arm and randomly inserting them into the genetic code segments belonging to the other robot arm that can reach the point.

Fitness evaluation

Based on the above encoding model, the sequence of tasks assigned to each robot arm can be obtained after decoding. However, due to the nature of multiple IKs, each task point may have more than one IK available and the final pose must be decided before calculating the fitness. This paper therefore develops a Greedy_Assignment function in algorithm 2 to provide high efficiency in finding the most suitable IK. The greedy strategy is to select the most favourable branching scheme in the current state, with the expectation that the result will be the best solution. ⁴² Although greedy search is also a local search method, combining it with the NSGA-III loops for global search will result in a fast and efficient hybrid search solution.

Algorithm 2.

Greedy_Assignment( $R_i$ ).

Input:   $R_i$ : the sub-task set assigned to the i-th robot. Begin  For each $r_{i,j}$ in $R_i$ , $r_{i,j}·\mathrm{SelectedIK}=$ 0. /* Initially, not selected */  For j = 1 to $\left|r_i\right|$   If j! =  $\left|r_i\right|$ Then $r_{i,\mathrm{next}}=r_{i,j+1}$ Else $r_{i,\mathrm{next}}=r_{i,1}$ .     Compute and select the minimal operation time for all IK combinations between $r_{i,j}$ and $r_{i,\mathrm{next}}$ ; and set $r_{i,j}·\mathrm{SelectedIK}$  = (the selected IK index).   Next $j$  End

The algorithm selects the expected IK by evaluating the operation time from a point $r_{i,j}$ and its next point $r_{i,j+1}$ . The combination of the IK solutions of these two points that reproduced the smallest joint rotation time is then adopted. This low-cost algorithm is a computationally efficient solution with a computation complexity of O(2m ^l ), where m is the number of task points and l is the number of IKs of each point. This is significantly less than the exhaustive search method, which has a computation complexity of O(m ^l ).

Please refer to Algorithm 4 for a detailed overview of the fitness function. The algorithm of Objective_Constraint_Evaluation develops the proposed greedy search, trajectory collision detection (as shown in Algorithm 3) and objective evaluation.

Algorithm 3.

$\mathrm{\phi }$ (R, $\mathrm{\delta }$ ).

Input:  R: the task assignment matrix.   $\mathrm{\delta }$ : the timing resolution to detect the collision. Output:   $\mathrm{\gamma }$ : Collision count. Begin   $\mathrm{\gamma }$  = 0  For i = 1 to n   For j = 1 to n, $j\ne i$    Compute $T{R}_i$ and $T{R}_j$ according to equation (6).     ${\mathrm{t}}_{\min }$  = min( $T{R}_i$ , $T{R}_j$ ).    For each k in $\{0,\cdots {\mathrm{t}}_{\min }/\mathrm{\delta }\}$      $\mathrm{\gamma }$ + = Collision_Detection(      the pose of the i-th robot at the k-th time tick,      the pose of the j-th robot at the k-th time tick).    Next k   Next j  Next i Return $\mathrm{\gamma }$ End

Algorithm 4.

Objective_Constraint_Evaluation(P, $\mathrm{\delta }$ ).

Input:  P: the population.   $\mathrm{\delta }$ : the timing resolution to evaluate the collision amount. Begin  For each p in P    Decode p to obtain the evolutionary assignment plan $R$ . /* As defined in equation (6)*/   For each row vector ${\mathrm{R}}_i=\{r_{i,1}\cdots r_{i,m}\}$ in $R$    Call Greedy_Assignment( ${\mathrm{R}}_i$ ) to assign best-fit IKs.   Next ${\mathrm{R}}_i$   Call $\mathrm{\phi }$ (R, $\mathrm{\delta }$ ) in Algorithm 3 to obtain the collision count and set it as the constraint value.   Compute objectives ${\mathrm{\phi }}_1,{\mathrm{\phi }}_2$ of p according to equation (6).  Next p End

Constraint handling

Since this optimization problem is highly constrained, we employ appropriate constraint handling rules ⁴³ in conjunction with the tournament selection method. The approach offers the advantage of not requiring an absolute quantification of constraint violation, but rather focusing on the relative merits of the two individuals.

Hamming-distance-based method

This paper proposes a Hamming-distance-based method to enhance the genotype diversity. Hamming distance is defined as the number of different gene codes between two chromosomes. By deleting individuals (chromosomes) in the population whose Hamming distance to the others is too low, the diversity can be further enhanced. The pseudo-code is shown in Algorithm 5.

Algorithm 5.

Hamming_Distance_Threshold( $P_1,P_2,\epsilon )$ .

Input:   $P_1$ : Parent population.   $P_2$ : Offspring population.   $\epsilon$ : Hamming-distance threshold. Output:  Modified $P_1$ . Begin   $P_{1+2}=P_1\cup P_2$ .   $\mathrm{Q}=\varnothing$ .  idxCenter = 1.  nKeptCount = 1.   While (nKeptCount < (the size of $P_1$ )) and (idxCenter < (the size of $P_{1+2}$ ))   next = idxCenter+ 1.  If $P_{1+2}\left(\mathrm{idxCenter}\right)P_{1+2}(\mathrm{next})<\mathrm{\epsilon }$ Then   Move $P_{1+2}\left(\mathrm{next}\right)$ to Q.  Else   idxCenter = next.   nKeptCount = nKeptCount + 1.  End  Next i  Copy $P_{1+2}(1:\mathrm{nKeptCount})$ to $P_1$ .  If $\mathrm{nKeptCount}<(\mathrm{the}\mathrm{size}\mathrm{of}{P}_1)$ Then    Randomly select ( $(\mathrm{the}\mathrm{size}\mathrm{of}{P}_1)-\mathrm{nKeptCount})$ elements of Q to $P_1$ .  End End

Finally, based on the individual genetic operations proposed above, the overall HD-NSGA-III pseudo-steps are shown in Algorithm 6.

Algorithm 6.

HD-NSGA-III.

Input:  (1) $P$ : the parent population.  (2) Genotype (Hamming) distance threshold $\epsilon$ .  (3) The stop criteria (e.g. evolutionary generation limits or sufficiently small objective values). Output:  Pareto solutions. Begin   Preprocess the desired poses and compute their inverse kinematics (IK) solutions.  Initialize the parent population $P$ .  Call Objective_Constraint_Evaluation( $P,\epsilon$ ).  Do    Perform the selection, crossover and mutation on $P$ to obtain the offspring $Q$ .   Perform repair operation on $Q$ .   Call Objective_Constraint_Evaluation( $Q,\epsilon$ ).    Perform non-domination sort on $P\cup Q$ to obtain non-dominated fronts F.    Suppose F =  $\{F_1,\cdots ,F_l\cdots ,F_{\mathrm{end}}\}$ under the condition that (the size of $F_1\cup \cdots F_{l-1}$ ) < (the size of $P$ ) and (the size of $F_1\cup \cdots F_l$ ) ≥ (the size of $P$ ).    If (the size of $F_1\cup \cdots F_l$ ) = (the size of $P$ ) Then $P=F_1\cup \cdots F_l$ ; Q =  $F_{l+1}\cup \cdots F_{\mathrm{end}}$ .   Else     Perform the objective reference points evaluation of NSGA-III on $F_l$ to obtain $F_{l\_\mathrm{reserved}}$ and $F_{l\_\mathrm{removed}}$ .     $P=F_1\cup \cdots F_{l\_\mathrm{reserved}}$ .     $Q=F_{l\_\mathrm{removed}}\cup \cdots F_{\mathrm{end}}$ .   End    If (the stop criteria are not met) Then Call Hamming_Distance_Threshold ( $P$ , $Q$ , $\mathrm{\epsilon })$ .   Else Return $F_1$ . /*Exiting the function. */  Loop End

Experiments and discussions

This paper employs two industrial robot models, FANUC ER-4iA and PUMA560, to investigate the CRTSA issues. The experiments are conducted in a hardware and software environment that includes a CPU i7-12700, Windows 10 and 32G RAM. The NSGA-III implementation was inherited from geatpy, ⁴⁴ and used the pyGeo ⁴⁵ vector geometry method to compute the intersection values between the geometries.

In the manufacturing task setup, the endpoint of the manipulator is the centre of the flange plane in the sixth joints of the robots, and its direction is perpendicular to the manufacturing task point. It is also worth noting that when machining tools are used, they can also be attached by bounding boxes for collision avoidance calculations.

As illustrated in Figure 6, in our proposed experiments, we initially examine the performance of two arms in executing 25, 50 and 75 task points. We randomly generate the required task points in the world coordinates in the range of $x\in [100,800]$ , $y\in [-350,350]$ and $z\in [0,50]$ .

Figure 6.

Example task points: (a) 25 points, (b) 50 points and (c) 75 points.

In this case, the robot arm base is placed as shown in Figure 7. This is a common way of placing two robots opposite to each other; the first robot is positioned at (0,0,0) with base orientation $(0°,0°,0°)$ , while the second robot is positioned at (900,0,0) with base orientation $(180°,0°,0°)$ .

Figure 7.

Placement of dual robots.

Experimental parameter settings

The experiment was conducted with a population size of 50, a crossover rate of 1.0, a mutation rate of 0.1, an evolutionary generation limit of 10,000 and a 10% Hamming-distance threshold. Each case runs 10 times to calculate the PMI.

Collision-free trajectory planning

There are numerous potential approaches for generating a trajectory that contains temporal information for a chosen path. In this paper, we choose the joint space planning method and use the fifth-degree polynomial smoothing method to perform the trajectory interpolation between the path points. ³⁷ This trajectory requires six boundary conditions, including the initial position, end-position, initial velocity, end-point velocity, initial acceleration and endpoint acceleration. Although these conditions can be specified as customized values according to practical needs, here we have set the initial and terminal velocities to 0, which is commonly used in many manufacturing tasks such as pick-and-place or spot welding processes. The two task points are set to be the initial and terminal positions, and the operation time is evaluated by assuming that the maximum average angular velocities of the joints are the same, and taking $\mathrm{\omega }=100\mathrm{ing}$ as the planning basis to calculate the operation time in equation (6), and then convert it back to the acceleration values. In practice, these values can also be specified by a look-up table according to the manufacturing requirements and the specifications of robot arms.

In conjunction with the above planning method, we generate all trajectory points at 20 Hz, which is the permissible trajectory specification for most industrial robots. ³⁷ However, if the corresponding joint angle difference between some trajectory points is too large, the collision detection may be distorted due to the low resolution of the trajectory. To address this, this paper additionally adds a timing resolution threshold $\mathrm{\delta }$ to the collision detection in Algorithm 3. To ensure that when the trajectory points are generated, the collision-free detecting rate can be higher than the generated resolution. Please note that this sublevel collision-free detection process is primarily used as a logical computation to evaluate the collision problem. In the event that the trajectory confirms that no collision occurs, the actual trajectory of the robot arm will still maintain a 20 Hz frequency to output the trajectory points.

Performance measurement index

For a multiobjective optimization problem, different solution schemes will explore their respective sets of Pareto solutions. In general, the Pareto front of each scheme will contain multiple solutions that may not be able to dominate each other. Consequently, it is necessary to use specially designed performance measurement indices (PMIs) to examine the relative merits of two sets of Pareto solutions derived from two distinct schemes.

In this paper, to assess the solution quality of HD-NSGA-III and the original NSGA-III for the CRTSA problem, we use two distinct PMIs, C-function ¹³ and Distribution Metric (DM), ⁴⁶ for comparison. First, the C-function is used to calculate the domination ratio of the two sets of Pareto solutions. This can be formulated as follows:

$$C(X\prime ,X″):=\frac{\left|\{a″\in X″;\exists a\prime \in X\prime :a\prime \preccurlyeq a″\}\right|}{|X″|},$$

(7)

where $X^{\prime }$ and $X^{″}$ are non-dominated solutions (in this case, Pareto solutions), $a^{\prime }\preccurlyeq a^{″}$ means that $a^{\prime }$ dominates or at least equal to $a^{″}$ ; when C( $X^{\prime },X^{″}$ ) = 1, all points of $X^{\prime }$ dominate or equal to all points of $X^{″}$ , and when C( $X^{\prime },X^{″}$ ) = 0, it means that none of the points of $X^{\prime }$ overruns over $X^{″}$ . Therefore, if C( $X^{\prime },X^{″}$ ) > C( $X^{″},X^{\prime }$ ), then $X^{\prime }$ has better quality solutions compared to $X^{″}$ .

In contrast to the C-function, DM measures the proximity of the objective to the minimum and assesses the distribution of Pareto solutions across the solution space. Its equation is as follows ⁴⁶ :

$$DM=\frac{1}{\left|S\right|}\sum _{h=1}^H\frac{{\sigma }_h/{\mu }_h}{U_h/\left|f_h\left(V_G\right)-f_h\left(V_B\right)\right|},$$

(8)

where $\left|S\right|$ is the number of non-dominated solutions, H is the number of objectives, $U_h$ is the difference between the largest and smallest objective values in the objective h, $f_h\left(s\right)$ is the function value of the h-th objective at current solution s, ${\sigma }_h=\frac{1}{\left|S\right|-2}\sum _{e=1}^{\left|S\right|-1}{({d^{\prime }}_e^h-\overline{{d^{\prime }}^h})}^2$ , ${\mu }_h=\frac{\sum _{e=1}^{\left|S\right|-1}{d^{\prime }}_e^h}{\left|S\right|-2}$ , ${d^{\prime }}^h$ is the distance to the objective h, $\overline{{d^{\prime }}^h}$ is the mean value of ${d^{\prime }}^h$ , and $V_B$ and $V_G$ are the Utopia and Nadir design points, which are composed of the best and worst values of objectives in the Pareto front, respectively. For further details, please refer to Zheng et al. ⁴⁶ In accordance with the definition, a smaller the value of $\mathrm{DM}$ accurate representation of the distributional nature of the Pareto front.

Results and comparison

The results of the comparison between HD-NSGA-III and NSGA-III using the C-function and DM are presented in Tables 3 and 4.

Table 3.

Comparison of optimization results in terms of C-function.

Case	FANUC (A: HD-NSGA-III, B: NSGA-III)			PUMA560 (A: HD-NSGA-III, B: NSGA-III)
	C(A,B) > C(B,A)	C(A,B) < C(B,A)	C(A,B) = C(B,A)	C(A,B) > C(B,A)	C(A,B) < C(B,A)	C(A,B) = C(B,A)
25	5	5	0	7	3	0
50	6	4	0	6	4	0
75	6	3	1	7	3	0

Table 4.

Comparison of optimization results with DM.

Case	FANUC ER-4iA		PUMA560
	$\mathrm{DM}(\mathrm{NSGA}-\mathrm{III})$	$\mathrm{DM}(\mathrm{HD}-\mathrm{NSGA}-\mathrm{III})$	$\mathrm{DM}(\mathrm{NSGA}-\mathrm{III})$	$\mathrm{DM}(\mathrm{HD}-\mathrm{NSGA}-\mathrm{III})$
25	0.628 (0.467)	0.289 (0.158)	0.517 (0.359)	0.302 (0.191)
50	0.668 (0.166)	0.543 (0.298)	0.458 (0.305)	0.386 (0.432)
75	2.256 (1.681)	1.038 (0.568)	0.327 (0.373)	0.190 (0.083)

As illustrated in Table 3, the two sets of Pareto solutions have comparable outcomes in some test instances. Both can converge to similar objective values, approaching to the minimum. However, HD-NSGA-III has demonstrated slightly better in most C-function test instances. Conversely, as illustrated in Table 4, in DM, HD-NSGA-III demonstrably outperforms NSGA-III in terms of quality. This indicates that HD-NSGA-III explores a better distribution of objectives in phenotypic space. Consequently, the subsequent experiments will be conducted using HD-NSGA-III.

Regarding the obtained assignment solutions, we will illustrate the details of the case that was executed with 25 points. Please refer to Table 5 for a detailed breakdown of the task points. Once the algorithm has been executed, the Pareto fronts are obtained as shown in Figure 8. In particular, the Pareto solutions are listed in Table 6, which describes the objectives and the sub-task sequence assignments are described. For example, the first solution assigns the task sequences (20, 7, 1, 2, 12, 3, 19, 25, 4, 9, 15, 21, 23, 18) and (8, 5, 16, 24, 13, 6, 22, 10, 14, 11, 17) to the two robots, with the objective values for the manufacturing time and the balance of task loading being 5.808 and 0.000006, respectively.

Table 5.

Task data used in the case of 25 task points.

No	$(\mathrm{x},\mathrm{y},\mathrm{z},\hat{\mathrm{x}},\hat{\mathrm{y}},\hat{\mathrm{z}})$	No	$(\mathrm{x},\mathrm{y},\mathrm{z},\hat{\mathrm{x}},\hat{\mathrm{y}},\hat{\mathrm{z}})$	No	$(\mathrm{x},\mathrm{y},\mathrm{z},\hat{\mathrm{x}},\hat{\mathrm{y}},\hat{\mathrm{z}})$
1	(357.23, −134.48, 18.85, 0, 0, 0.13)	10	(545.73, −310.78, 41.92, 0, 0, 0.24)	19	(481.79, 64.43, 35.14, 0, 0, −0.25)
2	(440.43, −128.51, 46.97, 0, 0, 0.24)	11	(645.22, 250.82, 41.40, 0, 0, 0.22)	20	(198.30, 15.04, 9.26, 0, 0, −0.44)
3	(401.12, −305.23, 29.88, 0, 0, −0.73)	12	(433.87, −106.32, 0.26, 0, 0, 0.37)	21	(292.70, 208.68, 28.03, 0, 0, −0.48)
4	(420.37, 167.04, 14.87, 0, 0, −0.17)	13	(617.96, −289.90, 7.95, 0, 0, −0.31)	22	(470.58, −282.55, 27.80, 0, 0, 0.40)
5	(646.11, 161.21, 48.04, 0, 0, 0.72)	14	(499.32, 252.20, 21.12, 0, 0, 0.36)	23	(251.96, 270.87, 27.60, 0, 0, 0.26)
6	(534.41, −228.30, 39.15, 0, 0, 0.50)	15	(189.51, 319.36, 20.05, 0, 0, 0.66)	24	(622.76, −280.76, 20.19, 0, 0, 0.72)
7	(170.17, −195.46, 28.33, 0, 0, 0.45)	16	(592.89, 344.96, 30.08, 0, 0, 0.33)	25	(456.41, 142.89, 16.54, 0, 0, 0.10)
8	(592.21, −15.79, 6.70, 0, 0, −0.35)	17	(621.58, 67.64, 22.76, 0, 0, −0.60)
9	(429.11, 236.44, 38.00, 0, 0, 0.40)	18	(304.12, −164.71, 8.99, 0, 0, −0.11)

Figure 8.

Pareto solutions in the case of 25 task points assigned to two robots: (a) FAUNC ER-4iA and (b) PUMA560.

Table 6.

Solutions of the case of 25 task points of HD-NSGA-III.

No	$({\varphi }_1,{\varphi }_2)$	FANUC ER-4iA #1	FANUC ER-4iA #2
1	(5.808, 0.000006)	20, 7, 1, 2, 12, 3, 19, 25, 4, 9, 15, 21, 23, 18	8, 5, 16, 24, 13, 6, 22, 10, 14, 11, 17
2	(4.983, 0.000013)	23, 15, 1, 9, 4, 25, 19, 3, 2, 12, 18, 7, 20, 21	22, 10, 13, 8, 14, 16, 11, 5, 24, 6, 17
3	(4.758, 0.000020)	23, 20, 21, 18, 7, 1, 2, 3, 12, 19, 25, 4, 9, 15	6, 24, 8, 13, 17, 16, 14, 11, 5, 22, 10
4	(4.339, 0.000049)	15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23, 21	10, 22, 6, 13, 24, 8, 11, 14, 16, 5, 17
5		23, 20, 7, 18, 1, 2, 12, 3, 19, 25, 4, 9, 15, 21	10, 22, 6, 24, 13, 8, 17, 11, 16, 14, 5
6		23, 20, 7, 18, 1, 2, 12, 3, 19, 25, 4, 9, 15, 21	10, 22, 6, 24, 13, 5, 14, 16, 11, 17, 8
7		21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23	10, 22, 6, 24, 13, 8, 5, 14, 16, 11, 17
8		21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23	10, 22, 6, 24, 13, 8, 11, 16, 14, 5, 17
9		23, 21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20	10, 22, 6, 24, 13, 17, 11, 16, 14, 5, 8
10		21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23	10, 22, 6, 24, 13, 8, 11, 16, 14, 5, 17
11		21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23	10, 22, 6, 24, 13, 8, 17, 11, 16, 14, 5
12		21, 15, 9, 4, 25, 19, 3, 12, 2, 1, 18, 7, 20, 23	17, 5, 14, 16, 11, 8, 13, 24, 6, 22, 10
13	(4.288, 0.000136)	23, 21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20	22, 10, 13, 17, 16, 14, 11, 5, 8, 24, 6
14		23, 20, 18, 7, 1, 12, 2, 3, 19, 4, 25, 9, 15, 21	22, 10, 13, 17, 16, 14, 11, 5, 8, 24, 6
15		23, 21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
16		23, 21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
17		23, 21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
18		21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20, 23	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
19		23, 20, 18, 7, 1, 12, 2, 3, 19, 4, 25, 9, 15, 21	22, 10, 13, 5, 11, 16, 14, 17, 8, 24, 6
20		21, 15, 9, 25, 4, 19, 3, 2, 12, 1, 7, 18, 20, 23	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
21		23, 20, 18, 7, 1, 12, 2, 3, 19, 4, 25, 9, 15, 21	22, 10, 13, 17, 14, 16, 11, 5, 8, 24, 6
22	(4.270, 0.000260)	23, 20, 18, 7, 1, 2, 3, 12, 19, 25, 4, 9, 15, 21	22, 6, 14, 11, 16, 5, 17, 8, 24, 13, 10
23	(4.097, 0.000483)	23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	17, 14, 16, 11, 5, 8, 13, 24, 6, 10, 22
24		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	5, 11, 14, 16, 17, 8, 24, 13, 6, 10, 22
25		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	8, 5, 11, 16, 14, 17, 13, 24, 6, 10, 22
26		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	17, 16, 14, 11, 5, 8, 13, 24, 6, 10, 22
27		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	17, 16, 14, 11, 5, 8, 24, 13, 6, 10, 22
28		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	8, 5, 11, 16, 14, 17, 24, 13, 6, 10, 22
29		15, 23, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	8, 5, 11, 14, 16, 17, 24, 13, 6, 10, 22
30		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	17, 16, 14, 11, 5, 8, 13, 24, 6, 10, 22
31		23, 15, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	5, 11, 16, 14, 17, 8, 24, 13, 6, 10, 22
32		15, 23, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 21	17, 16, 14, 11, 5, 8, 24, 13, 6, 10, 22
33	(4.072, 0.000690)	23, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4, 15	22, 24, 13, 14, 16, 11, 5, 17, 8, 10, 6
34	(4.034, 0.001372)	23, 15, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	6, 16, 14, 11, 5, 17, 8, 13, 24, 22, 10
35		23, 15, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	8, 17, 5, 11, 14, 16, 6, 13, 24, 22, 10
36		23, 15, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	6, 16, 14, 11, 5, 17, 8, 13, 24, 22, 10
37		23, 15, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	8, 16, 14, 11, 5, 17, 6, 13, 24, 22, 10
38		15, 23, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	6, 16, 14, 11, 5, 17, 8, 13, 24, 22, 10
39		15, 23, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	8, 17, 5, 11, 14, 16, 6, 13, 24, 22, 10
40		15, 23, 21, 20, 7, 18, 1, 12, 3, 2, 19, 25, 9, 4	8, 16, 14, 11, 5, 17, 6, 13, 24, 22, 10
41	(4.031, 0.001860)	23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 17, 14, 16, 11, 5, 8, 24
42		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 17, 14, 16, 11, 5, 8, 24
43		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 8, 5, 11, 16, 14, 17, 6, 24
44		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 5, 11, 16, 14, 17, 8, 24
45		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 5, 11, 16, 14, 17, 8, 24
46		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 8, 5, 11, 16, 14, 17, 6, 24
47		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 5, 11, 16, 14, 17, 8, 24
48		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 8, 5, 11, 16, 14, 17, 6, 24
49		23, 15, 4, 9, 25, 19, 2, 3, 12, 1, 18, 7, 20, 21	22, 10, 13, 6, 17, 14, 16, 11, 5, 8, 24
50	(3.945, 0.007814)	9, 25, 4, 19, 2, 3, 12, 1, 18, 7, 20, 21, 15, 23	8, 17, 5, 11, 16, 14, 22, 13, 10, 6, 24

Furthermore, the table displays 11 objective combinations located at the Pareto front. However, all the different phenotype combinations assigned to either the first or second robot are present in the 50 populations. Therefore, the diversity in the phenotype space is maintained to the greatest extent possible.

Expanded application for four cooperative robotic arms

The proposed method can be extended to set up additional robot arms and configurations for various applications. In this experiment, we take four FANUC ER-4iA robot arms with linear and diagonal configurations, as illustrated in Figure 9 and in Table 7. Due to the larger number of arms, we use 75 task points with more evolutionary generations (100,000 generations) for exploring the global solution. After the execution, the proposed HD-NSGA-III also successfully identifies several non-dominated solutions under the manufacturing conditions. The resulting Pareto fronts are presented in Figure 10. Figure 11 shows examples of Pareto solutions for task sequence assignments. As shown in the figure, the markers are the task points and the dashed lines are the manufacturing paths. The algorithm schedules the task sequences for the robot arms under collision-free conditions. Note that for the sake of the screen clarity, the robot pose changes between task points are not shown in the figures.

Figure 9.

Different configurations with four FANUC robots: (a) diagonal configure and (b) linear configure.

Table 7.

Robot locations in the configurations.

Robot no.	Diagonal configuration		Linear configuration
	Site location	Site direction	Site location	Site direction
1	$(0,0,0)$	$(0°,0°,0°)$	(0, 400, 0)	$(0°,0°,0°)$
2	$(900,0,0)$	$(180°,0°,0°)$	(0, −400, 0)	$(0°,0°,0°)$
3	$(450,450,0)$	$(-90°,0°,0°)$	(0, 1200, 0)	$(0°,0°,0°)$
4	$(450,-450,0)$	$(90,0°,0°)$	(0, −1200, 0)	$(0°,0°,0°)$

Figure 10.

Pareto solutions for 75 task points assigned to four FANUC robots: (a) Pareto front in the diagonal case and (b) Pareto front in the linear case.

Figure 11.

Example task sequence assignments in the two configurations: (a) example solution in the diagonal case and (b) example solution in the linear case.

The main reason why more evolutionary generations are needed than in the case of two robot arms is that the volume of the robot arms occupies a larger workspace. If more robot arms are arranged in a limited area, the operating space of the robot arms becomes too crowded, which increases the probability of collision and causes more constraint violations. Consequently, more generations are needed to traverse a greater number of combinations of task sequences during the actual execution.

Conclusions

This paper presents a transformation of the collaborative robotic task sequence assignment (CRTSA) problem into a multiobjective optimization problem with two objectives: task completion time and task load balancing. The optimal CRTSA problem contains a greater number of combinations within the solution space than the multi-traveling salesman (MTSP) problem. This is due to the existence of multiple IKs that can be taken to manipulate robots to reach task points. Moreover, the necessity for multiple robot arms to work together increases the likelihood of collisions between robotic motions, which is a significant and critical issue. The process of assigning automated task sequences to robots becomes increasingly complex. Therefore, this paper presents the development of a hybrid NSGA-III (HD-NSGA-III) approach, which integrates three distinct methods: a Hamming distance-based method, a geometric inference method and a greedy method. The hybridization approach ensures the quality of the solution.

To assess the efficacy and scalability of the proposed scheme, we have conducted simulated experiments on different scale sizes using two and four industrial robots. The experiments demonstrate that the proposed approach is effective in obtaining Pareto solutions that contain collision-free robotic trajectories. Furthermore, the validation results of different performance measurement indices (PMIs), including the C-function and the DM metric, demonstrate that the proposed HD-NSGA-III can obtain a superior solution quality for the optimal CRTSA problem.

Nevertheless, the findings of this study are subject to certain limitations if the practical trajectory planning has different requirements. One such example is a hybrid trajectory, which necessitates the inclusion of a constant-velocity linear segment with parabolic blends at the beginning and end of the segment. In this case, it is necessary to make adjustments to the trajectory generation implementation.

Further research could investigate the issue of stable robot control under heavy loads. This would require a detailed study of robot dynamics models during trajectory design. Additionally, it would be beneficial to investigate more flexible planning scenarios. One possible research direction is to explore the integration of additional robot types and robot site locations in task sequencing planning, provided that robot reconfiguration is a viable option.