A multi-objective dynamic scheduling decision method and system

Abstract

The application discloses a multi-target dynamic production scheduling decision method and system, and belongs to the technical field of flexible job shop production scheduling. The method comprises the following steps: collecting order, equipment and material data, and constructing a flexible job shop production scheduling model with the maximum completion time, weighted delay and equipment utilization rate as targets; solving the model by using a double-file cooperative evolution non-dominated sorting genetic algorithm, embedding a robust buffer before a key path process during decoding; deciding a Pareto non-dominated solution set by using an entropy weight-approximate ideal solution sorting method; quantifying disturbance severity based on right shift propagation simulation of a disjunctive graph, matching right shift repair, local re-optimization or global rolling rearrangement of a critical path according to the disturbance severity, and returning disturbance samples to a knowledge base to update a threshold value and a work hour standard deviation. The system comprises data acquisition, optimization engine, decision, disturbance monitoring, rescheduling, knowledge base and visualization modules. The application improves the quality of a Pareto frontier, the scheme has disturbance absorption capacity, and rescheduling is adaptively executed according to the severity.

Description

Technical Field

[0001] This invention belongs to the field of flexible workshop production scheduling technology, specifically relating to a multi-objective dynamic production scheduling decision-making method and system. Background Technology

[0002] With the widespread application of mixed-flow production (MJSP) in discrete manufacturing industries such as equipment manufacturing, machining, and electronic assembly, flexible job shops have become the mainstream workshop form in modern manufacturing. The Flexible Job Shop Scheduling Problem (FJSP) involves multiple candidate machines for each process. Therefore, solving it requires selecting suitable machines for each process and determining the start time of that process on the selected machines. In engineering implementation, the FJSP typically optimizes multiple objectives, including maximum completion time, weighted delay, and equipment utilization. Maximum completion time reflects the overall production cycle length, weighted delay reflects the deviation from customer delivery commitments, and equipment utilization reflects the efficiency of equipment resource use. These three factors are interdependent in most real-world workshop scenarios. The FJSP has been proven to be NP-hard. For medium to large-scale instances, it is difficult to obtain the optimal solution within an acceptable engineering timeframe using exact solution methods. Therefore, evolutionary computation methods are typically used to approximate the Pareto front within a reasonable timeframe. In evolutionary methods for solving multi-objective production scheduling problems, non-dominated sorting genetic algorithms, with their elite retention mechanism and crowding distance evaluation mechanism, have become a classic method in the field of multi-objective evolution and are widely used in engineering scheduling, parameter optimization and other fields.

[0003] In actual production operations, various disturbances frequently occur, primarily categorized into three types: equipment failure, urgent order insertion, and processing time deviations. Once a disturbance occurs, the original production scheduling plan becomes infeasible or no longer optimal at the execution level, requiring a corresponding rescheduling mechanism to repair or optimize the original plan to restore production order. The rescheduling mechanism should employ different response costs based on the severity of the disturbance: minor disturbances should be repaired at a low cost, while major disturbances should be rescheduled at a high cost, ensuring timely response while avoiding unnecessary global rescheduling that could cause fluctuations in on-site execution.

[0004] Existing multi-objective flexible job shop scheduling methods suffer from several technical problems when facing real-world issues such as multi-objective conflicts, frequent disturbances, and high timeliness requirements for decision-making. In the solution phase, existing methods often employ standard non-dominated sorting genetic algorithms. However, the crowding distance evaluation mechanism used in these algorithms lacks effective support for the front edge in target spaces with three or more objectives, resulting in a sparse distribution of the Pareto approximate front in the edge region, with relatively insufficient volume envelope and dispersion uniformity. The scheduling schemes derived by existing methods lack embedded robustness to absorb random fluctuations in shop floor processing time. Even slight delays in critical process processing time accumulate and amplify along the critical path, ultimately delaying the overall production cycle, leading to low stability at the execution level. Furthermore, in the Pareto front decision-making stage, existing methods often rely on manual weighting of each objective based on experience, followed by weighted summation to select a scheme. The decision results fluctuate with the subjective preferences of the decision-makers, lacking objectivity and traceability. Existing methods for handling workshop disturbances often employ either defining response rules separately for each type or uniformly performing global rescheduling. The former results in high maintenance costs for disturbance response rules, while the latter overreacts to small disturbances but fails to adequately cover large disturbances, making it difficult to balance timely response with the reasonable cost of rescheduling. Furthermore, existing methods mostly rely on static parameter configuration, failing to automatically adjust thresholds and time statistics parameters based on accumulated workshop operational experience to better reflect the actual statistical characteristics of the workshop, thus hindering the continuous optimization of the methods' effectiveness over time.

[0005] In summary, designing a multi-objective dynamic scheduling decision-making method and system to solve the aforementioned problems in existing technologies is of great significance. Summary of the Invention

[0006] To address the problems existing in the background art, the present invention provides a multi-objective dynamic scheduling decision-making method, comprising the following steps:

[0007] S1: Collect and preprocess order, equipment, material, and real-time status data;

[0008] S2: Construct a multi-objective flexible workshop scheduling model;

[0009] S3: Solve the model using a dual-file co-evolutionary non-dominated sorting genetic algorithm. During decoding, a robust buffer is embedded before the critical path process to obtain the Pareto non-dominated solution set.

[0010] S4: The initial execution plan is obtained by using the entropy weight-approximation ideal solution sorting method on the Pareto non-dominated solution set;

[0011] S5: Perform right-shift propagation simulation based on the disjunction graph corresponding to the current production schedule, and quantify the severity of the disturbance;

[0012] S6: Perform right-shift repair, local re-optimization of the critical path, or global rolling rearrangement according to the severity of the disturbance to obtain an updated solution;

[0013] S7: Output the update scheme and decision report, and return the disturbance sample to the disturbance statistical sample knowledge base.

[0014] Furthermore, step S1 includes:

[0015] S11: Read the workpiece number, process number, set of processable equipment for each process, standard processing time for each process on each processable equipment, release time of each workpiece, delivery date and lateness penalty weight of each workpiece through the manufacturing execution system; read the current availability status, remaining processing time and expected release time of each equipment from the programmable logic controller through the equipment IoT interface; read the current inventory and available time window of each material through the material management system;

[0016] S12: Map the data collected in step S11 to the memory data structure of the production scheduling solver according to preset fields; fill missing standard processing time, release time, remaining processing time, and inventory quantity numerical fields with historical averages or averages of similar processes; perform integrity checks on missing workpiece number fields, process sequence fields, and processable equipment set fields, and mark orders that fail the integrity check as pending manual confirmation; perform unified unit conversion on each field, unifying time fields to minutes and quantity fields according to the corresponding material's unit of measurement; obtain a production scheduling input dataset including workpiece set, equipment set, process set, processing time matrix, and constraint parameter list, as input for subsequent steps.

[0017] Furthermore, the multi-objective flexible workshop scheduling model constructed in step S2 uses a binary decision variable regarding whether the j-th process of workpiece i is assigned to the k-th machine for processing. and the actual start time of this process. The decision variables are defined as follows: minimizing the maximum completion time, minimizing the weighted delay, and maximizing equipment utilization. The optimization objectives are: process sequence constraints, unique equipment allocation constraints for processes, mutual exclusion constraints for processing on the same equipment, workpiece release time constraints, equipment machinability constraints, and material availability constraints. A value of 1 indicates allocation, and a value of 0 indicates no allocation. Let i be the actual start time of the j-th process of workpiece i; i is the workpiece number; j is the process number within the same workpiece; k is the equipment number; the weighted delay is the sum of the products of the delay value of each workpiece whose actual completion time is later than its delivery time and the corresponding delay penalty weight; the equipment utilization rate is the average of the proportion of the effective processing time of each equipment to the available time of that equipment.

[0018] The encoding and evolutionary operations of the dual-file co-evolutionary non-dominated sorting genetic algorithm in step S3 include:

[0019] S31: Using a double-layered chromosome as the individual code, the first segment is the process sorting segment, the length of which is the sum of the number of processes for each workpiece. Each gene position is filled with the corresponding workpiece number, and the number of times the same workpiece number appears in this segment is equal to the number of processes for that workpiece. The second segment is the equipment allocation segment, the length of which is the same as the first segment. Each gene position is filled with the sequence number of the equipment selected by the process corresponding to that position in its set of processable equipment.

[0020] S32: A hybrid strategy is used to generate the initial population. Some individuals in the initial population are generated by the global selection initialization method, while the remaining individuals are generated in a completely random manner.

[0021] S33: For the parent individual's process sorting segment, priority processes are maintained with crossover; for the equipment allocation segment, positional uniform crossover is adopted. Mutation is performed at a preset mutation rate. The process sorting segment mutation adopts the method of randomly swapping the positions of two different workpiece numbers. The equipment allocation segment mutation is performed on the gene bit corresponding to the critical path process identified by the decoding in step S34: when there are no less than two machines in the set of processable machines, the gene bit is changed to the sequence number of the candidate machine in the set of processable machines whose standard processing time is less than that of the currently selected machine; otherwise, the equipment allocation segment mutation is not performed.

[0022] Furthermore, step S3 further includes:

[0023] S34: Based on the individual's two-layer chromosome decoding, a preliminary production scheduling plan is obtained, and a directed acyclic disjunctive graph corresponding to the preliminary production scheduling plan is constructed; the earliest start time and earliest finish time of each process are calculated by forward traversal from the source node to the sink node in topological order; then, with the global maximum finish time of the preliminary production scheduling plan as the endpoint, the latest start time and latest finish time of each process are calculated by reverse traversal from the sink node to the source node in topological reverse order; the process with the difference between the latest start time and the earliest start time being zero is identified as the critical path process, and the corresponding process set is denoted as CP;

[0024] S35: For each critical process in the process set CP The key process is read from the disturbance statistical sample knowledge base. Standard deviation of processing time The key process is calculated according to the following formula. The robust buffer length that should be embedded :

[0025] ;

[0026] In the formula: Let be the robust buffer length that should be embedded before the j-th process of workpiece i when it is processed on equipment k. It is a robust multiple, dimensionless, and a positive real number; Let be the standard deviation of the processing time of the j-th operation on workpiece i on equipment k; i is the workpiece number; j is the operation number within the same workpiece; k is the equipment number. The j-th process of workpiece i is assigned to the process node processed by equipment k in the preliminary production scheduling plan;

[0027] The obtained robust buffer length Before the start time of the corresponding process is inserted as an additional delay, the earliest start time and earliest completion time of the direct successor node of the same workpiece, the direct successor node of the same equipment, and the recursive successor node of the aforementioned successor node are updated along the disjunction graph. The updated scheme is used as the input of step S34 again, and steps S34 and S35 are repeated until the critical path process set obtained by two consecutive iterations is the same or the number of iterations reaches the preset upper limit. The feasible production scheduling scheme corresponding to this individual is then output.

[0028] Furthermore, step S3 further includes:

[0029] S36: During the selection phase of each generation, the set obtained by merging the main evolutionary population with the current generation's offspring is divided into several non-dominant levels according to the non-dominant ranking. The non-dominant levels are then filled into the convergence file in order of increasing non-dominant level and decreasing crowding distance within the same non-dominant level. Until the converged file is reached. Fill the preset capacity;

[0030] S37: Normalize the objective value of each non-dominated solution in the merged set using the objective range to obtain a normalized objective vector; when the files are distributed... When the set is empty, the non-dominated solution with the largest Euclidean distance between the normalized target vector and the current ideal point in the merged set is taken as the dispersed file. The initial solution; calculate the non-dominated solutions in the merged set according to the following formula. Relative to the scattered archives Riesz2-Energy Contribution Value :

[0031] ;

[0032] In the formula: Non-dominated solution Relative to the scattered archives The Riesz2- energy contribution value is dimensionless. Non-dominated solution With the aforementioned scattered archives China-Africa Domination Solution The Euclidean distance over the normalized target vector is dimensionless. It is a preset minimum positive real number used to prevent division by zero when the distance is zero; Indicates taking and The larger of the two; For the distributed archive; i is the number of the non-dominated solution to be evaluated; l is the distributed archive. The numbers of the non-dominated solutions already included in the database; This indicates the location of the scattered archives. Unless the solution is dominated Summing the remaining non-dominated solutions;

[0033] according to The remaining non-dominated solutions in the merged set are filled into the dispersed file in ascending order. until the scattered archives Fill the preset capacity; every preset exchange interval algebra, from the converged file The non-dominated solution with the largest current congestion distance is selected for migration into the dispersed archive. and from the scattered archives The non-dominated solution with the smallest current non-dominated level value is selected and migrated into the convergence archive. The converged archive is then migrated. Non-dominated solutions from the scattered archives After removing from the middle, press From smallest to largest, select an undoed non-dominated solution from the merged set and add it to the scattered file. After the iteration terminates, the converged file will be... With the aforementioned scattered archives Take the union of the solutions, remove inferior solutions dominated by other solutions and duplicate solutions, and obtain the Pareto non-dominated solution set.

[0034] Furthermore, step S4 includes:

[0035] S41: Arrange the values ​​of each candidate solution in the Pareto non-dominated solution set for the three objectives of maximum completion time, weighted delay, and equipment utilization rate by row as candidate solutions and by column as objectives, to obtain an N-row, 3-column decision matrix, where N is the total number of candidate solutions in the Pareto non-dominated solution set; when N=1, directly issue the unique candidate solution in the Pareto non-dominated solution set as the initial execution solution to the workshop execution interface, and step S4 ends;

[0036] S42: Normalize the maximum completion time column and weighted delay column of the decision matrix according to cost-type indicators, and normalize the equipment utilization rate column according to benefit-type indicators to obtain a normalized decision matrix; when the range of a certain target column is zero, the normalized value of that column is uniformly set to a preset constant of 0.5; calculate the Shannon entropy of the q-th column of the normalized decision matrix according to the following formula. :

[0037] ;

[0038] In the formula: Let be the Shannon entropy of the q-th column of the normalized decision matrix, which is dimensionless; q is the target number, and q takes values ​​of 1, 2, and 3, which correspond to the three targets of maximum completion time, weighted delay, and equipment utilization rate, respectively; N is the total number of candidate solutions; The normalized value in the nth row and qth column of the normalized decision matrix is ​​the proportion of the sum of the normalized values ​​in the qth column, which is dimensionless; n is the candidate scheme number. The natural logarithm operator; This represents the summation over n=1 to n=N; when When it is zero, The item is treated as zero;

