An optimization method for the balancing problem of mixed flow U-shaped disassembly line considering operator differences

Abstract

An optimization method for the balance problem of a mixed-flow U-shaped disassembly line considering operator differences relates to the technical field of facility layout. The method comprises the following steps: establishing a mathematical model, initializing algorithm parameters, encoding and generating an initial population and its corresponding external archive, performing K-mean clustering on the external archive and then starting a Markov chain, selecting any individual in the current population, judging the state of the individual, performing a Q-learning operation, determining whether to generate a new population for subsequent steps based on the maximum Markov chain length, re-clustering the external archive, performing an annealing operation, judging whether conditions are met, and outputting a final solution if conditions are met. The present invention combines a simulated annealing algorithm with a machine learning algorithm, so that the algorithm can select the action that best suits the individual's characteristics according to the individual's state characteristics to perturb the individual, improve the individual's objective function value, increase the algorithm's computational efficiency, and reasonably divide task attributes to be more in line with actual conditions.

Description

Technical Field

[0001] The present invention relates to the technical field of facility layout, and in particular to an optimization method for a mixed flow U-shaped disassembly line balancing problem taking into account operator differences. Background Art

[0002] The balancing problem of mixed-flow assembly lines refers to the disassembly of multiple products on the same assembly line. However, in the actual operation of a human-robot collaborative assembly line, the specialized tasks of disassembling different products may require different types of robots, necessitating the consideration of introducing multiple types of robots into the assembly line. Furthermore, a single assembly line also has different levels of operators. Due to their skill level, entry-level operators are only able to complete common tasks and entry-level complex tasks. Senior operators possess higher operational skills and can complete common tasks, entry-level complex tasks, and advanced complex tasks. Therefore, the problem of optimizing the balancing of human-robot collaborative assembly lines involving multiple levels of workers and multiple types of robot operators needs to be addressed.

[0003] Incomplete disassembly refers to the situation where disassembly tasks are assigned without requiring the entire task to be disassembled, but all complex and specialized tasks, including their immediate predecessors, must be performed. Completing each task generates a certain profit for the enterprise, while also incurring certain production costs. Previous research on the profit of incomplete disassembly often prioritizes the completion of all tasks whose benefits exceed their costs. However, since performing more valuable tasks (tasks where the disassembly benefits exceed the costs) may require the opening of additional workstations or the purchase of more robots, which increases fixed investment, using a long-term profit objective function that includes fixed costs can effectively balance this problem. When the disassembly line's operating time is short, the benefits of performing more valuable tasks are less than the fixed costs of operating workstations or purchasing robots, so these valuable tasks are not performed. However, if the enterprise's operating time is long, the benefits of performing more valuable tasks exceed the fixed costs of operating workstations or purchasing robots, so more fixed costs can be invested to achieve greater total profit. Therefore, this problem has more practical research value.

[0004] Current research has not considered the fact that different products require specialized tasks performed by different robot models, nor has it considered the simultaneous integration of multiple levels of workers and multiple types of robots in human-robot collaborative assembly lines. Furthermore, previous research on profit-oriented assembly line balance has not factored changes in the number of workstations and operators caused by incomplete assembly into the long-term profitability of assembly line operations, failing to fully reflect actual production conditions. Summary of the Invention

[0005] In view of this, the main purpose of the present invention is to propose an optimization method for the balancing problem of a mixed flow U-shaped disassembly line taking into account operator differences, which can efficiently optimize a multi-stage disassembly line including mixed flow operators.

[0006] The technical solution of the present invention is a method for optimizing the balancing problem of a mixed flow U-shaped disassembly line taking into account differences among operators, comprising the following steps:

[0007] Step S1: With the goal of minimizing the number of workstations, minimizing the idle time equilibrium index, and maximizing long-term profit, a mathematical model of the U-shaped disassembly line balance problem considering the differences of mixed-flow operators is established;

[0008] Step S2: Initialize the algorithm parameters, generate the initial population and its corresponding external archives by encoding, perform K-mean clustering on the external archives, and start the Markov chain;

[0009] Step S3: Select any individual current in the current population. Based on the characteristics of the individual in the external archive, compare the individual current with a random external individual in the same cluster in the external archive, determine the state of the individual current, and execute the Q-learning algorithm based on the state of the individual current. Then update the Q table and use the improved Metropolis criterion until all individuals in the current population have executed the Q-learning algorithm.

[0010] Step S4: Determine whether the chain length is greater than the maximum Markov chain length. If not, return to step S3. If greater, update the external archive and generate a new population using the Pareto method and crowding distance strategy.

[0011] Step S5: Re-cluster the external files, perform annealing operation, and determine whether the current temperature is greater than the termination temperature. If not, return to the operation step of starting the Markov chain in step S2. If it is greater, end the algorithm and output the final solution.

[0012] The technical effects of the present invention are:

[0013] 1. By combining the simulated annealing algorithm with machine learning algorithms (K-means and Q-learning), the present invention enables the algorithm to select the action that best suits the individual's characteristics to perturb the individual based on the individual's state characteristics, including task sequence, operator sequence, and number of unexecuted tasks, thereby improving the individual's objective function value. Compared with the traditional incomplete mixed-flow disassembly line algorithm, which performs perturbations on each individual without identifying its characteristics, the algorithm's computational efficiency is increased, and multiple types of robots and multiple levels of workers are simultaneously taken into account in the disassembly line. The task attributes are also reasonably divided, which is more in line with the actual situation.

[0014] 2. For the profit-oriented disassembly line balance problem, previous studies did not take the changes in the number of workstations and operators caused by incomplete disassembly into account in the long-term profit of the disassembly line operation, and could not fully reflect the actual production situation. The objective function setting of the present invention has more practical application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments.

[0016] Figure 1 It is the overall flow chart of the present invention;

[0017] Figure 2 Schematic diagram of the task sequence in the present invention;

[0018] Figure 3 (a) is a schematic diagram of the merge priority relationship matrix of the merge task;

[0019] Figure 3 (b) is a schematic diagram of the two-stage encoding process;

[0020] Figure 4 Schematic diagram of the state division of the Q-learning algorithm in the present invention;

[0021] Figure 5 Schematic diagram of the operation of the mutation action of the Q-learning algorithm in the present invention;

[0022] Figure 6 Schematic diagram of the crossover operation of the Q-learning algorithm in the present invention;

[0023] Figure 7 (a) is a schematic diagram of the operation of changing whether a task executes an action in the Q-learning algorithm of the present invention;

[0024] Figure 7 (b) is a schematic diagram of the operation of changing the operator sequence action of the Q-learning algorithm in the present invention;

[0025] Figure 8 A priority relationship diagram for a large-scale instance of the present invention;

[0026] Figure 9 Graph showing the results of parameter orthogonal experiments on a large-scale example of the present invention. DETAILED DESCRIPTION

[0027] The present invention will be further described in detail below with reference to the embodiments and accompanying drawings.

[0028] To make the purpose, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. The described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0029] Example:

[0030] The specific process of the mixed flow U-shaped disassembly line balance problem optimization method proposed in the present invention considering operator differences is as follows: Figure 1 As shown, the following steps are included:

[0031] Step S1: With the goal of minimizing the number of workstations, minimizing the idle time balance index, and maximizing long-term profits, a mathematical model for the U-shaped disassembly line balancing problem considering the differences among mixed-flow operators is established. In this invention, in order to enable accurate quantification, it is assumed that the scenario meets the following conditions:

[0032] 1. The cycle time and the information of each disassembly task are known;

[0033] 2. Ignoring the turnaround time of operators in the workstation to dismantle the import and export side tasks;

[0034] 3. Unlimited supply of disassembled products during the operation of the disassembly line;

[0035] 4. There is no unexpected interruption during the operation of the disassembly line.

[0036] The objective function based on this is:

[0037] F=min[ f 1 ,f 2 ,-f 3] (1)

[0038] In formula (1), f1 is a function related to the number of workstations, as shown in formula (2):

[0039]

[0040] f2 is the idle time related function, as shown in formula (3):

[0041]

[0042] Formula (1) represents the three objective functions, namely, minimizing the number of workstations, minimizing the idle time balance index, and maximizing long-term profits. Formula (2) represents the objective function of minimizing the number of workstations. Reducing the number of workstations can reduce the length of the disassembly line operation and the floor space occupied. Formula (3) represents the objective function of minimizing the idle time balance index of the workstation, which aims to balance the load of the operators at each workstation and prevent excessive differences in the work intensity of the operators at each workstation due to load imbalance.

[0043] f3 is the long-term profit correlation function, as shown in formula (4):

[0044]

[0045] In formula (4), C d Represents the cost of each product, as shown in formula (5):

[0046]

[0047] Equation (4) represents the objective function for maximizing the long-term operating profit of the disassembly line. The cost includes fixed costs and operating costs. Fixed costs include the workstation startup cost and the robot purchase cost. The operating cost is expressed as Equation (5), which includes the operating costs of different operators performing different tasks and the standby costs of different operators during the idle time of each workstation. The disassembly line profit is expressed as the difference between the revenue of the disassembly line performing tasks and the operating cost, multiplied by the number of products disassembled per day, multiplied by the number of operating days, minus the fixed cost. The purchase cost of the robot is high, but the subsequent maintenance and operation costs are relatively low; the worker operator has no purchase cost, but the subsequent production cost is high. Therefore, the objective function aims to solve the corresponding disassembly plan that maximizes the profit based on different production times.

[0048] The constraints that the above objective function needs to satisfy are:

[0049]

[0050] Equations (6) and (7) represent incomplete disassembly constraints. A task can be executed or not, and can be executed at most once, but the predecessor task must be executed when the successor task is executed.

[0051]

[0052] Equations (8) and (9) represent the constraints on the number of operators in a workstation. When a task is assigned to a workstation, the workstation must be turned on, and there cannot be more than one operator in the workstation.

[0053]

[0054] Formula (10) represents the workstation takt time constraint. The total working time of each workstation is expressed as the sum of the time of all tasks assigned to the workstation and the time of the operator at the workstation to perform the task, and its value should not be greater than the takt time.

