Time-Impact Optimal Trajectory Planning for Wafer Handling Robotic Arms Based on Improved Snake Optimization Algorithm

1. Introduction

With the rapid advancement of the semiconductor industry, wafer-handling manipu-lator arms have been confronted with elevated demands in terms of precision, flexibility, and stability [1,2,3,4]. Consequently, reliable tracing of the anticipated trajectory during operation has become a key concern in the current research. When deliberating the trajectory planning of the manipulator arm, a choice must be made between the planning schemes in Cartesian space and Joint space, which are typically determined by specific task requirements, characteristics of the manipulator arm, and tracking accuracy, among other factors. For instance, trajectory planning in Cartesian space enables direct control of the position and orientation of the manipulator end effector, which is relatively more intuitive than that in Joint space, and can more effectively plan smooth motion trajectories [5,6,7]. From the perspective of advanced computational methods, Xiao et al. [8] successfully mitigated the shock and vibration experienced by a robot during its movement from the starting position to the target position in Cartesian space, thereby ensuring the accuracy and safety of point-to-point movements. In contrast, Zhang et al. [9] successfully planned two sets of trajectories using polynomial interpolation in Joint space and verified the model by using SimMechanics. However, owing to the insufficient robustness of the PID controller against external environmental factors, the simplification of the model significantly affected the accuracy of the experiment. Additionally, the improved quantum particle swarm optimization algorithm proposed by Luo et al. [10] performed outstandingly in the trajectory planning of the robot joint space, effectively enhancing the efficiency and accuracy of trajectory planning. Regarding the S-shaped speed curve, Wu et al. [11] investigated how to plan trajectories in Joint space of industrial robots that satisfy kinematic constraints and have continuous jerk. Nevertheless, trajectory planning in Joint space requires complex inverse kinematic solutions. Owing to the presence of joint coupling, the trajectory converted to Cartesian space might be distorted or deviate from the expected path. This constitutes a potential issue for manipulator arms employed in wafer-handling scenarios.

Trajectory planning curves are vital for motion control strategies. The common types include linear interpolation, B-spline curves, polynomial curves, T-shaped speed curves, and S-shaped speed curves. Polynomial interpolation is frequently employed for this purpose. By selecting the appropriate polynomial order, its adapt ability can be flexibly adjusted to match different numbers of data points and smooth the curve, thus effectively reducing the vibration and shock of the robotic arm during movement [12,13,14]. Although high-order polynomial interpolation can provide a smoother curve fit, it may lead to the Runge phenomenon with an excessive number of data points. In addition, they are highly sensitive to minor data point variations, potentially leading to numerical instability. Furthermore, Cong et al. [15] proposed the use of fifth-order Non-uniform B-spline curves in the trajectory planning and crafted a multi-objective trajectory optimization model employing the Non-dominated Sorting Genetic Algorithm-II (NSGA2). Regarding the T-shaped speed curve, Ayazbay et al. [16] applied a mini-mum time-impact planning algorithm to optimize the curve, targeting an enhanced dynamic performance and reduced mechanical wear. However, this curve lacked smoothness at the endpoints, necessitating additional processing to ensure the overall continuity and smoothness of the trajectory. In response to these limitations, researchers have proposed S-shaped speed curves [17,18] as trajectory planning curves, which are distinguished by smooth transitions at the velocity inflection points. This results in less impact force and noise compared to the T-shaped acceleration and deceleration phases. S-shaped speed curves are commonly used in applications requiring precise control and minimization of mechanical vibration.