[0039] Then calculate the objective entropy weight corresponding to the q-th target using the following formula. :

[0040] ;

[0041] In the formula: Let be the objective entropy weight corresponding to the q-th target, which is dimensionless; For the summation index, Values ​​are 1, 2, and 3; For the normalized decision matrix, the first... The Shannon entropy of the column is dimensionless; Indicates to =1 to Summing = 3;

[0042] S43: Multiply each element in the q-th column of the normalized decision matrix by its corresponding objective entropy weight. The weighted normalized decision matrix is ​​obtained. From each column of the weighted normalized decision matrix, the maximum value is selected to form the positive ideal solution, and the minimum value is selected to form the negative ideal solution. The relative closeness of the nth candidate solution to the positive ideal solution is calculated using the following formula. :

[0043] ;

[0044] In the formula: Let be the relative closeness of the nth candidate solution to the positive ideal solution, which is dimensionless; Let be the Euclidean distance between the nth row vector of the weighted normalized decision matrix and the positive ideal solution, which is dimensionless; Let be the Euclidean distance between the nth row vector of the weighted normalized decision matrix and the negative ideal solution, which is dimensionless;

[0045] The relative closeness The candidate solution with the largest value is sent to the workshop execution interface as the initial execution solution.

[0046] Furthermore, step S5 includes:

[0047] S51: Continuously receive disturbance events e through the disturbance alarm interface, and perform one of the following processing on the disjunction diagram corresponding to the currently executing production scheduling plan according to the event type:

[0048] In the first category, when the disturbance event e is a device failure, the equivalent node of the dissociative graph node corresponding to the currently being processed or the next process to be processed on the faulty device is delayed as the expected fault repair time.

[0049] The second type is when the disturbance event e is an urgent order insertion type. The urgent order is added to the disjunction graph as a set of new workpiece nodes. The release time of the first process of the new workpiece is set as the urgent order insertion time. The process sequence constraints, equipment manufacturability constraints and material availability constraints of the new workpiece are added to the disjunction graph.

[0050] The third type is when the disturbance event e is a processing time deviation type. The equivalent node of the disjunction graph node corresponding to the process whose deviation exceeds the preset tolerance is delayed to the difference between the actual processing time and the standard processing time of the process.

[0051] S52: Using the processed disjunctive graph as input, traverse forward from the source node to the sink node in topological order. For each node: take the new start time of the node as the larger of the maximum value of the new completion times of all its direct predecessor nodes and the node's release time constraint; take the new completion time of the node as the sum of the new start time, standard processing time, and equivalent node delay; after traversal, obtain the updated start time and completion time of each process, and record the maximum value of the updated completion times as... ;

[0052] S53: Calculate the severity of the disturbance event e using the following formula. :

[0053] ;

[0054] In the formula: Let be the severity of the disturbance event e, which is dimensionless; The global maximum completion time obtained after the right-shift propagation simulation in step S52; The global maximum completion time of the current production schedule before the occurrence of the disturbance event e; This indicates taking the larger of a and b; the superscript RS indicates the state after the right-shift propagation simulation; the superscript plan indicates the planned state before the disturbance.

[0055] The severity of the disturbance This serves as the basis for matching the response strategy in step S6.

[0056] Furthermore, steps S6 and S7 include:

[0057] S61: Determine the severity of the disturbance With the slight threshold maintained in the perturbation statistics sample knowledge base and severe threshold Comparison: When Perform rightward repair when Perform local re-optimization of the critical path when A global scroll reordering is performed at that time; among which, The threshold value is light and dimensionless. The severity threshold is dimensionless. When the cumulative number of samples in the disturbance statistics sample knowledge base is less than a preset sample threshold, and Take the preset initial threshold from the disturbance statistics sample knowledge base;

[0058] S62: The right shift repair maintains the equipment allocation of each process in the current production schedule and the processing order of processes on the same equipment unchanged, and uses the new start time of each process obtained in step S52 to cover the original start time as the output of the updated plan;

[0059] S63: The critical path local re-optimization uses the right-shift propagation simulation scheme obtained in step S52 as the initial value, identifies its critical path, and unlocks the start time and equipment allocation of the processes on the critical path, the direct successor processes of the same workpiece on the critical path, the direct successor processes of the same equipment on the critical path, and the adjacent processes related to equipment reassignment as decision variables. The start time and equipment allocation of the remaining processes are fixed as determined parameters. The dual-file co-evolutionary non-dominated sorting genetic algorithm of step S3 is called to perform small-scale re-optimization on the reduced subproblem. During the re-optimization process, conflict verification is performed on the fixed process time windows on the reassigned equipment, and the insertionable time window constraint ensures that the processing of the same equipment is mutually exclusive. Then, the Pareto non-dominated solution set obtained by re-optimization is decided according to step S4 to obtain the updated scheme.

[0060] S64: The global rolling rescheduling fixes the processes that have started but not yet completed at the current time on the originally assigned equipment for continued processing. The start time of the process is fixed to the current time, and the remaining processing time of the process is fixed to the remaining processing time of the process. The processes that have not yet started are set as a new set of processes to be scheduled. The updated equipment status, material status and remaining orders are used as new scheduling inputs. Steps S2, S3 and S4 are re-executed to obtain the updated scheme.

[0061] S71: Output the update scheme in the form of a two-dimensional Gantt chart in the workpiece-equipment dimension, and attach a decision report containing three target values, critical path identifiers and buffer embedding location identifiers;

[0062] S72: The characteristic parameters of this disturbance event and the type of response strategy adopted are fed back to the disturbance statistical sample knowledge base as a sample. The characteristic parameters include the disturbance event type, the affected process number, the affected equipment number, the equivalent node delay value, and the disturbance severity. When the cumulative number of samples in the perturbation statistical sample knowledge base reaches the preset sample threshold, the perturbation statistical sample knowledge base updates the mild threshold based on the preset quantiles of the relative increment distribution of the maximum completion time actually generated by the right-shift repair strategy in the cumulative samples. and the severity threshold Based on the deviation between the actual processing time and the standard processing time of each process on different equipment in the cumulative sample, the corresponding standard deviation of the processing time is updated on a rolling basis according to the workpiece-process-equipment combination. .

[0063] This invention also provides a multi-objective dynamic scheduling decision-making system, including a data acquisition and preprocessing module, a multi-objective optimization engine module, an automatic solution decision-making module, a disturbance monitoring and impact assessment module, an adaptive rescheduling decision-making module, a disturbance statistical sample knowledge base module, and a solution output and visualization module, wherein:

[0064] The data acquisition and preprocessing module communicates with the workshop manufacturing execution system, enterprise resource management system, material management system, equipment IoT interface and programmable logic controller via industrial Ethernet;

[0065] The multi-objective optimization engine module includes a model building submodule, a key process robust decoding submodule, and a dual-file co-evolution submodule connected in sequence;

[0066] The output of the automatic decision-making module is connected to the first input of the visualization module and the workshop execution interface;

[0067] The disturbance monitoring and impact assessment module communicates with the disturbance alarm interface and programmable logic controller via industrial Ethernet;

[0068] The output of the adaptive rescheduling decision module is connected to the second input of the scheme output and visualization module, the feedback input of the multi-objective optimization engine module, and the return input of the perturbation statistical sample knowledge base module, respectively.

[0069] The perturbation statistics sample knowledge base module provides the processing time standard deviation to the robust decoding submodule of key processes, and provides mild and severe thresholds to the adaptive rescheduling decision module;

[0070] The output of the solution and visualization module is connected to the workshop visualization equipment and the operator's tablet terminal.

[0071] In the preferred embodiment, the output of the data acquisition and preprocessing module is connected to the input of the multi-objective optimization engine module via the system's internal data bus; the output of the multi-objective optimization engine module is connected to the input of the automatic scheme decision module via the system's internal data bus; and the output of the disturbance monitoring and impact assessment module is connected to the input of the adaptive rescheduling decision module via the system's internal data bus.

[0072] The critical process robust decoding submodule includes a critical path identification unit and a robust buffer embedding unit. The output of the critical path identification unit is connected to the first input of the robust buffer embedding unit, and the second input of the robust buffer embedding unit is connected to the output of the perturbation statistical sample knowledge base module. The critical path identification unit is used to identify the set of critical path processes by topological forward traversal and topological reverse traversal. The robust buffer embedding unit is used to embed a robust buffer before the critical path processes identified by the critical path identification unit.

[0073] The dual-archive co-evolution submodule includes a convergent archive storage unit, a distributed archive storage unit, a Riesz2-energy calculation unit, and an elite migration control unit. The output of the Riesz2-energy calculation unit is connected to the input of the distributed archive storage unit, and the elite migration control unit is bidirectionally connected to both the convergent archive storage unit and the distributed archive storage unit.

[0074] The disturbance monitoring and impact assessment module includes a disturbance event acquisition unit, a node delay mapping unit, a disjunction graph right shift propagation simulation unit, and a disturbance severity calculation unit connected in sequence. The other input of the disjunction graph right shift propagation simulation unit is connected to the output of the multi-objective optimization engine module to read the disjunction graph structure data of the currently executing production scheduling scheme.

[0075] The adaptive rescheduling decision module includes a threshold matching unit and a right-shift repair execution unit, a critical path local re-optimization execution unit, and a global rolling reordering execution unit that can be selectively connected to the output of the threshold matching unit. The other input of the threshold matching unit is connected to the output of the perturbation statistics sample knowledge base module to read the mild threshold and the severe threshold.

[0076] The beneficial effects achieved by this invention are as follows:

[0077] This invention introduces an adaptive robust buffer embedding method for critical processes during the decoding stage. For each critical process identified through critical path identification, the corresponding processing time standard deviation is read from the disturbance statistical sample knowledge base according to the workpiece-process-equipment combination. Then, the buffer length is embedded before the start time of the critical process by using the product of a preset robustness factor and the processing time standard deviation as the buffer length. The subsequent nodes affected by the buffer are updated in a chain along the disjunctive graph. Based on Chebyshev's inequality, the original production scheduling scheme has the ability to absorb first-order random disturbances within a controllable probability range under the objective condition that there are random fluctuations in workshop processing time. This avoids the phenomenon that the processing delay of a single critical process accumulates and amplifies along the critical path, ultimately lengthening the overall production cycle. It also allows for a flexible trade-off between the compactness and stability of the scheme by adjusting the robustness factor value during engineering implementation, making it easy to configure an appropriate robustness factor for different workshop disturbance characteristics.

[0078] This invention introduces a dual-archive co-evolutionary approach based on Riesz2-energy during the selection phase. By maintaining convergent archives selected according to non-dominated level and crowding distance, and dispersed archives selected according to Riesz2-energy contribution value, and performing elite migration between the two archives at a preset exchange interval for generations, the population of each generation achieves synergistic improvement in the two mutually constraining dimensions of convergence and dispersion. Addressing the front sparsity problem in traditional non-dominated sorting genetic algorithms in a three-objective space, where the crowding distance index loses its effective support for the front edge, this invention leverages the mathematical property of Riesz2-energy to simultaneously maintain central and edge support solutions in multi-objective fronts. This results in simultaneous improvements in the Pareto approximate front in two complementary dimensions: volume envelope and approximation accuracy. The algorithm exhibits stability in multiple independent runs, ensuring reproducibility of the solution during engineering implementation and facilitating the reproduction of the solution across different production batches and operators.

[0079] This invention introduces a combination of entropy weighting and a method for ranking near-ideal solutions in the decision-making stage. By automatically assigning objective weights to each objective based on Shannon entropy, and then selecting the final solution based on the weighted relative proximity, it avoids the subjective bias caused by assigning weights to multiple objectives based on human experience. The ability of each objective to distinguish the decision is automatically quantified into objective entropy weights according to its dispersion in the candidate solution set. Then, the relative proximity of each candidate solution to the positive and negative ideal solutions is calculated according to the objective entropy weights. The selected solution is the one that is closest to the positive ideal solution and furthest from the negative ideal solution. This makes the decision-making process repeatable and traceable. This invention can issue and execute solutions according to unified objective decision-making rules in different workshops and different production batches, which facilitates standardized deployment and production process auditing during engineering implementation.

[0080] This invention introduces a disturbance severity quantification based on right-shift propagation simulation of the disjunctive graph, a graded response based on disturbance severity, and a closed-loop maintenance mechanism for disturbance sample feedback during the dynamic execution phase. Different types of disturbance events, such as equipment failures, urgent order insertions, and processing time deviations, are uniformly mapped to a single dimensionless index. Based on the relative relationship between this index and mild and severe thresholds, it automatically selects one of three response strategies with different costs: right-shift repair, local re-optimization of the critical path, or global rolling rearrangement. This avoids the problems of overreaction or underreaction caused by applying the same response to every type of disturbance. Furthermore, by feeding disturbance samples and response records back to the disturbance statistical sample knowledge base to continuously update the mild threshold, severe threshold, and processing time standard deviation, the method of this invention continuously approximates the actual statistical characteristics of the workshop as workshop operating experience accumulates. This allows the invention to work with initial parameters provided by optimized values ​​in the early stages of project operation, and its graded response judgment and robust buffer embedding gradually become more accurate over time. The beneficial effects of the method continuously improve with the accumulation of workshop operating history. Attached Figure Description

[0081] Figure 1 This is a comparison chart of the hypervolume index convergence curves of Example 1 and Comparative Examples 1 to 4 on 10 standard benchmark examples from Brandimarte MK01 to MK10.

[0082] Figure 2 This is a comparison chart of the distribution of supervolume indexes obtained from 30 independent runs of Example 1 and Comparative Examples 1 to 4 on the BrandimarteMK02 benchmark instance.

[0083] Figure 3 The robustness multiple under Examples 1 to 4 and supplementary values. For the maximum completion time A comparison curve of the trade-off between the disturbance absorption rate of the proposed scheme and the actual perturbation absorption rate.

[0084] Figure 4 This is a comparison chart of the comprehensive performance of Example 1 and Comparative Examples 1, 3, and 4 on five standard benchmark examples from Brandimarte MK01 to MK05. (a) is a comparison chart of the supervolume index, and (b) is a comparison chart of the IGD index.

[0085] Figure 5 This is a flowchart of a multi-objective dynamic scheduling decision-making method according to the present invention. Detailed Implementation