[0055]

[0056] Equations (11) to (14) represent task attribute constraints. Tasks with attributes must be executed. At the same time, primary complex tasks must be assigned to worker workstations, advanced complex tasks must be assigned to advanced worker workstations, Class A special tasks must be assigned to Type A robot workstations, and Class B special tasks must be assigned to Type B workstations. Considering equipment maintenance and operating costs, Class A special tasks are defined as complex tasks that need to be handled by robots, and Class B special tasks are defined as simple tasks that need to be handled by robots.

[0057]

[0058]

[0059] Equations (15) and (16) represent precedence constraints. When both the predecessor and successor tasks are on the import side, the successor task cannot be executed at the workstation before the predecessor task. When the predecessor task is assigned to the export side, the successor task must be assigned to the export side and cannot be executed at the workstation after the predecessor task. Considering the incomplete disassembly constraint, the successor task can be omitted when the predecessor task is executed.

[0060]

[0061] Equations (17) and (18) represent the workstation startup constraints. When a workstation is started, there must be tasks assigned to it, and the workstations must be started in sequence.

[0062] In formulas (1)-(18), W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i and j are task numbers respectively; n is the workstation number; m is the judgment variable for the import and export ends of the disassembly line, m=1 represents the import end, and m=2 represents the export end; k is the operator variable, k=1 represents the operator is a junior worker, k=2 represents the operator is a senior worker, k=3 represents the operator is a Class A robot, and k=4 represents the operator is a Class B robot; M is a positive number greater than the total number of tasks; K1 is the known set of primary complex tasks; K2 is the known set of advanced complex tasks; S1 is the known set of Class A special tasks; S2 is the known set of Class B special tasks; TP ijis the priority relationship matrix. If task i is the immediate predecessor of task j, then TP ij =1; CT is the takt time; is the disassembly time of task i when the operator is k; is the cost of task i when operator is k; r k is the standby cost when the number of operators at the workstation is l; P i is the disassembly benefit of task i; R k is the single acquisition cost of type k robot; V is the unit workstation startup cost; C d is the cost of disassembling each EOL (End of Life) product; T is the number of days the disassembly line operates; G is the total number of parts disassembled by the disassembly line per day; is a binary variable, Indicates that task i is assigned to port m of the nth workstation whose operator is k, otherwise x j mkn is a binary variable, x j mkn =1 means task j is assigned to port m of the nth workstation whose operator is k, otherwise x j mkn =0;z kn is a binary variable, z kn =1 means the nth workstation is open and the operator is k, otherwise z kn =0;

[0063] Step S2: Initialize the algorithm parameters, generate the initial population and its corresponding external archives by encoding, perform K-mean clustering on the external archives, and start the Markov chain;

[0064] The algorithm used in this paper is SA-KQL (Simulated annealing algorithm with K-means and Q-Learning). Its basic idea is to use K-means clustering method to assist in individual state division under the framework of simulated annealing, and adopt ε-greedy strategy to select the action to be executed according to the Q table.

[0065] Among them, the encoding process takes into account that the tasks in the U-shaped disassembly line need to be assigned to both the import and export sides, so the task sequence is usually encoded in a mixed way of positive and negative numbers. Positive numbers indicate that the task is assigned to the import side of the workstation, and negative numbers indicate that the task is assigned to the export side of the workstation. In terms of incomplete disassembly, previous studies usually adopted a two-layer encoding method, introducing a second layer of task execution sequence to indicate whether the task is executed. Figure 2As shown, the present invention introduces complex domain task sequence encoding, improving the original two-layer encoding to a single layer. The task execution is determined by judging whether the imaginary part of the task sequence is 0. In the figure, white tasks (i.e., imaginary task numbers) indicate unexecuted tasks, while colored tasks (i.e., real task numbers) indicate executed tasks.

[0066] Generally speaking, DLBP (Disassembly Line Balancing Problem) problems need to satisfy task priority constraints. In order to facilitate computer calculation, it is necessary to use Figure 3 The merge priority relationship matrix shown in (a) is used to represent the relationship between merged tasks. Based on this, the encoding process of the present invention adopts a two-stage encoding method. The first stage generates the initial task sequence. The specific operation process is as follows: tasks with no immediate predecessor or immediate successor are randomly selected for assignment. If the task has no immediate predecessor, the task is assigned to the entrance of the U-shaped disassembly line with a positive task number; if the task has no immediate successor, the task is assigned to the exit of the U-shaped disassembly line with a negative task number. The TP matrix is then updated, and the above operation is repeated until all tasks are assigned. The second stage selects unexecuted tasks and generates an operator sequence. The specific operation process is as follows: first, a set of tasks E that can be unexecuted is determined. Any tasks not exceeding E (including 0 tasks) are randomly selected. The selected tasks and their immediate successors are multiplied by the imaginary unit i in the task sequence to represent them as unexecuted tasks. Finally, considering the different operators at each workstation, an operator sequence consisting of random numbers [1, 4] is generated. The nth number in the sequence indicates that the operator of the nth workstation is k, and the operator represented by k is the same as the operator represented by k in Section 1.2. The specific operation process is as follows Figure 3 As shown in (b).

