Balancing Stochastic Mixed-Model Two-Sided Disassembly Line Using Multiobjective Genetic Flatworm Algorithm

The disassembly industry is still labor-intensive, and for the disassembly of products with similar assembly structures, achieving mixed disassembly of these products on the same disassembly line is a more economical solution and very popular in disassembly companies. When large-volume end-of-life (EOL) products are disassembled on the two-sided disassembly line, there is uncertainty in the disassembly time during the disassembly process due to the differences in the recovery quality of EOL products. Therefore, a stochastic mixed-integer programming model is formulated that considers workstation activation cost, workload smoothness index between workstations, and total hazard index considering the order of hazardous task removal. Then, a multi-objective genetic flatworm algorithm is developed in which two different mechanisms for new individual generation are embedded, namely the flatworm genetic operations and the regeneration operations. The effectiveness of the model and the efficiency of the algorithm are verified by the disassembly cases of the printer and the classical mixed-model on the two-sided disassembly line. Finally, the proposed model and algorithm are applied to the disassembly of multiple types of cars, which proves their practical industrial application value.


I. INTRODUCTION
With the rapid development of science and technology, the lifecycle of products is constantly shortened. Many products are rapidly eliminated due to outdated functions or scrapped, more and more waste is generated [1]. This is particularly prominent for China, as the manufacturing center of the whole world. The government proposed to comprehensively promote cleaner production, develop the circular economy, improve the efficiency of resource recycling, and build a green manufacturing system in ''made in China 2025'' [2]. To solve these problems, improve resource utilization and protect the environment, it is necessary to vigorously The associate editor coordinating the review of this manuscript and approving it for publication was Seyedali Mirjalili . implement sustainable waste utilization and manufacturing technologies. Remanufacturing is the manufacturing of endof-life (EOL) products as raw materials, not only can effectively save energy and resources, but also can greatly reduce production costs. Disassembly should be carried out precede remanufacturing by removing valuable parts from the EOL products. In recent years, the paced disassembly line has had a great efficiency advantage in handling a large number of EOL products to be disassembled.
When disassembly is performed, the parts that make up the EOL products not only have disassembly priorities among themselves, but some may also carry toxic residues. How to assign the disassembly tasks to the stations on the disassembly line in the proper order to achieve the optimal disassembly goals under these constraints is called the disassembly line VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ balancing problem (DLBP) [3], which was first proposed by Gupta et al. and has attracted a lot of attention from the disassembly industry. The layout of the disassembly line is mainly in the form of straight [4], U-shaped [5], parallel [6]. The workstations on these types of disassembly lines are arranged on one side and are particularly suitable for the disassembly of waste washing machines, microwave ovens, and other white and black household appliances, whose volumes are relatively small. For large-volume refrigerators and cars, the abovementioned disassembly lines are not as efficient and costeffective, which makes the two-sided disassembly line an efficient alternative. Kucukkoc [7] described the two-sided disassembly line problem and modeled it mathematically for the first time. Wang et al. solved the partial disassembly problem of EOL printers using the two-sided disassembly line [8]. Liang et al. [2], [9] studied the complex execution constraints and fuzzy operation time problems in the two-sided disassembly lines, respectively. All of these cases directly illustrate the efficiency and cost advantages of the two-sided disassembly line for dealing with large-volume EOL products.
DLBP is mainly divided into two main types according to the different research objectives. The first type (DLBP-I) is to determine the minimum number of workstations for a given cycle time, which is usually carried out during the design and installation stage for the disassembly industry. The second type (DLBP-II) is to minimize the cycle time with a fixed number of stations, which is usually optimized for the existing line to improve the disassembly flexibility [10]. Almost all research on DLP is devoted to Type-I and extended the optimization goals to other aspects based on this to make it more realistic for disassembly companies. Ren et al. [11] studied the total disassembly revenue and cost. Wang et al. [12], [13] studied the disassembly energy consumption, workload smoothness among the stations, and disassembly profit. EDIS [14] additionally optimized the goals of disassembly direction changes, total demand index, and hazard index. More details of the research objectives can be found in the recent literature review [15].
Different from assembly line balancing issues, parts of EOL products with disassembly may have physical defects or functional defects, or both, causing many uncertainties in the disassembly process, the most intuitive manifestation of which is that the disassembly time is not a constant value, but a certain amount of fluctuation. Alavidoost et al. [16], Zhang et al. [17], and Kalayci et al. [18] expressed the task disassembly time as a triangular fuzzy number, and [2] introduced the α-cut sets to enrich the diversity of fuzzy disassembly model cases. In recent studies, [6], [19]- [21] simulated the uncertainty disassembly time as a number that conforms to a normal random distribution so that the workload of each workstation does not exceed the cycle time constraint with different confidence levels.
DLBP was proved to be an NP-complete problem [22]. There are two main ways to solve it since the problem was proposed. At the beginning of the research, due to the small size of the DLBP disassembly task, the desired results can be obtained in a short time by using mixed-integer linear programming and branch-and-bound methods [23]. Li et al. [24] proposed a branch, bound, and remember method to solve the simple DLBP recently. With the increase of the scale of DLBP, the exact methods cannot give acceptable disassembly solutions in a reasonable time, and heuristic methods and metaheuristic algorithms were proposed [25]. Typical methods include genetic algorithm [7], ant colony optimization [26], particle swarm optimization [27], [28], simulated annealing algorithm [29], [30], artificial bee colony [31], [32] and variable neighborhood search [1]. Recently, some new algorithms, such as memetic algorithm [33], whale optimization algorithm [34], flower pollination algorithm [35], migrating birds optimization algorithm [36], non-dominated sorting genetic algorithm [37], etc., have also efficiently solved the DLBP in large-scale.
Some scholars have also studied the robot disassembly line [25], [30], but in practice, the disassembly industry is still labor-intensive, and disassembly shops usually have only a limited number of paced disassembly lines. For disassembly of EOL products that are similar in terms of assembly structure, a more economical solution is to achieve mixed disassembly of these products on the same disassembly line. Especially for large-volume of EOL products with different kinds of individual configurations, when disassembled on two-sided disassembly lines, the disassembly tasks may be different due to differences in part configurations, or the disassembly of the same part may lead to different operating times due to the different quality of EOL products. For example, various brands of decommissioned vehicles, some have sunroofs and some do not. Some hood bolts may be corroded and will result in longer disassembly times. However, to the best of the author's knowledge, only Mutlu et al. [33] have studied the balancing of mixed-model two-sided disassembly line problem (M-TDLBP), whose model not only ignores the uncertainty of the part removal time but also optimizes only the total number of opened stations.
In view of the fact that the mixed-mode two-sided line is very popular in real disassembly companies because of the great flexibility it exhibits in the face of disassembly of EOL products with similar structures, and that there is less related research work, this paper further complements and improves the research in this area. The main contributions are as follows.
1) A mixed-integer programming (MIP) model is formulated for the stochastic mixed-mode two-sided disassembly line problem (SM-TDLBP). The model considers both the cost of opened stations, the workload smoothness index among workstations, and the total hazard index considering the removal sequence of hazard tasks. 2) A multiobjective genetic flatworm algorithm (GFA) for solving SM-TDLBP is developed, in which two different mechanisms for generating new individuals, the flatworm genetic operation and the regeneration operation, are embedded to improve the efficiency of solving SM-TDLBP. 3) Providing rich case information for the SM-TDLBP study. Includes not only modifications of cases from existing literature, but also data and information on disassembly of actual mixed EOL vehicles. The rest of this paper is structured as follows. Section II describes the problem and formulates the MIP model of SM-TDLBP. Section III explains the proposed multiobjective genetic flatworm algorithm, and then Section IV gives computational examples to validate the performance of the proposed model and algorithm. The last section discusses the conclusions, including directions for further research. The two workstations directly facing each other are called a mated-station [7], and each workstation is equipped with the appropriate operator and disassembly tools as needed. The workstation on the left side can only disassemble those parts of the EOL product that are to the left, as is the workstation on the right. If multiple incoming EOL products with similar structure are disassembled in a mixed flow on the running paced two-sided disassembly line, then is called the mixedmodel two-sided disassembly line.
Each EOL model has its task disassembly precedence relationship (DPR), and all DPRs can be combined into a joint disassembly precedence diagram. It should be pointed out that there is no DPR conflict for the same disassembly task contained in different models.
The disassembly quantity for each EOL product during the scheduled disassembly planning period is generally different, and to simplify the problem formulation, it is assumed that the disassembly quantity of each EOL model is equal.

