Disassembly sequence planning based on a modified grey wolf optimizer

Disassembly sequence planning (DSP) can effectively increase the disassembly efficiency, shorten the disassembly cycle, reduce disassembly costs, and reduce environmental hazards of end-of-life (EOL) products, playing an important role in manufacturing industries. Thus, it is urgent to propose an approach to solve the DSP problem. DSP is a famous NP-hard combinatorial optimization problem. As the size of components increases, exact algorithms can hardly obtain the optimal disassembly sequence. Therefore, we propose a promising intelligence algorithm, modified grey wolf optimizer (MGWO), to solve the DSP problem. MGWO inherits the main idea of the hierarchy and hunting mechanism of the original grey wolf optimizer (GWO). Three new operators are designed in MGWO to ensure the feasibility of solutions under the complex constraint of disassembly precedence. The feasible solution generator (FSG) is designed to obtain feasible disassembly sequences, the neighborhood search operator (NSO) is developed to make wolves (solutions) self-evolving, and the guided search operator (GSO) is used to make the wolf group guided by three leaders of wolves. Two engineering cases are applied to validate the effectiveness of the proposed operators. Then, they and two real-world applications are used to compare the MGWO with other reported methods. The results demonstrate that MGWO can solve the DSP problem effectively.


Introduction
With the renewal of modern technology and the shortening of product life, an enormous amount of end-of-life (EOL) products has been generated. Minimizing the adverse environmental effect of EOL products and increasing their recovery rate become a global trend. Disassembly, a systematic method to divide products into components and subassemblies, is the first and often the most challenging process of remanufacturing EOL products [1]. A good disassembly sequence can effectively improve the disassembly efficiency, shorten the disassembly cycle, and reduce the disassembly cost. Hence, the key of the disassembly process is to obtain the optimal disassembly sequence.
Disassembly sequence planning (DSP) studies the disassembly sequence of components for a given EOL product to minimize disassembly time and cost, maintain stability during the disassembly process, and improve recovery efficiency. It has been proved to be an NP-hard problem [2]. The common methods for solving DSP include modeling, mathematical programming, and intelligent optimization algorithms. The modeling methods contain AND/OR graphs, Petri nets, and matrix-based models. Lambert [3] used an AND/OR graph to represent the disassembly process and designed a binary integer linear programming approach based on the AND/OR graph to solve DSP. Cappelli et al. [4] proposed an AND/OR graph of mechanical assemblies to represent the solution space of disassembly sequences. Rai et al. [5] combined a Petri net with heuristic search procedures to solve DSP. The Petri net guides the disassembly process through a token game. Kuo [6] proposed a Petri netbased analysis approach to solve DSP. The Petri net determined the optimal tradeoff between the cost and the environmental effectiveness of the disassembly processes. Li et al. [7] designed an interference matrix to represent the geometric constraints between components. The feasible disassembly sequence was generated based on the interference matrix. In addition to these modeling methods, many linear programming methods are also used to solve DSP. Zhu et al. [8] introduced a machinereadable disassembly information model and proposed a linear programming-based optimization method to solve the method. Ma et al. [9] designed an integer programming method to solve DSP in a parallel disassembly environment to maximize profit. Behdad et al. [10] developed a stochastic mixed-integer nonlinear programming method to solve DSP considering uncertain outcomes, such as time, cost, and the probability of causing damage. Kim and Lee [11] developed an integer programming model to represent disassembly sequences and proposed a branch and bound method for obtaining the lower and upper bounds of DSP. The limitation of modeling and mathematical programming is that they can only solve small-scale DSP problems. With the expansion of the scale of problems, the computational time of modeling and mathematical programming increases exponentially. The high computational time is unacceptable. Hence, the intelligent optimization algorithm is the most widely used method for solving DSP. Wang et al. [12] proposed a genetic algorithm (GA) to solve DSP. Their crossover and mutation operators can effectively help the algorithm to find the optimal disassembly sequence. Go et al. [13] designed an optimization model for DSP and presented a GA to solve it. Yeh [14] proposed a simplified swarm optimization (SSO) algorithm for solving DSP. His algorithm combined the precedence preservative operator, feasible solution generator, repetitive pairwise exchange, and self-adaptive parameter control procedures. Then, Yeh [15] modified the SSO algorithm by improving its update mechanism and revising its self-adaptive parameter control procedure to make it better solve DSP. Percoco and Diella [16] used an artificial bee colony (ABC) algorithm for DSP considering the disassembly cost and the environmental impact. Tian et al. [17] presented an ABC algorithm to solve DSP considering the stochastic cost and time of product disassembly. Additionally, many intelligent optimization algorithms, including GA [18][19][20], particle swarm optimization (PSO) [21,22], ant colony optimization (ACO) [23,24], teaching-learning-based optimization (TLBO) [25,26], and Bees algorithm (BA) [27,28], have also been designed for solving DSP. The limitation of intelligent optimization algorithms is that their accuracy is not exact. The intelligent optimization algorithms can only obtain an approximate optimal solution, but cannot guarantee to achieve an exact solution. Therefore, how to design an intelligent optimization algorithm to balance solution quality and computational efficiency has become the focus of our research. Besides, the parameter settings in most of the above algorithms need constant tunings, such as acceleration constant and inertia weight in PSO, the probability of crossover and mutation in GA, and pheromone evaporation rate in ACO. The characteristic makes it difficult for above algorithms to adapt to a variety of EOL products. Based on this, we design an algorithm with few parameters, which are adaptive and robust for various situations. It is an obvious advantage of our algorithm compared with previously reported methods.
Grey wolf optimizer (GWO) is a novel intelligence optimization method designed by Mirjalili et al. [29], inspired by mimicking the hunting mechanism and the social hierarchy of the wolf group. Due to its high convergence speed and simple implementation, GWO has been applied in many industrial fields. For example, Mohanty et al. [30] used GWO to implement a maximum power point tracking design for maximizing power extraction in a photovoltaic system. Jayakumar et al. [31] proposed GWO to allocate generation and heat outputs to the committed units in a power system operation. Their proposed method maximized power output while minimizing cost. Lu et al. [32] applied a multi-objective GWO to the welding shop scheduling. Besides, GWO was also used in engineering design [33], path planning [34], and power scheduling [35], etc. Unlike many meta-heuristic algorithms, GWO does not require many parameters to be tuned. Hence, it is very convenient to use GWO to solve optimization problems.
However, GWO is designed for continuous optimization problems, so it cannot be directly applied to the DSP problem. This paper develops a modified grey wolf optimizer (MGWO) to solve the DSP problem. The algorithm inherits the significant thought of the social hierarchy and hunting mechanism of the original GWO. The novelty of our GWO is as follows: (1) Three new operators, i.e., feasible solution generator (FSG), neighborhood search operator (NSO), and guided search operator (GSO), are designed for DSP. FSG is developed for generating feasible disassembly sequences. NSO and GSO are designed to update disassembly sequences and make the wolf group evolve in a good direction. (2) In the original GWO, the δ wolf is the third-best wolf. However, through experimental analysis, we found that if δ wolf is set as a new wolf generated by FSG, the performance of GWO is better for solving DSP. (3) In addition, in the original GWO, the first three wolves are preserved to the next generation without being guided by other wolves, but in our GWO, α, β, and δ wolves would exchange their prey information with each other.
The remaining sections of this paper are organized as follows. Section 2 gives the disassembly model for DSP. Section 3 describes the MGWO for obtaining the optimal disassembly sequence. Section 4 verifies the superiority of MGWO by utilizing two case studies and two real-world applications. Section 5 discusses the effectiveness and superiority of MGWO. Finally, S e c t i o n 6 p r o v i d e s t h e c o n c l u s i o n a n d f u t u r e researches.