Heuristic algorithms are increasingly employed in trajectory planning and typically address complex optimization problems when it is difficult to find an exact solution or when the cost of computing an exact solution is too high. Intelligent algorithms, such as genetic algorithms (GA), particle swarm optimization (PSO) algorithms, annealing algorithms (SA), whale optimization algorithms (WOA), and grey wolf optimization (GWO) algorithms provide novel solutions for addressing intricate trajectory planning problems. They have the characteristics of flexibility, robustness, and diversity, which can help users find the global optimal solutions. Among them, the PSO algorithm [19], which simulates the social behavior of bird flocks to achieve optimization, is widely used to optimize the use of robotic arms owing to its simplicity, ease of implementation, and strong ability for multi-modal optimization. However, the PSO algorithm cannot avoid problems of parameter sensitivity and convergence acceleration. Therefore, research teams have improved to varying degrees to overcome these shortcomings [20,21,22,23]. Lu and Wang [24] built an optimal trajectory planning method based on a B-spline and a WOA to improve the work efficiency of a robotic arm. However, they did not fully compare the performance differences between the WOA and other existing algorithms. By combining the back propagation neural network (BPNN) with the GA, the overall performance of the upper limb rehabilitation robotic arm trajectory planning was improved to some extent. However, the parameter adjustment after the combination was very complex and required a large number of experiments and professional knowledge for correction [25]. The SO algorithm excels outstanding in the realm of optimization owing to its exceptional flexibility and robustness. It is not only easy to integrate chaos theory and learning mechanisms to bolster search efficiency but also to maintain efficient exploration ability under poor initial conditions. In addition, it can rapidly adapt to the changes in dynamic environments and optimize search strategies in real time, thereby safeguarding the continuity and effectiveness of the optimization process. Li et al. used a SO algorithm to find optimal paths in dynamic environments. Although its performance has been improved compared with other intelligent algorithms, the improved algorithm still requires careful parameter tuning [26]. To address the shortcomings of the SO algorithm, recent literatures [27,28] begun to focus on solving its core mechanism to raise the overall performance of optimization problems. The ISO has also been gradually applied to path planning and to solve general constrained optimization problems.

Therefore, considering the shortcomings of the afore-mentioned studies and the main factors that should be prioritized during the actual operation of the wafer manipulator, this study proposes an improved control strategy for a manipulator that combines the single-segment start-stop S-shaped speed of the circular arc transition, an ISO algorithm, and a time-impact optimization method of trajectory planning. It is worth noting that it has more advantages in terms of the stability of the control effect or the convergence ability of the algorithm compared with other types of methods. And it can adapt to parameter optimization in different scenarios. In summary, the contributions of this study are as follows.

(1)

In the area of intelligent optimization methods, the following improvements have been made to the SO algorithm.

●

Tent Chaotic Map initialize populations, instead of randomly generating populations with ‘ rand ’ functions.

●

Dynamic learning factors, integrated with random perturbations, replace static learning factors.

●

A cosine annealing learning rate scheme is employed for the adaptive updating of individual positions.

●

The Bayesian method is utilized to optimize thresholds, surpassing the reliance on empirical thresholds.

(2)

In the area of trajectory planning, the following improvements have been made to the S-shaped speed curve.

●

The interpolated arc transition supplants the abrupt right-angle transition

●

The continuous single-segment start-stop mechanism supersedes the segmented start-stop mechanism.

2. Materials and Methods

2.1. The Structure Parameters of the Wafer-Handling Robotic Arms

Assuming that the task requirements can be realized in a plane, the system model is simplified when the robotic arm transports the wafer, as shown in Figure 1. Before planning the trajectory, it is necessary to establish the link coordinate system correctly. In general, the points along the Z_i axis are aligned with the i-th joint axis. Furthermore, the points along the X_i axis are aligned with a common vertical line. In addition, the Y_i axis is oriented according to the right-hand rule. In this study, the parameters of the linkage rods are expressed using the improved Denavit-Hartenberg (DH) method. The parameters of each joint and link are obtained using the link coordinate system as a reference, as shown in Table 1.

To facilitate the visual analysis of the subsequent trajectory planning and verify the correctness of the simplified model, this study utilizes the MATLAB robot toolbox developed by Peter Corke. It is a powerful robotic simulation and analysis tool that provides a range of capabilities for modeling, simulating, controlling, and visualizing robots. A robotic arm consistent with the simplified model is successfully constructed on the interface by following the building instructions in the MATLAB robot toolbox, as shown in Figure 2.

2.2. The Mathematical Model of Multiple Objective Function