B. STOCHASTIC DISASSEMBLY TIME
To facilitate the presentation of SM-TDLBP, the notations used (including indices, parameters, decision variables, and indicative variables) are summarized in the Appendix.
Assuming that all the task disassembly time in each model follows an independent normal distribution, and the mean value and variance are known, denoted as t im ∼ N (µ im , σ 2 im ). To handle stochastic task times, Özcan et al. [38] proposed a simplified treatment that uses the additive property of independent normal random distributions to determine the magnitude of each other. However, the maximum of two normal variables is not normally distributed, Chiang et al. [39] and Tang et al. [40] improved this model by considering the joint probability that the workstations on both sides of the same mated-station are completed within the cycle time. The difficulty of the stochastic mix-model on the two-sided disassembly line lies in the handling of unfinished tasks.
If unfinished tasks can be completed without affecting the subsequent mated-stations, then the completion probability can be represented by formula (1), and all mated-stations can be treated individually [39], [40].
If the disassembly line must be stopped in order to wait for the unfinished task to complete, then all subsequent matedstations will be affected, and formula (2) is adopted to handle the completion probability under this condition.
Let be the probability density function of the standard normal distribution. Clearly, the joint probability of any two tasks being assigned to each side of the same mated-station to be completed without exceeding the cycle time constraint is greater than the probability of them being assigned to the same station. 

