Learning to Handle Parameter Perturbations in Combinatorial Optimization: an Application to Facility Location

TL;DR

This paper introduces a machine learning approach to adapt solutions of facility location problems to new instances by predicting solution applicability and magnitude of change, enhancing optimization with historical data.

Contribution

It develops a novel method combining classifiers and regressors to incorporate past experience into combinatorial optimization, applicable beyond facility location problems.

Findings

Classifier predicts if reference solution applies to new instance

Regressor estimates the extent of solution change needed

Empirical results show improved adaptation to new problem instances

Abstract

We present an approach to couple the resolution of Combinatorial Optimization problems with methods from Machine Learning, applied to the single source, capacitated, facility location problem. Our study is framed in the context where a reference facility location optimization problem is given. Assuming there exist data for many variations of the reference problem (historical or simulated) along with their optimal solution, we study how one can exploit these to make predictions about an unseen new instance. We demonstrate how a classifier can be built from these data to determine whether the solution to the reference problem still applies to a new instance. In case the reference solution is partially applicable, we build a regressor indicating the magnitude of the expected change, and conversely how much of it can be kept for the new instance. This insight, derived from a priori information, is expressed via an additional constraint in the original mathematical programming formulation. We present an empirical evaluation and discuss the benefits, drawbacks and perspectives of such an approach. Although presented through the application to the facility location problem, the approach developed here is general and explores a new perspective on the exploitation of past experience in combinatorial optimization.

Tables

Table 1. Table 1 : Training Data: repartition between binary classification and regression

Reference	Classification	Regression	Total
A	1434	9283	10000
B	8056	5965	10000
C	7686	6179	10000

Table 2. Table 2 : Comparison of binary classifiers

Reference	Model	Accuracy	False Positive Errors
A	Extra Trees	0.756	0.106
	Neural Net	0.717	0.160
	Logistic	0.648	0.189
	Naive Bayes	0.641	0.191
B	Extra Trees	0.840	0.049
	Neural Net	0.857	0.083
	Naive Bayes	0.782	0.100
	Logistic	0.629	0.192
C	Extra Trees	0.864	0.033
	Neural Net	0.858	0.085
	Naive Bayes	0.776	0.104
	Logistic	0.616	0.199

Table 3. Table 3 : Comparison of regression models by Mean Squared Error

Reference	Linear Regression	Extra Trees	Neural Net
$A$	0.998	0.883	0.975
$B$	2.087	1.780	2.099
$C$	1.114	0.750	1.097

Equations

$$\underset{x,y}{\mathrm{minimize}}$$

$$\sum _{i\in I}{d}_i{x}_{ij}\le s_j{y}_j$$

$$\sum _{j\in J}{x}_{ij}=1$$

$$x_{ij}\in \{0,1\},y_j\in \{0,1\}$$

$$I=\{1,2,\dots ,N_C\},J=\{1,2,\dots ,N_F\},$$

$$\underset{x,y}{\mathrm{minimize}}$$

$$\sum _{r\in R_{ref}}{y}_r\ge \hat{\gamma }{N}_{ref},$$

Full text

1 Introduction

Solving combinatorial optimization problems can be highly challenging. Many problems have been shown to be computationally hard to solve (Garey and Johnson,, 1990), that is, no polynomial-time algorithms exist. Despite this intrinsic complexity, the state-of-the-art modeling and solving solutions can handle real-world instances via exact methods, heuristics or a mix of the two.

Our rationale goes as follows: given a combinatorial problem, we assume the existence of a reference problem instance $P_{ref}$, derived from the description of a system in its nominal state. This reference can be modeled via Mixed Integer Programming (MIP) and solved to its optimum $x_{ref}^{*}$ via state-of-the-art solvers based on branch-and-bound (BB) methods like CPLEX (IBM,, 2018). If a modification of the reference parameters occurs, the problem needs to be solved again. We operate as if we had to respond to this parameter perturbation (also referred to as disruption) under a tight optimization time budget, aiming at accelerating the descent towards the optimum. Given the modifications that have occurred in the past or that can be simulated a priori, one can build a dataset of resolutions. Through the example of the Facility Location Problem, we study how one can extract information from this dataset to facilitate the resolution of future instances. We cast this problem as a Supervised Learning problem to predict whether a modification will affect and to what extent the reference optimal solution. According to this prediction, additional constraints are added to the MIP formulation of the new instance. In order to evaluate the approach, the two formulations are fed to the solver and the performance are compared.

