Machine Learning Insides OptVerse AI Solver:
Design Principles and Applications

Data, i.e., the collection of MILP and SAT instances plays a fundamental role in developing advanced MILP and SAT solvers, primarily from three aspects. First, for tasks such as hyperparameter configuration and ML model training, a substantial and diverse set of realistic and independent identically distributed (i.i.d.) instances is necessary. Second, the evaluation of solvers requires as many instances as possible to identify potential issues and weaknesses through white-box testing. Finally, to inspire more efficient algorithms and promote research investigations, the research community keeps calling for challenging instances for better benchmarking and competitions.

However, the limited availability of real-world instances, due to labor-intensive data collection and proprietary issues, remains a critical challenge and often leads to sub-optimal decisions and biased assessments. This challenge motivate a suite of synthetic instance generation techniques, which fall into two categories. Some methods rely heavily on expert-designed formulations for specific problems, such as Traveling Salesman Problems (TSPs) or Set Covering problems. However, these methods cannot cover real-world applications where domain-specific expertise or access to the combinatorial structures is limited. Other methods construct general MILP and SAT instances by sampling from an encoding space that controls a few specific statistics. However, these methods often struggle to capture the rich features and the underlying combinatorial structures, resulting in an unsatisfactory alignment with real-world instances. Moreover, these works have scarcely investigated the ability of the generators to produce challenging MILP and SAT instances.

We argue that developing deep learning (DL)-based generators is a promising way to address the challenge of limited data availability for the following reasons. First, these generators can actively learn from real-world instances and generate new ones without relying on expert-designed formulations. The generated instances can simulate realistic scenarios, cover a wide range of cases, significantly enrich datasets, and facilitate the development of solvers at a relatively low cost. Second, learning-based generators have the capacity to capture instance features, enabling them to efficiently explore the combinatorial problem space to construct challenging instances. Finally, this approach shows promising technical prospects for understanding the problem space and learning representations. Although deep learning techniques have been widely applied to generate Boolean satisfiability (SAT) problems, but the previous work all ignore retaining the difficulty of generated SAT instances. Besides, the development of deep learning-based MILP instance generators remains a complete blank due to higher technical difficulties, i.e., it involves not only the intrinsic combinatorial structure preservation but also high-precision numerical prediction.

To enhance the data reliability of SAT instances, we conduct an in-depth analysis of the computational hardness degradation issue in SAT instance generation and propose HardSATGEN, which to our best knowledge, is the first generative model that generates instances maintaining similar computational hardness and capturing the global structural properties simultaneously. We adopt the one-to-one bipartite graph representation for SAT instances, i.e. the literal-clause graph (LCG), which is composed of a literal node set and a clause node set, and an edge therein indicates the occurrence of the literal in the clause. The methodology lies upon the node split-merge framework You et al. (2019), which first gradually splits clause nodes in the bipartite graph to form a forest of template and then enforces a graph neural network (GNN) to learn the reverse node merging process to generate new bipartite graphs. A concrete analysis of the hardness degradation issue for the one-stage node split-merge framework reveals that the root cause lies in the inherent homogeneity in the split-merge procedure and the difficulty of semantic formation of structures led by oversplit substructures. In addition, we observe that industrial formulae exhibit community structure correlated with hardness. Based on the analysis, we introduce a fine-grained control mechanism over the community structure and the unsatisfiable cores (of unsatisfiable instances) to the split-merge paradigm and propose a multi-stage generation pipeline involving scrambling the core, connecting within communities, and connecting cross the communities to guide the iterative split-merge procedure for better similarity in structure and computational hardness.

Specifically, as shown in Fig. 3, for an input formula, core detection on the formula and community detection on its variable-incidence graph (VIG) representation, which consists of variables as nodes and the edges indicating the co-occurrence of two variables in one clause, are applied to incorporate knowledge about community and core structures into the graphs. During the splitting phase, we retain the detected core and then remove cross-community connections and the remaining in-community connections to obtain split node pairs, which serve as the training data for two constructive phases, and the remained graph is saved as the templates. In training, two GNNs are trained to correspond to cross-community and in-community connections, learning the structural features and predicting whether the specific node pair should be merged, which is framed as a binary classification task. The trained GNNs are utilized in the inference phase to generate new instances, where the core scrambling operation involves random permutation of variables and clauses, as well as random flipping of literals, and the subsequent in-community and cross-community node merging phases utilize the learned models with set operations on the node sets to construct connections within and cross communities, respectively. For more details about the methodology, please refer to Li et al. (2023c).

To enhance the data reliability of MILP instances, we propose G2MILP, which to the best of our knowledge, is the first deep generative framework for MILP instances. We represent MILP instances as weighted bipartite graphs, where variables and constraints are vertices, and non-zero coefficients are edges. This graph representation enables us to use graph neural networks (GNNs) to effectively capture the essential features of MILP instances using graph neural networks (GNNs). In this way, we recast the original task as a graph generation problem. To accommodate various application scenarios, we consider two task settings for utilizing G2MILP in MILP instance generation: realistic MILP instance generation and hard MILP instance generation, as shown in Figure 4. We begin by focusing on realistic MILP instance generation, where our objective is to generate new MILP instances that closely resemble real-world instances in terms of their structures and computational hardness. Since generating the complex bipartite graphs from scratch can be computationally expensive and may destroy the intrinsic combinatorial structures of the problems, we propose a masked variational autoencoder (VAE) paradigm inspired by the masked autoencoder (MAE) and the variational autoencoder (VAE) theories. The model overview is in Figure 5. Specifically, this paradigm iteratively corrupts and replaces parts of the original graphs using sampled latent vectors. To implement the complicated generation steps, we design a decoder consisting of four modules that work cooperatively to determine multiple components of new instances, encompassing both structure and and numerical prediction tasks simultaneously. Subsequently, we work on hard MILP instance generation, where our goal is to construct challenging MILP instances. To achieve this, we propose a hardness-oriented iterative augmenting scheme. In each iteration, G2MILP generates a batch of new instances and stores the most difficult instances in a storage. We then fine-tune G2MILP using the instances in the storage as a training set, so that the model is specifically oriented towards generating challenging instances. For more details about this technique, please refer to Geng et al. (2023b); Wang et al. (2023b).

