Optimization method for solving mixed integer linear programming problem

Abstract

The invention discloses an optimization method for solving a mixed integer linear programming problem, and the method comprises the steps: employing a bipartite graph structure composed of variable nodes and constraint nodes for the modeling of the MILP problem, and the bipartite graph structure comprises variable node features, constraint node features and edge features; variables of the bipartite graph structure are divided into stable variables and unstable variables, a corresponding total loss function is determined, and the total loss function is determined based on a cross entropy loss function and a comparison loss function; carrying out label learning on the stable variables by adopting a cross entropy loss function, and carrying out comparative learning on the unstable variables by adopting a comparative loss function; based on the total loss function, training learning of the corresponding relation of the MILP problem and the solution is carried out to obtain a corresponding prediction model, and based on the prediction model, the MILP problem is solved. According to the method, prediction precision and solution feasibility can be considered in the MILP problem solving process, and good stability and practicability are shown in various types of actual problems.

Description

Technical Field

[0001] This invention relates to the field of computer artificial intelligence technology, and in particular to an optimization method for solving mixed integer linear programming problems. Background Technology

[0002] Mixed-integer linear programming (MILP) is a mathematical model widely used in optimizing decision-making processes for complex systems. It describes a problem through a set of linear objective functions and constraints, including both integer and continuous variables. Due to its powerful expressiveness and flexibility, MILP has become a key tool for modeling and solving many practical problems. Specifically, MILP has important applications in various fields such as manufacturing scheduling, logistics route optimization, power system management, and financial portfolio allocation.

[0003] However, since MILP problems are inherently nondeterministic polynomial-time (NP) difficult problems, traditional solution methods such as branch and bound, cutting plane methods, linear relaxation techniques, and heuristic search strategies still face significant challenges when running on high-performance computing hardware as the problem size or structural complexity increases. Specifically: on the one hand, as the problem size increases, the computation time of traditional algorithms grows exponentially, making it difficult to obtain high-quality feasible solutions within a reasonable timeframe; on the other hand, in large-scale, dynamically changing real-world application scenarios, the efficiency and stability of traditional methods are limited, failing to meet the dual requirements of real-time performance and accuracy; moreover, even with high-performance computing resources, traditional methods may still fail to achieve efficient solutions due to excessive consumption of computing resources.

[0004] In recent years, with the development of deep learning and graph neural networks, the industry has begun to explore new paradigms that combine machine learning methods with traditional MILP solvers to improve solution efficiency and quality. Among these, the end-to-end variable assignment prediction method can construct a prediction model based on GNNs, modeling the MILP problem as a bipartite graph structure. It utilizes the information transfer mechanism between nodes to capture the complex dependencies between variables and constraints, thereby directly predicting the values ​​of key variables, narrowing the search space, and guiding the solver to converge to a high-quality solution more quickly. However, this method still faces several significant challenges in practical applications that remain unresolved, including:

[0005] (1) Due to the complex interdependencies among variables in MILP problems, if the model makes an incorrect prediction of the key variables, it may lead to the entire solution becoming infeasible or even deviating significantly from the optimal solution.

[0006] (2) The quality and diversity of training data directly affect model performance. Current methods still struggle to train high-performance models in terms of sample construction and feature extraction.

[0007] (3) The existing prediction methods mentioned above use a fixed “trust region” or search range, that is, only allow variables to change within a small range near the initial predicted value. This fixed hyperparameter setting method lacks the ability to adapt to the problem structure and prediction uncertainty. When facing problems of different scales or types, it is easy to find a feasible solution if the search space is too small, or to waste computing resources if it is too large. Moreover, it has a weak ability to transfer across problem types and lacks adaptability to different constraint structures and variable distributions, which limits its widespread application in industrial scenarios.

[0008] In view of this, the present invention is hereby proposed. Summary of the Invention

[0009] The purpose of this invention is to provide an optimization method for solving mixed-integer linear programming problems, so as to balance the prediction accuracy and solution feasibility in the solution process of mixed-integer linear programming problems, and to solve the technical problems existing in the prior art.

[0010] The objective of this invention is achieved through the following technical solution:

[0011] An optimization method for solving mixed-integer linear programming problems includes:

[0012] The modeling of the mixed integer linear programming (MILP) problem adopts a bipartite graph structure composed of variable nodes and constraint nodes. The bipartite graph structure includes variable node features, constraint node features, and edge features.