C. OBJECTIVES AND CONSTRAINTS
A mathematical model for the stochastic two-sided partial disassembly line balancing problem was given by Chiang et al. [39] and Tang et al. [40]. Based on this, the MIP model of the stochastic mixed-model two-sided disassembly line balancing problem can be constructed as shown below.
Three objectives need to be optimized in SM-TDLBP. The workstation activation cost is expressed by formula (4), which is used to calculate the number of opened workstations to meet the demand of all disassembly EOL products. β is the sharing coefficient of disassembly equipment and tools in a mated-station, which is set to 0.8. For the workstation with only one side is utilized, the coefficient is 1. Therefore, the lower the starting cost of the workstation, the better. The weighted smoothness index is designed to measure the load deviation between mated-stations and it is calculated by the formula (5). The smaller the smoothness index means the smaller the load deviation between different mated-stations, and 0 is the ideal state. The total hazard index of EOL products is calculated by the formula (6), which is related not VOLUME 9, 2021 only to the number of each EOL product to disassembly, and the hazard degree of each task, but also related to the disassembly sequence of the hazardous tasks. (t s im + t f im )/2CT is the penalty coefficient indicating that tasks with a higher hazard level will also have a higher hazard index if the later they are removed. The lower hazard index is more preferred.
The following constraints need to be met when assigning disassembly tasks to workstations.
Among the constraints, formula (7) ensures that each disassembly task can only be assigned to exactly one station. Formula (8) indicates that the DPR constraint is satisfied. Formula (9) controls the disassembly completion time of any task on the station that does not exceed the cycle time. Formula (10) ensures that any two disassembly tasks with DPR must be satisfied when they are assigned to the same workstation or two different workstations in a mated-station. The formula adopts the fact that there is no overlap between the start disassembly time and the completion disassembly time between the two tasks. For tasks g and i, if there is no DPR between them, they can theoretically be disassembled in any order. When they are assigned to the same workstation, the order for disassembly between them should be specified. Formula (11) indicates that task i is operated before task g, and formula (12) indicates that task i is disassembled after task g. When they are assigned to different workstations, a positive infinity integer is used to relieve the constraints between each other. Formula (13) indicates that the mean disassembly time for any task is a nonnegative value, and 0 means that the disassembly task does not exist. Formula (14) indicates that the variance of disassembly time for any task is a nonnegative value, and 0 means that the disassembly time is deterministic.