Section 2 reviews the related literature on applying Machine Learning to Optimization problems. Section 3 provides a description of the Facility Location Problem and Section 4 introduces how the prediction problem can be formulated mathematically and linked to existing Machine Learning approaches. Experiments and numerical results are presented and discussed in Section 5. Section 6 discusses the advantages and drawbacks of the proposed method and present future research directions.

2 Machine Learning for Optimization

The remarkable results achieved in recent years by Machine Learning (ML), for example in the field of image recognition (Krizhevsky et al.,, 2012) or Reinforcement Learning (Mnih et al.,, 2015), sparked a lot of interest in the Operations Research community. Some results have emerged in applications such as guiding the branching process in branch-and-bound optimization algorithms, in the form of node exploration ordering (He et al.,, 2014) and approximating performant branching rules (Alvarez et al.,, 2017; Khalil et al.,, 2016; Lodi and Zarpellon,, 2017). In Di Liberto et al., (2016) one can find algorithmic ideas on how to handle the many existing heuristics already developed for branch-and-bound methods.

The traveling salesman problem has been addressed via an ad-hoc neural network architecture in a Reinforcement Learning context (Bello et al.,, 2017; Khalil et al.,, 2017; Deudon et al.,, 2018), aiming at learning to build incrementally a solution.

Another field of application is the interaction between existing optimization solvers and a priori knowledge on the problem. For instance, Kruber et al., (2017) find an approach to automatically handle MIP problem decompositions, by detecting if and what reformulation to apply within a dedicated software. On the same topic, Basso et al., (2018) bring evidence to why such an approach can be carried out, showing empirically that the role of an expert needed when using decomposition methods can be, at least in part, automated via learning algorithms. In Bonami et al., (2018), the authors present a method to select the best resolution method for a Mixed Integer Quadratic Programming problem from the different algorithms offered within the CPLEX framework.

Machine Learning has proven useful in producing a description of the optimal solution of a yet unseen instance. Fischetti and Fraccaro, (2018) describe a real-world application with the problem of a wind turbine park layout. Minimizing the complex turbulence induced by the positioning of wind turbines generates a series of difficult combinatorial problems. Where different configurations would need to be evaluated at a considerable computational cost, the authors approximate the solution values to these optimizations via ML, thanks to a dataset built a priori. On a similar line lies the work of Larsen et al., (2018). When the problem of choosing the optimal load planning for containers on freight trains cannot be solved online because of time constraints and insufficient information (tactical level), a set of offline cases can be collected and used to get a description of the solution of a new instance, at an aggregate level. Such description provides meaningful insights to decision makers in a real-time context.

Strictly connected with our work, Xavier et al., (2019) consider a number of ML techniques to extract information from previously solved instances of Unit Commitment problems and leverage such pieces of information to improve the MIP performance when solving similar instances again and again.

The interested reader is referred to Bengio et al., (2018) for a recent survey on the subject.

3 Constrained Facility Location Problems

Our application of choice will be that of the Capacitated Facility Location Problem (CFLP). Given a group of customers and a list of potential facilities that can serve them, one must open a certain amount of facilities and assign a unique facility to each of the customers. The goal of the problem is to satisfy the demand of the customers while minimizing the associated operational and start-up costs. The facilities are constrained in capacity, that is, the quantity of total service each can provide is limited. When solving this problem one determines which facilities must be opened and which customers they will serve.

3.1 Mathematical Programming Formulation

The CFLP can be written as a Binary Integer Program, with decision variables $y_{j},x_{ij}\in\{0,1\}$, indicating respectively whether facility $j$ is activated (also referred to as open) and whether customer $i$ is served by $j$.