Generalizing learning-based solvers to out-of-distribution (OOD) instances poses significant challenges for real-world applications. First, the solvers are often trained on a given dataset and need to process instances of varying sizes in real-world scenarios. The sizes of instances encountered in practice can be much larger than those in the training dataset. Second, the environment changes often introduce perturbations to instance structures. For example, in staff scheduling problems, emergencies may arise and impose additional restrictions that are not included in the previous MILP problems. We thus need to consider additional kinds of constraints representing these restrictions, leading to changes in the problem structures. Third, acquiring or generating more instances from other distributions or with larger sizes may not be feasible due to information security/privacy concerns or high data acquisition costs. For example, in a company’s job-shop scheduling problems, confidential business information, such as the production capacity and costs, can be inferred from the MILP instances. Thus, the company may provide only a few anonymous training instances to avoid privacy disclosure.

The solving time of the learning-based solvers under the aforementioned settings could be much longer compared to traditional heuristics, particularly when the training instances are few. Existing work on the generalization ability of CO solvers mainly focuses on the problem-specific approximation solvers (such as the neural solvers for routing problems) instead of the general exact MILP solvers. For the problem-specific solvers, we can get access to the problem type and then employ domain-randomization type approaches by generating a large number of instances from different distributions Wang et al. (2023a); Geisler et al. (2022). However, data generation has severe limitations when it comes to anonymous datasets where prior knowledge of the problem type is unavailable. Consequently, the problem-specific generation algorithms cannot be applied in such cases. Therefore, we thus focuses on the restrictive practical setting where training data only comes from limited training domains, and acquiring or generating new instances within a short time is intractable. This challenge significantly hinders the applicability of learning-based solvers to real-world scenarios.

The instance augmentation policy, denoted as $\Phi:\mathcal{G}\to\mathcal{G}$, maps an instance bipartite graph $G$ to the augmented instance graph $\Phi(G)$ by applying augmentation operators. However, two challenges for the instance augmentation remain to address. First, employing structural augmentation on a bipartite graph leads to a non-differentiable learning objective for the augmentation policy, making it intractable to train through backward propagation and gradient-based updates. Second, the augmentation policy receives significantly less training data (instances graphs) compared to that of the branching policy (branching samples), resulting in low training efficiency for the augmentation policy. The adversarial augmentation policy operates on MILP instance graphs $\{G_{k}\}_{k}$ to generate augmented instances $\{\Phi(G_{k})\}_{k}$, and each augmented instance can generate a sequence of branching samples.

Addressing the issues of non-differentiability and unsatisfactory training efficiency requires a sample-efficient training algorithm under sparse supervised signals. Thus, we formulate the instances augmentation process as a contextual bandit problem and adopt a contextual-bandit version of the proximal policy optimization (PPO) algorithm to train the augmentation policy Schulman et al. (2017).

To formulate the graph augmentation as a contextual bandit problem, we view the solver as the environment and the adversarial augmentation policy as the agent. The contextual bandit problem can be represented as the tuple $(\mathcal{G},\mathcal{A}^{\text{CB}},r^{\text{CB}})$. To distinguish the notations from the branching MDP, here we use the superscript CB to specify the contextual bandit formulation. Then we specify the state space ${\mathcal{G}}$, the action space $\mathcal{A}^{\text{CB}}$ and reward function $r^{\text{CB}}:{\mathcal{G}}\times\mathcal{A}^{\text{CB}}\to\mathbb{R}$ for a MILP instance as follows.

• •

  The state space $\mathcal{G}$, or context vector space, also known as the context vector space, is the set of instance graph representations of MILP instances.

• •

  An action $a^{\text{CB}}=(V^{\text{CB}},E^{\text{CB}},W^{\text{CB}})\in\mathcal{A}^{\text{CB}}$, consists of subsets of variables $V^{\text{CB}}$, edges $E^{\text{CB}}$, and constraints $W^{\text{CB}}$ to be masked. The agent applies action $a^{\text{CB}}$ to an instance graph $G$ to obtain an augmented instance $\Phi(G)$, on which branching samples are collected.

• •

  The reward The reward function $r^{\text{CB}}(G,a^{\text{CB}})$ is defined as the average loss of the branching policy $F(\pi_{\theta},\Phi(G),\mathcal{T})$ using the samples collected from the augmented graph $\Phi(G)$ derived from $G$.

