A distributed blocking flow shop scheduling optimization system based on a co-evolution algorithm

Abstract

The application provides a distributed blocking flow shop scheduling optimization system based on a co-evolution algorithm, and designs a knowledge-driven cross-region interactive learning mechanism assisted co-evolution algorithm to optimize the distributed blocking flow shop scheduling under the energy consumption constraint condition. A heuristic method based on the problem characteristics is designed to reduce the energy consumption and total delay in the production process. In order to improve the diversity and convergence rate of the population, the application designs a cross-region interactive learning mechanism guided by the optimal state value function of reinforcement learning, realizes comprehensive cooperation from the algorithm level, the parameter level and the individual level, and improves the accuracy and efficiency of the solution. In order to improve the solution quality, a poor solution repair strategy and an individual rebirth mechanism are provided. In addition, the energy consumption is reduced through the energy-saving operation on the critical path. Through comparison and verification on the continuous optimization problem test set and the examples composed of different numbers of workshops, workpieces and machines, the optimization effect of the optimization system designed by the application is better than that of other optimization systems.

Description

Technical Field

[0001] This invention belongs to the field of manufacturing production scheduling and route planning technology, and relates to a distributed congested flow shop scheduling optimization system based on a co-evolutionary algorithm. Background Technology

[0002] The rapid development of society has led to optimization problems in intelligent computing being applied to more fields such as continuous optimization, engineering applications, and intelligent transportation. Furthermore, as the complexity of optimization problems increases, different improved optimization systems have been proposed based on the original system framework. The main objective of optimization is to find the optimal solution from a set of candidate solutions, typically a minimization or maximization problem.

[0003] Many typical problems in complex process industrial systems can be abstracted into distributed congested flow shop scheduling problems. For example, in the steel industry, industrial waste treatment, and non-ferrous metal smelting processes, the selection of the factory, the allocation of machines, the machine operating speed, and the sequence of workpieces directly affect the enterprise's production efficiency, revenue, and energy consumption. Therefore, designing a distributed congested flow shop scheduling optimization system with high production efficiency and low energy loss has significant research significance and value.

[0004] Continuous optimization and distributed shop floor scheduling optimization problems in complex process industrial systems are often complex and multimodal. Intelligent optimization algorithms, due to their wide applicability, robustness, adaptability, self-learning capabilities, and independence from the nature of the problem, are widely used in optimization systems. Typical intelligent optimization algorithms include Differential Evolution (DE), Evolutionary Algorithm (EDA), and Genetic Algorithm (GA). DE guides evolutionary behavior through mutation operators and difference terms, including mutation, crossover, and selection. It evaluates each solution based on the objective function, uses a greedy strategy to select a better solution, and iteratively approaches the optimal solution. DE is widely used in solving complex optimization problems due to its good performance, few parameters, simple structure, and strong global search capability. EDA is a new evolutionary algorithm developed from the Genetic Algorithm, based on probability theory and statistics, involving mathematical statistics and intelligent computing. EDA constructs a probabilistic model based on selected dominant individuals, thereby generating high-quality solutions that replace inferior solutions in the original population, and can find the global optimum with a higher probability.

[0005] Co-evolutionary algorithms leverage the strengths of each algorithm through appropriate collaborative strategies, eliminating their individual limitations and even achieving complementary advantages. This algorithmic framework can formulate optimization schemes for potentially adaptive search regions, improving convergence speed or search accuracy to varying degrees. Reinforcement learning is used to learn sequential decision-making tasks. Agents optimize their behavior through continuous experimentation. Agents learn by interacting with dynamic environments and evaluate the current environment through reward feedback. In recent years, the combination of reinforcement learning and intelligent optimization algorithms has been successfully used to solve various optimization problems. Summary of the Invention

[0006] The purpose of this invention is to address the problems existing in the prior art. With the support of the "Innovation Star" project for postgraduate students funded by the Gansu Provincial Department of Education (Project No. 2023CXZX-476), this invention conducts relevant research and proposes technical solutions, specifically a distributed congested flow shop scheduling optimization system based on a co-evolutionary algorithm. This invention designs a cross-regional interactive learning mechanism driven by problem characteristics and guided by reinforcement learning, achieving a comprehensive collaborative optimization method at the algorithm, parameter, and individual levels. Figure 1 As shown, detailed explanations can be found in steps 1-3. Performance verification was conducted in a continuous optimization problem, and the system was applied to a distributed congested flow shop scheduling problem with energy constraints, demonstrating that the optimization system designed in this invention can improve production efficiency while reducing energy consumption.

[0007] A distributed congested flow shop scheduling optimization system based on a co-evolutionary algorithm is characterized by the following: a heuristic method based on problem characteristics that considers total energy consumption and total delay; a cross-regional interactive learning mechanism guided by the optimal state value function based on reinforcement learning; a suboptimal solution repair strategy and individual regeneration mechanism for specific phenomena of the optimization problem; and energy-saving operations based on the critical path. The specific interactive methods of the key ideas of this invention are as follows: Figure 2 As shown, the specific steps include 1-3:

[0008] Step 1: In the initialization method considering total delay and total energy consumption, the position and operating speed of each workpiece are randomly initialized. This method considers the effects of front-end delay, blocking time, and idle time, especially for the first position, using the workpiece with the minimum delay time. Based on the workpiece sequence, the factory with the maximum total delay is selected as the critical factory. By optimizing the total delay of the critical factory, the total delay of the distributed blocking flow shop with energy consumption constraints is optimized.

[0009] Step 2: The cross-regional interactive learning mechanism guided by the optimal state value function based on reinforcement learning, and the inferior solution repair strategy and individual regeneration mechanism for specific phenomena of the optimization problem, include steps A1-A7:

[0010] Step A1: Divide the area;