[0067] For the subsequent decoding process in the present invention, since different operators have different working task times, the number of tasks that can be assigned to the workstation cannot be determined when the workstation has not determined the operator; since different tasks need to be assigned to different operators to perform the work, the operator in the workstation cannot be determined when the task attributes in the workstation are not determined. In order to resolve this contradiction, the present invention proposes a decoding method for a variable operator sequence, the main idea of which is: when the task in the workstation conflicts with the operator k in the operator sequence, if there is an operator k_new that can execute all tasks within the takt time, then the operator sequence is changed and the operator of the workstation is changed to k_new; if all operators cannot fully execute the tasks in the workstation, that is, the workstation has tasks with multiple attributes, then operator k is still used as the operator of the current workstation. The specific operation is as follows: after opening the ws-th workstation, first determine the set of executable tasks I that meet the takt time constraints under the execution time of each operator k, check the set of tasks I performed by the operator ZL(ws) in the operator sequence ZL ZL(ws) Is there a task conflict? If not, update the task index i=i+size(I ZL(ws) ); If there is a task conflict, check the remaining I k Check if there is a task attribute conflict. If so, change the operator ZL(ws) to k_new and update the task index. If k_new does not exist, use the original operator to execute the task before the task with the task attribute conflict. In special cases, if the first task of the workstation is a task with a task attribute conflict and the original operator cannot be used, change the operator to enable the workstation to start smoothly. The specific operation is shown in the pseudo code in Table 1:

[0068] Table 1 Decoding operation pseudo code

[0069]

[0070]

[0071] In the solution process, especially for large-scale problems, the exact solver's solution time will increase exponentially, which makes it impossible to solve large-scale instances in a limited time. Therefore, in order to obtain better results, it is necessary to determine the SA-KQL algorithm parameters, including the initial population size pop_size, the initial temperature T initial , cooling coefficient q, maximum Markov chain length maxL, discount factor γ, learning rate α, the determination method mainly includes the following steps:

[0072] Step 1) Determine the value range of each algorithm parameter and perform multivariate variance analysis on each algorithm parameter;

[0073] Step 2) Select parameters with a p-value less than 0.05 in the multivariate analysis of variance to conduct an orthogonal experiment with hypervolume as the indicator, and use the maximum value of each parameter in the orthogonal experiment as the value of the parameter in the algorithm;

[0074] Step 3) The middle value of the range of parameters with a p-value not less than 0.05 in the multivariate variance analysis is used as its value in the algorithm.

[0075] Step S3: Select any individual current in the current population. Based on the characteristics of the individual in the external archive, compare the individual current with a random external individual in the same cluster in the external archive, determine the state of the individual current, and execute the Q-learning algorithm based on the state of the individual current. Then update the Q table and use the improved Metropolis criterion until all individuals in the current population have executed the Q-learning algorithm.

[0076] The three basic related factors of the Q-learning process are: state, action and reward. State information represents the changes caused by the agent's perception of environmental information and its own actions. In previous studies, the division design of the state space is usually only associated with the fitness value. However, since the basic problem characteristics of the combinatorial optimization problem are not unique, a state space that takes into account both the characteristics of the feasible solution sequence and the fitness value is proposed. In terms of fitness value, since this problem is a multi-objective optimization problem, it is impossible to simply explain the fitness value state of the individual by comparing the size of the target value. By observing previous studies on DLBP problems, it is found that the solution set of DLBP is usually divided into multiple clusters according to the optimization direction of each target value. Therefore, this study uses the K-means clustering method to divide the external solution set into K clusters, and then divides the individuals into the Kth cluster according to the original division method. n cluster, then the target value state of the individual is K n .

[0077] In order to eliminate the influence of different orders of magnitude of target values in calculating the Euclidean distance to the cluster center, the target data are standardized by using the Z score, i.e., the standardized score. The Z score calculation formula is shown in formula (24).

[0078]

[0079] Where, f a represents the ath objective function value of an individual, μ represents the mean of the objective function data set, and σ represents the standard deviation of the data set.

[0080] In terms of sequence resolution, the sequence of problems in this invention has three main characteristics: the order of executed tasks, the sequence of operators, and the tasks not executed. These three characteristics will all affect the final objective function value, so the sequence resolution should consider all three aspects simultaneously. On this basis, the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive is determined in terms of the three individual characteristics of the order of executed tasks, the sequence of operators, and the tasks not executed. After combining the similarities in each sequence characteristic, combined with the fitness value classification of the K clusters in the K-mean clustering of the individual current, the state of any individual current in the current population is obtained.

[0081] Among them, the similarity judgment of individual characteristics of task execution sequence is shown in formula (19):

[0082]

[0083] In formula (19), D hm_XL is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of the order of executing tasks; W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i is the task number, i = 1 when the i-th task or operator is the same, otherwise i = 0; XL_d extrnal [i] is the task execution sequence of a random external individual in the same cluster in the external file; XL_d current [i] is the task execution sequence of any individual current in the current population;