Geometric constraints of the disassembly process
In the space rectangular coordinate system, the disassembly directions of components are generally divided into 6 directions d k = {−x, −y, −z, +x, +y, +z}. In this paper, a decimal interference matrix IM 10+ is defined to denote the interference relationship between components along +x-, +y-, and +z-axes. The interference relationship between components along −x-, −y-, and -z-axes can be obtained by the inversion of IM 10+ .
where I ij denotes the interference relationship between component i and component j, when component i is disassembled in the positive direction of each axis.     Fig. 5 The pseudo-code of GSO Table 1 An example of executing GSO to obtain a new disassembly sequence Step

Changes of disassembly directions and tools
In disassembling EOL products, we found that changes in disassembly directions and tools increase extra expenses of time and cost and reduce the disassembly efficiency and flexibility. Hence, few changes of disassembly directions and tools benefit disassembly. This paper defines f d (x) to denote the fitness evaluation for changes of disassembly directions and f t (x) to represent the fitness evaluation for changes of disassembly tools.
where X = {x 1 , x 2 , …, x n }denotes the disassembly component sequence, d(x i ) and t(x i ) can be calculated as follows:

Changes of component materials
To facilitate the subsequent recycling of EOL products, we want to disassemble the components with the same materials in priority. Therefore, it is necessary to minimize the material changes of components during a disassembly process.
In this paper, f m (x) represents the fitness evaluation for material changes of the disassembly sequence.
where m(x i ) can be calculated as follows:

Stability of the disassembly process
The stability of disassembly sequences is used to evaluate the stability between connected components during the disassembly process. An unstable situation may cause damage to recyclable components, reduce disassembly efficiency, and increase costs. For this reason, a stability matrix S = (s ij ) n×n is defined to represent the stability between components. If component i supports component j, s ij = 1; otherwise, s ij = 0.
In this paper, f s (x) is defined to represent the fitness evaluation for component stability of the disassembly process.

Fitness function
A disassembly sequence firstly needs to satisfy geometric constraints. Then, according to the literature [18], a linear weighted expression considering the above four indicators is designed to measure the disassembly quality for the feasible sequence.
where w d , w t , w m , and w s are the weight factors of f d (X), f t (X), f m (X), f s (X), respectively. In subsequent experiments, all weight factors are set to 0.25. The optimal value is indicated in bold

Proposed MGWO for DSP
This section introduces the MGWO for the DSP problem in detail. The original GWO is described briefly, firstly. Afterward, the framework of MGWO is proposed. Finally, three critical operators are designed in MGWO for solving the combinatorial optimization problem.

Brief introduction of GWO
GWO is a new meta-heuristic algorithm characterized by presenting the social hierarchy and hunting mechanism of the wolf group. All wolves are divided into four hierarchies based on their fitness values. The best wolf is represented as α, the second wolf is denoted as β, and the third-best wolf is expressed as δ. The remaining wolves are described as ω.
During the search process, α, β, and δ wolves guide ω wolves. The prey is firstly encircled by grey wolves when hunting. Equation 13 is expressed to the model of the encircling behavior.
where D denotes the distance between the prey and a wolf, X and X p indicate the position vector of the wolf and the prey, respectively, and t represents the current iteration. Finally, the vector A and C are formulated using Eqs. (15) and (16).
where r 1 and r 2 are random vectors between 0 and 1, and a decreases linearly from 2 to 0 throughout the iterations. α, β, and δ wolves guide the wolf group to hunt for the prey. Nevertheless, the optimum (prey) position cannot be determined in an unknown and abstract solution space. Aiming at mathematically describing the hunting practices of wolves, the Material  algorithm assumes that α, β, and δ wolves are closer to the prey than ω wolves. Hence, the three best wolves are saved and oblige the other wolves to update their position toward the prey. Equations (17)(18)(19)(20)(21)(22)(23) calculate the hunting practices of the wolf group [29].
In summary, wolves are first divided into four kinds of hierarchies in GWO. The first three wolves lead the remaining wolves to search for the prey. According to Eq. (14), when |A| < 1, wolves will move between themselves and the prey. Otherwise, they attack the prey. Since the actual prey position is unknown, Eqs. (20)(21)(22) assumes the position of the prey is the same as that of α, β, and δ wolves. Finally, the optimal solution is obtained when the stop criterion is satisfied.

