Parallel batch processing machine scheduling method and system based on adaptive differential evolution algorithm

Abstract

The present invention belongs to the field of production scheduling, and specifically discloses a parallel batch processing machine scheduling method and system based on an adaptive differential evolution algorithm, which includes: S1 generates an initial population, and sets the number of iterations t=1; S2 randomly selects a target individual from the population; S3 selects the selection probability of a mutation operator, a scaling factor, and a crossover probability interval with equal probability or adaptively according to historical data according to whether the number of iterations t is within a preset learning period; S4 performs differential mutation and crossover on the target individual to obtain a test individual; S5 updates the individuals in the population based on the objective function value of the test individual; S6 repeats S2 to S5 until all individuals in the population are selected, and updates an external archive set based on the individuals in the population; S7 determines whether t reaches the maximum number of iterations: if so, the iteration ends, otherwise t=t+1 is set, and returns to S2. The present invention has the characteristics of quickly solving the multi-objective scheduling problem of parallel batch processing machines, improving machine utilization and reducing energy consumption.

Description

Technical Field

[0001] The present invention belongs to the field of production scheduling, and more specifically, relates to a parallel batch processing machine scheduling method and system based on an adaptive differential evolution algorithm. Background Art

[0002] Unlike traditional processing machines, batch processing machines can process multiple workpieces at the same time, which improves the production efficiency of enterprises. The decision-making of batch processing machine scheduling problems includes how to assign workpieces to machines, how to group batches, how to sort batches, etc., which is an NP-hard problem. Batch processing machine scheduling problems are widely present in the fields of semiconductor manufacturing, steel mills, 3D printing, etc. In the multi-variety, small-batch customized production environment of intelligent manufacturing, facing the characteristics of short delivery time, studying batch processing machine scheduling problems can effectively improve the production efficiency and competitiveness of enterprises, which has important theoretical and practical significance. The existing methods for batch processing machine scheduling problems mainly focus on single objectives, and often assume that workpieces have the same delivery time. In actual production, the JIT (Just-In-Time) production mode emphasizes the on-time completion of workpieces. Early completion and delayed completion lead to increased inventory costs and delay penalty costs for enterprises, respectively. At the same time, studying multi-objective optimization scheduling with integrated energy consumption has become one of the hot topics.

[0003] Therefore, traditional research on scheduling problems for single objectives and single batch processors can no longer meet the actual needs of production. It is urgent to study intelligent optimization algorithms for multi-objective scheduling problems of parallel batch processors, with the goal of simultaneously minimizing the maximum completion time, lead time / tardiness time and total energy consumption, thereby improving the company's production efficiency and on-time rate and reducing production costs and energy consumption.

[0004] Differential evolution (DE) algorithm is a population-based global optimization algorithm, which is widely regarded as one of the most commonly used algorithms for solving complex optimization problems. By combining differential mutation, crossover and selection operations in each generation, the population is led to the global optimum. It has the advantages of simplicity and high efficiency, and is a hot topic of concern and research by scholars in different fields. The mutation strategy and control parameters of traditional DE algorithms are often specified by users and remain unchanged throughout the evolution process, which is not conducive to balancing the exploration and development capabilities of the algorithm. Summary of the invention

[0005] In view of the above defects or improvement needs of the prior art, the present invention provides a parallel batch processor scheduling method and system based on an adaptive differential evolution algorithm, the purpose of which is to improve the differential evolution algorithm for the multi-objective scheduling problem of parallel batch processors and achieve a rapid solution to the multi-objective scheduling problem of parallel batch processors.

[0006] To achieve the above object, according to one aspect of the present invention, a parallel batch processor scheduling method based on an adaptive differential evolution algorithm is proposed, comprising the following steps:

[0007] Determine the decision variables and objective functions of the multi-objective scheduling problem of parallel batch processing machines and build a multi-objective scheduling model; use an adaptive differential evolution algorithm to solve the multi-objective scheduling model and obtain a workpiece production scheduling plan;

[0008] When using the adaptive differential evolution algorithm to solve the multi-objective scheduling model, the following steps are included:

[0009] S1. Generate an initial population, where each individual in the population represents a workpiece production scheduling plan; select some individuals in the initial population as the initial external archive set, and set the number of iterations t = 1;

[0010] S2, randomly select a target individual from the population;

[0011] S3, determine whether the number of iterations t is within the preset learning period lp: if t≤lp, then select each mutation operator, scaling factor interval and crossover probability interval with equal probability; otherwise, determine the selection probability of each mutation operator, scaling factor interval and crossover probability interval according to the performance of each mutation operator, scaling factor interval and crossover probability interval in the previous iteration;

[0012] S4, determining the mutation operator, scaling factor and crossover probability according to the selection probability of S3, and then performing differential mutation and crossover on the target individual to obtain the test individual;

[0013] S5. Based on the objective function value of the test individual, update the individuals in the population;

[0014] S6, repeating steps S2 to S5 until all individuals in the population of this generation have been selected at least once; then updating the external archive set based on the individuals in the current population;

[0015] S7. Determine whether the number of iterations t reaches the preset maximum number of iterations: if so, the iteration ends and the individual in the external archive set at this time is taken as the final solution; otherwise, set t=t+1 and return to step S2.

[0016] As further preferred, in step S5, based on the objective function value of the test individual, the individuals in the population are updated, comprising the following steps:

[0017] Decode the test individuals, calculate the corresponding objective function value, and update the reference point according to the objective function value; based on the reference point, use the Chebyshev aggregation function to select and update the individuals in the population;

[0018] Among them, the decoding method for individuals is: using heuristic rules to group and sort the workpieces assigned to the machine, and then using RLS rules to optimize the batch start time;

[0019] Reference point Z * =min{z j ,j=1,…,d},z j represents the global minimum of the jth objective function in the current population and test individuals, and d is the number of targets.

