Solving Recurrent MIPs with Semi-supervised Graph Neural Networks
We propose an ML-based model that automates and expedites the solution of
MIPs by predicting the values of variables. Our approach is motivated by the
observation that many problem instances share salient features and solution
structures since they differ only in few (time-varying) parameters. Examples
include transportation and routing problems where decisions need to be
re-optimized whenever commodity volumes or link costs change. Our method is the
first to exploit the sequential nature of the instances being solved
periodically, and can be trained with ``unlabeled'' instances, when exact
solutions are unavailable, in a semi-supervised setting. Also, we provide a
principled way of transforming the probabilistic predictions into integral
solutions. Using a battery of experiments with representative binary MIPs, we
show the gains of our model over other ML-based optimization approaches.
Konstantinos Benidis, Ugo Rosolia, Syama Rangapuram, George Iosifidis, Georgios Paschos
# Solving Recurrent MIPs with Semi-supervised Graph Neural Networks

###### Abstract

We propose an ML-based model that automates and expedites the solution of MIPs by predicting the values of variables. Our approach is motivated by the observation that many problem instances share salient features and solution structures since they differ only in few (time-varying) parameters. Examples include transportation and routing problems where decisions need to be re-optimized whenever commodity volumes or link costs change. Our method is the first to exploit the sequential nature of the instances being solved periodically, and can be trained with "unlabeled" instances, when exact solutions are unavailable, in a semi-supervised setting. Also, we provide a principled way of transforming the probabilistic predictions into integral solutions. Using a battery of experiments with representative binary MIPs, we show the gains of our model over other ML-based optimization approaches.

Graph Neural Networks, Discrete Optimization, Semi-supervised Learning, MIP

## 1 Introduction

The solution of mixed integer programming (MIP) problems is as important as it is challenging to achieve. Indeed, MIP is employed for optimizing real-world operations, decision processes and systems, including transportation [28, 34], facility location [25], production planning [62, 14], content distribution networks [29], all the way to structured predictions [30, 46] and neural network (NN) training [48]. However, despite the success and wide availability of MIP solvers, the fast and accurate optimization of such problems remains largely elusive. This has given rise to perpetual efforts for approximation algorithms [59] and heuristic solutions [37], and for improving the performance of solvers [54].

Nonetheless, an aspect that has received less attention is that several MIP instances share properties and common structures. In fact, practitioners, more often than not, need to solve repeatedly problems that differ only in few parameters, cf. [51, 11]. For example, routing and other transportation problems are re-solved periodically over the same graph whenever new demands or costs are obtained. This motivates the use of ML in order to explore correlations among problem properties and solution values, which in turn can be leveraged to expedite the solution of new instances, enabling their solution even in real-time as an increasing number of applications require [10]. The idea of ML-assisted optimization is not new per se; it has been successfully applied to configure solvers [44], and design heuristics, e.g., in Branch-and-Bound (BnB) techniques [31]. Importantly, recent studies focused on learning problem structures towards _predicting_ (some) variables [31, 19, 52, 39] and/or active constraints [10, 11], with promising results - we review them in Sec. 2.
In this work we make several steps towards the latter direction by developing an ML-based method for facilitating the solution of _binary_ MIPs which (might) exhibit temporal structure. For example, internet traffic follows a diurnal (and not random) pattern [47; 27]; and the same holds for a plenitude of transportation problems [49]. Motivated by these observations, we propose MIPnet, a Long Short-Term Memory (LSTM)-based probabilistic variable-prediction model [56] that operates on a permutation and scale-invariant embedding space created by a Graph Convolutional Network (GCN) [8]. The GCN exploits the intuitive bipartite graph MIP representation [31; 52], towards encoding the problems' salient properties without manual feature engineering.

Furthermore, unlike prior works, MIPnet employs _semi-supervised_ learning to augment the training data without computational costs, while still benefiting from the robustness of supervised learning when possible. And including this unsupervised loss provides additional gains. Namely, not all variables have the same impact on the MIP's objective and on constraint violation. Indeed, misprediction of some variables can render the problem infeasible, while others affect only the objective (and to different extent). Clearly, while a model using only supervised loss cannot discern such conditions, MIPnet can learn these effects with per-variable granularity, and mitigate their impact accordingly. Our goal is to predict a significant portion of the binary variables so as to facilitate the optimization of the remaining (binary and continuous) variables. The selection of these variables is based on a confidence metric following a principled Bayesian approach that enables tunable variable selection without (expensive) model re-training.

MIPnet is evaluated through a battery of experiments, using both real and synthetic traces, and in diverse real-world problems that include network routing, facility location and TSP, among others. These are representative problems in operational research, and have been used in prior works [10; 11; 52; 19] and references therein. The experiments reveal significant gains over these previous works in terms of variable prediction accuracy, while maintaining high constraint feasibility. We also demonstrate the importance of exploiting the temporal dimension of these problems (whenever it is prevalent), towards improving our predictions. Finally, we stress that we focus on binary MILPs which have applications in a vast range of real-world problems;2 yet, our method paves the road towards tackling general MIPs.

Footnote 2: E.g., 164 of 240 benchmarks in [33] are binary/integer LPs; while \(90\%\) of variables are binary in the rest.

In summary, the contributions of this work are:

* We propose MIPnet, a GCN-based probabilistic model that learns to solve MIPs and selects the most confident variables based on a principled Bayesian approach. MIPnet speeds up significantly the solution of _general binary_ MIPs while exploiting, for the first time, any temporal relation across the problem instances.
* MIPnet employs semi-supervised learning to benefit from labeled data when available, and remain operative when these are hard to obtain. This hybrid approach further allows the model to identify and prioritize variables that are pivotal for feasibility and optimality.
* In a series of experiments, MIPnet is proved to find accurate solutions to a range of different problems, consistently outperforming SotA benchmarks such as [10; 52].

**Paper Structure**. Sec. 2 reviews the related work and Sec. 3 introduces the model, training and inference approach. Sec. 4 presents the experimental results and Sec. 5 concludes the study. Model and dataset details, and additional experiments are provided as supplementary material.

## 2 Related work

Learning Configurations & HeuristicsA first thrust of works learn to tune (hyper) parameters of solvers [22; 23; 50; 44; 4; 6]. For example, [22; 23] leverage local search to identify good configurations and [50] developed a pertinent software library. These suggestions are oblivious to, and hence unable to benefit from, the problem structure. A different approach is to learn heuristics. Many works in this context learn how to select BnB variables [1; 21; 31; 2; 35] or BnB nodes [36; 38], and how to run primal heuristics [20]. Similar ideas are used in cutting planes [7] and in optimization over graphs [16; 40; 41]. On the other hand, [12; 45] (and references therein) learn optimization models via alternative decompositions and problem reformulations.

Identifying SubstructuresThe third thrust of related work focuses on learning and exploiting problem substructures. Early efforts include predicting _backdoor variables_[17; 18], i.e., instantiating key variables to increase the problem's tractability. Similarly, assuming the availability of training data, [51; 43] propose sampling methods for predicting active constraints. Along these lines, [10; 11] builds Optimal Classification Trees and Feedforward NNs to predict active constraints and variables. Similarly, [61] assigns variables using a \(k\)-Nearest Neighbours classifier. Finally, [57; 60] use Imitation Learning and RL for splitting ILPs. These works do not employ richer (and more promising) problem representations, nor they design bespoke NN architectures for the problems at hand.
Employing GCNsTo that end, [31; 52; 19] use GCNs where the MIPs are encoded with bipartite graphs. This approach is permutation-invariant, thus can yield generalizable learned policies. In [31] the GCN is used for variable selection in BnB; [52] learns also how to initialize (_Neural Diving_); and [19] uses cross-entropy loss to predict probability distributions for the binary variables so as to expedite branch exploration. On the other hand, [39] proposes a GNN-based unsupervised approach for graph-related problems, where the integral solution is recovered with derandomization through sequential decoding. This interesting approach yields feasible solutions with high probability, yet is not applicable to general MIPs.

Unlike these works (see overview [9]), \(\mathtt{MIPnet}\) considers the _temporal evolution_ of problem's parameters; employs semi-supervised learning to reduce the requirement for training datasets; and follows a Bayesian approach to variable selection through learned distributions with tunable confidence and infeasibility penalization. Appendix J includes further comparisons with prior work.

## 3 Semi-supervised Temporal GCN

### Problem Setting

Consider the binary MILP in standard form:3

Footnote 3: We do not specifically include equality constraints since they can be rewritten as two inequality constraints.

\[\underset{\mathbf{z}}{\text{minimize}} \mathbf{c}^{T}\mathbf{z}\] (1) subject to \[\mathbf{A}\mathbf{z}\preceq\mathbf{b},\]

where \(\mathbf{z}\!=\![\mathbf{z}^{(\mathbf{b})};\mathbf{z}^{(\mathbf{c})}]\), with \(\mathbf{z}^{(\mathbf{b})}\!\in\!\{0,1\}^{D_{x}^{(\mathbf{c})}}\), and \(\mathbf{z}^{(\mathbf{c})}\!\in\mathbb{R}^{D_{x}^{(\mathbf{c})}}\) being the binary and continuous variables, respectively, and \(D_{x}\!=\!D_{x}^{(\mathbf{b})}\!+\!D_{x}^{(\mathbf{c})}\). An instance of (1) is defined by the set of parameters \(\phi\!=\!\{\mathbf{c},\mathbf{A},\mathbf{b}\}\), with \(\mathbf{c}\!\in\!\mathbb{R}^{D_{x}}\), \(\mathbf{A}\!\in\!\mathbb{R}^{D_{x}\times D_{y}}\), \(\mathbf{b}\!\in\!\mathbb{R}^{D_{x}}\) (notation details in Appendix A).

We focus on applications where (1) needs to be solved recurrently with different parameters that follow a temporal structure. Consider a set of time series \(\mathcal{S}\) with \(\mathcal{T}_{s}\) timesteps, where each timestep corresponds to an instance with parameters \(\phi_{s,t}\), with \(s\in\mathcal{S}\) and \(t\in\mathcal{T}_{s}\). For the set of instances \(\{\phi_{s,t}\}_{s\in\mathcal{S},t\in\mathcal{T}_{s}}\), consider their optimal solutions \(\{\mathbf{z}_{s,t}^{\star}\}_{s\in\mathcal{S},t\in\mathcal{T}_{s}}\), with \(\mathcal{S}\subseteq\mathcal{S}\) and \(\mathcal{T}_{s}\subseteq\mathcal{T}_{s}\), i.e., we assume the availability of solution only for a subset of instances. Our goal is to train a model that learns a globally shared mapping from the problem parameters \(\phi_{s,t}\) to a probabilistic representation \(p(\mathbf{z}_{s,t}^{\star})\) of the optimal solution. Leveraging this mapping, we present a method to select a subset of variables for which the optimal assignment is known with high confidence - the size of the subset is a user-defined parameter. Then, we fix these variables in problem (1) and we solve the resulting lower dimensional sub-problem. As shown in the result section, for every new set of instances \(\{\phi_{s,t}^{\prime}\}_{s\in\mathcal{S}^{\prime},t\in\mathcal{T}_{s}^{\prime}}\), the proposed methods is able to expedite the computation of the optimal solution.

We elaborate next on the model architecture, training and inference procedures. For clarity of exposition we drop the time series and timestep subscripts unless necessary.

### Proposed Method

Inspired by [31; 19; 52], we represent an instance of (1) as a bipartite graph by having one node for each variable and constraint, and connecting a variable node to a constraint node iff the variable is active on that constraint. We use a GCN block to map the nodes to an embedding space, followed by an LSTM block to capture the temporal structure of the MILP instances \(\{\phi_{s,t}\}_{s\in\mathcal{S},t\in\mathcal{T}_{s}}\). The model is trained in a semi-supervised manner, exploiting both a supervised and an unsupervised component.

Two key properties of (1) are the scale and permutation invariance, i.e., the solution of the problem does not change if we scale (appropriately) the problem, or permute its parameters and variables. If a model does not satisfy these properties, it would need to be re-trained for each variation of a problem. Our architecture satisfies both scale invariance by introducing a normalization step and permutation invariance by allowing only symmetric functions among nodes.

Parameter NormalizationWe initially normalize all parameters \(\phi\) of (1) in order to make the training more stable and the model scale-invariant. This means that the model will learn to solve the problem with parameters drawn from a "base" distribution, and every scaled problem version can be solved just by normalizing the parameters. The normalization needs to ensure that the relative weights of the objective parameters and each constraint remain the same;
therefore we normalize them separately. The normalized parameters \(\forall i\leq D_{\text{c}},j\leq D_{\text{z}}\), are computed as follows:

\[a_{i,j}\!=\!\frac{a_{i,j}}{\|[\mathbf{a}_{i}^{T};b_{i}]\|_{p}},\ b_{i}\!=\!\frac{ b_{i}}{\|[\mathbf{a}_{i}^{T};b_{i}]\|_{p}},\ c_{j}\!=\!\frac{c_{j}}{\|\mathbf{c}\|_{p}}, \tag{2}\]

where \(\mathbf{a}_{i}\) is the \(i\)-th row of \(\mathbf{A}\). We use \(p=2\). For problem instances of different sizes the normalization can lead to a distribution shift. We discuss this in Appendix B.

Model Architecture.The model builds on the bipartite graph representation of MILPs. Formally, consider a graph \(G=(\mathcal{V},\mathcal{E},\mathbf{A}^{\text{(adj)}})\) defined by the set of nodes \(\mathcal{V}\), with \(|\mathcal{V}|=D_{\text{z}}+D_{\text{c}}\), the set of edges \(\mathcal{E}\), with \(|\mathbf{E}|\) equal to the number of nonzero entries in the constraint matrix \(\mathbf{A}\), and the graph adjacency matrix \(\mathbf{A}^{\text{(adj)}}=[\mathbf{I}_{D_{\text{z}}},\mathbf{A}^{T};\mathbf{ A},\mathbf{I}_{D_{\text{z}}}]\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}\), i.e., the (non-binary) adjacency matrix contains the coefficients from the constraint matrix \(\mathbf{A}\) and self-loops for all nodes given by the identity matrices. One set of \(D_{\text{z}}\) nodes in the bipartite graph corresponds to variables and the other set of \(D_{\text{c}}\) nodes to constraints. An edge \(e_{i,j}\in\mathcal{E}\) indicates that variable \(j\) appears in the \(i\)-th constraint.

