Optimization method for balance problem of mixed-flow U-shaped disassembly line considering operator difference

Abstract

A mixed flow U-shaped disassembly line balance problem optimization method considering operator differences relates to the technical field of facility layout, and comprises the following steps: establishing a mathematical model, initializing algorithm parameters, encoding to generate an initial population and an external file corresponding to the initial population, performing K-mean clustering on the external file, then opening a Markov chain, selecting any individual in the current population, and selecting the selected individual in the current population; judging the state of the individual, carrying out Q-learning operation, judging whether a new population is generated or not to carry out subsequent steps according to the maximum Markov chain length, carrying out re-clustering on external archives, executing annealing operation, judging whether conditions are met or not, and outputting a final solution if the conditions are met; according to the method, the simulated annealing algorithm and the machine learning algorithm are combined, so that the algorithm can select the action most suitable for the individual characteristics to disturb the individual according to the state characteristics of the individual, the objective function value of the individual is improved, the calculation efficiency of the algorithm is improved, the task attributes are reasonably divided, and the actual situation is better met.

Description

Technical Field

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

[0002] The problem of balancing a mixed-flow disassembly line refers to the disassembly of multiple products on the same disassembly line. However, in the actual operation process of a human-machine collaborative disassembly line, the special tasks of different disassembled products may need to be assigned to different types of robots for execution, which requires considering introducing multiple types of robots into the disassembly line. At the same time, there are also workers of different levels on a disassembly line. Due to the reason of skill proficiency, junior workers can only complete ordinary tasks and primary complex tasks; senior workers have higher operation skills and can complete three types of tasks: ordinary tasks, primary complex tasks, and advanced complex tasks. Therefore, it is necessary to solve the problem of optimizing the balance of a human-machine collaborative disassembly line including multi-level workers and multi-type robot operators.

[0003] Incomplete disassembly means that when assigning disassembly tasks, it is not necessary to disassemble all tasks, but all complex tasks, special tasks, and their preceding tasks must be executed. After each task is completed, the enterprise will obtain a certain income, and at the same time, a certain production cost will be generated. In previous studies on the profit of incomplete disassembly, the goal was often to complete all tasks with income greater than cost. However, since executing more valuable tasks (tasks with disassembly income greater than cost) may open more workstations or purchase more robots, increasing additional fixed investment, using the long-term profit including fixed costs as the objective function can effectively balance this problem: when the operation time of the disassembly line is short, the income brought by executing more valuable tasks is less than the fixed cost invested in opening workstations or purchasing robots, so these valuable tasks are not executed; if the enterprise's operation time is long, the income brought by executing more valuable tasks is greater than the fixed cost invested in opening workstations or purchasing robots, then more fixed costs can be invested to obtain more total profit. Therefore, this problem has more practical research value.

[0004] In current research, it does not consider that the special tasks of different products need to be completed by different types of robots, nor does it consider multi-level workers and multi-type robots simultaneously in the human-machine collaborative disassembly line. At the same time, for the problem of balancing a disassembly line with profit orientation, in previous studies, the change in the number of workstations and the number of operators caused by incomplete disassembly was not considered in the long-term profit of the disassembly line operation, and it could not fully reflect the actual production situation. Summary of the Invention

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

[0006] The technical solution of the present invention is an optimization method for the problem of balancing a mixed-flow U-shaped disassembly line considering operator differences, including the following steps:

[0007] Step S1: Establish a mathematical model for the problem of balancing a U-shaped disassembly line considering mixed-flow operator differences with the goals of minimizing the number of workstations, minimizing the idle time balance index, and maximizing the long-term profit;

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

[0009] Step S3: Select any individual current in the current population. Based on the individual characteristics of the external archive, compare this individual current with a randomly selected external individual external in the same cluster in the external archive, judge the state of this individual current, execute the Q-learning algorithm according to the state of this individual current, and then update the Q table and according to the improved Metropolis criterion until all individuals in the current population have executed the Q-learning algorithm;

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