We describe the training process for the augmentation policy network $\Phi_{\eta}$ with parameters $\eta$ under the contextual bandit formulation as follows. First, the augmentation policy network $\Phi_{\eta}$ takes as input an instance graph $G$ and outputs the masking probability of each node and edge. An action $a^{\text{CB}}$ is constructed by selecting the nodes and edges that are most likely to be masked based on the predicted probability and a predefined masking proportion. We denote the predicted probability of the selected action $a^{\text{CB}}$ by $P_{\eta}(a^{\text{CB}}|G)$. Second, we leverage a learnable state-value function $V_{\alpha}:\mathcal{G}\to\mathbb{R}$ with parameters $\alpha$ to calculate the variance-reduced estimation of the advantage function $\hat{A}=r^{\text{CB}}-V_{\alpha}(G)$. Third, to improve the training efficiency of the augmentation policy, we use a replay buffer $\mathcal{B}$ consisting of the tuples $(G,\hat{A},P_{\eta}(a^{\text{CB}}|G)$ obtained from the history augmentation policy for experience replay. We use $\lambda(\eta)=\frac{P_{\eta_{\text{current}}}(a^{\text{CB}}|G)}{P_{\eta_{\text{old}}}(a^{\text{CB}}|G)}$ to represent the likelihood ratio between the current policy and old policies sampled from the buffer. Finally, the training objective for the adversarial augmentation policy, which includes the policy network and state-value function, is achieved by minimizing the following equations:

	$\displaystyle L_{V}(\alpha)$	$\displaystyle=\mathbb{E}_{\mathcal{B}}\left[\hat{A}^{2}\right]\quad$		(1)
	$\displaystyle L_{\Phi}(\eta)$	$\displaystyle=\mathbb{E}_{\mathcal{B}}\left[\min(\lambda(\eta),\text{clip}(\lambda(\eta),1-\epsilon,1+\epsilon))\hat{A})\right],$		(2)

where $\epsilon$ is the clip ratio.

For more details about this technique, please refer to Liu et al. (2023).

To evaluate the performance in real-world applications, we deploy AdaSolver to a real-world Anonymous dataset, a widely-used subset of the MILPLIB library. The number of constraints in the Anonymous dataset ranges from 3,375 to 159,823, and the number of variables ranges from 1,613 to 92,261. The wide range of problem sizes as well as real-world perturbations pose significant challenges for the learning-based methods. The results in Table 1 demonstrate that AdaSolver is able to improve the performance of both IL- and RL-based solvers. AdaSolver-IL achieves the best PD gap while AdaSolver-RL achieves the best PD integral metric. The performance of AdaSolver sheds light on the generalizable learning-based solvers in real-world applications.

The Simplex method Dantzig and Wolfe (1960) is a pioneering method for solving large-scale LP problems. It starts with initial basis $\mathcal{B}^{(0)}$ and routinely pivots to a neighboring basis with improvement till reaching an optimal basis $\mathcal{B}^{*}$. The performance of the Simplex method largely depends on the choice of initial basis. Current strategies Bixby (1992), Gould and Reid (1989), Junior and Lins (2005), Ploskas et al. (2021), Galabova and Hall (2020) for selecting an advanced initial basis are mostly rule-based, requiring extensive expert knowledge and empirical study. Yet, most of them still fail to exhibit consistent improvement over the canonical initial basis, because they often rely on heuristics that do not generalize well. Nevertheless, the quest for universally effective strategies may not be necessary in practice, as we often encounter sets of LP problems that exhibit considerable similarities. For instance, in scenarios like a manufacturer’s daily production planning or an airport’s hourly flight scheduling, similar LP problems regularly occur. Therefore, it is natural to ask if a learning-based approach could be more effective, leveraging the correlation among LP problems and enhancing the efficiency of the Simplex method.

Denote the set of past solved Linear Programs (LPs) as $\mathcal{D}=\left\{[(P^{k}),(x^{k},s^{k})]\right\}_{k=1}^{K}$, where $(P^{k})$ is the $k$-th problem and $(x^{k},s^{k})$ is its optimal solution. Given $\mathcal{D}$ and a new LP $P^{test}$ from the same distribution, the goal is to predict an initial basis for the new LP. This basis should enable the Simplex algorithm to converge more rapidly, taking advantage of the inherent similarities within the set $\mathcal{D}$. We consider the following standard format of LP:

$$\begin{array}[]{cl}\min_{x\in\mathbb{R}^{n},s\in\mathbb{R}^{m}}&c^{T}x\\ \text{s.t.}&Ax=s\\ &\ell^{x}\leq x\leq u^{x}\\ &\ell^{s}\leq s\leq u^{s},\end{array}$$

(P)

where $\ell^{x},\ell^{s}$ can reach $-\infty$, and $u^{x},u^{s}$ can reach $\infty$.

This problem presents several challenges: (1) LP problem $P^{k}$ is of varying size. (2) For each LP problem, the number of potential bases is exponential w.r.t the problem size. (3) The selected basis has to be valid, i.e., the corresponding basis matrix is non-singular and the states of non-basic variables are consistent with their bounds.

We represent a linear programming (LP) problem as a bipartite graph, as suggested by Gasse et al. (2019a). This approach allows for the effective handling of LP problems of varying sizes and complexities. Furthermore, we employ a 2-layer Graph Convolutional Network (GCN), following the model proposed by Morris et al. (2019), to predict the optimal basis for an LP problem. The bipartite graph representation effectively captures the relationship between variables and constraints in an LP problem, transforming it into a structure that can be efficiently processed by the GCN. The 2-layer GCN then processes this graph to predict the status of each variable in the LP problem, classifying them into categories such as basic, non-basic at the upper bound, or non-basic at the lower bound. This prediction is essential for selecting an appropriate initial basis for the simplex method. In this way, we transform the complex task of initial basis selection into a more manageable classification task.

