Path planning of spot welding robot based on multi-objective grey wolf algorithm

During the welding process, the robot corresponds to a time series {t₁, t₂, ⋯ , t _i , ⋯ , t_n-1, t _n  } when it moves between welding joints {γ₁, γ₂, ⋯ , γ _i , ⋯ , γ _n  }. t _i represents the movement time between welding joint γ _i and welding joint γ_i+1. t _n represents the movement time between welding joint γ _n and γ₁. The welding time includes working time and movement time. Because the working time T _work used for welding is fixed. In order to simplify the problem, the welding time mainly considers the movement time between the welding joints. All welding joints are welded, the total welding time required is shown in Equation (2). T = ∑ i = 1 n - 1 t i + T work(2)

In order to optimize welding time, the movement time of the robot between two welding joints can be shortened as much as possible according to the kinematic and inverse kinematic constraints of the robot, i.e. angular velocity and angular acceleration constraints.

The motion trajectory of the robot between two welding joints is a quintic polynomial. It is difficult to calculate the path length q_i,i+1 between two points in Cartesian space. Therefore, the joint space is used to calculate the path length. The end of robot torch is approximated as the position of welding joint. In the joint space, according to the position and attitude of the two welding joints γ _i , γ_i+1 on the path, the angular displacement corresponding is obtained and denoted by θ _i . The velocity θ ̇ i and acceleration θ ̈ i of the joint are calculated by the angular displacement. The expression of the joint angular displacement between two welding joints is expressed as: θ i = σ 5 t i 5 + σ 4 t i 4 + σ 3 t i 3 + σ 2 t i 2 + σ 1 t i + σ 0(3)

Where, σ₅, σ₄, σ₃, σ₂, σ₁, σ₀ are the coefficient. t _i is the time.

The motion trajectory between welding joints γ _i , γ_i+1 can be divided into equal parts of p-segment in joint space. Then the angular displacement θ_i1, ⋯ , θ _ip of p points are calculated on the trajectory according to Equation (3). Through the robot kinematic, the Cartesian coordinate values corresponding to the p points are obtained. With the increase of the segment numbers, the sum of p-segment trajectories is approximated as the path length q_i,i+1 between two welding joints. q i , i + 1 = ∑ j = 1 p - 1 d j , j + 1(4)

Where d_j,j+1 is the Euclidean distance between the adjacent point position (a _j , b _j , c _j ) and point position (a_j+1, b_j+1, c_j+1) of the p-segment trajectory [9].

d i , i + 1 = ( a j - a j + 1 ) 2 + ( b j - b j + 1 ) 2 + ( c j - c j + 1 ) 2(5)

All welding joints on the welding path are calculated. The total length from the start welding joint to the end welding joint follows Equation (6). L = ∑ i = 1 n - 1 q i , i + 1 + q n , 1(6)

The welding time model Equation (2) and welding path length model Equation (6) established are used as objective functions of path planning. The multi-objective model of path planning for spot welding robot is established. The formula Equation (1) is expressed as follows Equations (7)∼(8). Minimize : f 1 ( x ) = L ( x )(7) Minimize : f 2 ( x ) = T ( x )(8)

Where x∈ { x₁, x₂, ⋯ , x _i , ⋯ , x _Δ } is the X set of all welding joints sequences. x _i is a sequence of welding joints {γ₁, γ₂, ⋯ , γ _i , ⋯ , γ _n  }.

The maximum angular velocity and maximum angular acceleration are the constraints. The constraint of angular velocity is: | θ ̇ | ⩽ θ ̇ max . θ ̇ is the speed of each robot joint in the trajectory. θ ̇ max is the maximum speed that the joint cannot exceed. The constraint of angular acceleration is: | θ ̇ | ⩽ θ ̇ max . θ ̇ is the acceleration of each robot joint in the trajectory. θ ̈ max is the maximum acceleration that the joint cannot exceed.

3.1 Grey wolf optimization algorithm

GWO is a new intelligence optimization algorithm proposed by Mirjalili S in 2015 [11]. During the iterative process, the algorithm updates all positions of grey wolves according to Equations (10). D = | ξ · x p ( g ) - x ( g ) |(9) x ( g + 1 ) = x p ( g ) - A · D(10)

Where, D is the distance between the grey wolf and the prey, g is the evolutionary algebra, x _p  (g) represents the position of prey when the algorithm iterates to the g-th generation, and x (g) represents the position of grey wolf in the g-th generation.