[0013] Based on the bipartite graph structure, variables are divided into stable variables with consistent values ​​across all solutions and unstable variables with discrepancies, and a corresponding total loss function is determined. The total loss function is based on the cross-entropy loss function and the contrastive loss function. Furthermore, the stable variables are labeled using the cross-entropy loss function, and the unstable variables are labeled using the contrastive loss function.

[0014] Based on the total loss function, the corresponding prediction model is obtained by training and learning the correspondence between the MILP problem and the solution, and the MILP problem is solved based on the prediction model.

[0015] The variable node features include:

[0016] Normalized objective function coefficients are used to reflect the importance of variables in the objective function;

[0017] The average value of the variable's coefficients across all constraints is used to provide the overall impact of the variable under all constraints.

[0018] The degree of a variable node is used to express how closely connected it is in the network;

[0019] The maximum and minimum values ​​of the coefficient of this variable in all constraints are used to describe the upper and lower limits of the variable's scope of influence, respectively.

[0020] Whether it is an integer variable is used to distinguish between continuous variables and integer variables using binary identifiers;

[0021] Position encoding is used to indicate the specific position of a variable in the original sequence.

[0022] The characteristics of the constraint node include:

[0023] The average value of all variable coefficients within the constraint is used to represent the strength of the constraint's influence on the overall variable set;

[0024] The degree of a constraint node indicates the number of variables associated with that constraint;

[0025] The normalized value of the right-hand side term is the constant value of the constraint right-hand side after standardization.

[0026] Constraint direction, used to describe the properties of the constraint condition.

[0027] The edge features include: the specific coefficients of the variables in the linear constraint, which directly reflect the actual magnitude of the variable's contribution to the constraint.

[0028] The calculation method for the total loss function includes:

[0029] ;

[0030] ;

[0031] ;

[0032] in, Let cross-entropy be the loss function. To compare the loss functions, α is the weighting coefficient. It is a variable The true label, The variables predicted by the model The probability of taking the value 1. It's a temperature parameter. It is a positive sample , It is a positive sample set. It is a negative sample , It is a negative sample set.

[0033] The collection process of the positive sample set includes:

[0034] For each MILP instance, a predetermined number of optimal solutions are collected using an advanced solver over a predetermined time period to form a positive sample set.

[0035] The collection process of the negative sample set includes:

[0036] Analyze the solutions in the positive sample set, determine the stable and unstable variables, and randomly perturb a predetermined number of variables in the unstable variables to make the solution infeasible. When a solution is infeasible, it is added to the negative sample set. Each positive sample is used to generate multiple infeasible solutions.

[0037] The random perturbation processing includes: randomly selecting 10% of the unstable variables, flipping the values ​​of the corresponding unstable variables, i.e., changing 0 to 1 and 1 to 0; and each positive sample is used to generate 10 infeasible solutions.

[0038] The training process of the prediction model includes:

[0039] The Adam optimizer was trained using adaptive moment estimation with a learning rate of 0.001, a batch size of 8, and 1000 epochs. An early stopping strategy was employed to prevent overfitting. This strategy involves dynamically monitoring the model performance during training and terminating training prematurely when preset conditions are met but the required number of epochs has not been reached.

[0040] The method also includes an adaptive search region adjustment process, which includes:

[0041] The MILP solver is used to verify whether the solution within the neighborhood is feasible. If the solution is not feasible, a heuristic algorithm is applied to extract the irreducible and infeasible subsystem (IIS) and identify the subset of constraints that cause the conflict. If the solution is feasible, the solver is used to continue solving the remaining variables.

[0042] Based on the constraints extracted from IIS, identify the key variables, remove their fixed state, and reconstruct and adjust the search area.

[0043] If adjusting the search area still does not work, the extracted IIS will continue to be identified and the set of unfixed variables will be gradually expanded according to the set strategy until the predicted solution is a feasible solution.

[0044] Compared with the prior art, the present invention provides an optimization method for solving mixed-integer linear programming problems, which can balance prediction accuracy and solution feasibility in the process of solving mixed-integer linear programming problems. Specifically, it can dynamically adjust the search area while maintaining efficient solution, provide tolerance for prediction errors, and show good stability and practicality in various types of practical problems. Attached Figure Description

[0045] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0046] Figure 1 This is a schematic diagram of the hybrid training strategy framework provided in an embodiment of the present invention;

[0047] Figure 2 This is a schematic diagram of the adaptive search region adjustment framework provided in an embodiment of the present invention. Detailed Implementation