[0084] When D hm _XL≥size(XL_d extrnal ) / 2, that is, when the similarity of the task execution sequence of individual current and external individual external is greater than half, the two are in a similar state, otherwise they are in a dissimilar state;

[0085] The similarity judgment of individual characteristics of operator sequences is shown in formula (20):

[0086]

[0087] In formula (20), D hm _ZL is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of operator sequence; W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i is the task number, i = 1 when the i-th task or operator is the same, otherwise i = 0; ZL extrnal [i] is the operator sequence of a random external individual external in the same cluster in the external file; ZL current [i] is the operator sequence of any individual current in the current population;

[0088] When D hm _ZL≥size(ZL_d extrnal ) / 2, that is, when the operator sequence similarity between individual current and external individual external is greater than half, the two are in a similar state, otherwise they are in a dissimilar state;

[0089] The similarity judgment of individual characteristics of unexecuted tasks is shown in formula (21):

[0090] ΔJ=|J external -J current | (21)

[0091] △J is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of unperformed tasks; J extrnal is the number of unexecuted tasks of a random external individual external in the same cluster in the external archive; J current is the number of unexecuted tasks of any individual current in the current population;

[0092] When △J≥J extrnal / 2, that is, when the number of unexecuted tasks in the individual current is more than half of the number of unexecuted tasks in the external individual external, the two are in a similar state, otherwise they are in a dissimilar state.

[0093] The fitness value and sequence similarity can be used to design a specific state space. The fitness value is classified into K states according to the K clusters. The three sequence similarities have two states - similar and dissimilar. In this case, the three sequences can be combined to obtain a total of 8 states. The combined 8 states are then multiplied by the K states obtained from the K cluster classification, resulting in a total of 8×K states. It can be seen that these 8×K states are all the states that the Q-learning algorithm can have, and the sequence of any individual current will belong to one of these 8×K states. The specific division process of the state is shown in the figure below. Figure 4 shown.

[0094] In terms of Q-learning action design, the actions (disturbance operations) in traditional DLBP optimization algorithms usually only consider perturbations to the task sequence. However, in the problem addressed in this invention, the objective function value is not only related to the task sequence, but also requires the design of perturbations to the non-executed tasks (the imaginary part of the task sequence) and the operator sequence.

[0095] Therefore, four actions are designed: mutation, crossover, changing whether the task is executed, and changing the operator sequence. Among them, the mutation action is designed as follows:

[0096] A mutation operation that can change the execution of the import and export side of the task. First, randomly select whether the mutation task needs to change the import and export side, then determine the mutation section of the task on that side, and finally randomly insert the task into the mutation section. The specific operation is as follows Figure 5As shown in the example, a task (such as Task 3 in the figure) is randomly selected with equal probability. A random number, rand, between 0 and 1 is then generated. If rand is less than 0.5, its import and export terminals are changed (i.e., their sign is positive or negative). If rand is greater than 0.5, the import and export terminals are not changed. The mutable range is determined by the immediately preceding and immediately succeeding tasks. After the mutated task, the task must still satisfy the precedence constraint.

[0097] Crossover action is designed to:

[0098] The crossover operation is a process in which the current individual learns from the selected external individual. The influence of the import and export sides also needs to be considered. Randomly select the import or export side of the current individual's task sequence, rearrange it according to the actual execution order of the external individual, and return it to the original sequence. The specific operation is as follows: Figure 6 As shown in the figure - first randomly select a positive sequence representing the import end, or a negative sequence representing the export end, and then randomly select a fragment from the randomly selected sequence and rearrange it according to the order in the external sequence.

[0099] The action to change whether a task is executed is designed as follows:

[0100] In order to consider the characteristics of the unexecuted task sequence, i.e., the impact of the imaginary part of the task in the sequence, a variable task execution action is designed. A task is randomly selected from the set of non-attribute tasks and its immediate predecessor tasks (i.e., tasks that can be unexecuted). If the task is executed in the original sequence, the task and its immediate predecessor tasks are multiplied by the imaginary unit i to convert them into unexecuted tasks. If the task is not executed, the imaginary units of the task and its immediate predecessor tasks are removed to convert them into executable tasks. The specific operation is as follows: Figure 7 As shown in (a).

[0101] Change the operator sequence action design to:

[0102] Considering the impact of the operator sequence on the individual fitness value, a mutation operation is designed for the operator sequence. Since the actual number of workstations in the decoding of this problem is less than the length of the operator sequence, in order to prevent invalid mutation, according to the value of the objective function f1_external of the number of workstations of the external archive individual external, one of the first f1_external operators in the current individual operator sequence is randomly selected for perturbation and changed to another operator. The specific operator is as follows Figure 7 (b) shown.

[0103] Furthermore, in order to increase the disturbance of the action on the original individual, the four actions of mutation, crossover, changing whether the task is executed, and changing the operator sequence are combined and executed simultaneously in pairs, resulting in a total of eight actions: mutation, crossover, changing whether the task is executed, changing the operator sequence, mutation and changing whether the task is executed at the same time, mutation and changing the operator sequence at the same time, crossover and changing whether the task is executed at the same time, and crossover and changing the operator sequence at the same time.