[0086] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. In addition, the forms of the various structures described in the following embodiments are merely illustrative. The present invention is not limited to the structures described in the following embodiments. All other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0087] This embodiment addresses the common problems of multi-objective conflict, frequent disturbances, and high timeliness requirements for decision-making in multi-variety, small-batch mixed-flow production in flexible operation workshops, providing an engineering implementation path. The method provided by this invention consists of the following main steps: data acquisition and preprocessing, construction of a multi-objective scheduling model, solution using a non-dominated sorting genetic algorithm based on dual-archive co-evolution, objective decision-making based on entropy weight and approximation of the ideal solution sorting method, disturbance severity quantification based on disjunctive graph right-shift propagation simulation, hierarchical rescheduling based on disturbance severity, and scheme output and sample feedback, denoted as steps S1 to S7 respectively. (Refer to...) Figure 5 The seven main steps can be executed sequentially from step S1 to step S4 at the beginning of the scheduling cycle to complete the initial scheme issuance, and steps S5 to S7 form a continuously operating disturbance response loop during the scheduling cycle.

[0088] In step S1, the input information required for production scheduling is collected from the workshop production data source and standardized preprocessing is completed. Step S1 consists of two sub-steps: step S11 and step S12.

[0089] In step S11, the Manufacturing Execution System (MES) reads the workpiece number, the process number, the set of processable equipment for each process, the standard processing time for each process on each processable equipment, the release time of each workpiece, the delivery date of each workpiece, and the lateness penalty weight. The IoT interface reads the current availability status, remaining processing time, and expected release time of each device from the Programmable Logic Controller (PLC). The Material Management System (MMS) reads the current inventory and available time window of each material. The MES is the authoritative data source for shop floor orders and process flows. The PLC reflects the real-time operating status of physical equipment. The MMS constrains material completeness. These three types of data together constitute the complete input of the flexible workshop scheduling model in the time, resource, and material dimensions.

[0090] In step S12, the multi-source data collected in step S11 is mapped to the memory data structure of the production scheduling solver according to preset fields. Considering the various realities of inconsistent formats, different units, and missing data in the field, step S12 performs differentiated cleaning on different types of fields. Specifically, for missing numerical fields such as standard processing time, release time, remaining processing time, and inventory, historical averages or averages of similar processes are used to fill in the missing fields to ensure the model is solvable. For missing fields such as workpiece number, process sequence, and set of processable equipment, no filling is performed; instead, an integrity check is conducted. Orders that fail the integrity check are marked as awaiting manual confirmation to avoid misscheduling caused by mistakenly filling in non-numerical fields based on historical averages. Then, the units of each field are uniformly converted, with time fields uniformly converted to min and quantity fields converted according to the measurement units of the corresponding materials. Finally, the processed data is organized into workpiece sets, equipment sets, process sets, processing time matrices, and constraint parameter lists, which serve as the production scheduling input dataset for subsequent steps.

[0091] It should be noted that "min" as a time unit is a commonly used granularity in workshop scheduling. Its alignment with machine cycle time and manual shift scheduling is simple and direct, facilitating subsequent analysis, modeling, and solving. Preferably, for fields with hours, seconds, or days as the original unit, integerization can be performed in step S12 to reduce the accumulation of floating-point errors. For material measurement unit conversion, it is preferable to perform the conversion according to the smallest measurement unit specified in the material process specification. For example, metal raw materials measured in kilograms, standard parts measured in pieces, and profiles measured in meters should be converted according to their respective smallest measurement units.

[0092] In step S2, a multi-objective flexible workshop scheduling model is constructed. The model uses a binary decision variable—whether the j-th process of workpiece i is assigned to the k-th machine—as a basis for decision. and the actual start time of this process. As a decision variable; among which A value of 1 indicates that the j-th process of workpiece i is assigned to the k-th machine for processing, while a value of 0 indicates that it is not assigned. Let be the actual start time (min) of the j-th process of workpiece i; i is the workpiece number; j is the process number within the same workpiece; k is the equipment number.

[0093] The optimization objectives of the model are to minimize the maximum completion time, minimize weighted delay, and maximize equipment utilization. Maximum completion time refers to the time when the last workpiece in the current batch is completed, reflecting the overall production cycle length. Weighted delay is the sum of the products of the actual completion time of each workpiece being later than its delivery date and the corresponding delay penalty weight, reflecting the degree of deviation from customer delivery commitments. Equipment utilization is the average ratio of the effective processing time of each piece of equipment to the available time of that equipment, reflecting the efficiency of equipment resource utilization. The smaller the first two objectives, the better; the larger the third objective, the better.

[0094] The hard constraints of the model include six categories. Category 1 is the process sequence constraint, requiring a strict sequential order between adjacent processes of the same workpiece, meaning the start time of the later process cannot be earlier than the completion time of the earlier process. Category 2 is the unique allocation constraint of process equipment, requiring each process to be assigned to one piece of equipment from its set of available equipment in each feasible solution of the model. Category 3 is the mutual exclusion constraint for processing on the same equipment, requiring that one process be processed on the same equipment at the same time, and that the processing time intervals of different processes on the same equipment do not overlap. Category 4 is the workpiece release time constraint, requiring that the start time of the first process of a workpiece cannot be earlier than the release time of the workpiece. Category 5 is the equipment manufacturability constraint, requiring each process to be assigned to equipment within the set of available equipment for that process. Category 6 is the material availability constraint, requiring that the start time of a process be within the available time window of the materials required for that process, and not exceed the current inventory of those materials.

[0095] The Flexible Job Shop Scheduling Problem (FJSP) differs from the traditional job shop scheduling problem in that each process has multiple candidate machines to choose from. Therefore, solving it requires determining both the machine allocation and the start time of each process, resulting in a higher decision-making dimension and greater difficulty. This invention incorporates three objectives into the optimization model and employs a multi-objective evolutionary method for solving it. This avoids the subjectivity inherent in manually selecting weights and reducing multiple objectives to a single objective. The Pareto front objectively reflects the trade-offs between the three objectives, allowing decision-makers to select the appropriate execution plan based on the actual preferences of the current production task.

[0096] In step S3, the model constructed in step S2 is solved using a non-dominated sorting genetic algorithm based on dual-archive co-evolution. The standard form of the non-dominated sorting genetic algorithm is abbreviated as NSGA-II. This invention introduces two improvements to NSGA-II: robust buffer embedding, a key process in the decoding stage, and dual-archive co-evolution in the selection stage. Step S3 consists of seven sub-steps, from S31 to S37.

[0097] In step S31, a double-layer chromosome is used as the encoding for each individual in the algorithm. The first segment of the double-layer chromosome is called the process sequencing segment, and its length is the sum of the number of processes for each workpiece. Each gene position in the process sequencing segment is filled with the corresponding workpiece number, and the number of times the same workpiece number appears in this segment is equal to the number of processes for that workpiece. The number of workpiece i that appears p times from left to right represents the relative priority of the p-th process of workpiece i during decoding, where p is the cumulative number of times the workpiece number appears in the process sequencing segment. The second segment of the double-layer chromosome is called the equipment allocation segment, and its length is the same as that of the process sequencing segment. Each gene position in the equipment allocation segment is filled with the sequence number of the equipment selected by the corresponding process in its set of processable equipment. The process sequencing segment determines the scheduling order, and the equipment allocation segment determines the resource allocation. Together, they constitute a complete description of the flexible job shop scheduling solution.

[0098] In step S32, a hybrid strategy is used to generate the initial population. Specifically, a portion of the individuals in the initial population are generated using a global selection initialization method. This method prioritizes equipment with lower cumulative allocated loads in the set of equipment available for processing for each process. The remaining individuals in the initial population are generated using a completely random method, where one piece of equipment is randomly selected with equal probability from the set of equipment available for processing for each process as its equipment number. The former focuses on the initial quality of the population, ensuring that the maximum completion time distribution of the initial population is biased towards a more favorable range; the latter focuses on the diversity of the population, preventing the evolutionary process from prematurely falling into local optima. Preferably, the proportion of the individuals in the hybrid strategy can be 50% to 80% of the initial population size.

[0099] In step S33, priority process crossover and positional uniform crossover are performed on the parent generation individuals according to a preset crossover rate. Priority process crossover is abbreviated as POX. The specific process of POX crossover is as follows: a subset of workpieces is randomly selected in parent generation 1, and all processes of the workpieces contained in the subset are copied to the same positions in the offspring according to their relative order in parent generation 1. The remaining positions are filled in sequentially according to the relative order of the remaining workpieces in parent generation 2, excluding the copied workpieces. Positional uniform crossover is used for the equipment allocation segment, and each gene position is copied from parent generation 1 or parent generation 2 with a probability of 0.5. Mutation is performed at a preset mutation rate. The mutation of the process order segment is carried out by randomly swapping the positions of two different workpiece numbers, introducing perturbation while maintaining the legality of the process order. The mutation of the equipment allocation segment is not performed equally on all gene positions, but only on the gene positions corresponding to the critical path processes identified by decoding in step S34, using the following rule: when there are at least two machines in the set of processable equipment for that process, the gene position is changed to the sequence number of the candidate machine whose standard processing time is less than that of the currently selected machine; otherwise, the mutation of the equipment allocation segment is not performed. The basis of this targeted mutation strategy is that the processing time of the critical path process directly determines the maximum completion time, and mutating its equipment allocation has a higher objective function improvement expectation than mutating non-critical path processes, thereby accelerating convergence.

[0100] In step S34, decoding is performed on each individual to obtain the preliminary production scheduling plan corresponding to that individual, and the set of processes contained in the critical path of the plan is identified. Specifically, based on the individual's double-layer chromosome, the processes are first arranged sequentially to the earliest available time window of the corresponding equipment according to a predetermined decoding rule to obtain the preliminary production scheduling plan; then, a directed acyclic disjunctive graph corresponding to the preliminary production scheduling plan is constructed. The disjunctive graph is a classic graph model of the job shop scheduling problem, where nodes correspond to processes, and arcs are divided into two categories: one is the process sequence arc, which connects two adjacent processes of the same workpiece, and the other is the resource occupancy arc, which connects two adjacent processing processes on the same equipment; the source node of the disjunctive graph corresponds to a virtual start event, and the sink node corresponds to a virtual end event.

[0101] On the disjunctive graph, traverse forward from the source node to the sink node in topological order. For each node, take the maximum value among the earliest completion times of its direct predecessor nodes as the earliest start time of that node, and add it to the standard processing time of the process on the allocated equipment to obtain the earliest completion time of that node. After the traversal, the earliest start time and earliest completion time of each process are obtained, and the largest earliest completion time is the global maximum completion time of the preliminary production scheduling plan. Then, using this global maximum completion time as the endpoint, traverse backward from the sink node to the source node in topological reverse order. For each node, take the minimum value among the latest start times of its direct successor nodes as the latest completion time of that node, and subtract it from the standard processing time of the process on the allocated equipment to obtain the latest start time of that node. Finally, identify the processes whose difference between the latest start time and the earliest start time is zero as critical path processes, and the corresponding process set is denoted as CP. The difference between the latest start time and the earliest start time is called the total float time of the process. A total float time of zero means that any delay in the start or end of the process will be directly passed on as a delay in the maximum completion time.

[0102] The critical path identification method combining forward and backward traversal is known as the Critical Path Method (CPM) in the field of engineering network planning. It has a time complexity of O(n+a) on deterministic networks, where n is the number of nodes and a is the number of arcs. It is also effective on disjunctive graphs of flexible job shops. On medium-sized disjunctive graphs, a single critical path identification can be completed in milliseconds, making it suitable for use in each generation of the decoding loop.

[0103] In step S35, for each key process in the process set CP Perform robust buffer embedding. (The key process...) The j-th process of workpiece i is assigned to equipment k in the preliminary production scheduling plan. The key process is retrieved from the disturbance statistics sample knowledge base. Standard deviation of processing time The key process is calculated according to the following formula. The length of the robust buffer that should be embedded.

[0104] ;(Formula 1);

[0105] In formula 1: Let be the robust buffer length (min) that should be embedded before the j-th process of workpiece i when it is processed on equipment k. It is a robust multiple, dimensionless, and a positive real number; Let be the standard deviation (in minutes) of the processing time of the j-th operation on equipment k for workpiece i; i is the workpiece number; j is the operation number within the same workpiece; k is the equipment number.

[0106] Formula 1 is based on Chebyshev's inequality in probability theory. This inequality, given by the Russian mathematician Chebyshev in the 19th century, states that for any random variable X and any positive number... X deviates from its mean by more than The probability of more than one standard deviation does not exceed Specifically, regarding the random fluctuations in processing time in this invention, When the value is 1, the probability that the processing time deviates from the mean by more than one standard deviation is no more than 1; when the value is 2, the probability is no more than 1 / 4; and when the value is 3, the probability is no more than 1 / 9. Therefore, the robustness factor... Standard deviation of processing time Before multiplying and embedding the buffer length into critical processes, first-order random fluctuations can be absorbed with a controllable probability, improving the robustness of the solution to time disturbances without significantly increasing the maximum completion time. Preferably, the robustness factor... A value of 1 can be chosen to achieve a better balance between probability constraints and scheme compactness; when the workshop time statistics show that the processing time tail is thick, It can be appropriately increased; when the workshop has high requirements for the compactness of the production cycle, It can be reduced appropriately.

[0107] Then the obtained robust buffer length Before the start time of the corresponding process is inserted as an additional delay, the earliest start time and earliest completion time of the direct successor nodes of the same workpiece, the direct successor nodes of the same equipment, and the recursive successor nodes of the aforementioned successor nodes are updated along the disjunction graph. The recursive successor nodes refer to all subsequent nodes obtained by continuously tracing back along the process sequence arc and resource occupancy arc of the disjunction graph from the aforementioned two types of direct successor nodes. The update logic ensures that the propagation of the buffer embedding along the time axis of the scheme is fully reflected, avoiding time conflicts caused by local insertion of buffers without corresponding right shift of subsequent nodes.

[0108] The updated solution is then used as input for step S34, and steps S34 and S35 are repeated until the critical path obtained from two consecutive iterations contains the same set of processes or the number of iterations reaches a preset upper limit. The feasible scheduling solution corresponding to this individual is then output. The necessity of repeated iterations lies in the fact that buffer embedding changes the start time of processes, potentially altering the critical path itself. Therefore, it is necessary to re-evaluate and embed buffers on the new critical path until the critical path stabilizes. Preferably, the iteration upper limit can be 2 to 3 times; in most instances, the critical path stabilizes after 2 iterations.