[TABLE]

where $i\in I$ are the customers and $j\in J$ the available facilities; $c_{ij}\geq 0$ is the cost of serving customer $i$ from $j$; $f_{j}\geq 0$ is the fixed cost for opening facility $j$; $d_{i}\geq 0$ is the demand of customer $i$; $s_{j}\geq 0$ is the capacity of facility $j$; $N_{F}$ and $N_{C}$ are the number of facilities and of customers, respectively.

The objective function (1) quantifies the total cost of a given assignment, composed of the fixed start-up costs $\sum_{j\in J}f_{j}y_{j}$ and $\sum_{i\in I}\sum_{j\in J}c_{ij}x_{ij}$, the operational costs. Constraints (2) ensure that the total demand of the customers assigned to $j$ does not exceed its capacity $s_{j}$. Constraints (3) ensure that each customer is served by exactly one facility. For convenience, we write an instance as $P=\left(\boldsymbol{c},\boldsymbol{f},\boldsymbol{d},\boldsymbol{s}\right)$.

In the various forms of the Facility Location Problem (FLP), the binary variables associated with facility activation and customer-to-facility assignment make the problem computationally challenging. Even the uncapacitated version, for instance, has been proven to be NP-hard (Cornuéjols et al.,, 1983). For a comprehensive treatment on exact formulations of the FLP, we refer the reader to the work of Klose and Drexl, (2005). Beyond mathematical programming and exact methods, heuristics for the FLP have received a fair amount of attention (Cornuéjols et al.,, 1983; Guastaroba and Speranza,, 2012, 2014).

3.2 Perturbations in a Reference Problem

We assume the existence of a reference instance $P_{ref}=\left(\boldsymbol{c}_{ref},\boldsymbol{f}_{ref},\boldsymbol{d}_{ref},\boldsymbol{s}_{ref}\right)$ where the parameters are considered to be in their nominal state. This problem has $\{\boldsymbol{x}_{ref}^{*},\boldsymbol{y}_{ref}^{*}\}$ as its reference solution. In case a disruption in the system occurs, for example an anomaly in the capacity, we wish to learn whether the reference solution, or part of it, is still applicable without running the full optimization. To that end, we wish to predict what parts of the solution need to be changed by providing hints to the solver, via an additional constraint. We shall write a new instance $P^{\prime}=\left(\boldsymbol{c}^{\prime},\boldsymbol{f}^{\prime},\boldsymbol{d}^{\prime},\boldsymbol{s}^{\prime}\right)$, which we see as a variation $P^{\prime}=P_{ref}+\Delta P$, where a disruption $\Delta P$ affected one or more of $P_{ref}$’s parameters. In practice, one would proceed as in the following example. The given (simple) $P_{ref}$ with capacities $\boldsymbol{s}_{ref}=(5000,5000,5000)$ is affected by a disruption and the derived $P^{\prime}$ has capacities $\boldsymbol{s}^{\prime}=(5000,5088,4304)$. The ML algorithm will take $\Delta P=\{0,88,-696\}$ as input and predict how $\{\boldsymbol{x^{\prime}}^{*},\boldsymbol{y^{\prime}}^{*}\}$ could differ from $\{\boldsymbol{x}_{ref}^{*},\boldsymbol{y}_{ref}^{*}\}$, before running the optimization for $P^{\prime}$ (see Section 4). In our case, the hint to the solver will specify how many (if any at all) of the facilities active in $P_{ref}$ will be still active, without specifying which of the three.

The FLP can be adapted to different domains. Consider for instance the airport system of a region. If a runway at a major airport has problems or needs maintenance, one could think of this as a disruption. One could then want to learn and predict how to react to the closure of a runway, how much of the traffic would need to be delayed, rescheduled or rerouted or, in general, how to be guided to manage the disruption.

4 Learning constraints

In this Section, we formalize the search for a constraint that will be used to accelerate the resolution of FLP instances facing disruptions. We formulate this as a statistical learning problem and describe how we couple it with MIP resolution.