[0020] As a further preferred embodiment, the heuristic rule adopted is the IUDERT rule, that is, the workpieces assigned to the machine are grouped into batches according to the IUD rule, and the processing batches on the machine are sorted according to the ERT rule; in the batching process, the batch completion time, batch arrival time and batch delivery time are considered at the same time.

[0021] As a further preferred method, the batch start time is optimized by using the RLS rule, including the following steps:

[0022] S521, generate a random variable q=random{0,1}; initially, set the batch index k=1;

[0023] S522, if q=1, go to step S523; otherwise, the kth batch B k The start time remains unchanged, set k=k+1, and go to step S521;

[0024] S523, calculate batch B k Processing time P k , arrival time R k , and batch B k The median of the delivery time of the workpiece

[0025] S524, if k=1, set batch B k Completion time CT k =R k +P k ; Otherwise, set CT k =max{R k ,CT k-1}+P k ;

[0026] S525, if Then move batch B to the right k Make its completion time equal to like Then continue to judge as follows:

[0027] If R k ≤CT k-1 , define continuous np The processing batch is Block p, and the set of workpieces in Block p is A; calculate the parameter value of Block p Among them, the collection Indicates that Block p satisfies c j >d j Conditional artifact set, collection Indicates that Block p satisfies c j =d j Conditional artifact set, collection Indicates that Block p satisfies c j <d j The set of artifacts of the condition; c j represents the completion time of the workpiece, d j Indicates the delivery time of the workpiece;

[0028] If Δ p <0: Block p does not move;

[0029] If Δ p ≥0: (1) Initially, k = 1. When batch B k The start time is ST k ≥R k ,CT k-1 ≤R k When , Block p is shifted left by one time unit and ST is updated k =ST k -1, CT k =CT k -1; otherwise, Block p does not move; (2) Set k = k + 1 and return to step (1) until k = n p ;

[0030] If R k >CT k-1 , then no movement is performed;

[0031] S526, set k=k+1, and return to step S521 until For batch B k The number of workpieces in .

[0032] As a further preferred method, based on the reference point, the individuals in the population are updated using the Chebyshev aggregation function, including the following steps:

[0033] Randomly select an individual Y from the current population and calculate the test individual according to the following formula and the Chebyshev aggregation function value g of individual Y tche (X|λ,z):

[0034]

[0035] Where: z j represents the global minimum of the jth objective function, represents the global maximum value of the jth objective function; λ j is the weight vector of the jth target; f j is the jth objective function value of the individual;

[0036] According to the Chebyshev aggregation function value, the experimental individuals Compared with individual Y: If is better than individual Y, then use Replace individual Y, otherwise, keep individual Y.

[0037] As a further preferred embodiment, in step S3, the selection probability of each mutation operator, scaling factor interval and crossover probability interval is determined according to the performance of each mutation operator, scaling factor interval and crossover probability interval in the previous iteration; including:

[0038] The selection probability pm of the oth mutation operator o,t The calculation is as follows:

[0039]

[0040]

[0041] in, represents the total number of test individuals generated by the mutation operator o during the iteration from t-lp to t-1 population. It represents the total number of experimental individuals generated by mutation operator o in the iteration process of population from t-lp to t-1 generations that successfully enter the next generation; O is the total number of candidate mutation operators; is a positive number less than 0.01;

[0042] The probability of selecting the vth scaling factor interval pf v,t The calculation is as follows:

[0043]

[0044]

[0045] in, represents the sum of the experimental individuals generated by the candidate scaling factor interval v during the learning period [t-lp,t-1], It represents the total number of experimental individuals generated by the candidate scaling factor interval v during the iteration process from t-lp to t-1 generations that successfully enter the next generation; V is the total number of candidate scaling factor intervals;

[0046] The selection probability pc of the vth crossover probability interval v,t The calculation is as follows:

[0047]

[0048]

[0049] Among them, the candidate crossover probability interval is the same as the scaling factor interval; represents the sum of the experimental individuals generated in the candidate crossover probability interval v during the [t-lp, t-1] learning period, It represents the total number of experimental individuals generated during the iteration of the candidate crossover probability interval v from t-lp to t-1 generations that successfully enter the next generation.

[0050] As a further preferred embodiment, in step S6, updating the external archive set based on individuals in the current population includes:

[0051] S61, the current t generation population N POP individuals, and N in the external archive set of generation t-1 EP Individuals merge;

[0052] S62, use fast non-dominated sorting to sort N POP +N EP Individuals were ranked hierarchically;

[0053] S63, put the individuals into the external archive set layer by layer; when all the individuals in a certain layer are put in, the number of individuals in the external archive set is greater than N EP When the HV index value is used to sort the individuals in this layer, the better individuals are selected in order to keep the number of individuals in the external archive set to N. EP .

[0054] As a further preferred method, when the adaptive differential evolution algorithm is used for solving, a mixture uniform design method is used in advance to generate N with the same number of individuals as the population. POP Initial weight vector set Each weight vector corresponds to an individual; calculate each weight vector λ i With other weight vector λ k The Euclidean distance between them is used to determine the weight variables λ i T neighborhood vectors of P ;

[0055] When performing differential mutation and crossover on the target individual in step S4, individuals are selected from the current population or the neighborhood corresponding to the target individual according to a random rule, and differential mutation is performed on the target individual based on the selected individuals.

[0056] As further preferred, the objective function includes minimizing maximum completion time, lead time / delay time, and total energy consumption.

[0057] According to another aspect of the present invention, a parallel batch processor scheduling system based on an adaptive differential evolution algorithm is provided, comprising a processor, wherein the processor is used to execute the above-mentioned parallel batch processor scheduling method based on an adaptive differential evolution algorithm.

[0058] In general, the above technical solution conceived by the present invention has the following technical advantages compared with the prior art:

[0059] 1. Aiming at the multi-objective scheduling problem of parallel batch processing machines, the present invention makes an adaptive improvement to the differential evolution algorithm. Specifically, due to the influence of mutation operators and control parameters on the random optimization of the differential evolution algorithm, the present invention sets a learning period and uses historical data to adaptively select mutation operators and control parameters, thereby avoiding the differential evolution algorithm from falling into local convergence and generating a better scheduling plan quickly and efficiently.