[0109] In step S36, the set obtained by merging the main evolutionary population with the current generation's offspring is divided into several non-dominated levels according to the non-dominated ranking. The dominance relationship between two solutions in the non-dominated ranking means that one solution is no worse than the other in all objectives and is strictly superior to the other in at least one objective. The set of solutions with non-dominated level 1 is the set of elite solutions in the merged set that do not dominate each other; the solutions with non-dominated level 1 are considered the Pareto non-dominated front of the current generation. The set of solutions with non-dominated level 2 is the set of non-dominated solutions remaining after removing solutions with non-dominated level 1 from the merged set, and so on. Then, the convergence file is filled in sequentially according to the non-dominated level from low to high, and within the same non-dominated level, according to the crowding distance from large to small. The crowding distance, given by Deb et al. in NSGA-II, is calculated for each solution in each objective dimension as the sum of the objective differences between its two neighboring solutions, and then normalized to reflect the sparsity of the local region where the solution resides. A larger crowding distance indicates that the solution's support for the frontier is more independent and therefore more worthy of preservation. The filling operation continues until the convergence archive is closed. Fill the preset capacity. Preferably, converge the file. The preset capacity is consistent with the population size, so that each generation can retain elite solutions of the same size as the population in the archive.

[0110] In step S37, the distributed archive is constructed and maintained. Specifically, the objective value of each non-dominated solution in the merged set is normalized using the range of each objective to obtain the normalized objective vector of that non-dominated solution. This range normalization maps the value of each objective dimension to the interval between 0 and 1, avoiding distortion in Euclidean distance calculation caused by inconsistent dimensions of different objectives. When the archives are scattered... When the set is empty, the non-dominated solution with the largest Euclidean distance between the normalized objective vector and the current ideal point in the merged set is first used as the dispersed file. The initial solution; the current ideal point refers to the virtual vector composed of the currently achievable optimal values ​​of each objective. In the distributed archives When not empty, calculate each non-dominated solution in the merged set according to the following formula. Relative to the scattered archives Riesz2-energy contribution value.

[0111] ;(Formula 2);

[0112] In formula 2: Non-dominated solution Relative to the scattered archives The Riesz2- energy contribution value is dimensionless. Non-dominated solution With the aforementioned scattered archives China-Africa Domination Solution The Euclidean distance over the normalized target vector is dimensionless. It is a preset minimum positive real number used to prevent division by zero when the distance is zero; Indicates taking and The larger of the two; For the distributed archive; i is the number of the non-dominated solution to be evaluated; l is the distributed archive. The numbers of the non-dominated solutions already included in the database; This indicates the location of the scattered archives. Unless the solution is dominated Summing the remaining non-dominated solutions; the superscript 2 indicates the squaring operation.

[0113] Formula 2 originates from the Rieszs-energy in potential energy theory, which makes the form on a compact manifold resemble... The finite set of points whose summation reaches its minimum value tends to a uniform distribution on the manifold as s approaches infinity. When s is 2, the summation form is concise and has a natural strong repulsion effect on boundary points, maintaining support solutions at the center and edges simultaneously on multi-objective fronts. This invention uses the energy contribution value... As a measure of non-dominated solutions An indicator of the degree of dispersion contribution; the smaller the energy contribution value, the more non-dominated the solution. The further away from existing solutions in the archive, the better it can serve as a support point for sparse regions. (Introduction) This avoids numerical instability where the denominator is zero in cases of repeated or near-repeated solutions. Preferably, Desirable Magnitude.

[0114] Then press Fill the remaining non-dominated solutions in the merged set into the distributed file in ascending order. until Fill the preset capacity. Preferably, distribute the files. Preset capacity and convergence file The preset capacity remains consistent.

[0115] To strengthen the collection of archives With scattered archives Information exchange between them avoids both becoming trapped in local preferences, and elite migration is performed algebraically at preset exchange intervals. Specifically, from the convergence archives... Select the non-dominated solution for the dispersed migration of the current congestion distance. and from scattered archives Select the non-dominated solution with the smallest current non-dominated level value and migrate it into the convergent archive. Relocation and collection of archives Non-dominated solutions from scattered archives After removing from the middle, press From smallest to largest, select an undoed non-dominated solution from the merged set and add it to the scattered file. The logical design basis for the elite migration is based on the convergence of archives. To scattered archives Relocate solutions that are crowded and have large distances, and supplement the scattered archives with solutions that have better convergence and are more sparsely distributed; from the scattered archives To collect files The solution with the smallest non-dominant rank value is removed, and solutions with better dispersion and competitive elite status are added to the convergence file; the combination of these two approaches synergistically enhances the diversity and eliteness of the solution set. Preferably, the preset exchange interval can be about 10 generations to balance the interaction frequency and algorithm stability.

[0116] When the iteration termination condition is met, the convergence file will be closed. With scattered archives The union of the solutions is used to remove inferior and duplicate solutions dominated by other solutions. The resulting solution set is the Pareto non-dominated solution set output to step S4. The iteration termination condition can be set to the current evolution generation reaching a preset maximum generation, or the convergence file. The Hypervolume metric (i.e., the hypervolume index) is defined as the relative growth rate of a function that does not exceed a preset threshold over several consecutive generations. The Hypervolume metric measures the target space volume of the Pareto front envelope and is a commonly used indicator for evaluating the convergence of multi-objective evolutionary algorithms.

[0117] In step S4, an objective decision is made on the Pareto non-dominated solution set output in step S3 to obtain the initial execution plan to be issued to the workshop. Step S4 consists of three sub-steps, S41 to S43.

[0118] In step S41, the values ​​of each candidate solution in the Pareto non-dominated solution set for the three objectives of maximum completion time, weighted delay, and equipment utilization are arranged by row as candidate solutions and by column as objectives, resulting in an N-row, 3-column decision matrix; where N is the total number of candidate solutions in the Pareto non-dominated solution set. When N equals 1, it indicates that the Pareto front contains only one candidate solution, and no further comparison is needed. This unique candidate solution is directly issued as the initial execution solution to the workshop execution interface, and step S4 ends; when N is greater than 1, proceed to step S42.

[0119] In step S42, range normalization is first performed on the decision matrix. Specifically, the maximum completion time column and the weighted delay column of the decision matrix are normalized according to cost-type indicators, and the equipment utilization rate column is normalized according to benefit-type indicators, resulting in a normalized decision matrix. The range normalization of the cost-type indicators is to subtract the difference of each element from the maximum value of the column and divide by the range of the column; the range normalization of the benefit-type indicators is to divide the difference between each element and the minimum value of the column and divide by the range of the column. Both methods ensure that the normalized values ​​fall within the range of 0 to 1, and both methods map the better target value to a larger normalized value. When the range of a target column is zero, the normalized value of that column is uniformly set to a preset constant of 0.5. The constant 0.5 is based on the fact that all candidate schemes in that column have the same value for that target. At this time, the target has no discriminatory power for decision-making. Setting the normalized value to the median of 0.5 can maintain the symmetry of the subsequent entropy weight calculation without bias in any direction.

[0120] Then, the Shannon entropy of the q-th column of the normalized decision matrix is ​​calculated according to the following formula.

[0121] ;(Formula 3);

[0122] In formula 3: Let be the Shannon entropy of the q-th column of the normalized decision matrix, which is dimensionless; q is the target number, where q takes the values ​​1, 2, and 3, corresponding to the three targets: maximum completion time, weighted delay, and equipment utilization rate, respectively; N is the total number of candidate solutions. The normalized value in the nth row and qth column of the normalized decision matrix is ​​the proportion of the sum of the normalized values ​​in the qth column, which is dimensionless; n is the candidate scheme number. The natural logarithm operator; This represents the summation over n=1 to n=N.

[0123] Formula 3 provides a measure of the uncertainty of a discrete random variable, called information entropy or Shannon entropy, which is defined as the negative log-expected value of the probability distribution. In the field of objective weighting, the values ​​of candidate solutions on a certain objective are considered as probabilistic events, and the resulting entropy value reflects the degree of dispersion of the objective across all candidate solutions. The normalized values ​​are then normalized column-wise to obtain a probability distribution. Substituting this into the Shannon entropy formula, we get... Pre-factor It is a normalized form of Shannon entropy, making The value of falls within the range of 0 to 1, which facilitates subsequent comparisons. When When it is zero, press This mathematical fact allows us to treat the corresponding terms as zero, thus preventing logarithmic divergence.

[0124] Then, the objective entropy weight corresponding to the q-th target is calculated according to the following formula.

[0125] ;(Formula 4);

[0126] In formula 4: Let be the objective entropy weight corresponding to the q-th target, which is dimensionless; For the summation index, Choose 1, 2, or 3; For the normalized decision matrix, the first... The Shannon entropy of the column is dimensionless; Indicates to =1 to =3, sum them up.

[0127] Formula 4 is based on, The larger the value, the more uniform the distribution of normalized values ​​among the candidate solutions in column q, and the weaker the ability of this objective to distinguish between different options in the decision-making process; therefore, it should be assigned a smaller weight. Conversely, a smaller value indicates a more uniform distribution of normalized values ​​among the candidate solutions in column q. The smaller the value, the more concentrated the distribution of the q-th column, and the stronger the ability of the objective to distinguish decisions, so it should be given a larger weight. That is, the distinguishing power of the target is quantified, and then the objective entropy weight is obtained by normalizing it according to the sum of the distinguishing powers of all targets. . satisfy The entropy weighting method is a representative of objective weighting methods based on Shannon entropy.

[0128] In step S43, the normalized decision matrix is ​​weighted according to entropy weights, and the relative closeness of each candidate solution to the positive and negative ideal solutions is recalculated. Specifically, each element in the q-th column of the normalized decision matrix is ​​multiplied by its corresponding objective entropy weight. Obtain the weighted normalized decision matrix; select the maximum value from each column of the weighted normalized decision matrix to form the positive ideal solution and select the minimum value to form the negative ideal solution; calculate the relative closeness of the nth candidate solution to the positive ideal solution according to the following formula.

[0129] ;(Formula 5);

[0130] In formula 5: Let be the relative closeness of the nth candidate solution to the positive ideal solution, which is dimensionless; The Euclidean distance between the nth row vector of the weighted normalized decision matrix and the positive ideal solution is dimensionless. The Euclidean distance between the nth row vector of the weighted normalized decision matrix and the negative ideal solution is dimensionless; the superscript plus sign identifies the positive ideal solution; the superscript minus sign identifies the negative ideal solution.

[0131] Formula 5 is the ranking method for approximating the ideal solution, abbreviated as TOPSIS. The core idea of ​​TOPSIS is that in multi-criteria decision-making, the optimal solution should simultaneously satisfy the conditions of being closest to the positive ideal solution and furthest from the negative ideal solution. The positive ideal solution consists of the best values ​​of each objective in the candidate set, while the negative ideal solution consists of the worst values ​​of each objective in the candidate set. This is achieved through relative proximity... This scalar measures both of the above requirements. The value ranges from 0 to 1. The closer it is to 1, the closer the candidate solution is to the positive ideal solution and the further it is from the negative ideal solution, and the better it is overall.

[0132] Finally, the relative closeness The candidate solution with the largest value is sent to the workshop execution interface as the initial execution solution, and step S4 ends. This invention adopts an objective decision-making method combining entropy weight and TOPSIS. Compared with the practice of directly selecting solutions from the Pareto front by humans based on experience, it avoids the introduction of subjective weights and makes the decision-making process repeatable and traceable.

[0133] In step S5, the impact of disturbance events occurring in real time during workshop execution is assessed and the severity of the disturbances is quantified, providing a basis for subsequent rescheduling strategy matching. Step S5 consists of three sub-steps, from step S51 to step S53.

[0134] In step S51, disturbance events e are continuously received through the disturbance alarm interface. The disturbance alarm interface can be triggered collaboratively by the workshop manufacturing execution system, programmable logic controller, and operator tablet terminal, covering three types of disturbance sources: equipment layer, order layer, and process layer. Based on the event type, the parsed graph corresponding to the currently executing production scheduling plan is processed in one of the following ways.

[0135] In category 1, when the disturbance event e is a device failure, the equivalent node delay of the disjunction graph node corresponding to the currently processing or next processing step on the failed device is set to the estimated fault recovery time (min). The estimated fault recovery time is provided by the programmable logic controller or maintenance system. During the device failure, the device cannot accept new processes; therefore, the delay of the failed node is set to the recovery time to reflect its rightward shift effect on the disjunction graph.

[0136] In the second category, when the disturbance event e is an urgent order insertion, the urgent order is added to the disjunction graph as a set of newly added workpiece nodes. The release time of the first operation of the newly added workpiece is set as the urgent order insertion time. The operation sequence constraints, equipment manufacturability constraints, and material availability constraints of the newly added workpiece are also added to the disjunction graph. The essence of an urgent order is a newly added workpiece and its operations. Compared to equating an urgent order with existing node delays, modeling an urgent order as a set of newly added nodes more accurately reflects its true structure in the scheduling problem.

[0137] In the third category, when the disturbance event e is a processing time deviation, the equivalent node of the disjunctive graph node corresponding to the process whose deviation exceeds the preset tolerance is delayed to the difference (min) between the actual processing time and the standard processing time of that process. The preset tolerance is used to filter random noise and avoid triggering an evaluation for every tiny fluctuation.

[0138] In step S52, a right-shift propagation simulation is performed using the processed disjunctive graph as input. Specifically, the topological order is traversed forward from the source node to the sink node, and the following operations are performed on each node: the new start time of the node is taken as the larger of the maximum value of the new completion times of all its direct predecessor nodes and the node's release time constraint; the new completion time of the node is taken as the sum of the new start time, standard processing time, and equivalent node delay. After the traversal is completed, the updated start time and completion time of each process are obtained, and the maximum value of the updated completion times is recorded as... The time complexity of the topological order traversal is O(number of nodes + number of arcs), which can be completed in milliseconds on medium-sized disjunctive graphs, making it suitable for real-time evaluation. The right-shift propagation simulation does not re-execute the optimization solution; it only shifts the delay node by node according to the existing scheduling order and device allocation, obtaining a measure of the instability of the scheme under the assumption of no reordering.

[0139] In step S53, the severity of the disturbance event e is calculated according to the following formula.

[0140] ;(Formula 6);

[0141] In formula 6: Let be the severity of the disturbance event e, which is dimensionless; The global maximum completion time (min) after the right-shift propagation simulation obtained in step S52; The global maximum completion time (min) of the current production schedule before the occurrence of disturbance event e; This indicates taking the larger of a and b; the superscript RS indicates the state after the right-shift propagation simulation; the superscript plan indicates the planned state before the disturbance.