[0011] Step S5: Re-cluster the external archive, perform annealing operation, judge whether the current temperature is greater than the termination temperature. If it is not greater, 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 and machine learning algorithms (K-means and Q-learning), the present invention enables the algorithm to select the action most suitable for the individual characteristics according to the state characteristics of the individual, including the task sequence, operator sequence, and number of non-executed tasks, to perturb the individual and improve the objective function value of the individual. Compared with the method of the traditional incomplete mixed-flow disassembly line algorithm that perturbs each individual without identifying its characteristics, the calculation efficiency of the algorithm is increased, multiple types of robots and multiple levels of workers are simultaneously considered in the disassembly line, and the task attributes are reasonably divided, which is more in line with the actual situation.

[0014] 2. For the profit-oriented disassembly line balancing problem, the changes in the number of workstations and operators caused by incomplete disassembly were not considered in the long-term profit of the disassembly line operation in previous studies, which could not fully reflect the actual production situation. However, the objective function of the present invention is set to be more valuable in practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the embodiments will be briefly introduced below.

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

[0017] Figure 2 It is the schematic diagram of the task sequence in the present invention;

[0018] Figure 3 In (a), it is the schematic diagram of the merge precedence relationship matrix of the merged tasks;

[0019] Figure 3 In (b), it is the schematic diagram of the two-stage coding process;

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

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

[0022] Figure 6 It is the schematic diagram of the operation of the crossover action of the Q-learning algorithm in the present invention;

[0023] Figure 7 In (a), it is the schematic diagram of the operation of the Q-learning algorithm in the present invention for changing whether a task is executed;

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

[0025] Figure 8 It is the precedence graph of the large-scale instance in the present invention;

[0026] Figure 9 It is the result graph of the parameter orthogonal experiment of the large-scale instance in the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

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

[0028] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. The described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0029] Example:

[0030] The specific process of an optimization method for the mixed-model U-shaped disassembly line balance problem considering operator differences proposed in the present invention is as Figure 1 shown, including the following steps:

[0031] Step S1: With the goals of minimizing the number of workstations, minimizing the idle time balance index, and maximizing the long-term profit, establish a mathematical model for the U-shaped disassembly line balance problem considering mixed-model operator differences. In the present invention, in order to be able to accurately quantify, it is assumed that the scenario satisfies the following conditions:

[0032] 1. The takt time is known, and the information of each disassembly task is known;

[0033] 2. Ignore the turning time of the operator for the tasks on the inlet and outlet sides of the disassembly within the workstation;

[0034] 3. The supply of disassembled products is unlimited during the operation of the disassembly line;

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

[0036] Based on this, the objective function is:

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

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

[0039]

[0040] f 2 is a function related to the idle time, specifically as shown in formula (3):

[0041]

[0042] Equation (1) represents that the three objective functions are respectively to minimize the number of workstations, minimize the idle time balance index, and maximize the long-term profit; Equation (2) represents the objective function of minimizing the number of workstations. Reducing the number of workstations can reduce the length and floor area of the disassembly line operation; Equation (3) represents the objective function of minimizing the idle time balance index of the workstations, aiming to balance the load of the operators at each workstation and prevent the workload difference of the operators at each workstation from being too large due to unbalanced load.

[0043] f 3 is a function related to long-term profit, specifically as shown in Equation (4):

[0044]

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

[0046]

[0047] Equation (4) represents the objective function of maximizing the long-term operating profit of the disassembly line. The cost includes fixed cost and operating cost. The fixed cost includes the cost of starting the workstation and the cost of purchasing robots; the operating cost is expressed as Equation (5), which includes the operating cost of different operators performing different tasks and the standby cost of the idle time of different operators at each workstation. The profit of the disassembly line 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 robots is relatively high, but the subsequent maintenance and operating costs are relatively low; there is no purchase cost for worker operators, but the subsequent production cost is relatively high. Therefore, this objective function aims to solve the corresponding disassembly plan that maximizes the revenue according to different production times.

[0048] The constraint conditions that the above objective functions need to satisfy are:

[0049]

[0050] Equations (6) and (7) represent the incomplete disassembly constraint. A task can be selected to be executed or not, and can be executed at most once, but when the subsequent task is executed, the preceding task must be executed.

[0051]

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

[0053]

[0054] Equation (10) represents the cycle time constraint of the workstation. The total working time of each workstation is expressed as the sum of the times of all tasks assigned to that workstation as executed by the operator at that workstation, and its value should not be greater than the cycle time.

[0055]