To solve the problem of path planning, the coding method of the grey wolf is changed from any value in the continuous domain to the positive integer sequence of 1∼n. The position vector of each grey wolf is expressed as a solution to the multi-objective path optimization of spot welding robot. Let the position of the i-th grey wolf be x _i  = (x_i1, x_i2, x_i3, ⋯ , x _in ), where x_i1, x_i2, x_i3, ⋯ , x _in represents a certain sequence of welding joints {γ₁, γ₂, ⋯ , γ _i , ⋯ , γ _n  }. The multi-objective function of the path planning is defined as minF (x) ={ f₁ (x) , f₂ (x) , ⋯ , f _m  (x) }.

The operation rules are redefined as follows.

Definition 1. The minus sign in Equations (10) is defined as the operator ⊗. It represents the intersection of prey position and grey wolf position after the action of ⊙.

Definition 2. The multiplication sign of Equations (10) is defined as operator ⊙. It represents the sequence adjusted by the probability coefficient to a new sequence.

According to the definition 1 and definition 2, Equation (9) is changed to Equations (11)∼(13). D α = | ( ξ 1 ⊗ x α ( g ) ) ⊙ x ( g ) |(11) D β = | ( ξ 2 ⊗ x β ( g ) ) ⊙ x ( g ) |(12) D δ = | ( ξ 3 ⊗ x δ ( g ) ) ⊙ x ( g ) |(13)

Where, ξ₁, ξ₂, ξ₃ are the coefficient. x _α (g), x _β (g), x _δ (g) represent positions of the three best wolves (α wolf, β wolf and δ wolf) in the pack. ξ = 2 r 1(14)

Where r₁ is a random number in the range of [0,1].

According to the definition 1 and definition 2, Equation (10) is changed to Equations (15)∼(17). x 1 ( g ) = x α ( g ) ⊙ ( A 1 ⊗ D α )(15) x 2 ( g ) = x β ( g ) ⊙ ( A 2 ⊗ D β )(16) x 3 ( g ) = x δ ( g ) ⊙ ( A 3 ⊗ D δ )(17)

Where A₁, A₂, A₃ are the convergence factor and r₂ is a random number in the range [0,1] in Equation (18). According to Equation (19), φ decreases linearly from 2 to 0 as the number of iterations increases, and Max_iter is the maximum number of iterations. A = 2 φ · r 2 - φ(18) φ = 2 - 2 g Max iter(19)

All positions x (g + 1) of grey wolves update according to the positions of α wolf, β wolf and δ wolf. Operator ⊕ in Equation (20) represents the superposition effect of the sequences. x ( g + 1 ) = x 1 ( g ) ⊕ x 2 ( g ) ⊕ x 3 ( g ) 3(20)