[0060] 2. The present invention adopts the IUDERT rule to group the workpieces, and specially considers the arrival time of the workpiece batches; at the same time, on the basis of grouping and batch sorting, a batch start time optimization method based on the RLS rule is proposed to improve the utilization rate of the machine and the on-time rate of the order, and reduce the energy consumption.

[0061] 3. Aiming at multi-objective problems, the present invention introduces the decomposition-based multi-objective optimization idea into the differential evolution algorithm, proposes a reference point, and uses an improved Chebyshev aggregation function based on the reference point to select solutions, thereby balancing the exploration and development capabilities of the multi-objective differential evolution algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is a flow chart of a parallel batch processor scheduling method based on an adaptive differential evolution algorithm according to an embodiment of the present invention;

[0063] Figure 2 This is a schematic diagram of an encoding method according to an embodiment of the present invention;

[0064] Figure 3 A schematic diagram of a variation method of an embodiment of the present invention;

[0065] Figure 4 A schematic diagram of a crossover method according to an embodiment of the present invention;

[0066] Figure 5 (a)-(c) are the target values ​​C of the MODE / D algorithm according to the embodiment of the present invention. max , the main effect value of parameter signal-to-noise ratio S / N when ET and TEC are respectively used as response columns;

[0067] Figure 6 (a)-(c) are the target values ​​C of the NSGA-II algorithm according to the embodiment of the present invention. max , the main effect value of parameter signal-to-noise ratio S / N when ET and TEC are respectively used as response columns;

[0068] Figure 7 (a)-(c) are the target values ​​C of the NSGA-III algorithm according to the embodiment of the present invention. max , the main effect value of parameter signal-to-noise ratio S / N when ET and TEC are respectively used as response columns;

[0069] Figure 8 (a)-(c) are the target values ​​C of the IUDERT-RLS decoding method according to the embodiment of the present invention. max , the main effect value of parameter signal-to-noise ratio S / N when ET and TEC are respectively used as response columns;

[0070] Fig. 9 (a)-(c) are box plots of C matrix values, NR values ​​and IGD values ​​decoded by FFERT-RLS and IUDERT-RLS of the four algorithms in the embodiments of the present invention, respectively. DETAILED DESCRIPTION

[0071] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0072] An embodiment of the present invention provides a parallel batch processor scheduling method based on an adaptive differential evolution algorithm.

[0073] The multi-objective scheduling problem of parallel batch processing machines can be described as assigning n different jobs J = {J1, J2, …, J n}arranged on m parallel batch processing machines M={M1,M2,…,M m}, in this embodiment, the maximum completion time C is minimized max , the lead time / tardiness (Earliness and Tardiness, ET) and total energy consumption (Total Energy Consumption, TEC) are the targets. j With different processing times p j , arrival time r j 、Sizes j and delivery date j . Different batch machines M iWith different capacity constraints Q i , different energy consumption per unit time l i The decision-making content includes allocating n workpieces to m machines, determining the batching method of the workpieces allocated to each machine, and the processing order and start time of the batch. The basic assumptions satisfied are: (1) The sum of the sizes of all workpieces in a batch cannot exceed the capacity constraint Q of the allocated batch processing machine i ; (2) When batch processor M i Processing a batch B ib The process cannot be interrupted until all workpieces in the batch are processed; (3) a batch B b The arrival time R b The processing time of a batch is determined by the arrival time of the latest workpiece in the batch. b =Equal to the processing time of the workpiece with the longest processing time in the batch; (4) The batch arrival time R is not met b All workpieces in the batch are not allowed to be processed.

[0074] Aiming at the multi-objective scheduling problem of the parallel batch processor mentioned above, a multi-objective scheduling model is constructed. The adaptive differential evolution algorithm (SMODE / D) is used to solve the multi-objective scheduling model. Figure 1 As shown, the following steps are included:

[0075] S1, generate the initial population and determine the neighborhood of individuals in the population; initialize the reference point and select some individuals in the initial population as the initial external archive set EP; set the number of iterations t = 1;

[0076] S2. Select target individuals from the population according to random rules

[0077] S3, determine whether the number of iterations t is within the preset learning period lp: if t≤lp, then select each mutation operator, scaling factor interval and crossover probability interval with equal probability; otherwise, determine the selection probability of each mutation operator, scaling factor interval and crossover probability interval based on historical data in previous iterations;

[0078] S4: Determine the mutation operator, scaling factor and crossover probability according to the selection probability of S3, and then perform differential mutation and crossover on the target individual to obtain the experimental individual.

[0079] S5. Decode the test individual and calculate the corresponding objective function value. Update the reference point Z according to the objective function value. * , and update the individuals in the population;

[0080] S6, repeating steps S2 to S5 until all individuals in the population of this generation have been selected at least once; then updating the external archive set based on the individuals in the current population;

[0081] S7. Determine whether the number of iterations t reaches the preset maximum number of iterations Maxt: if so, the iteration ends and the individual in the external archive set at this time is taken as the final solution; otherwise, set t=t+1 and return to step S2.

[0082] Furthermore, in step S1, generating an initial population includes:

[0083] S11. Randomly generate a sequence of n workpieces, and then randomly distribute the workpieces to m parallel batch processing machines to generate an initial population P0.

[0084] S12. For the obviously inferior solutions in the initial population, that is, the solutions that do not have workpieces assigned to a certain machine, use the following method to repair the solutions:

[0085] S121. Initialize the machine set for allocating workpieces The collection of machines that have no workpieces assigned to them

[0086] S122, determine whether i>m (i is the machine index, initially 1): if so, end; otherwise, go to the next step.

[0087] S123, Statistical Machine M i The number of workpieces n allocated i .

[0088] S124, if n i = 0, then the machine M i Put into the set UAJS, UAJS=UAJS∪{M i}; Otherwise, machine M i Put into the set AJS, AJS=AJS∪{M i}.