[0056] Equations (11) to (14) represent the 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, type A special tasks must be assigned to type A robot workstations, and type B special tasks must be assigned to type B workstations. Among them, considering equipment maintenance and operating costs, type A special tasks are defined as complex tasks that need to be processed by robots, and type B special tasks are defined as simple tasks that need to be processed by robots.

[0057]

[0058]

[0059] Equations (15) and (16) represent the precedence relationship constraints. When both the predecessor task and the successor task are on the import side, the successor task cannot be executed at a 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 a workstation after the predecessor task. At the same time, considering the incomplete disassembly constraint, when the predecessor task is executed, the successor task may not be executed.

[0060]

[0061] Equations (17) and (18) represent the workstation startup constraints. When a workstation starts up, tasks must be assigned to the workstation, and the workstations must start up in sequence.

[0062] In equations (1)-(18), W is the total number of disassembly tasks and the maximum number of workstations. Among them, the task set of product 1 is W 1 ={1, 2, 3…, w 1}, the task set of product 2 is W 2 ={w 1 +1, w 1 +2, w 1 +3…, w 1 +w 2}, W ∈ W 1 ∪W 2; i and j are respectively the task numbers; n is the workstation number; m is the determination variable at the inlet and outlet ends of the disassembly line, m = 1 represents the inlet end, and m = 2 represents the outlet 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 type-A robot, k = 4 represents the operator is a type-B robot; M is a positive number greater than the total number of tasks; K 1 is the known set of junior complex tasks; K 2 is the known set of senior complex tasks; S 1 is the known set of type-A special tasks; S 2 is the known set of type-B special tasks; TP ij is the precedence relation matrix. If task i is the immediate predecessor task of task j, then TP ij = 1; CT is the cycle time; is the disassembly time of task i when the operator is k; is the cost of task i when the operator is k; r k is the standby cost when the operator at the workstation is l; P i is the disassembly revenue of task i; R k is the single acquisition cost of a robot of type k; 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 operating days of the disassembly line; G is the total number of parts disassembled by the disassembly line per day; is a binary variable, indicating that task i is assigned to the m port of the nth workstation with the operator k, otherwise x j mkn is a binary variable, x j mkn = 1 indicates that task j is assigned to the m port of the nth workstation with the operator k, otherwise x j mkn = 0; z kn is a binary variable, z kn = 1 indicates that the nth workstation is turned on and the operator is k, otherwise z kn = 0;

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

[0064] The algorithm adopted in the present invention is SA-KQL (Simulated annealing algorithm with K-means and Q-Learning, a simulated annealing algorithm based on K-means clustering and Q-Learning). Its basic idea is that under the framework of simulated annealing, the K-means clustering method is used to assist in dividing the individual states, and the action to be executed is selected according to the Q-table using the ε-greedy strategy.

[0065] Among them, in the encoding process, considering that the tasks in the U-shaped disassembly line need to be assigned to both sides of the inlet and outlet, a mixed way of positive and negative numbers is usually adopted to encode the task sequence. Positive numbers indicate that the task is assigned to the inlet side of the workstation, and negative numbers indicate that the task is assigned to the outlet side of the workstation. In terms of incomplete disassembly, in previous studies, a double-layer encoding method was usually adopted, and a second-layer task execution sequence was introduced to indicate whether the task is executed. For example Figure 2 As shown, in the present invention, a complex-domain task sequence encoding is introduced, which improves the original double-layer encoding to a single layer. Whether the task is executed is indicated by judging whether the imaginary part of the task sequence is 0. The white tasks in the figure, that is, the imaginary task numbers, indicate the tasks that are not executed, and the colored tasks, that is, the real task numbers, indicate the tasks that are executed.