To ensure the predicted status of each variable is consistent with its bounds in the LP formulation, a knowledge-based masking technique Fan et al. (2022) is employed. This technique masks out probability entries based on the problem’s knowledge, ensuring that the resulting probabilities satisfy the feasibility of non-basic entries. During the inference phase, after generating a candidate basis from the predicted probabilities, the basis needs to be adjusted to ensure it is valid, i.e., the matrix $[A_{\mathcal{B}_{x}}\enspace-I^{m}_{\mathcal{B}_{s}}]$ is non-singular. This adjustment is inspired by a basis repair procedure, which involves factorization steps and, if necessary, removing columns that do not meet certain criteria to achieve a complete and successfully factored basis. Figure 6 shows the inference steps of the proposed method. Please refer to Fan et al. (2023b) for more details.

It is crucial that the LP problems solved using this method are similar in nature. This similarity ensures that the learned model can effectively generalize and apply its learned knowledge to new but related LP problems. Our method works best when applied to LP instances that share characteristics with those in the training set. In practice, the deployment steps could be:

1. 1.

   Initial Phase: In the initial phase, users upload LP problems which are solved completely using the default strategy to accumulate past solving experience. As LP problems are solved, their data (including the problems themselves and their solutions) are stored. This phase is critical for accumulating a diverse and representative training set, which is essential for training an effective GNN model.

2. 2.

   Training the GNN Model: When enough correlated data is accumulated, there is an option to start training the GNN model. Actually, it remains difficult to decide when to start training. In practice, we can divide the data into training and validation sets. The performance of trained GNN is closely monitored on the validation set. If the validation performance indicates a significant acceleration of the simplex algorithm, it implies the time to start training.

3. 3.

   Applying the Method to Future LPs: For future LP problems provided by the user, the trained GNN model is employed. The model accelerates the solving process by predicting an initial basis that is closer to the optimal solution, thereby reducing the number of iterations and overall computation time required for solving the LP.

Through these steps, the proposed method can be deployed to practical scenarios. Its effectiveness is evidenced by the speed enhancements observed in various benchmarks. Table 2 highlights the end-to-end time reduction achieved using our GNN model across different datasets: MIRP Papageorgiou et al. (2014), LIBSVM Zhu et al. (2003); Applegate et al. (2021), StochasticSupplyChain Castro and de la Lama-Zubirán (2020), a synthetic dataset (Generated) Bowly et al. (2020), and two Huawei supply chain datasets – SupplyChain (Small) and SupplyChain (Medium). The effectiveness of the GNN model is contrasted with baseline methods CA Bixby (1992), CA-MPC Ploskas et al. (2021), and CA-ANG Junior and Lins (2005), demonstrating its superior performance in accelerating LP problem solving by harnessing past-solving experiences and the GNN’s predictive capabilities.

Presolve simplifies input LP problems by equivalent transformations. It is widely recognized as one of the most crucial components in modern linear programming (LP) solvers Maros (2002); Elble (2010). In practice, we found LP problems from industry usually contain much redundancy, which could come from reasons like general purpose modeling systems and non-expert modelings, and it can severely decrease the efficiency and reliability of LP solvers to solve LPs Maros (2002). There are many types of redundancy, including multiple constraints (rows) that are linear dependent with each other, a single variable (column) whose value can be fixed in advance Mészáros and Suhl (2003), etc. Modern LP solvers integrate a rich set of presolvers to handle the redundancy from different aspects Andersen and Andersen (1995). Within OptVerse AI solver, we employ eighteen different presolvers to handle various redundancy from different aspects. An illustrate of the presolve process with a toy example can be found in Figure 7.

In practice, we found that presolve routines, i.e., a program determining a sequence of presolvers successively executed in the presolving phase, play critical roles in the efficiency of solving LPs. Previous research on designing high performance presolve routines is relatively limited. Thus, we empirically conclude three points in presolve routine design. That is, which presolvers to select, in what order to execute, and when to stop. Designing such presolve routines is challenging due to the enormous search space and the extensive requirement for expert knowledge. Moreover, hard-coded routines integrated in many modern solvers Forrest et al. (2022); Galabova (2023) cannot capture distinct patterns of different problems to achieve high performance.

In the OptVerse AI solver, we propose the first learning-based approach—that is, reinforcement learning for presolve (RL4Presolve)—to accelerate presolve in large-scale LPs. Specifically, there are three core components in RL4Presolve (see Figure 8):

1. 1.

   The Markov decision process (MDP) formulation. We formulate the task as a Markov decision process by defining : a) the presolve state as a 51-dimensional feature vector based on the activation condition of each presolver. b) the action space as the set of all possible presolver sequences. That is, $\mathcal{A}=\left\{\boldsymbol{a}\mid\boldsymbol{a}=(a_{i})_{i=1}^{n},\,a_{i}\in\mathcal{A}_{0},n\in\mathbb{N}\right\}$, where $\mathcal{A}_{0}$ is the set of all available presolvers in OptVerse. c) the environment transition as the problem changes after the presolver sequence is executed. d) the reward function as the executing time between two steps (i.e., $r=-t$). e) the discount factor $\gamma=1$, and thus the absolute value of the cumulative rewards equals to the total solving time of the input problem.

2. 2.

   The agent with adaptive action sequence. Compared to traditional reinforcement learning (RL) algorithms, we need to explicitly consider the cumulative decision time. Thus, select a single presolver for each decision can be relatively expensive. To tackle this problem, we propose a novel approach that replaces primitive presolvers with automatically generated presolver sequences at each step. Specifically, we employ an agent that learns a one-step conditional probability of executing presolver $j$ after presolver $i$ is executed at last step, which is motivated from combos in video games. Then, we can sample a sequence of presolvers from this learned probability matrix at each step. Intuitively, this is an automatic way to generate high-quality “combos” as each step.