[0142] Formula 6 defines the severity of the disturbance as the relative increment of the maximum completion time after the right-shift propagation simulation relative to the maximum completion time before the disturbance, and uses... The function takes a non-negative truncation. Its design principle is that the most direct negative impact of disturbances on the scheduling scheme is the delay of the maximum completion time; the dimensionless quantity of relative increment facilitates comparison across workshops and workpiece batches and the unified setting of thresholds; The function uses non-negative truncation to handle the rare case where the maximum completion time does not increase after right-shift propagation, thus avoiding logical errors caused by negative values ​​in threshold comparisons. The perturbation severity... Step S5 ends, serving as the basis for response strategy matching in step S6.

[0143] Through steps S51 to S53, disturbance events of different types and intensities are uniformly mapped to a single dimensionless index. This simplifies the matching rules for subsequent rescheduling strategies and decouples them from specific disturbance types. Compared to existing technologies that define response rules separately for each type of disturbance, the unified quantitative index of this invention reduces rule maintenance costs and facilitates the accumulation of samples by the knowledge base self-maintenance mechanism through a unified interface.

[0144] In step S6, based on the severity of the disturbance... Select and execute one of the three rescheduling strategies to obtain the updated solution. Step S6 consists of four sub-steps, from step S61 to step S64.

[0145] In step S61, the severity of the disturbance is... With the slight threshold maintained in the perturbation statistics sample knowledge base and severe threshold The comparison is performed, and the response strategy is matched according to the following rules. When Perform rightward repair when; Perform local re-optimization of the critical path when needed; A global scroll reordering is performed at that time. The threshold value is light and dimensionless. The severity threshold is dimensionless. When the cumulative number of samples in the disturbance statistics sample knowledge base is less than a preset sample threshold, and A preset initial threshold is taken from the disturbance statistical sample knowledge base. This preset initial threshold is used to address the engineering reality that the sample size is insufficient to support automatic maintenance during the initial stage of system operation. Preferably, the initial threshold... A value of approximately 0.10 can be used as the initial threshold. A value of approximately 0.30 can be taken, corresponding to the understanding from engineering experience that a delay of less than 10% is mildly absorbable and a delay of more than 30% requires rearrangement.

[0146] In step S62, the right-shift repair maintains the equipment allocation for each process in the current production schedule and the processing sequence of processes on the same equipment unchanged. The new start times for each process obtained in step S52 are used to overwrite the original start times, and this is output as the updated plan. The computational cost of the right-shift repair is approximately zero. Its engineering value lies in completing the plan repair for small disturbances at minimal cost, avoiding on-site execution fluctuations caused by unnecessary rearrangement.

[0147] In step S63, the local re-optimization of the critical path uses the right-shift propagation simulation scheme obtained in step S52 as the initial value, identifies its critical path, and unlocks the start time and equipment allocation of the processes on the critical path, the direct successor processes of the same workpiece on the critical path, the direct successor processes of the same equipment on the critical path, and the adjacent processes related to equipment reassignment as decision variables. The start time and equipment allocation of the remaining processes are fixed as determined parameters. The adjacent processes related to equipment reassignment refer to the fixed processes immediately adjacent to the process on the newly assigned equipment when the equipment allocation of the process on the critical path changes during re-optimization, which also need to be included in the decision variable set to avoid time conflicts on the new equipment during re-optimization. Then, the dual-file co-evolutionary non-dominated sorting genetic algorithm of step S3 is called to perform small-scale re-optimization on the reduced subproblem. During the re-optimization process, conflict verification is performed on the time windows of the fixed processes on the reassigned equipment, and the insertionable time window constraint ensures that the processing of the same equipment is mutually exclusive. Then, the Pareto non-dominated solution set obtained from the re-optimization is decided according to step S4 to obtain the updated scheme. Preferably, the maximum number of generations for the small-scale re-optimization can be taken as 1 / 3 to 1 / 5 of the maximum number of generations for the global optimization. The local re-optimization of the critical path achieves a balance between the solution scale and the extent of solution improvement, significantly reducing the computation time compared to global rearrangement and providing greater room for solution improvement compared to right-shift repair.

[0148] In step S64, the global rolling rescheduling fixes the processes that have started but not yet completed at the current time onto their original assigned equipment for continued processing. The start time of each process is fixed to the current time, and the remaining processing time is fixed to the remaining processing time of that process. The processes that have not yet started are treated as a new set of processes to be scheduled, and the updated equipment status, material status, and remaining orders are used as new scheduling inputs. Steps S2, S3, and S4 are then re-executed to obtain an updated plan. Global rolling rescheduling is suitable for situations where the disturbance is so severe that the overall structure of the original scheduling plan is no longer applicable. Its cost is a longer solution time, but under significant disturbances, it can fundamentally re-optimize the resource allocation for the entire remaining production period.

[0149] Through steps S61 to S64, the present invention forms a hierarchical response mechanism that matches the severity of disturbances with the response strategy. Small disturbances are responded to at a low cost, while large disturbances are responded to at a high cost, thus keeping the average rescheduling cost at a reasonable level overall.

[0150] In step S7, an updated plan and decision report are output, and the disturbance sample is fed back to the disturbance statistical sample knowledge base to support the self-maintenance of the knowledge base. Step S7 consists of two sub-steps, step S71 and step S72.

[0151] In step S71, the updated plan is output as a two-dimensional Gantt chart in the workpiece-equipment dimension, along with a decision report containing three target values, critical path identifiers, and buffer embedding location identifiers. The Gantt chart uses either the workpiece or the equipment as the vertical axis and time as the horizontal axis. Each rectangle represents the processing time interval of a process on a single piece of equipment, allowing on-site schedulers to intuitively grasp the overall picture. The target values ​​in the decision report reflect the specific level of the plan in terms of the three targets; the critical path identifier assists schedulers in identifying bottlenecks and monitoring the execution progress of critical processes; and the buffer embedding location identifiers clearly indicate which critical processes have reserved robust buffers, facilitating on-site understanding of why these areas are incompressible.

[0152] In step S72, the characteristic parameters of this disturbance event and the type of response strategy adopted are fed back as a sample to the disturbance statistical sample knowledge base; the characteristic parameters include the disturbance event type, the affected process number, the affected equipment number, the equivalent node delay value (min), and the disturbance severity. When the cumulative number of samples in the perturbation statistical sample knowledge base reaches the preset sample threshold, the perturbation statistical sample knowledge base updates the mild threshold based on the preset quantiles of the relative increment distribution of the maximum completion time actually generated by the right-shift repair strategy in the cumulative samples. and the severity threshold Based on the deviation between the actual processing time and the standard processing time of each process on different equipment in the cumulative sample, the corresponding standard deviation of the processing time is updated on a rolling basis according to the workpiece-process-equipment combination. The rolling update enables , and As the workshop's operational history continuously approximates the statistical characteristics of the true distribution, the threshold matching in step S61 and the robust buffer embedding in step S35 are continuously optimized over time. Preferably, the preset sample threshold can be approximately 30 to 50 samples to meet the minimum sample requirement for quantile estimation; the preset quantiles can be approximately 25% and 75% quantiles, respectively, to correspond to the mild and severe thresholds.

[0153] Through steps S71 and S72, this invention forms a closed loop with four interconnected stages: scheduling, execution, sample processing, and knowledge acquisition. This allows the effectiveness of the method to improve as workshop operating experience accumulates. In the initial stage of operation, the method relies on the initial threshold and the prior value of the standard deviation of working hours given by the optimized selection. As the samples accumulate above the preset threshold, the method switches to a self-maintaining knowledge base, whose threshold and standard deviation of working hours gradually reflect the true statistical characteristics of the workshop. This makes the subsequent rescheduling strategy matching more accurate and the buffer embedding more appropriate.

[0154] The system implementation of the present invention will be described in detail below. The system provided by the present invention is called a multi-objective dynamic scheduling decision system, which is used to implement the aforementioned method. The system includes seven primary modules: a data acquisition and preprocessing module, a multi-objective optimization engine module, an automatic scheme decision module, a disturbance monitoring and impact assessment module, an adaptive rescheduling decision module, a disturbance statistical sample knowledge base module, and a scheme output and visualization module. The seven primary modules are interconnected with the system's internal data bus via an industrial Ethernet, and correspond to the knowledge base self-maintenance mechanism in steps S1, S3, S4, S5, S6, and S72 of the method, as well as the hardware implementation of step S7.

[0155] The data acquisition and preprocessing module communicates with the workshop manufacturing execution system, enterprise resource management system, material management system, equipment IoT interface, and programmable logic controller (PLC) via industrial Ethernet, forming the system's data entry point. Internally, this module further includes a manufacturing execution system data acquisition unit, a PLC data acquisition unit, a material data acquisition unit, and a data standardization and mapping unit. The outputs of the first three acquisition units are connected to the input of the data standardization and mapping unit. This unit performs four functions: field mapping, missing data filling, unit unification, and integrity verification. Its output, serving as the output of the data acquisition and preprocessing module, is connected to the input of the multi-objective optimization engine module via the system's internal data bus. The core function of this module is to standardize and integrate multi-source heterogeneous production data, ensuring that the input received by the subsequent optimization engine is consistent in structure and dimensions.

[0156] The multi-objective optimization engine module includes a model building submodule, a key process robust decoding submodule, and a dual-file co-evolution submodule connected in sequence. The model building submodule constructs in-memory representations of decision variables, objective functions, and hard constraints based on the production scheduling input dataset output by the data acquisition and preprocessing module, undertaking the hardware implementation of method step S2. The key process robust decoding submodule further includes a critical path identification unit and a robust buffer embedding unit. The critical path identification unit is used to identify the set of processes contained in the critical path by topological forward traversal and topological reverse traversal, corresponding to the hardware implementation of method step S34. The robust buffer embedding unit is used to embed a robust buffer before the critical path processes identified by the critical path identification unit, corresponding to the hardware implementation of method step S35. The second input of the robust buffer embedding unit is connected to the output of the perturbation statistical sample knowledge base module for real-time reading of the processing time standard deviation. The dual-archive co-evolution submodule includes a convergent archive storage unit, a dispersed archive storage unit, a Riesz2-energy calculation unit, and an elite migration control unit. The output of the Riesz2-energy calculation unit is connected to the input of the dispersed archive storage unit, and is used to sort and fill the dispersed archives according to the energy contribution values ​​calculated by Formula 2. The elite migration control unit is bidirectionally connected to both the convergent archive storage unit and the dispersed archive storage unit, and is used to perform elite migration between archives at preset exchange intervals, corresponding to the hardware implementation of steps S36 and S37. The output of the multi-objective optimization engine module is connected to the input of the automatic scheme decision module through the system's internal data bus, and is responsible for sending the Pareto non-dominated solution set to the decision module.

[0157] The input of the automatic decision-making module is connected to the output of the multi-objective optimization engine module via the system's internal data bus. The output is connected to the first input of the solution output and visualization module and the workshop execution interface. The automatic decision-making module further includes a decision matrix construction unit, an entropy weight calculation unit, a proximity calculation unit, and an optimal solution selection unit. These four units are connected in series and respectively undertake the hardware implementation of steps S41, S42, and S43. The core function of this module is to convert the Pareto non-dominated solution set into a single executable initial execution plan and distribute it to the workshop execution interface, while simultaneously copying it to the solution output and visualization module for on-site display.

[0158] The disturbance monitoring and impact assessment module communicates with the disturbance alarm interface and programmable logic controller via industrial Ethernet, forming the system's entry point on the disturbance event dimension. This module further includes, in sequence, a disturbance event acquisition unit, a node delay mapping unit, a disjunction graph right shift propagation simulation unit, and a disturbance severity calculation unit. The disturbance event acquisition unit receives three types of events pushed by the disturbance alarm interface: equipment failure, urgent order insertion, and processing time deviation. The node delay mapping unit maps these events to equivalent node delays or new node sets on the disjunction graph according to the event type. The other input of the disjunction graph right shift propagation simulation unit is connected to the output of the multi-objective optimization engine module, used to read the disjunction graph structure data of the currently executing production scheduling plan and perform a right shift propagation simulation by forward traversing from the source node to the sink node in topological order. The disturbance severity calculation unit calculates the disturbance severity of disturbance event e according to Formula 6. The disturbance monitoring and impact assessment module is then output via the system's internal data bus and is used as the module's output. This corresponds to the hardware implementation of step S5 in the method.

[0159] The input of the adaptive rescheduling decision module receives the disturbance severity via the system's internal data bus. This module further includes a threshold matching unit and a right-shift repair execution unit, a critical path local re-optimization execution unit, and a global rolling reordering execution unit, which are selectively connected to the output of the threshold matching unit. Another input of the threshold matching unit is connected to the output of the perturbation statistical sample knowledge base module, used to read the mild threshold and the severe threshold. The selective connection means that the threshold matching unit enables one corresponding execution unit and closes the other two based on its judgment result, avoiding concurrent execution of the three response strategies and causing scheme conflicts. The right-shift repair execution unit is implemented independently in hardware by this module; the critical path local re-optimization execution unit and the global rolling reordering execution unit are connected to the multi-objective optimization engine module through a feedback path in hardware, reusing its computing resources. The output of the adaptive rescheduling decision module is divided into three paths: the first path connects to the scheme output and the second input of the visualization module, sending the updated scheme out for visualization; the second path feeds back to the feedback input of the multi-objective optimization engine module, triggering critical path local re-optimization or global rolling reordering; the third path connects to the return input of the perturbation statistical sample knowledge base module, returning the current response sample to the knowledge base. The adaptive rescheduling decision module corresponds to the hardware implementation of method step S6.

[0160] The disturbance statistical sample knowledge base module further includes a disturbance sample storage unit, a processing time statistical parameter calculation unit, and a threshold automatic update unit. The disturbance sample storage unit is responsible for storing and indexing the samples; the processing time statistical parameter calculation unit updates the corresponding processing time standard deviation according to the deviation between the actual processing time of each process on different equipment and the standard processing time, based on the workpiece-process-equipment combination. The automatic threshold update unit updates the mild and severe thresholds based on the preset quantiles of the relative increment distribution of the maximum completion time actually generated by the right-shift repair strategy in the samples. The module's external interface includes three outputs: providing the processing time standard deviation to the robust buffer embedding unit in the multi-objective optimization engine module; providing the mild and severe thresholds to the threshold matching unit in the adaptive rescheduling decision module; and providing the identification information required for the decision report to the scheme output and visualization module. It also includes one feedback input, receiving sample feedback from the adaptive rescheduling decision module. The perturbation statistics sample knowledge base module corresponds to the hardware implementation of the knowledge base self-maintenance mechanism in method step S72.