[0066] Generally speaking, problems of the DLBP (Disassembly Line Balancing Problem) type need to satisfy the task precedence relationship constraints. For the convenience of computer calculation, a combined precedence relationship matrix as shown in Figure 3 (a) is used to represent the relationship between the combined tasks. On this basis, the encoding process of the present invention is encoded by a two-stage encoding method. The first stage generates an initial task sequence. The specific operation process is as follows: Randomly select tasks without a direct predecessor or a direct successor for assignment. If the task has no direct predecessor, the task is assigned to the inlet side of the U-shaped disassembly line, and the task number is a positive number; if the task has no direct successor, the task is assigned to the outlet side of the U-shaped disassembly line, and the task number is a negative number. Then update the TP matrix, and repeat the above operation until all tasks are assigned. The second stage selects the tasks that are not executed and generates the operator sequence. The specific operation process is as follows: First, determine the set E of tasks that can be not executed, randomly select any number of tasks not exceeding the number of E (including 0 tasks), and multiply the selected tasks and their direct successors in the task sequence by the imaginary unit i, indicating that they are tasks that are not executed. Finally, considering that the operators of each workstation are different, an operator sequence composed of random numbers from [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 shown in Figure 3 (b).

[0067] For the subsequent decoding process of the present invention, since the working task times of different operators are different, the number of tasks that can be assigned to the workstation cannot be determined when the operator is not determined for the workstation; since different tasks need to be assigned to different operators for operation, the operator in the workstation cannot be determined when the task attributes in the workstation are not determined. To solve this contradiction, the present invention proposes a decoding method with a variable operator sequence, and its main idea is: when there is a conflict between the tasks in the workstation and the operator k in the operator sequence, if there is an operator k_new who can execute all tasks within the beat time, then change the operator sequence and change the operator of the workstation to k_new; if all operators are unable to fully execute the tasks in the workstation, that is, the tasks in the workstation have multiple attributes, then still use the operator k as the operator of the current workstation. The specific operations are as follows: After the ws-th workstation is started, first determine the set I of executable tasks that satisfy the beat time constraint under the execution times of each operator k , check whether there are task conflicts in the set I of tasks executed by the operator ZL(ws) in the operator sequence ZL ZL(ws) . If not, update the task index i = i + size(I ZL(ws) ); if there are task conflicts, then check whether there are no task attribute conflicts in the remaining I k . If so, change the operator ZL(ws) = k_new and update the task index; if there is no k_new, then still use the original operator to execute the tasks until the previous task of the task attribute conflict task. In special cases, if the first task of the workstation is the task attribute conflict task and the original operator cannot be used, then change the operator to enable the workstation to start smoothly. The specific operations are shown in the pseudo-code in Table 1:

[0068] Table 1 Pseudo-code for decoding operations

[0069]

[0070]

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

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

[0073] Step 2) Select the parameters with p-value less than 0.05 in the multivariate analysis of variance for the orthogonal experiment with hypervolume as the index, and take the maximum value of each parameter in the orthogonal experiment as the value of the parameter in the algorithm;

[0074] Step 3) Take the median of the value range of the parameters with p-value not less than 0.05 in the multivariate analysis of variance as their values in the algorithm.

[0075] Step S3: Select any individual current in the current population. Based on the characteristics of the external archive individuals, compare this individual current with a randomly selected external individual external in the same cluster in the external archive to determine the state of this individual current. According to the state of this individual current, execute the Q-learning algorithm, and then update the Q-table and according to the improved Metropolis criterion until all individuals in the current population have executed the Q-learning algorithm;

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

[0077] Among them, in order to eliminate the influence caused by different orders of magnitude of each objective value when calculating the Euclidean distance to the cluster center, the Z-score, that is, the standardized score method, is used to standardize each objective data, and the Z-score calculation formula is shown in Equation (24).

[0078]

[0079] In the formula, f a represents the a-th objective function value of the individual, μ represents the mean of the objective function data set, and σ represents the standard deviation of the data set.

[0080] In terms of the solution sequence, the sequence of problems in the present invention mainly has three characteristics: the order of task execution, the operator sequence, and the tasks not to be executed. All three characteristics will affect the final objective function value. Therefore, all three aspects should be considered simultaneously in the solution sequence. On this basis, the similarity between any individual current in the current population and a randomly selected external individual external in the same cluster in the external archive is judged in terms of the three individual characteristics of the order of task execution, the operator sequence, and the tasks not to be executed. After combining the similarity situations in each sequence characteristic aspect, it is classified according to the fitness value of the K clusters in the K-mean clustering of the individual current, and the state of any individual current in the current population is obtained;

[0081] Among them, the similarity judgment in terms of the individual characteristic of the order of task execution is shown in Equation (19):

[0082]