3. 3.

   Training algorithm with proximal policy optimization (PPO). Since evaluating the learned presolve routines is time-consuming, we employ the training efficient on-policy RL algorithm PPO to optimize the one-step conditional probabilities. Specifically, the probability ratio can be written as:

   $$r(\boldsymbol{\theta})=\frac{\pi_{\boldsymbol{\theta}}\left(a_{1}\mid\mathbf{s}\right)\prod_{i=2}^{n+1}\pi_{\boldsymbol{\theta}}\left(a_{i}\mid a_{i-1},\mathbf{s}\right)}{\pi_{\boldsymbol{\theta}_{\text{old}}}\left(a_{1}\mid\mathbf{s}\right)\prod_{i=2}^{n+1}\pi_{\boldsymbol{\theta}_{\text{old}}}\left(a_{i}\mid a_{i-1},\mathbf{s}\right)}=\prod_{i=1}^{n+1}\frac{\pi_{i}(\boldsymbol{\theta})}{\pi_{i}(\boldsymbol{\theta}_{\text{old}})},$$

   (3)

   where $\pi_{1}(\boldsymbol{\theta})=\pi_{\boldsymbol{\theta}}(a_{1}\mid\boldsymbol{s})$, $\pi_{i}(\boldsymbol{\theta})=\pi_{\boldsymbol{\theta}}(a_{i}\mid a_{i-1},\boldsymbol{s})$ for $i\geq 2$, and $a_{i}$ is the $i$th presolver in OptVerse.

Please refer to Kuang et al. (2023b) for more details.

Applying the RL agent directly to real-world applications is usually challenging due to the hardware constraints for high-end GPUs. Thus, we extract rules from learned policies for simple and efficient deployment to the solver. Specifically, we observe that RL4Presolve tends to generate highly similar action probabilities for instances from similar distributions, which might because that independent and identically distributed instances from real-world tasks are usually generated from similar models with different daily data instantiation. Thus, we can extract some new routines from learned policies that outperforms the default one. Specifically, we sample $20$ presolve routines on each task from learned policies and then replace the default routine in OptVerse with the best-performed one. We find this approach consistently improves the hard-coded presolve routines. We report the improved efficiency of solving large-scale LPs on four real-world benchmarks (three of which from Huawei’s supply chain) in Table 3.

Given a standard form of Mixed Integer Linear Programming (MILP), we drop all its integer constraints to obtain its linear programming (LP) relaxation, which takes the form of

$\displaystyle z_{\text{LP}}^{*}\triangleq\min_{\textbf{x}}\{\textbf{c}^{\top}\textbf{x}|\textbf{A}\textbf{x}\leq\textbf{b},\textbf{x}\in\mathbb{R}^{n}\}.$

(4)

Since the problem in (4) expands the feasible set, we have $z_{\text{LP}}^{*}\leq z^{*}$, where $z^{*}$ is the optimal solution to the given MILP. We denote any lower bound found via an LP relaxation by a dual bound. Given the LP relaxation in (4), cutting planes (cuts) are linear inequalities that are added to the LP relaxation in the attempt to tighten it without removing any integer feasible solutions of given MILP. Cuts generated by MILP solvers are added in successive rounds. Specifically, each round $k$ involves (1) solving the current LP relaxation, (2) generating a pool of candidate cuts $\mathcal{C}^{k}$, (3) selecting a subset $\mathcal{S}^{k}\subseteq\mathcal{C}^{k}$, (4) adding $\mathcal{S}^{k}$ to the current LP relaxation to obtain the next LP relaxation, (5) and proceeding to the next round. Adding all the generated cuts to the LP relaxation would maximally strengthen the LP relaxation and improve the lower bound at each round. However, adding too many cuts could lead to large models, which can increase the computational burden and present numerical instabilities (Wesselmann and Stuhl, 2012). Therefore, cut selection is proposed to select a proper subset of the candidate cuts, which is significant for improving the efficiency of solving MILPs (Tang et al., 2020).

Cut selection heavily depends on (P1) which cuts should be preferred, and (P2) how many cuts should be selected (Achterberg, 2007; Dey and Molinaro, 2018). Moreover, we observe from extensive empirical results that (P3) what order of selected cuts should be preferred significantly impacts the efficiency of solving MILPs as well. To improve the efficiency of MILP solvers, recent methods (Tang et al., 2020; Paulus et al., 2022; Huang et al., 2022) propose to learn cut selection policies via machine learning, especially reinforcement learning. They offer promising approaches to learn more effective heuristics by capturing underlying patterns among MILPs from specific applications. However, many existing learning-based methods (Tang et al., 2020; Paulus et al., 2022; Huang et al., 2022)—which learn a scoring function to measure cut quality and select a fixed ratio/number of cuts with high scores—suffer from two limitations. First, they learn which cuts should be preferred by learning a scoring function, neglecting the importance of learning the number and order of selected cuts (Dey and Molinaro, 2018). Second, they do not take into account the interaction among cuts when learning which cuts should be preferred, as they score each cut independently. As a result, they struggle to select cuts that complement each other nicely, which could severely hinder the efficiency of solving MILPs (Dey and Molinaro, 2018). Indeed, we empirically show that they tend to select many similar cuts with high scores.

In the OptVerse AI solver, we propose a novel hierarchical sequence model (HEM) to learn cut selection policies via reinforcement learning (Wang et al., 2023c), which is the first learning-based method that can tackle (P1)-(P3) simultaneously. Specifically, there are three core components in HEM (see Figure 9):