Framework of the proposed MGWO
The proposed MGWO inherits the main idea of the social hierarchy and hunting mechanism from the original GWO. Feasible solution generator (FSG), neighborhood search operator (NSO), and guided search operator (GSO) are designed in the algorithm. In the original GWO, the δ wolf is the third-best wolf. However, in the proposed MGWO, the δ wolf is a new wolf generated by FSG. In addition, in the original GWO, the first three wolves are preserved to the next generation without being guided by other wolves, but in the MGWO, α, β, and δ wolves would be guided by each other. The reason for the modification is explained in Section 4.1. The computation results verify that the modified algorithm outperforms the original algorithm in solving the DSP problem. The pseudocode of MGWO is shown in Fig. 1.  The optimal value is indicated in bold In the NS phase, α and β wolves with the best fitness values are identified first in the population, and the δ wolf is generated by FSG. Then, a neighborhood search is performed on each of the three wolves (disassembly sequences). First, a component in the disassembly sequence is randomly selected, and then its earliest insertion point a and latest insertion point b in the original disassembly sequence are determined. Finally, the component is inserted into a random position between [a, b] of the original disassembly sequence, and then the new disassembly sequence is generated.
In the GS phase, α, β, and δ wolves guide each wolf in the population to search the prey. The self-updating factor and the guided factors are calculated according to the fitness values of the chosen wolf and these three wolves. These self-adaptive factors determine the probability that the chosen wolf will be guided by the other three wolves.
A detailed description of FSG, NSO, and GSO will be introduced in the following sections.

Feasible solution generator
The constraint of disassembly precedence should be satisfied in the disassembly process. In our research, the constraint of Tool Fig. 11 Matrices, material, and tool sets of ball valve. a Interference matrix. b Support matrix. c Material set. d Tool set Fig. 12 Box plot of ball valve

Neighborhood search operator
NSO is used for neighborhood search on α, β, and δ wolves, making them update themselves. For each of these three disassembly sequences, a component is randomly selected to change its position or disassembly direction, making the sequence variation. If the fitness value of the disassembly sequence is improved after using NSO, the original sequence will be replaced by the new sequence. In the NS phase, the most critical thing is to find the earliest insertion point a and the latest insertion point b of the selected component. In this paper, the decimal interference matrix IM 10+ is used to find a and b. The pseudo-code of NSO is shown in Fig. 3.

Guided search operator
α, β, and δ wolves execute GWO to guide the wolf group to hunt prey. X represents a guided wolf in the population. The three wolves try to improve the fitness value of X wolf by sharing their positions to lead it to move toward them.
The new solution can be generated according to α, β, and δ wolves, and X wolves using GSO. Firstly, the lowest bound of the fitness value f lb is given according to the best disassembly sequence in the population. Then, four self-adaptive parameters, named self-updating factor p s , α-guiding factor p α , βguiding factor p β , and δ-guiding factor p δ , are used for modifying X solution. Each parameter can be calculated according to the fitness value of X wolf f X , the fitness value of α wolf f α , the fitness value of β wolf f β , the fitness value of δ wolf f δ and the lowest bound of the fitness value f lb . If the fitness value of a solution is lower, its guiding factor is higher. According to  The optimal value is indicated in bold Finally, the precedence preservative crossover (PPX) methodology [18] is used to ensure the precedence relationship between components. The PPX operator passes on the precedence relationship based on two disassembly sequences to a new sequence while ensuring no new precedence relationships are introduced. The pseudo-code of GSO is shown in Fig. 5.
In the GSO procedure, the new disassembly component sequence X new and its corresponding disassembly direction sequence D new are generated by choosing and setting the components one by one from X, α, β, and δ. Firstly, the probability p is randomly generated in the range of [0, 1], compared with p s , p α , p β , p δ . If p is less than p s , the leftmost component in X is pushed back to X new ; if p is less than p s + p α , the leftmost component in α is pushed back to X new ; if p is less than p s + p α + p β , the leftmost component in β is pushed back to X new ; otherwise, the leftmost component in δ is pushed back to X new . Then, the chosen component is removed from the four wolves. Meanwhile, the disassembly direction of the selected component is added to D new . Finally, a new solution is generated by adding components one by one, as above. The precedence relationship between components can be preserved by this method.