[0089] S125, if Select the energy consumption per unit time l in the set UAJS i The smallest machine M min The machine M with the most workpieces in the set AJS max .

[0090] S126, from machine M max Randomly select a max -1] is assigned to machine M min superior.

[0091] S127, update set UAJS = UAJS\M min, the set AJS=AJS∪M min .

[0092] S128, Update machine M min and M max The number of workpieces allocated, let i=i+1, and return to step S122.

[0093] S13, using heuristic rules to group and sort the workpieces assigned to the machine, and then using RLS rules to optimize the batch start time. The specific decoding method is described in steps S51 and S52, which will not be repeated here.

[0094] Furthermore, in step S1, the method for determining the neighborhood is as follows:

[0095] Mixture Uniform Design (MUD) is used to generate N with the same number of individuals as the population. POP The initial weight vector set Π={λ1,…,λ NP}. Then calculate each weight vector λ i (i=1,…,N P ) and other weight vectors λ k (k=1,…,N P ), select T vectors with a closer Euclidean distance as the weight variable λ i The neighborhood vector set B(i) = {λ [1] ,…,λ [T]}.

[0096] Furthermore, step S3 specifically includes:

[0097] (1) The steps for adaptively selecting mutation operators are as follows.

[0098] Set a set SET containing 4 mutation operators DE / rand / 1, DE / rand / 2, DE / best / 1, and DE / current-to-rand / 1.

[0099] When the number of iterations t≤lp, a mutation operator is selected from the set SET with an equal probability of 1 / 4 to perform a differential mutation operation. lp represents the length of the learning period, and in this embodiment, lp=0.3*G is set, and G represents the maximum number of iterations.

[0100] When t≥lp+1, the mutation operator is adaptively selected according to the probability of each mutation operator, and the probability of mutation operator o being selected is pm o,t It is calculated according to the following formula.

[0101]

[0102]

[0103] in, represents the total number of test individuals generated by the mutation operator o during the iteration from t-lp to t-1 population. It represents the total number of experimental individuals generated by the mutation operator o in the iteration process of the population from t-lp to t-1 generations that successfully enter the next generation. Set to 0.01 to prevent sm o,t is 0. The more experimental individuals generated by the mutation operator o successfully enter the next generation, the greater the probability of being selected.

[0104] In addition, due to the optimization of multiple objectives, the mutation operator DE / best / 1 Individuals cannot be optimally determined by the objective function value. Based on the Pareto non-dominated relationship, select from the external archive set EP Considering that there are also differences in convergence and distribution between individuals in the same dominance layer, the solution with the largest HV value is selected from EP as individual.

[0105] (2) The steps for adaptively selecting the shrinkage factor interval are as follows.

[0106] The set contains 4 candidate shrinkage factor intervals, namely F1: [0, 0.25), F2: [0.25, 0.5), F3: [0.5, 0.75), F4: [0.75, 1], forming the candidate set SETF.

[0107] When the number of iterations t≤lp, a candidate interval is randomly selected with equal probability in the candidate set, and the shrinkage factor F is generated to perform the differential mutation operation of the target individual.

[0108] When the number of iterations t≥lp+1, the algorithm adaptively selects candidate intervals to generate shrinkage factors based on the historical data accumulated in the previous lp generations. The probability pf of interval v being selected v,t It is calculated according to the following formula.

[0109]

[0110]

[0111] in, represents the total number of experimental individuals generated in the candidate interval v during the [t-lp, t-1] learning period, It represents the total number of experimental individuals generated in the candidate interval v during the iteration process of the population from t-lp to t-1 generations that successfully enter the next generation.

[0112] (3) The steps for adaptively selecting the crossover probability interval are as follows.

[0113] Set the candidate crossover probability interval to be the same as the shrinkage factor interval.

[0114] When the number of iterations t≤lp, a candidate interval is randomly selected with equal probability in the candidate set, and the crossover probability CR is generated to perform the crossover operation of the target individual.

[0115] When the number of iterations t≥lp+1, the algorithm adaptively selects candidate intervals to generate crossover probability based on the historical data accumulated in the previous lp generations. The probability pc of candidate interval v being selected v,t It is calculated according to the following formula.

[0116]

[0117]

[0118] in, represents the total number of experimental individuals generated in the candidate interval v during the [t-lp, t-1] learning period, It represents the total number of experimental individuals generated in the candidate interval v during the iteration process of the population from t-lp to t-1 generations that successfully enter the next generation.

[0119] Further, in step S4, after determining the mutation operator, scaling factor and crossover probability, differential mutation and crossover are performed on the target individual, including:

[0120] S41, such as Figure 3 As shown, the execution process of the differential mutation operator is as follows.

[0121] The differential mutation operation is performed on the artifact sorting layer. Here, the differential mutation process is explained by taking the mutation operator DE / rand / 1 as an example, as shown in the following formula.

[0122]

[0123] in, For variant individuals, represents different individuals (selected from the population or the neighborhood corresponding to the target individual according to random rules), and F represents the scaling factor.

[0124] Since the machine allocation layer of the workpiece is represented by integers, it is not possible to directly perform differential mutation operations, so the two-point mutation method is used. Randomly generate two random numbers p1 and p2 in [1,n], and p1≠p2, and exchange individuals The machine allocation genes at the corresponding positions realize the machine allocation layer variation of the workpiece.

[0125] S42, such as Figure 4 As shown, the execution process of the crossover operator is as follows.

[0126] Randomly generate n random numbers in the interval [0,1] to form a 1*n dimensional array randv. If the value of position j in the array randv is less than or equal to the crossover probability CR, the test individual The genetic inheritance variant individual V at this position i t The corresponding position gene; if it is greater than the crossover probability CR, the test individual The target individual for gene inheritance at this position The corresponding position gene. The crossover formula of the individual is as follows:

[0127]

[0128] Further, step S5 includes the following steps:

[0129] S51. Use heuristic rules to group and sort the workpieces assigned to the machine.