D. AN EXAMPLE
A stochastic mixed-model two-sided disassembly example is provided here, which is modified from Ozcan et al. [41].  it can be seen that a total of 16 tasks need to be completely disassembled on the joint disassembly precedence diagram, and Table 1 illustrates the disassembly order between tasks, disassembly direction, hazard degree, and disassembly time. The hazard degree is specified between 0 and 0.5 to simulate different impacts on the environment. The mean disassembly time is directly obtained from Ozcan et al. [41], and the task variance is randomly generated between 0 and 1/16 times the mean. There are two types of models A and B that need to be disassembled, and the tasks that need to be disassembled are not exactly the same in models A and B. For example, Model A does not have tasks 12, 9, and 4, while Model B does not have tasks 5, 3, and 1. Even for the same tasks of different models, their disassembly times are not exactly the same. For example, the disassembly time for task 6 on model A is (4,0.14), while it is (8,0.76) on model B. In the disassembly direction constraint, L means that the task can only be operated on the left station, R means the right, and E denotes that can be disassembled in the workstation on either side of the disassembly line.
Under the premise of satisfying the constraints of DPR, disassembly direction, and cycle time, all disassembly tasks of models A and B are assigned to the workstations on the line, and One of the disassembly schemes is shown in Fig. 1. The cycle time of this example is set to 18 time units, and the confidence level is 0.9.
In the disassembly scheme in Fig. 1, three mated-stations are utilized to complete all the disassembly tasks of models A and B. It is worth noting that there are two types of idle time in this disassembly scheme. The first type is the time left after the task is completed on the workstation, and the other is caused by the DPR constraint between tasks, which is called sequence-dependent idle time. For example, for model B, task 11 is the immediate successor of tasks 14 and 15, yet task 11 follows immediately after task 15 because task 15 takes longer to disassemble than task 14.

III. GENETIC FLATWORM ALGORITHM
The flatworm algorithm (FA) was first proposed and applied to the disassembly sequence planning problem by Tseng [42]. The algorithm mainly simulates the process of individual growth, splitting and regeneration to produce new individuals when the flatworms in nature encounter stimuli and retains the best individuals through a penalty mechanism to promote the population's continuous optimization until convergence. The process of generating new individuals through partial tissue regeneration is one of asexual reproduction. However, in addition to producing new individuals by tissue reproduction, most of them can also reproduce sexually through self-fertilization or hetero-fertilization to generate individuals with a richer gene pool. Among evolutionary intelligence algorithms, the genetic algorithm (GA) can particularly reflect the great hybridization advantages of sexual proliferation through the crossover and mutation operations between gene fragments. Therefore, a hybrid model of GA and FA is introduced in this paper, called the genetic flatworm algorithm (GFA), whose computational flow chart is detailed in Fig. 2, and each part of the proposed algorithm is illustrated in the next.
There are two main categories of parameter initialization: one is parameters related to each disassembly product itself, such as the fields in Table 1. The second is the parameters required by the algorithm, such as the total of iterations and population size, which are introduced in Section IV.

A. SOLUTION ENCODING AND DECODING
Efficient and reasonable individual encoding is not only a prerequisite for population construction but also the key to solving the SM-TDLBP, which needs to satisfy the following two requirements: firstly, it must ensure that the individual encoding is a feasible disassembly sequence, that is, it does not violate any DPR constraints. Secondly, it can perform the genetic operation and flatworm regeneration process stably, accurately and efficiently.
Task sequence-oriented encoding techniques are employed, where each gene of an individual represents a disassembly task with a length equal to the total number of tasks in the joint disassembly relationship graph. A feasible encoding scheme is depicted in Fig. 3. The number in the circle represents the disassembly task, the red color indicates that the disassembly task is hazardous to the environment, and the letter above the task indicates its initial preferred disassembly direction. It should be noted that the disassembly direction for type E in the coding scheme is not specified to be left or right, but is determined dynamically during decoding.
When decoding, Clark's approximation method [39], [40], [43] is used to approximate the probability of completion time. The detailed decoding procedures are as follows: Step 1: Start decoding and open the first mated-station.
Step 2: For each model, judge whether all tasks in the encoding individual have been assigned to the mated-station. If so, skip to step 6. Otherwise, go to step 3.
Step 3: Take the first task in the encoding individual as the active task to be assigned. If the cycle time is satisfied, go to step 4. Otherwise, skip to step 5.
Step 4: If the current task has a clear disassembly direction (L or R), assign it to the left or right side of the mated-station. If E, the current task can be randomly assigned to the left or right of the mated-station. After the assignment, the current task is removed from the encoding individual.
Step 5: Close the current mated-station and open a new one. Then skip to step 2 to continue.
Step 6: The decoding process is finished and the decoding result is output.
After decoding, the final actual disassembly direction of the individual is determined as shown below the task in Fig. 3, and the corresponding decoding result is shown in Fig. 1.
Next, a population containing a specified number of coded individuals can be created, and this population is then divided into two subpopulations, with subpopulation I performing the  flatworm genetic operations and subpopulation II performing the flatworm regeneration processes.