[0048] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the specific content of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments, which do not constitute a limitation of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0049] First, the following explanations are provided for the terms that may be used in this article:

[0050] The term "and / or" means that either or both can be achieved simultaneously. For example, X and / or Y means that it includes both "X" or "Y" as well as the three cases of "X and Y".

[0051] The terms "comprising," "including," "containing," "having," or other similar semantic descriptions should be interpreted as non-exclusive inclusion. For example, including a technical feature element (such as raw material, component, ingredient, carrier, dosage form, material, size, part, component, mechanism, device, step, process, method, reaction conditions, processing conditions, parameter, algorithm, signal, data, product or article of manufacture, etc.) should be interpreted as including not only the expressly listed technical feature element, but also other technical feature elements that are not expressly listed and are well-known in the art.

[0052] The term "composed of" excludes any technical features not expressly listed. When used in a claim, it closes the claim to exclude all technical features other than those expressly listed, except for associated conventional impurities. If the term appears only in a clause of a claim, it limits the claim to the elements expressly listed in that clause; elements recited in other clauses are not excluded from the overall claim.

[0053] Unless otherwise explicitly specified or limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; they can refer to mechanical connections or electrical connections; they can refer to direct connections or indirect connections through an intermediate medium; and they can refer to the internal connection between two components. Those skilled in the art can understand the specific meaning of the above terms in this document according to the specific circumstances.

[0054] When concentration, temperature, pressure, size, or other parameters are expressed as numerical ranges, such ranges should be understood to specifically disclose all ranges formed by any pairing of upper limits, lower limits, or preferred values ​​within that range, regardless of whether the range is explicitly stated; for example, if the numerical range "2 to 8" is stated, then that range should be interpreted to include ranges such as "2 to 7", "2 to 6", "5 to 7", "3 to 4 and 6 to 7", "3 to 5 and 7", "2 and 5 to 7", etc. Unless otherwise stated, the numerical ranges described herein include both their endpoints and all integers and fractions within that range.

[0055] The terms “center,” “longitudinal,” “lateral,” “length,” “width,” “thickness,” “upper,” “lower,” “front,” “back,” “left,” “right,” “vertical,” “horizontal,” “top,” “bottom,” “inner,” “outer,” “clockwise,” and “counterclockwise” indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience and simplification of description and do not imply that the device or component referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this document.

[0056] This invention, based on the intersection of artificial intelligence and operations research, provides a predictive optimization method based on Mixed Integer Linear Programming (MILP). Specifically, it employs a hybrid training strategy based on graph neural networks, combining contrastive learning and labeled supervised learning, and a method for dynamically adjusting the search space using irreducible infeasible subsystems. In particular, it provides a method and system for solving MILP that integrates a hybrid training strategy and an adaptive search region adjustment mechanism, which can effectively improve prediction accuracy and solution feasibility, and enhance the algorithm's robustness and generalization ability on complex problems.

[0057] Specifically, the optimization method for solving mixed-integer linear programming problems provided in this embodiment of the invention aims to at least improve the model's prediction accuracy and ensure the feasibility of the predicted solution. Specifically, this embodiment of the invention can significantly improve the model's prediction quality and the feasibility of the solution by introducing a structure-aware hybrid training strategy and a dynamic search region adjustment mechanism based on Irreducible Infeasible Subsystem (IIS) analysis, thereby enhancing the algorithm's practicality and robustness in complex industrial scenarios.

[0058] Furthermore, the optimization method for solving mixed-integer linear programming problems provided in this embodiment of the invention specifically employs a structure-aware hybrid training strategy for solving mixed-integer linear programming problems, which may include:

[0059] The modeling of the MILP problem adopts a bipartite graph structure composed of variable nodes and constraint nodes. The bipartite graph structure includes variable node features, constraint node features, and edge features.

[0060] Based on the bipartite graph structure, variables are divided into stable variables whose values ​​are consistent across all solutions and unstable variables whose values ​​differ across all solutions. A corresponding total loss function is then determined, based on the cross-entropy loss function and the contrastive loss function. Furthermore, the stable variables are labeled using the cross-entropy loss function, while the unstable variables are labeled using the contrastive loss function.

[0061] Based on the total loss function, the corresponding prediction model is obtained by training and learning the correspondence between the MILP problem and the solution, and the MILP problem is solved based on the prediction model.