[0104] In Q-learning, the cumulative reward of the action-state pair is used to maximize the benefit during the iteration process. The Q(state, action) value represents the Q value of the state corresponding to the action and is stored in the Q-table. The initial value of the Q-table is 0. The essence of QL is to continuously learn and update the Q-table during the iteration process. Its update formula is shown in Equation (25):

[0105]

[0106] Where state represents the current state, action represents the action selected in the current state, state' represents the next state, and action' represents the action selected in the next state. r represents the reward value, α represents the learning rate, whose value is positively correlated with the historical learning effect, and γ is the discount factor that represents the impact of the reward on the new state.

[0107] Due to the multi-objective nature of this problem, the reward and penalty value r is designed based on the Pareto method. When the objective function value after executing the action dominates the original individual, r = 2; when the two do not dominate each other, r = 1; otherwise, r = -1. This algorithm adopts a dynamic ε-greedy strategy. The ε value increases with the algorithm iteration, indicating that as the Q table continues to learn, the QL action selection should trust the Q table more. The dynamic ε-greedy strategy selects the action to be executed, which mainly includes the following steps:

[0108] Generate a random number rand, and when rand>ε, randomly select an action; otherwise, select the action with the maximum value in the Q table under this state. The calculation formula of ε is shown in formula (22):

[0109]

[0110] In formula (22), T is the current temperature value, T initial is the initial temperature, T final is the termination temperature.

[0111] In addition, since the present invention involves a multi-objective problem, the algorithm uses the improved Metropolis criterion to determine whether to accept a new solution when updating the population pop. The probability of accepting a new solution is shown in formula (23):

[0112]

[0113] In formula (19), P is the acceptance probability, represents the ath objective function value of the new solution, represents the ath objective function value of the current solution, and T is the current temperature value.

[0114] Step S4: Determine whether the chain length is greater than the maximum Markov chain length. If not, return to step S3. If greater, update the external archive and generate a new population using the Pareto method and crowding distance strategy.

[0115] Step S5: Re-cluster the external files, perform annealing operation, and determine whether the current temperature is greater than the termination temperature. If not, return to the operation step of starting the Markov chain in step S2. If it is greater, end the algorithm and output the final solution.

[0116] The following is a demonstration of the actual operation of this embodiment:

[0117] Program testing environment:

[0118] The simulation computing environment for the numerical experiments in this chapter is Windows 11 system, AMD Ryzen 76800H CPU, 3.20GHz, 16GB memory, and the Gurobi10.0.2 exact solver is used to run the solution.

[0119] Large-scale instance applications:

[0120] The U-shaped disassembly line is suitable for the disassembly of small and medium-sized electrical and electronic products. Therefore, this problem model is applied to a disassembly line that mixes monitors and laptops. In order to quantify this problem, the monitor is divided into 27 disassembly tasks and the laptop is divided into 29 disassembly tasks. The disassembly priority relationship is as follows: Figure 8 The specific disassembly information is shown in Table 2, where the order of disassembly time and cost is: junior worker, senior worker, type A robot, and type B robot. The cycle time is CT = 45s.

[0121] Table 2 Large-scale instance disassembly information

[0122]

[0123]

[0124] As the problem size increases, the exact solver's solution time increases exponentially, and the exact solver is no longer able to solve large-scale instances within a limited time. Therefore, in order to obtain better results, the SA-KQL algorithm parameters need to be determined. The hypervolume index (HV) is used as an evaluation index for the algorithm results. The higher the HV value, the better the solution convergence and solution quality. The calculation formula of the HV value is shown in Equation (26):

[0125]

[0126] In the formula, V(F(Xi),REF*) represents the hypervolume formed by the solution set and the reference point REF*. The reference point is determined by the worst value of each sub-objective in the different problems. The worst value of each objective is selected from a randomly generated initial population and rounded up to the nearest integer as the reference point for this example. In this example, the reference point [15, 6000, -100000] is selected for calculation.

[0127] Multivariate Analysis of Variance (MANOVA) can simultaneously consider the relationship between multiple variables and the interaction between each variable, thereby providing multi-dimensional statistical data analysis, which is of great significance in statistics. Parameters in SA-KQL: pop_size = {140, 160, 180}, maxL = {20, 22, 25}, T initial ={22, 25, 28}, q = {0.92, 0.93, 0.95}, α = {0.6, 0.7, 0.8}, γ = {0.2, 0.3, 0.4}; the termination temperature is fixed at 1, and 3 are shared. 6 Parameter combinations, calculate HV value 10 times for each combination and take the average value, a total of 10*3 6 =7290 times. The calculation results were subjected to multivariate variance analysis, and the analysis results are shown in Table 3:

[0128] Table 3 Results of multivariate analysis of variance

[0129]

[0130]

