Mixed-integer programming model and hybrid local search genetic algorithm for human–robot collaborative disassembly line balancing problem

ABSTRACT Human–robot collaborative technology maximises the advantages of the capabilities of humans and robots, and provides diverse operating scenarios for the remanufacturing industry. Accordingly, this paper proposes an innovative human–robot collaborative disassembly line balancing problem (HRC-DLBP). First, a mixed-integer programming (MIP) model is devised for the HRC-DLBP to minimise the number of workstations, smoothness index, and various costs. Second, a hybrid local search genetic algorithm (HLSGA) is developed to solve the proposed HRC-DLBP efficiently. According to the problem characteristics, a four-layer encoding and decoding strategy was constructed. The search mechanism of the local search operator was improved, and its search strategy was adjusted to suit the genetic algorithm structure better. Furthermore, the accuracy of the proposed MIP model and HLSGA is verified through two HRC-DLBP examples. Subsequently, three HRC-DLBP examples are used to prove that the HLSGA is superior to five other excellent algorithms. The case of the two-sided disassembly line problem reported in the literature is also solved using the HLSGA. The results are found to be significantly better than the reported outputs of the improved whale optimisation algorithm. Besides, HLSGA also outperforms the results reported in the literature in solving EOL state-oriented DLBP. Finally, the HLSGA is applied to a power battery disassembly problem, and several optimal allocation schemes are obtained.