[0130] The workpieces assigned to the machine by each individual in the population can be batched and sorted according to the FFERT (First Fit Earliest Ready Time) or IUDERT (Improved Updated Distance Earliest Ready Time) rule. It is preferred to use the IUDERT rule, that is, use the IUD rule to batch the workpieces assigned to the machine, and sort the processing batches on the machine according to the ERT rule, including:

[0131] S511, initialize machine M i Unbatch artifact collection UFJS i , the number of workpiece batches Distance Matrix Minimum distance d min =0, batch index k=1, machine index i=1.

[0132] S512, determine whether i>m: if so, end; otherwise, go to the next step.

[0133] S513, for machine M i k = 1 to Positions, unbatch Workpieces are put into a batch.

[0134] S514. Calculate the distance d(j,w) between two processing batches in sequence according to the distance function of the following formula, and update the distance matrix D.

[0135]

[0136] Among them, α, β, γ are parameters; pj 、p w is the processing time of the workpiece or batch; r j 、r w is the arrival time of the workpiece or batch; d j d w is the delivery date of the workpiece or batch; j 、s w is the size of the workpiece or batch; Q i For machine M i Initially, one workpiece is a batch, and the above parameters refer to the workpiece; after batching, if there are more than two workpieces in a batch, the above parameters refer to the batch.

[0137] S515, let the upper triangular matrix of the distance matrix D be D′. If all the distance values ​​in D′ are ∞, stop batching; otherwise, find the minimum distance value d min (If there are multiple d min , select workpieces or batches with small indicators), and set d min Position d rc The two related batches B r and B c Combine into a new batch B b .

[0138] S516, update batch completion time P b =max{P r ,P c}, batch arrival time R b =max{R r ,R c}, batch delivery date n2 represents the new batch B b contains the number of workpieces; Go to step S513. Wherein, P r , P c Batch B r , B c Processing time; R r , R c Batch B r , B c Arrival time.

[0139] S517, Machine M i The batches built above i = {B1, B2, …, B nib}, according to the ERT rules, the machine M i Sort the processing batches on the . Let i = i + 1, and go to step S512.

[0140] S52. Use RLS rules to optimize the batch start time with the goal of minimizing ET, including:

[0141] S521, generate a random variable q=random{0,1}; initially, set the batch index k=1;

[0142] S522, if q=1, go to step S523; otherwise, the kth batch B k The start time remains unchanged, set k=k+1, and go to step S521;

[0143] S523, B k n k Sort the workpieces in ascending order of due date and calculate: P k =max{p ρ ,ρ∈B k}, Among them, P k , R k Batch B k Processing time and arrival time; ρ 、r ρ Batch B k Processing time and arrival time of workpieces; Indicates batch B k The median of the workpiece delivery time;

[0144] S524, if k=1, set batch B k Completion time CT k =R k +P k ; Otherwise, set

[0145] CT k =max{R k ,CT k-1}+P k ;

[0146] S525, if Then move batch B to the right k To complete the work like Then continue to judge as follows:

[0147] If R k ≤CT k-1 , define continuous n p The processing batch is Block p, and the set of workpieces in Block p is A; calculate the parameter value of Block p Among them, the collection Indicates that Block p satisfies c j >d jConditional artifact set, collection Indicates that Block p satisfies c j =d j Conditional artifact set, collection Indicates that Block p satisfies c j <d j The set of artifacts of the condition; c j represents the completion time of the workpiece, d j Indicates the delivery time of the workpiece;

[0148] If Δ p <0: Block p does not move;

[0149] If Δ p ≥0: (1) Initially, k = 1. When batch B k The start time is ST k ≥R k ,CT k-1 ≤R k When , Block p is shifted left by one time unit and ST is updated k =ST k -1, CT k =CT k -1; otherwise, Block p does not move; (2) Set k = k + 1 and return to step (1) until k = n p ;

[0150] If R k >CT k-1 , then no movement is performed;

[0151] S526, set k=k+1, and return to step S521 until For batch B k The number of workpieces in .

[0152] S53, after decoding is completed through S51 and S52, the corresponding objective function value is calculated, and the reference point is updated according to the objective function value; specifically including:

[0153] Reference point Z * =min{z j ,j=1,…,d},z j represents the global minimum value of the jth objective function in the current population and test individuals, d is the number of targets, and in this embodiment, d=3.

[0154] S54. Experimental individuals using Chebyshev aggregation function Compare the target value of the solution with that of the original parent generation individual and retain the better solution; specifically:

[0155] S541. For the test individuals Randomly select an individual Y from the current population or neighborhood (determine whether it is the neighborhood of the test individual or the entire population according to the random number); calculate the test individual according to the following formula and the Chebyshev aggregation function value g of individual Y tche (X|λ,z):

[0156]

[0157] Where: z j represents the global minimum of the jth objective function, represents the global maximum value of the jth objective function; λ j is the weight vector of the jth target; f j is the jth objective function value of the individual.

[0158] S542, according to the Chebyshev aggregation function value, the test individuals Compared with individual Y: If is better than individual Y, then use Replace individual Y, otherwise, keep individual Y.

[0159] Further, step S6, using fast non-dominated sorting, HV index and elite strategy to update the external archive set EP, specifically includes the following steps:

[0160] S61, the current t generation population N POP individuals, and N in the external archive set of generation t-1 EP Individuals merge;

[0161] S62, use fast non-dominated sorting to sort N POP +N EP Individuals were ranked hierarchically;

[0162] S63, put the individuals into the external archive set layer by layer; when all the individuals in a certain layer are put in, the number of individuals in the external archive set is greater than N EP When the HV index value is used to sort the individuals in this layer, the better individuals are selected in order to keep the number of individuals in the external archive set to N. EP , ensuring that the concentrated population in the external archives has a better distribution.

[0163] Here is a specific example:

[0164] In order to verify the effectiveness of the method proposed in this invention, a problem experimental environment is designed as shown in Table 1. The experimental environment includes the processing time p of the workpiece j , arrival time r j , workpiece delivery time d j and workpiece sizej The batch machine capacity constraint Q i and energy consumption per unit time l i Parameters, the value ranges of each parameter and the generation formula are shown in Table 1.