B. THE GENETIC OPERATIONS
The genetic operations of flatworms in the subpopulation I mainly consists of crossover and mutation. All individuals in the subpopulation are involved in the crossover operation to ensure the diversity of feasible solutions. At the same time, in order to prevent falling into random search and to retain the good individuals, the individuals after crossover participate in mutation with adaptive mutation probability.

1) THE CROSSOVER
A technique of mapping information fragments between two gene loci is applied to the crossover operation. An example is illustrated in Fig. 4, and the specific procedures are described below.
Step 1: Two different individuals are randomly selected as parents in subpopulation I.
Step 2: Then two genetic loci are randomly selected as the crossover points. Gene locus 7 and 12 are selected in Fig. 4.
Step 3: Taking P1 as the parental individual and P2 as the reference individual, the gene fragments on the left and right sides of the crossover points are kept unchanged. Intercepting the gene fragment between the crossover points of P1, and {10, 8, 6, 3, 9, 7} is obtained, then finding the corresponding gene fragment of this fragment in P2 separately, ensuring that its position index is sorted from smallest to largest, and then {8, 6, 9, 10, 7, 3} is obtained. Swapping these two gene fragments to get the offspring 1, denoted as O1.
Step 4: Similarly, setting P2 as the parental individual and P1 as the reference individual, repeat step 3 will get O2.

2) THE MUTATION
The mutation probability P m is set to P m max − l/l max · (P m max − P m min ) to ensure that the mutation probability is larger at the start to get as many feasible solutions, while  the mutation probability gradually decreases as the iteration proceeds to ensure the stability of the solution and gradually converge to the optimal.
The mutation case of O1 obtained after the crossover is depicted in Fig. 5. First, a mutation task is randomly selected, assumed to be task 3. Then according to the DPR described in Table 1, the immediate predecessor and successor tasks are determined, which are task 6 and task 1 respectively. Finally, the position where task 3 can be inserted is after task 6 or 9 or 10, and before task 1. The mutation individual (M1) is obtained when is placed after task 9.

C. THE FLATWORM REGENERATION PROCESSES
The process of asexual proliferation of individual flatworms in nature is graphically depicted in Fig. 6. After the growth, a flatworm individual is cut and split into three parts: head, body, and tail. Then, after the regeneration process of tissue stem cells, the three flatworm tissues will grow into three complete flatworm individuals by their regenerative ability to grow the missing parts [42]. The flatworm individuals in the subpopulation II of the proposed algorithm can also simulate the above three processes.

1) THE GROWTH PROCESS
On the one hand, to simulate the principle of tissue stem cell proliferation with real flatworms, this paper assumes that the disassembly task on each gene locus is the basic unit of  growth. On the other hand, a fixed growth probability is set to 0.1 to avoid the infinite growth of the encoded individual.
An example of individual growth is depicted in Fig. 7. For an individual with 16 disassembly tasks, only two tasks are self-replicated separately, assuming that tasks 10 and 9 are selected randomly and growth occurs, then they are assigned to the reasonable locations of the original encoded individual according to the DPR constraint. Task 10 can be inserted after task 13 and before 7, while task 9 can only be inserted after 12 (the immediate predecessor tasks of 9 are 13 and 12, and task 12 is relatively backward in the encoding individual) and before task 7. The reproduced tasks maintain their original disassembly direction, and the grown individual becomes an infeasible disassembly sequence due to the self-replication tasks, but it makes organizational preparation for the next stage of splitting.

2) THE SPLITTING PROCESS
According to related theoretical studies, the grown flatworm can be cut to 1/279 the size of the body and still regenerate a complete flatworm [44]. If the encoded individual is cut too tiny in the splitting process, it will not be conducive to the preservation of excellent individuals in the regeneration process. Therefore, the probability of splitting is assumed to be equal to the probability of growth, and the splitting points are generated randomly. Fig. 8 depicts the splitting process. Two split points are randomly generated after tasks 10 and 9, and the complete individual is cut into three parts: head, body, and tail.