In recent years, the optimization of time, energy, impact, and other indicators has been shown to improve the operational efficiency and motion performance of robotic arms in the field of robot trajectory planning [29,30,31,32]. Considering the production background and optimization effect, a multi-segment point-to-point trajectory in Cartesian space is adopted to achieve multi-objective optimization of time and impact. First, according to the expected trajectory of the robotic arm, it must be discretized into several interpolation points in Cartesian space. Subsequently, several discrete points are connected using the trajectory function. The joint angle corresponding to the entire trajectory in the joint space can be calculated using the inverse kinematics equation of the robot. Finally, only the selected joint angle conforming to the motion trajectory is substituted into the forward kinematics of the robot to solve the cartesian coordinate value of the expected trajectory. It is assumed that the time at which the manipulator reaches each interpolation point is t_i . The time interval between adjacent interpolation points is h_i = t_i+1 – t_i , where i is the index of the interpolation points (i = 0,1...n). At the same time, through the preliminary analysis, it is clearly observed that there is an inverse correlation between time and impact. When the time of the robotic arm increases, the impact decreases. Therefore, it is impossible to simultaneously obtain the minimum values. Thus, considering the constraints of the geometric joints of the robotic arm, velocity, acceleration, and jerk, the primary goal of this study is to achieve an optimal balance between time and impact. Therefore, a Pareto optimal method is adopted in this study to achieve the tradeoff between optimization objectives. The time function is given by Equation (1). Because the impact quantity can not be directly measured, it is quantified by the acceleration, as shown in Equation (2). The constraints are expressed in Equations (3). The multi-objective optimization function is shown in Equation (4).

(1)

(2)

(3)

(4)

where F₁ , F_1min , and w₁ are the objective function, minimum value, and weight value of time, respectively. F₂ , F_2min , and w₂ are the objective function, minimum value, and weight of impact, respectively. h_i is the time interval between the i-th and i+1-th interpolation point. t_i and j_i denote the time node and jerk value of the i-th interpolation point. θ_max , v_max , a_max , and j_max are the maximum rotational angle, speed, acceleration, and jerk of the joint of the manipulator arm, respectively.

2.3. The Construction of Trajectory Planning

2.3.1. The Explanation of S-Shaped Speed Curve

The T-shaped speed curve is a widely used method for robotic arm control. Because of its simplicity and ease of control, it is extremely common in industrial applications. As shown in Figure 3, this curve may produce large impacts and vibrations at points A, B, C, and D, potentially causing the wafer to break. By contrast, the S-shaped speed curve uses nonlinear adjustments to increase and decrease the velocity, making the motion process smoother and effectively solving the discontinuity problem of the acceleration of the T-shaped speed curve. Therefore, the S-shaped speed curve is more suitable in scenarios requiring higher operational stability.

The implementation complexity of the S-shaped speed curve is higher than the T-shaped speed curve. However, it effectively reduces impact and improves the quality of motion. The S-shaped speed curve divides the entire motion process into seven stages: accelerating acceleration phase, uniform acceleration phase, decelerating acceleration phase, uniform speed phase, accelerating deceleration phase, uniform deceleration phase, and decelerating deceleration phase.

In this study, the kinematic characteristics of an S-shaped speed curve at different stages are described in detail. Specifically, the position, velocity, acceleration, and jerk functions of each stage are shown by Equations (5), (6), (7), and (8), respectively.

(5)

(6)

(7)

(8)

where t_k (k = 0, 1...7) represents the transition point time of part k, which indicates the specific point in time at which each stage begins or ends. τ_k (k = 0, 1...7) represents the local time period of part k, which refers to the difference from the starting point of the current part, called time zero, to the time point of the current part. The equation used τ_k = t_k - t_k-1 (k = 1, 2...7). s_k (k = 1, 2...6) represents the last shift of part k. v_k (k = 1, 2...6) represents the final speed of part k and v_s represents the initial speed. v_max is maximum speed. a_max is the maximum acceleration. j_max is the maximum acceleration.

2.3.2. The Improvement of S-Shaped Speed Curve