[0165] Table 1 Problem parameters of simulation experiments

[0166]

[0167] Set the number of small-scale workpieces n = 10, 20, and the number of machines m = 2, such as Figure 2 As shown in the figure, the number of medium-scale workpieces is n = 40, 60, 80, and the number of machines is m = 2, 3; the number of large-scale workpieces is n = 100, 150, 200, 250, and the number of machines is m = 3, 5. 16 groups of problem instances are generated. For each group of problem instances with different combinations of workpieces and machines, 10 groups of data are randomly generated for experiments, for a total of 160 groups of data.

[0168] Four groups of problems were selected from the 16 groups of problem examples, and the optimal values ​​of the performance parameters affecting the SMODE / D algorithm proposed in the present invention and the traditional NSGA-II, NSGA-III, and MODE / D algorithms were determined through Taguchi orthogonal experiments.

[0169] The main parameters that affect the performance of the MODE / D algorithm are: population number NP, scaling factor F and crossover probability CR. Through preliminary experiments, the value range of each parameter was determined to be NP = {50, 100, 150, 200}, F = {0.1, 0.3, 0.6, 0.9}, and CR = {0.1, 0.3, 0.6, 0.9}. Taking NP, F and CR as control factors, the orthogonal experimental combination is obtained as L 16 (4 3 ), use Minitab to generate an orthogonal experimental table, and set the target value C max ET and TEC were used as response columns for the signal-to-noise ratio experiment. The algorithm was run 10 times for each parameter combination in the orthogonal table, with a total of 480 experiments. The main effect values ​​of the signal-to-noise ratio of the three objectives of the MODE / D algorithm were obtained as follows: Figure 5 shown.

[0170] Since the optimal parameter combinations determined by the main effect values ​​of the signal-to-noise ratio obtained according to different targets as the response columns are not exactly the same, the normalization method is used to comprehensively determine the unique value of each parameter for the three targets.

[0171] For the SMODE / D algorithm, in order to make an effective comparison with the MODE / D algorithm, the population number NP is set to be consistent with that of the MODE / D algorithm.

[0172] For algorithms NSGA-II and NSGA-III, the parameters that affect algorithm performance are: population size NP, crossover probability p c and mutation probability p m The preliminary experiment determined that the range of each parameter is NP = {50, 100, 150, 200}, p m ={0.2,0.4,0.6,0.8}, p c ={0.2,0.4,0.6,0.8}. m and p c As the control factor, the main effects of the signal-to-noise ratio of algorithms NSGA-II and NSGA-III are obtained according to the same parameter adjustment process as the MODE / D algorithm. Figure 6 and Figure 7 The final parameter values ​​of the four algorithms are shown in Table 2.

[0173] Table 2 Experimental parameters of different algorithms

[0174]

[0175] The setting of coefficients α, β, and γ in the IUD distance formula will directly affect the batching of workpieces. In order to select the most appropriate coefficient combination, five levels are set for each coefficient, namely α = {0.2, 0.4, 0.6, 0.8, 1}, β = {0.2, 0.4, 0.6, 0.8, 1}, and γ = {0.2, 0.4, 0.6, 0.8, 1}. Taking 10 workpieces as an example, the coefficients α, β, and γ are used as parameters, and the number of parameter combinations in the orthogonal experimental table is L. 25 (5 3 ), the target value was used as the response column for the signal-to-noise ratio experiment. Each parameter combination algorithm in the orthogonal table was run 10 times, for a total of 250 experiments. The main effect value of the signal-to-noise ratio was obtained as follows Figure 8 According to the same normalization method as the algorithm parameters, the coefficients in the distance formula are α = 0.6, β = 0.2, and γ = 0.6.

[0176] The experiment was conducted using the experimental parameters shown in Table 2, and the C matrix, NR, and inverse generation distance IGD indicators were used to evaluate the performance of the algorithm. The larger the C matrix and NR values, the smaller the IGD value, and the better the performance of the algorithm.

[0177] When the batch grouping and batch sorting rules are FFERT, the C matrix values, NR and IGD values ​​obtained by the four algorithms are shown in Tables 3 and 4, respectively.

[0178] Table 3 C matrix values ​​obtained by different algorithms based on FFERT-RLS rule

[0179]

[0180] Table 4 NR and IGD values ​​obtained by different algorithms based on the FFERT-RLS rule

[0181]

[0182]

[0183] As can be seen from Table 3, the convergence of the four algorithms is ranked from worst to best as SMODE / D algorithm, MODE / D algorithm, NSGA-III algorithm, and NSGA-II algorithm. The reason why the convergence of NSGA-III algorithm is better than that of NSGA-II algorithm is that NSGA-II algorithm is good at processing low-dimensional optimization problems, while this embodiment optimizes three objectives at the same time, and the NSGA-III critical layer adopts the reference point method to select individuals, which can better maintain population diversity and convergence compared with NSGA-II. From the perspective of the C matrix value, the designed MODE / D algorithm and the improved SMODE / D algorithm are far superior to the NSGA-II algorithm and the NSGA-III algorithm as a whole. The MODE / D algorithm and the improved SMODE / D algorithm obtain the vast majority of non-dominated solutions for 75% of the instance problems. This is because the decomposition-based multi-objective evolutionary algorithm uses weight vectors to generate reference points, decomposes the multi-objective optimization problem into single-objective sub-problems, and can achieve better distribution and convergence effects on the Pareto frontier. In comparison, the C matrix value of the SMODE / D algorithm is better than that of the MODE / D algorithm. This is because SMODE / D adaptively selects the generation interval of the mutation operator and control parameters based on the historical data learning feedback of the search process, achieves the optimal adjustment of the mutation operator and control parameters, balances the search and development characteristics of the mutation operator and control parameters, and improves the algorithm's optimization ability.