3) THE REGENERATION PROCESS
The regeneration process mainly includes the following steps.
1) Make sure that the disassembly task is unique in each partial individual after the split. If a task is not unique, redundant tasks are randomly deleted. 2) For each segment of an individual, all the missing tasks are remembered according to the disassembly information. Then, these lost disassembly tasks are regenerated into one or more feasible disassembly sequences according to the DPR constraint. 3) Finally, the original fragments and the regenerated fragments are orderly spliced together to reconstruct a complete flatworm individual. Taking the regeneration of the tail fragment in Fig. 9 as an example. The missing disassembly task set is {6, 8,10,11,12, 13,14,15,16}, and during regeneration, the tasks with hazards can be prioritized for removal to accelerate convergence, the tasks in the set are sequenced as {15, 14,11,16,12, 13,8,6,10}. It is important to highlight that the regeneration process is considered to be instantaneous, regardless of the length of the original individual segment. In other words, there is no time difference in the regeneration process of three parts.

D. MULTIOBJECTIVE ELITE RETENTION METHOD
This paper minimizes the workstation activation cost, the disassembly line smoothness index and the total hazard index simultaneously, which has been proved in literature [2], [9] that there is a conflict among the objectives, that is, if one of the objectives is improved, it will cause the deterioration of the other objectives. Hence, SM-TDLBP is a multiobjective optimization problem.
To compare the quality of solutions, the concept of Paretodomination is introduced. Taking the minimization problem as an example, suppose there are two feasible solutions X 1 and X 2 , and f (X ) as the corresponding objective value. If f (X 1 ) ≤ f (X 2 ) is true for all objectives, and f (X 1 ) < f (X 2 ) is true in at least one objective, then X 2 is said to be dominated by X 1 or X 1 performs better than X 2 . If a solution set is not dominated by any other feasible solutions, then the solution set is called the Pareto optimal solution [45]. The Pareto elite selection and retention strategy is to store a certain number of excellent solutions in new populations.

IV. NUMERICAL EXPERIMENTS VERIFICATION
In this section, the efficiency of the proposed multiobjective GFA is demonstrated by different cases of mixed-model twosided disassembly line. All algorithms were coded in Matlab R2016a on Windows 7 with 3.6GHz CPU and 4GB RAM.

A. EXPERIMENTAL EVALUATION INDEX
In general, both heuristic and meta-heuristic optimization algorithms can only obtain a group of better approximate Pareto optimal solutions, and the quality evaluation of the solution set should take into account the convergence degree VOLUME 9, 2021 between the obtained Pareto optimal solutions with the true Pareto optimal solutions, the uniformity and extensiveness of the distribution in the target space, etc. The hypervolume (HV) metric is favored to evaluate the quality of the obtained Pareto optimal solutions, which was first proposed by Zitzler et al. [46], and HV denotes the volume of the hypercube enclosed by the individuals in the Pareto optimal solution set and the reference point in the target space [47]. According to the meaning of HV, if the Pareto optimal solution set PS 1 obtained by the algorithm of A 1 is better than PS 2 obtained by A 2 , then the HV index of PS 1 will be greater than that of PS 2 , and it can be further inferred that A 1 is more efficient than A 2 . Based on the advantages of HV, it is also embedded in the proposed GFA as an elite individual retention strategy with the Pareto domination together.

B. CASE ONE
Case one is taken from Wang et al. [8], who studied the disassembly of the single-model of a printer with 55 tasks on the two-sided disassembly line. The research assumes that the disassembly time of each task obeys a random normal distribution, and uses an improved multi-objective discrete flower pollination algorithm (IMOFPA) to obtain a set of Pareto optimal solutions. Obviously, this research is a special case of mixed-model disassembly. Meanwhile, the multiobjective GFA, FA, and IMOFPA are all population-based evolutionary algorithms, so the parameters, such as the population size and the maximum number of iterations, are kept consistent with those in IMOFPA and set to 100 and 200, P m max and P m min in GFA are set to 0.3 and 0.1, respectively. The confidence levels are set to 0.9, 0.95, and 0.975. Each algorithm is tested 5 times. Table 2 shows the average HV corresponding to the Pareto optimal solution at different cycle times and the minimum mean value for each objective. As can be seen from Table 2, the mean HV obtained by IMOFPA is 4.28% smaller than those obtained by GFA on average, except for 3 which are the same. The mean HV obtained by FA is all smaller than those obtained by GFA, with the smallest gap of 2.92% when CT = 105 and α = 0.9, and the largest gap is 15.2% when CT = 84 and α = 0.975, indicating the effectiveness of embedding genetic operations into FA. From the perspective of each objective, the results obtained by IMOFPA and FA are both greater than or equal to those of GFA. The results show that the performance of  the proposed GFA is better than FA and IMOFPA in solving SM-TDLBP.