Trajectory planning is a crucial step in the robotic motion control. When the wafer-handling robotic arm is in operation, it frequently encounters turning points where the direction of movement requires alteration. These turning points often serve as a source of abrupt jerk changes. To alleviate this problem, an S-shaped speed curve adopts a circular arc transition strategy to guarantee a smooth trajectory transition. This strategy not only mitigates the impact and vibration induced by the robotic arm but also contributes to shortening the execution time of the motion trajectory. Prior to constructing the S-shaped speed curve, it is imperative to determine the key interpolation points of the trajectory in Cartesian space, as indicated in Table 2. Simultaneously, the principal motion parameters must be set as presented in Table 3. v_ p₀, v_ p₁, and v_ p₂ are the velocities of the initial point, the circular arc transition section, and the terminal point, respectively. v_max , a_max , and j_max denote the maximum velocity, the maximum acceleration, and the maximum rate of change in the acceleration of an end effector in Cartesian space, respectively. r is the radius of the circular arc transition at the inflection point.

The trajectory curve shown in Figure 4. can be generated by applying the parameters listed in Table 3. to the S-shaped speed curve model. In this figure, the black lines represent the direct path without a circular transition, whereas the blue lines show a smooth path with a circular transition. It is noteworthy that the radius of the circular transition can be adjusted according to specific work requirements. Overall, the circular transition not only improves the smoothness and efficiency of the robotic arm's motion but also helps to protect the robotic arm and its payload, while enhancing the overall operational performance. This transition method allows the robotic arm to perform complex tasks more accurately and reliably.

By analyzing the data in Table 2, it is clear that this trajectory planning task involves three key nodes that the robotic arm must reach in sequence. S-shaped speed curve planning is a standard point-to-point start-stop motion scheme. Specifically, this task is represented by two movements from the starting point to the middle point and then from the middle point to the end point, with zero starting and ending speeds for each movement. Based on the data in Table 3, the velocity image of the scheme is shown in Figure 5. However, this scheme causes the manipulator to frequently start and stop, which not only increases the vibration during the motion process but also leads to time wastage in the acceleration and deceleration stages. To solve this problem, this study has proposed an improved scheme for the single-segment start-stop mechanism, optimizing the starting and stopping parts of the intermediate process into a smooth motion mode with constant speeds at both ends, as shown in Figure 6. When the velocity of the first segment of motion attains the maximum value of the set parameters, instead of reducing the velocity back to zero as in the standard scheme, the second segment of motion proceeds in a uniform transition manner. The specific parameters still refer to Table 3. By using this improved scheme, the problems existing in the conventional scheme can be effectively reduced. Simultaneously, the movement efficiency can be improved and the vibration can be reduced.

Comparing Figure 5. and 6, it can be seen that the improved motion scheme results in a smoother speed curve. Concurrently, it effectively reduces the vibration generated by the motor during starting and stopping. With the basic parameters unchanged, the unoptimized time schedule which does not include arc adjustment and smooth transition took 7.26 s, whereas the optimized time schedule, which includes arc adjustment and smooth transition, is shortened to 5.51 s. After the two parts are optimized, the total time is saved by 1.75 s, a 24.1% reduction compared to the original scheme.

2.4. Improved Snake Optimization Algorithm

2.4.1. Initialization Population of Tent Chaotic Map

The snake optimization (SO) algorithm is an intelligent optimization technique that simulates the specific mating behavior of snakes. It achieves an effective balance between global and local searches by using a unique search mechanism. In the initialization phase of the algorithm, individuals in the population are equally divided into male and female individuals. The composition of this initial population can be represented by Equation (9), where the numbers of male and female individuals are represented by Equations (10) and (11), respectively.

(9)

(10)

(11)

where X_i denotes the position of the i-th individual. The random function ‘ rand ’ generates a real random number between 0 and 1. The upper and lower bounds of X_i are ub and lb, respectively. N is the total population size. N_m and N_f are the numbers of individual males and females, respectively.

However, the use of the function ‘ rand ’ in Equation (9) of the SO algorithm may result in an uneven distribution of random numbers. Moreover, because this function is based on a deterministic algorithm, the same random sequence is generated when the same seed is given. This situation may lead to the premature convergence of the heuristic algorithm to the local optimal solution, making it challenging to find the global optimal solution. To solve this problem in population initialization, this study introduces a Tent Chaotic Map. It has excellent chaotic characteristics, which can effectively improve the uniformity of population distribution. In addition, compared with other Chaotic Maps, its implementation is more intuitive and simpler. The Tent Chaos Map can be expressed as Equation (12). And ‘ x_i ’ replaces ‘ rand ’ as shown in Equation (13).

(12)

(13)