[0184] The experimental results in Table 4 also fully demonstrate that the convergence and distribution of the SMODE / D algorithm are better than those of the other three algorithms. Among the 16 problem instances, the SMODE / D algorithm obtained the ideal Pareto * The vast majority of solutions in the solution set; the IGD value of the solution set obtained by the SMODE / D algorithm is much smaller than the IGD value of the solution set obtained by the NSGA-III algorithm and the NSGA-II algorithm, and is significantly smaller than the IGD value of the solution set obtained by the MODE / D algorithm.

[0185] When the batch grouping and batch sorting rules are FFERT, the C matrix values, NR and IGD values ​​obtained by the four algorithms are shown in Tables 5 and 6, respectively.

[0186] Table 5 C matrix values ​​obtained by different algorithms based on IUDERT-RLS rule

[0187]

[0188] Table 6 NR and IGD values ​​obtained by different algorithms based on the IUDERT-RLS rule

[0189]

[0190]

[0191] It can be seen from Tables 5 and 6 that the distribution and convergence of the solution set obtained by using the IUDERT-RLS rule for batch grouping and processing batch sorting are almost the same as the results of the FFERT-RLS rule compared with the other three algorithms. The SMODE / D algorithm of the present invention is better than the MODE / D algorithm, the MODE / D algorithm is better than the NSGA-III algorithm, and the NSGA-III algorithm is better than the NSGA-II algorithm. This is because in the process of solving the problem, the IUDERT-RLS and FFERT-RLS rules for the sorted workpieces assigned to the machine are only different in the batch grouping method, and the overall framework of the algorithm is the same, so the overall performance trend of the algorithm has not changed.

[0192] The Pareto solutions obtained by the four algorithms when the batch grouping and batch sorting rules are FFERT-RLS and IUDERT-RLS rules respectively are combined, and the ideal Pareto solution is obtained by non-dominated sorting. * The solution set is then used to calculate the C matrix value, NR value and IGD value. The results are shown in Tables 7 to 9 respectively.

[0193] Table 7 C matrix values ​​of different algorithms based on FFERT rule and IUDERT rule

[0194]

[0195] Table 8 NR values ​​of different algorithms based on FFERT rule and IUDERT rule

[0196]

[0197] Table 9 IGD values ​​of different algorithms based on FFERT rule and IUDERT rule

[0198]

[0199] It can be seen from Tables 7 to 9 that, compared with the FFERT-RLS batching rule, the C value, NR value, and IGD value obtained by the IUDERT-RLS rule batching and sorting algorithms have obvious advantages, that is, whether it is the convergence or distribution of the solution set, the IUDERT-RLS rule can find a better solution. This is because the IUDERT-RLS rule puts workpieces with similar processing time, arrival time, and delivery time in the same batch during the batching process, which can reduce the space waste within the batch caused by the large difference in the arrival time and processing time of the workpieces, which is conducive to minimizing the target value C max and TEC, and also reduce the difference in delivery time between different workpieces in a single batch. When the batch / block is moved left or right to the appropriate position, the difference in delivery time between different workpieces in a single batch is small, which is conducive to minimizing the target value ET.

[0200] like Fig. 9 As shown in the figure, the IGD values ​​obtained by algorithms NSGA-II, NSGA-III and MODE / D are within (1240.31, 2923.61), (1125.45, 2720.11), and (374.98, 606.36), respectively, while the IGD value obtained by the SMODE / D algorithm is within (0, 167.09), and the value range is much smaller than that of the comparison algorithm, indicating that the SMODE / D algorithm can obtain solutions with better convergence and distribution.

[0201] The box plots of the Pareto solutions obtained after the four algorithms were run 10 times using the FFERT-RLS batch method and the IUDERT-RLS batch method are shown in the figure below. Fig. 9 As shown in the figure, the IGD values ​​obtained by the four algorithms using the FFERT-RLS batch method are (1649.90, 3495.52), (1079.24, 3667.01), (1029.96, 2595.65), (893.59, 2810.41), and the IGD values ​​obtained by the IUDERT-RLS batch method are (0, 841.41), (0, 962.61), (0, 512.04), (0, 311.27), respectively. The results show that the range of IGD values ​​obtained by the IUDERT-RLS batch method of different algorithms is smaller than the IGD value obtained by FFERT-RLS, indicating that the IUDERT-RLS batch method can obtain solutions with better convergence and distribution.

[0202] It will be easily understood by those skilled in the art that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A parallel batch processor scheduling method based on adaptive differential evolution algorithm, characterized in that: The steps include: Determine the decision variables and objective functions of the multi-objective scheduling problem of parallel batch processing machines and build a multi-objective scheduling model; use an adaptive differential evolution algorithm to solve the multi-objective scheduling model and obtain a workpiece production scheduling plan; When using the adaptive differential evolution algorithm to solve the multi-objective scheduling model, the following steps are included: S1. Generate an initial population, where each individual in the population represents a workpiece production scheduling plan; select some individuals in the initial population as the initial external archive set, and set the number of iterations t = 1; S2, randomly select a target individual from the population; S3, determine whether the number of iterations t is within the preset learning period lp: if t≤lp, then select each mutation operator, scaling factor interval and crossover probability interval with equal probability; otherwise, determine the selection probability of each mutation operator, scaling factor interval and crossover probability interval according to the performance of each mutation operator, scaling factor interval and crossover probability interval in the previous iteration; S4, determining the mutation operator, scaling factor and crossover probability according to the selection probability of S3, and then performing differential mutation and crossover on the target individual to obtain the test individual; S5. Based on the objective function value of the test individual, update the individuals in the population; S6, repeating steps S2 to S5 until all individuals in the population of this generation have been selected at least once; then updating the external archive set based on the individuals in the current population; S7, judging whether the number of iterations t reaches the preset maximum number of iterations: if so, the iteration ends, and the individuals in the external archive set at this time are taken as the final solution; Otherwise, set t=t+1 and return to step S2.

2. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm according to claim 1, characterized in that: In step S5, based on the objective function value of the test individual, the individuals in the population are updated, including the following steps: Decode the test individuals, calculate the corresponding objective function value, and update the reference point according to the objective function value; based on the reference point, use the Chebyshev aggregation function to select and update the individuals in the population; Among them, the decoding method for individuals is: using heuristic rules to group and sort the workpieces assigned to the machine, and then using RLS rules to optimize the batch start time; Reference point Z * =min{z j ,j=1,…,d},z j It represents the global minimum of the jth objective function in the current population and test individuals, and d is the number of targets.

3. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm as claimed in claim 2, characterized in that: The heuristic rule adopted is the IUDERT rule, that is, the workpieces assigned to the machine are grouped into batches according to the IUD rule, and the processing batches on the machine are sorted according to the ERT rule; in the batch grouping process, the batch completion time, batch arrival time and batch delivery time are considered at the same time.

4. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm as claimed in claim 3, characterized in that: Using RLS rules to optimize batch start time includes the following steps: S521, generate a random variable q=random{0,1}; initially, set the batch index k=1; S522, if q=1, go to step S523; otherwise, the kth batch B k The start time remains unchanged, set k=k+1, and go to step S521; S523, calculate batch B k Processing time P k , arrival time R k , and batch B k The median of the delivery time of the workpiece S524, if k=1, set batch B k Completion time CT k =R k +P k ; Otherwise, set CT k =max{R k ,CT k-1 }+P k ; S525, if Then move batch B to the right k Make its completion time equal to like Then continue to judge as follows: If R k ≤CT k-1 , define continuous n p The processing batch is Block p, and the set of workpieces in Block p is A; calculate the parameter value of Block p Among them, the collection Indicates that Block p satisfies c j >d j Conditional artifact set, collection Indicates that Block p satisfies c j =d j Conditional artifact set, collection Indicates that Block p satisfies c j <d j The set of artifacts of the condition; c j represents the completion time of the workpiece, d j Indicates the delivery time of the workpiece; If Δ p <0: Block p does not move; If Δ p ≥0: (1) Initially, k = 1. When batch B k The start time of ST k ≥R k ,CT k-1 ≤R k When , Block p is shifted left by one time unit and ST is updated k =ST k -1, CT k =CT k -1; otherwise, Block p does not move; (2) Set k = k + 1 and return to step (1) until k = n p ; If R k >CT k-1 , then no movement is performed; S526, set k=k+1, and return to step S521 until For batch B k The number of workpieces in .

5. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm as claimed in claim 2, characterized in that: Based on the reference point, the Chebyshev aggregation function is used to update the individuals in the population, including the following steps: Randomly select an individual Y from the current population and calculate the test individual according to the following formula and the Chebyshev aggregation function value g of individual Y tche (X|λ,z): Where: z j represents the global minimum of the jth objective function, represents the global maximum value of the jth objective function; λ j is the weight vector of the jth target; f j is the jth objective function value of the individual; According to the Chebyshev aggregation function value, the experimental individuals Compared with individual Y: If is better than individual Y, then use Replace individual Y, otherwise, keep individual Y.

6. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm according to claim 1, characterized in that: In step S3, the selection probability of each mutation operator, scaling factor interval and crossover probability interval is determined according to the performance of each mutation operator, scaling factor interval and crossover probability interval in the previous iteration; including: The selection probability pm of the oth mutation operator o,t The calculation is as follows: in, represents the total number of test individuals generated by the mutation operator o during the iteration from t-lp to t-1 population. It represents the total number of experimental individuals generated by mutation operator o in the iteration process of population from t-lp to t-1 generations that successfully enter the next generation; O is the total number of candidate mutation operators; is a positive number less than 0.01; The probability of selecting the vth scaling factor interval pf v,t The calculation is as follows: in, represents the sum of the experimental individuals generated by the candidate scaling factor interval v during the learning period [t-lp,t-1], It represents the total number of experimental individuals generated by the candidate scaling factor interval v during the iteration process from t-lp to t-1 generations that successfully enter the next generation; V is the total number of candidate scaling factor intervals; The selection probability pc of the vth crossover probability interval v,t The calculation is as follows: Among them, the candidate crossover probability interval is the same as the scaling factor interval; represents the sum of the experimental individuals generated in the candidate crossover probability interval v during the [t-lp, t-1] learning period, It represents the total number of experimental individuals generated during the iteration of the candidate crossover probability interval v from t-lp to t-1 generations that successfully enter the next generation.

7. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm according to claim 1, characterized in that: In step S6, the external archive set is updated based on the individuals in the current population, including: S61, the current t generation population N POP individuals, and N in the external archive set of generation t-1 EP Individuals merge; S62, use fast non-dominated sorting to sort N POP +N EP Individuals were ranked hierarchically; S63, put the individuals into the external archive set layer by layer; when all the individuals in a certain layer are put in, the number of individuals in the external archive set is greater than N EP When the HV index value is used to sort the individuals in this layer, the better individuals are selected in order to keep the number of individuals in the external archive set to N. EP .

8. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm according to claim 1, characterized in that: When using the adaptive differential evolution algorithm to solve the problem, the mixture uniform design method is used in advance to generate N with the same number of individuals as the population. POP Initial weight vector set Each weight vector corresponds to an individual; calculate each weight vector λ i With other weight vector λ k The Euclidean distance between them is used to determine the weight variables λ i T neighborhood vectors of P ; When performing differential mutation and crossover on the target individual in step S4, individuals are selected from the current population or the neighborhood corresponding to the target individual according to a random rule, and differential mutation is performed on the target individual based on the selected individuals.

9. The parallel batch processor scheduling method based on the adaptive differential evolution algorithm according to any one of claims 1 to 8, characterized in that: The objective function includes minimizing the maximum completion time, lead time / tardiness time, and total energy consumption.

10. A parallel batch processing machine scheduling system based on adaptive differential evolution algorithm, characterized in that: It includes a processor, which is used to execute the parallel batch processor scheduling method based on the adaptive differential evolution algorithm as described in any one of claims 1 to 9.

Images