1. 1.

   Reinforcement Learning Formulation As shown in Figure 9, we formulate a MILP solver as the environment and our proposed HEM as the agent. We consider an MDP defined by the tuple $(\mathcal{S},\mathcal{A},r,f)$. Specifically, we specify the state space $\mathcal{S}$, the action space $\mathcal{A}$, the reward function $r:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$, the transition function $f$, and the terminal state in the following. (1) The state space $\mathcal{S}$. Since the current LP relaxation and the generated cuts contain the core information for cut selection, we define a state $s$ by $(M_{\text{LP}},\mathcal{C},\textbf{x}_{\text{LP}}^{*})$. Here $M_{\text{LP}}$ denotes the mathematical model of the current LP relaxation, $\mathcal{C}$ denotes the set of the candidate cuts, and $\textbf{x}_{\text{LP}}^{*}$ denotes the optimal solution of the LP relaxation. To encode the state information, we design thirteen features for each candidate cut based on the information of $(M_{\text{LP}},\mathcal{C},\textbf{x}_{\text{LP}}^{*})$. That is, we actually represent a state $s$ by a sequence of thirteen-dimensional feature vectors. (2) The action space $\mathcal{A}$. To take into account the ratio and order of selected cuts, we define the action space by all the ordered subsets of the candidate cuts $\mathcal{C}$. It can be challenging to explore the action space efficiently, as the cardinality of the action space can be extremely large due to its combinatorial structure. (3) The reward function $r$. To evaluate the impact of the added cuts on solving MILPs, we design the reward function by (i) measures collected at the end of solving LP relaxations such as the dual bound improvement, (ii) or end-of-run statistics, such as the solving time and the primal-dual gap integral. For the first, the reward $r(s,a)$ can be defined as the negative dual bound improvement at each step. For the second, the reward $r(s,a)$ can be defined as zero except for the last step $(s_{T},a_{T})$ in a trajectory, i.e., $r(s_{T},a_{T})$ is defined by the negative solving time or the negative primal-dual gap integral. (4) The transition function $f$. The transition function maps the current state $s$ and the action $a$ to the next state $s^{\prime}$, where $s^{\prime}$ represents the next LP relaxation generated by adding the selected cuts at the current LP relaxation. (5) The terminal state. There is no standard and unified criterion to determine when to terminate the cut separation procedure (Paulus et al., 2022). Suppose we set the cut separation rounds as $T$, then the solver environment terminates the cut separation after $T$ rounds.

2. 2.

   Hierarchical Sequence Model The policy network architecture of HEM is also illustrated in Figure 9. First, the higher-level policy learns the number of cuts that should be selected by predicting a proper ratio. Suppose the length of the state is $N$ and the predicted ratio is $k$, then the predicted number of cuts that should be selected is $\lfloor N*k\rfloor$, where $\lfloor\cdot\rfloor$ denotes the floor function. We define the higher-level policy by $\pi^{h}:\mathcal{S}\to\mathcal{P}(\left[0,1\right])$, where $\pi^{h}(\cdot|s)$ denotes the probability distribution over $\left[0,1\right]$ given the state $s$. Second, the lower-level policy learns to select an ordered subset with the size determined by the higher-level policy. We define the lower-level policy by $\pi^{l}:\mathcal{S}\times\left[0,1\right]\to\mathcal{P}(\mathcal{A})$, where $\pi^{l}(\cdot|s,k)$ denotes the probability distribution over the action space given the state $s$ and the ratio $k$. Specifically, we formulate the lower-level policy as a sequence model, which can capture the interaction among cuts. Finally, we derive the cut selection policy via the law of total probability, i.e., $\pi(a_{k}|s)=\mathbb{E}_{k\sim\pi^{h}(\cdot|s)}[\pi^{l}(a_{k}|s,k)],$ where $k$ denotes the given ratio and $a_{k}$ denotes the action. The policy is computed by an expectation, as $a_{k}$ cannot determine the ratio $k$. For example, suppose that $N=100$ and the length of $a_{k}$ is $10$, then the ratio $k$ can be any number in the interval $[0.1,0.11)$. Actually, we sample an action from the policy $\pi$ by first sampling a ratio $k$ from $\pi^{h}$ and then sampling an action from $\pi^{l}$ given the ratio.

3. 3.

   Training: Hierarchical Policy Gradient For the cut selection task, we aim to find $\theta$ that maximizes the expected reward over all trajectories

   $\displaystyle J(\theta)=\mathbb{E}_{s\sim\mu,a_{k}\sim\pi_{\theta}(\cdot|s)}[r(s,a_{k})],$

   (5)

   where $\theta=\left[\theta_{1},\theta_{2}\right]$ with $\left[\theta_{1},\theta_{2}\right]$ denoting the concatenation of the two vectors, $\pi_{\theta}(a_{k}|s)=\mathbb{E}_{k\sim\pi^{h}_{\theta_{1}}(\cdot|s)}[\pi^{l}_{\theta_{2}}(a_{k}|s,k)]$, and $\mu$ denotes the initial state distribution. To train the policy with a hierarchical structure, we derive a hierarchical policy gradient.

   Proposition 1

   Given the cut selection policy $\pi_{\theta}(a_{k}|s)=\mathbb{E}_{k\sim\pi^{h}_{\theta_{1}}(\cdot|s)}[\pi^{l}_{\theta_{2}}(a_{k}|s,k)]$ and the training objective (5), the hierarchical policy gradient takes the form of

   	$\displaystyle\nabla_{\theta_{1}}J(\left[\theta_{1},\theta_{2}\right])$	$\displaystyle=\mathbb{E}_{s\sim\mu,k\sim\pi^{h}_{\theta_{1}}(\cdot|s)}[\nabla_{\theta_{1}}\log(\pi^{h}_{\theta_{1}}(k|s))\mathbb{E}_{a_{k}\sim\pi^{l}_{\theta 2}(\cdot|s,k)}[r(s,a_{k})]],$
   	$\displaystyle\nabla_{\theta_{2}}J(\left[\theta_{1},\theta_{2}\right])$	$\displaystyle=\mathbb{E}_{s\sim\mu,k\sim\pi^{h}_{\theta_{1}}(\cdot|s),a_{k}\sim\pi^{l}_{\theta_{2}}(\cdot|s,k)}[\nabla_{\theta_{2}}\log\pi^{l}_{\theta_{2}}(a_{k}|s,k)r(s,a_{k})].$

   We use the derived hierarchical policy gradient to update the parameters of the higher-level and lower-level policies. We implement the training algorithm in a parallel manner that is closely related to the asynchronous advantage actor-critic (A3C) (Mnih et al., 2016).