[0062] Furthermore, the solution process for mixed-integer linear programming problems also includes a dynamic search region adjustment mechanism based on the analysis of irreducible infeasible subsystems (IIS) to achieve adaptive adjustment of the search region.

[0063] The structure-aware hybrid training strategy and the dynamic search region adjustment mechanism based on IIS analysis in the above embodiments of the present invention will be described in detail below.

[0064] (a) Hybrid training strategy for structure awareness

[0065] In this embodiment of the invention, the modeling of the MILP problem employs a bipartite graph structure composed of variable nodes and constraint nodes. This graph structure not only clearly expresses the relationship between variables and constraints in the MILP problem but also provides rich input information for the neural network. For example... Figure 1As shown, the hybrid training strategy includes a data processing part and a hybrid training part. The data processing part includes two processes: bipartite graph feature extraction and positive and negative sample collection. The hybrid training part includes two parts: contrastive learning and label learning. The hybrid training part is implemented based on the total loss function, which is based on the contrastive loss function of the contrastive learning part and the cross-entropy loss function of the label learning part. The implementation methods of each part of the hybrid training strategy will be explained in detail below.

[0066] (1) The following will provide a detailed explanation of the eigenvectors contained in each node and edge of the bipartite graph structure:

[0067] (11) Variable node characteristics

[0068] Each variable node is assigned a feature vector consisting of 18 dimensions, which may include:

[0069] Normalized objective function coefficients: used to reflect the importance of variables in the objective function;

[0070] The average of the coefficients of the variable across all constraints: used to provide the overall effect of the variable under all constraints;

[0071] Degree of a variable node: the number of constraints a variable participates in, used to express the degree of its connection in the network;

[0072] The maximum and minimum values ​​of the coefficient of this variable in all constraints: used to describe the upper and lower limits of the variable's scope of influence, respectively;

[0073] Whether it is an integer variable: This is used to distinguish between continuous variables and integer variables using binary identifiers;

[0074] Position encoding: Uses binary encoding, occupying 12 dimensions, to indicate the specific position of the variable in the original sequence.

[0075] (12) Constraint node features

[0076] Each constraint node has a 4-dimensional feature vector, which can specifically include:

[0077] The average value of all coefficients within the constraint: used to reflect the strength of the constraint's influence on the overall set of variables;

[0078] Degree of a constraint node: This indicates the number of variables associated with the constraint, which helps to understand the importance of the constraint in the network.

[0079] The normalized value of the right-hand side term b: It is the constant value of the right-hand side of the constraint after standardization, which facilitates cross-constraint comparison;

[0080] Constraint direction: This refers to the use of symbols such as "≤", "=", "≥" to precisely describe the nature of the constraint condition.

[0081] (13) Edge features

[0082] Each edge represents the form in which a specific variable exists within a certain linear constraint, represented by a single-dimensional eigenvector, i.e.:

[0083] The specific coefficients of a variable in a linear constraint: used to directly reflect the actual magnitude of the variable's contribution to the constraint.

[0084] (2) The process of determining the total loss function will be explained below;

[0085] Due to the complexity of MILP problems, not all correspondences between MILP solutions and their corresponding relationships can be learned. The goal of machine learning is to maximize the effective information learned. Therefore, to learn more information more fully, this invention proposes a novel hybrid training strategy based on the different types of variables in the solution and leveraging the advantages of different learning paradigms.

[0086] Based on the set of historical best or near-best solutions, the consistency of values ​​for each variable in the solutions is analyzed, and the variables are divided into "stable variables" and "unstable variables." Different training strategies are used for training and learning for the two types of variables, which may include:

[0087] (21) Label-supervised learning (i.e., label learning): For stable variables, label-supervised learning is used for training. Specifically, the cross-entropy loss function is used. Optimize the model to guide it to accurately learn the value patterns of stable variables in the optimal solution, thereby reducing the prediction error rate;

[0088] (22) Contrastive learning strategy: For unstable variables, a contrastive learning strategy is used for training. By constructing positive and negative sample pairs, a contrastive loss function is introduced. By narrowing the feature embedding distance of positive sample pairs and widening the feature embedding distance of negative sample pairs, the model is helped to learn effective association patterns between unstable variables, thereby enhancing the model's ability to distinguish between the good and bad combinations of variables.

[0089] Joint design of loss functions: combining the loss functions of the two classes of variables into the total loss function, i.e.:

[0090] ;

[0091] ;

[0092] ;