4.1 Learning to Bound the Changes to a Reference Solution

Given the information gained from past resolutions or simulated data, we investigate the possibility of predicting the number of facilities we expect to change. We want to predict that “only a proportion $\gamma=0.15$ of the facilities needs replanning”, and impose that as a constraint. Indeed, it is well known in the case of MIPs with binary variables (Fischetti and Lodi,, 2003) that it is easy to express the neighborhood of a feasible solution by a linear constraint and the resulting MIP is generally way easier to solve than the original problem. In addition, for FLP it is also well known that the critical choice is associated with the facilities to open while the assignment of clients to facilities is easier to deal with.

Then, the problem boils down to predicting the proportion of facilities opened in $\boldsymbol{y}_{ref}^{*}$ to be kept open in $\boldsymbol{y^{\prime}}^{*}$, the optimal solution of the perturbed instance $P^{\prime}$. This estimation is a scalar value $\widehat{\gamma}$ determined via a function $\Gamma(\Delta P)$, where $\Gamma:\mathbb{R}\rightarrow[0,1]$. Given the $\left\{\left(\Delta P_{k},\gamma_{k}\right)\right\}_{1\leq k\leq K}$ data from the $K$ past resolutions, this is a regression problem, which we can tackle with any Statistical Learning method (Hastie et al.,, 2009).

The drawback of such an approach is the potential introduction of bias in the solution of $P^{\prime}$ if the true optimum of $P^{\prime}$ is far from $\{\boldsymbol{x}_{ref}^{*},\boldsymbol{y}_{ref}^{*}\}$ and the ML fails to detect it or if it cuts the optimum from the admissible domain. In fact, in case the parameters $\left(\boldsymbol{c}^{\prime},\boldsymbol{f}^{\prime},\boldsymbol{d}^{\prime},\boldsymbol{s}^{\prime}\right)$ generating $P^{\prime}$ are very different from $P_{ref}$, a good model would impose no additional constraint, thus falling back to the full problem.

The information on $\widehat{\gamma}=\Gamma(\Delta P)$ can be modeled through the linear constraint (6), thus yielding the new problem

[TABLE]

where $R_{ref}=\{j\in J\mid y_{ref,j}^{*}=1\}$ is the set of indices of the facilities opened in the reference solution, $\widehat{\gamma}\in[0,1]$ is the parameter predicted via ML and $N_{ref}=|R_{ref}|$ is the number of facilities activated in the reference problem.

We do not attempt to identify directly which of the facilities need to change or stay as in $\{\boldsymbol{x}_{ref}^{*},\boldsymbol{y}_{ref}^{*}\}$, but rather to restrict the number of facilities that can be changed. If $\widehat{\gamma}\approx 0$, the disruption data conveys no information as to the new optimum and all the facility allocations should be left to the solver. In that case, it is likely that the variation on $P_{ref}$ is so strong that prior information is of no help and we need to solve the new instance as a full problem in its original formulation. Conversely, if $\widehat{\gamma}\approx 1$, we conclude that all of the reference facilities are to be left activated, without excluding further facilities from being activated.

4.2 Structuring the prediction problem

The overall problem of predicting $\widehat{\gamma}=\Gamma(\Delta P)$ is a regression problem. Previously solved instances $P_{k}$ provide values $(\Delta P_{k},\gamma_{k})$ to assemble a training set. Then, most ML methods will search for $\Gamma$ by minimizing a loss function defined as an expected value of individual losses over the training set, such as the least squares. Such methods can thus be sensitive to the training samples’ distribution and might overfit to the majority value. In the specific case of predicting the proportion of facilities to keep open, the distribution of observed values for samples $(\Delta P_{k},\gamma_{k})$ is strongly affected by the fact that in many cases, all the facilities in the reference solution should remain open (for details, see Section 5.1). To compensate for this imbalance, our predictor’s architecture features two levels.