[0011] The population is divided into regions based on the distance between individuals and elite solutions. These regions are further divided into three areas from farthest to closest: the region containing the closest 1 / 3 of individuals is defined as the strengthening region; the region containing the farthest 1 / 3 of individuals is defined as the decaying region; and the remaining regions at intermediate distances are defined as the stable region. The distance is defined as shown in Formula 1-2.

[0012] d(X_i,X_best)=‖X_i-X_best‖,i∈{1,2,…,NP} Formula 1

[0013] NP = NP_total / 3 (Formula 2)

[0014] Where X_best is the elite solution, i.e. the individual with the smallest fitness value, X_i is the vector of the i-th individual, NP is the number of individuals in each region, and NP_total is the total population size.

[0015] Step A2: Establish a reinforcement learning mechanism based on the ε-greedy policy and simulated annealing, specifically including sub-steps A2.1-A2.5:

[0016] Step A2.1: Considering that the evolutionary process is a continuous learning process, the discount factor γ emphasizes the contribution to the total reward within different iterations. Therefore, γ does not take a fixed value, but is updated adaptively;

[0017] γ = 0.5 × (0.9 - 0.6 × g / G_max) (Formula 3)

[0018] Where g is the current algebra and G_max is the maximum number of iterations;

[0019] Step A2.2: Approximately use the optimal state value function as the learning objective, and iteratively update it according to Formula 4;

[0020] Q_g←Q_g+α[R_(g+1)+γ*ΔQ_(g+1)] Formula 4

[0021] Where α is the learning rate, R_g is the reward for the g-th generation, and ΔQ_(g+1) is the increment of Q.

[0022] Step A2.3: R_g is set according to the distance, based on the distance between the individual and the optimal individual in each region as shown in Formula 5;

[0023] L(X_(i,k),X_(best,k))=‖X_(i,k)-X_(best,k)‖,k={1,2,3} Formula 5

[0024] Where X_(i,k) is the i-th individual in the k-th region, and X_(best,k) is the best individual in the k-th region;

[0025] Step A2.4: Reward or punish individuals based on their distance from the optimal solution to obtain R_g. A larger R_g indicates closer proximity to the optimal solution; a smaller R_g indicates further away from the optimal solution.

[0026] Step A2.5: Introduce the ε-greedy strategy and dynamically adjust the ε value using simulated annealing. The temperature is updated according to formula 6-7;

[0027] T_0 = T_max (Formula 6)

[0028] T_(g+1)=T_min+θ×(T_g-T_min) Formula 7

[0029] Where T_0 is the initial temperature, which reaches its maximum value, T_min is the minimum temperature, and θ∈[0,1] is the annealing factor. During the iteration process, the changes in the Q-value table are considered, and the ε-value is defined using the difference in Q-values ​​and the dynamic temperature value.

[0030] ε=exp((-(ΔQ_g)) / T_g) Formula 8

[0031] Where ΔQ_g is the increment of Q in the g-th generation, and T_g is the temperature value in the g-th generation;

[0032] Step A3: Implement a cross-regional dynamic resource allocation scheme;

[0033] The parameters for the three regions are updated according to different formulas, specifically including sub-steps A3.1-A3.4:

[0034] Step A3.1: The parameters in the mutation strategy for the enhanced region are updated according to formula 9-14;

[0035] F_(i,1)∈randn_i(μ_(F_1),0.1)

[0036] Cr_(i,1)∈randn_i(μ_(Cr_1),0.1) Formula 9

[0037] F_(i+1,1)'=c*F_(i+1,1)+(1-c)*F_(i,1)

[0038] Cr_(i+1,1)'=c*Cr_(i+1,1)+(1-c)*Cr_(i,1) Formula 10

[0039] c = 0.9 - 0.9 * 10^(-5) * g / G_max (Formula 11)

[0040] μ_(F_1)=(1-c)*μ_(F_1)+c*mean_A(ω_i*F_(i,1)) Formula 12

[0041] μ_(Cr_1)=(1-c)*μ_(Cr_1)+c*mean_A(ω_i*Cr_(i,1)) Formula 13

[0042] ω_i=(Δfit_i) / (∑(Δfit_k) Formula 14

[0043] Where F_i and Cr_i are the mutation factor and crossover factor, μ_F is the mean of F_i, μ_Cr is the mean of Cr_i, fit_i is the fitness value of the i-th individual, Δfit_i is the change in fitness value, and mean_A represents the arithmetic mean. F_i and Cr_i both follow a Gaussian distribution in this region; the change in fitness value is used as a weighting coefficient, and the arithmetic mean is used to calculate the mean.

[0044] Step A3.2: The parameters in the mutation strategy of the stable region are updated according to formula 15-19;

[0045] F_(i,2)∈randc_i(μ_(F_2),0.1)

[0046] Cr_(i,2)∈randn_i(μ_(Cr_2),0.1)) Formula 15

[0047] F_(i+1,2)'=c*F_(i+1,2)+(1-c)*F_(i,2)

[0048] Cr_(i+1,2)'=c*Cr_(i+1,2)+(1-c)*Cr_(i,2) Formula 16

[0049] μ_(F_2)=(1-c)*μ_(F_2)+c*mean_L(ω_i*F_(i,2)) Formula 17

[0050] μ_(Cr_2)=(1-c)*μ_(Cr_2)+c*mean_A(ω_i*Cr_(i,2)) Formula 18

[0051] mean_L(ω_i*F_(i,2))=(∑ω_i*F_(i,2)^2) / (∑ω_i*F_(i,2)) Formula 19

[0052] Where F_i follows a Cauchy distribution and the mean is the Lymer mean, mean_L(ω_i*F_(i,2)), and Cr_i follows a Gaussian distribution and is calculated using the arithmetic mean.