Introduction
The rapid development of the global economy and technology has accelerated the iteration of various industrial products (Bai, Zhou, and Sarkis 2023).Furthermore, the constantly changing demands have resulted in considerably shorter product scrapping cycles (Luo, Thevenin, and Dolgui 2022;Xie et al. 2022;Castañé et al. 2023).Recycling end-of-life (EOL) products to reduce resource wastage and environmental hazards has become a waste management consensus (Rhee, Jang, and Kim. 2021;Vieceli et al. 2021).Disassembly is a key step in recycling, reusing, and remanufacturing EOL products (Gebhardt et al. 2022).The disassembly line, as a means for efficiently disassembling EOL products, has been promoted and applied by resource recovery enterprises (Zhu et al. 2020;Bai and Zhang 2023).However, unbalanced operation time among stations decreases disassembly efficiency and increases disassembly costs (Feng and Che 2022;Bentaha et al. 2023;Schilling and hazardous components (such as lithium batteries and circuit boards) that can considerably harm the human body and endanger life (Wu et al. 2022).With the continuous upgrade of robotics technology, the use of industrial robots to replace manual labour has considerably attracted the interest of scholars (Liu et al. 2018;Zeng et al. 2022).However, the industrial application of robots also faces a difficulty.Although industrial robots have high disassembly efficiency, they cannot complete special and complex tasks independently (Xu et al. 2021).In this case, human-robot collaborative (HRC) disassembly can take into account not only the flexibility of human operators and efficiency of robots to maximise the role of the disassembly line but also the protection of physical and mental health as well as the well-being of people involved in manufacturing (Guo, Zhang, and Zhang 2023).As a result, the foregoing has become a research focus in intelligent manufacturing (Hjorth and Chrysostomou 2022;Jahanmahin et al. 2022).In a survey of HRC technology (Hanna et al. 2022;Saenz et al. 2018) classified the collaboration between humans robots into four types of work scenarios -coexistence, sequential cooperation, parallel cooperation, and collaboration -based on whether humans and robots share space, time, tasks, and goals.Leveraging these four types of work scenarios contributes to the development of a personalised HRC disassembly technology.However, such exploitation introduces many complex features to the DLBP.These complications include the means of dealing with the relationship between humans and robots, the relationship between humans and robots and EOL products, and the relationship between humans and robots and disassembly lines; the foregoing considerably increases the difficulty of solving the DLBP.Therefore, the introduction of HRC technology to the DLBP opens a new processing mode for the intelligent remanufacturing of EOL products.The study of this technology is not only important but also challenging.
However, to the best knowledge of the authors, previous studies have not integrated HRC techniques considering different collaboration scenarios into the disassembly line, as demonstrated by the literature review presented in the next section.With these points of interest, this study proposes an HRC-DLBP that considers different collaboration scenarios and formulates minimisation goals.The foregoing considers the number of workstations, smoothness index, and disassembly cost (including the robot energy consumption cost, robot disassembly cost, manual disassembly cost, fixed asset input cost, and additional disassembly cost for hazardous tasks).The contributions of this study to existing knowledge are as follows.
• A disassembly line considering different HRC scenarios is designed.• A mixed-integer programming (MIP) model for the HRC-DLBP is formulated to solve small-scale problems and verify the accuracy of the proposed metaheuristic algorithm.• A hybrid local search genetic algorithm (HLSGA) is developed; encoding and decoding are designed based on the features of the HRC-DLBP.Moreover, a genetic operator is effectively combined with a local search operator to enhance the optimisation ability of the algorithm.
• The accuracy and superiority of the proposed HLSGA are demonstrated by applying it to five test cases.• The proposed algorithm is applied to the disassembly of a scrap power battery module.
The remainder of this paper is organised as follows.Section 2 reviews previous research on DLBP.Section 3 describes the proposed HRC-DLBP and MIP model construction.The proposed HLSGA is introduced in Section 4. Section 5 presents the verification of the correctness of the MIP and the accuracy and superiority of the HLSGA through different case calculations.Section 6 discusses the application of the proposed algorithm to a disassembly example of a power battery module.Finally, Section 7 summarises the conclusions and future research directions.

Literature review
This section reviews and discusses reports in the literature to identify knowledge gaps in existing studies based on four aspects: disassembly operator, representation method of disassembly precedence relations, DLBP optimisation objective, and DLBP optimisation method.

Disassembly operator
Based on different operators, existing DLBPs can be classified into manual disassembly, robotic disassembly, and human-robot cooperative disassembly (Özceylan et al. 2019;Laili et al. 2020).Robotic disassembly compared with manual disassembly is attracting interest in the industry and academia because of its higher operating efficiency and lower operating costs.As shown in Table 1, Liu et al. first proposed the robot DLBP, which generated feasible disassembly sequences through spatial interference matrix.They also improved the discrete bee colony algorithm to provide solutions for robot disassembly lines (Liu et al. 2018;Liu et al. 2018).Subsequently, they included the disassembly sequence planning problem in their investigation to improve the disassembly efficiency (Liu et al. 2020).Fang et al. (2019) formulated a mathematical programming model to minimise the cycle time, energy consumption, peak energy consumption, and number of robots as well as developed an evolutionary algorithm to solve these problems.Based on the foregoing, Fang et al. (2020aFang et al. ( ,2020b) ) considered the robot resource constraints and interval processing time of a multi-robot DLBP.Çil, Mete, and Serin. (2020) presumed that the disassembly time of different robots to complete the same task varies.Accordingly, they developed a mixed-integer linear programming model and an intelligent algorithm based on ant colony optimisation to solve small-scale and large-scale problems, respectively.However, with respect to the current robotics technology, robotic disassembly is a highly customised operation method.If the residual value of the products to be disassembled is not high or the number of recovered products to be disassembled is unstable (the quantity is small), the recovery cycle of enterprise hardware input costs inevitably increases (Hanna et al. 2022).However, if the structure of the EOL product to be dismantled is extremely complex, the robot cannot complete the disassembly task, and the product is processed by flexible manual disassembly (Wu et al. 2022).Consequently, HRC disassembly has attracted the interest of researchers.In 2020, Xu et al. (2020) fully considered task failures based on the difficulty of part disassembly, assigning difficult disassembly tasks and easy tasks to humans and robots, respectively.They formulated four objectives: maximising profits, minimising energy consumption, reducing disassembly difficulty, and decreasing the number of workstations; they subsequently achieved these goals using artificial bee colony algorithm.In 2021, Xu et al. (2021) proposed an HRC-DLBP based on a safety strategy classification tree.They assumed that one worker and one robot were assigned to each workstation and collaborated to complete the disassembly task.Then, the proposed discrete bee colony algorithm was verified via a bearing coupling case study.Wu et al. (2022) assigned a worker or robot to a workstation according to task attributes and studied the HRC-DLBP for Tesla battery modules.

Representation method for disassembly precedence relations
Disassembly priorities are caused by the physical interference among the parts of scrap products, and priority relationship constraints must be satisfied to ensure the feasibility of the disassembly sequence.Presently, mainstream methods for representing priority relationships are broadly divided into two categories.The first category deals with AND/OR priorities.Among the typical considerations in this category are the AND/OR Graph (Homem de Mello and Sanderson 1990;Tian et al. 2019), transformed AND/OR Graph (TAOG) (Koc, Sabuncuoglu, and Erel 2009;Ren et al. 2018;Chen et al. 2022) and disassembly trees (DT) (Koc, Sabuncuoglu, and Erel 2009).This type of method can fully express the connection state of parts; however, it is extremely complicated when expressing the priority relationship of largescale products (Yin et al. 2022).In contrast, the second category ignores the OR relationship and only considers the AND relationship.It facilitates the expression of large-scale products, such as the task precedence diagram (TPD) (Koc, Sabuncuoglu, and Erel 2009) and part precedence diagram (PPD) (Zhang et al. 2017;Zhu, Zhang, and Wang 2018;Wang, Li, and Gao 2019;Yin, Zhang, and Jiang. 2021), The TPD only focuses on the priority relationship among disassembly tasks and cannot express the current state of disassembly information of parts.In contrast, the PPD can overcome this deficiency.Accordingly, this study focuses on the PPD.

DLBP optimisation objectives
In DLBP research, the following are typically targeted for optimisation: number of workstations, smoothing index, demand index, hazard index, and operating energy consumption of disassembly lines (Agrawal and Tiwari 2008;Gungor and Gupta 2001;Kalayci and Gupta 2013;Kalayci, Polat, and Gupta 2015;Kalaycılar, Azizoğlu, and Yeralan 2016;Mete et al. 2016;Altekin 2017;Ren et al. 2018;Özceylan et al. 2019;Fang et al. 2019;Wang, Li, and Gao 2019;Laili et al. 2020;Xu et al. 2020;Cevikcan, Aslan, and Yeni 2020;Edis 2021;Wang et al. 2021;Çil et al. 2022;Edis, Edis, and Ilgin 2022).The number of workstations determines the disassembly line length.Decreasing the number of workstations not only reduces the construction cost but also effectively reduces the area of the workshop to reduce the cost of industrial land (Wang et al. 2020;Liang et al. 2022;Guo et al. 2023).The smoothing index is the main method for balancing the load among the stations to improve disassembly efficiency.These two indicators are the main evaluation indices in the study of the DLBP (Zhu, Zhang, and Guan 2020;Wu et al. 2022).The demand and hazard indices require the early disassembly of high-value and hazardous parts to avoid the loss of disassembly profits (Zhang et al. 2021).The disassembly profit of the resource recovery company is also affected by the daily operating cost, which is an extremely expensive expense; it has been incorporated into the lean production management of the enterprise.The foregoing shows that in addition to the main objective, the disassembly cost must be considered more comprehensively.

DLBP optimisation method
The methods commonly used to solve DLBPs include exact and inexact methods (Edis, Edis, and Ilgin 2022).The advantage of the former over the latter is that it is more effective in solving small-scale problems and can be used to verify the accuracy of the inexact method.Typical exact methods include integer programming (Kalaycılar, Azizoğlu, and Yeralan 2016;Zhu et al. 2020), MIP (Altekin, Kandiller, and Ozdemirel 2008;Altekin and Akkan 2012;Zhang et al. 2021;Yin et al. 2022), dynamic programming (Koc, Sabuncuoglu, and Erel 2009), constraint programming (Çil et al. 2022), and chance constraints planning (Bentaha, Battaïa, andDolgui 2013, 2015).Mcgovern and Gupta (2007) have demonstrated that the DLBP is non-deterministic polynomial (NP)hard; consequently, exact methods cannot solve largescale problems.Furthermore, an accurate model is difficult to build when the DLBP constraints are extremely complex (Edis, Edis, and Ilgin 2022).
Inexact methods refer to heuristic and meta-heuristic algorithms, which can obtain a satisfactory solution within a certain period of time.Heuristic algorithms are typically designed with search strategies based on problem characteristics (Kalaycılar, Azizoğlu, and Yeralan 2016), which exclude random factors; given the input, these algorithms have fixed outputs.The metaheuristic algorithm has been widely studied; it is random and yields solutions with better diversity.Typical meta-heuristic algorithms include genetic algorithm (McGovern and Gupta2007), tabu search (Kalayci and Gupta 2014), simulated annealing algorithm (Kalayci, Gupta, and Nakashima 2012), ant colony optimisation algorithm (Kalayci and Gupta 2013), particle swarm optimisation algorithm (Kalayci and Gupta 2013), artificial fish swarm algorithm (Wang et al. 2017;Zhang et al. 2017), artificial bee colony algorithm (Kalayci and Gupta 2013;Wang et al. 2021), discrete flower pollination algorithm (Wang, Li, and Gao 2019).To combine the advantages of different algorithms, some hybrid algorithms have been developed and reported to exhibit improved optimisation performance.Some of these hybrids include genetic simulated annealing (Wang et al. 2021), hybrid group neighbourhood search algorithm (Zhu, Zhang, and Guan 2020), hybrid drive algorithm (Yin et al. 2022).The foregoing inspired the development of a new powerful hybrid algorithm for solving the HRC-DLBP-the HLSGA.

Research gap
The literature review indicates that the papers on HRC disassembly in the field of DLBP are few.Three studies on human-robot collaboration, listed in Table 1, have three deficiencies, as follows.
• The consideration of scene application is insufficient.Saenz et al. (2018) classified HRC scenarios into four categories: coexistence, sequential cooperation, parallel cooperation, and collaboration.The three studies did not fully consider collaboration scenarios.• The formulation of the optimisation objective is unreasonable.Only Wu et al. (2022) considered smooth metrics; however, they did not consider the energy consumption of robots.None of the three studies considered the disassembly line operating costs in all aspects.
• The development and application of solutions are insufficient.In terms of method development, the three studies neither established exact mathematical models to verify the accuracy of the proposed algorithm (decoding errors typically result in infeasible solution sequences) nor extensively compared the devised algorithm with existing ones to prove its superiority.In terms of method application, the object of the case study was a single small-scale EOL product that could not prove the effectiveness of both the meta-heuristic algorithm they proposed for solving large-scale problems and the proposed method for disassembling other EOL products.
Accordingly, the motivation of the present study is to fill the knowledge gap in the foregoing aspects.

Problem description
As shown in Figure 1, a worker and robot are assigned to a workstation.The human-robot collaboration can be divided into four scenarios: coexistence, sequential cooperation, parallel cooperation, and collaboration.
• Coexistence means that humans and robots do not share a workspace; hence, they are not in the same workstation although they coexist.• The sequential cooperation mentioned by Saenz et al. (2018) refers to situations in which humans and robots occupy the same space at different times.This is a collaboration scenario proposed for the case of a single workstation; it is frequently used in the study of disassembly sequence planning (Hanna et al. 2022).For the disassembly line, humans and robots typically only move within the station to avoid excessive travel time, which affects disassembly efficiency.Accordingly, in this study, the relationship among the operators of adjacent workstations is defined as sequential collaboration.
• Parallel cooperation means that people and robots are in the same workspace at the same time; operators have no contact.As observed in workstations 2 and 3, hazardous tasks are assigned to robots, and complex tasks are assigned to humans during operation.This cooperation scenario not only reduces the risk of injury to people but also increases space utilisation and shortens the duration of cyclic operation.• As shown in workstations 1 and 4, collaboration means that people and robots must share a working time (i.e. they must start and end simultaneously) when they share the workspace and collaborate to perform the same task.In this scenario, human-robot interaction (i.e.contact) is necessary.Moreover, human and robots in the collaborative scene do not interact when performing non-interactive disassembly tasks, but complete disassembly tasks independently.Thus, the tasks defined in this study for interactive disassembly must not be hazardous or complex.
Note that in this study, tasks that are both hazardous and complicated are defined as complex tasks since robots cannot complete complex tasks.The method of assigning tasks based on their attributes is referred to as the task attribute constraint.
The HRC-DLBP must also satisfy priority relationship and cycle time (CT) constraints.In a station, suppose a task is assigned to the human, while some of its immediately preceding tasks are assigned to the robot.In that case, the task cannot be performed unless the previous tasks assigned to the human have been completed and the tasks' immediately preceding disassembly tasks have also been performed.When the interactive cooperation task begins, the previous tasks of both operators must be complete.The traditional CT constraint means that the total disassembly time for each workstation must not exceed CT to maintain the efficient operation of the disassembly line.In contrast, a worker and robot present at the same time in one workstation are considered in this study; hence, the total disassembly time for each workstation cannot exceed 2CT.
The following are assumed in this study.
• The number of products to be disassembled is sufficient, and the task mode is complete disassembly.• The required time for completing parts is ignored.
• The tool replacement time of workers or robots is included in the disassembly time of parts.• All parts have no physical damage, and destructive disassembly is not allowed.• All workstations are equipped with auxiliary lifting equipment to overcome the restrictions in the disassembly direction of parts.
• Each worker has the same level of operating proficiency.

Mathematical model
The indexes, parameters and variables required for the MIP model are as follow: Based on the foregoing notations, the MIP formulation of the HRI-DLBP can be expressed as follows.Objective functions are: (2) (3) x im if tasks i is assigned to the workstation m, x im = 1; otherwise, x im = 0 y iw if tasks i is assigned to the operator w, y iw = 1; otherwise, y iw = 0 z ij if tasks i and j are assigned to the same operator at the same workstation and i is removed before j, z ij = 1; otherwise, j, if workstation m is tuined on, S m = 1; otherwise, S m = 0 P wm if operator w is assigned to the workstation m, P wm = 1; otherwise, The mathematical model formulated three objectives, as follows.The first is given by Equation (2); it aims to minimise the number of workstations required to optimise the line length.The second is given by Equation (3); it minimises the load smoothness among operators to increase disassembly efficiency.The third is given by Equation ( 4); it minimises the disassembly cost (including robot energy consumption cost, labour cost, fixed asset input cost, and additional disassembly cost) for hazardous tasks.The robot energy consumption cost is given by Equation ( 5), which calculates the standby and running energy consumptions of the robot.The labour cost includes the disassembly cost of the robot, and the worker cost is calculated using Equation ( 6).The fixed asset input cost includes the fixed cost of configuring a workstation and that of purchasing a disassembly robot, as calculated by Equation ( 7).Parts with hazardous properties can be detrimental to the environment and must be treated to become harmless.The additional processing cost of hazardous tasks is calculated using Equation ( 8).The following constraints are given according to the problem features and objective function.
The constraints are as follows: Task assignment constraints are represented by Equations ( 9)-( 13).Equation ( 9) indicates that all tasks must be assigned to workstations, and the same task can only be assigned to one station.Equations ( 10) and ( 11) restrict the assignment of non-interactive tasks to robots or workers, whereas robots and humans work together on interactive tasks.Finally, Equations ( 12) and ( 13) require that complex tasks must be assigned to humans, and hazardous tasks must be assigned to robots.
Time-dependent constraints are represented by Equations ( 14)-( 19).Equation ( 14) provides the upper and lower bounds of the number of workstations based on the sum of all task times (CT) and number of tasks.Equation (15) indicates that if task i is the predecessor of task j and tasks i and j are assigned to the same operator at the same station, then the start time of task j must be later than the end time of task i. Equations ( 16) and ( 17) stipulate that for any task, the start and end times must be in the CT of a workstation to avoid the operation of tasks across stations.Equation ( 18) restricts the start time of any predecessor task to a time later than the end time of its predecessor task.
The location-dependent constraints restricting the priority relationship between tasks i and j are represented by Equations ( 20)-( 24).Equation ( 20) establishes the relationship between tasks, operators, and workstations, and Equations ( 21)-( 24) ensure that tasks i and j satisfy the disassembly priority relationship.
The workstation allocation constraints are represented by Equations ( 25)-( 29).Equation (25) indicates that when task i is assigned to workstation m and task i is assigned to operator w, operator w must be assigned to workstation m.Equation ( 26) restricts the use of a workstation to a single pair of worker and robot.Equations ( 27) and ( 28) provide the lower and upper bounds of the workstations by task assignment.Equation ( 29) indicates that the workstations must be turned on sequentially.The binary variables are given by Equations ( 30)-(34).

Model linearisation
The solution to the nonlinear MIP model is extremely difficult to obtain, and linearising the model is necessary.Equations ( 3), ( 5), and ( 6) contain the product of the following two decision variables: x im and y iw .Let G imw = x im • y iw ; hence, Equations ( 3), ( 5), and ( 6) can be expressed as: Next, introduce the linearisation constraints as follows: Equations ( 38) and ( 39) indicate that the decision variable G imw must not exceed the decision variables x im and y iw , thus giving an upper bound for G imw .Equation ( 40) further restricts G imw to be equal to 1 when both x im and y iw take a value of 1.

Hybrid local search genetic algorithm
The non-dominated sorting genetic algorithm (NSGA-II) (Deb et al. 2002) exhibits excellent global optimisation characteristics for solving multi-objective problems.However, its local search ability is inadequate, and its search efficiency is low.The local search algorithm compared with NSGA-II has a strong local search ability; however, it is extremely greedy when searching locally or in a neighbourhood.Consequently, the algorithm facilely falls into local optima (Li, Kucukkoc, and Zhang 2019).
In view of the foregoing, this section presents the combination of the advantages of these two algorithms to develop a multi-objective HLSGA by introducing local search operators under the NSGA-II algorithm framework.The HLSGA flowchart is shown in Figure 2, and the optimising process of HRC-DLBP by HLSGA is as follows: Step 1: Algorithm parameter settings, including population size N pop , running time T r , crossover probability P c , mutation probability P m , local search probability P s and CT of workstation.
Step 2: Population initialisation.The proposed encoding method is used to create an initial population that satisfies the disassembly priority relation.Then, the proposed decoding method is used to assign tasks to the workstations, and the optimisation objective is calculated according to the assignment results.Finally, the Pareto method is used to screen non-dominated solutions and update external files N e .
Step 4: Perform crossover operation.Cross the best individuals from external files with the initial population to generate new individuals.
Step 5: Perform mutation operation.Mutate the best individuals from external files and the individuals generated from crossover operation to generate new individuals.
Step 6: Perform improved local search operation.To generate new individuals, utilise crowding distance to select the 20 best individuals from external files and perform a local search operation on them.
Step 7: Update solution set.After decoding all new individuals, the objective values of all individuals are obtained.
Step 8: Update external files N e .The Pareto method is applied to the new solution set to filter and obtain non-dominated solutions, which are subsequently used to update the external archive.
Step 9: Generate the next-generation parent population.Check the number (N o ) of individuals in the external file.If N o > N pop , perform a non-dominated sort on the individuals in the external file; otherwise, use the proposed encoding method to regenerate insufficient individuals.
Step 10: Termination condition check.The algorithm ends if the termination condition (running time) is satisfied.Otherwise, return to step 3 to continue execution.
Step 11: Output optimal solution set.
The coding and decoding operations, genetic operators, local search operators, and multi-objective evaluation methods of the proposed HLSGA are elaborated as follows.

Encoding and decoding
Because the HRC-DLBP is a discrete problem, integer number coding based on the task sequence is used in this study.As shown in Figure 3, the priority relationship diagram is converted into a priority relationship matrix to obtain the sequence code satisfying the priority relationship.Then, a group of feasible disassembly sequences is generated through matrix element transformation.The specific matrix element transformation method was described by Zhu et al. (2020).
The encoding operation generates a set of task sequences satisfying the precedence constraints, whereas the decoding operation assigns this set to operators in the workstation.This satisfies the beat time and operator's task attribute constraints described in Section 3.1.As shown in Figure 4, a four-layer decoding method is employed to ensure the uniqueness of the solution.The first level is the task sequence represented by a real number denoting the task index.The second layer is the workstation number corresponding to the task sequence; it is represented by a real number that denotes the workstation index.The third layer is the operator number corresponding to the task sequence and workstation number; it is represented by the operator-type index: 1 and 2 denote a worker and robot, respectively.The fourth layer is the cooperation mode corresponding to the workstation and operator numbers; 1 and 2 represent the collaboration and parallel cooperation, respectively.
Two heuristic decoding rules were established for objectives 2 and 3 of the formulated mathematical model.
• For objective 2, a current non-attribute task is preferentially assigned to the operator with more idle time to reduce the idle time.• For objective 3, when the two operators have the same remaining time, tasks are assigned to the robot first to reduce the standby energy consumption of the robot.
The foregoing can effectively reduce the disassembly cost because the disassembly cost per unit time of the robot is lower than that of the worker.
Notably, the workers and robots in the same station work at the same time.Therefore, in assigning tasks to two operators, recording the start and end times of each task is necessary to satisfy the task priority relationship and CT constraints.The specific decoding process is shown in Algorithm 1, In the algorithm, X represents the feasible disassembly task sequence; m represents the current workstation number; WTM represents the current working time node of the worker; RTM represents the remaining working time of the worker; WTR represents the current working time node of the robot; and RTR represents the remaining working time of the robot.Section 3.2 presents the meaning of the other parameters.In steps 1-11 of Algorithm 2, a hazardous task is assigned to the robot.However, determining whether the current task has a predecessor task is necessary; if such is the case, the current task can only be started after the predecessor task ends.In steps 12-23 of Algorithm 3, complex tasks are allocated to humans; the task priority relationship among different operators must also be observed.In steps 1-11 of Algorithm 3, interactive tasks are simultaneously assigned to humans and robots.The start and end times for the human and robot pair to complete the interactive tasks must be consistent.In steps 12-27 of Algorithm 3, the two heuristic decoding rules are formulated.

Algorithm 1 Decoding
Require: A disassembly sequence Ensure: Obtain a disassembly scheme that satisfies various constraints of HRC-DLBP. 1: Executing Algorithm 2 5: end if 6: if X(i) ∈ NI then 7: Executing Algorithm 3 8: end if 9: end for Algorithm 2 Decoding Require: Assign hazardous and complex tasks.Ensure: Hazardous and complex tasks are assigned to robots and workers, respectively.As shown in Figure 5, the task sequence in Figure 4 has been allocated to two workstations using the proposed four-layer encoding and heuristic decoding method.As workstations 1 and 2 have interactive tasks (6 and 1, respectively), no parallel collaboration mode exists for humans and robots in this disassembly scheme.Since interactive disassembly task 6 should begin and end simultaneously, the robot at the first workstation must wait for the worker to complete task 5 before commencing the operation.Complex tasks must be performed manually, so task 9 must be allocated to the human operator at the second workstation.Tasks (2,3,4,5,7,8,10) without specific attributes are assigned based on heuristic decoding rules, provided that they meet the priority relationship for disassembly.

Crossover and mutation
Task priorities must be considered in crossover and mutation operations to derive a feasible solution.As presented in this section, the crossover operation ensures that the solution sequence is feasible through the mapping crossover between two feasible sequences, whereas the mutation operation limits the interval of mutation points based on the priority relationship.The two crossing points are shown in Figure 6 (a), two two points, P 1 and P 1 , are randomly generated on the parent sequence.Using these two points as boundaries, a gene segment [6,5,3,7,4,9] can be obtained.Subsequently, a new gene segment [4,5,7,9,6,3] can be obtained by mapping the gene segment in the optimal solution sequence.Finally, a new solution with an excellent parent gene is obtained by replacing the original gene segment with a new gene segment.Single point variations are shown in Figure 6 (b).Task 5 is randomly selected as the variation point of the optimal solution using probability.After evaluating the priority relationship, task 5 must be disassembled after and before tasks 1 and 10, respectively; then, the variability interval of task 5 is [4,7,9,6,8].New solutions can be obtained by randomly disassembling the tasks into this interval.

Improved local search
The local search operator typically requires the fitness value calculation to provide feedback on the search effect.However, this search method is time-consuming when solving HRC-DLBP because the complex nature of the problem decelerates the decoding process.To improve the solution speed of the proposed HLSGA further, the stopping conditions and search methods of the local search are improved to adapt to the structural characteristics of the genetic algorithm.First, the Pareto solution set is filtered by crowding distance such that the number of solution sequences used in the local search operation does not exceed 20 to improve the computational efficiency.Second, the local search operation in the HLSGA does not calculate the fitness value but retains all the searched solution sequences and returns them to the external archives to adapt to the structural characteristics of the genetic algorithm.The improved local search is presented in Figure 7. Task 7 is randomly selected as the starting point for the search.After assessing the priority relationship, it is found that the movable range of task 7 is [4,1,9,10], and task 7 generates a new set of solutions through partial movement.The improved local search is similar to that of the single-point mutation operation.Although the search depth is enhanced, the efficiency of the algorithm is not reduced.

Multi-objective evaluation index
Hypervolume (HV) refers to the space volume formed by the pre-set reference points distributed in target space as the boundary and Pareto solution set obtained by the algorithm.To evaluate the convergence and diversity of the algorithm, HV metrics are typically used.The algorithm converges when the HV value tends to be stable.The higher the HV value, the larger the target space volume formed by the obtained solution; this proves that the diversity and quality of the Pareto solution set solved by the algorithm are better.The formula for calculating the HV index is as follows: where λ indicates the Lebesgue measure, for measure hypervolume; |S| indicates the number of solutions in the Pareto solution set; and v i indicates the hypervolume formed by the ith non-inferior solution and the reference point.

Verification of model and algorithm
This section presents the verification of the proposed model and algorithm through numerical tests in four parts.The first part verifies the accuracy and effectiveness of the MIP and HLSGA, and the second part confirms the superiority of the HLSGA in solving the HRC-DLBP at different scales.Because no test benchmark exists for the proposed HRC-DLBP in the existing knowledge base, the reported test cases (P10, P25, and P55) in the literature regarding these two parts are modified.Test cases P10 (a personal computer) and P25 (a mobile phone) are modified from Kucukkoc (2020), and P55 (a printer) is modified from Wang, Li, and Gao (2019).The priority relationship diagrams for these three cases are presented in Figs S1, S2 and S3 of Supplementary data (see Data availability statement).The third part computes and compares the two-sided DLBP (TDLBP) for an engine example (P34); its 34 parts are reported in the literature (Zhang et al. 2021).The fourth part again verifies the superiority of HLSGA with the EOL state-oriented DLBP example; its 55 parts are reported in the study by (Zhu et al. 2020).The algorithm program used in the research is developed using MATLAB R2020b; the operating system is Windows 11 Pro, and the central processing unit of the computer is Intel i9-10900K at 3.7 GHz with a 32-GB RAM.

Correctness and validity of model and algorithm
This section presents the use of two small-scale cases (P10 and P25) for verifying the accuracy and effectiveness of the MIP model and proposed HLSGA by comparing their calculation results.The MIP model was developed based on a free academic version of GUROBI 9.1.5.The known parameters required by MIP and HLSGA were referred to existing research (Wu et al. 2022).Specifically, these parameters include the additional cost of dealing with hazardous tasks, C h = 5 RMB; cost per unit time of robots, C rt = 5 RMB/h; cost per unit time of workers, C m = 37 RMB/h; the operation of robots power consumption, OE = 5 kW • h; and robot standby power consumption, SE = 0.2 kW • h.The CT of P10 and P25 are set to 62 and 29 s, respectively.This study establishes three objectives for the proposed HRC-DLBP.Objective 2 is a nonlinear (exponential) objective that prevents the exact solver from simultaneously solving multiple objectives.Because the real Pareto front contains the single-objective optimal solution, the single-objective optimal solution of the MIP model obtained using the GUROBI solver is compared with the HLSGA solution results.The parameter settings for the HLSGA are listed in Table S1 of Supplementary data.The algorithm program was independently run 10 times, have solved the Pareto solution set containing the single-objective optimal value with each run.Table 2 lists all the solutions obtained using the GUROBI solver and HLSGA.The table indicates that the two methods solve the P10 and P25 cases and obtain three solutions containing the single-objective optimal solution.In addition, when computing smallscale cases, there is little difference in the computation time between the MIP model and HLSGA.However, as the size of the problem increases, the time required for GUROBI to solve the MIP becomes significantly higher than HLSGA.For example, the solution of objective 3 for the P25 case took 59853.04s,which was much longer than HLSGA.
Figures 8 and 9 display the solutions containing the optimal single-objective solution obtained by two methods for solving the P10 and P25 cases, respectively.The different task categories in these two figures are marked using different colours.Consider Figure 8 as an example.The light purple, light yellow, and light blue colours represent the complex task (task 9), harmful task (task 7), and human-robot interaction task (tasks 1 and 6), respectively.According to the division of collaboration scenarios described in Section 3.1, the six solutions of the P10 case only include coexistence, sequential cooperation, and collaboration scenarios.However, the parallel cooperation scenario is excluded because each workstation has tasks that require human-robot interactive collaboration.All six solutions of the P25 case compared with those of the P10 case include four HRC scenarios.As shown in Figures 8 and 9, the solution results of both the MIP model and HLSGA satisfy the constraints of the HRC-DLBP proposed in Section 3.1.Notably, as shown in the schematic of HLSGA, the heuristic decoding rules proposed in Section 4.1 render the allocation scheme more reasonable, especially by assigning more tasks to the robot operator to increase the physical and mental wellbeing of the workers.Therefore, both the exact method and meta-heuristic algorithm proposed in this study are correct.

Verification of HLSGA with HRC-DLBP example
To verify the superiority of the HLSGA in solving the HRC-DLBP further, five powerful algorithms discrete flower pollination algorithm (FPA) (Wang, Li, and Gao 2019), NSGA-II (Deb et al. 2002), genetic simulated annealing algorithms (GASA) (Wang et al. 2021), multi-objective improved particle swarm optimisation (MIPSO) (Li et al. 2018) and goal-driven discrete cuckoo search (GDCS) (Li et al. 2018) are selected for comparative testing.The proposed HLSGA and five selected algorithms are used to calculate the three cases.The CT settings of P10, P25, and P55 are 62, 29, and 60 s, respectively.The different algorithm parameter settings are presented in Table S1 of Supplementary data.The parameters of the proposed HLSGA are determined after numerous calculations and tests, and those of the other algorithms are based on original papers published in the literature.
Each algorithm is independently run 10 times.Each running time is used as the algorithm stop condition; the solution times of P10, P25, and P55 are set to 1, 10, and 1000 s, respectively.The highest HV value of each algorithm (the best set of Pareto solutions) is selected for comparison with other algorithms.The results of each algorithm for P10, P25, and P55 are listed in Tables S2, S3 and S4 of Supplementary data, respectively.The Pareto solutions in the table are indicated by bold numbers.Figures 10(a,b,c) show the solution space distribution of each algorithm in Tables S2, S3 and S4 To evaluate the stability of the proposed HLSGA, the HV values of the results of 10 independent runs of the different algorithms were calculated.The reference points of the HV index for the P10, P25, and P55 cases are (3, 2000, 1.6e+5), (5, 1500, 3.5e+5), and (12, 5000, 6.5e+5), respectively.Note that the HV index of each algorithm is calculated using all the Pareto solutions obtained.After counting the HV calculation results 10 times, boxplots of P10, P25, and P55 cases are drawn, as shown in Figures 11 (a,b,c), respectively.The stability of the HLSGA is observed to be better than that of the other five algorithms for different problem sizes.Moreover, a higher HV value means that the solution obtained by the algorithm is closer to the real Pareto front; again, this proves that the HLSGA has better optimisation ability than the other algorithms.

Verification of HLSGA with two-sided DLBP example
To verify the solving performance of the proposed HLSGA further, an engine disassembly case (P34) on a TDLBP with 34 tasks is introduced, and the results are compared to those in the recently reported research  (Zhang et al. 2021) to verify the superiority of the HLSGA.The P34 case optimisation objectives include the number of workstations (f 1 ), smoothing index (f 2 ), demand index (f 3 ), and hazard index (f 4 ).After extensive testing, N pop = 135, P c = 0.9, P m = 0.1 and P s = 0.2 are set as the HLSGA parameters, and the algorithm deadline is 120 s; the calculation results are listed in Table 3. Eight Pareto solutions are obtained by the improved whale  :22,23,18,9,12,11,17→20→28,14,13,25→2,24,26,8,27,34 R:3,5,6,7,10,19,15→4,1,16→29→21,30,31,32 optimisation algorithm (IWOA), whereas 10 Pareto solutions are derived by the HLSGA.Specifically, the first and second solutions of the IWOA are dominated by the second and eighth solutions of the HLSGA.In addition, the IWOA does not obtain solutions with small values for f 2 , f 3 and f 4 , indicating that the solution distribution and solution depth obtained by the HLSGA are better than those of the IWOA.Therefore, the proposed HLSGA also outperformed the IWOA in solving the TDLBP.

Verification of HLSGA with EOL state-oriented DLBP example
In order to verify the superiority and versatility of the proposed algorithm to solve the classic DLBP, this subsection introduces a large-scale printer disassembly case with 55 parts, and compares the results with the existing literature.Zhu et al. (2020) formulated three optimisation objectives for the proposed EOL state oriented DLBP: number of workstations (F 1 ), idle time (F 2 ), and amount of disassembly resources (F 3 ).After extensive testing, N pop = 400, P c = 0.9, P m = 0.12 and P s = 0.25 are set as the HLSGA parameters, and the algorithm deadline is 1800 s.Table 4 lists the solution results (including objective value and solution time) of the basic migrating birds optimisation (MBO), hybrid migrating birds optimisation (HMBO) proposed by Zhu et al. (2020) and HLSGA proposed in this study.Table S5 (see Supplementary data) lists the disassembly schemes obtained by HLSGA.As can be seen from the Table 4, HLSGA and HMBO significantly outperform MBO in their optimisation capabilities.Furthermore, HLSGA obtains one more Pareto solution (7, 206.5450, 13) than HMBO despite taking the least time.This once again proves that the proposed HLSGA has strong superiority and versatility in solving similar discrete problems.

Instance application
Recently, the cascade utilisation of waste power batteries has become a research focus (Harper et al. 2019;Xu et al. 2022).The first step in the cascade utilisation of EOL products, such as waste power batteries, is disassembly.Accordingly, this section introduces the disassembly case of a Tesla power battery module; the HRC-DLBP proposed in this study is applied to this case.
In the case application, the proposed HLSGA is compared with FPA, NSGA-II, GASA, MIPSO, and GDCS.The known parameters of each algorithm are the same as those described in Section 5.1: CT = 200, and the algorithm stop time is 500 s; the other parameters are listed in Table S6 of Supplementary data.

Instance information
The disassembly case information of the Tesla power battery module is obtained from the report of Wu et al. (2022), The three-dimensional structure, shown in Figure S4 of Supplementary data, has 44 disassembly tasks.The names of the parts shown in the figure are summarised in Table S7 of Supplementary data.The table also lists the disassembly time of each task, task attributes, priority relationships, and collaborative disassembly task.Note that because the HRC disassembly could improve the disassembly efficiency, the disassembly time of the collaborative disassembly task was modified in this study to adapt to the proposed HRC-DLBP.

Results of instance application
Each algorithm is executed 10 times, stops after 500 s of runtime, and selects the best results for comparison.
Because each algorithm obtains dozens of Pareto solutions, 10 solutions are filtered using crowding distance, as listed in Table 5.When the results of each group of algorithms are combined for a unified comparison, the solution set of the HLSGA completely dominates the solutions of the other five algorithms.This means that after the solution space increases with the problem scale, the other algorithms become unable to derive better solutions in a short time.This proves that NSGA-II has better search ability over a short time than FPA, NSGA-II, GASA, MIPSO, and GDCS.
To evaluate the search and convergence abilities of the proposed HLSGA in various aspects, the changes in the HV index when each algorithm is run are tracked.The HV value is recorded every second, and an iterative curve is drawn, as shown in Figure 12.In the first 10 s, the HV values of the algorithms (other than the FPA) rapidly increase.Within 10-50 s, the HV values of NSGA-II, GASA, MIPSO, and GDCS gradually tended to stabilise.In contrast, the HLSGA continued to search the solution space and did not stabilise until 131 s.Based on the HV value, the proposed HLSGA could not fall into the local optimum compared with NSGA-II, GASA, MIPSO, and GDCS.In addition, the FPA tended to stabilise at 185 s; however, it obtained the least HV value among the algorithms.This indicates that the reference point and solution it derived result in the least HV.The FPA can be considered as the worst performing algorithm in solving the HRC-DLBP proposed in this study because it has inadequate search performance and the slowest convergence.Overall, the convergence speed and optimisation ability of the HLSGA are better than those of the other five algorithms.

Results analysis of disassembly scheme
Three schemes were screened from the obtained Pareto solution set considering the varying objective preferences of different enterprise managers.As shown in Figure 13, light purple, light yellow, and light blue indicate complex tasks, hazardous tasks, and tasks accomplished by the collaboration between humans and robots, respectively; this is similar to the description presented in Section 5.1.Moreover, all these schemes satisfy the set constraints of CT, task attribute, and HRC disassembly.Notably, the number of workstations in the three schemes is optimal.Scheme 1 is a compromise scheme; that is, objectives 2 and 3 are neither the best nor the worst, and the evaluation indicators in all aspects are relatively balanced.Objective 2 (smoothing index) of scheme 2 is optimal, and the difference among the idle times of the stations is the smallest, balancing the workload among the operators and not readily causing blockage of the disassembly line.In contrast to scheme 2, scheme 3 has the worst objective 2; however, it costs less to dismantle.This is because more tasks are allocated to the robot, and the energy consumption and daily maintenance costs of the robot are low, thereby reducing disassembly costs.Furthermore, according to the task allocation plan, the people and robots in the first two workstations of scheme 1 are in parallel cooperation, whereas those in the other workstations collaborate.In contrast to scheme 1, the people and robots in the first three workstations of schemes 2 and 3 are all in parallel cooperation, and those in the last four workstations collaborate.In addition, in the three schemes, the people and robots in different workstations are in the coexistence scenario, whereas those in adjacent workstations are in sequential cooperation.

Conclusion
This paper introduces different HRC scenarios into the field of DLBP for the first time.New applications of HRC technology in the intelligent remanufacturing of EOL products are also presented.This study proposes an innovative HRC-DLBP with the aim of realising four collaboration scenarios-coexistence, sequential cooperation, parallel cooperation, and collaboration.Based on this problem, an HRC disassembly line is designed.Tasks with hazardous and complex attributes are assigned to robots and humans, respectively.Tasks requiring interactive disassembly are assigned to humans and robots for completion.The cooperation mode between the two is determined according to the final assignment plan.Subsequently, an MIP model is formulated for the proposed HRC-DLBP to optimise the number of workstations, smoothing index, and various disassembly costs in the operation of the disassembly line.Two case tests demonstrated that the established MIP model could solve small-scale cases over a short time to obtain an optimal single-objective solution.
Furthermore, to solve this complex combinatorial optimisation problem efficiently, this study develops a new hybrid algorithm, HLSGA, combining local search and genetic algorithms.First, the coding and decoding strategy is redesigned according to the characteristics of the problem.The local search strategy is improved based on the priority relationship of parts, and the local search operator is adjusted to adapt to the structure of the genetic algorithm.Then, the accuracy of the proposed HLSGA is verified through two HRC-DLBP examples.Subsequently, through three HRC-DLBP examples, the HLSGA is proved to be superior to the other five excellent algorithms.Moreover, the comparisons with the published literature results proves that the proposed HLSGA also outperformed the IWOA and HMBO in solving the TDLBP and EOL state-oriented DLBP, respectively.Finally, the HLSGA is applied to a disassembly case of a power battery module; the results also demonstrate that the optimisation and convergence capabilities of the HLSGA are better than those of the other five algorithms reported in the literature.
This study has the following limitations and drawbacks.Firstly, the demand index was ignored in formulating the multiple objectives.Failure to dismantle parts with demand (value) early in the process may damage the parts, which could impact the revenue of the resource recovery company.Secondly, this study only focuses on some disassembly constraints between humans, robots, products, and disassembly lines under the four cooperation modes without considering human risk assessment and remedial measures in case of task failure.In actual disassembly production, HRC disassembly also involves many problems, such as the disassembly of multiple mixed products and scrap status of products.The limitations mentioned above will be resolved in future studies.
Fixed cost of configuring a robot C rt Given unit time cost of robot disassembly C mt Given unit time cost of manual disassembly C h Given additional cost of dealing with hazardous task OE Given energy consumption of robot SE Given standby energy consumption of robot ϕ A large positive number t i Disassembly time for task i CT Cycle time of the disassembly line TP ij Priority relationship matrix for disassembly tasks, TP ij = [a ij ] n×n ; if a ij = 1, task i is an immediately preceding task of task j Decision variables:

Figure 3 .
Figure 3. Priority relationship diagram is converted into the priority relationship matrix.

Figure 5 .
Figure 5.The operational results of the proposed four-layer encoding and heuristic decoding method.
, respectively.It is obvious that the solutions of HLSGA are all on the 9. Single-objective optimal schemes of GUROBI and HLSGA for P25.(a) Solution 1 of GUROBI (b) Solution 2 of GUROBI (c) Solution 3 of GUROBI (d) Solution 1 of HLSGA (e) Solution 2 of HLSGA (f) Solution 3 of HLSGA.Pareto frontier, and the results of the three cases solved by the proposed HLSGA are better distributed in solution space than those of the other five algorithms.Specifically, as indicated by Figure 10 (a), the ranking of the solving capabilities of the six algorithms according to the number of Pareto solutions obtained is HLSGA (six Pareto solutions) > NSGA-II (five Pareto solutions) = GDCS (five Pareto solutions) > FPA (three Pareto solutions) > MIPSO (three Pareto solutions) > GASA (zero Pareto solutions).Figure 10 (b) similarly indicates that the proposed HLSGA is better than the other algorithms in solving P25.Note that each algorithm obtained dozens of solutions after solving P55, and 10 solutions were screened by crowding distance for comparison.As indicated by Figure 10 (c), the solving ability of the HLSGA is better than those of the other algorithms.
* CM: Collaboration mode.* * Interactive tasks are assigned to robots and workers to complete cooperatively, and common tasks are assigned according to heuristic rules.

Table 2 .
Results of GUROBI and the HLSGA for solving P10 and P25 of the HRI-DLBP.

Table 3 .
Solutions obtained by HLSGA and IWOA for solving P34 of TDLBP.

Table 5 .
Solutions obtained by different algorithms for solving Tesla power battery module.