[0161] The solution output and visualization module further includes a Gantt chart generation unit, a decision report generation unit, and a visualization distribution unit. The Gantt chart generation unit generates a two-dimensional Gantt chart in the workpiece-equipment dimension based on the initial execution plan or updated plan. The decision report generation unit generates a decision report containing three target values, critical path identifiers, and buffer embedding position identifiers. The visualization distribution unit's output is connected to the workshop visualization equipment and the operator's tablet terminal via an industrial Ethernet connection, simultaneously distributing the Gantt chart and decision report to the large screen and the terminal. The workshop visualization equipment is typically a display screen in the workshop office area or near the workstation, facilitating all personnel to grasp the current scheduling status. The operator's tablet terminal is a portable terminal held by the on-site operator, enabling them to precisely grasp the execution rhythm based on the critical path and buffer embedding position identified in the decision report.

[0162] The system provided by this invention, at the hardware implementation level, can utilize an industrial-grade server or industrial control computer deployed on the workshop side to carry the computing units of steps S2 to S6. It can also use fieldbus or industrial Ethernet to form the external communication channels for the data acquisition and preprocessing module and the disturbance monitoring and impact assessment module. Furthermore, it can use a workshop large screen and operator tablet terminals to form the physical terminals for the solution output and visualization module. The input and output terminals between each module are connected in series via the system's internal data bus, while the external interfaces are implemented according to industrial Ethernet and fieldbus protocol specifications. This allows the system to be easily embedded into existing workshop information systems without significant modifications to existing upstream systems such as manufacturing execution systems and enterprise resource management systems.

[0163] Example 1 applies the method of the present invention to a flexible workshop scheduling scenario in a machinery manufacturing enterprise. This workshop undertakes mixed-flow production of multiple product types and small batches for machining. In this example, a planned 8-hour day shift is selected, with 10 workpieces to be processed, 6 machines, and a total of 55 processes. The workpiece release time covers the period from 0 to 120 minutes from the start of the shift, and the standard processing time for a single workpiece process is distributed between 10 and 45 minutes. The workpiece, equipment, and processing time matrix used in this example is based on the MK02 standard benchmark example.

[0164] In step S1 of this invention, the manufacturing execution system, the equipment IoT interface, and the material management system collaboratively collect production scheduling input information. In step S11, the workpiece number, the process number, the set of processable equipment for each process, the standard processing time for each process on each processable equipment, the release time, delivery date, and delay penalty weight for each workpiece are read. The current availability status, remaining processing time, and expected release time of the six equipment units are also read, along with the current inventory and available time windows of metal raw materials, cutting tools, and fixtures. In step S12, the multi-source data is mapped to the memory data structure of the production scheduling solver according to preset fields. For the very few missing standard processing time fields, the average value of similar processes is used to fill in the gaps. Fields with inconsistent dimensions are unified as min and piece, resulting in the production scheduling input dataset of this embodiment.

[0165] In step S2, a multi-objective flexible workshop scheduling model is constructed. In step S21, the workpiece is used as the basis for scheduling. The Is the process assigned to the first... Binary decision variables for processing with machine tools and the actual start time of this process These are the decision variables. In step S22, the three optimization objectives are to minimize the maximum completion time, minimize the weighted delay, and maximize the equipment utilization rate. In step S23, the six hard constraints are: process sequence constraint, unique allocation constraint of process equipment, mutual exclusion constraint of processing on the same equipment, workpiece release time constraint, equipment machinability constraint, and material availability constraint.

[0166] In step S3, a dual-file co-evolutionary non-dominated sorting genetic algorithm is executed. In step S31, a dual-layer chromosome is used as the individual encoding, with a process sorting segment length of 55 and a device allocation segment length of 55. In step S32, an initial population is generated using a hybrid strategy, with a population size of 100, and 60% of the initial population is globally selected as initial individuals. In step S33, priority processes with a crossover rate of 0.85 are used to maintain crossover and perform positional uniform crossover, while random exchange mutation of the process sorting segment and directed mutation of the critical path device allocation segment are performed with a mutation rate of 0.10.

[0167] In step S34, the earliest start time and earliest finish time of each process are calculated by traversing the topology from the source node to the sink node in a forward direction. Then, using the global maximum finish time as the endpoint, the latest start time and latest finish time of each process are calculated by traversing the topology from the sink node to the source node in a reverse direction. Processes whose difference between the latest start time and the earliest start time is zero are identified as the critical path process set CP. In step S35, each critical process in the process set CP is... Read the processing time standard deviation from the perturbation statistical sample knowledge base. Calculate the robust buffer length using Formula 1. And before embedding the start time of the corresponding process, update the earliest start time and earliest completion time of the affected nodes along the disjunction graph; in this embodiment, the robustness factor is... Set the value to 1.0, and set the maximum number of iterations to 3.

[0168] In step S36, the set obtained by merging the main evolutionary population with the current generation's offspring is divided into several non-dominated levels according to the non-dominated sorting. The non-dominated levels are then filled into the convergence file in order from low to high, and within the same non-dominated level, they are filled in order from large to small crowding distance. The convergence file size is set to 100. In step S37, the relative dispersion file for each non-dominated solution is calculated according to Formula 2. Riesz2-Energy Contribution Value ,according to Fill in the scattered files in ascending order of size. The capacity of the distributed archives is set to 100, with zero-distance protection. Pick Elite migration is performed every 10 generations; the maximum number of evolution generations is 200. After the iteration terminates, [the process will proceed as follows]. and By taking the union of the solutions and removing inferior and duplicate solutions dominated by other solutions, we obtain the Pareto non-dominated solution set in this embodiment.

[0169] In step S4, the entropy weight-approximation ideal solution ranking method is executed for decision-making. In step S41, the Pareto non-dominated solution set is organized into a decision matrix; in step S42, the objective entropy weight of each objective is calculated according to formulas 3 and 4; in step S43, the relative proximity of each candidate solution is calculated according to formula 5. The candidate scheme with the highest relative proximity value is selected as the initial execution scheme and sent to the workshop execution interface.

[0170] In step S5, workshop disturbance events are continuously monitored. In step S51, three types of disturbance events—equipment failure, urgent order insertion, and processing time deviation—are received through the disturbance alarm interface, and the disjunction graph corresponding to the currently executing production scheduling plan is processed according to the event type. In step S52, a right-shift propagation simulation is performed, traversing forward from the source node to the sink node in topological order. In step S53, disturbance events are quantified according to formula 6. Severity of disturbance .

[0171] In step S6, hierarchical rescheduling is performed. In step S61, the disturbance severity is... With mild threshold and severe threshold Comparison; This embodiment The initial value is 0.10. The initial value is set to 0.30, and the preset sample threshold is set to 40. In steps S62, S63, and S64, one of the three response strategies—right shift repair, local re-optimization of the critical path, and global rolling rearrangement—is implemented to obtain the updated solution.

[0172] In step S7, an update plan and decision report are output, and sample backflow is performed. In step S71, the update plan is output as a two-dimensional Gantt chart in the workpiece-equipment dimension to the workshop visualization equipment and operator tablet terminals, along with a decision report containing three target values, critical path identifiers, and buffer embedding location identifiers. In step S72, the characteristic parameters of this disturbance event and the type of response strategy adopted are used as a sample backflow to the disturbance statistical sample knowledge base, and the knowledge base updates the mild threshold based on the accumulated samples. Severe threshold and processing time standard deviation This concludes the implementation process of Example 1.

[0173] Example 2: Workshop environment, problem size, algorithm parameters divided by robustness factor in Example 2 compared to Example 1. The external aspects remain consistent, the difference being the robustness factor. Take 0.5. This embodiment is suitable for scenarios where the workshop has high requirements for tight production cycles and low tolerance for processing time disturbances, corresponding to formula 1. The engineering trade-off range where smaller values ​​result in a more compact solution but reduced robustness.

[0174] Example 3: Workshop environment, problem size, algorithm parameters divided by robustness factor in Example 3 compared to Example 1. The external aspects remain consistent, the difference being the robustness factor. Take 2.0. This embodiment is suitable for scenarios with large disturbances in workshop processing time and high requirements for solution stability. It corresponds to the buffering method in Formula 1, which embeds a buffer at twice the standard deviation. Based on Chebyshev's inequality, the deviation of single-process time from the mean can be reduced. The upper bound of the probability is reduced from 100% to the 25% engineering range.

[0175] Example 4: Workshop environment, problem size, algorithm parameters divided by robustness factor in Example 4 compared to Example 1. The external aspects remain consistent, the difference being the robustness factor. We take 3.0. This embodiment is applicable to scenarios where the processing time of critical processes has serious uncertainties and the highest robustness must be prioritized, corresponding to an engineering range where the upper bound of the Chebyshev inequality probability is further reduced to approximately 11.1%.

[0176] Example 5 applies the method of the present invention to the BrandimarteMK04 benchmark instance and its derived perturbation scenarios. In this example, there are 15 workpieces, 8 machines, and a total of 90 processes. The remaining algorithm parameters are consistent with Example 1, and the robustness factor is... 1.0 is used. This embodiment is used to verify the feasibility and scalability of the present invention in addressing production scheduling issues in medium-sized flexible workshops.

[0177] Comparative Example 1 maintains the same workshop environment, problem size, and problem data as Example 1. The difference lies in the solution phase, which employs the standard non-dominated sorting genetic algorithm NSGA-II proposed by Deb et al. in 2002, without introducing a dual-archive co-evolutionary mechanism or robust buffer embedding for critical processes. The decision-making phase utilizes a manually weighted three-objective summation decision, with weights of 0.33, 0.33, and 0.34 for the three objectives. This comparative example corresponds to the baseline scheme without any of the improved features of this invention.

[0178] Comparative Example 2, compared to Example 1, uses the standard NSGA-II algorithm in the solution stage but introduces the same key step as Example 1: robust buffer embedding operation; it does not introduce the dual-archive co-evolutionary mechanism; and in the decision-making stage, it uses the same entropy weight-approximation ideal solution sorting method as Example 1. This comparative example is used to isolate and evaluate the contribution of the dual-archive co-evolutionary mechanism to the Pareto approximate front quality while retaining robust buffer embedding.

[0179] Comparative Example 3, compared to Example 1, employs the same dual-archive co-evolutionary non-dominated sorting genetic algorithm in the solution phase but without introducing a robust buffer embedding operation for critical steps; the decision-making phase uses the same entropy-weighted approximation of the ideal solution sorting method as in Example 1. This comparative example is used to isolate and evaluate the contribution of robust buffer embedding for critical steps to the robustness of the scheme to perturbations while preserving dual-archive co-evolution.

[0180] Comparative Example 4, compared to Example 1, uses the SPEA2 algorithm as another mainstream multi-objective evolutionary algorithm baseline in the solution stage, without introducing a dual-file co-evolutionary mechanism or robust buffer embedding for key processes; the decision-making stage uses the same entropy weight-approximation of ideal solution sorting method as Example 1. This comparative example is used to demonstrate the relative advantages of the present invention compared to other mainstream multi-objective evolutionary algorithms.

[0181] Experiment 1 uses 10 standard benchmark flexible job shop scheduling instances (Brandimarte MK01 to MK10) as the underlying problem. The method in Example 1, as well as the methods in Comparative Examples 1, 2, 3, and 4, are each evolved 200 times. For each generation, the Hypervolume index value of the Pareto non-dominated solution set relative to the reference point is recorded. The reference point is taken as 1.1 times the worst value of each target in all examples and comparative examples, forming the upper right corner reference point for the 3D Hypervolume calculation. The Hypervolume index is a commonly used convergence metric for measuring the target space volume of the Pareto approximate front envelope. Each algorithm in each instance is run independently 30 times. The average of the 30 generational HV values ​​is used as the representative value of the algorithm for that generation in that instance. The average of the representative values ​​from the 10 benchmark instances is then used to obtain the overall mean of the algorithm for that generation. Simultaneously, the standard deviation of the 30 generational HV values ​​is used as the width of the corresponding error band.

[0182] The experimental results of Example 1 are as follows: Figure 1 As shown. Figure 1 The horizontal axis represents the number of generations, and the vertical axis represents the standardized Hypervolume index value. Figure 1 Five main curves are plotted to represent the intergenerational mean HV values ​​of Example 1 and Comparative Examples 1 to 4, respectively. The shaded areas on both sides of each main curve represent the ±1 standard deviation range of the HV index of the corresponding generation. The main curve corresponding to Example 1 is drawn with a solid red line, Comparative Example 1 with a dashed blue line, Comparative Example 2 with a dotted green line, Comparative Example 3 with a dotted orange line, and Comparative Example 4 with a solid purple line.

[0183] from Figure 1As can be seen, the red curve corresponding to Example 1 entered the rapid rise phase earliest within the first 80 generations, crossing the horizontal line where the Hypervolume index is 0.7 as early as around generation 55. The orange curve corresponding to Comparative Example 3 was next, and the curves corresponding to Comparative Example 2, Comparative Example 4, and Comparative Example 1 crossed the horizontal line later in that order. The final Hypervolume value at the end of 200 generations of evolution differed among the five methods. The average final Hypervolume value of Example 1 was approximately 0.85, that of Comparative Example 3 was approximately 0.82, that of Comparative Example 2 was approximately 0.78, that of Comparative Example 4 was approximately 0.76, and that of Comparative Example 1 was approximately 0.75. The width of the shaded bands on both sides of the five curves differed. The shaded band corresponding to Example 1 was the narrowest, and that of Comparative Example 1 was the widest, indicating that the variance of Example 1 under 30 independent runs was significantly smaller than that of the comparative examples.

[0184] The method of this invention outperforms the four comparative methods in terms of Pareto approximate front quality. The final Hypervolume improvement compared to Comparative Example 1 is approximately 10 percentage points, compared to Comparative Example 3 approximately 3 percentage points, compared to Comparative Example 2 approximately 7 percentage points, and compared to Comparative Example 4 approximately 9 percentage points. Introducing only the dual-archive co-evolution mechanism brings about a 7 percentage point improvement in Hypervolume for Comparative Example 3 compared to Comparative Example 1, and further introducing robust buffer embedding for key processes brings about another 3 percentage point improvement. The method of this invention exhibits the smallest variance in multiple independent repeated experiments, demonstrating the stability of the solution process. These conclusions prove that both the robust buffer embedding for key processes in step S35 and the dual-archive co-evolution in steps S36 to S37 of this invention substantially contribute to the Pareto approximate front quality, and the two have a synergistic effect that can be superimposed.