Case studies and industrial application
This section aims to demonstrate the effectiveness of MGWO.
Firstly, in Section 4.1, a vise [36] with 11 components and a ball valve [37] with 17 components are used to verify the effectiveness of the modified strategies in MGWO.
Secondly, in Section 4.2, two case studies are implemented to present the superiority of MGWO compared with five algorithms: GA [12], SSO [14], STLBO [25], IHS [36], and ACO [23]. All of these methods have been proven to be effective and efficient in solving the sequence planning problem. To ensure the fairness of comparison, we set the population size of all algorithms to 20 and the maximum number of termination iterations to 1000. Other parameter settings can be referred to in the corresponding literatures.
Finally, a real-world engineering case is used to verify the practicality of the MGWO in Section 4.3.
All algorithms are coded in Matlab 2017a and carried out on an Intel Core 3.3-GHz PC with 4-GB memory. In the subsection, a vise is used to test the effectiveness of the proposed strategies. The assembly view of the vise is shown in Fig. 6, and its corresponding interference and support matrices, material, and tool sets are offered in Fig. 7. MGWO is compared with its variations, including not using NSO (MGWO-1), not using GSO (MGWO-2), α, β, and δ wolves do not participate in the GS phase (MGWO-3), and δ wolf is the third-best wolf, not a newly generated one (MGWO-4).
We set the population size to 20 and the maximum number of termination iterations to 1000. After 30 independent runs, the comparison results of different strategies are shown in Table 2. The average computation time is the average time taken by various strategies to obtain the optimal solution. The box plot of best fitness values obtained by different strategies in 30 runs is shown in Fig. 8. It can be observed that MGWO, MGWO-1, and MGWO-3 can obtain the optimal solution, but the mean of fitness values (3.8667) of MGWO is better (lower) than that of other strategies. Figure 8 presents the box plot of fitness values of the vise in 30 independent runs, showing that the solutions obtained by MGWO are more concentrated than its variants. MGWO can obtain the optimal value (3.7500) with the probability of 63.33%. MGWO-1 and MGWO-3 can obtain the optimal value with the probability of 30.00%. Additionally, MGWO-2 and MGWO-4 cannot obtain the optimal value. Compared with MGWO-1, MGWO makes α, β, and δ wolves have more opportunities to get the optimal solution through self-updating using NSO. Compared with MGWO-2, MGWO allows each wolf to be guided by the best wolves instead of inefficiently searching for the prey by itself. Compared with MGWO-3, in MGWO, α, β, and δ wolves would exchange their prey information, enabling them to move together toward the prey. Compared with MGWO-4, δ wolf is a new one in MGWO, increasing the diversity of the population and avoids the algorithm from falling into the local optimum. Figure 9 presents the convergence curves of MGWO compared with four other strategies with the optimal solution in 30 independent runs. It reveals that MGWO has the highest probability of obtaining the optimal solution and converges faster to the optimal solution than other strategies. In Table 3, the Wilcoxon rank-sum test with a significant level of 0.05 is used to measure the significance between MGWO with other strategies. Significance value = 1 indicates the result obtained by MGWO is significantly different from the other strategies used in comparison; significance value = 0 represents the result obtained by MGWO and is not substantially different from the other strategies used in the comparison. Obviously, our proposed modified strategies are The optimal value is indicated in bold  Fig. 16 Box plot of ball valve effective, and the performance of MGWO is significantly better than its variants.