Please refer to Wang et al. (2023d) for more details.

Deployment and Performance Gain We deploy HEM to large-scale real-world production planning problems from Huawei. The results in Table 4 show that HEM significantly outperforms all the baselines in terms of the Time and PD integral. The results demonstrate the strong ability to enhance the Optverse solver with our proposed HEM in real-world applications.

Finding a good feasible solution in a short time is a very challenging task in solving Mixed Integer Linear Programming (MILP) problems. One of the most promising methods for leveraging machine learning to solve it is a Neural Diving heuristic Nair et al. (2020). Without loss of generality, we will restrict the description only to the binary MILPs. The main idea behind the Neural Diving method is to learn to generate correct partial assignments to selected variables of a given problem in order to solve corresponding sub-problems. Providing that Neural Diving correctly predicts these assignments the problem could be solved much faster.

A deep neural network is used to produce multiple partial assignments of the binary variables of the original MILP problem. The remaining unassigned variables define smaller “sub-problems” (sub-MILPs), which are solved using an off-the-shelf solver to obtain solutions. This method can be easily parallelized since every sub-MILP can be solved independently. Neural Diving can indeed help to quickly get good solutions, though it does not guarantee global optimality.

Following Gasse et al. (2019b), we encoded each MILP as a bi-partite graph with two sets of vertices, representing its variables and constraints, and gave it as input to a Graph Convolutional Neural Network (GCNN) that can predict the values of binary variables. Such a network can be trained in a supervised manner using the binary cross-entropy loss. We enhanced GCNN architecture with layer norms and skip connections to make the training more stable, utilized the primal integral Berthold (2013) for performance evaluation. We observed that performance of the neural diving method is mostly influenced by the following key factors: (1) by the number of parallel sub-MILPs, (2) by their sizes which are determined by the size of partial assignments, and, what is also very important, (3) by the choice of binary variables that constitute these partial assignments. While the first factor is only related to the available hardware resources and the second can be balanced by fine-tuning, the optimal choice of the variables is not so straight forward.

To make a proper choice of the variables, we must take into account the confidence of their prediction and, at the same time, also consider the CPU time which is spent on solving corresponding sub-MILPs. To tackle this trade-off, we developed our own method that is based on ranking the variables according to their ability to reduce the CPU time that would be required to solve the MILP if they were assigned. We proved that good rankings can be learnt off-line and used during the inference to generate faster sub-MILPs. Also we enhanced the GCNN training procedure by using auxiliary losses, one for pumping the feasibility of generated sub-MILPs and the other for better balance of 0 and 1 values in partial assignments. We performed tests on various synthetic and real life datasets which proved that it can provide 3 to 10x speedup compared to the baseline version of OptVerse AI solver (see Figure 10).

Since it is hard for practitioners to explicitly model the relationship between the solver’s performance and corresponding hyper-parameters, performance function of hyper-parameters is actually a black-box function. Hence, Bayesian Optimization (BO), which is a classic black-box optimization technique with high sample efficiency, is widely-used for tuning solvers.

Despite that the BO framework shows availability in tuning tasks for solvers, there still exist some issues in practice Cowen-Rivers et al. (2022). Firstly, Gaussian Process (GP), which is a common surrogate model in BO, typically requires two important modelling assumptions stated as data stationarity and homoscedasticity of the noise distribution. Unfortunately, even simple hyper-parameter tuning tasks can exhibit significant non-stationarity and heteroscedastic. Consequently, a distinct approach is required to process such non-ideal data in real scenario. Secondly, there are various acquisition functions trading off exploration and exploitation to the approximated surrogate model in general BO framework. However, various acquisition functions may result in opposing solutions when pursuing their own optima. Thus, we also need to tackle the conflicting acquisitions.

To address the aforementioned issues, we integrate Heteroscedastic and Evolutionary Bayesian Optimization (HEBO) algorithm which is also developed by ourselves in our parameter tuning framework Cowen-Rivers et al. (2022). This algorithm introduces an Kumaraswamy input warping which performs a transformation to the input data (i.e., hyper-parameter variables) and corrects the non-stationarity. Furthermore, to tackle the data heteroscedasticity, HEBO leverage ideas from warped GP E. et al. (2003) where the Box-Cox transformation is considered as a corrective mapping for non-Gaussian data. As for the conflicting acquisition functions, HEBO designs a multi-objective acquisition function that combines three different commonly-used acquisitions, i.e., Expected Improvement (EI), Probability of Improvement (PI), and Upper Confidence Bound (UCB). Formally, we solve