To obtain the Pareto set, the algorithm sets the external archive set A to retain the non-dominated solutions. When the number of non-dominated solutions exceeds the capacity of archive set A, “ K -nearest” neighbor method is used to filter poor solutions. In order to enhance the population diversity in the search process, a mechanism of shift-based density estimation (SDE) is introduced [18, 22]. When the density of the surrounding areas of individuals in population estimated, the positions of other individuals are shifted according to their convergence to reflect the relative distance of individuals to the Pareto front. f m S ( x l ) = { f ′ m ( x j ) , i f f ′ m ( x l ) < f ′ m ( x j ) f ′ m ( x l ) , o t h e r w i s e(21)

Where f m S ( x l ) is the objective function of shifted. All shifted objective functions are F S ( x l ) = { f 1 S ( x l ) , f 2 S ( x l ) , ⋯ , f m S ( x l ) } . are the normalized functions of f_m (x _l ) , f_m (x _j ) [23].

The search process includes the distribution information and convergence information of individuals. The density is estimated by the difference of objective function from the individual to the “K-nearest” neighbor. Given the density estimation, the individuals with poor convergence in the population are moved to the area with high crowding degree. These individuals are given a higher density value, which makes them easier to be eliminated during the iteration process.

The complete DeMOGWO for path planning of spot welding robot is described in Algorithm 1.

This process repeats until the termination criterion is met. Finally, the solutions in archive set A are outputted.

4.1 Performance comparison on the Mo-TSP benchmark problems

Generally speaking, Mo-TSP benchmark problems [25, 26] is used to test the performance of six algorithms. The mathematical definitions and other relevant parameters of Mo-TSP are incorporated in Table 1.

The parameter settings of other five algorithms are given as follows: NSGA-II refers to [16], MOPSO refers to [17], SPEA2SDE refers to [18], IBEA refers to [19], and MOEADD refers to [20]. The experimental environment and the iteration numbers of algorithms are the same as above.

The hypervolume (HV) is a metric, where there is no prior knowledge about the true Pareto front required in the metric calculations [27]. A higher HV value indicates the better convergence as well as diversity of the points in the solution space. Given a solution set SP, the HV value of SP is defined as the area covered by SP with respect to a set

of predefined reference points SR in the objective space. λ is the Lebesgue measure.HV (SP, SR) = λ (H (SP, SR)). Due to its good performance, the HV metric is a very popular quality metric.

Since the range of the Pareto front (PF) is unknown in Mo-TSP benchmark problems, the indicator HV is considered to evaluate the performance of each algorithm. The volume of objective space between the solution set and the reference point calculated by HV is better to use a larger value. A larger HV is preferable, indicating a better approximation to the true PF and better diversity.

Run Time is the time required for the algorithm to iterate the maximum numbers. It describes computing complexity of the algorithm.

Tables 2-4 show the HV values and RunTime values of compared algorithms for different Mo-TSP problems. From the HV results shown in Tables 2-4, a phenomenon observed from these results is that NSGA-II, MOPSO and MOEADD do not provide the good solution for the Mo-TSP benchmark problems.

Tables 2, 3, and 4 show that advantage of DeMOGWO algorithm on HV is greater than other algorithms. The higher the HV value represents the better comprehensive performance of the DeMOGWO algorithm. DeMOGWO is slightly more than NSGA-II, MOPSO and IBEA with computational complexity while less than SPEA2SDE and MOEADD, this is acceptable. Moreover, as the number of cities increases, the performance differences can be better observed. It shows that DeMOGWO algorithm has strong search ability. Therefore, DeMOGWO has better performance than other five evolutionary algorithms in Mo-TSP benchmark problems.

Next, the welding joints sequence planning of BIW workpieces is will solved. In practical, a welding robot is allocated about 30 to 50 welding joints during the welding process. The paper takes 47 welding joints as an example. Figure 1 shows the welding joints assignment of a BIW workpiece.

OPEN IN VIEWER

Fig. 1

Welding joints assignment of a BIW workpiece.

The ball in Fig. 1 indicates the welding joints that need to be welded. Three welding joints are marked in the figure. Table 5 shows the coordinate values and the attitude values of 47 welding joints while welding (x, y, z are the position coordinate values, Ψ₁, Ψ₂, Ψ₃ are the attitude values). The path planning of spot welding robot can be summarized as the following steps:

Robot D-H model establishment.

Due to the multi-objective problem of path planning, the true PF cannot be obtained. Therefore, the HV and RunTime indicators are used to quantitatively evaluate the quality of the optimized solution. The parameters of the six algorithms are the same to the Section 4.1.

Table 6 shows the HV and RunTime results of the six algorithms NSGA-II, MOPSO, SPEA2SDE, IBEA, MOEADD, and DeMOGWO for the multi-objective path planning. It can be concluded from Table 6 that the HV indicators from large to small are DeMOGWO, SPEA2SDE, IBEA, MOEADD, NSGA-II, MOPSO. The larger of HV, the convergence and diversity of the algorithm is better. And the RunTime of DeMOGWO has better performance than MOEADD and SPEA2SDE.

After the initialization of the DeMOGWO algorithm, the welding sequence covering all welding joints is randomly generated. Then, the welding length model and welding time model are taken as the objective functions. Finally, the Pareto solution set is obtained by the six algorithms, as shown in Fig. 2. f1 represents the welding length (unit: 10000 mm). f2 represents the welding time (unit: 100 s). Each point in the figure represents a specific candidate solution.

OPEN IN VIEWER

Fig. 2

Pareto solution set obtained by six algorithms.

The purpose of Pareto solution set obtained through multi-objective algorithm is to provide the optimal solution for the production process. The results are helpful to improve welding efficiency and increase production capacity. In actual production, technicians can set different weights for two goals. The optimal welding path corresponding to the weight can be obtained from the solution set. For example, when the weight of welding length and welding time is set as (0.5, 0.5), the optimal welding sequence is:

44–>45–>12–>5–>28–>22–>29–> 26–> 23–> 36–>7–>17–>18–>19–>14–>10–>38–>34–>24–>20–>8–>11–>27–>6–>9–>21–>13–>25–> 39–> 37–> 16–>15–>1–> –>2–> 4–> 33–>41–>40–> 42–> 35–> 30–>31–>3–>32–>43–>47–>46. Then the path length and welding time obtained are respectively 8122.415 mm and 99.31 s.

It can be seen from the Fig. 2 and Table 6, DeMOGWO has a good diversity. One possible reason is that the density estimation mechanism enables the algorithm to maintain diversity in the iterative process. Due to the grey wolf algorithm has a strong global search ability, it can jump out of the local optimal solution through the hunting mechanism. Therefore, the convergence of the DeMOGWO algorithm is also excellent. The DeMOGWO algorithm not only coordinates the global and local search ability, but also improves the method of elite solutions obtained. This makes the non-inferior solutions can be distributed evenly in the solution space.