The nodes of the graph are associated with a set of features \(\mathbf{U}=\mathbf{U}(\phi)\in\mathbb{R}^{|\mathcal{V}|\times D_{\text{u}}}\) derived by the problem parameters \(\phi\), where \(D_{u}\) is the feature dimension. The feature vector of each node is constructed by linearly combining in a permutation invariant way all the variable and constraint parameters that the node is related to. We refer the reader to Appendix C for a detailed description.

The graph representation of (1) motivates using a GCN [42]. A GCN with \(L\) layers is denoted:

\[\mathbf{X}^{(l+1)}\!=\!g(\mathbf{X}^{(l)},\mathbf{A}^{\text{(adj)}};\theta_{ g}^{(l)}),\ \forall\ l=0,\ldots,L-1, \tag{3}\]

with \(\mathbf{X}^{(0)}=\mathbf{U}\) and \(\theta_{g}^{(l)}\in\Theta\) the learnable parameters of each layer. The adjacency matrix defines the graph connectivity and determines how information is aggregated at each layer, while the number of layers define the number of hops away that a node gets information from. Note that the GCN is applied to each node in parallel and the resulting embedding is by construction permutation invariant. We follow the GCN propagation rule as defined in [42; 8], with a few (optional) modifications between each layer of the GCN as in [52]: (i) We include a nonlinearity after each layer, (ii) we include skip connection inputs and (iii) we apply layer normalization [5] in the output of each layer.

The output of the final layer \(\mathbf{X}=\mathbf{X}^{(L)}\in\mathbb{R}^{|\mathcal{V}|\times D_{\text{z}}}\) is the final embedding of each node which will be used as input to the next NN block of the model architecture. We use an LSTM network to capture the temporal evolution of (1). The inputs to the LSTM at time \(t\) are the embedded nodes \(\mathbf{X}\) as well as the previous network output \(\mathbf{h}_{s,t-1}\), i.e., \(\mathbf{h}_{s,t}=h(\mathbf{h}_{s,t-1},\mathbf{X}_{s,t-1};\theta_{h})\), with \(\theta_{h}\in\Theta\) the learnable parameters of the LSTM. The LSTM retains the permutation invariance since it acts on the feature dimension. The network output at time \(t\), i.e., \(\mathbf{h}_{s,t}\), is appropriately projected to the parameters of the selected probabilistic representation of each variable using a multilayer perceptron (MLP). In general the projection has the following form:

\[\psi_{s,t,j}=f(\mathbf{x}_{s,t,j},\mathbf{h}_{s,t};\theta_{f}), \tag{4}\]

where \(\theta_{f}\in\Theta\) are the learnable MLP parameters and \(\mathbf{x}_{s,t,j}=[\mathbf{X}_{s,t}]_{j}\) are the embedded features of the \(j\)-th variable node that acts as a skip connection, while the constraint nodes are discarded. Note that the mapping layer can apply any function in the feature dimension of \(\mathbf{h}_{s,t}\) but only symmetric functions across the node dimension. Here, \(\psi\) represents a generic set of parameters for the binary variables that may differ based on the selected model output. If (1) includes continuous variables the model outputs their values and not a distribution, i.e., \(\psi_{s,t}\!=\!\{p(\mathbf{a}_{s,t}^{(\mathbf{b})}),\hat{\mathbf{z}}_{s,t}^{( \mathbf{c})}\}\). In particular, we leverage two MLPs to project the output of the LSTM to the estimates of binary and continuous variables (for a discussion about the prediction of continuous variables please see Appendix D). The probabilistic representation of binary variables becomes concrete in Sec. 3.3.

### Training

A common characteristic of MILP applications is the difficulty in obtaining labels. Specifically, to obtain one label, one has to solve an instance of the MILP, which depending on the use case may take from minutes to hours or even days. It is therefore of high interest for the generality of MIPnet to be able to train with a semi-supervised loss, exploiting labeled data and benefiting from unlabeled data.

Supervised Setting.For the instances that have available labels (optimal solutions) we can train the model with an MLE approach. A natural choice for a probabilistic representation of binary variables is the Bernoulli distribution \(\text{Ber}(z;\pi)\)[52; 60], where \(\pi\) is the Bernoulli parameter. Given a Bernoulli representation of a binary variable it is straightforward to assign a value to the variable, but hard to quantify how accurate or confident this prediction is (see Appendix E for a discussion). A different approach was proposed in [52] that used the SelectiveNet [32] and learned
two sets of binary variables. One set indicates if a variable is going to be selected and the other its binary value (to be used if selected). Here, we introduce a principled way to define how accurate or confident is a prediction via its own variance. We leverage this notion of variance and follow a Bayesian approach by using the Beta distribution \(\text{Beta}(\pi;\alpha,\beta)\) to model the Bernoulli parameters, where \(\alpha\) and \(\beta\) are the parameters of the Beta distribution. This approach can be readily used as a statistically principled selection method, and is a key advantage of our method.

Now, we elaborate on the model output (4). The model yields the Beta distribution parameters \(\psi\!=\!\{\alpha,\beta\}\) (we ignore the continuous variables and drop subscripts), and for each variable with optimal value \(z^{\star}\) the likelihood is computed as \(\int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)\text{Ber}(z^{\star};\pi)d\pi\), i.e., we need to integrate over all values of \(\pi\). Unfortunately this integral does not have a closed-form solution; yet, we can approximate it efficiently using a quadrature method (Monte-Carlo sampling can be also used, but is not as efficient). We use the Clenshaw-Curtis quadrature method [15] where the function to be integrated is evaluated at the \(K\) roots of a Chebyshev polynomial and the integral is approximated by a weighted average of the integral's function \(I(\pi)\!=\!\text{Beta}(\pi;\alpha,\beta)\text{Ber}(z^{\star};\pi)\) in specific and predefined points \(\{\bar{\pi}_{k}\}_{k=0}^{K/2}\), i.e.,

\[\int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)\text{Ber}(z^{\star};\pi)d\pi\!=\! \int_{0}^{1}\!I(\pi)d\pi\approx\mathbf{w}^{T}I(\bar{\pi}),\]

where \(\mathbf{w}\in\mathbb{R}^{K/2+1}\) is independent of the function (precomputed) and \(I(\bar{\pi})=[I(\bar{\pi}_{0}),\ldots,I(\bar{\pi}_{K/2})]^{T}\). For details of the method we refer the reader to Appendix F.1. With the above formulation the negative log-likelihood (NLL) of all the instances becomes:

\[\ell_{\text{sup}}(\mathbf{z}^{\star},\psi)\approx-\sum_{s=1}^{|\mathcal{S}|} \sum_{t=1}^{|\mathcal{T}_{s}|}\sum_{j=1}^{D_{s}}\log\left(\mathbf{w}^{T}I( \bar{\pi}_{s,t,j})\right), \tag{5}\]

where \(\psi\!=\!\{\boldsymbol{\alpha}_{s,t},\boldsymbol{\beta}_{s,t}\}\).

In Appendix F.2 we propose an additional regularization term, while in Appendix F.3 we describe a weighted version of 5 for imbalanced labels.

Unsupervised Setting.Apart from (5) we consider an unsupervised loss that is optimized together and is essentially the objective of (1) with a penalty term for the constraint violations:

\[\ell_{\text{unsup}}=\sum_{s=1}^{|\mathcal{S}|}\sum_{t=1}^{|\mathcal{T}_{s}|} \mathbf{c}_{s,t}^{T}\hat{\mathbf{z}}_{s,t}+\lambda_{\text{c}}\|\left(\mathbf{ A}_{s,t}\hat{\mathbf{z}}_{s,t}-\mathbf{b}_{s,t}\right)_{+}\|^{2}, \tag{6}\]

where \(\lambda_{\text{c}}\geq 0\) and \((x)_{+}=\max(x,0)\). Here, \(\hat{\mathbf{z}}_{s,t}\) represents a prediction based on the learned Beta distribution that is computed as \(\hat{z}=\sigma\left(\frac{\alpha}{\alpha+\beta}\right)\), where \(\sigma(\cdot)\) is the sigmoid function and \(\frac{\alpha}{\alpha+\beta}\) is the mean of Beta. Applying the sigmoid to the mean gives values close to the limits that we use as a proxy for the true binary since exact rounding would block the gradient flow. Further, in the unsupervised loss we include the continuous variables that are directly given by the network output (4), i.e., \(\psi_{s,t}=\{\boldsymbol{\alpha}_{s,t},\boldsymbol{\beta}_{s,t},\hat{\mathbf{ z}}_{s,t}^{(c)}\}\) where \(\boldsymbol{\alpha}_{s,t},\boldsymbol{\beta}_{s,t}\in\mathbb{R}^{D^{(b)}_{s}}\) and \(\hat{\mathbf{z}}_{s,t}^{(c)}\in\mathbb{R}^{D^{(s)}_{s}}\).

The unsupervised loss, as explained earlier, allow us to work with instances for which labels are hard to obtain and, at the same time, provides the means to learn the impact of erroneous predictions (with per-variable granularity) on

Figure 1: MIPnet – At each time step \(t\), we construct the features \(\mathbf{U}_{s,t}\) of all nodes based on the MILP temporal parameters, feed them into a GCN that produces the final embedding \(\mathbf{X}_{s,t}\). The embeddings are fed to an LSTM whose output is projected to the set of parameters \(\psi_{s,t}\) and the corresponding loss \(\ell_{s,t}\) is computed.
the objective function and constraint violation. This, in turn, is proved an effective practical control for mitigating the effect of the (few) miss-predictions that MIPnet yields. Note also that the continuous variables are optimized jointly with the binary through this loss.

Putting the above together, the overall loss is given by

\[\ell=\ell_{\text{sup}}+\lambda\ell_{\text{unsup}}, \tag{7}\]

with \(\lambda\geq 0\). Note that if a label is not available, the supervised loss of that instance is masked and we compute only the unsupervised loss. Fig. 1 summarizes the architecture.

### Inference and Variable Selection

Our prediction problem is very challenging since it can return an infeasible solution even if it mispredicts only few variables. Hence, instead of targeting an all-or-nothing solution (as in [11], see Table 4), we apply a _variable selection_ method where we generate a recommended assignment for each variable associated with a notion of confidence. By fixing the value of high-confidence variables and calling an MILP solver on the reduced problem, we speed up the discovery of high-quality solutions, as shown also in [19; 52; 11; 61]. Differently to these works however, we take a formal Bayesian approach into this variable selection problem, which additionally, does not require retraining.

Given a trained model and a new instance \(\phi^{\prime}=\{\mathbf{c},\mathbf{A},\mathbf{b}\}\) the model outputs the parameters \(\psi^{\prime}=\{\boldsymbol{\alpha},\boldsymbol{\beta},\hat{\mathbf{z}}^{( \mathrm{c})}\}\) with \(\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathbb{R}^{D^{(\mathrm{b})}_{z}}\) and \(\hat{\mathbf{z}}^{(\mathrm{c})}\in\mathbb{R}^{D^{(\mathrm{c})}_{z}}\). For each of the \(D^{(\mathrm{b})}_{z}\) binary variables we compute the mean and the standard deviation of its Beta distribution, given by \(\mu_{j}=\frac{\alpha_{j}}{\alpha_{j}+\beta_{j}}\) and \(\sigma_{j}=\sqrt{\frac{\alpha_{j}\beta_{j}}{(\alpha_{j}+\beta_{j})(\alpha_{j} +\beta_{j}+1)}}\). Since a confident output has mean close to 0 or 1 and low variance, we define the score \(s_{j}=\min(\mu_{j},1-\mu_{j})+\gamma\sigma_{j}\), for \(j=1,\ldots,D^{(\mathrm{b})}_{z}\), and we select a desired percentage \(\rho\) of the variables with the lowest score. Parameter \(\gamma\geq 0\) is tuned on a validation set. For a given set of selected variables, we fix their values rounding their mean, and optimize the remaining binary and continuous variables.Following this approach we solve the remaining problem only once, although one has the option to take multiple samples from the learned distribution, solve the remaining problems multiple times and keep the best solution [52].

## 4 Experiments

We evaluate MIPnet using carefully-selected problems and datasets, as detailed in Appendix G. For the routing problem we used the dataset Geant[53] collecting real-world _temporal_ traffic data from Internet Service Providers; for the facility-loc we used the real-world topology from the nobel-germany dataset [53]; and for the caching problem we used actual requests from the MovieLens datasets [24]. We generated synthetic data for the tsp, energy-grid and revenue-max problems. We ran every experiment 3 times in order to evaluate model variance. Appendices H and I include evaluation details and additional experiments, respectively.

### Evaluating Accuracy & Feasibility

We assess the accuracy, optimality and feasibility of MIPnet against various state of the art methods. Namely, we compare our results with Neural Diving (N-div) [52], a GCN model that selects variables using the SelectiveNet approach [32]; the method proposed in [11], an approach that learns a classifier from problem parameters to complete solutions (based on the unique solutions appeared in the training set); and a variant of our model trained only with supervised loss, MIPnet-sup, to evaluate the importance of the unsupervised loss component. As a further ablation, we also evaluate our method when trained only with unsupervised loss which, however, did not achieve satisfactory performance (see Appendix I). As a final method we use the MIP solver SCIP [58]. For all problems expect for routing, SCIP finds an optimal solution for all instances within the \(15\)-minutes time limit. On the other hand for routing, SCIP finds an optimal solution for \(99.73\%\) of instances within the \(60\)-minutes time limit, which highlights the problem's complexity. These optimal solutions - and the suboptimal ones when optimals are not available - serve as _ground truth_, and the results of all other methods are presented relative to them.