Case 2: ball valve
In this case, to further measure the significance of our proposed strategies, we compare MGWO with MGWO-1, MGWO-2, MGWO-3, and MGWO-4 for another case, the ball valve [37] with 17 components. The assembly view of the ball valve is shown in Fig. 10, and its corresponding interference and support matrices, material, and tool sets are offered in Fig. 11. We set the population size to 20 and the maximum number of termination iterations to 1000. After 30 independent runs, the comparison results of different strategies are shown in Table 4. The average computation time is the average time taken by various strategies to obtain the optimal solution. The box plot of best fitness values obtained by different strategies in 30 runs is shown in Fig. 12. Table 4 shows that all strategies can obtain the optimal solution, but the mean of fitness values (3.7500) of MGWO is better than that of other strategies. Figure 12 presents the box plot of fitness values of the ball valve in 30 independent runs, showing that the solutions obtained by MGWO are more concentrated than those of other algorithms. MGWO, MGWO-1, MGWO-2, MGWO-3, and MGWO-4 can reach the optimal value (3.7500) with the probability of 100%, 70%, 13.33%, 43.33%, and 20%. Figure 13 shows the convergence curves of MGWO compared with other strategies with the optimal solution in 30 independent runs. It reveals that compared with MGWO-1, MGWO-2, MGWO-3, and MGWO-4, MGWO has a higher probability and faster convergence to obtain the optimal solution. Table 5 uses the Wilcoxon ranksum test with a significant level of 0.05 to measure the significance between MGWO with other strategies. It shows that MGWO outperforms other strategies in solving the disassembly sequence of the ball valve.

Case 1: vise
In this subsection, the vise introduced in Section 4.1.1 is used to test the performance of MGWO. Five classic algorithms, GA, ACO, SSO, STLBO, and IHS, are used to compare with MGWO. Like Section 4.1, we set the population size to 20 and the maximum number of termination iterations to 1000. After 30 independent runs, Table 6 presents the comparison results of the six algorithms. The average computation time is the average time taken by algorithms to get the optimal disassembly sequence. The box plot of the distribution of best fitness in 30 runs is shown in Fig. 14.
We can observe from Table 6 that only MGWO, GA, and ACO can obtain the optimal solution, and the mean of fitness values (3.8667) of MGWO is better than that of other algorithms. Figure 14 also presents that the solutions obtained by MGWO are better than other algorithms. MGWO, GA, and ACO can receive the optimal solution (3.7500) with the probability of 63.33%, 13.33%, and 16.67%. Additionally, other algorithms cannot get the optimal solution. The optimal disassembly component sequence is (7,6,9,1,11,10,8,3,4,2,5), and its corresponding disassembly direction sequence is (−y, −y, +z, +y, +z, +z, +z, +z, +z, +z, +z). Figure 15 presents the convergence curves of all algorithms with the optimal solution in 30 independent runs. It reveals that MGWO has the highest probability of obtaining the optimal solution and converges to the optimal solution with the least number of iterations in all algorithms. In Table 7, the Wilcoxon ranksum test with the significant level of 0.05 is used to test the

Case 2: ball valve
To further measure the effectiveness of MGWO, we compare the algorithm with GA, ACO, SSO, STLBO, and IHS for the ball valve introduced in Section 4.1.2. We set the population size to 20 and the maximum number of termination iterations to 1000. After 30 independent runs, the experimental results of all algorithms are presented in Table 8. The average computation time is the average time taken by algorithms to obtain the optimal solution. Table 8 shows that all algorithms can obtain the optimal solution except SSO, and the mean of fitness values (3.7500) of MGWO is better than that of other algorithms. Figure 16 presents the box plot of fitness values of the ball valve in 30 independent runs, showing that the solutions obtained by MGWO are more concentrated than other algorithms. MGWO, GA, ACO, IHS, and STLBO can achieve the optimal value (3.7500) with the probability of 100%, 16.67%, 20.00%, 3.33%, and 13.33%, respectively. In addition, SSO cannot reach optimal value. The optimal disassembly component sequence is (16,17,11,10,8,9,12,6,13,7,5,3,4,15,14,2,1), and its corresponding disassembly direction sequence is (+z, +z, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y, −y). Figure 17 shows the convergence curves of all algorithms with the optimal solution in 30 independent runs. It reveals that MGWO has the highest probability and the fastest convergence speed to obtain the optimal solution in all algorithms. Table 9 uses the Wilcoxon rank-sum test with a Disassemblealong -Z and then along +Y Disassemble along +Y