C. CASE TWO
For case two, the disassembly tasks, disassembly directions, and the DPR between tasks are obtained from three different scale benchmark instances of assembly lines: two small-scale instances, P16 is taken from Lee et al. [48] and P24 from Kim et al. [49], one medium-scale instance, P65 is taken from Lee et al. [48], two large-scale instances, P148 is taken from Bartholdi [50] and P205 is from Lee et al. [48].
The disassembly time of each task in each model follows a normal distribution. The model types and disassembly times are from Özcan et al. [41], where the operating time for each task in each model is assumed as the mean disassembly time. The variances are randomly generated using the method from Özcan et al. [38]. That is, the low variance of each disassembly task is generated between 0 and µ 2 im /4, and the high task time variance is between 0 and µ 2 im /2. The hazard degrees of the hazard tasks are obtained from Liang et al. [9].
Five multi-objective algorithms are chosen for comparison: particle swarm optimization (PSO), artificial bee VOLUME 9, 2021 Before the comparison test, it is necessary to determine the reasonable parameter values, which is time-consuming and challenging work. The orthogonal experimental design and analysis of variance are utilized to calibrate the parameters of each algorithm. Meanwhile, all comparison algorithms, like GFA, are based on population optimization, so the common parameters such as the population size and the maximum number of iterations are only needed to be calibrated once in GFA to ensure fairness. The inertia weight coefficient, individual learning and global learning coefficients of PSO are set to 0.6, 1, and 2, respectively. the number of employed bees in ABC is the same as the population size, and the honey abandonment limit is set at twice the population size. the crossover rate and mutation rate of NSGA2 and SPEA2 are set to 0.8 and 0.2, respectively. For each confidence level, each case with the same cycle time is tested 5 times under the low and high task time variance conditions, and a total of 2 × 2 × (6 + 6 + 5 + 7 + 10) × 6 × 5 = 4080 case datasets are obtained. Tables 3 and 4 show the mean HV obtained by each algorithm at two confidence levels, respectively. The * in the tables indicates that there is no feasible solution to the problem, and when α =0.9 and CT=326, only P65 has no feasible solution at high task time variance. When α =0.95, four infeasible solutions appear in Table 4. The reason is that as the confidence level increases, it increases the probability that the disassembly time of a single task exceeds the CT because of the effect of random task disassembly time. It can also be predicted that more infeasible solutions will appear when α is further increased.
It can also be seen in Tables 3 and 4 that in the smallscale instances, especially in P16, all algorithms can obtain more consistent Pareto optimal solutions. As the number of disassembly tasks increases, the algorithms start to show significant differences. Overall, the average HV obtained by GFA is 5.96%, 12.83%, 22.65%, 37.95% and 45.66% higher than those obtained by NSGA2, SPEA2, ABC, FA, and PSO, respectively. This verifies that the performance of the proposed algorithm and the comparison algorithms from  strong to weak is GFA>NSGA2>SPEA2>ABC>FA>PSO. The minimum mean values of the Pareto optimal solution obtained by GFA for each objective are given in Table 5.