[0185] Dispersed archives Based on the Riesz2 energy contribution value defined in Formula 2 The solution set is filtered by value from smallest to largest, so that the non-dominated solutions retained in the dispersed archive tend to be located at the support points in the front sparse region. This avoids the boundary sparsity phenomenon caused by crowding distance failure in the traditional NSGA-II in a 3-object space, thus improving the Hypervolume index on the front edge volume; converged archive The elite solutions selected from low to high non-dominated levels are retained, and elite migration is performed every 10 generations to synergistically improve the convergence and dispersion of the solution set. At the same time, the robust buffer of the key process is embedded in the decoding stage to provide each individual with a buffer proportional to the historical working time variance, so that the decoded solution can improve the absorption of first-order random fluctuations without significantly increasing the maximum completion time, and converge to the region with higher frontier quality in the search process.

[0186] Experiment 2: In Experiment 2, the method of Example 1 and the methods of Comparative Examples 1 to 4 were each run independently 30 times on the Brandimarte MK02 benchmark instance, with 200 generations per run. The Hypervolume index value at termination was recorded. A paired significance test was performed using the Wilcoxon rank-sum test on the Hypervolume index distributions of Example 1 and each comparative example. The significance level was set at... We set the value to 0.05. The Wilcoxon rank-sum test is a nonparametric two-sample rank test that makes relaxed assumptions about the shape of the population distribution and is suitable for comparing the significance of repeated experimental data in multi-objective evolutionary algorithms.

[0187] The experimental results of Example 2 are as follows Figure 2 As shown. Figure 2 The horizontal axis, from left to right, represents Example 1, Comparative Example 1, Comparative Example 2, Comparative Example 3, and Comparative Example 4, while the vertical axis represents the Hypervolume index value at the end. The bottom and top edges of each box plot represent the 25th and 75th percentiles, respectively. The center line of the box represents the median, the whiskers represent extreme values ​​within 1.5 interquartile ranges, and the plus signs represent outliers. The box corresponding to Example 1 is filled with red, Comparative Example 1 with blue, Comparative Example 2 with green, Comparative Example 3 with orange, and Comparative Example 4 with purple. Figure 2 The Wilcoxon rank-sum test results of Example 1 and four comparative examples are listed below. Value; considering the floating-point display precision of statistical software, extremely small. The values ​​in the text box are uniformly set to... The formal expression, each The value is followed by an asterisk indicating the significance level, where three asterisks indicate... The level of highly significant difference is less than 0.001.

[0188] from Figure 2 As can be seen, the box position corresponding to Example 1 is the highest, with a median of approximately 0.852. Comparative Example 3 is next, at approximately 0.821, Comparative Example 2 at approximately 0.778, Comparative Example 4 at approximately 0.765, and Comparative Example 1 is the lowest, at approximately 0.750. The box span, i.e., the interquartile range, corresponding to Example 1 is the smallest, reflecting the smallest variance across 30 independent runs. The four values ​​obtained from the Wilcoxon rank-sum test... The values ​​are all less than 0.0001 and are marked with three asterisks, indicating that the Hypervolume distribution of Example 1 is statistically significant relative to the four comparative examples at a significance level of 0.05.

[0189] The present invention demonstrates a statistically significant advantage over four comparative Hypervolumes in 30 independent runs, indicating that the advantage does not depend on a random initial population or stochastic process, but rather stems from a stable improvement in the algorithm mechanism itself. The minimum span of the corresponding box proves that the method of the present invention has the highest solution stability, facilitating repeatable solution results during engineering implementation.

[0190] Dispersed files in each independent run Filter by formula 2 and then... Value sorting, convergence archive Sort by non-dominated level and crowding distance, so that each independent run can stably retain two types of solutions covering the sparse region and elite region of the frontier from the merged set; the elite migration mechanism executed every 10 generations further reduces the trajectory divergence under different initial populations, so that the final frontier quality remains compactly distributed under multiple repetitions; this mechanism is the algorithmic basis for the significance results obtained by the Wilcoxon rank sum test.

[0191] Experimental Example 3 uses the robustness multiples corresponding to Examples 1, 2, 3, and 4, respectively. Choose from four values: 1.0, 0.5, 2.0, and 3.0, and supplement them. Take the middle values ​​of 0.25, 0.75, and 1.5, for a total of 7 values. Value points. The corresponding initial execution schemes are solved on the BrandimarteMK02 benchmark instance, and each value point is recorded. Deterministic maximum completion time of the solution under the given value Then, 1000 perturbation Monte Carlo simulations were performed on each scheme. The perturbation method was to add a mean of 0 and a standard deviation of 0 to the standard processing time of each process on the equipment allocated to that process. The Gaussian random perturbation is used to determine the proportion of samples in which the original process sequence remains feasible under the perturbation, which is then taken as the perturbation absorption rate. Simultaneously, Chebyshev's inequality is applied... Give the corresponding Theoretical perturbation delay lower bound for the value ;when When the value is less than or equal to 1, the lower bound of the theory is 0.

[0192] The experimental results of Example 3 are as follows: Figure 3 As shown. Figure 3 The horizontal axis represents the robustness factor. The value of is given by the left vertical axis, which represents the deterministic maximum completion time. The unit is min; the right vertical axis represents the perturbation absorption rate and the lower bound of Chebyshev's theory, which is dimensionless and takes values ​​between 0 and 1. Figure 3 The curve formed by the solid red line connecting the dots corresponds to the left vertical axis. Follow The trend of change is shown by the blue solid line connecting the squares, which corresponds to the measured perturbation absorption rate on the right vertical axis, and the green dashed line connecting the triangles, which corresponds to the lower bound of Chebyshev's theory on the right vertical axis. . Figure 3 Examples 2, 1, 3, and 4 are respectively in Five-pointed stars with gold filling and black borders at positions 0.5, 1.0, 2.0, and 3.0 are superimposed on the measured perturbation absorption rate curve. Each five-pointed star is accompanied by a label with a white background and black border indicating the embodiment number. Item 4 of the legend clearly identifies the five-pointed star markings corresponding to Embodiments 1 to 4. The position on the axis makes it easy to read the location of the four embodiments. The values ​​and corresponding perturbation absorption rates.

[0193] from Figure 3 It can be seen from this that Curve follows Monotonous rise, from At approximately 0.25, the value rose to approximately 196 minutes. At approximately 248 minutes with a time of 3.0, the overall increase was about 27%; the measured disturbance absorption rate increased with... Monotonous rise, from When taken at 0.25, it rises to approximately 0.44. The value is approximately 0.98 at 3.0; the measured disturbance absorption rate is consistently higher than the lower bound of Chebyshev's theory, consistent with the mathematical fact that the actual exceedance rate is less than the theoretical upper bound given by Chebyshev's inequality when the working time distribution is approximately normal; this can be read from the positions of the four pentagram markers. When taking 1.0 The increase is about 5% while the disturbance absorption rate is about 0.74, which is in the cost-benefit equilibrium region, corresponding to Example 1; When taken as 0.5 The elevation is approximately 3% while the disturbance absorption rate is approximately 0.58, corresponding to Example 2 for a compact scenario; When using version 2.0 The elevation is approximately 9% while the disturbance absorption rate reaches approximately 0.93, corresponding to Example 3 of a moderately robust scenario; When taking version 3.0 The elevation is about 15% while the disturbance absorption rate is about 0.98, corresponding to Example 4 of a highly robust scenario.

[0194] Robustness multiple Taking 1.0 as the preferred value will help maintain as much as possible. Achieving high disturbance absorption capacity while maintaining compactness; when robustness is a priority in the workshop, the size can be appropriately increased. The robust buffer length defined in Formula 1 and and The proportional design is adjustable during engineering implementation, making it easy to select the appropriate design based on the actual disturbance characteristics of the workshop. Values; the relative positions of the measured perturbation absorption rate and the lower bound of Chebyshev's theory are within 7. The consistent value across all points verifies the correctness of the design basis of Formula 1, which is based on Chebyshev's inequality.

[0195] Formula 1 will increase the robustness factor. Standard deviation of processing time The product of these factors, embedded as a buffer length before the critical process, is equivalent to pre-absorbing the first-order random fluctuations in the processing time of the critical process with a controllable probability. This probability is guaranteed by a strict upper bound given by Chebyshev's inequality; when the working time distribution shape is close to normal, the actual deviation from the mean exceeds... The probability of a disturbance being twice the standard deviation is less than the Chebyshev upper bound, which makes the measured perturbation absorption rate higher than the theoretical lower bound.

[0196] Experiment 4 uses five standard benchmark instances (Brandimarte MK01 to MK05) as the problem's underlying layer, solving them using the methods of Example 1, Comparative Examples 1, 3, and 4 respectively. Comparative Example 2 is omitted to reduce graph density because its comparison dimension partially overlaps with Comparative Example 3. Each instance is run independently 30 times using each method, and the mean Hypervolume index and the mean Inverse Generation Distance (IGD) index are recorded for each of the 30 runs. The IGD index is defined as the mean distance from each solution to the approximate front on the reference front; a smaller IGD value indicates a higher accuracy in approximating the reference front. The reference front in this experiment is obtained by combining the non-dominated solution sets obtained from 30 independent runs of each instance using the four methods, and then removing dominated solutions to obtain the union front.

[0197] The experimental results of Example 4 are as follows: Figure 4 As shown. Figure 4 It consists of two images, (a) and (b). Figure 4 In the middle (a), there is a bar chart of Hypervolume index grouping for four methods under five benchmark instances. The horizontal axis represents benchmark instances MK01 to MK05, and the vertical axis represents Hypervolume index values. The four colored bars represent Example 1, Comparative Example 1, Comparative Example 3, and Comparative Example 4, respectively. Figure 4 (b) is a grouped bar chart of the IGD indicator under the same layout, with the vertical axis representing the IGD indicator value. The smaller the IGD value, the better. Figure 4 The red bars represent Example 1, the blue bars represent Comparative Example 1, the orange bars represent Comparative Example 3, and the purple bars represent Comparative Example 4.

[0198] from Figure 4 It can be seen from this that Figure 4 In Figure (a), the red bar corresponding to Example 1 ranks first in all five benchmark examples. The relative improvement of Hypervolume compared to Comparative Example 1 is between 8% and 15%, the improvement compared to Comparative Example 3 is about 2% to 5%, and the improvement compared to Comparative Example 4 is about 8% to 14%. Figure 4 (b) shows that the red bar corresponding to Example 1 is the lowest among the five benchmark examples, with a reduction in IGD of 30% to 45% relative to Comparative Example 1, a reduction of about 15% to 25% relative to Comparative Example 3, and a reduction of about 30% to 40% relative to Comparative Example 4. From MK01 to MK05, as the example size gradually increases, the advantage of Example 1 relative to Comparative Example 1 and Comparative Example 4 generally expands, while the advantage relative to Comparative Example 3 remains stable.

[0199] The method of this invention outperforms four main comparative methods in comprehensive performance across multiple benchmark instances. The Hypervolume metric and the IGD metric jointly verify that the volume envelope and approximation accuracy of the Pareto approximation front of this invention are both superior. The two key techniques, robust buffer embedding and dual-file co-evolution, are applied to different stages of the solution process. The former applies to the decoding stage, and the latter applies to the selection stage. The synergy of the two techniques enhances the solver in both key aspects of the solution space search. The Hypervolume metric measures the target space volume of the Pareto front envelope, while the IGD metric measures the approximation accuracy between the approximation front and the reference front. The two techniques reflect the front quality from two complementary dimensions: volume and accuracy.

[0200] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A multi-objective dynamic scheduling decision-making method, characterized in that, Includes the following steps: S1: Collect and preprocess order, equipment, material, and real-time status data; S2: Construct a multi-objective flexible workshop scheduling model; S3: Solve the model using a dual-file co-evolutionary non-dominated sorting genetic algorithm. During decoding, a robust buffer is embedded before the critical path process to obtain the Pareto non-dominated solution set. S4: The initial execution plan is obtained by using the entropy weight-approximation ideal solution sorting method on the Pareto non-dominated solution set; S5: Perform right-shift propagation simulation based on the disjunction graph corresponding to the current production schedule, and quantify the severity of the disturbance; S6: Perform right-shift repair, local re-optimization of the critical path, or global rolling rearrangement according to the severity of the disturbance to obtain an updated solution; S7: Output the update scheme and decision report, and return the disturbance sample to the disturbance statistical sample knowledge base.

2. The multi-objective dynamic scheduling decision-making method according to claim 1, characterized in that, Step S1 includes: S11: Read the workpiece number, process number, set of processable equipment for each process, standard processing time for each process on each processable equipment, release time of each workpiece, delivery date and lateness penalty weight of each workpiece through the manufacturing execution system; read the current availability status, remaining processing time and expected release time of each equipment from the programmable logic controller through the equipment IoT interface; read the current inventory and available time window of each material through the material management system; S12: Map the data collected in step S11 to the memory data structure of the production scheduling solver according to preset fields; fill missing standard processing time, release time, remaining processing time, and inventory quantity numerical fields with historical averages or averages of similar processes; perform integrity checks on missing workpiece number fields, process sequence fields, and processable equipment set fields, and mark orders that fail the integrity check as pending manual confirmation; perform unified unit conversion on each field, unifying time fields to minutes and quantity fields according to the corresponding material's unit of measurement; obtain a production scheduling input dataset including workpiece set, equipment set, process set, processing time matrix, and constraint parameter list, as input for subsequent steps.