[0083] In Equation (19), D hm _XL is the similarity degree in the order of task execution between any individual current in the current population and a randomly selected external individual external in the same cluster in the external archive; W is the total number of disassembly tasks and the maximum number of workstations. Among them, the task set of Product 1 is W 1 ={1, 2, 3…, w 1}, the task set of Product 2 is W 2 ={w 1 +1, w 1 +2, w 1 +3…, w 1 +w 2}, W ∈ W 1 ∪W 2 ; i is the task number. When the i-th task or operator is the same, i = 1, otherwise i = 0; XL_d extrnal [i] is the task execution sequence of a randomly selected external individual external in the same cluster in the external archive; 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 ), that is, when the similarity degree of the task execution sequences of the individual current and the 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 in terms of the individual characteristic of the operator sequence is shown in Equation (20):

[0086]

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

[0088] When D hm _ZL ≥ size(ZL_d extrnal ), that is, when the similarity of the operator sequences 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;

[0089] The similarity judgment of the characteristics of task individuals not to be executed is shown in formula (21):

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

[0091] △J is the similarity in the tasks not to be executed between any individual current in the current population and a random external individual external in the same cluster in the external archive; J extrnal is the number of tasks not to be executed of a random external individual external in the same cluster in the external archive; J current is the number of tasks not to be executed of any individual current in the current population;

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

[0093] A specific state space can be designed through fitness values and sequence similarities. The fitness values have a total of K states classified into K clusters, and each of the three sequence similarities has 2 states - similar and dissimilar. In this case, combining the three sequences can result in a total of 8 states. Then, multiplying and combining these 8 states obtained after combination with the K states obtained from the K-cluster classification gives 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 possess, and the sequence of any individual current will belong to one of these 8×K states. The schematic diagram of the specific state division process is as follows Figure 4 shown.

[0094] In terms of the action design of Q-learning, the actions (perturbation operations) in the traditional DLBP optimization algorithm usually only consider the perturbation of the task sequence. However, in the problem addressed in the present invention, the objective function value is not only related to the task sequence, and perturbations to non-executing tasks (the imaginary part in the task sequence) and the operator sequence also need to be designed.

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

[0096] A mutation operation that can change the execution on the inlet and outlet sides of a task. First, randomly select whether to change the inlet and outlet sides of the mutated task, then determine the mutable section of the task on this side, and finally randomly insert the task into the mutable section. The specific operation is as follows Figure 5 shown - first, randomly select a task with equal probability (such as task 3 in the figure), then generate a random number rand between 0 and 1. If rand is less than 0.5, then change its inlet and outlet ends, that is, the positive and negative signs. If rand is greater than 0.5, then do not change the inlet and outlet ends. The mutable section is determined by the preceding task and the succeeding task, and the mutated task must still satisfy the precedence relation constraint.

[0097] The crossover action is designed as:

[0098] The crossover operation represents the process of the current individual learning from the selected external individual, and the influence of the inlet and outlet sides also needs to be considered. Randomly select a part of the task sequence on the inlet side or outlet side of the current individual, rearrange it according to the actual execution order of the external individual, and then return it to the original sequence. The specific operation is as follows Figure 6 shown - first, randomly select the positive number sequence representing the inlet end or the negative number sequence representing the outlet end, and then randomly select a segment from the randomly selected sequence and rearrange it according to the order in the external sequence.

[0099] The action of changing whether a task is executed is designed as:

[0100] To consider the influence of the task sequence features that are not executed, i.e., the tasks with imaginary parts in the sequence, a variable action of whether a task is executed is designed. Randomly select a task from the set of non-attribute tasks and its immediate predecessor tasks (i.e., tasks that can be not executed). If the task is executed in the original sequence, then the task and its subsequent tasks are multiplied by the imaginary unit i to be converted into tasks that are not executed. If the task is not executed, then the imaginary unit of the task and its immediate predecessor tasks is removed to be converted into executable tasks. The specific operation is as shown in Figure 7 (a) below.

[0101] The design of the operator sequence action is changed to:

[0102] Considering the influence of the operator sequence on the individual fitness value, a mutation operation on the operator sequence is designed. Since the actual number of workstations in the decoding of this problem is less than the length of the operator sequence, to prevent invalid mutations, according to the workstation number objective function f 1 _external of the external archive individual external, randomly select one of the first f 1 _external operators in the current individual operator sequence for perturbation and change it to another operator. The specific operator is as shown in Figure 7 (b) below.