[0093] in, For cross-entropy loss, To compare the losses, α is the weighting coefficient. It is a variable The true label, The variables predicted by the model The probability of taking the value 1. It's a temperature parameter. It is a positive sample , It is a negative sample .

[0094] The training process of the prediction model may include: training using the Adam (adaptive moment estimation) optimizer, setting the learning rate to 0.001, the batch size to 8, training for 1000 epochs, and adopting an early stopping strategy to prevent overfitting. The early stopping strategy is to dynamically monitor the model performance during the training process and terminate the training early when the preset conditions are met (such as the validation set performance failing to improve for several consecutive epochs) but the training epochs have not been reached.

[0095] (3) The collection process of the positive and negative sample sets used in the corresponding training and learning process will be described in detail below;

[0096] (31) Collection of positive sample set

[0097] For each MILP instance, a predetermined number of optimal or near-optimal solutions are collected using an advanced solver (such as Gurobi) over a predetermined time period t (e.g., 3600 seconds) to form a positive sample set. .

[0098] (32) Collection of negative sample sets

[0099] To support contrastive learning, a negative sample set needs to be generated. The corresponding negative sample set collection process may include:

[0100] First, it is necessary to identify the variable type and analyze the positive sample set. The solution in the solution determines the stable variable (whose value is consistent across all solutions, and whose index is ). ) and unstable variables (values ​​differ across all solutions, indexed as ) and unstable variables (values ​​differ across all solutions, indexed as ) );

[0101] Then, perturb some (partial or all) of the variables to make the solution infeasible, i.e., for each positive sample. For example, you can choose randomly. Take 10% of the variables, flip their values ​​(0 becomes 1, 1 becomes 0), then verify whether the perturbed solution is feasible. If not, add it to the negative sample set. ;