3. The multi-objective dynamic scheduling decision-making method according to claim 1, characterized in that: The multi-objective flexible workshop scheduling model constructed in step S2 uses a binary decision variable to determine whether the j-th process of workpiece i is assigned to the k-th machine. and the actual start time of this process. The decision variables are defined as follows: minimizing the maximum completion time, minimizing the weighted delay, and maximizing equipment utilization. The optimization objectives are: process sequence constraints, unique equipment allocation constraints for processes, mutual exclusion constraints for processing on the same equipment, workpiece release time constraints, equipment machinability constraints, and material availability constraints. A value of 1 indicates allocation, and a value of 0 indicates no allocation. Let i be the actual start time of the j-th process of workpiece i; i is the workpiece number; j is the process number within the same workpiece; k is the equipment number; the weighted delay is the sum of the products of the delay value of each workpiece whose actual completion time is later than its delivery time and the corresponding delay penalty weight; the equipment utilization rate is the average of the proportion of the effective processing time of each equipment to the available time of that equipment. The encoding and evolutionary operations of the dual-file co-evolutionary non-dominated sorting genetic algorithm in step S3 include: S31: Using a double-layered chromosome as the individual code, the first segment is the process sorting segment, the length of which is the sum of the number of processes for each workpiece. Each gene position is filled with the corresponding workpiece number, and the number of times the same workpiece number appears in this segment is equal to the number of processes for that workpiece. The second segment is the equipment allocation segment, the length of which is the same as the first segment. Each gene position is filled with the sequence number of the equipment selected by the process corresponding to that position in its set of processable equipment. S32: A hybrid strategy is used to generate the initial population. Some individuals in the initial population are generated by the global selection initialization method, while the remaining individuals are generated in a completely random manner. S33: For the parent individual's process sorting segment, priority processes are maintained with crossover; for the equipment allocation segment, positional uniform crossover is adopted. Mutation is performed at a preset mutation rate. The process sorting segment mutation adopts the method of randomly swapping the positions of two different workpiece numbers. The equipment allocation segment mutation is performed on the gene bit corresponding to the critical path process identified by the decoding in step S34: when there are no less than two machines in the set of processable machines, the gene bit is changed to the sequence number of the candidate machine in the set of processable machines whose standard processing time is less than that of the currently selected machine; otherwise, the equipment allocation segment mutation is not performed.

4. The multi-objective dynamic scheduling decision-making method according to claim 3, characterized in that, Step S3 further includes: S34: Based on the individual's two-layer chromosome decoding, a preliminary production scheduling plan is obtained, and a directed acyclic disjunctive graph corresponding to the preliminary production scheduling plan is constructed; the earliest start time and earliest finish time of each process are calculated by forward traversal from the source node to the sink node in topological order; then, with the global maximum finish time of the preliminary production scheduling plan as the endpoint, the latest start time and latest finish time of each process are calculated by reverse traversal from the sink node to the source node in topological reverse order; the process with the difference between the latest start time and the earliest start time being zero is identified as the critical path process, and the corresponding process set is denoted as CP; S35: For each critical process in the process set CP The key process is read from the disturbance statistical sample knowledge base. Standard deviation of processing time The key process is calculated according to the following formula. The robust buffer length that should be embedded : ； In the formula: Let be the robust buffer length that should be embedded before the j-th process of workpiece i when it is processed on equipment k. It is a robust multiple, dimensionless, and a positive real number; Let be the standard deviation of the processing time of the j-th operation on workpiece i on equipment k; i is the workpiece number; j is the operation number within the same workpiece; k is the equipment number. The j-th process of workpiece i is assigned to the process node processed by equipment k in the preliminary production scheduling plan; The obtained robust buffer length Before the start time of the corresponding process is inserted as an additional delay, the earliest start time and earliest completion time of the direct successor node of the same workpiece, the direct successor node of the same equipment, and the recursive successor node of the aforementioned successor node are updated along the disjunction graph. The updated scheme is used as the input of step S34 again, and steps S34 and S35 are repeated until the critical path process set obtained by two consecutive iterations is the same or the number of iterations reaches the preset upper limit. The feasible production scheduling scheme corresponding to this individual is then output.

5. The multi-objective dynamic scheduling decision-making method according to claim 3, characterized in that, Step S3 further includes: S36: During the selection phase of each generation, the set obtained by merging the main evolutionary population with the current generation's offspring is divided into several non-dominant levels according to the non-dominant ranking. The non-dominant levels are then filled into the convergence file in order of increasing non-dominant level and decreasing crowding distance within the same non-dominant level. Until the converged file is reached. Fill the preset capacity; S37: Normalize the objective value of each non-dominated solution in the merged set using the objective range to obtain a normalized objective vector; when the files are distributed... When the set is empty, the non-dominated solution with the largest Euclidean distance between the normalized target vector and the current ideal point in the merged set is taken as the dispersed file. The initial solution; calculate the non-dominated solutions in the merged set according to the following formula. Relative to the scattered archives Riesz2-Energy Contribution Value : ； In the formula: Non-dominated solution Relative to the scattered archives The Riesz2- energy contribution value is dimensionless. Non-dominated solution With the aforementioned scattered archives China-Africa Domination Solution The Euclidean distance over the normalized target vector is dimensionless. It is a preset minimum positive real number used to prevent division by zero when the distance is zero; Indicates taking and The larger of the two; For the distributed archive; i is the number of the non-dominated solution to be evaluated; l is the distributed archive. The numbers of the non-dominated solutions already included in the database; This indicates the location of the scattered archives. Unless the solution is dominated Summing the remaining non-dominated solutions; according to The remaining non-dominated solutions in the merged set are filled into the dispersed file in ascending order. until the scattered archives Fill the preset capacity; every preset exchange interval algebra, from the converged file The non-dominated solution with the largest current congestion distance is selected for migration into the dispersed archive. and from the scattered archives The non-dominated solution with the smallest current non-dominated level value is selected and migrated into the convergence archive. The converged archive is then moved in. Non-dominated solutions from the scattered archives After removing from the middle, press From smallest to largest, select an undoed non-dominated solution from the merged set and add it to the scattered file. After the iteration terminates, the converged file will be... With the aforementioned scattered archives Take the union of the solutions, remove inferior solutions dominated by other solutions and duplicate solutions, and obtain the Pareto non-dominated solution set.

6. The multi-objective dynamic scheduling decision-making method according to claim 1, characterized in that, Step S4 includes: S41: Arrange the values ​​of each candidate solution in the Pareto non-dominated solution set for the three objectives of maximum completion time, weighted delay, and equipment utilization rate by row as candidate solutions and by column as objectives, to obtain an N-row, 3-column decision matrix, where N is the total number of candidate solutions in the Pareto non-dominated solution set; when N=1, directly issue the unique candidate solution in the Pareto non-dominated solution set as the initial execution solution to the workshop execution interface, and step S4 ends; S42: Normalize the maximum completion time column and weighted delay column of the decision matrix according to cost-type indicators, and normalize the equipment utilization rate column according to benefit-type indicators to obtain a normalized decision matrix; when the range of a certain target column is zero, the normalized value of that column is uniformly set to a preset constant of 0.5; calculate the Shannon entropy of the q-th column of the normalized decision matrix according to the following formula. : ； In the formula: Let be the Shannon entropy of the q-th column of the normalized decision matrix, which is dimensionless; q is the target number, and q takes values ​​of 1, 2, and 3, which correspond to the three targets of maximum completion time, weighted delay, and equipment utilization rate, respectively; N is the total number of candidate solutions; The normalized value in the nth row and qth column of the normalized decision matrix is ​​the proportion of the sum of the normalized values ​​in the qth column, which is dimensionless; n is the candidate scheme number. The natural logarithm operator; This represents the summation over n=1 to n=N; when When it is zero, The item is treated as zero; Then calculate the objective entropy weight corresponding to the q-th target using the following formula. : ； In the formula: Let be the objective entropy weight corresponding to the q-th target, which is dimensionless; For the summation index, Values ​​are 1, 2, and 3; For the normalized decision matrix, the first... The Shannon entropy of the column is dimensionless; Indicates to =1 to Summing = 3; S43: Multiply each element in the q-th column of the normalized decision matrix by its corresponding objective entropy weight. The weighted normalized decision matrix is ​​obtained. From each column of the weighted normalized decision matrix, the maximum value is selected to form the positive ideal solution, and the minimum value is selected to form the negative ideal solution. The relative closeness of the nth candidate solution to the positive ideal solution is calculated using the following formula. : ； In the formula: Let be the relative closeness of the nth candidate solution to the positive ideal solution, which is dimensionless; Let be the Euclidean distance between the nth row vector of the weighted normalized decision matrix and the positive ideal solution, which is dimensionless; Let be the Euclidean distance between the nth row vector of the weighted normalized decision matrix and the negative ideal solution, which is dimensionless; The relative closeness The candidate solution with the largest value is sent to the workshop execution interface as the initial execution solution.

7. The multi-objective dynamic scheduling decision-making method according to claim 1, characterized in that, Step S5 includes: S51: Continuously receive disturbance events e through the disturbance alarm interface, and perform one of the following processing on the disjunction diagram corresponding to the currently executing production scheduling plan according to the event type: In the first category, when the disturbance event e is a device failure, the equivalent node of the dissociative graph node corresponding to the currently being processed or the next process to be processed on the faulty device is delayed as the expected fault repair time. The second type is when the disturbance event e is an urgent order insertion type. The urgent order is added to the disjunction graph as a set of new workpiece nodes. The release time of the first process of the new workpiece is set as the urgent order insertion time. The process sequence constraints, equipment manufacturability constraints and material availability constraints of the new workpiece are added to the disjunction graph. The third type is when the disturbance event e is a processing time deviation type. The equivalent node of the disjunction graph node corresponding to the process whose deviation exceeds the preset tolerance is delayed to the difference between the actual processing time and the standard processing time of the process. S52: Using the processed disjunctive graph as input, traverse forward from the source node to the sink node in topological order. For each node: take the new start time of the node as the larger of the maximum value of the new completion times of all its direct predecessor nodes and the node's release time constraint; take the new completion time of the node as the sum of the new start time, standard processing time, and equivalent node delay; after traversal, obtain the updated start time and completion time of each process, and record the maximum value of the updated completion times as... ; S53: Calculate the severity of the disturbance event e using the following formula. : ； In the formula: Let be the severity of the disturbance event e, which is dimensionless; The global maximum completion time obtained after the right-shift propagation simulation in step S52; The global maximum completion time of the current production schedule before the occurrence of the disturbance event e; This indicates taking the larger of a and b; the superscript RS indicates the state after the right-shift propagation simulation; the superscript plan indicates the planned state before the disturbance. The severity of the disturbance This serves as the basis for matching the response strategy in step S6.

8. The multi-objective dynamic scheduling decision-making method according to claim 7, characterized in that, Steps S6 and S7 include: S61: Determine the severity of the disturbance With the slight threshold maintained in the perturbation statistics sample knowledge base and severe threshold Comparison: When Perform rightward repair when Perform local re-optimization of the critical path when A global scroll reordering is performed at that time; among which, The threshold value is light and dimensionless. The severity threshold is dimensionless. When the cumulative number of samples in the disturbance statistics sample knowledge base is less than a preset sample threshold, and Take the preset initial threshold from the disturbance statistics sample knowledge base; S62: The right shift repair maintains the equipment allocation of each process in the current production schedule and the processing order of processes on the same equipment unchanged, and uses the new start time of each process obtained in step S52 to cover the original start time as the output of the updated plan; S63: The critical path local re-optimization uses the right-shift propagation simulation scheme obtained in step S52 as the initial value, identifies its critical path, and unlocks the start time and equipment allocation of the processes on the critical path, the direct successor processes of the same workpiece on the critical path, the direct successor processes of the same equipment on the critical path, and the adjacent processes related to equipment reassignment as decision variables. The start time and equipment allocation of the remaining processes are fixed as determined parameters. The dual-file co-evolutionary non-dominated sorting genetic algorithm of step S3 is called to perform small-scale re-optimization on the reduced subproblem. During the re-optimization process, conflict verification is performed on the fixed process time windows on the reassigned equipment, and the insertionable time window constraint ensures that the processing of the same equipment is mutually exclusive. Then, the Pareto non-dominated solution set obtained by re-optimization is decided according to step S4 to obtain the updated scheme. S64: The global rolling rescheduling fixes the processes that have started but not yet completed at the current time on the originally assigned equipment for continued processing. The start time of the process is fixed to the current time, and the remaining processing time of the process is fixed to the remaining processing time of the process. The processes that have not yet started are set as a new set of processes to be scheduled. The updated equipment status, material status and remaining orders are used as new scheduling inputs. Steps S2, S3 and S4 are re-executed to obtain the updated scheme. S71: Output the update scheme in the form of a two-dimensional Gantt chart in the workpiece-equipment dimension, along with a decision report containing three target values, critical path identifiers, and buffer embedding location identifiers; S72: The characteristic parameters of this disturbance event and the type of response strategy adopted are fed back to the disturbance statistical sample knowledge base as a sample. The characteristic parameters include the disturbance event type, the affected process number, the affected equipment number, the equivalent node delay value, and the disturbance severity. When the cumulative number of samples in the perturbation statistical sample knowledge base reaches the preset sample threshold, the perturbation statistical sample knowledge base updates the mild threshold based on the preset quantiles of the relative increment distribution of the maximum completion time actually generated by the right-shift repair strategy in the cumulative samples. and the severity threshold Based on the deviation between the actual processing time and the standard processing time of each process on different equipment in the cumulative sample, the corresponding standard deviation of the processing time is updated on a rolling basis according to the workpiece-process-equipment combination. .

9. A multi-objective dynamic production scheduling decision-making system, characterized in that, The system is used to implement the multi-objective dynamic scheduling decision-making method according to any one of claims 1 to 8. The system includes a data acquisition and preprocessing module, a multi-objective optimization engine module, an automatic scheme decision-making module, a disturbance monitoring and impact assessment module, an adaptive rescheduling decision-making module, a disturbance statistical sample knowledge base module, and a scheme output and visualization module, wherein: The data acquisition and preprocessing module communicates with the workshop manufacturing execution system, enterprise resource management system, material management system, equipment IoT interface and programmable logic controller via industrial Ethernet; The multi-objective optimization engine module includes a model building submodule, a key process robust decoding submodule, and a dual-file co-evolution submodule connected in sequence; The output of the automatic decision-making module is connected to the first input of the visualization module and the workshop execution interface; The disturbance monitoring and impact assessment module communicates with the disturbance alarm interface and programmable logic controller via industrial Ethernet; The output of the adaptive rescheduling decision module is connected to the second input of the scheme output and visualization module, the feedback input of the multi-objective optimization engine module, and the return input of the perturbation statistical sample knowledge base module, respectively. The perturbation statistics sample knowledge base module provides the processing time standard deviation to the robust decoding submodule of key processes, and provides mild and severe thresholds to the adaptive rescheduling decision module; The output of the solution and visualization module is connected to the workshop visualization equipment and the operator's tablet terminal.

Images