D. CASE APPLICATION
The proposed multi-objective GFA has been validated in the above mixed-model disassembly cases of different scales. In this part, the GFA will be applied to a car disassembly line in Chongqing, China. A total of three types of EOL cars are disassembled on the two-sided disassembly line. The three types of EOL cars differ in their part configurations, e.g., sporty utility vehicles generally have nerf bars and luggage racks, and some cars have sunroofs. There are 74 tasks in the joint disassembly precedence diagram. The disassembly time for each task is measured by a stopwatch and its mean time and unbiased estimated variance are obtained by statistical analysis techniques. The hazard degree is assessed by environmental professionals, the disassembly direction and Immediate successors of each task are detailed in Table 6.
The CT is set to 540 seconds, the population size of GFA is set to 120, and the maximum number of iterations is set to 300, P m max and P m min are set to 0.3 and 0.1. The algorithm runs 5 times independently at confidence levels of 0.9, 0.95, and 0.975, and their respective approximate Pareto frontiers are obtained, as shown in Fig. 11.
At the low confidence level, a total of 53 Pareto optimal solutions are obtained, of which the minimum workstation activation cost is 14.8, the minimum smoothness index is 58.83, and the minimum total hazard index is 1.71. When the confidence level is 0.95, a total of 36 Pareto optimal solutions are obtained, and the three minimum objectives are 15.4, 51.58 and 1.72, respectively. At the high confidence level, the three minimum objectives are 15.4, 62.46, and 1.73 in a total of 61 Pareto optimal solutions. On the whole, as the confidence level increases, the workstation activation cost gradually increases due to the fact that the disassembly time per task is continuously enlarged to ensure a sufficiently high confidence level. At the low confidence level, the maximum workstation activation cost is 19.8, while the maximum cost increases to 20.8 and 22.4 at confidence levels of 0.95 and 0.975, respectively. At the same time, the smoothness index and the total hazard index have also increased accordingly. Tables 6, 7 and 8 summarize the assignments schemes of disassembly tasks at the confidence levels, respectively of 0.9, 0.95 and 0.975, for the lowest workstation activation cost.
Ten disassembly mated-stations are utilized for each of the three confidence levels. The difference is that at the confidence level of 0.9, a total of 18 workstations in 10 matedstations are utilized, with the right workstations of 6# and 10# not being utilized. Under the confidence level of 0.95, a total of 19 workstations are utilized, of which the right workstation of 10# is not utilized. Also, it can be found that there is no disassembly task for model B on the right workstation of 8# and the left workstation of 10#, implying that when models A and C are transported to these workstations on the disassembly line, the operator will be in the processing state, while when model B is transported to these workstations, the operators will be in an idle state. Obviously, this will inevitably cause a higher smoothness index. In fact, this disassembly scheme has a smoothness index of 3.63 although the workstation activation cost is the lowest of 15.4, which makes the workload fluctuation between workstations to the maximum. Similarly, when the confidence level is 0.975, a total of 19 workstations are utilized, of which the right workstation of 9# is not utilized. The three objectives for this disassembly scheme are relatively balanced, at 15.4, 64.13, and 2.96, respectively. The shop floor decision makers and production managers can choose their preferred disassembly scheme according to their actual production needs.

V. CONCLUSION
This paper investigates the mixed disassembly of largevolume EOL products with similar assembly structures on the two-sided disassembly line. A stochastic mixed-integer programming model is formulated for the uncertainty in the disassembly process, which considers the workstation activation cost, the workload smoothness index between workstations, and the total hazard index considering the order of hazardous task removal. A multi-objective genetic flatworm algorithm is developed, and the effectiveness of the model and the efficiency of the algorithm are verified by different scales of disassembly cases of mixed-model on the two-sided disassembly lines. Finally, the proposed model and algorithm are applied to the disassembly of three types of cars, which is helpful to decision-makers to select the desired disassembly schemes.
Although GFA obtains better results than GA and FA, it comes at the cost of more time consuming because of its more complex evolutionary structure. According to statistics, GFA takes about 1.21% more computation time than GA or FA alone, about 7.74% more computation time for mediumscale problems, and about 19.32% more time for large-scale problems. The interaction effect between GA and the growth, splitting and regeneration processes of flatworms can be considered, and those operations that do not significantly improve the generation of new individuals can be canceled or merged to enhance the computational efficiency.
This study assumes that each EOL model has an equal disassembly quantity, while in practice, it is obvious that the number of different models will fluctuate due to recycling policies and cost constraints. The disassembly quantity and sequence of each EOL model will affect the balance of the disassembly line. In future research, the joint optimization of mixed-model disassembly sequencing and disassembly line balancing can be considered. In addition, how to efficiently complete the selective disassembly of specific parts on the disassembly line will also be very meaningful.

APPENDIX
Notations of the stochastic MIP model.

Indices:
i, g, h : disassembly task index, i, g, h ∈ I . j, q : mated-station index, j, q ∈ J . k : side index of the disassembly line, 1 represents the left side, and 2 represents the right side. set of tasks without any immediate predecessors. P(i) : immediate predecessors of disassembly task i. K (i) : disassembly direction of task i. α : the predetermined confidence level. β : the sharing coefficient of disassembly equipment and tools in a mated-station. P m : the dynamic mutation rate. P m max : the maximum of mutation rate. P m min : the minimum of mutation rate. l/l max : current iteration/maximum iteration.
Decision variables: t s im , t f im : start and completion time of task i for model m.
x ijk : 1, if disassembly task i is assigned to the k side station of the mated-station j; 0, otherwise. Indicator variables: ψ : an infinite positive integer. y ig : 1, if task i starts earlier than g in the same station; 0, otherwise. NPR : a pair of disassembly tasks without DPR are assigned to the same workstation orderly.