[0053] Step A3.3: The parameters in the mutation strategy of the decay region are updated according to formulas 20-25;

[0054] F_(i,3)∈randc_i(μ_(F_3),0.1)

[0055] Cr_(i,3)∈randc_i(μ_(Cr_3),0.1) Formula 20

[0056] F_(i+1,3)'=c*F_(i+1,3)+(1-c)*F_(i,3)

[0057] Cr_(i+1,3)'=c*Cr_(i+1,3)+(1-c)*Cr_(i,3) Formula 21

[0058] μ_(F_3)=(1-c)*μ_(F_3)+c*mean_L(ω_i*F_(i,3)) Formula 22

[0059] μ_(Cr_3)=(1-c)*μ_(Cr_3)+c*mean_l(ω_i*Cr_(i,3)) Formula 23

[0060] mean_L(ω_i*F_(i,3))=(∑ω_i*F_(i,3)^2) / (∑ω_i*F_(i,3)) Formula 24

[0061] mean_L(ω_i*Cr_(i,3))=(∑ω_i*Cr_(i,3)^2) / (∑ω_i*Cr_(i,3) Formula 25

[0062] Both F_i and Cr_i follow a Cauchy distribution, and the mean is calculated using the Lemo mean, which yields a larger value.

[0063] Step A3.4: Average fitness of individuals per generation in each region. and the average fitness after mutation and crossover. Calculate according to formulas 26 and 27 respectively;

[0064]

[0065]

[0066] Step A4: When the conditions for successful mutation are met as shown in Formula 28, execute the self-feedback strategy;

[0067]

[0068] If formula 28 is satisfied, the generation has successfully mutated, and parameter self-feedback is executed. The self-feedback design for the three regions is shown in formulas 10, 16, and 21. Otherwise, it is executed according to formulas 9, 15, and 20.

[0069] Step A5: Improve the individuals with small Q-values ​​and those that are too densely packed in the region in the reinforcement learning mechanism by implementing a suboptimal solution improvement strategy, specifically including sub-steps A5.1-A5.6:

[0070] Step A5.1: Improve the individuals with small Q values ​​and those that are too densely packed in the region in the reinforcement learning mechanism, and use them in the population evolution process;

[0071] The density of individuals is measured by the average distance between an individual and the centroid of a region, as well as the average fitness value. The centroid of a region is obtained according to Formula 29.

[0072] m_c=(∑X_i) / NP Formula 29

[0073] Step A5.2: The average distance from each individual to its respective region's centroid is calculated using formula 30;

[0074] d-=(∑‖X_i-m_c‖) / NP Formula 30

[0075] Step A5.3: The average fitness increment is defined according to Formula 31;

[0076] Δf=(∑((f(X_i)-f(m_c))) / NP (Formula 31)

[0077] Step A5.4: The density DD is calculated according to formula 32;

[0078] DD=Δf*(1+λ / d-) Formula 32

[0079] Where Δf represents the change in fitness value;

[0080] Step A5.5: The value of λ is calculated according to formula 33;

[0081] λ=(∑((U_d-L_d) / NP)^2)^1 / 2 Formula 33

[0082] Where D represents dimension, L_d represents lower bound, and U_d represents lower bound.

[0083] Step A5.6: Individuals with small Q values ​​and overly dense populations are defined as inferior solutions, denoted by I_i, and the repair strategy is Equation 34:

[0084] I_i=I_i+ζ*(m_c-I_i),i∈{1,2,…,NP} Formula 34

[0085] Where I_i is the i-th repaired solution, and ζ∈[0,1] is used to control the degree to which the inferior solution is close to the centroid.

[0086] Step A6: When the fitness value has not been updated for several consecutive generations, execute the individual regeneration mechanism, specifically including sub-steps A6.1-A6.3:

[0087] Step A6.1: Use the Cholesky decomposition method to generate a low-order triangular matrix M from the covariance matrix C;

[0088] C = M * M^T (Formula 35)

[0089] Step A6.2: Sample a 1×D matrix R using the standard normal distribution N(0,1);

[0090] Step A6.3: Generate a new individual according to formula 36;

[0091] X_new=μ+η(g)*M*R Formula 36

[0092] Where η(g) is the scaling factor, and η(g) > 1. Update η(g) according to formulas 37-40.

[0093] η(g)=η_min+((G_max-g) / (G_max))^β*(η_max-η_min) Formula 37

[0094] η_max=(U_d-L_d)*0.01 Formula 38

[0095] η_min=η_max / 10 (Formula 39)

[0096] β=1-e^(g / (0.15*G_max)) Formula 40

[0097] Where η_max and η_min represent the maximum and minimum allowed values, determined according to the size of the search space, with L_d representing the lower bound and U_d representing the upper bound, G_max being the maximum number of iterations, and β being the correction parameter for nonlinear adjustment of η.

[0098] Step A7: Once the development phase begins, perform a distance-based EDA local search, following formulas 41-42.

[0099] X_i(g+1)=N(μ_i,δ_i) Formula 41

[0100] Where N(μ_i,δ_i) is a Gaussian distribution with mean μ_i and variance δ_i;

[0101] δ_i=τ*‖X_best-X_i(g)‖ Formula 42

[0102] Where τ is the scaling factor.

[0103] Step 3: Energy-saving operations based on the critical path, proceeding according to steps B1-B3;

[0104] Step B1: Reducing the processing speed of some workpieces can reduce energy consumption;

[0105] Step B2: Reduce energy consumption and achieve energy saving by decreasing the speed of non-critical components on the critical path;

[0106] Step B3: If both total energy consumption and total delay decrease, or if either indicator decreases, then accept the solution; otherwise, discard the new solution and retain the original solution.

[0107] Compared with the prior art, the beneficial improvements of the present invention are:

[0108] (1) With the goal of reducing energy consumption and total delay in the production process, a heuristic method based on the characteristics of the problem was designed to achieve a comprehensive balance and bidirectional optimization of production efficiency and energy consumption indicators.

[0109] (2) In order to improve the diversity and convergence rate of the population, this invention designs a cross-regional interactive learning mechanism guided by the optimal state value function of reinforcement learning, which achieves comprehensive collaboration at the algorithm level, parameter level and individual level, thereby improving the accuracy and efficiency of the solution.

[0110] (3) To improve the quality of the solution, a suboptimal solution repair strategy and an individual regeneration mechanism were proposed. In addition, energy consumption was reduced by energy-saving operations on the critical path, which improved production efficiency and achieved the goal of energy saving. Attached Figure Description

[0111] Figure 1 This is a flowchart of the improved method proposed in this invention;

[0112] Figure 2 This is a design concept diagram of the present invention;

[0113] Figure 3 This is a box diagram of our optimization system and the comparative optimization system used in the continuous optimization problem;

[0114] Figure 4 It is an energy-saving operation diagram of the factory's critical path. Detailed Implementation

[0115] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0116] A distributed congested flow shop scheduling optimization system based on a co-evolutionary algorithm is characterized by comprising: heuristic steps based on problem characteristics considering total energy consumption and total delay; cross-regional interactive learning steps guided by the optimal state value function based on reinforcement learning; and inferior solution repair steps and individual regeneration steps for specific phenomena of the optimization problem; and energy-saving operations based on the critical path; specifically including the following steps:

[0117] Step 1: In the initialization step considering total delay and total energy consumption, the position and operation speed of each workpiece are randomly initialized. In this step, the effects of front-end delay, blocking time, and idle time are considered, especially for the first position, which uses the workpiece with the minimum delay time. According to the workpiece sequence, the factory with the maximum total delay is selected as the critical factory. The total delay of the distributed blocking flow shop with energy consumption constraints is optimized by optimizing the total delay of the critical factory.

[0118] Step 2: The cross-regional interactive learning step guided by the optimal state value function based on reinforcement learning, and the inferior solution repair step and individual regeneration step for specific phenomena of the optimization problem, including steps A1-A7:

[0119] Step A1: Divide the population into regions; divide the population into three regions based on the distance between individuals and elite solutions, from farthest to closest. The region containing the closest 1 / 3 of individuals is defined as the enhancement region, the region containing the farthest 1 / 3 of individuals is defined as the decay region, and the remaining regions at intermediate distances are defined as the stable region. The distance is defined as shown in Formula 1-2:

[0120] d(X_i,X_best)=‖X_i-X_best‖,i∈{1,2,…,NP} Formula 1

[0121] NP = NP_total / 3 (Formula 2)

[0122] Where X_best is the elite solution, i.e. the individual with the smallest fitness value, X_i is the vector of the i-th individual, NP is the number of individuals in each region, and NP_total is the total population size.

[0123] Step A2: Establish a reinforcement learning mechanism based on the ε-greedy policy and simulated annealing, specifically including sub-steps A2.1-A2.5:

[0124] Step A2.1: Considering that the evolutionary process is a continuous learning process, the discount factor γ emphasizes the contribution to the total reward within different iterations; therefore, γ does not take a fixed value, but is updated adaptively.

[0125] γ = 0.5 × (0.9 - 0.6 × g / G_max) (Formula 3)

[0126] Where g is the current algebra and G_max is the maximum number of iterations;

[0127] Step A2.2: Approximately use the optimal state value function as the learning objective, and iteratively update it according to Formula 4;

[0128] Q_g←Q_g+α[R_(g+1)+γ*ΔQ_(g+1)] Formula 4

[0129] Where α is the learning rate, R_g is the reward in the g-th generation, and ΔQ_(g+1) is the increment of Q;

[0130] Step A2.3: R_g is set according to the distance, based on the distance between the individual and the optimal individual in each region as shown in Formula 5;

[0131] L(X_(i,k),X_(best,k))=‖X_(i,k)-X_(best,k)‖,k={1,2,3} Formula 5

[0132] Where X_(i,k) is the i-th individual in the k-th region, and X_(best,k) is the best individual in the k-th region;

[0133] Step A2.4: Reward or punish individuals by comparing their distances to obtain R_g; the larger R_g is, the closer the individual is to the optimal solution; the smaller R_g is, the further the individual is from the optimal solution.

[0134] Step A2.5: Introduce the ε-greedy strategy and simulated annealing to dynamically adjust the ε value, and update the temperature according to formula 6-7;

[0135] T_0 = T_max (Formula 6)

[0136] T_(g+1)=T_min+θ×(T_g-T_min) Formula 7

[0137] Where T_0 is the initial temperature and reaches its maximum value, T_min is the minimum temperature, and θ∈[0,1] is the annealing factor; during the iteration process, the changes in the Q value table are considered, and the ε value is defined using the difference in Q values ​​and the dynamic temperature value;

[0138] ε=exp((-(ΔQ_g)) / T_g) Formula 8

[0139] Where ΔQ_g is the increment of Q in the g-th generation, and T_g is the temperature value in the g-th generation;

[0140] Step A3: Implement a cross-regional dynamic resource allocation scheme; the parameters for the three regions are updated according to different formulas, specifically including sub-steps A3.1-A3.4:

[0141] Step A3.1: The parameters in the mutation strategy for the enhanced region are updated according to formula 9-14;

[0142] F_(i,1)∈randn_i(μ_(F_1),0.1)

[0143] Cr_(i,1)∈randn_i(μ_(Cr_1),0.1) Formula 9

[0144] F_(i+1,1)'=c*F_(i+1,1)+(1-c)*F_(i,1)

[0145] Cr_(i+1,1)'=c*Cr_(i+1,1)+(1-c)*Cr_(i,1) Formula 10

[0146] c = 0.9 - 0.9 * 10^(-5) * g / G_max (Formula 11)

[0147] μ_(F_1)=(1-c)*μ_(F_1)+c*mean_A(ω_i*F_(i,1)) Formula 12

[0148] μ_(Cr_1)=(1-c)*μ_(Cr_1)+c*mean_A(ω_i*Cr_(i,1)) Formula 13

[0149] ω_i=(Δfit_i) / (∑(Δfit_k) Formula 14

[0150] Where F_i and Cr_i are the mutation factor and crossover factor, μ_F is the mean of F_i, μ_Cr is the mean of Cr_i, fit_i is the fitness value of the i-th individual, Δfit_i is the change in fitness value, and mean_A represents the arithmetic mean. F_i and Cr_i both follow a Gaussian distribution in this region; the change in fitness value is used as a weighting coefficient, and the arithmetic mean is used to calculate the mean.

[0151] Step A3.2: The parameters in the mutation strategy of the stable region are updated according to formula 15-19;

[0152] F_(i,2)∈randc_i(μ_(F_2),0.1)

[0153] Cr_(i,2)∈randn_i(μ_(Cr_2),0.1)) Formula 15

[0154] F_(i+1,2)'=c*F_(i+1,2)+(1-c)*F_(i,2)

[0155] Cr_(i+1,2)'=c*Cr_(i+1,2)+(1-c)*Cr_(i,2) Formula 16

[0156] μ_(F_2)=(1-c)*μ_(F_2)+c*mean_L(ω_i*F_(i,2)) Formula 17

[0157] μ_(Cr_2)=(1-c)*μ_(Cr_2)+c*mean_A(ω_i*Cr_(i,2)) Formula 18

[0158] mean_L(ω_i*F_(i,2))=(∑ω_i*F_(i,2)^2) / (∑ω_i*F_(i,2)) Formula 19

[0159] Where F_i follows a Cauchy distribution and its mean is the Lymer mean mean_L(ω_i*F_(i,2)), and Cr_i follows a Gaussian distribution and its arithmetic mean is used;

[0160] Step A3.3: The parameters in the mutation strategy of the decay region are updated according to formulas 20-25;

[0161] F_(i,3)∈randc_i(μ_(F_3),0.1)

[0162] Cr_(i,3)∈randc_i(μ_(Cr_3),0.1) Formula 20

[0163] F_(i+1,3)'=c*F_(i+1,3)+(1-c)*F_(i,3)

[0164] Cr_(i+1,3)'=c*Cr_(i+1,3)+(1-c)*Cr_(i,3) Formula 21

[0165] μ_(F_3)=(1-c)*μ_(F_3)+c*mean_L(ω_i*F_(i,3)) Formula 22

[0166] μ_(Cr_3)=(1-c)*μ_(Cr_3)+c*mean_l(ω_i*Cr_(i,3)) Formula 23

[0167] mean_L(ω_i*F_(i,3))=(∑ω_i*F_(i,3)^2) / (∑ω_i*F_(i,3)) Formula 24

[0168] mean_L(ω_i*Cr_(i,3))=(∑ω_i*Cr_(i,3)^2) / (∑ω_i*Cr_(i,3) Formula 25

[0169] Both F_i and Cr_i follow a Cauchy distribution, and the mean is calculated using the Lemo mean, which yields a larger value.

[0170] Step A3.4: Average fitness of individuals per generation in each region. and the average fitness after mutation and crossover. Calculate according to formulas 26 and 27 respectively;

[0171]

[0172]

[0173]

[0174] If formula 28 is satisfied, the generation has successfully mutated and parameter self-feedback is executed; otherwise, formulas 9, 15, and 20 are executed. The self-feedback design for the three regions is shown in formulas 10, 16, and 21.

[0175] Step A4: When the conditions for successful mutation are met, execute the self-feedback strategy;

[0176] Step A5: Improve the reinforcement learning mechanism by addressing individuals with small Q-values ​​and those that are too densely packed within a region, implementing a suboptimal solution improvement process, specifically including sub-steps A5.1-A5.6:

[0177] Step A5.1: Improve the individuals with small Q values ​​and those that are too densely packed in the region in the reinforcement learning mechanism, and use them in the population evolution process;

[0178] The density of individuals is measured by the average distance between an individual and the centroid of a region, as well as the average fitness value. The centroid of a region is obtained according to Formula 29.

[0179] m_c=(∑X_i) / NP Formula 29

[0180] Step A5.2: The average distance from each individual to its respective region's centroid is calculated using formula 30;

[0181] d-=(∑‖X_i-m_c‖) / NP Formula 30

[0182] Step A5.3: The average fitness increment is defined according to Formula 31;

[0183] Δf=(∑((f(X_i)-f(m_c))) / NP (Formula 31)

[0184] Step A5.4: The density DD is calculated according to formula 32;

[0185] DD=Δf*(1+λ / d-) Formula 32

[0186] Where Δf represents the change in fitness value;

[0187] Step A5.5: The value of λ is calculated according to formula 33;

[0188] λ=(∑((U_d-L_d) / NP)^2)^1 / 2 Formula 33

[0189] Where D represents the dimension, L_d represents the lower bound, and U_d represents the lower bound. Obviously, λ is adaptively adjusted with the search space.

[0190] Step A5.6: Individuals with small Q values ​​and overly dense populations are defined as inferior solutions, denoted by I_i, and the repair strategy is Equation 34:

[0191] I_i=I_i+ζ*(m_c-I_i),i∈{1,2,…,NP} Formula 34

[0192] Where I_i is the i-th repaired solution, and ζ∈[0,1] is used to control the degree to which the inferior solution is close to the centroid;

[0193] Step A6: When the fitness value has not been updated for several consecutive generations, execute the individual regeneration mechanism, specifically including sub-steps A6.1-A6.3:

[0194] Step A6.1: Use the Cholesky decomposition method to generate a low-order triangular matrix M from the covariance matrix C;

[0195] C = M * M^T (Formula 35)

[0196] Step A6.2: Sample a 1×D matrix R using the standard normal distribution N(0,1);

[0197] Step A6.3: Generate a new individual according to formula 36;

[0198] X_new=μ+η(g)*M*R Formula 36

[0199] Where η(g) is the scaling factor, η(g)>1, and η(g) is updated according to formula 37-40;

[0200] η(g)=η_min+((G_max-g) / (G_max))^β*(η_max-η_min) Formula 37

[0201] η_max=(U_d-L_d)*0.01 Formula 38

[0202] η_min=η_max / 10 (Formula 39)

[0203] β=1-e^(g / (0.15*G_max)) Formula 40

[0204] Where η_max and η_min represent the maximum and minimum allowed values, which are determined according to the size of the search space. L_d represents the lower bound, U_d represents the upper bound, G_max is the maximum number of iterations, and β is the correction parameter for nonlinear adjustment of η.

[0205] Step A7: Once the development phase begins, perform a distance-based EDA local search, following formulas 41-42.

[0206] X_i(g+1)=N(μ_i,δ_i) Formula 41

[0207] Where N(μ_i,δ_i) is a Gaussian distribution with mean μ_i and variance δ_i;

[0208] δ_i=τ*‖X_best-X_i(g)‖ Formula 42

[0209] Where τ is the scaling factor.

[0210] Step 3: Energy-saving operations based on the critical path, such as... Figure 4 As shown, different numbers represent different workpieces, and the specific steps are B1-B3.

[0211] Step B1: Reducing the processing speed of some workpieces can reduce energy consumption;

[0212] Step B2: Reduce the speed of non-critical components on the critical path to reduce energy consumption and achieve energy saving.

[0213] Step B3: If both total energy consumption and total delay decrease, or if either indicator decreases, then accept the solution; otherwise, discard the new solution and retain the original solution.

[0214] The optimized system proposed in this invention was verified in Example 1.

[0215] Step 1: Verify the performance of our optimized system compared to the comparative optimized systems in the table by comparing the mean and variance. The mean reflects the overall optimization performance, while the variance reflects the stability of the system. Evaluations were performed on 10-dimensional (10D), 30-dimensional (30D), 50-dimensional (50D), and 100-dimensional (100D) systems, as shown in Tables 1-4. The optimal results are marked in bold.

[0216] Table 1. Comparison of mean and variance of several typical optimization systems (10D)

[0217]

[0218] Table 2 Comparison of mean and variance of several typical optimization systems (30D)

[0219]

[0220]

[0221] Table 3. Comparison of mean and variance of several typical optimization systems (50D)

[0222]

[0223] Table 4. Comparison of mean and variance of several typical optimization systems (100D)

[0224]

[0225]

[0226] Step 2: The results above show that the optimization method proposed in this invention has better overall performance compared to the other four optimization methods.

[0227] Step 3: Box diagrams for optimization problems f3, f6, f14, and f26, such as... Figure 3 As shown, the horizontal axis represents each optimized system, and the vertical axis represents the error value. Figure 3 This reflects the stability effect of the optimized system. The stability of the optimization method proposed in this invention is optimal compared with other comparative optimization methods.

[0228] Step 4: Table 5 shows the significant differences between this optimization system and other optimization systems.

[0229] Table 5 Statistical test values ​​of eight optimization systems

[0230]

[0231] The optimized system proposed in this invention was verified in Example 2.

[0232] The proposed optimization system was validated against current classic optimizers on 720 instances with varying numbers of jobs, machines, and factories. The number of jobs and machines were set to {(20, 50, 100)*(5, 10, 20); 200*(10, 20); 500*20}, and the number of factories was set to 2, 3, 4, 5, 6, and 7. Non-dominated vector generation (ONVG) and the C metric were used to evaluate performance.

[0233] Table 6 Comparison of ONVG Indices for Four Optimization Systems

[0234]

[0235]

[0236] Table 6 shows that the result obtained by this optimization system is optimal.

[0237] Table 7 Nonparametric tests for the optimized system

[0238]

[0239] The results in Table 7 show that this optimized system is significantly different from other optimized systems.

[0240] In summary, the optimization system designed in this invention outperforms other optimization systems in this factory production instance.

[0241] The basic principles, main features, and advantages of the present invention have been described above in conjunction with the accompanying drawings. For those skilled in the art, modifications or variations can be made without departing from the principles of the present invention, and these improvements are also considered to be within the scope of protection of the present invention.

Claims

1. A distributed congested flow shop scheduling optimization system based on a co-evolutionary algorithm, characterized in that, include: Heuristic steps based on problem characteristics that consider total energy consumption and total delay, cross-regional interactive learning steps guided by the optimal state value function based on reinforcement learning, and inferior solution repair steps and individual regeneration steps for specific phenomena of optimization problems, and energy-saving operations based on the critical path. Specifically, the following steps are included: Step 1: In the initialization step that considers total delay and total energy consumption, the position and operating speed of each workpiece are randomly initialized. In this step, the effects of front-end delay, blocking time and idle time are considered. For the first position, the workpiece with the minimum delay time is selected. According to the workpiece sequence, the factory with the maximum total delay is selected as the critical factory. Step 2: The cross-regional interactive learning step guided by the optimal state value function based on reinforcement learning, and the inferior solution repair step and individual regeneration step for specific phenomena of the optimization problem, including steps A1-A7: Step A1: Divide the area; Step A2: Establish a reinforcement learning mechanism based on the ε-greedy policy and simulated annealing; Step A3: Implement a cross-regional dynamic resource allocation scheme; Step A4: When the conditions for successful mutation are met, execute the self-feedback strategy; Step A5: Improve the reinforcement learning mechanism by refining individuals with small Q-values ​​and those that are too densely packed within a region, and perform a suboptimal solution improvement process; Step A6: When the fitness value has not been updated for several consecutive generations, execute the individual regeneration mechanism; Step A7: Once the development phase begins, perform a distance-based local EDA search; Step 3: Energy-saving operations based on the critical path, proceeding according to steps B1-B3; Step B1: Reducing the processing speed of some workpieces can reduce energy consumption; Step B2: Reduce the speed of non-critical components on the critical path to reduce energy consumption and achieve energy saving. Step B3: If both total energy consumption and total delay decrease, or if either indicator decreases, then accept the solution; otherwise, discard the new solution and retain the original solution.

2. The distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that, In step 1: The factory with the largest total delay is designated as the critical factory, and the total delay of the distributed congested flow workshop with energy consumption constraints is optimized by optimizing the total delay of the critical factory.

3. The distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that, In sub-step A1 of step 2, regions are divided based on the distance between individuals in the population and elite solutions. These regions are divided into three areas from farthest to closest: the region containing the closest 1 / 3 of individuals is defined as the enhancement region, the region containing the farthest 1 / 3 of individuals is defined as the decay region, and the remaining regions at intermediate distances are defined as the stable region. The distance is defined as shown in Formula 1-2. d(X_i,X_best)=║X_i-X_bes║,i ∈ {1,2,…,NP}Formula 1 NP = NP_total / 3 (Formula 2) Where X_best is the elite solution, i.e. the individual with the smallest fitness value, X_i is the vector of the i-th individual, NP is the number of individuals in each region, and NP_total is the total population size.

4. A distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that, Sub-step A2 in step 2 specifically includes sub-steps C1-C5: Step C1: Considering that the evolutionary process is a continuous learning process, the discount factor γ emphasizes the contribution to the total reward within different iterations; therefore, γ does not take a fixed value, but is updated adaptively. γ = 0.5 × (0.9 - 0.6 × g / G_max) Formula 3; Where g is the current algebra and G_max is the maximum number of iterations; Step C2: Approximate the optimal state value function as the learning objective, and iteratively update it according to Formula 4; Q_g←Q_g+α[R_ (g+1)+γ* Δ Q_ (g+1)] Formula 4; Where α is the learning rate, R_g is the reward in the g-th generation, and ΔQ_(g+1) is the increment of Q; Step C3: R_g is set according to the distance, based on the distance between the individual and the optimal individual in each region as shown in Formula 5; L(X_ (i,k) ,X_ (best,k))=║X_ (i,k)-X_ (best,k) ║,k={1,2,3} Formula 5; Where X_(i,k) is the i-th individual in the k-th region, and X_(best,k) is the best individual in the k-th region; Step C4: Reward or punish individuals based on their distance to obtain R_g; the larger R_g is, the closer the individual is to the optimal solution; the smaller R_g is, the further the individual is from the optimal solution. Step C5: Introduce the ε_greedy strategy and simulated annealing to dynamically adjust the ε value, and update the temperature according to formula 6-7; Formula 6: T_0 = T_max; T_ (g+1)=T_min+θ×(T_g-T_min) Formula 7; Where T_0 is the initial temperature and reaches its maximum value, T_min is the minimum temperature, and θ∈ [0,1] is the annealing factor; during the iteration process, the change of the Q value table is considered, and the ε value is defined by the difference of the Q value and the dynamic temperature value; ε=exp((-( Δ Q_g)) / T_g) formula 8; Where ΔQ_g is the increment of Q in the g-th generation, and T_g is the temperature value in the g-th generation.

5. A distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that: In sub-step A3 of step 2, the parameters of the three regions are updated according to different formulas, specifically including sub-steps D1-D4: Step D1: The parameters in the mutation strategy for the enhanced region are updated according to Equation 9_14; F_ (i,1) ∈randn_i (μ_ (F_1 ) ,0.1) Cr_ (i,1) ∈randn_i (μ_ (Cr_1 ) ,0.1) Formula 9; F_ (i+1,1) '=c*F_ (i+1,1)+(1-c)*F_ (i,1) Cr_ (i+1,1) '=c*Cr_ (i+1,1)+(1-c)*Cr_ (i,1) Formula 10; c = 0.9 - 0.9 * 10^(-5) * g / G_max (Formula 11); μ_ (F_1)=(1-c)*μ_ (F_1)+c*mean_A (ω_i*F_ (i,1)) Formula 12; μ_ (Cr_1)=(1_c)*μ_ (Cr_1)+c*mean_A (ω_i*Cr_ (i,1)) Formula 13; ; Where F_i and Cr_i are the mutation factor and crossover factor, μ_F is the mean of F_i, μ_Cr is the mean of Cr_i, fit_i is the fitness value of the i-th individual, Δfit_i is the change in fitness value, and mean_A represents the arithmetic mean; F_i and Cr_i both follow a Gaussian distribution in this region; the change in fitness value is used as the weighting coefficient, and the arithmetic mean is used for the calculation of the mean; Step D2: The parameters in the mutation strategy of the stable region are updated according to formula 15-19; F_ (i,2) ∈randc_i (μ_ (F_2 ) ,0.1) Cr_ (i,2) ∈randn_i (μ_ (Cr_2 ) ,0.1)) Formula 15; F_ (i+1,2) '=c*F_ (i+1,2)+(1-c)*F_ (i,2) Cr_ (i+1,2) '=c*Cr_ (i+1,2)+(1-c)*Cr_ (i,2) Formula 16; μ_ (F_2)=(1-c)*μ_ (F_2)+c*mean_L (ω_i*F_ (i,2)) Formula 17; μ_ (Cr_2)=(1-c)*μ_ (Cr_2)+c*mean_A (ω_i*Cr_ (i,2)) Formula 18; mean_L ( ω_i*F_ (i,2))=( ∑ ω_i*F_ (i,2)^2) / ( ∑ ω_i*F_ (i,2)) Formula 19; Where F_i follows a Cauchy distribution and its mean is the Lymer mean mean_L ( ω_i*F_ (i,2)), and Cr_i follows a Gaussian distribution and its arithmetic mean is used; Step D3: The parameters in the mutation strategy of the decay region are updated according to formulas 20-25; F_ (i,3) ∈randc_i (μ_ (F_3 ) ,0.1) Cr_ (i,3) ∈randc_i (μ_ (Cr_3 ) ,0.1) Formula 20; F_ (i+1,3) '=c*F_ (i+1,3)+(1-c)*F_ (i,3) Cr_ (i+1,3) '=c*Cr_ (i+1,3)+(1-c)*Cr_ (i,3) Formula 21; μ_ (F_3)=(1-c)*μ_ (F_3)+c*mean_L (ω_i*F_ (i,3)) Formula 22; μ_ (Cr_3)=(1-c)*μ_ (Cr_3 )+c*mean_l (ω_i*Cr_ (i,3)) Formula 23; mean_L (ω_i*F_ (i,3))=(∑ω_i*F_ (i,3)^2) / (∑ω_i*F_ (i,3) ) Formula 24; mean_L (ω_i*Cr_ (i,3))=(∑ω_i*Cr_ (i,3)^2) / (∑ω_i*Cr_ (i,3)Formula 25; Both F_i and Cr_i follow a Cauchy distribution, and the mean is calculated using the Lemo mean, which yields a larger value. Step D4: In each region, the average fitness φ of individuals per generation and the average fitness φ' after mutation and crossover are calculated according to formulas 26-27 respectively; φ=(∑f(X_i)) / NP (Formula 26) φ '=( ∑f(U_i)) / NP formula 27; φ'<φ Formula 28; Where U_i is the individual vector of the i-th individual after mutation and crossover; If formula 28 is satisfied, the generation has successfully mutated and the parameter self-feedback is executed; otherwise, formulas 9, 15, and 20 are executed. The self-feedback design for the three regions is shown in formulas 10, 16, and 21.

6. A distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that: The sub-step A5 of step 2, which improves the inferior solution, specifically includes sub-steps E1-E6: Step E1: Improve the individuals with small Q values ​​and those that are too densely packed in the region in the reinforcement learning mechanism, and use them in the population evolution process; The density of individuals is measured by the average distance between an individual and the centroid of a region, as well as the average fitness value. The centroid of a region is obtained according to Formula 29. m_c=(∑X_i ) / NP Formula 29; Step E2: The average distance from each individual to its respective region's centroid is calculated using formula 30; d_=(∑║X_i-m_c║) / NP (Formula 30); Step E3: The average fitness increment is defined according to Formula 31; Δf=(∑((f(X_i )-f(m_c))) / NP Formula 31; Step E4: The density DD is calculated according to formula 32; DD = Δf * (1 + λ / d - ) Formula 32; Where Δf represents the change in fitness value; Step E5: The value of λ is calculated according to formula 33; λ=(∑((U_d-L_d) / NP)^2) ^1 / 2Formula 33; Where D represents the dimension, L_d represents the lower bound, and U_d represents the upper bound. Obviously, λ is adaptively adjusted with the search space. Step E6: Individuals with small Q values ​​and overly dense populations are defined as inferior solutions, denoted by I_i, and the repair strategy is Equation 34: I_i=I_i+ζ*(m_c-I_i ) ,i ∈ {1,2,… ,NP} Formula 34; Where I_i is the i-th repaired solution, and ζ∈ [0,1] is used to control the degree to which the inferior solution is close to the centroid.

7. A distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that: In sub-step A6 of step 2, regarding the phenomenon that the fitness value does not update for several consecutive generations during the evolutionary process, the specific methods for regenerating individuals include sub-steps F1-F3: Step F1: Use the Cholesky decomposition method to generate a low-order triangular matrix M from the covariance matrix C; Formula 35: C=M*M^T; Step F2: Sample a 1×D matrix R using the standard normal distribution N(0,1); Step F3: Generate a new individual according to formula 36; X_new=μ+η(g)*M*R Formula 36; Where η(g) is the scaling factor, η(g)>1, and η(g) is updated according to formula 37-40; η(g)=η_min+((G_max-g) / (G_max))^β*(η_max-η_min) Formula 37; η_max=(U_d-L_d)*0.01 Formula 38; η_min = η_max / 10 (Formula 39) β = 1 - e^(g / (0.15*G_max)) Formula 40; Where η_max and η_min represent the maximum and minimum allowed values, determined according to the size of the search space, with L_d representing the lower bound and U_d representing the upper bound, G_max being the maximum number of iterations, and β being the correction parameter for nonlinear adjustment of η.

8. A distributed congested flow shop scheduling optimization system based on co-evolutionary algorithm according to claim 1, characterized in that: In sub-step A7 of step 2, the local search is performed according to formulas 41-42; X_i (g+1)=N(μ_i,δ_i)Formula 41 Where N(μ_i,δ_i) is a Gaussian distribution with mean μ_i and variance δ_i; δ_i=τ*║X_best-X_i (g) ║Formula 42 Where τ is the scaling factor.

Images