Application 1: azimuth thruster propeller
Used to verify the practicality of MGWO, azimuth thruster propeller [36], a vital product of a Chinese boat company, is employed in this section. Azimuth thruster is the critical equipment of the marine dynamic positioning system. Propeller is used to generate thrust for the azimuth thruster. Many components of the azimuth thruster propeller can be recycled and reused at the end of life. A good disassembly sequence can dramatically improve the recovery efficiency of these components, so the DSP of the azimuth thruster propeller has significant research value. The azimuth thruster propeller has been modified for simplification to satisfy the requirement of disassembly modeling. The exploded view of the azimuth thruster propeller is shown in Fig. 18. The simplification of the disassembling crank is shown in Fig. 19, and the matrices of the azimuth thruster propeller are presented in Fig.  20.
The proposed MGWO is applied to the actual case, and GA, ACO, SSO, STLBO, and IHS are used to compare with MGWO. We set the population size to 20 and the maximum number of termination iterations to 1000. After 30 independent runs, the statistical results of all algorithms are summarized in Table 10, and their corresponding box plots are presented in Fig. 21. The mean of fitness values (5.5669) of MGWO is better (lower) than that of other algorithms. MGWO, GA, ACO, IHS, and STLBO can obtain the optimal solution with the probability of 76.67%, 20.00%, 50.00%, 30.00%, and 73.33%, respectively. In addition, SSO cannot obtain the optimal solution. It can be observed that MGWO performs better than the compared algorithms in terms of the probability of obtaining the optimal solution and the concentration of solutions in this case. The results show that MGWO is very suitable to solve the DSP problem. Figure 22 shows the convergence curves of all algorithms with the optimal solution in 30 independent runs. It demonstrates that MGWO has the highest probability and the fastest convergence speed to obtain the optimal solution in all algorithms. In Table 11, the Wilcoxon rank-sum test with a significant level of 0.05 is used to measure the significance between MGWO with other algorithms. It shows that MGWO performs significantly better than other algorithms except for STLBO. MGWO can reach the optimal value of 5.5000 with the corresponding disassembly component sequence (1,2,3,8,17,14,11,5,12,9,6,15,13,16,7,10,4) and the disassembly direction sequence (−y, −y, −y, −x, −z, +x, +z, +y, +y, +y, +y, +y, +y, +y, +y, +y, +y), which is  accorded with the actual industrial disassembly sequence. Thus, the real-world application demonstrates the validity and practicality of the proposed method.

Application 2: robot arm
In this section, we use a robot arm [38] to validate the practicability of MGWO further. The robot arm consists of 30 functional components. It is more complex than the azimuth thruster propeller. The assembly of the robot arm is shown in Fig. 23, and its corresponding interference and support matrices, material, and tool sets are provided in Fig. 24. We compare MGWO with GA, ACO, SSO, STLBO, and HIS on the robot arm. The population size is set to 20, and the maximum number of termination iterations is set to 1000. Table 12 gives the statistical results of all algorithms, and Fig. 25 provides their corresponding box plots. The mean of fitness values (9.6250) of MGWO is better than that of ACO, IHS, STLBO, and SSO and is slightly worse than that of GA. The optimal value (8.7500) can only be obtained by MGWO, and cannot be found by the other algorithms. Besides, the computation time of MGWO is shorter than that of all other algorithms. The experimental results demonstrate that MGWO outperforms the other algorithms in both solution quality and computational efficiency. Figure 26 shows the convergence curves of all algorithms with the optimal solution in 30 independent runs. It presents that only MGWO can converge to the optimal solution. In Table 13, the Wilcoxon rank-sum test with a significant level of 0.05 is used to measure the significance between MGWO with other algorithms. It indicates that MGWO significantly outperforms other algorithms except for GA. MGWO can obtain the optimal disassembly component sequence (26,30,29,28,14,3,27,5,25,19,13,12,11,9,7,4,6,20,18,10,17,16,22,21,15,2,23,1,8,24) and its corresponding disassembly direction sequence (−z, −z, −z, −z, +y, +x, +x, +x, +x, +z, +z, −x, −x, −x, −x, −x, −x, +z, +z, +z, +z, +z ,+y ,+y, +y, +y, +y, +y, +y, +z). The two sequences are accorded with the actual industrial disassembly sequence. Hence, the industrial application of the robot arm verifies the practicality of our proposed method.