[0102] In the above processing, multiple samples can be generated for each positive sample (let's say...). (1) infeasible solutions It is usually set to 10 to balance the sample size.

[0103] (4) Prediction model architecture based on hybrid training strategy

[0104] In this embodiment of the invention, a graph neural network is used to encode the bipartite graph structure of the MILP problem and establish a prediction model. The established prediction model may specifically include the following modules:

[0105] Input layer: Construct feature vectors for variable nodes and constraint nodes respectively. The input includes the upper and lower bounds of variables in the MILP model, whether they are integer variables, the right-hand side values ​​of constraints, and the sparse structure of the coefficient matrix.

[0106] Graph convolutional layer: It adopts a two-layer graph convolutional structure, first aggregating information from variable nodes to constraint nodes, and then passing it back from constraint nodes to variable nodes in reverse.

[0107] Activation function: ReLU activation function is used after each layer;

[0108] Output layer: The probability of each variable taking the value 1 is output through the MLP as the prediction result.

[0109] During the above processing, the calculation results of each loss function are backpropagated to guide the parameter updates of modules such as graph convolutional layers and MLPs, thereby optimizing the model.

[0110] (ii) Adaptive search region adjustment mechanism

[0111] The trust region search mechanism constructs a trust region based on the predicted probability during the prediction phase. Controlled The neighborhood is constructed by selecting the neighborhood with the highest prediction probability. Each variable is fixed at 1, representing the lowest predicted probability. One variable is fixed at 0; to increase fault tolerance, it is allowed to... The error of each variable, usually It is fixed as a constant.

[0112] An irreducible infeasible subsystem (IIS) is a subset of constraints that are mutually contradictory, for which no solution can satisfy all constraints, and which becomes feasible by removing one constraint.

[0113] In practical applications, if we consider the trust region search mechanism... If the variable is fixed as a constant, due to the complexity of the MILP problem and the characteristics of neural networks, it cannot be guaranteed that the prediction of variables for unseen instances will still be feasible, which greatly limits the application potential of this field. To solve this problem, the adaptive search region adjustment mechanism proposed in this embodiment of the invention refers to... Figure 2 As shown, it includes an inference phase and an adaptive neighborhood adjustment phase, wherein:

[0114] In the inference phase, the MILP model is analyzed using a graph structure, and the probability of each variable taking a specific value is output. The corresponding probability reflects the "confidence" of the variable in the optimal solution - the higher the value, the greater the probability that the variable will take the corresponding value. The output of the inference phase serves as the data basis for the adaptive neighborhood adjustment phase.

[0115] The specific implementation process of the corresponding adaptive neighborhood adjustment phase may include:

[0116] Feasibility verification: Verification using the MILP solver Is a solution within the neighborhood feasible?

[0117] Irreducible Infeasibility Subsystem (IIS) Extraction: This involves infeasibility assessment and IIS calculation, including: based on the feasibility test results, if the solution is infeasible, using a filtering algorithm or other heuristic algorithm to approximate the IIS and identify the subset of constraints that cause the conflict; if the solution is feasible, using a solver to continue solving the remaining variables.

[0118] Local variable unlocking mechanism: namely, variable correction and adjustment of dependency domain, including: locating key variables based on constraints involved in IIS, removing their fixed state, reconstructing the search area, and improving the feasibility of the solution;

[0119] The iterative correction strategy involves the solver solving the problem. If the solution is still not feasible after adjusting the search area, the IIS will continue to be identified and the set of unfixed variables will be gradually expanded according to the set strategy until the predicted solution is a feasible solution.

[0120] The above process mainly includes: first, identifying the core constraint set (IIS) that makes the predicted solution infeasible, and then gradually expanding the range of unfixed variables by adjusting the values ​​of these variables until a feasible predicted solution that satisfies all constraints is obtained.

[0121] The embodiments of the present invention have significant advantages and positive effects compared with the prior art in solving mixed-integer linear programming problems, specifically including:

[0122] High Accuracy: The hybrid training strategy significantly improves the accuracy of variable prediction through targeted training. For stable variables (variables whose values ​​remain consistent across the optimal solution set), cross-entropy loss is used for precise label learning, ensuring the model accurately captures the fixed patterns of these variables. For variable variables (variables whose values ​​change across the solution set, i.e., unstable variables), contrastive loss is used for contrastive learning to uncover complex dependencies between variables and enhance the model's understanding of the problem structure. This dual training mechanism enables the prediction model to exhibit higher accuracy across various MILP instances and reduces the generation of suboptimal solutions.

[0123] Strong Adaptability: The adaptive search region adjustment mechanism overcomes the limitations of fixed hyperparameters in traditional methods by dynamically adjusting the search region. Trust region search mechanisms rely on a fixed number of variables and a fixed range of perturbations, making them difficult to adapt to MILP instances of varying sizes and complexities, often leading to infeasible solutions or excessively large search spaces. In this embodiment, however, the irreducible infeasible subsystem (IIS) analysis is used. When a predicted solution is infeasible, the variables causing the conflict are accurately identified and released from fixed constraints, dynamically updating the search region. This adaptive mechanism ensures the method's robustness to diverse datasets.

[0124] High efficiency: The embodiments of this invention work in both the training and inference phases, enabling the model to outperform existing technologies in both solution speed and solution quality. Ultimately, it achieves a significant efficiency improvement compared to traditional SCIP solvers.

[0125] In summary, the embodiments of the present invention provide a novel MILP optimization method that can balance prediction accuracy and solution feasibility. It can dynamically adjust the search region while maintaining efficient solution, enhance the tolerance to prediction errors, and demonstrate good stability and practicality in various types of practical problems.

[0126] The embodiments of the present invention have also undergone corresponding practical tests in specific scenarios, and the specific application process will be described in detail below.

[0127] Furthermore, taking the factory order scheduling scenario as an example, the implementation process of the technical solution provided by the embodiments of the present invention is described in detail below, demonstrating its specific application in industrial optimization; in this scenario, machine learning methods are used to learn the relationship between similar MILP problems and solutions, thereby accelerating scheduling optimization.

[0128] (1) Application scenario description

[0129] In manufacturing, order scheduling is a crucial step in optimizing production efficiency and resource utilization. Taking a single factory as an example, its basic parameters (such as the number of machines, processing flow, and process dependencies) typically remain constant, while order data (such as order quantity, delivery deadlines, and product types) changes with the production cycle (e.g., monthly). Since scheduling problems across different cycles have similar constraint structures and variable relationships, they can be modeled as similar MILP problems. In the implementation of this invention, machine learning methods are used to learn the mapping relationship between these MILP problems and high-quality solutions, significantly improving the speed and quality of scheduling optimization. This invention uses the monthly order scheduling of a factory as an example to demonstrate how this method can be modeled, trained, and deployed to achieve efficient production scheduling.

[0130] (2) The process of MILP modeling

[0131] The order scheduling problem is modeled as a MILP problem, with the goal of minimizing production costs (including machine operating costs, overtime expenses, and late payment penalties) or maximizing production efficiency (such as on-time delivery rate). Possible modeling methods include:

[0132] Decision variables:

[0133] Binary variables: represent the decision of order allocation to machine (e.g., whether order i is allocated to machine j).

[0134] Continuous variables: represent the start time or completion time of a process.

[0135] Auxiliary variable: indicating whether overtime work or whether the work is postponed.

[0136] Objective function: Minimize total cost, which is a weighted sum of machine uptime, inventory costs, and delay penalties.

[0137] Constraints:

[0138] Resource constraints: Each machine can only process one process at a time, and the number of machines is fixed (e.g., 10 machines).

[0139] Process dependency: The processes in an order must be completed in a specified order (e.g., assembly after processing).

[0140] Time constraints: Orders must be completed before the delivery deadline, taking into account processing time and machine availability.

[0141] Integer constraint: Binary variables take the values ​​0 or 1.

[0142] To generate diverse MILP instances, a baseline problem is built based on real order data (e.g., a month containing 100 orders involving 500 processes), and the data is augmented in the following ways:

[0143] Randomly adjust order volume (80-120 orders), delivery time (±10% fluctuation), or product type ratio. Introduce additional constraints (such as periodic maintenance downtime for some machines) or relax constraints (such as allowing limited delays). Using the above methods, generate 400 MILP instances, each containing variables and constraints that fluctuate within a certain range to ensure coverage of the diversity of scheduling issues.

[0144] (3) Data generation

[0145] To train a machine learning model, a dataset of positive and negative samples needs to be generated:

[0146] Positive sample collection:

[0147] For each MILP instance, run the Gurobi solver for 3600 seconds (or longer) to collect 50 high-quality solutions that are as close as possible to the best known solution. Each solution is stored as a vector containing the assignments (0 or 1) of binary variables and the values ​​of continuous variables (such as process start time).

[0148] Negative sample generation:

[0149] Analyze the positive sample set to identify stable variables (values ​​consistent across all solutions, such as fixed machine allocation) and unstable variables (values ​​that change, such as process time). For each positive sample, randomly select 10% (at least one) of the unstable variables and flip their binary values ​​(0 becomes 1, 1 becomes 0). Verify whether the perturbed solution is infeasible (violation of constraints, such as machine overload or process conflict). Generate 10 negative samples for each positive sample.

[0150] (4) Model training

[0151] The specific steps for training a graphical convolutional network model using a hybrid training strategy are as follows:

[0152] Data partitioning:

[0153] Training set: 240 MILP instances.

[0154] Validation set: 60 instances.

[0155] Test set: 100 instances.

[0156] Graph Neural Network Architecture: Using the model architecture described in this embodiment of the invention, the input is a bipartite graph feature, and the output is the probability that each variable takes the value 1.

[0157] Training process: The specific process is as described above and will not be repeated here.

[0158] (5) Model Deployment

[0159] The trained model is deployed into the factory's production scheduling optimization program. The specific process is as follows:

[0160] System Integration:

[0161] Embed the model into an Enterprise Resource Planning (ERP) system or dedicated scheduling software, receiving new order data via an API interface. The software automatically converts the order data (order quantity, deadline, product type) into a bipartite graph format, including constraint nodes, variable nodes, and edge features. The model processes the bipartite graph and outputs the probability of each variable taking the value 1.

[0162] Optimize the process:

[0163] Based on the trust region search mechanism, variables are arranged in descending order of probability, with high-probability variables fixed at 1 and low-probability variables at 0, defining a trust region (i.e., a delta neighborhood). A solver is used to solve the MILP problem constrained by the trust region. If the solution is infeasible, an adaptive search region adjustment mechanism is triggered, calculating the irreducible infeasible subsystem (IIS), identifying conflicting variables, dynamically releasing fixed constraints, and resolving the problem, iterating a maximum of 5 times to ensure a balance between solution efficiency and feasibility. The final solution (including machine allocation and process time) is fed back to the production system to generate a scheduling plan.

[0164] The above-mentioned application in the solution process of mixed integer linear programming problems can balance prediction accuracy and solution feasibility, and demonstrates good stability and practicality in various types of practical problems.

[0165] Through the above description of the embodiments, those skilled in the art can clearly understand that the above embodiments can be implemented by software, or by using software plus necessary general-purpose hardware platforms. Based on this understanding, the technical solutions of the above embodiments can be embodied in the form of a software product, which can be stored in a non-volatile storage medium (such as a CD-ROM, USB flash drive, mobile hard drive, etc.), including several instructions to cause a computer device (such as a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments of the present invention.

[0166] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims. The information disclosed in the background section is intended only to enhance the understanding of the overall background technology of the present invention and should not be construed as an admission or implication in any way that such information constitutes prior art known to those skilled in the art.

Claims

1. An optimization method for solving mixed-integer linear programming problems, characterized in that, include: The modeling of the mixed integer linear programming (MILP) problem adopts a bipartite graph structure composed of variable nodes and constraint nodes. The bipartite graph structure includes variable node features, constraint node features, and edge features. Based on the bipartite graph structure, variables are divided into stable variables with consistent values ​​across all solutions and unstable variables with discrepancies, and a corresponding total loss function is determined. The total loss function is based on the cross-entropy loss function and the contrastive loss function. Furthermore, the stable variables are labeled using the cross-entropy loss function, and the unstable variables are labeled using the contrastive loss function. Based on the total loss function, the corresponding prediction model is obtained by training and learning the correspondence between the MILP problem and the solution, and the MILP problem is solved based on the prediction model.

2. The method according to claim 1, characterized in that, The variable node features include: Normalized objective function coefficients are used to reflect the importance of variables in the objective function; The average value of the variable's coefficients across all constraints is used to provide the overall impact of the variable under all constraints. The degree of a variable node is used to express how closely connected it is in the network; The maximum and minimum values ​​of the coefficient of this variable in all constraints are used to describe the upper and lower limits of the variable's scope of influence, respectively. Whether it is an integer variable is used to distinguish between continuous variables and integer variables using binary identifiers; Position encoding is used to indicate the specific position of a variable in the original sequence.

3. The method according to claim 1, characterized in that, The characteristics of the constraint node include: The average value of all variable coefficients within the constraint is used to represent the strength of the constraint's influence on the overall variable set; The degree of a constraint node indicates the number of variables associated with that constraint; The normalized value of the right-hand side term is the constant value of the constraint right-hand side after standardization. Constraint direction, used to describe the properties of the constraint condition.

4. The method according to claim 1, characterized in that, The edge features include: The specific coefficients of a variable in a linear constraint are used to directly reflect the actual magnitude of the variable's contribution to that constraint.

5. The method according to any one of claims 1 to 4, characterized in that, The calculation method for the total loss function includes: ； ； ； in, Let cross-entropy be the loss function. To compare the loss functions, α is the weighting coefficient. It is a variable The true label, The variables predicted by the model The probability of taking the value 1. It's a temperature parameter. It is a positive sample , It is a positive sample set. It is a negative sample , It is a negative sample set.

6. The method according to claim 5, characterized in that, The collection process of the positive sample set includes: For each MILP instance, a predetermined number of optimal solutions are collected using an advanced solver over a predetermined time period to form a positive sample set.

7. The method according to claim 5, characterized in that, The collection process of the negative sample set includes: Analyze the solutions in the positive sample set, determine the stable and unstable variables, and randomly perturb a predetermined number of variables in the unstable variables to make the solution infeasible. When a solution is infeasible, it is added to the negative sample set. Each positive sample is used to generate multiple infeasible solutions.

8. The method according to claim 7, characterized in that, The random perturbation processing includes: randomly selecting 10% of the unstable variables, flipping the values ​​of the corresponding unstable variables, i.e., changing 0 to 1 and 1 to 0; and each positive sample is used to generate 10 infeasible solutions.

9. The method according to claim 5, characterized in that, The training process of the prediction model includes: The Adam optimizer was trained using adaptive moment estimation with a learning rate of 0.001, a batch size of 8, and 1000 epochs. An early stopping strategy was employed to prevent overfitting. This strategy involves dynamically monitoring the model performance during training and terminating training prematurely when preset conditions are met but the required number of epochs has not been reached.

10. The method according to claim 5, characterized in that, The method also includes an adaptive search region adjustment process, which includes: The MILP solver is used to verify whether the solution within the neighborhood is feasible. If the solution is not feasible, a heuristic algorithm is applied to extract the irreducible and infeasible subsystem (IIS) and identify the subset of constraints that cause the conflict. If the solution is feasible, the solver is used to continue solving the remaining variables. Based on the constraints extracted from IIS, identify the key variables, remove their fixed state, and reconstruct and adjust the search area. If adjusting the search area still does not work, the extracted IIS will continue to be identified and the set of unfixed variables will be gradually expanded according to the set strategy until the predicted solution is a feasible solution.

Images