[0103] Furthermore, to increase the perturbation of the action to the original individual, the four actions of mutation, crossover, changing whether a task is executed, and changing the operator sequence are combined in pairs and executed simultaneously to obtain eight actions: mutation, crossover, changing whether a task is executed, changing the operator sequence, mutation and changing whether a task is executed simultaneously, mutation and changing the operator sequence simultaneously, crossover and changing whether a task is executed simultaneously, and crossover and changing the operator sequence simultaneously.

[0104] In Q-learning, the cumulative reward of the action-state pair is used to obtain the maximum benefit during the iterative process. The Q(state, action) value represents the Q value of the action corresponding to the state 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 iterative process, and its update formula is as shown in Equation (25):

[0105]

[0106] In the formula, 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, and its value is positively correlated with the historical learning effect. γ is the discount factor representing the influence of the reward on the new state.

[0107] Due to the multi-objective nature of this problem, the design of the reward and punishment value r refers to the Pareto method. When the objective function value after the action is executed dominates the original individual, r = 2; when neither of them dominates the other, r = 1; otherwise, r = -1. This algorithm adopts a dynamic ε-greedy strategy, and the value of ε increases continuously with the algorithm iteration, indicating that as the Q-table is continuously learned, the action selection of QL should trust the Q-table more. Among them, the actions selected by the dynamic ε-greedy strategy mainly include 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 in this state. Among them, the calculation formula of the ε value is shown in Equation (22):

[0109]

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

[0111] In addition, since the present invention relates to a multi-objective problem, an improved Metropolis criterion is adopted to judge whether to accept a new solution when updating the population pop. Among them, the probability of accepting a new solution is shown in Equation (23):

[0112]

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

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

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

[0116] The following is the effect display in the actual operation of this embodiment:

[0117] Program test environment:

[0118] The simulation calculation environment for the numerical experiments in this chapter is the Windows 11 system, AMD Ryzen 7 6800H CPU, 3.20 GHz, 16 GB of memory, and the Gurobi 10.0.2 exact solver is used for operation and solution.

[0119] Large-scale instance application:

[0120] The U-shaped disassembly line is suitable for the disassembly of medium and small-sized electrical and electronic products. Therefore, this problem model is applied to a disassembly line with a mixed flow of monitors and laptops. 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 Figure 8 shown, and the specific disassembly information is shown in Table 2. The disassembly time and cost sequence in the table are: junior workers, senior workers, type A robots, and type B robots. The cycle time is CT = 45s.

[0121] Table 2 Disassembly information of large-scale instances

[0122]

[0123]

[0124] As the problem scale increases, the solution time of the exact solver grows exponentially, and the exact solver can no longer solve large-scale instances within a limited time. Therefore, to obtain better results, it is necessary to determine the SA-KQL algorithm parameters. The hypervolume index (HV) is used as the evaluation index for the algorithm results. The higher the HV value, the better the convergence and solution quality of the solution. 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*, and the reference point is determined by the worst values of each sub-objective of different problems. The worst values of each objective in a randomly generated initial population are selected and rounded up as the reference point for this instance calculation. In this instance, the reference point [15, 6000, -100000] is selected for calculation.

[0127] Multivariate Analysis of Variance (MANOVA) can simultaneously consider the mutual relationships of multiple variables and the interactions of each independent variable, thereby providing multi-dimensional statistical data analysis, which is of great significance in statistics. The 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 a total of 3 6For each parameter combination, the HV value was calculated 10 times and the average value was taken, and a total of 10 * 3 6 = 7290 times. Multivariate analysis of variance was performed on the calculation results, 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 p-values of parameters T initial , q, maxL, and maxL * α are all less than 0.05, which proves that they have an impact on the algorithm performance. Since the solution space of large-scale problems is large and the set value of the population size is much smaller than the size of the solution space, it can no longer have an impact on the solution performance; the p-value of the single factor of the learning rate α is greater than 0.05, but the p-value of the interaction with maxL is less than 0.05. This is because the impact of α on the algorithm performance is non-linear and is not easily found in single-factor analysis, but becomes significant when considering the interaction. To select the best parameter combination, an orthogonal experiment was performed on the above 4 algorithm parameters, and the evaluation index was HV. The experimental results are as Figure 9 shown.