where x_i and x_i+1 stand for the i-th and i+1-th values of Tent Chaotic Map, respectively. X_i stands for the i-th individual position. Moreover, parameter a is the variable value, which is 0.5 in this study [33].

2.4.2. Dynamic Learning Factor with Random Noise

In heuristic algorithms, learning factors play a key role in regulating the search behavior, including the balance between exploration and exploitation. Exploration refers to the ability to search for new solutions randomly, whereas development refers to the ability to search locally for the best solution. The learning factors of the SO algorithm include the constant value C₁ of the food calculation equation, constant value C₂ of the food exploration stage, and constant value C₃ of the food development stage. The initial values of these constants are typically set to C₁ = 0.5, C₂ = 0.05, and C₃ = 2. The equations containing C₁, C₂, and C₃ are Equations (14), (15), (16), (17), and (18). The subscript ‘m’ and ‘f’ refers to the male and female snakes, respectively. The special mention of female snake expressions will be omitted in the following text.

(14)

(15)

(16)

(17)

(18)

where Q is the amount of food consumed. t denotes the current number of iterations. T denotes the maximum number of iterations. Temp is the iteration temperature. X_{i, m_1} , X_{i, m_2} , X_{i, m_3} , and X_{i, m_4} are the locations of food in the exploration phase of male snakes, food in the development phase of male snakes, combat mode in the development phase of male snakes, and mating mode in the development phase of male snakes, respectively. X_{rand,m_1} is the randomly selected male population. X_{best_ f} is the best position for the females in this population. A_m, F_m, and M_m are measured for the male's ability to find food, the male's ability to fight, and the male's ability to mate.

Because the learning factors C₁, C₂, and C₃ remain unchanged during the iteration, the algorithm may not adapt to the dynamic changes in the solution space. Hence, the delta value is introduced for the normal distribution of random tubs. This perturbation helps the SO algorithm eliminate the trap of local optimal solutions, thus enhancing its global search abilities. By dynamically adjusting these fixed learning factors C₁, C₂, and C₃, the SO algorithm can achieve an effective dynamic balance between global and local searches, thereby improving the performance and efficiency of the algorithm. The concrete expression of the random disturbance is shown in Equations (19) and (20). In Equation (19), alpha is used to store the parameter values for learning factors C₁, C₂, and C₃.

(19)

(20)

where sigma represents the amplitude of the perturbation. It is 0.1 in this paper [34]. alpha is limited to a range between 0.01 and 1. Random function ‘ randn ’ represents a random value with a standard normal distribution.

2.4.3. Cosine Annealing Learning Rate

In Equations (15), (16), (17), and (18), the strategy for updating the individual positions is relatively simple. Without considering adjusting the strategy as an iteration progress, the SO algorithm mainly relies on the current individual's position and velocity vector to determine the next individual's position. This method is more likely to cause the algorithm to become stuck in a locally optimal solution. To improve the update strategy of the original SO algorithm, this study introduced a cosine learning annealing rate. This strategy simulates the dynamic adjustment of the learning rate based on the cosine function, which can dynamically adjust the speed of the position updating based on the number of iterations. In the early stages of the iteration, the learning rate can be reduced rapidly to achieve fast convergence. In the later stages of the iteration, fine adjustments are made to optimize the precision of the solution by reducing the learning rate, as shown in Figure 7. The specific expression of the cosine learning annealing rate is shown in Equation (21).

(21)

where η_t is the current learning rate. η_min and η_max are the minimum and maximum values of the learning rate, which are 0.001 and 0.01 in this paper. T_cur and T_max are the current and total iterations, respectively. T_max is usually set to half the number of iterations of the main algorithm, which is 20 in this article [35].

The cosine annealing learning rate exhibits periodic change, which enables the SO algorithm to adopt diverse search strategies in different stages, thus enhancing the dynamic adjustment ability of the algorithm. This flexibility helps the algorithm to avoid premature convergence during the search process, thereby increasing the possibility of finding a global optimal solution. Therefore, by introducing the cosine annealing learning rate based on Equations (15), (16), (17), and (18), the improved Equations (22), (23), (24), and (25) are obtained.

(22)

(23)

(24)

(25)

2.4.4. Cosine Annealing Learning Rate

