2.1. Collaborative system model

As shown in Figure 1, the base coordinate system $\{\mathrm{O}_{\mathrm{B}}\}$ of the parallel robot is established at the center of its static platform, and its tool reference coordinate system $\{\mathrm{O}_{\mathrm{P}}\}$ is established at the center of its moving platform. The base coordinate system $\{\mathrm{O}'_{\mathrm{B}}\}$ of the serial robot is established at its base, and its tool coordinate system $\{\mathrm{O}_{\mathrm{T}}\}$ is established at its end-effector. The coordinate system $\{\mathrm{O}_{\mathrm{W}}\}$ represents the world coordinate system of the entire collaborative assembly system. For analytical convenience, the world coordinate system coincides with the base coordinate system of the parallel robot [Reference Romdhane24].

Figure 1.

Schematic diagram of a collaborative assembly system with serial-parallel robots.

2.2. Kinematic description

Let the transformation matrices from the base coordinate systems of the serial and parallel robots to the world coordinate system be ${_{B}^{W}{\boldsymbol{T}}}{}$ and $_{\mathrm{B}'}^{W}{\boldsymbol{T}}{}$ , respectively. Meanwhile, the transformation matrices from their corresponding tool coordinate systems to the base coordinate system are $_{P}^{B}{\boldsymbol{T}}{}$ and $_{T}^{\mathrm{B}'}{\boldsymbol{T}}{}$ , respectively. Besides, let $_{P}^{T}{\boldsymbol{T}}{}$ denote the transformation matrix between the end-tool coordinate system of the parallel robot and that of the serial robot. Then one can obtain:

(1) \begin{equation} ^{W}_{B}{\boldsymbol{T}}{_{P}^{B}{\boldsymbol{T}}}={_{\mathrm{B}'}^{W}{\boldsymbol{T}}}{_{T}^{\mathrm{B}'}{\boldsymbol{T}}} {_{P}^{T}{\boldsymbol{T}}} \end{equation}

Multiplying both sides of the Eq. (1) by $(_{B}^{W}{\boldsymbol{T}}{})^{-1}$ on the left yields:

(2) \begin{equation} {_{P}^{B}{\boldsymbol{T}}}={(_{B}^{W}{\boldsymbol{T}}})^{-1}{_{\mathrm{B}'}^{W}{\boldsymbol{T}}}{_{T}^{\mathrm{B}'}{\boldsymbol{T}}}{_{P}^{T}{\boldsymbol{T}}} \end{equation}

In a dual-robot collaborative assembly system with determined robot structure and installation layout, matrices $(_{B}^{W}{\boldsymbol{T}}{})^{-1}$ and $_{\mathrm{B}'}^{W}{\boldsymbol{T}}{}$ are constant. Matrix $_{P}^{T}{\boldsymbol{T}}{}$ is determined by the work task of system. Once the work task is defined, the corresponding joint trajectories of the robots can be solved, enabling coordinated motion within the collaborative system.

2.3. Workspace determination

The collaborative workspace of a dual-robot system is defined as the intersection of the workspaces of the individual robots involved [Reference Xu, Zhao, He, He and Wu25]. Therefore, to accurately define the collaborative workspace, the workspaces of the individual robots must be precisely characterized first. For the collaborative assembly system shown in Figure 1, the workspace of the serial robot is solved by using the Monte Carlo method [Reference Peidró, Reinoso, Gil, Marín and Payá26], and that of the parallel robot is solved by using the lattice method.

The specific steps are as follows: first, solve for the point cloud set $W_{\text{serial}}$ of the serial robot’s workspace under coordinate system $\{\mathrm{O}'_{\mathrm{B}}\}$ , then transform it to obtain the point cloud set $W'_{\text{serial}}$ under coordinate system $\{\mathrm{O}_{\mathrm{B}}\}$ . Substitute the points in $W'_{\text{serial}}$ into the parallel robot workspace solution algorithm, retain the point cloud data within the parallel robot workspace, and obtain the point cloud set $W_{\text{paral}}$ under the coordinate system $\{\mathrm{O}_{\mathrm{B}}\}$ . After coordinate transformation, the point cloud set $W$ under the world coordinate system $\{\mathrm{O}_{\mathrm{W}}\}$ is obtained, which represents the collaborative workspace of the system [Reference Rega, Marino, Pasquariello, Vitolo, Patalano, Zanella and Lanzotti27]. The algorithm flow is shown in Figure 2.

Figure 2.

Collaborative workspace solution flowchart.

This method leverages the simplicity of the forward kinematic solution for serial robots and the inverse kinematic solution for parallel robots. It thus avoids information loss when solving overlapping parts and ensures uniform point cloud density within the collaborative workspace.

3.1. Selection of decision variables

For the collaborative assembly system designed in this paper, the key parameters to be optimized are the radii of the moving and static platforms of the parallel robots, the distribution of the articulation points, the limiting lengths of the legs of the parallel robots, and the layout parameters between the series and the parallel dual robots. After the selection and necessary processing of these main parameters, the decision variables can be derived as:

(3) \begin{equation} \begin{array}{ll} {\boldsymbol{X}} & =\left(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6},x_{7},x_{8},x_{9}\right)^{\mathrm{T}}\\ & =\left(\tau _{R},\alpha _{\mathrm{P}},R_{\mathrm{B}},\alpha _{\mathrm{B}},\tau _{l},l_{\max },X_{\mathrm{rel}},Z_{\mathrm{rel}},\beta _{\mathrm{rel}}\right)^{\mathrm{T}} \end{array} \end{equation}

where the length is in mm, the angle is in °, and the physical significance of each component is:

$\tau _{R}$ : The ratio of the radius of the outer circle of the moving platform $R_{\mathrm{P}}$ to the radius of the outer circle of the stationary platform $R_{\mathrm{B}}$ ;

$\alpha _{\mathrm{P}}$ : Angle of the center of gravity of the same set of ball hinges of the moving platform;

$\alpha _{\mathrm{B}}$ : Static platform with the same set of hook hinge tensor circumcenter angle;

$\tau _{l}$ : Ratio of minimum leg length $l_{\min }$ to maximum leg length $l_{\max }$ ;

$X_{\mathrm{rel}},Z_{\mathrm{rel}},\beta _{\mathrm{rel}}$ : Define the pose of the tandem robot’s base frame relative to the world coordinate system (i.e., the parallel robot base frame). Specifically, $X_{\mathrm{rel}}$ and $Z_{\mathrm{rel}}$ are the translational offsets along the X and Z axes, respectively, and $\beta _{\mathrm{rel}}$ is the rotational offset (yaw angle) about the Y axis.

3.2. Determination of constraints

According to the design requirements, the constraints are given and classified into linear and non-linear constraints according to their nature.

1. (1) Linear constraints: Considering factors such as the manufacture and installation of parallel robots, the installation of serial robots, and the safety of collaboration between dual robots, the decision variable should satisfy the conditions described as follows:

   (4) \begin{equation} x_{i}\in \left[x_{i\min },\,x_{i\max }\right]\begin{array}{c} \begin{array}{c} \end{array} \end{array}\begin{array}{c} \left(i=1,2,\,\cdots, 9\right) \end{array} \end{equation}
   where $x_{i\min }$ and $x_{i\max }$ are the lower and upper bounds of the decision variables, respectively. In this research, they are assigned to:
   \begin{align*} {\boldsymbol{X}}_{\min } & =\left(0.1,10,200,10,0.5,600,-800,0,0\right)^{\mathrm{T}}\\[4pt]{\boldsymbol{X}}_{\max } & =\left(0.9,110,600,110,0.9,1000,-200,200,30\right)^{\mathrm{T}} \end{align*}

2. (2) Nonlinear constraints: Since the configuration of the serial robots in the collaborative assembly system has been determined, their workspace $W_{\text{serial}}$ has been determined, so the collaborative workspace $W({\boldsymbol{X}})$ of the collaborative system must satisfy:

   (5) \begin{equation} W\left({\boldsymbol{X}}\right)\subset W_{\text{serial}} \end{equation}

In addition, for the collaborative workspace $W({\boldsymbol{X}})$ , there are constraints on the parallel robot Hook hinge pendulum angle $\theta_{B}({\boldsymbol{X}})$ and ball hinge pendulum angle $\theta _{\mathrm{P}}({\boldsymbol{X}})$ in the collaborative system as follows:

(6) \begin{equation} \left\{\begin{array}{l} \theta_{B}({\boldsymbol{X}}) \leq \theta_{Bmax}\\[5pt]\theta_{P}({\boldsymbol{X}}) \leq \theta_{Pmax} \end{array}\right. \end{equation}

where $\theta_{Bmax}^{}{}{}$ and $\theta_{Pmax}^{}{}{}$ are the upper boundaries of the swing angle ranges of the hook and ball hinges, respectively, with $\theta_{Bmax}^{}{}{}$ smaller than 90° and $\theta_{Pmax}^{}{}{}$ smaller than 45° in general.

The swing angle $\theta \leq \theta_{\max }^{}{}{}$ constraint is enforced via a penalty function in the objective: $\kappa \cdot \max (0,$ $\theta -\theta _{\max })^{2}$ and all constraints are handled through this unified penalty framework. Besides $\kappa$ is an adaptive penalty coefficient (initialized to 10³ and scaled by 10 per infeasible iteration). This ensures gradient-driven constraint satisfaction, while also not interfering with the evolutionary mechanism of the subsequent algorithm.

3.3. Establishment of the objective function

In practical applications of collaborative assembly systems, different application scenarios have different requirements for the performance indicators of the system. This paper focuses on the flexibility of the assembly process, particularly the size of the collaborative system’s workspace and its mechanical properties, leading to the establishment of the following objective functions:

(1) Collaborative workspace volume function: Using the collaborative workspace volume calculation method shown in Figure 2, the workspace volume $V$ of the collaborative system corresponding to the decision variable is obtained. Define $f_{1}({\boldsymbol{X}})$ as the collaborative workspace volume function:

(7) \begin{equation} f_{1}\left({\boldsymbol{X}}\right)=-V \end{equation}

Therefore, the smaller the $f_{1}({\boldsymbol{X}})$ is, the larger the collaborative workspace volume of the collaborative system.

(2) Force transfer isotropic performance index function: When the configuration of the 6-SPS robot is determined, the length of each leg remains constant, forming a support chain with the hook and ball hinges at both ends, connecting the moving and static platforms into a stable structure. If a six-dimensional force vector acts on the moving platform, there will exist a corresponding reaction force on each leg, as depicted in Figure 3.

Figure 3.

Force analysis of parallel robot.

In the moving platform coordinate system $\{\mathrm{O}_{\mathrm{P}}\}$ , considering the static equilibrium of the moving platform and only the reaction force along the direction of the legs, the mathematical model based on screw theory is established:

(8) \begin{equation} \sum _{i=1}^{6}f_{i}{\boldsymbol{\$ }}_{i}={\boldsymbol{F}}+\in {\boldsymbol{M}} \end{equation}

where ${\boldsymbol{f}}=(f_{1},f_{2},f_{3},f_{4},f_{5},f_{6})^{\mathrm{T}}$ is the reaction force on the 6 legs, ${\boldsymbol{\$}}_{i}$ is the unit vector of the i-th leg axis to the fixed coordinate system, ${\boldsymbol{F}}=(F_{x},F_{y},F_{z})^{\mathrm{T}}$ is the principal vector of the external force, and ${\boldsymbol{M}}=(M_{x},M_{y},M_{z})^{\mathrm{T}}$ is the principal moment to the fixed coordinate origin.

Equation (8) can be rewritten in matrix form, dividing the matrix of first-order static influence coefficients ${\boldsymbol{G}}$ into a force operability matrix ${\boldsymbol{G}}_{f}^{F}$ and a torque operability matrix ${\boldsymbol{G}}_{f}^{M}$ :

(9) \begin{equation} \left\{\begin{array}{l} {\boldsymbol{F}}\\ {\boldsymbol{M}} \end{array}\right\}={\boldsymbol{Gf}}=\left[\begin{array}{l} {\boldsymbol{G}}_{f}^{F}\\[6pt] {\boldsymbol{G}}_{f}^{M} \end{array}\right]{\boldsymbol{f}} \end{equation}

where ${\boldsymbol{G}}_{f}^{F}$ determines the transfer relationship between the parallel robot leg driving force and the end output force, indicating the force transfer performance of the parallel robot. ${\boldsymbol{G}}_{f}^{M}$ determines the transfer relationship between the driving force of legs in parallel robot and the end output torque, indicating the torque transfer performance of the parallel robot.

According to the principle of virtual work, when maintaining a static equilibrium condition, the sum of the virtual work done by all the active forces acting on the particle system in any virtual displacement is equal to zero. Therefore, for the 6-SPS parallel robot, there is:

(10) \begin{equation} \left\{\begin{array}{ll} {\boldsymbol{F}} & {\boldsymbol{M}} \end{array}\right\}\delta {\boldsymbol{L}}={\boldsymbol{f}}^{\mathrm{T}}\delta {\boldsymbol{l}} \end{equation}

(11) \begin{equation} \delta {\boldsymbol{L}}\boldsymbol{=}{\boldsymbol{J}}\delta {\boldsymbol{l}} \end{equation}

where $\delta {\boldsymbol{L}}$ represents the virtual displacement vector at the end of the robot, $\delta {\boldsymbol{l}}$ represents the virtual displacement vector in the robot’s joint space, and ${\boldsymbol{J}}$ is the Jacobian matrix of the robot. Therefore, the relationship between the robot Jacobian matrix ${\boldsymbol{J}}$ and the robot first-order static influence coefficient matrix ${\boldsymbol{G}}^{-1}$ is as follows.

(12) \begin{equation} {\boldsymbol{J}}=\left({\boldsymbol{G}}^{-1}\right)^{\mathrm{T}} \end{equation}

The force transfer isotropic performance index $\varphi _{F}$ and the torque transfer isotropic performance index $\varphi _{M}$ can be defined, respectively:

(13) \begin{equation} \varphi _{F}=\frac{\lambda _{F\min }}{\lambda _{F\max }} \end{equation}

(14) \begin{equation} \varphi _{M}=\frac{\lambda _{M\min }}{\lambda _{M\max }} \end{equation}

where $\lambda _{F\min }$ and $\lambda _{F\max }$ are the minimum and maximum singular values obtained from the singular value decomposition of matrix ${\boldsymbol{G}}_{f}^{F}$ , respectively. The $\lambda _{M\min }$ and $\lambda _{M\max }$ are the minimum and maximum singular values obtained from the singular value decomposition of matrix ${\boldsymbol{G}}_{f}^{M}$ , respectively. Larger values of $\varphi _{F}$ and $\varphi _{M}$ indicate better isotropic performance for force and torque transmission, respectively.

The force isotropy index $\varphi =\lambda _{\min }/\lambda _{\max }$ quantifies the uniformity of force transmission capability in the workspace. In screw theory terms, the Jacobian ${\boldsymbol{J}}$ maps the actuator wrench $\boldsymbol{\tau }$ to the end-effector force/torque ${\boldsymbol{F}}$ via ${\boldsymbol{F}}={\boldsymbol{J}}^{\mathrm{T}}\boldsymbol{\tau }$ . The condition number $\text{cond}({\boldsymbol{J}})=\lambda _{\max }/\lambda _{\min }$ directly measures the dexterity of force transmission. Minimizing $\text{cond}({\boldsymbol{J}})$ , which is equivalent to maximizing $\varphi$ , ensures minimal directional bias, which is optimal for collaborative tasks requiring uniform load capacity [Reference Jia, Huang, Wang and Li23].

To evaluate the force transfer isotropic performance of the collaborative system under different decision variables ${\boldsymbol{X}}$ , the average force transfer isotropic performance index $\Gamma_{F}$ and average torque isotropic transfer performance index $\Gamma_{M}$ are defined within the corresponding collaborative workspace $W({\boldsymbol{X}})$ :

(15) \begin{equation} \Gamma_{F}=\frac{1}{N_{W}}\sum _{j=1}^{N_{W}}\varphi _{Fj} \end{equation}

(16) \begin{equation} \Gamma_{M}=\frac{1}{N_{W}}\sum _{j=1}^{N_{W}}\varphi _{Mj} \end{equation}

where $N_{W}$ is the number of points contained in the $W({\boldsymbol{X}})$ point set.

The force transfer isotropic performance index function $f_{2}({\boldsymbol{X}})$ is defined as:

(17) \begin{equation} f_{2}\left({\boldsymbol{X}}\right)=-\Gamma_{F} \end{equation}

with $-1\leq f_{2}({\boldsymbol{X}})\leq 0$ . A smaller $f_{2}({\boldsymbol{X}})$ indicates better force transfer isotropic performance of the system.

(3) Torque transfer isotropic performance index function: Based on Eq. (16), $f_{3}({\boldsymbol{X}})$ can be defined as the torque transfer isotropic performance index function.

(18) \begin{equation} f_{3}\left({\boldsymbol{X}}\right)=-\Gamma_{M} \end{equation}

where the inequality $-1\leq f_{3}({\boldsymbol{X}})\leq 0$ is satisfied, and when $f_{3}({\boldsymbol{X}})$ is smaller, the isotropic performance of the collaborative system for torque transmission is better.

(4) Comprehensive force transfer performance index function: The comprehensive force transfer performance index $k_{F}$ and comprehensive torque transfer performance index $k_{M}$ are defined as:

(19) \begin{equation} k_{F}=\prod _{i=1}^{3}\lambda _{Fi}\begin{array}{ll} & \left(i=1,2,3\right) \end{array} \end{equation}

(20) \begin{equation} k_{M}=\prod _{i=1}^{3}\lambda _{Mi}\begin{array}{ll} & \left(i=1,2,3\right) \end{array} \end{equation}

where $\lambda _{Fi}$ is the i-th singular value from the singular value decomposition of matrix ${\boldsymbol{G}}_{f}^{F}$ , and $\lambda _{Mi}$ is the i-th singular value from the singular value decomposition of matrix ${\boldsymbol{G}}_{f}^{M}$ . These indices reflect the overall force and torque transfer performance of the robot.

To evaluate the force transfer performance of the collaborative system under different decision variables ${\boldsymbol{X}}$ , the average comprehensive force transfer performance index $K_{F}$ and average comprehensive torque transfer performance index $K_{M}$ are defined within the corresponding collaborative workspace $W({\boldsymbol{X}})$ :

(21) \begin{equation} K_{F}=\frac{1}{N_{W}}\sum _{j=1}^{N_{W}}k_{Fj} \end{equation}

(22) \begin{equation} K_{M}=\frac{1}{N_{W}}\sum _{j=1}^{N_{W}}k_{Mj} \end{equation}

where $N_{W}$ is the number of points in the $W({\boldsymbol{X}})$ point set.

The comprehensive force transfer performance index function $f_{4}({\boldsymbol{X}})$ is defined as:

(23) \begin{equation} f_{4}\left({\boldsymbol{X}}\right)=-K_{F} \end{equation}

with $f_{4}({\boldsymbol{X}})\leq 0$ . A smaller $f_{4}({\boldsymbol{X}})$ indicates better comprehensive force transfer performance of the robot.

(5) Comprehensive torque transfer performance index function: Based on Eq. (22), the comprehensive torque transfer performance index function $f_{5}({\boldsymbol{X}})$ is defined as:

(24) \begin{equation} f_{5}\left({\boldsymbol{X}}\right)=-K_{M} \end{equation}

with $f_{5}({\boldsymbol{X}})\leq 0$ . A smaller $f_{5}({\boldsymbol{X}})$ indicates better comprehensive torque transfer performance of the robot.

The five objective functions have a complex mutual constraint relationship; consequently, they can hardly reach their own optimal solution simultaneously. From the application point of view, the goal is to find a set of equilibrium solutions balanced across each objective function, where no other solution can improve all objectives or one objective without worsening others.

Therefore, the multi-objective optimization model for the collaborative assembly system is formulated as:

(25) \begin{equation} \left\{\begin{array}{l@{\quad}l} \text{find} & {\boldsymbol{X}}\\[5pt] \min & {\boldsymbol{Y}}\left({\boldsymbol{X}}\right)=\left(f_{1}\left({\boldsymbol{X}}\right)\!,\ f_{2}\left({\boldsymbol{X}}\right)\!, \ f_{3}\left({\boldsymbol{X}}\right)\!,\ f_{4}\left({\boldsymbol{X}}\right)\!,\ f_{5}\left({\boldsymbol{X}}\right)\right)^{\mathrm{T}}\\[5pt] \mathrm{s}.\mathrm{t}. & x_{\min }\leq x_{i}\leq x_{\max}\,\, i=1,2,\cdots,9\\[5pt] & W\left({\boldsymbol{X}}\right)\in W_{\text{serial}}\\[5pt] & \begin{array}{l} \theta_{B}^{}{\left({\boldsymbol{X}}\right)}{}\leq \theta_{Bmax}^{}{}{}\\[5pt] \theta_{P}^{}{\left({\boldsymbol{X}}\right)}{}\leq \theta_{Pmax} \end{array} \end{array}\right. \end{equation}

Given the significant disparity in the magnitudes of the five objective functions (e.g., workspace volume on the order of 10⁷ and force parameters within [0, 10]), the raw values are normalized before performing non-dominated sorting and crowding distance calculation. This prevents the optimization from being biased towards the objectives with larger numerical ranges. For each generation and for each objective function $f_{i}$ , the values are normalized to the interval [0, 1] by using a min–max approach based on the current population’s minimum ( $f_{i,\min }$ ) and maximum ( $f_{i,\max }$ ) values:

(26) \begin{equation} f_{i,\text{norm}}=\frac{f_{i}-f_{i,\min }}{f_{i,\max }-f_{i,\min }} \end{equation}

The evolutionary operations (selection, crossover, and mutation) are then based on these normalized objective values.

4.1. Novel strategy in population initialization and genetic operators

In multi-objective optimization problems such as those formulated in Eq. (25), the selection of optimization algorithms critically determines robotic system performance. As a renowned multi-objective evolutionary algorithm (MOEA), the Non-dominated Sorting Genetic Algorithm II (NSGA-II) core operation involves ranking solutions into non-dominated fronts and using a crowding distance metric to promote diversity within each front, ultimately converging towards a Pareto-optimal set. This paper proposes an improved NSGA-II algorithm to resolve the multi-objective optimization problem in the serial-parallel robotic collaborative system. To enhance the global search capability of the algorithm, this paper improves the NSGA-II algorithm in two aspects: the selection of the initial population and the settings for crossover and mutation probabilities [Reference Deb, Pratap, Agarwal and Meyarivan28].

(1) Pure random initialization causes slow convergence in high-DOF systems, while pure HC risks local optima. The proposed hybrid strategy leverages the fast local search of HC and the boundary exploration of the Beta distribution. The HC method initializes the population based on Pareto theory, as illustrated in Figure 4. The detailed algorithmic process is as follows:

1. a) Generate a random point ${\boldsymbol{X}}_{\mathrm{cur}}$ within the search space $S_{\text{whole}}$ , and calculate the corresponding objective function value vector ${\boldsymbol{Y}}({\boldsymbol{X}}_{\mathrm{cur}})$ .

2. b) Generate a new point ${\boldsymbol{X}}_{\mathrm{new}}$ in the neighborhood $N({\boldsymbol{X}}_{\mathrm{cur}},d{\boldsymbol{X}})$ of point ${\boldsymbol{X}}_{\mathrm{cur}}$ , and compute the corresponding objective function value vector ${\boldsymbol{Y}}({\boldsymbol{X}}_{\mathrm{new}})$ .

3. c) Assess the relative superiority of the two points; if ${\boldsymbol{Y}}({\boldsymbol{X}}_{\mathrm{new}})$ dominates ${\boldsymbol{Y}}({\boldsymbol{X}}_{\mathrm{cur}})$ , meaning that ${\boldsymbol{X}}_{\mathrm{new}}$ is superior to ${\boldsymbol{X}}_{\mathrm{cur}}$ , then set ${\boldsymbol{X}}_{\mathrm{cur}}={\boldsymbol{X}}_{\mathrm{new}}$ ; otherwise, return to step 1 and repeat until the termination condition is met, ultimately returning the optimal solution $({\boldsymbol{X}}_{\mathrm{cur}},{\boldsymbol{Y}}({\boldsymbol{X}}_{\mathrm{cur}}))$ found.

Figure 4.

Flowchart of the Hill-Climbing algorithm for population initialization.

In the above process, $S_{\text{whole}}=\{{\boldsymbol{X}} \, | \, {\boldsymbol{X}}\in [{\boldsymbol{X}}_{\min },{\boldsymbol{X}}_{\max }]\}$ and $d{\boldsymbol{X}}$ represent the search radius. The generation methods for points ${\boldsymbol{X}}_{\mathrm{cur}}$ and ${\boldsymbol{X}}_{\mathrm{new}}$ are described by Eqs. (27) and (28), respectively. Here, $u$ is a random number uniformly distributed between [0, 1], and the termination condition is typically defined as the number of individual search iterations reaching a specified value $G_{\mathrm{HC}}$ .

(27) \begin{equation} {\boldsymbol{X}}={\boldsymbol{X}}_{\min }+u\left({\boldsymbol{X}}_{\max }-{\boldsymbol{X}}_{\min }\right) \end{equation}

(28) \begin{equation} {\boldsymbol{X}}_{\mathrm{new}}=\left({\boldsymbol{X}}_{\mathrm{cur}}-d{\boldsymbol{X}}\right)+2ud{\boldsymbol{X}} \end{equation}

(2) Mutation and Crossover Probabilities Based on the Beta Distribution:

In NSGA-II, the Beta distribution is ideal for adaptive mutation and crossover probabilities because its support is strictly [0, 1], eliminating truncation errors. Its two shape parameters enable real-time skewing of the probability mass toward either exploitation (elite fronts) or exploration (dominated fronts) based on crowding distance and rank. Thus, the Beta prior yields closed-form, fitness-driven moments that seamlessly update SBX and polynomial mutation rates without violating NSGA-II’s elitist preservation mechanism.

In genetic algorithms, the selection of crossover probability $P_{\mathrm{c}}$ and mutation probability $P_{\mathrm{m}}$ directly impacts algorithm performance. The HC algorithm can yield a high-quality population within a local scope; however, it may also guide the population into local optima. Thus, during the early stages of evolution, to preserve advantageous genes while avoiding local optima, a lower $P_{\mathrm{c}}$ and a higher $P_{\mathrm{m}}$ should be set. In the middle and later stages of evolution, an excessively high $P_{\mathrm{c}}$ can disrupt the elite individuals in the population, reducing population stability. Therefore, at this stage, it is advisable to decrease $P_{\mathrm{m}}$ while increasing $P_{\mathrm{c}}$ .

Therefore, the calculation formulas for the crossover probability $P_{\mathrm{c}i}$ and mutation probability $P_{\mathrm{m}i}$ for the population in generation $i$ are given as follows:

(29) \begin{equation} \left\{\begin{array}{l} P_{\mathrm{c}i}=\int _{0}^{0.9}f\left(x; \, t_{i}\right)dx\\[9pt] P_{\mathrm{m}i}=\int _{0}^{0.1}f\left(x; \, t_{i}\right)dx \end{array}\right. \end{equation}

(30) \begin{equation} f\left(x;\,t_{i}\right)=\frac{x^{{t_{i}}}\left(1-x\right)^{{t_{i}}}}{\int _{0}^{1}u^{{t_{i}}}\left(1-u\right)^{{t_{i}}}du} \end{equation}

where $t_{i}$ is a parameter derived from the number of iterations $N_{\mathrm{G}}$ set within the algorithm, with the relationship $t_{i}=i/N_{\mathrm{G}}$ .

1. Step 1: Utilize a random method to obtain a population $P_{0}$ with a number of individuals denoted by $N$ .

2. Step 2: Employ the aforementioned HC algorithm to search through all individuals within the population $P_{0}$ , thereby obtaining the parent population $P$ .

3. Step 3: Subject the population $P$ to genetic algorithm operations, including selection, crossover, and mutation, to generate a temporary offspring population $P_{t}$ . Combine $P$ with $P_{t}$ to form new population $P'$ .

4. Step 4: Conduct a rapid non-dominated sorting and crowding distance calculation for the individuals within the population $P'$ . Based on the elite preservation strategy, obtain the next generation population $P_{0}$ .

5. Step 5: Assess whether the termination conditions are met. If so, output the result $P_{0}$ . Otherwise, use the population $P_{0}$ as the next generation’s parent population $P$ , and proceed to Step 3.

The improved NSGA-II algorithm program flow is depicted in Figure 5.

Figure 5.

Flowchart of the improved NSAGA-II algorithm.

4.2. Algorithm performance index design

Let ${\boldsymbol{Y}}_{\mathrm{ref}}=(f_{\mathrm{ref}1},f_{\mathrm{ref}2},\cdots,f_{\text{ref }5})^{\mathrm{T}}$ be the objective function value corresponding to the custom reference point, and the objective function value vector of individual ${\boldsymbol{X}}_{i}$ is set as ${\boldsymbol{Y}}({\boldsymbol{X}}_{i})=(f_{1}({\boldsymbol{X}}_{i}),f_{2}({\boldsymbol{X}}_{i}),\cdots,f_{5}({\boldsymbol{X}}_{i}))^{\mathrm{T}}$ . $\tilde{{\boldsymbol{Y}}}_{\mathrm{ref}}$ and $\tilde{{\boldsymbol{Y}}}({\boldsymbol{X}}_{i})$ are obtained by normalizing ${\boldsymbol{Y}}_{\mathrm{ref}}$ and ${\boldsymbol{Y}}({\boldsymbol{X}}_{i})$ according to Eq. (31), and the distance $d_{fi}$ between the two is calculated, and the sum of all individuals in the population corresponding to $d_{fi}$ is calculated as $D_{f}$ according to Eq. (32), and $D_{f}$ is used as an indicator of convergence.

(31) \begin{equation} \tilde{f}_{i}=\frac{f_{i}-f_{\min }}{f_{\max }-f_{\min }}\left(i=1,2,\ldots,5\right) \end{equation}

(32) \begin{equation} D_{f}=\sum _{i=1}^{N}d_{fi}=\sum _{i=1}^{N}\sum _{j=1}^{M}\left| \tilde{f}_{j}\left({\boldsymbol{X}}_{i}\right)-\tilde{f}_{\text{ref }j}\right| \end{equation}

where $N$ is the number of individuals in the population and $M$ is the dimension of the objective function.

For solving the multi-objective minimum optimization problem, if the maximum value of each objective function is known, the vector composed of the maximum value of each objective function can be used as the reference point ${\boldsymbol{Y}}_{\mathrm{ref}}$ . In this case, the larger $D_{f}$ is, it means that the optimization result is farther away from the maximum value point, and the optimization effect is better; if the minimum value of each objective function is known, the vector composed of the minimum value of each objective function can be used as the reference point ${\boldsymbol{Y}}_{\mathrm{ref}}$ , and in this case, the smaller $D_{f}$ is, it means that the optimization result is closer to the minimum value point, and the optimization effect is better.

The convergence index $D_{f}$ can only indicate the convergence of the optimization results, and to analyze the uniform distribution of the optimization results on the Pareto front surface, the overall distribution index $C_{f}$ is introduced, which is calculated as in Eq. (33), where $C({\boldsymbol{X}}_{i})$ is the crowding degree corresponding to individual ${\boldsymbol{X}}_{i}$ . Comprehensive analysis shows that when the population convergence index $D_{f}$ is consistent, the larger $C_{f}$ indicates that the final optimization results are more evenly distributed on the Pareto front surface, and the optimization effect is better.

(33) \begin{equation} C_{f}=\sum _{i=1}^{N}C\left({\boldsymbol{X}}_{i}\right) \end{equation}

4.3. Analysis of optimal solution to multi-objective optimization problems

The multi-objective optimization problem represented by Eq. (24) is solved using the standard NSGA-II algorithm and the improved algorithm of this paper according to the parameter values listed in

Table I, and the results are shown in Figures 6 and 7.

Table I.

Optimization problem solving parameter settings.

Figure 6.

Convergence metric $D_{f}$ curve.

Figure 7.

Distribution metric $C_{f}$ curve.

The Pareto-optimal solutions presented below are the result of the optimization process employing the normalization technique described in Section 3.3 to ensure unbiased algorithm performance. All results are reported in their original physical units for clear interpretation.

It is critical to note that for complex multi-objective problems, the true Pareto-optimal set is typically unknown. The goal of MOEAs like NSGA-II is not to guarantee global optimality, which is computationally intractable for NP-hard problems, but to find a well-converged and diverse approximation of the Pareto front.

A comprehensive analysis of Figures 6 and 7 show that the $D_{f}$ metric and $C_{f}$ of the proposed improved NSGA-II algorithm stabilize faster than the standard NSGA-II algorithm after 100 generations of iterations (displaying one out of every eight points), indicating that the proposed algorithm achieves faster convergence.

Figure 8 illustrates the final optimization results (displaying one out of every eight individuals), representing the optimization outcomes for each performance metric in the collaborative assembly system. Figure 8(a)–(e) shows the distribution of the final-generation populations obtained by each algorithm across various metrics, depicting the projections of Pareto solution sets for the multi-objective optimization problem along five performance criteria.

Figure 8.

Distribution of each metric in the optimization results: a – Metric $V$ Distribution; b – Metric $\Gamma_{F}$ Distribution; c – Metric $\Gamma_{M}$ Distribution; d – Metric $K_{F}$ Distribution; e – Metric $K_{M}$ Distribution.

The improved NSGA-II method proposed in this paper demonstrates superior distribution and continuity of the Pareto front across all performance dimensions.

Table V lists a subset of the Pareto solutions obtained by the proposed algorithm for solving the multi-objective optimization problem in collaborative assembly systems. Depending on specific task requirements, appropriate performance parameters can be selected from these solutions to determine the parallel robot structure and layout parameters for designing the collaborative assembly system.

The quality of the obtained Pareto-optimal set is further evaluated using the hypervolume (HV) indicator, a standard metric that comprehensively measures both convergence and diversity. Given the high-dimensionality (5 dimensions) of the objective space, the reference point ${\boldsymbol{R}}$ was set to twice the Nadir point ( ${\boldsymbol{R}} \boldsymbol{=}\,\max (\text{Nadir point})\times 2$ ) to ensure it was strictly dominated by all solutions. The obtained HV value is 0.0902. While this absolute value may appear small, it is a direct consequence of the exponential scaling of volume in high-dimensional spaces (the curse of dimensionality) and is consistent with values reported in other many-objective optimization studies.

Table II.

The D-H parameters of the 6-DOF serial robot.