[0132] According to Figure 9 shown, the parameter values with better solution effects in the orthogonal experiment results were selected as the algorithm parameters. For the parameters with p-values greater than 0.05, the intermediate values within the parameter range were selected as the algorithm parameter values.

[0133] To better prove the superiority of the algorithm proposed in this study, it was compared with the more advanced algorithms in current human-robot collaboration research: multi-objective enhanced differential evolution algorithm (MEDE), multi-objective shuffled frog leaping algorithm (SFLA), and multi-objective modified teaching and learning optimization algorithm (MTLO) using the above large-scale examples. The results are shown in Table 4. The data in the table are the 10 final solutions sorted by crowding distance for each algorithm. The bold data are the optimal values obtained for the objectives, and the underlined data indicate that the solution is dominated by the solutions of other algorithms.

[0134] Comparison of Algorithm Solving Results in Table 4

[0135]

[0136]

[0137] It is found from the comparison results that from the perspective of single objective, the optimal value of the idle time balance index in large-scale instances is still 0, indicating that the combination of incomplete disassembly and multi-type operator U-shaped disassembly line can effectively improve the flexibility of the disassembly line. At the same time, for f 1 Only SA-KQL and MTLO obtained the optimal solution 6. For f 3 Only SA-KQL and SFLA obtained the optimal solution 286930. And for f 2 Only SA-KQL obtained the optimal solution 0, indicating that SA-KQL has strong depth search ability. From the multi-objective aspect, 8, 9, and 10 final solutions of the three comparison algorithms were dominated by other algorithms respectively, while no solution of the SA-KQL algorithm was dominated by other algorithms, indicating that the non-dominated solutions obtained by SA-KQL are closer to the true Pareto front. In summary, SA-KQL has better solving performance than classical algorithms.

[0138] As mentioned above, the above are only the preferred specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed in the embodiments of the present invention should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. An optimization method for the balance 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 of the U-shaped disassembly line balance problem considering the differences of mixed flow operators is established, and its objective function is: F=min[ f 1 ,f 2 ,-f 3] (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, j are the task numbers respectively; n is the workstation number; m is the judgment variable of 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 a known set of primary complex tasks; K2 is a known set of advanced complex tasks; S1 is a known set of class A special tasks; S2 is a 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 beat 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 purchase cost of a robot of type k; 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 every day; is a binary variable, indicates that task i is assigned to port m of the nth workstation whose operator is k, otherwise is a binary variable, indicates that task j is assigned to port m of the nth workstation whose operator is k, otherwise 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 algorithm parameters, encode and generate the initial population and its corresponding external archives, perform K-mean clustering on the external archives, and start the Markov chain; Step S3: Select any individual current in the current population, and compare the individual current with a random external individual external in the same cluster in the external archive based on the individual characteristics of the external archive, determine the state of the individual current, perform Q-learning operation according to the state of the individual current, and then update the Q table and perform Q-learning operation according to the improved Metropolis criterion until all individuals in the current population have undergone Q-learning operation; 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 through the Pareto method and crowding distance strategy. Step S5: Re-cluster the external archives, 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 greater, end the algorithm and output the final solution.

2. The method for optimizing the balance 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 balance 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 the parameters with p-value less than 0.05 in the multivariate variance analysis 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 analysis of variance is used as its value in the algorithm.

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

5. The method for optimizing the balance problem of a mixed flow U-shaped disassembly line considering operator differences according to claim 4, characterized in that: The method for determining the state 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 three individual characteristics: the order of executing tasks, the sequence of operators, and the tasks not executed. After combining the similarity of each sequence characteristic, combine the fitness value classification of the K clusters in the K-mean clustering of the individual current to obtain the state 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 archive; 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 a sequence of operators of a random external individual in the same cluster in the external archive; 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 the individual characteristics of the 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 balance 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 balance 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 executed, changing an operator sequence, mutation and changing whether a task is executed simultaneously, mutation and changing an operator sequence simultaneously, crossover and changing whether a task is executed simultaneously, and crossover and changing an operator sequence simultaneously.

8. The method for optimizing the balance 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 adopts 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, where the calculation formula of the ε value 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 balance 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 balance 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, and 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 beat time, the operator sequence is changed and the operator of the workstation is changed to k_new; if all operators cannot completely execute the tasks in the workstation, that is, the workstation has tasks with multiple attributes, then the operator k is still used as the operator of the current workstation.

Images