For a problem with \(D^{(b)}_{z}\) binary variables, we fix the \(\left\lceil\rho D^{(b)}_{z}\right\rceil\) most-confident ones, where \(\rho\in(0,1)\), and use a solver for the exact optimization of the remaining subproblem, subject to MIPnet's assignments. Table 1 summarizes the percentage of correct predictions on the assigned variables (_accuracy_) across all instances, while Table 2 presents the percentage of infeasible instances. Table 3 compares the optimality gap of the different approaches, which is computed using _only_ the set of feasible solutions (denoted with \(\mathcal{N}_{f}\)), and is defined as the average difference between
a method's objective (obj) and SCIP's optimal value (obj\({}^{*}\)), i.e., \(\texttt{opt\_gap}=\frac{1}{|\mathcal{N}_{f}|}\sum_{\texttt{obj}\in\mathcal{N}_{ f}}(\texttt{obj}-\texttt{obj}^{*})\). We repeat the experiments for different \(\rho\) values.4

Footnote 4: In MIPnet, we can decide a different portion of assigned variables without model retraining. This stands in stark contrast to NI-div where the variable selection percentage is predefined.

We observe that MIPnet matches or outperforms N-div in almost all cases, while it exhibits more consistent behaviour across runs. Compared to MIPnet-sup, it is either on par or better, especially in feasibility. A notable observation is that accuracy itself is not adequate to validate the performance of a model. For example, in the routing and tsp datasets we observe that accuracy values well-above \(90\%\) can lead to (almost) \(100\%\) infeasibility; while in facility-loc accuracy values \(>99\%\) can yield a significant portion of infeasible instances - especially for MIPnet-sup and N-div. This is associated with the complexity of the constraints and highlights the importance of using an unsupervised component in model training.

Since [11] is an approach that selects a complete solution, i.e., there is no option for selecting \(\rho\%\) of variables, we present the results of this method separately. As shown in Table 4, this all-or-nothing approach learns feasible solutions for tsp, facility-loc, and caching (no variance was observed across runs). In these problems the constraints are time-invariant (do not change across instances). Hence, optimal solutions from the training dataset are feasible for the test dataset, and indeed we see that [11] finds feasible solutions. Alas, these are often far from optimal, e.g., for tsp and caching the average gap is \(19.61\%\) and \(52.49\%\), respectively. On the other hand, for revenue-max where the constraints change across instances, optimal solutions from the training dataset are infeasible for new problem instances and the method returns infeasible variable assignments. That is, this approach does not generalize to solutions not seen in the training dataset. Finally, we highlight that we were able to run the approach from [11] only on a smaller version of the routing dataset, which alludes to potential implementation difficulties of this approach - see Appendix I.

\begin{table}
\begin{tabular}{c c|c c c c c c} \hline \hline \(\rho\) (\%) & Method & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline \multirow{3}{*}{30\%} & MIPnet & **97.0\(\pm\)0.13** & **100.00\(\pm\)0.00** & **98.12\(\pm\) 1.38** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & **97.0\(\pm\)0.04** & **100.98\(\pm\)0.00** & 94.56\(\pm\) 1.73 & 99.97\(\pm\) 0.05 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 84.0\(\pm\)0.72 & **100.00\(\pm\)0.00** & 93.43\(\pm\) 0.00 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \multirow{3}{*}{40\%} & MIPnet & **96.03\(\pm\)0.13** & **100.00\(\pm\)0.00** & **97.54\(\pm\) 1.73** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & 80.01\(\pm\)0.15 & 98.97\(\pm\) 1.69 & 93.97\(\pm\) 1.49 & 99.97\(\pm\) 0.00 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 87.17\(\pm\) 3.68 & 99.98\(\pm\) 0.02 & 90.69\(\pm\) 1.49 & 99.20\(\pm\) 1.31 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \multirow{3}{*}{50\%} & MIPnet & **94.46\(\pm\)0.09** & **100.00\(\pm\)0.00** & **96.97\(\pm\) 1.49** & **99.98\(\pm\) 0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & **94.46\(\pm\)0.11** & 99.79\(\pm\) 0.24 & 93.42\(\pm\) 1.49 & 99.94\(\pm\) 0.00 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 89.34\(\pm\) 0.11 & 99.91\(\pm\) 0.03 & 95.92\(\pm\) 0.54 & 99.61\(\pm\) 0.53 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \multirow{3}{*}{60\%} & MIPnet & 92.33\(\pm\)0.08 & **99.01\(\pm\) 90.03** & **96.38\(\pm\) 2.05** & 99.96\(\pm\) 0.01 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & **93.15\(\pm\) 4.11** & 99.71\(\pm\) 0.29 & **93.26\(\pm\) 1.35** & 99.86\(\pm\) 0.13 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 90.38\(\pm\) 2.55 & 99.91\(\pm\) 0.02 & 95.41\(\pm\) 2.61 & **99.96\(\pm\) 0.01** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \multirow{3}{*}{70\%} & MIPnet & 91.25\(\pm\) 0.07 & **93.84\(\pm\) 0.58** & 95.72\(\pm\) 2.85 & **99.89\(\pm\) 0.01** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & **91.86\(\pm\) 0.09** & 99.41\(\pm\) 0.08 & 93.92\(\pm\) 1.19 & 99.74\(\pm\) 0.14 & 100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 90.30\(\pm\) 0.92 & 99.64\(\pm\) 0.06 & 93.93\(\pm\) 0.84 & 99.64\(\pm\) 0.11 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Accuracy (mean \(\pm\) std)\(\%\) of MIPnet vs. N-div. Bold indicates the best method (higher values are better).

\begin{table}
\begin{tabular}{c c|c c c c c c} \hline \hline \(\rho\) (\%) & Method & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline \multirow{3}{*}{30\%} & MIPnet & **97.0\(\pm\)0.13** & **100.00\(\pm\)0.00** & **98.12\(\pm\) 1.38** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & 61.27\(\pm\)0.04** & **100.00\(\pm\)0.00** & 93.81\(\pm\) 0.13 & 99.97\(\pm\) 0.05 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 84.0\(\pm\)0.72 & **100.00\(\pm\)0.00** & 93.43\(\pm\) 0.00 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\ \hline \multirow{3}{*}{40\%} & MIPnet & **96.03\(\pm\)0.13** & **100.00\(\pm\)0.00** & **97.54\(\pm\) 1.73** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & NPnet-sup & 60.01\(\pm\)0.15 & 98.97\(\pm\) 1.06 & 93.97\(\pm\) 1.49 & 99.97\(\pm\) 0.04 & **100.00\(\pm\)0.00** & **100.00\(\pm\)0.00** \\  & -div & 87.17\(\pm\) 3.68 & 99.98\(\pm\) 0.02 & 90.69\(\pm\) 1.
### Evaluating the Temporal Component of MIPnet

In this section we explore the benefit of leveraging the temporal aspect of the problem instances. Figure 2 illustrates the differences in accuracy and infeasibility when MIPnet employs an MLP instead of an LSTM, for three problems with different temporal properties. Namely, the real-world demands in routing follow a diurnal pattern; the demands in facility-loc have a strong time-dimension; while the edge costs in tsp are generated by adding random perturbations in successive instances.

We observe that the LSTM-based MIPnet outperforms consistently, and most-often substantially, its MLP-based version both in terms of mean and variance values. In other words, all else being equal in the model and inference approach, the addition of LSTM has substantial gains, without requiring any compromises (MLP never outperforms LSTM). This finding underlines the importance of leveraging the temporal evolution pattern of the problem parameters, a hitherto overlooked aspect in all prior works.

### Trading Off Accuracy and Solution Speed via \(\rho\)

Finally, we explore how the percentage of assigned variables \(\rho\), affects the solution time. These experiments reveal also the gains of MIPnet in terms of solution speed, compared to off-the-shelf solvers. In detail, we denote with \(t_{p}\) the time the SCIP solver requires for solving the optimization problem, where \(p=1-\rho\) is the percentage of variables that are _not_ assigned by MIPnet (hence, need to be optimized by the solver). Figure 3 reports the normalized time defined as the ration \(t_{p}/t_{100}\), i.e., the ratio between the solver (SCIP) time when \(\rho\) variables are assigned using MIPnet, and the time for solving the entire problem with SCIP.5 As expected, the time decreases when a higher percentage of variables is assigned. For all problems an average of, at least, \(2\times\) speedup is obtained by setting \(\rho=0.7\), while for revenue-max, energy-grid, and caching we achieve a 2x speedup already for \(\rho=0.3\), i.e., when assigning only \(30\%\) of the variables, and a maximum \(5\times\) speed up. These experiments manifest the solution speed gains of MIPnet compared to a solver, and the importance of our variable selection threshold which indeed can be used to trade off accuracy with solution speed.

Footnote 5: These times are calibrated in order to account for the different computational environments we have used in the experiments.

## 5 Conclusions

Solving MIPs via ML methods can revolutionize the solution of large-scale NP-hard problems and avail MIP to new application domains. This potential is hampered by unprecedented challenges: the values of variables are inherently correlated and need to be jointly predicted; even tiny value assignment errors may lead to infeasibility or unbounded

\begin{table}
\begin{tabular}{c|c|c|c} \hline \hline  & **Accuracy** & **Infeasibility** & **Optimality gap** \\ \hline routing & - & - & - \\ facility-loc & 98.13 & 0.00 & 1.78 \\ tsp & 86.29 & 0.00 & 19.61 \\ revenue-max & 99.22 & 100.00 & - \\ caching & 99.91 & 0.00 & 52.49 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Accuracy, infeasibility and optimality gap of [11].