In the process of determining threshold values of the SO algorithm, parameter setting often depends on the experience and intuition of algorithm developers or needs to be determined by practical users through practice. This approach lacks universality and requires considerable time for experimentation and adjustments. If the parameters are not selected properly, the adaptability of the algorithm may be affected. To solve the problem in which the SO algorithm has difficulty adapting parameters in the iterative process, this study introduces a hyperparameter tuning method commonly used in machine learning, which is known as Bayesian optimization. The Gaussian process mainly includes Gaussian process regression and an acquisition function.

First, it establishes prior knowledge of the global behavior of the objective function through a Gaussian process. The mean vector and covariance matrix involved in the prior distribution are shown in Equations (26) and (27). Subsequently, using the observed objective function to obtain values at different input points, this prior knowledge is updated to form a posterior distribution. The posterior mean and variance are given by Equations (28) and (29). Then, the improved collection function is expected to consider the current best observations and the posterior distribution of the Gaussian process to determine the next evaluation point. The collection function used in this study is expected to improve equation shown in (30). In short, this method of guiding the search based on existing knowledge can avoid many useless attempts.

(26)

(27)

We assume that the above expression follows a multi-variate normal distribution, where μ₀(x) represents the mean function.∑₀(x,x') represents the kernel function, also known as the covariance function, which describes the similarity between points x and x' in the input space.

(28)

(29)

(30)

where K_x1:n,x1:n are the kernel matrix. σ² is the noise variance. f(x₁:n) is the observed vector. f *= max_{i = 1, ... n} is the currently observed best target value. f(x_i) is the function value.

2.4.5. ISO Algorithm Steps

The ISO algorithm is outlined in Algorithm 1, as shown in Table 4.

The flow chart of the ISO algorithm, which illustrates the step-by-step process, is shown in Figure 8.

3. Results

3.1. Ability Test and Analysis for ISO

3.1.1. Parameter Settings and Simulation Results

In this study, all algorithm experiments and simulations are performed using MATLAB 2022B software platform. To verify the convergence performance and accuracy of the ISO algorithm, five heuristic algorithms—SO [36], SA [37], PSO [19], GWO [38], and WOA [39]—are added for comparison. In addition, to ensure the rigor and accuracy of the experimental results, the CEC2005 benchmark test set is selected as the objective function, as listed in Table 5. Five unimodal functions (F₁-F₅) and three multi-modal functions (F₆-F₈) are listed in Table 5. Through an in-depth analysis of these benchmark functions, the performance of the ISO algorithm is evaluated to a certain extent.

Before simulating the benchmark test functions for the six heuristic algorithms, a brief introduction to their main parameter settings is necessary to ensure that researchers can obtain similar experimental results when conducting simulations. To ensure the fairness of experiments, the population size (N) of all algorithms is set to 50, the maximum iteration number (T_max) is set to 100, the dimension (Dim) is set to 30, and each algorithm is run 20 times independently. For the other key parameters within the algorithms, recommended values based on experience are used for setting, and the specific parameter values are listed in Table 6. After successfully configuring these parameters, a MATLAB script is used to start the six algorithms simultaneously. The experimental results are summarized in Table 7, and the iteration curves of the benchmark test functions are displayed in Figure 9.

3.1.2. Simulation Result Analysis

The results of the eight benchmark functions listed in Table 7 are comprehensively analyzed and compared. As shown in Figure 9(a)-(e), it is clearly observed that when the dimension (Dim) is 30, the ISO algorithm achieves an optimal fitness value of 0 for the F₁ and F₃ problems in the unimodal test function. Its convergence curve can rapidly jump out and continue to decline after falling into the local optimum for a short time. This indicates that the ISO algorithm can effectively escape from a local optimal solution and successfully obtain a global optimal solution. By contrast, to a certain extent, the other five algorithms fall into the local optimum. For the F₂ and F₄ problems, although the average value of the ISO algorithm is not zero, the variance is zero, indicating that the algorithm can obtain a global optimal solution to some extent. Its average precision is better than that of the second ranked GWO algorithm. In addition, the convergence speed of the ISO algorithm is faster than those of the other algorithms. For the F₅ problem, the performances of the ISO algorithm and the other five algorithms are similar, but its overall optimization precision is slightly higher, and the convergence curve is better, without falling into the local optimum too early. Overall, the optimization effect of the ISO algorithm on unimodal functions is significantly better than that of the other five algorithms. Further analyzing Figure 9(f)-(h), when the dimension (Dim) is 30 in the multi-modal test functions, the SO, GWO, WOA, and ISO algorithms can achieve the global optimal value in the F₆ and F₇ problems with good overall performance, but the ISO algorithm has a faster iteration speed. However, the performance of the six algorithms is not outstanding in the F₈ problem, but the average value and variance of the ISO algorithm are the best, with an optimization precision approximately twice that of the other algorithms. In summary, the ISO algorithm demonstrates higher optimization precision in both unimodal and multi-modal functions than the other algorithms. It also improves the avoidance of falling into the local optimum compared to the SO algorithm, showing the expected optimization effect.

3.1.3. Friedman Test

Further analysis of the results from 20 experimental groups using the Friedman test method reveal significant differences in the overall performance of the six algorithms across eight benchmark functions. It is a flexible and effective tool that does not depend on the specific distribution of data, making it particularly suitable for comparing the algorithm performance. When dimension (Dim) is set to 30, the average ranking of the mean and variance of each algorithm is shown in Table 8. A lower average ranking indicates that the algorithm performed better, and the results are more stable. According to the data in Table 8, the ISO algorithm ranks first in terms of the mean value and variance. Compared with the GWO algorithm, which ranked second, the average ranking score of the ISO algorithm's mean and variance are 38.9% higher. Combined with the SO algorithm, the ISO algorithm shows a significant improvement in the average ranking score of the mean and variance, reaching 60.8% and 63.4%, respectively. Considering the overall performance of the algorithm, the ISO algorithm exhibits better performance in terms of stability and convergence accuracy.

3.2. Time-Jerk Optimal Trajectory Planning Based on ISO

After elaborating on the ISO algorithm and the modified S-shaped speed curve in detail, a combination of the two is performed. The objective function is obtained in this study, as presented in Equation (4). Firstly, Equation (4) adopts a compromise algorithm by adjusting the values of w₁ and w₂ to determine the equilibrium point in the multi-objective optimization. Because the weights of time and impact are considered consistent in this study, both w₁ and w₂ are set to 0.5. Secondly, the time function F₁ and the impact function F₂ belong to different categories of physical quantities. When the weighted sum of their calculation results is taken directly, it leads to an imbalance between the time and impact. To address this issue, a minimum-maximum normalization method is employed in this study. It is typically used to eliminate the influence of dimensions among different features, making the data comparable on the same scale, as shown in Equations (31) and (32).

(31)

(32)

After the above operations, the domains of the values for F₁ and F₂ functions are restricted within the same range. Based on this, a set of excellent solutions is selected by solving a multi-objective function. Then, the optimal solution is identified using the Pareto optimization method. The Pareto judgment standard is as follows. If a new solution improves at least one objective without causing a performance decline in other objectives, then the new solution is considered to be superior to the original solution. When these Pareto optimal solutions are smoothly connected, a Pareto frontier is formed. Ultimately, this study aims to determine a balance point on this curve to achieve the best overall effect.

The specific parameter configuration of the simulation experiment can be found in Table 9 for the application of the ISO algorithm to the S-shaped speed curve optimization problem, and the kinematic parameters are set according to Industrial Robot ISO 9283:1998. In the simulation experiments, a_max and j_max are set as independent variables and their upper bound (ub) and lower bound (lb) are between 80 and 120. In the SO algorithm, the empirical values of threshold parameters threshold₁, threshold₂, and threshold₃ are 0.25, 0.6, and 0.6, respectively. To determine more reasonable parameter combinations, the ranges of these thresholds are unified and adjusted. The range of threshold₁ is limited to 0.05 to 0.45, while the ranges of threshold₂ and threshold₃ are adjusted to 0.4 to 0.8. Correspondingly, the search range of the threshold parameter (Bayesian threshold) in the Bayesian optimization process is also set to 0.05 to 0.8 to facilitate determining the optimal parameter configuration. The other parameters have the same meanings as those described earlier. The Bayesian optimization threshold simulation results are shown in Table 10, which includes the iteration, the actual observed objective value of the current iteration, the best observed value among all observations (BestSoFar_observed), the model-predicted best value (BestSoFar_estimated), and the best thresholds combination obtained in this generation. However, because there are no enough data to support the construction of the first-generation model of Bayesian optimization, Objective, BestSoFar_estimated, and BestSoFar_observed are not included in the calculation. In addition, a simulation diagram of the minimum objective value and function calculation times corresponding to Bayesian optimization is shown in Figure 10, where the initial value of the first generation is replaced by number 1. For the dual-objective optimization problem considered in this study, the distribution of the Pareto frontier is generally in the form of a curve, which is intuitively shown in Figure 11.

According to the data provided in Table 10, Bayesian optimization is used to select the second-generation simulation results with the minimum objective function value. This selection is based on minimization of the target value during the iteration. Figure 10. clearly shows the dynamic change in the target value with the iteration, making the change process clear. Further analysis of Figure 11, the minimum values of F₁ and F₂ respectively reach 5.6294 s and 115.937 mm/s. And the two minimum values are respectively located at the two ends of the Pareto frontier composed of the red dots. To optimize both the time and impact, an optimal advantage must be found on the Pareto frontier. A specific point is identified in Figure 11, which achieves an optimal balance between the time and impact. When a_max and j_max are 78.1476 mm/s² and 119.8756 mm/s³, the time and impact of the Pareto optimal value are 5.8028 s and 167.5427 mm/s, respectively.

Using the computed a_max and j_max, these data are reintegrated into the trajectory planning of the S-shaped speed curve. Through this process, it is possible to accurately determine the partial displacement of the trajectory planning in the directions of the X-axis and Y-axis, as well as the associated combined velocities, combined accelerations, and combined accelerations, as shown in Figure 12. These figures clearly show that the planned S-shaped speed curve exhibited smooth characteristics during movement, effectively avoiding sudden changes and shocks. This is extremely important for precision operations such as wafer handling, as it helps reduce the damage that can be caused by impact.

In addition, this study meticulously records the key parameters, particularly the joint angles, during the S-shaped speed curve movement. Two sets of joint angles consisting of several interpolation points are substituted into the forward kinematics formula, as shown in Equation (33). X and Y represent the coordinate values in the Cartesian coordinate system, respectively.

(33)

The actual motion trajectory depicted in the planar space is shown in Figure 13. By comparing Figure 4 and 13, the high consistency of the two trajectories can be clearly verified. This comparison clearly validates the accuracy and reliability of the ISO. This planning method not only elevates the accuracy of trajectory planning, but also ensures the safety and reliability of the motion.

4. Conclusions

This study proposes an improved snake optimization (ISO) algorithm. Meanwhile, a single-segment start-stop S-shaped speed curve with arc transition is constructed based on ISO, which is a method specifically designed for the rapid and stable operation requirements of the wafer-handling robotic arm. This method has conducted in-depth research and optimization on the key issues of time and impact within the context of robotic arm trajectory planning. Through a comparative analysis with the original S-shaped speed curve scheme under the same parameter settings, the new method improves the time efficiency by 24.1%, thereby sig nificantly reducing the time required for the robotic arm to complete the task. Furthermore, the performance of the ISO is compared with that of SO, SA, PSO, GWO, and WOA on a series of benchmark test functions. The results indicate that, in terms of the two key performance indicators, —average value and variance, —the ISO algorithm outperforms the second ranked GWO algorithm by a margin of 38.9%. These improved simulations fully demonstrate the effectiveness and superiority of the proposed method in addressing the trajectory planning problem for wafer-handling robotic arms.

In this study, the shortcomings of the SO algorithm have been extensively studied and improved. Although some improvements have been made in searching for the optimal solution, there is still a risk of falling into local optima when dealing with some uni modal and multi-modal problems. Additionally, although Bayesian optimization improves the search efficiency, it consumes more computer resources. Therefore, the follow-up work of this study will focus on further adjusting the search mechanism of the SO algorithm to enhance its convergence abilities. Simultaneously, we plan to explore the combination of Bayesian optimization and bandit-based methods, aiming to reduce computational costs while improving the optimization efficiency.