[0131] In MANOVA, when the p-value is less than 0.05, it is considered that the parameter has an impact on the algorithm performance. By observing the data in the table, the parameter T initialThe p-values of , q, maxL, and maxL*α are all less than 0.05, proving that they have an impact on the algorithm performance. Since the solution space for large-scale problems is large, the population size setting value is much smaller than the solution space size, so it has no effect on the solution performance; the learning rate α single factor p-value is greater than 0.05, but the p-value of the interaction with maxL is less than 0.05. This is because α has a nonlinear effect on the algorithm performance and is not easy to be discovered in a single factor analysis. The effect becomes significant when the interaction is considered. In order to select the best parameter combination, an orthogonal experiment was conducted on the above four algorithm parameters, and the evaluation index was HV. The experimental results are as follows: Figure 9 shown.

[0132] according to Figure 9 As shown in the figure, the parameter values with better solution effects in the orthogonal test results are selected as algorithm parameters. For parameters with p-value greater than 0.05, the middle value within the parameter range is selected as the algorithm parameter value.

[0133] To better demonstrate the superiority of the proposed algorithm, we compared it with several leading algorithms in current human-computer collaboration research: the multi-objective enhanced differential evolution algorithm (MEDE), the multi-objective shuffled frog leaping algorithm (SFLA), and the multi-objective modified teaching and learning optimization algorithm (MTLO). The results are shown in Table 4. The data in this table represent the 10 final solutions selected by sorting each algorithm using crowding distance. The bolded data represent the optimal values obtained for the objective, while the underlined data indicate that the solution was dominated by the solutions of other algorithms.

[0134] Table 4 Comparison of algorithm solution results

[0135]

[0136]

[0137] Comparative results show that, from a single-objective perspective, the optimal value of the idle time balance indicator remains 0 in large-scale instances, indicating that combining incomplete disassembly with a multi-operator U-shaped disassembly line effectively improves disassembly line flexibility. Furthermore, for f1, only SA-KQL and MTLO achieved an optimal solution of 6, for f3, only SA-KQL and SFLA achieved an optimal solution of 286,930, and for f2, only SA-KQL achieved an optimal solution of 0, demonstrating SA-KQL's strong deep search capabilities. From a multi-objective perspective, 8, 9, and 10 final solutions of the three compared algorithms, respectively, were dominated by other algorithms, while none of the SA-KQL algorithms' solutions were dominated by other algorithms. This indicates that the non-inferior solutions obtained by SA-KQL are closer to the true Pareto frontier. In summary, SA-KQL exhibits superior solution performance compared to classical algorithms.

[0138] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the embodiments of the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. An optimization method for the balancing problem of a mixed flow U-shaped disassembly line considering operator differences, characterized in that: The following steps are involved: Step S1: With the goal of minimizing the number of workstations, minimizing the idle time balance index, and maximizing long-term profits, a mathematical model for the U-shaped disassembly line balance problem considering the differences among mixed-flow operators is established. Its objective function is: F=min[f1,f2,-f3] (1) In formula (1), f1 is a function related to the number of workstations, as shown in formula (2): f2 is the idle time related function, as shown in formula (3): f3 is the long-term profit correlation function, as shown in formula (4): In formula (4), C d Represents the cost of each product, as shown in formula (5): The constraints that the objective function needs to satisfy are: In formulas (1)-(18), W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i and j are task numbers respectively; n is the workstation number; m is the judgment variable for the import and export ends of the disassembly line, m=1 represents the import end, and m=2 represents the export end; k is the operator variable, k=1 represents the operator is a junior worker, k=2 represents the operator is a senior worker, k=3 represents the operator is a Class A robot, and k=4 represents the operator is a Class B robot; M is a positive number greater than the total number of tasks; K1 is the known set of primary complex tasks; K2 is the known set of advanced complex tasks; S1 is the known set of Class A special tasks; S2 is the known set of Class B special tasks; TP ij is the priority relationship matrix. If task i is the immediate predecessor of task j, then TP ij =1; CT is the takt time; is the disassembly time of task i when the operator is k; is the cost of task i when operator is k; r k is the standby cost when the number of operators at the workstation is l; P i is the disassembly benefit of task i; R k is the single acquisition cost of type k robot; V is the unit workstation startup cost; C d is the cost of disassembling each EOL product; T is the number of days the disassembly line operates; G is the total number of parts disassembled by the disassembly line per day; is a binary variable, Indicates that task i is assigned to port m of the nth workstation whose operator is k, otherwise x j mkn is a binary variable, x j mkn =1 means task j is assigned to port m of the nth workstation whose operator is k, otherwise x j mkn =0;z kn is a binary variable, z kn =1 means the nth workstation is open and the operator is k, otherwise z kn =0; Step S2: Initialize the algorithm parameters, generate the initial population and its corresponding external archives by encoding, perform K-mean clustering on the external archives, and start the Markov chain; Step S3: Select any individual current in the current population. Based on the characteristics of the individual in the external archive, compare the individual current with a random external individual in the same cluster in the external archive, determine the state of the individual current, and perform Q-learning operations based on the state of the individual current. Then update the Q table and use the improved Metropolis criterion until all individuals in the current population have undergone Q-learning operations. Step S4: Determine whether the chain length is greater than the maximum Markov chain length. If not, return to step S3. If greater, update the external archive and generate a new population using the Pareto method and crowding distance strategy. Step S5: Re-cluster the external files, perform annealing operation, and determine whether the current temperature is greater than the termination temperature. If not, return to the operation step of starting the Markov chain in step S2. If it is greater, end the algorithm and output the final solution.

2. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The algorithm parameters in step S2 include the initial population size pop_size, the initial temperature T initial , cooling coefficient q, maximum Markov chain length maxL, discount factor γ, and learning rate α.

3. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 2, characterized in that: The method for determining the algorithm parameters in step S2 mainly includes the following steps: Step 1) Determine the value range of each algorithm parameter and perform multivariate variance analysis on each algorithm parameter; Step 2) Select parameters with a p-value less than 0.05 in the multivariate analysis of variance to conduct an orthogonal experiment with hypervolume as the indicator, and use the maximum value of each parameter in the orthogonal experiment as the value of the parameter in the algorithm; Step 3) The middle value of the range of parameters with a p-value not less than 0.05 in the multivariate variance analysis is used as its value in the algorithm.

4. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The individual features of the external archive in step S3 include the order of executed tasks, the sequence of operators, and the tasks not executed.

5. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 4, characterized in that: The method for determining the status of any individual current in the current population in step S3 is: Determine the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of the three individual characteristics of task execution sequence, operator sequence, and unexecuted tasks. After combining the similarities in each sequence characteristic, combine the fitness value classification of the K clusters in the K-mean clustering of the individual current to obtain the status of any individual current in the current population. Among them, the similarity judgment of individual characteristics of task execution sequence is shown in formula (19): In formula (19), D hm _XL is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of the order of executing tasks; W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i is the task number, i = 1 when the i-th task or operator is the same, otherwise i = 0; XL_d extrnal [i] is the task execution sequence of a random external individual in the same cluster in the external file; XL_d current [i] is the task execution sequence of any individual current in the current population; When D hm _XL≥size(XL_d extrnal ) / 2, that is, when the similarity of the task execution sequence of individual current and external individual external is greater than half, the two are in a similar state, otherwise they are in a dissimilar state; The similarity judgment of individual characteristics of operator sequences is shown in formula (20): In formula (20), D hm _ZL is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of operator sequence; W is the total number of disassembly tasks and the maximum number of workstations, where the task set of product 1 is W1 = {1, 2, 3…, w1}, and the task set of product 2 is W2 = {w1+1, w1+2, w1+3…, w1+w2}, W∈W1∪W2; i is the task number, i = 1 when the i-th task or operator is the same, otherwise i = 0; ZL extrnal [i] is the operator sequence of a random external individual external in the same cluster in the external file; ZL current [i] is the operator sequence of any individual current in the current population; When D hm _ZL≥size(ZL_d extrnal ) / 2, that is, when the operator sequence similarity between individual current and external individual external is greater than half, the two are in a similar state, otherwise they are in a dissimilar state; The similarity judgment of individual characteristics of unexecuted tasks is shown in formula (21): ΔJ=|J external -J current | (21) △J is the similarity between any individual current in the current population and a random external individual external in the same cluster in the external archive in terms of unperformed tasks; J extrnal is the number of unexecuted tasks of a random external individual external in the same cluster in the external archive; J current is the number of unexecuted tasks of any individual current in the current population; When △J≥J extrnal / 2, that is, when the number of unexecuted tasks in the individual current is more than half of the number of unexecuted tasks in the external individual external, the two are in a similar state, otherwise they are in a dissimilar state.

6. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 5, characterized in that: The number of state types of the Q-learning algorithm in step S3 is 8×K, where K is the number of fitness value classifications of K clusters in K-mean clustering.

7. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The actions performed by the Q-learning algorithm in step S3 include mutation, crossover, changing whether a task is to be executed, changing an operator sequence, mutation and changing whether a task is to be executed simultaneously, mutation and changing an operator sequence simultaneously, crossover and changing whether a task is to be executed simultaneously, and crossover and changing an operator sequence simultaneously.

8. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The Q-learning algorithm in step S3 uses a dynamic ε-greedy strategy to select the action to be executed, which mainly includes the following steps: Generate a random number rand, and when rand>ε, randomly select an action; otherwise, select the action with the maximum value in the Q table under this state. The calculation formula of ε is shown in formula (22): In formula (22), T is the current temperature value, T initial is the initial temperature, T final is the termination temperature.

9. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The probability of accepting a new solution in the improved Metropolis criterion in step S3 is shown in formula (23): In formula (19), P is the acceptance probability, represents the ath objective function value of the new solution, represents the ath objective function value of the current solution, and T is the current temperature value.

10. The method for optimizing the balancing problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 1, characterized in that: The decoding method is a decoding method for a variable operator sequence. The specific steps include: when the task in the workstation conflicts with the operator k in the operator sequence, if there is an operator k_new that can execute all tasks within the takt time, then the operator sequence is changed and the operator of the workstation is changed to k_new; if all operators cannot fully execute the tasks in the workstation, that is, the workstation has tasks with multiple attributes, then operator k is still used as the operator of the current workstation.

Images