$$\max\limits_{\mathbf{x}}\left(\bar{\alpha}^{\theta}_{EI}(\mathbf{x}|\mathcal{D}),\bar{\alpha}^{\theta}_{PI}(\mathbf{x}|\mathcal{D}),\bar{\alpha}^{\theta}_{UCB}(\mathbf{x}|\mathcal{D})\right)$$

(6)

where $\bar{\alpha}^{\theta}_{(\cdot)}(\mathbf{x}|\mathcal{D})$ denotes a single acquisition in which $\mathcal{D}$ and $\mathbf{x}$ are historical data (hyper-parameter configurations and corresponding solver performance) and newly sampled points. This comprehensive formulation avoids the hard selection of an optimal acquisition such that the solution will not be dominated by one single acquisition. Correspondingly, HEBO employs the non-dominated sorting genetic algorithm II (NSGA-II) for seeking Pareto-front solutions of the multi-objective acquisitions. With all designs mentioned above, HEBO significantly outperforms existing black-box optimisers on 108 machine learning hyperparameter tuning tasks comprising the Bayesmark benchmark, and win the championship of the NIPS 2020 Black-Box Optimisation challenge Cowen-Rivers et al. (2022).

To further improve sample efficiency on new-coming target tasks, one available approach is to use meta BO to transfer knowledge from related source tasks in between Bai et al. (2023). We also propose an end-to-end meta BO framework termed Transformer BO Maraval et al. (2023) for solver tuning. This framework uses a new transformer-based neural process to directly predict acquisition function values from observed data, without relying on traditional surrogate models like GP. After collecting enough observed data from source tasks, an improved reinforcement learning protocol with inductive bias is adopted to train the Transformer architecture, and then the well-trained meta neural process can be used for tuning target tasks more efficiently with prior knowledge from source tasks.

Differential Evolution (DE) Storn and Price (1997) is a powerful and versatile optimization algorithm that falls under the category of evolutionary algorithms, which are inspired by the process of natural selection and evolution. DE operates on a population of candidate solutions and iteratively refines them to find the optimal solution. It uses mutation, crossover, and selection operations to update the population. The main steps of the DE algorithm are as follows: (1) Initialization: Initialize a population of candidate solutions (vectors) randomly within predefined bounds. (2) Mutation: Create new candidate solutions by perturbing the current population using mutation strategies. In particular, the mutant vector $\mathbf{v}$ may be taken as $\mathbf{v}=\mathbf{a}+F(\mathbf{b}-\mathbf{c})$, where $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{c}$ are three random distinct vectors from the population. (3) Crossover: Change the coordinates of the mutated solution $\mathbf{v}$ by corresponding coordinates of the current vector $\mathbf{x}$ with the cross-over probability $C$ to create trial solutions. (4) Selection: Evaluate the objective function for the trial vector and select the candidates that produce better results, replacing the current population with the best candidates. (5) Termination Criteria: Determine when to stop the algorithm, usually based on a maximum number of iterations, a convergence threshold or the value of the fitness function to be optimized.

DE is known for its ability to efficiently explore the search space and converge to optimal solutions when tackling problems with non-linearity, high-dimensionality, and multimodality. Throughout various competitions organized under the IEEE Congress on Evolutionary Computation conference series, different variants of DE consistently achieve top rankings Das et al. (2016).

The efficiency of the classical implementation of DE depends on three crucial control parameters, namely, the scale factor, crossover rate, and population size. In particular, the scale factor and crossover rate play a pivotal role in determining the algorithm’s behaviour, enabling a fine balance between global and local search, often referred to as exploration and exploitation Črepinšek et al. (2013). A larger population size introduces more diversity into the population and tend to perform better in global search, as they have a higher chance of covering various areas in the search space, although at higher computational cost. Through adjusting the controlling strategies of aforementioned three parameters, some variants of DE can be acquired, which are summarized in Table 5.

To analyze the efficiency of algorithms, we choose a pool of testing functions. They include:

• •

  convex quadratic functions with different large number of dimensions;

• •

  polynomial curve fitting task;

• •

  non-convex functions: Rosenbrock’s, Rastrigin’s, Ackley’s, Griewank’s, Schaffer’s F6, HGBat, Schwefel, Weierstrass.

Based on our numerical tests for chosen pool of test problems we propose an advanced version of DE which we refer to LJADE. Its features are summarized below.

We use rand-to-p-best1 strategy, described in Zhang and Sanderson (2009b), that is expressed as follows:

$$\textbf{rand-to-p-best/1: }\mathbf{v}=\mathbf{a}+F(\mathbf{g_{p}}-\mathbf{a})+F(\mathbf{b}-\mathbf{c}),$$

(7)

where $\mathbf{g_{p}}$ is the random vector from a pool of the top $\mathbf{p}$ individuals in the current population with $\mathbf{p}\in(0,1]$ * 100 %. We also use the self-adaptation scheme for the mutation factor Zhang and Sanderson (2007) to independently select parameters that favour the survival of individuals in case they are not initialized optimally enough, and introduce the modified scheme for the crossover probability Peng et al. (2009) to balance the bias towards small values during self-adaptation process. Furthermore, the population size is dynamically reduced using a linear population size reduction method Tanabe and Fukunaga (2014) and an initialization scheme based on quasi-random Halton sequences Halton (1964) is leveraged, which is somewhat superior to random initialization scheme.

The underlying concept behind these modifications is that the algorithm can adapt its approach based on the rate at which the residual function is decreasing. This adaptability allows it to switch between more aggressive searches and finer local searches as needed.

The results of experiments with several instances from MIPLIB2017 benchmark are shown in Table 6. To optimize hyper-parameters, 1000 trials of tuning are used, while the NNI concurrency is 25. The population size of LJADE is set to 75.