\begin{table}
\begin{tabular}{c c|c c c c c c} \hline \hline \(\rho\) (\%) & Method & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline \multirow{3}{*}{30\%} & MIPnet & \(6.59\pm 0.15\) & \(0.07\pm 0.01\) & \(\mathbf{0.24\pm 0.02}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & MIPnet-sup & \(\mathbf{6.03\pm 0.63}\) & \(\mathbf{0.04\pm 0.05}\) & \(2.38\pm 0.00\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & S-div & & \(0.28\pm 0.23\) & \(0.25\pm 0.04\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \hline \multirow{3}{*}{40\%} & MIPnet & \(\mathbf{6.27\pm 0.00}\) & \(\mathbf{16.0\pm 0.01}\) & \(0.37\pm 0.07\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & MIPnet-sup & \(27.59\pm 0.00\) & \(0.20\pm 0.10\) & \(\mathbf{0.20\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & N-div & - & \(0.44\pm 0.35\) & \(\mathbf{0.26\pm 0.04}\) & \(0.08\pm 0.04\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \hline \multirow{3}{*}{50\%} & MIPnet & - & \(\mathbf{0.34\pm 0.03}\) & \(\mathbf{0.74\pm 0.19}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & MIPnet-sup & - & \(0.72\pm 0.42\) & \(0.65\pm 0.00\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & N-div & - & \(3.55\pm 4.88\) & \(5.70\pm 0.63\) & \(8.81\pm 1.245\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \hline \multirow{3}{*}{60\%} & MIPnet & - & \(\mathbf{0.73\pm 0.10}\) & \(1.49\pm 0.30\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & MIPnet-sup & - & \(2.15\pm 1.72\) & \(4.18\pm 0.00\) & \(0.01\pm 0.01\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & N-div & - & \(1.22\pm 0.55\) & \(\mathbf{0.48\pm 0.01}\) & \(0.01\pm 0.00\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \hline \multirow{3}{*}{70\%} & MIPnet & - & \(\mathbf{266\pm 0.72}\) & \(2.08\pm 0.17\) & \(\mathbf{0.01\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\  & MPnet-sup & - & \(5.80\pm 5.04\) & \(\mathbf{0.01\pm 0.01}\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \cline{1-1}  & S-div & - & \(5.83\pm 1.97\) & \(\mathbf{1.79\pm 0.21}\) & \(0.04\pm 0.03\) & \(\mathbf{0.00\pm 0.00}\) & \(\mathbf{0.00\pm 0.00}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Optimality gap (mean \(\pm\) std)\(\%\) of MIPnet vs. N-div. Bold indicates the best method (lower values are better).
suboptimality; training data are not readily available; and generalizing trained models is highly non-trivial. In this work, we propose MIPnet, a systematic framework for tackling general binary MILPs that may exhibit temporal structure, which: employs GCNs to automate embeddings and capture correlation among variables; uses an LSTM to account for, hitherto-overlooked, time dependencies across problem instances (a regularly-encountered operational condition); and includes a principled Bayesian assignment strategy with user-tunable confidence. We follow a semi-supervised approach so as to account for the lack of solved instances, a practical limitation when dealing with NP-hard problems, and in order, also, to identify key variables impacting feasibility and optimality. Our work fills (some of) the gaps in the literature and, we believe, contributes towards establishing ML as an off-the-shelf MIP technology, including for non-binary and nonlinear problems that were not addressed by the proposed approach.

Figure 3: Calibrated running time for various \(\rho\) values. \(t_{p}\) is the solution time when a percentage \(\rho\) of variables is assigned by MIPnet and \(p\!=\!1\!-\!\rho\) variables are assigned by the solver. \(t_{100}\) is the time when using only the solver.

Figure 2: Comparison of MIPnet wl and wo a temporal component.
## References

* [1] A. H. Land, and A. G. Doig. An Automatic Method of Solving Discrete Programming Problems. _Econometrica_, 28(3):497-520, 1960.
* [2] A. M. Alvarez. A Machine Learning-based Approximation of Strong Branching. _INFORMS Journal on Computing_, 29(1):185-195, 2017.
* [3] Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. _Advances in Neural Information Processing Systems (NeurIPS)_, 2019.
* [4] Carlos Ansotegui, Yuri Malitsky, Horst Samulowitz, Meinolf Sellmann, and Kevin Tierney. Model-based genetic algorithms for algorithm configuration. In _Twenty-Fourth International Joint Conference on Artificial Intelligence_, 2015.
* [5] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. _arXiv preprint arXiv:1607.06450_, 2016.
* [6] Maria-Florina Balcan, Travis Dick, Tuomas Sandholm, and Ellen Vitercik. Learning to branch. In _International Conference on Machine Learning (ICML)_, pages 344-353, 2018.
* [7] Maria-Florina F Balcan, Siddharth Prasad, Tuomas Sandholm, and Ellen Vitercik. Sample complexity of tree search configuration: Cutting planes and beyond. _Advances in Neural Information Processing Systems (NeurIPS)_, 2021.
* [8] Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al. Relational inductive biases, deep learning, and graph networks. _arXiv preprint arXiv:1806.01261_, 2018.
* [9] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost. Machine learning for combinatorial optimization: a methodological tour d'horizon. _European Journal of Operational Research_, 290(2):405-421, 2021.
* [10] Dimitris Bertsimas and Bartolomeo Stellato. Online Mixed-integer Optimization in Milliseconds. _arXiv preprint arXiv:1907.02206_, 2019.
* [11] Dimitris Bertsimas and Bartolomeo Stellato. The Voice of Optimization. _Machine Learning_, 110(2):249-277, 2021.
* [12] Pierre Bonami, Andrea Lodi, and Giulia Zarpellon. Learning a classification of mixed-integer quadratic programming problems. In _International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research_, pages 595-604. Springer, 2018.
* [13] Francesco Borrelli, Alberto Bemporad, and Manfred Morari. _Predictive control for linear and hybrid systems_. Cambridge University Press, 2017.
* [14] Zhi-Long Chen. Integrated production and outbound distribution scheduling: review and extensions. _Operations research_, 58(1):130-148, 2010.
* [15] Charles W Clenshaw and Alan R Curtis. A method for numerical integration on an automatic computer. _Numerische Mathematik_, 2(1):197-205, 1960.
* [16] Hanjun Dai, Bo Dai, and Le Song. Discriminative embeddings of latent variable models for structured data. In _International conference on machine learning_, pages 2702-2711, 2016.
* [17] Bistra Dilkina, Carla P Gomes, Yuri Malitsky, Ashish Sabharwal, and Meinolf Sellmann. Backdoors to combinatorial optimization: Feasibility and optimality. In _International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research_, pages 56-70. Springer, 2009.
* [18] Bistra Dilkina, Carla P Gomes, and Ashish Sabharwal. Backdoors in the context of learning. In _International Conference on Theory and Applications of Satisfiability Testing_, pages 73-79. Springer, 2009.
* [19] Jian-Ya Ding, Chao Zhang, Lei Shen, Shengyin Li, Bing Wang, Yinghui Xu, and Le Song. Accelerating primal solution findings for mixed integer programs based on solution prediction. In _Proceedings of the AAAI Conference on Artificial Intelligence_, volume 34, pages 1452-1459, 2020.
* [20] E. B. Khalil, B. Dilkina, G. L. Nemhauser, S. Ahmed, and Y. Shao. Learning to Run Heuristics in Tree Search. In _International Joint Conference on Artificial Intelligence (IJCAI)_, 2017.
* [21] E. B. Khalil, et al. Learning to Branch in Mixed Integer Programming. In _Proceedings of the AAAI conference on artificial intelligence_, page 724-731, 2016.
* [22] F. Hutter, et al. ParamILS: An Automatic Algorithm Configuration Framework. _Journal of Artificial Intelligence Research_, 36(1), 2009.
* [23] F. Hutter, et al. Sequential Model-based Optimization for General Algorithm Configuration. _Learning and intelligent optimization_, pages 507-523, 2011.
* [24] F. M. Harper, and J. A. Konstan. The MovieLens Datasets: History and Context. _ACM Transactions on Interactive Intelligent Systems_, 5:19:1-19:2, 2015.
* [25] R. Z. Farahani and M. Hekmatfar. Facility Location: Concepts, Models, Algorithms and Case Studies. _Springer_, 2009.
* [26] Nuno P. Fafsca, Vivek Dua, and Efstratios N. Pistikopoulos. _Multiparametric Linear and Quadratic Programming_, chapter 1, pages 1-23. John Wiley & Sons, Ltd, 2007. ISBN 9783527631216. doi: [https://doi.org/10.1002/9783527631216.ch1](https://doi.org/10.1002/9783527631216.ch1). URL [https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527631216.ch1](https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527631216.ch1).
* [27] G. Barlacchi, et al. A Multi-source Dataset of Urban Life in the City of Milan and the Province of Trentino. _Scientific Data_, 2, 2015.
* [28] G. Laporte. Fifty Years of Vehicle Routing. _Transportation Science_, pages 408-416, 2009.
* [29] G. S. Paschos, et al. Cache Optimization Models and Algorithms. _Found. Trends Commun. Inf. Theory_, 16(3-4):156-343, 2019.
* [30] G. Vivek, K. Srikumar, and D. Roth. On amortizing inference cost for structured prediction. _Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning_, pages 1114-1124, 2012.
* [31] Maxime Gasse, Didier Chetelat, Nicola Ferroni, Laurent Charlin, and Andrea Lodi. Exact combinatorial optimization with graph convolutional neural networks. _Advances in Neural Information Processing Systems (NeurIPS)_, 2019.
* [32] Yonatan Geifman and Ran El-Yaniv. Selectivenet: A deep neural network with an integrated reject option. In _International Conference on Machine Learning (ICML)_, pages 2151-2159, 2019.
* [33] Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bastubbe, Timo Berthold, Philipp M. Christophel, Kati Jardx, Thorsten Koch, Jeff Linderoth, Marco Lubbecke, Hans D. Mittelmann, Derya Ozyurt, Ted K. Ralphs, Domenico Salvagnin, and Yuji Shinano. MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library. _Mathematical Programming Computation_, 2021. doi: 10.1007/s12532-020-00194-3. URL [https://doi.org/10.1007/s12532-020-00194-3](https://doi.org/10.1007/s12532-020-00194-3).
* [34] B. Gopalakrishnan and E. L. Johnson. Airline Crew Scheduling: State-of-the-Art. _Annals of Operations Research_, 140(1):305-337, 2005.
* [35] Christoph Hansknecht, Imke Joormann, and Sebastian Stiller. Cuts, primal heuristics, and learning to branch for the time-dependent traveling salesman problem. _arXiv preprint arXiv:1805.01415_, 2018.
* [36] He He, Hal Daume III, and Jason M Eisner. Learning to search in branch and bound algorithms. _Advances in Neural Information Processing Systems (NeurIPS)_, 2014.
* [37] I. Boussaid, et al. A Survey on Optimization Metaheuristics. _Information Sciences_, 237:82-117, 2013.
* [38] J. Song, R. Lanka, A. Zhao, Y. Yue, and M. Ono. Learning to Search via Retrospective Imitation. In _arXiv:1804.00846_, 2018.
* [39] Nikolaos Karalias and Andres Loukas. Erdos goes neural: an unsupervised learning framework for combinatorial optimization on graphs. In _Proceedings of NeurIPS_, 2020.
* [40] Elias Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Le Song. Learning combinatorial optimization algorithms over graphs. _Advances in Neural Information Processing Systems (NeurIPS)_, 2017.
* [41] Minsu Kim, Jinkyoo Park, et al. Learning collaborative policies to solve np-hard routing problems. _Advances in Neural Information Processing Systems (NeurIPS)_, 2021.
* [42] Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. _arXiv preprint arXiv:1609.02907_, 2016.
* [43] Martin Klauco, Martin Kaluz, and Michal Kvasnica. Machine learning-based warm starting of active set methods in embedded model predictive control. _Engineering Applications of Artificial Intelligence_, 77:1-8, 2019.
* [44] Robert Kleinberg, Kevin Leyton-Brown, Brendan Lucier, and Devon Graham. Procrasinating with confidence: Near-optimal, anytime, adaptive algorithm configuration. _Advances in Neural Information Processing Systems (NeurIPS)_, 2019.
* [45] Markus Kruber, Marco E Lubbecke, and Axel Parmentier. Learning when to use a decomposition. In _International conference on AI and OR techniques in constraint programming for combinatorial optimization problems_, pages 202-210. Springer, 2017.
* [46] Alex Kulesza and Fernando Pereira. Structured learning with approximate inference. _Advances in Neural Information Processing Systems (NeurIPS)_, 2007.
* [47] Anukool Lakhina, Konstantina Papagiannaki, Mark Crovella, Christophe Diot, Eric D. Kolaczyk, and Nina Taft. Structural analysis of network traffic flows. In _Proceedings of ACM Sigmetrics_, 2004.
* [48] Cong Leng, Zesheng Dou, Hao Li, Shenghuo Zhu, and Rong Jin. Extremely low bit neural network: Squeeze the last bit out with ADMM. In _Proceedings of the AAAI Conference on Artificial Intelligence_, 2018.
* [49] M. C. Gonzalez, et al. Understanding Individual Human Mobility Patterns. _Nature_, 453:779-782, 2008.
* [50] M. Lopez-Ibanez, et al. The Trace Package: Iterated Racing for Automatic Algorithm Configuration. _Operations Research Perspectives_, pages 43-58, 2016.
* [51] Sidhant Misra, Line Roald, and Yeesian Ng. Learning for constrained optimization: Identifying optimal active constraint sets. _INFORMS Journal on Computing_, 34(1):463-480, 2022.
* [52] Vinod Nair, Sergey Bartunov, Felix Gimeno, Ingrid von Glehn, Pawel Lichocki, Ivan Lobov, Brendan O'Donoghue, Nicolas Sonnerat, Christian Tjandraatmadja, Pengming Wang, et al. Solving mixed integer programs using neural networks. _arXiv preprint arXiv:2012.13349_, 2020.
* [53] S. Orlowski, M. Pioro, A. Tomaszewski, and R. Wessaly. SNDlib 1.0-Survivable Network Design Library. In _Proceedings of the 3rd International Network Optimization Conference (INOC 2007), Spa, Belgium_, April. [http://sndlib.zib.de](http://sndlib.zib.de), extended version accepted in Networks, 2009.
* [54] R. E. Bixby. A Brief History of Linear and Mixed-Integer Programming Compucation. _Documenta Mathematica_, pages 107-121, 2010.
* [55] S. Chopra, I. Gilboa, and S. T. Sastry. Source Sink Flows with Capacity Installation in Batches. _Discrete Applied Mathematics_, 85:165-192, 1998.
* [56] S. Hochreiter, et al. Long short-term memory. _Neural Computation_, 9(8):1735-1780, 1997.
* [57] Jialin Song, Yisong Yue, Bistra Dilkina, et al. A general large neighborhood search framework for solving integer linear programs. _Advances in Neural Information Processing Systems (NeurIPS)_, 2020.
* [58] T. Achterberg. SCIP: Solving Constraint Integer Programs. _Mathematical Programming Computation_, 1:1-41, 2009.
* [59] David Williamson and David Shmoys. _The Design of Approximation Algorithms_. Cambridge University Press, 2011.
* [60] Yaoxin Wu, Wen Song, Zhiguang Cao, and Jie Zhang. Learning large neighborhood search policy for integer programming. _Advances in Neural Information Processing Systems (NeurIPS)_, 2021.
* [61] Alinson S Xavier, Feng Qiu, and Shabbir Ahmed. Learning to solve large-scale security-constrained unit commitment problems. _INFORMS Journal on Computing_, 33(2):739-756, 2021.
* [62] Y. Pochet, and L. A. Wolsey. Production Planning by Mixed Integer Programming. _Springer Science and Business Media_, 2006.
**Supplementary Materials**

## Appendix A Notation

Matrices, vectors and scalars are denoted by uppercase bold, lowercase bold and lowercase normal letters, i.e., \(\mathbf{X}\), \(\mathbf{x}\) and \(x\), respectively. The set of \(D\)-dimensional real numbers is defined by \(\mathbb{R}^{D}\), \(\mathbf{x}\in\mathbb{R}^{D}\) and \(\mathbf{z}\in\mathbb{Z}^{D}\) denote vectors of real and integer variables with dimension \(D\). \([\mathbf{x};\mathbf{z}]\) and \([\mathbf{x},\mathbf{z}]\) indicates column and row concatenation, respectively. \(\mathbf{X}^{T}\) denotes the transpose of matrix \(\mathbf{X}\) and \([\mathbf{X}]_{j}\) its \(j\)-th row. \(\mathbf{I}_{D}\) and \(\mathbf{0}_{D}\) denote the \(D\)-dimensional identity and zero matrix, respectively. The number of variables is denoted as \(D_{z}\) and the number of constraints as \(D_{\mathrm{c}}\). We use a hat to indicate a variable prediction, e.g., \(\hat{z}\). We use the subscript \(j\) as a variable index and \(i\) as a constraint index. We use uppercase calligraphic and Greek letters to define sets. We denote by \(\mathcal{S}\) a set of time series (indexed by \(s\)) and by \(\mathcal{T}_{s}=\{1,\ldots,|\mathcal{T}_{s}|\}\) the set of timesteps of time series \(s\) (indexed by \(t\)). \(\phi\) represents the set of parameters of an optimization problem and \(\Theta\) the set of parameters of a NN. The symbol \(\mathbbm{1}_{\text{condition}}\) denotes the indicator function with value 1 if condition is true and 0 otherwise. \(\text{KL}\left(p||q\right)\) denotes the KL-divergence of distributions \(p\) and \(q\).

## Appendix B Parameter Normalization & Distribution Shift

As mentioned in Section 3, the parameter normalization (see Eq. (2)) makes our model scale invariant, i.e., our trained model can be applied to new instances where the magnitudes of the constraint matrix and cost vector are entirely different from what was seen during training. More importantly, our trained model can also be used in solving MIP problems whose sizes (i.e., length of the cost vector \(D_{z}\)) are different from that of training instances. However, during test time, if the problem size grows (or shrinks), the parameters will become smaller (or larger), causing a distribution shift.

The distribution shift is caused since in training time with the normalization given in (2) (with \(p=2\)) the parameters \(a_{ij}\) and \(b_{i}\) are proportional to \(\frac{1}{\sqrt{D_{z}+1}}\), and \(c_{j}\) is proportional to \(\frac{1}{\sqrt{D_{z}}}\), whereas at test time these would be proportional to \(\frac{1}{\sqrt{D_{z}+1}}\) and \(\frac{1}{\sqrt{D_{z}}}\), respectively. This can be rectified by a further rescaling of parameters during test time.

More precisely, given a model that is trained on instances of size \(D_{z}\) and applied on test instances of size \(\bar{D}_{z}\), we normalize the test instance's parameters as

\[a_{i,j}\!=\!\sqrt{\frac{\bar{D}_{z}+1}{D_{z}+1}}\frac{a_{i,j}}{\|[\mathbf{a} _{i}^{T};b_{i}]\|_{2}},\quad b_{i}\!=\!\sqrt{\frac{\bar{D}_{z}+1}{D_{z}+1}} \frac{b_{i}}{\|[\mathbf{a}_{i}^{T};b_{i}]\|_{2}},\quad c_{j}\!=\!\sqrt{\frac{ \bar{D}_{z}}{D_{z}}}\frac{c_{j}}{\|\mathbf{c}\|_{2}}. \tag{8}\]

This reduces to Eq. (2) (with \(p=2\)), if \(\bar{D}_{z}=D_{z}\). With this normalization, the parameters \((a_{i,j},b_{i})\) and \(c_{j}\) during test time are also proportional to \(\frac{1}{\sqrt{D_{z}+1}}\) and \(\frac{1}{\sqrt{D_{z}}}\) respectively.

Given this normalization, our model can also be trained on MIP instances of different sizes. We just need to fix one size \(D_{z}\) as the reference size and for every instance of different size \(\bar{D}_{z}\), the parameters would have a further rescaling as in (8).

## Appendix C Graph Features

An integral part of MIPnet is the creation of the features of each node that form the matrix \(\mathbf{U}\in\mathbb{R}^{|\mathcal{V}|\times D_{u}}\), which is fed into the GCN (see Eq. (3)), where \(D_{u}\) is the feature dimension. We differentiate the feature construction between variable and constraint nodes.

The set of parameters that the \(j\)-th variable is associated with are: (i) its cost in the objective \(c_{j}\), (ii) all its parameters in the constraint matrix \(\mathbf{A}\), i.e., \(\mathbf{a}_{j}=[\mathbf{A}^{T}]_{j}\) and (iii) the right-hand side of the constraints \(\mathbf{b}\). Usually the constraint matrix \(\mathbf{A}\) is sparse, especially for large problems, which means that each variable participates in only a small subset of the constraints. We leverage this observation and we include as relevant features only the subset of constraints that each variable is included, i.e., for the \(j\)-th variable we include the constraint \(i\) if \(a_{i,j}\neq 0\). In particular, we calculate the maximum number of constraints that any variable appears, \(m_{c}\), and we keep that many constraint parameters for all variables (so all the feature dimensions are consistent).6
Now, for the \(j\)-th variable we construct \(m_{c}\) triplets \((a_{i,j},b_{i},c_{j})\), where \(i\in\mathcal{I}_{j}=\{i|a_{i,j}\neq 0\}\). If a variable is included in less that \(m_{c}\) constraints, i.e., if \(|\mathcal{I}_{j}|<m_{c}\), then we include \(m_{c}-|\mathcal{I}_{j}|\) zero triplets. Each triplet is mapped to the target feature dimension \(D_{u}\) using an MLP with parameters learned jointly across all variables:

\[\mathbf{V}_{j}=f^{\text{(var-map)}}(\mathbf{a}_{\mathcal{I}_{j},j},\mathbf{b} _{\mathcal{I}_{j}},c_{j})\!=\!\begin{bmatrix}a_{i_{1},j}&b_{i_{1}}&c_{j}\\ \vdots&\vdots&\vdots\\ a_{i_{m_{c}},j}&b_{i_{m_{c}}}&c_{j}\end{bmatrix}\mathbf{W}^{\text{(var-map)}},\]

where \(\mathbf{W}^{\text{(var-map)}}\in\mathbb{R}^{3\times D_{u}}\) is the learnable weight matrix of the MLP and \(\mathbf{V}_{j}\in\mathbb{R}^{m_{c}\times D_{u}}\). Finally, we aggregate over all rows using a symmetric function (e.g., average):

\[\mathbf{u}_{j}^{\text{(var)}}=f^{\text{(var-agg)}}(\mathbf{V}_{j}), \tag{9}\]

where \(\mathbf{u}_{j}^{\text{(var)}}\in\mathbb{R}^{D_{u}}\). Note that \(\mathbf{u}_{j}^{\text{(var)}}\) is variable and constraint permutation invariant. The complete feature matrix for the variable nodes is constructed by concatenating all the feature vectors \(\mathbf{u}_{j}^{\text{(var)}}\), i.e., \(\mathbf{U}^{\text{(var)}}=[\mathbf{u}_{1}^{\text{(var)}},\dots,\mathbf{u}_{ D_{x}}^{\text{(var)}}]\in\mathbb{R}^{D_{x}\times D_{u}}\).

The features of the constraints are constructed in a similar manner. The set of parameters that the \(i\)-th constraint is associated with are: (i) all its parameters in the constraint matrix \(\mathbf{A}\), i.e., \(\mathbf{a}_{i}=[\mathbf{A}]_{i}\) (ii) the right-hand side of the constraint \(b_{i}\) and (iii) the costs in the objective \(\mathbf{c}\) for the variables in the constraint. As in the variable nodes, we leverage the sparsity of \(\mathbf{A}\) and find the maximum number of variables in all the constraints, \(m_{v}\). Now, for the \(i\)-th constraint we create \(m_{v}\) triplets \((a_{i,j},b_{i},c_{j})\), where \(j\in\mathcal{J}_{i}=\{j:a_{i,j}\neq 0\}\). The triplets are projected to the target dimension \(D_{u}\) and then aggregated in the same way as the variable features:

\[\mathbf{C}_{i}=f^{\text{(con-map)}}(\mathbf{a}_{i,\mathcal{J}_{i }},b_{i},\mathbf{c}_{\mathcal{J}_{i}}) \!=\!\begin{bmatrix}a_{i,j_{1}}&b_{i}&c_{j_{1}}\\ \vdots&\vdots&\vdots\\ a_{i,j_{m_{v}}}&b_{i}&c_{j_{m_{v}}}\end{bmatrix}\mathbf{W}^{\text{(con-map)}},\] \[\mathbf{u}_{i}^{\text{(con)}}=f^{\text{(con-agg)}}(\mathbf{C}_{i}),\]

where \(\mathbf{W}^{\text{(con-map)}}\in\mathbb{R}^{3\times D_{u}}\) is the learnable weight matrix of the MLP, \(\mathbf{C}_{i}\in\mathbb{R}^{m_{v}\times D_{u}}\) and \(\mathbf{u}_{i}^{\text{(con)}}\in\mathbb{R}^{D_{u}}\). Note that \(\mathbf{u}_{i}^{\text{(con)}}\) is variable and constraint permutation invariant. The complete feature matrix for the constraint nodes is constructed by concatenating all the feature vectors \(\mathbf{u}_{i}^{\text{(con)}}\), i.e., \(\mathbf{U}^{\text{(con)}}=[\mathbf{u}_{1}^{\text{(con)}},\dots,\mathbf{u}_{D_ {c}}^{\text{(con)}}]\in\mathbb{R}^{D_{c}\times D_{u}}\).

It is possible to differentiate between the equality and inequality constraints and treat them separately. In that case, a different set of triplets is constructed for each constraint and a different MLP can be used to map the triplets to the feature dimension.

The final feature matrix \(\mathbf{U}\) is:

\[\mathbf{U}=\begin{bmatrix}\mathbf{U}^{\text{(var)}}\\ \mathbf{U}^{\text{(con)}}\end{bmatrix}\in\mathbb{R}^{|\mathcal{V}|\times D_{u}}. \tag{10}\]

## Appendix D MLP for Continuous Variables

Given an assignment of integer variables, we can re-write the MILP (1) as the following parametric linear program:

\[V(\mathbf{z}^{\text{(b)}})=\underset{\mathbf{z}}{\text{minimize}} (\mathbf{c}^{\text{(b)}})^{T}\mathbf{z}^{\text{(b)}}+(\mathbf{c}^{ \text{(c)}})^{T}\mathbf{z}^{\text{(c)}}\] (11) subject to \[\mathbf{A}^{\text{(b)}}\mathbf{z}^{\text{(c)}}\leq\mathbf{b}- \mathbf{A}^{\text{(b)}}\mathbf{z}^{\text{(b)}}.\]

where \(\mathbf{z}=[\mathbf{z}^{\text{(b)}};\mathbf{z}^{\text{(c)}}]\), \(\mathbf{A}=[\mathbf{A}^{\text{(b)}},\mathbf{A}^{\text{(c)}}]\), and \(\mathbf{c}=[\mathbf{c}^{\text{(b)}},\mathbf{c}^{\text{(c)}}]\). The above function \(V:\mathbb{R}^{D_{x}}\rightarrow\mathbb{R}\) is a piecewise convex function [26, 13], and therefore it can be approximated by an MLP with ReLU activation units. This fact motivated us to use an MLP to map the embedding from the GCN to the continuous variables. Notice that also a convex layer [3] can be used to approximate optimization programs. However, training convex layers is expensive as it requires to differentiate through the optimization program and in our testing it significantly increased the time needed for training. For this reason, we opted to use an MLP instead.

## Appendix E On the Bernoulli and Beta Distributions

The goal of our model is to learn the value of each binary variable as well as measure how reliable is each prediction. A natural approach to assess how reliable or risky is a decision is via the variance. Recall that a Bernoulli distribution,
defined as \(\text{Ber}(\pi)\), has a mean \(\pi\) and variance \(\pi(1-\pi)\). Thus, for a given mean, the variance is completely defined and fixed, i.e., we cannot have different variance values for a given mean value. This essentially renders the variance uninformative (given the mean) for the purpose of assessing how reliable is a binary variable prediction. On the other hand, in the Beta distribution a fixed mean does not correspond to a fixed variance value and therefore can be used (along with the mean) as a measure to select trustworthy variables. For example, when selecting the \(\rho\%\) most reliable variable predictions, one might prefer to include a variable with mean 0.2 and small variance (and fix it to the integer value 0) than a variable with mean 0.1 and very large variance, which indicates that the model has low confidence on that particular prediction although its mean is closer to 0.

## Appendix F Training

### Clenshaw-Curtis Quadrature

The Clenshaw-Curtis quadrature [15] is a method for numerical integration that is based on the expansion of the integrand in terms of Chebyshev polynomials, where the function \(f(x)\) to be integrated over the fixed interval \([-1,1]\) is evaluated at the \(K\) roots of a Chebyshev polynomial. Then, the integral can be approximated as:

\[\int_{-1}^{1}f(x)dx\approx\mathbf{w}^{\top}\mathbf{y}, \tag{12}\]

where

\[\mathbf{w}=\mathbf{D}^{\top}\mathbf{d}, \tag{13}\]

are the quadrature weights, with \(\mathbf{d}\in\mathbb{R}^{K/2+1}\) and \(\mathbf{D}\in\mathbb{R}^{(K/2+1)\times(K/2+1)}\). In particular, the \(k\)-th element of \(\mathbf{d}\) (with zero-based indexing) is given by

\[d_{k}=\begin{cases}1,&k=0,\\ 2/(1-(2k)^{2}),&k=1,\ldots,K/2-1,\\ 1/(1-K^{2}),&k=K/2,\end{cases} \tag{14}\]

while the \((m,k)\)-th element of \(\mathbf{D}\) is defined as

\[D_{mk}=\frac{2}{K}\cos\left(\frac{mk\pi}{K/2}\right)\times\begin{cases}1/2,&k= 0,K/2,\\ 1,&\text{otherwise}.\end{cases} \tag{15}\]

Notice that both \(\mathbf{D}\) and \(\mathbf{d}\) are independent of the function \(f(x)\) and therefore can be precomputed in \(O(K\log K)\). The information of the function is encoded in the vector \(\mathbf{y}\in\mathbb{R}^{K/2+1}\), where its \(k\)-th element, with \(k=0,\ldots,K/2\), can be computed as

\[y_{k}=f(\cos(k\pi/K))+f(-\cos(k\pi/K)). \tag{16}\]

Finally, if the function needs to be integrated in the \([0,1]\) interval, we can apply the simple change of variable \(x^{\prime}=2x-1\to x=(x^{\prime}+1)/2\) and the new integral becomes

\[\int_{-1}^{1}\frac{1}{2}f((x^{\prime}+1)/2)dx^{\prime}, \tag{17}\]

where we have also accounted for the Jacobian factor which is equal to \(1/2\). The integral in (17) has the right limits and the Clenshaw-Curtis quadrature can be applied accordingly.

### Regularization

Ideally, we would like the model to produce Beta distributions with mean close to 0.5 and high variance if it is not confident. We propose the following regularization term to reinforce this behaviour:

\[r(\alpha,\beta,z^{\star})= \int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)|z^{\star}-\pi|d\pi \cdot\text{KL}\left(\text{U}(0,1)||\text{Beta}(\pi,\alpha,\beta)\right)\] \[= \left((1-z^{\star})\int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)\pi d \pi+z^{\star}\int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)(1-\pi)d\pi\right)\] \[\times\int_{0}^{1}\log\left(\frac{1}{\text{Beta}(\pi;\alpha, \beta)}\right)d\pi \tag{18}\] \[= \left((1-z^{\star})\frac{\alpha}{\alpha+\beta}+z^{\star}\frac{ \beta}{\alpha+\beta}\right)\left(-\int_{0}^{1}\log\left(\frac{\pi^{\alpha-1}( 1-\pi)^{\beta-1}}{B(\alpha,\beta)}\right)d\pi\right)\] \[= \left(\frac{(1-z^{\star})\alpha+z^{\star}\beta}{\alpha+\beta} \right)\left(\alpha-1+\beta-1+\log\left(B(\alpha,\beta)\right)\right),\]
where \(\text{Uni}(0,1)\) is the uniform distribution in the unit interval and \(B(\alpha,\beta)\) the Beta function. The goal is to minimize the KL-divergence between the learned Beta distribution and the Uniform only when a prediction is wrong. The term \(\int_{0}^{1}\text{Beta}(\pi;\alpha,\beta)|z^{\star}-\pi|d\pi\) evaluates how far is the distribution from the true label \(z^{\star}\): if Beta _is close_ to the label and has small variance (confident and correct) then this term is very small and the regularization becomes negligible. If the distribution _is far_ from the label, the KL divergence regularizes Beta towards Uniform (pushes its mean towards 0.5) and increases its variance.

The overall supervised loss is given by:

\[\ell_{\text{sup-reg}}=\ell_{\text{sup}}(\mathbf{z}^{\star},\psi=\{\boldsymbol{ \alpha}_{s,t},\boldsymbol{\beta}_{s,t}\})+\lambda_{\text{reg}}r(\mathbf{z}^{ \star},\psi=\{\boldsymbol{\alpha}_{s,t},\boldsymbol{\beta}_{s,t}\}), \tag{19}\]

where \(\lambda_{\text{reg}}\geq 0\) is a regularization parameter.

### Weighted Loss

It is common in classification tasks to include weights in the loss when the classes are imbalanced or to trade off precision and recall.

Our problem can be interpreted as a classification task with \(2^{D^{\text{(b)}}_{z}}\) classes, one for each possible binary sequence of length \(D^{\text{(b)}}_{z}\), and each instance then belongs to a specific class. With this interpretation the number of classes explodes very fast and no meaningful technique can be applied since the samples required to have a decent representation of each class in unattainable (e.g., for \(D^{\text{(b)}}_{z}=100\) and assuming an average of only 10 samples per class we would require \(2^{100}\cdot 10\approx 10^{7}\) samples).

An interpretation that scales is to consider our problem as \(D^{\text{(b)}}_{z}\) classification tasks (one fore each variable) with 2 classes each (since we consider binary variables). With this approach, each instance is a combination of \(D^{\text{(b)}}_{z}\) classes. Now, each classification task has as many samples as the number of instances, which are split in a different way among the 2 classes.

Note that the above interpretations do not change the mathematical formulation of the problem or the loss but the way we think of the class weighting. Following the second approach, we can calculate the class representation percentage \(r_{j}\) of the \(j\)-th variable (i.e., each classification task) from the available instances (labels). For example, if in \(80\%\) of the labels \(z_{j}=0\) and in \(20\%\) of the labels \(z_{j}=1\), then \(r_{j}=0.2\). Then, the (supervised) loss in its generic form would be:

\[\ell(\mathbf{z}^{\star},\phi)=-\sum_{s=1}^{S}\sum_{t=1}^{T}\sum_{j=1}^{D_{s}} \frac{\log(\cdot)}{r_{j}^{z_{s,t,j}}(1-r_{j})^{1-z_{s,t,j}^{\star}}}. \tag{20}\]

The use (or not) of these weights for the computation of the supervised loss is a hyperparameter that is optimized for each dataset.

## Appendix G Datasets

We introduce the six types of problems which we employed in the evaluation of MIPnet. The results for the problems from Tables 5 and 7 can be found in Section 4, and the results from problems from Tables 6 and 8 are discussed in Section I.1 of the Appendix. These problems were selected in line with the evaluation studies in prior works, see [52; 19; 11] and references therein; they appear regularly in different application domains and, furthermore, constitute building blocks for a very wide range of optimization problems. In the sequel we describe the mathematical formulation of each problem, denoted as \(\mathbb{P}_{p}^{t}\) with \(t\in\mathcal{T}\) and \(p=1,\ldots,6\), to indicate that each formulation corresponds to a particular instance of the problem with time-varying parameters.78 We further identify the problems' complexity and class (type of problem), and we map them to real-world applications. Finally, we provide details about the datasets and instances we used in training, testing and validation - see also the summarized view in Tables 5-8.

Footnote 7: We have ignored the time series index \(s\) to avoid notational cluttering.

Footnote 8: For the problem formulations we have used \(t\) as a superscript in all the time-varying parameters (\(t\) is fixed per instance and it just indicates that a specific parameter changes across instances), while all the subscripts are instance-dependent indices.

### Multicommodity Network Design and Routing

We consider the joint unsplittable-routing and network design problem where the goal is to minimize the aggregate routing cost by deciding how much capacity to purchase for each link and which path to select for each commodity.
Formally, the network is modeled by a graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) and serves a set of \(\mathcal{K}\) commodities with time varying demands \(\boldsymbol{d}^{t}=\{d_{k}^{t}\}_{k\in\mathcal{K}}\), and predetermined sources and destinations. Specifically, our goal is to solve the following _binary problem:_

\[\mathbb{P}_{1}^{t}:\quad\underset{\boldsymbol{p},\boldsymbol{q}}{ \text{minimize}} \sum_{k\in\mathcal{K}}\sum_{l\in\mathcal{P}_{k}}v_{k,l}p_{k,l}d_{k}^ {t}+\sum_{e\in\mathcal{E}}\sum_{i\in\mathcal{I}}u_{i}q_{e,i}\] subject to \[\sum_{l\in\mathcal{P}_{k}}p_{k,l}=1, \forall k\in\mathcal{K},\] \[\sum_{k\in\mathcal{K}}\sum_{l\in\mathcal{P}_{k}}a_{e,k,l}p_{k,l}d _{k}^{t}\leq b_{e}+\sum_{i\in\mathcal{I}}c_{i}q_{e,i}, \forall e\in\mathcal{E},\] \[\boldsymbol{p}\in\{0,1\}^{|\mathcal{K}|},\ \boldsymbol{q}\in\{0,1\}^{| \mathcal{E}||\mathcal{I}|}.\]

Problem \(\mathbb{P}_{1}^{t}\) has the following binary _variables_:

* \(p_{k,l}\in\{0,1\}\): selection of path \(l\) for routing commodity \(k\in\mathcal{K}\),
* \(q_{e,i}\in\{0,1\}\): installment of \(i\)-type capacity at edge \(e\in\mathcal{E}\),

and the following _parameters_:

* \(\mathcal{P}_{k}\): set of eligible paths for commodity \(k\in\mathcal{K}\),
* \(v_{k,l}\in\mathbb{R}_{+}\): routing cost (per unit of traffic) of path \(l\in\mathcal{P}_{k}\),
* \(a_{e,l}^{k}\in\{0,1\}\): parameter indicating that edge \(e\) is contained in path \(l\in\mathcal{P}_{k}\),
* \(b_{e}\in\mathbb{R}_{+}\): initial capacity installed at edge \(e\in\mathcal{E}\),

\begin{table}
\begin{tabular}{c|c c c} \hline \hline  & routing-sm & revenue-max-sm \\ \hline \hline \# Integer variables & 1044 & 10000 \\ \hline \# Continuous variables & 0 & 0 \\ \hline \# Equality constraints & 390 & 0 \\ \hline \# Inequality constraints & 36 & 10 \\ \hline \hline Data & real & synthetic \\ \hline Fraction of non-zeros in \(\mathbf{A}\) & 0.0333 & 0.501 \\ \hline Fraction of non-zeros in \(\mathbf{b}\) & 1.0 & 1.0 \\ \hline Fraction of non-zeros in \(\mathbf{c}\) & 0.9169 & 1.0 \\ \hline Fraction of non-zeros in \(\mathbf{z}^{*}\) & 0.4298 & 0.0166 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Additional datasets problem sizes.

\begin{table}
\begin{tabular}{c|c c c c c c} \hline \hline  & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline \# Integer variables & 2812 & 1989 & 144 & 1000 & 30000 & 39942 \\ \hline \# Continuous variables & 0 & 0 & 500 & 0 & 0 & 0 \\ \hline \# Equality constraints & 450 & 50 & 36 & 1 & 0 & 0 \\ \hline \# Inequality constraints & 36 & 39 & 4082 & 1010 & 20 & 1 \\ \hline Data & real & mixed & synthetic & synthetic & synthetic & mixed \\ \hline Fraction of non-zeros in \(\mathbf{A}\) & 0.0374 & 0.0226 & 0.1414 & 0.6723 & 0.50 & 1.0 \\ \hline Fraction of non-zeros in \(\mathbf{b}\) & 1.0 & 0.5 & 0.8333 & 0.9876 & 1.0 & 1.0 \\ \hline Fraction of non-zeros in \(\mathbf{c}\) & 0.579 & 0.9804 & 0.9996 & 1.0 & 1.0 & 0.0187 \\ \hline Fraction of non-zeros in \(\mathbf{z}^{*}\) & 0.0001 & 0.0278 & 0.0833 & 0.7375 & 0.0044 & 0.0006 \\ \hline \end{tabular}
\end{table}
Table 5: Datasets problem sizes.
* \(u_{i}\in\mathbb{R}_{+}\): cost of buying \(i\)-th capacity installment,
* \(c_{i}\in\mathbb{R}_{+}\): capacity of \(i\)-th installment,
* \(d_{k}^{t}\in\mathbb{R}_{+}\): demand of commodity \(k\) at time \(t\).

**Complexity**. \(\mathbb{P}_{1}^{t}\) is known to be NP-hard even for a single commodity, i.e., \(|\mathcal{K}|=1\), cf. [55].

**Datasets**. The problem parameters, including the graph topology, were taken from the SNDlib database [53] that contains real-world traffic matrices and network topologies from different communication (and other) networks (e.g., from backbone ISP networks).9 In particular, for this problem we used the dataset _geant_ dataset, the size of which can be seen in Table 5. The capacity and installment costs are fixed over time, while the demand vector \(\mathbf{d}^{t}\) is being updated every \(15\) minutes over a period of \(4\) months. We note that we have scaled the demand \(\mathbf{d}^{t}\) and edge capacity by a factor of \(100\) and \(1/40\), respectively, to increase the complexity of the problem.10 The number of time series we used, and the number of instances for each time series can be found in Table 7, where the instances differ in vector \(\mathbf{d}^{t}\).

Footnote 9: SNDlib is a rich and well-known database which is regularly used in the evaluation of routing algorithms in communication networks.

Footnote 10: This scaling renders capacity purchase necessary in order to fulfill the demands.

### Facility Location

In the facility location problem, given a set \(\mathcal{I}\) of _facility_ locations and a set \(\mathcal{J}\) of _clients_, the goal is to decide which facilities to open and how to assign clients to the opened facilities. The objective is to minimize the total opening and assignment cost. More formally, we solve the following _binary problem_:

\[\mathbb{P}_{2}^{t}:\quad\text{minimize} \sum_{i\in\mathcal{I}}\sum_{j\in\mathcal{J}}c_{j,t}d_{j}^{t}z_{j, i}+\sum_{i\in\mathcal{I}}f_{i}x_{i}\] subject to \[\sum_{i\in\mathcal{I}}z_{j,i}=1, \forall j\in\mathcal{J},\] \[\sum_{j\in\mathcal{J}}z_{j,i}\leq 2|\mathcal{J}|x_{i}, \forall i\in\mathcal{I},\] \[\mathbf{z}\in\{0,1\}^{|\mathcal{I}||\mathcal{J}||\mathcal{I}|},\]

where the vector of binary variables \(\mathbf{z}=[z_{11},\dots,z_{|\mathcal{I}||\mathcal{J}|},x_{1},\dots,x_{| \mathcal{I}|}]\). In the above problem, we introduced the following _binary variables_:

* \(z_{j,i}\in\{0,1\}\): assignment of client \(j\in\mathcal{J}\) to facility \(i\in\mathcal{I}\),
* \(x_{i}\in\{0,1\}\): deployment of facility \(i\in\mathcal{I}\),

\begin{table}
\begin{tabular}{c|c c c c c c} \hline \hline  & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline \# of time series for training & 46 & 80 & 80 & 80 & 80 & 80 \\ \hline \# of time series for testing & 5 & 20 & 20 & 20 & 20 & 20 \\ \hline \# of time series for validation & 5 & 20 & 20 & 20 & 20 & 20 \\ \hline \# of instances per time series & 96 & 700 & 400 & 100 & 100 & 150 \\ \hline \hline \end{tabular}
\end{table}
Table 7: Datasets time series dimensions.

\begin{table}
\begin{tabular}{c|c c} \hline \hline  & routing-sm & revenue-max-sm \\ \hline \hline \# of time series for training & 96 & 80 \\ \hline \# of time series for testing & 18 & 20 \\ \hline \# of time series for validation & 10 & 20 \\ \hline \# of instances per time series & 96 & 100 \\ \hline \hline \end{tabular}
\end{table}
Table 8: Additional datasets time series dimensions.
and the following _parameters_:

* \(f_{i}\in\mathbb{R}_{+}\): cost for opening facility \(i\in\mathcal{I}\),
* \(c_{j,i}\in\mathbb{R}_{+}\): cost for associating client \(j\in\mathcal{J}\) to facility \(i\in\mathcal{I}\),
* \(d^{t}_{j}\in\mathbb{R}_{+}\): demand of client \(j\in\mathcal{J}\) at time \(t\).

**Complexity**. \(\mathbb{P}^{t}_{2}\) is the metric uncapacitated facility location problem, which is NP-hard and inapproximable below \(1.463\) ratio, with the currently-known algorithm achieving a solution within \(1.488\) ratio, see [67].

**Datasets**. The graph topology is borrowed from the _germany50_ dataset [53], the cost parameters \(f_{i}\) are randomly generated from a uniform distribution for each problem instance and the time-varying demand \(\boldsymbol{d}^{t}=\{d^{t}_{j}\}_{j\in\mathcal{J}}\) is defined as:

\[\boldsymbol{d}^{t+1}=\max\left\{0,\boldsymbol{A}\boldsymbol{d}^{t}+a_{1}\sin( t/\mathtt{period}_{1})+a_{2}\sin(t/\mathtt{period}_{2})+\mathbf{w}\right\},\]

for a Hurwitz matrix \(\mathbf{A}\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}\) constructed by randomly sampling the eigenvalues \(\lambda_{i}\sim\mathcal{U}(0.98,0.999)\) and the associated orthonormal eigenvectors \(\mathbf{v}_{1}\)11 for \(i\in\{1,\ldots,|\mathcal{J}|\}\). The random variables \(a_{1}\sim\mathcal{U}(1,5)\) and \(a_{2}\sim\mathcal{U}(2,10)\), and the random vector \(\mathbf{w}\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma}_{1})\), where the vector of zeros \(\mathbf{0}\in\mathbb{R}^{|\mathcal{J}|}\) and the covariance matrix \(\mathbf{\Sigma}_{1}\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}\) is a diagonal matrix of ones. In the above definition the \(\max\) operator is applied component-wise and guarantees nonnegative demand. Finally, the periodicity of the sine functions are \(\mathtt{period}_{1}=20\) and \(\mathtt{period}_{2}=70\).

Footnote 11: The eigenvectors are computed randomly sampling from a standard normal and then using the Gram–Schmidt procedure to compute a set of orthonormalized vectors.

### Travelling Salesman

We consider the Travelling Salesman Problem (TSP) under the Dantzig-Fulkerson-Johnson formulation [64]. In detail, we have a set \(\mathcal{N}\) of \(N\) cities, which we need to visit exactly once and return to the starting node. The cities are connected with a graph \(G=(\mathcal{N},\mathcal{E})\) which induces routing (or distance) costs \(c_{i,j},\forall(i,j)\in\mathcal{E}\). Our goal is to minimize the aggregate (travelling) cost of the route. In detail, the formulated _binary problem_ is:

\[\mathbb{P}^{t}_{3}:\quad\underset{\mathbf{z}}{\text{minimize}} \sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{N}}c^{t}_{i,j}z_{i,j}\] subject to \[\sum_{i\in\mathcal{N}}z_{i,j}=1, \forall j\in\mathcal{N},\] \[\sum_{j\in\mathcal{J}}z_{i,j}=1, \forall i\in\mathcal{N},\] \[\sum_{i\in\mathcal{S}}\sum_{j\in\mathcal{S}}z_{i,j}\leq|\mathcal{ S}|-1, \forall\mathcal{S}\subseteq\{1,\ldots,N\},\;|\mathcal{S}|\geq 2,\] \[\mathbf{z}\in\{0,1\}^{N^{2}}.\]

The problem includes the _binary variables_:

* \(z_{i,j}\in\{0,1\}\): inclusion of link \((i,j)\in\mathcal{E}\) or not (\(x_{ij}=0\)) in the route,

and _parameters_:

* \(c^{t}_{i,j}\in\mathbb{R}_{+}\): cost for traversing link \((i,j)\) at time \(t\).

**Complexity**. \(\mathbb{P}^{t}_{3}\) is NP-hard to approximate within any polynomial factor for general distance functions, while improved approximation ratios are available for certain restricted cases (metric spaces, etc.). We refer the reader to [63] for an up-to-date discussion on the complexity of TSP.

**Datasets**. We created instances with \(N=12\) nodes. The location of each node is randomly generated using a uniform distribution and we assumed that the graph is fully connected. The cost vector \(\boldsymbol{c}^{t}\in\mathbb{R}^{n}\) is given by the distance between each city pair and it changes at each time step as follows:

\[\boldsymbol{c}^{t+1}=\max\{0,\boldsymbol{c}^{t}+\mathbf{w}\},\]

where \(\mathbf{w}\sim\text{U}(d_{min},d_{max})^{n}\) and the \(\max\) operator is applied component-wise; where \(d_{min},d_{max}\) are positive constants.
### Revenue Maximization

We consider the problem of shipping commodities from a set \(\mathcal{N}\) of \(N=|\mathcal{N}|\) source nodes (each commodity corresponds to one node) towards their intended destinations over predetermined (overlapping) paths. Our goal is to maximize the revenue from delivering as many commodities as possible, while satisfying the time-varying capacity of each edge \(i\in\mathcal{I}\) that lies along the path of each commodity \(n\in\mathcal{N}\), i.e., \(i\geq n\). Formally, we consider the following _binary_ problem:

\[\mathbb{P}_{4t}:\quad\underset{\mathbf{z}}{\text{maximize}} \sum_{n\in\mathcal{N}}c_{n}^{t}z_{n}\] subject to \[\sum_{n\in\mathcal{N}}a_{i,n}z_{n}\leq b_{i}^{t},\quad\forall i\in \mathcal{I},\] \[\mathbf{z}\in\{0,1\}^{N}.\]

The problem includes the _binary variables_:

* \(z_{n}\in\{0,1\}\): transporting commodity \(n\) to its destination via a predetermined path, or not,

and _parameters_:

* \(c_{n}^{t}\in\mathbb{R}_{+}\): revenue when delivering commodity \(n\), at time \(t\),
* \(b_{i}^{t}\in\mathbb{R}_{+}\): transportation capacity of link \(i\in\mathcal{I}\), at time \(t\),
* \(a_{i,n}\in\{0,1\}\): indicating if link \(i\) lies in the predetermined path of commodity \(n\), i.e., \(i\geq n\).

**Complexity**. \(\mathbb{P}_{4}^{t}\) is NP-complete as it falls in the category of multi-dimensional Knapsack problems, see [66].

**Datasets.** We randomly generated from a uniform distribution the scalars \(a_{i,n}\) representing the resources needed to ship commodity \(n\) from the \(i\)-th facility. Furthermore, for each time series the time-varying capacity and revenue vectors are defined as follows:

\[\mathbf{c^{t+1}} =\mathbf{c^{t}}+a_{1}\sin(t/\mathtt{period}_{1})+a_{2}\sin(t/ \mathtt{period}_{2})+\mathbf{w}_{c},\] \[\mathbf{b^{t+1}} =\mathbf{b^{t}}+\mathbf{w}_{b},\]

where \(\mathbf{c^{t}}=[c_{1}^{t},\ldots,c_{N}^{t}]\), \(\mathbf{b_{t}}=[b_{1}^{t},\ldots,b_{|\mathcal{I}|}^{t}]\), the scalars \(a_{1}\), \(a_{2}\) are randomly generated for each time series, and the random vectors \(\mathbf{w}_{c}\in\mathbb{R}^{|\mathcal{I}||\mathcal{N}|}\) and \(\mathbf{w}_{b}\in\mathbb{R}^{|\mathcal{I}|}\). For further details about the number of time series and problem instances per time series please refer to Table 7.

### Energy Grid

We consider a microgrid energy-sharing problem, cf. [65], where a set \(\mathcal{F}\) of \(F=|\mathcal{F}|\) energy producers-consumers (or, _prosumers_) coordinate their energy consumption plan over several time periods. The prosumers have at their disposal a hierarchical energy storage system, where the primary batteries are already in place and the secondary batteries are deployed upon demand by paying additional cost. The goal is to optimize the overall microgrid operation, by solving

Figure 4: Microgrid energy storage optimization for a community of \(\mathcal{F}\) prosumers and a two-tier battery system.
the _mixed binary problem_:

\[\mathbb{P}_{5}^{t}:\underset{\mathbf{z}}{\text{maximize}} \sum_{f\in\mathcal{F}}c_{f}^{t}z_{f}^{(c)}-\sum_{i\in\mathcal{I}}p_{ i}^{t}z_{i}^{(b)}\] subject to \[\sum_{f\in\mathcal{F}}a_{nf}z_{f}^{(c)}\leq b_{n}^{t}+\sum_{i\in \mathcal{I}}d_{n,i}z_{i}^{(b)},\quad\forall n\in\mathcal{N},\] \[\sum_{f\in\mathcal{F}}z_{f}^{(c)}=1,\] \[\mathbf{z}^{(c)}\in\mathbb{R}^{F},\,\mathbf{z}^{(b)}\in\{0,1\}^{F},\]

where the last constraint enforces the prosumers \(\mathcal{F}\) to keep at least 1 (normalized) unit of energy for serving the needs of their community.

The main _variables_ of \(\mathbb{P}_{5}^{t}\) are:

* \(z_{f}^{(c)}\in\mathbb{R}\): continuous variable deciding how much energy prosumer \(f\in\mathcal{F}\) will store (if positive) or purchase (if negative),12 Footnote 12: Negative transfers release capacity that can be used for stored-energy transfers.
* \(z_{i}^{(b)}\in\{0,1\}\): decides whether to deploy or not (\(z_{i}=0\)) the secondary battery \(i\in\mathcal{I}\),

and the main _parameters_:

* \(a_{n,f}\in\mathbb{R}_{+}\): energy transfer loss coefficient when prosumer \(f\in\mathcal{F}\) transfers energy to battery \(n\in\mathcal{N}\). Larger values indicate larger losses,
* \(c_{f}^{t}\in\mathbb{R}_{+}\): selling or purchase (when \(z_{f}^{(c)}<\)0) price for the energy of prosumer \(f\in\mathcal{F}\), at time \(t\),
* \(b_{n}^{t}\in\mathbb{R}_{+}\): available capacity at primary battery \(n\in\mathcal{N}\) at time \(t\),
* \(d_{n,i}\in\mathbb{R}_{+}\): energy loss coefficient when transferring energy from primary battery \(n\in\mathcal{N}\) to secondary battery \(i\in\mathcal{I}\).

**Complexity**. \(\mathbb{P}_{5}^{t}\) is NP-hard as it generalizes the facility location problem [65].

**Datasets**. We created instances using a non-stationary formula for the time-varying vectors that exhibits temporal properties as follows:

\[\mathbf{c}^{t+1} =\mathbf{c}^{t}+a\sin(t/\texttt{period}_{1})+\mathbf{w}_{c},\] \[\mathbf{b}^{t+1} =\mathbf{b}^{t}+\mathbf{w}_{\mathbf{b}},\]

where the scalar \(a\) is randomly generated for each time series, and the random vectors \(\mathbf{w}_{c}\in\mathbb{R}^{|\mathcal{I}||\mathcal{P}|}\) and \(\mathbf{w}_{b}\in\mathbb{R}^{|\mathcal{I}|}\).

### Caching

Our final problem is a standard data caching problem where we wish to store the most popular files in a cache of limited capacity, see [29]. In particular, we consider a content library of \(\mathcal{N}=\{1,2,\ldots,N\}\) files and a single cache with \(C\) bytes capacity. Each file \(n\in\mathcal{N}\) has popularity \(p_{n}^{t}\) at time \(t\), which captures the expected requests for file \(n\) at that time, and size of \(q_{n}\) bytes. The goal is to store, at each time, those files that will be requested by the larger number of content viewers, i.e., to maximize the cache hits. This can be formalized with the following _binary problem_:

\[\mathbb{P}_{6}^{t}:\underset{\mathbf{z}}{\text{maximize}} \sum_{n=1}^{N}p_{n}^{t}x_{n}\] subject to \[\sum_{n=1}^{N}x_{n}q_{n}\leq C\] \[x_{n}\in\{0,1\},\,\,\,\forall n\in\mathcal{N}.\]

The _variables_ of \(\mathbb{P}_{6}^{t}\) are:

* \(x_{n}\in\{0,1\}\): binary variables that select a file to be cached or not,
and the main _parameters_:

* \(p_{n}^{t}\in\mathbb{R}_{+}\): normalized file popularity parameters at time \(t\),
* \(q_{n}\in\mathbb{R}_{+}\): size of file \(n\in\mathcal{N}\).

**Complexity**. \(\mathbb{P}_{6}^{t}\) is a standard Knapsack problem which is NP-Complete but can be solved in pseudo-polynomial time through dynamic programming.

**Datasets**. The file sizes were created randomly from a uniform distribution \(\mathcal{U}(1,10)\), while for the file popularity we used the standard dataset for request traces from Movielens [24], which we sliced in order to create the different time series.

## Appendix H Experimental Details

### Benchmarks

In our experiments we use as benchmark the N-div model [52] as it is the most related approach to our method MIPnet.

N-div uses a GCN based on the bipartite graph representation of the MILP, similar to our approach, while the output of the GCN is mapped to Bernoulli parameters using an MLP. The model is trained by minimizing the NLL (cross-entropy) modified based on the SelectiveNet approach [32]. The idea is to learn two sets of binary variables: one set indicates if a variable is going to be selected and the other its binary value (to be used if selected). Therefore, the overall loss can be seen as a weighted cross-entropy. Further, this loss assumes a predefined target selection percentage and therefore a different model is required for different variable selection percentages.

We implemented N-div by keeping the common structure of our model and changing the final layers and the loss. Further, we used the same set of features (see Appendix C) as in our method for a fair comparison. Note that in [52] it is not clearly mentioned what features are used for the N-div method. The authors only mention a set of solver-based features that can be used in their neural-branching method but are not applicable for N-div since no such features are available during inference.

For the method from [11], we use the implementation available online at github.com/bstellato/mlopt. This strategy learns a mapping from problem parameters to feasible variable assignments. First, the input data is categorized into strategies, i.e., feasible variables assignment. Then, the mapping is learned as a multi-class classification problem. Notice that this approach does not generalize when a class is not in the training dataset, i.e., when at test time all feasible assignments seen during training are note feasible for the new problem instance. In order to scale to large datasets, we modified their implementation so that it accepts sparse training data. Moreover, we made sure that their method receives the same train/validation/test split as our method.

### Hyperparameters

For all datasets we used the Adam optimizer for the gradient updates with \(10^{-5}\) weight-decay and clipped the gradient norm to 10. The learning rate was warmed-up linearly from \(1\times 10^{-4}\) to \(1\times 10^{-2}\) for the first \(500\) steps after which a cosine decay follows for the remaining time steps with a decay rate of 0.99.

For each dataset we performed a hyperparameter optimization on the model and loss parameters. The batch size and the number of training steps were tuned per dataset with values in \(\{8,12,16,32\}\) and \(\{8000,15000,20000\}\), respectively. All the loss hyperparameters were linearly warmed-up to an initial value and then linearly scaled up to

\begin{table}
\begin{tabular}{c|c c c} \hline \hline Schedule & \multicolumn{3}{c}{Parameter} \\ \hline  & \(\lambda\) & \(\lambda_{\text{avg}}\) & \(\lambda_{\text{c}}\) \\ \hline \hline warm-up steps & \{500, 1000, 1500, 2000\} & \{250,500\} & \{1000,1500\} \\ warm-up initial value & \(\{0.01,0.1\}\) & \(\{0.01,0.1\}\) & \(0.1\) \\ warm-up final value & \(\{0.1,1\}\) & \(\{0.1,1\}\) & \(\{1,10\}\) \\ final value & \(\{1,10,50\}\) & \(\{0.1,1,10,100\}\) & \(\{10,100\}\) \\ \hline \hline \end{tabular}
\end{table}
Table 9: Hyperparameters of the loss components.
a final target value in the final training step. The loss scheduling hyperparameter ranges are provided in Table 9. The model hyperparameter ranges are provided in Table 10.13

Footnote 13: Note that the number of units in GCN increases due to skip connections.

### Computational Environment

For running our experiments we used p3.8xlarge AWS EC2 instances with 4 Tesla V100 GPUs, 36 CPUs, and 244 GB of memory, and m5.8xlarge instances with 32 CPUs and 128 GB of memory.

## Appendix I Additional Results

### Additional Datasets

In Table 11 we provide results for accuracy, infeasibility, and optimality gap for the two additional datasets routing-sm and revenue-max-sm, which are smaller versions of routing and revenue-max-sm, respectively. Table 12 shows the accuracy, infeasibility, and optimality gap on these two additional datasets for the method from [11]. We notice that our method performs better than N-div and the method from [11] in almost all additional experiments. Notice that the method from [11] is able to find a feasible solution for routing-sm in more cases than MIPnet and N-div. However, this method is always infeasible in revenue-max, MIPnet and N-div always find a feasible solution.

### Unsupervised Loss

As we already showed in Section 4.1, training MIPnet only with supervised loss (MIPnet-sup) leads to decreased performance for reasons we discussed throughout the paper. Here, we complete the picture by including the performance of MIPnet-unsup, i.e., MIPnet trained only with unsupervised loss. Tables 13-15 show the accuracy, infeasibility and optimality gap, respectively. It is clear that the performance of MIPnet-unsup compared to MIPnet is significantly degraded in most datasets, with the sole exception of the caching.

It is important to highlight some key observations based on these results. First, the accuracy of MIPnet-unsup in the energy-grid dataset is less than 50% for all values of \(\rho\). This does not lead to significant infeasibilities since apparently the constraints of the problem are not that hard to satisfy, however it causes a significant optimality gap. This is a case where labels make a big difference in the model training process and in fact the performance of MIPnet-sup is in par with MIPnet, i.e., most of the learning is done from the labels. On the other end, we observe that in the

\begin{table}
\begin{tabular}{c|c c c c} \hline \hline Parameter & \multicolumn{4}{c}{Block} \\ \hline  & MLP (features) & GCN & LSTM & MLP (binary map) & MLP (continuous map) \\ \hline \hline \# of layers & 1 & \(\{2,3\}\) & \(\{1,2,3\}\) & 1 & 1 \\ \# of units & \(\{8,16\}\) & \(\{8,16\}\) & \(\{16,32\}\) & \(2\times D_{i}^{(b)}\) & \(D_{i}^{(c)}\) \\ nonlinearity & \(\{\text{ReLU},\text{No}\}\) & \(\{\text{ReLU},\text{No}\}\) & No & ReLU & No \\ \hline \hline \end{tabular}
\end{table}
Table 10: Hyperparameters of the model blocks.

\begin{table}
\begin{tabular}{c c|c|c c|c c} \hline \hline  & & \multicolumn{2}{c|}{**Accuracy**} & \multicolumn{2}{c|}{**Infeasibility**} & \multicolumn{2}{c}{**Optimality gap**} \\ \hline \(\rho\) & Method & routing-sm & revenue-max-sm & routing-sm & revenue-max-sm & routing-sm & revenue-max-sm \\ \hline \hline \multirow{3}{*}{30\%} & MIPnet & **99.18 \(\pm\) 0.17** & **100.00 \(\pm\) 0.00** & **7.88 \(\pm\) 0.63** & **0.00 \(\pm\) 0.00** & **0.35 \(\pm\) 0.01** & **0.12 \(\pm\) 0.03** \\  & N-div & 95.67 \(\pm\) 0.48 & **100.00 \(\pm\) 0.00** & 77.91 \(\pm\) 9.73 & **0.00 \(\pm\) 0.00** & **0.39 \(\pm\) 0.04** & 0.16 \(\pm\) 0.05 \\ \hline \multirow{3}{*}{40\%} & MIPnet & **98.29 \(\pm\) 0.22** & **100.00 \(\pm\) 0.00** & **35.63 \(\pm\) 21.02** & **0.00 \(\pm\) 0.00** & 0.95 \(\pm\) 0.37 & **0.07 \(\pm\) 0.02** \\  & N-div & 78.33 \(\pm\) 25.14 & **100.00 \(\pm\) 0.00** & 96.27 \(\pm\) 4.32 & **0.00 \(\pm\) 0.00** & **0.27 \(\pm\) 0.02** & 0.11 \(\pm\) 0.07 \\ \hline \multirow{3}{*}{50\%} & MIPnet & **97.32 \(\pm\) 0.40** & **100.00 \(\pm\) 0.00** & **65.07 \(\pm\) 29.85** & **0.00 \(\pm\) 0.00** & **1.48 \(\pm\) 0.89** & 0.05 \(\pm\) 0.02 \\  & N-div & 89.46 \(\pm\) 0.93 & **100.00 \(\pm\) 0.00** & 100.00 \(\pm\) 0.00 & **0.00 \(\pm\) 0.00** & - & **0.03 \(\pm\) 0.03** \\ \hline \multirow{3}{*}{60\%} & MIPnet & **95.19 \(\pm\) 0.42** & **100.00 \(\pm\) 0.00** & **100.00 \(\pm\) 0.00** & **0.00 \(\pm\) 0.00** & - & 0.03 \(\pm\) 0.03 \\  & N-div & 92.41 \(\pm\) 1.26 & **100.00 \(\pm\) 0.00** & **100.00 \(\pm\) 0.00** & **0.00 \(\pm\) 0.00** & - & **0.02 \(\pm\) 0.01** \\ \hline \multirow{3}{*}{70\%} & MIPnet & **92.79 \(\pm\) 0.43** & **100.00 \(\pm\) 0.00** & **100.00 \(\pm\) 0.00** & **0.00 \(\pm\) 0.00** & - & **0.03 \(\pm\) 0.04** \\  & N-div & 89.52 \(\pm\) 25.55 & **100.00 \(\pm\) 0.00** & **100.00 \(\pm\) 0.00** & **0.00 \(\pm\) 0.00** & - & **0.03 \(\pm\) 0.01** \\ \hline \hline \end{tabular}
\end{table}
Table 11: Accuracy, infeasibility and optimality gap (mean \(\pm\) std) of MIPnet vs. N-div. Bold indicates best method.
caching dataset the model is able to learn the optimal values without labels. Arguably this is an easy problem to learn since all the models achieve optimal performance. In such cases the overhead of collecting (many) labels might not be necessary. All the other datasets fall in the case where both MIPnet-sup and MIPnet-unsup underperform compared to MIPnet and a combination of supervised and unsupervised loss gives the best performance.

### Missing Labels

Finally, we examine how missing labels affect the performance of a model trained with and without unsupervised loss. For this experiment we used the tsp dataset. Figure 5 illustrates the accuracy of the models for different percentages of missing labels. In the presence of missing labels the model that is trained with both supervised and unsupervised loss clearly outperforms the one without unsupervised loss in all cases.

## Appendix J Discussion of Related Work

Table 16 summarizes the key differences of MIPnet from the most-related models and approaches. In the evaluation section we have presented detailed comparisons of MIPnet with [52], [32] and [11], while here we focus on qualitative

\begin{table}
\begin{tabular}{c|c c c c c} \hline \hline \(\rho\) (\(\%\)) & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline
30\% & \(95.46\pm 1.47\) & \(99.02\pm 1.12\) & \(96.45\pm 2.69\) & \(49.38\pm 1.36\) & \(100.00\pm 0.00\) & \(100.00\pm 0.00\) \\ \hline
40\% & \(94.13\pm 1.16\) & \(98.90\pm 1.25\) & \(95.83\pm 2.90\) & \(48.16\pm 0.75\) & \(100.00\pm 0.00\) & \(100.00\pm 0.00\) \\ \hline
50\% & \(92.60\pm 1.13\) & \(98.85\pm 1.26\) & \(95.28\pm 3.04\) & \(47.61\pm 0.77\) & \(99.99\pm 0.00\) & \(100.00\pm 0.00\) \\ \hline
60\% & \(91.24\pm 1.01\) & \(98.76\pm 1.21\) & \(94.76\pm 3.03\) & \(47.59\pm 0.86\) & \(99.99\pm 0.00\) & \(100.00\pm 0.00\) \\ \hline
70\% & \(89.49\pm 0.76\) & \(98.50\pm 1.02\) & \(94.32\pm 2.84\) & \(47.78\pm 0.79\) & \(99.99\pm 0.00\) & \(100.00\pm 0.00\) \\ \hline \end{tabular}
\end{table}
Table 13: Accuracy (mean \(\pm\) std)\(\%\) of MIPnet-unsup.

\begin{table}
\begin{tabular}{c|c c c c c c} \hline \hline \(\rho\) (\(\%\)) & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline
30\% & \(100.00\pm 0.00\) & \(6.92\pm 0.56\) & \(33.43\pm 47.07\) & \(0.02\pm 0.03\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
40\% & \(100.00\pm 0.00\) & \(13.87\pm 10.09\) & \(41.94\pm 42.22\) & \(0.06\pm 0.00\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
50\% & \(100.00\pm 0.00\) & \(29.89\pm 25.36\) & \(54.68\pm 39.66\) & \(0.11\pm 0.06\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
60\% & \(100.00\pm 0.00\) & \(47.28\pm 24.40\) & \(69.16\pm 40.00\) & \(0.15\pm 0.03\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
70\% & \(100.00\pm 0.00\) & \(72.68\pm 15.37\) & \(84.04\pm 21.92\) & \(0.17\pm 0.08\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline \end{tabular}
\end{table}
Table 15: Optimality gap (mean \(\pm\) std)\(\%\) of MIPnet-unsup.

\begin{table}
\begin{tabular}{c|c c c c c c} \hline \hline \(\rho\) (\(\%\)) & routing & facility-loc & tsp & energy-grid & revenue-max & caching \\ \hline \hline
30\% & - & \(8.83\pm 12.33\) & \(1.55\pm 1.27\) & \(8.59\pm 0.04\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
40\% & - & \(12.44\pm 17.31\) & \(3.11\pm 2.78\) & \(12.62\pm 0.36\) & \(0.00\pm 0.00\) & \(0.00\pm 0.00\) \\ \hline
50\% & - & \(15.08\pm 20.06\) & \(5.95\pm 5.38\) & \(17.16\pm 0.63\) & \(0.66\pm 0.01\) & \(0.00\pm 0.00\) \\ \hline
60\% & - & \(18.97\pm 24.36\) & \(11.06\pm 9.82\) & \(22.00\pm 0.95\) & \(1.23\pm 0.01\) & \(0.00\pm 0.00\) \\ \hline
70\% & - & \(32.20\pm 27.12\) & \(19.73\pm 17.64\) & \(27.41\pm 1.11\) & \(1.57\pm 0.07\) & \(0.00\pm 0.00\) \\ \hline \end{tabular}
\end{table}
Table 12: Accuracy, infeasibility and optimality gap of [11].
differences. In detail, we see that MIPnet is the only model that uses semi-supervised learning. This is very important since obtaining labels requires, in most cases, solving large-scale NP-hard problems, while the unsupervised loss component can improve the performance of the model further, by allowing us to identify the effect (in terms of optimality gap and infeasibility) of mispredictions at a per-variable granularity. These benefits are evident in the experiments. We note that other works, such as [52] and [19], identify the importance of unsupervised learning as well, but this approach is essesentially explored only in [39].

In particular, [39] proposes an unsupervised approach specifically for graph problems where a GNN is used to learn a distribution over the graph nodes, representing a solution. In order to produce an integral solution, the authors derandomize the continuous values using sequential decoding. The method comes with theoretical guarantees, giving good and feasible solutions with high probability. However, this approach is not applicable to general MIPs, while even for the targeted graph problems it is not trivial to include general constraints. Besides, this work does not benefit from the availability of labels, which, in certain operational environments are available.

On the other hand, in [52] the authors provide only a preliminary evaluation with training datasets that might include noisy labels (obtained from suboptimal solutions); while [19] performs an initial variable assignment and then uses neighborhood search to increase the training data without solving the problem exactly - which works under the assumption of locality, as the authors stress. On the contrary, MIPnet uses a _semi-supervised_ learning approach which, apart from allowing us to expand the training dataset, enables tuning the impact of each variable on the feasibility and objective value of the problem.

Another distinguishing feature of MIPnet is that, unlike all prior works, it accounts for the temporal structure across the different instances, which, as shown in the experiments, indeed enhances the model's performance. This aspect is crucial as, more often than not, the different instances that practitioners solve exhibit a temporal structure, namely problem parameters such as the commodity volumes, user requests, transportation costs or electricity prices, follow some type of diurnal pattern.

We also note that the related works in Table 16 can be separated to those that use GCN in order to automate feature embedding, and to those that do not follow this approach, as e.g., [10; 11]. GCNs seem to enhance the performance and provide the means to account for dependencies across the variables and, indirectly, among the constraints as well. Finally, MIPnet is the only work that uses a Bayesian approach and a tunable confidence for performing the assignment of variables, unlike, e.g., the threshold-based variable selection rule of [19]. We note also that [52] offers the option to select different percentage \(\rho\) of the binary variables, but this decision needs to be made in advance so as to train the model accordingly.

\begin{table}
\begin{tabular}{c|c c c c} \hline \hline  & **Use GCN** & **Vars Select Method** & **Temporal** & **Learning** \\ \hline \hline
[52] & Yes & Bernoulli & No & Superv. \\ \hline
[11] & No & No & No & Superv. \\ \hline
[19] & Yes & Threshold-based & No & Superv. \\ \hline
[10] & No & No & No & Superv. \\ \hline
[31] & Yes & Threshold-based & No & Superv. \\ \hline
[39] & Yes & No & No & Unsuperv. \\ \hline MIPnet & Yes & Bayesian (\& tunable thresh.) & Yes & Superv. + Unsuperv. \\ \hline \hline \end{tabular}
\end{table}
Table 16: Comparison of MIPnet with most related models.

Figure 5: Comparison of MIPnet w and wo unsupervised loss for different missing label percentages.