First, we fit a binary classifier to indicate whether the open facilities in $\{\boldsymbol{x}_{ref}^{*},\boldsymbol{y}_{ref}^{*}\}$ should all be kept open. If this is the case, then $\widehat{\gamma}:=1$ without any further computation. Second, if we cannot assert that all the reference active facilities need to be left active, we call a specific regression function. This $\Gamma(\Delta P)$ function is trained on all the data points except for those for which $\gamma_{k}=1$. Then, given a new problem $P^{\prime}$, we predict how much of the reference solution has to be changed, adding that as constraint (6) to the original problem. Figure 1 summarizes the $\widehat{\gamma}$ estimation process and the resolution of a new instance $P^{\prime}$.

5 Experimental Evaluation

Given the reference instance $P_{ref}$ and an instance $P^{\prime}$ derived by perturbing $P_{ref}$, our aim is to provide a good solution to $P^{\prime}$ within a small time budget.

Our $P_{ref}$ is taken from the CFLP instances of Guastaroba and Speranza, (2014)111Data available at http://or-brescia.unibs.it/instances/instances_sscflp. We test our method on three different reference cases, respectively $capa_{1}$, $capb_{1}$ and $capc_{1}$. For conciseness, we refer to these as $A,B,C$, and their corresponding MIP models as $P_{ref}^{A},P_{ref}^{B},P_{ref}^{C}$. In each case, we have $N_{F}=100$ facilities available that share the same capacity $s^{A}=s_{1}^{A}=s_{2}^{A}=\dots=s_{100}^{A}$.

5.1 Generation of Training Data

For each of the three instances, we produce three datasets of sizes $N_{A}=N_{B}=N_{C}=10000$. Given $P_{ref}$, we generate our training data by perturbing the capacities of the facilities. We randomly perturb between 5% and 50% of the facilities as follows.

Pick uniformly at random between 5% and 50% of the $N_{F}=100$ facilities. 2. 2.

For each facility $j$, apply a Gaussian noise to its capacity:

 $s_{j}\leftarrow s_{j}+\mathcal{N}(0,\sigma=0.2\times s_{j})$. 3. 3.

All the remaining data is left unchanged. This yields $\Delta P_{k}$. 4. 4.

Solve the new instance $P_{k}=P_{ref}+\Delta P_{k}$ with a MIP solver and get the number $Y_{k}$ of facilities that were in $x_{ref}^{*}$ and are also in $x_{k}^{*}$. 5. 5.

Add the new point $(\Delta P_{k},\gamma_{k}=Y_{k}/N_{F})$ to the dataset.

After generating the data, we can observe the distribution of $Y$, the number of original open facilities still open in the optimal solution of the randomized instances. At their nominal values, the optimal solutions for $A$, $B$ and $C$ have respectively 7, 11 and 11 facilities open. For instance, given the $k$-th point in our generated data from case $A$, if $Y_{k}=4$, this means that after the perturbation of the capacities, of the seven facilities open in the reference optimal solution, four were included in the optimal solution of the derived problem $P_{k}$.

Figure 2 presents three different scenarios. With instance $A$ we see that on average less than half of the 7 original facilities are kept open, that is, the perturbations affect the new optimal solutions. In the case of $B$, for most of the data points the original reference solution is still optimal after the perturbations, as seen in the right-most bar of Figure 2(b). Case $C$ is intermediate, with at least eight facilities out of eleven still open after the variations of the capacities.

5.2 Learning the Prediction Model

As mentioned in Section 4.2, we need a binary classifier and a regressor to account for sample distribution imbalance. The number of data for binary classification depends on the shape of the generated data, that is, on how many extreme cases we observe (rightmost column in histograms of Figure 2).

We split the training data into classification and regression data as follows. In case A, we have $717$ data points of value $7$, corresponding to the observations where all the facilities open in the reference solution are to be left open in the perturbed instances. Thus, we set these $717$ points aside and randomly sampled an equal amount of points among the other points, yielding a learning data set of $N_{Classification}A=1434$. The remaining $9283$ points form the regression learning data.

This process was repeated for all the three cases, yielding the distributions reported in Table 1.