Discussion
The above experimental results demonstrate that MGWO outperforms the other compared algorithms. This is detailed below. In general, MGWO outperforms the other algorithms in both convergence speed and solution quality.
In Section 4.1, Tables 2 and 3 show that MGWO obtains the optimal solution with the highest probability in all strategies. For the vise, MGWO receives the optimal solution with the probability of 63.33%, but the probability of other strategies obtaining the optimal solution is less than 30%. MGWO obtains the optimal solution with the probability of 100% for the ball valve, but the probability of other strategies obtaining the optimal solution is less than 70%. Figures 9 and 13 illustrate that the convergence speed of MGWO is faster than that of other strategies. These experimental results suggest that the effectiveness of NSO and GSO operators.
In Section 4.2, Tables 6 and 8 show that MGWO obtains the optimal solution with the highest probability in all algorithms. For the vise, MGWO obtains the optimal solution with the probability of 63.33%, but the probability of other algorithms obtaining the optimal solution is less than 16.67%. MGWO obtains the optimal solution with the probability of 100% for the ball     valve, but the probability of other strategies obtaining the optimal solution is less than 20.00%. Figures 15 and  17 illustrate that MGWO converges faster than other algorithms. The computation results verify that MGWO is more appropriate for solving the DSP problem than other compared algorithms. In Section 4.3, two real-world industrial applications prove that MGWO can obtain the disassembly sequence following practical production. MGWO receives the optimal solution with the probability of 76.67% for the azimuth thruster propeller, but the probability of other algorithms obtaining the optimal solution is less than 73.33%. Only MGWO can obtain the optimal solution for the robot arm, and other comparison algorithms cannot find it. Figure 22 illustrates that MGWO converges faster than other algorithms. Figure 26 shows that only MGWO can converge to the optimal solution.
The outstanding performance of MGWO for the DSP problem mostly depends on NSO and GSO operators. The two operators can well balance local search and global search. Additionally, the hierarchy and hunting mechanism allows the population to evolve in a good direction without falling into the local optimum. Another important outcome is that MGWO can be well applied to engineering cases, proving our proposed algorithm can effectively guide the disassembly practice.

Conclusions and future researches
In this paper, an MGWO algorithm is proposed for solving the DSP problem. Three key operators, FSG, NSO, and GSO, are designed for the problem. First, two engineering cases are employed to test the effectiveness of these operators. Then, they and two real-world applications are used to test the validity of MGWO. The experimental results verify that MGWO can solve the DSP problem effectively. The advantages of MGWO are as follows: & Three key operators are designed in this paper: FSG, NSO, and GSO. FSG generates the feasible disassembly sequence; NSO and GSO balance both the local search and global search. & The implementation of MGWO is simple. Moreover, all parameters of MGWO are self-adapted, so its performance is stable. & Due to the hierarchy and hunting mechanism of the wolf group, the convergence speed of MGWO is very fast.
Although MGWO shows a good performance for solving the DSP problem, the algorithm still has some limitations. One is that only three leading wolves execute the neighborhood search operator, and the remaining wolves obey the three wolves and lack self-exploration. Another limitation is that our mathematical model considers only one constraint of component interferences, but other constraints also need to be introduced into the model. Concerning future work, the first research direction is to design a novel self-search mechanism for remaining wolves to enhance their self-explore ability. The  Data availability All data generated or analyzed during this study are included in this paper.

Declarations
Ethical approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Competing interests The authors declare no competing interests.