In accordance to standard practice in ML (Hastie et al.,, 2009), at learning time the data set was randomly split into a training and a test partition. The predictive models were fit on the training data but evaluated on the test data.

5.2.1 Comparison of ML algorithms: Classification

For the binary classification task, we tested a set of binary classifiers among which Extremely Randomized Trees (Extra Trees) (Geurts et al.,, 2006), Neural Network classifier (Goodfellow et al.,, 2016), Logistic Regression classification and Naive Bayes classification (Hastie et al.,, 2009).

As can be seen across Table 2, Extra Trees emerged as the best performing in terms of accuracy and false positives errors, that is, over the total of the prediction, how many were predicted as false positives. Here, false positives are cases where the reference solution’s open facilities should not be all kept open in the new instance but the classifier prescribes otherwise. This could introduce bias, making this metric important for our task.

5.2.2 Comparison of ML algorithms: Regression

For the regression task we selected Extremely Randomized Trees out of a pool of models comprising Multiple Linear Regression, Extremely Randomized Trees and a Multi-layer Perceptron artificial Neural Network.

In Table 3 are the values of the Mean Squared Error (MSE) computed on the test partition. Extra Trees yielded the best performance of the tested methods. While Neural Networks have recently been the object of intense research and remarkable results, we think its poorer results can explained by the limited amount of training data ($N\leq 10000$) and the highly nonlinear relations between features and response variable, for which more data could have been needed.

Figures 3(a), 3(c) and 3(e), plot the predictions of $Y$ (number of facilities kept open) versus the true values on the test data. The ideal case of perfect predictions $\hat{Y}=Y$, is denoted by the red line. The green line is a simple linear regression of the form $\widehat{Y}=\hat{m}Y+\hat{b}$ that graphically renders the actual relation between $Y$ and $\widehat{Y}$.

We remark that the data are highly concentrated around some central values. The tails of the distribution are scarcely represented in our training data, and the learning function cannot generalize well, as one can see more prominently with instance $A$, where most of the data take values 2, 3 or 4.

5.3 Optimization Results

We report the performance of our method on the three references $A,B$ and $C$, generating three batches $D_{A},D_{B}$ and $D_{C}$ of $50$ new instances, unseen during training. We apply the same procedure followed to generate the learning data (Section 5.1).

To measure the effect of the ML bound with respect to resolution time, each instance in $D_{A}$ (but also $D_{B}$ and $D_{C}$) is optimized twice (with and without bound) with time limits $t_{lim}\in\{5,10,30,60,120\}$ seconds. At the end of the time-limited run, we record the objective value and the effective solving time (some instances are optimized before the limit is reached).

Our ML constraint cannot make the solution infeasible, as the total number of active facilities is not bounded, but it could cut off the optimal solution. For instance, one could predict all of the reference facilities to be open ($\hat{\gamma}=1$) while none of them should be. We run the optimization with 3,600 seconds of time limit, to measure the difference in objective values after a long run, with and without the additional constraint.

5.3.1 Experimental setup

For each iteration of our experiments, we proceed as follows. A new instance is generated as detailed in Section 5.1. The instance is then solved once for each time limit without the additional constraint. The objectives and effective times are recorded. The ML algorithm is run (see Section 4) and constraint (6) is introduced in the model’s formulation. The constrained model is then solved once for each time limit. To verify that the optimal solution was preserved by the introduced bound constraint, we run both models with and without constraint (6) for a full hour to check if any remarkable deviation has been caused by the ML.

When using a highly-optimized MIP solver like CPLEX, it is sometimes hard to evaluate the impact of model modifications (imposed by the user) because the complex algorithmic components within the solver might mitigate the effect of those modifications. For this reason, it is common practice to deactivate, for testing purposes, some of those CPLEX algorithmic ingredients.

In our special case, we are interested in evaluating the effect of the ML-predicted constraint (6) as a way to fast dive toward good solutions of a modified FLP instance with respect to the solution of a reference one. To achieve this goal, we have found useful to run CPLEX without the presolve feature and we first report the results of this configuration. We nonetheless performed all experiments and report the results with CPLEX default, i.e., with presolve activated. We show that the overall message provided by the experiments is in both cases very similar.

Finally, given the interest, especially in the academic community, of using non-commercial MIP solvers, we report in Appendix A the results obtained by replacing CPLEX with the most widely adopted of those solvers, namely SCIP (default version) (Gleixner et al.,, 2017).

5.3.2 Experimental results

Within $3,600$ seconds, most of the runs attain an optimality gap of the order of $0.05\%$ or manage to reach optimality. The bias introduced by the learned bound is on average around $0.03\%$ and at most of the order of $0.5\%$. That is, at the $3,600$ seconds time limit, the objective value of the ML-bounded run is on average $0.03\%$ higher than the objective for the regular run, but never greater than $0.5\%$.

Figures 4(a), 4(c) and 4(e) compare the values of the objective functions of the constrained and unconstrained resolutions at the time limit. Each point represents the same new, unseen, instance optimized twice: once with the ML bound and once without. The color coding helps determining the time limit imposed on this pair of runs. On the $x$ axis one can see the objective value without ML, on the $y$ axis the value with the ML bound. Points on the lower left corner correspond to runs with longer time limits, where the difference between the two approaches becomes more negligible. The points below the line are runs for which our method outperformed the original formulation.

In Figures 4(b), 4(d) and 4(f), the same information is presented focusing on the temporal dimension. On the $x$ axis is the effective resolution time, which can be inferior to the time limit. Along the $y$ axis is the objective value at the time limit. The points aligned at the bottom of Figures 4(d) and 4(f) are those for which the optimal solution was found before the time limit expired. For short resolution times, the ML-constrained formulation yields better objective function values than the original formulation without cutting the optimal solution which is eventually found given more computational time. On average, the blue points representing the objective of the ML-bounded model lie below the red ones of the unbounded model. In case $C$, our method allows for superior performances than those obtained simply via the solver. In case $A$ we are still on average improving the performances or at least keeping in line with the unbounded model. Here, the amplitude of the perturbations imply completely different solutions, making transfer difficult. We remark however that, on average, our approach avoids cutting the optimal solution and thus prevents negative transfer.

As noted above, we operated with the presolve feature switched off. For completeness, we run all of the experiments a second time with the presolve feature activated (Figure 5). The results reflect what was observed previously, although the presolve heuristics of CPLEX manage to reduce in part the effect of our approach in case $B$ (Figure 5 c and d).

6 Conclusions and Perspectives

We have shown that the existing or simulated data of a recurring optimization problem can be exploited to gain insight into the optimization process. We have used these data to fit a binary classifier and a regressor. We predict whether a new instance, derived from a reference instance, will share all or a part of its optimal solution with the solution of the reference instance. This piece of information, translated into an additional constraint, is given to a solver and allows to dive faster, on average, towards a good solution. Moreover, we empirically illustrated that this additional constraint preserves the optimal solution and thus prevents negative transfer.

Building upon these results, we plan on extending the reach of our approach by experimenting with other families of problems, where, potentially, different types of constraints should be considered. This will be the way of fully demonstrating the wide applicability of the approach to both solving repetition (with different data) and re-optimization.

Acknowledgements

The authors wish to acknowledge the support of the Canada Excellence Research Chair in Data Science for Real-Time Decision-Making at Polytechnique Montréal and the ISAE-SUPAERO foundation.

Appendix A Results with SCIP

The MIP solver SCIP, whose source code is accessible for academic purposes, has a considerable adoption in academic research, which is why it is interesting to test our approach in conjunction with such a solver. The experiments reported in Figure 6 are the same as those in Section 5.3, with the exception of the time limits of $5$ seconds; for such a short time limit, most of the runs did not yield a feasible solution and we did not report it. What observed previously with CPLEX is true here as well, although the effect of constraint (6) are amplified. The additional bound constraint allows the solver to dive much faster towards a good solution. Even at the 120 seconds mark, the gap between the two objective values is still considerable.